22 implicit none ;
private
24 #include <MOM_memory.h>
26 character(len=40) ::
mdl =
"seamount_initialization"
43 real,
dimension(G%isd:G%ied,G%jsd:G%jed), &
46 real,
intent(in) :: max_depth
50 real :: x, y, delta, lx, rlx, ly, rly
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.)
66 rlx = 0. ;
if (lx>0.) rlx = 1. / lx
67 rly = 0. ;
if (ly>0.) rly = 1. / ly
69 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
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) )
84 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
88 logical,
optional,
intent(in) :: just_read_params
93 real :: eta1d(szk_(g)+1)
96 real :: s_surf, s_range, s_ref, s_light, s_dense
98 character(len=20) :: verticalcoordinate
100 integer :: i, j, k, is, ie, js, je, nz
102 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
104 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
106 if (.not.just_read) &
107 call mom_mesg(
"MOM_initialization.F90, initialize_thickness_uniform: setting thickness")
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
209 real :: t_ref, t_light, t_dense, s_ref, s_light, s_dense, a1, frac_dense, k_frac, res_rat
211 character(len=20) :: verticalcoordinate, density_profile
213 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
215 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
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)
223 'Initial surface salinity', units=
'1e-3', default=34., do_not_log=just_read)
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)
233 select case ( coordinatemode(verticalcoordinate) )
234 case ( regridding_layer )
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
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
250 a1 = 2.0 * res_rat / (1.0 + res_rat)
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
260 if (just_read)
return
261 do j=js,je ;
do i=is,ie
264 xi1 = xi0 + gv%H_to_Z * h(i,j,k) / g%max_depth
265 select case ( trim(density_profile) )
268 s(i,j,k) = s_surf + ( 0.5 * s_range ) * (xi0 + xi1)
269 t(i,j,k) = t_surf + t_range * 0.5 * (xi0 + xi1)
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)
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)
278 call mom_error(fatal,
'Unknown value for "INITIAL_DENSITY_PROFILE"')