FMS  2025.04
Flexible Modeling System
grid2.F90
1 !***********************************************************************
2 !* Apache License 2.0
3 !*
4 !* This file is part of the GFDL Flexible Modeling System (FMS).
5 !*
6 !* Licensed under the Apache License, Version 2.0 (the "License");
7 !* you may not use this file except in compliance with the License.
8 !* You may obtain a copy of the License at
9 !*
10 !* http://www.apache.org/licenses/LICENSE-2.0
11 !*
12 !* FMS is distributed in the hope that it will be useful, but WITHOUT
13 !* WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied;
14 !* without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
15 !* PARTICULAR PURPOSE. See the License for the specific language
16 !* governing permissions and limitations under the License.
17 !***********************************************************************
18 !> @defgroup grid2_mod grid2_mod
19 !> @ingroup mosaic2
20 !> @brief Routines for grid calculations, using @ref fms2_io
21 
22 !> @addtogroup grid2_mod
23 !> @{
24 module grid2_mod
25 
26 use mpp_mod, only : mpp_root_pe, mpp_error, uppercase, lowercase, fatal, note
27 use constants_mod, only : pi, radius
28 use fms2_io_mod, only : get_global_attribute, read_data, global_att_exists, &
29  variable_exists, file_exists, open_file, close_file, get_variable_size, &
30  fmsnetcdffile_t
31 use fms_string_utils_mod, only: string
34 
35 ! the following two use statement are only needed for define_cube_mosaic
36 use mpp_domains_mod, only : domain2d, mpp_define_mosaic, mpp_get_compute_domain, &
38 use mosaic2_mod, only : get_mosaic_ncontacts, get_mosaic_contact
39 use platform_mod
40 
41 implicit none;private
42 
43 ! ==== public interfaces =====================================================
44 ! grid dimension inquiry subroutines
45 public :: get_great_circle_algorithm ! returns great_circle_algorithm
46 public :: get_grid_ntiles ! returns number of tiles
47 public :: get_grid_size ! returns horizontal sizes of the grid
48 ! grid geometry inquiry subroutines
49 public :: get_grid_cell_centers
50 public :: get_grid_cell_vertices
51 ! grid area inquiry subroutines
52 public :: get_grid_cell_area
53 public :: get_grid_comp_area
54 ! decompose cubed sphere domains -- probably does not belong here, but it should
55 ! be in some place available for component models
56 public :: define_cube_mosaic
57 public :: grid_init
58 public :: grid_end
59 ! ==== end of public interfaces ==============================================
60 
61 !> Gets the size of the grid for one or all tiles
62 !> @ingroup grid2_mod
63 interface get_grid_size
64  module procedure get_grid_size_for_all_tiles
65  module procedure get_grid_size_for_one_tile
66 end interface
67 
68 !> Gets arrays of global grid cell boundaries for given model component and
69 !! mosaic tile number
70 !> @ingroup grid2_mod
72  module procedure get_grid_cell_vertices_1d_r4
73  module procedure get_grid_cell_vertices_1d_r8
74  module procedure get_grid_cell_vertices_2d_r4
75  module procedure get_grid_cell_vertices_2d_r8
76  module procedure get_grid_cell_vertices_ug_r4
77  module procedure get_grid_cell_vertices_ug_r8
78 end interface
79 
80 !> Gets grid cell centers
81 !> @ingroup grid2_mod
83  module procedure get_grid_cell_centers_1d_r4
84  module procedure get_grid_cell_centers_1d_r8
85  module procedure get_grid_cell_centers_2d_r4
86  module procedure get_grid_cell_centers_2d_r8
87  module procedure get_grid_cell_centers_ug_r4
88  module procedure get_grid_cell_centers_ug_r8
89 end interface
90 
91 !> Finds area of a grid cell
92 !> @ingroup grid2_mod
94  module procedure get_grid_cell_area_sg_r4
95  module procedure get_grid_cell_area_sg_r8
96  module procedure get_grid_cell_area_ug_r4
97  module procedure get_grid_cell_area_ug_r8
98 end interface get_grid_cell_area
99 
100 !> Gets the area of a given component per grid cell
101 !> @ingroup grid2_mod
103  module procedure get_grid_comp_area_sg_r4
104  module procedure get_grid_comp_area_sg_r8
105  module procedure get_grid_comp_area_ug_r4
106  module procedure get_grid_comp_area_ug_r8
107 end interface get_grid_comp_area
108 
109 !> @addtogroup grid2_mod
110 !> @{
111 ! ==== module constants ======================================================
112 character(len=*), parameter :: &
113  module_name = 'grid2_mod'
114 
115 ! Include variable "version" to be written to log file.
116 #include<file_version.h>
117 
118 character(len=*), parameter :: &
119  grid_dir = 'INPUT/', & !< root directory for all grid files
120  grid_file = 'INPUT/grid_spec.nc' !< name of the grid spec file
121 
122 integer, parameter :: &
123  max_name = 256, & !< max length of the variable names
124  version_geolon_t = 0, & !< indicates gelon_t variable is present in grid_file
125  version_x_t = 1, & !< indicates x_t variable is present in grid_file
126  version_ocn_mosaic_file = 2, & !< indicates ocn_mosaic_file variable is present in grid_file
127  version_gridfiles = 3 !< indicates gridfiles variable is present in grid_file
128 
129 integer, parameter :: bufsize = 1048576 !< This is used to control memory usage in get_grid_comp_area
130  !! We may change this to a namelist variable is needed.
131 
132 ! ==== module variables ======================================================
133 integer :: grid_version = -1 !< Value to indicate what type of grid file is being read,
134  !! based on which variables are present
135 logical :: great_circle_algorithm = .false.
136 logical :: module_is_initialized = .false.
137 logical :: grid_spec_exists = .true.
138 type(fmsnetcdffile_t) :: gridfileobj
139 type(fmsnetcdffile_t), dimension(3) :: mosaic_fileobj
140 
141 contains
142 
143 !> @brief Initialize the grid2 module
144 subroutine grid_init
145  if (module_is_initialized) return
146  if (.not. file_exists(grid_file)) then
147  module_is_initialized = .true.
148  grid_spec_exists = .false.
149  return
150  endif
151  call open_grid_file(gridfileobj, grid_file)
152  great_circle_algorithm = get_great_circle_algorithm()
153  grid_version = get_grid_version(gridfileobj)
156  module_is_initialized = .true.
157 end subroutine grid_init
158 
159 !> @brief Shutdown the grid2 module
160 subroutine grid_end
161  if (.not. module_is_initialized) return
162  if (grid_spec_exists) then
164  call close_file(gridfileobj)
165  endif
166  module_is_initialized = .false.
167 end subroutine grid_end
168 
169 !> @brief Checks that the grid2 module was initialized propertly
170 subroutine init_checks(subroutine_name)
171  character(len=*), intent(in) :: subroutine_name !< Name of the subroutine calling this from
172 
173  if (.not. module_is_initialized) then
174  call mpp_error(fatal, "grid2_mod::"//trim(subroutine_name)//" is being called but grid2 was never initialized. "//&
175  "Please ensure that grid_init is called before calling "//trim(subroutine_name)//".")
176  endif
177 
178  if (.not. grid_spec_exists) then
179  call mpp_error(fatal, "grid2_mod::"//trim(subroutine_name)//" is being called, but "//trim(grid_file)//&
180  " does not exist, so grid2 was not initialized properly."//&
181  " Please ensure that "//trim(grid_file)//" is accessible.")
182  endif
183 end subroutine init_checks
184 
185 !> @brief Determine if we are using the great circle algorithm
186 !! @return Logical flag describing if we are using the great circlealgorithm
188  character(len=128) :: attvalue
189  logical :: get_great_circle_algorithm
190 
192  if (.not. grid_spec_exists) return
193  if (global_att_exists(gridfileobj, "great_circle_algorithm")) then
194  call get_global_attribute(gridfileobj, "great_circle_algorithm", attvalue)
195  if(trim(attvalue) == "TRUE") then
197  else if(trim(attvalue) .NE. "FALSE") then
198  call mpp_error(fatal, module_name//&
199  '/get_great_circle_algorithm value of global attribute "great_circle_algorthm" in file'// &
200  trim(grid_file)//' should be TRUE or FALSE')
201  endif
202  endif
203 end function get_great_circle_algorithm
204 
205 !> @brief Open a grid file
206 subroutine open_grid_file(myfileobj, myfilename)
207  type(fmsnetcdffile_t), intent(out) :: myfileobj !< File object of grid file
208  character(len=*), intent(in) :: myfilename!< Name of the grid file
209  if(.not. open_file(myfileobj, myfilename, 'read')) then
210  call mpp_error(fatal, 'grid2_mod(open_grid_file):Error in opening file '//trim(myfilename))
211  endif
212 end subroutine open_grid_file
213 
214 !> @brief Open a mosaic file
215 subroutine open_mosaic_file(mymosaicfileobj, component)
216  type(fmsnetcdffile_t), intent(out) :: mymosaicfileobj !< File object returned
217  character(len=3), intent(in) :: component !< Component (atm, lnd, etc.)
218 
219  character(len=FMS_PATH_LEN) :: mosaicfilename
220 
221  call read_data(gridfileobj,trim(lowercase(component))//'_mosaic_file', mosaicfilename)
222  call open_grid_file(mymosaicfileobj, grid_dir//trim(mosaicfilename))
223 end subroutine open_mosaic_file
224 
225 !> @brief Read a tile file name from a netcdf file
226 !! @return Name of the file as a string
227 function read_file_name(thisfileobj, filevar, level)
228  type(fmsnetcdffile_t), intent(in) :: thisfileobj !< File object of file
229  character(len=*), intent(in) :: filevar!< Variable containing file names
230  integer, intent(in) :: level !< Level of tile file
231  integer, dimension(2) :: file_list_size
232  character(len=FMS_PATH_LEN) :: read_file_name
233  character(len=FMS_PATH_LEN), dimension(:), allocatable :: file_names
234 
235  call get_variable_size(thisfileobj, filevar, file_list_size)
236  allocate(file_names(file_list_size(2)))
237  call read_data(thisfileobj, filevar, file_names)
238  read_file_name = file_names(level)
239  deallocate(file_names)
240 end function read_file_name
241 
242 !> @brief Get the grid version from a file object
243 !! @return An integer representation of the grid version
244 function get_grid_version(fileobj)
245  type(fmsnetcdffile_t), intent(in) :: fileobj !< File object of grid file
246  integer :: get_grid_version
247 
248  if(grid_version<0) then
249  if(variable_exists(fileobj, 'geolon_t')) then
251  else if(variable_exists(fileobj, 'x_T')) then
253  else if(variable_exists(fileobj, 'ocn_mosaic_file') ) then
255  else if(variable_exists(fileobj, 'gridfiles') ) then
257  else
258  call mpp_error(fatal, module_name//'/get_grid_version Can''t determine the version of the grid spec:'// &
259  & ' none of "x_T", "geolon_t", or "ocn_mosaic_file" exist in file "'//trim(grid_file)//'"')
260  endif
261  endif
262 end function get_grid_version
263 
264 !> @brief Assign the component mosaic files if grid_spec is Version 3
266  mosaic_fileobj(1) = gridfileobj
267  mosaic_fileobj(2) = gridfileobj
268  mosaic_fileobj(3) = gridfileobj
269 end subroutine assign_component_mosaics
270 
271 !> @brief Open the component mosaic files for atm, lnd, and ocn
273  if (variable_exists(gridfileobj, 'atm_mosaic_file')) call open_mosaic_file(mosaic_fileobj(1), 'atm')
274  if (variable_exists(gridfileobj, 'ocn_mosaic_file')) call open_mosaic_file(mosaic_fileobj(2), 'ocn')
275  if (variable_exists(gridfileobj, 'lnd_mosaic_file')) call open_mosaic_file(mosaic_fileobj(3), 'lnd')
276 end subroutine open_component_mosaics
277 
278 !> @brief Close the component mosaic files for atm, lnd, and ocn
280  if (variable_exists(gridfileobj, 'atm_mosaic_file')) call close_file(mosaic_fileobj(1))
281  if (variable_exists(gridfileobj, 'ocn_mosaic_file')) call close_file(mosaic_fileobj(2))
282  if (variable_exists(gridfileobj, 'lnd_mosaic_file')) call close_file(mosaic_fileobj(3))
283 end subroutine close_component_mosaics
284 
285 !> @brief Get the component number of a model component (atm, lnd, ocn)
286 !! @return Integer component number
287 function get_component_number(component)
288  character(len=*), intent(in) :: component !< Component model (atm, lnd, ocn)
289  integer :: get_component_number
290  select case(lowercase(component))
291  case('atm')
293  case('ocn')
295  case('lnd')
297  end select
298 end function get_component_number
299 
300 !> @brief returns number of tiles for a given component
301 subroutine get_grid_ntiles(component,ntiles)
302  character(len=*) :: component !< Component model (atm, lnd, ocn)
303  integer, intent(out) :: ntiles !< Number of tiles
304 
305  call init_checks("get_grid_ntiles")
306  select case (grid_version)
308  ntiles = 1
310  ntiles = get_mosaic_ntiles(mosaic_fileobj(get_component_number(trim(component))))
311  end select
312 end subroutine get_grid_ntiles
313 
314 !> @brief returns size of the grid for each of the tiles
315 subroutine get_grid_size_for_all_tiles(component,nx,ny)
316  character(len=*) :: component !< Component model (atm, lnd, ocn)
317  integer, intent(inout) :: nx(:),ny(:) !< Grid size in x and y
318 
319  ! local vars
320  integer :: siz(2) ! for the size of external fields
321  character(len=MAX_NAME) :: varname1
322 
323  call init_checks("get_grid_size")
324  varname1 = 'AREA_'//trim(uppercase(component))
325 
326  select case (grid_version)
328  call get_variable_size(gridfileobj, varname1, siz)
329  nx(1) = siz(1); ny(1)=siz(2)
330  case(version_ocn_mosaic_file, version_gridfiles) ! mosaic file
331  call get_mosaic_grid_sizes(mosaic_fileobj(get_component_number(trim(component))),nx,ny)
332  end select
333 end subroutine get_grid_size_for_all_tiles
334 
335 !> @brief returns size of the grid for one of the tiles
336 subroutine get_grid_size_for_one_tile(component,tile,nx,ny)
337  character(len=*) :: component !< Component model (atm, lnd, ocn)
338  integer, intent(in) :: tile !< Tile number
339  integer, intent(inout) :: nx,ny !< Grid size in x and y
340 
341  ! local vars
342  integer, allocatable :: nnx(:), nny(:)
343  integer :: ntiles
344 
345  call init_checks("get_grid_size")
346  call get_grid_ntiles(component, ntiles)
347  if(tile>0.and.tile<=ntiles) then
348  allocate(nnx(ntiles),nny(ntiles))
349  call get_grid_size_for_all_tiles(component,nnx,nny)
350  nx = nnx(tile); ny = nny(tile)
351  deallocate(nnx,nny)
352  else
353  call mpp_error(fatal, 'get_grid_size'//&
354  'requested tile index '//trim(string(tile))//' is out of bounds (1:'//trim(string(ntiles))//')')
355  endif
356 end subroutine get_grid_size_for_one_tile
357 
358 !> @brief given a model component, a layout, and (optionally) a halo size, returns a
359 !! domain for current processor
360 subroutine define_cube_mosaic(component, domain, layout, halo, maskmap)
361  character(len=*) , intent(in) :: component !< Component model (atm, lnd, ocn)
362  type(domain2d) , intent(inout) :: domain !< Domain
363  integer , intent(in) :: layout(2) !< Layout
364  integer, optional, intent(in) :: halo !< Halo
365  logical, optional, intent(in) :: maskmap(:,:,:) !< Maskmap
366 
367  ! ---- local vars
368  integer :: ntiles ! number of tiles
369  integer :: ncontacts ! number of contacts between mosaic tiles
370  integer :: n
371  integer :: ng, pe_pos, npes ! halo size
372  integer, allocatable :: nlon(:), nlat(:), global_indices(:,:)
373  integer, allocatable :: pe_start(:), pe_end(:), layout_2d(:,:)
374  integer, allocatable :: tile1(:),tile2(:)
375  integer, allocatable :: is1(:),ie1(:),js1(:),je1(:)
376  integer, allocatable :: is2(:),ie2(:),js2(:),je2(:)
377 
378  call init_checks("define_cube_mosaic")
379  call get_grid_ntiles(component,ntiles)
380  allocate(nlon(ntiles), nlat(ntiles))
381  allocate(global_indices(4,ntiles))
382  allocate(pe_start(ntiles),pe_end(ntiles))
383  allocate(layout_2d(2,ntiles))
384  call get_grid_size(component,nlon,nlat)
385 
386  pe_pos = mpp_root_pe()
387  do n = 1, ntiles
388  global_indices(:,n) = (/ 1, nlon(n), 1, nlat(n) /)
389  layout_2d(:,n) = layout
390  if(present(maskmap)) then
391  npes = count(maskmap(:,:,n))
392  else
393  npes = layout(1)*layout(2)
394  endif
395  pe_start(n) = pe_pos
396  pe_end(n) = pe_pos + npes - 1
397  pe_pos = pe_end(n) + 1
398  enddo
399 
400  ! get the contact information from mosaic file
401  ncontacts = get_mosaic_ncontacts(mosaic_fileobj(get_component_number(trim(component))))
402  allocate(tile1(ncontacts),tile2(ncontacts))
403  allocate(is1(ncontacts),ie1(ncontacts),js1(ncontacts),je1(ncontacts))
404  allocate(is2(ncontacts),ie2(ncontacts),js2(ncontacts),je2(ncontacts))
405  call get_mosaic_contact(mosaic_fileobj(get_component_number(trim(component))), tile1, tile2, &
406  is1, ie1, js1, je1, is2, ie2, js2, je2)
407 
408  ng = 0
409  if(present(halo)) ng = halo
410  ! create the domain2d variable
411  call mpp_define_mosaic ( global_indices, layout_2d, domain, &
412  ntiles, ncontacts, tile1, tile2, &
413  is1, ie1, js1, je1, &
414  is2, ie2, js2, je2, &
415  pe_start=pe_start, pe_end=pe_end, symmetry=.true., &
416  shalo = ng, nhalo = ng, whalo = ng, ehalo = ng, &
417  maskmap = maskmap, &
418  name = trim(component)//'Cubic-Sphere Grid' )
419 
420  deallocate(nlon,nlat,global_indices,pe_start,pe_end,layout_2d)
421  deallocate(tile1,tile2)
422  deallocate(is1,ie1,js1,je1)
423  deallocate(is2,ie2,js2,je2)
424 end subroutine define_cube_mosaic
425 
426 #include "grid2_r4.fh"
427 #include "grid2_r8.fh"
428 
429 end module grid2_mod
430 !> @}
431 ! close documentation grouping
Close a netcdf or domain file opened with open_file or open_virtual_file.
Definition: fms2_io.F90:165
Opens a given netcdf or domain file.
Definition: fms2_io.F90:121
Read data from a defined field in a file.
Definition: fms2_io.F90:291
character(:) function, allocatable, public string(v, fmt)
Converts a number or a Boolean value to a string.
integer grid_version
Value to indicate what type of grid file is being read, based on which variables are present.
Definition: grid2.F90:133
integer, parameter version_gridfiles
indicates gridfiles variable is present in grid_file
Definition: grid2.F90:122
subroutine, public grid_end
Shutdown the grid2 module.
Definition: grid2.F90:161
subroutine, public grid_init
Initialize the grid2 module.
Definition: grid2.F90:145
subroutine open_mosaic_file(mymosaicfileobj, component)
Open a mosaic file.
Definition: grid2.F90:216
integer function get_component_number(component)
Get the component number of a model component (atm, lnd, ocn)
Definition: grid2.F90:288
integer, parameter version_ocn_mosaic_file
indicates ocn_mosaic_file variable is present in grid_file
Definition: grid2.F90:122
integer, parameter version_geolon_t
indicates gelon_t variable is present in grid_file
Definition: grid2.F90:122
character(len=fms_path_len) function read_file_name(thisfileobj, filevar, level)
Read a tile file name from a netcdf file.
Definition: grid2.F90:228
subroutine open_component_mosaics
Open the component mosaic files for atm, lnd, and ocn.
Definition: grid2.F90:273
subroutine, public get_grid_ntiles(component, ntiles)
returns number of tiles for a given component
Definition: grid2.F90:302
subroutine assign_component_mosaics
Assign the component mosaic files if grid_spec is Version 3.
Definition: grid2.F90:266
character(len= *), parameter grid_dir
root directory for all grid files
Definition: grid2.F90:118
integer, parameter version_x_t
indicates x_t variable is present in grid_file
Definition: grid2.F90:122
subroutine open_grid_file(myfileobj, myfilename)
Open a grid file.
Definition: grid2.F90:207
subroutine close_component_mosaics
Close the component mosaic files for atm, lnd, and ocn.
Definition: grid2.F90:280
integer, parameter bufsize
This is used to control memory usage in get_grid_comp_area We may change this to a namelist variable ...
Definition: grid2.F90:129
subroutine init_checks(subroutine_name)
Checks that the grid2 module was initialized propertly.
Definition: grid2.F90:171
subroutine get_grid_size_for_one_tile(component, tile, nx, ny)
returns size of the grid for one of the tiles
Definition: grid2.F90:337
integer, parameter max_name
max length of the variable names
Definition: grid2.F90:122
character(len= *), parameter grid_file
name of the grid spec file
Definition: grid2.F90:118
integer function get_grid_version(fileobj)
Get the grid version from a file object.
Definition: grid2.F90:245
subroutine, public define_cube_mosaic(component, domain, layout, halo, maskmap)
given a model component, a layout, and (optionally) a halo size, returns a domain for current process...
Definition: grid2.F90:361
logical function, public get_great_circle_algorithm()
Determine if we are using the great circle algorithm.
Definition: grid2.F90:188
subroutine get_grid_size_for_all_tiles(component, nx, ny)
returns size of the grid for each of the tiles
Definition: grid2.F90:316
Finds area of a grid cell.
Definition: grid2.F90:93
Gets grid cell centers.
Definition: grid2.F90:82
Gets arrays of global grid cell boundaries for given model component and mosaic tile number.
Definition: grid2.F90:71
Gets the area of a given component per grid cell.
Definition: grid2.F90:102
Gets the size of the grid for one or all tiles.
Definition: grid2.F90:63
integer function, public get_mosaic_xgrid_size(fileobj)
return exchange grid size of mosaic xgrid file.
Definition: mosaic2.F90:118
subroutine, public get_mosaic_contact(fileobj, tile1, tile2, istart1, iend1, jstart1, jend1, istart2, iend2, jstart2, jend2)
Get contact information from mosaic_file Example usage: call get_mosaic_contact(mosaic_file,...
Definition: mosaic2.F90:222
integer function, public get_mosaic_ntiles(fileobj)
Get number of tiles in the mosaic_file.
Definition: mosaic2.F90:134
integer function, public get_mosaic_ncontacts(fileobj)
Get number of contacts in the mosaic_file.
Definition: mosaic2.F90:150
subroutine, public get_mosaic_grid_sizes(fileobj, nx, ny)
Get grid size of each tile from mosaic_file.
Definition: mosaic2.F90:172
subroutine mpp_define_mosaic(global_indices, layout, domain, num_tile, num_contact, tile1, tile2, istart1, iend1, jstart1, jend1, istart2, iend2, jstart2, jend2, pe_start, pe_end, pelist, whalo, ehalo, shalo, nhalo, xextent, yextent, maskmap, name, memory_size, symmetry, xflags, yflags, tile_id)
Defines a domain for mosaic tile grids.
These routines retrieve the axis specifications associated with the compute domains....
These routines retrieve the axis specifications associated with the global domains....
Passes data from a structured grid to an unstructured grid Example usage:
The domain2D type contains all the necessary information to define the global, compute and data domai...
Domain information for managing data on unstructured grids.
Error handler.
Definition: mpp.F90:381