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 ADVU(DRHOX,ADVUA,DTI2) C VERSION(11/16/90) C INCLUDE 'comdeck' SAVE C DIMENSION DRHOX(IM,JM,KB),XFLUX(IM,JM,KB),YFLUX(IM,JM,KB), . ADVUA(IM,JM),CURV4(IM,JM,KB) . ,XFLUXWAVE(IM,JM,KB),YFLUXWAVE(IM,JM,KB) !RMarsooli EQUIVALENCE (A,XFLUX),(C,YFLUX),(VH,CURV4) c!!!!!!!!!!!!!!!!!!!RMarsooli DO J=1,JM DO I=1,IM FSMRAD(I,J)=FSM(I,J)*DHMWAD(I,J)*FSMW(I,J) !IF(DT(I,J).LT.1.0.OR.WVHT(I,J).LT.0.30) FSMRAD(I,J)=0.0 IF(DT(I,J).LT.1.0) FSMRAD(I,J)=0.0 ENDDO ENDDO DO 95 K=1,KB DO 95 J=1,JM DO 95 I=1,IM XFLUXWAVE(I,J,K)=0.0 !RMarsooli YFLUXWAVE(I,J,K)=0.0 !RMarsooli UF (I,J,K)=0.0 CURV4(I,J,K)=0.0 XFLUX(I,J,K)=0.0 95 YFLUX(I,J,K)=0.0 IF (ADVECT.EQ.'NON-LINEAR') THEN DO 60 J=2,JMM1 DO 60 K=1,KBM1 DO 60 I=2,IMM1 CURV4(I,J,K)= . +.125*((V(I,J+1,K)+V(I,J,K))* . ((H2(I+1,J)*FSM(I+1,J)+H2(I,J)*FSM(I,J))/ . (FSM(I+1,J)+FSM(I,J)+1.E-30)- . (H2(I,J)*FSM(I,J)+H2(I-1,J)*FSM(I-1,J))/ . (FSM(I,J)+FSM(I-1,J)+1.E-30)) . -(U(I+1,J,K)+U(I,J,K))* . ((H1(I,J+1)*FSM(I,J+1)+H1(I,J)*FSM(I,J))/ . (FSM(I,J+1)+FSM(I,J)+1.E-30)- . (H1(I,J)*FSM(I,J)+H1(I,J-1)*FSM(I,J-1))/ . (FSM(I,J)+FSM(I,J-1)+1.E-30)) ) . /ART(I,J) . +.25*COR(I,J) 60 CONTINUE ELSE DO 62 J=2,JMM1 DO 62 K=1,KBM1 DO 62 I=2,IMM1 62 CURV4(I,J,K)=0.25*COR(I,J) END IF DO 65 J=1,JM DO 65 I=1,IM 65 CURV42D(I,J)=0.0 DO 70 J=1,JM DO 70 K=1,KBM1 DO 70 I=1,IM 70 CURV42D(I,J)=CURV42D(I,J)+CURV4(I,J,K)*DZ(K) C C-------- HORIZONTAL ADVECTION ----------------------------------------- IF (ADVECT.EQ.'NON-LINEAR') THEN DO 100 K=1,KBM1 DO 100 J=2,JM DO 100 I=2,IMM1 100 XFLUX(I,J,K)=0.25*(XMFL3D(I,J,K)+XMFL3D(I+1,J,K)) . *(U(I,J,K)+U(I+1,J,K)) cC UPWINDING, RMarsooli c UFACE=(DUMWAD(I,J)*U(I,J,K)+DUMWAD(I+1,J)*U(I+1,J,K)) c . /(DUMWAD(I,J)+DUMWAD(I+1,J)+SMALLVAL) c CONVECTEMP=H2(I,J)*DT(I,J) c IF(UFACE.GE.0.) THEN c XFLUX(I,J,K)=CONVECTEMP*U(I,J,K)*UFACE c ELSE c XFLUX(I,J,K)=CONVECTEMP*U(I+1,J,K)*UFACE c ENDIF c 100 CONTINUE DO 120 K=1,KBM1 DO 120 J=2,JMM1 DO 120 I=2,IMM1 120 YFLUX(I,J,K)=0.25*(YMFL3D(I-1,J,K)+YMFL3D(I,J,K)) . *(U(I,J-1,K)+U(I,J,K)) cC UPWINDING, RMarsooli c VFACE=(DVMWAD(I,J)*V(I,J,K)+DVMWAD(I-1,J)*V(I-1,J,K)) c . /(DVMWAD(I,J)+DVMWAD(I-1,J)+SMALLVAL) c CONVECTEMP=( (DHMWAD(I,J)*H1(I,J)+DHMWAD(I-1,J)*H1(I-1,J) c . +DHMWAD(I,J-1)*H1(I,J-1)+DHMWAD(I-1,J-1)*H1(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(VFACE.GE.0.) THEN c YFLUX(I,J,K)=CONVECTEMP*U(I,J-1,K)*VFACE c ELSE c YFLUX(I,J,K)=CONVECTEMP*U(I,J,K)*VFACE c ENDIF c 120 CONTINUE END IF C C---------WAVE RADIATION STRESS, RMarsooli, 2015 DO 21 K=1,KBM1 DO 21 J=2,JM DO 21 I=2,IMM1 21 XFLUXWAVE(I,J,K)=RAMP*DT(I,J)*Sxx(I,J,K) *H2(I,J) DO 22 K=1,KBM1 DO 22 J=2,JMM1 DO 22 I=2,IMM1 22 YFLUXWAVE(I,J,K)=RAMP*0.125 . * ( (DT(I,J)+DT(I-1,J))*(Sxy(I,J,K)+Sxy(I-1,J,K)) . +(DT(I,J-1)+DT(I-1,J-1))*(Sxy(I,J-1,K)+Sxy(I-1,J-1,K)) ) . *0.25*(H1(I,J)+H1(I-1,J)+H1(I,J-1)+H1(I-1,J-1)) C C-------- HORIZONTAL DIFFUSION ----------------------------------------- DO 110 K=1,KBM1 DO 110 J=1,JM DO 110 I=1,IM 110 PROD(I,J,K)=0.0 DO 500 J=2,JM DO 500 K=1,KBM1 DO 500 I=2,IM C C-------- THIS IS USED IN advv.f AS WELL AS HERE ----------------------- PROD(I,J,K)=.25*(DT(I,J)+DT(I-1,J)+DT(I,J-1)+DT(I-1,J-1)) . *(AAM(I,J,K)+AAM(I-1,J,K)+AAM(I,J-1,K)+AAM(I-1,J-1,K)) 500 CONTINUE C DO 700 J=2,JMM1 DO 700 K=1,KBM1 DO 700 I=2,IMM1 XFLUX(I,J,K)=XFLUX(I,J,K) . -DT(I,J)*AAM(I,J,K)*2.*(UB(I+1,J,K)-UB(I,J,K))*H2(I,J)/H1(I,J) YFLUX(I,J,K)=YFLUX(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*(H1(I,J)+H1(I-1,J)+H1(I,J-1)+H1(I-1,J-1)) . *DUM(I,J)*DUM(I,J-1) * DVM(I,J)*DVM(I-1,J) ! FOR THIN DAM . *DUMWAD(I,J)*DUMWAD(I,J-1) * DVMWAD(I,J)*DVMWAD(I-1,J) ! NG09302011, not in original NYHOPS c . *DUM(I,J)*DUM(I,J-1) 700 CONTINUE C C-------- DO VERTICAL ADVECTION ---------------------------------------- C-------- ADD HORIZONTAL ADVECTION ------------------------------------- IF (ADVECT.EQ.'NON-LINEAR') THEN DO 130 J=1,JM DO 130 I=1,IM 130 UF(I,J,1)=0.0 DO 140 K=2,KBM1 DO 140 J=2,JMM1 DO 140 I=2,IMM1 140 UF(I,J,K)=.25*(W(I,J,K)+W(I-1,J,K))*(U(I,J,K)+U(I,J,K-1)) DO 160 J=1,JM DO 160 I=1,IM 160 UF(I,J,KB)=0.0 DO 145 J=2,JMM1 DO 145 K=1,KBM1 DO 145 I=2,IMM1 145 UF(I,J,K)=DZR(K)*(UF(I,J,K)-UF(I,J,K+1))*ARU(I,J) END IF C DO 146 K=1,KBM1 DO 146 J=2,JMM1 DO 146 I=2,IMM1 146 UF(I,J,K)=UF(I,J,K) . +(XFLUX(I,J,K)-XFLUX(I-1,J,K)) . +(XFLUXWAVE(I,J,K)-XFLUXWAVE(I-1,J,K)) !RMarsooli . *FSMRAD(I,J)*FSMRAD(I-1,J) . +(YFLUX(I,J+1,K)-YFLUX(I,J,K)) . +(YFLUXWAVE(I,J+1,K)-YFLUXWAVE(I,J,K)) !RMarsooli . *FSMRAD(I,J)*FSMRAD(I-1,J) . *FSMRAD(I,J+1)*FSMRAD(I-1,J+1) . *FSMRAD(I,J-1)*FSMRAD(I-1,J-1) C C-------- SAVE HORIZONTAL ADVECTION & DIFFUSION FLUXES ----------------- C-------------------- FOR EXTERNAL MODE -------------------------------- DO 790 J=1,JM DO 790 I=1,IM 790 ADVUA(I,J)=0.0 DO 800 J=2,JMM1 DO 800 K=1,KBM1 DO 800 I=2,IMM1 800 ADVUA(I,J)=ADVUA(I,J)+DZ(K)*UF(I,J,K) C DO 310 N=1,NUMEBC IE=IETA(N) JE=JETA(N) IC=ICON(N) JC=JCON(N) IF(FSM(IE+1,JE).EQ.0.0.AND.JE.EQ.JC) THEN CURV42D(IE,JE)=CURV42D(IE-1,JE) ADVUA (IE,JE)=0.0 DO 320 K=1,KBM1 CURV4(IE,JE,K)=CURV4(IE-1,JE,K) 320 UF (IE,JE,K)=0.0 ELSE IF(FSM(IE-1,JE).EQ.0.0.AND.JE.EQ.JC) THEN CURV42D(IE,JE)=CURV42D(IE+1,JE) ADVUA (IE+1,JE)=0.0 DO 330 K=1,KBM1 CURV4(IE,JE,K)=CURV4(IE+1,JE,K) 330 UF (IE+1,JE,K)=0.0 ENDIF 310 CONTINUE CNG2012 CORNERS FIX - PATCH FOR CORIOLIS PROBLEM DO N=1,NUMECR Do 360 K = 1, KBM1 IF(NCHARD(N).EQ.1) THEN ! Bottom Left Corner V(ICHARD(N),JCHARD(N)+1,K) = + .5*(V(ICHARD(N)+1,JCHARD(N)+1,K)+V(ICHARD(N),JCHARD(N)+2,K)) ELSEIF(NCHARD(N).EQ.2) THEN ! Bottom Right Corner V(ICHARD(N),JCHARD(N)+1,K) = + .5*(V(ICHARD(N)-1,JCHARD(N)+1,K)+V(ICHARD(N),JCHARD(N)+2,K)) ELSEIF(NCHARD(N).EQ.3) THEN ! Top Right Corner V(ICHARD(N),JCHARD(N),K) = + .5*(V(ICHARD(N)-1,JCHARD(N),K)+V(ICHARD(N),JCHARD(N)-1,K)) ELSE ! Top Left Corner V(ICHARD(N),JCHARD(N),K) = + .5*(V(ICHARD(N)+1,JCHARD(N),K)+V(ICHARD(N),JCHARD(N)-1,K)) ENDIF 360 CONTINUE ENDDO C C-------- -FVD + GDEG/H1 + BAROCLINIC TERM ----------------------------- DO 150 J=2,JMM1 DO 150 K=1,KBM1 DO 150 I=3,IMM1 150 UF(I,J,K)=UF(I,J,K) . -ARU(I,J)*(CURV4(I,J,K)*DT(I,J)*(V(I,J+1,K)+V(I,J,K)) . +CURV4(I-1,J,K)*DT(I-1,J)*(V(I-1,J+1,K)+V(I-1,J,K))) . +GRAV*.25*(DT(I,J)+DT(I-1,J)) . *.5*(EGF(I,J)-EGF(I-1,J)+EGB(I,J)-EGB(I-1,J)) . *(H2(I,J)+H2(I-1,J)) . +GRAV*.25*(DT(I,J)+DT(I-1,J)) . *(DATUM(I,J)-DATUM(I-1,J))*RAMP*(H2(I,J)+H2(I-1,J)) . +DRHOX(I,J,K) cqa Add Atmospheric Pressure Terms by Quamrul QA 4/20/00 * +ARU(I,J)*(DT(I,J)+DT(I-1,J))*(PATM(I,J)-PATM(I-1,J)) CNG * /RHO0/(H1(I,J)+H1(I-1,J)) * *RAMP*PATMOPT/RHO0/(H1(I,J)+H1(I-1,J)) !WETLAND VEGETATION DRAG FORCE-RMarsooli_Jan2015 . -(DHMWAD(I,J)*DVEG3D(I,J,K)+DHMWAD(I-1,J)*DVEG3D(I-1,J,K)) . /(DHMWAD(I,J)+DHMWAD(I-1,J)+SMALLVAL)*U(I,J,K) !SIDE WALL FRICTION FORCE-RMarsooli-Jan2015 . -DWALLU3D(I,J,K)*U(I,J,K) C C-------- STEP FORWARD IN TIME ----------------------------------------- DO 190 J=2,JMM1 DO 190 K=1,KBM1 DO 190 I=3,IMM1 190 UF(I,J,K)= . ( (H(I,J)+ETB(I,J)+H(I-1,J)+ETB(I-1,J))*ARU(I,J)*UB(I,J,K) . -2.*DTI2*UF(I,J,K) !WETLAND VEGETATION INERTIA FORCE-RMarsooli_Jan2015 . +2.*DTI2*UF(I,J,K)*(1.0-1.0/ . (1.0-0.5*(IVEG3D(I,J,K)+IVEG3D(I-1,J,K)))) ) . / ((H(I,J)+ETF(I,J)+H(I-1,J)+ETF(I-1,J))*ARU(I,J)) RETURN END