23 implicit none ;
private
25 #include <MOM_memory.h>
27 public thickness_diffuse, thickness_diffuse_init, thickness_diffuse_end
29 public thickness_diffuse_get_kh
39 real :: khth_slope_cff
46 logical :: thickness_diffuse
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(szib_(g), szjb_(g))
143 real,
dimension(:,:),
pointer :: cg1 => null()
144 logical :: use_varmix, resoln_scaled, use_stored_slopes, khth_use_ebt_struct, use_visbeck
145 logical :: use_qg_leith
146 integer :: i, j, k, is, ie, js, je, nz
147 real :: hu(szi_(g), szj_(g))
148 real :: hv(szi_(g), szj_(g))
149 real :: kh_u_lay(szi_(g), szj_(g))
150 real :: kh_v_lay(szi_(g), szj_(g))
152 if (.not.
associated(cs))
call mom_error(fatal,
"MOM_thickness_diffuse:"// &
153 "Module must be initialized before it is used.")
155 if ((.not.cs%thickness_diffuse) .or. &
156 .not.( cs%Khth > 0.0 .or.
associated(varmix) .or.
associated(meke) ) )
return
158 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
159 h_neglect = gv%H_subroundoff
161 if (
associated(meke))
then
162 if (
associated(meke%GM_src))
then
163 do j=js,je ;
do i=is,ie ; meke%GM_src(i,j) = 0. ;
enddo ;
enddo
167 use_varmix = .false. ; resoln_scaled = .false. ; use_stored_slopes = .false.
168 khth_use_ebt_struct = .false. ; use_visbeck = .false. ; use_qg_leith = .false.
170 if (
associated(varmix))
then
171 use_varmix = varmix%use_variable_mixing .and. (cs%KHTH_Slope_Cff > 0.)
172 resoln_scaled = varmix%Resoln_scaled_KhTh
173 use_stored_slopes = varmix%use_stored_slopes
174 khth_use_ebt_struct = varmix%khth_use_ebt_struct
175 use_visbeck = varmix%use_Visbeck
176 use_qg_leith = varmix%use_QG_Leith_GM
177 if (
associated(varmix%cg1)) cg1 => varmix%cg1
184 do j=js,je ;
do i=is-1,ie
185 kh_u_cfl(i,j) = (0.25*cs%max_Khth_CFL) / &
186 (dt*(g%IdxCu(i,j)*g%IdxCu(i,j) + g%IdyCu(i,j)*g%IdyCu(i,j)))
189 do j=js-1,je ;
do i=is,ie
190 kh_v_cfl(i,j) = (0.25*cs%max_Khth_CFL) / &
191 (dt*(g%IdxCv(i,j)*g%IdxCv(i,j) + g%IdyCv(i,j)*g%IdyCv(i,j)))
195 call find_eta(h, tv, g, gv, us, e, halo_size=1)
203 do j=js,je;
do i=is-1,ie
204 khth_loc_u(i,j) = cs%Khth
209 if (use_visbeck)
then
210 do j=js,je ;
do i=is-1,ie
211 khth_loc_u(i,j) = khth_loc_u(i,j) + cs%KHTH_Slope_Cff*varmix%L2u(i,j)*varmix%SN_u(i,j)
216 if (
associated(meke))
then ;
if (
associated(meke%Kh))
then
218 if (cs%MEKE_GEOMETRIC)
then
219 do j=js,je ;
do i=is-1,ie
220 khth_loc_u(i,j) = khth_loc_u(i,j) + &
221 g%mask2dCu(i,j) * cs%MEKE_GEOMETRIC_alpha * 0.5*(meke%MEKE(i,j)+meke%MEKE(i+1,j)) / &
222 (varmix%SN_u(i,j) + cs%MEKE_GEOMETRIC_epsilon)
225 do j=js,je ;
do i=is-1,ie
226 khth_loc_u(i,j) = khth_loc_u(i,j) + meke%KhTh_fac*sqrt(meke%Kh(i,j)*meke%Kh(i+1,j))
231 if (resoln_scaled)
then
233 do j=js,je;
do i=is-1,ie
234 khth_loc_u(i,j) = khth_loc_u(i,j) * varmix%Res_fn_u(i,j)
238 if (cs%Khth_Max > 0)
then
240 do j=js,je;
do i=is-1,ie
241 khth_loc_u(i,j) = max(cs%Khth_min, min(khth_loc_u(i,j),cs%Khth_Max))
245 do j=js,je;
do i=is-1,ie
246 khth_loc_u(i,j) = max(cs%Khth_min, khth_loc_u(i,j))
250 do j=js,je;
do i=is-1,ie
251 kh_u(i,j,1) = min(kh_u_cfl(i,j), khth_loc_u(i,j))
254 if (khth_use_ebt_struct)
then
256 do k=2,nz+1 ;
do j=js,je ;
do i=is-1,ie
257 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) )
258 enddo ;
enddo ;
enddo
261 do k=2,nz+1 ;
do j=js,je ;
do i=is-1,ie
262 kh_u(i,j,k) = kh_u(i,j,1)
263 enddo ;
enddo ;
enddo
268 if (use_qg_leith)
then
269 do k=1,nz ;
do j=js,je ;
do i=is-1,ie
270 kh_u(i,j,k) = varmix%KH_u_QG(i,j,k)
271 enddo ;
enddo ;
enddo
276 if (cs%use_GME_thickness_diffuse)
then
277 do k=1,nz+1 ;
do j=js,je ;
do i=is-1,ie
278 cs%KH_u_GME(i,j,k) = kh_u(i,j,k)
279 enddo ;
enddo ;
enddo
283 do j=js-1,je ;
do i=is,ie
284 khth_loc(i,j) = cs%Khth
289 if (use_visbeck)
then
290 do j=js-1,je ;
do i=is,ie
291 khth_loc(i,j) = khth_loc(i,j) + cs%KHTH_Slope_Cff*varmix%L2v(i,j)*varmix%SN_v(i,j)
295 if (
associated(meke))
then ;
if (
associated(meke%Kh))
then
297 if (cs%MEKE_GEOMETRIC)
then
298 do j=js-1,je ;
do i=is,ie
299 khth_loc(i,j) = khth_loc(i,j) + &
300 g%mask2dCv(i,j) * cs%MEKE_GEOMETRIC_alpha * 0.5*(meke%MEKE(i,j)+meke%MEKE(i,j+1)) / &
301 (varmix%SN_v(i,j) + cs%MEKE_GEOMETRIC_epsilon)
304 do j=js-1,je ;
do i=is,ie
305 khth_loc(i,j) = khth_loc(i,j) + meke%KhTh_fac*sqrt(meke%Kh(i,j)*meke%Kh(i,j+1))
310 if (resoln_scaled)
then
312 do j=js-1,je ;
do i=is,ie
313 khth_loc(i,j) = khth_loc(i,j) * varmix%Res_fn_v(i,j)
317 if (cs%Khth_Max > 0)
then
319 do j=js-1,je ;
do i=is,ie
320 khth_loc(i,j) = max(cs%Khth_min, min(khth_loc(i,j),cs%Khth_Max))
324 do j=js-1,je ;
do i=is,ie
325 khth_loc(i,j) = max(cs%Khth_min, khth_loc(i,j))
329 if (cs%max_Khth_CFL > 0.0)
then
331 do j=js-1,je ;
do i=is,ie
332 kh_v(i,j,1) = min(kh_v_cfl(i,j), khth_loc(i,j))
336 if (khth_use_ebt_struct)
then
338 do k=2,nz+1 ;
do j=js-1,je ;
do i=is,ie
339 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) )
340 enddo ;
enddo ;
enddo
343 do k=2,nz+1 ;
do j=js-1,je ;
do i=is,ie
344 kh_v(i,j,k) = kh_v(i,j,1)
345 enddo ;
enddo ;
enddo
350 if (use_qg_leith)
then
351 do k=1,nz ;
do j=js-1,je ;
do i=is,ie
352 kh_v(i,j,k) = varmix%KH_v_QG(i,j,k)
353 enddo ;
enddo ;
enddo
358 if (cs%use_GME_thickness_diffuse)
then
359 do k=1,nz+1 ;
do j=js-1,je ;
do i=is,ie
360 cs%KH_v_GME(i,j,k) = kh_v(i,j,k)
361 enddo ;
enddo ;
enddo
364 if (
associated(meke))
then ;
if (
associated(meke%Kh))
then
366 if (cs%MEKE_GEOMETRIC)
then
367 do j=js,je ;
do i=is,ie
368 meke%Kh(i,j) = cs%MEKE_GEOMETRIC_alpha * meke%MEKE(i,j) / &
369 (0.25*(varmix%SN_u(i,j)+varmix%SN_u(i-1,j)+varmix%SN_v(i,j)+varmix%SN_v(i,j-1)) + &
370 cs%MEKE_GEOMETRIC_epsilon)
377 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
379 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
382 if (cs%detangle_interfaces)
then
383 call add_detangling_kh(h, e, kh_u, kh_v, kh_u_cfl, kh_v_cfl, tv, dt, g, gv, us, &
384 cs, int_slope_u, int_slope_v)
388 call uvchksum(
"Kh_[uv]", kh_u, kh_v, g%HI,haloshift=0)
389 call uvchksum(
"int_slope_[uv]", int_slope_u, int_slope_v, g%HI, haloshift=0)
390 call hchksum(h,
"thickness_diffuse_1 h", g%HI, haloshift=1, scale=gv%H_to_m)
391 call hchksum(e,
"thickness_diffuse_1 e", g%HI, haloshift=1, scale=us%Z_to_m)
392 if (use_stored_slopes)
then
393 call uvchksum(
"VarMix%slope_[xy]", varmix%slope_x, varmix%slope_y, &
396 if (
associated(tv%eqn_of_state))
then
397 call hchksum(tv%T,
"thickness_diffuse T", g%HI, haloshift=1)
398 call hchksum(tv%S,
"thickness_diffuse S", g%HI, haloshift=1)
403 if (use_stored_slopes)
then
404 call thickness_diffuse_full(h, e, kh_u, kh_v, tv, uhd, vhd, cg1, dt, g, gv, us, meke, cs, &
405 int_slope_u, int_slope_v, varmix%slope_x, varmix%slope_y)
407 call thickness_diffuse_full(h, e, kh_u, kh_v, tv, uhd, vhd, cg1, dt, g, gv, us, meke, cs, &
408 int_slope_u, int_slope_v)
411 if (
associated(meke) .AND.
associated(varmix))
then
412 if (
associated(meke%Rd_dx_h) .and.
associated(varmix%Rd_dx_h))
then
414 do j=js,je ;
do i=is,ie
415 meke%Rd_dx_h(i,j) = varmix%Rd_dx_h(i,j)
421 if (query_averaging_enabled(cs%diag))
then
422 if (cs%id_uhGM > 0)
call post_data(cs%id_uhGM, uhd, cs%diag)
423 if (cs%id_vhGM > 0)
call post_data(cs%id_vhGM, vhd, cs%diag)
424 if (cs%id_GMwork > 0)
call post_data(cs%id_GMwork, cs%GMwork, cs%diag)
425 if (cs%id_KH_u > 0)
call post_data(cs%id_KH_u, kh_u, cs%diag)
426 if (cs%id_KH_v > 0)
call post_data(cs%id_KH_v, kh_v, cs%diag)
427 if (cs%id_KH_u1 > 0)
call post_data(cs%id_KH_u1, kh_u(:,:,1), cs%diag)
428 if (cs%id_KH_v1 > 0)
call post_data(cs%id_KH_v1, kh_v(:,:,1), cs%diag)
435 if (cs%id_KH_t > 0 .or. cs%id_KH_t1 > 0 .or. cs%Use_KH_in_MEKE)
then
439 do j=js,je ;
do i=is-1,ie
440 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)
441 if (hu(i,j) /= 0.0) hu(i,j) = 1.0
442 kh_u_lay(i,j) = 0.5*(kh_u(i,j,k)+kh_u(i,j,k+1))
444 do j=js-1,je ;
do i=is,ie
445 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)
446 if (hv(i,j) /= 0.0) hv(i,j) = 1.0
447 kh_v_lay(i,j) = 0.5*(kh_v(i,j,k)+kh_v(i,j,k+1))
450 do j=js,je ;
do i=is,ie
451 kh_t(i,j,k) = ((hu(i-1,j)*kh_u_lay(i-1,j)+hu(i,j)*kh_u_lay(i,j)) &
452 +(hv(i,j-1)*kh_v_lay(i,j-1)+hv(i,j)*kh_v_lay(i,j))) &
453 / (hu(i-1,j)+hu(i,j)+hv(i,j-1)+hv(i,j)+h_neglect)
457 if (cs%Use_KH_in_MEKE)
then
458 meke%Kh_diff(:,:) = 0.0
460 do j=js,je ;
do i=is,ie
461 meke%Kh_diff(i,j) = meke%Kh_diff(i,j) + kh_t(i,j,k) * h(i,j,k)
465 do j=js,je ;
do i=is,ie
466 meke%Kh_diff(i,j) = meke%Kh_diff(i,j) / max(1.0,g%bathyT(i,j))
470 if (cs%id_KH_t > 0)
call post_data(cs%id_KH_t, kh_t, cs%diag)
471 if (cs%id_KH_t1 > 0)
call post_data(cs%id_KH_t1, kh_t(:,:,1), cs%diag)
478 do j=js,je ;
do i=is-1,ie
479 uhtr(i,j,k) = uhtr(i,j,k) + uhd(i,j,k)*dt
480 if (
associated(cdp%uhGM)) cdp%uhGM(i,j,k) = uhd(i,j,k)
482 do j=js-1,je ;
do i=is,ie
483 vhtr(i,j,k) = vhtr(i,j,k) + vhd(i,j,k)*dt
484 if (
associated(cdp%vhGM)) cdp%vhGM(i,j,k) = vhd(i,j,k)
486 do j=js,je ;
do i=is,ie
487 h(i,j,k) = h(i,j,k) - dt * g%IareaT(i,j) * &
488 ((uhd(i,j,k) - uhd(i-1,j,k)) + (vhd(i,j,k) - vhd(i,j-1,k)))
489 if (h(i,j,k) < gv%Angstrom_H) h(i,j,k) = gv%Angstrom_H
496 call diag_update_remap_grids(cs%diag)
499 call uvchksum(
"thickness_diffuse [uv]hD", uhd, vhd, &
500 g%HI, haloshift=0, scale=gv%H_to_m)
501 call uvchksum(
"thickness_diffuse [uv]htr", uhtr, vhtr, &
502 g%HI, haloshift=0, scale=gv%H_to_m)
503 call hchksum(h,
"thickness_diffuse h", g%HI, haloshift=0, scale=gv%H_to_m)
506 end subroutine thickness_diffuse
511 subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV, US, MEKE, &
512 CS, int_slope_u, int_slope_v, slope_x, slope_y)
516 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
517 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1),
intent(in) :: e
518 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)+1),
intent(in) :: Kh_u
520 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)+1),
intent(in) :: Kh_v
523 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(out) :: uhD
525 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(out) :: vhD
527 real,
dimension(:,:),
pointer :: cg1
528 real,
intent(in) :: dt
531 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)+1),
optional,
intent(in) :: int_slope_u
535 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)+1),
optional,
intent(in) :: int_slope_v
539 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)+1),
optional,
intent(in) :: slope_x
540 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)+1),
optional,
intent(in) :: slope_y
542 real,
dimension(SZI_(G), SZJ_(G), SZK_(G)) :: &
553 real,
dimension(SZI_(G), SZJB_(G), SZK_(G)+1) :: &
557 real,
dimension(SZIB_(G), SZJ_(G), SZK_(G)+1) :: &
561 real,
dimension(SZI_(G), SZJ_(G), SZK_(G)+1) :: &
564 real,
dimension(SZIB_(G)) :: &
567 real,
dimension(SZI_(G)) :: &
570 real :: uhtot(SZIB_(G), SZJ_(G))
571 real :: vhtot(SZI_(G), SZJB_(G))
572 real,
dimension(SZIB_(G)) :: &
576 real,
dimension(SZI_(G)) :: &
580 real :: Work_u(SZIB_(G), SZJ_(G))
581 real :: Work_v(SZI_(G), SZJB_(G))
590 real :: drdi_u(SZIB_(G), SZK_(G)+1)
591 real :: drdj_v(SZI_(G), SZK_(G)+1)
592 real :: drdkDe_u(SZIB_(G),SZK_(G)+1)
594 real :: drdkDe_v(SZI_(G),SZK_(G)+1)
596 real :: hg2A, hg2B, hg2L, hg2R
597 real :: haA, haB, haL, haR
599 real :: wtA, wtB, wtL, wtR
603 real :: c2_h_u(SZIB_(G), SZK_(G)+1)
604 real :: c2_h_v(SZI_(G), SZK_(G)+1)
605 real :: hN2_u(SZIB_(G), SZK_(G)+1)
606 real :: hN2_v(SZI_(G), SZK_(G)+1)
609 real :: Sfn_unlim_u(SZIB_(G), SZK_(G)+1)
610 real :: Sfn_unlim_v(SZI_(G), SZK_(G)+1)
611 real :: slope2_Ratio_u(SZIB_(G), SZK_(G)+1)
612 real :: slope2_Ratio_v(SZI_(G), SZK_(G)+1)
635 real,
dimension(SZIB_(G), SZJ_(G), SZK_(G)+1) :: diag_sfn_x, diag_sfn_unlim_x
636 real,
dimension(SZI_(G), SZJB_(G), SZK_(G)+1) :: diag_sfn_y, diag_sfn_unlim_y
637 logical :: present_int_slope_u, present_int_slope_v
638 logical :: present_slope_x, present_slope_y, calc_derivatives
639 integer :: is, ie, js, je, nz, IsdB
641 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke ; isdb = g%IsdB
644 i_slope_max2 = 1.0 / (cs%slope_max**2)
645 g_scale = gv%g_Earth*us%L_to_m**2*us%s_to_T**2 * gv%H_to_m
646 h_neglect = gv%H_subroundoff ; h_neglect2 = h_neglect**2
647 dz_neglect = gv%H_subroundoff*gv%H_to_Z
648 g_rho0 = gv%g_Earth*us%L_to_m**2*us%s_to_T**2 / gv%Rho0
649 n2_floor = cs%N2_floor*us%Z_to_m**2
651 use_eos =
associated(tv%eqn_of_state)
652 present_int_slope_u =
PRESENT(int_slope_u)
653 present_int_slope_v =
PRESENT(int_slope_v)
654 present_slope_x =
PRESENT(slope_x)
655 present_slope_y =
PRESENT(slope_y)
657 nk_linear = max(gv%nkml, 1)
659 slope_x_pe(:,:,:) = 0.0
660 slope_y_pe(:,:,:) = 0.0
661 hn2_x_pe(:,:,:) = 0.0
662 hn2_y_pe(:,:,:) = 0.0
665 if (
associated(meke)) find_work =
associated(meke%GM_src)
666 find_work = (
associated(cs%GMwork) .or. find_work)
669 call vert_fill_ts(h, tv%T, tv%S, cs%kappa_smooth*dt, t, s, g, gv, 1, larger_h_denom=.true.)
672 if (cs%use_FGNV_streamfn .and. .not.
associated(cg1))
call mom_error(fatal, &
673 "cg1 must be associated when using FGNV streamfunction.")
680 do j=js-1,je+1 ;
do i=is-1,ie+1
681 h_avail_rsum(i,j,1) = 0.0
684 h_avail(i,j,1) = max(i4dt*g%areaT(i,j)*(h(i,j,1)-gv%Angstrom_H),0.0)
685 h_avail_rsum(i,j,2) = h_avail(i,j,1)
687 pres(i,j,2) = pres(i,j,1) + gv%H_to_Pa*h(i,j,1)
691 do k=2,nz ;
do i=is-1,ie+1
692 h_avail(i,j,k) = max(i4dt*g%areaT(i,j)*(h(i,j,k)-gv%Angstrom_H),0.0)
693 h_avail_rsum(i,j,k+1) = h_avail_rsum(i,j,k) + h_avail(i,j,k)
694 h_frac(i,j,k) = 0.0 ;
if (h_avail(i,j,k) > 0.0) &
695 h_frac(i,j,k) = h_avail(i,j,k) / h_avail_rsum(i,j,k+1)
696 pres(i,j,k+1) = pres(i,j,k) + gv%H_to_Pa*h(i,j,k)
700 do j=js,je ;
do i=is-1,ie
701 uhtot(i,j) = 0.0 ; work_u(i,j) = 0.0
702 diag_sfn_x(i,j,1) = 0.0 ; diag_sfn_unlim_x(i,j,1) = 0.0
703 diag_sfn_x(i,j,nz+1) = 0.0 ; diag_sfn_unlim_x(i,j,nz+1) = 0.0
706 do j=js-1,je ;
do i=is,ie
707 vhtot(i,j) = 0.0 ; work_v(i,j) = 0.0
708 diag_sfn_y(i,j,1) = 0.0 ; diag_sfn_unlim_y(i,j,1) = 0.0
709 diag_sfn_y(i,j,nz+1) = 0.0 ; diag_sfn_unlim_y(i,j,nz+1) = 0.0
727 do i=is-1,ie ; hn2_u(i,1) = 0. ; hn2_u(i,nz+1) = 0. ;
enddo
729 if (find_work .and. .not.(use_eos))
then
730 drdia = 0.0 ; drdib = 0.0
731 drdkl = gv%Rlay(k)-gv%Rlay(k-1) ; drdkr = drdkl
734 calc_derivatives = use_eos .and. (k >= nk_linear) .and. &
735 (find_work .or. .not. present_slope_x .or. cs%use_FGNV_streamfn)
738 if (calc_derivatives)
then
740 pres_u(i) = 0.5*(pres(i,j,k) + pres(i+1,j,k))
741 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)))
742 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)))
745 drho_ds_u, (is-isdb+1)-1, ie-is+2, tv%eqn_of_state)
749 if (calc_derivatives)
then
751 drdia = drho_dt_u(i) * (t(i+1,j,k-1)-t(i,j,k-1)) + &
752 drho_ds_u(i) * (s(i+1,j,k-1)-s(i,j,k-1))
753 drdib = drho_dt_u(i) * (t(i+1,j,k)-t(i,j,k)) + &
754 drho_ds_u(i) * (s(i+1,j,k)-s(i,j,k))
757 drdkl = (drho_dt_u(i) * (t(i,j,k)-t(i,j,k-1)) + &
758 drho_ds_u(i) * (s(i,j,k)-s(i,j,k-1)))
759 drdkr = (drho_dt_u(i) * (t(i+1,j,k)-t(i+1,j,k-1)) + &
760 drho_ds_u(i) * (s(i+1,j,k)-s(i+1,j,k-1)))
761 drdkde_u(i,k) = drdkr * e(i+1,j,k) - drdkl * e(i,j,k)
762 elseif (find_work)
then
763 drdkde_u(i,k) = drdkr * e(i+1,j,k) - drdkl * e(i,j,k)
766 if (find_work) drdi_u(i,k) = drdib
768 if (k > nk_linear)
then
770 if (cs%use_FGNV_streamfn .or. find_work .or. .not.present_slope_x)
then
771 hg2l = h(i,j,k-1)*h(i,j,k) + h_neglect2
772 hg2r = h(i+1,j,k-1)*h(i+1,j,k) + h_neglect2
773 hal = 0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect
774 har = 0.5*(h(i+1,j,k-1) + h(i+1,j,k)) + h_neglect
775 if (gv%Boussinesq)
then
776 dzal = hal * gv%H_to_Z ; dzar = har * gv%H_to_Z
778 dzal = 0.5*(e(i,j,k-1) - e(i,j,k+1)) + dz_neglect
779 dzar = 0.5*(e(i+1,j,k-1) - e(i+1,j,k+1)) + dz_neglect
783 wtl = hg2l*(har*dzar) ; wtr = hg2r*(hal*dzal)
785 drdz = (wtl * drdkl + wtr * drdkr) / (dzal*wtl + dzar*wtr)
789 hg2a = h(i,j,k-1)*h(i+1,j,k-1) + h_neglect2
790 hg2b = h(i,j,k)*h(i+1,j,k) + h_neglect2
791 haa = 0.5*(h(i,j,k-1) + h(i+1,j,k-1)) + h_neglect
792 hab = 0.5*(h(i,j,k) + h(i+1,j,k)) + h_neglect
795 hn2_u(i,k) = (0.5 * gv%H_to_Z * ( hg2a / haa + hg2b / hab )) * &
796 max(drdz*g_rho0, n2_floor)
798 if (present_slope_x)
then
799 slope = slope_x(i,j,k)
800 slope2_ratio_u(i,k) = slope**2 * i_slope_max2
804 wta = hg2a*hab ; wtb = hg2b*haa
806 drdx = ((wta * drdia + wtb * drdib) / (wta + wtb) - &
807 drdz * (e(i,j,k)-e(i+1,j,k))) * g%IdxCu(i,j)
811 mag_grad2 = drdx**2 + (us%m_to_Z*drdz)**2
812 if (mag_grad2 > 0.0)
then
813 slope = drdx / sqrt(mag_grad2)
814 slope2_ratio_u(i,k) = slope**2 * i_slope_max2
817 slope2_ratio_u(i,k) = 1.0e20
823 if (present_int_slope_u)
then
824 slope = (1.0 - int_slope_u(i,j,k)) * slope + &
825 int_slope_u(i,j,k) * us%Z_to_m*((e(i+1,j,k)-e(i,j,k)) * g%IdxCu(i,j))
826 slope2_ratio_u(i,k) = (1.0 - int_slope_u(i,j,k)) * slope2_ratio_u(i,k)
829 slope_x_pe(i,j,k) = min(slope,cs%slope_max)
830 hn2_x_pe(i,j,k) = hn2_u(i,k) * us%m_to_Z
831 if (cs%id_slope_x > 0) cs%diagSlopeX(i,j,k) = slope
834 sfn_unlim_u(i,k) = -((kh_u(i,j,k)*g%dy_Cu(i,j))*us%m_to_Z*slope)
838 if (sfn_unlim_u(i,k) > 0.0)
then
839 if (e(i,j,k) < e(i+1,j,nz+1))
then
840 sfn_unlim_u(i,k) = 0.0
842 elseif (e(i+1,j,nz+1) > e(i,j,k+1))
then
845 sfn_unlim_u(i,k) = sfn_unlim_u(i,k) * ((e(i,j,k) - e(i+1,j,nz+1)) / &
846 ((e(i,j,k) - e(i,j,k+1)) + dz_neglect))
849 if (e(i+1,j,k) < e(i,j,nz+1))
then ; sfn_unlim_u(i,k) = 0.0
850 elseif (e(i,j,nz+1) > e(i+1,j,k+1))
then
851 sfn_unlim_u(i,k) = sfn_unlim_u(i,k) * ((e(i+1,j,k) - e(i,j,nz+1)) / &
852 ((e(i+1,j,k) - e(i+1,j,k+1)) + dz_neglect))
857 if (present_slope_x)
then
858 slope = slope_x(i,j,k)
860 slope = us%Z_to_m*((e(i,j,k)-e(i+1,j,k))*g%IdxCu(i,j)) * g%mask2dCu(i,j)
862 if (cs%id_slope_x > 0) cs%diagSlopeX(i,j,k) = slope
863 sfn_unlim_u(i,k) = ((kh_u(i,j,k)*g%dy_Cu(i,j))*us%m_to_Z*slope)
864 hn2_u(i,k) = us%L_to_m**2*us%s_to_T**2*gv%g_prime(k)
867 hn2_u(i,k) = n2_floor * dz_neglect
868 sfn_unlim_u(i,k) = 0.
870 if (cs%id_sfn_unlim_x>0) diag_sfn_unlim_x(i,j,k) = sfn_unlim_u(i,k)
874 if (cs%use_FGNV_streamfn)
then
875 do k=1,nz ;
do i=is-1,ie ;
if (g%mask2dCu(i,j)>0.)
then
876 h_harm = max( h_neglect, &
877 2. * h(i,j,k) * h(i+1,j,k) / ( ( h(i,j,k) + h(i+1,j,k) ) + h_neglect ) )
878 c2_h_u(i,k) = cs%FGNV_scale * ( 0.5*( cg1(i,j) + cg1(i+1,j) ) )**2 / (gv%H_to_Z*h_harm)
879 endif ;
enddo ;
enddo
883 if (g%mask2dCu(i,j)>0.)
then
884 sfn_unlim_u(i,:) = ( 1. + cs%FGNV_scale ) * sfn_unlim_u(i,:)
885 call streamfn_solver(nz, c2_h_u(i,:), hn2_u(i,:), sfn_unlim_u(i,:))
887 sfn_unlim_u(i,:) = 0.
894 if (k > nk_linear)
then
897 if (uhtot(i,j) <= 0.0)
then
899 sfn_safe = uhtot(i,j) * (1.0 - h_frac(i,j,k)) * gv%H_to_Z
901 sfn_safe = uhtot(i,j) * (1.0 - h_frac(i+1,j,k)) * gv%H_to_Z
905 sfn_est = (sfn_unlim_u(i,k) + slope2_ratio_u(i,k)*sfn_safe) / (1.0 + slope2_ratio_u(i,k))
907 sfn_est = sfn_unlim_u(i,k)
912 sfn_in_h = min(max(sfn_est * gv%Z_to_H, -h_avail_rsum(i,j,k)), h_avail_rsum(i+1,j,k))
916 uhd(i,j,k) = max(min((sfn_in_h - uhtot(i,j)), h_avail(i,j,k)), &
919 if (cs%id_sfn_x>0) diag_sfn_x(i,j,k) = diag_sfn_x(i,j,k+1) + uhd(i,j,k)
931 if (uhtot(i,j) <= 0.0)
then
932 uhd(i,j,k) = -uhtot(i,j) * h_frac(i,j,k)
934 uhd(i,j,k) = -uhtot(i,j) * h_frac(i+1,j,k)
937 if (cs%id_sfn_x>0) diag_sfn_x(i,j,k) = diag_sfn_x(i,j,k+1) + uhd(i,j,k)
945 uhtot(i,j) = uhtot(i,j) + uhd(i,j,k)
954 work_u(i,j) = work_u(i,j) + g_scale * &
955 ( uhtot(i,j) * drdkde_u(i,k) - &
956 (uhd(i,j,k) * drdi_u(i,k)) * 0.25 * &
957 ((e(i,j,k) + e(i,j,k+1)) + (e(i+1,j,k) + e(i+1,j,k+1))) )
980 if (find_work .and. .not.(use_eos))
then
981 drdja = 0.0 ; drdjb = 0.0
982 drdkl = gv%Rlay(k)-gv%Rlay(k-1) ; drdkr = drdkl
985 calc_derivatives = use_eos .and. (k >= nk_linear) .and. &
986 (find_work .or. .not. present_slope_y)
988 if (calc_derivatives)
then
990 pres_v(i) = 0.5*(pres(i,j,k) + pres(i,j+1,k))
991 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)))
992 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)))
995 drho_ds_v, is, ie-is+1, tv%eqn_of_state)
998 if (calc_derivatives)
then
1000 drdja = drho_dt_v(i) * (t(i,j+1,k-1)-t(i,j,k-1)) + &
1001 drho_ds_v(i) * (s(i,j+1,k-1)-s(i,j,k-1))
1002 drdjb = drho_dt_v(i) * (t(i,j+1,k)-t(i,j,k)) + &
1003 drho_ds_v(i) * (s(i,j+1,k)-s(i,j,k))
1006 drdkl = (drho_dt_v(i) * (t(i,j,k)-t(i,j,k-1)) + &
1007 drho_ds_v(i) * (s(i,j,k)-s(i,j,k-1)))
1008 drdkr = (drho_dt_v(i) * (t(i,j+1,k)-t(i,j+1,k-1)) + &
1009 drho_ds_v(i) * (s(i,j+1,k)-s(i,j+1,k-1)))
1010 drdkde_v(i,k) = drdkr * e(i,j+1,k) - drdkl * e(i,j,k)
1011 elseif (find_work)
then
1012 drdkde_v(i,k) = drdkr * e(i,j+1,k) - drdkl * e(i,j,k)
1015 if (find_work) drdj_v(i,k) = drdjb
1017 if (k > nk_linear)
then
1019 if (cs%use_FGNV_streamfn .or. find_work .or. .not. present_slope_y)
then
1020 hg2l = h(i,j,k-1)*h(i,j,k) + h_neglect2
1021 hg2r = h(i,j+1,k-1)*h(i,j+1,k) + h_neglect2
1022 hal = 0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect
1023 har = 0.5*(h(i,j+1,k-1) + h(i,j+1,k)) + h_neglect
1024 if (gv%Boussinesq)
then
1025 dzal = hal * gv%H_to_Z ; dzar = har * gv%H_to_Z
1027 dzal = 0.5*(e(i,j,k-1) - e(i,j,k+1)) + dz_neglect
1028 dzar = 0.5*(e(i,j+1,k-1) - e(i,j+1,k+1)) + dz_neglect
1032 wtl = hg2l*(har*dzar) ; wtr = hg2r*(hal*dzal)
1034 drdz = (wtl * drdkl + wtr * drdkr) / (dzal*wtl + dzar*wtr)
1038 hg2a = h(i,j,k-1)*h(i,j+1,k-1) + h_neglect2
1039 hg2b = h(i,j,k)*h(i,j+1,k) + h_neglect2
1040 haa = 0.5*(h(i,j,k-1) + h(i,j+1,k-1)) + h_neglect
1041 hab = 0.5*(h(i,j,k) + h(i,j+1,k)) + h_neglect
1044 hn2_v(i,k) = (0.5 * gv%H_to_Z * ( hg2a / haa + hg2b / hab )) * &
1045 max(drdz*g_rho0, n2_floor)
1047 if (present_slope_y)
then
1048 slope = slope_y(i,j,k)
1049 slope2_ratio_v(i,k) = slope**2 * i_slope_max2
1053 wta = hg2a*hab ; wtb = hg2b*haa
1055 drdy = ((wta * drdja + wtb * drdjb) / (wta + wtb) - &
1056 drdz * (e(i,j,k)-e(i,j+1,k))) * g%IdyCv(i,j)
1060 mag_grad2 = drdy**2 + (us%m_to_Z*drdz)**2
1061 if (mag_grad2 > 0.0)
then
1062 slope = drdy / sqrt(mag_grad2)
1063 slope2_ratio_v(i,k) = slope**2 * i_slope_max2
1066 slope2_ratio_v(i,k) = 1.0e20
1072 if (present_int_slope_v)
then
1073 slope = (1.0 - int_slope_v(i,j,k)) * slope + &
1074 int_slope_v(i,j,k) * us%Z_to_m*((e(i,j+1,k)-e(i,j,k)) * g%IdyCv(i,j))
1075 slope2_ratio_v(i,k) = (1.0 - int_slope_v(i,j,k)) * slope2_ratio_v(i,k)
1078 slope_y_pe(i,j,k) = min(slope,cs%slope_max)
1079 hn2_y_pe(i,j,k) = hn2_v(i,k) * us%m_to_Z
1080 if (cs%id_slope_y > 0) cs%diagSlopeY(i,j,k) = slope
1083 sfn_unlim_v(i,k) = -((kh_v(i,j,k)*g%dx_Cv(i,j))*us%m_to_Z*slope)
1087 if (sfn_unlim_v(i,k) > 0.0)
then
1088 if (e(i,j,k) < e(i,j+1,nz+1))
then
1089 sfn_unlim_v(i,k) = 0.0
1091 elseif (e(i,j+1,nz+1) > e(i,j,k+1))
then
1094 sfn_unlim_v(i,k) = sfn_unlim_v(i,k) * ((e(i,j,k) - e(i,j+1,nz+1)) / &
1095 ((e(i,j,k) - e(i,j,k+1)) + dz_neglect))
1098 if (e(i,j+1,k) < e(i,j,nz+1))
then ; sfn_unlim_v(i,k) = 0.0
1099 elseif (e(i,j,nz+1) > e(i,j+1,k+1))
then
1100 sfn_unlim_v(i,k) = sfn_unlim_v(i,k) * ((e(i,j+1,k) - e(i,j,nz+1)) / &
1101 ((e(i,j+1,k) - e(i,j+1,k+1)) + dz_neglect))
1106 if (present_slope_y)
then
1107 slope = slope_y(i,j,k)
1109 slope = us%Z_to_m*((e(i,j,k)-e(i,j+1,k))*g%IdyCv(i,j)) * g%mask2dCv(i,j)
1111 if (cs%id_slope_y > 0) cs%diagSlopeY(i,j,k) = slope
1112 sfn_unlim_v(i,k) = ((kh_v(i,j,k)*g%dx_Cv(i,j))*us%m_to_Z*slope)
1113 hn2_v(i,k) = us%L_to_m**2*us%s_to_T**2*gv%g_prime(k)
1116 hn2_v(i,k) = n2_floor * dz_neglect
1117 sfn_unlim_v(i,k) = 0.
1119 if (cs%id_sfn_unlim_y>0) diag_sfn_unlim_y(i,j,k) = sfn_unlim_v(i,k)
1123 if (cs%use_FGNV_streamfn)
then
1124 do k=1,nz ;
do i=is,ie ;
if (g%mask2dCv(i,j)>0.)
then
1125 h_harm = max( h_neglect, &
1126 2. * h(i,j,k) * h(i,j+1,k) / ( ( h(i,j,k) + h(i,j+1,k) ) + h_neglect ) )
1127 c2_h_v(i,k) = cs%FGNV_scale * ( 0.5*( cg1(i,j) + cg1(i,j+1) ) )**2 / (gv%H_to_Z*h_harm)
1128 endif ;
enddo ;
enddo
1132 if (g%mask2dCv(i,j)>0.)
then
1133 sfn_unlim_v(i,:) = ( 1. + cs%FGNV_scale ) * sfn_unlim_v(i,:)
1134 call streamfn_solver(nz, c2_h_v(i,:), hn2_v(i,:), sfn_unlim_v(i,:))
1136 sfn_unlim_v(i,:) = 0.
1143 if (k > nk_linear)
then
1146 if (vhtot(i,j) <= 0.0)
then
1148 sfn_safe = vhtot(i,j) * (1.0 - h_frac(i,j,k)) * gv%H_to_Z
1150 sfn_safe = vhtot(i,j) * (1.0 - h_frac(i,j+1,k)) * gv%H_to_Z
1154 sfn_est = (sfn_unlim_v(i,k) + slope2_ratio_v(i,k)*sfn_safe) / (1.0 + slope2_ratio_v(i,k))
1156 sfn_est = sfn_unlim_v(i,k)
1161 sfn_in_h = min(max(sfn_est * gv%Z_to_H, -h_avail_rsum(i,j,k)), h_avail_rsum(i,j+1,k))
1165 vhd(i,j,k) = max(min((sfn_in_h - vhtot(i,j)), h_avail(i,j,k)), &
1168 if (cs%id_sfn_y>0) diag_sfn_y(i,j,k) = diag_sfn_y(i,j,k+1) + vhd(i,j,k)
1180 if (vhtot(i,j) <= 0.0)
then
1181 vhd(i,j,k) = -vhtot(i,j) * h_frac(i,j,k)
1183 vhd(i,j,k) = -vhtot(i,j) * h_frac(i,j+1,k)
1186 if (cs%id_sfn_y>0) diag_sfn_y(i,j,k) = diag_sfn_y(i,j,k+1) + vhd(i,j,k)
1194 vhtot(i,j) = vhtot(i,j) + vhd(i,j,k)
1203 work_v(i,j) = work_v(i,j) + g_scale * &
1204 ( vhtot(i,j) * drdkde_v(i,k) - &
1205 (vhd(i,j,k) * drdj_v(i,k)) * 0.25 * &
1206 ((e(i,j,k) + e(i,j,k+1)) + (e(i,j+1,k) + e(i,j+1,k+1))) )
1214 if (.not.find_work .or. .not.(use_eos))
then
1215 do j=js,je ;
do i=is-1,ie ; uhd(i,j,1) = -uhtot(i,j) ;
enddo ;
enddo
1216 do j=js-1,je ;
do i=is,ie ; vhd(i,j,1) = -vhtot(i,j) ;
enddo ;
enddo
1222 pres_u(i) = 0.5*(pres(i,j,1) + pres(i+1,j,1))
1223 t_u(i) = 0.5*(t(i,j,1) + t(i+1,j,1))
1224 s_u(i) = 0.5*(s(i,j,1) + s(i+1,j,1))
1227 drho_ds_u, (is-isdb+1)-1, ie-is+2, tv%eqn_of_state)
1230 uhd(i,j,1) = -uhtot(i,j)
1233 drdib = drho_dt_u(i) * (t(i+1,j,1)-t(i,j,1)) + &
1234 drho_ds_u(i) * (s(i+1,j,1)-s(i,j,1))
1236 work_u(i,j) = work_u(i,j) + g_scale * &
1237 ( (uhd(i,j,1) * drdib) * 0.25 * &
1238 ((e(i,j,1) + e(i,j,2)) + (e(i+1,j,1) + e(i+1,j,2))) )
1247 pres_v(i) = 0.5*(pres(i,j,1) + pres(i,j+1,1))
1248 t_v(i) = 0.5*(t(i,j,1) + t(i,j+1,1))
1249 s_v(i) = 0.5*(s(i,j,1) + s(i,j+1,1))
1252 drho_ds_v, is, ie-is+1, tv%eqn_of_state)
1255 vhd(i,j,1) = -vhtot(i,j)
1258 drdjb = drho_dt_v(i) * (t(i,j+1,1)-t(i,j,1)) + &
1259 drho_ds_v(i) * (s(i,j+1,1)-s(i,j,1))
1261 work_v(i,j) = work_v(i,j) - g_scale * &
1262 ( (vhd(i,j,1) * drdjb) * 0.25 * &
1263 ((e(i,j,1) + e(i,j,2)) + (e(i,j+1,1) + e(i,j+1,2))) )
1270 if (find_work)
then ;
do j=js,je ;
do i=is,ie
1272 work_h = 0.5 * g%IareaT(i,j) * &
1273 ((work_u(i-1,j) + work_u(i,j)) + (work_v(i,j-1) + work_v(i,j)))
1274 pe_release_h = -0.25*(kh_u(i,j,k)*(slope_x_pe(i,j,k)**2) * hn2_x_pe(i,j,k) + &
1275 kh_u(i-1,j,k)*(slope_x_pe(i-1,j,k)**2) * hn2_x_pe(i-1,j,k) + &
1276 kh_v(i,j,k)*(slope_y_pe(i,j,k)**2) * hn2_y_pe(i,j,k) + &
1277 kh_v(i,j-1,k)*(slope_y_pe(i,j-1,k)**2) * hn2_y_pe(i,j-1,k))
1278 if (
associated(cs%GMwork)) cs%GMwork(i,j) = work_h
1279 if (
associated(meke))
then ;
if (
associated(meke%GM_src))
then
1280 if (cs%GM_src_alt)
then
1281 meke%GM_src(i,j) = meke%GM_src(i,j) + pe_release_h
1283 meke%GM_src(i,j) = meke%GM_src(i,j) + work_h
1287 enddo ;
enddo ;
endif
1289 if (cs%id_slope_x > 0)
call post_data(cs%id_slope_x, cs%diagSlopeX, cs%diag)
1290 if (cs%id_slope_y > 0)
call post_data(cs%id_slope_y, cs%diagSlopeY, cs%diag)
1291 if (cs%id_sfn_x > 0)
call post_data(cs%id_sfn_x, diag_sfn_x, cs%diag)
1292 if (cs%id_sfn_y > 0)
call post_data(cs%id_sfn_y, diag_sfn_y, cs%diag)
1293 if (cs%id_sfn_unlim_x > 0)
call post_data(cs%id_sfn_unlim_x, diag_sfn_unlim_x, cs%diag)
1294 if (cs%id_sfn_unlim_y > 0)
call post_data(cs%id_sfn_unlim_y, diag_sfn_unlim_y, cs%diag)
1296 end subroutine thickness_diffuse_full
1299 subroutine streamfn_solver(nk, c2_h, hN2, sfn)
1300 integer,
intent(in) :: nk
1301 real,
dimension(nk),
intent(in) :: c2_h
1302 real,
dimension(nk+1),
intent(in) :: hN2
1303 real,
dimension(nk+1),
intent(inout) :: sfn
1309 real :: b_denom, beta, d1, c1(nk)
1312 b_denom = hn2(2) + c2_h(1)
1313 beta = 1.0 / ( b_denom + c2_h(2) )
1315 sfn(2) = ( beta * hn2(2) )*sfn(2)
1317 c1(k-1) = beta * c2_h(k-1)
1318 b_denom = hn2(k) + d1*c2_h(k-1)
1319 beta = 1.0 / (b_denom + c2_h(k))
1321 sfn(k) = beta * (hn2(k)*sfn(k) + c2_h(k-1)*sfn(k-1))
1323 c1(nk) = beta * c2_h(nk)
1326 sfn(k) = sfn(k) + c1(k)*sfn(k+1)
1329 end subroutine streamfn_solver
1332 subroutine add_detangling_kh(h, e, Kh_u, Kh_v, KH_u_CFL, KH_v_CFL, tv, dt, G, GV, US, CS, &
1333 int_slope_u, int_slope_v)
1337 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
1338 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1),
intent(in) :: e
1339 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)+1),
intent(inout) :: Kh_u
1341 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)+1),
intent(inout) :: Kh_v
1343 real,
dimension(SZIB_(G),SZJ_(G)),
intent(in) :: Kh_u_CFL
1345 real,
dimension(SZI_(G),SZJB_(G)),
intent(in) :: Kh_v_CFL
1348 real,
intent(in) :: dt
1350 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)+1),
intent(inout) :: int_slope_u
1354 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)+1),
intent(inout) :: int_slope_v
1359 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
1362 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: &
1365 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)) :: &
1368 real,
dimension(SZI_(G),SZJ_(G)) :: &
1401 real :: denom, I_denom
1410 real,
dimension(SZIB_(G),SZK_(G)+1) :: &
1433 real,
dimension(SZIB_(G)) :: &
1435 logical,
dimension(SZIB_(G)) :: &
1437 integer :: i, j, k, n, ish, jsh, is, ie, js, je, nz, k_top
1438 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1440 k_top = gv%nk_rho_varies + 1
1441 h_neglect = gv%H_subroundoff
1446 if (cs%detangle_time > dt) kh_scale = 0.5 * dt / cs%detangle_time
1448 do j=js-1,je+1 ;
do i=is-1,ie+1
1449 de_top(i,j,k_top) = 0.0 ; de_bot(i,j) = 0.0
1451 do k=k_top+1,nz ;
do j=js-1,je+1 ;
do i=is-1,ie+1
1452 de_top(i,j,k) = de_top(i,j,k-1) + h(i,j,k-1)
1453 enddo ;
enddo ;
enddo
1455 do j=js,je ;
do i=is-1,ie
1456 kh_lay_u(i,j,nz) = 0.0 ; kh_lay_u(i,j,k_top) = 0.0
1458 do j=js-1,je ;
do i=is,ie
1459 kh_lay_v(i,j,nz) = 0.0 ; kh_lay_v(i,j,k_top) = 0.0
1462 do k=nz-1,k_top+1,-1
1464 do j=js-1,je+1 ;
do i=is-1,ie+1
1465 de_bot(i,j) = de_bot(i,j) + h(i,j,k+1)
1468 do j=js,je ;
do i=is-1,ie ;
if (g%mask2dCu(i,j) > 0.0)
then
1469 if (h(i,j,k) > h(i+1,j,k))
then
1471 h1 = max( h(i+1,j,k), h2 - min(de_bot(i+1,j), de_top(i+1,j,k)) )
1474 h1 = max( h(i,j,k), h2 - min(de_bot(i,j), de_top(i,j,k)) )
1476 jag_rat = (h2 - h1)**2 / (h2 + h1 + h_neglect)**2
1477 kh_lay_u(i,j,k) = (kh_scale * kh_u_cfl(i,j)) * jag_rat**2
1478 endif ;
enddo ;
enddo
1480 do j=js-1,je ;
do i=is,ie ;
if (g%mask2dCv(i,j) > 0.0)
then
1481 if (h(i,j,k) > h(i,j+1,k))
then
1483 h1 = max( h(i,j+1,k), h2 - min(de_bot(i,j+1), de_top(i,j+1,k)) )
1486 h1 = max( h(i,j,k), h2 - min(de_bot(i,j), de_top(i,j,k)) )
1488 jag_rat = (h2 - h1)**2 / (h2 + h1 + h_neglect)**2
1489 kh_lay_v(i,j,k) = (kh_scale * kh_v_cfl(i,j)) * jag_rat**2
1490 endif ;
enddo ;
enddo
1495 i_4t = us%Z_to_m*kh_scale / (4.0*dt)
1498 if (n==1)
then ; jsh = js ; ish = is-1
1499 else ; jsh = js-1 ; ish = is ;
endif
1506 do_i(i) = (g%mask2dCu(i,j) > 0.0)
1507 kh_max_max(i) = kh_u_cfl(i,j)
1509 do k=1,nz+1 ;
do i=ish,ie
1510 kh_bg(i,k) = kh_u(i,j,k) ; kh(i,k) = kh_bg(i,k)
1511 kh_min_max_p(i,k) = kh_bg(i,k) ; kh_min_max_m(i,k) = kh_bg(i,k)
1512 kh_detangle(i,k) = 0.0
1516 do_i(i) = (g%mask2dCv(i,j) > 0.0) ; kh_max_max(i) = kh_v_cfl(i,j)
1518 do k=1,nz+1 ;
do i=ish,ie
1519 kh_bg(i,k) = kh_v(i,j,k) ; kh(i,k) = kh_bg(i,k)
1520 kh_min_max_p(i,k) = kh_bg(i,k) ; kh_min_max_m(i,k) = kh_bg(i,k)
1521 kh_detangle(i,k) = 0.0
1526 do k=k_top,nz ;
do i=ish,ie ;
if (do_i(i))
then
1529 denom = ((g%IareaT(i+1,j) + g%IareaT(i,j))*g%dy_Cu(i,j))
1533 dh = i_4t * ((e(i+1,j,k) - e(i+1,j,k+1)) - &
1534 (e(i,j,k) - e(i,j,k+1))) / denom
1538 sign = 1.0*us%Z_to_m ;
if (dh < 0) sign = -1.0*us%Z_to_m
1539 sl_k = sign * (e(i+1,j,k)-e(i,j,k)) * g%IdxCu(i,j)
1540 sl_kp1 = sign * (e(i+1,j,k+1)-e(i,j,k+1)) * g%IdxCu(i,j)
1545 denom = (sl_k**2 + sl_kp1**2)
1546 wt1 = 0.5 ; wt2 = 0.5
1547 if (denom > 0.0)
then
1548 wt1 = sl_k**2 / denom ; wt2 = sl_kp1**2 / denom
1550 kh_detangle(i,k) = kh_detangle(i,k) + wt1*kh_lay_u(i,j,k)
1551 kh_detangle(i,k+1) = kh_detangle(i,k+1) + wt2*kh_lay_u(i,j,k)
1554 denom = ((g%IareaT(i,j+1) + g%IareaT(i,j))*g%dx_Cv(i,j))
1556 dh = i_4t * ((e(i,j+1,k) - e(i,j+1,k+1)) - &
1557 (e(i,j,k) - e(i,j,k+1))) / denom
1561 sign = 1.0*us%Z_to_m ;
if (dh < 0) sign = -1.0*us%Z_to_m
1562 sl_k = sign * (e(i,j+1,k)-e(i,j,k)) * g%IdyCv(i,j)
1563 sl_kp1 = sign * (e(i,j+1,k+1)-e(i,j,k+1)) * g%IdyCv(i,j)
1568 denom = (sl_k**2 + sl_kp1**2)
1569 wt1 = 0.5 ; wt2 = 0.5
1570 if (denom > 0.0)
then
1571 wt1 = sl_k**2 / denom ; wt2 = sl_kp1**2 / denom
1573 kh_detangle(i,k) = kh_detangle(i,k) + wt1*kh_lay_v(i,j,k)
1574 kh_detangle(i,k+1) = kh_detangle(i,k+1) + wt2*kh_lay_v(i,j,k)
1577 if (adh == 0.0)
then
1578 kh_min_m(i,k+1) = 1.0 ; kh0_min_m(i,k+1) = 0.0
1579 kh_max_m(i,k+1) = 1.0 ; kh0_max_m(i,k+1) = 0.0
1580 kh_min_p(i,k) = 1.0 ; kh0_min_p(i,k) = 0.0
1581 kh_max_p(i,k) = 1.0 ; kh0_max_p(i,k) = 0.0
1582 elseif (adh > 0.0)
then
1583 if (sl_k <= sl_kp1)
then
1586 kh_min_m(i,k+1) = 1.0 ; kh0_min_m(i,k+1) = 0.0
1587 kh_max_m(i,k+1) = 1.0 ; kh0_max_m(i,k+1) = 0.0
1588 kh_min_p(i,k) = 1.0 ; kh0_min_p(i,k) = 0.0
1589 kh_max_p(i,k) = 1.0 ; kh0_max_p(i,k) = 0.0
1590 elseif (sl_k <= 0.0)
then
1591 i_sl = -1.0 / sl_kp1
1593 irsl = 1e9 ;
if (rsl > 1e-9) irsl = 1.0/rsl
1596 if (kh_max_max(i) > 0) &
1597 fn_r = min(sqrt(rsl), rsl + (adh * i_sl) / kh_max_max(i))
1599 kh_min_m(i,k+1) = fn_r ; kh0_min_m(i,k+1) = 0.0
1600 kh_max_m(i,k+1) = rsl ; kh0_max_m(i,k+1) = adh * i_sl
1601 kh_min_p(i,k) = irsl ; kh0_min_p(i,k) = -adh * (i_sl*irsl)
1602 kh_max_p(i,k) = 1.0/(fn_r + 1.0e-30) ; kh0_max_p(i,k) = 0.0
1603 elseif (sl_kp1 < 0.0)
then
1604 i_sl_k = 1e18 ;
if (sl_k > 1e-18) i_sl_k = 1.0 / sl_k
1605 i_sl_kp1 = 1e18 ;
if (-sl_kp1 > 1e-18) i_sl_kp1 = -1.0 / sl_kp1
1607 kh_min_m(i,k+1) = 0.0 ; kh0_min_m(i,k+1) = 0.0
1608 kh_max_m(i,k+1) = - sl_k*i_sl_kp1 ; kh0_max_m(i,k+1) = adh*i_sl_kp1
1609 kh_min_p(i,k) = 0.0 ; kh0_min_p(i,k) = 0.0
1610 kh_max_p(i,k) = sl_kp1*i_sl_k ; kh0_max_p(i,k) = adh*i_sl_k
1614 kh_max = adh / (sl_k - sl_kp1)
1615 kh_min_max_p(i,k) = max(kh_min_max_p(i,k), kh_max)
1616 kh_min_max_m(i,k+1) = max(kh_min_max_m(i,k+1), kh_max)
1620 irsl = 1e9 ;
if (rsl > 1e-9) irsl = 1.0/rsl
1624 if (kh_max_max(i) > 0) &
1625 fn_r = min(sqrt(rsl), rsl + (adh * i_sl) / kh_max_max(i))
1627 kh_min_m(i,k+1) = irsl ; kh0_min_m(i,k+1) = -adh * (i_sl*irsl)
1628 kh_max_m(i,k+1) = 1.0/(fn_r + 1.0e-30) ; kh0_max_m(i,k+1) = 0.0
1629 kh_min_p(i,k) = fn_r ; kh0_min_p(i,k) = 0.0
1630 kh_max_p(i,k) = rsl ; kh0_max_p(i,k) = adh * i_sl
1633 endif ;
enddo ;
enddo
1635 do k=k_top,nz+1,nz+1-k_top ;
do i=ish,ie ;
if (do_i(i))
then
1637 kh_min_m(i,k) = 0.0 ; kh0_min_m(i,k) = 0.0
1638 kh_max_m(i,k) = 0.0 ; kh0_max_m(i,k) = 0.0
1639 kh_min_p(i,k) = 0.0 ; kh0_min_p(i,k) = 0.0
1640 kh_max_p(i,k) = 0.0 ; kh0_max_p(i,k) = 0.0
1641 kh_min_max_p(i,k) = kh_bg(i,k)
1642 kh_min_max_m(i,k) = kh_bg(i,k)
1643 endif ;
enddo ;
enddo
1653 do k=k_top+1,nz ;
do i=ish,ie ;
if (do_i(i))
then
1654 kh(i,k) = max(kh_bg(i,k), kh_detangle(i,k), &
1655 min(kh_min_m(i,k)*kh(i,k-1) + kh0_min_m(i,k), kh(i,k-1)))
1657 if (kh0_max_m(i,k) > kh_bg(i,k)) kh(i,k) = min(kh(i,k), kh0_max_m(i,k))
1658 if (kh0_max_p(i,k) > kh_bg(i,k)) kh(i,k) = min(kh(i,k), kh0_max_p(i,k))
1659 endif ;
enddo ;
enddo
1661 do k=nz,k_top+1,-1 ;
do i=ish,ie ;
if (do_i(i))
then
1662 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)))
1664 kh_max = max(kh_min_max_p(i,k), kh_max_p(i,k)*kh(i,k+1) + kh0_max_p(i,k))
1665 kh(i,k) = min(kh(i,k), kh_max)
1666 endif ;
enddo ;
enddo
1671 do k=k_top+1,nz ;
do i=ish,ie ;
if (do_i(i))
then
1672 kh_max = max(kh_min_max_m(i,k), kh_max_m(i,k)*kh(i,k-1) + kh0_max_m(i,k))
1673 if (kh(i,k) > kh_max) kh(i,k) = kh_max
1674 endif ;
enddo ;
enddo
1728 do k=k_top+1,nz ;
do i=ish,ie
1729 if (kh(i,k) > kh_u(i,j,k))
then
1730 dkh = (kh(i,k) - kh_u(i,j,k))
1731 int_slope_u(i,j,k) = dkh / kh(i,k)
1732 kh_u(i,j,k) = kh(i,k)
1736 do k=k_top+1,nz ;
do i=ish,ie
1737 if (kh(i,k) > kh_v(i,j,k))
then
1738 dkh = kh(i,k) - kh_v(i,j,k)
1739 int_slope_v(i,j,k) = dkh / kh(i,k)
1740 kh_v(i,j,k) = kh(i,k)
1748 end subroutine add_detangling_kh
1751 subroutine thickness_diffuse_init(Time, G, GV, US, param_file, diag, CDp, CS)
1752 type(time_type),
intent(in) :: time
1757 type(
diag_ctrl),
target,
intent(inout) :: diag
1762 #include "version_variable.h"
1763 character(len=40) :: mdl =
"MOM_thickness_diffuse"
1764 real :: omega, strat_floor, flux_to_kg_per_s
1766 if (
associated(cs))
then
1767 call mom_error(warning, &
1768 "Thickness_diffuse_init called with an associated control structure.")
1770 else ;
allocate(cs) ;
endif
1776 call get_param(param_file, mdl,
"THICKNESSDIFFUSE", cs%thickness_diffuse, &
1777 "If true, interface heights are diffused with a "//&
1778 "coefficient of KHTH.", default=.false.)
1779 call get_param(param_file, mdl,
"KHTH", cs%Khth, &
1780 "The background horizontal thickness diffusivity.", &
1781 units =
"m2 s-1", default=0.0)
1782 call get_param(param_file, mdl,
"KHTH_SLOPE_CFF", cs%KHTH_Slope_Cff, &
1783 "The nondimensional coefficient in the Visbeck formula "//&
1784 "for the interface depth diffusivity", units=
"nondim", &
1786 call get_param(param_file, mdl,
"KHTH_MIN", cs%KHTH_Min, &
1787 "The minimum horizontal thickness diffusivity.", &
1788 units =
"m2 s-1", default=0.0)
1789 call get_param(param_file, mdl,
"KHTH_MAX", cs%KHTH_Max, &
1790 "The maximum horizontal thickness diffusivity.", &
1791 units =
"m2 s-1", default=0.0)
1792 call get_param(param_file, mdl,
"KHTH_MAX_CFL", cs%max_Khth_CFL, &
1793 "The maximum value of the local diffusive CFL ratio that "//&
1794 "is permitted for the thickness diffusivity. 1.0 is the "//&
1795 "marginally unstable value in a pure layered model, but "//&
1796 "much smaller numbers (e.g. 0.1) seem to work better for "//&
1797 "ALE-based models.", units =
"nondimensional", default=0.8)
1803 if (cs%max_Khth_CFL < 0.0) cs%max_Khth_CFL = 0.0
1804 call get_param(param_file, mdl,
"DETANGLE_INTERFACES", cs%detangle_interfaces, &
1805 "If defined add 3-d structured enhanced interface height "//&
1806 "diffusivities to horizontally smooth jagged layers.", &
1808 cs%detangle_time = 0.0
1809 if (cs%detangle_interfaces) &
1810 call get_param(param_file, mdl,
"DETANGLE_TIMESCALE", cs%detangle_time, &
1811 "A timescale over which maximally jagged grid-scale "//&
1812 "thickness variations are suppressed. This must be "//&
1813 "longer than DT, or 0 to use DT.", units =
"s", default=0.0)
1814 call get_param(param_file, mdl,
"KHTH_SLOPE_MAX", cs%slope_max, &
1815 "A slope beyond which the calculated isopycnal slope is "//&
1816 "not reliable and is scaled away.", units=
"nondim", default=0.01)
1817 call get_param(param_file, mdl,
"KD_SMOOTH", cs%kappa_smooth, &
1818 "A diapycnal diffusivity that is used to interpolate "//&
1819 "more sensible values of T & S into thin layers.", &
1820 default=1.0e-6, scale=us%m_to_Z**2)
1821 call get_param(param_file, mdl,
"KHTH_USE_FGNV_STREAMFUNCTION", cs%use_FGNV_streamfn, &
1822 "If true, use the streamfunction formulation of "//&
1823 "Ferrari et al., 2010, which effectively emphasizes "//&
1824 "graver vertical modes by smoothing in the vertical.", &
1826 call get_param(param_file, mdl,
"FGNV_FILTER_SCALE", cs%FGNV_scale, &
1827 "A coefficient scaling the vertical smoothing term in the "//&
1828 "Ferrari et al., 2010, streamfunction formulation.", &
1829 default=1., do_not_log=.not.cs%use_FGNV_streamfn)
1830 call get_param(param_file, mdl,
"FGNV_C_MIN", cs%FGNV_c_min, &
1831 "A minium wave speed used in the Ferrari et al., 2010, "//&
1832 "streamfunction formulation.", &
1833 default=0., units=
"m s-1", do_not_log=.not.cs%use_FGNV_streamfn)
1834 call get_param(param_file, mdl,
"FGNV_STRAT_FLOOR", strat_floor, &
1835 "A floor for Brunt-Vasaila frequency in the Ferrari et al., 2010, "//&
1836 "streamfunction formulation, expressed as a fraction of planetary "//&
1837 "rotation, OMEGA. This should be tiny but non-zero to avoid degeneracy.", &
1838 default=1.e-15, units=
"nondim", do_not_log=.not.cs%use_FGNV_streamfn)
1839 call get_param(param_file, mdl,
"OMEGA",omega, &
1840 "The rotation rate of the earth.", units=
"s-1", &
1841 default=7.2921e-5, do_not_log=.not.cs%use_FGNV_streamfn)
1842 if (cs%use_FGNV_streamfn) cs%N2_floor = (strat_floor*omega)**2
1843 call get_param(param_file, mdl,
"DEBUG", cs%debug, &
1844 "If true, write out verbose debugging data.", &
1845 default=.false., debuggingparam=.true.)
1847 call get_param(param_file, mdl,
"MEKE_GM_SRC_ALT", cs%GM_src_alt, &
1848 "If true, use the GM energy conversion form S^2*N^2*kappa rather \n"//&
1849 "than the streamfunction for the GM source term.", default=.false.)
1850 call get_param(param_file, mdl,
"MEKE_GEOMETRIC", cs%MEKE_GEOMETRIC, &
1851 "If true, uses the GM coefficient formulation \n"//&
1852 "from the GEOMETRIC framework (Marshall et al., 2012).", default=.false.)
1853 if (cs%MEKE_GEOMETRIC)
then
1855 call get_param(param_file, mdl,
"MEKE_GEOMETRIC_EPSILON", cs%MEKE_GEOMETRIC_epsilon, &
1856 "Minimum Eady growth rate used in the calculation of \n"//&
1857 "GEOMETRIC thickness diffusivity.", units=
"s-1", default=1.0e-7)
1859 call get_param(param_file, mdl,
"MEKE_GEOMETRIC_ALPHA", cs%MEKE_GEOMETRIC_alpha, &
1860 "The nondimensional coefficient governing the efficiency of the GEOMETRIC \n"//&
1861 "thickness diffusion.", units=
"nondim", default=0.05)
1864 call get_param(param_file, mdl,
"USE_KH_IN_MEKE", cs%Use_KH_in_MEKE, &
1865 "If true, uses the thickness diffusivity calculated here to diffuse \n"//&
1866 "MEKE.", default=.false.)
1868 call get_param(param_file, mdl,
"USE_GME", cs%use_GME_thickness_diffuse, &
1869 "If true, use the GM+E backscatter scheme in association \n"//&
1870 "with the Gent and McWilliams parameterization.", default=.false.)
1872 if (cs%use_GME_thickness_diffuse)
then
1873 call safe_alloc_ptr(cs%KH_u_GME,g%IsdB,g%IedB,g%jsd,g%jed,g%ke+1)
1874 call safe_alloc_ptr(cs%KH_v_GME,g%isd,g%ied,g%JsdB,g%JedB,g%ke+1)
1877 if (gv%Boussinesq)
then ; flux_to_kg_per_s = gv%Rho0
1878 else ; flux_to_kg_per_s = 1. ;
endif
1880 cs%id_uhGM = register_diag_field(
'ocean_model',
'uhGM', diag%axesCuL, time, &
1881 'Time Mean Diffusive Zonal Thickness Flux',
'kg s-1', &
1882 y_cell_method=
'sum', v_extensive=.true., conversion=flux_to_kg_per_s)
1883 if (cs%id_uhGM > 0)
call safe_alloc_ptr(cdp%uhGM,g%IsdB,g%IedB,g%jsd,g%jed,g%ke)
1884 cs%id_vhGM = register_diag_field(
'ocean_model',
'vhGM', diag%axesCvL, time, &
1885 'Time Mean Diffusive Meridional Thickness Flux',
'kg s-1', &
1886 x_cell_method=
'sum', v_extensive=.true., conversion=flux_to_kg_per_s)
1887 if (cs%id_vhGM > 0)
call safe_alloc_ptr(cdp%vhGM,g%isd,g%ied,g%JsdB,g%JedB,g%ke)
1889 cs%id_GMwork = register_diag_field(
'ocean_model',
'GMwork', diag%axesT1, time, &
1890 'Integrated Tendency of Ocean Mesoscale Eddy KE from Parameterized Eddy Advection', &
1891 'W m-2', cmor_field_name=
'tnkebto', &
1892 cmor_long_name=
'Integrated Tendency of Ocean Mesoscale Eddy KE from Parameterized Eddy Advection',&
1893 cmor_standard_name=
'tendency_of_ocean_eddy_kinetic_energy_content_due_to_parameterized_eddy_advection')
1894 if (cs%id_GMwork > 0)
call safe_alloc_ptr(cs%GMwork,g%isd,g%ied,g%jsd,g%jed)
1896 cs%id_KH_u = register_diag_field(
'ocean_model',
'KHTH_u', diag%axesCui, time, &
1897 'Parameterized mesoscale eddy advection diffusivity at U-point',
'm2 s-1')
1898 cs%id_KH_v = register_diag_field(
'ocean_model',
'KHTH_v', diag%axesCvi, time, &
1899 'Parameterized mesoscale eddy advection diffusivity at V-point',
'm2 s-1')
1900 cs%id_KH_t = register_diag_field(
'ocean_model',
'KHTH_t', diag%axesTL, time, &
1901 'Ocean Tracer Diffusivity due to Parameterized Mesoscale Advection',
'm2 s-1',&
1902 cmor_field_name=
'diftrblo', &
1903 cmor_long_name=
'Ocean Tracer Diffusivity due to Parameterized Mesoscale Advection', &
1904 cmor_units=
'm2 s-1', &
1905 cmor_standard_name=
'ocean_tracer_diffusivity_due_to_parameterized_mesoscale_advection')
1907 cs%id_KH_u1 = register_diag_field(
'ocean_model',
'KHTH_u1', diag%axesCu1, time, &
1908 'Parameterized mesoscale eddy advection diffusivity at U-points (2-D)',
'm2 s-1')
1909 cs%id_KH_v1 = register_diag_field(
'ocean_model',
'KHTH_v1', diag%axesCv1, time, &
1910 'Parameterized mesoscale eddy advection diffusivity at V-points (2-D)',
'm2 s-1')
1911 cs%id_KH_t1 = register_diag_field(
'ocean_model',
'KHTH_t1', diag%axesT1, time,&
1912 'Parameterized mesoscale eddy advection diffusivity at T-points (2-D)',
'm2 s-1')
1914 cs%id_slope_x = register_diag_field(
'ocean_model',
'neutral_slope_x', diag%axesCui, time, &
1915 'Zonal slope of neutral surface',
'nondim')
1916 if (cs%id_slope_x > 0)
call safe_alloc_ptr(cs%diagSlopeX,g%IsdB,g%IedB,g%jsd,g%jed,g%ke+1)
1917 cs%id_slope_y = register_diag_field(
'ocean_model',
'neutral_slope_y', diag%axesCvi, time, &
1918 'Meridional slope of neutral surface',
'nondim')
1919 if (cs%id_slope_y > 0)
call safe_alloc_ptr(cs%diagSlopeY,g%isd,g%ied,g%JsdB,g%JedB,g%ke+1)
1920 cs%id_sfn_x = register_diag_field(
'ocean_model',
'GM_sfn_x', diag%axesCui, time, &
1921 'Parameterized Zonal Overturning Streamfunction',
'm3 s-1')
1922 cs%id_sfn_y = register_diag_field(
'ocean_model',
'GM_sfn_y', diag%axesCvi, time, &
1923 'Parameterized Meridional Overturning Streamfunction',
'm3 s-1')
1924 cs%id_sfn_unlim_x = register_diag_field(
'ocean_model',
'GM_sfn_unlim_x', diag%axesCui, time, &
1925 'Parameterized Zonal Overturning Streamfunction before limiting/smoothing', &
1926 'm3 s-1', conversion=us%Z_to_m)
1927 cs%id_sfn_unlim_y = register_diag_field(
'ocean_model',
'GM_sfn_unlim_y', diag%axesCvi, time, &
1928 'Parameterized Meridional Overturning Streamfunction before limiting/smoothing', &
1929 'm3 s-1', conversion=us%Z_to_m)
1931 end subroutine thickness_diffuse_init
1934 subroutine thickness_diffuse_get_kh(CS, KH_u_GME, KH_v_GME, G)
1938 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)+1),
intent(inout) :: kh_u_gme
1940 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)+1),
intent(inout) :: kh_v_gme
1945 do k=1,g%ke+1 ;
do j = g%jsc, g%jec ;
do i = g%isc-1, g%iec
1946 kh_u_gme(i,j,k) = cs%KH_u_GME(i,j,k)
1947 enddo ;
enddo ;
enddo
1949 do k=1,g%ke+1 ;
do j = g%jsc-1, g%jec ;
do i = g%isc, g%iec
1950 kh_v_gme(i,j,k) = cs%KH_v_GME(i,j,k)
1951 enddo ;
enddo ;
enddo
1953 end subroutine thickness_diffuse_get_kh
1956 subroutine thickness_diffuse_end(CS)
1959 if (
associated(cs))
deallocate(cs)
1960 end subroutine thickness_diffuse_end