Merge pull request #35 from tugraskan/jeffs

12/23 - changes to MUSLE C factor and some carbon updates

Author: odav

.vs/cmake.db

@@ -115,7 +115,7 @@ subroutine cal_parmchg_read
             case ("cha")
               cal_upd(i)%num_elem = sp_ob%chan
             case ("swq")
-              cal_upd(i)%num_elem = sp_ob%chan
+              cal_upd(i)%num_elem = db_mx%ch_nut
             case ("res")
               cal_upd(i)%num_elem = sp_ob%res
             case ("sdc")

src/cal_parmchg_read.f90

@@ -292,13 +292,8 @@ subroutine cbn_zhang2
       !!calculate c/n dynamics for each soil layer
       !!===========================================
       do k = 1, soil(j)%nly
-        if (k == 1) then
-          !10 cm / 1000 = 0.01m; 1 ha = 10000 m2; ton/m3; * 1000 --> final unit is kg/ha; rock fraction is considered
-          sol_mass = (10) / 1000.* 10000. * soil(j)%phys(k)%bd * 1000. * (1. - soil(j)%phys(k)%rock / 100.)            
-        else
-          sol_mass = (soil(j)%phys(k)%d - soil(j)%phys(k-1)%d) / 1000. * 10000. * soil(j)%phys(k)%bd * 1000. * &
-                    (1- soil(j)%phys(k)%rock / 100.)
-        end if        
+        !! mm / 1000 * 10000 m2 / ha * ton/m3 * 1000 kg/ha -> kg/ha; rock fraction is considered
+        sol_mass = 10000. * soil(j)%phys(k)%thick * soil(j)%phys(k)%bd * (1 - soil(j)%phys(k)%rock / 100.)
 
         ! if k = 1, then using temperature, soil moisture in layer 2 to calculate decomposition factor
         if (k == 1) then

src/cbn_zhang2.f90

@@ -66,11 +66,10 @@ subroutine ero_cfactor
         end if
       else
         !! new method using residue and biomass cover
-        rsd_sumfac = 0.
         grnd_sumfac = 0.
+        rsd_sumfac = pldb(idp)%rsd_pctcov * (soil1(j)%rsd(1)%m +1.) / 1000.
         do ipl = 1, pcom(j)%npl
           idp = pcom(j)%plcur(ipl)%idplt
-          rsd_sumfac = rsd_sumfac + pldb(idp)%rsd_pctcov * (soil1(j)%rsd(1)%m +1.) / 1000.
           if (pl_mass(j)%ab_gr(ipl)%m > 1.e-6) then
             ab_gr_t = pl_mass(j)%ab_gr(ipl)%m / 1000.
             grnd_sumfac = grnd_sumfac + 100. * pldb(idp)%usle_c / ab_gr_t
@@ -97,25 +96,30 @@ subroutine ero_cfactor
         grnd_covfact = max (0., grnd_covfact)
 
         !! ***jga
-        grnd_covfact = 1.34 + 0.225 * log(pldb(idp)%usle_c)
-        grnd_covfact = amin1 (1., grnd_covfact)
-        grnd_covfact = max (0., grnd_covfact)
+        !grnd_covfact = 1.34 + 0.225 * log(pldb(idp)%usle_c)
+        !grnd_covfact = amin1 (1., grnd_covfact)
+        !grnd_covfact = max (0., grnd_covfact)
         c = Max(1.e-10, rsd_covfact * can_covfact * grnd_covfact)
 
         !! newer method using residue and biomass cover
-        !rsd_sumfac = 0.
-        !grnd_sumfac = 0.
-        !can_covfact = 10000.
-        !do ipl = 1, pcom(j)%npl
-        !  ab_gr_t = pl_mass(j)%ab_gr(ipl)%m / 1000.
-        !  grnd_sumfac = grnd_sumfac + ab_gr_t
-        !  rsd_sumfac = rsd_sumfac + (soil1(j)%rsd(1)%m +1.) / 1000.
-        !  can_covfact = amin1 (can_covfact, pcom(j)%plg(ipl)%cht)
-        !end do
+        rsd_sumfac = (soil1(j)%rsd(1)%m +1.) / 1000.
+        grnd_covfact = 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))
+          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)
-        !bio_covfac = 1. - grnd_covfact * exp(-0.1 * can_covfact)
-        !c = Max(1.e-10, rsd_covfact * bio_covfac)
+        rsd_covfact = exp(-0.75 * 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)
 
         !! erosion output variables
         ero_output(j)%ero_d%c = c

src/ero_cfactor.f90

@@ -413,7 +413,7 @@ subroutine hyd_connect
     do while (idone == 0)
         do i = 1, sp_ob%objs
 
-        if (iord > 1000) then        
+        if (iord > 5000) then        
           open (9002,file="looping.con",recl = 8000)
           write (9002, *) "LOOPING.CON CHECKING INFINITE LOOPS"
           do iob = 1, sp_ob%objs

src/hyd_connect.f90

@@ -164,7 +164,7 @@ subroutine nut_orgnc2
           end if
           sol_thick = 0.
           sol_thick = soil(j)%phys(k)%d-soil(j)%phys(k-1)%d
-          soil1(j)%tot(1)%c = soil1(j)%str(k)%c + soil1(j)%meta(k)%c + soil1(j)%hp(k)%c + soil1(j)%hs(k)%c 
+          soil1(j)%tot(1)%c = soil1(j)%hp(k)%c + soil1(j)%hs(k)%c 
           Y1 = soil1(j)%microb(k)%c + VBC
           VBC=0.
           IF(Y1>=.01)THEN

src/nut_orgnc2.f90

@@ -51,14 +51,14 @@ subroutine res_nutrient (iob)
       phosk = amin1 (phosk, 1.)
       phosk = max (phosk, 0.)
 
-      !! remove nutrients from reservoir by settling
+      !! remove nutrients from reservoir by settling - exclude soluble nutrients
       !! other part of equation 29.1.3 in SWAT manual
-      wbody%solp = wbody%solp * (1. - phosk)
+      !wbody%solp = wbody%solp * (1. - phosk)
       wbody%sedp = wbody%sedp * (1. - phosk)
       wbody%orgn = wbody%orgn * (1. - nitrok)
-      wbody%no3 = wbody%no3 * (1. - nitrok)
-      wbody%nh3 = wbody%nh3 * (1. - nitrok)
-      wbody%no2 = wbody%no2 * (1. - nitrok)
+      !wbody%no3 = wbody%no3 * (1. - nitrok)
+      !wbody%nh3 = wbody%nh3 * (1. - nitrok)
+      !wbody%no2 = wbody%no2 * (1. - nitrok)
 
       !! calculate chlorophyll-a and water clarity
       chlaco = 0.

src/res_nutrient.f90

@@ -49,7 +49,7 @@ subroutine soil_nutcarb_write(out_freq)
 
       !! sum the output for the entire soil profile
       do j = 1, sp_ob%hru
-        iob = sp_ob1%hru + j - 1   !!!!!! added for new output write  !!nbs
+        iob = sp_ob1%hru + j - 1
         soil1(j)%tot_org = soil_org_z
         soil_prof_hact = soil_org_z
         soil_prof_hsta = soil_org_z