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>
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 [R Z3 T-3 ~> W m-2].
171 diag_tke_mke, & !< The resolved KE source of TKE [R Z3 T-3 ~> W m-2].
172 diag_tke_conv, & !< The convective source of TKE [R 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
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)
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) :: &
726 real :: TKE_left_min, TKE_left_max
727 real :: Kddt_h_max, Kddt_h_min
734 real :: vstar_unit_scale
736 logical :: convectively_stable
737 logical :: sfc_connected
739 logical :: sfc_disconnect
749 real :: MLD_guess, MLD_found
764 logical :: OBL_converged
767 real :: Surface_Scale
769 logical :: debug=.false.
772 real :: dPE_debug, mixing_debug, taux2, tauy2
773 real,
dimension(20) :: TKE_left_itt, PE_chg_itt, Kddt_h_itt, dPEa_dKd_itt, MKE_src_itt
774 real,
dimension(SZK_(GV)) :: mech_TKE_k, conv_PErel_k, nstar_k
775 integer,
dimension(SZK_(GV)) :: num_itts
777 integer :: k, nz, itt, max_itt
781 if (.not.
associated(cs))
call mom_error(fatal,
"energetic_PBL: "//&
782 "Module must be initialized before it is used.")
784 debug = .false. ;
if (
allocated(ecd%dT_expect) .or.
allocated(ecd%dS_expect)) debug = .true.
786 h_neglect = gv%H_subroundoff
789 dt__diag = dt ;
if (
present(dt_diag)) dt__diag = dt_diag
790 i_dtdiag = 1.0 / dt__diag
794 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)
795 vstar_unit_scale = us%m_to_Z * us%T_to_s
806 do k=1,nz+1 ; kd(k) = 0.0 ;
enddo
810 dmass = gv%H_to_RZ * h(k)
811 dpres = us%L_to_Z**2 * gv%g_Earth * dmass
812 dt_to_dpe(k) = (dmass * (pres_z(k) + 0.5*dpres)) * dsv_dt(k)
813 ds_to_dpe(k) = (dmass * (pres_z(k) + 0.5*dpres)) * dsv_ds(k)
814 dt_to_dcolht(k) = dmass * dsv_dt(k)
815 ds_to_dcolht(k) = dmass * dsv_ds(k)
817 pres_z(k+1) = pres_z(k) + dpres
821 h_sum = h_neglect ;
do k=1,nz ; h_sum = h_sum + h(k) ;
enddo
822 i_hs = 0.0 ;
if (h_sum > 0.0) i_hs = 1.0 / h_sum
827 hb_hs(k) = h_bot * i_hs
830 mld_output = h(1)*gv%H_to_Z
834 max_mld = 0.0 ;
do k=1,nz ; max_mld = max_mld + h(k)*gv%H_to_Z ;
enddo
839 if (mld_guess <= min_mld) mld_guess = 0.5 * (min_mld + max_mld)
842 obl_converged = .false.
843 do obl_it=1,cs%Max_MLD_Its
845 if (.not. obl_converged)
then
847 if (.not.cs%Use_MLD_iteration) obl_converged = .true.
849 if (debug)
then ; mech_tke_k(:) = 0.0 ; conv_perel_k(:) = 0.0 ;
endif
852 mld_output = h(1)*gv%H_to_Z
853 sfc_connected = .true.
858 h=h, u_h=u, v_h=v, waves=waves)
859 call find_mstar(cs, us, b_flux, u_star, u_star_mean, mld_guess, absf, &
860 mstar_total, langmuir_number=la, convect_langmuir_number=lamod,&
863 call find_mstar(cs, us, b_flux, u_star, u_star_mean, mld_guess, absf, mstar_total)
868 mech_tke = (dt*mstar_total*gv%Rho0) * u_star**3
870 mech_tke = mstar_total * (dt*gv%Rho0* u_star**3)
873 if (cs%TKE_diagnostics)
then
874 ecd%dTKE_conv = 0.0 ; ecd%dTKE_mixing = 0.0
875 ecd%dTKE_MKE = 0.0 ; ecd%dTKE_mech_decay = 0.0 ; ecd%dTKE_conv_decay = 0.0
877 ecd%dTKE_wind = mech_tke * i_dtdiag
878 if (tke_forcing(1) <= 0.0)
then
879 ecd%dTKE_forcing = max(-mech_tke, tke_forcing(1)) * i_dtdiag
882 ecd%dTKE_forcing = cs%nstar*tke_forcing(1) * i_dtdiag
887 if (tke_forcing(1) <= 0.0)
then
888 mech_tke = mech_tke + tke_forcing(1)
889 if (mech_tke < 0.0) mech_tke = 0.0
892 conv_perel = tke_forcing(1)
897 do k=1,nz+1 ; mixvel(k) = 0.0 ; mixlen(k) = 0.0 ;
enddo
900 if ((.not.cs%Use_MLD_iteration) .or. &
901 (cs%transLay_scale >= 1.0) .or. (cs%transLay_scale < 0.0) )
then
903 mixlen_shape(k) = 1.0
905 elseif (mld_guess <= 0.0)
then
906 if (cs%transLay_scale > 0.0)
then ;
do k=1,nz+1
907 mixlen_shape(k) = cs%transLay_scale
908 enddo ;
else ;
do k=1,nz+1
909 mixlen_shape(k) = 1.0
914 i_mld = 1.0 / mld_guess
916 mixlen_shape(1) = 1.0
918 h_rsum = h_rsum + h(k-1)*gv%H_to_Z
919 if (cs%MixLenExponent==2.0)
then
920 mixlen_shape(k) = cs%transLay_scale + (1.0 - cs%transLay_scale) * &
921 (max(0.0, (mld_guess - h_rsum)*i_mld) )**2
923 mixlen_shape(k) = cs%transLay_scale + (1.0 - cs%transLay_scale) * &
924 (max(0.0, (mld_guess - h_rsum)*i_mld) )**cs%MixLenExponent
929 kd(1) = 0.0 ; kddt_h(1) = 0.0
931 dt_to_dpe_a(1) = dt_to_dpe(1) ; dt_to_dcolht_a(1) = dt_to_dcolht(1)
932 ds_to_dpe_a(1) = ds_to_dpe(1) ; ds_to_dcolht_a(1) = ds_to_dcolht(1)
934 htot = h(1) ; uhtot = u(1)*h(1) ; vhtot = v(1)*h(1)
937 mech_tke_k(1) = mech_tke ; conv_perel_k(1) = conv_perel
938 nstar_k(:) = 0.0 ; nstar_k(1) = cs%nstar ; num_itts(:) = -1
949 idecay_len_tke = (cs%TKE_decay * absf / u_star) * gv%H_to_Z
951 if (idecay_len_tke > 0.0) exp_kh = exp(-h(k-1)*idecay_len_tke)
952 if (cs%TKE_diagnostics) &
953 ecd%dTKE_mech_decay = ecd%dTKE_mech_decay + (exp_kh-1.0) * mech_tke * i_dtdiag
954 mech_tke = mech_tke * exp_kh
958 if (tke_forcing(k) > 0.0)
then
959 conv_perel = conv_perel + tke_forcing(k)
960 if (cs%TKE_diagnostics) &
961 ecd%dTKE_forcing = ecd%dTKE_forcing + cs%nstar*tke_forcing(k) * i_dtdiag
965 mech_tke_k(k) = mech_tke ; conv_perel_k(k) = conv_perel
970 if (cs%nstar * conv_perel > 0.0)
then
974 nstar_fc = cs%nstar * conv_perel / (conv_perel + 0.2 * &
975 sqrt(0.5 * dt * gv%Rho0 * (absf*(htot*gv%H_to_Z))**3 * conv_perel))
978 if (debug) nstar_k(k) = nstar_fc
980 tot_tke = mech_tke + nstar_fc * conv_perel
984 if (tke_forcing(k) < 0.0)
then
985 if (tke_forcing(k) + tot_tke < 0.0)
then
987 if (cs%TKE_diagnostics)
then
988 ecd%dTKE_mixing = ecd%dTKE_mixing + tot_tke * i_dtdiag
989 ecd%dTKE_forcing = ecd%dTKE_forcing - tot_tke * i_dtdiag
991 ecd%dTKE_conv_decay = ecd%dTKE_conv_decay + (cs%nstar-nstar_fc) * conv_perel * i_dtdiag
993 tot_tke = 0.0 ; mech_tke = 0.0 ; conv_perel = 0.0
996 tke_reduc = (tot_tke + tke_forcing(k)) / tot_tke
997 if (cs%TKE_diagnostics)
then
998 ecd%dTKE_mixing = ecd%dTKE_mixing - tke_forcing(k) * i_dtdiag
999 ecd%dTKE_forcing = ecd%dTKE_forcing + tke_forcing(k) * i_dtdiag
1000 ecd%dTKE_conv_decay = ecd%dTKE_conv_decay + &
1001 (1.0-tke_reduc)*(cs%nstar-nstar_fc) * conv_perel * i_dtdiag
1003 tot_tke = tke_reduc*tot_tke
1004 mech_tke = tke_reduc*mech_tke
1005 conv_perel = tke_reduc*conv_perel
1010 if (cs%orig_PE_calc)
then
1012 dte_t2 = 0.0 ; dse_t2 = 0.0
1014 dte_t2 = kddt_h(k-1) * ((t0(k-2) - t0(k-1)) + dte(k-2))
1015 dse_t2 = kddt_h(k-1) * ((s0(k-2) - s0(k-1)) + dse(k-2))
1018 dt_h = (gv%Z_to_H**2*dt) / max(0.5*(h(k-1)+h(k)), 1e-15*h_sum)
1026 convectively_stable = ( 0.0 <= &
1027 ( (dt_to_dcolht(k) + dt_to_dcolht(k-1) ) * (t0(k-1)-t0(k)) + &
1028 (ds_to_dcolht(k) + ds_to_dcolht(k-1) ) * (s0(k-1)-s0(k)) ) )
1030 if ((mech_tke + conv_perel) <= 0.0 .and. convectively_stable)
then
1032 tot_tke = 0.0 ; mech_tke = 0.0 ; conv_perel = 0.0
1033 kd(k) = 0.0 ; kddt_h(k) = 0.0
1034 sfc_disconnect = .true.
1044 if (cs%orig_PE_calc)
then
1045 dte(k-1) = b1 * ( dte_t2 )
1046 dse(k-1) = b1 * ( dse_t2 )
1050 dt_to_dpe_a(k) = dt_to_dpe(k)
1051 ds_to_dpe_a(k) = ds_to_dpe(k)
1052 dt_to_dcolht_a(k) = dt_to_dcolht(k)
1053 ds_to_dcolht_a(k) = ds_to_dcolht(k)
1056 sfc_disconnect = .false.
1060 if (cs%orig_PE_calc)
then
1062 dt_km1_t2 = (t0(k)-t0(k-1))
1063 ds_km1_t2 = (s0(k)-s0(k-1))
1065 dt_km1_t2 = (t0(k)-t0(k-1)) - &
1066 (kddt_h(k-1) / hp_a) * ((t0(k-2) - t0(k-1)) + dte(k-2))
1067 ds_km1_t2 = (s0(k)-s0(k-1)) - &
1068 (kddt_h(k-1) / hp_a) * ((s0(k-2) - s0(k-1)) + dse(k-2))
1070 dte_term = dte_t2 + hp_a * (t0(k-1)-t0(k))
1071 dse_term = dse_t2 + hp_a * (s0(k-1)-s0(k))
1074 th_a(k-1) = h(k-1) * t0(k-1) ; sh_a(k-1) = h(k-1) * s0(k-1)
1076 th_a(k-1) = h(k-1) * t0(k-1) + kddt_h(k-1) * te(k-2)
1077 sh_a(k-1) = h(k-1) * s0(k-1) + kddt_h(k-1) * se(k-2)
1079 th_b(k) = h(k) * t0(k) ; sh_b(k) = h(k) * s0(k)
1087 if ((cs%MKE_to_TKE_effic > 0.0) .and. (htot*h(k) > 0.0))
then
1090 dmke_max = (us%L_to_Z**2*gv%H_to_RZ * cs%MKE_to_TKE_effic) * 0.5 * &
1091 (h(k) / ((htot + h(k))*htot)) * &
1092 ((uhtot-u(k)*htot)**2 + (vhtot-v(k)*htot)**2)
1095 mke2_hharm = (htot + h(k) + 2.0*h_neglect) / &
1096 ((htot+h_neglect) * (h(k)+h_neglect))
1105 h_tt = htot + h_tt_min
1106 tke_here = mech_tke + cs%wstar_ustar_coef*conv_perel
1107 if (tke_here > 0.0)
then
1109 vstar = cs%vstar_scale_fac * vstar_unit_scale * (i_dtrho*tke_here)**c1_3
1111 surface_scale = max(0.05, 1.0 - htot/mld_guess)
1112 vstar = cs%vstar_scale_fac * surface_scale * (cs%vstar_surf_fac*u_star + &
1113 vstar_unit_scale * (cs%wstar_ustar_coef*conv_perel*i_dtrho)**c1_3)
1115 hbs_here = gv%H_to_Z * min(hb_hs(k), mixlen_shape(k))
1116 mixlen(k) = max(cs%min_mix_len, ((h_tt*hbs_here)*vstar) / &
1117 ((cs%Ekman_scale_coef * absf) * (h_tt*hbs_here) + vstar))
1120 if (.not.cs%Use_MLD_iteration)
then
1121 kd_guess0 = vstar * cs%vonKar * ((h_tt*hbs_here)*vstar) / &
1122 ((cs%Ekman_scale_coef * absf) * (h_tt*hbs_here) + vstar)
1124 kd_guess0 = vstar * cs%vonKar * mixlen(k)
1127 vstar = 0.0 ; kd_guess0 = 0.0
1130 kddt_h_g0 = kd_guess0 * dt_h
1132 if (cs%orig_PE_calc)
then
1134 dt_km1_t2, ds_km1_t2, dt_to_dpe(k), ds_to_dpe(k), &
1135 dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), &
1136 pres_z(k), dt_to_dcolht(k), ds_to_dcolht(k), &
1137 dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
1138 pe_chg=pe_chg_g0, dpec_dkd=dpea_dkd_g0, dpe_max=pe_chg_max, &
1139 dpec_dkd_0=dpec_dkd_kd0 )
1142 th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
1143 dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), dt_to_dpe(k), ds_to_dpe(k), &
1144 pres_z(k), dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
1145 dt_to_dcolht(k), ds_to_dcolht(k), &
1146 pe_chg=pe_chg_g0, dpec_dkd=dpea_dkd_g0, dpe_max=pe_chg_max, &
1147 dpec_dkd_0=dpec_dkd_kd0 )
1150 mke_src = dmke_max*(1.0 - exp(-kddt_h_g0 * mke2_hharm))
1153 if ((pe_chg_g0 < 0.0) .or. ((vstar == 0.0) .and. (dpec_dkd_kd0 < 0.0)))
then
1155 if (pe_chg_max <= 0.0)
then
1157 tke_here = mech_tke + cs%wstar_ustar_coef*(conv_perel-pe_chg_max)
1158 if (tke_here > 0.0)
then
1160 vstar = cs%vstar_scale_fac * vstar_unit_scale * (i_dtrho*tke_here)**c1_3
1162 surface_scale = max(0.05, 1. - htot/mld_guess)
1163 vstar = cs%vstar_scale_fac * surface_scale * (cs%vstar_surf_fac*u_star + &
1164 vstar_unit_scale * (cs%wstar_ustar_coef*conv_perel*i_dtrho)**c1_3)
1166 hbs_here = gv%H_to_Z * min(hb_hs(k), mixlen_shape(k))
1167 mixlen(k) = max(cs%min_mix_len, ((h_tt*hbs_here)*vstar) / &
1168 ((cs%Ekman_scale_coef * absf) * (h_tt*hbs_here) + vstar))
1169 if (.not.cs%Use_MLD_iteration)
then
1172 kd(k) = vstar * cs%vonKar * ((h_tt*hbs_here)*vstar) / &
1173 ((cs%Ekman_scale_coef * absf) * (h_tt*hbs_here) + vstar)
1175 kd(k) = vstar * cs%vonKar * mixlen(k)
1178 vstar = 0.0 ; kd(k) = 0.0
1182 if (cs%orig_PE_calc)
then
1184 dt_km1_t2, ds_km1_t2, dt_to_dpe(k), ds_to_dpe(k), &
1185 dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), &
1186 pres_z(k), dt_to_dcolht(k), ds_to_dcolht(k), &
1187 dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
1191 th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
1192 dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), dt_to_dpe(k), ds_to_dpe(k), &
1193 pres_z(k), dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
1194 dt_to_dcolht(k), ds_to_dcolht(k), &
1198 if (dpe_conv > 0.0)
then
1199 kd(k) = kd_guess0 ; dpe_conv = pe_chg_g0
1201 mke_src = dmke_max*(1.0 - exp(-(kd(k)*dt_h) * mke2_hharm))
1205 kd(k) = kd_guess0 ; dpe_conv = pe_chg_g0
1208 conv_perel = conv_perel - dpe_conv
1209 mech_tke = mech_tke + mke_src
1210 if (cs%TKE_diagnostics)
then
1211 ecd%dTKE_conv = ecd%dTKE_conv - cs%nstar*dpe_conv * i_dtdiag
1212 ecd%dTKE_MKE = ecd%dTKE_MKE + mke_src * i_dtdiag
1214 if (sfc_connected)
then
1215 mld_output = mld_output + gv%H_to_Z * h(k)
1218 kddt_h(k) = kd(k) * dt_h
1219 elseif (tot_tke + (mke_src - pe_chg_g0) >= 0.0)
then
1223 kddt_h(k) = kddt_h_g0
1226 tot_tke = tot_tke + mke_src
1228 if (tot_tke > 0.0) tke_reduc = (tot_tke - pe_chg_g0) / tot_tke
1229 if (cs%TKE_diagnostics)
then
1230 ecd%dTKE_mixing = ecd%dTKE_mixing - pe_chg_g0 * i_dtdiag
1231 ecd%dTKE_MKE = ecd%dTKE_MKE + mke_src * i_dtdiag
1232 ecd%dTKE_conv_decay = ecd%dTKE_conv_decay + &
1233 (1.0-tke_reduc)*(cs%nstar-nstar_fc) * conv_perel * i_dtdiag
1235 tot_tke = tke_reduc*tot_tke
1236 mech_tke = tke_reduc*(mech_tke + mke_src)
1237 conv_perel = tke_reduc*conv_perel
1238 if (sfc_connected)
then
1239 mld_output = mld_output + gv%H_to_Z * h(k)
1242 elseif (tot_tke == 0.0)
then
1244 kd(k) = 0.0 ; kddt_h(k) = 0.0
1245 tot_tke = 0.0 ; conv_perel = 0.0 ; mech_tke = 0.0
1246 sfc_disconnect = .true.
1250 kddt_h_max = kddt_h_g0 ; kddt_h_min = 0.0
1251 tke_left_max = tot_tke + (mke_src - pe_chg_g0)
1252 tke_left_min = tot_tke
1256 kddt_h_guess = tot_tke * kddt_h_max / max( pe_chg_g0 - mke_src, &
1257 kddt_h_max * (dpec_dkd_kd0 - dmke_max * mke2_hharm) )
1264 tke_left_itt(:) = 0.0 ; dpea_dkd_itt(:) = 0.0 ; pe_chg_itt(:) = 0.0
1265 mke_src_itt(:) = 0.0 ; kddt_h_itt(:) = 0.0
1268 if (cs%orig_PE_calc)
then
1270 dt_km1_t2, ds_km1_t2, dt_to_dpe(k), ds_to_dpe(k), &
1271 dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), &
1272 pres_z(k), dt_to_dcolht(k), ds_to_dcolht(k), &
1273 dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
1274 pe_chg=pe_chg, dpec_dkd=dpec_dkd )
1277 th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
1278 dt_to_dpe_a(k-1), ds_to_dpe_a(k-1), dt_to_dpe(k), ds_to_dpe(k), &
1279 pres_z(k), dt_to_dcolht_a(k-1), ds_to_dcolht_a(k-1), &
1280 dt_to_dcolht(k), ds_to_dcolht(k), &
1283 mke_src = dmke_max * (1.0 - exp(-mke2_hharm * kddt_h_guess))
1284 dmke_src_dk = dmke_max * mke2_hharm * exp(-mke2_hharm * kddt_h_guess)
1286 tke_left = tot_tke + (mke_src - pe_chg)
1287 if (debug .and. itt<=20)
then
1288 kddt_h_itt(itt) = kddt_h_guess ; mke_src_itt(itt) = mke_src
1289 pe_chg_itt(itt) = pe_chg ; dpea_dkd_itt(itt) = dpec_dkd
1290 tke_left_itt(itt) = tke_left
1294 if (tke_left >= 0.0)
then
1295 kddt_h_min = kddt_h_guess ; tke_left_min = tke_left
1297 kddt_h_max = kddt_h_guess ; tke_left_max = tke_left
1303 if (dpec_dkd - dmke_src_dk <= 0.0)
then
1306 dkddt_h_newt = tke_left / (dpec_dkd - dmke_src_dk)
1307 kddt_h_newt = kddt_h_guess + dkddt_h_newt
1308 if ((kddt_h_newt > kddt_h_max) .or. (kddt_h_newt < kddt_h_min)) &
1313 kddt_h_next = kddt_h_guess + dkddt_h_newt
1314 dkddt_h = dkddt_h_newt
1316 kddt_h_next = (tke_left_max * kddt_h_min - kddt_h_max * tke_left_min) / &
1317 (tke_left_max - tke_left_min)
1318 dkddt_h = kddt_h_next - kddt_h_guess
1321 if ((abs(dkddt_h) < 1e-9*kddt_h_guess) .or. (itt==max_itt))
then
1323 if (debug) num_itts(k) = itt
1326 kddt_h_guess = kddt_h_next
1329 kd(k) = kddt_h_guess / dt_h ; kddt_h(k) = kd(k) * dt_h
1332 if (cs%TKE_diagnostics)
then
1333 ecd%dTKE_mixing = ecd%dTKE_mixing - (tot_tke + mke_src) * i_dtdiag
1334 ecd%dTKE_MKE = ecd%dTKE_MKE + mke_src * i_dtdiag
1335 ecd%dTKE_conv_decay = ecd%dTKE_conv_decay + &
1336 (cs%nstar-nstar_fc) * conv_perel * i_dtdiag
1339 if (sfc_connected) mld_output = mld_output + &
1340 (pe_chg / (pe_chg_g0)) * gv%H_to_Z * h(k)
1342 tot_tke = 0.0 ; mech_tke = 0.0 ; conv_perel = 0.0
1343 sfc_disconnect = .true.
1346 kddt_h(k) = kd(k) * dt_h
1349 b1 = 1.0 / (hp_a + kddt_h(k))
1350 c1(k) = kddt_h(k) * b1
1351 if (cs%orig_PE_calc)
then
1352 dte(k-1) = b1 * ( kddt_h(k)*(t0(k)-t0(k-1)) + dte_t2 )
1353 dse(k-1) = b1 * ( kddt_h(k)*(s0(k)-s0(k-1)) + dse_t2 )
1356 hp_a = h(k) + (hp_a * b1) * kddt_h(k)
1357 dt_to_dpe_a(k) = dt_to_dpe(k) + c1(k)*dt_to_dpe_a(k-1)
1358 ds_to_dpe_a(k) = ds_to_dpe(k) + c1(k)*ds_to_dpe_a(k-1)
1359 dt_to_dcolht_a(k) = dt_to_dcolht(k) + c1(k)*dt_to_dcolht_a(k-1)
1360 ds_to_dcolht_a(k) = ds_to_dcolht(k) + c1(k)*ds_to_dcolht_a(k-1)
1365 if (sfc_disconnect)
then
1370 sfc_connected = .false.
1372 uhtot = uhtot + u(k)*h(k)
1373 vhtot = vhtot + v(k)*h(k)
1379 te(1) = b1*(h(1)*t0(1))
1380 se(1) = b1*(h(1)*s0(1))
1382 te(k-1) = b1 * (h(k-1) * t0(k-1) + kddt_h(k-1) * te(k-2))
1383 se(k-1) = b1 * (h(k-1) * s0(k-1) + kddt_h(k-1) * se(k-2))
1392 te(nz) = b1 * (h(nz) * t0(nz) + kddt_h(nz) * te(nz-1))
1393 se(nz) = b1 * (h(nz) * s0(nz) + kddt_h(nz) * se(nz-1))
1394 ecd%dT_expect(nz) = te(nz) - t0(nz) ; ecd%dS_expect(nz) = se(nz) - s0(nz)
1396 te(k) = te(k) + c1(k+1)*te(k+1)
1397 se(k) = se(k) + c1(k+1)*se(k+1)
1398 ecd%dT_expect(k) = te(k) - t0(k) ; ecd%dS_expect(k) = se(k) - s0(k)
1403 dpe_debug = dpe_debug + (dt_to_dpe(k) * (te(k) - t0(k)) + &
1404 ds_to_dpe(k) * (se(k) - s0(k)))
1406 mixing_debug = dpe_debug * i_dtdiag
1417 mld_found = mld_output
1418 if (mld_found - cs%MLD_tol > mld_guess)
then
1420 elseif (abs(mld_guess - mld_found) < cs%MLD_tol)
then
1421 obl_converged = .true.
1428 mld_guess = 0.5*(min_mld + max_mld)
1432 ecd%LA = la ; ecd%LAmod = lamod ; ecd%mstar = mstar_total ; ecd%mstar_LT = mstar_lt
1434 ecd%LA = 0.0 ; ecd%LAmod = 0.0 ; ecd%mstar = mstar_total ; ecd%mstar_LT = 0.0
1443 subroutine find_pe_chg(Kddt_h0, dKddt_h, hp_a, hp_b, Th_a, Sh_a, Th_b, Sh_b, &
1444 dT_to_dPE_a, dS_to_dPE_a, dT_to_dPE_b, dS_to_dPE_b, &
1445 pres_Z, dT_to_dColHt_a, dS_to_dColHt_a, dT_to_dColHt_b, dS_to_dColHt_b, &
1446 PE_chg, dPEc_dKd, dPE_max, dPEc_dKd_0, PE_ColHt_cor)
1447 real,
intent(in) :: Kddt_h0
1450 real,
intent(in) :: dKddt_h
1453 real,
intent(in) :: hp_a
1457 real,
intent(in) :: hp_b
1461 real,
intent(in) :: Th_a
1464 real,
intent(in) :: Sh_a
1467 real,
intent(in) :: Th_b
1470 real,
intent(in) :: Sh_b
1473 real,
intent(in) :: dT_to_dPE_a
1477 real,
intent(in) :: dS_to_dPE_a
1481 real,
intent(in) :: dT_to_dPE_b
1485 real,
intent(in) :: dS_to_dPE_b
1489 real,
intent(in) :: pres_Z
1492 real,
intent(in) :: dT_to_dColHt_a
1496 real,
intent(in) :: dS_to_dColHt_a
1500 real,
intent(in) :: dT_to_dColHt_b
1504 real,
intent(in) :: dS_to_dColHt_b
1509 real,
optional,
intent(out) :: PE_chg
1511 real,
optional,
intent(out) :: dPEc_dKd
1513 real,
optional,
intent(out) :: dPE_max
1516 real,
optional,
intent(out) :: dPEc_dKd_0
1518 real,
optional,
intent(out) :: PE_ColHt_cor
1542 bdt1 = hp_a * hp_b + kddt_h0 * hps
1543 dt_c = hp_a * th_b - hp_b * th_a
1544 ds_c = hp_a * sh_b - hp_b * sh_a
1545 pec_core = hp_b * (dt_to_dpe_a * dt_c + ds_to_dpe_a * ds_c) - &
1546 hp_a * (dt_to_dpe_b * dt_c + ds_to_dpe_b * ds_c)
1547 colht_core = hp_b * (dt_to_dcolht_a * dt_c + ds_to_dcolht_a * ds_c) - &
1548 hp_a * (dt_to_dcolht_b * dt_c + ds_to_dcolht_b * ds_c)
1550 if (
present(pe_chg))
then
1553 y1_3 = dkddt_h / (bdt1 * (bdt1 + dkddt_h * hps))
1554 pe_chg = pec_core * y1_3
1555 colht_chg = colht_core * y1_3
1556 if (colht_chg < 0.0) pe_chg = pe_chg - pres_z * colht_chg
1557 if (
present(pe_colht_cor)) pe_colht_cor = -pres_z * min(colht_chg, 0.0)
1558 elseif (
present(pe_colht_cor))
then
1559 y1_3 = dkddt_h / (bdt1 * (bdt1 + dkddt_h * hps))
1560 pe_colht_cor = -pres_z * min(colht_core * y1_3, 0.0)
1563 if (
present(dpec_dkd))
then
1565 y1_4 = 1.0 / (bdt1 + dkddt_h * hps)**2
1566 dpec_dkd = pec_core * y1_4
1567 colht_chg = colht_core * y1_4
1568 if (colht_chg < 0.0) dpec_dkd = dpec_dkd - pres_z * colht_chg
1571 if (
present(dpe_max))
then
1573 y1_3 = 1.0 / (bdt1 * hps)
1574 dpe_max = pec_core * y1_3
1575 colht_chg = colht_core * y1_3
1576 if (colht_chg < 0.0) dpe_max = dpe_max - pres_z * colht_chg
1579 if (
present(dpec_dkd_0))
then
1581 y1_4 = 1.0 / bdt1**2
1582 dpec_dkd_0 = pec_core * y1_4
1583 colht_chg = colht_core * y1_4
1584 if (colht_chg < 0.0) dpec_dkd_0 = dpec_dkd_0 - pres_z * colht_chg
1593 dT_km1_t2, dS_km1_t2, dT_to_dPE_k, dS_to_dPE_k, &
1594 dT_to_dPEa, dS_to_dPEa, pres_Z, dT_to_dColHt_k, &
1595 dS_to_dColHt_k, dT_to_dColHta, dS_to_dColHta, &
1596 PE_chg, dPEc_dKd, dPE_max, dPEc_dKd_0)
1597 real,
intent(in) :: Kddt_h
1600 real,
intent(in) :: h_k
1601 real,
intent(in) :: b_den_1
1605 real,
intent(in) :: dTe_term
1607 real,
intent(in) :: dSe_term
1609 real,
intent(in) :: dT_km1_t2
1611 real,
intent(in) :: dS_km1_t2
1613 real,
intent(in) :: pres_Z
1616 real,
intent(in) :: dT_to_dPE_k
1620 real,
intent(in) :: dS_to_dPE_k
1624 real,
intent(in) :: dT_to_dPEa
1628 real,
intent(in) :: dS_to_dPEa
1632 real,
intent(in) :: dT_to_dColHt_k
1636 real,
intent(in) :: dS_to_dColHt_k
1640 real,
intent(in) :: dT_to_dColHta
1644 real,
intent(in) :: dS_to_dColHta
1649 real,
optional,
intent(out) :: PE_chg
1651 real,
optional,
intent(out) :: dPEc_dKd
1653 real,
optional,
intent(out) :: dPE_max
1656 real,
optional,
intent(out) :: dPEc_dKd_0
1673 real :: dT_k, dT_km1
1674 real :: dS_k, dS_km1
1677 real :: ddT_k_dKd, ddT_km1_dKd
1678 real :: ddS_k_dKd, ddS_km1_dKd
1680 b1 = 1.0 / (b_den_1 + kddt_h)
1688 i_kr_denom = 1.0 / (h_k*b_den_1 + (b_den_1 + h_k)*kddt_h)
1690 dt_k = (kddt_h*i_kr_denom) * dte_term
1691 ds_k = (kddt_h*i_kr_denom) * dse_term
1693 if (
present(pe_chg))
then
1696 dt_km1 = b1kd * ( dt_k + dt_km1_t2 )
1697 ds_km1 = b1kd * ( ds_k + ds_km1_t2 )
1698 pe_chg = (dt_to_dpe_k * dt_k + dt_to_dpea * dt_km1) + &
1699 (ds_to_dpe_k * ds_k + ds_to_dpea * ds_km1)
1700 colht_chg = (dt_to_dcolht_k * dt_k + dt_to_dcolhta * dt_km1) + &
1701 (ds_to_dcolht_k * ds_k + ds_to_dcolhta * ds_km1)
1702 if (colht_chg < 0.0) pe_chg = pe_chg - pres_z * colht_chg
1705 if (
present(dpec_dkd))
then
1707 dkr_dkd = (h_k*b_den_1) * i_kr_denom**2
1709 ddt_k_dkd = dkr_dkd * dte_term
1710 dds_k_dkd = dkr_dkd * dse_term
1711 ddt_km1_dkd = (b1**2 * b_den_1) * ( dt_k + dt_km1_t2 ) + b1kd * ddt_k_dkd
1712 dds_km1_dkd = (b1**2 * b_den_1) * ( ds_k + ds_km1_t2 ) + b1kd * dds_k_dkd
1715 dpec_dkd = (dt_to_dpe_k * ddt_k_dkd + dt_to_dpea * ddt_km1_dkd) + &
1716 (ds_to_dpe_k * dds_k_dkd + ds_to_dpea * dds_km1_dkd)
1717 dcolht_dkd = (dt_to_dcolht_k * ddt_k_dkd + dt_to_dcolhta * ddt_km1_dkd) + &
1718 (ds_to_dcolht_k * dds_k_dkd + ds_to_dcolhta * dds_km1_dkd)
1719 if (dcolht_dkd < 0.0) dpec_dkd = dpec_dkd - pres_z * dcolht_dkd
1722 if (
present(dpe_max))
then
1724 dpe_max = (dt_to_dpea * dt_km1_t2 + ds_to_dpea * ds_km1_t2) + &
1725 ((dt_to_dpe_k + dt_to_dpea) * dte_term + &
1726 (ds_to_dpe_k + ds_to_dpea) * dse_term) / (b_den_1 + h_k)
1727 dcolht_max = (dt_to_dcolhta * dt_km1_t2 + ds_to_dcolhta * ds_km1_t2) + &
1728 ((dt_to_dcolht_k + dt_to_dcolhta) * dte_term + &
1729 (ds_to_dcolht_k + ds_to_dcolhta) * dse_term) / (b_den_1 + h_k)
1730 if (dcolht_max < 0.0) dpe_max = dpe_max - pres_z*dcolht_max
1733 if (
present(dpec_dkd_0))
then
1735 dpec_dkd_0 = (dt_to_dpea * dt_km1_t2 + ds_to_dpea * ds_km1_t2) / (b_den_1) + &
1736 (dt_to_dpe_k * dte_term + ds_to_dpe_k * dse_term) / (h_k*b_den_1)
1737 dcolht_dkd = (dt_to_dcolhta * dt_km1_t2 + ds_to_dcolhta * ds_km1_t2) / (b_den_1) + &
1738 (dt_to_dcolht_k * dte_term + ds_to_dcolht_k * dse_term) / (h_k*b_den_1)
1739 if (dcolht_dkd < 0.0) dpec_dkd_0 = dpec_dkd_0 - pres_z*dcolht_dkd
1745 subroutine find_mstar(CS, US, Buoyancy_Flux, UStar, UStar_Mean,&
1746 BLD, Abs_Coriolis, MStar, Langmuir_Number,&
1747 MStar_LT, Convect_Langmuir_Number)
1750 real,
intent(in) :: UStar
1751 real,
intent(in) :: UStar_Mean
1752 real,
intent(in) :: Abs_Coriolis
1753 real,
intent(in) :: Buoyancy_Flux
1754 real,
intent(in) :: BLD
1755 real,
intent(out) :: Mstar
1756 real,
optional,
intent(in) :: Langmuir_Number
1757 real,
optional,
intent(out) :: MStar_LT
1758 real,
optional,
intent(out) :: Convect_Langmuir_number
1762 real :: MSCR_term1, MSCR_term2
1763 real :: MStar_Conv_Red
1764 real :: MStar_S, MStar_N
1771 mstar = cs%Fixed_MStar
1775 if (cs%answers_2018)
then
1777 mstar_s = cs%MStar_coef*sqrt(max(0.0,buoyancy_flux) / ustar**2 / &
1778 (abs_coriolis + 1.e-10*us%T_to_s) )
1780 mstar_n = cs%C_Ek * log( max( 1., ustar / (abs_coriolis + 1.e-10*us%T_to_s) / bld ) )
1783 mstar_s = cs%MSTAR_COEF*sqrt(max(0.0, buoyancy_flux) / (ustar**2 * max(abs_coriolis, 1.e-20*us%T_to_s)))
1786 if (ustar > abs_coriolis * bld) mstar_n = cs%C_EK * log(ustar / (abs_coriolis * bld))
1790 mstar = max(mstar_s, min(1.25, mstar_n))
1791 if (cs%MStar_Cap > 0.0) mstar = min( cs%MStar_Cap,mstar )
1793 if (cs%answers_2018)
then
1794 mstar_n = cs%RH18_MStar_cn1 * ( 1.0 - 1.0 / ( 1. + cs%RH18_MStar_cn2 * &
1795 exp( cs%RH18_mstar_CN3 * bld * abs_coriolis / ustar) ) )
1797 msn_term = cs%RH18_MStar_cn2 * exp( cs%RH18_mstar_CN3 * bld * abs_coriolis / ustar)
1798 mstar_n = (cs%RH18_MStar_cn1 * msn_term) / ( 1. + msn_term)
1800 mstar_s = cs%RH18_MStar_CS1 * ( max(0.0, buoyancy_flux)**2 * bld / &
1801 ( ustar**5 * max(abs_coriolis,1.e-20*us%T_to_s) ) )**cs%RH18_mstar_cs2
1802 mstar = mstar_n + mstar_s
1806 if (cs%answers_2018)
then
1807 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) / &
1808 ( (-min(0.0,buoyancy_flux) + 1.e-10*us%T_to_s**3*us%m_to_Z**2) + &
1809 2.0 *mstar * ustar**3 / bld )
1811 mscr_term1 = -bld * min(0.0, buoyancy_flux)
1812 mscr_term2 = 2.0*mstar * ustar**3
1813 if ( abs(mscr_term2) > 0.0)
then
1814 mstar_conv_red = ((1.-cs%mstar_convect_coef) * mscr_term1 + mscr_term2) / (mscr_term1 + mscr_term2)
1816 mstar_conv_red = 1.-cs%mstar_convect_coef
1821 mstar = mstar * mstar_conv_red
1823 if (
present(langmuir_number))
then
1825 call mstar_langmuir(cs, us, abs_coriolis, buoyancy_flux, ustar, bld, langmuir_number, mstar, &
1826 mstar_lt, convect_langmuir_number)
1832 subroutine mstar_langmuir(CS, US, Abs_Coriolis, Buoyancy_Flux, UStar, BLD, Langmuir_Number, &
1833 Mstar, MStar_LT, Convect_Langmuir_Number)
1836 real,
intent(in) :: Abs_Coriolis
1837 real,
intent(in) :: Buoyancy_Flux
1838 real,
intent(in) :: UStar
1839 real,
intent(in) :: BLD
1840 real,
intent(inout) :: Mstar
1841 real,
intent(in) :: Langmuir_Number
1842 real,
intent(out) :: MStar_LT
1843 real,
intent(out) :: Convect_Langmuir_number
1846 real,
parameter :: Max_ratio = 1.0e16
1847 real :: enhance_mstar
1848 real :: mstar_LT_add
1854 real :: Ekman_Obukhov
1856 real :: MLD_Obukhov_stab
1857 real :: Ekman_Obukhov_stab
1858 real :: MLD_Obukhov_un
1859 real :: Ekman_Obukhov_un
1862 enhance_mstar = 1.0 ; mstar_lt_add = 0.0
1866 if (cs%answers_2018)
then
1867 il_ekman = abs_coriolis / ustar
1868 il_obukhov = buoyancy_flux*cs%vonkar / ustar**3
1869 ekman_obukhov_stab = abs(max(0., il_obukhov / (il_ekman + 1.e-10*us%Z_to_m)))
1870 ekman_obukhov_un = abs(min(0., il_obukhov / (il_ekman + 1.e-10*us%Z_to_m)))
1871 mld_obukhov_stab = abs(max(0., bld*il_obukhov))
1872 mld_obukhov_un = abs(min(0., bld*il_obukhov))
1873 mld_ekman = abs( bld*il_ekman )
1875 ekman_obukhov = max_ratio ; mld_obukhov = max_ratio ; mld_ekman = max_ratio
1876 i_f = 0.0 ;
if (abs(abs_coriolis) > 0.0) i_f = 1.0 / abs_coriolis
1877 i_ustar = 0.0 ;
if (abs(ustar) > 0.0) i_ustar = 1.0 / ustar
1878 if (abs(buoyancy_flux*cs%vonkar) < max_ratio*(abs_coriolis * ustar**2)) &
1879 ekman_obukhov = abs(buoyancy_flux*cs%vonkar) * (i_f * i_ustar**2)
1880 if (abs(bld*buoyancy_flux*cs%vonkar) < max_ratio*ustar**3) &
1881 mld_obukhov = abs(bld*buoyancy_flux*cs%vonkar) * i_ustar**3
1882 if (bld*abs_coriolis < max_ratio*ustar) &
1883 mld_ekman = bld*abs_coriolis * i_ustar
1885 if (buoyancy_flux > 0.0)
then
1886 ekman_obukhov_stab = ekman_obukhov ; ekman_obukhov_un = 0.0
1887 mld_obukhov_stab = mld_obukhov ; mld_obukhov_un = 0.0
1889 ekman_obukhov_un = ekman_obukhov ; ekman_obukhov_stab = 0.0
1890 mld_obukhov_un = mld_obukhov ; mld_obukhov_stab = 0.0
1897 convect_langmuir_number = langmuir_number * &
1898 ( (1.0 + max(-0.5, cs%LaC_MLDoEK * mld_ekman)) + &
1899 ((cs%LaC_EKoOB_stab * ekman_obukhov_stab + cs%LaC_EKoOB_un * ekman_obukhov_un) + &
1900 (cs%LaC_MLDoOB_stab * mld_obukhov_stab + cs%LaC_MLDoOB_un * mld_obukhov_un)) )
1904 enhance_mstar = min(cs%Max_Enhance_M, &
1905 (1. + cs%LT_ENHANCE_COEF * convect_langmuir_number**cs%LT_ENHANCE_EXP) )
1908 mstar_lt_add = cs%LT_ENHANCE_COEF * convect_langmuir_number**cs%LT_ENHANCE_EXP
1912 mstar_lt = (enhance_mstar - 1.0)*mstar + mstar_lt_add
1913 mstar = mstar*enhance_mstar + mstar_lt_add
1923 real,
dimension(SZI_(G),SZJ_(G)),
intent(out) :: mld
1924 real,
optional,
intent(in) :: m_to_mld_units
1930 scale = us%Z_to_m ;
if (
present(m_to_mld_units)) scale = scale * m_to_mld_units
1932 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
1933 mld(i,j) = scale*cs%ML_Depth(i,j)
1941 type(time_type),
target,
intent(in) :: time
1946 type(
diag_ctrl),
target,
intent(inout) :: diag
1951 # include "version_variable.h"
1952 character(len=40) :: mdl =
"MOM_energetic_PBL"
1953 character(len=20) :: tmpstr
1954 real :: omega_frac_dflt
1955 real :: r_z3_t3_to_kg_s3
1956 integer :: isd, ied, jsd, jed
1957 integer :: mstar_mode, lt_enhance, wt_mode
1958 logical :: default_2018_answers
1959 logical :: use_temperature, use_omega
1960 logical :: use_la_windsea
1961 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1963 if (
associated(cs))
then
1964 call mom_error(warning,
"mixedlayer_init called with an associated"//&
1965 "associated control structure.")
1967 else ;
allocate(cs) ;
endif
1977 call get_param(param_file, mdl,
"OMEGA", cs%omega, &
1978 "The rotation rate of the earth.", units=
"s-1", &
1979 default=7.2921e-5, scale=us%T_to_S)
1980 call get_param(param_file, mdl,
"ML_USE_OMEGA", use_omega, &
1981 "If true, use the absolute rotation rate instead of the "//&
1982 "vertical component of rotation when setting the decay "//&
1983 "scale for turbulence.", default=.false., do_not_log=.true.)
1984 omega_frac_dflt = 0.0
1986 call mom_error(warning,
"ML_USE_OMEGA is depricated; use ML_OMEGA_FRAC=1.0 instead.")
1987 omega_frac_dflt = 1.0
1989 call get_param(param_file, mdl,
"ML_OMEGA_FRAC", cs%omega_frac, &
1990 "When setting the decay scale for turbulence, use this "//&
1991 "fraction of the absolute rotation rate blended with the "//&
1992 "local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).", &
1993 units=
"nondim", default=omega_frac_dflt)
1994 call get_param(param_file, mdl,
"EKMAN_SCALE_COEF", cs%Ekman_scale_coef, &
1995 "A nondimensional scaling factor controlling the inhibition "//&
1996 "of the diffusive length scale by rotation. Making this larger "//&
1997 "decreases the PBL diffusivity.", units=
"nondim", default=1.0)
1998 call get_param(param_file, mdl,
"DEFAULT_2018_ANSWERS", default_2018_answers, &
1999 "This sets the default value for the various _2018_ANSWERS parameters.", &
2001 call get_param(param_file, mdl,
"EPBL_2018_ANSWERS", cs%answers_2018, &
2002 "If true, use the order of arithmetic and expressions that recover the "//&
2003 "answers from the end of 2018. Otherwise, use updated and more robust "//&
2004 "forms of the same expressions.", default=default_2018_answers)
2007 call get_param(param_file, mdl,
"EPBL_ORIGINAL_PE_CALC", cs%orig_PE_calc, &
2008 "If true, the ePBL code uses the original form of the "//&
2009 "potential energy change code. Otherwise, the newer "//&
2010 "version that can work with successive increments to the "//&
2011 "diffusivity in upward or downward passes is used.", default=.true.)
2013 call get_param(param_file, mdl,
"MKE_TO_TKE_EFFIC", cs%MKE_to_TKE_effic, &
2014 "The efficiency with which mean kinetic energy released "//&
2015 "by mechanically forced entrainment of the mixed layer "//&
2016 "is converted to turbulent kinetic energy.", units=
"nondim", &
2018 call get_param(param_file, mdl,
"TKE_DECAY", cs%TKE_decay, &
2019 "TKE_DECAY relates the vertical rate of decay of the "//&
2020 "TKE available for mechanical entrainment to the natural "//&
2021 "Ekman depth.", units=
"nondim", default=2.5)
2025 call get_param(param_file, mdl,
"EPBL_MSTAR_SCHEME", tmpstr, &
2026 "EPBL_MSTAR_SCHEME selects the method for setting mstar. Valid values are: \n"//&
2027 "\t CONSTANT - Use a fixed mstar given by MSTAR \n"//&
2028 "\t OM4 - Use L_Ekman/L_Obukhov in the sabilizing limit, as in OM4 \n"//&
2029 "\t REICHL_H18 - Use the scheme documented in Reichl & Hallberg, 2018.", &
2031 call get_param(param_file, mdl,
"MSTAR_MODE", mstar_mode, default=-1)
2032 if (mstar_mode == 0)
then
2034 call mom_error(warning,
"Use EPBL_MSTAR_SCHEME = CONSTANT instead of the archaic MSTAR_MODE = 0.")
2035 elseif (mstar_mode == 1)
then
2036 call mom_error(fatal,
"You are using a legacy mstar mode in ePBL that has been phased out. "//&
2037 "If you need to use this setting please report this error. Also use "//&
2038 "EPBL_MSTAR_SCHEME to specify the scheme for mstar.")
2039 elseif (mstar_mode == 2)
then
2041 call mom_error(warning,
"Use EPBL_MSTAR_SCHEME = OM4 instead of the archaic MSTAR_MODE = 2.")
2042 elseif (mstar_mode == 3)
then
2044 call mom_error(warning,
"Use EPBL_MSTAR_SCHEME = REICHL_H18 instead of the archaic MSTAR_MODE = 3.")
2045 elseif (mstar_mode > 3)
then
2046 call mom_error(fatal,
"An unrecognized value of the obsolete parameter MSTAR_MODE was specified.")
2048 call log_param(param_file, mdl,
"EPBL_MSTAR_SCHEME", tmpstr, &
2049 "EPBL_MSTAR_SCHEME selects the method for setting mstar. Valid values are: \n"//&
2050 "\t CONSTANT - Use a fixed mstar given by MSTAR \n"//&
2051 "\t OM4 - Use L_Ekman/L_Obukhov in the sabilizing limit, as in OM4 \n"//&
2052 "\t REICHL_H18 - Use the scheme documented in Reichl & Hallberg, 2018.", &
2055 select case (tmpstr)
2063 call mom_mesg(
'energetic_PBL_init: EPBL_MSTAR_SCHEME ="'//trim(tmpstr)//
'"', 0)
2064 call mom_error(fatal,
"energetic_PBL_init: Unrecognized setting "// &
2065 "EPBL_MSTAR_SCHEME = "//trim(tmpstr)//
" found in input file.")
2068 call get_param(param_file, mdl,
"MSTAR", cs%fixed_mstar, &
2069 "The ratio of the friction velocity cubed to the TKE input to the "//&
2070 "mixed layer. This option is used if EPBL_MSTAR_SCHEME = CONSTANT.", &
2071 units=
"nondim", default=1.2, do_not_log=(cs%mstar_scheme/=
use_fixed_mstar))
2072 call get_param(param_file, mdl,
"MSTAR_CAP", cs%mstar_cap, &
2073 "If this value is positive, it sets the maximum value of mstar "//&
2074 "allowed in ePBL. (This is not used if EPBL_MSTAR_SCHEME = CONSTANT).", &
2075 units=
"nondim", default=-1.0, do_not_log=(cs%mstar_scheme==
use_fixed_mstar))
2077 call get_param(param_file, mdl,
"MSTAR2_COEF1", cs%MSTAR_COEF, &
2078 "Coefficient in computing mstar when rotation and stabilizing "//&
2079 "effects are both important (used if EPBL_MSTAR_SCHEME = OM4).", &
2080 units=
"nondim", default=0.3, do_not_log=(cs%mstar_scheme/=
mstar_from_ekman))
2081 call get_param(param_file, mdl,
"MSTAR2_COEF2", cs%C_EK, &
2082 "Coefficient in computing mstar when only rotation limits "// &
2083 "the total mixing (used if EPBL_MSTAR_SCHEME = OM4)", &
2084 units=
"nondim", default=0.085, do_not_log=(cs%mstar_scheme/=
mstar_from_ekman))
2086 call get_param(param_file, mdl,
"RH18_MSTAR_CN1", cs%RH18_mstar_cn1,&
2087 "MSTAR_N coefficient 1 (outter-most coefficient for fit). "//&
2088 "The value of 0.275 is given in RH18. Increasing this "//&
2089 "coefficient increases MSTAR for all values of Hf/ust, but more "//&
2090 "effectively at low values (weakly developed OSBLs).", &
2091 units=
"nondim", default=0.275, do_not_log=(cs%mstar_scheme/=
mstar_from_rh18))
2092 call get_param(param_file, mdl,
"RH18_MSTAR_CN2", cs%RH18_mstar_cn2,&
2093 "MSTAR_N coefficient 2 (coefficient outside of exponential decay). "//&
2094 "The value of 8.0 is given in RH18. Increasing this coefficient "//&
2095 "increases MSTAR for all values of HF/ust, with a much more even "//&
2096 "effect across a wide range of Hf/ust than CN1.", &
2097 units=
"nondim", default=8.0, do_not_log=(cs%mstar_scheme/=
mstar_from_rh18))
2098 call get_param(param_file, mdl,
"RH18_MSTAR_CN3", cs%RH18_mstar_CN3,&
2099 "MSTAR_N coefficient 3 (exponential decay coefficient). "//&
2100 "The value of -5.0 is given in RH18. Increasing this increases how "//&
2101 "quickly the value of MSTAR decreases as Hf/ust increases.", &
2102 units=
"nondim", default=-5.0, do_not_log=(cs%mstar_scheme/=
mstar_from_rh18))
2103 call get_param(param_file, mdl,
"RH18_MSTAR_CS1", cs%RH18_mstar_cs1,&
2104 "MSTAR_S coefficient for RH18 in stabilizing limit. "//&
2105 "The value of 0.2 is given in RH18 and increasing it increases "//&
2106 "MSTAR in the presence of a stabilizing surface buoyancy flux.", &
2107 units=
"nondim", default=0.2, do_not_log=(cs%mstar_scheme/=
mstar_from_rh18))
2108 call get_param(param_file, mdl,
"RH18_MSTAR_CS2", cs%RH18_mstar_cs2,&
2109 "MSTAR_S exponent for RH18 in stabilizing limit. "//&
2110 "The value of 0.4 is given in RH18 and increasing it increases MSTAR "//&
2111 "exponentially in the presence of a stabilizing surface buoyancy flux.", &
2112 units=
"nondim", default=0.4, do_not_log=(cs%mstar_scheme/=
mstar_from_rh18))
2116 call get_param(param_file, mdl,
"NSTAR", cs%nstar, &
2117 "The portion of the buoyant potential energy imparted by "//&
2118 "surface fluxes that is available to drive entrainment "//&
2119 "at the base of mixed layer when that energy is positive.", &
2120 units=
"nondim", default=0.2)
2121 call get_param(param_file, mdl,
"MSTAR_CONV_ADJ", cs%mstar_convect_coef, &
2122 "Coefficient used for reducing mstar during convection "//&
2123 "due to reduction of stable density gradient.", &
2124 units=
"nondim", default=0.0)
2128 call get_param(param_file, mdl,
"USE_MLD_ITERATION", cs%Use_MLD_iteration, &
2129 "A logical that specifies whether or not to use the "//&
2130 "distance to the bottom of the actively turbulent boundary "//&
2131 "layer to help set the EPBL length scale.", default=.false.)
2132 call get_param(param_file, mdl,
"EPBL_TRANSITION_SCALE", cs%transLay_scale, &
2133 "A scale for the mixing length in the transition layer "//&
2134 "at the edge of the boundary layer as a fraction of the "//&
2135 "boundary layer thickness.", units=
"nondim", default=0.1)
2136 if ( cs%Use_MLD_iteration .and. abs(cs%transLay_scale-0.5) >= 0.5)
then
2137 call mom_error(fatal,
"If flag USE_MLD_ITERATION is true, then "//&
2138 "EPBL_TRANSITION should be greater than 0 and less than 1.")
2141 call get_param(param_file, mdl,
"MLD_ITERATION_GUESS", cs%MLD_ITERATION_GUESS, &
2142 "A logical that specifies whether or not to use the "//&
2143 "previous timestep MLD as a first guess in the MLD iteration. "//&
2144 "The default is false to facilitate reproducibility.", default=.false.)
2145 call get_param(param_file, mdl,
"EPBL_MLD_TOLERANCE", cs%MLD_tol, &
2146 "The tolerance for the iteratively determined mixed "//&
2147 "layer depth. This is only used with USE_MLD_ITERATION.", &
2148 units=
"meter", default=1.0, scale=us%m_to_Z)
2149 call get_param(param_file, mdl,
"EPBL_MLD_MAX_ITS", cs%max_MLD_its, &
2150 "The maximum number of iterations that can be used to find a self-consistent "//&
2151 "mixed layer depth. For now, due to the use of bisection, the maximum number "//&
2152 "iteractions needed is set by Depth/2^MAX_ITS < EPBL_MLD_TOLERANCE.", &
2153 default=20, do_not_log=.not.cs%Use_MLD_iteration)
2154 if (.not.cs%Use_MLD_iteration) cs%Max_MLD_Its = 1
2155 call get_param(param_file, mdl,
"EPBL_MIN_MIX_LEN", cs%min_mix_len, &
2156 "The minimum mixing length scale that will be used "//&
2157 "by ePBL. The default (0) does not set a minimum.", &
2158 units=
"meter", default=0.0, scale=us%m_to_Z)
2160 call get_param(param_file, mdl,
"MIX_LEN_EXPONENT", cs%MixLenExponent, &
2161 "The exponent applied to the ratio of the distance to the MLD "//&
2162 "and the MLD depth which determines the shape of the mixing length. "//&
2163 "This is only used if USE_MLD_ITERATION is True.", &
2164 units=
"nondim", default=2.0)
2167 call get_param(param_file, mdl,
"EPBL_VEL_SCALE_SCHEME", tmpstr, &
2168 "Selects the method for translating TKE into turbulent velocities. "//&
2169 "Valid values are: \n"//&
2170 "\t CUBE_ROOT_TKE - A constant times the cube root of remaining TKE. \n"//&
2171 "\t REICHL_H18 - Use the scheme based on a combination of w* and v* as \n"//&
2172 "\t documented in Reichl & Hallberg, 2018.", &
2174 call get_param(param_file, mdl,
"EPBL_VEL_SCALE_MODE", wt_mode, default=-1)
2175 if (wt_mode == 0)
then
2177 call mom_error(warning,
"Use EPBL_VEL_SCALE_SCHEME = CUBE_ROOT_TKE instead of the archaic EPBL_VEL_SCALE_MODE = 0.")
2178 elseif (wt_mode == 1)
then
2180 call mom_error(warning,
"Use EPBL_VEL_SCALE_SCHEME = REICHL_H18 instead of the archaic EPBL_VEL_SCALE_MODE = 1.")
2181 elseif (wt_mode >= 2)
then
2182 call mom_error(fatal,
"An unrecognized value of the obsolete parameter EPBL_VEL_SCALE_MODE was specified.")
2184 call log_param(param_file, mdl,
"EPBL_VEL_SCALE_SCHEME", tmpstr, &
2185 "Selects the method for translating TKE into turbulent velocities. "//&
2186 "Valid values are: \n"//&
2187 "\t CUBE_ROOT_TKE - A constant times the cube root of remaining TKE. \n"//&
2188 "\t REICHL_H18 - Use the scheme based on a combination of w* and v* as \n"//&
2189 "\t documented in Reichl & Hallberg, 2018.", &
2192 select case (tmpstr)
2198 call mom_mesg(
'energetic_PBL_init: EPBL_VEL_SCALE_SCHEME ="'//trim(tmpstr)//
'"', 0)
2199 call mom_error(fatal,
"energetic_PBL_init: Unrecognized setting "// &
2200 "EPBL_VEL_SCALE_SCHEME = "//trim(tmpstr)//
" found in input file.")
2203 call get_param(param_file, mdl,
"WSTAR_USTAR_COEF", cs%wstar_ustar_coef, &
2204 "A ratio relating the efficiency with which convectively "//&
2205 "released energy is converted to a turbulent velocity, "//&
2206 "relative to mechanically forced TKE. Making this larger "//&
2207 "increases the BL diffusivity", units=
"nondim", default=1.0)
2208 call get_param(param_file, mdl,
"EPBL_VEL_SCALE_FACTOR", cs%vstar_scale_fac, &
2209 "An overall nondimensional scaling factor for wT. "//&
2210 "Making this larger increases the PBL diffusivity.", &
2211 units=
"nondim", default=1.0)
2212 call get_param(param_file, mdl,
"VSTAR_SURF_FAC", cs%vstar_surf_fac,&
2213 "The proportionality times ustar to set vstar at the surface.", &
2214 units=
"nondim", default=1.2)
2217 call get_param(param_file, mdl,
"USE_LA_LI2016", use_la_windsea, &
2218 "A logical to use the Li et al. 2016 (submitted) formula to "//&
2219 "determine the Langmuir number.", units=
"nondim", default=.false.)
2221 if (use_la_windsea)
then
2224 call get_param(param_file, mdl,
"EPBL_LT", cs%USE_LT, &
2225 "A logical to use a LT parameterization.", &
2226 units=
"nondim", default=.false.)
2229 call get_param(param_file, mdl,
"EPBL_LANGMUIR_SCHEME", tmpstr, &
2230 "EPBL_LANGMUIR_SCHEME selects the method for including Langmuir turbulence. "//&
2231 "Valid values are: \n"//&
2232 "\t NONE - Do not do any extra mixing due to Langmuir turbulence \n"//&
2233 "\t RESCALE - Use a multiplicative rescaling of mstar to account for Langmuir turbulence \n"//&
2234 "\t ADDITIVE - Add a Langmuir turblence contribution to mstar to other contributions", &
2236 call get_param(param_file, mdl,
"LT_ENHANCE", lt_enhance, default=-1)
2237 if (lt_enhance == 0)
then
2239 call mom_error(warning,
"Use EPBL_LANGMUIR_SCHEME = NONE instead of the archaic LT_ENHANCE = 0.")
2240 elseif (lt_enhance == 1)
then
2241 call mom_error(fatal,
"You are using a legacy LT_ENHANCE mode in ePBL that has been phased out. "//&
2242 "If you need to use this setting please report this error. Also use "//&
2243 "EPBL_LANGMUIR_SCHEME to specify the scheme for mstar.")
2244 elseif (lt_enhance == 2)
then
2246 call mom_error(warning,
"Use EPBL_LANGMUIR_SCHEME = RESCALE instead of the archaic LT_ENHANCE = 2.")
2247 elseif (lt_enhance == 3)
then
2249 call mom_error(warning,
"Use EPBL_LANGMUIR_SCHEME = ADDITIVE instead of the archaic LT_ENHANCE = 3.")
2250 elseif (lt_enhance > 3)
then
2251 call mom_error(fatal,
"An unrecognized value of the obsolete parameter LT_ENHANCE was specified.")
2253 call log_param(param_file, mdl,
"EPBL_LANGMUIR_SCHEME", tmpstr, &
2254 "EPBL_LANGMUIR_SCHEME selects the method for including Langmuir turbulence. "//&
2255 "Valid values are: \n"//&
2256 "\t NONE - Do not do any extra mixing due to Langmuir turbulence \n"//&
2257 "\t RESCALE - Use a multiplicative rescaling of mstar to account for Langmuir turbulence \n"//&
2258 "\t ADDITIVE - Add a Langmuir turblence contribution to mstar to other contributions", &
2261 select case (tmpstr)
2269 call mom_mesg(
'energetic_PBL_init: EPBL_LANGMUIR_SCHEME ="'//trim(tmpstr)//
'"', 0)
2270 call mom_error(fatal,
"energetic_PBL_init: Unrecognized setting "// &
2271 "EPBL_LANGMUIR_SCHEME = "//trim(tmpstr)//
" found in input file.")
2274 call get_param(param_file, mdl,
"LT_ENHANCE_COEF", cs%LT_ENHANCE_COEF, &
2275 "Coefficient for Langmuir enhancement of mstar", &
2276 units=
"nondim", default=0.447, do_not_log=(cs%LT_enhance_form==
no_langmuir))
2277 call get_param(param_file, mdl,
"LT_ENHANCE_EXP", cs%LT_ENHANCE_EXP, &
2278 "Exponent for Langmuir enhancementt of mstar", &
2279 units=
"nondim", default=-1.33, do_not_log=(cs%LT_enhance_form==
no_langmuir))
2280 call get_param(param_file, mdl,
"LT_MOD_LAC1", cs%LaC_MLDoEK, &
2281 "Coefficient for modification of Langmuir number due to "//&
2282 "MLD approaching Ekman depth.", &
2283 units=
"nondim", default=-0.87, do_not_log=(cs%LT_enhance_form==
no_langmuir))
2284 call get_param(param_file, mdl,
"LT_MOD_LAC2", cs%LaC_MLDoOB_stab, &
2285 "Coefficient for modification of Langmuir number due to "//&
2286 "MLD approaching stable Obukhov depth.", &
2287 units=
"nondim", default=0.0, do_not_log=(cs%LT_enhance_form==
no_langmuir))
2288 call get_param(param_file, mdl,
"LT_MOD_LAC3", cs%LaC_MLDoOB_un, &
2289 "Coefficient for modification of Langmuir number due to "//&
2290 "MLD approaching unstable Obukhov depth.", &
2291 units=
"nondim", default=0.0, do_not_log=(cs%LT_enhance_form==
no_langmuir))
2292 call get_param(param_file, mdl,
"LT_MOD_LAC4", cs%Lac_EKoOB_stab, &
2293 "Coefficient for modification of Langmuir number due to "//&
2294 "ratio of Ekman to stable Obukhov depth.", &
2295 units=
"nondim", default=0.95, do_not_log=(cs%LT_enhance_form==
no_langmuir))
2296 call get_param(param_file, mdl,
"LT_MOD_LAC5", cs%Lac_EKoOB_un, &
2297 "Coefficient for modification of Langmuir number due to "//&
2298 "ratio of Ekman to unstable Obukhov depth.", &
2299 units=
"nondim", default=0.95, do_not_log=(cs%LT_enhance_form==
no_langmuir))
2305 cs%ustar_min = 2e-4*cs%omega*(gv%Angstrom_Z + gv%H_to_Z*gv%H_subroundoff)
2306 call log_param(param_file, mdl,
"EPBL_USTAR_MIN", cs%ustar_min*us%Z_to_m*us%s_to_T, &
2307 "The (tiny) minimum friction velocity used within the "//&
2308 "ePBL code, derived from OMEGA and ANGSTROM.", units=
"m s-1")
2312 r_z3_t3_to_kg_s3 = us%R_to_kg_m3 * us%Z_to_m**3 * us%s_to_T**3
2313 cs%id_ML_depth = register_diag_field(
'ocean_model',
'ePBL_h_ML', diag%axesT1, &
2314 time,
'Surface boundary layer depth',
'm', conversion=us%Z_to_m, &
2315 cmor_long_name=
'Ocean Mixed Layer Thickness Defined by Mixing Scheme')
2316 cs%id_TKE_wind = register_diag_field(
'ocean_model',
'ePBL_TKE_wind', diag%axesT1, &
2317 time,
'Wind-stirring source of mixed layer TKE',
'm3 s-3', conversion=r_z3_t3_to_kg_s3)
2318 cs%id_TKE_MKE = register_diag_field(
'ocean_model',
'ePBL_TKE_MKE', diag%axesT1, &
2319 time,
'Mean kinetic energy source of mixed layer TKE',
'm3 s-3', conversion=r_z3_t3_to_kg_s3)
2320 cs%id_TKE_conv = register_diag_field(
'ocean_model',
'ePBL_TKE_conv', diag%axesT1, &
2321 time,
'Convective source of mixed layer TKE',
'm3 s-3', conversion=r_z3_t3_to_kg_s3)
2322 cs%id_TKE_forcing = register_diag_field(
'ocean_model',
'ePBL_TKE_forcing', diag%axesT1, &
2323 time,
'TKE consumed by mixing surface forcing or penetrative shortwave radation '//&
2324 'through model layers',
'm3 s-3', conversion=r_z3_t3_to_kg_s3)
2325 cs%id_TKE_mixing = register_diag_field(
'ocean_model',
'ePBL_TKE_mixing', diag%axesT1, &
2326 time,
'TKE consumed by mixing that deepens the mixed layer',
'm3 s-3', conversion=r_z3_t3_to_kg_s3)
2327 cs%id_TKE_mech_decay = register_diag_field(
'ocean_model',
'ePBL_TKE_mech_decay', diag%axesT1, &
2328 time,
'Mechanical energy decay sink of mixed layer TKE',
'm3 s-3', conversion=r_z3_t3_to_kg_s3)
2329 cs%id_TKE_conv_decay = register_diag_field(
'ocean_model',
'ePBL_TKE_conv_decay', diag%axesT1, &
2330 time,
'Convective energy decay sink of mixed layer TKE',
'm3 s-3', conversion=r_z3_t3_to_kg_s3)
2331 cs%id_Mixing_Length = register_diag_field(
'ocean_model',
'Mixing_Length', diag%axesTi, &
2332 time,
'Mixing Length that is used',
'm', conversion=us%Z_to_m)
2333 cs%id_Velocity_Scale = register_diag_field(
'ocean_model',
'Velocity_Scale', diag%axesTi, &
2334 time,
'Velocity Scale that is used.',
'm s-1', conversion=us%Z_to_m*us%s_to_T)
2335 cs%id_MSTAR_mix = register_diag_field(
'ocean_model',
'MSTAR', diag%axesT1, &
2336 time,
'Total mstar that is used.',
'nondim')
2339 cs%id_LA = register_diag_field(
'ocean_model',
'LA', diag%axesT1, &
2340 time,
'Langmuir number.',
'nondim')
2341 cs%id_LA_mod = register_diag_field(
'ocean_model',
'LA_MOD', diag%axesT1, &
2342 time,
'Modified Langmuir number.',
'nondim')
2343 cs%id_MSTAR_LT = register_diag_field(
'ocean_model',
'MSTAR_LT', diag%axesT1, &
2344 time,
'Increase in mstar due to Langmuir Turbulence.',
'nondim')
2347 call get_param(param_file, mdl,
"ENABLE_THERMODYNAMICS", use_temperature, &
2348 "If true, temperature and salinity are used as state "//&
2349 "variables.", default=.true.)
2351 if (max(cs%id_TKE_wind, cs%id_TKE_MKE, cs%id_TKE_conv, &
2352 cs%id_TKE_mixing, cs%id_TKE_mech_decay, cs%id_TKE_forcing, &
2353 cs%id_TKE_conv_decay) > 0)
then
2354 call safe_alloc_alloc(cs%diag_TKE_wind, isd, ied, jsd, jed)
2355 call safe_alloc_alloc(cs%diag_TKE_MKE, isd, ied, jsd, jed)
2356 call safe_alloc_alloc(cs%diag_TKE_conv, isd, ied, jsd, jed)
2357 call safe_alloc_alloc(cs%diag_TKE_forcing, isd, ied, jsd, jed)
2358 call safe_alloc_alloc(cs%diag_TKE_mixing, isd, ied, jsd, jed)
2359 call safe_alloc_alloc(cs%diag_TKE_mech_decay, isd, ied, jsd, jed)
2360 call safe_alloc_alloc(cs%diag_TKE_conv_decay, isd, ied, jsd, jed)
2362 cs%TKE_diagnostics = .true.
2364 if (cs%id_Velocity_Scale>0)
call safe_alloc_alloc(cs%Velocity_Scale, isd, ied, jsd, jed, gv%ke+1)
2365 if (cs%id_Mixing_Length>0)
call safe_alloc_alloc(cs%Mixing_Length, isd, ied, jsd, jed, gv%ke+1)
2367 call safe_alloc_alloc(cs%ML_depth, isd, ied, jsd, jed)
2368 if (max(cs%id_mstar_mix, cs%id_LA, cs%id_LA_mod, cs%id_MSTAR_LT ) >0)
then
2369 call safe_alloc_alloc(cs%Mstar_mix, isd, ied, jsd, jed)
2370 call safe_alloc_alloc(cs%LA, isd, ied, jsd, jed)
2371 call safe_alloc_alloc(cs%LA_MOD, isd, ied, jsd, jed)
2372 call safe_alloc_alloc(cs%MSTAR_LT, isd, ied, jsd, jed)
2382 if (.not.
associated(cs))
return
2384 if (
allocated(cs%ML_depth))
deallocate(cs%ML_depth)
2385 if (
allocated(cs%LA))
deallocate(cs%LA)
2386 if (
allocated(cs%LA_MOD))
deallocate(cs%LA_MOD)
2387 if (
allocated(cs%MSTAR_MIX))
deallocate(cs%MSTAR_MIX)
2388 if (
allocated(cs%MSTAR_LT))
deallocate(cs%MSTAR_LT)
2389 if (
allocated(cs%diag_TKE_wind))
deallocate(cs%diag_TKE_wind)
2390 if (
allocated(cs%diag_TKE_MKE))
deallocate(cs%diag_TKE_MKE)
2391 if (
allocated(cs%diag_TKE_conv))
deallocate(cs%diag_TKE_conv)
2392 if (
allocated(cs%diag_TKE_forcing))
deallocate(cs%diag_TKE_forcing)
2393 if (
allocated(cs%diag_TKE_mixing))
deallocate(cs%diag_TKE_mixing)
2394 if (
allocated(cs%diag_TKE_mech_decay))
deallocate(cs%diag_TKE_mech_decay)
2395 if (
allocated(cs%diag_TKE_conv_decay))
deallocate(cs%diag_TKE_conv_decay)
2396 if (
allocated(cs%Mixing_Length))
deallocate(cs%Mixing_Length)
2397 if (
allocated(cs%Velocity_Scale))
deallocate(cs%Velocity_Scale)