Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

+Non-Boussinesq revisions to wave_interface #432

Merged
merged 2 commits into from
Aug 6, 2023
Merged

+Non-Boussinesq revisions to wave_interface #432

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
39fd048
+Non-Boussinesq revisions to wave_interface
Hallberg-NOAA Aug 3, 2023
27abbfa
Merge branch 'dev/gfdl' into nonBous_wave_interface
marshallward Aug 6, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
26 changes: 20 additions & 6 deletions src/core/MOM.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,7 @@ module MOM
use MOM_error_handler, only : MOM_set_verbosity, callTree_showQuery
use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint
use MOM_file_parser, only : read_param, get_param, log_version, param_file_type
use MOM_forcing_type, only : forcing, mech_forcing
use MOM_forcing_type, only : forcing, mech_forcing, find_ustar
use MOM_forcing_type, only : MOM_forcing_chksum, MOM_mech_forcing_chksum
use MOM_get_input, only : Get_MOM_Input, directories
use MOM_io, only : MOM_io_init, vardesc, var_desc
Expand Down Expand Up @@ -91,7 +91,7 @@ module MOM
use MOM_grid, only : set_first_direction, rescale_grid_bathymetry
use MOM_hor_index, only : hor_index_type, hor_index_init
use MOM_hor_index, only : rotate_hor_index
use MOM_interface_heights, only : find_eta, calc_derived_thermo
use MOM_interface_heights, only : find_eta, calc_derived_thermo, thickness_to_dz
use MOM_interface_filter, only : interface_filter, interface_filter_init, interface_filter_end
use MOM_interface_filter, only : interface_filter_CS
use MOM_lateral_mixing_coeffs, only : calc_slope_functions, VarMix_init, VarMix_end
Expand Down Expand Up @@ -544,8 +544,13 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS
! the end of a stepping cycle (whatever that may mean).
logical :: therm_reset ! If true, reset running sums of thermodynamic quantities.
real :: cycle_time ! The length of the coupled time-stepping cycle [T ~> s].
real, dimension(SZI_(CS%G),SZJ_(CS%G)) :: &
U_star ! The wind friction velocity, calculated using the Boussinesq reference density or
! the time-evolving surface density in non-Boussinesq mode [Z T-1 ~> m s-1]
real, dimension(SZI_(CS%G),SZJ_(CS%G)) :: &
ssh ! sea surface height, which may be based on eta_av [Z ~> m]
real, dimension(SZI_(CS%G),SZJ_(CS%G),SZK_(CS%GV)) :: &
dz ! Vertical distance across layers [Z ~> m]

real, dimension(:,:,:), pointer :: &
u => NULL(), & ! u : zonal velocity component [L T-1 ~> m s-1]
Expand Down Expand Up @@ -672,13 +677,18 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS

dt = time_interval / real(n_max)
dt_therm = dt ; ntstep = 1

if (CS%UseWaves .and. associated(fluxes%ustar)) &
call pass_var(fluxes%ustar, G%Domain, clock=id_clock_pass, halo=1)
if (CS%UseWaves .and. associated(fluxes%tau_mag)) &
call pass_var(fluxes%tau_mag, G%Domain, clock=id_clock_pass, halo=1)

if (associated(fluxes%p_surf)) p_surf => fluxes%p_surf
CS%tv%p_surf => NULL()
if (CS%use_p_surf_in_EOS .and. associated(fluxes%p_surf)) then
CS%tv%p_surf => fluxes%p_surf
if (allocated(CS%tv%SpV_avg)) call pass_var(fluxes%p_surf, G%Domain, clock=id_clock_pass)
endif
if (CS%UseWaves) call pass_var(fluxes%ustar, G%Domain, clock=id_clock_pass)
endif

if (therm_reset) then
Expand Down Expand Up @@ -722,12 +732,16 @@ subroutine step_MOM(forces_in, fluxes_in, sfc_state, Time_start, time_int_in, CS
if (CS%UseWaves) then
! Update wave information, which is presently kept static over each call to step_mom
call enable_averages(time_interval, Time_start + real_to_time(US%T_to_s*time_interval), CS%diag)
call Update_Stokes_Drift(G, GV, US, Waves, h, forces%ustar, time_interval, do_dyn)
call find_ustar(forces, CS%tv, U_star, G, GV, US, halo=1)
call thickness_to_dz(h, CS%tv, dz, G, GV, US, halo_size=1)
call Update_Stokes_Drift(G, GV, US, Waves, dz, U_star, time_interval, do_dyn)
call disable_averaging(CS%diag)
endif
else ! not do_dyn.
if (CS%UseWaves) then ! Diagnostics are not enabled in this call.
call Update_Stokes_Drift(G, GV, US, Waves, h, fluxes%ustar, time_interval, do_dyn)
call find_ustar(fluxes, CS%tv, U_star, G, GV, US, halo=1)
call thickness_to_dz(h, CS%tv, dz, G, GV, US, halo_size=1)
call Update_Stokes_Drift(G, GV, US, Waves, dz, U_star, time_interval, do_dyn)
endif
endif

Expand Down Expand Up @@ -3261,7 +3275,7 @@ subroutine finish_MOM_initialization(Time, dirs, CS)
! Write initial conditions
if (CS%write_IC) then
allocate(restart_CSp_tmp)
restart_CSP_tmp = CS%restart_CS
restart_CSp_tmp = CS%restart_CS
call restart_registry_lock(restart_CSp_tmp, unlocked=.true.)
allocate(z_interface(SZI_(G),SZJ_(G),SZK_(GV)+1))
call find_eta(CS%h, CS%tv, G, GV, US, z_interface, dZref=G%Z_ref)
Expand Down
6 changes: 4 additions & 2 deletions src/core/MOM_dynamics_split_RK2.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -481,7 +481,8 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
Use_Stokes_PGF = associated(Waves)
if (Use_Stokes_PGF) Use_Stokes_PGF = Waves%Stokes_PGF
if (Use_Stokes_PGF) then
call Stokes_PGF(G, GV, h, u, v, CS%PFu_Stokes, CS%PFv_Stokes, Waves)
call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1)
call Stokes_PGF(G, GV, US, dz, u, v, CS%PFu_Stokes, CS%PFv_Stokes, Waves)

! We are adding Stokes_PGF to hydrostatic PGF here. The diag PFu/PFv
! will therefore report the sum total PGF and we avoid other
Expand Down Expand Up @@ -748,7 +749,8 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
Use_Stokes_PGF = associated(Waves)
if (Use_Stokes_PGF) Use_Stokes_PGF = Waves%Stokes_PGF
if (Use_Stokes_PGF) then
call Stokes_PGF(G, GV, h, u, v, CS%PFu_Stokes, CS%PFv_Stokes, Waves)
call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1)
call Stokes_PGF(G, GV, US, dz, u, v, CS%PFu_Stokes, CS%PFv_Stokes, Waves)
if (.not.Waves%Passive_Stokes_PGF) then
do k=1,nz
do j=js,je ; do I=Isq,Ieq
Expand Down
19 changes: 12 additions & 7 deletions src/parameterizations/vertical/MOM_CVMix_KPP.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@ module MOM_CVMix_KPP
use MOM_file_parser, only : get_param, log_param, log_version, param_file_type
use MOM_file_parser, only : openParameterBlock, closeParameterBlock
use MOM_grid, only : ocean_grid_type, isPointInCell
use MOM_interface_heights, only : thickness_to_dz
use MOM_unit_scaling, only : unit_scale_type
use MOM_variables, only : thermo_var_ptrs
use MOM_verticalGrid, only : verticalGrid_type
Expand Down Expand Up @@ -604,7 +605,7 @@ subroutine KPP_calculate(CS, G, GV, US, h, tv, uStar, buoyFlux, Kt, Ks, Kv, &
type(ocean_grid_type), intent(in) :: G !< Ocean grid
type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer/level thicknesses [H ~> m or kg m-2]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure.
real, dimension(SZI_(G),SZJ_(G)), intent(in) :: uStar !< Surface friction velocity [Z T-1 ~> m s-1]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: buoyFlux !< Surface buoyancy flux [L2 T-3 ~> m2 s-3]
Expand Down Expand Up @@ -863,14 +864,14 @@ subroutine KPP_calculate(CS, G, GV, US, h, tv, uStar, buoyFlux, Kt, Ks, Kv, &
Kt(i,j,k) = Kt(i,j,k) + GV%m2_s_to_HZ_T * Kdiffusivity(k,1)
Ks(i,j,k) = Ks(i,j,k) + GV%m2_s_to_HZ_T * Kdiffusivity(k,2)
Kv(i,j,k) = Kv(i,j,k) + GV%m2_s_to_HZ_T * Kviscosity(k)
if (CS%Stokes_Mixing) Waves%KvS(i,j,k) = GV%H_to_Z*Kv(i,j,k)
if (CS%Stokes_Mixing) Waves%KvS(i,j,k) = Kv(i,j,k)
enddo
else ! KPP replaces prior diffusivity when former is non-zero
do k=1, GV%ke+1
if (Kdiffusivity(k,1) /= 0.) Kt(i,j,k) = GV%m2_s_to_HZ_T * Kdiffusivity(k,1)
if (Kdiffusivity(k,2) /= 0.) Ks(i,j,k) = GV%m2_s_to_HZ_T * Kdiffusivity(k,2)
if (Kviscosity(k) /= 0.) Kv(i,j,k) = GV%m2_s_to_HZ_T * Kviscosity(k)
if (CS%Stokes_Mixing) Waves%KvS(i,j,k) = GV%H_to_Z*Kv(i,j,k)
if (CS%Stokes_Mixing) Waves%KvS(i,j,k) = Kv(i,j,k)
enddo
endif
endif
Expand Down Expand Up @@ -912,7 +913,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl
type(ocean_grid_type), intent(inout) :: G !< Ocean grid
type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer/level thicknesses [H ~> m or kg m-2]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: Temp !< potential/cons temp [C ~> degC]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: Salt !< Salinity [S ~> ppt]
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: u !< Velocity i-component [L T-1 ~> m s-1]
Expand All @@ -924,6 +925,8 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl
real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: lamult !< Langmuir enhancement factor [nondim]

! Local variables
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: dz ! Height change across layers [Z ~> m]

! Variables for passing to CVMix routines, often in MKS units
real, dimension( GV%ke ) :: Ws_1d ! Profile of vertical velocity scale for scalars in MKS units [m s-1]
real, dimension( GV%ke ) :: deltaRho ! delta Rho in numerator of Bulk Ri number [R ~> kg m-3]
Expand All @@ -940,7 +943,6 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl
real :: Coriolis ! Coriolis parameter at tracer points in MKS units [s-1]
real :: KPP_OBL_depth ! Boundary layer depth calculated by CVMix_kpp_compute_OBL_depth in MKS units [m]


! Variables for EOS calculations
real, dimension( 3*GV%ke ) :: rho_1D ! A column of densities [R ~> kg m-3]
real, dimension( 3*GV%ke ) :: pres_1D ! A column of pressures [R L2 T-2 ~> Pa]
Expand Down Expand Up @@ -996,6 +998,9 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl
GoRho = US%Z_to_m*US%s_to_T**2 * GoRho_Z_L2
buoy_scale = US%L_to_m**2*US%s_to_T**3

! Find the vertical distances across layers.
call thickness_to_dz(h, tv, dz, G, GV, US)

! loop over horizontal points on processor
!$OMP parallel do default(none) private(surfFricVel, iFaceHeight, hcorr, dh, cellHeight, &
!$OMP surfBuoyFlux, U_H, V_H, Coriolis, pRef, SLdepth_0d, vt2_1d, &
Expand All @@ -1005,7 +1010,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl
!$OMP Temp_1D, salt_1D, surfBuoyFlux2, MLD_guess, LA, rho_1D, &
!$OMP deltarho, N2_1d, ws_1d, LangEnhVT2,KPP_OBL_depth, z_cell, &
!$OMP z_inter, OBL_depth, BulkRi_1d, zBottomMinusOffset) &
!$OMP shared(G, GV, CS, US, uStar, h, buoy_scale, buoyFlux, &
!$OMP shared(G, GV, CS, US, uStar, h, dz, buoy_scale, buoyFlux, &
!$OMP Temp, Salt, waves, tv, GoRho, GoRho_Z_L2, u, v, lamult)
do j = G%jsc, G%jec
do i = G%isc, G%iec ; if (G%mask2dT(i,j) > 0.0) then
Expand Down Expand Up @@ -1127,7 +1132,7 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl
if ( (CS%LT_K_ENHANCEMENT .or. CS%LT_VT2_ENHANCEMENT) .and. .not. present(lamult)) then
MLD_guess = max( CS%MLD_guess_min, abs(CS%OBLdepthprev(i,j) ) )
call get_Langmuir_Number(LA, G, GV, US, MLD_guess, uStar(i,j), i, j, &
H=H(i,j,:), U_H=U_H, V_H=V_H, WAVES=WAVES)
dz=dz(i,j,:), U_H=U_H, V_H=V_H, WAVES=WAVES)
CS%La_SL(i,j) = LA
endif

Expand Down
26 changes: 16 additions & 10 deletions src/parameterizations/vertical/MOM_energetic_PBL.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@ module MOM_energetic_PBL
use MOM_file_parser, only : get_param, log_param, log_version, param_file_type
use MOM_forcing_type, only : forcing
use MOM_grid, only : ocean_grid_type
use MOM_interface_heights, only : thickness_to_dz
use MOM_string_functions, only : uppercase
use MOM_unit_scaling, only : unit_scale_type
use MOM_variables, only : thermo_var_ptrs
Expand Down Expand Up @@ -309,6 +310,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS
! Local variables
real, dimension(SZI_(G),SZK_(GV)) :: &
h_2d, & ! A 2-d slice of the layer thickness [H ~> m or kg m-2].
dz_2d, & ! A 2-d slice of the vertical distance across layers [Z ~> m].
T_2d, & ! A 2-d slice of the layer temperatures [C ~> degC].
S_2d, & ! A 2-d slice of the layer salinities [S ~> ppt].
TKE_forced_2d, & ! A 2-d slice of TKE_forced [R Z3 T-2 ~> J m-2].
Expand All @@ -320,6 +322,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS
Kd_2d ! A 2-d version of the diapycnal diffusivity [Z2 T-1 ~> m2 s-1].
real, dimension(SZK_(GV)) :: &
h, & ! The layer thickness [H ~> m or kg m-2].
dz, & ! The vertical distance across layers [Z ~> m].
T0, & ! The initial layer temperatures [C ~> degC].
S0, & ! The initial layer salinities [S ~> ppt].
dSV_dT_1d, & ! The partial derivatives of specific volume with temperature [R-1 C-1 ~> m3 kg-1 degC-1].
Expand Down Expand Up @@ -362,7 +365,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS

! Zero out diagnostics before accumulation.
if (CS%TKE_diagnostics) then
!!OMP parallel do default(none) shared(is,ie,js,je,CS)
!!OMP parallel do default(none) shared(is,ie,js,je,CS)
do j=js,je ; do i=is,ie
CS%diag_TKE_wind(i,j) = 0.0 ; CS%diag_TKE_MKE(i,j) = 0.0
CS%diag_TKE_conv(i,j) = 0.0 ; CS%diag_TKE_forcing(i,j) = 0.0
Expand All @@ -373,8 +376,8 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS
! if (CS%id_Mixing_Length>0) CS%Mixing_Length(:,:,:) = 0.0
! if (CS%id_Velocity_Scale>0) CS%Velocity_Scale(:,:,:) = 0.0

!!OMP parallel do default(private) shared(js,je,nz,is,ie,h_3d,u_3d,v_3d,tv,dt, &
!!OMP CS,G,GV,US,fluxes,TKE_forced,dSV_dT,dSV_dS,Kd_int)
!!OMP parallel do default(private) shared(js,je,nz,is,ie,h_3d,u_3d,v_3d,tv,dt, &
!!OMP CS,G,GV,US,fluxes,TKE_forced,dSV_dT,dSV_dS,Kd_int)
do j=js,je
! Copy the thicknesses and other fields to 2-d arrays.
do k=1,nz ; do i=is,ie
Expand All @@ -383,6 +386,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS
TKE_forced_2d(i,k) = TKE_forced(i,j,k)
dSV_dT_2d(i,k) = dSV_dT(i,j,k) ; dSV_dS_2d(i,k) = dSV_dS(i,j,k)
enddo ; enddo
call thickness_to_dz(h_3d, tv, dz_2d, j, G, GV)

! Determine the initial mech_TKE and conv_PErel, including the energy required
! to mix surface heating through the topmost cell, the energy released by mixing
Expand All @@ -394,7 +398,8 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS

! Copy the thicknesses and other fields to 1-d arrays.
do k=1,nz
h(k) = h_2d(i,k) + GV%H_subroundoff ; u(k) = u_2d(i,k) ; v(k) = v_2d(i,k)
h(k) = h_2d(i,k) + GV%H_subroundoff ; dz(k) = dz_2d(i,k) + GV%dZ_subroundoff
u(k) = u_2d(i,k) ; v(k) = v_2d(i,k)
T0(k) = T_2d(i,k) ; S0(k) = S_2d(i,k) ; TKE_forcing(k) = TKE_forced_2d(i,k)
dSV_dT_1d(k) = dSV_dT_2d(i,k) ; dSV_dS_1d(k) = dSV_dS_2d(i,k)
enddo
Expand All @@ -421,15 +426,15 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS

! Perhaps provide a first guess for MLD based on a stored previous value.
MLD_io = -1.0
if (CS%MLD_iteration_guess .and. (CS%ML_Depth(i,j) > 0.0)) MLD_io = CS%ML_Depth(i,j)
if (CS%MLD_iteration_guess .and. (CS%ML_depth(i,j) > 0.0)) MLD_io = CS%ML_depth(i,j)

if (stoch_CS%pert_epbl) then ! stochastics are active
call ePBL_column(h, u, v, T0, S0, dSV_dT_1d, dSV_dS_1d, TKE_forcing, B_flux, absf, &
call ePBL_column(h, dz, u, v, T0, S0, dSV_dT_1d, dSV_dS_1d, TKE_forcing, B_flux, absf, &
u_star, u_star_mean, dt, MLD_io, Kd, mixvel, mixlen, GV, &
US, CS, eCD, Waves, G, i, j, &
TKE_gen_stoch=stoch_CS%epbl1_wts(i,j), TKE_diss_stoch=stoch_CS%epbl2_wts(i,j))
else
call ePBL_column(h, u, v, T0, S0, dSV_dT_1d, dSV_dS_1d, TKE_forcing, B_flux, absf, &
call ePBL_column(h, dz, u, v, T0, S0, dSV_dT_1d, dSV_dS_1d, TKE_forcing, B_flux, absf, &
u_star, u_star_mean, dt, MLD_io, Kd, mixvel, mixlen, GV, &
US, CS, eCD, Waves, G, i, j)
endif
Expand Down Expand Up @@ -499,12 +504,13 @@ end subroutine energetic_PBL

!> This subroutine determines the diffusivities from the integrated energetics
!! mixed layer model for a single column of water.
subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, absf, &
subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, absf, &
u_star, u_star_mean, dt, MLD_io, Kd, mixvel, mixlen, GV, US, CS, eCD, &
Waves, G, i, j, TKE_gen_stoch, TKE_diss_stoch)
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
real, dimension(SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
real, dimension(SZK_(GV)), intent(in) :: dz !< The vertical distance across layers [Z ~> m].
real, dimension(SZK_(GV)), intent(in) :: u !< Zonal velocities interpolated to h points
!! [L T-1 ~> m s-1].
real, dimension(SZK_(GV)), intent(in) :: v !< Zonal velocities interpolated to h points
Expand Down Expand Up @@ -828,7 +834,7 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs
!/ Here we get MStar, which is the ratio of convective TKE driven mixing to UStar**3
MLD_guess_z = GV%H_to_Z*MLD_guess ! Convert MLD from thickness to height coordinates for these calls
if (CS%Use_LT) then
call get_Langmuir_Number(LA, G, GV, US, abs(MLD_guess_z), u_star_mean, i, j, h, Waves, &
call get_Langmuir_Number(LA, G, GV, US, abs(MLD_guess_z), u_star_mean, i, j, dz, Waves, &
U_H=u, V_H=v)
call find_mstar(CS, US, B_flux, u_star, u_star_Mean, MLD_guess_z, absf, &
MStar_total, Langmuir_Number=La, Convect_Langmuir_Number=LAmod,&
Expand Down Expand Up @@ -1931,7 +1937,7 @@ subroutine energetic_PBL_get_MLD(CS, MLD, G, US, m_to_MLD_units)
scale = 1.0 ; if (present(m_to_MLD_units)) scale = US%Z_to_m * m_to_MLD_units

do j=G%jsc,G%jec ; do i=G%isc,G%iec
MLD(i,j) = scale*CS%ML_Depth(i,j)
MLD(i,j) = scale*CS%ML_depth(i,j)
enddo ; enddo

end subroutine energetic_PBL_get_MLD
Expand Down
Loading