Line data Source code
1 : !> \file mo_mrm_write.f90
2 : !> \brief \copybrief mo_mrm_write
3 : !> \details \copydetails mo_mrm_write
4 :
5 : !> \brief write of discharge and restart files
6 : !> \details This module contains the subroutines for writing the discharge files and optionally the restart files.
7 : !> \authors Stephan Thober
8 : !> \date Aug 2015
9 : !> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
10 : !! mHM is released under the LGPLv3+ license \license_note
11 : !> \ingroup f_mrm
12 : module mo_mrm_write
13 :
14 : use mo_kind, only : i4, dp
15 : use mo_mrm_write_fluxes_states, only : OutputDataset
16 : use mo_mrm_global_variables, only : output_time_reference_mrm
17 : use mo_message, only : message, error_message
18 :
19 : implicit none
20 :
21 : public :: mrm_write
22 : public :: mrm_write_optinamelist
23 : public :: mrm_write_optifile
24 :
25 : private
26 :
27 : contains
28 :
29 : ! ------------------------------------------------------------------
30 :
31 : ! NAME
32 : ! mrm_write
33 :
34 : ! PURPOSE
35 : !> \brief write discharge and restart files
36 :
37 : !> \details First, this subroutine calls the writing or restart files that only
38 : !> succeeds if it happens after the write of mHM restart files because
39 : !> mHM restart files must exist. Second, simulated discharge is aggregated to the daily
40 : !> scale and then written to file jointly with observed discharge
41 :
42 : ! HISTORY
43 : !> \authors Juliane Mai, Rohini Kumar & Stephan Thober
44 :
45 : !> \date Aug 2015
46 :
47 : ! Modifications:
48 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
49 : ! Pallav Shrestha, Husain Najafi Mar 2022 - refactoring for measured timestep output
50 :
51 13 : subroutine mrm_write
52 :
53 : use mo_common_constants, only : nodata_dp
54 : use mo_common_mhm_mrm_variables, only : evalPer, mrm_coupling_mode, nTstepDay, simPer, warmingDays
55 : use mo_common_variables, only : mrmFileRestartOut, domainMeta, write_restart
56 : use mo_mrm_global_variables, only : domain_mrm, &
57 : gauge, mRM_runoff, nGaugesTotal, nMeasPerDay
58 : use mo_mrm_restart, only : mrm_write_restart
59 :
60 : implicit none
61 :
62 : integer(i4) :: domainID, iDomain
63 :
64 : integer(i4) :: iDay, iSubDay, iS, iE
65 :
66 : integer(i4) :: maxDailyTimeSteps, maxMeasTimeSteps
67 :
68 : integer(i4) :: tt
69 :
70 : integer(i4) :: gg
71 :
72 : integer(i4) :: nTimeSteps
73 :
74 13 : real(dp), dimension(:, :), allocatable :: d_Qmod, subd_Qmod ! Sim discharge at daily and subdaily time step
75 13 : real(dp), dimension(:, :), allocatable :: d_Qobs, subd_Qobs ! Obs discharge at daily and subdaily time step
76 :
77 : ! between simulated and measured time scale
78 : integer(i4) :: factor
79 :
80 : ! simulated Timesteps per Day
81 : integer(i4) :: TPD_sim
82 :
83 : ! observed Timesteps per Day
84 : integer(i4) :: TPD_obs
85 :
86 : ! --------------------------------------------------------------------------
87 : ! WRITE CONFIG FILE
88 : ! --------------------------------------------------------------------------
89 0 : if (mrm_coupling_mode .eq. 0_i4) call write_configfile()
90 :
91 : ! --------------------------------------------------------------------------
92 : ! WRITE RESTART
93 : ! --------------------------------------------------------------------------
94 13 : if (write_restart) then
95 38 : do iDomain = 1, domainMeta%nDomains
96 25 : domainID = domainMeta%indices(iDomain)
97 38 : if (domainMeta%doRouting(iDomain)) call mrm_write_restart(iDomain, domainID, mrmFileRestartOut)
98 : end do
99 : end if
100 :
101 : ! --------------------------------------------------------------------------
102 : ! STORE DAILY DISCHARGE TIMESERIES OF EACH GAUGING STATION
103 : ! FOR SIMULATIONS DURING THE EVALUATION PERIOD
104 : !
105 : ! **** AT DAILY TIME STEPS ****
106 : ! Note:: Observed Q are stored only for the evaluation period and not for
107 : ! the warming days
108 : ! --------------------------------------------------------------------------
109 :
110 13 : if (nGaugesTotal > 0 .and. nMeasPerDay > 0) then
111 : ! copy time resolution to local variables
112 13 : TPD_sim = nTstepDay
113 13 : TPD_obs = nMeasPerDay
114 :
115 : ! check if modelled timestep is an integer multiple of measured timesteps
116 13 : if (modulo(TPD_sim, TPD_obs) .eq. 0) then
117 13 : factor = TPD_sim / TPD_obs
118 : else
119 0 : call error_message(' Error: Number of modelled datapoints is no multiple of measured datapoints per day')
120 : end if
121 :
122 :
123 51 : maxDailyTimeSteps = maxval(evalPer(1 : domainMeta%nDomains)%julEnd - evalPer(1 : domainMeta%nDomains)%julStart + 1)
124 38 : maxMeasTimeSteps = maxval(evalPer(1 : domainMeta%nDomains)%julEnd - evalPer(1 : domainMeta%nDomains)%julStart + 1) * TPD_obs
125 52 : allocate(d_Qmod (maxDailyTimeSteps, nGaugesTotal))
126 39 : allocate(d_Qobs (maxDailyTimeSteps, nGaugesTotal))
127 52 : allocate(subd_Qmod ( maxMeasTimeSteps, nGaugesTotal))
128 39 : allocate(subd_Qobs ( maxMeasTimeSteps, nGaugesTotal))
129 12101 : d_Qmod = nodata_dp
130 12101 : d_Qobs = nodata_dp
131 12101 : subd_Qmod = nodata_dp
132 12101 : subd_Qobs = nodata_dp
133 :
134 :
135 : ! loop over domains
136 38 : do iDomain = 1, domainMeta%nDomains
137 38 : if (domainMeta%doRouting(iDomain)) then
138 20 : domainID = domainMeta%indices(iDomain)
139 :
140 : ! Convert simulated values to daily
141 20 : nTimeSteps = (simPer(iDomain)%julEnd - simPer(iDomain)%julStart + 1) * nTstepDay
142 20 : iDay = 0
143 8786 : do tt = warmingDays(iDomain) * nTstepDay + 1, nTimeSteps, nTstepDay
144 8766 : iS = tt
145 8766 : iE = tt + nTstepDay - 1
146 8766 : iDay = iDay + 1
147 : ! over gauges
148 18282 : do gg = 1, domain_mrm(iDomain)%nGauges
149 : ! simulation
150 18992 : d_Qmod(iDay, domain_mrm(iDomain)%gaugeIndexList(gg)) = &
151 265158 : sum(mRM_runoff(iS : iE, domain_mrm(iDomain)%gaugeIndexList(gg))) / real(nTstepDay, dp)
152 : end do
153 : end do
154 :
155 20 : dailycheck: if (nMeasPerDay > 1) then
156 : ! Convert observed values to daily
157 0 : nTimeSteps = (simPer(iDomain)%julEnd - simPer(iDomain)%julStart + 1) * nMeasPerDay
158 0 : iDay = 0
159 0 : do tt = 1, nTimeSteps, nMeasPerDay
160 0 : iS = tt
161 0 : iE = tt + nMeasPerDay - 1
162 0 : iDay = iDay + 1
163 : ! over gauges
164 0 : do gg = 1, domain_mrm(iDomain)%nGauges
165 : ! when -9999 value/s are present in current day, daily value remain -9999.
166 0 : if (.not.(any(gauge%Q(iS : iE, domain_mrm(iDomain)%gaugeIndexList(gg)) == nodata_dp))) then
167 : ! observation
168 0 : d_Qobs(iDay, domain_mrm(iDomain)%gaugeIndexList(gg)) = &
169 0 : sum( gauge%Q(iS : iE, domain_mrm(iDomain)%gaugeIndexList(gg))) / real(nMeasPerDay, dp)
170 : end if
171 : end do
172 : end do
173 : else
174 : ! observed values are already at daily (nMeasPerDay = 1) and stored for evalper
175 : ! over gauges
176 42 : do gg = 1, domain_mrm(iDomain)%nGauges
177 : ! observation
178 11363 : d_Qobs(:, domain_mrm(iDomain)%gaugeIndexList(gg)) = gauge%Q(:, domain_mrm(iDomain)%gaugeIndexList(gg))
179 : end do
180 : end if dailycheck
181 :
182 :
183 20 : subdailycheck: if (nMeasPerDay > 1) then
184 :
185 : ! Convert simulated values to subdaily
186 0 : nTimeSteps = (simPer(iDomain)%julEnd - simPer(iDomain)%julStart + 1) * nTstepDay
187 0 : iSubDay = 0
188 0 : do tt = warmingDays(iDomain) * nTstepDay + 1, nTimeSteps, factor
189 0 : iS = tt
190 0 : iE = tt + factor - 1
191 0 : iSubDay = iSubDay + 1
192 : ! over gauges
193 0 : do gg = 1, domain_mrm(iDomain)%nGauges
194 : ! simulation
195 0 : subd_Qmod(iSubDay, domain_mrm(iDomain)%gaugeIndexList(gg)) = &
196 0 : sum(mRM_runoff(iS : iE, domain_mrm(iDomain)%gaugeIndexList(gg))) / real(factor, dp)
197 : end do
198 : end do
199 :
200 : ! Convert observed values to subdaily
201 :
202 : ! observed values are already at subdaily (nMeasPerDay) and stored for evalper
203 : ! over gauges
204 0 : do gg = 1, domain_mrm(iDomain)%nGauges
205 : ! observation
206 0 : subd_Qobs(:, domain_mrm(iDomain)%gaugeIndexList(gg)) = gauge%Q(:, domain_mrm(iDomain)%gaugeIndexList(gg))
207 : end do
208 :
209 : end if subdailycheck
210 :
211 : end if
212 : end do
213 :
214 :
215 : ! write in an ASCII file ! OBS[nModeling_days X nGauges_total] , SIM[nModeling_days X nGauges_total]
216 : ! ToDo: is this if statement reasonable
217 13 : if (allocated(gauge%Q)) call write_daily_obs_sim_discharge(d_Qobs(:, :), d_Qmod(:, :))
218 :
219 : ! write in an ASCII file ! OBS[nMeasPerDay X nGauges_total] , SIM[nMeasPerDay X nGauges_total]
220 13 : if (nMeasPerDay > 1 .and. allocated(gauge%Q)) call write_subdaily_obs_sim_discharge(subd_Qobs(:, :), subd_Qmod(:, :), factor)
221 : ! The subdaily routine is only called if subdaily Q data is provided
222 :
223 : ! free space
224 13 : deallocate(d_Qmod, d_Qobs, subd_Qmod, subd_Qobs)
225 : end if
226 13 : end subroutine mrm_write
227 : ! ------------------------------------------------------------------
228 :
229 : ! NAME
230 : ! write_configfile
231 :
232 : ! PURPOSE
233 : !> \brief This modules writes the results of the configuration into an ASCII-file
234 :
235 : !> \details TODO: add description
236 :
237 : ! HISTORY
238 : !> \authors Christoph Schneider
239 :
240 : !> \date May 2013
241 :
242 : ! Modifications:
243 : ! Juliane Mai May 2013 - module version and documentation
244 : ! Stephan Thober Jun 2014 - bug fix in L11 config print out
245 : ! Stephan Thober Jun 2014 - updated read_restart
246 : ! Rohini, Luis Jul 2015 - updated version, L1 level prints
247 : ! Stephan Thober Sep 2015 - updated write of stream network
248 : ! Stephan Thober Nov 2016 - adapted write to selected case for routing process
249 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
250 :
251 0 : Subroutine write_configfile
252 :
253 13 : use mo_common_constants, only : nodata_dp
254 : use mo_common_file, only : file_config, uconfig
255 : use mo_common_mHM_mRM_variables, only : LCyearId, SimPer, evalPer, mrm_coupling_mode, read_restart, &
256 : resolutionRouting, timeStep, warmPer
257 : use mo_common_variables, only : LC_year_end, LC_year_start, LCfilename, &
258 : dirConfigOut, dirLCover, dirMorpho, dirOut, mrmFileRestartOut, global_parameters, &
259 : global_parameters_name, level0, level1, domainMeta, nLCoverScene, processMatrix, &
260 : resolutionHydrology, write_restart
261 : use mo_kind, only : dp, i4
262 : use mo_mrm_file, only : version
263 : use mo_mrm_global_variables, only : InflowGauge, L11_L1_Id, L11_fromN, L11_label, &
264 : L11_length, L11_netPerm, L11_rOrder, L11_slope, L11_toN, L1_L11_Id, domain_mrm, &
265 : dirGauges, dirTotalRunoff, gauge, level11, nGaugesTotal, nInflowGaugesTotal
266 : use mo_string_utils, only : num2str
267 : use mo_utils, only : ge
268 : use mo_os, only : check_path_isdir
269 :
270 : implicit none
271 :
272 : character(256) :: fName
273 :
274 : integer(i4) :: i, iDomain, domainID, j
275 :
276 : integer(i4) :: err
277 :
278 :
279 0 : fName = trim(adjustl(dirConfigOut)) // trim(adjustl(file_config))
280 0 : call message()
281 0 : call message(' Log-file written to ', trim(fName))
282 : !checking whether the directory exists where the file shall be created or opened
283 0 : call check_path_isdir(trim(adjustl(dirConfigOut)), raise=.true.)
284 0 : open(uconfig, file = fName, status = 'unknown', action = 'write', iostat = err)
285 0 : if (err .ne. 0) then
286 0 : call error_message(' Problems while creating File. Error-Code ', num2str(err))
287 : end if
288 0 : write(uconfig, 200)
289 0 : write(uconfig, 100) 'mRM-UFZ v-' // trim(version)
290 0 : write(uconfig, 100) 'S. Thober, L. Samaniego & R. Kumar, UFZ'
291 0 : write(uconfig, 200)
292 0 : write(uconfig, 100)
293 0 : write(uconfig, 201) ' M A I N mRM C O N F I G U R A T I O N I N F O R M A T I O N '
294 0 : write(uconfig, 100)
295 0 : write(uconfig, 103) 'Number of domains ', domainMeta%nDomains
296 0 : write(uconfig, 103) 'Total No. of gauges ', nGaugesTotal
297 0 : write(uconfig, 103) 'Time Step [h] ', timeStep
298 0 : do iDomain = 1, domainMeta%nDomains
299 0 : domainID = domainMeta%indices(iDomain)
300 0 : write(uconfig, 103) 'Total No. of nodes ', level11(iDomain)%nCells
301 0 : write(uconfig, 103) 'No. of cells L0 ', level0(domainMeta%L0DataFrom(iDomain))%nCells
302 0 : write(uconfig, 103) 'No. of cells L1 ', level1(iDomain)%nCells
303 0 : write(uconfig, 103) 'No. of cells L11 ', level11(iDomain)%nCells
304 :
305 : ! select case (iFlag_cordinate_sys)
306 : ! case (0)
307 0 : write(uconfig, 301) 'domain ', domainID, ' Hydrology Resolution [m] ', resolutionHydrology(iDomain)
308 0 : write(uconfig, 301) 'domain ', domainID, ' Routing Resolution [m] ', resolutionRouting(iDomain)
309 : ! case(1)
310 : ! write(uconfig, 302) 'domain ',domainID, ' Hydrology Resolution [o] ', resolutionHydrology(iDomain)
311 : ! write(uconfig, 302) 'domain ',domainID, ' Routing Resolution [o] ', resolutionRouting(iDomain)
312 : ! end select
313 : end do
314 0 : write(uconfig, 126) 'Flag READ restart ', read_restart
315 0 : write(uconfig, 126) 'Flag WRITE restart ', write_restart
316 : !
317 : !******************
318 : ! Model Run period
319 : !******************
320 0 : do iDomain = 1, domainMeta%nDomains
321 0 : domainID = domainMeta%indices(iDomain)
322 0 : write(uconfig, 115) ' Model Run Periods for domain ', num2str(domainID)
323 : write(uconfig, 116) &
324 0 : 'From To', &
325 0 : ' Day Month Year Day Month Year'
326 : write(uconfig, 117) &
327 0 : 'Warming Period (1) ', &
328 0 : warmPer(iDomain)%dStart, warmPer(iDomain)%mStart, warmPer(iDomain)%yStart, &
329 0 : warmPer(iDomain)%dEnd, warmPer(iDomain)%mEnd, warmPer(iDomain)%yEnd
330 : write(uconfig, 117) &
331 0 : 'Evaluation Period (2) ', &
332 0 : evalPer(iDomain)%dStart, evalPer(iDomain)%mStart, evalPer(iDomain)%yStart, &
333 0 : evalPer(iDomain)%dEnd, evalPer(iDomain)%mEnd, evalPer(iDomain)%yEnd
334 : write(uconfig, 117) &
335 0 : 'Simulation Period (1)+(2) ', &
336 0 : SimPer(iDomain)%dStart, SimPer(iDomain)%mStart, SimPer(iDomain)%yStart, &
337 0 : SimPer(iDomain)%dEnd, SimPer(iDomain)%mEnd, SimPer(iDomain)%yEnd
338 : end do
339 :
340 : !*********************************
341 : ! Model Land Cover Observations
342 : !*********************************
343 0 : if (processMatrix(8, 1) .eq. 1) then
344 0 : do iDomain = 1, domainMeta%nDomains
345 0 : domainID = domainMeta%indices(iDomain)
346 0 : write(uconfig, 118) ' Land Cover Observations for domain ', num2str(domainID)
347 0 : write(uconfig, 119) ' Start Year', ' End Year', ' Land cover scene', 'Land Cover File'
348 0 : do i = 1, nLCoverScene
349 0 : write(uconfig, 120) LC_year_start(i), LC_year_end(i), &
350 0 : LCyearId(max(evalPer(iDomain)%yStart, LC_year_start(i)), iDomain), trim(LCfilename(i))
351 : end do
352 : end do
353 : end if
354 : !*********************************
355 : ! Initial Parameter Ranges
356 : !*********************************
357 0 : write(uconfig, 121) ' Initial Transfer Function Parameter Ranges (gammas) '
358 : !
359 : ! Transfer functions
360 : write(uconfig, 122) &
361 0 : ' i', ' min', ' max', ' current', &
362 0 : ' name'
363 0 : do i = 1, size(global_parameters, 1)
364 : write(uconfig, 123) &
365 0 : i, global_parameters(i, 1), global_parameters(i, 2), global_parameters(i, 3), &
366 0 : trim(adjustl(global_parameters_name(i)))
367 : end do
368 : ! domain runoff data
369 0 : write(uconfig, 202) ' domain Runoff Data '
370 0 : write(uconfig, 107) ' Gauge No.', ' domain Id', ' Qmax[m3/s]', ' Qmin[m3/s]'
371 0 : do i = 1, nGaugesTotal
372 0 : if(any(gauge%Q(:, i) > nodata_dp)) then
373 0 : write(uconfig, 108) i, gauge%domainId(i), maxval(gauge%Q(:, i), gauge%Q(:, i) > nodata_dp), &
374 0 : minval(gauge%Q(:, i), gauge%Q(:, i) > nodata_dp)
375 : else
376 0 : write(uconfig, 108) i, gauge%domainId(i), nodata_dp, nodata_dp
377 : end if
378 : end do
379 : ! inflow gauge data
380 0 : if (nInflowGaugesTotal .GT. 0) then
381 0 : write(uconfig, 202) ' domain Inflow Data '
382 0 : write(uconfig, 107) ' Gauge No.', ' domain Id', ' Qmax[m3/s]', ' Qmin[m3/s]'
383 0 : do i = 1, nInflowGaugesTotal
384 0 : if(all(InflowGauge%Q(:, i) > nodata_dp)) then
385 0 : write(uconfig, 108) i, InflowGauge%domainId(i), maxval(InflowGauge%Q(:, i), InflowGauge%Q(:, i) > nodata_dp), &
386 0 : minval(InflowGauge%Q(:, i), InflowGauge%Q(:, i) > nodata_dp)
387 : else
388 0 : write(uconfig, 108) i, InflowGauge%domainId(i), nodata_dp, nodata_dp
389 : end if
390 : end do
391 : end if
392 : ! domain config
393 0 : write(uconfig, 218) 'domain-wise Configuration'
394 0 : do iDomain = 1, domainMeta%nDomains
395 0 : domainID = domainMeta%indices(iDomain)
396 0 : write(uconfig, 103) 'domain No. ', domainID, &
397 0 : 'No. of gauges ', domain_mrm(iDomain)%nGauges
398 :
399 0 : write(uconfig, 222) 'Directory list'
400 :
401 0 : write(uconfig, 224) 'Directory to morphological input ', dirMorpho(iDomain)
402 0 : write(uconfig, 224) 'Directory to land cover input ', dirLCover(iDomain)
403 0 : write(uconfig, 224) 'Directory to gauging station input ', dirGauges(iDomain)
404 0 : if (mrm_coupling_mode .eq. 0) then
405 0 : write(uconfig, 224) 'Directory to simulated runoff input ', dirTotalRunoff(iDomain)
406 : end if
407 0 : write(uconfig, 224) 'Directory to write output by default ', dirOut(iDomain)
408 0 : write(uconfig, 224) 'File to write mRM output when restarted ', mrmFileRestartOut(iDomain)
409 :
410 0 : write(uconfig, 102) 'River Network (Routing level)'
411 0 : write(uconfig, 100) 'Label 0 = intermediate draining cell '
412 0 : write(uconfig, 100) 'Label 1 = headwater cell '
413 0 : write(uconfig, 100) 'Label 2 = sink cell '
414 :
415 0 : if (processMatrix(8, 1) .eq. 1_i4) then
416 0 : write(uconfig, 104) ' Overall', &
417 0 : ' From', &
418 0 : ' To', &
419 0 : ' Routing', &
420 0 : ' Label', &
421 0 : ' Length', &
422 0 : ' Mean', &
423 0 : ' Link', &
424 0 : ' Routing', &
425 0 : ' Routing', &
426 0 : ' Sequence', &
427 0 : ' ', &
428 0 : ' ', &
429 0 : ' Slope'
430 : !
431 0 : write(uconfig, 105) ' Id', &
432 0 : ' Node', &
433 0 : ' Node', &
434 0 : '', &
435 0 : '', &
436 0 : ' [km]', &
437 0 : ' [o/oo]'
438 : !
439 0 : do j = level11(iDomain)%iStart, level11(iDomain)%iEnd - 1
440 0 : i = L11_netPerm(j) + level11(iDomain)%iStart - 1 ! adjust permutation for multi-domain option
441 0 : write(uconfig, 106) i, L11_fromN(i), L11_toN(i), L11_rOrder(i), L11_label(i), &
442 0 : L11_length(i) / 1000.0_dp, L11_slope(i) * 1.0e3_dp
443 : end do
444 :
445 0 : else if (processMatrix(8, 1) .eq. 2_i4) then
446 0 : write(uconfig, 134) ' Overall', &
447 0 : ' From', &
448 0 : ' To', &
449 0 : ' Routing', &
450 0 : ' Label', &
451 0 : ' Link', &
452 0 : ' Routing', &
453 0 : ' Routing', &
454 0 : ' Sequence', &
455 0 : ' '
456 : !
457 0 : write(uconfig, 135) ' Id', &
458 0 : ' Node', &
459 0 : ' Node', &
460 0 : '', &
461 0 : ''
462 : !
463 0 : do j = level11(iDomain)%iStart, level11(iDomain)%iEnd - 1
464 0 : i = L11_netPerm(j) + level11(iDomain)%iStart - 1 ! adjust permutation for multi-domain option
465 0 : write(uconfig, 136) i, L11_fromN(i), L11_toN(i), L11_rOrder(i), L11_label(i)
466 : end do
467 : end if
468 : ! draining node at L11
469 0 : write(uconfig, 109) ' Overall', ' domain', &
470 0 : ' Cell', ' Routing', &
471 0 : ' Id', ' Node Id'
472 0 : do i = level11(iDomain)%Id(1), level11(iDomain)%Id(level11(iDomain)%nCells)
473 0 : write(uconfig, 110) i + level11(iDomain)%iStart - 1, i
474 : end do
475 :
476 : ! L1 level information
477 0 : write(uconfig, 111) ' Modeling', ' Routing', ' Effective', &
478 0 : ' Cell', ' Cell Id', ' Area', &
479 0 : ' Id', ' [-]', ' [km2]'
480 0 : if (ge(resolutionRouting(iDomain), resolutionHydrology(iDomain))) then
481 0 : do i = level1(iDomain)%Id(1), level1(iDomain)%Id(level1(iDomain)%nCells)
482 0 : write(uconfig, 113) i + level1(iDomain)%iStart - 1, L1_L11_Id (i + level1(iDomain)%iStart - 1), &
483 0 : level1(iDomain)%CellArea(i) * 1.E-6_dp
484 0 : domainID = domainMeta%indices(iDomain)
485 : end do
486 : else
487 0 : do i = level11(iDomain)%Id(1), level11(iDomain)%Id(level11(iDomain)%nCells)
488 0 : write(uconfig, 110) i + level11(iDomain)%iStart - 1, L11_L1_Id (i + level11(iDomain)%iStart - 1)
489 : end do
490 : end if
491 0 : write(uconfig, 114) ' Total[km2]', sum(level1(iDomain)%CellArea) * 1.E-6_dp
492 : end do
493 :
494 0 : write(uconfig, *)
495 0 : close(uconfig)
496 :
497 : !! Formats
498 : 100 format (a80)
499 : 102 format (/ 30('-') / a30 / 30('-'))
500 : 103 format (a20, 10x, i10)
501 : 104 format (/ 75('-') / 5a10, 5x, 2a10 / 5a10, 5x, 2a10)
502 : 105 format (5a10, 5x, 2a10 / 75('-'))
503 : 106 format (5i10, 5x, 2f10.3)
504 : 107 format (2a10, 2a15)
505 : 108 format (2i10, 2f15.3)
506 : !
507 : 109 format (/ 20('-') / 2a10 / 2a10 / 2a10 / 20('-'))
508 : 110 format (2i10)
509 : !
510 : 111 format (/ 30('-') / 3a10 / 3a10 / 3a10 / 30('-'))
511 : 113 format ( 2i10, 1f10.3 )
512 : 114 format (30('-') / a15, 5x, 1f10.3 /)
513 : !
514 : 115 format (/61('-')/ a50, a10 /61('-'))
515 : 116 format (39x, a22 / 25x, a36)
516 : 117 format (3(a25, 6(i6)))
517 : !
518 : 118 format (/50('-')/ a40, a10 /50('-'))
519 : 119 format (a10, a10, a20, a20/)
520 : 120 format (i10, i10, 10x, i10, a20)
521 : !
522 : 121 format (/55('-')/ a55 /55('-'))
523 : 122 format (a10, 3a15, a35)
524 : 123 format (i10, 3f15.3, a35)
525 : !
526 : 126 format (a30, 9x, L1)
527 : !
528 : 134 format (/ 50('-') / 5a10 / 5a10)
529 : 135 format (5a10 / 50('-'))
530 : 136 format (5i10)
531 : !
532 : 200 format (80('-'))
533 : 201 format (a80)
534 : 202 format (/50('-')/ a50 /50('-'))
535 : !
536 : 218 format (/ 80('-')/ 26x, a24, 26x, /80('-'))
537 : 222 format (/80('-')/ 26x, a21 /80('-'))
538 : 224 format (a40, 5x, a256)
539 :
540 : 301 format (a7, i2, a32, f15.0)
541 : ! 302 format (a7, i2, a32,es20.8)
542 0 : end Subroutine write_configfile
543 :
544 : ! ------------------------------------------------------------------
545 :
546 : ! NAME
547 : ! write_daily_obs_sim_discharge
548 :
549 : ! PURPOSE
550 : !> \brief Write a file for the daily observed and simulated discharge timeseries
551 : !> during the evaluation period for each gauging station
552 :
553 : !> \details Write a file for the daily observed and simulated discharge timeseries
554 : !> during the evaluation period for each gauging station
555 :
556 : ! INTENT(IN)
557 : !> \param[in] "real(dp), dimension(:, :) :: Qobs" daily time series of observed dischargedims = (nModeling_days
558 : !> , nGauges_total)
559 : !> \param[in] "real(dp), dimension(:, :) :: Qsim" daily time series of modeled dischargedims = (nModeling_days ,
560 : !> nGauges_total)
561 :
562 : ! HISTORY
563 : !> \authors Rohini Kumar
564 :
565 : !> \date August 2013
566 :
567 : ! Modifications:
568 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
569 :
570 26 : subroutine write_daily_obs_sim_discharge(Qobs, Qsim)
571 :
572 0 : use mo_common_constants, only : nodata_dp
573 : use mo_common_mhm_mrm_variables, only : evalPer
574 : use mo_common_variables, only : dirOut, domainMeta
575 : use mo_errormeasures, only : kge, nse
576 : use mo_julian, only : dec2date
577 : use mo_mrm_file, only : file_daily_discharge, ncfile_discharge, udaily_discharge
578 : use mo_mrm_global_variables, only : domain_mrm, gauge
579 : use mo_string_utils, only : num2str
580 : use mo_utils, only : ge
581 : use mo_netcdf, only : NcDataset, NcDimension, NcVariable
582 :
583 : implicit none
584 :
585 : ! daily time series of observed dischargedims = (nModeling_days , nGauges_total)
586 : real(dp), dimension(:, :), intent(in) :: Qobs
587 :
588 : ! daily time series of modeled dischargedims = (nModeling_days , nGauges_total)
589 : real(dp), dimension(:, :), intent(in) :: Qsim
590 :
591 : character(256) :: fName, formHeader, formData, dummy
592 :
593 : integer(i4) :: domainID, iDomain, gg, tt, err
594 :
595 : integer(i4) :: igauge_start, igauge_end
596 :
597 : integer(i4) :: day, month, year
598 :
599 : integer(i4) :: tlength
600 :
601 : ! time axis
602 13 : integer(i4), allocatable, dimension(:) :: taxis
603 :
604 13 : real(dp) :: newTime
605 :
606 : ! nc related variables
607 : type(NcDataset) :: nc_out
608 : type(NcDimension) :: dim, dim_bnd
609 : type(NcVariable) :: var
610 :
611 : ! initalize igauge_start
612 13 : igauge_start = 1
613 :
614 : ! domain loop
615 38 : do iDomain = 1, domainMeta%nDomains
616 25 : domainID = domainMeta%indices(iDomain)
617 25 : if(domain_mrm(iDomain)%nGauges .lt. 1) cycle
618 :
619 : ! estimate igauge_end
620 19 : igauge_end = igauge_start + domain_mrm(iDomain)%nGauges - 1
621 :
622 : ! check the existance of file
623 19 : fName = trim(adjustl(dirOut(iDomain))) // trim(adjustl(file_daily_discharge))
624 19 : open(udaily_discharge, file = trim(fName), status = 'unknown', action = 'write', iostat = err)
625 19 : if(err .ne. 0) then
626 0 : call error_message(' IOError while openening "', trim(fName), '". Error-Code ', num2str(err))
627 : end if
628 :
629 : ! header
630 19 : write(formHeader, *) '( 4a8, ', domain_mrm(iDomain)%nGauges, '(2X, a5, i10.10, 2X, a5, i10.10) )'
631 19 : write(udaily_discharge, formHeader) 'No', 'Day', 'Mon', 'Year', &
632 22 : ('Qobs_', gauge%gaugeId(gg), &
633 60 : 'Qsim_', gauge%gaugeId(gg), gg = igauge_start, igauge_end)
634 :
635 : ! form data
636 19 : write(formData, *) '( 4I8, ', domain_mrm(iDomain)%nGauges, '(2X, f15.7 , 2X, f15.7 ) )'
637 :
638 : ! write data
639 19 : newTime = real(evalPer(iDomain)%julStart, dp) - 0.5_dp
640 :
641 8420 : do tt = 1, (evalPer(iDomain)%julEnd - evalPer(iDomain)%julStart + 1)
642 8401 : call dec2date(newTime, yy = year, mm = month, dd = day)
643 17897 : write(udaily_discharge, formData) tt, day, month, year, (Qobs(tt, gg), Qsim(tt, gg), gg = igauge_start, igauge_end)
644 8420 : newTime = newTime + 1.0_dp
645 : end do
646 :
647 : ! close file
648 19 : close(udaily_discharge)
649 :
650 : ! ======================================================================
651 : ! write netcdf file
652 : ! ======================================================================
653 19 : fName = trim(adjustl(dirOut(iDomain))) // trim(adjustl(ncfile_discharge))
654 19 : nc_out = NcDataset(trim(fName), "w")
655 19 : tlength = evalPer(iDomain)%julEnd - evalPer(iDomain)%julStart + 1
656 : ! write time
657 57 : allocate(taxis(tlength))
658 :
659 : ! tt is dependent on the unit of the time axis and is set to hours in mRM
660 : select case( output_time_reference_mrm)
661 : case(0)
662 8420 : forall(tt = 1 : tlength) taxis(tt) = (tt-1) * 24
663 : case(1)
664 0 : forall(tt = 1 : tlength) taxis(tt) = tt * 24 - 12
665 : case(2)
666 19 : forall(tt = 1 : tlength) taxis(tt) = tt * 24
667 : end select
668 :
669 19 : call dec2date(real(evalPer(iDomain)%julStart, dp) - 0.5_dp, yy = year, mm = month, dd = day)
670 19 : dim = nc_out%setDimension("time", tlength)
671 38 : var = nc_out%setVariable("time", "i32", [dim])
672 19 : call var%setData(taxis)
673 : call var%setAttribute( &
674 : "units", &
675 : 'hours since '//trim(num2str(year))//'-'//trim(num2str(month, '(i2.2)'))//'-'//trim(num2str(day, '(i2.2)'))//' 00:00:00' &
676 19 : )
677 19 : call var%setAttribute("long_name", "time in hours")
678 19 : call var%setAttribute("bounds", "time_bnds")
679 19 : call var%setAttribute("axis", "T")
680 19 : dim_bnd = nc_out%setDimension("bnds", 2)
681 57 : var = nc_out%setVariable("time_bnds", "i32", [dim_bnd, dim])
682 8420 : do tt = 1, tlength
683 25203 : call var%setData((tt - 1) * 24, (/1, tt/))
684 25222 : call var%setData(tt * 24, (/2, tt/))
685 : end do
686 19 : deallocate(taxis)
687 : ! write gauges
688 41 : do gg = igauge_start, igauge_end
689 44 : var = nc_out%setVariable('Qsim_' // trim(num2str(gauge%gaugeID(gg), '(i10.10)')), "f64", [dim])
690 22 : call var%setFillValue(nodata_dp)
691 22 : call var%setData(Qsim(1 : tlength, gg))
692 22 : call var%setAttribute("units", "m3 s-1")
693 22 : call var%setAttribute("long_name", 'simulated discharge at gauge ' // trim(num2str(gauge%gaugeID(gg), '(i10.10)')))
694 22 : call var%setAttribute("missing_value", nodata_dp)
695 : ! write observed discharge at that gauge
696 44 : var = nc_out%setVariable('Qobs_' // trim(num2str(gauge%gaugeID(gg), '(i10.10)')), "f64", [dim])
697 22 : call var%setFillValue(nodata_dp)
698 22 : call var%setData(Qobs(1 : tlength, gg))
699 22 : call var%setAttribute("units", "m3 s-1")
700 22 : call var%setAttribute("long_name", 'observed discharge at gauge ' // trim(num2str(gauge%gaugeID(gg), '(i10.10)')))
701 41 : call var%setAttribute("missing_value", nodata_dp)
702 : end do
703 19 : call nc_out%close()
704 :
705 : ! ======================================================================
706 : ! screen output
707 : ! ======================================================================
708 :
709 : ! if ( nMeasPerDay == 1_i4 ) then ! only print daily stats for daily Qobs
710 :
711 19 : call message()
712 19 : write(dummy, '(I3)') domainID
713 19 : call message(' OUTPUT: saved daily discharge file for domain ', trim(adjustl(dummy)))
714 19 : call message(' to ', trim(fname))
715 41 : do gg = igauge_start, igauge_end
716 11362 : if (count(ge(Qobs(:, gg), 0.0_dp)) > 1) then
717 : call message(' KGE of daily discharge (gauge #', trim(adjustl(num2str(gg))), '): ', &
718 10611 : trim(adjustl(num2str(kge(Qobs(:, gg), Qsim(:, gg), mask = (ge(Qobs(:, gg), 0.0_dp)))))))
719 : call message(' NSE of daily discharge (gauge #', trim(adjustl(num2str(gg))), '): ', &
720 10611 : trim(adjustl(num2str(nse(Qobs(:, gg), Qsim(:, gg), mask = (ge(Qobs(:, gg), 0.0_dp)))))))
721 : end if
722 : end do
723 :
724 : ! end if
725 :
726 : ! update igauge_start
727 19 : igauge_start = igauge_end + 1
728 : !
729 : end do
730 : !
731 13 : end subroutine write_daily_obs_sim_discharge
732 :
733 :
734 : ! ------------------------------------------------------------------
735 :
736 : ! NAME
737 : ! write_subdaily_obs_sim_discharge
738 :
739 : ! PURPOSE
740 : !> \brief Write a file for the simulated discharge timeseries
741 : !> during the evaluation period for each gauging station
742 :
743 : !> \details Write a file for the simulated discharge timeseries
744 : !> during the evaluation period for each gauging station
745 :
746 : ! INTENT(IN)
747 : !> \param[in] "real(dp), dimension(:, :) :: Qobs" time series of observed discharge dims = (nMeasTimeSteps ,
748 : !> nGauges_total)
749 : !> \param[in] "real(dp), dimension(:, :) :: Qsim" time series of modeled discharge dims = (nMeasTimeSteps ,
750 : !> nGauges_total)
751 : !> \param[in] "integer(i4), :: factor" ratio of modelled time steps per day to observation time steps per day
752 :
753 : ! HISTORY
754 : !> \authors Rohini Kumar
755 :
756 : !> \date August 2013
757 :
758 : ! Modifications:
759 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
760 : ! Pallav Shrestha Jul 2021 - ported for printing out simulations at hourly time step
761 : ! Pallav Shrestha, Husain Najafi Mar 2022 - refactoring for subdaily timestep output
762 :
763 0 : subroutine write_subdaily_obs_sim_discharge(Qobs, Qsim, factor)
764 :
765 13 : use mo_common_constants, only : nodata_dp
766 : use mo_common_mhm_mrm_variables, only : evalPer
767 : use mo_common_variables, only : dirOut, domainMeta
768 : use mo_errormeasures, only : kge, nse
769 : use mo_julian, only : dec2date
770 : use mo_mrm_file, only : ncfile_subdaily_discharge, file_subdaily_discharge, &
771 : usubdaily_discharge
772 : use mo_mrm_global_variables, only : domain_mrm, gauge, nMeasPerDay
773 : use mo_string_utils, only : num2str
774 : use mo_utils, only : ge
775 : use mo_netcdf, only : NcDataset, NcDimension, NcVariable
776 :
777 : implicit none
778 :
779 : ! time series of observed discharge. dims = (nMeasTimeSteps , nGauges_total)
780 : real(dp), dimension(:, :), intent(in) :: Qobs
781 :
782 : ! time series of modeled discharge. dims = (nMeasTimeSteps , nGauges_total)
783 : real(dp), dimension(:, :), intent(in) :: Qsim
784 :
785 : ! ratio of modelled time steps per day to observation time steps per day
786 : integer(i4), intent(in) :: factor
787 :
788 : character(256) :: fName, formHeader, formData, dummy
789 :
790 : integer(i4) :: domainID, iDomain, gg, tt, err
791 :
792 : integer(i4) :: igauge_start, igauge_end
793 :
794 : integer(i4) :: hour, day, month, year
795 :
796 : integer(i4) :: tlength
797 :
798 : ! time axis
799 0 : integer(i4), allocatable, dimension(:) :: taxis
800 :
801 0 : real(dp) :: newTime
802 :
803 : ! nc related variables
804 : type(NcDataset) :: nc_out
805 : type(NcDimension) :: dim, dim_bnd
806 : type(NcVariable) :: var
807 :
808 : ! use minutes if needed
809 : logical :: use_minutes
810 : integer(i4) :: time_unit_factor
811 :
812 : ! initalize igauge_start
813 0 : igauge_start = 1
814 :
815 : ! domain loop
816 0 : do iDomain = 1, domainMeta%nDomains
817 0 : domainID = domainMeta%indices(iDomain)
818 0 : if(domain_mrm(iDomain)%nGauges .lt. 1) cycle
819 :
820 : ! estimate igauge_end
821 0 : igauge_end = igauge_start + domain_mrm(iDomain)%nGauges - 1
822 :
823 :
824 : ! ======================================================================
825 : ! write text file
826 : ! ======================================================================
827 :
828 : ! check the existance of file
829 0 : fName = trim(adjustl(dirOut(iDomain))) // trim(adjustl(file_subdaily_discharge))
830 0 : open(usubdaily_discharge, file = trim(fName), status = 'unknown', action = 'write', iostat = err)
831 0 : if(err .ne. 0) then
832 0 : call error_message(' IOError while openening "', trim(fName), '". Error-Code ', num2str(err))
833 : end if
834 :
835 : ! header
836 0 : write(formHeader, *) '( 5a8, ', domain_mrm(iDomain)%nGauges, '(2X, a5, i10.10, 2X, a5, i10.10) )'
837 0 : write(usubdaily_discharge, formHeader) 'No', 'Hour', 'Day', 'Mon', 'Year', &
838 0 : ('Qobs_', gauge%gaugeId(gg), &
839 0 : 'Qsim_', gauge%gaugeId(gg), gg = igauge_start, igauge_end)
840 :
841 : ! form data
842 0 : write(formData, *) '( 5I8, ', domain_mrm(iDomain)%nGauges, '(2X, f15.7 , 2X, f15.7 ) )'
843 :
844 : ! write data
845 0 : newTime = real(evalPer(iDomain)%julStart, dp) - 0.5_dp
846 :
847 0 : do tt = 1, (evalPer(iDomain)%julEnd - evalPer(iDomain)%julStart + 1) * nMeasPerDay
848 0 : call dec2date(newTime, yy = year, mm = month, dd = day, hh = hour)
849 0 : write(usubdaily_discharge, formData) tt, hour, day, month, year, (Qobs(tt, gg), Qsim(tt, gg), gg = igauge_start, igauge_end)
850 0 : newTime = newTime + 1.0_dp / real(nMeasPerDay, dp)
851 : end do
852 :
853 : ! close file
854 0 : close(usubdaily_discharge)
855 :
856 : ! ======================================================================
857 : ! write netcdf file
858 : ! ======================================================================
859 0 : fName = trim(adjustl(dirOut(iDomain))) // trim(adjustl(ncfile_subdaily_discharge))
860 0 : nc_out = NcDataset(trim(fName), "w")
861 0 : tlength = (evalPer(iDomain)%julEnd - evalPer(iDomain)%julStart + 1) * nMeasPerDay
862 : ! write time
863 0 : allocate(taxis(tlength))
864 :
865 0 : use_minutes = .false.
866 0 : time_unit_factor = 1
867 0 : if ( mod(factor, 2) == 1 ) then
868 0 : use_minutes = .true.
869 0 : time_unit_factor = 60
870 : end if
871 :
872 : ! tt is dependent on the unit of the time axis and is set to hours in mRM
873 : select case( output_time_reference_mrm)
874 : case(0)
875 0 : forall(tt = 1 : tlength) taxis(tt) = (tt-1) * factor * time_unit_factor
876 : case(1)
877 0 : forall(tt = 1 : tlength) taxis(tt) = tt * factor - factor * time_unit_factor / 2
878 : case(2)
879 0 : forall(tt = 1 : tlength) taxis(tt) = tt * factor * time_unit_factor
880 : end select
881 :
882 0 : call dec2date(real(evalPer(iDomain)%julStart, dp) - 0.5_dp, yy = year, mm = month, dd = day, hh = hour)
883 0 : dim = nc_out%setDimension("time", tlength)
884 0 : var = nc_out%setVariable("time", "i32", [dim])
885 0 : call var%setData(taxis)
886 0 : if (use_minutes) then
887 : call var%setAttribute( &
888 : "units", &
889 : 'minutes since '//trim(num2str(year))//'-'//trim(num2str(month, '(i2.2)'))//'-'//trim(num2str(day, '(i2.2)'))//' '// &
890 : trim(num2str(hour, '(i2.2)'))//':00:00' &
891 0 : )
892 0 : call var%setAttribute("long_name", "time in minutes")
893 : else
894 : call var%setAttribute( &
895 : "units", &
896 : 'hours since '//trim(num2str(year))//'-'//trim(num2str(month, '(i2.2)'))//'-'//trim(num2str(day, '(i2.2)'))//' '// &
897 : trim(num2str(hour, '(i2.2)'))//':00:00' &
898 0 : )
899 0 : call var%setAttribute("long_name", "time in hours")
900 : end if
901 0 : call var%setAttribute("bounds", "time_bnds")
902 0 : call var%setAttribute("axis", "T")
903 0 : dim_bnd = nc_out%setDimension("bnds", 2)
904 0 : var = nc_out%setVariable("time_bnds", "i32", [dim_bnd, dim])
905 0 : do tt = 1, tlength
906 0 : call var%setData((tt - 1) * factor * time_unit_factor, (/1, tt/))
907 0 : call var%setData(tt * factor * time_unit_factor, (/2, tt/))
908 : end do
909 0 : deallocate(taxis)
910 : ! write gauges
911 0 : do gg = igauge_start, igauge_end
912 0 : var = nc_out%setVariable('Qsim_' // trim(num2str(gauge%gaugeID(gg), '(i10.10)')), "f64", [dim])
913 0 : call var%setFillValue(nodata_dp)
914 0 : call var%setData(Qsim(1 : tlength, gg))
915 0 : call var%setAttribute("units", "m3 s-1")
916 0 : call var%setAttribute("long_name", 'simulated discharge at gauge ' // trim(num2str(gauge%gaugeID(gg), '(i10.10)')))
917 0 : call var%setAttribute("missing_value", nodata_dp)
918 : ! write observed discharge at that gauge
919 0 : var = nc_out%setVariable('Qobs_' // trim(num2str(gauge%gaugeID(gg), '(i10.10)')), "f64", [dim])
920 0 : call var%setFillValue(nodata_dp)
921 0 : call var%setData(Qobs(1 : tlength, gg))
922 0 : call var%setAttribute("units", "m3 s-1")
923 0 : call var%setAttribute("long_name", 'observed discharge at gauge ' // trim(num2str(gauge%gaugeID(gg), '(i10.10)')))
924 0 : call var%setAttribute("missing_value", nodata_dp)
925 : end do
926 0 : call nc_out%close()
927 :
928 : ! ======================================================================
929 : ! screen output
930 : ! ======================================================================
931 0 : call message()
932 0 : write(dummy, '(I3)') domainID
933 0 : call message(' OUTPUT: saved subdaily discharge file for domain ', trim(adjustl(dummy)))
934 0 : call message(' to ', trim(fname))
935 0 : do gg = igauge_start, igauge_end
936 0 : if (count(ge(Qobs(:, gg), 0.0_dp)) > 1) then
937 : call message(' KGE of subdaily discharge (gauge #', trim(adjustl(num2str(gg))), '): ', &
938 0 : trim(adjustl(num2str(kge(Qobs(:, gg), Qsim(:, gg), mask = (ge(Qobs(:, gg), 0.0_dp)))))))
939 : call message(' NSE of subdaily discharge (gauge #', trim(adjustl(num2str(gg))), '): ', &
940 0 : trim(adjustl(num2str(nse(Qobs(:, gg), Qsim(:, gg), mask = (ge(Qobs(:, gg), 0.0_dp)))))))
941 : end if
942 : end do
943 :
944 : ! update igauge_start
945 0 : igauge_start = igauge_end + 1
946 : !
947 : end do
948 : !
949 0 : end subroutine write_subdaily_obs_sim_discharge
950 :
951 : ! ------------------------------------------------------------------
952 :
953 : ! NAME
954 : ! mrm_write_optifile
955 :
956 : ! PURPOSE
957 : !> \brief Write briefly final optimization results.
958 :
959 : !> \details Write overall best objective function and the best optimized parameter set to a file_opti.
960 :
961 : ! INTENT(IN)
962 : !> \param[in] "real(dp) :: best_OF" best objective function value as returnedby the
963 : !> optimization routine
964 : !> \param[in] "real(dp), dimension(:) :: best_paramSet" best associated global parameter setCalled only
965 : !> when optimize is .TRUE.
966 : !> \param[in] "character(len = *), dimension(:) :: param_names"
967 :
968 : ! HISTORY
969 : !> \authors David Schaefer
970 :
971 : !> \date July 2013
972 :
973 : ! Modifications:
974 : ! Rohini Kumar Aug 2013 - change in structure of the code including call statements
975 : ! Juliane Mai Oct 2013 - clear parameter names added
976 : ! - double precision written
977 : ! Stephan Thober Oct 2015 - ported to mRM
978 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
979 :
980 0 : subroutine mrm_write_optifile(best_OF, best_paramSet, param_names)
981 :
982 0 : use mo_common_mhm_mrm_file, only : file_opti, uopti
983 : use mo_common_variables, only : dirConfigOut
984 : use mo_string_utils, only : num2str
985 :
986 : implicit none
987 :
988 : ! best objective function value as returnedby the optimization routine
989 : real(dp), intent(in) :: best_OF
990 :
991 : ! best associated global parameter setCalled only when optimize is .TRUE.
992 : real(dp), dimension(:), intent(in) :: best_paramSet
993 :
994 : character(len = *), dimension(:), intent(in) :: param_names
995 :
996 : character(256) :: fName, formHeader, formParams
997 :
998 : integer(i4) :: ii, err, n_params
999 :
1000 :
1001 : ! number of parameters
1002 0 : n_params = size(best_paramSet)
1003 :
1004 : ! open file
1005 0 : fName = trim(adjustl(dirConfigOut)) // trim(adjustl(file_opti))
1006 0 : open(uopti, file = fName, status = 'unknown', action = 'write', iostat = err, recl = (n_params + 1) * 40)
1007 0 : if(err .ne. 0) then
1008 0 : call error_message(' IOError while openening "', trim(fName), '" Error-Code ', num2str(err))
1009 : end if
1010 :
1011 : ! header
1012 0 : write(formHeader, *) '(a40,', n_params, 'a40)'
1013 : ! len(param_names(1))=256 but only 39 characters taken here
1014 : ! write(uopti, formHeader) 'OF', (trim(adjustl(param_names(ii))), ii=1, n_params)
1015 0 : write(uopti, formHeader) 'OF', (trim(adjustl(param_names(ii)(1 : 39))), ii = 1, n_params)
1016 :
1017 : ! output
1018 0 : write(formParams, *) '( es40.14, ', n_params, '(es40.14) )'
1019 0 : write(uopti, formParams) best_OF, (best_paramSet(ii), ii = 1, n_params)
1020 :
1021 : ! close file
1022 0 : close(uopti)
1023 :
1024 : ! screen output
1025 0 : call message()
1026 0 : call message(' Optimized parameters written to ', trim(fName))
1027 :
1028 0 : end subroutine mrm_write_optifile
1029 :
1030 : ! ------------------------------------------------------------------
1031 :
1032 : ! NAME
1033 : ! mrm_write_optinamelist
1034 :
1035 : ! PURPOSE
1036 : !> \brief Write final, optimized parameter set in a namelist format.
1037 :
1038 : !> \details Write final, optimized parameter set in a namelist format.
1039 :
1040 : ! INTENT(IN)
1041 : !> \param[in] "real(dp), dimension(:, :) :: parameters" (min, max, opti)
1042 : !> \param[in] "logical, dimension(size(parameters, 1)) :: maskpara" .true. if parameter was
1043 : !> calibrated
1044 : !> \param[in] "character(len = *), dimension(size(parameters, 1)) :: parameters_name" clear names of parameters
1045 :
1046 : ! HISTORY
1047 : !> \authors Juliane Mai
1048 :
1049 : !> \date Dec 2013
1050 :
1051 : ! Modifications:
1052 : ! Stephan Thober Oct 2015 - adapted to mRM
1053 : ! Stephan Thober Nov 2016 - adapt header to routing process
1054 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
1055 :
1056 0 : subroutine mrm_write_optinamelist(parameters, maskpara, parameters_name)
1057 :
1058 0 : use mo_common_mhm_mrm_file, only : file_opti_nml, uopti_nml
1059 : use mo_common_variables, only : dirConfigOut, processMatrix
1060 : use mo_string_utils, only : num2str
1061 :
1062 : implicit none
1063 :
1064 : ! (min, max, opti)
1065 : real(dp), dimension(:, :), intent(in) :: parameters
1066 :
1067 : ! .true. if parameter was calibrated
1068 : logical, dimension(size(parameters, 1)), intent(in) :: maskpara
1069 :
1070 : ! clear names of parameters
1071 : character(len = *), dimension(size(parameters, 1)), intent(in) :: parameters_name
1072 :
1073 : character(256) :: fName
1074 :
1075 : character(3) :: flag
1076 :
1077 : integer(i4) :: err
1078 :
1079 : integer(i4) :: iPar
1080 :
1081 :
1082 : ! open file
1083 0 : fName = trim(adjustl(dirConfigOut)) // trim(adjustl(file_opti_nml))
1084 0 : open(uopti_nml, file = fName, status = 'unknown', action = 'write', iostat = err)
1085 0 : if(err .ne. 0) then
1086 0 : call message (' IOError while openening "', trim(fName), '" Error-Code ', num2str(err))
1087 : end if
1088 :
1089 0 : write(uopti_nml, *) '!global_parameters'
1090 0 : write(uopti_nml, *) '!PARAMETER lower_bound upper_bound value FLAG SCALING'
1091 :
1092 0 : write(uopti_nml, *) '! ', trim(adjustl('routing'))
1093 :
1094 0 : if (processMatrix(8, 1) .eq. 1_i4) write(uopti_nml, *) '&routing1'
1095 0 : if (ProcessMatrix(8, 1) .eq. 2_i4) write(uopti_nml, *) '&routing2'
1096 0 : if (ProcessMatrix(8, 1) .eq. 3_i4) write(uopti_nml, *) '&routing3'
1097 :
1098 0 : do iPar = 1, size(parameters, 1)
1099 0 : if (maskpara(iPar)) then
1100 0 : flag = ' 1 '
1101 : else
1102 0 : flag = ' 0 '
1103 : end if
1104 0 : write(uopti_nml, *) trim(adjustl(parameters_name(iPar))), ' = ', &
1105 0 : parameters(iPar, 1), ' , ', &
1106 0 : parameters(iPar, 2), ' , ', &
1107 0 : parameters(iPar, 3), ' , ', &
1108 0 : flag, ', 1 '
1109 : end do
1110 :
1111 0 : write(uopti_nml, *) '/'
1112 0 : write(uopti_nml, *) ' '
1113 :
1114 : ! close file
1115 0 : close(uopti_nml)
1116 :
1117 : ! screen output
1118 0 : call message()
1119 0 : call message(' Optimized parameters written in namelist format to ', trim(fName))
1120 :
1121 0 : end subroutine mrm_write_optinamelist
1122 :
1123 :
1124 : end module mo_mrm_write
|
