nminrl.f

水文模型的原始代码

!!    ~ ~ ~ SUBROUTINES/FUNCTIONS CALLED ~ ~ ~
!!    Intrinsic: Max, Exp, Sqrt, Min, Abs

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

      use parm

      integer :: j, k, kk
      real :: rmn1, rmp, xx, csf, rwn, hmn, hmp, r4, cnr, cnrf, cpr
      real :: cprf, ca, decr, rdc, wdn, cdg, sut

      j = 0
      j = ihru


      do k = 1, sol_nly(j)
 
        kk =0 
        if (k == 1) then
          kk = 2
        else
          kk = k
        end if

        !! mineralization can occur only if temp above 0 deg
        if (sol_tmp(kk,j) > 0.) then
          !! compute soil water factor
          sut = 0.
          sut = .1 + .9 * Sqrt(sol_st(kk,j) / sol_fc(kk,j))
          sut = Min(1., sut)
          sut = Max(.05, sut)

          !!compute soil temperature factor
          xx = 0.
          cdg = 0.
          xx = sol_tmp(kk,j)
          cdg = .9 * xx / (xx + Exp(9.93 - .312 * xx)) + .1
          cdg = Max(.1, cdg)

          !! compute combined factor
          xx = 0.
          csf = 0.
          xx = cdg * sut
          if (xx < 0.) xx = 0.
          if (xx > 1.e6) xx = 1.e6
          csf = Sqrt(xx)

          !! compute flow from active to stable pools
          rwn = 0.
          rwn = .1e-4 * (sol_aorgn(k,j) * (1. / nactfr - 1.) -          &
     &                                                    sol_orgn(k,j))
          if (rwn > 0.) then
            rwn = Min(rwn, sol_aorgn(k,j))
          else
            rwn = -(Min(Abs(rwn), sol_orgn(k,j)))
          endif
          sol_orgn(k,j) = Max(1.e-6, sol_orgn(k,j) + rwn)
          sol_aorgn(k,j) = Max(1.e-6, sol_aorgn(k,j) - rwn)

          !! compute humus mineralization on active organic n
          hmn = 0.
          hmn = cmn * csf * sol_aorgn(k,j)
          hmn = Min(hmn, sol_aorgn(k,j))
          !! compute humus mineralization on active organic p
          xx = 0.
          hmp = 0.
          xx = sol_orgn(k,j) + sol_aorgn(k,j)
          if (xx > 1.e-6) then
            hmp = 1.4 * hmn * sol_orgp(k,j) / xx
          else
            hmp = 0.
          end if
          hmp = Min(hmp, sol_orgp(k,j))
          !! move mineralized nutrients between pools
          sol_aorgn(k,j) = Max(1.e-6, sol_aorgn(k,j) - hmn)
          sol_no3(k,j) = sol_no3(k,j) + hmn
          sol_orgp(k,j) = sol_orgp(k,j) - hmp
          sol_solp(k,j) = sol_solp(k,j) + hmp

          !! compute residue decomp and mineralization of 
          !! fresh organic n and p (upper two layers only)
          rmn1 = 0.
          rmp = 0.
          if (k <= 2) then
            r4 = 0.
            r4 = .58 * sol_rsd(k,j)

            if (sol_fon(k,j) + sol_no3(k,j) > 1.e-4) then
              cnr = 0.
              cnr = r4 / (sol_fon(k,j) + sol_no3(k,j))
              if (cnr > 500.) cnr = 500.
              cnrf = 0.
              cnrf = Exp(-.693 * (cnr - 25.) / 25.)
            else
              cnrf = 1.
            end if

            if (sol_fop(k,j) + sol_solp(k,j) > 1.e-4) then
              cpr = 0.
              cpr = r4 / (sol_fop(k,j) + sol_solp(k,j))
              if (cpr > 5000.) cpr = 5000.
              cprf = 0.
              cprf = Exp(-.693 * (cpr - 200.) / 200.)
            else
              cprf = 1.
            end if

            ca = 0.
            decr = 0.
            rdc = 0.
            ca = Min(cnrf, cprf, 1.)
            decr = rsdco_pl(idplt(nro(j),icr(j),j)) * ca * csf
            decr = Max(.01, decr)
            decr = Min(decr, 1.)
            sol_rsd(k,j) = amax1(1.e-6,sol_rsd(k,j))
            rdc = decr * sol_rsd(k,j)
            sol_rsd(k,j) = sol_rsd(k,j) - rdc
            if (sol_rsd(k,j) < 0.) sol_rsd(k,j) = 0.
            rmn1 = decr * sol_fon(k,j)
            sol_fop(k,j) = amax1(1.e-6,sol_fop(k,j))
            rmp = decr * sol_fop(k,j)

            sol_fop(k,j) = sol_fop(k,j) - rmp
            sol_fon(k,j) = amax1(1.e-6,sol_fon(k,j))
            sol_fon(k,j) = sol_fon(k,j) - rmn1
            sol_no3(k,j) = sol_no3(k,j) + .8 * rmn1
            sol_aorgn(k,j) = sol_aorgn(k,j) + .2 * rmn1
            sol_solp(k,j) = sol_solp(k,j) + .8 * rmp
            sol_orgp(k,j) = sol_orgp(k,j) + .2 * rmp
          end if

          !! compute denitrification 
          wdn = 0.
          if (sut >= sdnco) then
            wdn = sol_no3(k,j) * (1. - Exp(-cdn * cdg * sol_cbn(k,j)))
          else
            wdn = 0.
          end if
          sol_no3(k,j) = sol_no3(k,j) - wdn

          !! summary calculations
          if (curyr > nyskip) then
            wshd_hmn = wshd_hmn + hmn * hru_dafr(j)
            wshd_rwn = wshd_rwn + rwn * hru_dafr(j)
            wshd_hmp = wshd_hmp + hmp * hru_dafr(j)
            wshd_rmn = wshd_rmn + rmn1 * hru_dafr(j)
            wshd_rmp = wshd_rmp + rmp * hru_dafr(j)
            wshd_dnit = wshd_dnit + wdn * hru_dafr(j)
            hmntl = hmntl + hmn
            rwntl = rwntl + rwn
            hmptl = hmptl + hmp
            rmn2tl = rmn2tl + rmn1
            rmptl = rmptl + rmp
            wdntl = wdntl + wdn
          end if
        end if
      end do


      return
      end