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 1219 of file MOM_grid_initialize.F90.

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

Referenced by initialize_masks().

Here is the caller graph for this function:

◆ 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 956 of file MOM_grid_initialize.F90.

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

Referenced by set_grid_metrics_mercator().

Here is the caller graph for this function:

◆ 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 926 of file MOM_grid_initialize.F90.

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

References dx_di().

Referenced by set_grid_metrics_mercator().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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 940 of file MOM_grid_initialize.F90.

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

References dy_dj().

Referenced by set_grid_metrics_mercator().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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 1091 of file MOM_grid_initialize.F90.

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

Referenced by ds_di(), and set_grid_metrics_mercator().

Here is the caller graph for this function:

◆ 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 1113 of file MOM_grid_initialize.F90.

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

Referenced by ds_dj(), and set_grid_metrics_mercator().

Here is the caller graph for this function:

◆ 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 1185 of file MOM_grid_initialize.F90.

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

Referenced by set_grid_metrics_from_mosaic().

Here is the caller graph for this function:

◆ 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 979 of file MOM_grid_initialize.F90.

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

Referenced by set_grid_metrics_mercator().

Here is the caller graph for this function:

◆ 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 116 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, haloshift=halo, scale=l_to_m)
129 
130  call uvchksum(trim(parent)//': dxC[uv]', g%dxCu, g%dyCv, g%HI, haloshift=halo, scale=l_to_m)
131 
132  call uvchksum(trim(parent)//': dxC[uv]', g%dyCu, g%dxCv, g%HI, haloshift=halo, scale=l_to_m)
133 
134  call bchksum_pair(trim(parent)//': dxB[uv]', g%dxBu, g%dyBu, g%HI, haloshift=halo, scale=l_to_m)
135 
136  call hchksum_pair(trim(parent)//': Id[xy]T', g%IdxT, g%IdyT, g%HI, haloshift=halo, scale=m_to_l)
137 
138  call uvchksum(trim(parent)//': Id[xy]C[uv]', g%IdxCu, g%IdyCv, g%HI, haloshift=halo, scale=m_to_l)
139 
140  call uvchksum(trim(parent)//': Id[xy]C[uv]', g%IdyCu, g%IdxCv, g%HI, haloshift=halo, scale=m_to_l)
141 
142  call bchksum_pair(trim(parent)//': Id[xy]B[uv]', g%IdxBu, g%IdyBu, g%HI, haloshift=halo, scale=m_to_l)
143 
144  call hchksum(g%areaT, trim(parent)//': areaT',g%HI, haloshift=halo, scale=l_to_m**2)
145  call bchksum(g%areaBu, trim(parent)//': areaBu',g%HI, haloshift=halo, scale=l_to_m**2)
146 
147  call hchksum(g%IareaT, trim(parent)//': IareaT',g%HI, haloshift=halo, scale=m_to_l**2)
148  call bchksum(g%IareaBu, trim(parent)//': IareaBu',g%HI, haloshift=halo, scale=m_to_l**2)
149 
150  call hchksum(g%geoLonT,trim(parent)//': geoLonT',g%HI, haloshift=halo)
151  call hchksum(g%geoLatT,trim(parent)//': geoLatT',g%HI, haloshift=halo)
152 
153  call bchksum(g%geoLonBu, trim(parent)//': geoLonBu',g%HI, haloshift=halo)
154  call bchksum(g%geoLatBu, trim(parent)//': geoLatBu',g%HI, haloshift=halo)
155 
156  call uvchksum(trim(parent)//': geoLonC[uv]', g%geoLonCu, g%geoLonCv, g%HI, haloshift=halo)
157 
158  call uvchksum(trim(parent)//': geoLatC[uv]', g%geoLatCu, g%geoLatCv, g%HI, haloshift=halo)
159 

Referenced by set_grid_metrics().

Here is the caller graph for this function:

◆ 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 1235 of file MOM_grid_initialize.F90.

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

References adcroft_reciprocal(), mom_error_handler::calltree_enter(), and mom_error_handler::calltree_leave().

Referenced by mom_fixed_initialization::mom_initialize_fixed().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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 1102 of file MOM_grid_initialize.F90.

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

Referenced by set_grid_metrics_mercator().

Here is the caller graph for this function:

◆ 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 1142 of file MOM_grid_initialize.F90.

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

Referenced by set_grid_metrics_mercator().

Here is the caller graph for this function:

◆ 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 63 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()")

References mom_error_handler::calltree_enter(), mom_error_handler::calltree_leave(), grid_metrics_chksum(), mom_dyn_horgrid::set_derived_dyn_horgrid(), set_grid_metrics_cartesian(), set_grid_metrics_from_mosaic(), set_grid_metrics_mercator(), and set_grid_metrics_spherical().

Referenced by mom_oda_driver_mod::init_oda(), mom_ice_shelf::initialize_ice_shelf(), and mom_fixed_initialization::mom_initialize_fixed().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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 418 of file MOM_grid_initialize.F90.

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

References mom_error_handler::calltree_enter(), and mom_error_handler::calltree_leave().

Referenced by set_grid_metrics().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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 166 of file MOM_grid_initialize.F90.

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

References mom_error_handler::calltree_enter(), mom_error_handler::calltree_leave(), and extrapolate_metric().

Referenced by set_grid_metrics().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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 695 of file MOM_grid_initialize.F90.

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

References mom_error_handler::calltree_enter(), mom_error_handler::calltree_leave(), dl(), ds_di(), ds_dj(), dx_di(), dy_dj(), find_root(), int_di_dx(), and int_dj_dy().

Referenced by set_grid_metrics().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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 556 of file MOM_grid_initialize.F90.

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

References mom_error_handler::calltree_enter(), and mom_error_handler::calltree_leave().

Referenced by set_grid_metrics().

Here is the call graph for this function:
Here is the caller graph for this function: