c
C**** SUBROUTINE GXKESO is called from group 13 OF GREX3 by setting
C the value ascribed to GROUND in the COVAL statement, and
C is entered only when the patch name begins with the
C character 'KESO'. It is also activated by TURMOD(KLMODL) or
C TURMOD(KEMODL) option.
C
C.... The arguments VIST and LEN1 are the integer names used in the
C SATELLITE, indicating which whole-field store will be used for
C the turbulent kinematic viscosity and length scale of phase 1
C respectively. The argument VISL refers to the laminar kinematic
C viscosity which is required for the low-Reynolds-number extension
C to the k-eps model.
C
C.... The library cases 151,154,174, 191 (k-l model), 152,175,192,193
C (k-e model) make use of it.
C
C.... Following subroutines are used:
C
C fn1(y,a) y = a
C fn7(y,x,a,b,c) y = (a+b*x)/(1.+c*x)
C fn15(y,x1,x2,a,b1) y = (a+b1*x1)/x2
C fn21(y,x1,x2,a,b) y = a+b*x1*x2
C fn26(y,x) y = y*x
C fn27(y,x) y = y/x
C fn28(y,x,a) y = a/x
C fn31(y,x1,x2,a,b1,b2) y = a*(x1**b1)*(x2**b2)
C fn34(y,x,a) y = y+a*x
C fn37(y,x,a) y = y*x**a
C fn56(y,x1,x2,x3,a) y = a*x1*x2/x3
C fn54(y,x1,x2,a) y = y+a*x1/x2
C fn53(y,x1,x2,a) y = y+a*x1*x2
C fn63(y,x,a) y = a*x**4
C fn68(y,x1,x2,x3,a,b) y = (a+b*x1)/(x2*x3)
C fn76(Y,X1,X2) y = y*x1/x2
C fn77(Y,X1,X2) y = y*x1*x2
C fn78(Y,X1,X2,X3,A) y = y+a*x1*x2/x3
C
SUBROUTINE GXKESO(VIST,LEN1,VISL,FIXVAL)
INCLUDE 'farray'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'parear'
COMMON /IGE/IXF,IXL,IYF,IYL,IREG,NZSTEP,IGR,ISC,IRUN,IZSTEP,ITHYD,
1 ISWEEP,ISTEP,INDVAR,VAL,CO,NDIREC,WALDIS,PATGEO,IGES20(6)
COMMON/IDATA/ NX,NY,NZ,LUPR1,IGFILL(116)
1 /RDATA/ TINY,RD1(17),ENUL,RD2(66)
1 /LDEBUG/LDBG1(35),TSTGNK,LDBG2(8)
C 1 /TSKEM/ GKTDKP,LBKP,LBKT,LBET,LBVOSQ,LBOMEG
COMMON/TSKEMI/ LBKP,LBKT,LBET,LBVOSQ,LBOMEG
1 /TSKEMR/ GKTDKP
1 /LRNTM1/L0WDIS,L0FMU,L0FONE,L0FTWO,L0REYN,L0REYT,L0UD1,
1 L0UD2,L0UD3,L0UD4
1 /LRNTM2/LBWDIS,LBFMU,LBFONE,LBFTWO,LBREYN,LBREYT,LBEPKE
1 /MMKKE/ LBFOMG,L0FOMG
COMMON/RDA3/ PRNDTL(150) /RDA4/PRT(150) /RDA6/VARMIN(150)
1 /RDA7/ VARMAX(150) /IDA2/LITER(150)
COMMON/GENI/ NXNY,IGFIL(59)
COMMON/NAMFN/NAMFUN,NAMSUB
COMMON/KEPRES/ SUMPK,SUMPE
LOGICAL LDBG1,TSTGNK,LDBG2,ZRGRN3,SLD,ISACTIVEIJ,ISACTIVECELLS
INTEGER VAL,CO,WALDIS,PATGEO,VIST,VISL,COGRN,VALGRN
CHARACTER*6 NAMFUN,NAMSUB
C
NAMSUB = 'GXKESO'
COGRN=ISC-1 ; VALGRN=ISC-12
C.... Set addresses when additional 3D-storage is provided:
IF(LBFMU.NE.0) L0FMU = L0F(LBFMU)
IF(LBFONE.NE.0) L0FONE= L0F(LBFONE)
IF(LBFTWO.NE.0) L0FTWO= L0F(LBFTWO)
IF(LBREYN.NE.0) L0REYN= L0F(LBREYN)
IF(LBREYT.NE.0) L0REYT= L0F(LBREYT)
IF(LBFOMG.NE.0) L0FOMG= L0F(LBFOMG)
IF(COGRN.EQ.4) THEN
C.... Coefficient part of linearized dissipation-of-turbulent-kinetic-
C energy source
C.... Modification for 2-layer model
IF(INDVAR.EQ.KE) THEN
IF(IENUTA.EQ.3.OR.IENUTA.EQ.4) THEN
CALL GXREYT(-L0REYT,KE,EP,VISL)
CALL GXREYN(-L0REYN,KE,-L0WDIS,VISL)
CALL GXLRDF(1,L0FMU,L0FONE,L0FTWO,L0REYT,L0REYN)
ELSEIF(IENUTA.EQ.8) THEN
CALL GXREYN(-L0REYN,KE,-L0WDIS,VISL)
CALL GXFMU (L0FMU,L0REYN)
CALL GXFEPS(L0FTWO,L0REYN)
ELSEIF(IENUTA.EQ.11) THEN
CALL FN68(-L0REYT,KE,LBOMEG,VISL,0.,1.0)
CALL FN7(-L0FMU,-L0REYT,1./40.,1./6.,1./6.)
CALL FN7(-L0FONE,-L0REYT,0.1,1./2.7,1./2.7)
CALL FN27(-L0FONE,-L0FMU)
CALL FN63(-L0FTWO,-L0REYT,1./4096.)
CALL FN7(-L0FTWO,-L0FTWO,5./18.,1.,1.)
ENDIF
ENDIF
C.... Apply linearization for the sources of KE and EP according to the
C KELIN setting
IF(KELIN.EQ.0) THEN
C>>>>>>>>>>>>>>>>>>>>>>>> KELIN = 0 <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
IF(INDVAR.EQ.KE) THEN
CONST= CD/CMU
ELSEIF(INDVAR.EQ.EP) THEN
CONST= C2E*CD/CMU
ELSEIF(INDVAR.EQ.LBVOSQ) THEN
CONST= 2.0*(C2E-1.0)*CD/CMU
ELSEIF(INDVAR.EQ.LBOMEG) THEN
CONST= C2E/(CMU*CMU)
IF(STORE(EP)) CALL FN21(EP,LBOMEG,KE,0.,CMUCD)
ENDIF
CALL FN31(CO,VIST,LEN1,CONST,1.0,-2.0)
if(tstgnk) then
write(14,*) 'debug from GXKESO for indvar= ',indvar
call prn('vist',vist)
call prn('len1',len1)
call prn('co ',co)
endif
IF(IENUTA.EQ.3.OR.IENUTA.EQ.4) THEN
IF(INDVAR.EQ.KE) THEN
CALL FN27(CO,-L0FMU)
ELSEIF(INDVAR.EQ.EP) THEN
CALL FN76(CO,-L0FTWO,-L0FMU)
ENDIF
C.... 2-layer model & low-re k-omega model
ELSEIF(IENUTA.EQ.8.OR.IENUTA.EQ.11) THEN
IF(INDVAR.EQ.KE) CALL FN76(CO,-L0FTWO,-L0FMU)
IF(INDVAR.EQ.LBOMEG) CALL FN27(CO,-L0FMU)
ELSEIF(IENUTA.EQ.12) THEN
CALL FN27(CO,-L0FOMG)
ENDIF
ELSEIF(KELIN.EQ.1) THEN
C>>>>>>>>>>>>>>>>>>>>>>>> KELIN = 1 <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
IF(INDVAR.EQ.KE) THEN
IF(IENUTA.EQ.3.OR.IENUTA.EQ.4) THEN
CALL FN56(CO,KE,-L0FMU,VIST,C2E*CMUCD)
IF(GEN1.NE.0) CALL FN78(CO,VIST,GEN1,KE,0.5)
C.... 2-layer model
ELSEIF(IENUTA.EQ.8) THEN
CALL FN56(CO,KE,-L0FMU,VIST,C2E*CMUCD)
CALL FN26(CO,-L0FTWO)
IF(GEN1.NE.0) CALL FN78(CO,VIST,GEN1,KE,0.5)
ELSEIF(IENUTA.EQ.12) THEN
CALL FN56(CO,KE,-L0FOMG,VIST,C2E*CMUCD)
IF(GEN1.NE.0) CALL FN78(CO,VIST,GEN1,KE,0.5)
ELSE
IF(GEN1.NE.0) THEN
CALL FN56(CO,VIST,GEN1,KE,0.5)
IF(IENUTA.EQ.13) CALL FN26(CO,-L0FOMG) ! KL K-E turbulence model
ELSE
CALL FN1(CO,0.0)
ENDIF
CALL FN54(CO,KE,VIST,C2E*CMUCD)
ENDIF
CALL FN0(EASP7,CO)
ELSEIF(INDVAR.EQ.EP) THEN
CALL FN28(CO,VIST, (2.0*C2E-1.0)*CMUCD)
CALL FN26(CO,KE)
IF(IENUTA.EQ.3.OR.IENUTA.EQ.4) THEN
CALL FN77(CO,-L0FTWO,-L0FMU)
ELSEIF(IENUTA.EQ.12) THEN
CALL FN26(CO,-L0FOMG)
ENDIF
ENDIF
ELSEIF(KELIN.EQ.2) THEN
C>>>>>>>>>>>>>>>>>>>>>>>> KELIN = 2 <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
IF(INDVAR.EQ.KE) THEN
CALL FN15(CO,KE,VIST,0.0,CMUCD)
C.... 2-layer model
IF(IENUTA.EQ.3.OR.IENUTA.EQ.4) THEN
CALL FN26(CO,-L0FMU)
ELSEIF(IENUTA.EQ.8) THEN
CALL FN77(CO,-L0FMU,-L0FTWO)
ELSEIF(IENUTA.EQ.12) THEN
CALL FN26(CO,-L0FOMG)
ENDIF
ELSEIF(INDVAR.EQ.EP) THEN
CALL FN15(CO,KE,VIST,0.0,C2E*CMUCD)
IF(IENUTA.EQ.3.OR.IENUTA.EQ.4) THEN
CALL FN77(CO,-L0FMU,-L0FTWO)
C... MMK high-Reynolds-number K-E turbulence model
ELSEIF(IENUTA.EQ.12) THEN
CALL FN26(CO,-L0FOMG)
ENDIF
ENDIF
ELSEIF(KELIN.EQ.3) THEN
C>>>>>>>>>>>>>>>>>>>>>>>> KELIN = 3 <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
L0EPKE= L0F(LBEPKE)
C.... Set LITER=1 as sign that a change has been made
LITER(LBEPKE)= 1
L0CO= L0F(CO)
IF(INDVAR.EQ.KE) THEN
L0KE= L0F(KE)
FACTOR= 1.5
IF(SOLVE(EP)) THEN
L0EP= L0F(EP)
DO 300 I=1,NXNY
IF(SLD(I)) THEN
F(L0EPKE+I)=0.0
ELSE
EPKE= F(L0EP+I)/(F(L0KE+I)+TINY)
F(L0EPKE+I)= AMAX1(VARMIN(LBEPKE),
1 AMIN1(VARMAX(LBEPKE),EPKE))
ENDIF
300 CONTINUE
ELSE
L0LEN= L0F(LEN1)
DO 301 I=1,NXNY
IF(SLD(I)) THEN
F(L0EPKE+I)=0.0
ELSE
EPKE= CD*F(L0KE+I)**0.5/(F(L0LEN+I)+TINY)
F(L0EPKE+I)= AMAX1(VARMIN(LBEPKE),
1 AMIN1(VARMAX(LBEPKE),EPKE))
ENDIF
301 CONTINUE
ENDIF
ELSEIF(INDVAR.EQ.EP) THEN
FACTOR= 1.3333*C2E
ENDIF
C.... Co = 0.0 + EP/KE * FACTOR
CALL FN2(CO,LBEPKE,0.0,FACTOR)
cc call prn('co ',co)
C.... KELIN= 3 for low-Re
IF((IENUTA.EQ.3.OR.IENUTA.EQ.4.OR.IENUTA.EQ.8) .AND.
1 (INDVAR.EQ.EP)) CALL FN26(CO,-L0FTWO)
ENDIF
C.... 2-layer model
IF(INDVAR.EQ.EP .AND. IENUTA.EQ.8) THEN
L0CO= L0F(CO)
CONST= 0.8*FIXVAL
DO 101 J=1,NXNY
IF(F(L0REYN+J).LT.350.0) F(L0CO+J)= CONST
101 CONTINUE
ENDIF
ELSEIF(VALGRN.EQ.4) THEN
C.... Value part of linearized dissipation-of-turbulent-kinetic-energy
C source
IF(KELIN.EQ.0) THEN
C>>>>>>>>>>>>>>>>>>>>>>>> KELIN = 0 <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
IF(INDVAR.EQ.KE) THEN
CONST= CMU/CD
IF(GEN1.NE.0) THEN
CALL FN31(VAL,GEN1,LEN1,CONST,1.0,2.0)
ELSE
CALL FN1(VAL,0.0)
ENDIF
IF(IENUTA.EQ.3.OR.IENUTA.EQ.4) CALL FN26(VAL,-L0FMU)
C.... 2-layer K-E model & low-re k-omega model
IF(IENUTA.EQ.8.OR.IENUTA.EQ.11)
1 CALL FN76(VAL,-L0FMU,-L0FTWO)
IF(IENUTA.EQ.12.OR.IENUTA.EQ.13) CALL FN26(VAL,-L0FOMG)
ELSEIF(INDVAR.EQ.EP) THEN
IF(GEN1.NE.0) THEN
CALL FN21(VAL,GEN1,VIST,0.0,C1E/C2E)
ELSE
CALL FN1(VAL,0.0)
ENDIF
IF(IENUTA.EQ.3.OR.IENUTA.EQ.4)
1 CALL FN76(VAL,-L0FONE,-L0FTWO)
IF(IENUTA.EQ.13) CALL FN26(VAL,-L0FOMG)
ELSEIF(INDVAR.EQ.LBVOSQ) THEN
IF(GEN1.NE.0) THEN
CALL FN31(VAL,GEN1,LEN1,CMU/CD,1.0,2.0)
ELSE
CALL FN1(VAL,0.0)
ENDIF
CALL FN26(VAL,LBVOSQ)
CALL FN27(VAL,KE)
CALL FN25(VAL,(CMU/CD)*((C1E-1.0)/(C2E-1.0))**2)
ELSEIF(INDVAR.EQ.LBOMEG) THEN
IF(GEN1.NE.0) THEN
CALL FN31(VAL,GEN1,LEN1,C1E/C2E,1.0,2.0)
ELSE
CALL FN1(VAL,0.0)
ENDIF
CALL FN26(VAL,LBOMEG)
CALL FN27(VAL,KE)
CALL FN25(VAL,CMU*CMU)
IF(IENUTA.EQ.11) CALL FN77(VAL,-L0FMU,-L0FONE)
ENDIF
ELSEIF(KELIN.EQ.1) THEN
C>>>>>>>>>>>>>>>>>>>>>>>> KELIN = 1 <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
IF(INDVAR.EQ.KE) THEN
CALL FN56(VAL,KE,KE,VIST, (C2E-1.0)*CMUCD)
IF(IENUTA.EQ.3.OR.IENUTA.EQ.4) THEN
CALL FN26(VAL,-L0FMU)
ELSEIF(IENUTA.EQ.8) THEN
C.... 2-layer K-E model
CALL FN77(VAL,-L0FMU,-L0FTWO)
ELSEIF(IENUTA.EQ.12) THEN
CALL FN26(VAL,-L0FOMG)
ENDIF
IF(GEN1.NE.0) CALL FN53(VAL,VIST,GEN1,1.5)
CALL FN27(VAL,EASP7)
ELSEIF(INDVAR.EQ.EP) THEN
IF(GEN1.NE.0) THEN
CALL FN21(VAL,GEN1,VIST,0.0,C1E/(2.0*C2E-1.0))
ELSE
CALL FN1(VAL,0.0)
ENDIF
IF(IENUTA.EQ.3 .OR. IENUTA.EQ.4)
1 CALL FN76(VAL,-L0FONE,-L0FTWO)
CALL FN34(VAL,EP, (C2E-1.0)/(2.0*C2E-1.0))
ENDIF
ELSEIF(KELIN.EQ.2) THEN
C>>>>>>>>>>>>>>>>>>>>>>>> KELIN = 2 <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
IF(INDVAR.EQ.KE) THEN
IF(GEN1.NE.0) THEN
CALL FN15(VAL,GEN1,KE,0.0,1.0/CMUCD)
ELSE
CALL FN1(VAL,0.0)
ENDIF
CALL FN37(VAL,VIST,2.0)
IF(IENUTA.EQ.3.OR.IENUTA.EQ.4) THEN
C.... Low-Re K-E models
CALL FN27(VAL,-L0FMU)
ELSEIF(IENUTA.EQ.8) THEN
C.... 2-layer K-E model
CALL FN27(VAL,-L0FMU)
CALL FN27(VAL,-L0FTWO)
ELSEIF(IENUTA.EQ.12) THEN
CALL FN27(VAL,-L0FOMG)
ENDIF
ELSEIF(INDVAR.EQ.EP) THEN
IF(GEN1.NE.0) THEN
CALL FN21(VAL,VIST,GEN1,0.0,C1E/C2E)
ELSE
CALL FN1(VAL,0.0)
ENDIF
IF(IENUTA.EQ.3.OR.IENUTA.EQ.4)
1 CALL FN76(VAL,-L0FONE,-L0FTWO)
ENDIF
ELSEIF(KELIN.EQ.3) THEN
C>>>>>>>>>>>>>>>>>>>>>>>> KELIN = 3 <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
L0EPKE= L0F(LBEPKE); L0VAR = L0F(INDVAR)
L0VAL = L0F(VAL); L0CO=L0F(CO)
L0SU = L0F(LSU); L0AP=L0F(LAP)
IF(GEN1.NE.0) L0GEN= L0F(GEN1)
L0VIST= L0F(VIST)
L0EA1 = L0F(EASP1)
L0M1 = L0F(LM1)
IF(INDVAR.EQ.KE) THEN
FACTOR= 0.333
IF(IZ==1) SUMPK=0.0
C.... Note that VAL is used only for the EP part of the source
C the generation part is put directly into SU.
IF(IENUTA.EQ.3.OR.IENUTA.EQ.4) THEN
C.... Low-Re K-E models
IF(GEN1.EQ.0) THEN
DO 312 I= 1,NXNY
312 F(L0VAL+I)= FACTOR*F(L0VAR+I)
ELSE
DO 302 I= 1,NXNY
IF(F(L0AP+I)>FIXVAL) THEN
F(L0VAL+I)=0.; F(L0CO+I)=0.0; SOURCE=0.0
ELSE
F(L0VAL+I)= FACTOR*F(L0VAR+I)
SOURCE = F(L0VIST+I)*F(L0GEN+I)*F(L0M1+I)
F(L0SU+I) = F(L0SU+I) + SOURCE
IF(ISACTIVEIJ(I,IZ)) SUMPK=SUMPK+SOURCE
ENDIF
302 CONTINUE
ENDIF
ELSE
IF(GEN1.EQ.0) THEN
DO 313 IX= IXF,IXL
IADX= (IX-1)*NY
DO 313 IY= IYF,IYL
IJ= IY+IADX
IF(ZRGRN3(IX,IY,IZSTEP,IJ)) THEN
F(L0VAL+IJ)= F(L0VAR+IJ)
ELSE
F(L0VAL+IJ)= FACTOR*F(L0VAR+IJ)
ENDIF
313 CONTINUE
ELSE
DO 303 IX= IXF,IXL
IADX= (IX-1)*NY
DO 303 IY= IYF,IYL
IJ= IY+IADX
IF(ZRGRN3(IX,IY,IZSTEP,IJ)) THEN
F(L0VAL+IJ)= F(L0VAR+IJ)
ELSE
IF(F(L0AP+IJ)>FIXVAL) THEN
F(L0VAL+IJ)=0.; F(L0CO+IJ)=0; SOURCE=0.
ELSE
F(L0VAL+IJ)= FACTOR*F(L0VAR+IJ)
SOURCE = F(L0VIST+IJ)*F(L0GEN+IJ)*F(L0M1+IJ)
IF(IENUTA.EQ.13) SOURCE = SOURCE*F(L0FOMG+IJ)
F(L0SU+IJ) = F(L0SU+IJ) + SOURCE
IF(ISACTIVECELLS(IX,IY,IZ)) SUMPK=SUMPK+SOURCE
ENDIF
ENDIF
303 CONTINUE
ENDIF
ENDIF
IF(IZ==NZ.AND.MIMD.AND.NPROC>1) CALL GLSUM(SUMPK)
ELSEIF(INDVAR.EQ.EP) THEN
FACTOR= 0.25
IF(IZ==1) SUMPE=0.0
IF(IENUTA.EQ.3.OR.IENUTA.EQ.4) THEN
C.... Low-Re K-E models
IF(GEN1.EQ.0) THEN
DO 314 I= 1,NXNY
314 F(L0VAL+I)= FACTOR*F(L0VAR+I)
ELSE
DO 304 I= 1,NXNY
IF(F(L0AP+I)>FIXVAL) THEN
F(L0VAL+I)=0.; F(L0CO+I)=0.; SOURCE=0.
ELSE
F(L0VAL+I)= FACTOR*F(L0VAR+I)
SOURCE = F(L0VIST+I)*F(L0GEN+I)*F(L0M1+I)
F(L0SU+I)=F(L0SU+I)+C1E*F(L0FONE+I)*F(L0EPKE+I)*
1 SOURCE
IF(ISACTIVEIJ(I,IZ)) SUMPE=SUMPE+
1 C1E*F(L0FONE+I)*F(L0EPKE+I)*SOURCE
ENDIF
304 CONTINUE
ENDIF
ELSE
IF(GEN1.EQ.0) THEN
DO 315 I= 1,NXNY
315 F(L0VAL+I)= FACTOR*F(L0VAR+I)
ELSE
DO 305 I= 1,NXNY
IF(F(L0AP+I)>FIXVAL) THEN
F(L0VAL+I)=0.; F(L0CO+I)=0.; SOURCE=0.
ELSE
F(L0VAL+I)= FACTOR*F(L0VAR+I)
SOURCE = F(L0VIST+I)*F(L0GEN+I)*F(L0M1+I)
IF(IENUTA.EQ.13) SOURCE = SOURCE*F(L0FOMG+I)
F(L0SU+I) = F(L0SU+I) + C1E*F(L0EPKE+I)*SOURCE
IF(ISACTIVEIJ(I,IZ)) SUMPE=SUMPE+
1 C1E*F(L0EPKE+I)*SOURCE
ENDIF
305 CONTINUE
ENDIF
ENDIF
IF(IZ==NZ.AND.MIMD.AND.NPROC>1) CALL GLSUM(SUMPE)
ENDIF
ENDIF
C.... 2-layer K-E model
IF(INDVAR.EQ.EP .AND. IENUTA.EQ.8) THEN
L0KE=L0F(KE) ; L0LEN=L0F(LEN1) ; L0VAL=L0F(VAL)
DO 201 I= 1,NXNY
IF(F(L0REYN+I).LT.350.0) F(L0VAL+I)=
1 CD*ABS(F(L0KE+I))**1.5*F(L0FTWO+I)/(F(L0LEN+I)+TINY)
201 CONTINUE
ENDIF
C.... For high-Re models protect wall-function settings:
IF(.NOT.(IENUTA.EQ.3 .OR. IENUTA.EQ.4)) THEN
IF(INDVAR.EQ.KE .AND. KELIN.NE.3) THEN
L0VAR=L0F(INDVAR) ; L0VAL=L0F(VAL)
C.... Do not modify KE source for high-Re models near walls with GRND3
C wall-function
DO 202 IX= IXF,IXL
IADX= (IX-1)*NY
DO 202 IY= IYF,IYL
IJ= IY+IADX
IF(ZRGRN3(IX,IY,IZSTEP,IJ)) F(L0VAL+IJ)= F(L0VAR+IJ)
202 CONTINUE
ENDIF
ENDIF
ENDIF
NAMSUB = 'gxkeso'
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
C**** SUBROUTINE GXLESO is called from group 13, GREX3, by
C setting the value ascribed to GROUND in the COVAL statement,
C and is entered only when the patch name begins with the
C character 'LESO'.
C
C.... The dummy DIREC is the direction-index, which is set to 'SOUTH'
C in GREX3.
C
C.... The library case 975 makes use of it.
C
SUBROUTINE GXLESO(DIREC)
INCLUDE 'farray'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
COMMON /IGE/IXF,IXL,IYF,IYL,IREG,NZSTEP,IGR,ISC,IRUN,IZSTEP,ITHYD,
1 ISWEEP,ISTEP,INDVAR,VAL,CO,NDIREC,WALDIS,PATGEO,IGES20(6)
INTEGER VAL,CO,WALDIS,PATGEO
CHARACTER*(*) DIREC
COMMON /NAMFN/NAMFUN,NAMSUB
CHARACTER*6 NAMFUN,NAMSUB
C
NAMSUB = 'GXLESO'
IF(ISC.EQ.13) THEN
C.... Value set to + or - const * INTFRC * mean abs(dwdy)
C for the two-fluid model of turbulence.
IF(MOD(INDVAR,2).EQ.0) THEN
C fn10(y,x1,x2,a,b1,b2) y = a+b1*x1+b2*x2
C fn40(y) y = abs(y)
C fnav(y,x,direc) This subroutine puts the average of x and
C X(DIREC) into Y, where DIREC is 'NORTH','SOUTH','EAST',
C 'WEST','HIGH' or 'LOW'.
CALL FN10(GEN1,V1,V2,0.0,ELSOA,-ELSOA)
CALL FN40(GEN1)
CALL FNAV(VAL,GEN1,DIREC)
ELSE
CALL FNAV(VAL,GEN1,DIREC)
ENDIF
ENDIF
NAMSUB = 'gxleso'
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
C**** SUBROUTINE GXSHSO is called from group 13 of GREX3, by setting
C values ascribed to GROUND in the COVAL statement in SATELLITE.
C The subroutine is entered only when the patch name begins with
C the characters 'SHSO'.
C
C.... The dummy INTFRC and LEN1 are the integer names used in the
C SATELLITE, indicating which whole-field store will be used
C for the inter-phase friction coefficient and the length-scale
C of phase 1 in response to the command STORE(CFIPS) and STORE(EL1)
C respectively.
C
C.... The library case 975 (VAL=GRND5) exemplifies the use of the
C subroutine.
C
SUBROUTINE GXSHSO(INTFRC,LEN1)
INCLUDE 'farray'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
COMMON /IGE/IXF,IXL,IYF,IYL,IREG,NZSTEP,IGR,ISC,IRUN,IZSTEP,ITHYD,
1 ISWEEP,ISTEP,INDVAR,VAL,CO,NDIREC,WALDIS,PATGEO,IGES20(6)
INTEGER VAL,CO,WALDIS,PATGEO
COMMON /NAMFN/NAMFUN,NAMSUB
CHARACTER*6 NAMFUN,NAMSUB
C
NAMSUB = 'GXSHSO'
IF(ISC.EQ.17) THEN
IF(INDVAR.EQ.V1) THEN
C.... Multiply interphase-friction coefficient by length scale
C fn21(y,x1,x2,a,b) y = a + b*x1*x2
CALL FN21(EASP1,INTFRC,LEN1,0.0,1.0)
C.... Get average value at v location
CALL FNAV(EASP1,EASP1,'NORTH')
C.... Arithmetic-mean w-gradient times SHSOA
C fn10(y,x1,x2,a,b1,b2) y = a+b1*x1+b2*x2
CALL FN10(VAL,W1,W2,0.0,0.5*SHSOA,0.5*SHSOA)
C fn103(y,x,idir) This subroutine puts x(next)-x into y, where
C 'next' is the next value in the north, south,
C east or west direction, indicated by
C IDIR=1,2,3 or 4.
CALL FN103(VAL,VAL,1)
C fn56(y,x1,x2,x3,a) y = a*x1*x2/x3
CALL FN56(VAL,VAL,EASP1,DYG2D,1.0)
C fn40(y) y = abs(y)
CALL FN40(VAL)
C.... Store for use for second-phase value
C fn0(y,x) y = x
CALL FN0(EASP1,VAL)
C fn2(y,x,a,b) y = a+b*x
ELSEIF(INDVAR.EQ.V2) THEN
CALL FN2(VAL,EASP1,0.0,-1.0)
ENDIF
ENDIF
NAMSUB = 'gxshso'
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C
LOGICAL FUNCTION ZRGRN3(IX,IY,IZ,IJ)
INCLUDE 'farray'
INCLUDE 'patnos'
COMMON /GENI/NXNY,IGFILL(59)
COMMON/IDATA/NX,NY,NZ,LUPR1,IDFIL1(116)
COMMON /INTDMN/IDMN,IDUM(8)
COMMON /LRNTM3/L0UTAU,NMWALL,L0WALL,L0DSKN,IVPRST
LOGICAL WALPRE,GRNDNO
ZRGRN3= .FALSE.
IF(NMWALL.GT.0) THEN
IF(WALPRE(IJ)) THEN
C.... Wall-type patch is present in the cell, so loop over wall-patches
C to check for GRND3 coefficient for KE
DO 10 IWL= 1,NMWALL
IPT= NINT(F(L0WALL+IWL))
IF(IPT.LT.0) THEN
CALL GETIZS(-IPT,IZF,IZL)
IF(IZF.LE.IZ .AND. IZ.LE.IZL) THEN
CALL GETIXS(-IPT,IXF,IXL)
IF(IXF.LE.IX .AND. IX.LE.IXL) THEN
CALL GETIYS(-IPT,IYF,IYL)
IF(IYF.LE.IY .AND. IY.LE.IYL) THEN
ZRGRN3= .TRUE.
RETURN
ENDIF
ENDIF
ENDIF
ENDIF
10 CONTINUE
ENDIF
ENDIF
C.... Check surrounding cells to see if any EGWF GRND3 cells
DO IFACE=1,6
IZZ=IZ; IJJ=IJ
IF(IFACE==1) THEN ! check to North
IF(IY==NY) CYCLE
IJJ=IJ+1
ELSEIF(IFACE==2) THEN ! South
IF(IY==1) CYCLE
IJJ=IJ-1
ELSEIF(IFACE==3) THEN ! East
IF(IX==NX) CYCLE
IJJ=IJ+NY
ELSEIF(IFACE==4) THEN ! West
IF(IX==1) CYCLE
IJJ=IJ-NY
ELSEIF(IFACE==5) THEN ! High
IF(IZ==NZ) CYCLE
IZZ=IZ+1
ELSEIF(IFACE==6) THEN ! Low
IF(IZ==1) CYCLE
IZZ=IZ-1
ENDIF
L0PAT=L0PATNO(IDMN)+(IZZ-1)*NXNY
IPAT=NINT(F(L0PAT+IJJ)) ! retrieve number of parent patch
IF(IPAT>0) THEN
GWALLCO=F(L0PATWLCO(IDMN)+IPAT)
IF(GRNDNO(3,GWALLCO)) THEN
ZRGRN3=.TRUE.
RETURN
ENDIF
ENDIF
ENDDO
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c