Merge pull request #14 from crazyzlj/2012rev689

Update to 2012rev689

VERSIONS

@@ -1,3 +1,3 @@
 VERSION_MAJOR 2012
-VERSION_MINOR 688
+VERSION_MINOR 689
 VERSION_PATCH

src/CMakeLists.txt

@@ -35,11 +35,16 @@ endif()
 
 ## Special sources that have fixed length of 72
 set(LEN72_SRCS "grow" "tran")
+## Special sources that have long length of 132 but retained the suffix of "*.f"
+set(LENLONG_SRCS "subbasin")
 foreach(F77FILE ${F77SRCS})
     get_filename_component(CORENAME ${F77FILE} NAME_WE)
     list(FIND LEN72_SRCS ${CORENAME} _FOUND_LEN72)
+    list(FIND LENLONG_SRCS ${CORENAME} _FOUND_LENLONG)
     if(${_FOUND_LEN72} GREATER -1)
         set_source_files_properties(${F77FILE} PROPERTIES COMPILE_FLAGS ${FLAG_L72})
+    elseif(${_FOUND_LENLONG} GREATER -1)
+        set_source_files_properties(${F77FILE} PROPERTIES COMPILE_FLAGS ${FLAG_LLONG})
     else()
         set_source_files_properties(${F77FILE} PROPERTIES COMPILE_FLAGS ${FLAG_L80})
     endif()

src/allocate_parms.f

@@ -60,6 +60,13 @@ subroutine allocate_parms
       mstdo = 113
       motot = 1200    !!! 100 year limit
 
+!! ppet arrays for tropical growth
+      allocate (ppet(mhru))
+      do ihru = 1, mhru
+        allocate (ppet(ihru)%precip(ppet(ihru)%ndays))
+        allocate (ppet(ihru)%pet(ppet(ihru)%ndays))
+      end do
+
 !! Srini 11_1_22
       allocate (tmp_win1(mhru))
       allocate (tmp_win2(mhru))

src/autoirr.f

@@ -83,11 +83,6 @@ subroutine autoirr
       j = 0
       j = ihru
 
-!!!! Srin's irrigation source by each application changes
-      irrsc(j) = irr_sca(j)
-      irrno(j) = irr_noa(j)
-!!!! Srin's irrigation source by each application changes
-
       if ((wstrs_id(j) == 3) .or.
      & (wstrs_id(j) == 1 .and. strsw(j) < auto_wstr(j)) .or.          
      & (wstrs_id(j)==2.and.sol_sumfc(j)-sol_sw(j)>auto_wstr(j))) then
@@ -106,7 +101,7 @@ subroutine autoirr
         select case (irrsc(j))
           case (3)   !! shallow aquifer source
             do k = 1, nhru
-              if (hru_sub(k) == irrno(j)) then
+              if (k == irrno(j)) then
                 cnv = 0.
                 cnv = hru_ha(k) * 10.
                 vmma = vmma + shallst(k) * cnv
@@ -120,7 +115,7 @@ subroutine autoirr
 
           case (4)   !! deep aquifer source
             do k = 1, nhru
-              if (hru_sub(k) == irrno(j)) then
+              if (k == irrno(j)) then
                 cnv = 0.
                 cnv = hru_ha(k) * 10.
                 vmma = vmma + deepst(k) * cnv

src/header.f

@@ -93,7 +93,7 @@ subroutine header
      &         "    WTMPdegc","     Salt1  ","     Salt2  ",
      &         "     Salt3  ","     Salt4  ","     Salt5  ",
      &         "     Salt6  ","     Salt7  ","     Salt8  ",
-     &         "     Salt9  ","     Salt10 ","       SAR  ",
+     &         "     Salt9  ","     Salt10  ","      SAR  ",
      &         "        EC  "/)
 
 !!    numbers printed to VB interface reach output file

src/main.f

@@ -48,7 +48,7 @@ program main
       use parm
       implicit none
       !prog = "SWAT Dec 1   VER 2022/Merge Rev 663/Rev 687"
-      prog = "SWAT Jan 20   VER 2023/Rev 688"
+      prog = "SWAT May 10   VER 2023/Rev 689"
       write (*,1000)
  1000 format(1x,"               SWAT2022               ",/,             
      &          "               Rev. 688               ",/,             

src/modparm.f

@@ -2,6 +2,23 @@ module parm
       integer icalen
       real*8 :: prf_bsn
 
+      type precip_pet_moving_average                  !!for tropical plant growth  
+        integer :: trop = 0                           !!      |1=tropical growth ?moisture driven
+        integer :: peren = 0                          !!      |0=annual crop; 1=perennial
+        integer :: mce = 0                            !!      |my current element -day in the p/pet arrays 
+        integer :: mon_seas = 0                       !!      |0=not monsoon season; 1=in monsoon
+        integer :: ndays_mon = 0                      !!      |number of days in the monsoon period 
+        integer :: curday_mon = 0                     !!      |current day into the monsoon period
+        integer :: ndays = 30                         !!      |number of days for precip/pet moving average
+        real :: trig = 0.5                            !!mm/mm |precip/pet ratio to trigger plant/restart
+        real :: precip_sum                            !!mm    |sum of precip during moving average period
+        real :: pet_sum                               !!mm    |sum of pet during moving average period
+        real, dimension (:), allocatable :: precip    !!mm    |precip dimensioned by ndays 
+        real, dimension (:), allocatable :: pet       !!mm    |pet dimensioned by ndays 
+        real :: rto = 0                               !!      |sum of precip/sum of pet
+      end type precip_pet_moving_average
+      type (precip_pet_moving_average), dimension(:),allocatable :: ppet  !dimensioned by subbasin
+
 !!    srin - co2 (EPA)
       real*8 :: co2_x2, co2_x
 
@@ -11,7 +28,7 @@ module parm
      & tmp_spr2, tmp_fal1, tmp_fal2
 
       real*8 :: wtmp
-
+      real*8, dimension (12) :: pcpmm        
       real*8, dimension (:), allocatable :: alph_e
       real*8, dimension (:), allocatable :: co_p, surlag, cdn, nperco
       real*8, dimension (:), allocatable :: cmn, phoskd, psp, sdnco

src/readfcst.f

@@ -57,7 +57,7 @@ subroutine readfcst
       use parm
 
       character (len=80) :: titldum
-      real*8, dimension (12) :: pcpmm, pcpd
+      real*8, dimension (12) :: pcpd
       integer :: mon, mdays, j, fcstregtot
 
 

src/readsub.f

@@ -454,6 +454,20 @@ subroutine readsub
 !!read in weather generator parameter values
       call readwgn
       plaps(i) = plaps(i) / pcpdays(i)
+        
+      !!initialize ppet for tropical growth
+      !! use average monthly December precip and pet from wgn - assume Jan 1 start
+      mo_ppet = 12
+      ppet(ihru)%precip_sum = 0.
+      ppet(ihru)%pet_sum = 0.
+      do inext = 1, ppet(ihru)%ndays
+        !! set previous precip and pet to December ave daily precip
+        ppet(ihru)%precip(inext) = pcpmm(12) / 30.
+        ppet(ihru)%pet(inext) = pcpmm(12) / 30.
+        ppet(ihru)%precip_sum = ppet(ihru)%precip_sum + ppet(ihru)%precip(inext)
+        ppet(ihru)%pet_sum = ppet(ihru)%pet_sum + ppet(ihru)%pet(inext)
+      end do
+
 !!read in subbasin impoundment parameter values
       call readpnd
 !!read in subbasin water use parameter values

src/readwgn.f

@@ -133,7 +133,7 @@ subroutine readwgn
       character (len=80) :: titldum
       real*8 :: xx, lattan, x1, x2, x3, tav, tmin, tmax, rain_yrs
       real*8 :: summx_t, summn_t, summm_p, sum, rnm2, r6, xlv, pcp
-      real*8, dimension (12) :: rainhhmx, rain_hhsm, pcpmm, pcpd
+      real*8, dimension (12) :: rainhhmx, rain_hhsm, pcpd
       real*8 :: tmpsoil, sffc, rndm1, dl
       integer :: mon, mdays, j, m1, nda, xrnd
 
@@ -346,6 +346,23 @@ subroutine readwgn
         end do
         ffc(ihru) = sffc
         dormhr(ihru) = dl
+        !! set initial precip/pet values for tropical plant growth
+        !! compute potential et with Preistley-Taylor Method
+        mdays = ndays(13) - ndays(12)
+        tav = (tmpmx(12,i) + tmpmn(12,i)) / 2.
+        tk = tav  + 273.
+        alb = .15   !tropical rainforests (0.05-0.15)
+        d = EXP(21.255 - 5304. / tk) * 5304. / tk ** 2
+        gma = d / (d +.68)
+        ho = 23.9 * solarav(12,i) * (1. - alb) / 58.3
+        aph = 1.28
+        pet_dec = aph * ho * gma * 30.
+        ppet(j)%precip = pcpmm(12) / mdays
+        ppet(j)%pet = pet_dec / mdays
+        ppet(j)%precip_sum = pcpmm(12)
+        ppet(j)%pet_sum = pet_dec
+        ppet(j)%rto = ppet(j)%precip_sum / ppet(j)%pet_sum
+
       end do
 
       close (114)

src/rteinit.f

@@ -124,8 +124,8 @@ subroutine rteinit
         else 
          if (bio_e(idplt(idum)) > 1.e-6) then
           subfr_nowtr(isb) = subfr_nowtr(isb) + hru_fr(idum)
-	   end if
-	  end if
+         end if
+       end if
       end do
 
 !!    read in areas associated with .fig record files

src/sched_mgt.f

@@ -61,6 +61,34 @@ subroutine sched_mgt
      &        sol_sumsolp(j)
             end if
 
+          case (18)  !! monsoon planting for tropical growth
+            igro(j) = 1
+            lai_init = mgt5op(nop(j),j)
+            bio_init = mgt6op(nop(j),j)
+            ppet(j)%trop = 1
+            ppet(j)%trig = mgt7op(nop(j),j)
+            bio_targ(j) = mgt8op(nop(j),j) * 1000.
+            cnop = mgt9op(nop(j),j)
+            ppet(j)%ndays_mon = mgt2iop(nop(j),j)
+            ppet(j)%peren = mgt3iop(nop(j),j)
+            if (curyr_mat(j) == 0) igrotree(j) = 1
+
+            idplt(j) = mgt1iop(nop(j),j)
+            ppet(j)%mon_seas = 1
+            ppet(j)%curday_mon = 0
+
+            if (mgt4op(nop(j),j) < 700.) mgt4op(nop(j),j) = 1700.
+!            if (mgt4op(nop(j),j) > 5000.) mgt4op(nop(j),j) = 5000.
+            phu_plt(j) = mgt4op(nop(j),j)
+
+          case (19)  !! reset phenology during monsoon season for tropical perennial growth
+            ppet(j)%mon_seas = 1
+            ppet(j)%trop = 1
+            ppet(j)%peren = mgt3iop(nop(j),j)
+            ppet(j)%curday_mon = 0
+            ppet(j)%trig = mgt7op(nop(j),j)
+            ppet(j)%ndays_mon = mgt2iop(nop(j),j)
+
 
           case (2)  !! irrigation operation
             irr_sc(ihru) = mgt2iop(nop(j),j)     !!NUBZ
@@ -234,7 +262,7 @@ subroutine sched_mgt
           if (imgt ==1) then
             write (143, 1010) subnum(j), hruno(j), iyr, i_mo, 
      *      iida, hru_km(j), "        ",
-     *      "SCHED AUTORR", phubase(j), phuacc(j), sol_sw(j), bio_ms(j),
+     *      "SCHED AUTOIRR", phubase(j), phuacc(j), sol_sw(j), bio_ms(j),
      *      sol_rsd(1,j), sol_sumno3(j),sol_sumsolp(j)
 1010  format (a5,1x,a4,3i6,1x,e10.5,1x,2a15,7f10.2)
           end if

src/stdaa.f

@@ -533,8 +533,7 @@ subroutine stdaa
      & 'AUTOPkh ',t84,'MIXEF',t90,'PRECmm',t97,'SURQGENmm',t109,        
      & 'GWQmm',t118,'ETmm',t125,'SEDth ',t132,'NO3kgh ',t140,           
      & 'ORGNkgh ',t148,'BIOMth',t157,'YLDth',t164,'SURQmm')
- !1900 format (i7,i4,3x,a16,3x,e8.3,17f10.2) 
- 1900 format (i7,i4,3x,a8,3x,e8.3,17f8.2) 
+ 1900 format (i7,i4,3x,a8,3x,e8.3,17f8.2)
  2000 format (///,t17,'AVE MONTHLY BASIN VALUES',/t20,'SNOW',t46,       
      &   'WATER',t66,'SED',/t3,'MON',t11,'RAIN',t20,'FALL',t27,'SURF Q',
      &    t37,'LAT Q',t46,'YIELD',t58,'ET',t64,'YIELD',t75,'PET',/t11,  

src/subbasin.f

@@ -197,6 +197,79 @@ subroutine subbasin
         !! calculate soil temperature for soil layers
         call solt
 
+        !! compute evapotranspiration
+        call etpot
+!        if (pot_vol(j) < 1.e-6) call etact
+        call etact
+
+        !! Update CMI and Precip minus PET 30 day moving sum
+        ppet(j)%mce = ppet(j)%mce + 1
+        if (ppet(j)%mce > ppet(j)%ndays) ppet(j)%mce = 1
+        ppet_mce = ppet(j)%mce
+        !! calculate climatic moisture index - cumulative p/pet
+        !! subtract the 30 day previous and add the current day precip/pet
+        ppet(j)%precip_sum = ppet(j)%precip_sum + precipday - ppet(j)%precip(ppet_mce)
+        ppet(j)%pet_sum = ppet(j)%pet_sum + pet_day - ppet(j)%pet(ppet_mce)
+        ppet(j)%rto = ppet(j)%precip_sum / ppet(j)%pet_sum
+        ppet(j)%precip(ppet_mce) = precipday
+        ppet(j)%pet(ppet_mce) = pet_day
+        if (ppet(j)%trop > 0) then
+        if (ppet(j)%mon_seas > 0) then
+          ppet(j)%curday_mon = ppet(j)%curday_mon + 1
+          !! plant if last day of monsoon and not planted
+          if (ppet(j)%curday_mon == ppet(j)%ndays_mon) then
+            ppet(j)%mon_seas = 0
+            if (ppet(j)%peren == 0) then
+              !! annual planting
+              call plantop 
+              if (imgt == 1) then
+                write (143, 1000) subnum(j), hruno(j), iyr, i_mo, iida,   
+     &          hru_km(j),cpnm(idplt(j))," PLANT", phubase(j), phuacc(j), 
+     &          sol_sw(j),bio_ms(j), sol_rsd(1,j),sol_sumno3(j),
+     &          sol_sumsolp(j)
+              end if
+            else
+              !! perennial phenology reset
+              igro(j) = 1
+              idorm(j) = 0
+              phuacc(j) = 0.
+              if (imgt == 1) then
+                write (143, 1000) subnum(j), hruno(j), iyr, i_mo, iida,   
+     &          hru_km(j),cpnm(idplt(j))," PHENO-RESET", phubase(j), phuacc(j), 
+     &          sol_sw(j),bio_ms(j), sol_rsd(1,j),sol_sumno3(j),
+     &          sol_sumsolp(j)
+              end if
+            end if
+          end if
+        !else
+            !! check if ppet ratio is exceeded
+            if (ppet(j)%rto > ppet(j)%trig) then
+            ppet(j)%mon_seas = 0
+            if (ppet(j)%peren == 0) then
+              !! annual planting
+              call plantop 
+              if (imgt == 1) then
+                write (143, 1000) subnum(j), hruno(j), iyr, i_mo, iida,   
+     &          hru_km(j),cpnm(idplt(j))," PLANT", phubase(j), phuacc(j), 
+     &          sol_sw(j),bio_ms(j), sol_rsd(1,j),sol_sumno3(j),
+     &          sol_sumsolp(j)
+              end if
+            else
+              !! perennial phenology reset
+              igro(j) = 1
+              idorm(j) = 0
+              phuacc(j) = 0.
+              if (imgt == 1) then
+                write (143, 1000) subnum(j), hruno(j), iyr, i_mo, iida,   
+     &          hru_km(j),cpnm(idplt(j))," PHENO-RESET", phubase(j), phuacc(j), 
+     &          sol_sw(j),bio_ms(j), sol_rsd(1,j),sol_sumno3(j),
+     &          sol_sumsolp(j)
+              end if
+            end if
+          end if
+        end if
+        end if
+
         call surface
 
        !! compute effective rainfall (amount that percs into soil)
@@ -207,16 +280,16 @@ subroutine subbasin
         !! perform management operations
         if (yr_skip(j) == 0) call operatn
 
+!!!! Srin's irrigation source by each application changes
+      irrsc(j) = irr_sca(j)
+      irrno(j) = irr_noa(j)
+!!!! Srin's irrigation source by each application changes
+
         if (irrsc(j) > 2) call autoirr       
 
         !! perform soil water routing
         call percmain
 
-        !! compute evapotranspiration
-        call etpot
-!        if (pot_vol(j) < 1.e-6) call etact
-        call etact
-
         !! compute water table depth using climate drivers
         call wattable
 
@@ -251,8 +324,9 @@ subroutine subbasin
         !! compute crop growth
         call plantmod
 
-        !! check for dormancy
-        if (igro(j) == 1) call dormant
+        !! check for dormancy - not for tropical plants
+        if (igro(j) == 1 .and. ppet(j)%trop == 0) call dormant
+
         !! compute actual ET for day in HRU
         etday = ep_day + es_day + canev
 

src/writed.f

@@ -183,6 +183,5 @@ subroutine writed
 5001  format(2i5,500f12.4)
 5100  format(1x,a5,a4,1x,i4,1x,i3,1x,250(e16.4,1x))
 5200  format(i7,i9,i6,i5,1x,e9.4,f12.3,f7.1,f14.3)
-!!6200  format(i5,13f7.2,2f5.2,1x,5f8.2)
 6200  format(i5,15f8.2,1x,4f8.2)
       end