! FILE writgen.f ! ! SUBROUTINES IN THIS FILE CALLED FROM ! subroutine wr_daily main ! subroutine wr_month(mo1) main ! subroutine wr_annual main module output use utilities, only : & dp, & log_debug, & log_info, & log_error, & identifier_max_length, & output_max_variables, & path_max_length, & output_label_max_length implicit none integer, parameter :: nih = 7 integer :: nvrch = 18 integer :: nvsub = 30 integer :: nns = 30 ! integer, save :: ieapu ! area covered by snow [m2] real(dp), save :: area_tot_snow ! average snow depth [mm] real(dp), save :: depth_ave_snow ! area covered by glaciers [m2] real(dp), save :: area_tot_glacier ! real(dp), save :: depth_ave_glacier ! average glacier depth [mm], averaged over the entire catchment ! annual SUM of precipitation in subbasin real(dp), save, dimension(:), allocatable :: sbpy ! annual SUBBASIN outputs (dif. components) real(dp), save, dimension(:, :), allocatable :: sysub ! total SUBBASIN outputs = SUM(sysub) real(dp), save, dimension(:, :), allocatable :: stsub integer, save :: giscounter = 1; ! annual REACH outputs (dif. components) real(dp), save, dimension(:, :), allocatable :: syrch ! total REACH outputs = SUM(syrch) real(dp), save, dimension(:, :), allocatable :: strch integer :: nsb = 30 ! Directory to write output to character(len=path_max_length) :: output_dir = "output" ! Output specs namelist path character(len=path_max_length) :: output_specs_path = "input/output.nml" ! Write out interval in timesteps: daily D, monthly M, annual Y character(len=1) :: output_write_interval = "M" ! Float csv output format character(len=identifier_max_length) :: output_float_format = "f12.4" ! Format of space index in csv output character(len=identifier_max_length) :: output_space_index_format = "i6" ! Default file format if none given in output file nml character(len=identifier_max_length) :: output_default_format = "csv" ! Log levels and master logfile character(len=path_max_length) :: master_logfile = "swim.log" character(len=identifier_max_length) :: & log_stderr_level = "warning", & log_stdout_level = "info", & log_file_level = "info" ! Output file attributes (mostly as integer indeces of corresponding module variables) type output_file ! Suffix after <space>_<time>_<name>.<format> character(len=identifier_max_length) :: name ! Requested variables as columns character(len=identifier_max_length) :: variables(output_max_variables) ! Indeces of dimensions integer :: format, space, time ! Number and indeces of variables in output_requested_vars integer :: nvars, variable_ix(output_max_variables) ! ID (and record for binary output) integer :: file_id, current_record = 1 end type ! Output files array type(output_file), dimension(:), allocatable :: output_files ! Variables requested and registered in the code character(len=identifier_max_length), dimension(output_max_variables) :: & output_requested_vars = "", & output_registered_vars = "", & output_hydrotope_vars = "", & output_subbasin_vars = "" ! Dimensions character(len=identifier_max_length), parameter :: & output_space_dim(5) = (/"hydrotope ", "subbasin ", "catchment ", & "hydrotope_label", "subbasin_label "/), & output_time_dim(3) = (/"daily ", "monthly", "annual "/), & output_format_dim(2) = (/"csv", "bin"/) ! Counters integer :: & output_id_hydrotope_counter = 0, & output_id_subbasin_counter = 0, & output_store_day_counter = 1, & output_ndays_per_month = 0, & output_store_month_counter = 1, & output_ndays_per_year = 0, & output_ndays = 0, & output_nmonths = 0, & output_nyears = 0, & output_nvars = 0, & output_ncatchments = 0, & output_nsubbasins = 0, & output_nhydrotopes = 0, & output_nfiles = 0, & output_specs_fid = 0 ! Storage arrays for daily raw variable values (vars, space, time) real, allocatable, dimension(:, :, :) :: & output_storage_hydrotope, & output_storage_subbasin ! Space and time aggregated storage arrays (vars, space, time) real, allocatable, dimension(:, :, :) :: & output_agg_hydrotope_monthly, & output_agg_hydrotope_annual, & output_agg_subbasin_daily, & output_agg_subbasin_monthly, & output_agg_subbasin_annual, & output_agg_catchment_daily, & output_agg_catchment_monthly, & output_agg_catchment_annual integer, allocatable :: & ! Requested output convenience matrix (vars, space, time) output_is_requested(:, :, :), & ! Variable hydrotope/subbasin index in output_requested_vars, 0 if not output_hydrotope_var_ix(:), & output_subbasin_var_ix(:), & ! Indeces of hydrotope raw storage for hydrotope_daily output output_hydrotope_requested_ix(:), & ! Indeces variables in storage arrays (nfiles, nvars) output_storage_index(:, :), & ! Hydrotope/subbasin indeces output_hydrotope_id(:), & output_subbasin_id(:), & output_catchment_id(:), & ! Number of hydrotopes in each subbasin to unsparse 2D hydrotope arrays output_subbasin_nhydrotopes(:), & ! Start index of hydrotope to unsparse 2D hydrotope arrays output_subbasin_hydrotope_stix(:), & ! Hydrotope-subbasin-catchment maps output_hydrotope_subbasin_ix(:), & output_subbasin_catchment_ix(:), & ! Labelled values indeces in storage output_hydrotope_label_ix(:), & output_subbasin_label_ix(:) ! Unit areas for weighted averages real, allocatable, dimension(:) :: & output_hydrotope_subbasin_share, & output_subbasin_catchment_share, & output_catchment_basin_share ! Time indeces for csv output character(len=10), allocatable :: output_day_ix(:) character(len=7), allocatable :: output_month_ix(:) ! Hydrotope/subbasin labels character(len=output_label_max_length), allocatable, dimension(:) :: & output_hydrotope_label, & output_subbasin_label ! Temporal aggregation per variable: average if true, sum if false logical, allocatable :: output_var_state(:) namelist / OUTPUT_PARAMETERS / & nns, & nsb, & nvrch, & nvsub, & output_dir, & output_write_interval, & output_float_format, & output_space_index_format, & output_default_format, & master_logfile, & log_stderr_level, & log_stdout_level, & log_file_level contains subroutine output_allocate(mb, mch) integer, intent(in) :: mb, mch allocate(sbpy(mb)) allocate(strch(18, mch)) allocate(stsub(30, mb)) allocate(syrch(18, mch)) allocate(sysub(30, mb)) sbpy = 0. strch = 0. stsub = 0. syrch = 0. sysub = 0. end subroutine output_allocate subroutine dealloc_output deallocate(sbpy) deallocate(strch) deallocate(stsub) deallocate(syrch) deallocate(sysub) end subroutine dealloc_output integer function output_open_file(file_name) result(file_id) ! Abstract output file opener (so far only for formatted output) use utilities, only : open_file character(len=*), intent(in) :: file_name file_id = open_file(trim(output_dir)//"/"//trim(file_name), "w") end function output_open_file !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! PARAMETERS & VARIABLES ! ! susb(30, j) = monthly SUBBASIN outputs (dif. components) sysub ! sysub(30, j) = annual SUBBASIN outputs (dif. components) ! stsub(30, j) = total SUBBASIN outputs = SUM(sysub) ! sub(30) = daily BASIN outputs (dif. components, weighted sums) ! smm(30) = monthly BASIN outputs (dif. components) ! smy(30) = annual BASIN outputs (dif. components) ! sm(30) = average annual BASIN outputs (dif. components) ! srch(18, ir) = monthly REACH outputs (dif. components) ! syrch(18, ir)= annual REACH outputs (dif. components) ! strch(18, ir)= total REACH outputs = SUM(syrch) ! sbp(j) = monthly SUM of precipitation in subbasin ! smq(j) = monthly SUM of surface runoff for subbasin ! smsq(j) = monthly SUM of subsurface runoff for subbasin ! sym(j) = monthly SUM of sediment yield for subbasins ! syq(j) = SUM(smq) annual SUM surface runoff for subbasins ! sysq(j) = SUM(smsq) annual SUM sub-surface runoff for subbasins ! syy(j) = SUM(sym) annual SUM sed. yield for subbasins ! sq(j) = SUM(syq) total SUM (whole period) surf. runoff in subb. ! ssq(j) = SUM(sysq) total SUM (whole period) sub-surf. runoff in subb. ! sy(j) = SUM(syy) total SUM (whole period) sed. yield in subb. ! DESCRIPTION of BASIN OUTPUTS sub() (Analogue - for susb, smm, smy): ! 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. ! sub(5) = SUM(tmn*flu) average "weighted" min temp. ! sub(6) = SUM(ra*flu) radiation ! 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 q, mm ! sub(16) = SUM(revap*flu) revap from g-w to soil prof., mm ! sub(17) = SUM(gwseep*flu) g-w percolation ! sub(18) = SUM(gwchrg*flu) g-w recharge ! sub(19) = SUM(xswimd*flu) soil water ! sub(20) = SUM(wysb*flu) wysb=qi+ssf+gwq(j)-qtl, water yield, mm ! sub(21) = SUM(yd/(100*da*flu) sed yield ! sub(22) = SUM(yon*flu) org. N loss with sed ! sub(23) = SUM(yph*flu) org P loss with sed ! sub(24) = SUM(ysp*flu) sol P ! sub(25) = SUM(xyno3*flu) no3 in sur. ! sub(26) = SUM(xssfn*flu) no3 in sub sur ! sub(27) = SUM(xpercn*flu) no3 leached ! sub(28) = SUM(xuno3*flu) N uptake by plants ! Analogue - for susb, smm, smy: ! SUB(1)= DAILY PRECIPITATION ! SMM(1)= TOTAL MONTHLY PRECIPITATION ! SMY(1)= TOTAL YEARLY PRECIPITATION ! DESCRIPTION OF monthly REACH OUTPUTS ! srch(1, ir) = SUM(varoute(2, )/86400.) surface flow, inflow ! srch(2, ir) = SUM(xxqd/86400.) surface flow, outflow ! srch(3, ir) = SUM(varoute(3, )) sediment in ! srch(4, ir) = SUM(yd) sediment out ! srch(5, ir) = SUM(sedcon) sediment conc. ! srch(6, ir) = SUM(varoute(4, )) organic N in ! srch(7, ir) = SUM(yon*dart(ihout)*100) organic N out ! srch(8, ir) = SUM(varoute(5, 2)) organic P in ! srch(9, ir) = SUM(yph*dart(ihout)*100) organic P out ! srch(10, ir) = SUM(evp/86400.) evaporation (not active!) ! srch(11, ir) = SUM(tlc/86400. transmission losses (not active!) ! srch(12, ir) = SUM(rl/86400. seepage (not active!) ! srch(13, ir) = SUM(diver/86400.) diversion (not active!) ! srch(14, ir) = SUM(rflow/86400.) return flow (not active!) ! srch(15, ir) = SUM(varoute(6, 2)) nitrate N in ! srch(16, ir) = SUM(xxnit) nitrate N out ! srch(17, ir) = SUM(varoute(7, 2)*100) soluble P in ! srch(18, ir) = SUM(xysp*dart(ihout)*100) soluble P out ! Analogue - for syrch = SUM(srch) - annual sums !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ subroutine output_nashsutcliffe_efficiency(qo, qs, icd, istyr, iy) !**** PURPOSE: THIS SUBROUTINE COMPUTES CRITERIA OF FIT: ! difference in calulated water balance, relative difference**2, ! Nash&Sutcliffe Efficiency & LOG-Nash&Sutcliffe Efficiency !**** CALLED IN: MAIN !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! PARAMETERS & VARIABLES ! ! >>>>> PARAMETERS & VARIABLES in TITLE ! qo(1:inn) = observed water discharge, m3/sec. ! qs(1:inn) = simulated water discharge, m3/sec. ! inn = number of days ! icd = code: 1 - one year, 2 - whole period ! >>>>> COMMON PARAMETERS ! istyr = starting year ! iy = current year ! >>>>> ! >>>>> STATIC PARAMETERS ! akk(20) = accumulated values, internal ! crdif = difference in calulated water balance ! crdifp = difference in calulated water balance, % ! dif = difference between qs and qo ! eff = Nash & Sutcliffe Efficiency for variables qo and qs ! efflog = Nash & Sutcliffe Efficiency for log(qo) and log(qs) ! f00 = local par ! f11 = local par ! ik = local par ! im = local par ! qmid = local par ! qmid2 = local par ! reldif2 = relative difference **2 ! xdif = xqo - xqs ! xqo = log(qo) ! xqs = log(qs) ! >>>>> !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !**** Include common parameters & descriptions integer, intent(in) :: istyr, icd, iy real(dp), dimension(:), intent(in) :: qo, qs integer ik, im real(dp) efflog, f00, f11, qmid, qmid2, reldif2, xdif, xqo, xqs, akk(20), crdif, crdifp, dif, eff eff = 0. efflog = 0. do im = 1, 20 akk(im) = 0. end do do ik = 1, size(qo) dif = qs(ik) - qo (ik) qmid = 0.5 * (abs(qo(ik)) + qs(ik)) qmid2 = qmid * qmid akk(1) = akk(1) + 1. akk(2) = akk(2) + qo(ik) akk(3) = akk(3) + qo(ik) * qo(ik) akk(4) = akk(4) + dif akk(5) = akk(5) + dif * dif if (qmid2 .gt. 0.001) then akk(6) = akk(6) + dif * dif / qmid2 else endif if (qo(ik) .gt. 0. .and. qs(ik) .gt. 0.) then xqo = log(qo(ik)) xqs = log(qs(ik)) xdif = xqo - xqs akk(7) = akk(7) + 1. akk(8) = akk(8) + xqo akk(9) = akk(9) + xqo * xqo akk(10) = akk(10) + xdif * xdif endif end do crdif = akk(4) if (akk(2) .ne. 0) then crdifp = akk(4) * 100 / akk(2) else crdifp = 0. endif reldif2 = akk(6) if (akk(1) .gt. 0.) then f00 = akk(3) - akk(2) * akk(2) / akk(1) if (f00 .gt. 0.0001) eff = (f00 - akk(5)) / f00 endif if (akk(7) .gt. 0.) then f11 = akk(9) - akk(8) * akk(8) / akk(7) if (f11 .gt. 0.0001) efflog = (f11 - akk(10)) / f11 endif if (icd .eq. 1) then call log_debug("output_nashsutcliffe_efficiency", & 'CRITERIA of fit in (year): Abs. difference, Difference %, Rel.dif**2, L-Efficiency, Efficiency:', & i1=istyr+iy-1, reals=(/crdif, crdifp, reldif2, efflog, eff/)) else call log_debug("output_nashsutcliffe_efficiency", & 'CRITERIA of fit TOTAL: Abs. difference, Difference %, Rel.dif**2, L-Efficiency, Efficiency:', & reals=(/crdif, crdifp, reldif2, efflog, eff/)) endif return end subroutine output_nashsutcliffe_efficiency subroutine output_initialise(hydrotope_id, hydrotope_area, hydrotope_subbasin, & subbasin_id, subbasin_nhydrotopes, subbasin_catchment, catchment_id, & hydrotope_input_fid, subbasin_input_fid) ! Process user input before variables are registered use utilities, only: string_index integer, dimension(:), intent(in) :: & ! Hydrotope IDs and subbasin mapping (must be same length) hydrotope_id, hydrotope_subbasin, & ! Subbasin IDs and catchment mapping (must be same length) subbasin_id, subbasin_nhydrotopes, subbasin_catchment, & !Catchment IDs catchment_id ! Hydrotope area for weighted averages real, dimension(:), intent(in) :: hydrotope_area ! File ids to read output labels integer, intent(in) :: hydrotope_input_fid integer, intent(in) :: subbasin_input_fid integer s ! Hydrotope arrays output_nhydrotopes = size(hydrotope_id) allocate(output_hydrotope_id(output_nhydrotopes)) allocate(output_hydrotope_subbasin_ix(output_nhydrotopes)) allocate(output_hydrotope_subbasin_share(output_nhydrotopes)) output_hydrotope_id = hydrotope_id ! Subbasin arrays output_nsubbasins = size(subbasin_id) allocate(output_subbasin_id(output_nsubbasins)) allocate(output_subbasin_catchment_ix(output_nsubbasins)) ! initialised with 2 because if catchment ids not read with subcatch = False ! 2 is the index of the default catchment ID, 1 being the total catchment output_subbasin_catchment_ix = 2 allocate(output_subbasin_catchment_share(output_nsubbasins)) output_subbasin_id = subbasin_id allocate(output_subbasin_nhydrotopes(output_nsubbasins)) output_subbasin_nhydrotopes = subbasin_nhydrotopes allocate(output_subbasin_hydrotope_stix(output_nsubbasins)) output_subbasin_hydrotope_stix = 0 ! Make cumsum out of sparse hydrotope numbers for quick start index do s = 1, output_nsubbasins-1 output_subbasin_hydrotope_stix(s+1) = & output_subbasin_hydrotope_stix(s) + subbasin_nhydrotopes(s) end do ! Catchment arrays output_ncatchments = size(catchment_id) + 1 ! one extra for catchment total allocate(output_catchment_id(output_ncatchments)) allocate(output_catchment_basin_share(output_ncatchments)) output_catchment_id(1) = 0 output_catchment_id(2:output_ncatchments) = catchment_id allocate(output_storage_index(output_nfiles, output_nvars)) output_storage_index = 0 call output_initialise_unit_shares(hydrotope_subbasin, subbasin_catchment, hydrotope_area) call output_initialise_time ! Fill convenience matrix call output_initialise_is_requested call output_initialise_storage(hydrotope_input_fid, subbasin_input_fid) call output_allocate(output_nsubbasins, output_nsubbasins) end subroutine output_initialise subroutine output_initialise_unit_shares( & hydrotope_subbasin, subbasin_catchment, hydrotope_area) ! Get array indeces of subbasin and catchment maps and create areal weights integer, dimension(:), intent(in) :: hydrotope_subbasin, subbasin_catchment real, dimension(:), intent(in) :: hydrotope_area integer h, s, c output_catchment_basin_share = 0 output_subbasin_catchment_share = 0 ! indeces and area sums do s = 1, output_nsubbasins do h = 1, output_nhydrotopes if (hydrotope_subbasin(h) == output_subbasin_id(s)) then output_hydrotope_subbasin_ix(h) = s output_subbasin_catchment_share(s) = & output_subbasin_catchment_share(s) + hydrotope_area(h) end if end do end do do c = 2, output_ncatchments do s = 1, output_nsubbasins if (subbasin_catchment(s) == output_catchment_id(c)) then output_subbasin_catchment_ix(s) = c output_catchment_basin_share(c) = output_catchment_basin_share(c) & + output_subbasin_catchment_share(s) end if end do end do ! Make shares output_hydrotope_subbasin_share = hydrotope_area / & output_subbasin_catchment_share(output_hydrotope_subbasin_ix) output_subbasin_catchment_share = output_subbasin_catchment_share & / output_catchment_basin_share(output_subbasin_catchment_ix) output_catchment_basin_share(1) = sum(hydrotope_area) ! entire catchment output_catchment_basin_share = output_catchment_basin_share & / output_catchment_basin_share(1) end subroutine output_initialise_unit_shares subroutine output_initialise_user_input ! Read and interpret output namelist parameters use utilities, only: & extend_unique_string, & string_index, & open_file, & log_create, & master_log use input, only: get_config_fid integer ios, i, v character(len=identifier_max_length) :: & name, format, variables(output_max_variables), space, time namelist / file / & name, format, variables, space, time ! General options from central config nml read(get_config_fid(), nml=OUTPUT_PARAMETERS) ! Create master log file master_log = log_create(trim(output_dir)//"/"//trim(master_logfile), & log_stderr_level, log_stdout_level, log_file_level) ! count requested output files output_specs_fid = open_file(output_specs_path) output_nfiles = -1 ios = 0 do while (ios /= -1) output_nfiles = output_nfiles + 1 read(output_specs_fid, nml=file, iostat=ios) end do rewind(output_specs_fid) if (output_nfiles == 0) & call log_error("output_initialise_user_input", & "No output file definition found in "//output_specs_path) ! Allocate arrays of nfiles length allocate(output_files(output_nfiles)) do i = 1, output_nfiles name = '' format = output_default_format variables = '' space = '' time = '' read(output_specs_fid, nml=file) ! Check input and store if (trim(name) == '') & call log_error("output_initialise_user_input", "No name given for output file.") output_files(i)%name = name output_files(i)%nvars = string_index('', variables) - 1 if (output_files(i)%nvars == 0) & call log_error("output_initialise_user_input", & "No variables define for output file "//trim(name)) output_files(i)%variables = variables output_files(i)%space = string_index(space, output_space_dim) if (output_files(i)%space == 0) call log_error("output_initialise_user_input", & trim(space)//" is not a valid output space domain in file "//trim(name)) output_files(i)%time = string_index(time, output_time_dim) if (output_files(i)%time == 0) call log_error("output_initialise_user_input", & trim(time)//" is not a valid output time domain in file "//trim(name)) output_files(i)%format = string_index(format, output_format_dim) if (output_files(i)%format == 0) call log_error("output_initialise_user_input", & trim(format)//" is not a valid output format in file "//trim(name)) ! Store unique variables call extend_unique_string(variables, output_requested_vars) end do ! Get variable indeces for each file%variables in output_requested_vars do i = 1, output_nfiles do v = 1, output_files(i)%nvars output_files(i)%variable_ix(v) = string_index(output_files(i)%variables(v), & output_requested_vars) end do end do ! Count all requested variables output_nvars = string_index("", output_requested_vars) - 1 allocate(output_var_state(output_nvars)) output_var_state = .True. end subroutine output_initialise_user_input subroutine output_initialise_is_requested type(output_file) :: file integer i, v, sp ! Allocate without labelled space units as they are just copied allocate(output_is_requested(output_nvars, 3, size(output_time_dim))) output_is_requested = 0 ! Marke 1 where needed in (var, space, time) matrix do i = 1, output_nfiles file = output_files(i) do v = 1, file%nvars sp = file%space ! Make labelled count the same as unlabelled if (sp > 3) sp = sp - 3 output_is_requested(file%variable_ix(v), sp, file%time) = 1 end do end do end subroutine output_initialise_is_requested integer function output_register_hydrotope_var(name, state) ! Register hydrotope variable name and optionally "state" (default .True.). ! Returns output ID if variable requested for output, otherwise 0. ! If state = .False. temporal aggregation is sum not average. character(*), intent(in) :: name logical, intent(in), optional :: state logical s s = .True. if (present(state)) s = state call output_register_var(name, s, output_id_hydrotope_counter, & output_hydrotope_vars, output_register_hydrotope_var) end function output_register_hydrotope_var integer function output_register_subbasin_var(name, state) ! Register subbasin variable name and optionally "state" (default .True.). ! Returns output ID if variable requested for output, otherwise 0. ! If state = .False. temporal aggregation is sum not average. character(*), intent(in) :: name logical, intent(in), optional :: state logical s s = .True. if (present(state)) s = state call output_register_var(name, s, output_id_subbasin_counter, & output_subbasin_vars, output_register_subbasin_var) end function output_register_subbasin_var subroutine output_register_var(name, state, counter, vars, id) ! Abstract subroutine of register_hydrotope/subbasin_var functions. use utilities, only: string_index character(*), intent(in) :: name logical, intent(in) :: state integer, intent(inout) :: counter character(*), intent(out) :: vars(:) integer, intent(out) :: id integer vi ! make sure name hasnt been registered already if (string_index(name, vars) > 0) & call log_error("output_register_var", trim(name)//" has already been registered.") ! is variable requested? vi = string_index(name, output_requested_vars) output_registered_vars(string_index("", output_registered_vars)) = name ! increase counter if requested if (vi > 0) then counter = counter + 1 vars(counter) = name id = counter output_var_state(vi) = state else id = 0 end if end subroutine output_register_var subroutine output_initialise_storage(hydrotope_input_fid, subbasin_input_fid) ! Initialisation storage after all variables have been registered registration. integer, intent(in) :: hydrotope_input_fid, subbasin_input_fid integer :: nd, nm, ny integer :: nspacetime(3, 3) call output_initialise_labelled(hydrotope_input_fid, subbasin_input_fid) call output_check_requested_vars call output_open_files call output_initialise_storage_index nd = output_ndays nm = output_nmonths ny = output_nyears ! Raw daily stores allocate(output_storage_hydrotope(output_id_hydrotope_counter, output_nhydrotopes, nd)) allocate(output_storage_subbasin(output_id_subbasin_counter, output_nsubbasins, nd)) ! Total variables requested for (space, time), maxval because [csv, bin] = 1 nspacetime = sum(output_is_requested, dim=1) allocate(output_agg_hydrotope_monthly(nspacetime(1, 2), output_nhydrotopes, nm)) allocate(output_agg_hydrotope_annual(nspacetime(1, 3), output_nhydrotopes, ny)) allocate(output_agg_subbasin_daily(nspacetime(2, 1), output_nsubbasins, nd)) allocate(output_agg_subbasin_monthly(nspacetime(2, 2), output_nsubbasins, nm)) allocate(output_agg_subbasin_annual(nspacetime(2, 3), output_nsubbasins, ny)) allocate(output_agg_catchment_daily(nspacetime(3, 1), output_ncatchments, nd)) allocate(output_agg_catchment_monthly(nspacetime(3, 2), output_ncatchments, nm)) allocate(output_agg_catchment_annual(nspacetime(3, 3), output_ncatchments, ny)) ! set all to 0 output_storage_hydrotope = 0 output_storage_subbasin = 0 output_agg_hydrotope_monthly = 0 output_agg_hydrotope_annual = 0 output_agg_subbasin_daily = 0 output_agg_subbasin_monthly = 0 output_agg_subbasin_annual = 0 output_agg_catchment_daily = 0 output_agg_catchment_monthly = 0 output_agg_catchment_annual = 0 end subroutine output_initialise_storage subroutine output_initialise_time ! decide time dimension/index output_ndays = 1 output_nmonths = 1 output_nyears = 1 if (output_write_interval == "M") then output_ndays = 31 else if (output_write_interval == "Y") then output_ndays = 366 output_nmonths = 12 else if (output_write_interval /= "D") then call log_error("output_initialise_time", & "output_write_interval must be one of D, M, Y not "//trim(output_write_interval)) end if allocate(output_day_ix(output_ndays)) allocate(output_month_ix(output_nmonths)) end subroutine output_initialise_time subroutine output_initialise_storage_index use utilities, only: string_index, indeces_in_array integer fi, sp, i, cumsum(output_nvars) type(output_file) :: fl do fi = 1, output_nfiles fl = output_files(fi) sp = fl%space if (sp > 3) sp = sp - 3 cumsum = 0 do i = 1, output_nvars cumsum(i) = sum(output_is_requested(:i, sp, fl%time)) end do output_storage_index(fi, :fl%nvars) = cumsum(fl%variable_ix(:fl%nvars)) end do ! hydrotope/subbasin storage indeces allocate(output_hydrotope_var_ix(output_nvars)) output_hydrotope_var_ix = 0 allocate(output_subbasin_var_ix(output_nvars)) output_subbasin_var_ix = 0 allocate(output_hydrotope_requested_ix(sum(output_is_requested(:, 1, 1)))) output_hydrotope_requested_ix = 0 ! Hydrotope/subbasin indeces output_hydrotope_requested_ix = indeces_in_array( & pack(output_requested_vars(:output_nvars), output_is_requested(:, 1, 1) == 1), & output_hydrotope_vars) output_hydrotope_var_ix = indeces_in_array( & output_requested_vars(:output_nvars), output_hydrotope_vars) output_subbasin_var_ix = indeces_in_array(& output_requested_vars(:output_nvars), output_subbasin_vars) end subroutine output_initialise_storage_index subroutine output_initialise_labelled(hydrotope_input_fid, subbasin_input_fid) use utilities, only: string_index !use input, only: read_string_column use input, only: read_string_column, has_column integer, intent(in) :: hydrotope_input_fid, subbasin_input_fid integer i, ucount logical ir ir = any(output_files(:)%space == 4) if (ir .and. has_column("output_label", hydrotope_input_fid)) then allocate(output_hydrotope_label(output_nhydrotopes)) output_hydrotope_label = '' call read_string_column(hydrotope_input_fid, 'output_label', output_hydrotope_label) ucount = 0 do i = 1, output_nhydrotopes if (trim(output_hydrotope_label(i)) /= '') ucount = ucount + 1 end do allocate(output_hydrotope_label_ix(ucount)) ucount = 0 do i = 1, output_nhydrotopes if (trim(output_hydrotope_label(i)) /= '') then ucount = ucount + 1 output_hydrotope_label_ix(ucount) = i end if end do else if (ir) then call log_error("output_initialise_labelled", & "hydrotope_label output requested but output_label column not found.") end if ir = any(output_files(:)%space == 5) if (ir .and. has_column("output_label", subbasin_input_fid)) then allocate(output_subbasin_label(output_nsubbasins)) output_subbasin_label = '' call read_string_column(subbasin_input_fid, 'output_label', output_subbasin_label) ucount = 0 do i = 1, output_nsubbasins if (trim(output_subbasin_label(i)) /= '') ucount = ucount + 1 end do allocate(output_subbasin_label_ix(ucount)) ucount = 0 do i = 1, output_nsubbasins if (trim(output_subbasin_label(i)) /= '') then ucount = ucount + 1 output_subbasin_label_ix(ucount) = i end if end do else if (ir) then call log_error("output_initialise_labelled", & "subbasin_label output requested but output_label column not found.") end if end subroutine output_initialise_labelled subroutine output_open_files ! Open all needed output files and write csv headers use utilities, only : open_file integer fi, fid, rl type(output_file) :: fl character(len=path_max_length) :: fname integer, parameter :: nbytes_per_real = 4 do fi = 1, output_nfiles fl = output_files(fi) fname = trim(output_space_dim(fl%space))//"_"//& trim(output_time_dim(fl%time))//"_"//trim(fl%name)// & "."//trim(output_format_dim(fl%format)) if (fl%format == 1) then fid = output_open_file(fname) else if (fl%format == 2) then rl = output_nhydrotopes if (fl%space == 2) rl = output_nsubbasins if (fl%space == 3) rl = output_ncatchments if (fl%space == 4) rl = size(output_hydrotope_label_ix) if (fl%space == 5) rl = size(output_subbasin_label_ix) rl = rl * nbytes_per_real * fl%nvars ! TODO implement binary open file in open_file open(fid, file=fname, form="unformatted", access="direct", recl=rl) end if output_files(fi)%file_id = fid ! Write header if needed if (fl%format == 1) & call output_write_csv_header(fid, output_space_dim(fl%space), & output_requested_vars(fl%variable_ix(:fl%nvars))) end do end subroutine output_open_files subroutine output_write_csv_header(fid, space, variables) ! Write csv header <time> <space> <variables> integer, intent(in) :: fid character(len=identifier_max_length), intent(in) :: space, variables(:) character(len=identifier_max_length) :: fmt integer i, ffl ! Format lengths as int read(output_float_format(2:index(output_float_format, ".") - 1), *) ffl ! time write(fid, "(a)", advance="no") "time, " ! space write(fid, "(a)", advance="no") trim(space)//"," ! variables do i = 1, size(variables) - 1 ! Make sure variable name fits into the size of real format write(fmt, "('(a',i4,')')") max(ffl, len(trim(variables(i)))+1) + 1 write(fid, trim(fmt), advance="no") trim(variables(i))//"," end do ! Last column write(fmt, "('(a',i4,')')") max(ffl, len(trim(variables(size(variables))))+1) write(fid, fmt) trim(variables(size(variables))) end subroutine output_write_csv_header subroutine output_check_requested_vars ! Make sure all requested variables have been registered. use utilities, only: string_index character(len=identifier_max_length) vs integer i, ih, is do i = 1, output_nvars vs = output_requested_vars(i) ih = string_index(vs, output_hydrotope_vars) is = string_index(vs, output_subbasin_vars) ! Is requested variable registered? if (ih + is == 0) call log_error("output_check_requested_vars", & "Requested output variable invalid: "//vs) ! Is subbasin variable requested at hydrotope level? if (ih == 0 .and. maxval(output_is_requested(i, 1, :)) > 0) then call log_error("output_check_requested_vars", & "Requested variable at hydrotope level is not a hydrotope variable: "//vs) end if end do end subroutine output_check_requested_vars subroutine output_store_hydrotope_values(id, array) ! Store hydrotope values at the end of timestep. integer, intent(in) :: id real, dimension(:, :), intent(in) :: array integer s, ns, nh if (id > 0) then ! Unsparse hydrotope data do s = 1, output_nsubbasins ns = output_subbasin_hydrotope_stix(s) nh = output_subbasin_nhydrotopes(s) output_storage_hydrotope(id, ns:ns+nh-1, output_store_day_counter) = array(s, :nh) end do end if end subroutine output_store_hydrotope_values subroutine output_store_subbasin_values(id, array) ! Store subbasin values at the end of timestep. integer, intent(in) :: id real, dimension(:), intent(in) :: array if (id > 0) output_storage_subbasin(id, :, output_store_day_counter) = array end subroutine output_store_subbasin_values subroutine output_store_hydrotope_value(id, subbasinix, hydrotopeix, value) ! Store hydrotope values at the end of timestep. integer, intent(in) :: id ! Subbasin and hydrotope (within subbasin) indeces integer, intent(in) :: subbasinix, hydrotopeix real, intent(in) :: value integer ns if (id > 0) then ns = output_subbasin_hydrotope_stix(subbasinix) output_storage_hydrotope(id, ns+hydrotopeix, output_store_day_counter) = value end if end subroutine output_store_hydrotope_value subroutine output_store_subbasin_value(id, subbasinix, value) ! Store a single subbasin value at the end of timestep. integer, intent(in) :: id, subbasinix real, intent(in) :: value if (id > 0) output_storage_subbasin(id, subbasinix, output_store_day_counter) = value end subroutine output_store_subbasin_value function output_hydrotope_to_subbasin(hydrotope_values) result(subbasin_values) ! Average hydrotope values to subbasins real, intent(in) :: hydrotope_values(output_nhydrotopes) real :: subbasin_values(output_nsubbasins) integer h, hsix real share subbasin_values = 0 do h = 1, output_nhydrotopes hsix = output_hydrotope_subbasin_ix(h) share = hydrotope_values(h) * output_hydrotope_subbasin_share(h) subbasin_values(hsix) = subbasin_values(hsix) + share end do end function output_hydrotope_to_subbasin function output_hydrotope_to_catchment(hydrotope_values) result(catchment_values) ! Average hydrotope values to catchment real, intent(in) :: hydrotope_values(output_nhydrotopes) real :: catchment_values(output_ncatchments) integer h, hsix, scix real share catchment_values = 0 do h = 1, output_nhydrotopes hsix = output_hydrotope_subbasin_ix(h) scix = output_subbasin_catchment_ix(hsix) ! Share of hydrotope in catchment share = hydrotope_values(h) * output_hydrotope_subbasin_share(h) & * output_subbasin_catchment_share(hsix) catchment_values(scix) = catchment_values(scix) + share ! Entire basin catchment_values(1) = & catchment_values(1) + share * output_catchment_basin_share(scix) end do end function output_hydrotope_to_catchment function output_subbasin_to_catchment(subbasin_values) result(catchment_values) ! Average subbasin values to catchment real, intent(in) :: subbasin_values(output_nsubbasins) real :: catchment_values(output_ncatchments) integer s, scix real share catchment_values = 0 do s = 1, output_nsubbasins scix = output_subbasin_catchment_ix(s) share = subbasin_values(s) * output_subbasin_catchment_share(s) catchment_values(scix) = catchment_values(scix) + share catchment_values(1) = & catchment_values(1) + share * output_catchment_basin_share(scix) end do end function output_subbasin_to_catchment subroutine output_day(year, month, day) ! Aggregate space/time and write to csv/bin, if interval is daily. integer, intent(in) :: year, month, day integer i, dc, mc integer :: sti(3, 3), stb(3, 3) real :: hydval(output_nhydrotopes), subval(output_nsubbasins), & catval(output_ncatchments) dc = output_store_day_counter mc = output_store_month_counter ! Format day label for csv output write(output_day_ix(dc), "(i0.4, '-', i0.2, '-', i0.2)") year, month, day ! aggregate day to month do i = 1, output_nvars ! 0/1 (space, time) for current variable stb = output_is_requested(i, :, :) ! indeces in aggregation arrays of current variable sti = sum(output_is_requested(:i, :, :), dim=1) ! current hydrotope/subbasin variables if (output_hydrotope_var_ix(i) > 0) then hydval = output_storage_hydrotope(output_hydrotope_var_ix(i), :, dc) subval = output_hydrotope_to_subbasin(hydval) catval = output_hydrotope_to_catchment(hydval) end if if (output_subbasin_var_ix(i) > 0) then hydval = 0 subval = output_storage_subbasin(output_subbasin_var_ix(i), :, dc) catval = output_subbasin_to_catchment(subval) end if ! Add daily values to aggregation arrays as needed if (stb(1, 2) == 1) output_agg_hydrotope_monthly(sti(1, 2), :, mc) = & output_agg_hydrotope_monthly(sti(1, 2), :, mc) + hydval if (stb(1, 3) == 1) output_agg_hydrotope_annual(sti(1, 3), :, 1) = & output_agg_hydrotope_annual(sti(1, 3), :, 1) + hydval if (stb(2, 1) == 1) output_agg_subbasin_daily(sti(2, 1), :, dc) = & output_agg_subbasin_daily(sti(2, 1), :, dc) + subval if (stb(2, 2) == 1) output_agg_subbasin_monthly(sti(2, 2), :, mc) = & output_agg_subbasin_monthly(sti(2, 2), :, mc) + subval if (stb(2, 3) == 1) output_agg_subbasin_annual(sti(2, 3), :, 1) = & output_agg_subbasin_annual(sti(2, 3), :, 1) + subval if (stb(3, 1) == 1) output_agg_catchment_daily(sti(3, 1), :, dc) = & output_agg_catchment_daily(sti(3, 1), :, dc) + catval if (stb(3, 2) == 1) output_agg_catchment_monthly(sti(3, 2), :, mc) = & output_agg_catchment_monthly(sti(3, 2), :, mc) + catval if (stb(3, 3) == 1) output_agg_catchment_annual(sti(3, 3), :, 1) = & output_agg_catchment_annual(sti(3, 3), :, 1) + catval end do ! Increase/reset day counter if (output_write_interval == "D") then call output_write_daily output_store_day_counter = 0 end if output_store_day_counter = output_store_day_counter + 1 output_ndays_per_month = output_ndays_per_month + 1 end subroutine output_day subroutine output_month(year, month) ! Take monthly averages and write to csv/bin, if interval is monthly. integer, intent(in) :: year, month integer :: sti(3, 3), stb(3, 3) integer i, nd ! Format month label for csv output write(output_month_ix(output_store_month_counter), "(i0.4, '-', i0.2)") year, month ! Take average of monthly values as needed nd = output_ndays_per_month do i = 1, output_nvars ! 0/1 (space, time) for current variable stb = output_is_requested(i, :, :) ! indeces in aggregation arrays of current variable sti = sum(output_is_requested(:i, :, :), dim=1) if (output_var_state(i)) then if (stb(1, 2) == 1) output_agg_hydrotope_monthly(sti(1, 2), :, :) = & output_agg_hydrotope_monthly(sti(1, 2), :, :) / nd if (stb(2, 2) == 1) output_agg_subbasin_monthly(sti(2, 2), :, :) = & output_agg_subbasin_monthly(sti(2, 2), :, :) / nd if (stb(3, 2) == 1) output_agg_catchment_monthly(sti(3, 2), :, :) = & output_agg_catchment_monthly(sti(3, 2), :, :) / nd end if end do ! Increase/reset counters and write out if necessary if (output_write_interval == "D") then call output_write_monthly else if (output_write_interval == "M") then call output_write_daily call output_write_monthly output_store_day_counter = 1 else if (output_write_interval == "Y") then output_store_month_counter = output_store_month_counter + 1 end if output_ndays_per_year = output_ndays_per_year + output_ndays_per_month output_ndays_per_month = 0 end subroutine output_month subroutine output_year ! Take annual averages and write to csv/bin, if interval is annual. integer :: sti(3, 3), stb(3, 3) integer i, nd ! Take average of annual values as needed nd = output_ndays_per_year do i = 1, output_nvars ! 0/1 (space, time) for current variable stb = output_is_requested(i, :, :) ! indeces in aggregation arrays of current variable sti = sum(output_is_requested(:i, :, :), dim=1) if (output_var_state(i)) then if (stb(1, 3) == 1) output_agg_hydrotope_annual(sti(1, 3), :, :) = & output_agg_hydrotope_annual(sti(1, 3), :, :) / nd if (stb(2, 3) == 1) output_agg_subbasin_annual(sti(2, 3), :, :) = & output_agg_subbasin_annual(sti(2, 3), :, :) / nd if (stb(3, 3) == 1) output_agg_catchment_annual(sti(3, 3), :, :) = & output_agg_catchment_annual(sti(3, 3), :, :) / nd end if end do ! Write out before reset counters if (output_write_interval == "Y") then call output_write_daily call output_write_monthly output_store_day_counter = 1 output_store_month_counter = 1 end if call output_write_annual output_ndays_per_year = 0 end subroutine output_year subroutine output_write_daily ! Write out storage/agg arrays as requested and reset to 0 integer nd ! For varying year and month lengths nd = output_store_day_counter ! Undo day increase of previous output_day call if (output_write_interval /= "D") nd = nd - 1 call output_write_time(1, output_day_ix(:nd), & output_storage_hydrotope(output_hydrotope_requested_ix, :, :nd), & output_agg_subbasin_daily(:, :, :nd), & output_agg_catchment_daily(:, :, :nd)) ! Set all storage arrays to 0 output_storage_hydrotope = 0 output_storage_subbasin = 0 output_agg_subbasin_daily = 0 output_agg_catchment_daily = 0 end subroutine output_write_daily subroutine output_write_monthly ! Write out storage/agg monthly arrays CSV output call output_write_time(2, output_month_ix, output_agg_hydrotope_monthly, & output_agg_subbasin_monthly, output_agg_catchment_monthly) ! Set all storage arrays to 0 output_agg_hydrotope_monthly = 0 output_agg_subbasin_monthly = 0 output_agg_catchment_monthly = 0 end subroutine output_write_monthly subroutine output_write_annual ! Write out storage/agg annual arrays CSV output call output_write_time(3, (/output_month_ix(1)(:4)/), output_agg_hydrotope_annual, & output_agg_subbasin_annual, output_agg_catchment_annual) ! Set all storage arrays to 0 output_agg_hydrotope_annual = 0 output_agg_subbasin_annual = 0 output_agg_catchment_annual = 0 end subroutine output_write_annual subroutine output_write_time(timeix, time_label, hydrotope, subbasin, catchment) ! Abstraction of write out storage/agg arrays integer, intent(in) :: timeix character(len=*), dimension(:), intent(in) :: time_label real, dimension(:, :, :), intent(in) :: hydrotope, subbasin, catchment type(output_file) :: fl integer i, t, nv, ntm, osix(output_nvars) ntm = size(time_label) do i = 1, output_nfiles fl = output_files(i) ! Only continue with files of the timeix frequency if (fl%time /= timeix) cycle ! Abbreviations nv = fl%nvars osix = output_storage_index(i, :) ! CSV output if (fl%format == 1) then select case (fl%space) case (1) call output_write_space_time_csv(fl%file_id, hydrotope(osix(:nv), :, :), & output_hydrotope_id, time_label) case (2) call output_write_space_time_csv(fl%file_id, subbasin(osix(:nv), :, :), & output_subbasin_id, time_label) case (3) call output_write_space_time_csv(fl%file_id, catchment(osix(:nv), :, :), & output_catchment_id, time_label) case (4) call output_write_space_time_csv(fl%file_id, hydrotope(osix(:nv), & output_hydrotope_label_ix, :), time_label=time_label, & space_label=output_hydrotope_label(output_hydrotope_label_ix)) case (5) call output_write_space_time_csv(fl%file_id, subbasin(osix(:nv), & output_subbasin_label_ix, :), time_label=time_label, & space_label=output_subbasin_label(output_subbasin_label_ix)) end select ! Binary output else if (fl%format == 2) then select case (fl%space) case (1) do t = 1, ntm call output_array_to_bin(fl%file_id, hydrotope(osix(:nv), :, t), fl%current_record) output_files(i)%current_record = output_files(i)%current_record + 1 end do case (2) do t = 1, ntm call output_array_to_bin(fl%file_id, subbasin(osix(:nv), :, t), fl%current_record) output_files(i)%current_record = output_files(i)%current_record + 1 end do case (3) do t = 1, ntm call output_array_to_bin(fl%file_id, catchment(osix(:nv), :, t), fl%current_record) output_files(i)%current_record = output_files(i)%current_record + 1 end do case (4) do t = 1, ntm call output_array_to_bin(fl%file_id, hydrotope(osix(:nv), & output_hydrotope_label_ix, t), fl%current_record) output_files(i)%current_record = output_files(i)%current_record + 1 end do case (5) do t = 1, ntm call output_array_to_bin(fl%file_id, subbasin(osix(:nv), & output_subbasin_label_ix, t), fl%current_record) output_files(i)%current_record = output_files(i)%current_record + 1 end do end select end if end do end subroutine output_write_time subroutine output_close integer fi ! Close files do fi = 1, output_nfiles close(output_files(fi)%file_id) end do ! TODO are deallocate arrays really necessary? end subroutine output_close subroutine output_write_space_time_csv(file_id, array, space_ix, time_label, & space_label) ! Write a 3D (vars, space, time) array with time, space index to tidy csv integer, intent(in) :: file_id integer, intent(in), optional :: space_ix(:) character(len=*), intent(in) :: time_label(:) character(len=*), intent(in), optional :: space_label(:) real, intent(in) :: array(:, :, :) integer it if (present(space_ix)) then do it = 1, size(time_label) call output_array_to_csv(file_id, array(:, :, it), space_ix, time_label(it:it)) end do else if (present(space_label)) then do it = 1, size(time_label) call output_array_to_csv(file_id, array(:, :, it), & time_label=time_label(it:it), space_label=space_label) end do else call log_error("output_write_space_time_csv", & "space_ix or space_label have to be given.") end if end subroutine output_write_space_time_csv subroutine output_array_to_csv(file_id, array, space_ix, time_label, space_label) ! Write 2D real array to csv in file_id with a space and time index ! Both space_ix and time_label must be either the same length as array or of length 1 ! to repeat over all rows. integer, intent(in) :: file_id real, intent(in) :: array(:, :) integer, intent(in), optional :: space_ix(:) character(len=*), intent(in) :: time_label(:) character(len=*), intent(in), optional :: space_label(:) character(len=identifier_max_length) :: tsfmt integer :: sh(2), i, ii, lens, lent ! Check input if (rank(array) /= 2) call log_error("output_array_to_csv", & "output_array_to_csv array is not 2D.") sh = shape(array) lens = 1 if (present(space_ix)) then lens = size(space_ix) if (lens /= sh(2) .and. lens /= 1) call log_error("output_array_to_csv", & "space_ix not the right length.") else if (present(space_label)) then lens = size(space_label) if (lens /= sh(2) .and. lens /= 1) & call log_error("output_array_to_csv", "space_label not the right length.") else call log_error("output_array_to_csv", "space_ix or space_label have to be given.") end if lent = size(time_label) if (lent /= sh(2) .and. lent /= 1) & call log_error("output_array_to_csv", "time_label not the right length or 1.") ! Time, space format tsfmt = "(a, ',', "//trim(output_space_index_format)//", ',')" do i = 1, sh(2) ! Time, space indeces if (present(space_ix)) then write(file_id, tsfmt, advance="no") time_label(min(i, lent)), & space_ix(min(i, lens)) else write(file_id, "(a, ', ', a, ', ')", advance="no") time_label(min(i, lent)), & trim(space_label(min(i, lens))) end if do ii = 1, sh(1)-1 write(file_id, "("//trim(output_float_format)//", ',')", advance="no") & array(ii, i) end do ! Last column without , write(file_id, "("//trim(output_float_format)//")") array(sh(1), i) end do end subroutine output_array_to_csv subroutine output_array_to_bin(file_id, array, record) ! Write array to unformatted (binary) output in file_id integer, intent(in) :: file_id, record real, intent(in) :: array(:, :) write(file_id, rec=record) array end subroutine output_array_to_bin subroutine output_print_variables ! Print registered subbasin/hydrotope output variables use utilities, only : string_index write(*, "(A)") output_registered_vars(:string_index("", output_registered_vars)-1) stop end subroutine output_print_variables end module output