6 use mom_cpu_clock,
only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, clock_routine
22 implicit none ;
private
24 #include <MOM_memory.h>
26 public energetic_pbl, energetic_pbl_init, energetic_pbl_end
27 public energetic_pbl_get_mld
52 logical :: use_mld_iteration=.false.
53 logical :: mld_iteration_guess=.false.
55 integer :: max_mld_its
57 real :: mixlenexponent
60 real :: mke_to_tke_effic
65 real :: ekman_scale_coef
68 real :: translay_scale
82 real :: wstar_ustar_coef
86 real :: vstar_surf_fac
88 real :: vstar_scale_fac
93 integer :: mstar_scheme
94 logical :: mstar_flatcap=.true.
112 real :: mstar_coef = 0.3
115 real :: rh18_mstar_cn1
119 real :: rh18_mstar_cn2
123 real :: rh18_mstar_cn3
126 real :: rh18_mstar_cs1
128 real :: rh18_mstar_cs2
132 real :: mstar_convect_coef
135 logical :: use_lt = .false.
136 integer :: lt_enhance_form
137 real :: lt_enhance_coef
138 real :: lt_enhance_exp
141 real :: lac_mldoob_stab
143 real :: lac_ekoob_stab
145 real :: lac_mldoob_un
149 real :: max_enhance_m = 5.
152 type(time_type),
pointer :: time=>null()
154 logical :: tke_diagnostics = .false.
155 logical :: answers_2018
158 logical :: orig_pe_calc
165 real,
allocatable,
dimension(:,:) :: &
169 real,
allocatable,
dimension(:,:) :: &
170 diag_tke_wind, & !< The wind source of TKE [kg m-3 Z3 T-3 ~> W m-2].
171 diag_tke_mke, & !< The resolved KE source of TKE [kg m-3 Z3 T-3 ~> W m-2].
172 diag_tke_conv, & !< The convective source of TKE [kg m-3 Z3 T-3 ~> W m-2].
173 diag_tke_forcing, & !< The TKE sink required to mix surface penetrating shortwave heating
175 diag_tke_mech_decay, &
176 diag_tke_conv_decay, &
184 real,
allocatable,
dimension(:,:,:) :: &
185 velocity_scale, & !< The velocity scale used in getting Kd [Z T-1 ~> m s-1]
188 integer :: id_ml_depth = -1, id_tke_wind = -1, id_tke_mixing = -1
189 integer :: id_tke_mke = -1, id_tke_conv = -1, id_tke_forcing = -1
190 integer :: id_tke_mech_decay = -1, id_tke_conv_decay = -1
191 integer :: id_mixing_length = -1, id_velocity_scale = -1
192 integer :: id_mstar_mix = -1, id_la_mod = -1, id_la = -1, id_mstar_lt = -1
197 integer,
parameter :: use_fixed_mstar = 0
198 integer,
parameter :: mstar_from_ekman = 2
200 integer,
parameter :: mstar_from_rh18 = 3
201 integer,
parameter :: no_langmuir = 0
202 integer,
parameter :: langmuir_rescale = 2
204 integer,
parameter :: langmuir_add = 3
206 integer,
parameter :: wt_from_croot_tke = 0
208 integer,
parameter :: wt_from_rh18 = 1
211 character*(20),
parameter :: constant_string =
"CONSTANT"
212 character*(20),
parameter :: om4_string =
"OM4"
213 character*(20),
parameter :: rh18_string =
"REICHL_H18"
214 character*(20),
parameter :: root_tke_string =
"CUBE_ROOT_TKE"
215 character*(20),
parameter :: none_string =
"NONE"
216 character*(20),
parameter :: rescaled_string =
"RESCALE"
217 character*(20),
parameter :: additive_string =
"ADDITIVE"
223 real :: dtke_conv, dtke_forcing, dtke_wind, dtke_mixing
224 real :: dtke_mke, dtke_mech_decay, dtke_conv_decay
230 real,
allocatable,
dimension(:) :: dt_expect
231 real,
allocatable,
dimension(:) :: ds_expect
240 subroutine energetic_pbl(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS, &
241 dSV_dT, dSV_dS, TKE_forced, buoy_flux, dt_diag, last_call, &
242 dT_expected, dS_expected, Waves )
246 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
247 intent(inout) :: h_3d
248 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
251 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
254 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
258 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
261 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
262 intent(in) :: tke_forced
268 type(
forcing),
intent(inout) :: fluxes
271 real,
intent(in) :: dt
272 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
273 intent(out) :: kd_int
277 real,
dimension(SZI_(G),SZJ_(G)), &
278 intent(in) :: buoy_flux
279 real,
optional,
intent(in) :: dt_diag
281 logical,
optional,
intent(in) :: last_call
285 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
286 optional,
intent(out) :: dt_expected
289 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
290 optional,
intent(out) :: ds_expected
294 optional,
pointer :: waves
319 real,
dimension(SZI_(G),SZK_(GV)) :: &
328 real,
dimension(SZI_(G),SZK_(GV)+1) :: &
330 real,
dimension(SZK_(GV)) :: &
339 real,
dimension(SZK_(GV)+1) :: &
354 logical :: write_diags
355 logical :: reset_diags
357 logical :: debug=.false.
360 integer :: i, j, k, is, ie, js, je, nz
362 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
364 if (.not.
associated(cs))
call mom_error(fatal,
"energetic_PBL: "//&
365 "Module must be initialized before it is used.")
367 if (.not.
associated(tv%eqn_of_state))
call mom_error(fatal, &
368 "energetic_PBL: Temperature, salinity and an equation of state "//&
370 if (.NOT.
associated(fluxes%ustar))
call mom_error(fatal, &
371 "energetic_PBL: No surface TKE fluxes (ustar) defined in mixedlayer!")
372 debug = .false. ;
if (
present(dt_expected) .or.
present(ds_expected)) debug = .true.
374 if (debug)
allocate(ecd%dT_expect(nz), ecd%dS_expect(nz))
376 h_neglect = gv%H_subroundoff
378 dt__diag = dt ;
if (
present(dt_diag)) dt__diag = dt_diag
379 write_diags = .true. ;
if (
present(last_call)) write_diags = last_call
384 if (
present(dt_diag) .and. write_diags .and. (dt__diag > dt)) &
385 reset_diags = .false.
387 if (reset_diags)
then
388 if (cs%TKE_diagnostics)
then
390 do j=js,je ;
do i=is,ie
391 cs%diag_TKE_wind(i,j) = 0.0 ; cs%diag_TKE_MKE(i,j) = 0.0
392 cs%diag_TKE_conv(i,j) = 0.0 ; cs%diag_TKE_forcing(i,j) = 0.0
393 cs%diag_TKE_mixing(i,j) = 0.0 ; cs%diag_TKE_mech_decay(i,j) = 0.0
394 cs%diag_TKE_conv_decay(i,j) = 0.0
406 do k=1,nz ;
do i=is,ie
407 h_2d(i,k) = h_3d(i,j,k) ; u_2d(i,k) = u_3d(i,j,k) ; v_2d(i,k) = v_3d(i,j,k)
408 t_2d(i,k) = tv%T(i,j,k) ; s_2d(i,k) = tv%S(i,j,k)
409 tke_forced_2d(i,k) = tke_forced(i,j,k)
410 dsv_dt_2d(i,k) = dsv_dt(i,j,k) ; dsv_ds_2d(i,k) = dsv_ds(i,j,k)
419 do i=is,ie ;
if (g%mask2dT(i,j) > 0.5)
then
423 h(k) = h_2d(i,k) + gv%H_subroundoff ; u(k) = u_2d(i,k) ; v(k) = v_2d(i,k)
424 t0(k) = t_2d(i,k) ; s0(k) = s_2d(i,k) ; tke_forcing(k) = tke_forced_2d(i,k)
425 dsv_dt_1d(k) = dsv_dt_2d(i,k) ; dsv_ds_1d(k) = dsv_ds_2d(i,k)
427 do k=1,nz+1 ; kd(k) = 0.0 ;
enddo
430 u_star = fluxes%ustar(i,j)
431 u_star_mean = fluxes%ustar_gustless(i,j)
432 b_flux = buoy_flux(i,j)
433 if (
associated(fluxes%ustar_shelf) .and.
associated(fluxes%frac_shelf_h))
then
434 if (fluxes%frac_shelf_h(i,j) > 0.0) &
435 u_star = (1.0 - fluxes%frac_shelf_h(i,j)) * u_star + &
436 fluxes%frac_shelf_h(i,j) * fluxes%ustar_shelf(i,j)
438 if (u_star < cs%ustar_min) u_star = cs%ustar_min
439 if (cs%omega_frac >= 1.0)
then
442 absf = 0.25*((abs(g%CoriolisBu(i,j)) + abs(g%CoriolisBu(i-1,j-1))) + &
443 (abs(g%CoriolisBu(i,j-1)) + abs(g%CoriolisBu(i-1,j))))
444 if (cs%omega_frac > 0.0) &
445 absf = sqrt(cs%omega_frac*4.0*cs%omega**2 + (1.0-cs%omega_frac)*absf**2)
450 if (cs%MLD_iteration_guess .and. (cs%ML_Depth(i,j) > 0.0)) mld_io = cs%ML_Depth(i,j)
452 call epbl_column(h, u, v, t0, s0, dsv_dt_1d, dsv_ds_1d, tke_forcing, b_flux, absf, &
453 u_star, u_star_mean, dt, mld_io, kd, mixvel, mixlen, gv, &
454 us, cs, ecd, dt_diag=dt_diag, waves=waves, g=g, i=i, j=j)
461 cs%ML_depth(i,j) = mld_io
463 if (
present(dt_expected))
then
464 do k=1,nz ; dt_expected(i,j,k) = ecd%dT_expect(k) ;
enddo
466 if (
present(ds_expected))
then
467 do k=1,nz ; ds_expected(i,j,k) = ecd%dS_expect(k) ;
enddo
470 if (cs%TKE_diagnostics)
then
471 cs%diag_TKE_MKE(i,j) = cs%diag_TKE_MKE(i,j) + ecd%dTKE_MKE
472 cs%diag_TKE_conv(i,j) = cs%diag_TKE_conv(i,j) + ecd%dTKE_conv
473 cs%diag_TKE_forcing(i,j) = cs%diag_TKE_forcing(i,j) + ecd%dTKE_forcing
474 cs%diag_TKE_wind(i,j) = cs%diag_TKE_wind(i,j) + ecd%dTKE_wind
475 cs%diag_TKE_mixing(i,j) = cs%diag_TKE_mixing(i,j) + ecd%dTKE_mixing
476 cs%diag_TKE_mech_decay(i,j) = cs%diag_TKE_mech_decay(i,j) + ecd%dTKE_mech_decay
477 cs%diag_TKE_conv_decay(i,j) = cs%diag_TKE_conv_decay(i,j) + ecd%dTKE_conv_decay
481 if (cs%id_Mixing_Length>0)
then ;
do k=1,nz
482 cs%Mixing_Length(i,j,k) = mixlen(k)
484 if (cs%id_Velocity_Scale>0)
then ;
do k=1,nz
485 cs%Velocity_Scale(i,j,k) = mixvel(k)
487 if (
allocated(cs%mstar_mix)) cs%mstar_mix(i,j) = ecd%mstar
488 if (
allocated(cs%mstar_lt)) cs%mstar_lt(i,j) = ecd%mstar_LT
489 if (
allocated(cs%La)) cs%La(i,j) = ecd%LA
490 if (
allocated(cs%La_mod)) cs%La_mod(i,j) = ecd%LAmod
493 do k=1,nz+1 ; kd_2d(i,k) = 0. ;
enddo
494 cs%ML_depth(i,j) = 0.0
496 if (
present(dt_expected))
then
497 do k=1,nz ; dt_expected(i,j,k) = 0.0 ;
enddo
499 if (
present(ds_expected))
then
500 do k=1,nz ; ds_expected(i,j,k) = 0.0 ;
enddo
504 do k=1,nz+1 ;
do i=is,ie ; kd_int(i,j,k) = kd_2d(i,k) ;
enddo ;
enddo
508 if (write_diags)
then
509 if (cs%id_ML_depth > 0)
call post_data(cs%id_ML_depth, cs%ML_depth, cs%diag)
510 if (cs%id_TKE_wind > 0)
call post_data(cs%id_TKE_wind, cs%diag_TKE_wind, cs%diag)
511 if (cs%id_TKE_MKE > 0)
call post_data(cs%id_TKE_MKE, cs%diag_TKE_MKE, cs%diag)
512 if (cs%id_TKE_conv > 0)
call post_data(cs%id_TKE_conv, cs%diag_TKE_conv, cs%diag)
513 if (cs%id_TKE_forcing > 0)
call post_data(cs%id_TKE_forcing, cs%diag_TKE_forcing, cs%diag)
514 if (cs%id_TKE_mixing > 0)
call post_data(cs%id_TKE_mixing, cs%diag_TKE_mixing, cs%diag)
515 if (cs%id_TKE_mech_decay > 0) &
516 call post_data(cs%id_TKE_mech_decay, cs%diag_TKE_mech_decay, cs%diag)
517 if (cs%id_TKE_conv_decay > 0) &
518 call post_data(cs%id_TKE_conv_decay, cs%diag_TKE_conv_decay, cs%diag)
519 if (cs%id_Mixing_Length > 0)
call post_data(cs%id_Mixing_Length, cs%Mixing_Length, cs%diag)
520 if (cs%id_Velocity_Scale >0)
call post_data(cs%id_Velocity_Scale, cs%Velocity_Scale, cs%diag)
521 if (cs%id_MSTAR_MIX > 0)
call post_data(cs%id_MSTAR_MIX, cs%MSTAR_MIX, cs%diag)
522 if (cs%id_LA > 0)
call post_data(cs%id_LA, cs%LA, cs%diag)
523 if (cs%id_LA_MOD > 0)
call post_data(cs%id_LA_MOD, cs%LA_MOD, cs%diag)
524 if (cs%id_MSTAR_LT > 0)
call post_data(cs%id_MSTAR_LT, cs%MSTAR_LT, cs%diag)
527 if (debug)
deallocate(ecd%dT_expect, ecd%dS_expect)
529 end subroutine energetic_pbl
535 subroutine epbl_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, absf, &
536 u_star, u_star_mean, dt, MLD_io, Kd, mixvel, mixlen, GV, US, CS, eCD, &
537 dt_diag, Waves, G, i, j)
540 real,
dimension(SZK_(GV)),
intent(in) :: h
541 real,
dimension(SZK_(GV)),
intent(in) :: u
543 real,
dimension(SZK_(GV)),
intent(in) :: v
545 real,
dimension(SZK_(GV)),
intent(in) :: T0
546 real,
dimension(SZK_(GV)),
intent(in) :: S0
548 real,
dimension(SZK_(GV)),
intent(in) :: dSV_dT
551 real,
dimension(SZK_(GV)),
intent(in) :: dSV_dS
553 real,
dimension(SZK_(GV)),
intent(in) :: TKE_forcing
556 real,
intent(in) :: B_flux
557 real,
intent(in) :: absf
558 real,
intent(in) :: u_star
559 real,
intent(in) :: u_star_mean
561 real,
intent(inout) :: MLD_io
563 real,
intent(in) :: dt
564 real,
dimension(SZK_(GV)+1), &
567 real,
dimension(SZK_(GV)+1), &
568 intent(out) :: mixvel
570 real,
dimension(SZK_(GV)+1), &
571 intent(out) :: mixlen
575 real,
optional,
intent(in) :: dt_diag
578 optional,
pointer :: Waves
580 optional,
intent(inout) :: G
581 integer,
optional,
intent(in) :: i
582 integer,
optional,
intent(in) :: j
598 real,
dimension(SZK_(GV)+1) :: &
611 real :: Idecay_len_TKE
614 real,
dimension(SZK_(GV)) :: &
648 real,
dimension(SZK_(GV)+1) :: &
724 real :: TKE_left_min, TKE_left_max
725 real :: Kddt_h_max, Kddt_h_min
732 real :: vstar_unit_scale
734 logical :: convectively_stable
735 logical :: sfc_connected
737 logical :: sfc_disconnect
747 real :: MLD_guess, MLD_found
762 logical :: OBL_converged
765 real :: Surface_Scale
767 logical :: debug=.false.
770 real :: dPE_debug, mixing_debug, taux2, tauy2
771 real,
dimension(20) :: TKE_left_itt, PE_chg_itt, Kddt_h_itt, dPEa_dKd_itt, MKE_src_itt
772 real,
dimension(SZK_(GV)) :: mech_TKE_k, conv_PErel_k, nstar_k
773 integer,
dimension(SZK_(GV)) :: num_itts
775 integer :: k, nz, itt, max_itt
779 if (.not.
associated(cs))
call mom_error(fatal,
"energetic_PBL: "//&
780 "Module must be initialized before it is used.")
782 debug = .false. ;
if (
allocated(ecd%dT_expect) .or.
allocated(ecd%dS_expect)) debug = .true.
784 h_neglect = gv%H_subroundoff
787 dt__diag = dt ;
if (
present(dt_diag)) dt__diag = dt_diag
788 i_dtdiag = 1.0 / dt__diag
792 i_dtrho = 0.0 ;
if (dt*gv%Rho0 > 0.0) i_dtrho = (us%Z_to_m**3*us%s_to_T**3) / (dt*gv%Rho0)
793 vstar_unit_scale = us%m_to_Z * us%T_to_s
804 do k=1,nz+1 ; kd(k) = 0.0 ;
enddo
808 dmass = us%m_to_Z * gv%H_to_kg_m2 * h(k)
809 dpres = us%L_to_Z**2 * gv%g_Earth * dmass
810 dt_to_dpe(k) = (dmass * (pres_z(k) + 0.5*dpres)) * dsv_dt(k)
811 ds_to_dpe(k) = (dmass * (pres_z(k) + 0.5*dpres)) * dsv_ds(k)
812 dt_to_dcolht(k) = dmass * dsv_dt(k)
813 ds_to_dcolht(k) = dmass * dsv_ds(k)
815 pres_z(k+1) = pres_z(k) + dpres
819 h_sum = h_neglect ;
do k=1,nz ; h_sum = h_sum + h(k) ;
enddo
820 i_hs = 0.0 ;
if (h_sum > 0.0) i_hs = 1.0 / h_sum
825 hb_hs(k) = h_bot * i_hs
828 mld_output = h(1)*gv%H_to_Z
832 max_mld = 0.0 ;
do k=1,nz ; max_mld = max_mld + h(k)*gv%H_to_Z ;
enddo
837 if (mld_guess <= min_mld) mld_guess = 0.5 * (min_mld + max_mld)
840 obl_converged = .false.
841 do obl_it=1,cs%Max_MLD_Its
843 if (.not. obl_converged)
then
845 if (.not.cs%Use_MLD_iteration) obl_converged = .true.
847 if (debug)
then ; mech_tke_k(:) = 0.0 ; conv_perel_k(:) = 0.0 ;
endif
850 mld_output = h(1)*gv%H_to_Z
851 sfc_connected = .true.
855 call get_langmuir_number(la, g, gv, us, abs(mld_guess), u_star_mean, i, j, &
856 h=h, u_h=u, v_h=v, waves=waves)
857 call find_mstar(cs, us, b_flux, u_star, u_star_mean, mld_guess, absf, &
858 mstar_total, langmuir_number=la, convect_langmuir_number=lamod,&
861 call find_mstar(cs, us, b_flux, u_star, u_star_mean, mld_guess, absf, mstar_total)
865 if ((cs%answers_2018) .and. (cs%mstar_scheme==use_fixed_mstar))
then
866 mech_tke = (dt*mstar_total*gv%Rho0) * u_star**3
868 mech_tke = mstar_total * (dt*gv%Rho0* u_star**3)
871 if (cs%TKE_diagnostics)
then
872 ecd%dTKE_conv = 0.0 ; ecd%dTKE_mixing = 0.0
873 ecd%dTKE_MKE = 0.0 ; ecd%dTKE_mech_decay = 0.0 ; ecd%dTKE_conv_decay = 0.0
875 ecd%dTKE_wind = mech_tke * i_dtdiag
876 if (tke_forcing(1) <= 0.0)
then
877 ecd%dTKE_forcing = max(-mech_tke, tke_forcing(1)) * i_dtdiag
880 ecd%dTKE_forcing = cs%nstar*tke_forcing(1) * i_dtdiag
885 if (tke_forcing(1) <= 0.0)
then
886 mech_tke = mech_tke + tke_forcing(1)
887 if (mech_tke < 0.0) mech_tke = 0.0
890 conv_perel = tke_forcing(1)
895 do k=1,nz+1 ; mixvel(k) = 0.0 ; mixlen(k) = 0.0 ;
enddo
898 if ((.not.cs%Use_MLD_iteration) .or. &
899 (cs%transLay_scale >= 1.0) .or. (cs%transLay_scale < 0.0) )
then
901 mixlen_shape(k) = 1.0
903 elseif (mld_guess <= 0.0)
then
904 if (cs%transLay_scale > 0.0)
then ;
do k=1,nz+1
905 mixlen_shape(k) = cs%transLay_scale
906 enddo ;
else ;
do k=1,nz+1
907 mixlen_shape(k) = 1.0
912 i_mld = 1.0 / mld_guess
914 mixlen_shape(1) = 1.0
916 h_rsum = h_rsum + h(k-1)*gv%H_to_Z
917 if (cs%MixLenExponent==2.0)
then
918 mixlen_shape(k) = cs%transLay_scale + (1.0 - cs%transLay_scale) * &
919 (max(0.0, (mld_guess - h_rsum)*i_mld) )**2
921 mixlen_shape(k) = cs%transLay_scale + (1.0 - cs%transLay_scale) * &
922 (max(0.0, (mld_guess - h_rsum)*i_mld) )**cs%MixLenExponent
927 kd(1) = 0.0 ; kddt_h(1) = 0.0
929 dt_to_dpe_a(1) = dt_to_dpe(1) ; dt_to_dcolht_a(1) = dt_to_dcolht(1)
930 ds_to_dpe_a(1) = ds_to_dpe(1) ; ds_to_dcolht_a(1) = ds_to_dcolht(1)
932 htot = h(1) ; uhtot = u(1)*h(1) ; vhtot = v(1)*h(1)
935 mech_tke_k(1) = mech_tke ; conv_perel_k(1) = conv_perel
936 nstar_k(:) = 0.0 ; nstar_k(1) = cs%nstar ; num_itts(:) = -1
947 idecay_len_tke = (cs%TKE_decay * absf / u_star) * gv%H_to_Z
949 if (idecay_len_tke > 0.0) exp_kh = exp(-h(k-1)*idecay_len_tke)
950 if (cs%TKE_diagnostics) &
951 ecd%dTKE_mech_decay = ecd%dTKE_mech_decay + (exp_kh-1.0) * mech_tke * i_dtdiag
952 mech_tke = mech_tke * exp_kh
956 if (tke_forcing(k) > 0.0)
then
957 conv_perel = conv_perel + tke_forcing(k)
958 if (cs%TKE_diagnostics) &
959 ecd%dTKE_forcing = ecd%dTKE_forcing + cs%nstar*tke_forcing(k) * i_dtdiag
963 mech_tke_k(k) = mech_tke ; conv_perel_k(k) = conv_perel
968 if (cs%nstar * conv_perel > 0.0)
then
972 nstar_fc = cs%nstar * conv_perel / (conv_perel + 0.2 * &
973 sqrt(0.5 * dt * gv%Rho0 * (absf*(htot*gv%H_to_Z))**3 * conv_perel))
976 if (debug) nstar_k(k) = nstar_fc
978 tot_tke = mech_tke + nstar_fc * conv_perel
982 if (tke_forcing(k) < 0.0)
then
983 if (tke_forcing(k) + tot_tke < 0.0)
then
985 if (cs%TKE_diagnostics)
then
986 ecd%dTKE_mixing = ecd%dTKE_mixing + tot_tke * i_dtdiag
987 ecd%dTKE_forcing = ecd%dTKE_forcing - tot_tke * i_dtdiag
989 ecd%dTKE_conv_decay = ecd%dTKE_conv_decay + (cs%nstar-nstar_fc) * conv_perel * i_dtdiag
991 tot_tke = 0.0 ; mech_tke = 0.0 ; conv_perel = 0.0
994 tke_reduc = (tot_tke + tke_forcing(k)) / tot_tke
995 if (cs%TKE_diagnostics)
then
996 ecd%dTKE_mixing = ecd%dTKE_mixing - tke_forcing(k) * i_dtdiag
997 ecd%dTKE_forcing = ecd%dTKE_forcing + tke_forcing(k) * i_dtdiag
998 ecd%dTKE_conv_decay = ecd%dTKE_conv_decay + &
999 (1.0-tke_reduc)*(cs%nstar-nstar_fc) * conv_perel * i_dtdiag
1001 tot_tke = tke_reduc*tot_tke
1002 mech_tke = tke_reduc*mech_tke
1003 conv_perel = tke_reduc*conv_perel
1008 if (cs%orig_PE_calc)
then
1010 dte_t2 = 0.0 ; dse_t2 = 0.0
1012 dte_t2 = kddt_h(k-1) * ((t0(k-2) - t0(k-1)) + dte(k-2))
1013 dse_t2 = kddt_h(k-1) * ((s0(k-2) - s0(k-1)) + dse(k-2))
1016 dt_h = (gv%Z_to_H**2*dt) / max(0.5*(h(k-1)+h(k)), 1e-15*h_sum)
1024 convectively_stable = ( 0.0 <= &
1025 ( (dt_to_dcolht(k) + dt_to_dcolht(k-1) ) * (t0(k-1)-t0(k)) + &
1026 (ds_to_dcolht(k) + ds_to_dcolht(k-1) ) * (s0(k-1)-s0(k)) ) )
1028 if ((mech_tke + conv_perel) <= 0.0 .and. convectively_stable)
then
1030 tot_tke = 0.0 ; mech_tke = 0.0 ; conv_perel = 0.0
1031 kd(k) = 0.0 ; kddt_h(k) = 0.0
1032 sfc_disconnect = .true.
1042 if (cs%orig_PE_calc)
then
1043 dte(k-1) = b1 * ( dte_t2 )
1044 dse(k-1) = b1 * ( dse_t2 )
1048 dt_to_dpe_a(k) = dt_to_dpe(k)
1049 ds_to_dpe_a(k) = ds_to_dpe(k)
1050 dt_to_dcolht_a(k) = dt_to_dcolht(k)
1051 ds_to_dcolht_a(k) = ds_to_dcolht(k)
1054 sfc_disconnect = .false.
1058 if (cs%orig_PE_calc)
then
1060 dt_km1_t2 = (t0(k)-t0(k-1))
1061 ds_km1_t2 = (s0(k)-s0(k-1))
1063 dt_km1_t2 = (t0(k)-t0(k-1)) - &
1064 (kddt_h(k-1) / hp_a) * ((t0(k-2) - t0(k-1)) + dte(k-2))
1065 ds_km1_t2 = (s0(k)-s0(k-1)) - &
1066 (kddt_h(k-1) / hp_a) * ((s0(k-2) - s0(k-1)) + dse(k-2))
1068 dte_term = dte_t2 + hp_a * (t0(k-1)-t0(k))
1069 dse_term = dse_t2 + hp_a * (s0(k-1)-s0(k))
1072 th_a(k-1) = h(k-1) * t0(k-1) ; sh_a(k-1) = h(k-1) * s0(k-1)
1074 th_a(k-1) = h(k-1) * t0(k-1) + kddt_h(k-1) * te(k-2)
1075 sh_a(k-1) = h(k-1) * s0(k-1) + kddt_h(k-1) * se(k-2)
1077 th_b(k) = h(k) * t0(k) ; sh_b(k) = h(k) * s0(k)
1085 if ((cs%MKE_to_TKE_effic > 0.0) .and. (htot*h(k) > 0.0))
then
1088 dmke_max = (us%m_to_Z**3*us%T_to_s**2)*(gv%H_to_kg_m2 * cs%MKE_to_TKE_effic) * 0.5 * &
1089 (h(k) / ((htot + h(k))*htot)) * &
1090 ((uhtot-u(k)*htot)**2 + (vhtot-v(k)*htot)**2)
1093 mke2_hharm = (htot + h(k) + 2.0*h_neglect) / &
1094 ((htot+h_neglect) * (h(k)+h_neglect))
1103 h_tt = htot + h_tt_min
1104 tke_here = mech_tke + cs%wstar_ustar_coef*conv_perel
1105 if (tke_here > 0.0)
then
1106 if (cs%wT_scheme==wt_from_croot_tke)
then
1107 vstar = cs%vstar_scale_fac * vstar_unit_scale * (i_dtrho*tke_here)**c1_3
1108 elseif (cs%wT_scheme==wt_from_rh18)
then
1109 surface_scale = max(0.05, 1.0 - htot/mld_guess)
1110 vstar = cs%vstar_scale_fac * surface_scale * (cs%vstar_surf_fac*u_star + &
1111 vstar_unit_scale * (cs%wstar_ustar_coef*conv_perel*i_dtrho)**c1_3)
1113 hbs_here = gv%H_to_Z * min(hb_hs(k), mixlen_shape(k))
1114 mixlen(k) = max(cs%min_mix_len, ((h_tt*hbs_here)*vstar) / &
1115 ((cs%Ekman_scale_coef * absf) * (h_tt*hbs_here) + vstar))
1118 if (.not.cs%Use_MLD_iteration)
then
1119 kd_guess0 = vstar * cs%vonKar * ((h_tt*hbs_here)*vstar) / &
1120 ((cs%Ekman_scale_coef * absf) * (h_tt*hbs_here) + vstar)
1122 kd_guess0 = vstar * cs%vonKar * mixlen(k)
1125 vstar = 0.0 ; kd_guess0 = 0.0
1128 kddt_h_g0 = kd_guess0 * dt_h
1130 if (cs%orig_PE_calc)
then
1131 call find_pe_chg_orig(kddt_h_g0, h(k), hp_a, dte_term, dse_term, &
1132 dt_km1_t2, ds_km1_t2, dt_to_dpe(k), ds_to_dpe(k), &
1133 dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), &
1134 pres_z(k), dt_to_dcolht(k), ds_to_dcolht(k), &
1135 dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
1136 pe_chg=pe_chg_g0, dpec_dkd=dpea_dkd_g0, dpe_max=pe_chg_max, &
1137 dpec_dkd_0=dpec_dkd_kd0 )
1139 call find_pe_chg(0.0, kddt_h_g0, hp_a, h(k), &
1140 th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
1141 dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), dt_to_dpe(k), ds_to_dpe(k), &
1142 pres_z(k), dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
1143 dt_to_dcolht(k), ds_to_dcolht(k), &
1144 pe_chg=pe_chg_g0, dpec_dkd=dpea_dkd_g0, dpe_max=pe_chg_max, &
1145 dpec_dkd_0=dpec_dkd_kd0 )
1148 mke_src = dmke_max*(1.0 - exp(-kddt_h_g0 * mke2_hharm))
1151 if ((pe_chg_g0 < 0.0) .or. ((vstar == 0.0) .and. (dpec_dkd_kd0 < 0.0)))
then
1153 if (pe_chg_max <= 0.0)
then
1155 tke_here = mech_tke + cs%wstar_ustar_coef*(conv_perel-pe_chg_max)
1156 if (tke_here > 0.0)
then
1157 if (cs%wT_scheme==wt_from_croot_tke)
then
1158 vstar = cs%vstar_scale_fac * vstar_unit_scale * (i_dtrho*tke_here)**c1_3
1159 elseif (cs%wT_scheme==wt_from_rh18)
then
1160 surface_scale = max(0.05, 1. - htot/mld_guess)
1161 vstar = cs%vstar_scale_fac * surface_scale * (cs%vstar_surf_fac*u_star + &
1162 vstar_unit_scale * (cs%wstar_ustar_coef*conv_perel*i_dtrho)**c1_3)
1164 hbs_here = gv%H_to_Z * min(hb_hs(k), mixlen_shape(k))
1165 mixlen(k) = max(cs%min_mix_len, ((h_tt*hbs_here)*vstar) / &
1166 ((cs%Ekman_scale_coef * absf) * (h_tt*hbs_here) + vstar))
1167 if (.not.cs%Use_MLD_iteration)
then
1170 kd(k) = vstar * cs%vonKar * ((h_tt*hbs_here)*vstar) / &
1171 ((cs%Ekman_scale_coef * absf) * (h_tt*hbs_here) + vstar)
1173 kd(k) = vstar * cs%vonKar * mixlen(k)
1176 vstar = 0.0 ; kd(k) = 0.0
1180 if (cs%orig_PE_calc)
then
1181 call find_pe_chg_orig(kd(k)*dt_h, h(k), hp_a, dte_term, dse_term, &
1182 dt_km1_t2, ds_km1_t2, dt_to_dpe(k), ds_to_dpe(k), &
1183 dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), &
1184 pres_z(k), dt_to_dcolht(k), ds_to_dcolht(k), &
1185 dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
1188 call find_pe_chg(0.0, kd(k)*dt_h, hp_a, h(k), &
1189 th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
1190 dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), dt_to_dpe(k), ds_to_dpe(k), &
1191 pres_z(k), dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
1192 dt_to_dcolht(k), ds_to_dcolht(k), &
1196 if (dpe_conv > 0.0)
then
1197 kd(k) = kd_guess0 ; dpe_conv = pe_chg_g0
1199 mke_src = dmke_max*(1.0 - exp(-(kd(k)*dt_h) * mke2_hharm))
1203 kd(k) = kd_guess0 ; dpe_conv = pe_chg_g0
1206 conv_perel = conv_perel - dpe_conv
1207 mech_tke = mech_tke + mke_src
1208 if (cs%TKE_diagnostics)
then
1209 ecd%dTKE_conv = ecd%dTKE_conv - cs%nstar*dpe_conv * i_dtdiag
1210 ecd%dTKE_MKE = ecd%dTKE_MKE + mke_src * i_dtdiag
1212 if (sfc_connected)
then
1213 mld_output = mld_output + gv%H_to_Z * h(k)
1216 kddt_h(k) = kd(k) * dt_h
1217 elseif (tot_tke + (mke_src - pe_chg_g0) >= 0.0)
then
1221 kddt_h(k) = kddt_h_g0
1224 tot_tke = tot_tke + mke_src
1226 if (tot_tke > 0.0) tke_reduc = (tot_tke - pe_chg_g0) / tot_tke
1227 if (cs%TKE_diagnostics)
then
1228 ecd%dTKE_mixing = ecd%dTKE_mixing - pe_chg_g0 * i_dtdiag
1229 ecd%dTKE_MKE = ecd%dTKE_MKE + mke_src * i_dtdiag
1230 ecd%dTKE_conv_decay = ecd%dTKE_conv_decay + &
1231 (1.0-tke_reduc)*(cs%nstar-nstar_fc) * conv_perel * i_dtdiag
1233 tot_tke = tke_reduc*tot_tke
1234 mech_tke = tke_reduc*(mech_tke + mke_src)
1235 conv_perel = tke_reduc*conv_perel
1236 if (sfc_connected)
then
1237 mld_output = mld_output + gv%H_to_Z * h(k)
1240 elseif (tot_tke == 0.0)
then
1242 kd(k) = 0.0 ; kddt_h(k) = 0.0
1243 tot_tke = 0.0 ; conv_perel = 0.0 ; mech_tke = 0.0
1244 sfc_disconnect = .true.
1248 kddt_h_max = kddt_h_g0 ; kddt_h_min = 0.0
1249 tke_left_max = tot_tke + (mke_src - pe_chg_g0)
1250 tke_left_min = tot_tke
1254 kddt_h_guess = tot_tke * kddt_h_max / max( pe_chg_g0 - mke_src, &
1255 kddt_h_max * (dpec_dkd_kd0 - dmke_max * mke2_hharm) )
1262 tke_left_itt(:) = 0.0 ; dpea_dkd_itt(:) = 0.0 ; pe_chg_itt(:) = 0.0
1263 mke_src_itt(:) = 0.0 ; kddt_h_itt(:) = 0.0
1266 if (cs%orig_PE_calc)
then
1267 call find_pe_chg_orig(kddt_h_guess, h(k), hp_a, dte_term, dse_term, &
1268 dt_km1_t2, ds_km1_t2, dt_to_dpe(k), ds_to_dpe(k), &
1269 dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), &
1270 pres_z(k), dt_to_dcolht(k), ds_to_dcolht(k), &
1271 dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
1272 pe_chg=pe_chg, dpec_dkd=dpec_dkd )
1274 call find_pe_chg(0.0, kddt_h_guess, hp_a, h(k), &
1275 th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
1276 dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), dt_to_dpe(k), ds_to_dpe(k), &
1277 pres_z(k), dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
1278 dt_to_dcolht(k), ds_to_dcolht(k), &
1281 mke_src = dmke_max * (1.0 - exp(-mke2_hharm * kddt_h_guess))
1282 dmke_src_dk = dmke_max * mke2_hharm * exp(-mke2_hharm * kddt_h_guess)
1284 tke_left = tot_tke + (mke_src - pe_chg)
1285 if (debug .and. itt<=20)
then
1286 kddt_h_itt(itt) = kddt_h_guess ; mke_src_itt(itt) = mke_src
1287 pe_chg_itt(itt) = pe_chg ; dpea_dkd_itt(itt) = dpec_dkd
1288 tke_left_itt(itt) = tke_left
1292 if (tke_left >= 0.0)
then
1293 kddt_h_min = kddt_h_guess ; tke_left_min = tke_left
1295 kddt_h_max = kddt_h_guess ; tke_left_max = tke_left
1301 if (dpec_dkd - dmke_src_dk <= 0.0)
then
1304 dkddt_h_newt = tke_left / (dpec_dkd - dmke_src_dk)
1305 kddt_h_newt = kddt_h_guess + dkddt_h_newt
1306 if ((kddt_h_newt > kddt_h_max) .or. (kddt_h_newt < kddt_h_min)) &
1311 kddt_h_next = kddt_h_guess + dkddt_h_newt
1312 dkddt_h = dkddt_h_newt
1314 kddt_h_next = (tke_left_max * kddt_h_min - kddt_h_max * tke_left_min) / &
1315 (tke_left_max - tke_left_min)
1316 dkddt_h = kddt_h_next - kddt_h_guess
1319 if ((abs(dkddt_h) < 1e-9*kddt_h_guess) .or. (itt==max_itt))
then
1321 if (debug) num_itts(k) = itt
1324 kddt_h_guess = kddt_h_next
1327 kd(k) = kddt_h_guess / dt_h ; kddt_h(k) = kd(k) * dt_h
1330 if (cs%TKE_diagnostics)
then
1331 ecd%dTKE_mixing = ecd%dTKE_mixing - (tot_tke + mke_src) * i_dtdiag
1332 ecd%dTKE_MKE = ecd%dTKE_MKE + mke_src * i_dtdiag
1333 ecd%dTKE_conv_decay = ecd%dTKE_conv_decay + &
1334 (cs%nstar-nstar_fc) * conv_perel * i_dtdiag
1337 if (sfc_connected) mld_output = mld_output + &
1338 (pe_chg / (pe_chg_g0)) * gv%H_to_Z * h(k)
1340 tot_tke = 0.0 ; mech_tke = 0.0 ; conv_perel = 0.0
1341 sfc_disconnect = .true.
1344 kddt_h(k) = kd(k) * dt_h
1347 b1 = 1.0 / (hp_a + kddt_h(k))
1348 c1(k) = kddt_h(k) * b1
1349 if (cs%orig_PE_calc)
then
1350 dte(k-1) = b1 * ( kddt_h(k)*(t0(k)-t0(k-1)) + dte_t2 )
1351 dse(k-1) = b1 * ( kddt_h(k)*(s0(k)-s0(k-1)) + dse_t2 )
1354 hp_a = h(k) + (hp_a * b1) * kddt_h(k)
1355 dt_to_dpe_a(k) = dt_to_dpe(k) + c1(k)*dt_to_dpe_a(k-1)
1356 ds_to_dpe_a(k) = ds_to_dpe(k) + c1(k)*ds_to_dpe_a(k-1)
1357 dt_to_dcolht_a(k) = dt_to_dcolht(k) + c1(k)*dt_to_dcolht_a(k-1)
1358 ds_to_dcolht_a(k) = ds_to_dcolht(k) + c1(k)*ds_to_dcolht_a(k-1)
1363 if (sfc_disconnect)
then
1368 sfc_connected = .false.
1370 uhtot = uhtot + u(k)*h(k)
1371 vhtot = vhtot + v(k)*h(k)
1377 te(1) = b1*(h(1)*t0(1))
1378 se(1) = b1*(h(1)*s0(1))
1380 te(k-1) = b1 * (h(k-1) * t0(k-1) + kddt_h(k-1) * te(k-2))
1381 se(k-1) = b1 * (h(k-1) * s0(k-1) + kddt_h(k-1) * se(k-2))
1390 te(nz) = b1 * (h(nz) * t0(nz) + kddt_h(nz) * te(nz-1))
1391 se(nz) = b1 * (h(nz) * s0(nz) + kddt_h(nz) * se(nz-1))
1392 ecd%dT_expect(nz) = te(nz) - t0(nz) ; ecd%dS_expect(nz) = se(nz) - s0(nz)
1394 te(k) = te(k) + c1(k+1)*te(k+1)
1395 se(k) = se(k) + c1(k+1)*se(k+1)
1396 ecd%dT_expect(k) = te(k) - t0(k) ; ecd%dS_expect(k) = se(k) - s0(k)
1401 dpe_debug = dpe_debug + (dt_to_dpe(k) * (te(k) - t0(k)) + &
1402 ds_to_dpe(k) * (se(k) - s0(k)))
1404 mixing_debug = dpe_debug * i_dtdiag
1415 mld_found = mld_output
1416 if (mld_found - cs%MLD_tol > mld_guess)
then
1418 elseif (abs(mld_guess - mld_found) < cs%MLD_tol)
then
1419 obl_converged = .true.
1426 mld_guess = 0.5*(min_mld + max_mld)
1430 ecd%LA = la ; ecd%LAmod = lamod ; ecd%mstar = mstar_total ; ecd%mstar_LT = mstar_lt
1432 ecd%LA = 0.0 ; ecd%LAmod = 0.0 ; ecd%mstar = mstar_total ; ecd%mstar_LT = 0.0
1437 end subroutine epbl_column
1441 subroutine find_pe_chg(Kddt_h0, dKddt_h, hp_a, hp_b, Th_a, Sh_a, Th_b, Sh_b, &
1442 dT_to_dPE_a, dS_to_dPE_a, dT_to_dPE_b, dS_to_dPE_b, &
1443 pres_Z, dT_to_dColHt_a, dS_to_dColHt_a, dT_to_dColHt_b, dS_to_dColHt_b, &
1444 PE_chg, dPEc_dKd, dPE_max, dPEc_dKd_0, ColHt_cor)
1445 real,
intent(in) :: Kddt_h0
1448 real,
intent(in) :: dKddt_h
1451 real,
intent(in) :: hp_a
1455 real,
intent(in) :: hp_b
1459 real,
intent(in) :: Th_a
1462 real,
intent(in) :: Sh_a
1465 real,
intent(in) :: Th_b
1468 real,
intent(in) :: Sh_b
1471 real,
intent(in) :: dT_to_dPE_a
1475 real,
intent(in) :: dS_to_dPE_a
1479 real,
intent(in) :: dT_to_dPE_b
1483 real,
intent(in) :: dS_to_dPE_b
1487 real,
intent(in) :: pres_Z
1490 real,
intent(in) :: dT_to_dColHt_a
1494 real,
intent(in) :: dS_to_dColHt_a
1498 real,
intent(in) :: dT_to_dColHt_b
1502 real,
intent(in) :: dS_to_dColHt_b
1507 real,
optional,
intent(out) :: PE_chg
1509 real,
optional,
intent(out) :: dPEc_dKd
1511 real,
optional,
intent(out) :: dPE_max
1514 real,
optional,
intent(out) :: dPEc_dKd_0
1516 real,
optional,
intent(out) :: ColHt_cor
1540 bdt1 = hp_a * hp_b + kddt_h0 * hps
1541 dt_c = hp_a * th_b - hp_b * th_a
1542 ds_c = hp_a * sh_b - hp_b * sh_a
1543 pec_core = hp_b * (dt_to_dpe_a * dt_c + ds_to_dpe_a * ds_c) - &
1544 hp_a * (dt_to_dpe_b * dt_c + ds_to_dpe_b * ds_c)
1545 colht_core = hp_b * (dt_to_dcolht_a * dt_c + ds_to_dcolht_a * ds_c) - &
1546 hp_a * (dt_to_dcolht_b * dt_c + ds_to_dcolht_b * ds_c)
1548 if (
present(pe_chg))
then
1551 y1_3 = dkddt_h / (bdt1 * (bdt1 + dkddt_h * hps))
1552 pe_chg = pec_core * y1_3
1553 colht_chg = colht_core * y1_3
1554 if (colht_chg < 0.0) pe_chg = pe_chg - pres_z * colht_chg
1555 if (
present(colht_cor)) colht_cor = -pres_z * min(colht_chg, 0.0)
1556 elseif (
present(colht_cor))
then
1557 y1_3 = dkddt_h / (bdt1 * (bdt1 + dkddt_h * hps))
1558 colht_cor = -pres_z * min(colht_core * y1_3, 0.0)
1561 if (
present(dpec_dkd))
then
1563 y1_4 = 1.0 / (bdt1 + dkddt_h * hps)**2
1564 dpec_dkd = pec_core * y1_4
1565 colht_chg = colht_core * y1_4
1566 if (colht_chg < 0.0) dpec_dkd = dpec_dkd - pres_z * colht_chg
1569 if (
present(dpe_max))
then
1571 y1_3 = 1.0 / (bdt1 * hps)
1572 dpe_max = pec_core * y1_3
1573 colht_chg = colht_core * y1_3
1574 if (colht_chg < 0.0) dpe_max = dpe_max - pres_z * colht_chg
1577 if (
present(dpec_dkd_0))
then
1579 y1_4 = 1.0 / bdt1**2
1580 dpec_dkd_0 = pec_core * y1_4
1581 colht_chg = colht_core * y1_4
1582 if (colht_chg < 0.0) dpec_dkd_0 = dpec_dkd_0 - pres_z * colht_chg
1585 end subroutine find_pe_chg
1590 subroutine find_pe_chg_orig(Kddt_h, h_k, b_den_1, dTe_term, dSe_term, &
1591 dT_km1_t2, dS_km1_t2, dT_to_dPE_k, dS_to_dPE_k, &
1592 dT_to_dPEa, dS_to_dPEa, pres_Z, dT_to_dColHt_k, &
1593 dS_to_dColHt_k, dT_to_dColHta, dS_to_dColHta, &
1594 PE_chg, dPEc_dKd, dPE_max, dPEc_dKd_0)
1595 real,
intent(in) :: Kddt_h
1598 real,
intent(in) :: h_k
1599 real,
intent(in) :: b_den_1
1603 real,
intent(in) :: dTe_term
1605 real,
intent(in) :: dSe_term
1607 real,
intent(in) :: dT_km1_t2
1609 real,
intent(in) :: dS_km1_t2
1611 real,
intent(in) :: pres_Z
1614 real,
intent(in) :: dT_to_dPE_k
1618 real,
intent(in) :: dS_to_dPE_k
1622 real,
intent(in) :: dT_to_dPEa
1626 real,
intent(in) :: dS_to_dPEa
1630 real,
intent(in) :: dT_to_dColHt_k
1634 real,
intent(in) :: dS_to_dColHt_k
1638 real,
intent(in) :: dT_to_dColHta
1642 real,
intent(in) :: dS_to_dColHta
1647 real,
optional,
intent(out) :: PE_chg
1649 real,
optional,
intent(out) :: dPEc_dKd
1651 real,
optional,
intent(out) :: dPE_max
1654 real,
optional,
intent(out) :: dPEc_dKd_0
1671 real :: dT_k, dT_km1
1672 real :: dS_k, dS_km1
1675 real :: ddT_k_dKd, ddT_km1_dKd
1676 real :: ddS_k_dKd, ddS_km1_dKd
1678 b1 = 1.0 / (b_den_1 + kddt_h)
1686 i_kr_denom = 1.0 / (h_k*b_den_1 + (b_den_1 + h_k)*kddt_h)
1688 dt_k = (kddt_h*i_kr_denom) * dte_term
1689 ds_k = (kddt_h*i_kr_denom) * dse_term
1691 if (
present(pe_chg))
then
1694 dt_km1 = b1kd * ( dt_k + dt_km1_t2 )
1695 ds_km1 = b1kd * ( ds_k + ds_km1_t2 )
1696 pe_chg = (dt_to_dpe_k * dt_k + dt_to_dpea * dt_km1) + &
1697 (ds_to_dpe_k * ds_k + ds_to_dpea * ds_km1)
1698 colht_chg = (dt_to_dcolht_k * dt_k + dt_to_dcolhta * dt_km1) + &
1699 (ds_to_dcolht_k * ds_k + ds_to_dcolhta * ds_km1)
1700 if (colht_chg < 0.0) pe_chg = pe_chg - pres_z * colht_chg
1703 if (
present(dpec_dkd))
then
1705 dkr_dkd = (h_k*b_den_1) * i_kr_denom**2
1707 ddt_k_dkd = dkr_dkd * dte_term
1708 dds_k_dkd = dkr_dkd * dse_term
1709 ddt_km1_dkd = (b1**2 * b_den_1) * ( dt_k + dt_km1_t2 ) + b1kd * ddt_k_dkd
1710 dds_km1_dkd = (b1**2 * b_den_1) * ( ds_k + ds_km1_t2 ) + b1kd * dds_k_dkd
1713 dpec_dkd = (dt_to_dpe_k * ddt_k_dkd + dt_to_dpea * ddt_km1_dkd) + &
1714 (ds_to_dpe_k * dds_k_dkd + ds_to_dpea * dds_km1_dkd)
1715 dcolht_dkd = (dt_to_dcolht_k * ddt_k_dkd + dt_to_dcolhta * ddt_km1_dkd) + &
1716 (ds_to_dcolht_k * dds_k_dkd + ds_to_dcolhta * dds_km1_dkd)
1717 if (dcolht_dkd < 0.0) dpec_dkd = dpec_dkd - pres_z * dcolht_dkd
1720 if (
present(dpe_max))
then
1722 dpe_max = (dt_to_dpea * dt_km1_t2 + ds_to_dpea * ds_km1_t2) + &
1723 ((dt_to_dpe_k + dt_to_dpea) * dte_term + &
1724 (ds_to_dpe_k + ds_to_dpea) * dse_term) / (b_den_1 + h_k)
1725 dcolht_max = (dt_to_dcolhta * dt_km1_t2 + ds_to_dcolhta * ds_km1_t2) + &
1726 ((dt_to_dcolht_k + dt_to_dcolhta) * dte_term + &
1727 (ds_to_dcolht_k + ds_to_dcolhta) * dse_term) / (b_den_1 + h_k)
1728 if (dcolht_max < 0.0) dpe_max = dpe_max - pres_z*dcolht_max
1731 if (
present(dpec_dkd_0))
then
1733 dpec_dkd_0 = (dt_to_dpea * dt_km1_t2 + ds_to_dpea * ds_km1_t2) / (b_den_1) + &
1734 (dt_to_dpe_k * dte_term + ds_to_dpe_k * dse_term) / (h_k*b_den_1)
1735 dcolht_dkd = (dt_to_dcolhta * dt_km1_t2 + ds_to_dcolhta * ds_km1_t2) / (b_den_1) + &
1736 (dt_to_dcolht_k * dte_term + ds_to_dcolht_k * dse_term) / (h_k*b_den_1)
1737 if (dcolht_dkd < 0.0) dpec_dkd_0 = dpec_dkd_0 - pres_z*dcolht_dkd
1740 end subroutine find_pe_chg_orig
1743 subroutine find_mstar(CS, US, Buoyancy_Flux, UStar, UStar_Mean,&
1744 BLD, Abs_Coriolis, MStar, Langmuir_Number,&
1745 MStar_LT, Convect_Langmuir_Number)
1748 real,
intent(in) :: UStar
1749 real,
intent(in) :: UStar_Mean
1750 real,
intent(in) :: Abs_Coriolis
1751 real,
intent(in) :: Buoyancy_Flux
1752 real,
intent(in) :: BLD
1753 real,
intent(out) :: Mstar
1754 real,
optional,
intent(in) :: Langmuir_Number
1755 real,
optional,
intent(out) :: MStar_LT
1756 real,
optional,
intent(out) :: Convect_Langmuir_number
1760 real :: MSCR_term1, MSCR_term2
1761 real :: MStar_Conv_Red
1762 real :: MStar_S, MStar_N
1768 if (cs%mstar_scheme == use_fixed_mstar)
then
1769 mstar = cs%Fixed_MStar
1771 elseif (cs%mstar_scheme == mstar_from_ekman)
then
1773 if (cs%answers_2018)
then
1775 mstar_s = cs%MStar_coef*sqrt(max(0.0,buoyancy_flux) / ustar**2 / &
1776 (abs_coriolis + 1.e-10*us%T_to_s) )
1778 mstar_n = cs%C_Ek * log( max( 1., ustar / (abs_coriolis + 1.e-10*us%T_to_s) / bld ) )
1781 mstar_s = cs%MSTAR_COEF*sqrt(max(0.0, buoyancy_flux) / (ustar**2 * max(abs_coriolis, 1.e-20*us%T_to_s)))
1784 if (ustar > abs_coriolis * bld) mstar_n = cs%C_EK * log(ustar / (abs_coriolis * bld))
1788 mstar = max(mstar_s, min(1.25, mstar_n))
1789 if (cs%MStar_Cap > 0.0) mstar = min( cs%MStar_Cap,mstar )
1790 elseif ( cs%mstar_scheme == mstar_from_rh18 )
then
1791 if (cs%answers_2018)
then
1792 mstar_n = cs%RH18_MStar_cn1 * ( 1.0 - 1.0 / ( 1. + cs%RH18_MStar_cn2 * &
1793 exp( cs%RH18_mstar_CN3 * bld * abs_coriolis / ustar) ) )
1795 msn_term = cs%RH18_MStar_cn2 * exp( cs%RH18_mstar_CN3 * bld * abs_coriolis / ustar)
1796 mstar_n = (cs%RH18_MStar_cn1 * msn_term) / ( 1. + msn_term)
1798 mstar_s = cs%RH18_MStar_CS1 * ( max(0.0, buoyancy_flux)**2 * bld / &
1799 ( ustar**5 * max(abs_coriolis,1.e-20*us%T_to_s) ) )**cs%RH18_mstar_cs2
1800 mstar = mstar_n + mstar_s
1804 if (cs%answers_2018)
then
1805 mstar_conv_red = 1. - cs%MStar_Convect_coef * (-min(0.0,buoyancy_flux) + 1.e-10*us%T_to_s**3*us%m_to_Z**2) / &
1806 ( (-min(0.0,buoyancy_flux) + 1.e-10*us%T_to_s**3*us%m_to_Z**2) + &
1807 2.0 *mstar * ustar**3 / bld )
1809 mscr_term1 = -bld * min(0.0, buoyancy_flux)
1810 mscr_term2 = 2.0*mstar * ustar**3
1811 mstar_conv_red = ((1.-cs%mstar_convect_coef) * mscr_term1 + mscr_term2) / (mscr_term1 + mscr_term2)
1815 mstar = mstar * mstar_conv_red
1817 if (
present(langmuir_number))
then
1819 call mstar_langmuir(cs, us, abs_coriolis, buoyancy_flux, ustar, bld, langmuir_number, mstar, &
1820 mstar_lt, convect_langmuir_number)
1823 end subroutine find_mstar
1826 subroutine mstar_langmuir(CS, US, abs_Coriolis, buoyancy_flux, ustar, BLD, Langmuir_Number, &
1827 mstar, mstar_LT, Convect_Langmuir_Number)
1830 real,
intent(in) :: Abs_Coriolis
1831 real,
intent(in) :: Buoyancy_Flux
1832 real,
intent(in) :: UStar
1833 real,
intent(in) :: BLD
1834 real,
intent(inout) :: Mstar
1835 real,
intent(in) :: Langmuir_Number
1836 real,
intent(out) :: MStar_LT
1837 real,
intent(out) :: Convect_Langmuir_number
1840 real,
parameter :: Max_ratio = 1.0e16
1841 real :: enhance_mstar
1842 real :: mstar_LT_add
1848 real :: Ekman_Obukhov
1850 real :: MLD_Obukhov_stab
1851 real :: Ekman_Obukhov_stab
1852 real :: MLD_Obukhov_un
1853 real :: Ekman_Obukhov_un
1856 enhance_mstar = 1.0 ; mstar_lt_add = 0.0
1858 if (cs%LT_Enhance_Form /= no_langmuir)
then
1860 if (cs%answers_2018)
then
1861 il_ekman = abs_coriolis / ustar
1862 il_obukhov = buoyancy_flux*cs%vonkar / ustar**3
1863 ekman_obukhov_stab = abs(max(0., il_obukhov / (il_ekman + 1.e-10*us%Z_to_m)))
1864 ekman_obukhov_un = abs(min(0., il_obukhov / (il_ekman + 1.e-10*us%Z_to_m)))
1865 mld_obukhov_stab = abs(max(0., bld*il_obukhov))
1866 mld_obukhov_un = abs(min(0., bld*il_obukhov))
1867 mld_ekman = abs( bld*il_ekman )
1869 ekman_obukhov = max_ratio ; mld_obukhov = max_ratio ; mld_ekman = max_ratio
1870 i_f = 0.0 ;
if (abs(abs_coriolis) > 0.0) i_f = 1.0 / abs_coriolis
1871 i_ustar = 0.0 ;
if (abs(ustar) > 0.0) i_ustar = 1.0 / ustar
1872 if (abs(buoyancy_flux*cs%vonkar) < max_ratio*(abs_coriolis * ustar**2)) &
1873 ekman_obukhov = abs(buoyancy_flux*cs%vonkar) * (i_f * i_ustar**2)
1874 if (abs(bld*buoyancy_flux*cs%vonkar) < max_ratio*ustar**3) &
1875 mld_obukhov = abs(bld*buoyancy_flux*cs%vonkar) * i_ustar**3
1876 if (bld*abs_coriolis < max_ratio*ustar) &
1877 mld_ekman = bld*abs_coriolis * i_ustar
1879 if (buoyancy_flux > 0.0)
then
1880 ekman_obukhov_stab = ekman_obukhov ; ekman_obukhov_un = 0.0
1881 mld_obukhov_stab = mld_obukhov ; mld_obukhov_un = 0.0
1883 ekman_obukhov_un = ekman_obukhov ; ekman_obukhov_stab = 0.0
1884 mld_obukhov_un = mld_obukhov ; mld_obukhov_stab = 0.0
1891 convect_langmuir_number = langmuir_number * &
1892 ( (1.0 + max(-0.5, cs%LaC_MLDoEK * mld_ekman)) + &
1893 ((cs%LaC_EKoOB_stab * ekman_obukhov_stab + cs%LaC_EKoOB_un * ekman_obukhov_un) + &
1894 (cs%LaC_MLDoOB_stab * mld_obukhov_stab + cs%LaC_MLDoOB_un * mld_obukhov_un)) )
1896 if (cs%LT_Enhance_Form == langmuir_rescale)
then
1898 enhance_mstar = min(cs%Max_Enhance_M, &
1899 (1. + cs%LT_ENHANCE_COEF * convect_langmuir_number**cs%LT_ENHANCE_EXP) )
1900 elseif (cs%LT_ENHANCE_Form == langmuir_add)
then
1902 mstar_lt_add = cs%LT_ENHANCE_COEF * convect_langmuir_number**cs%LT_ENHANCE_EXP
1906 mstar_lt = (enhance_mstar - 1.0)*mstar + mstar_lt_add
1907 mstar = mstar*enhance_mstar + mstar_lt_add
1909 end subroutine mstar_langmuir
1913 subroutine energetic_pbl_get_mld(CS, MLD, G, US, m_to_MLD_units)
1917 real,
dimension(SZI_(G),SZJ_(G)),
intent(out) :: mld
1918 real,
optional,
intent(in) :: m_to_mld_units
1924 scale = us%Z_to_m ;
if (
present(m_to_mld_units)) scale = scale * m_to_mld_units
1926 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
1927 mld(i,j) = scale*cs%ML_Depth(i,j)
1930 end subroutine energetic_pbl_get_mld
1934 subroutine energetic_pbl_init(Time, G, GV, US, param_file, diag, CS)
1935 type(time_type),
target,
intent(in) :: time
1940 type(
diag_ctrl),
target,
intent(inout) :: diag
1945 # include "version_variable.h"
1946 character(len=40) :: mdl =
"MOM_energetic_PBL"
1947 character(len=20) :: tmpstr
1948 real :: omega_frac_dflt
1949 real :: z3_t3_to_m3_s3
1950 integer :: isd, ied, jsd, jed
1951 integer :: mstar_mode, lt_enhance, wt_mode
1952 logical :: default_2018_answers
1953 logical :: use_temperature, use_omega
1954 logical :: use_la_windsea
1955 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1957 if (
associated(cs))
then
1958 call mom_error(warning,
"mixedlayer_init called with an associated"//&
1959 "associated control structure.")
1961 else ;
allocate(cs) ;
endif
1971 call get_param(param_file, mdl,
"OMEGA", cs%omega, &
1972 "The rotation rate of the earth.", units=
"s-1", &
1973 default=7.2921e-5, scale=us%T_to_S)
1974 call get_param(param_file, mdl,
"ML_USE_OMEGA", use_omega, &
1975 "If true, use the absolute rotation rate instead of the "//&
1976 "vertical component of rotation when setting the decay "//&
1977 "scale for turbulence.", default=.false., do_not_log=.true.)
1978 omega_frac_dflt = 0.0
1980 call mom_error(warning,
"ML_USE_OMEGA is depricated; use ML_OMEGA_FRAC=1.0 instead.")
1981 omega_frac_dflt = 1.0
1983 call get_param(param_file, mdl,
"ML_OMEGA_FRAC", cs%omega_frac, &
1984 "When setting the decay scale for turbulence, use this "//&
1985 "fraction of the absolute rotation rate blended with the "//&
1986 "local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).", &
1987 units=
"nondim", default=omega_frac_dflt)
1988 call get_param(param_file, mdl,
"EKMAN_SCALE_COEF", cs%Ekman_scale_coef, &
1989 "A nondimensional scaling factor controlling the inhibition "//&
1990 "of the diffusive length scale by rotation. Making this larger "//&
1991 "decreases the PBL diffusivity.", units=
"nondim", default=1.0)
1992 call get_param(param_file, mdl,
"DEFAULT_2018_ANSWERS", default_2018_answers, &
1993 "This sets the default value for the various _2018_ANSWERS parameters.", &
1995 call get_param(param_file, mdl,
"EPBL_2018_ANSWERS", cs%answers_2018, &
1996 "If true, use the order of arithmetic and expressions that recover the "//&
1997 "answers from the end of 2018. Otherwise, use updated and more robust "//&
1998 "forms of the same expressions.", default=default_2018_answers)
2001 call get_param(param_file, mdl,
"EPBL_ORIGINAL_PE_CALC", cs%orig_PE_calc, &
2002 "If true, the ePBL code uses the original form of the "//&
2003 "potential energy change code. Otherwise, the newer "//&
2004 "version that can work with successive increments to the "//&
2005 "diffusivity in upward or downward passes is used.", default=.true.)
2007 call get_param(param_file, mdl,
"MKE_TO_TKE_EFFIC", cs%MKE_to_TKE_effic, &
2008 "The efficiency with which mean kinetic energy released "//&
2009 "by mechanically forced entrainment of the mixed layer "//&
2010 "is converted to turbulent kinetic energy.", units=
"nondim", &
2012 call get_param(param_file, mdl,
"TKE_DECAY", cs%TKE_decay, &
2013 "TKE_DECAY relates the vertical rate of decay of the "//&
2014 "TKE available for mechanical entrainment to the natural "//&
2015 "Ekman depth.", units=
"nondim", default=2.5)
2019 call get_param(param_file, mdl,
"EPBL_MSTAR_SCHEME", tmpstr, &
2020 "EPBL_MSTAR_SCHEME selects the method for setting mstar. Valid values are: \n"//&
2021 "\t CONSTANT - Use a fixed mstar given by MSTAR \n"//&
2022 "\t OM4 - Use L_Ekman/L_Obukhov in the sabilizing limit, as in OM4 \n"//&
2023 "\t REICHL_H18 - Use the scheme documented in Reichl & Hallberg, 2018.", &
2024 default=constant_string, do_not_log=.true.)
2025 call get_param(param_file, mdl,
"MSTAR_MODE", mstar_mode, default=-1)
2026 if (mstar_mode == 0)
then
2027 tmpstr = constant_string
2028 call mom_error(warning,
"Use EPBL_MSTAR_SCHEME = CONSTANT instead of the archaic MSTAR_MODE = 0.")
2029 elseif (mstar_mode == 1)
then
2030 call mom_error(fatal,
"You are using a legacy mstar mode in ePBL that has been phased out. "//&
2031 "If you need to use this setting please report this error. Also use "//&
2032 "EPBL_MSTAR_SCHEME to specify the scheme for mstar.")
2033 elseif (mstar_mode == 2)
then
2035 call mom_error(warning,
"Use EPBL_MSTAR_SCHEME = OM4 instead of the archaic MSTAR_MODE = 2.")
2036 elseif (mstar_mode == 3)
then
2037 tmpstr = rh18_string
2038 call mom_error(warning,
"Use EPBL_MSTAR_SCHEME = REICHL_H18 instead of the archaic MSTAR_MODE = 3.")
2039 elseif (mstar_mode > 3)
then
2040 call mom_error(fatal,
"An unrecognized value of the obsolete parameter MSTAR_MODE was specified.")
2042 call log_param(param_file, mdl,
"EPBL_MSTAR_SCHEME", tmpstr, &
2043 "EPBL_MSTAR_SCHEME selects the method for setting mstar. Valid values are: \n"//&
2044 "\t CONSTANT - Use a fixed mstar given by MSTAR \n"//&
2045 "\t OM4 - Use L_Ekman/L_Obukhov in the sabilizing limit, as in OM4 \n"//&
2046 "\t REICHL_H18 - Use the scheme documented in Reichl & Hallberg, 2018.", &
2047 default=constant_string)
2048 tmpstr = uppercase(tmpstr)
2049 select case (tmpstr)
2050 case (constant_string)
2051 cs%mstar_Scheme = use_fixed_mstar
2053 cs%mstar_Scheme = mstar_from_ekman
2055 cs%mstar_Scheme = mstar_from_rh18
2057 call mom_mesg(
'energetic_PBL_init: EPBL_MSTAR_SCHEME ="'//trim(tmpstr)//
'"', 0)
2058 call mom_error(fatal,
"energetic_PBL_init: Unrecognized setting "// &
2059 "EPBL_MSTAR_SCHEME = "//trim(tmpstr)//
" found in input file.")
2062 call get_param(param_file, mdl,
"MSTAR", cs%fixed_mstar, &
2063 "The ratio of the friction velocity cubed to the TKE input to the "//&
2064 "mixed layer. This option is used if EPBL_MSTAR_SCHEME = CONSTANT.", &
2065 units=
"nondim", default=1.2, do_not_log=(cs%mstar_scheme/=use_fixed_mstar))
2066 call get_param(param_file, mdl,
"MSTAR_CAP", cs%mstar_cap, &
2067 "If this value is positive, it sets the maximum value of mstar "//&
2068 "allowed in ePBL. (This is not used if EPBL_MSTAR_SCHEME = CONSTANT).", &
2069 units=
"nondim", default=-1.0, do_not_log=(cs%mstar_scheme==use_fixed_mstar))
2071 call get_param(param_file, mdl,
"MSTAR2_COEF1", cs%MSTAR_COEF, &
2072 "Coefficient in computing mstar when rotation and stabilizing "//&
2073 "effects are both important (used if EPBL_MSTAR_SCHEME = OM4).", &
2074 units=
"nondim", default=0.3, do_not_log=(cs%mstar_scheme/=mstar_from_ekman))
2075 call get_param(param_file, mdl,
"MSTAR2_COEF2", cs%C_EK, &
2076 "Coefficient in computing mstar when only rotation limits "// &
2077 "the total mixing (used if EPBL_MSTAR_SCHEME = OM4)", &
2078 units=
"nondim", default=0.085, do_not_log=(cs%mstar_scheme/=mstar_from_ekman))
2080 call get_param(param_file, mdl,
"RH18_MSTAR_CN1", cs%RH18_mstar_cn1,&
2081 "MSTAR_N coefficient 1 (outter-most coefficient for fit). "//&
2082 "The value of 0.275 is given in RH18. Increasing this "//&
2083 "coefficient increases MSTAR for all values of Hf/ust, but more "//&
2084 "effectively at low values (weakly developed OSBLs).", &
2085 units=
"nondim", default=0.275, do_not_log=(cs%mstar_scheme/=mstar_from_rh18))
2086 call get_param(param_file, mdl,
"RH18_MSTAR_CN2", cs%RH18_mstar_cn2,&
2087 "MSTAR_N coefficient 2 (coefficient outside of exponential decay). "//&
2088 "The value of 8.0 is given in RH18. Increasing this coefficient "//&
2089 "increases MSTAR for all values of HF/ust, with a much more even "//&
2090 "effect across a wide range of Hf/ust than CN1.", &
2091 units=
"nondim", default=8.0, do_not_log=(cs%mstar_scheme/=mstar_from_rh18))
2092 call get_param(param_file, mdl,
"RH18_MSTAR_CN3", cs%RH18_mstar_CN3,&
2093 "MSTAR_N coefficient 3 (exponential decay coefficient). "//&
2094 "The value of -5.0 is given in RH18. Increasing this increases how "//&
2095 "quickly the value of MSTAR decreases as Hf/ust increases.", &
2096 units=
"nondim", default=-5.0, do_not_log=(cs%mstar_scheme/=mstar_from_rh18))
2097 call get_param(param_file, mdl,
"RH18_MSTAR_CS1", cs%RH18_mstar_cs1,&
2098 "MSTAR_S coefficient for RH18 in stabilizing limit. "//&
2099 "The value of 0.2 is given in RH18 and increasing it increases "//&
2100 "MSTAR in the presence of a stabilizing surface buoyancy flux.", &
2101 units=
"nondim", default=0.2, do_not_log=(cs%mstar_scheme/=mstar_from_rh18))
2102 call get_param(param_file, mdl,
"RH18_MSTAR_CS2", cs%RH18_mstar_cs2,&
2103 "MSTAR_S exponent for RH18 in stabilizing limit. "//&
2104 "The value of 0.4 is given in RH18 and increasing it increases MSTAR "//&
2105 "exponentially in the presence of a stabilizing surface buoyancy flux.", &
2106 units=
"nondim", default=0.4, do_not_log=(cs%mstar_scheme/=mstar_from_rh18))
2110 call get_param(param_file, mdl,
"NSTAR", cs%nstar, &
2111 "The portion of the buoyant potential energy imparted by "//&
2112 "surface fluxes that is available to drive entrainment "//&
2113 "at the base of mixed layer when that energy is positive.", &
2114 units=
"nondim", default=0.2)
2115 call get_param(param_file, mdl,
"MSTAR_CONV_ADJ", cs%mstar_convect_coef, &
2116 "Coefficient used for reducing mstar during convection "//&
2117 "due to reduction of stable density gradient.", &
2118 units=
"nondim", default=0.0)
2122 call get_param(param_file, mdl,
"USE_MLD_ITERATION", cs%Use_MLD_iteration, &
2123 "A logical that specifies whether or not to use the "//&
2124 "distance to the bottom of the actively turbulent boundary "//&
2125 "layer to help set the EPBL length scale.", default=.false.)
2126 call get_param(param_file, mdl,
"EPBL_TRANSITION_SCALE", cs%transLay_scale, &
2127 "A scale for the mixing length in the transition layer "//&
2128 "at the edge of the boundary layer as a fraction of the "//&
2129 "boundary layer thickness.", units=
"nondim", default=0.1)
2130 if ( cs%Use_MLD_iteration .and. abs(cs%transLay_scale-0.5) >= 0.5)
then
2131 call mom_error(fatal,
"If flag USE_MLD_ITERATION is true, then "//&
2132 "EPBL_TRANSITION should be greater than 0 and less than 1.")
2135 call get_param(param_file, mdl,
"MLD_ITERATION_GUESS", cs%MLD_ITERATION_GUESS, &
2136 "A logical that specifies whether or not to use the "//&
2137 "previous timestep MLD as a first guess in the MLD iteration. "//&
2138 "The default is false to facilitate reproducibility.", default=.false.)
2139 call get_param(param_file, mdl,
"EPBL_MLD_TOLERANCE", cs%MLD_tol, &
2140 "The tolerance for the iteratively determined mixed "//&
2141 "layer depth. This is only used with USE_MLD_ITERATION.", &
2142 units=
"meter", default=1.0, scale=us%m_to_Z)
2143 call get_param(param_file, mdl,
"EPBL_MLD_MAX_ITS", cs%max_MLD_its, &
2144 "The maximum number of iterations that can be used to find a self-consistent "//&
2145 "mixed layer depth. For now, due to the use of bisection, the maximum number "//&
2146 "iteractions needed is set by Depth/2^MAX_ITS < EPBL_MLD_TOLERANCE.", &
2147 default=20, do_not_log=.not.cs%Use_MLD_iteration)
2148 if (.not.cs%Use_MLD_iteration) cs%Max_MLD_Its = 1
2149 call get_param(param_file, mdl,
"EPBL_MIN_MIX_LEN", cs%min_mix_len, &
2150 "The minimum mixing length scale that will be used "//&
2151 "by ePBL. The default (0) does not set a minimum.", &
2152 units=
"meter", default=0.0, scale=us%m_to_Z)
2154 call get_param(param_file, mdl,
"MIX_LEN_EXPONENT", cs%MixLenExponent, &
2155 "The exponent applied to the ratio of the distance to the MLD "//&
2156 "and the MLD depth which determines the shape of the mixing length. "//&
2157 "This is only used if USE_MLD_ITERATION is True.", &
2158 units=
"nondim", default=2.0)
2161 call get_param(param_file, mdl,
"EPBL_VEL_SCALE_SCHEME", tmpstr, &
2162 "Selects the method for translating TKE into turbulent velocities. "//&
2163 "Valid values are: \n"//&
2164 "\t CUBE_ROOT_TKE - A constant times the cube root of remaining TKE. \n"//&
2165 "\t REICHL_H18 - Use the scheme based on a combination of w* and v* as \n"//&
2166 "\t documented in Reichl & Hallberg, 2018.", &
2167 default=root_tke_string, do_not_log=.true.)
2168 call get_param(param_file, mdl,
"EPBL_VEL_SCALE_MODE", wt_mode, default=-1)
2169 if (wt_mode == 0)
then
2170 tmpstr = root_tke_string
2171 call mom_error(warning,
"Use EPBL_VEL_SCALE_SCHEME = CUBE_ROOT_TKE instead of the archaic EPBL_VEL_SCALE_MODE = 0.")
2172 elseif (wt_mode == 1)
then
2173 tmpstr = rh18_string
2174 call mom_error(warning,
"Use EPBL_VEL_SCALE_SCHEME = REICHL_H18 instead of the archaic EPBL_VEL_SCALE_MODE = 1.")
2175 elseif (wt_mode >= 2)
then
2176 call mom_error(fatal,
"An unrecognized value of the obsolete parameter EPBL_VEL_SCALE_MODE was specified.")
2178 call log_param(param_file, mdl,
"EPBL_VEL_SCALE_SCHEME", tmpstr, &
2179 "Selects the method for translating TKE into turbulent velocities. "//&
2180 "Valid values are: \n"//&
2181 "\t CUBE_ROOT_TKE - A constant times the cube root of remaining TKE. \n"//&
2182 "\t REICHL_H18 - Use the scheme based on a combination of w* and v* as \n"//&
2183 "\t documented in Reichl & Hallberg, 2018.", &
2184 default=root_tke_string)
2185 tmpstr = uppercase(tmpstr)
2186 select case (tmpstr)
2187 case (root_tke_string)
2188 cs%wT_scheme = wt_from_croot_tke
2190 cs%wT_scheme = wt_from_rh18
2192 call mom_mesg(
'energetic_PBL_init: EPBL_VEL_SCALE_SCHEME ="'//trim(tmpstr)//
'"', 0)
2193 call mom_error(fatal,
"energetic_PBL_init: Unrecognized setting "// &
2194 "EPBL_VEL_SCALE_SCHEME = "//trim(tmpstr)//
" found in input file.")
2197 call get_param(param_file, mdl,
"WSTAR_USTAR_COEF", cs%wstar_ustar_coef, &
2198 "A ratio relating the efficiency with which convectively "//&
2199 "released energy is converted to a turbulent velocity, "//&
2200 "relative to mechanically forced TKE. Making this larger "//&
2201 "increases the BL diffusivity", units=
"nondim", default=1.0)
2202 call get_param(param_file, mdl,
"EPBL_VEL_SCALE_FACTOR", cs%vstar_scale_fac, &
2203 "An overall nondimensional scaling factor for wT. "//&
2204 "Making this larger increases the PBL diffusivity.", &
2205 units=
"nondim", default=1.0)
2206 call get_param(param_file, mdl,
"VSTAR_SURF_FAC", cs%vstar_surf_fac,&
2207 "The proportionality times ustar to set vstar at the surface.", &
2208 units=
"nondim", default=1.2)
2211 call get_param(param_file, mdl,
"USE_LA_LI2016", use_la_windsea, &
2212 "A logical to use the Li et al. 2016 (submitted) formula to "//&
2213 "determine the Langmuir number.", units=
"nondim", default=.false.)
2215 if (use_la_windsea)
then
2218 call get_param(param_file, mdl,
"EPBL_LT", cs%USE_LT, &
2219 "A logical to use a LT parameterization.", &
2220 units=
"nondim", default=.false.)
2223 call get_param(param_file, mdl,
"EPBL_LANGMUIR_SCHEME", tmpstr, &
2224 "EPBL_LANGMUIR_SCHEME selects the method for including Langmuir turbulence. "//&
2225 "Valid values are: \n"//&
2226 "\t NONE - Do not do any extra mixing due to Langmuir turbulence \n"//&
2227 "\t RESCALE - Use a multiplicative rescaling of mstar to account for Langmuir turbulence \n"//&
2228 "\t ADDITIVE - Add a Langmuir turblence contribution to mstar to other contributions", &
2229 default=none_string, do_not_log=.true.)
2230 call get_param(param_file, mdl,
"LT_ENHANCE", lt_enhance, default=-1)
2231 if (lt_enhance == 0)
then
2232 tmpstr = none_string
2233 call mom_error(warning,
"Use EPBL_LANGMUIR_SCHEME = NONE instead of the archaic LT_ENHANCE = 0.")
2234 elseif (lt_enhance == 1)
then
2235 call mom_error(fatal,
"You are using a legacy LT_ENHANCE mode in ePBL that has been phased out. "//&
2236 "If you need to use this setting please report this error. Also use "//&
2237 "EPBL_LANGMUIR_SCHEME to specify the scheme for mstar.")
2238 elseif (lt_enhance == 2)
then
2239 tmpstr = rescaled_string
2240 call mom_error(warning,
"Use EPBL_LANGMUIR_SCHEME = RESCALE instead of the archaic LT_ENHANCE = 2.")
2241 elseif (lt_enhance == 3)
then
2242 tmpstr = additive_string
2243 call mom_error(warning,
"Use EPBL_LANGMUIR_SCHEME = ADDITIVE instead of the archaic LT_ENHANCE = 3.")
2244 elseif (lt_enhance > 3)
then
2245 call mom_error(fatal,
"An unrecognized value of the obsolete parameter LT_ENHANCE was specified.")
2247 call log_param(param_file, mdl,
"EPBL_LANGMUIR_SCHEME", tmpstr, &
2248 "EPBL_LANGMUIR_SCHEME selects the method for including Langmuir turbulence. "//&
2249 "Valid values are: \n"//&
2250 "\t NONE - Do not do any extra mixing due to Langmuir turbulence \n"//&
2251 "\t RESCALE - Use a multiplicative rescaling of mstar to account for Langmuir turbulence \n"//&
2252 "\t ADDITIVE - Add a Langmuir turblence contribution to mstar to other contributions", &
2253 default=none_string)
2254 tmpstr = uppercase(tmpstr)
2255 select case (tmpstr)
2257 cs%LT_enhance_form = no_langmuir
2258 case (rescaled_string)
2259 cs%LT_enhance_form = langmuir_rescale
2260 case (additive_string)
2261 cs%LT_enhance_form = langmuir_add
2263 call mom_mesg(
'energetic_PBL_init: EPBL_LANGMUIR_SCHEME ="'//trim(tmpstr)//
'"', 0)
2264 call mom_error(fatal,
"energetic_PBL_init: Unrecognized setting "// &
2265 "EPBL_LANGMUIR_SCHEME = "//trim(tmpstr)//
" found in input file.")
2268 call get_param(param_file, mdl,
"LT_ENHANCE_COEF", cs%LT_ENHANCE_COEF, &
2269 "Coefficient for Langmuir enhancement of mstar", &
2270 units=
"nondim", default=0.447, do_not_log=(cs%LT_enhance_form==no_langmuir))
2271 call get_param(param_file, mdl,
"LT_ENHANCE_EXP", cs%LT_ENHANCE_EXP, &
2272 "Exponent for Langmuir enhancementt of mstar", &
2273 units=
"nondim", default=-1.33, do_not_log=(cs%LT_enhance_form==no_langmuir))
2274 call get_param(param_file, mdl,
"LT_MOD_LAC1", cs%LaC_MLDoEK, &
2275 "Coefficient for modification of Langmuir number due to "//&
2276 "MLD approaching Ekman depth.", &
2277 units=
"nondim", default=-0.87, do_not_log=(cs%LT_enhance_form==no_langmuir))
2278 call get_param(param_file, mdl,
"LT_MOD_LAC2", cs%LaC_MLDoOB_stab, &
2279 "Coefficient for modification of Langmuir number due to "//&
2280 "MLD approaching stable Obukhov depth.", &
2281 units=
"nondim", default=0.0, do_not_log=(cs%LT_enhance_form==no_langmuir))
2282 call get_param(param_file, mdl,
"LT_MOD_LAC3", cs%LaC_MLDoOB_un, &
2283 "Coefficient for modification of Langmuir number due to "//&
2284 "MLD approaching unstable Obukhov depth.", &
2285 units=
"nondim", default=0.0, do_not_log=(cs%LT_enhance_form==no_langmuir))
2286 call get_param(param_file, mdl,
"LT_MOD_LAC4", cs%Lac_EKoOB_stab, &
2287 "Coefficient for modification of Langmuir number due to "//&
2288 "ratio of Ekman to stable Obukhov depth.", &
2289 units=
"nondim", default=0.95, do_not_log=(cs%LT_enhance_form==no_langmuir))
2290 call get_param(param_file, mdl,
"LT_MOD_LAC5", cs%Lac_EKoOB_un, &
2291 "Coefficient for modification of Langmuir number due to "//&
2292 "ratio of Ekman to unstable Obukhov depth.", &
2293 units=
"nondim", default=0.95, do_not_log=(cs%LT_enhance_form==no_langmuir))
2299 cs%ustar_min = 2e-4*cs%omega*(gv%Angstrom_Z + gv%H_to_Z*gv%H_subroundoff)
2300 call log_param(param_file, mdl,
"EPBL_USTAR_MIN", cs%ustar_min*us%Z_to_m*us%s_to_T, &
2301 "The (tiny) minimum friction velocity used within the "//&
2302 "ePBL code, derived from OMEGA and ANGSTROM.", units=
"m s-1")
2306 z3_t3_to_m3_s3 = us%Z_to_m**3 * us%s_to_T**3
2307 cs%id_ML_depth = register_diag_field(
'ocean_model',
'ePBL_h_ML', diag%axesT1, &
2308 time,
'Surface boundary layer depth',
'm', conversion=us%Z_to_m, &
2309 cmor_long_name=
'Ocean Mixed Layer Thickness Defined by Mixing Scheme')
2310 cs%id_TKE_wind = register_diag_field(
'ocean_model',
'ePBL_TKE_wind', diag%axesT1, &
2311 time,
'Wind-stirring source of mixed layer TKE',
'm3 s-3', conversion=z3_t3_to_m3_s3)
2312 cs%id_TKE_MKE = register_diag_field(
'ocean_model',
'ePBL_TKE_MKE', diag%axesT1, &
2313 time,
'Mean kinetic energy source of mixed layer TKE',
'm3 s-3', conversion=z3_t3_to_m3_s3)
2314 cs%id_TKE_conv = register_diag_field(
'ocean_model',
'ePBL_TKE_conv', diag%axesT1, &
2315 time,
'Convective source of mixed layer TKE',
'm3 s-3', conversion=z3_t3_to_m3_s3)
2316 cs%id_TKE_forcing = register_diag_field(
'ocean_model',
'ePBL_TKE_forcing', diag%axesT1, &
2317 time,
'TKE consumed by mixing surface forcing or penetrative shortwave radation '//&
2318 'through model layers',
'm3 s-3', conversion=z3_t3_to_m3_s3)
2319 cs%id_TKE_mixing = register_diag_field(
'ocean_model',
'ePBL_TKE_mixing', diag%axesT1, &
2320 time,
'TKE consumed by mixing that deepens the mixed layer',
'm3 s-3', conversion=z3_t3_to_m3_s3)
2321 cs%id_TKE_mech_decay = register_diag_field(
'ocean_model',
'ePBL_TKE_mech_decay', diag%axesT1, &
2322 time,
'Mechanical energy decay sink of mixed layer TKE',
'm3 s-3', conversion=z3_t3_to_m3_s3)
2323 cs%id_TKE_conv_decay = register_diag_field(
'ocean_model',
'ePBL_TKE_conv_decay', diag%axesT1, &
2324 time,
'Convective energy decay sink of mixed layer TKE',
'm3 s-3', conversion=z3_t3_to_m3_s3)
2325 cs%id_Mixing_Length = register_diag_field(
'ocean_model',
'Mixing_Length', diag%axesTi, &
2326 time,
'Mixing Length that is used',
'm', conversion=us%Z_to_m)
2327 cs%id_Velocity_Scale = register_diag_field(
'ocean_model',
'Velocity_Scale', diag%axesTi, &
2328 time,
'Velocity Scale that is used.',
'm s-1', conversion=us%Z_to_m*us%s_to_T)
2329 cs%id_MSTAR_mix = register_diag_field(
'ocean_model',
'MSTAR', diag%axesT1, &
2330 time,
'Total mstar that is used.',
'nondim')
2333 cs%id_LA = register_diag_field(
'ocean_model',
'LA', diag%axesT1, &
2334 time,
'Langmuir number.',
'nondim')
2335 cs%id_LA_mod = register_diag_field(
'ocean_model',
'LA_MOD', diag%axesT1, &
2336 time,
'Modified Langmuir number.',
'nondim')
2337 cs%id_MSTAR_LT = register_diag_field(
'ocean_model',
'MSTAR_LT', diag%axesT1, &
2338 time,
'Increase in mstar due to Langmuir Turbulence.',
'nondim')
2341 call get_param(param_file, mdl,
"ENABLE_THERMODYNAMICS", use_temperature, &
2342 "If true, temperature and salinity are used as state "//&
2343 "variables.", default=.true.)
2345 if (max(cs%id_TKE_wind, cs%id_TKE_MKE, cs%id_TKE_conv, &
2346 cs%id_TKE_mixing, cs%id_TKE_mech_decay, cs%id_TKE_forcing, &
2347 cs%id_TKE_conv_decay) > 0)
then
2348 call safe_alloc_alloc(cs%diag_TKE_wind, isd, ied, jsd, jed)
2349 call safe_alloc_alloc(cs%diag_TKE_MKE, isd, ied, jsd, jed)
2350 call safe_alloc_alloc(cs%diag_TKE_conv, isd, ied, jsd, jed)
2351 call safe_alloc_alloc(cs%diag_TKE_forcing, isd, ied, jsd, jed)
2352 call safe_alloc_alloc(cs%diag_TKE_mixing, isd, ied, jsd, jed)
2353 call safe_alloc_alloc(cs%diag_TKE_mech_decay, isd, ied, jsd, jed)
2354 call safe_alloc_alloc(cs%diag_TKE_conv_decay, isd, ied, jsd, jed)
2356 cs%TKE_diagnostics = .true.
2358 if (cs%id_Velocity_Scale>0)
call safe_alloc_alloc(cs%Velocity_Scale, isd, ied, jsd, jed, gv%ke+1)
2359 if (cs%id_Mixing_Length>0)
call safe_alloc_alloc(cs%Mixing_Length, isd, ied, jsd, jed, gv%ke+1)
2361 call safe_alloc_alloc(cs%ML_depth, isd, ied, jsd, jed)
2362 if (max(cs%id_mstar_mix, cs%id_LA, cs%id_LA_mod, cs%id_MSTAR_LT ) >0)
then
2363 call safe_alloc_alloc(cs%Mstar_mix, isd, ied, jsd, jed)
2364 call safe_alloc_alloc(cs%LA, isd, ied, jsd, jed)
2365 call safe_alloc_alloc(cs%LA_MOD, isd, ied, jsd, jed)
2366 call safe_alloc_alloc(cs%MSTAR_LT, isd, ied, jsd, jed)
2369 end subroutine energetic_pbl_init
2372 subroutine energetic_pbl_end(CS)
2376 if (.not.
associated(cs))
return
2378 if (
allocated(cs%ML_depth))
deallocate(cs%ML_depth)
2379 if (
allocated(cs%LA))
deallocate(cs%LA)
2380 if (
allocated(cs%LA_MOD))
deallocate(cs%LA_MOD)
2381 if (
allocated(cs%MSTAR_MIX))
deallocate(cs%MSTAR_MIX)
2382 if (
allocated(cs%MSTAR_LT))
deallocate(cs%MSTAR_LT)
2383 if (
allocated(cs%diag_TKE_wind))
deallocate(cs%diag_TKE_wind)
2384 if (
allocated(cs%diag_TKE_MKE))
deallocate(cs%diag_TKE_MKE)
2385 if (
allocated(cs%diag_TKE_conv))
deallocate(cs%diag_TKE_conv)
2386 if (
allocated(cs%diag_TKE_forcing))
deallocate(cs%diag_TKE_forcing)
2387 if (
allocated(cs%diag_TKE_mixing))
deallocate(cs%diag_TKE_mixing)
2388 if (
allocated(cs%diag_TKE_mech_decay))
deallocate(cs%diag_TKE_mech_decay)
2389 if (
allocated(cs%diag_TKE_conv_decay))
deallocate(cs%diag_TKE_conv_decay)
2390 if (
allocated(cs%Mixing_Length))
deallocate(cs%Mixing_Length)
2391 if (
allocated(cs%Velocity_Scale))
deallocate(cs%Velocity_Scale)
2395 end subroutine energetic_pbl_end