18 implicit none ;
private
20 #include <MOM_memory.h>
22 public phillips_initialize_thickness
23 public phillips_initialize_velocity
24 public phillips_initialize_sponges
25 public phillips_initialize_topography
33 #include "version_variable.h"
38 subroutine phillips_initialize_thickness(h, G, GV, US, param_file, just_read_params)
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))
117 end subroutine phillips_initialize_thickness
120 subroutine phillips_initialize_velocity(u, v, G, GV, US, param_file, just_read_params)
124 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
126 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
130 logical,
optional,
intent(in) :: just_read_params
133 real :: jet_width, jet_height, x_2, y_2
134 real :: velocity_amplitude, pi
135 integer :: i, j, k, is, ie, js, je, nz, m
137 character(len=40) :: mdl =
"Phillips_initialize_velocity"
138 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
140 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
142 if (.not.just_read)
call log_version(param_file, mdl, version)
143 call get_param(param_file, mdl,
"VELOCITY_IC_PERTURB_AMP", velocity_amplitude, &
144 "The magnitude of the initial velocity perturbation.", &
145 units=
"m s-1", default=0.001, do_not_log=just_read)
146 call get_param(param_file, mdl,
"JET_WIDTH", jet_width, &
147 "The width of the zonal-mean jet.", units=
"km", &
148 fail_if_missing=.not.just_read, do_not_log=just_read)
149 call get_param(param_file, mdl,
"JET_HEIGHT", jet_height, &
150 "The interface height scale associated with the "//&
151 "zonal-mean jet.", units=
"m", scale=us%m_to_Z, &
152 fail_if_missing=.not.just_read, do_not_log=just_read)
154 if (just_read)
return
162 do k=nz-1,1 ;
do j=js,je ;
do i=is-1,ie
163 y_2 = g%geoLatCu(i,j) - g%south_lat - 0.5*g%len_lat
169 u(i,j,k) = u(i,j,k+1) + (1e-3 * (jet_height / jet_width) * &
170 (sech(y_2 / jet_width))**2 ) * &
171 (2.0 * us%L_to_m**2*us%s_to_T**2*gv%g_prime(k+1) * us%T_to_s / (g%CoriolisBu(i,j) + g%CoriolisBu(i,j-1)))
172 enddo ;
enddo ;
enddo
174 do k=1,nz ;
do j=js,je ;
do i=is-1,ie
175 y_2 = (g%geoLatCu(i,j) - g%south_lat - 0.5*g%len_lat) / g%len_lat
176 x_2 = (g%geoLonCu(i,j) - g%west_lon - 0.5*g%len_lon) / g%len_lon
177 if (g%geoLonCu(i,j) == g%west_lon)
then
181 x_2 = ((g%west_lon + g%len_lon*real(g%ieg-(g%isg-1))/real(g%Domain%niglobal)) - &
182 g%west_lon - 0.5*g%len_lon) / g%len_lon
184 u(i,j,k) = u(i,j,k) + velocity_amplitude * ((real(k)-0.5)/real(nz)) * &
185 (0.5 - abs(2.0*x_2) + 0.1*abs(cos(10.0*pi*x_2)) - abs(sin(5.0*pi*y_2)))
187 u(i,j,k) = u(i,j,k) + 0.2*velocity_amplitude * ((real(k)-0.5)/real(nz)) * &
188 cos(2.0*m*pi*x_2 + 2*m) * cos(6.0*pi*y_2)
190 enddo ;
enddo ;
enddo
192 end subroutine phillips_initialize_velocity
197 subroutine phillips_initialize_sponges(G, GV, US, tv, param_file, CSp, h)
212 real,
intent(in),
dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h
215 real :: eta0(szk_(g)+1)
216 real :: eta(szi_(g),szj_(g),szk_(g)+1)
217 real :: temp(szi_(g),szj_(g),szk_(g))
218 real :: idamp(szi_(g),szj_(g))
219 real :: eta_im(szj_(g),szk_(g)+1)
220 real :: idamp_im(szj_(g))
227 character(len=40) :: mdl =
"Phillips_initialize_sponges"
229 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz
230 logical,
save :: first_call = .true.
232 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
233 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
235 eta(:,:,:) = 0.0 ; temp(:,:,:) = 0.0 ; idamp(:,:) = 0.0
236 eta_im(:,:) = 0.0 ; idamp_im(:) = 0.0
238 if (first_call)
call log_version(param_file, mdl, version)
240 call get_param(param_file, mdl,
"HALF_STRAT_DEPTH", half_strat, &
241 "The fractional depth where the stratificaiton is centered.", &
242 units=
"nondim", default = 0.5)
243 call get_param(param_file, mdl,
"SPONGE_RATE", damp_rate, &
244 "The rate at which the zonal-mean sponges damp.", units=
"s-1", &
245 default = 1.0/(10.0*86400.0))
247 call get_param(param_file, mdl,
"JET_WIDTH", jet_width, &
248 "The width of the zonal-mean jet.", units=
"km", &
249 fail_if_missing=.true.)
250 call get_param(param_file, mdl,
"JET_HEIGHT", jet_height, &
251 "The interface height scale associated with the "//&
252 "zonal-mean jet.", units=
"m", scale=us%m_to_Z, &
253 fail_if_missing=.true.)
255 half_depth = g%max_depth*half_strat
256 eta0(1) = 0.0 ; eta0(nz+1) = -g%max_depth
257 do k=2,1+nz/2 ; eta0(k) = -half_depth*(2.0*(k-1)/real(nz)) ;
enddo
259 eta0(k) = -g%max_depth - 2.0*(g%max_depth-half_depth) * ((k-(nz+1))/real(nz))
263 idamp_im(j) = damp_rate
264 eta_im(j,1) = 0.0 ; eta_im(j,nz+1) = -g%max_depth
266 do k=2,nz ;
do j=js,je
267 y_2 = g%geoLatT(is,j) - g%south_lat - 0.5*g%len_lat
268 eta_im(j,k) = eta0(k) + jet_height * tanh(y_2 / jet_width)
270 if (eta_im(j,k) > 0.0) eta_im(j,k) = 0.0
271 if (eta_im(j,k) < -g%max_depth) eta_im(j,k) = -g%max_depth
274 call initialize_sponge(idamp, eta, g, param_file, csp, gv, idamp_im, eta_im)
276 end subroutine phillips_initialize_sponges
280 real,
intent(in) :: x
284 if (abs(x) > 228.)
then
287 sech = 2.0 / (exp(x) + exp(-x))
292 subroutine phillips_initialize_topography(D, G, param_file, max_depth, US)
294 real,
dimension(G%isd:G%ied,G%jsd:G%jed), &
297 real,
intent(in) :: max_depth
302 real :: pi, htop, wtop, ltop, offset, dist
303 real :: x1, x2, x3, x4, y1, y2
304 integer :: i,j,is,ie,js,je
305 character(len=40) :: mdl =
"Phillips_initialize_topography"
307 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
310 m_to_z = 1.0 ;
if (
present(us)) m_to_z = us%m_to_Z
312 call get_param(param_file, mdl,
"PHILLIPS_HTOP", htop, &
313 "The maximum height of the topography.", units=
"m", scale=m_to_z, &
314 fail_if_missing=.true.)
317 ltop = 0.25*g%len_lon
318 offset = 0.1*g%len_lat
319 dist = 0.333*g%len_lon
322 y1=g%south_lat+0.5*g%len_lat+offset-0.5*wtop; y2=y1+wtop
323 x1=g%west_lon+0.1*g%len_lon; x2=x1+ltop; x3=x1+dist; x4=x3+3.0/2.0*ltop
325 do i=is,ie ;
do j=js,je
327 if (g%geoLonT(i,j)>x1 .and. g%geoLonT(i,j)<x2)
then
328 d(i,j) = htop*sin(pi*(g%geoLonT(i,j)-x1)/(x2-x1))**2
329 if (g%geoLatT(i,j)>y1 .and. g%geoLatT(i,j)<y2)
then
330 d(i,j) = d(i,j)*(1-sin(pi*(g%geoLatT(i,j)-y1)/(y2-y1))**2)
332 elseif (g%geoLonT(i,j)>x3 .and. g%geoLonT(i,j)<x4 .and. &
333 g%geoLatT(i,j)>y1 .and. g%geoLatT(i,j)<y2)
then
334 d(i,j) = 2.0/3.0*htop*sin(pi*(g%geoLonT(i,j)-x3)/(x4-x3))**2 &
335 *sin(pi*(g%geoLatT(i,j)-y1)/(y2-y1))**2
337 d(i,j) = max_depth - d(i,j)
340 end subroutine phillips_initialize_topography