! Code Saturne User Subroutines
! Swirl Inlet Velocity Profile

! This program is free software (including commercial use) and distributed under GNU GPL v3 license.

! Function calculates local velocity components at the boundary mesh face
! Program is only compatible with Linux (determining executable path)

! [!] Temperature column may present in the table only if non-compressible model is used with temperature scalar.
! For compressible model or with enthalpy scalar setting temperature ! will have no effect or will induce errors.
! In most cases you just don't need to set the temperature profile.

! Input file "SwirlProfile.csv" should be placed in DATA directory and
! has strict CSV format with features listed below.
! 1. First line contents coordinates of inlet center with "Cnt" identifier in first column
!    (format is Cnt,x,y,z).
! 2. Second line defines the name of inlet, case sensitive (format is InletMnem,Inlet).
! 3. Third line defines multiplier for normal velocity component (format is VelNrmMul,1.0).
! 4. Fourth line contents table header. Each column has parameter identifier:
!    Rad => Radius [m]
!    VelNrm => Axial velocity [m/s]
!    VelTan => Tangential velocity [m/s]
!    VelRad => Radial velocity [m/s]
!    Temp => Temperature [C]
!    TurbEng => Turbulent kinetic energy
!    TurbDiss => Turbulent energy dissipation
!    Identifiers are cese-sensitive.
!    Column headers may have extra characters like units.
! 5. Next and other lines contents parameter values without extra characters (except spaces).
! 6. Supported field (value) separators are "," and ";",
!    new line magic numbers expected in Win/DOS or Linux formats (0D, 0A, 0A0D, 0D0A).
! 7. Decimal separator must be "." (point), no other decimal separators allowed.
! 8. There may be empty lines (without separators), but empty fields and comment lines are not allowed!
! 9. There may be additional trailing seperetors at line ends. It's a standard side-effect
!    during table conversion to CSV format.

! Common parameters and structures
module tdcCommon
integer,parameter::nDim=3,& ! Number of dimentions
                   ! Direction identifiers
                   dirNrm=1,& ! Normal
                   dirTang=2,& ! Tangential
                   dirRad=3,& ! Radial
                   dirX=1,& ! Axis X
                   dirY=2,& ! Axis Y
                   dirZ=3 ! Axis Z
double precision,parameter::eps=1.0d-50 ! Relative real tolerance
end module tdcCommon

! Variant parameters and structures
module tdcSwirl
use tdcCommon ! Common parameters
integer,parameter::& ! Array dimensions
                   nSwirlMdlTabNodeMax=10000,& ! Nodes in table
                   ! Identifiers
                   ! Table file reading mode
                   idSwirlFileTabReadModeCnt=1,& ! Center coordinates reading
                   idSwirlFileTabReadModePar=2,& ! Parameters reading
                   idSwirlFileTabReadModeHeader=3,& ! Header reading
                   idSwirlFileTabReadModeValue=4 ! Values reading
! Coordinate system vector components
type SwirlMdlCoordVct
 double precision & ! Vector components [m]
                  nrm(nDim),& ! Normal
                  rad(nDim),& ! Radial
                  tang(nDim) ! Tangential
end type SwirlMdlCoordVct
! Table node parameters
type SwirlMdlTabNode
 double precision rad,& ! Node radial position [m]
                  ! Tabulated parameter values
                  velNrm,& ! Velocity normal [m/s]
                  velTang,& ! Velocity tangential [m/s]
                  velRad,& ! Velocity radial [m/s]
                  temp,& ! Temperature [C]
                  turbEng,& ! Turbulence kinetic energy
                  turbDiss ! Turbulence energy dissipation
end type SwirlMdlTabNode
! Tabulated values
type SwirlMdlTab
 integer nNode,& ! Number of table nodes
         indNodeSel,& ! Index of current (selected) node
         ! Parameter column indexes
         indRadCol,& ! Radius
         indVelNrmCol,& ! Velocity normal
         indVelTangCol,& ! Velocity tangential
         indVelRadCol,& ! Velocity radial
         indTempCol,& ! Temperature
         indTurbEngCol,& ! Turbulence kinetic energy
         indTurbDissCol ! Turbulence energy dissipation
 type (SwirlMdlTabNode) node(nSwirlMdlTabNodeMax) ! Node paramenetrs
end type SwirlMdlTab
! Local velocity parameters in stationary coordinate system
type SwirlMdlVelLocalStat
 double precision comp(nDim) ! Components [m/s]
end type SwirlMdlVelLocalStat
! Local velocity parameters in rotating coordinate system
type SwirlMdlVelLocalRotate
 double precision nrmMul ! Normal velocity multiplier
 double precision comp(nDim) ! Components [m/s]
end type SwirlMdlVelLocalRotate
! Local velocity parameters (at current radius)
type SwirlMdlVelLocal
 double precision mag ! Velocity magnitude [m/s]
 type (SwirlMdlVelLocalStat) stat ! Stationary (main) coordinate system
 type (SwirlMdlVelLocalRotate) rotate ! Rotating coordinate system
end type SwirlMdlVelLocal
! Velocity parameters
type SwirlMdlVel
 type (SwirlMdlVelLocal) local ! Local velocity parameters (at current radius)
end type SwirlMdlVel
! Inlet centroid coordinates
type SwirlMdlCnt
 double precision coord(nDim) ! Coordinates in stationary system [m]
end type SwirlMdlCnt
! Mesh face on inlet surface for velocity calculation
type SwirlMdlFace
 double precision nrmVct(nDim) ! Normal vector, lenght is a face area
end type SwirlMdlFace
! Point on inlet surface for velocity calculation
type SwirlMdlPoint
 double precision coord(nDim),& ! Coordinates in stationary system [m]
                  rad ! Radius relative to inlet centroid [m]
 type (SwirlMdlCoordVct) rcVct ! Rotating coordinate system vector components in stationary system
end type SwirlMdlPoint
! Swirl inlet model parameters
type SwirlMdl
 character(255) inletMnem ! Mnemonic name of inlet to set velocity profile
 type (SwirlMdlCnt) cnt ! Inlet centroid coordinates
 type (SwirlMdlFace) face ! Face on inlet surface for velocity calculation
 type (SwirlMdlPoint) point ! Point on inlet surface for velocity calculation
 type (SwirlMdlTab) tab ! Tabulated velocity components (function of radius)
 type (SwirlMdlVel) vel ! Velocity parameters
end type SwirlMdl
! Swirl inlet parameters
type Swirl
 type (SwirlMdl) mdl ! Model
end type Swirl
type (Swirl) prog ! Program parameters
end module tdcSwirl

! Print simple text
subroutine TextPrint(text)
implicit none
character(*) text
write (*,*) text
end subroutine TextPrint

! Real value equal or greater than zero
double precision function GeZeroR(arg)
use tdcCommon
implicit none
double precision arg
GeZeroR=max(arg,0.0d0) ! Function value
end function GeZeroR

! Real value equal or greater than zero
double precision function GtZeroR(arg)
use tdcCommon
implicit none
double precision arg
GtZeroR=max(arg,eps) ! Function value
end function GtZeroR

! Non-zero real value
double precision function NonZeroR(arg)
use tdcCommon
implicit none
double precision arg
if ((arg.le.-eps).or.(arg.ge.eps)) then ! Allowed value ranges
 NonZeroR=arg ! Function value
else ! Absolute value is less than tolerance
 if (arg.lt.0) then ! Negative value
  NonZeroR=-eps ! Function value
 else ! Positive value
  NonZeroR=eps ! Function value
 end if
end if
end function NonZeroR

! Vector length (assuming size is nDim)
double precision function VctLen(vct)
use tdcCommon
implicit none
double precision vct(*)
integer iCoord
double precision,external::GeZeroR
VctLen=0 ! Initialise vector length
do iCoord=1,nDim ! Coordinate loop
 VctLen=VctLen+vct(iCoord)**2 ! Vector length quadrature sum [m2]
end do
VctLen=GeZeroR(VctLen)**0.5d0 ! Vector length [m]
end function VctLen

! Vector normalize (assuming size is nDim)
subroutine VctNorm(vctIn,vct)
use tdcCommon
implicit none
double precision vctIn(*),vct(*)
integer iCoord
double precision length
double precision,external::VctLen,GtZeroR
length=VctLen(vctIn) ! Vector length
do iCoord=1,nDim ! Coordinate loop
 vct(iCoord)=vctIn(iCoord)/GtZeroR(length) ! Normalized vector component [m]
end do
end subroutine VctNorm

! Vector product (assuming size is nDim)
subroutine VctMulVct(vctPri,vctSec,vctMul)
use tdcCommon
implicit none
double precision vctPri(*),vctSec(*),vctMul(*)
! Components of resulting vector
vctMul(dirX)=vctPri(dirY)*vctSec(dirZ)-vctPri(dirZ)*vctSec(dirY) ! Component X
vctMul(dirY)=vctPri(dirZ)*vctSec(dirX)-vctPri(dirX)*vctSec(dirZ) ! Component Y
vctMul(dirZ)=vctPri(dirX)*vctSec(dirY)-vctPri(dirY)*vctSec(dirX) ! Component Z
end subroutine VctMulVct

! Convert string to number
double precision function StrToNum(str)
implicit none
character(*) str
integer idState
double precision num
num=0 ! Initialise number
read (str,*,iostat=idState) num ! Read number from string
StrToNum=num ! Function value
end function StrToNum

! Convert number to string
character(*) function NumToStr(num,fmt)
implicit none
double precision num
character(*) fmt
character(255) str
str="" ! String
write (str,fmt) num ! Write number to string
NumToStr=adjustl(str) ! Function value
end function NumToStr

! Get program executabe (binary) directory
! Program directory not used in library version of this program
logical function ProgDirGet(progDirName)
implicit none
character(*) progDirName
integer indSlashPos
character(255) progFileName
! Initialisation
ProgDirGet=.false. ! Function value
! Determining program directory
! call getarg(0,progFileName) ! Full program file name (with path)
indSlashPos=index(progFileName,'/',.true.) ! Leading slash position
if (indSlashPos.le.0) then ! Slash not found in program path
 call TextPrint("[!] Error parsing Swirl program directory name!")
 return
end if
if (indSlashPos.gt.len(progDirName)) then ! Program path is too long
 call TextPrint("[!] Swirl profile directory name is too long!")
 return
end if
progDirName=progFileName(1:indSlashPos) ! Program directory name
ProgDirGet=.true. ! Function value
end function ProgDirGet

! Read tabulated velocity components (rotating coordinate system)
! Subroutine implies that the file "SwirlProfile.csv" has a CSV format
logical function SwirlFileTabRead(mdl)
use tdcSwirl
implicit none
type (SwirlMdl) mdl
integer,parameter::nSep=2,nLf=2
integer iChar,iSep,iLf
integer indFile,indField,indValRow,idState,idMode
integer(1) bt
integer(1) lf(nLf)
logical isSuccess,isFileExist,isSep,isLf,isLfPrv,doFileRead
character(1) chr
character(1) sep(nSep)
character(255) progDirName,fileName,fieldStr,parMnem
logical,external::ProgDirGet
double precision,external::StrToNum
! Initialisation
SwirlFileTabRead=.false. ! Function value
mdl%inletMnem="" ! Inlet mnemonic name
mdl%vel%local%rotate%nrmMul=-1 ! Normal velocity multiplier
mdl%tab%indVelNrmCol=-1 ! Normal velocity column index
mdl%tab%indVelTangCol=-1 ! Tangential velocity column index
mdl%tab%indVelRadCol=-1 ! Radial velocity column index
mdl%tab%indTempCol=-1 ! Temperature column index
mdl%tab%indTurbEngCol=-1 ! Turbulence kinetic energy column index
mdl%tab%indTurbDissCol=-1 ! Turbulence energy dissipation column index
! Separators
sep(1)="," ! First separator
sep(2)=";" ! Second separator
! Line feed marks
lf(1)=13 ! First mark (0Dh)
lf(2)=10 ! Second mark (0Ah)
! Reading file
! isSuccess=ProgDirGet(progDirName) ! Program executable directory
! fileName=trim(progDirName)//"SwirlProfile.csv" ! File name
fileName="SwirlProfile.csv" ! Base file name
inquire(file=fileName,exist=isFileExist) ! Check if file exist
if (.not.isFileExist) then ! File not found
 call TextPrint("[!] File ["//trim(fileName)//"] not found!")
 return
end if
indFile=1 ! File index
open (unit=indFile,file=fileName,action='read',form='unformatted',access='stream',iostat=idState) ! Open file
if (idState.ne.0) then ! Error opening file
 call TextPrint("[!] Error opening file ["//trim(fileName)//"]!")
 return
end if
doFileRead=.true. ! File read loop indicator
idMode=idSwirlFileTabReadModeCnt ! Read mode identifier
indField=1 ! Field index
fieldStr="" ! Current field string
parMnem="" ! Parameter mnemonic name
indValRow=1 ! Value row index
isLfPrv=.true. ! Line feed flag
do while (doFileRead) ! File read loop
 read (indFile,iostat=idState) bt ! Read current byte from file
 if (idState.lt.0) then ! Check end of file
  ! doFileRead=.false. ! Read loop indicator
  close (indFile) ! Close file
  exit ! Exit loop
 end if
 if (idState.gt.0) then ! Error reading file
  call TextPrint("[!] Error reading Swirl profile!")
  close (indFile) ! Close file
  return
 end if
 chr=char(bt) ! Current character
 ! Check for field separator
 isSep=.false. ! Separator flag
 do iSep=1,nSep ! Separator compare loop
  if (chr.eq.sep(iSep)) then ! Current character is a separator
   isSep=.true. ! Separator flag
   exit ! Exit loop
  end if
 end do
 ! Check for line feed
 isLf=.false. ! Line feed flag
 if (.not.isSep) then ! Current character is not a separator
  do iLf=1,nLf ! Line feed compare loop
   if (bt.eq.lf(iLf)) then ! Current character is a line feed
    isLf=.true. ! Line feed flag
    exit ! Exit loop
   end if
  end do
 end if
 ! Process current character or field
 if ((.not.isSep).and.(.not.isLf)) then ! Regular character
  fieldStr=trim(fieldStr)//chr ! Add character to current field
 else ! Separator or line feed
  if ((fieldStr.ne."").and.(.not.isLfPrv)) then ! Field is not empty
   select case (idMode) ! Select mode
   case (idSwirlFileTabReadModeCnt) ! Center coordinates
    if (indField.eq.1) then ! First field
     if (index(fieldStr,"Cnt").le.0) then ! Magic string "Cnt" not found
      call TextPrint("[!] Incorrect first field in Swirl file, keyword [Cnt] missing!")
      close (indFile) ! Close file
      return
     end if
    else ! Center coordinate field
     if (indField.le.(nDim+1)) then ! Coordinate field
      mdl%cnt%coord(indField-1)=StrToNum(fieldStr) ! Set current center coordinate [m]
     else ! Surplus (trailing) field
      ! Just read trailing fields to move file pointer to next line
     end if
    end if
   case (idSwirlFileTabReadModePar) ! Parameter
    if (indField.eq.1) then ! First (parameter name) field
     parMnem=fieldStr ! Parameter mnemonic name
    else ! Intermediate (parameter value) field
     if (len_trim(fieldStr).gt.0) then ! Valid parameter field
      if (index(parMnem,"InletMnem").gt.0) then ! Inlet mnemonic name
       mdl%inletMnem=fieldStr ! Set inlet mnemonic name
      end if
      if (index(parMnem,"VelNrmMul").gt.0) then ! Normal velocity multiplier
       mdl%vel%local%rotate%nrmMul=StrToNum(fieldStr) ! Set velocity multiplier
       if (mdl%vel%local%rotate%nrmMul.lt.(0.01)) then ! Incorrect multiplier
        call TextPrint("[!] Incorrect multiplier for normal velocity in Swirl parameters!")
        close (indFile) ! Close file
        return
       end if
      end if
      ! parMnem="" ! Clear parameter mnemonic name
     end if
    end if
   case (idSwirlFileTabReadModeHeader) ! Header
    if (indField.eq.1) then ! First header field (radius)
     if (index(fieldStr,"Rad").le.0) then ! Magic string "Rad" not found
      call TextPrint("[!] Incorrect first field in Swirl table, keyword [Rad] missing!")
      close (indFile) ! Close file
      return
     end if
     mdl%tab%indRadCol=1 ! Set first column as radius
    else ! Intermeduate header field (parameter name)
     ! Set matching parameters column indexes
     if (index(fieldStr,"VelNrm").gt.0) mdl%tab%indVelNrmCol=indField ! Normal velocity column index
     if (index(fieldStr,"VelTan").gt.0) mdl%tab%indVelTangCol=indField ! Tangential velocity column index
     if (index(fieldStr,"VelRad").gt.0) mdl%tab%indVelRadCol=indField ! Radial velocity column index
     if (index(fieldStr,"Temp").gt.0) mdl%tab%indTempCol=indField ! Temperature column index
     if (index(fieldStr,"TurbEng").gt.0) mdl%tab%indTurbEngCol=indField ! Turbulence kinetic energy column index
     if (index(fieldStr,"TurbDiss").gt.0) mdl%tab%indTurbDissCol=indField ! Turbulence energy dissipation column index
    end if    
   case (idSwirlFileTabReadModeValue) ! Values
    ! Set matching parameter value
    if (indField.eq.mdl%tab%indRadCol) mdl%tab%node(indValRow)%rad=StrToNum(fieldStr) ! Radius [m]
    if (indField.eq.mdl%tab%indVelNrmCol) mdl%tab%node(indValRow)%velNrm=StrToNum(fieldStr) ! Normal velocity [m/s]
    if (indField.eq.mdl%tab%indVelTangCol) mdl%tab%node(indValRow)%velTang=StrToNum(fieldStr) ! Tangential velocity [m/s]
    if (indField.eq.mdl%tab%indVelRadCol) mdl%tab%node(indValRow)%velRad=StrToNum(fieldStr) ! Radial velocity [m/s]
    if (indField.eq.mdl%tab%indTempCol) mdl%tab%node(indValRow)%temp=StrToNum(fieldStr) ! Temperature [C]
    if (indField.eq.mdl%tab%indTurbEngCol) mdl%tab%node(indValRow)%turbEng=StrToNum(fieldStr) ! Turbulence kinetic energy
    if (indField.eq.mdl%tab%indTurbDissCol) mdl%tab%node(indValRow)%turbDiss=StrToNum(fieldStr) ! Turbulence energy dissipation
   end select
   indField=indField+1 ! Field index
   fieldStr="" ! Initialize field string
   if (indField.gt.100) then ! Unrealistic number of columns (incorrect file)
    call TextPrint("[!] Incorrect Swirl file, more than 100 columns found!")
    close (indFile) ! Close file
    return
   end if
  end if
 end if
 if (isLf.and.(.not.isLfPrv)) then ! Line feed
  select case (idMode) ! Select read mode
  case (idSwirlFileTabReadModeHeader) ! Header mode
   idMode=idSwirlFileTabReadModeValue ! Change mode to values
  case (idSwirlFileTabReadModePar) ! Parameters mode
   if (index(parMnem,"VelNrmMul").gt.0) idMode=idSwirlFileTabReadModeHeader ! Change mode to header
  case (idSwirlFileTabReadModeCnt) ! Center coordinates mode
   idMode=idSwirlFileTabReadModePar ! Change mode to parameters
  case (idSwirlFileTabReadModeValue) ! Parameter value mode
   mdl%tab%nNode=indValRow ! Number of table nodes
   if (indField.gt.2) indValRow=indValRow+1 ! Parameter value row index (only non-empty rows)
   if (indValRow.gt.nSwirlMdlTabNodeMax) then ! Table does not fit in program array
    call TextPrint("[!] Swirl table contains too many rows abd cannot be read!")
    close (indFile) ! Close file
    return
   end if
  end select
  indField=1 ! Initialise field index
 end if
 isLfPrv=isLf ! Line feed flag
end do
close (indFile) ! Close file
! Check parameters
if (len_trim(mdl%inletMnem).eq.0) then ! Inlet name not provided
 call TextPrint("[!] Inlet mnemonic name for velocity swirl is empty!")
 return
end if
if (mdl%vel%local%rotate%nrmMul.lt.0) then ! Incorrect normal velocity multiplier
 mdl%vel%local%rotate%nrmMul=0 ! Set zero multiplier
 call TextPrint("[!] Multiplier for normal velocity in Swirl parameters not defined!")
 return
end if
SwirlFileTabRead=.true. ! Function value
end function SwirlFileTabRead

! Print summary info
subroutine SwirlIfcSummaryPrint(mdl)
use tdcSwirl
implicit none
type (SwirlMdl) mdl
integer iNode
integer nNode
double precision velNrmMin,velNrmMax,velTangMin,velTangMax,velRadMin,velRadMax
character(255),external::NumToStr
! Initialisation
! Normal velocity
velNrmMin=1.0d10 ! Minimum [m/s]
velNrmMax=-1.0d10 ! Maximum [m/s]
! Tangential velocity
velTangMin=1.0d10 ! Minimum [m/s]
velTangMax=-1.0d10 ! Maximum [m/s]
! Radial velocity
velRadMin=1.0d10 ! Minimum [m/s]
velRadMax=-1.0d10 ! Maximum [m/s]
! Calculate velocity component ranges
nNode=mdl%tab%nNode ! Number of table nodes
do iNode=1,nNode ! Table node cycle
 ! Correct velocity ranges
 ! Normal velocity
 if (velNrmMin.gt.mdl%tab%node(iNode)%velNrm) velNrmMin=mdl%tab%node(iNode)%velNrm ! Minimum [m/s]
 if (velNrmMax.lt.mdl%tab%node(iNode)%velNrm) velNrmMax=mdl%tab%node(iNode)%velNrm ! Maximum [m/s]
 ! Tangential velocity
 if (velTangMin.gt.mdl%tab%node(iNode)%velTang) velTangMin=mdl%tab%node(iNode)%velTang ! Minimum [m/s]
 if (velTangMax.lt.mdl%tab%node(iNode)%velTang) velTangMax=mdl%tab%node(iNode)%velTang ! Maximum [m/s]
 ! Radial velocity
 if (velRadMin.gt.mdl%tab%node(iNode)%velRad) velRadMin=mdl%tab%node(iNode)%velRad ! Minimum [m/s]
 if (velRadMax.lt.mdl%tab%node(iNode)%velRad) velRadMax=mdl%tab%node(iNode)%velRad ! Maximum [m/s]
end do
! Print summary info
call TextPrint("") ! Empty line
call TextPrint("[=] Swirl velocity inlet parameters summary")
call TextPrint("===========================================")
call TextPrint("[i] Inlet name ["//trim(mdl%inletMnem)//"], Normal velocity coefficient ["&
               //trim(NumToStr(mdl%vel%local%rotate%nrmMul,"(f10.5)"))//"]")
call TextPrint("[i] Inlet center ["&
               //trim(NumToStr(mdl%cnt%coord(dirX),"(f10.3)"))//" m, "&
               //trim(NumToStr(mdl%cnt%coord(dirY),"(f10.3)"))//" m, "&
               //trim(NumToStr(mdl%cnt%coord(dirZ),"(f10.3)"))//" m], Radius ["&
               //trim(NumToStr(mdl%tab%node(1)%rad,"(f10.3)"))//"..."&
               //trim(NumToStr(mdl%tab%node(nNode)%rad,"(f10.3)"))//" m]")
call TextPrint("[i] Vel Nrm ["//trim(NumToStr(velNrmMin,"(f10.2)"))//&
               "..."//trim(NumToStr(velNrmMax,"(f10.2)"))//&
               " m/s], Tang ["//trim(NumToStr(velTangMin,"(f10.2)"))//&
               "..."//trim(NumToStr(velTangMax,"(f10.2)"))//&
               " m/s], Rad ["//trim(NumToStr(velRadMin,"(f10.2)"))//&
               "..."//trim(NumToStr(velRadMax,"(f10.2)"))//" m/s]")
call TextPrint("===========================================")
call TextPrint("") ! Empty line
end subroutine SwirlIfcSummaryPrint

! Select table node index by relative (radial) point coordinate
logical function SwirlMdlTabSelect(mdl)
use tdcSwirl
implicit none
type (SwirlMdl) mdl
integer iNode
double precision radMin,radMax
! Initialisation
SwirlMdlTabSelect=.false. ! Function value
mdl%tab%indNodeSel=0 ! Selected table node
! Table search
radMin=mdl%tab%node(1)%rad ! Table radius minimum [m]
radMax=mdl%tab%node(mdl%tab%nNode)%rad ! Table raduis maximum [m]
if (mdl%point%rad.le.radMin) then ! Point radius is lower than table minimum
 mdl%tab%indNodeSel=1 ! Selected table node
else ! Point radius is greater than table minimum
 if (mdl%point%rad.ge.radMax) then ! Point radius is greater than table maximum
  mdl%tab%indNodeSel=mdl%tab%nNode ! Selected table node
 else ! Point radius is inside table range
  do iNode=1,mdl%tab%nNode-1 ! Table node loop
   if ((mdl%point%rad.ge.mdl%tab%node(iNode)%rad).and.&
       (mdl%point%rad.le.mdl%tab%node(iNode+1)%rad)) then ! Point is between current nodes
    if ((mdl%point%rad-mdl%tab%node(iNode)%rad).lt.&
        (mdl%tab%node(iNode+1)%rad-mdl%point%rad)) then ! Nearest node is lower
     mdl%tab%indNodeSel=iNode ! Selected table node
    else ! Nearest node is upper
     mdl%tab%indNodeSel=iNode+1 ! Selected table node
    end if
    exit ! Exit loop   
   end if
  end do
 end if
end if
if (mdl%tab%indNodeSel.le.0) return ! Check node selection
SwirlMdlTabSelect=.true. ! Function value
end function SwirlMdlTabSelect

! Calculate point radial coordinate and locar rotating axes vectors
logical function SwirlMdlPointCoordGet(mdl)
use tdcSwirl
implicit none
type (SwirlMdl) mdl
integer iCoord
double precision,external::GeZeroR,GtZeroR,VctLen
! Initialisation
SwirlMdlPointCoordGet=.false. ! Function value
! Calculate coordinates
! Radial coordinate
mdl%point%rad=0 ! Point radial coordinate [m]
do iCoord=1,nDim ! Coordinate loop
 mdl%point%rcVct%rad(iCoord)=mdl%point%coord(iCoord)-mdl%cnt%coord(iCoord) ! Radial vector component [m]
 mdl%point%rad=mdl%point%rad+mdl%point%rcVct%rad(iCoord)**2 ! Radius quadrature sum [m2]
end do
mdl%point%rad=GeZeroR(mdl%point%rad)**0.5d0 ! Point radial coordinate [m]
! Radial vector (from point to center)
do iCoord=1,nDim ! Coordinate loop
 mdl%point%rcVct%rad(iCoord)=mdl%point%rcVct%rad(iCoord)/GtZeroR(mdl%point%rad) ! Radius quadrature sum [m2]
end do
! Normal vector
call VctNorm(mdl%face%nrmVct,mdl%point%rcVct%nrm) ! Just normalize face normal vector
! Tangential vector
call VctMulVct(mdl%point%rcVct%nrm,mdl%point%rcVct%rad,mdl%point%rcVct%tang) ! Vector product of normal and radial vectors
SwirlMdlPointCoordGet=.true. ! Function value
end function SwirlMdlPointCoordGet

! Calculate velocity vector components in stationary coordinate system
logical function SwirlMdlVelLocalStatGet(mdl)
use tdcSwirl
implicit none
type (SwirlMdl) mdl
integer iCoord
integer indTabNode
double precision velNrm
double precision,external::GeZeroR,GtZeroR
! Initialisation
SwirlMdlVelLocalStatGet=.false. ! Function value
mdl%vel%local%rotate%comp(:)=0 ! Velocity components in rotating system
! Get vector components from table
indTabNode=mdl%tab%indNodeSel ! Index of selected table row (current radius)
if (indTabNode.gt.0) then ! Table row (node) selected
 ! Velocity components [m/s]
 if (mdl%tab%indVelNrmCol.gt.0) mdl%vel%local%rotate%comp(dirNrm)=mdl%tab%node(indTabNode)%velNrm ! Normal
 if (mdl%tab%indVelTangCol.gt.0) mdl%vel%local%rotate%comp(dirTang)=mdl%tab%node(indTabNode)%velTang ! Tangential
 if (mdl%tab%indVelRadCol.gt.0) mdl%vel%local%rotate%comp(dirRad)=mdl%tab%node(indTabNode)%velRad ! Radial
end if
velNrm=mdl%vel%local%rotate%nrmMul*mdl%vel%local%rotate%comp(dirNrm) ! Corrected normal velocity [m/s]
! Calculate vector components [m/s]
! Component along axis X
mdl%vel%local%stat%comp(dirX)=mdl%point%rcVct%nrm(dirX)*velNrm+&
                              mdl%point%rcVct%tang(dirX)*mdl%vel%local%rotate%comp(dirTang)+&
                              mdl%point%rcVct%rad(dirX)*mdl%vel%local%rotate%comp(dirRad)
! Component along axis Y
mdl%vel%local%stat%comp(dirY)=mdl%point%rcVct%nrm(dirY)*velNrm+&
                              mdl%point%rcVct%tang(dirY)*mdl%vel%local%rotate%comp(dirTang)+&
                              mdl%point%rcVct%rad(dirY)*mdl%vel%local%rotate%comp(dirRad)
! Component along axis Z
mdl%vel%local%stat%comp(dirZ)=mdl%point%rcVct%nrm(dirZ)*velNrm+&
                              mdl%point%rcVct%tang(dirZ)*mdl%vel%local%rotate%comp(dirTang)+&
                              mdl%point%rcVct%rad(dirZ)*mdl%vel%local%rotate%comp(dirRad)
SwirlMdlVelLocalStatGet=.true. ! Function value
end function SwirlMdlVelLocalStatGet

! Calculate local velocity components in global (stationary domain) coordinate system
! Main input variables
! 1. Tabulated velocity components (axial, tangential, radial) in rotating system
! 2. Inlet centroid coords (stationary system)
! 3. Current point coords (stationary system)
logical function SwirlMdlPointVelCompGet(mdl)
use tdcSwirl
implicit none
type (SwirlMdl) mdl
logical isSuccess
logical,external::SwirlMdlPointCoordGet,SwirlMdlTabSelect,SwirlMdlVelLocalStatGet
! Initialize
SwirlMdlPointVelCompGet=.false. ! Function value
! Prepare coordinates and find table values
isSuccess=SwirlMdlPointCoordGet(mdl) ! Determine local rotating coordinate system directional vectors
if (.not.isSuccess) return ! Check local coordinates directions determination
isSuccess=SwirlMdlTabSelect(mdl) ! Determine radial position and select values from table
if (.not.isSuccess) return ! Check table index determination
! Convert normal (axial), tangential and radial velocity components to stationary system
isSuccess=SwirlMdlVelLocalStatGet(mdl) ! Get velocity components in stationary coordinate system
if (.not.isSuccess) return ! Check coordinate transformation
SwirlMdlPointVelCompGet=.true. ! Function value
end function SwirlMdlPointVelCompGet

! Set Inlet boundary condition for Code Saturne model
! This function is based on example routine (Base, EDF S. A.)
! distributed under GNU GPL v2 (or later) license
subroutine cs_f_user_boundary_conditions(nvar,nscal,icodcl,itrifb,itypfb,izfppp,dt,rcodcl)
! Code Saturne modules
use paramx
use numvar
use optcal
use cstphy
use cstnum
use entsor
use parall
use period
use ihmpre
use ppppar
use ppthch
use coincl
use cpincl
use ppincl
use ppcpfu
use atincl
use atsoil
use ctincl
use cs_fuel_incl
use mesh
use field
use turbomachinery
use iso_c_binding
use cs_c_bindings
! Variant modules
use tdcSwirl,nDimSwirl=>nDim
implicit none
! Arguments
integer nvar,nscal
integer itrifb(nfabor),itypfb(nfabor),izfppp(nfabor)
integer icodcl(nfabor,nvarcl)
double precision dt(ncelet)
double precision rcodcl(nfabor,nvarcl,3)
! Local variables
integer iCoord,iFace,iCell
integer indFace,indCell,indTabNode,nFace
integer,allocatable,dimension(:)::faceDirList
logical isSuccess
logical,save::isStart=.true.,isFileSuccess=.true.
double precision temp,turbEng,turbDiss,turbEngCoef
logical,external::SwirlFileTabRead,SwirlMdlPointVelCompGet
! Initialization
if (isStart) then ! Program start
 isFileSuccess=SwirlFileTabRead(prog%mdl) ! Read initial data file (CSV)
 if (isFileSuccess) then ! File read successfully
  call SwirlIfcSummaryPrint(prog%mdl) ! Print summary info
 else ! Error reading file
  call TextPrint("") ! Empty line
  call TextPrint("[!] Error reading Swirl parameter file!")
  call TextPrint("[!] Program will retain default inlet parameters!")
  call TextPrint("") ! Empty line
 end if
end if
if (.not.isFileSuccess) return ! Check initial data file reading
allocate(faceDirList(nfabor)) ! Allocate face directing list (faces on inlet boundary)
! Set velocity components
call getfbr(prog%mdl%inletMnem,nFace,faceDirList) ! Get directing list for boundary faces
do iFace=1,nFace ! Inlet face loop
 indFace=faceDirList(iFace) ! Face index in common (mesh) array
 indCell=ifabor(indFace) ! Adjacent element (cell) index in common (mesh) array
 ! Point coordinates (cell center) and face normal vector
 do iCoord=1,nDim ! Coordinate loop
  prog%mdl%point%coord(iCoord)=cdgfbo(iCoord,indFace) ! Coordinate of point at cell center [m]
  prog%mdl%face%nrmVct(iCoord)=-surfbo(iCoord,indFace) ! Component of face normal vector (norm is face area)
 end do
 isSuccess=SwirlMdlPointVelCompGet(prog%mdl) ! Calculate velocity components in stationary coordinates
 if (.not.isSuccess) continue ! Select next face if velocity components not calculated
 indTabNode=prog%mdl%tab%indNodeSel ! Index of selected table row (current radius)
 if (indTabNode.gt.0) then ! Table row selected
  temp=prog%mdl%tab%node(indTabNode)%temp ! Temperature [C]
  turbEng=prog%mdl%tab%node(indTabNode)%turbEng ! Turbulence kinetic energy
  turbDiss=prog%mdl%tab%node(indTabNode)%turbDiss ! Turbulence energy dissipation
 else ! Table row not selected
  temp=0 ! Temperature [C]
  turbEng=0 ! Turbulence kinetic energy
  turbDiss=0 ! Turbulence energy dissipation
 end if
 ! Inlet parameters
 itypfb(indFace)=ientre ! Type of boundary condition (inlet)
 ! Velocity components [m/s]
 rcodcl(indFace,iu,1)=prog%mdl%vel%local%stat%comp(dirX) ! Axis X
 rcodcl(indFace,iv,1)=prog%mdl%vel%local%stat%comp(dirY) ! Axis Y
 rcodcl(indFace,iw,1)=prog%mdl%vel%local%stat%comp(dirZ) ! Axis Z
 if (prog%mdl%tab%indTempCol.gt.0) then ! Temperature defined in parameter table
  if (iscalt.gt.0) then ! Temperature scalar is present in the model
   rcodcl(indFace,isca(iscalt),1)=temp ! Temperature [C]
  end if
 end if  
 if ((prog%mdl%tab%indTurbEngCol.gt.0).and.(prog%mdl%tab%indTurbDissCol.gt.0)) then ! Turbulence parameters defined in table
  ! Turbulence model in Code Saturne identified by variables "iturb" (two digit integer,
  ! model and submodel) and "itytur" (first digit of "iturb" variable)
  turbEngCoef=2.0d0/3.0d0 ! Turbulent kinetic energy coefficient
  select case (itytur) ! Select turbulence model
  case (2) ! K-epsilon models
   rcodcl(indFace,ik,1)=turbEng ! Turbulent kinetic energy [J/kg]
   rcodcl(indFace,iep,1)=turbDiss ! Turbulent energy dissipation rate [J/(kg*s)]
  case (3) ! Reynolds stress models
   ! Turbulent stress components
   ! Normal components
   rcodcl(indFace,ir11,1)=turbEngCoef*turbEng ! Axis X
   rcodcl(indFace,ir22,1)=turbEngCoef*turbEng ! Axis Y
   rcodcl(indFace,ir33,1)=turbEngCoef*turbEng ! Axis Z
   ! Shear components
   rcodcl(indFace,ir12,1)=0.0d0 ! Plane XY
   rcodcl(indFace,ir13,1)=0.0d0 ! Plane XZ
   rcodcl(indFace,ir23,1)=0.0d0 ! Plane YZ
   rcodcl(indFace,iep,1)=turbDiss ! Turbulent energy dissipation rate [J/(kg*s)]
   select case (iturb) ! Select turbulent submodel
   case (32) ! Submodel EBRSM
    rcodcl(indFace,ial,1)=1.0d0 ! Parameter Alpha for EBRSM
   end select
  case (5) ! V2F model
   rcodcl(indFace,ik,1)=turbEng ! Turbulent kinetic energy [J/kg]
   rcodcl(indFace,iep,1)=turbDiss ! Turbulent energy dissipation rate [J/(kg*s)]
   rcodcl(indFace,iphi,1)=turbEngCoef ! Parameter Phi for V2F model
   select case (iturb) ! Select turbulent submodel
   case (50) ! Submodel V2F Phi
    rcodcl(indFace,ifb,1)=0.0d0 ! Parameter f for V2F Phi model
   case (51) ! Submodel V2F BL-V2-K
    rcodcl(indFace,ial,1)=0.0d0 ! Parameter Alpha for V2F BL-V2-K model
   end select
  case (6) ! K-Omega SST Model
   ! Parameter Cmu defined by user or main program
   rcodcl(indFace,ik,1)=turbEng ! Turbulent kinetic energy [J/kg]
   rcodcl(indFace,iomg,1)=turbDiss/cmu/turbEng ! Specific turbulent energy dissipation rate
  case (7) ! Spalart-Allmaras Model
   rcodcl(indFace,inusa,1)=(cmu*turbEng**2)/turbDiss ! Specific parameter for Spalart-Allmaras model
  case default ! Turbulence model not supported in this plugin
   if ((prog%mdl%tab%indTurbEngCol.gt.0).or.&
       (prog%mdl%tab%indTurbDissCol.gt.0)) then ! Turbulence parameters present in table
    if (isStart) call TextPrint("[!] Warning! Turbulence model not supported in Swirl plugin!")
   end if
  end select
 end if
end do
! Free resources
deallocate(faceDirList)  ! Deallocate temporary array for boundary faces
isStart=.false. ! Program start (initialization) indicator
end subroutine cs_f_user_boundary_conditions