Update to 2012Rev670

Author: crazyzlj

VERSIONS

@@ -1,3 +1,3 @@
 VERSION_MAJOR 2012
-VERSION_MINOR 664
+VERSION_MINOR 670
 VERSION_PATCH

src/HQDAV.f90

@@ -0,0 +1,55 @@
+      SUBROUTINE HQDAV(A,CBW,QQ,SSS,ZCH,ZX,CHW,FPW,jrch)
+!     adopted from APEX1501 by Jaehak Jeong 2017
+!     THIS SUBPROGRAM COMPUTES FLOW AREA AND DEPTH GIVEN RATE in a reach
+    
+      USE PARM
+      
+      real*8, intent (in out) :: A, ZX, CHW, FPW
+      real*8, intent (in) :: CBW, QQ, SSS, ZCH
+      integer, intent (in) :: jrch
+      
+      real*8:: RFPW, RFPX
+      ZX=.5*ZCH 
+      RFPW = ch_w(2,jrch) * 4. !width of floodplain
+      RFPX = SQRT(ch_s(2,jrch)) * RFPW / ch_n(1,jrch)
+
+      DO IT=1,10
+          IF(QQ>QCAP(jrch))THEN
+              ZX=MAX(ZX,ZCH)
+              ZZ=ZX-ZCH
+              !COMPUTE CH FLOW ABOVE QCAP
+              ACH=CHXA(jrch)+ZZ*ch_w(2,jrch)
+              R=ACH/CHXP(jrch)
+              QCH=ACH*R**.66667*RCHX(jrch)
+              CHW=ch_w(2,jrch)
+              !COMPUTE FP FLOW
+              AFP=ZZ*(RFPW-ch_w(2,jrch))
+              QFP=AFP*ZZ**.66667*RFPX/RFPW
+              Q=QCH+QFP
+              A=ACH+AFP
+              FPW=RFPW
+              NBCF=1
+          ELSE
+              X1=ZX*RCSS(jrch)
+              A=ZX*(CBW+X1)
+              P=CBW+2.*SSS*ZX
+              Q=A**1.66667*RCHX(jrch)/P**.66667
+              CHW=CBW+2.*X1
+              FPW=0.
+              NBCX=1     
+          END IF
+          FU=Q-QQ
+          X6=MAX(1.,QQ)
+          IF(ABS(FU/X6)<.001)EXIT
+          IF(IT==1)THEN
+              DFQ=-.1*ZX
+          ELSE
+              DFDZ=(FU-FU1)/(ZX-ZX1)
+              DFQ=FU/DFDZ
+          END IF
+          FU1=FU
+          ZX1=ZX
+          ZX=ZX-DFQ
+      END DO
+      RETURN
+      END

src/NCsed_leach.f90

@@ -54,9 +54,9 @@ subroutine orgncswat2(iwave)
 
       integer, intent (in) :: iwave
       integer :: j
-      real :: xx, wt1, er, conc
-      real :: sol_mass, QBC, VBC, YBC, YOC, YW, TOT, YEW, X1, PRMT_21, PRMT_44
-      real :: DK,  V, X3, CO, CS, perc_clyr, latc_clyr  
+      real*8 :: xx, wt1, er, conc
+      real*8 :: sol_mass, QBC, VBC, YBC, YOC, YW, TOT, YEW, X1, PRMT_21, PRMT_44
+      real*8 :: DK,  V, X3, CO, CS, perc_clyr, latc_clyr  
       integer :: k
       latc_clyr = 0.
 

src/addh.f

@@ -69,7 +69,6 @@ subroutine addh
 
 !! add hydrograph points (hourly time step)
       if (ievent > 0) then
-   !     do kk = 1, 24
         do kk = 1, nstep  ! modified for urban modeling by J.Jeong 4/15/2008
           if (hhvaroute(2,inum1,kk) + hhvaroute(2,inum2,kk) > 0.1) then
             hhvaroute(1,ihout,kk) = (hhvaroute(1,inum1,kk) *            
@@ -79,12 +78,16 @@ subroutine addh
           end if
         end do
         do ii = 2, mvaro
-   !       do kk = 1, 24
           do kk = 1, nstep  ! modified for urban modeling by J.Jeong 4/15/2008
             hhvaroute(ii,ihout,kk) = hhvaroute(ii,inum1,kk) +           
      *                                            hhvaroute(ii,inum2,kk)
           end do
         end do
+        
+        DO K = 1, nstep
+          QHY(K,ihout,IHX(1))=QHY(K,inum1,IHX(1))+QHY(K,inum2,IHX(1)) !flood routing jaehak 2017
+        END DO
+
       endif
 
       do ii = 29, mvaro

src/albedo.f

@@ -36,7 +36,7 @@ subroutine albedo
       use parm
 
       integer :: j
-      real :: cej, eaj
+      real*8 :: cej, eaj
 
       j = 0
       j = ihru

src/allocate_parms.f

@@ -1173,6 +1173,9 @@ subroutine allocate_parms
       allocate (pnd_cla(mhru))
       allocate (pnd_sag(mhru))
       allocate (pnd_lag(mhru))
+      
+      allocate (twlpnd(mhru))     !!srini pond/wet infiltration to shallow gw storage
+      allocate (twlwet(mhru))     !!srini pond/wet infiltration to shallow gw storage
 
       allocate (pnd_solp(mhru))
       allocate (pnd_solpg(mhru))
@@ -1740,8 +1743,8 @@ subroutine allocate_parms
 !!    LID simulations
 !!    Common variable
 !!    van Genuchten equation's coefficients
-      allocate(lid_vgcl,lid_vgcm,lid_qsurf_total,
-     & lid_farea_sum)
+!      allocate(lid_vgcl,lid_vgcm,lid_qsurf_total,
+!     & lid_farea_sum)
       allocate(lid_cuminf_last(mhru,4),lid_sw_last(mhru,4),
      & interval_last(mhru,4),lid_f_last(mhru,4),lid_cumr_last(mhru,4),
      & lid_str_last(mhru,4),lid_farea(mhru,4),lid_qsurf(mhru,4),
@@ -1881,6 +1884,10 @@ subroutine allocate_parms
        tillage_factor = 0.
       !! By Zhang for C/N cycling
       !! ============================
+       
+       !FLOOD ROUTING
+       allocate(QHY(nstep+1,mhyd,4), NHY(4*msub))
+       allocate(RCHX(msub),RCSS(msub),QCAP(msub),CHXA(msub),CHXP(msub))
       	  
       call zero0
       call zero1

src/alph.f

@@ -70,7 +70,7 @@ subroutine alph(iwave)
 
       integer, intent (in) :: iwave
       integer :: j, k, kk, jj
-      real :: ab, ajp, preceff, rainsum
+      real*8 :: ab, ajp, preceff, rainsum
 
       j = 0
       j = ihru

src/anfert.f

@@ -167,9 +167,9 @@ subroutine anfert
 
       use parm
 
-      real, parameter :: rtoaf = 0.50
+      real*8, parameter :: rtoaf = 0.50
       integer :: j, ly, ifrt
-      real :: tsno3, tpno3, dwfert, xx, targn, tfp
+      real*8 :: tsno3, tpno3, dwfert, xx, targn, tfp
 
       j = 0
       j = ihru

src/apex_day.f

@@ -178,12 +178,12 @@ subroutine apex_day
 
       if (ievent > 0) then
         do ii = 1, nstep
-          hhvaroute(2,ihout,ii) = flodaya(inum1) / real(nstep)
-          hhvaroute(3,ihout,ii) = seddaya(inum1) / real(nstep)
-          hhvaroute(4,ihout,ii) = orgndaya(inum1) / real(nstep)
-          hhvaroute(5,ihout,ii) = orgpdaya(inum1) / real(nstep)
-          hhvaroute(6,ihout,ii) = no3daya(inum1) / real(nstep)
-          hhvaroute(7,ihout,ii) = minpdaya(inum1) / real(nstep)
+          hhvaroute(2,ihout,ii) = flodaya(inum1) / dfloat(nstep)
+          hhvaroute(3,ihout,ii) = seddaya(inum1) / dfloat(nstep)
+          hhvaroute(4,ihout,ii) = orgndaya(inum1) / dfloat(nstep)
+          hhvaroute(5,ihout,ii) = orgpdaya(inum1) / dfloat(nstep)
+          hhvaroute(6,ihout,ii) = no3daya(inum1) / dfloat(nstep)
+          hhvaroute(7,ihout,ii) = minpdaya(inum1) / dfloat(nstep)
 
         end do
       end if

src/apply.f

@@ -64,7 +64,7 @@ subroutine apply
       use parm
 
       integer :: j, kk, k, jj
-      real :: xx, gc
+      real*8 :: xx, gc
 
       j = 0
       j = ihru