18 implicit none ;
private
20 #include <MOM_memory.h>
33 #include "version_variable.h"
42 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
46 logical,
optional,
intent(in) :: just_read_params
49 real :: eta0(szk_(g)+1)
50 real :: eta_im(szj_(g),szk_(g)+1)
51 real :: eta1d(szk_(g)+1)
58 character(len=40) :: mdl =
"Phillips_initialize_thickness"
59 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz
61 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
62 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
66 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
68 if (.not.just_read)
call log_version(param_file, mdl, version)
69 call get_param(param_file, mdl,
"HALF_STRAT_DEPTH", half_strat, &
70 "The fractional depth where the stratification is centered.", &
71 units=
"nondim", default = 0.5, do_not_log=just_read)
72 call get_param(param_file, mdl,
"JET_WIDTH", jet_width, &
73 "The width of the zonal-mean jet.", units=
"km", &
74 fail_if_missing=.not.just_read, do_not_log=just_read)
75 call get_param(param_file, mdl,
"JET_HEIGHT", jet_height, &
76 "The interface height scale associated with the "//&
77 "zonal-mean jet.", units=
"m", scale=us%m_to_Z, &
78 fail_if_missing=.not.just_read, do_not_log=just_read)
82 half_depth = g%max_depth*half_strat
83 eta0(1) = 0.0 ; eta0(nz+1) = -g%max_depth
84 do k=2,1+nz/2 ; eta0(k) = -half_depth*(2.0*(k-1)/real(nz)) ;
enddo
86 eta0(k) = -g%max_depth - 2.0*(g%max_depth-half_depth) * ((k-(nz+1))/real(nz))
90 eta_im(j,1) = 0.0 ; eta_im(j,nz+1) = -g%max_depth
92 do k=2,nz ;
do j=js,je
93 y_2 = g%geoLatT(is,j) - g%south_lat - 0.5*g%len_lat
94 eta_im(j,k) = eta0(k) + jet_height * tanh(y_2 / jet_width)
96 if (eta_im(j,k) > 0.0) eta_im(j,k) = 0.0
97 if (eta_im(j,k) < -g%max_depth) eta_im(j,k) = -g%max_depth
100 do j=js,je ;
do i=is,ie
105 eta1d(nz+1) = -g%bathyT(i,j)
107 eta1d(k) = eta_im(j,k)
108 if (eta1d(k) < (eta1d(k+1) + gv%Angstrom_Z))
then
109 eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
110 h(i,j,k) = gv%Angstrom_H
112 h(i,j,k) = gv%Z_to_H * (eta1d(k) - eta1d(k+1))
123 real,
dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
125 real,
dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
130 logical,
optional,
intent(in) :: just_read_params
137 real :: velocity_amplitude
139 integer :: i, j, k, is, ie, js, je, nz, m
141 character(len=40) :: mdl =
"Phillips_initialize_velocity"
142 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
144 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
146 if (.not.just_read)
call log_version(param_file, mdl, version)
147 call get_param(param_file, mdl,
"VELOCITY_IC_PERTURB_AMP", velocity_amplitude, &
148 "The magnitude of the initial velocity perturbation.", &
149 units=
"m s-1", default=0.001, scale=us%m_s_to_L_T, do_not_log=just_read)
150 call get_param(param_file, mdl,
"JET_WIDTH", jet_width, &
151 "The width of the zonal-mean jet.", units=
"km", &
152 fail_if_missing=.not.just_read, do_not_log=just_read)
153 call get_param(param_file, mdl,
"JET_HEIGHT", jet_height, &
154 "The interface height scale associated with the "//&
155 "zonal-mean jet.", units=
"m", scale=us%m_to_Z, &
156 fail_if_missing=.not.just_read, do_not_log=just_read)
158 if (just_read)
return
166 do k=nz-1,1 ;
do j=js,je ;
do i=is-1,ie
167 y_2 = g%geoLatCu(i,j) - g%south_lat - 0.5*g%len_lat
173 u(i,j,k) = u(i,j,k+1) + (1e-3 * (jet_height / (us%m_to_L*jet_width)) * &
174 (
sech(y_2 / jet_width))**2 ) * &
175 (2.0 * gv%g_prime(k+1) / (g%CoriolisBu(i,j) + g%CoriolisBu(i,j-1)))
176 enddo ;
enddo ;
enddo
178 do k=1,nz ;
do j=js,je ;
do i=is-1,ie
179 y_2 = (g%geoLatCu(i,j) - g%south_lat - 0.5*g%len_lat) / g%len_lat
180 x_2 = (g%geoLonCu(i,j) - g%west_lon - 0.5*g%len_lon) / g%len_lon
181 if (g%geoLonCu(i,j) == g%west_lon)
then
185 x_2 = ((g%west_lon + g%len_lon*real(g%ieg-(g%isg-1))/real(g%Domain%niglobal)) - &
186 g%west_lon - 0.5*g%len_lon) / g%len_lon
188 u(i,j,k) = u(i,j,k) + velocity_amplitude * ((real(k)-0.5)/real(nz)) * &
189 (0.5 - abs(2.0*x_2) + 0.1*abs(cos(10.0*pi*x_2)) - abs(sin(5.0*pi*y_2)))
191 u(i,j,k) = u(i,j,k) + 0.2*velocity_amplitude * ((real(k)-0.5)/real(nz)) * &
192 cos(2.0*m*pi*x_2 + 2*m) * cos(6.0*pi*y_2)
194 enddo ;
enddo ;
enddo
216 real,
intent(in),
dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h
219 real :: eta0(szk_(g)+1)
220 real :: eta(szi_(g),szj_(g),szk_(g)+1)
221 real :: temp(szi_(g),szj_(g),szk_(g))
222 real :: idamp(szi_(g),szj_(g))
223 real :: eta_im(szj_(g),szk_(g)+1)
224 real :: idamp_im(szj_(g))
231 character(len=40) :: mdl =
"Phillips_initialize_sponges"
233 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz
234 logical,
save :: first_call = .true.
236 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
237 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
239 eta(:,:,:) = 0.0 ; temp(:,:,:) = 0.0 ; idamp(:,:) = 0.0
240 eta_im(:,:) = 0.0 ; idamp_im(:) = 0.0
242 if (first_call)
call log_version(param_file, mdl, version)
244 call get_param(param_file, mdl,
"HALF_STRAT_DEPTH", half_strat, &
245 "The fractional depth where the stratificaiton is centered.", &
246 units=
"nondim", default = 0.5)
247 call get_param(param_file, mdl,
"SPONGE_RATE", damp_rate, &
248 "The rate at which the zonal-mean sponges damp.", units=
"s-1", &
249 default = 1.0/(10.0*86400.0))
251 call get_param(param_file, mdl,
"JET_WIDTH", jet_width, &
252 "The width of the zonal-mean jet.", units=
"km", &
253 fail_if_missing=.true.)
254 call get_param(param_file, mdl,
"JET_HEIGHT", jet_height, &
255 "The interface height scale associated with the "//&
256 "zonal-mean jet.", units=
"m", scale=us%m_to_Z, &
257 fail_if_missing=.true.)
259 half_depth = g%max_depth*half_strat
260 eta0(1) = 0.0 ; eta0(nz+1) = -g%max_depth
261 do k=2,1+nz/2 ; eta0(k) = -half_depth*(2.0*(k-1)/real(nz)) ;
enddo
263 eta0(k) = -g%max_depth - 2.0*(g%max_depth-half_depth) * ((k-(nz+1))/real(nz))
267 idamp_im(j) = damp_rate
268 eta_im(j,1) = 0.0 ; eta_im(j,nz+1) = -g%max_depth
270 do k=2,nz ;
do j=js,je
271 y_2 = g%geoLatT(is,j) - g%south_lat - 0.5*g%len_lat
272 eta_im(j,k) = eta0(k) + jet_height * tanh(y_2 / jet_width)
274 if (eta_im(j,k) > 0.0) eta_im(j,k) = 0.0
275 if (eta_im(j,k) < -g%max_depth) eta_im(j,k) = -g%max_depth
284 real,
intent(in) :: x
288 if (abs(x) > 228.)
then
291 sech = 2.0 / (exp(x) + exp(-x))
298 real,
dimension(G%isd:G%ied,G%jsd:G%jed), &
301 real,
intent(in) :: max_depth
306 real :: pi, htop, wtop, ltop, offset, dist
307 real :: x1, x2, x3, x4, y1, y2
308 integer :: i,j,is,ie,js,je
309 character(len=40) :: mdl =
"Phillips_initialize_topography"
311 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
314 m_to_z = 1.0 ;
if (
present(us)) m_to_z = us%m_to_Z
316 call get_param(param_file, mdl,
"PHILLIPS_HTOP", htop, &
317 "The maximum height of the topography.", units=
"m", scale=m_to_z, &
318 fail_if_missing=.true.)
321 ltop = 0.25*g%len_lon
322 offset = 0.1*g%len_lat
323 dist = 0.333*g%len_lon
326 y1=g%south_lat+0.5*g%len_lat+offset-0.5*wtop; y2=y1+wtop
327 x1=g%west_lon+0.1*g%len_lon; x2=x1+ltop; x3=x1+dist; x4=x3+3.0/2.0*ltop
329 do i=is,ie ;
do j=js,je
331 if (g%geoLonT(i,j)>x1 .and. g%geoLonT(i,j)<x2)
then
332 d(i,j) = htop*sin(pi*(g%geoLonT(i,j)-x1)/(x2-x1))**2
333 if (g%geoLatT(i,j)>y1 .and. g%geoLatT(i,j)<y2)
then
334 d(i,j) = d(i,j)*(1-sin(pi*(g%geoLatT(i,j)-y1)/(y2-y1))**2)
336 elseif (g%geoLonT(i,j)>x3 .and. g%geoLonT(i,j)<x4 .and. &
337 g%geoLatT(i,j)>y1 .and. g%geoLatT(i,j)<y2)
then
338 d(i,j) = 2.0/3.0*htop*sin(pi*(g%geoLonT(i,j)-x3)/(x4-x3))**2 &
339 *sin(pi*(g%geoLatT(i,j)-y1)/(y2-y1))**2
341 d(i,j) = max_depth - d(i,j)