MOM6
benchmark_initialization.F90
Go to the documentation of this file.
1 !> Initialization for the "bench mark" configuration
3 
4 ! This file is part of MOM6. See LICENSE.md for the license.
5 
10 use mom_get_input, only : directories
11 use mom_grid, only : ocean_grid_type
17 
18 implicit none ; private
19 
20 #include <MOM_memory.h>
21 
25 
26 ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
27 ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
28 ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
29 ! vary with the Boussinesq approximation, the Boussinesq variant is given first.
30 
31 contains
32 
33 !> This subroutine sets up the benchmark test case topography.
34 subroutine benchmark_initialize_topography(D, G, param_file, max_depth, US)
35  type(dyn_horgrid_type), intent(in) :: g !< The dynamic horizontal grid type
36  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
37  intent(out) :: d !< Ocean bottom depth in m or [Z ~> m] if US is present
38  type(param_file_type), intent(in) :: param_file !< Parameter file structure
39  real, intent(in) :: max_depth !< Maximum model depth in the units of D
40  type(unit_scale_type), optional, intent(in) :: us !< A dimensional unit scaling type
41 
42  ! Local variables
43  real :: min_depth ! The minimum and maximum depths [Z ~> m].
44  real :: pi ! 3.1415926... calculated as 4*atan(1)
45  real :: d0 ! A constant to make the maximum !
46  ! basin depth MAXIMUM_DEPTH. !
47  real :: m_to_z ! A dimensional rescaling factor.
48  real :: x, y
49  ! This include declares and sets the variable "version".
50 # include "version_variable.h"
51  character(len=40) :: mdl = "benchmark_initialize_topography" ! This subroutine's name.
52  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
53  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
54  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
55 
56  call mom_mesg(" benchmark_initialization.F90, benchmark_initialize_topography: setting topography", 5)
57 
58  m_to_z = 1.0 ; if (present(us)) m_to_z = us%m_to_Z
59 
60  call log_version(param_file, mdl, version, "")
61  call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, &
62  "The minimum depth of the ocean.", units="m", default=0.0, scale=m_to_z)
63 
64  pi = 4.0*atan(1.0)
65  d0 = max_depth / 0.5
66 
67 ! Calculate the depth of the bottom.
68  do j=js,je ; do i=is,ie
69  x = (g%geoLonT(i,j)-g%west_lon) / g%len_lon
70  y = (g%geoLatT(i,j)-g%south_lat) / g%len_lat
71 ! This sets topography that has a reentrant channel to the south.
72  d(i,j) = -d0 * ( y*(1.0 + 0.6*cos(4.0*pi*x)) &
73  + 0.75*exp(-6.0*y) &
74  + 0.05*cos(10.0*pi*x) - 0.7 )
75  if (d(i,j) > max_depth) d(i,j) = max_depth
76  if (d(i,j) < min_depth) d(i,j) = 0.
77  enddo ; enddo
78 
80 
81 !> Initializes layer thicknesses for the benchmark test case,
82 !! by finding the depths of interfaces in a specified latitude-dependent
83 !! temperature profile with an exponentially decaying thermocline on top of a
84 !! linear stratification.
85 subroutine benchmark_initialize_thickness(h, G, GV, US, param_file, eqn_of_state, &
86  P_ref, just_read_params)
87  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
88  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
89  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
90  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
91  intent(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2].
92  type(param_file_type), intent(in) :: param_file !< A structure indicating the open file
93  !! to parse for model parameter values.
94  type(eos_type), pointer :: eqn_of_state !< integer that selects the
95  !! equation of state.
96  real, intent(in) :: p_ref !< The coordinate-density
97  !! reference pressure [Pa].
98  logical, optional, intent(in) :: just_read_params !< If present and true, this call will
99  !! only read parameters without changing h.
100  ! Local variables
101  real :: e0(szk_(gv)+1) ! The resting interface heights, in depth units [Z ~> m],
102  ! usually negative because it is positive upward.
103  real :: e_pert(szk_(gv)+1) ! Interface height perturbations, positive upward,
104  ! in depth units [Z ~> m].
105  real :: eta1d(szk_(gv)+1) ! Interface height relative to the sea surface
106  ! positive upward, in depth units [Z ~> m].
107  real :: sst ! The initial sea surface temperature [degC].
108  real :: t_int ! The initial temperature of an interface [degC].
109  real :: ml_depth ! The specified initial mixed layer depth, in depth units [Z ~> m].
110  real :: thermocline_scale ! The e-folding scale of the thermocline, in depth units [Z ~> m].
111  real, dimension(SZK_(GV)) :: &
112  t0, pres, s0, & ! drho
113  rho_guess, & ! Potential density at T0 & S0 [R ~> kg m-3].
114  drho_dt, & ! Derivative of density with temperature [R degC-1 ~> kg m-3 degC-1].
115  drho_ds ! Derivative of density with salinity [R ppt-1 ~> kg m-3 ppt-1].
116  real :: a_exp ! The fraction of the overall stratification that is exponential.
117  real :: i_ts, i_md ! Inverse lengthscales [Z-1 ~> m-1].
118  real :: t_frac ! A ratio of the interface temperature to the range
119  ! between SST and the bottom temperature.
120  real :: err, derr_dz ! The error between the profile's temperature and the
121  ! interface temperature for a given z and its derivative.
122  real :: pi, z
123  logical :: just_read
124  ! This include declares and sets the variable "version".
125 # include "version_variable.h"
126  character(len=40) :: mdl = "benchmark_initialize_thickness" ! This subroutine's name.
127  integer :: i, j, k, k1, is, ie, js, je, nz, itt
128 
129  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
130 
131  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params
132  if (.not.just_read) call log_version(param_file, mdl, version, "")
133  call get_param(param_file, mdl, "BENCHMARK_ML_DEPTH_IC", ml_depth, &
134  "Initial mixed layer depth in the benchmark test case.", &
135  units='m', default=50.0, scale=us%m_to_Z, do_not_log=just_read)
136  call get_param(param_file, mdl, "BENCHMARK_THERMOCLINE_SCALE", thermocline_scale, &
137  "Initial thermocline depth scale in the benchmark test case.", &
138  default=500.0, units="m", scale=us%m_to_Z, do_not_log=just_read)
139 
140  if (just_read) return ! This subroutine has no run-time parameters.
141 
142  call mom_mesg(" benchmark_initialization.F90, benchmark_initialize_thickness: setting thickness", 5)
143 
144  k1 = gv%nk_rho_varies + 1
145 
146  a_exp = 0.9
147 
148 ! This block calculates T0(k) for the purpose of diagnosing where the
149 ! interfaces will be found.
150  do k=1,nz
151  pres(k) = p_ref ; s0(k) = 35.0
152  enddo
153  t0(k1) = 29.0
154  call calculate_density(t0(k1), s0(k1), pres(k1), rho_guess(k1), eqn_of_state, scale=us%kg_m3_to_R)
155  call calculate_density_derivs(t0, s0, pres, drho_dt, drho_ds, k1, 1, eqn_of_state, scale=us%kg_m3_to_R)
156 
157 ! A first guess of the layers' temperatures.
158  do k=1,nz
159  t0(k) = t0(k1) + (gv%Rlay(k) - rho_guess(k1)) / drho_dt(k1)
160  enddo
161 
162 ! Refine the guesses for each layer.
163  do itt=1,6
164  call calculate_density(t0, s0, pres, rho_guess, 1, nz, eqn_of_state, scale=us%kg_m3_to_R)
165  call calculate_density_derivs(t0, s0, pres, drho_dt, drho_ds, 1, nz, eqn_of_state, scale=us%kg_m3_to_R)
166  do k=1,nz
167  t0(k) = t0(k) + (gv%Rlay(k) - rho_guess(k)) / drho_dt(k)
168  enddo
169  enddo
170 
171  pi = 4.0*atan(1.0)
172  i_ts = 1.0 / thermocline_scale
173  i_md = 1.0 / g%max_depth
174  do j=js,je ; do i=is,ie
175  sst = 0.5*(t0(k1)+t0(nz)) - 0.9*0.5*(t0(k1)-t0(nz)) * &
176  cos(pi*(g%geoLatT(i,j)-g%south_lat)/(g%len_lat))
177 
178  do k=1,nz ; e_pert(k) = 0.0 ; enddo
179 
180  ! This sets the initial thickness (in [H ~> m or kg m-2]) of the layers. The thicknesses
181  ! are set to insure that: 1. each layer is at least Gv%Angstrom_m thick, and
182  ! 2. the interfaces are where they should be based on the resting depths and interface
183  ! height perturbations, as long at this doesn't interfere with 1.
184  eta1d(nz+1) = -g%bathyT(i,j)
185 
186  do k=nz,2,-1
187  t_int = 0.5*(t0(k) + t0(k-1))
188  t_frac = (t_int - t0(nz)) / (sst - t0(nz))
189  ! Find the z such that T_frac = a exp(z/thermocline_scale) + (1-a) (z+D)/D
190  z = 0.0
191  do itt=1,6
192  err = a_exp * exp(z*i_ts) + (1.0 - a_exp) * (z*i_md + 1.0) - t_frac
193  derr_dz = a_exp * i_ts * exp(z*i_ts) + (1.0 - a_exp) * i_md
194  z = z - err / derr_dz
195  enddo
196  e0(k) = z
197 ! e0(K) = -ML_depth + thermocline_scale * log((T_int - T0(nz)) / (SST - T0(nz)))
198 
199  eta1d(k) = e0(k) + e_pert(k)
200 
201  if (eta1d(k) > -ml_depth) eta1d(k) = -ml_depth
202 
203  if (eta1d(k) < eta1d(k+1) + gv%Angstrom_Z) &
204  eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
205 
206  h(i,j,k) = max(gv%Z_to_H * (eta1d(k) - eta1d(k+1)), gv%Angstrom_H)
207  enddo
208  h(i,j,1) = max(gv%Z_to_H * (0.0 - eta1d(2)), gv%Angstrom_H)
209 
210  enddo ; enddo
211 
212 end subroutine benchmark_initialize_thickness
213 
214 !> Initializes layer temperatures and salinities for benchmark
215 subroutine benchmark_init_temperature_salinity(T, S, G, GV, US, param_file, &
216  eqn_of_state, P_Ref, just_read_params)
217  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
218  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
219  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: t !< The potential temperature
220  !! that is being initialized.
221  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: s !< The salinity that is being
222  !! initialized.
223  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
224  type(param_file_type), intent(in) :: param_file !< A structure indicating the
225  !! open file to parse for
226  !! model parameter values.
227  type(eos_type), pointer :: eqn_of_state !< integer that selects the
228  !! equation of state.
229  real, intent(in) :: p_ref !< The coordinate-density
230  !! reference pressure [Pa].
231  logical, optional, intent(in) :: just_read_params !< If present and true, this call will
232  !! only read parameters without changing h.
233  ! Local variables
234  real :: t0(szk_(g)), s0(szk_(g))
235  real :: pres(szk_(g)) ! Reference pressure [kg m-3].
236  real :: drho_dt(szk_(g)) ! Derivative of density with temperature [R degC-1 ~> kg m-3 degC-1].
237  real :: drho_ds(szk_(g)) ! Derivative of density with salinity [R ppt-1 ~> kg m-3 ppt-1].
238  real :: rho_guess(szk_(g)) ! Potential density at T0 & S0 [R ~> kg m-3].
239  real :: pi ! 3.1415926... calculated as 4*atan(1)
240  real :: sst ! The initial sea surface temperature [degC].
241  real :: lat
242  logical :: just_read ! If true, just read parameters but set nothing.
243  character(len=40) :: mdl = "benchmark_init_temperature_salinity" ! This subroutine's name.
244  integer :: i, j, k, k1, is, ie, js, je, nz, itt
245 
246  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
247 
248  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params
249 
250  if (just_read) return ! All run-time parameters have been read, so return.
251 
252  k1 = gv%nk_rho_varies + 1
253 
254  do k=1,nz
255  pres(k) = p_ref ; s0(k) = 35.0
256  enddo
257 
258  t0(k1) = 29.0
259  call calculate_density(t0(k1),s0(k1),pres(k1),rho_guess(k1),eqn_of_state, scale=us%kg_m3_to_R)
260  call calculate_density_derivs(t0,s0,pres,drho_dt,drho_ds,k1,1,eqn_of_state, scale=us%kg_m3_to_R)
261 
262 ! A first guess of the layers' temperatures. !
263  do k=1,nz
264  t0(k) = t0(k1) + (gv%Rlay(k) - rho_guess(k1)) / drho_dt(k1)
265  enddo
266 
267 ! Refine the guesses for each layer. !
268  do itt = 1,6
269  call calculate_density(t0,s0,pres,rho_guess,1,nz,eqn_of_state, scale=us%kg_m3_to_R)
270  call calculate_density_derivs(t0,s0,pres,drho_dt,drho_ds,1,nz,eqn_of_state, scale=us%kg_m3_to_R)
271  do k=1,nz
272  t0(k) = t0(k) + (gv%Rlay(k) - rho_guess(k)) / drho_dt(k)
273  enddo
274  enddo
275 
276  do k=1,nz ; do i=is,ie ; do j=js,je
277  t(i,j,k) = t0(k)
278  s(i,j,k) = s0(k)
279  enddo ; enddo ; enddo
280  pi = 4.0*atan(1.0)
281  do i=is,ie ; do j=js,je
282  sst = 0.5*(t0(k1)+t0(nz)) - 0.9*0.5*(t0(k1)-t0(nz)) * &
283  cos(pi*(g%geoLatT(i,j)-g%south_lat)/(g%len_lat))
284  do k=1,k1-1
285  t(i,j,k) = sst
286  enddo
287  enddo ; enddo
288 
290 
291 end module benchmark_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_file_parser::log_version
An overloaded interface to log version information about modules.
Definition: MOM_file_parser.F90:109
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
benchmark_initialization::benchmark_initialize_topography
subroutine, public benchmark_initialize_topography(D, G, param_file, max_depth, US)
This subroutine sets up the benchmark test case topography.
Definition: benchmark_initialization.F90:35
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
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
mom_eos::eos_type
A control structure for the equation of state.
Definition: MOM_EOS.F90:91
benchmark_initialization::benchmark_init_temperature_salinity
subroutine, public benchmark_init_temperature_salinity(T, S, G, GV, US, param_file, eqn_of_state, P_Ref, just_read_params)
Initializes layer temperatures and salinities for benchmark.
Definition: benchmark_initialization.F90:217
mom_verticalgrid::verticalgrid_type
Describes the vertical ocean grid, including unit conversion factors.
Definition: MOM_verticalGrid.F90:24
mom_variables
Provides transparent structures with groups of MOM6 variables and supporting routines.
Definition: MOM_variables.F90:2
benchmark_initialization
Initialization for the "bench mark" configuration.
Definition: benchmark_initialization.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
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
benchmark_initialization::benchmark_initialize_thickness
subroutine, public benchmark_initialize_thickness(h, G, GV, US, param_file, eqn_of_state, P_ref, just_read_params)
Initializes layer thicknesses for the benchmark test case, by finding the depths of interfaces in a s...
Definition: benchmark_initialization.F90:87
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
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
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
mom_eos::calculate_density
Calculates density of sea water from T, S and P.
Definition: MOM_EOS.F90:60