23 implicit none ;
private
25 #include <MOM_memory.h>
39 real :: khth_slope_cff
47 logical :: use_fgnv_streamfn
56 logical :: detangle_interfaces
64 logical :: use_gme_thickness_diffuse
66 logical :: meke_geometric
68 real :: meke_geometric_alpha
70 real :: meke_geometric_epsilon
72 logical :: use_kh_in_meke
76 real,
pointer :: gmwork(:,:) => null()
77 real,
pointer :: diagslopex(:,:,:) => null()
78 real,
pointer :: diagslopey(:,:,:) => null()
80 real,
dimension(:,:,:),
pointer :: &
86 integer :: id_uhgm = -1, id_vhgm = -1, id_gmwork = -1
87 integer :: id_kh_u = -1, id_kh_v = -1, id_kh_t = -1
88 integer :: id_kh_u1 = -1, id_kh_v1 = -1, id_kh_t1 = -1
89 integer :: id_slope_x = -1, id_slope_y = -1
90 integer :: id_sfn_unlim_x = -1, id_sfn_unlim_y = -1, id_sfn_x = -1, id_sfn_y = -1
99 subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp, CS)
103 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(inout) :: h
104 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(inout) :: uhtr
106 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(inout) :: vhtr
109 real,
intent(in) :: dt
115 real :: e(szi_(g), szj_(g), szk_(g)+1)
117 real :: uhd(szib_(g), szj_(g), szk_(g))
118 real :: vhd(szi_(g), szjb_(g), szk_(g))
120 real,
dimension(SZIB_(G), SZJ_(G), SZK_(G)+1) :: &
126 real,
dimension(SZI_(G), SZJB_(G), SZK_(G)+1) :: &
132 real,
dimension(SZI_(G), SZJ_(G), SZK_(G)) :: &
135 real,
dimension(SZIB_(G), SZJ_(G)) :: &
137 real,
dimension(SZI_(G), SZJB_(G)) :: &
139 real :: khth_loc_u(szib_(g), szj_(g))
140 real :: khth_loc_v(szi_(g), szjb_(g))
141 real :: khth_loc(szib_(g), szjb_(g))
144 real,
dimension(:,:),
pointer :: cg1 => null()
145 logical :: use_varmix, resoln_scaled, depth_scaled, use_stored_slopes, khth_use_ebt_struct, use_visbeck
146 logical :: use_qg_leith
147 integer :: i, j, k, is, ie, js, je, nz
148 real :: hu(szi_(g), szj_(g))
149 real :: hv(szi_(g), szj_(g))
150 real :: kh_u_lay(szi_(g), szj_(g))
151 real :: kh_v_lay(szi_(g), szj_(g))
153 if (.not.
associated(cs))
call mom_error(fatal,
"MOM_thickness_diffuse:"// &
154 "Module must be initialized before it is used.")
156 if ((.not.cs%thickness_diffuse) .or. &
157 .not.( cs%Khth > 0.0 .or.
associated(varmix) .or.
associated(meke) ) )
return
159 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
160 h_neglect = gv%H_subroundoff
162 if (
associated(meke))
then
163 if (
associated(meke%GM_src))
then
164 do j=js,je ;
do i=is,ie ; meke%GM_src(i,j) = 0. ;
enddo ;
enddo
168 use_varmix = .false. ; resoln_scaled = .false. ; use_stored_slopes = .false.
169 khth_use_ebt_struct = .false. ; use_visbeck = .false. ; use_qg_leith = .false.
170 depth_scaled = .false.
172 if (
associated(varmix))
then
173 use_varmix = varmix%use_variable_mixing .and. (cs%KHTH_Slope_Cff > 0.)
174 resoln_scaled = varmix%Resoln_scaled_KhTh
175 depth_scaled = varmix%Depth_scaled_KhTh
176 use_stored_slopes = varmix%use_stored_slopes
177 khth_use_ebt_struct = varmix%khth_use_ebt_struct
178 use_visbeck = varmix%use_Visbeck
179 use_qg_leith = varmix%use_QG_Leith_GM
180 if (
associated(varmix%cg1)) cg1 => varmix%cg1
187 do j=js,je ;
do i=is-1,ie
188 kh_u_cfl(i,j) = (0.25*cs%max_Khth_CFL) / &
189 (dt * (g%IdxCu(i,j)*g%IdxCu(i,j) + g%IdyCu(i,j)*g%IdyCu(i,j)))
192 do j=js-1,je ;
do i=is,ie
193 kh_v_cfl(i,j) = (0.25*cs%max_Khth_CFL) / &
194 (dt * (g%IdxCv(i,j)*g%IdxCv(i,j) + g%IdyCv(i,j)*g%IdyCv(i,j)))
198 call find_eta(h, tv, g, gv, us, e, halo_size=1)
207 do j=js,je;
do i=is-1,ie
208 khth_loc_u(i,j) = cs%Khth
212 if (use_visbeck)
then
214 do j=js,je ;
do i=is-1,ie
215 khth_loc_u(i,j) = khth_loc_u(i,j) + &
216 cs%KHTH_Slope_Cff*varmix%L2u(i,j) * varmix%SN_u(i,j)
221 if (
associated(meke))
then ;
if (
associated(meke%Kh))
then
222 if (cs%MEKE_GEOMETRIC)
then
224 do j=js,je ;
do i=is-1,ie
225 khth_loc_u(i,j) = khth_loc_u(i,j) + g%mask2dCu(i,j) * cs%MEKE_GEOMETRIC_alpha * &
226 0.5*(meke%MEKE(i,j)+meke%MEKE(i+1,j)) / &
227 (varmix%SN_u(i,j) + cs%MEKE_GEOMETRIC_epsilon)
230 do j=js,je ;
do i=is-1,ie
231 khth_loc_u(i,j) = khth_loc_u(i,j) + meke%KhTh_fac*sqrt(meke%Kh(i,j)*meke%Kh(i+1,j))
236 if (resoln_scaled)
then
238 do j=js,je;
do i=is-1,ie
239 khth_loc_u(i,j) = khth_loc_u(i,j) * varmix%Res_fn_u(i,j)
243 if (depth_scaled)
then
245 do j=js,je;
do i=is-1,ie
246 khth_loc_u(i,j) = khth_loc_u(i,j) * varmix%Depth_fn_u(i,j)
250 if (cs%Khth_Max > 0)
then
252 do j=js,je;
do i=is-1,ie
253 khth_loc_u(i,j) = max(cs%Khth_Min, min(khth_loc_u(i,j), cs%Khth_Max))
257 do j=js,je;
do i=is-1,ie
258 khth_loc_u(i,j) = max(cs%Khth_Min, khth_loc_u(i,j))
262 do j=js,je;
do i=is-1,ie
263 kh_u(i,j,1) = min(kh_u_cfl(i,j), khth_loc_u(i,j))
266 if (khth_use_ebt_struct)
then
268 do k=2,nz+1 ;
do j=js,je ;
do i=is-1,ie
269 kh_u(i,j,k) = kh_u(i,j,1) * 0.5 * ( varmix%ebt_struct(i,j,k-1) + varmix%ebt_struct(i+1,j,k-1) )
270 enddo ;
enddo ;
enddo
273 do k=2,nz+1 ;
do j=js,je ;
do i=is-1,ie
274 kh_u(i,j,k) = kh_u(i,j,1)
275 enddo ;
enddo ;
enddo
279 if (use_qg_leith)
then
281 do k=1,nz ;
do j=js,je ;
do i=is-1,ie
282 kh_u(i,j,k) = varmix%KH_u_QG(i,j,k)
283 enddo ;
enddo ;
enddo
287 if (cs%use_GME_thickness_diffuse)
then
289 do k=1,nz+1 ;
do j=js,je ;
do i=is-1,ie
290 cs%KH_u_GME(i,j,k) = kh_u(i,j,k)
291 enddo ;
enddo ;
enddo
295 do j=js-1,je ;
do i=is,ie
296 khth_loc_v(i,j) = cs%Khth
300 if (use_visbeck)
then
302 do j=js-1,je ;
do i=is,ie
303 khth_loc_v(i,j) = khth_loc_v(i,j) + cs%KHTH_Slope_Cff*varmix%L2v(i,j)*varmix%SN_v(i,j)
307 if (
associated(meke))
then ;
if (
associated(meke%Kh))
then
308 if (cs%MEKE_GEOMETRIC)
then
310 do j=js-1,je ;
do i=is,ie
311 khth_loc(i,j) = khth_loc(i,j) + g%mask2dCv(i,j) * cs%MEKE_GEOMETRIC_alpha * &
312 0.5*(meke%MEKE(i,j)+meke%MEKE(i,j+1)) / &
313 (varmix%SN_v(i,j) + cs%MEKE_GEOMETRIC_epsilon)
316 do j=js-1,je ;
do i=is,ie
317 khth_loc_v(i,j) = khth_loc_v(i,j) + meke%KhTh_fac*sqrt(meke%Kh(i,j)*meke%Kh(i,j+1))
322 if (resoln_scaled)
then
324 do j=js-1,je;
do i=is,ie
325 khth_loc_v(i,j) = khth_loc_v(i,j) * varmix%Res_fn_v(i,j)
329 if (depth_scaled)
then
331 do j=js-1,je ;
do i=is,ie
332 khth_loc_v(i,j) = khth_loc_v(i,j) * varmix%Depth_fn_v(i,j)
336 if (cs%Khth_Max > 0)
then
338 do j=js-1,je ;
do i=is,ie
339 khth_loc_v(i,j) = max(cs%Khth_Min, min(khth_loc_v(i,j), cs%Khth_Max))
343 do j=js-1,je ;
do i=is,ie
344 khth_loc_v(i,j) = max(cs%Khth_Min, khth_loc_v(i,j))
348 if (cs%max_Khth_CFL > 0.0)
then
350 do j=js-1,je ;
do i=is,ie
351 kh_v(i,j,1) = min(kh_v_cfl(i,j), khth_loc_v(i,j))
355 if (khth_use_ebt_struct)
then
357 do k=2,nz+1 ;
do j=js-1,je ;
do i=is,ie
358 kh_v(i,j,k) = kh_v(i,j,1) * 0.5 * ( varmix%ebt_struct(i,j,k-1) + varmix%ebt_struct(i,j+1,k-1) )
359 enddo ;
enddo ;
enddo
362 do k=2,nz+1 ;
do j=js-1,je ;
do i=is,ie
363 kh_v(i,j,k) = kh_v(i,j,1)
364 enddo ;
enddo ;
enddo
368 if (use_qg_leith)
then
370 do k=1,nz ;
do j=js-1,je ;
do i=is,ie
371 kh_v(i,j,k) = varmix%KH_v_QG(i,j,k)
372 enddo ;
enddo ;
enddo
376 if (cs%use_GME_thickness_diffuse)
then
378 do k=1,nz+1 ;
do j=js-1,je ;
do i=is,ie
379 cs%KH_v_GME(i,j,k) = kh_v(i,j,k)
380 enddo ;
enddo ;
enddo
383 if (
associated(meke))
then ;
if (
associated(meke%Kh))
then
384 if (cs%MEKE_GEOMETRIC)
then
386 do j=js,je ;
do i=is,ie
388 meke%Kh(i,j) = cs%MEKE_GEOMETRIC_alpha * meke%MEKE(i,j) / &
389 (0.25*(varmix%SN_u(i,j)+varmix%SN_u(i-1,j)+varmix%SN_v(i,j)+varmix%SN_v(i,j-1)) + &
390 cs%MEKE_GEOMETRIC_epsilon)
397 do k=1,nz+1 ;
do j=js,je ;
do i=is-1,ie ; int_slope_u(i,j,k) = 0.0 ;
enddo ;
enddo ;
enddo
399 do k=1,nz+1 ;
do j=js-1,je ;
do i=is,ie ; int_slope_v(i,j,k) = 0.0 ;
enddo ;
enddo ;
enddo
402 if (cs%detangle_interfaces)
then
403 call add_detangling_kh(h, e, kh_u, kh_v, kh_u_cfl, kh_v_cfl, tv, dt, g, gv, us, &
404 cs, int_slope_u, int_slope_v)
408 call uvchksum(
"Kh_[uv]", kh_u, kh_v, g%HI, haloshift=0, scale=us%L_to_m**2*us%s_to_T)
409 call uvchksum(
"int_slope_[uv]", int_slope_u, int_slope_v, g%HI, haloshift=0)
410 call hchksum(h,
"thickness_diffuse_1 h", g%HI, haloshift=1, scale=gv%H_to_m)
411 call hchksum(e,
"thickness_diffuse_1 e", g%HI, haloshift=1, scale=us%Z_to_m)
412 if (use_stored_slopes)
then
413 call uvchksum(
"VarMix%slope_[xy]", varmix%slope_x, varmix%slope_y, &
416 if (
associated(tv%eqn_of_state))
then
417 call hchksum(tv%T,
"thickness_diffuse T", g%HI, haloshift=1)
418 call hchksum(tv%S,
"thickness_diffuse S", g%HI, haloshift=1)
423 if (use_stored_slopes)
then
424 call thickness_diffuse_full(h, e, kh_u, kh_v, tv, uhd, vhd, cg1, dt, g, gv, us, meke, cs, &
425 int_slope_u, int_slope_v, varmix%slope_x, varmix%slope_y)
427 call thickness_diffuse_full(h, e, kh_u, kh_v, tv, uhd, vhd, cg1, dt, g, gv, us, meke, cs, &
428 int_slope_u, int_slope_v)
431 if (
associated(meke) .AND.
associated(varmix))
then
432 if (
associated(meke%Rd_dx_h) .and.
associated(varmix%Rd_dx_h))
then
434 do j=js,je ;
do i=is,ie
435 meke%Rd_dx_h(i,j) = varmix%Rd_dx_h(i,j)
441 if (query_averaging_enabled(cs%diag))
then
442 if (cs%id_uhGM > 0)
call post_data(cs%id_uhGM, uhd, cs%diag)
443 if (cs%id_vhGM > 0)
call post_data(cs%id_vhGM, vhd, cs%diag)
444 if (cs%id_GMwork > 0)
call post_data(cs%id_GMwork, cs%GMwork, cs%diag)
445 if (cs%id_KH_u > 0)
call post_data(cs%id_KH_u, kh_u, cs%diag)
446 if (cs%id_KH_v > 0)
call post_data(cs%id_KH_v, kh_v, cs%diag)
447 if (cs%id_KH_u1 > 0)
call post_data(cs%id_KH_u1, kh_u(:,:,1), cs%diag)
448 if (cs%id_KH_v1 > 0)
call post_data(cs%id_KH_v1, kh_v(:,:,1), cs%diag)
455 if (cs%id_KH_t > 0 .or. cs%id_KH_t1 > 0 .or. cs%Use_KH_in_MEKE)
then
459 do j=js,je ;
do i=is-1,ie
460 hu(i,j) = 2.0*h(i,j,k)*h(i+1,j,k)/(h(i,j,k)+h(i+1,j,k)+h_neglect)
461 if (hu(i,j) /= 0.0) hu(i,j) = 1.0
462 kh_u_lay(i,j) = 0.5*(kh_u(i,j,k)+kh_u(i,j,k+1))
464 do j=js-1,je ;
do i=is,ie
465 hv(i,j) = 2.0*h(i,j,k)*h(i,j+1,k)/(h(i,j,k)+h(i,j+1,k)+h_neglect)
466 if (hv(i,j) /= 0.0) hv(i,j) = 1.0
467 kh_v_lay(i,j) = 0.5*(kh_v(i,j,k)+kh_v(i,j,k+1))
470 do j=js,je ;
do i=is,ie
471 kh_t(i,j,k) = ((hu(i-1,j)*kh_u_lay(i-1,j)+hu(i,j)*kh_u_lay(i,j)) &
472 +(hv(i,j-1)*kh_v_lay(i,j-1)+hv(i,j)*kh_v_lay(i,j))) &
473 / (hu(i-1,j)+hu(i,j)+hv(i,j-1)+hv(i,j)+h_neglect)
477 if (cs%Use_KH_in_MEKE)
then
478 meke%Kh_diff(:,:) = 0.0
480 do j=js,je ;
do i=is,ie
481 meke%Kh_diff(i,j) = meke%Kh_diff(i,j) + kh_t(i,j,k) * h(i,j,k)
485 do j=js,je ;
do i=is,ie
486 meke%Kh_diff(i,j) = meke%Kh_diff(i,j) / max(1.0,g%bathyT(i,j))
490 if (cs%id_KH_t > 0)
call post_data(cs%id_KH_t, kh_t, cs%diag)
491 if (cs%id_KH_t1 > 0)
call post_data(cs%id_KH_t1, kh_t(:,:,1), cs%diag)
498 do j=js,je ;
do i=is-1,ie
499 uhtr(i,j,k) = uhtr(i,j,k) + uhd(i,j,k) * dt
500 if (
associated(cdp%uhGM)) cdp%uhGM(i,j,k) = uhd(i,j,k)
502 do j=js-1,je ;
do i=is,ie
503 vhtr(i,j,k) = vhtr(i,j,k) + vhd(i,j,k) * dt
504 if (
associated(cdp%vhGM)) cdp%vhGM(i,j,k) = vhd(i,j,k)
506 do j=js,je ;
do i=is,ie
507 h(i,j,k) = h(i,j,k) - dt * g%IareaT(i,j) * &
508 ((uhd(i,j,k) - uhd(i-1,j,k)) + (vhd(i,j,k) - vhd(i,j-1,k)))
509 if (h(i,j,k) < gv%Angstrom_H) h(i,j,k) = gv%Angstrom_H
519 call uvchksum(
"thickness_diffuse [uv]hD", uhd, vhd, &
520 g%HI, haloshift=0, scale=gv%H_to_m*us%L_to_m**2*us%s_to_T)
521 call uvchksum(
"thickness_diffuse [uv]htr", uhtr, vhtr, &
522 g%HI, haloshift=0, scale=us%L_to_m**2*gv%H_to_m)
523 call hchksum(h,
"thickness_diffuse h", g%HI, haloshift=0, scale=gv%H_to_m)
531 subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV, US, MEKE, &
532 CS, int_slope_u, int_slope_v, slope_x, slope_y)
536 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
537 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1),
intent(in) :: e
538 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)+1),
intent(in) :: Kh_u
540 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)+1),
intent(in) :: Kh_v
543 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(out) :: uhD
545 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(out) :: vhD
547 real,
dimension(:,:),
pointer :: cg1
548 real,
intent(in) :: dt
551 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)+1),
optional,
intent(in) :: int_slope_u
555 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)+1),
optional,
intent(in) :: int_slope_v
559 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)+1),
optional,
intent(in) :: slope_x
560 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)+1),
optional,
intent(in) :: slope_y
562 real,
dimension(SZI_(G), SZJ_(G), SZK_(G)) :: &
571 real,
dimension(SZI_(G), SZJB_(G), SZK_(G)+1) :: &
575 real,
dimension(SZIB_(G), SZJ_(G), SZK_(G)+1) :: &
579 real,
dimension(SZI_(G), SZJ_(G), SZK_(G)+1) :: &
582 real,
dimension(SZIB_(G)) :: &
585 real,
dimension(SZI_(G)) :: &
588 real :: uhtot(SZIB_(G), SZJ_(G))
589 real :: vhtot(SZI_(G), SZJB_(G))
590 real,
dimension(SZIB_(G)) :: &
594 real,
dimension(SZI_(G)) :: &
598 real :: Work_u(SZIB_(G), SZJ_(G))
599 real :: Work_v(SZI_(G), SZJB_(G))
608 real :: drdi_u(SZIB_(G), SZK_(G)+1)
609 real :: drdj_v(SZI_(G), SZK_(G)+1)
610 real :: drdkDe_u(SZIB_(G),SZK_(G)+1)
612 real :: drdkDe_v(SZI_(G),SZK_(G)+1)
614 real :: hg2A, hg2B, hg2L, hg2R
615 real :: haA, haB, haL, haR
617 real :: wtA, wtB, wtL, wtR
621 real :: c2_h_u(SZIB_(G), SZK_(G)+1)
622 real :: c2_h_v(SZI_(G), SZK_(G)+1)
623 real :: hN2_u(SZIB_(G), SZK_(G)+1)
624 real :: hN2_v(SZI_(G), SZK_(G)+1)
627 real :: Sfn_unlim_u(SZIB_(G), SZK_(G)+1)
628 real :: Sfn_unlim_v(SZI_(G), SZK_(G)+1)
629 real :: slope2_Ratio_u(SZIB_(G), SZK_(G)+1)
630 real :: slope2_Ratio_v(SZI_(G), SZK_(G)+1)
653 real,
dimension(SZIB_(G), SZJ_(G), SZK_(G)+1) :: diag_sfn_x, diag_sfn_unlim_x
654 real,
dimension(SZI_(G), SZJB_(G), SZK_(G)+1) :: diag_sfn_y, diag_sfn_unlim_y
655 logical :: present_int_slope_u, present_int_slope_v
656 logical :: present_slope_x, present_slope_y, calc_derivatives
657 integer :: is, ie, js, je, nz, IsdB
659 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke ; isdb = g%IsdB
662 i_slope_max2 = 1.0 / (cs%slope_max**2)
663 g_scale = gv%g_Earth * gv%H_to_Z
665 h_neglect = gv%H_subroundoff ; h_neglect2 = h_neglect**2
666 dz_neglect = gv%H_subroundoff*gv%H_to_Z
667 g_rho0 = gv%g_Earth / gv%Rho0
668 n2_floor = cs%N2_floor*us%Z_to_L**2
670 use_eos =
associated(tv%eqn_of_state)
671 present_int_slope_u =
PRESENT(int_slope_u)
672 present_int_slope_v =
PRESENT(int_slope_v)
673 present_slope_x =
PRESENT(slope_x)
674 present_slope_y =
PRESENT(slope_y)
676 nk_linear = max(gv%nkml, 1)
678 slope_x_pe(:,:,:) = 0.0
679 slope_y_pe(:,:,:) = 0.0
680 hn2_x_pe(:,:,:) = 0.0
681 hn2_y_pe(:,:,:) = 0.0
684 if (
associated(meke)) find_work =
associated(meke%GM_src)
685 find_work = (
associated(cs%GMwork) .or. find_work)
688 call vert_fill_ts(h, tv%T, tv%S, cs%kappa_smooth*dt, t, s, g, gv, 1, larger_h_denom=.true.)
691 if (cs%use_FGNV_streamfn .and. .not.
associated(cg1))
call mom_error(fatal, &
692 "cg1 must be associated when using FGNV streamfunction.")
699 do j=js-1,je+1 ;
do i=is-1,ie+1
700 h_avail_rsum(i,j,1) = 0.0
703 h_avail(i,j,1) = max(i4dt*g%areaT(i,j)*(h(i,j,1)-gv%Angstrom_H),0.0)
704 h_avail_rsum(i,j,2) = h_avail(i,j,1)
706 pres(i,j,2) = pres(i,j,1) + gv%H_to_Pa*h(i,j,1)
710 do k=2,nz ;
do i=is-1,ie+1
711 h_avail(i,j,k) = max(i4dt*g%areaT(i,j)*(h(i,j,k)-gv%Angstrom_H),0.0)
712 h_avail_rsum(i,j,k+1) = h_avail_rsum(i,j,k) + h_avail(i,j,k)
713 h_frac(i,j,k) = 0.0 ;
if (h_avail(i,j,k) > 0.0) &
714 h_frac(i,j,k) = h_avail(i,j,k) / h_avail_rsum(i,j,k+1)
715 pres(i,j,k+1) = pres(i,j,k) + gv%H_to_Pa*h(i,j,k)
719 do j=js,je ;
do i=is-1,ie
720 uhtot(i,j) = 0.0 ; work_u(i,j) = 0.0
721 diag_sfn_x(i,j,1) = 0.0 ; diag_sfn_unlim_x(i,j,1) = 0.0
722 diag_sfn_x(i,j,nz+1) = 0.0 ; diag_sfn_unlim_x(i,j,nz+1) = 0.0
725 do j=js-1,je ;
do i=is,ie
726 vhtot(i,j) = 0.0 ; work_v(i,j) = 0.0
727 diag_sfn_y(i,j,1) = 0.0 ; diag_sfn_unlim_y(i,j,1) = 0.0
728 diag_sfn_y(i,j,nz+1) = 0.0 ; diag_sfn_unlim_y(i,j,nz+1) = 0.0
746 do i=is-1,ie ; hn2_u(i,1) = 0. ; hn2_u(i,nz+1) = 0. ;
enddo
748 if (find_work .and. .not.(use_eos))
then
749 drdia = 0.0 ; drdib = 0.0
750 drdkl = gv%Rlay(k) - gv%Rlay(k-1) ; drdkr = drdkl
753 calc_derivatives = use_eos .and. (k >= nk_linear) .and. &
754 (find_work .or. .not. present_slope_x .or. cs%use_FGNV_streamfn)
757 if (calc_derivatives)
then
759 pres_u(i) = 0.5*(pres(i,j,k) + pres(i+1,j,k))
760 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)))
761 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)))
764 drho_ds_u, (is-isdb+1)-1, ie-is+2, tv%eqn_of_state, scale=us%kg_m3_to_R)
768 if (calc_derivatives)
then
770 drdia = drho_dt_u(i) * (t(i+1,j,k-1)-t(i,j,k-1)) + &
771 drho_ds_u(i) * (s(i+1,j,k-1)-s(i,j,k-1))
772 drdib = drho_dt_u(i) * (t(i+1,j,k)-t(i,j,k)) + &
773 drho_ds_u(i) * (s(i+1,j,k)-s(i,j,k))
776 drdkl = (drho_dt_u(i) * (t(i,j,k)-t(i,j,k-1)) + &
777 drho_ds_u(i) * (s(i,j,k)-s(i,j,k-1)))
778 drdkr = (drho_dt_u(i) * (t(i+1,j,k)-t(i+1,j,k-1)) + &
779 drho_ds_u(i) * (s(i+1,j,k)-s(i+1,j,k-1)))
780 drdkde_u(i,k) = drdkr * e(i+1,j,k) - drdkl * e(i,j,k)
781 elseif (find_work)
then
782 drdkde_u(i,k) = drdkr * e(i+1,j,k) - drdkl * e(i,j,k)
785 if (find_work) drdi_u(i,k) = drdib
787 if (k > nk_linear)
then
789 if (cs%use_FGNV_streamfn .or. find_work .or. .not.present_slope_x)
then
790 hg2l = h(i,j,k-1)*h(i,j,k) + h_neglect2
791 hg2r = h(i+1,j,k-1)*h(i+1,j,k) + h_neglect2
792 hal = 0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect
793 har = 0.5*(h(i+1,j,k-1) + h(i+1,j,k)) + h_neglect
794 if (gv%Boussinesq)
then
795 dzal = hal * gv%H_to_Z ; dzar = har * gv%H_to_Z
797 dzal = 0.5*(e(i,j,k-1) - e(i,j,k+1)) + dz_neglect
798 dzar = 0.5*(e(i+1,j,k-1) - e(i+1,j,k+1)) + dz_neglect
802 wtl = hg2l*(har*dzar) ; wtr = hg2r*(hal*dzal)
804 drdz = (wtl * drdkl + wtr * drdkr) / (dzal*wtl + dzar*wtr)
808 hg2a = h(i,j,k-1)*h(i+1,j,k-1) + h_neglect2
809 hg2b = h(i,j,k)*h(i+1,j,k) + h_neglect2
810 haa = 0.5*(h(i,j,k-1) + h(i+1,j,k-1)) + h_neglect
811 hab = 0.5*(h(i,j,k) + h(i+1,j,k)) + h_neglect
814 hn2_u(i,k) = (0.5 * gv%H_to_Z * ( hg2a / haa + hg2b / hab )) * &
815 max(drdz*g_rho0, n2_floor)
817 if (present_slope_x)
then
818 slope = slope_x(i,j,k)
819 slope2_ratio_u(i,k) = slope**2 * i_slope_max2
823 wta = hg2a*hab ; wtb = hg2b*haa
825 drdx = ((wta * drdia + wtb * drdib) / (wta + wtb) - &
826 drdz * (e(i,j,k)-e(i+1,j,k))) * g%IdxCu(i,j)
830 mag_grad2 = drdx**2 + (us%L_to_Z*drdz)**2
831 if (mag_grad2 > 0.0)
then
832 slope = drdx / sqrt(mag_grad2)
833 slope2_ratio_u(i,k) = slope**2 * i_slope_max2
836 slope2_ratio_u(i,k) = 1.0e20
842 if (present_int_slope_u)
then
843 slope = (1.0 - int_slope_u(i,j,k)) * slope + &
844 int_slope_u(i,j,k) * us%Z_to_L*((e(i+1,j,k)-e(i,j,k)) * g%IdxCu(i,j))
845 slope2_ratio_u(i,k) = (1.0 - int_slope_u(i,j,k)) * slope2_ratio_u(i,k)
848 slope_x_pe(i,j,k) = min(slope,cs%slope_max)
849 hn2_x_pe(i,j,k) = hn2_u(i,k)
850 if (cs%id_slope_x > 0) cs%diagSlopeX(i,j,k) = slope
853 sfn_unlim_u(i,k) = -((kh_u(i,j,k)*g%dy_Cu(i,j))*us%L_to_Z*slope)
857 if (sfn_unlim_u(i,k) > 0.0)
then
858 if (e(i,j,k) < e(i+1,j,nz+1))
then
859 sfn_unlim_u(i,k) = 0.0
861 elseif (e(i+1,j,nz+1) > e(i,j,k+1))
then
864 sfn_unlim_u(i,k) = sfn_unlim_u(i,k) * ((e(i,j,k) - e(i+1,j,nz+1)) / &
865 ((e(i,j,k) - e(i,j,k+1)) + dz_neglect))
868 if (e(i+1,j,k) < e(i,j,nz+1))
then ; sfn_unlim_u(i,k) = 0.0
869 elseif (e(i,j,nz+1) > e(i+1,j,k+1))
then
870 sfn_unlim_u(i,k) = sfn_unlim_u(i,k) * ((e(i+1,j,k) - e(i,j,nz+1)) / &
871 ((e(i+1,j,k) - e(i+1,j,k+1)) + dz_neglect))
876 if (present_slope_x)
then
877 slope = slope_x(i,j,k)
879 slope = us%Z_to_L*((e(i,j,k)-e(i+1,j,k))*g%IdxCu(i,j)) * g%mask2dCu(i,j)
881 if (cs%id_slope_x > 0) cs%diagSlopeX(i,j,k) = slope
882 sfn_unlim_u(i,k) = ((kh_u(i,j,k)*g%dy_Cu(i,j))*us%L_to_Z*slope)
883 hn2_u(i,k) = gv%g_prime(k)
886 hn2_u(i,k) = n2_floor * dz_neglect
887 sfn_unlim_u(i,k) = 0.
889 if (cs%id_sfn_unlim_x>0) diag_sfn_unlim_x(i,j,k) = sfn_unlim_u(i,k)
893 if (cs%use_FGNV_streamfn)
then
894 do k=1,nz ;
do i=is-1,ie ;
if (g%mask2dCu(i,j)>0.)
then
895 h_harm = max( h_neglect, &
896 2. * h(i,j,k) * h(i+1,j,k) / ( ( h(i,j,k) + h(i+1,j,k) ) + h_neglect ) )
897 c2_h_u(i,k) = cs%FGNV_scale * &
898 ( 0.5*( cg1(i,j) + cg1(i+1,j) ) )**2 / (gv%H_to_Z*h_harm)
899 endif ;
enddo ;
enddo
903 if (g%mask2dCu(i,j)>0.)
then
904 sfn_unlim_u(i,:) = ( 1. + cs%FGNV_scale ) * sfn_unlim_u(i,:)
907 sfn_unlim_u(i,:) = 0.
914 if (k > nk_linear)
then
917 if (uhtot(i,j) <= 0.0)
then
919 sfn_safe = uhtot(i,j) * (1.0 - h_frac(i,j,k)) * gv%H_to_Z
921 sfn_safe = uhtot(i,j) * (1.0 - h_frac(i+1,j,k)) * gv%H_to_Z
925 sfn_est = (sfn_unlim_u(i,k) + slope2_ratio_u(i,k)*sfn_safe) / (1.0 + slope2_ratio_u(i,k))
927 sfn_est = sfn_unlim_u(i,k)
932 sfn_in_h = min(max(sfn_est * gv%Z_to_H, -h_avail_rsum(i,j,k)), h_avail_rsum(i+1,j,k))
936 uhd(i,j,k) = max(min((sfn_in_h - uhtot(i,j)), h_avail(i,j,k)), &
939 if (cs%id_sfn_x>0) diag_sfn_x(i,j,k) = diag_sfn_x(i,j,k+1) + uhd(i,j,k)
951 if (uhtot(i,j) <= 0.0)
then
952 uhd(i,j,k) = -uhtot(i,j) * h_frac(i,j,k)
954 uhd(i,j,k) = -uhtot(i,j) * h_frac(i+1,j,k)
957 if (cs%id_sfn_x>0) diag_sfn_x(i,j,k) = diag_sfn_x(i,j,k+1) + uhd(i,j,k)
965 uhtot(i,j) = uhtot(i,j) + uhd(i,j,k)
974 work_u(i,j) = work_u(i,j) + g_scale * &
975 ( uhtot(i,j) * drdkde_u(i,k) - &
976 (uhd(i,j,k) * drdi_u(i,k)) * 0.25 * &
977 ((e(i,j,k) + e(i,j,k+1)) + (e(i+1,j,k) + e(i+1,j,k+1))) )
1000 if (find_work .and. .not.(use_eos))
then
1001 drdja = 0.0 ; drdjb = 0.0
1002 drdkl = gv%Rlay(k) - gv%Rlay(k-1) ; drdkr = drdkl
1005 calc_derivatives = use_eos .and. (k >= nk_linear) .and. &
1006 (find_work .or. .not. present_slope_y)
1008 if (calc_derivatives)
then
1010 pres_v(i) = 0.5*(pres(i,j,k) + pres(i,j+1,k))
1011 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)))
1012 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)))
1015 drho_ds_v, is, ie-is+1, tv%eqn_of_state, scale=us%kg_m3_to_R)
1018 if (calc_derivatives)
then
1020 drdja = drho_dt_v(i) * (t(i,j+1,k-1)-t(i,j,k-1)) + &
1021 drho_ds_v(i) * (s(i,j+1,k-1)-s(i,j,k-1))
1022 drdjb = drho_dt_v(i) * (t(i,j+1,k)-t(i,j,k)) + &
1023 drho_ds_v(i) * (s(i,j+1,k)-s(i,j,k))
1026 drdkl = (drho_dt_v(i) * (t(i,j,k)-t(i,j,k-1)) + &
1027 drho_ds_v(i) * (s(i,j,k)-s(i,j,k-1)))
1028 drdkr = (drho_dt_v(i) * (t(i,j+1,k)-t(i,j+1,k-1)) + &
1029 drho_ds_v(i) * (s(i,j+1,k)-s(i,j+1,k-1)))
1030 drdkde_v(i,k) = drdkr * e(i,j+1,k) - drdkl * e(i,j,k)
1031 elseif (find_work)
then
1032 drdkde_v(i,k) = drdkr * e(i,j+1,k) - drdkl * e(i,j,k)
1035 if (find_work) drdj_v(i,k) = drdjb
1037 if (k > nk_linear)
then
1039 if (cs%use_FGNV_streamfn .or. find_work .or. .not. present_slope_y)
then
1040 hg2l = h(i,j,k-1)*h(i,j,k) + h_neglect2
1041 hg2r = h(i,j+1,k-1)*h(i,j+1,k) + h_neglect2
1042 hal = 0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect
1043 har = 0.5*(h(i,j+1,k-1) + h(i,j+1,k)) + h_neglect
1044 if (gv%Boussinesq)
then
1045 dzal = hal * gv%H_to_Z ; dzar = har * gv%H_to_Z
1047 dzal = 0.5*(e(i,j,k-1) - e(i,j,k+1)) + dz_neglect
1048 dzar = 0.5*(e(i,j+1,k-1) - e(i,j+1,k+1)) + dz_neglect
1052 wtl = hg2l*(har*dzar) ; wtr = hg2r*(hal*dzal)
1054 drdz = (wtl * drdkl + wtr * drdkr) / (dzal*wtl + dzar*wtr)
1058 hg2a = h(i,j,k-1)*h(i,j+1,k-1) + h_neglect2
1059 hg2b = h(i,j,k)*h(i,j+1,k) + h_neglect2
1060 haa = 0.5*(h(i,j,k-1) + h(i,j+1,k-1)) + h_neglect
1061 hab = 0.5*(h(i,j,k) + h(i,j+1,k)) + h_neglect
1064 hn2_v(i,k) = (0.5 * gv%H_to_Z * ( hg2a / haa + hg2b / hab )) * &
1065 max(drdz*g_rho0, n2_floor)
1067 if (present_slope_y)
then
1068 slope = slope_y(i,j,k)
1069 slope2_ratio_v(i,k) = slope**2 * i_slope_max2
1073 wta = hg2a*hab ; wtb = hg2b*haa
1075 drdy = ((wta * drdja + wtb * drdjb) / (wta + wtb) - &
1076 drdz * (e(i,j,k)-e(i,j+1,k))) * g%IdyCv(i,j)
1080 mag_grad2 = drdy**2 + (us%L_to_Z*drdz)**2
1081 if (mag_grad2 > 0.0)
then
1082 slope = drdy / sqrt(mag_grad2)
1083 slope2_ratio_v(i,k) = slope**2 * i_slope_max2
1086 slope2_ratio_v(i,k) = 1.0e20
1092 if (present_int_slope_v)
then
1093 slope = (1.0 - int_slope_v(i,j,k)) * slope + &
1094 int_slope_v(i,j,k) * us%Z_to_L*((e(i,j+1,k)-e(i,j,k)) * g%IdyCv(i,j))
1095 slope2_ratio_v(i,k) = (1.0 - int_slope_v(i,j,k)) * slope2_ratio_v(i,k)
1098 slope_y_pe(i,j,k) = min(slope,cs%slope_max)
1099 hn2_y_pe(i,j,k) = hn2_v(i,k)
1100 if (cs%id_slope_y > 0) cs%diagSlopeY(i,j,k) = slope
1103 sfn_unlim_v(i,k) = -((kh_v(i,j,k)*g%dx_Cv(i,j))*us%L_to_Z*slope)
1107 if (sfn_unlim_v(i,k) > 0.0)
then
1108 if (e(i,j,k) < e(i,j+1,nz+1))
then
1109 sfn_unlim_v(i,k) = 0.0
1111 elseif (e(i,j+1,nz+1) > e(i,j,k+1))
then
1114 sfn_unlim_v(i,k) = sfn_unlim_v(i,k) * ((e(i,j,k) - e(i,j+1,nz+1)) / &
1115 ((e(i,j,k) - e(i,j,k+1)) + dz_neglect))
1118 if (e(i,j+1,k) < e(i,j,nz+1))
then ; sfn_unlim_v(i,k) = 0.0
1119 elseif (e(i,j,nz+1) > e(i,j+1,k+1))
then
1120 sfn_unlim_v(i,k) = sfn_unlim_v(i,k) * ((e(i,j+1,k) - e(i,j,nz+1)) / &
1121 ((e(i,j+1,k) - e(i,j+1,k+1)) + dz_neglect))
1126 if (present_slope_y)
then
1127 slope = slope_y(i,j,k)
1129 slope = us%Z_to_L*((e(i,j,k)-e(i,j+1,k))*g%IdyCv(i,j)) * g%mask2dCv(i,j)
1131 if (cs%id_slope_y > 0) cs%diagSlopeY(i,j,k) = slope
1132 sfn_unlim_v(i,k) = ((kh_v(i,j,k)*g%dx_Cv(i,j))*us%L_to_Z*slope)
1133 hn2_v(i,k) = gv%g_prime(k)
1136 hn2_v(i,k) = n2_floor * dz_neglect
1137 sfn_unlim_v(i,k) = 0.
1139 if (cs%id_sfn_unlim_y>0) diag_sfn_unlim_y(i,j,k) = sfn_unlim_v(i,k)
1143 if (cs%use_FGNV_streamfn)
then
1144 do k=1,nz ;
do i=is,ie ;
if (g%mask2dCv(i,j)>0.)
then
1145 h_harm = max( h_neglect, &
1146 2. * h(i,j,k) * h(i,j+1,k) / ( ( h(i,j,k) + h(i,j+1,k) ) + h_neglect ) )
1147 c2_h_v(i,k) = cs%FGNV_scale * &
1148 ( 0.5*( cg1(i,j) + cg1(i,j+1) ) )**2 / (gv%H_to_Z*h_harm)
1149 endif ;
enddo ;
enddo
1153 if (g%mask2dCv(i,j)>0.)
then
1154 sfn_unlim_v(i,:) = ( 1. + cs%FGNV_scale ) * sfn_unlim_v(i,:)
1157 sfn_unlim_v(i,:) = 0.
1164 if (k > nk_linear)
then
1167 if (vhtot(i,j) <= 0.0)
then
1169 sfn_safe = vhtot(i,j) * (1.0 - h_frac(i,j,k)) * gv%H_to_Z
1171 sfn_safe = vhtot(i,j) * (1.0 - h_frac(i,j+1,k)) * gv%H_to_Z
1175 sfn_est = (sfn_unlim_v(i,k) + slope2_ratio_v(i,k)*sfn_safe) / (1.0 + slope2_ratio_v(i,k))
1177 sfn_est = sfn_unlim_v(i,k)
1182 sfn_in_h = min(max(sfn_est * gv%Z_to_H, -h_avail_rsum(i,j,k)), h_avail_rsum(i,j+1,k))
1186 vhd(i,j,k) = max(min((sfn_in_h - vhtot(i,j)), h_avail(i,j,k)), -h_avail(i,j+1,k))
1188 if (cs%id_sfn_y>0) diag_sfn_y(i,j,k) = diag_sfn_y(i,j,k+1) + vhd(i,j,k)
1200 if (vhtot(i,j) <= 0.0)
then
1201 vhd(i,j,k) = -vhtot(i,j) * h_frac(i,j,k)
1203 vhd(i,j,k) = -vhtot(i,j) * h_frac(i,j+1,k)
1206 if (cs%id_sfn_y>0) diag_sfn_y(i,j,k) = diag_sfn_y(i,j,k+1) + vhd(i,j,k)
1214 vhtot(i,j) = vhtot(i,j) + vhd(i,j,k)
1223 work_v(i,j) = work_v(i,j) + g_scale * &
1224 ( vhtot(i,j) * drdkde_v(i,k) - &
1225 (vhd(i,j,k) * drdj_v(i,k)) * 0.25 * &
1226 ((e(i,j,k) + e(i,j,k+1)) + (e(i,j+1,k) + e(i,j+1,k+1))) )
1234 if (.not.find_work .or. .not.(use_eos))
then
1235 do j=js,je ;
do i=is-1,ie ; uhd(i,j,1) = -uhtot(i,j) ;
enddo ;
enddo
1236 do j=js-1,je ;
do i=is,ie ; vhd(i,j,1) = -vhtot(i,j) ;
enddo ;
enddo
1242 pres_u(i) = 0.5*(pres(i,j,1) + pres(i+1,j,1))
1243 t_u(i) = 0.5*(t(i,j,1) + t(i+1,j,1))
1244 s_u(i) = 0.5*(s(i,j,1) + s(i+1,j,1))
1247 drho_ds_u, (is-isdb+1)-1, ie-is+2, tv%eqn_of_state, scale=us%kg_m3_to_R)
1250 uhd(i,j,1) = -uhtot(i,j)
1253 drdib = drho_dt_u(i) * (t(i+1,j,1)-t(i,j,1)) + &
1254 drho_ds_u(i) * (s(i+1,j,1)-s(i,j,1))
1256 work_u(i,j) = work_u(i,j) + g_scale * &
1257 ( (uhd(i,j,1) * drdib) * 0.25 * &
1258 ((e(i,j,1) + e(i,j,2)) + (e(i+1,j,1) + e(i+1,j,2))) )
1267 pres_v(i) = 0.5*(pres(i,j,1) + pres(i,j+1,1))
1268 t_v(i) = 0.5*(t(i,j,1) + t(i,j+1,1))
1269 s_v(i) = 0.5*(s(i,j,1) + s(i,j+1,1))
1272 drho_ds_v, is, ie-is+1, tv%eqn_of_state, scale=us%kg_m3_to_R)
1275 vhd(i,j,1) = -vhtot(i,j)
1278 drdjb = drho_dt_v(i) * (t(i,j+1,1)-t(i,j,1)) + &
1279 drho_ds_v(i) * (s(i,j+1,1)-s(i,j,1))
1281 work_v(i,j) = work_v(i,j) - g_scale * &
1282 ( (vhd(i,j,1) * drdjb) * 0.25 * &
1283 ((e(i,j,1) + e(i,j,2)) + (e(i,j+1,1) + e(i,j+1,2))) )
1288 if (find_work)
then ;
do j=js,je ;
do i=is,ie
1290 work_h = 0.5 * g%IareaT(i,j) * &
1291 ((work_u(i-1,j) + work_u(i,j)) + (work_v(i,j-1) + work_v(i,j)))
1292 if (
associated(cs%GMwork)) cs%GMwork(i,j) = work_h
1293 if (
associated(meke) .and. .not.cs%GM_src_alt)
then ;
if (
associated(meke%GM_src))
then
1294 meke%GM_src(i,j) = meke%GM_src(i,j) + work_h
1296 enddo ;
enddo ;
endif
1298 if (find_work .and. cs%GM_src_alt .and.
associated(meke))
then ;
if (
associated(meke%GM_src))
then
1299 do j=js,je ;
do i=is,ie ;
do k=nz,1,-1
1300 pe_release_h = -0.25*(kh_u(i,j,k)*(slope_x_pe(i,j,k)**2) * hn2_x_pe(i,j,k) + &
1301 kh_u(i-1,j,k)*(slope_x_pe(i-1,j,k)**2) * hn2_x_pe(i-1,j,k) + &
1302 kh_v(i,j,k)*(slope_y_pe(i,j,k)**2) * hn2_y_pe(i,j,k) + &
1303 kh_v(i,j-1,k)*(slope_y_pe(i,j-1,k)**2) * hn2_y_pe(i,j-1,k))
1304 meke%GM_src(i,j) = meke%GM_src(i,j) + us%L_to_Z**2 * gv%Rho0 * pe_release_h
1305 enddo ;
enddo ;
enddo
1308 if (cs%id_slope_x > 0)
call post_data(cs%id_slope_x, cs%diagSlopeX, cs%diag)
1309 if (cs%id_slope_y > 0)
call post_data(cs%id_slope_y, cs%diagSlopeY, cs%diag)
1310 if (cs%id_sfn_x > 0)
call post_data(cs%id_sfn_x, diag_sfn_x, cs%diag)
1311 if (cs%id_sfn_y > 0)
call post_data(cs%id_sfn_y, diag_sfn_y, cs%diag)
1312 if (cs%id_sfn_unlim_x > 0)
call post_data(cs%id_sfn_unlim_x, diag_sfn_unlim_x, cs%diag)
1313 if (cs%id_sfn_unlim_y > 0)
call post_data(cs%id_sfn_unlim_y, diag_sfn_unlim_y, cs%diag)
1319 integer,
intent(in) :: nk
1320 real,
dimension(nk),
intent(in) :: c2_h
1321 real,
dimension(nk+1),
intent(in) :: hN2
1322 real,
dimension(nk+1),
intent(inout) :: sfn
1328 real :: b_denom, beta, d1, c1(nk)
1331 b_denom = hn2(2) + c2_h(1)
1332 beta = 1.0 / ( b_denom + c2_h(2) )
1334 sfn(2) = ( beta * hn2(2) )*sfn(2)
1336 c1(k-1) = beta * c2_h(k-1)
1337 b_denom = hn2(k) + d1*c2_h(k-1)
1338 beta = 1.0 / (b_denom + c2_h(k))
1340 sfn(k) = beta * (hn2(k)*sfn(k) + c2_h(k-1)*sfn(k-1))
1342 c1(nk) = beta * c2_h(nk)
1345 sfn(k) = sfn(k) + c1(k)*sfn(k+1)
1351 subroutine add_detangling_kh(h, e, Kh_u, Kh_v, KH_u_CFL, KH_v_CFL, tv, dt, G, GV, US, CS, &
1352 int_slope_u, int_slope_v)
1356 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
1357 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1),
intent(in) :: e
1358 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)+1),
intent(inout) :: Kh_u
1360 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)+1),
intent(inout) :: Kh_v
1362 real,
dimension(SZIB_(G),SZJ_(G)),
intent(in) :: Kh_u_CFL
1364 real,
dimension(SZI_(G),SZJB_(G)),
intent(in) :: Kh_v_CFL
1367 real,
intent(in) :: dt
1369 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)+1),
intent(inout) :: int_slope_u
1373 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)+1),
intent(inout) :: int_slope_v
1378 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
1381 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: &
1384 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)) :: &
1387 real,
dimension(SZI_(G),SZJ_(G)) :: &
1420 real :: denom, I_denom
1429 real,
dimension(SZIB_(G),SZK_(G)+1) :: &
1452 real,
dimension(SZIB_(G)) :: &
1454 logical,
dimension(SZIB_(G)) :: &
1456 integer :: i, j, k, n, ish, jsh, is, ie, js, je, nz, k_top
1457 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1459 k_top = gv%nk_rho_varies + 1
1460 h_neglect = gv%H_subroundoff
1465 if (cs%detangle_time > dt) kh_scale = 0.5 * dt / cs%detangle_time
1467 do j=js-1,je+1 ;
do i=is-1,ie+1
1468 de_top(i,j,k_top) = 0.0 ; de_bot(i,j) = 0.0
1470 do k=k_top+1,nz ;
do j=js-1,je+1 ;
do i=is-1,ie+1
1471 de_top(i,j,k) = de_top(i,j,k-1) + h(i,j,k-1)
1472 enddo ;
enddo ;
enddo
1474 do j=js,je ;
do i=is-1,ie
1475 kh_lay_u(i,j,nz) = 0.0 ; kh_lay_u(i,j,k_top) = 0.0
1477 do j=js-1,je ;
do i=is,ie
1478 kh_lay_v(i,j,nz) = 0.0 ; kh_lay_v(i,j,k_top) = 0.0
1481 do k=nz-1,k_top+1,-1
1483 do j=js-1,je+1 ;
do i=is-1,ie+1
1484 de_bot(i,j) = de_bot(i,j) + h(i,j,k+1)
1487 do j=js,je ;
do i=is-1,ie ;
if (g%mask2dCu(i,j) > 0.0)
then
1488 if (h(i,j,k) > h(i+1,j,k))
then
1490 h1 = max( h(i+1,j,k), h2 - min(de_bot(i+1,j), de_top(i+1,j,k)) )
1493 h1 = max( h(i,j,k), h2 - min(de_bot(i,j), de_top(i,j,k)) )
1495 jag_rat = (h2 - h1)**2 / (h2 + h1 + h_neglect)**2
1496 kh_lay_u(i,j,k) = (kh_scale * kh_u_cfl(i,j)) * jag_rat**2
1497 endif ;
enddo ;
enddo
1499 do j=js-1,je ;
do i=is,ie ;
if (g%mask2dCv(i,j) > 0.0)
then
1500 if (h(i,j,k) > h(i,j+1,k))
then
1502 h1 = max( h(i,j+1,k), h2 - min(de_bot(i,j+1), de_top(i,j+1,k)) )
1505 h1 = max( h(i,j,k), h2 - min(de_bot(i,j), de_top(i,j,k)) )
1507 jag_rat = (h2 - h1)**2 / (h2 + h1 + h_neglect)**2
1508 kh_lay_v(i,j,k) = (kh_scale * kh_v_cfl(i,j)) * jag_rat**2
1509 endif ;
enddo ;
enddo
1514 i_4t = kh_scale / (4.0 * dt)
1517 if (n==1)
then ; jsh = js ; ish = is-1
1518 else ; jsh = js-1 ; ish = is ;
endif
1525 do_i(i) = (g%mask2dCu(i,j) > 0.0)
1526 kh_max_max(i) = kh_u_cfl(i,j)
1528 do k=1,nz+1 ;
do i=ish,ie
1529 kh_bg(i,k) = kh_u(i,j,k) ; kh(i,k) = kh_bg(i,k)
1530 kh_min_max_p(i,k) = kh_bg(i,k) ; kh_min_max_m(i,k) = kh_bg(i,k)
1531 kh_detangle(i,k) = 0.0
1535 do_i(i) = (g%mask2dCv(i,j) > 0.0) ; kh_max_max(i) = kh_v_cfl(i,j)
1537 do k=1,nz+1 ;
do i=ish,ie
1538 kh_bg(i,k) = kh_v(i,j,k) ; kh(i,k) = kh_bg(i,k)
1539 kh_min_max_p(i,k) = kh_bg(i,k) ; kh_min_max_m(i,k) = kh_bg(i,k)
1540 kh_detangle(i,k) = 0.0
1545 do k=k_top,nz ;
do i=ish,ie ;
if (do_i(i))
then
1548 denom = ((g%IareaT(i+1,j) + g%IareaT(i,j)) * g%dy_Cu(i,j))
1552 dh = i_4t * ((e(i+1,j,k) - e(i+1,j,k+1)) - &
1553 (e(i,j,k) - e(i,j,k+1))) / denom
1557 sign = 1.0 ;
if (dh < 0) sign = -1.0
1558 sl_k = sign * (e(i+1,j,k)-e(i,j,k)) * g%IdxCu(i,j)
1559 sl_kp1 = sign * (e(i+1,j,k+1)-e(i,j,k+1)) * g%IdxCu(i,j)
1564 denom = (sl_k**2 + sl_kp1**2)
1565 wt1 = 0.5 ; wt2 = 0.5
1566 if (denom > 0.0)
then
1567 wt1 = sl_k**2 / denom ; wt2 = sl_kp1**2 / denom
1569 kh_detangle(i,k) = kh_detangle(i,k) + wt1*kh_lay_u(i,j,k)
1570 kh_detangle(i,k+1) = kh_detangle(i,k+1) + wt2*kh_lay_u(i,j,k)
1573 denom = ((g%IareaT(i,j+1) + g%IareaT(i,j)) * g%dx_Cv(i,j))
1575 dh = i_4t * ((e(i,j+1,k) - e(i,j+1,k+1)) - &
1576 (e(i,j,k) - e(i,j,k+1))) / denom
1580 sign = 1.0 ;
if (dh < 0) sign = -1.0
1581 sl_k = sign * (e(i,j+1,k)-e(i,j,k)) * g%IdyCv(i,j)
1582 sl_kp1 = sign * (e(i,j+1,k+1)-e(i,j,k+1)) * g%IdyCv(i,j)
1587 denom = (sl_k**2 + sl_kp1**2)
1588 wt1 = 0.5 ; wt2 = 0.5
1589 if (denom > 0.0)
then
1590 wt1 = sl_k**2 / denom ; wt2 = sl_kp1**2 / denom
1592 kh_detangle(i,k) = kh_detangle(i,k) + wt1*kh_lay_v(i,j,k)
1593 kh_detangle(i,k+1) = kh_detangle(i,k+1) + wt2*kh_lay_v(i,j,k)
1596 if (adh == 0.0)
then
1597 kh_min_m(i,k+1) = 1.0 ; kh0_min_m(i,k+1) = 0.0
1598 kh_max_m(i,k+1) = 1.0 ; kh0_max_m(i,k+1) = 0.0
1599 kh_min_p(i,k) = 1.0 ; kh0_min_p(i,k) = 0.0
1600 kh_max_p(i,k) = 1.0 ; kh0_max_p(i,k) = 0.0
1601 elseif (adh > 0.0)
then
1602 if (sl_k <= sl_kp1)
then
1605 kh_min_m(i,k+1) = 1.0 ; kh0_min_m(i,k+1) = 0.0
1606 kh_max_m(i,k+1) = 1.0 ; kh0_max_m(i,k+1) = 0.0
1607 kh_min_p(i,k) = 1.0 ; kh0_min_p(i,k) = 0.0
1608 kh_max_p(i,k) = 1.0 ; kh0_max_p(i,k) = 0.0
1609 elseif (sl_k <= 0.0)
then
1610 i_sl = -1.0 / sl_kp1
1612 irsl = 1e9 ;
if (rsl > 1e-9) irsl = 1.0/rsl
1615 if (kh_max_max(i) > 0) &
1616 fn_r = min(sqrt(rsl), rsl + (adh * i_sl) / (kh_max_max(i)))
1618 kh_min_m(i,k+1) = fn_r ; kh0_min_m(i,k+1) = 0.0
1619 kh_max_m(i,k+1) = rsl ; kh0_max_m(i,k+1) = adh * i_sl
1620 kh_min_p(i,k) = irsl ; kh0_min_p(i,k) = -adh * (i_sl*irsl)
1621 kh_max_p(i,k) = 1.0/(fn_r + 1.0e-30) ; kh0_max_p(i,k) = 0.0
1622 elseif (sl_kp1 < 0.0)
then
1623 i_sl_k = 1e18*us%Z_to_L ;
if (sl_k > 1e-18*us%L_to_Z) i_sl_k = 1.0 / sl_k
1624 i_sl_kp1 = 1e18*us%Z_to_L ;
if (-sl_kp1 > 1e-18*us%L_to_Z) i_sl_kp1 = -1.0 / sl_kp1
1626 kh_min_m(i,k+1) = 0.0 ; kh0_min_m(i,k+1) = 0.0
1627 kh_max_m(i,k+1) = - sl_k*i_sl_kp1 ; kh0_max_m(i,k+1) = adh*i_sl_kp1
1628 kh_min_p(i,k) = 0.0 ; kh0_min_p(i,k) = 0.0
1629 kh_max_p(i,k) = sl_kp1*i_sl_k ; kh0_max_p(i,k) = adh*i_sl_k
1633 kh_max = adh / (sl_k - sl_kp1)
1634 kh_min_max_p(i,k) = max(kh_min_max_p(i,k), kh_max)
1635 kh_min_max_m(i,k+1) = max(kh_min_max_m(i,k+1), kh_max)
1639 irsl = 1e9 ;
if (rsl > 1e-9) irsl = 1.0/rsl
1643 if (kh_max_max(i) > 0) &
1644 fn_r = min(sqrt(rsl), rsl + (adh * i_sl) / kh_max_max(i))
1646 kh_min_m(i,k+1) = irsl ; kh0_min_m(i,k+1) = -adh * (i_sl*irsl)
1647 kh_max_m(i,k+1) = 1.0/(fn_r + 1.0e-30) ; kh0_max_m(i,k+1) = 0.0
1648 kh_min_p(i,k) = fn_r ; kh0_min_p(i,k) = 0.0
1649 kh_max_p(i,k) = rsl ; kh0_max_p(i,k) = adh * i_sl
1652 endif ;
enddo ;
enddo
1654 do k=k_top,nz+1,nz+1-k_top ;
do i=ish,ie ;
if (do_i(i))
then
1656 kh_min_m(i,k) = 0.0 ; kh0_min_m(i,k) = 0.0
1657 kh_max_m(i,k) = 0.0 ; kh0_max_m(i,k) = 0.0
1658 kh_min_p(i,k) = 0.0 ; kh0_min_p(i,k) = 0.0
1659 kh_max_p(i,k) = 0.0 ; kh0_max_p(i,k) = 0.0
1660 kh_min_max_p(i,k) = kh_bg(i,k)
1661 kh_min_max_m(i,k) = kh_bg(i,k)
1662 endif ;
enddo ;
enddo
1672 do k=k_top+1,nz ;
do i=ish,ie ;
if (do_i(i))
then
1673 kh(i,k) = max(kh_bg(i,k), kh_detangle(i,k), &
1674 min(kh_min_m(i,k)*kh(i,k-1) + kh0_min_m(i,k), kh(i,k-1)))
1676 if (kh0_max_m(i,k) > kh_bg(i,k)) kh(i,k) = min(kh(i,k), kh0_max_m(i,k))
1677 if (kh0_max_p(i,k) > kh_bg(i,k)) kh(i,k) = min(kh(i,k), kh0_max_p(i,k))
1678 endif ;
enddo ;
enddo
1680 do k=nz,k_top+1,-1 ;
do i=ish,ie ;
if (do_i(i))
then
1681 kh(i,k) = max(kh(i,k), min(kh_min_p(i,k)*kh(i,k+1) + kh0_min_p(i,k), kh(i,k+1)))
1683 kh_max = max(kh_min_max_p(i,k), kh_max_p(i,k)*kh(i,k+1) + kh0_max_p(i,k))
1684 kh(i,k) = min(kh(i,k), kh_max)
1685 endif ;
enddo ;
enddo
1690 do k=k_top+1,nz ;
do i=ish,ie ;
if (do_i(i))
then
1691 kh_max = max(kh_min_max_m(i,k), kh_max_m(i,k)*kh(i,k-1) + kh0_max_m(i,k))
1692 if (kh(i,k) > kh_max) kh(i,k) = kh_max
1693 endif ;
enddo ;
enddo
1745 do k=k_top+1,nz ;
do i=ish,ie
1746 if (kh(i,k) > kh_u(i,j,k))
then
1747 dkh = (kh(i,k) - kh_u(i,j,k))
1748 int_slope_u(i,j,k) = dkh / kh(i,k)
1749 kh_u(i,j,k) = kh(i,k)
1753 do k=k_top+1,nz ;
do i=ish,ie
1754 if (kh(i,k) > kh_v(i,j,k))
then
1755 dkh = kh(i,k) - kh_v(i,j,k)
1756 int_slope_v(i,j,k) = dkh / kh(i,k)
1757 kh_v(i,j,k) = kh(i,k)
1769 type(time_type),
intent(in) :: time
1774 type(
diag_ctrl),
target,
intent(inout) :: diag
1779 #include "version_variable.h"
1780 character(len=40) :: mdl =
"MOM_thickness_diffuse"
1784 if (
associated(cs))
then
1786 "Thickness_diffuse_init called with an associated control structure.")
1788 else ;
allocate(cs) ;
endif
1794 call get_param(param_file, mdl,
"THICKNESSDIFFUSE", cs%thickness_diffuse, &
1795 "If true, interface heights are diffused with a "//&
1796 "coefficient of KHTH.", default=.false.)
1797 call get_param(param_file, mdl,
"KHTH", cs%Khth, &
1798 "The background horizontal thickness diffusivity.", &
1799 default=0.0, units=
"m2 s-1", scale=us%m_to_L**2*us%T_to_s)
1800 call get_param(param_file, mdl,
"KHTH_SLOPE_CFF", cs%KHTH_Slope_Cff, &
1801 "The nondimensional coefficient in the Visbeck formula "//&
1802 "for the interface depth diffusivity", units=
"nondim", &
1804 call get_param(param_file, mdl,
"KHTH_MIN", cs%KHTH_Min, &
1805 "The minimum horizontal thickness diffusivity.", &
1806 default=0.0, units=
"m2 s-1", scale=us%m_to_L**2*us%T_to_s)
1807 call get_param(param_file, mdl,
"KHTH_MAX", cs%KHTH_Max, &
1808 "The maximum horizontal thickness diffusivity.", &
1809 default=0.0, units=
"m2 s-1", scale=us%m_to_L**2*us%T_to_s)
1810 call get_param(param_file, mdl,
"KHTH_MAX_CFL", cs%max_Khth_CFL, &
1811 "The maximum value of the local diffusive CFL ratio that "//&
1812 "is permitted for the thickness diffusivity. 1.0 is the "//&
1813 "marginally unstable value in a pure layered model, but "//&
1814 "much smaller numbers (e.g. 0.1) seem to work better for "//&
1815 "ALE-based models.", units =
"nondimensional", default=0.8)
1821 if (cs%max_Khth_CFL < 0.0) cs%max_Khth_CFL = 0.0
1822 call get_param(param_file, mdl,
"DETANGLE_INTERFACES", cs%detangle_interfaces, &
1823 "If defined add 3-d structured enhanced interface height "//&
1824 "diffusivities to horizontally smooth jagged layers.", &
1826 cs%detangle_time = 0.0
1827 if (cs%detangle_interfaces) &
1828 call get_param(param_file, mdl,
"DETANGLE_TIMESCALE", cs%detangle_time, &
1829 "A timescale over which maximally jagged grid-scale "//&
1830 "thickness variations are suppressed. This must be "//&
1831 "longer than DT, or 0 to use DT.", units=
"s", default=0.0, scale=us%s_to_T)
1832 call get_param(param_file, mdl,
"KHTH_SLOPE_MAX", cs%slope_max, &
1833 "A slope beyond which the calculated isopycnal slope is "//&
1834 "not reliable and is scaled away.", units=
"nondim", default=0.01)
1835 call get_param(param_file, mdl,
"KD_SMOOTH", cs%kappa_smooth, &
1836 "A diapycnal diffusivity that is used to interpolate "//&
1837 "more sensible values of T & S into thin layers.", &
1838 default=1.0e-6, scale=us%m_to_Z**2*us%T_to_s)
1839 call get_param(param_file, mdl,
"KHTH_USE_FGNV_STREAMFUNCTION", cs%use_FGNV_streamfn, &
1840 "If true, use the streamfunction formulation of "//&
1841 "Ferrari et al., 2010, which effectively emphasizes "//&
1842 "graver vertical modes by smoothing in the vertical.", &
1844 call get_param(param_file, mdl,
"FGNV_FILTER_SCALE", cs%FGNV_scale, &
1845 "A coefficient scaling the vertical smoothing term in the "//&
1846 "Ferrari et al., 2010, streamfunction formulation.", &
1847 default=1., do_not_log=.not.cs%use_FGNV_streamfn)
1848 call get_param(param_file, mdl,
"FGNV_C_MIN", cs%FGNV_c_min, &
1849 "A minium wave speed used in the Ferrari et al., 2010, "//&
1850 "streamfunction formulation.", &
1851 default=0., units=
"m s-1", scale=us%m_s_to_L_T, do_not_log=.not.cs%use_FGNV_streamfn)
1852 call get_param(param_file, mdl,
"FGNV_STRAT_FLOOR", strat_floor, &
1853 "A floor for Brunt-Vasaila frequency in the Ferrari et al., 2010, "//&
1854 "streamfunction formulation, expressed as a fraction of planetary "//&
1855 "rotation, OMEGA. This should be tiny but non-zero to avoid degeneracy.", &
1856 default=1.e-15, units=
"nondim", do_not_log=.not.cs%use_FGNV_streamfn)
1857 call get_param(param_file, mdl,
"OMEGA", omega, &
1858 "The rotation rate of the earth.", &
1859 default=7.2921e-5, units=
"s-1", scale=us%T_to_s, do_not_log=.not.cs%use_FGNV_streamfn)
1860 if (cs%use_FGNV_streamfn) cs%N2_floor = (strat_floor*omega)**2
1861 call get_param(param_file, mdl,
"DEBUG", cs%debug, &
1862 "If true, write out verbose debugging data.", &
1863 default=.false., debuggingparam=.true.)
1865 call get_param(param_file, mdl,
"MEKE_GM_SRC_ALT", cs%GM_src_alt, &
1866 "If true, use the GM energy conversion form S^2*N^2*kappa rather \n"//&
1867 "than the streamfunction for the GM source term.", default=.false.)
1868 call get_param(param_file, mdl,
"MEKE_GEOMETRIC", cs%MEKE_GEOMETRIC, &
1869 "If true, uses the GM coefficient formulation \n"//&
1870 "from the GEOMETRIC framework (Marshall et al., 2012).", default=.false.)
1871 if (cs%MEKE_GEOMETRIC)
then
1873 call get_param(param_file, mdl,
"MEKE_GEOMETRIC_EPSILON", cs%MEKE_GEOMETRIC_epsilon, &
1874 "Minimum Eady growth rate used in the calculation of \n"//&
1875 "GEOMETRIC thickness diffusivity.", units=
"s-1", default=1.0e-7, scale=us%T_to_s)
1877 call get_param(param_file, mdl,
"MEKE_GEOMETRIC_ALPHA", cs%MEKE_GEOMETRIC_alpha, &
1878 "The nondimensional coefficient governing the efficiency of the GEOMETRIC \n"//&
1879 "thickness diffusion.", units=
"nondim", default=0.05)
1882 call get_param(param_file, mdl,
"USE_KH_IN_MEKE", cs%Use_KH_in_MEKE, &
1883 "If true, uses the thickness diffusivity calculated here to diffuse \n"//&
1884 "MEKE.", default=.false.)
1886 call get_param(param_file, mdl,
"USE_GME", cs%use_GME_thickness_diffuse, &
1887 "If true, use the GM+E backscatter scheme in association \n"//&
1888 "with the Gent and McWilliams parameterization.", default=.false.)
1890 if (cs%use_GME_thickness_diffuse)
then
1891 call safe_alloc_ptr(cs%KH_u_GME,g%IsdB,g%IedB,g%jsd,g%jed,g%ke+1)
1892 call safe_alloc_ptr(cs%KH_v_GME,g%isd,g%ied,g%JsdB,g%JedB,g%ke+1)
1895 cs%id_uhGM = register_diag_field(
'ocean_model',
'uhGM', diag%axesCuL, time, &
1896 'Time Mean Diffusive Zonal Thickness Flux', &
1897 'kg s-1', conversion=gv%H_to_kg_m2*us%L_to_m**2*us%s_to_T, &
1898 y_cell_method=
'sum', v_extensive=.true.)
1899 if (cs%id_uhGM > 0)
call safe_alloc_ptr(cdp%uhGM,g%IsdB,g%IedB,g%jsd,g%jed,g%ke)
1900 cs%id_vhGM = register_diag_field(
'ocean_model',
'vhGM', diag%axesCvL, time, &
1901 'Time Mean Diffusive Meridional Thickness Flux', &
1902 'kg s-1', conversion=gv%H_to_kg_m2*us%L_to_m**2*us%s_to_T, &
1903 x_cell_method=
'sum', v_extensive=.true.)
1904 if (cs%id_vhGM > 0)
call safe_alloc_ptr(cdp%vhGM,g%isd,g%ied,g%JsdB,g%JedB,g%ke)
1906 cs%id_GMwork = register_diag_field(
'ocean_model',
'GMwork', diag%axesT1, time, &
1907 'Integrated Tendency of Ocean Mesoscale Eddy KE from Parameterized Eddy Advection', &
1908 'W m-2', conversion=us%R_to_kg_m3*us%Z_to_m*us%L_to_m**2*us%s_to_T**3, cmor_field_name=
'tnkebto', &
1909 cmor_long_name=
'Integrated Tendency of Ocean Mesoscale Eddy KE from Parameterized Eddy Advection', &
1910 cmor_standard_name=
'tendency_of_ocean_eddy_kinetic_energy_content_due_to_parameterized_eddy_advection')
1911 if (cs%id_GMwork > 0)
call safe_alloc_ptr(cs%GMwork,g%isd,g%ied,g%jsd,g%jed)
1913 cs%id_KH_u = register_diag_field(
'ocean_model',
'KHTH_u', diag%axesCui, time, &
1914 'Parameterized mesoscale eddy advection diffusivity at U-point', &
1915 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
1916 cs%id_KH_v = register_diag_field(
'ocean_model',
'KHTH_v', diag%axesCvi, time, &
1917 'Parameterized mesoscale eddy advection diffusivity at V-point', &
1918 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
1919 cs%id_KH_t = register_diag_field(
'ocean_model',
'KHTH_t', diag%axesTL, time, &
1920 'Ocean Tracer Diffusivity due to Parameterized Mesoscale Advection', &
1921 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T, &
1922 cmor_field_name=
'diftrblo', &
1923 cmor_long_name=
'Ocean Tracer Diffusivity due to Parameterized Mesoscale Advection', &
1924 cmor_units=
'm2 s-1', &
1925 cmor_standard_name=
'ocean_tracer_diffusivity_due_to_parameterized_mesoscale_advection')
1927 cs%id_KH_u1 = register_diag_field(
'ocean_model',
'KHTH_u1', diag%axesCu1, time, &
1928 'Parameterized mesoscale eddy advection diffusivity at U-points (2-D)', &
1929 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
1930 cs%id_KH_v1 = register_diag_field(
'ocean_model',
'KHTH_v1', diag%axesCv1, time, &
1931 'Parameterized mesoscale eddy advection diffusivity at V-points (2-D)', &
1932 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
1933 cs%id_KH_t1 = register_diag_field(
'ocean_model',
'KHTH_t1', diag%axesT1, time, &
1934 'Parameterized mesoscale eddy advection diffusivity at T-points (2-D)', &
1935 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
1937 cs%id_slope_x = register_diag_field(
'ocean_model',
'neutral_slope_x', diag%axesCui, time, &
1938 'Zonal slope of neutral surface',
'nondim')
1939 if (cs%id_slope_x > 0)
call safe_alloc_ptr(cs%diagSlopeX,g%IsdB,g%IedB,g%jsd,g%jed,g%ke+1)
1940 cs%id_slope_y = register_diag_field(
'ocean_model',
'neutral_slope_y', diag%axesCvi, time, &
1941 'Meridional slope of neutral surface',
'nondim')
1942 if (cs%id_slope_y > 0)
call safe_alloc_ptr(cs%diagSlopeY,g%isd,g%ied,g%JsdB,g%JedB,g%ke+1)
1943 cs%id_sfn_x = register_diag_field(
'ocean_model',
'GM_sfn_x', diag%axesCui, time, &
1944 'Parameterized Zonal Overturning Streamfunction', &
1945 'm3 s-1', conversion=gv%H_to_m*us%L_to_m**2*us%s_to_T)
1946 cs%id_sfn_y = register_diag_field(
'ocean_model',
'GM_sfn_y', diag%axesCvi, time, &
1947 'Parameterized Meridional Overturning Streamfunction', &
1948 'm3 s-1', conversion=gv%H_to_m*us%L_to_m**2*us%s_to_T)
1949 cs%id_sfn_unlim_x = register_diag_field(
'ocean_model',
'GM_sfn_unlim_x', diag%axesCui, time, &
1950 'Parameterized Zonal Overturning Streamfunction before limiting/smoothing', &
1951 'm3 s-1', conversion=us%Z_to_m*us%L_to_m**2*us%s_to_T)
1952 cs%id_sfn_unlim_y = register_diag_field(
'ocean_model',
'GM_sfn_unlim_y', diag%axesCvi, time, &
1953 'Parameterized Meridional Overturning Streamfunction before limiting/smoothing', &
1954 'm3 s-1', conversion=us%Z_to_m*us%L_to_m**2*us%s_to_T)
1963 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)+1),
intent(inout) :: kh_u_gme
1965 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)+1),
intent(inout) :: kh_v_gme
1970 do k=1,g%ke+1 ;
do j = g%jsc, g%jec ;
do i = g%isc-1, g%iec
1971 kh_u_gme(i,j,k) = cs%KH_u_GME(i,j,k)
1972 enddo ;
enddo ;
enddo
1974 do k=1,g%ke+1 ;
do j = g%jsc-1, g%jec ;
do i = g%isc, g%iec
1975 kh_v_gme(i,j,k) = cs%KH_v_GME(i,j,k)
1976 enddo ;
enddo ;
enddo
1984 if (
associated(cs))
deallocate(cs)