virtual.f

水文模型的原始代码

!!    sub_sw(:)   |mm H2O        |amount of water in soil in subbasin
!!    sub_tran(:) |mm H2O        |transmission losses for day in subbasin
!!    sub_wyld(:) |mm H2O        |water yield on day in subbasin
!!    sub_yorgn(:)|kg N/ha       |organic N in surface runoff on day in subbasin
!!    sub_yorgp(:)|kg P/ha       |organic P in surface runoff on day in subbasin
!!    varoute(1,:) |deg C        |temperature
!!    varoute(2,:) |m^3 H2O      |water
!!    varoute(3,:) |metric tons  |sediment or suspended solid load
!!    varoute(4,:) |kg N         |organic nitrogen
!!    varoute(5,:) |kg P         |organic phosphorus
!!    varoute(6,:) |kg N         |nitrate
!!    varoute(7,:) |kg P         |mineral phosphorus
!!    varoute(11,:)|mg pst       |pesticide in solution
!!    varoute(12,:)|mg pst       |pesticide sorbed to sediment
!!    varoute(13,:)|kg           |chlorophyll-a
!!    varoute(16,:)|kg           |carbonaceous biological oxygen demand
!!    varoute(17,:)|kg           |dissolved oxygen
!!    varoute(18,:)|# cfu/100ml  |persistent bacteria
!!    varoute(19,:)|# cfu/100ml  |less persistent bacteria
!!    wcklsp(:)   |
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ LOCAL DEFINITIONS ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    cnv         |none          |conversion factor (mm/ha => m^3)
!!    difflw      |mm H2O        |difference in total loading and surface runoff
!!                               |loading
!!    ii          |none          |counter
!!    j           |none          |HRU number
!!    kk          |none          |counter
!!    sb          |none          |subbasin number
!!    sub_ha      |ha            |area of subbasin in hectares
!!    sub_hwyld(:)|mm H2O        |water yield from subbasin during hour
!!    tothhqd     |mm H2O        |sum of hourly surface runoff for day
!!    wtmp        |deg C         |temperature of water
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ SUBROUTINES/FUNCTIONS CALLED ~ ~ ~
!!    SWAT: hruday, impndday, subday
!!    SWAT: alph, pkq, ysed, enrsb, pesty, orgn, psed
!!    SWAT: Tair

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

      use parm

      integer :: j, sb, kk, ii
      real :: cnv, sub_ha, wtmp, tothhqd, difflw
      real :: sub_hwyld(24)

      j = 0
      sb = 0
      j = ihru
      sb = inum1

      cnv = 0.
      cnv = hru_ha(j) * 10.

!! write daily HRU output
      if (iprint == 1 .and. curyr > nyskip) call hruday
      if (iprint == 1 .and. curyr > nyskip) call impndday

!! sum HRU results for subbasin 
      if (sb > 0) then
      !! subbasin averages: water
        sub_subp(sb) = sub_subp(sb) + subp(j) * hru_dafr(j)
        sub_snom(sb) = sub_snom(sb) + snomlt * hru_dafr(j)
        sub_pet(sb) = sub_pet(sb) + pet_day * hru_dafr(j)
        sub_etday(sb) = sub_etday(sb) + etday * hru_dafr(j)
        sub_sumfc(sb) = sub_sumfc(sb) + sol_sumfc(j) * hru_dafr(j)
        sub_sw(sb) = sub_sw(sb) + sol_sw(j) * hru_dafr(j)
        sub_sep(sb) = sub_sep(sb) + sepbtm(j) * hru_dafr(j)
        sub_qd(sb) = sub_qd(sb) + qday * hru_dafr(j)
        sub_gwq(sb) = sub_gwq(sb) + gw_q(j) * hru_dafr(j)
        sub_wyld(sb) = sub_wyld(sb) + qdr(j) * hru_dafr(j)
        sub_latq(sb) = sub_latq(sb) + latq(j) * hru_dafr(j)

      !! subbasin averages: hourly water
        if (ievent > 2) then
          do ii = 1, 24
            sub_hhqd(sb,ii) = sub_hhqd(sb,ii) + hhqday(ii) * hru_dafr(j)
            !! equation 2.3.13 in SWAT manual
            sub_hhwtmp(sb,ii) = sub_hhwtmp(sb,ii) + hhqday(ii) *        &
     &                hru_fr(j) * (5.0 + 0.75 * Tair(ii,j))
          end do
        end if

      !! subbasin averages: sediment
        sub_sedy(sb) = sub_sedy(sb) + sedyld(j)

        surqno3(j) = amax1(1.e-6,surqno3(j))
        latno3(j) = amax1(1.e-6,latno3(j))
        no3gw(j) = amax1(1.e-6,no3gw(j))
        surqsolp(j) = amax1(1.e-6,surqsolp(j))
        minpgw(j) = amax1(1.e-6,minpgw(j))
        sedorgn(j) = amax1(1.e-6,sedorgn(j))
        sedorgp(j) = amax1(1.e-6,sedorgp(j))
        sedminpa(j) = amax1(1.e-6,sedminpa(j))
        sedminps(j) = amax1(1.e-6,sedminps(j))
        

      !! subbasin averages: nutrients
        if (latno3(j) < 1.e-6) latno3(j) = 0.0
        sub_no3(sb) = sub_no3(sb) + surqno3(j) * hru_dafr(j)
        sub_latno3(sb) = sub_latno3(sb) + latno3(j) * hru_dafr(j)
        sub_gwno3(sb) = sub_gwno3(sb) + no3gw(j) * hru_dafr(j)
        sub_solp(sb) = sub_solp(sb) + surqsolp(j) * hru_dafr(j)
        sub_gwsolp(sb) = sub_gwsolp(sb) + minpgw(j) * hru_dafr(j)
        sub_yorgn(sb) = sub_yorgn(sb) + sedorgn(j) * hru_dafr(j)
        sub_yorgp(sb) = sub_yorgp(sb) + sedorgp(j) * hru_dafr(j)
        sub_sedpa(sb) = sub_sedpa(sb) + sedminpa(j) * hru_dafr(j)
        sub_sedps(sb) = sub_sedps(sb) + sedminps(j) * hru_dafr(j)

      !! subbasin averages: pesticides
        if (irtpest > 0) then
          sub_solpst(sb) = sub_solpst(sb) + hrupstd(irtpest,1,j) +      &
     &                                              hrupstd(irtpest,4,j)
          sub_sorpst(sb) = sub_sorpst(sb) + hrupstd(irtpest,2,j)
        end if

      !! subbasin averages: bacteria
        sub_bactp(sb) = sub_bactp(sb) + (bactrop + bactsedp)            &
     &                                                     * hru_dafr(j)
        sub_bactlp(sb) = sub_bactlp(sb) + (bactrolp + bactsedlp)        &
     &                                                     * hru_dafr(j)

      !! subbasin averages: water quality indicators
        sub_chl(sb) = sub_chl(sb) + chl_a(j) * (qday * qdfr * cnv)      &
     &                                                           * 1.e-6
        sub_cbod(sb) = sub_cbod(sb) + cbodu(j) * (qdr(j) * qdfr * cnv)  &
     &                                                           * 1.e-3
        sub_dox(sb) = sub_dox(sb) + (doxq(j) * (qdr(j) * qdfr * cnv) +  &
     &               soxy * (qdr(j) * (1. - qdfr) * cnv)) * 1.e-3

      !! subbasin averages: water temperature
      !! Stefan and Preudhomme. 1993.  Stream temperature estimation
      !! from air temperature.  Water Res. Bull. p. 27-45
        wtmp = 0.
        wtmp = 5.0 + 0.75 * tmpav(j)
        sub_wtmp(sb) = sub_wtmp(sb) + wtmp * qdr(j) * hru_dafr(j)

      !! subbasin averages used in subbasin sediment calculations
        wcklsp(sb) = wcklsp(sb) + cklsp(j) * hru_dafr(j)
        sub_precip(sb) = sub_precip(sb) + precipday * hru_dafr(j)
        sub_surfq(sb) = sub_surfq(sb) + surfq(j) * hru_dafr(j)
        sub_tran(sb) = sub_tran(sb) + tloss * hru_dafr(j)
        sub_bd(sb) = sub_bd(sb) + sol_bd(1,j) * hru_dafr(j)
        sub_orgn(sb) = sub_orgn(sb) + (sol_orgn(1,j) +                  &
     &                      sol_aorgn(1,j) + sol_fon(1,j)) * hru_dafr(j)
       ! do kk = 1, mp
       ! sub_pst(kk,sb) = sub_pst(kk,sb) + sol_pst(k,j,1) * hru_dafr(j)
       ! end do
      end if
!! end subbasin summarization calculations

!! perform subbasin level operations after processing last HRU in subbasin
      if (iihru == hrutot(sb)) then
        sub_ha = 0.
        sub_ha = da_ha * sub_fr(sb)

        !! calculate subbasin average values for weighted parameters
        if (subfr_nowtr(sb) > 1.e-6) then
            sub_snom(sb) = sub_snom(sb) / subfr_nowtr(sb)
            sub_sumfc(sb) = sub_sumfc(sb) / subfr_nowtr(sb)
            sub_sw(sb) = sub_sw(sb) / subfr_nowtr(sb)
            sub_sep(sb) = sub_sep(sb) / subfr_nowtr(sb)
            sub_qd(sb) = sub_qd(sb) / subfr_nowtr(sb)
            sub_gwq(sb) = sub_gwq(sb) / subfr_nowtr(sb)
            sub_wyld(sb) = sub_wyld(sb) / subfr_nowtr(sb)
            sub_latq(sb) = sub_latq(sb) / subfr_nowtr(sb)
        else
            sub_snom(sb) = 0.0
            sub_sumfc(sb) = 0.0
            sub_sw(sb) = 0.0 
            sub_sep(sb) = 0.0
            sub_qd(sb) = 0.0
            sub_gwq(sb) = 0.0
            sub_wyld(sb) = 0.0
            sub_latq(sb) = 0.0
        end if

        sub_subp(sb) = sub_subp(sb) / sub_fr(sb)
        sub_pet(sb) = sub_pet(sb) / sub_fr(sb)
        sub_etday(sb) = sub_etday(sb) / sub_fr(sb)

        if (ievent > 2) then
          do ii = 1, 24
            sub_hhqd(sb,ii) = sub_hhqd(sb,ii) / sub_fr(sb)
            if (sub_hhqd(sb,ii) > 0.1) then
            sub_hhwtmp(sb,ii) = sub_hhwtmp(sb,ii) / sub_hhqd(sb,ii)
            else
            sub_hhwtmp(sb,ii) = 0.
            end if
          end do
   
        hqd = 0.
        do ii = 1, 24
          do ib = 1, itb(sb)
            hqd(ib+ii-1) = hqd(ib+ii-1) + sub_hhqd(sb,ii) * uh(sb,ib)
          end do
        end do
    
        do ii = 1, 24
          sub_hhqd(sb,ii) = hqd(ii) + hqdsave(ii)
        end do

        do ii = 1, itb(sb)
          hqdsave(ii) = hqdsave(ii+24) + hqd(ii+24)
        end do 
        end if