MOM6
seamount_initialization.F90
Go to the documentation of this file.
1 !> Configures the model for the idealized seamount test case.
3 
4 ! This file is part of MOM6. See LICENSE.md for the license.
5 
6 use mom_domains, only : sum_across_pes
10 use mom_get_input, only : directories
11 use mom_grid, only : ocean_grid_type
21 
22 implicit none ; private
23 
24 #include <MOM_memory.h>
25 
26 character(len=40) :: mdl = "seamount_initialization" !< This module's name.
27 
28 ! The following routines are visible to the outside world
32 
33 ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
34 ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
35 ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
36 ! vary with the Boussinesq approximation, the Boussinesq variant is given first.
37 
38 contains
39 
40 !> Initialization of topography.
41 subroutine seamount_initialize_topography( D, G, param_file, max_depth )
42  type(dyn_horgrid_type), intent(in) :: g !< The dynamic horizontal grid type
43  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
44  intent(out) :: d !< Ocean bottom depth in the units of depth_max
45  type(param_file_type), intent(in) :: param_file !< Parameter file structure
46  real, intent(in) :: max_depth !< Maximum ocean depth in arbitrary units
47 
48  ! Local variables
49  integer :: i, j
50  real :: x, y, delta, lx, rlx, ly, rly
51 
52  call get_param(param_file, mdl,"SEAMOUNT_DELTA",delta, &
53  "Non-dimensional height of seamount.", &
54  units="non-dim", default=0.5)
55  call get_param(param_file, mdl,"SEAMOUNT_X_LENGTH_SCALE",lx, &
56  "Length scale of seamount in x-direction. "//&
57  "Set to zero make topography uniform in the x-direction.", &
58  units="Same as x,y", default=20.)
59  call get_param(param_file, mdl,"SEAMOUNT_Y_LENGTH_SCALE",ly, &
60  "Length scale of seamount in y-direction. "//&
61  "Set to zero make topography uniform in the y-direction.", &
62  units="Same as x,y", default=0.)
63 
64  lx = lx / g%len_lon
65  ly = ly / g%len_lat
66  rlx = 0. ; if (lx>0.) rlx = 1. / lx
67  rly = 0. ; if (ly>0.) rly = 1. / ly
68 
69  do j=g%jsc,g%jec ; do i=g%isc,g%iec
70  ! Compute normalized zonal coordinates (x,y=0 at center of domain)
71  x = ( g%geoLonT(i,j) - g%west_lon ) / g%len_lon - 0.5
72  y = ( g%geoLatT(i,j) - g%south_lat ) / g%len_lat - 0.5
73  d(i,j) = g%max_depth * ( 1.0 - delta * exp(-(rlx*x)**2 -(rly*y)**2) )
74  enddo ; enddo
75 
76 end subroutine seamount_initialize_topography
77 
78 !> Initialization of thicknesses.
79 !! This subroutine initializes the layer thicknesses to be uniform.
80 subroutine seamount_initialize_thickness ( h, G, GV, US, param_file, just_read_params)
81  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
82  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
83  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
84  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
85  intent(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2].
86  type(param_file_type), intent(in) :: param_file !< A structure indicating the open file
87  !! to parse for model parameter values.
88  logical, optional, intent(in) :: just_read_params !< If present and true, this call will
89  !! only read parameters without changing h.
90 
91  real :: e0(szk_(g)+1) ! The resting interface heights [Z ~> m], usually
92  ! negative because it is positive upward.
93  real :: eta1d(szk_(g)+1)! Interface height relative to the sea surface
94  ! positive upward [Z ~> m].
95  real :: min_thickness ! The minimum layer thicknesses [Z ~> m].
96  real :: s_surf, s_range, s_ref, s_light, s_dense ! Various salinities [ppt].
97  real :: eta_ic_quanta ! The granularity of quantization of intial interface heights [Z-1 ~> m-1].
98  character(len=20) :: verticalcoordinate
99  logical :: just_read ! If true, just read parameters but set nothing.
100  integer :: i, j, k, is, ie, js, je, nz
101 
102  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
103 
104  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params
105 
106  if (.not.just_read) &
107  call mom_mesg("MOM_initialization.F90, initialize_thickness_uniform: setting thickness")
108 
109  call get_param(param_file, mdl,"MIN_THICKNESS",min_thickness, &
110  'Minimum thickness for layer', &
111  units='m', default=1.0e-3, do_not_log=just_read, scale=us%m_to_Z)
112  call get_param(param_file, mdl,"REGRIDDING_COORDINATE_MODE",verticalcoordinate, &
113  default=default_coordinate_mode, do_not_log=just_read)
114 
115  ! WARNING: this routine specifies the interface heights so that the last layer
116  ! is vanished, even at maximum depth. In order to have a uniform
117  ! layer distribution, use this line of code within the loop:
118  ! e0(k) = -G%max_depth * real(k-1) / real(nz)
119  ! To obtain a thickness distribution where the last layer is
120  ! vanished and the other thicknesses uniformly distributed, use:
121  ! e0(k) = -G%max_depth * real(k-1) / real(nz-1)
122  !do k=1,nz+1
123  ! e0(k) = -G%max_depth * real(k-1) / real(nz)
124  !enddo
125 
126  select case ( coordinatemode(verticalcoordinate) )
127 
128  case ( regridding_layer, regridding_rho ) ! Initial thicknesses for isopycnal coordinates
129  call get_param(param_file, mdl,"INITIAL_SSS", s_surf, default=34., do_not_log=.true.)
130  call get_param(param_file, mdl,"INITIAL_S_RANGE", s_range, default=2., do_not_log=.true.)
131  call get_param(param_file, mdl, "S_REF", s_ref, default=35.0, do_not_log=.true.)
132  call get_param(param_file, mdl, "TS_RANGE_S_LIGHT", s_light, default = s_ref, do_not_log=.true.)
133  call get_param(param_file, mdl, "TS_RANGE_S_DENSE", s_dense, default = s_ref, do_not_log=.true.)
134  call get_param(param_file, mdl, "INTERFACE_IC_QUANTA", eta_ic_quanta, &
135  "The granularity of initial interface height values "//&
136  "per meter, to avoid sensivity to order-of-arithmetic changes.", &
137  default=2048.0, units="m-1", scale=us%Z_to_m, do_not_log=just_read)
138  if (just_read) return ! All run-time parameters have been read, so return.
139 
140  do k=1,nz+1
141  ! Salinity of layer k is S_light + (k-1)/(nz-1) * (S_dense - S_light)
142  ! Salinity of interface K is S_light + (K-3/2)/(nz-1) * (S_dense - S_light)
143  ! Salinity at depth z should be S(z) = S_surf - S_range * z/max_depth
144  ! Equating: S_surf - S_range * z/max_depth = S_light + (K-3/2)/(nz-1) * (S_dense - S_light)
145  ! Equating: - S_range * z/max_depth = S_light - S_surf + (K-3/2)/(nz-1) * (S_dense - S_light)
146  ! Equating: z/max_depth = - ( S_light - S_surf + (K-3/2)/(nz-1) * (S_dense - S_light) ) / S_range
147  e0(k) = - g%max_depth * ( ( s_light - s_surf ) + ( s_dense - s_light ) * &
148  ( (real(k)-1.5) / real(nz-1) ) ) / s_range
149  ! Force round numbers ... the above expression has irrational factors ...
150  if (eta_ic_quanta > 0.0) &
151  e0(k) = nint(eta_ic_quanta*e0(k)) / eta_ic_quanta
152  e0(k) = min(real(1-k)*gv%Angstrom_Z, e0(k)) ! Bound by surface
153  e0(k) = max(-g%max_depth, e0(k)) ! Bound by bottom
154  enddo
155  do j=js,je ; do i=is,ie
156  eta1d(nz+1) = -g%bathyT(i,j)
157  do k=nz,1,-1
158  eta1d(k) = e0(k)
159  if (eta1d(k) < (eta1d(k+1) + gv%Angstrom_Z)) then
160  eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
161  h(i,j,k) = gv%Angstrom_H
162  else
163  h(i,j,k) = gv%Z_to_H * (eta1d(k) - eta1d(k+1))
164  endif
165  enddo
166  enddo ; enddo
167 
168  case ( regridding_zstar ) ! Initial thicknesses for z coordinates
169  if (just_read) return ! All run-time parameters have been read, so return.
170  do j=js,je ; do i=is,ie
171  eta1d(nz+1) = -g%bathyT(i,j)
172  do k=nz,1,-1
173  eta1d(k) = -g%max_depth * real(k-1) / real(nz)
174  if (eta1d(k) < (eta1d(k+1) + min_thickness)) then
175  eta1d(k) = eta1d(k+1) + min_thickness
176  h(i,j,k) = gv%Z_to_H * min_thickness
177  else
178  h(i,j,k) = gv%Z_to_H * (eta1d(k) - eta1d(k+1))
179  endif
180  enddo
181  enddo ; enddo
182 
183  case ( regridding_sigma ) ! Initial thicknesses for sigma coordinates
184  if (just_read) return ! All run-time parameters have been read, so return.
185  do j=js,je ; do i=is,ie
186  h(i,j,:) = gv%Z_to_H * g%bathyT(i,j) / dfloat(nz)
187  enddo ; enddo
188 
189 end select
190 
191 end subroutine seamount_initialize_thickness
192 
193 !> Initial values for temperature and salinity
194 subroutine seamount_initialize_temperature_salinity ( T, S, h, G, GV, param_file, &
195  eqn_of_state, just_read_params)
196  type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
197  type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
198  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: t !< Potential temperature [degC]
199  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: s !< Salinity [ppt]
200  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
201  type(param_file_type), intent(in) :: param_file !< Parameter file structure
202  type(eos_type), pointer :: eqn_of_state !< Equation of state structure
203  logical, optional, intent(in) :: just_read_params !< If present and true, this call will
204  !! only read parameters without changing h.
205 
206  ! Local variables
207  integer :: i, j, k, is, ie, js, je, nz, k_light
208  real :: xi0, xi1, dxi, r, s_surf, t_surf, s_range, t_range
209  real :: t_ref, t_light, t_dense, s_ref, s_light, s_dense, a1, frac_dense, k_frac, res_rat
210  logical :: just_read ! If true, just read parameters but set nothing.
211  character(len=20) :: verticalcoordinate, density_profile
212 
213  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
214 
215  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params
216 
217  call get_param(param_file, mdl, "REGRIDDING_COORDINATE_MODE", verticalcoordinate, &
218  default=default_coordinate_mode, do_not_log=just_read)
219  call get_param(param_file, mdl,"INITIAL_DENSITY_PROFILE", density_profile, &
220  'Initial profile shape. Valid values are "linear", "parabolic" '//&
221  'and "exponential".', default='linear', do_not_log=just_read)
222  call get_param(param_file, mdl,"INITIAL_SSS", s_surf, &
223  'Initial surface salinity', units='1e-3', default=34., do_not_log=just_read)
224  call get_param(param_file, mdl,"INITIAL_SST", t_surf, &
225  'Initial surface temperature', units='C', default=0., do_not_log=just_read)
226  call get_param(param_file, mdl,"INITIAL_S_RANGE", s_range, &
227  'Initial salinity range (bottom - surface)', units='1e-3', &
228  default=2., do_not_log=just_read)
229  call get_param(param_file, mdl,"INITIAL_T_RANGE", t_range, &
230  'Initial temperature range (bottom - surface)', units='C', &
231  default=0., do_not_log=just_read)
232 
233  select case ( coordinatemode(verticalcoordinate) )
234  case ( regridding_layer ) ! Initial thicknesses for layer isopycnal coordinates
235  ! These parameters are used in MOM_fixed_initialization.F90 when CONFIG_COORD="ts_range"
236  call get_param(param_file, mdl, "T_REF", t_ref, default=10.0, do_not_log=.true.)
237  call get_param(param_file, mdl, "TS_RANGE_T_LIGHT", t_light, default=t_ref, do_not_log=.true.)
238  call get_param(param_file, mdl, "TS_RANGE_T_DENSE", t_dense, default=t_ref, do_not_log=.true.)
239  call get_param(param_file, mdl, "S_REF", s_ref, default=35.0, do_not_log=.true.)
240  call get_param(param_file, mdl, "TS_RANGE_S_LIGHT", s_light, default = s_ref, do_not_log=.true.)
241  call get_param(param_file, mdl, "TS_RANGE_S_DENSE", s_dense, default = s_ref, do_not_log=.true.)
242  call get_param(param_file, mdl, "TS_RANGE_RESOLN_RATIO", res_rat, default=1.0, do_not_log=.true.)
243  if (just_read) return ! All run-time parameters have been read, so return.
244 
245  ! Emulate the T,S used in the "ts_range" coordinate configuration code
246  k_light = gv%nk_rho_varies + 1
247  do j=js,je ; do i=is,ie
248  t(i,j,k_light) = t_light ; s(i,j,k_light) = s_light
249  enddo ; enddo
250  a1 = 2.0 * res_rat / (1.0 + res_rat)
251  do k=k_light+1,nz
252  k_frac = real(k-k_light)/real(nz-k_light)
253  frac_dense = a1 * k_frac + (1.0 - a1) * k_frac**2
254  do j=js,je ; do i=is,ie
255  t(i,j,k) = frac_dense * (t_dense - t_light) + t_light
256  s(i,j,k) = frac_dense * (s_dense - s_light) + s_light
257  enddo ; enddo
258  enddo
259  case ( regridding_sigma, regridding_zstar, regridding_rho ) ! All other coordinate use FV initialization
260  if (just_read) return ! All run-time parameters have been read, so return.
261  do j=js,je ; do i=is,ie
262  xi0 = 0.0
263  do k = 1,nz
264  xi1 = xi0 + gv%H_to_Z * h(i,j,k) / g%max_depth
265  select case ( trim(density_profile) )
266  case ('linear')
267  !S(i,j,k) = S_surf + S_range * 0.5 * (xi0 + xi1)
268  s(i,j,k) = s_surf + ( 0.5 * s_range ) * (xi0 + xi1) ! Coded this way to reproduce old hard-coded answers
269  t(i,j,k) = t_surf + t_range * 0.5 * (xi0 + xi1)
270  case ('parabolic')
271  s(i,j,k) = s_surf + s_range * (2.0 / 3.0) * (xi1**3 - xi0**3) / (xi1 - xi0)
272  t(i,j,k) = t_surf + t_range * (2.0 / 3.0) * (xi1**3 - xi0**3) / (xi1 - xi0)
273  case ('exponential')
274  r = 0.8 ! small values give sharp profiles
275  s(i,j,k) = s_surf + s_range * (exp(xi1/r)-exp(xi0/r)) / (xi1 - xi0)
276  t(i,j,k) = t_surf + t_range * (exp(xi1/r)-exp(xi0/r)) / (xi1 - xi0)
277  case default
278  call mom_error(fatal, 'Unknown value for "INITIAL_DENSITY_PROFILE"')
279  end select
280  xi0 = xi1
281  enddo
282  enddo ; enddo
283  end select
284 
286 
287 end module seamount_initialization
mom_verticalgrid
Provides a transparent vertical ocean grid type and supporting routines.
Definition: MOM_verticalGrid.F90:2
mom_sponge::set_up_sponge_field
subroutine, public set_up_sponge_field(sp_val, f_ptr, G, nlay, CS, sp_val_i_mean)
This subroutine stores the reference profile for the variable whose address is given by f_ptr....
Definition: MOM_sponge.F90:214
mom_variables::thermo_var_ptrs
Pointers to an assortment of thermodynamic fields that may be available, including potential temperat...
Definition: MOM_variables.F90:78
mom_get_input::directories
Container for paths and parameter file names.
Definition: MOM_get_input.F90:20
mom_dyn_horgrid
Contains a shareable dynamic type for describing horizontal grids and metric data and utilty routines...
Definition: MOM_dyn_horgrid.F90:3
mom_error_handler::mom_mesg
subroutine, public mom_mesg(message, verb, all_print)
This provides a convenient interface for writing an informative comment.
Definition: MOM_error_handler.F90:53
mom_eos
Provides subroutines for quantities specific to the equation of state.
Definition: MOM_EOS.F90:2
mom_tracer_registry
This module contains the tracer_registry_type and the subroutines that handle registration of tracers...
Definition: MOM_tracer_registry.F90:5
mom_file_parser::param_file_type
A structure that can be parsed to read and document run-time parameters.
Definition: MOM_file_parser.F90:54
mom_file_parser::get_param
An overloaded interface to read and log the values of various types of parameters.
Definition: MOM_file_parser.F90:102
seamount_initialization::seamount_initialize_thickness
subroutine, public seamount_initialize_thickness(h, G, GV, US, param_file, just_read_params)
Initialization of thicknesses. This subroutine initializes the layer thicknesses to be uniform.
Definition: seamount_initialization.F90:81
regrid_consts::regridding_layer
integer, parameter regridding_layer
Layer mode identifier.
Definition: regrid_consts.F90:13
mom_get_input
Reads the only Fortran name list needed to boot-strap the model.
Definition: MOM_get_input.F90:6
mom_unit_scaling::unit_scale_type
Describes various unit conversion factors.
Definition: MOM_unit_scaling.F90:14
regrid_consts
Contains constants for interpreting input parameters that control regridding.
Definition: regrid_consts.F90:2
mom_eos::eos_type
A control structure for the equation of state.
Definition: MOM_EOS.F90:91
seamount_initialization::seamount_initialize_temperature_salinity
subroutine, public seamount_initialize_temperature_salinity(T, S, h, G, GV, param_file, eqn_of_state, just_read_params)
Initial values for temperature and salinity.
Definition: seamount_initialization.F90:196
seamount_initialization::seamount_initialize_topography
subroutine, public seamount_initialize_topography(D, G, param_file, max_depth)
Initialization of topography.
Definition: seamount_initialization.F90:42
regrid_consts::regridding_sigma
integer, parameter regridding_sigma
Sigma coordinates identifier.
Definition: regrid_consts.F90:16
mom_verticalgrid::verticalgrid_type
Describes the vertical ocean grid, including unit conversion factors.
Definition: MOM_verticalGrid.F90:24
mom_domains
Describes the decomposed MOM domain and has routines for communications across PEs.
Definition: MOM_domains.F90:2
mom_variables
Provides transparent structures with groups of MOM6 variables and supporting routines.
Definition: MOM_variables.F90:2
mom_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
mom_sponge
Implements sponge regions in isopycnal mode.
Definition: MOM_sponge.F90:2
regrid_consts::regridding_zstar
integer, parameter regridding_zstar
z* coordinates identifier
Definition: regrid_consts.F90:14
mom_eos::calculate_density_derivs
Calculate the derivatives of density with temperature and salinity from T, S, and P.
Definition: MOM_EOS.F90:70
mom_tracer_registry::tracer_registry_type
Type to carry basic tracer information.
Definition: MOM_tracer_registry.F90:138
mom_grid
Provides the ocean grid type.
Definition: MOM_grid.F90:2
mom_sponge::sponge_cs
This control structure holds memory and parameters for the MOM_sponge module.
Definition: MOM_sponge.F90:41
mom_unit_scaling
Provides a transparent unit rescaling type to facilitate dimensional consistency testing.
Definition: MOM_unit_scaling.F90:2
mom_error_handler::is_root_pe
logical function, public is_root_pe()
This returns .true. if the current PE is the root PE.
Definition: MOM_error_handler.F90:44
regrid_consts::coordinatemode
integer function coordinatemode(string)
Parse a string parameter specifying the coordinate mode and return the appropriate enumerated integer...
Definition: regrid_consts.F90:54
regrid_consts::default_coordinate_mode
character(len= *), parameter default_coordinate_mode
Default coordinate mode.
Definition: regrid_consts.F90:35
mom_error_handler::mom_error
subroutine, public mom_error(level, message, all_print)
This provides a convenient interface for writing an mpp_error message with run-time filter based on a...
Definition: MOM_error_handler.F90:72
regrid_consts::regridding_rho
integer, parameter regridding_rho
Density coordinates identifier.
Definition: regrid_consts.F90:15
mom_sponge::initialize_sponge
subroutine, public initialize_sponge(Iresttime, int_height, G, param_file, CS, GV, Iresttime_i_mean, int_height_i_mean)
This subroutine determines the number of points which are within sponges in this computational domain...
Definition: MOM_sponge.F90:90
seamount_initialization
Configures the model for the idealized seamount test case.
Definition: seamount_initialization.F90:2
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2
mom_dyn_horgrid::dyn_horgrid_type
Describes the horizontal ocean grid with only dynamic memory arrays.
Definition: MOM_dyn_horgrid.F90:23
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:26
seamount_initialization::mdl
character(len=40) mdl
This module's name.
Definition: seamount_initialization.F90:26
mom_eos::calculate_density
Calculates density of sea water from T, S and P.
Definition: MOM_EOS.F90:60