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 VEGFORCE C REZA MARSOOLI, 2015 C CALCULATES THE VEGETATION DRAG & INERTIA FORCES IN MOMENTUM EQS. C VEGETATION STEMS ARE ASSUMED TO BE CYLINDRICAL STANDING ELEMENTS C DRAG TERM=0.5*CD*N*B*U C INERTIA TERM=CM*N*A C CD: VEGETATION DRAG COEFFICIENT (-) C N : VEGETATION DENSITY (STEMS/M^2) C B : STEM DIAMETER (M) C U : FLOW VELOCITY ACTING ON STEMS (M/S) C CM: VEGETATION INERTIA COEFFICIENT (-), =2.0 C A : HORIZONTAL AREA OF STEMS (M^2) INCLUDE 'comdeck' SAVE C !FIND THE ACTUAL HEIGHT OF THE STEM (DUE TO FLEXIBILITY) C !TSUJIMOTO ET AL. (1996) METHOD (CANTILEVER BEAM THEORY) C !EQ. (10) IN ZHAO AND CHEN (2013) C IF SIMULATING FLEXIBLE VEG (E.G. SEAGRASS) OPEN THE LINES IN DO LOOP 5 DO 5 II=1,NUMVEG I=IDXVEG(II) J=IDYVEG(II) HVEGA(I,J)=HVEG(I,J) c !VEGETAL STRESS 'TAUVEG' c TAUVEG=0.0 c DO K=1,KBM1 c TAUVEG=TAUVEG-(DVEG3D(I,J,K)/ART(I,J)/(D(I,J)+SMALLVAL)) c . *(DZ(K)*D(I,J)) c . *SQRT( 0.25*((U(I,J,K)+U(I+1,J,K))**2+ c . (V(I,J,K)+V(I,J+1,K))**2) )*1020. c TAUVEG=AMAX1(1.0E-15,TAUVEG) c ENDDO c TEMPH=EXP(-4.66*NVEG(I,J)*STIFF(I,J) c . /(TAUVEG*HVEG(I,J)**2)) c HVEGA(I,J)=AMIN1( HVEG(I,J) , HVEG(I,J)*(1.-0.89*TEMPH) ) c HVEGA(I,J)=AMAX1(HVEGA(I,J),0.0) 5 CONTINUE C INITIALIZE THE VEGEATTION DRAG AND INERTIA TERMS DO 10 II=1,NUMVEG I=IDXVEG(II) J=IDYVEG(II) DO 10 K=1,KBM1 DVEG3D(I,J,K)=0.0 IVEG3D(I,J,K)=0.0 10 CONTINUE DO 20 II=1,NUMVEG I=IDXVEG(II) J=IDYVEG(II) DVEG2D(I,J)=0.0 IVEG2D(I,J)=0.0 20 CONTINUE C !CALCULATE VEG. DRAG AND INERTIA FORCES DO 30 II=1,NUMVEG I=IDXVEG(II) J=IDYVEG(II) DO 30 K=1,KBM1 C !CALCULATE SOLID VOLUME FRACTION VEGPHI=NVEG(I,J)*0.785*BVEG(I,J)**2 C !FIND VELOCITY INTEGRATED OVER STEM HEIGHT VELAVEG=0.0 VELAVCOUNT=SMALLVAL DO KK=1,KBM1 IF( HVEGA(I,J).GT.(-Z(KB)+Z(KK+1))*D(I,J) ) THEN VELAVCOUNT=VELAVCOUNT+1.0 VELAVEG=VELAVEG+SQRT(0.25*((U(I,J,KK)+U(I+1,J,KK))**2+ . (V(I,J,KK)+V(I,J+1,KK))**2) ) !may not be very accurate ENDIF ENDDO VELAVEG=VELAVEG/VELAVCOUNT C !CALCULATE PLANT REYNOLDS NUMBER REYNP=SMALLVAL+(VELAVEG*BVEG(I,J)/0.000001) C !USE LOCAL VELOCITY INSTEAD OF DEPTH INTEGRATED OVER STEM C REYNP=0.0000000001+SQRT(0.25*((U(I,J,K)+U(I+1,J,K))**2+ C . (V(I,J,K)+V(I,J+1,K))**2) )*BVEG(I,J)/0.000001 c IF(REYNP.LT.100.0) GOTO 30 !for model stability, can be changed to smaller value C !CALCULATE VEGETATION DRAG COEFFICIENT IF(VEGCDM.EQ.'CONSTANT') THEN TEMPCD=CDVEG(I,J) ELSEIF(VEGCDM.EQ.' WHITE') THEN TEMPCD=1.+10./REYNP**0.66667 TEMPCD=AMIN1(TEMPCD,10.) !LIMIT DRAG COEFFICIENT TO 10.0 ELSEIF(VEGCDM.EQ.'TAN&NEPF') THEN !FIND COEFFICIENTS ALPHA0 AND ALPHA1 IF(VEGPHI.LE.0.091) THEN TEMPA0=181.8*VEGPHI+3.30 ELSE TEMPA0=84.0 ENDIF TEMPA1=0.46+3.8*VEGPHI TEMPCD=2.*(TEMPA1+TEMPA0/ABS(REYNP)) TEMPCD=AMIN1(TEMPCD,5.) !LIMIT DRAG COEFFICIENT ELSE WRITE(*,*) 'THE VEGETATION DRAG FORMULA IS NOT SUPPORTED' STOP ENDIF IF( HVEGA(I,J).GT.(-Z(KB)+Z(K+1))*D(I,J) ) THEN C !WETLAND VEGETATION DRAG FORCE DVEG3D(I,J,K)=-0.5*TEMPCD*NVEG(I,J)*BVEG(I,J)* . SQRT( 0.25*((U(I,J,K)+U(I+1,J,K))**2+ . (V(I,J,K)+V(I,J+1,K))**2) )*ART(I,J)*D(I,J) DVEG3D(I,J,K)=DVEG3D(I,J,K)*FSM(I,J)*DHMWAD(I,J) c . *DUMWAD(I,J)*DUMWAD(I+1,J)*DVMWAD(I,J)*DVMWAD(I,J+1) C !WETLAND VEGETATION INERTIA FORCE IVEG3D(I,J,K)=-CDINER*NVEG(I,J)*0.785*BVEG(I,J)**2* . ART(I,J)*D(I,J) IVEG3D(I,J,K)=IVEG3D(I,J,K)*FSM(I,J)*DHMWAD(I,J) c . *DUMWAD(I,J)*DUMWAD(I+1,J)*DVMWAD(I,J)*DVMWAD(I,J+1) ENDIF 30 CONTINUE DO J=2,JM-1 DO I=2,IM-1 IF(FSM(I,J).EQ.0) GOTO 50 IF(DHMWAD(I,J).LT.1.0) THEN !DRY CELL CELLWET=DHMWAD(I+1,J)+DHMWAD(I,J+1)+ . DHMWAD(I-1,J)+DHMWAD(I,J-1) IF(CELLWET.GT.0.0) THEN DO K=1,KBM1 DVEG3D(I,J,K)=(DVEG3D(I+1,J,K)*DHMWAD(I+1,J)+ . DVEG3D(I-1,J,K)*DHMWAD(I-1,J)+ . DVEG3D(I,J+1,K)*DHMWAD(I,J+1)+ . DVEG3D(I,J-1,K)*DHMWAD(I,J-1))/CELLWET IVEG3D(I,J,K)=(IVEG3D(I+1,J,K)*DHMWAD(I+1,J)+ . IVEG3D(I-1,J,K)*DHMWAD(I-1,J)+ . IVEG3D(I,J+1,K)*DHMWAD(I,J+1)+ . IVEG3D(I,J-1,K)*DHMWAD(I,J-1))/CELLWET ENDDO ENDIF ENDIF 50 CONTINUE ENDDO ENDDO c Depth-averaged drag and inertia forces DO 40 II=1,NUMVEG I=IDXVEG(II) J=IDYVEG(II) DO 40 K=1,KBM1 DVEG2D(I,J)=DVEG2D(I,J)+DVEG3D(I,J,K)*DZ(K) IVEG2D(I,J)=IVEG2D(I,J)+IVEG3D(I,J,K)*DZ(K) 40 CONTINUE RETURN END