! FILE eros.f ! ! SUBROUTINES IN THIS FILE CALLED FROM ! subroutine ecklsp(j, je, k, n) hydrotop ! subroutine ysed(j) subbasin ! subroutine enrsb(j) subbasin ! subroutine orgnsed(j) subbasin ! subroutine psed(j) subbasin module erosion use utilities, only : dp implicit none ! weighted aver. st.org.N in lay 1 in subb., SUM(anors), g / t real(dp), save :: xnorg ! weighted aver. org.P in layer 1 in subb., SUM(porg), g / t real(dp), save :: xporg ! weighted aver. st.org.N in lay 1 in subb., SUM(anors), kg / ha real(dp), save :: xnorgp ! combined c, k, ls, p factor for subbasin, calc in subbasin real(dp), save :: xcklsp ! weighted aver. total P (porg + pms + pma) in subb., kg / ha real(dp), save :: xpsedp ! enrichment ratio for subbasin real(dp), save :: er ! xnorg * er, g / t: N org. in I layer for subbasin, corrected for enrichment real(dp), save :: conn ! xporg * er, g / t: P org. in I layer in subbasin, g / t real(dp), save :: cpp ! USLE erosion control practice factor P real(dp), save, dimension(:), allocatable :: ecp ! USLE slope length / slope steepness factor real(dp), save, dimension(:), allocatable :: sl ! org. N loss with erosion (calc in orgnsed), kg / ha real(dp), save, dimension(:), allocatable :: yone ! org. P loss with erosion, kg / ha real(dp), save, dimension(:), allocatable :: yphe ! particle size distribution, calc in subbasin real(dp), save, dimension(:, :), allocatable :: parsz ! delivery ratios (part. size distr. for sediments) real(dp), save, dimension(:, :), allocatable :: pct ! combined c, k, ls, p factor real(dp), save, dimension(:, :), allocatable :: cklsp ! minimum value of C factor for water erosion, readcrp real(dp), save, dimension(:), allocatable :: cvm ! correction coef. for chnnel USLE K factor real(dp), save :: chxkc0 = 1. ! correction coef. for chnnel USLE C factor real(dp), save :: chcc0 = 1. ! sediment conc in return flow, ppm real(dp), save, dimension(:), allocatable :: css ! channel USLE K factor real(dp), save, dimension(:), allocatable :: chxk ! channel USLE C factor real(dp), save, dimension(:), allocatable :: chc integer, save :: nsz = 5 real(dp), save, dimension(5) :: dia = (/ 14.14, 3.16, 1.41, 5.48, 22.36 /) namelist / EROSION_PARAMETERS / & chcc0, & chxkc0 contains subroutine erosion_initialise(mb, mcrdb, meap, subbasin_input_file_id, da, flu, stp, slope_length, tc, tp6, tp5, al, chl, chnn, chd, chss) use input, only : get_config_fid integer, intent(in) :: mb, mcrdb, meap, subbasin_input_file_id real(dp), intent(in) :: da real(dp), intent(in) :: flu(:) real(dp), intent(in) :: stp(:) real(dp), intent(in) :: slope_length(:) real(dp), intent(in) :: tc(:) real(dp), intent(in) :: tp6(:) real(dp), intent(in) :: tp5(:) real(dp), intent(out) :: al(:) real(dp), intent(in) :: chl(:, :) real(dp), intent(in) :: chnn(:) real(dp), intent(in) :: chd(:) real(dp), intent(in) :: chss(:) integer j real(dp) vm, a2, xm, sxc read(get_config_fid(), EROSION_PARAMETERS) call erosion_allocate(mb, mcrdb, meap) call erosion_read_input(subbasin_input_file_id) !**** CALC subbasin parameters sl(), al(), css() do j = 1, mb xm = .6 * (1. - exp(- 35.835 * stp(j))) sl(j) = (slope_length(j) / 22.127) ** xm * (65.41 * stp(j) * stp(j) + 4.56 * stp(j) + .065) a2 = - log10(tp5(j) / tp6(j)) / 1.0792 al(j) = tp6(j) * (tc(j) / 6.) ** a2 * (da / 3.6) * flu(j) / tp5(j) if (chl(2, j) .ne. 0.) then sxc = sqrt(chss(j)) / chnn(j) vm = chd(j) ** .6667 * sxc end if css(j) = css(j) * 1.e-6 end do end subroutine erosion_initialise subroutine erosion_read_input(subbasin_input_file_id) use input, only : read_real_column integer, intent(in) :: subbasin_input_file_id call read_real_column(subbasin_input_file_id, "css", css, 0.5_dp) call read_real_column(subbasin_input_file_id, "ecp", ecp, 0.5_dp) call read_real_column(subbasin_input_file_id, "chxk", chxk, 0.28_dp) call read_real_column(subbasin_input_file_id, "chc", chc, 1.0_dp) end subroutine erosion_read_input subroutine erosion_allocate(mb, mcrdb, meap) integer, intent(in) :: mb, mcrdb, meap ! NOTE: Something during the init process sets `mch = mb`, so use mb as we ! we have reordered some intialization code ! allocate(chxk(mch)) allocate(chc(mb)) allocate(chxk(mb)) allocate(cklsp(mb, meap)) allocate(css(mb)) allocate(cvm(mcrdb)) allocate(ecp(mb)) allocate(parsz(5, mb)) allocate(pct(5, mb)) allocate(sl(mb)) allocate(yone(mb)) allocate(yphe(mb)) chc = 1.0 * chcc0 chxk = 0.28 * chxkc0 cklsp = 0. css = 0.5 cvm = 0. ecp = 0.5 parsz = 0. pct = 0. sl = 0. yone = 0. yphe = 0. end subroutine erosion_allocate subroutine dealloc_erosion deallocate(chc) deallocate(chxk) deallocate(cklsp) deallocate(css) deallocate(cvm) deallocate(ecp) deallocate(parsz) deallocate(pct) deallocate(sl) deallocate(yone) deallocate(yphe) end subroutine dealloc_erosion subroutine erosion_cklsp_factor(j, je, k, cva, ek, igro, nucr, is_cropland, is_natural_vegetation, is_forest) !**** PURPOSE: THIS SUBROUTINE ESTIMATES COMBINED CKLSP factor ! FOR WATER EROSION !**** CALLED IN: HYDROTOP !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! PARAMETERS & VARIABLES ! ! >>>>> COMMON PARAMETERS & VARIABLES ! cklsp(j, je) = combined c, k, ls, p factor ! cva(j, je) = land cover, kg/ha, calc in crpmd ! cvm(icr) = minimum value of C factor for water erosion, readcrp ! dm(j, je) = total biomass, kg/ha ! ecp(j) = P factor, readsub ! ek(k) = USLE soil K factor, read in readsol ! ida = current day ! ieros = switch code to print from eros() ! iersb = number of subbasin to print from eros(), if ieros = 1 ! igro(j, je) = vegetation index for cropland, = 1 if planted ! nucr(j, je) = crop number ! sl(j) = USLE slope length/slope steepness factor ! >>>>> ! >>>>> STATIC PARAMETER ! c = C (land cover) factor ! >>>>> !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ real(dp), dimension(:, :), intent(in) :: cva real(dp), dimension(:), intent(in) :: ek integer, dimension(:, :), intent(in) :: igro integer, dimension(:, :), intent(in) :: nucr logical, intent(in) :: is_cropland, is_natural_vegetation, is_forest integer j, je, k real(dp) c c = 0. if (is_natural_vegetation) c = 0.1 if (is_forest) c = 0.45 if (is_cropland) then if (igro(j, je) .eq. 1) then c = exp((- .2231 - cvm(nucr(j, je))) * & exp(- .00115 * cva(j, je)) + cvm(nucr(j, je))) else c = 0.8 endif endif !**** CALC combined c, k, ls, p factor cklsp(j, je) = c * ek(k) * ecp(j) * sl(j) return end subroutine erosion_cklsp_factor !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& subroutine erosion_soil_loss(aff, pr, xqd, yd) !**** PURPOSE: THIS SUBROUTINE ESTIMATES DAILY SOIL LOSS CAUSED BY WATER ! EROSION USING THE MODIFIED UNIVERSAL SOIL LOSS EQUATION !**** CALLED IN: SUBBASIN !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! PARAMETERS & VARIABLES ! ! >>>>> COMMON PARAMETERS & VARIABLES ! aff = 1000. * da * flu(j), calc in subbasin ! ida = current day ! ieros = switch code to print from eros() ! iersb = number of subbasin to print from eros(), if ieros = 1 ! pr = peak runoff rate, m3/sec., calc in peakq ! precip = precipitation, mm ! sl(j) = USLE slope length/slope steepness factor ! xcklsp = combined c, k, ls, p factor for subbasin, calc in subbasin ! xqd = surface runoff for subbasin, mm, calc in subbasin ! yd = daily soil loss from subbasin caused by water erosion, t ! >>>>> !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !**** Include common parameters real(dp), intent(in) :: aff real(dp), intent(in) :: pr real(dp), intent(in) :: xqd real(dp), intent(inout) :: yd yd = 11.8 * (xqd * pr * aff) ** .56 * xcklsp return end subroutine erosion_soil_loss subroutine erosion_enritchment_ratio(j, da, da9, flu, pr, precip, rp, xqd, yd) !**** PURPOSE: THIS SUBROUTINE CALCULATES ENRICHMENT RATIO FOR SUBBASIN !**** CALLED IN: SUBBASIN !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! PARAMETERS & VARIABLES ! ! >>>>> COMMON PARAMETERS & VARIABLES ! da = area of the basin, km2 ! da9 = 100. * da = basin area in ha, from readbas ! er = enrichment ratio for subbasin ! flu(j) = fraction of subbasin area in the basin ! parsz(5, j) = particle size distribution, calc in subbasin ! pct(5, j) = delivery ratios (part. size distr. for sediments) ! pr = peak runoff rate, m3/sec., in peakq ! precip = precipitation, mm ! rp = alpha for rainfall, the fraction of total rainfall ! occuring during 0.5h ! xqd = daily surface runoff, mm, calc in subbasin ! yd = daily soil loss, t, calc in ysed ! >>>>> ! >>>>> STATIC PARAMETERS ! bet = local par ! cy = local par ! dia = local par ! dr = local par ! durf = local par ! jj = local par ! nsz = local par ! rep = local par ! rinf = local par ! rp1 = local par ! tot = local par ! x1 = local par ! x2 = local par ! >>>>> !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !**** Include common parameters real(dp), intent(in) :: da real(dp), intent(in) :: da9 real(dp), dimension(:), intent(in) :: flu real(dp), intent(in) :: pr real(dp), intent(in) :: precip real(dp), intent(in) :: rp real(dp), intent(in) :: xqd real(dp), intent(in) :: yd integer j, jj real(dp) bet, cy, dr, durf, rep, rinf, rp1, tot, x1, x2 rp1 = - 2. * precip *log(1. - rp) durf = 4.605 * precip / rp1 rinf = (precip - xqd) / durf rep = (rp1 - rinf) * da * flu(j) / 3.6 dr = (pr / (rep + 1.e-10)) ** .56 bet =log(dr) / 4.47 tot = 0. do jj = 1, nsz pct(jj, j) = parsz(jj, j) * exp(bet * dia(jj)) tot = tot + pct(jj, j) end do do jj = 1, nsz pct(jj, j) = pct(jj, j) / tot end do cy = 100000. * yd / (da9 * xqd + 1.e-6) x2 = - log10(dr) / 2.699 x1 = 1. / .25 ** x2 !**** CALC ENRICHMENT RATIO FOR SUBBASIN if (cy > 1.e-6) then er = .78 * cy ** (- .2468) else er = 0. endif if (er .lt. 1.) er = 1. if (er .gt. 3.5) er = 3.5 return end subroutine erosion_enritchment_ratio !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& subroutine erosion_organic_nitrate_loss(j, da9, yd, yon) !**** PURPOSE: COMPUTES organic N loss with erosion !**** CALLED IN: SUBBASIN !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! PARAMETERS & VARIABLES ! ! >>>>> COMMON PARAMETERS & VARIABLES ! conn = xnorg * er, g/t ! da9 = 100. * da = basin area in ha, from readbas ! er = enrichment ratio, from enrsb ! xnorg = N org. in I layer for subbasin, g/t ! xnorgp = N org. in I layer for subbasin, kg/ha ! yd = daily soil loss (erosion), in t, calc in ysed ! yon = org N loss with erosion, kg/ha ! yone(j) = org N loss with erosion, kg/ha ! >>>>> !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !**** Include common parameters real(dp), intent(in) :: da9 real(dp), intent(in) :: yd real(dp), intent(inout) :: yon integer j !**** CALC org N loss with erosion conn = xnorg * er yon = .001 * conn * yd / da9 !**** Correction: AnjaH !!! if (yon.lt.xnorgp) yon = 0. if (yon .gt. xnorgp) yon = xnorgp if (yon .le. 0) yon = 0. yone(j) = yon return end subroutine erosion_organic_nitrate_loss !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& subroutine erosion_phosphorus_loss(j, da9, yd, yph) !**** PURPOSE: COMPUTES P loss with erosion !**** CALLED IN: SUBBASIN !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! PARAMETERS & VARIABLES ! ! >>>>> COMMON PARAMETERS & VARIABLES ! cpp = xporg * er, g/t ! da9 = 100. * da = basin area in ha, from readbas ! er = enrichment ration, from enrsb ! xporg = P org. in I layer in subbasin, g/t ! xpsedp = SUM(porg+pms+pma) in subbasin, kg/ha ! yd = daily soil loss, in t, calc in ysed ! yph = P org. loss with erosion, kg/ha ! yphe(j) = P org. loss with erosion, kg/ha ! >>>>> !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !**** Include common parameters real(dp), intent(in) :: da9 real(dp), intent(in) :: yd real(dp), intent(inout) :: yph integer j !**** CALC P org. loss with erosion cpp = xporg * er yph = .001 * cpp * yd / da9 !**** Correction: AnjaH !!! if (yph.lt.xporgp) yph = 0. if (yph .gt. xpsedp) yph = xpsedp if (yph .le. 0) yph = 0. yphe(j) = yph return end subroutine erosion_phosphorus_loss end module erosion