MOM_fixed_initialization.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Initializes fixed aspects of the model, such as horizontal grid metrics,
6!! topography and Coriolis.
8
9use mom_debugging, only : hchksum, qchksum, uvchksum
10use mom_domains, only : pass_var
12use mom_error_handler, only : mom_mesg, mom_error, fatal, warning, is_root_pe
14use mom_file_parser, only : get_param, read_param, log_param, param_file_type
15use mom_file_parser, only : log_version
16use mom_io, only : slasher
32
47
48implicit none ; private
49
50public mom_initialize_fixed, mom_initialize_rotation, mom_initialize_topography
51
52contains
53
54! -----------------------------------------------------------------------------
55!> MOM_initialize_fixed sets up time-invariant quantities related to MOM6's
56!! horizontal grid, bathymetry, and the Coriolis parameter.
57subroutine mom_initialize_fixed(G, US, OBC, PF)
58 type(dyn_horgrid_type), intent(inout) :: g !< The ocean's grid structure.
59 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
60 type(ocean_obc_type), pointer :: obc !< Open boundary structure.
61 type(param_file_type), intent(in) :: pf !< A structure indicating the open file
62 !! to parse for model parameter values.
63
64 ! Local variables
65 character(len=200) :: config
66 logical :: obc_projection_bug, open_corners, enable_bugs
67 logical :: read_porous_file, read_meansl_file
68 character(len=40) :: mdl = "MOM_fixed_initialization" ! This module's name.
69 integer :: i, j
70 logical :: debug
71 ! This include declares and sets the variable "version".
72# include "version_variable.h"
73
74 call calltree_enter("MOM_initialize_fixed(), MOM_fixed_initialization.F90")
75 call get_param(pf, mdl, "DEBUG", debug, default=.false.)
76
77 ! Set up the parameters of the physical domain (i.e. the grid), G
78 call set_grid_metrics(g, pf, us)
79
80 ! Read time mean sea level from file
81 call get_param(pf, mdl, "READ_MEAN_SEA_LEVEL", read_meansl_file, &
82 "If true, use a 2D map for time mean sea level, which is used to calculate "// &
83 "time mean ocean total thickness.", default=.false.)
84 if (read_meansl_file) &
85 call set_meansl_from_file(g%meanSL, g, pf, us)
86
87 ! Set up the bottom depth, G%bathyT either analytically or from file
88 ! This also sets G%max_depth based on the input parameter MAXIMUM_DEPTH,
89 ! or, if absent, is diagnosed as G%max_depth = max( G%D(:,:) )
90 call mom_initialize_topography(g%bathyT, g%max_depth, g, pf, us, meansl=g%meanSL)
91
92 ! To initialize masks, the bathymetry in halo regions must be filled in
93 call pass_var(g%bathyT, g%Domain)
94
95 ! Determine the position of any open boundaries and create OBC
96 call open_boundary_config(g, us, pf, obc)
97
98 ! Make bathymetry (if OBC_PROJECTION_BUG) and masks consistent with open boundaries.
99 if (associated(obc)) then
100 call get_param(pf, mdl, "ENABLE_BUGS_BY_DEFAULT", enable_bugs, &
101 default=.true., do_not_log=.true.) ! This is logged from MOM.F90.
102 call get_param(pf, mdl, "OBC_PROJECTION_BUG", obc_projection_bug, &
103 "If false, use only interior ocean points at OBCs to specify several "//&
104 "calculations at OBC points, and it avoids applying a land mask at the "//&
105 "bay-like intersection of orthogonal OBC segments. Otherwise the "//&
106 "calculation of terms like the potential vorticity used in the barotropic "//&
107 "solver relies on bathymetry or other fields being projected outward across "//&
108 "OBCs. This option changes answers for some configurations that use OBCs.", &
109 default=enable_bugs)
110 open_corners = .not.obc_projection_bug
111
112 if (obc_projection_bug .and. read_meansl_file) &
113 ! OBC_projection_bug modifies bathyT outside of the open boundaries, so meanSL would have to be
114 ! modified as well.
115 call mom_error(fatal, "MOM_initialize_fixed: To read mean sea level file, "//&
116 "OBC_PROJECTION_BUG needs to be False.")
117
118 ! This call sets masks that prohibit flow over any point interpreted as land
119 if (obc_projection_bug) &
120 call open_boundary_impose_normal_slope(obc, g, g%bathyT)
121 call initialize_masks(g, pf, us, obc_dir_u=obc%segnum_u, obc_dir_v=obc%segnum_v, &
122 open_corner_obcs=open_corners)
123 ! Make OBC mask consistent with land mask
124 call open_boundary_impose_land_mask(obc, g, g%areaCu, g%areaCv, us)
125 else
126 call initialize_masks(g, pf, us)
127 endif
128
129 if (debug) then
130 call hchksum(g%bathyT, 'MOM_initialize_fixed: depth ', g%HI, haloshift=1, unscale=us%Z_to_m)
131 call hchksum(g%mask2dT, 'MOM_initialize_fixed: mask2dT ', g%HI)
132 call uvchksum('MOM_initialize_fixed: mask2dC[uv]', g%mask2dCu, &
133 g%mask2dCv, g%HI)
134 call qchksum(g%mask2dBu, 'MOM_initialize_fixed: mask2dBu ', g%HI)
135 endif
136
137 ! Set up other fixed quantities
138 ! Parameters below are logged under "module MOM_fixed_initialization".
139 call log_version(pf, mdl, version, "")
140 ! Modulate geometric scales according to geography.
141 call get_param(pf, mdl, "CHANNEL_CONFIG", config, &
142 "A parameter that determines which set of channels are \n"//&
143 "restricted to specific widths. Options are:\n"//&
144 " \t none - All channels have the grid width.\n"//&
145 " \t global_1deg - Sets 16 specific channels appropriate \n"//&
146 " \t\t for a 1-degree model, as used in CM2G.\n"//&
147 " \t list - Read the channel locations and widths from a \n"//&
148 " \t\t text file, like MOM_channel_list in the MOM_SIS \n"//&
149 " \t\t test case.\n"//&
150 " \t file - Read open face widths everywhere from a \n"//&
151 " \t\t NetCDF file on the model grid.", &
152 default="none")
153 select case ( trim(config) )
154 case ("none")
155 case ("list") ; call reset_face_lengths_list(g, pf, us)
156 case ("file") ; call reset_face_lengths_file(g, pf, us)
157 case ("global_1deg") ; call reset_face_lengths_named(g, pf, trim(config), us)
158 case default ; call mom_error(fatal, "MOM_initialize_fixed: "// &
159 "Unrecognized channel configuration "//trim(config))
160 end select
161
162 ! This call sets the topography at velocity points.
163 if (g%bathymetry_at_vel) then
164 call get_param(pf, mdl, "VELOCITY_DEPTH_CONFIG", config, &
165 "A string that determines how the topography is set at "//&
166 "velocity points. This may be 'min' or 'max'.", &
167 default="max")
168 select case ( trim(config) )
169 case ("max") ; call set_velocity_depth_max(g)
170 case ("min") ; call set_velocity_depth_min(g)
171 case default ; call mom_error(fatal, "MOM_initialize_fixed: "// &
172 "Unrecognized velocity depth configuration "//trim(config))
173 end select
174 endif
175
176 ! Read sub-grid scale topography parameters at velocity points used for porous barrier calculation
177 ! TODO: The following routine call may eventually be merged as one of the CHANNEL_CONFIG options
178 call get_param(pf, mdl, "SUBGRID_TOPO_AT_VEL", read_porous_file, &
179 "If true, use variables from TOPO_AT_VEL_FILE as parameters for porous barrier.", &
180 default=.false.)
181 if (read_porous_file) &
182 call set_subgrid_topo_at_vel_from_file(g, pf, us)
183
184 ! Calculate the value of the Coriolis parameter at the latitude !
185 ! of the q grid points [T-1 ~> s-1].
186 call mom_initialize_rotation(g%CoriolisBu, g, pf, us=us)
187 ! Calculate the components of grad f (beta)
188 call mom_calculate_grad_coriolis(g%dF_dx, g%dF_dy, g, us=us)
189 ! Calculate the square of the Coriolis parameter
190 do i=g%IsdB,g%IedB ; do j=g%JsdB,g%JedB
191 g%Coriolis2Bu(i,j) = g%CoriolisBu(i,j)**2
192 enddo ; enddo
193
194 if (debug) then
195 call qchksum(g%CoriolisBu, "MOM_initialize_fixed: f ", g%HI, unscale=us%s_to_T)
196 call qchksum(g%Coriolis2Bu, "MOM_initialize_fixed: f2 ", g%HI, unscale=us%s_to_T**2)
197 call hchksum(g%dF_dx, "MOM_initialize_fixed: dF_dx ", g%HI, unscale=us%m_to_L*us%s_to_T)
198 call hchksum(g%dF_dy, "MOM_initialize_fixed: dF_dy ", g%HI, unscale=us%m_to_L*us%s_to_T)
199 endif
200
201 call initialize_grid_rotation_angle(g, pf)
202
203 ! Compute global integrals of grid values for later use in scalar diagnostics !
204 call compute_global_grid_integrals(g, us=us)
205
206 call calltree_leave('MOM_initialize_fixed()')
207
208end subroutine mom_initialize_fixed
209
210!> MOM_initialize_topography makes the appropriate call to set up the bathymetry in units of [Z ~> m].
211subroutine mom_initialize_topography(D, max_depth, G, PF, US, meanSL)
212 type(dyn_horgrid_type), intent(in) :: g !< The dynamic horizontal grid type
213 real, dimension(G%isd:G%ied,G%jsd:G%jed), &
214 intent(out) :: d !< Ocean bottom depth [Z ~> m]
215 type(param_file_type), intent(in) :: pf !< Parameter file structure
216 real, intent(out) :: max_depth !< Maximum depth or geometric thickness,
217 !! with meanSL present, of model [Z ~> m]
218 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
219 real, dimension(G%isd:G%ied,G%jsd:G%jed), &
220 optional, intent(in) :: meansl !< Mean sea level [Z ~> m]
221
222 ! This subroutine makes the appropriate call to set up the bottom depth.
223 ! This is a separate subroutine so that it can be made public and shared with
224 ! the ice-sheet code or other components.
225
226 ! Local variables
227 real :: max_depth_default = -1.e9 ! Default value of MAXIMUM_DEPTH parameter [m]
228 character(len=40) :: mdl = "MOM_initialize_topography" ! This subroutine's name.
229 character(len=200) :: config
230 real, dimension(G%isd:G%ied, G%jsd:G%jed) :: d_meansl ! depth (positive below meanSL) referenced
231 ! to meanSL. A temporary field used to diagnose maximum
232 ! static column thickness. D_meanSL = D + meanSL [Z ~> m].
233 integer :: i, j
234
235 call get_param(pf, mdl, "TOPO_CONFIG", config, &
236 "This specifies how bathymetry is specified: \n"//&
237 " \t file - read bathymetric information from the file \n"//&
238 " \t\t specified by (TOPO_FILE).\n"//&
239 " \t flat - flat bottom set to MAXIMUM_DEPTH. \n"//&
240 " \t bowl - an analytically specified bowl-shaped basin \n"//&
241 " \t\t ranging between MAXIMUM_DEPTH and MINIMUM_DEPTH. \n"//&
242 " \t spoon - a similar shape to 'bowl', but with an vertical \n"//&
243 " \t\t wall at the southern face. \n"//&
244 " \t halfpipe - a zonally uniform channel with a half-sine \n"//&
245 " \t\t profile in the meridional direction. \n"//&
246 " \t bbuilder - build topography from list of functions. \n"//&
247 " \t benchmark - use the benchmark test case topography. \n"//&
248 " \t Neverworld - use the Neverworld test case topography. \n"//&
249 " \t DOME - use a slope and channel configuration for the \n"//&
250 " \t\t DOME sill-overflow test case. \n"//&
251 " \t ISOMIP - use a slope and channel configuration for the \n"//&
252 " \t\t ISOMIP test case. \n"//&
253 " \t DOME2D - use a shelf and slope configuration for the \n"//&
254 " \t\t DOME2D gravity current/overflow test case. \n"//&
255 " \t Kelvin - flat but with rotated land mask.\n"//&
256 " \t seamount - Gaussian bump for spontaneous motion test case.\n"//&
257 " \t dumbbell - Sloshing channel with reservoirs on both ends.\n"//&
258 " \t shelfwave - exponential slope for shelfwave test case.\n"//&
259 " \t Phillips - ACC-like idealized topography used in the Phillips config.\n"//&
260 " \t dense - Denmark Strait-like dense water formation and overflow.\n"//&
261 " \t USER - call a user modified routine.", &
262 fail_if_missing=.true.)
263 call get_param(pf, mdl, "MAXIMUM_DEPTH", max_depth, units="m", default=max_depth_default, &
264 scale=us%m_to_Z, do_not_log=.true.)
265 select case ( trim(config) )
266 case ("file"); call initialize_topography_from_file(d, g, pf, us)
267 case ("flat"); call initialize_topography_named(d, g, pf, config, max_depth, us)
268 case ("spoon"); call initialize_topography_named(d, g, pf, config, max_depth, us)
269 case ("bowl"); call initialize_topography_named(d, g, pf, config, max_depth, us)
270 case ("halfpipe"); call initialize_topography_named(d, g, pf, config, max_depth, us)
271 case ("DOME"); call dome_initialize_topography(d, g, pf, max_depth, us)
272 case ("ISOMIP"); call isomip_initialize_topography(d, g, pf, max_depth, us)
273 case ("bbuilder"); call basin_builder_topography(d, g, pf, max_depth)
274 case ("benchmark"); call benchmark_initialize_topography(d, g, pf, max_depth, us)
275 case ("Neverworld","Neverland"); call neverworld_initialize_topography(d, g, pf, max_depth)
276 case ("DOME2D"); call dome2d_initialize_topography(d, g, pf, max_depth)
277 case ("Kelvin"); call kelvin_initialize_topography(d, g, pf, max_depth, us)
278 case ("sloshing"); call sloshing_initialize_topography(d, g, pf, max_depth)
279 case ("seamount"); call seamount_initialize_topography(d, g, pf, max_depth)
280 case ("dumbbell"); call dumbbell_initialize_topography(d, g, pf, max_depth)
281 case ("shelfwave"); call shelfwave_initialize_topography(d, g, pf, max_depth, us)
282 case ("Phillips"); call phillips_initialize_topography(d, g, pf, max_depth, us)
283 case ("dense"); call dense_water_initialize_topography(d, g, pf, max_depth)
284 case ("USER"); call user_initialize_topography(d, g, pf, max_depth, us)
285 case default ; call mom_error(fatal,"MOM_initialize_topography: "// &
286 "Unrecognized topography setup '"//trim(config)//"'")
287 end select
288 if (max_depth /= max_depth_default * us%m_to_Z) then
289 call log_param(pf, mdl, "MAXIMUM_DEPTH", max_depth, &
290 "The maximum depth of the ocean.", units="m", unscale=us%Z_to_m)
291 if (trim(config) /= "DOME") then
292 call limit_topography(d, g, pf, max_depth, us)
293 endif
294 else
295 if (present(meansl)) then
296 d_meansl(:,:) = 0.0
297 do j=g%jsc,g%jec ; do i=g%isc,g%iec ; d_meansl(i,j) = d(i,j) + meansl(i,j) ; enddo ; enddo
298 max_depth = diagnosemaximumdepth(d_meansl, g)
299 else
300 max_depth = diagnosemaximumdepth(d, g)
301 endif
302 call log_param(pf, mdl, "!MAXIMUM_DEPTH", max_depth, &
303 "The (diagnosed) maximum depth of the ocean.", &
304 units="m", unscale=us%Z_to_m, like_default=.true.)
305 if (trim(config) /= "DOME") then
306 ! MAXIMUM_DEPTH is not set and topography does not need to be trimmed by its maximum depth.
307 call limit_topography(d, g, pf, -max_depth_default * us%m_to_Z, us)
308 endif
309 endif
310
311end subroutine mom_initialize_topography
312