Line data Source code
1 : !> \file mo_grid.f90
2 : !> \brief \copybrief mo_grid
3 : !> \details \copydetails mo_grid
4 :
5 : !> \brief gridding tools
6 : !> \details Common tools to deal with grids in mHM.
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_grid
13 : use mo_kind, only : dp, i4
14 :
15 : IMPLICIT NONE
16 :
17 : PRIVATE
18 :
19 : PUBLIC :: init_lowres_level, set_domain_indices, L0_grid_setup, &
20 : mapCoordinates, geoCoordinates
21 : contains
22 : ! ------------------------------------------------------------------
23 :
24 : ! NAME
25 : ! init_lowres_level
26 :
27 : ! PURPOSE
28 : !> \brief Level-1 variable initialization
29 :
30 : !> \details following tasks are performed for L1 datasets
31 : !> - cell id & numbering
32 : !> - mask creation
33 : !> - storage of cell cordinates (row and coloum id)
34 : !> - sorage of four corner L0 cordinates
35 : !> If a variable is added or removed here, then it also has to
36 : !> be added or removed in the subroutine config_variables_set in
37 : !> module mo_restart and in the subroutine set_config in module
38 : !> mo_set_netcdf_restart
39 :
40 : ! INTENT(IN)
41 : !> \param[in] "type(Grid) :: highres"
42 : !> \param[in] "real(dp) :: target_resolution"
43 :
44 : ! INTENT(INOUT)
45 : !> \param[inout] "type(Grid) :: lowres"
46 :
47 : ! INTENT(INOUT), OPTIONAL
48 : !> \param[inout] "type(GridRemapper), optional :: highres_lowres_remap"
49 :
50 : ! HISTORY
51 : !> \authors Rohini Kumar
52 :
53 : !> \date Jan 2013
54 :
55 : ! Modifications:
56 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
57 :
58 75 : subroutine init_lowres_level(highres, target_resolution, lowres, highres_lowres_remap)
59 :
60 : use mo_common_constants, only : nodata_dp, nodata_i4
61 : use mo_common_types, only : Grid, GridRemapper
62 :
63 : implicit none
64 :
65 : type(Grid), target, intent(in) :: highres
66 :
67 : real(dp), intent(in) :: target_resolution
68 :
69 : type(Grid), target, intent(inout) :: lowres
70 :
71 : type(GridRemapper), intent(inout), optional :: highres_lowres_remap
72 :
73 75 : real(dp), dimension(:, :), allocatable :: areaCell0_2D
74 :
75 75 : real(dp) :: cellFactor
76 :
77 : integer(i4) :: iup, idown
78 :
79 : integer(i4) :: jl, jr
80 :
81 : integer(i4) :: i, j, k, ic, jc
82 :
83 : ! STEPS ::
84 :
85 : !--------------------------------------------------------
86 : ! 1) Estimate each variable locally for a given domain
87 : ! 2) Pad each variable to its corresponding global one
88 : !--------------------------------------------------------
89 : ! grid properties
90 75 : if (.not. allocated(lowres%mask)) then
91 : call calculate_grid_properties(highres%nrows, highres%ncols, &
92 : highres%xllcorner, highres%yllcorner, highres%cellsize, &
93 : target_resolution, &
94 : lowres%nrows, lowres%ncols, &
95 75 : lowres%xllcorner, lowres%yllcorner, lowres%cellsize)
96 : ! cellfactor = leve1-1 / level-0
97 75 : cellFactor = anint(lowres%cellsize / highres%cellsize)
98 :
99 : ! allocation and initalization of mask at level-1
100 300 : allocate(lowres%mask(lowres%nrows, lowres%ncols))
101 5457 : lowres%mask(:, :) = .FALSE.
102 :
103 : ! create mask at level-1
104 31323 : do j = 1, highres%ncols
105 31248 : jc = ceiling(real(j, dp) / cellFactor)
106 8961627 : do i = 1, highres%nrows
107 8930304 : if (.NOT. highres%mask(i, j)) cycle
108 3337755 : ic = ceiling(real(i, dp) / cellFactor)
109 8961552 : lowres%mask(ic, jc) = .TRUE.
110 : end do
111 : end do
112 :
113 : ! estimate ncells and initalize related variables
114 5457 : lowres%nCells = count(lowres%mask)
115 : ! allocate and initalize cell1 related variables
116 225 : allocate(lowres%Id (lowres%nCells))
117 8655 : lowres%Id = (/ (k, k = 1, lowres%nCells) /)
118 : end if
119 :
120 75 : if (present(highres_lowres_remap)) then
121 : ! cellfactor = leve1-1 / level-0, set again in case not yet initialized
122 49 : cellFactor = lowres%cellsize / highres%cellsize
123 :
124 : ! lowres additional properties
125 196 : allocate(areaCell0_2D(highres%nrows, highres%ncols))
126 49 : areaCell0_2D(:, :) = UNPACK(highres%CellArea, highres%mask, nodata_dp)
127 :
128 49 : if (.not. allocated(lowres%CellCoor)) then
129 147 : allocate(lowres%CellCoor (lowres%nCells, 2))
130 147 : allocate(lowres%CellArea (lowres%nCells))
131 : end if
132 :
133 147 : allocate(highres_lowres_remap%lower_bound(lowres%nCells))
134 147 : allocate(highres_lowres_remap%upper_bound(lowres%nCells))
135 147 : allocate(highres_lowres_remap%left_bound (lowres%nCells))
136 147 : allocate(highres_lowres_remap%right_bound(lowres%nCells))
137 147 : allocate(highres_lowres_remap%n_subcells (lowres%nCells))
138 196 : allocate(highres_lowres_remap%lowres_id_on_highres (highres%nrows, highres%ncols))
139 5849569 : highres_lowres_remap%lowres_id_on_highres = nodata_i4
140 :
141 49 : highres_lowres_remap%high_res_grid => highres
142 49 : highres_lowres_remap%low_res_grid => lowres
143 :
144 49 : k = 0
145 519 : do jc = 1, lowres%ncols
146 3859 : do ic = 1, lowres%nrows
147 3340 : if (.NOT. lowres%mask(ic, jc)) cycle
148 1981 : k = k + 1
149 :
150 1981 : lowres%CellCoor(k, 1) = ic
151 1981 : lowres%CellCoor(k, 2) = jc
152 :
153 : ! coord. of all corners -> of finer scale level-0
154 1981 : iup = (ic - 1) * nint(cellFactor, i4) + 1
155 1981 : idown = ic * nint(cellFactor, i4)
156 1981 : jl = (jc - 1) * nint(cellFactor, i4) + 1
157 1981 : jr = jc * nint(cellFactor, i4)
158 :
159 : ! constrain the range of up, down, left, and right boundaries
160 1981 : if(iup < 1) iup = 1
161 1981 : if(idown > highres%nrows) idown = highres%nrows
162 1981 : if(jl < 1) jl = 1
163 1981 : if(jr > highres%ncols) jr = highres%ncols
164 :
165 1981 : highres_lowres_remap%upper_bound (k) = iup
166 1981 : highres_lowres_remap%lower_bound (k) = idown
167 1981 : highres_lowres_remap%left_bound (k) = jl
168 1981 : highres_lowres_remap%right_bound(k) = jr
169 :
170 : ! effective area [km2] & total no. of L0 cells within a given L1 cell
171 3706453 : lowres%CellArea(k) = sum(areacell0_2D(iup : idown, jl : jr), highres%mask(iup : idown, jl : jr))
172 3706453 : highres_lowres_remap%n_subcells(k) = count(highres%mask(iup : idown, jl : jr))
173 : ! Delimitation of level-11 cells on level-0
174 3706923 : highres_lowres_remap%lowres_id_on_highres(iup : idown, jl : jr) = k
175 : end do
176 : end do
177 :
178 : ! free space
179 49 : deallocate(areaCell0_2D)
180 :
181 : end if
182 :
183 75 : end subroutine init_lowres_level
184 :
185 : ! NAME
186 : ! set_domain_indices
187 :
188 : ! PURPOSE
189 : !> \brief TODO: add description
190 :
191 : !> \details TODO: add description
192 :
193 : ! INTENT(INOUT)
194 : !> \param[inout] "type(Grid), dimension(:) :: grids"
195 :
196 : ! HISTORY
197 : !> \authors Robert Schweppe
198 :
199 : !> \date Jun 2018
200 :
201 : ! Modifications:
202 : ! Stephan Thober, Aug 2019 - added optional indices for L0 data because L0 data can be shared among domains
203 :
204 162 : subroutine set_domain_indices(grids, indices)
205 :
206 75 : use mo_common_types, only: Grid
207 :
208 : implicit none
209 :
210 : type(Grid), intent(inout), dimension(:) :: grids
211 : integer(i4), intent(in), dimension(:), optional :: indices
212 :
213 : integer(i4) :: iDomain
214 :
215 :
216 234 : do iDomain = 1, size(grids)
217 : ! Saving indices of mask and packed data
218 153 : if(iDomain .eq. 1_i4) then
219 81 : grids(iDomain)%iStart = 1_i4
220 : else
221 72 : if (present(indices)) then
222 24 : grids(iDomain)%iStart = grids(indices(iDomain - 1))%iEnd + 1_i4
223 : else
224 48 : grids(iDomain)%iStart = grids(iDomain - 1_i4)%iEnd + 1_i4
225 : end if
226 : end if
227 234 : grids(iDomain)%iEnd = grids(iDomain)%iStart + grids(iDomain)%nCells - 1_i4
228 : end do
229 :
230 81 : end subroutine set_domain_indices
231 :
232 : ! ------------------------------------------------------------------
233 :
234 : ! NAME
235 : ! L0_grid_setup
236 :
237 : ! PURPOSE
238 : !> \brief level 0 variable initialization
239 :
240 : !> \details following tasks are performed for L0 data sets
241 : !> - cell id & numbering
242 : !> - storage of cell cordinates (row and coloum id)
243 : !> - empirical dist. of terrain slope
244 : !> - flag to determine the presence of a particular soil id
245 : !> in this configuration of the model run
246 : !> If a variable is added or removed here, then it also has to
247 : !> be added or removed in the subroutine config_variables_set in
248 : !> module mo_restart and in the subroutine set_config in module
249 : !> mo_set_netcdf_restart
250 :
251 : ! INTENT(INOUT)
252 : !> \param[inout] "type(Grid) :: new_grid"
253 :
254 : ! HISTORY
255 : !> \authors Rohini Kumar
256 :
257 : !> \date Jan 2013
258 :
259 : ! Modifications:
260 : ! Rohini Kumar & Matthias Cuntz May 2014 - cell area calulation based on a regular lat-lon grid or
261 : ! on a regular X-Y coordinate system
262 : ! Matthias Cuntz May 2014 - changed empirical distribution function so that doubles get the same value
263 : ! Matthias Zink & Matthias Cuntz Feb 2016 - code speed up due to reformulation of CDF calculation
264 : ! Rohini Kumar Mar 2016 - changes for handling multiple soil database options
265 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
266 :
267 22 : subroutine L0_grid_setup(new_grid)
268 :
269 81 : use mo_common_types, only: Grid
270 : use mo_common_variables, only : iFlag_cordinate_sys
271 : use mo_constants, only : RadiusEarth_dp, TWOPI_dp
272 :
273 : implicit none
274 :
275 : type(Grid), intent(inout) :: new_grid
276 :
277 22 : real(dp), dimension(:, :), allocatable :: areaCell_2D
278 :
279 : integer(i4) :: i, j, k
280 :
281 22 : real(dp) :: rdum, degree_to_radian, degree_to_metre
282 :
283 : ! STEPS ::
284 :
285 :
286 : !--------------------------------------------------------
287 : ! 1) Estimate each variable locally for a given domain
288 : ! 2) Pad each variable to its corresponding global one
289 : !--------------------------------------------------------
290 : ! level-0 information
291 2612662 : new_grid%nCells = count(new_grid%mask)
292 :
293 66 : allocate(new_grid%CellCoor(new_grid%nCells, 2))
294 66 : allocate(new_grid%Id(new_grid%nCells))
295 66 : allocate(new_grid%CellArea(new_grid%nCells))
296 88 : allocate(areaCell_2D(new_grid%nrows, new_grid%ncols))
297 :
298 2918894 : new_grid%Id = (/ (k, k = 1, new_grid%nCells) /)
299 :
300 : !------------------------------------------------
301 : ! start looping for cell cordinates and ids
302 : !------------------------------------------------
303 22 : k = 0
304 9142 : do j = 1, new_grid%ncols
305 2612662 : do i = 1, new_grid%nrows
306 2603520 : if (.NOT. new_grid%mask(i, j)) cycle
307 972950 : k = k + 1
308 972950 : new_grid%cellCoor(k, 1) = i
309 2612640 : new_grid%cellCoor(k, 2) = j
310 : end do
311 : end do
312 :
313 : ! ESTIMATE AREA [m2]
314 :
315 : ! regular X-Y coordinate system
316 22 : if(iFlag_cordinate_sys .eq. 0) then
317 972972 : new_grid%CellArea(:) = new_grid%cellsize * new_grid%cellsize
318 :
319 : ! regular lat-lon coordinate system
320 0 : else if(iFlag_cordinate_sys .eq. 1) then
321 :
322 0 : degree_to_radian = TWOPI_dp / 360.0_dp
323 0 : degree_to_metre = RadiusEarth_dp * TWOPI_dp / 360.0_dp
324 0 : do i = new_grid%ncols, 1, -1
325 0 : j = new_grid%ncols - i + 1
326 : ! get latitude in degrees
327 0 : rdum = new_grid%yllcorner + (real(j, dp) - 0.5_dp) * new_grid%cellsize
328 : ! convert to radians
329 0 : rdum = rdum * degree_to_radian
330 : ! AREA [m2]
331 0 : areaCell_2D(:, i) = (new_grid%cellsize * cos(rdum) * degree_to_metre) * (new_grid%cellsize * degree_to_metre)
332 : end do
333 0 : new_grid%CellArea(:) = pack(areaCell_2D(:, :), new_grid%mask)
334 :
335 : end if
336 :
337 : ! free space
338 22 : deallocate(areaCell_2D)
339 :
340 22 : end subroutine L0_grid_setup
341 :
342 :
343 : !------------------------------------------------------------------
344 : ! NAME
345 : ! mapCoordinates
346 :
347 : ! PURPOSE
348 : !> \brief Generate map coordinates
349 :
350 : !> \details Generate map coordinate arrays for given domain and level
351 :
352 : ! INTENT(IN)
353 : !> \param[in] "type(Grid) :: level" -> grid reference
354 :
355 : ! INTENT(OUT)
356 : !> \param[out] "real(dp), dimension(:) :: x, y"
357 : !> \param[out] "real(dp), dimension(:) :: x, y"
358 :
359 : ! HISTORY
360 : !> \authors Matthias Zink
361 :
362 : !> \date Apr 2013
363 :
364 : ! Modifications:
365 : ! Stephan Thober Nov 2013 - removed fproj dependency
366 : ! David Schaefer Jun 2015 - refactored the former subroutine CoordSystem
367 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
368 :
369 :
370 19 : subroutine mapCoordinates(level, y, x)
371 :
372 22 : use mo_common_types, only: Grid
373 :
374 : implicit none
375 :
376 : ! -> grid reference
377 : type(Grid), intent(in) :: level
378 :
379 : real(dp), intent(out), allocatable, dimension(:) :: x, y
380 :
381 : integer(i4) :: ii, ncols, nrows
382 :
383 19 : real(dp) :: cellsize
384 :
385 :
386 19 : cellsize = level%cellsize
387 19 : nrows = level%nrows
388 19 : ncols = level%ncols
389 :
390 95 : allocate(x(nrows), y(ncols))
391 :
392 19 : x(1) = level%xllcorner + 0.5_dp * cellsize
393 124 : do ii = 2, nrows
394 124 : x(ii) = x(ii - 1) + cellsize
395 : end do
396 :
397 : ! inverse for Panoply, ncview display
398 19 : y(ncols) = level%yllcorner + 0.5_dp * cellsize
399 181 : do ii = ncols - 1, 1, -1
400 181 : y(ii) = y(ii + 1) + cellsize
401 : end do
402 :
403 19 : end subroutine mapCoordinates
404 :
405 : !------------------------------------------------------------------
406 : ! NAME
407 : ! geoCoordinates
408 :
409 : ! PURPOSE
410 : !> \brief Generate geographic coordinates
411 :
412 : !> \details Generate geographic coordinate arrays for given domain and level
413 :
414 : ! INTENT(IN)
415 : !> \param[in] "type(Grid) :: level" -> grid reference
416 :
417 : ! INTENT(OUT)
418 : !> \param[out] "real(dp), dimension(:, :) :: lat, lon"
419 : !> \param[out] "real(dp), dimension(:, :) :: lat, lon"
420 :
421 : ! HISTORY
422 : !> \authors Matthias Zink
423 :
424 : !> \date Apr 2013
425 :
426 : ! Modifications:
427 : ! Stephan Thober Nov 2013 - removed fproj dependency
428 : ! David Schaefer Jun 2015 - refactored the former subroutine CoordSystem
429 : ! Stephan Thober Sep 2015 - using mask to unpack coordinates
430 : ! Stephan Thober Oct 2015 - writing full lat/lon again
431 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
432 :
433 19 : subroutine geoCoordinates(level, lat, lon)
434 :
435 19 : use mo_common_types, only: Grid
436 :
437 : implicit none
438 :
439 : ! -> grid reference
440 : type(Grid), intent(in) :: level
441 :
442 : real(dp), intent(out), allocatable, dimension(:, :) :: lat, lon
443 :
444 :
445 1549 : lat = level%y
446 1549 : lon = level%x
447 :
448 19 : end subroutine geoCoordinates
449 :
450 : ! ------------------------------------------------------------------
451 :
452 : ! NAME
453 : ! calculate_grid_properties
454 :
455 : ! PURPOSE
456 : !> \brief Calculates basic grid properties at a required coarser level using
457 : !> information of a given finer level.
458 : !> Calculates basic grid properties at a required coarser level (e.g., L11) using
459 : !> information of a given finer level (e.g., L0). Basic grid properties such as
460 : !> nrows, ncols, xllcorner, yllcorner cellsize are estimated in this
461 : !> routine.
462 :
463 : !> \details TODO: add description
464 :
465 : ! INTENT(IN)
466 : !> \param[in] "integer(i4) :: nrowsIn" no. of rows at an input level
467 : !> \param[in] "integer(i4) :: ncolsIn" no. of cols at an input level
468 : !> \param[in] "real(dp) :: xllcornerIn" xllcorner at an input level
469 : !> \param[in] "real(dp) :: yllcornerIn" yllcorner at an input level
470 : !> \param[in] "real(dp) :: cellsizeIn" cell size at an input level
471 : !> \param[in] "real(dp) :: aimingResolution" resolution of an output level
472 :
473 : ! INTENT(OUT)
474 : !> \param[out] "integer(i4) :: nrowsOut" no. of rows at an output level
475 : !> \param[out] "integer(i4) :: ncolsOut" no. of cols at an output level
476 : !> \param[out] "real(dp) :: xllcornerOut" xllcorner at an output level
477 : !> \param[out] "real(dp) :: yllcornerOut" yllcorner at an output level
478 : !> \param[out] "real(dp) :: cellsizeOut" cell size at an output level
479 :
480 : ! HISTORY
481 : !> \authors Matthias Zink & Rohini Kumar
482 :
483 : !> \date Feb 2013
484 :
485 : ! Modifications:
486 : ! R. Kumar Sep 2013 - documentation added according to the template
487 : ! Robert Schweppe Jun 2018 - refactoring and reformatting
488 :
489 75 : subroutine calculate_grid_properties(nrowsIn, ncolsIn, xllcornerIn, yllcornerIn, cellsizeIn, aimingResolution, &
490 : nrowsOut, ncolsOut, xllcornerOut, yllcornerOut, cellsizeOut)
491 :
492 19 : use mo_message, only : error_message
493 : use mo_string_utils, only : num2str
494 :
495 : implicit none
496 :
497 : ! no. of rows at an input level
498 : integer(i4), intent(in) :: nrowsIn
499 :
500 : ! no. of cols at an input level
501 : integer(i4), intent(in) :: ncolsIn
502 :
503 : ! xllcorner at an input level
504 : real(dp), intent(in) :: xllcornerIn
505 :
506 : ! yllcorner at an input level
507 : real(dp), intent(in) :: yllcornerIn
508 :
509 : ! cell size at an input level
510 : real(dp), intent(in) :: cellsizeIn
511 :
512 : ! resolution of an output level
513 : real(dp), intent(in) :: aimingResolution
514 :
515 : ! no. of rows at an output level
516 : integer(i4), intent(out) :: nrowsOut
517 :
518 : ! no. of cols at an output level
519 : integer(i4), intent(out) :: ncolsOut
520 :
521 : ! xllcorner at an output level
522 : real(dp), intent(out) :: xllcornerOut
523 :
524 : ! yllcorner at an output level
525 : real(dp), intent(out) :: yllcornerOut
526 :
527 : ! cell size at an output level
528 : real(dp), intent(out) :: cellsizeOut
529 :
530 75 : real(dp) :: cellFactor, rounded
531 : integer(i4) :: rounded_int
532 :
533 :
534 75 : cellFactor = aimingResolution / cellsizeIn
535 75 : rounded = anint(cellFactor)
536 75 : rounded_int = nint(cellFactor)
537 :
538 75 : if (abs(rounded - cellFactor) > 1.e-7_dp) then
539 : call error_message( &
540 : '***ERROR: Two resolutions size do not confirm: ', &
541 : trim(adjustl(num2str(nint(AimingResolution)))), &
542 0 : trim(adjustl(num2str(nint(cellsizeIn)))))
543 : end if
544 :
545 75 : cellsizeOut = aimingResolution
546 75 : ncolsOut = nint(real(ncolsIn, dp) / cellFactor)
547 75 : nrowsOut = nint(real(nrowsIn, dp) / cellFactor)
548 :
549 : ! if we rounded down, but now we would miss cells, add rows and/or cols
550 75 : if ( ncolsOut * rounded_int < ncolsIn ) ncolsOut = ncolsOut + 1_i4
551 75 : if ( nrowsOut * rounded_int < nrowsIn ) nrowsOut = nrowsOut + 1_i4
552 :
553 75 : xllcornerOut = xllcornerIn + real(ncolsIn, dp) * aimingResolution / rounded - real(ncolsOut, dp) * cellsizeOut
554 75 : yllcornerOut = yllcornerIn + real(nrowsIn, dp) * aimingResolution / rounded - real(nrowsOut, dp) * cellsizeOut
555 :
556 75 : end subroutine calculate_grid_properties
557 :
558 : end module mo_grid
|
