clicon.f

水文模型的原始代码

          if (slrsim == 2) then
            hru_ra(k) = rabsb
            hru_rmx(k) = rmxbsb
            dayl(k) = daylbsb
            npcp(k) = npcpbsb
            do ii = 1, 24
              frad(k,ii) = fradbsb(ii)
            end do
          end if
          if (wndsim == 2) u10(k) = u10bsb
        else
          if (pcpsim == 2) call pgen(k)
          if (tmpsim == 2) then
            call weatgn(k)
            call tgen(k)
          end if
          if (slrsim == 2) then
            call clgen(k)
            call slrgen(k)
          end if
          if (rhsim == 2) call rhgen(k)
          if (ipet == 1) then
            if (wndsim == 2) call wndgen(k)
          end if
          !! set subbasin generated values
          inum3sprev = 0
          tmxbsb = 0.
          tmnbsb = 0.
          rbsb = 0.
          rhdbsb = 0.
          rabsb = 0.
          rmxbsb = 0.
          daylbsb = 0.
          npcpbsb = 0.
          u10bsb = 0.
          fradbsb = 0.
          inum3sprev = hru_sub(k)
          tmxbsb = tmx(k)
          tmnbsb = tmn(k)
          rbsb = subp(k)
          if (ievent > 0) then
            rhrbsb = 0.
            rstpbsb = 0.
            do l = 1, 24
              rhrbsb(l) = hhsubp(k,l)
            end do
            do l = 1, nstep
              rstpbsb(l) = rainsub(k,l)
            end do
          end if
          rhdbsb = rhd(k)
          rabsb = hru_ra(k)
          rmxbsb = hru_rmx(k)
          daylbsb = dayl(k)
          npcpbsb = npcp(k)
          u10bsb = u10(k)
          do ii = 1, 24
            fradbsb(ii) = frad(k,ii)
          end do
        end if
        tmpav(k) = (tmx(k) + tmn(k)) / 2.
      end do

!! Pinue adjustments!!
c	written by Ann van Griensven
	if (iclb.eq.9) then
      do k = 1, nhru
	call pinuesamp(psamp,isamp(1))	
        subp(k) = subp(k) *                                             &
     &	  (1. + (stprain(k,1)+stprain(k,2)* psamp)/ 100.)
        if (subp(k) < 0.) subp(k) = 0.
        if (nstep > 0) then
          do ii = 1, nstep
            rainsub(k,ii) = rainsub(k,ii) *                             &
     &                 (1. + (stprain(k,1)+stprain(k,2)* psamp) / 100.)
            if (rainsub(k,ii) < 0.) rainsub(k,ii) = 0.
          end do
          do ii = 1, 24
            hhsubp(k,ii) = hhsubp(k,ii) *                               &
     &         (1.+(stprain(k,1)+stprain(k,2)* psamp) / 100.)
            if (hhsubp(k,ii) < 0.) hhsubp(k,ii) = 0.
          end do
        end if
	call pinuesamp(psamp,isamp(2))	
        tmx(k) = tmx(k) + (stptemp(k,1)+stptemp(k,2)* psamp)
        tmn(k) = tmn(k) +  (stptemp(k,1)+stptemp(k,2)* psamp)
        tmpav(k) = tmpav(k) +  (stptemp(k,1)+stptemp(k,2)* psamp)
	call pinuesamp(psamp,isamp(3))	
        hru_ra(k) = hru_ra(k) +  (stprad(k,1)+stprad(k,2)* psamp)
        hru_ra(k) = Max(0.,hru_ra(k))
	call pinuesamp(psamp,isamp(4))	
        rhd(k) = rhd(k) +  (stprhd(k,1)+stprhd(k,2)* psamp)
        rhd(k) = Max(0.01,rhd(k))
        rhd(k) = Min(0.99,rhd(k))
      end do
	end if
!! Climate Change Adjustments !!
      do k = 1, nhru
        subp(k) = subp(k) * (1. + rfinc(hru_sub(k),i_mo) / 100.)
        if (subp(k) < 0.) subp(k) = 0.
        if (nstep > 0) then
          do ii = 1, nstep
            rainsub(k,ii) = rainsub(k,ii) *                             &
     &                              (1. + rfinc(hru_sub(k),i_mo) / 100.)
            if (rainsub(k,ii) < 0.) rainsub(k,ii) = 0.
          end do
          do ii = 1, 24
            hhsubp(k,ii) = hhsubp(k,ii) *                               &
     &                              (1. + rfinc(hru_sub(k),i_mo) / 100.)
            if (hhsubp(k,ii) < 0.) hhsubp(k,ii) = 0.
          end do
        end if
        tmx(k) = tmx(k) + tmpinc(hru_sub(k),i_mo)
        tmn(k) = tmn(k) + tmpinc(hru_sub(k),i_mo)
        tmpav(k) = tmpav(k) + tmpinc(hru_sub(k),i_mo)
        hru_ra(k) = hru_ra(k) + radinc(hru_sub(k),i_mo)
        hru_ra(k) = Max(0.,hru_ra(k))
        rhd(k) = rhd(k) + huminc(hru_sub(k),i_mo)
        rhd(k) = Max(0.01,rhd(k))
        rhd(k) = Min(0.99,rhd(k))
      end do

!! Elevation Adjustments !!
      do k = 1, nhru
        if (elevb(1,hru_sub(k)) > 0. .and.                              &
     &                                elevb_fr(1,hru_sub(k)) > 0.) then
        !! determine temperature and precipitation for individual bands
        ratio = 0.
        do ib = 1, 10
          if (elevb_fr(ib,hru_sub(k)) < 0.) exit
          tdif = 0.
          pdif = 0.
          if (tmpsim == 1) then
            tdif = (elevb(ib,hru_sub(k)) -                              &
     &      Real(elevt(itgage(hru_sub(k))))) * tlaps(hru_sub(k)) / 1000.
          else
            tdif = (elevb(ib,hru_sub(k)) - welev(hru_sub(k)))           &
     &                                       * tlaps(hru_sub(k)) / 1000.
          end if
          if (pcpsim == 1) then
            pdif = (elevb(ib,hru_sub(k)) -                              &
     &      Real(elevp(irgage(hru_sub(k))))) * plaps(hru_sub(k)) / 1000.
          else
            pdif = (elevb(ib,hru_sub(k)) - welev(hru_sub(k)))           &
     &                                       * plaps(hru_sub(k)) / 1000.
          end if
          tavband(ib,k) = tmpav(k) + tdif
          tmxband(ib,k) = tmx(k) + tdif
          tmnband(ib,k) = tmn(k) + tdif
          if (subp(k) > 0.01) then
            pcpband(ib,k) = subp(k) + pdif
            if (pcpband(ib,k) < 0.) pcpband(ib,k) = 0.
          end if
          ratio = ratio + pdif * elevb_fr(ib,hru_sub(k))
        end do
        !! determine fraction change in precipitation for HRU
        if (subp(k) >= 0.01) then
          ratio = ratio / subp(k)
        else
          ratio = 0.
        end if
        !! determine new overall temperature and precipitation values
        !! for HRU
        tmpav(k) = 0.
        tmx(k) = 0.
        tmn(k) = 0.
        subp(k) = 0.
        do ib = 1, 10
          if (elevb_fr(ib,hru_sub(k)) < 0.) exit
          tmpav(k) = tmpav(k) + tavband(ib,k) * elevb_fr(ib,hru_sub(k))
          tmx(k) = tmx(k) + tmxband(ib,k) * elevb_fr(ib,hru_sub(k))
          tmn(k) = tmn(k) + tmnband(ib,k) * elevb_fr(ib,hru_sub(k))
          subp(k) = subp(k) + pcpband(ib,k) * elevb_fr(ib,hru_sub(k))
        end do
        if (nstep > 0) then
          do ii = 1, nstep
            if (rainsub(k,ii) > 0.01) then
              rainsub(k,ii) = rainsub(k,ii) + ratio * rainsub(k,ii)
              if (rainsub(k,ii) < 0.) rainsub(k,ii) = 0.
            end if
          end do
          do ii = 1, 24
            if (hhsubp(k,ii) > 0.01) then
              hhsubp(k,ii) = hhsubp(k,ii) + ratio * hhsubp(k,ii)
              if (hhsubp(k,ii) < 0.) hhsubp(k,ii) = 0.
            end if
          end do
        end if
        end if
      end do

      if (nstep > 0) then
        deallocate (rhrbsb)
        deallocate (rstpbsb)
      end if

      return
 5000 format (i4,i3,f5.1)
 5100 format (7x,f5.1)
      end