hhnoqual.f

水文模型的原始代码

      subroutine hhnoqual

!!    ~ ~ ~ PURPOSE ~ ~ ~
!!    this subroutine performs in-stream nutrient calculations. No trans-
!!    formations are calculated

!!    ~ ~ ~ INCOMING VARIABLES ~ ~ ~
!!    name             |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    ai0              |ug chla/mg alg|ratio of chlorophyll-a to algal biomass
!!    algae(:)         |mg alg/L      |algal biomass concentration in reach
!!    ammonian(:)      |mg N/L        |ammonia concentration in reach
!!    chlora(:)        |mg chl-a/L    |chlorophyll-a concentration in reach
!!    disolvp(:)       |mg P/L        |dissolved P concentration in reach
!!    hhvaroute(2,:,:) |m^3 H2O       |water
!!    hhvaroute(4,:,:) |kg N          |organic nitrogen
!!    hhvaroute(5,:,:) |kg P          |organic posphorus
!!    hhvaroute(6,:,:) |kg N          |nitrate
!!    hhvaroute(7,:,:) |kg P          |soluble phosphorus
!!    hhvaroute(13,:,:)|kg            |chlorophyll-a
!!    hhvaroute(14,:,:)|kg N          |ammonium
!!    hhvaroute(15,:,:)|kg N          |nitrite
!!    hhvaroute(16,:,:)|kg            |carbonaceous biological oxygen demand
!!    hhvaroute(17,:,:)|kg O2         |dissolved oxygen
!!    hrchwtr(ii)      |m^3 H2O       |water stored in reach at beginning of day
!!    hrtwtr(:)        |m^3 H2O       |flow out of reach
!!    inum1            |none          |reach number
!!    inum2            |none          |inflow hydrograph storage location number
!!    nitraten(:)      |mg N/L        |nitrate concentration in reach
!!    nitriten(:)      |mg N/L        |nitrite concentration in reach
!!    organicn(:)      |mg N/L        |organic nitrogen concentration in reach
!!    organicp(:)      |mg P/L        |organic phosphorus concentration in reach
!!    rch_cbod(:)      |mg O2/L       |carbonaceous biochemical oxygen demand in
!!                                    |reach 
!!    rch_dox(:)       |mg O2/L       |dissolved oxygen concentration in reach
!!    rnum1            |none          |fraction of overland flow
!!    tmpav(:)         |deg C         |average air temperature on current day
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ OUTGOING VARIABLES ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    algae(:)    |mg alg/L      |algal biomass concentration in reach
!!    ammonian(:) |mg N/L        |ammonia concentration in reach
!!    chlora(:)   |mg chl-a/L    |chlorophyll-a concentration in reach
!!    disolvp(:)  |mg P/L        |dissolved phosphorus concentration in reach
!!    nitraten(:) |mg N/L        |nitrate concentration in reach
!!    nitriten(:) |mg N/L        |nitrite concentration in reach
!!    organicn(:) |mg N/L        |organic nitrogen concentration in reach
!!    organicp(:) |mg P/L        |organic phosphorus concentration in reach
!!    rch_cbod(:) |mg O2/L       |carbonaceous biochemical oxygen demand in
!!                               |reach
!!    rch_dox(:)  |mg O2/L       |dissolved oxygen concentration in reach
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ LOCAL DEFINITIONS ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    algcon      |mg alg/L      |initial algal biomass concentration in reach
!!    algin       |mg alg/L      |algal biomass concentration in inflow
!!    ammoin      |mg N/L        |ammonium N concentration in inflow
!!    cbodcon     |mg/L          |initial carbonaceous biological oxygen demand
!!                               |concentration in reach
!!    cbodin      |mg/L          |carbonaceous biological oxygen demand 
!!                               |concentration in inflow
!!    chlin       |mg chl-a/L    |chlorophyll-a concentration in inflow
!!    disoxin     |mg O2/L       |dissolved oxygen concentration in inflow
!!    dispin      |mg P/L        |soluble P concentration in inflow
!!    jrch        |none          |reach number
!!    nh3con      |mg N/L        |initial ammonia concentration in reach
!!    nitratin    |mg N/L        |nitrate concentration in inflow
!!    nitritin    |mg N/L        |nitrite concentration in inflow
!!    no2con      |mg N/L        |initial nitrite concentration in reach
!!    no3con      |mg N/L        |initial nitrate concentration in reach
!!    o2con       |mg O2/L       |initial dissolved oxygen concentration in 
!!                               |reach
!!    orgncon     |mg N/L        |initial organic N concentration in reach
!!    orgnin      |mg N/L        |organic N concentration in inflow
!!    orgpcon     |mg P/L        |initial organic P concentration in reach
!!    orgpin      |mg P/L        |organic P concentration in inflow
!!    solpcon     |mg P/L        |initial soluble P concentration in reach
!!    wtmp        |deg C         |temperature of water in reach
!!    wtrin       |m^3 H2O       |water flowing into reach on day
!!    wtrtot      |m^3 H2O       |inflow + storage water
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ SUBROUTINES/FUNCTIONS CALLED ~ ~ ~
!!    Intrinsic: Log, Exp, Min
!!    SWAT: Theta

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

      use parm

      integer :: jrch, ii
      real :: wtrin, chlin, algin, orgnin, ammoin, nitratin, nitritin
      real :: orgpin, dispin, cbodin, disoxin, wtmp
      real :: orgpcon, solpcon, cbodcon, o2con, wtrtot
      real :: algcon, orgncon, nh3con, no2con, no3con

      jrch = 0
      jrch = inum1

!! hourly loop
      do ii = 1, 24
       !! initialize water flowing into reach
       wtrin = 0.
       wtrin = hhvaroute(2,inum2,ii) * (1. - rnum1)

       if (hrtwtr(ii) / 3600. > 0.01) then
!! concentrations
         !! initialize inflow concentrations
         chlin = 0.
         algin = 0.
         orgnin = 0.
         ammoin = 0.
         nitritin = 0.
         nitratin = 0.
         orgpin = 0.
         dispin = 0.
         cbodin = 0.
         disoxin = 0.
         if (wtrin > 0.001) then
         chlin = 1000. * hhvaroute(13,inum2,ii) * (1. - rnum1) / wtrin
         algin = 1000. * chlin / ai0        !! QUAL2E equation III-1
         orgnin = 1000. * hhvaroute(4,inum2,ii) * (1. - rnum1) / wtrin
         ammoin = 1000. * hhvaroute(14,inum2,ii) * (1. - rnum1) / wtrin
         nitritin = 1000. * hhvaroute(15,inum2,ii) * (1. - rnum1) /     &
     &                                                             wtrin
         nitratin = 1000. * hhvaroute(6,inum2,ii) * (1. - rnum1) / wtrin
         orgpin = 1000. * hhvaroute(5,inum2,ii) * (1. - rnum1) / wtrin
         dispin = 1000. * hhvaroute(7,inum2,ii) * (1. - rnum1) / wtrin
         cbodin = 1000. * hhvaroute(16,inum2,ii) * (1. - rnum1) / wtrin
         disoxin= 1000. * hhvaroute(17,inum2,ii) * (1. - rnum1) / wtrin
         end if

         !! initialize concentration of nutrient in reach
         wtrtot = 0.
         algcon = 0.
         orgncon = 0.
         nh3con = 0.
         no2con = 0.
         no3con = 0.
         orgpcon = 0.
         solpcon = 0.
         cbodcon = 0.
         o2con = 0.
         wtrtot = wtrin + hrchwtr(ii)
         if (ii == 1) then
         algcon = (algin * wtrin + algae(jrch) * hrchwtr(ii)) / wtrtot
         orgncon = (orgnin * wtrin + organicn(jrch) * hrchwtr(ii))      &
     &                                                          / wtrtot
         nh3con = (ammoin * wtrin + ammonian(jrch) * hrchwtr(ii))       &
     &                                                          / wtrtot
         no2con = (nitritin * wtrin + nitriten(jrch) * hrchwtr(ii))     &
     &                                                          / wtrtot
         no3con = (nitratin * wtrin + nitraten(jrch) * hrchwtr(ii))     &
     &                                                          / wtrtot
         orgpcon = (orgpin * wtrin + organicp(jrch) * hrchwtr(ii))      &
     &                                                          / wtrtot
         solpcon = (dispin * wtrin + disolvp(jrch) * hrchwtr(ii))       &
     &                                                          / wtrtot
         cbodcon = (cbodin * wtrin + rch_cbod(jrch) * hrchwtr(ii))      &
     &                                                          / wtrtot
         o2con = (disoxin * wtrin + rch_dox(jrch) * hrchwtr(ii))        &
     &                                                          / wtrtot
         else
         algcon = (algin * wtrin + halgae(ii-1) * hrchwtr(ii)) / wtrtot
         orgncon = (orgnin * wtrin + horgn(ii-1) * hrchwtr(ii)) / wtrtot
         nh3con = (ammoin * wtrin + hnh4(ii-1) * hrchwtr(ii)) / wtrtot
         no2con = (nitritin * wtrin + hno2(ii-1) * hrchwtr(ii)) / wtrtot
         no3con = (nitratin * wtrin + hno3(ii-1) * hrchwtr(ii)) / wtrtot
         orgpcon = (orgpin * wtrin + horgp(ii-1) * hrchwtr(ii)) / wtrtot
         solpcon = (dispin * wtrin + hsolp(ii-1) * hrchwtr(ii)) / wtrtot
         cbodcon = (cbodin * wtrin + hbod(ii-1) * hrchwtr(ii)) / wtrtot
         o2con = (disoxin * wtrin + hdisox(ii-1) * hrchwtr(ii)) / wtrtot
         end if

         !! calculate temperature in stream
         !! Stefan and Preudhomme. 1993.  Stream temperature estimation 
         !! from air temperature.  Water Res. Bull. p. 27-45
         !! SWAT manual equation 2.3.13
         wtmp = 0.
         wtmp = 5.0 + 0.75 * tmpav(jrch)
         if (wtmp <= 0.) wtmp = 0.1

!! end initialize concentrations

         !! calculate algal biomass concentration at end of hour
         !! (phytoplanktonic algae)
         halgae(ii) = 0.
         halgae(ii) = algcon
         if (halgae(ii) < 0.) halgae(ii) = 0.

         !! calculate chlorophyll-a concentration at end of hour
         hchla(ii) = 0.
         hchla(ii) = halgae(ii) * ai0 / 1000.

!! oxygen calculations
         !! calculate carbonaceous biological oxygen demand at end of hour
         hbod(ii) = 0.
         hbod(ii) = cbodcon
         if (hbod(ii) < 0.) hbod(ii) = 0.

         !! calculate dissolved oxygen concentration if reach at 
         !! end of hour
         hdisox(ii) = 0.
         hdisox(ii) = o2con 
         if (hdisox(ii) < 0.) hdisox(ii) = 0.
!! end oxygen calculations

!! nitrogen calculations
         !! calculate organic N concentration at end of hour
         horgn(ii) = 0.
         horgn(ii) = orgncon 
         if (horgn(ii) < 0.) horgn(ii) = 0.

        !! calculate ammonia nitrogen concentration at end of hour
        hnh4(ii) = 0.
        hnh4(ii) = nh3con 
        if (hnh4(ii) < 0.) hnh4(ii) = 0.

        !! calculate concentration of nitrite at end of hour
        hno2(ii) = 0.
        hno2(ii) = no2con 
        if (hno2(ii) < 0.) hno2(ii) = 0.

        !! calculate nitrate concentration at end of hour
        hno3(ii) = 0.
        hno3(ii) = no3con 
        if (hno3(ii) < 0.) hno3(ii) = 0.
!! end nitrogen calculations

!! phosphorus calculations
        !! calculate organic phosphorus concentration at end of hour
        horgp(ii) = 0.
        horgp(ii) = orgpcon 
        if (horgp(ii) < 0.) horgp(ii) = 0.

        !! calculate dissolved phosphorus concentration at end of hour
        hsolp(ii) = 0.
        hsolp(ii) = solpcon 
        if (hsolp(ii) < 0.) hsolp(ii) = 0.
!! end phosphorus calculations

      else
        !! all water quality variables set to zero when no flow
        algin = 0.0
        chlin = 0.0
        orgnin = 0.0
        ammoin = 0.0
        nitritin = 0.0
        nitratin = 0.0
        orgpin = 0.0
        dispin = 0.0
        cbodin = 0.0
        disoxin = 0.0
        halgae(ii) = 0.0
        hchla(ii) = 0.0
        horgn(ii) = 0.0
        hnh4(ii) = 0.0
        hno2(ii) = 0.0
        hno3(ii) = 0.0
        horgp(ii) = 0.0
        hsolp(ii) = 0.0
        hbod(ii) = 0.0
        hdisox(ii) = 0.0
      endif
      end do
!! end hourly loop

!! set end of day concentrations
      algae(jrch) = halgae(24)
      chlora(jrch) = hchla(24)
      organicn(jrch) = horgn(24)
      ammonian(jrch) = hnh4(24)
      nitriten(jrch) = hno2(24)
      nitraten(jrch) = hno3(24)
      organicp(jrch) = horgp(24)
      disolvp(jrch) = hsolp(24)
      rch_cbod(jrch) = hbod(24)
      rch_dox(jrch) = hdisox(24)

      return
      end