wetlan.f

水文模型的原始代码

        xx = wet_fr(j) * hru_ha(j) * hru_dafr(j)
        wet_solp(j) = wet_solp(j) + (surqsolp(j) + sedminpa(j)) * xx
        wet_psed(j) = wet_psed(j) + sedminps(j) * xx
        wet_orgp(j) = wet_orgp(j) + sedorgp(j) * xx
        wet_solpg(j) = wet_solpg(j) + minpgw(j) * xx
        wet_orgn(j) = wet_orgn(j) + sedorgn(j) * xx
        wet_no3(j) = wet_no3(j) + surqno3(j) * xx
        wet_no3s(j) = wet_no3s(j) + latno3(j) * xx
        wet_no3g(j) = wet_no3g(j) + no3gw(j) * xx

        !! remove nutrients entering wetlands from HRU loadings
        xx = 0.
        xx = 1. - wet_fr(j)
        sedorgn(j) = sedorgn(j) * xx
        surqno3(j) = surqno3(j) * xx
        latno3(j) = latno3(j) * xx
        no3gw(j) = no3gw(j) * xx
        sedorgp(j) = sedorgp(j) * xx
        sedminpa(j) = sedminpa(j) * xx
        sedminps(j) = sedminps(j) * xx
        surqsolp(j) = surqsolp(j) * xx
        minpgw(j) = minpgw(j) * xx


        !! new water volume for day 
        wet_vol(j) = wet_vol(j) - wetsep - wetev + wetpcp + wetflwi


        if (wet_vol(j) < 0.001) then

          !! check for volume deficit in wetland
          !! reduce seepage so that the wetland volume is zero
          wetsep = wetsep + wet_vol(j)
          wet_vol(j) = 0.
          !! if seepage is less than the volume deficit, take the 
          !! remainder from evaporation
          if (wetsep < 0.) then
            wetev = wetev + wetsep
            wetsep = 0.
          end if
          wet_sed(j) = 0.
          wet_solp(j) = 0.
          wet_psed(j) = 0.
          wet_orgp(j) = 0.
          wet_solpg(j) = 0.
          wet_orgn(j) = 0.
          wet_no3(j) = 0.
          wet_no3s(j) = 0.
          wet_no3g(j) = 0.
          wet_chla(j) = 0.
          wet_seci(j) = 0.

        else
 
          !! compute new sediment concentration
          wet_sed(j) = (sed * vol + wetsedi) / wet_vol(j)

          !! compute outflow if wetland water volume > 0
          if (wet_vol(j) <= wet_nvol(j)) then
            wetflwo = 0.
          else
            if (wet_vol(j) <= wet_mxvol(j)) then
              wetflwo = (wet_vol(j) - wet_nvol(j)) / 10.
              wet_vol(j) = wet_vol(j) - wetflwo
            else
              wetflwo = wet_vol(j) - wet_mxvol(j)
              wet_vol(j) = wet_mxvol(j)
            end if
          end if
          qdr(j) = qdr(j) + wetflwo / cnv

          !! compute sediment settling
          if (wet_sed(j) > wet_nsed(j)) then
            wet_sed(j) = (wet_sed(j) - wet_nsed(j)) * sed_stl(j) +      &
     &                                                       wet_nsed(j)
          end if

          !! compute sediment leaving wetland
          wetsedo = wet_sed(j) * wetflwo
          sedyld(j) = sedyld(j) + wetsedo

          !! net change in amount of sediment in wetland for day
          wetsedc = vol * sed + wetsedi - wetsedo - wet_sed(j) *        &
     &                                                        wet_vol(j)
          !! determine settling rate for nutrients
          !! part of equation 29.1.3 in SWAT manual
          if (i_mo >= ipnd1(j) .and. i_mo <= ipnd2(j)) then
            iseas = 1
          else
            iseas = 2
          endif
          phosk = 0.
          nitrok = 0.
          phosk = psetlw(iseas,j) * wetsa * 10000. / wet_vol(j)
          phosk = Min(phosk, 1.)
          nitrok = nsetlw(iseas,j) * wetsa * 10000. / wet_vol(j)
          nitrok = Min(nitrok, 1.)

          !! remove nutrients by settling
          !! other part of equation 29.1.3 in SWAT manual
          wet_solp(j) = wet_solp(j) * (1. - phosk)
          wet_psed(j) = wet_psed(j) * (1. - phosk)
          wet_orgp(j) = wet_orgp(j) * (1. - phosk)
          wet_solpg(j) = wet_solpg(j) * (1. - phosk)
          wet_orgn(j) = wet_orgn(j) * (1. - nitrok)
          wet_no3(j) = wet_no3(j) * (1. - nitrok)
          wet_no3s(j) = wet_no3s(j) * (1. - nitrok)
          wet_no3g(j) = wet_no3g(j) * (1. - nitrok)

          if (wet_vol(j) < 1.e-6) wet_vol(j) = 0.0
          if (wet_orgn(j) < 1.e-6) wet_orgn(j) = 0.0
          if (wet_no3(j) < 1.e-6) wet_no3(j) = 0.0
          if (wet_no3s(j) < 1.e-6) wet_no3s(j) = 0.0
          if (wet_no3g(j) < 1.e-6) wet_no3g(j) = 0.0
          if (wet_orgp(j) < 1.e-6) wet_orgp(j) = 0.0
          if (wet_psed(j) < 1.e-6) wet_psed(j) = 0.0
          if (wet_solp(j) < 1.e-6) wet_solp(j) = 0.0
          if (wet_solpg(j) < 1.e-6) wet_solpg(j) = 0.0

          tpco = 0.
          tpco = 1.e+6 * (wet_solp(j) + wet_orgp(j) + wet_psed(j) +     &
     &                            wet_solpg(j)) / (wet_vol(j) + wetflwo)
          chlaco = 0.
          wet_chla(j) = 0.
          wet_seci(j) = 0.
          if (tpco > 1.e-4) then
            !! equation 29.1.6 in SWAT manual
            chlaco = chlaw(j) * 0.551 * (tpco**0.76)
            wet_chla(j) = chlaco * (wet_vol(j) + wetflwo) * 1.e-6
          endif
          if (chlaco > 1.e-4) then
            !! equation 29.1.8 in SWAT manual
            wet_seci(j) = secciw(j) * 6.35 * (chlaco**(-0.473))
          endif

          !! compute nutrients leaving wetland
          yy = 0.
          yy = wetflwo / (wet_vol(j) + wetflwo)
          sedorgn(j) = sedorgn(j) + wet_orgn(j) * yy / hru_ha(j)
          surqno3(j) = surqno3(j) + wet_no3(j) * yy / hru_ha(j)
          latno3(j) = latno3(j) + wet_no3s(j) * yy / hru_ha(j)
          no3gw(j) = no3gw(j) + wet_no3g(j) * yy / hru_ha(j)
          sedorgp(j) = sedorgp(j) + wet_orgp(j) * yy / hru_ha(j)
          sedminps(j) = sedminps(j) + wet_psed(j) * yy / hru_ha(j)
          surqsolp(j) = surqsolp(j) + wet_solp(j) * yy / hru_ha(j)
          minpgw(j) = minpgw(j) + wet_solpg(j) * yy / hru_ha(j)

          !!update nutrient pools in wetlands
          wet_orgn(j) = wet_orgn(j) * (1. - yy)
          wet_no3(j) = wet_no3(j) * (1. - yy)
          wet_no3s(j) = wet_no3s(j) * (1. - yy)
          wet_no3g(j) = wet_no3g(j) * (1. - yy)
          wet_orgp(j) = wet_orgp(j) * (1. - yy)
          wet_psed(j) = wet_psed(j) * (1. - yy)
          wet_solp(j) = wet_solp(j) * (1. - yy)
          wet_solpg(j) = wet_solpg(j) * (1. - yy)
          wet_chla(j) = wet_chla(j) * (1. - yy)


        end if

        !! add impoundment seepage to shallow aquifer 
        shallst(j) = shallst(j) + wetsep / cnv

        !! compute seepage depth for HRU water balance
        twlwet = wetsep / cnv

      end if

      if (qdr(j) < 0.) qdr(j) = 0.
      if (sedyld(j) < 0.) sedyld(j) = 0.

      return
      end