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, 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 logical, intent(in) :: staggered_in_x !< True is the x-axis location is at u or q points
390 logical, intent(in) :: staggered_in_y !< True is the y-axis location is at v or q points
391 real, dimension(:,:,:), pointer :: mask !< A mask for the field [nondim].
392 real, dimension(:,:,:), intent(in) :: field(:,:,:) !< The diagnostic field to be remapped [A]
393 real, dimension(:,:,:), intent(out) :: remapped_field !< Field remapped to new coordinate [A]
394
395 ! Local variables
396 integer :: isdf, jsdf !< The starting i- and j-indices in memory for field
397
398 call assert(remap_cs%initialized, 'diag_remap_do_remap: remap_cs not initialized.')
399 call assert(size(field, 3) == size(h, 3), &
400 'diag_remap_do_remap: Remap field and thickness z-axes do not match.')
401
402 isdf = g%isd ; if (staggered_in_x) isdf = g%IsdB
403 jsdf = g%jsd ; if (staggered_in_y) jsdf = g%JsdB
404
405 if (associated(mask)) then
406 call do_remap(remap_cs, g, gv, us, isdf, jsdf, h, staggered_in_x, staggered_in_y, &
407 field, remapped_field, mask(:,:,1))
408 else
409 call do_remap(remap_cs, g, gv, us, isdf, jsdf, h, staggered_in_x, staggered_in_y, &
410 field, remapped_field)
411 endif
412
413end subroutine diag_remap_do_remap
414
415!> The internal routine to remap a diagnostic field to an alternative vertical grid.
416subroutine do_remap(remap_cs, G, GV, US, isdf, jsdf, h, staggered_in_x, staggered_in_y, &
417 field, remapped_field, mask)
418 type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coordinate control structure
419 type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
420 type(verticalgrid_type), intent(in) :: GV !< ocean vertical grid structure
421 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
422 integer, intent(in) :: isdf !< The starting i-index in memory for field
423 integer, intent(in) :: jsdf !< The starting j-index in memory for field
424 real, dimension(G%isd:,G%jsd:,:), &
425 intent(in) :: h !< The current thicknesses [H ~> m or kg m-2] or [Z ~> m],
426 !! depending on the value of remap_CS%Z_based_coord
427 logical, intent(in) :: staggered_in_x !< True is the x-axis location is at u or q points
428 logical, intent(in) :: staggered_in_y !< True is the y-axis location is at v or q points
429 real, dimension(isdf:,jsdf:,:), &
430 intent(in) :: field !< The diagnostic field to be remapped [A]
431 real, dimension(isdf:,jsdf:,:), &
432 intent(out) :: remapped_field !< Field remapped to new coordinate [A]
433 real, dimension(isdf:,jsdf:), &
434 optional, intent(in) :: mask !< A mask for the field [nondim]
435
436 ! Local variables
437 real, dimension(remap_cs%nz) :: h_dest ! Destination thicknesses [H ~> m or kg m-2] or [Z ~> m]
438 real, dimension(size(h,3)) :: h_src ! A column of source thicknesses [H ~> m or kg m-2] or [Z ~> m]
439 integer :: nz_src, nz_dest ! The number of layers on the native and remapped grids
440 integer :: i, j ! Grid index
441
442 nz_src = size(field,3)
443 nz_dest = remap_cs%nz
444 remapped_field(:,:,:) = 0.
445
446 if (staggered_in_x .and. .not. staggered_in_y) then
447 ! U-points
448 if (present(mask)) then
449 do j=g%jsc,g%jec ; do i=g%IscB,g%IecB ; if (mask(i,j) > 0.) then
450 h_src(:) = 0.5 * (h(i,j,:) + h(i+1,j,:))
451 h_dest(:) = 0.5 * (remap_cs%h(i,j,:) + remap_cs%h(i+1,j,:))
452 call remapping_core_h(remap_cs%remap_cs, nz_src, h_src(:), field(i,j,:), &
453 nz_dest, h_dest(:), remapped_field(i,j,:))
454 endif ; enddo ; enddo
455 else
456 do j=g%jsc,g%jec ; do i=g%IscB,g%IecB
457 h_src(:) = 0.5 * (h(i,j,:) + h(i+1,j,:))
458 h_dest(:) = 0.5 * (remap_cs%h(i,j,:) + remap_cs%h(i+1,j,:))
459 call remapping_core_h(remap_cs%remap_cs, nz_src, h_src(:), field(i,j,:), &
460 nz_dest, h_dest(:), remapped_field(i,j,:))
461 enddo ; enddo
462 endif
463 elseif (staggered_in_y .and. .not. staggered_in_x) then
464 ! V-points
465 if (present(mask)) then
466 do j=g%jscB,g%jecB ; do i=g%isc,g%iec ; if (mask(i,j) > 0.) then
467 h_src(:) = 0.5 * (h(i,j,:) + h(i,j+1,:))
468 h_dest(:) = 0.5 * (remap_cs%h(i,j,:) + remap_cs%h(i,j+1,:))
469 call remapping_core_h(remap_cs%remap_cs, nz_src, h_src(:), field(i,j,:), &
470 nz_dest, h_dest(:), remapped_field(i,j,:))
471 endif ; enddo ; enddo
472 else
473 do j=g%jscB,g%jecB ; do i=g%isc,g%iec
474 h_src(:) = 0.5 * (h(i,j,:) + h(i,j+1,:))
475 h_dest(:) = 0.5 * (remap_cs%h(i,j,:) + remap_cs%h(i,j+1,:))
476 call remapping_core_h(remap_cs%remap_cs, nz_src, h_src(:), field(i,j,:), &
477 nz_dest, h_dest(:), remapped_field(i,j,:))
478 enddo ; enddo
479 endif
480 elseif ((.not. staggered_in_x) .and. (.not. staggered_in_y)) then
481 ! H-points
482 if (present(mask)) then
483 do j=g%jsc,g%jec ; do i=g%isc,g%iec ; if (mask(i,j) > 0.) then
484 call remapping_core_h(remap_cs%remap_cs, nz_src, h(i,j,:), field(i,j,:), &
485 nz_dest, remap_cs%h(i,j,:), remapped_field(i,j,:))
486 endif ; enddo ; enddo
487 else
488 do j=g%jsc,g%jec ; do i=g%isc,g%iec
489 call remapping_core_h(remap_cs%remap_cs, nz_src, h(i,j,:), field(i,j,:), &
490 nz_dest, remap_cs%h(i,j,:), remapped_field(i,j,:))
491 enddo ; enddo
492 endif
493 else
494 call assert(.false., 'diag_remap_do_remap: Unsupported axis combination')
495 endif
496
497end subroutine do_remap
498
499!> Calculate masks for target grid
500subroutine diag_remap_calc_hmask(remap_cs, G, mask)
501 type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coordinate control structure
502 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
503 real, dimension(G%isd:,G%jsd:,:), &
504 intent(out) :: mask !< h-point mask for target grid [nondim]
505
506 ! Local variables
507 real, dimension(remap_cs%nz) :: h_dest ! Destination thicknesses [H ~> m or kg m-2] or [Z ~> m]
508 integer :: i, j, k
509 logical :: mask_vanished_layers
510 real :: h_tot ! Sum of all thicknesses [H ~> m or kg m-2] or [Z ~> m]
511 real :: h_err ! An estimate of a negligible thickness [H ~> m or kg m-2] or [Z ~> m]
512
513 call assert(remap_cs%initialized, 'diag_remap_calc_hmask: remap_cs not initialized.')
514
515 ! Only z*-like diagnostic coordinates should have a 3d mask
516 mask_vanished_layers = (remap_cs%vertical_coord == coordinatemode('ZSTAR'))
517 mask(:,:,:) = 0.
518
519 do j=g%jsc-1,g%jec+1 ; do i=g%isc-1,g%iec+1
520 if (g%mask2dT(i,j)>0.) then
521 if (mask_vanished_layers) then
522 h_dest(:) = remap_cs%h(i,j,:)
523 h_tot = 0.
524 h_err = 0.
525 do k=1, remap_cs%nz
526 h_tot = h_tot + h_dest(k)
527 ! This is an overestimate of how thick a vanished layer might be, that
528 ! appears due to round-off.
529 h_err = h_err + epsilon(h_tot) * h_tot
530 ! Mask out vanished layers
531 if (h_dest(k)<=8.*h_err) then
532 mask(i,j,k) = 0.
533 else
534 mask(i,j,k) = 1.
535 endif
536 enddo
537 else ! all layers might contain data
538 mask(i,j,:) = 1.
539 endif
540 endif
541 enddo ; enddo
542
543end subroutine diag_remap_calc_hmask
544
545!> Vertically re-grid an already vertically-integrated diagnostic field to alternative vertical grid.
546subroutine vertically_reintegrate_diag_field(remap_cs, G, h, h_target, staggered_in_x, staggered_in_y, &
547 mask, field, reintegrated_field)
548 type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coordinate control structure
549 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
550 real, dimension(:,:,:), intent(in) :: h !< The thicknesses of the source grid [H ~> m or kg m-2] or [Z ~> m]
551 real, dimension(:,:,:), intent(in) :: h_target !< The thicknesses of the target grid [H ~> m or kg m-2] or [Z ~> m]
552 logical, intent(in) :: staggered_in_x !< True is the x-axis location is at u or q points
553 logical, intent(in) :: staggered_in_y !< True is the y-axis location is at v or q points
554 real, dimension(:,:,:), pointer :: mask !< A mask for the field [nondim]. Note that because this
555 !! is a pointer it retains its declared indexing conventions.
556 real, dimension(:,:,:), intent(in) :: field !< The diagnostic field to be remapped [A]
557 real, dimension(:,:,:), intent(out) :: reintegrated_field !< Field argument remapped to alternative coordinate [A]
558
559 ! Local variables
560 integer :: isdf, jsdf !< The starting i- and j-indices in memory for field
561
562 call assert(remap_cs%initialized, 'vertically_reintegrate_diag_field: remap_cs not initialized.')
563 call assert(size(field, 3) == size(h, 3), &
564 'vertically_reintegrate_diag_field: Remap field and thickness z-axes do not match.')
565
566 isdf = g%isd ; if (staggered_in_x) isdf = g%IsdB
567 jsdf = g%jsd ; if (staggered_in_y) jsdf = g%JsdB
568
569 if (associated(mask)) then
570 call vertically_reintegrate_field(remap_cs, g, isdf, jsdf, h, h_target, staggered_in_x, staggered_in_y, &
571 field, reintegrated_field, mask(:,:,1))
572 else
573 call vertically_reintegrate_field(remap_cs, g, isdf, jsdf, h, h_target, staggered_in_x, staggered_in_y, &
574 field, reintegrated_field)
575 endif
576
578
579!> The internal routine to vertically re-grid an already vertically-integrated diagnostic field to
580!! an alternative vertical grid.
581subroutine vertically_reintegrate_field(remap_cs, G, isdf, jsdf, h, h_target, staggered_in_x, staggered_in_y, &
582 field, reintegrated_field, mask)
583 type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coordinate control structure
584 type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
585 integer, intent(in) :: isdf !< The starting i-index in memory for field
586 integer, intent(in) :: jsdf !< The starting j-index in memory for field
587 real, dimension(G%isd:,G%jsd:,:), &
588 intent(in) :: h !< The thicknesses of the source grid [H ~> m or kg m-2] or [Z ~> m]
589 real, dimension(G%isd:,G%jsd:,:), &
590 intent(in) :: h_target !< The thicknesses of the target grid [H ~> m or kg m-2] or [Z ~> m]
591 logical, intent(in) :: staggered_in_x !< True is the x-axis location is at u or q points
592 logical, intent(in) :: staggered_in_y !< True is the y-axis location is at v or q points
593 real, dimension(isdf:,jsdf:,:), &
594 intent(in) :: field !< The diagnostic field to be remapped [A]
595 real, dimension(isdf:,jsdf:,:), &
596 intent(out) :: reintegrated_field !< Field argument remapped to alternative coordinate [A]
597 real, dimension(isdf:,jsdf:), &
598 optional, intent(in) :: mask !< A mask for the field [nondim]
599
600 ! Local variables
601 real, dimension(remap_cs%nz) :: h_dest ! Destination thicknesses [H ~> m or kg m-2] or [Z ~> m]
602 real, dimension(size(h,3)) :: h_src ! A column of source thicknesses [H ~> m or kg m-2] or [Z ~> m]
603 integer :: nz_src, nz_dest ! The number of layers on the native and remapped grids
604 integer :: i, j ! Grid index
605
606 nz_src = size(field,3)
607 nz_dest = remap_cs%nz
608 reintegrated_field(:,:,:) = 0.
609
610 if (staggered_in_x .and. .not. staggered_in_y) then
611 ! U-points
612 if (present(mask)) then
613 do j=g%jsc,g%jec ; do i=g%IscB,g%IecB ; if (mask(i,j) > 0.0) then
614 h_src(:) = 0.5 * (h(i,j,:) + h(i+1,j,:))
615 h_dest(:) = 0.5 * (h_target(i,j,:) + h_target(i+1,j,:))
616 call reintegrate_column(nz_src, h_src, field(i,j,:), &
617 nz_dest, h_dest, reintegrated_field(i,j,:))
618 endif ; enddo ; enddo
619 else
620 do j=g%jsc,g%jec ; do i=g%IscB,g%IecB
621 h_src(:) = 0.5 * (h(i,j,:) + h(i+1,j,:))
622 h_dest(:) = 0.5 * (h_target(i,j,:) + h_target(i+1,j,:))
623 call reintegrate_column(nz_src, h_src, field(i,j,:), &
624 nz_dest, h_dest, reintegrated_field(i,j,:))
625 enddo ; enddo
626 endif
627 elseif (staggered_in_y .and. .not. staggered_in_x) then
628 ! V-points
629 if (present(mask)) then
630 do j=g%jscB,g%jecB ; do i=g%isc,g%iec ; if (mask(i,j) > 0.0) then
631 h_src(:) = 0.5 * (h(i,j,:) + h(i,j+1,:))
632 h_dest(:) = 0.5 * (h_target(i,j,:) + h_target(i,j+1,:))
633 call reintegrate_column(nz_src, h_src, field(i,j,:), &
634 nz_dest, h_dest, reintegrated_field(i,j,:))
635 endif ; enddo ; enddo
636 else
637 do j=g%jscB,g%jecB ; do i=g%isc,g%iec
638 h_src(:) = 0.5 * (h(i,j,:) + h(i,j+1,:))
639 h_dest(:) = 0.5 * (h_target(i,j,:) + h_target(i,j+1,:))
640 call reintegrate_column(nz_src, h_src, field(i,j,:), &
641 nz_dest, h_dest, reintegrated_field(i,j,:))
642 enddo ; enddo
643 endif
644 elseif ((.not. staggered_in_x) .and. (.not. staggered_in_y)) then
645 ! H-points
646 if (present(mask)) then
647 do j=g%jsc,g%jec ; do i=g%isc,g%iec ; if (mask(i,j) > 0.0) then
648 call reintegrate_column(nz_src, h(i,j,:), field(i,j,:), &
649 nz_dest, h_target(i,j,:), reintegrated_field(i,j,:))
650 endif ; enddo ; enddo
651 else
652 do j=g%jsc,g%jec ; do i=g%isc,g%iec
653 call reintegrate_column(nz_src, h(i,j,:), field(i,j,:), &
654 nz_dest, h_target(i,j,:), reintegrated_field(i,j,:))
655 enddo ; enddo
656 endif
657 else
658 call assert(.false., 'vertically_reintegrate_diag_field: Q point remapping is not coded yet.')
659 endif
660
661end subroutine vertically_reintegrate_field
662
663!> Vertically interpolate diagnostic field to alternative vertical grid.
664subroutine vertically_interpolate_diag_field(remap_cs, G, h, staggered_in_x, staggered_in_y, &
665 mask, field, interpolated_field)
666 type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coordinate control structure
667 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
668 real, dimension(:,:,:), intent(in) :: h !< The current thicknesses [H ~> m or kg m-2] or [Z ~> m],
669 !! depending on the value of remap_cs%Z_based_coord
670 logical, intent(in) :: staggered_in_x !< True is the x-axis location is at u or q points
671 logical, intent(in) :: staggered_in_y !< True is the y-axis location is at v or q points
672 real, dimension(:,:,:), pointer :: mask !< A mask for the field [nondim]. Note that because this
673 !! is a pointer it retains its declared indexing conventions.
674 real, dimension(:,:,:), intent(in) :: field !< The diagnostic field to be remapped [A]
675 real, dimension(:,:,:), intent(inout) :: interpolated_field !< Field argument remapped to alternative coordinate [A]
676
677 ! Local variables
678 integer :: isdf, jsdf !< The starting i- and j-indices in memory for field
679
680 call assert(remap_cs%initialized, 'vertically_interpolate_diag_field: remap_cs not initialized.')
681 call assert(size(field, 3) == size(h, 3)+1, &
682 'vertically_interpolate_diag_field: Remap field and thickness z-axes do not match.')
683
684 isdf = g%isd ; if (staggered_in_x) isdf = g%IsdB
685 jsdf = g%jsd ; if (staggered_in_y) jsdf = g%JsdB
686
687 if (associated(mask)) then
688 call vertically_interpolate_field(remap_cs, g, isdf, jsdf, h, staggered_in_x, staggered_in_y, &
689 field, interpolated_field, mask(:,:,1))
690 else
691 call vertically_interpolate_field(remap_cs, g, isdf, jsdf, h, staggered_in_x, staggered_in_y, &
692 field, interpolated_field)
693 endif
694
696
697!> Internal routine to vertically interpolate a diagnostic field to an alternative vertical grid.
698subroutine vertically_interpolate_field(remap_cs, G, isdf, jsdf, h, staggered_in_x, staggered_in_y, &
699 field, interpolated_field, mask)
700 type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coordinate control structure
701 type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
702 integer, intent(in) :: isdf !< The starting i-index in memory for field
703 integer, intent(in) :: jsdf !< The starting j-index in memory for field
704 real, dimension(G%isd:,G%jsd:,:), &
705 intent(in) :: h !< The current thicknesses [H ~> m or kg m-2] or [Z ~> m],
706 !! depending on the value of remap_cs%Z_based_coord
707 logical, intent(in) :: staggered_in_x !< True is the x-axis location is at u or q points
708 logical, intent(in) :: staggered_in_y !< True is the y-axis location is at v or q points
709 real, dimension(isdf:,jsdf:,:), &
710 intent(in) :: field !< The diagnostic field to be remapped [A]
711 real, dimension(isdf:,jsdf:,:), &
712 intent(out) :: interpolated_field !< Field argument remapped to alternative coordinate [A]
713 real, dimension(isdf:,jsdf:), &
714 optional, intent(in) :: mask !< A mask for the field [nondim]
715
716 ! Local variables
717 real, dimension(remap_cs%nz) :: h_dest ! Destination thicknesses [H ~> m or kg m-2] or [Z ~> m]
718 real, dimension(size(h,3)) :: h_src ! A column of source thicknesses [H ~> m or kg m-2] or [Z ~> m]
719 integer :: nz_src, nz_dest ! The number of layers on the native and remapped grids
720 integer :: i, j !< Grid index
721
722 interpolated_field(:,:,:) = 0.
723
724 nz_src = size(h,3)
725 nz_dest = remap_cs%nz
726
727 if (staggered_in_x .and. .not. staggered_in_y) then
728 ! U-points
729 if (present(mask)) then
730 do j=g%jsc,g%jec ; do i=g%IscB,g%IecB ; if (mask(i,j) > 0.0) then
731 h_src(:) = 0.5 * (h(i,j,:) + h(i+1,j,:))
732 h_dest(:) = 0.5 * (remap_cs%h(i,j,:) + remap_cs%h(i+1,j,:))
733 call interpolate_column(nz_src, h_src, field(i,j,:), &
734 nz_dest, h_dest, interpolated_field(i,j,:), .true.)
735 endif ; enddo ; enddo
736 else
737 do j=g%jsc,g%jec ; do i=g%IscB,g%IecB
738 h_src(:) = 0.5 * (h(i,j,:) + h(i+1,j,:))
739 h_dest(:) = 0.5 * (remap_cs%h(i,j,:) + remap_cs%h(i+1,j,:))
740 call interpolate_column(nz_src, h_src, field(i,j,:), &
741 nz_dest, h_dest, interpolated_field(i,j,:), .true.)
742 enddo ; enddo
743 endif
744 elseif (staggered_in_y .and. .not. staggered_in_x) then
745 ! V-points
746 if (present(mask)) then
747 do j=g%jscB,g%jecB ; do i=g%isc,g%iec ; if (mask(i,j) > 0.0) then
748 h_src(:) = 0.5 * (h(i,j,:) + h(i,j+1,:))
749 h_dest(:) = 0.5 * (remap_cs%h(i,j,:) + remap_cs%h(i,j+1,:))
750 call interpolate_column(nz_src, h_src, field(i,j,:), &
751 nz_dest, h_dest, interpolated_field(i,j,:), .true.)
752 endif ; enddo ; enddo
753 else
754 do j=g%jscB,g%jecB ; do i=g%isc,g%iec
755 h_src(:) = 0.5 * (h(i,j,:) + h(i,j+1,:))
756 h_dest(:) = 0.5 * (remap_cs%h(i,j,:) + remap_cs%h(i,j+1,:))
757 call interpolate_column(nz_src, h_src, field(i,j,:), &
758 nz_dest, h_dest, interpolated_field(i,j,:), .true.)
759 enddo ; enddo
760 endif
761 elseif ((.not. staggered_in_x) .and. (.not. staggered_in_y)) then
762 ! H-points
763 if (present(mask)) then
764 do j=g%jsc,g%jec ; do i=g%isc,g%iec ; if (mask(i,j) > 0.0) then
765 call interpolate_column(nz_src, h(i,j,:), field(i,j,:), &
766 nz_dest, remap_cs%h(i,j,:), interpolated_field(i,j,:), .true.)
767 endif ; enddo ; enddo
768 else
769 do j=g%jsc,g%jec ; do i=g%isc,g%iec
770 call interpolate_column(nz_src, h(i,j,:), field(i,j,:), &
771 nz_dest, remap_cs%h(i,j,:), interpolated_field(i,j,:), .true.)
772 enddo ; enddo
773 endif
774 else
775 call assert(.false., 'vertically_interpolate_diag_field: Q point remapping is not coded yet.')
776 endif
777
778end subroutine vertically_interpolate_field
779
780!> Horizontally average a diagnostic field
781subroutine horizontally_average_diag_field(G, GV, h, staggered_in_x, staggered_in_y, &
782 is_layer, is_extensive, &
783 field, averaged_field, &
784 averaged_mask)
785 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
786 type(verticalgrid_type), intent(in) :: gv !< The ocean vertical grid structure
787 real, dimension(:,:,:), intent(in) :: h !< The current thicknesses [H ~> m or kg m-2]
788 logical, intent(in) :: staggered_in_x !< True if the x-axis location is at u or q points
789 logical, intent(in) :: staggered_in_y !< True if the y-axis location is at v or q points
790 logical, intent(in) :: is_layer !< True if the z-axis location is at h points
791 logical, intent(in) :: is_extensive !< True if the z-direction is spatially integrated (over layers)
792 real, dimension(:,:,:), intent(in) :: field !< The diagnostic field to be remapped [A]
793 real, dimension(:), intent(out) :: averaged_field !< Field argument horizontally averaged [A]
794 logical, dimension(:), intent(out) :: averaged_mask !< Mask for horizontally averaged field [nondim]
795
796 ! Local variables
797 integer :: isdf, jsdf !< The starting i- and j-indices in memory for field
798
799 isdf = g%isd ; if (staggered_in_x) isdf = g%IsdB
800 jsdf = g%jsd ; if (staggered_in_y) jsdf = g%JsdB
801
802 call horizontally_average_field(g, gv, isdf, jsdf, h, staggered_in_x, staggered_in_y, &
803 is_layer, is_extensive, field, averaged_field, averaged_mask)
804
806
807!> Horizontally average a diagnostic field
808subroutine horizontally_average_field(G, GV, isdf, jsdf, h, staggered_in_x, staggered_in_y, &
809 is_layer, is_extensive, field, averaged_field, averaged_mask)
810 type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
811 type(verticalgrid_type), intent(in) :: GV !< The ocean vertical grid structure
812 integer, intent(in) :: isdf !< The starting i-index in memory for field
813 integer, intent(in) :: jsdf !< The starting j-index in memory for field
814 real, dimension(G%isd:,G%jsd:,:), &
815 intent(in) :: h !< The current thicknesses [H ~> m or kg m-2]
816 logical, intent(in) :: staggered_in_x !< True if the x-axis location is at u or q points
817 logical, intent(in) :: staggered_in_y !< True if the y-axis location is at v or q points
818 logical, intent(in) :: is_layer !< True if the z-axis location is at h points
819 logical, intent(in) :: is_extensive !< True if the z-direction is spatially integrated (over layers)
820 real, dimension(isdf:,jsdf:,:), &
821 intent(in) :: field !< The diagnostic field to be remapped [A]
822 real, dimension(:), intent(out) :: averaged_field !< Field argument horizontally averaged [A]
823 logical, dimension(:), intent(out) :: averaged_mask !< Mask for horizontally averaged field [nondim]
824
825 ! Local variables
826 real :: volume(G%isc:G%iec, G%jsc:G%jec, size(field,3)) ! The area [L2 ~> m2], volume [L2 m ~> m3]
827 ! or mass [L2 kg m-2 ~> kg] of each cell.
828 real :: stuff(G%isc:G%iec, G%jsc:G%jec, size(field,3)) ! The area, volume or mass-weighted integral of the
829 ! field being averaged in each cell, in [L2 a ~> m2 A],
830 ! [L2 m a ~> m3 A] or [L2 kg m-2 A ~> kg A],
831 ! depending on the weighting for the averages and whether the
832 ! model makes the Boussinesq approximation.
833 real, dimension(size(field, 3)) :: vol_sum ! The global sum of the areas [m2], volumes [m3] or mass [kg]
834 ! in the cells that used in the weighted averages.
835 real, dimension(size(field, 3)) :: stuff_sum ! The global sum of the weighted field in all cells, in
836 ! [A m2], [A m3] or [A kg]
837 type(efp_type), dimension(2*size(field,3)) :: sums_EFP ! Sums of volume or stuff by layer
838 real :: height ! An average thickness attributed to an velocity point [H ~> m or kg m-2]
839 integer :: i, j, k, nz
840
841 nz = size(field, 3)
842
843 ! TODO: These averages could potentially be modified to use the function in
844 ! the MOM_spatial_means module.
845 ! NOTE: Reproducible sums must be computed in the original MKS units
846
847 if (staggered_in_x .and. .not. staggered_in_y) then
848 if (is_layer) then
849 ! U-points
850 do k=1,nz
851 vol_sum(k) = 0.
852 stuff_sum(k) = 0.
853 if (is_extensive) then
854 do j=g%jsc, g%jec ; do i=g%isc, g%iec
855 volume(i,j,k) = g%areaCu(i,j) * g%mask2dCu(i,j)
856 stuff(i,j,k) = volume(i,j,k) * field(i,j,k)
857 enddo ; enddo
858 else ! Intensive
859 do j=g%jsc, g%jec ; do i=g%isc, g%iec
860 height = 0.5 * (h(i,j,k) + h(i+1,j,k))
861 volume(i,j,k) = g%areaCu(i,j) * (gv%H_to_MKS * height) * g%mask2dCu(i,j)
862 stuff(i,j,k) = volume(i,j,k) * field(i,j,k)
863 enddo ; enddo
864 endif
865 enddo
866 else ! Interface
867 do k=1,nz
868 do j=g%jsc, g%jec ; do i=g%isc, g%iec
869 volume(i,j,k) = g%areaCu(i,j) * g%mask2dCu(i,j)
870 stuff(i,j,k) = volume(i,j,k) * field(i,j,k)
871 enddo ; enddo
872 enddo
873 endif
874 elseif (staggered_in_y .and. .not. staggered_in_x) then
875 if (is_layer) then
876 ! V-points
877 do k=1,nz
878 if (is_extensive) then
879 do j=g%jsc, g%jec ; do i=g%isc, g%iec
880 volume(i,j,k) = g%areaCv(i,j) * g%mask2dCv(i,j)
881 stuff(i,j,k) = volume(i,j,k) * field(i,j,k)
882 enddo ; enddo
883 else ! Intensive
884 do j=g%jsc, g%jec ; do i=g%isc, g%iec
885 height = 0.5 * (h(i,j,k) + h(i,j+1,k))
886 volume(i,j,k) = g%areaCv(i,j) * (gv%H_to_MKS * height) * g%mask2dCv(i,j)
887 stuff(i,j,k) = volume(i,j,k) * field(i,j,k)
888 enddo ; enddo
889 endif
890 enddo
891 else ! Interface
892 do k=1,nz
893 do j=g%jsc, g%jec ; do i=g%isc, g%iec
894 volume(i,j,k) = g%areaCv(i,j) * g%mask2dCv(i,j)
895 stuff(i,j,k) = volume(i,j,k) * field(i,j,k)
896 enddo ; enddo
897 enddo
898 endif
899 elseif ((.not. staggered_in_x) .and. (.not. staggered_in_y)) then
900 if (is_layer) then
901 ! H-points
902 do k=1,nz
903 if (is_extensive) then
904 do j=g%jsc, g%jec ; do i=g%isc, g%iec
905 if (h(i,j,k) > 0.) then
906 volume(i,j,k) = g%areaT(i,j) * g%mask2dT(i,j)
907 stuff(i,j,k) = volume(i,j,k) * field(i,j,k)
908 else
909 volume(i,j,k) = 0.
910 stuff(i,j,k) = 0.
911 endif
912 enddo ; enddo
913 else ! Intensive
914 do j=g%jsc, g%jec ; do i=g%isc, g%iec
915 volume(i,j,k) = g%areaT(i,j) * (gv%H_to_MKS * h(i,j,k)) * g%mask2dT(i,j)
916 stuff(i,j,k) = volume(i,j,k) * field(i,j,k)
917 enddo ; enddo
918 endif
919 enddo
920 else ! Interface
921 do k=1,nz
922 do j=g%jsc, g%jec ; do i=g%isc, g%iec
923 volume(i,j,k) = g%areaT(i,j) * g%mask2dT(i,j)
924 stuff(i,j,k) = volume(i,j,k) * field(i,j,k)
925 enddo ; enddo
926 enddo
927 endif
928 else
929 call assert(.false., 'horizontally_average_diag_field: Q point averaging is not coded yet.')
930 endif
931
932 ! Packing the sums into a single array with a single call to sum across PEs saves reduces
933 ! the costs of communication.
934 do k=1,nz
935 sums_efp(2*k-1) = reproducing_sum_efp(volume(:,:,k), only_on_pe=.true., unscale=g%US%L_to_m**2)
936 sums_efp(2*k) = reproducing_sum_efp(stuff(:,:,k), only_on_pe=.true., unscale=g%US%L_to_m**2)
937 enddo
938 call efp_sum_across_pes(sums_efp, 2*nz)
939 do k=1,nz
940 vol_sum(k) = efp_to_real(sums_efp(2*k-1))
941 stuff_sum(k) = efp_to_real(sums_efp(2*k))
942 enddo
943
944 averaged_mask(:) = .true.
945 do k=1,nz
946 if (vol_sum(k) > 0.) then
947 averaged_field(k) = stuff_sum(k) / vol_sum(k)
948 else
949 averaged_field(k) = 0.
950 averaged_mask(k) = .false.
951 endif
952 enddo
953
954end subroutine horizontally_average_field
955
956end module mom_diag_remap