MOM_diag_remap.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 runtime remapping of diagnostics to z star, sigma and
6!! rho vertical coordinates.
7!!
8!! The diag_remap_ctrl type represents a remapping of diagnostics to a particular
9!! vertical coordinate. The module is used by the diag mediator module in the
10!! following way:
11!! 1. diag_remap_init() is called to initialize a diag_remap_ctrl instance.
12!! 2. diag_remap_configure_axes() is called to read the configuration file and set up the
13!! vertical coordinate / axes definitions.
14!! 3. diag_remap_get_axes_info() returns information needed for the diag mediator to
15!! define new axes for the remapped diagnostics.
16!! 4. diag_remap_update() is called periodically (whenever h, T or S change) to either
17!! create or update the target remapping grids.
18!! 5. diag_remap_do_remap() is called from within a diag post() to do the remapping before
19!! the diagnostic is written out.
20
21
22! NOTE: In the following functions, the fields are initially passed using 1-based
23! indexing, which are then passed to separate private internal routines that shift
24! the indexing to use the same indexing conventions used elsewhere in the MOM6 code.
25!
26! * diag_remap_do_remap, which calls do_remap
27! * vertically_reintegrate_diag_field, which calls vertically_reintegrate_field
28! * vertically_interpolate_diag_field, which calls vertically_interpolate_field
29! * horizontally_average_diag_field, which calls horizontally_average_field
30
31
32module mom_diag_remap
33
34use mom_coms, only : reproducing_sum_efp, efp_to_real
35use mom_coms, only : efp_type, assignment(=), efp_sum_across_pes
36use mom_error_handler, only : mom_error, fatal, assert, warning
38use mom_diag_manager_infra,only : mom_diag_axis_init
39use mom_file_parser, only : get_param, log_param, param_file_type
41use mom_grid, only : ocean_grid_type
44use mom_eos, only : eos_type
55
56
57implicit none ; private
58
59#include "MOM_memory.h"
60
61public diag_remap_ctrl
70
71!> Represents remapping of diagnostics to a particular vertical coordinate.
72!!
73!! There is one of these types for each vertical coordinate. The vertical axes
74!! of a diagnostic will reference an instance of this type indicating how (or
75!! if) the diagnostic should be vertically remapped when being posted.
76type :: diag_remap_ctrl
77 logical :: configured = .false. !< Whether vertical coordinate has been configured
78 logical :: initialized = .false. !< Whether remappping initialized
79 logical :: used = .false. !< Whether this coordinate actually gets used.
80 integer :: vertical_coord = 0 !< The vertical coordinate that we remap to
81 character(len=10) :: vertical_coord_name ='' !< The coordinate name as understood by ALE
82 logical :: z_based_coord = .false. !< If true, this coordinate is based on remapping of
83 !! geometric distances across layers (in [Z ~> m]) rather
84 !! than layer thicknesses (in [H ~> m or kg m-2]). This
85 !! distinction only matters in non-Boussinesq mode.
86 character(len=16) :: diag_coord_name = '' !< A name for the purpose of run-time parameters
87 character(len=8) :: diag_module_suffix = '' !< The suffix for the module to appear in diag_table
88 type(remapping_cs) :: remap_cs !< Remapping control structure use for this axes
89 type(regridding_cs) :: regrid_cs !< Regridding control structure that defines the coordinates for this axes
90 integer :: nz = 0 !< Number of vertical levels used for remapping
91 real, dimension(:,:,:), allocatable :: h !< Remap grid thicknesses in [H ~> m or kg m-2] or
92 !! vertical extents in [Z ~> m], depending on the setting of Z_based_coord.
93 real, dimension(:,:,:), allocatable :: h_extensive !< Remap grid thicknesses in [H ~> m or kg m-2] or
94 !! vertical extents in [Z ~> m] for remapping extensive variables
95 integer :: interface_axes_id = 0 !< Vertical axes id for remapping at interfaces
96 integer :: layer_axes_id = 0 !< Vertical axes id for remapping on layers
97 logical :: om4_remap_via_sub_cells !< Use the OM4-era ramap_via_sub_cells
98 integer :: answer_date !< The vintage of the order of arithmetic and expressions
99 !! to use for remapping. Values below 20190101 recover
100 !! the answers from 2018, while higher values use more
101 !! robust forms of the same remapping expressions.
102
103end type diag_remap_ctrl
104
105contains
106
107!> Initialize a diagnostic remapping type with the given vertical coordinate.
108subroutine diag_remap_init(remap_cs, coord_tuple, om4_remap_via_sub_cells, answer_date, GV)
109 type(diag_remap_ctrl), intent(inout) :: remap_cs !< Diag remapping control structure
110 character(len=*), intent(in) :: coord_tuple !< A string in form of
111 !! MODULE_SUFFIX PARAMETER_SUFFIX COORDINATE_NAME
112 logical, intent(in) :: om4_remap_via_sub_cells !< Use the OM4-era ramap_via_sub_cells
113 integer, intent(in) :: answer_date !< The vintage of the order of arithmetic and expressions
114 !! to use for remapping. Values below 20190101 recover
115 !! the answers from 2018, while higher values use more
116 !! robust forms of the same remapping expressions.
117 type(verticalgrid_type), intent(in) :: gv !< The ocean vertical grid structure, used here to evaluate
118 !! whether the model is in non-Boussinesq mode.
119
120 remap_cs%diag_module_suffix = trim(extractword(coord_tuple, 1))
121 remap_cs%diag_coord_name = trim(extractword(coord_tuple, 2))
122 remap_cs%vertical_coord_name = trim(extractword(coord_tuple, 3))
123 remap_cs%vertical_coord = coordinatemode(remap_cs%vertical_coord_name)
124 remap_cs%Z_based_coord = .false.
125 if (.not.(gv%Boussinesq .or. gv%semi_Boussinesq) .and. &
126 ((remap_cs%vertical_coord == coordinatemode('ZSTAR')) .or. &
127 (remap_cs%vertical_coord == coordinatemode('SIGMA')) .or. &
128 (remap_cs%vertical_coord == coordinatemode('RHO'))) ) &
129 remap_cs%Z_based_coord = .true.
130
131 remap_cs%configured = .false.
132 remap_cs%initialized = .false.
133 remap_cs%used = .false.
134 remap_cs%om4_remap_via_sub_cells = om4_remap_via_sub_cells
135 remap_cs%answer_date = answer_date
136 remap_cs%nz = 0
137
138end subroutine diag_remap_init
139
140!> De-init a diagnostic remapping type.
141!! Free allocated memory.
142subroutine diag_remap_end(remap_cs)
143 type(diag_remap_ctrl), intent(inout) :: remap_cs !< Diag remapping control structure
144
145 if (allocated(remap_cs%h)) deallocate(remap_cs%h)
146
147 remap_cs%configured = .false.
148 remap_cs%initialized = .false.
149 remap_cs%used = .false.
150 remap_cs%nz = 0
151
152end subroutine diag_remap_end
153
154!> Inform that all diagnostics have been registered.
155!! If _set_active() has not been called on the remapping control structure
156!! will be disabled. This saves time in the case that a vertical coordinate was
157!! configured but no diagnostics which use the coordinate appeared in the
158!! diag_table.
159subroutine diag_remap_diag_registration_closed(remap_cs)
160 type(diag_remap_ctrl), intent(inout) :: remap_cs !< Diag remapping control structure
161
162 if (.not. remap_cs%used) then
163 call diag_remap_end(remap_cs)
164 call end_regridding(remap_cs%regrid_cs)
165 endif
166
168
169!> Indicate that this remapping type is actually used by the diag manager.
170!! If this is never called then the type will be disabled to save time.
171!! See further explanation with diag_remap_registration_closed.
172subroutine diag_remap_set_active(remap_cs)
173 type(diag_remap_ctrl), intent(inout) :: remap_cs !< Diag remapping control structure
174
175 remap_cs%used = .true.
176
177end subroutine diag_remap_set_active
178
179!> Configure the vertical axes for a diagnostic remapping control structure.
180!! Reads a configuration parameters to determine coordinate generation.
181subroutine diag_remap_configure_axes(remap_cs, G, GV, US, param_file)
182 type(diag_remap_ctrl), intent(inout) :: remap_cs !< Diag remap control structure
183 type(ocean_grid_type), intent(in) :: g !< The ocean's grid type
184 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
185 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
186 type(param_file_type), intent(in) :: param_file !< Parameter file structure
187
188 ! Local variables
189 character(len=40) :: mod = "MOM_diag_remap" ! This module's name.
190 character(len=8) :: units
191 character(len=34) :: longname
192 real, allocatable, dimension(:) :: &
193 interfaces, & ! Numerical values for interface vertical coordinates, in unscaled units
194 ! that might be [m], [kg m-3] or [nondim], depending on the coordinate.
195 layers ! Numerical values for layer vertical coordinates, in unscaled units
196 ! that might be [m], [kg m-3] or [nondim], depending on the coordinate.
197
198 call initialize_regridding(remap_cs%regrid_cs, g, gv, us, gv%max_depth, param_file, mod, &
199 trim(remap_cs%vertical_coord_name), "DIAG_COORD", trim(remap_cs%diag_coord_name))
200 call set_regrid_params(remap_cs%regrid_cs, min_thickness=0., integrate_downward_for_e=.false.)
201
202 remap_cs%nz = get_regrid_size(remap_cs%regrid_cs)
203
204 if (remap_cs%vertical_coord == coordinatemode('SIGMA')) then
205 units = 'nondim'
206 longname = 'Fraction'
207 elseif (remap_cs%vertical_coord == coordinatemode('RHO')) then
208 units = 'kg m-3'
209 longname = 'Target Potential Density'
210 else
211 units = 'meters'
212 longname = 'Depth'
213 endif
214
215 ! Make axes objects
216 allocate(interfaces(remap_cs%nz+1))
217 allocate(layers(remap_cs%nz))
218
219 interfaces(:) = getcoordinateinterfaces(remap_cs%regrid_cs, undo_scaling=.true.)
220 layers(:) = 0.5 * ( interfaces(1:remap_cs%nz) + interfaces(2:remap_cs%nz+1) )
221
222 remap_cs%interface_axes_id = mom_diag_axis_init(lowercase(trim(remap_cs%diag_coord_name))//'_i', &
223 interfaces, trim(units), 'z', &
224 trim(longname)//' at interface', direction=-1)
225 remap_cs%layer_axes_id = mom_diag_axis_init(lowercase(trim(remap_cs%diag_coord_name))//'_l', &
226 layers, trim(units), 'z', &
227 trim(longname)//' at cell center', direction=-1, &
228 edges=remap_cs%interface_axes_id)
229
230 ! Axes have now been configured.
231 remap_cs%configured = .true.
232
233 deallocate(interfaces)
234 deallocate(layers)
235
236end subroutine diag_remap_configure_axes
237
238!> Get layer and interface axes ids for this coordinate
239!! Needed when defining axes groups.
240subroutine diag_remap_get_axes_info(remap_cs, nz, id_layer, id_interface)
241 type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coordinate control structure
242 integer, intent(out) :: nz !< Number of vertical levels for the coordinate
243 integer, intent(out) :: id_layer !< 1D-axes id for layer points
244 integer, intent(out) :: id_interface !< 1D-axes id for interface points
245
246 nz = remap_cs%nz
247 id_layer = remap_cs%layer_axes_id
248 id_interface = remap_cs%interface_axes_id
249
250end subroutine diag_remap_get_axes_info
251
252
253!> Whether or not the axes for this vertical coordinated has been configured.
254!! Configuration is complete when diag_remap_configure_axes() has been
255!! successfully called.
256function diag_remap_axes_configured(remap_cs)
257 type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coordinate control structure
259
260 diag_remap_axes_configured = remap_cs%configured
261
262end function
263
264!> Build/update target vertical grids for diagnostic remapping.
265!! \note The target grids need to be updated whenever sea surface
266!! height or layer thicknesses changes. In the case of density-based
267!! coordinates then technically we should also regenerate the
268!! target grid whenever T/S change.
269subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state, h_target)
270 type(diag_remap_ctrl), intent(inout) :: remap_cs !< Diagnostic coordinate control structure
271 type(ocean_grid_type), pointer :: g !< The ocean's grid type
272 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
273 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
274 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
275 intent(in) :: h !< New thickness in [H ~> m or kg m-2] or [Z ~> m], depending
276 !! on the value of remap_cs%Z_based_coord
277 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
278 intent(in) :: t !< New temperatures [C ~> degC]
279 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
280 intent(in) :: s !< New salinities [S ~> ppt]
281 type(eos_type), intent(in) :: eqn_of_state !< A pointer to the equation of state
282 real, dimension(SZI_(G),SZJ_(G),remap_cs%nz), &
283 intent(inout) :: h_target !< The new diagnostic thicknesses in [H ~> m or kg m-2]
284 !! or [Z ~> m], depending on the value of remap_cs%Z_based_coord
285
286 ! Local variables
287 real, dimension(remap_cs%nz + 1) :: zinterfaces ! Interface positions [H ~> m or kg m-2] or [Z ~> m]
288 real :: h_neglect, h_neglect_edge ! Negligible thicknesses [H ~> m or kg m-2] or [Z ~> m]
289 real :: bottom_depth(szi_(g),szj_(g)) ! The depth of the bathymetry in [H ~> m or kg m-2] or [Z ~> m]
290 real :: h_tot(szi_(g),szj_(g)) ! The total thickness of the water column [H ~> m or kg m-2] or [Z ~> m]
291 real :: z_unit_scale ! A conversion factor from Z-units the internal work units in this routine,
292 ! in units of [H Z-1 ~> 1 or kg m-3] or [nondim], depending on remap_cs%Z_based_coord.
293 integer :: i, j, k, is, ie, js, je, nz
294
295 ! Note that coordinateMode('LAYER') is never 'configured' so will always return here.
296 if (.not. remap_cs%configured) return
297
298 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
299
300 ! Set the bottom depth and negligible thicknesses used in the coordinate remapping in the right units.
301 if (remap_cs%Z_based_coord) then
302 h_neglect = set_dz_neglect(gv, us, remap_cs%answer_date, h_neglect_edge)
303 z_unit_scale = 1.0
304 do j=js-1,je+1 ; do i=is-1,ie+1
305 bottom_depth(i,j) = g%bathyT(i,j) + g%Z_ref
306 enddo ; enddo
307 else
308 h_neglect = set_h_neglect(gv, remap_cs%answer_date, h_neglect_edge)
309 z_unit_scale = gv%Z_to_H ! This branch is not used in fully non-Boussinesq mode.
310 do j=js-1,je+1 ; do i=is-1,ie+1
311 bottom_depth(i,j) = gv%Z_to_H * (g%bathyT(i,j) + g%Z_ref)
312 enddo ; enddo
313 endif
314
315 if (.not. remap_cs%initialized) then
316 ! Initialize remapping and regridding on the first call
317 call initialize_remapping(remap_cs%remap_cs, 'PPM_IH4', boundary_extrapolation=.false., &
318 om4_remap_via_sub_cells=remap_cs%om4_remap_via_sub_cells, &
319 answer_date=remap_cs%answer_date, &
320 h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)
321 remap_cs%initialized = .true.
322 endif
323
324 ! Calculate the total thickness of the water column, if it is needed,
325 if ((remap_cs%vertical_coord == coordinatemode('ZSTAR')) .or. &
326 (remap_cs%vertical_coord == coordinatemode('SIGMA'))) then
327 if (remap_cs%answer_date >= 20240201) then
328 ! Avoid using sum to have a specific order for the vertical sums.
329 ! For some compilers, the explicit expression gives the same answers as the sum function.
330 h_tot(:,:) = 0.0
331 do k=1,gv%ke ; do j=js-1,je+1 ; do i=is-1,ie+1
332 h_tot(i,j) = h_tot(i,j) + h(i,j,k)
333 enddo ; enddo ; enddo
334 else
335 do j=js-1,je+1 ; do i=is-1,ie+1
336 h_tot(i,j) = sum(h(i,j,:))
337 enddo ; enddo
338 endif
339 endif
340
341 ! Calculate remapping thicknesses for different target grids based on
342 ! nominal/target interface locations. This happens for every call on the
343 ! assumption that h, T, S has changed.
344 h_target(:,:,:) = 0.0
345
346 nz = remap_cs%nz
347 if (remap_cs%vertical_coord == coordinatemode('ZSTAR')) then
348 do j=js-1,je+1 ; do i=is-1,ie+1 ; if (g%mask2dT(i,j) > 0.0) then
349 ! This function call can work with the last 4 arguments all in units of [Z ~> m] or [H ~> kg m-2].
350 call build_zstar_column(get_zlike_cs(remap_cs%regrid_cs), &
351 bottom_depth(i,j), h_tot(i,j), zinterfaces, zscale=z_unit_scale)
352 do k=1,nz ; h_target(i,j,k) = zinterfaces(k) - zinterfaces(k+1) ; enddo
353 endif ; enddo ; enddo
354 elseif (remap_cs%vertical_coord == coordinatemode('SIGMA')) then
355 do j=js-1, je+1 ; do i=is-1,ie+1 ; if (g%mask2dT(i,j) > 0.0) then
356 ! This function call can work with the last 3 arguments all in units of [Z ~> m] or [H ~> kg m-2].
357 call build_sigma_column(get_sigma_cs(remap_cs%regrid_cs), &
358 bottom_depth(i,j), h_tot(i,j), zinterfaces)
359 do k=1,nz ; h_target(i,j,k) = zinterfaces(k) - zinterfaces(k+1) ; enddo
360 endif ; enddo ; enddo
361 elseif (remap_cs%vertical_coord == coordinatemode('RHO')) then
362 do j=js-1,je+1 ; do i=is-1,ie+1 ; if (g%mask2dT(i,j) > 0.0) then
363 ! This function call can work with 5 arguments in units of [Z ~> m] or [H ~> kg m-2].
364 call build_rho_column(get_rho_cs(remap_cs%regrid_cs), gv%ke, &
365 bottom_depth(i,j), h(i,j,:), t(i,j,:), s(i,j,:), &
366 eqn_of_state, zinterfaces, h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)
367 do k=1,nz ; h_target(i,j,k) = zinterfaces(k) - zinterfaces(k+1) ; enddo
368 endif ; enddo ; enddo
369 elseif (remap_cs%vertical_coord == coordinatemode('HYCOM1')) then
370 call mom_error(fatal,"diag_remap_update: HYCOM1 coordinate not coded for diagnostics yet!")
371! do j=js-1,je+1 ; do i=is-1,ie+1 ; if (G%mask2dT(i,j) > 0.0) then
372! call build_hycom1_column(remap_cs%regrid_cs, nz, &
373! bottom_depth(i,j), h_tot(i,j), zInterfaces)
374! do k=1,nz ; h_target(i,j,k) = zInterfaces(K) - zInterfaces(K+1) ; enddo
375! endif ; enddo ; enddo
376 endif
377
378end subroutine diag_remap_update
379
380!> Remap diagnostic field to alternative vertical grid.
381subroutine diag_remap_do_remap(remap_cs, G, GV, US, h, OBC_u, OBC_v, staggered_in_x, staggered_in_y, &
382 mask, field, remapped_field)
383 type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coordinate control structure
384 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
385 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
386 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
387 real, dimension(:,:,:), intent(in) :: h !< The current thicknesses [H ~> m or kg m-2] or [Z ~> m],
388 !! depending on the value of remap_CS%Z_based_coord
389 integer, dimension(:,:), intent(in) :: obc_u !< An array that indicates the presence and direction
390 !! of any open boundary conditions at u-points,
391 !! with a value of 0 for no OBC, 1 for an Eastern OBC
392 !! or -1 for a Western OBC
393 integer, dimension(:,:), intent(in) :: obc_v !< An array that indicates the presence and direction
394 !! of any open boundary conditions at v-points,
395 !! with a value of 0 for no OBC, 1 for a Northern OBC
396 !! or -1 for a Southern OBC
397 logical, intent(in) :: staggered_in_x !< True is the x-axis location is at u or q points
398 logical, intent(in) :: staggered_in_y !< True is the y-axis location is at v or q points
399 real, dimension(:,:,:), pointer :: mask !< A mask for the field [nondim].
400 real, dimension(:,:,:), intent(in) :: field(:,:,:) !< The diagnostic field to be remapped [A]
401 real, dimension(:,:,:), intent(out) :: remapped_field !< Field remapped to new coordinate [A]
402
403 ! Local variables
404 integer :: isdf, jsdf !< The starting i- and j-indices in memory for field
405
406 call assert(remap_cs%initialized, 'diag_remap_do_remap: remap_cs not initialized.')
407 call assert(size(field, 3) == size(h, 3), &
408 'diag_remap_do_remap: Remap field and thickness z-axes do not match.')
409
410 isdf = g%isd ; if (staggered_in_x) isdf = g%IsdB
411 jsdf = g%jsd ; if (staggered_in_y) jsdf = g%JsdB
412
413 if (associated(mask)) then
414 call do_remap(remap_cs, g, gv, us, isdf, jsdf, h, obc_u, obc_v, staggered_in_x, staggered_in_y, &
415 field, remapped_field, mask(:,:,1))
416 else
417 call do_remap(remap_cs, g, gv, us, isdf, jsdf, h, obc_u, obc_v, staggered_in_x, staggered_in_y, &
418 field, remapped_field)
419 endif
420
421end subroutine diag_remap_do_remap
422
423!> The internal routine to remap a diagnostic field to an alternative vertical grid.
424subroutine do_remap(remap_cs, G, GV, US, isdf, jsdf, h, OBC_u, OBC_v, staggered_in_x, staggered_in_y, &
425 field, remapped_field, mask)
426 type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coordinate control structure
427 type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
428 type(verticalgrid_type), intent(in) :: GV !< ocean vertical grid structure
429 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
430 integer, intent(in) :: isdf !< The starting i-index in memory for field
431 integer, intent(in) :: jsdf !< The starting j-index in memory for field
432 real, dimension(G%isd:,G%jsd:,:), &
433 intent(in) :: h !< The current thicknesses [H ~> m or kg m-2] or [Z ~> m],
434 !! depending on the value of remap_CS%Z_based_coord
435 integer, dimension(G%IsdB:,G%jsd:), &
436 intent(in) :: OBC_u !< An array that indicates the presence and direction
437 !! of any open boundary conditions at u-points,
438 !! with a value of 0 for no OBC, 1 for an Eastern OBC
439 !! or -1 for a Western OBC
440 integer, dimension(G%isd:,G%JsdB:), &
441 intent(in) :: OBC_v !< An array that indicates the presence and direction
442 !! of any open boundary conditions at v-points,
443 !! with a value of 0 for no OBC, 1 for a Northern OBC
444 !! or -1 for a Southern OBC
445 logical, intent(in) :: staggered_in_x !< True is the x-axis location is at u or q points
446 logical, intent(in) :: staggered_in_y !< True is the y-axis location is at v or q points
447 real, dimension(isdf:,jsdf:,:), &
448 intent(in) :: field !< The diagnostic field to be remapped [A]
449 real, dimension(isdf:,jsdf:,:), &
450 intent(out) :: remapped_field !< Field remapped to new coordinate [A]
451 real, dimension(isdf:,jsdf:), &
452 optional, intent(in) :: mask !< A mask for the field [nondim]
453
454 ! Local variables
455 real, dimension(remap_cs%nz) :: h_dest ! Destination thicknesses [H ~> m or kg m-2] or [Z ~> m]
456 real, dimension(size(h,3)) :: h_src ! A column of source thicknesses [H ~> m or kg m-2] or [Z ~> m]
457 integer :: nz_src, nz_dest ! The number of layers on the native and remapped grids
458 integer :: i, j ! Grid index
459
460 nz_src = size(field,3)
461 nz_dest = remap_cs%nz
462 remapped_field(:,:,:) = 0.
463
464 if (staggered_in_x .and. .not. staggered_in_y) then
465 ! U-points
466 if (present(mask)) then
467 do j=g%jsc,g%jec ; do i=g%IscB,g%IecB ; if (mask(i,j) > 0.) then
468 if (obc_u(i,j) == 0) then ! This is not an OBC face.
469 h_src(:) = 0.5*(h(i,j,:) + h(i+1,j,:))
470 h_dest(:) = 0.5*(remap_cs%h(i,j,:) + remap_cs%h(i+1,j,:))
471 elseif (obc_u(i,j) < 0) then ! This is a western OBC face.
472 h_src(:) = h(i+1,j,:)
473 h_dest(:) = remap_cs%h(i+1,j,:)
474 else ! (OBC_u(I,j) > 0) ! This is a eastern OBC face.
475 h_src(:) = h(i,j,:)
476 h_dest(:) = remap_cs%h(i,j,:)
477 endif
478 call remapping_core_h(remap_cs%remap_cs, nz_src, h_src(:), field(i,j,:), &
479 nz_dest, h_dest(:), remapped_field(i,j,:))
480 endif ; enddo ; enddo
481 else
482 do j=g%jsc,g%jec ; do i=g%IscB,g%IecB
483 if (obc_u(i,j) == 0) then ! This is not an OBC face.
484 h_src(:) = 0.5*(h(i,j,:) + h(i+1,j,:))
485 h_dest(:) = 0.5*(remap_cs%h(i,j,:) + remap_cs%h(i+1,j,:))
486 elseif (obc_u(i,j) < 0) then ! This is a western OBC face.
487 h_src(:) = h(i+1,j,:)
488 h_dest(:) = remap_cs%h(i+1,j,:)
489 else ! (OBC_u(I,j) > 0) ! This is a eastern OBC face.
490 h_src(:) = h(i,j,:)
491 h_dest(:) = remap_cs%h(i,j,:)
492 endif
493 call remapping_core_h(remap_cs%remap_cs, nz_src, h_src(:), field(i,j,:), &
494 nz_dest, h_dest(:), remapped_field(i,j,:))
495 enddo ; enddo
496 endif
497 elseif (staggered_in_y .and. .not. staggered_in_x) then
498 ! V-points
499 if (present(mask)) then
500 do j=g%jscB,g%jecB ; do i=g%isc,g%iec ; if (mask(i,j) > 0.) then
501 if (obc_v(i,j) == 0) then ! This is not an OBC face.
502 h_src(:) = 0.5*(h(i,j,:) + h(i,j+1,:))
503 h_dest(:) = 0.5*(remap_cs%h(i,j,:) + remap_cs%h(i,j+1,:))
504 elseif (obc_v(i,j) < 0) then ! This is a southern OBC face
505 h_src(:) = h(i,j+1,:)
506 h_dest(:) = remap_cs%h(i,j+1,:)
507 else ! (OBC_v(i,J) > 0) ! This is a northern OBC face
508 h_src(:) = h(i,j,:)
509 h_dest(:) = remap_cs%h(i,j,:)
510 endif
511 call remapping_core_h(remap_cs%remap_cs, nz_src, h_src(:), field(i,j,:), &
512 nz_dest, h_dest(:), remapped_field(i,j,:))
513 endif ; enddo ; enddo
514 else
515 do j=g%jscB,g%jecB ; do i=g%isc,g%iec
516 if (obc_v(i,j) == 0) then ! This is not an OBC face.
517 h_src(:) = 0.5*(h(i,j,:) + h(i,j+1,:))
518 h_dest(:) = 0.5*(remap_cs%h(i,j,:) + remap_cs%h(i,j+1,:))
519 elseif (obc_v(i,j) < 0) then ! This is a southern OBC face
520 h_src(:) = h(i,j+1,:)
521 h_dest(:) = remap_cs%h(i,j+1,:)
522 else ! (OBC_v(i,J) > 0) ! This is a northern OBC face
523 h_src(:) = h(i,j,:)
524 h_dest(:) = remap_cs%h(i,j,:)
525 endif
526 call remapping_core_h(remap_cs%remap_cs, nz_src, h_src(:), field(i,j,:), &
527 nz_dest, h_dest(:), remapped_field(i,j,:))
528 enddo ; enddo
529 endif
530 elseif ((.not. staggered_in_x) .and. (.not. staggered_in_y)) then
531 ! H-points
532 if (present(mask)) then
533 do j=g%jsc,g%jec ; do i=g%isc,g%iec ; if (mask(i,j) > 0.) then
534 call remapping_core_h(remap_cs%remap_cs, nz_src, h(i,j,:), field(i,j,:), &
535 nz_dest, remap_cs%h(i,j,:), remapped_field(i,j,:))
536 endif ; enddo ; enddo
537 else
538 do j=g%jsc,g%jec ; do i=g%isc,g%iec
539 call remapping_core_h(remap_cs%remap_cs, nz_src, h(i,j,:), field(i,j,:), &
540 nz_dest, remap_cs%h(i,j,:), remapped_field(i,j,:))
541 enddo ; enddo
542 endif
543 else
544 call assert(.false., 'diag_remap_do_remap: Unsupported axis combination')
545 endif
546
547end subroutine do_remap
548
549!> Calculate masks for target grid
550subroutine diag_remap_calc_hmask(remap_cs, G, mask)
551 type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coordinate control structure
552 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
553 real, dimension(G%isd:,G%jsd:,:), &
554 intent(out) :: mask !< h-point mask for target grid [nondim]
555
556 ! Local variables
557 real, dimension(remap_cs%nz) :: h_dest ! Destination thicknesses [H ~> m or kg m-2] or [Z ~> m]
558 integer :: i, j, k
559 logical :: mask_vanished_layers
560 real :: h_tot ! Sum of all thicknesses [H ~> m or kg m-2] or [Z ~> m]
561 real :: h_err ! An estimate of a negligible thickness [H ~> m or kg m-2] or [Z ~> m]
562
563 call assert(remap_cs%initialized, 'diag_remap_calc_hmask: remap_cs not initialized.')
564
565 ! Only z*-like diagnostic coordinates should have a 3d mask
566 mask_vanished_layers = (remap_cs%vertical_coord == coordinatemode('ZSTAR'))
567 mask(:,:,:) = 0.
568
569 do j=g%jsc-1,g%jec+1 ; do i=g%isc-1,g%iec+1
570 if (g%mask2dT(i,j)>0.) then
571 if (mask_vanished_layers) then
572 h_dest(:) = remap_cs%h(i,j,:)
573 h_tot = 0.
574 h_err = 0.
575 do k=1, remap_cs%nz
576 h_tot = h_tot + h_dest(k)
577 ! This is an overestimate of how thick a vanished layer might be, that
578 ! appears due to round-off.
579 h_err = h_err + epsilon(h_tot) * h_tot
580 ! Mask out vanished layers
581 if (h_dest(k)<=8.*h_err) then
582 mask(i,j,k) = 0.
583 else
584 mask(i,j,k) = 1.
585 endif
586 enddo
587 else ! all layers might contain data
588 mask(i,j,:) = 1.
589 endif
590 endif
591 enddo ; enddo
592
593end subroutine diag_remap_calc_hmask
594
595!> Vertically re-grid an already vertically-integrated diagnostic field to alternative vertical grid.
596subroutine vertically_reintegrate_diag_field(remap_cs, G, h, h_target, OBC_u, OBC_v, &
597 staggered_in_x, staggered_in_y, mask, field, reintegrated_field)
598 type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coordinate control structure
599 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
600 real, dimension(:,:,:), intent(in) :: h !< The thicknesses of the source grid [H ~> m or kg m-2] or [Z ~> m]
601 real, dimension(:,:,:), intent(in) :: h_target !< The thicknesses of the target grid [H ~> m or kg m-2] or [Z ~> m]
602 integer, dimension(:,:), intent(in) :: obc_u !< An array that indicates the presence and direction
603 !! of any open boundary conditions at u-points,
604 !! with a value of 0 for no OBC, 1 for an Eastern OBC
605 !! or -1 for a Western OBC
606 integer, dimension(:,:), intent(in) :: obc_v !< An array that indicates the presence and direction
607 !! of any open boundary conditions at v-points,
608 !! with a value of 0 for no OBC, 1 for a Northern OBC
609 !! or -1 for a Southern OBC
610 logical, intent(in) :: staggered_in_x !< True is the x-axis location is at u or q points
611 logical, intent(in) :: staggered_in_y !< True is the y-axis location is at v or q points
612 real, dimension(:,:,:), pointer :: mask !< A mask for the field [nondim]. Note that because this
613 !! is a pointer it retains its declared indexing conventions.
614 real, dimension(:,:,:), intent(in) :: field !< The diagnostic field to be remapped [A]
615 real, dimension(:,:,:), intent(out) :: reintegrated_field !< Field argument remapped to alternative coordinate [A]
616
617 ! Local variables
618 integer :: isdf, jsdf !< The starting i- and j-indices in memory for field
619
620 call assert(remap_cs%initialized, 'vertically_reintegrate_diag_field: remap_cs not initialized.')
621 call assert(size(field, 3) == size(h, 3), &
622 'vertically_reintegrate_diag_field: Remap field and thickness z-axes do not match.')
623
624 isdf = g%isd ; if (staggered_in_x) isdf = g%IsdB
625 jsdf = g%jsd ; if (staggered_in_y) jsdf = g%JsdB
626
627 if (associated(mask)) then
628 call vertically_reintegrate_field(remap_cs, g, isdf, jsdf, h, h_target, obc_u, obc_v, &
629 staggered_in_x, staggered_in_y, field, reintegrated_field, mask(:,:,1))
630 else
631 call vertically_reintegrate_field(remap_cs, g, isdf, jsdf, h, h_target, obc_u, obc_v, &
632 staggered_in_x, staggered_in_y, field, reintegrated_field)
633 endif
634
636
637!> The internal routine to vertically re-grid an already vertically-integrated diagnostic field to
638!! an alternative vertical grid.
639subroutine vertically_reintegrate_field(remap_cs, G, isdf, jsdf, h, h_target, OBC_u, OBC_v, &
640 staggered_in_x, staggered_in_y, field, reintegrated_field, mask)
641 type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coordinate control structure
642 type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
643 integer, intent(in) :: isdf !< The starting i-index in memory for field
644 integer, intent(in) :: jsdf !< The starting j-index in memory for field
645 real, dimension(G%isd:,G%jsd:,:), &
646 intent(in) :: h !< The thicknesses of the source grid [H ~> m or kg m-2] or [Z ~> m]
647 real, dimension(G%isd:,G%jsd:,:), &
648 intent(in) :: h_target !< The thicknesses of the target grid [H ~> m or kg m-2] or [Z ~> m]
649 integer, dimension(G%IsdB:,G%jsd:), &
650 intent(in) :: OBC_u !< An array that indicates the presence and direction
651 !! of any open boundary conditions at u-points,
652 !! with a value of 0 for no OBC, 1 for an Eastern OBC
653 !! or -1 for a Western OBC
654 integer, dimension(G%isd:,G%JsdB:), &
655 intent(in) :: OBC_v !< An array that indicates the presence and direction
656 !! of any open boundary conditions at v-points,
657 !! with a value of 0 for no OBC, 1 for a Northern OBC
658 !! or -1 for a Southern OBC
659 logical, intent(in) :: staggered_in_x !< True is the x-axis location is at u or q points
660 logical, intent(in) :: staggered_in_y !< True is the y-axis location is at v or q points
661 real, dimension(isdf:,jsdf:,:), &
662 intent(in) :: field !< The diagnostic field to be remapped [A]
663 real, dimension(isdf:,jsdf:,:), &
664 intent(out) :: reintegrated_field !< Field argument remapped to alternative coordinate [A]
665 real, dimension(isdf:,jsdf:), &
666 optional, intent(in) :: mask !< A mask for the field [nondim]
667
668 ! Local variables
669 real, dimension(remap_cs%nz) :: h_dest ! Destination thicknesses [H ~> m or kg m-2] or [Z ~> m]
670 real, dimension(size(h,3)) :: h_src ! A column of source thicknesses [H ~> m or kg m-2] or [Z ~> m]
671 integer :: nz_src, nz_dest ! The number of layers on the native and remapped grids
672 integer :: i, j ! Grid index
673
674 nz_src = size(field,3)
675 nz_dest = remap_cs%nz
676 reintegrated_field(:,:,:) = 0.
677
678 if (staggered_in_x .and. .not. staggered_in_y) then
679 ! U-points
680 if (present(mask)) then
681 do j=g%jsc,g%jec ; do i=g%IscB,g%IecB ; if (mask(i,j) > 0.0) then
682 if (obc_u(i,j) == 0) then ! This is not an OBC face.
683 h_src(:) = 0.5*(h(i,j,:) + h(i+1,j,:))
684 h_dest(:) = 0.5*(h_target(i,j,:) + h_target(i+1,j,:))
685 elseif (obc_u(i,j) < 0) then ! This is a western OBC face
686 h_src(:) = h(i+1,j,:)
687 h_dest(:) = h_target(i+1,j,:)
688 else ! (OBC_u(I,j) > 0) ! This is an eastern OBC face
689 h_src(:) = h(i,j,:)
690 h_dest(:) = h_target(i,j,:)
691 endif
692 call reintegrate_column(nz_src, h_src, field(i,j,:), &
693 nz_dest, h_dest, reintegrated_field(i,j,:))
694 endif ; enddo ; enddo
695 else
696 do j=g%jsc,g%jec ; do i=g%IscB,g%IecB
697 if (obc_u(i,j) == 0) then ! This is not an OBC face.
698 h_src(:) = 0.5*(h(i,j,:) + h(i+1,j,:))
699 h_dest(:) = 0.5*(h_target(i,j,:) + h_target(i+1,j,:))
700 elseif (obc_u(i,j) < 0) then ! This is a western OBC face
701 h_src(:) = h(i+1,j,:)
702 h_dest(:) = h_target(i+1,j,:)
703 else ! (OBC_u(I,j) > 0) ! This is an eastern OBC face
704 h_src(:) = h(i,j,:)
705 h_dest(:) = h_target(i,j,:)
706 endif
707 call reintegrate_column(nz_src, h_src, field(i,j,:), &
708 nz_dest, h_dest, reintegrated_field(i,j,:))
709 enddo ; enddo
710 endif
711 elseif (staggered_in_y .and. .not. staggered_in_x) then
712 ! V-points
713 if (present(mask)) then
714 do j=g%jscB,g%jecB ; do i=g%isc,g%iec ; if (mask(i,j) > 0.0) then
715 if (obc_v(i,j) == 0) then ! This is not an OBC face.
716 h_src(:) = 0.5*(h(i,j,:) + h(i,j+1,:))
717 h_dest(:) = 0.5*(h_target(i,j,:) + h_target(i,j+1,:))
718 elseif (obc_v(i,j) < 0) then ! This is a southern OBC face
719 h_src(:) = h(i,j+1,:)
720 h_dest(:) = h_target(i,j+1,:)
721 else ! (OBC_v(i,J) > 0) ! This is a northern OBC face
722 h_src(:) = h(i,j,:)
723 h_dest(:) = h_target(i,j,:)
724 endif
725 call reintegrate_column(nz_src, h_src, field(i,j,:), &
726 nz_dest, h_dest, reintegrated_field(i,j,:))
727 endif ; enddo ; enddo
728 else
729 do j=g%jscB,g%jecB ; do i=g%isc,g%iec
730 if (obc_v(i,j) == 0) then ! This is not an OBC face.
731 h_src(:) = 0.5*(h(i,j,:) + h(i,j+1,:))
732 h_dest(:) = 0.5*(h_target(i,j,:) + h_target(i,j+1,:))
733 elseif (obc_v(i,j) < 0) then ! This is a southern OBC face
734 h_src(:) = h(i,j+1,:)
735 h_dest(:) = h_target(i,j+1,:)
736 else ! (OBC_v(i,J) > 0) ! This is a northern OBC face
737 h_src(:) = h(i,j,:)
738 h_dest(:) = h_target(i,j,:)
739 endif
740 call reintegrate_column(nz_src, h_src, field(i,j,:), &
741 nz_dest, h_dest, reintegrated_field(i,j,:))
742 enddo ; enddo
743 endif
744 elseif ((.not. staggered_in_x) .and. (.not. staggered_in_y)) then
745 ! H-points
746 if (present(mask)) then
747 do j=g%jsc,g%jec ; do i=g%isc,g%iec ; if (mask(i,j) > 0.0) then
748 call reintegrate_column(nz_src, h(i,j,:), field(i,j,:), &
749 nz_dest, h_target(i,j,:), reintegrated_field(i,j,:))
750 endif ; enddo ; enddo
751 else
752 do j=g%jsc,g%jec ; do i=g%isc,g%iec
753 call reintegrate_column(nz_src, h(i,j,:), field(i,j,:), &
754 nz_dest, h_target(i,j,:), reintegrated_field(i,j,:))
755 enddo ; enddo
756 endif
757 else
758 call assert(.false., 'vertically_reintegrate_diag_field: Q point remapping is not coded yet.')
759 endif
760
761end subroutine vertically_reintegrate_field
762
763!> Vertically interpolate diagnostic field to alternative vertical grid.
764subroutine vertically_interpolate_diag_field(remap_cs, G, h, OBC_u, OBC_v, staggered_in_x, staggered_in_y, &
765 mask, field, interpolated_field)
766 type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coordinate control structure
767 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
768 real, dimension(:,:,:), intent(in) :: h !< The current thicknesses [H ~> m or kg m-2] or [Z ~> m],
769 !! depending on the value of remap_cs%Z_based_coord
770 integer, dimension(:,:), intent(in) :: obc_u !< An array that indicates the presence and direction
771 !! of any open boundary conditions at u-points,
772 !! with a value of 0 for no OBC, 1 for an Eastern OBC
773 !! or -1 for a Western OBC
774 integer, dimension(:,:), intent(in) :: obc_v !< An array that indicates the presence and direction
775 !! of any open boundary conditions at v-points,
776 !! with a value of 0 for no OBC, 1 for a Northern OBC
777 !! or -1 for a Southern OBC
778 logical, intent(in) :: staggered_in_x !< True is the x-axis location is at u or q points
779 logical, intent(in) :: staggered_in_y !< True is the y-axis location is at v or q points
780 real, dimension(:,:,:), pointer :: mask !< A mask for the field [nondim]. Note that because this
781 !! is a pointer it retains its declared indexing conventions.
782 real, dimension(:,:,:), intent(in) :: field !< The diagnostic field to be remapped [A]
783 real, dimension(:,:,:), intent(inout) :: interpolated_field !< Field argument remapped to alternative coordinate [A]
784
785 ! Local variables
786 integer :: isdf, jsdf !< The starting i- and j-indices in memory for field
787
788 call assert(remap_cs%initialized, 'vertically_interpolate_diag_field: remap_cs not initialized.')
789 call assert(size(field, 3) == size(h, 3)+1, &
790 'vertically_interpolate_diag_field: Remap field and thickness z-axes do not match.')
791
792 isdf = g%isd ; if (staggered_in_x) isdf = g%IsdB
793 jsdf = g%jsd ; if (staggered_in_y) jsdf = g%JsdB
794
795 if (associated(mask)) then
796 call vertically_interpolate_field(remap_cs, g, isdf, jsdf, h, obc_u, obc_v, staggered_in_x, staggered_in_y, &
797 field, interpolated_field, mask(:,:,1))
798 else
799 call vertically_interpolate_field(remap_cs, g, isdf, jsdf, h, obc_u, obc_v, staggered_in_x, staggered_in_y, &
800 field, interpolated_field)
801 endif
802
804
805!> Internal routine to vertically interpolate a diagnostic field to an alternative vertical grid.
806subroutine vertically_interpolate_field(remap_cs, G, isdf, jsdf, h, OBC_u, OBC_v, &
807 staggered_in_x, staggered_in_y, field, interpolated_field, mask)
808 type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coordinate control structure
809 type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
810 integer, intent(in) :: isdf !< The starting i-index in memory for field
811 integer, intent(in) :: jsdf !< The starting j-index in memory for field
812 real, dimension(G%isd:,G%jsd:,:), &
813 intent(in) :: h !< The current thicknesses [H ~> m or kg m-2] or [Z ~> m],
814 !! depending on the value of remap_cs%Z_based_coord
815 integer, dimension(G%IsdB:,G%jsd:), &
816 intent(in) :: OBC_u !< An array that indicates the presence and direction
817 !! of any open boundary conditions at u-points,
818 !! with a value of 0 for no OBC, 1 for an Eastern OBC
819 !! or -1 for a Western OBC
820 integer, dimension(G%isd:,G%JsdB:), &
821 intent(in) :: OBC_v !< An array that indicates the presence and direction
822 !! of any open boundary conditions at v-points,
823 !! with a value of 0 for no OBC, 1 for a Northern OBC
824 !! or -1 for a Southern OBC
825 logical, intent(in) :: staggered_in_x !< True is the x-axis location is at u or q points
826 logical, intent(in) :: staggered_in_y !< True is the y-axis location is at v or q points
827 real, dimension(isdf:,jsdf:,:), &
828 intent(in) :: field !< The diagnostic field to be remapped [A]
829 real, dimension(isdf:,jsdf:,:), &
830 intent(out) :: interpolated_field !< Field argument remapped to alternative coordinate [A]
831 real, dimension(isdf:,jsdf:), &
832 optional, intent(in) :: mask !< A mask for the field [nondim]
833
834 ! Local variables
835 real, dimension(remap_cs%nz) :: h_dest ! Destination thicknesses [H ~> m or kg m-2] or [Z ~> m]
836 real, dimension(size(h,3)) :: h_src ! A column of source thicknesses [H ~> m or kg m-2] or [Z ~> m]
837 integer :: nz_src, nz_dest ! The number of layers on the native and remapped grids
838 integer :: i, j !< Grid index
839
840 interpolated_field(:,:,:) = 0.
841
842 nz_src = size(h,3)
843 nz_dest = remap_cs%nz
844
845 if (staggered_in_x .and. .not. staggered_in_y) then
846 ! U-points
847 if (present(mask)) then
848 do j=g%jsc,g%jec ; do i=g%IscB,g%IecB ; if (mask(i,j) > 0.0) then
849 if (obc_u(i,j) == 0) then ! This is not an OBC face.
850 h_src(:) = 0.5*(h(i,j,:) + h(i+1,j,:))
851 h_dest(:) = 0.5*(remap_cs%h(i,j,:) + remap_cs%h(i+1,j,:))
852 elseif (obc_u(i,j) < 0) then ! This is a western OBC face.
853 h_src(:) = h(i+1,j,:)
854 h_dest(:) = remap_cs%h(i+1,j,:)
855 else ! (OBC_u(I,j) > 0) ! This is a eastern OBC face.
856 h_src(:) = h(i,j,:)
857 h_dest(:) = remap_cs%h(i,j,:)
858 endif
859 call interpolate_column(nz_src, h_src, field(i,j,:), &
860 nz_dest, h_dest, interpolated_field(i,j,:), .true.)
861 endif ; enddo ; enddo
862 else
863 do j=g%jsc,g%jec ; do i=g%IscB,g%IecB
864 if (obc_u(i,j) == 0) then ! This is not an OBC face.
865 h_src(:) = 0.5*(h(i,j,:) + h(i+1,j,:))
866 h_dest(:) = 0.5*(remap_cs%h(i,j,:) + remap_cs%h(i+1,j,:))
867 elseif (obc_u(i,j) < 0) then ! This is a western OBC face.
868 h_src(:) = h(i+1,j,:)
869 h_dest(:) = remap_cs%h(i+1,j,:)
870 else ! (OBC_u(I,j) > 0) ! This is a eastern OBC face.
871 h_src(:) = h(i,j,:)
872 h_dest(:) = remap_cs%h(i,j,:)
873 endif
874 call interpolate_column(nz_src, h_src, field(i,j,:), &
875 nz_dest, h_dest, interpolated_field(i,j,:), .true.)
876 enddo ; enddo
877 endif
878 elseif (staggered_in_y .and. .not. staggered_in_x) then
879 ! V-points
880 if (present(mask)) then
881 do j=g%jscB,g%jecB ; do i=g%isc,g%iec ; if (mask(i,j) > 0.0) then
882 if (obc_v(i,j) == 0) then ! This is not an OBC face.
883 h_src(:) = 0.5*(h(i,j,:) + h(i,j+1,:))
884 h_dest(:) = 0.5*(remap_cs%h(i,j,:) + remap_cs%h(i,j+1,:))
885 elseif (obc_v(i,j) < 0) then ! This is a southern OBC face
886 h_src(:) = h(i,j+1,:)
887 h_dest(:) = remap_cs%h(i,j+1,:)
888 else ! (OBC_v(i,J) > 0) ! This is a northern OBC face
889 h_src(:) = h(i,j,:)
890 h_dest(:) = remap_cs%h(i,j,:)
891 endif
892 call interpolate_column(nz_src, h_src, field(i,j,:), &
893 nz_dest, h_dest, interpolated_field(i,j,:), .true.)
894 endif ; enddo ; enddo
895 else
896 do j=g%jscB,g%jecB ; do i=g%isc,g%iec
897 if (obc_v(i,j) == 0) then ! This is not an OBC face.
898 h_src(:) = 0.5*(h(i,j,:) + h(i,j+1,:))
899 h_dest(:) = 0.5*(remap_cs%h(i,j,:) + remap_cs%h(i,j+1,:))
900 elseif (obc_v(i,j) < 0) then ! This is a southern OBC face
901 h_src(:) = h(i,j+1,:)
902 h_dest(:) = remap_cs%h(i,j+1,:)
903 else ! (OBC_v(i,J) > 0) ! This is a northern OBC face
904 h_src(:) = h(i,j,:)
905 h_dest(:) = remap_cs%h(i,j,:)
906 endif
907 call interpolate_column(nz_src, h_src, field(i,j,:), &
908 nz_dest, h_dest, interpolated_field(i,j,:), .true.)
909 enddo ; enddo
910 endif
911 elseif ((.not. staggered_in_x) .and. (.not. staggered_in_y)) then
912 ! H-points
913 if (present(mask)) then
914 do j=g%jsc,g%jec ; do i=g%isc,g%iec ; if (mask(i,j) > 0.0) then
915 call interpolate_column(nz_src, h(i,j,:), field(i,j,:), &
916 nz_dest, remap_cs%h(i,j,:), interpolated_field(i,j,:), .true.)
917 endif ; enddo ; enddo
918 else
919 do j=g%jsc,g%jec ; do i=g%isc,g%iec
920 call interpolate_column(nz_src, h(i,j,:), field(i,j,:), &
921 nz_dest, remap_cs%h(i,j,:), interpolated_field(i,j,:), .true.)
922 enddo ; enddo
923 endif
924 else
925 call assert(.false., 'vertically_interpolate_diag_field: Q point remapping is not coded yet.')
926 endif
927
928end subroutine vertically_interpolate_field
929
930!> Horizontally average a diagnostic field
931subroutine horizontally_average_diag_field(G, GV, h, staggered_in_x, staggered_in_y, &
932 is_layer, is_extensive, &
933 field, averaged_field, &
934 averaged_mask)
935 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
936 type(verticalgrid_type), intent(in) :: gv !< The ocean vertical grid structure
937 real, dimension(:,:,:), intent(in) :: h !< The current thicknesses [H ~> m or kg m-2]
938 logical, intent(in) :: staggered_in_x !< True if the x-axis location is at u or q points
939 logical, intent(in) :: staggered_in_y !< True if the y-axis location is at v or q points
940 logical, intent(in) :: is_layer !< True if the z-axis location is at h points
941 logical, intent(in) :: is_extensive !< True if the z-direction is spatially integrated (over layers)
942 real, dimension(:,:,:), intent(in) :: field !< The diagnostic field to be remapped [A]
943 real, dimension(:), intent(out) :: averaged_field !< Field argument horizontally averaged [A]
944 logical, dimension(:), intent(out) :: averaged_mask !< Mask for horizontally averaged field [nondim]
945
946 ! Local variables
947 integer :: isdf, jsdf !< The starting i- and j-indices in memory for field
948
949 isdf = g%isd ; if (staggered_in_x) isdf = g%IsdB
950 jsdf = g%jsd ; if (staggered_in_y) jsdf = g%JsdB
951
952 call horizontally_average_field(g, gv, isdf, jsdf, h, staggered_in_x, staggered_in_y, &
953 is_layer, is_extensive, field, averaged_field, averaged_mask)
954
956
957!> Horizontally average a diagnostic field
958subroutine horizontally_average_field(G, GV, isdf, jsdf, h, staggered_in_x, staggered_in_y, &
959 is_layer, is_extensive, field, averaged_field, averaged_mask)
960 type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
961 type(verticalgrid_type), intent(in) :: GV !< The ocean vertical grid structure
962 integer, intent(in) :: isdf !< The starting i-index in memory for field
963 integer, intent(in) :: jsdf !< The starting j-index in memory for field
964 real, dimension(G%isd:,G%jsd:,:), &
965 intent(in) :: h !< The current thicknesses [H ~> m or kg m-2]
966 logical, intent(in) :: staggered_in_x !< True if the x-axis location is at u or q points
967 logical, intent(in) :: staggered_in_y !< True if the y-axis location is at v or q points
968 logical, intent(in) :: is_layer !< True if the z-axis location is at h points
969 logical, intent(in) :: is_extensive !< True if the z-direction is spatially integrated (over layers)
970 real, dimension(isdf:,jsdf:,:), &
971 intent(in) :: field !< The diagnostic field to be remapped [A]
972 real, dimension(:), intent(out) :: averaged_field !< Field argument horizontally averaged [A]
973 logical, dimension(:), intent(out) :: averaged_mask !< Mask for horizontally averaged field [nondim]
974
975 ! Local variables
976 real :: volume(G%isc:G%iec, G%jsc:G%jec, size(field,3)) ! The area [L2 ~> m2], volume [L2 m ~> m3]
977 ! or mass [L2 kg m-2 ~> kg] of each cell.
978 real :: stuff(G%isc:G%iec, G%jsc:G%jec, size(field,3)) ! The area, volume or mass-weighted integral of the
979 ! field being averaged in each cell, in [L2 a ~> m2 A],
980 ! [L2 m a ~> m3 A] or [L2 kg m-2 A ~> kg A],
981 ! depending on the weighting for the averages and whether the
982 ! model makes the Boussinesq approximation.
983 real, dimension(size(field, 3)) :: vol_sum ! The global sum of the areas [m2], volumes [m3] or mass [kg]
984 ! in the cells that used in the weighted averages.
985 real, dimension(size(field, 3)) :: stuff_sum ! The global sum of the weighted field in all cells, in
986 ! [A m2], [A m3] or [A kg]
987 type(efp_type), dimension(2*size(field,3)) :: sums_EFP ! Sums of volume or stuff by layer
988 real :: height ! An average thickness attributed to an velocity point [H ~> m or kg m-2]
989 integer :: i, j, k, nz
990
991 nz = size(field, 3)
992
993 ! TODO: These averages could potentially be modified to use the function in
994 ! the MOM_spatial_means module.
995 ! NOTE: Reproducible sums must be computed in the original MKS units
996
997 if (staggered_in_x .and. .not. staggered_in_y) then
998 if (is_layer) then
999 ! U-points
1000 do k=1,nz
1001 vol_sum(k) = 0.
1002 stuff_sum(k) = 0.
1003 if (is_extensive) then
1004 do j=g%jsc, g%jec ; do i=g%isc, g%iec
1005 volume(i,j,k) = g%areaCu(i,j) * g%mask2dCu(i,j)
1006 stuff(i,j,k) = volume(i,j,k) * field(i,j,k)
1007 enddo ; enddo
1008 else ! Intensive
1009 do j=g%jsc, g%jec ; do i=g%isc, g%iec
1010 height = 0.5 * (h(i,j,k) + h(i+1,j,k))
1011 volume(i,j,k) = g%areaCu(i,j) * (gv%H_to_MKS * height) * g%mask2dCu(i,j)
1012 stuff(i,j,k) = volume(i,j,k) * field(i,j,k)
1013 enddo ; enddo
1014 endif
1015 enddo
1016 else ! Interface
1017 do k=1,nz
1018 do j=g%jsc, g%jec ; do i=g%isc, g%iec
1019 volume(i,j,k) = g%areaCu(i,j) * g%mask2dCu(i,j)
1020 stuff(i,j,k) = volume(i,j,k) * field(i,j,k)
1021 enddo ; enddo
1022 enddo
1023 endif
1024 elseif (staggered_in_y .and. .not. staggered_in_x) then
1025 if (is_layer) then
1026 ! V-points
1027 do k=1,nz
1028 if (is_extensive) then
1029 do j=g%jsc, g%jec ; do i=g%isc, g%iec
1030 volume(i,j,k) = g%areaCv(i,j) * g%mask2dCv(i,j)
1031 stuff(i,j,k) = volume(i,j,k) * field(i,j,k)
1032 enddo ; enddo
1033 else ! Intensive
1034 do j=g%jsc, g%jec ; do i=g%isc, g%iec
1035 height = 0.5 * (h(i,j,k) + h(i,j+1,k))
1036 volume(i,j,k) = g%areaCv(i,j) * (gv%H_to_MKS * height) * g%mask2dCv(i,j)
1037 stuff(i,j,k) = volume(i,j,k) * field(i,j,k)
1038 enddo ; enddo
1039 endif
1040 enddo
1041 else ! Interface
1042 do k=1,nz
1043 do j=g%jsc, g%jec ; do i=g%isc, g%iec
1044 volume(i,j,k) = g%areaCv(i,j) * g%mask2dCv(i,j)
1045 stuff(i,j,k) = volume(i,j,k) * field(i,j,k)
1046 enddo ; enddo
1047 enddo
1048 endif
1049 elseif ((.not. staggered_in_x) .and. (.not. staggered_in_y)) then
1050 if (is_layer) then
1051 ! H-points
1052 do k=1,nz
1053 if (is_extensive) then
1054 do j=g%jsc, g%jec ; do i=g%isc, g%iec
1055 if (h(i,j,k) > 0.) then
1056 volume(i,j,k) = g%areaT(i,j) * g%mask2dT(i,j)
1057 stuff(i,j,k) = volume(i,j,k) * field(i,j,k)
1058 else
1059 volume(i,j,k) = 0.
1060 stuff(i,j,k) = 0.
1061 endif
1062 enddo ; enddo
1063 else ! Intensive
1064 do j=g%jsc, g%jec ; do i=g%isc, g%iec
1065 volume(i,j,k) = g%areaT(i,j) * (gv%H_to_MKS * h(i,j,k)) * g%mask2dT(i,j)
1066 stuff(i,j,k) = volume(i,j,k) * field(i,j,k)
1067 enddo ; enddo
1068 endif
1069 enddo
1070 else ! Interface
1071 do k=1,nz
1072 do j=g%jsc, g%jec ; do i=g%isc, g%iec
1073 volume(i,j,k) = g%areaT(i,j) * g%mask2dT(i,j)
1074 stuff(i,j,k) = volume(i,j,k) * field(i,j,k)
1075 enddo ; enddo
1076 enddo
1077 endif
1078 else
1079 call assert(.false., 'horizontally_average_diag_field: Q point averaging is not coded yet.')
1080 endif
1081
1082 ! Packing the sums into a single array with a single call to sum across PEs saves reduces
1083 ! the costs of communication.
1084 do k=1,nz
1085 sums_efp(2*k-1) = reproducing_sum_efp(volume(:,:,k), only_on_pe=.true., unscale=g%US%L_to_m**2)
1086 sums_efp(2*k) = reproducing_sum_efp(stuff(:,:,k), only_on_pe=.true., unscale=g%US%L_to_m**2)
1087 enddo
1088 call efp_sum_across_pes(sums_efp, 2*nz)
1089 do k=1,nz
1090 vol_sum(k) = efp_to_real(sums_efp(2*k-1))
1091 stuff_sum(k) = efp_to_real(sums_efp(2*k))
1092 enddo
1093
1094 averaged_mask(:) = .true.
1095 do k=1,nz
1096 if (vol_sum(k) > 0.) then
1097 averaged_field(k) = stuff_sum(k) / vol_sum(k)
1098 else
1099 averaged_field(k) = 0.
1100 averaged_mask(k) = .false.
1101 endif
1102 enddo
1103
1104end subroutine horizontally_average_field
1105
1106end module mom_diag_remap