12 implicit none ;
private
14 #include <MOM_memory.h>
16 public calc_isoneutral_slopes, vert_fill_ts
26 subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, &
27 slope_x, slope_y, N2_u, N2_v, halo)
31 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
32 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1),
intent(in) :: e
36 real,
intent(in) :: dt_kappa_smooth
38 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)+1),
intent(inout) :: slope_x
39 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)+1),
intent(inout) :: slope_y
40 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)+1), &
41 optional,
intent(inout) :: n2_u
43 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)+1), &
44 optional,
intent(inout) :: n2_v
46 integer,
optional,
intent(in) :: halo
51 real,
dimension(SZI_(G), SZJ_(G), SZK_(G)) :: &
57 real,
dimension(SZI_(G), SZJ_(G), SZK_(G)+1) :: &
59 real,
dimension(SZIB_(G)) :: &
62 real,
dimension(SZI_(G)) :: &
65 real,
dimension(SZIB_(G)) :: &
69 real,
dimension(SZI_(G)) :: &
79 real :: haa, hab, hal, har
81 real :: wta, wtb, wtl, wtr
95 real :: g_rho0, n2, dzn2, h_x(szib_(g)), h_y(szi_(g))
102 logical :: present_n2_u, present_n2_v
103 integer :: is, ie, js, je, nz, isdb
106 if (
present(halo))
then
107 is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo
109 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
111 nz = g%ke ; isdb = g%IsdB
113 h_neglect = gv%H_subroundoff ; h_neglect2 = h_neglect**2
114 z_to_l = us%Z_to_m ; h_to_z = gv%H_to_Z
118 l_to_z = 1.0 / z_to_l
119 dz_neglect = gv%H_subroundoff * h_to_z
121 use_eos =
associated(tv%eqn_of_state)
123 present_n2_u =
PRESENT(n2_u)
124 present_n2_v =
PRESENT(n2_v)
125 g_rho0 = (us%L_to_Z*us%L_to_m*l_to_z*us%s_to_T**2*gv%g_Earth) / gv%Rho0
126 if (present_n2_u)
then
127 do j=js,je ;
do i=is-1,ie
132 if (present_n2_v)
then
133 do j=js-1,je ;
do i=is,ie
140 if (
present(halo))
then
141 call vert_fill_ts(h, tv%T, tv%S, dt_kappa_smooth, t, s, g, gv, halo+1)
143 call vert_fill_ts(h, tv%T, tv%S, dt_kappa_smooth, t, s, g, gv, 1)
149 do j=js-1,je+1 ;
do i=is-1,ie+1
151 pres(i,j,2) = pres(i,j,1) + gv%H_to_Pa*h(i,j,1)
155 do k=2,nz ;
do i=is-1,ie+1
156 pres(i,j,k+1) = pres(i,j,k) + gv%H_to_Pa*h(i,j,k)
167 do j=js,je ;
do k=nz,2,-1
168 if (.not.(use_eos))
then
169 drdia = 0.0 ; drdib = 0.0
170 drdkl = gv%Rlay(k)-gv%Rlay(k-1) ; drdkr = gv%Rlay(k)-gv%Rlay(k-1)
176 pres_u(i) = 0.5*(pres(i,j,k) + pres(i+1,j,k))
177 t_u(i) = 0.25*((t(i,j,k) + t(i+1,j,k)) + (t(i,j,k-1) + t(i+1,j,k-1)))
178 s_u(i) = 0.25*((s(i,j,k) + s(i+1,j,k)) + (s(i,j,k-1) + s(i+1,j,k-1)))
181 drho_ds_u, (is-isdb+1)-1, ie-is+2, tv%eqn_of_state)
187 drdia = drho_dt_u(i) * (t(i+1,j,k-1)-t(i,j,k-1)) + &
188 drho_ds_u(i) * (s(i+1,j,k-1)-s(i,j,k-1))
189 drdib = drho_dt_u(i) * (t(i+1,j,k)-t(i,j,k)) + &
190 drho_ds_u(i) * (s(i+1,j,k)-s(i,j,k))
193 drdkl = (drho_dt_u(i) * (t(i,j,k)-t(i,j,k-1)) + &
194 drho_ds_u(i) * (s(i,j,k)-s(i,j,k-1)))
195 drdkr = (drho_dt_u(i) * (t(i+1,j,k)-t(i+1,j,k-1)) + &
196 drho_ds_u(i) * (s(i+1,j,k)-s(i+1,j,k-1)))
201 hg2a = h(i,j,k-1)*h(i+1,j,k-1) + h_neglect2
202 hg2b = h(i,j,k)*h(i+1,j,k) + h_neglect2
203 hg2l = h(i,j,k-1)*h(i,j,k) + h_neglect2
204 hg2r = h(i+1,j,k-1)*h(i+1,j,k) + h_neglect2
205 haa = 0.5*(h(i,j,k-1) + h(i+1,j,k-1))
206 hab = 0.5*(h(i,j,k) + h(i+1,j,k)) + h_neglect
207 hal = 0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect
208 har = 0.5*(h(i+1,j,k-1) + h(i+1,j,k)) + h_neglect
209 if (gv%Boussinesq)
then
210 dzal = hal * h_to_z ; dzar = har * h_to_z
212 dzal = 0.5*(e(i,j,k-1) - e(i,j,k+1)) + dz_neglect
213 dzar = 0.5*(e(i+1,j,k-1) - e(i+1,j,k+1)) + dz_neglect
217 wta = hg2a*hab ; wtb = hg2b*haa
218 wtl = hg2l*(har*dzar) ; wtr = hg2r*(hal*dzal)
220 drdz = (wtl * drdkl + wtr * drdkr) / (dzal*wtl + dzar*wtr)
225 drdx = ((wta * drdia + wtb * drdib) / (wta + wtb) - &
226 drdz * (e(i,j,k)-e(i+1,j,k))) * g%IdxCu(i,j)
230 mag_grad2 = drdx**2 + (l_to_z*drdz)**2
231 if (mag_grad2 > 0.0)
then
232 slope_x(i,j,k) = drdx / sqrt(mag_grad2)
237 if (present_n2_u) n2_u(i,j,k) = g_rho0 * drdz * g%mask2dCu(i,j)
240 slope_x(i,j,k) = (z_to_l*(e(i,j,k)-e(i+1,j,k))) * g%IdxCu(i,j)
254 do j=js-1,je ;
do k=nz,2,-1
255 if (.not.(use_eos))
then
256 drdja = 0.0 ; drdjb = 0.0
257 drdkl = gv%Rlay(k)-gv%Rlay(k-1) ; drdkr = gv%Rlay(k)-gv%Rlay(k-1)
262 pres_v(i) = 0.5*(pres(i,j,k) + pres(i,j+1,k))
263 t_v(i) = 0.25*((t(i,j,k) + t(i,j+1,k)) + (t(i,j,k-1) + t(i,j+1,k-1)))
264 s_v(i) = 0.25*((s(i,j,k) + s(i,j+1,k)) + (s(i,j,k-1) + s(i,j+1,k-1)))
267 drho_ds_v, is, ie-is+1, tv%eqn_of_state)
272 drdja = drho_dt_v(i) * (t(i,j+1,k-1)-t(i,j,k-1)) + &
273 drho_ds_v(i) * (s(i,j+1,k-1)-s(i,j,k-1))
274 drdjb = drho_dt_v(i) * (t(i,j+1,k)-t(i,j,k)) + &
275 drho_ds_v(i) * (s(i,j+1,k)-s(i,j,k))
278 drdkl = (drho_dt_v(i) * (t(i,j,k)-t(i,j,k-1)) + &
279 drho_ds_v(i) * (s(i,j,k)-s(i,j,k-1)))
280 drdkr = (drho_dt_v(i) * (t(i,j+1,k)-t(i,j+1,k-1)) + &
281 drho_ds_v(i) * (s(i,j+1,k)-s(i,j+1,k-1)))
285 hg2a = h(i,j,k-1)*h(i,j+1,k-1) + h_neglect2
286 hg2b = h(i,j,k)*h(i,j+1,k) + h_neglect2
287 hg2l = h(i,j,k-1)*h(i,j,k) + h_neglect2
288 hg2r = h(i,j+1,k-1)*h(i,j+1,k) + h_neglect2
289 haa = 0.5*(h(i,j,k-1) + h(i,j+1,k-1)) + h_neglect
290 hab = 0.5*(h(i,j,k) + h(i,j+1,k)) + h_neglect
291 hal = 0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect
292 har = 0.5*(h(i,j+1,k-1) + h(i,j+1,k)) + h_neglect
293 if (gv%Boussinesq)
then
294 dzal = hal * h_to_z ; dzar = har * h_to_z
296 dzal = 0.5*(e(i,j,k-1) - e(i,j,k+1)) + dz_neglect
297 dzar = 0.5*(e(i,j+1,k-1) - e(i,j+1,k+1)) + dz_neglect
301 wta = hg2a*hab ; wtb = hg2b*haa
302 wtl = hg2l*(har*dzar) ; wtr = hg2r*(hal*dzal)
304 drdz = (wtl * drdkl + wtr * drdkr) / (dzal*wtl + dzar*wtr)
309 drdy = ((wta * drdja + wtb * drdjb) / (wta + wtb) - &
310 drdz * (e(i,j,k)-e(i,j+1,k))) * g%IdyCv(i,j)
314 mag_grad2 = drdy**2 + (l_to_z*drdz)**2
315 if (mag_grad2 > 0.0)
then
316 slope_y(i,j,k) = drdy / sqrt(mag_grad2)
321 if (present_n2_v) n2_v(i,j,k) = g_rho0 * drdz * g%mask2dCv(i,j)
324 slope_y(i,j,k) = (z_to_l*(e(i,j,k)-e(i,j+1,k))) * g%IdyCv(i,j)
330 end subroutine calc_isoneutral_slopes
334 subroutine vert_fill_ts(h, T_in, S_in, kappa_dt, T_f, S_f, G, GV, halo_here, larger_h_denom)
337 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
338 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: t_in
339 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: s_in
340 real,
intent(in) :: kappa_dt
342 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(out) :: t_f
343 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(out) :: s_f
344 integer,
optional,
intent(in) :: halo_here
346 logical,
optional,
intent(in) :: larger_h_denom
352 real :: ent(szi_(g),szk_(g)+1)
354 real :: b1(szi_(g)), d1(szi_(g))
355 real :: c1(szi_(g),szk_(g))
363 integer :: i, j, k, is, ie, js, je, nz, halo
365 halo=0 ;
if (
present(halo_here)) halo = max(halo_here,0)
367 is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo ; nz = gv%ke
369 h_neglect = gv%H_subroundoff
370 kap_dt_x2 = (2.0*kappa_dt)*gv%Z_to_H**2
372 if (
present(larger_h_denom))
then
373 if (larger_h_denom) h0 = 1.0e-16*sqrt(kappa_dt)*gv%Z_to_H
376 if (kap_dt_x2 <= 0.0)
then
378 do k=1,nz ;
do j=js,je ;
do i=is,ie
379 t_f(i,j,k) = t_in(i,j,k) ; s_f(i,j,k) = s_in(i,j,k)
380 enddo ;
enddo ;
enddo
385 ent(i,2) = kap_dt_x2 / ((h(i,j,1)+h(i,j,2)) + h0)
386 h_tr = h(i,j,1) + h_neglect
387 b1(i) = 1.0 / (h_tr + ent(i,2))
389 t_f(i,j,1) = (b1(i)*h_tr)*t_in(i,j,1)
390 s_f(i,j,1) = (b1(i)*h_tr)*s_in(i,j,1)
392 do k=2,nz-1 ;
do i=is,ie
393 ent(i,k+1) = kap_dt_x2 / ((h(i,j,k)+h(i,j,k+1)) + h0)
394 h_tr = h(i,j,k) + h_neglect
395 c1(i,k) = ent(i,k) * b1(i)
396 b1(i) = 1.0 / ((h_tr + d1(i)*ent(i,k)) + ent(i,k+1))
397 d1(i) = b1(i) * (h_tr + d1(i)*ent(i,k))
398 t_f(i,j,k) = b1(i) * (h_tr*t_in(i,j,k) + ent(i,k)*t_f(i,j,k-1))
399 s_f(i,j,k) = b1(i) * (h_tr*s_in(i,j,k) + ent(i,k)*s_f(i,j,k-1))
402 c1(i,nz) = ent(i,nz) * b1(i)
403 h_tr = h(i,j,nz) + h_neglect
404 b1(i) = 1.0 / (h_tr + d1(i)*ent(i,nz))
405 t_f(i,j,nz) = b1(i) * (h_tr*t_in(i,j,nz) + ent(i,nz)*t_f(i,j,nz-1))
406 s_f(i,j,nz) = b1(i) * (h_tr*s_in(i,j,nz) + ent(i,nz)*s_f(i,j,nz-1))
408 do k=nz-1,1,-1 ;
do i=is,ie
409 t_f(i,j,k) = t_f(i,j,k) + c1(i,k+1)*t_f(i,j,k+1)
410 s_f(i,j,k) = s_f(i,j,k) + c1(i,k+1)*s_f(i,j,k+1)
415 end subroutine vert_fill_ts