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 EXTRNL(ADVUA,ADVVA,TRNU,TRNV,DTE2,DTI2) cpl SUBROUTINE EXTRNL(ADVUA,ADVVA,TRNU,TRNV,DTE2,BFRIC,DTI2) C VERSION(11/16/90) C INCLUDE 'comdeck' SAVE REAL*8 QEPC C C----------------------------------------------------------------------- C CALCULATE FREE SURFACE HEIGHT C----------------------------------------------------------------------- C DIMENSION ADVUA(IM,JM),ADVVA(IM,JM),TRNU(IM,JM),TRNV(IM,JM) C DO 400 J=2,JMM1 DO 400 I=2,IM 400 FLUXUA(I,J)=.25*(D(I,J)+D(I-1,J))*(H2(I,J)+H2(I-1,J))*UA(I,J) DO 405 J=2,JM DO 405 I=2,IMM1 405 FLUXVA(I,J)=.25*(D(I,J)+D(I,J-1))*(H1(I,J)+H1(I,J-1))*VA(I,J) C cwad DO J=2,JMM1 cwad DO I=2,IMM1 cwad IF(D(I,J).LE.WETMIN) THEN cwad ELF(I,J)=EL(I,J)-0.5*DTE2* cwad . (FLUXUA(I+1,J)-FLUXUA(I,J)+FLUXVA(I,J+1)-FLUXVA(I,J))/ cwad & ART(I,J) cwad ELSE cwad ELF(I,J)=ELB(I,J)-DTE2* cwad . (FLUXUA(I+1,J)-FLUXUA(I,J)+FLUXVA(I,J+1)-FLUXVA(I,J))/ cwad & ART(I,J) cwad ENDIF cwad ENDDO cwad ENDDO DO 410 J=2,JMM1 DO 410 I=2,IMM1 410 ELF(I,J)=ELB(I,J)-DTE2* . (FLUXUA(I+1,J)-FLUXUA(I,J)+FLUXVA(I,J+1)-FLUXVA(I,J))/ART(I,J) CNG12162010 INCLUDE RAIN/EVAP FLUX . -DTE2*(QEVAP2D(I,J)-QPREC2D(I,J))*RAMP*FSM(I,J) C----------------------------------------------------------------------- C IMPOSE MASS FLUX BOUNDARY CONDITIONS C----------------------------------------------------------------------- DO 406 N=1,NUMDBC ID=IDD(N) JD=JDD(N) ELF(ID,JD)=ELF(ID,JD)+ DTE2*QDIFF(N)*RAMP /ART(ID,JD) 406 CONTINUE C C----------------------------------------------------------------------- CNG04232008 IMPOSE ELEVATION POTENTIAL INTERIOR CONDITIONS C----------------------------------------------------------------------- IF (OPTEPC(1:8).EQ.'PCLAMPED') THEN DO N=1,NUMEPC I=IEPC(N) J=JEPC(N) IF (EPCOR(N).NE.-99.AND.D(I,J).GT.WETMIN) THEN CPHL=DUM(I,J)*DTE2/(0.5*(H1(I,J)+H1(I-1,J)))* + SQRT(GRAV*0.5*(D(I,J)+D(I-1,J))) CPHR=DUM(I+1,J)*DTE2/(0.5*(H1(I,J)+H1(I+1,J)))* + SQRT(GRAV*0.5*(D(I,J)+D(I+1,J))) CPHU=DVM(I,J+1)*DTE2/(0.5*(H2(I,J)+H2(I,J+1)))* + SQRT(GRAV*0.5*(D(I,J)+D(I,J+1))) CPHD=DVM(I,J)*DTE2/(0.5*(H2(I,J)+H2(I,J-1)))* + SQRT(GRAV*0.5*(D(I,J)+D(I,J-1))) ELFTEMP=(ELB(I,J)+CPHL*(2.*EL(I-1,J)-ELB(I,J))+DTE2* + (EPCOR(N)-ELB(I,J))/TNUDGEPC(N))/(1.+CPHL) + + (ELB(I,J)+CPHR*(2.*EL(I+1,J)-ELB(I,J))+DTE2* + (EPCOR(N)-ELB(I,J))/TNUDGEPC(N))/(1.+CPHR) + + (ELB(I,J)+CPHU*(2.*EL(I,J+1)-ELB(I,J))+DTE2* + (EPCOR(N)-ELB(I,J))/TNUDGEPC(N))/(1.+CPHU) + + (ELB(I,J)+CPHD*(2.*EL(I,J-1)-ELB(I,J))+DTE2* + (EPCOR(N)-ELB(I,J))/TNUDGEPC(N))/(1.+CPHD) ELFTEMP=0.25*ELFTEMP QEPC=RAMP*ART(I,J)*(ELFTEMP-ELF(I,J))/DTE2 ! Q is positive incoming (if ELFTEMP>ELF) QEPCNUDF(N)=QEPCNUDF(N)+QEPC*ISPI ! AVERAGE OVER THE EXTERNAL MODE ELF(I,J)=ELF(I,J) + DTE2*QEPC/ART(I,J) ELSE QEPCNUDF(N)=0.0 ! ALLOW NO EPC FLUX ENDIF ENDDO ELSE ! FOR NUDGING DO N=1,NUMEPC I=IEPC(N) J=JEPC(N) IF (EPCOR(N).NE.-99.AND.D(I,J).GT.WETMIN) THEN CNG Note: I tested the one below as well, (epcor applied on a dry CNG element, if positive). It also worked, but it was not as smooth as the CNG above (only allowing lateral, not EPCOR-related wetting). CNG IF ( (EPCOR(N).NE.-99).AND. CNG . (D(I,J).GT.WETMIN.OR.EPCOR(N).GT.ELB(I,J)) ) THEN QEPC=RAMP*ART(I,J)*(EPCOR(N)-ELB(I,J))/TNUDGEPC(N) ! Q is positive incoming (if EPCOR(N)>ELB) QEPCNUDF(N)=QEPCNUDF(N)+QEPC*ISPI ! AVERAGE OVER THE EXTERNAL MODE ELF(I,J)=ELF(I,J) + DTE2*QEPC/ART(I,J) ELSE QEPCNUDF(N)=0.0 ! ALLOW NO EPC FLUX ENDIF ENDDO ENDIF C CALL BCOND(1,DTI2,0) C C Introduced Inverted Reid and Bodine (IRANDB) and Optimized Clamped C and MIXED TYPE BOUNDARY CONDITIONS In Barotropic or external Mode C C MODIFIED BY Quamrul QA 01/03/99 C cqa IF (BCTYPE.EQ.'RANDB '.OR.BCTYPE.EQ.'SHULMAN') THEN cqa CALL SHULMAN cqa ELSE IF (BCTYPE.EQ.'PCLAMP ') THEN cqa CALL PCLAMP cqa ENDIF C Modified by Quamrul QA adopted from JKL 8/23/00 IF (BCTYPE.EQ.'OCLAMP '.OR.BCTYPE.EQ.'IRANDB '.OR. + BCTYPE.EQ.'MIXED ') CALL OCLAMP(CLAM) IF (BCTYPE.EQ.'PCLAMP ') CALL PCLAMP C CALL ADVAVE(ADVUA,ADVVA) C DO 95 J=1,JM DO 95 I=1,IM 95 CURV42D(I,J)=.25*COR(I,J) C C----------------------------------------------------------------------- C COMPUTE BOTTOM FRICTION AND CURVATURE TERMS C IF RUN IS BAROTROPIC MODE ONLY C----------------------------------------------------------------------- C IF(TOR.NE.'BAROTROPIC') GO TO 5000 C DO 60 J=2,JMM1 DO 60 I=2,IMM1 CURV42D(I,J)= CURV42D(I,J) . +.125*((VA(I,J+1)+VA(I,J))* . ((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)) . -(UA(I+1,J)+UA(I,J))* . ((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) 60 CONTINUE C DO 100 J=2,JMM1 DO 100 I=2,IMM1 CNG04232008 TXBOT(I,J)=-BFRIC TXBOT(I,J)=-0.5*(CBC(I,J)+CBC(I-1,J)) . *SQRT(UAB(I,J)**2+(.25*(VAB(I,J) . +VAB(I,J+1)+VAB(I-1,J)+VAB(I-1,J+1)))**2)*UAB(I,J) TXICE(I,J)=0.5*(AREAICE(I,J)/100.*VARWIF(I,J) ! NG FOR ICE . +AREAICE(I-1,J)/100.*VARWIF(I-1,J)) ! NG FOR ICE . *SQRT(UAB(I,J)**2+(.25*(VAB(I,J) . +VAB(I,J+1)+VAB(I-1,J)+VAB(I-1,J+1)))**2)*UAB(I,J) 100 CONTINUE DO 102 J=2,JMM1 DO 102 I=2,IMM1 CNG04232008 TYBOT(I,J)=-BFRIC TYBOT(I,J)=-0.5*(CBC(I,J)+CBC(I,J-1)) . *SQRT((.25*(UAB(I,J)+UAB(I+1,J) . +UAB(I,J-1)+UAB(I+1,J-1)))**2+VAB(I,J)**2)*VAB(I,J) TYICE(I,J)=-0.5*(AREAICE(I,J)/100.*VARWIF(I,J) ! NG FOR ICE . +AREAICE(I,J-1)/100.*VARWIF(I,J-1)) ! NG FOR ICE . *SQRT((.25*(UAB(I,J)+UAB(I+1,J) . +UAB(I,J-1)+UAB(I+1,J-1)))**2+VAB(I,J)**2)*VAB(I,J) 102 CONTINUE CNG INCLUDE SIMILARLY TO 3D TREATMENT CNG WHETHER THIS IS A GOOD IDEA NEEDS FURTHER INVESTIGATION! 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) TYBOT(IE,JE)=TYBOT(IE-1,JE) ! NG same as internal bcond, cosmetic ELSE IF(FSM(IE-1,JE).EQ.0.0.AND.JE.EQ.JC) THEN CURV42D(IE,JE)=CURV42D(IE+1,JE) TYBOT(IE,JE)=TYBOT(IE+1,JE) ! NG same as internal bcond, cosmetic ELSE IF(FSM(IE,JE+1).EQ.0.0.AND.IE.EQ.IC) THEN CURV42D(IE,JE)=CURV42D(IE,JE-1) TXBOT(IE,JE)=TXBOT(IE,JE-1) ! NG same as internal bcond, cosmetic ELSE IF(FSM(IE,JE-1).EQ.0.0.AND.IE.EQ.IC) THEN CURV42D(IE,JE)=CURV42D(IE,JE+1) TXBOT(IE,JE)=TXBOT(IE,JE+1) ! NG same as internal bcond, cosmetic ENDIF 310 CONTINUE 5000 CONTINUE C C----------------------------------------------------------------------- C CALCULATE U COMPONENT OF VELOCITY C----------------------------------------------------------------------- C CNG2012 CORNERS FIX - NO NEED TO HARDWIRE ANYMORE DO N=1,NUMECR IF(NCHARD(N).EQ.1) THEN ! Bottom Left Corner Va(ICHARD(N),JCHARD(N)+1) = + .5*(Va(ICHARD(N)+1,JCHARD(N)+1)+Va(ICHARD(N),JCHARD(N)+2)) Ua(ICHARD(N)+1,JCHARD(N)) = + .5*(Ua(ICHARD(N)+1,JCHARD(N)+1)+Ua(ICHARD(N)+2,JCHARD(N))) ELSEIF(NCHARD(N).EQ.2) THEN ! Bottom Right Corner Va(ICHARD(N),JCHARD(N)+1) = + .5*(Va(ICHARD(N)-1,JCHARD(N)+1)+Va(ICHARD(N),JCHARD(N)+2)) Ua(ICHARD(N),JCHARD(N)) = + .5*(Ua(ICHARD(N)-1,JCHARD(N))+Ua(ICHARD(N),JCHARD(N)+1)) ELSEIF(NCHARD(N).EQ.3) THEN ! Top Right Corner Va(ICHARD(N),JCHARD(N)) = + .5*(Va(ICHARD(N)-1,JCHARD(N))+Va(ICHARD(N),JCHARD(N)-1)) Ua(ICHARD(N),JCHARD(N)) = + .5*(Ua(ICHARD(N)-1,JCHARD(N))+Ua(ICHARD(N),JCHARD(N)-1)) ELSE ! Top Left Corner Va(ICHARD(N),JCHARD(N)) = + .5*(Va(ICHARD(N)+1,JCHARD(N))+Va(ICHARD(N),JCHARD(N)-1)) Ua(ICHARD(N)+1,JCHARD(N)) = + .5*(Ua(ICHARD(N)+2,JCHARD(N))+Ua(ICHARD(N)+1,JCHARD(N)-1)) ENDIF ENDDO DO 420 J=2,JMM1 DO 420 I=3,IMM1 420 UAF(I,J)=ADVUA(I,J) . -ARU(I,J)*( CURV42D(I,J)*D(I,J)*(VA(I,J+1)+VA(I,J)) . +CURV42D(I-1,J)*D(I-1,J)*(VA(I-1,J+1)+VA(I-1,J)) ) . +.25*GRAV*(H2(I,J)+H2(I-1,J))*(D(I,J)+D(I-1,J)) . *( (1.-2.*THETA)*(EL(I,J)-EL(I-1,J)) . +THETA*(ELB(I,J)-ELB(I-1,J)+ELF(I,J)-ELF(I-1,J)) ) . +.25*GRAV*(H2(I,J)+H2(I-1,J))*(D(I,J)+D(I-1,J)) . *(DATUM(I,J)-DATUM(I-1,J))*RAMP . +RAMP*TRNU(I,J) CNG . -ARU(I,J)*(-.5*(TXSURF(I,J)+TXSURF(I-1,J))+TXBOT(I,J) ) . -ARU(I,J)* . ( -.5*((100.-AREAICE(I,J))/100.*TXSURF(I,J)+ ! NG FOR ICE . (100.-AREAICE(I-1,J))/100.*TXSURF(I-1,J) ) ! NG FOR ICE . -TXICE(I,J) ! NG FOR ICE . +TXBOT(I,J) ) cqa Add Atmospheric Pressure Terms by Quamrul QA 4/20/00 * +ARU(I,J)*(D(I,J)+D(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)*DVEG2D(I,J)+DHMWAD(I-1,J)*DVEG2D(I-1,J)) . /(DHMWAD(I,J)+DHMWAD(I-1,J)+SMALLVAL)*UA(I,J) !SIDE WALL FRICTION FORCE-RMarsooli_Jan2015 . -DWALLU2D(I,J)*UA(I,J) DO 425 J=2,JMM1 DO 425 I=3,IMM1 425 UAF(I,J)= . ( (H(I,J)+ELB(I,J)+H(I-1,J)+ELB(I-1,J))*ARU(I,J)*UAB(I,J) . -4.*DTE*UAF(I,J) !WETLAND VEGETATION INERTIA FORCE-RMarsooli_Jan2015 . +4.*DTE*UAF(I,J)*(1.0-1.0/ . (1.0-0.5*(IVEG2D(I,J)+IVEG2D(I-1,J)))) ) . /((H(I,J)+ELF(I,J)+H(I-1,J)+ELF(I-1,J))*ARU(I,J)) C C----------------------------------------------------------------------- C CALCULATE V COMPONENT OF VELOCITY C----------------------------------------------------------------------- C DO 430 J=3,JMM1 DO 430 I=2,IMM1 430 VAF(I,J)=ADVVA(I,J) . +ARV(I,J)*( CURV42D(I,J)*D(I,J)*(UA(I+1,J)+UA(I,J)) . +CURV42D(I,J-1)*D(I,J-1)*(UA(I+1,J-1)+UA(I,J-1)) ) . +.25*GRAV*(H1(I,J)+H1(I,J-1))*(D(I,J)+D(I,J-1)) . *( (1.-2.*THETA)*(EL(I,J)-EL(I,J-1)) . +THETA*(ELB(I,J)-ELB(I,J-1)+ELF(I,J)-ELF(I,J-1)) ) . +.25*GRAV*(H1(I,J)+H1(I,J-1))*(D(I,J)+D(I,J-1)) . *(DATUM(I,J)-DATUM(I,J-1))*RAMP . +RAMP*TRNV(I,J) CNG . + ARV(I,J)*( .5*(TYSURF(I,J)+TYSURF(I,J-1))-TYBOT(I,J) ) . + ARV(I,J)* . ( .5*((100.-AREAICE(I,J))/100.*TYSURF(I,J)+ ! NG FOR ICE . (100.-AREAICE(I,J-1))/100.*TYSURF(I,J-1) ) ! NG FOR ICE . +TYICE(I,J) ! NG FOR ICE . -TYBOT(I,J) ) cqa Add Atmospheric Pressure Terms by Quamrul QA 4/20/00 * +ARV(I,J)*(D(I,J)+D(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)*DVEG2D(I,J)+DHMWAD(I,J-1)*DVEG2D(I,J-1)) . /(DHMWAD(I,J)+DHMWAD(I,J-1)+SMALLVAL)*VA(I,J) !SIDE WALL FRICTION FORCE-RMarsooli_Jan2015 . -DWALLV2D(I,J)*VA(I,J) DO 435 J=3,JMM1 DO 435 I=2,IMM1 435 VAF(I,J)= . ( (H(I,J)+ELB(I,J)+H(I,J-1)+ELB(I,J-1))*VAB(I,J)*ARV(I,J) . -4.*DTE*VAF(I,J) !WETLAND VEGETATION INERTIA FORCE-RMarsooli_Jan2015 . +4.*DTE*VAF(I,J)*(1.0-1.0/ . (1.0-0.5*(IVEG2D(I,J)+IVEG2D(I,J-1)))) ) . /((H(I,J)+ELF(I,J)+H(I,J-1)+ELF(I,J-1))*ARV(I,J)) C CJKL IF(BCTYPE.EQ.'RANDB ') CALL RANDB IF(BCTYPE.EQ.'MIXED ') CALL MIXED(1) CJKL CALL BCOND(2,DTI2,0) C chli print *,'returning from EXTRNL' RETURN END