C
C General function of the subroutine
C ----------------------------------
C
C GXHOL can be used for the simulation of free-surface flows,
C provided that there is no "overturning" of the surface, so
C that the surface height is a single-valued function of the
C other two coordinate directions.
C
C The "height" may be measured in the x, y or z direction,
C the selection being made by setting IHOLA to 1, 2 or 3 in the
C Q1 file, respectively. Negative IHOLAs mean that the ''height''
C is measured with reference to one of the 'high' boundaries of
C the flow domain.
C
C Blockages may not be present in the flow domain.
C
C Grids may be cartesian, polar-cylindrical or body-fitted.
C
C Functions of particular parts of the subroutine
C -----------------------------------------------
C
C The following Groups are visited:-
C
C * Group 1, in order to allocate storage and to make certain
C once-for-all settings
C
C * Group 8, in order to compute the column-wise inflows and
C outflows of liquid, from which the heights of
C liquid in the columns are computed.
C
C * Group 9, in order to compute the volume fraction of liquid
C from which material properties are deduced.
C
C * Group 13, in order to specify outflow boundary conditions
C for flows where the body force is normal to the
C main flow direction.
C
C * Group 19, in order to initialize and update various arrays
C which are used in other parts of this subroutine.
C Also updated are material properties.
C
C For more detailed description of HOL, users are referred to the
C PHOENICS Encyclopaedia accessible via POLIS.
C
C Definitions of FNs are given where they appear in the code. The
C FUNCTion entry of the Encyclopaedia provides full explanations.
C
C***********************************************************************
c
SUBROUTINE GXHOL
INCLUDE 'patcmn'
INCLUDE 'farray'
INCLUDE 'satear'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'grdear'
C
PARAMETER (NLG=20, NIG=20, NRG=100, NCG=10)
COMMON/LGRND/LG(NLG)/IGRND/IG(NIG)/RGRND/RG(NRG)/CGRND/CG(NCG)
COMMON/GENI/NXNY,IGFIL0(1),NXNYST,IGFIL1(52),IPRPS,IGFIL2(4)
COMMON/NAMFN/NAMFUN,NAMSUB
LOGICAL LG,VISITB,VISITX,VISITY,VISITZ,XDIREC,YDIREC,ZDIREC
CHARACTER*6 NAMFUN,NAMSUB
CHARACTER*4 CG
SAVE VISITX,VISITY,VISITZ,VISITB,XDIREC,YDIREC,ZDIREC,IFIRST,N,
1 IFREQU,GRELAX,RLXFAC,L0FLUX,L0MSIN,L0MOUT,LNOLDP,LNVFOL,ICHK,
1 RHOLIQ,RHOGAS,L0TMOL,L0MOLO,L0LMOL,L0INOU,L0TLIQ,L0VOLL,GGDT,
1 L0VFOL,L0VOFH,LCVOL,IZP1,L0EASP
C
NAMSUB='GXHOL'
if(flag) call banner(1,namsub,011116)
IXL=IABS(IXL)
C>>>>>>>>>>>>>>>>>>>>> Group 1 ------- Section 1 <<<<<<<<<<<<<<<<<<<<<
C Storage allocations for the slab-wise variables.
C
C Note that RHO1=GRND10 is being used as the activator of HOL
C Set local variables with meaningful names to values which
C have been set in Q1.
C.......................................................................
IF(IGR.EQ.1.AND.ISC.EQ.1) THEN
IF(IPRPS.NE.0) THEN
RLXFAC=MIN(1.0,ABS(DTFALS(P2)))
XDIREC=IABS(IHOLA).EQ.1
YDIREC=IABS(IHOLA).EQ.2
ZDIREC=IABS(IHOLA).EQ.3
RUPLIM = 1.0 - RLOLIM
IF(.NOT.STEADY) THEN
IFIRST = 1
IFREQU = 1
RLXFAC=1.0
ELSE
IFIRST=-1; IFREQU=-1
CALL GETSDI('HOL','IFIRST',IFIRST)
IF(IFIRST.EQ.-1) IFIRST=IDISPB
CALL GETSDI('HOL','IFREQU',IFREQU)
IF(IFREQU.EQ.-1) IFREQU=IDISPC
IF(IFIRST.EQ.0) IFIRST=5
IF(IFREQU.EQ.0) IFREQU=5
ENDIF
GRELAX = RLXFAC
CALL WRITST
WRITE(LUPR1,'('' HOL Control Parameters'')')
WRITE(LUPR1,'('' First update sweep:'',I4)') IFIRST
WRITE(LUPR1,'('' Update frequency: '',I4)') IFREQU
WRITE(LUPR1,'('' Relaxation factor: '',1PE10.2)') GRELAX
CALL WRITST
ELSE
WRITE(LUPR1,*)'PRPS must be stored in Q1 for hol'
CALL SET_ERR(243,'Error. PRPS must be stored in Q1 for hol.',
1 1)
ENDIF
ELSEIF(IGR.EQ.1.AND.ISC.EQ.3) THEN
C>>>>>>>>>>>>>>>>>>>>> Group 1 ------- Section 3 <<<<<<<<<<<<<<<<<<<<<
C Provision of slab-wise storage
C TMOL - Total Mass Of Liquid in a cell column
C MOLO - TMOL at the previous time step/isweep
C LMOL - Lower Mass Of Liquid in a cell column
C ie TMOL below the interface cell
C INOU - Rate of change of TMOL through IN- and OUT-flows
C VOLL - Volume of cell of the lower IZ slab
C TLIQ - Total amount of liquid in a cell column possible
IF(XDIREC) THEN
N = NY*NZ
ELSEIF(YDIREC) THEN
N = NX*NZ
ELSEIF(ZDIREC) THEN
N = NY*NX
ELSE
WRITE(LUPR1,*) 'HOL direction not set in Q1'
CALL SET_ERR(244,'Error. HOL direction not set in Q1.',1)
ENDIF
CALL GXMAKE(L0TMOL,N,'TMOL')
CALL GXMAKE(L0MOLO,N,'MOLO')
CALL GXMAKE(L0LMOL,N,'LMOL')
CALL GXMAKE(L0INOU,N,'INOU')
CALL GXMAKE(L0TLIQ,N,'TLIQ')
CALL GXMAKE(L0EASP,NY*NX,'EASP')
IF(ZDIREC) CALL GXMAKE(L0VOLL,N,'VOLL')
C>>>>>>>>>>>>>>>>>>>>> Group 1 ------- Section 2 <<<<<<<<<<<<<<<<<<<<<
C Variable initializations.
C.......................................................................
ELSEIF(IGR.EQ.1.AND.ISC.EQ.2) THEN
L0FLUX = L0F(LD5)
L0MSIN = L0F(LMIN1)
L0MOUT = L0F(LMOUT1)
LNOLDP = LBNAME('OLDP')
LNVFOL = LBNAME('VFOL')
RHOLIQ = F(INDPRTB(IPRPSA,0)+1)
RHOGAS = F(INDPRTB(IPRPSB,0)+1)
CALL WRITST
WRITE(LUPR1,'('' HOL densities'')')
CALL WRIT2R('RHOliq ',RHOLIQ,'RHOgas ',RHOGAS)
CALL WRITST
C>>>>>>>>>>>>>>>>>>>>> Group 8 ------- Section 7 <<<<<<<<<<<<<<<<<<<<<
C Compute volumetric source for GALA. Note that an isentropic gas
C flow is assumed.
C.......................................................................
ELSEIF(IGR.EQ.8.AND.ISC.EQ.7) THEN
IF(STORE(LNOLDP)) THEN
L0SU = L0F(LSU)
L0DEN1=L0F(DEN1)
L0OLDN=L0F(OLD(DEN1))
DO 701 IY = 1 , NY
DO 701 IX = 1 , NX
I = IY+(IX-1)*NY
VOLSOR=( 1.0-F(L0VFOL+I) ) *
1 ( 1.0/F(L0DEN1+I) - 1.0/F(L0OLDN+I) ) / DT
F(L0SU+I) = F(L0SU+I) + VOLSOR
701 CONTINUE
ENDIF
C>>>>>>>>>>>>>>>>>>>>> Group 8 ------- Section 8 <<<<<<<<<<<<<<<<<<<<<
C Compute inflows & outflows of liquid to a cell column
C at the current slab, so as to be able later to obtain
C the total mass of liquid in the cell column
C
C F(L0FLUX+...=the algebraic volumetric convective flux
c calculated inside EARTH
C
C Logicals VISITX, VISITY and VISITZ are used to ensure
c that the sum is calculated only once for each cell in
c each slab. ICHK serves a similar purpose
C
C INDVAR= R1 means that the 1st-phase convection fluxes
C are being computed
C.......................................................................
ELSEIF(IGR.EQ.8.AND.ISC.EQ.8) THEN
IF(INDVAR.NE.R1.OR.ICHK.GT.NZ) RETURN
C NDIRCA=1 means that the north direction is in question.
IF(NDIRCA.EQ.1 .AND. VISITY) THEN
DO 801 IY = 1,NY-1
CDIR$ IVDEP
DO 801 IX = 1,NX
I = IY+(IX-1)*NY
C Find the "upwind" volume fraction of liquid and
C multiply it by the algebraic flux times the density
C times (for unsteady flows) the time step.
IF(F(L0FLUX+I) .GE. 0.0) THEN
YFLUXN = GGDT*RHOLIQ*F(L0VFOL+I )*F(L0FLUX+I)
ELSE
YFLUXN = GGDT*RHOLIQ*F(L0VFOL+I+ 1)*F(L0FLUX+I)
ENDIF
C Add and subtract YFLUXN appropriately to and from
C the contents of the local and the neighbouring cells.
IF(ZDIREC) THEN
F(L0INOU+I ) = F(L0INOU+I ) - YFLUXN
F(L0INOU+I+ 1) = F(L0INOU+I+1) + YFLUXN
ELSEIF(XDIREC) THEN
ID=IY+(IZ-1)*NY
F(L0INOU+ID ) = F(L0INOU+ID ) - YFLUXN
F(L0INOU+ID+1) = F(L0INOU+ID+1) + YFLUXN
ENDIF
C This is all that needs to be done if solids are not present.
801 CONTINUE
C The next statement ensures that only one visit is made to
C the above sequence in each sweep.
VISITY = .FALSE.
C NDIRCA = 3 means that the east direction is in question
ELSEIF(NDIRCA.EQ.3.AND.VISITX) THEN
C Here is done for the x-direction what is done above for
C the y-direction
CDIR$ IVDEP
DO 802 IX = 1,NX-1
DO 802 IY = 1,NY
I = IY+(IX-1)*NY
IF(F(L0FLUX+I) .GE. 0.0) THEN
XFLUXE = GGDT*RHOLIQ*F(L0VFOL+I )*F(L0FLUX+I)
ELSE
XFLUXE = GGDT*RHOLIQ*F(L0VFOL+I+NY)*F(L0FLUX+I)
ENDIF
IF(ZDIREC) THEN
F(L0INOU+I ) = F(L0INOU+I ) - XFLUXE
F(L0INOU+I+NY) = F(L0INOU+I+NY) + XFLUXE
ELSEIF(YDIREC) THEN
ID = IX+(IZ-1)*NX
F(L0INOU+ID ) = F(L0INOU+ID ) - XFLUXE
F(L0INOU+ID+ 1) = F(L0INOU+ID+ 1) + XFLUXE
ENDIF
802 CONTINUE
VISITX = .FALSE.
C NDIRCA=5 means that the high direction is in question
ELSEIF(VISITZ .AND. NDIRCA.EQ.5 .AND.
1 .NOT.ZDIREC .and. iz.lt.nz) THEN
C Here is done for the z-direction what is done above
C for the x- and y-directions.
DO 803 IX = 1,NX
CDIR$ IVDEP
DO 803 IY = 1,NY
I = IY+(IX-1)*NY
IF(F(L0FLUX+I) .GE. 0.0) THEN
ZFLUXH = GGDT*RHOLIQ*F(L0VFOL+I)*F(L0FLUX+I)
ELSE
ZFLUXH = GGDT*RHOLIQ*F(L0VOFH+I)*F(L0FLUX+I)
ENDIF
IF(YDIREC) THEN
ID = IX+(IZ-1)*NX
F(L0INOU+ID ) = F(L0INOU+ID ) - ZFLUXH
F(L0INOU+ID+NX) = F(L0INOU+ID+NX) + ZFLUXH
ELSE
ID = IY+(IZ-1)*NY
F(L0INOU+ID ) = F(L0INOU+ID ) - ZFLUXH
F(L0INOU+ID+NY) = F(L0INOU+ID+NY) + ZFLUXH
ENDIF
803 CONTINUE
VISITZ = .FALSE.
ENDIF
C Add in inflow mass fluxes
IF(VISITB) THEN
DO 804 IR=1,NUMPAT
CALL GETPTC(NAMPAT(IR),GTYPE,IXF,IXL,IYF,IYL,IZF,IZL,JF,JL)
IF(NINT(GTYPE).GT.15.AND.NINT(GTYPE).NE.23) GO TO 804
CALL GETCOV(NAMPAT(IR),LNVFOL,GCO,GVAL)
IF(GVAL.GT.0.0) THEN
CONST=GGDT*(1.0-RHOGAS*GVAL)/(1.0-RHOGAS/RHOLIQ)
IF(ZDIREC) THEN
CALL FN34(-L0INOU,LMIN1,CONST)
ELSE
DO 805 IX = IXF,IXL
CDIR$ IVDEP
DO 805 IY = IYF,IYL
I = IY+(IX-1)*NY
IF(YDIREC) THEN
ID = IX+(IZ-1)*NX
F(L0INOU+ID)=F(L0INOU+ID)+F(L0MSIN+I)*CONST
ELSE
ID = IY+(IZ-1)*NY
F(L0INOU+ID)=F(L0INOU+ID)+F(L0MSIN+I)*CONST
ENDIF
805 CONTINUE
ENDIF
ENDIF
c.... reset ixl and iyl to ensure whole-slab print-out
IXL=NX
IYL=NY
804 CONTINUE
C subtract outflow mass fluxes
DO 806 IX = 1,NX
CDIR$ IVDEP
DO 806 IY = 1,NY
I =IY+(IX-1)*NY
OMASS = GGDT*F(L0MOUT+I)*F(L0VFOL+I)
1 /((1.0-F(L0VFOL+I))*RHOGAS/RHOLIQ+F(L0VFOL+I))
IF(ZDIREC) THEN
F(L0INOU+I) = F(L0INOU+I) - OMASS
ELSEIF(YDIREC) THEN
ID = IX+(IZ-1)*NX
F(L0INOU+ID) = F(L0INOU+ID) - OMASS
ELSE
ID = IY+(IZ-1)*NY
F(L0INOU+ID) = F(L0INOU+ID) - OMASS
ENDIF
806 CONTINUE
VISITB=.FALSE.
ENDIF
C>>>>>>>>>>>>>>>>>>>>> GROUP 9 ------- SECTION 1 <<<<<<<<<<<<<<<<<<<<<
c coding provided for computing, slab by slab, the vfols.
c call to fn34 - lmol(of current slab)=lmol(of last slab)
c +rholiq*vol(of last slab)
c call to fn67 - calculate volume fraction of liquid from
c local density
c izp1 ensures that lmol is calculated once for each slab.
c for transient flow: both ifirst and ifrequ are set to 1.
c for steady problem: IDISPB=number of sweep at which vfol
c undergoes first update.
c IDISPC=number of sweep after which a
c new VFOL field is calculated( >=2 ).
c definition of fns
c fn1 (y,a): y = a
c fn34(y,x,a): y = y + a*x
c fn67(y,x1,x2,x3,b1,b2): y = (b1*x1 + b2*x2)/x3
c fn22(y,a): y = amax1(y,a)
c fn23(y,a): y = amin1(y,a)
c fn2 (y,x,a,b): y = a + b*x
c fn59(y,x): y = 0 where x = 0
c fn26(y,x): y = y * x
C.......................................................................
ELSEIF(IGR.EQ.9.AND.ISC.EQ.1) THEN
c height measured in z direction
IF(ZDIREC) THEN
IF(IZSTEP.EQ.1) THEN
c at the bottom, the mass of liquid in the cell is set to zero
CALL FN1 (-L0LMOL, 0.0)
IZP1 = 2
ELSEIF(IZP1.EQ.IZSTEP) THEN
IZP1=IZSTEP+1
c lmol added up slab by slab
CALL FN34(-L0LMOL,-L0VOLL,RHOLIQ)
ENDIF
c store lower slab cell volume
CALL FN0(-L0VOLL,VOL)
ENDIF
IF(ISWEEP.GT.FSWEEP.AND.ISWEEP.GT.IFIRST.AND.
1 MOD(ISWEEP-1,IFREQU).EQ.0) THEN
ccc vfol is calculated from (tmol - lmol) / (vol * rholiq)
LCVOL = L0F(VOL)
IF(ZDIREC) THEN
CONT=1.0/RHOLIQ
CALL FN22(-L0TMOL,0.0)
IF(IHOLA.LT.0) THEN
CALL FN34(-L0TMOL,-L0TLIQ,-1.0)
CALL FN23(-L0TMOL,0.0)
CALL FN34(-L0TMOL,-L0TLIQ,1.0)
ENDIF
CALL FN10(LNVFOL,-L0TLIQ,-L0TMOL,0.0, 1.0, -1.0)
CALL FN40(LNVFOL)
CALL FN67(LNVFOL, LNVFOL,-L0LMOL,VOL,CONT,-CONT)
ELSEIF(YDIREC) THEN
DO 9101 IX = 1,NX
ID=IX+(IZ-1)*NX
F(L0LMOL+ID) = 0.
IF(F(L0TMOL+ID).LT.0.0) F(L0TMOL+ID)=0.0
IF(IHOLA.LT.0.AND.F(L0TMOL+ID).GT.F(L0TLIQ+ID))
1 F(L0TMOL+ID)=F(L0TLIQ+ID)
DO 9101 IY = 1,NY
I =IY+(IX-1)*NY
IF(IY.GT.1) F(L0LMOL+ID)=
1 F(L0LMOL+ID)+RHOLIQ*F(LCVOL+I- 1)
F(L0VFOL+I)=(ABS(F(L0TLIQ+ID)-F(L0TMOL+ID))-
1 F(L0LMOL+ID))/(RHOLIQ*F(LCVOL+I))
9101 CONTINUE
ELSE
DO 9102 IY = 1,NY
ID=IY+(IZ-1)*NY
F(L0LMOL+ID) = 0.
IF(F(L0TMOL+ID).LT.0.0) F(L0TMOL+ID)=0.0
IF(IHOLA.LT.0.AND.F(L0TMOL+ID).GT.F(L0TLIQ+ID))
1 F(L0TMOL+ID)=F(L0TLIQ+ID)
CDIR$ IVDEP
DO 9102 IX = 1,NX
I =IY+(IX-1)*NY
IF(IX.GT.1) F(L0LMOL+ID)=
1 F(L0LMOL+ID)+RHOLIQ*F(LCVOL+I-NY)
F(L0VFOL+I)=(ABS(F(L0TLIQ+ID)-F(L0TMOL+ID))-
1 F(L0LMOL+ID))/(RHOLIQ*F(LCVOL+I))
9102 CONTINUE
ENDIF
ccc call fn67(lnvfol,-l0tmol,-l0lmol,vol,1.0/rholiq,-1.0/rholiq)
c limiters ensuring vfol stay within 0 and 1
c fn22(y,a) .... y = amax1(y,a)
c fn23(y,a) .... y = amin1(y,a)
CALL FN22(LNVFOL,0.0)
CALL FN23(LNVFOL,1.0)
IF(IHOLA.LT.0) CALL FN2(LNVFOL,LNVFOL,1.0,-1.0)
ENDIF
c
c Note that fluid properties are calculated in Section 3, Group 19.
C.......................................................................
C>>>>>>>>>>>>>>>>>>>>> GROUP 13 ------- SECTION 1 <<<<<<<<<<<<<<<<<<<<<
c the following entry deals with the situation where an
c outflow boundary involves both fluids.
C.......................................................................
ELSEIF(IGR.EQ.13.AND.ISC.EQ.1) THEN
IF(INDVAR.EQ.P1) CALL FN0(CO,DEN1)
C>>>>>>>>>>>>>>>>>>>>> GROUP 13 ------- SECTION 12 <<<<<<<<<<<<<<<<<<<<<
c the following type of exit boundary treatment is usually
c required for specifying mass outflow when gravity acts at an
c angle to the main local flow direction. ship simulation is a
c case in point.
c definition of fn
c fn21(y,x1,x2,a,b) : y = a + b*x1*x2
C.......................................................................
ELSEIF(IGR.EQ.13.AND.ISC.EQ.12) THEN
IF(INDVAR.EQ.P1) THEN
CALL GETPTC(NPATCH,TYPE,IXF,IXL,IYF,IYL,J,J,J,J)
ITYPE=TYPE+0.1
c outflow through 'e' boundary
IF(ITYPE.EQ.2) THEN
INDVEL=WEST(U1)
ELSEIF(ITYPE.EQ.3) THEN
INDVEL=U1
ELSEIF(ITYPE.EQ.4) THEN
INDVEL=SOUTH(V1)
ELSEIF(ITYPE.EQ.5) THEN
INDVEL=V1
ELSEIF(ITYPE.EQ.6) THEN
INDVEL=LOW(W1)
ELSEIF(ITYPE.EQ.7) THEN
INDVEL=W1
ENDIF
CALL FN21(VAL,INDVEL,DEN1,0.0,-1.0)
ENDIF
C>>>>>>>>>>>>>>>>>>>>> GROUP 19 ------- SECTION 1 <<<<<<<<<<<<<<<<<<<<<
C * ------------------- SECTION 1 ---- Start of time step.
ELSEIF(IGR.EQ.19.AND.ISC.EQ.1) THEN
ICHK = 0
GGDT = DT
IF(STEADY) GGDT=1.0
IF(ISTEP.EQ.FSTEP .AND. ZDIREC) THEN
C Both "old mass of liquid" and tliq are set to zero
CALL FN1(-L0MOLO,0.0)
CALL FN1(-L0TLIQ,0.0)
ENDIF
C>>>>>>>>>>>>>>>>>>>>> Group 19 ------- Section 2 <<<<<<<<<<<<<<<<<<<<<
C * ------------------- SECTION 2 ---- Start of sweep.
ELSEIF(IGR.EQ.19.AND.ISC.EQ.2) THEN
IF(ICHK.Ge.NZ) RETURN
C Clear content of INOU at start of sweep, for ZDIREC only
IF(ZDIREC) THEN
IF(ISWEEP.LT.LSWEEP) CALL FN1(-L0INOU,0.0)
C Give MOLO value of TMOL for the previous time step
IF(ISWEEP.EQ.1.AND.ISTEP.GT.FSTEP) CALL FN0(-L0MOLO,-L0TMOL)
C MOLO after IFIRST sweeps and then every IFREQU sweeps
IF(STEADY.AND.ISWEEP.GT.FSWEEP.AND.ISWEEP.GT.IFIRST.
1 AND.MOD(ISWEEP-1,IFREQU).EQ.0)
1 CALL FN0(-L0MOLO,-L0TMOL)
ENDIF
IF(STEADY) THEN
IF(ISWEEP.GE.IFIRST.AND.MOD(ISWEEP,IFREQU).EQ.0) THEN
RLXFAC = GRELAX
ELSE
RLXFAC = 0.0
ENDIF
ENDIF
C>>>>>>>>>>>>>>>>>>>>> Group 19 ------- Section 3 <<<<<<<<<<<<<<<<<<<<<
C * ------------------- SECTION 3 ---- Start of iz slab.
C
ELSEIF(IGR.EQ.19.AND.ISC.EQ.3) THEN
IF(ICHK.GE.NZ) RETURN
C Store old P1
IF(STORE(LNOLDP).AND.ISWEEP.EQ.1) CALL FN0(LNOLDP,P1)
C Calculate fluid properties
CALL FNPRPS(IPRPS,LNVFOL,RLOLIM,RUPLIM,IPRPSB,IPRPSA)
C Set logicals to .TRUE., to ensure that the appropriate
C code segments are visited once for this slab
VISITX=.TRUE.
VISITY=.TRUE.
VISITZ=.TRUE.
VISITB=.TRUE.
C Find F-array pointer for the volume fraction of liquid
L0VFOL = L0F(LNVFOL)
C Store old values of VFOL, for use in under-relaxation
IF(LNVFOL.NE.0) CALL FN0(-L0EASP,LNVFOL)
IF(ZDIREC) RETURN
IF(NZ.GT.1) L0VOFH = L0F(HIGH(LNVFOL))
IF(YDIREC) THEN
CDIR$ IVDEP
DO 1931 IX = 1,NX
F(L0INOU+IX) = 0.0
IF(IZ.LT.NZ) F(L0INOU+IX+IZ*NX) = 0.0
1931 CONTINUE
ELSEIF(XDIREC) THEN
CDIR$ IVDEP
DO 1932 IY = 1,NY
F(L0INOU+IY) = 0.0
IF(IZ.LT.NZ) F(L0INOU+IY+IZ*NY) = 0.0
1932 CONTINUE
ENDIF
C>>>>>>>>>>>>>>>>>>>>> Group 19 ------- Section 4 start of iteration loop
C.......................................................................
ELSEIF(IGR.EQ.19.AND.ISC.EQ.4) THEN
C Total amount of liquid in each z-directed cell column
C at the beginning of the calculation
ICHK=ICHK+1
LCVOL = L0F(VOL)
IF(ISTEP.GT.FSTEP.OR.ISWEEP.GT.FSWEEP.OR.ICHK.GT.NZ) RETURN
IF(ZDIREC) THEN
CALL FN53(-L0MOLO,LNVFOL,VOL,RHOLIQ)
C TMOL is MOLO for 1st time step
CALL FN0 (-L0TMOL,-L0MOLO)
IF(IHOLA.LT.0) CALL FN34(-L0TLIQ,VOL,RHOLIQ)
ELSEIF(YDIREC) THEN
DO 1941 IX = 1,NX
ID=IX+(IZ-1)*NX
F(L0TMOL+ID) = 0.
F(L0TLIQ+ID) = 0.
DO 1941 IY = 1,NY
I =IY+(IX-1)*NY
IF(IHOLA.LT.0) F(L0TLIQ+ID)=F(L0TLIQ+ID)+RHOLIQ*F(LCVOL+I)
F(L0TMOL+ID) = F(L0TMOL+ID)+RHOLIQ*F(LCVOL+I)*F(L0VFOL+I)
1941 CONTINUE
ELSE
CDIR$ IVDEP
DO 1942 IY=1,NY
ID=IY+(IZ-1)*NY
F(L0TMOL+ID) = 0.0
F(L0TLIQ+ID) = 0.0
DO 1942 IX = 1,NX
I = IY+(IX-1)*NY
IF(IHOLA.LT.0) F(L0TLIQ+ID)=F(L0TLIQ+ID)+RHOLIQ*F(LCVOL+I)
F(L0TMOL+ID) = F(L0TMOL+ID)+RHOLIQ*F(LCVOL+I)*F(L0VFOL+I)
1942 CONTINUE
ENDIF
C>>>>>>>>>>>>>>>>>>>>> Group 19 ------- Section 6 finish of iz slab
C Tmol computed as molo + inou
C Under-relax vfol
C.......................................................................
ELSEIF(IGR.EQ.19.AND.ISC.EQ.6) THEN
IF(ICHK.GT.NZ) RETURN
IF(YDIREC) THEN
CDIR$ IVDEP
DO 1961 IX = 1 , NX
ID=IX+(IZ-1)*NX
IF(ISWEEP.EQ.IFIRST) THEN
F(L0MOLO+ID) = F(L0TMOL+ID)
ELSEIF(STEADY.AND.ISWEEP.GE.FSWEEP.AND.ISWEEP.GE.IFIRST.
1 AND.MOD(ISWEEP,IFREQU).EQ.0) THEN
F(L0MOLO+ID) = F(L0TMOL+ID)
ENDIF
F(L0TMOL+ID) = F(L0MOLO+ID) + RLXFAC*F(L0INOU+ID)
1961 CONTINUE
ELSEIF(XDIREC) THEN
CDIR$ IVDEP
DO 1962 IY = 1 , NY
ID=IY+(IZ-1)*NY
IF(ISWEEP.EQ.IFIRST) THEN
F(L0MOLO+ID) = F(L0TMOL+ID)
ELSEIF(STEADY.AND.ISWEEP.GE.FSWEEP.AND.ISWEEP.GE.IFIRST.
1 AND.MOD(ISWEEP,IFREQU).EQ.0) THEN
F(L0MOLO+ID) = F(L0TMOL+ID)
ENDIF
F(L0TMOL+ID) = F(L0MOLO+ID) + RLXFAC*F(L0INOU+ID)
1962 CONTINUE
ENDIF
IF(DTFALS(LNVFOL).LT.0.0) THEN
RELAX=-DTFALS(LNVFOL)
CALL FN10(LNVFOL,LNVFOL,-L0EASP,0.0,RELAX,1.0-RELAX)
ENDIF
ELSEIF(IGR.EQ.19.AND.ISC.EQ.7) THEN
C>>>>>>>>>>>>>>>>>>>>> Group 19 ------- Section 7 finish of sweep
C Update tmol for height-in-z-direction = .true.
C fn10(y,x1,x2,a,b1,b2) : y = a + b1*x1 + b2*x2
C.......................................................................
ICHK=0
IF(ZDIREC) THEN
CALL FN10(-L0TMOL,-L0MOLO,-L0INOU,0.0,1.0,RLXFAC)
IF(STORE(LBNAME('WAVE')))
1 CALL FN15(ANYZ(LBNAME('WAVE'),1),-L0TMOL,AHIGH,0.0,1.0/RHOLIQ)
ENDIF
ENDIF
NAMSUB='gxhol'
if(flag) call banner(2,namsub,0)
END
c