MOM_grid.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Provides the ocean grid type
6module mom_grid
7
9use mom_domains, only : mom_domain_type, get_domain_extent, compute_block_extent
10use mom_domains, only : get_global_shape, deallocate_mom_domain
11use mom_error_handler, only : mom_error, mom_mesg, fatal
12use mom_file_parser, only : get_param, log_param, log_version, param_file_type
14
15implicit none ; private
16
17#include <MOM_memory.h>
18
21
22! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
23! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
24! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
25! vary with the Boussinesq approximation, the Boussinesq variant is given first.
26
27!> Ocean grid type. See mom_grid for details.
28type, public :: ocean_grid_type
29 type(mom_domain_type), pointer :: domain => null() !< Ocean model domain
30 type(mom_domain_type), pointer :: domain_aux => null() !< A non-symmetric auxiliary domain type.
31 type(hor_index_type) :: hi !< Horizontal index ranges
32 type(hor_index_type), allocatable :: hid(:) !< Horizontal index ranges for downsampling
33
34 integer :: isc !< The start i-index of cell centers within the computational domain
35 integer :: iec !< The end i-index of cell centers within the computational domain
36 integer :: jsc !< The start j-index of cell centers within the computational domain
37 integer :: jec !< The end j-index of cell centers within the computational domain
38
39 integer :: isd !< The start i-index of cell centers within the data domain
40 integer :: ied !< The end i-index of cell centers within the data domain
41 integer :: jsd !< The start j-index of cell centers within the data domain
42 integer :: jed !< The end j-index of cell centers within the data domain
43
44 integer :: isg !< The start i-index of cell centers within the global domain
45 integer :: ieg !< The end i-index of cell centers within the global domain
46 integer :: jsg !< The start j-index of cell centers within the global domain
47 integer :: jeg !< The end j-index of cell centers within the global domain
48
49 integer :: iscb !< The start i-index of cell vertices within the computational domain
50 integer :: iecb !< The end i-index of cell vertices within the computational domain
51 integer :: jscb !< The start j-index of cell vertices within the computational domain
52 integer :: jecb !< The end j-index of cell vertices within the computational domain
53
54 integer :: isdb !< The start i-index of cell vertices within the data domain
55 integer :: iedb !< The end i-index of cell vertices within the data domain
56 integer :: jsdb !< The start j-index of cell vertices within the data domain
57 integer :: jedb !< The end j-index of cell vertices within the data domain
58
59 integer :: isgb !< The start i-index of cell vertices within the global domain
60 integer :: iegb !< The end i-index of cell vertices within the global domain
61 integer :: jsgb !< The start j-index of cell vertices within the global domain
62 integer :: jegb !< The end j-index of cell vertices within the global domain
63
64 integer :: isd_global !< The value of isd in the global index space (decomposition invariant).
65 integer :: jsd_global !< The value of isd in the global index space (decomposition invariant).
66 integer :: idg_offset !< The offset between the corresponding global and local i-indices.
67 integer :: jdg_offset !< The offset between the corresponding global and local j-indices.
68 integer :: ke !< The number of layers in the vertical.
69 logical :: symmetric !< True if symmetric memory is used.
70 logical :: nonblocking_updates !< If true, non-blocking halo updates are
71 !! allowed. The default is .false. (for now).
72 integer :: first_direction !< An integer that indicates which direction is
73 !! to be updated first in directionally split
74 !! parts of the calculation. This can be altered
75 !! during the course of the run via calls to
76 !! set_first_direction.
77
78 real allocable_, dimension(NIMEM_,NJMEM_) :: &
79 mask2dt, & !< 0 for land points and 1 for ocean points on the h-grid [nondim].
80 geolatt, & !< The geographic latitude at tracer (h) points [degrees_N] or [km] or [m]
81 geolont, & !< The geographic longitude at tracer (h) points [degrees_E] or [km] or [m]
82 dxt, & !< dxT is delta x at h points [L ~> m].
83 idxt, & !< 1/dxT [L-1 ~> m-1].
84 dyt, & !< dyT is delta y at h points [L ~> m].
85 idyt, & !< IdyT is 1/dyT [L-1 ~> m-1].
86 areat, & !< The area of an h-cell [L2 ~> m2].
87 iareat, & !< 1/areaT [L-2 ~> m-2].
88 sin_rot, & !< The sine of the angular rotation between the local model grid's northward
89 !! and the true northward directions [nondim].
90 cos_rot !< The cosine of the angular rotation between the local model grid's northward
91 !! and the true northward directions [nondim].
92
93 real allocable_, dimension(NIMEMB_PTR_,NJMEM_) :: &
94 mask2dcu, & !< 0 for boundary points and 1 for ocean points on the u grid [nondim].
95 obcmaskcu, & !< 0 for boundary or OBC points and 1 for ocean points on the u grid [nondim].
96 geolatcu, & !< The geographic latitude at u points [degrees_N] or [km] or [m]
97 geoloncu, & !< The geographic longitude at u points [degrees_E] or [km] or [m].
98 dxcu, & !< dxCu is delta x at u points [L ~> m].
99 idxcu, & !< 1/dxCu [L-1 ~> m-1].
100 idxcu_obcmask, & !< 1/dxCu or 0 at boundary or OBC points [L-1 ~> m-1].
101 dycu, & !< dyCu is delta y at u points [L ~> m].
102 idycu, & !< 1/dyCu [L-1 ~> m-1].
103 dy_cu, & !< The unblocked lengths of the u-faces of the h-cell [L ~> m].
104 iareacu, & !< The masked inverse areas of u-grid cells [L-2 ~> m-2].
105 areacu !< The areas of the u-grid cells [L2 ~> m2].
106
107 real allocable_, dimension(NIMEM_,NJMEMB_PTR_) :: &
108 mask2dcv, & !< 0 for boundary points and 1 for ocean points on the v grid [nondim].
109 obcmaskcv, & !< 0 for boundary or OBC points and 1 for ocean points on the v grid [nondim].
110 geolatcv, & !< The geographic latitude at v points [degrees_N] or [km] or [m]
111 geoloncv, & !< The geographic longitude at v points [degrees_E] or [km] or [m].
112 dxcv, & !< dxCv is delta x at v points [L ~> m].
113 idxcv, & !< 1/dxCv [L-1 ~> m-1].
114 dycv, & !< dyCv is delta y at v points [L ~> m].
115 idycv, & !< 1/dyCv [L-1 ~> m-1].
116 idycv_obcmask, & !< 1/dxCv or 0 at boundary or OBC points [L-1 ~> m-1].
117 dx_cv, & !< The unblocked lengths of the v-faces of the h-cell [L ~> m].
118 iareacv, & !< The masked inverse areas of v-grid cells [L-2 ~> m-2].
119 areacv !< The areas of the v-grid cells [L2 ~> m2].
120
121 real allocable_, dimension(NIMEMB_PTR_,NJMEM_) :: &
122 porous_dminu, & !< minimum topographic height (deepest) of U-face [Z ~> m]
123 porous_dmaxu, & !< maximum topographic height (shallowest) of U-face [Z ~> m]
124 porous_davgu !< average topographic height of U-face [Z ~> m]
125
126 real allocable_, dimension(NIMEM_,NJMEMB_PTR_) :: &
127 porous_dminv, & !< minimum topographic height (deepest) of V-face [Z ~> m]
128 porous_dmaxv, & !< maximum topographic height (shallowest) of V-face [Z ~> m]
129 porous_davgv !< average topographic height of V-face [Z ~> m]
130
131 real allocable_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: &
132 mask2dbu, & !< 0 for boundary points and 1 for ocean points on the q grid [nondim].
133 geolatbu, & !< The geographic latitude at q points [degrees_N] or [km] or [m]
134 geolonbu, & !< The geographic longitude at q points [degrees_E] or [km] or [m].
135 dxbu, & !< dxBu is delta x at q points [L ~> m].
136 idxbu, & !< 1/dxBu [L-1 ~> m-1].
137 dybu, & !< dyBu is delta y at q points [L ~> m].
138 idybu, & !< 1/dyBu [L-1 ~> m-1].
139 areabu, & !< areaBu is the area of a q-cell [L2 ~> m2]
140 iareabu !< IareaBu = 1/areaBu [L-2 ~> m-2].
141
142 real, pointer, dimension(:) :: &
143 gridlatt => null(), & !< The latitude of T points for the purpose of labeling the output axes,
144 !! often in units of [degrees_N] or [km] or [m] or [gridpoints].
145 !! On many grids this is the same as geoLatT.
146 gridlatb => null() !< The latitude of B points for the purpose of labeling the output axes,
147 !! often in units of [degrees_N] or [km] or [m] or [gridpoints].
148 !! On many grids this is the same as geoLatBu.
149 real, pointer, dimension(:) :: &
150 gridlont => null(), & !< The longitude of T points for the purpose of labeling the output axes,
151 !! often in units of [degrees_E] or [km] or [m] or [gridpoints].
152 !! On many grids this is the same as geoLonT.
153 gridlonb => null() !< The longitude of B points for the purpose of labeling the output axes,
154 !! often in units of [degrees_E] or [km] or [m] or [gridpoints].
155 !! On many grids this is the same as geoLonBu.
156 character(len=40) :: &
157 ! Except on a Cartesian grid, these are usually some variant of "degrees".
158 x_axis_units, & !< The units that are used in labeling the x coordinate axes.
159 y_axis_units, & !< The units that are used in labeling the y coordinate axes.
160 ! These are internally generated names, including "m", "km", "deg_E" and "deg_N".
161 x_ax_unit_short, & !< A short description of the x-axis units for documenting parameter units
162 y_ax_unit_short !< A short description of the y-axis units for documenting parameter units
163
164 real allocable_, dimension(NIMEM_,NJMEM_) :: &
165 bathyt !< Ocean bottom depth, referenced to Z_ref at tracer points. bathyT is in
166 !! depth units and positive *below* Z_ref [Z ~> m].
167 real allocable_, dimension(NIMEM_,NJMEM_) :: &
168 meansl !< Spatially varying time mean sea level, referenced to Z_ref at tracer points.
169 !! meanSL is in height units and positive *above* Z_ref. It is used
170 !! a) as the height where p = p_atm or zero;
171 !! b) to calculate time mean thickness of the water column, where
172 !! mean thickness = max(meanSL + bathyT, 0.0).
173 !! meanSL is 2D for the consideration of a domain with spatically varying mean
174 !! height, e.g. the Great Lakes system [Z ~> m].
175 real :: z_ref !< A reference value for all geometric height fields, such as bathyT [Z ~> m].
176
177 logical :: bathymetry_at_vel !< If true, there are separate values for the
178 !! basin depths at velocity points. Otherwise the effects of
179 !! of topography are entirely determined from thickness points.
180 real allocable_, dimension(NIMEMB_PTR_,NJMEM_) :: &
181 dblock_u, & !< Topographic depths at u-points at which the flow is blocked [Z ~> m].
182 dopen_u !< Topographic depths at u-points at which the flow is open at width dy_Cu [Z ~> m].
183 real allocable_, dimension(NIMEM_,NJMEMB_PTR_) :: &
184 dblock_v, & !< Topographic depths at v-points at which the flow is blocked [Z ~> m].
185 dopen_v !< Topographic depths at v-points at which the flow is open at width dx_Cv [Z ~> m].
186 real allocable_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: &
187 coriolisbu, & !< The Coriolis parameter at corner points [T-1 ~> s-1].
188 coriolis2bu !< The square of the Coriolis parameter at corner points [T-2 ~> s-2].
189 real allocable_, dimension(NIMEM_,NJMEM_) :: &
190 df_dx, & !< Derivative d/dx f (Coriolis parameter) at h-points [T-1 L-1 ~> s-1 m-1].
191 df_dy !< Derivative d/dy f (Coriolis parameter) at h-points [T-1 L-1 ~> s-1 m-1].
192
193 ! These variables are global sums that are useful for 1-d diagnostics.
194 real :: areat_global !< Global sum of h-cell area [L2 ~> m2]
195 real :: iareat_global !< Global sum of inverse h-cell area (1/areaT_global) [L-2 ~> m-2].
196
197 type(unit_scale_type), pointer :: us => null() !< A dimensional unit scaling type
198
199
200 ! These variables are for block structures.
201 integer :: nblocks !< The number of sub-PE blocks on this PE
202 type(hor_index_type), pointer :: block(:) => null() !< Index ranges for each block
203
204 ! These parameters are run-time parameters that are used during some
205 ! initialization routines (but not all)
206 real :: grid_unit_to_l !< A factor that converts a the geoLat and geoLon variables and related
207 !! variables like len_lat and len_lon into rescaled horizontal distance
208 !! units on a Cartesian grid, in [L km ~> 1000] or [L m-1 ~> 1] or
209 !! is 0 for a non-Cartesian grid.
210 real :: south_lat !< The latitude (or y-coordinate) of the first v-line [degrees_N] or [km] or [m]
211 real :: west_lon !< The longitude (or x-coordinate) of the first u-line [degrees_E] or [km] or [m]
212 real :: len_lat !< The latitudinal (or y-coord) extent of physical domain [degrees_N] or [km] or [m]
213 real :: len_lon !< The longitudinal (or x-coord) extent of physical domain [degrees_E] or [km] or [m]
214 real :: rad_earth_l !< The radius of the planet in rescaled units [L ~> m]
215 real :: max_depth !< The maximum depth of the ocean in depth units [Z ~> m]
216end type ocean_grid_type
217
218contains
219
220!> MOM_grid_init initializes the ocean grid array sizes and grid memory.
221subroutine mom_grid_init(G, param_file, US, HI, global_indexing, bathymetry_at_vel)
222 type(ocean_grid_type), intent(inout) :: g !< The horizontal grid type
223 type(param_file_type), intent(in) :: param_file !< Parameter file handle
224 type(unit_scale_type), optional, pointer :: us !< A dimensional unit scaling type
225 type(hor_index_type), &
226 optional, intent(in) :: hi !< A hor_index_type for array extents
227 logical, optional, intent(in) :: global_indexing !< If true use global index
228 !! values instead of having the data domain on each
229 !! processor start at 1.
230 logical, optional, intent(in) :: bathymetry_at_vel !< If true, there are
231 !! separate values for the ocean bottom depths at
232 !! velocity points. Otherwise the effects of topography
233 !! are entirely determined from thickness points.
234
235 ! Local variables
236 real :: mean_sealev_scale ! A scaling factor for the reference height variable [1] or [Z m-1 ~> 1]
237 integer :: isd, ied, jsd, jed
238 integer :: isdb, iedb, jsdb, jedb
239 integer :: ied_max, jed_max
240 integer :: niblock, njblock, nihalo, njhalo, nblocks, n, i, j
241 logical :: local_indexing ! If false use global index values instead of having
242 ! the data domain on each processor start at 1.
243 ! This include declares and sets the variable "version".
244# include "version_variable.h"
245
246 integer, allocatable, dimension(:) :: ibegin, iend, jbegin, jend
247 character(len=40) :: mod_nm = "MOM_grid" ! This module's name.
248
249 mean_sealev_scale = 1.0 ; if (associated(g%US)) mean_sealev_scale = g%US%m_to_Z
250
251 ! Read all relevant parameters and write them to the model log.
252 call get_param(param_file, mod_nm, "REFERENCE_HEIGHT", g%Z_ref, &
253 units="m", default=0.0, scale=mean_sealev_scale, do_not_log=.true.)
254 call log_version(param_file, mod_nm, version, &
255 "Parameters providing information about the lateral grid.", &
256 log_to_all=.true., layout=.true., all_default=(g%Z_ref==0.0))
257
258 call get_param(param_file, mod_nm, "NIBLOCK", niblock, "The number of blocks "// &
259 "in the x-direction on each processor (for openmp).", default=1, &
260 layoutparam=.true.)
261 call get_param(param_file, mod_nm, "NJBLOCK", njblock, "The number of blocks "// &
262 "in the y-direction on each processor (for openmp).", default=1, &
263 layoutparam=.true.)
264 if (present(us)) then ; if (associated(us)) g%US => us ; endif
265
266 call get_param(param_file, mod_nm, "REFERENCE_HEIGHT", g%Z_ref, &
267 "A reference value for geometric height fields, such as bathyT.", &
268 units="m", default=0.0, scale=mean_sealev_scale)
269
270 if (present(hi)) then
271 g%HI = hi
272
273 g%isc = hi%isc ; g%iec = hi%iec ; g%jsc = hi%jsc ; g%jec = hi%jec
274 g%isd = hi%isd ; g%ied = hi%ied ; g%jsd = hi%jsd ; g%jed = hi%jed
275 g%isg = hi%isg ; g%ieg = hi%ieg ; g%jsg = hi%jsg ; g%jeg = hi%jeg
276
277 g%IscB = hi%IscB ; g%IecB = hi%IecB ; g%JscB = hi%JscB ; g%JecB = hi%JecB
278 g%IsdB = hi%IsdB ; g%IedB = hi%IedB ; g%JsdB = hi%JsdB ; g%JedB = hi%JedB
279 g%IsgB = hi%IsgB ; g%IegB = hi%IegB ; g%JsgB = hi%JsgB ; g%JegB = hi%JegB
280
281 g%idg_offset = hi%idg_offset ; g%jdg_offset = hi%jdg_offset
282 g%isd_global = g%isd + hi%idg_offset ; g%jsd_global = g%jsd + hi%jdg_offset
283 g%symmetric = hi%symmetric
284 else
285 local_indexing = .true.
286 if (present(global_indexing)) local_indexing = .not.global_indexing
287 call hor_index_init(g%Domain, g%HI, param_file, &
288 local_indexing=local_indexing)
289
290 ! get_domain_extent ensures that domains start at 1 for compatibility between
291 ! static and dynamically allocated arrays, unless global_indexing is true.
292 call get_domain_extent(g%Domain, g%isc, g%iec, g%jsc, g%jec, &
293 g%isd, g%ied, g%jsd, g%jed, &
294 g%isg, g%ieg, g%jsg, g%jeg, &
295 g%idg_offset, g%jdg_offset, g%symmetric, &
296 local_indexing=local_indexing)
297 g%isd_global = g%isd+g%idg_offset ; g%jsd_global = g%jsd+g%jdg_offset
298 endif
299
300 g%nonblocking_updates = g%Domain%nonblocking_updates
301
302 ! Set array sizes for fields that are discretized at tracer cell boundaries.
303 g%IscB = g%isc ; g%JscB = g%jsc
304 g%IsdB = g%isd ; g%JsdB = g%jsd
305 g%IsgB = g%isg ; g%JsgB = g%jsg
306 if (g%symmetric) then
307 g%IscB = g%isc-1 ; g%JscB = g%jsc-1
308 g%IsdB = g%isd-1 ; g%JsdB = g%jsd-1
309 g%IsgB = g%isg-1 ; g%JsgB = g%jsg-1
310 endif
311 g%IecB = g%iec ; g%JecB = g%jec
312 g%IedB = g%ied ; g%JedB = g%jed
313 g%IegB = g%ieg ; g%JegB = g%jeg
314
315 call mom_mesg(" MOM_grid.F90, MOM_grid_init: allocating metrics", 5)
316
317 call allocate_metrics(g)
318
319 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
320 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
321
322 g%bathymetry_at_vel = .false.
323 if (present(bathymetry_at_vel)) g%bathymetry_at_vel = bathymetry_at_vel
324 if (g%bathymetry_at_vel) then
325 alloc_(g%Dblock_u(isdb:iedb, jsd:jed)) ; g%Dblock_u(:,:) = -g%Z_ref
326 alloc_(g%Dopen_u(isdb:iedb, jsd:jed)) ; g%Dopen_u(:,:) = -g%Z_ref
327 alloc_(g%Dblock_v(isd:ied, jsdb:jedb)) ; g%Dblock_v(:,:) = -g%Z_ref
328 alloc_(g%Dopen_v(isd:ied, jsdb:jedb)) ; g%Dopen_v(:,:) = -g%Z_ref
329 endif
330
331! setup block indices.
332 nihalo = g%Domain%nihalo
333 njhalo = g%Domain%njhalo
334 nblocks = niblock * njblock
335 if (nblocks < 1) call mom_error(fatal, "MOM_grid_init: " // &
336 "nblocks(=NI_BLOCK*NJ_BLOCK) must be no less than 1")
337
338 allocate(ibegin(niblock), iend(niblock), jbegin(njblock), jend(njblock))
339 call compute_block_extent(g%HI%isc,g%HI%iec,niblock,ibegin,iend)
340 call compute_block_extent(g%HI%jsc,g%HI%jec,njblock,jbegin,jend)
341 !-- make sure the last block is the largest.
342 do i = 1, niblock-1
343 if (iend(i)-ibegin(i) > iend(niblock)-ibegin(niblock) ) call mom_error(fatal, &
344 "MOM_grid_init: the last block size in x-direction is not the largest")
345 enddo
346 do j = 1, njblock-1
347 if (jend(j)-jbegin(j) > jend(njblock)-jbegin(njblock) ) call mom_error(fatal, &
348 "MOM_grid_init: the last block size in y-direction is not the largest")
349 enddo
350
351 g%nblocks = nblocks
352 allocate(g%Block(nblocks))
353 ied_max = 1 ; jed_max = 1
354 do n = 1,nblocks
355 ! Copy all information from the array index type describing the local grid.
356 g%Block(n) = g%HI
357
358 i = mod((n-1), niblock) + 1
359 j = (n-1)/niblock + 1
360 !--- isd and jsd are always 1 for each block to permit array reuse.
361 g%Block(n)%isd = 1 ; g%Block(n)%jsd = 1
362 g%Block(n)%isc = g%Block(n)%isd+nihalo
363 g%Block(n)%jsc = g%Block(n)%jsd+njhalo
364 g%Block(n)%iec = g%Block(n)%isc + iend(i) - ibegin(i)
365 g%Block(n)%jec = g%Block(n)%jsc + jend(j) - jbegin(j)
366 g%Block(n)%ied = g%Block(n)%iec + nihalo
367 g%Block(n)%jed = g%Block(n)%jec + njhalo
368 g%Block(n)%IscB = g%Block(n)%isc ; g%Block(n)%IecB = g%Block(n)%iec
369 g%Block(n)%JscB = g%Block(n)%jsc ; g%Block(n)%JecB = g%Block(n)%jec
370 ! For symmetric memory domains, the first block will have the extra point
371 ! at the lower boundary of its computational domain.
372 if (g%symmetric) then
373 if (i==1) g%Block(n)%IscB = g%Block(n)%IscB-1
374 if (j==1) g%Block(n)%JscB = g%Block(n)%JscB-1
375 endif
376 g%Block(n)%IsdB = g%Block(n)%isd ; g%Block(n)%IedB = g%Block(n)%ied
377 g%Block(n)%JsdB = g%Block(n)%jsd ; g%Block(n)%JedB = g%Block(n)%jed
378 !--- For symmetric memory domain, every block will have an extra point
379 !--- at the lower boundary of its data domain.
380 if (g%symmetric) then
381 g%Block(n)%IsdB = g%Block(n)%IsdB-1
382 g%Block(n)%JsdB = g%Block(n)%JsdB-1
383 endif
384 g%Block(n)%idg_offset = (ibegin(i) - g%Block(n)%isc) + g%HI%idg_offset
385 g%Block(n)%jdg_offset = (jbegin(j) - g%Block(n)%jsc) + g%HI%jdg_offset
386 ! Find the largest values of ied and jed so that all blocks will have the
387 ! same size in memory.
388 ied_max = max(ied_max, g%Block(n)%ied)
389 jed_max = max(jed_max, g%Block(n)%jed)
390 enddo
391
392 ! Reset all of the data domain sizes to match the largest for array reuse,
393 ! recalling that all block have isd=jed=1 for array reuse.
394 do n = 1,nblocks
395 g%Block(n)%ied = ied_max ; g%Block(n)%IedB = ied_max
396 g%Block(n)%jed = jed_max ; g%Block(n)%JedB = jed_max
397 enddo
398
399 !-- do some bounds error checking
400 if ( g%block(nblocks)%ied+g%block(nblocks)%idg_offset > g%HI%ied + g%HI%idg_offset ) &
401 call mom_error(fatal, "MOM_grid_init: G%ied_bk > G%ied")
402 if ( g%block(nblocks)%jed+g%block(nblocks)%jdg_offset > g%HI%jed + g%HI%jdg_offset ) &
403 call mom_error(fatal, "MOM_grid_init: G%jed_bk > G%jed")
404
405end subroutine mom_grid_init
406
407!> set_derived_metrics calculates metric terms that are derived from other metrics.
408subroutine set_derived_metrics(G, US)
409 type(ocean_grid_type), intent(inout) :: g !< The horizontal grid structure
410 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
411! Various inverse grid spacings and derived areas are calculated within this
412! subroutine.
413 integer :: i, j, isd, ied, jsd, jed
414 integer :: isdb, iedb, jsdb, jedb
415
416 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
417 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
418
419 do j=jsd,jed ; do i=isd,ied
420 if (g%dxT(i,j) < 0.0) g%dxT(i,j) = 0.0
421 if (g%dyT(i,j) < 0.0) g%dyT(i,j) = 0.0
422 g%IdxT(i,j) = adcroft_reciprocal(g%dxT(i,j))
423 g%IdyT(i,j) = adcroft_reciprocal(g%dyT(i,j))
424 g%IareaT(i,j) = adcroft_reciprocal(g%areaT(i,j))
425 enddo ; enddo
426
427 do j=jsd,jed ; do i=isdb,iedb
428 if (g%dxCu(i,j) < 0.0) g%dxCu(i,j) = 0.0
429 if (g%dyCu(i,j) < 0.0) g%dyCu(i,j) = 0.0
430 g%IdxCu(i,j) = adcroft_reciprocal(g%dxCu(i,j))
431 g%IdyCu(i,j) = adcroft_reciprocal(g%dyCu(i,j))
432 g%IdxCu_OBCmask(i,j) = g%OBCmaskCu(i,j) * g%IdxCu(i,j) ! This may be reset if masks are reset.
433 enddo ; enddo
434
435 do j=jsdb,jedb ; do i=isd,ied
436 if (g%dxCv(i,j) < 0.0) g%dxCv(i,j) = 0.0
437 if (g%dyCv(i,j) < 0.0) g%dyCv(i,j) = 0.0
438 g%IdxCv(i,j) = adcroft_reciprocal(g%dxCv(i,j))
439 g%IdyCv(i,j) = adcroft_reciprocal(g%dyCv(i,j))
440 g%IdyCv_OBCmask(i,j) = g%OBCmaskCv(i,j) * g%IdyCv(i,j) ! This may be reset if masks are reset.
441 enddo ; enddo
442
443 do j=jsdb,jedb ; do i=isdb,iedb
444 if (g%dxBu(i,j) < 0.0) g%dxBu(i,j) = 0.0
445 if (g%dyBu(i,j) < 0.0) g%dyBu(i,j) = 0.0
446
447 g%IdxBu(i,j) = adcroft_reciprocal(g%dxBu(i,j))
448 g%IdyBu(i,j) = adcroft_reciprocal(g%dyBu(i,j))
449 ! areaBu has usually been set to a positive area elsewhere.
450 if (g%areaBu(i,j) <= 0.0) g%areaBu(i,j) = g%dxBu(i,j) * g%dyBu(i,j)
451 g%IareaBu(i,j) = adcroft_reciprocal(g%areaBu(i,j))
452 enddo ; enddo
453end subroutine set_derived_metrics
454
455!> Adcroft_reciprocal(x) = 1/x for |x|>0 or 0 for x=0.
456function adcroft_reciprocal(val) result(I_val)
457 real, intent(in) :: val !< The value being inverted [A].
458 real :: i_val !< The Adcroft reciprocal of val [A-1].
459
460 i_val = 0.0 ; if (val /= 0.0) i_val = 1.0/val
461end function adcroft_reciprocal
462
463!> Returns true if the coordinates (x,y) are within the h-cell (i,j)
464logical function ispointincell(G, i, j, x, y)
465 type(ocean_grid_type), intent(in) :: g !< Grid type
466 integer, intent(in) :: i !< i index of cell to test
467 integer, intent(in) :: j !< j index of cell to test
468 real, intent(in) :: x !< x coordinate of point [degrees_E]
469 real, intent(in) :: y !< y coordinate of point [degrees_N]
470 ! Local variables
471 real :: xne, xnw, xse, xsw ! Longitudes of cell corners [degrees_E]
472 real :: yne, ynw, yse, ysw ! Latitudes of cell corners [degrees_N]
473 real :: l0, l1, l2, l3 ! Crossed products of differences in position [degrees_E degrees_N]
474 real :: p0, p1, p2, p3 ! Trinary unitary values reflecting the signs of the crossed products [nondim]
475 ispointincell = .false.
476 xne = g%geoLonBu(i ,j ) ; yne = g%geoLatBu(i ,j )
477 xnw = g%geoLonBu(i-1,j ) ; ynw = g%geoLatBu(i-1,j )
478 xse = g%geoLonBu(i ,j-1) ; yse = g%geoLatBu(i ,j-1)
479 xsw = g%geoLonBu(i-1,j-1) ; ysw = g%geoLatBu(i-1,j-1)
480 ! This is a crude calculation that assumes a geographic coordinate system
481 if (x<min(xne,xnw,xse,xsw) .or. x>max(xne,xnw,xse,xsw) .or. &
482 y<min(yne,ynw,yse,ysw) .or. y>max(yne,ynw,yse,ysw) ) then
483 return ! Avoid the more complicated calculation
484 endif
485 l0 = (x-xsw)*(yse-ysw) - (y-ysw)*(xse-xsw)
486 l1 = (x-xse)*(yne-yse) - (y-yse)*(xne-xse)
487 l2 = (x-xne)*(ynw-yne) - (y-yne)*(xnw-xne)
488 l3 = (x-xnw)*(ysw-ynw) - (y-ynw)*(xsw-xnw)
489
490 p0 = sign(1., l0) ; if (l0 == 0.) p0=0.
491 p1 = sign(1., l1) ; if (l1 == 0.) p1=0.
492 p2 = sign(1., l2) ; if (l2 == 0.) p2=0.
493 p3 = sign(1., l3) ; if (l3 == 0.) p3=0.
494
495 if ( (abs(p0)+abs(p2)) + (abs(p1)+abs(p3)) == abs((p0+p2) + (p1+p3)) ) then
496 ispointincell=.true.
497 endif
498end function ispointincell
499
500!> Store an integer indicating which direction to work on first.
501subroutine set_first_direction(G, y_first)
502 type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure
503 integer, intent(in) :: y_first !< The first direction to store
504
505 g%first_direction = y_first
506end subroutine set_first_direction
507
508!> Return global shape of horizontal grid
509subroutine get_global_grid_size(G, niglobal, njglobal)
510 type(ocean_grid_type), intent(inout) :: g !< The horizontal grid type
511 integer, intent(out) :: niglobal !< i-index global size of grid
512 integer, intent(out) :: njglobal !< j-index global size of grid
513
514 call get_global_shape(g%domain, niglobal, njglobal)
515
516end subroutine get_global_grid_size
517
518!> Allocate memory used by the ocean_grid_type and related structures.
519subroutine allocate_metrics(G)
520 type(ocean_grid_type), intent(inout) :: G !< The horizontal grid type
521 integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, isg, ieg, jsg, jeg
522
523 ! This subroutine allocates the lateral elements of the ocean_grid_type that
524 ! are always used and zeros them out.
525
526 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
527 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
528 isg = g%isg ; ieg = g%ieg ; jsg = g%jsg ; jeg = g%jeg
529
530 alloc_(g%dxT(isd:ied,jsd:jed)) ; g%dxT(:,:) = 0.0
531 alloc_(g%dxCu(isdb:iedb,jsd:jed)) ; g%dxCu(:,:) = 0.0
532 alloc_(g%dxCv(isd:ied,jsdb:jedb)) ; g%dxCv(:,:) = 0.0
533 alloc_(g%dxBu(isdb:iedb,jsdb:jedb)) ; g%dxBu(:,:) = 0.0
534 alloc_(g%IdxT(isd:ied,jsd:jed)) ; g%IdxT(:,:) = 0.0
535 alloc_(g%IdxCu(isdb:iedb,jsd:jed)) ; g%IdxCu(:,:) = 0.0
536 alloc_(g%IdxCu_OBCmask(isdb:iedb,jsd:jed)) ; g%IdxCu_OBCmask(:,:) = 0.0
537 alloc_(g%IdxCv(isd:ied,jsdb:jedb)) ; g%IdxCv(:,:) = 0.0
538 alloc_(g%IdxBu(isdb:iedb,jsdb:jedb)) ; g%IdxBu(:,:) = 0.0
539
540 alloc_(g%dyT(isd:ied,jsd:jed)) ; g%dyT(:,:) = 0.0
541 alloc_(g%dyCu(isdb:iedb,jsd:jed)) ; g%dyCu(:,:) = 0.0
542 alloc_(g%dyCv(isd:ied,jsdb:jedb)) ; g%dyCv(:,:) = 0.0
543 alloc_(g%dyBu(isdb:iedb,jsdb:jedb)) ; g%dyBu(:,:) = 0.0
544 alloc_(g%IdyT(isd:ied,jsd:jed)) ; g%IdyT(:,:) = 0.0
545 alloc_(g%IdyCu(isdb:iedb,jsd:jed)) ; g%IdyCu(:,:) = 0.0
546 alloc_(g%IdyCv(isd:ied,jsdb:jedb)) ; g%IdyCv(:,:) = 0.0
547 alloc_(g%IdyCv_OBCmask(isd:ied,jsdb:jedb)) ; g%IdyCv_OBCmask(:,:) = 0.0
548 alloc_(g%IdyBu(isdb:iedb,jsdb:jedb)) ; g%IdyBu(:,:) = 0.0
549
550 alloc_(g%areaT(isd:ied,jsd:jed)) ; g%areaT(:,:) = 0.0
551 alloc_(g%IareaT(isd:ied,jsd:jed)) ; g%IareaT(:,:) = 0.0
552 alloc_(g%areaBu(isdb:iedb,jsdb:jedb)) ; g%areaBu(:,:) = 0.0
553 alloc_(g%IareaBu(isdb:iedb,jsdb:jedb)) ; g%IareaBu(:,:) = 0.0
554
555 alloc_(g%mask2dT(isd:ied,jsd:jed)) ; g%mask2dT(:,:) = 0.0
556 alloc_(g%mask2dCu(isdb:iedb,jsd:jed)) ; g%mask2dCu(:,:) = 0.0
557 alloc_(g%OBCmaskCu(isdb:iedb,jsd:jed)) ; g%OBCmaskCu(:,:) = 0.0
558 alloc_(g%mask2dCv(isd:ied,jsdb:jedb)) ; g%mask2dCv(:,:) = 0.0
559 alloc_(g%OBCmaskCv(isd:ied,jsdb:jedb)) ; g%OBCmaskCv(:,:) = 0.0
560 alloc_(g%mask2dBu(isdb:iedb,jsdb:jedb)) ; g%mask2dBu(:,:) = 0.0
561 alloc_(g%geoLatT(isd:ied,jsd:jed)) ; g%geoLatT(:,:) = 0.0
562 alloc_(g%geoLatCu(isdb:iedb,jsd:jed)) ; g%geoLatCu(:,:) = 0.0
563 alloc_(g%geoLatCv(isd:ied,jsdb:jedb)) ; g%geoLatCv(:,:) = 0.0
564 alloc_(g%geoLatBu(isdb:iedb,jsdb:jedb)) ; g%geoLatBu(:,:) = 0.0
565 alloc_(g%geoLonT(isd:ied,jsd:jed)) ; g%geoLonT(:,:) = 0.0
566 alloc_(g%geoLonCu(isdb:iedb,jsd:jed)) ; g%geoLonCu(:,:) = 0.0
567 alloc_(g%geoLonCv(isd:ied,jsdb:jedb)) ; g%geoLonCv(:,:) = 0.0
568 alloc_(g%geoLonBu(isdb:iedb,jsdb:jedb)) ; g%geoLonBu(:,:) = 0.0
569
570 alloc_(g%dx_Cv(isd:ied,jsdb:jedb)) ; g%dx_Cv(:,:) = 0.0
571 alloc_(g%dy_Cu(isdb:iedb,jsd:jed)) ; g%dy_Cu(:,:) = 0.0
572
573 alloc_(g%porous_DminU(isdb:iedb,jsd:jed)) ; g%porous_DminU(:,:) = 0.0
574 alloc_(g%porous_DmaxU(isdb:iedb,jsd:jed)) ; g%porous_DmaxU(:,:) = 0.0
575 alloc_(g%porous_DavgU(isdb:iedb,jsd:jed)) ; g%porous_DavgU(:,:) = 0.0
576
577 alloc_(g%porous_DminV(isd:ied,jsdb:jedb)) ; g%porous_DminV(:,:) = 0.0
578 alloc_(g%porous_DmaxV(isd:ied,jsdb:jedb)) ; g%porous_DmaxV(:,:) = 0.0
579 alloc_(g%porous_DavgV(isd:ied,jsdb:jedb)) ; g%porous_DavgV(:,:) = 0.0
580
581 alloc_(g%areaCu(isdb:iedb,jsd:jed)) ; g%areaCu(:,:) = 0.0
582 alloc_(g%areaCv(isd:ied,jsdb:jedb)) ; g%areaCv(:,:) = 0.0
583 alloc_(g%IareaCu(isdb:iedb,jsd:jed)) ; g%IareaCu(:,:) = 0.0
584 alloc_(g%IareaCv(isd:ied,jsdb:jedb)) ; g%IareaCv(:,:) = 0.0
585
586 alloc_(g%bathyT(isd:ied, jsd:jed)) ; g%bathyT(:,:) = -g%Z_ref
587 alloc_(g%meanSL(isd:ied, jsd:jed)) ; g%meanSL(:,:) = g%Z_ref
588 alloc_(g%CoriolisBu(isdb:iedb, jsdb:jedb)) ; g%CoriolisBu(:,:) = 0.0
589 alloc_(g%Coriolis2Bu(isdb:iedb, jsdb:jedb)) ; g%Coriolis2Bu(:,:) = 0.0
590 alloc_(g%dF_dx(isd:ied, jsd:jed)) ; g%dF_dx(:,:) = 0.0
591 alloc_(g%dF_dy(isd:ied, jsd:jed)) ; g%dF_dy(:,:) = 0.0
592
593 alloc_(g%sin_rot(isd:ied,jsd:jed)) ; g%sin_rot(:,:) = 0.0
594 alloc_(g%cos_rot(isd:ied,jsd:jed)) ; g%cos_rot(:,:) = 1.0
595
596 allocate(g%gridLonT(isg:ieg), source=0.0)
597 allocate(g%gridLonB(g%IsgB:g%IegB), source=0.0)
598 allocate(g%gridLatT(jsg:jeg), source=0.0)
599 allocate(g%gridLatB(g%JsgB:g%JegB), source=0.0)
600
601end subroutine allocate_metrics
602
603!> Release memory used by the ocean_grid_type and related structures.
604subroutine mom_grid_end(G)
605 type(ocean_grid_type), intent(inout) :: g !< The horizontal grid type
606
607 deallocate(g%Block)
608
609 if (g%bathymetry_at_vel) then
610 dealloc_(g%Dblock_u) ; dealloc_(g%Dopen_u)
611 dealloc_(g%Dblock_v) ; dealloc_(g%Dopen_v)
612 endif
613
614 dealloc_(g%dxT) ; dealloc_(g%dxCu) ; dealloc_(g%dxCv) ; dealloc_(g%dxBu)
615 dealloc_(g%IdxT) ; dealloc_(g%IdxCu) ; dealloc_(g%IdxCv) ; dealloc_(g%IdxBu)
616
617 dealloc_(g%dyT) ; dealloc_(g%dyCu) ; dealloc_(g%dyCv) ; dealloc_(g%dyBu)
618 dealloc_(g%IdyT) ; dealloc_(g%IdyCu) ; dealloc_(g%IdyCv) ; dealloc_(g%IdyBu)
619
620 dealloc_(g%IdxCu_OBCmask) ; dealloc_(g%IdyCv_OBCmask)
621
622 dealloc_(g%areaT) ; dealloc_(g%IareaT)
623 dealloc_(g%areaBu) ; dealloc_(g%IareaBu)
624 dealloc_(g%areaCu) ; dealloc_(g%IareaCu)
625 dealloc_(g%areaCv) ; dealloc_(g%IareaCv)
626
627 dealloc_(g%mask2dT) ; dealloc_(g%mask2dCu) ; dealloc_(g%OBCmaskCu)
628 dealloc_(g%mask2dCv) ; dealloc_(g%OBCmaskCv) ; dealloc_(g%mask2dBu)
629
630 dealloc_(g%geoLatT) ; dealloc_(g%geoLatCu)
631 dealloc_(g%geoLatCv) ; dealloc_(g%geoLatBu)
632 dealloc_(g%geoLonT) ; dealloc_(g%geoLonCu)
633 dealloc_(g%geoLonCv) ; dealloc_(g%geoLonBu)
634
635 dealloc_(g%dx_Cv) ; dealloc_(g%dy_Cu)
636
637 dealloc_(g%bathyT) ; dealloc_(g%meanSL)
638 dealloc_(g%CoriolisBu) ; dealloc_(g%Coriolis2Bu)
639 dealloc_(g%dF_dx) ; dealloc_(g%dF_dy)
640 dealloc_(g%sin_rot) ; dealloc_(g%cos_rot)
641
642 dealloc_(g%porous_DminU) ; dealloc_(g%porous_DmaxU) ; dealloc_(g%porous_DavgU)
643 dealloc_(g%porous_DminV) ; dealloc_(g%porous_DmaxV) ; dealloc_(g%porous_DavgV)
644
645 deallocate(g%gridLonT) ; deallocate(g%gridLatT)
646 deallocate(g%gridLonB) ; deallocate(g%gridLatB)
647
648 ! The cursory flag avoids doing any deallocation of memory in the underlying
649 ! infrastructure to avoid problems due to shared pointers.
650 call deallocate_mom_domain(g%Domain, cursory=.true.)
651
652end subroutine mom_grid_end
653
654!> \namespace mom_grid
655!!
656!! Grid metrics and their inverses are labelled according to their staggered location on a Arakawa C (or B) grid.
657!! - Metrics centered on h- or T-points are labelled T, e.g. dxT is the distance across the cell in the x-direction.
658!! - Metrics centered on u-points are labelled Cu (C-grid u location). e.g. dyCu is the y-distance between
659!! two corners of a T-cell.
660!! - Metrics centered on v-points are labelled Cv (C-grid v location). e.g. dyCv is the y-distance between two -points.
661!! - Metrics centered on q-points are labelled Bu (B-grid u,v location). e.g. areaBu is the area centered on a q-point.
662!!
663!! \image html Grid_metrics.png "The labelling of distances (grid metrics) at various staggered
664!! location on an T-cell and around a q-point."
665!!
666!! Areas centered at T-, u-, v- and q- points are `areaT`, `areaCu`, `areaCv` and `areaBu` respectively.
667!!
668!! The reciprocal of metrics are pre-calculated and also stored in the ocean_grid_type with a I prepended to the name.
669!! For example, `1./areaT` is called `IareaT`, and `1./dyCv` is `IdyCv`.
670!!
671!! Geographic latitude and longitude (or model coordinates if not on a sphere) are stored in
672!! `geoLatT`, `geoLonT` for T-points.
673!! u-, v- and q- point coordinates are follow same pattern of replacing T with Cu, Cv and Bu respectively.
674!!
675!! Each location also has a 2D mask indicating whether the entire column is land or ocean.
676!! `mask2dT` is 1 if the column is wet or 0 if the T-cell is land.
677!! `mask2dCu` is 1 if both neighboring columns are ocean, and 0 if either is land.
678!! `OBCmasku` is 1 if both neighboring columns are ocean, and 0 if either is land of if this is OBC point.
679
680end module mom_grid