10 use mom_cpu_clock,
only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
11 use mom_cpu_clock,
only : clock_module_driver, clock_module, clock_routine
15 use mom_diabatic_aux,
only : make_frazil, adjust_salt, insert_brine, differential_diffuse_t_s, tridiagts
16 use mom_diabatic_aux,
only : find_uv_at_h, diagnosemldbydensitydifference, applyboundaryfluxesinout
22 use mom_diag_mediator,
only : diag_copy_diag_to_storage, diag_copy_storage_to_diag
38 use mom_eos,
only : calculate_specific_vol_derivs
39 use mom_error_handler,
only : mom_error, fatal, warning, calltree_showquery,mom_mesg
43 use mom_forcing_type,
only : calculatebuoyancyflux2d, forcing_singlepointprint
54 use mom_cvmix_kpp,
only : kpp_nonlocaltransport_temp, kpp_nonlocaltransport_saln
63 use mom_time_manager,
only : time_type, real_to_time,
operator(-),
operator(<=)
74 implicit none ;
private
76 #include <MOM_memory.h>
79 public diabatic_driver_init
80 public diabatic_driver_end
81 public extract_diabatic_member
83 public adiabatic_driver_init
94 logical :: use_legacy_diabatic
97 logical :: bulkmixedlayer
99 logical :: use_energetic_pbl
104 logical :: use_kappa_shear
106 logical :: use_cvmix_shear
108 logical :: use_cvmix_ddiff
109 logical :: use_tidal_mixing
110 logical :: use_cvmix_conv
112 logical :: use_sponge
116 logical :: use_geothermal
117 logical :: use_int_tides
119 logical :: epbl_is_additive
123 real :: uniform_test_cg
125 logical :: usealealgorithm
128 logical :: aggregate_fw_forcing
138 logical :: massless_match_targets
140 logical :: mix_boundary_tracers
151 real :: minimum_forcing_depth = 0.001
153 real :: evap_cfl_limit = 0.8
155 integer :: halo_ts_diff = 0
157 logical :: usekpp = .false.
158 logical :: salt_reject_below_ml
159 logical :: kppispassive
161 logical :: debugconservation
162 logical :: tracer_tridiag
164 logical :: debug_energy_req
166 real :: mlddensitydifference
171 integer :: id_cg1 = -1
172 integer,
allocatable,
dimension(:) :: id_cn
173 integer :: id_wd = -1, id_ea = -1, id_eb = -1
174 integer :: id_dudt_dia = -1, id_dvdt_dia = -1, id_ea_s = -1, id_eb_s = -1
175 integer :: id_ea_t = -1, id_eb_t = -1
176 integer :: id_kd_heat = -1, id_kd_salt = -1, id_kd_interface = -1, id_kd_epbl = -1
177 integer :: id_tdif = -1, id_tadv = -1, id_sdif = -1, id_sadv = -1
178 integer :: id_mld_003 = -1, id_mld_0125 = -1, id_mld_user = -1, id_mlotstsq = -1
179 integer :: id_submln2 = -1, id_brine_lay = -1
182 integer :: id_u_predia = -1, id_v_predia = -1, id_h_predia = -1
183 integer :: id_t_predia = -1, id_s_predia = -1, id_e_predia = -1
185 integer :: id_diabatic_diff_temp_tend = -1
186 integer :: id_diabatic_diff_saln_tend = -1
187 integer :: id_diabatic_diff_heat_tend = -1
188 integer :: id_diabatic_diff_salt_tend = -1
189 integer :: id_diabatic_diff_heat_tend_2d = -1
190 integer :: id_diabatic_diff_salt_tend_2d = -1
191 integer :: id_diabatic_diff_h= -1
193 integer :: id_boundary_forcing_h = -1
194 integer :: id_boundary_forcing_h_tendency = -1
195 integer :: id_boundary_forcing_temp_tend = -1
196 integer :: id_boundary_forcing_saln_tend = -1
197 integer :: id_boundary_forcing_heat_tend = -1
198 integer :: id_boundary_forcing_salt_tend = -1
199 integer :: id_boundary_forcing_heat_tend_2d = -1
200 integer :: id_boundary_forcing_salt_tend_2d = -1
202 integer :: id_frazil_h = -1
203 integer :: id_frazil_temp_tend = -1
204 integer :: id_frazil_heat_tend = -1
205 integer :: id_frazil_heat_tend_2d = -1
208 logical :: diabatic_diff_tendency_diag = .false.
209 logical :: boundary_forcing_tendency_diag = .false.
210 logical :: frazil_tendency_diag = .false.
211 real,
allocatable,
dimension(:,:,:) :: frazil_heat_diag
212 real,
allocatable,
dimension(:,:,:) :: frazil_temp_diag
229 type(
kpp_cs),
pointer :: kpp_csp => null()
234 type(group_pass_type) :: pass_hold_eb_ea
235 type(group_pass_type) :: pass_kv
238 real,
allocatable,
dimension(:,:,:) :: kpp_nltheat
239 real,
allocatable,
dimension(:,:,:) :: kpp_nltscalar
240 real,
allocatable,
dimension(:,:,:) :: kpp_buoy_flux
241 real,
allocatable,
dimension(:,:) :: kpp_temp_flux
242 real,
allocatable,
dimension(:,:) :: kpp_salt_flux
244 type(time_type),
pointer :: time
248 integer :: id_clock_entrain, id_clock_mixedlayer, id_clock_set_diffusivity
249 integer :: id_clock_tracers, id_clock_tridiag, id_clock_pass, id_clock_sponge
250 integer :: id_clock_geothermal, id_clock_differential_diff, id_clock_remap
251 integer :: id_clock_kpp
257 subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, &
258 G, GV, US, CS, WAVES)
262 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(inout) :: u
263 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(inout) :: v
264 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(inout) :: h
267 real,
dimension(:,:),
pointer :: hml
268 type(
forcing),
intent(inout) :: fluxes
275 real,
intent(in) :: dt
276 type(time_type),
intent(in) :: time_end
281 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: &
283 real,
dimension(SZI_(G),SZJ_(G),CS%nMode) :: &
285 real,
dimension(SZI_(G),SZJ_(G),G%ke) :: temp_diag
286 real,
dimension(SZI_(G),SZJ_(G)) :: tke_itidal_input_test
288 integer :: i, j, k, m, is, ie, js, je, nz
289 logical :: showcalltree
291 if (g%ke == 1)
return
293 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
295 if (.not.
associated(cs))
call mom_error(fatal,
"MOM_diabatic_driver: "// &
296 "Module must be initialized before it is used.")
297 if (dt == 0.0)
call mom_error(fatal,
"MOM_diabatic_driver: "// &
298 "diabatic was called with a zero length timestep.")
299 if (dt < 0.0)
call mom_error(fatal,
"MOM_diabatic_driver: "// &
300 "diabatic was called with a negative timestep.")
302 showcalltree = calltree_showquery()
306 if (cs%id_u_predia > 0)
call post_data(cs%id_u_predia, u, cs%diag)
307 if (cs%id_v_predia > 0)
call post_data(cs%id_v_predia, v, cs%diag)
308 if (cs%id_h_predia > 0)
call post_data(cs%id_h_predia, h, cs%diag)
309 if (cs%id_T_predia > 0)
call post_data(cs%id_T_predia, tv%T, cs%diag)
310 if (cs%id_S_predia > 0)
call post_data(cs%id_S_predia, tv%S, cs%diag)
311 if (cs%id_e_predia > 0)
then
312 call find_eta(h, tv, g, gv, us, eta, eta_to_m=1.0)
313 call post_data(cs%id_e_predia, eta, cs%diag)
316 dt_in_t = dt * us%s_to_T
319 call mom_forcing_chksum(
"Start of diabatic", fluxes, g, us, haloshift=0)
321 if (cs%debugConservation)
call mom_state_stats(
'Start of diabatic', u, v, h, tv%T, tv%S, g)
323 if (cs%debug_energy_req) &
324 call diapyc_energy_req_test(h, dt_in_t, tv, g, gv, us, cs%diapyc_en_rec_CSp)
326 call cpu_clock_begin(id_clock_set_diffusivity)
327 call set_bbl_tke(u, v, h, fluxes, visc, g, gv, us, cs%set_diff_CSp)
328 call cpu_clock_end(id_clock_set_diffusivity)
333 if (
associated(tv%T) .AND.
associated(tv%frazil))
then
335 call enable_averaging(0.5*dt, time_end - real_to_time(0.5*dt), cs%diag)
336 if (cs%frazil_tendency_diag)
then
337 do k=1,nz ;
do j=js,je ;
do i=is,ie
338 temp_diag(i,j,k) = tv%T(i,j,k)
339 enddo ;
enddo ;
enddo
342 if (
associated(fluxes%p_surf_full))
then
343 call make_frazil(h, tv, g, gv, cs%diabatic_aux_CSp, fluxes%p_surf_full, halo=cs%halo_TS_diff)
345 call make_frazil(h, tv, g, gv, cs%diabatic_aux_CSp, halo=cs%halo_TS_diff)
347 if (showcalltree)
call calltree_waypoint(
"done with 1st make_frazil (diabatic)")
349 if (cs%frazil_tendency_diag)
then
350 call diagnose_frazil_tendency(tv, h, temp_diag, 0.5*dt, g, gv, cs)
351 if (cs%id_frazil_h > 0)
call post_data(cs%id_frazil_h, h, cs%diag)
353 call disable_averaging(cs%diag)
355 if (cs%debugConservation)
call mom_state_stats(
'1st make_frazil', u, v, h, tv%T, tv%S, g)
358 if (cs%use_int_tides)
then
360 call set_int_tide_input(u, v, h, tv, fluxes, cs%int_tide_input, dt, g, gv, us, &
361 cs%int_tide_input_CSp)
363 if (cs%uniform_test_cg > 0.0)
then
364 do m=1,cs%nMode ; cn_igw(:,:,m) = cs%uniform_test_cg ;
enddo
366 call wave_speeds(h, tv, g, gv, us, cs%nMode, cn_igw, full_halos=.true.)
369 call propagate_int_tide(h, tv, cn_igw, cs%int_tide_input%TKE_itidal_input, cs%int_tide_input%tideamp, &
370 cs%int_tide_input%Nb, dt, g, gv, us, cs%int_tide_CSp)
371 if (showcalltree)
call calltree_waypoint(
"done with propagate_int_tide (diabatic)")
375 if (cs%useALEalgorithm .and. cs%use_legacy_diabatic)
then
376 call diabatic_ale_legacy(u, v, h, tv, hml, fluxes, visc, adp, cdp, dt, time_end, &
377 g, gv, us, cs, waves)
378 elseif (cs%useALEalgorithm)
then
379 call diabatic_ale(u, v, h, tv, hml, fluxes, visc, adp, cdp, dt, time_end, &
380 g, gv, us, cs, waves)
382 call layered_diabatic(u, v, h, tv, hml, fluxes, visc, adp, cdp, dt, time_end, &
383 g, gv, us, cs, waves)
388 call cpu_clock_begin(id_clock_pass)
389 if (
associated(visc%Kv_shear)) &
390 call pass_var(visc%Kv_shear, g%Domain, to_all+omit_corners, halo=1)
391 call cpu_clock_end(id_clock_pass)
393 call disable_averaging(cs%diag)
397 if (
associated(tv%T) .AND.
associated(tv%frazil))
then
398 call enable_averaging(0.5*dt, time_end, cs%diag)
399 if (cs%frazil_tendency_diag)
then
400 do k=1,nz ;
do j=js,je ;
do i=is,ie
401 temp_diag(i,j,k) = tv%T(i,j,k)
402 enddo ;
enddo ;
enddo
405 if (
associated(fluxes%p_surf_full))
then
406 call make_frazil(h, tv, g, gv, cs%diabatic_aux_CSp, fluxes%p_surf_full)
408 call make_frazil(h, tv, g, gv, cs%diabatic_aux_CSp)
411 if (cs%frazil_tendency_diag)
then
412 call diagnose_frazil_tendency(tv, h, temp_diag, 0.5*dt, g, gv, cs)
413 if (cs%id_frazil_h > 0 )
call post_data(cs%id_frazil_h, h, cs%diag)
416 if (showcalltree)
call calltree_waypoint(
"done with 2nd make_frazil (diabatic)")
417 if (cs%debugConservation)
call mom_state_stats(
'2nd make_frazil', u, v, h, tv%T, tv%S, g)
418 call disable_averaging(cs%diag)
424 call enable_averaging(dt, time_end, cs%diag)
425 if (cs%id_MLD_003 > 0 .or. cs%id_subMLN2 > 0 .or. cs%id_mlotstsq > 0)
then
426 call diagnosemldbydensitydifference(cs%id_MLD_003, h, tv, 0.03, g, gv, us, cs%diag, &
427 id_n2subml=cs%id_subMLN2, id_mldsq=cs%id_mlotstsq, dz_subml=cs%dz_subML_N2)
429 if (cs%id_MLD_0125 > 0)
then
430 call diagnosemldbydensitydifference(cs%id_MLD_0125, h, tv, 0.125, g, gv, us, cs%diag)
432 if (cs%id_MLD_user > 0)
then
433 call diagnosemldbydensitydifference(cs%id_MLD_user, h, tv, cs%MLDdensityDifference, g, gv, us, cs%diag)
435 if (cs%use_int_tides)
then
436 if (cs%id_cg1 > 0)
call post_data(cs%id_cg1, cn_igw(:,:,1),cs%diag)
437 do m=1,cs%nMode ;
if (cs%id_cn(m) > 0)
call post_data(cs%id_cn(m), cn_igw(:,:,m), cs%diag) ;
enddo
439 call disable_averaging(cs%diag)
441 if (cs%debugConservation)
call mom_state_stats(
'leaving diabatic', u, v, h, tv%T, tv%S, g)
443 end subroutine diabatic
449 subroutine diabatic_ale_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, &
450 G, GV, US, CS, WAVES)
454 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(inout) :: u
455 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(inout) :: v
456 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(inout) :: h
459 real,
dimension(:,:),
pointer :: Hml
460 type(
forcing),
intent(inout) :: fluxes
467 real,
intent(in) :: dt
468 type(time_type),
intent(in) :: Time_end
473 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
493 real,
dimension(SZI_(G),SZJ_(G)) :: &
496 real,
dimension(SZI_(G),SZJ_(G),G%ke) :: h_diag
497 real,
dimension(SZI_(G),SZJ_(G),G%ke) :: temp_diag
498 real,
dimension(SZI_(G),SZJ_(G),G%ke) :: saln_diag
499 real,
dimension(SZI_(G),SZJ_(G)) :: tendency_2d
503 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
508 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1),
target :: &
518 logical :: in_boundary(SZI_(G))
538 real :: htot(SZIB_(G))
539 real :: b1(SZIB_(G)), d1(SZIB_(G))
540 real :: c1(SZIB_(G),SZK_(G))
547 logical :: showCallTree
548 integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb, m, halo
553 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
554 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
555 nkmb = gv%nk_rho_varies
556 h_neglect = gv%H_subroundoff ; h_neglect2 = h_neglect*h_neglect
557 kd_heat(:,:,:) = 0.0 ; kd_salt(:,:,:) = 0.0
559 showcalltree = calltree_showquery()
560 if (showcalltree)
call calltree_enter(
"diabatic_ALE_legacy(), MOM_diabatic_driver.F90")
563 dt_in_t = dt * us%s_to_T
566 call enable_averaging(dt, time_end, cs%diag)
568 if (cs%use_geothermal)
then
569 halo = cs%halo_TS_diff
571 do k=1,nz ;
do j=js-halo,je+halo ;
do i=is-halo,ie+halo
572 h_orig(i,j,k) = h(i,j,k) ; eatr(i,j,k) = 0.0 ; ebtr(i,j,k) = 0.0
573 enddo ;
enddo ;
enddo
576 if (cs%use_geothermal)
then
577 call cpu_clock_begin(id_clock_geothermal)
578 call geothermal(h, tv, dt, eatr, ebtr, g, gv, cs%geothermal_CSp, halo=cs%halo_TS_diff)
579 call cpu_clock_end(id_clock_geothermal)
580 if (showcalltree)
call calltree_waypoint(
"geothermal (diabatic)")
581 if (cs%debugConservation)
call mom_state_stats(
'geothermal', u, v, h, tv%T, tv%S, g)
586 call diag_update_remap_grids(cs%diag)
591 if (
associated(cs%optics)) &
592 call set_pen_shortwave(cs%optics, fluxes, g, gv, cs%diabatic_aux_CSp, cs%opacity_CSp, cs%tracer_flow_CSp)
594 if (cs%debug)
call mom_state_chksum(
"before find_uv_at_h", u, v, h, g, gv, haloshift=0)
596 if (cs%use_kappa_shear .or. cs%use_CVMix_shear)
then
597 if (cs%use_geothermal)
then
598 call find_uv_at_h(u, v, h_orig, u_h, v_h, g, gv, eatr, ebtr)
600 call hchksum(eatr,
"after find_uv_at_h eatr",g%HI, scale=gv%H_to_m)
601 call hchksum(ebtr,
"after find_uv_at_h ebtr",g%HI, scale=gv%H_to_m)
604 call find_uv_at_h(u, v, h, u_h, v_h, g, gv)
606 if (showcalltree)
call calltree_waypoint(
"done with find_uv_at_h (diabatic)")
609 call cpu_clock_begin(id_clock_set_diffusivity)
612 call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, cs%optics, visc, dt_in_t, g, gv, us, &
613 cs%set_diff_CSp, kd_lay, kd_int)
614 call cpu_clock_end(id_clock_set_diffusivity)
615 if (showcalltree)
call calltree_waypoint(
"done with set_diffusivity (diabatic)")
619 if (.not.cs%use_legacy_diabatic .or. cs%useKPP)
then
621 do k=1,nz+1 ;
do j=js,je ;
do i=is,ie
622 kd_salt(i,j,k) = kd_int(i,j,k)
623 kd_heat(i,j,k) = kd_int(i,j,k)
624 enddo ;
enddo ;
enddo
626 if (
associated(visc%Kd_extra_S))
then
628 do k=1,nz+1 ;
do j=js,je ;
do i=is,ie
629 kd_salt(i,j,k) = kd_salt(i,j,k) + visc%Kd_extra_S(i,j,k)
630 enddo ;
enddo ;
enddo
632 if (
associated(visc%Kd_extra_T))
then
634 do k=1,nz+1 ;
do j=js,je ;
do i=is,ie
635 kd_heat(i,j,k) = kd_heat(i,j,k) + visc%Kd_extra_T(i,j,k)
636 enddo ;
enddo ;
enddo
641 call mom_state_chksum(
"after set_diffusivity ", u, v, h, g, gv, haloshift=0)
642 call mom_forcing_chksum(
"after set_diffusivity ", fluxes, g, us, haloshift=0)
643 call mom_thermovar_chksum(
"after set_diffusivity ", tv, g)
644 call hchksum(kd_int,
"after set_diffusivity Kd_Int", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
645 call hchksum(kd_heat,
"after set_diffusivity Kd_heat", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
646 call hchksum(kd_salt,
"after set_diffusivity Kd_salt", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
650 call cpu_clock_begin(id_clock_kpp)
652 if (.not.cs%use_legacy_diabatic)
then
653 do k=1,nz+1 ;
do j=js,je ;
do i=is,ie
654 visc%Kv_shear(i,j,k) = visc%Kv_shear(i,j,k) + visc%Kv_slow(i,j,k)
655 enddo ;
enddo ;
enddo
663 call calculatebuoyancyflux2d(g, gv, us, fluxes, cs%optics, h, tv%T, tv%S, tv, &
664 cs%KPP_buoy_flux, cs%KPP_temp_flux, cs%KPP_salt_flux)
667 call kpp_compute_bld(cs%KPP_CSp, g, gv, us, h, tv%T, tv%S, u, v, tv%eqn_of_state, &
668 fluxes%ustar, cs%KPP_buoy_flux, waves=waves)
670 call kpp_calculate(cs%KPP_CSp, g, gv, us, h, fluxes%ustar, cs%KPP_buoy_flux, kd_heat, &
671 kd_salt, visc%Kv_shear, cs%KPP_NLTheat, cs%KPP_NLTscalar, waves=waves)
673 if (
associated(hml))
then
675 call kpp_get_bld(cs%KPP_CSp, hml(:,:), g)
677 call pass_var(hml, g%domain, halo=1)
679 if (
associated(visc%MLD)) visc%MLD(:,:) = hml(:,:)
682 if (cs%use_legacy_diabatic .and. .not.cs%KPPisPassive)
then
684 do k=1,nz+1 ;
do j=js,je ;
do i=is,ie
685 kd_int(i,j,k) = min( kd_salt(i,j,k), kd_heat(i,j,k) )
686 enddo ;
enddo ;
enddo
687 if (
associated(visc%Kd_extra_S))
then
689 do k=1,nz+1 ;
do j=js,je ;
do i=is,ie
690 visc%Kd_extra_S(i,j,k) = (kd_salt(i,j,k) - kd_int(i,j,k))
691 enddo ;
enddo ;
enddo
693 if (
associated(visc%Kd_extra_T))
then
695 do k=1,nz+1 ;
do j=js,je ;
do i=is,ie
696 visc%Kd_extra_T(i,j,k) = (kd_heat(i,j,k) - kd_int(i,j,k))
697 enddo ;
enddo ;
enddo
701 call cpu_clock_end(id_clock_kpp)
702 if (showcalltree)
call calltree_waypoint(
"done with KPP_calculate (diabatic)")
705 call mom_forcing_chksum(
"after KPP", fluxes, g, us, haloshift=0)
706 call mom_thermovar_chksum(
"after KPP", tv, g)
707 call hchksum(kd_heat,
"after KPP Kd_heat", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
708 call hchksum(kd_salt,
"after KPP Kd_salt", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
715 call cpu_clock_begin(id_clock_kpp)
717 call hchksum(cs%KPP_temp_flux,
"before KPP_applyNLT netHeat", g%HI, haloshift=0, scale=gv%H_to_m)
718 call hchksum(cs%KPP_salt_flux,
"before KPP_applyNLT netSalt", g%HI, haloshift=0, scale=gv%H_to_m)
719 call hchksum(cs%KPP_NLTheat,
"before KPP_applyNLT NLTheat", g%HI, haloshift=0)
720 call hchksum(cs%KPP_NLTscalar,
"before KPP_applyNLT NLTscalar", g%HI, haloshift=0)
724 call kpp_nonlocaltransport_temp(cs%KPP_CSp, g, gv, h, cs%KPP_NLTheat, cs%KPP_temp_flux, dt, tv%T, tv%C_p)
725 call kpp_nonlocaltransport_saln(cs%KPP_CSp, g, gv, h, cs%KPP_NLTscalar, cs%KPP_salt_flux, dt, tv%S)
726 call cpu_clock_end(id_clock_kpp)
727 if (showcalltree)
call calltree_waypoint(
"done with KPP_applyNonLocalTransport (diabatic)")
728 if (cs%debugConservation)
call mom_state_stats(
'KPP_applyNonLocalTransport', u, v, h, tv%T, tv%S, g)
732 call mom_forcing_chksum(
"after KPP_applyNLT ", fluxes, g, us, haloshift=0)
733 call mom_thermovar_chksum(
"after KPP_applyNLT ", tv, g)
739 if (
associated(visc%Kd_extra_T) .and.
associated(visc%Kd_extra_S) .and.
associated(tv%T) .and. &
740 (cs%use_legacy_diabatic .or. .not.cs%use_CVMix_ddiff))
then
742 call cpu_clock_begin(id_clock_differential_diff)
743 call differential_diffuse_t_s(h, tv, visc, dt_in_t, g, gv)
744 call cpu_clock_end(id_clock_differential_diff)
746 if (showcalltree)
call calltree_waypoint(
"done with differential_diffuse_T_S (diabatic)")
747 if (cs%debugConservation)
call mom_state_stats(
'differential_diffuse_T_S', u, v, h, tv%T, tv%S, g)
751 if (.not. cs%useKPP)
then
753 do k=2,nz ;
do j=js,je ;
do i=is,ie
754 kd_heat(i,j,k) = kd_heat(i,j,k) + visc%Kd_extra_T(i,j,k)
755 kd_salt(i,j,k) = kd_salt(i,j,k) + visc%Kd_extra_S(i,j,k)
756 enddo ;
enddo ;
enddo
762 if (cs%use_CVMix_conv)
then
763 call calculate_cvmix_conv(h, tv, g, gv, us, cs%CVMix_conv_csp, hml)
765 if (cs%use_legacy_diabatic)
then
767 do k=1,nz+1 ;
do j=js,je ;
do i=is,ie
768 kd_int(i,j,k) = kd_int(i,j,k) + cs%CVMix_conv_csp%kd_conv(i,j,k)
769 visc%Kv_slow(i,j,k) = visc%Kv_slow(i,j,k) + cs%CVMix_conv_csp%kv_conv(i,j,k)
770 enddo ;
enddo ;
enddo
773 do k=1,nz+1 ;
do j=js,je ;
do i=is,ie
774 kd_heat(i,j,k) = kd_heat(i,j,k) + cs%CVMix_conv_csp%kd_conv(i,j,k)
775 kd_salt(i,j,k) = kd_salt(i,j,k) + cs%CVMix_conv_csp%kd_conv(i,j,k)
777 visc%Kv_shear(i,j,k) = visc%Kv_shear(i,j,k) + cs%CVMix_conv_csp%kv_conv(i,j,k)
779 visc%Kv_slow(i,j,k) = visc%Kv_slow(i,j,k) + cs%CVMix_conv_csp%kv_conv(i,j,k)
781 enddo ;
enddo ;
enddo
786 if (cs%use_legacy_diabatic)
then
787 do j=js,je ;
do i=is,ie
791 do k=2,nz ;
do j=js,je ;
do i=is,ie
792 hval=1.0/(h_neglect + 0.5*(h(i,j,k-1) + h(i,j,k)))
793 ea_s(i,j,k) = (gv%Z_to_H**2) * dt_in_t * hval * kd_int(i,j,k)
794 eb_s(i,j,k-1) = ea_s(i,j,k)
795 ea_t(i,j,k-1) = ea_s(i,j,k-1) ; eb_t(i,j,k-1) = eb_s(i,j,k-1)
796 enddo ;
enddo ;
enddo
797 do j=js,je ;
do i=is,ie
799 ea_t(i,j,nz) = ea_s(i,j,nz) ; eb_t(i,j,nz) = eb_s(i,j,nz)
801 if (showcalltree)
call calltree_waypoint(
"done setting ea,eb from Kd_int (diabatic)")
805 call mom_forcing_chksum(
"after calc_entrain ", fluxes, g, us, haloshift=0)
806 call mom_thermovar_chksum(
"after calc_entrain ", tv, g)
808 call hchksum(ea_s,
"after calc_entrain ea_s", g%HI, haloshift=0, scale=gv%H_to_m)
809 call hchksum(eb_s,
"after calc_entrain eb_s", g%HI, haloshift=0, scale=gv%H_to_m)
813 if (cs%boundary_forcing_tendency_diag)
then
814 do k=1,nz ;
do j=js,je ;
do i=is,ie
815 h_diag(i,j,k) = h(i,j,k)
816 temp_diag(i,j,k) = tv%T(i,j,k)
817 saln_diag(i,j,k) = tv%S(i,j,k)
818 enddo ;
enddo ;
enddo
822 call cpu_clock_begin(id_clock_remap)
825 do k=1,nz ;
do j=js,je ;
do i=is,ie
826 h_prebound(i,j,k) = h(i,j,k)
827 enddo ;
enddo ;
enddo
828 if (cs%use_energetic_PBL)
then
830 skinbuoyflux(:,:) = 0.0
831 call applyboundaryfluxesinout(cs%diabatic_aux_CSp, g, gv, us, dt, fluxes, cs%optics, &
832 optics_nbands(cs%optics), h, tv, cs%aggregate_FW_forcing, cs%evap_CFL_limit, &
833 cs%minimum_forcing_depth, ctke, dsv_dt, dsv_ds, skinbuoyflux=skinbuoyflux)
836 call hchksum(ea_t,
"after applyBoundaryFluxes ea_t",g%HI,haloshift=0, scale=gv%H_to_m)
837 call hchksum(eb_t,
"after applyBoundaryFluxes eb_t",g%HI,haloshift=0, scale=gv%H_to_m)
838 call hchksum(ea_s,
"after applyBoundaryFluxes ea_s",g%HI,haloshift=0, scale=gv%H_to_m)
839 call hchksum(eb_s,
"after applyBoundaryFluxes eb_s",g%HI,haloshift=0, scale=gv%H_to_m)
840 call hchksum(ctke,
"after applyBoundaryFluxes cTKE",g%HI,haloshift=0)
841 call hchksum(dsv_dt,
"after applyBoundaryFluxes dSV_dT",g%HI,haloshift=0)
842 call hchksum(dsv_ds,
"after applyBoundaryFluxes dSV_dS",g%HI,haloshift=0)
845 call find_uv_at_h(u, v, h, u_h, v_h, g, gv)
846 call energetic_pbl(h, u_h, v_h, tv, fluxes, dt_in_t, kd_epbl, g, gv, us, &
847 cs%energetic_PBL_CSp, dsv_dt, dsv_ds, ctke, skinbuoyflux, waves=waves)
849 if (
associated(hml))
then
850 call energetic_pbl_get_mld(cs%energetic_PBL_CSp, hml(:,:), g, us)
851 call pass_var(hml, g%domain, halo=1)
853 if (
associated(visc%MLD)) visc%MLD(:,:) = hml(:,:)
854 elseif (
associated(visc%MLD))
then
855 call energetic_pbl_get_mld(cs%energetic_PBL_CSp, visc%MLD, g, us)
856 call pass_var(visc%MLD, g%domain, halo=1)
860 do k=2,nz ;
do j=js,je ;
do i=is,ie
862 if (cs%ePBL_is_additive)
then
863 kd_add_here = kd_epbl(i,j,k)
864 visc%Kv_shear(i,j,k) = visc%Kv_shear(i,j,k) + kd_epbl(i,j,k)
866 kd_add_here = max(kd_epbl(i,j,k) - visc%Kd_shear(i,j,k), 0.0)
867 visc%Kv_shear(i,j,k) = max(visc%Kv_shear(i,j,k), kd_epbl(i,j,k))
870 if (cs%use_legacy_diabatic)
then
871 ent_int = kd_add_here * (gv%Z_to_H**2 * dt_in_t) / &
872 (0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect)
873 eb_s(i,j,k-1) = eb_s(i,j,k-1) + ent_int
874 ea_s(i,j,k) = ea_s(i,j,k) + ent_int
875 kd_int(i,j,k) = kd_int(i,j,k) + kd_add_here
878 kd_heat(i,j,k) = kd_heat(i,j,k) + kd_int(i,j,k)
879 kd_salt(i,j,k) = kd_salt(i,j,k) + kd_int(i,j,k)
881 kd_heat(i,j,k) = kd_heat(i,j,k) + kd_add_here
882 kd_salt(i,j,k) = kd_salt(i,j,k) + kd_add_here
885 enddo ;
enddo ;
enddo
888 call hchksum(ea_t,
"after ePBL ea_t",g%HI,haloshift=0, scale=gv%H_to_m)
889 call hchksum(eb_t,
"after ePBL eb_t",g%HI,haloshift=0, scale=gv%H_to_m)
890 call hchksum(ea_s,
"after ePBL ea_s",g%HI,haloshift=0, scale=gv%H_to_m)
891 call hchksum(eb_s,
"after ePBL eb_s",g%HI,haloshift=0, scale=gv%H_to_m)
892 call hchksum(kd_epbl,
"after ePBL Kd_ePBL", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
896 call applyboundaryfluxesinout(cs%diabatic_aux_CSp, g, gv, us, dt, fluxes, cs%optics, &
897 optics_nbands(cs%optics), h, tv, cs%aggregate_FW_forcing, &
898 cs%evap_CFL_limit, cs%minimum_forcing_depth)
905 if (cs%boundary_forcing_tendency_diag)
then
906 call diagnose_boundary_forcing_tendency(tv, h, temp_diag, saln_diag, h_diag, dt, g, gv, cs)
907 if (cs%id_boundary_forcing_h > 0)
call post_data(cs%id_boundary_forcing_h, h, cs%diag, alt_h = h_diag)
910 call diag_update_remap_grids(cs%diag)
911 call cpu_clock_end(id_clock_remap)
913 call mom_forcing_chksum(
"after applyBoundaryFluxes ", fluxes, g, us, haloshift=0)
914 call mom_thermovar_chksum(
"after applyBoundaryFluxes ", tv, g)
915 call mom_state_chksum(
"after applyBoundaryFluxes ", u, v, h, g, gv, haloshift=0)
917 if (showcalltree)
call calltree_waypoint(
"done with applyBoundaryFluxes (diabatic)")
918 if (cs%debugConservation)
call mom_state_stats(
'applyBoundaryFluxes', u, v, h, tv%T, tv%S, g)
926 if (cs%use_legacy_diabatic)
then
930 hold(i,j,1) = h(i,j,1)
932 hold(i,j,nz) = h(i,j,nz)
934 if (h(i,j,1) <= 0.0) h(i,j,1) = gv%Angstrom_H
935 if (h(i,j,nz) <= 0.0) h(i,j,nz) = gv%Angstrom_H
937 do k=2,nz-1 ;
do i=is,ie
938 hold(i,j,k) = h(i,j,k)
941 if (h(i,j,k) <= 0.0) h(i,j,k) = gv%Angstrom_H
945 call diag_update_remap_grids(cs%diag)
950 call mom_forcing_chksum(
"after negative check ", fluxes, g, us, haloshift=0)
951 call mom_thermovar_chksum(
"after negative check ", tv, g)
953 if (showcalltree)
call calltree_waypoint(
"done with h=ea-eb (diabatic)")
954 if (cs%debugConservation)
call mom_state_stats(
'h=ea-eb', u, v, h, tv%T, tv%S, g)
957 if (
associated(tv%T))
then
960 call hchksum(ea_t,
"before triDiagTS ea_t ",g%HI,haloshift=0, scale=gv%H_to_m)
961 call hchksum(eb_t,
"before triDiagTS eb_t ",g%HI,haloshift=0, scale=gv%H_to_m)
962 call hchksum(ea_s,
"before triDiagTS ea_s ",g%HI,haloshift=0, scale=gv%H_to_m)
963 call hchksum(eb_s,
"before triDiagTS eb_s ",g%HI,haloshift=0, scale=gv%H_to_m)
966 call cpu_clock_begin(id_clock_tridiag)
973 if (
associated(tv%S) .and.
associated(tv%salt_deficit)) &
974 call adjust_salt(h, tv, g, gv, cs%diabatic_aux_CSp)
976 if (cs%diabatic_diff_tendency_diag)
then
977 do k=1,nz ;
do j=js,je ;
do i=is,ie
978 temp_diag(i,j,k) = tv%T(i,j,k)
979 saln_diag(i,j,k) = tv%S(i,j,k)
980 enddo ;
enddo ;
enddo
983 if (cs%use_legacy_diabatic)
then
985 do k=1,nz ;
do j=js,je ;
do i=is,ie
986 ea_t(i,j,k) = ea_s(i,j,k) ; eb_t(i,j,k) = eb_s(i,j,k)
987 enddo ;
enddo ;
enddo
988 if (cs%tracer_tridiag)
then
989 call tracer_vertdiff(hold, ea_t, eb_t, dt, tv%T, g, gv)
990 call tracer_vertdiff(hold, ea_s, eb_s, dt, tv%S, g, gv)
992 call tridiagts(g, gv, is, ie, js, je, hold, ea_s, eb_s, tv%T, tv%S)
999 if (cs%diabatic_diff_tendency_diag)
then
1000 call diagnose_diabatic_diff_tendency(tv, hold, temp_diag, saln_diag, dt, g, gv, cs)
1001 if (cs%id_diabatic_diff_h > 0)
call post_data(cs%id_diabatic_diff_h, hold, cs%diag, alt_h = hold)
1006 do j=js,je ;
do i=is,ie
1007 ea_t(i,j,1) = 0.; ea_s(i,j,1) = 0.
1011 do k=2,nz ;
do j=js,je ;
do i=is,ie
1012 hval = 1.0 / (h_neglect + 0.5*(h(i,j,k-1) + h(i,j,k)))
1013 ea_t(i,j,k) = (gv%Z_to_H**2) * dt_in_t * hval * kd_heat(i,j,k)
1014 eb_t(i,j,k-1) = ea_t(i,j,k)
1015 ea_s(i,j,k) = (gv%Z_to_H**2) * dt_in_t * hval * kd_salt(i,j,k)
1016 eb_s(i,j,k-1) = ea_s(i,j,k)
1017 enddo ;
enddo ;
enddo
1018 do j=js,je ;
do i=is,ie
1019 eb_t(i,j,nz) = 0. ; eb_s(i,j,nz) = 0.
1021 if (showcalltree)
call calltree_waypoint(
"done setting ea_t,ea_s,eb_t,eb_s from Kd_heat" //&
1022 "and Kd_salt (diabatic)")
1025 call tracer_vertdiff(h, ea_t, eb_t, dt, tv%T, g, gv)
1026 call tracer_vertdiff(h, ea_s, eb_s, dt, tv%S, g, gv)
1029 if (cs%diabatic_diff_tendency_diag) &
1030 call diagnose_diabatic_diff_tendency(tv, h, temp_diag, saln_diag, dt, g, gv, cs)
1033 call cpu_clock_end(id_clock_tridiag)
1035 if (showcalltree)
call calltree_waypoint(
"done with triDiagTS (diabatic)")
1039 if (cs%debugConservation)
call mom_state_stats(
'triDiagTS', u, v, h, tv%T, tv%S, g)
1043 call mom_thermovar_chksum(
"after mixed layer ", tv, g)
1048 if (cs%id_dudt_dia > 0 .or. cs%id_dvdt_dia > 0) &
1050 call pass_var(h, g%domain, to_west+to_south+omit_corners, halo=1)
1051 call diag_update_remap_grids(cs%diag)
1055 if ((cs%id_Tdif > 0) .or. (cs%id_Tadv > 0))
then
1056 do j=js,je ;
do i=is,ie
1057 tdif_flx(i,j,1) = 0.0 ; tdif_flx(i,j,nz+1) = 0.0
1058 tadv_flx(i,j,1) = 0.0 ; tadv_flx(i,j,nz+1) = 0.0
1061 do k=2,nz ;
do j=js,je ;
do i=is,ie
1062 tdif_flx(i,j,k) = (idt * 0.5*(ea_t(i,j,k) + eb_t(i,j,k-1))) * &
1063 (tv%T(i,j,k-1) - tv%T(i,j,k))
1064 tadv_flx(i,j,k) = (idt * (ea_t(i,j,k) - eb_t(i,j,k-1))) * &
1065 0.5*(tv%T(i,j,k-1) + tv%T(i,j,k))
1066 enddo ;
enddo ;
enddo
1068 if ((cs%id_Sdif > 0) .or. (cs%id_Sadv > 0))
then
1069 do j=js,je ;
do i=is,ie
1070 sdif_flx(i,j,1) = 0.0 ; sdif_flx(i,j,nz+1) = 0.0
1071 sadv_flx(i,j,1) = 0.0 ; sadv_flx(i,j,nz+1) = 0.0
1074 do k=2,nz ;
do j=js,je ;
do i=is,ie
1075 sdif_flx(i,j,k) = (idt * 0.5*(ea_s(i,j,k) + eb_s(i,j,k-1))) * &
1076 (tv%S(i,j,k-1) - tv%S(i,j,k))
1077 sadv_flx(i,j,k) = (idt * (ea_s(i,j,k) - eb_s(i,j,k-1))) * &
1078 0.5*(tv%S(i,j,k-1) + tv%S(i,j,k))
1079 enddo ;
enddo ;
enddo
1083 call cpu_clock_begin(id_clock_tracers)
1085 if (cs%mix_boundary_tracers)
then
1086 tr_ea_bbl = gv%Z_to_H * sqrt(dt_in_t*cs%Kd_BBL_tr)
1090 ebtr(i,j,nz) = eb_s(i,j,nz)
1092 in_boundary(i) = (g%mask2dT(i,j) > 0.0)
1094 do k=nz,2,-1 ;
do i=is,ie
1095 if (in_boundary(i))
then
1096 htot(i) = htot(i) + h(i,j,k)
1105 add_ent = ((dt_in_t * cs%Kd_min_tr) * gv%Z_to_H**2) * &
1106 ((h(i,j,k-1)+h(i,j,k)+h_neglect) / &
1107 (h(i,j,k-1)*h(i,j,k)+h_neglect2)) - &
1108 0.5*(ea_s(i,j,k) + eb_s(i,j,k-1))
1109 if (htot(i) < tr_ea_bbl)
then
1110 add_ent = max(0.0, add_ent, &
1111 (tr_ea_bbl - htot(i)) - min(ea_s(i,j,k),eb_s(i,j,k-1)))
1112 elseif (add_ent < 0.0)
then
1113 add_ent = 0.0 ; in_boundary(i) = .false.
1116 ebtr(i,j,k-1) = eb_s(i,j,k-1) + add_ent
1117 eatr(i,j,k) = ea_s(i,j,k) + add_ent
1119 ebtr(i,j,k-1) = eb_s(i,j,k-1) ; eatr(i,j,k) = ea_s(i,j,k)
1122 if (
associated(visc%Kd_extra_S))
then ;
if (visc%Kd_extra_S(i,j,k) > 0.0)
then
1123 if (cs%use_legacy_diabatic)
then
1124 add_ent = ((dt_in_t * visc%Kd_extra_S(i,j,k)) * gv%Z_to_H**2) / &
1125 (0.25 * ((h(i,j,k-1) + h(i,j,k)) + (hold(i,j,k-1) + hold(i,j,k))) + h_neglect)
1127 add_ent = ((dt_in_t * visc%Kd_extra_S(i,j,k)) * gv%Z_to_H**2) / &
1128 (0.5 * (h(i,j,k-1) + h(i,j,k)) + &
1131 ebtr(i,j,k-1) = ebtr(i,j,k-1) + add_ent
1132 eatr(i,j,k) = eatr(i,j,k) + add_ent
1135 do i=is,ie ; eatr(i,j,1) = ea_s(i,j,1) ;
enddo
1141 call call_tracer_column_fns(h_prebound, h, ea_s, eb_s, fluxes, hml, dt, g, gv, tv, &
1142 cs%optics, cs%tracer_flow_CSp, cs%debug, &
1143 evap_cfl_limit = cs%evap_CFL_limit, &
1144 minimum_forcing_depth = cs%minimum_forcing_depth)
1146 elseif (
associated(visc%Kd_extra_S))
then
1148 do j=js,je ;
do i=is,ie
1149 ebtr(i,j,nz) = eb_s(i,j,nz) ; eatr(i,j,1) = ea_s(i,j,1)
1152 do k=nz,2,-1 ;
do j=js,je ;
do i=is,ie
1153 if (visc%Kd_extra_S(i,j,k) > 0.0)
then
1154 if (cs%use_legacy_diabatic)
then
1155 add_ent = ((dt_in_t * visc%Kd_extra_S(i,j,k)) * gv%Z_to_H**2) / &
1156 (0.25 * ((h(i,j,k-1) + h(i,j,k)) + (hold(i,j,k-1) + hold(i,j,k))) + h_neglect)
1158 add_ent = ((dt_in_t * visc%Kd_extra_S(i,j,k)) * gv%Z_to_H**2) / &
1159 (0.5 * (h(i,j,k-1) + h(i,j,k)) + &
1165 ebtr(i,j,k-1) = eb_s(i,j,k-1) + add_ent
1166 eatr(i,j,k) = ea_s(i,j,k) + add_ent
1167 enddo ;
enddo ;
enddo
1170 call call_tracer_column_fns(h_prebound, h, eatr, ebtr, fluxes, hml, dt, g, gv, tv, &
1171 cs%optics, cs%tracer_flow_CSp, cs%debug,&
1172 evap_cfl_limit = cs%evap_CFL_limit, &
1173 minimum_forcing_depth = cs%minimum_forcing_depth)
1176 call call_tracer_column_fns(h_prebound, h, eatr, ebtr, fluxes, hml, dt, g, gv, tv, &
1177 cs%optics, cs%tracer_flow_CSp, cs%debug, &
1178 evap_cfl_limit = cs%evap_CFL_limit, &
1179 minimum_forcing_depth = cs%minimum_forcing_depth)
1182 call cpu_clock_end(id_clock_tracers)
1185 if (cs%use_sponge)
then
1186 call cpu_clock_begin(id_clock_sponge)
1187 if (
associated(cs%ALE_sponge_CSp))
then
1188 call apply_ale_sponge(h, dt, g, gv, us, cs%ALE_sponge_CSp, cs%Time)
1191 call cpu_clock_end(id_clock_sponge)
1194 call mom_thermovar_chksum(
"apply_sponge ", tv, g)
1198 call disable_averaging(cs%diag)
1200 call enable_averaging(dt, time_end, cs%diag)
1202 if (cs%id_Kd_interface > 0)
call post_data(cs%id_Kd_interface, kd_int, cs%diag)
1203 if (cs%id_Kd_heat > 0)
call post_data(cs%id_Kd_heat, kd_heat, cs%diag)
1204 if (cs%id_Kd_salt > 0)
call post_data(cs%id_Kd_salt, kd_salt, cs%diag)
1205 if (cs%id_Kd_ePBL > 0)
call post_data(cs%id_Kd_ePBL, kd_epbl, cs%diag)
1207 if (cs%id_ea > 0)
call post_data(cs%id_ea, ea_s, cs%diag)
1208 if (cs%id_eb > 0)
call post_data(cs%id_eb, eb_s, cs%diag)
1209 if (cs%id_ea_t > 0)
call post_data(cs%id_ea_t, ea_t, cs%diag)
1210 if (cs%id_eb_t > 0)
call post_data(cs%id_eb_t, eb_t, cs%diag)
1211 if (cs%id_ea_s > 0)
call post_data(cs%id_ea_s, ea_s, cs%diag)
1212 if (cs%id_eb_s > 0)
call post_data(cs%id_eb_s, eb_s, cs%diag)
1214 if (cs%id_dudt_dia > 0)
call post_data(cs%id_dudt_dia, adp%du_dt_dia, cs%diag)
1215 if (cs%id_dvdt_dia > 0)
call post_data(cs%id_dvdt_dia, adp%dv_dt_dia, cs%diag)
1216 if (cs%id_wd > 0)
call post_data(cs%id_wd, cdp%diapyc_vel, cs%diag)
1218 if (cs%id_Tdif > 0)
call post_data(cs%id_Tdif, tdif_flx, cs%diag)
1219 if (cs%id_Tadv > 0)
call post_data(cs%id_Tadv, tadv_flx, cs%diag)
1220 if (cs%id_Sdif > 0)
call post_data(cs%id_Sdif, sdif_flx, cs%diag)
1221 if (cs%id_Sadv > 0)
call post_data(cs%id_Sadv, sadv_flx, cs%diag)
1223 call disable_averaging(cs%diag)
1225 if (showcalltree)
call calltree_leave(
"diabatic_ALE_legacy()")
1227 end subroutine diabatic_ale_legacy
1232 subroutine diabatic_ale(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, &
1233 G, GV, US, CS, Waves)
1237 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(inout) :: u
1238 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(inout) :: v
1239 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(inout) :: h
1242 real,
dimension(:,:),
pointer :: Hml
1243 type(
forcing),
intent(inout) :: fluxes
1250 real,
intent(in) :: dt
1251 type(time_type),
intent(in) :: Time_end
1256 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
1276 real,
dimension(SZI_(G),SZJ_(G)) :: &
1279 real,
dimension(SZI_(G),SZJ_(G),G%ke) :: h_diag
1280 real,
dimension(SZI_(G),SZJ_(G),G%ke) :: temp_diag
1281 real,
dimension(SZI_(G),SZJ_(G),G%ke) :: saln_diag
1282 real,
dimension(SZI_(G),SZJ_(G)) :: tendency_2d
1286 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
1291 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1),
target :: &
1301 logical :: in_boundary(SZI_(G))
1321 real :: htot(SZIB_(G))
1322 real :: b1(SZIB_(G)), d1(SZIB_(G))
1323 real :: c1(SZIB_(G),SZK_(G))
1330 logical :: showCallTree
1331 integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb, m, halo
1336 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1337 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
1338 nkmb = gv%nk_rho_varies
1339 h_neglect = gv%H_subroundoff ; h_neglect2 = h_neglect*h_neglect
1340 kd_heat(:,:,:) = 0.0 ; kd_salt(:,:,:) = 0.0
1342 showcalltree = calltree_showquery()
1343 if (showcalltree)
call calltree_enter(
"diabatic_ALE(), MOM_diabatic_driver.F90")
1345 if (.not. (cs%useALEalgorithm))
call mom_error(fatal,
"MOM_diabatic_driver: "// &
1346 "The ALE algorithm must be enabled when using MOM_diabatic_driver.")
1348 dt_in_t = dt * us%s_to_T
1351 call enable_averaging(dt, time_end, cs%diag)
1353 if (cs%use_geothermal)
then
1354 halo = cs%halo_TS_diff
1356 do k=1,nz ;
do j=js-halo,je+halo ;
do i=is-halo,ie+halo
1357 h_orig(i,j,k) = h(i,j,k) ; eatr(i,j,k) = 0.0 ; ebtr(i,j,k) = 0.0
1358 enddo ;
enddo ;
enddo
1361 if (cs%use_geothermal)
then
1362 call cpu_clock_begin(id_clock_geothermal)
1363 call geothermal(h, tv, dt, eatr, ebtr, g, gv, cs%geothermal_CSp, halo=cs%halo_TS_diff)
1364 call cpu_clock_end(id_clock_geothermal)
1365 if (showcalltree)
call calltree_waypoint(
"geothermal (diabatic)")
1366 if (cs%debugConservation)
call mom_state_stats(
'geothermal', u, v, h, tv%T, tv%S, g)
1371 call diag_update_remap_grids(cs%diag)
1376 if (
associated(cs%optics)) &
1377 call set_pen_shortwave(cs%optics, fluxes, g, gv, cs%diabatic_aux_CSp, cs%opacity_CSp, cs%tracer_flow_CSp)
1379 if (cs%debug)
call mom_state_chksum(
"before find_uv_at_h", u, v, h, g, gv, haloshift=0)
1381 if (cs%use_kappa_shear .or. cs%use_CVMix_shear)
then
1382 if (cs%use_geothermal)
then
1383 call find_uv_at_h(u, v, h_orig, u_h, v_h, g, gv, eatr, ebtr)
1385 call hchksum(eatr,
"after find_uv_at_h eatr",g%HI, scale=gv%H_to_m)
1386 call hchksum(ebtr,
"after find_uv_at_h ebtr",g%HI, scale=gv%H_to_m)
1389 call find_uv_at_h(u, v, h, u_h, v_h, g, gv)
1391 if (showcalltree)
call calltree_waypoint(
"done with find_uv_at_h (diabatic)")
1394 call cpu_clock_begin(id_clock_set_diffusivity)
1397 call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, cs%optics, visc, dt_in_t, g, gv, us, &
1398 cs%set_diff_CSp, kd_lay, kd_int)
1399 call cpu_clock_end(id_clock_set_diffusivity)
1400 if (showcalltree)
call calltree_waypoint(
"done with set_diffusivity (diabatic)")
1403 call mom_state_chksum(
"after set_diffusivity ", u, v, h, g, gv, haloshift=0)
1404 call mom_forcing_chksum(
"after set_diffusivity ", fluxes, g, us, haloshift=0)
1405 call mom_thermovar_chksum(
"after set_diffusivity ", tv, g)
1406 call hchksum(kd_int,
"after set_diffusivity Kd_Int", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
1412 do k=1,nz+1 ;
do j=js,je ;
do i=is,ie
1413 kd_salt(i,j,k) = kd_int(i,j,k)
1414 kd_heat(i,j,k) = kd_int(i,j,k)
1415 enddo ;
enddo ;
enddo
1417 if (
associated(visc%Kd_extra_S))
then
1419 do k=1,nz+1 ;
do j=js,je ;
do i=is,ie
1420 kd_salt(i,j,k) = kd_salt(i,j,k) + visc%Kd_extra_S(i,j,k)
1421 enddo ;
enddo ;
enddo
1423 if (
associated(visc%Kd_extra_T))
then
1425 do k=1,nz+1 ;
do j=js,je ;
do i=is,ie
1426 kd_heat(i,j,k) = kd_heat(i,j,k) + visc%Kd_extra_T(i,j,k)
1427 enddo ;
enddo ;
enddo
1431 call hchksum(kd_heat,
"after set_diffusivity Kd_heat", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
1432 call hchksum(kd_salt,
"after set_diffusivity Kd_salt", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
1436 call cpu_clock_begin(id_clock_kpp)
1438 do k=1,nz+1 ;
do j=js,je ;
do i=is,ie
1439 visc%Kv_shear(i,j,k) = visc%Kv_shear(i,j,k) + visc%Kv_slow(i,j,k)
1440 enddo ;
enddo ;
enddo
1447 call calculatebuoyancyflux2d(g, gv, us, fluxes, cs%optics, h, tv%T, tv%S, tv, &
1448 cs%KPP_buoy_flux, cs%KPP_temp_flux, cs%KPP_salt_flux)
1451 call kpp_compute_bld(cs%KPP_CSp, g, gv, us, h, tv%T, tv%S, u, v, tv%eqn_of_state, &
1452 fluxes%ustar, cs%KPP_buoy_flux, waves=waves)
1454 call kpp_calculate(cs%KPP_CSp, g, gv, us, h, fluxes%ustar, cs%KPP_buoy_flux, kd_heat, &
1455 kd_salt, visc%Kv_shear, cs%KPP_NLTheat, cs%KPP_NLTscalar, waves=waves)
1457 if (
associated(hml))
then
1459 call kpp_get_bld(cs%KPP_CSp, hml(:,:), g)
1461 call pass_var(hml, g%domain, halo=1)
1463 if (
associated(visc%MLD)) visc%MLD(:,:) = hml(:,:)
1466 call cpu_clock_end(id_clock_kpp)
1467 if (showcalltree)
call calltree_waypoint(
"done with KPP_calculate (diabatic)")
1470 call mom_forcing_chksum(
"after KPP", fluxes, g, us, haloshift=0)
1471 call mom_thermovar_chksum(
"after KPP", tv, g)
1472 call hchksum(kd_heat,
"after KPP Kd_heat", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
1473 call hchksum(kd_salt,
"after KPP Kd_salt", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
1480 call cpu_clock_begin(id_clock_kpp)
1482 call hchksum(cs%KPP_temp_flux,
"before KPP_applyNLT netHeat", g%HI, haloshift=0, scale=gv%H_to_m)
1483 call hchksum(cs%KPP_salt_flux,
"before KPP_applyNLT netSalt", g%HI, haloshift=0, scale=gv%H_to_m)
1484 call hchksum(cs%KPP_NLTheat,
"before KPP_applyNLT NLTheat", g%HI, haloshift=0)
1485 call hchksum(cs%KPP_NLTscalar,
"before KPP_applyNLT NLTscalar", g%HI, haloshift=0)
1489 call kpp_nonlocaltransport_temp(cs%KPP_CSp, g, gv, h, cs%KPP_NLTheat, cs%KPP_temp_flux, dt, tv%T, tv%C_p)
1490 call kpp_nonlocaltransport_saln(cs%KPP_CSp, g, gv, h, cs%KPP_NLTscalar, cs%KPP_salt_flux, dt, tv%S)
1491 call cpu_clock_end(id_clock_kpp)
1492 if (showcalltree)
call calltree_waypoint(
"done with KPP_applyNonLocalTransport (diabatic)")
1493 if (cs%debugConservation)
call mom_state_stats(
'KPP_applyNonLocalTransport', u, v, h, tv%T, tv%S, g)
1497 call mom_forcing_chksum(
"after KPP_applyNLT ", fluxes, g, us, haloshift=0)
1498 call mom_thermovar_chksum(
"after KPP_applyNLT ", tv, g)
1504 if (
associated(visc%Kd_extra_T) .and.
associated(visc%Kd_extra_S) .and.
associated(tv%T) .and. &
1505 (.not.cs%use_CVMix_ddiff))
then
1507 call cpu_clock_begin(id_clock_differential_diff)
1508 call differential_diffuse_t_s(h, tv, visc, dt_in_t, g, gv)
1509 call cpu_clock_end(id_clock_differential_diff)
1511 if (showcalltree)
call calltree_waypoint(
"done with differential_diffuse_T_S (diabatic)")
1512 if (cs%debugConservation)
call mom_state_stats(
'differential_diffuse_T_S', u, v, h, tv%T, tv%S, g)
1516 if (.not. cs%useKPP)
then
1518 do k=2,nz ;
do j=js,je ;
do i=is,ie
1519 kd_heat(i,j,k) = kd_heat(i,j,k) + visc%Kd_extra_T(i,j,k)
1520 kd_salt(i,j,k) = kd_salt(i,j,k) + visc%Kd_extra_S(i,j,k)
1521 enddo ;
enddo ;
enddo
1527 if (cs%use_CVMix_conv)
then
1528 call calculate_cvmix_conv(h, tv, g, gv, us, cs%CVMix_conv_csp, hml)
1531 do k=1,nz+1 ;
do j=js,je ;
do i=is,ie
1532 kd_heat(i,j,k) = kd_heat(i,j,k) + cs%CVMix_conv_csp%kd_conv(i,j,k)
1533 kd_salt(i,j,k) = kd_salt(i,j,k) + cs%CVMix_conv_csp%kd_conv(i,j,k)
1535 visc%Kv_shear(i,j,k) = visc%Kv_shear(i,j,k) + cs%CVMix_conv_csp%kv_conv(i,j,k)
1537 visc%Kv_slow(i,j,k) = visc%Kv_slow(i,j,k) + cs%CVMix_conv_csp%kv_conv(i,j,k)
1539 enddo ;
enddo ;
enddo
1543 if (cs%boundary_forcing_tendency_diag)
then
1544 do k=1,nz ;
do j=js,je ;
do i=is,ie
1545 h_diag(i,j,k) = h(i,j,k)
1546 temp_diag(i,j,k) = tv%T(i,j,k)
1547 saln_diag(i,j,k) = tv%S(i,j,k)
1548 enddo ;
enddo ;
enddo
1552 call cpu_clock_begin(id_clock_remap)
1555 do k=1,nz ;
do j=js,je ;
do i=is,ie
1556 h_prebound(i,j,k) = h(i,j,k)
1557 enddo ;
enddo ;
enddo
1558 if (cs%use_energetic_PBL)
then
1560 skinbuoyflux(:,:) = 0.0
1561 call applyboundaryfluxesinout(cs%diabatic_aux_CSp, g, gv, us, dt, fluxes, cs%optics, &
1562 optics_nbands(cs%optics), h, tv, cs%aggregate_FW_forcing, cs%evap_CFL_limit, &
1563 cs%minimum_forcing_depth, ctke, dsv_dt, dsv_ds, skinbuoyflux=skinbuoyflux)
1566 call hchksum(ea_t,
"after applyBoundaryFluxes ea_t",g%HI,haloshift=0, scale=gv%H_to_m)
1567 call hchksum(eb_t,
"after applyBoundaryFluxes eb_t",g%HI,haloshift=0, scale=gv%H_to_m)
1568 call hchksum(ea_s,
"after applyBoundaryFluxes ea_s",g%HI,haloshift=0, scale=gv%H_to_m)
1569 call hchksum(eb_s,
"after applyBoundaryFluxes eb_s",g%HI,haloshift=0, scale=gv%H_to_m)
1570 call hchksum(ctke,
"after applyBoundaryFluxes cTKE",g%HI,haloshift=0)
1571 call hchksum(dsv_dt,
"after applyBoundaryFluxes dSV_dT",g%HI,haloshift=0)
1572 call hchksum(dsv_ds,
"after applyBoundaryFluxes dSV_dS",g%HI,haloshift=0)
1575 call find_uv_at_h(u, v, h, u_h, v_h, g, gv)
1576 call energetic_pbl(h, u_h, v_h, tv, fluxes, dt_in_t, kd_epbl, g, gv, us, &
1577 cs%energetic_PBL_CSp, dsv_dt, dsv_ds, ctke, skinbuoyflux, waves=waves)
1579 if (
associated(hml))
then
1580 call energetic_pbl_get_mld(cs%energetic_PBL_CSp, hml(:,:), g, us)
1581 call pass_var(hml, g%domain, halo=1)
1583 if (
associated(visc%MLD)) visc%MLD(:,:) = hml(:,:)
1584 elseif (
associated(visc%MLD))
then
1585 call energetic_pbl_get_mld(cs%energetic_PBL_CSp, visc%MLD, g, us)
1586 call pass_var(visc%MLD, g%domain, halo=1)
1590 do k=2,nz ;
do j=js,je ;
do i=is,ie
1592 if (cs%ePBL_is_additive)
then
1593 kd_add_here = kd_epbl(i,j,k)
1594 visc%Kv_shear(i,j,k) = visc%Kv_shear(i,j,k) + kd_epbl(i,j,k)
1596 kd_add_here = max(kd_epbl(i,j,k) - visc%Kd_shear(i,j,k), 0.0)
1597 visc%Kv_shear(i,j,k) = max(visc%Kv_shear(i,j,k), kd_epbl(i,j,k))
1600 kd_heat(i,j,k) = kd_heat(i,j,k) + kd_add_here
1601 kd_salt(i,j,k) = kd_salt(i,j,k) + kd_add_here
1603 enddo ;
enddo ;
enddo
1606 call hchksum(ea_t,
"after ePBL ea_t",g%HI,haloshift=0, scale=gv%H_to_m)
1607 call hchksum(eb_t,
"after ePBL eb_t",g%HI,haloshift=0, scale=gv%H_to_m)
1608 call hchksum(ea_s,
"after ePBL ea_s",g%HI,haloshift=0, scale=gv%H_to_m)
1609 call hchksum(eb_s,
"after ePBL eb_s",g%HI,haloshift=0, scale=gv%H_to_m)
1610 call hchksum(kd_epbl,
"after ePBL Kd_ePBL", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
1614 call applyboundaryfluxesinout(cs%diabatic_aux_CSp, g, gv, us, dt, fluxes, cs%optics, &
1615 optics_nbands(cs%optics), h, tv, cs%aggregate_FW_forcing, &
1616 cs%evap_CFL_limit, cs%minimum_forcing_depth)
1623 if (cs%boundary_forcing_tendency_diag)
then
1624 call diagnose_boundary_forcing_tendency(tv, h, temp_diag, saln_diag, h_diag, dt, g, gv, cs)
1625 if (cs%id_boundary_forcing_h > 0)
call post_data(cs%id_boundary_forcing_h, h, cs%diag, alt_h = h_diag)
1628 call diag_update_remap_grids(cs%diag)
1629 call cpu_clock_end(id_clock_remap)
1631 call mom_forcing_chksum(
"after applyBoundaryFluxes ", fluxes, g, us, haloshift=0)
1632 call mom_thermovar_chksum(
"after applyBoundaryFluxes ", tv, g)
1633 call mom_state_chksum(
"after applyBoundaryFluxes ", u, v, h, g, gv, haloshift=0)
1635 if (showcalltree)
call calltree_waypoint(
"done with applyBoundaryFluxes (diabatic)")
1636 if (cs%debugConservation)
call mom_state_stats(
'applyBoundaryFluxes', u, v, h, tv%T, tv%S, g)
1638 if (showcalltree)
call calltree_waypoint(
"done with h=ea-eb (diabatic)")
1639 if (cs%debugConservation)
call mom_state_stats(
'h=ea-eb', u, v, h, tv%T, tv%S, g)
1642 if (
associated(tv%T))
then
1645 call hchksum(ea_t,
"before triDiagTS ea_t ",g%HI,haloshift=0, scale=gv%H_to_m)
1646 call hchksum(eb_t,
"before triDiagTS eb_t ",g%HI,haloshift=0, scale=gv%H_to_m)
1647 call hchksum(ea_s,
"before triDiagTS ea_s ",g%HI,haloshift=0, scale=gv%H_to_m)
1648 call hchksum(eb_s,
"before triDiagTS eb_s ",g%HI,haloshift=0, scale=gv%H_to_m)
1651 call cpu_clock_begin(id_clock_tridiag)
1658 if (
associated(tv%S) .and.
associated(tv%salt_deficit)) &
1659 call adjust_salt(h, tv, g, gv, cs%diabatic_aux_CSp)
1661 if (cs%diabatic_diff_tendency_diag)
then
1662 do k=1,nz ;
do j=js,je ;
do i=is,ie
1663 temp_diag(i,j,k) = tv%T(i,j,k)
1664 saln_diag(i,j,k) = tv%S(i,j,k)
1665 enddo ;
enddo ;
enddo
1671 do j=js,je ;
do i=is,ie
1672 ea_t(i,j,1) = 0.; ea_s(i,j,1) = 0.
1676 do k=2,nz ;
do j=js,je ;
do i=is,ie
1677 hval = 1.0 / (h_neglect + 0.5*(h(i,j,k-1) + h(i,j,k)))
1678 ea_t(i,j,k) = (gv%Z_to_H**2) * dt_in_t * hval * kd_heat(i,j,k)
1679 eb_t(i,j,k-1) = ea_t(i,j,k)
1680 ea_s(i,j,k) = (gv%Z_to_H**2) * dt_in_t * hval * kd_salt(i,j,k)
1681 eb_s(i,j,k-1) = ea_s(i,j,k)
1682 enddo ;
enddo ;
enddo
1683 do j=js,je ;
do i=is,ie
1684 eb_t(i,j,nz) = 0. ; eb_s(i,j,nz) = 0.
1686 if (showcalltree)
call calltree_waypoint(
"done setting ea_t,ea_s,eb_t,eb_s from Kd_heat" //&
1687 "and Kd_salt (diabatic)")
1693 ea_t(i,js-1,k) = 0.0 ; eb_t(i,js-1,k) = 0.0
1694 ea_s(i,js-1,k) = 0.0 ; eb_s(i,js-1,k) = 0.0
1695 ea_t(i,je+1,k) = 0.0 ; eb_t(i,je+1,k) = 0.0
1696 ea_s(i,je+1,k) = 0.0 ; eb_s(i,je+1,k) = 0.0
1699 ea_t(is-1,j,k) = 0.0 ; eb_t(is-1,j,k) = 0.0
1700 ea_s(is-1,j,k) = 0.0 ; eb_s(is-1,j,k) = 0.0
1701 ea_t(ie+1,j,k) = 0.0 ; eb_t(ie+1,j,k) = 0.0
1702 ea_s(ie+1,j,k) = 0.0 ; eb_s(ie+1,j,k) = 0.0
1707 call tracer_vertdiff(h, ea_t, eb_t, dt, tv%T, g, gv)
1708 call tracer_vertdiff(h, ea_s, eb_s, dt, tv%S, g, gv)
1712 if (cs%diabatic_diff_tendency_diag)
then
1713 call diagnose_diabatic_diff_tendency(tv, h, temp_diag, saln_diag, dt, g, gv, cs)
1715 call cpu_clock_end(id_clock_tridiag)
1717 if (showcalltree)
call calltree_waypoint(
"done with triDiagTS (diabatic)")
1721 if (cs%debugConservation)
call mom_state_stats(
'triDiagTS', u, v, h, tv%T, tv%S, g)
1725 call mom_thermovar_chksum(
"after mixed layer ", tv, g)
1730 call diag_update_remap_grids(cs%diag)
1734 if ((cs%id_Tdif > 0) .or. (cs%id_Tadv > 0))
then
1735 do j=js,je ;
do i=is,ie
1736 tdif_flx(i,j,1) = 0.0 ; tdif_flx(i,j,nz+1) = 0.0
1737 tadv_flx(i,j,1) = 0.0 ; tadv_flx(i,j,nz+1) = 0.0
1740 do k=2,nz ;
do j=js,je ;
do i=is,ie
1741 tdif_flx(i,j,k) = (idt * 0.5*(ea_t(i,j,k) + eb_t(i,j,k-1))) * &
1742 (tv%T(i,j,k-1) - tv%T(i,j,k))
1743 tadv_flx(i,j,k) = (idt * (ea_t(i,j,k) - eb_t(i,j,k-1))) * &
1744 0.5*(tv%T(i,j,k-1) + tv%T(i,j,k))
1745 enddo ;
enddo ;
enddo
1747 if ((cs%id_Sdif > 0) .or. (cs%id_Sadv > 0))
then
1748 do j=js,je ;
do i=is,ie
1749 sdif_flx(i,j,1) = 0.0 ; sdif_flx(i,j,nz+1) = 0.0
1750 sadv_flx(i,j,1) = 0.0 ; sadv_flx(i,j,nz+1) = 0.0
1753 do k=2,nz ;
do j=js,je ;
do i=is,ie
1754 sdif_flx(i,j,k) = (idt * 0.5*(ea_s(i,j,k) + eb_s(i,j,k-1))) * &
1755 (tv%S(i,j,k-1) - tv%S(i,j,k))
1756 sadv_flx(i,j,k) = (idt * (ea_s(i,j,k) - eb_s(i,j,k-1))) * &
1757 0.5*(tv%S(i,j,k-1) + tv%S(i,j,k))
1758 enddo ;
enddo ;
enddo
1762 call cpu_clock_begin(id_clock_tracers)
1764 if (cs%mix_boundary_tracers)
then
1765 tr_ea_bbl = gv%Z_to_H * sqrt(dt_in_t*cs%Kd_BBL_tr)
1769 ebtr(i,j,nz) = eb_s(i,j,nz)
1771 in_boundary(i) = (g%mask2dT(i,j) > 0.0)
1773 do k=nz,2,-1 ;
do i=is,ie
1774 if (in_boundary(i))
then
1775 htot(i) = htot(i) + h(i,j,k)
1784 add_ent = ((dt_in_t * cs%Kd_min_tr) * gv%Z_to_H**2) * &
1785 ((h(i,j,k-1)+h(i,j,k)+h_neglect) / &
1786 (h(i,j,k-1)*h(i,j,k)+h_neglect2)) - &
1787 0.5*(ea_s(i,j,k) + eb_s(i,j,k-1))
1788 if (htot(i) < tr_ea_bbl)
then
1789 add_ent = max(0.0, add_ent, &
1790 (tr_ea_bbl - htot(i)) - min(ea_s(i,j,k),eb_s(i,j,k-1)))
1791 elseif (add_ent < 0.0)
then
1792 add_ent = 0.0 ; in_boundary(i) = .false.
1795 ebtr(i,j,k-1) = eb_s(i,j,k-1) + add_ent
1796 eatr(i,j,k) = ea_s(i,j,k) + add_ent
1798 ebtr(i,j,k-1) = eb_s(i,j,k-1) ; eatr(i,j,k) = ea_s(i,j,k)
1801 if (
associated(visc%Kd_extra_S))
then ;
if (visc%Kd_extra_S(i,j,k) > 0.0)
then
1802 add_ent = ((dt_in_t * visc%Kd_extra_S(i,j,k)) * gv%Z_to_H**2) / &
1803 (0.5 * (h(i,j,k-1) + h(i,j,k)) + &
1805 ebtr(i,j,k-1) = ebtr(i,j,k-1) + add_ent
1806 eatr(i,j,k) = eatr(i,j,k) + add_ent
1809 do i=is,ie ; eatr(i,j,1) = ea_s(i,j,1) ;
enddo
1815 call call_tracer_column_fns(h_prebound, h, ea_s, eb_s, fluxes, hml, dt, g, gv, tv, &
1816 cs%optics, cs%tracer_flow_CSp, cs%debug, &
1817 evap_cfl_limit = cs%evap_CFL_limit, &
1818 minimum_forcing_depth = cs%minimum_forcing_depth)
1820 elseif (
associated(visc%Kd_extra_S))
then
1822 do j=js,je ;
do i=is,ie
1823 ebtr(i,j,nz) = eb_s(i,j,nz) ; eatr(i,j,1) = ea_s(i,j,1)
1826 do k=nz,2,-1 ;
do j=js,je ;
do i=is,ie
1827 if (visc%Kd_extra_S(i,j,k) > 0.0)
then
1828 add_ent = ((dt_in_t * visc%Kd_extra_S(i,j,k)) * gv%Z_to_H**2) / &
1829 (0.5 * (h(i,j,k-1) + h(i,j,k)) + &
1834 ebtr(i,j,k-1) = eb_s(i,j,k-1) + add_ent
1835 eatr(i,j,k) = ea_s(i,j,k) + add_ent
1836 enddo ;
enddo ;
enddo
1839 call call_tracer_column_fns(h_prebound, h, eatr, ebtr, fluxes, hml, dt, g, gv, tv, &
1840 cs%optics, cs%tracer_flow_CSp, cs%debug,&
1841 evap_cfl_limit = cs%evap_CFL_limit, &
1842 minimum_forcing_depth = cs%minimum_forcing_depth)
1845 call call_tracer_column_fns(h_prebound, h, eatr, ebtr, fluxes, hml, dt, g, gv, tv, &
1846 cs%optics, cs%tracer_flow_CSp, cs%debug, &
1847 evap_cfl_limit = cs%evap_CFL_limit, &
1848 minimum_forcing_depth = cs%minimum_forcing_depth)
1851 call cpu_clock_end(id_clock_tracers)
1854 if (cs%use_sponge)
then
1855 call cpu_clock_begin(id_clock_sponge)
1856 if (
associated(cs%ALE_sponge_CSp))
then
1857 call apply_ale_sponge(h, dt, g, gv, us, cs%ALE_sponge_CSp, cs%Time)
1860 call cpu_clock_end(id_clock_sponge)
1863 call mom_thermovar_chksum(
"apply_sponge ", tv, g)
1867 call cpu_clock_begin(id_clock_pass)
1868 if (g%symmetric)
then ; dir_flag = to_all+omit_corners
1869 else ; dir_flag = to_west+to_south+omit_corners ;
endif
1874 call do_group_pass(cs%pass_hold_eb_ea, g%Domain)
1876 if (
associated(visc%Kv_slow)) &
1877 call pass_var(visc%Kv_slow, g%Domain, to_all+omit_corners, halo=1)
1878 call cpu_clock_end(id_clock_pass)
1880 call disable_averaging(cs%diag)
1883 call enable_averaging(dt, time_end, cs%diag)
1885 if (cs%id_Kd_interface > 0)
call post_data(cs%id_Kd_interface, kd_int, cs%diag)
1886 if (cs%id_Kd_heat > 0)
call post_data(cs%id_Kd_heat, kd_heat, cs%diag)
1887 if (cs%id_Kd_salt > 0)
call post_data(cs%id_Kd_salt, kd_salt, cs%diag)
1888 if (cs%id_Kd_ePBL > 0)
call post_data(cs%id_Kd_ePBL, kd_epbl, cs%diag)
1890 if (cs%id_ea_t > 0)
call post_data(cs%id_ea_t, ea_t, cs%diag)
1891 if (cs%id_eb_t > 0)
call post_data(cs%id_eb_t, eb_t, cs%diag)
1892 if (cs%id_ea_s > 0)
call post_data(cs%id_ea_s, ea_s, cs%diag)
1893 if (cs%id_eb_s > 0)
call post_data(cs%id_eb_s, eb_s, cs%diag)
1895 if (cs%id_dudt_dia > 0)
call post_data(cs%id_dudt_dia, adp%du_dt_dia, cs%diag)
1896 if (cs%id_dvdt_dia > 0)
call post_data(cs%id_dvdt_dia, adp%dv_dt_dia, cs%diag)
1898 if (cs%id_Tdif > 0)
call post_data(cs%id_Tdif, tdif_flx, cs%diag)
1899 if (cs%id_Tadv > 0)
call post_data(cs%id_Tadv, tadv_flx, cs%diag)
1900 if (cs%id_Sdif > 0)
call post_data(cs%id_Sdif, sdif_flx, cs%diag)
1901 if (cs%id_Sadv > 0)
call post_data(cs%id_Sadv, sadv_flx, cs%diag)
1903 call disable_averaging(cs%diag)
1905 if (showcalltree)
call calltree_leave(
"diabatic_ALE()")
1907 end subroutine diabatic_ale
1911 subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, &
1912 G, GV, US, CS, WAVES)
1916 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(inout) :: u
1917 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(inout) :: v
1918 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(inout) :: h
1921 real,
dimension(:,:),
pointer :: Hml
1922 type(
forcing),
intent(inout) :: fluxes
1929 real,
intent(in) :: dt
1930 type(time_type),
intent(in) :: Time_end
1934 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
1950 real,
dimension(SZI_(G),SZJ_(G)) :: &
1953 real,
dimension(SZI_(G),SZJ_(G),G%ke) :: h_diag
1954 real,
dimension(SZI_(G),SZJ_(G),G%ke) :: temp_diag
1955 real,
dimension(SZI_(G),SZJ_(G),G%ke) :: saln_diag
1956 real,
dimension(SZI_(G),SZJ_(G)) :: tendency_2d
1960 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
target :: &
1966 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1),
target :: &
1977 real,
pointer,
dimension(:,:,:) :: &
1983 integer :: kb(SZI_(G),SZJ_(G))
1986 real :: p_ref_cv(SZI_(G))
1990 logical :: in_boundary(SZI_(G))
2010 real :: htot(SZIB_(G))
2011 real :: b1(SZIB_(G)), d1(SZIB_(G))
2012 real :: c1(SZIB_(G),SZK_(G))
2020 logical :: showCallTree
2021 integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb, m, halo
2026 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
2027 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
2028 nkmb = gv%nk_rho_varies
2029 h_neglect = gv%H_subroundoff ; h_neglect2 = h_neglect*h_neglect
2030 kd_heat(:,:,:) = 0.0 ; kd_salt(:,:,:) = 0.0
2033 showcalltree = calltree_showquery()
2034 if (showcalltree)
call calltree_enter(
"layered_diabatic(), MOM_diabatic_driver.F90")
2037 eaml => eatr ; ebml => ebtr
2038 dt_in_t = dt * us%s_to_T
2041 call enable_averaging(dt, time_end, cs%diag)
2043 if ((cs%ML_mix_first > 0.0) .or. cs%use_geothermal)
then
2044 halo = cs%halo_TS_diff
2046 do k=1,nz ;
do j=js-halo,je+halo ;
do i=is-halo,ie+halo
2047 h_orig(i,j,k) = h(i,j,k) ; eaml(i,j,k) = 0.0 ; ebml(i,j,k) = 0.0
2048 enddo ;
enddo ;
enddo
2051 if (cs%use_geothermal)
then
2052 call cpu_clock_begin(id_clock_geothermal)
2053 call geothermal(h, tv, dt, eaml, ebml, g, gv, cs%geothermal_CSp, halo=cs%halo_TS_diff)
2054 call cpu_clock_end(id_clock_geothermal)
2055 if (showcalltree)
call calltree_waypoint(
"geothermal (diabatic)")
2056 if (cs%debugConservation)
call mom_state_stats(
'geothermal', u, v, h, tv%T, tv%S, g)
2061 call diag_update_remap_grids(cs%diag)
2066 if (
associated(cs%optics)) &
2067 call set_pen_shortwave(cs%optics, fluxes, g, gv, cs%diabatic_aux_CSp, cs%opacity_CSp, cs%tracer_flow_CSp)
2069 if (cs%bulkmixedlayer)
then
2070 if (cs%debug)
call mom_forcing_chksum(
"Before mixedlayer", fluxes, g, us, haloshift=0)
2072 if (cs%ML_mix_first > 0.0)
then
2080 call find_uv_at_h(u, v, h, u_h, v_h, g, gv)
2082 call cpu_clock_begin(id_clock_mixedlayer)
2083 if (cs%ML_mix_first < 1.0)
then
2085 call bulkmixedlayer(h, u_h, v_h, tv, fluxes, dt_in_t*cs%ML_mix_first, &
2086 eaml,ebml, g, gv, us, cs%bulkmixedlayer_CSp, cs%optics, &
2087 hml, cs%aggregate_FW_forcing, dt_in_t, last_call=.false.)
2088 if (cs%salt_reject_below_ML) &
2089 call insert_brine(h, tv, g, gv, fluxes, nkmb, cs%diabatic_aux_CSp, &
2090 dt*cs%ML_mix_first, cs%id_brine_lay)
2093 call bulkmixedlayer(h, u_h, v_h, tv, fluxes, dt_in_t, eaml, ebml, &
2094 g, gv, us, cs%bulkmixedlayer_CSp, cs%optics, &
2095 hml, cs%aggregate_FW_forcing, dt_in_t, last_call=.true.)
2104 if (
associated(tv%S) .and.
associated(tv%salt_deficit)) &
2105 call adjust_salt(h, tv, g, gv, cs%diabatic_aux_CSp)
2106 call cpu_clock_end(id_clock_mixedlayer)
2109 call mom_forcing_chksum(
"After mixedlayer", fluxes, g, us, haloshift=0)
2111 if (showcalltree)
call calltree_waypoint(
"done with 1st bulkmixedlayer (diabatic)")
2112 if (cs%debugConservation)
call mom_state_stats(
'1st bulkmixedlayer', u, v, h, tv%T, tv%S, g)
2118 if (cs%use_kappa_shear .or. cs%use_CVMix_shear)
then
2119 if ((cs%ML_mix_first > 0.0) .or. cs%use_geothermal)
then
2120 call find_uv_at_h(u, v, h_orig, u_h, v_h, g, gv, eaml, ebml)
2122 call hchksum(eaml,
"after find_uv_at_h eaml",g%HI, scale=gv%H_to_m)
2123 call hchksum(ebml,
"after find_uv_at_h ebml",g%HI, scale=gv%H_to_m)
2126 call find_uv_at_h(u, v, h, u_h, v_h, g, gv)
2128 if (showcalltree)
call calltree_waypoint(
"done with find_uv_at_h (diabatic)")
2131 call cpu_clock_begin(id_clock_set_diffusivity)
2134 if ((cs%halo_TS_diff > 0) .and. (cs%ML_mix_first > 0.0))
then
2135 if (
associated(tv%T))
call pass_var(tv%T, g%Domain, halo=cs%halo_TS_diff, complete=.false.)
2136 if (
associated(tv%T))
call pass_var(tv%S, g%Domain, halo=cs%halo_TS_diff, complete=.false.)
2137 call pass_var(h, g%domain, halo=cs%halo_TS_diff, complete=.true.)
2139 call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, cs%optics, visc, dt_in_t, g, gv, us, &
2140 cs%set_diff_CSp, kd_lay, kd_int)
2141 call cpu_clock_end(id_clock_set_diffusivity)
2142 if (showcalltree)
call calltree_waypoint(
"done with set_diffusivity (diabatic)")
2145 call mom_state_chksum(
"after set_diffusivity ", u, v, h, g, gv, haloshift=0)
2146 call mom_forcing_chksum(
"after set_diffusivity ", fluxes, g, us, haloshift=0)
2147 call mom_thermovar_chksum(
"after set_diffusivity ", tv, g)
2148 call hchksum(kd_lay,
"after set_diffusivity Kd_lay", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
2149 call hchksum(kd_int,
"after set_diffusivity Kd_Int", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
2154 call cpu_clock_begin(id_clock_kpp)
2160 call calculatebuoyancyflux2d(g, gv, us, fluxes, cs%optics, h, tv%T, tv%S, tv, &
2161 cs%KPP_buoy_flux, cs%KPP_temp_flux, cs%KPP_salt_flux)
2167 do k=1,nz+1 ;
do j=js,je ;
do i=is,ie
2168 kd_salt(i,j,k) = kd_int(i,j,k)
2169 kd_heat(i,j,k) = kd_int(i,j,k)
2170 enddo ;
enddo ;
enddo
2172 if (
associated(visc%Kd_extra_S))
then
2174 do k=1,nz+1 ;
do j=js,je ;
do i=is,ie
2175 kd_salt(i,j,k) = kd_salt(i,j,k) + visc%Kd_extra_S(i,j,k)
2176 enddo ;
enddo ;
enddo
2178 if (
associated(visc%Kd_extra_T))
then
2180 do k=1,nz+1 ;
do j=js,je ;
do i=is,ie
2181 kd_heat(i,j,k) = kd_heat(i,j,k) + visc%Kd_extra_T(i,j,k)
2182 enddo ;
enddo ;
enddo
2185 call kpp_compute_bld(cs%KPP_CSp, g, gv, us, h, tv%T, tv%S, u, v, tv%eqn_of_state, &
2186 fluxes%ustar, cs%KPP_buoy_flux, waves=waves)
2188 call kpp_calculate(cs%KPP_CSp, g, gv, us, h, fluxes%ustar, cs%KPP_buoy_flux, kd_heat, &
2189 kd_salt, visc%Kv_shear, cs%KPP_NLTheat, cs%KPP_NLTscalar, waves=waves)
2191 if (
associated(hml))
then
2192 call kpp_get_bld(cs%KPP_CSp, hml(:,:), g)
2193 call pass_var(hml, g%domain, halo=1)
2195 if (
associated(visc%MLD)) visc%MLD(:,:) = hml(:,:)
2198 if (.not. cs%KPPisPassive)
then
2200 do k=1,nz+1 ;
do j=js,je ;
do i=is,ie
2201 kd_int(i,j,k) = min( kd_salt(i,j,k), kd_heat(i,j,k) )
2202 enddo ;
enddo ;
enddo
2203 if (
associated(visc%Kd_extra_S))
then
2205 do k=1,nz+1 ;
do j=js,je ;
do i=is,ie
2206 visc%Kd_extra_S(i,j,k) = (kd_salt(i,j,k) - kd_int(i,j,k))
2207 enddo ;
enddo ;
enddo
2209 if (
associated(visc%Kd_extra_T))
then
2211 do k=1,nz+1 ;
do j=js,je ;
do i=is,ie
2212 visc%Kd_extra_T(i,j,k) = (kd_heat(i,j,k) - kd_int(i,j,k))
2213 enddo ;
enddo ;
enddo
2217 call cpu_clock_end(id_clock_kpp)
2218 if (showcalltree)
call calltree_waypoint(
"done with KPP_calculate (diabatic)")
2221 call mom_forcing_chksum(
"after KPP", fluxes, g, us, haloshift=0)
2222 call mom_thermovar_chksum(
"after KPP", tv, g)
2223 call hchksum(kd_lay,
"after KPP Kd_lay", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
2224 call hchksum(kd_int,
"after KPP Kd_Int", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
2230 if (cs%use_CVMix_conv)
then
2231 call calculate_cvmix_conv(h, tv, g, gv, us, cs%CVMix_conv_csp, hml)
2233 do k=1,nz+1 ;
do j=js,je ;
do i=is,ie
2234 kd_int(i,j,k) = kd_int(i,j,k) + cs%CVMix_conv_csp%kd_conv(i,j,k)
2235 visc%Kv_slow(i,j,k) = visc%Kv_slow(i,j,k) + cs%CVMix_conv_csp%kv_conv(i,j,k)
2236 enddo ;
enddo ;
enddo
2240 call cpu_clock_begin(id_clock_kpp)
2242 call hchksum(cs%KPP_temp_flux,
"before KPP_applyNLT netHeat", g%HI, haloshift=0, scale=gv%H_to_m)
2243 call hchksum(cs%KPP_salt_flux,
"before KPP_applyNLT netSalt", g%HI, haloshift=0, scale=gv%H_to_m)
2244 call hchksum(cs%KPP_NLTheat,
"before KPP_applyNLT NLTheat", g%HI, haloshift=0)
2245 call hchksum(cs%KPP_NLTscalar,
"before KPP_applyNLT NLTscalar", g%HI, haloshift=0)
2249 call kpp_nonlocaltransport_temp(cs%KPP_CSp, g, gv, h, cs%KPP_NLTheat, cs%KPP_temp_flux, dt, tv%T, tv%C_p)
2250 call kpp_nonlocaltransport_saln(cs%KPP_CSp, g, gv, h, cs%KPP_NLTscalar, cs%KPP_salt_flux, dt, tv%S)
2251 call cpu_clock_end(id_clock_kpp)
2252 if (showcalltree)
call calltree_waypoint(
"done with KPP_applyNonLocalTransport (diabatic)")
2253 if (cs%debugConservation)
call mom_state_stats(
'KPP_applyNonLocalTransport', u, v, h, tv%T, tv%S, g)
2257 call mom_forcing_chksum(
"after KPP_applyNLT ", fluxes, g, us, haloshift=0)
2258 call mom_thermovar_chksum(
"after KPP_applyNLT ", tv, g)
2264 if (
associated(visc%Kd_extra_T) .and.
associated(visc%Kd_extra_S) .and.
associated(tv%T))
then
2266 call cpu_clock_begin(id_clock_differential_diff)
2267 call differential_diffuse_t_s(h, tv, visc, dt_in_t, g, gv)
2268 call cpu_clock_end(id_clock_differential_diff)
2269 if (showcalltree)
call calltree_waypoint(
"done with differential_diffuse_T_S (diabatic)")
2270 if (cs%debugConservation)
call mom_state_stats(
'differential_diffuse_T_S', u, v, h, tv%T, tv%S, g)
2274 if (.not. cs%useKPP)
then
2276 do k=2,nz ;
do j=js,je ;
do i=is,ie
2277 kd_heat(i,j,k) = kd_heat(i,j,k) + visc%Kd_extra_T(i,j,k)
2278 kd_salt(i,j,k) = kd_salt(i,j,k) + visc%Kd_extra_S(i,j,k)
2279 enddo ;
enddo ;
enddo
2289 call cpu_clock_begin(id_clock_entrain)
2292 call entrainment_diffusive(h, tv, fluxes, dt_in_t, g, gv, us, cs%entrain_diffusive_CSp, &
2293 ea, eb, kb, kd_lay=kd_lay, kd_int=kd_int)
2294 call cpu_clock_end(id_clock_entrain)
2295 if (showcalltree)
call calltree_waypoint(
"done with Entrainment_diffusive (diabatic)")
2298 call mom_forcing_chksum(
"after calc_entrain ", fluxes, g, us, haloshift=0)
2299 call mom_thermovar_chksum(
"after calc_entrain ", tv, g)
2301 call hchksum(ea,
"after calc_entrain ea", g%HI, haloshift=0, scale=gv%H_to_m)
2302 call hchksum(eb,
"after calc_entrain eb", g%HI, haloshift=0, scale=gv%H_to_m)
2306 if (cs%boundary_forcing_tendency_diag)
then
2307 do k=1,nz ;
do j=js,je ;
do i=is,ie
2308 h_diag(i,j,k) = h(i,j,k)
2309 temp_diag(i,j,k) = tv%T(i,j,k)
2310 saln_diag(i,j,k) = tv%S(i,j,k)
2311 enddo ;
enddo ;
enddo
2325 hold(i,j,1) = h(i,j,1)
2326 h(i,j,1) = h(i,j,1) + (eb(i,j,1) - ea(i,j,2))
2327 hold(i,j,nz) = h(i,j,nz)
2328 h(i,j,nz) = h(i,j,nz) + (ea(i,j,nz) - eb(i,j,nz-1))
2329 if (h(i,j,1) <= 0.0)
then
2330 h(i,j,1) = gv%Angstrom_H
2332 if (h(i,j,nz) <= 0.0)
then
2333 h(i,j,nz) = gv%Angstrom_H
2336 do k=2,nz-1 ;
do i=is,ie
2337 hold(i,j,k) = h(i,j,k)
2338 h(i,j,k) = h(i,j,k) + ((ea(i,j,k) - eb(i,j,k-1)) + &
2339 (eb(i,j,k) - ea(i,j,k+1)))
2340 if (h(i,j,k) <= 0.0)
then
2341 h(i,j,k) = gv%Angstrom_H
2346 call diag_update_remap_grids(cs%diag)
2349 call mom_state_chksum(
"after negative check ", u, v, h, g, gv, haloshift=0)
2350 call mom_forcing_chksum(
"after negative check ", fluxes, g, us, haloshift=0)
2351 call mom_thermovar_chksum(
"after negative check ", tv, g)
2353 if (showcalltree)
call calltree_waypoint(
"done with h=ea-eb (diabatic)")
2354 if (cs%debugConservation)
call mom_state_stats(
'h=ea-eb', u, v, h, tv%T, tv%S, g)
2360 if (cs%bulkmixedlayer)
then
2362 if (
associated(tv%T))
then
2363 call cpu_clock_begin(id_clock_tridiag)
2368 if (cs%massless_match_targets)
then
2372 h_tr = hold(i,j,1) + h_neglect
2373 b1(i) = 1.0 / (h_tr + eb(i,j,1))
2374 d1(i) = h_tr * b1(i)
2375 tv%T(i,j,1) = b1(i) * (h_tr*tv%T(i,j,1))
2376 tv%S(i,j,1) = b1(i) * (h_tr*tv%S(i,j,1))
2378 do k=2,nkmb ;
do i=is,ie
2379 c1(i,k) = eb(i,j,k-1) * b1(i)
2380 h_tr = hold(i,j,k) + h_neglect
2381 b_denom_1 = h_tr + d1(i)*ea(i,j,k)
2382 b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
2383 if (k<nkmb) d1(i) = b_denom_1 * b1(i)
2384 tv%T(i,j,k) = b1(i) * (h_tr*tv%T(i,j,k) + ea(i,j,k)*tv%T(i,j,k-1))
2385 tv%S(i,j,k) = b1(i) * (h_tr*tv%S(i,j,k) + ea(i,j,k)*tv%S(i,j,k-1))
2388 do k=nkmb+1,nz ;
do i=is,ie
2389 if (k == kb(i,j))
then
2390 c1(i,k) = eb(i,j,k-1) * b1(i)
2391 d1(i) = (((eb(i,j,nkmb)-eb(i,j,k-1)) + hold(i,j,nkmb) + h_neglect) + &
2392 d1(i)*ea(i,j,nkmb)) * b1(i)
2393 h_tr = hold(i,j,k) + h_neglect
2394 b_denom_1 = h_tr + d1(i)*ea(i,j,k)
2395 b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
2396 d1(i) = b_denom_1 * b1(i)
2397 tv%T(i,j,k) = b1(i) * (h_tr*tv%T(i,j,k) + ea(i,j,k)*tv%T(i,j,nkmb))
2398 tv%S(i,j,k) = b1(i) * (h_tr*tv%S(i,j,k) + ea(i,j,k)*tv%S(i,j,nkmb))
2399 elseif (k > kb(i,j))
then
2400 c1(i,k) = eb(i,j,k-1) * b1(i)
2401 h_tr = hold(i,j,k) + h_neglect
2402 b_denom_1 = h_tr + d1(i)*ea(i,j,k)
2403 b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
2404 d1(i) = b_denom_1 * b1(i)
2405 tv%T(i,j,k) = b1(i) * (h_tr*tv%T(i,j,k) + ea(i,j,k)*tv%T(i,j,k-1))
2406 tv%S(i,j,k) = b1(i) * (h_tr*tv%S(i,j,k) + ea(i,j,k)*tv%S(i,j,k-1))
2407 elseif (eb(i,j,k) < eb(i,j,k-1))
then
2412 tv%T(i,j,nkmb) = tv%T(i,j,nkmb) + b1(i) * (eb(i,j,k-1) - eb(i,j,k)) * tv%T(i,j,k)
2413 tv%S(i,j,nkmb) = tv%S(i,j,nkmb) + b1(i) * (eb(i,j,k-1) - eb(i,j,k)) * tv%S(i,j,k)
2417 do k=nz-1,nkmb,-1 ;
do i=is,ie
2418 if (k >= kb(i,j))
then
2419 tv%T(i,j,k) = tv%T(i,j,k) + c1(i,k+1)*tv%T(i,j,k+1)
2420 tv%S(i,j,k) = tv%S(i,j,k) + c1(i,k+1)*tv%S(i,j,k+1)
2423 do i=is,ie ;
if (kb(i,j) <= nz)
then
2424 tv%T(i,j,nkmb) = tv%T(i,j,nkmb) + c1(i,kb(i,j))*tv%T(i,j,kb(i,j))
2425 tv%S(i,j,nkmb) = tv%S(i,j,nkmb) + c1(i,kb(i,j))*tv%S(i,j,kb(i,j))
2427 do k=nkmb-1,1,-1 ;
do i=is,ie
2428 tv%T(i,j,k) = tv%T(i,j,k) + c1(i,k+1)*tv%T(i,j,k+1)
2429 tv%S(i,j,k) = tv%S(i,j,k) + c1(i,k+1)*tv%S(i,j,k+1)
2436 if (cs%tracer_tridiag)
then
2437 call tracer_vertdiff(hold, ea, eb, dt, tv%T, g, gv)
2438 call tracer_vertdiff(hold, ea, eb, dt, tv%S, g, gv)
2440 call tridiagts(g, gv, is, ie, js, je, hold, ea, eb, tv%T, tv%S)
2443 call cpu_clock_end(id_clock_tridiag)
2446 if (cs%debugConservation)
call mom_state_stats(
'BML tridiag', u, v, h, tv%T, tv%S, g)
2448 if ((cs%ML_mix_first > 0.0) .or. cs%use_geothermal)
then
2452 do k=1,nz ;
do j=js,je ;
do i=is,ie
2453 hold(i,j,k) = h_orig(i,j,k)
2454 ea(i,j,k) = ea(i,j,k) + eaml(i,j,k)
2455 eb(i,j,k) = eb(i,j,k) + ebml(i,j,k)
2456 enddo ;
enddo ;
enddo
2458 call hchksum(ea,
"after ea = ea + eaml",g%HI,haloshift=0, scale=gv%H_to_m)
2459 call hchksum(eb,
"after eb = eb + ebml",g%HI,haloshift=0, scale=gv%H_to_m)
2463 if (cs%ML_mix_first < 1.0)
then
2472 call find_uv_at_h(u, v, hold, u_h, v_h, g, gv, ea, eb)
2473 if (cs%debug)
call mom_state_chksum(
"find_uv_at_h1 ", u, v, h, g, gv, haloshift=0)
2475 dt_mix = min(dt_in_t, dt_in_t*(1.0 - cs%ML_mix_first))
2476 call cpu_clock_begin(id_clock_mixedlayer)
2478 call bulkmixedlayer(h, u_h, v_h, tv, fluxes, dt_mix, ea, eb, &
2479 g, gv, us, cs%bulkmixedlayer_CSp, cs%optics, &
2480 hml, cs%aggregate_FW_forcing, dt_in_t, last_call=.true.)
2482 if (cs%salt_reject_below_ML) &
2483 call insert_brine(h, tv, g, gv, fluxes, nkmb, cs%diabatic_aux_CSp, us%T_to_s*dt_mix, &
2492 if (
associated(tv%S) .and.
associated(tv%salt_deficit)) &
2493 call adjust_salt(h, tv, g, gv, cs%diabatic_aux_CSp)
2495 call cpu_clock_end(id_clock_mixedlayer)
2496 if (showcalltree)
call calltree_waypoint(
"done with 2nd bulkmixedlayer (diabatic)")
2497 if (cs%debugConservation)
call mom_state_stats(
'2nd bulkmixedlayer', u, v, h, tv%T, tv%S, g)
2503 if (
associated(tv%T))
then
2506 call hchksum(ea,
"before triDiagTS ea ",g%HI,haloshift=0, scale=gv%H_to_m)
2507 call hchksum(eb,
"before triDiagTS eb ",g%HI,haloshift=0, scale=gv%H_to_m)
2509 call cpu_clock_begin(id_clock_tridiag)
2517 if (
associated(tv%S) .and.
associated(tv%salt_deficit)) &
2518 call adjust_salt(h, tv, g, gv, cs%diabatic_aux_CSp)
2520 if (cs%diabatic_diff_tendency_diag)
then
2521 do k=1,nz ;
do j=js,je ;
do i=is,ie
2522 temp_diag(i,j,k) = tv%T(i,j,k)
2523 saln_diag(i,j,k) = tv%S(i,j,k)
2524 enddo ;
enddo ;
enddo
2528 if (cs%tracer_tridiag)
then
2529 call tracer_vertdiff(hold, ea, eb, dt, tv%T, g, gv)
2530 call tracer_vertdiff(hold, ea, eb, dt, tv%S, g, gv)
2532 call tridiagts(g, gv, is, ie, js, je, hold, ea, eb, tv%T, tv%S)
2538 if (cs%diabatic_diff_tendency_diag)
then
2539 call diagnose_diabatic_diff_tendency(tv, hold, temp_diag, saln_diag, dt, g, gv, cs)
2540 if (cs%id_diabatic_diff_h > 0)
call post_data(cs%id_diabatic_diff_h, hold, cs%diag, alt_h = hold)
2543 call cpu_clock_end(id_clock_tridiag)
2544 if (showcalltree)
call calltree_waypoint(
"done with triDiagTS (diabatic)")
2547 if (cs%debugConservation)
call mom_state_stats(
'triDiagTS', u, v, h, tv%T, tv%S, g)
2553 call mom_thermovar_chksum(
"after mixed layer ", tv, g)
2554 call hchksum(ea,
"after mixed layer ea", g%HI, scale=gv%H_to_m)
2555 call hchksum(eb,
"after mixed layer eb", g%HI, scale=gv%H_to_m)
2558 call cpu_clock_begin(id_clock_remap)
2559 call regularize_layers(h, tv, dt, ea, eb, g, gv, cs%regularize_layers_CSp)
2560 call cpu_clock_end(id_clock_remap)
2561 if (showcalltree)
call calltree_waypoint(
"done with regularize_layers (diabatic)")
2562 if (cs%debugConservation)
call mom_state_stats(
'regularize_layers', u, v, h, tv%T, tv%S, g)
2566 if (cs%id_dudt_dia > 0 .or. cs%id_dvdt_dia > 0) &
2568 call pass_var(h, g%domain, to_west+to_south+omit_corners, halo=1)
2569 call diag_update_remap_grids(cs%diag)
2573 if ((cs%id_Tdif > 0) .or. (cs%id_Tadv > 0))
then
2574 do j=js,je ;
do i=is,ie
2575 tdif_flx(i,j,1) = 0.0 ; tdif_flx(i,j,nz+1) = 0.0
2576 tadv_flx(i,j,1) = 0.0 ; tadv_flx(i,j,nz+1) = 0.0
2579 do k=2,nz ;
do j=js,je ;
do i=is,ie
2580 tdif_flx(i,j,k) = (idt * 0.5*(ea(i,j,k) + eb(i,j,k-1))) * &
2581 (tv%T(i,j,k-1) - tv%T(i,j,k))
2582 tadv_flx(i,j,k) = (idt * (ea(i,j,k) - eb(i,j,k-1))) * &
2583 0.5*(tv%T(i,j,k-1) + tv%T(i,j,k))
2584 enddo ;
enddo ;
enddo
2586 if ((cs%id_Sdif > 0) .or. (cs%id_Sadv > 0))
then
2587 do j=js,je ;
do i=is,ie
2588 sdif_flx(i,j,1) = 0.0 ; sdif_flx(i,j,nz+1) = 0.0
2589 sadv_flx(i,j,1) = 0.0 ; sadv_flx(i,j,nz+1) = 0.0
2592 do k=2,nz ;
do j=js,je ;
do i=is,ie
2593 sdif_flx(i,j,k) = (idt * 0.5*(ea(i,j,k) + eb(i,j,k-1))) * &
2594 (tv%S(i,j,k-1) - tv%S(i,j,k))
2595 sadv_flx(i,j,k) = (idt * (ea(i,j,k) - eb(i,j,k-1))) * &
2596 0.5*(tv%S(i,j,k-1) + tv%S(i,j,k))
2597 enddo ;
enddo ;
enddo
2601 call cpu_clock_begin(id_clock_tracers)
2602 if (cs%mix_boundary_tracers)
then
2603 tr_ea_bbl = gv%Z_to_H * sqrt(dt_in_t*cs%Kd_BBL_tr)
2607 ebtr(i,j,nz) = eb(i,j,nz)
2609 in_boundary(i) = (g%mask2dT(i,j) > 0.0)
2611 do k=nz,2,-1 ;
do i=is,ie
2612 if (in_boundary(i))
then
2613 htot(i) = htot(i) + h(i,j,k)
2622 add_ent = ((dt_in_t * cs%Kd_min_tr) * gv%Z_to_H**2) * &
2623 ((h(i,j,k-1)+h(i,j,k)+h_neglect) / &
2624 (h(i,j,k-1)*h(i,j,k)+h_neglect2)) - &
2625 0.5*(ea(i,j,k) + eb(i,j,k-1))
2626 if (htot(i) < tr_ea_bbl)
then
2627 add_ent = max(0.0, add_ent, &
2628 (tr_ea_bbl - htot(i)) - min(ea(i,j,k),eb(i,j,k-1)))
2629 elseif (add_ent < 0.0)
then
2630 add_ent = 0.0 ; in_boundary(i) = .false.
2633 ebtr(i,j,k-1) = eb(i,j,k-1) + add_ent
2634 eatr(i,j,k) = ea(i,j,k) + add_ent
2636 ebtr(i,j,k-1) = eb(i,j,k-1) ; eatr(i,j,k) = ea(i,j,k)
2638 if (
associated(visc%Kd_extra_S))
then ;
if (visc%Kd_extra_S(i,j,k) > 0.0)
then
2639 add_ent = ((dt_in_t * visc%Kd_extra_S(i,j,k)) * gv%Z_to_H**2) / &
2640 (0.25 * ((h(i,j,k-1) + h(i,j,k)) + (hold(i,j,k-1) + hold(i,j,k))) + &
2642 ebtr(i,j,k-1) = ebtr(i,j,k-1) + add_ent
2643 eatr(i,j,k) = eatr(i,j,k) + add_ent
2646 do i=is,ie ; eatr(i,j,1) = ea(i,j,1) ;
enddo
2650 call call_tracer_column_fns(hold, h, eatr, ebtr, fluxes, hml, dt, g, gv, tv, &
2651 cs%optics, cs%tracer_flow_CSp, cs%debug)
2653 elseif (
associated(visc%Kd_extra_S))
then
2655 do j=js,je ;
do i=is,ie
2656 ebtr(i,j,nz) = eb(i,j,nz) ; eatr(i,j,1) = ea(i,j,1)
2659 do k=nz,2,-1 ;
do j=js,je ;
do i=is,ie
2660 if (visc%Kd_extra_S(i,j,k) > 0.0)
then
2661 add_ent = ((dt_in_t * visc%Kd_extra_S(i,j,k)) * gv%Z_to_H**2) / &
2662 (0.25 * ((h(i,j,k-1) + h(i,j,k)) + (hold(i,j,k-1) + hold(i,j,k))) + &
2667 ebtr(i,j,k-1) = eb(i,j,k-1) + add_ent
2668 eatr(i,j,k) = ea(i,j,k) + add_ent
2669 enddo ;
enddo ;
enddo
2671 call call_tracer_column_fns(hold, h, eatr, ebtr, fluxes, hml, dt, g, gv, tv, &
2672 cs%optics, cs%tracer_flow_CSp, cs%debug)
2675 call call_tracer_column_fns(hold, h, ea, eb, fluxes, hml, dt, g, gv, tv, &
2676 cs%optics, cs%tracer_flow_CSp, cs%debug)
2680 call cpu_clock_end(id_clock_tracers)
2683 if (cs%use_sponge)
then
2684 call cpu_clock_begin(id_clock_sponge)
2686 if (cs%bulkmixedlayer .and.
associated(tv%eqn_of_state))
then
2687 do i=is,ie ; p_ref_cv(i) = tv%P_Ref ;
enddo
2691 is, ie-is+1, tv%eqn_of_state)
2693 call apply_sponge(h, dt, g, gv, ea, eb, cs%sponge_CSp, rcv_ml)
2695 call apply_sponge(h, dt, g, gv, ea, eb, cs%sponge_CSp)
2697 call cpu_clock_end(id_clock_sponge)
2700 call mom_thermovar_chksum(
"apply_sponge ", tv, g)
2705 if (
associated(cdp%diapyc_vel))
then
2708 do k=2,nz ;
do i=is,ie
2709 cdp%diapyc_vel(i,j,k) = idt * (gv%H_to_m * (ea(i,j,k) - eb(i,j,k-1)))
2712 cdp%diapyc_vel(i,j,1) = 0.0
2713 cdp%diapyc_vel(i,j,nz+1) = 0.0
2721 if (cs%bulkmixedlayer)
then
2723 call hchksum(ea,
"before net flux rearrangement ea",g%HI, scale=gv%H_to_m)
2724 call hchksum(eb,
"before net flux rearrangement eb",g%HI, scale=gv%H_to_m)
2728 do k=2,gv%nkml ;
do i=is,ie
2729 net_ent = ea(i,j,k) - eb(i,j,k-1)
2730 ea(i,j,k) = max(net_ent, 0.0)
2731 eb(i,j,k-1) = max(-net_ent, 0.0)
2735 call hchksum(ea,
"after net flux rearrangement ea",g%HI, scale=gv%H_to_m)
2736 call hchksum(eb,
"after net flux rearrangement eb",g%HI, scale=gv%H_to_m)
2744 hold(i,js-1,k) = gv%Angstrom_H ; ea(i,js-1,k) = 0.0 ; eb(i,js-1,k) = 0.0
2745 hold(i,je+1,k) = gv%Angstrom_H ; ea(i,je+1,k) = 0.0 ; eb(i,je+1,k) = 0.0
2748 hold(is-1,j,k) = gv%Angstrom_H ; ea(is-1,j,k) = 0.0 ; eb(is-1,j,k) = 0.0
2749 hold(ie+1,j,k) = gv%Angstrom_H ; ea(ie+1,j,k) = 0.0 ; eb(ie+1,j,k) = 0.0
2753 call cpu_clock_begin(id_clock_pass)
2754 if (g%symmetric)
then ; dir_flag = to_all+omit_corners
2755 else ; dir_flag = to_west+to_south+omit_corners ;
endif
2759 call do_group_pass(cs%pass_hold_eb_ea, g%Domain)
2760 call cpu_clock_end(id_clock_pass)
2767 call hchksum(ea,
"before u/v tridiag ea",g%HI, scale=gv%H_to_m)
2768 call hchksum(eb,
"before u/v tridiag eb",g%HI, scale=gv%H_to_m)
2769 call hchksum(hold,
"before u/v tridiag hold",g%HI, scale=gv%H_to_m)
2771 call cpu_clock_begin(id_clock_tridiag)
2775 if (
associated(adp%du_dt_dia)) adp%du_dt_dia(i,j,1) = u(i,j,1)
2776 hval = (hold(i,j,1) + hold(i+1,j,1)) + (ea(i,j,1) + ea(i+1,j,1)) + h_neglect
2777 b1(i) = 1.0 / (hval + (eb(i,j,1) + eb(i+1,j,1)))
2778 d1(i) = hval * b1(i)
2779 u(i,j,1) = b1(i) * (hval * u(i,j,1))
2781 do k=2,nz ;
do i=isq,ieq
2782 if (
associated(adp%du_dt_dia)) adp%du_dt_dia(i,j,k) = u(i,j,k)
2783 c1(i,k) = (eb(i,j,k-1)+eb(i+1,j,k-1)) * b1(i)
2784 eaval = ea(i,j,k) + ea(i+1,j,k)
2785 hval = hold(i,j,k) + hold(i+1,j,k) + h_neglect
2786 b1(i) = 1.0 / ((eb(i,j,k) + eb(i+1,j,k)) + (hval + d1(i)*eaval))
2787 d1(i) = (hval + d1(i)*eaval) * b1(i)
2788 u(i,j,k) = (hval*u(i,j,k) + eaval*u(i,j,k-1))*b1(i)
2790 do k=nz-1,1,-1 ;
do i=isq,ieq
2791 u(i,j,k) = u(i,j,k) + c1(i,k+1)*u(i,j,k+1)
2792 if (
associated(adp%du_dt_dia)) &
2793 adp%du_dt_dia(i,j,k) = (u(i,j,k) - adp%du_dt_dia(i,j,k)) * idt
2795 if (
associated(adp%du_dt_dia))
then
2797 adp%du_dt_dia(i,j,nz) = (u(i,j,nz)-adp%du_dt_dia(i,j,nz)) * idt
2802 call mom_state_chksum(
"aft 1st loop tridiag ", u, v, h, g, gv, haloshift=0)
2807 if (
associated(adp%dv_dt_dia)) adp%dv_dt_dia(i,j,1) = v(i,j,1)
2808 hval = (hold(i,j,1) + hold(i,j+1,1)) + (ea(i,j,1) + ea(i,j+1,1)) + h_neglect
2809 b1(i) = 1.0 / (hval + (eb(i,j,1) + eb(i,j+1,1)))
2810 d1(i) = hval * b1(i)
2811 v(i,j,1) = b1(i) * (hval * v(i,j,1))
2813 do k=2,nz ;
do i=is,ie
2814 if (
associated(adp%dv_dt_dia)) adp%dv_dt_dia(i,j,k) = v(i,j,k)
2815 c1(i,k) = (eb(i,j,k-1)+eb(i,j+1,k-1)) * b1(i)
2816 eaval = ea(i,j,k) + ea(i,j+1,k)
2817 hval = hold(i,j,k) + hold(i,j+1,k) + h_neglect
2818 b1(i) = 1.0 / ((eb(i,j,k) + eb(i,j+1,k)) + (hval + d1(i)*eaval))
2819 d1(i) = (hval + d1(i)*eaval) * b1(i)
2820 v(i,j,k) = (hval*v(i,j,k) + eaval*v(i,j,k-1))*b1(i)
2822 do k=nz-1,1,-1 ;
do i=is,ie
2823 v(i,j,k) = v(i,j,k) + c1(i,k+1)*v(i,j,k+1)
2824 if (
associated(adp%dv_dt_dia)) &
2825 adp%dv_dt_dia(i,j,k) = (v(i,j,k) - adp%dv_dt_dia(i,j,k)) * idt
2827 if (
associated(adp%dv_dt_dia))
then
2829 adp%dv_dt_dia(i,j,nz) = (v(i,j,nz)-adp%dv_dt_dia(i,j,nz)) * idt
2833 call cpu_clock_end(id_clock_tridiag)
2838 call disable_averaging(cs%diag)
2840 call enable_averaging(dt, time_end, cs%diag)
2842 if (cs%id_Kd_interface > 0)
call post_data(cs%id_Kd_interface, kd_int, cs%diag)
2843 if (cs%id_Kd_heat > 0)
call post_data(cs%id_Kd_heat, kd_heat, cs%diag)
2844 if (cs%id_Kd_salt > 0)
call post_data(cs%id_Kd_salt, kd_salt, cs%diag)
2845 if (cs%id_Kd_ePBL > 0)
call post_data(cs%id_Kd_ePBL, kd_epbl, cs%diag)
2847 if (cs%id_ea > 0)
call post_data(cs%id_ea, ea, cs%diag)
2848 if (cs%id_eb > 0)
call post_data(cs%id_eb, eb, cs%diag)
2850 if (cs%id_dudt_dia > 0)
call post_data(cs%id_dudt_dia, adp%du_dt_dia, cs%diag)
2851 if (cs%id_dvdt_dia > 0)
call post_data(cs%id_dvdt_dia, adp%dv_dt_dia, cs%diag)
2852 if (cs%id_wd > 0)
call post_data(cs%id_wd, cdp%diapyc_vel, cs%diag)
2854 if (cs%id_Tdif > 0)
call post_data(cs%id_Tdif, tdif_flx, cs%diag)
2855 if (cs%id_Tadv > 0)
call post_data(cs%id_Tadv, tadv_flx, cs%diag)
2856 if (cs%id_Sdif > 0)
call post_data(cs%id_Sdif, sdif_flx, cs%diag)
2857 if (cs%id_Sadv > 0)
call post_data(cs%id_Sadv, sadv_flx, cs%diag)
2859 call disable_averaging(cs%diag)
2861 if (showcalltree)
call calltree_leave(
"layered_diabatic()")
2863 end subroutine layered_diabatic
2867 subroutine extract_diabatic_member(CS, opacity_CSp, optics_CSp, &
2868 evap_CFL_limit, minimum_forcing_depth, diabatic_aux_CSp)
2871 type(
opacity_cs),
optional,
pointer :: opacity_csp
2872 type(
optics_type),
optional,
pointer :: optics_csp
2873 real,
optional,
intent( out) :: evap_cfl_limit
2875 real,
optional,
intent( out) :: minimum_forcing_depth
2881 if (
present(opacity_csp)) opacity_csp => cs%opacity_CSp
2882 if (
present(optics_csp)) optics_csp => cs%optics
2883 if (
present(diabatic_aux_csp)) diabatic_aux_csp => cs%diabatic_aux_CSp
2886 if (
present(evap_cfl_limit)) evap_cfl_limit = cs%evap_CFL_limit
2887 if (
present(minimum_forcing_depth)) minimum_forcing_depth = cs%minimum_forcing_depth
2889 end subroutine extract_diabatic_member
2892 subroutine adiabatic(h, tv, fluxes, dt, G, GV, CS)
2894 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
2897 type(
forcing),
intent(inout) :: fluxes
2898 real,
intent(in) :: dt
2902 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: zeros
2906 call call_tracer_column_fns(h, h, zeros, zeros, fluxes, zeros(:,:,1), dt, g, gv, tv, &
2907 cs%optics, cs%tracer_flow_CSp, cs%debug)
2909 end subroutine adiabatic
2915 subroutine diagnose_diabatic_diff_tendency(tv, h, temp_old, saln_old, dt, G, GV, CS)
2919 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
2920 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: temp_old
2921 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: saln_old
2922 real,
intent(in) :: dt
2926 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: work_3d
2927 real,
dimension(SZI_(G),SZJ_(G)) :: work_2d
2929 real :: ppt2mks = 0.001
2930 integer :: i, j, k, is, ie, js, je, nz
2932 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
2933 idt = 0.0 ;
if (dt > 0.0) idt = 1. / dt
2934 work_3d(:,:,:) = 0.0
2939 do k=1,nz ;
do j=js,je ;
do i=is,ie
2940 work_3d(i,j,k) = (tv%T(i,j,k)-temp_old(i,j,k))*idt
2941 enddo ;
enddo ;
enddo
2942 if (cs%id_diabatic_diff_temp_tend > 0)
then
2943 call post_data(cs%id_diabatic_diff_temp_tend, work_3d, cs%diag, alt_h = h)
2947 if (cs%id_diabatic_diff_heat_tend > 0 .or. cs%id_diabatic_diff_heat_tend_2d > 0)
then
2948 do k=1,nz ;
do j=js,je ;
do i=is,ie
2949 work_3d(i,j,k) = h(i,j,k) * gv%H_to_kg_m2 * tv%C_p * work_3d(i,j,k)
2950 enddo ;
enddo ;
enddo
2951 if (cs%id_diabatic_diff_heat_tend > 0)
then
2952 call post_data(cs%id_diabatic_diff_heat_tend, work_3d, cs%diag, alt_h = h)
2954 if (cs%id_diabatic_diff_heat_tend_2d > 0)
then
2955 do j=js,je ;
do i=is,ie
2958 work_2d(i,j) = work_2d(i,j) + work_3d(i,j,k)
2961 call post_data(cs%id_diabatic_diff_heat_tend_2d, work_2d, cs%diag)
2966 if (cs%id_diabatic_diff_saln_tend > 0)
then
2967 do k=1,nz ;
do j=js,je ;
do i=is,ie
2968 work_3d(i,j,k) = (tv%S(i,j,k)-saln_old(i,j,k))*idt
2969 enddo ;
enddo ;
enddo
2970 call post_data(cs%id_diabatic_diff_saln_tend, work_3d, cs%diag, alt_h = h)
2974 if (cs%id_diabatic_diff_salt_tend > 0 .or. cs%id_diabatic_diff_salt_tend_2d > 0)
then
2975 do k=1,nz ;
do j=js,je ;
do i=is,ie
2976 work_3d(i,j,k) = h(i,j,k) * gv%H_to_kg_m2 * ppt2mks * work_3d(i,j,k)
2977 enddo ;
enddo ;
enddo
2978 if (cs%id_diabatic_diff_salt_tend > 0)
then
2979 call post_data(cs%id_diabatic_diff_salt_tend, work_3d, cs%diag, alt_h = h)
2981 if (cs%id_diabatic_diff_salt_tend_2d > 0)
then
2982 do j=js,je ;
do i=is,ie
2985 work_2d(i,j) = work_2d(i,j) + work_3d(i,j,k)
2988 call post_data(cs%id_diabatic_diff_salt_tend_2d, work_2d, cs%diag)
2992 end subroutine diagnose_diabatic_diff_tendency
2999 subroutine diagnose_boundary_forcing_tendency(tv, h, temp_old, saln_old, h_old, &
3004 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
3006 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
3007 intent(in) :: temp_old
3008 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
3009 intent(in) :: saln_old
3010 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
3012 real,
intent(in) :: dt
3016 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: work_3d
3017 real,
dimension(SZI_(G),SZJ_(G)) :: work_2d
3019 real :: ppt2mks = 0.001
3020 integer :: i, j, k, is, ie, js, je, nz
3022 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
3023 idt = 0.0 ;
if (dt > 0.0) idt = 1. / dt
3024 work_3d(:,:,:) = 0.0
3028 if (cs%id_boundary_forcing_h_tendency > 0)
then
3029 do k=1,nz ;
do j=js,je ;
do i=is,ie
3030 work_3d(i,j,k) = (h(i,j,k) - h_old(i,j,k))*idt
3031 enddo ;
enddo ;
enddo
3032 call post_data(cs%id_boundary_forcing_h_tendency, work_3d, cs%diag, alt_h = h_old)
3036 if (cs%id_boundary_forcing_temp_tend > 0)
then
3037 do k=1,nz ;
do j=js,je ;
do i=is,ie
3038 work_3d(i,j,k) = (tv%T(i,j,k)-temp_old(i,j,k))*idt
3039 enddo ;
enddo ;
enddo
3040 call post_data(cs%id_boundary_forcing_temp_tend, work_3d, cs%diag, alt_h = h_old)
3044 if (cs%id_boundary_forcing_heat_tend > 0 .or. cs%id_boundary_forcing_heat_tend_2d > 0)
then
3045 do k=1,nz ;
do j=js,je ;
do i=is,ie
3046 work_3d(i,j,k) = gv%H_to_kg_m2 * tv%C_p * idt * (h(i,j,k) * tv%T(i,j,k) - h_old(i,j,k) * temp_old(i,j,k))
3047 enddo ;
enddo ;
enddo
3048 if (cs%id_boundary_forcing_heat_tend > 0)
then
3049 call post_data(cs%id_boundary_forcing_heat_tend, work_3d, cs%diag, alt_h = h_old)
3051 if (cs%id_boundary_forcing_heat_tend_2d > 0)
then
3052 do j=js,je ;
do i=is,ie
3055 work_2d(i,j) = work_2d(i,j) + work_3d(i,j,k)
3058 call post_data(cs%id_boundary_forcing_heat_tend_2d, work_2d, cs%diag)
3063 if (cs%id_boundary_forcing_saln_tend > 0)
then
3064 do k=1,nz ;
do j=js,je ;
do i=is,ie
3065 work_3d(i,j,k) = (tv%S(i,j,k)-saln_old(i,j,k))*idt
3066 enddo ;
enddo ;
enddo
3067 call post_data(cs%id_boundary_forcing_saln_tend, work_3d, cs%diag, alt_h = h_old)
3071 if (cs%id_boundary_forcing_salt_tend > 0 .or. cs%id_boundary_forcing_salt_tend_2d > 0)
then
3072 do k=1,nz ;
do j=js,je ;
do i=is,ie
3073 work_3d(i,j,k) = gv%H_to_kg_m2 * ppt2mks * idt * (h(i,j,k) * tv%S(i,j,k) - h_old(i,j,k) * saln_old(i,j,k))
3074 enddo ;
enddo ;
enddo
3075 if (cs%id_boundary_forcing_salt_tend > 0)
then
3076 call post_data(cs%id_boundary_forcing_salt_tend, work_3d, cs%diag, alt_h = h_old)
3078 if (cs%id_boundary_forcing_salt_tend_2d > 0)
then
3079 do j=js,je ;
do i=is,ie
3082 work_2d(i,j) = work_2d(i,j) + work_3d(i,j,k)
3085 call post_data(cs%id_boundary_forcing_salt_tend_2d, work_2d, cs%diag)
3089 end subroutine diagnose_boundary_forcing_tendency
3096 subroutine diagnose_frazil_tendency(tv, h, temp_old, dt, G, GV, CS)
3101 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
3102 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: temp_old
3103 real,
intent(in) :: dt
3105 real,
dimension(SZI_(G),SZJ_(G)) :: work_2d
3107 integer :: i, j, k, is, ie, js, je, nz
3109 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
3110 idt = 0.0 ;
if (dt > 0.0) idt = 1. / dt
3113 if (cs%id_frazil_temp_tend > 0)
then
3114 do k=1,nz ;
do j=js,je ;
do i=is,ie
3115 cs%frazil_temp_diag(i,j,k) = idt * (tv%T(i,j,k)-temp_old(i,j,k))
3116 enddo ;
enddo ;
enddo
3117 call post_data(cs%id_frazil_temp_tend, cs%frazil_temp_diag(:,:,:), cs%diag)
3121 if (cs%id_frazil_heat_tend > 0 .or. cs%id_frazil_heat_tend_2d > 0)
then
3122 do k=1,nz ;
do j=js,je ;
do i=is,ie
3123 cs%frazil_heat_diag(i,j,k) = gv%H_to_kg_m2 * tv%C_p * h(i,j,k) * idt * (tv%T(i,j,k)-temp_old(i,j,k))
3124 enddo ;
enddo ;
enddo
3125 if (cs%id_frazil_heat_tend > 0)
call post_data(cs%id_frazil_heat_tend, cs%frazil_heat_diag(:,:,:), cs%diag)
3129 if (cs%id_frazil_heat_tend_2d > 0)
then
3130 do j=js,je ;
do i=is,ie
3133 work_2d(i,j) = work_2d(i,j) + cs%frazil_heat_diag(i,j,k)
3136 call post_data(cs%id_frazil_heat_tend_2d, work_2d, cs%diag)
3140 end subroutine diagnose_frazil_tendency
3146 subroutine adiabatic_driver_init(Time, G, param_file, diag, CS, &
3148 type(time_type),
intent(in) :: time
3151 type(
diag_ctrl),
target,
intent(inout) :: diag
3157 #include "version_variable.h"
3158 character(len=40) :: mdl =
"MOM_diabatic_driver"
3160 if (
associated(cs))
then
3161 call mom_error(warning,
"adiabatic_driver_init called with an "// &
3162 "associated control structure.")
3164 else ;
allocate(cs) ;
endif
3167 if (
associated(tracer_flow_csp)) cs%tracer_flow_CSp => tracer_flow_csp
3171 "The following parameters are used for diabatic processes.")
3173 end subroutine adiabatic_driver_init
3177 subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, diag, &
3178 ADp, CDp, CS, tracer_flow_CSp, sponge_CSp, &
3180 type(time_type),
target :: time
3185 logical,
intent(in) :: usealealgorithm
3186 type(
diag_ctrl),
target,
intent(inout) :: diag
3198 logical :: use_temperature, differentialdiffusion
3201 #include "version_variable.h"
3202 character(len=40) :: mdl =
"MOM_diabatic_driver"
3203 character(len=48) :: thickness_units
3204 character(len=40) :: var_name
3205 character(len=160) :: var_descript
3206 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb, nz, nbands, m
3207 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = g%ke
3208 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
3210 if (
associated(cs))
then
3211 call mom_error(warning,
"diabatic_driver_init called with an "// &
3212 "associated control structure.")
3221 if (
associated(tracer_flow_csp)) cs%tracer_flow_CSp => tracer_flow_csp
3222 if (
associated(sponge_csp)) cs%sponge_CSp => sponge_csp
3223 if (
associated(ale_sponge_csp)) cs%ALE_sponge_CSp => ale_sponge_csp
3225 cs%useALEalgorithm = usealealgorithm
3226 cs%bulkmixedlayer = (gv%nkml > 0)
3230 "The following parameters are used for diabatic processes.")
3231 call get_param(param_file, mdl,
"USE_LEGACY_DIABATIC_DRIVER", cs%use_legacy_diabatic, &
3232 "If true, use a legacy version of the diabatic subroutine. "//&
3233 "This is temporary and is needed to avoid change in answers.", &
3235 call get_param(param_file, mdl,
"SPONGE", cs%use_sponge, &
3236 "If true, sponges may be applied anywhere in the domain. "//&
3237 "The exact location and properties of those sponges are "//&
3238 "specified via calls to initialize_sponge and possibly "//&
3239 "set_up_sponge_field.", default=.false.)
3240 call get_param(param_file, mdl,
"ENABLE_THERMODYNAMICS", use_temperature, &
3241 "If true, temperature and salinity are used as state "//&
3242 "variables.", default=.true.)
3243 call get_param(param_file, mdl,
"ENERGETICS_SFC_PBL", cs%use_energetic_PBL, &
3244 "If true, use an implied energetics planetary boundary "//&
3245 "layer scheme to determine the diffusivity and viscosity "//&
3246 "in the surface boundary layer.", default=.false.)
3247 call get_param(param_file, mdl,
"EPBL_IS_ADDITIVE", cs%ePBL_is_additive, &
3248 "If true, the diffusivity from ePBL is added to all "//&
3249 "other diffusivities. Otherwise, the larger of kappa-shear "//&
3250 "and ePBL diffusivities are used.", default=.true.)
3251 call get_param(param_file, mdl,
"DOUBLE_DIFFUSION", differentialdiffusion, &
3252 "If true, apply parameterization of double-diffusion.", &
3254 call get_param(param_file, mdl,
"USE_KPP", cs%use_KPP, &
3255 "If true, turns on the [CVMix] KPP scheme of Large et al., 1994, "//&
3256 "to calculate diffusivities and non-local transport in the OBL.", &
3257 default=.false., do_not_log=.true.)
3258 cs%use_CVMix_ddiff = cvmix_ddiff_is_used(param_file)
3260 if (cs%use_CVMix_ddiff .and. differentialdiffusion)
then
3261 call mom_error(fatal,
'diabatic_driver_init: '// &
3262 'Multiple double-diffusion options selected (DOUBLE_DIFFUSION and'//&
3263 'USE_CVMIX_DDIFF), please disable all but one option to proceed.')
3266 cs%use_kappa_shear = kappa_shear_is_used(param_file)
3267 cs%use_CVMix_shear = cvmix_shear_is_used(param_file)
3269 if (cs%bulkmixedlayer)
then
3270 call get_param(param_file, mdl,
"ML_MIX_FIRST", cs%ML_mix_first, &
3271 "The fraction of the mixed layer mixing that is applied "//&
3272 "before interior diapycnal mixing. 0 by default.", &
3273 units=
"nondim", default=0.0)
3274 call get_param(param_file, mdl,
"NKBL", cs%nkbl, default=2, do_not_log=.true.)
3276 cs%ML_mix_first = 0.0
3278 if (use_temperature)
then
3279 call get_param(param_file, mdl,
"DO_GEOTHERMAL", cs%use_geothermal, &
3280 "If true, apply geothermal heating.", default=.false.)
3282 cs%use_geothermal = .false.
3284 call get_param(param_file, mdl,
"INTERNAL_TIDES", cs%use_int_tides, &
3285 "If true, use the code that advances a separate set of "//&
3286 "equations for the internal tide energy density.", default=.false.)
3288 if (cs%use_int_tides)
then
3289 call get_param(param_file, mdl,
"INTERNAL_TIDE_MODES", cs%nMode, &
3290 "The number of distinct internal tide modes "//&
3291 "that will be calculated.", default=1, do_not_log=.true.)
3292 call get_param(param_file, mdl,
"UNIFORM_TEST_CG", cs%uniform_test_cg, &
3293 "If positive, a uniform group velocity of internal tide for test case", &
3294 default=-1., units=
"m s-1")
3297 call get_param(param_file, mdl,
"MASSLESS_MATCH_TARGETS", &
3298 cs%massless_match_targets, &
3299 "If true, the temperature and salinity of massless layers "//&
3300 "are kept consistent with their target densities. "//&
3301 "Otherwise the properties of massless layers evolve "//&
3302 "diffusively to match massive neighboring layers.", &
3305 call get_param(param_file, mdl,
"AGGREGATE_FW_FORCING", cs%aggregate_FW_forcing, &
3306 "If true, the net incoming and outgoing fresh water fluxes are combined "//&
3307 "and applied as either incoming or outgoing depending on the sign of the net. "//&
3308 "If false, the net incoming fresh water flux is added to the model and "//&
3309 "thereafter the net outgoing is removed from the topmost non-vanished "//&
3310 "layers of the updated state.", default=.true.)
3312 call get_param(param_file, mdl,
"DEBUG", cs%debug, &
3313 "If true, write out verbose debugging data.", &
3314 default=.false., debuggingparam=.true.)
3315 call get_param(param_file, mdl,
"DEBUG_CONSERVATION", cs%debugConservation, &
3316 "If true, monitor conservation and extrema.", &
3317 default=.false., debuggingparam=.true.)
3319 call get_param(param_file, mdl,
"DEBUG_ENERGY_REQ", cs%debug_energy_req, &
3320 "If true, debug the energy requirements.", default=.false., do_not_log=.true.)
3321 call get_param(param_file, mdl,
"MIX_BOUNDARY_TRACERS", cs%mix_boundary_tracers, &
3322 "If true, mix the passive tracers in massless layers at "//&
3323 "the bottom into the interior as though a diffusivity of "//&
3324 "KD_MIN_TR were operating.", default=.true.)
3326 if (cs%mix_boundary_tracers)
then
3327 call get_param(param_file, mdl,
"KD", kd, fail_if_missing=.true.)
3328 call get_param(param_file, mdl,
"KD_MIN_TR", cs%Kd_min_tr, &
3329 "A minimal diffusivity that should always be applied to "//&
3330 "tracers, especially in massless layers near the bottom. "//&
3331 "The default is 0.1*KD.", units=
"m2 s-1", default=0.1*kd, scale=us%m2_s_to_Z2_T)
3332 call get_param(param_file, mdl,
"KD_BBL_TR", cs%Kd_BBL_tr, &
3333 "A bottom boundary layer tracer diffusivity that will "//&
3334 "allow for explicitly specified bottom fluxes. The "//&
3335 "entrainment at the bottom is at least sqrt(Kd_BBL_tr*dt) "//&
3336 "over the same distance.", units=
"m2 s-1", default=0., scale=us%m2_s_to_Z2_T)
3339 call get_param(param_file, mdl,
"TRACER_TRIDIAG", cs%tracer_tridiag, &
3340 "If true, use the passive tracer tridiagonal solver for T and S", &
3343 call get_param(param_file, mdl,
"MINIMUM_FORCING_DEPTH", cs%minimum_forcing_depth, &
3344 "The smallest depth over which forcing can be applied. This "//&
3345 "only takes effect when near-surface layers become thin "//&
3346 "relative to this scale, in which case the forcing tendencies "//&
3347 "scaled down by distributing the forcing over this depth scale.", &
3348 units=
"m", default=0.001)
3349 call get_param(param_file, mdl,
"EVAP_CFL_LIMIT", cs%evap_CFL_limit, &
3350 "The largest fraction of a layer than can be lost to forcing "//&
3351 "(e.g. evaporation, sea-ice formation) in one time-step. The unused "//&
3352 "mass loss is passed down through the column.", &
3353 units=
"nondim", default=0.8)
3357 if (gv%Boussinesq)
then ; thickness_units =
"m"
3358 else ; thickness_units =
"kg m-2" ;
endif
3360 cs%id_ea_t = register_diag_field(
'ocean_model',
'ea_t',diag%axesTL,time, &
3361 'Layer (heat) entrainment from above per timestep',
'm')
3362 cs%id_eb_t = register_diag_field(
'ocean_model',
'eb_t',diag%axesTL,time, &
3363 'Layer (heat) entrainment from below per timestep',
'm')
3364 cs%id_ea_s = register_diag_field(
'ocean_model',
'ea_s',diag%axesTL,time, &
3365 'Layer (salt) entrainment from above per timestep',
'm')
3366 cs%id_eb_s = register_diag_field(
'ocean_model',
'eb_s',diag%axesTL,time, &
3367 'Layer (salt) entrainment from below per timestep',
'm')
3369 cs%id_ea = register_diag_field(
'ocean_model',
'ea',diag%axesTL,time, &
3370 'Layer entrainment from above per timestep',
'm')
3371 cs%id_eb = register_diag_field(
'ocean_model',
'eb',diag%axesTL,time, &
3372 'Layer entrainment from below per timestep',
'm')
3373 cs%id_wd = register_diag_field(
'ocean_model',
'wd',diag%axesTi,time, &
3374 'Diapycnal velocity',
'm s-1')
3375 if (cs%id_wd > 0)
call safe_alloc_ptr(cdp%diapyc_vel,isd,ied,jsd,jed,nz+1)
3377 cs%id_dudt_dia = register_diag_field(
'ocean_model',
'dudt_dia',diag%axesCuL,time, &
3378 'Zonal Acceleration from Diapycnal Mixing',
'm s-2')
3379 cs%id_dvdt_dia = register_diag_field(
'ocean_model',
'dvdt_dia',diag%axesCvL,time, &
3380 'Meridional Acceleration from Diapycnal Mixing',
'm s-2')
3382 if (cs%use_int_tides)
then
3383 cs%id_cg1 = register_diag_field(
'ocean_model',
'cn1', diag%axesT1, &
3384 time,
'First baroclinic mode (eigen) speed',
'm s-1')
3385 allocate(cs%id_cn(cs%nMode)) ; cs%id_cn(:) = -1
3387 write(var_name,
'("cn_mode",i1)') m
3388 write(var_descript,
'("Baroclinic (eigen) speed of mode ",i1)') m
3389 cs%id_cn(m) = register_diag_field(
'ocean_model',var_name, diag%axesT1, &
3390 time, var_descript,
'm s-1')
3391 call mom_mesg(
"Registering "//trim(var_name)//
", Described as: "//var_descript, 5)
3395 if (use_temperature)
then
3396 cs%id_Tdif = register_diag_field(
'ocean_model',
"Tflx_dia_diff",diag%axesTi, &
3397 time,
"Diffusive diapycnal temperature flux across interfaces", &
3399 cs%id_Tadv = register_diag_field(
'ocean_model',
"Tflx_dia_adv",diag%axesTi, &
3400 time,
"Advective diapycnal temperature flux across interfaces", &
3402 cs%id_Sdif = register_diag_field(
'ocean_model',
"Sflx_dia_diff",diag%axesTi, &
3403 time,
"Diffusive diapycnal salnity flux across interfaces", &
3405 cs%id_Sadv = register_diag_field(
'ocean_model',
"Sflx_dia_adv",diag%axesTi, &
3406 time,
"Advective diapycnal salnity flux across interfaces", &
3408 cs%id_MLD_003 = register_diag_field(
'ocean_model',
'MLD_003', diag%axesT1, time, &
3409 'Mixed layer depth (delta rho = 0.03)',
'm', conversion=us%Z_to_m, &
3410 cmor_field_name=
'mlotst', cmor_long_name=
'Ocean Mixed Layer Thickness Defined by Sigma T', &
3411 cmor_standard_name=
'ocean_mixed_layer_thickness_defined_by_sigma_t')
3412 cs%id_mlotstsq = register_diag_field(
'ocean_model',
'mlotstsq',diag%axesT1, time, &
3413 long_name=
'Square of Ocean Mixed Layer Thickness Defined by Sigma T', &
3414 standard_name=
'square_of_ocean_mixed_layer_thickness_defined_by_sigma_t', &
3415 units=
'm2', conversion=us%Z_to_m**2)
3416 cs%id_MLD_0125 = register_diag_field(
'ocean_model',
'MLD_0125',diag%axesT1,time, &
3417 'Mixed layer depth (delta rho = 0.125)',
'm', conversion=us%Z_to_m)
3418 cs%id_subMLN2 = register_diag_field(
'ocean_model',
'subML_N2',diag%axesT1,time, &
3419 'Squared buoyancy frequency below mixed layer',
's-2', conversion=us%s_to_T**2)
3420 cs%id_MLD_user = register_diag_field(
'ocean_model',
'MLD_user',diag%axesT1,time, &
3421 'Mixed layer depth (used defined)',
'm', conversion=us%Z_to_m)
3423 call get_param(param_file, mdl,
"DIAG_MLD_DENSITY_DIFF", cs%MLDdensityDifference, &
3424 "The density difference used to determine a diagnostic mixed "//&
3425 "layer depth, MLD_user, following the definition of Levitus 1982. "//&
3426 "The MLD is the depth at which the density is larger than the "//&
3427 "surface density by the specified amount.", units=
'kg/m3', default=0.1)
3428 call get_param(param_file, mdl,
"DIAG_DEPTH_SUBML_N2", cs%dz_subML_N2, &
3429 "The distance over which to calculate a diagnostic of the "//&
3430 "stratification at the base of the mixed layer.", &
3431 units=
'm', default=50.0, scale=us%m_to_Z)
3433 if (cs%id_dudt_dia > 0)
call safe_alloc_ptr(adp%du_dt_dia,isdb,iedb,jsd,jed,nz)
3434 if (cs%id_dvdt_dia > 0)
call safe_alloc_ptr(adp%dv_dt_dia,isd,ied,jsdb,jedb,nz)
3437 cs%id_u_predia = register_diag_field(
'ocean_model',
'u_predia', diag%axesCuL, time, &
3438 'Zonal velocity before diabatic forcing',
'm s-1')
3439 cs%id_v_predia = register_diag_field(
'ocean_model',
'v_predia', diag%axesCvL, time, &
3440 'Meridional velocity before diabatic forcing',
'm s-1')
3441 cs%id_h_predia = register_diag_field(
'ocean_model',
'h_predia', diag%axesTL, time, &
3442 'Layer Thickness before diabatic forcing', thickness_units, v_extensive=.true.)
3443 cs%id_e_predia = register_diag_field(
'ocean_model',
'e_predia', diag%axesTi, time, &
3444 'Interface Heights before diabatic forcing',
'm')
3445 if (use_temperature)
then
3446 cs%id_T_predia = register_diag_field(
'ocean_model',
'temp_predia', diag%axesTL, time, &
3447 'Potential Temperature',
'degC')
3448 cs%id_S_predia = register_diag_field(
'ocean_model',
'salt_predia', diag%axesTL, time, &
3454 cs%id_Kd_interface = register_diag_field(
'ocean_model',
'Kd_interface', diag%axesTi, time, &
3455 'Total diapycnal diffusivity at interfaces',
'm2 s-1', conversion=us%Z2_T_to_m2_s)
3456 if (cs%use_energetic_PBL)
then
3457 cs%id_Kd_ePBL = register_diag_field(
'ocean_model',
'Kd_ePBL', diag%axesTi, time, &
3458 'ePBL diapycnal diffusivity at interfaces',
'm2 s-1', conversion=us%Z2_T_to_m2_s)
3461 cs%id_Kd_heat = register_diag_field(
'ocean_model',
'Kd_heat', diag%axesTi, time, &
3462 'Total diapycnal diffusivity for heat at interfaces',
'm2 s-1', conversion=us%Z2_T_to_m2_s, &
3463 cmor_field_name=
'difvho', &
3464 cmor_standard_name=
'ocean_vertical_heat_diffusivity', &
3465 cmor_long_name=
'Ocean vertical heat diffusivity')
3466 cs%id_Kd_salt = register_diag_field(
'ocean_model',
'Kd_salt', diag%axesTi, time, &
3467 'Total diapycnal diffusivity for salt at interfaces',
'm2 s-1', conversion=us%Z2_T_to_m2_s, &
3468 cmor_field_name=
'difvso', &
3469 cmor_standard_name=
'ocean_vertical_salt_diffusivity', &
3470 cmor_long_name=
'Ocean vertical salt diffusivity')
3474 cs%useKPP = kpp_init(param_file, g, gv, us, diag, time, cs%KPP_CSp, passive=cs%KPPisPassive)
3476 allocate( cs%KPP_NLTheat(isd:ied,jsd:jed,nz+1) ) ; cs%KPP_NLTheat(:,:,:) = 0.
3477 allocate( cs%KPP_NLTscalar(isd:ied,jsd:jed,nz+1) ) ; cs%KPP_NLTscalar(:,:,:) = 0.
3480 allocate( cs%KPP_buoy_flux(isd:ied,jsd:jed,nz+1) ) ; cs%KPP_buoy_flux(:,:,:) = 0.
3481 allocate( cs%KPP_temp_flux(isd:ied,jsd:jed) ) ; cs%KPP_temp_flux(:,:) = 0.
3482 allocate( cs%KPP_salt_flux(isd:ied,jsd:jed) ) ; cs%KPP_salt_flux(:,:) = 0.
3485 if (cs%useKPP .and. differentialdiffusion)
then
3486 call mom_error(fatal,
'diabatic_driver_init: '// &
3487 'DOUBLE_DIFFUSION (old method) does not work with KPP. Please'//&
3488 'set DOUBLE_DIFFUSION=False and USE_CVMIX_DDIFF=True.')
3491 call get_param(param_file, mdl,
"SALT_REJECT_BELOW_ML", cs%salt_reject_below_ML, &
3492 "If true, place salt from brine rejection below the mixed layer, "// &
3493 "into the first non-vanished layer for which the column remains stable", &
3496 if (cs%salt_reject_below_ML)
then
3497 cs%id_brine_lay = register_diag_field(
'ocean_model',
'brine_layer',diag%axesT1,time, &
3498 'Brine insertion layer',
'none')
3505 cs%id_diabatic_diff_h = register_diag_field(
'ocean_model',
'diabatic_diff_h', diag%axesTL, time, &
3506 long_name =
'Cell thickness used during diabatic diffusion', units=
'm', v_extensive=.true.)
3507 if (cs%useALEalgorithm)
then
3508 cs%id_diabatic_diff_temp_tend = register_diag_field(
'ocean_model', &
3509 'diabatic_diff_temp_tendency', diag%axesTL, time, &
3510 'Diabatic diffusion temperature tendency',
'degC s-1')
3511 if (cs%id_diabatic_diff_temp_tend > 0)
then
3512 cs%diabatic_diff_tendency_diag = .true.
3515 cs%id_diabatic_diff_saln_tend = register_diag_field(
'ocean_model',&
3516 'diabatic_diff_saln_tendency', diag%axesTL, time, &
3517 'Diabatic diffusion salinity tendency',
'psu s-1')
3518 if (cs%id_diabatic_diff_saln_tend > 0)
then
3519 cs%diabatic_diff_tendency_diag = .true.
3522 cs%id_diabatic_diff_heat_tend = register_diag_field(
'ocean_model', &
3523 'diabatic_heat_tendency', diag%axesTL, time, &
3524 'Diabatic diffusion heat tendency', &
3525 'W m-2',cmor_field_name=
'opottempdiff', &
3526 cmor_standard_name=
'tendency_of_sea_water_potential_temperature_expressed_as_heat_content_'// &
3527 'due_to_parameterized_dianeutral_mixing', &
3528 cmor_long_name=
'Tendency of sea water potential temperature expressed as heat content '// &
3529 'due to parameterized dianeutral mixing',&
3531 if (cs%id_diabatic_diff_heat_tend > 0)
then
3532 cs%diabatic_diff_tendency_diag = .true.
3535 cs%id_diabatic_diff_salt_tend = register_diag_field(
'ocean_model', &
3536 'diabatic_salt_tendency', diag%axesTL, time, &
3537 'Diabatic diffusion of salt tendency', &
3538 'kg m-2 s-1',cmor_field_name=
'osaltdiff', &
3539 cmor_standard_name=
'tendency_of_sea_water_salinity_expressed_as_salt_content_'// &
3540 'due_to_parameterized_dianeutral_mixing', &
3541 cmor_long_name=
'Tendency of sea water salinity expressed as salt content '// &
3542 'due to parameterized dianeutral mixing', &
3544 if (cs%id_diabatic_diff_salt_tend > 0)
then
3545 cs%diabatic_diff_tendency_diag = .true.
3549 cs%id_diabatic_diff_heat_tend_2d = register_diag_field(
'ocean_model', &
3550 'diabatic_heat_tendency_2d', diag%axesT1, time, &
3551 'Depth integrated diabatic diffusion heat tendency', &
3552 'W m-2',cmor_field_name=
'opottempdiff_2d', &
3553 cmor_standard_name=
'tendency_of_sea_water_potential_temperature_expressed_as_heat_content_'//&
3554 'due_to_parameterized_dianeutral_mixing_depth_integrated', &
3555 cmor_long_name=
'Tendency of sea water potential temperature expressed as heat content '//&
3556 'due to parameterized dianeutral mixing depth integrated')
3557 if (cs%id_diabatic_diff_heat_tend_2d > 0)
then
3558 cs%diabatic_diff_tendency_diag = .true.
3562 cs%id_diabatic_diff_salt_tend_2d = register_diag_field(
'ocean_model', &
3563 'diabatic_salt_tendency_2d', diag%axesT1, time, &
3564 'Depth integrated diabatic diffusion salt tendency', &
3565 'kg m-2 s-1',cmor_field_name=
'osaltdiff_2d', &
3566 cmor_standard_name=
'tendency_of_sea_water_salinity_expressed_as_salt_content_'// &
3567 'due_to_parameterized_dianeutral_mixing_depth_integrated', &
3568 cmor_long_name=
'Tendency of sea water salinity expressed as salt content '// &
3569 'due to parameterized dianeutral mixing depth integrated')
3570 if (cs%id_diabatic_diff_salt_tend_2d > 0)
then
3571 cs%diabatic_diff_tendency_diag = .true.
3577 cs%id_boundary_forcing_h = register_diag_field(
'ocean_model',
'boundary_forcing_h', diag%axesTL, time, &
3578 long_name =
'Cell thickness after applying boundary forcing', units=
'm', v_extensive=.true.)
3579 cs%id_boundary_forcing_h_tendency = register_diag_field(
'ocean_model', &
3580 'boundary_forcing_h_tendency', diag%axesTL, time, &
3581 'Cell thickness tendency due to boundary forcing',
'm s-1', &
3582 v_extensive = .true.)
3583 if (cs%id_boundary_forcing_h_tendency > 0)
then
3584 cs%boundary_forcing_tendency_diag = .true.
3587 cs%id_boundary_forcing_temp_tend = register_diag_field(
'ocean_model',&
3588 'boundary_forcing_temp_tendency', diag%axesTL, time, &
3589 'Boundary forcing temperature tendency',
'degC s-1')
3590 if (cs%id_boundary_forcing_temp_tend > 0)
then
3591 cs%boundary_forcing_tendency_diag = .true.
3594 cs%id_boundary_forcing_saln_tend = register_diag_field(
'ocean_model',&
3595 'boundary_forcing_saln_tendency', diag%axesTL, time, &
3596 'Boundary forcing saln tendency',
'psu s-1')
3597 if (cs%id_boundary_forcing_saln_tend > 0)
then
3598 cs%boundary_forcing_tendency_diag = .true.
3601 cs%id_boundary_forcing_heat_tend = register_diag_field(
'ocean_model',&
3602 'boundary_forcing_heat_tendency', diag%axesTL, time, &
3603 'Boundary forcing heat tendency',
'W m-2', &
3604 v_extensive = .true.)
3605 if (cs%id_boundary_forcing_heat_tend > 0)
then
3606 cs%boundary_forcing_tendency_diag = .true.
3609 cs%id_boundary_forcing_salt_tend = register_diag_field(
'ocean_model',&
3610 'boundary_forcing_salt_tendency', diag%axesTL, time, &
3611 'Boundary forcing salt tendency',
'kg m-2 s-1', &
3612 v_extensive = .true.)
3613 if (cs%id_boundary_forcing_salt_tend > 0)
then
3614 cs%boundary_forcing_tendency_diag = .true.
3618 cs%id_boundary_forcing_heat_tend_2d = register_diag_field(
'ocean_model',&
3619 'boundary_forcing_heat_tendency_2d', diag%axesT1, time, &
3620 'Depth integrated boundary forcing of ocean heat',
'W m-2')
3621 if (cs%id_boundary_forcing_heat_tend_2d > 0)
then
3622 cs%boundary_forcing_tendency_diag = .true.
3626 cs%id_boundary_forcing_salt_tend_2d = register_diag_field(
'ocean_model',&
3627 'boundary_forcing_salt_tendency_2d', diag%axesT1, time, &
3628 'Depth integrated boundary forcing of ocean salt',
'kg m-2 s-1')
3629 if (cs%id_boundary_forcing_salt_tend_2d > 0)
then
3630 cs%boundary_forcing_tendency_diag = .true.
3635 cs%id_frazil_h = register_diag_field(
'ocean_model',
'frazil_h', diag%axesTL, time, &
3636 long_name =
'Cell Thickness', standard_name=
'cell_thickness', units=
'm', v_extensive=.true.)
3639 cs%id_frazil_temp_tend = register_diag_field(
'ocean_model',&
3640 'frazil_temp_tendency', diag%axesTL, time, &
3641 'Temperature tendency due to frazil formation',
'degC s-1')
3642 if (cs%id_frazil_temp_tend > 0)
then
3643 cs%frazil_tendency_diag = .true.
3647 cs%id_frazil_heat_tend = register_diag_field(
'ocean_model',&
3648 'frazil_heat_tendency', diag%axesTL, time, &
3649 'Heat tendency due to frazil formation',
'W m-2', v_extensive = .true.)
3650 if (cs%id_frazil_heat_tend > 0)
then
3651 cs%frazil_tendency_diag = .true.
3655 cs%id_frazil_heat_tend_2d = register_diag_field(
'ocean_model',&
3656 'frazil_heat_tendency_2d', diag%axesT1, time, &
3657 'Depth integrated heat tendency due to frazil formation',
'W m-2')
3658 if (cs%id_frazil_heat_tend_2d > 0)
then
3659 cs%frazil_tendency_diag = .true.
3662 if (cs%frazil_tendency_diag)
then
3663 allocate(cs%frazil_temp_diag(isd:ied,jsd:jed,nz) ) ; cs%frazil_temp_diag(:,:,:) = 0.
3664 allocate(cs%frazil_heat_diag(isd:ied,jsd:jed,nz) ) ; cs%frazil_heat_diag(:,:,:) = 0.
3668 cs%use_tidal_mixing = tidal_mixing_init(time, g, gv, us, param_file, diag, &
3669 cs%tidal_mixing_CSp)
3673 cs%use_CVMix_conv = cvmix_conv_init(time, g, gv, us, param_file, diag, cs%CVMix_conv_csp)
3675 call entrain_diffusive_init(time, g, gv, us, param_file, diag, cs%entrain_diffusive_CSp)
3678 if (cs%use_geothermal) &
3679 call geothermal_init(time, g, param_file, diag, cs%geothermal_CSp)
3682 if (cs%use_int_tides)
then
3683 call int_tide_input_init(time, g, gv, us, param_file, diag, cs%int_tide_input_CSp, &
3685 call internal_tides_init(time, g, gv, us, param_file, diag, cs%int_tide_CSp)
3689 call set_diffusivity_init(time, g, gv, us, param_file, diag, cs%set_diff_CSp, &
3690 cs%int_tide_CSp, cs%tidal_mixing_CSp, cs%halo_TS_diff)
3693 id_clock_entrain = cpu_clock_id(
'(Ocean diabatic entrain)', grain=clock_module)
3694 if (cs%bulkmixedlayer) &
3695 id_clock_mixedlayer = cpu_clock_id(
'(Ocean mixed layer)', grain=clock_module)
3696 id_clock_remap = cpu_clock_id(
'(Ocean vert remap)', grain=clock_module)
3697 if (cs%use_geothermal) &
3698 id_clock_geothermal = cpu_clock_id(
'(Ocean geothermal)', grain=clock_routine)
3699 id_clock_set_diffusivity = cpu_clock_id(
'(Ocean set_diffusivity)', grain=clock_module)
3700 id_clock_kpp = cpu_clock_id(
'(Ocean KPP)', grain=clock_module)
3701 id_clock_tracers = cpu_clock_id(
'(Ocean tracer_columns)', grain=clock_module_driver+5)
3702 if (cs%use_sponge) &
3703 id_clock_sponge = cpu_clock_id(
'(Ocean sponges)', grain=clock_module)
3704 id_clock_tridiag = cpu_clock_id(
'(Ocean diabatic tridiag)', grain=clock_routine)
3705 id_clock_pass = cpu_clock_id(
'(Ocean diabatic message passing)', grain=clock_routine)
3706 id_clock_differential_diff = -1 ;
if (differentialdiffusion) &
3707 id_clock_differential_diff = cpu_clock_id(
'(Ocean differential diffusion)', grain=clock_routine)
3710 call diabatic_aux_init(time, g, gv, us, param_file, diag, cs%diabatic_aux_CSp, &
3711 cs%useALEalgorithm, cs%use_energetic_PBL)
3714 if (cs%bulkmixedlayer) &
3715 call bulkmixedlayer_init(time, g, gv, us, param_file, diag, cs%bulkmixedlayer_CSp)
3716 if (cs%use_energetic_PBL) &
3717 call energetic_pbl_init(time, g, gv, us, param_file, diag, cs%energetic_PBL_CSp)
3719 call regularize_layers_init(time, g, gv, param_file, diag, cs%regularize_layers_CSp)
3721 if (cs%debug_energy_req) &
3722 call diapyc_energy_req_init(time, g, gv, us, param_file, diag, cs%diapyc_en_rec_CSp)
3725 if (use_temperature)
then
3726 call get_param(param_file, mdl,
"PEN_SW_NBANDS", nbands, default=1)
3727 if (nbands > 0)
then
3729 call opacity_init(time, g, gv, us, param_file, diag, cs%opacity_CSp, cs%optics)
3734 call diag_grid_storage_init(cs%diag_grids_prev, g, diag)
3736 end subroutine diabatic_driver_init
3740 subroutine diabatic_driver_end(CS)
3743 if (.not.
associated(cs))
return
3745 call diabatic_aux_end(cs%diabatic_aux_CSp)
3747 call entrain_diffusive_end(cs%entrain_diffusive_CSp)
3748 call set_diffusivity_end(cs%set_diff_CSp)
3751 deallocate( cs%KPP_buoy_flux )
3752 deallocate( cs%KPP_temp_flux )
3753 deallocate( cs%KPP_salt_flux )
3756 deallocate( cs%KPP_NLTheat )
3757 deallocate( cs%KPP_NLTscalar )
3758 call kpp_end(cs%KPP_CSp)
3761 if (cs%use_tidal_mixing)
call tidal_mixing_end(cs%tidal_mixing_CSp)
3763 if (cs%use_CVMix_conv)
call cvmix_conv_end(cs%CVMix_conv_csp)
3765 if (cs%use_energetic_PBL) &
3766 call energetic_pbl_end(cs%energetic_PBL_CSp)
3767 if (cs%debug_energy_req) &
3768 call diapyc_energy_req_end(cs%diapyc_en_rec_CSp)
3770 if (
associated(cs%optics))
then
3771 call opacity_end(cs%opacity_CSp, cs%optics)
3772 deallocate(cs%optics)
3783 end subroutine diabatic_driver_end