Line data Source code
1 : !> \file mo_common_mHM_mRM_read_config.f90
2 : !> \brief \copybrief mo_common_mhm_mrm_read_config
3 : !> \details \copydetails mo_common_mhm_mrm_read_config
4 :
5 : !> \brief Reading of main model configurations.
6 : !> \details This routine reads the configurations of common program parts
7 : !> \authors Matthias Zink
8 : !> \date Dec 2012
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_common
12 : MODULE mo_common_mHM_mRM_read_config
13 :
14 : use mo_kind, only : i4, dp
15 : use mo_message, only : message, error_message
16 :
17 : IMPLICIT NONE
18 :
19 : PRIVATE
20 :
21 : PUBLIC :: common_mHM_mRM_read_config, common_check_resolution, check_optimization_settings
22 :
23 : CONTAINS
24 :
25 :
26 : !> \brief Read main configurations for common parts
27 : !> \changelog
28 : !! - Robert Schweppe Dec 2017
29 : !! - based on mhm_read_config
30 : !! - Stephan Thober Jan 2022
31 : !! - added nTStepForcingDay
32 : !> \authors Matthias Zink
33 : !> \date Dec 2012
34 14 : subroutine common_mHM_mRM_read_config(file_namelist, unamelist)
35 :
36 : use mo_common_constants, only : maxNoDomains, nodata_i4
37 : use mo_common_mHM_mRM_variables, only : LCyearId, dds_r, mhmFileRestartIn, mrmFileRestartIn, evalPer,&
38 : mcmc_error_params, mcmc_opti, nIterations, &
39 : nTStepDay, opti_function, opti_method, optimize, optimize_restart, &
40 : read_restart, mrm_read_river_network, resolutionRouting, sa_temp, &
41 : sce_ngs, sce_npg, sce_nps, seed, &
42 : simPer, timestep, warmPer, warmingDays, read_old_style_restart_bounds, &
43 : restart_reset_fluxes_states
44 : use mo_common_read_config, only : set_land_cover_scenes_id
45 : use mo_common_types, only: period
46 : use mo_common_variables, only : LCfilename, domainMeta, processMatrix
47 : use mo_julian, only : caldat, julday
48 : use mo_nml, only : close_nml, open_nml, position_nml
49 : use mo_string_utils, only : num2str
50 :
51 : implicit none
52 :
53 : character(*), intent(in) :: file_namelist !< namelist file name
54 : integer, intent(in) :: unamelist !< unit to open namelist file
55 :
56 : integer(i4) :: jday
57 :
58 : integer(i4) :: domainID, iDomain
59 :
60 : integer(i4), dimension(maxNoDomains) :: warming_Days
61 :
62 : type(period), dimension(maxNoDomains) :: eval_Per
63 :
64 714 : real(dp), dimension(maxNoDomains) :: resolution_Routing
65 :
66 : character(256), dimension(maxNoDomains) :: mhm_file_RestartIn
67 : character(256), dimension(maxNoDomains) :: mrm_file_RestartIn
68 :
69 :
70 : ! namelist spatial & temporal resolution, otmization information
71 : namelist /mainconfig_mhm_mrm/ timestep, resolution_Routing, optimize, &
72 : optimize_restart, opti_method, opti_function, &
73 : read_restart, mrm_read_river_network, read_old_style_restart_bounds, restart_reset_fluxes_states, &
74 : mhm_file_RestartIn, mrm_file_RestartIn
75 : ! namelist for optimization settings
76 : namelist /Optimization/ nIterations, seed, dds_r, sa_temp, sce_ngs, &
77 : sce_npg, sce_nps, mcmc_opti, mcmc_error_params
78 : ! namelist for time settings
79 : namelist /time_periods/ warming_Days, eval_Per
80 :
81 : ! set default values for optional arguments
82 14 : mrm_read_river_network = .false.
83 14 : read_old_style_restart_bounds = .false.
84 14 : restart_reset_fluxes_states = .false.
85 :
86 : !===============================================================
87 : ! Read namelist main directories
88 : !===============================================================
89 14 : call open_nml(file_namelist, unamelist, quiet = .true.)
90 :
91 : !===============================================================
92 : ! Read namelist specifying the model configuration
93 : !===============================================================
94 14 : call position_nml('mainconfig_mhm_mrm', unamelist)
95 14 : read(unamelist, nml = mainconfig_mhm_mrm)
96 : ! consistency between read_restart and mrm_read_river_network
97 14 : if (read_restart) then
98 1 : if (.not. mrm_read_river_network) then
99 1 : call message('***WARNING: mrm_read_river_network is set to .true. because read_restart is .true.')
100 : end if
101 1 : mrm_read_river_network = .true.
102 : end if
103 :
104 42 : allocate(resolutionRouting(domainMeta%nDomains))
105 42 : allocate(mhmFileRestartIn(domainMeta%nDomains))
106 28 : allocate(mrmFileRestartIn(domainMeta%nDomains))
107 40 : do iDomain = 1, domainMeta%nDomains
108 26 : domainID = domainMeta%indices(iDomain)
109 26 : mhmFileRestartIn(iDomain) = mhm_file_RestartIn(domainID)
110 26 : mrmFileRestartIn(iDomain) = mrm_file_RestartIn(domainID)
111 40 : resolutionRouting(iDomain) = resolution_Routing(domainID)
112 : end do
113 :
114 : ! check for optimize and read restart
115 14 : if ((read_restart) .and. (optimize)) then
116 0 : call message()
117 0 : call error_message('***ERROR: cannot read states from restart file when optimizing')
118 : end if
119 :
120 40 : do iDomain = 1, domainMeta%nDomains
121 26 : domainID = domainMeta%indices(iDomain)
122 40 : if (processMatrix(8, 1) > 0 .and. domainMeta%optidata(iDomain) > 1 .and. optimize) then
123 5 : domainMeta%doRouting(iDomain) = .FALSE.
124 5 : call message('Warning: although defined in namelist, routing is switched off for domain', trim(num2str(domainID)))
125 5 : call message(' since the calibration of Q is not possible with the chosen opti input')
126 : end if
127 : end do
128 :
129 : !===============================================================
130 : ! INIT !!! (merged from mo_startup and mo_mrm_read_config)
131 : !===============================================================
132 : ! transformation of time units & constants
133 14 : if (mod(24, timeStep) > 0) then
134 0 : call error_message('mo_startup: timeStep must be a divisor of 24: ', num2str(timeStep))
135 : end if
136 14 : nTStepDay = 24_i4 / timeStep ! # of time steps per day
137 :
138 : ! allocate time periods
139 42 : allocate(simPer(domainMeta%nDomains))
140 28 : allocate(evalPer(domainMeta%nDomains))
141 42 : allocate(warmingDays(domainMeta%nDomains))
142 28 : allocate(warmPer(domainMeta%nDomains))
143 :
144 : !===============================================================
145 : ! read simulation time periods incl. warming days
146 : !===============================================================
147 14 : call position_nml('time_periods', unamelist)
148 14 : read(unamelist, nml = time_periods)
149 40 : do iDomain = 1, domainMeta%nDomains
150 26 : domainID = domainMeta%indices(iDomain)
151 26 : warmingDays(iDomain) = warming_Days(domainID)
152 : ! this will be a procedure subroutine
153 : ! therefore inout first, in second
154 40 : call period_copy_period_data(evalPer(iDomain), eval_Per(domainID))
155 : end do
156 : ! evalPer = eval_Per(1 : domainMeta%nDomains)
157 :
158 : !===============================================================
159 : ! determine simulation time period incl. warming days for each
160 : ! domain
161 : !===============================================================
162 40 : do iDomain = 1, domainMeta%nDomains
163 : ! julian days for evaluation period
164 26 : jday = julday(dd = evalPer(iDomain)%dStart, mm = evalPer(iDomain)%mStart, yy = evalPer(iDomain)%yStart)
165 26 : evalPer(iDomain)%julStart = jday
166 :
167 26 : jday = julday(dd = evalPer(iDomain)%dEnd, mm = evalPer(iDomain)%mEnd, yy = evalPer(iDomain)%yEnd)
168 26 : evalPer(iDomain)%julEnd = jday
169 :
170 : ! determine warming period
171 26 : warmPer(iDomain)%julStart = evalPer(iDomain)%julStart - warmingDays(iDomain)
172 26 : warmPer(iDomain)%julEnd = evalPer(iDomain)%julStart - 1
173 :
174 0 : call caldat(warmPer(iDomain)%julStart, dd = warmPer(iDomain)%dStart, mm = warmPer(iDomain)%mStart, &
175 26 : yy = warmPer(iDomain)%yStart)
176 0 : call caldat(warmPer(iDomain)%julEnd, dd = warmPer(iDomain)%dEnd, mm = warmPer(iDomain)%mEnd, &
177 26 : yy = warmPer(iDomain)%yEnd)
178 :
179 : ! simulation Period = warming Period + evaluation Period
180 26 : simPer(iDomain)%dStart = warmPer(iDomain)%dStart
181 26 : simPer(iDomain)%mStart = warmPer(iDomain)%mStart
182 26 : simPer(iDomain)%yStart = warmPer(iDomain)%yStart
183 26 : simPer(iDomain)%julStart = warmPer(iDomain)%julStart
184 26 : simPer(iDomain)%dEnd = evalPer(iDomain)%dEnd
185 26 : simPer(iDomain)%mEnd = evalPer(iDomain)%mEnd
186 26 : simPer(iDomain)%yEnd = evalPer(iDomain)%yEnd
187 92 : simPer(iDomain)%julEnd = evalPer(iDomain)%julEnd
188 : end do
189 :
190 : call set_land_cover_scenes_id(simPer, LCyearId)
191 :
192 : !===============================================================
193 : ! Settings for Optimization
194 : !===============================================================
195 : ! namelist for Optimization settings
196 14 : call position_nml('Optimization', unamelist)
197 14 : read(unamelist, nml = Optimization)
198 : ! checking of settings and default value initialization moved to new subroutine
199 : ! because global_parameters need to be set, which is not the case right now
200 14 : call close_nml(unamelist)
201 :
202 14 : end subroutine common_mHM_mRM_read_config
203 :
204 :
205 : !> \brief check optimization settings
206 : !> \authors Robert Schweppe
207 : !> \date Jun 2018
208 14 : subroutine check_optimization_settings
209 :
210 14 : use mo_common_mHM_mRM_variables, only : dds_r, nIterations, sce_ngs, sce_npg, sce_nps
211 : use mo_common_variables, only : global_parameters
212 :
213 : implicit none
214 :
215 : integer(i4) :: n_true_pars
216 :
217 :
218 : ! check and set default values
219 14 : if (nIterations .le. 0_i4) then
220 0 : call error_message('Number of iterations for Optimization (nIterations) must be greater than zero')
221 : end if
222 14 : if (dds_r .lt. 0.0_dp .or. dds_r .gt. 1.0_dp) then
223 0 : call error_message('dds_r must be between 0.0 and 1.0')
224 : end if
225 14 : if (sce_ngs .lt. 1_i4) then
226 0 : call error_message('number of complexes in SCE (sce_ngs) must be at least 1')
227 : end if
228 : ! number of points in each complex: default = 2n+1
229 14 : if (sce_npg .lt. 0_i4) then
230 757 : n_true_pars = count(nint(global_parameters(:, 4)) .eq. 1)
231 14 : sce_npg = 2 * n_true_pars + 1_i4
232 : end if
233 : ! number of points in each sub-complex: default = n+1
234 14 : if (sce_nps .lt. 0_i4) then
235 757 : n_true_pars = count(nint(global_parameters(:, 4)) .eq. 1)
236 14 : sce_nps = n_true_pars + 1_i4
237 : end if
238 14 : if (sce_npg .lt. sce_nps) then
239 0 : call error_message('number of points per complex (sce_npg) must be greater or', raise=.false.)
240 0 : call error_message('equal number of points per sub-complex (sce_nps)')
241 : end if
242 :
243 14 : end subroutine check_optimization_settings
244 :
245 :
246 : !> \brief check resolution
247 : !> \authors Robert Schweppe
248 : !> \date Jun 2018
249 28 : subroutine common_check_resolution(do_message, allow_subgrid_routing)
250 :
251 14 : use mo_common_mHM_mRM_variables, only : resolutionRouting
252 : use mo_common_variables, only : domainMeta, resolutionHydrology
253 : use mo_string_utils, only : num2str
254 :
255 : implicit none
256 :
257 : logical, intent(in) :: do_message !< flag to print messages
258 : logical, intent(in) :: allow_subgrid_routing !< flag to allow subgrid routing
259 :
260 : integer(i4) :: iDomain, domainID
261 :
262 : ! conversion factor L11 to L1
263 28 : real(dp) :: cellFactorRbyH
264 :
265 :
266 : !===============================================================
267 : ! check matching of resolutions: hydrology, forcing and routing
268 : !===============================================================
269 80 : do iDomain = 1, domainMeta%nDomains
270 52 : domainID = domainMeta%indices(iDomain)
271 52 : cellFactorRbyH = resolutionRouting(iDomain) / resolutionHydrology(iDomain)
272 52 : if (do_message) then
273 26 : call message()
274 26 : call message('domain ', trim(adjustl(num2str(domainID))), ': ')
275 : call message('resolution Hydrology (domain ', trim(adjustl(num2str(domainID))), ') = ', &
276 26 : trim(adjustl(num2str(resolutionHydrology(iDomain)))))
277 : call message('resolution Routing (domain ', trim(adjustl(num2str(domainID))), ') = ', &
278 26 : trim(adjustl(num2str(resolutionRouting(iDomain)))))
279 : end if
280 : !
281 80 : if(nint(cellFactorRbyH * 100.0_dp) .eq. 100) then
282 50 : if (do_message) then
283 25 : call message()
284 25 : call message('Resolution of routing and hydrological modeling are equal!')
285 : end if
286 :
287 2 : else if ((nint(cellFactorRbyH * 100.0_dp) .gt. 100) .and. .not.allow_subgrid_routing) then
288 1 : if(nint(mod(cellFactorRbyH, 2.0_dp) * 100.0_dp) .ne. 0) then
289 0 : call error_message('***ERROR: Resolution of routing is not a multiple of hydrological model resolution!', raise=.false.)
290 0 : call error_message(' FILE: mhm.nml, namelist: mainconfig, variable: resolutionRouting')
291 : end if
292 : !
293 1 : if (do_message) then
294 1 : call message()
295 : call message('Resolution of routing is bigger than hydrological model resolution by ', &
296 1 : trim(adjustl(num2str(nint(cellFactorRbyH)))), ' times !')
297 : end if
298 : end if
299 : !
300 : end do
301 :
302 28 : end subroutine common_check_resolution
303 :
304 :
305 : ! ToDo: make this a procedure of period
306 : !> \brief copy period data
307 26 : subroutine period_copy_period_data(toPeriod, fromPeriod)
308 28 : use mo_common_types, only: period
309 : type(period), intent(inout) :: toPeriod !< copy to this period
310 : type(period), intent(in) :: fromPeriod !< copy from this period
311 :
312 26 : toPeriod%dStart = fromPeriod%dStart ! first day
313 26 : toPeriod%mStart = fromPeriod%mStart ! first month
314 26 : toPeriod%yStart = fromPeriod%yStart ! first year
315 26 : toPeriod%dEnd = fromPeriod%dEnd ! last day
316 26 : toPeriod%mEnd = fromPeriod%mEnd ! last month
317 26 : toPeriod%yEnd = fromPeriod%yEnd ! last year
318 26 : toPeriod%julStart = 0 ! first julian day
319 26 : toPeriod%julEnd = 0 ! last julian day
320 26 : toPeriod%nObs = 0 ! total number of observations
321 :
322 26 : end subroutine period_copy_period_data
323 :
324 : END MODULE mo_common_mHM_mRM_read_config
|
