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 the ocean grid type | |
| 6 | module MOM_grid | |
| 7 | ||
| 8 | use MOM_hor_index, only : hor_index_type, hor_index_init | |
| 9 | use MOM_domains, only : MOM_domain_type, get_domain_extent, compute_block_extent | |
| 10 | use MOM_domains, only : get_global_shape, deallocate_MOM_domain | |
| 11 | use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL | |
| 12 | use MOM_file_parser, only : get_param, log_param, log_version, param_file_type | |
| 13 | use MOM_unit_scaling, only : unit_scale_type | |
| 14 | ||
| 15 | implicit none ; private | |
| 16 | ||
| 17 | #include <MOM_memory.h> | |
| 18 | ||
| 19 | public MOM_grid_init, MOM_grid_end, set_derived_metrics, set_first_direction | |
| 20 | public isPointInCell, hor_index_type, get_global_grid_size | |
| 21 | ||
| 22 | ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional | |
| 23 | ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with | |
| 24 | ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units | |
| 25 | ! vary with the Boussinesq approximation, the Boussinesq variant is given first. | |
| 26 | ||
| 27 | !> Ocean grid type. See mom_grid for details. | |
| 28 | type, public :: ocean_grid_type | |
| 29 | type(MOM_domain_type), pointer :: Domain => NULL() !< Ocean model domain | |
| 30 | type(MOM_domain_type), pointer :: Domain_aux => NULL() !< A non-symmetric auxiliary domain type. | |
| 31 | type(hor_index_type) :: HI !< Horizontal index ranges | |
| 32 | type(hor_index_type) :: HId2 !< Horizontal index ranges for level-2-downsampling | |
| 33 | ||
| 34 | integer :: isc !< The start i-index of cell centers within the computational domain | |
| 35 | integer :: iec !< The end i-index of cell centers within the computational domain | |
| 36 | integer :: jsc !< The start j-index of cell centers within the computational domain | |
| 37 | integer :: jec !< The end j-index of cell centers within the computational domain | |
| 38 | ||
| 39 | integer :: isd !< The start i-index of cell centers within the data domain | |
| 40 | integer :: ied !< The end i-index of cell centers within the data domain | |
| 41 | integer :: jsd !< The start j-index of cell centers within the data domain | |
| 42 | integer :: jed !< The end j-index of cell centers within the data domain | |
| 43 | ||
| 44 | integer :: isg !< The start i-index of cell centers within the global domain | |
| 45 | integer :: ieg !< The end i-index of cell centers within the global domain | |
| 46 | integer :: jsg !< The start j-index of cell centers within the global domain | |
| 47 | integer :: jeg !< The end j-index of cell centers within the global domain | |
| 48 | ||
| 49 | integer :: IscB !< The start i-index of cell vertices within the computational domain | |
| 50 | integer :: IecB !< The end i-index of cell vertices within the computational domain | |
| 51 | integer :: JscB !< The start j-index of cell vertices within the computational domain | |
| 52 | integer :: JecB !< The end j-index of cell vertices within the computational domain | |
| 53 | ||
| 54 | integer :: IsdB !< The start i-index of cell vertices within the data domain | |
| 55 | integer :: IedB !< The end i-index of cell vertices within the data domain | |
| 56 | integer :: JsdB !< The start j-index of cell vertices within the data domain | |
| 57 | integer :: JedB !< The end j-index of cell vertices within the data domain | |
| 58 | ||
| 59 | integer :: IsgB !< The start i-index of cell vertices within the global domain | |
| 60 | integer :: IegB !< The end i-index of cell vertices within the global domain | |
| 61 | integer :: JsgB !< The start j-index of cell vertices within the global domain | |
| 62 | integer :: JegB !< The end j-index of cell vertices within the global domain | |
| 63 | ||
| 64 | integer :: isd_global !< The value of isd in the global index space (decomposition invariant). | |
| 65 | integer :: jsd_global !< The value of isd in the global index space (decomposition invariant). | |
| 66 | integer :: idg_offset !< The offset between the corresponding global and local i-indices. | |
| 67 | integer :: jdg_offset !< The offset between the corresponding global and local j-indices. | |
| 68 | integer :: ke !< The number of layers in the vertical. | |
| 69 | logical :: symmetric !< True if symmetric memory is used. | |
| 70 | logical :: nonblocking_updates !< If true, non-blocking halo updates are | |
| 71 | !! allowed. The default is .false. (for now). | |
| 72 | integer :: first_direction !< An integer that indicates which direction is | |
| 73 | !! to be updated first in directionally split | |
| 74 | !! parts of the calculation. This can be altered | |
| 75 | !! during the course of the run via calls to | |
| 76 | !! set_first_direction. | |
| 77 | ||
| 78 | real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: & | |
| 79 | mask2dT, & !< 0 for land points and 1 for ocean points on the h-grid [nondim]. | |
| 80 | geoLatT, & !< The geographic latitude at tracer (h) points [degrees_N] or [km] or [m] | |
| 81 | geoLonT, & !< The geographic longitude at tracer (h) points [degrees_E] or [km] or [m] | |
| 82 | dxT, & !< dxT is delta x at h points [L ~> m]. | |
| 83 | IdxT, & !< 1/dxT [L-1 ~> m-1]. | |
| 84 | dyT, & !< dyT is delta y at h points [L ~> m]. | |
| 85 | IdyT, & !< IdyT is 1/dyT [L-1 ~> m-1]. | |
| 86 | areaT, & !< The area of an h-cell [L2 ~> m2]. | |
| 87 | IareaT, & !< 1/areaT [L-2 ~> m-2]. | |
| 88 | sin_rot, & !< The sine of the angular rotation between the local model grid's northward | |
| 89 | !! and the true northward directions [nondim]. | |
| 90 | cos_rot !< The cosine of the angular rotation between the local model grid's northward | |
| 91 | !! and the true northward directions [nondim]. | |
| 92 | ||
| 93 | real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_) :: & | |
| 94 | mask2dCu, & !< 0 for boundary points and 1 for ocean points on the u grid [nondim]. | |
| 95 | OBCmaskCu, & !< 0 for boundary or OBC points and 1 for ocean points on the u grid [nondim]. | |
| 96 | geoLatCu, & !< The geographic latitude at u points [degrees_N] or [km] or [m] | |
| 97 | geoLonCu, & !< The geographic longitude at u points [degrees_E] or [km] or [m]. | |
| 98 | dxCu, & !< dxCu is delta x at u points [L ~> m]. | |
| 99 | IdxCu, & !< 1/dxCu [L-1 ~> m-1]. | |
| 100 | IdxCu_OBCmask, & !< 1/dxCu or 0 at boundary or OBC points [L-1 ~> m-1]. | |
| 101 | dyCu, & !< dyCu is delta y at u points [L ~> m]. | |
| 102 | IdyCu, & !< 1/dyCu [L-1 ~> m-1]. | |
| 103 | dy_Cu, & !< The unblocked lengths of the u-faces of the h-cell [L ~> m]. | |
| 104 | IareaCu, & !< The masked inverse areas of u-grid cells [L-2 ~> m-2]. | |
| 105 | areaCu !< The areas of the u-grid cells [L2 ~> m2]. | |
| 106 | ||
| 107 | real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_) :: & | |
| 108 | mask2dCv, & !< 0 for boundary points and 1 for ocean points on the v grid [nondim]. | |
| 109 | OBCmaskCv, & !< 0 for boundary or OBC points and 1 for ocean points on the v grid [nondim]. | |
| 110 | geoLatCv, & !< The geographic latitude at v points [degrees_N] or [km] or [m] | |
| 111 | geoLonCv, & !< The geographic longitude at v points [degrees_E] or [km] or [m]. | |
| 112 | dxCv, & !< dxCv is delta x at v points [L ~> m]. | |
| 113 | IdxCv, & !< 1/dxCv [L-1 ~> m-1]. | |
| 114 | dyCv, & !< dyCv is delta y at v points [L ~> m]. | |
| 115 | IdyCv, & !< 1/dyCv [L-1 ~> m-1]. | |
| 116 | IdyCv_OBCmask, & !< 1/dxCv or 0 at boundary or OBC points [L-1 ~> m-1]. | |
| 117 | dx_Cv, & !< The unblocked lengths of the v-faces of the h-cell [L ~> m]. | |
| 118 | IareaCv, & !< The masked inverse areas of v-grid cells [L-2 ~> m-2]. | |
| 119 | areaCv !< The areas of the v-grid cells [L2 ~> m2]. | |
| 120 | ||
| 121 | real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_) :: & | |
| 122 | porous_DminU, & !< minimum topographic height (deepest) of U-face [Z ~> m] | |
| 123 | porous_DmaxU, & !< maximum topographic height (shallowest) of U-face [Z ~> m] | |
| 124 | porous_DavgU !< average topographic height of U-face [Z ~> m] | |
| 125 | ||
| 126 | real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_) :: & | |
| 127 | porous_DminV, & !< minimum topographic height (deepest) of V-face [Z ~> m] | |
| 128 | porous_DmaxV, & !< maximum topographic height (shallowest) of V-face [Z ~> m] | |
| 129 | porous_DavgV !< average topographic height of V-face [Z ~> m] | |
| 130 | ||
| 131 | real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: & | |
| 132 | mask2dBu, & !< 0 for boundary points and 1 for ocean points on the q grid [nondim]. | |
| 133 | geoLatBu, & !< The geographic latitude at q points [degrees_N] or [km] or [m] | |
| 134 | geoLonBu, & !< The geographic longitude at q points [degrees_E] or [km] or [m]. | |
| 135 | dxBu, & !< dxBu is delta x at q points [L ~> m]. | |
| 136 | IdxBu, & !< 1/dxBu [L-1 ~> m-1]. | |
| 137 | dyBu, & !< dyBu is delta y at q points [L ~> m]. | |
| 138 | IdyBu, & !< 1/dyBu [L-1 ~> m-1]. | |
| 139 | areaBu, & !< areaBu is the area of a q-cell [L2 ~> m2] | |
| 140 | IareaBu !< IareaBu = 1/areaBu [L-2 ~> m-2]. | |
| 141 | ||
| 142 | real, pointer, dimension(:) :: & | |
| 143 | gridLatT => NULL(), & !< The latitude of T points for the purpose of labeling the output axes, | |
| 144 | !! often in units of [degrees_N] or [km] or [m] or [gridpoints]. | |
| 145 | !! On many grids this is the same as geoLatT. | |
| 146 | gridLatB => NULL() !< The latitude of B points for the purpose of labeling the output axes, | |
| 147 | !! often in units of [degrees_N] or [km] or [m] or [gridpoints]. | |
| 148 | !! On many grids this is the same as geoLatBu. | |
| 149 | real, pointer, dimension(:) :: & | |
| 150 | gridLonT => NULL(), & !< The longitude of T points for the purpose of labeling the output axes, | |
| 151 | !! often in units of [degrees_E] or [km] or [m] or [gridpoints]. | |
| 152 | !! On many grids this is the same as geoLonT. | |
| 153 | gridLonB => NULL() !< The longitude of B points for the purpose of labeling the output axes, | |
| 154 | !! often in units of [degrees_E] or [km] or [m] or [gridpoints]. | |
| 155 | !! On many grids this is the same as geoLonBu. | |
| 156 | character(len=40) :: & | |
| 157 | ! Except on a Cartesian grid, these are usually some variant of "degrees". | |
| 158 | x_axis_units, & !< The units that are used in labeling the x coordinate axes. | |
| 159 | y_axis_units, & !< The units that are used in labeling the y coordinate axes. | |
| 160 | ! These are internally generated names, including "m", "km", "deg_E" and "deg_N". | |
| 161 | x_ax_unit_short, & !< A short description of the x-axis units for documenting parameter units | |
| 162 | y_ax_unit_short !< A short description of the y-axis units for documenting parameter units | |
| 163 | ||
| 164 | real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: & | |
| 165 | bathyT !< Ocean bottom depth, referenced to Z_ref at tracer points. bathyT is in | |
| 166 | !! depth units and positive *below* Z_ref [Z ~> m]. | |
| 167 | real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: & | |
| 168 | meanSL !< Spatially varying time mean sea level, referenced to Z_ref at tracer points. | |
| 169 | !! meanSL is in height units and positive *above* Z_ref. It is used | |
| 170 | !! a) as the height where p = p_atm or zero; | |
| 171 | !! b) to calculate time mean thickness of the water column, where | |
| 172 | !! mean thickness = max(meanSL + bathyT, 0.0). | |
| 173 | !! meanSL is 2D for the consideration of a domain with spatically varying mean | |
| 174 | !! height, e.g. the Great Lakes system [Z ~> m]. | |
| 175 | real :: Z_ref !< A reference value for all geometric height fields, such as bathyT [Z ~> m]. | |
| 176 | ||
| 177 | logical :: bathymetry_at_vel !< If true, there are separate values for the | |
| 178 | !! basin depths at velocity points. Otherwise the effects of | |
| 179 | !! of topography are entirely determined from thickness points. | |
| 180 | real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_) :: & | |
| 181 | Dblock_u, & !< Topographic depths at u-points at which the flow is blocked [Z ~> m]. | |
| 182 | Dopen_u !< Topographic depths at u-points at which the flow is open at width dy_Cu [Z ~> m]. | |
| 183 | real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_) :: & | |
| 184 | Dblock_v, & !< Topographic depths at v-points at which the flow is blocked [Z ~> m]. | |
| 185 | Dopen_v !< Topographic depths at v-points at which the flow is open at width dx_Cv [Z ~> m]. | |
| 186 | real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: & | |
| 187 | CoriolisBu, & !< The Coriolis parameter at corner points [T-1 ~> s-1]. | |
| 188 | Coriolis2Bu !< The square of the Coriolis parameter at corner points [T-2 ~> s-2]. | |
| 189 | real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: & | |
| 190 | df_dx, & !< Derivative d/dx f (Coriolis parameter) at h-points [T-1 L-1 ~> s-1 m-1]. | |
| 191 | df_dy !< Derivative d/dy f (Coriolis parameter) at h-points [T-1 L-1 ~> s-1 m-1]. | |
| 192 | ||
| 193 | ! These variables are global sums that are useful for 1-d diagnostics. | |
| 194 | real :: areaT_global !< Global sum of h-cell area [L2 ~> m2] | |
| 195 | real :: IareaT_global !< Global sum of inverse h-cell area (1/areaT_global) [L-2 ~> m-2]. | |
| 196 | ||
| 197 | type(unit_scale_type), pointer :: US => NULL() !< A dimensional unit scaling type | |
| 198 | ||
| 199 | ||
| 200 | ! These variables are for block structures. | |
| 201 | integer :: nblocks !< The number of sub-PE blocks on this PE | |
| 202 | type(hor_index_type), pointer :: Block(:) => NULL() !< Index ranges for each block | |
| 203 | ||
| 204 | ! These parameters are run-time parameters that are used during some | |
| 205 | ! initialization routines (but not all) | |
| 206 | real :: grid_unit_to_L !< A factor that converts a the geoLat and geoLon variables and related | |
| 207 | !! variables like len_lat and len_lon into rescaled horizontal distance | |
| 208 | !! units on a Cartesian grid, in [L km ~> 1000] or [L m-1 ~> 1] or | |
| 209 | !! is 0 for a non-Cartesian grid. | |
| 210 | real :: south_lat !< The latitude (or y-coordinate) of the first v-line [degrees_N] or [km] or [m] | |
| 211 | real :: west_lon !< The longitude (or x-coordinate) of the first u-line [degrees_E] or [km] or [m] | |
| 212 | real :: len_lat !< The latitudinal (or y-coord) extent of physical domain [degrees_N] or [km] or [m] | |
| 213 | real :: len_lon !< The longitudinal (or x-coord) extent of physical domain [degrees_E] or [km] or [m] | |
| 214 | real :: Rad_Earth_L !< The radius of the planet in rescaled units [L ~> m] | |
| 215 | real :: max_depth !< The maximum depth of the ocean in depth units [Z ~> m] | |
| 216 | end type ocean_grid_type | |
| 217 | ||
| 218 | contains | |
| 219 | ||
| 220 | !> MOM_grid_init initializes the ocean grid array sizes and grid memory. | |
| 221 | 1 | subroutine MOM_grid_init(G, param_file, US, HI, global_indexing, bathymetry_at_vel) |
| 222 | type(ocean_grid_type), intent(inout) :: G !< The horizontal grid type | |
| 223 | type(param_file_type), intent(in) :: param_file !< Parameter file handle | |
| 224 | type(unit_scale_type), optional, pointer :: US !< A dimensional unit scaling type | |
| 225 | type(hor_index_type), & | |
| 226 | optional, intent(in) :: HI !< A hor_index_type for array extents | |
| 227 | logical, optional, intent(in) :: global_indexing !< If true use global index | |
| 228 | !! values instead of having the data domain on each | |
| 229 | !! processor start at 1. | |
| 230 | logical, optional, intent(in) :: bathymetry_at_vel !< If true, there are | |
| 231 | !! separate values for the ocean bottom depths at | |
| 232 | !! velocity points. Otherwise the effects of topography | |
| 233 | !! are entirely determined from thickness points. | |
| 234 | ||
| 235 | ! Local variables | |
| 236 | real :: mean_SeaLev_scale ! A scaling factor for the reference height variable [1] or [Z m-1 ~> 1] | |
| 237 | integer :: isd, ied, jsd, jed | |
| 238 | integer :: IsdB, IedB, JsdB, JedB | |
| 239 | integer :: ied_max, jed_max | |
| 240 | integer :: niblock, njblock, nihalo, njhalo, nblocks, n, i, j | |
| 241 | logical :: local_indexing ! If false use global index values instead of having | |
| 242 | ! the data domain on each processor start at 1. | |
| 243 | ! This include declares and sets the variable "version". | |
| 244 | # include "version_variable.h" | |
| 245 | ||
| 246 | 1 | integer, allocatable, dimension(:) :: ibegin, iend, jbegin, jend |
| 247 | character(len=40) :: mod_nm = "MOM_grid" ! This module's name. | |
| 248 | ||
| 249 | 1 | mean_SeaLev_scale = 1.0 ; if (associated(G%US)) mean_SeaLev_scale = G%US%m_to_Z |
| 250 | ||
| 251 | ! Read all relevant parameters and write them to the model log. | |
| 252 | call get_param(param_file, mod_nm, "REFERENCE_HEIGHT", G%Z_ref, & | |
| 253 | 1 | units="m", default=0.0, scale=mean_SeaLev_scale, do_not_log=.true.) |
| 254 | call log_version(param_file, mod_nm, version, & | |
| 255 | "Parameters providing information about the lateral grid.", & | |
| 256 | 1 | log_to_all=.true., layout=.true., all_default=(G%Z_ref==0.0)) |
| 257 | ||
| 258 | call get_param(param_file, mod_nm, "NIBLOCK", niblock, "The number of blocks "// & | |
| 259 | "in the x-direction on each processor (for openmp).", default=1, & | |
| 260 | 1 | layoutParam=.true.) |
| 261 | call get_param(param_file, mod_nm, "NJBLOCK", njblock, "The number of blocks "// & | |
| 262 | "in the y-direction on each processor (for openmp).", default=1, & | |
| 263 | 1 | layoutParam=.true.) |
| 264 | 1 | if (present(US)) then ; if (associated(US)) G%US => US ; endif |
| 265 | ||
| 266 | call get_param(param_file, mod_nm, "REFERENCE_HEIGHT", G%Z_ref, & | |
| 267 | "A reference value for geometric height fields, such as bathyT.", & | |
| 268 | 1 | units="m", default=0.0, scale=mean_SeaLev_scale) |
| 269 | ||
| 270 | 1 | if (present(HI)) then |
| 271 | 1 | G%HI = HI |
| 272 | ||
| 273 | 1 | G%isc = HI%isc ; G%iec = HI%iec ; G%jsc = HI%jsc ; G%jec = HI%jec |
| 274 | 1 | G%isd = HI%isd ; G%ied = HI%ied ; G%jsd = HI%jsd ; G%jed = HI%jed |
| 275 | 1 | G%isg = HI%isg ; G%ieg = HI%ieg ; G%jsg = HI%jsg ; G%jeg = HI%jeg |
| 276 | ||
| 277 | 1 | G%IscB = HI%IscB ; G%IecB = HI%IecB ; G%JscB = HI%JscB ; G%JecB = HI%JecB |
| 278 | 1 | G%IsdB = HI%IsdB ; G%IedB = HI%IedB ; G%JsdB = HI%JsdB ; G%JedB = HI%JedB |
| 279 | 1 | G%IsgB = HI%IsgB ; G%IegB = HI%IegB ; G%JsgB = HI%JsgB ; G%JegB = HI%JegB |
| 280 | ||
| 281 | 1 | G%idg_offset = HI%idg_offset ; G%jdg_offset = HI%jdg_offset |
| 282 | 1 | G%isd_global = G%isd + HI%idg_offset ; G%jsd_global = G%jsd + HI%jdg_offset |
| 283 | 1 | G%symmetric = HI%symmetric |
| 284 | else | |
| 285 | 0 | local_indexing = .true. |
| 286 | 0 | if (present(global_indexing)) local_indexing = .not.global_indexing |
| 287 | call hor_index_init(G%Domain, G%HI, param_file, & | |
| 288 | 0 | local_indexing=local_indexing) |
| 289 | ||
| 290 | ! get_domain_extent ensures that domains start at 1 for compatibility between | |
| 291 | ! static and dynamically allocated arrays, unless global_indexing is true. | |
| 292 | call get_domain_extent(G%Domain, G%isc, G%iec, G%jsc, G%jec, & | |
| 293 | G%isd, G%ied, G%jsd, G%jed, & | |
| 294 | G%isg, G%ieg, G%jsg, G%jeg, & | |
| 295 | G%idg_offset, G%jdg_offset, G%symmetric, & | |
| 296 | 0 | local_indexing=local_indexing) |
| 297 | 0 | G%isd_global = G%isd+G%idg_offset ; G%jsd_global = G%jsd+G%jdg_offset |
| 298 | endif | |
| 299 | ||
| 300 | 1 | G%nonblocking_updates = G%Domain%nonblocking_updates |
| 301 | ||
| 302 | ! Set array sizes for fields that are discretized at tracer cell boundaries. | |
| 303 | 1 | G%IscB = G%isc ; G%JscB = G%jsc |
| 304 | 1 | G%IsdB = G%isd ; G%JsdB = G%jsd |
| 305 | 1 | G%IsgB = G%isg ; G%JsgB = G%jsg |
| 306 | 1 | if (G%symmetric) then |
| 307 | 1 | G%IscB = G%isc-1 ; G%JscB = G%jsc-1 |
| 308 | 1 | G%IsdB = G%isd-1 ; G%JsdB = G%jsd-1 |
| 309 | 1 | G%IsgB = G%isg-1 ; G%JsgB = G%jsg-1 |
| 310 | endif | |
| 311 | 1 | G%IecB = G%iec ; G%JecB = G%jec |
| 312 | 1 | G%IedB = G%ied ; G%JedB = G%jed |
| 313 | 1 | G%IegB = G%ieg ; G%JegB = G%jeg |
| 314 | ||
| 315 | 1 | call MOM_mesg(" MOM_grid.F90, MOM_grid_init: allocating metrics", 5) |
| 316 | ||
| 317 | 1 | call allocate_metrics(G) |
| 318 | ||
| 319 | 1 | isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed |
| 320 | 1 | IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB |
| 321 | ||
| 322 | 1 | G%bathymetry_at_vel = .false. |
| 323 | 1 | if (present(bathymetry_at_vel)) G%bathymetry_at_vel = bathymetry_at_vel |
| 324 | 1 | if (G%bathymetry_at_vel) then |
| 325 | 0 | ALLOC_(G%Dblock_u(IsdB:IedB, jsd:jed)) ; G%Dblock_u(:,:) = -G%Z_ref |
| 326 | 0 | ALLOC_(G%Dopen_u(IsdB:IedB, jsd:jed)) ; G%Dopen_u(:,:) = -G%Z_ref |
| 327 | 0 | ALLOC_(G%Dblock_v(isd:ied, JsdB:JedB)) ; G%Dblock_v(:,:) = -G%Z_ref |
| 328 | 0 | ALLOC_(G%Dopen_v(isd:ied, JsdB:JedB)) ; G%Dopen_v(:,:) = -G%Z_ref |
| 329 | endif | |
| 330 | ||
| 331 | ! setup block indices. | |
| 332 | 1 | nihalo = G%Domain%nihalo |
| 333 | 1 | njhalo = G%Domain%njhalo |
| 334 | 1 | nblocks = niblock * njblock |
| 335 | 1 | if (nblocks < 1) call MOM_error(FATAL, "MOM_grid_init: " // & |
| 336 | 0 | "nblocks(=NI_BLOCK*NJ_BLOCK) must be no less than 1") |
| 337 | ||
| 338 | 1 | allocate(ibegin(niblock), iend(niblock), jbegin(njblock), jend(njblock)) |
| 339 | 1 | call compute_block_extent(G%HI%isc,G%HI%iec,niblock,ibegin,iend) |
| 340 | 1 | call compute_block_extent(G%HI%jsc,G%HI%jec,njblock,jbegin,jend) |
| 341 | !-- make sure the last block is the largest. | |
| 342 | 1 | do i = 1, niblock-1 |
| 343 | 0 | if (iend(i)-ibegin(i) > iend(niblock)-ibegin(niblock) ) call MOM_error(FATAL, & |
| 344 | 1 | "MOM_grid_init: the last block size in x-direction is not the largest") |
| 345 | enddo | |
| 346 | 1 | do j = 1, njblock-1 |
| 347 | 0 | if (jend(j)-jbegin(j) > jend(njblock)-jbegin(njblock) ) call MOM_error(FATAL, & |
| 348 | 1 | "MOM_grid_init: the last block size in y-direction is not the largest") |
| 349 | enddo | |
| 350 | ||
| 351 | 1 | G%nblocks = nblocks |
| 352 | 1 | allocate(G%Block(nblocks)) |
| 353 | 1 | ied_max = 1 ; jed_max = 1 |
| 354 | 2 | do n = 1,nblocks |
| 355 | ! Copy all information from the array index type describing the local grid. | |
| 356 | 1 | G%Block(n) = G%HI |
| 357 | ||
| 358 | 1 | i = mod((n-1), niblock) + 1 |
| 359 | 1 | j = (n-1)/niblock + 1 |
| 360 | !--- isd and jsd are always 1 for each block to permit array reuse. | |
| 361 | 1 | G%Block(n)%isd = 1 ; G%Block(n)%jsd = 1 |
| 362 | 1 | G%Block(n)%isc = G%Block(n)%isd+nihalo |
| 363 | 1 | G%Block(n)%jsc = G%Block(n)%jsd+njhalo |
| 364 | 1 | G%Block(n)%iec = G%Block(n)%isc + iend(i) - ibegin(i) |
| 365 | 1 | G%Block(n)%jec = G%Block(n)%jsc + jend(j) - jbegin(j) |
| 366 | 1 | G%Block(n)%ied = G%Block(n)%iec + nihalo |
| 367 | 1 | G%Block(n)%jed = G%Block(n)%jec + njhalo |
| 368 | 1 | G%Block(n)%IscB = G%Block(n)%isc ; G%Block(n)%IecB = G%Block(n)%iec |
| 369 | 1 | G%Block(n)%JscB = G%Block(n)%jsc ; G%Block(n)%JecB = G%Block(n)%jec |
| 370 | ! For symmetric memory domains, the first block will have the extra point | |
| 371 | ! at the lower boundary of its computational domain. | |
| 372 | 1 | if (G%symmetric) then |
| 373 | 1 | if (i==1) G%Block(n)%IscB = G%Block(n)%IscB-1 |
| 374 | 1 | if (j==1) G%Block(n)%JscB = G%Block(n)%JscB-1 |
| 375 | endif | |
| 376 | 1 | G%Block(n)%IsdB = G%Block(n)%isd ; G%Block(n)%IedB = G%Block(n)%ied |
| 377 | 1 | G%Block(n)%JsdB = G%Block(n)%jsd ; G%Block(n)%JedB = G%Block(n)%jed |
| 378 | !--- For symmetric memory domain, every block will have an extra point | |
| 379 | !--- at the lower boundary of its data domain. | |
| 380 | 1 | if (G%symmetric) then |
| 381 | 1 | G%Block(n)%IsdB = G%Block(n)%IsdB-1 |
| 382 | 1 | G%Block(n)%JsdB = G%Block(n)%JsdB-1 |
| 383 | endif | |
| 384 | 1 | G%Block(n)%idg_offset = (ibegin(i) - G%Block(n)%isc) + G%HI%idg_offset |
| 385 | 1 | G%Block(n)%jdg_offset = (jbegin(j) - G%Block(n)%jsc) + G%HI%jdg_offset |
| 386 | ! Find the largest values of ied and jed so that all blocks will have the | |
| 387 | ! same size in memory. | |
| 388 | 1 | ied_max = max(ied_max, G%Block(n)%ied) |
| 389 | 2 | jed_max = max(jed_max, G%Block(n)%jed) |
| 390 | enddo | |
| 391 | ||
| 392 | ! Reset all of the data domain sizes to match the largest for array reuse, | |
| 393 | ! recalling that all block have isd=jed=1 for array reuse. | |
| 394 | 2 | do n = 1,nblocks |
| 395 | 1 | G%Block(n)%ied = ied_max ; G%Block(n)%IedB = ied_max |
| 396 | 2 | G%Block(n)%jed = jed_max ; G%Block(n)%JedB = jed_max |
| 397 | enddo | |
| 398 | ||
| 399 | !-- do some bounds error checking | |
| 400 | 1 | if ( G%block(nblocks)%ied+G%block(nblocks)%idg_offset > G%HI%ied + G%HI%idg_offset ) & |
| 401 | 0 | call MOM_error(FATAL, "MOM_grid_init: G%ied_bk > G%ied") |
| 402 | 1 | if ( G%block(nblocks)%jed+G%block(nblocks)%jdg_offset > G%HI%jed + G%HI%jdg_offset ) & |
| 403 | 0 | call MOM_error(FATAL, "MOM_grid_init: G%jed_bk > G%jed") |
| 404 | ||
| 405 | call get_domain_extent(G%Domain, G%HId2%isc, G%HId2%iec, G%HId2%jsc, G%HId2%jec, & | |
| 406 | G%HId2%isd, G%HId2%ied, G%HId2%jsd, G%HId2%jed, & | |
| 407 | 1 | G%HId2%isg, G%HId2%ieg, G%HId2%jsg, G%HId2%jeg, coarsen=2) |
| 408 | ||
| 409 | ! Set array sizes for fields that are discretized at tracer cell boundaries. | |
| 410 | 1 | G%HId2%IscB = G%HId2%isc ; G%HId2%JscB = G%HId2%jsc |
| 411 | 1 | G%HId2%IsdB = G%HId2%isd ; G%HId2%JsdB = G%HId2%jsd |
| 412 | 1 | G%HId2%IsgB = G%HId2%isg ; G%HId2%JsgB = G%HId2%jsg |
| 413 | 1 | if (G%symmetric) then |
| 414 | 1 | G%HId2%IscB = G%HId2%isc-1 ; G%HId2%JscB = G%HId2%jsc-1 |
| 415 | 1 | G%HId2%IsdB = G%HId2%isd-1 ; G%HId2%JsdB = G%HId2%jsd-1 |
| 416 | 1 | G%HId2%IsgB = G%HId2%isg-1 ; G%HId2%JsgB = G%HId2%jsg-1 |
| 417 | endif | |
| 418 | 1 | G%HId2%IecB = G%HId2%iec ; G%HId2%JecB = G%HId2%jec |
| 419 | 1 | G%HId2%IedB = G%HId2%ied ; G%HId2%JedB = G%HId2%jed |
| 420 | 1 | G%HId2%IegB = G%HId2%ieg ; G%HId2%JegB = G%HId2%jeg |
| 421 | ||
| 422 | 1 | end subroutine MOM_grid_init |
| 423 | ||
| 424 | !> set_derived_metrics calculates metric terms that are derived from other metrics. | |
| 425 | 1 | subroutine set_derived_metrics(G, US) |
| 426 | type(ocean_grid_type), intent(inout) :: G !< The horizontal grid structure | |
| 427 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 428 | ! Various inverse grid spacings and derived areas are calculated within this | |
| 429 | ! subroutine. | |
| 430 | integer :: i, j, isd, ied, jsd, jed | |
| 431 | integer :: IsdB, IedB, JsdB, JedB | |
| 432 | ||
| 433 | 1 | isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed |
| 434 | 1 | IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB |
| 435 | ||
| 436 | 8773 | do j=jsd,jed ; do i=isd,ied |
| 437 | 8704 | if (G%dxT(i,j) < 0.0) G%dxT(i,j) = 0.0 |
| 438 | 8704 | if (G%dyT(i,j) < 0.0) G%dyT(i,j) = 0.0 |
| 439 | 8704 | G%IdxT(i,j) = Adcroft_reciprocal(G%dxT(i,j)) |
| 440 | 8704 | G%IdyT(i,j) = Adcroft_reciprocal(G%dyT(i,j)) |
| 441 | 8772 | G%IareaT(i,j) = Adcroft_reciprocal(G%areaT(i,j)) |
| 442 | enddo ; enddo | |
| 443 | ||
| 444 | 8841 | do j=jsd,jed ; do I=IsdB,IedB |
| 445 | 8772 | if (G%dxCu(I,j) < 0.0) G%dxCu(I,j) = 0.0 |
| 446 | 8772 | if (G%dyCu(I,j) < 0.0) G%dyCu(I,j) = 0.0 |
| 447 | 8772 | G%IdxCu(I,j) = Adcroft_reciprocal(G%dxCu(I,j)) |
| 448 | 8772 | G%IdyCu(I,j) = Adcroft_reciprocal(G%dyCu(I,j)) |
| 449 | 8840 | G%IdxCu_OBCmask(I,j) = G%OBCmaskCu(I,j) * G%IdxCu(I,j) ! This may be reset if masks are reset. |
| 450 | enddo ; enddo | |
| 451 | ||
| 452 | 8902 | do J=JsdB,JedB ; do i=isd,ied |
| 453 | 8832 | if (G%dxCv(i,J) < 0.0) G%dxCv(i,J) = 0.0 |
| 454 | 8832 | if (G%dyCv(i,J) < 0.0) G%dyCv(i,J) = 0.0 |
| 455 | 8832 | G%IdxCv(i,J) = Adcroft_reciprocal(G%dxCv(i,J)) |
| 456 | 8832 | G%IdyCv(i,J) = Adcroft_reciprocal(G%dyCv(i,J)) |
| 457 | 8901 | G%IdyCv_OBCmask(i,J) = G%OBCmaskCv(i,J) * G%IdyCv(i,J) ! This may be reset if masks are reset. |
| 458 | enddo ; enddo | |
| 459 | ||
| 460 | 8971 | do J=JsdB,JedB ; do I=IsdB,IedB |
| 461 | 8901 | if (G%dxBu(I,J) < 0.0) G%dxBu(I,J) = 0.0 |
| 462 | 8901 | if (G%dyBu(I,J) < 0.0) G%dyBu(I,J) = 0.0 |
| 463 | ||
| 464 | 8901 | G%IdxBu(I,J) = Adcroft_reciprocal(G%dxBu(I,J)) |
| 465 | 8901 | G%IdyBu(I,J) = Adcroft_reciprocal(G%dyBu(I,J)) |
| 466 | ! areaBu has usually been set to a positive area elsewhere. | |
| 467 | 8901 | if (G%areaBu(I,J) <= 0.0) G%areaBu(I,J) = G%dxBu(I,J) * G%dyBu(I,J) |
| 468 | 8970 | G%IareaBu(I,J) = Adcroft_reciprocal(G%areaBu(I,J)) |
| 469 | enddo ; enddo | |
| 470 | 1 | end subroutine set_derived_metrics |
| 471 | ||
| 472 | !> Adcroft_reciprocal(x) = 1/x for |x|>0 or 0 for x=0. | |
| 473 | 88023 | function Adcroft_reciprocal(val) result(I_val) |
| 474 | real, intent(in) :: val !< The value being inverted [A]. | |
| 475 | real :: I_val !< The Adcroft reciprocal of val [A-1]. | |
| 476 | ||
| 477 | 88023 | I_val = 0.0 ; if (val /= 0.0) I_val = 1.0/val |
| 478 | 88023 | end function Adcroft_reciprocal |
| 479 | ||
| 480 | !> Returns true if the coordinates (x,y) are within the h-cell (i,j) | |
| 481 | 0 | logical function isPointInCell(G, i, j, x, y) |
| 482 | type(ocean_grid_type), intent(in) :: G !< Grid type | |
| 483 | integer, intent(in) :: i !< i index of cell to test | |
| 484 | integer, intent(in) :: j !< j index of cell to test | |
| 485 | real, intent(in) :: x !< x coordinate of point [degrees_E] | |
| 486 | real, intent(in) :: y !< y coordinate of point [degrees_N] | |
| 487 | ! Local variables | |
| 488 | real :: xNE, xNW, xSE, xSW ! Longitudes of cell corners [degrees_E] | |
| 489 | real :: yNE, yNW, ySE, ySW ! Latitudes of cell corners [degrees_N] | |
| 490 | real :: l0, l1, l2, l3 ! Crossed products of differences in position [degrees_E degrees_N] | |
| 491 | real :: p0, p1, p2, p3 ! Trinary unitary values reflecting the signs of the crossed products [nondim] | |
| 492 | 0 | isPointInCell = .false. |
| 493 | 0 | xNE = G%geoLonBu(i ,j ) ; yNE = G%geoLatBu(i ,j ) |
| 494 | 0 | xNW = G%geoLonBu(i-1,j ) ; yNW = G%geoLatBu(i-1,j ) |
| 495 | 0 | xSE = G%geoLonBu(i ,j-1) ; ySE = G%geoLatBu(i ,j-1) |
| 496 | 0 | xSW = G%geoLonBu(i-1,j-1) ; ySW = G%geoLatBu(i-1,j-1) |
| 497 | ! This is a crude calculation that assumes a geographic coordinate system | |
| 498 | if (x<min(xNE,xNW,xSE,xSW) .or. x>max(xNE,xNW,xSE,xSW) .or. & | |
| 499 | 0 | y<min(yNE,yNW,ySE,ySW) .or. y>max(yNE,yNW,ySE,ySW) ) then |
| 500 | 0 | return ! Avoid the more complicated calculation |
| 501 | endif | |
| 502 | 0 | l0 = (x-xSW)*(ySE-ySW) - (y-ySW)*(xSE-xSW) |
| 503 | 0 | l1 = (x-xSE)*(yNE-ySE) - (y-ySE)*(xNE-xSE) |
| 504 | 0 | l2 = (x-xNE)*(yNW-yNE) - (y-yNE)*(xNW-xNE) |
| 505 | 0 | l3 = (x-xNW)*(ySW-yNW) - (y-yNW)*(xSW-xNW) |
| 506 | ||
| 507 | 0 | p0 = sign(1., l0) ; if (l0 == 0.) p0=0. |
| 508 | 0 | p1 = sign(1., l1) ; if (l1 == 0.) p1=0. |
| 509 | 0 | p2 = sign(1., l2) ; if (l2 == 0.) p2=0. |
| 510 | 0 | p3 = sign(1., l3) ; if (l3 == 0.) p3=0. |
| 511 | ||
| 512 | 0 | if ( (abs(p0)+abs(p2)) + (abs(p1)+abs(p3)) == abs((p0+p2) + (p1+p3)) ) then |
| 513 | 0 | isPointInCell=.true. |
| 514 | endif | |
| 515 | 0 | end function isPointInCell |
| 516 | ||
| 517 | !> Store an integer indicating which direction to work on first. | |
| 518 | 1 | subroutine set_first_direction(G, y_first) |
| 519 | type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure | |
| 520 | integer, intent(in) :: y_first !< The first direction to store | |
| 521 | ||
| 522 | 1 | G%first_direction = y_first |
| 523 | 1 | end subroutine set_first_direction |
| 524 | ||
| 525 | !> Return global shape of horizontal grid | |
| 526 | 0 | subroutine get_global_grid_size(G, niglobal, njglobal) |
| 527 | type(ocean_grid_type), intent(inout) :: G !< The horizontal grid type | |
| 528 | integer, intent(out) :: niglobal !< i-index global size of grid | |
| 529 | integer, intent(out) :: njglobal !< j-index global size of grid | |
| 530 | ||
| 531 | 0 | call get_global_shape(G%domain, niglobal, njglobal) |
| 532 | ||
| 533 | 0 | end subroutine get_global_grid_size |
| 534 | ||
| 535 | !> Allocate memory used by the ocean_grid_type and related structures. | |
| 536 | 1 | subroutine allocate_metrics(G) |
| 537 | type(ocean_grid_type), intent(inout) :: G !< The horizontal grid type | |
| 538 | integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, isg, ieg, jsg, jeg | |
| 539 | ||
| 540 | ! This subroutine allocates the lateral elements of the ocean_grid_type that | |
| 541 | ! are always used and zeros them out. | |
| 542 | ||
| 543 | 1 | isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed |
| 544 | 1 | IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB |
| 545 | 1 | isg = G%isg ; ieg = G%ieg ; jsg = G%jsg ; jeg = G%jeg |
| 546 | ||
| 547 | 8773 | ALLOC_(G%dxT(isd:ied,jsd:jed)) ; G%dxT(:,:) = 0.0 |
| 548 | 8841 | ALLOC_(G%dxCu(IsdB:IedB,jsd:jed)) ; G%dxCu(:,:) = 0.0 |
| 549 | 8902 | ALLOC_(G%dxCv(isd:ied,JsdB:JedB)) ; G%dxCv(:,:) = 0.0 |
| 550 | 8971 | ALLOC_(G%dxBu(IsdB:IedB,JsdB:JedB)) ; G%dxBu(:,:) = 0.0 |
| 551 | 8773 | ALLOC_(G%IdxT(isd:ied,jsd:jed)) ; G%IdxT(:,:) = 0.0 |
| 552 | 8841 | ALLOC_(G%IdxCu(IsdB:IedB,jsd:jed)) ; G%IdxCu(:,:) = 0.0 |
| 553 | 8841 | ALLOC_(G%IdxCu_OBCmask(IsdB:IedB,jsd:jed)) ; G%IdxCu_OBCmask(:,:) = 0.0 |
| 554 | 8902 | ALLOC_(G%IdxCv(isd:ied,JsdB:JedB)) ; G%IdxCv(:,:) = 0.0 |
| 555 | 8971 | ALLOC_(G%IdxBu(IsdB:IedB,JsdB:JedB)) ; G%IdxBu(:,:) = 0.0 |
| 556 | ||
| 557 | 8773 | ALLOC_(G%dyT(isd:ied,jsd:jed)) ; G%dyT(:,:) = 0.0 |
| 558 | 8841 | ALLOC_(G%dyCu(IsdB:IedB,jsd:jed)) ; G%dyCu(:,:) = 0.0 |
| 559 | 8902 | ALLOC_(G%dyCv(isd:ied,JsdB:JedB)) ; G%dyCv(:,:) = 0.0 |
| 560 | 8971 | ALLOC_(G%dyBu(IsdB:IedB,JsdB:JedB)) ; G%dyBu(:,:) = 0.0 |
| 561 | 8773 | ALLOC_(G%IdyT(isd:ied,jsd:jed)) ; G%IdyT(:,:) = 0.0 |
| 562 | 8841 | ALLOC_(G%IdyCu(IsdB:IedB,jsd:jed)) ; G%IdyCu(:,:) = 0.0 |
| 563 | 8902 | ALLOC_(G%IdyCv(isd:ied,JsdB:JedB)) ; G%IdyCv(:,:) = 0.0 |
| 564 | 8902 | ALLOC_(G%IdyCv_OBCmask(isd:ied,JsdB:JedB)) ; G%IdyCv_OBCmask(:,:) = 0.0 |
| 565 | 8971 | ALLOC_(G%IdyBu(IsdB:IedB,JsdB:JedB)) ; G%IdyBu(:,:) = 0.0 |
| 566 | ||
| 567 | 8773 | ALLOC_(G%areaT(isd:ied,jsd:jed)) ; G%areaT(:,:) = 0.0 |
| 568 | 8773 | ALLOC_(G%IareaT(isd:ied,jsd:jed)) ; G%IareaT(:,:) = 0.0 |
| 569 | 8971 | ALLOC_(G%areaBu(IsdB:IedB,JsdB:JedB)) ; G%areaBu(:,:) = 0.0 |
| 570 | 8971 | ALLOC_(G%IareaBu(IsdB:IedB,JsdB:JedB)) ; G%IareaBu(:,:) = 0.0 |
| 571 | ||
| 572 | 8773 | ALLOC_(G%mask2dT(isd:ied,jsd:jed)) ; G%mask2dT(:,:) = 0.0 |
| 573 | 8841 | ALLOC_(G%mask2dCu(IsdB:IedB,jsd:jed)) ; G%mask2dCu(:,:) = 0.0 |
| 574 | 8841 | ALLOC_(G%OBCmaskCu(IsdB:IedB,jsd:jed)) ; G%OBCmaskCu(:,:) = 0.0 |
| 575 | 8902 | ALLOC_(G%mask2dCv(isd:ied,JsdB:JedB)) ; G%mask2dCv(:,:) = 0.0 |
| 576 | 8902 | ALLOC_(G%OBCmaskCv(isd:ied,JsdB:JedB)) ; G%OBCmaskCv(:,:) = 0.0 |
| 577 | 8971 | ALLOC_(G%mask2dBu(IsdB:IedB,JsdB:JedB)) ; G%mask2dBu(:,:) = 0.0 |
| 578 | 8773 | ALLOC_(G%geoLatT(isd:ied,jsd:jed)) ; G%geoLatT(:,:) = 0.0 |
| 579 | 8841 | ALLOC_(G%geoLatCu(IsdB:IedB,jsd:jed)) ; G%geoLatCu(:,:) = 0.0 |
| 580 | 8902 | ALLOC_(G%geoLatCv(isd:ied,JsdB:JedB)) ; G%geoLatCv(:,:) = 0.0 |
| 581 | 8971 | ALLOC_(G%geoLatBu(IsdB:IedB,JsdB:JedB)) ; G%geoLatBu(:,:) = 0.0 |
| 582 | 8773 | ALLOC_(G%geoLonT(isd:ied,jsd:jed)) ; G%geoLonT(:,:) = 0.0 |
| 583 | 8841 | ALLOC_(G%geoLonCu(IsdB:IedB,jsd:jed)) ; G%geoLonCu(:,:) = 0.0 |
| 584 | 8902 | ALLOC_(G%geoLonCv(isd:ied,JsdB:JedB)) ; G%geoLonCv(:,:) = 0.0 |
| 585 | 8971 | ALLOC_(G%geoLonBu(IsdB:IedB,JsdB:JedB)) ; G%geoLonBu(:,:) = 0.0 |
| 586 | ||
| 587 | 8902 | ALLOC_(G%dx_Cv(isd:ied,JsdB:JedB)) ; G%dx_Cv(:,:) = 0.0 |
| 588 | 8841 | ALLOC_(G%dy_Cu(IsdB:IedB,jsd:jed)) ; G%dy_Cu(:,:) = 0.0 |
| 589 | ||
| 590 | 8841 | ALLOC_(G%porous_DminU(IsdB:IedB,jsd:jed)) ; G%porous_DminU(:,:) = 0.0 |
| 591 | 8841 | ALLOC_(G%porous_DmaxU(IsdB:IedB,jsd:jed)) ; G%porous_DmaxU(:,:) = 0.0 |
| 592 | 8841 | ALLOC_(G%porous_DavgU(IsdB:IedB,jsd:jed)) ; G%porous_DavgU(:,:) = 0.0 |
| 593 | ||
| 594 | 8902 | ALLOC_(G%porous_DminV(isd:ied,JsdB:JedB)) ; G%porous_DminV(:,:) = 0.0 |
| 595 | 8902 | ALLOC_(G%porous_DmaxV(isd:ied,JsdB:JedB)) ; G%porous_DmaxV(:,:) = 0.0 |
| 596 | 8902 | ALLOC_(G%porous_DavgV(isd:ied,JsdB:JedB)) ; G%porous_DavgV(:,:) = 0.0 |
| 597 | ||
| 598 | 8841 | ALLOC_(G%areaCu(IsdB:IedB,jsd:jed)) ; G%areaCu(:,:) = 0.0 |
| 599 | 8902 | ALLOC_(G%areaCv(isd:ied,JsdB:JedB)) ; G%areaCv(:,:) = 0.0 |
| 600 | 8841 | ALLOC_(G%IareaCu(IsdB:IedB,jsd:jed)) ; G%IareaCu(:,:) = 0.0 |
| 601 | 8902 | ALLOC_(G%IareaCv(isd:ied,JsdB:JedB)) ; G%IareaCv(:,:) = 0.0 |
| 602 | ||
| 603 | 8773 | ALLOC_(G%bathyT(isd:ied, jsd:jed)) ; G%bathyT(:,:) = -G%Z_ref |
| 604 | 8773 | ALLOC_(G%meanSL(isd:ied, jsd:jed)) ; G%meanSL(:,:) = G%Z_ref |
| 605 | 8971 | ALLOC_(G%CoriolisBu(IsdB:IedB, JsdB:JedB)) ; G%CoriolisBu(:,:) = 0.0 |
| 606 | 8971 | ALLOC_(G%Coriolis2Bu(IsdB:IedB, JsdB:JedB)) ; G%Coriolis2Bu(:,:) = 0.0 |
| 607 | 8773 | ALLOC_(G%dF_dx(isd:ied, jsd:jed)) ; G%dF_dx(:,:) = 0.0 |
| 608 | 8773 | ALLOC_(G%dF_dy(isd:ied, jsd:jed)) ; G%dF_dy(:,:) = 0.0 |
| 609 | ||
| 610 | 8773 | ALLOC_(G%sin_rot(isd:ied,jsd:jed)) ; G%sin_rot(:,:) = 0.0 |
| 611 | 8773 | ALLOC_(G%cos_rot(isd:ied,jsd:jed)) ; G%cos_rot(:,:) = 1.0 |
| 612 | ||
| 613 | 121 | allocate(G%gridLonT(isg:ieg), source=0.0) |
| 614 | 122 | allocate(G%gridLonB(G%IsgB:G%IegB), source=0.0) |
| 615 | 61 | allocate(G%gridLatT(jsg:jeg), source=0.0) |
| 616 | 62 | allocate(G%gridLatB(G%JsgB:G%JegB), source=0.0) |
| 617 | 1 | end subroutine allocate_metrics |
| 618 | ||
| 619 | !> Release memory used by the ocean_grid_type and related structures. | |
| 620 | 1 | subroutine MOM_grid_end(G) |
| 621 | type(ocean_grid_type), intent(inout) :: G !< The horizontal grid type | |
| 622 | ||
| 623 | 1 | deallocate(G%Block) |
| 624 | ||
| 625 | 1 | if (G%bathymetry_at_vel) then |
| 626 | 0 | DEALLOC_(G%Dblock_u) ; DEALLOC_(G%Dopen_u) |
| 627 | 0 | DEALLOC_(G%Dblock_v) ; DEALLOC_(G%Dopen_v) |
| 628 | endif | |
| 629 | ||
| 630 | 1 | DEALLOC_(G%dxT) ; DEALLOC_(G%dxCu) ; DEALLOC_(G%dxCv) ; DEALLOC_(G%dxBu) |
| 631 | 1 | DEALLOC_(G%IdxT) ; DEALLOC_(G%IdxCu) ; DEALLOC_(G%IdxCv) ; DEALLOC_(G%IdxBu) |
| 632 | ||
| 633 | 1 | DEALLOC_(G%dyT) ; DEALLOC_(G%dyCu) ; DEALLOC_(G%dyCv) ; DEALLOC_(G%dyBu) |
| 634 | 1 | DEALLOC_(G%IdyT) ; DEALLOC_(G%IdyCu) ; DEALLOC_(G%IdyCv) ; DEALLOC_(G%IdyBu) |
| 635 | ||
| 636 | 1 | DEALLOC_(G%IdxCu_OBCmask) ; DEALLOC_(G%IdyCv_OBCmask) |
| 637 | ||
| 638 | 1 | DEALLOC_(G%areaT) ; DEALLOC_(G%IareaT) |
| 639 | 1 | DEALLOC_(G%areaBu) ; DEALLOC_(G%IareaBu) |
| 640 | 1 | DEALLOC_(G%areaCu) ; DEALLOC_(G%IareaCu) |
| 641 | 1 | DEALLOC_(G%areaCv) ; DEALLOC_(G%IareaCv) |
| 642 | ||
| 643 | 1 | DEALLOC_(G%mask2dT) ; DEALLOC_(G%mask2dCu) ; DEALLOC_(G%OBCmaskCu) |
| 644 | 1 | DEALLOC_(G%mask2dCv) ; DEALLOC_(G%OBCmaskCv) ; DEALLOC_(G%mask2dBu) |
| 645 | ||
| 646 | 1 | DEALLOC_(G%geoLatT) ; DEALLOC_(G%geoLatCu) |
| 647 | 1 | DEALLOC_(G%geoLatCv) ; DEALLOC_(G%geoLatBu) |
| 648 | 1 | DEALLOC_(G%geoLonT) ; DEALLOC_(G%geoLonCu) |
| 649 | 1 | DEALLOC_(G%geoLonCv) ; DEALLOC_(G%geoLonBu) |
| 650 | ||
| 651 | 1 | DEALLOC_(G%dx_Cv) ; DEALLOC_(G%dy_Cu) |
| 652 | ||
| 653 | 1 | DEALLOC_(G%bathyT) ; DEALLOC_(G%meanSL) |
| 654 | 1 | DEALLOC_(G%CoriolisBu) ; DEALLOC_(G%Coriolis2Bu) |
| 655 | 1 | DEALLOC_(G%dF_dx) ; DEALLOC_(G%dF_dy) |
| 656 | 1 | DEALLOC_(G%sin_rot) ; DEALLOC_(G%cos_rot) |
| 657 | ||
| 658 | 1 | DEALLOC_(G%porous_DminU) ; DEALLOC_(G%porous_DmaxU) ; DEALLOC_(G%porous_DavgU) |
| 659 | 1 | DEALLOC_(G%porous_DminV) ; DEALLOC_(G%porous_DmaxV) ; DEALLOC_(G%porous_DavgV) |
| 660 | ||
| 661 | 1 | deallocate(G%gridLonT) ; deallocate(G%gridLatT) |
| 662 | 1 | deallocate(G%gridLonB) ; deallocate(G%gridLatB) |
| 663 | ||
| 664 | ! The cursory flag avoids doing any deallocation of memory in the underlying | |
| 665 | ! infrastructure to avoid problems due to shared pointers. | |
| 666 | 1 | call deallocate_MOM_domain(G%Domain, cursory=.true.) |
| 667 | ||
| 668 | 1 | end subroutine MOM_grid_end |
| 669 | ||
| 670 | !> \namespace mom_grid | |
| 671 | !! | |
| 672 | !! Grid metrics and their inverses are labelled according to their staggered location on a Arakawa C (or B) grid. | |
| 673 | !! - Metrics centered on h- or T-points are labelled T, e.g. dxT is the distance across the cell in the x-direction. | |
| 674 | !! - Metrics centered on u-points are labelled Cu (C-grid u location). e.g. dyCu is the y-distance between | |
| 675 | !! two corners of a T-cell. | |
| 676 | !! - Metrics centered on v-points are labelled Cv (C-grid v location). e.g. dyCv is the y-distance between two -points. | |
| 677 | !! - Metrics centered on q-points are labelled Bu (B-grid u,v location). e.g. areaBu is the area centered on a q-point. | |
| 678 | !! | |
| 679 | !! \image html Grid_metrics.png "The labelling of distances (grid metrics) at various staggered | |
| 680 | !! location on an T-cell and around a q-point." | |
| 681 | !! | |
| 682 | !! Areas centered at T-, u-, v- and q- points are `areaT`, `areaCu`, `areaCv` and `areaBu` respectively. | |
| 683 | !! | |
| 684 | !! The reciprocal of metrics are pre-calculated and also stored in the ocean_grid_type with a I prepended to the name. | |
| 685 | !! For example, `1./areaT` is called `IareaT`, and `1./dyCv` is `IdyCv`. | |
| 686 | !! | |
| 687 | !! Geographic latitude and longitude (or model coordinates if not on a sphere) are stored in | |
| 688 | !! `geoLatT`, `geoLonT` for T-points. | |
| 689 | !! u-, v- and q- point coordinates are follow same pattern of replacing T with Cu, Cv and Bu respectively. | |
| 690 | !! | |
| 691 | !! Each location also has a 2D mask indicating whether the entire column is land or ocean. | |
| 692 | !! `mask2dT` is 1 if the column is wet or 0 if the T-cell is land. | |
| 693 | !! `mask2dCu` is 1 if both neighboring columns are ocean, and 0 if either is land. | |
| 694 | !! `OBCmasku` is 1 if both neighboring columns are ocean, and 0 if either is land of if this is OBC point. | |
| 695 | ||
| 696 | 0 | end module MOM_grid |