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) :: hid2 !< Horizontal index ranges for level-2-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
405 call get_domain_extent(g%Domain, g%HId2%isc, g%HId2%iec, g%HId2%jsc, g%HId2%jec, &
406 g%HId2%isd, g%HId2%ied, g%HId2%jsd, g%HId2%jed, &
407 g%HId2%isg, g%HId2%ieg, g%HId2%jsg, g%HId2%jeg, coarsen=2)
408
409 ! Set array sizes for fields that are discretized at tracer cell boundaries.
410 g%HId2%IscB = g%HId2%isc ; g%HId2%JscB = g%HId2%jsc
411 g%HId2%IsdB = g%HId2%isd ; g%HId2%JsdB = g%HId2%jsd
412 g%HId2%IsgB = g%HId2%isg ; g%HId2%JsgB = g%HId2%jsg
413 if (g%symmetric) then
414 g%HId2%IscB = g%HId2%isc-1 ; g%HId2%JscB = g%HId2%jsc-1
415 g%HId2%IsdB = g%HId2%isd-1 ; g%HId2%JsdB = g%HId2%jsd-1
416 g%HId2%IsgB = g%HId2%isg-1 ; g%HId2%JsgB = g%HId2%jsg-1
417 endif
418 g%HId2%IecB = g%HId2%iec ; g%HId2%JecB = g%HId2%jec
419 g%HId2%IedB = g%HId2%ied ; g%HId2%JedB = g%HId2%jed
420 g%HId2%IegB = g%HId2%ieg ; g%HId2%JegB = g%HId2%jeg
421
422end subroutine mom_grid_init
423
424!> set_derived_metrics calculates metric terms that are derived from other metrics.
425subroutine set_derived_metrics(G, US)
426 type(ocean_grid_type), intent(inout) :: g !< The horizontal grid structure
427 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
428! Various inverse grid spacings and derived areas are calculated within this
429! subroutine.
430 integer :: i, j, isd, ied, jsd, jed
431 integer :: isdb, iedb, jsdb, jedb
432
433 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
434 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
435
436 do j=jsd,jed ; do i=isd,ied
437 if (g%dxT(i,j) < 0.0) g%dxT(i,j) = 0.0
438 if (g%dyT(i,j) < 0.0) g%dyT(i,j) = 0.0
439 g%IdxT(i,j) = adcroft_reciprocal(g%dxT(i,j))
440 g%IdyT(i,j) = adcroft_reciprocal(g%dyT(i,j))
441 g%IareaT(i,j) = adcroft_reciprocal(g%areaT(i,j))
442 enddo ; enddo
443
444 do j=jsd,jed ; do i=isdb,iedb
445 if (g%dxCu(i,j) < 0.0) g%dxCu(i,j) = 0.0
446 if (g%dyCu(i,j) < 0.0) g%dyCu(i,j) = 0.0
447 g%IdxCu(i,j) = adcroft_reciprocal(g%dxCu(i,j))
448 g%IdyCu(i,j) = adcroft_reciprocal(g%dyCu(i,j))
449 g%IdxCu_OBCmask(i,j) = g%OBCmaskCu(i,j) * g%IdxCu(i,j) ! This may be reset if masks are reset.
450 enddo ; enddo
451
452 do j=jsdb,jedb ; do i=isd,ied
453 if (g%dxCv(i,j) < 0.0) g%dxCv(i,j) = 0.0
454 if (g%dyCv(i,j) < 0.0) g%dyCv(i,j) = 0.0
455 g%IdxCv(i,j) = adcroft_reciprocal(g%dxCv(i,j))
456 g%IdyCv(i,j) = adcroft_reciprocal(g%dyCv(i,j))
457 g%IdyCv_OBCmask(i,j) = g%OBCmaskCv(i,j) * g%IdyCv(i,j) ! This may be reset if masks are reset.
458 enddo ; enddo
459
460 do j=jsdb,jedb ; do i=isdb,iedb
461 if (g%dxBu(i,j) < 0.0) g%dxBu(i,j) = 0.0
462 if (g%dyBu(i,j) < 0.0) g%dyBu(i,j) = 0.0
463
464 g%IdxBu(i,j) = adcroft_reciprocal(g%dxBu(i,j))
465 g%IdyBu(i,j) = adcroft_reciprocal(g%dyBu(i,j))
466 ! areaBu has usually been set to a positive area elsewhere.
467 if (g%areaBu(i,j) <= 0.0) g%areaBu(i,j) = g%dxBu(i,j) * g%dyBu(i,j)
468 g%IareaBu(i,j) = adcroft_reciprocal(g%areaBu(i,j))
469 enddo ; enddo
470end subroutine set_derived_metrics
471
472!> Adcroft_reciprocal(x) = 1/x for |x|>0 or 0 for x=0.
473function adcroft_reciprocal(val) result(I_val)
474 real, intent(in) :: val !< The value being inverted [A].
475 real :: i_val !< The Adcroft reciprocal of val [A-1].
476
477 i_val = 0.0 ; if (val /= 0.0) i_val = 1.0/val
478end function adcroft_reciprocal
479
480!> Returns true if the coordinates (x,y) are within the h-cell (i,j)
481logical function ispointincell(G, i, j, x, y)
482 type(ocean_grid_type), intent(in) :: g !< Grid type
483 integer, intent(in) :: i !< i index of cell to test
484 integer, intent(in) :: j !< j index of cell to test
485 real, intent(in) :: x !< x coordinate of point [degrees_E]
486 real, intent(in) :: y !< y coordinate of point [degrees_N]
487 ! Local variables
488 real :: xne, xnw, xse, xsw ! Longitudes of cell corners [degrees_E]
489 real :: yne, ynw, yse, ysw ! Latitudes of cell corners [degrees_N]
490 real :: l0, l1, l2, l3 ! Crossed products of differences in position [degrees_E degrees_N]
491 real :: p0, p1, p2, p3 ! Trinary unitary values reflecting the signs of the crossed products [nondim]
492 ispointincell = .false.
493 xne = g%geoLonBu(i ,j ) ; yne = g%geoLatBu(i ,j )
494 xnw = g%geoLonBu(i-1,j ) ; ynw = g%geoLatBu(i-1,j )
495 xse = g%geoLonBu(i ,j-1) ; yse = g%geoLatBu(i ,j-1)
496 xsw = g%geoLonBu(i-1,j-1) ; ysw = g%geoLatBu(i-1,j-1)
497 ! This is a crude calculation that assumes a geographic coordinate system
498 if (x<min(xne,xnw,xse,xsw) .or. x>max(xne,xnw,xse,xsw) .or. &
499 y<min(yne,ynw,yse,ysw) .or. y>max(yne,ynw,yse,ysw) ) then
500 return ! Avoid the more complicated calculation
501 endif
502 l0 = (x-xsw)*(yse-ysw) - (y-ysw)*(xse-xsw)
503 l1 = (x-xse)*(yne-yse) - (y-yse)*(xne-xse)
504 l2 = (x-xne)*(ynw-yne) - (y-yne)*(xnw-xne)
505 l3 = (x-xnw)*(ysw-ynw) - (y-ynw)*(xsw-xnw)
506
507 p0 = sign(1., l0) ; if (l0 == 0.) p0=0.
508 p1 = sign(1., l1) ; if (l1 == 0.) p1=0.
509 p2 = sign(1., l2) ; if (l2 == 0.) p2=0.
510 p3 = sign(1., l3) ; if (l3 == 0.) p3=0.
511
512 if ( (abs(p0)+abs(p2)) + (abs(p1)+abs(p3)) == abs((p0+p2) + (p1+p3)) ) then
513 ispointincell=.true.
514 endif
515end function ispointincell
516
517!> Store an integer indicating which direction to work on first.
518subroutine set_first_direction(G, y_first)
519 type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure
520 integer, intent(in) :: y_first !< The first direction to store
521
522 g%first_direction = y_first
523end subroutine set_first_direction
524
525!> Return global shape of horizontal grid
526subroutine get_global_grid_size(G, niglobal, njglobal)
527 type(ocean_grid_type), intent(inout) :: g !< The horizontal grid type
528 integer, intent(out) :: niglobal !< i-index global size of grid
529 integer, intent(out) :: njglobal !< j-index global size of grid
530
531 call get_global_shape(g%domain, niglobal, njglobal)
532
533end subroutine get_global_grid_size
534
535!> Allocate memory used by the ocean_grid_type and related structures.
536subroutine allocate_metrics(G)
537 type(ocean_grid_type), intent(inout) :: G !< The horizontal grid type
538 integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, isg, ieg, jsg, jeg
539
540 ! This subroutine allocates the lateral elements of the ocean_grid_type that
541 ! are always used and zeros them out.
542
543 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
544 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
545 isg = g%isg ; ieg = g%ieg ; jsg = g%jsg ; jeg = g%jeg
546
547 alloc_(g%dxT(isd:ied,jsd:jed)) ; g%dxT(:,:) = 0.0
548 alloc_(g%dxCu(isdb:iedb,jsd:jed)) ; g%dxCu(:,:) = 0.0
549 alloc_(g%dxCv(isd:ied,jsdb:jedb)) ; g%dxCv(:,:) = 0.0
550 alloc_(g%dxBu(isdb:iedb,jsdb:jedb)) ; g%dxBu(:,:) = 0.0
551 alloc_(g%IdxT(isd:ied,jsd:jed)) ; g%IdxT(:,:) = 0.0
552 alloc_(g%IdxCu(isdb:iedb,jsd:jed)) ; g%IdxCu(:,:) = 0.0
553 alloc_(g%IdxCu_OBCmask(isdb:iedb,jsd:jed)) ; g%IdxCu_OBCmask(:,:) = 0.0
554 alloc_(g%IdxCv(isd:ied,jsdb:jedb)) ; g%IdxCv(:,:) = 0.0
555 alloc_(g%IdxBu(isdb:iedb,jsdb:jedb)) ; g%IdxBu(:,:) = 0.0
556
557 alloc_(g%dyT(isd:ied,jsd:jed)) ; g%dyT(:,:) = 0.0
558 alloc_(g%dyCu(isdb:iedb,jsd:jed)) ; g%dyCu(:,:) = 0.0
559 alloc_(g%dyCv(isd:ied,jsdb:jedb)) ; g%dyCv(:,:) = 0.0
560 alloc_(g%dyBu(isdb:iedb,jsdb:jedb)) ; g%dyBu(:,:) = 0.0
561 alloc_(g%IdyT(isd:ied,jsd:jed)) ; g%IdyT(:,:) = 0.0
562 alloc_(g%IdyCu(isdb:iedb,jsd:jed)) ; g%IdyCu(:,:) = 0.0
563 alloc_(g%IdyCv(isd:ied,jsdb:jedb)) ; g%IdyCv(:,:) = 0.0
564 alloc_(g%IdyCv_OBCmask(isd:ied,jsdb:jedb)) ; g%IdyCv_OBCmask(:,:) = 0.0
565 alloc_(g%IdyBu(isdb:iedb,jsdb:jedb)) ; g%IdyBu(:,:) = 0.0
566
567 alloc_(g%areaT(isd:ied,jsd:jed)) ; g%areaT(:,:) = 0.0
568 alloc_(g%IareaT(isd:ied,jsd:jed)) ; g%IareaT(:,:) = 0.0
569 alloc_(g%areaBu(isdb:iedb,jsdb:jedb)) ; g%areaBu(:,:) = 0.0
570 alloc_(g%IareaBu(isdb:iedb,jsdb:jedb)) ; g%IareaBu(:,:) = 0.0
571
572 alloc_(g%mask2dT(isd:ied,jsd:jed)) ; g%mask2dT(:,:) = 0.0
573 alloc_(g%mask2dCu(isdb:iedb,jsd:jed)) ; g%mask2dCu(:,:) = 0.0
574 alloc_(g%OBCmaskCu(isdb:iedb,jsd:jed)) ; g%OBCmaskCu(:,:) = 0.0
575 alloc_(g%mask2dCv(isd:ied,jsdb:jedb)) ; g%mask2dCv(:,:) = 0.0
576 alloc_(g%OBCmaskCv(isd:ied,jsdb:jedb)) ; g%OBCmaskCv(:,:) = 0.0
577 alloc_(g%mask2dBu(isdb:iedb,jsdb:jedb)) ; g%mask2dBu(:,:) = 0.0
578 alloc_(g%geoLatT(isd:ied,jsd:jed)) ; g%geoLatT(:,:) = 0.0
579 alloc_(g%geoLatCu(isdb:iedb,jsd:jed)) ; g%geoLatCu(:,:) = 0.0
580 alloc_(g%geoLatCv(isd:ied,jsdb:jedb)) ; g%geoLatCv(:,:) = 0.0
581 alloc_(g%geoLatBu(isdb:iedb,jsdb:jedb)) ; g%geoLatBu(:,:) = 0.0
582 alloc_(g%geoLonT(isd:ied,jsd:jed)) ; g%geoLonT(:,:) = 0.0
583 alloc_(g%geoLonCu(isdb:iedb,jsd:jed)) ; g%geoLonCu(:,:) = 0.0
584 alloc_(g%geoLonCv(isd:ied,jsdb:jedb)) ; g%geoLonCv(:,:) = 0.0
585 alloc_(g%geoLonBu(isdb:iedb,jsdb:jedb)) ; g%geoLonBu(:,:) = 0.0
586
587 alloc_(g%dx_Cv(isd:ied,jsdb:jedb)) ; g%dx_Cv(:,:) = 0.0
588 alloc_(g%dy_Cu(isdb:iedb,jsd:jed)) ; g%dy_Cu(:,:) = 0.0
589
590 alloc_(g%porous_DminU(isdb:iedb,jsd:jed)) ; g%porous_DminU(:,:) = 0.0
591 alloc_(g%porous_DmaxU(isdb:iedb,jsd:jed)) ; g%porous_DmaxU(:,:) = 0.0
592 alloc_(g%porous_DavgU(isdb:iedb,jsd:jed)) ; g%porous_DavgU(:,:) = 0.0
593
594 alloc_(g%porous_DminV(isd:ied,jsdb:jedb)) ; g%porous_DminV(:,:) = 0.0
595 alloc_(g%porous_DmaxV(isd:ied,jsdb:jedb)) ; g%porous_DmaxV(:,:) = 0.0
596 alloc_(g%porous_DavgV(isd:ied,jsdb:jedb)) ; g%porous_DavgV(:,:) = 0.0
597
598 alloc_(g%areaCu(isdb:iedb,jsd:jed)) ; g%areaCu(:,:) = 0.0
599 alloc_(g%areaCv(isd:ied,jsdb:jedb)) ; g%areaCv(:,:) = 0.0
600 alloc_(g%IareaCu(isdb:iedb,jsd:jed)) ; g%IareaCu(:,:) = 0.0
601 alloc_(g%IareaCv(isd:ied,jsdb:jedb)) ; g%IareaCv(:,:) = 0.0
602
603 alloc_(g%bathyT(isd:ied, jsd:jed)) ; g%bathyT(:,:) = -g%Z_ref
604 alloc_(g%meanSL(isd:ied, jsd:jed)) ; g%meanSL(:,:) = g%Z_ref
605 alloc_(g%CoriolisBu(isdb:iedb, jsdb:jedb)) ; g%CoriolisBu(:,:) = 0.0
606 alloc_(g%Coriolis2Bu(isdb:iedb, jsdb:jedb)) ; g%Coriolis2Bu(:,:) = 0.0
607 alloc_(g%dF_dx(isd:ied, jsd:jed)) ; g%dF_dx(:,:) = 0.0
608 alloc_(g%dF_dy(isd:ied, jsd:jed)) ; g%dF_dy(:,:) = 0.0
609
610 alloc_(g%sin_rot(isd:ied,jsd:jed)) ; g%sin_rot(:,:) = 0.0
611 alloc_(g%cos_rot(isd:ied,jsd:jed)) ; g%cos_rot(:,:) = 1.0
612
613 allocate(g%gridLonT(isg:ieg), source=0.0)
614 allocate(g%gridLonB(g%IsgB:g%IegB), source=0.0)
615 allocate(g%gridLatT(jsg:jeg), source=0.0)
616 allocate(g%gridLatB(g%JsgB:g%JegB), source=0.0)
617
618end subroutine allocate_metrics
619
620!> Release memory used by the ocean_grid_type and related structures.
621subroutine mom_grid_end(G)
622 type(ocean_grid_type), intent(inout) :: g !< The horizontal grid type
623
624 deallocate(g%Block)
625
626 if (g%bathymetry_at_vel) then
627 dealloc_(g%Dblock_u) ; dealloc_(g%Dopen_u)
628 dealloc_(g%Dblock_v) ; dealloc_(g%Dopen_v)
629 endif
630
631 dealloc_(g%dxT) ; dealloc_(g%dxCu) ; dealloc_(g%dxCv) ; dealloc_(g%dxBu)
632 dealloc_(g%IdxT) ; dealloc_(g%IdxCu) ; dealloc_(g%IdxCv) ; dealloc_(g%IdxBu)
633
634 dealloc_(g%dyT) ; dealloc_(g%dyCu) ; dealloc_(g%dyCv) ; dealloc_(g%dyBu)
635 dealloc_(g%IdyT) ; dealloc_(g%IdyCu) ; dealloc_(g%IdyCv) ; dealloc_(g%IdyBu)
636
637 dealloc_(g%IdxCu_OBCmask) ; dealloc_(g%IdyCv_OBCmask)
638
639 dealloc_(g%areaT) ; dealloc_(g%IareaT)
640 dealloc_(g%areaBu) ; dealloc_(g%IareaBu)
641 dealloc_(g%areaCu) ; dealloc_(g%IareaCu)
642 dealloc_(g%areaCv) ; dealloc_(g%IareaCv)
643
644 dealloc_(g%mask2dT) ; dealloc_(g%mask2dCu) ; dealloc_(g%OBCmaskCu)
645 dealloc_(g%mask2dCv) ; dealloc_(g%OBCmaskCv) ; dealloc_(g%mask2dBu)
646
647 dealloc_(g%geoLatT) ; dealloc_(g%geoLatCu)
648 dealloc_(g%geoLatCv) ; dealloc_(g%geoLatBu)
649 dealloc_(g%geoLonT) ; dealloc_(g%geoLonCu)
650 dealloc_(g%geoLonCv) ; dealloc_(g%geoLonBu)
651
652 dealloc_(g%dx_Cv) ; dealloc_(g%dy_Cu)
653
654 dealloc_(g%bathyT) ; dealloc_(g%meanSL)
655 dealloc_(g%CoriolisBu) ; dealloc_(g%Coriolis2Bu)
656 dealloc_(g%dF_dx) ; dealloc_(g%dF_dy)
657 dealloc_(g%sin_rot) ; dealloc_(g%cos_rot)
658
659 dealloc_(g%porous_DminU) ; dealloc_(g%porous_DmaxU) ; dealloc_(g%porous_DavgU)
660 dealloc_(g%porous_DminV) ; dealloc_(g%porous_DmaxV) ; dealloc_(g%porous_DavgV)
661
662 deallocate(g%gridLonT) ; deallocate(g%gridLatT)
663 deallocate(g%gridLonB) ; deallocate(g%gridLatB)
664
665 ! The cursory flag avoids doing any deallocation of memory in the underlying
666 ! infrastructure to avoid problems due to shared pointers.
667 call deallocate_mom_domain(g%Domain, cursory=.true.)
668
669end subroutine mom_grid_end
670
671!> \namespace mom_grid
672!!
673!! Grid metrics and their inverses are labelled according to their staggered location on a Arakawa C (or B) grid.
674!! - Metrics centered on h- or T-points are labelled T, e.g. dxT is the distance across the cell in the x-direction.
675!! - Metrics centered on u-points are labelled Cu (C-grid u location). e.g. dyCu is the y-distance between
676!! two corners of a T-cell.
677!! - Metrics centered on v-points are labelled Cv (C-grid v location). e.g. dyCv is the y-distance between two -points.
678!! - Metrics centered on q-points are labelled Bu (B-grid u,v location). e.g. areaBu is the area centered on a q-point.
679!!
680!! \image html Grid_metrics.png "The labelling of distances (grid metrics) at various staggered
681!! location on an T-cell and around a q-point."
682!!
683!! Areas centered at T-, u-, v- and q- points are `areaT`, `areaCu`, `areaCv` and `areaBu` respectively.
684!!
685!! The reciprocal of metrics are pre-calculated and also stored in the ocean_grid_type with a I prepended to the name.
686!! For example, `1./areaT` is called `IareaT`, and `1./dyCv` is `IdyCv`.
687!!
688!! Geographic latitude and longitude (or model coordinates if not on a sphere) are stored in
689!! `geoLatT`, `geoLonT` for T-points.
690!! u-, v- and q- point coordinates are follow same pattern of replacing T with Cu, Cv and Bu respectively.
691!!
692!! Each location also has a 2D mask indicating whether the entire column is land or ocean.
693!! `mask2dT` is 1 if the column is wet or 0 if the T-cell is land.
694!! `mask2dCu` is 1 if both neighboring columns are ocean, and 0 if either is land.
695!! `OBCmasku` is 1 if both neighboring columns are ocean, and 0 if either is land of if this is OBC point.
696
697end module mom_grid