C ********************************************************************** C * * C * SOFTWARE LICENSING * C * * C * This program is free software; you can redistribute * C * it and/or modify it under the terms of the GNU * C * General Public License as published by the Free * C * Software Foundation, either Version 2 of the * C * license, or (at your option) any later version. * C * * C * This program is distributed in the hope that it * C * will be useful, but without any warranty; without * C * even the implied warranty of merchantability or * C * fitness for a particular purpose. See the GNU * C * General Public License for more details. * C * * C * A copy of the GNU General Public License is * C * available at http://www.gnu.org/copyleft/gpl.html * C * or by writing to the Free Software Foundation, Inc.,* C * 59 Temple Place - Suite 330, Boston, MA 02111, USA. * C * * C ********************************************************************** SUBROUTINE ADVT(FB,F,FMEAN,DTI2,FF,FDIF2,TSDIS,TSBDRY) C INCLUDE 'comdeck' SAVE REAL TSQEPC(EPCM,KBM1) C C----------------------------------------------------------------------- C THIS SUBROUTINE INTEGRATES CONSERVATIVE CONSTITUENT EQUATIONS C----------------------------------------------------------------------- C DIMENSION FB(IM,JM,KB),F(IM,JM,KB),FF(IM,JM,KB) DIMENSION FMEAN(IM,JM,KB),FBSTART(IM,JM,KB),ETA(IM,JM) DIMENSION XMFLUX(IM,JM,KB),YMFLUX(IM,JM,KB),ZZFLUX(IM,JM,KB) DIMENSION XFLUX(IM,JM,KB),YFLUX(IM,JM,KB),ZFLUX(IM,JM,KB) DIMENSION FDIF2(DBCM,KBM1),TSDIS(QBCM),TSBDRY(EBCM,KBM1) EQUIVALENCE (XFLUX,A),(YFLUX,C),(ZFLUX,PROD) EQUIVALENCE (FBSTART,VH),(ETA,VHP) C IF(SCHEME .EQ. 'CENTRAL ')THEN DO 9 J=1,JM DO 9 I=1,IM F (I,J,KB)=F (I,J,KBM1) 9 FB(I,J,KB)=FB(I,J,KBM1) C C-------- DO ADVECTION FLUXES AS CENTER DIFFERENCES -------------------- DO 22 K=1,KBM1 DO 22 J=2,JM DO 22 I=2,IM XFLUX(I,J,K)=0.5*(F(I-1,J,K)+F(I,J,K)) * XMFL3D(I,J,K) 22 YFLUX(I,J,K)=0.5*(F(I,J-1,K)+F(I,J,K)) * YMFL3D(I,J,K) C C DO 236 N=1,NUMEBC IE=IETA(N) JE=JETA(N) IC=ICON(N) JC=JCON(N) C IF(FSM(IE+1,JE).EQ.0.0.AND.JE.EQ.JC) THEN C C-------- EAST SIDE ---------------------------------------------------- DO 21 K=1,KBM1 UUP =0.5*(XMFL3D(IE,JE,K)+ABS(XMFL3D(IE,JE,K))) UDOWN=0.5*(XMFL3D(IE,JE,K)-ABS(XMFL3D(IE,JE,K))) 21 XFLUX(IE,JE,K)=UUP*FB(IC,JC,K)+UDOWN*FB(IE,JE,K) GO TO 236 C ELSE IF (FSM(IE-1,JE).EQ.0.0.AND.JE.EQ.JC) THEN C C-------- WEST SIDE ---------------------------------------------------- DO 222 K=1,KBM1 UUP =0.5*(XMFL3D(IC,JC,K)+ABS(XMFL3D(IC,JC,K))) UDOWN=0.5*(XMFL3D(IC,JC,K)-ABS(XMFL3D(IC,JC,K))) 222 XFLUX(IC,JC,K)=UUP*FB(IE,JE,K)+UDOWN*FB(IC,JC,K) GO TO 236 C ELSE IF(FSM(IE,JE+1).EQ.0.0.AND.IE.EQ.IC) THEN C C-------- NORTH SIDE --------------------------------------------------- DO 23 K=1,KBM1 VUP =0.5*(YMFL3D(IE,JE,K)+ABS(YMFL3D(IE,JE,K))) VDOWN=0.5*(YMFL3D(IE,JE,K)-ABS(YMFL3D(IE,JE,K))) 23 YFLUX(IE,JE,K)=VUP*FB(IC,JC,K)+VDOWN*FB(IE,JE,K) GO TO 236 C ELSE IF(FSM(IE,JE-1).EQ.0.0.AND.IE.EQ.IC) THEN C C-------- SOUTH SIDE --------------------------------------------------- DO 24 K=1,KBM1 VUP =0.5*(YMFL3D(IC,JC,K)+ABS(YMFL3D(IC,JC,K))) VDOWN=0.5*(YMFL3D(IC,JC,K)-ABS(YMFL3D(IC,JC,K))) 24 YFLUX(IC,JC,K)=VUP*FB(IE,JE,K)+VDOWN*FB(IC,JC,K) C C-------- DONE --------------------------------------------------------- END IF 236 CONTINUE C DO 121 J=2,JMM1 DO 121 I=2,IMM1 121 FF(I,J,1)=-.5*DZR(1)*(F(I,J,1)+F(I,J,2))*W(I,J,2)*ART(I,J) DO 130 J=2,JMM1 DO 130 K=2,KBM1 DO 130 I=2,IMM1 130 FF(I,J,K)=.5*DZR(K)*((F(I,J,K-1)+F(I,J,K))*W(I,J,K) . -(F(I,J,K)+F(I,J,K+1))*W(I,J,K+1))*ART(I,J) C C************* ADJUST FOR RIVER/WALL FLUXES ********************* DO 50 N=1,NUMQBC ID=IQD(N) JD=JQD(N) IC=IQC(N) JC=JQC(N) IF(JD.EQ.JC) THEN IF(IC.GT.ID) THEN DO 60 K=1,KBM1 60 XFLUX(IC,JC,K)=F(IC,JC,K)*XMFL3D(IC,JC,K) ELSE DO 70 K=1,KBM1 70 XFLUX(ID,JD,K)=F(IC,JC,K)*XMFL3D(ID,JD,K) ENDIF ELSE IF(JC.GT.JD) THEN DO 80 K=1,KBM1 80 YFLUX(IC,JC,K)=F(IC,JC,K)*YMFL3D(IC,JC,K) ELSE DO 90 K=1,KBM1 90 YFLUX(ID,JD,K)=F(IC,JC,K)*YMFL3D(ID,JD,K) ENDIF ENDIF 50 CONTINUE C C-------- ADD NET ADVECTION FLUXES & THEN STEP FORWARD IN TIME -------- DO 141 J=2,JMM1 DO 141 K=1,KBM1 DO 141 I=2,IMM1 FF(I,J,K)=FF(I,J,K) . +XFLUX(I+1,J,K)-XFLUX(I,J,K) . +YFLUX(I,J+1,K)-YFLUX(I,J,K) 141 FF(I,J,K)=(FB(I,J,K)*(H(I,J)+ETB(I,J))*ART(I,J)-DTI2*FF(I,J,K)) . /((H(I,J)+ETF(I,J))*ART(I,J)) C ELSE C DO 10 J=1,JM DO 10 I=1,IM ETA(I,J)=ETB(I,J) 10 FB(I,J,KB)=FB(I,J,KBM1) DO 11 K=1,KB DO 11 J=1,JM DO 11 I=1,IM FBSTART(I,J,K)=FB(I,J,K) XMFLUX(I,J,K)=XMFL3D(I,J,K) YMFLUX(I,J,K)=YMFL3D(I,J,K) ZZFLUX(I,J,K)=W(I,J,K) 11 CONTINUE C----------------------------------------------------------------------- C SMOLARKIEWICZ'S SCHEME C----------------------------------------------------------------------- LOOP=0 5050 LOOP= LOOP+1 C C-------- DO ADVECTION FLUXES, FIRST AS UPWIND ------------- DO 20 K=1,KBM1 DO 20 J=2,JM DO 20 I=2,IM XFLUX(I,J,K)=AMAX1(0.,XMFLUX(I,J,K))*FBSTART(I-1,J,K)+ . AMIN1(0.,XMFLUX(I,J,K))*FBSTART(I,J,K) YFLUX(I,J,K)=AMAX1(0.,YMFLUX(I,J,K))*FBSTART(I,J-1,K)+ . AMIN1(0.,YMFLUX(I,J,K))*FBSTART(I,J,K) 20 CONTINUE DO 120 J=2,JMM1 DO 120 I=2,IMM1 ZFLUX(I,J,1)=0.0 120 ZFLUX(I,J,KB)=0.0 DO 125 J=2,JMM1 DO 125 K=2,KBM1 DO 125 I=2,IMM1 ZFLUX(I,J,K)=AMAX1(0.,ZZFLUX(I,J,K))*FBSTART(I,J,K)+ . AMIN1(0.,ZZFLUX(I,J,K))*FBSTART(I,J,K-1) 125 ZFLUX(I,J,K)=ZFLUX(I,J,K)*ART(I,J) C C----- ADD NET ADVECTIVE FLUXES & STEP FORWARD IN TIME DO 140 J=2,JMM1 DO 140 K=1,KBM1 DO 140 I=2,IMM1 FF(I,J,K)= (XFLUX(I+1,J,K)-XFLUX(I,J,K)) . +(YFLUX(I,J+1,K)-YFLUX(I,J,K)) . +DZR(K)*(ZFLUX(I,J,K) -ZFLUX(I,J,K+1)) 140 FF(I,J,K)=(FBSTART(I,J,K)*(H(I,J)+ETA(I,J))*ART(I,J)- . DTI2*FF(I,J,K))/((H(I,J)+ETF(I,J))*ART(I,J)) C C----- ANTIDIFFUSION VELOCITY ----------------------------------------- C IF(SCHEME .EQ. 'UPWIND ' .OR. LOOP .EQ. 2) GO TO 5051 C CNG09282011 CALL ANTIDIF(XMFLUX,YMFLUX,ZZFLUX,FB,FF,DTI2,TSDIS,TSBDRY) CALL ANTIDIF(XMFLUX,YMFLUX,ZZFLUX,FB,FF,DTI2,TSDIS,TSBDRY,FDIF2) C DO 146 J=1,JM DO 146 K=1,KB DO 146 I=1,IM ETA(I,J)=ETF(I,J) FBSTART(I,J,K)=FF(I,J,K) 146 CONTINUE GO TO 5050 C END IF C----------------------------------------------------------------------- C SMOLARKIEWICZ'S SCHEME END C----------------------------------------------------------------------- 5051 CONTINUE C C-------- ADD DIFFUSIVE FLUXES (HORIZONTAL) ---------------------------- C DO 30 K=1,KB ! SHALLOW SEAS C DO 30 J=1,JM C DO 30 I=1,IM C 30 FB(I,J,K)=FB(I,J,K)-FMEAN(I,J,K) DO 40 K=1,KBM1 DO 40 J=2,JM DO 40 I=2,IM CNG05192010WAD XMFLUX(I,J,K)=.5*(AAM(I,J,K)+AAM(I-1,J,K)) CNG05192010WAD YMFLUX(I,J,K)=.5*(AAM(I,J,K)+AAM(I,J-1,K)) XMFLUX(I,J,K)=.5*(AAM(I,J,K)+AAM(I-1,J,K))*DUMWAD(I,J) ! Not in original NYHOPS W&D YMFLUX(I,J,K)=.5*(AAM(I,J,K)+AAM(I,J-1,K))*DVMWAD(I,J) ! Not in original NYHOPS W&D 40 CONTINUE DO 100 J=2,JM DO 100 K=1,KBM1 DO 100 I=2,IM XFLUX(I,J,K)= . -XMFLUX(I,J,K)/HPRNU*(H(I,J)+ETB(I,J)+H(I-1,J)+ETB(I-1,J)) .*(FB(I,J,K)-FB(I-1,J,K))*DUM(I,J)*DUMWAD(I,J)/(H1(I,J)+H1(I-1,J)) ! Not in original NYHOPS W&D . *0.5*(H2(I,J)+H2(I-1,J)) YFLUX(I,J,K)= . -YMFLUX(I,J,K)/HPRNU*(H(I,J)+ETB(I,J)+H(I,J-1)+ETB(I,J-1)) .*(FB(I,J,K)-FB(I,J-1,K))*DVM(I,J)*DVMWAD(I,J)/(H2(I,J)+H2(I,J-1)) ! Not in original NYHOPS W&D . *0.5*(H1(I,J)+H1(I,J-1)) 100 CONTINUE C DO 110 K=1,KB ! SHALLOW SEAS C DO 110 J=1,JM C DO 110 I=1,IM C 110 FB(I,J,K)=FB(I,J,K)+FMEAN(I,J,K) C C------------- ADJUST FOR RIVER/WALL FLUXES --------------- DO 51 N=1,NUMQBC ID=IQD(N) JD=JQD(N) IC=IQC(N) JC=JQC(N) IF(JD.EQ.JC) THEN IF(IC.GT.ID) THEN DO 61 K=1,KBM1 61 XFLUX(IC,JC,K)=0.0 ELSE DO 71 K=1,KBM1 71 XFLUX(ID,JD,K)=0.0 ENDIF ELSE IF(JC.GT.JD) THEN DO 81 K=1,KBM1 81 YFLUX(IC,JC,K)=0.0 ELSE DO 91 K=1,KBM1 91 YFLUX(ID,JD,K)=0.0 ENDIF ENDIF 51 CONTINUE C C----- ADD NET HORIZONTAL FLUXES & STEP FORWARD IN TIME (CENTER DIFF.) DO 145 J=2,JMM1 DO 145 K=1,KBM1 DO 145 I=2,IMM1 145 FF(I,J,K)=FF(I,J,K)-DTI2*( . +(XFLUX(I+1,J,K)-XFLUX(I,J,K)) . +(YFLUX(I,J+1,K)-YFLUX(I,J,K))) . /((H(I,J)+ETF(I,J))*ART(I,J)) C----------------------------------------------------------------------- C IMPOSE TRACER FLUX BOUNDARY CONDITIONS C----------------------------------------------------------------------- DO 150 N=1,NUMDBC ID=IDD(N) JD=JDD(N) DO 150 K=1,KBM1 IF(H(ID,JD)+ETF(ID,JD).LE.WETMIN.AND.QDIFF(N).GT.0.) THEN ! IF DIFFUSER ON DRY CELL (NG04132011) FF(ID,JD,K)=FDIF2(N,K) ! SET TO DIFFUSER'S SCALAR ELSE FF(ID,JD,K)=FF(ID,JD,K)+DTI2*FDIF2(N,K)*QDIFF(N)*RAMP* . VDDIST(N,K)/100./DZ(K)/((H(ID,JD)+ETF(ID,JD))*ART(ID,JD)) ENDIF 150 CONTINUE C----------------------------------------------------------------------- CNG04232008 IMPOSE ELEVATION POTENTIAL INTERIOR ADJUSTMENT C----------------------------------------------------------------------- !write (*,*)time,ff(iepc(1),jepc(1),10),qepcnudg(1) !write (*,*)time,ff(iepc(4),jepc(4),10),qepcnudg(4) DO N=1,NUMEPC I=IEPC(N) J=JEPC(N) DO K=1,KBM1 FF(I,J,K)=FF(I,J,K)+DTI2*F(I,J,K)*QEPCNUDG(N)/ . ((H(I,J)+ETF(I,J))*ART(I,J)) ENDDO ENDDO C----------------------------------------------------------------------- C RETURN END