graze.f

水文模型的原始代码

!!    SWAT: Erfc

!!    ~ ~ ~ ~ ~ ~ END SPECIFICATIONS ~ ~ ~ ~ ~ ~

      use parm

      integer :: j, l, it
      real :: dmi, dmii, gc, gc1, swf, frt_t, xx

      j = 0
      j = ihru

!! if HRU currently not grazed, check to see if it is time
!! to initialize grazing
      if (igrz(j) == 0) then
        if (igraz(nro(j),ngr(j),j) > 0 .and.                            &
     &                              iida >= igraz(nro(j),ngr(j),j)) then
          igrz(j) = 1
          ndeat(j) = 1
        else if (phuacc(j) > phug(nro(j),ngr(j),j)) then
          igrz(j) = 1
          ndeat(j) = 1
        else
          return
        end if
      else
        !! if not first day of grazing increment total days of grazing by one
        ndeat(j) = ndeat(j) + 1
      end if

!! graze only if adequate biomass in HRU
      if (bio_ms(j) > bio_min(j)) then

        !! determine new biomass in HRU
        dmi = 0.
        dmi = bio_ms(j)
        bio_ms(j) = bio_ms(j) - bio_eat(nro(j),ngr(j),j)
        if (bio_ms(j) < bio_min(j)) bio_ms(j) = bio_min(j)

        !! adjust nutrient content of biomass
        plantn(j) = plantn(j) - (dmi - bio_ms(j)) * pltfr_n(j)
        plantp(j) = plantp(j) - (dmi - bio_ms(j)) * pltfr_p(j)
        if (plantn(j) < 0.) plantn(j) = 0.
        if (plantp(j) < 0.) plantp(j) = 0.

        !! remove trampled biomass and add to residue
        dmii = 0.
        dmii = bio_ms(j)
        bio_ms(j) = bio_ms(j) - bio_trmp(nro(j),ngr(j),j)
        if (bio_ms(j) < bio_min(j))  then
          sol_rsd(1,j) = sol_rsd(1,j) + dmii - bio_min(j)
          bio_ms(j) = bio_min(j)
        else
          sol_rsd(1,j) = sol_rsd(1,j) + bio_trmp(nro(j),ngr(j),j)
        endif
        sol_rsd(1,j) = Max(sol_rsd(1,j),0.)
        bio_ms(j) = Max(bio_ms(j),0.)

        !! adjust nutrient content of residue and biomass for
        !! trampling
        plantn(j) = plantn(j) - (dmii - bio_ms(j)) * pltfr_n(j)
        plantp(j) = plantp(j) - (dmii - bio_ms(j)) * pltfr_p(j)
        if (plantn(j) < 0.) plantn(j) = 0.
        if (plantp(j) < 0.) plantp(j) = 0.
        if (dmii - bio_ms(j) > 0.) then
          sol_fon(1,j) = (dmii - bio_ms(j)) * pltfr_n(j) + sol_fon(1,j)
          sol_fop(1,j) = (dmii - bio_ms(j)) * pltfr_p(j) + sol_fop(1,j) 
        end if


        !! apply manure
        it = 0
        it = manure_id(nro(j),ngr(j),j)
        if (manure_kg(nro(j),ngr(j),j) > 0.) then
          l = 1

          sol_no3(l,j) = sol_no3(l,j) + manure_kg(nro(j),ngr(j),j) *    &
     &                 (1. - fnh3n(it)) * fminn(it)
          sol_fon(l,j) = sol_fon(l,j) + manure_kg(nro(j),ngr(j),j) *    &
     &                 forgn(it)
          sol_nh3(l,j) = sol_nh3(l,j) + manure_kg(nro(j),ngr(j),j) *    &
     &                 fnh3n(it) * fminn(it)
          sol_solp(l,j) = sol_solp(l,j) + manure_kg(nro(j),ngr(j),j) *  &
     &                 fminp(it)
          sol_fop(l,j) = sol_fop(l,j) + manure_kg(nro(j),ngr(j),j) *    &
     &                 forgp(it)

!! add bacteria - #cfu/g * t(manure)/ha * 1.e6 g/t * ha/10,000 m^2 = 100.
!! calculate ground cover
          gc = 0.
          gc = (1.99532 - Erfc(1.333 * laiday(j) - 2.)) / 2.1
          if (gc < 0.) gc = 0.

          gc1 = 0.
          gc1 = 1. - gc

          swf = .15

          frt_t = 0.
          frt_t = bact_swf * manure_kg(nro(j),ngr(j),j) / 1000.

          bactp_plt(j) = gc * bactpdb(it) * frt_t * 100. + bactp_plt(j)
          bactlp_plt(j) = gc * bactlpdb(it) * frt_t * 100.+bactlp_plt(j)

          bactpq(j) = gc1 * bactpdb(it)  * frt_t * 100. + bactpq(j)
          bactpq(j) = bactkddb(it) * bactpq(j)

          bactps(j) = gc1 * bactpdb(it) * frt_t * 100. + bactps(j)
          bactps(j) = (1. - bactkddb(it)) * bactps(j)

          bactlpq(j) = gc1 * bactlpdb(it) * frt_t * 100. + bactlpq(j)     
          bactlpq(j) = bactkddb(it) * bactlpq(j)

          bactlps(j) = gc1 * bactlpdb(it) * frt_t * 100. + bactlps(j)
          bactlps(j) = (1. - bactkddb(it)) * bactlps(j)

        endif

        !! reset leaf area index and fraction of growing season
        if (dmi > 1.) then
          laiday(j) = laiday(j) * bio_ms(j) / dmi
          phuacc(j) = phuacc(j) * bio_ms(j) / dmi
        else
          laiday(j) = 0.05
          phuacc(j) = 0.
        endif


        !! summary calculations
        grazn = grazn + manure_kg(nro(j),ngr(j),j) *                    &
     &               (fminn(it) + forgn(it))
        grazp = grazp + manure_kg(nro(j),ngr(j),j) *                    &
     &               (fminp(it) + forgp(it))
        tgrazn(j) = tgrazn(j) + grazn
        tgrazp(j) = tgrazp(j) + grazp

        if (curyr > nyskip) then
          wshd_ftotn = wshd_ftotn + manure_kg(nro(j),ngr(j),j) *        &
     &               hru_dafr(j) * (fminn(it) + forgn(it))
          wshd_forgn = wshd_forgn + manure_kg(nro(j),ngr(j),j) *        &
     &               hru_dafr(j) * forgn(it)
          wshd_fno3 = wshd_fno3 + manure_kg(nro(j),ngr(j),j) *          &
     &               hru_dafr(j) * fminn(it) * (1. - fnh3n(it))
          wshd_fnh3 = wshd_fnh3 + manure_kg(nro(j),ngr(j),j) *          &
     &               hru_dafr(j) * fminn(it) * fnh3n(it)
          wshd_ftotp = wshd_ftotp + manure_kg(nro(j),ngr(j),j) *        &
     &               hru_dafr(j) * (fminp(it) + forgp(it))
          wshd_fminp = wshd_fminp + manure_kg(nro(j),ngr(j),j) *        &
     &               hru_dafr(j) * fminp(it)
          wshd_forgp = wshd_forgp + manure_kg(nro(j),ngr(j),j) *        &
     &               hru_dafr(j) * forgp(it)
          yldkg(nro(j),1,j) = yldkg(nro(j),1,j) + (dmi - bio_ms(j))
          !yldkg(nro(j),icr(j),j) = yldkg(nro(j),icr(j),j) + (dmi - bio_ms(j))
        end if
      end if

!! check to set if grazing period is over
      if (ndeat(j) == grz_days(nro(j),ngr(j),j)) then
        igrz(j) = 0
        ndeat(j) = 0
        ngr(j) = ngr(j) + 1
      end if


      return
      end