virtual.f

水文模型的原始代码

        sub_no3(sb) = sub_no3(sb) / sub_fr(sb)
        sub_latno3(sb) = sub_latno3(sb) / sub_fr(sb)
        sub_gwno3(sb) = sub_gwno3(sb) / sub_fr(sb)
        sub_solp(sb) = sub_solp(sb) / sub_fr(sb)
        sub_gwsolp(sb) = sub_gwsolp(sb) / sub_fr(sb)
        sub_yorgn(sb) = sub_yorgn(sb) / sub_fr(sb)
        sub_yorgp(sb) = sub_yorgp(sb) / sub_fr(sb)
        sub_sedpa(sb) = sub_sedpa(sb) / sub_fr(sb)
        sub_sedps(sb) = sub_sedps(sb) / sub_fr(sb)

        sub_bactp(sb) = sub_bactp(sb) / sub_fr(sb)
        sub_bactlp(sb) = sub_bactlp(sb) / sub_fr(sb)

        if (sub_wyld(sb) > 0.1) then
          sub_wtmp(sb) = sub_wtmp(sb) / sub_wyld(sb)
        else
          sub_wtmp(sb) = 0.0
        end if

        wcklsp(sb) = wcklsp(sb) / sub_fr(sb)
        sub_precip(sb) = sub_precip(sb) / sub_fr(sb)
        sub_surfq(sb) = sub_surfq(sb) / sub_fr(sb)
        sub_tran(sb) = sub_tran(sb) / sub_fr(sb)
        sub_bd(sb) = sub_bd(sb) / sub_fr(sb)
        sub_orgn(sb) = sub_orgn(sb) / sub_fr(sb)
        sub_orgp(sb) = sub_orgp(sb) / sub_fr(sb)
        sub_minp(sb) = sub_minp(sb) / sub_fr(sb)
        sub_minpa(sb) = sub_minpa(sb) / sub_fr(sb)
        sub_minps(sb) = sub_minps(sb) / sub_fr(sb)
       ! do kk = 1, mp
       ! sub_pst(kk,sb) = sub_pst(kk,sb) / sub_fr(sb)
       ! end do

        !!! compute subbasin sediment and organic chemical loadings
        !if (sub_precip(sb) > 1.e-6 .and. (sub_qd(sb)+sub_tran(sb)) > 1.e-6) then
        !  !! calculate subbasin value for al5
        !  call alph(sb)
        !  !! calculate subbasin value for peakr
        !  call pkq(sb)
        !  !! calculate subbasin value for sedyld
        !  call ysed(sb)
        !  if (sedyld(j) > 1.e-6 .and. peakr > 1.e-6) then
        !    !! calculate subbasin value for enratio
        !    call enrsb(sb)
        !    !! in order to calculate sorbed pesticide in runoff for 
        !    !! subbasin, new arrays need to be made to store pesticide
        !    !! info at the subbasin level (pst_sed, sol_pst etc)
        !    !! I did not modify pesty to accomodate subbasin calculations
        !    !! beyond what was already in there
        !    call pesty(sb)
        !    !! total amounts of nutrients in subbasin are calculated after
        !    !! removing nutrients in HRU sediment loop calculations, so
        !    !! subbasin values will always be different than sum of HRU
        !    !! values
        !    !! calculate subbasin value for sedorgn
        !    call orgn(sb)
        !    !! calculate subbasin value for sedorgp and sedminp
        !    call psed(sb)
        !  end if
        !end if
        !!!compare sub_sedy(sb) (HRU average) to sedyld (subbasin calc.)
        !!!compare sub_yorgn(sb) (HRU average) to sedorgn (subbasin calc.)
        !!!compare sub_yorgp(sb) (HRU average) to sedorgp (subbasin calc.)
        !!!compare sub_sedp(sb) (HRU average) to sedminp (subbasin calc.)

        !! assign reach loadings for subbasin
         !! zero out hydrograph storage locations
         do ii = 1, mvaro
           varoute(ii,ihout) = 0.
         end do

         !! set values for different routing variables
         !! storage locations set to zero are not currently used
         varoute(1,ihout) = sub_wtmp(sb)                 !!wtmp
         varoute(2,ihout) = sub_wyld(sb) * sub_ha * 10.  !!qdr
         varoute(3,ihout) = sub_sedy(sb)                 !!sedyld
         varoute(4,ihout) = sub_yorgn(sb) * sub_ha       !!sedorgn
         varoute(5,ihout) = (sub_yorgp(sb) + sub_sedps(sb)) * sub_ha  
                                                  !!sedorgp & sedminps
         varoute(6,ihout) = (sub_no3(sb) + sub_latno3(sb) +             &
     &                  sub_gwno3(sb)) * sub_ha   !!surqno3 & latno3 & no3gw
         varoute(7,ihout) = (sub_solp(sb) + sub_gwsolp(sb) +            &
     &                  sub_sedpa(sb)) * sub_ha   !!surqsolp & minpgw & sedminpa
         varoute(8,ihout) = 0.
         varoute(9,ihout) = 0.
         varoute(10,ihout) = 0.
         varoute(11,ihout) = sub_solpst(sb)              !!sol pst
         varoute(12,ihout) = sub_sorpst(sb)              !!sorb pst
         varoute(13,ihout) = sub_chl(sb)                 !!chl_a
         varoute(14,ihout) = 0.                          !! NH3
         varoute(15,ihout) = 0.                          !! NO2
         varoute(16,ihout) = sub_cbod(sb)                !!cbodu
         varoute(17,ihout) = sub_dox(sb)                 !!doxq & soxy
         if (varoute(2,ihout) > .1) then
          varoute(18,ihout) = sub_bactp(sb) * sub_ha / varoute(2,ihout)
          varoute(19,ihout) = sub_bactlp(sb) * sub_ha / varoute(2,ihout)
         end if
         varoute(20,ihout) = 0.                          !! cmetal #1
         varoute(21,ihout) = 0.                          !! cmetal #2
         varoute(22,ihout) = 0.                          !! cmetal #3
         !! varoute array has space for 33 different routing components

         !! sum variables for hyd.out
         do ii = 1, 6
           shyd(ii,ihout) = shyd(ii,ihout) + varoute(ii+1,ihout)
         end do
         shyd(7,ihout) = shyd(7,ihout) + varoute(11,ihout)

        !! sub-daily calculations
        if (ievent > 2) then
          !! determine water loading other than surface runoff for day
          tothhqd = 0.
          difflw = 0.
          do ii = 1, 24
            tothhqd = tothhqd + sub_hhqd(sb,ii)
          end do
          difflw = sub_wyld(sb) - tothhqd
          if (difflw < 0.) difflw = 0.
         
          !! assume water loadings other than surface runoff (eg groundwater,
          !! etc) are evenly distributed over 24 hr period
          sub_hwyld = 0.
          do ii = 1, 24
            sub_hwyld(ii) = sub_hhqd(sb,ii) + difflw / 24.
          end do

          !! assign reach loadings for subbasin
          !! zero out hydrograph storage locations
          do ii = 1, mvaro
            do kk = 1, 24
              hhvaroute(ii,ihout,kk) = 0.
            end do
          end do

          !! set values for different routing variables
          !! storage locations set to zero are not currently used
          do ii = 1, 24
            ratio = 0.
            if (sub_wyld(sb) > 1.e-3)                                   &
     &                              ratio = sub_hwyld(ii) / sub_wyld(sb)
            if (sub_hwyld(ii) > 0.) then
              hhvaroute(1,ihout,ii) = sub_hhwtmp(sb,ii)            !!wtmp
              hhvaroute(2,ihout,ii) = sub_hwyld(ii) * sub_ha * 10. !!water
              hhvaroute(3,ihout,ii) = varoute(3,ihout) * ratio     !!sedyld
              hhvaroute(4,ihout,ii) = varoute(4,ihout) * ratio     !!sedorgn
              hhvaroute(5,ihout,ii) = varoute(5,ihout) * ratio     !!sedorgp
              hhvaroute(6,ihout,ii) = varoute(6,ihout) * ratio     !!no3
              hhvaroute(7,ihout,ii) = varoute(7,ihout) * ratio     !!minp
              hhvaroute(8,ihout,ii) = 0.
              hhvaroute(9,ihout,ii) = 0.
              hhvaroute(10,ihout,ii) = 0.
              hhvaroute(11,ihout,ii) = varoute(11,ihout) * ratio   !!sol pst
              hhvaroute(12,ihout,ii) = varoute(12,ihout) * ratio   !!sorb pst
              hhvaroute(13,ihout,ii) = varoute(13,ihout) * ratio   !!chl_a
              hhvaroute(14,ihout,ii) = 0.                          !! NH3
              hhvaroute(15,ihout,ii) = 0.                          !! NO2
              hhvaroute(16,ihout,ii) = varoute(16,ihout) * ratio   !!cbodu
              hhvaroute(17,ihout,ii) = varoute(17,ihout) * ratio   !!doxq & soxy
              hhvaroute(18,ihout,ii) = varoute(18,ihout) * ratio   !!bactp
              hhvaroute(19,ihout,ii) = varoute(19,ihout) * ratio   !!bactlp
              hhvaroute(20,ihout,ii) = 0.                          !!cmetal#1
              hhvaroute(21,ihout,ii) = 0.                          !!cmetal#2
              hhvaroute(22,ihout,ii) = 0.                          !!cmetal#3
            end if
          end do
        end if

        !! summary calculations
        if (curyr > nyskip) then
         submono(1,sb) = submono(1,sb) + sub_subp(sb)
         submono(2,sb) = submono(2,sb) + sub_snom(sb)
         submono(3,sb) = submono(3,sb) + sub_qd(sb)
         submono(4,sb) = submono(4,sb) + sub_wyld(sb)
         submono(5,sb) = submono(5,sb) + sub_pet(sb)
         submono(6,sb) = submono(6,sb) + sub_etday(sb)
         submono(7,sb) = submono(7,sb) + sub_sedy(sb) / sub_ha
         submono(8,sb) = submono(8,sb) + sub_yorgn(sb)
         submono(9,sb) = submono(9,sb) + sub_yorgp(sb)
         submono(10,sb) = submono(10,sb) + sub_no3(sb)
         submono(11,sb) = submono(11,sb) + sub_solp(sb)
         submono(12,sb) = submono(12,sb) + sub_gwq(sb)
         submono(13,sb) = submono(13,sb) + sub_sep(sb)
         submono(14,sb) = submono(14,sb) + sub_sedpa(sb) + sub_sedps(sb)
         submono(15,sb) = submono(15,sb) + sub_latq(sb)
         submono(16,sb) = submono(16,sb) + sub_latno3(sb)

          if (iprint == 1) call subday
        end if
      endif

!! initialize irrigation/overland flow variables
!! these variable must be initialized here because irrigation/overland
!! flow is performed in different command loops and the water used in
!! irrigation will not be reported in output files if the 
!! variables are initialized in the regular place
      aird(j) = 0.
      shallirr(j) = 0.
      deepirr(j) = 0.
      ovrlnd(j) = 0.
      potflwi(j) = 0.
      potsedi(j) = 0.
!! end of day calculations
      tmpavp(j) = 0.
      tmpavp(j) = tmpav(j)
      

      return   
      end