14 use regrid_consts,
only : coordinatemode, default_coordinate_mode
18 implicit none ;
private
20 character(len=40) :: mdl =
"adjustment_initialization"
22 #include <MOM_memory.h>
24 public adjustment_initialize_thickness
25 public adjustment_initialize_temperature_salinity
35 subroutine adjustment_initialize_thickness ( h, G, GV, US, param_file, just_read_params)
39 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
43 logical,
optional,
intent(in) :: just_read_params
48 real :: eta1d(szk_(g)+1)
50 real :: x, y, yy, delta_s_strat, dsdz, delta_s, s_ref
51 real :: min_thickness, adjustment_width, adjustment_delta, adjustment_deltas
52 real :: front_wave_amp, front_wave_length, front_wave_asym
53 real :: target_values(szk_(g)+1)
55 character(len=20) :: verticalcoordinate
57 #include "version_variable.h"
58 integer :: i, j, k, is, ie, js, je, nz
60 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
62 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
65 call mom_mesg(
"initialize_thickness_uniform: setting thickness")
68 if (.not.just_read)
call log_version(param_file, mdl, version,
"")
69 call get_param(param_file, mdl,
"S_REF",s_ref,fail_if_missing=.true.,do_not_log=.true.)
70 call get_param(param_file, mdl,
"MIN_THICKNESS",min_thickness,
'Minimum layer thickness', &
71 units=
'm', default=1.0e-3, do_not_log=just_read, scale=us%m_to_Z)
74 call get_param(param_file, mdl,
"REGRIDDING_COORDINATE_MODE",verticalcoordinate, &
75 default=default_coordinate_mode, do_not_log=just_read)
76 call get_param(param_file, mdl,
"ADJUSTMENT_WIDTH",adjustment_width, &
77 "Width of frontal zone", &
78 units=
"same as x,y", fail_if_missing=.not.just_read, do_not_log=just_read)
79 call get_param(param_file, mdl,
"DELTA_S_STRAT",delta_s_strat, &
80 "Top-to-bottom salinity difference of stratification", &
81 units=
"1e-3", fail_if_missing=.not.just_read, do_not_log=just_read)
82 call get_param(param_file, mdl,
"ADJUSTMENT_DELTAS",adjustment_deltas, &
83 "Salinity difference across front", &
84 units=
"1e-3", fail_if_missing=.not.just_read, do_not_log=just_read)
85 call get_param(param_file, mdl,
"FRONT_WAVE_AMP",front_wave_amp, &
86 "Amplitude of trans-frontal wave perturbation", &
87 units=
"same as x,y",default=0., do_not_log=just_read)
88 call get_param(param_file, mdl,
"FRONT_WAVE_LENGTH",front_wave_length, &
89 "Wave-length of trans-frontal wave perturbation", &
90 units=
"same as x,y",default=0., do_not_log=just_read)
91 call get_param(param_file, mdl,
"FRONT_WAVE_ASYM",front_wave_asym, &
92 "Amplitude of frontal asymmetric perturbation", &
93 default=0., do_not_log=just_read)
105 dsdz = -delta_s_strat / g%max_depth
107 select case ( coordinatemode(verticalcoordinate) )
109 case ( regridding_layer, regridding_rho )
110 if (delta_s_strat /= 0.)
then
112 adjustment_delta = (adjustment_deltas / delta_s_strat) * g%max_depth
114 e0(k) = adjustment_delta - (g%max_depth + 2*adjustment_delta) * (real(k-1) / real(nz))
117 adjustment_delta = 2.*g%max_depth
119 e0(k) = -g%max_depth * (real(k-1) / real(nz))
122 target_values(1) = gv%Rlay(1) + 0.5*(gv%Rlay(1)-gv%Rlay(2))
123 target_values(nz+1) = gv%Rlay(nz) + 0.5*(gv%Rlay(nz)-gv%Rlay(nz-1))
125 target_values(k) = target_values(k-1) + ( gv%Rlay(nz) - gv%Rlay(1) ) / (nz-1)
127 target_values(:) = target_values(:) - 1000.
128 do j=js,je ;
do i=is,ie
129 if (front_wave_length /= 0.)
then
130 y = ( 0.125 + g%geoLatT(i,j) / front_wave_length ) * ( 4. * acos(0.) )
131 yy = 2. * ( g%geoLatT(i,j) - 0.5 * g%len_lat ) / adjustment_width
132 yy = min(1.0, yy); yy = max(-1.0, yy)
133 yy = yy * 2. * acos( 0. )
134 y = front_wave_amp*sin(y) + front_wave_asym*sin(yy)
138 x = ( ( g%geoLonT(i,j) - 0.5 * g%len_lon ) + y ) / adjustment_width
139 x = min(1.0, x); x = max(-1.0, x)
141 delta_s = adjustment_deltas * 0.5 * (1. - sin( x ) )
144 eta1d(k) = ( target_values(k) - ( s_ref + delta_s ) ) / dsdz
146 eta1d(k) = e0(k) - (0.5*adjustment_delta) * sin( x )
148 eta1d(k) = max( eta1d(k), -g%max_depth )
149 eta1d(k) = min( eta1d(k), 0. )
151 eta1d(1) = 0.; eta1d(nz+1) = -g%max_depth
153 if (eta1d(k) > 0.)
then
154 eta1d(k) = max( eta1d(k+1) + min_thickness, 0. )
155 h(i,j,k) = gv%Z_to_H * max( eta1d(k) - eta1d(k+1), min_thickness )
156 elseif (eta1d(k) <= (eta1d(k+1) + min_thickness))
then
157 eta1d(k) = eta1d(k+1) + min_thickness
158 h(i,j,k) = gv%Z_to_H * min_thickness
160 h(i,j,k) = gv%Z_to_H * (eta1d(k) - eta1d(k+1))
165 case ( regridding_zstar, regridding_sigma )
167 eta1d(k) = -g%max_depth * (real(k-1) / real(nz))
168 eta1d(k) = max(min(eta1d(k), 0.), -g%max_depth)
170 do j=js,je ;
do i=is,ie
172 h(i,j,k) = gv%Z_to_H * (eta1d(k) - eta1d(k+1))
177 call mom_error(fatal,
"adjustment_initialize_thickness: "// &
178 "Unrecognized i.c. setup - set ADJUSTMENT_IC")
182 end subroutine adjustment_initialize_thickness
185 subroutine adjustment_initialize_temperature_salinity(T, S, h, G, GV, param_file, &
186 eqn_of_state, just_read_params)
189 real,
dimension(SZI_(G),SZJ_(G), SZK_(G)),
intent(out) :: t
190 real,
dimension(SZI_(G),SZJ_(G), SZK_(G)),
intent(out) :: s
191 real,
dimension(SZI_(G),SZJ_(G), SZK_(G)),
intent(in) :: h
194 type(
eos_type),
pointer :: eqn_of_state
195 logical,
optional,
intent(in) :: just_read_params
198 integer :: i, j, k, is, ie, js, je, nz
200 integer :: index_bay_z
203 real :: s_range, t_range
205 real :: xi0, xi1, dsdz, delta_s, delta_s_strat
206 real :: adjustment_width, adjustment_deltas
207 real :: front_wave_amp, front_wave_length, front_wave_asym
208 real :: eta1d(szk_(g)+1)
210 character(len=20) :: verticalcoordinate
212 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
214 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
217 call get_param(param_file, mdl,
"S_REF",s_ref,
'Reference salinity', units=
'1e-3', &
218 fail_if_missing=.not.just_read, do_not_log=just_read)
219 call get_param(param_file, mdl,
"T_REF",t_ref,
'Reference temperature', units=
'C', &
220 fail_if_missing=.not.just_read, do_not_log=just_read)
221 call get_param(param_file, mdl,
"S_RANGE",s_range,
'Initial salinity range', units=
'1e-3', &
222 default=2.0, do_not_log=just_read)
223 call get_param(param_file, mdl,
"T_RANGE",t_range,
'Initial temperature range', units=
'C', &
224 default=0.0, do_not_log=just_read)
226 call get_param(param_file, mdl,
"REGRIDDING_COORDINATE_MODE",verticalcoordinate, &
227 default=default_coordinate_mode, do_not_log=just_read)
228 call get_param(param_file, mdl,
"ADJUSTMENT_WIDTH", adjustment_width, &
229 fail_if_missing=.not.just_read, do_not_log=.true.)
230 call get_param(param_file, mdl,
"ADJUSTMENT_DELTAS", adjustment_deltas, &
231 fail_if_missing=.not.just_read, do_not_log=.true.)
232 call get_param(param_file, mdl,
"DELTA_S_STRAT", delta_s_strat, &
233 fail_if_missing=.not.just_read, do_not_log=.true.)
234 call get_param(param_file, mdl,
"FRONT_WAVE_AMP", front_wave_amp, default=0., &
236 call get_param(param_file, mdl,
"FRONT_WAVE_LENGTH",front_wave_length, &
237 default=0., do_not_log=.true.)
238 call get_param(param_file, mdl,
"FRONT_WAVE_ASYM", front_wave_asym, default=0., &
241 if (just_read)
return
247 select case ( coordinatemode(verticalcoordinate) )
249 case ( regridding_zstar, regridding_sigma )
250 dsdz = -delta_s_strat / g%max_depth
251 do j=js,je ;
do i=is,ie
252 eta1d(nz+1) = -g%bathyT(i,j)
254 eta1d(k) = eta1d(k+1) + h(i,j,k)*gv%H_to_Z
256 if (front_wave_length /= 0.)
then
257 y = ( 0.125 + g%geoLatT(i,j) / front_wave_length ) * ( 4. * acos(0.) )
258 yy = 2. * ( g%geoLatT(i,j) - 0.5 * g%len_lat ) / front_wave_length
259 yy = min(1.0, yy); yy = max(-1.0, yy)
260 yy = yy * 2. * acos( 0. )
261 y = front_wave_amp*sin(y) + front_wave_asym*sin(yy)
265 x = ( ( g%geoLonT(i,j) - 0.5 * g%len_lon ) + y ) / adjustment_width
266 x = min(1.0, x); x = max(-1.0, x)
268 delta_s = adjustment_deltas * 0.5 * (1. - sin( x ) )
270 s(i,j,k) = s_ref + delta_s + 0.5 * ( eta1d(k)+eta1d(k+1) ) * dsdz
271 x = abs(s(i,j,k) - 0.5*real(nz-1)/real(nz)*s_range)/s_range*real(2*nz)
279 case ( regridding_layer, regridding_rho )
281 s(:,:,k) = s_ref + s_range * ( (real(k)-0.5) / real( nz ) )
288 call mom_error(fatal,
"adjustment_initialize_temperature_salinity: "// &
289 "Unrecognized i.c. setup - set ADJUSTMENT_IC")
293 end subroutine adjustment_initialize_temperature_salinity