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