subroutine ch_watqual4 !! ~ ~ ~ PURPOSE ~ ~ ~ !! this subroutine performs in-stream nutrient transformations and water !! quality calculations use channel_module use hydrograph_module use climate_module use channel_data_module use sd_channel_module use water_body_module implicit none real :: tday, wtmp, fll, gra real :: lambda, fnn, fpp, algi, fl_1, xx, yy, zz, ww real :: cordo, f1, algcon real :: thgra = 1.047, thrho = 1.047, thrs1 = 1.024 real :: thrs2 = 1.074, thrs3 = 1.074, thrs4 = 1.024, thrs5 = 1.024 real :: thbc1 = 1.083, thbc2 = 1.047, thbc3 = 1.047, thbc4 = 1.047 real :: thrk1 = 1.047, thrk2 = 1.024, thrk3 = 1.024, thrk4 = 1.060 real :: soxy !mg O2/L |saturation concetration of dissolved oxygen real :: theta, rs2_s, rs3_s, rk4_s real :: disoxin, dispin, ammoin, cinn real :: algin real :: factk, wq_semianalyt, alg_m1, alg_m, alg_m2 real :: alg_no3_m, alg_nh4_m, alg_p_m, alg_set, algcon_out real :: cbodo, cbodoin, rk1_k, wq_k2m, rk1_m , rk3_k, factm, & bc1_k,bc3_k, rs4_k, bc3_m real :: rk2_m, rk2_k, alg_m_o2, bc2_k, bc1_m, bc2_m, bc4_k, & bc4_m, rs5_k integer :: iwgn jrch = isdch !! calculate flow duration tday = rttime / 24.0 tday = amin1 (1., tday) rt_delt = 1. !! use maximum daily flow depth !rchdep = 0. !do istep = 1, time%step ! rchdep = Max (rchdep, flo_dep(istep)) !end do !! 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 = 5.0 + 0.75 * wst(iwst)%weat%tave if (wtmp <= 0.) wtmp = 0.1 ht2%temp = wtmp !! benthic sources/losses in mg rs2_s = Theta(ch_nut(jnut)%rs2,thrs2,wtmp) * ben_area !ch_hyd(jhyd)%l *ch_hyd(jhyd)%w * rt_delt rs3_s = Theta(ch_nut(jnut)%rs3,thrs3,wtmp) * ben_area !ch_hyd(jhyd)%l *ch_hyd(jhyd)%w * rt_delt rk4_s = Theta(ch_nut(jnut)%rk4,thrk4,wtmp) * ben_area !ch_hyd(jhyd)%l *ch_hyd(jhyd)%w * rt_delt !! ht3 = concentration of incoming nutrients !rchdep = ch_rcurv(jrch)%in2%dep if (ht3%flo > 0. .and. rchdep > 0.) then disoxin = ht3%dox - rk4_s / (ht3%flo * 1000.) !m3*1000 l/m3 = liters disoxin = max (0., disoxin) dispin = ht3%solp + rs2_s / (ht3%flo * 1000.) ammoin = ht3%nh3 + rs3_s / (ht3%flo * 1000.) !! calculate effective concentration of available nitrogen QUAL2E equation III-15 cinn = ch_stor(jrch)%nh3 + ch_stor(jrch)%no3 !! calculate saturation concentration for dissolved oxygen !! QUAL2E section 3.6.1 equation III-29 ww = -139.34410 + (1.575701e05 / (wtmp + 273.15)) xx = 6.642308e07 / ((wtmp + 273.15)**2) yy = 1.243800e10 / ((wtmp + 273.15)**3) zz = 8.621949e11 / ((wtmp + 273.15)**4) soxy = Exp(ww - xx + yy - zz) if (soxy < 1.e-6) soxy = 0. !! end initialize concentrations !! O2 impact calculations !! calculate nitrification rate correction factor for low oxygen QUAL2E equation III-21 if (ht3%dox < 0.001) ht3%dox = 0.001 if (ht3%dox > 30.) ht3%dox = 30. cordo = 1.0 - Exp(-0.6 * ht3%dox) !! end O2 impact calculations !! algal growth !! calculate chlorophyll-a concentration at end of day QUAL2E equation III-1 algcon = 1000. * ht3%chla / ch_nut(jnut)%ai0 algin = 1000. * ht3%chla / ch_nut(jnut)%ai0 !! calculate light extinction coefficient (algal self shading) QUAL2E equation III-12 if (ch_nut(jnut)%ai0 * algcon > 1.e-6) then lambda = ch_nut(jnut)%lambda0 + (ch_nut(jnut)%lambda1 * ch_nut(jnut)%ai0 * algcon) & + ch_nut(jnut)%lambda2 * (ch_nut(jnut)%ai0 * algcon) ** (.66667) else lambda = ch_nut(jnut)%lambda0 endif if (lambda > ch_nut(jnut)%lambda0) lambda = ch_nut(jnut)%lambda0 !! calculate algal growth limitation factors for nitrogen !! and phosphorus QUAL2E equations III-13 & III-14 fnn = cinn / (cinn + ch_nut(jnut)%k_n) fpp = ch_stor(jrch)%solp / (ch_stor(jrch)%solp + ch_nut(jnut)%k_p) !! calculate daylight average, photosynthetically active, light intensity QUAL2E equation III-8 !! Light Averaging Option # 2 iwgn = wst(iwst)%wco%wgn if (wgn_pms(iwgn)%daylth > 0.) then algi = wst(iwst)%weat%solrad * ch_nut(jnut)%tfact / wgn_pms(iwgn)%daylth else algi = 0.00001 end if !! calculate growth attenuation factor for light, based on !! daylight average light intensity QUAL2E equation III-7b fl_1 = (1. / (lambda * rchdep)) * Log((ch_nut(jnut)%k_l + algi) / & (ch_nut(jnut)%k_l + algi * (Exp(-lambda * rchdep)))) fll = 0.92 * (wgn_pms(iwgn)%daylth / 24.) * fl_1 !! calculcate local algal growth rate if (algcon < 5000.) then select case (ch_nut(jnut)%igropt) case (1) !! multiplicative QUAL2E equation III-3a gra = ch_nut(jnut)%mumax * fll * fnn * fpp case (2) !! limiting nutrient QUAL2E equation III-3b gra = ch_nut(jnut)%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 = ch_nut(jnut)%mumax * fll * 2. / ((1. / fnn) + (1. / fpp)) else gra = 0. endif end select end if !! calculate algal biomass concentration at end of day (phytoplanktonic algae) QUAL2E equation III-2 factk = Theta(gra,thgra,wtmp) - Theta(ch_nut(jnut)%rhoq, thrho, wtmp) algcon = 1000. * ht3%cla / ch_nut(jnut)%ai0 alg_m1 = wq_semianalyt (tday, rt_delt, 0., factk, algcon, algin) alg_m = wq_semianalyt (tday, rt_delt, 0., factk, algcon, algin) alg_m2 = alg_m - alg_m1 !! calculate fraction of algal nitrogen uptake from ammonia pool QUAL2E equation III-18 f1 = ch_nut(jnut)%p_n * ch_stor(jrch)%nh3 / (ch_nut(jnut)%p_n * ch_stor(jrch)%nh3 + & (1. - ch_nut(jnut)%p_n) * ch_stor(jrch)%no3 + 1.e-6) alg_no3_m = -alg_m * (1. - f1) * ch_nut(jnut)%ai1 alg_nh4_m = -alg_m * f1 * ch_nut(jnut)%ai1 alg_P_m = -alg_m * ch_nut(jnut)%ai2 alg_set = 0. if (rchdep > 0.001) alg_set = Theta (ch_nut(jnut)%rs1, thrs1, wtmp) / rchdep algcon_out = wq_semianalyt (tday, rt_delt, alg_m, -alg_set, algcon, algin) if (algcon_out < 1.e-6) algcon_out = 0. if (algcon_out > 5000.) algcon_out = 5000. !! calculate chlorophyll-a concentration at end of day QUAL2E equation III-1 ht2%chla = algcon_out * ch_nut(jnut)%ai0 / 1000. !! end algal growth !! oxygen calculations !! calculate carbonaceous biological oxygen demand at end of day QUAL2E section 3.5 equation III-26 !! adjust rk1 to m-term and BOD & O2 mass availability cbodo = min (ch_stor(jrch)%cbod, ch_stor(jrch)%dox) cbodoin = min (ht3%cbod, ht3%dox) rk1_k = -Theta (ch_nut(jnut)%rk1, thrk1,wtmp) rk1_m = wq_k2m (tday, rt_delt, rk1_k, ch_stor(jrch)%cbod, ht3%cbod) !! calculate corresponding m-term rk3_k=0. if (rchdep > 0.001) rk3_k = -Theta (ch_nut(jnut)%rk3, thrk3, wtmp) / rchdep factm = rk1_m factk = rk3_k ht2%cbod = wq_semianalyt (tday, rt_delt, factm, factk, ch_stor(jrch)%cbod, ht3%cbod) !! nitrogen calculations !! calculate organic N concentration at end of day !! QUAL2E section 3.3.1 equation III-16 bc1_k = Theta(ch_nut(jnut)%bc1,thbc1,wtmp) bc3_k = Theta(ch_nut(jnut)%bc3,thbc3,wtmp) bc1_k = bc1_k * 2. bc3_k = bc3_k * 2. rs4_k = 0. if (rchdep > 0.001) rs4_k = Theta (ch_nut(jnut)%rs4, thrs4, wtmp) / rchdep bc3_m = wq_k2m (tday, rt_delt, -bc3_k, ch_stor(jrch)%orgn, ht3%orgn) factk = -rs4_k factm = bc3_m ht2%orgn = wq_semianalyt (tday, rt_delt, factm, factk, ch_stor(jrch)%orgn, ht3%orgn) if (ht2%orgn <0.) ht2%orgn = 0. !! calculate dissolved oxygen concentration if reach at end of day QUAL2E section 3.6 equation III-28 rk2_m = Theta (ch_nut(jnut)%rk2, thrk2, wtmp) * soxy rk2_k = Theta (ch_nut(jnut)%rk2, thrk2, wtmp) alg_m_o2 = ch_nut(jnut)%ai4 * alg_m2 + ch_nut(jnut)%ai3 * alg_m1 factk = - rk2_k bc2_k = -Theta (ch_nut(jnut)%bc2, thbc2, wtmp) bc1_m = wq_k2m (tday, rt_delt, factk, ch_stor(jrch)%nh3, ammoin) bc2_m = wq_k2m (tday, rt_delt, bc2_k, ch_stor(jrch)%no2, ht3%no2) factm = rk1_m + rk2_m - rs4_k + bc1_m * ch_nut(jnut)%ai5 + bc2_m * ch_nut(jnut)%ai6 ht2%dox = wq_semianalyt (tday, rt_delt, factm, factk, ch_stor(jrch)%dox, ht3%dox) if (ht2%dox <0.) ht2%dox = 0. !! end oxygen calculations !! calculate ammonia nitrogen concentration at end of day QUAL2E section 3.3.2 equation III-17 factk = -bc1_k factm = bc1_m - bc3_m ht2%nh3 = wq_semianalyt (tday, rt_delt, factm, 0., ch_stor(jrch)%nh3, ammoin) if (ht2%nh3 < 1.e-6) ht2%nh3 = 0. !! calculate concentration of nitrite at end of day QUAL2E section 3.3.3 equation III-19 factm = -bc1_m + bc2_m ht2%no2 = wq_semianalyt (tday, rt_delt, factm, 0., ch_stor(jrch)%no2, ht3%no2) if (ht2%no2 < 1.e-6) ht2%no2 = 0. !! calculate nitrate concentration at end of day QUAL2E section 3.3.4 equation III-20 factk = 0. factm = -bc2_m ht2%no3 = wq_semianalyt (tday, rt_delt, factm, 0., ch_stor(jrch)%no3, ht3%no3) if (ht2%no3 < 1.e-6) ht3%no3 = 0. !! end nitrogen calculations !! phosphorus calculations !! calculate organic phosphorus concentration at end of day QUAL2E section 3.3.6 equation III-24 bc4_k = Theta (ch_nut(jnut)%bc4, thbc4,wtmp) bc4_m = wq_k2m (tday, rt_delt, -bc4_k, ch_stor(jrch)%sedp, ht3%sedp) rs5_k = 0. if (rchdep > 0.001) rs5_k = Theta (ch_nut(jnut)%rs5, thrs5, wtmp) / rchdep factk = -rs5_k factm = bc4_m ht2%sedp = wq_semianalyt (tday, rt_delt, factm, factk, ch_stor(jrch)%sedp, ht3%sedp) if (ht2%sedp < 1.e-6) ht2%sedp = 0. !! calculate dissolved phosphorus concentration at end of day QUAL2E section 3.4.2 equation III-25 factk = 0. factm = -bc4_m + ch_nut(jnut)%ai2 * alg_m ht2%solp = wq_semianalyt (tday, rt_delt, factm, 0., ch_stor(jrch)%solp, dispin) if (ht2%solp < 1.e-6) ht2%solp = 0. !! end phosphorus calculations else !! all water quality variables set to zero when no flow ht2 = hz ch_stor(jrch) = hz endif return end subroutine ch_watqual4