Line data Source code
1 : !> \file mo_common_restart.f90
2 : !> \brief \copybrief mo_common_read_data
3 : !> \details \copydetails mo_common_read_data
4 :
5 : !> \brief common restart tools
6 : !> \details Routines to deal with grid infos for restart files
7 : !> \authors Robert Schweppe
8 : !> \date Jun 2018
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_restart
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 :: write_grid_info
22 : PUBLIC :: read_grid_info ! read restart files for configuration from a given path
23 : PUBLIC :: read_nLAI_and_check_dims
24 :
25 :
26 : !> \brief check consistency of two given items
27 : INTERFACE check_consistency_element
28 : MODULE PROCEDURE check_consistency_element_i4, check_consistency_element_dp
29 : end interface check_consistency_element
30 :
31 :
32 : CONTAINS
33 :
34 :
35 : !> \brief write restart files for each domain
36 : !> \details write restart files for each domain. For each domain
37 : !! three restart files are written. These are xxx_states.nc,
38 : !! xxx_L11_config.nc, and xxx_config.nc (xxx being the three digit
39 : !! domain index). If a variable is added here, it should also be added
40 : !! in the read restart routines below.
41 : !> \changelog
42 : !! - Stephan Thober Aug 2015
43 : !! - moved write of routing states to mRM
44 : !! - David Schaefer Nov 2015
45 : !! - mo_netcdf
46 : !! - Stephan Thober Nov 2016
47 : !! - moved processMatrix to common variables
48 : !! - Zink M. Demirel C. Mar 2017
49 : !! - Added Jarvis soil water stress function at SM process(3)
50 : !! - Robert Schweppe Feb 2018
51 : !! - Removed all L0 references
52 : !! - Robert Schweppe Jun 2018
53 : !! - refactoring and reformatting
54 : !! - Stephan Thober May 2019
55 : !! - where statement for gnu73 to translate level0 mask
56 : !> \authors Stephan Thober
57 : !> \date Jun 2014
58 92 : subroutine write_grid_info(grid_in, level_name, nc)
59 :
60 : use mo_common_constants, only : nodata_dp, nodata_i4
61 : use mo_common_types, only: Grid
62 : use mo_netcdf, only : NcDataset, NcDimension, NcVariable
63 :
64 : implicit none
65 :
66 : !> level to be written
67 : type(Grid), intent(in) :: grid_in
68 :
69 : !> level_id
70 : character(*), intent(in) :: level_name
71 :
72 : !> NcDataset to write information to
73 : type(NcDataset), intent(inout) :: nc
74 :
75 : ! dummy for gnu73
76 92 : integer(i4), allocatable :: dummy(:, :)
77 :
78 : type(NcDimension) :: rows, cols
79 :
80 : type(NcVariable) :: var
81 :
82 :
83 92 : rows = nc%setDimension("nrows" // trim(level_name), grid_in%nrows)
84 92 : cols = nc%setDimension("ncols" // trim(level_name), grid_in%ncols)
85 :
86 : ! now set everything related to the grid
87 276 : var = nc%setVariable("L" // trim(level_name) // "_domain_mask", "i32", (/rows, cols/))
88 92 : call var%setFillValue(nodata_i4)
89 : ! transform from logical to i32
90 : ! ST: where statement is used because gnu73 does not properly translate with merge
91 368 : allocate(dummy(size(grid_in%mask, 1), size(grid_in%mask, 2)))
92 4231207 : dummy = 0_i4
93 4231115 : where(grid_in%mask) dummy = 1_i4
94 92 : call var%setData(dummy)
95 92 : deallocate(dummy)
96 92 : call var%setAttribute("long_name", "Mask at level " // trim(level_name))
97 :
98 276 : var = nc%setVariable("L" // trim(level_name) // "_domain_lat", "f64", (/rows, cols/))
99 92 : call var%setFillValue(nodata_dp)
100 92 : call var%setData(grid_in%y)
101 92 : call var%setAttribute("long_name", "Latitude at level " // trim(level_name))
102 :
103 276 : var = nc%setVariable("L" // trim(level_name) // "_domain_lon", "f64", (/rows, cols/))
104 92 : call var%setFillValue(nodata_dp)
105 92 : call var%setData(grid_in%x)
106 92 : call var%setAttribute("long_name", "Longitude at level " // trim(level_name))
107 :
108 276 : var = nc%setVariable("L" // trim(level_name) // "_domain_cellarea", "f64", (/rows, cols/))
109 92 : call var%setFillValue(nodata_dp)
110 1575971 : call var%setData(unpack(grid_in%CellArea * 1.0E-6_dp, grid_in%mask, nodata_dp))
111 92 : call var%setAttribute("long_name", "Cell area at level " // trim(level_name))
112 :
113 92 : call nc%setAttribute("xllcorner_L" // trim(level_name), grid_in%xllcorner)
114 92 : call nc%setAttribute("yllcorner_L" // trim(level_name), grid_in%yllcorner)
115 92 : call nc%setAttribute("cellsize_L" // trim(level_name), grid_in%cellsize)
116 92 : call nc%setAttribute("nrows_L" // trim(level_name), grid_in%nrows)
117 92 : call nc%setAttribute("ncols_L" // trim(level_name), grid_in%ncols)
118 92 : call nc%setAttribute("nCells_L" // trim(level_name), grid_in%nCells)
119 :
120 92 : end subroutine write_grid_info
121 :
122 :
123 : !> \brief reads configuration apart from Level 11 configuration from a restart directory
124 : !> \details read configuration variables from a given restart
125 : !> directory and initializes all configuration variables,
126 : !> that are initialized in the subroutine initialise,
127 : !> contained in module mo_startup.
128 : !> \changelog
129 : !! - David Schaefer Nov 2015
130 : !! - mo_netcdf
131 : !! - Zink M. Demirel C. Mar 2017
132 : !! - Added Jarvis soil water stress function at SM process(3)
133 : !! - Robert Schweppe Feb 2018
134 : !! - Removed all L0 references
135 : !! - Robert Schweppe Jun 2018
136 : !! - refactoring and reformatting
137 : !! - Stephan Thober May 2019
138 : !! - added allocation check for mask and cellArea because cellArea needs to be read by mRM, but mask is created before by mHM
139 : !> \authors Stephan Thober
140 : !> \date Apr 2013
141 4 : subroutine read_grid_info(InFile, level_name, new_grid)
142 :
143 92 : use mo_common_types, only: Grid
144 : use mo_netcdf, only : NcDataset, NcVariable
145 :
146 : implicit none
147 :
148 : !> Input Path including trailing slash
149 : character(256), intent(in) :: InFile
150 :
151 : !> level_name (id)
152 : character(*), intent(in) :: level_name
153 :
154 : !> grid to save information to
155 : type(Grid), intent(inout) :: new_grid
156 :
157 : ! dummy, 2 dimension I4
158 4 : integer(i4), dimension(:, :), allocatable :: dummyI2
159 :
160 : ! dummy, 2 dimension DP
161 4 : real(dp), dimension(:, :), allocatable :: dummyD2
162 :
163 : character(256) :: Fname
164 :
165 : type(NcDataset) :: nc
166 :
167 : type(NcVariable) :: var
168 :
169 : integer(i4) :: k
170 :
171 :
172 : ! read config
173 4 : fname = trim(InFile)
174 4 : call message(' Reading config from ', trim(adjustl(Fname)), ' ...')
175 :
176 4 : nc = NcDataset(fname, "r")
177 :
178 : ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
179 : ! Read L1 variables <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
180 : ! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
181 : ! read the grid properties
182 4 : call nc%getAttribute("xllcorner_L" // trim(level_name), new_grid%xllcorner)
183 4 : call nc%getAttribute("yllcorner_L" // trim(level_name), new_grid%yllcorner)
184 4 : call nc%getAttribute("nrows_L" // trim(level_name), new_grid%nrows)
185 4 : call nc%getAttribute("ncols_L" // trim(level_name), new_grid%ncols)
186 4 : call nc%getAttribute("cellsize_L" // trim(level_name), new_grid%cellsize)
187 4 : call nc%getAttribute("nCells_L" // trim(level_name), new_grid%nCells)
188 :
189 13 : if (.not. allocated(new_grid%mask)) allocate(new_grid%mask(new_grid%nrows, new_grid%ncols))
190 13 : if (.not. allocated(new_grid%x)) allocate(new_grid%x(new_grid%nrows, new_grid%ncols))
191 13 : if (.not. allocated(new_grid%y)) allocate(new_grid%y(new_grid%nrows, new_grid%ncols))
192 : ! read L1 mask
193 4 : var = nc%getVariable("L" // trim(level_name) // "_domain_mask")
194 : ! read integer
195 4 : call var%getData(dummyI2)
196 : ! transform to logical
197 249830 : new_grid%mask = (dummyI2 .eq. 1_i4)
198 :
199 4 : var = nc%getVariable("L" // trim(level_name) // "_domain_lat")
200 4 : call var%getData(new_grid%y)
201 :
202 4 : var = nc%getVariable("L" // trim(level_name) // "_domain_lon")
203 4 : call var%getData(new_grid%x)
204 :
205 4 : var = nc%getVariable("L" // trim(level_name) // "_domain_cellarea")
206 4 : call var%getData(dummyD2)
207 171594 : if (.not. allocated(new_grid%CellArea)) new_grid%CellArea = pack(dummyD2 / 1.0E-6_dp, new_grid%mask)
208 : ! new_grid%CellArea = pack(dummyD2 / 1.0E-6_dp, new_grid%mask)
209 :
210 4 : call nc%close()
211 :
212 279482 : new_grid%Id = (/ (k, k = 1, new_grid%nCells) /)
213 :
214 4 : end subroutine read_grid_info
215 :
216 :
217 : !> \brief read nubmer of LAI time steps and check dimension configurations read from restart file
218 : !> \author Robert Schweppe
219 : !> \date Aug 2019
220 : !> \author Sebastian Mueller
221 : !> \date Feb 2023
222 1 : subroutine read_nLAI_and_check_dims(iDomain, InFile)
223 :
224 4 : use mo_mpr_global_variables, only: nLAI, LAIBoundaries ! may read from restart
225 : use mo_netcdf, only : NcDataset, NcVariable, NcDimension
226 : use mo_common_constants, only : soilHorizonsVarName, landCoverPeriodsVarName, LAIVarName
227 : use mo_common_mHM_mRM_variables, only: read_old_style_restart_bounds
228 :
229 : implicit none
230 :
231 : !> domain counter (not ID)
232 : integer(i4), intent(in) :: iDomain
233 : !> Input Path including trailing slash
234 : character(256), intent(in) :: InFile
235 :
236 : character(256) :: fname
237 : type(NcDataset) :: nc
238 : type(NcVariable) :: var
239 : type(NcDimension) :: nc_dim
240 :
241 : integer(i4) :: nSoilHorizons_temp, nLAIs_temp, nLandCoverPeriods_temp
242 1 : real(dp), dimension(:), allocatable :: landCoverPeriodBoundaries_temp, soilHorizonBoundaries_temp, &
243 1 : LAIBoundaries_temp
244 :
245 : ! dummy, 2 dimension
246 1 : real(dp), dimension(:, :), allocatable :: dummyD2, dummyD2_tmp
247 :
248 : integer(i4) :: ii
249 :
250 : ! read config
251 1 : fname = trim(InFile)
252 1 : call message(' Reading and checking LAI, land-cover and soil-horizons from ', trim(adjustl(Fname)), ' ...')
253 :
254 1 : nc = NcDataset(fname, "r")
255 :
256 : ! get the dimensions
257 1 : var = nc%getVariable(trim(soilHorizonsVarName)//'_bnds')
258 1 : call var%getData(dummyD2_tmp)
259 : if (allocated(dummyD2)) deallocate(dummyD2)
260 1 : if ( read_old_style_restart_bounds ) then
261 0 : allocate(dummyD2(size(dummyD2_tmp,2), size(dummyD2_tmp,1)))
262 0 : dummyD2 = transpose(dummyD2_tmp)
263 : else
264 4 : allocate(dummyD2(size(dummyD2_tmp,1), size(dummyD2_tmp,2)))
265 8 : dummyD2 = dummyD2_tmp
266 : end if
267 1 : deallocate(dummyD2_tmp)
268 1 : nSoilHorizons_temp = size(dummyD2, 2)
269 3 : allocate(soilHorizonBoundaries_temp(nSoilHorizons_temp+1))
270 3 : soilHorizonBoundaries_temp(1:nSoilHorizons_temp) = dummyD2(1, :)
271 1 : soilHorizonBoundaries_temp(nSoilHorizons_temp+1) = dummyD2(2, nSoilHorizons_temp)
272 :
273 : ! get the landcover dimension
274 1 : var = nc%getVariable(trim(landCoverPeriodsVarName)//'_bnds')
275 1 : call var%getData(dummyD2_tmp)
276 1 : if (allocated(dummyD2)) deallocate(dummyD2)
277 1 : if ( read_old_style_restart_bounds ) then
278 0 : allocate(dummyD2(size(dummyD2_tmp,2), size(dummyD2_tmp,1)))
279 0 : dummyD2 = transpose(dummyD2_tmp)
280 : else
281 4 : allocate(dummyD2(size(dummyD2_tmp,1), size(dummyD2_tmp,2)))
282 8 : dummyD2 = dummyD2_tmp
283 : end if
284 1 : deallocate(dummyD2_tmp)
285 1 : nLandCoverPeriods_temp = size(dummyD2, 2)
286 3 : allocate(landCoverPeriodBoundaries_temp(nLandCoverPeriods_temp+1))
287 3 : landCoverPeriodBoundaries_temp(1:nLandCoverPeriods_temp) = dummyD2(1, :)
288 1 : landCoverPeriodBoundaries_temp(nLandCoverPeriods_temp+1) = dummyD2(2, nLandCoverPeriods_temp)
289 :
290 : ! get the LAI dimension
291 1 : if (nc%hasVariable(trim(LAIVarName)//'_bnds')) then
292 1 : var = nc%getVariable(trim(LAIVarName)//'_bnds')
293 1 : call var%getData(dummyD2_tmp)
294 1 : if (allocated(dummyD2)) deallocate(dummyD2)
295 1 : if ( read_old_style_restart_bounds ) then
296 0 : allocate(dummyD2(size(dummyD2_tmp,2), size(dummyD2_tmp,1)))
297 0 : dummyD2 = transpose(dummyD2_tmp)
298 : else
299 4 : allocate(dummyD2(size(dummyD2_tmp,1), size(dummyD2_tmp,2)))
300 38 : dummyD2 = dummyD2_tmp
301 : end if
302 1 : deallocate(dummyD2_tmp)
303 1 : nLAIs_temp = size(dummyD2, 2)
304 3 : allocate(LAIBoundaries_temp(nLAIs_temp+1))
305 13 : LAIBoundaries_temp(1:nLAIs_temp) = dummyD2(1, :)
306 1 : LAIBoundaries_temp(nLAIs_temp+1) = dummyD2(2, nLAIs_temp)
307 0 : else if (nc%hasDimension('L1_LAITimesteps')) then
308 0 : nc_dim = nc%getDimension('L1_LAITimesteps')
309 0 : nLAIs_temp = nc_dim%getLength()
310 0 : allocate(LAIBoundaries_temp(nLAIs_temp+1))
311 0 : LAIBoundaries_temp = [(ii, ii=1, nLAIs_temp+1)]
312 : else
313 0 : call error_message('***ERROR: no LAI information in restart file for reading')
314 : end if
315 :
316 : ! check LAI for consistency on all domains (-1 indicates first reading)
317 1 : if (nLAI == -1_i4) then
318 1 : nLAI = nLAIs_temp
319 3 : allocate(LAIBoundaries(nLAI + 1_i4))
320 15 : LAIBoundaries = LAIBoundaries_temp
321 : end if
322 :
323 : call check_dimension_consistency(iDomain, nSoilHorizons_temp, soilHorizonBoundaries_temp, &
324 1 : nLAIs_temp, LAIBoundaries_temp, nLandCoverPeriods_temp, landCoverPeriodBoundaries_temp)
325 :
326 1 : call nc%close()
327 :
328 1 : end subroutine read_nLAI_and_check_dims
329 :
330 : !> \brief checks dimension configurations read from restart file
331 : !> \authors Robert Schweppe
332 : !> \date Aug 2019
333 1 : subroutine check_dimension_consistency(iDomain, nSoilHorizons_temp, soilHorizonBoundaries_temp, &
334 1 : nLAIs_temp, LAIBoundaries_temp, nLandCoverPeriods_temp, landCoverPeriodBoundaries_temp)
335 :
336 1 : use mo_mpr_global_variables, only: nSoilHorizons_mHM, HorizonDepth_mHM, nLAI, LAIBoundaries ! may read from restart
337 : use mo_common_variables, only: nLCoverScene, LC_year_start, LC_year_end ! read from nml
338 : use mo_string_utils, only: compress, num2str
339 :
340 : integer(i4), intent(in) :: iDomain
341 :
342 : integer(i4), intent(in) :: nSoilHorizons_temp, nLAIs_temp, nLandCoverPeriods_temp
343 : real(dp), dimension(:), intent(inout) :: landCoverPeriodBoundaries_temp, soilHorizonBoundaries_temp, &
344 : LAIBoundaries_temp
345 : character(256) :: errorString
346 :
347 : integer(i4) :: k
348 :
349 : ! compare local to global
350 1 : call check_consistency_element(nLCoverScene, nLandCoverPeriods_temp, 'number of land cover periods', iDomain)
351 1 : call check_consistency_element(nSoilHorizons_mHM, nSoilHorizons_temp, 'number of soil horizons', iDomain)
352 1 : call check_consistency_element(nLAI, nLAIs_temp, 'number of LAI timesteps', iDomain)
353 :
354 : ! now check the boundaries
355 3 : do k=1, nLCoverScene
356 2 : errorString = compress(trim(num2str(k)))//'th land cover boundary'
357 3 : call check_consistency_element(real(LC_year_start(k), dp), landCoverPeriodBoundaries_temp(k), errorString, iDomain)
358 : end do
359 1 : errorString = 'last land cover boundary (with 1 year added due to real/int conversion) '
360 0 : call check_consistency_element(real(LC_year_end(nLCoverScene) + 1_i4, dp), &
361 1 : landCoverPeriodBoundaries_temp(nLCoverScene+1), errorString, iDomain)
362 :
363 : ! last soil horizon is spatially variable, so this is not checked yet
364 : ! first soil horizon 0 and not contained in HorizonDepth_mHM, so skip that, too
365 2 : do k=2, nSoilHorizons_mHM
366 1 : errorString = compress(trim(num2str(k)))//'th soil horizon boundary'
367 2 : call check_consistency_element(HorizonDepth_mHM(k-1), soilHorizonBoundaries_temp(k), errorString, iDomain)
368 : end do
369 :
370 14 : do k=1, nLAI+1
371 13 : errorString = compress(trim(num2str(k)))//'th LAI period boundary'
372 14 : call check_consistency_element(LAIBoundaries(k), LAIBoundaries_temp(k), errorString, iDomain)
373 : end do
374 :
375 1 : end subroutine check_dimension_consistency
376 :
377 : !> \copydoc check_consistency_element
378 17 : subroutine check_consistency_element_dp(item1, item2, name, iDomain)
379 1 : use mo_utils, only: ne
380 : use mo_string_utils, only: compress, num2str
381 :
382 : real(dp), intent(in) :: item1, item2
383 : character(*), intent(in) :: name
384 : integer(i4), intent(in) :: iDomain
385 :
386 17 : if (ne(item1, item2)) then
387 : call error_message( &
388 : 'The ', trim(name),&
389 : ' as set in the configuration file (', &
390 : compress(trim(num2str(item1))), &
391 : ') does not conform with domain ', &
392 0 : compress(trim(num2str(iDomain))), ' (', compress(trim(num2str(item2))), ').')
393 : end if
394 17 : end subroutine check_consistency_element_dp
395 :
396 : !> \copydoc check_consistency_element
397 3 : subroutine check_consistency_element_i4(item1, item2, name, iDomain)
398 17 : use mo_string_utils, only: compress, num2str
399 :
400 : integer(i4), intent(in) :: item1, item2
401 : character(*), intent(in) :: name
402 : integer(i4), intent(in) :: iDomain
403 :
404 3 : if (item1 /= item2) then
405 : call error_message( &
406 : 'The ', trim(name),&
407 : ' as set in the configuration file (', &
408 : compress(trim(num2str(item1))), &
409 : ') does not conform with domain ', &
410 0 : compress(trim(num2str(iDomain))), ' (', compress(trim(num2str(item2))), ').')
411 : end if
412 3 : end subroutine check_consistency_element_i4
413 :
414 : end module mo_common_restart
|
