MOM_spatial_means.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!> Functions and routines to take area, volume, mass-weighted, layerwise, zonal or meridional means
7
8use mom_coms, only : efp_type, operator(+), operator(-), assignment(=)
9use mom_coms, only : efp_to_real, real_to_efp, efp_sum_across_pes
10use mom_coms, only : reproducing_sum, reproducing_sum_efp, efp_to_real
12use mom_coms, only : max_across_pes, min_across_pes
13use mom_error_handler, only : mom_error, note, warning, fatal, is_root_pe
14use mom_file_parser, only : get_param, log_version, param_file_type
15use mom_grid, only : ocean_grid_type
18
19implicit none ; private
20
21#include <MOM_memory.h>
22
29
30! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
31! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
32! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
33! vary with the Boussinesq approximation, the Boussinesq variant is given first.
34! The functions in this module work with variables with arbitrary units, in which case the
35! arbitrary rescaled units are indicated with [A ~> a], while the unscaled units are just [a].
36
37contains
38
39!> Return the global area mean of a variable, perhaps with a change of units. This uses reproducing sums.
40function global_area_mean(var, G, scale, tmp_scale, unscale)
41 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
42 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: var !< The variable to average in arbitrary units [a],
43 !! or arbitrary rescaled units [A ~> a] if unscale
44 !! or tmp_scale is present
45 !! arbitrary, possibly rescaled units [A ~> a]
46 real, optional, intent(in) :: scale !< A rescaling factor for the variable [a A-1 ~> 1]
47 !! that converts it back to unscaled (e.g., mks)
48 !! units to enable the use of the reproducing sums
49 real, optional, intent(in) :: tmp_scale !< A temporary rescaling factor for the variable
50 !! that is reversed in the return value [a A-1 ~> 1],
51 !! or [b B-1 ~> 1] if unscale is also present.
52 real, optional, intent(in) :: unscale !< A rescaling factor for the variable [a A-1 ~> 1]
53 !! that converts it back to unscaled (e.g., mks)
54 !! units to enable the use of the reproducing sums, or
55 !! a factor converting between rescaled units if
56 !! tmp_scale is also present [B A-1 ~> b a-1].
57 !! Here scale and unscale are synonymous, but unscale
58 !! is preferred and takes precedence.
59 real :: global_area_mean ! The mean of the variable in arbitrary unscaled units [a] or scaled units [A ~> a]
60 ! or [B ~> b], depending on which optional arguments are provided
61
62 ! Local variables
63 ! In the following comments, [A ~> a] is used to indicate the arbitrary, possibly rescaled units of the
64 ! input array while [a] indicates the unscaled (e.g., mks) units that can be used with the reproducing sums.
65 ! [A ~> a] and [B ~> b] are the same units unless tmp_scale and unscale are both present.
66 real :: tmpforsumming(szi_(g),szj_(g)) ! An unscaled cell integral in [a L2 ~> a m2] or a
67 ! scaled cell integral in [A L2 ~> a m2] or [B L2 ~> b m2]
68 real :: scalefac ! A scaling factor for the variable that is not reversed [a A-1 ~> 1] or [B A-1 ~> b a-1] or [1]
69 real :: temp_scale ! A temporary scaling factor [a A-1 ~> 1] or [b B-1 ~> 1] or [1]
70 integer :: i, j, is, ie, js, je
71 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
72
73 temp_scale = 1.0 ; if (present(tmp_scale)) temp_scale = tmp_scale
74
75 scalefac = 1.0
76 if (present(unscale)) then ; scalefac = unscale
77 elseif (present(scale)) then ; scalefac = scale ; endif
78
79 tmpforsumming(:,:) = 0.
80 do j=js,je ; do i=is,ie
81 tmpforsumming(i,j) = var(i,j) * (scalefac * g%areaT(i,j) * g%mask2dT(i,j))
82 enddo ; enddo
83
84 global_area_mean = reproducing_sum(tmpforsumming, unscale=temp_scale*g%US%L_to_m**2) * g%IareaT_global
85
86end function global_area_mean
87
88!> Return the global area mean of a variable. This uses reproducing sums.
89function global_area_mean_v(var, G, tmp_scale)
90 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
91 real, dimension(SZI_(G),SZJB_(G)), intent(in) :: var !< The variable to average in arbitrary
92 !! units [a], or arbitrary rescaled units
93 !! [A ~> a] if tmp_scale is present
94 real, optional, intent(in) :: tmp_scale !< A temporary rescaling factor for the
95 !! variable that converts it back to unscaled
96 !! (e.g., mks) units to enable the use of the
97 !! reproducing sums [a A-1 ~> 1], but is reversed
98 !! before output so that the return value has
99 !! the same units as var
100
101 real :: global_area_mean_v ! The mean of the variable in the same arbitrary units as var [A ~> a]
102
103 ! Local variables
104 ! In the following comments, [A ~> a] is used to indicate the arbitrary, possibly rescaled units of the
105 ! input array while [a] indicates the unscaled (e.g., mks) units that can be used with the reproducing sums
106 ! [A ~> a] and [B ~> b] are the same unless tmp_scale and unscale are both present.
107 real, dimension(SZI_(G),SZJ_(G)) :: tmpforsumming ! An unscaled cell integral [A L2 ~> a m2]
108 real :: temp_scale ! A temporary scaling factor [a A-1 ~> 1] or [1]
109 integer :: i, j, is, ie, js, je, isb, ieb, jsb, jeb
110
111 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
112 isb = g%iscB ; ieb = g%iecB ; jsb = g%jscB ; jeb = g%jecB
113
114 temp_scale = 1.0 ; if (present(tmp_scale)) temp_scale = tmp_scale
115
116 tmpforsumming(:,:) = 0.
117 do j=js,je ; do i=is,ie
118 tmpforsumming(i,j) = g%areaT(i,j) * &
119 (var(i,j) * g%mask2dCv(i,j) + var(i,j-1) * g%mask2dCv(i,j-1)) / &
120 max(1.e-20, g%mask2dCv(i,j)+g%mask2dCv(i,j-1))
121 enddo ; enddo
122 global_area_mean_v = reproducing_sum(tmpforsumming, unscale=g%US%L_to_m**2*temp_scale) * g%IareaT_global
123
124end function global_area_mean_v
125
126!> Return the global area mean of a variable on U grid. This uses reproducing sums.
127function global_area_mean_u(var, G, tmp_scale)
128 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
129 real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: var !< The variable to average in arbitrary
130 !! units [a], or arbitrary rescaled units
131 !! [A ~> a] if tmp_scale is present
132 real, optional, intent(in) :: tmp_scale !< A temporary rescaling factor for the
133 !! variable that converts it back to unscaled
134 !! (e.g., mks) units to enable the use of the
135 !! reproducing sums [a A-1 ~> 1], but is reversed
136 !! before output so that the return value has
137 !! the same units as var
138 real :: global_area_mean_u ! The mean of the variable in the same arbitrary units as var [A ~> a]
139
140 ! Local variables
141 ! In the following comments, [A ~> a] is used to indicate the arbitrary, possibly rescaled units of the
142 ! input array while [a] indicates the unscaled (e.g., mks) units that can be used with the reproducing sums
143 real, dimension(SZI_(G),SZJ_(G)) :: tmpforsumming ! An unscaled cell integral [A L2 ~> a m2]
144 real :: temp_scale ! A temporary scaling factor [a A-1 ~> 1] or [1]
145 integer :: i, j, is, ie, js, je, isb, ieb, jsb, jeb
146
147 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
148 isb = g%iscB ; ieb = g%iecB ; jsb = g%jscB ; jeb = g%jecB
149
150 temp_scale = 1.0 ; if (present(tmp_scale)) temp_scale = tmp_scale
151
152 tmpforsumming(:,:) = 0.
153 do j=js,je ; do i=is,ie
154 tmpforsumming(i,j) = g%areaT(i,j) * &
155 (var(i,j) * g%mask2dCu(i,j) + var(i-1,j) * g%mask2dCu(i-1,j)) / &
156 max(1.e-20, g%mask2dCu(i,j)+g%mask2dCu(i-1,j))
157 enddo ; enddo
158 global_area_mean_u = reproducing_sum(tmpforsumming, unscale=g%US%L_to_m**2*temp_scale) * g%IareaT_global
159
160end function global_area_mean_u
161
162!> Return the global area integral of a variable using reproducing sums, perhaps with a change of units.
163!! By default the integral uses the masked area from the grid, but an alternate could be used instead.
164!! The presence of the optional tmp_scale argument determines whether the returned value is in scaled
165!! (if it is present) or unscaled units for both the variable itself and for the area in the integral.
166function global_area_integral(var, G, scale, area, tmp_scale, unscale)
167 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
168 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: var !< The variable to integrate in arbitrary units [a],
169 !! or arbitrary rescaled units [A ~> a] if unscale
170 !! or tmp_scale is present
171 real, optional, intent(in) :: scale !< A rescaling factor for the variable [a A-1 ~> 1]
172 !! that converts it back to unscaled (e.g., mks)
173 !! units to enable the use of the reproducing sums
174 real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: area !< The alternate area to use, including
175 !! any required masking [L2 ~> m2].
176 real, optional, intent(in) :: tmp_scale !< A temporary rescaling factor for the variable
177 !! that is reversed in the return value [a A-1 ~> 1],
178 !! or [b B-1 ~> 1] if unscale is also present.
179 real, optional, intent(in) :: unscale !< A rescaling factor for the variable [a A-1 ~> 1]
180 !! that converts it back to unscaled (e.g., mks)
181 !! units to enable the use of the reproducing sums, or
182 !! a factor converting between rescaled units if
183 !! tmp_scale is also present [B A-1 ~> b a-1].
184 !! Here scale and unscale are synonymous, but unscale
185 !! is preferred and takes precedence if both are present.
186 real :: global_area_integral !< The returned area integral, usually in the units of var times an area,
187 !! [a m2] or [A L2 ~> a m2] or [B L2 ~> b m2], depending on which optional
188 !! arguments are provided
189
190 ! Local variables
191 ! In the following comments, [A ~> a] is used to indicate the arbitrary, possibly rescaled units of the
192 ! input array while [a] indicates the unscaled (e.g., mks) units that can be used with the reproducing sums.
193 ! [A ~> a] and [B ~> b] are the same units unless tmp_scale and unscale are both present.
194 real, dimension(SZI_(G),SZJ_(G)) :: tmpforsumming ! An unscaled cell integral in [a m2] or
195 ! a scaled cell integral in [B L2 ~> b m2] or other units as indicated below
196 real :: scalefac ! An overall scaling factor for the areas and variable, in units of [a m2 A-1 L-2 ~> 1]
197 ! or [1] or [B m2 A-1 L-2 ~> b a-1] or [B A-1 ~> b a-1] depending on which
198 ! optional arguments are present.
199 !_______________________________________________________________________________________________
200 ! Table of units of scalefac and tmpForSumming, depending on the presence of optional arguments |
201 !_______________________________________________________________________________________________|
202 ! present(tmp_scale) | present(unscale) | scalefac units | tmpForSumming units |
203 !____________________|__________________|_________________________|_____________________________!
204 ! True | True | [B A-1 ~> b a-1] | [B L2 ~> b m2] |
205 ! True | False | [1] | [A L2 ~> a m2] |
206 ! False | True | [a m2 A-1 L-2 ~> b a-1] | [a m2] |
207 ! False | False | [m2 L-2 ~> 1] | [a m2] |
208 !____________________|__________________|_________________________|_____________________________!
209 real :: temp_scale ! A temporary scaling factor [a m2 L-2 A-1 ~> 1] or [b m2 L-2 B-1 ~> 1] or [1]
210 integer :: i, j, is, ie, js, je
211 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
212
213 if (present(tmp_scale)) then
214 temp_scale = g%US%L_to_m**2 * tmp_scale ! Units of [a m2 A-1 L-2 ~> 1] or [b m2 B-1 L-2 ~> 1]
215 scalefac = 1.0
216 else
217 temp_scale = 1.0
218 scalefac = g%US%L_to_m**2
219 endif
220 if (present(unscale)) then ; scalefac = scalefac * unscale
221 elseif (present(scale)) then ; scalefac = scalefac * scale ; endif
222
223 tmpforsumming(:,:) = 0.
224 if (present(area)) then
225 do j=js,je ; do i=is,ie
226 tmpforsumming(i,j) = var(i,j) * (scalefac * area(i,j))
227 enddo ; enddo
228 else
229 do j=js,je ; do i=is,ie
230 tmpforsumming(i,j) = var(i,j) * (scalefac * g%areaT(i,j) * g%mask2dT(i,j))
231 enddo ; enddo
232 endif
233
234 global_area_integral = reproducing_sum(tmpforsumming, unscale=temp_scale)
235
236end function global_area_integral
237
238!> Return the layerwise global thickness-weighted mean of a variable. This uses reproducing sums.
239function global_layer_mean(var, h, G, GV, scale, tmp_scale, unscale)
240 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
241 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
242 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: var !< The variable to average in arbitrary units [a],
243 !! or arbitrary rescaled units [A ~> a] if unscale
244 !! or tmp_scale is present
245 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
246 real, optional, intent(in) :: scale !< A rescaling factor for the variable [a A-1 ~> 1]
247 !! that converts it back to unscaled (e.g., mks)
248 !! units to enable the use of the reproducing sums
249 real, optional, intent(in) :: tmp_scale !< A temporary rescaling factor for the
250 !! variable for use in the reproducing sums
251 !! that is reversed in the return value [a A-1 ~> 1],
252 !! or [b B-1 ~> 1] if unscale is also present.
253 real, optional, intent(in) :: unscale !< A rescaling factor for the variable [a A-1 ~> 1]
254 !! that converts it back to unscaled (e.g., mks)
255 !! units to enable the use of the reproducing sums, or
256 !! a factor converting between rescaled units if
257 !! tmp_scale is also present [B A-1 ~> b a-1].
258 !! Here scale and unscale are synonymous, but unscale
259 !! is preferred and takes precedence.
260 real, dimension(SZK_(GV)) :: global_layer_mean !< The mean of the variable in the arbitrary scaled [A ~> a]
261 !! or [B ~> b] or unscaled [a] units of var, depending on which
262 !! optional arguments are provided
263
264 ! Local variables
265 ! In the following comments, [A ~> a] is used to indicate the arbitrary, possibly rescaled units of the
266 ! input array while [a] indicates the unscaled (e.g., mks) units that can be used with the reproducing sums
267 real :: tmpforsumming(g%isc:g%iec,g%jsc:g%jec,szk_(gv)) ! An unscaled cell integral in [L2 a m ~> a m3] or
268 ! [L2 a kg m-2 ~> a kg] or a scaled cell integral in
269 ! [L2 B m ~> b m3] or [L2 B m ~> b m3] or other units
270 ! as indicated the table below.
271 real :: weight(g%isc:g%iec,g%jsc:g%jec,szk_(gv)) ! The volume or mass of each cell, depending on whether
272 ! the model is Boussinesq, used as a weight [L2 m ~> m3]
273 ! or [L2 kg m-2 ~> kg]
274 real :: scalefac ! A scaling factor for the variable [a A-1 ~> 1] or [B A-1 ~> b a-1] or [1]
275 !__________________________________________________________________________________________________
276 ! Units of weight, scalefac and tmpForSumming, depending on the presence of optional arguments |
277 !_________________________________________________________________________________________________|
278 ! Boussinesq | tmp_scale | unscale | weight units | scalefac units | tmpForSumming units |
279 ! | present | present | | | |
280 !____________|___________|_________|___________________|__________________|_______________________!
281 ! True | True | True | [L2 m ~> m3] | [B A-1 ~> b a-1] | [B L2 m ~> b m3] |
282 ! True | True | False | [L2 m ~> m3] | [1] | [A L2 m ~> a m3] |
283 ! True | False | True | [L2 m ~> m3] | [a A-1 ~> 1] | [L2 a m ~> a m3] |
284 ! True | False | False | [L2 m ~> m3] | [1] | [L2 a m ~> a m3] |
285 ! False | True | True | [L2 kg m-2 ~> kg] | [B A-1 ~> b a-1] | [B L2 kg m-2 ~> b kg] |
286 ! False | True | False | [L2 kg m-2 ~> kg] | [1] | [A L2 kg m-2 ~> a kg] |
287 ! False | False | True | [L2 kg m-2 ~> kg] | [a A-1 ~> 1] | [L2 a kg m-2 ~> a kg] |
288 ! False | False | False | [L2 kg m-2 ~> kg] | [1] | [L2 a kg m-2 ~> a kg] |
289 !____________|___________|_________|___________________|__________________|_______________________!
290 type(efp_type) :: laysums(2*szk_(gv)) ! A vector of sums with heterogeneous meanings, with the first
291 ! half being the tracer integrals in [b m3] or [b kg] and the
292 ! second half being the summed weights in [m3] or [kg]
293 real :: global_temp_scalar ! The global integral of the tracer over all
294 ! layers [L2 a m ~> a m3] or [L2 a kg m-2 ~> a kg]
295 real :: global_weight_scalar ! The global integral of the volume or mass over all
296 ! layers [L2 m ~> m3] or [L2 kg m-2 ~> kg]
297 real :: temp_scale ! A temporary scaling factor [a A-1 ~> 1] or [b B-1 ~> 1] or [1]
298 integer :: i, j, k, is, ie, js, je, nz
299 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
300
301 temp_scale = 1.0 ; if (present(tmp_scale)) temp_scale = tmp_scale
302
303 scalefac = 1.0
304 if (present(unscale)) then ; scalefac = unscale
305 elseif (present(scale)) then ; scalefac = scale ; endif
306 tmpforsumming(:,:,:) = 0. ; weight(:,:,:) = 0.
307
308 do k=1,nz ; do j=js,je ; do i=is,ie
309 weight(i,j,k) = (gv%H_to_MKS * h(i,j,k)) * (g%areaT(i,j) * g%mask2dT(i,j))
310 tmpforsumming(i,j,k) = scalefac * var(i,j,k) * weight(i,j,k)
311 enddo ; enddo ; enddo
312
313 global_temp_scalar = reproducing_sum(tmpforsumming, efp_lay_sums=laysums(1:nz), only_on_pe=.true., &
314 unscale=temp_scale*g%US%L_to_m**2)
315 global_weight_scalar = reproducing_sum(weight, efp_lay_sums=laysums(nz+1:2*nz), only_on_pe=.true., &
316 unscale=g%US%L_to_m**2)
317 call efp_sum_across_pes(laysums, 2*nz)
318
319 ! Note that temp_scale appears in the denominator here because the variables returned via the
320 ! EFP_lay_sums arguments to reproducing sums stay in unscaled mks units.
321 do k=1,nz
322 global_layer_mean(k) = efp_to_real(laysums(k)) / (temp_scale*efp_to_real(laysums(nz+k)))
323 enddo
324
325end function global_layer_mean
326
327!> Find the global thickness-weighted mean of a variable. This uses reproducing sums.
328function global_volume_mean(var, h, G, GV, scale, tmp_scale, unscale)
329 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
330 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
331 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
332 intent(in) :: var !< The variable to average in arbitrary units [a],
333 !! or arbitrary rescaled units [A ~> a] if unscale
334 !! or tmp_scale is present
335 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
336 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
337 real, optional, intent(in) :: scale !< A rescaling factor for the variable [a A-1 ~> 1]
338 !! that converts it back to unscaled (e.g., mks)
339 !! units to enable the use of the reproducing sums
340 real, optional, intent(in) :: tmp_scale !< A temporary rescaling factor for the
341 !! variable that is reversed in the return value [a A-1 ~> 1],
342 !! or [b B-1 ~> 1] if unscale is also present.
343 real, optional, intent(in) :: unscale !< A rescaling factor for the variable [a A-1 ~> 1]
344 !! that converts it back to unscaled (e.g., mks)
345 !! units to enable the use of the reproducing sums, or
346 !! a factor converting between rescaled units if
347 !! tmp_scale is also present [B A-1 ~> b a-1].
348 !! Here scale and unscale are synonymous, but unscale
349 !! is preferred and takes precedence if both are present.
350 real :: global_volume_mean !< The thickness-weighted average of var in the arbitrary scaled [A ~> a] or [B ~> b] or
351 !! unscaled [a] units of var, depending on which optional arguments are provided
352
353 ! Local variables
354 ! In the following comments, [A ~> a] is used to indicate the arbitrary, possibly rescaled units of the
355 ! input array while [a] indicates the unscaled (e.g., mks) units that can be used with the reproducing sums
356 ! [A ~> a] and [B ~> b] are the same units unless tmp_scale and unscale are both present.
357 real :: temp_scale ! A temporary scaling factor [a A-1 ~> 1] or [b B-1 ~> 1] or [1]
358 real :: scalefac ! A scaling factor for the variable [a A-1 ~> 1] or [B A-1 ~> b a-1] or [1]
359 real :: weight_here ! The volume or mass of a grid cell [L2 m ~> m3] or [L2 kg m-2 ~> kg]
360 real, dimension(SZI_(G),SZJ_(G)) :: tmpforsumming ! The volume or mass integral of the variable in a column
361 ! [B L2 m ~> b m3] or [B L2 kg m-2 ~> b kg] or
362 ! [L2 a m ~> a m3] or [L2 a kg m-2 ~> a kg]
363 real, dimension(SZI_(G),SZJ_(G)) :: sum_weight ! The volume or mass of each column of water
364 ! [L2 m ~> m3] or [L2 kg m-2 ~> kg]
365 integer :: i, j, k, is, ie, js, je, nz
366 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
367
368 temp_scale = 1.0 ; if (present(tmp_scale)) temp_scale = tmp_scale
369
370 scalefac = 1.0
371 if (present(unscale)) then ; scalefac = unscale
372 elseif (present(scale)) then ; scalefac = scale ; endif
373 tmpforsumming(:,:) = 0. ; sum_weight(:,:) = 0.
374
375 do k=1,nz ; do j=js,je ; do i=is,ie
376 weight_here = (gv%H_to_MKS * h(i,j,k)) * (g%areaT(i,j) * g%mask2dT(i,j))
377 tmpforsumming(i,j) = tmpforsumming(i,j) + scalefac * var(i,j,k) * weight_here
378 sum_weight(i,j) = sum_weight(i,j) + weight_here
379 enddo ; enddo ; enddo
380 global_volume_mean = (reproducing_sum(tmpforsumming, unscale=temp_scale*g%US%L_to_m**2)) / &
381 (reproducing_sum(sum_weight, unscale=g%US%L_to_m**2))
382
383end function global_volume_mean
384
385
386!> Find the global mass-weighted integral of a variable. The presence of the optional tmp_scale
387!! argument determines whether the returned value is in scaled (if it is present) or unscaled units
388!! for both the variable itself and for the mass in the integral. This function uses reproducing sums.
389function global_mass_integral(h, G, GV, var, on_PE_only, scale, tmp_scale, unscale)
390 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
391 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
392 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
393 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
394 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
395 optional, intent(in) :: var !< The variable to integrate in arbitrary units [a],
396 !! or arbitrary rescaled units [A ~> a] if unscale
397 !! or tmp_scale is present
398 logical, optional, intent(in) :: on_pe_only !< If present and true, the sum is only done
399 !! on the local PE, and it is _not_ order invariant.
400 real, optional, intent(in) :: scale !< A rescaling factor for the variable [a A-1 ~> 1]
401 !! that converts it back to unscaled (e.g., mks)
402 !! units to enable the use of the reproducing sums
403 real, optional, intent(in) :: tmp_scale !< A temporary rescaling factor for the variable
404 !! that is reversed in the return value [a A-1 ~> 1],
405 !! or [b B-1 ~> 1] if unscale is also present.
406 real, optional, intent(in) :: unscale !< A rescaling factor for the variable [a A-1 ~> 1]
407 !! that converts it back to unscaled (e.g., mks)
408 !! units to enable the use of the reproducing sums, or
409 !! a factor converting between rescaled units if
410 !! tmp_scale is also present [B A-1 ~> b a-1].
411 !! Here scale and unscale are synonymous, but unscale
412 !! is preferred and takes precedence if both are present.
413 real :: global_mass_integral !< The mass-weighted integral of var (or 1) in kg times the arbitrary
414 !! units of var [kg a] or in [R Z L2 A ~> kg a] if tmp_scale is present
415 !! or [R Z L2 B ~> kg b] if both unscale and tmp_scale are present
416
417 ! Local variables
418 ! In the following comments, [A ~> a] is used to indicate the arbitrary, possibly rescaled units of the
419 ! input array while [a] indicates the unscaled (e.g., mks) units that can be used with the reproducing sums
420 ! [A ~> a] and [B ~> b] are the same units unless tmp_scale and unscale are both present.
421 real :: tmpforsumming(szi_(g),szj_(g)) ! The mass-weighted integral of the variable in a column in
422 ! [kg a] or [kg] or if tmp_scale is present in [B R Z L2 ~> kg b] or
423 ! [A R Z L2 !> kg m] or [R Z L2 ~> kg]
424 real :: scalefac ! An overall scaling factor for the cell mass and variable in [a kg A-1 R-1 Z-1 L-2 ~> 1]
425 ! or [kg R-1 Z-1 L-2 ~> 1] or [1] or [B A-1 ~> b a-1] if tmp_scale is present.
426 real :: temp_scale ! A temporary scaling factor [1] or if tmp_scale is present this could be in
427 ! [kg a R-1 Z-1 L-2 A-1 ~> 1] or [kg b R-1 Z-1 L-2 B-1 ~> 1] or [kg R-1 Z-1 L-2 ~> 1]
428 !_______________________________________________________________________________________
429 ! Units of scalefac and tmpForSumming, depending on the presence of optional arguments |
430 !______________________________________________________________________________________|
431 ! var | tmp_scale | unscale | scalefac units | tmpForSumming units |
432 ! present | present | present | | |
433 !_________|___________|_________|_____________________________|________________________!
434 ! True | True | True | [B A-1 ~> b a-1] | [B R Z L2 ~> b kg] |
435 ! True | True | False | [1] | [A R Z L2 ~> a kg] |
436 ! True | False | True | [a kg A-1 R-1 Z-1 L-2 ~> 1] | [a kg] |
437 ! True | False | False | [kg R-1 Z-1 L-2 ~> 1] | [a kg] |
438 ! False | True | either | [1] | [R Z L2 ~> kg] |
439 ! False | False | either | [kg R-1 Z-1 L-2 ~> 1] | [kg] |
440 !_________|___________|_________|_____________________________|________________________!
441 logical :: global_sum ! If true do the sum globally, but if false only do the sum on the current PE.
442 integer :: i, j, k, is, ie, js, je, nz
443 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
444
445 if (present(tmp_scale)) then
446 temp_scale = g%US%RZL2_to_kg * tmp_scale
447 if (.not.present(var)) temp_scale = g%US%RZL2_to_kg
448 scalefac = 1.0
449 else
450 temp_scale = 1.0
451 scalefac = g%US%RZL2_to_kg
452 endif
453 if (present(var)) then
454 if (present(unscale)) then ; scalefac = scalefac * unscale
455 elseif (present(scale)) then ; scalefac = scalefac * scale ; endif
456 endif
457
458 tmpforsumming(:,:) = 0.0
459 if (present(var)) then
460 do k=1,nz ; do j=js,je ; do i=is,ie
461 tmpforsumming(i,j) = tmpforsumming(i,j) + var(i,j,k) * &
462 ((gv%H_to_RZ * h(i,j,k)) * (scalefac*g%areaT(i,j) * g%mask2dT(i,j)))
463 enddo ; enddo ; enddo
464 else
465 do k=1,nz ; do j=js,je ; do i=is,ie
466 tmpforsumming(i,j) = tmpforsumming(i,j) + &
467 ((gv%H_to_RZ * h(i,j,k)) * (scalefac*g%areaT(i,j) * g%mask2dT(i,j)))
468 enddo ; enddo ; enddo
469 endif
470 global_sum = .true. ; if (present(on_pe_only)) global_sum = .not.on_pe_only
471 if (global_sum) then
472 global_mass_integral = reproducing_sum(tmpforsumming, unscale=temp_scale)
473 else
475 do j=js,je ; do i=is,ie
476 global_mass_integral = global_mass_integral + tmpforsumming(i,j)
477 enddo ; enddo
478 endif
479
480end function global_mass_integral
481
482!> Find the global mass-weighted order invariant integral of a variable in mks units,
483!! returning the value as an EFP_type. This uses reproducing sums.
484function global_mass_int_efp(h, G, GV, var, on_PE_only, scale, unscale)
485 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
486 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
487 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
488 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
489 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
490 optional, intent(in) :: var !< The variable to integrate in arbitrary units [a],
491 !! or arbitrary rescaled units [A ~> a] if unscale
492 !! is present
493 logical, optional, intent(in) :: on_pe_only !< If present and true, the sum is only done
494 !! on the local PE, but it is still order invariant.
495 real, optional, intent(in) :: scale !< A rescaling factor for the variable [a A-1 ~> 1]
496 !! that converts it back to unscaled (e.g., mks)
497 !! units to enable the use of the reproducing sums
498 real, optional, intent(in) :: unscale !< A rescaling factor for the variable [a A-1 ~> 1]
499 !! that converts it back to unscaled (e.g., mks)
500 !! units to enable the use of the reproducing sums.
501 !! Here scale and unscale are synonymous, but unscale
502 !! is preferred and takes precedence if both are present.
503 type(efp_type) :: global_mass_int_efp !< The mass-weighted integral of var (or 1) in
504 !! kg times the arbitrary units of var [kg a]
505
506 ! Local variables
507 ! In the following comments, [A ~> a] is used to indicate the arbitrary, possibly rescaled units of the
508 ! input array while [a] indicates the unscaled (e.g., mks) units that can be used with the reproducing sums
509 real :: tmpforsum(szi_(g),szj_(g)) ! The mass-weighted integral of the variable in a column [kg a] or [kg]
510 real :: scalefac ! An overall scaling factor for the cell mass and variable [a kg A-1 H-1 L-2 ~> kg m-3 or 1]
511 integer :: i, j, k, is, ie, js, je, nz, isr, ier, jsr, jer
512
513 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
514 isr = is - (g%isd-1) ; ier = ie - (g%isd-1) ; jsr = js - (g%jsd-1) ; jer = je - (g%jsd-1)
515
516 scalefac = gv%H_to_kg_m2 * g%US%L_to_m**2
517 if (present(unscale)) then ; scalefac = unscale * scalefac
518 elseif (present(scale)) then ; scalefac = scale * scalefac ; endif
519
520 tmpforsum(:,:) = 0.0
521 if (present(var)) then
522 do k=1,nz ; do j=js,je ; do i=is,ie
523 tmpforsum(i,j) = tmpforsum(i,j) + var(i,j,k) * &
524 ((scalefac * h(i,j,k)) * (g%areaT(i,j) * g%mask2dT(i,j)))
525 enddo ; enddo ; enddo
526 else
527 do k=1,nz ; do j=js,je ; do i=is,ie
528 tmpforsum(i,j) = tmpforsum(i,j) + &
529 ((scalefac * h(i,j,k)) * (g%areaT(i,j) * g%mask2dT(i,j)))
530 enddo ; enddo ; enddo
531 endif
532
533 global_mass_int_efp = reproducing_sum_efp(tmpforsum, isr, ier, jsr, jer, only_on_pe=on_pe_only)
534
535end function global_mass_int_efp
536
537
538!> Determine the global mean of a field along rows of constant i, returning it
539!! in a 1-d array using the local indexing. This uses reproducing sums.
540subroutine global_i_mean(array, i_mean, G, mask, scale, tmp_scale, unscale)
541 type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure
542 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: array !< The variable to integrate in arbitrary units [a],
543 !! or arbitrary rescaled units [A ~> a] if unscale
544 !! is present
545 real, dimension(SZJ_(G)), intent(out) :: i_mean !< Global mean of array along its i-axis [a] or [A ~> a]
546 real, dimension(SZI_(G),SZJ_(G)), &
547 optional, intent(in) :: mask !< An array used for weighting the i-mean [nondim]
548 real, optional, intent(in) :: scale !< A rescaling factor for the output variable [a A-1 ~> 1]
549 !! that converts it back to unscaled (e.g., mks)
550 !! units to enable the use of the reproducing sums
551 real, optional, intent(in) :: tmp_scale !< A rescaling factor for the internal
552 !! calculations that is removed from the output [a A-1 ~> 1]
553 real, optional, intent(in) :: unscale !< A rescaling factor for the variable [a A-1 ~> 1]
554 !! that converts it back to unscaled (e.g., mks)
555 !! units to enable the use of the reproducing sums.
556 !! Here scale and unscale are synonymous, but unscale
557 !! is preferred and takes precedence if both are present.
558
559 ! Local variables
560 ! In the following comments, [A ~> a] is used to indicate the arbitrary, possibly rescaled units of the
561 ! input array while [a] indicates the unscaled (e.g., mks) units that can be used with the reproducing sums
562 type(efp_type), allocatable, dimension(:) :: asum ! The masked sum of the variable in each row [a]
563 type(efp_type), allocatable, dimension(:) :: mask_sum ! The sum of the mask values in each row [nondim]
564 real :: scalefac ! A scaling factor for the variable [a A-1 ~> 1]
565 real :: rescale ! A factor for redoing any internal rescaling before output [A a-1 ~> 1]
566 real :: mask_sum_r ! The sum of the mask values in a row [nondim]
567 integer :: is, ie, js, je, idg_off, jdg_off
568 integer :: i, j
569
570 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
571 idg_off = g%idg_offset ; jdg_off = g%jdg_offset
572
573 scalefac = 1.0
574 if (present(unscale)) then ; scalefac = unscale
575 elseif (present(scale)) then ; scalefac = scale ; endif
576
577 rescale = 1.0
578 if (present(tmp_scale)) then ; if (tmp_scale /= 0.0) then
579 scalefac = scalefac * tmp_scale ; rescale = 1.0 / tmp_scale
580 endif ; endif
581 call reset_efp_overflow_error()
582
583 allocate(asum(g%jsg:g%jeg))
584 if (present(mask)) then
585 allocate(mask_sum(g%jsg:g%jeg))
586
587 do j=g%jsg,g%jeg
588 asum(j) = real_to_efp(0.0) ; mask_sum(j) = real_to_efp(0.0)
589 enddo
590
591 do j=js,je ; do i=is,ie
592 asum(j+jdg_off) = asum(j+jdg_off) + real_to_efp(scalefac*array(i,j)*mask(i,j))
593 mask_sum(j+jdg_off) = mask_sum(j+jdg_off) + real_to_efp(mask(i,j))
594 enddo ; enddo
595
596 if (query_efp_overflow_error()) call mom_error(fatal, &
597 "global_i_mean overflow error occurred before sums across PEs.")
598
599 call efp_sum_across_pes(asum(g%jsg:g%jeg), g%jeg-g%jsg+1)
600 call efp_sum_across_pes(mask_sum(g%jsg:g%jeg), g%jeg-g%jsg+1)
601
602 if (query_efp_overflow_error()) call mom_error(fatal, &
603 "global_i_mean overflow error occurred during sums across PEs.")
604
605 do j=js,je
606 mask_sum_r = efp_to_real(mask_sum(j+jdg_off))
607 if (mask_sum_r == 0.0 ) then ; i_mean(j) = 0.0 ; else
608 i_mean(j) = efp_to_real(asum(j+jdg_off)) / mask_sum_r
609 endif
610 enddo
611
612 deallocate(mask_sum)
613 else
614 do j=g%jsg,g%jeg ; asum(j) = real_to_efp(0.0) ; enddo
615
616 do j=js,je ; do i=is,ie
617 asum(j+jdg_off) = asum(j+jdg_off) + real_to_efp(scalefac*array(i,j))
618 enddo ; enddo
619
620 if (query_efp_overflow_error()) call mom_error(fatal, &
621 "global_i_mean overflow error occurred before sum across PEs.")
622
623 call efp_sum_across_pes(asum(g%jsg:g%jeg), g%jeg-g%jsg+1)
624
625 if (query_efp_overflow_error()) call mom_error(fatal, &
626 "global_i_mean overflow error occurred during sum across PEs.")
627
628 do j=js,je
629 i_mean(j) = efp_to_real(asum(j+jdg_off)) / real(g%ieg-g%isg+1)
630 enddo
631 endif
632
633 if (rescale /= 1.0) then ; do j=js,je ; i_mean(j) = rescale*i_mean(j) ; enddo ; endif
634
635 deallocate(asum)
636
637end subroutine global_i_mean
638
639!> Determine the global mean of a field along rows of constant j, returning it
640!! in a 1-d array using the local indexing. This uses reproducing sums.
641subroutine global_j_mean(array, j_mean, G, mask, scale, tmp_scale, unscale)
642 type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure
643 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: array !< The variable to integrate in arbitrary units [a],
644 !! or arbitrary rescaled units [A ~> a] if unscale
645 !! is present
646 real, dimension(SZI_(G)), intent(out) :: j_mean !< Global mean of array along its j-axis [a] or [A ~> a]
647 real, dimension(SZI_(G),SZJ_(G)), &
648 optional, intent(in) :: mask !< An array used for weighting the j-mean [nondim]
649 real, optional, intent(in) :: scale !< A rescaling factor for the output variable [a A-1 ~> 1]
650 !! that converts it back to unscaled (e.g., mks)
651 !! units to enable the use of the reproducing sums
652 real, optional, intent(in) :: tmp_scale !< A rescaling factor for the internal
653 !! calculations that is removed from the output [a A-1 ~> 1]
654 real, optional, intent(in) :: unscale !< A rescaling factor for the variable [a A-1 ~> 1]
655 !! that converts it back to unscaled (e.g., mks)
656 !! units to enable the use of the reproducing sums.
657 !! Here scale and unscale are synonymous, but unscale
658 !! is preferred and takes precedence if both are present.
659
660 ! Local variables
661 ! In the following comments, [A ~> a] is used to indicate the arbitrary, possibly rescaled units of the
662 ! input array while [a] indicates the unscaled (e.g., mks) units that can be used with the reproducing sums
663 type(efp_type), allocatable, dimension(:) :: asum ! The masked sum of the variable in each row [a]
664 type(efp_type), allocatable, dimension(:) :: mask_sum ! The sum of the mask values in each row [nondim]
665 real :: mask_sum_r ! The sum of the mask values in a row [nondim]
666 real :: scalefac ! A scaling factor for the variable [a A-1 ~> 1]
667 real :: rescale ! A factor for redoing any internal rescaling before output [A a-1 ~> 1]
668 integer :: is, ie, js, je, idg_off, jdg_off
669 integer :: i, j
670
671 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
672 idg_off = g%idg_offset ; jdg_off = g%jdg_offset
673
674 scalefac = 1.0
675 if (present(unscale)) then ; scalefac = unscale
676 elseif (present(scale)) then ; scalefac = scale ; endif
677
678 rescale = 1.0
679 if (present(tmp_scale)) then ; if (tmp_scale /= 0.0) then
680 scalefac = scalefac * tmp_scale ; rescale = 1.0 / tmp_scale
681 endif ; endif
682 call reset_efp_overflow_error()
683
684 allocate(asum(g%isg:g%ieg))
685 if (present(mask)) then
686 allocate (mask_sum(g%isg:g%ieg))
687
688 do i=g%isg,g%ieg
689 asum(i) = real_to_efp(0.0) ; mask_sum(i) = real_to_efp(0.0)
690 enddo
691
692 do i=is,ie ; do j=js,je
693 asum(i+idg_off) = asum(i+idg_off) + real_to_efp(scalefac*array(i,j)*mask(i,j))
694 mask_sum(i+idg_off) = mask_sum(i+idg_off) + real_to_efp(mask(i,j))
695 enddo ; enddo
696
697 if (query_efp_overflow_error()) call mom_error(fatal, &
698 "global_j_mean overflow error occurred before sums across PEs.")
699
700 call efp_sum_across_pes(asum(g%isg:g%ieg), g%ieg-g%isg+1)
701 call efp_sum_across_pes(mask_sum(g%isg:g%ieg), g%ieg-g%isg+1)
702
703 if (query_efp_overflow_error()) call mom_error(fatal, &
704 "global_j_mean overflow error occurred during sums across PEs.")
705
706 do i=is,ie
707 mask_sum_r = efp_to_real(mask_sum(i+idg_off))
708 if (mask_sum_r == 0.0 ) then ; j_mean(i) = 0.0 ; else
709 j_mean(i) = efp_to_real(asum(i+idg_off)) / mask_sum_r
710 endif
711 enddo
712
713 deallocate(mask_sum)
714 else
715 do i=g%isg,g%ieg ; asum(i) = real_to_efp(0.0) ; enddo
716
717 do i=is,ie ; do j=js,je
718 asum(i+idg_off) = asum(i+idg_off) + real_to_efp(scalefac*array(i,j))
719 enddo ; enddo
720
721 if (query_efp_overflow_error()) call mom_error(fatal, &
722 "global_j_mean overflow error occurred before sum across PEs.")
723
724 call efp_sum_across_pes(asum(g%isg:g%ieg), g%ieg-g%isg+1)
725
726 if (query_efp_overflow_error()) call mom_error(fatal, &
727 "global_j_mean overflow error occurred during sum across PEs.")
728
729 do i=is,ie
730 j_mean(i) = efp_to_real(asum(i+idg_off)) / real(g%jeg-g%jsg+1)
731 enddo
732 endif
733
734 if (rescale /= 1.0) then ; do i=is,ie ; j_mean(i) = rescale*j_mean(i) ; enddo ; endif
735
736 deallocate(asum)
737
738end subroutine global_j_mean
739
740!> Adjust 2d array such that area mean is zero without moving the zero contour
741subroutine adjust_area_mean_to_zero(array, G, scaling, unit_scale, unscale)
742 type(ocean_grid_type), intent(in) :: g !< Grid structure
743 real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: array !< 2D array to be adjusted in arbitrary units [a],
744 !! or arbitrary rescaled units [A ~> a] if unscale
745 !! is present
746 real, optional, intent(out) :: scaling !< The scaling factor used [nondim]
747 real, optional, intent(in) :: unit_scale !< A rescaling factor for the variable [a A-1 ~> 1]
748 !! that converts it back to unscaled (e.g., mks)
749 !! units to enable the use of the reproducing sums
750 real, optional, intent(in) :: unscale !< A rescaling factor for the variable [a A-1 ~> 1]
751 !! that converts it back to unscaled (e.g., mks)
752 !! units to enable the use of the reproducing sums.
753 !! Here unit_scale and unscale are synonymous, but unscale
754 !! is preferred and takes precedence if both are present.
755 ! Local variables
756 ! In the following comments, [A ~> a] is used to indicate the arbitrary, possibly rescaled units of the
757 ! input array while [a] indicates the unscaled (e.g., mks) units that can be used with the reproducing sums
758 real :: posvals(g%isc:g%iec, g%jsc:g%jec) ! The positive values in a cell or 0 [A ~> a]
759 real :: negvals(g%isc:g%iec, g%jsc:g%jec) ! The negative values in a cell or 0 [A ~> a]
760 real :: areaxposvals(g%isc:g%iec, g%jsc:g%jec) ! The cell area integral of the positive values [L2 A ~> m2 a]
761 real :: areaxnegvals(g%isc:g%iec, g%jsc:g%jec) ! The cell area integral of the negative values [L2 A ~> m2 a]
762 type(efp_type), dimension(2) :: areaint_efp ! An EFP version integral of the values on the current PE [m2 a]
763 real :: scalefac ! A scaling factor for the variable [a A-1 ~> 1]
764 real :: areaintposvals, areaintnegvals ! The global area integral of the positive and negative values [m2 a]
765 real :: posscale, negscale ! The scaling factor to apply to positive or negative values [nondim]
766 integer :: i,j
767
768 scalefac = 1.0
769 if (present(unscale)) then ; scalefac = unscale
770 elseif (present(unit_scale)) then ; scalefac = unit_scale ; endif
771
772 ! areaXposVals(:,:) = 0. ! This zeros out halo points.
773 ! areaXnegVals(:,:) = 0. ! This zeros out halo points.
774
775 do j=g%jsc,g%jec ; do i=g%isc,g%iec
776 posvals(i,j) = max(0., array(i,j))
777 areaxposvals(i,j) = g%areaT(i,j) * posvals(i,j)
778 negvals(i,j) = min(0., array(i,j))
779 areaxnegvals(i,j) = g%areaT(i,j) * negvals(i,j)
780 enddo ; enddo
781
782 ! Combining the sums like this avoids separate blocking global sums.
783 areaint_efp(1) = reproducing_sum_efp( areaxposvals, only_on_pe=.true., unscale=scalefac*g%US%L_to_m**2 )
784 areaint_efp(2) = reproducing_sum_efp( areaxnegvals, only_on_pe=.true., unscale=scalefac*g%US%L_to_m**2 )
785 call efp_sum_across_pes(areaint_efp, 2)
786 areaintposvals = efp_to_real( areaint_efp(1) )
787 areaintnegvals = efp_to_real( areaint_efp(2) )
788
789 posscale = 0.0 ; negscale = 0.0
790 if ((areaintposvals>0.).and.(areaintnegvals<0.)) then ! Only adjust if possible
791 if (areaintposvals>-areaintnegvals) then ! Scale down positive values
792 posscale = - areaintnegvals / areaintposvals
793 do j=g%jsc,g%jec ; do i=g%isc,g%iec
794 array(i,j) = (posscale * posvals(i,j)) + negvals(i,j)
795 enddo ; enddo
796 elseif (areaintposvals<-areaintnegvals) then ! Scale down negative values
797 negscale = - areaintposvals / areaintnegvals
798 do j=g%jsc,g%jec ; do i=g%isc,g%iec
799 array(i,j) = posvals(i,j) + (negscale * negvals(i,j))
800 enddo ; enddo
801 endif
802 endif
803 if (present(scaling)) scaling = posscale - negscale
804
805end subroutine adjust_area_mean_to_zero
806
807
808!> Find the global maximum and minimum of a tracer array and return the locations of the extrema.
809!! When there multiple cells with the same extreme values, the reported locations are from the
810!! uppermost layer where they occur, and then from the logically northernmost and then eastermost
811!! such location on the unrotated version of the grid within that layer. Only ocean points (as
812!! indicated by a positive value of G%mask2dT) are evaluated, and if there are no ocean points
813!! anywhere in the domain, the reported extrema and their locations are all returned as 0.
814subroutine array_global_min_max(tr_array, G, nk, g_min, g_max, &
815 xgmin, ygmin, zgmin, xgmax, ygmax, zgmax, unscale)
816 integer, intent(in) :: nk !< The number of vertical levels
817 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
818 real, dimension(SZI_(G),SZJ_(G),nk), intent(in) :: tr_array !< The tracer array to search for
819 !! extrema in arbitrary concentration units [CU ~> conc]
820 real, intent(out) :: g_min !< The global minimum of tr_array, either in
821 !! the same units as tr_array [CU ~> conc] or in
822 !! unscaled units if unscale is present [conc]
823 real, intent(out) :: g_max !< The global maximum of tr_array, either in
824 !! the same units as tr_array [CU ~> conc] or in
825 !! unscaled units if unscale is present [conc]
826 real, optional, intent(out) :: xgmin !< The x-position of the global minimum in the
827 !! units of G%geoLonT, often [degrees_E] or [km] or [m]
828 real, optional, intent(out) :: ygmin !< The y-position of the global minimum in the
829 !! units of G%geoLatT, often [degrees_N] or [km] or [m]
830 real, optional, intent(out) :: zgmin !< The z-position of the global minimum [layer]
831 real, optional, intent(out) :: xgmax !< The x-position of the global maximum in the
832 !! units of G%geoLonT, often [degrees_E] or [km] or [m]
833 real, optional, intent(out) :: ygmax !< The y-position of the global maximum in the
834 !! units of G%geoLatT, often [degrees_N] or [km] or [m]
835 real, optional, intent(out) :: zgmax !< The z-position of the global maximum [layer]
836 real, optional, intent(in) :: unscale !< A factor to use to undo any scaling of
837 !! the input tracer array [conc CU-1 ~> 1]
838
839 ! Local variables
840 real :: tmax, tmin ! Maximum and minimum tracer values, in the same units as tr_array [CU ~> conc]
841 integer :: ijk_min_max(2) ! Integers encoding the global grid positions of the global minimum and maximum values
842 real :: xyz_min_max(6) ! A single array with the x-, y- and z-positions of the minimum and
843 ! maximum values in units that vary between the array elements [various]
844 logical :: valid_pe ! True if there are any valid points on the local PE.
845 logical :: find_location ! If true, report the locations of the extrema
846 integer :: ijk_loc_max ! An integer encoding the global grid position of the maximum tracer value on this PE
847 integer :: ijk_loc_min ! An integer encoding the global grid position of the minimum tracer value on this PE
848 integer :: ijk_loc_here ! An integer encoding the global grid position of the current grid point
849 integer :: itmax, jtmax, ktmax, itmin, jtmin, ktmin
850 integer :: i, j, k, isc, iec, jsc, jec
851
852 isc = g%isc ; iec = g%iec ; jsc = g%jsc ; jec = g%jec
853
854 find_location = (present(xgmin) .or. present(ygmin) .or. present(zgmin) .or. &
855 present(xgmax) .or. present(ygmax) .or. present(zgmax))
856
857 ! The initial values set here are never used if there are any valid points.
858 tmax = -huge(tmax) ; tmin = huge(tmin)
859
860 if (find_location) then
861 ! Find the maximum and minimum tracer values on this PE and their locations.
862 valid_pe = .false.
863 itmax = 0 ; jtmax = 0 ; ktmax = 0 ; ijk_loc_max = 0
864 itmin = 0 ; jtmin = 0 ; ktmin = 0 ; ijk_loc_min = 0
865 do k=1,nk ; do j=jsc,jec ; do i=isc,iec ; if (g%mask2dT(i,j) > 0.0) then
866 valid_pe = .true.
867 if (tr_array(i,j,k) > tmax) then
868 tmax = tr_array(i,j,k)
869 itmax = i ; jtmax = j ; ktmax = k
870 ijk_loc_max = ijk_loc(i, j, k, nk, g%HI)
871 elseif ((tr_array(i,j,k) == tmax) .and. (k <= ktmax)) then
872 ijk_loc_here = ijk_loc(i, j, k, nk, g%HI)
873 if (ijk_loc_here > ijk_loc_max) then
874 itmax = i ; jtmax = j ; ktmax = k
875 ijk_loc_max = ijk_loc_here
876 endif
877 endif
878 if (tr_array(i,j,k) < tmin) then
879 tmin = tr_array(i,j,k)
880 itmin = i ; jtmin = j ; ktmin = k
881 ijk_loc_min = ijk_loc(i, j, k, nk, g%HI)
882 elseif ((tr_array(i,j,k) == tmin) .and. (k <= ktmin)) then
883 ijk_loc_here = ijk_loc(i, j, k, nk, g%HI)
884 if (ijk_loc_here > ijk_loc_min) then
885 itmin = i ; jtmin = j ; ktmin = k
886 ijk_loc_min = ijk_loc_here
887 endif
888 endif
889 endif ; enddo ; enddo ; enddo
890 else
891 ! Only the maximum and minimum values are needed, and not their positions.
892 do k=1,nk ; do j=jsc,jec ; do i=isc,iec ; if (g%mask2dT(i,j) > 0.0) then
893 if (tr_array(i,j,k) > tmax) tmax = tr_array(i,j,k)
894 if (tr_array(i,j,k) < tmin) tmin = tr_array(i,j,k)
895 endif ; enddo ; enddo ; enddo
896 endif
897
898 ! Find the global maximum and minimum tracer values.
899 g_max = tmax ; g_min = tmin
900 call max_across_pes(g_max)
901 call min_across_pes(g_min)
902
903 if (find_location) then
904 if (g_max < g_min) then
905 ! This only occurs if there are no unmasked points anywhere in the domain.
906 xyz_min_max(:) = 0.0
907 else
908 ! Find the global indices of the maximum and minimum locations. This can
909 ! occur on multiple PEs.
910 ijk_min_max(1:2) = 0
911 if (valid_pe) then
912 if (g_min == tmin) ijk_min_max(1) = ijk_loc_min
913 if (g_max == tmax) ijk_min_max(2) = ijk_loc_max
914 endif
915 ! If MOM6 supported taking maxima on arrays of integers, these could be combined as:
916 ! call max_across_PEs(ijk_min_max, 2)
917 call max_across_pes(ijk_min_max(1))
918 call max_across_pes(ijk_min_max(2))
919
920 ! Set the positions of the extrema if they occur on this PE. This will only
921 ! occur on a single PE.
922 xyz_min_max(1:6) = -huge(xyz_min_max) ! These huge negative values are never selected by max_across_PEs.
923 if (valid_pe) then
924 if (ijk_min_max(1) == ijk_loc_min) then
925 xyz_min_max(1) = g%geoLonT(itmin,jtmin)
926 xyz_min_max(2) = g%geoLatT(itmin,jtmin)
927 xyz_min_max(3) = real(ktmin)
928 endif
929 if (ijk_min_max(2) == ijk_loc_max) then
930 xyz_min_max(4) = g%geoLonT(itmax,jtmax)
931 xyz_min_max(5) = g%geoLatT(itmax,jtmax)
932 xyz_min_max(6) = real(ktmax)
933 endif
934 endif
935
936 call max_across_pes(xyz_min_max, 6)
937 endif
938
939 if (present(xgmin)) xgmin = xyz_min_max(1)
940 if (present(ygmin)) ygmin = xyz_min_max(2)
941 if (present(zgmin)) zgmin = xyz_min_max(3)
942 if (present(xgmax)) xgmax = xyz_min_max(4)
943 if (present(ygmax)) ygmax = xyz_min_max(5)
944 if (present(zgmax)) zgmax = xyz_min_max(6)
945 endif
946
947 if (g_max < g_min) then
948 ! There are no unmasked points anywhere in the domain.
949 g_max = 0.0 ; g_min = 0.0
950 endif
951
952 if (present(unscale)) then
953 ! Rescale g_min and g_max, perhaps changing their units from [CU ~> conc] to [conc]
954 g_max = unscale * g_max
955 g_min = unscale * g_min
956 endif
957
958end subroutine array_global_min_max
959
960! Return a positive integer encoding the rotationally invariant global position of a tracer cell
961function ijk_loc(i, j, k, nk, HI)
962 integer, intent(in) :: i !< Local i-index
963 integer, intent(in) :: j !< Local j-index
964 integer, intent(in) :: k !< Local k-index
965 integer, intent(in) :: nk !< Range of k-index, used to pick out a low-k position.
966 type(hor_index_type), intent(in) :: hi !< Horizontal index ranges
967 integer :: ijk_loc ! An integer encoding the cell position in the global grid.
968
969 ! Local variables
970 integer :: ig, jg ! Global index values with a global computational domain start value of 1.
971 integer :: ij_loc ! The encoding of the horizontal position
972 integer :: qturns ! The number of counter-clockwise quarter turns of the grid that have to be undone
973
974 ! These global i-grid positions run from 1 to HI%niglobal, and analogously for jg.
975 ig = i + hi%idg_offset + (1 - hi%isg)
976 jg = j + hi%jdg_offset + (1 - hi%jsg)
977
978 ! Compensate for the rotation of the model grid to give a rotationally invariant encoding.
979 qturns = modulo(hi%turns, 4)
980 if (qturns == 0) then
981 ij_loc = ig + hi%niglobal * jg
982 elseif (qturns == 1) then
983 ij_loc = jg + hi%njglobal * ((hi%niglobal+1)-ig)
984 elseif (qturns == 2) then
985 ij_loc = ((hi%niglobal+1)-ig) + hi%niglobal * ((hi%njglobal+1)-jg)
986 elseif (qturns == 3) then
987 ij_loc = ((hi%njglobal+1)-jg) + hi%njglobal * ig
988 endif
989
990 ijk_loc = ij_loc + (hi%niglobal*hi%njglobal) * (nk-k)
991
992end function ijk_loc
993
994
995end module mom_spatial_means