23 implicit none ;
private
25 #include <MOM_memory.h>
27 character(len=40) ::
mdl =
"dumbbell_initialization"
44 real,
dimension(G%isd:G%ied,G%jsd:G%jed), &
47 real,
intent(in) :: max_depth
51 real :: x, y, delta, dblen, dbfrac
54 'Lateral Length scale for dumbbell.',&
55 units=
'k', default=600., do_not_log=.false.)
56 call get_param(param_file,
mdl,
"DUMBBELL_FRACTION",dbfrac, &
57 'Meridional fraction for narrow part of dumbbell.',&
58 units=
'nondim', default=0.5, do_not_log=.false.)
60 if (g%x_axis_units ==
'm')
then
64 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
66 x = ( g%geoLonT(i,j) ) / dblen
67 y = ( g%geoLatT(i,j) ) / g%len_lat
69 if ((x>=-0.25 .and. x<=0.25) .and. (y <= -0.5*dbfrac .or. y >= 0.5*dbfrac))
then
81 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
85 logical,
optional,
intent(in) :: just_read_params
90 real :: eta1d(szk_(g)+1)
93 real :: s_surf, s_range, s_ref, s_light, s_dense
96 # include "version_variable.h"
97 character(len=20) :: verticalcoordinate
99 integer :: i, j, k, is, ie, js, je, nz
101 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
103 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
105 if (.not.just_read) &
106 call mom_mesg(
"MOM_initialization.F90, initialize_thickness_uniform: setting thickness")
108 if (.not.just_read)
call log_version(param_file,
mdl, version,
"")
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)
126 select case ( coordinatemode(verticalcoordinate) )
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
147 e0(k) = - g%max_depth * ( ( s_light - s_surf ) + ( s_dense - s_light ) * &
148 ( (real(k)-1.5) / real(nz-1) ) ) / s_range
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))
153 e0(k) = max(-g%max_depth, e0(k))
155 do j=js,je ;
do i=is,ie
156 eta1d(nz+1) = -g%bathyT(i,j)
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
163 h(i,j,k) = gv%Z_to_H * (eta1d(k) - eta1d(k+1))
168 case ( regridding_zstar )
169 if (just_read)
return
170 do j=js,je ;
do i=is,ie
171 eta1d(nz+1) = -g%bathyT(i,j)
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
178 h(i,j,k) = gv%Z_to_H * (eta1d(k) - eta1d(k+1))
184 if (just_read)
return
185 do j=js,je ;
do i=is,ie
186 h(i,j,:) = gv%Z_to_H * g%bathyT(i,j) / dfloat(nz)
195 eqn_of_state, just_read_params)
198 real,
dimension(SZI_(G),SZJ_(G), SZK_(G)),
intent(out) :: t
199 real,
dimension(SZI_(G),SZJ_(G), SZK_(G)),
intent(out) :: s
200 real,
dimension(SZI_(G),SZJ_(G), SZK_(G)),
intent(in) :: h
202 type(
eos_type),
pointer :: eqn_of_state
203 logical,
optional,
intent(in) :: just_read_params
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
210 real :: t_ref, t_light, t_dense, s_ref, s_light, s_dense, a1, frac_dense, k_frac, res_rat
212 character(len=20) :: verticalcoordinate, density_profile
214 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
216 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
220 call get_param(param_file,
mdl,
"REGRIDDING_COORDINATE_MODE", verticalcoordinate, &
221 default=default_coordinate_mode, do_not_log=just_read)
222 call get_param(param_file,
mdl,
"INITIAL_DENSITY_PROFILE", density_profile, &
223 'Initial profile shape. Valid values are "linear", "parabolic" '// &
224 'and "exponential".', default=
'linear', do_not_log=just_read)
225 call get_param(param_file,
mdl,
"DUMBBELL_SREF", s_surf, &
226 'DUMBBELL REFERENCE SALINITY', units=
'1e-3', default=34., do_not_log=just_read)
227 call get_param(param_file,
mdl,
"DUMBBELL_S_RANGE", s_range, &
228 'DUMBBELL salinity range (right-left)', units=
'1e-3', &
229 default=2., do_not_log=just_read)
231 'Lateral Length scale for dumbbell ',&
232 units=
'k', default=600., do_not_log=just_read)
234 if (g%x_axis_units ==
'm')
then
241 x = ( g%geoLonT(i,j) ) / dblen
247 s(i,j,k)=s_surf + 0.5*s_range
252 s(i,j,k)=s_surf - 0.5*s_range
268 logical,
intent(in) :: use_ale
272 real :: sponge_time_scale
274 real,
dimension(SZI_(G),SZJ_(G)) :: idamp
275 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h, t, s
276 real,
dimension(SZK_(GV)+1) :: e0, eta1d
278 integer :: i, j, k, nz
279 real :: x, zi, zmid, dist, min_thickness, dblen
280 real :: mld, s_ref, s_range, s_dense, t_ref, sill_height
282 'Lateral Length scale for dumbbell ',&
283 units=
'k', default=600., do_not_log=.true.)
285 if (g%x_axis_units ==
'm')
then
291 call get_param(param_file,
mdl,
"DUMBBELL_SPONGE_TIME_SCALE", sponge_time_scale, &
292 "The time scale in the reservoir for restoring. If zero, the sponge is disabled.", &
293 units=
"s", default=0.)
294 call get_param(param_file,
mdl,
"DUMBBELL_SREF", s_ref, do_not_log=.true.)
295 call get_param(param_file,
mdl,
"DUMBBELL_S_RANGE", s_range, do_not_log=.true.)
296 call get_param(param_file,
mdl,
"MIN_THICKNESS", min_thickness, &
297 'Minimum thickness for layer',&
298 units=
'm', default=1.0e-3, do_not_log=.true., scale=us%m_to_Z)
301 if (sponge_time_scale <= 0.)
return
308 if (g%mask2dT(i,j) > 0.)
then
310 x = (g%geoLonT(i,j) ) / dblen
311 if (x > 0.25 .or. x < -0.25)
then
313 idamp(i,j) = 1. / sponge_time_scale
321 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
322 eta1d(nz+1) = -g%bathyT(i,j)
324 eta1d(k) = -g%max_depth * real(k-1) / real(nz)
325 if (eta1d(k) < (eta1d(k+1) + min_thickness))
then
326 eta1d(k) = eta1d(k+1) + min_thickness
327 h(i,j,k) = gv%Z_to_H * min_thickness
329 h(i,j,k) = gv%Z_to_H * (eta1d(k) - eta1d(k+1))
340 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
342 x = ( g%geoLonT(i,j) ) / dblen
345 s(i,j,k)=s_ref + 0.5*s_range
350 s(i,j,k)=s_ref - 0.5*s_range