← back to index

src/initialization/MOM_fixed_initialization.F90

portedportable, not yet portedexecuted, not portableexecutable, not hit by this run

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Initializes fixed aspects of the model, such as horizontal grid metrics,
6!! topography and Coriolis.
7module MOM_fixed_initialization
8
9use MOM_debugging, only : hchksum, qchksum, uvchksum
10use MOM_domains, only : pass_var
11use MOM_dyn_horgrid, only : dyn_horgrid_type
12use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, WARNING, is_root_pe
13use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint
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
17use MOM_grid_initialize, only : initialize_masks, set_grid_metrics
18use MOM_open_boundary, only : ocean_OBC_type
19use MOM_open_boundary, only : open_boundary_config, open_boundary_query
20use MOM_open_boundary, only : open_boundary_impose_normal_slope
21use MOM_open_boundary, only : open_boundary_impose_land_mask
22use MOM_shared_initialization, only : MOM_initialize_rotation, MOM_calculate_grad_Coriolis
23use MOM_shared_initialization, only : initialize_topography_from_file, apply_topography_edits_from_file
24use MOM_shared_initialization, only : initialize_topography_named, limit_topography, diagnoseMaximumDepth
25use MOM_shared_initialization, only : set_rotation_planetary, set_rotation_beta_plane, initialize_grid_rotation_angle
26use MOM_shared_initialization, only : reset_face_lengths_named, reset_face_lengths_file, reset_face_lengths_list
27use MOM_shared_initialization, only : read_face_length_list, set_velocity_depth_max, set_velocity_depth_min
28use MOM_shared_initialization, only : set_subgrid_topo_at_vel_from_file
29use MOM_shared_initialization, only : compute_global_grid_integrals
30use MOM_shared_initialization, only : set_meanSL_from_file
31use MOM_unit_scaling, only : unit_scale_type
32
33use user_initialization, only : user_initialize_topography
34use DOME_initialization, only : DOME_initialize_topography
35use ISOMIP_initialization, only : ISOMIP_initialize_topography
36use basin_builder, only : basin_builder_topography
37use benchmark_initialization, only : benchmark_initialize_topography
38use Neverworld_initialization, only : Neverworld_initialize_topography
39use DOME2d_initialization, only : DOME2d_initialize_topography
40use Kelvin_initialization, only : Kelvin_initialize_topography
41use sloshing_initialization, only : sloshing_initialize_topography
42use seamount_initialization, only : seamount_initialize_topography
43use dumbbell_initialization, only : dumbbell_initialize_topography
44use shelfwave_initialization, only : shelfwave_initialize_topography
45use Phillips_initialization, only : Phillips_initialize_topography
46use dense_water_initialization, only : dense_water_initialize_topography
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.
571subroutine 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
741 call callTree_enter("MOM_initialize_fixed(), MOM_fixed_initialization.F90")
751 call get_param(PF, mdl, "DEBUG", debug, default=.false.)
76
77 ! Set up the parameters of the physical domain (i.e. the grid), G
781 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 "// &
831 "time mean ocean total thickness.", default=.False.)
841 if (read_meanSL_file) &
850 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(:,:) )
901 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
931 call pass_var(G%bathyT, G%Domain)
94
95 ! Determine the position of any open boundaries and create OBC
961 call open_boundary_config(G, US, PF, OBC)
97
98 ! Make bathymetry (if OBC_PROJECTION_BUG) and masks consistent with open boundaries.
991 if (associated(OBC)) then
100 call get_param(PF, mdl, "ENABLE_BUGS_BY_DEFAULT", enable_bugs, &
1010 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.", &
1090 default=enable_bugs)
1100 open_corners = .not.OBC_projection_bug
111
1120 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, "//&
1160 "OBC_PROJECTION_BUG needs to be False.")
117
118 ! This call sets masks that prohibit flow over any point interpreted as land
1190 if (OBC_projection_bug) &
1200 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, &
1220 open_corner_OBCs=open_corners)
123 ! Make OBC mask consistent with land mask
1240 call open_boundary_impose_land_mask(OBC, G, G%areaCu, G%areaCv, US)
125 else
1261 call initialize_masks(G, PF, US)
127 endif
128
1291 if (debug) then
1300 call hchksum(G%bathyT, 'MOM_initialize_fixed: depth ', G%HI, haloshift=1, unscale=US%Z_to_m)
1310 call hchksum(G%mask2dT, 'MOM_initialize_fixed: mask2dT ', G%HI)
132 call uvchksum('MOM_initialize_fixed: mask2dC[uv]', G%mask2dCu, &
1330 G%mask2dCv, G%HI)
1340 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".
1391 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.", &
1521 default="none")
1533 select case ( trim(config) )
154 case ("none")
1550 case ("list") ; call reset_face_lengths_list(G, PF, US)
1560 case ("file") ; call reset_face_lengths_file(G, PF, US)
1570 case ("global_1deg") ; call reset_face_lengths_named(G, PF, trim(config), US)
158 case default ; call MOM_error(FATAL, "MOM_initialize_fixed: "// &
1592 "Unrecognized channel configuration "//trim(config))
160 end select
161
162 ! This call sets the topography at velocity points.
1631 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'.", &
1670 default="max")
1680 select case ( trim(config) )
1690 case ("max") ; call set_velocity_depth_max(G)
1700 case ("min") ; call set_velocity_depth_min(G)
171 case default ; call MOM_error(FATAL, "MOM_initialize_fixed: "// &
1720 "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.", &
1801 default=.False.)
1811 if (read_porous_file) &
1820 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].
1861 call MOM_initialize_rotation(G%CoriolisBu, G, PF, US=US)
187 ! Calculate the components of grad f (beta)
1881 call MOM_calculate_grad_Coriolis(G%dF_dx, G%dF_dy, G, US=US)
189 ! Calculate the square of the Coriolis parameter
1909031 do I=G%IsdB,G%IedB ; do J=G%JsdB,G%JedB
1919030 G%Coriolis2Bu(I,J) = G%CoriolisBu(I,J)**2
192 enddo ; enddo
193
1941 if (debug) then
1950 call qchksum(G%CoriolisBu, "MOM_initialize_fixed: f ", G%HI, unscale=US%s_to_T)
1960 call qchksum(G%Coriolis2Bu, "MOM_initialize_fixed: f2 ", G%HI, unscale=US%s_to_T**2)
1970 call hchksum(G%dF_dx, "MOM_initialize_fixed: dF_dx ", G%HI, unscale=US%m_to_L*US%s_to_T)
1980 call hchksum(G%dF_dy, "MOM_initialize_fixed: dF_dy ", G%HI, unscale=US%m_to_L*US%s_to_T)
199 endif
200
2011 call initialize_grid_rotation_angle(G, PF)
202
203 ! Compute global integrals of grid values for later use in scalar diagnostics !
2041 call compute_global_grid_integrals(G, US=US)
205
2061 call callTree_leave('MOM_initialize_fixed()')
207
2081end subroutine MOM_initialize_fixed
209
210!> MOM_initialize_topography makes the appropriate call to set up the bathymetry in units of [Z ~> m].
2111subroutine 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
2302 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.", &
2621 fail_if_missing=.true.)
263 call get_param(PF, mdl, "MAXIMUM_DEPTH", max_depth, units="m", default=max_depth_default, &
2641 scale=US%m_to_Z, do_not_log=.true.)
2652 select case ( trim(config) )
2660 case ("file"); call initialize_topography_from_file(D, G, PF, US)
2670 case ("flat"); call initialize_topography_named(D, G, PF, config, max_depth, US)
2680 case ("spoon"); call initialize_topography_named(D, G, PF, config, max_depth, US)
2690 case ("bowl"); call initialize_topography_named(D, G, PF, config, max_depth, US)
2700 case ("halfpipe"); call initialize_topography_named(D, G, PF, config, max_depth, US)
2710 case ("DOME"); call DOME_initialize_topography(D, G, PF, max_depth, US)
2720 case ("ISOMIP"); call ISOMIP_initialize_topography(D, G, PF, max_depth, US)
2730 case ("bbuilder"); call basin_builder_topography(D, G, PF, max_depth)
2741 case ("benchmark"); call benchmark_initialize_topography(D, G, PF, max_depth, US)
2750 case ("Neverworld","Neverland"); call Neverworld_initialize_topography(D, G, PF, max_depth)
2760 case ("DOME2D"); call DOME2d_initialize_topography(D, G, PF, max_depth)
2770 case ("Kelvin"); call Kelvin_initialize_topography(D, G, PF, max_depth, US)
2780 case ("sloshing"); call sloshing_initialize_topography(D, G, PF, max_depth)
2790 case ("seamount"); call seamount_initialize_topography(D, G, PF, max_depth)
2800 case ("dumbbell"); call dumbbell_initialize_topography(D, G, PF, max_depth)
2810 case ("shelfwave"); call shelfwave_initialize_topography(D, G, PF, max_depth, US)
2820 case ("Phillips"); call Phillips_initialize_topography(D, G, PF, max_depth, US)
2830 case ("dense"); call dense_water_initialize_topography(D, G, PF, max_depth)
2840 case ("USER"); call user_initialize_topography(D, G, PF, max_depth, US)
285 case default ; call MOM_error(FATAL,"MOM_initialize_topography: "// &
2862 "Unrecognized topography setup '"//trim(config)//"'")
287 end select
2881 if (max_depth /= max_depth_default * US%m_to_Z) then
289 call log_param(PF, mdl, "MAXIMUM_DEPTH", max_depth, &
2901 "The maximum depth of the ocean.", units="m", unscale=US%Z_to_m)
2911 if (trim(config) /= "DOME") then
2921 call limit_topography(D, G, PF, max_depth, US)
293 endif
294 else
2950 if (present(meanSL)) then
2960 D_meanSL(:,:) = 0.0
2970 do j=G%jsc,G%jec ; do i=G%isc,G%iec ; D_meanSL(i,j) = D(i,j) + meanSL(i,j) ; enddo ; enddo
2980 max_depth = diagnoseMaximumDepth(D_meanSL, G)
299 else
3000 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.", &
3040 units="m", unscale=US%Z_to_m, like_default=.true.)
3050 if (trim(config) /= "DOME") then
306 ! MAXIMUM_DEPTH is not set and topography does not need to be trimmed by its maximum depth.
3070 call limit_topography(D, G, PF, -max_depth_default * US%m_to_Z, US)
308 endif
309 endif
310
3111end subroutine MOM_initialize_topography
312
313end module MOM_fixed_initialization