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