Merge pull request #36 from tugraskan/jeffs

Refactor modules for plant class handling, erosion, sediment, & biomass.

src/basin_module.f90

@@ -111,7 +111,7 @@ module basin_module
         real :: tlaps = 6.5         !! deg C/km     |temperature lapse rate: deg C per km of elevation difference
         real :: nfixmx = 20.0       !! max daily n-fixation (kg/ha)
         real :: decr_min = 0.01     !! minimum daily residue decay
-        real :: rsd_covco = 0.30    !! residue cover factor for computing frac of cover         
+        real :: rsd_covco = 0.75    !! residue cover factor for computing frac of cover         
         real :: urb_init_abst = 1.  !! maximum initial abstraction for urban areas when using Green and Ampt
         real :: petco_pmpt = 100.0  !! PET adjustment (%) for Penman-Montieth and Preiestly-Taylor methods
         real :: uhalpha = 1.0       !! alpha coeff for est unit hydrograph using gamma func

src/basin_prm_default.f90

@@ -37,7 +37,7 @@ subroutine basin_prm_default
          !if (bsn_prm%dorm_hr < 1.e-6) bsn_prm%dorm_hr = -1.        !! time threshold used to define dormant (hrs)
          if (bsn_prm%nfixmx < 1.e-6) bsn_prm%nfixmx = 20.0        !! max daily n-fixation (kg/ha)
          if (bsn_prm%decr_min < 1.e-6) bsn_prm%decr_min = 0.01    !! 
-         if (bsn_prm%rsd_covco < 1.e-6) bsn_prm%rsd_covco = 0.30  !! residue cover factor for computing frac of cover         
+         if (bsn_prm%rsd_covco < 1.e-6) bsn_prm%rsd_covco = 0.75  !! residue cover factor for computing C factor equation         
          if (bsn_prm%urb_init_abst < 1.e-6) bsn_prm%urb_init_abst = 0.  !! PET adjustment (%) for Penman-Montieth and Preiestly-Taylor methods
          if (bsn_prm%petco_pmpt < 0.5 .and. bsn_prm%petco_pmpt > 0.) bsn_prm%petco_pmpt = 0.0   !! reservoir sediment settling coeff
          bsn_prm%petco_pmpt = (100. + bsn_prm%petco_pmpt) / 100.    !! convert to fraction

src/cal_conditions.f90

@@ -7,6 +7,7 @@ subroutine cal_conditions
       use hru_module, only : hru
       use soil_module
       use plant_module
+      use plant_data_module
       use time_module
       use climate_module, only : pcp, tmp
 
@@ -74,6 +75,10 @@ subroutine cal_conditions
                   if (pl_find == "n") cond_met = "n"
                   exit
               end do
+            case ("pl_class")
+              if (cal_upd(ichg_par)%cond(ic)%targc /= pl_class(ielem)) then
+                cond_met = "n"
+              end if
             case ("landuse")    !for hru
               if (cal_upd(ichg_par)%cond(ic)%targc /= hru(ielem)%land_use_mgt_c) then 
                 cond_met = "n"

src/cal_parm_select.f90

@@ -34,7 +34,8 @@ subroutine cal_parm_select (ielem, ly, chg_parm, chg_typ, chg_val, absmin, absma
       use hydrograph_module
       use pesticide_data_module
       use plant_module
-        use gwflow_module
+      use plant_data_module
+      use gwflow_module
 
       implicit none
 
@@ -80,12 +81,16 @@ subroutine cal_parm_select (ielem, ly, chg_parm, chg_typ, chg_val, absmin, absma
         hru(ielem)%lumv%usle_p = chg_par (hru(ielem)%lumv%usle_p,         &
                           chg_typ, chg_val, absmin, absmax)
 
+      case ("usle_c")
+        pldb(ielem)%usle_c = chg_par (pldb(ielem)%usle_c,         &
+                          chg_typ, chg_val, absmin, absmax)
+
       case ("ovn")
         hru(ielem)%luse%ovn = chg_par (hru(ielem)%luse%ovn,               &
                           chg_typ, chg_val, absmin, absmax)
 
       case ("elev")
-        hru(ielem)%topo%elev = chg_par (hru(ielem)%topo%elev,             &
+        hru(ielem)%topo%elev = chg_par (hru(ielem)%topo%elev,             & 
                           chg_typ, chg_val, absmin, absmax)
 
       case ("slope")
@@ -429,8 +434,8 @@ subroutine cal_parm_select (ielem, ly, chg_parm, chg_typ, chg_val, absmin, absma
         bsn_prm%n_updis = chg_par(bsn_prm%n_updis,                      &
                          chg_typ, chg_val, absmin, absmax)
 
-      case ("p_updis")
-        bsn_prm%p_updis = chg_par(bsn_prm%p_updis,                      &
+      case ("rsd_covco")
+        bsn_prm%rsd_covco = chg_par(bsn_prm%rsd_covco,                      &
                          chg_typ, chg_val, absmin, absmax)
 
       case ("dorm_hr")

src/cal_parmchg_read.f90

@@ -102,6 +102,8 @@ subroutine cal_parmchg_read
               cal_upd(i)%num_elem = db_mx%dtbl_res
             case ("plt")
               cal_upd(i)%num_elem = sp_ob%hru
+            case ("pl_class")
+              cal_upd(i)%num_elem = db_mx%plantparm
             case ("lyr")
               cal_upd(i)%num_elem = sp_ob%hru
             case ("sol")

src/ero_cfactor.f90

@@ -103,23 +103,24 @@ subroutine ero_cfactor
 
         !! newer method using residue and biomass cover
         rsd_sumfac = (soil1(j)%rsd(1)%m +1.) / 1000.
-        grnd_covfact = 0.
+        grnd_sumfac = 0.
         can_covfact = 10000.
         do ipl = 1, pcom(j)%npl
           idp = pcom(j)%plcur(ipl)%idplt
           ab_gr_t = pl_mass(j)%ab_gr(ipl)%m / 1000.
-          !grnd_sumfac = grnd_sumfac + ab_gr_t
-          grnd_covfact = grnd_covfact + pldb(idp)%usle_c * ab_gr_t / (ab_gr_t + exp(1.175 - 1.748 * ab_gr_t))
+          grnd_sumfac = grnd_sumfac + ab_gr_t
+          !! grnd_covfact = grnd_covfact + pldb(idp)%usle_c * ab_gr_t / (ab_gr_t + exp(1.175 - 1.748 * ab_gr_t))
           can_covfact = amin1 (can_covfact, pcom(j)%plg(ipl)%cht)
         end do
         !grnd_covfact = grnd_sumfac / (grnd_sumfac + exp(1.175 - 1.748 * grnd_sumfac))
-        rsd_covfact = exp(-0.75 * rsd_sumfac)
+        rsd_covfact = exp(-bsn_prm%rsd_covco * rsd_sumfac)
 
         can_frcov = amin1 (1., pcom(j)%lai_sum / 3.)
         can_covfact = 1. - can_frcov * Exp(-.328 * pcom(j)%cht_mx)
 
-        bio_covfac = 1. - grnd_covfact * exp(-0.1 * can_covfact)
-        c = Max(1.e-10, rsd_covfact * grnd_covfact * bio_covfac)
+        grnd_covfact = exp(-pldb(idp)%usle_c * grnd_sumfac)
+        !! bio_covfac = 1. - grnd_covfact * exp(-0.1 * can_covfact)
+        c = Max(1.e-10, rsd_covfact * grnd_covfact)  ! * can_covfact)
 
         !! erosion output variables
         ero_output(j)%ero_d%c = c

src/mgt_harvbiomass.f90

@@ -40,11 +40,10 @@ subroutine mgt_harvbiomass (jj, iplant, iharvop)
       harv_seed = hi_tot * pl_mass(j)%seed(ipl)
       harv_leaf = hi_tot * pl_mass(j)%leaf(ipl)
       harv_stem = hi_tot * pl_mass(j)%stem(ipl)
-      pl_yield = harv_seed + harv_leaf
-      pl_yield = pl_yield + harv_stem
+      pl_yield = harv_seed + harv_leaf + harv_stem
 
       !! apply pest stress to harvest index - mass lost due to pests - don't add to residue
-      pl_yield = (1. - pcom(j)%plcur(ipl)%pest_stress) * (1. - harveff) * pl_yield
+      pl_yield = (1. - pcom(j)%plcur(ipl)%pest_stress) * pl_yield
       !! add plant carbon for printing
       hrc_d(j)%plant_surf_c = hrc_d(j)%plant_surf_c + pl_yield%c
       hpc_d(j)%harv_abgr_c = hpc_d(j)%harv_abgr_c + pl_yield%c

src/nut_nminrl.f90

@@ -32,7 +32,7 @@ subroutine nut_nminrl
       use septic_data_module
       use basin_module
       use organic_mineral_mass_module
-      use hru_module, only : rsdco_plcom, i_sep, ihru, ipl, isep 
+      use hru_module, only : rsdco_plcom, i_sep, ihru, isep 
       use soil_module
       use plant_module
       use output_landscape_module, only : hnb_d
@@ -43,7 +43,7 @@ subroutine nut_nminrl
       integer :: k = 0      !none          |counter (soil layer)
       integer :: kk = 0     !none          |soil layer used to compute soil water and
                             !              |soil temperature factors
-      integer :: idp = 0
+      !integer :: idp = 0
       real :: rmn1 = 0.     !kg N/ha       |amount of nitrogen moving from fresh organic
                             !              |to nitrate(80%) and active organic(20%)
                             !              |pools in layer
@@ -62,7 +62,7 @@ subroutine nut_nminrl
       real :: cprf = 0.     !              |carbon phosphorus ratio factor
       real :: ca = 0.       !              |
       real :: decr = 0.     !              |
-      real :: rdc = 0.      !              |
+      !real :: rdc = 0.      !              |
       real :: wdn = 0.      !kg N/ha       |amount of nitrogen lost from nitrate pool in
                             !              |layer due to denitrification
       real :: cdg = 0.      !none          |soil temperature factor
@@ -77,6 +77,8 @@ subroutine nut_nminrl
       hnb_d(j)%org_lab_p = 0.
       hnb_d(j)%act_sta_n = 0.
       hnb_d(j)%denit = 0.
+      hnb_d(j)%rsd_nitorg_n = 0.
+      hnb_d(j)%rsd_laborg_p = 0.
 
       !! compute humus mineralization of organic soil pools 
       do k = 1, soil(j)%nly
@@ -147,33 +149,36 @@ subroutine nut_nminrl
             cnrf = 1.
           end if
 
-            if (soil1(j)%rsd(k)%p > 1.e-4) then
-              cpr = soil1(j)%rsd(k)%c / soil1(j)%rsd(k)%p
-              if (cpr > 5000.) cpr = 5000.
-              cprf = Exp(-.693 * (cpr - 200.) / 200.)
-            else
-              cprf = 1.
-            end if
+          if (soil1(j)%rsd(k)%p > 1.e-4) then
+            cpr = soil1(j)%rsd(k)%c / soil1(j)%rsd(k)%p
+            if (cpr > 5000.) cpr = 5000.
+            cprf = Exp(-.693 * (cpr - 200.) / 200.)
+          else
+            cprf = 1.
+          end if
 
-            ca = Min(cnrf, cprf, 1.)
+          ca = Min(cnrf, cprf, 1.)
 
-            !! compute root and incorporated residue decomposition
-            !! all plant residue in soil is mixed - don't track individual plant residue in soil
+          !! compute root and incorporated residue decomposition
+          !! all plant residue in soil is mixed - don't track individual plant residue in soil
 
-            if (pcom(j)%npl > 0) then
-              decr = rsdco_plcom(j) / pcom(j)%npl * ca * csf
-            else
-              decr = 0.05
-            end if
-            decr = Max(bsn_prm%decr_min, decr)
-            decr = Min(decr, 1.)
-            decomp = decr * soil1(j)%rsd(k)
-            soil1(j)%rsd(k) = soil1(j)%rsd(k) - decomp
-            soil1(j)%mn(k)%no3 = soil1(j)%mn(k)%no3 + .8 * decomp%n
-            soil1(j)%hact(k)%n = soil1(j)%hact(k)%n + .2 * decomp%n
-            soil1(j)%mp(k)%lab = soil1(j)%mp(k)%lab + .8 * decomp%p
-            soil1(j)%hsta(k)%p = soil1(j)%hsta(k)%p + .2 * decomp%p
-
+          if (pcom(j)%npl > 0) then
+            decr = rsdco_plcom(j) / pcom(j)%npl * ca * csf
+          else
+            decr = 0.05
+          end if
+          decr = Max(bsn_prm%decr_min, decr)
+          decr = Min(decr, 1.)
+          decomp = decr * soil1(j)%rsd(k)
+          soil1(j)%rsd(k) = soil1(j)%rsd(k) - decomp
+          soil1(j)%mn(k)%no3 = soil1(j)%mn(k)%no3 + .8 * decomp%n
+          soil1(j)%hact(k)%n = soil1(j)%hact(k)%n + .2 * decomp%n
+          soil1(j)%mp(k)%lab = soil1(j)%mp(k)%lab + .8 * decomp%p
+          soil1(j)%hsta(k)%p = soil1(j)%hsta(k)%p + .2 * decomp%p
+
+          hnb_d(j)%rsd_nitorg_n = hnb_d(j)%rsd_nitorg_n + .8 * decomp%n
+          hnb_d(j)%rsd_laborg_p = hnb_d(j)%rsd_laborg_p + .8 * decomp%p
+
           !!  compute denitrification
           wdn = 0.   
           if (i_sep(j) /= k .or. sep(isep)%opt  /= 1) then

src/pl_dormant.f90

@@ -47,7 +47,7 @@ subroutine pl_dormant
           end if
           lai_drop = max (0., lai_drop)
           lai_drop = amin1 (1., lai_drop)
-          leaf_drop%m = lai_drop * pl_mass(j)%leaf(ipl)%m
+          leaf_drop%m = rto * lai_drop * pl_mass(j)%leaf(ipl)%m
           leaf_drop%n = leaf_drop%m * pcom(j)%plm(ipl)%n_fr
           leaf_drop%n = max (0., leaf_drop%n)
           leaf_drop%p = leaf_drop%m * pcom(j)%plm(ipl)%p_fr

src/pl_fert.f90

@@ -19,7 +19,7 @@ subroutine pl_fert (ifrt, frt_kg, fertop)
 
       implicit none 
 
-      real, parameter :: rtof=0.5         !none          |weighting factor used to partition the 
+      real :: rtof             !none          |weighting factor used to partition the 
                                           !              |organic N & P concentration of septic effluent
                                           !              |between the fresh organic and the stable organic pools
       integer :: j = 0                    !none          |hru counter
@@ -36,6 +36,7 @@ subroutine pl_fert (ifrt, frt_kg, fertop)
 
       j = ihru
 
+      rtof = 0.5
       !! calculate c:n ratio for manure applications for SWAT-C
       if (bsn_cc%cswat == 2) then
         org_frt%m = frt_kg

src/plant_data_module.f90

@@ -3,7 +3,7 @@ module plant_data_module
       implicit none
 
       character(len=40), dimension (:), allocatable :: plts_bsn     !none      |plant names simulated in current run
-
+      character(len=25), dimension(:), allocatable :: pl_class      !none      |plant class - row crop, tree, grass, etc
       type plant_db
         character(len=40) :: plantnm = ""  !none              |crop name
         character(len=18) :: typ = ""    !none              |plant category

src/plant_parm_read.f90

@@ -6,9 +6,10 @@ subroutine plant_parm_read
 
       implicit none 
 
-      integer :: ic = 0               !none       |counter
-      character (len=80) :: titldum = ""!           |title of file
-      character (len=80) :: header = "" !           |header of file
+      integer :: ic = 0                   !none       |plant counter
+      character (len=80) :: titldum = ""  !           |title of file
+      character (len=80) :: header = ""   !           |header of file
+      character (len=80) :: plclass = ""  !           |plant class - row crop, close grown, grass, tree, etc
       integer :: eof = 0              !           |end of file
       integer :: imax = 0             !none       |determine max number for array (imax) and total number in file
       integer :: mpl = 0              !           | 
@@ -23,6 +24,7 @@ subroutine plant_parm_read
       if (.not. i_exist .or. in_parmdb%plants_plt == " null") then
         allocate (pldb(0:0))
         allocate (plcp(0:0))
+        allocate (pl_class(0:0))
       else
       do
         open (104,file=in_parmdb%plants_plt)
@@ -37,15 +39,25 @@ subroutine plant_parm_read
           end do
         allocate (pldb(0:imax))
         allocate (plcp(0:imax))
+        allocate (pl_class(0:imax))
 
         rewind (104)
         read (104,*,iostat=eof) titldum
         if (eof < 0) exit
         read (104,*,iostat=eof) header
         if (eof < 0) exit
 
+        read (104,*,iostat=eof) plclass
+        if (plclass /= "nam1") then
+          backspace (105)
+        end if
+
         do ic = 1, imax
-          read (104,*,iostat=eof) pldb(ic)
+          if (plclass /= "nam1") then
+            read (104,*,iostat=eof) pldb(ic)
+          else
+            read (104,*,iostat=eof) pldb(ic), pl_class(ic)
+          end if
           if (eof < 0) exit
           pldb(ic)%mat_yrs = Max (1, pldb(ic)%mat_yrs)
         end do