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