From b913ba79ece0af0a0399bc1635badc299408f945 Mon Sep 17 00:00:00 2001 From: Brandon Reichl Date: Mon, 19 Dec 2016 13:43:52 -0500 Subject: [PATCH 1/9] Working on adding mstar as function of MLD, Ekman Depth, and Obhukov depth --- .../lateral/MOM_hor_visc.F90 | 4 +- .../vertical/MOM_diabatic_driver.F90 | 10 +- .../vertical/MOM_energetic_PBL.F90 | 122 ++++++++++++++---- 3 files changed, 103 insertions(+), 33 deletions(-) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 1a6e11d8fb..e06fae5a7c 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -207,7 +207,7 @@ module MOM_hor_visc contains subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, OBC) - type(ocean_grid_type), intent(in) :: G + type(ocean_grid_type), intent(inout) :: G type(verticalGrid_type), intent(in) :: GV real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v @@ -703,7 +703,7 @@ end subroutine horizontal_viscosity subroutine hor_visc_init(Time, G, param_file, diag, CS) type(time_type), intent(in) :: Time - type(ocean_grid_type), intent(in) :: G + type(ocean_grid_type), intent(inout) :: G type(param_file_type), intent(in) :: param_file type(diag_ctrl), target, intent(inout) :: diag type(hor_visc_CS), pointer :: CS diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 0ffc550558..24c7140fd3 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -742,6 +742,8 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, GV, CS) if (CS%use_energetic_PBL) then call applyBoundaryFluxesInOut(CS%diabatic_aux_CSp, G, GV, dt, fluxes, CS%optics, & h, tv, CS%aggregate_FW_forcing, cTKE, dSV_dT, dSV_dS) + call calculateBuoyancyFlux2d(G, GV, fluxes, CS%optics, h, tv%T, tv%S, tv, & + CS%KPP_buoy_flux, CS%KPP_temp_flux, CS%KPP_salt_flux) if (CS%debug) then call hchksum(ea, "after applyBoundaryFluxes ea",G%HI,haloshift=0) @@ -753,7 +755,7 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, GV, CS) call find_uv_at_h(u, v, h, u_h, v_h, G, GV) call energetic_PBL(h, u_h, v_h, tv, fluxes, dt, Kd_ePBL, G, GV, & - CS%energetic_PBL_CSp, dSV_dT, dSV_dS, cTKE) + CS%energetic_PBL_CSp, dSV_dT, dSV_dS, cTKE, CS%KPP_Buoy_Flux) ! If visc%MLD exists, copy the ePBL's MLD into it if (associated(visc%MLD)) then @@ -2029,13 +2031,13 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, ! CS%useKPP is set to True if KPP-scheme is to be used, False otherwise. ! KPP_init() allocated CS%KPP_Csp and also sets CS%KPPisPassive CS%useKPP = KPP_init(param_file, G, diag, Time, CS%KPP_CSp, passive=CS%KPPisPassive) - if (CS%useKPP) then + if (CS%useKPP .or. CS%use_energetic_PBL) then allocate( CS%KPP_NLTheat(isd:ied,jsd:jed,nz+1) ) ; CS%KPP_NLTheat(:,:,:) = 0. allocate( CS%KPP_NLTscalar(isd:ied,jsd:jed,nz+1) ) ; CS%KPP_NLTscalar(:,:,:) = 0. allocate( CS%KPP_buoy_flux(isd:ied,jsd:jed,nz+1) ) ; CS%KPP_buoy_flux(:,:,:) = 0. allocate( CS%KPP_temp_flux(isd:ied,jsd:jed) ) ; CS%KPP_temp_flux(:,:) = 0. allocate( CS%KPP_salt_flux(isd:ied,jsd:jed) ) ; CS%KPP_salt_flux(:,:) = 0. - endif + endif call get_param(param_file, mod, "SALT_REJECT_BELOW_ML", CS%salt_reject_below_ML, & "If true, place salt from brine rejection below the mixed layer,\n"// & @@ -2272,7 +2274,7 @@ subroutine diabatic_driver_end(CS) call entrain_diffusive_end(CS%entrain_diffusive_CSp) call set_diffusivity_end(CS%set_diff_CSp) - if (CS%useKPP) then + if (CS%useKPP .or. CS%use_energetic_PBL) then deallocate( CS%KPP_NLTheat ) deallocate( CS%KPP_NLTscalar ) deallocate( CS%KPP_buoy_flux ) diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index 70fb49f728..b2308dcb0d 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -87,12 +87,22 @@ module MOM_energetic_PBL public energetic_PBL_get_MLD type, public :: energetic_PBL_CS ; private + logical :: MSTARFIXED real :: mstar ! The ratio of the friction velocity cubed to the - ! TKE input to the mixed layer, nondimensional. + ! TKE available to drive entrainment, nondimensional. + ! This quantity is the vertically integrated + ! shear production minus the vertically integrated + ! dissipation of TKE produced by shear. real :: nstar ! The fraction of the TKE input to the mixed layer ! available to drive entrainment, nondim. + ! This quantity is the vertically integrated + ! buoyancy production minus the vertically integrated + ! dissipation of TKE produced by buoyancy. real :: TKE_decay ! The ratio of the natural Ekman depth to the TKE ! decay scale, nondimensional. + ! This process is generally not included in models + ! and therefore cannot be found from (for example) + ! tuning against K-Epsilon. real :: MKE_to_TKE_effic ! The efficiency with which mean kinetic energy ! released by mechanically forced entrainment of ! the mixed layer is converted to TKE, nondim. @@ -178,8 +188,8 @@ module MOM_energetic_PBL contains subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & - dSV_dT, dSV_dS, TKE_forced, dt_diag, last_call, & - dT_expected, dS_expected) + dSV_dT, dSV_dS, TKE_forced, Buoy_Flux, dt_diag, last_call, & + dT_expected, dS_expected ) type(ocean_grid_type), intent(inout) :: G type(verticalGrid_type), intent(in) :: GV real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h_3d @@ -190,6 +200,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & real, intent(in) :: dt real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(out) :: Kd_int type(energetic_PBL_CS), pointer :: CS + real, dimension(SZI_(G),SZJ_(G)), intent(in) :: Buoy_Flux real, optional, intent(in) :: dt_diag logical, optional, intent(in) :: last_call real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & @@ -448,6 +459,8 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & integer, save :: NOTCONVERGED! !-End BGR iteration parameters----------------------------------------- real :: N2_dissipation + real :: B_STAR + real :: Stab, StabL, mstar_effective type(cvmix_global_params_type) :: CVMIX_gravity logical :: debug=.false. ! Change this hard-coded value for debugging. ! The following arrays are used only for debugging purposes. @@ -504,6 +517,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & if (CS%Mixing_Diagnostics) then CS%Mixing_Length(:,:,:) = 0.0 CS%Velocity_Scale(:,:,:) = 0.0 + !CS%Mstar_Used(:,:) = 0.0 endif !!OMP end parallel endif @@ -571,36 +585,42 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & fluxes%frac_shelf_h(i,j) * fluxes%ustar_shelf(i,j) endif if (U_Star < CS%ustar_min) U_Star = CS%ustar_min + ! Computing u10 based on u_star and COARE 3.5 relationships call ust_2_u10_coare3p5(U_STAR*sqrt(GV%Rho0/1.225),U_10) + ! Computing B_Star, or the Buoyancy flux over the friction velocity. + B_Star = min(0.0,-buoy_Flux(i,j)/U_Star) if (CS%omega_frac >= 1.0) then ; absf(i) = 2.0*CS%omega else absf(i) = 0.25*((abs(G%CoriolisBu(I,J)) + abs(G%CoriolisBu(I-1,J-1))) + & (abs(G%CoriolisBu(I,J-1)) + abs(G%CoriolisBu(I-1,J)))) if (CS%omega_frac > 0.0) & - absf = sqrt(CS%omega_frac*4.0*CS%omega**2 + (1.0-CS%omega_frac)*absf**2) + absf(i) = sqrt(CS%omega_frac*4.0*CS%omega**2 + (1.0-CS%omega_frac)*absf(i)**2) endif - mech_TKE(i) = (dt*CS%mstar*GV%Rho0)*((U_Star**3)) - conv_PErel(i) = 0.0 + if (CS%mstarFixed) then + mech_TKE(i) = (dt*CS%mstar*GV%Rho0)*((U_Star**3)) + conv_PErel(i) = 0.0 + + if (CS%TKE_diagnostics) then + CS%diag_TKE_wind(i,j) = CS%diag_TKE_wind(i,j) + mech_TKE(i) * IdtdR0 + if (TKE_forced(i,j,1) <= 0.0) then + CS%diag_TKE_forcing(i,j) = CS%diag_TKE_forcing(i,j) + & + max(-mech_TKE(i), TKE_forced(i,j,1)) * IdtdR0 + ! CS%diag_TKE_unbalanced_forcing(i,j) = CS%diag_TKE_unbalanced_forcing(i,j) + & + ! min(0.0, TKE_forced(i,j,1) + mech_TKE(i)) * IdtdR0 + else + CS%diag_TKE_forcing(i,j) = CS%diag_TKE_forcing(i,j) + CS%nstar*TKE_forced(i,j,1) * IdtdR0 + endif + endif - if (CS%TKE_diagnostics) then - CS%diag_TKE_wind(i,j) = CS%diag_TKE_wind(i,j) + mech_TKE(i) * IdtdR0 if (TKE_forced(i,j,1) <= 0.0) then - CS%diag_TKE_forcing(i,j) = CS%diag_TKE_forcing(i,j) + & - max(-mech_TKE(i), TKE_forced(i,j,1)) * IdtdR0 - ! CS%diag_TKE_unbalanced_forcing(i,j) = CS%diag_TKE_unbalanced_forcing(i,j) + & - ! min(0.0, TKE_forced(i,j,1) + mech_TKE(i)) * IdtdR0 + mech_TKE(i) = mech_TKE(i) + TKE_forced(i,j,1) + if (mech_TKE(i) < 0.0) mech_TKE(i) = 0.0 else - CS%diag_TKE_forcing(i,j) = CS%diag_TKE_forcing(i,j) + CS%nstar*TKE_forced(i,j,1) * IdtdR0 + conv_PErel(i) = conv_PErel(i) + TKE_forced(i,j,1) endif - endif - if (TKE_forced(i,j,1) <= 0.0) then - mech_TKE(i) = mech_TKE(i) + TKE_forced(i,j,1) - if (mech_TKE(i) < 0.0) mech_TKE(i) = 0.0 - else - conv_PErel(i) = conv_PErel(i) + TKE_forced(i,j,1) endif ! endif ; enddo @@ -636,7 +656,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & ! Store the initial mechanical TKE and convectively released PE to ! enable multiple iterations. - mech_TKE_top(i) = mech_TKE(i) ; conv_PErel_top(i) = conv_PErel(i) + !mech_TKE_top(i) = mech_TKE(i) ; conv_PErel_top(i) = conv_PErel(i) !/The following lines are for the iteration over MLD !{ @@ -668,13 +688,57 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & !CS%ML_depth2(i,j) = h(i,1)*GV%H_to_m sfc_connected(i) = .true. - ! Multiply mech TKE at surface by enhancement (w'/ust) cubed. - ! This has not been confirmed to reproduce LES Langmuir simulations, - ! but is consistent if we assume we are trying to capture the enhancement - ! of the TKE flux/ (turbulent velocity scale cubed) knowing that the - ! enhancement factor is that of the turbulent velocity scale. - mech_TKE(i) = mech_TKE_top(i)*ENHANCE_V**(3.) ; conv_PErel(i) = conv_PErel_top(i) + if (.not.CS%MstarFixed) then + ! Note this value of mech_TKE(i) now must be iterated over, so it is moved here + + ! First solve for the TKE to PE length scale + ! Change the (2.0) to a runtime parameter + StabL = -u_star**2 / ( VonKar * ( B_STAR + (2.0) * u_star * absf(i))) + Stab = MLD_guess / StabL + ! Limit Stab to regime where 3rd order fit is reasonable + Stab = max(-0.79,min(5.0,Stab)) + Mstar_effective = 0.005007927814101 * Stab**4 + & + -0.059699214537634 * Stab**3 + & + 0.259989170346221 * Stab**2 + & + 0.706210570634930 * Stab + & + 0.367256544478824 + if ((Stab).ge.0.0 .and. STAB .le.0.5)then + Mstar_effective=Mstar_effective*(1.+.2*((0.5-STAB)/0.5)) + elseif ((Stab).lt.0.0 .and. STAB .gt.-5.0) then + Mstar_effective=Mstar_effective*(1.-.3*((5.0+STAB)/5.0)) + endif + mech_TKE(i) = (dt*Mstar_effective*GV%Rho0)*((U_Star**3)) + conv_PErel(i) = 0.0 + print*,'---' + print*,Stab,Mstar_effective + print*,MLD_guess,-u_star**2/b_star,-2.*u_star/absf(i) + if (CS%TKE_diagnostics) then + CS%diag_TKE_wind(i,j) = CS%diag_TKE_wind(i,j) + mech_TKE(i) * IdtdR0 + if (TKE_forced(i,j,1) <= 0.0) then + CS%diag_TKE_forcing(i,j) = CS%diag_TKE_forcing(i,j) + & + max(-mech_TKE(i), TKE_forced(i,j,1)) * IdtdR0 + ! CS%diag_TKE_unbalanced_forcing(i,j) = CS%diag_TKE_unbalanced_forcing(i,j) + & + ! min(0.0, TKE_forced(i,j,1) + mech_TKE(i)) * IdtdR0 + else + CS%diag_TKE_forcing(i,j) = CS%diag_TKE_forcing(i,j) + CS%nstar*TKE_forced(i,j,1) * IdtdR0 + endif + endif + + if (TKE_forced(i,j,1) <= 0.0) then + mech_TKE(i) = mech_TKE(i) + TKE_forced(i,j,1) + if (mech_TKE(i) < 0.0) mech_TKE(i) = 0.0 + else + conv_PErel(i) = conv_PErel(i) + TKE_forced(i,j,1) + endif + else + ! Multiply mech TKE at surface by enhancement (w'/ust) cubed. + ! This has not been confirmed to reproduce LES Langmuir simulations, + ! but is consistent if we assume we are trying to capture the enhancement + ! of the TKE flux/ (turbulent velocity scale cubed) knowing that the + ! enhancement factor is that of the turbulent velocity scale. + mech_TKE(i) = mech_TKE_top(i)*ENHANCE_V**(3.) ; conv_PErel(i) = conv_PErel_top(i) + endif if (CS%TKE_diagnostics) then dTKE_conv = 0.0 ; dTKE_forcing = 0.0 ; dTKE_mixing = 0.0 @@ -705,7 +769,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & do K=2,nz+1 h_rsum = h_rsum + h(i,k-1) MixLen_shape(K) = CS%transLay_scale + (1.0 - CS%transLay_scale) * & - (max(0.0, (MLD_guess - h_rsum)*I_MLD) )**2 !BR prefer 1 or 2 for exponent? + (max(0.0, (MLD_guess - h_rsum)*I_MLD) )**1.0 !BR prefer 1 or 2 for exponent? enddo endif @@ -1739,6 +1803,10 @@ subroutine energetic_PBL_init(Time, G, GV, param_file, diag, CS) call get_param(param_file, mod, "MSTAR", CS%mstar, & "The ratio of the friction velocity cubed to the TKE \n"//& "input to the mixed layer.", "units=nondim", default=1.2) + call get_param(param_file, mod, "MSTAR_FIXED", CS%MstarFixed, & + "True to use the fixed value of mstar, false to depend \n"//& + "on the Obhukov length and Ekman length.","units=nondim",& + default=.true.) call get_param(param_file, mod, "NSTAR", CS%nstar, & "The portion of the buoyant potential energy imparted by \n"//& "surface fluxes that is available to drive entrainment \n"//& From a6b846f7adc0fdf33b8312174b65d3a107e20124 Mon Sep 17 00:00:00 2001 From: Brandon Reichl Date: Tue, 3 Jan 2017 17:19:19 -0500 Subject: [PATCH 2/9] Updates to ePBL before merging from dev/master. --- .../vertical/MOM_energetic_PBL.F90 | 161 ++++++++++-------- 1 file changed, 92 insertions(+), 69 deletions(-) diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index b2308dcb0d..f37dca8bc4 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -87,7 +87,6 @@ module MOM_energetic_PBL public energetic_PBL_get_MLD type, public :: energetic_PBL_CS ; private - logical :: MSTARFIXED real :: mstar ! The ratio of the friction velocity cubed to the ! TKE available to drive entrainment, nondimensional. ! This quantity is the vertically integrated @@ -100,9 +99,6 @@ module MOM_energetic_PBL ! dissipation of TKE produced by buoyancy. real :: TKE_decay ! The ratio of the natural Ekman depth to the TKE ! decay scale, nondimensional. - ! This process is generally not included in models - ! and therefore cannot be found from (for example) - ! tuning against K-Epsilon. real :: MKE_to_TKE_effic ! The efficiency with which mean kinetic energy ! released by mechanically forced entrainment of ! the mixed layer is converted to TKE, nondim. @@ -146,7 +142,15 @@ module MOM_energetic_PBL ! applies to in subsequent revisions of this code. ! "_Neg" and "_Pos" refer to which scale is applied as a ! function of negative or positive local buoyancy. + real :: MSTAR_SLOPE ! Slope of the function which relates the shear production + ! to the mixing layer depth, Ekman depth, and Monin-Obukhov + ! depth. + real :: MSTAR_XINT ! Value where MSTAR function transitions from linear + ! to decay toward MSTAR->0 at fully developed Ekman depth. + real :: MSTAR_AT_XINT ! Intercept value of MSTAR at value where function + ! changes to linear transition. type(time_type), pointer :: Time ! A pointer to the ocean model's clock. + logical :: Use_Mstar_Fixed = .true. ! A logical to revert to a fixed m* logical :: TKE_diagnostics = .false. logical :: Use_LT_LiFunction = .false. logical :: orig_PE_calc = .true. @@ -248,6 +252,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & ! salinity, in m3 kg-1 ppt-1. ! (in) TKE_forced - The forcing requirements to homogenize the forcing ! that has been applied to each layer through each layer, in J m-2. +! (in) Buoy_Flux - The surface buoyancy flux. in m2/s3. ! (in,opt) dt_diag - The diagnostic time step, which may be less than dt ! if there are two callse to mixedlayer, in s. ! (in,opt) last_call - if true, this is the last call to mixedlayer in the @@ -459,8 +464,19 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & integer, save :: NOTCONVERGED! !-End BGR iteration parameters----------------------------------------- real :: N2_dissipation - real :: B_STAR - real :: Stab, StabL, mstar_effective + real :: B_STAR ! Buoyancy flux over ustar + real :: STAB_SCALE ! Composite of Stabilizing length scales: + ! Ekman scale and Monin-Obukhov scale. + real :: C_MO = 1. ! Constant in STAB_SCALE for Monin-Obukhov + real :: C_EK = 2. ! Constant in STAB_SCALE for Ekman length + real :: MLD_over_STAB ! Mixing layer depth divided by STAB_SCALE + real :: MSTAR_N = -2. ! Exponent in hyperolic decay at negative values of MLD_over_STAB + real :: MSTAR_A ! MSTAR_A and MSTAR_B are coefficients in hyperbolic decay which + real :: MSTAR_B ! are computed to match the value and slope of the linear fit at the + ! value of mld_over_stab where the transition is set. + real :: MSTAR_MIX ! The value of mstar (Proportionality of TKE to drive mixing to ustar + ! cubed) which is computed as a function of latitude, boundary layer depth, + ! and the Monin-Obhukov depth. type(cvmix_global_params_type) :: CVMIX_gravity logical :: debug=.false. ! Change this hard-coded value for debugging. ! The following arrays are used only for debugging purposes. @@ -517,7 +533,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & if (CS%Mixing_Diagnostics) then CS%Mixing_Length(:,:,:) = 0.0 CS%Velocity_Scale(:,:,:) = 0.0 - !CS%Mstar_Used(:,:) = 0.0 + !CS%Mstar_Used(:,:) = 0.0 ! Should add diagnostic Mstar_Used endif !!OMP end parallel endif @@ -589,7 +605,9 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & call ust_2_u10_coare3p5(U_STAR*sqrt(GV%Rho0/1.225),U_10) ! Computing B_Star, or the Buoyancy flux over the friction velocity. B_Star = min(0.0,-buoy_Flux(i,j)/U_Star) - + ! Computing stability scale which correlates with TKE for mixing, where + ! TKE for mixing = TKE production minus TKE dissipation + Stab_Scale = -u_star**2 / ( VonKar * ( C_MO * B_STAR + C_EK * u_star * absf(i))) if (CS%omega_frac >= 1.0) then ; absf(i) = 2.0*CS%omega else absf(i) = 0.25*((abs(G%CoriolisBu(I,J)) + abs(G%CoriolisBu(I-1,J-1))) + & @@ -598,7 +616,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & absf(i) = sqrt(CS%omega_frac*4.0*CS%omega**2 + (1.0-CS%omega_frac)*absf(i)**2) endif - if (CS%mstarFixed) then + if (CS%Use_Mstar_Fixed) then mech_TKE(i) = (dt*CS%mstar*GV%Rho0)*((U_Star**3)) conv_PErel(i) = 0.0 @@ -656,7 +674,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & ! Store the initial mechanical TKE and convectively released PE to ! enable multiple iterations. - !mech_TKE_top(i) = mech_TKE(i) ; conv_PErel_top(i) = conv_PErel(i) + mech_TKE_top(i) = mech_TKE(i) ; conv_PErel_top(i) = conv_PErel(i) !/The following lines are for the iteration over MLD !{ @@ -673,7 +691,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & MLD_guess=CS%ML_Depth2(i,j) else !Otherwise guess middle of water column. - MLD_guess = 0.5 * (min_MLD+max_MLD) + MLD_guess = min(abs(stab_scale),0.5 * (min_MLD+max_MLD)) endif ! Iterate up to MAX_OBL_IT times to determine a converged EPBL depth. @@ -689,30 +707,22 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & sfc_connected(i) = .true. - if (.not.CS%MstarFixed) then + if (.not.CS%Use_Mstar_Fixed) then ! Note this value of mech_TKE(i) now must be iterated over, so it is moved here ! First solve for the TKE to PE length scale ! Change the (2.0) to a runtime parameter - StabL = -u_star**2 / ( VonKar * ( B_STAR + (2.0) * u_star * absf(i))) - Stab = MLD_guess / StabL - ! Limit Stab to regime where 3rd order fit is reasonable - Stab = max(-0.79,min(5.0,Stab)) - Mstar_effective = 0.005007927814101 * Stab**4 + & - -0.059699214537634 * Stab**3 + & - 0.259989170346221 * Stab**2 + & - 0.706210570634930 * Stab + & - 0.367256544478824 - if ((Stab).ge.0.0 .and. STAB .le.0.5)then - Mstar_effective=Mstar_effective*(1.+.2*((0.5-STAB)/0.5)) - elseif ((Stab).lt.0.0 .and. STAB .gt.-5.0) then - Mstar_effective=Mstar_effective*(1.-.3*((5.0+STAB)/5.0)) + MLD_over_Stab = 1.0 * MLD_guess / Stab_Scale + MSTAR_A = CS%MSTAR_AT_XINT**(1./MSTAR_N) + MSTAR_B = CS%MSTAR_SLOPE / (MSTAR_N*MSTAR_A**(MSTAR_N-1.)) + if ((MLD_over_Stab) .ge. CS%MSTAR_XINT) then + MSTAR_mix = CS%MSTAR_SLOPE*(MLD_over_Stab-CS%MSTAR_XINT)+CS%MSTAR_AT_XINT + else + MSTAR_mix = (MSTAR_B*(MLD_over_Stab-CS%MSTAR_XINT)+MSTAR_A)**(MSTAR_N) endif - mech_TKE(i) = (dt*Mstar_effective*GV%Rho0)*((U_Star**3)) + !print*,MLD_guess,STAB_SCALE,MSTAR_MIX + mech_TKE(i) = (dt*MSTAR_mix*GV%Rho0)*((U_Star**3)) conv_PErel(i) = 0.0 - print*,'---' - print*,Stab,Mstar_effective - print*,MLD_guess,-u_star**2/b_star,-2.*u_star/absf(i) if (CS%TKE_diagnostics) then CS%diag_TKE_wind(i,j) = CS%diag_TKE_wind(i,j) + mech_TKE(i) * IdtdR0 if (TKE_forced(i,j,1) <= 0.0) then @@ -1271,47 +1281,48 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & !1. Check if guess was too shallow !Adding -1 m as cushion, will help avoid ! non-convergence flag when nearly converged. - elseif (MLD_FOUND-CS%MLD_tol > MLD_guess) then - !elseif (MLD_FOUND > MLD_guess) then - !/ Guess was too shallow, set new minimum guess - min_MLD = MLD_guess - FIRST_OBL = .false. !Break OBL loop - !2. Check if Guess minus found MLD - ! is less than thickness of level (=converged) - ! - We could try to add a more precise - ! value for found MLD, but seems difficult to - ! to contrain beyond within a level. - elseif ((MLD_guess-MLD_FOUND) < max(CS%MLD_tol,h(i,k-1)*GV%H_to_m)) then - ! Converged. Exit iteration. - FIRST_OBL = .false.!Break OBL loop - OBL_CONVERGED = .true.!Break convergence loop - ! Testing Output, comment for use. - !print*,'Converged--------' - !print*,MLD_FOUND,MLD_guess - !/ - if (OBL_IT_STATS) then !Compute iteration statistics - MAXIT = max(MAXIT,obl_it) - MINIT = min(MINIT,obl_it) - SUMIT = SUMIT+obl_it - NUMIT = NUMIT+1 - print*,MAXIT,MINIT,SUMIT/NUMIT - endif - !BGR this can be where MLD is stored for next time... - CS%ML_Depth2(i,j) = MLD_guess - !/ - !2. If not, guess was too deep - else - !Guess was too deep, set new maximum guess - max_MLD = MLD_guess !We know this guess was too deep - FIRST_OBL = .false.!Break OBL loop - endif - endif - enddo + else + FIRST_OBL = .false. + endif + endif + enddo + MLD_FOUND=CS%ML_DEPTH(i,j) + if (MLD_FOUND-CS%MLD_tol > MLD_guess) then + !elseif (MLD_FOUND > MLD_guess) then + !/ Guess was too shallow, set new minimum guess + min_MLD = MLD_guess + !2. Check if Guess minus found MLD + ! is less than thickness of level (=converged) + ! - We could try to add a more precise + ! value for found MLD, but seems difficult to + ! to contrain beyond within a level. + elseif ((MLD_guess-MLD_FOUND) < max(CS%MLD_tol,h(i,k-1)*GV%H_to_m)) then + ! Converged. Exit iteration. + OBL_CONVERGED = .true.!Break convergence loop + ! Testing Output, comment for use. + !print*,'Converged--------' + !print*,MLD_FOUND,MLD_guess + !/ + if (OBL_IT_STATS) then !Compute iteration statistics + MAXIT = max(MAXIT,obl_it) + MINIT = min(MINIT,obl_it) + SUMIT = SUMIT+obl_it + NUMIT = NUMIT+1 + print*,MAXIT,MINIT,SUMIT/NUMIT + endif + !BGR this can be where MLD is stored for next time... + CS%ML_Depth2(i,j) = MLD_guess + !/ + !2. If not, guess was too deep + else + !Guess was too deep, set new maximum guess + max_MLD = MLD_guess !We know this guess was too deep + endif ! For next pass, guess average of minimum and maximum values. - MLD_guess = min_MLD*0.5 + max_MLD*0.5 + MLD_guess = MLD_FOUND!min_MLD*0.5 + max_MLD*0.5 ITresult(obl_it) = MLD_FOUND endif ; enddo ! Iteration loop for converged boundary layer thickness. - + print*,mld_guess,CS%ml_depth(i,j) if (.not.OBL_CONVERGED) then !/Temp output, warn that EPBL didn't converge !/Print guess/found for every iteration step @@ -1803,8 +1814,20 @@ subroutine energetic_PBL_init(Time, G, GV, param_file, diag, CS) call get_param(param_file, mod, "MSTAR", CS%mstar, & "The ratio of the friction velocity cubed to the TKE \n"//& "input to the mixed layer.", "units=nondim", default=1.2) - call get_param(param_file, mod, "MSTAR_FIXED", CS%MstarFixed, & - "True to use the fixed value of mstar, false to depend \n"//& + call get_param(param_file, mod, "MSTAR_SLOPE", CS%mstar_slope, & + "The ratio of the friction velocity cubed to the TKE \n"//& + "input to the mixed layer (used if MSTAR_FIXED=false).",& + "units=nondim", default=1.2) + call get_param(param_file, mod, "MSTAR_XINT", CS%mstar_xint, & + "The ratio of the friction velocity cubed to the TKE \n"//& + "input to the mixed layer (used if MSTAR_FIXED=false).",& + "units=nondim", default=1.2) + call get_param(param_file, mod, "MSTAR_AT_XINT", CS%mstar_at_xint, & + "The ratio of the friction velocity cubed to the TKE \n"//& + "input to the mixed layer (used if MSTAR_FIXED=false).",& + "units=nondim", default=1.2) + call get_param(param_file, mod, "MSTAR_FIXED", CS%Use_Mstar_Fixed, & + "True to use a fixed value of mstar, if false mstar depends \n"//& "on the Obhukov length and Ekman length.","units=nondim",& default=.true.) call get_param(param_file, mod, "NSTAR", CS%nstar, & From e4c8d572f3e465367c366edf5781313f2e591c79 Mon Sep 17 00:00:00 2001 From: Brandon Reichl Date: Tue, 3 Jan 2017 18:08:18 -0500 Subject: [PATCH 3/9] New EPBL version where MSTAR no longer must be treated as a constant. - Adds pass of buoyancy flux to EPBL (needed to compute MO length scale) - Adds calculation of variable MSTAR to EPBL (option is added to revert to use constant) - MSTAR calculation is based on MLD and composite Ekman/Monin-Obuhkov length scale. - Relationship between the length scale and MSTAR is determined to be roughly linear and optional inputs for line parameters are added (with defaults based on fit to K-epsilon simulations). - Reworks EPBL MLD iteration to use MLD value determined by EPBL rather than simple value based on vertical TKE decay (this may require more robust testing but it is a necessary change since MLD must now be estimated more precisely). - Adds input for exponent in EPBL mixing length (default 1, was 2) --- .../lateral/MOM_hor_visc.F90 | 2 +- .../vertical/MOM_energetic_PBL.F90 | 126 +++++++++--------- 2 files changed, 67 insertions(+), 61 deletions(-) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index a3dd791b3d..5edf85af79 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -207,7 +207,7 @@ module MOM_hor_visc contains subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, CS, OBC) - type(ocean_grid_type), intent(inout) :: G + type(ocean_grid_type), intent(in) :: G type(verticalGrid_type), intent(in) :: GV real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index f37dca8bc4..9abf0f2149 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -97,6 +97,9 @@ module MOM_energetic_PBL ! This quantity is the vertically integrated ! buoyancy production minus the vertically integrated ! dissipation of TKE produced by buoyancy. + real :: MixLenExponent ! Exponent in the mixing length shape-function. + ! 1 is law-of-the-wall at top and bottom, + ! 2 is more KPP like. real :: TKE_decay ! The ratio of the natural Ekman depth to the TKE ! decay scale, nondimensional. real :: MKE_to_TKE_effic ! The efficiency with which mean kinetic energy @@ -450,7 +453,8 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & ! e.g. M=12 for DEPTH=4000m and DZ=1m real, dimension(SZK_(GV)+1) :: Vstar_Used, & ! 1D arrays used to store Mixing_Length_Used ! Vstar and Mixing_Length - + logical :: OLD_MLD_METHOD = .false. ! Can be set to compute MLD based on depth + ! TKE expires instead of using TKE threshold !/BGR - remaining variables are related to tracking iteration statistics. logical :: OBL_IT_STATS=.false. ! Flag for computing OBL iteration statistics REAL :: ITguess(20), ITresult(20),ITmax(20),ITmin(20) ! Flag for storing guess/result @@ -690,7 +694,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & !If prev value is present use for guess. MLD_guess=CS%ML_Depth2(i,j) else - !Otherwise guess middle of water column. + !Otherwise guess middle of water column (or stab_scale if smaller). MLD_guess = min(abs(stab_scale),0.5 * (min_MLD+max_MLD)) endif @@ -708,19 +712,20 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & sfc_connected(i) = .true. if (.not.CS%Use_Mstar_Fixed) then - ! Note this value of mech_TKE(i) now must be iterated over, so it is moved here - + ! Note the value of mech_TKE(i) now must be iterated over, so it is moved here ! First solve for the TKE to PE length scale - ! Change the (2.0) to a runtime parameter - MLD_over_Stab = 1.0 * MLD_guess / Stab_Scale + MLD_over_Stab = MLD_guess / Stab_Scale + !Fitting coefficients for hyperbolic decay as MLD -> Ekman depth MSTAR_A = CS%MSTAR_AT_XINT**(1./MSTAR_N) MSTAR_B = CS%MSTAR_SLOPE / (MSTAR_N*MSTAR_A**(MSTAR_N-1.)) if ((MLD_over_Stab) .ge. CS%MSTAR_XINT) then - MSTAR_mix = CS%MSTAR_SLOPE*(MLD_over_Stab-CS%MSTAR_XINT)+CS%MSTAR_AT_XINT + !Within linear regime + MSTAR_mix = CS%MSTAR_SLOPE*(MLD_over_Stab-CS%MSTAR_XINT)+CS%MSTAR_AT_XINT else - MSTAR_mix = (MSTAR_B*(MLD_over_Stab-CS%MSTAR_XINT)+MSTAR_A)**(MSTAR_N) + !Within hypterbolic decay regime + MSTAR_mix = (MSTAR_B*(MLD_over_Stab-CS%MSTAR_XINT)+MSTAR_A)**(MSTAR_N) endif - !print*,MLD_guess,STAB_SCALE,MSTAR_MIX + !Reset mech_tke and conv_perel values (based on new mstar) mech_TKE(i) = (dt*MSTAR_mix*GV%Rho0)*((U_Star**3)) conv_PErel(i) = 0.0 if (CS%TKE_diagnostics) then @@ -742,11 +747,6 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & conv_PErel(i) = conv_PErel(i) + TKE_forced(i,j,1) endif else - ! Multiply mech TKE at surface by enhancement (w'/ust) cubed. - ! This has not been confirmed to reproduce LES Langmuir simulations, - ! but is consistent if we assume we are trying to capture the enhancement - ! of the TKE flux/ (turbulent velocity scale cubed) knowing that the - ! enhancement factor is that of the turbulent velocity scale. mech_TKE(i) = mech_TKE_top(i)*ENHANCE_V**(3.) ; conv_PErel(i) = conv_PErel_top(i) endif @@ -779,7 +779,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & do K=2,nz+1 h_rsum = h_rsum + h(i,k-1) MixLen_shape(K) = CS%transLay_scale + (1.0 - CS%transLay_scale) * & - (max(0.0, (MLD_guess - h_rsum)*I_MLD) )**1.0 !BR prefer 1 or 2 for exponent? + (max(0.0, (MLD_guess - h_rsum)*I_MLD) )**CS%MixLenExponent enddo endif @@ -1265,64 +1265,66 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & mixing_debug = dPE_debug * IdtdR0 endif k = nz ! This is here to allow a breakpoint to be set. - - !/BGR: The following lines are used for the iteration + !/BGR: The following lines are used for the iteration + ! note the iteration has been altered to use the value predicted by + ! the TKE threshold (ML_DEPTH). This is because the MSTAR + ! is now dependent on the ML, and therefore the ML needs to be estimated + ! more precisely than the grid spacing. ITmax(obl_it) = max_MLD ! Track max } ITmin(obl_it) = min_MLD ! Track min } For debug purpose ITguess(obl_it) = MLD_guess ! Track guess } - MLD_FOUND=0.0 ; FIRST_OBL=.true. - ! MLD_FOUND=CS%ML_depth2(i,J) - do k=2,nz - if (FIRST_OBL) then !Breaks when OBL found - if (Vstar_Used(k) > 1.e-10 .and. k < nz) then - !If within OBL, keep integrating to find OBL - !/ Add thickness of level above to OBL depth. - MLD_FOUND = MLD_FOUND+h(i,k-1)*GV%H_to_m - !1. Check if guess was too shallow - !Adding -1 m as cushion, will help avoid - ! non-convergence flag when nearly converged. - else - FIRST_OBL = .false. - endif - endif - enddo - MLD_FOUND=CS%ML_DEPTH(i,j) - if (MLD_FOUND-CS%MLD_tol > MLD_guess) then - !elseif (MLD_FOUND > MLD_guess) then - !/ Guess was too shallow, set new minimum guess - min_MLD = MLD_guess - !2. Check if Guess minus found MLD - ! is less than thickness of level (=converged) - ! - We could try to add a more precise - ! value for found MLD, but seems difficult to - ! to contrain beyond within a level. - elseif ((MLD_guess-MLD_FOUND) < max(CS%MLD_tol,h(i,k-1)*GV%H_to_m)) then - ! Converged. Exit iteration. - OBL_CONVERGED = .true.!Break convergence loop - ! Testing Output, comment for use. - !print*,'Converged--------' - !print*,MLD_FOUND,MLD_guess - !/ - if (OBL_IT_STATS) then !Compute iteration statistics + if (OLD_MLD_METHOD) then + MLD_FOUND=0.0 ; FIRST_OBL=.true. + MLD_FOUND=CS%ML_depth2(i,J) + do k=2,nz + if (FIRST_OBL) then !Breaks when OBL found + if (Vstar_Used(k) > 1.e-10 .and. k < nz) then + !If within OBL, keep integrating to find OBL + !/ Add thickness of level above to OBL depth. + MLD_FOUND = MLD_FOUND+h(i,k-1)*GV%H_to_m + else + FIRST_OBL = .false. + endif + endif + enddo + else + MLD_FOUND=CS%ML_DEPTH(i,j) + endif + if (MLD_FOUND-CS%MLD_tol > MLD_guess) then + !elseif (MLD_FOUND > MLD_guess) then + !/ Guess was too shallow, set new minimum guess + min_MLD = MLD_guess + !2. Check if Guess minus found MLD + ! is less than thickness of level (=converged) + ! - We could try to add a more precise + ! value for found MLD, but seems difficult to + ! to contrain beyond within a level. + elseif ((MLD_guess-MLD_FOUND) < max(CS%MLD_tol,h(i,k-1)*GV%H_to_m)) then + ! Converged. Exit iteration. + OBL_CONVERGED = .true.!Break convergence loop + ! Testing Output, comment for use. + !print*,'Converged--------' + !print*,MLD_FOUND,MLD_guess + !/ + if (OBL_IT_STATS) then !Compute iteration statistics MAXIT = max(MAXIT,obl_it) MINIT = min(MINIT,obl_it) SUMIT = SUMIT+obl_it NUMIT = NUMIT+1 - print*,MAXIT,MINIT,SUMIT/NUMIT - endif - !BGR this can be where MLD is stored for next time... - CS%ML_Depth2(i,j) = MLD_guess - !/ - !2. If not, guess was too deep - else - !Guess was too deep, set new maximum guess + !print*,MAXIT,MINIT,SUMIT/NUMIT + endif + !BGR this can be where MLD is stored for next time... + CS%ML_Depth2(i,j) = MLD_guess + !/ + !2. If not, guess was too deep + else + !Guess was too deep, set new maximum guess max_MLD = MLD_guess !We know this guess was too deep endif ! For next pass, guess average of minimum and maximum values. MLD_guess = MLD_FOUND!min_MLD*0.5 + max_MLD*0.5 ITresult(obl_it) = MLD_FOUND endif ; enddo ! Iteration loop for converged boundary layer thickness. - print*,mld_guess,CS%ml_depth(i,j) if (.not.OBL_CONVERGED) then !/Temp output, warn that EPBL didn't converge !/Print guess/found for every iteration step @@ -1814,6 +1816,10 @@ subroutine energetic_PBL_init(Time, G, GV, param_file, diag, CS) call get_param(param_file, mod, "MSTAR", CS%mstar, & "The ratio of the friction velocity cubed to the TKE \n"//& "input to the mixed layer.", "units=nondim", default=1.2) + call get_param(param_file, mod, "MIX_LEN_EXPONENT", CS%MixLenExponent, & + "The exponent applied to the ratio of the distance to the MLD \n"//& + "and the MLD depth which determines the shape of the mixing length.",& + "units=nondim", default=1.0) call get_param(param_file, mod, "MSTAR_SLOPE", CS%mstar_slope, & "The ratio of the friction velocity cubed to the TKE \n"//& "input to the mixed layer (used if MSTAR_FIXED=false).",& From f16cd209862b652f12205fae465174bfe01238b1 Mon Sep 17 00:00:00 2001 From: Brandon Reichl Date: Tue, 3 Jan 2017 18:12:33 -0500 Subject: [PATCH 4/9] Defaults fixed for input linear MSTAR relationship parameters --- src/parameterizations/vertical/MOM_energetic_PBL.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index 9abf0f2149..7e886573ff 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -1823,15 +1823,15 @@ subroutine energetic_PBL_init(Time, G, GV, param_file, diag, CS) call get_param(param_file, mod, "MSTAR_SLOPE", CS%mstar_slope, & "The ratio of the friction velocity cubed to the TKE \n"//& "input to the mixed layer (used if MSTAR_FIXED=false).",& - "units=nondim", default=1.2) + "units=nondim", default=1.0) call get_param(param_file, mod, "MSTAR_XINT", CS%mstar_xint, & "The ratio of the friction velocity cubed to the TKE \n"//& "input to the mixed layer (used if MSTAR_FIXED=false).",& - "units=nondim", default=1.2) + "units=nondim", default=-0.25) call get_param(param_file, mod, "MSTAR_AT_XINT", CS%mstar_at_xint, & "The ratio of the friction velocity cubed to the TKE \n"//& "input to the mixed layer (used if MSTAR_FIXED=false).",& - "units=nondim", default=1.2) + "units=nondim", default=0.15) call get_param(param_file, mod, "MSTAR_FIXED", CS%Use_Mstar_Fixed, & "True to use a fixed value of mstar, if false mstar depends \n"//& "on the Obhukov length and Ekman length.","units=nondim",& From 0458b7ca5c8954dbc64a7541b7a49cdac7a3cee1 Mon Sep 17 00:00:00 2001 From: Brandon Reichl Date: Wed, 4 Jan 2017 13:54:47 -0500 Subject: [PATCH 5/9] Changed ePBL defaults. --- .../vertical/MOM_energetic_PBL.F90 | 23 +++++++++++-------- 1 file changed, 13 insertions(+), 10 deletions(-) diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index 7e886573ff..80a23c7a84 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -158,6 +158,7 @@ module MOM_energetic_PBL logical :: Use_LT_LiFunction = .false. logical :: orig_PE_calc = .true. logical :: Use_MLD_iteration=.false. ! False to use old ePBL method. + logical :: Orig_MLD_iteration=.false. ! False to use old MLD value logical :: MLD_iteration_guess=.false. ! False to default to guessing half the ! ocean depth for the iteration. logical :: Mixing_Diagnostics = .false. ! Will be true when outputing mixing @@ -453,8 +454,6 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & ! e.g. M=12 for DEPTH=4000m and DZ=1m real, dimension(SZK_(GV)+1) :: Vstar_Used, & ! 1D arrays used to store Mixing_Length_Used ! Vstar and Mixing_Length - logical :: OLD_MLD_METHOD = .false. ! Can be set to compute MLD based on depth - ! TKE expires instead of using TKE threshold !/BGR - remaining variables are related to tracking iteration statistics. logical :: OBL_IT_STATS=.false. ! Flag for computing OBL iteration statistics REAL :: ITguess(20), ITresult(20),ITmax(20),ITmin(20) ! Flag for storing guess/result @@ -1273,7 +1272,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & ITmax(obl_it) = max_MLD ! Track max } ITmin(obl_it) = min_MLD ! Track min } For debug purpose ITguess(obl_it) = MLD_guess ! Track guess } - if (OLD_MLD_METHOD) then + if (CS%Orig_MLD_iteration) then MLD_FOUND=0.0 ; FIRST_OBL=.true. MLD_FOUND=CS%ML_depth2(i,J) do k=2,nz @@ -1821,20 +1820,20 @@ subroutine energetic_PBL_init(Time, G, GV, param_file, diag, CS) "and the MLD depth which determines the shape of the mixing length.",& "units=nondim", default=1.0) call get_param(param_file, mod, "MSTAR_SLOPE", CS%mstar_slope, & - "The ratio of the friction velocity cubed to the TKE \n"//& - "input to the mixed layer (used if MSTAR_FIXED=false).",& + "The slope of the linear relationship between mstar \n"//& + "and the length scale ratio (used if MSTAR_FIXED=false).",& "units=nondim", default=1.0) call get_param(param_file, mod, "MSTAR_XINT", CS%mstar_xint, & - "The ratio of the friction velocity cubed to the TKE \n"//& - "input to the mixed layer (used if MSTAR_FIXED=false).",& + "The value of the length scale ratio where the mstar \n"//& + "is linear above (used if MSTAR_FIXED=false).",& "units=nondim", default=-0.25) call get_param(param_file, mod, "MSTAR_AT_XINT", CS%mstar_at_xint, & - "The ratio of the friction velocity cubed to the TKE \n"//& - "input to the mixed layer (used if MSTAR_FIXED=false).",& + "The value of mstar at MSTAR_XINT \n"//& + "(used if MSTAR_FIXED=false).",& "units=nondim", default=0.15) call get_param(param_file, mod, "MSTAR_FIXED", CS%Use_Mstar_Fixed, & "True to use a fixed value of mstar, if false mstar depends \n"//& - "on the Obhukov length and Ekman length.","units=nondim",& + "on the composite Obhukov length and Ekman length.","units=nondim",& default=.true.) call get_param(param_file, mod, "NSTAR", CS%nstar, & "The portion of the buoyant potential energy imparted by \n"//& @@ -1888,6 +1887,10 @@ subroutine energetic_PBL_init(Time, G, GV, param_file, diag, CS) "A logical that specifies whether or not to use the \n"//& "distance to the bottom of the actively turblent boundary \n"//& "layer to help set the EPBL length scale.", default=.false.) + call get_param(param_file, mod, "ORIG_MLD_ITERATION", CS%ORIG_MLD_ITERATION, & + "A logical that specifies whether or not to use the \n"//& + "old method for determining MLD depth in iteration, which \n"//& + "is limited to resolution.", default=.true.) call get_param(param_file, mod, "MLD_ITERATION_GUESS", CS%MLD_ITERATION_GUESS, & "A logical that specifies whether or not to use the \n"//& "previous timestep MLD as a first guess in the MLD iteration.\n"//& From 24bfbdbc2aede5979f438842081cccec41500bc4 Mon Sep 17 00:00:00 2001 From: Brandon Reichl Date: Wed, 4 Jan 2017 14:55:57 -0500 Subject: [PATCH 6/9] Change default for Mixing Length exponent to 2.0 --- src/parameterizations/vertical/MOM_energetic_PBL.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index 80a23c7a84..f6dca489e5 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -1818,7 +1818,7 @@ subroutine energetic_PBL_init(Time, G, GV, param_file, diag, CS) call get_param(param_file, mod, "MIX_LEN_EXPONENT", CS%MixLenExponent, & "The exponent applied to the ratio of the distance to the MLD \n"//& "and the MLD depth which determines the shape of the mixing length.",& - "units=nondim", default=1.0) + "units=nondim", default=2.0) call get_param(param_file, mod, "MSTAR_SLOPE", CS%mstar_slope, & "The slope of the linear relationship between mstar \n"//& "and the length scale ratio (used if MSTAR_FIXED=false).",& From f314853e6b24be7dd4cdf9d54bab1307269256e0 Mon Sep 17 00:00:00 2001 From: Brandon Reichl Date: Wed, 4 Jan 2017 16:52:38 -0500 Subject: [PATCH 7/9] Restore orignal MLD guess and next MLD iteration guess to as-in previous ePBL version. --- src/parameterizations/vertical/MOM_energetic_PBL.F90 | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index f6dca489e5..aa95bda0a2 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -694,7 +694,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & MLD_guess=CS%ML_Depth2(i,j) else !Otherwise guess middle of water column (or stab_scale if smaller). - MLD_guess = min(abs(stab_scale),0.5 * (min_MLD+max_MLD)) + MLD_guess = 0.5 * (min_MLD+max_MLD) endif ! Iterate up to MAX_OBL_IT times to determine a converged EPBL depth. @@ -1321,7 +1321,11 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & max_MLD = MLD_guess !We know this guess was too deep endif ! For next pass, guess average of minimum and maximum values. - MLD_guess = MLD_FOUND!min_MLD*0.5 + max_MLD*0.5 + if (CS%Orig_MLD_iteration) then + MLD_guess = min_MLD*0.5 + max_MLD*0.5 + else + MLD_guess = MLD_Found + endif ITresult(obl_it) = MLD_FOUND endif ; enddo ! Iteration loop for converged boundary layer thickness. if (.not.OBL_CONVERGED) then From 9577d2744656bc53a55e81f168307ad310f403bb Mon Sep 17 00:00:00 2001 From: Brandon Reichl Date: Thu, 5 Jan 2017 12:51:56 -0500 Subject: [PATCH 8/9] Fix issues to reproduce old answers. --- .../vertical/MOM_energetic_PBL.F90 | 126 ++++++++++-------- 1 file changed, 67 insertions(+), 59 deletions(-) diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index aa95bda0a2..bc77093fbb 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -88,8 +88,8 @@ module MOM_energetic_PBL type, public :: energetic_PBL_CS ; private real :: mstar ! The ratio of the friction velocity cubed to the - ! TKE available to drive entrainment, nondimensional. - ! This quantity is the vertically integrated + ! TKE available to drive entrainment, nondimensional. + ! This quantity is the vertically integrated ! shear production minus the vertically integrated ! dissipation of TKE produced by shear. real :: nstar ! The fraction of the TKE input to the mixed layer @@ -134,23 +134,23 @@ module MOM_energetic_PBL real :: min_mix_len ! The minimum mixing length scale that will be ! used by ePBL, in m. The default (0) does not ! set a minimum. - real :: N2_Dissipation_Scale_Neg + real :: N2_Dissipation_Scale_Neg real :: N2_Dissipation_Scale_Pos ! A nondimensional scaling factor controlling the - ! loss of TKE due to enhanced dissipation in the presence - ! of stratification. This dissipation is applied to the - ! available TKE which includes both that generated at the + ! loss of TKE due to enhanced dissipation in the presence + ! of stratification. This dissipation is applied to the + ! available TKE which includes both that generated at the ! surface and that generated at depth. It may be important ! to distinguish which TKE flavor that this dissipation ! applies to in subsequent revisions of this code. - ! "_Neg" and "_Pos" refer to which scale is applied as a + ! "_Neg" and "_Pos" refer to which scale is applied as a ! function of negative or positive local buoyancy. real :: MSTAR_SLOPE ! Slope of the function which relates the shear production - ! to the mixing layer depth, Ekman depth, and Monin-Obukhov + ! to the mixing layer depth, Ekman depth, and Monin-Obukhov ! depth. real :: MSTAR_XINT ! Value where MSTAR function transitions from linear ! to decay toward MSTAR->0 at fully developed Ekman depth. - real :: MSTAR_AT_XINT ! Intercept value of MSTAR at value where function + real :: MSTAR_AT_XINT ! Intercept value of MSTAR at value where function ! changes to linear transition. type(time_type), pointer :: Time ! A pointer to the ocean model's clock. logical :: Use_Mstar_Fixed = .true. ! A logical to revert to a fixed m* @@ -188,7 +188,7 @@ module MOM_energetic_PBL integer :: id_TKE_mech_decay = -1, id_TKE_conv_decay = -1 integer :: id_Hsfc_used = -1 integer :: id_Mixing_Length = -1, id_Velocity_Scale = -1 - integer :: id_OSBL = -1, id_LT_Enhancement = -1 + integer :: id_OSBL = -1, id_LT_Enhancement = -1 end type energetic_PBL_CS integer :: num_msg = 0, max_msg = 2 @@ -358,7 +358,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & real :: U_star ! The surface friction velocity, in m s-1. real :: U_10 ! An approximate (neutral) wind speed backed out of the friction velocity real :: vstar ! An in-situ turbulent velocity, in m s-1. - real :: Enhance_V ! An enhancement factor for vstar, based here on Langmuir impact. + real :: Enhance_V ! An enhancement factor for vstar, based here on Langmuir impact. real :: hbs_here ! The local minimum of hb_hs and MixLen_shape, times a ! conversion factor from H to M, in m H-1. real :: nstar_FC ! The fraction of conv_PErel that can be converted to mixing, nondim. @@ -475,12 +475,12 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & real :: MLD_over_STAB ! Mixing layer depth divided by STAB_SCALE real :: MSTAR_N = -2. ! Exponent in hyperolic decay at negative values of MLD_over_STAB real :: MSTAR_A ! MSTAR_A and MSTAR_B are coefficients in hyperbolic decay which - real :: MSTAR_B ! are computed to match the value and slope of the linear fit at the + real :: MSTAR_B ! are computed to match the value and slope of the linear fit at the ! value of mld_over_stab where the transition is set. real :: MSTAR_MIX ! The value of mstar (Proportionality of TKE to drive mixing to ustar ! cubed) which is computed as a function of latitude, boundary layer depth, ! and the Monin-Obhukov depth. - type(cvmix_global_params_type) :: CVMIX_gravity + type(cvmix_global_params_type) :: CVMIX_gravity logical :: debug=.false. ! Change this hard-coded value for debugging. ! The following arrays are used only for debugging purposes. real :: dPE_debug, mixing_debug @@ -606,7 +606,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & if (U_Star < CS%ustar_min) U_Star = CS%ustar_min ! Computing u10 based on u_star and COARE 3.5 relationships call ust_2_u10_coare3p5(U_STAR*sqrt(GV%Rho0/1.225),U_10) - ! Computing B_Star, or the Buoyancy flux over the friction velocity. + ! Computing B_Star, or the Buoyancy flux over the friction velocity. B_Star = min(0.0,-buoy_Flux(i,j)/U_Star) ! Computing stability scale which correlates with TKE for mixing, where ! TKE for mixing = TKE production minus TKE dissipation @@ -745,7 +745,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & else conv_PErel(i) = conv_PErel(i) + TKE_forced(i,j,1) endif - else + else mech_TKE(i) = mech_TKE_top(i)*ENHANCE_V**(3.) ; conv_PErel(i) = conv_PErel_top(i) endif @@ -777,8 +777,13 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & MixLen_shape(1) = 1.0 do K=2,nz+1 h_rsum = h_rsum + h(i,k-1) - MixLen_shape(K) = CS%transLay_scale + (1.0 - CS%transLay_scale) * & - (max(0.0, (MLD_guess - h_rsum)*I_MLD) )**CS%MixLenExponent + if (CS%MixLenExponent==2.0)then + MixLen_shape(K) = CS%transLay_scale + (1.0 - CS%transLay_scale) * & + (max(0.0, (MLD_guess - h_rsum)*I_MLD) )**2!CS%MixLenExponent + else + MixLen_shape(K) = CS%transLay_scale + (1.0 - CS%transLay_scale) * & + (max(0.0, (MLD_guess - h_rsum)*I_MLD) )**CS%MixLenExponent + endif enddo endif @@ -1004,7 +1009,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & !Negative buoyancy (increases PE) N2_dissipation = 1.+CS%N2_DISSIPATION_SCALE_NEG else - !Positive buoyancy (decreases PE) + !Positive buoyancy (decreases PE) N2_dissipation = 1.+CS%N2_DISSIPATION_SCALE_POS endif @@ -1100,7 +1105,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & ! There is not enough energy to support the mixing, so reduce the ! diffusivity to what can be supported. Kddt_h_max = Kddt_h_g0 ; Kddt_h_min = 0.0 - TKE_left_max = tot_TKE + (MKE_src - N2_DISSIPATION*PE_chg_g0) ; + TKE_left_max = tot_TKE + (MKE_src - N2_DISSIPATION*PE_chg_g0) ; TKE_left_min = tot_TKE ! As a starting guess, take the minimum of a false position estimate @@ -1264,66 +1269,69 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & mixing_debug = dPE_debug * IdtdR0 endif k = nz ! This is here to allow a breakpoint to be set. - !/BGR: The following lines are used for the iteration + !/BGR + ! The following lines are used for the iteration ! note the iteration has been altered to use the value predicted by ! the TKE threshold (ML_DEPTH). This is because the MSTAR ! is now dependent on the ML, and therefore the ML needs to be estimated - ! more precisely than the grid spacing. + ! more precisely than the grid spacing. + !/ ITmax(obl_it) = max_MLD ! Track max } ITmin(obl_it) = min_MLD ! Track min } For debug purpose - ITguess(obl_it) = MLD_guess ! Track guess } + ITguess(obl_it) = MLD_guess ! Track guess } + !/ + MLD_FOUND=0.0 ; FIRST_OBL=.true. if (CS%Orig_MLD_iteration) then - MLD_FOUND=0.0 ; FIRST_OBL=.true. - MLD_FOUND=CS%ML_depth2(i,J) + !This is how the iteration was original conducted do k=2,nz if (FIRST_OBL) then !Breaks when OBL found if (Vstar_Used(k) > 1.e-10 .and. k < nz) then - !If within OBL, keep integrating to find OBL - !/ Add thickness of level above to OBL depth. MLD_FOUND = MLD_FOUND+h(i,k-1)*GV%H_to_m else FIRST_OBL = .false. + if (MLD_FOUND-CS%MLD_tol > MLD_guess) then + min_MLD = MLD_guess + elseif ((MLD_guess-MLD_FOUND) < max(CS%MLD_tol,h(i,k-1)*GV%H_to_m)) then + OBL_CONVERGED = .true.!Break convergence loop + if (OBL_IT_STATS) then !Compute iteration statistics + MAXIT = max(MAXIT,obl_it) + MINIT = min(MINIT,obl_it) + SUMIT = SUMIT+obl_it + NUMIT = NUMIT+1 + print*,MAXIT,MINIT,SUMIT/NUMIT + endif + CS%ML_Depth2(i,j) = MLD_guess + else + max_MLD = MLD_guess !We know this guess was too deep + endif endif endif enddo else + !New method uses ML_DEPTH as computed in ePBL routine MLD_FOUND=CS%ML_DEPTH(i,j) - endif - if (MLD_FOUND-CS%MLD_tol > MLD_guess) then - !elseif (MLD_FOUND > MLD_guess) then - !/ Guess was too shallow, set new minimum guess - min_MLD = MLD_guess - !2. Check if Guess minus found MLD - ! is less than thickness of level (=converged) - ! - We could try to add a more precise - ! value for found MLD, but seems difficult to - ! to contrain beyond within a level. - elseif ((MLD_guess-MLD_FOUND) < max(CS%MLD_tol,h(i,k-1)*GV%H_to_m)) then - ! Converged. Exit iteration. - OBL_CONVERGED = .true.!Break convergence loop - ! Testing Output, comment for use. - !print*,'Converged--------' - !print*,MLD_FOUND,MLD_guess - !/ - if (OBL_IT_STATS) then !Compute iteration statistics - MAXIT = max(MAXIT,obl_it) - MINIT = min(MINIT,obl_it) - SUMIT = SUMIT+obl_it - NUMIT = NUMIT+1 - !print*,MAXIT,MINIT,SUMIT/NUMIT + if (MLD_FOUND-CS%MLD_tol > MLD_guess) then + min_MLD = MLD_guess + elseif ((MLD_guess-MLD_FOUND) < (CS%MLD_tol)) then + OBL_CONVERGED = .true.!Break convergence loop + if (OBL_IT_STATS) then !Compute iteration statistics + MAXIT = max(MAXIT,obl_it) + MINIT = min(MINIT,obl_it) + SUMIT = SUMIT+obl_it + NUMIT = NUMIT+1 + print*,MAXIT,MINIT,SUMIT/NUMIT + endif + CS%ML_Depth2(i,j) = MLD_guess + else + max_MLD = MLD_guess !We know this guess was too deep endif - !BGR this can be where MLD is stored for next time... - CS%ML_Depth2(i,j) = MLD_guess - !/ - !2. If not, guess was too deep - else - !Guess was too deep, set new maximum guess - max_MLD = MLD_guess !We know this guess was too deep endif ! For next pass, guess average of minimum and maximum values. if (CS%Orig_MLD_iteration) then MLD_guess = min_MLD*0.5 + max_MLD*0.5 else + !No longer use binary search, rather use the value found + !as the next guess in new method. MLD_guess = MLD_Found endif ITresult(obl_it) = MLD_FOUND @@ -1748,7 +1756,7 @@ subroutine ust_2_u10_coare3p5(USTair,U10) real :: z0sm, z0, z0rough, u10a, alpha, CD integer :: CT - ! Uses empirical formula for z0 to convert ustar_air to u10 based on the + ! Uses empirical formula for z0 to convert ustar_air to u10 based on the ! COARE 3.5 paper (Edson et al., 2013) !alpha=m*U10+b !Note in Edson et al. 2013, eq. 13 m is given as 0.017. However, @@ -1773,8 +1781,8 @@ subroutine ust_2_u10_coare3p5(USTair,U10) ! and converged rapidly (e.g. 2 cycles) ! for ustar=0.0001:0.0001:10. if (CT>20) then - u10 = USTair/sqrt(0.0015) ! I don't expect to get here, but just - ! in case it will output a reasonable value. + u10 = USTair/sqrt(0.0015) ! I don't expect to get here, but just + ! in case it will output a reasonable value. exit endif enddo From 5545872bf524b876985b3f7100631a1436ecac1e Mon Sep 17 00:00:00 2001 From: Brandon Reichl Date: Thu, 5 Jan 2017 16:54:09 -0500 Subject: [PATCH 9/9] Tweaked issues with iteration and changed defaults as appropriate. --- .../vertical/MOM_energetic_PBL.F90 | 20 +++++++------------ 1 file changed, 7 insertions(+), 13 deletions(-) diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index bc77093fbb..2cb16ae789 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -608,9 +608,6 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & call ust_2_u10_coare3p5(U_STAR*sqrt(GV%Rho0/1.225),U_10) ! Computing B_Star, or the Buoyancy flux over the friction velocity. B_Star = min(0.0,-buoy_Flux(i,j)/U_Star) - ! Computing stability scale which correlates with TKE for mixing, where - ! TKE for mixing = TKE production minus TKE dissipation - Stab_Scale = -u_star**2 / ( VonKar * ( C_MO * B_STAR + C_EK * u_star * absf(i))) if (CS%omega_frac >= 1.0) then ; absf(i) = 2.0*CS%omega else absf(i) = 0.25*((abs(G%CoriolisBu(I,J)) + abs(G%CoriolisBu(I-1,J-1))) + & @@ -618,7 +615,10 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & if (CS%omega_frac > 0.0) & absf(i) = sqrt(CS%omega_frac*4.0*CS%omega**2 + (1.0-CS%omega_frac)*absf(i)**2) endif - + ! Computing stability scale which correlates with TKE for mixing, where + ! TKE for mixing = TKE production minus TKE dissipation + Stab_Scale = -u_star**2 / ( VonKar * ( C_MO * B_STAR + C_EK * u_star * absf(i))) + if (CS%Use_Mstar_Fixed) then mech_TKE(i) = (dt*CS%mstar*GV%Rho0)*((U_Star**3)) conv_PErel(i) = 0.0 @@ -1312,7 +1312,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & MLD_FOUND=CS%ML_DEPTH(i,j) if (MLD_FOUND-CS%MLD_tol > MLD_guess) then min_MLD = MLD_guess - elseif ((MLD_guess-MLD_FOUND) < (CS%MLD_tol)) then + elseif (abs(MLD_guess-MLD_FOUND) < (CS%MLD_tol)) then OBL_CONVERGED = .true.!Break convergence loop if (OBL_IT_STATS) then !Compute iteration statistics MAXIT = max(MAXIT,obl_it) @@ -1327,13 +1327,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, CS, & endif endif ! For next pass, guess average of minimum and maximum values. - if (CS%Orig_MLD_iteration) then - MLD_guess = min_MLD*0.5 + max_MLD*0.5 - else - !No longer use binary search, rather use the value found - !as the next guess in new method. - MLD_guess = MLD_Found - endif + MLD_guess = min_MLD*0.5 + max_MLD*0.5 ITresult(obl_it) = MLD_FOUND endif ; enddo ! Iteration loop for converged boundary layer thickness. if (.not.OBL_CONVERGED) then @@ -1842,7 +1836,7 @@ subroutine energetic_PBL_init(Time, G, GV, param_file, diag, CS) call get_param(param_file, mod, "MSTAR_AT_XINT", CS%mstar_at_xint, & "The value of mstar at MSTAR_XINT \n"//& "(used if MSTAR_FIXED=false).",& - "units=nondim", default=0.15) + "units=nondim", default=0.13) call get_param(param_file, mod, "MSTAR_FIXED", CS%Use_Mstar_Fixed, & "True to use a fixed value of mstar, if false mstar depends \n"//& "on the composite Obhukov length and Ekman length.","units=nondim",&