Line data Source code
1 : !> \dir meteo
2 : !> \brief \copybrief f_meteo
3 : !> \details \copydetails f_meteo
4 :
5 : !> \defgroup f_meteo Meteo - Fortran modules
6 : !> \brief Core modules to deal with meteorological forcings.
7 : !> \details These modules provide the meteo handler, the spatial and temporal remapping algorithms and helper routines.
8 :
9 : !> \file mo_meteo_handler.f90
10 : !> \brief \copybrief mo_meteo_handler
11 : !> \details \copydetails mo_meteo_handler
12 :
13 : !> \brief Class for the meteo handler
14 : !> \details Handler for meteorological forcings in mHM.
15 : !! Is independent of global variables and provides 3 methods to access forcings:
16 : !! - get_corrected_pet : get the modified pet for mHM for the current timestep
17 : !! - get_temp : get the temporal disaggregated temperature for the current timestep
18 : !! - get_prec : get the temporal disaggregated percipitation for the current timestep
19 : !!
20 : !> \version 0.1
21 : !> \authors Sebastian Mueller
22 : !> \date Mar 2023
23 : !> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
24 : !! mHM is released under the LGPLv3+ license \license_note
25 : !> \ingroup f_meteo
26 : module mo_meteo_handler
27 :
28 : USE mo_kind, ONLY : i4, dp
29 : USE mo_constants, ONLY : YearMonths
30 : use mo_common_types, only: Grid, period
31 : use mo_message, only : message, error_message
32 : use mo_coupling_type, only : couple_cfg_type
33 : use mo_sentinel, only : set_sentinel, check_sentinel
34 : use mo_datetime, only : datetime, timedelta, zero_delta, one_hour, one_day
35 :
36 : implicit none
37 :
38 : private
39 :
40 : !> \class meteo_handler_type
41 : !> \brief This is a handler for the meteorological forcings
42 : !> \details This class provides all procedures to handle meteorological forcings from file or interfaces.
43 : type, public :: meteo_handler_type
44 : ! config settings
45 : character(256) :: dir_nml_name = 'directories_mHM' !< namelist name in mhm namelist
46 : character(256) :: weight_nml_name = 'nightDayRatio' !< namelist name in mhm namelist
47 : !> Input nCols and nRows of binary meteo and LAI files are in header file
48 : CHARACTER(256) :: file_meteo_header = 'header.txt'
49 : !> Unit for meteo header file
50 : INTEGER :: umeteo_header = 50
51 : !> .FALSE. to only warn about bound (lower, upper) violations in meteo files, default = .TRUE. - raise an error
52 : logical :: bound_error
53 : integer(i4) :: pet_case !< process case for PET (processCase(5))
54 : integer(i4) :: riv_temp_case !< process case for river temperature (processCase(11))
55 : type(period), public :: readPer !< start and end dates of read period
56 : integer(i4), public :: nTstepForcingDay !< Number of forcing intervals per day
57 : !> flag wether forcings are given at hourly timestep
58 : logical, public :: is_hourly_forcing
59 : integer(i4), public :: timeStep !< [h] simulation time step (= TS) in [h]
60 : integer(i4), public :: nTstepDay !< Number of time intervals per day
61 : real(dp), public :: nTstepDay_dp !< Number of time intervals per day as real
62 : ! -------------------------------------------------------------------
63 : ! level 2 description
64 : ! -------------------------------------------------------------------
65 : type(Grid), dimension(:), allocatable, public :: level2 !< Reference of the metereological variables
66 : ! -------------------------------------------------------------------
67 : ! INPUT variables for configuration of mHM
68 : ! -------------------------------------------------------------------
69 : integer(i4), dimension(:), allocatable, public :: timeStep_model_inputs !< frequency for reading meteo input
70 : logical, public :: read_meteo_weights !< read weights for tavg and pet
71 : character(256), public :: inputFormat_meteo_forcings !< format of meteo input data (nc)
72 : ! ------------------------------------------------------------------
73 : ! DIRECTORIES
74 : ! ------------------------------------------------------------------
75 : ! has the dimension of nDomains
76 : character(256), dimension(:), allocatable, public :: dir_meteo_header !< Directory where the meteo header file is located
77 : character(256), dimension(:), allocatable, public :: dirPrecipitation !< Directory where precipitation files are located
78 : character(256), dimension(:), allocatable, public :: dirTemperature !< Directory where temperature files are located
79 : character(256), dimension(:), allocatable, public :: dirMinTemperature !< Directory where minimum temp. files are located
80 : character(256), dimension(:), allocatable, public :: dirMaxTemperature !< Directory where maximum temp. files are located
81 : character(256), dimension(:), allocatable, public :: dirNetRadiation !< Directory where abs. vap. pressure files are located
82 : character(256), dimension(:), allocatable, public :: dirabsVapPressure !< Directory where abs. vap. pressure files are located
83 : character(256), dimension(:), allocatable, public :: dirwindspeed !< Directory where windspeed files are located
84 : character(256), dimension(:), allocatable, public :: dirReferenceET !< Directory where reference-ET files are located
85 : ! riv-temp releated
86 : character(256), dimension(:), allocatable, public :: dirRadiation !< Directory where short/long-wave rad. files are located
87 : ! -------------------------------------------------------------------
88 : ! L1 DOMAIN description
89 : ! -------------------------------------------------------------------
90 : ! Forcings
91 : ! dim1 = number grid cells L1
92 : ! dim2 = number of meteorological time steps
93 : ! dim2 = month of year
94 : real(dp), public, dimension(:, :, :), allocatable :: L1_temp_weights !< hourly temperature weights for daily values
95 : real(dp), public, dimension(:, :, :), allocatable :: L1_pet_weights !< hourly pet weights for daily values
96 : real(dp), public, dimension(:, :, :), allocatable :: L1_pre_weights !< hourly pre weights for daily values
97 : real(dp), public, dimension(:, :), allocatable :: L1_pre !< [mm] Precipitation
98 : real(dp), public, dimension(:, :), allocatable :: L1_temp !< [degC] Air temperature
99 : real(dp), public, dimension(:, :), allocatable :: L1_pet !< [mm TS-1] Potential evapotranspiration
100 : real(dp), public, dimension(:, :), allocatable :: L1_tmin !< [degC] minimum daily air temperature
101 : real(dp), public, dimension(:, :), allocatable :: L1_tmax !< [degC] maximum daily air temperature
102 : real(dp), public, dimension(:, :), allocatable :: L1_netrad !< [W m2] net radiation
103 : real(dp), public, dimension(:, :), allocatable :: L1_absvappress !< [Pa] absolute vapour pressure
104 : real(dp), public, dimension(:, :), allocatable :: L1_windspeed !< [m s-1] windspeed
105 : ! riv-temp related
106 : real(dp), public, dimension(:, :), allocatable :: L1_ssrd !< [W m2] short wave radiation
107 : real(dp), public, dimension(:, :), allocatable :: L1_strd !< [W m2] long wave radiation
108 : real(dp), public, dimension(:, :), allocatable :: L1_tann !< [degC] annual mean air temperature
109 : ! -------------------------------------------------------------------
110 : ! Monthly day/night variation of Meteorological variables
111 : ! for temporal disaggregation
112 : ! -------------------------------------------------------------------
113 : ! dim1 = number of months in a year
114 : real(dp), public, dimension(int(YearMonths, i4)) :: fday_prec !< [-] Day ratio precipitation < 1
115 : real(dp), public, dimension(int(YearMonths, i4)) :: fnight_prec !< [-] Night ratio precipitation < 1
116 : real(dp), public, dimension(int(YearMonths, i4)) :: fday_pet !< [-] Day ratio PET < 1
117 : real(dp), public, dimension(int(YearMonths, i4)) :: fnight_pet !< [-] Night ratio PET < 1
118 : real(dp), public, dimension(int(YearMonths, i4)) :: fday_temp !< [-] Day factor mean temp
119 : real(dp), public, dimension(int(YearMonths, i4)) :: fnight_temp !< [-] Night factor mean temp
120 : real(dp), public, dimension(int(YearMonths, i4)) :: fday_ssrd !< [-] Day factor short-wave rad.
121 : real(dp), public, dimension(int(YearMonths, i4)) :: fnight_ssrd !< [-] Night factor short-wave rad.
122 : real(dp), public, dimension(int(YearMonths, i4)) :: fday_strd !< [-] Day factor long-wave rad.
123 : real(dp), public, dimension(int(YearMonths, i4)) :: fnight_strd !< [-] Night factor long-wave rad.
124 :
125 : ! -------------------------------------------------------------------
126 : ! current mHM array indizes
127 : ! -------------------------------------------------------------------
128 : !> start index of meteo variables
129 : integer(i4) :: s_meteo
130 : !> end index of meteo variables
131 : integer(i4) :: e_meteo
132 : !> index of meteo time-step
133 : integer(i4) :: iMeteoTS
134 : !> start index of level-1 variables
135 : integer(i4) :: s1
136 : !> end index of level-1 variables
137 : integer(i4) :: e1
138 : !> number of domains
139 : integer(i4) :: nDomains
140 : integer(i4), dimension(:), allocatable :: indices !< indices
141 : integer(i4), dimension(:), allocatable :: L0DataFrom !< index of associated level-0 domain
142 : !> current julian time
143 : real(dp) :: time
144 :
145 : ! -------------------------------------------------------------------
146 : ! coupling settings
147 : ! -------------------------------------------------------------------
148 : type(couple_cfg_type), public :: couple_cfg !< coupling configuration class
149 : type(timedelta) :: couple_step_delta !< timedelta for the coupling meteo time-step
150 : type(datetime) :: couple_pre_time !< current time from coupling for pre
151 : type(datetime) :: couple_temp_time !< current time from coupling for temp
152 : type(datetime) :: couple_pet_time !< current time from coupling for pet
153 : type(datetime) :: couple_tmin_time !< current time from coupling for tmin
154 : type(datetime) :: couple_tmax_time !< current time from coupling for tmax
155 : type(datetime) :: couple_netrad_time !< current time from coupling for netrad
156 : type(datetime) :: couple_absvappress_time !< current time from coupling for absvappress
157 : type(datetime) :: couple_windspeed_time !< current time from coupling for windspeed
158 : type(datetime) :: couple_ssrd_time !< current time from coupling for ssrd
159 : type(datetime) :: couple_strd_time !< current time from coupling for strd
160 : type(datetime) :: couple_tann_time !< current time from coupling for tann
161 : logical, public :: couple_pre !< coupling config for pre
162 : logical, public :: couple_temp !< coupling config for temp
163 : logical, public :: couple_pet !< coupling config for pet
164 : logical, public :: couple_tmin !< coupling config for tmin
165 : logical, public :: couple_tmax !< coupling config for tmax
166 : logical, public :: couple_netrad !< coupling config for netrad
167 : logical, public :: couple_absvappress !< coupling config for absvappress
168 : logical, public :: couple_windspeed !< coupling config for windspeed
169 : logical, public :: couple_ssrd !< coupling config for ssrd
170 : logical, public :: couple_strd !< coupling config for strd
171 : logical, public :: couple_tann !< coupling config for tann
172 : logical, public :: couple_all !< flag to indicated that all meteo-data is coming from the coupler
173 : logical, public :: couple_is_hourly !< flag to indicated hourly data from coupler
174 : contains
175 : !> \copydoc mo_meteo_handler::clean_up
176 : procedure :: clean_up !< \see mo_meteo_handler::clean_up
177 : !> \copydoc mo_meteo_handler::config
178 : procedure :: config !< \see mo_meteo_handler::config
179 : !> \copydoc mo_meteo_handler::single_read
180 : procedure :: single_read !< \see mo_meteo_handler::single_read
181 : !> \copydoc mo_meteo_handler::init_level2
182 : procedure :: init_level2 !< \see mo_meteo_handler::init_level2
183 : !> \copydoc mo_meteo_handler::prepare_data
184 : procedure :: prepare_data !< \see mo_meteo_handler::prepare_data
185 : !> \copydoc mo_meteo_handler::update_timestep
186 : procedure :: update_timestep !< \see mo_meteo_handler::update_timestep
187 : !> \copydoc mo_meteo_handler::get_corrected_pet
188 : procedure :: get_corrected_pet !< \see mo_meteo_handler::get_corrected_pet
189 : !> \copydoc mo_meteo_handler::get_temp
190 : procedure :: get_temp !< \see mo_meteo_handler::get_temp
191 : !> \copydoc mo_meteo_handler::get_prec
192 : procedure :: get_prec !< \see mo_meteo_handler::get_prec
193 : !> \copydoc mo_meteo_handler::get_ssrd
194 : procedure :: get_ssrd !< \see mo_meteo_handler::get_ssrd
195 : !> \copydoc mo_meteo_handler::get_strd
196 : procedure :: get_strd !< \see mo_meteo_handler::get_strd
197 : !> \copydoc mo_meteo_handler::get_tann
198 : procedure :: get_tann !< \see mo_meteo_handler::get_tann
199 : !> \copydoc mo_meteo_handler::set_meteo
200 : procedure :: set_meteo !< \see mo_meteo_handler::set_meteo
201 : end type meteo_handler_type
202 :
203 : contains
204 :
205 : !> \brief clean up
206 14 : subroutine clean_up(self)
207 : implicit none
208 :
209 : class(meteo_handler_type), intent(inout) :: self
210 :
211 14 : if ( allocated(self%indices) ) deallocate(self%indices)
212 14 : if ( allocated(self%L0DataFrom) ) deallocate(self%L0DataFrom)
213 14 : if ( allocated(self%timeStep_model_inputs) ) deallocate(self%timeStep_model_inputs)
214 14 : if ( allocated(self%dir_meteo_header) ) deallocate(self%dir_meteo_header)
215 14 : if ( allocated(self%dirPrecipitation) ) deallocate(self%dirPrecipitation)
216 14 : if ( allocated(self%dirTemperature) ) deallocate(self%dirTemperature)
217 14 : if ( allocated(self%dirMinTemperature) ) deallocate(self%dirMinTemperature)
218 14 : if ( allocated(self%dirMaxTemperature) ) deallocate(self%dirMaxTemperature)
219 14 : if ( allocated(self%dirNetRadiation) ) deallocate(self%dirNetRadiation)
220 14 : if ( allocated(self%dirabsVapPressure) ) deallocate(self%dirabsVapPressure)
221 14 : if ( allocated(self%dirwindspeed) ) deallocate(self%dirwindspeed)
222 14 : if ( allocated(self%dirReferenceET) ) deallocate(self%dirReferenceET)
223 14 : if ( allocated(self%dirRadiation) ) deallocate(self%dirRadiation)
224 40 : if ( allocated(self%level2) ) deallocate(self%level2)
225 14 : if ( allocated(self%L1_temp_weights) ) deallocate(self%L1_temp_weights)
226 14 : if ( allocated(self%L1_pet_weights) ) deallocate(self%L1_pet_weights)
227 14 : if ( allocated(self%L1_pre_weights) ) deallocate(self%L1_pre_weights)
228 14 : if ( allocated(self%L1_pre) ) deallocate(self%L1_pre)
229 14 : if ( allocated(self%L1_temp) ) deallocate(self%L1_temp)
230 14 : if ( allocated(self%L1_pet) ) deallocate(self%L1_pet)
231 14 : if ( allocated(self%L1_tmin) ) deallocate(self%L1_tmin)
232 14 : if ( allocated(self%L1_tmax) ) deallocate(self%L1_tmax)
233 14 : if ( allocated(self%L1_netrad) ) deallocate(self%L1_netrad)
234 14 : if ( allocated(self%L1_absvappress) ) deallocate(self%L1_absvappress)
235 14 : if ( allocated(self%L1_windspeed) ) deallocate(self%L1_windspeed)
236 14 : if ( allocated(self%L1_ssrd) ) deallocate(self%L1_ssrd)
237 14 : if ( allocated(self%L1_strd) ) deallocate(self%L1_strd)
238 14 : if ( allocated(self%L1_tann) ) deallocate(self%L1_tann)
239 :
240 14 : end subroutine clean_up
241 :
242 : !> \brief configure the \ref meteo_handler_type class from the mhm namelist
243 14 : subroutine config(self, file_namelist, unamelist, optimize, domainMeta, processMatrix, timeStep, couple_cfg)
244 :
245 14 : use mo_common_constants, only : maxNoDomains, nodata_i4
246 : use mo_common_types, only : domain_meta
247 : use mo_nml, only : close_nml, open_nml, position_nml
248 : use mo_common_variables, only : nProcesses
249 :
250 : implicit none
251 :
252 : class(meteo_handler_type), intent(inout) :: self
253 : character(*), intent(in) :: file_namelist !< mhm namelist file
254 : integer, intent(in) :: unamelist !< unit to open namelist file
255 : logical, intent(in) :: optimize !< Optimization flag
256 : type(domain_meta), intent(in) :: domainMeta !< domain general description
257 : integer(i4), dimension(nProcesses, 3), intent(in) :: processMatrix !< Info about which process runs in which option
258 : integer(i4), intent(in) :: timeStep !< [h] simulation time step (= TS) in [h]
259 : type(couple_cfg_type), intent(in) :: couple_cfg !< coupling configuration class
260 :
261 : integer(i4), dimension(maxNoDomains) :: time_step_model_inputs
262 : character(256), dimension(maxNoDomains) :: dir_meteo_header
263 : character(256), dimension(maxNoDomains) :: dir_Precipitation
264 : character(256), dimension(maxNoDomains) :: dir_Temperature
265 : character(256), dimension(maxNoDomains) :: dir_MinTemperature
266 : character(256), dimension(maxNoDomains) :: dir_MaxTemperature
267 : character(256), dimension(maxNoDomains) :: dir_NetRadiation
268 : character(256), dimension(maxNoDomains) :: dir_windspeed
269 : character(256), dimension(maxNoDomains) :: dir_absVapPressure
270 : character(256), dimension(maxNoDomains) :: dir_ReferenceET
271 : character(256), dimension(maxNoDomains) :: dir_Radiation ! riv-temp related
272 : character(256) :: inputFormat_meteo_forcings
273 :
274 : logical :: read_meteo_weights, bound_error
275 182 : real(dp), dimension(int(YearMonths, i4)) :: fnight_prec
276 182 : real(dp), dimension(int(YearMonths, i4)) :: fnight_pet
277 182 : real(dp), dimension(int(YearMonths, i4)) :: fnight_temp
278 182 : real(dp), dimension(int(YearMonths, i4)) :: fnight_ssrd
279 182 : real(dp), dimension(int(YearMonths, i4)) :: fnight_strd
280 :
281 : integer(i4) :: domainID, iDomain
282 :
283 : ! namelist directories
284 : namelist /directories_mHM/ &
285 : inputFormat_meteo_forcings, &
286 : bound_error, &
287 : dir_meteo_header, &
288 : dir_Precipitation, &
289 : dir_Temperature, &
290 : dir_ReferenceET, &
291 : dir_MinTemperature, &
292 : dir_MaxTemperature, &
293 : dir_absVapPressure, &
294 : dir_windspeed, &
295 : dir_NetRadiation, &
296 : dir_Radiation, &
297 : time_step_model_inputs
298 :
299 : ! namelist for night-day ratio of precipitation, referenceET and temperature
300 : namelist /nightDayRatio/ &
301 : read_meteo_weights, &
302 : fnight_prec, &
303 : fnight_pet, &
304 : fnight_temp, &
305 : fnight_ssrd, &
306 : fnight_strd
307 :
308 : ! store coupling config
309 14 : self%couple_cfg = couple_cfg
310 :
311 : ! store needed domain meta infos
312 14 : self%nDomains = domainMeta%nDomains
313 42 : allocate(self%indices(self%nDomains))
314 42 : allocate(self%L0DataFrom(self%nDomains))
315 40 : self%indices(:) = domainMeta%indices(:)
316 40 : self%L0DataFrom(:) = domainMeta%L0DataFrom(:)
317 :
318 : ! # init of number of forcing timesteps, will be set when reading forcings
319 14 : self%nTStepForcingDay = nodata_i4
320 :
321 : ! store important process cases
322 14 : self%pet_case = processMatrix(5,1)
323 14 : self%riv_temp_case = processMatrix(11,1)
324 : ! store time-stepping info
325 14 : self%timeStep = timeStep
326 14 : self%nTStepDay = 24_i4 / timeStep ! # of time steps per day
327 14 : self%nTstepDay_dp = real(self%nTStepDay, dp)
328 :
329 : ! allocate variables
330 42 : allocate(self%dir_meteo_header(self%nDomains))
331 42 : allocate(self%dirPrecipitation(self%nDomains))
332 42 : allocate(self%dirTemperature(self%nDomains))
333 42 : allocate(self%dirwindspeed(self%nDomains))
334 42 : allocate(self%dirabsVapPressure(self%nDomains))
335 42 : allocate(self%dirReferenceET(self%nDomains))
336 42 : allocate(self%dirMinTemperature(self%nDomains))
337 42 : allocate(self%dirMaxTemperature(self%nDomains))
338 42 : allocate(self%dirNetRadiation(self%nDomains))
339 42 : allocate(self%dirRadiation(self%nDomains))
340 : ! allocate time periods
341 42 : allocate(self%timestep_model_inputs(self%nDomains))
342 :
343 : ! open the namelist file
344 14 : call open_nml(file_namelist, unamelist, quiet=.true.)
345 :
346 : !===============================================================
347 : ! Read namelist main directories
348 : !===============================================================
349 714 : call set_sentinel(dir_meteo_header) ! set sentinal to check reading
350 14 : inputFormat_meteo_forcings = "nc"
351 14 : bound_error = .TRUE.
352 14 : call position_nml(self%dir_nml_name, unamelist)
353 14 : read(unamelist, nml = directories_mHM)
354 :
355 14 : self%bound_error = bound_error
356 14 : self%inputFormat_meteo_forcings = inputFormat_meteo_forcings
357 :
358 40 : do iDomain = 1, self%nDomains
359 26 : domainID = self%indices(iDomain)
360 26 : self%timestep_model_inputs(iDomain) = time_step_model_inputs(domainID)
361 26 : self%dirPrecipitation(iDomain) = dir_Precipitation(domainID)
362 26 : self%dirTemperature(iDomain) = dir_Temperature(domainID)
363 26 : self%dirReferenceET(iDomain) = dir_ReferenceET(domainID)
364 26 : self%dirMinTemperature(iDomain) = dir_MinTemperature(domainID)
365 26 : self%dirMaxTemperature(iDomain) = dir_MaxTemperature(domainID)
366 26 : self%dirNetRadiation(iDomain) = dir_NetRadiation(domainID)
367 26 : self%dirwindspeed(iDomain) = dir_windspeed(domainID)
368 26 : self%dirabsVapPressure(iDomain) = dir_absVapPressure(domainID)
369 : ! riv-temp related
370 26 : self%dirRadiation(iDomain) = dir_Radiation(domainID)
371 : ! meteo header directory (if not given, use precipitation dir)
372 40 : if (check_sentinel(dir_meteo_header(domainID))) then
373 26 : self%dir_meteo_header(iDomain) = self%dirPrecipitation(iDomain)
374 : else
375 0 : self%dir_meteo_header(iDomain) = dir_meteo_header(domainID)
376 : end if
377 : end do
378 :
379 : ! consistency check for timestep_model_inputs
380 54 : if (any(self%timestep_model_inputs .ne. 0) .and. .not. all(self%timestep_model_inputs .ne. 0)) then
381 0 : call error_message('***ERROR: timestep_model_inputs either have to be all zero or all non-zero')
382 : end if
383 : ! check for optimzation and timestep_model_inputs options
384 35 : if (optimize .and. (any(self%timestep_model_inputs .ne. 0))) then
385 0 : call error_message('***ERROR: optimize and chunk read is switched on! (set timestep_model_inputs to zero)')
386 : end if
387 :
388 : !===============================================================
389 : ! Read night-day ratios and pan evaporation
390 : !===============================================================
391 : ! default values for long/shortwave rad.
392 14 : fnight_ssrd = 0.0_dp
393 182 : fnight_strd = 0.45_dp
394 :
395 : ! namelist for night-day ratio of precipitation, referenceET and temperature
396 14 : call position_nml(self%weight_nml_name, unamelist)
397 14 : read(unamelist, nml = nightDayRatio)
398 14 : self%read_meteo_weights = read_meteo_weights
399 182 : self%fnight_prec = fnight_prec
400 182 : self%fnight_pet = fnight_pet
401 182 : self%fnight_temp = fnight_temp
402 182 : self%fnight_ssrd = fnight_ssrd
403 182 : self%fnight_strd = fnight_strd
404 182 : self%fday_prec = 1.0_dp - self%fnight_prec
405 182 : self%fday_pet = 1.0_dp - self%fnight_pet
406 182 : self%fday_temp = -1.0_dp * self%fnight_temp
407 182 : self%fday_ssrd = 1.0_dp - self%fnight_ssrd
408 182 : self%fday_strd = 1.0_dp - self%fnight_strd
409 :
410 : ! TODO-RIV-TEMP:
411 : ! - add short- and long-wave raidiation weights (nc files)
412 :
413 : ! closing the namelist file
414 14 : call close_nml(unamelist)
415 :
416 : ! check coupling configuration matching process cases
417 14 : self%couple_is_hourly = .false.
418 14 : self%couple_all = .false.
419 14 : self%couple_pre = self%couple_cfg%active() .and. self%couple_cfg%meteo_expect_pre
420 14 : self%couple_temp = self%couple_cfg%active() .and. self%couple_cfg%meteo_expect_temp
421 14 : self%couple_pet = self%couple_cfg%active() .and. self%couple_cfg%meteo_expect_pet
422 14 : self%couple_tmin = self%couple_cfg%active() .and. self%couple_cfg%meteo_expect_tmin
423 14 : self%couple_tmax = self%couple_cfg%active() .and. self%couple_cfg%meteo_expect_tmax
424 14 : self%couple_netrad = self%couple_cfg%active() .and. self%couple_cfg%meteo_expect_netrad
425 14 : self%couple_absvappress = self%couple_cfg%active() .and. self%couple_cfg%meteo_expect_absvappress
426 14 : self%couple_windspeed = self%couple_cfg%active() .and. self%couple_cfg%meteo_expect_windspeed
427 14 : self%couple_ssrd = self%couple_cfg%active() .and. self%couple_cfg%meteo_expect_ssrd
428 14 : self%couple_strd = self%couple_cfg%active() .and. self%couple_cfg%meteo_expect_strd
429 14 : self%couple_tann = self%couple_cfg%active() .and. self%couple_cfg%meteo_expect_tann
430 14 : if (self%couple_cfg%active()) then
431 0 : self%couple_step_delta = timedelta(hours=self%couple_cfg%meteo_timestep)
432 0 : self%couple_is_hourly = self%couple_step_delta == one_hour()
433 : ! default init values for coupling times: 0001-01-01
434 0 : if (self%couple_cfg%meteo_expect_pre) self%couple_pre_time = datetime()
435 0 : if (self%couple_cfg%meteo_expect_temp) self%couple_temp_time = datetime()
436 : ! PET releated
437 0 : if (self%couple_cfg%meteo_expect_pet) self%couple_pet_time = datetime()
438 0 : if (self%couple_cfg%meteo_expect_tmin) self%couple_tmin_time = datetime()
439 0 : if (self%couple_cfg%meteo_expect_tmax) self%couple_tmax_time = datetime()
440 0 : if (self%couple_cfg%meteo_expect_netrad) self%couple_netrad_time = datetime()
441 0 : if (self%couple_cfg%meteo_expect_absvappress) self%couple_absvappress_time = datetime()
442 0 : if (self%couple_cfg%meteo_expect_windspeed) self%couple_windspeed_time = datetime()
443 : ! RIV-TEMP releated
444 0 : if (self%couple_cfg%meteo_expect_ssrd) self%couple_ssrd_time = datetime()
445 0 : if (self%couple_cfg%meteo_expect_strd) self%couple_strd_time = datetime()
446 0 : if (self%couple_cfg%meteo_expect_tann) self%couple_tann_time = datetime()
447 : ! PET related meteo
448 0 : self%couple_all = self%couple_cfg%meteo_expect_pre .and. self%couple_cfg%meteo_expect_temp
449 0 : select case (self%pet_case)
450 : case(-1 : 0) ! pet is input
451 0 : self%couple_all = self%couple_all .and. self%couple_cfg%meteo_expect_pet
452 0 : if (self%couple_cfg%meteo_expect_tmin) call error_message("Coupling: tmin expected but not needed for PET.")
453 0 : if (self%couple_cfg%meteo_expect_tmax) call error_message("Coupling: tmax expected but not needed for PET.")
454 0 : if (self%couple_cfg%meteo_expect_netrad) call error_message("Coupling: netrad expected but not needed for PET.")
455 0 : if (self%couple_cfg%meteo_expect_absvappress) call error_message("Coupling: absvappress expected but not needed for PET.")
456 0 : if (self%couple_cfg%meteo_expect_windspeed) call error_message("Coupling: windspeed expected but not needed for PET.")
457 :
458 : case(1) ! Hargreaves-Samani formulation (input: minimum and maximum Temperature)
459 0 : self%couple_all = self%couple_all .and. self%couple_cfg%meteo_expect_tmin
460 0 : self%couple_all = self%couple_all .and. self%couple_cfg%meteo_expect_tmax
461 0 : if (self%couple_cfg%meteo_expect_pet) call error_message("Coupling: pet expected but not needed for PET.")
462 0 : if (self%couple_cfg%meteo_expect_netrad) call error_message("Coupling: netrad expected but not needed for PET.")
463 0 : if (self%couple_cfg%meteo_expect_absvappress) call error_message("Coupling: absvappress expected but not needed for PET.")
464 0 : if (self%couple_cfg%meteo_expect_windspeed) call error_message("Coupling: windspeed expected but not needed for PET.")
465 :
466 : case(2) ! Priestley-Taylor formulation (input: net radiation)
467 0 : self%couple_all = self%couple_all .and. self%couple_cfg%meteo_expect_netrad
468 0 : if (self%couple_cfg%meteo_expect_pet) call error_message("Coupling: pet expected but not needed for PET.")
469 0 : if (self%couple_cfg%meteo_expect_tmin) call error_message("Coupling: tmin expected but not needed for PET.")
470 0 : if (self%couple_cfg%meteo_expect_tmax) call error_message("Coupling: tmax expected but not needed for PET.")
471 0 : if (self%couple_cfg%meteo_expect_absvappress) call error_message("Coupling: absvappress expected but not needed for PET.")
472 0 : if (self%couple_cfg%meteo_expect_windspeed) call error_message("Coupling: windspeed expected but not needed for PET.")
473 :
474 : case(3) ! Penman-Monteith formulation (input: net radiationm absulute vapour pressure, windspeed)
475 0 : self%couple_all = self%couple_all .and. self%couple_cfg%meteo_expect_netrad
476 0 : self%couple_all = self%couple_all .and. self%couple_cfg%meteo_expect_absvappress
477 0 : self%couple_all = self%couple_all .and. self%couple_cfg%meteo_expect_windspeed
478 0 : if (self%couple_cfg%meteo_expect_pet) call error_message("Coupling: pet expected but not needed for PET.")
479 0 : if (self%couple_cfg%meteo_expect_tmin) call error_message("Coupling: tmin expected but not needed for PET.")
480 0 : if (self%couple_cfg%meteo_expect_tmax) call error_message("Coupling: tmax expected but not needed for PET.")
481 : end select
482 : ! river temperature related meteo
483 0 : if ( self%riv_temp_case == 0 ) then
484 0 : if (self%couple_cfg%meteo_expect_ssrd) call error_message("Coupling: ssrd expected but river temperature not activated.")
485 0 : if (self%couple_cfg%meteo_expect_strd) call error_message("Coupling: strd expected but river temperature not activated.")
486 0 : if (self%couple_cfg%meteo_expect_tann) call error_message("Coupling: tann expected but river temperature not activated.")
487 : else
488 0 : self%couple_all = self%couple_all .and. self%couple_cfg%meteo_expect_ssrd
489 0 : self%couple_all = self%couple_all .and. self%couple_cfg%meteo_expect_strd
490 0 : self%couple_all = self%couple_all .and. self%couple_cfg%meteo_expect_tann
491 : end if
492 : end if
493 14 : end subroutine config
494 :
495 : !> \brief whether meteo data should be read completely at the begining
496 : !> \return True if meteo data is retrieved with a single read
497 2020989 : logical function single_read(self, iDomain)
498 : implicit none
499 : class(meteo_handler_type), intent(in) :: self
500 : integer(i4), intent(in) :: iDomain !< current domain
501 2020989 : single_read = self%timeStep_model_inputs(iDomain) == 0_i4
502 14 : end function single_read
503 :
504 : !> \brief Initialize meteo data and level-2 grid
505 28 : subroutine init_level2(self, level0, level1)
506 :
507 2020989 : use mo_grid, only : set_domain_indices
508 : use mo_common_types, only : grid
509 : use mo_grid, only : init_lowres_level
510 : use mo_read_spatial_data, only : read_header_ascii
511 : use mo_string_utils, only : num2str
512 :
513 : implicit none
514 :
515 : class(meteo_handler_type), intent(inout) :: self
516 : !> grid information at level-0
517 : type(Grid), dimension(:), intent(in) :: level0
518 : !> grid information at level-1 if all meteo data is coupled
519 : type(Grid), dimension(:), intent(in) :: level1
520 :
521 : ! header info
522 : integer(i4) :: nrows2, ncols2
523 14 : real(dp) :: xllcorner2, yllcorner2
524 14 : real(dp) :: cellsize2, nodata_dummy
525 : ! file name
526 : character(256) :: fName
527 : ! looping variable
528 : integer(i4) :: iDomain
529 :
530 : ! create level-2 info
531 68 : allocate(self%level2(self%nDomains))
532 :
533 : ! we don't need level 2 if all meteo data comes from the coupler
534 14 : if (self%couple_all) then
535 0 : self%level2(:) = level1(:)
536 : return
537 : end if
538 :
539 40 : do iDomain = 1, self%nDomains
540 : ! read header
541 26 : fName = trim(adjustl(self%dir_meteo_header(iDomain))) // trim(adjustl(self%file_meteo_header))
542 : call read_header_ascii(trim(fName), self%umeteo_header, &
543 26 : nrows2, ncols2, xllcorner2, yllcorner2, cellsize2, nodata_dummy)
544 : ! check grid compatibility
545 26 : call init_lowres_level(level0(self%L0DataFrom(iDomain)), cellsize2, self%level2(iDomain))
546 : ! check
547 0 : if ((ncols2 .ne. self%level2(iDomain)%ncols) .or. &
548 : (nrows2 .ne. self%level2(iDomain)%nrows) .or. &
549 : (abs(xllcorner2 - self%level2(iDomain)%xllcorner) .gt. tiny(1.0_dp)) .or. &
550 26 : (abs(yllcorner2 - self%level2(iDomain)%yllcorner) .gt. tiny(1.0_dp)) .or. &
551 40 : (abs(cellsize2 - self%level2(iDomain)%cellsize) .gt. tiny(1.0_dp))) then
552 : call error_message(' ***ERROR: subroutine L2_variable_init: size mismatch in grid file for level2 in domain ', &
553 0 : trim(adjustl(num2str(iDomain))), '!', raise=.false.)
554 0 : call error_message(' Expected to have following properties (based on L0):', raise=.false.)
555 0 : call error_message('... rows: ', trim(adjustl(num2str(self%level2(iDomain)%nrows))), ', ', raise=.false.)
556 0 : call error_message('... cols: ', trim(adjustl(num2str(self%level2(iDomain)%ncols))), ', ', raise=.false.)
557 0 : call error_message('... cellsize: ', trim(adjustl(num2str(self%level2(iDomain)%cellsize))), ', ', raise=.false.)
558 0 : call error_message('... xllcorner:', trim(adjustl(num2str(self%level2(iDomain)%xllcorner))), ', ', raise=.false.)
559 0 : call error_message('... yllcorner:', trim(adjustl(num2str(self%level2(iDomain)%yllcorner))), ', ', raise=.false.)
560 0 : call error_message(' Provided (in precipitation file):', raise=.false.)
561 0 : call error_message('... rows: ', trim(adjustl(num2str(nrows2))), ', ', raise=.false.)
562 0 : call error_message('... cols: ', trim(adjustl(num2str(ncols2))), ', ', raise=.false.)
563 0 : call error_message('... cellsize: ', trim(adjustl(num2str(cellsize2))), ', ', raise=.false.)
564 0 : call error_message('... xllcorner:', trim(adjustl(num2str(xllcorner2))), ', ', raise=.false.)
565 0 : call error_message('... yllcorner:', trim(adjustl(num2str(yllcorner2))), ', ')
566 : end if
567 : end do
568 :
569 : ! set indices
570 14 : call set_domain_indices(self%level2)
571 :
572 14 : end subroutine init_level2
573 :
574 : !> \brief update the current time-step of the \ref meteo_handler_type class
575 4040544 : subroutine update_timestep(self, tt, time, iDomain, level1, simPer)
576 :
577 : implicit none
578 :
579 : class(meteo_handler_type), intent(inout) :: self
580 : integer(i4), intent(in) :: tt !< current time step
581 : real(dp), intent(in) :: time !< current decimal Julian day
582 : integer(i4), intent(in) :: iDomain !< current domain
583 : !> grid information at hydrologic level
584 : type(Grid), dimension(:), intent(in) :: level1
585 : !> warmPer + evalPer
586 : type(period), dimension(:), intent(in) :: simPer
587 :
588 : ! store current indizes and time
589 : ! only needed for the "get_<var>" methods
590 2020272 : self%s1 = level1(iDomain)%iStart
591 2020272 : self%e1 = level1(iDomain)%iEnd
592 2020272 : self%time = time
593 :
594 : ! time increment is done right after call to mrm (and initially before looping)
595 2020272 : if (self%single_read(iDomain) .or. self%couple_all) then
596 : ! whole meteorology is already read or all meteo is coupled
597 :
598 : ! set start and end of meteo position
599 1954752 : self%s_meteo = level1(iDomain)%iStart
600 1954752 : self%e_meteo = level1(iDomain)%iEnd
601 :
602 : ! time step for meteorological variable (daily values)
603 : ! iMeteoTS = ceiling(real(tt, dp) / real(nTstepDay, dp))
604 1954752 : if (self%couple_all) then
605 0 : self%iMeteoTS = 1_i4
606 : else
607 1954752 : self%iMeteoTS = ceiling(real(tt, dp) / real(nint( 24._dp / real(self%nTstepForcingDay, dp)), dp))
608 : end if
609 : else
610 : ! read chunk of meteorological forcings data (reading, upscaling/downscaling)
611 65520 : call self%prepare_data(tt, iDomain, level1, simPer)
612 : ! set start and end of meteo position
613 65520 : self%s_meteo = 1
614 65520 : self%e_meteo = level1(iDomain)%iEnd - level1(iDomain)%iStart + 1
615 : ! time step for meteorological variable (daily values)
616 : self%iMeteoTS = ceiling(real(tt, dp) / real(nint( 24._dp / real(self%nTstepForcingDay, dp)), dp)) &
617 65520 : - (self%readPer%julStart - simPer(iDomain)%julStart)
618 : end if
619 :
620 14 : end subroutine update_timestep
621 :
622 : !> \brief Prepare meteorological forcings data for a given variable
623 : !> \details Prepare meteorological forcings data for a given variable.
624 : !! Internally this subroutine calls another routine meteo_wrapper
625 : !! for different meterological variables
626 : !> \changelog
627 : !! - Matthias Zink, Jun 2013
628 : !! - addded NetCDf reader
629 : !! - Rohini Kumar, Aug 2013
630 : !! - name changed "inputFormat" to inputFormat_meteo_forcings
631 : !! - Matthias Zink, Feb 2014
632 : !! - added read in for different PET processes (process 5)
633 : !! - Stephan Thober, Jun 2014
634 : !! - add chunk_config for chunk read, copied L2 initialization to mo_startup
635 : !! - Stephan Thober, Nov 2016
636 : !! - moved processMatrix to common variables
637 : !! - Stephan Thober, Jan 2017
638 : !! - added subroutine for meteo_weights
639 : !! - Robert Schweppe Jun 2018
640 : !! - refactoring and reformatting
641 : !! - Sebastian Müller Mar 2023
642 : !! - converted routine to meteo-handler method
643 : !> \authors Rohini Kumar
644 : !> \date Jan 2013
645 131082 : subroutine prepare_data(self, tt, iDomain, level1, simPer)
646 :
647 2020272 : use mo_string_utils, only : num2str
648 : use mo_timer, only : timer_get, timer_start, timer_stop
649 : use mo_meteo_helper, only : meteo_forcings_wrapper, meteo_weights_wrapper, chunk_config
650 :
651 : implicit none
652 :
653 : class(meteo_handler_type), intent(inout) :: self
654 : integer(i4), intent(in) :: tt !< current timestep
655 : integer(i4), intent(in) :: iDomain !< Domain number
656 : !> grid information at hydrologic level
657 : type(Grid), dimension(:), intent(in) :: level1
658 : !> warmPer + evalPer
659 : type(period), dimension(:), intent(in) :: simPer
660 :
661 : ! indicate whether data should be read
662 : logical :: read_flag
663 : integer(i4) :: domainID ! current domain ID
664 :
665 : ! allocate arrays only once if they are coupled
666 65541 : if (self%couple_pre .and. .not. allocated(self%L1_pre)) allocate(self%L1_pre(level1(iDomain)%nCells, 1))
667 65541 : if (self%couple_temp .and. .not. allocated(self%L1_temp)) allocate(self%L1_temp(level1(iDomain)%nCells, 1))
668 65541 : if (self%couple_pet .and. .not. allocated(self%L1_pet)) allocate(self%L1_pet(level1(iDomain)%nCells, 1))
669 65541 : if (self%couple_tmin .and. .not. allocated(self%L1_tmin)) allocate(self%L1_tmin(level1(iDomain)%nCells, 1))
670 65541 : if (self%couple_tmax .and. .not. allocated(self%L1_tmax)) allocate(self%L1_tmax(level1(iDomain)%nCells, 1))
671 65541 : if (self%couple_netrad .and. .not. allocated(self%L1_netrad)) allocate(self%L1_netrad(level1(iDomain)%nCells, 1))
672 65541 : if (self%couple_absvappress .and. .not. allocated(self%L1_absvappress)) allocate(self%L1_absvappress(level1(iDomain)%nCells, 1))
673 65541 : if (self%couple_windspeed .and. .not. allocated(self%L1_windspeed)) allocate(self%L1_windspeed(level1(iDomain)%nCells, 1))
674 65541 : if (self%couple_ssrd .and. .not. allocated(self%L1_ssrd)) allocate(self%L1_ssrd(level1(iDomain)%nCells, 1))
675 65541 : if (self%couple_strd .and. .not. allocated(self%L1_strd)) allocate(self%L1_strd(level1(iDomain)%nCells, 1))
676 65541 : if (self%couple_tann .and. .not. allocated(self%L1_tann)) allocate(self%L1_tann(level1(iDomain)%nCells, 1))
677 :
678 65541 : domainID = self%indices(iDomain)
679 :
680 : ! configuration of chunk_read
681 65541 : call chunk_config(iDomain, tt, self%nTstepDay, simPer, self%timestep, self%timeStep_model_inputs, read_flag, self%readPer)
682 :
683 : ! only read, if read_flag is true
684 65541 : if (read_flag) then
685 :
686 : ! read weights for hourly disaggregation
687 111 : if (tt .eq. 1) then
688 : ! TODO-RIV-TEMP: No NC files for weights for radiation at the moment
689 26 : if (self%single_read(iDomain)) call message(' read meteo weights for tavg ...')
690 0 : call meteo_weights_wrapper(iDomain, self%read_meteo_weights, self%dirTemperature(iDomain), &
691 26 : self%L1_temp_weights, level1=level1, level2=self%level2, ncvarName = 'tavg_weight')
692 :
693 26 : if (self%single_read(iDomain)) call message(' read meteo weights for pet ...')
694 0 : call meteo_weights_wrapper(iDomain, self%read_meteo_weights, self%dirReferenceET(iDomain), &
695 26 : self%L1_pet_weights, level1=level1, level2=self%level2, ncvarName = 'pet_weight')
696 :
697 26 : if (self%single_read(iDomain)) call message(' read meteo weights for pre ...')
698 0 : call meteo_weights_wrapper(iDomain, self%read_meteo_weights, self%dirPrecipitation(iDomain), &
699 26 : self%L1_pre_weights, level1=level1, level2=self%level2, ncvarName = 'pre_weight')
700 : end if
701 :
702 : ! free L1 variables if chunk read is activated
703 111 : if (self%timeStep_model_inputs(iDomain) .ne. 0) then
704 90 : if (.not. self%couple_pre .and. allocated(self%L1_pre)) deallocate(self%L1_pre)
705 90 : if (.not. self%couple_temp .and. allocated(self%L1_temp)) deallocate(self%L1_temp)
706 90 : if (.not. self%couple_pet .and. allocated(self%L1_pet)) deallocate(self%L1_pet)
707 90 : if (.not. self%couple_tmin .and. allocated(self%L1_tmin)) deallocate(self%L1_tmin)
708 90 : if (.not. self%couple_tmax .and. allocated(self%L1_tmax)) deallocate(self%L1_tmax)
709 90 : if (.not. self%couple_netrad .and. allocated(self%L1_netrad)) deallocate(self%L1_netrad)
710 90 : if (.not. self%couple_absvappress .and. allocated(self%L1_absvappress)) deallocate(self%L1_absvappress)
711 90 : if (.not. self%couple_windspeed .and. allocated(self%L1_windspeed)) deallocate(self%L1_windspeed)
712 : end if
713 :
714 : ! Domain characteristics and read meteo header
715 111 : if (self%single_read(iDomain)) then
716 21 : call message(' Reading meteorological forcings for Domain: ', trim(adjustl(num2str(domainID))), ' ...')
717 21 : call timer_start(1)
718 : end if
719 :
720 : ! precipitation
721 111 : if (.not. self%couple_pre) then
722 111 : if (self%single_read(iDomain)) call message(' read precipitation ...')
723 : ! upper bound: 1825 mm/d in La Réunion 7-8 Jan 1966
724 0 : call meteo_forcings_wrapper(iDomain, self%dirPrecipitation(iDomain), self%inputFormat_meteo_forcings, &
725 : dataOut1=self%L1_pre, &
726 0 : readPer=self%readPer, nTstepForcingDay=self%nTstepForcingDay, level1=level1, level2=self%level2, &
727 111 : lower = 0.0_dp, upper = 2000._dp, ncvarName = 'pre', bound_error=self%bound_error)
728 : end if
729 :
730 : ! temperature
731 111 : if (.not. self%couple_temp) then
732 111 : if (self%single_read(iDomain)) call message(' read temperature ...')
733 0 : call meteo_forcings_wrapper(iDomain, self%dirTemperature(iDomain), self%inputFormat_meteo_forcings, &
734 : dataOut1=self%L1_temp, &
735 0 : readPer=self%readPer, nTstepForcingDay=self%nTstepForcingDay, level1=level1, level2=self%level2, &
736 111 : lower = -100._dp, upper = 100._dp, ncvarName = 'tavg', bound_error=self%bound_error)
737 : end if
738 :
739 : ! read input for PET (process 5) depending on specified option
740 : ! 0 - input, 1 - Hargreaves-Samani, 2 - Priestley-Taylor, 3 - Penman-Monteith
741 201 : select case (self%pet_case)
742 : case(-1 : 0) ! pet is input
743 90 : if (.not. self%couple_pet) then
744 90 : if (self%single_read(iDomain)) call message(' read pet ...')
745 0 : call meteo_forcings_wrapper(iDomain, self%dirReferenceET(iDomain), self%inputFormat_meteo_forcings, &
746 : dataOut1=self%L1_pet, &
747 0 : readPer=self%readPer, nTstepForcingDay=self%nTstepForcingDay, level1=level1, level2=self%level2, &
748 90 : lower = 0.0_dp, upper = 1000._dp, ncvarName = 'pet', bound_error=self%bound_error)
749 : end if
750 : ! allocate PET and dummies for mhm_call
751 90 : if ((iDomain.eq.self%nDomains) .OR. (self%timeStep_model_inputs(iDomain) .NE. 0)) then
752 78 : allocate(self%L1_tmin(1, 1))
753 78 : allocate(self%L1_tmax(1, 1))
754 78 : allocate(self%L1_netrad(1, 1))
755 78 : allocate(self%L1_absvappress(1, 1))
756 78 : allocate(self%L1_windspeed(1, 1))
757 : end if
758 :
759 : case(1) ! Hargreaves-Samani formulation (input: minimum and maximum Temperature)
760 2 : if (.not. self%couple_tmin) then
761 2 : if (self%single_read(iDomain)) call message(' read min. temperature ...')
762 0 : call meteo_forcings_wrapper(iDomain, self%dirMinTemperature(iDomain), self%inputFormat_meteo_forcings, &
763 : dataOut1=self%L1_tmin, &
764 0 : readPer=self%readPer, nTstepForcingDay=self%nTstepForcingDay, level1=level1, level2=self%level2, &
765 2 : lower = -100.0_dp, upper = 100._dp, ncvarName = 'tmin', bound_error=self%bound_error)
766 : end if
767 2 : if (.not. self%couple_tmax) then
768 2 : if (self%single_read(iDomain)) call message(' read max. temperature ...')
769 0 : call meteo_forcings_wrapper(iDomain, self%dirMaxTemperature(iDomain), self%inputFormat_meteo_forcings, &
770 : dataOut1=self%L1_tmax, &
771 0 : readPer=self%readPer, nTstepForcingDay=self%nTstepForcingDay, level1=level1, level2=self%level2, &
772 2 : lower = -100.0_dp, upper = 100._dp, ncvarName = 'tmax', bound_error=self%bound_error)
773 : end if
774 : ! allocate PET and dummies for mhm_call
775 2 : if ((iDomain .eq. self%nDomains) .OR. (self%timeStep_model_inputs(iDomain) .NE. 0)) then
776 8 : allocate(self%L1_pet (size(self%L1_tmax, dim = 1), size(self%L1_tmax, dim = 2)))
777 2 : allocate(self%L1_netrad(1, 1))
778 2 : allocate(self%L1_absvappress(1, 1))
779 2 : allocate(self%L1_windspeed(1, 1))
780 : end if
781 :
782 : case(2) ! Priestley-Taylor formulation (input: net radiation)
783 18 : if (.not. self%couple_netrad) then
784 18 : if (self%single_read(iDomain)) call message(' read net radiation ...')
785 0 : call meteo_forcings_wrapper(iDomain, self%dirNetRadiation(iDomain), self%inputFormat_meteo_forcings, &
786 : dataOut1=self%L1_netrad, &
787 0 : readPer=self%readPer, nTstepForcingDay=self%nTstepForcingDay, level1=level1, level2=self%level2, &
788 18 : lower = -500.0_dp, upper = 1500._dp, ncvarName = 'net_rad', bound_error=self%bound_error)
789 : end if
790 : ! allocate PET and dummies for mhm_call
791 18 : if ((iDomain .eq. self%nDomains) .OR. (self%timeStep_model_inputs(iDomain) .NE. 0)) then
792 72 : allocate(self%L1_pet (size(self%L1_netrad, dim = 1), size(self%L1_netrad, dim = 2)))
793 18 : allocate(self%L1_tmin(1, 1))
794 18 : allocate(self%L1_tmax(1, 1))
795 18 : allocate(self%L1_absvappress(1, 1))
796 18 : allocate(self%L1_windspeed(1, 1))
797 : end if
798 :
799 : case(3) ! Penman-Monteith formulation (input: net radiationm absulute vapour pressure, windspeed)
800 1 : if (.not. self%couple_netrad) then
801 1 : if (self%single_read(iDomain)) call message(' read net radiation ...')
802 0 : call meteo_forcings_wrapper(iDomain, self%dirNetRadiation(iDomain), self%inputFormat_meteo_forcings, &
803 : dataOut1=self%L1_netrad, &
804 0 : readPer=self%readPer, nTstepForcingDay=self%nTstepForcingDay, level1=level1, level2=self%level2, &
805 1 : lower = -500.0_dp, upper = 1500._dp, ncvarName = 'net_rad', bound_error=self%bound_error)
806 : end if
807 1 : if (.not. self%couple_absvappress) then
808 1 : if (self%single_read(iDomain)) call message(' read absolute vapour pressure ...')
809 0 : call meteo_forcings_wrapper(iDomain, self%dirabsVapPressure(iDomain), self%inputFormat_meteo_forcings, &
810 : dataOut1=self%L1_absvappress, &
811 0 : readPer=self%readPer, nTstepForcingDay=self%nTstepForcingDay, level1=level1, level2=self%level2, &
812 1 : lower = 0.0_dp, upper = 15000.0_dp, ncvarName = 'eabs', bound_error=self%bound_error)
813 : end if
814 1 : if (.not. self%couple_windspeed) then
815 1 : if (self%single_read(iDomain)) call message(' read windspeed ...')
816 0 : call meteo_forcings_wrapper(iDomain, self%dirwindspeed(iDomain), self%inputFormat_meteo_forcings, &
817 : dataOut1=self%L1_windspeed, &
818 0 : readPer=self%readPer, nTstepForcingDay=self%nTstepForcingDay, level1=level1, level2=self%level2, &
819 1 : lower = 0.0_dp, upper = 250.0_dp, ncvarName = 'windspeed', bound_error=self%bound_error)
820 : end if
821 : ! allocate PET and dummies for mhm_call
822 112 : if ((iDomain.eq.self%nDomains) .OR. (self%timeStep_model_inputs(iDomain) .NE. 0)) then
823 4 : allocate(self%L1_pet (size(self%L1_absvappress, dim = 1), size(self%L1_absvappress, dim = 2)))
824 1 : allocate(self%L1_tmin(1, 1))
825 1 : allocate(self%L1_tmax(1, 1))
826 : end if
827 : end select
828 :
829 : ! long/short-wave radiation and annual mean temperature for river-temperature routing
830 111 : if ( self%riv_temp_case .ne. 0 ) then
831 : ! free L1 variables if chunk read is activated
832 18 : if (self%timeStep_model_inputs(iDomain) .ne. 0) then
833 18 : if (.not. self%couple_ssrd .and. allocated(self%L1_ssrd)) deallocate(self%L1_ssrd)
834 18 : if (.not. self%couple_strd .and. allocated(self%L1_strd)) deallocate(self%L1_strd)
835 18 : if (.not. self%couple_tann .and. allocated(self%L1_tann)) deallocate(self%L1_tann)
836 : end if
837 18 : if (.not. self%couple_ssrd) then
838 18 : if (self%single_read(iDomain)) call message(' read short-wave radiation ...')
839 : call meteo_forcings_wrapper( &
840 0 : iDomain, self%dirRadiation(iDomain), self%inputFormat_meteo_forcings, &
841 : dataOut1=self%L1_ssrd, &
842 0 : readPer=self%readPer, nTstepForcingDay=self%nTstepForcingDay, level1=level1, level2=self%level2, &
843 18 : lower = 0.0_dp, upper = 1500._dp, ncvarName = 'ssrd', bound_error=self%bound_error)
844 : end if
845 18 : if (.not. self%couple_strd) then
846 18 : if (self%single_read(iDomain)) call message(' read long-wave radiation ...')
847 : call meteo_forcings_wrapper( &
848 0 : iDomain, self%dirRadiation(iDomain), self%inputFormat_meteo_forcings, &
849 : dataOut1=self%L1_strd, &
850 0 : readPer=self%readPer, nTstepForcingDay=self%nTstepForcingDay, level1=level1, level2=self%level2, &
851 18 : lower = 0.0_dp, upper = 1500._dp, ncvarName = 'strd', bound_error=self%bound_error)
852 : end if
853 18 : if (.not. self%couple_tann) then
854 18 : if (self%single_read(iDomain)) call message(' read annual mean temperature ...')
855 : call meteo_forcings_wrapper( &
856 0 : iDomain, self%dirTemperature(iDomain), self%inputFormat_meteo_forcings, &
857 : dataOut1=self%L1_tann, &
858 0 : readPer=self%readPer, nTstepForcingDay=self%nTstepForcingDay, level1=level1, level2=self%level2, &
859 18 : lower = -100.0_dp, upper = 100._dp, ncvarName = 'tann', bound_error=self%bound_error)
860 : end if
861 : end if
862 :
863 111 : if (self%single_read(iDomain)) then
864 21 : call timer_stop(1)
865 21 : call message(' in ', trim(num2str(timer_get(1), '(F9.3)')), ' seconds.')
866 : end if
867 : end if
868 :
869 65541 : if (tt .eq. 1) then
870 : ! set hourly flags (once at begining)
871 26 : if (self%couple_all) then
872 0 : self%nTstepForcingDay = int(one_day() / self%couple_step_delta, i4)
873 0 : self%is_hourly_forcing = self%couple_is_hourly
874 : else
875 26 : self%is_hourly_forcing = (self%nTstepForcingDay .eq. 24_i4)
876 : end if
877 :
878 : ! check hourly data for PET (once at begining)
879 28 : select case (self%pet_case)
880 : case(1) ! Hargreaves-Samani formulation
881 2 : if (self%couple_temp .and. self%couple_tmin .and. self%couple_tmax) then
882 0 : if (self%couple_is_hourly) call error_message("Coupling: PET - Hargreaves-Samani needs daily data. Got hourly.")
883 2 : else if (self%couple_temp .or. self%couple_tmin .or. self%couple_tmax) then
884 0 : if (self%couple_is_hourly) call error_message("Coupling: PET - Hargreaves-Samani needs daily data. Got hourly.")
885 0 : if (self%is_hourly_forcing) call error_message("Meteo: PET - Hargreaves-Samani needs daily data. Got hourly.")
886 : else
887 2 : if (self%is_hourly_forcing) call error_message("Meteo: PET - Hargreaves-Samani needs daily data. Got hourly.")
888 : end if
889 :
890 : case(2) ! Priestley-Taylor formulation
891 1 : if (self%couple_temp .and. self%couple_netrad) then
892 0 : if (self%couple_is_hourly) call error_message("Coupling: PET - Priestley-Taylor needs daily data. Got hourly.")
893 1 : else if (self%couple_temp .or. self%couple_netrad) then
894 0 : if (self%couple_is_hourly) call error_message("Coupling: PET - Priestley-Taylor needs daily data. Got hourly.")
895 0 : if (self%is_hourly_forcing) call error_message("Meteo: PET - Priestley-Taylor needs daily data. Got hourly.")
896 : else
897 1 : if (self%is_hourly_forcing) call error_message("Meteo: PET - Priestley-Taylor needs daily data. Got hourly.")
898 : end if
899 :
900 : case(3) ! Penman-Monteith formulation
901 26 : if (self%couple_temp .and. self%couple_netrad .and. self%couple_absvappress .and. self%couple_windspeed) then
902 0 : if (self%couple_is_hourly) call error_message("Coupling: PET - Penman-Monteith needs daily data. Got hourly.")
903 1 : else if (self%couple_temp .or. self%couple_netrad .or. self%couple_absvappress .or. self%couple_windspeed) then
904 0 : if (self%couple_is_hourly) call error_message("Coupling: PET - Penman-Monteith needs daily data. Got hourly.")
905 0 : if (self%is_hourly_forcing) call error_message("Meteo: PET - Penman-Monteith needs daily data. Got hourly.")
906 : else
907 1 : if (self%is_hourly_forcing) call error_message("Meteo: PET - Penman-Monteith needs daily data. Got hourly.")
908 : end if
909 : end select
910 : end if
911 :
912 65541 : end subroutine prepare_data
913 :
914 : !> \brief get corrected PET for the current timestep and domain
915 2020272 : subroutine get_corrected_pet(self, pet_calc, &
916 2020272 : petLAIcorFactorL1, fAsp, HarSamCoeff, latitude, PrieTayAlpha, aeroResist, surfResist)
917 :
918 65541 : use mo_mhm_constants, only : HarSamConst
919 : use mo_julian, only : date2dec, dec2date
920 : use mo_meteo_temporal_tools, only : temporal_disagg_meteo_weights, temporal_disagg_flux_daynight
921 : use mo_string_utils, only : num2str
922 : use mo_pet, only : pet_hargreaves, pet_penman, pet_priestly
923 :
924 : implicit none
925 :
926 : class(meteo_handler_type), intent(inout) :: self
927 : !> [mm TS-1] estimated PET (if PET is input = corrected values (fAsp*PET))
928 : real(dp), dimension(:), intent(inout) :: pet_calc
929 : !> PET correction factor based on LAI at level 1
930 : real(dp), dimension(:), intent(in) :: petLAIcorFactorL1
931 : !> [1] PET correction for Aspect at level 1
932 : real(dp), dimension(:), intent(in) :: fAsp
933 : !> [1] PET Hargreaves Samani coefficient at level 1
934 : real(dp), dimension(:), intent(in) :: HarSamCoeff
935 : !> latitude on level 1
936 : real(dp), dimension(:), intent(in) :: latitude
937 : !> [1] PET Priestley Taylor coefficient at level 1
938 : real(dp), dimension(:), intent(in) :: PrieTayAlpha
939 : !> [s m-1] PET aerodynamical resitance at level 1
940 : real(dp), dimension(:), intent(in) :: aeroResist
941 : !> [s m-1] PET bulk surface resitance at level 1
942 : real(dp), dimension(:), intent(in) :: surfResist
943 :
944 : ! pet in [mm d-1]
945 2020272 : real(dp) :: pet
946 : logical :: isday, is_hourly
947 : integer(i4) :: year, month, day, hour
948 : type(datetime) :: curr_dt
949 : type(timedelta) :: meteo_time_delta
950 :
951 : ! doy of the year [1-365 or 1-366]
952 : integer(i4) :: doy
953 :
954 : ! number of L1 cells
955 : integer(i4) :: nCells1
956 : ! cell index
957 : integer(i4) :: k, i, s1
958 : ! individual meteo time-steps for all variables due to possible coupling
959 : integer(i4) :: mTS_pet, mTS_temp, mTS_tmin, mTS_tmax, mTS_rn, mTS_avp, mTS_ws
960 :
961 : ! date and month of this timestep
962 2020272 : call dec2date(self%time, yy = year, mm = month, dd = day, hh = hour)
963 2020272 : curr_dt = datetime(year, month, day, hour)
964 2020272 : doy = curr_dt%doy()
965 :
966 2020272 : nCells1 = self%e1 - self%s1 + 1
967 2020272 : s1 = self%s1
968 :
969 2020272 : if (self%couple_pet) then
970 0 : meteo_time_delta = curr_dt - self%couple_pet_time
971 : ! check that the PET from the interface has the correct time-stamp
972 0 : if (meteo_time_delta < zero_delta() .or. meteo_time_delta >= self%couple_step_delta) &
973 0 : call error_message("meteo_handler: PET was expected from coupler, but has a wrong time-stamp.")
974 0 : mTS_pet = 1_i4
975 0 : is_hourly = self%couple_is_hourly
976 : else
977 2020272 : mTS_pet = self%iMeteoTS
978 2020272 : is_hourly = self%is_hourly_forcing ! not needed to set with other variables
979 : end if
980 :
981 2020272 : if (self%couple_temp) then
982 0 : meteo_time_delta = curr_dt - self%couple_temp_time
983 : ! check that the temperature from the interface has the correct time-stamp
984 0 : if (meteo_time_delta < zero_delta() .or. meteo_time_delta >= self%couple_step_delta) &
985 0 : call error_message("meteo_handler: temperature was expected from coupler, but has a wrong time-stamp.")
986 0 : mTS_temp = 1_i4
987 0 : is_hourly = .false.
988 : else
989 2020272 : mTS_temp = self%iMeteoTS
990 : end if
991 :
992 2020272 : if (self%couple_tmin) then
993 0 : meteo_time_delta = curr_dt - self%couple_tmin_time
994 : ! check that the tmin from the interface has the correct time-stamp
995 0 : if (meteo_time_delta < zero_delta() .or. meteo_time_delta >= self%couple_step_delta) &
996 0 : call error_message("meteo_handler: minimum temperature was expected from coupler, but has a wrong time-stamp.")
997 0 : mTS_tmin = 1_i4
998 0 : is_hourly = .false.
999 : else
1000 2020272 : mTS_tmin = self%iMeteoTS
1001 : end if
1002 :
1003 2020272 : if (self%couple_tmax) then
1004 0 : meteo_time_delta = curr_dt - self%couple_tmax_time
1005 : ! check that the tmax from the interface has the correct time-stamp
1006 0 : if (meteo_time_delta < zero_delta() .or. meteo_time_delta >= self%couple_step_delta) &
1007 0 : call error_message("meteo_handler: maximum temperature was expected from coupler, but has a wrong time-stamp.")
1008 0 : mTS_tmax = 1_i4
1009 0 : is_hourly = .false.
1010 : else
1011 2020272 : mTS_tmax = self%iMeteoTS
1012 : end if
1013 :
1014 2020272 : if (self%couple_netrad) then
1015 0 : meteo_time_delta = curr_dt - self%couple_netrad_time
1016 : ! check that the netrad from the interface has the correct time-stamp
1017 0 : if (meteo_time_delta < zero_delta() .or. meteo_time_delta >= self%couple_step_delta) &
1018 0 : call error_message("meteo_handler: net-radiation was expected from coupler, but has a wrong time-stamp.")
1019 0 : mTS_rn = 1_i4
1020 0 : is_hourly = .false.
1021 : else
1022 2020272 : mTS_rn = self%iMeteoTS
1023 : end if
1024 :
1025 2020272 : if (self%couple_absvappress) then
1026 0 : meteo_time_delta = curr_dt - self%couple_absvappress_time
1027 : ! check that the absvappress from the interface has the correct time-stamp
1028 0 : if (meteo_time_delta < zero_delta() .or. meteo_time_delta >= self%couple_step_delta) &
1029 0 : call error_message("meteo_handler: Abs. vapour pressure was expected from coupler, but has a wrong time-stamp.")
1030 0 : mTS_avp = 1_i4
1031 0 : is_hourly = .false.
1032 : else
1033 2020272 : mTS_avp = self%iMeteoTS
1034 : end if
1035 :
1036 2020272 : if (self%couple_windspeed) then
1037 0 : meteo_time_delta = curr_dt - self%couple_windspeed_time
1038 : ! check that the windspeed from the interface has the correct time-stamp
1039 0 : if (meteo_time_delta < zero_delta() .or. meteo_time_delta >= self%couple_step_delta) &
1040 0 : call error_message("meteo_handler: windspeed was expected from coupler, but has a wrong time-stamp.")
1041 0 : mTS_ws = 1_i4
1042 0 : is_hourly = .false.
1043 : else
1044 2020272 : mTS_ws = self%iMeteoTS
1045 : end if
1046 :
1047 : ! flag for day or night depending on hours of the day
1048 2020272 : isday = (hour .gt. 6) .AND. (hour .le. 18)
1049 :
1050 : !$OMP parallel default(shared) &
1051 : !$OMP private(k, pet, i)
1052 : !$OMP do SCHEDULE(STATIC)
1053 78667320 : do k = 1, nCells1
1054 :
1055 : ! correct index on concatenated arrays
1056 76647048 : i = self%s_meteo - 1 + k
1057 :
1058 : ! PET calculation
1059 77538120 : select case (self%pet_case)
1060 : case(-1) ! PET is input ! correct pet for every day only once at the first time step
1061 891072 : pet = petLAIcorFactorL1(k) * self%L1_pet(i, mTS_pet)
1062 :
1063 : case(0) ! PET is input ! correct pet for every day only once at the first time step
1064 71743704 : pet = fAsp(k) * self%L1_pet(i, mTS_pet)
1065 :
1066 : case(1) ! Hargreaves-Samani
1067 : ! estimate day of the year (doy) for approximation of the extraterrestrial radiation
1068 893520 : if (self%L1_tmax(i, mTS_tmax) .lt. self%L1_tmin(i, mTS_tmin)) &
1069 : call message('WARNING: tmax smaller than tmin at doy ', &
1070 0 : num2str(doy), ' in year ', num2str(year), ' at cell', num2str(k), '!')
1071 893520 : pet = fAsp(k) * pet_hargreaves( &
1072 893520 : HarSamCoeff=HarSamCoeff(k), &
1073 : HarSamConst=HarSamConst, &
1074 0 : tavg=self%L1_temp(i, mTS_temp), &
1075 0 : tmax=self%L1_tmax(i, mTS_tmax), &
1076 0 : tmin=self%L1_tmin(i, mTS_tmin), &
1077 893520 : latitude=latitude(k), &
1078 2680560 : doy=doy)
1079 :
1080 : case(2) ! Priestley-Taylor
1081 : ! Priestley Taylor is not defined for values netrad < 0.0_dp
1082 : pet = pet_priestly( &
1083 445536 : PrieTayParam=PrieTayAlpha(k), &
1084 0 : Rn=max(self%L1_netrad(i, mTS_rn), 0.0_dp), &
1085 891072 : tavg=self%L1_temp(i, mTS_temp))
1086 :
1087 : case(3) ! Penman-Monteith
1088 : pet = pet_penman( &
1089 0 : net_rad=max(self%L1_netrad(i, mTS_rn), 0.0_dp), &
1090 0 : tavg=self%L1_temp(i, mTS_temp), &
1091 0 : act_vap_pressure=self%L1_absvappress(i, mTS_avp) / 1000.0_dp, &
1092 2673216 : aerodyn_resistance=aeroResist(k) / self%L1_windspeed(i, mTS_ws), &
1093 0 : bulksurface_resistance=surfResist(k), &
1094 : a_s=1.0_dp, &
1095 79320264 : a_sh=1.0_dp)
1096 : end select
1097 :
1098 : ! temporal disaggreagtion of forcing variables
1099 78667320 : if (self%is_hourly_forcing) then
1100 0 : pet_calc(k) = pet
1101 : else
1102 76647048 : if (self%read_meteo_weights) then
1103 : ! all meteo forcings are disaggregated with given weights
1104 : call temporal_disagg_meteo_weights( &
1105 : meteo_val_day=pet, &
1106 0 : meteo_val_weights=self%L1_pet_weights(s1 - 1 + k, month, hour + 1), &
1107 0 : meteo_val=pet_calc(k))
1108 : else
1109 : ! all meteo forcings are disaggregated with day-night correction values
1110 : call temporal_disagg_flux_daynight( &
1111 : isday=isday, &
1112 : ntimesteps_day=self%nTstepDay_dp, &
1113 : meteo_val_day=pet, &
1114 76647048 : fday_meteo_val=self%fday_pet(month), &
1115 : fnight_meteo_val=self%fnight_pet(month), &
1116 153294096 : meteo_val=pet_calc(k))
1117 : end if
1118 : end if
1119 : end do
1120 : !$OMP end do
1121 : !$OMP end parallel
1122 :
1123 2020272 : end subroutine get_corrected_pet
1124 :
1125 : !> \brief get surface temperature for the current timestep and domain
1126 2020272 : subroutine get_temp(self, temp_calc)
1127 :
1128 2020272 : use mo_julian, only : dec2date
1129 : use mo_meteo_temporal_tools, only : temporal_disagg_meteo_weights, temporal_disagg_state_daynight
1130 : use mo_constants, only : T0_dp ! 273.15 - Celcius <-> Kelvin [K]
1131 :
1132 : implicit none
1133 :
1134 : class(meteo_handler_type), intent(inout) :: self
1135 : !> [degC] temperature for current time step
1136 : real(dp), dimension(:), intent(inout) :: temp_calc
1137 :
1138 : logical :: isday, is_hourly
1139 : integer(i4) :: year, month, day, hour
1140 : type(datetime) :: curr_dt
1141 : type(timedelta) :: meteo_time_delta
1142 :
1143 : ! number of L1 cells
1144 : integer(i4) :: nCells1
1145 : ! cell index
1146 : integer(i4) :: k, i, s1, mTS
1147 :
1148 : ! date and month of this timestep
1149 2020272 : call dec2date(self%time, yy=year, mm=month, dd=day, hh=hour)
1150 :
1151 2020272 : nCells1 = self%e1 - self%s1 + 1
1152 2020272 : s1 = self%s1
1153 2020272 : if (self%couple_temp) then
1154 0 : curr_dt = datetime(year, month, day, hour)
1155 0 : meteo_time_delta = curr_dt - self%couple_temp_time
1156 : ! check that the temperature from the interface has the correct time-stamp
1157 0 : if (meteo_time_delta < zero_delta() .or. meteo_time_delta >= self%couple_step_delta) &
1158 0 : call error_message("meteo_handler: temperature was expected from coupler, but has a wrong time-stamp.")
1159 0 : mTS = 1_i4
1160 0 : is_hourly = self%couple_is_hourly
1161 : else
1162 2020272 : mTS = self%iMeteoTS
1163 2020272 : is_hourly = self%is_hourly_forcing
1164 : end if
1165 :
1166 : ! shortcut hourly data
1167 2020272 : if (is_hourly) then
1168 0 : temp_calc(:) = self%L1_temp(self%s_meteo : self%e_meteo, mTS)
1169 : return
1170 : end if
1171 :
1172 : ! flag for day or night depending on hours of the day
1173 2020272 : isday = (hour .gt. 6) .AND. (hour .le. 18)
1174 :
1175 : !$OMP parallel default(shared) &
1176 : !$OMP private(k, i)
1177 : !$OMP do SCHEDULE(STATIC)
1178 78667320 : do k = 1, nCells1
1179 :
1180 : ! correct index on concatenated arrays
1181 76647048 : i = self%s_meteo - 1 + k
1182 :
1183 78667320 : if (self%read_meteo_weights) then
1184 : ! all meteo forcings are disaggregated with given weights
1185 : call temporal_disagg_meteo_weights( &
1186 0 : meteo_val_day=self%L1_temp(i, mTS), &
1187 0 : meteo_val_weights=self%L1_temp_weights(s1 - 1 + k, month, hour + 1), &
1188 0 : meteo_val=temp_calc(k), &
1189 0 : weights_correction=T0_dp)
1190 : else
1191 : ! all meteo forcings are disaggregated with day-night correction values
1192 : call temporal_disagg_state_daynight( &
1193 : isday=isday, &
1194 : ntimesteps_day=self%nTstepDay_dp, &
1195 0 : meteo_val_day=self%L1_temp(i, mTS), &
1196 76647048 : fday_meteo_val=self%fday_temp(month), &
1197 : fnight_meteo_val=self%fnight_temp(month), &
1198 76647048 : meteo_val=temp_calc(k), &
1199 229941144 : add_correction=.true.)
1200 : end if
1201 : end do
1202 : !$OMP end do
1203 : !$OMP end parallel
1204 :
1205 4040544 : end subroutine get_temp
1206 :
1207 : !> \brief get precipitation for the current timestep and domain
1208 2020272 : subroutine get_prec(self, prec_calc)
1209 :
1210 2020272 : use mo_julian, only : dec2date
1211 : use mo_meteo_temporal_tools, only : temporal_disagg_meteo_weights, temporal_disagg_flux_daynight
1212 :
1213 : implicit none
1214 :
1215 : class(meteo_handler_type), intent(inout) :: self
1216 : !> [mm TS-1] precipitation for current time step
1217 : real(dp), dimension(:), intent(inout) :: prec_calc
1218 :
1219 : logical :: isday, is_hourly
1220 : integer(i4) :: year, month, day, hour
1221 : type(datetime) :: curr_dt
1222 : type(timedelta) :: meteo_time_delta
1223 :
1224 : ! number of L1 cells
1225 : integer(i4) :: nCells1
1226 : ! cell index
1227 : integer(i4) :: k, i, s1, mTS
1228 :
1229 : ! date and month of this timestep
1230 2020272 : call dec2date(self%time, yy=year, mm=month, dd=day, hh=hour)
1231 :
1232 2020272 : nCells1 = self%e1 - self%s1 + 1
1233 2020272 : s1 = self%s1
1234 2020272 : if (self%couple_pre) then
1235 0 : curr_dt = datetime(year, month, day, hour)
1236 0 : meteo_time_delta = curr_dt - self%couple_pre_time
1237 : ! check that the precipitation from the interface has the correct time-stamp
1238 0 : if (meteo_time_delta < zero_delta() .or. meteo_time_delta >= self%couple_step_delta) &
1239 0 : call error_message("meteo_handler: precipitation was expected from coupler, but has a wrong time-stamp.")
1240 0 : mTS = 1_i4
1241 0 : is_hourly = self%couple_is_hourly
1242 : else
1243 2020272 : mTS = self%iMeteoTS
1244 2020272 : is_hourly = self%is_hourly_forcing
1245 : end if
1246 :
1247 : ! shortcut hourly data
1248 2020272 : if (is_hourly) then
1249 0 : prec_calc(:) = self%L1_pre(self%s_meteo : self%e_meteo, mTS)
1250 : return
1251 : end if
1252 :
1253 : ! flag for day or night depending on hours of the day
1254 2020272 : isday = (hour .gt. 6) .AND. (hour .le. 18)
1255 :
1256 : !$OMP parallel default(shared) &
1257 : !$OMP private(k, i)
1258 : !$OMP do SCHEDULE(STATIC)
1259 78667320 : do k = 1, nCells1
1260 :
1261 : ! correct index on concatenated arrays
1262 76647048 : i = self%s_meteo - 1 + k
1263 :
1264 : ! temporal disaggreagtion of forcing variables
1265 78667320 : if (self%read_meteo_weights) then
1266 : ! all meteo forcings are disaggregated with given weights
1267 : call temporal_disagg_meteo_weights( &
1268 0 : meteo_val_day=self%L1_pre(i, mTS), &
1269 0 : meteo_val_weights=self%L1_pre_weights(s1 - 1 + k, month, hour + 1), &
1270 0 : meteo_val=prec_calc(k))
1271 : else
1272 : ! all meteo forcings are disaggregated with day-night correction values
1273 : call temporal_disagg_flux_daynight( &
1274 : isday=isday, &
1275 : ntimesteps_day=self%nTstepDay_dp, &
1276 0 : meteo_val_day=self%L1_pre(i, mTS), &
1277 76647048 : fday_meteo_val=self%fday_prec(month), &
1278 : fnight_meteo_val=self%fnight_prec(month), &
1279 153294096 : meteo_val=prec_calc(k))
1280 : end if
1281 : end do
1282 : !$OMP end do
1283 : !$OMP end parallel
1284 :
1285 4040544 : end subroutine get_prec
1286 :
1287 : !> \brief get surface short-wave (solar) radiation downwards for the current timestep and domain
1288 13104 : subroutine get_ssrd(self, ssrd_calc)
1289 :
1290 2020272 : use mo_julian, only : dec2date
1291 : use mo_meteo_temporal_tools, only : temporal_disagg_state_daynight
1292 :
1293 : implicit none
1294 :
1295 : class(meteo_handler_type), intent(inout) :: self
1296 : !> [W m2] surface short-wave (solar) radiation downwards for current time step
1297 : real(dp), dimension(:), intent(inout) :: ssrd_calc
1298 :
1299 : logical :: isday, is_hourly
1300 : integer(i4) :: year, month, day, hour
1301 : type(datetime) :: curr_dt
1302 : type(timedelta) :: meteo_time_delta
1303 :
1304 : ! number of L1 cells
1305 : integer(i4) :: nCells1
1306 : ! cell index
1307 : integer(i4) :: k, i, s1, mTS
1308 :
1309 : ! date and month of this timestep
1310 13104 : call dec2date(self%time, yy=year, mm=month, dd=day, hh=hour)
1311 :
1312 13104 : nCells1 = self%e1 - self%s1 + 1
1313 13104 : s1 = self%s1
1314 13104 : if (self%couple_ssrd) then
1315 0 : curr_dt = datetime(year, month, day, hour)
1316 0 : meteo_time_delta = curr_dt - self%couple_ssrd_time
1317 : ! check that the ssrd from the interface has the correct time-stamp
1318 0 : if (meteo_time_delta < zero_delta() .or. meteo_time_delta >= self%couple_step_delta) &
1319 0 : call error_message("meteo_handler: ssrd was expected from coupler, but has a wrong time-stamp.")
1320 0 : mTS = 1_i4
1321 0 : is_hourly = self%couple_is_hourly
1322 : else
1323 13104 : mTS = self%iMeteoTS
1324 13104 : is_hourly = self%is_hourly_forcing
1325 : end if
1326 :
1327 : ! shortcut hourly data
1328 13104 : if (is_hourly) then
1329 0 : ssrd_calc(:) = self%L1_ssrd(self%s_meteo : self%e_meteo, mTS)
1330 : return
1331 : end if
1332 :
1333 : ! flag for day or night depending on hours of the day
1334 13104 : isday = (hour .gt. 6) .AND. (hour .le. 18)
1335 :
1336 : !$OMP parallel default(shared) &
1337 : !$OMP private(k, i)
1338 : !$OMP do SCHEDULE(STATIC)
1339 458640 : do k = 1, nCells1
1340 : ! correct index on concatenated arrays
1341 445536 : i = self%s_meteo - 1 + k
1342 : ! TODO-RIV-TEMP: add weights for ssrd
1343 : ! temporal disaggreagtion of forcing variables
1344 : call temporal_disagg_state_daynight( &
1345 : isday=isday, &
1346 : ntimesteps_day=self%nTstepDay_dp, &
1347 0 : meteo_val_day=self%L1_ssrd(i, mTS), &
1348 445536 : fday_meteo_val=self%fday_ssrd(month), &
1349 : fnight_meteo_val=self%fnight_ssrd(month), &
1350 445536 : meteo_val=ssrd_calc(k) &
1351 1349712 : )
1352 : end do
1353 : !$OMP end do
1354 : !$OMP end parallel
1355 :
1356 26208 : end subroutine get_ssrd
1357 :
1358 : !> \brief get surface long-wave (thermal) radiation downwards for the current timestep and domain
1359 13104 : subroutine get_strd(self, strd_calc)
1360 :
1361 13104 : use mo_julian, only : dec2date
1362 : use mo_meteo_temporal_tools, only : temporal_disagg_state_daynight
1363 :
1364 : implicit none
1365 :
1366 : class(meteo_handler_type), intent(inout) :: self
1367 : !> [W m2] surface long-wave (thermal) radiation downwards for current time step
1368 : real(dp), dimension(:), intent(inout) :: strd_calc
1369 :
1370 : logical :: isday, is_hourly
1371 : integer(i4) :: year, month, day, hour
1372 : type(datetime) :: curr_dt
1373 : type(timedelta) :: meteo_time_delta
1374 :
1375 : ! number of L1 cells
1376 : integer(i4) :: nCells1
1377 : ! cell index
1378 : integer(i4) :: k, i, s1, mTS
1379 :
1380 : ! date and month of this timestep
1381 13104 : call dec2date(self%time, yy=year, mm=month, dd=day, hh=hour)
1382 :
1383 13104 : nCells1 = self%e1 - self%s1 + 1
1384 13104 : s1 = self%s1
1385 13104 : if (self%couple_strd) then
1386 0 : curr_dt = datetime(year, month, day, hour)
1387 0 : meteo_time_delta = curr_dt - self%couple_strd_time
1388 : ! check that the strd from the interface has the correct time-stamp
1389 0 : if (meteo_time_delta < zero_delta() .or. meteo_time_delta >= self%couple_step_delta) &
1390 0 : call error_message("meteo_handler: strd was expected from coupler, but has a wrong time-stamp.")
1391 0 : mTS = 1_i4
1392 0 : is_hourly = self%couple_is_hourly
1393 : else
1394 13104 : mTS = self%iMeteoTS
1395 13104 : is_hourly = self%is_hourly_forcing
1396 : end if
1397 :
1398 : ! shortcut hourly data
1399 13104 : if (is_hourly) then
1400 0 : strd_calc(:) = self%L1_strd(self%s_meteo : self%e_meteo, mTS)
1401 : return
1402 : end if
1403 :
1404 : ! flag for day or night depending on hours of the day
1405 13104 : isday = (hour .gt. 6) .AND. (hour .le. 18)
1406 :
1407 : !$OMP parallel default(shared) &
1408 : !$OMP private(k, i)
1409 : !$OMP do SCHEDULE(STATIC)
1410 458640 : do k = 1, nCells1
1411 : ! correct index on concatenated arrays
1412 445536 : i = self%s_meteo - 1 + k
1413 : ! TODO-RIV-TEMP: add weights for strd
1414 : ! temporal disaggreagtion of forcing variables
1415 : call temporal_disagg_state_daynight( &
1416 : isday=isday, &
1417 : ntimesteps_day=self%nTstepDay_dp, &
1418 0 : meteo_val_day=self%L1_strd(i, mTS), &
1419 445536 : fday_meteo_val=self%fday_strd(month), &
1420 : fnight_meteo_val=self%fnight_strd(month), &
1421 445536 : meteo_val=strd_calc(k) &
1422 1349712 : )
1423 : end do
1424 : !$OMP end do
1425 : !$OMP end parallel
1426 :
1427 26208 : end subroutine get_strd
1428 :
1429 : !> \brief get annual mean surface temperature for the current timestep and domain
1430 13104 : subroutine get_tann(self, tann_calc)
1431 13104 : use mo_julian, only : dec2date
1432 :
1433 : implicit none
1434 :
1435 : class(meteo_handler_type), intent(inout) :: self
1436 : !> [degC] annual mean air temperature
1437 : real(dp), dimension(:), intent(inout) :: tann_calc
1438 :
1439 : type(datetime) :: curr_dt
1440 : type(timedelta) :: meteo_time_delta
1441 : integer(i4) :: year, month, day, hour, mTS
1442 :
1443 13104 : if (self%couple_tann) then
1444 : ! date and month of this timestep
1445 0 : call dec2date(self%time, yy=year, mm=month, dd=day, hh=hour)
1446 0 : curr_dt = datetime(year, month, day, hour)
1447 0 : meteo_time_delta = curr_dt - self%couple_tann_time
1448 : ! check that the tann from the interface has the correct time-stamp
1449 0 : if (meteo_time_delta < zero_delta() .or. meteo_time_delta >= self%couple_step_delta) &
1450 0 : call error_message("meteo_handler: tann was expected from coupler, but has a wrong time-stamp.")
1451 0 : mTS = 1_i4
1452 : else
1453 13104 : mTS = self%iMeteoTS
1454 : end if
1455 :
1456 : ! annual temperature is not disaggregated
1457 458640 : tann_calc(:) = self%L1_tann(self%s_meteo : self%e_meteo, mTS)
1458 :
1459 13104 : end subroutine get_tann
1460 :
1461 : !> \brief set meteo_data from coupling
1462 0 : subroutine set_meteo(self, &
1463 : year, month, day, hour, &
1464 0 : pre, &
1465 0 : temp, &
1466 0 : pet, &
1467 0 : tmin, &
1468 0 : tmax, &
1469 0 : netrad, &
1470 0 : absvappress, &
1471 0 : windspeed, &
1472 0 : ssrd, &
1473 0 : strd, &
1474 0 : tann)
1475 : implicit none
1476 : class(meteo_handler_type), intent(inout) :: self
1477 : integer(i4), intent(in) :: year
1478 : integer(i4), intent(in) :: month
1479 : integer(i4), intent(in) :: day
1480 : integer(i4), intent(in), optional :: hour
1481 : real(dp), dimension(:), optional, intent(in) :: pre
1482 : real(dp), dimension(:), optional, intent(in) :: temp
1483 : real(dp), dimension(:), optional, intent(in) :: pet
1484 : real(dp), dimension(:), optional, intent(in) :: tmin
1485 : real(dp), dimension(:), optional, intent(in) :: tmax
1486 : real(dp), dimension(:), optional, intent(in) :: netrad
1487 : real(dp), dimension(:), optional, intent(in) :: absvappress
1488 : real(dp), dimension(:), optional, intent(in) :: windspeed
1489 : real(dp), dimension(:), optional, intent(in) :: ssrd
1490 : real(dp), dimension(:), optional, intent(in) :: strd
1491 : real(dp), dimension(:), optional, intent(in) :: tann
1492 :
1493 : integer(i4) :: hour_
1494 : type(datetime) :: input_time
1495 :
1496 0 : if (.not. self%couple_cfg%active()) &
1497 0 : call error_message("meteo_handler%set_meteo: coupling was not activated.")
1498 :
1499 : ! determine input time
1500 0 : hour_ = -1_i4
1501 0 : if (present(hour)) then
1502 0 : hour_ = hour
1503 0 : else if (self%couple_cfg%meteo_timestep == 24_i4) then
1504 0 : hour_ = 0_i4
1505 : end if
1506 0 : if (hour_ == -1_i4) &
1507 0 : call error_message("meteo_handler%set_meteo: hour for the meteo date needs to be given if the timestep is not daily.")
1508 0 : input_time = datetime(year, month, day, hour_)
1509 :
1510 : ! fix input time, if reference point is at the end of the time interval
1511 0 : if (self%couple_cfg%meteo_time_ref_endpoint) input_time = input_time - self%couple_step_delta
1512 :
1513 : ! check if input time matches the required time step
1514 0 : if (mod(input_time%hour, self%couple_cfg%meteo_timestep) /= 0) &
1515 0 : call error_message("meteo_handler%set_meteo: given time doesn't match couple timestep: ", input_time%str())
1516 :
1517 : ! precipitation
1518 0 : if (present(pre)) then
1519 0 : if (.not. self%couple_pre) &
1520 0 : call error_message("meteo_handler%set_meteo: precipitation was not set to be coupled.")
1521 0 : self%couple_pre_time = input_time
1522 0 : self%L1_pre(:, 1_i4) = pre(:)
1523 : end if
1524 :
1525 : ! temperature
1526 0 : if (present(temp)) then
1527 0 : if (.not. self%couple_temp) &
1528 0 : call error_message("meteo_handler%set_meteo: avg. temperature was not set to be coupled.")
1529 0 : self%couple_temp_time = input_time
1530 0 : self%L1_temp(:, 1_i4) = temp(:)
1531 : end if
1532 :
1533 : ! PET
1534 0 : if (present(pet)) then
1535 0 : if (.not. self%couple_pet) &
1536 0 : call error_message("meteo_handler%set_meteo: PET was not set to be coupled.")
1537 0 : self%couple_pet_time = input_time
1538 0 : self%L1_pet(:, 1_i4) = pet(:)
1539 : end if
1540 :
1541 : ! tmin
1542 0 : if (present(tmin)) then
1543 0 : if (.not. self%couple_tmin) &
1544 0 : call error_message("meteo_handler%set_meteo: tmin was not set to be coupled.")
1545 0 : self%couple_tmin_time = input_time
1546 0 : self%L1_tmin(:, 1_i4) = tmin(:)
1547 : end if
1548 :
1549 : ! tmax
1550 0 : if (present(tmax)) then
1551 0 : if (.not. self%couple_tmax) &
1552 0 : call error_message("meteo_handler%set_meteo: tmax was not set to be coupled.")
1553 0 : self%couple_tmax_time = input_time
1554 0 : self%L1_tmax(:, 1_i4) = tmax(:)
1555 : end if
1556 :
1557 : ! netrad
1558 0 : if (present(netrad)) then
1559 0 : if (.not. self%couple_netrad) &
1560 0 : call error_message("meteo_handler%set_meteo: netrad was not set to be coupled.")
1561 0 : self%couple_netrad_time = input_time
1562 0 : self%L1_netrad(:, 1_i4) = netrad(:)
1563 : end if
1564 :
1565 : ! absvappress
1566 0 : if (present(absvappress)) then
1567 0 : if (.not. self%couple_absvappress) &
1568 0 : call error_message("meteo_handler%set_meteo: absvappress was not set to be coupled.")
1569 0 : self%couple_absvappress_time = input_time
1570 0 : self%L1_absvappress(:, 1_i4) = absvappress(:)
1571 : end if
1572 :
1573 : ! windspeed
1574 0 : if (present(windspeed)) then
1575 0 : if (.not. self%couple_windspeed) &
1576 0 : call error_message("meteo_handler%set_meteo: windspeed was not set to be coupled.")
1577 0 : self%couple_windspeed_time = input_time
1578 0 : self%L1_windspeed(:, 1_i4) = windspeed(:)
1579 : end if
1580 :
1581 : ! ssrd
1582 0 : if (present(ssrd)) then
1583 0 : if (.not. self%couple_ssrd) &
1584 0 : call error_message("meteo_handler%set_meteo: ssrd was not set to be coupled.")
1585 0 : self%couple_ssrd_time = input_time
1586 0 : self%L1_ssrd(:, 1_i4) = ssrd(:)
1587 : end if
1588 :
1589 : ! strd
1590 0 : if (present(strd)) then
1591 0 : if (.not. self%couple_strd) &
1592 0 : call error_message("meteo_handler%set_meteo: strd was not set to be coupled.")
1593 0 : self%couple_strd_time = input_time
1594 0 : self%L1_strd(:, 1_i4) = strd(:)
1595 : end if
1596 :
1597 : ! tann
1598 0 : if (present(tann)) then
1599 0 : if (.not. self%couple_tann) &
1600 0 : call error_message("meteo_handler%set_meteo: tann was not set to be coupled.")
1601 0 : self%couple_tann_time = input_time
1602 0 : self%L1_tann(:, 1_i4) = tann(:)
1603 : end if
1604 :
1605 13104 : end subroutine set_meteo
1606 :
1607 0 : end module mo_meteo_handler
|
