MOM6
Data Types | Functions/Subroutines
mom_grid_initialize Module Reference

Detailed Description

Initializes horizontal grid.

The metric terms have the form Dzp, IDzp, or DXDYp, where z can be X or Y, and p can be q, u, v, or h. z describes the direction of the metric, while p describes the location. IDzp is the inverse of Dzp, while DXDYp is the product of DXp and DYp except that areaT is calculated analytically from the latitudes and longitudes of the surrounding q points.

On a sphere, a variety of grids can be implemented by defining analytic expressions for dx_di, dy_dj (where x and y are latitude and longitude, and i and j are grid indices) and the expressions for the integrals of their inverses in the four subroutines dy_dj, Int_dj_dy, dx_di, and Int_di_dx.

initialize_masks sets up land masks based on the depth field. The one argument is the minimum ocean depth. Depths that are less than this are interpreted as land points.

Data Types

type  gps
 Global positioning system (aka container for information to describe the grid) More...
 

Functions/Subroutines

subroutine, public set_grid_metrics (G, param_file, US)
 set_grid_metrics is used to set the primary values in the model's horizontal grid. The bathymetry, land-sea mask and any restricted channel widths are not known yet, so these are set later. More...
 
subroutine grid_metrics_chksum (parent, G, US)
 grid_metrics_chksum performs a set of checksums on metrics on the grid for debugging. More...
 
subroutine set_grid_metrics_from_mosaic (G, param_file, US)
 Sets the grid metrics from a mosaic file. More...
 
subroutine set_grid_metrics_cartesian (G, param_file, US)
 Calculate the values of the metric terms for a Cartesian grid that might be used and save them in arrays. More...
 
subroutine set_grid_metrics_spherical (G, param_file, US)
 Calculate the values of the metric terms that might be used and save them in arrays. More...
 
subroutine set_grid_metrics_mercator (G, param_file, US)
 Calculate the values of the metric terms that might be used and save them in arrays. More...
 
real function ds_di (x, y, GP)
 This function returns the grid spacing in the logical x direction. More...
 
real function ds_dj (x, y, GP)
 This function returns the grid spacing in the logical y direction. More...
 
real function dl (x1, x2, y1, y2)
 This function returns the contribution from the line integral along one of the four sides of a cell face to the area of a cell, assuming that the sides follow a linear path in latitude and longitude (i.e., on a Mercator grid). More...
 
real function find_root (fn, dy_df, GP, fnval, y1, ymin, ymax, ittmax)
 This subroutine finds and returns the value of y at which the monotonically increasing function fn takes the value fnval, also returning in ittmax the number of iterations of Newton's method that were used to polish the root. More...
 
real function dx_di (x, GP)
 This function calculates and returns the value of dx/di, where x is the longitude in Radians, and i is the integral north-south grid index. More...
 
real function int_di_dx (x, GP)
 This function calculates and returns the integral of the inverse of dx/di to the point x, in radians. More...
 
real function dy_dj (y, GP)
 This subroutine calculates and returns the value of dy/dj, where y is the latitude in Radians, and j is the integral north-south grid index. More...
 
real function int_dj_dy (y, GP)
 This subroutine calculates and returns the integral of the inverse of dy/dj to the point y, in radians. More...
 
subroutine extrapolate_metric (var, jh, missing)
 Extrapolates missing metric data into all the halo regions. More...
 
real function, public adcroft_reciprocal (val)
 This function implements Adcroft's rule for reciprocals, namely that Adcroft_Inv(x) = 1/x for |x|>0 or 0 for x=0. More...
 
subroutine, public initialize_masks (G, PF, US)
 Initializes the grid masks and any metrics that come with masks already applied. More...
 

Function/Subroutine Documentation

◆ adcroft_reciprocal()

real function, public mom_grid_initialize::adcroft_reciprocal ( real, intent(in)  val)

This function implements Adcroft's rule for reciprocals, namely that Adcroft_Inv(x) = 1/x for |x|>0 or 0 for x=0.

Parameters
[in]valThe value being inverted.
Returns
The Adcroft reciprocal of val.

Definition at line 1220 of file MOM_grid_initialize.F90.

1221  real, intent(in) :: val !< The value being inverted.
1222  real :: I_val !< The Adcroft reciprocal of val.
1223 
1224  i_val = 0.0
1225  if (val /= 0.0) i_val = 1.0/val

◆ dl()

real function mom_grid_initialize::dl ( real, intent(in)  x1,
real, intent(in)  x2,
real, intent(in)  y1,
real, intent(in)  y2 
)
private

This function returns the contribution from the line integral along one of the four sides of a cell face to the area of a cell, assuming that the sides follow a linear path in latitude and longitude (i.e., on a Mercator grid).

Parameters
[in]x1Segment starting longitude, in degrees E.
[in]x2Segment ending longitude, in degrees E.
[in]y1Segment ending latitude, in degrees N.
[in]y2Segment ending latitude, in degrees N.

Definition at line 957 of file MOM_grid_initialize.F90.

958  real, intent(in) :: x1 !< Segment starting longitude, in degrees E.
959  real, intent(in) :: x2 !< Segment ending longitude, in degrees E.
960  real, intent(in) :: y1 !< Segment ending latitude, in degrees N.
961  real, intent(in) :: y2 !< Segment ending latitude, in degrees N.
962  ! Local variables
963  real :: dL
964  real :: r, dy
965 
966  dy = y2 - y1
967 
968  if (abs(dy) > 2.5e-8) then
969  r = ((1.0 - cos(dy))*cos(y1) + sin(dy)*sin(y1)) / dy
970  else
971  r = (0.5*dy*cos(y1) + sin(y1))
972  endif
973  dl = r * (x2 - x1)
974 

◆ ds_di()

real function mom_grid_initialize::ds_di ( real, intent(in)  x,
real, intent(in)  y,
type(gps), intent(in)  GP 
)
private

This function returns the grid spacing in the logical x direction.

Parameters
[in]xThe longitude in question
[in]yThe latitude in question
[in]gpA structure of grid parameters

Definition at line 927 of file MOM_grid_initialize.F90.

928  real, intent(in) :: x !< The longitude in question
929  real, intent(in) :: y !< The latitude in question
930  type(GPS), intent(in) :: GP !< A structure of grid parameters
931  real :: ds_di
932  ! Local variables
933 
934  ds_di = gp%Rad_Earth * cos(y) * dx_di(x,gp)
935  ! In general, this might be...
936  ! ds_di = GP%Rad_Earth * sqrt( cos(y)*cos(y) * dx_di(x,y,GP)*dx_di(x,y,GP) + &
937  ! dy_di(x,y,GP)*dy_di(x,y,GP))

◆ ds_dj()

real function mom_grid_initialize::ds_dj ( real, intent(in)  x,
real, intent(in)  y,
type(gps), intent(in)  GP 
)
private

This function returns the grid spacing in the logical y direction.

Parameters
[in]xThe longitude in question
[in]yThe latitude in question
[in]gpA structure of grid parameters

Definition at line 941 of file MOM_grid_initialize.F90.

942  real, intent(in) :: x !< The longitude in question
943  real, intent(in) :: y !< The latitude in question
944  type(GPS), intent(in) :: GP !< A structure of grid parameters
945  ! Local variables
946  real :: ds_dj
947 
948  ds_dj = gp%Rad_Earth * dy_dj(y,gp)
949  ! In general, this might be...
950  ! ds_dj = GP%Rad_Earth * sqrt( cos(y)*cos(y) * dx_dj(x,y,GP)*dx_dj(x,y,GP) + &
951  ! dy_dj(x,y,GP)*dy_dj(x,y,GP))

◆ dx_di()

real function mom_grid_initialize::dx_di ( real, intent(in)  x,
type(gps), intent(in)  GP 
)
private

This function calculates and returns the value of dx/di, where x is the longitude in Radians, and i is the integral north-south grid index.

Parameters
[in]xThe longitude in question
[in]gpA structure of grid parameters

Definition at line 1092 of file MOM_grid_initialize.F90.

1093  real, intent(in) :: x !< The longitude in question
1094  type(GPS), intent(in) :: GP !< A structure of grid parameters
1095  real :: dx_di
1096 
1097  dx_di = (gp%len_lon * 4.0*atan(1.0)) / (180.0 * gp%niglobal)
1098 

◆ dy_dj()

real function mom_grid_initialize::dy_dj ( real, intent(in)  y,
type(gps), intent(in)  GP 
)
private

This subroutine calculates and returns the value of dy/dj, where y is the latitude in Radians, and j is the integral north-south grid index.

Parameters
[in]yThe latitude in question
[in]gpA structure of grid parameters

Definition at line 1114 of file MOM_grid_initialize.F90.

1115  real, intent(in) :: y !< The latitude in question
1116  type(GPS), intent(in) :: GP !< A structure of grid parameters
1117  real :: dy_dj
1118  ! Local variables
1119  real :: PI ! 3.1415926... calculated as 4*atan(1)
1120  real :: C0 ! The constant that converts the nominal y-spacing in
1121  ! gridpoints to the nominal spacing in Radians.
1122  real :: y_eq_enhance ! The latitude in radians within which the resolution
1123  ! is enhanced.
1124  pi = 4.0*atan(1.0)
1125  if (gp%isotropic) then
1126  c0 = (gp%len_lon * pi) / (180.0 * gp%niglobal)
1127  y_eq_enhance = pi*abs(gp%lat_eq_enhance)/180.0
1128  if (abs(y) < y_eq_enhance) then
1129  dy_dj = c0 * (cos(y) / (1.0 + 0.5*cos(y) * (gp%lat_enhance_factor - 1.0) * &
1130  (1.0+cos(pi*y/y_eq_enhance)) ))
1131  else
1132  dy_dj = c0 * cos(y)
1133  endif
1134  else
1135  c0 = (gp%len_lat * pi) / (180.0 * gp%njglobal)
1136  dy_dj = c0
1137  endif
1138 

◆ extrapolate_metric()

subroutine mom_grid_initialize::extrapolate_metric ( real, dimension(:,:), intent(inout)  var,
integer, intent(in)  jh,
real, intent(in), optional  missing 
)
private

Extrapolates missing metric data into all the halo regions.

Parameters
[in,out]varThe array in which to fill in halos
[in]jhThe size of the halos to be filled
[in]missingThe missing data fill value, 0 by default.

Definition at line 1186 of file MOM_grid_initialize.F90.

1187  real, dimension(:,:), intent(inout) :: var !< The array in which to fill in halos
1188  integer, intent(in) :: jh !< The size of the halos to be filled
1189  real, optional, intent(in) :: missing !< The missing data fill value, 0 by default.
1190  ! Local variables
1191  real :: badval
1192  integer :: i,j
1193 
1194  badval = 0.0 ; if (present(missing)) badval = missing
1195 
1196  ! Fill in southern halo by extrapolating from the computational domain
1197  do j=lbound(var,2)+jh,lbound(var,2),-1 ; do i=lbound(var,1),ubound(var,1)
1198  if (var(i,j)==badval) var(i,j) = 2.0*var(i,j+1)-var(i,j+2)
1199  enddo ; enddo
1200 
1201  ! Fill in northern halo by extrapolating from the computational domain
1202  do j=ubound(var,2)-jh,ubound(var,2) ; do i=lbound(var,1),ubound(var,1)
1203  if (var(i,j)==badval) var(i,j) = 2.0*var(i,j-1)-var(i,j-2)
1204  enddo ; enddo
1205 
1206  ! Fill in western halo by extrapolating from the computational domain
1207  do j=lbound(var,2),ubound(var,2) ; do i=lbound(var,1)+jh,lbound(var,1),-1
1208  if (var(i,j)==badval) var(i,j) = 2.0*var(i+1,j)-var(i+2,j)
1209  enddo ; enddo
1210 
1211  ! Fill in eastern halo by extrapolating from the computational domain
1212  do j=lbound(var,2),ubound(var,2) ; do i=ubound(var,1)-jh,ubound(var,1)
1213  if (var(i,j)==badval) var(i,j) = 2.0*var(i-1,j)-var(i-2,j)
1214  enddo ; enddo
1215 

◆ find_root()

real function mom_grid_initialize::find_root ( real, external  fn,
real, external  dy_df,
type(gps), intent(in)  GP,
real, intent(in)  fnval,
real, intent(in)  y1,
real, intent(in)  ymin,
real, intent(in)  ymax,
integer, intent(out)  ittmax 
)
private

This subroutine finds and returns the value of y at which the monotonically increasing function fn takes the value fnval, also returning in ittmax the number of iterations of Newton's method that were used to polish the root.

Returns
The value of y where fn(y) = fnval that will be returned
Parameters
fnThe external function whose root is being sought
dy_dfThe inverse of the derivative of that function
[in]gpA structure of grid parameters
[in]fnvalThe value of fn being sought
[in]y1A first guess for y
[in]yminThe minimum permitted value of y
[in]ymaxThe maximum permitted value of y
[out]ittmaxThe number of iterations used to polish the root

Definition at line 980 of file MOM_grid_initialize.F90.

981  real :: find_root !< The value of y where fn(y) = fnval that will be returned
982  real, external :: fn !< The external function whose root is being sought
983  real, external :: dy_df !< The inverse of the derivative of that function
984  type(GPS), intent(in) :: GP !< A structure of grid parameters
985  real, intent(in) :: fnval !< The value of fn being sought
986  real, intent(in) :: y1 !< A first guess for y
987  real, intent(in) :: ymin !< The minimum permitted value of y
988  real, intent(in) :: ymax !< The maximum permitted value of y
989  integer, intent(out) :: ittmax !< The number of iterations used to polish the root
990  ! Local variables
991  real :: y, y_next
992  real :: ybot, ytop, fnbot, fntop
993  integer :: itt
994  character(len=256) :: warnmesg
995 
996  real :: dy_dfn, dy, fny
997 
998 ! Bracket the root. Do not use the bounding values because the value at the
999 ! function at the bounds could be infinite, as is the case for the Mercator
1000 ! grid recursion relation. (I.e., this is a search on an open interval.)
1001  ybot = y1
1002  fnbot = fn(ybot,gp) - fnval ; itt = 0
1003  do while (fnbot > 0.0)
1004  if ((ybot - 2.0*dy_df(ybot,gp)) < (0.5*(ybot+ymin))) then
1005  ! Go twice as far as the secant method would normally go.
1006  ybot = ybot - 2.0*dy_df(ybot,gp)
1007  else ! But stay within the open interval!
1008  ybot = 0.5*(ybot+ymin) ; itt = itt + 1
1009  endif
1010  fnbot = fn(ybot,gp) - fnval
1011 
1012  if ((itt > 50) .and. (fnbot > 0.0)) then
1013  write(warnmesg, '("PE ",I2," unable to find bottom bound for grid function. &
1014  &x = ",ES10.4,", xmax = ",ES10.4,", fn = ",ES10.4,", dfn_dx = ",ES10.4,&
1015  &", seeking fn = ",ES10.4," - fn = ",ES10.4,".")') &
1016  pe_here(),ybot,ymin,fn(ybot,gp),dy_df(ybot,gp),fnval, fnbot
1017  call mom_error(fatal,warnmesg)
1018  endif
1019  enddo
1020 
1021  ytop = y1
1022  fntop = fn(ytop,gp) - fnval ; itt = 0
1023  do while (fntop < 0.0)
1024  if ((ytop + 2.0*dy_df(ytop,gp)) < (0.5*(ytop+ymax))) then
1025  ! Go twice as far as the secant method would normally go.
1026  ytop = ytop + 2.0*dy_df(ytop,gp)
1027  else ! But stay within the open interval!
1028  ytop = 0.5*(ytop+ymax) ; itt = itt + 1
1029  endif
1030  fntop = fn(ytop,gp) - fnval
1031 
1032  if ((itt > 50) .and. (fntop < 0.0)) then
1033  write(warnmesg, '("PE ",I2," unable to find top bound for grid function. &
1034  &x = ",ES10.4,", xmax = ",ES10.4,", fn = ",ES10.4,", dfn_dx = ",ES10.4, &
1035  &", seeking fn = ",ES10.4," - fn = ",ES10.4,".")') &
1036  pe_here(),ytop,ymax,fn(ytop,gp),dy_df(ytop,gp),fnval,fntop
1037  call mom_error(fatal,warnmesg)
1038  endif
1039  enddo
1040 
1041  ! Find the root using a bracketed variant of Newton's method, starting
1042  ! with a false-positon method first guess.
1043  if ((fntop < 0.0) .or. (fnbot > 0.0) .or. (ytop < ybot)) then
1044  write(warnmesg, '("PE ",I2," find_root failed to bracket function. y = ",&
1045  &2ES10.4,", fn = ",2ES10.4,".")') pe_here(),ybot,ytop,fnbot,fntop
1046  call mom_error(fatal, warnmesg)
1047  endif
1048 
1049  if (fntop == 0.0) then ; y = ytop ; fny = fntop
1050  elseif (fnbot == 0.0) then ; y = ybot ; fny = fnbot
1051  else
1052  y = (ybot*fntop - ytop*fnbot) / (fntop - fnbot)
1053  fny = fn(y,gp) - fnval
1054  if (fny < 0.0) then ; fnbot = fny ; ybot = y
1055  else ; fntop = fny ; ytop = y ; endif
1056  endif
1057 
1058  do itt=1,50
1059  dy_dfn = dy_df(y,gp)
1060 
1061  dy = -1.0* fny * dy_dfn
1062  y_next = y + dy
1063  if ((y_next >= ytop) .or. (y_next <= ybot)) then
1064  ! The Newton's method estimate has escaped bracketing, so use the
1065  ! false-position method instead. The complicated test is to properly
1066  ! handle the case where the iteration is down to roundoff level differences.
1067  y_next = y
1068  if (abs(fntop - fnbot) > epsilon(y) * (abs(fntop) + abs(fnbot))) &
1069  y_next = (ybot*fntop - ytop*fnbot) / (fntop - fnbot)
1070  endif
1071 
1072  dy = y_next - y
1073  if (abs(dy) < (2.0*epsilon(y)*(abs(y) + abs(y_next)) + 1.0e-20)) then
1074  y = y_next ; exit
1075  endif
1076  y = y_next
1077 
1078  fny = fn(y,gp) - fnval
1079  if (fny > 0.0) then ; ytop = y ; fntop = fny
1080  elseif (fny < 0.0) then ; ybot = y ; fnbot = fny
1081  else ; exit ; endif
1082 
1083  enddo
1084  if (abs(y) < 1e-12) y = 0.0
1085 
1086  ittmax = itt
1087  find_root = y

◆ grid_metrics_chksum()

subroutine mom_grid_initialize::grid_metrics_chksum ( character(len=*), intent(in)  parent,
type(dyn_horgrid_type), intent(in)  G,
type(unit_scale_type), intent(in), optional  US 
)
private

grid_metrics_chksum performs a set of checksums on metrics on the grid for debugging.

Parameters
[in]parentA string identifying the caller
[in]gThe dynamic horizontal grid type
[in]usA dimensional unit scaling type

Definition at line 115 of file MOM_grid_initialize.F90.

116  character(len=*), intent(in) :: parent !< A string identifying the caller
117  type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type
118  type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type
119 
120  real :: m_to_L ! A unit conversion factor [L m-1 ~> nondim]
121  real :: L_to_m ! A unit conversion factor [m L-1 ~> nondim]
122  integer :: halo
123  m_to_l = 1.0 ; if (present(us)) m_to_l = us%m_to_L
124  l_to_m = 1.0 ; if (present(us)) l_to_m = us%L_to_m
125 
126  halo = min(g%ied-g%iec, g%jed-g%jec, 1)
127 
128  call hchksum_pair(trim(parent)//': d[xy]T', g%dxT, g%dyT, g%HI, &
129  haloshift=halo, scale=l_to_m, scalar_pair=.true.)
130 
131  call uvchksum(trim(parent)//': dxC[uv]', g%dxCu, g%dyCv, g%HI, haloshift=halo, scale=l_to_m)
132 
133  call uvchksum(trim(parent)//': dxC[uv]', g%dyCu, g%dxCv, g%HI, haloshift=halo, scale=l_to_m)
134 
135  call bchksum_pair(trim(parent)//': dxB[uv]', g%dxBu, g%dyBu, g%HI, haloshift=halo, scale=l_to_m)
136 
137  call hchksum_pair(trim(parent)//': Id[xy]T', g%IdxT, g%IdyT, g%HI, &
138  haloshift=halo, scale=m_to_l, scalar_pair=.true.)
139 
140  call uvchksum(trim(parent)//': Id[xy]C[uv]', g%IdxCu, g%IdyCv, g%HI, haloshift=halo, scale=m_to_l)
141 
142  call uvchksum(trim(parent)//': Id[xy]C[uv]', g%IdyCu, g%IdxCv, g%HI, haloshift=halo, scale=m_to_l)
143 
144  call bchksum_pair(trim(parent)//': Id[xy]B[uv]', g%IdxBu, g%IdyBu, g%HI, haloshift=halo, scale=m_to_l)
145 
146  call hchksum(g%areaT, trim(parent)//': areaT',g%HI, haloshift=halo, scale=l_to_m**2)
147  call bchksum(g%areaBu, trim(parent)//': areaBu',g%HI, haloshift=halo, scale=l_to_m**2)
148 
149  call hchksum(g%IareaT, trim(parent)//': IareaT',g%HI, haloshift=halo, scale=m_to_l**2)
150  call bchksum(g%IareaBu, trim(parent)//': IareaBu',g%HI, haloshift=halo, scale=m_to_l**2)
151 
152  call hchksum(g%geoLonT,trim(parent)//': geoLonT',g%HI, haloshift=halo)
153  call hchksum(g%geoLatT,trim(parent)//': geoLatT',g%HI, haloshift=halo)
154 
155  call bchksum(g%geoLonBu, trim(parent)//': geoLonBu',g%HI, haloshift=halo)
156  call bchksum(g%geoLatBu, trim(parent)//': geoLatBu',g%HI, haloshift=halo)
157 
158  call uvchksum(trim(parent)//': geoLonC[uv]', g%geoLonCu, g%geoLonCv, g%HI, haloshift=halo)
159 
160  call uvchksum(trim(parent)//': geoLatC[uv]', g%geoLatCu, g%geoLatCv, g%HI, haloshift=halo)
161 

◆ initialize_masks()

subroutine, public mom_grid_initialize::initialize_masks ( type(dyn_horgrid_type), intent(inout)  G,
type(param_file_type), intent(in)  PF,
type(unit_scale_type), intent(in), optional  US 
)

Initializes the grid masks and any metrics that come with masks already applied.

Initialize_masks sets mask2dT, mask2dCu, mask2dCv, and mask2dBu to mask out flow over any points which are shallower than Dmin and permit an appropriate treatment of the boundary conditions. mask2dCu and mask2dCv are 0.0 at any points adjacent to a land point. mask2dBu is 0.0 at any land or boundary point. For points in the interior, mask2dCu, mask2dCv, and mask2dBu are all 1.0.

Parameters
[in,out]gThe dynamic horizontal grid type
[in]pfParameter file structure
[in]usA dimensional unit scaling type

Definition at line 1236 of file MOM_grid_initialize.F90.

1237  type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type
1238  type(param_file_type), intent(in) :: PF !< Parameter file structure
1239  type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type
1240  ! Local variables
1241  real :: m_to_Z_scale ! A unit conversion factor from m to Z.
1242  real :: m_to_L ! A unit conversion factor [L m-1 ~> nondim]
1243  real :: Dmin ! The depth for masking in the same units as G%bathyT [Z ~> m].
1244  real :: min_depth ! The minimum ocean depth in the same units as G%bathyT [Z ~> m].
1245  real :: mask_depth ! The depth shallower than which to mask a point as land [Z ~> m].
1246  character(len=40) :: mdl = "MOM_grid_init initialize_masks"
1247  integer :: i, j
1248 
1249  call calltree_enter("initialize_masks(), MOM_grid_initialize.F90")
1250  m_to_z_scale = 1.0 ; if (present(us)) m_to_z_scale = us%m_to_Z
1251  m_to_l = 1.0 ; if (present(us)) m_to_l = us%m_to_L
1252 
1253  call get_param(pf, mdl, "MINIMUM_DEPTH", min_depth, &
1254  "If MASKING_DEPTH is unspecified, then anything shallower than "//&
1255  "MINIMUM_DEPTH is assumed to be land and all fluxes are masked out. "//&
1256  "If MASKING_DEPTH is specified, then all depths shallower than "//&
1257  "MINIMUM_DEPTH but deeper than MASKING_DEPTH are rounded to MINIMUM_DEPTH.", &
1258  units="m", default=0.0, scale=m_to_z_scale)
1259  call get_param(pf, mdl, "MASKING_DEPTH", mask_depth, &
1260  "The depth below which to mask points as land points, for which all "//&
1261  "fluxes are zeroed out. MASKING_DEPTH is ignored if negative.", &
1262  units="m", default=-9999.0, scale=m_to_z_scale)
1263 
1264  dmin = min_depth
1265  if (mask_depth>=0.) dmin = mask_depth
1266 
1267  g%mask2dCu(:,:) = 0.0 ; g%mask2dCv(:,:) = 0.0 ; g%mask2dBu(:,:) = 0.0
1268 
1269  ! Construct the h-point or T-point mask
1270  do j=g%jsd,g%jed ; do i=g%isd,g%ied
1271  if (g%bathyT(i,j) <= dmin) then
1272  g%mask2dT(i,j) = 0.0
1273  else
1274  g%mask2dT(i,j) = 1.0
1275  endif
1276  enddo ; enddo
1277 
1278  do j=g%jsd,g%jed ; do i=g%isd,g%ied-1
1279  if ((g%bathyT(i,j) <= dmin) .or. (g%bathyT(i+1,j) <= dmin)) then
1280  g%mask2dCu(i,j) = 0.0
1281  else
1282  g%mask2dCu(i,j) = 1.0
1283  endif
1284  enddo ; enddo
1285 
1286  do j=g%jsd,g%jed-1 ; do i=g%isd,g%ied
1287  if ((g%bathyT(i,j) <= dmin) .or. (g%bathyT(i,j+1) <= dmin)) then
1288  g%mask2dCv(i,j) = 0.0
1289  else
1290  g%mask2dCv(i,j) = 1.0
1291  endif
1292  enddo ; enddo
1293 
1294  do j=g%jsd,g%jed-1 ; do i=g%isd,g%ied-1
1295  if ((g%bathyT(i+1,j) <= dmin) .or. (g%bathyT(i+1,j+1) <= dmin) .or. &
1296  (g%bathyT(i,j) <= dmin) .or. (g%bathyT(i,j+1) <= dmin)) then
1297  g%mask2dBu(i,j) = 0.0
1298  else
1299  g%mask2dBu(i,j) = 1.0
1300  endif
1301  enddo ; enddo
1302 
1303  call pass_var(g%mask2dBu, g%Domain, position=corner)
1304  call pass_vector(g%mask2dCu, g%mask2dCv, g%Domain, to_all+scalar_pair, cgrid_ne)
1305 
1306  do j=g%jsd,g%jed ; do i=g%IsdB,g%IedB
1307  g%dy_Cu(i,j) = g%mask2dCu(i,j) * g%dyCu(i,j)
1308  g%areaCu(i,j) = g%dxCu(i,j) * g%dy_Cu(i,j)
1309  g%IareaCu(i,j) = g%mask2dCu(i,j) * adcroft_reciprocal(g%areaCu(i,j))
1310  enddo ; enddo
1311 
1312  do j=g%JsdB,g%JedB ; do i=g%isd,g%ied
1313  g%dx_Cv(i,j) = g%mask2dCv(i,j) * g%dxCv(i,j)
1314  g%areaCv(i,j) = g%dyCv(i,j) * g%dx_Cv(i,j)
1315  g%IareaCv(i,j) = g%mask2dCv(i,j) * adcroft_reciprocal(g%areaCv(i,j))
1316  enddo ; enddo
1317 
1318  call calltree_leave("initialize_masks()")

◆ int_di_dx()

real function mom_grid_initialize::int_di_dx ( real, intent(in)  x,
type(gps), intent(in)  GP 
)
private

This function calculates and returns the integral of the inverse of dx/di to the point x, in radians.

Parameters
[in]xThe longitude in question
[in]gpA structure of grid parameters

Definition at line 1103 of file MOM_grid_initialize.F90.

1104  real, intent(in) :: x !< The longitude in question
1105  type(GPS), intent(in) :: GP !< A structure of grid parameters
1106  real :: Int_di_dx
1107 
1108  int_di_dx = x * ((180.0 * gp%niglobal) / (gp%len_lon * 4.0*atan(1.0)))
1109 

◆ int_dj_dy()

real function mom_grid_initialize::int_dj_dy ( real, intent(in)  y,
type(gps), intent(in)  GP 
)
private

This subroutine calculates and returns the integral of the inverse of dy/dj to the point y, in radians.

Parameters
[in]yThe latitude in question
[in]gpA structure of grid parameters

Definition at line 1143 of file MOM_grid_initialize.F90.

1144  real, intent(in) :: y !< The latitude in question
1145  type(GPS), intent(in) :: GP !< A structure of grid parameters
1146  real :: Int_dj_dy
1147  ! Local variables
1148  real :: I_C0 = 0.0 ! The inverse of the constant that converts the
1149  ! nominal spacing in gridpoints to the nominal
1150  ! spacing in Radians.
1151  real :: PI ! 3.1415926... calculated as 4*atan(1)
1152  real :: y_eq_enhance ! The latitude in radians from
1153  ! from the equator within which the
1154  ! meridional grid spacing is enhanced by
1155  ! a factor of GP%lat_enhance_factor.
1156  real :: r
1157 
1158  pi = 4.0*atan(1.0)
1159  if (gp%isotropic) then
1160  i_c0 = (180.0 * gp%niglobal) / (gp%len_lon * pi)
1161  y_eq_enhance = pi*abs(gp%lat_eq_enhance)/180.0
1162 
1163  if (y >= 0.0) then
1164  r = i_c0 * log((1.0 + sin(y))/cos(y))
1165  else
1166  r = -1.0 * i_c0 * log((1.0 - sin(y))/cos(y))
1167  endif
1168 
1169  if (y >= y_eq_enhance) then
1170  r = r + i_c0*0.5*(gp%lat_enhance_factor - 1.0)*y_eq_enhance
1171  elseif (y <= -y_eq_enhance) then
1172  r = r - i_c0*0.5*(gp%lat_enhance_factor - 1.0)*y_eq_enhance
1173  else
1174  r = r + i_c0*0.5*(gp%lat_enhance_factor - 1.0) * &
1175  (y + (y_eq_enhance/pi)*sin(pi*y/y_eq_enhance))
1176  endif
1177  else
1178  i_c0 = (180.0 * gp%njglobal) / (gp%len_lat * pi)
1179  r = i_c0 * y
1180  endif
1181 
1182  int_dj_dy = r

◆ set_grid_metrics()

subroutine, public mom_grid_initialize::set_grid_metrics ( type(dyn_horgrid_type), intent(inout)  G,
type(param_file_type), intent(in)  param_file,
type(unit_scale_type), intent(in), optional  US 
)

set_grid_metrics is used to set the primary values in the model's horizontal grid. The bathymetry, land-sea mask and any restricted channel widths are not known yet, so these are set later.

Parameters
[in,out]gThe dynamic horizontal grid type
[in]param_fileParameter file structure
[in]usA dimensional unit scaling type

Definition at line 62 of file MOM_grid_initialize.F90.

63  type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type
64  type(param_file_type), intent(in) :: param_file !< Parameter file structure
65  type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type
66  ! Local variables
67 ! This include declares and sets the variable "version".
68 #include "version_variable.h"
69  logical :: debug
70  character(len=256) :: config
71 
72  call calltree_enter("set_grid_metrics(), MOM_grid_initialize.F90")
73  call log_version(param_file, "MOM_grid_init", version, "")
74  call get_param(param_file, "MOM_grid_init", "GRID_CONFIG", config, &
75  "A character string that determines the method for "//&
76  "defining the horizontal grid. Current options are: \n"//&
77  " \t mosaic - read the grid from a mosaic (supergrid) \n"//&
78  " \t file set by GRID_FILE.\n"//&
79  " \t cartesian - use a (flat) Cartesian grid.\n"//&
80  " \t spherical - use a simple spherical grid.\n"//&
81  " \t mercator - use a Mercator spherical grid.", &
82  fail_if_missing=.true.)
83  call get_param(param_file, "MOM_grid_init", "DEBUG", debug, &
84  "If true, write out verbose debugging data.", &
85  default=.false., debuggingparam=.true.)
86 
87  ! These are defaults that may be changed in the next select block.
88  g%x_axis_units = "degrees_east" ; g%y_axis_units = "degrees_north"
89  select case (trim(config))
90  case ("mosaic"); call set_grid_metrics_from_mosaic(g, param_file, us)
91  case ("cartesian"); call set_grid_metrics_cartesian(g, param_file, us)
92  case ("spherical"); call set_grid_metrics_spherical(g, param_file, us)
93  case ("mercator"); call set_grid_metrics_mercator(g, param_file, us)
94  case ("file"); call mom_error(fatal, "MOM_grid_init: set_grid_metrics "//&
95  'GRID_CONFIG "file" is no longer a supported option. Use a '//&
96  'mosaic file ("mosaic") or one of the analytic forms instead.')
97  case default ; call mom_error(fatal, "MOM_grid_init: set_grid_metrics "//&
98  "Unrecognized grid configuration "//trim(config))
99  end select
100 
101  ! Calculate derived metrics (i.e. reciprocals and products)
102  call calltree_enter("set_derived_metrics(), MOM_grid_initialize.F90")
103  call set_derived_dyn_horgrid(g, us)
104  call calltree_leave("set_derived_metrics()")
105 
106  if (debug) call grid_metrics_chksum('MOM_grid_init/set_grid_metrics', g, us)
107 
108  call calltree_leave("set_grid_metrics()")

◆ set_grid_metrics_cartesian()

subroutine mom_grid_initialize::set_grid_metrics_cartesian ( type(dyn_horgrid_type), intent(inout)  G,
type(param_file_type), intent(in)  param_file,
type(unit_scale_type), intent(in), optional  US 
)
private

Calculate the values of the metric terms for a Cartesian grid that might be used and save them in arrays.

Within this subroutine, the x- and y- grid spacings and their inverses and the cell areas centered on h, q, u, and v points are calculated, as are the geographic locations of each of these 4 sets of points.

Parameters
[in,out]gThe dynamic horizontal grid type
[in]param_fileParameter file structure
[in]usA dimensional unit scaling type

Definition at line 419 of file MOM_grid_initialize.F90.

420  type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type
421  type(param_file_type), intent(in) :: param_file !< Parameter file structure
422  type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type
423  ! Local variables
424  integer :: i, j, isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, I1off, J1off
425  integer :: niglobal, njglobal
426  real :: grid_latT(G%jsd:G%jed), grid_latB(G%JsdB:G%JedB)
427  real :: grid_lonT(G%isd:G%ied), grid_lonB(G%IsdB:G%IedB)
428  real :: dx_everywhere, dy_everywhere ! Grid spacings [m].
429  real :: I_dx, I_dy ! Inverse grid spacings [m-1].
430  real :: PI
431  real :: m_to_L ! A unit conversion factor [L m-1 ~> nondim]
432  real :: L_to_m ! A unit conversion factor [m L-1 ~> nondim]
433  character(len=80) :: units_temp
434  character(len=48) :: mdl = "MOM_grid_init set_grid_metrics_cartesian"
435 
436  niglobal = g%Domain%niglobal ; njglobal = g%Domain%njglobal
437  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
438  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
439  i1off = g%idg_offset ; j1off = g%jdg_offset
440 
441  call calltree_enter("set_grid_metrics_cartesian(), MOM_grid_initialize.F90")
442 
443  m_to_l = 1.0 ; if (present(us)) m_to_l = us%m_to_L
444  l_to_m = 1.0 ; if (present(us)) l_to_m = us%L_to_m
445  pi = 4.0*atan(1.0)
446 
447  call get_param(param_file, mdl, "AXIS_UNITS", units_temp, &
448  "The units for the Cartesian axes. Valid entries are: \n"//&
449  " \t degrees - degrees of latitude and longitude \n"//&
450  " \t m - meters \n \t k - kilometers", default="degrees")
451  call get_param(param_file, mdl, "SOUTHLAT", g%south_lat, &
452  "The southern latitude of the domain or the equivalent "//&
453  "starting value for the y-axis.", units=units_temp, &
454  fail_if_missing=.true.)
455  call get_param(param_file, mdl, "LENLAT", g%len_lat, &
456  "The latitudinal or y-direction length of the domain.", &
457  units=units_temp, fail_if_missing=.true.)
458  call get_param(param_file, mdl, "WESTLON", g%west_lon, &
459  "The western longitude of the domain or the equivalent "//&
460  "starting value for the x-axis.", units=units_temp, &
461  default=0.0)
462  call get_param(param_file, mdl, "LENLON", g%len_lon, &
463  "The longitudinal or x-direction length of the domain.", &
464  units=units_temp, fail_if_missing=.true.)
465  call get_param(param_file, mdl, "RAD_EARTH", g%Rad_Earth, &
466  "The radius of the Earth.", units="m", default=6.378e6)
467 
468  if (units_temp(1:1) == 'k') then
469  g%x_axis_units = "kilometers" ; g%y_axis_units = "kilometers"
470  elseif (units_temp(1:1) == 'm') then
471  g%x_axis_units = "meters" ; g%y_axis_units = "meters"
472  endif
473  call log_param(param_file, mdl, "explicit AXIS_UNITS", g%x_axis_units)
474 
475  ! Note that the dynamic grid always uses symmetric memory for the global
476  ! arrays G%gridLatB and G%gridLonB.
477  do j=g%jsg-1,g%jeg
478  g%gridLatB(j) = g%south_lat + g%len_lat*real(j-(g%jsg-1))/real(njglobal)
479  enddo
480  do j=g%jsg,g%jeg
481  g%gridLatT(j) = g%south_lat + g%len_lat*(real(j-g%jsg)+0.5)/real(njglobal)
482  enddo
483  do i=g%isg-1,g%ieg
484  g%gridLonB(i) = g%west_lon + g%len_lon*real(i-(g%isg-1))/real(niglobal)
485  enddo
486  do i=g%isg,g%ieg
487  g%gridLonT(i) = g%west_lon + g%len_lon*(real(i-g%isg)+0.5)/real(niglobal)
488  enddo
489 
490  do j=jsdb,jedb
491  grid_latb(j) = g%south_lat + g%len_lat*real(j+j1off-(g%jsg-1))/real(njglobal)
492  enddo
493  do j=jsd,jed
494  grid_latt(j) = g%south_lat + g%len_lat*(real(j+j1off-g%jsg)+0.5)/real(njglobal)
495  enddo
496  do i=isdb,iedb
497  grid_lonb(i) = g%west_lon + g%len_lon*real(i+i1off-(g%isg-1))/real(niglobal)
498  enddo
499  do i=isd,ied
500  grid_lont(i) = g%west_lon + g%len_lon*(real(i+i1off-g%isg)+0.5)/real(niglobal)
501  enddo
502 
503  if (units_temp(1:1) == 'k') then ! Axes are measured in km.
504  dx_everywhere = 1000.0 * g%len_lon / (real(niglobal))
505  dy_everywhere = 1000.0 * g%len_lat / (real(njglobal))
506  elseif (units_temp(1:1) == 'm') then ! Axes are measured in m.
507  dx_everywhere = g%len_lon / (real(niglobal))
508  dy_everywhere = g%len_lat / (real(njglobal))
509  else ! Axes are measured in degrees of latitude and longitude.
510  dx_everywhere = g%Rad_Earth * g%len_lon * pi / (180.0 * niglobal)
511  dy_everywhere = g%Rad_Earth * g%len_lat * pi / (180.0 * njglobal)
512  endif
513 
514  i_dx = 1.0 / dx_everywhere ; i_dy = 1.0 / dy_everywhere
515 
516  do j=jsdb,jedb ; do i=isdb,iedb
517  g%geoLonBu(i,j) = grid_lonb(i) ; g%geoLatBu(i,j) = grid_latb(j)
518 
519  g%dxBu(i,j) = m_to_l*dx_everywhere ; g%IdxBu(i,j) = l_to_m*i_dx
520  g%dyBu(i,j) = m_to_l*dy_everywhere ; g%IdyBu(i,j) = l_to_m*i_dy
521  g%areaBu(i,j) = m_to_l**2*dx_everywhere * dy_everywhere ; g%IareaBu(i,j) = l_to_m**2*i_dx * i_dy
522  enddo ; enddo
523 
524  do j=jsd,jed ; do i=isd,ied
525  g%geoLonT(i,j) = grid_lont(i) ; g%geoLatT(i,j) = grid_latt(j)
526  g%dxT(i,j) = m_to_l*dx_everywhere ; g%IdxT(i,j) = l_to_m*i_dx
527  g%dyT(i,j) = m_to_l*dy_everywhere ; g%IdyT(i,j) = l_to_m*i_dy
528  g%areaT(i,j) = m_to_l**2*dx_everywhere * dy_everywhere ; g%IareaT(i,j) = l_to_m**2*i_dx * i_dy
529  enddo ; enddo
530 
531  do j=jsd,jed ; do i=isdb,iedb
532  g%geoLonCu(i,j) = grid_lonb(i) ; g%geoLatCu(i,j) = grid_latt(j)
533 
534  g%dxCu(i,j) = m_to_l*dx_everywhere ; g%IdxCu(i,j) = l_to_m*i_dx
535  g%dyCu(i,j) = m_to_l*dy_everywhere ; g%IdyCu(i,j) = l_to_m*i_dy
536  enddo ; enddo
537 
538  do j=jsdb,jedb ; do i=isd,ied
539  g%geoLonCv(i,j) = grid_lont(i) ; g%geoLatCv(i,j) = grid_latb(j)
540 
541  g%dxCv(i,j) = m_to_l*dx_everywhere ; g%IdxCv(i,j) = l_to_m*i_dx
542  g%dyCv(i,j) = m_to_l*dy_everywhere ; g%IdyCv(i,j) = l_to_m*i_dy
543  enddo ; enddo
544 
545  call calltree_leave("set_grid_metrics_cartesian()")

◆ set_grid_metrics_from_mosaic()

subroutine mom_grid_initialize::set_grid_metrics_from_mosaic ( type(dyn_horgrid_type), intent(inout)  G,
type(param_file_type), intent(in)  param_file,
type(unit_scale_type), intent(in), optional  US 
)
private

Sets the grid metrics from a mosaic file.

Parameters
[in,out]gThe dynamic horizontal grid type
[in]param_fileParameter file structure
[in]usA dimensional unit scaling type

Definition at line 167 of file MOM_grid_initialize.F90.

168  type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type
169  type(param_file_type), intent(in) :: param_file !< Parameter file structure
170  type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type
171  ! Local variables
172  real, dimension(G%isd :G%ied ,G%jsd :G%jed ) :: tempH1, tempH2, tempH3, tempH4
173  real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB) :: tempQ1, tempQ2, tempQ3, tempQ4
174  real, dimension(G%IsdB:G%IedB,G%jsd :G%jed ) :: tempE1, tempE2
175  real, dimension(G%isd :G%ied ,G%JsdB:G%JedB) :: tempN1, tempN2
176  ! These arrays are a holdover from earlier code in which the arrays in G were
177  ! macros and may have had reduced dimensions.
178  real, dimension(G%isd :G%ied ,G%jsd :G%jed ) :: dxT, dyT, areaT
179  real, dimension(G%IsdB:G%IedB,G%jsd :G%jed ) :: dxCu, dyCu
180  real, dimension(G%isd :G%ied ,G%JsdB:G%JedB) :: dxCv, dyCv
181  real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB) :: dxBu, dyBu, areaBu
182  ! This are symmetric arrays, corresponding to the data in the mosaic file
183  real, dimension(2*G%isd-2:2*G%ied+1,2*G%jsd-2:2*G%jed+1) :: tmpT
184  real, dimension(2*G%isd-3:2*G%ied+1,2*G%jsd-2:2*G%jed+1) :: tmpU
185  real, dimension(2*G%isd-2:2*G%ied+1,2*G%jsd-3:2*G%jed+1) :: tmpV
186  real, dimension(2*G%isd-3:2*G%ied+1,2*G%jsd-3:2*G%jed+1) :: tmpZ
187  real, dimension(:,:), allocatable :: tmpGlbl
188  real :: m_to_L ! A unit conversion factor [L m-1 ~> nondim]
189  character(len=200) :: filename, grid_file, inputdir
190  character(len=64) :: mdl = "MOM_grid_init set_grid_metrics_from_mosaic"
191  integer :: err=0, ni, nj, global_indices(4)
192  type(MOM_domain_type) :: SGdom ! Supergrid domain
193  logical :: lon_bug ! If true use an older buggy answer in the tripolar longitude.
194  integer :: i, j, i2, j2
195  integer :: npei,npej
196  integer, dimension(:), allocatable :: exni,exnj
197  integer :: start(4), nread(4)
198 
199  call calltree_enter("set_grid_metrics_from_mosaic(), MOM_grid_initialize.F90")
200 
201  m_to_l = 1.0 ; if (present(us)) m_to_l = us%m_to_L
202  call get_param(param_file, mdl, "GRID_FILE", grid_file, &
203  "Name of the file from which to read horizontal grid data.", &
204  fail_if_missing=.true.)
205  call get_param(param_file, mdl, "USE_TRIPOLAR_GEOLONB_BUG", lon_bug, &
206  "If true, use older code that incorrectly sets the longitude "//&
207  "in some points along the tripolar fold to be off by 360 degrees.", &
208  default=.false.)
209  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
210  inputdir = slasher(inputdir)
211  filename = trim(adjustl(inputdir)) // trim(adjustl(grid_file))
212  call log_param(param_file, mdl, "INPUTDIR/GRID_FILE", filename)
213  if (.not.file_exists(filename)) &
214  call mom_error(fatal," set_grid_metrics_from_mosaic: Unable to open "//&
215  trim(filename))
216 
217  ! Initialize everything to 0.
218  dxcu(:,:) = 0.0 ; dycu(:,:) = 0.0
219  dxcv(:,:) = 0.0 ; dycv(:,:) = 0.0
220  dxbu(:,:) = 0.0 ; dybu(:,:) = 0.0 ; areabu(:,:) = 0.0
221 
222  !<MISSING CODE TO READ REFINEMENT LEVEL>
223  ni = 2*(g%iec-g%isc+1) ! i size of supergrid
224  nj = 2*(g%jec-g%jsc+1) ! j size of supergrid
225 
226  ! Define a domain for the supergrid (SGdom)
227  npei = g%domain%layout(1) ; npej = g%domain%layout(2)
228  allocate(exni(npei)) ; allocate(exnj(npej))
229  call mpp_get_domain_extents(g%domain%mpp_domain, exni, exnj)
230  allocate(sgdom%mpp_domain)
231  sgdom%nihalo = 2*g%domain%nihalo+1
232  sgdom%njhalo = 2*g%domain%njhalo+1
233  sgdom%niglobal = 2*g%domain%niglobal
234  sgdom%njglobal = 2*g%domain%njglobal
235  sgdom%layout(:) = g%domain%layout(:)
236  sgdom%io_layout(:) = g%domain%io_layout(:)
237  global_indices(1) = 1+sgdom%nihalo
238  global_indices(2) = sgdom%niglobal+sgdom%nihalo
239  global_indices(3) = 1+sgdom%njhalo
240  global_indices(4) = sgdom%njglobal+sgdom%njhalo
241  exni(:) = 2*exni(:) ; exnj(:) = 2*exnj(:)
242  if (associated(g%domain%maskmap)) then
243  call mom_define_domain(global_indices, sgdom%layout, sgdom%mpp_domain, &
244  xflags=g%domain%X_FLAGS, yflags=g%domain%Y_FLAGS, &
245  xhalo=sgdom%nihalo, yhalo=sgdom%njhalo, &
246  xextent=exni,yextent=exnj, &
247  symmetry=.true., name="MOM_MOSAIC", maskmap=g%domain%maskmap)
248  else
249  call mom_define_domain(global_indices, sgdom%layout, sgdom%mpp_domain, &
250  xflags=g%domain%X_FLAGS, yflags=g%domain%Y_FLAGS, &
251  xhalo=sgdom%nihalo, yhalo=sgdom%njhalo, &
252  xextent=exni,yextent=exnj, &
253  symmetry=.true., name="MOM_MOSAIC")
254  endif
255 
256  call mom_define_io_domain(sgdom%mpp_domain, sgdom%io_layout)
257  deallocate(exni)
258  deallocate(exnj)
259 
260  ! Read X from the supergrid
261  tmpz(:,:) = 999.
262  call mom_read_data(filename, 'x', tmpz, sgdom, position=corner)
263 
264  if (lon_bug) then
265  call pass_var(tmpz, sgdom, position=corner)
266  else
267  call pass_var(tmpz, sgdom, position=corner, inner_halo=0)
268  endif
269  call extrapolate_metric(tmpz, 2*(g%jsc-g%jsd)+2, missing=999.)
270  do j=g%jsd,g%jed ; do i=g%isd,g%ied ; i2 = 2*i ; j2 = 2*j
271  g%geoLonT(i,j) = tmpz(i2-1,j2-1)
272  enddo ; enddo
273  do j=g%JsdB,g%JedB ; do i=g%IsdB,g%IedB ; i2 = 2*i ; j2 = 2*j
274  g%geoLonBu(i,j) = tmpz(i2,j2)
275  enddo ; enddo
276  do j=g%jsd,g%jed ; do i=g%IsdB,g%IedB ; i2 = 2*i ; j2 = 2*j
277  g%geoLonCu(i,j) = tmpz(i2,j2-1)
278  enddo ; enddo
279  do j=g%JsdB,g%JedB ; do i=g%isd,g%ied ; i2 = 2*i ; j2 = 2*j
280  g%geoLonCv(i,j) = tmpz(i2-1,j2)
281  enddo ; enddo
282  ! For some reason, this messes up the solution...
283  ! call pass_var(G%geoLonBu, G%domain, position=CORNER)
284 
285  ! Read Y from the supergrid
286  tmpz(:,:) = 999.
287  call mom_read_data(filename, 'y', tmpz, sgdom, position=corner)
288 
289  call pass_var(tmpz, sgdom, position=corner)
290  call extrapolate_metric(tmpz, 2*(g%jsc-g%jsd)+2, missing=999.)
291  do j=g%jsd,g%jed ; do i=g%isd,g%ied ; i2 = 2*i ; j2 = 2*j
292  g%geoLatT(i,j) = tmpz(i2-1,j2-1)
293  enddo ; enddo
294  do j=g%JsdB,g%JedB ; do i=g%IsdB,g%IedB ; i2 = 2*i ; j2 = 2*j
295  g%geoLatBu(i,j) = tmpz(i2,j2)
296  enddo ; enddo
297  do j=g%jsd,g%jed ; do i=g%IsdB,g%IedB ; i2 = 2*i ; j2 = 2*j
298  g%geoLatCu(i,j) = tmpz(i2,j2-1)
299  enddo ; enddo
300  do j=g%JsdB,g%JedB ; do i=g%isd,g%ied ; i2 = 2*i ; j2 = 2*j
301  g%geoLatCv(i,j) = tmpz(i2-1,j2)
302  enddo ; enddo
303 
304  ! Read DX,DY from the supergrid
305  tmpu(:,:) = 0. ; tmpv(:,:) = 0.
306  call mom_read_data(filename,'dx',tmpv,sgdom,position=north_face)
307  call mom_read_data(filename,'dy',tmpu,sgdom,position=east_face)
308  call pass_vector(tmpu, tmpv, sgdom, to_all+scalar_pair, cgrid_ne)
309  call extrapolate_metric(tmpv, 2*(g%jsc-g%jsd)+2, missing=0.)
310  call extrapolate_metric(tmpu, 2*(g%jsc-g%jsd)+2, missing=0.)
311 
312  do j=g%jsd,g%jed ; do i=g%isd,g%ied ; i2 = 2*i ; j2 = 2*j
313  dxt(i,j) = tmpv(i2-1,j2-1) + tmpv(i2,j2-1)
314  dyt(i,j) = tmpu(i2-1,j2-1) + tmpu(i2-1,j2)
315  enddo ; enddo
316 
317  do j=g%jsd,g%jed ; do i=g%IsdB,g%IedB ; i2 = 2*i ; j2 = 2*j
318  dxcu(i,j) = tmpv(i2,j2-1) + tmpv(i2+1,j2-1)
319  dycu(i,j) = tmpu(i2,j2-1) + tmpu(i2,j2)
320  enddo ; enddo
321 
322  do j=g%JsdB,g%JedB ; do i=g%isd,g%ied ; i2 = 2*i ; j2 = 2*j
323  dxcv(i,j) = tmpv(i2-1,j2) + tmpv(i2,j2)
324  dycv(i,j) = tmpu(i2-1,j2) + tmpu(i2-1,j2+1)
325  enddo ; enddo
326 
327  do j=g%JsdB,g%JedB ; do i=g%IsdB,g%IedB ; i2 = 2*i ; j2 = 2*j
328  dxbu(i,j) = tmpv(i2,j2) + tmpv(i2+1,j2)
329  dybu(i,j) = tmpu(i2,j2) + tmpu(i2,j2+1)
330  enddo ; enddo
331 
332  ! Read AREA from the supergrid
333  tmpt(:,:) = 0.
334  call mom_read_data(filename, 'area', tmpt, sgdom)
335  call pass_var(tmpt, sgdom)
336  call extrapolate_metric(tmpt, 2*(g%jsc-g%jsd)+2, missing=0.)
337 
338  do j=g%jsd,g%jed ; do i=g%isd,g%ied ; i2 = 2*i ; j2 = 2*j
339  areat(i,j) = (tmpt(i2-1,j2-1) + tmpt(i2,j2)) + &
340  (tmpt(i2-1,j2) + tmpt(i2,j2-1))
341  enddo ; enddo
342  do j=g%JsdB,g%JedB ; do i=g%IsdB,g%IedB ; i2 = 2*i ; j2 = 2*j
343  areabu(i,j) = (tmpt(i2,j2) + tmpt(i2+1,j2+1)) + &
344  (tmpt(i2,j2+1) + tmpt(i2+1,j2))
345  enddo ; enddo
346 
347  ni=sgdom%niglobal
348  nj=sgdom%njglobal
349  call mpp_deallocate_domain(sgdom%mpp_domain)
350  deallocate(sgdom%mpp_domain)
351 
352  call pass_vector(dycu, dxcv, g%Domain, to_all+scalar_pair, cgrid_ne)
353  call pass_vector(dxcu, dycv, g%Domain, to_all+scalar_pair, cgrid_ne)
354  call pass_vector(dxbu, dybu, g%Domain, to_all+scalar_pair, bgrid_ne)
355  call pass_var(areat, g%Domain)
356  call pass_var(areabu, g%Domain, position=corner)
357 
358  do i=g%isd,g%ied ; do j=g%jsd,g%jed
359  g%dxT(i,j) = m_to_l*dxt(i,j) ; g%dyT(i,j) = m_to_l*dyt(i,j) ; g%areaT(i,j) = m_to_l**2*areat(i,j)
360  enddo ; enddo
361  do i=g%IsdB,g%IedB ; do j=g%jsd,g%jed
362  g%dxCu(i,j) = m_to_l*dxcu(i,j) ; g%dyCu(i,j) = m_to_l*dycu(i,j)
363  enddo ; enddo
364  do i=g%isd,g%ied ; do j=g%JsdB,g%JedB
365  g%dxCv(i,j) = m_to_l*dxcv(i,j) ; g%dyCv(i,j) = m_to_l*dycv(i,j)
366  enddo ; enddo
367  do i=g%IsdB,g%IedB ; do j=g%JsdB,g%JedB
368  g%dxBu(i,j) = m_to_l*dxbu(i,j) ; g%dyBu(i,j) = m_to_l*dybu(i,j) ; g%areaBu(i,j) = m_to_l**2*areabu(i,j)
369  enddo ; enddo
370 
371  ! Construct axes for diagnostic output (only necessary because "ferret" uses
372  ! broken convention for interpretting netCDF files).
373  start(:) = 1 ; nread(:) = 1
374  start(2) = 2 ; nread(1) = ni+1 ; nread(2) = 2
375  allocate( tmpglbl(ni+1,2) )
376  if (is_root_pe()) &
377  call read_data(filename, "x", tmpglbl, start, nread, no_domain=.true.)
378  call broadcast(tmpglbl, 2*(ni+1), root_pe())
379 
380  ! I don't know why the second axis is 1 or 2 here. -RWH
381  do i=g%isg,g%ieg
382  g%gridLonT(i) = tmpglbl(2*(i-g%isg)+2,2)
383  enddo
384  ! Note that the dynamic grid always uses symmetric memory for the global
385  ! arrays G%gridLatB and G%gridLonB.
386  do i=g%isg-1,g%ieg
387  g%gridLonB(i) = tmpglbl(2*(i-g%isg)+3,1)
388  enddo
389  deallocate( tmpglbl )
390 
391  allocate( tmpglbl(1, nj+1) )
392  start(:) = 1 ; nread(:) = 1
393  start(1) = int(ni/4)+1 ; nread(2) = nj+1
394  if (is_root_pe()) &
395  call read_data(filename, "y", tmpglbl, start, nread, no_domain=.true.)
396  call broadcast(tmpglbl, nj+1, root_pe())
397 
398  do j=g%jsg,g%jeg
399  g%gridLatT(j) = tmpglbl(1,2*(j-g%jsg)+2)
400  enddo
401  do j=g%jsg-1,g%jeg
402  g%gridLatB(j) = tmpglbl(1,2*(j-g%jsg)+3)
403  enddo
404  deallocate( tmpglbl )
405 
406  call calltree_leave("set_grid_metrics_from_mosaic()")

◆ set_grid_metrics_mercator()

subroutine mom_grid_initialize::set_grid_metrics_mercator ( type(dyn_horgrid_type), intent(inout)  G,
type(param_file_type), intent(in)  param_file,
type(unit_scale_type), intent(in), optional  US 
)
private

Calculate the values of the metric terms that might be used and save them in arrays.

Within this subroutine, the x- and y- grid spacings and their inverses and the cell areas centered on h, q, u, and v points are calculated, as are the geographic locations of each of these 4 sets of points.

Parameters
[in,out]gThe dynamic horizontal grid type
[in]param_fileParameter file structure
[in]usA dimensional unit scaling type

Definition at line 696 of file MOM_grid_initialize.F90.

697  type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type
698  type(param_file_type), intent(in) :: param_file !< Parameter file structure
699  type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type
700  ! Local variables
701  integer :: i, j, isd, ied, jsd, jed
702  integer :: I_off, J_off
703  type(GPS) :: GP
704  character(len=128) :: warnmesg
705  character(len=48) :: mdl = "MOM_grid_init set_grid_metrics_mercator"
706  real :: PI, PI_2! PI = 3.1415926... as 4*atan(1), PI_2 = (PI) /2.0
707  real :: y_q, y_h, jd, x_q, x_h, id
708  real, dimension(G%isd:G%ied,G%jsd:G%jed) :: &
709  xh, yh ! Latitude and longitude of h points in radians.
710  real, dimension(G%IsdB:G%IedB,G%jsd:G%jed) :: &
711  xu, yu ! Latitude and longitude of u points in radians.
712  real, dimension(G%isd:G%ied,G%JsdB:G%JedB) :: &
713  xv, yv ! Latitude and longitude of v points in radians.
714  real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB) :: &
715  xq, yq ! Latitude and longitude of q points in radians.
716  real :: fnRef ! fnRef is the value of Int_dj_dy or
717  ! Int_dj_dy at a latitude or longitude that is
718  real :: jRef, iRef ! being set to be at grid index jRef or iRef.
719  real :: m_to_L ! A unit conversion factor [L m-1 ~> nondim]
720  integer :: itt1, itt2
721  logical :: debug = .false., simple_area = .true.
722  integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, IsdB, IedB, JsdB, JedB
723 
724  ! All of the metric terms should be defined over the domain from
725  ! isd to ied. Outside of the physical domain, both the metrics
726  ! and their inverses may be set to zero.
727  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
728  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
729  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
730  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
731  i_off = g%idg_offset ; j_off = g%jdg_offset
732 
733  gp%niglobal = g%Domain%niglobal
734  gp%njglobal = g%Domain%njglobal
735 
736  call calltree_enter("set_grid_metrics_mercator(), MOM_grid_initialize.F90")
737 
738  m_to_l = 1.0 ; if (present(us)) m_to_l = us%m_to_L
739  ! Calculate the values of the metric terms that might be used
740  ! and save them in arrays.
741  pi = 4.0*atan(1.0) ; pi_2 = 0.5*pi
742 
743  call get_param(param_file, mdl, "SOUTHLAT", gp%south_lat, &
744  "The southern latitude of the domain.", units="degrees", &
745  fail_if_missing=.true.)
746  call get_param(param_file, mdl, "LENLAT", gp%len_lat, &
747  "The latitudinal length of the domain.", units="degrees", &
748  fail_if_missing=.true.)
749  call get_param(param_file, mdl, "WESTLON", gp%west_lon, &
750  "The western longitude of the domain.", units="degrees", &
751  default=0.0)
752  call get_param(param_file, mdl, "LENLON", gp%len_lon, &
753  "The longitudinal length of the domain.", units="degrees", &
754  fail_if_missing=.true.)
755  call get_param(param_file, mdl, "RAD_EARTH", gp%Rad_Earth, &
756  "The radius of the Earth.", units="m", default=6.378e6)
757  g%south_lat = gp%south_lat ; g%len_lat = gp%len_lat
758  g%west_lon = gp%west_lon ; g%len_lon = gp%len_lon
759  g%Rad_Earth = gp%Rad_Earth
760  call get_param(param_file, mdl, "ISOTROPIC", gp%isotropic, &
761  "If true, an isotropic grid on a sphere (also known as "//&
762  "a Mercator grid) is used. With an isotropic grid, the "//&
763  "meridional extent of the domain (LENLAT), the zonal "//&
764  "extent (LENLON), and the number of grid points in each "//&
765  "direction are _not_ independent. In MOM the meridional "//&
766  "extent is determined to fit the zonal extent and the "//&
767  "number of grid points, while grid is perfectly isotropic.", &
768  default=.false.)
769  call get_param(param_file, mdl, "EQUATOR_REFERENCE", gp%equator_reference, &
770  "If true, the grid is defined to have the equator at the "//&
771  "nearest q- or h- grid point to (-LOWLAT*NJGLOBAL/LENLAT).", &
772  default=.true.)
773  call get_param(param_file, mdl, "LAT_ENHANCE_FACTOR", gp%Lat_enhance_factor, &
774  "The amount by which the meridional resolution is "//&
775  "enhanced within LAT_EQ_ENHANCE of the equator.", &
776  units="nondim", default=1.0)
777  call get_param(param_file, mdl, "LAT_EQ_ENHANCE", gp%Lat_eq_enhance, &
778  "The latitude range to the north and south of the equator "//&
779  "over which the resolution is enhanced.", units="degrees", &
780  default=0.0)
781 
782  ! With an isotropic grid, the north-south extent of the domain,
783  ! the east-west extent, and the number of grid points in each
784  ! direction are _not_ independent. Here the north-south extent
785  ! will be determined to fit the east-west extent and the number of
786  ! grid points. The grid is perfectly isotropic.
787  if (gp%equator_reference) then
788  ! With the following expression, the equator will always be placed
789  ! on either h or q points, in a position consistent with the ratio
790  ! GP%south_lat to GP%len_lat.
791  jref = (g%jsg-1) + 0.5*floor(gp%njglobal*((-1.0*gp%south_lat*2.0)/gp%len_lat)+0.5)
792  fnref = int_dj_dy(0.0, gp)
793  else
794  ! The following line sets the reference latitude GP%south_lat at j=js-1 (or -2?)
795  jref = (g%jsg-1)
796  fnref = int_dj_dy((gp%south_lat*pi/180.0), gp)
797  endif
798 
799  ! These calculations no longer depend on the the order in which they
800  ! are performed because they all use the same (poor) starting guess and
801  ! iterate to convergence.
802  ! Note that the dynamic grid always uses symmetric memory for the global
803  ! arrays G%gridLatB and G%gridLonB.
804  do j=g%jsg-1,g%jeg
805  jd = fnref + (j - jref)
806  y_q = find_root(int_dj_dy, dy_dj, gp, jd, 0.0, -1.0*pi_2, pi_2, itt2)
807  g%gridLatB(j) = y_q*180.0/pi
808  ! if (is_root_pe()) &
809  ! write(*, '("J, y_q = ",I4,ES14.4," itts = ",I4)') j, y_q, itt2
810  enddo
811  do j=g%jsg,g%jeg
812  jd = fnref + (j - jref) - 0.5
813  y_h = find_root(int_dj_dy, dy_dj, gp, jd, 0.0, -1.0*pi_2, pi_2, itt1)
814  g%gridLatT(j) = y_h*180.0/pi
815  ! if (is_root_pe()) &
816  ! write(*, '("j, y_h = ",I4,ES14.4," itts = ",I4)') j, y_h, itt1
817  enddo
818  do j=jsdb+j_off,jedb+j_off
819  jd = fnref + (j - jref)
820  y_q = find_root(int_dj_dy, dy_dj, gp, jd, 0.0, -1.0*pi_2, pi_2, itt2)
821  do i=isdb,iedb ; yq(i,j-j_off) = y_q ; enddo
822  do i=isd,ied ; yv(i,j-j_off) = y_q ; enddo
823  enddo
824  do j=jsd+j_off,jed+j_off
825  jd = fnref + (j - jref) - 0.5
826  y_h = find_root(int_dj_dy, dy_dj, gp, jd, 0.0, -1.0*pi_2, pi_2, itt1)
827  if ((j >= jsd+j_off) .and. (j <= jed+j_off)) then
828  do i=isd,ied ; yh(i,j-j_off) = y_h ; enddo
829  do i=isdb,iedb ; yu(i,j-j_off) = y_h ; enddo
830  endif
831  enddo
832 
833  ! Determine the longitudes of the various points.
834 
835  ! These two lines place the western edge of the domain at GP%west_lon.
836  iref = (g%isg-1) + gp%niglobal
837  fnref = int_di_dx(((gp%west_lon+gp%len_lon)*pi/180.0), gp)
838 
839  ! These calculations no longer depend on the the order in which they
840  ! are performed because they all use the same (poor) starting guess and
841  ! iterate to convergence.
842  do i=g%isg-1,g%ieg
843  id = fnref + (i - iref)
844  x_q = find_root(int_di_dx, dx_di, gp, id, 0.0, -4.0*pi, 4.0*pi, itt2)
845  g%gridLonB(i) = x_q*180.0/pi
846  enddo
847  do i=g%isg,g%ieg
848  id = fnref + (i - iref) - 0.5
849  x_h = find_root(int_di_dx, dx_di, gp, id, 0.0, -4.0*pi, 4.0*pi, itt1)
850  g%gridLonT(i) = x_h*180.0/pi
851  enddo
852  do i=isdb+i_off,iedb+i_off
853  id = fnref + (i - iref)
854  x_q = find_root(int_di_dx, dx_di, gp, id, 0.0, -4.0*pi, 4.0*pi, itt2)
855  do j=jsdb,jedb ; xq(i-i_off,j) = x_q ; enddo
856  do j=jsd,jed ; xu(i-i_off,j) = x_q ; enddo
857  enddo
858  do i=isd+i_off,ied+i_off
859  id = fnref + (i - iref) - 0.5
860  x_h = find_root(int_di_dx, dx_di, gp, id, 0.0, -4.0*pi, 4.0*pi, itt1)
861  do j=jsd,jed ; xh(i-i_off,j) = x_h ; enddo
862  do j=jsdb,jedb ; xv(i-i_off,j) = x_h ; enddo
863  enddo
864 
865  do j=jsdb,jedb ; do i=isdb,iedb
866  g%geoLonBu(i,j) = xq(i,j)*180.0/pi
867  g%geoLatBu(i,j) = yq(i,j)*180.0/pi
868  g%dxBu(i,j) = m_to_l*ds_di(xq(i,j), yq(i,j), gp)
869  g%dyBu(i,j) = m_to_l*ds_dj(xq(i,j), yq(i,j), gp)
870 
871  g%areaBu(i,j) = g%dxBu(i,j) * g%dyBu(i,j)
872  g%IareaBu(i,j) = 1.0 / (g%areaBu(i,j))
873  enddo ; enddo
874 
875  do j=jsd,jed ; do i=isd,ied
876  g%geoLonT(i,j) = xh(i,j)*180.0/pi
877  g%geoLatT(i,j) = yh(i,j)*180.0/pi
878  g%dxT(i,j) = m_to_l*ds_di(xh(i,j), yh(i,j), gp)
879  g%dyT(i,j) = m_to_l*ds_dj(xh(i,j), yh(i,j), gp)
880 
881  g%areaT(i,j) = g%dxT(i,j)*g%dyT(i,j)
882  g%IareaT(i,j) = 1.0 / (g%areaT(i,j))
883  enddo ; enddo
884 
885  do j=jsd,jed ; do i=isdb,iedb
886  g%geoLonCu(i,j) = xu(i,j)*180.0/pi
887  g%geoLatCu(i,j) = yu(i,j)*180.0/pi
888  g%dxCu(i,j) = m_to_l*ds_di(xu(i,j), yu(i,j), gp)
889  g%dyCu(i,j) = m_to_l*ds_dj(xu(i,j), yu(i,j), gp)
890  enddo ; enddo
891 
892  do j=jsdb,jedb ; do i=isd,ied
893  g%geoLonCv(i,j) = xv(i,j)*180.0/pi
894  g%geoLatCv(i,j) = yv(i,j)*180.0/pi
895  g%dxCv(i,j) = m_to_l*ds_di(xv(i,j), yv(i,j), gp)
896  g%dyCv(i,j) = m_to_l*ds_dj(xv(i,j), yv(i,j), gp)
897  enddo ; enddo
898 
899  if (.not.simple_area) then
900  do j=jsdb+1,jed ; do i=isdb+1,ied
901  g%areaT(i,j) = m_to_l**2*gp%Rad_Earth**2 * &
902  (dl(xq(i-1,j-1),xq(i-1,j),yq(i-1,j-1),yq(i-1,j)) + &
903  (dl(xq(i-1,j),xq(i,j),yq(i-1,j),yq(i,j)) + &
904  (dl(xq(i,j),xq(i,j-1),yq(i,j),yq(i,j-1)) + &
905  dl(xq(i,j-1),xq(i-1,j-1),yq(i,j-1),yq(i-1,j-1)))))
906  enddo ;enddo
907  if ((isdb == isd) .or. (jsdb == jsq)) then
908  ! Fill in row and column 1 to calculate the area in the southernmost
909  ! and westernmost land cells when we are not using symmetric memory.
910  ! The pass_var call updates these values if they are not land cells.
911  g%areaT(isd+1,jsd) = g%areaT(isd+1,jsd+1)
912  do j=jsd,jed ; g%areaT(isd,j) = g%areaT(isd+1,j) ; enddo
913  do i=isd,ied ; g%areaT(i,jsd) = g%areaT(i,jsd+1) ; enddo
914  ! Now replace the data in the halos, if value values exist.
915  call pass_var(g%areaT,g%Domain)
916  endif
917  do j=jsd,jed ; do i=isd,ied
918  g%IareaT(i,j) = 1.0 / (g%areaT(i,j))
919  enddo ; enddo
920  endif
921 
922  call calltree_leave("set_grid_metrics_mercator()")

◆ set_grid_metrics_spherical()

subroutine mom_grid_initialize::set_grid_metrics_spherical ( type(dyn_horgrid_type), intent(inout)  G,
type(param_file_type), intent(in)  param_file,
type(unit_scale_type), intent(in), optional  US 
)
private

Calculate the values of the metric terms that might be used and save them in arrays.

Within this subroutine, the x- and y- grid spacings and their inverses and the cell areas centered on h, q, u, and v points are calculated, as are the geographic locations of each of these 4 sets of points.

Parameters
[in,out]gThe dynamic horizontal grid type
[in]param_fileParameter file structure
[in]usA dimensional unit scaling type

Definition at line 557 of file MOM_grid_initialize.F90.

558  type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type
559  type(param_file_type), intent(in) :: param_file !< Parameter file structure
560  type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type
561  ! Local variables
562  real :: PI, PI_180! PI = 3.1415926... as 4*atan(1)
563  integer :: i, j, isd, ied, jsd, jed
564  integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, IsdB, IedB, JsdB, JedB
565  integer :: i_offset, j_offset
566  real :: grid_latT(G%jsd:G%jed), grid_latB(G%JsdB:G%JedB)
567  real :: grid_lonT(G%isd:G%ied), grid_lonB(G%IsdB:G%IedB)
568  real :: dLon,dLat,latitude,longitude,dL_di
569  real :: m_to_L ! A unit conversion factor [L m-1 ~> nondim]
570  character(len=48) :: mdl = "MOM_grid_init set_grid_metrics_spherical"
571 
572  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
573  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
574  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
575  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
576  i_offset = g%idg_offset ; j_offset = g%jdg_offset
577 
578  call calltree_enter("set_grid_metrics_spherical(), MOM_grid_initialize.F90")
579  m_to_l = 1.0 ; if (present(us)) m_to_l = us%m_to_L
580 
581 ! Calculate the values of the metric terms that might be used
582 ! and save them in arrays.
583  pi = 4.0*atan(1.0); pi_180 = atan(1.0)/45.
584 
585  call get_param(param_file, mdl, "SOUTHLAT", g%south_lat, &
586  "The southern latitude of the domain.", units="degrees", &
587  fail_if_missing=.true.)
588  call get_param(param_file, mdl, "LENLAT", g%len_lat, &
589  "The latitudinal length of the domain.", units="degrees", &
590  fail_if_missing=.true.)
591  call get_param(param_file, mdl, "WESTLON", g%west_lon, &
592  "The western longitude of the domain.", units="degrees", &
593  default=0.0)
594  call get_param(param_file, mdl, "LENLON", g%len_lon, &
595  "The longitudinal length of the domain.", units="degrees", &
596  fail_if_missing=.true.)
597  call get_param(param_file, mdl, "RAD_EARTH", g%Rad_Earth, &
598  "The radius of the Earth.", units="m", default=6.378e6)
599 
600  dlon = g%len_lon/g%Domain%niglobal
601  dlat = g%len_lat/g%Domain%njglobal
602 
603  ! Note that the dynamic grid always uses symmetric memory for the global
604  ! arrays G%gridLatB and G%gridLonB.
605  do j=g%jsg-1,g%jeg
606  latitude = g%south_lat + dlat*(real(j-(g%jsg-1)))
607  g%gridLatB(j) = min(max(latitude,-90.),90.)
608  enddo
609  do j=g%jsg,g%jeg
610  latitude = g%south_lat + dlat*(real(j-g%jsg)+0.5)
611  g%gridLatT(j) = min(max(latitude,-90.),90.)
612  enddo
613  do i=g%isg-1,g%ieg
614  g%gridLonB(i) = g%west_lon + dlon*(real(i-(g%isg-1)))
615  enddo
616  do i=g%isg,g%ieg
617  g%gridLonT(i) = g%west_lon + dlon*(real(i-g%isg)+0.5)
618  enddo
619 
620  do j=jsdb,jedb
621  latitude = g%south_lat + dlat* real(j+j_offset-(g%jsg-1))
622  grid_latb(j) = min(max(latitude,-90.),90.)
623  enddo
624  do j=jsd,jed
625  latitude = g%south_lat + dlat*(real(j+j_offset-g%jsg)+0.5)
626  grid_latt(j) = min(max(latitude,-90.),90.)
627  enddo
628  do i=isdb,iedb
629  grid_lonb(i) = g%west_lon + dlon*real(i+i_offset-(g%isg-1))
630  enddo
631  do i=isd,ied
632  grid_lont(i) = g%west_lon + dlon*(real(i+i_offset-g%isg)+0.5)
633  enddo
634 
635  dl_di = (g%len_lon * 4.0*atan(1.0)) / (180.0 * g%Domain%niglobal)
636  do j=jsdb,jedb ; do i=isdb,iedb
637  g%geoLonBu(i,j) = grid_lonb(i)
638  g%geoLatBu(i,j) = grid_latb(j)
639 
640  ! The following line is needed to reproduce the solution from
641  ! set_grid_metrics_mercator when used to generate a simple spherical grid.
642  g%dxBu(i,j) = m_to_l*g%Rad_Earth * cos( g%geoLatBu(i,j)*pi_180 ) * dl_di
643 ! G%dxBu(I,J) = m_to_L*G%Rad_Earth * dLon*PI_180 * COS( G%geoLatBu(I,J)*PI_180 )
644  g%dyBu(i,j) = m_to_l*g%Rad_Earth * dlat*pi_180
645  g%areaBu(i,j) = g%dxBu(i,j) * g%dyBu(i,j)
646  enddo ; enddo
647 
648  do j=jsdb,jedb ; do i=isd,ied
649  g%geoLonCv(i,j) = grid_lont(i)
650  g%geoLatCv(i,j) = grid_latb(j)
651 
652  ! The following line is needed to reproduce the solution from
653  ! set_grid_metrics_mercator when used to generate a simple spherical grid.
654  g%dxCv(i,j) = m_to_l*g%Rad_Earth * cos( g%geoLatCv(i,j)*pi_180 ) * dl_di
655 ! G%dxCv(i,J) = m_to_L*G%Rad_Earth * (dLon*PI_180) * COS( G%geoLatCv(i,J)*PI_180 )
656  g%dyCv(i,j) = m_to_l*g%Rad_Earth * dlat*pi_180
657  enddo ; enddo
658 
659  do j=jsd,jed ; do i=isdb,iedb
660  g%geoLonCu(i,j) = grid_lonb(i)
661  g%geoLatCu(i,j) = grid_latt(j)
662 
663  ! The following line is needed to reproduce the solution from
664  ! set_grid_metrics_mercator when used to generate a simple spherical grid.
665  g%dxCu(i,j) = m_to_l*g%Rad_Earth * cos( g%geoLatCu(i,j)*pi_180 ) * dl_di
666 ! G%dxCu(I,j) = m_to_L*G%Rad_Earth * dLon*PI_180 * COS( latitude )
667  g%dyCu(i,j) = m_to_l*g%Rad_Earth * dlat*pi_180
668  enddo ; enddo
669 
670  do j=jsd,jed ; do i=isd,ied
671  g%geoLonT(i,j) = grid_lont(i)
672  g%geoLatT(i,j) = grid_latt(j)
673 
674  ! The following line is needed to reproduce the solution from
675  ! set_grid_metrics_mercator when used to generate a simple spherical grid.
676  g%dxT(i,j) = m_to_l*g%Rad_Earth * cos( g%geoLatT(i,j)*pi_180 ) * dl_di
677 ! G%dxT(i,j) = G%Rad_Earth * dLon*PI_180 * COS( latitude )
678  g%dyT(i,j) = m_to_l*g%Rad_Earth * dlat*pi_180
679 
680 ! latitude = G%geoLatCv(i,J)*PI_180 ! In radians
681 ! dL_di = G%geoLatCv(i,max(jsd,J-1))*PI_180 ! In radians
682 ! G%areaT(i,j) = m_to_L**2 * Rad_Earth**2*dLon*dLat*ABS(SIN(latitude)-SIN(dL_di))
683  g%areaT(i,j) = g%dxT(i,j) * g%dyT(i,j)
684  enddo ; enddo
685 
686  call calltree_leave("set_grid_metrics_spherical()")