cC.... File Name ...... GXGENK.FOR........... 120416
C File includes the following subroutines and functions:
C GXGENF, GXSQUR, GFTRM1, GFTRM2, GFDUDX, GFDUDY, GFDUDZ, GFDVDX,
C GFDVDY, GFDVDZ, GFDWDX, GFDWDY, GFDWDZ, AVRGVL, GXKEGB, AXIBFC.
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C.... GXGENF calculates derivatives of flow velocities and Generation
C for a current slab. It is called twice; from Gr.1 Sec.2 to
C define auxiliary variables; and from Gr.19 Sec.4 to make the
C calculations. If IGEN0=0 when GXGENF is called from Gr.19 Sec.4,
C it calculates and stores only DUDX,...,DWDZ. This call is normally
C done for Reynolds stress turbulence model.
C
C.... The dummy IPH is the phase index, which indicates that the Generation
C is calculated for the first-phase (IPH=1) or the second-phase (IPH=2)
C The IGEN0 is the slab-wise address of the Generation storage.
C
C.... There is one /LSG/-logicals temporary in use:
C LSG53=T activates the use of cell-centered averaged velocities
C for both 'in-line' and 'off-line' derivatives.
C.... To check calculations of the generation function the user may
C activate DEBUG 3D-storage of averaged velocities, or velocity
C derivatives by storing appropriate variables in Q1 (see details
C below).
c
SUBROUTINE GXGENF(IPH,IGEN0)
INCLUDE 'farray'
INCLUDE 'satear'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'grdear'
COMMON /CMNGN/ IUG0(2),IVG0(2),IWG0(2),IUIZM0(2),IUIZP0(2),
1 IVIZM0(2),IVIZP0(2),IWIZM0(2),IWIZP0(2),
1 IDUDI0,IDUDJ0,IDUDK0,IDVDI0,IDVDJ0,IDVDK0,
1 IDWDI0,IDWDJ0,IDWDK0
1 /GENI/ IG1(2),NXNYST,IG2(6),NFM,IG3(23),ILTLS,IG4(26)
1 /GENL/ LGL1(14),debgtz,LGL2(43),STRGN1,STRGN2
1 /NAMFN/NAMFUN,NAMSUB
1 /BFICRT/IUCRT(2),IVCRT(2),IWCRT(2),ITMP(6)
1 /GENFL/ LDUDX,LDUDY,LDUDZ,LDVDX,LDVDY,LDVDZ,LDWDX,LDWDY,
1 LDWDZ,COLVEL,AXISYM /FNDRL/NXNE1,NYNE1,NZNE1
1 /LBDFDL/IDUDX,IDUDY,IDUDZ,IDVDX,IDVDY,IDVDZ,IDWDX,IDWDY,
1 IDWDZ,IDSDX,IDSDY,IDSDZ,IDU2X,IDU2Y,IDU2Z,IDV2X,
1 IDV2Y,IDV2Z,IDW2X,IDW2Y,IDW2Z
1 /RSTM2/ J0DUDX,J0DUDY,J0DUDZ,J0DVDX,J0DVDY,J0DVDZ,J0DWDX,
1 J0DWDY,J0DWDZ,JRSTF(45)
1 /F0/ LB1(109),KZXCY,LB2(194)
1 /LSG/ LSGF1(52),LSG53,LSGF2(47)
LOGICAL LSGF1,LSG53,LSGF2
LOGICAL NXNE1,NYNE1,NZNE1,COLVEL,AXISYM,LDUDX,LDUDY,LDUDZ,LDVDX,
1 LDVDY,LDVDZ,LDWDX,LDWDY,LDWDZ,LSTOU,LSTOV,LSTOW,AVEVEL,
1 DVDL3D,LDOGEN,LGL1,STRGN1,STRGN2,LDBAVR,dbgloc,LGL2,debgtz
CHARACTER*6 NAMFUN,NAMSUB
DIMENSION IU1C(2),IV1C(2),IW1C(2)
SAVE ITMP10,ITMP20,ITMP30,IU1C,IV1C,IW1C
SAVE LSTOU,LSTOV,LSTOW,AVEVEL,DVDL3D,LDBAVR
C
NAMSUB = 'GXGENF'
if(indvar>0) then
dbgloc=debgtz.and.dbgphi(indvar)
else
dbgloc=.false.
endif
if(flag.or.dbgloc) then
call banner(1,namsub,200918)
call writ2i('igr ',igr,'isc ',isc)
endif
IF(IGR==1 .AND. ISC==2) THEN ! Call from Gr.1; Sec.2 to prepare
AVEVEL=LSG53; AXISYM=.NOT.CARTES
IF(BFC .AND. NX==1) CALL AXIBFC(AXISYM) ! axisymmetry
! (if Y-Z grid is axisymmetric, additional term V/R appears as
! part of dU/dX):
C.... Define velocities to use for Generation calculation:
IF(BFC) THEN
IUG0(1)=IUCRT(1); IVG0(1)=IVCRT(1); IWG0(1)=IWCRT(1)
IUG0(2)=IUCRT(2); IVG0(2)=IVCRT(2); IWG0(2)=IWCRT(2)
C.... Temporary slab-storage for DUDI,DUDJ,...,DWDK:
IDUDI0=L0F(LD4); IDVDI0=L0F(LD7); IDWDI0=L0F(LD10)
IDUDJ0=L0F(LD5); IDVDJ0=L0F(LD8); IDWDJ0=L0F(LD11)
IDUDK0=L0F(LD6); IDVDK0=L0F(LD9); IDWDK0=L0F(LDSTAG)
ELSEIF(CCM) THEN
IUG0(1)=LBNAME('UC1'); IUG0(2)=LBNAME('UC2')
IVG0(1)=LBNAME('VC1'); IVG0(2)=LBNAME('VC2')
IWG0(1)=LBNAME('WC1'); IWG0(2)=LBNAME('WC2')
ELSE
C.... For staggered solver on Cartesian geometry there are two options:
C 1. AVEVEL=F (default) Use face-centered velocities U1,...,W1 for
C the 'in-line' derivatives (dU/dX,...) and averaged
C cell-centered velocities for the 'off-line' (dU/dY..).
C 2. AVEVEL=T Use averaged cell-centered velocities for the both
C 'in-line' and 'off-line' derivatives.
C Storage for averaged velocities at the current and high slabs is
C allocated in EARTH. LD4,LD5 and LD6 are used as temporary storage
C for low slab.
LSTOU=STORE(3); LSTOV=STORE(5); LSTOW=STORE(7)
IUG0(1)=IUIZM0(1); IVG0(1)=IVIZM0(1); IWG0(1)=IWIZM0(1)
IUIZM0(1)=L0F(LD4); IVIZM0(1)=L0F(LD5); IWIZM0(1)=L0F(LD6)
C.... Debug 3D-storage of averaged velocities:
IU1C(1)= LBNAME('U1C'); IV1C(1)= LBNAME('V1C')
IW1C(1)= LBNAME('W1C')
IU1C(2)= LBNAME('U2C'); IV1C(2)= LBNAME('V2C')
IW1C(2)= LBNAME('W2C')
LDBAVR= IU1C(1)/=0 .OR. IV1C(1)/=0 .OR. IW1C(1)/=0 .OR.
1 IU1C(2)/=0 .OR. IV1C(2)/=0 .OR. IW1C(2)/=0
IUG0(2)=0; IVG0(2)=0; IWG0(2)=0
IF(.NOT.ONEPHS) THEN
IUG0(2)=IUIZM0(2); IVG0(2)=IVIZM0(2); IWG0(2)=IWIZM0(2)
ENDIF
ENDIF
IF(BFC.AND.KZXCY/=0) CALL XCYVEL(IGR,IZ,IUSL0,IVSL0,IWSL0)
C.... COLVEL indicates that velocities are cell-centered
COLVEL= BFC .OR. CCM
C.... Find zero initial addresses for colocated velocities:
IF(COLVEL) THEN
IF(IUG0(1)/=0) IUG0(1)= L0F(IUG0(1))
IF(IVG0(1)/=0) IVG0(1)= L0F(IVG0(1))
IF(IWG0(1)/=0) IWG0(1)= L0F(IWG0(1))
IF(IUG0(2)/=0) IUG0(2)= L0F(IUG0(2))
IF(IVG0(2)/=0) IVG0(2)= L0F(IVG0(2))
IF(IWG0(2)/=0) IWG0(2)= L0F(IWG0(2))
LSTOU= IUG0(1)/=0
LSTOV= IVG0(1)/=0
LSTOW= IWG0(1)/=0
ENDIF
C.... Define auxiliary logicals:
LDUDX=.FALSE.; LDUDY=.FALSE.; LDUDZ=.FALSE.
LDVDX=.FALSE.; LDVDY=.FALSE.; LDVDZ=.FALSE.
LDWDX=.FALSE.; LDWDY=.FALSE.; LDWDZ=.FALSE.
IF(.NOT.PARAB) THEN
LDUDX= LSTOU .AND. (DUDX .OR. GENK .OR. RSTM) .AND. NXNE1
LDUDX= LDUDX .OR. (AXISYM .AND. LSTOV)
LDUDY= LSTOU .AND. (DUDY .OR. GENK .OR. RSTM) .AND. NYNE1
LDUDZ= LSTOU .AND. (DUDZ .OR. GENK .OR. RSTM) .AND. NZNE1
LDVDX= LSTOV .AND. (DVDX .OR. GENK .OR. RSTM) .AND. NXNE1
LDVDY= LSTOV .AND. (DVDY .OR. GENK .OR. RSTM) .AND. NYNE1
LDVDZ= LSTOV .AND. (DVDZ .OR. GENK .OR. RSTM) .AND. NZNE1
LDWDX= LSTOW .AND. (DWDX .OR. GENK .OR. RSTM) .AND. NXNE1
LDWDY= LSTOW .AND. (DWDY .OR. GENK .OR. RSTM) .AND. NYNE1
LDWDZ= LSTOW .AND. (DWDZ .OR. GENK .OR. RSTM) .AND. NZNE1
ELSE
LDWDX= LSTOW .AND. (DWDX .OR. GENK .OR. RSTM) .AND. NXNE1
LDWDY= LSTOW .AND. (DWDY .OR. GENK .OR. RSTM) .AND. NYNE1
ENDIF
C.... Define addresses for temporary storage (LD1,LD2,LD3 are in use):
ITMP10=L0F(LD1); ITMP20=L0F(LD2); ITMP30=L0F(LD3)
C.... Store velocity derivatives for the dump into RESULT file:
DVDL3D= IDUDX/=0 .OR. IDUDY/=0 .OR. IDUDZ/=0 .OR.
1 IDVDX/=0 .OR. IDVDY/=0 .OR. IDVDZ/=0 .OR.
1 IDWDX/=0 .OR. IDWDY/=0 .OR. IDWDZ/=0 .OR.
1 IDU2X/=0 .OR. IDU2Y/=0 .OR. IDU2Z/=0 .OR.
1 IDV2X/=0 .OR. IDV2Y/=0 .OR. IDV2Z/=0 .OR.
1 IDW2X/=0 .OR. IDW2Y/=0 .OR. IDW2Z/=0
C.... For Reynolds stress model DUDX,...,DWDZ are put into the storage
C provided in UST191 subroutine:
DVDL3D= DVDL3D .AND. .NOT.RSTM
ELSE ! Call from Gr.19; Sec.4 to calculate Generation or DUDX,...,etc.
C-----------------------------------------------------------------------
C.... If GEN1 (or GEN2) is in 3D-storage don't calculate it again when
C Gr.19 Sec.4 is entered for variables solved whole-field:
IF( (IPH==1.AND.STRGN1 .OR. IPH==2.AND.STRGN2) .AND.
1 INDVAR>10 .AND. INDVAR/=ILTLS ) RETURN
C.... If IGEN0/=0 calculate generation term, otherwise only DUDX,...
LDOGEN= IGEN0/=0
IF(GCV) THEN
CALL BFGENB(IGEN0)
RETURN
ENDIF
C.... Calculate Generation for MB-FGE technique
IF(CCM .AND. NUMBLK>0) THEN
CALL UNGENB(IPH-1,IGEN0)
RETURN
ENDIF
C.... Recalculate geometric coefficients for moving BFC-grids:
IF(MOVBFC) CALL IJKXYZ(IZSTEP)
IF(.NOT.COLVEL) THEN
C.... Average staggered velocities:
IF(LSTOU) THEN
IF(NXNE1 .AND. (LDUDX.OR.LDUDY.OR.LDUDZ)) THEN
IF(LDUDZ) THEN
IF(IZSTEP>1 ) THEN
CALL SHINX2(IUIZM0(IPH),IUG0(IPH),IUG0(IPH),IUIZP0(IPH))
ELSE
CALL AVRGVL(4,IUG0(IPH),IPH,IZSTEP,0)
ENDIF
IF(IZSTEP1) THEN
CALL SHINX2(IVIZM0(IPH),IVG0(IPH),IVG0(IPH),IVIZP0(IPH))
ELSE
CALL AVRGVL(2,IVG0(IPH),IPH,IZSTEP,0)
ENDIF
IF(IZSTEP1 ) THEN
CALL SHINX2(IWIZM0(IPH),IWG0(IPH),IWG0(IPH),IWIZP0(IPH))
ELSE
c w average for current slab
CALL AVRGVL(6,IWG0(IPH),IPH,IZSTEP,0)
ENDIF
c w average for high slab
IF(IZSTEP
SUBROUTINE GFTRM1(IGEN0,A1,ID10,A2,ID20,A3,ID30)
INCLUDE 'farray'
COMMON /GENI/IG1(2),NXNYST,IG2(57)
I1M0= ID10-IGEN0
I2M0= ID20-IGEN0
I3M0= ID30-IGEN0
DO I= IGEN0+1,IGEN0+NXNYST
F(I)= 2.*(A1*F(I1M0+I)**2 + A2*F(I2M0+I)**2 + A3*F(I3M0+I)**2)
ENDDO
END
C-----------------------------------------------------------------------
c
SUBROUTINE GFTRM2(IGEN0,A1,ID10,A2,ID20)
INCLUDE 'farray'
COMMON /GENI/IG1(2),NXNYST,IG2(57)
LOGICAL QGE
I1M0= ID10-IGEN0
I2M0= ID20-IGEN0
DO I= IGEN0+1,IGEN0+NXNYST
IF(QGE(ABS(A1*F(I1M0+I)),1.E10).OR.QGE(ABS(A2*F(I2M0+I)),1.E10))
1 THEN
F(I)=1.E10
ELSE
F(I)= F(I) + ( A1*F(I1M0+I) + A2*F(I2M0+I) )**2
ENDIF
ENDDO
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C
c
SUBROUTINE GFDUDX(IPH,IZZ,IDUDX0,AVEVEL)
INCLUDE 'cmndmn'
INCLUDE 'farray'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'd_earth/pbcstr'
COMMON /CMNGN/IUG0(2),IVG0(2),IWG0(2),IC1(12),IDUDI0,IDUDJ0,
` IDUDK0,IC2(6) /GENFL/LGF(9),COLVEL,AXISYM
1 /FNDRI/IDIDX0,IDJDX0,IDKDX0,IDF(8),IJPW(6),IJPS(6),IDSG(6),
1 IDR(6),JDR(6),KDR(6) /FNDRL/NXNE1,NYNE1,NZNE1
1 /GENI/ IG1(2),NXNYST,IG2(6),NFM,IG3(50)
1 /RDATA/TINY,RDT1(16),RINNER,RDT2(67)
1 /IDATA/NX,NY,NZ,IDT1(117) /LDATA/LDT1(19),BFC,LDT2(64)
LOGICAL LGF,COLVEL,AXISYM,SLD,AVEVEL,LDT1,BFC,LDT2,NXNE1,NYNE1,
1 NZNE1
COMMON/F0/LB1(29),KXYDX,LB2(274)
COMMON/GEODMN0/I3DAEX,I3DANY,I3DAHZ,I3DVOL,I3DDXG,I3DDYG,I3DDZG,
1 I3DDX,I3DDY,I3DDZ,I2DAWB,I2DASB,I2DALB
COMMON/CCTDBL/LCUTDBL
LOGICAL LCUTDBL
COMMON/INTDMN/IDMN,INTDM1(8)
LOGICAL LNB2CMM,LNB2CPP,BLK2CM,BLK2CP
LOGICAL LSLDR,LSLDIR,BLKM,BLKP,LNB2CM,LNB2CP,LCUTCL
C
IADS= (IZZ-1)*NFM
C.... Main term dU/dX:
IF(NXNE1) THEN
IF(COLVEL) THEN
IF(BFC) THEN
IADW= (IZZ-1)*NXNYST
IDX=IDIDX0+IADW; JDX=IDJDX0+IADW; KDX=IDKDX0+IADW
DO 10 IJ= 1,NXNYST
F(IDUDX0+IJ)= F(IDUDI0+IJ)*F(IDX+IJ)+F(IDUDJ0+IJ)*F(JDX+IJ)
1 + F(IDUDK0+IJ)*F(KDX+IJ)
10 CONTINUE
ELSE
CALL FNDFDX(IUG0(IPH)+IADS,IDUDX0,IZZ,.FALSE.)
ENDIF
ELSEIF(LCUTDBL) THEN
IDIR=3; IDIP=4
IADW=(IZZ-1)*NXNYST; IJ=0
DO IX=1,NX
DO IY=1,NY
IJ=IJ+1
IF(SLD(IJ)) THEN
F(IDUDX0+IJ)=0.
ELSE
IJK=IJ+IADW
IF(IX>1) THEN
IJKNBM=IJK-NY; LNB2CM=LCUTCL(IJKNBM)
IF(LNB2CM) THEN
IPBCNB=NINT(F(NPBIJK(IDMN)+IJKNBM))
ITNBM=ISHPB(IDMN,IPBCNB)
LNB2CM=F(PBC_CCTYP+ITNBM)==2..AND.
1 F(I3DAEX+IJKNBM)<=0..AND.F(PBC_3RDAR(3)+ITNBM)>0.
ENDIF
IF(IX>2.AND..NOT.LNB2CM) THEN
IJKNBM=IJKNBM-NY; LNB2CMM=LCUTCL(IJKNBM)
IF(LNB2CMM) THEN
IPBCNB=NINT(F(NPBIJK(IDMN)+IJKNBM))
ITNBMM=ISHPB(IDMN,IPBCNB)
LNB2CMM=F(PBC_CCTYP+ITNBMM)==2..AND.
1 F(I3DAEX+IJKNBM-NY)<=0..AND.F(PBC_3RDAR(3)+ITNBMM)>0.
ENDIF
ELSE
LNB2CMM=.FALSE.
ENDIF
BLK2CM=F(I3DAEX+IJKNBM)<=0..AND..NOT.LNB2CM.OR.
1 F(I3DAEX+IJKNBM-NY)<=0..AND..NOT.LNB2CMM
ELSE
LNB2CM=.FALSE.; LNB2CMM=.FALSE.; BLK2CM=.FALSE.
ENDIF
c
IF(IX0.
ENDIF
IF(IX0.
ENDIF
ELSE
LNB2CPP=.FALSE.
ENDIF
BLK2CP=F(I3DAEX+IJK)<=0..AND..NOT.LNB2CP.OR.
1 F(I3DAEX+IJKNBP)<=0..AND..NOT.LNB2CPP
ELSE
LNB2CP=.FALSE.; LNB2CPP=.FALSE.; BLK2CP=.FALSE.
ENDIF
c
BLKM=LSLDR(IJ,IDIP).AND..NOT.LNB2CM.OR.BLK2CM
BLKP=LSLDR(IJ,IDIR).AND..NOT.LNB2CP.OR.BLK2CP
IF(BLKM .AND. BLKP) THEN
F(IDUDX0+IJ)=0.
ELSE
IDXC=KXYDX+IJ
IF(AVEVEL) THEN
IFC=IUG0(IPH)+IJ
DSC=0.5*F(IDXC); FC=F(IFC)
IF(.NOT.BLKM) THEN
IF(LNB2CM) THEN
DSM=F(PBC_3RDDST(3)+ITNBM)
c FM=F(PBC_3RDVEL(3)+ITNBM)
ELSE
DSM=0.5*F(IDXC+IJPS(IDIP))
c FM=F(IFC+IJPS(IDIP))
ENDIF
FM=F(IFC+IJPS(IDIP))
ENDIF
RDSM= 1./(DSM+DSC+TINY)
c
IF(.NOT.BLKP) THEN
IF(LNB2CP) THEN
DSP=F(PBC_3RDDST(4)+ITNBP)
c FP=F(PBC_3RDVEL(4)+ITNBP)
ELSE
DSP=0.5*F(IDXC+IJPS(IDIR))
c FP=F(IFC+IJPS(IDIR))
ENDIF
FP=F(IFC+IJPS(IDIR))
ENDIF
RDSP=1./(DSC+DSP+TINY)
c
IF(BLKM) THEN
FWLM=FC-(FP-FC)*DSC*RDSP
FWLP=(FC*DSP+FP*DSC)*RDSP
F(IDUDX0+IJ)=(FWLP-FWLM)/(F(IDXC)+TINY)
ELSEIF(BLKP) THEN
FWLM=(FM*DSC+FC*DSM)*RDSM
FWLP=FC-(FM-FC)*DSC*RDSM
F(IDUDX0+IJ)=(FWLP-FWLM)/(F(IDXC)+TINY)
ELSE !blkm & !blkp
FWLM=(FM*DSC+FC*DSM)*RDSM
FWLP=(FC*DSP+FP*DSC)*RDSP
F(IDUDX0+IJ)=(FWLP-FWLM)/(F(IDXC)+TINY)
ENDIF
ELSE ! not.avevel
IUC=L0F(3+IPH-1)+IJ
IF(BLKM) THEN
IF(LSLDIR(IJ+IADW+IJPW(IDIR),IDIR).AND.
1 .NOT.LNB2CPP) THEN
F(IDUDX0+IJ)=0.
ELSE
IF(LNB2CPP) THEN
FP=F(PBC_3RDVEL(4)+ITNBPP)
ELSE
FP=F(IUC+IJPS(IDIR))
ENDIF
F(IDUDX0+IJ)=(FP-F(IUC))/(DLL(IJ,IZZ,IDIR)+TINY)
ENDIF
ELSEIF(BLKP) THEN
IF(LSLDIR(IJ+IADW+IJPW(IDIP),IDIP).AND.
1 .NOT.LNB2CMM) THEN
F(IDUDX0+IJ)=0.
ELSE
IUM=IUC+IJPS(IDIP)
IF(LNB2CMM) THEN
FM=F(PBC_3RDVEL(3)+ITNBMM)
ELSE
FM=F(IUM+IJPS(IDIP))
ENDIF
F(IDUDX0+IJ)=(F(IUM)-FM)/DLL(IJ,IZZ,IDIP)
ENDIF
ELSE
IF(LNB2CP) THEN
FP=F(PBC_3RDVEL(4)+ITNBP)
ELSE
FP=F(IUC)
ENDIF
IF(LNB2CM) THEN
FM=F(PBC_3RDVEL(3)+ITNBM)
ELSE
FM=F(IUC+IJPS(IDIP))
ENDIF
F(IDUDX0+IJ)=(FP-FM)/F(IDXC)
ENDIF
ENDIF ! blkm or blkp
ENDIF ! blkm & blkp
ENDIF ! sld
ENDDO
ENDDO
ELSE
C.... Treatment of staggered velocities on non-BFC grids:
IF(AVEVEL) THEN
CALL FNDFDX(IUG0(IPH),IDUDX0,IZZ,.FALSE.)
ELSE
CALL AVDVEL(3,4,L0F(3+IPH-1),IDUDX0,IZZ)
ENDIF
ENDIF
ELSE
CALL ZERNUM(IDUDX0,NXNYST)
ENDIF
C.... Additional term for polar-cylindrical geometry:
IF(AXISYM .AND. IVG0(IPH)/=0) THEN
IVC0= ITWO(IVG0(IPH)+IADS,IVG0(IPH),COLVEL)
IF(BFC) THEN
C.... For BFC this term appears only if NX=1:
DO 20 IY= 1,NY
IF(.NOT.SLD(IY)) THEN
IDUDX= IDUDX0+IY
XCP = COORD1(1,IY,IZZ,1)
YCP = COORD1(1,IY,IZZ,2)
RCP = SQRT(XCP*XCP+YCP*YCP) + RINNER
F(IDUDX)= F(IDUDX) + F(IVC0+IY)/(RCP + TINY)
ENDIF
20 CONTINUE
ELSE
IR0= L0F(LXYR)
DO 30 IJ= 1,NXNYST
IF(.NOT.SLD(IJ)) THEN
IDUDX = IDUDX0+IJ
F(IDUDX)= F(IDUDX) + F(IVC0+IJ)/(F(IR0+IJ) + TINY)
ENDIF
30 CONTINUE
ENDIF
ENDIF
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C
SUBROUTINE GFDUDY(IPH,IZZ,IDUDY0)
INCLUDE 'cmndmn'
INCLUDE 'farray'
INCLUDE 'cmnrot'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'd_earth/pbcstr'
COMMON /CMNGN/IUG0(2),IVG0(2),IWG0(2),IC1(12),IDUDI0,IDUDJ0,
1 IDUDK0,IC2(6) /GENFL/LGF(9),COLVEL,AXISYM
1 /FNDRI/IDF1(3),IDIDY0,IDJDY0,IDKDY0,IDF2(5),IJPW(6),
1 IJPS(6),IDSG(6),IDR(6),JDR(6),KDR(6)
1 /GENI/ IG1(2),NXNYST,IG2(6),NFM,IG3(6),NXNYNZ,IG4(43)
1 /RDATA/TINY,RDT1(16),RINNER,RDT2(43),FIXVAL,RDT3(23)
1 /IDATA/NX,NY,NZ,IDT1(117) /LDATA/LDT1(19),BFC,LDT2(64)
COMMON /LROT/LROTOR
LOGICAL LGF,COLVEL,AXISYM,SLD,LDT1,BFC,LDT2,DOIT,LROTOR
COMMON/F0/IFIL1(30),KXYDY,IFIL2(39),KAP,IFIL3(233)
COMMON/GEODMN0/I3DAEX,I3DANY,I3DAHZ,I3DVOL,I3DDXG,I3DDYG,I3DDZG,
1 I3DDX,I3DDY,I3DDZ,I2DAWB,I2DASB,I2DALB
COMMON/CCTDBL/LCUTDBL
LOGICAL LCUTDBL,LCUTCL
COMMON/INTDMN/IDMN,INTDM1(8)
LOGICAL LSLDR,BLKM,BLKP,LNB2CM,LNB2CP,BLK2CM,BLK2CP
C
C.... Main term dU/dY:
IUC0= ITWO(IUG0(IPH)+(IZZ-1)*NFM,IUG0(IPH),COLVEL)
IF(BFC) THEN
IADW= (IZZ-1)*NXNYST
IDY=IDIDY0+IADW; JDY=IDJDY0+IADW; KDY=IDKDY0+IADW
DO 10 IJ= 1,NXNYST
F(IDUDY0+IJ)= F(IDUDI0+IJ)*F(IDY+IJ)+F(IDUDJ0+IJ)*F(JDY+IJ)
1 + F(IDUDK0+IJ)*F(KDY+IJ)
10 CONTINUE
ELSEIF(LCUTDBL) THEN
IADZ=(IZZ-1)*NXNYST
DO IJ=1,NXNYST
IF(SLD(IJ)) THEN
F(IDUDY0+IJ)=0.
ELSE
IDIR=1; IDIP=2
IJK=IJ+IADZ; IY=1+MOD(IJ-1,NY)
IF(IY>1) THEN
IJKNBM=IJK-1; LNB2CM=LCUTCL(IJKNBM)
IF(LNB2CM) THEN
IPBCNB=NINT(F(NPBIJK(IDMN)+IJKNBM))
ITNBM=ISHPB(IDMN,IPBCNB)
LNB2CM=F(PBC_CCTYP+ITNBM)==2..AND.F(I3DANY+IJK-1)<=0.
1 .AND.F(PBC_3RDAR(1)+ITNBM)>0..AND.F(PBC_3RDAR(3)+ITNBM)>0.
ENDIF
BLK2CM=F(I3DANY+IJK-1)<=0..AND..NOT.LNB2CM
ELSE
LNB2CM=.FALSE.; BLK2CM=.FALSE.
ENDIF
c
IF(IY0..AND.F(PBC_3RDAR(3)+ITNBP)>0.
ENDIF
BLK2CP=F(I3DANY+IJK)<=0..AND..NOT.LNB2CP
ELSE
LNB2CP=.FALSE.; BLK2CP=.FALSE.
ENDIF
c
BLKM=LSLDR(IJ,IDIP).AND..NOT.LNB2CM.OR.BLK2CM
BLKP=LSLDR(IJ,IDIR).AND..NOT.LNB2CP.OR.BLK2CP
c
IF(BLKM.AND.BLKP) THEN
F(IDUDY0+IJ)=0.
ELSE
IDYC=KXYDY+IJ; IFC=IUC0+IJ
DSC=0.5*F(IDYC); FC=F(IFC)
c
IF(.NOT.BLKM) THEN
IF(LNB2CM) THEN
DSM=F(PBC_3RDDST(1)+ITNBM)
IF(NINT(F(ITNBM+PBC_IX))0) THEN
IY= 1+MOD(IJ-1,NY)
F(IDUDY)= F(IDUDY)-(F(IUC0+IJ)+ANGV(IROTC)*F(IR0+IJ))/
1 (F(IR0+IJ) + TINY)
ELSE
F(IDUDY)= F(IDUDY)-F(IUC0+IJ)/(F(IR0+IJ) + TINY)
ENDIF
ELSE
F(IDUDY)= F(IDUDY) - F(IUC0+IJ)/(F(IR0+IJ) + TINY)
ENDIF
ENDIF
30 CONTINUE
ENDIF
ENDIF
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C
c
SUBROUTINE GFDUDZ(IPH,IZZ,IDUDZ0)
INCLUDE 'cmndmn'
INCLUDE 'farray'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'd_earth/pbcstr'
COMMON /CMNGN/IUG0(2),IVG0(2),IWG0(2),IUIZM0(2),IUIZP0(2),IC1(8),
1 IDUDI0,IDUDJ0,IDUDK0,IC2(6)
1 /FNDRI/IDF1(6),IDIDZ0,IDJDZ0,IDKDZ0,IDF2(38)
1 /GENI/ IG1(2),NXNYST,IG2(6),NFM,IG3(50)
1 /LDATA/LDT1(19),BFC,LDT2(64) /GENFL/LGF(9),COLVEL,AXISYM
LOGICAL LGF,COLVEL,AXISYM,LDT1,BFC,LDT2
COMMON/RDATA/TINY,RD1(84)
COMMON/IDATA/NX,NY,NZ,IDT1(117)
COMMON/F0/LB1(16),KZDZ,LB2(53),KAP,LB3(233)
COMMON/GEODMN0/I3DAEX,I3DANY,I3DAHZ,I3DVOL,I3DDXG,I3DDYG,I3DDZG,
1 I3DDX,I3DDY,I3DDZ,I2DAWB,I2DASB,I2DALB
COMMON/CCTDBL/LCUTDBL
LOGICAL LCUTDBL,LCUTCL
COMMON/INTDMN/IDMN,INTDM1(8)
LOGICAL SLD,LSLDR,BLKM,BLKP,LNB2CM,LNB2CP,BLK2CM,BLK2CP
C
IF(COLVEL) THEN
IF(BFC) THEN
IADW= (IZZ-1)*NXNYST
IDZ=IDIDZ0+IADW; JDZ=IDJDZ0+IADW; KDZ=IDKDZ0+IADW
DO 10 IJ= 1,NXNYST
F(IDUDZ0+IJ)= F(IDUDI0+IJ)*F(IDZ+IJ)+F(IDUDJ0+IJ)*F(JDZ+IJ)
1 + F(IDUDK0+IJ)*F(KDZ+IJ)
10 CONTINUE
ELSE
CALL FNDFDZ(IUG0(IPH)+(IZZ-1)*NFM,IDUDZ0,IZZ,.FALSE.)
ENDIF
ELSEIF(LCUTDBL) THEN
IFIZM0=IUIZM0(IPH); IFIZ0=IUG0(IPH); IFIZP0=IUIZP0(IPH)
IDZC=KZDZ+IZZ; DKDZ=1./(F(IDZC)+TINY)
DSC=.5*F(IDZC); DSM=.5*F(IDZC-1); DSP=.5*F(IDZC+1)
RDSM=1./(DSM+DSC+TINY); RDSP=1./(DSC+DSP+TINY)
IADZ=(IZZ-1)*NXNYST
DO IJ=1,NXNYST
IF(SLD(IJ)) THEN
F(IDUDZ0+IJ)=0.0
ELSE
IJK=IJ+IADZ
IDIR=5; IDIP=6
c
IF(IZZ>1) THEN
IJKNBM=IJK-NXNYST; LNB2CM=LCUTCL(IJKNBM)
IF(LNB2CM) THEN
IPBCNB=NINT(F(NPBIJK(IDMN)+IJKNBM))
ITNBM=ISHPB(IDMN,IPBCNB)
LNB2CM=F(PBC_CCTYP+ITNBM)==2..AND.
1 F(I3DAHZ+IJK-NXNYST)<=0..AND.
1 F(PBC_3RDAR(5)+ITNBM)>0..AND.F(PBC_3RDAR(3)+ITNBM)>0.
ENDIF
BLK2CM=F(I3DAHZ+IJK-NXNYST)<=0..AND..NOT.LNB2CM
ELSE
LNB2CM=.FALSE.; BLK2CM=.FALSE.
ENDIF
c
IF(IZZ0..AND.F(PBC_3RDAR(3)+ITNBP)>0.
ENDIF
BLK2CP=F(I3DAHZ+IJK)<=0..AND..NOT.LNB2CP
ELSE
LNB2CP=.FALSE.; BLK2CP=.FALSE.
ENDIF
c
BLKM=LSLDR(IJ,IDIP).AND..NOT.LNB2CM.OR.BLK2CM
BLKP=LSLDR(IJ,IDIR).AND..NOT.LNB2CP.OR.BLK2CP
c
IF(BLKM.AND.BLKP) THEN
F(IDUDZ0+IJ)=0.0
ELSE
FC=F(IFIZ0+IJ)
c
IF(.NOT.BLKM) THEN
IF(LNB2CM) THEN
DSM=F(PBC_3RDDST(5)+ITNBM)
IF(NINT(F(ITNBM+PBC_IZ))
SUBROUTINE GFDVDX(IPH,IZZ,IDVDX0)
INCLUDE 'cmndmn'
INCLUDE 'farray'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'd_earth/pbcstr'
COMMON /CMNGN/IUG0(2),IVG0(2),IWG0(2),IC1(15),IDVDI0,IDVDJ0,
1 IDVDK0,IC2(3)
COMMON/GENI/NXNY,NXM1NY,NXNYST,NDIR,KDUMM,NXST,NYST,NXNYFM,IZ,
1 NFM,IG2(50)
COMMON/LDATA/LDT1(7),XCYCLE,LDT2(11),BFC,LDT3(64)
COMMON/GENFL/LGF(9),COLVEL,AXISYM
COMMON/FNDRI/IDIDX0,IDJDX0,IDKDX0,IDF(8),IJPW(6),IJPS(6),IDSG(6),
1 IDR(6),JDR(6),KDR(6)
LOGICAL LGF,COLVEL,AXISYM,LDT1,XCYCLE,LDT2,BFC,LDT3
COMMON/IDATA/NX,NY,NZ,IDT1(117)
LOGICAL XCYCZ,NEZ,SLD
COMMON/F0/IFIL1(29),KXYDX,IFIL2(274)
COMMON/RDATA/TINY,RD1(84)
COMMON/GEODMN0/I3DAEX,I3DANY,I3DAHZ,I3DVOL,I3DDXG,I3DDYG,I3DDZG,
1 I3DDX,I3DDY,I3DDZ,I2DAWB,I2DASB,I2DALB
COMMON/CCTDBL/LCUTDBL
LOGICAL LCUTDBL, LCUTCL
COMMON/INTDMN/IDMN,INTDM1(8)
LOGICAL LSLDR,BLKM,BLKP,LNB2CM,LNB2CP,BLK2CM,BLK2CP
C
IF(BFC) THEN
IADW= (IZZ-1)*NXNYST
IDX=IDIDX0+IADW; JDX=IDJDX0+IADW; KDX=IDKDX0+IADW
DO 10 IJ= 1,NXNYST
F(IDVDX0+IJ)= F(IDVDI0+IJ)*F(IDX+IJ) + F(IDVDJ0+IJ)*F(JDX+IJ)
1 + F(IDVDK0+IJ)*F(KDX+IJ)
10 CONTINUE
ELSEIF(LCUTDBL) THEN
IF0=IVG0(IPH)
IADZ=(IZZ-1)*NXNYST
DO IJ=1,NXNYST
IF(SLD(IJ)) THEN
F(IDVDX0+IJ)=0.
ELSE
IDIR=3; IDIP=4
IJK=IJ+IADZ; IX=1+(IJ-1)/NY
IF(IX>1) THEN
IJKNBM=IJK-NY; LNB2CM=LCUTCL(IJKNBM)
IF(LNB2CM) THEN
IPBCNB=NINT(F(NPBIJK(IDMN)+IJKNBM))
ITNBM=ISHPB(IDMN,IPBCNB)
LNB2CM=F(PBC_CCTYP+ITNBM)==2..AND.F(I3DAEX+IJK-NY)<=0.
1 .AND.F(PBC_3RDAR(3)+ITNBM)>0..AND.F(PBC_3RDAR(1)+ITNBM)>0.
ENDIF
BLK2CM=F(I3DAEX+IJK-NY)<=0..AND..NOT.LNB2CM
ELSE
LNB2CM=.FALSE.; BLK2CM=.FALSE.
ENDIF
c
IF(IX0..AND.F(PBC_3RDAR(1)+ITNBP)>0.
ENDIF
BLK2CP=F(I3DAEX+IJK)<=0..AND..NOT.LNB2CP
ELSE
LNB2CP=.FALSE.; BLK2CP=.FALSE.
ENDIF
c
BLKM=LSLDR(IJ,IDIP).AND..NOT.LNB2CM.OR.BLK2CM
BLKP=LSLDR(IJ,IDIR).AND..NOT.LNB2CP.OR.BLK2CP
IF(BLKM.AND.BLKP) THEN
F(IDVDX0+IJ)=0.
ELSE
IDXC=KXYDX+IJ; IFC=IF0+IJ
DSC=0.5*F(IDXC); FC=F(IFC)
c
IF(.NOT.BLKM) THEN
IF(LNB2CM) THEN
DSM=F(PBC_3RDDST(3)+ITNBM)
IF(NINT(F(ITNBM+PBC_IY))
SUBROUTINE GFDVDY(IPH,IZZ,IDVDY0,AVEVEL)
INCLUDE 'cmndmn'
INCLUDE 'farray'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'd_earth/pbcstr'
COMMON /CMNGN/IUG0(2),IVG0(2),IWG0(2),IC1(15),IDVDI0,IDVDJ0,
1 IDVDK0,IC2(3) /GENI/IG1(2),NXNYST,IG2(6),NFM,IG3(50)
1 /LDATA/LDT1(19),BFC,LDT2(64) /GENFL/LGF(9),COLVEL,AXISYM
1 /FNDRI/IDF1(3),IDIDY0,IDJDY0,IDKDY0,IDF2(5),IJPW(6),IJPS(6)
1 ,IDSG(6),IDR(6),JDR(6),KDR(6)
LOGICAL LGF,COLVEL,AXISYM,AVEVEL,LDT1,BFC,LDT2
COMMON/IDATA/NX,NY,NZ,IDT1(117)
COMMON/RDATA/TINY,RD1(84)
COMMON/F0/LB1(30),KXYDY,LB2(273)
COMMON/GEODMN0/I3DAEX,I3DANY,I3DAHZ,I3DVOL,I3DDXG,I3DDYG,I3DDZG,
1 I3DDX,I3DDY,I3DDZ,I2DAWB,I2DASB,I2DALB
COMMON/CCTDBL/LCUTDBL
LOGICAL LCUTDBL
COMMON/INTDMN/IDMN,INTDM1(8)
LOGICAL LNB2CMM,LNB2CPP,BLK2CM,BLK2CP
LOGICAL SLD,LSLDR,LSLDIR,BLKM,BLKP,LNB2CM,LNB2CP,LCUTCL
C
IF(COLVEL) THEN
IF(BFC) THEN
IADW= (IZZ-1)*NXNYST
IDY=IDIDY0+IADW; JDY=IDJDY0+IADW; KDY=IDKDY0+IADW
DO 10 IJ= 1,NXNYST
F(IDVDY0+IJ)= F(IDVDI0+IJ)*F(IDY+IJ)+F(IDVDJ0+IJ)*F(JDY+IJ)
1 + F(IDVDK0+IJ)*F(KDY+IJ)
10 CONTINUE
ELSE
CALL FNDFDY(IVG0(IPH)+(IZZ-1)*NFM,IDVDY0,IZZ,.FALSE.)
ENDIF
ELSEIF(LCUTDBL) THEN
IDIR=1; IDIP=2
IADW=(IZZ-1)*NXNYST; IJ=0
DO IX=1,NX
DO IY=1,NY
IJ=IJ+1
IF(SLD(IJ)) THEN
F(IDVDY0+IJ)=0.
ELSE
IJK=IJ+IADW
IF(IY>1) THEN
IJKNBM=IJK-1; LNB2CM=LCUTCL(IJKNBM)
IF(LNB2CM) THEN
IPBCNB=NINT(F(NPBIJK(IDMN)+IJKNBM))
ITNBM=ISHPB(IDMN,IPBCNB)
LNB2CM=F(PBC_CCTYP+ITNBM)==2..AND.
1 F(I3DANY+IJKNBM)<=0..AND.F(PBC_3RDAR(1)+ITNBM)>0.
ENDIF
IF(IY>2.AND..NOT.LNB2CM) THEN
IJKNBM=IJKNBM-1; LNB2CMM=LCUTCL(IJKNBM)
IF(LNB2CMM) THEN
IPBCNB=NINT(F(NPBIJK(IDMN)+IJKNBM))
ITNBMM=ISHPB(IDMN,IPBCNB)
LNB2CMM=F(PBC_CCTYP+ITNBMM)==2..AND.
1 F(I3DANY+IJK-2)<=0..AND.F(PBC_3RDAR(1)+ITNBMM)>0.
ENDIF
ELSE
LNB2CMM=.FALSE.
ENDIF
BLK2CM=F(I3DANY+IJKNBM)<=0..AND..NOT.LNB2CM.OR.
1 F(I3DANY+IJKNBM-1)<=0..AND..NOT.LNB2CMM
ELSE
LNB2CM=.FALSE.; LNB2CMM=.FALSE.; BLK2CM=.FALSE.
ENDIF
c
IF(IY0.
ENDIF
IF(IY0.
ENDIF
ELSE
LNB2CPP=.FALSE.
ENDIF
BLK2CP=F(I3DANY+IJK)<=0..AND..NOT.LNB2CP.OR.
1 F(I3DANY+IJKNBP)<=0..AND..NOT.LNB2CPP
ELSE
LNB2CP=.FALSE.; LNB2CPP=.FALSE.; BLK2CP=.FALSE.
ENDIF
c
BLKM=LSLDR(IJ,IDIP).AND..NOT.LNB2CM.OR.BLK2CM
BLKP=LSLDR(IJ,IDIR).AND..NOT.LNB2CP.OR.BLK2CP
IF(BLKM.AND.BLKP) THEN
F(IDVDY0+IJ)=0.
ELSE
IDYC=KXYDY+IJ
IF(AVEVEL) THEN
IFC=IVG0(IPH)+IJ
DSC=0.5*F(IDYC); FC=F(IFC)
IF(.NOT.BLKM) THEN
IF(LNB2CM) THEN
DSM=F(PBC_3RDDST(1)+ITNBM)
c FM=F(PBC_3RDVEL(1)+ITNBM)
ELSE
DSM=0.5*F(IDYC+IJPS(IDIP))
c FM=F(IFC+IJPS(IDIP))
ENDIF
FM=F(IFC+IJPS(IDIP))
ENDIF
RDSM=1./(DSM+DSC+TINY)
c
IF(.NOT.BLKP) THEN
IF(LNB2CP) THEN
DSP=F(PBC_3RDDST(2)+ITNBP)
c FP=F(PBC_3RDVEL(2)+ITNBP)
ELSE
DSP=0.5*F(IDYC+IJPS(IDIR))
c FP=F(IFC+IJPS(IDIR))
ENDIF
FP=F(IFC+IJPS(IDIR))
ENDIF
RDSP=1./(DSC+DSP+TINY)
c
IF(BLKM) THEN
FWLM=FC-(FP-FC)*DSC*RDSP
FWLP=(FC*DSP+FP*DSC)*RDSP
F(IDVDY0+IJ)=(FWLP-FWLM)/(F(IDYC)+TINY)
ELSEIF(BLKP) THEN
FWLM=(FM*DSC+FC*DSM)*RDSM
FWLP=FC-(FM-FC)*DSC*RDSM
F(IDVDY0+IJ)=(FWLP-FWLM)/(F(IDYC)+TINY)
ELSE !blkm & !blkp
FWLM=(FM*DSC+FC*DSM)*RDSM
FWLP=(FC*DSP+FP*DSC)*RDSP
F(IDVDY0+IJ)=(FWLP-FWLM)/(F(IDYC)+TINY)
ENDIF
ELSE ! not.avevel
IVC=L0F(5+IPH-1)+IJ
IF(BLKM) THEN
IF(LSLDIR(IJ+IADW+IJPW(IDIR),IDIR).AND.
1 .NOT.LNB2CPP) THEN
F(IDVDY0+IJ)=0.
ELSE
IF(LNB2CPP) THEN
FP=F(PBC_3RDVEL(2)+ITNBPP)
ELSE
FP=F(IVC+IJPS(IDIR))
ENDIF
F(IDVDY0+IJ)=(FP-F(IVC))/(DLL(IJ,IZZ,IDIR)+TINY)
ENDIF
ELSEIF(BLKP) THEN
IF(LSLDIR(IJ+IADW+IJPW(IDIP),IDIP).AND.
1 .NOT.LNB2CMM) THEN
F(IDVDY0+IJ)=0.
ELSE
IVM=IVC+IJPS(IDIP)
IF(LNB2CMM) THEN
FM=F(PBC_3RDVEL(1)+ITNBMM)
ELSE
FM=F(IVM+IJPS(IDIP))
ENDIF
F(IDVDY0+IJ)=(F(IVM)-FM)/DLL(IJ,IZZ,IDIP)
ENDIF
ELSE
IF(LNB2CP) THEN
FP=F(PBC_3RDVEL(2)+ITNBP)
ELSE
FP=F(IVC)
ENDIF
IF(LNB2CM) THEN
FM=F(PBC_3RDVEL(1)+ITNBM)
ELSE
FM=F(IVC+IJPS(IDIP))
ENDIF
F(IDVDY0+IJ)=(FP-FM)/F(IDYC)
ENDIF
ENDIF ! blkm or blkp
ENDIF ! blkm & blkp
ENDIF ! sld
ENDDO
ENDDO
ELSE
C.... Treatment of staggered velocities:
IF(AVEVEL) THEN
CALL FNDFDY(IVG0(IPH),IDVDY0,IZZ,.FALSE.)
ELSE
CALL AVDVEL(1,2,L0F(5+IPH-1),IDVDY0,IZZ)
ENDIF
ENDIF
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C
c
SUBROUTINE GFDVDZ(IPH,IZZ,IDVDZ0)
INCLUDE 'cmndmn'
INCLUDE 'farray'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'd_earth/pbcstr'
COMMON /CMNGN/IUG0(2),IVG0(2),IWG0(2),IC1(4),IVIZM0(2),IVIZP0(2),
1 IC2(7),IDVDI0,IDVDJ0,IDVDK0,IC3(3)
1 /FNDRI/IDF1(6),IDIDZ0,IDJDZ0,IDKDZ0,IDF2(38)
1 /GENI/ IG1(2),NXNYST,IG2(6),NFM,IG3(50)
1 /LDATA/LDT1(19),BFC,LDT2(64) /GENFL/LGF(9),COLVEL,AXISYM
LOGICAL LGF,COLVEL,AXISYM,LDT1,BFC,LDT2
COMMON/RDATA/TINY,RD1(84)
COMMON/IDATA/NX,NY,NZ,IDT1(117)
COMMON/F0/LB1(16),KZDZ,LB2(287)
COMMON/GEODMN0/I3DAEX,I3DANY,I3DAHZ,I3DVOL,I3DDXG,I3DDYG,I3DDZG,
1 I3DDX,I3DDY,I3DDZ,I2DAWB,I2DASB,I2DALB
COMMON/CCTDBL/LCUTDBL
LOGICAL LCUTDBL, LCUTCL
COMMON/INTDMN/IDMN,INTDM1(8)
LOGICAL SLD,LSLDR,BLKM,BLKP,LNB2CT,LNB2CM,LNB2CP,BLK2CM,BLK2CP
C
IF(COLVEL) THEN
IF(BFC) THEN
IADW= (IZZ-1)*NXNYST
IDZ=IDIDZ0+IADW; JDZ=IDJDZ0+IADW; KDZ=IDKDZ0+IADW
DO 10 IJ= 1,NXNYST
F(IDVDZ0+IJ)= F(IDVDI0+IJ)*F(IDZ+IJ)+F(IDVDJ0+IJ)*F(JDZ+IJ)
1 + F(IDVDK0+IJ)*F(KDZ+IJ)
10 CONTINUE
ELSE
CALL FNDFDZ(IVG0(IPH)+(IZZ-1)*NFM,IDVDZ0,IZZ,.FALSE.)
ENDIF
ELSEIF(LCUTDBL) THEN
IFIZM0=IVIZM0(IPH); IFIZ0=IVG0(IPH); IFIZP0=IVIZP0(IPH)
IDZC=KZDZ+IZZ; DKDZ=1./(F(IDZC)+TINY)
DSC=.5*F(IDZC); DSM=.5*F(IDZC-1); DSP=.5*F(IDZC+1)
RDSM=1./(DSM+DSC+TINY); RDSP=1./(DSC+DSP+TINY)
IADZ=(IZZ-1)*NXNYST
DO IJ=1,NXNYST
IF(SLD(IJ)) THEN
F(IDVDZ0+IJ)=0.
ELSE
IJK=IJ+IADZ
IDIR=5; IDIP=6
c
IF(IZZ>1) THEN
IJKNBM=IJK-NXNYST; LNB2CM=LCUTCL(IJKNBM)
IF(LNB2CM) THEN
IPBCNB=NINT(F(NPBIJK(IDMN)+IJKNBM))
ITNBM=ISHPB(IDMN,IPBCNB)
LNB2CT=F(PBC_CCTYP+ITNBM)==2..AND.
1 F(I3DAHZ+IJK-NXNYST)<=0..AND.
1 F(PBC_3RDAR(5)+ITNBM)>0..AND.F(PBC_3RDAR(1)+ITNBM)>0.
ENDIF
BLK2CM=F(I3DAHZ+IJK-NXNYST)<=0..AND..NOT.LNB2CM
ELSE
LNB2CM=.FALSE.; BLK2CM=.FALSE.
ENDIF
c
IF(IZZ0..AND.F(PBC_3RDAR(1)+ITNBP)>0.
ENDIF
BLK2CP=F(I3DAHZ+IJK)<=0..AND..NOT.LNB2CP
ELSE
LNB2CP=.FALSE.; BLK2CP=.FALSE.
ENDIF
c
BLKM=LSLDR(IJ,IDIP).AND..NOT.LNB2CM.OR.BLK2CM
BLKP=LSLDR(IJ,IDIR).AND..NOT.LNB2CP.OR.BLK2CP
c
IF(BLKM.AND.BLKP) THEN
F(IDVDZ0+IJ)=0.
ELSE
FC=F(IFIZ0+IJ)
IF(.NOT.BLKM) THEN
IF(LNB2CM) THEN
DSM=F(PBC_3RDDST(5)+ITNBM)
IF(NINT(F(ITNBM+PBC_IY))
SUBROUTINE GFDWDX(IPH,IZZ,IDWDX0)
INCLUDE 'cmndmn'
INCLUDE 'farray'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'd_earth/pbcstr'
COMMON /CMNGN/IUG0(2),IVG0(2),IWG0(2),IC(18),IDWDI0,IDWDJ0,IDWDK0
1 /GENI/ IG1(2),NXNYST,IG2(6),NFM,IG3(50)
1 /LDATA/LDT1(19),BFC,LDT2(64) /GENFL/LGF(9),COLVEL,AXISYM
c 1 /FNDRI/IDIDX0,IDJDX0,IDKDX0,IDF(44)
COMMON /FNDRI/IDIDX0,IDJDX0,IDKDX0,IDF(8),IJPW(6),IJPS(6),IDSG(6),
1 IDR(6),JDR(6),KDR(6)
LOGICAL LGF,COLVEL,AXISYM,LDT1,BFC,LDT2
COMMON/RDATA/TINY,RD1(84)
COMMON/IDATA/NX,NY,NZ,IDT1(117)
COMMON/F0/IFIL1(29),KXYDX,IFIL4(274)
COMMON/GEODMN0/I3DAEX,I3DANY,I3DAHZ,I3DVOL,I3DDXG,I3DDYG,I3DDZG,
1 I3DDX,I3DDY,I3DDZ,I2DAWB,I2DASB,I2DALB
COMMON/CCTDBL/LCUTDBL
LOGICAL LCUTDBL, LCUTCL
COMMON/INTDMN/IDMN,INTDM1(8)
LOGICAL SLD,LSLDR,BLKM,BLKP,LNB2CM,LNB2CP,BLK2CM,BLK2CP
C
IF(BFC) THEN
IADW= (IZZ-1)*NXNYST
IDX=IDIDX0+IADW; JDX=IDJDX0+IADW; KDX=IDKDX0+IADW
DO 10 IJ= 1,NXNYST
F(IDWDX0+IJ)= F(IDWDI0+IJ)*F(IDX+IJ) + F(IDWDJ0+IJ)*F(JDX+IJ)
1 + F(IDWDK0+IJ)*F(KDX+IJ)
10 CONTINUE
ELSEIF(LCUTDBL) THEN
IF0=IWG0(IPH)
IADZ=(IZZ-1)*NXNYST
DO IJ=1,NXNYST
IF(SLD(IJ)) THEN
F(IDWDX0+IJ)=0.
ELSE
IDIR=3; IDIP=4
IJK=IJ+IADZ; IX=1+(IJ-1)/NY
IF(IX>1) THEN
IJKNBM=IJK-NY; LNB2CM=LCUTCL(IJKNBM)
IF(LNB2CM) THEN
IPBCNB=NINT(F(NPBIJK(IDMN)+IJKNBM))
ITNBM=ISHPB(IDMN,IPBCNB)
LNB2CM=F(PBC_CCTYP+ITNBM)==2..AND.F(I3DAEX+IJK-NY)<=0.
1 .AND.F(PBC_3RDAR(5)+ITNBM)>0..AND.F(PBC_3RDAR(1)+ITNBM)>0.
ENDIF
BLK2CM=F(I3DAEX+IJK-NY)<=0..AND..NOT.LNB2CM
ELSE
LNB2CM=.FALSE.; BLK2CM=.FALSE.
ENDIF
c
IF(IX0..AND.F(PBC_3RDAR(1)+ITNBP)>0.
ENDIF
BLK2CP=F(I3DAEX+IJK)<=0..AND..NOT.LNB2CP
ELSE
LNB2CP=.FALSE.; BLK2CP=.FALSE.
ENDIF
c
BLKM=LSLDR(IJ,IDIP).AND..NOT.LNB2CM.OR.BLK2CM
BLKP=LSLDR(IJ,IDIR).AND..NOT.LNB2CP.OR.BLK2CP
IF(BLKM.AND.BLKP) THEN
F(IDWDX0+IJ)=0.
ELSE
IDXC=KXYDX+IJ; IFC=IF0+IJ
DSC=0.5*F(IDXC); FC=F(IFC)
c
IF(.NOT.BLKM) THEN
IF(LNB2CM) THEN
DSM=F(PBC_3RDDST(3)+ITNBM)
IF(NINT(F(ITNBM+PBC_IZ))
SUBROUTINE GFDWDY(IPH,IZZ,IDWDY0)
INCLUDE 'cmndmn'
INCLUDE 'farray'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'd_earth/pbcstr'
COMMON /CMNGN/IUG0(2),IVG0(2),IWG0(2),IC(18),IDWDI0,IDWDJ0,IDWDK0
1 /GENI/ IG1(2),NXNYST,IG2(6),NFM,IG3(50)
1 /LDATA/LDT1(19),BFC,LDT2(64) /GENFL/LGF(9),COLVEL,AXISYM
c 1 /FNDRI/IDF1(3),IDIDY0,IDJDY0,IDKDY0,IDF2(41)
COMMON /FNDRI/IDF1(3),IDIDY0,IDJDY0,IDKDY0,IDF2(5),
1 IJPW(6),IJPS(6),IDSG(6),IDR(6),JDR(6),KDR(6)
LOGICAL LGF,COLVEL,AXISYM,LDT1,BFC,LDT2
COMMON/RDATA/TINY,RD1(84)
COMMON/IDATA/NX,NY,NZ,IDT1(117)
COMMON/F0/LB1(30),KXYDY,LB2(273)
COMMON/GEODMN0/I3DAEX,I3DANY,I3DAHZ,I3DVOL,I3DDXG,I3DDYG,I3DDZG,
1 I3DDX,I3DDY,I3DDZ,I2DAWB,I2DASB,I2DALB
COMMON/CCTDBL/LCUTDBL
LOGICAL LCUTDBL, LCUTCL
COMMON/INTDMN/IDMN,INTDM1(8)
LOGICAL SLD,LSLDR,BLKM,BLKP,LNB2CM,LNB2CP,BLK2CM,BLK2CP
C
IF(BFC) THEN
IADW= (IZZ-1)*NXNYST
IDY=IDIDY0+IADW; JDY=IDJDY0+IADW; KDY=IDKDY0+IADW
DO 10 IJ= 1,NXNYST
F(IDWDY0+IJ)= F(IDWDI0+IJ)*F(IDY+IJ)+F(IDWDJ0+IJ)*F(JDY+IJ)
1 + F(IDWDK0+IJ)*F(KDY+IJ)
10 CONTINUE
ELSEIF(LCUTDBL) THEN
IF0=IWG0(IPH)
IADZ=(IZZ-1)*NXNYST
DO IJ=1,NXNYST
IF(SLD(IJ)) THEN
F(IDWDY0+IJ)=0.
ELSE
IDIR=1; IDIP=2
IJK=IJ+IADZ; IY=1+MOD(IJ-1,NY)
IF(IY>1) THEN
IJKNBM=IJK-1; LNB2CM=LCUTCL(IJKNBM)
IF(LNB2CM) THEN
IPBCNB=NINT(F(NPBIJK(IDMN)+IJKNBM))
ITNBM=ISHPB(IDMN,IPBCNB)
LNB2CM=F(PBC_CCTYP+ITNBM)==2..AND.F(I3DANY+IJK-1)<=0.
1 .AND.F(PBC_3RDAR(1)+ITNBM)>0..AND.F(PBC_3RDAR(5)+ITNBM)>0.
ENDIF
BLK2CM=F(I3DANY+IJK-1)<=0..AND..NOT.LNB2CM
ELSE
LNB2CM=.FALSE.; BLK2CM=.FALSE.
ENDIF
c
IF(IY0..AND.F(PBC_3RDAR(5)+ITNBP)>0.
ENDIF
BLK2CP=F(I3DANY+IJK)<=0..AND..NOT.LNB2CP
ELSE
LNB2CP=.FALSE.; BLK2CP=.FALSE.
ENDIF
c
BLKM=LSLDR(IJ,IDIP).AND..NOT.LNB2CM.OR.BLK2CM
BLKP=LSLDR(IJ,IDIR).AND..NOT.LNB2CP.OR.BLK2CP
c
IF(BLKM.AND.BLKP) THEN
F(IDWDY0+IJ)=0.
ELSE
IDYC=KXYDY+IJ; IFC=IF0+IJ
DSC=0.5*F(IDYC); FC=F(IFC)
c
IF(.NOT.BLKM) THEN
IF(LNB2CM) THEN
DSM=F(PBC_3RDDST(1)+ITNBM)
IF(NINT(F(ITNBM+PBC_IZ))
SUBROUTINE GFDWDZ(IPH,IZZ,IDWDZ0,AVEVEL)
INCLUDE 'cmndmn'
INCLUDE 'farray'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'd_earth/pbcstr'
COMMON /CMNGN/IUG0(2),IVG0(2),IWG0(2),IC1(8),IWIZM0(2),IWIZP0(2),
1 IC2(6),IDWDI0,IDWDJ0,IDWDK0
1 /FNDRI/IDF1(6),IDIDZ0,IDJDZ0,IDKDZ0,IDF2(2),IJPW(6),IJPS(6)
1 ,IDSG(6),IDR(6),JDR(6),KDR(6)
1 /GENI/ IG1(2),NXNYST,IG2(6),NFM,IG3(50)
1 /LDATA/LDT1(19),BFC,LDT2(64) /GENFL/LGF(9),COLVEL,AXISYM
LOGICAL AVEVEL,LGF,COLVEL,AXISYM,LDT1,BFC,LDT2
COMMON /IDATA/NX,NY,NZ,IDT1(117)
COMMON/RDATA/TINY,RD1(84)
c COMMON/F0/LB1(16),KZDZ,LB2(287)
COMMON/F0/LB1(16),KZDZ,LB2(148),KXYDZ,LB3(138)
COMMON/GEODMN0/I3DAEX,I3DANY,I3DAHZ,I3DVOL,I3DDXG,I3DDYG,I3DDZG,
1 I3DDX,I3DDY,I3DDZ,I2DAWB,I2DASB,I2DALB
COMMON/CCTDBL/LCUTDBL
LOGICAL LCUTDBL
COMMON/INTDMN/IDMN,INTDM1(8)
LOGICAL LNB2CMM,LNB2CPP,BLK2CM,BLK2CP
LOGICAL SLD,LSLDR,LSLDIR,BLKM,BLKP,LNB2CM,LNB2CP,LCUTCL
C
IF(COLVEL) THEN
IF(BFC) THEN
IADW= (IZZ-1)*NXNYST
IDZ=IDIDZ0+IADW; JDZ=IDJDZ0+IADW; KDZ=IDKDZ0+IADW
DO 10 IJ= 1,NXNYST
F(IDWDZ0+IJ)= F(IDWDI0+IJ)*F(IDZ+IJ)+F(IDWDJ0+IJ)*F(JDZ+IJ)
1 + F(IDWDK0+IJ)*F(KDZ+IJ)
10 CONTINUE
ELSE
CALL FNDFDZ(IWG0(IPH)+(IZZ-1)*NFM,IDWDZ0,IZZ,.FALSE.)
ENDIF
ELSEIF(LCUTDBL) THEN
IDIR=5; IDIP=6
IADW=(IZZ-1)*NXNYST; IJ=0
DO IX=1,NX
DO IY=1,NY
IJ=IJ+1
IF(SLD(IJ)) THEN
F(IDWDZ0+IJ)= 0.0
ELSE
IJK=IJ+IADW
IF(IZZ>1) THEN
IJKNBM=IJK-NXNYST; LNB2CM=LCUTCL(IJKNBM)
IF(LNB2CM) THEN
IPBCNB=NINT(F(NPBIJK(IDMN)+IJKNBM))
ITNBM=ISHPB(IDMN,IPBCNB)
LNB2CM=F(PBC_CCTYP+ITNBM)==2..AND.
1 F(I3DAHZ+IJKNBM)<=0..AND.F(PBC_3RDAR(5)+ITNBM)>0.
ENDIF
IF(IZZ>2.AND..NOT.LNB2CM) THEN
IJKNBM=IJKNBM-NXNYST; LNB2CMM=LCUTCL(IJKNBM)
IF(LNB2CMM) THEN
IPBCNB=NINT(F(NPBIJK(IDMN)+IJKNBM))
ITNBMM=ISHPB(IDMN,IPBCNB)
LNB2CMM=F(PBC_CCTYP+ITNBMM)==2..AND.
1 F(I3DAHZ+IJKNBM-NXNYST)<=0..AND.F(PBC_3RDAR(5)+ITNBMM)>0.
ENDIF
ELSE
LNB2CMM=.FALSE.
ENDIF
BLK2CM=F(I3DAHZ+IJKNBM)<=0..AND..NOT.LNB2CM.OR.
1 F(I3DAHZ+IJKNBM-NXNYST)<=0..AND..NOT.LNB2CMM
ELSE
LNB2CM=.FALSE.; LNB2CMM=.FALSE.; BLK2CM=.FALSE.
ENDIF
c
IF(IZZ0.
ENDIF
IF(IZZ0.
ENDIF
ELSE
LNB2CPP=.FALSE.
ENDIF
BLK2CP=F(I3DAHZ+IJK)<=0..AND..NOT.LNB2CP.OR.
1 F(I3DAHZ+IJKNBP)<=0..AND..NOT.LNB2CPP
ELSE
LNB2CP=.FALSE.; LNB2CPP=.FALSE.; BLK2CP=.FALSE.
ENDIF
BLKM=LSLDR(IJ,IDIP).AND..NOT.LNB2CM.OR.BLK2CM
BLKP=LSLDR(IJ,IDIR).AND..NOT.LNB2CP.OR.BLK2CP
IF(BLKM.AND.BLKP) THEN
F(IDWDZ0+IJ)= 0.0
ELSE
IDZC=KZDZ+IZZ
IF(AVEVEL) THEN
IFC=IWG0(IPH)+IJ
DSC=0.5*F(IDZC); FC=F(IFC)
IF(.NOT.BLKM) THEN
IF(LNB2CM) THEN
DSM=F(PBC_3RDDST(3)+ITNBM)
ELSE
DSM=.5*F(IDZC-1)
ENDIF
FM=F(IWIZM0(IPH)+IJ)
ENDIF
RDSM= 1./(DSM+DSC+TINY)
c
IF(.NOT.BLKP) THEN
IF(LNB2CP) THEN
DSP=F(PBC_3RDDST(4)+ITNBP)
ELSE
DSP=.5*F(IDZC+1)
ENDIF
FP=F(IWIZP0(IPH)+IJ)
ENDIF
RDSP=1./(DSC+DSP+TINY)
c
IF(BLKM) THEN
FWLM=FC-(FP-FC)*DSC*RDSP
FWLP=(FC*DSP+FP*DSC)*RDSP
F(IDWDZ0+IJ)=(FWLP-FWLM)/(F(IDZC)+TINY)
ELSEIF(BLKP) THEN
FWLM=(FM*DSC+FC*DSM)*RDSM
FWLP=FC-(FM-FC)*DSC*RDSM
F(IDWDZ0+IJ)=(FWLP-FWLM)/(F(IDZC)+TINY)
ELSE !blkm & !blkp
FWLM=(FM*DSC+FC*DSM)*RDSM
FWLP=(FC*DSP+FP*DSC)*RDSP
F(IDWDZ0+IJ)=(FWLP-FWLM)/(F(IDZC)+TINY)
ENDIF
ELSE ! not.avevel
IWC=L0F(7+IPH-1)+IJ
IF(BLKM) THEN
IF(LSLDIR(IJ+IADW+IJPW(IDIR),IDIR).AND.
1 .NOT.LNB2CPP) THEN
F(IDWDZ0+IJ)=0.
ELSE
IF(LNB2CPP) THEN
FP=F(PBC_3RDVEL(6)+ITNBPP)
ELSE
FP=F(IWC+IJPS(IDIR))
ENDIF
F(IDWDZ0+IJ)=(FP-F(IWC))/(DLL(IJ,IZZ,IDIR)+TINY)
ENDIF
ELSEIF(BLKP) THEN
IF(LSLDIR(IJ+IADW+IJPW(IDIP),IDIP).AND.
1 .NOT.LNB2CMM) THEN
F(IDWDZ0+IJ)=0.
ELSE
IWM=IWC+IJPS(IDIP)
IF(LNB2CMM) THEN
FM=F(PBC_3RDVEL(5)+ITNBMM)
ELSE
FM=F(IWM+IJPS(IDIP))
ENDIF
F(IDWDZ0+IJ)=(F(IWM)-FM)/DLL(IJ,IZZ,IDIP)
ENDIF
ELSE
IF(LNB2CP) THEN
FP=F(PBC_3RDVEL(6)+ITNBP)
ELSE
FP=F(IWC)
ENDIF
IF(LNB2CM) THEN
FM=F(PBC_3RDVEL(5)+ITNBM)
ELSE
FM=F(IWC+IJPS(IDIP))
ENDIF
c F(IDWDZ0+IJ)=(FP-FM)/F(IDZC)
F(IDWDZ0+IJ)=(FP-FM)/F(KXYDZ+IJ)
ENDIF
ENDIF ! blkm or blkp
ENDIF ! blkm & blkp
ENDIF ! sld
ENDDO ! iy
ENDDO ! ix
ELSE
C.... Treatment of staggered velocities on non-BFC grids:
IF(AVEVEL) THEN
CALL FUVWDZ(IWIZM0(IPH),IWG0(IPH),IWIZP0(IPH),IDWDZ0,IZZ)
ELSE
CALL AVDVEL(5,6,L0F(7+IPH-1),IDWDZ0,IZZ)
ENDIF
ENDIF
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C.... AVRGVL averages staggered velocities U1,...,W2. Result is stored
C in temporary storage IUAV0. The later is used in calculations of
C the generation function for staggered solver on non-BFC geometry.
C The user might print-out averaged velocities by storing U1C,V1C,
C etc in Q1 (see GXGENF for details).
C.... IDIP is the direction index; which can be either IDIP=4 for U1;
C or IDIP=2 for V1; or IDIP=6 for W1.
C.... ISLAB is slab index; ISLAB=0 -averaging is for the current slab;
C ISLAB=1 -averaging is for the HIGH-slab.
C.... AVRGVL is called from GXGENF.
C
c
SUBROUTINE AVRGVL(IDIP,IUAV0,IPH,IZZ,ISLAB)
INCLUDE 'cmndmn'
INCLUDE 'farray'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'd_earth/pbcstr'
COMMON /IDATA/NX,NY,NZ,IDT1(117) /F0/LB1(109),KZXCY,LB2(194)
1 /GENI/ NXNY,NXM1NY,NXNYST,IG1(6),NFM,IG2(50)
1 /FNDRI/IDF(11),IJPW(6),IJPS(6),IDSG(6),IDR(6),JDR(6),KDR(6)
1 /FACES/L0FACE,L0FACZ /FACDIR/IBND(6),IFCP(6),IWLP(6)
LOGICAL SLD,LSLDR,BLKM,BLKP,XCYCZ,NEZ
COMMON/CCTDBL/LCUTDBL
LOGICAL LCUTDBL, LCUTCL
COMMON/INTDMN/IDMN,INTDM1(8)
COMMON/GEODMN0/I3DAEX,I3DANY,I3DAHZ,I3DVOL,I3DDXG,I3DDYG,I3DDZG,
1 I3DDX,I3DDY,I3DDZ,I2DAWB,I2DASB,I2DALB
LOGICAL LNB2CM,LNB2CP,BLK2CM,BLK2CP,LDO1,LDO2
C.... Preliminaries
IDIR=IDIP-1; L0FZST=L0FACZ
L0FACZ= L0FACZ + ISLAB*NXNYST
IF(IDIR==1) THEN
IUST0= L0F(V1+IPH-1) + ISLAB*NFM
ELSEIF(IDIR==3) THEN
IUST0= L0F(U1+IPH-1) + ISLAB*NFM
ELSE
IUST0= L0F(W1+IPH-1) + ISLAB*NFM
ENDIF
XCYCZ= IDIR==3 .AND. KZXCY/=0
IF(XCYCZ) XCYCZ= NEZ(F(KZXCY+IZZ))
IF(IDIR==1) THEN
L0AREA=I3DANY; ISH=1
ELSEIF(IDIR==3) THEN
L0AREA=I3DAEX; ISH=NY
ELSEIF(IDIR==5) THEN
L0AREA=I3DAHZ; ISH=NXNY
ENDIF
C.... Loop over the current slab
DO 20 IX= 1,NX
IADX= (IX-1)*NY
IF(XCYCZ .AND. IX==1) IJPS(4)= NXM1NY
DO 10 IY= 1,NY
IJ= IY+IADX
IF( SLD(IJ) ) THEN
F(IUAV0+IJ)= 0.0
ELSE
IUSTC= IUST0+IJ
BLKM = LSLDR(IJ,IDIP)
BLKP = LSLDR(IJ,IDIR)
IF(LCUTDBL) THEN
IJK=IJ+(IZZ-1)*NXNYST
IF(IDIR==1) THEN
LDO1=IY>1; LDO2=IY1; LDO2=IX1; LDO2=IZZ0.
ENDIF
BLK2CM=F(L0AREA+IJK-ISH)<=0..AND..NOT.LNB2CM
ELSE
LNB2CM=.FALSE.; BLK2CM=.FALSE.
ENDIF
c
IF(LDO2) THEN
IJKNBP=IJK+NY; LNB2CP=LCUTCL(IJKNBP)
IF(LNB2CP) THEN
IPBCNB=NINT(F(NPBIJK(IDMN)+IJKNBP))
ITNBP=ISHPB(IDMN,IPBCNB)
LNB2CP=F(PBC_CCTYP+ITNBP)==2..AND.
1 F(L0AREA+IJK)<=0..AND.F(PBC_3RDAR(IDIP)+ITNBP)>0.
ENDIF
BLK2CP=F(L0AREA+IJK)<=0..AND..NOT.LNB2CP
ELSE
LNB2CP=.FALSE.; BLK2CP=.FALSE.
ENDIF
c
BLKM=BLKM.OR.BLK2CM
BLKP=BLKP.OR.BLK2CP
c
IF(BLKM.AND.BLKP) THEN
F(IUAV0+IJ)=0.
ELSEIF(BLKM) THEN
IF(LNB2CP) THEN
F(IUAV0+IJ)=F(PBC_3RDVEL(IDIP)+ITNBP)
ELSE
F(IUAV0+IJ)= F(IUSTC)
ENDIF
ELSEIF(BLKP) THEN
IF(LNB2CM) THEN
F(IUAV0+IJ)=F(PBC_3RDVEL(IDIR)+ITNBM)
ELSE
F(IUAV0+IJ)=F(IUSTC+IJPS(IDIP))
ENDIF
ELSE
IF(LNB2CP) THEN
FP=F(PBC_3RDVEL(IDIP)+ITNBP)
ELSE
FP=F(IUSTC)
ENDIF
c
IF(LNB2CM) THEN
FM=F(PBC_3RDVEL(IDIR)+ITNBM)
ELSE
FM=F(IUSTC+IJPS(IDIP))
ENDIF
F(IUAV0+IJ)=0.5*(FP+FM)
ENDIF
ELSE
IF(BLKM .AND. BLKP) THEN
F(IUAV0+IJ)= 0.0
ELSEIF(BLKM) THEN
F(IUAV0+IJ)= F(IUSTC)
ELSEIF(BLKP) THEN
F(IUAV0+IJ)= F(IUSTC+IJPS(IDIP))
ELSE
F(IUAV0+IJ)= 0.5*( F(IUSTC) + F(IUSTC+IJPS(IDIP)) )
ENDIF
ENDIF
ENDIF
10 CONTINUE
IF(XCYCZ .AND. IX==1) IJPS(4)= -NY
20 CONTINUE
L0FACZ= L0FZST
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C.... GXKEGB is called from group 1, group 19 section 4, and group 13 of
C GREX3 by setting the coefficient & value to GRND4 in the COVAL
C statements for KE and EP (or OMEG) for the patch name KEBUOY.
C See the 'K-Epsilon turbulence model' and 'GRAVitational' entries
C in the PHOENICS Encyclopaedia for further details.
C
C.... BUOYA, BUOYB and BUOYC signify the resolutes of the gravitational
C acceleration in the cartesian coordinate direction XC, YC and ZC.
C
C.... The coefficient GCEB determines the magnitude of the buoyancy
C production term in the EP (or OMEG) equation. There is no "standard"
C value for GCEB, although it is recommended that it should be close to
C zero for stably-stratified flows, and close to unity for unstably-
C stratified flows. The default in PHOENICS-VR is to compute GCEB
C from the function: GCEB = tanh (v/U) where v is the velocity
C component parallel to the gravity vector, and U is the velocity
C component perpendicular to the gravity vector.
C This default setting is GCEB=1.0, but this may be overwritten by
C setting: SPEDAT(SET,KEBUOY,GCEB,R,0.2), say, in the Q1 file.
C
c
SUBROUTINE GXKEGB
INCLUDE 'farray'
INCLUDE 'satear'
INCLUDE 'grdear'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'grdbfc'
COMMON /INTDMN/IDMN,NUMDMN,NXMAX,NYMAX,NZMAX,NXYMAX,NXYZMX,NCLTOT
COMMON /GENI/ IG1(2),NXNYST,IG2(30),ILTLS,IG3(15),ITEM1,ITEM2,
1 IG4(4),IPRPS,IG5(4)
1 /F0/ LB1(109),KZXCY,LB2(194)
1 /FNDRL/NXNE1,NYNE1,NZNE1 /NAMFN/NAMFUN,NAMSUB
1 /GENL/ LGL1(14),debgtz,LGL2(45)
1 /FNDRI/IDF(11),IJPW(6),IJPS(6),IDSG(6),IDR(6),JDR(6),KDR(6)
1 /LRNTM2/LBWDIS,LBFMU,LBFONE,LBFTWO,LBREYN,LBREYT,LBEPKE
COMMON/TSKEMI/ LBKP,LBKT,LBET,LBVOSQ,LBOMEG
C!!! MRM 15.11.05 Introduce Automatic C3e function in GXKEGB
COMMON/BFICRT/IUCRT(2),IVCRT(2),IWCRT(2),ITMP(6)
DIMENSION GV(3),VF(3),VPG(3),VNG(3)
LOGICAL NXNE1,NYNE1,NZNE1,GRN,BOUSS,SLDEF,LGL1,debgtz,LGL2,
1 dbgloc,solu,solv,solw,VARC3E,XCYCZ,BLKM,BLKP,LSLDR,NEZ
CHARACTER*6 NAMFUN,NAMSUB
SAVE BOUSS,JSCAL,JDRDX,JDRDY,JDRDZ,JGENB,L0DRDX,L0DRDY,L0DRDZ,
C!!! MRM 15.11.05 Introduce Automatic C3e function in GXKEGB
1 L0GENB,L0SCAL,GCEB,RECPRT,SLDEF,L0C3EB0,SOLU,SOLV,SOLW,
1 JC3EB,VARC3E,C3ERLX
COMMON/USPGLOB/USPG, USPDBG, USPSCHEMA,UseVTK
LOGICAL(4) USPG,USPDBG,UseVTK
INTEGER(4) USPSCHEMA
C
NAMSUB='GXKEGB'
DBGLOC=DEBGTZ
IF(INDVAR>0) DBGLOC=DBGLOC.AND.DBGPHI(INDVAR)
IF(FLAG.OR.DBGLOC) CALL BANNER(1,NAMSUB,200918)
IF(CCM) THEN
CALL UNKEGB
RETURN
ENDIF
IF(IGR==1 .AND. ISC==2) THEN
C-----------------------------------------------------------------------
C CALL FROM GR.1; SEC.2 TO PREPARE CALCULATIONS
C-----------------------------------------------------------------------
C.... THE BOUSSINESQ APPROXIMATION IS EMPLOYED WHEN RHO1 IS CONSTANT.
BOUSS= .NOT.GRN(RHO1)
SLDEF= .FALSE.
IF(BOUSS) THEN
IF(SOLVE(ITEM1)) THEN
JSCAL= ITEM1
SLDEF= .TRUE.
ELSEIF(SOLVE(H1)) THEN
JSCAL= H1
ELSE
CALL WRYT40('GXKEGB REQUIRES TEM1 OR H1 TO BE SOLVED ')
CALL WRYT40('FOR THE CURRENT GRAVITY &/OR RHO1 MODELS')
CALL WRIT40('GXKEGB REQUIRES TEM1 OR H1 TO BE SOLVED ')
CALL WRIT40('FOR THE CURRENT GRAVITY &/OR RHO1 MODELS')
CALL SET_ERR(222,'ERROR. SEE RESULT FILE.',1)
C CALL EAROUT(1)
ENDIF
RECPRT= 1./PRT(JSCAL)
ELSE
IF(DEN1<=0) THEN
CALL WRIT40('GXKEGB REQUIRES STORE(DEN1) IN Q1! ')
CALL WRYT40('GXKEGB REQUIRES STORE(DEN1) IN Q1! ')
CALL SET_ERR(223,'ERROR. GXKEGB REQUIRES STORE(DEN1) IN Q1.'
* ,1)
C CALL EAROUT(1)
ENDIF
JSCAL = DEN1
RECPRT= 1./PRT(H1)
ENDIF
IF(USPG)THEN
JGENB=LBNAME('GENB')
IF(JGENB==0) CALL AddStored('GENB')
GCEB=1.0
CALL GETSDR('KEBUOY','GCEB',GCEB)
VARC3E=.FALSE.
IF(GCEB<0.0) VARC3E=.TRUE.
IF(VARC3E) THEN
JC3EB=LBNAME('C3EB')
IF(JC3EB==0) CALL AddStored('C3EB')
ENDIF
RETURN
ENDIF
C.... ALLOCATE STORAGE FOR SCALAR DERIVATIVES (LD1,LD2,LD3 ARE IN USE):
JDRDX=LBNAME('DRDX'); JDRDY=LBNAME('DRDY')
JDRDZ=LBNAME('DRDZ'); JGENB=LBNAME('GENB')
L0DRDX=L0F(LD1); L0DRDY=L0F(LD2); L0DRDZ=L0F(LD3)
C.... IF A 3D STORE HAS NOT BE PROVIDED BY STORE(GENB) IN Q1, CREATE
C A SLAB-WISE STORE AND INITIALISE CONTENTS TO ZERO.
IF(JGENB==0) CALL GXMAKE0(L0GENB,NXYMAX,'GENB')
C
GCEB=1.0
CALL GETSDR('KEBUOY','GCEB',GCEB)
C!!! MRM 19.10.05 INTRODUCE HENKES C3E FUNCTION
VARC3E=.FALSE.
IF(GCEB<0.0) VARC3E=.TRUE.
IF(VARC3E) THEN
JC3EB=LBNAME('C3EB')
IF(JC3EB==0) CALL GXMAKE0(L0C3EB0,NXYZMX,'C3EB')
C3ERLX=0.3
CALL GETSDR('KEBUOY','C3ERLX',C3ERLX)
SOLU=SOLVE(U1); SOLV=SOLVE(V1); SOLW=SOLVE(W1)
ENDIF
CALL WRITBL
CALL WRITST
CALL WRIT40('EP buoyancy production coefficient C3EB ')
IF(VARC3E) THEN
CALL WRIT40('Variable C3EB=tanh(V/U) ')
CALL WRIT40('Linear under-relaxation factor for C3EB:')
CALL WRIT1R('KE-C3RLX',C3ERLX)
ELSE
CALL WRIT1R('K-E GCEB',GCEB)
CALL WRIT40('The coefficient C3EB may be changed from')
CALL WRIT40('the Q1 file by setting: ')
CALL WRIT40('SPEDAT(SET,KEBUOY,GCEB,R,value)')
CALL WRIT40('A negative value of GCEB activates a ')
CALL WRIT40('variable C3EB=tanh(V/U). ')
ENDIF
CALL WRITST
ELSEIF(IGR==19 .AND. ISC==4) THEN
C-----------------------------------------------------------------------
C Call from Gr.19; Sec.4 to calculate Gb term.
C-----------------------------------------------------------------------
C.... If DRDX,...,GENB are stored in Q1, get slab-wise addresses:
IF(JDRDX/=0) L0DRDX= L0F(JDRDX)
IF(JDRDY/=0) L0DRDY= L0F(JDRDY)
IF(JDRDZ/=0) L0DRDZ= L0F(JDRDZ)
IF(JGENB/=0) L0GENB= L0F(JGENB)
C.... Calculate derivatives of JSCAL
L0SCAL= L0F(JSCAL)
CALL ZERNM3(L0DRDX,L0DRDY,L0DRDZ,NXNYST)
IF(NXNE1) CALL FNDFDX(L0SCAL,L0DRDX,IZSTEP,SLDEF)
IF(NYNE1) CALL FNDFDY(L0SCAL,L0DRDY,IZSTEP,SLDEF)
IF(NZNE1.AND..NOT.PARAB) CALL FNDFDZ(L0SCAL,L0DRDZ,IZSTEP,SLDEF)
IF(CARTES .OR. BFC) THEN
C... BFC or Cartesian geometries:
DO 10 IJ= 1,NXNYST
GDRDX= F(L0DRDX+IJ)
GDRDY= F(L0DRDY+IJ)
GDRDZ= F(L0DRDZ+IJ)
F(L0GENB+IJ)= BUOYA*GDRDX + BUOYB*GDRDY + BUOYC*GDRDZ
10 CONTINUE
ELSE
C.... Polar-cylindrical geometry:
L0XG= L0F(LXXG)
DO IX= 1,NX
IADX = (IX-1)*NY
ANGLE= F(L0XG+IX)
SINA = SIN(ANGLE); COSA = COS(ANGLE)
DO IY= 1,NY
IJ= IY+IADX
GDRDX= F(L0DRDX+IJ)
GDRDY= F(L0DRDY+IJ)
GDRDZ= F(L0DRDZ+IJ)
F(L0GENB+IJ)= BUOYA*( COSA*GDRDX+SINA*GDRDY)
1 + BUOYB*(-SINA*GDRDX+COSA*GDRDY) + BUOYC*GDRDZ
ENDDO
ENDDO
ENDIF
C.... The volumetric production of K due to buoyancy forces
L0ENUT= L0F(VIST)
DO IJ= 1,NXNYST
F(L0GENB+IJ)= -F(L0ENUT+IJ)*F(L0GENB+IJ)*RECPRT
ENDDO
IF(BOUSS) THEN
C.... Put expansion coefficient into L0GENB:
C.... FN46(Y,X,A,B) ........... Y = Y*(A+B*X):
CALL FN46(-L0GENB,LDVO1,0.0,-1.0)
C.... Divide by Cp when enthalpy is the solved for variable
C.... FN27(Y,A) ............... Y = Y / X
IF(JSCAL==H1) CALL FN27(-L0GENB,LCP1)
ELSE
CALL FN27(-L0GENB,DEN1)
ENDIF
C
C!!! MRM 15.11.05 Introduce Automatic C3e function in GXKEGB
C.... Compute C3E from velocity function
IF(VARC3E) THEN
C.... Unit gravity vector
CALL VECTOR(GV,BUOYA,BUOYB,BUOYC)
CALL NORM(GV)
IF(BFC) THEN
IF(SOLU) L0U1=L0F(IUCRT(1))
IF(SOLV) L0V1=L0F(IVCRT(1))
IF(SOLW) L0W1=L0F(IWCRT(1))
ELSE
IF(SOLU) L0U1=L0F(U1)
IF(SOLV) L0V1=L0F(V1)
IF(SOLW) L0W1=L0F(W1)
ENDIF
IF(JC3EB/=0) THEN
L0C3EB = L0F(JC3EB)
ELSE
L0C3EB=L0C3EB0+(IZ-1)*NXNYST
ENDIF
L0XG= L0F(LXXG)
XCYCZ= SOLU .AND. KZXCY/=0
IF(XCYCZ) XCYCZ= NEZ(F(KZXCY+IZ))
DO IX= 1,NX
IADX = (IX-1)*NY
IF(XCYCZ .AND. IX==1) IJPS(4)= NXNYST-NX
DO IY= 1,NY
IJ= IY+IADX
C.... Cartesian velocity vector
CALL VECTOR(VF,0.0,0.0,0.0)
IF(SOLW) THEN
IDIR=5; IDIP=6
BLKM = LSLDR(IJ,IDIP); BLKP = LSLDR(IJ,IDIR)
IF(BLKM .AND. BLKP) THEN
VF(3)=0.0
ELSEIF(BLKM) THEN
VF(3)=F(L0W1+IJ)
ELSEIF(BLKP) THEN
VF(3)=F(L0W1+IJ+IJPS(IDIP))
ELSE
VF(3)=0.5*(F(L0W1+IJ)+F(L0W1+IJ+IJPS(IDIP)))
ENDIF
ENDIF
IF(SOLU) THEN
IDIR=3; IDIP=4
BLKM = LSLDR(IJ,IDIP); BLKP = LSLDR(IJ,IDIR)
IF(BLKM .AND. BLKP) THEN
VF(1)=0.0
ELSEIF(BLKM) THEN
VF(1)=F(L0U1+IJ)
ELSEIF(BLKP) THEN
VF(1)=F(L0U1+IJ+IJPS(IDIP))
ELSE
VF(1)=0.5*(F(L0U1+IJ)+F(L0U1+IJ+IJPS(IDIP)))
ENDIF
ENDIF
IF(SOLV) THEN
IDIR=1; IDIP=2
BLKM = LSLDR(IJ,IDIP); BLKP = LSLDR(IJ,IDIR)
IF(BLKM .AND. BLKP) THEN
VF(2)=0.0
ELSEIF(BLKM) THEN
VF(2)=F(L0V1+IJ)
ELSEIF(BLKP) THEN
VF(2)=F(L0V1+IJ+IJPS(IDIP))
ELSE
VF(2)=0.5*(F(L0V1+IJ)+F(L0V1+IJ+IJPS(IDIP)))
ENDIF
ENDIF
IF(.NOT.CARTES) THEN
UVEL=VF(1); VVEL=VF(2)
ANGLE= F(L0XG+IX)
SINA = SIN(ANGLE); COSA = COS(ANGLE)
IF(SOLU) VF(1)= COSA*UVEL + SINA*VVEL
IF(SOLV) VF(2)= -SINA*UVEL + COSA*VVEL
ENDIF
C.... velocity vectors parallel and normal to g
VFCPG = DOT(VF,GV)
CALL MULVEC(GV,VFCPG,VPG)
CALL SUBVEC(VF,VPG,VNG)
VFCNG = VECMAG(VNG)
VRAT = ABS(VFCPG/(VFCNG + TINY))
C3ENEW = TANH(VRAT)
F(L0C3EB+IJ)=C3ENEW*C3ERLX+(1.-C3ERLX)*F(L0C3EB+IJ)
ENDDO
IF(XCYCZ .AND. IX==1) IJPS(4)= -NY
ENDDO
ENDIF
C
ELSEIF(IGR==13 .AND. ISC==5) THEN
C-----------------------------------------------------------------------
C Call from Gr.13; Sec.5 to process CO=GRND4.
C-----------------------------------------------------------------------
L0CO=L0F(CO); L0KE=L0F(KE)
IF(INDVAR==KE) THEN
DO 40 IJ= 1,NXNYST
40 F(L0CO+IJ)= AMAX1( FIXFLU, -F(L0GENB+IJ)/(F(L0KE+IJ)+TINY) )
ELSEIF(INDVAR==EP.OR.INDVAR==LBOMEG) THEN
C!!! MRM 15.11.05 Introduce Automatic C3e function in GXKEGB
IF(JC3EB/=0) THEN
L0C3EB = L0F(JC3EB)
ELSE
L0C3EB=L0C3EB0+(IZ-1)*NXNYST
ENDIF
DO 50 IJ= 1,NXNYST
C!!! MRM 15.11.05 Introduce Automatic C3e function in GXKEGB
GC3EB=GCEB
IF(VARC3E) GC3EB = F(L0C3EB+IJ)
50 F(L0CO+IJ)=AMAX1(FIXFLU,-GC3EB*F(L0GENB+IJ)/(F(L0KE+IJ)+TINY))
ENDIF
ELSEIF(IGR==13 .AND. ISC==16) THEN
C-----------------------------------------------------------------------
C Call from Gr.13; Sec.16 to process VAL=GRND4.
C-----------------------------------------------------------------------
L0VAL= L0F(VAL)
RFF = 1./FIXFLU
IF(INDVAR==KE) THEN
DO 60 IJ= 1,NXNYST
60 F(L0VAL+IJ)= AMAX1( 0.0, RFF*F(L0GENB+IJ) )
ELSEIF(INDVAR==EP.OR.INDVAR==LBOMEG) THEN
L0KE=L0F(KE)
IF(INDVAR==EP) THEN
L0EP=L0F(EP)
IF(LBEPKE/=0) L0EPKE=L0F(LBEPKE)
ELSE
L0OMEG=L0F(LBOMEG)
ENDIF
C!!! MRM 15.11.05 Introduce Automatic C3e function in GXKEGB
IF(JC3EB/=0) THEN
L0C3EB = L0F(JC3EB)
ELSE
L0C3EB=L0C3EB0+(IZ-1)*NXNYST
ENDIF
DO 70 IJ= 1,NXNYST
IF(INDVAR==EP) THEN
GEPDK= F(L0EP+IJ)/(F(L0KE+IJ)+TINY)
IF(LBEPKE/=0) F(L0EPKE+IJ) = GEPDK
ELSE
GEPDK= F(L0OMEG+IJ)/(F(L0KE+IJ)+TINY)
ENDIF
C!!! MRM 15.11.05 Introduce Automatic C3e function in GXKEGB
GC3EB=GCEB
IF(VARC3E) GC3EB = F(L0C3EB+IJ)
F(L0VAL+IJ)= AMAX1( RFF*GC3EB*GEPDK*F(L0GENB+IJ), 0.0 )
70 CONTINUE
ENDIF
ENDIF
namsub='gxkegb'
if(flag.or.dbgloc) call banner(2,namsub,200918)
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C.... AXIBFC checks the BFC-grid and returns AXISYM= T if the grid is
C axisymmetric.
c
SUBROUTINE AXIBFC(AXISYM)
COMMON /IDATA/ NX,NY,NZ,IDT1(117) /RDATA/TINY,RDT1(84)
1 /XYZCRN/X1,X2,X3,X4,X5,X6,X7,X8,Y1,Y2,Y3,Y4,Y5,Y6,Y7,Y8,
1 Z1,Z2,Z3,Z4,Z5,Z6,Z7,Z8
LOGICAL AXISYM
C
IF(NY<=2) THEN
AXISYM= .FALSE.
ELSE
AXISYM= .TRUE.
DO 10 IZZ= 1,NZ
CALL CELCOR(1,2,IZZ)
A1= ABS(2.*(X4-X1)/(Y1+Y4+TINY))
A2= ABS(2.*(X2-X3)/(Y2+Y3+TINY))
A3= ABS(2.*(X5-X8)/(Y5+Y8+TINY))
A4= ABS(2.*(X6-X7)/(Y6+Y7+TINY))
AA1= 0.25*(A1+A2+A3+A4)
CALL CELCOR(1,NY,IZZ)
A1= ABS(2.*(X4-X1)/(Y1+Y4+TINY))
A2= ABS(2.*(X2-X3)/(Y2+Y3+TINY))
A3= ABS(2.*(X5-X8)/(Y5+Y8+TINY))
A4= ABS(2.*(X6-X7)/(Y6+Y7+TINY))
AA2= 0.25*(A1+A2+A3+A4)
AXISYM= AXISYM .AND. ABS(AA1-AA2)<0.01*AA2
10 CONTINUE
ENDIF
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c