18 implicit none ;
private
20 character(len=40) ::
mdl =
"adjustment_initialization"
22 #include <MOM_memory.h>
39 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
43 logical,
optional,
intent(in) :: just_read_params
48 real :: eta1d(szk_(g)+1)
53 real :: delta_s_strat, dsdz, delta_s, s_ref
54 real :: min_thickness, adjustment_width, adjustment_delta
55 real :: adjustment_deltas
56 real :: front_wave_amp, front_wave_length, front_wave_asym
57 real :: target_values(szk_(g)+1)
59 character(len=20) :: verticalcoordinate
61 #include "version_variable.h"
62 integer :: i, j, k, is, ie, js, je, nz
64 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
66 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
69 call mom_mesg(
"initialize_thickness_uniform: setting thickness")
73 call get_param(param_file,
mdl,
"S_REF",s_ref,fail_if_missing=.true.,do_not_log=.true.)
74 call get_param(param_file,
mdl,
"MIN_THICKNESS",min_thickness,
'Minimum layer thickness', &
75 units=
'm', default=1.0e-3, do_not_log=just_read, scale=us%m_to_Z)
78 call get_param(param_file,
mdl,
"REGRIDDING_COORDINATE_MODE",verticalcoordinate, &
79 default=default_coordinate_mode, do_not_log=just_read)
80 call get_param(param_file,
mdl,
"ADJUSTMENT_WIDTH",adjustment_width, &
81 "Width of frontal zone", &
82 units=
"same as x,y", fail_if_missing=.not.just_read, do_not_log=just_read)
83 call get_param(param_file,
mdl,
"DELTA_S_STRAT",delta_s_strat, &
84 "Top-to-bottom salinity difference of stratification", &
85 units=
"1e-3", fail_if_missing=.not.just_read, do_not_log=just_read)
86 call get_param(param_file,
mdl,
"ADJUSTMENT_DELTAS",adjustment_deltas, &
87 "Salinity difference across front", &
88 units=
"1e-3", fail_if_missing=.not.just_read, do_not_log=just_read)
89 call get_param(param_file,
mdl,
"FRONT_WAVE_AMP",front_wave_amp, &
90 "Amplitude of trans-frontal wave perturbation", &
91 units=
"same as x,y",default=0., do_not_log=just_read)
92 call get_param(param_file,
mdl,
"FRONT_WAVE_LENGTH",front_wave_length, &
93 "Wave-length of trans-frontal wave perturbation", &
94 units=
"same as x,y",default=0., do_not_log=just_read)
95 call get_param(param_file,
mdl,
"FRONT_WAVE_ASYM",front_wave_asym, &
96 "Amplitude of frontal asymmetric perturbation", &
97 default=0., do_not_log=just_read)
109 dsdz = -delta_s_strat / g%max_depth
111 select case ( coordinatemode(verticalcoordinate) )
114 drho_ds = 1.0 * us%kg_m3_to_R
115 if (delta_s_strat /= 0.)
then
117 adjustment_delta = (adjustment_deltas / delta_s_strat) * g%max_depth
119 e0(k) = adjustment_delta - (g%max_depth + 2*adjustment_delta) * (real(k-1) / real(nz))
122 adjustment_delta = 2.*g%max_depth
124 e0(k) = -g%max_depth * (real(k-1) / real(nz))
127 target_values(1) = ( gv%Rlay(1) + 0.5*(gv%Rlay(1)-gv%Rlay(2)) )
128 target_values(nz+1) = ( gv%Rlay(nz) + 0.5*(gv%Rlay(nz)-gv%Rlay(nz-1)) )
130 target_values(k) = target_values(k-1) + ( gv%Rlay(nz) - gv%Rlay(1) ) / (nz-1)
132 target_values(:) = target_values(:) - 1000.*us%kg_m3_to_R
133 do j=js,je ;
do i=is,ie
134 if (front_wave_length /= 0.)
then
135 y = ( 0.125 + g%geoLatT(i,j) / front_wave_length ) * ( 4. * acos(0.) )
136 yy = 2. * ( g%geoLatT(i,j) - 0.5 * g%len_lat ) / adjustment_width
137 yy = min(1.0, yy); yy = max(-1.0, yy)
138 yy = yy * 2. * acos( 0. )
139 y = front_wave_amp*sin(y) + front_wave_asym*sin(yy)
143 x = ( ( g%geoLonT(i,j) - 0.5 * g%len_lon ) + y ) / adjustment_width
144 x = min(1.0, x); x = max(-1.0, x)
146 delta_s = adjustment_deltas * 0.5 * (1. - sin( x ) )
148 if (drho_ds*dsdz /= 0.)
then
149 eta1d(k) = ( target_values(k) - drho_ds*( s_ref + delta_s ) ) / (drho_ds*dsdz)
151 eta1d(k) = e0(k) - (0.5*adjustment_delta) * sin( x )
153 eta1d(k) = max( eta1d(k), -g%max_depth )
154 eta1d(k) = min( eta1d(k), 0. )
156 eta1d(1) = 0.; eta1d(nz+1) = -g%max_depth
158 if (eta1d(k) > 0.)
then
159 eta1d(k) = max( eta1d(k+1) + min_thickness, 0. )
160 h(i,j,k) = gv%Z_to_H * max( eta1d(k) - eta1d(k+1), min_thickness )
161 elseif (eta1d(k) <= (eta1d(k+1) + min_thickness))
then
162 eta1d(k) = eta1d(k+1) + min_thickness
163 h(i,j,k) = gv%Z_to_H * min_thickness
165 h(i,j,k) = gv%Z_to_H * (eta1d(k) - eta1d(k+1))
172 eta1d(k) = -g%max_depth * (real(k-1) / real(nz))
173 eta1d(k) = max(min(eta1d(k), 0.), -g%max_depth)
175 do j=js,je ;
do i=is,ie
177 h(i,j,k) = gv%Z_to_H * (eta1d(k) - eta1d(k+1))
182 call mom_error(fatal,
"adjustment_initialize_thickness: "// &
183 "Unrecognized i.c. setup - set ADJUSTMENT_IC")
191 eqn_of_state, just_read_params)
194 real,
dimension(SZI_(G),SZJ_(G), SZK_(G)),
intent(out) :: t
195 real,
dimension(SZI_(G),SZJ_(G), SZK_(G)),
intent(out) :: s
196 real,
dimension(SZI_(G),SZJ_(G), SZK_(G)),
intent(in) :: h
199 type(
eos_type),
pointer :: eqn_of_state
200 logical,
optional,
intent(in) :: just_read_params
203 integer :: i, j, k, is, ie, js, je, nz
205 integer :: index_bay_z
208 real :: s_range, t_range
210 real :: xi0, xi1, dsdz, delta_s, delta_s_strat
211 real :: adjustment_width, adjustment_deltas
212 real :: front_wave_amp, front_wave_length, front_wave_asym
213 real :: eta1d(szk_(g)+1)
215 character(len=20) :: verticalcoordinate
217 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
219 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
222 call get_param(param_file,
mdl,
"S_REF",s_ref,
'Reference salinity', units=
'1e-3', &
223 fail_if_missing=.not.just_read, do_not_log=just_read)
224 call get_param(param_file,
mdl,
"T_REF",t_ref,
'Reference temperature', units=
'C', &
225 fail_if_missing=.not.just_read, do_not_log=just_read)
226 call get_param(param_file,
mdl,
"S_RANGE",s_range,
'Initial salinity range', units=
'1e-3', &
227 default=2.0, do_not_log=just_read)
228 call get_param(param_file,
mdl,
"T_RANGE",t_range,
'Initial temperature range', units=
'C', &
229 default=0.0, do_not_log=just_read)
231 call get_param(param_file,
mdl,
"REGRIDDING_COORDINATE_MODE",verticalcoordinate, &
232 default=default_coordinate_mode, do_not_log=just_read)
233 call get_param(param_file,
mdl,
"ADJUSTMENT_WIDTH", adjustment_width, &
234 fail_if_missing=.not.just_read, do_not_log=.true.)
235 call get_param(param_file,
mdl,
"ADJUSTMENT_DELTAS", adjustment_deltas, &
236 fail_if_missing=.not.just_read, do_not_log=.true.)
237 call get_param(param_file,
mdl,
"DELTA_S_STRAT", delta_s_strat, &
238 fail_if_missing=.not.just_read, do_not_log=.true.)
239 call get_param(param_file,
mdl,
"FRONT_WAVE_AMP", front_wave_amp, default=0., &
241 call get_param(param_file,
mdl,
"FRONT_WAVE_LENGTH",front_wave_length, &
242 default=0., do_not_log=.true.)
243 call get_param(param_file,
mdl,
"FRONT_WAVE_ASYM", front_wave_asym, default=0., &
246 if (just_read)
return
252 select case ( coordinatemode(verticalcoordinate) )
255 dsdz = -delta_s_strat / g%max_depth
256 do j=js,je ;
do i=is,ie
257 eta1d(nz+1) = -g%bathyT(i,j)
259 eta1d(k) = eta1d(k+1) + h(i,j,k)*gv%H_to_Z
261 if (front_wave_length /= 0.)
then
262 y = ( 0.125 + g%geoLatT(i,j) / front_wave_length ) * ( 4. * acos(0.) )
263 yy = 2. * ( g%geoLatT(i,j) - 0.5 * g%len_lat ) / front_wave_length
264 yy = min(1.0, yy); yy = max(-1.0, yy)
265 yy = yy * 2. * acos( 0. )
266 y = front_wave_amp*sin(y) + front_wave_asym*sin(yy)
270 x = ( ( g%geoLonT(i,j) - 0.5 * g%len_lon ) + y ) / adjustment_width
271 x = min(1.0, x); x = max(-1.0, x)
273 delta_s = adjustment_deltas * 0.5 * (1. - sin( x ) )
275 s(i,j,k) = s_ref + delta_s + 0.5 * ( eta1d(k)+eta1d(k+1) ) * dsdz
276 x = abs(s(i,j,k) - 0.5*real(nz-1)/real(nz)*s_range)/s_range*real(2*nz)
286 s(:,:,k) = s_ref + s_range * ( (real(k)-0.5) / real( nz ) )
293 call mom_error(fatal,
"adjustment_initialize_temperature_salinity: "// &
294 "Unrecognized i.c. setup - set ADJUSTMENT_IC")