! FILE readsol.f ! ! SUBROUTINES IN THIS FILE CALLED FROM ! subroutine readsol main module soil use utilities, only : dp, path_max_length, log_debug, log_info, log_warn, log_error implicit none ! correction factor for saturated conductivity (all soils) real(dp), save, dimension(:), allocatable :: bsn_sccor ! correction factor for curve number (0.25 - 1.25) real(dp), save, dimension(:), allocatable :: bsn_cncor ! fraction of field capacity real(dp), save, dimension(:), allocatable :: ffc ! ml = max number of soil layers integer :: ml = 10 ! msdb = max number of soils in database (used in readsol), counted during runtime integer, save :: msdb = 0 ! total number of soils used (ms), counted during runtime integer, save :: ms = 0 ! parameter for calc of precip alpha factor in alpha real(dp) :: ab = 0.02083 ! parameter used to calc P availability index in pcycle real(dp) :: psp = 0.5 ! active N pool fraction real(dp) :: rtn = 0.15 real(dp), save, dimension(17) :: bsc = & (/ - 8.96847, 0., - 0.028212, 19.52348, 0.00018107, & - 0.0094125, - 8.395215, 0., 0.077718, 0., 0.0000173, 0.02733, & 0.001434, - 0.0000035, 0., - 0.00298, - 0.019492 /) ! random number, used in subbasin to calc alpha real(dp), save :: rn1 ! number of soil layers, calc in subbasin, cycle 100 integer, save :: nn ! alpha for rainfall, the fraction of total rainfall real(dp), save :: r1 ! alpha for rainfall, the fraction of total rainfall real(dp), save :: rp ! current CN in hydrotope real(dp), save :: cn ! daily surface runoff, mm real(dp), save :: qd ! water excess above FC, mm real(dp), save :: su ! percolation, calc in perc for 4mm slugs, recalc here, mm real(dp), save :: sep ! lateral subsurface flow, mm, calc in perc for 4mm layers real(dp), save :: prk ! surface runoff for subbasin, mm, calc in subbasin real(dp), save :: xqd ! peak runoff rate, m3 / sec. real(dp), save :: pr ! index for channel parameters chl(, ), chw(, ), chk(, ) integer, save :: iv ! runoff volume = xqd * da * flu(j) * 1000, m3, calc in subbasin real(dp), save :: vo ! flow duration, h real(dp), save :: dur ! remember old xqd, mm, to compare with new one in subbasin real(dp), save :: q1 ! code for saturated cond.: 0 - read, 1 - calc integer, save :: isc ! switch code for CN: 0: CN dif for soils (standard method) integer, save :: icn ! init. water content coef., later stin() = fc() * stinco real(dp), save :: stinco = 0.9 ! init. CN for cropland, cond 1 real(dp), save :: cnum1 = 2. ! init. CN for cropland, cond 3 real(dp), save :: cnum3 = 2. ! coefs for each subbasin real(dp), save, dimension(:), allocatable :: cncor ! soil erodibility correction factor real(dp), save :: ekc0 = 1. ! coefs for each subbasin real(dp), save, dimension(:), allocatable :: sccor real(dp), save :: sccor0 = 1. ! average annual air temp, degree C, used in solt() real(dp), save, dimension(:), allocatable :: avt ! annual 1/2 amplitude in daily average temp, calc readwet real(dp), save, dimension(:), allocatable :: amp ! sum of overland and channel concentration times, hours real(dp), save, dimension(:), allocatable :: tc ! fun(tc, tp5, tp6, flu, da), calc in readsub real(dp), save, dimension(:), allocatable :: al ! initial land cover, kg / ha real(dp), save, dimension(:), allocatable :: cv ! monthly prob. of rainy day, calc readwet real(dp), save, dimension(:, :), allocatable :: wft ! wi(m, j) = f(wim), used in alpha() for estim of precip. alpha factor real(dp), save, dimension(:, :), allocatable :: wi ! water storage in a layer, recalc here, mm real(dp), save, dimension(:, :, :), allocatable :: ste ! daily ave temp at the bottom of each layer, degree C real(dp), save, dimension(:, :, :), allocatable :: te ! bare soil surface temp, degree real(dp), save, dimension(:, :), allocatable :: te0 ! lateral subsurface flow, sum of prk for layer, mm real(dp), save, dimension(:, :, :), allocatable :: flate ! percolation, sum of sep for layer, mm real(dp), save, dimension(:, :, :), allocatable :: poe ! retention factor, corresponding cn1 real(dp), save, dimension(:, :), allocatable :: smx ! shape parameters for calc. of retention real(dp), save, dimension(:, :, :), allocatable :: wf ! return flow travel time, calc. hrflowtt, h real(dp), save, dimension(:, :, :), allocatable :: hrtc ! fun(field cap), calc in subbasin from hsumfc(j, jea) real(dp), save, dimension(:, :, :), allocatable :: hwss ! position of the soil (storing the soil number) integer, save, dimension(:), allocatable :: snum ! number of soil layers integer, save, dimension(:), allocatable :: ns integer, save, dimension(:), allocatable :: nsolgr real(dp), save, dimension(:), allocatable :: ek real(dp), save, dimension(:), allocatable :: bd real(dp), save, dimension(:), allocatable :: abd real(dp), save, dimension(:), allocatable :: swin real(dp), save, dimension(:), allocatable :: sumul real(dp), save, dimension(:), allocatable :: sumfc ! number of layers for arable soils integer, save, dimension(:), allocatable :: nsa ! silt content, % real(dp), save, dimension(:, :), allocatable :: silt ! clay content, % real(dp), save, dimension(:, :), allocatable :: clay ! sand content, % real(dp), save, dimension(:, :), allocatable :: sand ! particle size distribution real(dp), save, dimension(:, :), allocatable :: psz real(dp), save, dimension(:, :), allocatable :: z real(dp), save, dimension(:, :), allocatable :: cla real(dp), save, dimension(:, :), allocatable :: sil real(dp), save, dimension(:, :), allocatable :: san real(dp), save, dimension(:, :), allocatable :: bden real(dp), save, dimension(:, :), allocatable :: poros real(dp), save, dimension(:, :), allocatable :: awc real(dp), save, dimension(:, :), allocatable :: fc real(dp), save, dimension(:, :), allocatable :: cbn real(dp), save, dimension(:, :), allocatable :: wn real(dp), save, dimension(:, :), allocatable :: sc real(dp), save, dimension(:, :), allocatable :: wno3 real(dp), save, dimension(:, :), allocatable :: ap real(dp), save, dimension(:, :), allocatable :: wmn real(dp), save, dimension(:, :), allocatable :: wpo real(dp), save, dimension(:, :), allocatable :: pmn real(dp), save, dimension(:, :), allocatable :: op real(dp), save, dimension(:, :), allocatable :: hum real(dp), save, dimension(:, :), allocatable :: wp real(dp), save, dimension(:, :), allocatable :: up real(dp), save, dimension(:, :), allocatable :: por real(dp), save, dimension(:, :), allocatable :: stin real(dp), save, dimension(:, :), allocatable :: ul real(dp), save, dimension(:, :), allocatable :: hk ! Curve Numbers for hydrotope j, jea real(dp), save, dimension(:, :), allocatable :: cn2 integer, save, dimension(4) :: k7 = (/ 227, 57, 929, 37/) ! accumulated soil water content in the catchment at the end of the year, mm real(dp), save :: soil_acc_mm real(dp), save :: amt = 4. character(len=path_max_length), dimension(:), allocatable :: soil_filenames character(len=path_max_length) :: soil_input_dir = 'input/soils' integer :: soil_input_file_id ! Output variable ids integer :: & surface_runoff_output_id = 0, & subsurface_runoff_output_id = 0, & percolation_output_id = 0, & soil_water_index_output_id = 0 character(len=path_max_length) :: soil_input_file = "soil.csv" namelist / SOIL_PARAMETERS / & ab, & cnum1, & cnum3, & ekc0, & icn, & ml, & psp, & rtn, & sccor0, & stinco, & soil_input_dir contains subroutine soil_initialise(mb, meap, mstruc, neap, bSubcatch, da, slope_length, stp, ovn, chs, flu, chn, chl, cn2a, cn2b, cn2c, cn2d) use input, only : input_open_file, input_count_rows, get_config_fid use output, only : output_register_hydrotope_var use utilities, only : log_error, random_n integer, intent(in) :: mb, meap integer, dimension(:, :, :), intent(in) :: mstruc integer, intent(inout) :: neap(:) logical, intent(in) :: bSubcatch real(dp), intent(in) :: da real(dp), intent(in) :: slope_length(:) real(dp), intent(in) :: stp(:) real(dp), intent(in) :: ovn(mb) real(dp), intent(in) :: chs(mb) real(dp), intent(in) :: flu(mb) real(dp), intent(in) :: chn(mb) real(dp), intent(in) :: chl(2, mb) real(dp), dimension(:), intent(in) :: cn2a, cn2b, cn2c, cn2d integer j real(dp) ot, ct read(get_config_fid(), SOIL_PARAMETERS) if (bSubcatch .AND. icn == 1) call log_error("soil_initialise", & "Inconsistency in .bsn input file: "// & "Switch icn is set to 1 and subcatchment calibration is switched ON. "// & "Either set icn = 0 or deactivate subcatchment calibration.") if (isc.eq.0) call log_debug("soil_initialise", 'Codes: SC - read, isc =', int=isc) if (isc.eq.1) call log_debug("soil_initialise", 'Codes: SC - calculated, isc =', int=isc) if (icn.eq.0) call log_debug("soil_initialise", 'Codes: CN - standard, icn =', int=icn) if (icn.eq.1) call log_debug("soil_initialise", 'Codes: CN = f(cnum1, cnum3), icn =', int=icn) surface_runoff_output_id = output_register_hydrotope_var("surface_runoff", .False.) subsurface_runoff_output_id = output_register_hydrotope_var("subsurface_runoff", .False.) percolation_output_id = output_register_hydrotope_var("percolation", .False.) soil_water_index_output_id = output_register_hydrotope_var("soil_water_index", .False.) soil_input_file_id = input_open_file(soil_input_file) ! **** count number of soils used in simulation (soil.cio) msdb = input_count_rows(soil_input_file_id, .true.) ! if soil parms read from single soil files then no header ms = msdb call soil_allocate(mb, meap) call soil_read_input(mb, mstruc, neap, cn2a, cn2b, cn2c, cn2d) !**** COMPUTE TIME OF CONCENTRATION for basins tc() do j = 1, mb ot = .0556 * (slope_length(j) * ovn(j)) ** .6 / stp(j) ** .3 ct = .62 * chl(1, j) * chn(j) ** .75 / ((da * flu(j)) ** .125 * chs(j) ** .375) tc(j) = ot + ct end do rn1 = random_n(k7) end subroutine soil_initialise subroutine soil_allocate(mb, meap) integer, intent(in) :: mb, meap allocate(abd(ms)) allocate(al(mb)) allocate(amp(mb)) allocate(ap(ml, ms)) allocate(avt(mb)) allocate(awc(ml, ms)) allocate(bd(ms)) allocate(bden(ml, ms)) allocate(cbn(ml, ms)) allocate(cla(ml, ms)) allocate(clay(1, ms)) allocate(cn2(mb, meap)) allocate(cv(mb)) allocate(ek(ms)) allocate(fc(ml, ms)) allocate(ffc(mb)) allocate(flate(mb, meap, ml)) allocate(hk(ml, ms)) allocate(hrtc(mb, meap, ml)) allocate(hum(ml, ms)) allocate(hwss(2, mb, meap)) allocate(ns(ms)) ! can be deleted if nslayers is fully implemented allocate(nsa(ms)) allocate(nsolgr(ms)) allocate(op(ml, ms)) allocate(pmn(ml, ms)) allocate(poe(mb, meap, ml)) allocate(por(ml, ms)) allocate(poros(ml, ms)) allocate(psz(5, ms)) allocate(san(ml, ms)) allocate(sand(1, ms)) allocate(sc(ml, ms)) allocate(sccor(mb)) allocate(sil(ml, ms)) allocate(silt(1, ms)) allocate(smx(mb, meap)) allocate(snum(ms)) allocate(soil_filenames(msdb)) allocate(ste(mb, meap, ml)) allocate(stin(ml, ms)) allocate(sumfc(ms)) allocate(sumul(ms)) allocate(swin(ms)) allocate(tc(mb)) allocate(te(mb, meap, ml)) allocate(te0(mb, meap)) allocate(ul(ml, ms)) allocate(up(ml, ms)) allocate(wf(2, mb, meap)) allocate(wft(12, mb)) allocate(wi(12, mb)) allocate(wmn(ml, ms)) allocate(wn(ml, ms)) allocate(wno3(ml, ms)) allocate(wp(ml, ms)) allocate(wpo(ml, ms)) allocate(z(ml, ms)) allocate(cncor(mb)) cncor = 0. abd = 0. abd = 0. al = 0. amp = 0. ap = 0. avt = 0. awc = 0. bd = 0. bden = 0. cbn = 0. cla = 0. clay = 0. cn2 = 0. cv = 0.01 ek = 0. fc = 0. ffc = 0. flate = 0. flate = 0. hk = 0. hrtc = 0. hrtc = 0. hum = 0. hwss = 0. ns = 0 nsa = 0 nsolgr = 0 op = 0. pmn = 0. poe = 0. por = 0. poros = 0. psz = 0. psz = 0. san = 0. sand = 0. sc = 0. sccor = sccor0 sil = 0. silt = 0. smx = 0. snum = 0 ste = 0. ste = 0. stin = 0. sumfc = 0. sumfc = 0. sumul = 0. sumul = 0. swin = 0. tc = 0. te = 0. te0 = 0. ul = 0. up = 0. wf = 0. wft = 0. wi = 0. wmn = 0. wn = 0. wno3 = 0. wp = 0. wpo = 0. z = 0. end subroutine soil_allocate subroutine soil_allocate_subcatch(n_subcatch) integer, intent(in) :: n_subcatch allocate(bsn_sccor(n_subcatch)) bsn_sccor = 0. allocate(bsn_cncor(n_subcatch)) bsn_cncor = 0. end subroutine soil_allocate_subcatch subroutine dealloc_soil deallocate(abd) deallocate(al) deallocate(amp) deallocate(ap) deallocate(avt) deallocate(awc) deallocate(bd) deallocate(bden) deallocate(cbn) deallocate(cla) deallocate(clay) deallocate(cn2) deallocate(cv) deallocate(ek) deallocate(fc) deallocate(ffc) deallocate(flate) deallocate(hk) deallocate(hrtc) deallocate(hum) deallocate(hwss) deallocate(ns) deallocate(nsa) deallocate(op) deallocate(pmn) deallocate(poe) deallocate(por) deallocate(poros) deallocate(psz) deallocate(san) deallocate(sand) deallocate(sc) deallocate(sil) deallocate(silt) deallocate(smx) deallocate(soil_filenames) deallocate(ste) deallocate(stin) deallocate(sumfc) deallocate(sumul) deallocate(swin) deallocate(tc) deallocate(te) deallocate(te0) deallocate(ul) deallocate(up) deallocate(wf) deallocate(wft) deallocate(wi) deallocate(wmn) deallocate(wn) deallocate(wno3) deallocate(wp) deallocate(wpo) deallocate(z) end subroutine dealloc_soil subroutine soil_curve_number(cnn, j, je, hsumfc, hsumul) !**** PURPOSE: TO SET CURVE NUMBER PARAMETERS !**** CALLED IN: HYDROTOP !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! PARAMETERS & VARIABLES ! ! >>>>> COMMON PARAMETERS & VARIABLES ! cnn = Curve Number = cn2(k, n), from readsol, with title ! cnum1 = init. CN for cropland, cond 1 ! cnum3 = init. CN for cropland, cond 3 ! hsumfc(j, je) = sum of field capacity in soil, calc in subbasin, mm ! hsumul(j, je) = sum of upper limit water content in soil, ! calc in subbasin, mm ! icn = switch code for CN: 0 - CN dif for soils, ! 1 - CN=const from cnum1, cnum3 ! icurn = switch code to print from curn() ! icursb = number of subbasin to print from curn(), if icurn = 1 ! ida = current day ! smx(j, je) = retention factor, corresponding cn1 ! wf(2, j, je) = shape parameters for calc. of retention ! >>>>> ! >>>>> STATIC PARAMETERS ! c2 = local par ! cn1 = CN, conditions 1 ! cn3 = CN, conditions 3 ! s3 = local par ! yy = local par ! zz = local par ! >>>>> !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !**** Include common parameters real(dp), dimension(:, :), intent(in) :: hsumfc real(dp), dimension(:, :), intent(in) :: hsumul integer, intent(in) :: j, je ! cnn = cn2(k, n) real(dp), intent(in) :: cnn real(dp) :: c2, cn1, cn3, s3, yy, zz !**** CALC cn1, cn3 if (icn .eq. 0) then c2 = cnn * cnn cn1 = - 16.911 + 1.3481 * cnn - .013793 * c2 + .00011772 * c2 * cnn cn3 = cnn * exp(.006729 * (100. - cnn)) else cn1 = cnum1 cn3 = cnum3 endif !**** CALC retention factor smx() smx(j, je) = 254. * (100. / cn1 - 1.) s3 = 254. * (100. / cn3 - 1.) !**** CALC shape parameters for CN method wf(1, ) and wf(2, ) yy = hsumfc(j, je) / (1. - s3 / smx(j, je)) - hsumfc(j, je) if (yy .ne. 0.) then zz =log(yy) wf(2, j, je) = (zz -log(hsumul(j, je) / (1. - 2.54 / smx(j, je)) & - hsumul(j, je))) / (hsumul(j, je) - hsumfc(j, je)) wf(1, j, je) = zz + wf(2, j, je) * hsumfc(j, je) else wf(2, j, je) = 0. wf(1, j, je) = 0. endif return end subroutine soil_curve_number subroutine soil_curve_number_runoff(j, je, alai, blai, canstor, igro, nucr, preinf, cnmx) !**** PURPOSE: THIS SUBROUTINE COMPUTES DAILY RUNOFF GIVEN DAILY PRECIPITATION ! AND SNOW MELT USING A MODIFIED SCS CURVE NUMBER APPROACH !**** CALLED IN: HYDROTOP !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! PARAMETERS & VARIABLES ! ! >>>>> COMMON PARAMETERS & VARIABLES ! alai(j, je) = Leaf Area Index (LAI) ! blai(icr) = max LAI for crop ! canmax(n) = canopy maximum storage for interception, mm, calc in init ! canstor(j, je) = canopy water storage, mm ! cn = Curve Number, current ! icurn = switch code to print from curn() ! icursb = number of subbasin to print from curn(), if icurn = 1 ! ida = current day ! igro(j, je) = vegetation index, =1 if vegetation is growing ! nn = number of soil layers, calc in subbasin, cycle 100 ! nucr(j, je) = crop number (database) ! precip = precipitation, mm, read in readcli ! preinf(j, je) = precipitation adjusted for canopy storage, mm ! qd = daily surface runoff, mm ! smx(j, je) = retention coef, calc in curno ! ste(j, je, l) = water storage in a layer, mm, calc in hydrotop & purk ! te(j, je, l) = soil temperature, degree C, calc in solt ! wf(2, j, je) = shape parameters eq.6, calc in curno ! >>>>> ! >>>>> STATIC PARAMETERS ! bb = local par ! canmxl = local par ! l = local par ! pb = local par ! r2 = local par ! sum = soil water content in all layers ! xx = local par ! xx1 = local par ! xx3 = local par ! xx4 = local par ! >>>>> !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !**** Include common parameters real(dp), dimension(:, :), intent(in) :: alai real(dp), dimension(:), intent(in) :: blai real(dp), dimension(:, :), intent(inout) :: canstor integer, dimension(:, :), intent(in) :: igro integer, dimension(:, :), intent(in) :: nucr real(dp), dimension(:, :), intent(inout) :: preinf real(dp), intent(in) :: cnmx integer j, je, l real(dp) bb, canmxl, pb, r2, sum, xx, xx1, xx3, xx4 sum = 0. do l = 1, nn sum = sum + ste(j, je, l) end do !**** Canopy interception xx1 = 0. canmxl = 0. xx1 = preinf(j, je) !**** CALC canopy storage if (igro(j, je) .ge. 1) then canmxl = cnmx * alai(j, je) / blai(nucr(j, je)) else canmxl = 0. endif xx3 = preinf(j, je) - canmxl if (xx3 < 0.) then canstor(j, je) = xx1 else canstor(j, je) = canmxl endif preinf(j, je) = preinf(j, je) - canstor(j, je) xx4 = preinf(j, je) - canstor(j, je) xx = wf(1, j, je) - wf(2, j, je) * sum if (xx .lt. - 20.) xx = - 20. if (xx .gt. 20.) xx = 20. r2 = smx(j, je) * (1. - sum / (sum + exp(xx))) if (te(j, je, 2) .le. 0.) r2 = smx(j, je) * (1. - exp(- .000862 * r2)) cn = 25400. / (r2 + 254.) r2 = 25400. / cn - 254. bb = .2 * r2 pb = xx4 - bb !**** CALC daily surface runoff qd if (pb .gt. 0.) then qd = pb * pb / (xx4 + .8 * r2) else qd = 0. end if return end subroutine soil_curve_number_runoff !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& subroutine soil_curve_number_peak_runoff(j) !**** PURPOSE: THIS SUBROUTINE COMPUTES THE PEAK RUNOFF RATE ! USING A MODIFICATION OF THE RATIONAL FORMULA !**** CALLED IN: SUBBASIN !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! PARAMETERS & VARIABLES ! ! >>>>> COMMON PARAMETERS & VARIABLES ! al(j) = fun(tc, tp5, tp6, flu, da), calc in readsub ! pr = peak runoff rate, m3/sec. ! r1 = alpha for rainfall, the fraction of total rainfall ! occuring during 0.5h ! tc(j) = time of concentration, hours, calc in readsub ! xqd = surface runoff, mm, calc in volq ! >>>>> !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !**** Include common parameters integer j pr = r1 * al(j) * xqd / tc(j) return end subroutine soil_curve_number_peak_runoff !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& subroutine soil_curve_transmission_losses(j, chk, chl, chw) !**** PURPOSE: THIS SUBROUTINE COMPUTES CHANNEL TRANSMISSION LOSSES ! and recalculates surface runoff xqd and peak runoff rate pr !**** CALLED IN: SUBBASIN !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! PARAMETERS & VARIABLES ! ! >>>>> COMMON PARAMETERS & VARIABLES ! chk(2, j) = effective hydraulic conductivity of main channel, mm/h ! chl(2, j) = main channel length, km ! chw(2, j) = average width of main channel, m ! dur = flow duration, h ! ida = current day ! itran = switch code to print from tran() ! iv = 1 in subbasin (reach 1 as outlet) ! pr = peak runoff rate, m3/sec. ! q1 = remember old xqd, mm, to compare with new one in subbasin ! vo = runoff volume = xqd * da * flu(j) * 1000, m3, ! calc in subbasin ! xqd = surface runoff, mm, recalculation ! >>>>> ! >>>>> STATIC PARAMETERS ! a = local par ! ak = local par ! axw = local par ! b = local par ! bxw = local par ! pr1 = local par ! pxw = local par ! vol = local par ! xx = local par ! yy = local par ! zz = local par ! >>>>> !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !**** Include common parameters real(dp), dimension(:, :), intent(in) :: chk real(dp), dimension(:, :), intent(in) :: chl real(dp), dimension(:, :), intent(in) :: chw integer j real(dp) a, ak, axw, b, bxw, pr1, pxw, vol, xx, yy, zz ! q1 = xqd yy = vo / (xqd * 1000.) a = - .0001831 * chk(iv, j) * dur xx = .2649 * chk(iv, j) * dur / vo !**** RECALC xqd, pr ! iv = 1 (from subbasin) if (xx .lt. 1.) then ak = - 1.09 *log(1. - xx) b = exp(- ak) if (.not. ((1.-b) .le. 1.e-20)) then zz = - 2.04 * ak * chw(iv, j) * chl(iv, j) if (zz .ge. - 30.) then bxw = exp(zz) axw = (a / (1. - b)) * (1. - bxw) pxw = (- axw / bxw) vol = vo / 1233.5 xqd = 0. if (vol .gt. pxw) xqd = axw + bxw * vol xqd = xqd * 1.234 / yy pr1 = 35.3 * pr pr = 12.1 * axw / dur - (1. - bxw) * vol + bxw * pr1 pr = pr / 35.3 if (pr .lt. 0.) pr = 0. else xqd = 0. end if end if end if return end subroutine soil_curve_transmission_losses subroutine soil_curve_number_alpha(j, mo, precip, sml) !**** PURPOSE: THIS SUBROUTINE COMPUTES PRECIP ALPHA FACTOR: ! ALPHA is A DIMENSIONLESS PARAMETER THAT EXPRESSES ! THE FRACTION OF TOTAL RAINFALL THAT OCCURS DURING 0.5 H ! readwet: read wim() = monthly max .5h rain for period of record (mm) ! readwet: calc wi() = f(wim) = to be used in alpha() ! for estim of precip. alpha factor ! alpha: calc r1 = alpha for rainfall, the fraction of total rainfall ! occuring in 0.5h ! peakq: calc pr = peak runoff rate using r1 !**** CALLED IN: SUBBASIN !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! PARAMETERS & VARIABLES ! ! >>>>> COMMON PARAMETERS & VARIABLES !block ab = .02083 !block k7 = parameters for random number generator ! mo = current month ! precip = precipitation in subbasin, mm ! r1 = alpha for rainfall, the fraction of total rainfall ! occuring during 0.5h ! rn1 = random value, calc in readsub ! rp = alpha for rainfall, the fraction of total rainfall ! occuring during 0.5h ! sml = snow melt, calc in snom(), mm ! wi(m, j) = f(wim), wim = monthly max .5h rain for period of record, mm ! >>>>> ! >>>>> STATIC PARAMETERS ! aii = local parameter ! ajp = local parameter ! ei = local parameter ! >>>>> !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !**** Include common parameters use utilities, only : gamma_distribution integer, intent(in) :: mo real(dp), intent(in) :: precip real(dp), intent(in) :: sml integer j real(dp) aii, ajp, ei ! real(dp) gammad ei = precip - sml ! sl begin ei = MAX(ei, 0.0_dp) ! sl end aii = (1. - ab) / (wi(mo, j) - ab) if (ei .ge. 25.) then ajp = 1. - exp(- 125. / ei) else ajp = 1. end if r1 = gamma_distribution(rn1, aii, k7) r1 = (ei * (ab + r1 * (ajp - ab)) + sml * ab) / precip ! sl begin if (r1 >= 1.) r1 = 0.99999 ! sl end rp = r1 return end subroutine soil_curve_number_alpha subroutine soil_read_input(mb, mstruc, neap, cn2a, cn2b, cn2c, cn2d) !**** PURPOSE: THIS SUBROUTINE READS SOIL DATA !**** CALLED IN: MAIN !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! PARAMETERS & VARIABLES ! ! >>>>> COMMON PARAMETERS & VARIABLES ! abd(k) = sum of porosity for soil type, mm, used in solt ! alai0(n) = initial Leaf Area Index for land use type ! ap(l, k) = labile (soluble) phosphorus, g/t, then trans. to kg/ha ! asc(1:17) = empirical coef. (estimated here) to calc sc(l, k) ! awc(l, k) = available water capacity, % , /100. to transfer to m/m ! bd(k) = bulk density of the upper (!) soil layer, g/cm3, ! calc from porosity ! bden(l, k) = bulk density, input, g/cm3 ! bsc(1:17) = empirical coef to calc sc(l, k) ! cbn(l, k) = organic carbon content, % ! cla(l, k) = clay content, % ! clay(1, k) = clay content, % ! cn2(k, n, j) = Curve Numbers for soil k, land use n, and subbasin ! cn2a(n) = Curve Numbers for soil group A ! cn2b(n) = Curve Numbers for soil group B ! cn2c(n) = Curve Numbers for soil group C ! cn2d(n) = Curve Numbers for soil group D ! cv(j) = initial land cover, kg/ha ! ek(k) = soil erodibility factor ! fc(l, k) = field capacity, %, then trans. to mm ! hk(l, k) = beta coef. to estimate hydr. cond., used in perc ! hum(l, k) = humus content, kg/ha ! isc = code for saturated cond.: 0 - read, 1 - calc ! nn = number of soil layers ! ns(k) = number of soil layers ! nsa(k) = number of layers for arable soils !block nsolgr(k) = number of soils group for soil type, 1=A, 2=B, 3=C, 4=D ! op(l, k) = stable min. P, estimated from pmn(l, k), coef = 4., kg/ha ! pmn(l, k) = act. min. P, estimated from ap(l, k), psp=0.5, kg/ha ! por(l, k) = porosity, m/m (calc) ! poros(l, k) = porosity, % (input) !block psp = parameter = 0.5 ! psz(5, k) = particle size distribution !block rtn = parameter = .15 ! san(l, k) = sand content, % ! sand(1, k) = sand content, % ! sc(l, k) = saturated conductivity, mm/h, calc, if isc = 1 ! sccor = correction factor for saturated conductivity (all soils) ! sil(l, k) = silt content, % ! silt(1, k) = silt content, % ! snam(k) = soil name ! stin(l, k) = water content in soil layer, mm ! stinco = init. water content coef, from readbas ! sumfc(k) = sum of field capacity in soil, mm ! sumul(k) = sum of upper limit water content in soil, mm ! swin(k) = soil water content in soil profile, mm ! ul(l, k) = upper limit water content in layer, mm ! up(l, k) = upper limit water content, m/m (used only here to calc fc) ! wmn(l, k) = active org. N, estim. from wn(l, k), rtn=.15, kg/ha ! wn(l, k) = stable organic N content, %, ! recalc from % to g/t OR ! estim. from C using C:N=10. and 10000. coef from % to g/t, ! then trans. to kg/ha ! wno3(l, k) = NO3-N content (kg/ha), estimated in g/t, trans to kg/ha ! wp(l, k) = wilting point, in m/m (used here and in perc) ! wpo(l, k) = org. P, estimated from wn(), coef=.125, kg/ha ! z(l, k) = depth to bottom of layer, mm ! >>>>> ! >>>>> STATIC PARAMETERS ! ccf = correction coef for up and por by silt content ! dg = layer thickness ! i = cycle par ! j = par ! jj = par ! k = soil type number ! ka = par ! l = par ! n = par ! numsb = actual number of soils in database ! sfile = soil file name ! tex = sum of clay, silt and sand ! titldum = text ! wt1 = coef. to transform from g/t to kg/ha ! xx = par ! xz = par ! zz = par ! >>>>> !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !**** Include common parameters use utilities, only : open_file, check_range, to_string use input, only : read_string_column integer, intent(in) :: mb integer, dimension(:, :, :), intent(in) :: mstruc integer, dimension(:), intent(in) :: neap real(dp), dimension(:), intent(in) :: cn2a, cn2b, cn2c, cn2d ! soil counter integer :: k ! layer counter integer :: l integer :: jj, n integer :: j, jea character(len = 100) :: titldum, snam, sn character(len=path_max_length) :: soildat real(dp) :: dg, wt1, xx, xz, zz real(dp), dimension(17):: asc integer :: fid call log_debug('soil_read_input', 'SOIL PARAMETERS:') call log_debug('soil_read_input', 'Number of soils in soil.csv = ', int=msdb) !**** Reading Soil parameters from soilXX.dat files call read_string_column(soil_input_file_id, "soil_file", soil_filenames, "") do k = 1, msdb if (k .gt. msdb .or. k .eq. 0) then call log_error("soil_read_input", 'Soil type (k > msdb) = max soil number, or 0', & i1=k, i2=msdb) end if soildat = soil_filenames(k) ! open soil parameter files fid = open_file(trim(soil_input_dir)//'/'//soildat) !**** READ SOIL HYDRAULIC PROPERTIES ! THE SOIL IS DIVIDED VERTICALLY INTO UP TO 10 LAYERS ! LAYER THICKNESS CAN BE VARIABLE ! ns - number of layers, nsa - number of layers for arable soils read(fid, *, ERR=55) snum(k), ns(k), nsa(k), nsolgr(k) nn = ns(k) if (ns(k) .gt. ml .or. ns(k) .eq. 0 ) then call log_error("soil_read_input", 'Number of layers for (soil > ml) = max number, or 0', i1=k, i2=ml) end if if (nsa(k) .gt. ns(k) .or. nsa(k) .eq. 0 ) then call log_error("soil_read_input", 'Number of arrable layers grater than number of layers (soil) or 0', i1=k, ints=(/nsa(k), ns(k)/)) end if read(fid, "(a)") snam read(fid, "(a)") titldum read(fid, "(a)") titldum read(fid, *, ERR=55) (z(l, k), l = 1, nn) read(fid, *, ERR=55) (cla(l, k), l = 1, nn) read(fid, *, ERR=55) (sil(l, k), l = 1, nn) read(fid, *, ERR=55) (san(l, k), l = 1, nn) read(fid, *, ERR=55) (bden(l, k), l = 1, nn) read(fid, *, ERR=55) (poros(l, k), l = 1, nn) read(fid, *, ERR=55) (awc(l, k), l = 1, nn) read(fid, *, ERR=55) (fc(l, k), l = 1, nn) read(fid, *, ERR=55) (cbn(l, k), l = 1, nn) read(fid, *, ERR=55) (wn(l, k), l = 1, nn) read(fid, *, ERR=55) (sc(l, k), l = 1, nn) read(fid, *, ERR=55) ek(k) close(fid) if (k .gt. ms .or. k .eq. 0) call log_error("soil_read_input", & 'Soil type number > ms = max soil number, or 0', int=k) sn = to_string(int=k) call check_range(z(:, k), 'depth of soil '//trim(sn), (/0.001, 1e5/)) call check_range(cla(:, k) + sil(:, k) + san(:, k), 'texture of soil '//trim(sn), & (/99., 101./)) call check_range(bden(:, k), 'bulk density of soil '//trim(sn), (/0.001, 100./)) call check_range(poros(:, k), 'porosity of soil '//trim(sn), (/0.01, 100./)) call check_range(awc(:, k), 'water capacity of soil '//trim(sn), (/0.01, 100./)) call check_range(fc(:, k), 'field capacity of soil '//trim(sn), (/0.01, 100./)) call check_range(sc(:, k), 'saturated conductivity of soil '//trim(sn), (/0.01, 1e6/)) call check_range((/ek(k)/), "erodibility of soil "//trim(sn), (/0.01, 1./), index=k) !**** ESTIMATE CARBON FOR LOWER LAYERS IF AVAILABLE ONLY FOR UPPER LAYER if (cbn(3, k) .le. 0) then xx = z(2, k) do l = 3, nn dg = (z(l, k) - xx) if (cbn(l, k) .eq. 0.) cbn(l, k) = cbn(l - 1, k) * exp(- .001 * dg) xx = z(l, k) end do end if !**** CALCULATE INITIAL NUTRIENT CONTENT: wno3, wn, wmn, wpo, ap, pmn, op, hum ! INITIALIZE sc() if isc = 0 ! Units transformation: ! g/t --> kg/ha *wt1=bden()*dg/100., bden() - bulk density ! % --> g/t *10000. ! % --> kg/ha *wt1*10000. ! kg/ha --> g/t 1/wt1=100./dg/bden() ! mg/kg = g/t xx = 0. do l = 1, nn if (isc .eq. 1) sc(l, k) = 0. dg = (z(l, k) - xx) wt1 = bden(l, k) * dg / 100. if (wno3(l, k) .le. 0.) then wno3(l, k) = 10. * exp(- z(l, k) / 1000.) wno3(l, k) = wno3(l, k) * wt1 endif wn(l, k) = 10000. * wn(l, k) if (wn(l, k) .le. 0.) wn(l, k) = 10000. * cbn(l, k) / 10. wn(l, k) = wn(l, k) * wt1 wmn(l, k) = wn(l, k) * rtn wn(l, k) = wn(l, k) * (1. - rtn) wpo(l, k) = .125 * wn(l, k) if (ap(1, k) .le. 0.) ap(1, k) = 5. ap(l, k) = ap(1, k) * wt1 pmn(l, k) = ap(l, k) * (1. - psp) / psp op(l, k) = 4. * pmn(l, k) hum(l, k) = cbn(l, k) * wt1 * 172. xx = z(l, k) end do !**** CALCULATE PARTICLE SIZE DISTRIBUTION psz silt(1, k) = sil(1, k) / 100. clay(1, k) = cla(1, k) / 100. sand(1, k) = san(1, k) / 100. psz(1, k) = (1. - clay(1, k)) ** 2.49 * sand(1, k) psz(2, k) = .13 * silt(1, k) psz(3, k) = .20 * clay(1, k) if (clay(1, k) .le. .5) then if (clay(1, k) .lt. .25) then psz(4, k) = 2. * clay(1, k) else psz(4, k) = .28 * (clay(1, k) - .25) + .5 end if else psz(4, k) = .57 end if psz(5, k) = 1. - psz(1, k) - psz(2, k) - psz(3, k) - psz(4, k) !**** CALCULATE/RECALCULATE HYDRAULIC PARAMETERS wp, up, por, awc & ! CALCULATE SATURATED CONDUCTIVITY sc IN CASE isc=1 ! CORRECTION of up & por by silt content: Fred Hattermann 2002, ccf=.05 ! Removed 2011 from Tobias Vetter do l = 1, nn ! wp(l, k) = 0.4 * cla(l, k) * bden(l, k) / 100. wp(l, k) = fc(l, k) / 100 - awc(l, k) / 100 ! ccf = 0.05 ! up(l, k) = wp(l, k) + awc(l, k)/100. + ccf * sil(1, k) / 100. up(l, k) = fc(l, k) / 100 ! por(l, k) = 1. - bden(l, k) / 2.65 + ccf * sil(1, k) / 100. por(l, k) = poros(l, k) / 100 if (up(l, k) .ge. por(l, k)) then up(l, k) = por(l, k) - .05 wp(l, k) = up(l, k) - awc(l, k) / 100. if (wp(l, k) .le. 0.) then up(l, k) = por(l, k) * .75 wp(l, k) = por(l, k) * .25 end if end if awc(l, k) = up(l, k) - wp(l, k) if (isc .eq. 1) then asc(1) = 1. asc(2) = san(l, k) asc(3) = cla(l, k) asc(4) = poros(l, k) / 100. asc(5) = san(l, k) * san(l, k) asc(6) = cla(l, k) * cla(l, k) asc(7) = poros(l, k) * poros(l, k) / 10000. asc(8) = san(l, k) * cla(l, k) asc(9) = san(l, k) * poros(l, k) / 100. asc(10) = cla(l, k) * poros(l, k) / 100. asc(11) = asc(5) * cla(l, k) asc(12) = asc(6) * poros(l, k) / 100. asc(13) = asc(5) * poros(l, k) / 100. asc(14) = asc(6) * san(l, k) asc(15) = asc(7) * cla(l, k) asc(16) = asc(5) * asc(7) asc(17) = asc(6) * asc(7) do jj = 1, 17 sc(l, k) = sc(l, k) + asc(jj) * bsc(jj) end do sc(l, k) = exp(sc(l, k)) * 10. end if end do ek(k) = ek(k) * ekc0 end do ! loop over number of soil files !*********************************************************** END OF SOIL CYCLE !**** DETERMINE DRAINAGE COEFS ul, sumul, fc, sumfc, hk, wss & ! WRITE SELECTED SOIL PARAMETERS do k = 1, ms if (z(1, k) .gt. 0.) then nn = ns(k) bd(k) = 2.65 * (1. - por(1, k)) abd(k) = 0. zz = 0. xx = 0. swin(k) = 0. do l = 1, nn dg = z(l, k) - xx xz = por(l, k) * dg * 10. abd(k) = abd(k) + xz ul(l, k) = (por(l, k) - wp(l, k)) * dg sumul(k) = sumul(k) + ul(l, k) fc(l, k) = dg * (up(l, k) - wp(l, k)) sumfc(k) = sumfc(k) + fc(l, k) stin(l, k) = fc(l, k) * stinco hk(l, k) = - 2.655 / log10(fc(l, k) / ul(l, k)) xx = z(l, k) swin(k) = swin(k) + stin(l, k) end do abd(k) = abd(k) / (10. * z(nn, k)) end if end do !**** DEFINE CURVE NUMBERS ! cn2(k, n)- CN, conditions 2 do j = 1, mb do jea = 1, neap(j) n = mstruc(j, jea, 1) k = mstruc(j, jea, 2) select case(nsolgr(k)) case(1) cn2(j, jea) = cn2a(n) case(2) cn2(j, jea) = cn2b(n) case(3) cn2(j, jea) = cn2c(n) case(4) cn2(j, jea) = cn2d(n) end select end do end do cv = 0.01 call log_info('soil_read_input', 'Soil parameters read for N soils: ', int=k) return 55 continue call log_error('soil_read_input', "Failed to read soil file for soil", int=k) end subroutine soil_read_input subroutine soil_process(j, je, k, rain, swe) !**** PURPOSE: THIS SUBROUTINE DIVIDES EACH LAYER'S FLOW INTO 4 MM SLUGS ! and CALLS PERC to CALC PERCOLATION and LATERAL SUBSURFACE FLOW !**** CALLED IN: HYDROTOP !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! PARAMETERS & VARIABLES ! ! >>>>> COMMON PARAMETERS & VARIABLES ! fc(l, k) = field capacity, mm, calc readsol ! flate(j, je, l) = lateral subsurface flow, sum of prk for layer, mm ! ida = current day ! ipehd = number of hydrotope to print from perc(), if iperc = 1 ! iperc = switch code to print from perc() ! ipesb = number of subbasin to print from perc(), if iperc = 1 ! nn = number of soil layers, calc in subbasin ! poe(j, je, l) = percolation, sum of sep for layer, mm ! prk = lateral subsurface flow, calc in perc for 4mm slugs, mm ! rain = preinf(j, je)-qd, mm ! sep = percolation, calc in perc for 4mm slugs, recalc here, mm ! ste(j, je, l) = water storage in a layer, recalc here, mm ! su = water excess above FC, mm ! >>>>> ! >>>>> STATIC PARAMETERS ! add = local par ! amt = local par ! j1 = local par ! j2 = local par ! l = local par ! n1 = local par ! n2 = local par ! n3 = local par ! sum = local par ! tot = local par ! vvv = local par ! vvvs = local par ! >>>>> !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !**** Include common parameters real(dp), intent(in) :: rain real(dp), dimension(:, :), intent(in) :: swe integer j, je, k, j1, j2, l, n1, n2, n3 real(dp) river_route_add, sum, tot, vvv dimension vvv(11) sum = 0 n1 = 1 n3 = 1 !**** Calc vvv() n2 = 0 prk = 0. do l = 1, nn flate(j, je, l) = 0. poe(j, je, l) = 0. end do vvv(1) = ste(j, je, 1) - fc(1, k) if (vvv(1) .ge. 0.) then ste(j, je, 1) = fc(1, k) else vvv(1) = 0. end if vvv(1) = vvv(1) + rain if (vvv(1) .gt. 0.) then n1 = 1 n2 = 1 end if do l = 2, nn vvv(l) = ste(j, je, l) - fc(l, k) if (vvv(l) .gt. 0.) then ste(j, je, l) = fc(l, k) if (n2 .gt. 0) CYCLE n2 = 1 n1 = l else vvv(l) = 0. end if end do vvv(nn + 1) = 0. tot = 0. river_route_add = 0. !#### CALL PERC to calculate percolation and subsurface flow do while (n2 .ne. 0) n2 = 0 do j1 = n1, nn sum = 0. j2 = j1 + 1 if (j1 .eq. nn) j2 = nn if (vvv(j1) .gt. 0.) then if (n2 .le. 0) then n2 = 1 n3 = j1 end if su = vvv(j1) if (su .gt. amt) su = amt ste(j, je, j1) = ste(j, je, j1) + su sep = 0. ste(j, je, j1) = ste(j, je, j1) - sep sum = sep if (su .ne. sep) then if (ste(j, je, j1) .gt. fc(j1, k)) then call soil_percolation(j, je, k, j1, j2, swe) ste(j, je, j1) = ste(j, je, j1) - sep - prk sum = sum + sep river_route_add = river_route_add + prk flate(j, je, j1) = flate(j, je, j1) + prk poe(j, je, j1) = poe(j, je, j1) + sep endif end if vvv(j1) = vvv(j1) - su j2 = j1 + 1 vvv(j2) = vvv(j2) + sum end if end do tot = tot + sum n1 = n3 end do sep = tot prk = river_route_add return end subroutine soil_process !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& subroutine soil_percolation(j, je, k, j1, j2, swe) !**** PURPOSE: THIS SUBROUTINE COMPUTES PERCOLATION AND LATERAL SUBSURFACE ! FLOW FROM A SOIL LAYER WHEN FIELD CAPACITY IS EXCEEDED. !**** METHOD: SWRRB !**** CALLED IN: PURK !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! PARAMETERS & VARIABLES !inp j1 - layer, from purk !inp j2 - layer, from purk ! ! >>>>> COMMON PARAMETERS & VARIABLES ! hk(l, k) = beta coef. to estimate hydr. cond., used in perc ! hrtc(j, je, l) = return flow travel time, calc. hrflowtt, h ! hwss(2, j, je) = fun(field cap), calc in subbasin from hsumfc(j, jea) ! ida = current day ! ipehd = number of hydrotope to print from perc(), if iperc = 1 ! iperc = switch code to print from perc() ! ipesb = number of subbasin to print from perc(), if iperc = 1 ! prk = lateral subsurface flow, calc in perc for 4mm slugs, mm ! sc(l, k) = saturated conductivity, mm/h, calc, if isc = 1 ! sep = percolation, calc in perc for 4mm slugs, recalc here, mm ! ste(j, je, l) = water storage in a layer, recalc here, mm ! su = water excess above FC ! swe(j, je) = soil water content, mm ! te(j, je, l) = daily average temp at the bottom of each layer, degree C ! ul(l, k) = upper limit water content in layer, mm ! >>>>> ! >>>>> STATIC PARAMETERS ! adjf = local par ! cr = local par ! fx = local par ! iconfig = local par ! rtf = local par ! rtw = local par ! stu = local par ! stu1 = local par ! stu2 = local par ! stz = local par ! sup = local par ! sup1 = local par ! xx = local par ! zz = local par ! >>>>> !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !**** Include common parameters real(dp), dimension(:, :), intent(in) :: swe integer j, je, k, j1, j2 real(dp) cr, fx, rtf, rtw, stu, stu1, stu2, stz, sup, sup1, xx, zz !**** IF TEMP OF LAYER IS BELOW 0 C - NO WATER FLOW if (te(j, je, j1) .le. 0.) then sep = 0. prk = 0. rtf = 0. sup1 = 0. rtw = 0. else sup = su - sep sup1 = sup !**** CALC SEEPAGE TO NEXT LAYER xx = hwss(1, j, je) - hwss(2, j, je) * swe(j, je) if (xx .lt. - 20.) xx = - 20. rtw = 10. * (1. - swe(j, je) / (swe(j, je) + exp(xx))) rtf = hrtc(j, je, j1) * rtw if (rtf .gt. 0.) rtf = 1. - exp(- 1. / rtf) stz = ste(j, je, j1) / ul(j1, k) if (stz .ge. 1.) then fx = 1. else fx = stz ** hk(j1, k) end if stu = ste(j, je, j2) / ul(j2, k) stu1 = ul(j1, k) / ste(j, je, j1) stu2 = ste(j, je, j1) / ul(j1, k) if (stu .ge. 1.) then sep = 0. prk = rtf * sup else cr = sqrt(1. - stu) ! SL *** zz = 24. * fx * cr * sc(j1, k) * sccor(j) / sup if (zz .gt. 10.) then sep = sup prk = 0. else sep = sup * (1. - exp(- zz)) sup = sup - sep prk = rtf * sup endif end if end if !**** Write if (j .eq. 1 .and. je .eq. 1) then if (sep <= 0.) sep = .000001 end if return end subroutine soil_percolation subroutine soil_temperature(zz, j, je, k, bcv, ida, mo, pit, preinf, swe, tmax_tmp, tmin_tmp, tx_tmp) !**** PURPOSE: THIS SUBROUTINE ESTIMATES DAILY AVERAGE TEMPERATURE ! AT THE BOTTOM OF EACH SOIL LAYER ! =f(air temp, precip, residue, snow cover) !**** CALLED IN: HYDROTOP !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! PARAMETERS & VARIABLES ! ! >>>>> COMMON PARAMETERS & VARIABLES ! abd(k) = sum of porosity for soil type, mm ! amp(j) = annual 1/2 amplitude in daily average temp, calc readwet, ! degree C ! avt(j) = average annual air temperature, calc readwet, degree C ! bcv(j, je) = lag factor for residue and snow effects on surf.temp. ! ida = current day ! isolt = switch code to print from solt() ! isosb = number of subbasin to print from solt(), if isolt = 1 ! iy = current year ! mo = current month ! nn = number of soil layers, calc in subbasin, cycle 100 !block pit = 58.13 ! precip = precipitation, mm, calc in readcli ! swe(j, je) = soil water, mm, calc in hydrotop (previous day) ! te(j, je, l) = daily ave temp at the bottom of each layer, degree C ! te0(j, je) = bare soil surface temp, degree C ! tmn(j) = daily min temp. for subbasin, readcli, degree C ! tmx(j) = daily max temp. for subbasin, readcli, degree C ! tx(j) = daily aver. temperature, readcli, degree C ! wft(m, j) = monthly prob. of rainy day, calc readwet ! z(l, k) = soil depth, read in readsol, mm ! >>>>> ! >>>>> STATIC PARAMETERS ! alx = local par ! b = local par ! dd = local par ! mdp = max damping depth for the soil ! dt = local par ! f = local par ! l = local par ! st0 = bare soil surfac temp ! ta = local par ! wc = local par ! ww = local par ! xi = local par ! xx = local par ! zd = local par ! >>>>> !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !**** Include common parameters real(dp), dimension(:, :), intent(in) :: bcv integer, intent(in) :: ida integer, intent(in) :: mo real(dp), intent(in) :: pit real(dp), dimension(:, :), intent(in) :: preinf real(dp), dimension(:, :), intent(in) :: swe real(dp), intent(in) :: tmax_tmp, tmin_tmp, tx_tmp integer j, je, k, l real(dp) zz, alx, b, dd, mdp, dt, f, st0, ta, wc, ww, xi, xx, zd xi = ida alx = (xi - 200.) / pit f = abd(k) / (abd(k) + 686. * exp(- 5.63 * abd(k))) mdp = 1000. + 2500. * f ww = .356 - .144 * abd(k) b =log(500. / mdp) wc = swe(j, je) / (ww * z(nn, k)) f = exp(b * ((1. - wc) / (1. + wc)) ** 2.) dd = f * mdp !**** CALC st0 - bare soil surfac temp if (preinf(j, je) .le. 0.) then st0 = wft(mo, j) * (tmax_tmp - tx_tmp) + tx_tmp else st0 = wft(mo, j) * (tx_tmp - tmin_tmp) + tmin_tmp end if ta = avt(j) + amp(j) * cos(alx) xx = bcv(j, je) * te0(j, je) + (1. - bcv(j, je)) * st0 dt = xx - ta te0(j, je) = st0 zz = 2. * dd xx = 0. !**** CALC te() do l = 1, nn zd = - z(l, k) / dd te(j, je, l) = avt(j) + (amp(j) * cos(alx + zd) + dt) * exp(zd) end do return end subroutine soil_temperature end module soil