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 ADVV(DRHOY,ADVVA,DTI2) C VERSION(11/16/90) C INCLUDE 'comdeck' SAVE C DIMENSION DRHOY(IM,JM,KB),XFLUX(IM,JM,KB),YFLUX(IM,JM,KB), . ADVVA(IM,JM),CURV4(IM,JM,KB) . ,XFLUXWAVE(IM,JM,KB),YFLUXWAVE(IM,JM,KB) !RMarsooli EQUIVALENCE (A,XFLUX),(C,YFLUX),(VH,CURV4) C DO 10 K=1,KB DO 10 J=1,JM DO 10 I=1,IM XFLUXWAVE(I,J,K)=0.0 !RMarsooli YFLUXWAVE(I,J,K)=0.0 !RMarsooli VF (I,J,K)=0.0 XFLUX(I,J,K)=0.0 10 YFLUX(I,J,K)=0.0 C C-------- HORIZONTAL ADVECTION ----------------------------------------- IF (ADVECT.EQ.'NON-LINEAR') THEN DO 20 K=1,KBM1 DO 20 J=2,JMM1 DO 20 I=2,IM 20 XFLUX(I,J,K)=0.25*(XMFL3D(I,J-1,K)+XMFL3D(I,J,K)) . *(V(I-1,J,K)+V(I,J,K)) cC UPWINDING, RMarsooli c UFACE=(DUMWAD(I,J)*U(I,J,K)+DUMWAD(I,J-1)*U(I,J-1,K)) c . /(DUMWAD(I,J)+DUMWAD(I,J-1)+SMALLVAL) c CONVECTEMP=( (DHMWAD(I,J)*H2(I,J)+DHMWAD(I-1,J)*H2(I-1,J) c . +DHMWAD(I,J-1)*H2(I,J-1)+DHMWAD(I-1,J-1)*H2(I-1,J-1)) c . /(DHMWAD(I,J)+DHMWAD(I-1,J)+DHMWAD(I,J-1) c . +DHMWAD(I-1,J-1)+SMALLVAL) ) c . * ( (DHMWAD(I,J)*DT(I,J)+DHMWAD(I-1,J)*DT(I-1,J) c . +DHMWAD(I,J-1)*DT(I,J-1)+DHMWAD(I-1,J-1)*DT(I-1,J-1)) c . /(DHMWAD(I,J)+DHMWAD(I-1,J)+DHMWAD(I,J-1) c . +DHMWAD(I-1,J-1)+SMALLVAL) ) c IF(UFACE.GE.0.) THEN c XFLUX(I,J,K)=CONVECTEMP*V(I-1,J,K)*UFACE c ELSE c XFLUX(I,J,K)=CONVECTEMP*V(I,J,K)*UFACE c ENDIF c 20 CONTINUE DO 30 K=1,KBM1 DO 30 J=2,JMM1 DO 30 I=2,IMM1 30 YFLUX(I,J,K)=0.25*(YMFL3D(I,J,K)+YMFL3D(I,J+1,K)) . *(V(I,J,K)+V(I,J+1,K)) cC UPWINDING, RMarsoolic VFACE=0.5*(V(I,J,K)+V(I,J+1,K)) c VFACE=(DVMWAD(I,J)*V(I,J,K)+DVMWAD(I,J+1)*V(I,J+1,K)) c . /(DVMWAD(I,J)+DVMWAD(I,J+1)+SMALLVAL) c CONVECTEMP=H1(I,J)*DT(I,J) c IF(VFACE.GE.0.) THEN c YFLUX(I,J,K)=CONVECTEMP*V(I,J,K)*VFACE c ELSE c YFLUX(I,J,K)=CONVECTEMP*V(I,J+1,K)*VFACE c ENDIF c 30 CONTINUE END IF C C---------WAVE RADIATION STRESS, RMarsooli, 2015 DO 21 K=1,KBM1 DO 21 J=2,JMM1 DO 21 I=2,IM 21 XFLUXWAVE(I,J,K)=RAMP*0.125 . * ( (DT(I,J)+DT(I,J-1))*(Sxy(I,J,K)+Sxy(I,J-1,K)) . +(DT(I-1,J)+DT(I-1,J-1))*(Sxy(I-1,J,K)+Sxy(I-1,J-1,K)) ) . *0.25*(H2(I,J)+H2(I-1,J)+H2(I,J-1)+H2(I-1,J-1)) DO 22 K=1,KBM1 DO 22 J=2,JMM1 DO 22 I=2,IMM1 22 YFLUXWAVE(I,J,K)=RAMP*DT(I,J)*Syy(I,J,K) *H1(I,J) C C-------- HORIZONTAL DIFFUSION ----------------------------------------- C-------- NOTE: PROD = DT*AAM LEFT OVER FROM ADVU --------------------- DO 40 J=2,JMM1 DO 40 K=1,KBM1 DO 40 I=2,IMM1 XFLUX(I,J,K)=XFLUX(I,J,K) . -PROD(I,J,K)*((UB(I,J,K)-UB(I,J-1,K)) . /(H2(I,J)+H2(I-1,J)+H2(I,J-1)+H2(I-1,J-1)) . +(VB(I,J,K)-VB(I-1,J,K)) . /(H1(I,J)+H1(I-1,J)+H1(I,J-1)+H1(I-1,J-1))) . *0.25*(H2(I,J)+H2(I-1,J)+H2(I,J-1)+H2(I-1,J-1)) . *DVM(I,J)*DVM(I-1,J)*DUM(I,J)*DUM(I,J-1) ! FOR THIN DAM . *DVMWAD(I,J)*DVMWAD(I-1,J) * DUMWAD(I,J)*DUMWAD(I,J-1) ! NG09302011 Not in original NYHOPS c . *DVM(I,J)*DVM(I-1,J) YFLUX(I,J,K)=YFLUX(I,J,K) . -DT(I,J)*AAM(I,J,K)*2.0*(VB(I,J+1,K)-VB(I,J,K))*H1(I,J)/H2(I,J) 40 CONTINUE C C-------- DO VERTICAL ADVECTION ---------------------------------------- C-------- ADD HORIZONTAL ADVECTION ------------------------------------- IF (ADVECT.EQ.'NON-LINEAR') THEN DO 50 J=1,JM DO 50 I=1,IM 50 VF(I,J,1)=0.0 DO 60 K=2,KBM1 DO 60 J=3,JMM1 DO 60 I=2,IMM1 60 VF(I,J,K)=.25*(W(I,J,K)+W(I,J-1,K))*(V(I,J,K)+V(I,J,K-1)) DO 70 J=1,JM DO 70 I=1,IM 70 VF(I,J,KB)=0.0 DO 80 J=3,JMM1 DO 80 K=1,KBM1 DO 80 I=2,IMM1 80 VF(I,J,K)=DZR(K)*(VF(I,J,K)-VF(I,J,K+1))*ARV(I,J) END IF C DO 90 K=1,KBM1 DO 90 J=3,JMM1 DO 90 I=2,IMM1 90 VF(I,J,K)=VF(I,J,K) . +(XFLUX(I+1,J,K)-XFLUX(I,J,K)) . +(XFLUXWAVE(I+1,J,K)-XFLUXWAVE(I,J,K)) !RMarsooli . *FSMRAD(I,J)*FSMRAD(I,J-1) . *FSMRAD(I+1,J)*FSMRAD(I+1,J-1) . *FSMRAD(I-1,J)*FSMRAD(I-1,J-1) . +(YFLUX(I,J,K)-YFLUX(I,J-1,K)) . +(YFLUXWAVE(I,J,K)-YFLUXWAVE(I,J-1,K)) !RMarsooli . *FSMRAD(I,J)*FSMRAD(I,J-1) C C-------- SAVE HORIZONTAL ADVECTION & DIFFUSION FLUXES ----------------- C------------------- FOR EXTERNAL MODE --------------------------------- DO 790 J=1,JM DO 790 I=1,IM 790 ADVVA(I,J)=0.0 DO 800 J=3,JMM1 DO 800 K=1,KBM1 DO 800 I=2,IMM1 800 ADVVA(I,J)=ADVVA(I,J)+DZ(K)*VF(I,J,K) C DO 305 N=1,NUMEBC IE=IETA(N) JE=JETA(N) IC=ICON(N) JC=JCON(N) IF(FSM(IE,JE+1).EQ.0.0.AND.IE.EQ.IC) THEN CURV42D(IE,JE)=CURV42D(IE,JE-1) ADVVA(IE,JE)=0.0 DO 310 K=1,KBM1 CURV4(IE,JE,K)=CURV4(IE,JE-1,K) 310 VF(IE,JE,K)=0.0 ELSE IF(FSM(IE,JE-1).EQ.0.0.AND.IE.EQ.IC) THEN CURV42D(IE,JE)=CURV42D(IE,JE+1) ADVVA(IE,JE+1)=0.0 DO 315 K=1,KBM1 CURV4(IE,JE,K)=CURV4(IE,JE+1,K) 315 VF(IE,JE+1,K)=0.0 ENDIF 305 CONTINUE CNG2012 CORNERS FIX - PATCH FOR CORIOLIS PROBLEM DO N=1,NUMECR DO 270 K=1,KBM1 IF(NCHARD(N).EQ.1) THEN ! Bottom Left Corner U(ICHARD(N)+1,JCHARD(N),K) = + .5*(U(ICHARD(N)+1,JCHARD(N)+1,K)+U(ICHARD(N)+2,JCHARD(N),K)) ELSEIF(NCHARD(N).EQ.2) THEN ! Bottom Right Corner U(ICHARD(N),JCHARD(N),K) = + .5*(U(ICHARD(N)-1,JCHARD(N),K)+U(ICHARD(N),JCHARD(N)+1,K)) ELSEIF(NCHARD(N).EQ.3) THEN ! Top Right Corner U(ICHARD(N),JCHARD(N),K) = + .5*(U(ICHARD(N)-1,JCHARD(N),K)+U(ICHARD(N),JCHARD(N)-1,K)) ELSE ! Top Left Corner U(ICHARD(N)+1,JCHARD(N),K) = + .5*(U(ICHARD(N)+2,JCHARD(N),K)+U(ICHARD(N)+1,JCHARD(N)-1,K)) ENDIF 270 CONTINUE ENDDO C C--------- +FUD + GDEG/H2 + BAROCLINIC TERM ---------------------------- DO 340 J=3,JMM1 DO 340 K=1,KBM1 DO 340 I=2,IMM1 340 VF(I,J,K)=VF(I,J,K) . +ARV(I,J)*(CURV4(I,J,K)*DT(I,J)*(U(I+1,J,K)+U(I,J,K)) . +CURV4(I,J-1,K)*DT(I,J-1)*(U(I+1,J-1,K)+U(I,J-1,K))) . +GRAV*.25*(DT(I,J)+DT(I,J-1)) . *.5*(EGF(I,J)-EGF(I,J-1)+EGB(I,J)-EGB(I,J-1)) . *(H1(I,J)+H1(I,J-1)) . +GRAV*.25*(DT(I,J)+DT(I,J-1)) . *(DATUM(I,J)-DATUM(I,J-1))*RAMP*(H1(I,J)+H1(I,J-1)) . +DRHOY(I,J,K) cqa Add Atmospheric Pressure Terms by Quamrul QA 4/20/00 * +ARV(I,J)*(DT(I,J)+DT(I,J-1))*(PATM(I,J)-PATM(I,J-1)) CNG * /RHO0/(H2(I,J)+H2(I,J-1)) * *RAMP*PATMOPT/RHO0/(H2(I,J)+H2(I,J-1)) !WETLAND VEGETATION DRAG FORCE-RMarsooli_Jan2015 . -(DHMWAD(I,J)*DVEG3D(I,J,K)+DHMWAD(I,J-1)*DVEG3D(I,J-1,K)) . /(DHMWAD(I,J)+DHMWAD(I,J-1)+SMALLVAL)*V(I,J,K) !SIDE WALL FRICTION FORCE-RMarsooli_Jan2015 . -DWALLV3D(I,J,K)*V(I,J,K) C C-------- STEP FORWARD IN TIME ----------------------------------------- DO 390 J=3,JMM1 DO 390 K=1,KBM1 DO 390 I=2,IMM1 390 VF(I,J,K)= . ((H(I,J)+ETB(I,J)+H(I,J-1)+ETB(I,J-1))*ARV(I,J)*VB(I,J,K) . -2.*DTI2*VF(I,J,K) !WETLAND VEGETATION INERTIA FORCE-RMarsooli_Jan2015 . +2.*DTI2*VF(I,J,K)*(1.0-1.0/ . (1.0-0.5*(IVEG3D(I,J,K)+IVEG3D(I,J-1,K)))) ) . /((H(I,J)+ETF(I,J)+H(I,J-1)+ETF(I,J-1))*ARV(I,J)) RETURN END