pothole.f

水文模型的原始代码

      sumo = 0.

!! conversion factors
      cnv = 0.
      cnv = 10. * hru_ha(j)

!! when water is impounding
      if (imp_trig(nro(j),nrelease(j),ipot(j)) == 0) then

        !! update volume of water in pothole
        pot_vol(ipot(j)) = pot_vol(ipot(j)) + qday * Abs(pot_fr(j)) *   &
     &     cnv
        potflwi(ipot(j)) = potflwi(ipot(j)) + qday * Abs(pot_fr(j)) *   &
     &     cnv
        qday = qday * (1. - Abs(pot_fr(j)))

        !! update sediment in pothole  
        pot_sed(ipot(j)) = pot_sed(ipot(j)) + sedyld(j) * Abs(pot_fr(j)) 
        potsedi(ipot(j)) = potsedi(ipot(j)) + sedyld(j) * Abs(pot_fr(j))
        sedyld(j) = sedyld(j) * (1. - Abs(pot_fr(j)))

        !! update nitrate in pothole
        pot_no3(ipot(j)) = pot_no3(ipot(j)) + surqno3(j) *              &
     &     Abs(pot_fr(j)) * hru_ha(j)
        surqno3(j) = surqno3(j) * (1. - Abs(pot_fr(j)))


        !! adjust water balance if rice or pothole in HRU
        !! i.e. if ipot(j) == j

        if (pot_vol(j) > 0.) then

          !! compute surface area assuming a cone shape (m^2)
          potsa(j) = pi * (3. * pot_vol(j) / (pi * hru_slp(j)))**.66666
          potsa(j) = potsa(j) / 10000.                    !convert to ha
          potsa(j) = Min(potsa(j), hru_ha(j))

          !! compute rainfall on impounded water and add to current volume
          potpcp = subp(j) * potsa(j) * 10.
          pot_vol(j) = pot_vol(j) + potpcp

          !! check if max volume is exceeded
          if (pot_vol(j) > pot_volx(j)) then
            spillo = 0.
            sedloss = 0.
            no3loss = 0.
            spillo = pot_vol(j) - pot_volx(j)
            sumo = sumo + spillo
            pot_vol(j) = pot_volx(j)
            qday = qday + spillo / cnv

            sedloss = pot_sed(j) *  spillo / pot_vol(j)
            sedloss = Min(sedloss, pot_sed(j))
            pot_sed(j) = pot_sed(j) - sedloss
            potsedo = potsedo + sedloss
            sedyld(j) = sedyld(j) + sedloss

            no3loss = pot_no3(j) *  spillo / pot_vol(j)
            no3loss = Min(no3loss, pot_no3(j))
            pot_no3(j) = pot_no3(j) - no3loss
            surqno3(j) = surqno3(j) + no3loss / hru_ha(j)
          endif

          !! limit seepage into soil if profile is near field capacity
          yy = 0.
          if (sol_sw(j) / sol_sumfc(j) < .5) then
            yy = 1.
          elseif (sol_sw(j) / sol_sumfc(j) < 1.) then
            yy = 1. - sol_sw(j) / sol_sumfc(j)
          endif

          !! calculate seepage into soil
          potsep = yy * sol_k(1,j) * potsa(j) * 240.
          potsep = Min(potsep, pot_vol(j))
          pot_vol(j) = pot_vol(j) - potsep
          sol_st(1,j) = sol_st(1,j) + potsep / cnv

          !! redistribute water so that no layer exceeds maximum storage
          do ly = 1, sol_nly(j)
            dg = 0.
            stmax = 0.
            excess = 0.
            if (ly == 1) then
              dg = sol_z(ly,j)
            else
              dg = sol_z(ly,j) - sol_z(ly-1,j)
            end if
            stmax = sol_por(ly,j) * dg
            if (sol_st(ly,j) <= stmax) exit
            excess = stmax - sol_st(ly,j)
            sol_st(ly,j) = stmax
            if (ly + 1 <= sol_nly(j)) then
              sol_st(ly+1,j) = sol_st(ly+1,j) + excess
            end if
          end do

          !! recompute total soil water
          sol_sw(j) = 0.
          do ly = 1, sol_nly(j)
            sol_sw(j) = sol_sw(j) + sol_st(ly,j)
          end do

          !! compute evaporation from water surface
          if (laiday(j) < evlai) then
            potev = (1. - laiday(j) / evlai) * pet_day
            potev = 10. * evpot(j) * potev * potsa(j)          !!units mm => m^3
            potev = Min(potev, pot_vol(j))
            pot_vol(j) = pot_vol(j) - potev
          endif

          !! Check date for release/impounding water on rice fields
          if (iida == irelease(nro(j),nrelease(j),j) .and.              &
     &                        imp_trig(nro(j),nrelease(j),j) == 1) then
            qday = qday + pot_vol(j) / cnv
            sumo = sumo + pot_vol(j)
            pot_vol(j) = 0.
            potsedo = potsedo + pot_sed(j)
            sedyld(j) = sedyld(j) + pot_sed(j)
            pot_sed(j) = 0.
            surqno3(j) = surqno3(j) + pot_no3(j) / hru_ha(j)
            pot_no3(j) = 0.
          else
            tileo = Min(pot_tile(j), pot_vol(j))
            sumo = sumo + tileo
            pot_vol(j) = pot_vol(j) - tileo
            qday = qday + tileo / cnv
            if (pot_vol(j) > 1.) then
              sedloss = 0.
              no3loss = 0.
              sedloss = pot_sed(j) *  tileo / pot_vol(j)
              pot_sed(j) = pot_sed(j) - sedloss
              potsedo = potsedo + sedloss
              sedyld(j) = sedyld(j) + sedloss
              no3loss = pot_no3(j) *  tileo / pot_vol(j)
              pot_no3(j) = pot_no3(j) - no3loss
              surqno3(j) = surqno3(j) + no3loss / hru_ha(j)
            else
              sedyld(j) = 0.
              surqno3(j) = 0.
            endif
          endif

          !! calculate settling in the pothole
          sedsetl = 0.
          sedsetl = (pot_sed(j) - (pot_nsed(j) * pot_vol(j) / 1.e6)) *  &
     &       sed_stl(j)
          sedsetl = Min(sedsetl, pot_sed(j))
          pot_sed(j) = pot_sed(j) - sedsetl
          pot_sed(j) = Max(0., pot_sed(j))
          pot_no3(j) = pot_no3(j) * (1. - pot_no3l(j))
        endif
      endif


        potpcpmm = potpcp / cnv
        potevmm = potev / cnv
        potsepmm = potsep / cnv
        potflwo = sumo / cnv

!! summary calculations
      if (curyr > nyskip) then
        potmm = 0.
        potmm = pot_vol(j) / cnv

        spadyo = spadyo + potflwo * hru_dafr(j)
        spadyev = spadyev + potevmm * hru_dafr(j)
        spadysp = spadysp + potsepmm * hru_dafr(j)
        spadyrfv = spadyrfv + potpcpmm * hru_dafr(j)
      end if

      return
1000  format (1x,i4,2x,9(f8.2,2x))
      end