portedportable, not yet portedexecuted, not portableexecutable, not hit by this run
| 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 | |
| 6 | module MOM_spatial_means | |
| 7 | ||
| 8 | use MOM_coms, only : EFP_type, operator(+), operator(-), assignment(=) | |
| 9 | use MOM_coms, only : EFP_to_real, real_to_EFP, EFP_sum_across_PEs | |
| 10 | use MOM_coms, only : reproducing_sum, reproducing_sum_EFP, EFP_to_real | |
| 11 | use MOM_coms, only : query_EFP_overflow_error, reset_EFP_overflow_error | |
| 12 | use MOM_coms, only : max_across_PEs, min_across_PEs | |
| 13 | use MOM_error_handler, only : MOM_error, NOTE, WARNING, FATAL, is_root_pe | |
| 14 | use MOM_file_parser, only : get_param, log_version, param_file_type | |
| 15 | use MOM_grid, only : ocean_grid_type | |
| 16 | use MOM_hor_index, only : hor_index_type | |
| 17 | use MOM_verticalGrid, only : verticalGrid_type | |
| 18 | ||
| 19 | implicit none ; private | |
| 20 | ||
| 21 | #include <MOM_memory.h> | |
| 22 | ||
| 23 | public :: global_i_mean, global_j_mean | |
| 24 | public :: global_area_mean, global_area_mean_u, global_area_mean_v, global_layer_mean | |
| 25 | public :: global_area_integral | |
| 26 | public :: global_volume_mean, global_mass_integral, global_mass_int_EFP | |
| 27 | public :: adjust_area_mean_to_zero | |
| 28 | public :: array_global_min_max | |
| 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 | ||
| 37 | contains | |
| 38 | ||
| 39 | !> Return the global area mean of a variable, perhaps with a change of units. This uses reproducing sums. | |
| 40 | 12 | function 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 | 24 | 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 | 12 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec |
| 72 | ||
| 73 | 12 | temp_scale = 1.0 ; if (present(tmp_scale)) temp_scale = tmp_scale |
| 74 | ||
| 75 | 12 | scalefac = 1.0 |
| 76 | 12 | if (present(unscale)) then ; scalefac = unscale |
| 77 | 12 | elseif (present(scale)) then ; scalefac = scale ; endif |
| 78 | ||
| 79 | 105276 | tmpForSumming(:,:) = 0. |
| 80 | 87132 | do j=js,je ; do i=is,ie |
| 81 | 87120 | tmpForSumming(i,j) = var(i,j) * (scalefac * G%areaT(i,j) * G%mask2dT(i,j)) |
| 82 | enddo ; enddo | |
| 83 | ||
| 84 | 12 | global_area_mean = reproducing_sum(tmpForSumming, unscale=temp_scale*G%US%L_to_m**2) * G%IareaT_global |
| 85 | ||
| 86 | 12 | end function global_area_mean |
| 87 | ||
| 88 | !> Return the global area mean of a variable. This uses reproducing sums. | |
| 89 | 0 | function 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 | 0 | 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 | 0 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec |
| 112 | 0 | isB = G%iscB ; ieB = G%iecB ; jsB = G%jscB ; jeB = G%jecB |
| 113 | ||
| 114 | 0 | temp_scale = 1.0 ; if (present(tmp_scale)) temp_scale = tmp_scale |
| 115 | ||
| 116 | 0 | tmpForSumming(:,:) = 0. |
| 117 | 0 | 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 | 0 | max(1.e-20, G%mask2dCv(i,J)+G%mask2dCv(i,J-1)) |
| 121 | enddo ; enddo | |
| 122 | 0 | global_area_mean_v = reproducing_sum(tmpForSumming, unscale=G%US%L_to_m**2*temp_scale) * G%IareaT_global |
| 123 | ||
| 124 | 0 | end function global_area_mean_v |
| 125 | ||
| 126 | !> Return the global area mean of a variable on U grid. This uses reproducing sums. | |
| 127 | 0 | function 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 | 0 | 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 | 0 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec |
| 148 | 0 | isB = G%iscB ; ieB = G%iecB ; jsB = G%jscB ; jeB = G%jecB |
| 149 | ||
| 150 | 0 | temp_scale = 1.0 ; if (present(tmp_scale)) temp_scale = tmp_scale |
| 151 | ||
| 152 | 0 | tmpForSumming(:,:) = 0. |
| 153 | 0 | 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 | 0 | max(1.e-20, G%mask2dCu(I,j)+G%mask2dCu(I-1,j)) |
| 157 | enddo ; enddo | |
| 158 | 0 | global_area_mean_u = reproducing_sum(tmpForSumming, unscale=G%US%L_to_m**2*temp_scale) * G%IareaT_global |
| 159 | ||
| 160 | 0 | end 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. | |
| 166 | 0 | function 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 | 0 | 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 | 0 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec |
| 212 | ||
| 213 | 0 | if (present(tmp_scale)) then |
| 214 | 0 | 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 | 0 | scalefac = 1.0 |
| 216 | else | |
| 217 | 0 | temp_scale = 1.0 |
| 218 | 0 | scalefac = G%US%L_to_m**2 |
| 219 | endif | |
| 220 | 0 | if (present(unscale)) then ; scalefac = scalefac * unscale |
| 221 | 0 | elseif (present(scale)) then ; scalefac = scalefac * scale ; endif |
| 222 | ||
| 223 | 0 | tmpForSumming(:,:) = 0. |
| 224 | 0 | if (present(area)) then |
| 225 | 0 | do j=js,je ; do i=is,ie |
| 226 | 0 | tmpForSumming(i,j) = var(i,j) * (scalefac * area(i,j)) |
| 227 | enddo ; enddo | |
| 228 | else | |
| 229 | 0 | do j=js,je ; do i=is,ie |
| 230 | 0 | tmpForSumming(i,j) = var(i,j) * (scalefac * G%areaT(i,j) * G%mask2dT(i,j)) |
| 231 | enddo ; enddo | |
| 232 | endif | |
| 233 | ||
| 234 | 0 | global_area_integral = reproducing_sum(tmpForSumming, unscale=temp_scale) |
| 235 | ||
| 236 | 0 | end function global_area_integral |
| 237 | ||
| 238 | !> Return the layerwise global thickness-weighted mean of a variable. This uses reproducing sums. | |
| 239 | 0 | function 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 | 0 | 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 | 0 | 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 | 0 | 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 | 0 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke |
| 300 | ||
| 301 | 0 | temp_scale = 1.0 ; if (present(tmp_scale)) temp_scale = tmp_scale |
| 302 | ||
| 303 | 0 | scalefac = 1.0 |
| 304 | 0 | if (present(unscale)) then ; scalefac = unscale |
| 305 | 0 | elseif (present(scale)) then ; scalefac = scale ; endif |
| 306 | 0 | tmpForSumming(:,:,:) = 0. ; weight(:,:,:) = 0. |
| 307 | ||
| 308 | 0 | do k=1,nz ; do j=js,je ; do i=is,ie |
| 309 | 0 | weight(i,j,k) = (GV%H_to_MKS * h(i,j,k)) * (G%areaT(i,j) * G%mask2dT(i,j)) |
| 310 | 0 | 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 | 0 | 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 | 0 | unscale=G%US%L_to_m**2) |
| 317 | 0 | 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 | 0 | do k=1,nz |
| 322 | 0 | global_layer_mean(k) = EFP_to_real(laysums(k)) / (temp_scale*EFP_to_real(laysums(nz+k))) |
| 323 | enddo | |
| 324 | ||
| 325 | end function global_layer_mean | |
| 326 | ||
| 327 | !> Find the global thickness-weighted mean of a variable. This uses reproducing sums. | |
| 328 | 0 | function 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 | 0 | 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 | 0 | 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 | 0 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke |
| 367 | ||
| 368 | 0 | temp_scale = 1.0 ; if (present(tmp_scale)) temp_scale = tmp_scale |
| 369 | ||
| 370 | 0 | scalefac = 1.0 |
| 371 | 0 | if (present(unscale)) then ; scalefac = unscale |
| 372 | 0 | elseif (present(scale)) then ; scalefac = scale ; endif |
| 373 | 0 | tmpForSumming(:,:) = 0. ; sum_weight(:,:) = 0. |
| 374 | ||
| 375 | 0 | do k=1,nz ; do j=js,je ; do i=is,ie |
| 376 | 0 | weight_here = (GV%H_to_MKS * h(i,j,k)) * (G%areaT(i,j) * G%mask2dT(i,j)) |
| 377 | 0 | tmpForSumming(i,j) = tmpForSumming(i,j) + scalefac * var(i,j,k) * weight_here |
| 378 | 0 | 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 | 0 | (reproducing_sum(sum_weight, unscale=G%US%L_to_m**2)) |
| 382 | ||
| 383 | 0 | end 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. | |
| 389 | 0 | function 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 | 0 | 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 | 0 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke |
| 444 | ||
| 445 | 0 | if (present(tmp_scale)) then |
| 446 | 0 | temp_scale = G%US%RZL2_to_kg * tmp_scale |
| 447 | 0 | if (.not.present(var)) temp_scale = G%US%RZL2_to_kg |
| 448 | 0 | scalefac = 1.0 |
| 449 | else | |
| 450 | 0 | temp_scale = 1.0 |
| 451 | 0 | scalefac = G%US%RZL2_to_kg |
| 452 | endif | |
| 453 | 0 | if (present(var)) then |
| 454 | 0 | if (present(unscale)) then ; scalefac = scalefac * unscale |
| 455 | 0 | elseif (present(scale)) then ; scalefac = scalefac * scale ; endif |
| 456 | endif | |
| 457 | ||
| 458 | 0 | tmpForSumming(:,:) = 0.0 |
| 459 | 0 | if (present(var)) then |
| 460 | 0 | do k=1,nz ; do j=js,je ; do i=is,ie |
| 461 | tmpForSumming(i,j) = tmpForSumming(i,j) + var(i,j,k) * & | |
| 462 | 0 | ((GV%H_to_RZ * h(i,j,k)) * (scalefac*G%areaT(i,j) * G%mask2dT(i,j))) |
| 463 | enddo ; enddo ; enddo | |
| 464 | else | |
| 465 | 0 | do k=1,nz ; do j=js,je ; do i=is,ie |
| 466 | tmpForSumming(i,j) = tmpForSumming(i,j) + & | |
| 467 | 0 | ((GV%H_to_RZ * h(i,j,k)) * (scalefac*G%areaT(i,j) * G%mask2dT(i,j))) |
| 468 | enddo ; enddo ; enddo | |
| 469 | endif | |
| 470 | 0 | global_sum = .true. ; if (present(on_PE_only)) global_sum = .not.on_PE_only |
| 471 | 0 | if (global_sum) then |
| 472 | 0 | global_mass_integral = reproducing_sum(tmpForSumming, unscale=temp_scale) |
| 473 | else | |
| 474 | 0 | global_mass_integral = 0.0 |
| 475 | 0 | do j=js,je ; do i=is,ie |
| 476 | 0 | global_mass_integral = global_mass_integral + tmpForSumming(i,j) |
| 477 | enddo ; enddo | |
| 478 | endif | |
| 479 | ||
| 480 | 0 | end 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. | |
| 484 | 3 | function 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 | 6 | 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 | 3 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke |
| 514 | 3 | isr = is - (G%isd-1) ; ier = ie - (G%isd-1) ; jsr = js - (G%jsd-1) ; jer = je - (G%jsd-1) |
| 515 | ||
| 516 | 3 | scalefac = GV%H_to_kg_m2 * G%US%L_to_m**2 |
| 517 | 3 | if (present(unscale)) then ; scalefac = unscale * scalefac |
| 518 | 3 | elseif (present(scale)) then ; scalefac = scale * scalefac ; endif |
| 519 | ||
| 520 | 26319 | tmpForSum(:,:) = 0.0 |
| 521 | 3 | if (present(var)) then |
| 522 | 1633728 | do k=1,nz ; do j=js,je ; do i=is,ie |
| 523 | tmpForSum(i,j) = tmpForSum(i,j) + var(i,j,k) * & | |
| 524 | 1633500 | ((scalefac * h(i,j,k)) * (G%areaT(i,j) * G%mask2dT(i,j))) |
| 525 | enddo ; enddo ; enddo | |
| 526 | else | |
| 527 | 0 | do k=1,nz ; do j=js,je ; do i=is,ie |
| 528 | tmpForSum(i,j) = tmpForSum(i,j) + & | |
| 529 | 0 | ((scalefac * h(i,j,k)) * (G%areaT(i,j) * G%mask2dT(i,j))) |
| 530 | enddo ; enddo ; enddo | |
| 531 | endif | |
| 532 | ||
| 533 | 3 | global_mass_int_EFP = reproducing_sum_EFP(tmpForSum, isr, ier, jsr, jer, only_on_PE=on_PE_only) |
| 534 | ||
| 535 | 6 | end 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. | |
| 540 | 0 | subroutine 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 | 0 | type(EFP_type), allocatable, dimension(:) :: asum ! The masked sum of the variable in each row [a] |
| 563 | 0 | 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 | 0 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec |
| 571 | 0 | idg_off = G%idg_offset ; jdg_off = G%jdg_offset |
| 572 | ||
| 573 | 0 | scalefac = 1.0 |
| 574 | 0 | if (present(unscale)) then ; scalefac = unscale |
| 575 | 0 | elseif (present(scale)) then ; scalefac = scale ; endif |
| 576 | ||
| 577 | 0 | rescale = 1.0 |
| 578 | 0 | if (present(tmp_scale)) then ; if (tmp_scale /= 0.0) then |
| 579 | 0 | scalefac = scalefac * tmp_scale ; rescale = 1.0 / tmp_scale |
| 580 | endif ; endif | |
| 581 | 0 | call reset_EFP_overflow_error() |
| 582 | ||
| 583 | 0 | allocate(asum(G%jsg:G%jeg)) |
| 584 | 0 | if (present(mask)) then |
| 585 | 0 | allocate(mask_sum(G%jsg:G%jeg)) |
| 586 | ||
| 587 | 0 | do j=G%jsg,G%jeg |
| 588 | 0 | asum(j) = real_to_EFP(0.0) ; mask_sum(j) = real_to_EFP(0.0) |
| 589 | enddo | |
| 590 | ||
| 591 | 0 | do j=js,je ; do i=is,ie |
| 592 | 0 | asum(j+jdg_off) = asum(j+jdg_off) + real_to_EFP(scalefac*array(i,j)*mask(i,j)) |
| 593 | 0 | mask_sum(j+jdg_off) = mask_sum(j+jdg_off) + real_to_EFP(mask(i,j)) |
| 594 | enddo ; enddo | |
| 595 | ||
| 596 | 0 | if (query_EFP_overflow_error()) call MOM_error(FATAL, & |
| 597 | 0 | "global_i_mean overflow error occurred before sums across PEs.") |
| 598 | ||
| 599 | 0 | call EFP_sum_across_PEs(asum(G%jsg:G%jeg), G%jeg-G%jsg+1) |
| 600 | 0 | call EFP_sum_across_PEs(mask_sum(G%jsg:G%jeg), G%jeg-G%jsg+1) |
| 601 | ||
| 602 | 0 | if (query_EFP_overflow_error()) call MOM_error(FATAL, & |
| 603 | 0 | "global_i_mean overflow error occurred during sums across PEs.") |
| 604 | ||
| 605 | 0 | do j=js,je |
| 606 | 0 | mask_sum_r = EFP_to_real(mask_sum(j+jdg_off)) |
| 607 | 0 | if (mask_sum_r == 0.0 ) then ; i_mean(j) = 0.0 ; else |
| 608 | 0 | i_mean(j) = EFP_to_real(asum(j+jdg_off)) / mask_sum_r |
| 609 | endif | |
| 610 | enddo | |
| 611 | ||
| 612 | 0 | deallocate(mask_sum) |
| 613 | else | |
| 614 | 0 | do j=G%jsg,G%jeg ; asum(j) = real_to_EFP(0.0) ; enddo |
| 615 | ||
| 616 | 0 | do j=js,je ; do i=is,ie |
| 617 | 0 | asum(j+jdg_off) = asum(j+jdg_off) + real_to_EFP(scalefac*array(i,j)) |
| 618 | enddo ; enddo | |
| 619 | ||
| 620 | 0 | if (query_EFP_overflow_error()) call MOM_error(FATAL, & |
| 621 | 0 | "global_i_mean overflow error occurred before sum across PEs.") |
| 622 | ||
| 623 | 0 | call EFP_sum_across_PEs(asum(G%jsg:G%jeg), G%jeg-G%jsg+1) |
| 624 | ||
| 625 | 0 | if (query_EFP_overflow_error()) call MOM_error(FATAL, & |
| 626 | 0 | "global_i_mean overflow error occurred during sum across PEs.") |
| 627 | ||
| 628 | 0 | do j=js,je |
| 629 | 0 | i_mean(j) = EFP_to_real(asum(j+jdg_off)) / real(G%ieg-G%isg+1) |
| 630 | enddo | |
| 631 | endif | |
| 632 | ||
| 633 | 0 | if (rescale /= 1.0) then ; do j=js,je ; i_mean(j) = rescale*i_mean(j) ; enddo ; endif |
| 634 | ||
| 635 | 0 | deallocate(asum) |
| 636 | ||
| 637 | 0 | end 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. | |
| 641 | 0 | subroutine 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 | 0 | type(EFP_type), allocatable, dimension(:) :: asum ! The masked sum of the variable in each row [a] |
| 664 | 0 | 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 | 0 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec |
| 672 | 0 | idg_off = G%idg_offset ; jdg_off = G%jdg_offset |
| 673 | ||
| 674 | 0 | scalefac = 1.0 |
| 675 | 0 | if (present(unscale)) then ; scalefac = unscale |
| 676 | 0 | elseif (present(scale)) then ; scalefac = scale ; endif |
| 677 | ||
| 678 | 0 | rescale = 1.0 |
| 679 | 0 | if (present(tmp_scale)) then ; if (tmp_scale /= 0.0) then |
| 680 | 0 | scalefac = scalefac * tmp_scale ; rescale = 1.0 / tmp_scale |
| 681 | endif ; endif | |
| 682 | 0 | call reset_EFP_overflow_error() |
| 683 | ||
| 684 | 0 | allocate(asum(G%isg:G%ieg)) |
| 685 | 0 | if (present(mask)) then |
| 686 | 0 | allocate (mask_sum(G%isg:G%ieg)) |
| 687 | ||
| 688 | 0 | do i=G%isg,G%ieg |
| 689 | 0 | asum(i) = real_to_EFP(0.0) ; mask_sum(i) = real_to_EFP(0.0) |
| 690 | enddo | |
| 691 | ||
| 692 | 0 | do i=is,ie ; do j=js,je |
| 693 | 0 | asum(i+idg_off) = asum(i+idg_off) + real_to_EFP(scalefac*array(i,j)*mask(i,j)) |
| 694 | 0 | mask_sum(i+idg_off) = mask_sum(i+idg_off) + real_to_EFP(mask(i,j)) |
| 695 | enddo ; enddo | |
| 696 | ||
| 697 | 0 | if (query_EFP_overflow_error()) call MOM_error(FATAL, & |
| 698 | 0 | "global_j_mean overflow error occurred before sums across PEs.") |
| 699 | ||
| 700 | 0 | call EFP_sum_across_PEs(asum(G%isg:G%ieg), G%ieg-G%isg+1) |
| 701 | 0 | call EFP_sum_across_PEs(mask_sum(G%isg:G%ieg), G%ieg-G%isg+1) |
| 702 | ||
| 703 | 0 | if (query_EFP_overflow_error()) call MOM_error(FATAL, & |
| 704 | 0 | "global_j_mean overflow error occurred during sums across PEs.") |
| 705 | ||
| 706 | 0 | do i=is,ie |
| 707 | 0 | mask_sum_r = EFP_to_real(mask_sum(i+idg_off)) |
| 708 | 0 | if (mask_sum_r == 0.0 ) then ; j_mean(i) = 0.0 ; else |
| 709 | 0 | j_mean(i) = EFP_to_real(asum(i+idg_off)) / mask_sum_r |
| 710 | endif | |
| 711 | enddo | |
| 712 | ||
| 713 | 0 | deallocate(mask_sum) |
| 714 | else | |
| 715 | 0 | do i=G%isg,G%ieg ; asum(i) = real_to_EFP(0.0) ; enddo |
| 716 | ||
| 717 | 0 | do i=is,ie ; do j=js,je |
| 718 | 0 | asum(i+idg_off) = asum(i+idg_off) + real_to_EFP(scalefac*array(i,j)) |
| 719 | enddo ; enddo | |
| 720 | ||
| 721 | 0 | if (query_EFP_overflow_error()) call MOM_error(FATAL, & |
| 722 | 0 | "global_j_mean overflow error occurred before sum across PEs.") |
| 723 | ||
| 724 | 0 | call EFP_sum_across_PEs(asum(G%isg:G%ieg), G%ieg-G%isg+1) |
| 725 | ||
| 726 | 0 | if (query_EFP_overflow_error()) call MOM_error(FATAL, & |
| 727 | 0 | "global_j_mean overflow error occurred during sum across PEs.") |
| 728 | ||
| 729 | 0 | do i=is,ie |
| 730 | 0 | j_mean(i) = EFP_to_real(asum(i+idg_off)) / real(G%jeg-G%jsg+1) |
| 731 | enddo | |
| 732 | endif | |
| 733 | ||
| 734 | 0 | if (rescale /= 1.0) then ; do i=is,ie ; j_mean(i) = rescale*j_mean(i) ; enddo ; endif |
| 735 | ||
| 736 | 0 | deallocate(asum) |
| 737 | ||
| 738 | 0 | end subroutine global_j_mean |
| 739 | ||
| 740 | !> Adjust 2d array such that area mean is zero without moving the zero contour | |
| 741 | 0 | subroutine 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 | 0 | real :: posVals(G%isc:G%iec, G%jsc:G%jec) ! The positive values in a cell or 0 [A ~> a] |
| 759 | 0 | real :: negVals(G%isc:G%iec, G%jsc:G%jec) ! The negative values in a cell or 0 [A ~> a] |
| 760 | 0 | real :: areaXposVals(G%isc:G%iec, G%jsc:G%jec) ! The cell area integral of the positive values [L2 A ~> m2 a] |
| 761 | 0 | 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 | 0 | scalefac = 1.0 |
| 769 | 0 | if (present(unscale)) then ; scalefac = unscale |
| 770 | 0 | 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 | 0 | do j=G%jsc,G%jec ; do i=G%isc,G%iec |
| 776 | 0 | posVals(i,j) = max(0., array(i,j)) |
| 777 | 0 | areaXposVals(i,j) = G%areaT(i,j) * posVals(i,j) |
| 778 | 0 | negVals(i,j) = min(0., array(i,j)) |
| 779 | 0 | 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 | 0 | areaInt_EFP(1) = reproducing_sum_EFP( areaXposVals, only_on_PE=.true., unscale=scalefac*G%US%L_to_m**2 ) |
| 784 | 0 | areaInt_EFP(2) = reproducing_sum_EFP( areaXnegVals, only_on_PE=.true., unscale=scalefac*G%US%L_to_m**2 ) |
| 785 | 0 | call EFP_sum_across_PEs(areaInt_EFP, 2) |
| 786 | 0 | areaIntPosVals = EFP_to_real( areaInt_EFP(1) ) |
| 787 | 0 | areaIntNegVals = EFP_to_real( areaInt_EFP(2) ) |
| 788 | ||
| 789 | 0 | posScale = 0.0 ; negScale = 0.0 |
| 790 | 0 | if ((areaIntPosVals>0.).and.(areaIntNegVals<0.)) then ! Only adjust if possible |
| 791 | 0 | if (areaIntPosVals>-areaIntNegVals) then ! Scale down positive values |
| 792 | 0 | posScale = - areaIntNegVals / areaIntPosVals |
| 793 | 0 | do j=G%jsc,G%jec ; do i=G%isc,G%iec |
| 794 | 0 | array(i,j) = (posScale * posVals(i,j)) + negVals(i,j) |
| 795 | enddo ; enddo | |
| 796 | 0 | elseif (areaIntPosVals<-areaIntNegVals) then ! Scale down negative values |
| 797 | 0 | negScale = - areaIntPosVals / areaIntNegVals |
| 798 | 0 | do j=G%jsc,G%jec ; do i=G%isc,G%iec |
| 799 | 0 | array(i,j) = posVals(i,j) + (negScale * negVals(i,j)) |
| 800 | enddo ; enddo | |
| 801 | endif | |
| 802 | endif | |
| 803 | 0 | if (present(scaling)) scaling = posScale - negScale |
| 804 | ||
| 805 | 0 | end 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. | |
| 814 | 6 | subroutine 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 | 6 | 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 | 6 | 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 | 6 | tmax = -huge(tmax) ; tmin = huge(tmin) |
| 859 | ||
| 860 | 6 | if (find_location) then |
| 861 | ! Find the maximum and minimum tracer values on this PE and their locations. | |
| 862 | 6 | valid_PE = .false. |
| 863 | 6 | itmax = 0 ; jtmax = 0 ; ktmax = 0 ; ijk_loc_max = 0 |
| 864 | 6 | itmin = 0 ; jtmin = 0 ; ktmin = 0 ; ijk_loc_min = 0 |
| 865 | 3267456 | do k=1,nk ; do j=jsc,jec ; do i=isc,iec ; if (G%mask2dT(i,j) > 0.0) then |
| 866 | 2256300 | valid_PE = .true. |
| 867 | 2256300 | if (tr_array(i,j,k) > tmax) then |
| 868 | 2071 | tmax = tr_array(i,j,k) |
| 869 | 2071 | itmax = i ; jtmax = j ; ktmax = k |
| 870 | 2071 | ijk_loc_max = ijk_loc(i, j, k, nk, G%HI) |
| 871 | 2254229 | elseif ((tr_array(i,j,k) == tmax) .and. (k <= ktmax)) then |
| 872 | 5587 | ijk_loc_here = ijk_loc(i, j, k, nk, G%HI) |
| 873 | 5587 | if (ijk_loc_here > ijk_loc_max) then |
| 874 | 5587 | itmax = i ; jtmax = j ; ktmax = k |
| 875 | 5587 | ijk_loc_max = ijk_loc_here |
| 876 | endif | |
| 877 | endif | |
| 878 | 2256300 | if (tr_array(i,j,k) < tmin) then |
| 879 | 1596 | tmin = tr_array(i,j,k) |
| 880 | 1596 | itmin = i ; jtmin = j ; ktmin = k |
| 881 | 1596 | ijk_loc_min = ijk_loc(i, j, k, nk, G%HI) |
| 882 | 2254704 | elseif ((tr_array(i,j,k) == tmin) .and. (k <= ktmin)) then |
| 883 | 5571 | ijk_loc_here = ijk_loc(i, j, k, nk, G%HI) |
| 884 | 5571 | if (ijk_loc_here > ijk_loc_min) then |
| 885 | 5571 | itmin = i ; jtmin = j ; ktmin = k |
| 886 | 5571 | 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 | 0 | do k=1,nk ; do j=jsc,jec ; do i=isc,iec ; if (G%mask2dT(i,j) > 0.0) then |
| 893 | 0 | if (tr_array(i,j,k) > tmax) tmax = tr_array(i,j,k) |
| 894 | 0 | 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 | 6 | g_max = tmax ; g_min = tmin |
| 900 | 6 | call max_across_PEs(g_max) |
| 901 | 6 | call min_across_PEs(g_min) |
| 902 | ||
| 903 | 6 | if (find_location) then |
| 904 | 6 | if (g_max < g_min) then |
| 905 | ! This only occurs if there are no unmasked points anywhere in the domain. | |
| 906 | 0 | 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 | 6 | ijk_min_max(1:2) = 0 |
| 911 | 6 | if (valid_PE) then |
| 912 | 6 | if (g_min == tmin) ijk_min_max(1) = ijk_loc_min |
| 913 | 6 | 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 | 6 | call max_across_PEs(ijk_min_max(1)) |
| 918 | 6 | 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 | 42 | xyz_min_max(1:6) = -huge(xyz_min_max) ! These huge negative values are never selected by max_across_PEs. |
| 923 | 6 | if (valid_PE) then |
| 924 | 6 | if (ijk_min_max(1) == ijk_loc_min) then |
| 925 | 6 | xyz_min_max(1) = G%geoLonT(itmin,jtmin) |
| 926 | 6 | xyz_min_max(2) = G%geoLatT(itmin,jtmin) |
| 927 | 6 | xyz_min_max(3) = real(ktmin) |
| 928 | endif | |
| 929 | 6 | if (ijk_min_max(2) == ijk_loc_max) then |
| 930 | 6 | xyz_min_max(4) = G%geoLonT(itmax,jtmax) |
| 931 | 6 | xyz_min_max(5) = G%geoLatT(itmax,jtmax) |
| 932 | 6 | xyz_min_max(6) = real(ktmax) |
| 933 | endif | |
| 934 | endif | |
| 935 | ||
| 936 | 6 | call max_across_PEs(xyz_min_max, 6) |
| 937 | endif | |
| 938 | ||
| 939 | 6 | if (present(xgmin)) xgmin = xyz_min_max(1) |
| 940 | 6 | if (present(ygmin)) ygmin = xyz_min_max(2) |
| 941 | 6 | if (present(zgmin)) zgmin = xyz_min_max(3) |
| 942 | 6 | if (present(xgmax)) xgmax = xyz_min_max(4) |
| 943 | 6 | if (present(ygmax)) ygmax = xyz_min_max(5) |
| 944 | 6 | if (present(zgmax)) zgmax = xyz_min_max(6) |
| 945 | endif | |
| 946 | ||
| 947 | 6 | if (g_max < g_min) then |
| 948 | ! There are no unmasked points anywhere in the domain. | |
| 949 | 0 | g_max = 0.0 ; g_min = 0.0 |
| 950 | endif | |
| 951 | ||
| 952 | 6 | if (present(unscale)) then |
| 953 | ! Rescale g_min and g_max, perhaps changing their units from [CU ~> conc] to [conc] | |
| 954 | 6 | g_max = unscale * g_max |
| 955 | 6 | g_min = unscale * g_min |
| 956 | endif | |
| 957 | ||
| 958 | 6 | end subroutine array_global_min_max |
| 959 | ||
| 960 | ! Return a positive integer encoding the rotationally invariant global position of a tracer cell | |
| 961 | 14825 | function 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 | 14825 | ig = i + HI%idg_offset + (1 - HI%isg) |
| 976 | 14825 | 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 | 14825 | qturns = modulo(HI%turns, 4) |
| 980 | 14825 | if (qturns == 0) then |
| 981 | 14825 | ij_loc = ig + HI%niglobal * jg |
| 982 | 0 | elseif (qturns == 1) then |
| 983 | 0 | ij_loc = jg + HI%njglobal * ((HI%niglobal+1)-ig) |
| 984 | 0 | elseif (qturns == 2) then |
| 985 | 0 | ij_loc = ((HI%niglobal+1)-ig) + HI%niglobal * ((HI%njglobal+1)-jg) |
| 986 | 0 | elseif (qturns == 3) then |
| 987 | 0 | ij_loc = ((HI%njglobal+1)-jg) + HI%njglobal * ig |
| 988 | endif | |
| 989 | ||
| 990 | 14825 | ijk_loc = ij_loc + (HI%niglobal*HI%njglobal) * (nk-k) |
| 991 | ||
| 992 | 14825 | end function ijk_loc |
| 993 | ||
| 994 | ||
| 995 | end module MOM_spatial_means |