routres.f

水文模型的原始代码

      use parm

      integer :: jres, k, ii
      real :: sepmm, resorgpc, ressolpc, sedcon, resorgnc, resno3c
      real :: resno2c, resnh3c

      jres = 0
      jres = inum1

      !! initialize variables for reservoir daily simulation
      call resinit

      if (iyr > iyres(jres) .or.                                        &
     &              (i_mo >= mores(jres) .and. iyr == iyres(jres))) then

        !! Adjust Reservoir Storage for Irrigation Diversions
        call irr_res

        !! perform reservoir water/sediment balance
        call res

        !! perform reservoir nutrient balance
        call resnut

        !! perform reservoir pesticide transformations
        call lakeq

        !! add reservoir seepage to shallow aquifer convert from m^3 to mm 
        if (ressep > 0.) then
          sepmm = 0.
          sepmm = ressep / (da_ha * sub_fr(res_sub(jres)) * 10.)
          do k = 1, nhru
            if (hru_sub(k) == res_sub(jres)) then
              shallst(k) = shallst(k) + sepmm
            end if
          end do
        end if

        !! set values for routing variables
        varoute(1,ihout) = 0.           !!undefined
        varoute(2,ihout) = resflwo
        varoute(3,ihout) = ressedo
        varoute(4,ihout) = resorgno
        varoute(5,ihout) = resorgpo
        varoute(6,ihout) = resno3o
        varoute(7,ihout) = ressolpo
        varoute(8,ihout) = 0.           !!undefined
        varoute(9,ihout) = 0.           !!undefined
        varoute(10,ihout) = 0.          !!undefined
        varoute(11,ihout) = solpesto
        varoute(12,ihout) = sorpesto
        varoute(13,ihout) = reschlao
        varoute(14,ihout) = resnh3o
        varoute(15,ihout) = resno2o
        varoute(16,ihout) = 0.          !!CBOD
        varoute(17,ihout) = 0.          !!dis O2
        varoute(18,ihout) = varoute(18,inum2)  !!persistent bact
        varoute(19,ihout) = varoute(19,inum2)  !!less persistent bact
        varoute(20,ihout) = varoute(20,inum2)  !!conservative metal #1
        varoute(21,ihout) = varoute(21,inum2)  !!conservative metal #2
        varoute(22,ihout) = varoute(22,inum2)  !!conservative metal #3

        if (ievent > 2) then
          do ii = 1, 24
            hhvaroute(1,ihout,ii) = 0.           !!undefined
            hhvaroute(2,ihout,ii) = resflwo / 24.
            hhvaroute(3,ihout,ii) = ressedo / 24.
            hhvaroute(4,ihout,ii) = resorgno / 24.
            hhvaroute(5,ihout,ii) = resorgpo / 24.
            hhvaroute(6,ihout,ii) = resno3o / 24.
            hhvaroute(7,ihout,ii) = ressolpo / 24.
            hhvaroute(8,ihout,ii) = 0.           !!undefined
            hhvaroute(9,ihout,ii) = 0.           !!undefined
            hhvaroute(10,ihout,ii) = 0.          !!undefined
            hhvaroute(11,ihout,ii) = solpesto / 24.
            hhvaroute(12,ihout,ii) = sorpesto / 24.
            hhvaroute(13,ihout,ii) = reschlao / 24.
            hhvaroute(14,ihout,ii) = resnh3o / 24.
            hhvaroute(15,ihout,ii) = resno2o / 24.
            hhvaroute(16,ihout,ii) = 0.          !!CBOD
            hhvaroute(17,ihout,ii) = 0.          !!dis O2
            hhvaroute(18,ihout,ii) = varoute(18,inum2) / 24. !!persistent bact
            hhvaroute(19,ihout,ii) = varoute(19,inum2) / 24. !!less persist bact
            hhvaroute(20,ihout,ii) = varoute(20,inum2) / 24. !!cons metal #1
            hhvaroute(21,ihout,ii) = varoute(21,inum2) / 24. !!cons metal #2
            hhvaroute(22,ihout,ii) = varoute(22,inum2) / 24. !!cons metal #3
          end do
        end if

        !! summarization calculations
        if (curyr > nyskip) then
          !!calculate concentrations
          resorgnc = 0.
          resnh3c = 0.
          resno3c = 0.
          resno2c = 0.
          resorgpc = 0.
          ressolpc = 0.
          sedcon = 0.
          resorgnc = res_orgn(jres) / (res_vol(jres)+.1) * 1000.
          resno3c = res_no3(jres) / (res_vol(jres)+.1) * 1000.
          resno2c = res_no2(jres) / (res_vol(jres)+.1) * 1000.
          resnh3c = res_nh3(jres) / (res_vol(jres)+.1) * 1000.
          resorgpc = res_orgp(jres) / (res_vol(jres)+.1) * 1000.
          ressolpc = res_solp(jres) / (res_vol(jres)+.1) * 1000.
          sedcon = res_sed(jres) * 1.e6

          resoutm(1,jres) = resoutm(1,jres) + resflwi / 86400.
          resoutm(2,jres) = resoutm(2,jres) + resflwo / 86400.
          resoutm(3,jres) = resoutm(3,jres) + ressedi
          resoutm(4,jres) = resoutm(4,jres) + ressedo
          resoutm(5,jres) = resoutm(5,jres) + sedcon
          resoutm(6,jres) = resoutm(6,jres) + respesti
          resoutm(7,jres) = resoutm(7,jres) + reactw
          resoutm(8,jres) = resoutm(8,jres) + volatpst
          resoutm(9,jres) = resoutm(8,jres) + setlpst
          resoutm(10,jres) = resoutm(10,jres) + resuspst
          resoutm(11,jres) = resoutm(11,jres) - difus
          resoutm(12,jres) = resoutm(12,jres) + reactb
          resoutm(13,jres) = resoutm(13,jres) + bury
          resoutm(14,jres) = resoutm(14,jres) + solpesto + sorpesto
          resoutm(15,jres) = resoutm(15,jres) + lkpst_conc(jres)
          resoutm(16,jres) = resoutm(16,jres) + lkspst_conc(jres)
          resoutm(17,jres) = resoutm(17,jres) + resev
          resoutm(18,jres) = resoutm(18,jres) + ressep
          resoutm(19,jres) = resoutm(19,jres) + respcp
          resoutm(20,jres) = resoutm(20,jres) + resflwi
          resoutm(21,jres) = resoutm(21,jres) + resflwo
          resoutm(22,jres) = resoutm(22,jres) + varoute(4,inum2)
          resoutm(23,jres) = resoutm(23,jres) + resorgno
          resoutm(24,jres) = resoutm(24,jres) + varoute(5,inum2)
          resoutm(25,jres) = resoutm(25,jres) + resorgpo
          resoutm(26,jres) = resoutm(26,jres) + varoute(6,inum2)
          resoutm(27,jres) = resoutm(27,jres) + resno3o
          resoutm(28,jres) = resoutm(28,jres) + varoute(15,inum2)
          resoutm(29,jres) = resoutm(29,jres) + resno2o
          resoutm(30,jres) = resoutm(30,jres) + varoute(14,inum2)
          resoutm(31,jres) = resoutm(31,jres) + resnh3o
          resoutm(32,jres) = resoutm(32,jres) + varoute(7,inum2)
          resoutm(33,jres) = resoutm(33,jres) + ressolpo
          resoutm(34,jres) = resoutm(34,jres) + varoute(13,inum2)
          resoutm(35,jres) = resoutm(35,jres) + reschlao
          resoutm(36,jres) = resoutm(36,jres) + resorgpc
          resoutm(37,jres) = resoutm(37,jres) + ressolpc
          resoutm(38,jres) = resoutm(38,jres) + resorgnc
          resoutm(39,jres) = resoutm(39,jres) + resno3c
          resoutm(40,jres) = resoutm(40,jres) + resno2c
          resoutm(41,jres) = resoutm(41,jres) + resnh3c
          wshddayo(11) = wshddayo(11) + ressedc
          wshddayo(34) = wshddayo(34) + resflwi - resflwo
        end if
      else
        !! reservoir has not been constructed yet
        do ii = 1, mvaro
          varoute(ii,ihout) = varoute(ii,inum2)
        end do
      end if

      if (iprint == 1 .and. curyr > nyskip) then
        if (iscen == 1) then
        write (8,5000) jres, iida, res_vol(jres), resflwi / 86400.,     &
     &    (resflwo / 86400.), respcp, resev, ressep, ressedi, ressedo,  &
     &    sedcon, varoute(4,inum2), resorgno, resorgnc,                 &
     &    varoute(5,inum2), resorgpo, resorgpc, varoute(6,inum2),       &
     &    resno3o, resno3c, varoute(15,inum2), resno2o, resno2c,        &
     &    varoute(14,inum2), resnh3o, resnh3c, varoute(7,inum2),        &
     &    ressolpo, ressolpc, varoute(13,inum2), reschlao,              &
     &    res_seci(jres), respesti, reactw, volatpst, setlpst, resuspst,&
     &    difus, reactb, bury, solpesto + sorpesto, lkpst_conc(jres),   &
     &    lkspst_conc(jres)
        else if (isproj == 1) then
        write (22,5000) jres, iida, res_vol(jres), resflwi / 86400.,    &
     &    (resflwo / 86400.), respcp, resev, ressep, ressedi, ressedo,  &
     &    sedcon, varoute(4,inum2), resorgno, resorgnc,                 &
     &    varoute(5,inum2), resorgpo, resorgpc, varoute(6,inum2),       &
     &    resno3o, resno3c, varoute(15,inum2), resno2o, resno2c,        &
     &    varoute(14,inum2), resnh3o, resnh3c, varoute(7,inum2),        &
     &    ressolpo, ressolpc, varoute(13,inum2), reschlao,              &
     &    res_seci(jres), respesti, reactw, volatpst, setlpst, resuspst,&
     &    difus, reactb, bury, solpesto + sorpesto, lkpst_conc(jres),   &
     &    lkspst_conc(jres)
        endif
      endif


      return
 5000 format ('RES   ',i8,1x,i4,41e12.4)
      end