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 | !> Initializes horizontal grid | |
| 6 | module MOM_grid_initialize | |
| 7 | ||
| 8 | use MOM_checksums, only : hchksum, Bchksum, uvchksum, hchksum_pair, Bchksum_pair | |
| 9 | use MOM_domains, only : pass_var, pass_vector, pe_here, root_PE, broadcast | |
| 10 | use MOM_domains, only : AGRID, BGRID_NE, CGRID_NE, To_All, Scalar_Pair | |
| 11 | use MOM_domains, only : To_North, To_South, To_East, To_West | |
| 12 | use MOM_domains, only : MOM_domain_type, clone_MOM_domain, deallocate_MOM_domain | |
| 13 | use MOM_dyn_horgrid, only : dyn_horgrid_type, set_derived_dyn_horgrid | |
| 14 | use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, is_root_pe | |
| 15 | use MOM_error_handler, only : callTree_enter, callTree_leave | |
| 16 | use MOM_file_parser, only : get_param, log_param, log_version, param_file_type | |
| 17 | use MOM_io, only : MOM_read_data, slasher, file_exists, stdout | |
| 18 | use MOM_io, only : CORNER, NORTH_FACE, EAST_FACE | |
| 19 | use MOM_unit_scaling, only : unit_scale_type | |
| 20 | ||
| 21 | implicit none ; private | |
| 22 | ||
| 23 | public set_grid_metrics, initialize_masks, Adcroft_reciprocal | |
| 24 | ||
| 25 | ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional | |
| 26 | ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with | |
| 27 | ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units | |
| 28 | ! vary with the Boussinesq approximation, the Boussinesq variant is given first. | |
| 29 | ||
| 30 | !> Global positioning system (aka container for information to describe the grid) | |
| 31 | type, private :: GPS ; private | |
| 32 | real :: len_lon !< The longitudinal or x-direction length of the domain [degrees_E] or [km] or [m]. | |
| 33 | real :: len_lat !< The latitudinal or y-direction length of the domain [degrees_N] or [km] or [m]. | |
| 34 | real :: west_lon !< The western longitude of the domain or the equivalent | |
| 35 | !! starting value for the x-axis [degrees_E] or [km] or [m]. | |
| 36 | real :: south_lat !< The southern latitude of the domain or the equivalent | |
| 37 | !! starting value for the y-axis [degrees_N] or [km] or [m]. | |
| 38 | real :: Rad_Earth_L !< The radius of the Earth in rescaled units [L ~> m] | |
| 39 | real :: Lat_enhance_factor !< The amount by which the meridional resolution | |
| 40 | !! is enhanced within LAT_EQ_ENHANCE of the equator [nondim] | |
| 41 | real :: Lat_eq_enhance !< The latitude range to the north and south of the equator | |
| 42 | !! over which the resolution is enhanced [degrees_N] | |
| 43 | logical :: isotropic !< If true, an isotropic grid on a sphere (also known as a Mercator grid) | |
| 44 | !! is used. With an isotropic grid, the meridional extent of the domain | |
| 45 | !! (LENLAT), the zonal extent (LENLON), and the number of grid points in each | |
| 46 | !! direction are _not_ independent. In MOM the meridional extent is determined | |
| 47 | !! to fit the zonal extent and the number of grid points, while grid is | |
| 48 | !! perfectly isotropic. | |
| 49 | logical :: equator_reference !< If true, the grid is defined to have the equator at the | |
| 50 | !! nearest q- or h- grid point to (-LOWLAT*NJGLOBAL/LENLAT). | |
| 51 | integer :: niglobal !< The number of i-points in the global grid computational domain | |
| 52 | integer :: njglobal !< The number of j-points in the global grid computational domain | |
| 53 | end type GPS | |
| 54 | ||
| 55 | contains | |
| 56 | ||
| 57 | !> set_grid_metrics is used to set the primary values in the model's horizontal | |
| 58 | !! grid. The bathymetry, land-sea mask and any restricted channel widths are | |
| 59 | !! not known yet, so these are set later. | |
| 60 | 1 | subroutine set_grid_metrics(G, param_file, US) |
| 61 | type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type | |
| 62 | type(param_file_type), intent(in) :: param_file !< Parameter file structure | |
| 63 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 64 | ||
| 65 | ! Local variables | |
| 66 | ! This include declares and sets the variable "version". | |
| 67 | # include "version_variable.h" | |
| 68 | logical :: debug | |
| 69 | character(len=256) :: config | |
| 70 | ||
| 71 | 1 | call callTree_enter("set_grid_metrics(), MOM_grid_initialize.F90") |
| 72 | 1 | call log_version(param_file, "MOM_grid_init", version, "") |
| 73 | call get_param(param_file, "MOM_grid_init", "GRID_CONFIG", config, & | |
| 74 | "A character string that determines the method for "//& | |
| 75 | "defining the horizontal grid. Current options are: \n"//& | |
| 76 | " \t mosaic - read the grid from a mosaic (supergrid) \n"//& | |
| 77 | " \t file set by GRID_FILE.\n"//& | |
| 78 | " \t cartesian - use a (flat) Cartesian grid.\n"//& | |
| 79 | " \t spherical - use a simple spherical grid.\n"//& | |
| 80 | " \t mercator - use a Mercator spherical grid.", & | |
| 81 | 1 | fail_if_missing=.true.) |
| 82 | call get_param(param_file, "MOM_grid_init", "DEBUG", debug, & | |
| 83 | "If true, write out verbose debugging data.", & | |
| 84 | 1 | default=.false., debuggingParam=.true.) |
| 85 | ||
| 86 | ! These are defaults that may be changed in the next select block. | |
| 87 | 1 | G%x_axis_units = "degrees_east" ; G%y_axis_units = "degrees_north" |
| 88 | 1 | G%x_ax_unit_short = "degrees_E" ; G%y_ax_unit_short = "degrees_N" |
| 89 | 1 | G%grid_unit_to_L = 0.0 |
| 90 | 1 | G%Rad_Earth_L = -1.0*US%m_to_L ; G%len_lat = 0.0 ; G%len_lon = 0.0 |
| 91 | 2 | select case (trim(config)) |
| 92 | 0 | case ("mosaic"); call set_grid_metrics_from_mosaic(G, param_file, US) |
| 93 | 0 | case ("cartesian"); call set_grid_metrics_cartesian(G, param_file, US) |
| 94 | 1 | case ("spherical"); call set_grid_metrics_spherical(G, param_file, US) |
| 95 | 0 | case ("mercator"); call set_grid_metrics_mercator(G, param_file, US) |
| 96 | case ("file"); call MOM_error(FATAL, "MOM_grid_init: set_grid_metrics "//& | |
| 97 | 'GRID_CONFIG "file" is no longer a supported option. Use a '//& | |
| 98 | 0 | 'mosaic file ("mosaic") or one of the analytic forms instead.') |
| 99 | case default ; call MOM_error(FATAL, "MOM_grid_init: set_grid_metrics "//& | |
| 100 | 2 | "Unrecognized grid configuration "//trim(config)) |
| 101 | end select | |
| 102 | 1 | if (G%Rad_Earth_L <= 0.0) then |
| 103 | ! The grid metrics were set with an option that does not explicitly initialize Rad_Earth. | |
| 104 | call get_param(param_file, "MOM_grid_init", "RAD_EARTH", G%Rad_Earth_L, & | |
| 105 | 0 | "The radius of the Earth.", units="m", default=6.378e6, scale=US%m_to_L) |
| 106 | endif | |
| 107 | ||
| 108 | ! Calculate derived metrics (i.e. reciprocals and products) | |
| 109 | 1 | call callTree_enter("set_derived_metrics(), MOM_grid_initialize.F90") |
| 110 | 1 | call set_derived_dyn_horgrid(G, US) |
| 111 | 1 | call callTree_leave("set_derived_metrics()") |
| 112 | ||
| 113 | 1 | if (debug) call grid_metrics_chksum('MOM_grid_init/set_grid_metrics', G, US) |
| 114 | ||
| 115 | 1 | call callTree_leave("set_grid_metrics()") |
| 116 | 1 | end subroutine set_grid_metrics |
| 117 | ||
| 118 | ! ------------------------------------------------------------------------------ | |
| 119 | ||
| 120 | !> grid_metrics_chksum performs a set of checksums on metrics on the grid for | |
| 121 | !! debugging. | |
| 122 | 0 | subroutine grid_metrics_chksum(parent, G, US) |
| 123 | character(len=*), intent(in) :: parent !< A string identifying the caller | |
| 124 | type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type | |
| 125 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 126 | ||
| 127 | integer :: halo | |
| 128 | ||
| 129 | 0 | halo = min(G%ied-G%iec, G%jed-G%jec, 1) |
| 130 | ||
| 131 | call hchksum_pair(trim(parent)//': d[xy]T', G%dxT, G%dyT, G%HI, & | |
| 132 | 0 | haloshift=halo, unscale=US%L_to_m, scalar_pair=.true.) |
| 133 | ||
| 134 | 0 | call uvchksum(trim(parent)//': dxC[uv]', G%dxCu, G%dyCv, G%HI, haloshift=halo, unscale=US%L_to_m) |
| 135 | ||
| 136 | 0 | call uvchksum(trim(parent)//': dxC[uv]', G%dyCu, G%dxCv, G%HI, haloshift=halo, unscale=US%L_to_m) |
| 137 | ||
| 138 | 0 | call Bchksum_pair(trim(parent)//': dxB[uv]', G%dxBu, G%dyBu, G%HI, haloshift=halo, unscale=US%L_to_m) |
| 139 | ||
| 140 | call hchksum_pair(trim(parent)//': Id[xy]T', G%IdxT, G%IdyT, G%HI, & | |
| 141 | 0 | haloshift=halo, unscale=US%m_to_L, scalar_pair=.true.) |
| 142 | ||
| 143 | 0 | call uvchksum(trim(parent)//': Id[xy]C[uv]', G%IdxCu, G%IdyCv, G%HI, haloshift=halo, unscale=US%m_to_L) |
| 144 | ||
| 145 | 0 | call uvchksum(trim(parent)//': Id[xy]C[uv]', G%IdyCu, G%IdxCv, G%HI, haloshift=halo, unscale=US%m_to_L) |
| 146 | ||
| 147 | 0 | call Bchksum_pair(trim(parent)//': Id[xy]B[uv]', G%IdxBu, G%IdyBu, G%HI, haloshift=halo, unscale=US%m_to_L) |
| 148 | ||
| 149 | 0 | call hchksum(G%areaT, trim(parent)//': areaT',G%HI, haloshift=halo, unscale=US%L_to_m**2) |
| 150 | 0 | call Bchksum(G%areaBu, trim(parent)//': areaBu',G%HI, haloshift=halo, unscale=US%L_to_m**2) |
| 151 | ||
| 152 | 0 | call hchksum(G%IareaT, trim(parent)//': IareaT',G%HI, haloshift=halo, unscale=US%m_to_L**2) |
| 153 | 0 | call Bchksum(G%IareaBu, trim(parent)//': IareaBu',G%HI, haloshift=halo, unscale=US%m_to_L**2) |
| 154 | ||
| 155 | 0 | call hchksum(G%geoLonT,trim(parent)//': geoLonT',G%HI, haloshift=halo) |
| 156 | 0 | call hchksum(G%geoLatT,trim(parent)//': geoLatT',G%HI, haloshift=halo) |
| 157 | ||
| 158 | 0 | call Bchksum(G%geoLonBu, trim(parent)//': geoLonBu',G%HI, haloshift=halo) |
| 159 | 0 | call Bchksum(G%geoLatBu, trim(parent)//': geoLatBu',G%HI, haloshift=halo) |
| 160 | ||
| 161 | 0 | call uvchksum(trim(parent)//': geoLonC[uv]', G%geoLonCu, G%geoLonCv, G%HI, haloshift=halo) |
| 162 | ||
| 163 | 0 | call uvchksum(trim(parent)//': geoLatC[uv]', G%geoLatCu, G%geoLatCv, G%HI, haloshift=halo) |
| 164 | ||
| 165 | 0 | end subroutine grid_metrics_chksum |
| 166 | ||
| 167 | ! ------------------------------------------------------------------------------ | |
| 168 | ||
| 169 | !> Sets the grid metrics from a mosaic file. | |
| 170 | 0 | subroutine set_grid_metrics_from_mosaic(G, param_file, US) |
| 171 | type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type | |
| 172 | type(param_file_type), intent(in) :: param_file !< Parameter file structure | |
| 173 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 174 | ||
| 175 | ! Local variables | |
| 176 | ! These are symmetric arrays, corresponding to the data in the mosaic file | |
| 177 | 0 | real, dimension(2*G%isd-2:2*G%ied+1,2*G%jsd-2:2*G%jed+1) :: tmpT ! Areas [L2 ~> m2] |
| 178 | 0 | real, dimension(2*G%isd-3:2*G%ied+1,2*G%jsd-2:2*G%jed+1) :: tmpU ! East face supergrid spacing [L ~> m] |
| 179 | 0 | real, dimension(2*G%isd-2:2*G%ied+1,2*G%jsd-3:2*G%jed+1) :: tmpV ! North face supergrid spacing [L ~> m] |
| 180 | 0 | real, dimension(2*G%isd-3:2*G%ied+1,2*G%jsd-3:2*G%jed+1) :: tmpZ ! Corner latitudes [degrees_N] or |
| 181 | ! longitudes [degrees_E] | |
| 182 | 0 | real, dimension(:,:), allocatable :: tmpGlbl ! A global array of axis labels [degrees_N] or [km] or [m] |
| 183 | character(len=200) :: filename, grid_file, inputdir | |
| 184 | character(len=64) :: mdl = "MOM_grid_init set_grid_metrics_from_mosaic" | |
| 185 | type(MOM_domain_type), pointer :: SGdom => NULL() ! Supergrid domain | |
| 186 | logical :: lon_bug ! If true use an older buggy answer in the tripolar longitude. | |
| 187 | integer :: i, j, i2, j2, ni, nj | |
| 188 | integer :: start(4), nread(4) | |
| 189 | ||
| 190 | 0 | call callTree_enter("set_grid_metrics_from_mosaic(), MOM_grid_initialize.F90") |
| 191 | ||
| 192 | call get_param(param_file, mdl, "GRID_FILE", grid_file, & | |
| 193 | "Name of the file from which to read horizontal grid data.", & | |
| 194 | 0 | fail_if_missing=.true.) |
| 195 | call get_param(param_file, mdl, "USE_TRIPOLAR_GEOLONB_BUG", lon_bug, & | |
| 196 | "If true, use older code that incorrectly sets the longitude "//& | |
| 197 | "in some points along the tripolar fold to be off by 360 degrees.", & | |
| 198 | 0 | default=.false.) |
| 199 | 0 | call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".") |
| 200 | 0 | inputdir = slasher(inputdir) |
| 201 | 0 | filename = trim(adjustl(inputdir)) // trim(adjustl(grid_file)) |
| 202 | 0 | call log_param(param_file, mdl, "INPUTDIR/GRID_FILE", filename) |
| 203 | 0 | if (.not.file_exists(filename)) & |
| 204 | call MOM_error(FATAL," set_grid_metrics_from_mosaic: Unable to open "//& | |
| 205 | 0 | trim(filename)) |
| 206 | ||
| 207 | !<MISSING CODE TO READ REFINEMENT LEVEL> | |
| 208 | ||
| 209 | call clone_MOM_domain(G%domain, SGdom, symmetric=.true., domain_name="MOM_MOSAIC", & | |
| 210 | 0 | refine=2, extra_halo=1) |
| 211 | ||
| 212 | ! Read X from the supergrid | |
| 213 | 0 | tmpZ(:,:) = 999. |
| 214 | 0 | call MOM_read_data(filename, 'x', tmpZ, SGdom, position=CORNER) |
| 215 | ||
| 216 | 0 | if (lon_bug) then |
| 217 | 0 | call pass_var(tmpZ, SGdom, position=CORNER) |
| 218 | else | |
| 219 | 0 | call pass_var(tmpZ, SGdom, position=CORNER, inner_halo=0) |
| 220 | endif | |
| 221 | 0 | call extrapolate_metric(tmpZ, 2*(G%jsc-G%jsd)+2, missing=999.) |
| 222 | 0 | do j=G%jsd,G%jed ; do i=G%isd,G%ied ; i2 = 2*i ; j2 = 2*j |
| 223 | 0 | G%geoLonT(i,j) = tmpZ(i2-1,j2-1) |
| 224 | enddo ; enddo | |
| 225 | 0 | do J=G%JsdB,G%JedB ; do I=G%IsdB,G%IedB ; i2 = 2*I ; j2 = 2*J |
| 226 | 0 | G%geoLonBu(I,J) = tmpZ(i2,j2) |
| 227 | enddo ; enddo | |
| 228 | 0 | do j=G%jsd,G%jed ; do I=G%IsdB,G%IedB ; i2 = 2*i ; j2 = 2*j |
| 229 | 0 | G%geoLonCu(I,j) = tmpZ(i2,j2-1) |
| 230 | enddo ; enddo | |
| 231 | 0 | do J=G%JsdB,G%JedB ; do i=G%isd,G%ied ; i2 = 2*i ; j2 = 2*J |
| 232 | 0 | G%geoLonCv(i,J) = tmpZ(i2-1,j2) |
| 233 | enddo ; enddo | |
| 234 | ! For some reason, this messes up the solution... | |
| 235 | ! call pass_var(G%geoLonBu, G%domain, position=CORNER) | |
| 236 | ||
| 237 | ! Read Y from the supergrid | |
| 238 | 0 | tmpZ(:,:) = 999. |
| 239 | 0 | call MOM_read_data(filename, 'y', tmpZ, SGdom, position=CORNER) |
| 240 | ||
| 241 | 0 | call pass_var(tmpZ, SGdom, position=CORNER) |
| 242 | 0 | call extrapolate_metric(tmpZ, 2*(G%jsc-G%jsd)+2, missing=999.) |
| 243 | 0 | do j=G%jsd,G%jed ; do i=G%isd,G%ied ; i2 = 2*i ; j2 = 2*j |
| 244 | 0 | G%geoLatT(i,j) = tmpZ(i2-1,j2-1) |
| 245 | enddo ; enddo | |
| 246 | 0 | do J=G%JsdB,G%JedB ; do I=G%IsdB,G%IedB ; i2 = 2*I ; j2 = 2*J |
| 247 | 0 | G%geoLatBu(I,J) = tmpZ(i2,j2) |
| 248 | enddo ; enddo | |
| 249 | 0 | do j=G%jsd,G%jed ; do I=G%IsdB,G%IedB ; i2 = 2*i ; j2 = 2*j |
| 250 | 0 | G%geoLatCu(I,j) = tmpZ(i2,j2-1) |
| 251 | enddo ; enddo | |
| 252 | 0 | do J=G%JsdB,G%JedB ; do i=G%isd,G%ied ; i2 = 2*i ; j2 = 2*J |
| 253 | 0 | G%geoLatCv(i,J) = tmpZ(i2-1,j2) |
| 254 | enddo ; enddo | |
| 255 | ||
| 256 | ! This routine could be modified to support the use of a mosaic using Cartesian grid coordinates, | |
| 257 | ! in which case the values of G%x_axis_units, G%y_axis_units and G%grid_unit_to_L would need to be | |
| 258 | ! reset appropriately here, but this option has not yet been implemented, and the grid coordinates | |
| 259 | ! are assumed to be degrees of longitude and latitude. | |
| 260 | ||
| 261 | ! Read DX,DY from the supergrid | |
| 262 | 0 | tmpU(:,:) = 0. ; tmpV(:,:) = 0. |
| 263 | 0 | call MOM_read_data(filename, 'dx', tmpV, SGdom, position=NORTH_FACE, scale=US%m_to_L) |
| 264 | 0 | call MOM_read_data(filename, 'dy', tmpU, SGdom, position=EAST_FACE, scale=US%m_to_L) |
| 265 | 0 | call pass_vector(tmpU, tmpV, SGdom, To_All+Scalar_Pair, CGRID_NE) |
| 266 | 0 | call extrapolate_metric(tmpV, 2*(G%jsc-G%jsd)+2, missing=0.) |
| 267 | 0 | call extrapolate_metric(tmpU, 2*(G%jsc-G%jsd)+2, missing=0.) |
| 268 | ||
| 269 | 0 | do j=G%jsd,G%jed ; do i=G%isd,G%ied ; i2 = 2*i ; j2 = 2*j |
| 270 | 0 | G%dxT(i,j) = tmpV(i2-1,j2-1) + tmpV(i2,j2-1) |
| 271 | 0 | G%dyT(i,j) = tmpU(i2-1,j2-1) + tmpU(i2-1,j2) |
| 272 | enddo ; enddo | |
| 273 | ||
| 274 | 0 | do j=G%jsd,G%jed ; do I=G%IsdB,G%IedB ; i2 = 2*i ; j2 = 2*j |
| 275 | 0 | G%dxCu(I,j) = tmpV(i2,j2-1) + tmpV(i2+1,j2-1) |
| 276 | 0 | G%dyCu(I,j) = tmpU(i2,j2-1) + tmpU(i2,j2) |
| 277 | enddo ; enddo | |
| 278 | ||
| 279 | 0 | do J=G%JsdB,G%JedB ; do i=G%isd,G%ied ; i2 = 2*i ; j2 = 2*j |
| 280 | 0 | G%dxCv(i,J) = tmpV(i2-1,j2) + tmpV(i2,j2) |
| 281 | 0 | G%dyCv(i,J) = tmpU(i2-1,j2) + tmpU(i2-1,j2+1) |
| 282 | enddo ; enddo | |
| 283 | ||
| 284 | 0 | do J=G%JsdB,G%JedB ; do I=G%IsdB,G%IedB ; i2 = 2*i ; j2 = 2*j |
| 285 | 0 | G%dxBu(I,J) = tmpV(i2,j2) + tmpV(i2+1,j2) |
| 286 | 0 | G%dyBu(I,J) = tmpU(i2,j2) + tmpU(i2,j2+1) |
| 287 | enddo ; enddo | |
| 288 | ||
| 289 | ! Read AREA from the supergrid | |
| 290 | 0 | tmpT(:,:) = 0. |
| 291 | 0 | call MOM_read_data(filename, 'area', tmpT, SGdom, scale=US%m_to_L**2) |
| 292 | 0 | call pass_var(tmpT, SGdom) |
| 293 | 0 | call extrapolate_metric(tmpT, 2*(G%jsc-G%jsd)+2, missing=0.) |
| 294 | ||
| 295 | 0 | do j=G%jsd,G%jed ; do i=G%isd,G%ied ; i2 = 2*i ; j2 = 2*j |
| 296 | G%areaT(i,j) = (tmpT(i2-1,j2-1) + tmpT(i2,j2)) + & | |
| 297 | 0 | (tmpT(i2-1,j2) + tmpT(i2,j2-1)) |
| 298 | enddo ; enddo | |
| 299 | 0 | do J=G%JsdB,G%JedB ; do I=G%IsdB,G%IedB ; i2 = 2*i ; j2 = 2*j |
| 300 | G%areaBu(I,J) = (tmpT(i2,j2) + tmpT(i2+1,j2+1)) + & | |
| 301 | 0 | (tmpT(i2,j2+1) + tmpT(i2+1,j2)) |
| 302 | enddo ; enddo | |
| 303 | ||
| 304 | 0 | ni = SGdom%niglobal |
| 305 | 0 | nj = SGdom%njglobal |
| 306 | 0 | call deallocate_MOM_domain(SGdom) |
| 307 | ||
| 308 | 0 | call pass_vector(G%dyCu, G%dxCv, G%Domain, To_All+Scalar_Pair, CGRID_NE) |
| 309 | 0 | call pass_vector(G%dxCu, G%dyCv, G%Domain, To_All+Scalar_Pair, CGRID_NE) |
| 310 | 0 | call pass_vector(G%dxBu, G%dyBu, G%Domain, To_All+Scalar_Pair, BGRID_NE) |
| 311 | 0 | call pass_var(G%areaT, G%Domain) |
| 312 | 0 | call pass_var(G%areaBu, G%Domain, position=CORNER) |
| 313 | ||
| 314 | ! Construct axes for diagnostic output (only necessary because "ferret" uses | |
| 315 | ! broken convention for interpretting netCDF files). | |
| 316 | 0 | start(:) = 1 ; nread(:) = 1 |
| 317 | 0 | start(2) = 2 ; nread(1) = ni+1 ; nread(2) = 2 |
| 318 | 0 | allocate( tmpGlbl(ni+1,2) ) |
| 319 | 0 | if (is_root_PE()) & |
| 320 | call MOM_read_data(filename, "x", tmpGlbl, start, nread, & | |
| 321 | 0 | no_domain=.TRUE., turns=G%HI%turns) |
| 322 | 0 | call broadcast(tmpGlbl, 2*(ni+1), root_PE()) |
| 323 | ||
| 324 | ! I don't know why the second axis is 1 or 2 here. -RWH | |
| 325 | 0 | do i=G%isg,G%ieg |
| 326 | 0 | G%gridLonT(i) = tmpGlbl(2*(i-G%isg)+2,2) |
| 327 | enddo | |
| 328 | ! Note that the dynamic grid always uses symmetric memory for the global | |
| 329 | ! arrays G%gridLatB and G%gridLonB. | |
| 330 | 0 | do I=G%isg-1,G%ieg |
| 331 | 0 | G%gridLonB(I) = tmpGlbl(2*(I-G%isg)+3,1) |
| 332 | enddo | |
| 333 | 0 | deallocate( tmpGlbl ) |
| 334 | ||
| 335 | 0 | allocate( tmpGlbl(1, nj+1) ) |
| 336 | 0 | start(:) = 1 ; nread(:) = 1 |
| 337 | 0 | start(1) = int(ni/4)+1 ; nread(2) = nj+1 |
| 338 | 0 | if (is_root_PE()) & |
| 339 | call MOM_read_data(filename, "y", tmpGlbl, start, nread, & | |
| 340 | 0 | no_domain=.TRUE., turns=G%HI%turns) |
| 341 | 0 | call broadcast(tmpGlbl, nj+1, root_PE()) |
| 342 | ||
| 343 | 0 | do j=G%jsg,G%jeg |
| 344 | 0 | G%gridLatT(j) = tmpGlbl(1,2*(j-G%jsg)+2) |
| 345 | enddo | |
| 346 | 0 | do J=G%jsg-1,G%jeg |
| 347 | 0 | G%gridLatB(J) = tmpGlbl(1,2*(j-G%jsg)+3) |
| 348 | enddo | |
| 349 | 0 | deallocate( tmpGlbl ) |
| 350 | ||
| 351 | 0 | call callTree_leave("set_grid_metrics_from_mosaic()") |
| 352 | 0 | end subroutine set_grid_metrics_from_mosaic |
| 353 | ||
| 354 | ||
| 355 | ! ------------------------------------------------------------------------------ | |
| 356 | ||
| 357 | !> Calculate the values of the metric terms for a Cartesian grid that | |
| 358 | !! might be used and save them in arrays. | |
| 359 | !! | |
| 360 | !! Within this subroutine, the x- and y- grid spacings and their | |
| 361 | !! inverses and the cell areas centered on h, q, u, and v points are | |
| 362 | !! calculated, as are the geographic locations of each of these 4 | |
| 363 | !! sets of points. | |
| 364 | 0 | subroutine set_grid_metrics_cartesian(G, param_file, US) |
| 365 | type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type | |
| 366 | type(param_file_type), intent(in) :: param_file !< Parameter file structure | |
| 367 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 368 | ! Local variables | |
| 369 | integer :: i, j, isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, I1off, J1off | |
| 370 | integer :: niglobal, njglobal | |
| 371 | 0 | real :: grid_latT(G%jsd:G%jed), grid_latB(G%JsdB:G%JedB) ! Axis labels [degrees_N] or [km] or [m] |
| 372 | 0 | real :: grid_lonT(G%isd:G%ied), grid_lonB(G%IsdB:G%IedB) ! Axis labels [degrees_E] or [km] or [m] |
| 373 | real :: dx_everywhere, dy_everywhere ! Grid spacings [L ~> m]. | |
| 374 | real :: I_dx, I_dy ! Inverse grid spacings [L-1 ~> m-1]. | |
| 375 | real :: PI ! The ratio of the circumference of a circle to its diameter [nondim] | |
| 376 | character(len=80) :: units_temp | |
| 377 | character(len=48) :: mdl = "MOM_grid_init set_grid_metrics_cartesian" | |
| 378 | ||
| 379 | 0 | niglobal = G%Domain%niglobal ; njglobal = G%Domain%njglobal |
| 380 | 0 | isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed |
| 381 | 0 | IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB |
| 382 | 0 | I1off = G%idg_offset ; J1off = G%jdg_offset |
| 383 | ||
| 384 | 0 | call callTree_enter("set_grid_metrics_cartesian(), MOM_grid_initialize.F90") |
| 385 | ||
| 386 | 0 | PI = 4.0*atan(1.0) |
| 387 | ||
| 388 | call get_param(param_file, mdl, "AXIS_UNITS", units_temp, & | |
| 389 | "The units for the Cartesian axes. Valid entries are: \n"//& | |
| 390 | " \t degrees - degrees of latitude and longitude \n"//& | |
| 391 | " \t m or meter(s) - meters \n"//& | |
| 392 | 0 | " \t k or km or kilometer(s) - kilometers", default="degrees") |
| 393 | 0 | if (trim(units_temp) == "k") units_temp = "km" |
| 394 | ||
| 395 | call get_param(param_file, mdl, "SOUTHLAT", G%south_lat, & | |
| 396 | "The southern latitude of the domain or the equivalent "//& | |
| 397 | "starting value for the y-axis.", units=units_temp, & | |
| 398 | 0 | fail_if_missing=.true.) |
| 399 | call get_param(param_file, mdl, "LENLAT", G%len_lat, & | |
| 400 | "The latitudinal or y-direction length of the domain.", & | |
| 401 | 0 | units=units_temp, fail_if_missing=.true.) |
| 402 | call get_param(param_file, mdl, "WESTLON", G%west_lon, & | |
| 403 | "The western longitude of the domain or the equivalent "//& | |
| 404 | "starting value for the x-axis.", units=units_temp, & | |
| 405 | 0 | default=0.0) |
| 406 | call get_param(param_file, mdl, "LENLON", G%len_lon, & | |
| 407 | "The longitudinal or x-direction length of the domain.", & | |
| 408 | 0 | units=units_temp, fail_if_missing=.true.) |
| 409 | call get_param(param_file, mdl, "RAD_EARTH", G%Rad_Earth_L, & | |
| 410 | 0 | "The radius of the Earth.", units="m", default=6.378e6, scale=US%m_to_L) |
| 411 | ||
| 412 | 0 | if (units_temp(1:1) == 'k') then |
| 413 | 0 | G%x_axis_units = "kilometers" ; G%y_axis_units = "kilometers" |
| 414 | 0 | G%x_ax_unit_short = "km" ; G%y_ax_unit_short = "km" |
| 415 | 0 | elseif (units_temp(1:1) == 'm') then |
| 416 | 0 | G%x_axis_units = "meters" ; G%y_axis_units = "meters" |
| 417 | 0 | G%x_ax_unit_short = "m" ; G%y_ax_unit_short = "m" |
| 418 | endif | |
| 419 | 0 | call log_param(param_file, mdl, "explicit AXIS_UNITS", G%x_axis_units) |
| 420 | ||
| 421 | ! Note that the dynamic grid always uses symmetric memory for the global | |
| 422 | ! arrays G%gridLatB and G%gridLonB. | |
| 423 | 0 | do J=G%jsg-1,G%jeg |
| 424 | 0 | G%gridLatB(j) = G%south_lat + G%len_lat*REAL(J-(G%jsg-1))/REAL(njglobal) |
| 425 | enddo | |
| 426 | 0 | do j=G%jsg,G%jeg |
| 427 | 0 | G%gridLatT(j) = G%south_lat + G%len_lat*(REAL(j-G%jsg)+0.5)/REAL(njglobal) |
| 428 | enddo | |
| 429 | 0 | do I=G%isg-1,G%ieg |
| 430 | 0 | G%gridLonB(i) = G%west_lon + G%len_lon*REAL(I-(G%isg-1))/REAL(niglobal) |
| 431 | enddo | |
| 432 | 0 | do i=G%isg,G%ieg |
| 433 | 0 | G%gridLonT(i) = G%west_lon + G%len_lon*(REAL(i-G%isg)+0.5)/REAL(niglobal) |
| 434 | enddo | |
| 435 | ||
| 436 | 0 | do J=JsdB,JedB |
| 437 | 0 | grid_latB(J) = G%south_lat + G%len_lat*REAL(J+J1off-(G%jsg-1))/REAL(njglobal) |
| 438 | enddo | |
| 439 | 0 | do j=jsd,jed |
| 440 | 0 | grid_latT(J) = G%south_lat + G%len_lat*(REAL(j+J1off-G%jsg)+0.5)/REAL(njglobal) |
| 441 | enddo | |
| 442 | 0 | do I=IsdB,IedB |
| 443 | 0 | grid_lonB(I) = G%west_lon + G%len_lon*REAL(i+I1off-(G%isg-1))/REAL(niglobal) |
| 444 | enddo | |
| 445 | 0 | do i=isd,ied |
| 446 | 0 | grid_lonT(i) = G%west_lon + G%len_lon*(REAL(i+I1off-G%isg)+0.5)/REAL(niglobal) |
| 447 | enddo | |
| 448 | ||
| 449 | 0 | if (units_temp(1:1) == 'k') then ! Axes are measured in km. |
| 450 | 0 | G%grid_unit_to_L = 1000.0*US%m_to_L |
| 451 | 0 | dx_everywhere = 1000.0*US%m_to_L * G%len_lon / (REAL(niglobal)) |
| 452 | 0 | dy_everywhere = 1000.0*US%m_to_L * G%len_lat / (REAL(njglobal)) |
| 453 | 0 | elseif (units_temp(1:1) == 'm') then ! Axes are measured in m. |
| 454 | 0 | G%grid_unit_to_L = US%m_to_L |
| 455 | 0 | dx_everywhere = US%m_to_L*G%len_lon / (REAL(niglobal)) |
| 456 | 0 | dy_everywhere = US%m_to_L*G%len_lat / (REAL(njglobal)) |
| 457 | else ! Axes are measured in degrees of latitude and longitude. | |
| 458 | 0 | dx_everywhere = G%Rad_Earth_L * G%len_lon * PI / (180.0 * niglobal) |
| 459 | 0 | dy_everywhere = G%Rad_Earth_L * G%len_lat * PI / (180.0 * njglobal) |
| 460 | endif | |
| 461 | ||
| 462 | 0 | I_dx = 1.0 / dx_everywhere ; I_dy = 1.0 / dy_everywhere |
| 463 | ||
| 464 | 0 | do J=JsdB,JedB ; do I=IsdB,IedB |
| 465 | 0 | G%geoLonBu(I,J) = grid_lonB(I) ; G%geoLatBu(I,J) = grid_latB(J) |
| 466 | ||
| 467 | 0 | G%dxBu(I,J) = dx_everywhere ; G%IdxBu(I,J) = I_dx |
| 468 | 0 | G%dyBu(I,J) = dy_everywhere ; G%IdyBu(I,J) = I_dy |
| 469 | 0 | G%areaBu(I,J) = dx_everywhere * dy_everywhere ; G%IareaBu(I,J) = I_dx * I_dy |
| 470 | enddo ; enddo | |
| 471 | ||
| 472 | 0 | do j=jsd,jed ; do i=isd,ied |
| 473 | 0 | G%geoLonT(i,j) = grid_lonT(i) ; G%geoLatT(i,j) = grid_LatT(j) |
| 474 | 0 | G%dxT(i,j) = dx_everywhere ; G%IdxT(i,j) = I_dx |
| 475 | 0 | G%dyT(i,j) = dy_everywhere ; G%IdyT(i,j) = I_dy |
| 476 | 0 | G%areaT(i,j) = dx_everywhere * dy_everywhere ; G%IareaT(i,j) = I_dx * I_dy |
| 477 | enddo ; enddo | |
| 478 | ||
| 479 | 0 | do j=jsd,jed ; do I=IsdB,IedB |
| 480 | 0 | G%geoLonCu(I,j) = grid_lonB(I) ; G%geoLatCu(I,j) = grid_LatT(j) |
| 481 | ||
| 482 | 0 | G%dxCu(I,j) = dx_everywhere ; G%IdxCu(I,j) = I_dx |
| 483 | 0 | G%dyCu(I,j) = dy_everywhere ; G%IdyCu(I,j) = I_dy |
| 484 | enddo ; enddo | |
| 485 | ||
| 486 | 0 | do J=JsdB,JedB ; do i=isd,ied |
| 487 | 0 | G%geoLonCv(i,J) = grid_lonT(i) ; G%geoLatCv(i,J) = grid_latB(J) |
| 488 | ||
| 489 | 0 | G%dxCv(i,J) = dx_everywhere ; G%IdxCv(i,J) = I_dx |
| 490 | 0 | G%dyCv(i,J) = dy_everywhere ; G%IdyCv(i,J) = I_dy |
| 491 | enddo ; enddo | |
| 492 | ||
| 493 | 0 | call callTree_leave("set_grid_metrics_cartesian()") |
| 494 | 0 | end subroutine set_grid_metrics_cartesian |
| 495 | ||
| 496 | ! ------------------------------------------------------------------------------ | |
| 497 | ||
| 498 | !> Calculate the values of the metric terms that might be used | |
| 499 | !! and save them in arrays. | |
| 500 | !! | |
| 501 | !! Within this subroutine, the x- and y- grid spacings and their | |
| 502 | !! inverses and the cell areas centered on h, q, u, and v points are | |
| 503 | !! calculated, as are the geographic locations of each of these 4 | |
| 504 | !! sets of points. | |
| 505 | 1 | subroutine set_grid_metrics_spherical(G, param_file, US) |
| 506 | type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type | |
| 507 | type(param_file_type), intent(in) :: param_file !< Parameter file structure | |
| 508 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 509 | ! Local variables | |
| 510 | real :: PI ! PI = 3.1415926... as 4*atan(1) [nondim] | |
| 511 | real :: PI_180 ! The conversion factor from degrees to radians [radians degree-1] | |
| 512 | integer :: i, j, isd, ied, jsd, jed | |
| 513 | integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, IsdB, IedB, JsdB, JedB | |
| 514 | integer :: i_offset, j_offset | |
| 515 | 2 | real :: grid_latT(G%jsd:G%jed), grid_latB(G%JsdB:G%JedB) ! Axis labels [degrees_N] |
| 516 | 2 | real :: grid_lonT(G%isd:G%ied), grid_lonB(G%IsdB:G%IedB) ! Axis labels [degrees_E] |
| 517 | real :: dLon ! The change in longitude between successive grid points [degrees_E] | |
| 518 | real :: dLat ! The change in latitude between successive grid points [degrees_N] | |
| 519 | real :: dL_di ! dLon rescaled from degrees to radians [radians] | |
| 520 | real :: latitude ! The latitude of a grid point [degrees_N] | |
| 521 | character(len=48) :: mdl = "MOM_grid_init set_grid_metrics_spherical" | |
| 522 | ||
| 523 | 1 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec |
| 524 | 1 | isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed |
| 525 | 1 | Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB |
| 526 | 1 | IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB |
| 527 | 1 | i_offset = G%idg_offset ; j_offset = G%jdg_offset |
| 528 | ||
| 529 | 1 | call callTree_enter("set_grid_metrics_spherical(), MOM_grid_initialize.F90") |
| 530 | ||
| 531 | ! Calculate the values of the metric terms that might be used | |
| 532 | ! and save them in arrays. | |
| 533 | 1 | PI = 4.0*atan(1.0) ; PI_180 = atan(1.0)/45. |
| 534 | ||
| 535 | call get_param(param_file, mdl, "SOUTHLAT", G%south_lat, & | |
| 536 | "The southern latitude of the domain.", units="degrees_N", & | |
| 537 | 1 | fail_if_missing=.true.) |
| 538 | call get_param(param_file, mdl, "LENLAT", G%len_lat, & | |
| 539 | "The latitudinal length of the domain.", units="degrees_N", & | |
| 540 | 1 | fail_if_missing=.true.) |
| 541 | call get_param(param_file, mdl, "WESTLON", G%west_lon, & | |
| 542 | "The western longitude of the domain.", units="degrees_E", & | |
| 543 | 1 | default=0.0) |
| 544 | call get_param(param_file, mdl, "LENLON", G%len_lon, & | |
| 545 | "The longitudinal length of the domain.", units="degrees_E", & | |
| 546 | 1 | fail_if_missing=.true.) |
| 547 | call get_param(param_file, mdl, "RAD_EARTH", G%Rad_Earth_L, & | |
| 548 | 1 | "The radius of the Earth.", units="m", default=6.378e6, scale=US%m_to_L) |
| 549 | ||
| 550 | 1 | dLon = G%len_lon/G%Domain%niglobal |
| 551 | 1 | dLat = G%len_lat/G%Domain%njglobal |
| 552 | ||
| 553 | ! Note that the dynamic grid always uses symmetric memory for the global | |
| 554 | ! arrays G%gridLatB and G%gridLonB. | |
| 555 | 62 | do j=G%jsg-1,G%jeg |
| 556 | 61 | latitude = G%south_lat + dLat*(REAL(J-(G%jsg-1))) |
| 557 | 62 | G%gridLatB(J) = MIN(MAX(latitude,-90.),90.) |
| 558 | enddo | |
| 559 | 61 | do j=G%jsg,G%jeg |
| 560 | 60 | latitude = G%south_lat + dLat*(REAL(j-G%jsg)+0.5) |
| 561 | 61 | G%gridLatT(j) = MIN(MAX(latitude,-90.),90.) |
| 562 | enddo | |
| 563 | 122 | do i=G%isg-1,G%ieg |
| 564 | 122 | G%gridLonB(I) = G%west_lon + dLon*(REAL(I-(G%isg-1))) |
| 565 | enddo | |
| 566 | 121 | do i=G%isg,G%ieg |
| 567 | 121 | G%gridLonT(i) = G%west_lon + dLon*(REAL(i-G%isg)+0.5) |
| 568 | enddo | |
| 569 | ||
| 570 | 70 | do J=JsdB,JedB |
| 571 | 69 | latitude = G%south_lat + dLat* REAL(J+J_offset-(G%jsg-1)) |
| 572 | 70 | grid_LatB(J) = MIN(MAX(latitude,-90.),90.) |
| 573 | enddo | |
| 574 | 69 | do j=jsd,jed |
| 575 | 68 | latitude = G%south_lat + dLat*(REAL(j+J_offset-G%jsg)+0.5) |
| 576 | 69 | grid_LatT(j) = MIN(MAX(latitude,-90.),90.) |
| 577 | enddo | |
| 578 | 130 | do I=IsdB,IedB |
| 579 | 130 | grid_LonB(I) = G%west_lon + dLon*REAL(I+I_offset-(G%isg-1)) |
| 580 | enddo | |
| 581 | 129 | do i=isd,ied |
| 582 | 129 | grid_LonT(i) = G%west_lon + dLon*(REAL(i+I_offset-G%isg)+0.5) |
| 583 | enddo | |
| 584 | ||
| 585 | 1 | dL_di = (G%len_lon * 4.0*atan(1.0)) / (180.0 * G%Domain%niglobal) |
| 586 | 8971 | do J=JsdB,JedB ; do I=IsdB,IedB |
| 587 | 8901 | G%geoLonBu(I,J) = grid_lonB(I) |
| 588 | 8901 | G%geoLatBu(I,J) = grid_latB(J) |
| 589 | ||
| 590 | ! The following line is needed to reproduce the solution from | |
| 591 | ! set_grid_metrics_mercator when used to generate a simple spherical grid. | |
| 592 | 8901 | G%dxBu(I,J) = G%Rad_Earth_L * COS( G%geoLatBu(I,J)*PI_180 ) * dL_di |
| 593 | ! G%dxBu(I,J) = G%Rad_Earth_L * dLon*PI_180 * COS( G%geoLatBu(I,J)*PI_180 ) | |
| 594 | 8901 | G%dyBu(I,J) = G%Rad_Earth_L * dLat*PI_180 |
| 595 | 8970 | G%areaBu(I,J) = G%dxBu(I,J) * G%dyBu(I,J) |
| 596 | enddo ; enddo | |
| 597 | ||
| 598 | 8902 | do J=JsdB,JedB ; do i=isd,ied |
| 599 | 8832 | G%geoLonCv(i,J) = grid_LonT(i) |
| 600 | 8832 | G%geoLatCv(i,J) = grid_latB(J) |
| 601 | ||
| 602 | ! The following line is needed to reproduce the solution from | |
| 603 | ! set_grid_metrics_mercator when used to generate a simple spherical grid. | |
| 604 | 8832 | G%dxCv(i,J) = G%Rad_Earth_L * COS( G%geoLatCv(i,J)*PI_180 ) * dL_di |
| 605 | ! G%dxCv(i,J) = G%Rad_Earth_L * (dLon*PI_180) * COS( G%geoLatCv(i,J)*PI_180 ) | |
| 606 | 8901 | G%dyCv(i,J) = G%Rad_Earth_L * dLat*PI_180 |
| 607 | enddo ; enddo | |
| 608 | ||
| 609 | 8841 | do j=jsd,jed ; do I=IsdB,IedB |
| 610 | 8772 | G%geoLonCu(I,j) = grid_lonB(I) |
| 611 | 8772 | G%geoLatCu(I,j) = grid_LatT(j) |
| 612 | ||
| 613 | ! The following line is needed to reproduce the solution from | |
| 614 | ! set_grid_metrics_mercator when used to generate a simple spherical grid. | |
| 615 | 8772 | G%dxCu(I,j) = G%Rad_Earth_L * COS( G%geoLatCu(I,j)*PI_180 ) * dL_di |
| 616 | ! G%dxCu(I,j) = G%Rad_Earth_L * dLon*PI_180 * COS( latitude ) | |
| 617 | 8840 | G%dyCu(I,j) = G%Rad_Earth_L * dLat*PI_180 |
| 618 | enddo ; enddo | |
| 619 | ||
| 620 | 8773 | do j=jsd,jed ; do i=isd,ied |
| 621 | 8704 | G%geoLonT(i,j) = grid_LonT(i) |
| 622 | 8704 | G%geoLatT(i,j) = grid_LatT(j) |
| 623 | ||
| 624 | ! The following line is needed to reproduce the solution from | |
| 625 | ! set_grid_metrics_mercator when used to generate a simple spherical grid. | |
| 626 | 8704 | G%dxT(i,j) = G%Rad_Earth_L * COS( G%geoLatT(i,j)*PI_180 ) * dL_di |
| 627 | ! G%dxT(i,j) = G%Rad_Earth_L * dLon*PI_180 * COS( latitude ) | |
| 628 | 8704 | G%dyT(i,j) = G%Rad_Earth_L * dLat*PI_180 |
| 629 | ||
| 630 | ! latitude = G%geoLatCv(i,J)*PI_180 ! In radians | |
| 631 | ! dL_di = G%geoLatCv(i,max(jsd,J-1))*PI_180 ! In radians | |
| 632 | ! G%areaT(i,j) = Rad_Earth_L**2*dLon*dLat*ABS(SIN(latitude)-SIN(dL_di)) | |
| 633 | 8772 | G%areaT(i,j) = G%dxT(i,j) * G%dyT(i,j) |
| 634 | enddo ; enddo | |
| 635 | ||
| 636 | 1 | call callTree_leave("set_grid_metrics_spherical()") |
| 637 | 1 | end subroutine set_grid_metrics_spherical |
| 638 | ||
| 639 | !> Calculate the values of the metric terms that might be used | |
| 640 | !! and save them in arrays. | |
| 641 | !! | |
| 642 | !! Within this subroutine, the x- and y- grid spacings and their | |
| 643 | !! inverses and the cell areas centered on h, q, u, and v points are | |
| 644 | !! calculated, as are the geographic locations of each of these 4 | |
| 645 | !! sets of points. | |
| 646 | 0 | subroutine set_grid_metrics_mercator(G, param_file, US) |
| 647 | type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type | |
| 648 | type(param_file_type), intent(in) :: param_file !< Parameter file structure | |
| 649 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 650 | ! Local variables | |
| 651 | integer :: i, j, isd, ied, jsd, jed | |
| 652 | integer :: I_off, J_off | |
| 653 | type(GPS) :: GP | |
| 654 | character(len=48) :: mdl = "MOM_grid_init set_grid_metrics_mercator" | |
| 655 | real :: PI, PI_2 ! PI = 3.1415926... as 4*atan(1), PI_2 = (PI) /2.0 [nondim] | |
| 656 | real :: y_q, y_h ! Latitudes of a point [radians] | |
| 657 | real :: id ! The i-grid space positions whose longitude is being sought [gridpoints] | |
| 658 | real :: jd ! The j-grid space positions whose latitude is being sought [gridpoints] | |
| 659 | real :: x_q, x_h ! Longitudes of a point [radians] | |
| 660 | real, dimension(G%isd:G%ied,G%jsd:G%jed) :: & | |
| 661 | 0 | xh, yh ! Latitude and longitude of h points in radians [radians] |
| 662 | real, dimension(G%IsdB:G%IedB,G%jsd:G%jed) :: & | |
| 663 | 0 | xu, yu ! Latitude and longitude of u points in radians [radians] |
| 664 | real, dimension(G%isd:G%ied,G%JsdB:G%JedB) :: & | |
| 665 | 0 | xv, yv ! Latitude and longitude of v points in radians [radians] |
| 666 | real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB) :: & | |
| 667 | 0 | xq, yq ! Latitude and longitude of q points in radians [radians] |
| 668 | real :: fnRef ! fnRef is the value of Int_dj_dy or | |
| 669 | ! Int_dj_dy at a latitude or longitude that is | |
| 670 | ! being set to be at grid index jRef or iRef [gridpoints] | |
| 671 | real :: jRef, iRef ! The grid index at which fnRef is evaluated [gridpoints] | |
| 672 | integer :: itt1, itt2 | |
| 673 | logical, parameter :: simple_area = .true. | |
| 674 | integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, IsdB, IedB, JsdB, JedB | |
| 675 | ||
| 676 | ! All of the metric terms should be defined over the domain from | |
| 677 | ! isd to ied. Outside of the physical domain, both the metrics | |
| 678 | ! and their inverses may be set to zero. | |
| 679 | 0 | is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec |
| 680 | 0 | isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed |
| 681 | 0 | Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB |
| 682 | 0 | IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB |
| 683 | 0 | I_off = G%idg_offset ; J_off = G%jdg_offset |
| 684 | ||
| 685 | 0 | GP%niglobal = G%Domain%niglobal |
| 686 | 0 | GP%njglobal = G%Domain%njglobal |
| 687 | ||
| 688 | 0 | call callTree_enter("set_grid_metrics_mercator(), MOM_grid_initialize.F90") |
| 689 | ||
| 690 | ! Calculate the values of the metric terms that might be used | |
| 691 | ! and save them in arrays. | |
| 692 | 0 | PI = 4.0*atan(1.0) ; PI_2 = 0.5*PI |
| 693 | ||
| 694 | call get_param(param_file, mdl, "SOUTHLAT", GP%south_lat, & | |
| 695 | "The southern latitude of the domain.", units="degrees_N", & | |
| 696 | 0 | fail_if_missing=.true.) |
| 697 | call get_param(param_file, mdl, "LENLAT", GP%len_lat, & | |
| 698 | "The latitudinal length of the domain.", units="degrees_N", & | |
| 699 | 0 | fail_if_missing=.true.) |
| 700 | call get_param(param_file, mdl, "WESTLON", GP%west_lon, & | |
| 701 | "The western longitude of the domain.", units="degrees_E", & | |
| 702 | 0 | default=0.0) |
| 703 | call get_param(param_file, mdl, "LENLON", GP%len_lon, & | |
| 704 | "The longitudinal length of the domain.", units="degrees_E", & | |
| 705 | 0 | fail_if_missing=.true.) |
| 706 | call get_param(param_file, mdl, "RAD_EARTH", GP%Rad_Earth_L, & | |
| 707 | 0 | "The radius of the Earth.", units="m", default=6.378e6, scale=US%m_to_L) |
| 708 | 0 | G%south_lat = GP%south_lat ; G%len_lat = GP%len_lat |
| 709 | 0 | G%west_lon = GP%west_lon ; G%len_lon = GP%len_lon |
| 710 | 0 | G%Rad_Earth_L = GP%Rad_Earth_L |
| 711 | ||
| 712 | call get_param(param_file, mdl, "ISOTROPIC", GP%isotropic, & | |
| 713 | "If true, an isotropic grid on a sphere (also known as "//& | |
| 714 | "a Mercator grid) is used. With an isotropic grid, the "//& | |
| 715 | "meridional extent of the domain (LENLAT), the zonal "//& | |
| 716 | "extent (LENLON), and the number of grid points in each "//& | |
| 717 | "direction are _not_ independent. In MOM the meridional "//& | |
| 718 | "extent is determined to fit the zonal extent and the "//& | |
| 719 | "number of grid points, while grid is perfectly isotropic.", & | |
| 720 | 0 | default=.false.) |
| 721 | call get_param(param_file, mdl, "EQUATOR_REFERENCE", GP%equator_reference, & | |
| 722 | "If true, the grid is defined to have the equator at the "//& | |
| 723 | "nearest q- or h- grid point to (-LOWLAT*NJGLOBAL/LENLAT).", & | |
| 724 | 0 | default=.true.) |
| 725 | call get_param(param_file, mdl, "LAT_ENHANCE_FACTOR", GP%Lat_enhance_factor, & | |
| 726 | "The amount by which the meridional resolution is "//& | |
| 727 | "enhanced within LAT_EQ_ENHANCE of the equator.", & | |
| 728 | 0 | units="nondim", default=1.0) |
| 729 | call get_param(param_file, mdl, "LAT_EQ_ENHANCE", GP%Lat_eq_enhance, & | |
| 730 | "The latitude range to the north and south of the equator "//& | |
| 731 | "over which the resolution is enhanced.", units="degrees_N", & | |
| 732 | 0 | default=0.0) |
| 733 | ||
| 734 | ! With an isotropic grid, the north-south extent of the domain, | |
| 735 | ! the east-west extent, and the number of grid points in each | |
| 736 | ! direction are _not_ independent. Here the north-south extent | |
| 737 | ! will be determined to fit the east-west extent and the number of | |
| 738 | ! grid points. The grid is perfectly isotropic. | |
| 739 | 0 | if (GP%equator_reference) then |
| 740 | ! With the following expression, the equator will always be placed | |
| 741 | ! on either h or q points, in a position consistent with the ratio | |
| 742 | ! GP%south_lat to GP%len_lat. | |
| 743 | 0 | jRef = (G%jsg-1) + 0.5*FLOOR(GP%njglobal*((-1.0*GP%south_lat*2.0)/GP%len_lat)+0.5) |
| 744 | 0 | fnRef = Int_dj_dy(0.0, GP) |
| 745 | else | |
| 746 | ! The following line sets the reference latitude GP%south_lat at j=js-1 (or -2?) | |
| 747 | 0 | jRef = (G%jsg-1) |
| 748 | 0 | fnRef = Int_dj_dy((GP%south_lat*PI/180.0), GP) |
| 749 | endif | |
| 750 | ||
| 751 | ! These calculations no longer depend on the order in which they | |
| 752 | ! are performed because they all use the same (poor) starting guess and | |
| 753 | ! iterate to convergence. | |
| 754 | ! Note that the dynamic grid always uses symmetric memory for the global | |
| 755 | ! arrays G%gridLatB and G%gridLonB. | |
| 756 | 0 | do J=G%jsg-1,G%jeg |
| 757 | 0 | jd = fnRef + (J - jRef) |
| 758 | 0 | y_q = find_root(Int_dj_dy, dy_dj, GP, jd, 0.0, -1.0*PI_2, PI_2, itt2) |
| 759 | 0 | G%gridLatB(J) = y_q*180.0/PI |
| 760 | ! if (is_root_pe()) & | |
| 761 | ! write(stdout, '("J, y_q = ",I0,", ",ES14.4," itts = ",I0)') j, y_q, itt2 | |
| 762 | enddo | |
| 763 | 0 | do j=G%jsg,G%jeg |
| 764 | 0 | jd = fnRef + (j - jRef) - 0.5 |
| 765 | 0 | y_h = find_root(Int_dj_dy, dy_dj, GP, jd, 0.0, -1.0*PI_2, PI_2, itt1) |
| 766 | 0 | G%gridLatT(j) = y_h*180.0/PI |
| 767 | ! if (is_root_pe()) & | |
| 768 | ! write(stdout, '("j, y_h = ",I0,", ",ES14.4," itts = ",I0)') j, y_h, itt1 | |
| 769 | enddo | |
| 770 | 0 | do J=JsdB+J_off,JedB+J_off |
| 771 | 0 | jd = fnRef + (J - jRef) |
| 772 | 0 | y_q = find_root(Int_dj_dy, dy_dj, GP, jd, 0.0, -1.0*PI_2, PI_2, itt2) |
| 773 | 0 | do I=IsdB,IedB ; yq(I,J-J_off) = y_q ; enddo |
| 774 | 0 | do i=isd,ied ; yv(i,J-J_off) = y_q ; enddo |
| 775 | enddo | |
| 776 | 0 | do j=jsd+J_off,jed+J_off |
| 777 | 0 | jd = fnRef + (j - jRef) - 0.5 |
| 778 | 0 | y_h = find_root(Int_dj_dy, dy_dj, GP, jd, 0.0, -1.0*PI_2, PI_2, itt1) |
| 779 | 0 | if ((j >= jsd+J_off) .and. (j <= jed+J_off)) then |
| 780 | 0 | do i=isd,ied ; yh(i,j-J_off) = y_h ; enddo |
| 781 | 0 | do I=IsdB,IedB ; yu(I,j-J_off) = y_h ; enddo |
| 782 | endif | |
| 783 | enddo | |
| 784 | ||
| 785 | ! Determine the longitudes of the various points. | |
| 786 | ||
| 787 | ! These two lines place the western edge of the domain at GP%west_lon. | |
| 788 | 0 | iRef = (G%isg-1) + GP%niglobal |
| 789 | 0 | fnRef = Int_di_dx(((GP%west_lon+GP%len_lon)*PI/180.0), GP) |
| 790 | ||
| 791 | ! These calculations no longer depend on the order in which they | |
| 792 | ! are performed because they all use the same (poor) starting guess and | |
| 793 | ! iterate to convergence. | |
| 794 | 0 | do I=G%isg-1,G%ieg |
| 795 | 0 | id = fnRef + (I - iRef) |
| 796 | 0 | x_q = find_root(Int_di_dx, dx_di, GP, id, 0.0, -4.0*PI, 4.0*PI, itt2) |
| 797 | 0 | G%gridLonB(I) = x_q*180.0/PI |
| 798 | enddo | |
| 799 | 0 | do i=G%isg,G%ieg |
| 800 | 0 | id = fnRef + (i - iRef) - 0.5 |
| 801 | 0 | x_h = find_root(Int_di_dx, dx_di, GP, id, 0.0, -4.0*PI, 4.0*PI, itt1) |
| 802 | 0 | G%gridLonT(i) = x_h*180.0/PI |
| 803 | enddo | |
| 804 | 0 | do I=IsdB+I_off,IedB+I_off |
| 805 | 0 | id = fnRef + (I - iRef) |
| 806 | 0 | x_q = find_root(Int_di_dx, dx_di, GP, id, 0.0, -4.0*PI, 4.0*PI, itt2) |
| 807 | 0 | do J=JsdB,JedB ; xq(I-I_off,J) = x_q ; enddo |
| 808 | 0 | do j=jsd,jed ; xu(I-I_off,j) = x_q ; enddo |
| 809 | enddo | |
| 810 | 0 | do i=isd+I_off,ied+I_off |
| 811 | 0 | id = fnRef + (i - iRef) - 0.5 |
| 812 | 0 | x_h = find_root(Int_di_dx, dx_di, GP, id, 0.0, -4.0*PI, 4.0*PI, itt1) |
| 813 | 0 | do j=jsd,jed ; xh(i-I_off,j) = x_h ; enddo |
| 814 | 0 | do J=JsdB,JedB ; xv(i-I_off,J) = x_h ; enddo |
| 815 | enddo | |
| 816 | ||
| 817 | 0 | do J=JsdB,JedB ; do I=IsdB,IedB |
| 818 | 0 | G%geoLonBu(I,J) = xq(I,J)*180.0/PI |
| 819 | 0 | G%geoLatBu(I,J) = yq(I,J)*180.0/PI |
| 820 | 0 | G%dxBu(I,J) = ds_di(xq(I,J), yq(I,J), GP) |
| 821 | 0 | G%dyBu(I,J) = ds_dj(xq(I,J), yq(I,J), GP) |
| 822 | ||
| 823 | 0 | G%areaBu(I,J) = G%dxBu(I,J) * G%dyBu(I,J) |
| 824 | 0 | G%IareaBu(I,J) = 1.0 / (G%areaBu(I,J)) |
| 825 | enddo ; enddo | |
| 826 | ||
| 827 | 0 | do j=jsd,jed ; do i=isd,ied |
| 828 | 0 | G%geoLonT(i,j) = xh(i,j)*180.0/PI |
| 829 | 0 | G%geoLatT(i,j) = yh(i,j)*180.0/PI |
| 830 | 0 | G%dxT(i,j) = ds_di(xh(i,j), yh(i,j), GP) |
| 831 | 0 | G%dyT(i,j) = ds_dj(xh(i,j), yh(i,j), GP) |
| 832 | ||
| 833 | 0 | G%areaT(i,j) = G%dxT(i,j)*G%dyT(i,j) |
| 834 | 0 | G%IareaT(i,j) = 1.0 / (G%areaT(i,j)) |
| 835 | enddo ; enddo | |
| 836 | ||
| 837 | 0 | do j=jsd,jed ; do I=IsdB,IedB |
| 838 | 0 | G%geoLonCu(I,j) = xu(I,j)*180.0/PI |
| 839 | 0 | G%geoLatCu(I,j) = yu(I,j)*180.0/PI |
| 840 | 0 | G%dxCu(I,j) = ds_di(xu(I,j), yu(I,j), GP) |
| 841 | 0 | G%dyCu(I,j) = ds_dj(xu(I,j), yu(I,j), GP) |
| 842 | enddo ; enddo | |
| 843 | ||
| 844 | 0 | do J=JsdB,JedB ; do i=isd,ied |
| 845 | 0 | G%geoLonCv(i,J) = xv(i,J)*180.0/PI |
| 846 | 0 | G%geoLatCv(i,J) = yv(i,J)*180.0/PI |
| 847 | 0 | G%dxCv(i,J) = ds_di(xv(i,J), yv(i,J), GP) |
| 848 | 0 | G%dyCv(i,J) = ds_dj(xv(i,J), yv(i,J), GP) |
| 849 | enddo ; enddo | |
| 850 | ||
| 851 | if (.not.simple_area) then | |
| 852 | do j=JsdB+1,jed ; do i=IsdB+1,ied | |
| 853 | G%areaT(I,J) = GP%Rad_Earth_L**2 * & | |
| 854 | (dL(xq(I-1,J-1),xq(I-1,J),yq(I-1,J-1),yq(I-1,J)) + & | |
| 855 | (dL(xq(I-1,J),xq(I,J),yq(I-1,J),yq(I,J)) + & | |
| 856 | (dL(xq(I,J),xq(I,J-1),yq(I,J),yq(I,J-1)) + & | |
| 857 | dL(xq(I,J-1),xq(I-1,J-1),yq(I,J-1),yq(I-1,J-1))))) | |
| 858 | enddo ; enddo | |
| 859 | if ((IsdB == isd) .or. (JsdB == jsq)) then | |
| 860 | ! Fill in row and column 1 to calculate the area in the southernmost | |
| 861 | ! and westernmost land cells when we are not using symmetric memory. | |
| 862 | ! The pass_var call updates these values if they are not land cells. | |
| 863 | G%areaT(isd+1,jsd) = G%areaT(isd+1,jsd+1) | |
| 864 | do j=jsd,jed ; G%areaT(isd,j) = G%areaT(isd+1,j) ; enddo | |
| 865 | do i=isd,ied ; G%areaT(i,jsd) = G%areaT(i,jsd+1) ; enddo | |
| 866 | ! Now replace the data in the halos, if value values exist. | |
| 867 | call pass_var(G%areaT,G%Domain) | |
| 868 | endif | |
| 869 | do j=jsd,jed ; do i=isd,ied | |
| 870 | G%IareaT(i,j) = 1.0 / (G%areaT(i,j)) | |
| 871 | enddo ; enddo | |
| 872 | endif | |
| 873 | ||
| 874 | 0 | call callTree_leave("set_grid_metrics_mercator()") |
| 875 | 0 | end subroutine set_grid_metrics_mercator |
| 876 | ||
| 877 | ||
| 878 | !> This function returns the grid spacing in the logical x direction in [L ~> m]. | |
| 879 | 0 | function ds_di(x, y, GP) |
| 880 | real, intent(in) :: x !< The longitude in question [radians] | |
| 881 | real, intent(in) :: y !< The latitude in question [radians] | |
| 882 | type(GPS), intent(in) :: GP !< A structure of grid parameters | |
| 883 | ||
| 884 | real :: ds_di ! The returned grid spacing [L ~> m] | |
| 885 | ||
| 886 | 0 | ds_di = GP%Rad_Earth_L * cos(y) * dx_di(x,GP) |
| 887 | ! In general, this might be... | |
| 888 | ! ds_di = GP%Rad_Earth_L * sqrt( cos(y)*cos(y) * dx_di(x,y,GP)*dx_di(x,y,GP) + & | |
| 889 | ! dy_di(x,y,GP)*dy_di(x,y,GP)) | |
| 890 | 0 | end function ds_di |
| 891 | ||
| 892 | !> This function returns the grid spacing in the logical y direction in [L ~> m]. | |
| 893 | 0 | function ds_dj(x, y, GP) |
| 894 | real, intent(in) :: x !< The longitude in question [radians] | |
| 895 | real, intent(in) :: y !< The latitude in question [radians] | |
| 896 | type(GPS), intent(in) :: GP !< A structure of grid parameters | |
| 897 | ||
| 898 | real :: ds_dj ! The returned grid spacing [L ~> m] | |
| 899 | ||
| 900 | 0 | ds_dj = GP%Rad_Earth_L * dy_dj(y,GP) |
| 901 | ! In general, this might be... | |
| 902 | ! ds_dj = GP%Rad_Earth_L * sqrt( cos(y)*cos(y) * dx_dj(x,y,GP)*dx_dj(x,y,GP) + & | |
| 903 | ! dy_dj(x,y,GP)*dy_dj(x,y,GP)) | |
| 904 | 0 | end function ds_dj |
| 905 | ||
| 906 | !> This function returns the contribution from the line integral along one of the four sides of a | |
| 907 | !! cell face to the area of a cell, in [radians2], assuming that the sides follow a linear path in | |
| 908 | !! latitude and longitude (i.e., on a Mercator grid). | |
| 909 | 0 | function dL(x1, x2, y1, y2) |
| 910 | real, intent(in) :: x1 !< Segment starting longitude [radians] | |
| 911 | real, intent(in) :: x2 !< Segment ending longitude [radians] | |
| 912 | real, intent(in) :: y1 !< Segment starting latitude [radians] | |
| 913 | real, intent(in) :: y2 !< Segment ending latitude [radians] | |
| 914 | ! Local variables | |
| 915 | real :: dL ! A contribution to the spanned area the surface of the sphere [radian2] | |
| 916 | real :: r ! A contribution from the range of latitudes, including trigonometric factors [radians] | |
| 917 | real :: dy ! The spanned range of latitudes [radians] | |
| 918 | ||
| 919 | 0 | dy = y2 - y1 |
| 920 | ||
| 921 | 0 | if (ABS(dy) > 2.5e-8) then |
| 922 | 0 | r = ((1.0 - cos(dy))*cos(y1) + sin(dy)*sin(y1)) / dy |
| 923 | else | |
| 924 | 0 | r = (0.5*dy*cos(y1) + sin(y1)) |
| 925 | endif | |
| 926 | 0 | dL = r * (x2 - x1) |
| 927 | ||
| 928 | 0 | end function dL |
| 929 | ||
| 930 | !> This subroutine finds and returns the value of y at which the monotonically increasing | |
| 931 | !! function fn takes the value fnval, also returning in ittmax the number of iterations of | |
| 932 | !! Newton's method that were used to polish the root. | |
| 933 | 0 | function find_root( fn, dy_df, GP, fnval, y1, ymin, ymax, ittmax) |
| 934 | real :: find_root !< The value of y where fn(y) = fnval that will be returned [radians] | |
| 935 | real, external :: fn !< The external function whose root is being sought [gridpoints] | |
| 936 | real, external :: dy_df !< The inverse of the derivative of that function [radian gridpoint-1] | |
| 937 | type(GPS), intent(in) :: GP !< A structure of grid parameters | |
| 938 | real, intent(in) :: fnval !< The value of fn being sought [gridpoints] | |
| 939 | real, intent(in) :: y1 !< A first guess for y [radians] | |
| 940 | real, intent(in) :: ymin !< The minimum permitted value of y [radians] | |
| 941 | real, intent(in) :: ymax !< The maximum permitted value of y [radians] | |
| 942 | integer, intent(out) :: ittmax !< The number of iterations used to polish the root | |
| 943 | ! Local variables | |
| 944 | real :: y, y_next ! Successive guesses at the root position [radians] | |
| 945 | real :: ybot, ytop ! Brackets bounding the root [radians] | |
| 946 | real :: fnbot, fntop ! Values of fn at the bounding values of y [gridpoints] | |
| 947 | real :: dy_dfn ! The inverse of the local derivative of fn with y [radian gridpoint-1] | |
| 948 | real :: dy ! The jump to the next guess of y [radians] | |
| 949 | real :: fny ! The difference between fn(y) and the target value [gridpoints] | |
| 950 | integer :: itt | |
| 951 | character(len=256) :: warnmesg | |
| 952 | ||
| 953 | ! Bracket the root. Do not use the bounding values because the value at the | |
| 954 | ! function at the bounds could be infinite, as is the case for the Mercator | |
| 955 | ! grid recursion relation. (I.e., this is a search on an open interval.) | |
| 956 | 0 | ybot = y1 |
| 957 | 0 | fnbot = fn(ybot,GP) - fnval ; itt = 0 |
| 958 | 0 | do while (fnbot > 0.0) |
| 959 | 0 | if ((ybot - 2.0*dy_df(ybot,GP)) < (0.5*(ybot+ymin))) then |
| 960 | ! Go twice as far as the secant method would normally go. | |
| 961 | 0 | ybot = ybot - 2.0*dy_df(ybot,GP) |
| 962 | else ! But stay within the open interval! | |
| 963 | 0 | ybot = 0.5*(ybot+ymin) ; itt = itt + 1 |
| 964 | endif | |
| 965 | 0 | fnbot = fn(ybot,GP) - fnval |
| 966 | ||
| 967 | 0 | if ((itt > 50) .and. (fnbot > 0.0)) then |
| 968 | write(warnmesg, '("PE ",I0," unable to find bottom bound for grid function. & | |
| 969 | &x = ",ES10.4,", xmax = ",ES10.4,", fn = ",ES10.4,", dfn_dx = ",ES10.4,& | |
| 970 | &", seeking fn = ",ES10.4," - fn = ",ES10.4,".")') & | |
| 971 | 0 | pe_here(),ybot,ymin,fn(ybot,GP),dy_df(ybot,GP),fnval, fnbot |
| 972 | 0 | call MOM_error(FATAL,warnmesg) |
| 973 | endif | |
| 974 | enddo | |
| 975 | ||
| 976 | 0 | ytop = y1 |
| 977 | 0 | fntop = fn(ytop,GP) - fnval ; itt = 0 |
| 978 | 0 | do while (fntop < 0.0) |
| 979 | 0 | if ((ytop + 2.0*dy_df(ytop,GP)) < (0.5*(ytop+ymax))) then |
| 980 | ! Go twice as far as the secant method would normally go. | |
| 981 | 0 | ytop = ytop + 2.0*dy_df(ytop,GP) |
| 982 | else ! But stay within the open interval! | |
| 983 | 0 | ytop = 0.5*(ytop+ymax) ; itt = itt + 1 |
| 984 | endif | |
| 985 | 0 | fntop = fn(ytop,GP) - fnval |
| 986 | ||
| 987 | 0 | if ((itt > 50) .and. (fntop < 0.0)) then |
| 988 | write(warnmesg, '("PE ",I0," unable to find top bound for grid function. & | |
| 989 | &x = ",ES10.4,", xmax = ",ES10.4,", fn = ",ES10.4,", dfn_dx = ",ES10.4, & | |
| 990 | &", seeking fn = ",ES10.4," - fn = ",ES10.4,".")') & | |
| 991 | 0 | pe_here(),ytop,ymax,fn(ytop,GP),dy_df(ytop,GP),fnval,fntop |
| 992 | 0 | call MOM_error(FATAL,warnmesg) |
| 993 | endif | |
| 994 | enddo | |
| 995 | ||
| 996 | ! Find the root using a bracketed variant of Newton's method, starting | |
| 997 | ! with a false-positon method first guess. | |
| 998 | 0 | if ((fntop < 0.0) .or. (fnbot > 0.0) .or. (ytop < ybot)) then |
| 999 | write(warnmesg, '("PE ",I0," find_root failed to bracket function. y = ",& | |
| 1000 | 0 | &2ES10.4,", fn = ",2ES10.4,".")') pe_here(),ybot,ytop,fnbot,fntop |
| 1001 | 0 | call MOM_error(FATAL, warnmesg) |
| 1002 | endif | |
| 1003 | ||
| 1004 | 0 | if (fntop == 0.0) then ; y = ytop ; fny = fntop |
| 1005 | 0 | elseif (fnbot == 0.0) then ; y = ybot ; fny = fnbot |
| 1006 | else | |
| 1007 | 0 | y = (ybot*fntop - ytop*fnbot) / (fntop - fnbot) |
| 1008 | 0 | fny = fn(y,GP) - fnval |
| 1009 | 0 | if (fny < 0.0) then ; fnbot = fny ; ybot = y |
| 1010 | 0 | else ; fntop = fny ; ytop = y ; endif |
| 1011 | endif | |
| 1012 | ||
| 1013 | 0 | do itt=1,50 |
| 1014 | 0 | dy_dfn = dy_df(y,GP) |
| 1015 | ||
| 1016 | 0 | dy = -1.0* fny * dy_dfn |
| 1017 | 0 | y_next = y + dy |
| 1018 | 0 | if ((y_next >= ytop) .or. (y_next <= ybot)) then |
| 1019 | ! The Newton's method estimate has escaped bracketing, so use the | |
| 1020 | ! false-position method instead. The complicated test is to properly | |
| 1021 | ! handle the case where the iteration is down to roundoff level differences. | |
| 1022 | 0 | y_next = y |
| 1023 | 0 | if (abs(fntop - fnbot) > EPSILON(y) * (abs(fntop) + abs(fnbot))) & |
| 1024 | 0 | y_next = (ybot*fntop - ytop*fnbot) / (fntop - fnbot) |
| 1025 | endif | |
| 1026 | ||
| 1027 | 0 | dy = y_next - y |
| 1028 | 0 | if (ABS(dy) < (2.0*EPSILON(y)*(ABS(y) + ABS(y_next)) + 1.0e-20)) then |
| 1029 | 0 | y = y_next ; exit |
| 1030 | endif | |
| 1031 | 0 | y = y_next |
| 1032 | ||
| 1033 | 0 | fny = fn(y,GP) - fnval |
| 1034 | 0 | if (fny > 0.0) then ; ytop = y ; fntop = fny |
| 1035 | 0 | elseif (fny < 0.0) then ; ybot = y ; fnbot = fny |
| 1036 | 0 | else ; exit ; endif |
| 1037 | ||
| 1038 | enddo | |
| 1039 | 0 | if (ABS(y) < 1e-12) y = 0.0 |
| 1040 | ||
| 1041 | 0 | ittmax = itt |
| 1042 | 0 | find_root = y |
| 1043 | 0 | end function find_root |
| 1044 | ||
| 1045 | !> This function calculates and returns the value of dx/di in [radian gridpoint-1], | |
| 1046 | !! where x is the longitude in Radians, and i is the integral east-west grid index. | |
| 1047 | 0 | function dx_di(x, GP) |
| 1048 | real, intent(in) :: x !< The longitude in question [radians] | |
| 1049 | type(GPS), intent(in) :: GP !< A structure of grid parameters | |
| 1050 | real :: dx_di ! The derivative of zonal position with the grid index [radian gridpoint-1] | |
| 1051 | ||
| 1052 | 0 | dx_di = (GP%len_lon * 4.0*atan(1.0)) / (180.0 * GP%niglobal) |
| 1053 | ||
| 1054 | 0 | end function dx_di |
| 1055 | ||
| 1056 | !> This function calculates and returns the integral of the inverse | |
| 1057 | !! of dx/di to the point x, in radians [gridpoints] | |
| 1058 | 0 | function Int_di_dx(x, GP) |
| 1059 | real, intent(in) :: x !< The longitude in question [radians] | |
| 1060 | type(GPS), intent(in) :: GP !< A structure of grid parameters | |
| 1061 | real :: Int_di_dx ! A position in the global i-index space [gridpoints] | |
| 1062 | ||
| 1063 | 0 | Int_di_dx = x * ((180.0 * GP%niglobal) / (GP%len_lon * 4.0*atan(1.0))) |
| 1064 | ||
| 1065 | 0 | end function Int_di_dx |
| 1066 | ||
| 1067 | !> This subroutine calculates and returns the value of dy/dj in [radian gridpoint-1], | |
| 1068 | !! where y is the latitude in Radians, and j is the integral north-south grid index. | |
| 1069 | 0 | function dy_dj(y, GP) |
| 1070 | real, intent(in) :: y !< The latitude in question [radians] | |
| 1071 | type(GPS), intent(in) :: GP !< A structure of grid parameters | |
| 1072 | real :: dy_dj ! The derivative of meridional position with the grid index [radian gridpoint-1] | |
| 1073 | ! Local variables | |
| 1074 | real :: PI ! 3.1415926... calculated as 4*atan(1) [nondim] | |
| 1075 | real :: C0 ! The constant that converts the nominal y-spacing in | |
| 1076 | ! gridpoints to the nominal spacing in Radians [radian gridpoint-1] | |
| 1077 | real :: y_eq_enhance ! The latitude in radians within which the resolution | |
| 1078 | ! is enhanced [radians] | |
| 1079 | 0 | PI = 4.0*atan(1.0) |
| 1080 | 0 | if (GP%isotropic) then |
| 1081 | 0 | C0 = (GP%len_lon * PI) / (180.0 * GP%niglobal) |
| 1082 | 0 | y_eq_enhance = PI*abs(GP%lat_eq_enhance)/180.0 |
| 1083 | 0 | if (ABS(y) < y_eq_enhance) then |
| 1084 | dy_dj = C0 * (cos(y) / (1.0 + 0.5*cos(y) * (GP%lat_enhance_factor - 1.0) * & | |
| 1085 | 0 | (1.0+cos(PI*y/y_eq_enhance)) )) |
| 1086 | else | |
| 1087 | 0 | dy_dj = C0 * cos(y) |
| 1088 | endif | |
| 1089 | else | |
| 1090 | 0 | C0 = (GP%len_lat * PI) / (180.0 * GP%njglobal) |
| 1091 | 0 | dy_dj = C0 |
| 1092 | endif | |
| 1093 | ||
| 1094 | 0 | end function dy_dj |
| 1095 | ||
| 1096 | !> This subroutine calculates and returns the integral of the inverse | |
| 1097 | !! of dy/dj to the point y in radians [gridpoints] | |
| 1098 | 0 | function Int_dj_dy(y, GP) |
| 1099 | real, intent(in) :: y !< The latitude in question [radians] | |
| 1100 | type(GPS), intent(in) :: GP !< A structure of grid parameters | |
| 1101 | real :: Int_dj_dy ! The grid position of latitude y [gridpoints] | |
| 1102 | ! Local variables | |
| 1103 | real :: I_C0 ! The inverse of the constant that converts the | |
| 1104 | ! nominal spacing in gridpoints to the nominal | |
| 1105 | ! spacing in Radians [gridpoint radian-1] | |
| 1106 | real :: PI ! 3.1415926... calculated as 4*atan(1) [nondim] | |
| 1107 | real :: y_eq_enhance ! The latitude in radians from from the equator within which the meridional | |
| 1108 | ! grid spacing is enhanced by a factor of GP%lat_enhance_factor [radians] | |
| 1109 | real :: r ! The y grid position in the global index space [gridpoints] | |
| 1110 | ||
| 1111 | 0 | PI = 4.0*atan(1.0) |
| 1112 | 0 | if (GP%isotropic) then |
| 1113 | 0 | I_C0 = (180.0 * GP%niglobal) / (GP%len_lon * PI) |
| 1114 | 0 | y_eq_enhance = PI*ABS(GP%lat_eq_enhance)/180.0 |
| 1115 | ||
| 1116 | 0 | if (y >= 0.0) then |
| 1117 | 0 | r = I_C0 * log((1.0 + sin(y))/cos(y)) |
| 1118 | else | |
| 1119 | 0 | r = -1.0 * I_C0 * log((1.0 - sin(y))/cos(y)) |
| 1120 | endif | |
| 1121 | ||
| 1122 | 0 | if (y >= y_eq_enhance) then |
| 1123 | 0 | r = r + I_C0*0.5*(GP%lat_enhance_factor - 1.0)*y_eq_enhance |
| 1124 | 0 | elseif (y <= -y_eq_enhance) then |
| 1125 | 0 | r = r - I_C0*0.5*(GP%lat_enhance_factor - 1.0)*y_eq_enhance |
| 1126 | else | |
| 1127 | r = r + I_C0*0.5*(GP%lat_enhance_factor - 1.0) * & | |
| 1128 | 0 | (y + (y_eq_enhance/PI)*sin(PI*y/y_eq_enhance)) |
| 1129 | endif | |
| 1130 | else | |
| 1131 | 0 | I_C0 = (180.0 * GP%njglobal) / (GP%len_lat * PI) |
| 1132 | 0 | r = I_C0 * y |
| 1133 | endif | |
| 1134 | ||
| 1135 | 0 | Int_dj_dy = r |
| 1136 | 0 | end function Int_dj_dy |
| 1137 | ||
| 1138 | !> Extrapolates missing metric data into all the halo regions. | |
| 1139 | 0 | subroutine extrapolate_metric(var, jh, missing) |
| 1140 | real, dimension(:,:), intent(inout) :: var !< The array in which to fill in halos in arbitrary units [A] | |
| 1141 | integer, intent(in) :: jh !< The size of the halos to be filled | |
| 1142 | real, optional, intent(in) :: missing !< The missing data fill value, 0 by default [A] | |
| 1143 | ! Local variables | |
| 1144 | real :: badval ! A bad data value [A] | |
| 1145 | integer :: i, j | |
| 1146 | ||
| 1147 | 0 | badval = 0.0 ; if (present(missing)) badval = missing |
| 1148 | ||
| 1149 | ! Fill in southern halo by extrapolating from the computational domain | |
| 1150 | 0 | do j=lbound(var,2)+jh,lbound(var,2),-1 ; do i=lbound(var,1),ubound(var,1) |
| 1151 | 0 | if (var(i,j)==badval) var(i,j) = 2.0*var(i,j+1)-var(i,j+2) |
| 1152 | enddo ; enddo | |
| 1153 | ||
| 1154 | ! Fill in northern halo by extrapolating from the computational domain | |
| 1155 | 0 | do j=ubound(var,2)-jh,ubound(var,2) ; do i=lbound(var,1),ubound(var,1) |
| 1156 | 0 | if (var(i,j)==badval) var(i,j) = 2.0*var(i,j-1)-var(i,j-2) |
| 1157 | enddo ; enddo | |
| 1158 | ||
| 1159 | ! Fill in western halo by extrapolating from the computational domain | |
| 1160 | 0 | do j=lbound(var,2),ubound(var,2) ; do i=lbound(var,1)+jh,lbound(var,1),-1 |
| 1161 | 0 | if (var(i,j)==badval) var(i,j) = 2.0*var(i+1,j)-var(i+2,j) |
| 1162 | enddo ; enddo | |
| 1163 | ||
| 1164 | ! Fill in eastern halo by extrapolating from the computational domain | |
| 1165 | 0 | do j=lbound(var,2),ubound(var,2) ; do i=ubound(var,1)-jh,ubound(var,1) |
| 1166 | 0 | if (var(i,j)==badval) var(i,j) = 2.0*var(i-1,j)-var(i-2,j) |
| 1167 | enddo ; enddo | |
| 1168 | ||
| 1169 | 0 | end subroutine extrapolate_metric |
| 1170 | ||
| 1171 | !> This function implements Adcroft's rule for reciprocals, namely that | |
| 1172 | !! Adcroft_Inv(x) = 1/x for |x|>0 or 0 for x=0. | |
| 1173 | 17604 | function Adcroft_reciprocal(val) result(I_val) |
| 1174 | real, intent(in) :: val !< The value being inverted in arbitrary units [A] | |
| 1175 | real :: I_val !< The Adcroft reciprocal of val [A-1] | |
| 1176 | ||
| 1177 | 17604 | I_val = 0.0 |
| 1178 | 17604 | if (val /= 0.0) I_val = 1.0/val |
| 1179 | 17604 | end function Adcroft_reciprocal |
| 1180 | ||
| 1181 | !> Initializes the grid masks and any metrics that come with masks already applied. | |
| 1182 | !! | |
| 1183 | !! Initialize_masks sets mask2dT, mask2dCu, mask2dCv, and mask2dBu to mask out | |
| 1184 | !! flow over any points which are shallower than Dmask and permit an | |
| 1185 | !! appropriate treatment of the boundary conditions. mask2dCu and mask2dCv | |
| 1186 | !! are 0.0 at any points adjacent to a land point. mask2dBu is 0.0 at | |
| 1187 | !! any land or boundary point. For points in the ocean interior or at open boundary | |
| 1188 | !! condition points, mask2dCu, mask2dCv, and mask2dBu are all 1.0. | |
| 1189 | 0 | subroutine initialize_masks(G, PF, US, OBC_dir_u, OBC_dir_v, open_corner_OBCs, maskT) |
| 1190 | type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type | |
| 1191 | type(param_file_type), intent(in) :: PF !< Parameter file structure | |
| 1192 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 1193 | integer, dimension(G%IsdB:G%IedB,G%jsd:G%jed), & | |
| 1194 | optional, intent(in) :: OBC_dir_u !< Trinary values that indicate whether there | |
| 1195 | !! is an open boundary condition at zonal velocity | |
| 1196 | !! faces and their orientation, with 0 for no OBC, | |
| 1197 | !! a positive value for an Eastern OBC and | |
| 1198 | !! a negative value for a Western OBC. | |
| 1199 | integer, dimension(G%isd:G%ied,G%JsdB:G%JedB), & | |
| 1200 | optional, intent(in) :: OBC_dir_v !< Trinary values that indicate whether there | |
| 1201 | !! is an open boundary condition at zonal velocity | |
| 1202 | !! faces and their orientation, with 0 for no OBC, | |
| 1203 | !! a positive value for a Northern OBC and | |
| 1204 | !! a negative value for a Southern OBC. | |
| 1205 | logical, optional, intent(in) :: open_corner_OBCs !< If present and true, the bay-like corner | |
| 1206 | !! between two orthogonal open boundary segments is open, | |
| 1207 | !! otherwise it is closed. | |
| 1208 | real, dimension(G%isd:G%ied,G%jsd:G%jed), & | |
| 1209 | optional, intent(in) :: maskT !< If present, this array is used to set the | |
| 1210 | !! the mask at tracer points instead of using the | |
| 1211 | !! bathymetry to determine the masks [nondim] | |
| 1212 | ||
| 1213 | ! Local variables | |
| 1214 | real :: Dmask ! The depth for masking in the same units as G%bathyT [Z ~> m]. | |
| 1215 | real :: min_depth ! The minimum ocean depth in the same units as G%bathyT [Z ~> m]. | |
| 1216 | real :: mask_depth ! The depth shallower than which to mask a point as land [Z ~> m]. | |
| 1217 | logical :: open_corners ! If true, the bay-like corner between two orthogonal open boundary segments is open | |
| 1218 | character(len=40) :: mdl = "MOM_grid_init initialize_masks" | |
| 1219 | integer :: i, j | |
| 1220 | ||
| 1221 | 1 | call callTree_enter("initialize_masks(), MOM_grid_initialize.F90") |
| 1222 | ||
| 1223 | call get_param(PF, mdl, "MINIMUM_DEPTH", min_depth, & | |
| 1224 | "If MASKING_DEPTH is unspecified, then anything shallower than "//& | |
| 1225 | "MINIMUM_DEPTH is assumed to be land and all fluxes are masked out. "//& | |
| 1226 | "If MASKING_DEPTH is specified, then all depths shallower than "//& | |
| 1227 | "MINIMUM_DEPTH but deeper than MASKING_DEPTH are rounded to MINIMUM_DEPTH.", & | |
| 1228 | 1 | units="m", default=0.0, scale=US%m_to_Z) |
| 1229 | call get_param(PF, mdl, "MASKING_DEPTH", mask_depth, & | |
| 1230 | "The depth below which to mask points as land points, for which all "//& | |
| 1231 | "fluxes are zeroed out. MASKING_DEPTH is ignored if it has the special "//& | |
| 1232 | "default value.", & | |
| 1233 | 1 | units="m", default=-9999.0, scale=US%m_to_Z, do_not_log=present(maskT)) |
| 1234 | ||
| 1235 | 1 | Dmask = mask_depth |
| 1236 | 1 | if (mask_depth == -9999.0*US%m_to_Z) Dmask = min_depth |
| 1237 | ||
| 1238 | 1 | open_corners = .false. ; if (present(open_corner_OBCs)) open_corners = open_corner_OBCs |
| 1239 | ||
| 1240 | 35484 | G%mask2dT(:,:) = 0.0 ; G%mask2dCu(:,:) = 0.0 ; G%mask2dCv(:,:) = 0.0 ; G%mask2dBu(:,:) = 0.0 |
| 1241 | ||
| 1242 | ! Construct the h-point or T-point mask | |
| 1243 | 1 | if (present(maskT)) then |
| 1244 | 0 | do j=G%jsd,G%jed ; do i=G%isd,G%ied |
| 1245 | 0 | G%mask2dT(i,j) = max(min(maskT(i,j), 1.0), 0.0) |
| 1246 | enddo ; enddo | |
| 1247 | else | |
| 1248 | 8773 | do j=G%jsd,G%jed ; do i=G%isd,G%ied |
| 1249 | 8772 | if (G%bathyT(i,j) <= Dmask) then |
| 1250 | 3534 | G%mask2dT(i,j) = 0.0 |
| 1251 | else | |
| 1252 | 5170 | G%mask2dT(i,j) = 1.0 |
| 1253 | endif | |
| 1254 | enddo ; enddo | |
| 1255 | endif | |
| 1256 | ||
| 1257 | 1 | call pass_var(G%mask2dT, G%Domain) |
| 1258 | ||
| 1259 | 8705 | do j=G%jsd,G%jed ; do I=G%isd,G%ied-1 |
| 1260 | 8704 | G%mask2dCu(I,j) = G%mask2dT(i,j) * G%mask2dT(i+1,j) |
| 1261 | enddo ; enddo | |
| 1262 | ||
| 1263 | 1 | if (present(OBC_dir_u)) then |
| 1264 | 0 | do j=G%jsd,G%jed ; do I=G%isd,G%ied-1 |
| 1265 | 0 | if (OBC_dir_u(I,j) > 0) G%mask2dCu(I,j) = G%mask2dT(i,j) |
| 1266 | 0 | if (OBC_dir_u(I,j) < 0) G%mask2dCu(I,j) = G%mask2dT(i+1,j) |
| 1267 | enddo ; enddo | |
| 1268 | endif | |
| 1269 | ||
| 1270 | 8705 | do j=G%jsd,G%jed ; do I=G%isd,G%ied-1 |
| 1271 | ! This mask may be revised later after the open boundary positions are specified. | |
| 1272 | 8704 | G%OBCmaskCu(I,j) = G%mask2dCu(I,j) |
| 1273 | enddo ; enddo | |
| 1274 | ||
| 1275 | 8644 | do J=G%jsd,G%jed-1 ; do i=G%isd,G%ied |
| 1276 | 8643 | G%mask2dCv(i,J) = G%mask2dT(i,j) * G%mask2dT(i,j+1) |
| 1277 | enddo ; enddo | |
| 1278 | ||
| 1279 | 1 | if (present(OBC_dir_v)) then |
| 1280 | 0 | do J=G%jsd,G%jed-1 ; do i=G%isd,G%ied |
| 1281 | 0 | if (OBC_dir_v(i,J) > 0) G%mask2dCv(i,J) = G%mask2dT(i,j) |
| 1282 | 0 | if (OBC_dir_v(i,J) < 0) G%mask2dCv(i,J) = G%mask2dT(i,j+1) |
| 1283 | enddo ; enddo | |
| 1284 | endif | |
| 1285 | ||
| 1286 | 8644 | do J=G%jsd,G%jed-1 ; do i=G%isd,G%ied |
| 1287 | ! This mask may be revised later after the open boundary positions are specified. | |
| 1288 | 8643 | G%OBCmaskCv(i,J) = G%mask2dCv(i,J) |
| 1289 | enddo ; enddo | |
| 1290 | ||
| 1291 | ! The mask at the vertex can be determined from the masks at the faces. | |
| 1292 | ! This works at interior ocean points or at convex OBC points. | |
| 1293 | 8577 | do J=G%jsd,G%jed-1 ; do I=G%isd,G%ied-1 |
| 1294 | 8576 | G%mask2dBu(I,J) = (G%mask2dCu(I,j) * G%mask2dCu(I,j+1)) * (G%mask2dCv(i,J) * G%mask2dCv(i+1,J)) |
| 1295 | enddo ; enddo | |
| 1296 | ||
| 1297 | ! This block resets masks at the vertices when there are OBCs. The right logic is that if there | |
| 1298 | ! are 2 or more unmasked OBCs, this point should be open, but to recreate the previous answers, | |
| 1299 | 1 | if (present(OBC_dir_u)) then |
| 1300 | 0 | do J=G%jsd,G%jed-1 ; do I=G%isd,G%ied-1 |
| 1301 | ! These are conditions to set open vertex points on a straight north-south coastline | |
| 1302 | 0 | if ((G%mask2dCu(I,j) * OBC_dir_u(I,j)) * (G%mask2dCu(I,j+1) * OBC_dir_u(I,j+1)) > 0.) & |
| 1303 | 0 | G%mask2dBu(I,J) = 1.0 |
| 1304 | enddo ; enddo | |
| 1305 | endif | |
| 1306 | 1 | if (present(OBC_dir_v)) then |
| 1307 | 0 | do J=G%jsd,G%jed-1 ; do I=G%isd,G%ied-1 |
| 1308 | ! These are conditions to set open vertex points on a straight east-west coastline | |
| 1309 | 0 | if ((G%mask2dCv(i,J) * OBC_dir_v(i,J)) * (G%mask2dCv(i+1,J) * OBC_dir_v(i+1,J)) > 0.) & |
| 1310 | 0 | G%mask2dBu(I,J) = 1.0 |
| 1311 | enddo ; enddo | |
| 1312 | endif | |
| 1313 | 1 | if (open_corners .and. present(OBC_dir_u) .and. present(OBC_dir_v)) then |
| 1314 | 0 | do J=G%jsd,G%jed-1 ; do I=G%isd,G%ied-1 |
| 1315 | ! These are the 4 conditions to set an open point in a concave (bay-like) corner | |
| 1316 | 0 | if ((G%mask2dCu(I,j+1) * OBC_dir_u(I,j+1) < 0.) .and. (G%mask2dCv(i+1,J) * OBC_dir_v(i+1,J) < 0.)) & |
| 1317 | 0 | G%mask2dBu(I,J) = 1.0 ! Southwestern corner |
| 1318 | 0 | if ((G%mask2dCu(I,j+1) * OBC_dir_u(I,j+1) > 0.) .and. (G%mask2dCv(i,J) * OBC_dir_v(i,J) < 0.)) & |
| 1319 | 0 | G%mask2dBu(I,J) = 1.0 ! Southeastern corner |
| 1320 | 0 | if ((G%mask2dCu(I,j) * OBC_dir_u(I,j) < 0.) .and. (G%mask2dCv(i+1,J) * OBC_dir_v(i+1,J) > 0.)) & |
| 1321 | 0 | G%mask2dBu(I,J) = 1.0 ! Northwestern corner |
| 1322 | 0 | if ((G%mask2dCu(I,j) * OBC_dir_u(I,j) > 0.) .and. (G%mask2dCv(i,J) * OBC_dir_v(i,J) > 0.)) & |
| 1323 | 0 | G%mask2dBu(I,J) = 1.0 ! Northeastern corner |
| 1324 | enddo ; enddo | |
| 1325 | endif | |
| 1326 | ||
| 1327 | 1 | call pass_var(G%mask2dBu, G%Domain, position=CORNER) |
| 1328 | 1 | call pass_vector(G%mask2dCu, G%mask2dCv, G%Domain, To_All+Scalar_Pair, CGRID_NE) |
| 1329 | ||
| 1330 | 8841 | do j=G%jsd,G%jed ; do I=G%IsdB,G%IedB |
| 1331 | ! This open face length may be revised later. | |
| 1332 | 8772 | G%dy_Cu(I,j) = G%mask2dCu(I,j) * G%dyCu(I,j) |
| 1333 | 8772 | G%IdxCu_OBCmask(I,j) = G%OBCmaskCu(I,j) * G%IdxCu(I,j) |
| 1334 | 8772 | G%areaCu(I,j) = G%dxCu(I,j) * G%dy_Cu(I,j) |
| 1335 | 8840 | G%IareaCu(I,j) = G%mask2dCu(I,j) * Adcroft_reciprocal(G%areaCu(I,j)) |
| 1336 | enddo ; enddo | |
| 1337 | ||
| 1338 | 8902 | do J=G%JsdB,G%JedB ; do i=G%isd,G%ied |
| 1339 | ! This open face length may be revised later. | |
| 1340 | 8832 | G%dx_Cv(i,J) = G%mask2dCv(i,J) * G%dxCv(i,J) |
| 1341 | 8832 | G%IdyCv_OBCmask(i,J) = G%OBCmaskCv(i,J) * G%IdyCv(i,J) |
| 1342 | 8832 | G%areaCv(i,J) = G%dyCv(i,J) * G%dx_Cv(i,J) |
| 1343 | 8901 | G%IareaCv(i,J) = G%mask2dCv(i,J) * Adcroft_reciprocal(G%areaCv(i,J)) |
| 1344 | enddo ; enddo | |
| 1345 | ||
| 1346 | 1 | call callTree_leave("initialize_masks()") |
| 1347 | 1 | end subroutine initialize_masks |
| 1348 | ||
| 1349 | !> \namespace mom_grid_initialize | |
| 1350 | !! | |
| 1351 | !! The metric terms have the form Dzp, IDzp, or DXDYp, where z can | |
| 1352 | !! be X or Y, and p can be q, u, v, or h. z describes the direction | |
| 1353 | !! of the metric, while p describes the location. IDzp is the | |
| 1354 | !! inverse of Dzp, while DXDYp is the product of DXp and DYp except | |
| 1355 | !! that areaT is calculated analytically from the latitudes and | |
| 1356 | !! longitudes of the surrounding q points. | |
| 1357 | !! | |
| 1358 | !! On a sphere, a variety of grids can be implemented by defining | |
| 1359 | !! analytic expressions for dx_di, dy_dj (where x and y are latitude | |
| 1360 | !! and longitude, and i and j are grid indices) and the expressions | |
| 1361 | !! for the integrals of their inverses in the four subroutines | |
| 1362 | !! dy_dj, Int_dj_dy, dx_di, and Int_di_dx. | |
| 1363 | !! | |
| 1364 | !! initialize_masks sets up land masks based on the depth field. | |
| 1365 | !! The one argument is the minimum ocean depth. Depths that are | |
| 1366 | !! less than this are interpreted as land points. | |
| 1367 | ||
| 1368 | end module MOM_grid_initialize |