MOM_dyn_horgrid.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!> Contains a shareable dynamic type for describing horizontal grids and metric data
6!! and utilty routines that work on this type.
8
9use mom_array_transform, only : rotate_array, rotate_array_pair
10use mom_domains, only : mom_domain_type, deallocate_mom_domain
11use mom_error_handler, only : mom_error, mom_mesg, fatal, warning
14
15implicit none ; private
16
19
20! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
21! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
22! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
23! vary with the Boussinesq approximation, the Boussinesq variant is given first.
24
25!> Describes the horizontal ocean grid with only dynamic memory arrays
26type, public :: dyn_horgrid_type
27 type(mom_domain_type), pointer :: domain => null() !< Ocean model domain
28 type(mom_domain_type), pointer :: domain_aux => null() !< A non-symmetric auxiliary domain type.
29 type(hor_index_type) :: hi !< Horizontal index ranges
30
31 integer :: isc !< The start i-index of cell centers within the computational domain
32 integer :: iec !< The end i-index of cell centers within the computational domain
33 integer :: jsc !< The start j-index of cell centers within the computational domain
34 integer :: jec !< The end j-index of cell centers within the computational domain
35
36 integer :: isd !< The start i-index of cell centers within the data domain
37 integer :: ied !< The end i-index of cell centers within the data domain
38 integer :: jsd !< The start j-index of cell centers within the data domain
39 integer :: jed !< The end j-index of cell centers within the data domain
40
41 integer :: isg !< The start i-index of cell centers within the global domain
42 integer :: ieg !< The end i-index of cell centers within the global domain
43 integer :: jsg !< The start j-index of cell centers within the global domain
44 integer :: jeg !< The end j-index of cell centers within the global domain
45
46 integer :: iscb !< The start i-index of cell vertices within the computational domain
47 integer :: iecb !< The end i-index of cell vertices within the computational domain
48 integer :: jscb !< The start j-index of cell vertices within the computational domain
49 integer :: jecb !< The end j-index of cell vertices within the computational domain
50
51 integer :: isdb !< The start i-index of cell vertices within the data domain
52 integer :: iedb !< The end i-index of cell vertices within the data domain
53 integer :: jsdb !< The start j-index of cell vertices within the data domain
54 integer :: jedb !< The end j-index of cell vertices within the data domain
55
56 integer :: isgb !< The start i-index of cell vertices within the global domain
57 integer :: iegb !< The end i-index of cell vertices within the global domain
58 integer :: jsgb !< The start j-index of cell vertices within the global domain
59 integer :: jegb !< The end j-index of cell vertices within the global domain
60
61 integer :: isd_global !< The value of isd in the global index space (decompoistion invariant).
62 integer :: jsd_global !< The value of isd in the global index space (decompoistion invariant).
63 integer :: idg_offset !< The offset between the corresponding global and local i-indices.
64 integer :: jdg_offset !< The offset between the corresponding global and local j-indices.
65 logical :: symmetric !< True if symmetric memory is used.
66
67 logical :: nonblocking_updates !< If true, non-blocking halo updates are
68 !! allowed. The default is .false. (for now).
69 integer :: first_direction !< An integer that indicates which direction is to be updated first in
70 !! directionally split parts of the calculation. This can be altered
71 !! during the course of the run via calls to set_first_direction.
72
73 real, allocatable, dimension(:,:) :: &
74 mask2dt, & !< 0 for land points and 1 for ocean points on the h-grid [nondim].
75 geolatt, & !< The geographic latitude at q points [degrees of latitude] or [m].
76 geolont, & !< The geographic longitude at q points [degrees of longitude] or [m].
77 dxt, & !< dxT is delta x at h points [L ~> m].
78 idxt, & !< 1/dxT [L-1 ~> m-1].
79 dyt, & !< dyT is delta y at h points [L ~> m].
80 idyt, & !< IdyT is 1/dyT [L-1 ~> m-1].
81 areat, & !< The area of an h-cell [L2 ~> m2].
82 iareat !< 1/areaT [L-2 ~> m-2].
83 real, allocatable, dimension(:,:) :: sin_rot
84 !< The sine of the angular rotation between the local model grid's northward
85 !! and the true northward directions [nondim].
86 real, allocatable, dimension(:,:) :: cos_rot
87 !< The cosine of the angular rotation between the local model grid's northward
88 !! and the true northward directions [nondim].
89
90 real, allocatable, dimension(:,:) :: &
91 mask2dcu, & !< 0 for boundary points and 1 for ocean points on the u grid [nondim].
92 obcmaskcu, & !< 0 for boundary or OBC points and 1 for ocean points on the u grid [nondim].
93 geolatcu, & !< The geographic latitude at u points [degrees of latitude] or [m].
94 geoloncu, & !< The geographic longitude at u points [degrees of longitude] or [m].
95 dxcu, & !< dxCu is delta x at u points [L ~> m].
96 idxcu, & !< 1/dxCu [L-1 ~> m-1].
97 idxcu_obcmask, & !< 1/dxCu or 0 at boundary or OBC points [L-1 ~> m-1].
98 dycu, & !< dyCu is delta y at u points [L ~> m].
99 idycu, & !< 1/dyCu [L-1 ~> m-1].
100 dy_cu, & !< The unblocked lengths of the u-faces of the h-cell [L ~> m].
101 iareacu, & !< The masked inverse areas of u-grid cells [L-2 ~> m-2].
102 areacu !< The areas of the u-grid cells [L2 ~> m2].
103
104 real, allocatable, dimension(:,:) :: &
105 mask2dcv, & !< 0 for boundary points and 1 for ocean points on the v grid [nondim].
106 obcmaskcv, & !< 0 for boundary or OBC points and 1 for ocean points on the v grid [nondim].
107 geolatcv, & !< The geographic latitude at v points [degrees of latitude] or [m].
108 geoloncv, & !< The geographic longitude at v points [degrees of longitude] or [m].
109 dxcv, & !< dxCv is delta x at v points [L ~> m].
110 idxcv, & !< 1/dxCv [L-1 ~> m-1].
111 dycv, & !< dyCv is delta y at v points [L ~> m].
112 idycv, & !< 1/dyCv [L-1 ~> m-1].
113 idycv_obcmask, & !< 1/dxCv or 0 at boundary or OBC points [L-1 ~> m-1].
114 dx_cv, & !< The unblocked lengths of the v-faces of the h-cell [L ~> m].
115 iareacv, & !< The masked inverse areas of v-grid cells [L-2 ~> m-2].
116 areacv !< The areas of the v-grid cells [L2 ~> m2].
117
118 real, allocatable, dimension(:,:) :: &
119 porous_dminu, & !< minimum topographic height (deepest) of U-face [Z ~> m]
120 porous_dmaxu, & !< maximum topographic height (shallowest) of U-face [Z ~> m]
121 porous_davgu !< average topographic height of U-face [Z ~> m]
122
123 real, allocatable, dimension(:,:) :: &
124 porous_dminv, & !< minimum topographic height (deepest) of V-face [Z ~> m]
125 porous_dmaxv, & !< maximum topographic height (shallowest) of V-face [Z ~> m]
126 porous_davgv !< average topographic height of V-face [Z ~> m]
127
128 real, allocatable, dimension(:,:) :: &
129 mask2dbu, & !< 0 for boundary points and 1 for ocean points on the q grid [nondim].
130 geolatbu, & !< The geographic latitude at q points [degrees of latitude] or [m].
131 geolonbu, & !< The geographic longitude at q points [degrees of longitude] or [m].
132 dxbu, & !< dxBu is delta x at q points [L ~> m].
133 idxbu, & !< 1/dxBu [L-1 ~> m-1].
134 dybu, & !< dyBu is delta y at q points [L ~> m].
135 idybu, & !< 1/dyBu [L-1 ~> m-1].
136 areabu, & !< areaBu is the area of a q-cell [L ~> m]
137 iareabu !< IareaBu = 1/areaBu [L-2 ~> m-2].
138
139 real, pointer, dimension(:) :: gridlatt => null()
140 !< The latitude of T points for the purpose of labeling the output axes,
141 !! often in units of [degrees_N] or [km] or [m] or [gridpoints].
142 !! On many grids this is the same as geoLatT.
143 real, pointer, dimension(:) :: gridlatb => null()
144 !< The latitude of B points for the purpose of labeling the output axes,
145 !! often in units of [degrees_N] or [km] or [m] or [gridpoints].
146 !! On many grids this is the same as geoLatBu.
147 real, pointer, dimension(:) :: gridlont => null()
148 !< The longitude of T points for the purpose of labeling the output axes,
149 !! often in units of [degrees_E] or [km] or [m] or [gridpoints].
150 !! On many grids this is the same as geoLonT.
151 real, pointer, dimension(:) :: gridlonb => null()
152 !< The longitude of B points for the purpose of labeling the output axes,
153 !! often in units of [degrees_E] or [km] or [m] or [gridpoints].
154 !! On many grids this is the same as geoLonBu.
155 character(len=40) :: &
156 ! Except on a Cartesian grid, these are usually some variant of "degrees".
157 x_axis_units, & !< The units that are used in labeling the x coordinate axes.
158 y_axis_units, & !< The units that are used in labeling the y coordinate axes.
159 ! These are internally generated names, including "m", "km", "deg_E" and "deg_N".
160 x_ax_unit_short, & !< A short description of the x-axis units for documenting parameter units
161 y_ax_unit_short !< A short description of the y-axis units for documenting parameter units
162
163 real, allocatable, dimension(:,:) :: &
164 bathyt !< Ocean bottom depth, referenced to a zero reference height at tracer points.
165 !! bathyT is in depth units and positive *below* the reference height [Z ~> m].
166 real, allocatable, dimension(:,:) :: &
167 meansl !< Spatially varying time mean sea level, referenced to a zero reference height
168 !! at tracer points. meanSL is in height units and positive *above* zero. It is used
169 !! a) as the height where p = p_atm or zero;
170 !! b) to calculate time mean thickness of the water column, where
171 !! mean thickness = max(meanSL + bathyT, 0.0).
172 !! meanSL is 2D for the consideration of a domain with spatically varying mean
173 !! height, e.g. the Great Lakes system [Z ~> m].
174
175 logical :: bathymetry_at_vel !< If true, there are separate values for the
176 !! basin depths at velocity points. Otherwise the effects of
177 !! of topography are entirely determined from thickness points.
178 real, allocatable, dimension(:,:) :: &
179 dblock_u, & !< Topographic depths at u-points at which the flow is blocked [Z ~> m].
180 dopen_u !< Topographic depths at u-points at which the flow is open at width dy_Cu [Z ~> m].
181 real, allocatable, dimension(:,:) :: &
182 dblock_v, & !< Topographic depths at v-points at which the flow is blocked [Z ~> m].
183 dopen_v !< Topographic depths at v-points at which the flow is open at width dx_Cv [Z ~> m].
184 real, allocatable, dimension(:,:) :: &
185 coriolisbu, & !< The Coriolis parameter at corner points [T-1 ~> s-1].
186 coriolis2bu !< The square of the Coriolis parameter at corner points [T-2 ~> s-2].
187 real, allocatable, dimension(:,:) :: &
188 df_dx, & !< Derivative d/dx f (Coriolis parameter) at h-points [T-1 L-1 ~> s-1 m-1].
189 df_dy !< Derivative d/dy f (Coriolis parameter) at h-points [T-1 L-1 ~> s-1 m-1].
190
191 ! These variables are global sums that are useful for 1-d diagnostics.
192 real :: areat_global !< Global sum of h-cell area [L2 ~> m2]
193 real :: iareat_global !< Global sum of inverse h-cell area (1/areaT_global) [L-2 ~> m-2]
194
195 ! These parameters are run-time parameters that are used during some
196 ! initialization routines (but not all)
197 real :: grid_unit_to_l !< A factor that converts a the geoLat and geoLon variables and related
198 !! variables like len_lat and len_lon into rescaled horizontal distance
199 !! units on a Cartesian grid, in [L km ~> 1000] or [L m-1 ~> 1] or
200 !! is 0 for a non-Cartesian grid.
201 real :: south_lat !< The latitude (or y-coordinate) of the first v-line [degrees_N] or [km] or [m]
202 real :: west_lon !< The longitude (or x-coordinate) of the first u-line [degrees_E] or [km] or [m]
203 real :: len_lat !< The latitudinal (or y-coord) extent of physical domain [degrees_N] or [km] or [m]
204 real :: len_lon !< The longitudinal (or x-coord) extent of physical domain [degrees_E] or [km] or [m]
205 real :: rad_earth_l !< The radius of the planet in rescaled units [L ~> m]
206 real :: max_depth !< The maximum depth of the ocean [Z ~> m]
207end type dyn_horgrid_type
208
209contains
210
211!---------------------------------------------------------------------
212!> Allocate memory used by the dyn_horgrid_type and related structures.
213subroutine create_dyn_horgrid(G, HI, bathymetry_at_vel)
214 type(dyn_horgrid_type), pointer, intent(inout) :: g !< A pointer to the dynamic horizontal grid type
215 type(hor_index_type), intent(in) :: hi !< A hor_index_type for array extents
216 logical, optional, intent(in) :: bathymetry_at_vel !< If true, there are
217 !! separate values for the basin depths at velocity
218 !! points. Otherwise the effects of topography are
219 !! entirely determined from thickness points.
220 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb, isg, ieg, jsg, jeg
221
222 ! This subroutine allocates the lateral elements of the dyn_horgrid_type that
223 ! are always used and zeros them out.
224
225 if (associated(g)) then
226 call mom_error(warning, "create_dyn_horgrid called with an associated horgrid_type.")
227 else
228 allocate(g)
229 endif
230
231 g%HI = hi
232
233 g%isc = hi%isc ; g%iec = hi%iec ; g%jsc = hi%jsc ; g%jec = hi%jec
234 g%isd = hi%isd ; g%ied = hi%ied ; g%jsd = hi%jsd ; g%jed = hi%jed
235 g%isg = hi%isg ; g%ieg = hi%ieg ; g%jsg = hi%jsg ; g%jeg = hi%jeg
236
237 g%IscB = hi%IscB ; g%IecB = hi%IecB ; g%JscB = hi%JscB ; g%JecB = hi%JecB
238 g%IsdB = hi%IsdB ; g%IedB = hi%IedB ; g%JsdB = hi%JsdB ; g%JedB = hi%JedB
239 g%IsgB = hi%IsgB ; g%IegB = hi%IegB ; g%JsgB = hi%JsgB ; g%JegB = hi%JegB
240
241 g%idg_offset = hi%idg_offset ; g%jdg_offset = hi%jdg_offset
242 g%isd_global = g%isd + hi%idg_offset ; g%jsd_global = g%jsd + hi%jdg_offset
243 g%symmetric = hi%symmetric
244
245 g%bathymetry_at_vel = .false.
246 if (present(bathymetry_at_vel)) g%bathymetry_at_vel = bathymetry_at_vel
247
248 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
249 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
250 isg = g%isg ; ieg = g%ieg ; jsg = g%jsg ; jeg = g%jeg
251
252 allocate(g%dxT(isd:ied,jsd:jed), source=0.0)
253 allocate(g%dxCu(isdb:iedb,jsd:jed), source=0.0)
254 allocate(g%dxCv(isd:ied,jsdb:jedb), source=0.0)
255 allocate(g%dxBu(isdb:iedb,jsdb:jedb), source=0.0)
256 allocate(g%IdxT(isd:ied,jsd:jed), source=0.0)
257 allocate(g%IdxCu(isdb:iedb,jsd:jed), source=0.0)
258 allocate(g%IdxCu_OBCmask(isdb:iedb,jsd:jed), source=0.0)
259 allocate(g%IdxCv(isd:ied,jsdb:jedb), source=0.0)
260 allocate(g%IdxBu(isdb:iedb,jsdb:jedb), source=0.0)
261
262 allocate(g%dyT(isd:ied,jsd:jed), source=0.0)
263 allocate(g%dyCu(isdb:iedb,jsd:jed), source=0.0)
264 allocate(g%dyCv(isd:ied,jsdb:jedb), source=0.0)
265 allocate(g%dyBu(isdb:iedb,jsdb:jedb), source=0.0)
266 allocate(g%IdyT(isd:ied,jsd:jed), source=0.0)
267 allocate(g%IdyCu(isdb:iedb,jsd:jed), source=0.0)
268 allocate(g%IdyCv(isd:ied,jsdb:jedb), source=0.0)
269 allocate(g%IdyCv_OBCmask(isd:ied,jsdb:jedb), source=0.0)
270 allocate(g%IdyBu(isdb:iedb,jsdb:jedb), source=0.0)
271
272 allocate(g%areaT(isd:ied,jsd:jed), source=0.0)
273 allocate(g%IareaT(isd:ied,jsd:jed), source=0.0)
274 allocate(g%areaBu(isdb:iedb,jsdb:jedb), source=0.0)
275 allocate(g%IareaBu(isdb:iedb,jsdb:jedb), source=0.0)
276
277 allocate(g%mask2dT(isd:ied,jsd:jed), source=0.0)
278 allocate(g%mask2dCu(isdb:iedb,jsd:jed), source=0.0)
279 allocate(g%mask2dCv(isd:ied,jsdb:jedb), source=0.0)
280 allocate(g%mask2dBu(isdb:iedb,jsdb:jedb), source=0.0)
281 allocate(g%OBCmaskCu(isdb:iedb,jsd:jed), source=0.0)
282 allocate(g%OBCmaskCv(isd:ied,jsdb:jedb), source=0.0)
283 allocate(g%geoLatT(isd:ied,jsd:jed), source=0.0)
284 allocate(g%geoLatCu(isdb:iedb,jsd:jed), source=0.0)
285 allocate(g%geoLatCv(isd:ied,jsdb:jedb), source=0.0)
286 allocate(g%geoLatBu(isdb:iedb,jsdb:jedb), source=0.0)
287 allocate(g%geoLonT(isd:ied,jsd:jed), source=0.0)
288 allocate(g%geoLonCu(isdb:iedb,jsd:jed), source=0.0)
289 allocate(g%geoLonCv(isd:ied,jsdb:jedb), source=0.0)
290 allocate(g%geoLonBu(isdb:iedb,jsdb:jedb), source=0.0)
291
292 allocate(g%dx_Cv(isd:ied,jsdb:jedb), source=0.0)
293 allocate(g%dy_Cu(isdb:iedb,jsd:jed), source=0.0)
294
295 allocate(g%areaCu(isdb:iedb,jsd:jed), source=0.0)
296 allocate(g%areaCv(isd:ied,jsdb:jedb), source=0.0)
297 allocate(g%IareaCu(isdb:iedb,jsd:jed), source=0.0)
298 allocate(g%IareaCv(isd:ied,jsdb:jedb), source=0.0)
299
300 allocate(g%porous_DminU(isdb:iedb,jsd:jed), source=0.0)
301 allocate(g%porous_DmaxU(isdb:iedb,jsd:jed), source=0.0)
302 allocate(g%porous_DavgU(isdb:iedb,jsd:jed), source=0.0)
303
304 allocate(g%porous_DminV(isd:ied,jsdb:jedb), source=0.0)
305 allocate(g%porous_DmaxV(isd:ied,jsdb:jedb), source=0.0)
306 allocate(g%porous_DavgV(isd:ied,jsdb:jedb), source=0.0)
307
308 allocate(g%bathyT(isd:ied, jsd:jed), source=0.0)
309 allocate(g%meanSL(isd:ied, jsd:jed), source=0.0)
310 allocate(g%CoriolisBu(isdb:iedb, jsdb:jedb), source=0.0)
311 allocate(g%Coriolis2Bu(isdb:iedb, jsdb:jedb), source=0.0)
312 allocate(g%dF_dx(isd:ied, jsd:jed), source=0.0)
313 allocate(g%dF_dy(isd:ied, jsd:jed), source=0.0)
314
315 allocate(g%sin_rot(isd:ied,jsd:jed), source=0.0)
316 allocate(g%cos_rot(isd:ied,jsd:jed), source=1.0)
317
318 if (g%bathymetry_at_vel) then
319 allocate(g%Dblock_u(isdb:iedb, jsd:jed), source=0.0)
320 allocate(g%Dopen_u(isdb:iedb, jsd:jed), source=0.0)
321 allocate(g%Dblock_v(isd:ied, jsdb:jedb), source=0.0)
322 allocate(g%Dopen_v(isd:ied, jsdb:jedb), source=0.0)
323 endif
324
325 ! gridLonB and gridLatB are used as edge values in some cases, so they
326 ! always need to use symmetric memory allcoations.
327 allocate(g%gridLonT(isg:ieg), source=0.0)
328 allocate(g%gridLonB(isg-1:ieg), source=0.0)
329 allocate(g%gridLatT(jsg:jeg), source=0.0)
330 allocate(g%gridLatB(jsg-1:jeg), source=0.0)
331
332end subroutine create_dyn_horgrid
333
334
335!> Copy the rotated contents of one horizontal grid type into another. The input
336!! and output grid type arguments can not use the same object.
337subroutine rotate_dyn_horgrid(G_in, G, US, turns)
338 type(dyn_horgrid_type), intent(in) :: g_in !< The input horizontal grid type
339 type(dyn_horgrid_type), intent(inout) :: g !< An output rotated horizontal grid type
340 !! that has already been allocated, but whose
341 !! contents are largely replaced here.
342 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
343 integer, intent(in) :: turns !< Number of quarter turns
344
345 ! Center point
346 call rotate_array(g_in%geoLonT, turns, g%geoLonT)
347 call rotate_array(g_in%geoLatT, turns, g%geoLatT)
348 call rotate_array_pair(g_in%dxT, g_in%dyT, turns, g%dxT, g%dyT)
349 call rotate_array(g_in%areaT, turns, g%areaT)
350 call rotate_array(g_in%bathyT, turns, g%bathyT)
351 call rotate_array(g_in%meanSL, turns, g%meanSL)
352
353 call rotate_array_pair(g_in%df_dx, g_in%df_dy, turns, g%df_dx, g%df_dy)
354 call rotate_array(g_in%sin_rot, turns, g%sin_rot)
355 call rotate_array(g_in%cos_rot, turns, g%cos_rot)
356 call rotate_array(g_in%mask2dT, turns, g%mask2dT)
357
358 ! Face points
359 call rotate_array_pair(g_in%geoLonCu, g_in%geoLonCv, turns, g%geoLonCu, g%geoLonCv)
360 call rotate_array_pair(g_in%geoLatCu, g_in%geoLatCv, turns, g%geoLatCu, g%geoLatCv)
361 call rotate_array_pair(g_in%dxCu, g_in%dyCv, turns, g%dxCu, g%dyCv)
362 call rotate_array_pair(g_in%dxCv, g_in%dyCu, turns, g%dxCv, g%dyCu)
363 call rotate_array_pair(g_in%dx_Cv, g_in%dy_Cu, turns, g%dx_Cv, g%dy_Cu)
364
365 call rotate_array_pair(g_in%mask2dCu, g_in%mask2dCv, turns, g%mask2dCu, g%mask2dCv)
366 call rotate_array_pair(g_in%OBCmaskCu, g_in%OBCmaskCv, turns, g%OBCmaskCu, g%OBCmaskCv)
367 call rotate_array_pair(g_in%areaCu, g_in%areaCv, turns, g%areaCu, g%areaCv)
368 call rotate_array_pair(g_in%IareaCu, g_in%IareaCv, turns, g%IareaCu, g%IareaCv)
369
370 call rotate_array_pair(g_in%porous_DminU, g_in%porous_DminV, &
371 turns, g%porous_DminU, g%porous_DminV)
372 call rotate_array_pair(g_in%porous_DmaxU, g_in%porous_DmaxV, &
373 turns, g%porous_DmaxU, g%porous_DmaxV)
374 call rotate_array_pair(g_in%porous_DavgU, g_in%porous_DavgV, &
375 turns, g%porous_DavgU, g%porous_DavgV)
376
377
378 ! Vertex point
379 call rotate_array(g_in%geoLonBu, turns, g%geoLonBu)
380 call rotate_array(g_in%geoLatBu, turns, g%geoLatBu)
381 call rotate_array_pair(g_in%dxBu, g_in%dyBu, turns, g%dxBu, g%dyBu)
382 call rotate_array(g_in%areaBu, turns, g%areaBu)
383 call rotate_array(g_in%CoriolisBu, turns, g%CoriolisBu)
384 call rotate_array(g_in%Coriolis2Bu, turns, g%Coriolis2Bu)
385 call rotate_array(g_in%mask2dBu, turns, g%mask2dBu)
386
387 ! Topography at the cell faces
388 g%bathymetry_at_vel = g_in%bathymetry_at_vel
389 if (g%bathymetry_at_vel) then
390 call rotate_array_pair(g_in%Dblock_u, g_in%Dblock_v, turns, g%Dblock_u, g%Dblock_v)
391 call rotate_array_pair(g_in%Dopen_u, g_in%Dopen_v, turns, g%Dopen_u, g%Dopen_v)
392 endif
393
394 ! Nominal grid axes
395 ! TODO: We should not assign lat values to the lon axis, and vice versa.
396 ! We temporarily copy lat <-> lon since several components still expect
397 ! lat and lon sizes to match the first and second dimension sizes.
398 ! But we ought to instead leave them unchanged and adjust the references to
399 ! these axes.
400 if (modulo(turns, 2) /= 0) then
401 g%gridLonT(:) = g_in%gridLatT(g_in%jeg:g_in%jsg:-1)
402 g%gridLatT(:) = g_in%gridLonT(:)
403 g%gridLonB(:) = g_in%gridLatB(g_in%jeg:(g_in%jsg-1):-1)
404 g%gridLatB(:) = g_in%gridLonB(:)
405 else
406 g%gridLonT(:) = g_in%gridLonT(:)
407 g%gridLatT(:) = g_in%gridLatT(:)
408 g%gridLonB(:) = g_in%gridLonB(:)
409 g%gridLatB(:) = g_in%gridLatB(:)
410 endif
411
412 g%x_axis_units = g_in%y_axis_units
413 g%y_axis_units = g_in%x_axis_units
414 g%x_ax_unit_short = g_in%y_ax_unit_short
415 g%y_ax_unit_short = g_in%x_ax_unit_short
416 g%south_lat = g_in%south_lat
417 g%west_lon = g_in%west_lon
418 g%len_lat = g_in%len_lat
419 g%len_lon = g_in%len_lon
420
421 ! Rotation-invariant fields
422 g%grid_unit_to_L = g_in%grid_unit_to_L
423 g%areaT_global = g_in%areaT_global
424 g%IareaT_global = g_in%IareaT_global
425 g%Rad_Earth_L = g_in%Rad_Earth_L
426 g%max_depth = g_in%max_depth
427
428 call set_derived_dyn_horgrid(g, us)
429end subroutine rotate_dyn_horgrid
430
431
432!> rescale_dyn_horgrid_bathymetry permits a change in the internal units for the bathymetry on the
433!! grid, both rescaling the depths and recording the new internal depth units.
434subroutine rescale_dyn_horgrid_bathymetry(G, m_in_new_units)
435 type(dyn_horgrid_type), intent(inout) :: g !< The dynamic horizontal grid type
436 real, intent(in) :: m_in_new_units !< The new internal representation of 1 m depth [m Z-1 ~> 1]
437
438 ! Local variables
439 real :: rescale ! The inverse of m_in_new_units, used in rescaling bathymetry [Z m-1 ~> 1]
440 integer :: i, j, isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
441
442 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
443 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
444
445 if (m_in_new_units == 1.0) return
446 if (m_in_new_units < 0.0) &
447 call mom_error(fatal, "rescale_dyn_horgrid_bathymetry: Negative depth units are not permitted.")
448 if (m_in_new_units == 0.0) &
449 call mom_error(fatal, "rescale_dyn_horgrid_bathymetry: Zero depth units are not permitted.")
450
451 rescale = 1.0 / m_in_new_units
452 do j=jsd,jed ; do i=isd,ied
453 g%bathyT(i,j) = rescale*g%bathyT(i,j)
454 g%meanSL(i,j) = rescale*g%meanSL(i,j)
455 enddo ; enddo
456 if (g%bathymetry_at_vel) then ; do j=jsd,jed ; do i=isdb,iedb
457 g%Dblock_u(i,j) = rescale*g%Dblock_u(i,j) ; g%Dopen_u(i,j) = rescale*g%Dopen_u(i,j)
458 enddo ; enddo ; endif
459 if (g%bathymetry_at_vel) then ; do j=jsdb,jedb ; do i=isd,ied
460 g%Dblock_v(i,j) = rescale*g%Dblock_v(i,j) ; g%Dopen_v(i,j) = rescale*g%Dopen_v(i,j)
461 enddo ; enddo ; endif
462 g%max_depth = rescale*g%max_depth
463
465
466!> set_derived_dyn_horgrid calculates metric terms that are derived from other metrics.
467subroutine set_derived_dyn_horgrid(G, US)
468 type(dyn_horgrid_type), intent(inout) :: g !< The dynamic horizontal grid type
469 type(unit_scale_type), optional, intent(in) :: us !< A dimensional unit scaling type
470! Various inverse grid spacings and derived areas are calculated within this
471! subroutine.
472 integer :: i, j, isd, ied, jsd, jed
473 integer :: isdb, iedb, jsdb, jedb
474
475 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
476 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
477
478 do j=jsd,jed ; do i=isd,ied
479 if (g%dxT(i,j) < 0.0) g%dxT(i,j) = 0.0
480 if (g%dyT(i,j) < 0.0) g%dyT(i,j) = 0.0
481 g%IdxT(i,j) = adcroft_reciprocal(g%dxT(i,j))
482 g%IdyT(i,j) = adcroft_reciprocal(g%dyT(i,j))
483 g%IareaT(i,j) = adcroft_reciprocal(g%areaT(i,j))
484 enddo ; enddo
485
486 do j=jsd,jed ; do i=isdb,iedb
487 if (g%dxCu(i,j) < 0.0) g%dxCu(i,j) = 0.0
488 if (g%dyCu(i,j) < 0.0) g%dyCu(i,j) = 0.0
489 g%IdxCu(i,j) = adcroft_reciprocal(g%dxCu(i,j))
490 g%IdyCu(i,j) = adcroft_reciprocal(g%dyCu(i,j))
491 g%IdxCu_OBCmask(i,j) = g%OBCmaskCu(i,j) * g%IdxCu(i,j) ! This may be reset when the masks are set.
492 enddo ; enddo
493
494 do j=jsdb,jedb ; do i=isd,ied
495 if (g%dxCv(i,j) < 0.0) g%dxCv(i,j) = 0.0
496 if (g%dyCv(i,j) < 0.0) g%dyCv(i,j) = 0.0
497 g%IdxCv(i,j) = adcroft_reciprocal(g%dxCv(i,j))
498 g%IdyCv(i,j) = adcroft_reciprocal(g%dyCv(i,j))
499 g%IdyCv_OBCmask(i,j) = g%OBCmaskCv(i,j) * g%IdyCv(i,j) ! This may be reset when the masks are set.
500 enddo ; enddo
501
502 do j=jsdb,jedb ; do i=isdb,iedb
503 if (g%dxBu(i,j) < 0.0) g%dxBu(i,j) = 0.0
504 if (g%dyBu(i,j) < 0.0) g%dyBu(i,j) = 0.0
505
506 g%IdxBu(i,j) = adcroft_reciprocal(g%dxBu(i,j))
507 g%IdyBu(i,j) = adcroft_reciprocal(g%dyBu(i,j))
508 ! areaBu has usually been set to a positive area elsewhere.
509 if (g%areaBu(i,j) <= 0.0) g%areaBu(i,j) = g%dxBu(i,j) * g%dyBu(i,j)
510 g%IareaBu(i,j) = adcroft_reciprocal(g%areaBu(i,j))
511 enddo ; enddo
512
513end subroutine set_derived_dyn_horgrid
514
515!> Adcroft_reciprocal(x) = 1/x for |x|>0 or 0 for x=0.
516function adcroft_reciprocal(val) result(I_val)
517 real, intent(in) :: val !< The value being inverted in arbitrary units [A ~> a]
518 real :: i_val !< The Adcroft reciprocal of val [A-1 ~> a-1].
519
520 i_val = 0.0 ; if (val /= 0.0) i_val = 1.0/val
521end function adcroft_reciprocal
522
523!---------------------------------------------------------------------
524!> Release memory used by the dyn_horgrid_type and related structures.
525subroutine destroy_dyn_horgrid(G)
526 type(dyn_horgrid_type), pointer :: g !< The dynamic horizontal grid type
527
528 if (.not.associated(g)) then
529 call mom_error(fatal, "destroy_dyn_horgrid called with an unassociated horgrid_type.")
530 endif
531
532 deallocate(g%dxT) ; deallocate(g%dxCu) ; deallocate(g%dxCv) ; deallocate(g%dxBu)
533 deallocate(g%IdxT) ; deallocate(g%IdxCu) ; deallocate(g%IdxCv) ; deallocate(g%IdxBu)
534
535 deallocate(g%dyT) ; deallocate(g%dyCu) ; deallocate(g%dyCv) ; deallocate(g%dyBu)
536 deallocate(g%IdyT) ; deallocate(g%IdyCu) ; deallocate(g%IdyCv) ; deallocate(g%IdyBu)
537
538 deallocate(g%areaT) ; deallocate(g%IareaT)
539 deallocate(g%areaBu) ; deallocate(g%IareaBu)
540 deallocate(g%areaCu) ; deallocate(g%IareaCu)
541 deallocate(g%areaCv) ; deallocate(g%IareaCv)
542
543 deallocate(g%mask2dT) ; deallocate(g%mask2dCu) ; deallocate(g%OBCmaskCu)
544 deallocate(g%mask2dCv) ; deallocate(g%OBCmaskCv) ; deallocate(g%mask2dBu)
545 deallocate(g%IdxCu_OBCmask) ; deallocate(g%IdyCv_OBCmask)
546
547 deallocate(g%geoLatT) ; deallocate(g%geoLatCu)
548 deallocate(g%geoLatCv) ; deallocate(g%geoLatBu)
549 deallocate(g%geoLonT) ; deallocate(g%geoLonCu)
550 deallocate(g%geoLonCv) ; deallocate(g%geoLonBu)
551
552 deallocate(g%dx_Cv) ; deallocate(g%dy_Cu)
553
554 deallocate(g%porous_DminU) ; deallocate(g%porous_DmaxU) ; deallocate(g%porous_DavgU)
555 deallocate(g%porous_DminV) ; deallocate(g%porous_DmaxV) ; deallocate(g%porous_DavgV)
556
557 deallocate(g%bathyT) ; deallocate(g%meanSL)
558 deallocate(g%CoriolisBu) ; deallocate(g%Coriolis2Bu)
559 deallocate(g%dF_dx) ; deallocate(g%dF_dy)
560 deallocate(g%sin_rot) ; deallocate(g%cos_rot)
561
562 if (allocated(g%Dblock_u)) deallocate(g%Dblock_u)
563 if (allocated(g%Dopen_u)) deallocate(g%Dopen_u)
564 if (allocated(g%Dblock_v)) deallocate(g%Dblock_v)
565 if (allocated(g%Dopen_v)) deallocate(g%Dopen_v)
566
567 deallocate(g%gridLonT) ; deallocate(g%gridLatT)
568 deallocate(g%gridLonB) ; deallocate(g%gridLatB)
569
570 ! CS%debug is required to validate Domain_aux, so use allocation test
571 if (associated(g%Domain_aux)) call deallocate_mom_domain(g%Domain_aux)
572
573 call deallocate_mom_domain(g%Domain)
574
575 deallocate(g)
576
577end subroutine destroy_dyn_horgrid
578
579end module mom_dyn_horgrid