watqual2.f

水文模型的原始代码

!       write(104, *) cordo, 'cordo'
	  	if (o2con.le.0.001) o2con=0.001
	if (o2con.gt.30.) o2con=30.
        cordo = 1.0 - Exp(-0.6 * o2con)

        !! modify ammonia and nitrite oxidation rates to account for
        !! low oxygen
        bc1mod = 0.
        bc2mod = 0.
        bc1mod = bc1(jrch) * cordo
        bc2mod = bc2(jrch) * cordo
!! end O2 impact calculations

         
C	tday is the calculation time step = 1 day

         tday = 1.0

!! algal growth
         !! calculate light extinction coefficient 
         !! (algal self shading) QUAL2E equation III-12
         if (ai0 * algcon > 1.e-6) then
           lambda = lambda0 + (lambda1 * ai0 * algcon) + lambda2 *      &
     &                                        (ai0 * algcon) ** (.66667)
         else
           lambda = lambda0
         endif

         !! calculate algal growth limitation factors for nitrogen
         !! and phosphorus QUAL2E equations III-13 & III-14
         fnn = 0.
         fpp = 0.
         fnn = cinn / (cinn + k_n)
         fpp = solpcon / (solpcon + k_p)

         !! calculate daylight average, photosynthetically active,
         !! light intensity QUAL2E equation III-8
         !! Light Averaging Option # 2
         algi = 0.
         if (dayl(hru1(jrch)) > 0.) then
           algi = hru_ra(hru1(jrch)) * tfact / dayl(hru1(jrch))
         else
           algi = 0.00001
         end if

         !! calculate growth attenuation factor for light, based on
         !! daylight average light intensity QUAL2E equation III-7b
         fl_1 = 0.
         fll = 0.
         fl_1 = (1. / (lambda * rchdep)) *                              &
     &        Log((k_l + algi) / (k_l + algi * (Exp(-lambda * rchdep))))
         fll = 0.92 * (dayl(hru1(jrch)) / 24.) * fl_1

         !! calculcate local algal growth rate
         gra = 0.
         select case (igropt)
           case (1)
             !! multiplicative QUAL2E equation III-3a
             gra = mumax * fll * fnn * fpp
           case (2)
             !! limiting nutrient QUAL2E equation III-3b
             gra = mumax * fll * Min(fnn, fpp)
           case (3)
             !! harmonic mean QUAL2E equation III-3c
             if (fnn > 1.e-6 .and. fpp > 1.e-6) then
               gra = mumax * fll * 2. / ((1. / fnn) + (1. / fpp))
             else
               gra = 0.
             endif
         end select

    
         !! calculate algal biomass concentration at end of day
         !! (phytoplanktonic algae)
         !! QUAL2E equation III-2
         dalgae = 0.
	 setl=min(1.,Theta(rs1(jrch),thrs1,wtmp)/ rchdep)
         dalgae = algcon + (Theta(gra,thgra,wtmp) * algcon -            &
     &    Theta(rhoq,thrho,wtmp) * algcon - setl * algcon) * tday
         if (dalgae < 0.00001) algae(jrch) = 0.00001

         !! calculate chlorophyll-a concentration at end of day
         !! QUAL2E equation III-1

         dchla = 0.
         dchla = dalgae* ai0 / 1000.
!! end algal growth 

!! oxygen calculations
         !! calculate carbonaceous biological oxygen demand at end
         !! of day QUAL2E section 3.5 equation III-26
         yy = 0.
         zz = 0.
         yy = Theta(rk1(jrch),thrk1,wtmp) * cbodcon
         zz = Theta(rk3(jrch),thrk3,wtmp) * cbodcon
         dbod = 0.
         dbod = cbodcon - (yy + zz) * tday
         if (dbod < 0.00001) dbod = 0.00001

         !! calculate dissolved oxygen concentration if reach at 
         !! end of day QUAL2E section 3.6 equation III-28
  
      

	   uu = 0.
         vv = 0.
         ww = 0.
         xx = 0.
         yy = 0.
         zz = 0.
	   hh=Theta(rk2(jrch),thrk2,wtmp) 
         uu = Theta(rk2(jrch),thrk2,wtmp) * (soxy - o2con)
         if (algcon.gt.0.001) then
	                vv = (ai3 * Theta(gra,thgra,wtmp) - ai4           
     *                              * Theta(rhoq,thrho,wtmp)) * algcon
  
	   else
	   algcon=0.001
	   end if
         ww = Theta(rk1(jrch),thrk1,wtmp) * cbodcon

         if(rchdep.gt.0.001) xx = Theta(rk4(jrch),thrk4,wtmp)           &
     &	   / (rchdep * 1000.)
         if (nh3con.gt.0.001) then
	   yy = ai5 * Theta(bc1mod,thbc1,wtmp) * nh3con
         else
	   nh3con=0.001
	   end if
	   if (no2con.gt.0.001) then
	   zz = ai6 * Theta(bc2mod,thbc2,wtmp) * no2con
         else
	   n02con=0.001
	   end if
         ddisox = o2con + (uu + vv - ww - xx - yy - zz) * tday
		o2proc=o2con-ddisox
         if (ddisox < 0.00001) ddisox = 0.00001

 
!! end oxygen calculations

!! nitrogen calculations
         !! calculate organic N concentration at end of day
         !! QUAL2E section 3.3.1 equation III-16
         xx = 0.
         yy = 0.
         zz = 0.
         xx = ai1 * Theta(rhoq,thrho,wtmp) * algcon
         yy = Theta(bc3(jrch),thbc3,wtmp) * orgncon
         zz = Theta(rs4(jrch),thrs4,wtmp) * orgncon
         dorgn = 0.
         dorgn = orgncon + (xx - yy - zz) * tday
         if (dorgn < 0.00001) dorgn = 0.00001


        !! calculate fraction of algal nitrogen uptake from ammonia
        !! pool QUAL2E equation III-18
        f1 = 0.
        f1 = p_n * nh3con / (p_n * nh3con + (1. - p_n) * no3con +       &
     &                                                            1.e-6)

        !! calculate ammonia nitrogen concentration at end of day
        !! QUAL2E section 3.3.2 equation III-17
        ww = 0.
        xx = 0.
        yy = 0.
        zz = 0.
        ww = Theta(bc3(jrch),thbc3,wtmp) * orgncon
        xx = Theta(bc1mod,thbc1,wtmp) * nh3con
        yy = Theta(rs3(jrch),thrs3,wtmp) / (rchdep * 1000.)
        zz = f1 * ai1 * algcon * Theta(gra,thgra,wtmp)
        dnh4  = 0.
        dnh4  = nh3con + (ww - xx + yy - zz) * tday
        if (dnh4  < 1.e-6) dnh4  = 0.

         !! calculate concentration of nitrite at end of day
        !! QUAL2E section 3.3.3 equation III-19
        yy = 0.
        zz = 0.
        yy = Theta(bc1mod,thbc1,wtmp) * nh3con
        zz = Theta(bc2mod,thbc2,wtmp) * no2con
        dno2 = 0.
        dno2 = no2con + (yy - zz) * tday
        if (dno2 < 1.e-6) dno2 = 0.



        !! calculate nitrate concentration at end of day
        !! QUAL2E section 3.3.4 equation III-20
        yy = 0.
        zz = 0.
        yy = Theta(bc2mod,thbc2,wtmp) * no2con
        zz = (1. - f1) * ai1 * algcon * Theta(gra,thgra,wtmp)
        dno3  = 0.
        dno3  = no3con + (yy - zz) * tday
        if (dno3 < 1.e-6) dno3  = 0.
!! end nitrogen calculations

!! phosphorus calculations
        !! calculate organic phosphorus concentration at end of
        !! day QUAL2E section 3.3.6 equation III-24
        xx = 0.
        yy = 0.
        zz = 0.
        xx = ai2 * Theta(rhoq,thrho,wtmp) * algcon
        yy = Theta(bc4(jrch),thbc4,wtmp) * orgpcon
        zz = Theta(rs5(jrch),thrs5,wtmp) * orgpcon
        dorgp = 0.
        dorgp= orgpcon + (xx - yy - zz) * tday
        if (dorgp < 1.e-6) dorgp = 0.

        !! calculate dissolved phosphorus concentration at end
        !! of day QUAL2E section 3.4.2 equation III-25
        xx = 0.
        yy = 0.
        zz = 0.
        xx = Theta(bc4(jrch),thbc4,wtmp) * orgpcon
        yy = Theta(rs2(jrch),thrs2,wtmp) / (rchdep * 1000.)
        zz = ai2 * Theta(gra,thgra,wtmp) * algcon
        dsolp  = 0.
        dsolp  = solpcon + (xx + yy - zz) * tday
        if (dsolp  < 1.e-6) dsolp  = 0.




!! end phosphorus calculations






         wtrtot = wtrin + rchwtr
         !! initialize inflow concentrations
         chlin = 0.
         algin = 0.
         orgnin = 0.
         ammoin = 0.
         nitritin = 0.
         nitratin = 0.
         orgpin = 0.
         dispin = 0.
         cbodin = 0.
         disoxin = 0.
         cinn = 0.

         if (wtrin > 0.001) then
         chlin = 1000. * varoute(13,inum2) * (1. - rnum1) / wtrin
         algin = 1000. * chlin / ai0        !! QUAL2E equation III-1
         orgnin = 1000. * varoute(4,inum2) * (1. - rnum1) / wtrin
         ammoin = 1000. * varoute(14,inum2) * (1. - rnum1) / wtrin
         nitritin = 1000. * varoute(15,inum2) * (1. - rnum1) / wtrin
         nitratin = 1000. * varoute(6,inum2) * (1. - rnum1) / wtrin
         orgpin = 1000. * varoute(5,inum2) * (1. - rnum1) / wtrin
         dispin = 1000. * varoute(7,inum2) * (1. - rnum1) / wtrin
         cbodin = 1000. * varoute(16,inum2) * (1. - rnum1) / wtrin
         disoxin= 1000. * varoute(17,inum2) * (1. - rnum1) / wtrin
	   heatin = varoute(1,inum2) * (1. - rnum1)    
         end if


	   wattemp(jrch) =(heatin * wtrin + wtmp * rchwtr) / wtrtot
         algae(jrch) = (algin * wtrin + dalgae * rchwtr) / wtrtot

         organicn(jrch) = (orgnin * wtrin + dorgn * rchwtr) / wtrtot
         ammonian(jrch) = (ammoin * wtrin +  dnh4 * rchwtr) / wtrtot
         nitriten(jrch) = (nitritin * wtrin + dno2  * rchwtr) / wtrtot
         nitraten(jrch) = (nitratin * wtrin + dno3  * rchwtr) / wtrtot
         organicp(jrch) = (orgpin * wtrin +  dorgp * rchwtr) / wtrtot
         disolvp(jrch) = (dispin * wtrin +  dsolp * rchwtr) / wtrtot
         rch_cbod(jrch) = (cbodin * wtrin + dbod * rchwtr) / wtrtot
         rch_dox(jrch) = (disoxin * wtrin +  ddisox * rchwtr) / wtrtot

         !! calculate chlorophyll-a concentration at end of day
         !! QUAL2E equation III-1
         chlora(jrch) = 0.
         chlora(jrch) = algae(jrch) * ai0 / 1000.

      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
        algae(jrch) = 0.0
        chlora(jrch) = 0.0
        organicn(jrch) = 0.0
        ammonian(jrch) = 0.0
        nitriten(jrch) = 0.0
        nitraten(jrch) = 0.0
        organicp(jrch) = 0.0
        disolvp(jrch) = 0.0
        rch_cbod(jrch) = 0.0
        rch_dox(jrch) = 0.0
	   dalgae = 0.0
         dchla = 0.0
         dorgn = 0.0
         dnh4 = 0.0
         dno2 = 0.0
         dno3 = 0.0
         dorgp= 0.0
         dsolp = 0.0
         dbod = 0.0
         ddisox = 0.0   
        soxy = 0.0
      endif




      return
      end