cGxrstm.for
c
c
c
C.... UTURSO is called from Gr.8 Sec.12 of GXRSTM
C
SUBROUTINE UTURSO
INCLUDE 'farray'
INCLUDE 'satear'
INCLUDE 'grdear'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'rsmcmn'
COMMON/INDAUX/L0ISL,L0IST,L0SL,L0ST,NSTO,NSOL,L0SLRS,L0TTRS,
1 IFL(12)
COMMON/NAMFN/NAMFUN,NAMSUB /GENI/IGFIL1(49),ITEM1,IGFIL2(10)
COMMON/RSTMCM/L0GTRS,L0PRTR
LOGICAL GTURB
CHARACTER*6 NAMFUN,NAMSUB,NM13*3
DATA G2D3/0.66666667/
C
IF(JTURML.LT.1) RETURN
NAMSUB= 'UTURSO'
J0AP= L0F(LAP)
J0SU= L0F(LSU)
IF(GTURB(INDVAR)) THEN
NM13 = NAME(INDVAR)(1:3)
J0CEDK= J0NXY
J0PHMS= J0NXY2
DO 10 I= 1,NXNY
F(J0CEDK+I)= F(L0GCRS+ ISL(INDVAR))*F(J0EPDK+I)
F(J0PHMS+I)= F(J0DEN1+I)*F(J0VOL+I)
F(J0AP +I)= F(J0AP+I) + F(J0PHMS+I)*F(J0CEDK+I)
10 CONTINUE
C----------------------------------------------------------- KE
IF(INDVAR.EQ.KE) THEN
DO 20 I= 1,NXNY
F(J0SU+I)= F(J0SU+I) + F(J0PHMS+I)*(F(J0PK+I)
1 - F(J0EPDK+I)*F(J0KE+I))
20 CONTINUE
C----------------------------------------------------------- EP
ELSEIF(INDVAR.EQ.EP) THEN
IF(IRSMHM.EQ.3) THEN
GCE1D1= 2.*GCE1/GC1
GCE2D1= 2.*GCE2/GC1
ELSE
GCE1D1= GCE1/GC1
GCE2D1= GCE2/GC1
ENDIF
DO 30 I= 1,NXNY
F(J0SU+I) = F(J0SU+I) + F(J0PHMS+I) *
1 (F(J0CEDK+I)*(GCE1D1*F(J0PK+I)-GCE2D1*F(J0EPDK+I)*F(J0KE+I)))
30 CONTINUE
C----------------------------------------------------------- U2RS
ELSEIF(INDVAR.EQ.JU2RS) THEN
DO 40 I= 1,NXNY
F(J0SU+I)= F(J0SU+I) + F(J0PHMS+I) *
1 ( F(J0CEDK+I) * F(J0KE+I)*(G2D3 - F(J0U2DK+I) )
1 + F(J0PU2 +I) + F(J0R2U2+I) - G2D3*F(J0EPDK+I)*F(J0KE+I) )
40 CONTINUE
IF(IRSMHM.EQ.3) THEN
DO 41 I= 1,NXNY
F(J0SU+I)= F(J0SU+I) + F(J0PHMS+I) * ( F(J0R3U2+I) +
1 F(J0R4U2+I) + F(J0R5U2+I) - GC1ST*F(J0PK+I)*F(J0BU2+I) )
41 CONTINUE
ENDIF
C----------------------------------------------------------- V2RS
ELSEIF(INDVAR.EQ.JV2RS) THEN
DO 50 I=1,NXNY
F(J0SU+I) = F(J0SU+I) + F(J0PHMS+I) *
1 ( F(J0CEDK+I) * F(J0KE+I)*( G2D3 - F(J0V2DK+I) )
1 + F(J0PV2 +I) + F(J0R2V2+I) - G2D3*F(J0EPDK+I)*F(J0KE+I) )
50 CONTINUE
IF(IRSMHM.EQ.3) THEN
DO 51 I= 1,NXNY
F(J0SU+I)= F(J0SU+I) + F(J0PHMS+I) * ( F(J0R3V2+I) +
1 F(J0R4V2+I) + F(J0R5V2+I) - GC1ST*F(J0PK+I)*F(J0BV2+I) )
51 CONTINUE
ENDIF
C----------------------------------------------------------- W2RS
ELSEIF(INDVAR.EQ.JW2RS) THEN
DO 60 I= 1,NXNY
F(J0SU+I)= F(J0SU+I) + F(J0PHMS+I) *
1 ( F(J0CEDK+I) * F(J0KE+I)*( G2D3 - F(J0W2DK+I) )
1 + F(J0PW2 +I) + F(J0R2W2+I) - G2D3*F(J0EPDK+I)*F(J0KE+I) )
60 CONTINUE
IF(IRSMHM.EQ.3) THEN
DO 61 I= 1,NXNY
F(J0SU+I)= F(J0SU+I) + F(J0PHMS+I) * ( F(J0R3W2+I) +
1 F(J0R4W2+I) + F(J0R5W2+I) - GC1ST*F(J0PK+I)*F(J0BW2+I) )
61 CONTINUE
ENDIF
C----------------------------------------------------------- DFRS
ELSEIF(INDVAR.EQ.JDFRS) THEN
DO 70 I= 1,NXNY
F(J0SU+I)= F(J0SU+I) + F(J0PHMS+I) *
1 ( - F(J0CEDK+I) * F(J0DFRS+I)
1 + F(J0PV2 +I) - F(J0PU2+I) + F(J0R2V2+I) - F(J0R2U2+I) )
70 CONTINUE
C----------------------------------------------------------- UVRS
ELSEIF(INDVAR.EQ.JUVRS) THEN
DO 80 I= 1,NXNY
F(J0SU+I)= F(J0SU+I) + F(J0PHMS+I) *
1 (-F(J0CEDK+I)*F(J0UVDK+I)*F(J0KE+I)+F(J0PUV+I)+F(J0R2UV+I))
80 CONTINUE
IF(IRSMHM.EQ.3) THEN
DO 81 I= 1,NXNY
F(J0SU+I) = F(J0SU+I) + F(J0PHMS+I) * ( F(J0R3UV+I) +
1 F(J0R4UV+I) + F(J0R5UV+I) - GC1ST*F(J0PK+I)*F(J0BUV+I) )
81 CONTINUE
ENDIF
C----------------------------------------------------------- UWRS
ELSEIF(INDVAR.EQ.JUWRS) THEN
DO 90 I= 1,NXNY
F(J0SU+I)= F(J0SU+I) + F(J0PHMS+I) *
1 (-F(J0CEDK+I)*F(J0UWDK+I)*F(J0KE+I)+F(J0PUW+I)+F(J0R2UW+I))
90 CONTINUE
IF(IRSMHM.EQ.3) THEN
DO 91 I= 1,NXNY
F(J0SU+I) = F(J0SU+I) + F(J0PHMS+I) * ( F(J0R3UW+I) +
1 F(J0R4UW+I) + F(J0R5UW+I) - GC1ST*F(J0PK+I)*F(J0BUW+I) )
91 CONTINUE
ENDIF
C----------------------------------------------------------- VWRS
ELSEIF(INDVAR.EQ.JVWRS) THEN
DO 100 I= 1,NXNY
F(J0SU+I)= F(J0SU+I) + F(J0PHMS+I) *
1 (-F(J0CEDK+I)*F(J0VWDK+I)*F(J0KE+I)+F(J0PVW+I)+F(J0R2VW+I))
100 CONTINUE
IF(IRSMHM.EQ.3) THEN
DO 101 I= 1,NXNY
F(J0SU+I) = F(J0SU+I) + F(J0PHMS+I) * ( F(J0R3VW+I) +
1 F(J0R4VW+I) + F(J0R5VW+I) - GC1ST*F(J0PK+I)*F(J0BVW+I) )
101 CONTINUE
ENDIF
C----------------------------------------------------------- USCR
ELSEIF(NM13.EQ.'USC' .OR. NM13.EQ.'UTR') THEN
J0USCR= L0F(INDVAR)
DO 110 I= 1,NXNY
F(J0SU+I)= F(J0SU+I) + F(J0PHMS+I) *
1 (-F(J0CEDK+I)*F(J0USCR+I)
1 +F(J0PUS1+I)+(1.-GC2T)*F(J0PUS2+I))
110 CONTINUE
C----------------------------------------------------------- VSCR
ELSEIF(NM13.EQ.'VSC' .OR. NM13.EQ.'VTR') THEN
J0VSCR= L0F(INDVAR)
DO 120 I= 1,NXNY
F(J0SU+I) = F(J0SU+I) + F(J0PHMS+I) *
1 (-F(J0CEDK+I)*F(J0VSCR+I)
1 +F(J0PVS1+I)+(1.-GC2T)*F(J0PVS2+I))
120 CONTINUE
C----------------------------------------------------------- WSCR
ELSEIF(NM13.EQ.'WSC' .OR. NM13.EQ.'WTR') THEN
J0WSCR= L0F(INDVAR)
DO 130 I= 1,NXNY
F(J0SU+I) = F(J0SU+I) + F(J0PHMS+I) *
1 (-F(J0CEDK+I)*F(J0WSCR+I)
1 +F(J0PWS1+I)+(1.-GC2T)*F(J0PWS2+I))
130 CONTINUE
ENDIF
ELSEIF(JTURML.GT.1) THEN
IF(INDVAR.EQ.W1) THEN
C----------------------------------------------- W1
IF(.NOT.(PARAB.AND.IZ.EQ.1 .OR. NZ.EQ.1)) THEN
JPZ= ITWO(0,JUMPZ,PARAB)
DO 140 I= 1,NXNY
GAH= F(J0AH+I)
GAL= GAH
F(J0SU+I) = F(J0SU+I) - F(J0DEN1+I) *
1 (GAH*F(J0W2RS+I+JPZ) - GAL*F(J0W2RS+I-JUMPZ+JPZ))
140 CONTINUE
ENDIF
IF(NX.NE.1) THEN
NXL= ITWO(NX,NX-1,XCYCLE)
JPZ= ITWO(0,JUMPZ,PARAB.OR.NZ.EQ.1)
DO 150 IX= 1,NXL
I= (IX-1)*NY
IPLUS= MOD(I+NY,NXNY)
DO 150 IY= 1,NY
I= I + 1
IPLUS= IPLUS + 1
GAE = F(J0AE+I)
GUWE= 0.25*(F(J0UWRS+I) +F(J0UWRS+IPLUS)
1 +F(J0UWRS+I+JPZ)+F(J0UWRS+IPLUS+JPZ))
SUSUB= F(J0DEN1+I)*GUWE*GAE
F(J0SU+I) = F(J0SU+I) - SUSUB
F(J0SU+IPLUS)= F(J0SU+IPLUS) + SUSUB
150 CONTINUE
ENDIF
IF(NY.NE.1) THEN
JPZ= ITWO(0,JUMPZ,PARAB.OR.NZ.EQ.1)
DO 160 IX= 1,NX
I= (IX-1)*NY
DO 160 IY= 1,NY-1
I= I + 1
GAN = F(J0AN+I)
GVWN= 0.25*(F(J0VWRS+I) +F(J0VWRS+I+1)
1 +F(J0VWRS+I+JPZ)+F(J0VWRS+I+1+JPZ))
SUSUB= F(J0DEN1+I)*GVWN*GAN
F(J0SU+I) = F(J0SU+I) - SUSUB
F(J0SU+I+1)= F(J0SU+I+1) + SUSUB
160 CONTINUE
ENDIF
C-------------------------------------- Sources for V1
ELSEIF(INDVAR.EQ.V1) THEN
IF(NY.NE.1) THEN
DO 170 IX= 1,NX
I = (IX-1)*NY
GAN = F(J0VOL+I+1)/F(J0DY+I+1)
GV2N= F(J0V2RS+I+1)
DO 170 IY= 1,NY-1
I = I + 1
GAS = GAN
GV2S = GV2N
GAN = F(J0VOL+I+1)/F(J0DY+I+1)
GV2N = F(J0V2RS+I+1)
GSUAD= - F(J0DEN1+I)*(GAN*GV2N - GAS*GV2S)
IF(.NOT.CARTES) THEN
GVOL = 0.5*(F(J0VOL+I)+F(J0VOL+I+1))
GR = 0.5*(F(J0R+I)+F(J0R+I+JUMPY))
GU2RS= 0.5*(F(J0U2RS+I)+F(J0U2RS+I+1))
GSUAD= GSUAD+F(J0DEN1+I)*GVOL/GR*GU2RS
ENDIF
F(J0SU+I)= F(J0SU+I) + GSUAD
170 CONTINUE
ENDIF
IF(NZ.NE.1 .AND. .NOT.(PARAB.AND.IZ.EQ.1)) THEN
JPZ = ITWO(0,JUMPZ,PARAB)
GADH= ITWO(1,0,PARAB.OR.IZ.NE.NZ)
GADL= ITWO(1,0,PARAB.OR.IZ.NE.1)
DO 180 IX= 1,NX
I= (IX-1)*NY
DO 180 IY= 1,NY-1
I = I+1
GAH = 0.5*(F(J0AH+I)+F(J0AH+I+1))
GAL = GAH
GVWH= 0.25*(F(J0VWRS+I)+F(J0VWRS+I+1)
1 +F(J0VWRS+I+JPZ)+F(J0VWRS+I+1+JPZ))
GVWL= 0.25*(F(J0VWRS+I-JUMPZ)+F(J0VWRS+I+1-JUMPZ)
1 +F(J0VWRS+I-JUMPZ+JPZ)+F(J0VWRS+I+1-JUMPZ+JPZ))
F(J0SU+I)= F(J0SU+I)-F(J0DEN1+I)
1 *(GADH*GAH*GVWH-GADL*GAL*GVWL)
180 CONTINUE
ENDIF
IF(NX.NE.1) THEN
NXL= ITWO(NX,NX-1,XCYCLE)
DO 190 IX= 1,NXL
I= (IX-1)*NY
IPLUS= MOD(I+NY,NXNY)
DO 190 IY= 1,NY-1
I = I + 1
IPLUS= IPLUS + 1
GAE = 0.5*(F(J0AE+I)+F(J0AE+I+1))
GUVE= 0.25*(F(J0UVRS+I) +F(J0UVRS+IPLUS)
1 +F(J0UVRS+I+1)+F(J0UVRS+IPLUS+1))
SUSUB= F(J0DEN1+I)*GAE*GUVE
F(J0SU+I) = F(J0SU+I) - SUSUB
F(J0SU+IPLUS)= F(J0SU+IPLUS)+ SUSUB
190 CONTINUE
ENDIF
C-------------------------------------- Sources for U1
ELSEIF(INDVAR.EQ.U1) THEN
IF(NX.NE.1) THEN
NXL= ITWO(NX,NX-1,XCYCLE)
DO 200 IX= 1,NXL
I= (IX-1)*NY
IPLUS= MOD(I+NY,NXNY)
DO 200 IY= 1,NY
I= I + 1
IPLUS= IPLUS + 1
GAE = F(J0VOL+IPLUS)/F(J0DX+IPLUS)
GAW = F(J0VOL+I)/F(J0DX+I)
GU2E= F(J0U2RS+IPLUS)
GU2W= F(J0U2RS+I)
F(J0SU+I)= F(J0SU+I)-F(J0DEN1+I)*(GAE*GU2E-GAW*GU2W)
200 CONTINUE
ENDIF
IF(NZ.NE.1 .AND. .NOT.(PARAB.AND.IZ.EQ.1)) THEN
NXL = ITWO(NX,NX-1,XCYCLE.OR.NX.EQ.1)
GADH= ITWO(1,0,PARAB.OR.IZ.NE.NZ)
GADL= ITWO(1,0,PARAB.OR.IZ.NE.1)
JPZ = ITWO(0,JUMPZ,PARAB)
DO 210 IX= 1,NXL
I= (IX-1)*NY
IPLUS= MOD(I+NY,NXNY)
DO 210 IY= 1,NY
I = I + 1
IPLUS= IPLUS + 1
GAH = 0.5*(F(J0AH+I)+F(J0AH+IPLUS))
GAL = GAH
GUWH = 0.25*(F(J0UWRS+I)+F(J0UWRS+IPLUS)
1 + F(J0UWRS+I+JPZ)+F(J0UWRS+IPLUS+JPZ))
GUWL = 0.25*(F(J0UWRS+I-JUMPZ)+F(J0UWRS+IPLUS-JUMPZ)
1 + F(J0UWRS+I+JPZ-JUMPZ)+F(J0UWRS+IPLUS+JPZ-JUMPZ))
F(J0SU+I)= F(J0SU+I)-F(J0DEN1+I)
1 *(GADH*GAH*GUWH-GADL*GAL*GUWL)
210 CONTINUE
ENDIF
IF(NY.NE.1) THEN
NXL= ITWO(NX,NX-1,XCYCLE.OR.NX.EQ.1)
DO 220 IX=1,NXL
I= (IX-1)*NY
IPLUS= MOD(I+NY,NXNY)
DO 220 IY= 1,NY-1
I= I + 1
IPLUS= IPLUS + 1
GAN = 0.5*(F(J0AN+I)+F(J0AN+IPLUS))
GUVN = 0.25*(F(J0UVRS+I) + F(J0UVRS+I+1) +
1 F(J0UVRS+IPLUS) + F(J0UVRS+IPLUS+1))
SUSUB= F(J0DEN1+I)*GAN*GUVN
F(J0SU+I) = F(J0SU+I) - SUSUB
F(J0SU+I+1)= F(J0SU+I+1)+ SUSUB
220 CONTINUE
ENDIF
IF(.NOT.CARTES) THEN
NXL= ITWO(NX,NX-1,XCYCLE.OR.NX.EQ.1)
DO 230 IX= 1,NXL
I= (IX-1)*NY
IPLUS= MOD(I+NY,NXNY)
DO 230 IY= 1,NY
I = I + 1
IPLUS= IPLUS + 1
GVOL = 0.5*(F(J0VOL+I)+F(J0VOL+IPLUS))
GR = 0.5*(F(J0R+I)+F(J0R+IPLUS))
GUVRS= 0.5*(F(J0UVRS+I)+F(J0UVRS+IPLUS))
F(J0SU+I)= F(J0SU+I) - F(J0DEN1+I)*GVOL/GR*GUVRS
230 CONTINUE
ENDIF
C---------------------------------- Sources for H1,TEM1 & scalars
ELSEIF( (INDVAR.EQ.H1.OR.INDVAR.EQ.ITEM1.OR.NAME(INDVAR)(1:2).
1 EQ.'SC') .AND. JSCAML.GE.2) THEN
C.... Allow for specific heat for TEM1, but for now just use cp,p
IF(INDVAR.EQ.ITEM1) THEN
CALL FN0(-J0NXY,LCP1)
ELSE
CALL FN1(-J0NXY,1.0)
ENDIF
IF(NX.NE.1) THEN
J0USCR= J0SCRS(INDVAR,J0UTRS,'U')
NXL = ITWO(NX,NX-1,XCYCLE)
DO 240 IX= 1,NXL
I = (IX-1)*NY
IPLUS= MOD(I+NY,NXNY)
DO 240 IY=1,NY
I = I + 1
IPLUS= IPLUS + 1
GAE = F(J0AE+I)
GUSE = 0.5*(F(J0USCR+I)+F(J0USCR+IPLUS))
GCP = F(J0NXY+I)
SUSUB= F(J0DEN1+I)*GUSE*GAE*GCP
F(J0SU+I) = F(J0SU+I) - SUSUB
F(J0SU+IPLUS)= F(J0SU+IPLUS) + SUSUB
240 CONTINUE
ENDIF
IF(NY.NE.1) THEN
J0VSCR= J0SCRS(INDVAR,J0VTRS,'V')
DO 250 IX= 1,NX
I= (IX-1)*NY
DO 250 IY=1,NY-1
I = I + 1
GAN = F(J0AN+I)
GVSN= 0.5*(F(J0VSCR+I)+F(J0VSCR+I+1))
GCP = F(J0NXY+I)
SUSUB= F(J0DEN1+I)*GVSN*GAN*GCP
F(J0SU+I) = F(J0SU+I) - SUSUB
F(J0SU+I+1)= F(J0SU+I+1)+ SUSUB
250 CONTINUE
ENDIF
IF(NZ.NE.1 .AND. .NOT.(PARAB.AND.IZ.EQ.1)) THEN
J0WSCR= J0SCRS(INDVAR,J0WTRS,'W')
GADH = ITWO(1,0,PARAB.OR.IZ.NE.NZ)
GADL = ITWO(1,0,PARAB.OR.IZ.NE.1)
JPZ = ITWO(0,JUMPZ,PARAB)
DO 260 I= 1,NXNY
GAH = F(J0AH+I)
GAL = GAH
GWSH= 0.5*(F(J0WSCR+I)+F(J0WSCR+I+JPZ))
GWSL= 0.5*(F(J0WSCR+I-JUMPZ)+F(J0WSCR+I+JPZ-JUMPZ))
GCP = F(J0NXY+I)
F(J0SU+I)= F(J0SU+I)-F(J0DEN1+I)*(GADH*GAH*GWSH-
1 GADL*GAL*GWSL)*GCP
260 CONTINUE
ENDIF
ENDIF
ENDIF
C.... Common part:
IF(SOLVE(INDVAR)) THEN
CALL SUMAIF(SUMABS,J0SU,J0AP,FIXVAL)
SUMABS= SUMABS/(RESREF(INDVAR)+TINY)
F(L0SLRS+ISL(INDVAR))= SUMABS
IF(.NOT.(PARAB .OR. ITHYD.LT.LITHYD)) THEN
IF(IZSTEP.LE.1) F(L0GTRS+ISL(INDVAR)) = 0.0
F(L0GTRS+ISL(INDVAR)) = F(L0GTRS+ISL(INDVAR)) + SUMABS
F(L0TTRS+ISL(INDVAR)) = F(L0GTRS+ISL(INDVAR))
ENDIF
ENDIF
NAMSUB = 'uturso'
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C.... URSVIS is called from Gr.9, Sec.5 of GXRSTM.
C
SUBROUTINE URSVIS
INCLUDE 'farray'
INCLUDE 'satear'
INCLUDE 'grdear'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'rsmcmn'
COMMON /NAMFN/NAMFUN,NAMSUB
CHARACTER*6 NAMFUN,NAMSUB
LOGICAL QEQ,SLD
C
NAMSUB= 'URSVIS'
IF(QEQ(ENUT,GRND))THEN
J0VT= L0F(AUX(VIST))
IF(JTURML.GT.0) THEN
DO 10 I= 1,NXNY
IF(SLD(I)) THEN
F(J0VT+I)= 0.0
ELSE
F(J0EPDK+I)=F(J0EP+I)/(F(J0KE+I)+TINY)
F(J0VT+I)= GCS*F(J0KE+I)/F(J0EPDK+I)
ENDIF
10 CONTINUE
ELSE
DO 20 I= 1,NXNY
20 F(J0VT+I)= 0.0
ENDIF
ENDIF
NAMSUB= 'ursvis'
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C
SUBROUTINE URSLEN
INCLUDE 'farray'
INCLUDE 'satear'
INCLUDE 'grdear'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'rsmcmn'
COMMON /NAMFN/NAMFUN,NAMSUB
CHARACTER*6 NAMFUN,NAMSUB
LOGICAL QEQ,SLD
C
NAMSUB= 'URSLEN'
IF(QEQ(EL1,GRND)) THEN
J0L1= L0F(AUX(LEN1))
IF(JTURML.GT.0) THEN
ELMN= VARMIN(LEN1)
ELMX= VARMAX(LEN1)
DO 10 I= 1,NXNY
F(J0KE+I)= AMAX1(F(J0KE+I),TINY)
F(J0EP+I)= AMAX1(F(J0EP+I),TINY)
IF(SLD(I)) THEN
F(J0L1+I)= 0.0
ELSE
F(J0L1+I)= AMAX1(ELMN,AMIN1(ELMX,
1 CD*F(J0KE+I)**1.5/F(J0EP+I)))
ENDIF
10 CONTINUE
ELSE
DO 20 I= 1,NXNY
20 F(J0L1+I)=0.0
ENDIF
ENDIF
NAMSUB= 'urslen'
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C.... UDFMOD is called from Gr.8, Sec. 9 of GXRSTM.
C
SUBROUTINE UDFMOD
INCLUDE 'farray'
INCLUDE 'satear'
INCLUDE 'grdear'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'rsmcmn'
COMMON /F0/KF01(70),KAP,KSU,KF073(232)
1 /GENI/IGFIL1(49),ITEM1,IGFIL2(10) /NAMFN/NAMFUN,NAMSUB
LOGICAL GU2OV2,GTURB
CHARACTER*6 NAMFUN,NAMSUB
C
NAMSUB= 'UDFMOD'
IF(JTURML.GT.1) THEN
IF(NDIREC.EQ.1) THEN
IF(GTURB(INDVAR)) THEN
J0DIFN= L0F(LAN)
J0DIFS= L0F(LAS)
DO 10 IX= 1,NX
I= (IX-1)*NY
DO 10 IY= 1,NY-1
I= I+1
GMOD = 0.5*(F(J0V2DK+I) + F(J0V2DK+I+JUMPY))
F(J0DIFN+I)= F(J0DIFN+I)*GMOD*F(L0GCRT+ ISL(INDVAR))
10 CONTINUE
ENDIF
ELSEIF(NDIREC.EQ.3) THEN
IF(GTURB(INDVAR)) THEN
J0DIFE= L0F(LAE)
NXL= ITWO(NX,NX-1,XCYCLE)
I = 0
DO 20 IX= 1,NXL
IPLUS= MOD(I+NY,NXNY)
DO 20 IY= 1,NY
I = I + 1
IPLUS= IPLUS + 1
GMOD = 0.5*(F(J0U2DK+I)+F(J0U2DK+IPLUS))
F(J0DIFE+I)= F(J0DIFE+I)*GMOD*F(L0GCRT+ ISL(INDVAR))
20 CONTINUE
ENDIF
ELSEIF(NDIREC.EQ.5) THEN
IF(GTURB(INDVAR)) THEN
J0DIFH= L0F(LD11)
DO 30 I= 1,NXNY
GMOD= F(J0W2DK+I)
F(J0DIFH+I)= F(J0DIFH+I)*GMOD*F(L0GCRT+ ISL(INDVAR))
30 CONTINUE
ENDIF
ENDIF
ENDIF
IF(.NOT.CARTES .AND. NDIREC.EQ.1) THEN
IF(INDVAR.EQ.V1 .OR. GTURB(INDVAR)) THEN
J0VAR = L0F(INDVAR)
J0DIFN= L0F(LAN)
IF(INDVAR.EQ.JWTRS .OR. NAME(INDVAR)(1:3).EQ.'WSC') THEN
DO 40 I= 1,NXNY
40 F(J0NXY+I)= 0.0
ELSEIF(INDVAR.EQ.V1) THEN
DO 50 I= 1,NXNY
50 F(J0NXY+I)= 1.0
ELSEIF(INDVAR.EQ.JUVRS) THEN
DO 60 I= 1,NXNY
60 F(J0NXY+I)= 4.*F(J0U2DK+I)/F(J0V2DK+I)
ELSEIF(INDVAR.EQ.JUWRS.OR.INDVAR.EQ.JVWRS.OR.NAME(INDVAR)(2:3)
1 .EQ.'SC'.OR.NAME(INDVAR)(2:4).EQ.'TRS') THEN
DO 70 I= 1,NXNY
70 F(J0NXY+I)= F(J0U2DK+I)/F(J0V2DK+I)*F(L0GCRT+ ISL(INDVAR))
ELSEIF(INDVAR.EQ.JU2RS) THEN
DO 80 I= 1,NXNY
GU2DV2= F(J0U2DK+I)/F(J0V2DK+I)
F(J0NXY+I)= -2.*(1.-GU2DV2)
80 CONTINUE
ELSEIF(INDVAR.EQ.JV2RS) THEN
DO 90 I= 1,NXNY
GU2DV2= F(J0U2DK+I)/F(J0V2DK+I)
F(J0NXY+I)= 2.*GU2DV2*(1.-GU2DV2)
90 CONTINUE
ELSE
RETURN
ENDIF
JJ0R= ITWO(J0RV,J0R,INDVAR.EQ.V1)
JJNY= ITWO(NY-1,NY,INDVAR.EQ.V1)
GRADP = TINY
GRADN = F(JJ0R+1)
GU2OV2= INDVAR.EQ.JU2RS .OR. INDVAR.EQ.JV2RS
DO 110 IY = 1,JJNY
GRADS= GRADP
GRADP= GRADN
GRADN= F(JJ0R+IY+1)
APDAN= GRADN/GRADP - 1.
APDAS= GRADS/GRADP - 1.
IF(IY.EQ.1) APDAS= 0.0
I= IY - NY
DO 100 IX = 1,NX
I= I + NY
APAD= F(J0DIFN+I)*APDAN + F(J0DIFN+I-1)*APDAS
C.... Provide protection against -ve ap for u2rs & v2rs by not
C multiplying ap by f(j0nxy+i); since solving a correction
C equation this is permissible, although in future it may
C be better to adopt a more consistent treatment.
SUAD= - F(J0VAR+I)*APAD*F(J0NXY+I)
IF(.NOT.GU2OV2) APAD= APAD*F(J0NXY+I)
F(KAP+I)= F(KAP+I) + APAD
F(KSU+I)= F(KSU+I) + SUAD
100 CONTINUE
110 CONTINUE
ENDIF
ENDIF
C.... Eddy-diffusivity option for H1; cmucd=taudke*taudke
C for consistency with RSTM constants
IF((INDVAR.EQ.H1.OR.INDVAR.EQ.ITEM1.OR.NAME(INDVAR)(1:2).EQ.'SC')
1 .AND.JSCAML.EQ.1) THEN
IF(NDIREC.EQ.1) THEN
CALL FN25(LAN,F(L0GCRT+ ISL(INDVAR)))
CALL FN25(LAS,F(L0GCRT+ ISL(INDVAR)))
ELSEIF(NDIREC.EQ.3) THEN
CALL FN25(LAE,F(L0GCRT+ ISL(INDVAR)))
CALL FN25(LAW,F(L0GCRT+ ISL(INDVAR)))
ELSEIF(NDIREC.EQ.5) THEN
CALL FN25(LD11,F(L0GCRT+ ISL(INDVAR)))
ENDIF
ENDIF
NAMSUB = 'udfmod'
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C.... UCNMOD is called from Gr.9, Sec.10 of UCNMOD.
C
SUBROUTINE UCNMOD
INCLUDE 'farray'
INCLUDE 'satear'
INCLUDE 'grdear'
INCLUDE 'grdloc'
INCLUDE 'satgrd'
INCLUDE 'rsmcmn'
COMMON /F0/KF01(71),KSU,KF073(232) /NAMFN/NAMFUN,NAMSUB
LOGICAL GTURB
CHARACTER*6 NAMFUN,NAMSUB
CHARACTER*4 NMPH
C
NAMSUB= 'UCNMOD'
IF(GTURB(INDVAR).AND..NOT.CARTES.AND.STORE(U1)) THEN
IF(NX.EQ.1.AND.NDIREC.EQ.1.OR.NX.NE.1.AND.NDIREC.EQ.3) THEN
IF(INDVAR.EQ.JU2RS) THEN
DO 10 I= 1,NXNY
10 F(J0NXY+I)= -2.*F(J0UVRS+I)
ELSEIF(INDVAR.EQ.JV2RS) THEN
DO 20 I= 1,NXNY
20 F(J0NXY+I)= 2.*F(J0UVRS+I)
ELSEIF(INDVAR.EQ.JW2RS) THEN
RETURN
ELSEIF(INDVAR.EQ.JDFRS) THEN
DO 30 I= 1,NXNY
30 F(J0NXY+I)= 4.*F(J0UVRS+I)
ELSEIF(INDVAR.EQ.KE) THEN
RETURN
ELSEIF(INDVAR.EQ.EP) THEN
RETURN
ELSEIF(INDVAR.EQ.JUVRS) THEN
DO 40 I= 1,NXNY
40 F(J0NXY+I)= F(J0U2RS+I)-F(J0V2RS+I)
ELSEIF(INDVAR.EQ.JUWRS) THEN
DO 50 I= 1,NXNY
50 F(J0NXY+I)= -F(J0VWRS+I)
ELSEIF(INDVAR.EQ.JVWRS) THEN
DO 60 I= 1,NXNY
60 F(J0NXY+I)= F(J0UWRS+I)
ELSEIF(INDVAR.EQ.JUTRS) THEN
DO 70 I= 1,NXNY
70 F(J0NXY+I)= -F(J0VTRS+I)
ELSEIF(INDVAR.EQ.JVTRS) THEN
DO 80 I= 1,NXNY
80 F(J0NXY+I)= F(J0UTRS+I)
ELSEIF(INDVAR.EQ.JWTRS.OR.NAME(INDVAR)(1:3).EQ.'WSC') THEN
RETURN
ELSEIF(NAME(INDVAR)(1:3).EQ.'USC') THEN
NMPH = 'VSC'//NAME(INDVAR)(4:4)
J0VSCR= L0F(LBNAME(NMPH))
DO 90 I= 1,NXNY
90 F(J0NXY+I)= -F(J0VSCR+I)
ELSEIF(NAME(INDVAR)(1:3).EQ.'VSC') THEN
NMPH = 'USC'//NAME(INDVAR)(4:4)
J0USCR= L0F(LBNAME(NMPH))
DO 100 I= 1,NXNY
100 F(J0NXY+I)= F(J0USCR+I)
ENDIF
ENDIF
IF(NX.EQ.1.AND.NDIREC.EQ.1) THEN
DO 110 I= 1,NXNY
F(KSU+I)= F(KSU+I)+F(J0DEN1+I)*F(J0VOL+I)*F(J0NXY+I)*
1 F(J0U1+I)/F(J0R+I)
110 CONTINUE
ELSEIF(NDIREC.EQ.3) THEN
J0CNE= L0F(LD8)
J0CNW= L0F(LD7)
NXL =ITWO(NX,NX-1,XCYCLE)
DO 120 IX= 1,NXL
I= (IX-1)*NY
IPLUS= (MOD(IX+NX,NX))*NY
DO 120 IY= 1,NY
I= I + 1
IPLUS= IPLUS + 1
GFIPE= 0.5*(F(J0NXY+I)+F(J0NXY+IPLUS))
GDX = F(J0DX+I)
GRP = F(J0R+I)
TERM = GFIPE*GDX/GRP
F(J0CNE+I) = F(J0CNE+I) - TERM
F(J0CNW+IPLUS)= F(J0CNW+IPLUS)+ TERM
120 CONTINUE
ENDIF
ENDIF
NAMSUB= 'ucnmod'
END
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c