! FILE subbasin.f ! ! SUBROUTINES IN THIS FILE CALLED FROM ! subroutine subbasin(icode, ihout, inum1, inum2) main ! subroutine hrflowtt(j, je, k) subbasin module subbasin use utilities, only : dp, path_max_length, log_error implicit none real(dp), parameter :: clt = 57.296 ! Inflow Hydrograph integer, save, dimension(:), allocatable :: inum2s ! total number of subbasins integer, save :: mb ! mhyd = max number of hydrograph nodes integer, save :: mhyd ! mch = max number of channels integer, save :: mch ! array containing subbasin list integer, save, dimension(:), allocatable :: subs ! number of HRUs in subbasin integer, save, dimension(:), allocatable :: neap ! Hyd. Storage Location integer, save, dimension(:), allocatable :: ihouts ! Reach No. integer, save, dimension(:), allocatable :: inum1s ! Subbasin-level driving data file character(len=path_max_length) :: climate_input_file = "climate.csv" ! 1000. * da * flu(j), calc in subbasin real(dp), save :: aff ! precipitation, mm, read in readcli real(dp), save :: precip ! weighted aver. Cur.Num. in subbasin (SUM(cn)) real(dp), save :: sumcn ! weighted aver. surface flow in subbasin, SUM(qd), mm real(dp), save :: xqi ! weighted aver. subsurface flow in subbasin, SUM(ssf), mm real(dp), save :: xssf ! subbasin percolation, mm real(dp), save :: xsep ! weighted aver. pot. evapotr. in subbasin, SUM(eo), mm real(dp), save :: xeo ! weighted aver. soil wat. index in subbasin, SUM(swimd) real(dp), save :: xswind ! weighted aver. N loss with qd in subbasin, SUM(yno3), kg / ha real(dp), save :: xyno3 ! weighted aver. N loss with ssf in subbasin, SUM(ssfn), mm real(dp), save :: xssfn ! weighted aver. N leaching to g - w in subb., kg / ha real(dp), save :: xpercn ! transmission losses in subbasin, mm real(dp), save :: qtl ! water yield in subbasin = xqi + xssf + gwq(j) - qtl, mm real(dp), save :: wysb ! real(dp), save, dimension(30) :: sub = 0. character(len=path_max_length), dimension(:), allocatable :: station_ids ! Subbasin ids read from the subbasin.csv integer, dimension(:), allocatable :: subbasin_ids ! array storing an ID for each subbasin indicating the catchment ID integer, save, dimension(:), allocatable :: subcatch_id ! indeces of catchment ids in the catchment_id array integer, save, dimension(:), allocatable :: subcatch_idx ! water discharge in the basin outlet, m3/sec real(dp), dimension(:, :), allocatable :: obs_discharge ! subbasin numbers of observed discharge integer, save, dimension(100) :: obssb ! number of observed discharge columns, default 1 (only outlet) integer, save :: nqobs = 1 logical, save :: bRunoffDat = .false. character(len=path_max_length) :: discharge_input_file = "discharge.csv" integer :: discharge_input_file_id ! Counter for offset in discharge.csv column, as leap year logic etc is in outer loops integer :: discharge_offset = 1 ! fraction of subbasin area in the basin real(dp), save, dimension(:), allocatable :: flu ! SUM(xqd) - surface runoff for subbasin, mm real(dp), save, dimension(:), allocatable :: smq ! SUM of subsurface runoff for subbasin real(dp), save, dimension(:), allocatable :: smsq ! monthly SUM of precipitation in subbasin real(dp), save, dimension(:), allocatable :: sbp ! average daily temp., degree C real(dp), save, dimension(:), allocatable :: tx ! daily precipitation in the subbasin, mm real(dp), save, dimension(:), allocatable :: subp ! monthly SUBBASIN outputs (dif. components) sysub (see subbasin.f) real(dp), save, dimension(:, :), allocatable :: susb ! subbasin outputs for routing, andlogue to varoute() real(dp), save, dimension(:, :), allocatable :: sda ! code to switch between routing subroutines integer, save, dimension(:), allocatable :: icodes ! drainage area for subbasin, km2 real(dp), save, dimension(:), allocatable :: dart ! hydrograph storage location for each subbasin integer, save, dimension(:), allocatable :: subouthyd ! weighted aver. soluble P leaching in subbasin, kg / ha real(dp), save :: xysp ! daily soil loss from subbasin caused by water erosion, t real(dp), save :: yd ! real(dp), save, dimension(366) :: runs ! slope length (m) for time of concentration and erosion USLE slope length real(dp), save, dimension(:), allocatable :: slope_length ! average slope steepness for subbasin, m / m real(dp), save, dimension(:), allocatable :: stp ! subbasin area [m2] real(dp), save, dimension(:), allocatable :: sbar ! only used when CaMaFLood is activated real(dp), save, dimension(:), allocatable :: runoff_mm ! weighted aver. snow evap in subbasin, SUM(snoev), mm real(dp), save :: xsnoev ! weighted aver. org.P in layer 1 in subb., SUM(porg), kg / ha real(dp), save :: xporgp ! weighted aver. total P (porg + pms + pma) in subb., g / t real(dp), save :: xpsed ! climate.csv column positions integer :: tmean_column_ix, tmin_column_ix, tmax_column_ix, & precipitation_column_ix, radiation_column_ix, humidity_column_ix = 0 integer :: subbasin_climate_file_id integer :: river_runoff_output_id = 0 character(len=path_max_length) :: subbasin_input_file = "subbasin.csv" ! (12, :) average monthly max temp, degree C real(dp), save, dimension(12) :: & obmx = (/31.44, 32.71, 34.51, 34.73, 33.93, 31.46, 28.66, 29.07, 30.76, 30.56, 29.98, 30.27/) ! (12, :) average monthly min temp, degree C real(dp), save, dimension(12) :: & obmn = (/7.51, 8.72, 10.34, 11.32, 11.82, 10.66, 10.43, 10.64, 9.41, 10.52, 9.34, 7.89/) ! (12, :) monthly mean event of daily rainfall, mm real(dp), save, dimension(12) :: & rst = (/0.2, 0.2, 0.9, 1.3, 4.2, 7.6, 10.7, 10.3, 7.4, 4.3, 1.3, 0.3/) ! monthly max .5h rain for period of record, mm real(dp), save, dimension(12) :: & wim = (/7.0, 7.0, 9.0, 9.0, 10.0, 10.0, 10.0, 10.0, 9.0, 9.0, 7.0, 7.0/) ! wheather generator data from wgen.dat real(dp), save, dimension(:), allocatable :: tp5, tp6, tpnyr ! default values for the weather generator variables real(dp) :: & tp5_default = 15., & tp6_default = 22., & tpnyr_default = 54. ! prw(1, m, j) = monthly probability of wet day after dry day real(dp), save, dimension(2, 12) :: & prw = reshape((/0.03, 0.59, 0.06, 0.43, 0.09, 0.48, 0.14, 0.57, 0.21, 0.56, 0.32, & 0.5, 0.39, 0.41, 0.43, 0.51, 0.31, 0.42, 0.12, 0.49, 0.08, 0.56, 0.05, 0.63/), & shape=(/2, 12/)) integer :: subbasin_input_file_id integer :: subbasin_routing_input_file_id = 0 character(len=path_max_length) :: subbasin_routing_input_file = "subbasin_routing.csv" namelist / SUBBASIN_PARAMETERS / & subbasin_routing_input_file, & climate_input_file, & bRunoffDat, & discharge_input_file, & tp5_default, & tp6_default, & tpnyr_default, & obmx, & obmn, & wim, & prw, & rst contains subroutine subbasin_initialise use input, only : & header_column_index, & input_open_file, & input_count_rows, & get_config_fid use output, only: output_register_subbasin_var use utilities, only: hydrograph_storage_location #ifndef with_netcdf integer sfid #endif integer i read(get_config_fid(), SUBBASIN_PARAMETERS) !**** OPEN subbasin-level driving data (if not using NetCDF) #ifndef with_netcdf subbasin_climate_file_id = input_open_file(climate_input_file) sfid = subbasin_climate_file_id tmean_column_ix = header_column_index('tmean', sfid, required=.True.) tmin_column_ix = header_column_index('tmin', sfid, required=.True.) tmax_column_ix = header_column_index('tmax', sfid, required=.True.) precipitation_column_ix = header_column_index('precipitation', sfid, required=.True.) ! Optional columns radiation_column_ix = header_column_index('radiation', sfid) humidity_column_ix = header_column_index('humidity', sfid) #endif if (bRunoffDat) then discharge_input_file_id = input_open_file(discharge_input_file) obssb(:) = 1 end if river_runoff_output_id = output_register_subbasin_var("river_runoff") subbasin_input_file_id = input_open_file(subbasin_input_file) mb = input_count_rows(subbasin_input_file_id, .true.) mch = mb allocate(neap(mb)) ! number of HRUs in subbasin neap = 0 allocate(sbar(mb)) ! subbasin area [m2] sbar = 0. call subbasin_read_routing call subbasin_allocate call subbasin_read_input if (bRunoffDat) call subbasin_allocate_obs_discharge ! subbasin hydrograph storage location for routed Q output do i = 1, mb subs(i) = i end do ! Identify hydrograph storage location for each subbasin subouthyd = hydrograph_storage_location(icodes, ihouts, inum1s, mhyd, mb, subs) ! array subs not used thereafter... deallocate(subs) end subroutine subbasin_initialise subroutine subbasin_read_input use input, only : read_integer_column, read_string_column, read_real_column call read_integer_column(subbasin_input_file_id, "subbasin_id", subbasin_ids) call read_string_column(subbasin_input_file_id, "station_id", station_ids, "") call read_integer_column(subbasin_input_file_id, "catchment_id", subcatch_id, 1) call read_real_column(subbasin_input_file_id, "sl", slope_length) call read_real_column(subbasin_input_file_id, "stp", stp, range=(/0., 3./), closed="n") end subroutine subbasin_read_input subroutine subbasin_read_routing use input, only : & read_integer_column, & read_string_column, & input_open_file, & input_count_rows subbasin_routing_input_file_id = input_open_file(subbasin_routing_input_file) mhyd = input_count_rows(subbasin_routing_input_file_id, .true.) ! Allocate routing node variables allocate(icodes(mhyd)) allocate(ihouts(mhyd)) allocate(inum1s(mhyd)) allocate(inum2s(mhyd)) icodes = 0 ihouts = 0 inum1s = 0 inum2s = 0 ! **** READ 3 - routing structure file & count number of hydrograph nodes call read_integer_column(subbasin_routing_input_file_id, "icd", icodes) call read_integer_column(subbasin_routing_input_file_id, "iht", ihouts) call read_integer_column(subbasin_routing_input_file_id, "inm1", inum1s) call read_integer_column(subbasin_routing_input_file_id, "inm2", inum2s) end subroutine subbasin_read_routing subroutine subbasin_allocate allocate(dart(mhyd)) allocate(flu(mb)) allocate(runoff_mm(mb)) allocate(sbp(mb)) allocate(sda(10, mb)) allocate(smq(mb)) allocate(smsq(mb)) allocate(station_ids(mb)) allocate(subbasin_ids(mb)) allocate(subouthyd(mb)) allocate(subp(mb)) allocate(subs(mb)) allocate(susb(30, mb)) allocate(tx(mb)) allocate(slope_length(mb)) allocate(tp5(mb)) ! 10 year frequency of .5h rainfall (mm) allocate(tp6(mb)) ! 10 year frequency of .6h rainfall (mm) allocate(tpnyr(mb)) ! number of years of record max .5h rainfall allocate(subcatch_id(mb)) allocate(subcatch_idx(mb)) allocate(stp(mb)) stp = 0. subcatch_id = 1 subcatch_idx = 1 tp5 = 0. tp6 = 0. tpnyr = 0. slope_length = 0. dart = 0. flu = 0. precip = 0. runoff_mm = 0. sbp = 0. sda = 0. smq = 0. smsq = 0. subp = 0. susb = 0. tx = 0. susb = 0. subouthyd = 0 subs = 0 tp5 = tp5_default tp6 = tp6_default tpnyr = tpnyr_default end subroutine subbasin_allocate subroutine subbasin_close deallocate(dart) deallocate(flu) deallocate(icodes) deallocate(ihouts) deallocate(inum1s) deallocate(inum2s) deallocate(neap) deallocate(sbar) deallocate(sbp) deallocate(sda) deallocate(smq) deallocate(smsq) deallocate(station_ids) deallocate(subbasin_ids) deallocate(subp) deallocate(susb) deallocate(tx) deallocate(subcatch_id) deallocate(subcatch_idx) deallocate(tp6) deallocate(tpnyr) deallocate(slope_length) deallocate(stp) close(discharge_input_file_id) end subroutine subbasin_close subroutine subbasin_initialise_area(da) real(dp), intent(in) :: da integer id !**** COMPUTE DRAINAGE AREA FOR EACH SUBBASIN do id = 1, mhyd if (icodes(id) .eq. 0) exit if (icodes(id) .eq. 1) then dart(ihouts(id)) = da * flu(inum1s(id)) end if if (icodes(id) .eq. 2) then dart(ihouts(id)) = dart(inum2s(id)) end if if (icodes(id) .eq. 3) then dart(ihouts(id)) = dart(inum2s(id)) end if if (icodes(id) .eq. 5) then dart(ihouts(id)) = dart(inum1s(id)) + dart(inum2s(id)) end if if (icodes(id) .eq. 6) then dart(ihouts(id)) = 0. end if end do end subroutine subbasin_initialise_area subroutine subbasin_initialise_subbasin(canev, sml, xcklsp, xet, xnorg, xnorgp, xporg, xpsedp, xqd, yon, yph) !**** PURPOSE: this subroutine initialises subbasin variables !**** CALLED IN: subbasin !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ real(dp), intent(out) :: canev real(dp), intent(out) :: sml real(dp), intent(out) :: xcklsp real(dp), intent(out) :: xet real(dp), intent(out) :: xnorg real(dp), intent(out) :: xnorgp real(dp), intent(out) :: xporg real(dp), intent(out) :: xpsedp real(dp), intent(out) :: xqd real(dp), intent(out) :: yon real(dp), intent(out) :: yph sumcn = 0. xqd = 0. xqi = 0. xssf = 0. qtl = 0. yd = 0. yon = 0. yph = 0. xsep = 0. xeo = 0. xet = 0. xsnoev = 0. xswind = 0. xyno3 = 0. xysp = 0. xssfn = 0. xpercn = 0. xnorg = 0. xnorgp = 0. xporg = 0. xporgp = 0. xpsed = 0. xpsedp = 0. xcklsp = 0. sml = 0. canev = 0. end subroutine subbasin_initialise_subbasin subroutine subbasin_read_climate(humi, mb, ra, subp, tmn, tmx, tx) !**** PURPOSE: this subroutine reads climate data ! & hydrological data from runoff.dat !**** CALLED IN: MAIN !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! PARAMETERS & VARIABLES ! ! >>>>> COMMON PARAMETERS & VARIABLES ! daylen(j) = day length in subbasin ! humi(j) = air humidity in the subbasin, % ! ida = current day in the year (Julian day) ! mb = number of subbasins !block pit = parameter for estimation of the day length ! precip = precipitation in the current subbasin, mm ! ra(j) = radiation, J/cm^2 ! subp(j) = daily precipitation in the subbasin, mm ! tmn(j) = daily min temp in the subbasin, degree C ! tmx(j) = daily max temp in the subbasin, degree C ! tx(j) = daily aver temp in the subbasin, degree C ! ylc(j) = cos(lat()/clt), lat() - latitude, clt=57.296, ! used to calc rmx in evap ! yls(j) = sin(lat()/clt), lat() - latitude, clt=57.296, ! used to calc rmx in evap ! >>>>> ! >>>>> STATIC PARAMETERS ! ch = parameter for estim. day length ! h = parameter for estim. day length ! id = day ! iiyr = year ! imn = month ! k = count parameter ! sd = parameter for estim. day length ! xi = day (real(dp) number) ! >>>>> !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ use input, only : read_real_column real(dp), dimension(:), intent(inout) :: humi integer, intent(inout) :: mb real(dp), dimension(:), intent(inout) :: ra real(dp), dimension(:), intent(inout) :: subp real(dp), dimension(:), intent(inout) :: tmn real(dp), dimension(:), intent(inout) :: tmx real(dp), dimension(:), intent(inout) :: tx integer k, fid real(dp), dimension(mb) :: mx, mn !**** READ CLIMATE DATA & MEASURED DISCHARGE (current day) fid = subbasin_climate_file_id call read_real_column(fid, array=tx, index=tmean_column_ix) ! no skip call read_real_column(fid, array=tmn, index=tmin_column_ix, skip=-mb) call read_real_column(fid, array=tmx, index=tmax_column_ix, skip=-mb) call read_real_column(fid, array=subp, index=precipitation_column_ix, skip=-mb) if (radiation_column_ix > 0) then call read_real_column(fid, array=ra, index=radiation_column_ix, skip=-mb) else ra = 1. end if if (humidity_column_ix > 0) then call read_real_column(fid, array=humi, index=humidity_column_ix, skip=-mb) else humi = - 999.9 end if ! In order to avoid errors if Tmax < Tmin the following conversion is done do k = 1, mb mx(k) = max(tmx(k), tmn(k)) mn(k) = min(tmx(k), tmn(k)) tmx(k) = mx(k) tmn(k) = mn(k) end do end subroutine subbasin_read_climate subroutine subbasin_allocate_obs_discharge use input, only : has_column integer k nqobs = 0 do k = 1, size(station_ids) if (len_trim(station_ids(k)) > 0 .and. \ has_column(station_ids(k), discharge_input_file_id)) then nqobs = nqobs + 1 end if end do allocate(obs_discharge(366, nqobs)) end subroutine subbasin_allocate_obs_discharge subroutine runsubbasin(ihout, inum1, bSubcatch, da, da9, daycounter, ida, ieap, iy, iyr, mo, nbyr, nd) !**** PURPOSE: THIS SUBROUTINE PERFORMS SUBBASIN OPERATIONS: ! 1. ESTABLISHES INITIAL CONDITIONS IN HYDROTOPES, ! 2. CALLS HYDROTOP TO CALCULATE HYDROTOP PROCESSES, ! 3. AGGREGATES HYDROTOP OUTPUTS, ! 4. CALCULATES EROSION IN SUBBASIN, ! 5. CALCULATES RETENTION OF NUTRIENTS, ! 6. CALCULATES VARIABLES FOR OUTPUT, and ! 7. SETS OUTPUTS FROM SUBBASIN FOR ROUTING !**** CALLED IN: MAIN ! icode = 1, ihout = inum1 = inum2 = j !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! PARAMETERS & VARIABLES ! ! >>>>> INDICES ! i, ida = DAY ! j = SUBBASIN NR ! jea, je = HYDROTOPE NR ! k = SOIL NR ! n = LAND USE NR ! l = SOIL LAYER NR ! icr, iv = CROP or VEGETATION NR ! ih = HYDR. STORAGE LOCATON NR ! io = CROP OPERATION NR ! if = CROP FERTILISATION NR ! >>>>> ! ! >>>>> COMMON PARAMETERS & VARIABLES ! af = 1000. * da, km2 ! aff = af * flu(j) = 1000 * da * flu(j), unit transf. coef. ! alai(j, je) = leaf area index ! alai0(n) = initial leaf area index ! ano3(j, je, l) = wno3(l, k) N-NO3 in sub j, HRU je, layer l, kg/ha ! anora(j, je, l) = wmn(l, k) active organic N in sub j, HRU je, layer l, ! kg/ha ! anors(j, je, l) = wn(l, k) stable organic N in sub j, HRU je, layer l, ! kg/ha ! ap(l, k) = labile (soluble) phosphorus, kg/ha (initialisation) ! bd(k) = bulk density of the upper (!) soil layer, ! calc from porosity, g/cm3 ! canev = canopy evaporation, mm ! cklsp(j, je) = combined c, k, ls, p factor ! cn = current CN in hydrotope ! css(j) = sediment conc in return flow, ppm ! cv(j) = initial land cover, kg/ha ! cva(j, je) = vegetation cover, kg/ha ! da = basin area, km2 ! dart(j) = drainage area for subbasin, km2 ! degNgrw = N degradation in groundwater, (1/day) ! degNsub = N degradation in subsurface flow, (1/day) ! degNsur = N degradation in surface flow, (1/day) ! degPsur = P degradation in surface flow, (1/day) ! dm(j, je) = total biomass, kg/ha ! dur = flow duration, h ! eo = potential evapotranspiration, mm ! eopot = potential evapotranspiration, mm ! ep = plant transpirationm, mm ! es = soil evaporation, mm ! evasum(j, je) = SUM(ep + es + canev), mm ! fc(l, k) = field capacity, mm ! flu(j) = fraction of subbasin area in the basin ! frar(j, je) = fractional area of HYDROTOPE in SUBBASIN ! g(j, je) = fraction of heat unit accumulated ! gwchrg = groundwater recharge, mm ! gwq(j) = groundwaterw contribution to stream, mm ! gwseep = groundwater seepage, mm ! 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 ! hwss(2, j, je) = fun(field cap), calc i subbasin from hsumfc(j, je) ! ida = current day ! gis_m = switch code to write monthly results (water & crops) for GIS output ! gis_y = switch code to write annual results (water & crops) for GIS output ! gis_ave = switch code to write average annual results (water & crops) for GIS output ! ih1, ih2, ih3, ih4 = hydrotopes for HYDROTOPE PRINTOUT ! isb1, isb2, isb3, isb4 = subbasins for HYDROTOPE PRINTOUT ! isu1, isu2, isu3 = subbasins for SUBBASIN PRINTOUT ! iv = index for channel parameters chl(, ), chw(, ), chk(, ) ! iy = current year as counter (1, ..., nbyr) ! iyr = current year ! mstruc(j, je, 2) = HRU structure vector to define k, n ! ndpri = day to write crop yield for GIS output ! neap(j) = number of HYDROTOPEs in subbasin, calc readbas ! nn = number of soil layers ! ns(k) = number of soil layers ! nsa(k) = number of layers for arable soils ! nucr(j, je) = crop number (database) ! op(l, k) = stable min. P, estimated from pmn(), kg/ha ! parsz(5, j) = particle size distribution, calc in subbasin ! percn = N leaching to g-w, kg/h ! plab(j, je, l) = ap(l, k) labile (soluable) P in sub j, HRU je, ! layer l, kg/ha ! pma(j, je, l) = pmn(l, k) active mineral P in sub j, HRU je, ! layer l, kg/ha ! pmn(l, k) = act. min. P, estimated in readsol, kg/ha ! pms(j, je, l) = op(l, k) stable mineral P in sub j, HRU je, ! layer l, kg/ha ! porg(j, je, l) = wpo(l, k) organic P in sub j, HRU je, layer l, kg/ha ! pr = peak runoff rate, m3/sec. ! gwrsum(j, je) = annual sum of percolation for GIS output, SUM(sep), mm ! precip = precipitation in subbasin, mm ! presum(j, je) = annual sum of precip for GIS output, mm ! psz(5, k) = particle size distribution ! q1 = daily surface runoff = xqd before call tran, mm ! qd = daily surface runoff in hydrotope, mm ! qi = daily surface runoff in hydrotope, mm ! qtl = transmission losses in subbasin, mm ! ra(j) = radiation in J/cm^2 ! rd(j, je) = root depth, mm ! retNgrw = N retention in groundwater, days ! retNsub = N retention in subsurface flow, days ! retNsur = N retention in surface flow, days ! retPsur = P retention in surface flow, days ! revap = revaporation from groundwater, mm ! revapst(j) = transmission losses in the reach, calc in route, m3 ! rsd(j, je, 2) = crop residue, kg/ha ! runsum(j, je) = annual sum of surface+subsurface runoff for GIS output, ! SUM(qd+ssf), mm ! sbp(j) = precipitation in subbasin, mm ! sda(1:10, j) = subbasin outputs for routing, analogue to varoute() ! varoute(1:10, ihout) = sda(1:10, j) ! sep = percolation to g-w, mm ! sml = snow melt, calc in snom, mm ! smq(j) = SUM(xqd) - surface runoff for subbasin, mm ! smsq(j) = SUM(xssf) - subsurface runoff for subbasin, mm ! sno(j) = accumulated value of snow cover, mm ! snoa(j, je) = snow water content in sub j, HRU je, mm ! snoev = snow evaporation, mm ! snow = precipitation as snow, mm ! ssf = subsurface flow in HYDROTOPE, mm ! ssfn = N loss with subsurface flow, kg/ha ! ste(j, je, l) = water content in sub j, HRU je, layer l, mm ! stin(l, k) = soil water content in soil layer, mm ! sub(1:28) = basin outputs (weighted sums) daily, see desc. below ! sub(1) = SUM(precip) precipitation, mm ! sub(2) = SUM(precip), if tx()<0 snowfall, mm ! sub(3) = SUM(sml) snowmelt(mm) ! sub(4) = SUM((tmx*flu) average "weighted" max temp., degree C ! sub(5) = SUM(tmn*flu) average "weighted" min temp., degree C ! sub(6) = SUM(ra*flu) radiation, J/cm^2 ! sub(7) = SUM(sumcn*flu) curve number ! sub(8) = SUM(xqi*flu) surface runoff, mm ! sub(9) = SUM(xssf*flu) subsurface runoff, mm ! sub(10) = SUM(qtl*flu) channel transm. losses, mm ! sub(11) = SUM(xsep*flu) percolation, mm, ! sub(12) = SUM(xeo*flu) pot evapotr., mm ! sub(13) = SUM(xet*flu) evapotranspiration, mm ! sub(14) = SUM(snoev*flu) snow evaporation, mm ! sub(15) = SUM(gwq*flu) gr. water flow, mm ! sub(16) = SUM(revap*flu) revap from g-w to soil prof., mm ! sub(17) = SUM(gwseep*flu) g-w percolation, mm ! sub(18) = SUM(gwchrg*flu) g-w recharge, mm ! sub(19) = SUM(xswimd*flu) soil water index for basin ! sub(20) = SUM(wysb*flu) wysb=qi+ssf+gwq(j)-qtl, wat yld, mm ! sub(21) = SUM(yd/(100*da*flu) sed yield, t/ha ! sub(22) = SUM(yon*flu) org. N loss with sed, kg/ha ! sub(23) = SUM(yph*flu) org. P loss with sed, kg/ha ! sub(24) = SUM(ysp*flu) soluble P, kg/ha ! sub(25) = SUM(xyno3*flu) NO3-N in surface runoff, kg/ha ! sub(26) = SUM(xssfn*flu) NO3-N in subsurface runoff, kg/ha ! sub(27) = SUM(xpercn*flu) NO3-N leached to g-w, kg/ha ! sub(28) = SUM(xuno3*flu) NO3-N uptake by plants , kg/ha ! subp(j) = daily precipitation in subbasin, mm (from readcli) ! sumcn = weighted aver. Cur.Num. in subbasin (SUM(cn)) ! susb(30, j) = subbasin outputs daily, analogue to sub() above ! swe(j, je) = water content in sub j, HRU je ! swin(k) = soil water content in soil profile, mm ! swind = soil water index for hydrotope: ! swind=swe(j, je)/sumfc(k), 1/1 ! sym(j) = SUM(yd) = sum of sediment yield for subbasins, t ! tmn(j) = daily min temp in the subbasin, degree C ! tmx(j) = daily max temp in the subbasin, degree C ! tsnfall = threshold temperature for snow fall ! tx(j) = daily aver temp in the subbasin, degree C ! ul(l, k) = upper limit water content in layer, mm ! varoute(*, ih) = subbasin outputs for routing: ! Name Units Definition ! varoute(2, ih) |(m^3) |surface flow ! varoute(3, ih) |(tons) |sediment ! varoute(4, ih) |(kg) |organic N ! varoute(5, ih) |(kg) |organic P ! varoute(6, ih) |(kg) |nitrate N ! varoute(7, ih) |(kg) |soluble P ! varoute(8, ih) |(m^3) |subsurface + g-w flow ! vo = surface runoff, m3 ! wmn(l, k) = active org. N content, kg/ha ! wn(l, k) = stable organic N content, kg/ha ! wno3(l, k) = NO3-N content, kg/ha ! wpo(l, k) = organic P content, kg/ha ! wysb = water yield in subbasin = xqi+xssf+gwq(j)-qtl, mm ! xcklsp = weighted aver. comb C*K*LS*P in subbasin, SUM(cklsp) ! xeo = weighted aver. pot. evapotr. in subbasin, SUM(eo), mm ! xet = weighted aver. act. evapotr. in subbasin, ! SUM(ep+es), mm ! xnorg = weighted aver. st.org.N in lay 1 in subb., ! SUM(anors), g/t ! xnorgp = weighted aver. st.org.N in lay 1 in subb., ! SUM(anors), kg/ha ! xpercn = weighted aver. N leaching to g-w in subb., kg/ha ! xporg = weighted aver. org.P in layer 1 in subb., ! SUM(porg), g/t ! xporgp = weighted aver. org.P in layer 1 in subb., ! SUM(porg), kg/ha ! xpsed = weighted aver. total P (porg+pms+pma) in subb., g/t ! xpsedp = weighted aver. total P (porg+pms+pma) in subb., kg/ha ! xqd = weighted aver. surface flow in subbasin, SUM(qd), mm ! xqi = weighted aver. surface flow in subbasin, SUM(qd), mm ! xsep = weighted aver. percolation in subbasin, SUM(sep), mm ! xsnoev = weighted aver. snow evap in subbasin, SUM(snoev), mm ! xssf = weighted aver. subsurface flow in subbasin, ! SUM(ssf), mm ! xssfn = weighted aver. N loss with ssf in subbasin, ! SUM(ssfn), mm ! xswind = weighted aver. soil wat. index in subbasin, SUM(swimd) ! xwysb = water yield for subbasin, m3 ! xxswind = weighted aver. soil wat. index in the basin ! xyno3 = weighted aver. N loss with qd in subbasin, ! SUM(yno3), kg/ha ! xysp = weighted aver. soluble P leaching in subbasin, kg/ha ! yd = daily soil loss caused by water erosion, t ! yno3 = N loss with surface flow, kg/ha ! yon = org N loss with erosion, kg/ha ! yph = org P loss with erosion, kg/ha ! ysp = soluble P leaching, kg/ha ! >>>>> ! >>>>> STATIC PARAMETERS ! ii = local par ! im = local par ! j = subbasin nr ! jea = hydrotope nr ! jj = local par ! k = soil nr ! l = soil layer nr ! n = land use nr ! wabad = daily water balance ! >>>>> ! ! >>>>> UNITS TRANSFORMATION: ! g/t --> kg/ha *wt1=bden*dg/100., bden=bulk dens., dg=layer thickness ! % --> g/t *10000. ! % --> kg/ha *wt1*10000. ! kg/ha --> g/t 1/wt1 = 100./dg/bden() ! mg/kg = g/t ! mm --> m3 *(dart()*1000) ! kg/ha --> kg *(dart()*100) ! >>>>> !~~~~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ use crop, only : crop_yield_output, cva, ndgro, ndpri, nucr, ylda use erosion, only : & cklsp, & css, & erosion_enritchment_ratio, & erosion_organic_nitrate_loss, & erosion_phosphorus_loss, & erosion_soil_loss, & parsz, & xcklsp, & xnorg, & xnorgp, & xporg, & xpsedp use evapotranspiration, only : & canev, & eo, & eopot, & es, & et, & eta_output_id, & etp_output_id, & ra, & snoev, & soil_evaporation_output_id, & tmn use groundwater, only : & additionalGwUptake, & baseflow_output_id, & bff, & groundwater_process, & groundwater_recharge_output_id, & gwchrg, & gwq, & gwseep, & revap, & xet use hydrotope, only : & frar, & hsumfc, & hsumul, & hydrotope_process, & mstruc, & precipitation_output_id, & qi, & soil_water_content_output_id, & swe, & swind, & sym, & tmax_output_id, & tmean_output_id, & tmin_output_id, & tmpNgrw, & tmpNsub, & tmpNsur, & tmpPsur, & water_yield_output_id use landuse, only : & landuse_is_cropland, & landuse_is_forest, & landuse_is_natural_vegetation use management, only : & TSubbasin, & bWAM_Module, & bwamhydrograph, & management_is_transfer_subbasin, & management_subbasin_pointer, & management_transfer_out use nutrient, only : & ano3, & anora, & anors, & degNgrw, & degNsub, & degNsur, & degPsur, & percn, & plab, & pma, & pms, & porg, & retNgrw, & retNsub, & retNsur, & retPsur, & ssf, & ssfn, & yno3, & yon, & yph, & ysp use output, only : & area_tot_glacier, & area_tot_snow, & depth_ave_glacier, & depth_ave_snow, & output_store_hydrotope_value, & output_store_subbasin_value, & output_store_subbasin_values use reservoir, only : xwysb, xxswind use river, only : chk, chl, chw, pet_day, varoute use snow, only : & bSnowModule, & gla, & glac_acc_mm, & glac_acc_mm0, & precip_elev_cor, & precipe, & psnow, & sml, & smle, & sno, & snoa, & snowVal, & snow_acc_mm, & snow_acc_mm0, & snow_degree_day_melting, & snow_depth_weq_output_id, & snow_initialise_subbasin, & snow_process, & snowfall_weq_output_id, & tmax, & tmelt, & tmin, & tmit, & tmx, & tsnfall, & vsn, & xprecip, & xprecip_elev_cor, & xsml, & xsnow, & xtmax, & xtmin, & xtmit, & xvsn use soil, only : & ap, & bd, & cn, & cv, & dur, & fc, & hrtc, & hwss, & iv, & ms, & nn, & ns, & nsa, & op, & percolation_output_id, & pmn, & pr, & psz, & q1, & qd, & rp, & sc, & sccor, & sep, & snum, & soil_acc_mm, & soil_curve_number_alpha, & soil_curve_number_peak_runoff, & soil_curve_transmission_losses, & soil_water_index_output_id, & ste, & stin, & subsurface_runoff_output_id, & surface_runoff_output_id, & swin, & ul, & vo, & wmn, & wn use soil, only : & wno3, & wpo, & xqd, & z use vegetation, only : & alai, & almn, & dm, & ep, & g, & gwqLastday, & rd, & transpiration_output_id, & tsav, & wsav integer, intent(in) :: ihout, inum1 logical, intent(in) :: bSubcatch real(dp), intent(in) :: da real(dp), intent(in) :: da9 integer, intent(in) :: daycounter integer, intent(in) :: ida integer, intent(inout) :: ieap integer, intent(in) :: iy integer, intent(in) :: iyr integer, intent(in) :: mo integer, intent(in) :: nbyr integer, intent(in) :: nd integer :: ii, j, jea, jj, k, l, n, wet real(dp) :: wabad real(dp) :: pcp_tmp = 0. ! temporary precipitation values for output real(dp) :: snow_tmp = 0. ! temporary snow values for output real(dp) :: pcp_tot = 0. ! total preciptiation (rain + snow) for GIS output ! pointer on subbasin #### WATER MANAGEMENT MODULE #### TYPE (TSubbasin), POINTER :: pS !#### INITIALIZATION: call initsub call subbasin_initialise_subbasin(canev, sml, xcklsp, xet, xnorg, xnorgp, xporg, xpsedp, xqd, yon, yph) j = inum1 aff = da * 1000. * flu(j) !**** COMPUTE TOTALS for precip: precip = subp(j) !########################### !#### SNOW MODULE #### !########################### if (bSnowModule) then call snow_initialise_subbasin else susb(1, j) = susb(1, j) + precip sub(1) = sub(1) + precip * flu(j) sbp(j) = sbp(j) + precip !#### COMPUTE snow, recalc precip ! ATTN: Snow accumulation, melt, and snow evaporation ! are calculated at hydrotope level! snowVal = 0. if (tx(j) .le. tsnfall(j)) then snowVal = precip sub(2) = sub(2) + precip * flu(j) call output_store_subbasin_value(snowfall_weq_output_id, j, precip) susb(2, j) = susb(2, j) + precip precip = 0. end if end if !( bSnowModule ) !########################### !*********************************************************** START OF CYCLE 100 ! MAIN CYCLE IN SUBBASIN: 100 - HYDROTOPE caclculation: ! 1) initial conditions, ! 2) call HYDROTOP, ! 3) compute weighted averages for subbasin. do jea = 1, neap(j) n = mstruc(j, jea, 1) ! land use type k = mstruc(j, jea, 2) ! soil number wet = mstruc(j, jea, 3) ! wetland 0 / 1 !########################### !#### SNOW MODULE #### !########################### !#### CALL SNOM: COMPUTE snow & snowmelt, recalc precip ! ATTN: in snom() precip and sno() are recalculated if (bSnowModule) then call snow_process(j, jea, ida, mstruc, precip, tmn, tx) else if (snowVal > 0.) snoa(j, jea) = snoa(j, jea) + snowVal if (tmx(j) > tmelt(j) .AND. snoa(j, jea) > 1E-4 ) call snow_degree_day_melting(j, jea) end if !########################### !**** ESTABLISH number of soil layers: different for arable and non-arable soils nn = ns(k) if (landuse_is_cropland(n) ) nn = nsa(k) !**** CALC hsumul(j, jea), hsumfc(j, jea), hwss(2, j, jea) and hrtc(j, jea, l) call subbasin_assign_hsum(fc, hsumfc, hsumul, hwss, nn, snum, ul) !#### CALL HRFLOWTT call subbasin_flow_travel_time(j, jea, k, bSubcatch, bff, hrtc, nn, sc, sccor, z) !**** INITIAL CONDITIONS for water & nutrients if (iy == 1 .and. ida == 1) call subbasin_init_water_nutrients(n, nn, alai, almn, ano3, anora, anors, ap, cv, cva, dm, g, nucr, op, plab, pma, pmn, pms, porg, rd, sno, snoa, ste, stin, swe, swin, wmn, wn, wno3, wpo) !**** INITIALIZATION of forest - in the beginning of every year ! ATTN! Forest is init. every year, as no forest management is included if (ida .eq. 1) then if (landuse_is_forest(n) ) then dm(j, jea) = 200000. rd(j, jea) = 2000. end if end if !write(*,*) 'Inside subbasin.f95 ::: before olai(1, 1) = ', olai(1,1) !#### CALL HYDROTOPE - calculate processes in HYDROTOPE call hydrotope_process(j, jea, k, n, wet, bSubcatch, dart, daycounter, flu, frar, ida, iy, iyr, mo, nbyr, precip, sbar, tx) !write(*,*) 'Inside subbasin.f95 ::: olai(1, 1) = ', olai(1,1) ! Subroutine wstress is called after ep was calculated in evap.f90 ! This usually decreases ep but et is not updated thereafter. ! Important for hydrotope output! et = ep + es !**** COMPUTE WEIGHTED AVERAGES ! ATTN: xnorgp, xporgp - in kg/ha ! ATTN: xnorg, xporg - coef. (10./bd(k)) to tranf. kg/ha to g/t (I layer!) ! ATTN: xet - not from et, et is not recalc after cropmd sumcn = sumcn + cn * frar(j, jea) xqd = xqd + qd * frar(j, jea) xqi = xqi + qi * frar(j, jea) xssf = xssf + ssf * frar(j, jea) xsep = xsep + sep * frar(j, jea) xswind = xswind + swind * frar(j, jea) xeo = xeo + eopot * frar(j, jea) xet = xet + (ep + es + canev) * frar(j, jea) xsnoev = xsnoev + snoev * frar(j, jea) sno(j) = sno(j) + snoa(j, jea) * frar(j, jea) ! summarise total area covered by snow if (snoa(j, jea) > 1E-4 ) then area_tot_snow = area_tot_snow + frar(j, jea) * sbar(j) depth_ave_snow = depth_ave_snow + snoa(j, jea) * frar(j, jea) * sbar(j) / (da * 10 ** 6) call output_store_hydrotope_value(snow_depth_weq_output_id, j, jea, snoa(j, jea)) if (ida == nd) snow_acc_mm = snow_acc_mm + snoa(j, jea) * frar(j, jea) * sbar(j) / (da * 10 ** 6) - snow_acc_mm0 end if if (ida == nd) soil_acc_mm = soil_acc_mm + swe(j, jea) * frar(j, jea) * sbar(j) / (da * 10 ** 6) !########################### !#### SNOW MODULE #### !########################### if (bSnowModule) then if (gla(j, jea) > 1E-4 ) then area_tot_glacier = area_tot_glacier + frar(j, jea) * sbar(j) depth_ave_glacier = depth_ave_glacier + gla(j, jea) * frar(j, jea) * sbar(j) / (da * 10 ** 6) if (ida == nd) glac_acc_mm = glac_acc_mm + gla(j, jea) * frar(j, jea) * sbar(j) / (da * 10 ** 6) - glac_acc_mm0 end if xtmit = xtmit + tmit * frar(j, jea) xtmax = xtmax + tmax * frar(j, jea) xtmin = xtmin + tmin * frar(j, jea) xprecip = xprecip + precipe * frar(j, jea) ! weighted subbasin rainfall or rainfall + snow melt (not total precipitation), mm xprecip_elev_cor = xprecip_elev_cor + precip_elev_cor * frar(j, jea) ! weighted subbasin elevation - based corrected precipitation, mm xsml = xsml + smle * frar(j, jea) xsnow = xsnow + snowVal * frar(j, jea) ! weighted subbasin snow from snowfall, mm xvsn = xvsn + vsn * frar(j, jea) ! weighted subbasin from snow pack, mm end if !########################### xyno3 = xyno3 + yno3 * frar(j, jea) xssfn = xssfn + ssfn * frar(j, jea) xpercn = xpercn + percn * frar(j, jea) xysp = xysp + ysp * frar(j, jea) xcklsp = xcklsp + cklsp(j, jea) * frar(j, jea) xnorg = xnorg + anors(j, jea, 1) * (10. / bd(k)) * frar(j, jea) xnorgp = xnorgp + anors(j, jea, 1) * frar(j, jea) xporg = xporg + porg(j, jea, 1) * (10. / bd(k)) * frar(j, jea) xporgp = xporgp + porg(j, jea, 1) * frar(j, jea) xpsed = xpsed + (porg(j, jea, 1) + pms(j, jea, 1) + pma(j, jea, 1)) * (10. / bd(k)) * frar(j, jea) xpsedp = xpsedp + (porg(j, jea, 1) + pms(j, jea, 1) + pma(j, jea, 1)) * frar(j, jea) !########################### !#### SNOW MODULE #### !########################### if (bSnowModule) then pcp_tmp = precip_elev_cor pcp_tot = precipe + snowVal else pcp_tmp = precip pcp_tot = subp(j) end if !########################### call output_store_hydrotope_value(subsurface_runoff_output_id, j, jea, ssf) call output_store_hydrotope_value(percolation_output_id, j, jea, sep) call output_store_hydrotope_value(transpiration_output_id, j, jea, ep) call output_store_hydrotope_value(soil_evaporation_output_id, j, jea, es) call output_store_hydrotope_value(soil_water_index_output_id, j, jea, swind) call output_store_hydrotope_value(soil_water_content_output_id, j, jea, swe(j, jea)) !**** CALC variables for GIS output if (bSnowModule) then call output_store_hydrotope_value(precipitation_output_id, j, jea, precipe + snowVal) else call output_store_hydrotope_value(precipitation_output_id, j, jea, subp(j)) end if call output_store_hydrotope_value(surface_runoff_output_id, j, jea, qd) call output_store_hydrotope_value(groundwater_recharge_output_id, j, jea, sep) call output_store_hydrotope_value(eta_output_id, j, jea, ep + es + canev) call output_store_hydrotope_value(etp_output_id, j, jea, eo) !#### CALL CROP_GIS if (iy > 1 .AND. ida == ndpri) call crop_yield_output(j, jea, k, ieap, ms, ndgro, tsav, wsav, ylda, landuse_is_cropland(n)) !**** CALC parsz() for subbasins do jj = 1, 5 parsz(jj, j) = parsz(jj, j) + psz(jj, k) * frar(j, jea) end do end do !*********************************************************** END OF CYCLE 100 pet_day(j) = xeo ! for transmission losses !########################### !#### SNOW MODULE #### !########################### if (bSnowModule) then !#### recalculate precipitation and air temperature for the whole basin precip = xprecip tx(j) = xtmit tmx(j) = xtmax tmn(j) = xtmin sml = xsml psnow = xsnow xprecip = 0. xtmit = 0. xtmax = 0. xtmin = 0. xsml = 0. xsnow = 0. end if !########################### !#### CALL alpha to calc r1 - alpha to be used in peakq if (precip > 0.) then call soil_curve_number_alpha(j, mo, precip, sml) endif !**** COMPUTE PEAK RUNOFF RATE, TRANSMISSION LOSSES & EROSION for subbasin !#### CALL PEAKQ, TRAN, YSED ! xssf (mm), css (ppm) aff=1000*da*flu(j), yd (t) if (xqd .ne. 0.) then call soil_curve_number_peak_runoff(j) if (xqd .gt. 0. .and. pr .gt. 0.) then iv = 1 vo = xqd * da * flu(j) * 1000. dur = vo / (pr * 3600.) if (dur .gt. 24.) dur = 24. ! SL BEGIN: transmission losses are calculated in the routing process using subroutine tran_river ! SL only recalculate peak runoff rate (pr) in function tran(j) q1 = xqd ! SL call soil_curve_transmission_losses(j, chk, chl, chw) !qtl = q1 - xqd xqd = q1 ! if ( qtl.lt.0. ) then ! xqd = q1 ! qtl = 0. ! end if ! revapst(j) = revapst(j) + qtl ! SL END end if call erosion_soil_loss(aff, pr, xqd, yd) end if yd = yd + xssf * aff * css(j) !**** COMPUTE ENRICHMENT RATIO, ORG N and P LOSS WITH EROSION !#### CALL ENRSB, ORGNSED, PSED if (qd > 0. .AND. pr > 0. .AND. precip > 0.) then call erosion_enritchment_ratio(j, da, da9, flu, pr, precip, rp, xqd, yd) call erosion_organic_nitrate_loss(j, da9, yd, yon) call erosion_phosphorus_loss(j, da9, yd, yph) end if !**** COMPUTE GROUND WATER CONTRIBUTION TO FLOW & WATER YIELD !#### CALL GWMOD ! NOTE that revap from shallow aquifer is add to xet ! call gwmod(j) call groundwater_process(j, xsep, xet, xeo) wysb = xqi + xssf + gwq(j) ! - qtl ! water yield, transmission losses are calculated in the routing process runoff_mm(j) = wysb call output_store_subbasin_value(river_runoff_output_id, j, wysb) !**** COMPUTE RETENTION of NUTRIENTS in SUBBASIN ! Retention of N (xyno3, xssfn, xpercn) and soluble P (xysp) ! Method: from Fred Hattermann 2003 ! Updated 11.03.2008 from Shaochun's file xyno3 = (1. - exp(- 1. / retNsur - degNsur)) * xyno3 / (1. + retNsur * degNsur) + tmpNsur(j) * exp(- 1. / retNsur - degNsur) xssfn = (1. - exp(- 1. / retNsub - degNsub)) * xssfn / (1. + retNsub * degNsub) + tmpNsub(j) * exp(- 1. / retNsub - degNsub) xpercn = (1. - exp(- 1. / retNgrw - degNgrw)) * xpercn / (1. + retNgrw * degNgrw) + tmpNgrw(j) * exp(- 1. / retNgrw - degNgrw) xysp = (1. - exp(- 1. / retPsur - degPsur)) * xysp / (1. + retPsur * degPsur) + tmpPsur(j) * exp(- 1. / retPsur - degPsur) tmpNsur(j) = xyno3 tmpNsub(j) = xssfn tmpNgrw(j) = xpercn tmpPsur(j) = xysp !**** COMPUTE TOTALS for water: sub(1...13) !########################### !#### SNOW MODULE #### !########################### if (bSnowModule) then susb(1, j) = susb(1, j) + xprecip_elev_cor ! precip sub(1) = sub(1) + xprecip_elev_cor * flu(j) ! precip * flu(j) sbp(j) = sbp(j) + xprecip_elev_cor ! precip !precipe sub(2) = sub(2) + psnow * flu(j) call output_store_subbasin_value(snowfall_weq_output_id, j, psnow) susb(2, j) = susb(2, j) + psnow end if !########################### sub(3) = sub(3) + sml * flu(j) susb(3, j) = susb(3, j) + sml sub(4) = sub(4) + tmx(j) * flu(j) call output_store_subbasin_value(tmax_output_id, j, tmx(j)) susb(4, j) = susb(4, j) + tmx(j) sub(5) = sub(5) + tmn(j) * flu(j) call output_store_subbasin_value(tmin_output_id, j, tmn(j)) susb(5, j) = susb(5, j) + tmn(j) sub(29) = sub(29) + tx(j) * flu(j) ! / get_days_in_month(mo, iyr) call output_store_subbasin_value(tmean_output_id, j, tx(j)) susb(29, j) = susb(29, j) + tx(j) ! / get_days_in_month(mo, iyr) sub(6) = sub(6) + ra(j) * flu(j) susb(6, j) = susb(6, j) + ra(j) sub(7) = sub(7) + sumcn * flu(j) susb(8, j) = susb(8, j) + xqd sub(8) = sub(8) + xqd * flu(j) smq(j) = smq(j) + xqd susb(9, j) = susb(9, j) + xssf sub(9) = sub(9) + xssf * flu(j) smsq(j) = smsq(j) + xssf ! SL Transmission losses are calculated during the routing process... ! susb(10, j) = susb(10, j) + qtl ! sub(10) = sub(10) + qtl * flu(j) ! revapst(j) = revapst(j) + qtl susb(11, j) = susb(11, j) + xsep sub(11) = sub(11) + xsep * flu(j) susb(12, j) = susb(12, j) + xeo sub(12) = sub(12) + xeo * flu(j) susb(13, j) = susb(13, j) + xet sub(13) = sub(13) + xet * flu(j) !**** COMPUTE TOTALS for ground water: sub(15...19) susb(15, j) = susb(15, j) + gwq(j) sub(15) = sub(15) + gwq(j) * flu(j) call output_store_subbasin_value(baseflow_output_id, j, gwq(j)) susb(16, j) = susb(16, j) + revap sub(16) = sub(16) + revap * flu(j) susb(17, j) = susb(17, j) + gwseep sub(17) = sub(17) + gwseep * flu(j) susb(18, j) = susb(18, j) + gwchrg sub(18) = sub(18) + gwchrg * flu(j) susb(19, j) = susb(19, j) + xswind sub(19) = sub(19) + xswind * flu(j) xxswind = xxswind + xswind * flu(j) !**** Calc daily water balance wabad = sub(1) - sub(8) - sub(9) - sub(11) - sub(13) !**** COMPUTE TOTALS for water yield (wysb) & sediments: sub(20...25) susb(20, j) = susb(20, j) + wysb sub(20) = sub(20) + wysb * flu(j) call output_store_subbasin_value(water_yield_output_id, j, wysb) sub(21) = sub(21) + yd / (100. * da * flu(j)) sym(j) = sym(j) + yd susb(21, j) = susb(21, j) + yd / (100. * da * flu(j)) sub(22) = sub(22) + yon * flu(j) susb(22, j) = susb(22, j) + yon sub(23) = sub(23) + yph * flu(j) susb(23, j) = susb(23, j) + yph sub(24) = sub(24) + ysp * flu(j) susb(24, j) = susb(24, j) + ysp !**** COMPUTE SUBBASIN OUTPUTS FOR ROUTING: sda(), varoute() ! ATTN: sda(6, j) = sum of 3 fluxes after retention, new version ! ATTN: coef (dart()*1000) to transform from mm to m3 (42, 48) ! ATTN: coef (dart()*100) to transform kg/ha to kg (44-47) ! ATTN: wysb(mm)*dart(km2)*1000 = wysb*dart*1000 (m3) ! ATTN: sda(2, j), sda(8, j) in m3 ! ATTN: xyno3 in kg/ha, sda(6) & varoute(6, .) in kg sda(1, j) = precip sda(2, j) = xqi * dart(ihout) * 1000. sda(3, j) = yd sda(4, j) = yon * dart(ihout) * 100 sda(5, j) = yph * dart(ihout) * 100 sda(6, j) = (xyno3 + xssfn + xpercn) * dart(ihout) * 100 sda(7, j) = xysp * dart(ihout) * 100 sda(8, j) = (wysb - xqi) * dart(ihout) * 1000. xwysb = xwysb + wysb * dart(ihout) * 1000. !cc riparian zone take water form GW flow as demanded gwqLastday(j) = sda(8, j) * 0.9 sda(8, j) = sda(8, j) - additionalGwUptake(j) !################################# !#### WATER MANAGEMENT MODULE #### !################################# if (bWAM_Module) then ! check if water transfers take place in current subbasin if (management_is_transfer_subbasin(j) ) then pS => management_subbasin_pointer(j) !------------------------------------------------------------- ! Add supply from water user(s) to subbasins' hydrograph storage location (=j) ! If the subbasin is a reservoir, this step will be overwritten by the reservoir module !------------------------------------------------------------- ! if external inputs from water user(s) are greater than 0 ! add supply of previous day to subbasins' hydrograph storage location (surface runoff component) if (daycounter > 1) then if (pS % inflow(daycounter - 1) > 0. ) & sda(2, j) = sda(2, j) + pS % inflow(daycounter - 1) * 86400. ! convert input from m3 / s to m3 / d end if !------------------------------------------------------------- !------------------------------------------------------------- ! Withdraw water from subbasin outlet ! but only if subbasin is a headwater .and. not a reservoir. ! otherwise withdawals are taken during the routing process in subroutine 'add' ! Moreover, the water users %supplied and %tr_losses are calculated. !------------------------------------------------------------- ! Check if current subbasin is a headwater...in this case, ! bWAMHydrograph is .true. at the array position of the subbasin number ! ATTENTION: The values of sda(2, j) and sda(8, j) may be changed! if (bWAMHydrograph(j) ) then call management_transfer_out(j, sda(2, j), sda(8, j), daycounter, ida, iyr) end if !------------------------------------------------------------- pS % Q_act(daycounter) = real(sda(2, j) / 86400. + sda(8, j) / 86400., 4) end if end if !################################# do ii = 1, 8 varoute(ii, ihout) = sda(ii, j) end do !########################### !#### SNOW MODULE #### !########################### if (bSnowModule) then pcp_tmp = xprecip_elev_cor snow_tmp = psnow else pcp_tmp = precip snow_tmp = snowVal end if !########################### return !****************************************************************************** contains ! Subroutines inside subroutine subbasis subroutine subbasin_assign_hsum(fc, hsumfc, hsumul, hwss, nn, snum, ul) real(dp), dimension(:, :), intent(in) :: fc real(dp), dimension(:, :), intent(inout) :: hsumfc real(dp), dimension(:, :), intent(inout) :: hsumul real(dp), dimension(:, :, :), intent(inout) :: hwss integer, intent(in) :: nn integer, dimension(:), intent(in) :: snum real(dp), dimension(:, :), intent(in) :: ul ! Correction from Fred Hattermann ! Reason: nn is defined for HRU, & sumul(k) and sumfc(k) are not appropriate hsumul(j, jea) = 0. hsumfc(j, jea) = 0. do l = 1, nn hsumul(j, jea) = hsumul(j, jea) + ul(l, k) hsumfc(j, jea) = hsumfc(j, jea) + fc(l, k) end do if (hsumfc(j, jea) .lt. 0.001) then call log_error("subbasin_assign_hsum", "hsumfc less than 0.001, is this soil in soil.cio? (subbasin, hru, landuse, soil)", & ints=(/j, jea, n, snum(k)/)) end if hwss(1, j, jea) = - 1.2677 +log(hsumfc(j, jea)) hwss(2, j, jea) = 1.6234 / hsumfc(j, jea) end subroutine subbasin_assign_hsum !****************************************************************************** !****************************************************************************** subroutine subbasin_init_water_nutrients(n, nn, alai, almn, ano3, anora, anors, ap, cv, cva, dm, g, nucr, op, plab, pma, pmn, pms, porg, rd, sno, snoa, ste, stin, swe, swin, wmn, wn, wno3, wpo) use landuse, only : LULC, landuse_index, landuse_is_cropland, landuse_is_natural_vegetation real(dp), dimension(:, :), intent(inout) :: alai real(dp), dimension(:), intent(in) :: almn real(dp), dimension(:, :, :), intent(inout) :: ano3 real(dp), dimension(:, :, :), intent(inout) :: anora real(dp), dimension(:, :, :), intent(inout) :: anors real(dp), dimension(:, :), intent(in) :: ap real(dp), dimension(:), intent(in) :: cv real(dp), dimension(:, :), intent(inout) :: cva real(dp), dimension(:, :), intent(inout) :: dm real(dp), dimension(:, :), intent(inout) :: g integer, dimension(:, :), intent(inout) :: nucr real(dp), dimension(:, :), intent(in) :: op real(dp), dimension(:, :, :), intent(inout) :: plab real(dp), dimension(:, :, :), intent(inout) :: pma real(dp), dimension(:, :), intent(in) :: pmn real(dp), dimension(:, :, :), intent(inout) :: pms real(dp), dimension(:, :, :), intent(inout) :: porg real(dp), dimension(:, :), intent(inout) :: rd real(dp), dimension(:), intent(in) :: sno real(dp), dimension(:, :), intent(inout) :: snoa real(dp), dimension(:, :, :), intent(inout) :: ste real(dp), dimension(:, :), intent(in) :: stin real(dp), dimension(:, :), intent(inout) :: swe real(dp), dimension(:), intent(in) :: swin real(dp), dimension(:, :), intent(in) :: wmn real(dp), dimension(:, :), intent(in) :: wn real(dp), dimension(:, :), intent(in) :: wno3 real(dp), dimension(:, :), intent(in) :: wpo ! assign initial conditions for water & nutrients at the first ! day of simulation ! land use type integer, intent(in) :: n ! number of soil arable layers integer, intent(in) :: nn integer :: l = 0 do l = 1, nn ste(j, jea, l) = stin(l, k) end do swe(j, jea) = swin(k) !########################### !#### SNOW MODULE #### !########################### if (.NOT. bSnowModule) then snoa(j, jea) = sno(j) end if !########################### !**** Initial conditions for nutrients do l = 1, nn ano3(j, jea, l) = wno3(l, k) anora(j, jea, l) = wmn(l, k) anors(j, jea, l) = wn(l, k) plab(j, jea, l) = ap(l, k) porg(j, jea, l) = wpo(l, k) pma(j, jea, l) = pmn(l, k) pms(j, jea, l) = op(l, k) end do !**** Correction of initial values ano3(), plab() for noncropland if (landuse_is_cropland(n) ) then do l = 1, nn ano3(j, jea, l) = ano3(j, jea, l) * .1 plab(j, jea, l) = plab(j, jea, l) * .1 end do end if !**** INITIALIZATION of crop parameters g(j, jea) = 0. !alai(j, jea)=alai0(get_lu_index(n)) ! set LAI to minimum LAI as indicated in crop.dat if (LULC % veg_code(landuse_index(n)) > 0 ) then alai(j, jea) = almn(LULC % veg_code(landuse_index(n))) else alai(j, jea) = 0. end if !**** Correction of initial values ano3(), plab() for noncropland if (landuse_is_cropland(n) ) then cva(j, jea) = cv(j) else cva(j, jea) = 0. nucr(j, jea) = 0 ! crop number (crop.dat) dm(j, jea) = 0. ! total biomass rd(j, jea) = 0. ! root depth end if if (landuse_is_natural_vegetation(n) ) then dm(j, jea) = 3. rd(j, jea) = 300. end if end subroutine subbasin_init_water_nutrients !****************************************************************************** !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& end subroutine runsubbasin !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& subroutine subbasin_flow_travel_time(j, je, k, bSubcatch, bff, hrtc, nn, sc, sccor, z) !**** PURPOSE: THIS SUBROUTINE CALCULATES RETURN FLOW TRAVEL TIMES !**** CALLED IN: SUBBASIN !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! PARAMETERS & VARIABLES ! ! >>>>> COMMON PARAMETERS & VARIABLES ! bff = baseflow factor for basin, from readbas ! hrtc(j, je, l) = return flow travel time, calculated here, h ! nn = number of soil layers, calc in readsol ! sc(l, k) = sat. cond., input in readsol, mm/h ! z(l, k) = depth, input in readsol, mm ! >>>>> ! >>>>> STATIC PARAMETERS ! dg = local par ! i1 = local par ! ii = local par ! la = local par ! sumk = local par ! xx = local par ! yy = local par ! zk = local par ! >>>>> !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !**** Include common parameters logical, intent(in) :: bSubcatch real(dp), dimension(:), intent(inout) :: bff real(dp), dimension(:, :, :), intent(inout) :: hrtc integer, intent(in) :: nn real(dp), dimension(:, :), intent(in) :: sc real(dp), dimension(:), intent(in) :: sccor real(dp), dimension(:, :), intent(in) :: z integer j, je, k, i1, ii, la real(dp) dg, sumk, xx, yy, zk sumk = 0. yy = 0. zk = 0. do la = 1, nn ii = nn - la + 1 i1 = ii - 1 if (i1 .ne. 0) then xx = z(i1, k) else xx = 0. end if dg = z(ii, k) - xx yy = yy + dg ! SL begin ! sumk = sumk + dg / sc(ii, k) ! zk = yy / sumk ! zk = sc(ii, k) / zk if (bSubcatch) then sumk = sumk + dg / (sc(ii, k) * sccor(j)) zk = yy / sumk zk = (sc(ii, k) * sccor(j)) / zk else sumk = sumk + dg / sc(ii, k) zk = yy / sumk zk = sc(ii, k) / zk end if ! SL end xx = - 0.5447 + 0.01757 * zk if (xx .lt. - 20.) xx = - 20. if (bff(j) .le. 0.) bff(j) = .01 hrtc(j, je, ii) = 10. * (1. - zk / (zk + exp(xx))) / bff(j) end do return end subroutine subbasin_flow_travel_time subroutine subbasin_initialise_weather_gen(amp, avt, ffc, lat, nc, wft, wi) !**** PURPOSE: THIS SUBROUTINE READS MONTHLY STATISTICAL WEATHER PARAMETERS ! for the basin from wgen.dat !**** CALLED IN: readsub !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! PARAMETERS & VARIABLES ! ! >>>>> COMMON PARAMETERS & VARIABLES ! ATTN: Input parameters read from wgen.dat are BASIN parameters; ! Dimension (j) is optional, if subbasin statistics will be availbale ! amp(j) = annual amplitude in daily aver. temperature, degree C ! avt(j) = average annual air temp, degree C, used in solt() ! daylmn(j) = min day length, h ! ffc(j) = field capacity, m/m (not used) !block nc(m) = number of days passed in the beginning of month ! obmn(m, j) = average monthly min temp, degree C ! obmx(m, j) = average monthly max temp, degree C ! prw(1, m, j) = monthly probability of wet day after dry day ! prw(2, m, j) = monthly probability of wet day after wet day ! r(8) = vector for output in readwet ! rsm(m) = monthly max .5h rain for period of record, mm, ! smoothed wim() ! rsmm(m) = monthly number of rainy days ! rsmy(m) = monthly rainfall, mm ! rst(m, j) = monthly mean event of daily rainfall, mm ! tpnyr(j) = number of years of record max .5h rainfall ! tp5(j) = 10 year frequency of .5h rainfall (mm) ! tp6(j) = 10 year frequency of .6h rainfall (mm) ! wft(m, j) = monthly prob. of rainy day ! wi(m, j) = f(wim), used in alpha() for estim of precip. alpha factor ! wim(m) = monthly max .5h rain for period of record, mm ! ylc(j) = cos(lat()/clt), lat() - lat, clt=57.296, for rmx in evap ! yls(j) = sin(lat()/clt), lat() - lat, clt=57.296, for rmx in evap ! (convert degrees to radians (2pi/360=1/57.296) ) ! >>>>> STATIC PARAMETERS ! i = subbasin number (IN TITLE) ! ch = interm. parameter ! f = interm. parameter ! h = interm. parameter ! ii = cycle parameter ! j = cycle parameter ! mon = cycle parameter ! r25 = interm. parameter ! tas = interm. parameter to calc amp() ! tav = interm. parameter, monthly mean temp ! tbb = interm. parameter to calc amp() ! titldum = text ! xm = interm. parameter ! xx = interm. parameter ! xy2 = interm. parameter ! ytn = interm. parameter use input, only : get_config_fid use output, only : output_open_file real(dp), dimension(:), intent(inout) :: amp real(dp), dimension(:), intent(inout) :: avt real(dp), dimension(:), intent(inout) :: ffc real(dp), dimension(:), intent(in) :: lat integer, dimension(13), intent(in) :: nc real(dp), dimension(:, :), intent(inout) :: wft real(dp), dimension(:, :), intent(inout) :: wi integer i, mon real(dp) f, r25, tas, tav, tbb, xm, xx, xy2 real(dp), dimension(12) :: r, rsm, rsmm, rsmy !**** CALCULATION of WEATHER GENERATOR PARAMETERS, step 2 ! wft() - used in solt, wi() - used in alpha & peakq rsm = 0. rsmm = 0. rsmy = 0. rsm(1) = (wim(12) + wim(1) + wim(2)) / 3. do mon = 2, 11 rsm(mon) = (wim(mon - 1) + wim(mon) + wim(mon + 1)) / 3. end do rsm(12) = (wim(11) + wim(12) + wim(1)) / 3. tbb = 0. tas = 100. do i = 1, mb r = 0. do mon = 1, 12 xm = nc(mon + 1) - nc(mon) rsmm(mon) = xm * prw(1, mon) / & (1. - prw(2, mon) + prw(1, mon)) if (rsmm(mon) .le. 0.) rsmm(mon) = .001 r25 = rst(mon) rsmy(mon) = rsmm(mon) * r25 wft(mon, i) = rsmm(mon) / xm xy2 = .5 / tpnyr(i) f = xy2 / rsmm(mon) wi(mon, i) = - rsm(mon) /log(f) wi(mon, i) = 1. - exp(- wi(mon, i) / r25) if (wi(mon, i) .lt. .1) wi(mon, i) = .1 if (wi(mon, i) .gt. .95) wi(mon, i) = .95 r(1) = r(1) + obmx(mon) r(2) = r(2) + obmn(mon) r(8) = r(8) + rsmy(mon) tav = (obmx(mon) + obmn(mon)) / 2. if (tav .gt. tbb) tbb = tav if (tav .lt. tas) tas = tav end do !**** CALCULATION of WEATHER GENERATOR PARAMETERS, step 3 ! avt(), amp() - used in solt xx = lat(i) / clt avt(i) = (r(1) + r(2)) / 2. / 12. amp(i) = (tbb - tas) / 2. xx = r(8) ffc(i) = xx / (xx + exp(9.043 - .002135 * xx)) end do end subroutine subbasin_initialise_weather_gen end module subbasin