From 245bf6faef4970e829e9c7e18bd4cf9786f8590e Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Tue, 18 Jul 2023 05:27:01 -0400 Subject: [PATCH] +*Add and use find_ustar Added the new public interface find_ustar to extract the friction velocity from either a forcing type argument, or a mech_forcing_type argument, either directly or from tau_mag, and in non-Boussinesq mode by using the time-evolving surface specific volume. Find_ustar is an overloaded interface to find_ustar_fluxes or find_ustar_mech_forcing, which are the same but for the type of one of their arguments. For now, the subroutines bulkmixedlayer, mixedlayer_restrajt_OM4, mixedlayer_restrat_Bodner and mixedlayer_restrat_BML are calling find_ustar to avoid code duplication during the transition to work in fully non-Boussinesq mode, but it will eventually be used in about another half dozen other places. All Boussinesq answers are bitwise identical, but non-Boussinesq answers will change and become less dependent on the Boussinesq reference density, and there is a new publicly visible interface wrapping two subroutines. --- src/core/MOM_forcing_type.F90 | 141 +++++++++++++++++- .../lateral/MOM_mixed_layer_restrat.F90 | 36 +++-- .../vertical/MOM_bulk_mixed_layer.F90 | 19 ++- 3 files changed, 182 insertions(+), 14 deletions(-) diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index a59c33d525..bdedccb1bd 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -28,7 +28,7 @@ module MOM_forcing_type public extractFluxes1d, extractFluxes2d, optics_type public MOM_forcing_chksum, MOM_mech_forcing_chksum -public calculateBuoyancyFlux1d, calculateBuoyancyFlux2d +public calculateBuoyancyFlux1d, calculateBuoyancyFlux2d, find_ustar public forcing_accumulate, fluxes_accumulate public forcing_SinglePointPrint, mech_forcing_diags, forcing_diagnostics public register_forcing_type_diags, allocate_forcing_type, deallocate_forcing_type @@ -52,6 +52,12 @@ module MOM_forcing_type module procedure allocate_mech_forcing_from_ref end interface allocate_mech_forcing +!> Determine the friction velocity from a forcing type or a mechanical forcing type. +interface find_ustar + module procedure find_ustar_fluxes + module procedure find_ustar_mech_forcing +end interface find_ustar + ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units @@ -1072,6 +1078,139 @@ subroutine calculateBuoyancyFlux2d(G, GV, US, fluxes, optics, h, Temp, Salt, tv, end subroutine calculateBuoyancyFlux2d +!> Determine the friction velocity from the contenxts of a forcing type, perhaps +!! using the evolving surface density. +subroutine find_ustar_fluxes(fluxes, tv, U_star, G, GV, US, halo, H_T_units) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(forcing), intent(in) :: fluxes !< Surface fluxes container + type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any + !! available thermodynamic fields. + real, dimension(SZI_(G),SZJ_(G)), & + intent(out) :: U_star !< The surface friction velocity [Z T-1 ~> m s-1] or + !! [H T-1 ~> m s-1 or kg m-2 s-1], depending on H_T_units. + integer, optional, intent(in) :: halo !< The extra halo size to fill in, 0 by default + logical, optional, intent(in) :: H_T_units !< If present and true, return U_star in units + !! of [H T-1 ~> m s-1 or kg m-2 s-1] + + ! Local variables + real :: I_rho ! The inverse of the reference density times a ratio of scaling + ! factors [Z L-1 R-1 ~> m3 kg-1] or in some semi-Boussinesq cases + ! the rescaled reference density [H2 Z-1 L-1 R-1 ~> m3 kg-1 or kg m-3] + logical :: Z_T_units ! If true, U_star is returned in units of [Z T-1 ~> m s-1], otherwise it is + ! returned in [H T-1 ~> m s-1 or kg m-2 s-1] + integer :: i, j, k, is, ie, js, je, hs + + hs = 0 ; if (present(halo)) hs = max(halo, 0) + is = G%isc - hs ; ie = G%iec + hs ; js = G%jsc - hs ; je = G%jec + hs + + Z_T_units = .true. ; if (present(H_T_units)) Z_T_units = .not.H_T_units + + if (.not.(associated(fluxes%ustar) .or. associated(fluxes%tau_mag))) & + call MOM_error(FATAL, "find_ustar_fluxes requires that either ustar or tau_mag be associated.") + + if (associated(fluxes%ustar) .and. (GV%Boussinesq .or. .not.associated(fluxes%tau_mag))) then + if (Z_T_units) then + do j=js,je ; do i=is,ie + U_star(i,j) = fluxes%ustar(i,j) + enddo ; enddo + else + do j=js,je ; do i=is,ie + U_star(i,j) = GV%Z_to_H * fluxes%ustar(i,j) + enddo ; enddo + endif + elseif (allocated(tv%SpV_avg)) then + if (tv%valid_SpV_halo < 0) call MOM_error(FATAL, & + "find_ustar_fluxes called in non-Boussinesq mode with invalid values of SpV_avg.") + if (tv%valid_SpV_halo < hs) call MOM_error(FATAL, & + "find_ustar_fluxes called in non-Boussinesq mode with insufficient valid values of SpV_avg.") + if (Z_T_units) then + do j=js,je ; do i=is,ie + U_star(i,j) = sqrt(US%L_to_Z*fluxes%tau_mag(i,j) * tv%SpV_avg(i,j,1)) + enddo ; enddo + else + do j=js,je ; do i=is,ie + U_star(i,j) = GV%RZ_to_H * sqrt(US%L_to_Z*fluxes%tau_mag(i,j) / tv%SpV_avg(i,j,1)) + enddo ; enddo + endif + else + I_rho = US%L_to_Z * GV%Z_to_H * GV%RZ_to_H + if (Z_T_units) I_rho = US%L_to_Z * GV%H_to_Z * GV%RZ_to_H ! == US%L_to_Z / GV%Rho0 + do j=js,je ; do i=is,ie + U_star(i,j) = sqrt(fluxes%tau_mag(i,j) * I_rho) + enddo ; enddo + endif + +end subroutine find_ustar_fluxes + + +!> Determine the friction velocity from the contenxts of a forcing type, perhaps +!! using the evolving surface density. +subroutine find_ustar_mech_forcing(forces, tv, U_star, G, GV, US, halo, H_T_units) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(mech_forcing), intent(in) :: forces !< Surface forces container + type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any + !! available thermodynamic fields. + real, dimension(SZI_(G),SZJ_(G)), & + intent(out) :: U_star !< The surface friction velocity [Z T-1 ~> m s-1] + integer, optional, intent(in) :: halo !< The extra halo size to fill in, 0 by default + logical, optional, intent(in) :: H_T_units !< If present and true, return U_star in units + !! of [H T-1 ~> m s-1 or kg m-2 s-1] + + ! Local variables + real :: I_rho ! The inverse of the reference density times a ratio of scaling + ! factors [Z L-1 R-1 ~> m3 kg-1] or in some semi-Boussinesq cases + ! the rescaled reference density [H2 Z-1 L-1 R-1 ~> m3 kg-1 or kg m-3] + logical :: Z_T_units ! If true, U_star is returned in units of [Z T-1 ~> m s-1], otherwise it is + ! returned in [H T-1 ~> m s-1 or kg m-2 s-1] + integer :: i, j, k, is, ie, js, je, hs + + hs = 0 ; if (present(halo)) hs = max(halo, 0) + is = G%isc - hs ; ie = G%iec + hs ; js = G%jsc - hs ; je = G%jec + hs + + Z_T_units = .true. ; if (present(H_T_units)) Z_T_units = .not.H_T_units + + if (.not.(associated(forces%ustar) .or. associated(forces%tau_mag))) & + call MOM_error(FATAL, "find_ustar_mech requires that either ustar or tau_mag be associated.") + + if (associated(forces%ustar) .and. (GV%Boussinesq .or. .not.associated(forces%tau_mag))) then + if (Z_T_units) then + do j=js,je ; do i=is,ie + U_star(i,j) = forces%ustar(i,j) + enddo ; enddo + else + do j=js,je ; do i=is,ie + U_star(i,j) = GV%Z_to_H * forces%ustar(i,j) + enddo ; enddo + endif + elseif (allocated(tv%SpV_avg)) then + if (tv%valid_SpV_halo < 0) call MOM_error(FATAL, & + "find_ustar_mech called in non-Boussinesq mode with invalid values of SpV_avg.") + if (tv%valid_SpV_halo < hs) call MOM_error(FATAL, & + "find_ustar_mech called in non-Boussinesq mode with insufficient valid values of SpV_avg.") + if (Z_T_units) then + do j=js,je ; do i=is,ie + U_star(i,j) = sqrt(US%L_to_Z*forces%tau_mag(i,j) * tv%SpV_avg(i,j,1)) + enddo ; enddo + else + do j=js,je ; do i=is,ie + U_star(i,j) = GV%RZ_to_H * sqrt(US%L_to_Z*forces%tau_mag(i,j) / tv%SpV_avg(i,j,1)) + enddo ; enddo + endif + else + I_rho = US%L_to_Z * GV%Z_to_H * GV%RZ_to_H + if (Z_T_units) I_rho = US%L_to_Z * GV%H_to_Z * GV%RZ_to_H ! == US%L_to_Z / GV%Rho0 + do j=js,je ; do i=is,ie + U_star(i,j) = sqrt(forces%tau_mag(i,j) * I_rho) + enddo ; enddo + endif + +end subroutine find_ustar_mech_forcing + + !> Write out chksums for thermodynamic fluxes. subroutine MOM_forcing_chksum(mesg, fluxes, G, US, haloshift) character(len=*), intent(in) :: mesg !< message diff --git a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 index 206773ecb0..5b7ec60dee 100644 --- a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 +++ b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 @@ -11,7 +11,7 @@ module MOM_mixed_layer_restrat use MOM_error_handler, only : MOM_error, FATAL, WARNING use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_file_parser, only : openParameterBlock, closeParameterBlock -use MOM_forcing_type, only : mech_forcing +use MOM_forcing_type, only : mech_forcing, find_ustar use MOM_grid, only : ocean_grid_type use MOM_hor_index, only : hor_index_type use MOM_lateral_mixing_coeffs, only : VarMix_CS @@ -184,6 +184,9 @@ subroutine mixedlayer_restrat_OM4(h, uhtr, vhtr, tv, forces, dt, MLD_in, VarMix, h_avail ! The volume available for diffusion out of each face of each ! sublayer of the mixed layer, divided by dt [H L2 T-1 ~> m3 s-1 or kg s-1]. real, dimension(SZI_(G),SZJ_(G)) :: & + U_star_2d, & ! 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] MLD_fast, & ! Mixed layer depth actually used in MLE restratification parameterization [H ~> m or kg m-2] htot_fast, & ! The sum of the thicknesses of layers in the mixed layer [H ~> m or kg m-2] Rml_av_fast, & ! g_Rho0 times the average mixed layer density [L2 Z-1 T-2 ~> m s-2] @@ -254,6 +257,9 @@ subroutine mixedlayer_restrat_OM4(h, uhtr, vhtr, tv, forces, dt, MLD_in, VarMix, call MOM_error(FATAL, "mixedlayer_restrat_OM4: "// & "The resolution argument, Rd/dx, was not associated.") + ! Extract the friction velocity from the forcing type. + call find_ustar(forces, tv, U_star_2d, G, GV, US, halo=1) + if (CS%MLE_density_diff > 0.) then ! We need to calculate a mixed layer depth, MLD. !! TODO: use derivatives and mid-MLD pressure. Currently this is sigma-0. -AJA pRef_MLD(:) = 0. @@ -408,7 +414,7 @@ subroutine mixedlayer_restrat_OM4(h, uhtr, vhtr, tv, forces, dt, MLD_in, VarMix, if (CS%debug) then call hchksum(h,'mixed_layer_restrat: h', G%HI, haloshift=1, scale=GV%H_to_m) - call hchksum(forces%ustar,'mixed_layer_restrat: u*', G%HI, haloshift=1, scale=US%Z_to_m*US%s_to_T) + call hchksum(U_star_2d, 'mixed_layer_restrat: u*', G%HI, haloshift=1, scale=US%Z_to_m*US%s_to_T) call hchksum(MLD_fast,'mixed_layer_restrat: MLD', G%HI, haloshift=1, scale=GV%H_to_m) call hchksum(Rml_av_fast,'mixed_layer_restrat: rml', G%HI, haloshift=1, & scale=US%m_to_Z*US%L_T_to_m_s**2) @@ -421,7 +427,7 @@ subroutine mixedlayer_restrat_OM4(h, uhtr, vhtr, tv, forces, dt, MLD_in, VarMix, ! U - Component !$OMP do do j=js,je ; do I=is-1,ie - u_star = max(CS%ustar_min, 0.5*(forces%ustar(i,j) + forces%ustar(i+1,j))) + u_star = max(CS%ustar_min, 0.5*(U_star_2d(i,j) + U_star_2d(i+1,j))) absf = 0.5*(abs(G%CoriolisBu(I,J-1)) + abs(G%CoriolisBu(I,J))) ! If needed, res_scaling_fac = min( ds, L_d ) / l_f @@ -508,7 +514,7 @@ subroutine mixedlayer_restrat_OM4(h, uhtr, vhtr, tv, forces, dt, MLD_in, VarMix, ! V- component !$OMP do do J=js-1,je ; do i=is,ie - u_star = max(CS%ustar_min, 0.5*(forces%ustar(i,j) + forces%ustar(i,j+1))) + u_star = max(CS%ustar_min, 0.5*(U_star_2d(i,j) + U_star_2d(i,j+1))) absf = 0.5*(abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J))) ! If needed, res_scaling_fac = min( ds, L_d ) / l_f @@ -711,6 +717,9 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d wpup ! Turbulent vertical momentum [ ????? ~> m2 s-2] real :: uDml_diag(SZIB_(G),SZJ_(G)) ! A 2D copy of uDml for diagnostics [H L2 T-1 ~> m3 s-1 or kg s-1] real :: vDml_diag(SZI_(G),SZJB_(G)) ! A 2D copy of vDml for diagnostics [H L2 T-1 ~> m3 s-1 or kg s-1] + real :: U_star_2d(SZI_(G),SZJ_(G)) ! 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 :: covTS(SZI_(G)) ! SGS TS covariance in Stanley param; currently 0 [degC ppt] real :: varS(SZI_(G)) ! SGS S variance in Stanley param; currently 0 [ppt2] real :: dmu(SZK_(GV)) ! Change in mu(z) across layer k [nondim] @@ -762,12 +771,15 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d call pass_var(bflux, G%domain, halo=1) + ! Extract the friction velocity from the forcing type. + call find_ustar(forces, tv, U_star_2d, G, GV, US, halo=1) + if (CS%debug) then call hchksum(h,'mixed_Bodner: h', G%HI, haloshift=1, scale=GV%H_to_m) call hchksum(BLD, 'mle_Bodner: BLD in', G%HI, haloshift=1, scale=US%Z_to_m) if (associated(bflux)) & call hchksum(bflux, 'mle_Bodner: bflux', G%HI, haloshift=1, scale=US%Z_to_m**2*US%s_to_T**3) - call hchksum(forces%ustar,'mle_Bodner: u*', G%HI, haloshift=1, scale=US%Z_to_m*US%s_to_T) + call hchksum(U_star_2d, 'mle_Bodner: u*', G%HI, haloshift=1, scale=US%Z_to_m*US%s_to_T) call hchksum(CS%MLD_filtered, 'mle_Bodner: MLD_filtered 1', & G%HI, haloshift=1, scale=US%Z_to_m) call hchksum(CS%MLD_filtered_slow,'mle_Bodner: MLD_filtered_slow 1', & @@ -793,7 +805,7 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d do j = js-1, je+1 ; do i = is-1, ie+1 w_star3 = max(0., -bflux(i,j)) * BLD(i,j) & ! (this line in Z3 T-3 ~> m3 s-3) * ( ( US%Z_to_m * US%s_to_T )**3 ) ! m3 s-3 - u_star3 = ( US%Z_to_m * US%s_to_T * forces%ustar(i,j) )**3 ! m3 s-3 + u_star3 = ( US%Z_to_m * US%s_to_T * U_star_2d(i,j) )**3 ! m3 s-3 wpup(i,j) = max( CS%min_wstar2, & ! The max() avoids division by zero later ( CS%mstar * u_star3 + CS%nstar * w_star3 )**two_thirds ) & ! (this line m2 s-2) * ( ( US%m_to_Z * US%T_to_s )**2 ) ! Z2 T-2 ~> m2 s-2 @@ -965,7 +977,7 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d ! Offer diagnostic fields for averaging. if (query_averaging_enabled(CS%diag)) then - if (CS%id_ustar > 0) call post_data(CS%id_ustar, forces%ustar, CS%diag) + if (CS%id_ustar > 0) call post_data(CS%id_ustar, U_star_2d, CS%diag) if (CS%id_bflux > 0) call post_data(CS%id_bflux, bflux, CS%diag) if (CS%id_wpup > 0) call post_data(CS%id_wpup, wpup, CS%diag) if (CS%id_Rml > 0) call post_data(CS%id_Rml, buoy_av, CS%diag) @@ -1053,6 +1065,9 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS) h_avail ! The volume available for diffusion out of each face of each ! sublayer of the mixed layer, divided by dt [H L2 T-1 ~> m3 s-1 or kg s-1]. real, dimension(SZI_(G),SZJ_(G)) :: & + U_star_2d, & ! 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] htot, & ! The sum of the thicknesses of layers in the mixed layer [H ~> m or kg m-2] Rml_av ! g_Rho0 times the average mixed layer density [L2 Z-1 T-2 ~> m s-2] real :: g_Rho0 ! G_Earth/Rho0 [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1] @@ -1110,6 +1125,9 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS) if (CS%use_Stanley_ML) call MOM_error(FATAL, "mixedlayer_restrat_BML: "// & "The Stanley parameterization is not available with the BML.") + ! Extract the friction velocity from the forcing type. + call find_ustar(forces, tv, U_star_2d, G, GV, US, halo=1) + ! Fix this later for nkml >= 3. p0(:) = 0.0 @@ -1145,7 +1163,7 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS) do j=js,je ; do I=is-1,ie h_vel = 0.5*(htot(i,j) + htot(i+1,j)) * GV%H_to_Z - u_star = max(CS%ustar_min, 0.5*(forces%ustar(i,j) + forces%ustar(i+1,j))) + u_star = max(CS%ustar_min, 0.5*(U_star_2d(i,j) + U_star_2d(i+1,j))) absf = 0.5*(abs(G%CoriolisBu(I,J-1)) + abs(G%CoriolisBu(I,J))) @@ -1196,7 +1214,7 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS) do J=js-1,je ; do i=is,ie h_vel = 0.5*(htot(i,j) + htot(i,j+1)) * GV%H_to_Z - u_star = max(CS%ustar_min, 0.5*(forces%ustar(i,j) + forces%ustar(i,j+1))) + u_star = max(CS%ustar_min, 0.5*(U_star_2d(i,j) + U_star_2d(i,j+1))) absf = 0.5*(abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J))) diff --git a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 index 66e2dfa6b2..ceba8dad1a 100644 --- a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 +++ b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 @@ -9,7 +9,7 @@ module MOM_bulk_mixed_layer use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type use MOM_error_handler, only : MOM_error, FATAL, WARNING use MOM_file_parser, only : get_param, log_param, log_version, param_file_type -use MOM_forcing_type, only : extractFluxes1d, forcing +use MOM_forcing_type, only : extractFluxes1d, forcing, find_ustar use MOM_grid, only : ocean_grid_type use MOM_opacity, only : absorbRemainingSW, optics_type, extract_optics_slice use MOM_unit_scaling, only : unit_scale_type @@ -235,6 +235,9 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, C ksort ! The sorted k-index that each original layer goes to. real, dimension(SZI_(G),SZJ_(G)) :: & h_miss ! The summed absolute mismatch [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G)) :: & + U_star_2d ! 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_(G)) :: & TKE, & ! The turbulent kinetic energy available for mixing over a ! time step [Z L2 T-2 ~> m3 s-2]. @@ -412,6 +415,9 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, C max_BL_det(:) = -1 EOSdom(:) = EOS_domain(G%HI) + ! Extract the friction velocity from the forcing type. + call find_ustar(fluxes, tv, U_star_2d, G, GV, US) + !$OMP parallel default(shared) firstprivate(dKE_CA,cTKE,h_CA,max_BL_det,p_ref,p_ref_cv) & !$OMP private(h,u,v,h_orig,eps,T,S,opacity_band,d_ea,d_eb,R0,Rcv,ksort, & !$OMP dR0_dT,dR0_dS,dRcv_dT,dRcv_dS,htot,Ttot,Stot,TKE,Conv_en, & @@ -513,7 +519,7 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, C ! First the TKE at the depth of free convection that is available ! to drive mixing is calculated. - call find_starting_TKE(htot, h_CA, fluxes, Conv_En, cTKE, dKE_FC, dKE_CA, & + call find_starting_TKE(htot, h_CA, fluxes, U_star_2d, Conv_En, cTKE, dKE_FC, dKE_CA, & TKE, TKE_river, Idecay_len_TKE, cMKE, dt, Idt_diag, & j, ksort, G, GV, US, CS) @@ -1252,7 +1258,7 @@ end subroutine mixedlayer_convection !> This subroutine determines the TKE available at the depth of free !! convection to drive mechanical entrainment. -subroutine find_starting_TKE(htot, h_CA, fluxes, Conv_En, cTKE, dKE_FC, dKE_CA, & +subroutine find_starting_TKE(htot, h_CA, fluxes, U_star_2d, Conv_En, cTKE, dKE_FC, dKE_CA, & TKE, TKE_river, Idecay_len_TKE, cMKE, dt, Idt_diag, & j, ksort, G, GV, US, CS) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. @@ -1265,6 +1271,10 @@ subroutine find_starting_TKE(htot, h_CA, fluxes, Conv_En, cTKE, dKE_FC, dKE_CA, type(forcing), intent(in) :: fluxes !< A structure containing pointers to any !! possible forcing fields. Unused fields !! have NULL pointers. + real, dimension(SZI_(G),SZJ_(G)), intent(in) :: U_star_2d !< 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_(G)), intent(inout) :: Conv_En !< The buoyant turbulent kinetic energy source !! due to free convection [Z L2 T-2 ~> m3 s-2]. real, dimension(SZI_(G)), intent(in) :: dKE_FC !< The vertically integrated change in @@ -1325,7 +1335,8 @@ subroutine find_starting_TKE(htot, h_CA, fluxes, Conv_En, cTKE, dKE_FC, dKE_CA, if (CS%omega_frac >= 1.0) absf = 2.0*CS%omega do i=is,ie - U_star = fluxes%ustar(i,j) + U_star = U_star_2d(i,j) + if (associated(fluxes%ustar_shelf) .and. associated(fluxes%frac_shelf_h)) then if (fluxes%frac_shelf_h(i,j) > 0.0) & U_star = (1.0 - fluxes%frac_shelf_h(i,j)) * U_star + &