MOM6
Functions/Subroutines
mom_shared_initialization Module Reference

Detailed Description

Code that initializes fixed aspects of the model grid, such as horizontal grid metrics, topography and Coriolis, and can be shared between components.

Functions/Subroutines

subroutine, public mom_shared_init_init (PF)
 MOM_shared_init_init just writes the code version. More...
 
subroutine, public mom_initialize_rotation (f, G, PF, US)
 MOM_initialize_rotation makes the appropriate call to set up the Coriolis parameter. More...
 
subroutine, public mom_calculate_grad_coriolis (dF_dx, dF_dy, G, US)
 Calculates the components of grad f (Coriolis parameter) More...
 
real function, public diagnosemaximumdepth (D, G)
 Return the global maximum ocean bottom depth in the same units as the input depth. More...
 
subroutine, public initialize_topography_from_file (D, G, param_file, US)
 Read gridded depths from file. More...
 
subroutine, public apply_topography_edits_from_file (D, G, param_file, US)
 Applies a list of topography overrides read from a netcdf file. More...
 
subroutine, public initialize_topography_named (D, G, param_file, topog_config, max_depth, US)
 initialize the bathymetry based on one of several named idealized configurations More...
 
subroutine, public limit_topography (D, G, param_file, max_depth, US)
 limit_topography ensures that min_depth < D(x,y) < max_depth More...
 
subroutine, public set_rotation_planetary (f, G, param_file, US)
 This subroutine sets up the Coriolis parameter for a sphere. More...
 
subroutine, public set_rotation_beta_plane (f, G, param_file, US)
 This subroutine sets up the Coriolis parameter for a beta-plane or f-plane. More...
 
subroutine, public initialize_grid_rotation_angle (G, PF)
 initialize_grid_rotation_angle initializes the arrays with the sine and cosine of the angle between logical north on the grid and true north. More...
 
real function modulo_around_point (x, xc, Lx)
 Return the modulo value of x in an interval [xc-(Lx/2) xc+(Lx/2)] If Lx<=0, then it returns x without applying modulo arithmetic. More...
 
subroutine, public reset_face_lengths_named (G, param_file, name, US)
 This subroutine sets the open face lengths at selected points to restrict passages to their observed widths based on a named set of sizes. More...
 
subroutine, public reset_face_lengths_file (G, param_file, US)
 This subroutine sets the open face lengths at selected points to restrict passages to their observed widths from a arrays read from a file. More...
 
subroutine, public reset_face_lengths_list (G, param_file, US)
 This subroutine sets the open face lengths at selected points to restrict passages to their observed widths from a list read from a file. More...
 
subroutine, public read_face_length_list (iounit, filename, num_lines, lines)
 This subroutine reads and counts the non-blank lines in the face length list file, after removing comments. More...
 
subroutine, public set_velocity_depth_max (G)
 Set the bathymetry at velocity points to be the maximum of the depths at the neighoring tracer points. More...
 
subroutine, public set_velocity_depth_min (G)
 Set the bathymetry at velocity points to be the minimum of the depths at the neighoring tracer points. More...
 
subroutine, public compute_global_grid_integrals (G, US)
 Pre-compute global integrals of grid quantities (like masked ocean area) for later use in reporting diagnostics. More...
 
subroutine, public write_ocean_geometry_file (G, param_file, directory, geom_file, US)
 Write out a file describing the topography, Coriolis parameter, grid locations and various other fixed fields from the grid. More...
 

Function/Subroutine Documentation

◆ apply_topography_edits_from_file()

subroutine, public mom_shared_initialization::apply_topography_edits_from_file ( real, dimension(g%isd:g%ied,g%jsd:g%jed), intent(inout)  D,
type(dyn_horgrid_type), intent(in)  G,
type(param_file_type), intent(in)  param_file,
type(unit_scale_type), intent(in), optional  US 
)

Applies a list of topography overrides read from a netcdf file.

Parameters
[in]gThe dynamic horizontal grid type
[in,out]dOcean bottom depth in m or Z if US is present
[in]param_fileParameter file structure
[in]usA dimensional unit scaling type

Definition at line 186 of file MOM_shared_initialization.F90.

186  type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type
187  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
188  intent(inout) :: D !< Ocean bottom depth in m or Z if US is present
189  type(param_file_type), intent(in) :: param_file !< Parameter file structure
190  type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type
191 
192  ! Local variables
193  real :: m_to_Z ! A dimensional rescaling factor.
194  character(len=200) :: topo_edits_file, inputdir ! Strings for file/path
195  character(len=40) :: mdl = "apply_topography_edits_from_file" ! This subroutine's name.
196  integer :: n_edits, n, ashape(5), i, j, ncid, id, ncstatus, iid, jid, zid
197  integer, dimension(:), allocatable :: ig, jg
198  real, dimension(:), allocatable :: new_depth
199 
200  call calltree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")
201 
202  m_to_z = 1.0 ; if (present(us)) m_to_z = us%m_to_Z
203 
204  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
205  inputdir = slasher(inputdir)
206  call get_param(param_file, mdl, "TOPO_EDITS_FILE", topo_edits_file, &
207  "The file from which to read a list of i,j,z topography overrides.", &
208  default="")
209 
210  if (len_trim(topo_edits_file)==0) return
211 
212  topo_edits_file = trim(inputdir)//trim(topo_edits_file)
213  if (.not.file_exists(topo_edits_file, g%Domain)) call mom_error(fatal, &
214  'initialize_topography_from_file: Unable to open '//trim(topo_edits_file))
215 
216  ncstatus = nf90_open(trim(topo_edits_file), nf90_nowrite, ncid)
217  if (ncstatus /= nf90_noerr) call mom_error(fatal, 'apply_topography_edits_from_file: '//&
218  'Failed to open '//trim(topo_edits_file))
219 
220  ! Get nEdits
221  ncstatus = nf90_inq_dimid(ncid, 'nEdits', id)
222  if (ncstatus /= nf90_noerr) call mom_error(fatal, 'apply_topography_edits_from_file: '//&
223  'Failed to inq_dimid nEdits for '//trim(topo_edits_file))
224  ncstatus = nf90_inquire_dimension(ncid, id, len=n_edits)
225  if (ncstatus /= nf90_noerr) call mom_error(fatal, 'apply_topography_edits_from_file: '//&
226  'Failed to inquire_dimension nEdits for '//trim(topo_edits_file))
227 
228  ! Read ni
229  ncstatus = nf90_inq_varid(ncid, 'ni', id)
230  if (ncstatus /= nf90_noerr) call mom_error(fatal, 'apply_topography_edits_from_file: '//&
231  'Failed to inq_varid ni for '//trim(topo_edits_file))
232  ncstatus = nf90_get_var(ncid, id, i)
233  if (ncstatus /= nf90_noerr) call mom_error(fatal, 'apply_topography_edits_from_file: '//&
234  'Failed to get_var ni for '//trim(topo_edits_file))
235  if (i /= g%ieg) call mom_error(fatal, 'apply_topography_edits_from_file: '//&
236  'Incompatible i-dimension of grid in '//trim(topo_edits_file))
237 
238  ! Read nj
239  ncstatus = nf90_inq_varid(ncid, 'nj', id)
240  if (ncstatus /= nf90_noerr) call mom_error(fatal, 'apply_topography_edits_from_file: '//&
241  'Failed to inq_varid nj for '//trim(topo_edits_file))
242  ncstatus = nf90_get_var(ncid, id, j)
243  if (ncstatus /= nf90_noerr) call mom_error(fatal, 'apply_topography_edits_from_file: '//&
244  'Failed to get_var nj for '//trim(topo_edits_file))
245  if (j /= g%jeg) call mom_error(fatal, 'apply_topography_edits_from_file: '//&
246  'Incompatible j-dimension of grid in '//trim(topo_edits_file))
247 
248  ! Read iEdit
249  ncstatus = nf90_inq_varid(ncid, 'iEdit', id)
250  if (ncstatus /= nf90_noerr) call mom_error(fatal, 'apply_topography_edits_from_file: '//&
251  'Failed to inq_varid iEdit for '//trim(topo_edits_file))
252  allocate(ig(n_edits))
253  ncstatus = nf90_get_var(ncid, id, ig)
254  if (ncstatus /= nf90_noerr) call mom_error(fatal, 'apply_topography_edits_from_file: '//&
255  'Failed to get_var iEdit for '//trim(topo_edits_file))
256 
257  ! Read jEdit
258  ncstatus = nf90_inq_varid(ncid, 'jEdit', id)
259  if (ncstatus /= nf90_noerr) call mom_error(fatal, 'apply_topography_edits_from_file: '//&
260  'Failed to inq_varid jEdit for '//trim(topo_edits_file))
261  allocate(jg(n_edits))
262  ncstatus = nf90_get_var(ncid, id, jg)
263  if (ncstatus /= nf90_noerr) call mom_error(fatal, 'apply_topography_edits_from_file: '//&
264  'Failed to get_var jEdit for '//trim(topo_edits_file))
265 
266  ! Read zEdit
267  ncstatus = nf90_inq_varid(ncid, 'zEdit', id)
268  if (ncstatus /= nf90_noerr) call mom_error(fatal, 'apply_topography_edits_from_file: '//&
269  'Failed to inq_varid zEdit for '//trim(topo_edits_file))
270  allocate(new_depth(n_edits))
271  ncstatus = nf90_get_var(ncid, id, new_depth)
272  if (ncstatus /= nf90_noerr) call mom_error(fatal, 'apply_topography_edits_from_file: '//&
273  'Failed to get_var zEdit for '//trim(topo_edits_file))
274 
275  ! Close file
276  ncstatus = nf90_close(ncid)
277  if (ncstatus /= nf90_noerr) call mom_error(fatal, 'apply_topography_edits_from_file: '//&
278  'Failed to close '//trim(topo_edits_file))
279 
280  do n = 1, n_edits
281  i = ig(n) - g%isd_global + 2 ! +1 for python indexing and +1 for ig-isd_global+1
282  j = jg(n) - g%jsd_global + 2
283  if (i>=g%isc .and. i<=g%iec .and. j>=g%jsc .and. j<=g%jec) then
284  if (new_depth(n)/=0.) then
285  write(*,'(a,3i5,f8.2,a,f8.2,2i4)') &
286  'Ocean topography edit: ',n,ig(n),jg(n),d(i,j)/m_to_z,'->',abs(new_depth(n)),i,j
287  d(i,j) = abs(m_to_z*new_depth(n)) ! Allows for height-file edits (i.e. converts negatives)
288  else
289  call mom_error(fatal, ' apply_topography_edits_from_file: '//&
290  "A zero depth edit would change the land mask and is not allowed in"//trim(topo_edits_file))
291  endif
292  endif
293  enddo
294 
295  deallocate( ig, jg, new_depth )
296 
297  call calltree_leave(trim(mdl)//'()')

◆ compute_global_grid_integrals()

subroutine, public mom_shared_initialization::compute_global_grid_integrals ( type(dyn_horgrid_type), intent(inout)  G,
type(unit_scale_type), intent(in), optional  US 
)

Pre-compute global integrals of grid quantities (like masked ocean area) for later use in reporting diagnostics.

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

Definition at line 1149 of file MOM_shared_initialization.F90.

1149  type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid
1150  type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type
1151 
1152  ! Local variables
1153  real, dimension(G%isc:G%iec, G%jsc:G%jec) :: tmpForSumming
1154  real :: area_scale ! A scaling factor for area into MKS units
1155  integer :: i,j
1156 
1157  area_scale = 1.0 ; if (present(us)) area_scale = us%L_to_m**2
1158 
1159  tmpforsumming(:,:) = 0.
1160  g%areaT_global = 0.0 ; g%IareaT_global = 0.0
1161  do j=g%jsc,g%jec ; do i=g%isc,g%iec
1162  tmpforsumming(i,j) = area_scale*g%areaT(i,j) * g%mask2dT(i,j)
1163  enddo ; enddo
1164  g%areaT_global = reproducing_sum(tmpforsumming)
1165 
1166  if (g%areaT_global == 0.0) &
1167  call mom_error(fatal, "compute_global_grid_integrals: "//&
1168  "zero ocean area (check topography?)")
1169 
1170  g%IareaT_global = 1.0 / (g%areaT_global)
mom_coms::reproducing_sum
Find an accurate and order-invariant sum of a distributed 2d or 3d field.
Definition: MOM_coms.F90:53

◆ diagnosemaximumdepth()

real function, public mom_shared_initialization::diagnosemaximumdepth ( real, dimension(g%isd:g%ied,g%jsd:g%jed), intent(in)  D,
type(dyn_horgrid_type), intent(in)  G 
)

Return the global maximum ocean bottom depth in the same units as the input depth.

Parameters
[in]gThe dynamic horizontal grid type
[in]dOcean bottom depth in m or Z
Returns
The global maximum ocean bottom depth in m or Z

Definition at line 125 of file MOM_shared_initialization.F90.

125  type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type
126  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
127  intent(in) :: D !< Ocean bottom depth in m or Z
128  real :: diagnoseMaximumDepth !< The global maximum ocean bottom depth in m or Z
129  ! Local variables
130  integer :: i,j
131  diagnosemaximumdepth = d(g%isc,g%jsc)
132  do j=g%jsc, g%jec ; do i=g%isc, g%iec
133  diagnosemaximumdepth = max(diagnosemaximumdepth,d(i,j))
134  enddo ; enddo
135  call max_across_pes(diagnosemaximumdepth)

◆ initialize_grid_rotation_angle()

subroutine, public mom_shared_initialization::initialize_grid_rotation_angle ( type(dyn_horgrid_type), intent(inout)  G,
type(param_file_type), intent(in)  PF 
)

initialize_grid_rotation_angle initializes the arrays with the sine and cosine of the angle between logical north on the grid and true north.

Parameters
[in,out]gThe dynamic horizontal grid
[in]pfA structure indicating the open file to parse for model parameter values.

Definition at line 548 of file MOM_shared_initialization.F90.

548  type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid
549  type(param_file_type), intent(in) :: PF !< A structure indicating the open file
550  !! to parse for model parameter values.
551 
552  real :: angle, lon_scale
553  real :: len_lon ! The periodic range of longitudes, usually 360 degrees.
554  real :: pi_720deg ! One quarter the conversion factor from degrees to radians.
555  real :: lonB(2,2) ! The longitude of a point, shifted to have about the same value.
556  character(len=40) :: mdl = "initialize_grid_rotation_angle" ! This subroutine's name.
557  logical :: use_bugs
558  integer :: i, j, m, n
559 
560  call get_param(pf, mdl, "GRID_ROTATION_ANGLE_BUGS", use_bugs, &
561  "If true, use an older algorithm to calculate the sine and "//&
562  "cosines needed rotate between grid-oriented directions and "//&
563  "true north and east. Differences arise at the tripolar fold.", &
564  default=.false.)
565 
566  if (use_bugs) then
567  do j=g%jsc,g%jec ; do i=g%isc,g%iec
568  lon_scale = cos((g%geoLatBu(i-1,j-1) + g%geoLatBu(i,j-1 ) + &
569  g%geoLatBu(i-1,j) + g%geoLatBu(i,j)) * atan(1.0)/180)
570  angle = atan2((g%geoLonBu(i-1,j) + g%geoLonBu(i,j) - &
571  g%geoLonBu(i-1,j-1) - g%geoLonBu(i,j-1))*lon_scale, &
572  g%geoLatBu(i-1,j) + g%geoLatBu(i,j) - &
573  g%geoLatBu(i-1,j-1) - g%geoLatBu(i,j-1) )
574  g%sin_rot(i,j) = sin(angle) ! angle is the clockwise angle from lat/lon to ocean
575  g%cos_rot(i,j) = cos(angle) ! grid (e.g. angle of ocean "north" from true north)
576  enddo ; enddo
577 
578  ! This is not right at a tripolar or cubed-sphere fold.
579  call pass_var(g%cos_rot, g%Domain)
580  call pass_var(g%sin_rot, g%Domain)
581  else
582  pi_720deg = atan(1.0) / 180.0
583  len_lon = 360.0 ; if (g%len_lon > 0.0) len_lon = g%len_lon
584  do j=g%jsc,g%jec ; do i=g%isc,g%iec
585  do n=1,2 ; do m=1,2
586  lonb(m,n) = modulo_around_point(g%geoLonBu(i+m-2,j+n-2), g%geoLonT(i,j), len_lon)
587  enddo ; enddo
588  lon_scale = cos(pi_720deg*((g%geoLatBu(i-1,j-1) + g%geoLatBu(i,j)) + &
589  (g%geoLatBu(i,j-1) + g%geoLatBu(i-1,j)) ) )
590  angle = atan2(lon_scale*((lonb(1,2) - lonb(2,1)) + (lonb(2,2) - lonb(1,1))), &
591  (g%geoLatBu(i-1,j) - g%geoLatBu(i,j-1)) + &
592  (g%geoLatBu(i,j) - g%geoLatBu(i-1,j-1)) )
593  g%sin_rot(i,j) = sin(angle) ! angle is the clockwise angle from lat/lon to ocean
594  g%cos_rot(i,j) = cos(angle) ! grid (e.g. angle of ocean "north" from true north)
595  enddo ; enddo
596 
597  call pass_vector(g%cos_rot, g%sin_rot, g%Domain, stagger=agrid)
598  endif
599 

◆ initialize_topography_from_file()

subroutine, public mom_shared_initialization::initialize_topography_from_file ( real, dimension(g%isd:g%ied,g%jsd:g%jed), intent(out)  D,
type(dyn_horgrid_type), intent(in)  G,
type(param_file_type), intent(in)  param_file,
type(unit_scale_type), intent(in), optional  US 
)

Read gridded depths from file.

Parameters
[in]gThe dynamic horizontal grid type
[out]dOcean bottom depth in m or Z if US is present
[in]param_fileParameter file structure
[in]usA dimensional unit scaling type

Definition at line 141 of file MOM_shared_initialization.F90.

141  type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type
142  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
143  intent(out) :: D !< Ocean bottom depth in m or Z if US is present
144  type(param_file_type), intent(in) :: param_file !< Parameter file structure
145  type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type
146  ! Local variables
147  real :: m_to_Z ! A dimensional rescaling factor.
148  character(len=200) :: filename, topo_file, inputdir ! Strings for file/path
149  character(len=200) :: topo_varname ! Variable name in file
150  character(len=40) :: mdl = "initialize_topography_from_file" ! This subroutine's name.
151 
152  call calltree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")
153 
154  m_to_z = 1.0 ; if (present(us)) m_to_z = us%m_to_Z
155 
156  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
157  inputdir = slasher(inputdir)
158  call get_param(param_file, mdl, "TOPO_FILE", topo_file, &
159  "The file from which the bathymetry is read.", &
160  default="topog.nc")
161  call get_param(param_file, mdl, "TOPO_VARNAME", topo_varname, &
162  "The name of the bathymetry variable in TOPO_FILE.", &
163  default="depth")
164 
165  filename = trim(inputdir)//trim(topo_file)
166  call log_param(param_file, mdl, "INPUTDIR/TOPO_FILE", filename)
167 
168  if (.not.file_exists(filename, g%Domain)) call mom_error(fatal, &
169  " initialize_topography_from_file: Unable to open "//trim(filename))
170 
171  d(:,:) = -9.e30*m_to_z ! Initializing to a very large negative depth (tall mountains) everywhere
172  ! before reading from a file should do nothing. However, in the instance of
173  ! masked-out PEs, halo regions are not updated when a processor does not
174  ! exist. We need to ensure the depth in masked-out PEs appears to be that
175  ! of land so this line does that in the halo regions. For non-masked PEs
176  ! the halo region is filled properly with a later pass_var().
177  call mom_read_data(filename, trim(topo_varname), d, g%Domain, scale=m_to_z)
178 
179  call apply_topography_edits_from_file(d, g, param_file, us)
180 
181  call calltree_leave(trim(mdl)//'()')

◆ initialize_topography_named()

subroutine, public mom_shared_initialization::initialize_topography_named ( real, dimension(g%isd:g%ied,g%jsd:g%jed), intent(out)  D,
type(dyn_horgrid_type), intent(in)  G,
type(param_file_type), intent(in)  param_file,
character(len=*), intent(in)  topog_config,
real, intent(in)  max_depth,
type(unit_scale_type), intent(in), optional  US 
)

initialize the bathymetry based on one of several named idealized configurations

Parameters
[in]gThe dynamic horizontal grid type
[out]dOcean bottom depth in m or Z if US is present
[in]param_fileParameter file structure
[in]topog_configThe name of an idealized topographic configuration
[in]max_depthMaximum depth of model in the units of D
[in]usA dimensional unit scaling type

Definition at line 302 of file MOM_shared_initialization.F90.

302  type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type
303  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
304  intent(out) :: D !< Ocean bottom depth in m or Z if US is present
305  type(param_file_type), intent(in) :: param_file !< Parameter file structure
306  character(len=*), intent(in) :: topog_config !< The name of an idealized
307  !! topographic configuration
308  real, intent(in) :: max_depth !< Maximum depth of model in the units of D
309  type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type
310 
311  ! This subroutine places the bottom depth in m into D(:,:), shaped according to the named config.
312 
313  ! Local variables
314  real :: m_to_Z ! A dimensional rescaling factor.
315  real :: min_depth ! The minimum depth [Z ~> m].
316  real :: PI ! 3.1415926... calculated as 4*atan(1)
317  real :: D0 ! A constant to make the maximum basin depth MAXIMUM_DEPTH.
318  real :: expdecay ! A decay scale of associated with the sloping boundaries [m].
319  real :: Dedge ! The depth [Z ~> m], at the basin edge
320 ! real :: south_lat, west_lon, len_lon, len_lat, Rad_earth
321  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
322  character(len=40) :: mdl = "initialize_topography_named" ! This subroutine's name.
323  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
324  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
325 
326  call calltree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")
327  call mom_mesg(" MOM_shared_initialization.F90, initialize_topography_named: "//&
328  "TOPO_CONFIG = "//trim(topog_config), 5)
329 
330  m_to_z = 1.0 ; if (present(us)) m_to_z = us%m_to_Z
331 
332  call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, &
333  "The minimum depth of the ocean.", units="m", default=0.0, scale=m_to_z)
334  if (max_depth<=0.) call mom_error(fatal,"initialize_topography_named: "// &
335  "MAXIMUM_DEPTH has a non-sensical value! Was it set?")
336 
337  if (trim(topog_config) /= "flat") then
338  call get_param(param_file, mdl, "EDGE_DEPTH", dedge, &
339  "The depth at the edge of one of the named topographies.", &
340  units="m", default=100.0, scale=m_to_z)
341 ! call get_param(param_file, mdl, "SOUTHLAT", south_lat, &
342 ! "The southern latitude of the domain.", units="degrees", &
343 ! fail_if_missing=.true.)
344 ! call get_param(param_file, mdl, "LENLAT", len_lat, &
345 ! "The latitudinal length of the domain.", units="degrees", &
346 ! fail_if_missing=.true.)
347 ! call get_param(param_file, mdl, "WESTLON", west_lon, &
348 ! "The western longitude of the domain.", units="degrees", &
349 ! default=0.0)
350 ! call get_param(param_file, mdl, "LENLON", len_lon, &
351 ! "The longitudinal length of the domain.", units="degrees", &
352 ! fail_if_missing=.true.)
353 ! call get_param(param_file, mdl, "RAD_EARTH", Rad_Earth, &
354 ! "The radius of the Earth.", units="m", default=6.378e6)
355  call get_param(param_file, mdl, "TOPOG_SLOPE_SCALE", expdecay, &
356  "The exponential decay scale used in defining some of "//&
357  "the named topographies.", units="m", default=400000.0)
358  endif
359 
360 
361  pi = 4.0*atan(1.0)
362 
363  if (trim(topog_config) == "flat") then
364  do i=is,ie ; do j=js,je ; d(i,j) = max_depth ; enddo ; enddo
365  elseif (trim(topog_config) == "spoon") then
366  d0 = (max_depth - dedge) / &
367  ((1.0 - exp(-0.5*g%len_lat*g%Rad_earth*pi/(180.0 *expdecay))) * &
368  (1.0 - exp(-0.5*g%len_lat*g%Rad_earth*pi/(180.0 *expdecay))))
369  do i=is,ie ; do j=js,je
370  ! This sets a bowl shaped (sort of) bottom topography, with a !
371  ! maximum depth of max_depth. !
372  d(i,j) = dedge + d0 * &
373  (sin(pi * (g%geoLonT(i,j) - (g%west_lon)) / g%len_lon) * &
374  (1.0 - exp((g%geoLatT(i,j) - (g%south_lat+g%len_lat))*g%Rad_earth*pi / &
375  (180.0*expdecay)) ))
376  enddo ; enddo
377  elseif (trim(topog_config) == "bowl") then
378  d0 = (max_depth - dedge) / &
379  ((1.0 - exp(-0.5*g%len_lat*g%Rad_earth*pi/(180.0 *expdecay))) * &
380  (1.0 - exp(-0.5*g%len_lat*g%Rad_earth*pi/(180.0 *expdecay))))
381 
382  ! This sets a bowl shaped (sort of) bottom topography, with a
383  ! maximum depth of max_depth.
384  do i=is,ie ; do j=js,je
385  d(i,j) = dedge + d0 * &
386  (sin(pi * (g%geoLonT(i,j) - g%west_lon) / g%len_lon) * &
387  ((1.0 - exp(-(g%geoLatT(i,j) - g%south_lat)*g%Rad_Earth*pi/ &
388  (180.0*expdecay))) * &
389  (1.0 - exp((g%geoLatT(i,j) - (g%south_lat+g%len_lat))* &
390  g%Rad_Earth*pi/(180.0*expdecay)))))
391  enddo ; enddo
392  elseif (trim(topog_config) == "halfpipe") then
393  d0 = max_depth - dedge
394  do i=is,ie ; do j=js,je
395  d(i,j) = dedge + d0 * abs(sin(pi*(g%geoLatT(i,j) - g%south_lat)/g%len_lat))
396  enddo ; enddo
397  else
398  call mom_error(fatal,"initialize_topography_named: "// &
399  "Unrecognized topography name "//trim(topog_config))
400  endif
401 
402  ! This is here just for safety. Hopefully it doesn't do anything.
403  do i=is,ie ; do j=js,je
404  if (d(i,j) > max_depth) d(i,j) = max_depth
405  if (d(i,j) < min_depth) d(i,j) = 0.5*min_depth
406  enddo ; enddo
407 
408  call calltree_leave(trim(mdl)//'()')

◆ limit_topography()

subroutine, public mom_shared_initialization::limit_topography ( real, dimension(g%isd:g%ied,g%jsd:g%jed), intent(inout)  D,
type(dyn_horgrid_type), intent(in)  G,
type(param_file_type), intent(in)  param_file,
real, intent(in)  max_depth,
type(unit_scale_type), intent(in), optional  US 
)

limit_topography ensures that min_depth < D(x,y) < max_depth

Parameters
[in]gThe dynamic horizontal grid type
[in,out]dOcean bottom depth in m or Z if US is present
[in]param_fileParameter file structure
[in]max_depthMaximum depth of model in the units of D
[in]usA dimensional unit scaling type

Definition at line 415 of file MOM_shared_initialization.F90.

415  type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type
416  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
417  intent(inout) :: D !< Ocean bottom depth in m or Z if US is present
418  type(param_file_type), intent(in) :: param_file !< Parameter file structure
419  real, intent(in) :: max_depth !< Maximum depth of model in the units of D
420  type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type
421 
422  ! Local variables
423  real :: m_to_Z ! A dimensional rescaling factor.
424  integer :: i, j
425  character(len=40) :: mdl = "limit_topography" ! This subroutine's name.
426  real :: min_depth, mask_depth
427 
428  call calltree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")
429 
430  m_to_z = 1.0 ; if (present(us)) m_to_z = us%m_to_Z
431 
432  call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, &
433  "If MASKING_DEPTH is unspecified, then anything shallower than "//&
434  "MINIMUM_DEPTH is assumed to be land and all fluxes are masked out. "//&
435  "If MASKING_DEPTH is specified, then all depths shallower than "//&
436  "MINIMUM_DEPTH but deeper than MASKING_DEPTH are rounded to MINIMUM_DEPTH.", &
437  units="m", default=0.0, scale=m_to_z)
438  call get_param(param_file, mdl, "MASKING_DEPTH", mask_depth, &
439  "The depth below which to mask the ocean as land.", &
440  units="m", default=-9999.0, scale=m_to_z, do_not_log=.true.)
441 
442 ! Make sure that min_depth < D(x,y) < max_depth
443  if (mask_depth < -9990.*m_to_z) then
444  do j=g%jsd,g%jed ; do i=g%isd,g%ied
445  d(i,j) = min( max( d(i,j), 0.5*min_depth ), max_depth )
446  enddo ; enddo
447  else
448  do j=g%jsd,g%jed ; do i=g%isd,g%ied
449  if (d(i,j)>0.) then
450  d(i,j) = min( max( d(i,j), min_depth ), max_depth )
451  else
452  d(i,j) = 0.
453  endif
454  enddo ; enddo
455  endif
456 
457  call calltree_leave(trim(mdl)//'()')

◆ modulo_around_point()

real function mom_shared_initialization::modulo_around_point ( real, intent(in)  x,
real, intent(in)  xc,
real, intent(in)  Lx 
)
private

Return the modulo value of x in an interval [xc-(Lx/2) xc+(Lx/2)] If Lx<=0, then it returns x without applying modulo arithmetic.

Parameters
[in]xValue to which to apply modulo arithmetic
[in]xcCenter of modulo range
[in]lxModulo range width
Returns
x shifted by an integer multiple of Lx to be close to xc.

Definition at line 606 of file MOM_shared_initialization.F90.

606  real, intent(in) :: x !< Value to which to apply modulo arithmetic
607  real, intent(in) :: xc !< Center of modulo range
608  real, intent(in) :: Lx !< Modulo range width
609  real :: x_mod !< x shifted by an integer multiple of Lx to be close to xc.
610 
611  if (lx > 0.0) then
612  x_mod = modulo(x - (xc - 0.5*lx), lx) + (xc - 0.5*lx)
613  else
614  x_mod = x
615  endif

◆ mom_calculate_grad_coriolis()

subroutine, public mom_shared_initialization::mom_calculate_grad_coriolis ( real, dimension(g%isd:g%ied,g%jsd:g%jed), intent(out)  dF_dx,
real, dimension(g%isd:g%ied,g%jsd:g%jed), intent(out)  dF_dy,
type(dyn_horgrid_type), intent(inout)  G,
type(unit_scale_type), intent(in), optional  US 
)

Calculates the components of grad f (Coriolis parameter)

Parameters
[in,out]gThe dynamic horizontal grid type
[out]df_dxx-component of grad f [T-1 L-1 ~> s-1 m-1]
[out]df_dyy-component of grad f [T-1 L-1 ~> s-1 m-1]
[in]usA dimensional unit scaling type

Definition at line 91 of file MOM_shared_initialization.F90.

91  type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid type
92  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
93  intent(out) :: dF_dx !< x-component of grad f [T-1 L-1 ~> s-1 m-1]
94  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
95  intent(out) :: dF_dy !< y-component of grad f [T-1 L-1 ~> s-1 m-1]
96  type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type
97  ! Local variables
98  integer :: i,j
99  real :: m_to_L ! A unit conversion factor [L m-1 ~> nondim]
100  real :: f1, f2
101 
102  m_to_l = 1.0 ; if (present(us)) m_to_l = us%m_to_L
103 
104  if ((lbound(g%CoriolisBu,1) > g%isc-1) .or. &
105  (lbound(g%CoriolisBu,2) > g%isc-1)) then
106  ! The gradient of the Coriolis parameter can not be calculated with this grid.
107  df_dx(:,:) = 0.0 ; df_dy(:,:) = 0.0
108  return
109  endif
110 
111  do j=g%jsc, g%jec ; do i=g%isc, g%iec
112  f1 = 0.5*( g%CoriolisBu(i,j) + g%CoriolisBu(i,j-1) )
113  f2 = 0.5*( g%CoriolisBu(i-1,j) + g%CoriolisBu(i-1,j-1) )
114  df_dx(i,j) = g%IdxT(i,j) * ( f1 - f2 )
115  f1 = 0.5*( g%CoriolisBu(i,j) + g%CoriolisBu(i-1,j) )
116  f2 = 0.5*( g%CoriolisBu(i,j-1) + g%CoriolisBu(i-1,j-1) )
117  df_dy(i,j) = g%IdyT(i,j) * ( f1 - f2 )
118  enddo ; enddo
119  call pass_vector(df_dx, df_dy, g%Domain, stagger=agrid)
120 

◆ mom_initialize_rotation()

subroutine, public mom_shared_initialization::mom_initialize_rotation ( real, dimension(g%isdb:g%iedb,g%jsdb:g%jedb), intent(out)  f,
type(dyn_horgrid_type), intent(in)  G,
type(param_file_type), intent(in)  PF,
type(unit_scale_type), intent(in), optional  US 
)

MOM_initialize_rotation makes the appropriate call to set up the Coriolis parameter.

Parameters
[in]gThe dynamic horizontal grid type
[out]fThe Coriolis parameter [T-1 ~> s-1]
[in]pfParameter file structure
[in]usA dimensional unit scaling type

Definition at line 58 of file MOM_shared_initialization.F90.

58  type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type
59  real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), intent(out) :: f !< The Coriolis parameter [T-1 ~> s-1]
60  type(param_file_type), intent(in) :: PF !< Parameter file structure
61  type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type
62 
63 ! This subroutine makes the appropriate call to set up the Coriolis parameter.
64 ! This is a separate subroutine so that it can be made public and shared with
65 ! the ice-sheet code or other components.
66 ! Set up the Coriolis parameter, f, either analytically or from file.
67  character(len=40) :: mdl = "MOM_initialize_rotation" ! This subroutine's name.
68  character(len=200) :: config
69 
70  call calltree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")
71  call get_param(pf, mdl, "ROTATION", config, &
72  "This specifies how the Coriolis parameter is specified: \n"//&
73  " \t 2omegasinlat - Use twice the planetary rotation rate \n"//&
74  " \t\t times the sine of latitude.\n"//&
75  " \t betaplane - Use a beta-plane or f-plane.\n"//&
76  " \t USER - call a user modified routine.", &
77  default="2omegasinlat")
78  select case (trim(config))
79  case ("2omegasinlat"); call set_rotation_planetary(f, g, pf, us)
80  case ("beta"); call set_rotation_beta_plane(f, g, pf, us)
81  case ("betaplane"); call set_rotation_beta_plane(f, g, pf, us)
82  !case ("nonrotating") ! Note from AJA: Missing case?
83  case default ; call mom_error(fatal,"MOM_initialize: "// &
84  "Unrecognized rotation setup "//trim(config))
85  end select
86  call calltree_leave(trim(mdl)//'()')

◆ mom_shared_init_init()

subroutine, public mom_shared_initialization::mom_shared_init_init ( type(param_file_type), intent(in)  PF)

MOM_shared_init_init just writes the code version.

Parameters
[in]pfA structure indicating the open file to parse for model parameter values.

Definition at line 43 of file MOM_shared_initialization.F90.

43  type(param_file_type), intent(in) :: PF !< A structure indicating the open file
44  !! to parse for model parameter values.
45 
46  character(len=40) :: mdl = "MOM_shared_initialization" ! This module's name.
47 
48 ! This include declares and sets the variable "version".
49 #include "version_variable.h"
50  call log_version(pf, mdl, version, &
51  "Sharable code to initialize time-invariant fields, like bathymetry and Coriolis parameters.")
52 

◆ read_face_length_list()

subroutine, public mom_shared_initialization::read_face_length_list ( integer, intent(in)  iounit,
character(len=*), intent(in)  filename,
integer, intent(out)  num_lines,
character(len=120), dimension(:), pointer  lines 
)

This subroutine reads and counts the non-blank lines in the face length list file, after removing comments.

Parameters
[in]iounitAn open I/O unit number for the file
[in]filenameThe name of the face-length file to read
[out]num_linesThe number of non-blank lines in the file
linesThe non-blank lines, after removing comments

Definition at line 1055 of file MOM_shared_initialization.F90.

1055  integer, intent(in) :: iounit !< An open I/O unit number for the file
1056  character(len=*), intent(in) :: filename !< The name of the face-length file to read
1057  integer, intent(out) :: num_lines !< The number of non-blank lines in the file
1058  character(len=120), dimension(:), pointer :: lines !< The non-blank lines, after removing comments
1059 
1060  ! This subroutine reads and counts the non-blank lines in the face length
1061  ! list file, after removing comments.
1062  character(len=120) :: line, line_up
1063  logical :: found_u, found_v
1064  integer :: isu, isv, icom, verbose
1065  integer :: last
1066 
1067  num_lines = 0
1068 
1069  if (iounit <= 0) return
1070  rewind(iounit)
1071  do while(.true.)
1072  read(iounit, '(a)', end=8, err=9) line
1073  last = len_trim(line)
1074  ! Eliminate either F90 or C comments from the line.
1075  icom = index(line(:last), "!") ; if (icom > 0) last = icom-1
1076  icom = index(line(:last), "/*") ; if (icom > 0) last = icom-1
1077  if (last < 1) cycle
1078 
1079  ! Detect keywords
1080  line_up = uppercase(line)
1081  found_u = .false.; found_v = .false.
1082  isu = index(line_up(:last), "U_WIDTH" ); if (isu > 0) found_u = .true.
1083  isv = index(line_up(:last), "V_WIDTH" ); if (isv > 0) found_v = .true.
1084 
1085  if (found_u .and. found_v) call mom_error(fatal, &
1086  "read_face_length_list : both U_WIDTH and V_WIDTH found when "//&
1087  "reading the line "//trim(line(:last))//" in file "//trim(filename))
1088  if (found_u .or. found_v) then
1089  num_lines = num_lines + 1
1090  if (associated(lines)) then
1091  lines(num_lines) = line(1:last)
1092  endif
1093  endif
1094  enddo ! while (.true.)
1095 
1096 8 continue
1097  return
1098 
1099 9 call mom_error(fatal, "read_face_length_list : "//&
1100  "Error while reading file "//trim(filename))
1101 

◆ reset_face_lengths_file()

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

This subroutine sets the open face lengths at selected points to restrict passages to their observed widths from a arrays read from a file.

Parameters
[in,out]gThe dynamic horizontal grid
[in]param_fileA structure to parse for run-time parameters
[in]usA dimensional unit scaling type

Definition at line 750 of file MOM_shared_initialization.F90.

750  type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid
751  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
752  type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type
753 
754  ! Local variables
755  character(len=40) :: mdl = "reset_face_lengths_file" ! This subroutine's name.
756  character(len=256) :: mesg ! Message for error messages.
757  character(len=200) :: filename, chan_file, inputdir ! Strings for file/path
758  real :: m_to_L ! A unit conversion factor [L m-1 ~> nondim]
759  real :: L_to_m ! A unit conversion factor [m L-1 ~> nondim]
760  integer :: i, j, isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
761  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
762  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
763  ! These checks apply regardless of the chosen option.
764 
765  call calltree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")
766  m_to_l = 1.0 ; if (present(us)) m_to_l = us%m_to_L
767  l_to_m = 1.0 ; if (present(us)) l_to_m = us%L_to_m
768 
769  call get_param(param_file, mdl, "CHANNEL_WIDTH_FILE", chan_file, &
770  "The file from which the list of narrowed channels is read.", &
771  default="ocean_geometry.nc")
772  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
773  inputdir = slasher(inputdir)
774  filename = trim(inputdir)//trim(chan_file)
775  call log_param(param_file, mdl, "INPUTDIR/CHANNEL_WIDTH_FILE", filename)
776 
777  if (is_root_pe()) then ; if (.not.file_exists(filename)) &
778  call mom_error(fatal," reset_face_lengths_file: Unable to open "//&
779  trim(filename))
780  endif
781 
782  call mom_read_vector(filename, "dyCuo", "dxCvo", g%dy_Cu, g%dx_Cv, g%Domain, scale=m_to_l)
783  call pass_vector(g%dy_Cu, g%dx_Cv, g%Domain, to_all+scalar_pair, cgrid_ne)
784 
785  do j=jsd,jed ; do i=isdb,iedb
786  if (l_to_m*g%dy_Cu(i,j) > l_to_m*g%dyCu(i,j)) then
787  write(mesg,'("dy_Cu of ",ES11.4," exceeds unrestricted width of ",ES11.4,&
788  &" by ",ES11.4," at lon/lat of ", ES11.4, ES11.4)') &
789  l_to_m*g%dy_Cu(i,j), l_to_m*g%dyCu(i,j), l_to_m*g%dy_Cu(i,j)-l_to_m*g%dyCu(i,j), &
790  g%geoLonCu(i,j), g%geoLatCu(i,j)
791  call mom_error(fatal,"reset_face_lengths_file "//mesg)
792  endif
793  g%areaCu(i,j) = g%dxCu(i,j) * g%dy_Cu(i,j)
794  g%IareaCu(i,j) = 0.0
795  if (g%areaCu(i,j) > 0.0) g%IareaCu(i,j) = g%mask2dCu(i,j) / (g%areaCu(i,j))
796  enddo ; enddo
797 
798  do j=jsdb,jedb ; do i=isd,ied
799  if (l_to_m*g%dx_Cv(i,j) > l_to_m*g%dxCv(i,j)) then
800  write(mesg,'("dx_Cv of ",ES11.4," exceeds unrestricted width of ",ES11.4,&
801  &" by ",ES11.4, " at lon/lat of ", ES11.4, ES11.4)') &
802  l_to_m*g%dx_Cv(i,j), l_to_m*g%dxCv(i,j), l_to_m*g%dx_Cv(i,j)-l_to_m*g%dxCv(i,j), &
803  g%geoLonCv(i,j), g%geoLatCv(i,j)
804 
805  call mom_error(fatal,"reset_face_lengths_file "//mesg)
806  endif
807  g%areaCv(i,j) = g%dyCv(i,j) * g%dx_Cv(i,j)
808  g%IareaCv(i,j) = 0.0
809  if (g%areaCv(i,j) > 0.0) g%IareaCv(i,j) = g%mask2dCv(i,j) / (g%areaCv(i,j))
810  enddo ; enddo
811 
812  call calltree_leave(trim(mdl)//'()')

◆ reset_face_lengths_list()

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

This subroutine sets the open face lengths at selected points to restrict passages to their observed widths from a list read from a file.

Parameters
[in,out]gThe dynamic horizontal grid
[in]param_fileA structure to parse for run-time parameters
[in]usA dimensional unit scaling type

Definition at line 820 of file MOM_shared_initialization.F90.

820  type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid
821  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
822  type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type
823 
824  ! Local variables
825  character(len=120), pointer, dimension(:) :: lines => null()
826  character(len=120) :: line
827  character(len=200) :: filename, chan_file, inputdir, mesg ! Strings for file/path
828  character(len=40) :: mdl = "reset_face_lengths_list" ! This subroutine's name.
829  real, pointer, dimension(:,:) :: &
830  u_lat => null(), u_lon => null(), v_lat => null(), v_lon => null()
831  real, pointer, dimension(:) :: &
832  u_width => null(), v_width => null()
833  real :: m_to_L ! A unit conversion factor [L m-1 ~> nondim]
834  real :: L_to_m ! A unit conversion factor [m L-1 ~> nondim]
835  real :: lat, lon ! The latitude and longitude of a point.
836  real :: len_lon ! The periodic range of longitudes, usually 360 degrees.
837  real :: len_lat ! The range of latitudes, usually 180 degrees.
838  real :: lon_p, lon_m ! The longitude of a point shifted by 360 degrees.
839  logical :: check_360 ! If true, check for longitudes that are shifted by
840  ! +/- 360 degrees from the specified range of values.
841  logical :: found_u, found_v
842  logical :: unit_in_use
843  integer :: ios, iounit, isu, isv
844  integer :: last, num_lines, nl_read, ln, npt, u_pt, v_pt
845  integer :: i, j, isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
846  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
847  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
848 
849  call calltree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")
850  m_to_l = 1.0 ; if (present(us)) m_to_l = us%m_to_L
851  l_to_m = 1.0 ; if (present(us)) l_to_m = us%L_to_m
852 
853  call get_param(param_file, mdl, "CHANNEL_LIST_FILE", chan_file, &
854  "The file from which the list of narrowed channels is read.", &
855  default="MOM_channel_list")
856  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
857  inputdir = slasher(inputdir)
858  filename = trim(inputdir)//trim(chan_file)
859  call log_param(param_file, mdl, "INPUTDIR/CHANNEL_LIST_FILE", filename)
860  call get_param(param_file, mdl, "CHANNEL_LIST_360_LON_CHECK", check_360, &
861  "If true, the channel configuration list works for any "//&
862  "longitudes in the range of -360 to 360.", default=.true.)
863 
864  if (is_root_pe()) then
865  ! Open the input file.
866  if (.not.file_exists(filename)) call mom_error(fatal, &
867  " reset_face_lengths_list: Unable to open "//trim(filename))
868 
869  ! Find an unused unit number.
870  do iounit=10,512
871  INQUIRE(iounit,opened=unit_in_use) ; if (.not.unit_in_use) exit
872  enddo
873  if (iounit >= 512) call mom_error(fatal, &
874  "reset_face_lengths_list: No unused file unit could be found.")
875 
876  ! Open the parameter file.
877  open(iounit, file=trim(filename), access='SEQUENTIAL', &
878  form='FORMATTED', action='READ', position='REWIND', iostat=ios)
879  if (ios /= 0) call mom_error(fatal, &
880  "reset_face_lengths_list: Error opening "//trim(filename))
881 
882  ! Count the number of u_width and v_width entries.
883  call read_face_length_list(iounit, filename, num_lines, lines)
884  endif
885 
886  len_lon = 360.0 ; if (g%len_lon > 0.0) len_lon = g%len_lon
887  len_lat = 180.0 ; if (g%len_lat > 0.0) len_lat = g%len_lat
888  ! Broadcast the number of lines and allocate the required space.
889  call broadcast(num_lines, root_pe())
890  u_pt = 0 ; v_pt = 0
891  if (num_lines > 0) then
892  allocate (lines(num_lines))
893  if (num_lines > 0) then
894  allocate(u_lat(2,num_lines)) ; u_lat(:,:) = -1e34
895  allocate(u_lon(2,num_lines)) ; u_lon(:,:) = -1e34
896  allocate(u_width(num_lines)) ; u_width(:) = -1e34
897 
898  allocate(v_lat(2,num_lines)) ; v_lat(:,:) = -1e34
899  allocate(v_lon(2,num_lines)) ; v_lon(:,:) = -1e34
900  allocate(v_width(num_lines)) ; v_width(:) = -1e34
901  endif
902 
903  ! Actually read the lines.
904  if (is_root_pe()) then
905  call read_face_length_list(iounit, filename, nl_read, lines)
906  if (nl_read /= num_lines) &
907  call mom_error(fatal, 'reset_face_lengths_list : Found different '// &
908  'number of valid lines on second reading of '//trim(filename))
909  close(iounit) ; iounit = -1
910  endif
911 
912  ! Broadcast the lines.
913  call broadcast(lines, 120, root_pe())
914 
915  ! Populate the u_width, etc., data.
916  do ln=1,num_lines
917  line = lines(ln)
918  ! Detect keywords
919  found_u = .false.; found_v = .false.
920  isu = index(uppercase(line), "U_WIDTH" ); if (isu > 0) found_u = .true.
921  isv = index(uppercase(line), "V_WIDTH" ); if (isv > 0) found_v = .true.
922 
923  ! Store and check the relevant values.
924  if (found_u) then
925  u_pt = u_pt + 1
926  read(line(isu+8:),*) u_lon(1:2,u_pt), u_lat(1:2,u_pt), u_width(u_pt)
927  if (is_root_pe()) then
928  if (check_360) then
929  if ((abs(u_lon(1,u_pt)) > len_lon) .or. (abs(u_lon(2,u_pt)) > len_lon)) &
930  call mom_error(warning, "reset_face_lengths_list : Out-of-bounds "//&
931  "u-longitude found when reading line "//trim(line)//" from file "//&
932  trim(filename))
933  if ((abs(u_lat(1,u_pt)) > len_lat) .or. (abs(u_lat(2,u_pt)) > len_lat)) &
934  call mom_error(warning, "reset_face_lengths_list : Out-of-bounds "//&
935  "u-latitude found when reading line "//trim(line)//" from file "//&
936  trim(filename))
937  endif
938  if (u_lat(1,u_pt) > u_lat(2,u_pt)) &
939  call mom_error(warning, "reset_face_lengths_list : Out-of-order "//&
940  "u-face latitudes found when reading line "//trim(line)//" from file "//&
941  trim(filename))
942  if (u_lon(1,u_pt) > u_lon(2,u_pt)) &
943  call mom_error(warning, "reset_face_lengths_list : Out-of-order "//&
944  "u-face longitudes found when reading line "//trim(line)//" from file "//&
945  trim(filename))
946  if (u_width(u_pt) < 0.0) &
947  call mom_error(warning, "reset_face_lengths_list : Negative "//&
948  "u-width found when reading line "//trim(line)//" from file "//&
949  trim(filename))
950  endif
951  elseif (found_v) then
952  v_pt = v_pt + 1
953  read(line(isv+8:),*) v_lon(1:2,v_pt), v_lat(1:2,v_pt), v_width(v_pt)
954  if (is_root_pe()) then
955  if (check_360) then
956  if ((abs(v_lon(1,v_pt)) > len_lon) .or. (abs(v_lon(2,v_pt)) > len_lon)) &
957  call mom_error(warning, "reset_face_lengths_list : Out-of-bounds "//&
958  "v-longitude found when reading line "//trim(line)//" from file "//&
959  trim(filename))
960  if ((abs(v_lat(1,v_pt)) > len_lat) .or. (abs(v_lat(2,v_pt)) > len_lat)) &
961  call mom_error(warning, "reset_face_lengths_list : Out-of-bounds "//&
962  "v-latitude found when reading line "//trim(line)//" from file "//&
963  trim(filename))
964  endif
965  if (v_lat(1,v_pt) > v_lat(2,v_pt)) &
966  call mom_error(warning, "reset_face_lengths_list : Out-of-order "//&
967  "v-face latitudes found when reading line "//trim(line)//" from file "//&
968  trim(filename))
969  if (v_lon(1,v_pt) > v_lon(2,v_pt)) &
970  call mom_error(warning, "reset_face_lengths_list : Out-of-order "//&
971  "v-face longitudes found when reading line "//trim(line)//" from file "//&
972  trim(filename))
973  if (v_width(v_pt) < 0.0) &
974  call mom_error(warning, "reset_face_lengths_list : Negative "//&
975  "v-width found when reading line "//trim(line)//" from file "//&
976  trim(filename))
977  endif
978  endif
979  enddo
980 
981  deallocate(lines)
982  endif
983 
984  do j=jsd,jed ; do i=isdb,iedb
985  lat = g%geoLatCu(i,j) ; lon = g%geoLonCu(i,j)
986  if (check_360) then ; lon_p = lon+len_lon ; lon_m = lon-len_lon
987  else ; lon_p = lon ; lon_m = lon ; endif
988 
989  do npt=1,u_pt
990  if (((lat >= u_lat(1,npt)) .and. (lat <= u_lat(2,npt))) .and. &
991  (((lon >= u_lon(1,npt)) .and. (lon <= u_lon(2,npt))) .or. &
992  ((lon_p >= u_lon(1,npt)) .and. (lon_p <= u_lon(2,npt))) .or. &
993  ((lon_m >= u_lon(1,npt)) .and. (lon_m <= u_lon(2,npt)))) ) then
994 
995  g%dy_Cu(i,j) = g%mask2dCu(i,j) * m_to_l*min(l_to_m*g%dyCu(i,j), max(u_width(npt), 0.0))
996  if (j>=g%jsc .and. j<=g%jec .and. i>=g%isc .and. i<=g%iec) then ! Limit messages/checking to compute domain
997  if ( g%mask2dCu(i,j) == 0.0 ) then
998  write(*,'(A,2F8.2,A,4F8.2,A)') "read_face_lengths_list : G%mask2dCu=0 at ",lat,lon," (",&
999  u_lat(1,npt), u_lat(2,npt), u_lon(1,npt), u_lon(2,npt),") so grid metric is unmodified."
1000  else
1001  write(*,'(A,2F8.2,A,4F8.2,A5,F9.2,A1)') &
1002  "read_face_lengths_list : Modifying dy_Cu gridpoint at ",lat,lon," (",&
1003  u_lat(1,npt), u_lat(2,npt), u_lon(1,npt), u_lon(2,npt),") to ",l_to_m*g%dy_Cu(i,j),"m"
1004  endif
1005  endif
1006  endif
1007  enddo
1008 
1009  g%areaCu(i,j) = g%dxCu(i,j) * g%dy_Cu(i,j)
1010  g%IareaCu(i,j) = 0.0
1011  if (g%areaCu(i,j) > 0.0) g%IareaCu(i,j) = g%mask2dCu(i,j) / (g%areaCu(i,j))
1012  enddo ; enddo
1013 
1014  do j=jsdb,jedb ; do i=isd,ied
1015  lat = g%geoLatCv(i,j) ; lon = g%geoLonCv(i,j)
1016  if (check_360) then ; lon_p = lon+len_lon ; lon_m = lon-len_lon
1017  else ; lon_p = lon ; lon_m = lon ; endif
1018 
1019  do npt=1,v_pt
1020  if (((lat >= v_lat(1,npt)) .and. (lat <= v_lat(2,npt))) .and. &
1021  (((lon >= v_lon(1,npt)) .and. (lon <= v_lon(2,npt))) .or. &
1022  ((lon_p >= v_lon(1,npt)) .and. (lon_p <= v_lon(2,npt))) .or. &
1023  ((lon_m >= v_lon(1,npt)) .and. (lon_m <= v_lon(2,npt)))) ) then
1024  g%dx_Cv(i,j) = g%mask2dCv(i,j) * m_to_l*min(l_to_m*g%dxCv(i,j), max(v_width(npt), 0.0))
1025  if (i>=g%isc .and. i<=g%iec .and. j>=g%jsc .and. j<=g%jec) then ! Limit messages/checking to compute domain
1026  if ( g%mask2dCv(i,j) == 0.0 ) then
1027  write(*,'(A,2F8.2,A,4F8.2,A)') "read_face_lengths_list : G%mask2dCv=0 at ",lat,lon," (",&
1028  v_lat(1,npt), v_lat(2,npt), v_lon(1,npt), v_lon(2,npt),") so grid metric is unmodified."
1029  else
1030  write(*,'(A,2F8.2,A,4F8.2,A5,F9.2,A1)') &
1031  "read_face_lengths_list : Modifying dx_Cv gridpoint at ",lat,lon," (",&
1032  v_lat(1,npt), v_lat(2,npt), v_lon(1,npt), v_lon(2,npt),") to ",l_to_m*g%dx_Cv(i,j),"m"
1033  endif
1034  endif
1035  endif
1036  enddo
1037 
1038  g%areaCv(i,j) = g%dyCv(i,j) * g%dx_Cv(i,j)
1039  g%IareaCv(i,j) = 0.0
1040  if (g%areaCv(i,j) > 0.0) g%IareaCv(i,j) = g%mask2dCv(i,j) / (g%areaCv(i,j))
1041  enddo ; enddo
1042 
1043  if (num_lines > 0) then
1044  deallocate(u_lat) ; deallocate(u_lon) ; deallocate(u_width)
1045  deallocate(v_lat) ; deallocate(v_lon) ; deallocate(v_width)
1046  endif
1047 
1048  call calltree_leave(trim(mdl)//'()')

◆ reset_face_lengths_named()

subroutine, public mom_shared_initialization::reset_face_lengths_named ( type(dyn_horgrid_type), intent(inout)  G,
type(param_file_type), intent(in)  param_file,
character(len=*), intent(in)  name,
type(unit_scale_type), intent(in), optional  US 
)

This subroutine sets the open face lengths at selected points to restrict passages to their observed widths based on a named set of sizes.

Parameters
[in,out]gThe dynamic horizontal grid
[in]param_fileA structure to parse for run-time parameters
[in]nameThe name for the set of face lengths. Only "global_1deg" is currently implemented.
[in]usA dimensional unit scaling type

Definition at line 622 of file MOM_shared_initialization.F90.

622  type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid
623  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
624  character(len=*), intent(in) :: name !< The name for the set of face lengths. Only "global_1deg"
625  !! is currently implemented.
626  type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type
627 
628  ! Local variables
629  character(len=256) :: mesg ! Message for error messages.
630  real :: m_to_L ! A unit conversion factor [L m-1 ~> nondim]
631  real :: L_to_m ! A unit conversion factor [m L-1 ~> nondim]
632  real :: dx_2 = -1.0, dy_2 = -1.0
633  real :: pi_180
634  integer :: option = -1
635  integer :: i, j, isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
636  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
637  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
638  pi_180 = (4.0*atan(1.0))/180.0
639 
640  select case ( trim(name) )
641  case ("global_1deg") ; option = 1 ; dx_2 = 0.5*1.0
642  case default ; call mom_error(fatal, "reset_face_lengths_named: "//&
643  "Unrecognized channel configuration name "//trim(name))
644  end select
645 
646  m_to_l = 1.0 ; if (present(us)) m_to_l = us%m_to_L
647  l_to_m = 1.0 ; if (present(us)) l_to_m = us%L_to_m
648 
649  if (option==1) then ! 1-degree settings.
650  do j=jsd,jed ; do i=isdb,iedb ! Change any u-face lengths within this loop.
651  dy_2 = dx_2 * g%dyCu(i,j)*g%IdxCu(i,j) * cos(pi_180 * g%geoLatCu(i,j))
652 
653  if ((abs(g%geoLatCu(i,j)-35.5) < dy_2) .and. (g%geoLonCu(i,j) < -4.5) .and. &
654  (g%geoLonCu(i,j) > -6.5)) &
655  g%dy_Cu(i,j) = g%mask2dCu(i,j)*12000.0*m_to_l ! Gibraltar
656 
657  if ((abs(g%geoLatCu(i,j)-12.5) < dy_2) .and. (abs(g%geoLonCu(i,j)-43.0) < dx_2)) &
658  g%dy_Cu(i,j) = g%mask2dCu(i,j)*10000.0*m_to_l ! Red Sea
659 
660  if ((abs(g%geoLatCu(i,j)-40.5) < dy_2) .and. (abs(g%geoLonCu(i,j)-26.0) < dx_2)) &
661  g%dy_Cu(i,j) = g%mask2dCu(i,j)*5000.0*m_to_l ! Dardanelles
662 
663  if ((abs(g%geoLatCu(i,j)-41.5) < dy_2) .and. (abs(g%geoLonCu(i,j)+220.0) < dx_2)) &
664  g%dy_Cu(i,j) = g%mask2dCu(i,j)*35000.0*m_to_l ! Tsugaru strait at 140.0e
665 
666  if ((abs(g%geoLatCu(i,j)-45.5) < dy_2) .and. (abs(g%geoLonCu(i,j)+217.5) < 0.9)) &
667  g%dy_Cu(i,j) = g%mask2dCu(i,j)*15000.0*m_to_l ! Betw Hokkaido and Sakhalin at 217&218 = 142e
668 
669  ! Greater care needs to be taken in the tripolar region.
670  if ((abs(g%geoLatCu(i,j)-80.84) < 0.2) .and. (abs(g%geoLonCu(i,j)+64.9) < 0.8)) &
671  g%dy_Cu(i,j) = g%mask2dCu(i,j)*38000.0*m_to_l ! Smith Sound in Canadian Arch - tripolar region
672 
673  enddo ; enddo
674 
675  do j=jsdb,jedb ; do i=isd,ied ! Change any v-face lengths within this loop.
676  dy_2 = dx_2 * g%dyCv(i,j)*g%IdxCv(i,j) * cos(pi_180 * g%geoLatCv(i,j))
677  if ((abs(g%geoLatCv(i,j)-41.0) < dy_2) .and. (abs(g%geoLonCv(i,j)-28.5) < dx_2)) &
678  g%dx_Cv(i,j) = g%mask2dCv(i,j)*2500.0*m_to_l ! Bosporus - should be 1000.0 m wide.
679 
680  if ((abs(g%geoLatCv(i,j)-13.0) < dy_2) .and. (abs(g%geoLonCv(i,j)-42.5) < dx_2)) &
681  g%dx_Cv(i,j) = g%mask2dCv(i,j)*10000.0*m_to_l ! Red Sea
682 
683  if ((abs(g%geoLatCv(i,j)+2.8) < 0.8) .and. (abs(g%geoLonCv(i,j)+241.5) < dx_2)) &
684  g%dx_Cv(i,j) = g%mask2dCv(i,j)*40000.0*m_to_l ! Makassar Straits at 241.5 W = 118.5 E
685 
686  if ((abs(g%geoLatCv(i,j)-0.56) < 0.5) .and. (abs(g%geoLonCv(i,j)+240.5) < dx_2)) &
687  g%dx_Cv(i,j) = g%mask2dCv(i,j)*80000.0*m_to_l ! entry to Makassar Straits at 240.5 W = 119.5 E
688 
689  if ((abs(g%geoLatCv(i,j)-0.19) < 0.5) .and. (abs(g%geoLonCv(i,j)+230.5) < dx_2)) &
690  g%dx_Cv(i,j) = g%mask2dCv(i,j)*25000.0*m_to_l ! Channel betw N Guinea and Halmahara 230.5 W = 129.5 E
691 
692  if ((abs(g%geoLatCv(i,j)-0.19) < 0.5) .and. (abs(g%geoLonCv(i,j)+229.5) < dx_2)) &
693  g%dx_Cv(i,j) = g%mask2dCv(i,j)*25000.0*m_to_l ! Channel betw N Guinea and Halmahara 229.5 W = 130.5 E
694 
695  if ((abs(g%geoLatCv(i,j)-0.0) < 0.25) .and. (abs(g%geoLonCv(i,j)+228.5) < dx_2)) &
696  g%dx_Cv(i,j) = g%mask2dCv(i,j)*25000.0*m_to_l ! Channel betw N Guinea and Halmahara 228.5 W = 131.5 E
697 
698  if ((abs(g%geoLatCv(i,j)+8.5) < 0.5) .and. (abs(g%geoLonCv(i,j)+244.5) < dx_2)) &
699  g%dx_Cv(i,j) = g%mask2dCv(i,j)*20000.0*m_to_l ! Lombok Straits at 244.5 W = 115.5 E
700 
701  if ((abs(g%geoLatCv(i,j)+8.5) < 0.5) .and. (abs(g%geoLonCv(i,j)+235.5) < dx_2)) &
702  g%dx_Cv(i,j) = g%mask2dCv(i,j)*20000.0*m_to_l ! Timor Straits at 235.5 W = 124.5 E
703 
704  if ((abs(g%geoLatCv(i,j)-52.5) < dy_2) .and. (abs(g%geoLonCv(i,j)+218.5) < dx_2)) &
705  g%dx_Cv(i,j) = g%mask2dCv(i,j)*2500.0*m_to_l ! Russia and Sakhalin Straits at 218.5 W = 141.5 E
706 
707  ! Greater care needs to be taken in the tripolar region.
708  if ((abs(g%geoLatCv(i,j)-76.8) < 0.06) .and. (abs(g%geoLonCv(i,j)+88.7) < dx_2)) &
709  g%dx_Cv(i,j) = g%mask2dCv(i,j)*8400.0*m_to_l ! Jones Sound in Canadian Arch - tripolar region
710 
711  enddo ; enddo
712  endif
713 
714  ! These checks apply regardless of the chosen option.
715 
716  do j=jsd,jed ; do i=isdb,iedb
717  if (l_to_m*g%dy_Cu(i,j) > l_to_m*g%dyCu(i,j)) then
718  write(mesg,'("dy_Cu of ",ES11.4," exceeds unrestricted width of ",ES11.4,&
719  &" by ",ES11.4," at lon/lat of ", ES11.4, ES11.4)') &
720  l_to_m*g%dy_Cu(i,j), l_to_m*g%dyCu(i,j), l_to_m*g%dy_Cu(i,j)-l_to_m*g%dyCu(i,j), &
721  g%geoLonCu(i,j), g%geoLatCu(i,j)
722  call mom_error(fatal,"reset_face_lengths_named "//mesg)
723  endif
724  g%areaCu(i,j) = g%dxCu(i,j) * g%dy_Cu(i,j)
725  g%IareaCu(i,j) = 0.0
726  if (g%areaCu(i,j) > 0.0) g%IareaCu(i,j) = g%mask2dCu(i,j) / (g%areaCu(i,j))
727  enddo ; enddo
728 
729  do j=jsdb,jedb ; do i=isd,ied
730  if (l_to_m*g%dx_Cv(i,j) > l_to_m*g%dxCv(i,j)) then
731  write(mesg,'("dx_Cv of ",ES11.4," exceeds unrestricted width of ",ES11.4,&
732  &" by ",ES11.4, " at lon/lat of ", ES11.4, ES11.4)') &
733  l_to_m*g%dx_Cv(i,j), l_to_m*g%dxCv(i,j), l_to_m*g%dx_Cv(i,j)-l_to_m*g%dxCv(i,j), &
734  g%geoLonCv(i,j), g%geoLatCv(i,j)
735 
736  call mom_error(fatal,"reset_face_lengths_named "//mesg)
737  endif
738  g%areaCv(i,j) = g%dyCv(i,j) * g%dx_Cv(i,j)
739  g%IareaCv(i,j) = 0.0
740  if (g%areaCv(i,j) > 0.0) g%IareaCv(i,j) = g%mask2dCv(i,j) / (g%areaCv(i,j))
741  enddo ; enddo
742 

◆ set_rotation_beta_plane()

subroutine, public mom_shared_initialization::set_rotation_beta_plane ( real, dimension(g%isdb:g%iedb,g%jsdb:g%jedb), intent(out)  f,
type(dyn_horgrid_type), intent(in)  G,
type(param_file_type), intent(in)  param_file,
type(unit_scale_type), intent(in), optional  US 
)

This subroutine sets up the Coriolis parameter for a beta-plane or f-plane.

Parameters
[in]gThe dynamic horizontal grid
[out]fCoriolis parameter (vertical component) [T-1 ~> s-1]
[in]param_fileA structure to parse for run-time parameters
[in]usA dimensional unit scaling type

Definition at line 497 of file MOM_shared_initialization.F90.

497  type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid
498  real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), &
499  intent(out) :: f !< Coriolis parameter (vertical component) [T-1 ~> s-1]
500  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
501  type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type
502 
503 ! This subroutine sets up the Coriolis parameter for a beta-plane
504  integer :: I, J
505  real :: f_0 ! The reference value of the Coriolis parameter [T-1 ~> s-1]
506  real :: beta ! The meridional gradient of the Coriolis parameter [T-1 m-1 ~> s-1 m-1]
507  real :: y_scl, Rad_Earth
508  real :: T_to_s ! A time unit conversion factor
509  real :: PI
510  character(len=40) :: mdl = "set_rotation_beta_plane" ! This subroutine's name.
511  character(len=200) :: axis_units
512 
513  call calltree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")
514 
515  t_to_s = 1.0 ; if (present(us)) t_to_s = us%T_to_s
516 
517  call get_param(param_file, mdl, "F_0", f_0, &
518  "The reference value of the Coriolis parameter with the "//&
519  "betaplane option.", units="s-1", default=0.0, scale=t_to_s)
520  call get_param(param_file, mdl, "BETA", beta, &
521  "The northward gradient of the Coriolis parameter with "//&
522  "the betaplane option.", units="m-1 s-1", default=0.0, scale=t_to_s)
523  call get_param(param_file, mdl, "AXIS_UNITS", axis_units, default="degrees")
524 
525  pi = 4.0*atan(1.0)
526  select case (axis_units(1:1))
527  case ("d")
528  call get_param(param_file, mdl, "RAD_EARTH", rad_earth, &
529  "The radius of the Earth.", units="m", default=6.378e6)
530  y_scl = rad_earth/pi
531  case ("k"); y_scl = 1.e3
532  case ("m"); y_scl = 1.
533  case ("c"); y_scl = 1.e-2
534  case default ; call mom_error(fatal, &
535  " set_rotation_beta_plane: unknown AXIS_UNITS = "//trim(axis_units))
536  end select
537 
538  do i=g%IsdB,g%IedB ; do j=g%JsdB,g%JedB
539  f(i,j) = f_0 + beta * ( g%geoLatBu(i,j) * y_scl )
540  enddo ; enddo
541 
542  call calltree_leave(trim(mdl)//'()')

◆ set_rotation_planetary()

subroutine, public mom_shared_initialization::set_rotation_planetary ( real, dimension(g%isdb:g%iedb,g%jsdb:g%jedb), intent(out)  f,
type(dyn_horgrid_type), intent(in)  G,
type(param_file_type), intent(in)  param_file,
type(unit_scale_type), intent(in), optional  US 
)

This subroutine sets up the Coriolis parameter for a sphere.

Parameters
[in]gThe dynamic horizontal grid
[out]fCoriolis parameter (vertical component) [T-1 ~> s-1]
[in]param_fileA structure to parse for run-time parameters
[in]usA dimensional unit scaling type

Definition at line 464 of file MOM_shared_initialization.F90.

464  type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid
465  real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB), &
466  intent(out) :: f !< Coriolis parameter (vertical component) [T-1 ~> s-1]
467  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
468  type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type
469 
470 ! This subroutine sets up the Coriolis parameter for a sphere
471  character(len=30) :: mdl = "set_rotation_planetary" ! This subroutine's name.
472  integer :: I, J
473  real :: PI
474  real :: omega ! The planetary rotation rate [T-1 ~> s-1]
475  real :: T_to_s ! A time unit conversion factor
476 
477  call calltree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")
478 
479  t_to_s = 1.0 ; if (present(us)) t_to_s = us%T_to_s
480 
481  call get_param(param_file, "set_rotation_planetary", "OMEGA", omega, &
482  "The rotation rate of the earth.", units="s-1", &
483  default=7.2921e-5, scale=t_to_s)
484  pi = 4.0*atan(1.0)
485 
486  do i=g%IsdB,g%IedB ; do j=g%JsdB,g%JedB
487  f(i,j) = ( 2.0 * omega ) * sin( ( pi * g%geoLatBu(i,j) ) / 180.)
488  enddo ; enddo
489 
490  call calltree_leave(trim(mdl)//'()')

◆ set_velocity_depth_max()

subroutine, public mom_shared_initialization::set_velocity_depth_max ( type(dyn_horgrid_type), intent(inout)  G)

Set the bathymetry at velocity points to be the maximum of the depths at the neighoring tracer points.

Parameters
[in,out]gThe dynamic horizontal grid

Definition at line 1109 of file MOM_shared_initialization.F90.

1109  type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid
1110  ! This subroutine sets the 4 bottom depths at velocity points to be the
1111  ! maximum of the adjacent depths.
1112  integer :: i, j
1113 
1114  do i=g%isd,g%ied-1 ; do j=g%jsd,g%jed
1115  g%Dblock_u(i,j) = g%mask2dCu(i,j) * max(g%bathyT(i,j), g%bathyT(i+1,j))
1116  g%Dopen_u(i,j) = g%Dblock_u(i,j)
1117  enddo ; enddo
1118  do i=g%isd,g%ied ; do j=g%jsd,g%jed-1
1119  g%Dblock_v(i,j) = g%mask2dCv(i,j) * max(g%bathyT(i,j), g%bathyT(i,j+1))
1120  g%Dopen_v(i,j) = g%Dblock_v(i,j)
1121  enddo ; enddo

◆ set_velocity_depth_min()

subroutine, public mom_shared_initialization::set_velocity_depth_min ( type(dyn_horgrid_type), intent(inout)  G)

Set the bathymetry at velocity points to be the minimum of the depths at the neighoring tracer points.

Parameters
[in,out]gThe dynamic horizontal grid

Definition at line 1129 of file MOM_shared_initialization.F90.

1129  type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid
1130  ! This subroutine sets the 4 bottom depths at velocity points to be the
1131  ! minimum of the adjacent depths.
1132  integer :: i, j
1133 
1134  do i=g%isd,g%ied-1 ; do j=g%jsd,g%jed
1135  g%Dblock_u(i,j) = g%mask2dCu(i,j) * min(g%bathyT(i,j), g%bathyT(i+1,j))
1136  g%Dopen_u(i,j) = g%Dblock_u(i,j)
1137  enddo ; enddo
1138  do i=g%isd,g%ied ; do j=g%jsd,g%jed-1
1139  g%Dblock_v(i,j) = g%mask2dCv(i,j) * min(g%bathyT(i,j), g%bathyT(i,j+1))
1140  g%Dopen_v(i,j) = g%Dblock_v(i,j)
1141  enddo ; enddo

◆ write_ocean_geometry_file()

subroutine, public mom_shared_initialization::write_ocean_geometry_file ( type(dyn_horgrid_type), intent(inout)  G,
type(param_file_type), intent(in)  param_file,
character(len=*), intent(in)  directory,
character(len=*), intent(in), optional  geom_file,
type(unit_scale_type), intent(in), optional  US 
)

Write out a file describing the topography, Coriolis parameter, grid locations and various other fixed fields from the grid.

Parameters
[in,out]gThe dynamic horizontal grid
[in]param_fileParameter file structure
[in]directoryThe directory into which to place the geometry file.
[in]geom_fileIf present, the name of the geometry file (otherwise the file is "ocean_geometry")
[in]usA dimensional unit scaling type

Definition at line 1178 of file MOM_shared_initialization.F90.

1178  type(dyn_horgrid_type), intent(inout) :: G !< The dynamic horizontal grid
1179  type(param_file_type), intent(in) :: param_file !< Parameter file structure
1180  character(len=*), intent(in) :: directory !< The directory into which to place the geometry file.
1181  character(len=*), optional, intent(in) :: geom_file !< If present, the name of the geometry file
1182  !! (otherwise the file is "ocean_geometry")
1183  type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type
1184 
1185  ! Local variables.
1186  character(len=240) :: filepath
1187  character(len=40) :: mdl = "write_ocean_geometry_file"
1188  integer, parameter :: nFlds=23
1189  type(vardesc) :: vars(nFlds)
1190  type(fieldtype) :: fields(nFlds)
1191  real :: Z_to_m_scale ! A unit conversion factor from Z to m.
1192  real :: s_to_T_scale ! A unit conversion factor from T-1 to s-1.
1193  real :: L_to_m_scale ! A unit conversion factor from L to m.
1194  integer :: unit
1195  integer :: file_threading
1196  integer :: nFlds_used
1197  integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq
1198  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
1199  logical :: multiple_files
1200  real, dimension(G%isd :G%ied ,G%jsd :G%jed ) :: out_h
1201  real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB) :: out_q
1202  real, dimension(G%IsdB:G%IedB,G%jsd :G%jed ) :: out_u
1203  real, dimension(G%isd :G%ied ,G%JsdB:G%JedB) :: out_v
1204 
1205  call calltree_enter('write_ocean_geometry_file()')
1206 
1207  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1208  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
1209  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1210  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
1211 
1212  z_to_m_scale = 1.0 ; if (present(us)) z_to_m_scale = us%Z_to_m
1213  s_to_t_scale = 1.0 ; if (present(us)) s_to_t_scale = us%s_to_T
1214  l_to_m_scale = 1.0 ; if (present(us)) l_to_m_scale = us%L_to_m
1215 
1216 ! vardesc is a structure defined in MOM_io.F90. The elements of
1217 ! this structure, in order, are:
1218 ! (1) the variable name for the NetCDF file
1219 ! (2) the variable's long name
1220 ! (3) a character indicating the horizontal grid, which may be '1' (column),
1221 ! 'h', 'q', 'u', or 'v', for the corresponding C-grid variable
1222 ! (4) a character indicating the vertical grid, which may be 'L' (layer),
1223 ! 'i' (interface), or '1' (no vertical location)
1224 ! (5) a character indicating the time levels of the field, which may be
1225 ! 's' (snap-shot), 'p' (periodic), or '1' (no time variation)
1226 ! (6) the variable's units
1227  vars(1) = var_desc("geolatb","degree","latitude at corner (Bu) points",'q','1','1')
1228  vars(2) = var_desc("geolonb","degree","longitude at corner (Bu) points",'q','1','1')
1229  vars(3) = var_desc("geolat","degree", "latitude at tracer (T) points", 'h','1','1')
1230  vars(4) = var_desc("geolon","degree","longitude at tracer (T) points",'h','1','1')
1231  vars(5) = var_desc("D","meter","Basin Depth",'h','1','1')
1232  vars(6) = var_desc("f","s-1","Coriolis Parameter",'q','1','1')
1233  vars(7) = var_desc("dxCv","m","Zonal grid spacing at v points",'v','1','1')
1234  vars(8) = var_desc("dyCu","m","Meridional grid spacing at u points",'u','1','1')
1235  vars(9) = var_desc("dxCu","m","Zonal grid spacing at u points",'u','1','1')
1236  vars(10)= var_desc("dyCv","m","Meridional grid spacing at v points",'v','1','1')
1237  vars(11)= var_desc("dxT","m","Zonal grid spacing at h points",'h','1','1')
1238  vars(12)= var_desc("dyT","m","Meridional grid spacing at h points",'h','1','1')
1239  vars(13)= var_desc("dxBu","m","Zonal grid spacing at q points",'q','1','1')
1240  vars(14)= var_desc("dyBu","m","Meridional grid spacing at q points",'q','1','1')
1241  vars(15)= var_desc("Ah","m2","Area of h cells",'h','1','1')
1242  vars(16)= var_desc("Aq","m2","Area of q cells",'q','1','1')
1243 
1244  vars(17)= var_desc("dxCvo","m","Open zonal grid spacing at v points",'v','1','1')
1245  vars(18)= var_desc("dyCuo","m","Open meridional grid spacing at u points",'u','1','1')
1246  vars(19)= var_desc("wet", "nondim", "land or ocean?", 'h','1','1')
1247 
1248  vars(20) = var_desc("Dblock_u","m","Blocked depth at u points",'u','1','1')
1249  vars(21) = var_desc("Dopen_u","m","Open depth at u points",'u','1','1')
1250  vars(22) = var_desc("Dblock_v","m","Blocked depth at v points",'v','1','1')
1251  vars(23) = var_desc("Dopen_v","m","Open depth at v points",'v','1','1')
1252 
1253 
1254  nflds_used = 19 ; if (g%bathymetry_at_vel) nflds_used = 23
1255 
1256  if (present(geom_file)) then
1257  filepath = trim(directory) // trim(geom_file)
1258  else
1259  filepath = trim(directory) // "ocean_geometry"
1260  endif
1261 
1262  out_h(:,:) = 0.0
1263  out_u(:,:) = 0.0
1264  out_v(:,:) = 0.0
1265  out_q(:,:) = 0.0
1266 
1267  call get_param(param_file, mdl, "PARALLEL_RESTARTFILES", multiple_files, &
1268  "If true, each processor writes its own restart file, "//&
1269  "otherwise a single restart file is generated", &
1270  default=.false.)
1271  file_threading = single_file
1272  if (multiple_files) file_threading = multiple
1273 
1274  call create_file(unit, trim(filepath), vars, nflds_used, fields, &
1275  file_threading, dg=g)
1276 
1277  do j=jsq,jeq; do i=isq,ieq; out_q(i,j) = g%geoLatBu(i,j); enddo ; enddo
1278  call write_field(unit, fields(1), g%Domain%mpp_domain, out_q)
1279  do j=jsq,jeq; do i=isq,ieq; out_q(i,j) = g%geoLonBu(i,j); enddo ; enddo
1280  call write_field(unit, fields(2), g%Domain%mpp_domain, out_q)
1281  call write_field(unit, fields(3), g%Domain%mpp_domain, g%geoLatT)
1282  call write_field(unit, fields(4), g%Domain%mpp_domain, g%geoLonT)
1283 
1284  do j=js,je ; do i=is,ie ; out_h(i,j) = z_to_m_scale*g%bathyT(i,j) ; enddo ; enddo
1285  call write_field(unit, fields(5), g%Domain%mpp_domain, out_h)
1286  do j=jsq,jeq ; do i=isq,ieq ; out_q(i,j) = s_to_t_scale*g%CoriolisBu(i,j) ; enddo ; enddo
1287  call write_field(unit, fields(6), g%Domain%mpp_domain, out_q)
1288 
1289  ! I think that all of these copies are holdovers from a much earlier
1290  ! ancestor code in which many of the metrics were macros that could have
1291  ! had reduced dimensions, and that they are no longer needed in MOM6. -RWH
1292  do j=jsq,jeq ; do i=is,ie ; out_v(i,j) = l_to_m_scale*g%dxCv(i,j) ; enddo ; enddo
1293  call write_field(unit, fields(7), g%Domain%mpp_domain, out_v)
1294  do j=js,je ; do i=isq,ieq ; out_u(i,j) = l_to_m_scale*g%dyCu(i,j) ; enddo ; enddo
1295  call write_field(unit, fields(8), g%Domain%mpp_domain, out_u)
1296 
1297  do j=js,je ; do i=isq,ieq ; out_u(i,j) = l_to_m_scale*g%dxCu(i,j) ; enddo ; enddo
1298  call write_field(unit, fields(9), g%Domain%mpp_domain, out_u)
1299  do j=jsq,jeq ; do i=is,ie ; out_v(i,j) = l_to_m_scale*g%dyCv(i,j) ; enddo ; enddo
1300  call write_field(unit, fields(10), g%Domain%mpp_domain, out_v)
1301 
1302  do j=js,je ; do i=is,ie ; out_h(i,j) = l_to_m_scale*g%dxT(i,j); enddo ; enddo
1303  call write_field(unit, fields(11), g%Domain%mpp_domain, out_h)
1304  do j=js,je ; do i=is,ie ; out_h(i,j) = l_to_m_scale*g%dyT(i,j) ; enddo ; enddo
1305  call write_field(unit, fields(12), g%Domain%mpp_domain, out_h)
1306 
1307  do j=jsq,jeq ; do i=isq,ieq ; out_q(i,j) = l_to_m_scale*g%dxBu(i,j) ; enddo ; enddo
1308  call write_field(unit, fields(13), g%Domain%mpp_domain, out_q)
1309  do j=jsq,jeq ; do i=isq,ieq ; out_q(i,j) = l_to_m_scale*g%dyBu(i,j) ; enddo ; enddo
1310  call write_field(unit, fields(14), g%Domain%mpp_domain, out_q)
1311 
1312  do j=js,je ; do i=is,ie ; out_h(i,j) = g%areaT(i,j) ; enddo ; enddo
1313  call write_field(unit, fields(15), g%Domain%mpp_domain, out_h)
1314  do j=jsq,jeq ; do i=isq,ieq ; out_q(i,j) = g%areaBu(i,j) ; enddo ; enddo
1315  call write_field(unit, fields(16), g%Domain%mpp_domain, out_q)
1316 
1317  do j=jsq,jeq ; do i=is,ie ; out_v(i,j) = l_to_m_scale*g%dx_Cv(i,j) ; enddo ; enddo
1318  call write_field(unit, fields(17), g%Domain%mpp_domain, out_v)
1319  do j=js,je ; do i=isq,ieq ; out_u(i,j) = l_to_m_scale*g%dy_Cu(i,j) ; enddo ; enddo
1320  call write_field(unit, fields(18), g%Domain%mpp_domain, out_u)
1321  call write_field(unit, fields(19), g%Domain%mpp_domain, g%mask2dT)
1322 
1323  if (g%bathymetry_at_vel) then
1324  do j=js,je ; do i=isq,ieq ; out_u(i,j) = z_to_m_scale*g%Dblock_u(i,j) ; enddo ; enddo
1325  call write_field(unit, fields(20), g%Domain%mpp_domain, out_u)
1326  do j=js,je ; do i=isq,ieq ; out_u(i,j) = z_to_m_scale*g%Dopen_u(i,j) ; enddo ; enddo
1327  call write_field(unit, fields(21), g%Domain%mpp_domain, out_u)
1328  do j=jsq,jeq ; do i=is,ie ; out_v(i,j) = z_to_m_scale*g%Dblock_v(i,j) ; enddo ; enddo
1329  call write_field(unit, fields(22), g%Domain%mpp_domain, out_v)
1330  do j=jsq,jeq ; do i=is,ie ; out_v(i,j) = z_to_m_scale*g%Dopen_v(i,j) ; enddo ; enddo
1331  call write_field(unit, fields(23), g%Domain%mpp_domain, out_v)
1332  endif
1333 
1334  call close_file(unit)
1335 
1336  call calltree_leave('write_ocean_geometry_file()')