diff --git a/config_src/drivers/solo_driver/MOM_surface_forcing.F90 b/config_src/drivers/solo_driver/MOM_surface_forcing.F90 index 10b5f377fa..c1e125be83 100644 --- a/config_src/drivers/solo_driver/MOM_surface_forcing.F90 +++ b/config_src/drivers/solo_driver/MOM_surface_forcing.F90 @@ -106,7 +106,7 @@ module MOM_surface_forcing real :: gyres_taux_cos_amp !< The amplitude of cosine wind stress gyres [R L Z T-1 ~> Pa], if WIND_CONFIG=='gyres' real :: gyres_taux_n_pis !< The number of sine lobes in the basin if WIND_CONFIG=='gyres' integer :: answer_date !< This 8-digit integer gives the approximate date with which the order - !! of arithmetic and and expressions were added to the code. + !! of arithmetic and expressions were added to the code. !! Dates before 20190101 use original answers. !! Dates after 20190101 use a form of the gyre wind stresses that are !! rotationally invariant and more likely to be the same between compilers. @@ -161,8 +161,8 @@ module MOM_surface_forcing character(len=200) :: salinityrestore_file = '' !< The file from which to read the sea surface !! salinity to restore toward - character(len=80) :: stress_x_var = '' !< X-windstress variable name in the input file - character(len=80) :: stress_y_var = '' !< Y-windstress variable name in the input file + character(len=80) :: stress_x_var = '' !< X-wind stress variable name in the input file + character(len=80) :: stress_y_var = '' !< Y-wind stress variable name in the input file character(len=80) :: ustar_var = '' !< ustar variable name in the input file character(len=80) :: LW_var = '' !< longwave heat flux variable name in the input file character(len=80) :: SW_var = '' !< shortwave heat flux variable name in the input file @@ -447,6 +447,8 @@ subroutine wind_forcing_2gyre(sfc_state, forces, day, G, US, CS) forces%tauy(i,J) = 0.0 enddo ; enddo + if (associated(forces%ustar)) call stresses_to_ustar(forces, G, US, CS) + call callTree_leave("wind_forcing_2gyre") end subroutine wind_forcing_2gyre @@ -484,6 +486,8 @@ subroutine wind_forcing_1gyre(sfc_state, forces, day, G, US, CS) forces%tauy(i,J) = 0.0 enddo ; enddo + if (associated(forces%ustar)) call stresses_to_ustar(forces, G, US, CS) + call callTree_leave("wind_forcing_1gyre") end subroutine wind_forcing_1gyre @@ -499,8 +503,6 @@ subroutine wind_forcing_gyres(sfc_state, forces, day, G, US, CS) !! a previous surface_forcing_init call ! Local variables real :: PI ! A common irrational number, 3.1415926535... [nondim] - real :: I_rho ! The inverse of the reference density times a ratio of scaling - ! factors [Z L-1 R-1 ~> m3 kg-1] real :: y ! The latitude relative to the south normalized by the domain extent [nondim] integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq @@ -530,12 +532,7 @@ subroutine wind_forcing_gyres(sfc_state, forces, day, G, US, CS) forces%taux(i-1,j)*forces%taux(i-1,j) + forces%taux(i,j)*forces%taux(i,j)))/CS%Rho0) ) enddo ; enddo else - I_rho = US%L_to_Z / CS%Rho0 - do j=js,je ; do i=is,ie - forces%ustar(i,j) = sqrt( (CS%gust_const + & - sqrt(0.5*((forces%tauy(i,J-1)**2 + forces%tauy(i,J)**2) + & - (forces%taux(I-1,j)**2 + forces%taux(I,j)**2))) ) * I_rho ) - enddo ; enddo + call stresses_to_ustar(forces, G, US, CS) endif call callTree_leave("wind_forcing_gyres") @@ -558,8 +555,6 @@ subroutine Neverworld_wind_forcing(sfc_state, forces, day, G, US, CS) real :: Pa_to_RLZ_T2 ! A unit conversion factor from Pa to the internal units ! for wind stresses [R Z L T-2 Pa-1 ~> 1] real :: PI ! A common irrational number, 3.1415926535... [nondim] - real :: I_rho ! The inverse of the reference density times a ratio of scaling - ! factors [Z L-1 R-1 ~> m3 kg-1] real :: y ! The latitude relative to the south normalized by the domain extent [nondim] real :: tau_max ! The magnitude of the wind stress [R Z L T-2 ~> Pa] real :: off ! An offset in the relative latitude [nondim] @@ -602,14 +597,7 @@ subroutine Neverworld_wind_forcing(sfc_state, forces, day, G, US, CS) enddo ; enddo ! Set the surface friction velocity, in units of [Z T-1 ~> m s-1]. ustar is always positive. - if (associated(forces%ustar)) then - I_rho = US%L_to_Z / CS%Rho0 - do j=js,je ; do i=is,ie - forces%ustar(i,j) = sqrt( (CS%gust_const + & - sqrt(0.5*((forces%tauy(i,J-1)**2 + forces%tauy(i,J)**2) + & - (forces%taux(I-1,j)**2 + forces%taux(I,j)**2))) ) * I_rho ) - enddo ; enddo - endif + if (associated(forces%ustar)) call stresses_to_ustar(forces, G, US, CS) end subroutine Neverworld_wind_forcing @@ -625,8 +613,6 @@ subroutine scurve_wind_forcing(sfc_state, forces, day, G, US, CS) !! a previous surface_forcing_init call ! Local variables integer :: i, j, kseg - real :: I_rho ! The inverse of the reference density times a ratio of scaling - ! factors [Z L-1 R-1 ~> m3 kg-1] real :: y_curve ! The latitude relative to the southern end of a curve segment [degreesN] real :: L_curve ! The latitudinal extent of a curve segment [degreesN] ! real :: ydata(7) = (/ -70., -45., -15., 0., 15., 45., 70. /) @@ -657,14 +643,7 @@ subroutine scurve_wind_forcing(sfc_state, forces, day, G, US, CS) enddo ; enddo ! Set the surface friction velocity, in units of [Z T-1 ~> m s-1]. ustar is always positive. - if (associated(forces%ustar)) then - I_rho = US%L_to_Z / CS%Rho0 - do j=G%jsc,G%jec ; do i=G%isc,G%iec - forces%ustar(i,j) = sqrt( (CS%gust_const + & - sqrt(0.5*((forces%tauy(i,J-1)**2 + forces%tauy(i,J)**2) + & - (forces%taux(I-1,j)**2 + forces%taux(I,j)**2))) ) * I_rho ) - enddo ; enddo - endif + if (associated(forces%ustar)) call stresses_to_ustar(forces, G, US, CS) end subroutine scurve_wind_forcing @@ -892,6 +871,37 @@ subroutine wind_forcing_by_data_override(sfc_state, forces, day, G, US, CS) call callTree_leave("wind_forcing_by_data_override") end subroutine wind_forcing_by_data_override +!> Translate the wind stresses into the friction velocity, including effects of background gustiness. +subroutine stresses_to_ustar(forces, G, US, CS) + type(mech_forcing), intent(inout) :: forces !< A structure with the driving mechanical forces + type(ocean_grid_type), intent(in) :: G !< Grid structure. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(surface_forcing_CS), pointer :: CS !< pointer to control structure returned by + !! a previous surface_forcing_init call + ! 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] + integer :: i, j, is, ie, js, je + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + + I_rho = US%L_to_Z / CS%Rho0 + + if (CS%read_gust_2d) then + do j=js,je ; do i=is,ie + forces%ustar(i,j) = sqrt( (CS%gust(i,j) + & + sqrt(0.5*((forces%tauy(i,j-1)**2 + forces%tauy(i,j)**2) + & + (forces%taux(i-1,j)**2 + forces%taux(i,j)**2))) ) * I_rho ) + enddo ; enddo + else + do j=js,je ; do i=is,ie + forces%ustar(i,j) = sqrt( (CS%gust_const + & + sqrt(0.5*((forces%tauy(i,J-1)**2 + forces%tauy(i,J)**2) + & + (forces%taux(I-1,j)**2 + forces%taux(I,j)**2))) ) * I_rho ) + enddo ; enddo + endif + +end subroutine stresses_to_ustar !> Specifies zero surface buoyancy fluxes from input files. subroutine buoyancy_forcing_from_files(sfc_state, fluxes, day, dt, G, US, CS) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index f6c714b3be..267a162b00 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1982,10 +1982,9 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & default=.false.) CS%tv%T_is_conT = use_conT_absS ; CS%tv%S_is_absS = use_conT_absS call get_param(param_file, "MOM", "ADIABATIC", CS%adiabatic, & - "There are no diapycnal mass fluxes if ADIABATIC is "//& - "true. This assumes that KD = KDML = 0.0 and that "//& - "there is no buoyancy forcing, but makes the model "//& - "faster by eliminating subroutine calls.", default=.false.) + "There are no diapycnal mass fluxes if ADIABATIC is true. "//& + "This assumes that KD = 0.0 and that there is no buoyancy forcing, "//& + "but makes the model faster by eliminating subroutine calls.", default=.false.) call get_param(param_file, "MOM", "DO_DYNAMICS", CS%do_dynamics, & "If False, skips the dynamics calls that update u & v, as well as "//& "the gravity wave adjustment to h. This may be a fragile feature, "//& diff --git a/src/core/MOM_porous_barriers.F90 b/src/core/MOM_porous_barriers.F90 index cac8a44bf0..0e48cf07fd 100644 --- a/src/core/MOM_porous_barriers.F90 +++ b/src/core/MOM_porous_barriers.F90 @@ -284,7 +284,7 @@ subroutine calc_eta_at_uv(eta_u, eta_v, interp, dmask, h, tv, G, GV, US, eta_bt) real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: eta_bt !< optional barotropic variable !! used to dilate the layer thicknesses !! [H ~> m or kg m-2]. - real, intent(in) :: dmask !< The depth shaller than which + real, intent(in) :: dmask !< The depth shallower than which !! porous barrier is not applied [Z ~> m] integer, intent(in) :: interp !< eta interpolation method real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(out) :: eta_u !< Layer interface heights at u points [Z ~> m] @@ -456,13 +456,13 @@ subroutine porous_barriers_init(Time, US, param_file, diag, CS) ! The sign needs to be inverted to be consistent with the sign convention of Davg_[UV] CS%mask_depth = -CS%mask_depth call get_param(param_file, mdl, "PORBAR_ETA_INTERP", interp_method, & - "A string describing the method that decicdes how the "//& + "A string describing the method that decides how the "//& "interface heights at the velocity points are calculated. "//& "Valid values are:\n"//& "\t MAX (the default) - maximum of the adjacent cells \n"//& "\t MIN - minimum of the adjacent cells \n"//& "\t ARITHMETIC - arithmetic mean of the adjacent cells \n"//& - "\t HARMOINIC - harmonic mean of the adjacent cells \n", & + "\t HARMONIC - harmonic mean of the adjacent cells \n", & default=ETA_INTERP_MAX_STRING) select case (interp_method) case (ETA_INTERP_MAX_STRING) ; CS%eta_interp = ETA_INTERP_MAX diff --git a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 index 864669a217..a8efa4cf12 100644 --- a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 +++ b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 @@ -48,6 +48,7 @@ module MOM_mixed_layer_restrat logical :: MLE_use_PBL_MLD !< If true, use the MLD provided by the PBL parameterization. !! if false, MLE will calculate a MLD based on a density difference !! based on the parameter MLE_DENSITY_DIFF. + real :: vonKar !< The von Karman constant as used for mixed layer viscosity [nomdim] real :: MLE_MLD_decay_time !< Time-scale to use in a running-mean when MLD is retreating [T ~> s]. real :: MLE_MLD_decay_time2 !< Time-scale to use in a running-mean when filtered MLD is retreating [T ~> s]. real :: MLE_density_diff !< Density difference used in detecting mixed-layer depth [R ~> kg m-3]. @@ -189,6 +190,8 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var real :: dh ! Portion of the layer thickness that is in the mixed layer [H ~> m or kg m-2] real :: res_scaling_fac ! The resolution-dependent scaling factor [nondim] real :: I_LFront ! The inverse of the frontal length scale [L-1 ~> m-1] + real :: vonKar_x_pi2 ! A scaling constant that is approximately the von Karman constant times + ! pi squared [nondim] logical :: line_is_empty, keep_going, res_upscale integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz @@ -200,6 +203,8 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var covTS(:)=0.0 !!Functionality not implemented yet; in future, should be passed in tv varS(:)=0.0 + vonKar_x_pi2 = CS%vonKar * 9.8696 + if (.not.associated(tv%eqn_of_state)) call MOM_error(FATAL, "MOM_mixedlayer_restrat: "// & "An equation of state must be used with this module.") if (.not. allocated(VarMix%Rd_dx_h) .and. CS%front_length > 0.) & @@ -380,11 +385,10 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var ( sqrt( 0.5 * ( G%dxCu(I,j)**2 + G%dyCu(I,j)**2 ) ) * I_LFront ) & * min( 1., 0.5*( VarMix%Rd_dx_h(i,j) + VarMix%Rd_dx_h(i+1,j) ) ) - ! peak ML visc: u_star * 0.41 * (h_ml*u_star)/(absf*h_ml + 4.0*u_star) + ! peak ML visc: u_star * von_Karman * (h_ml*u_star)/(absf*h_ml + 4.0*u_star) ! momentum mixing rate: pi^2*visc/h_ml^2 - ! 0.41 is the von Karmen constant, 9.8696 = pi^2. h_vel = 0.5*((htot_fast(i,j) + htot_fast(i+1,j)) + h_neglect) * GV%H_to_Z - mom_mixrate = (0.41*9.8696)*u_star**2 / & + mom_mixrate = vonKar_x_pi2*u_star**2 / & (absf*h_vel**2 + 4.0*(h_vel+dz_neglect)*u_star) timescale = 0.0625 * (absf + 2.0*mom_mixrate) / (absf**2 + mom_mixrate**2) timescale = timescale * CS%ml_restrat_coef @@ -393,7 +397,7 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var (Rml_av_fast(i+1,j)-Rml_av_fast(i,j)) * (h_vel**2 * GV%Z_to_H) ! As above but using the slow filtered MLD h_vel = 0.5*((htot_slow(i,j) + htot_slow(i+1,j)) + h_neglect) * GV%H_to_Z - mom_mixrate = (0.41*9.8696)*u_star**2 / & + mom_mixrate = vonKar_x_pi2*u_star**2 / & (absf*h_vel**2 + 4.0*(h_vel+dz_neglect)*u_star) timescale = 0.0625 * (absf + 2.0*mom_mixrate) / (absf**2 + mom_mixrate**2) timescale = timescale * CS%ml_restrat_coef2 @@ -456,11 +460,10 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var ( sqrt( 0.5 * ( (G%dxCv(i,J))**2 + (G%dyCv(i,J))**2 ) ) * I_LFront ) & * min( 1., 0.5*( VarMix%Rd_dx_h(i,j) + VarMix%Rd_dx_h(i,j+1) ) ) - ! peak ML visc: u_star * 0.41 * (h_ml*u_star)/(absf*h_ml + 4.0*u_star) + ! peak ML visc: u_star * von_Karman * (h_ml*u_star)/(absf*h_ml + 4.0*u_star) ! momentum mixing rate: pi^2*visc/h_ml^2 - ! 0.41 is the von Karmen constant, 9.8696 = pi^2. h_vel = 0.5*((htot_fast(i,j) + htot_fast(i,j+1)) + h_neglect) * GV%H_to_Z - mom_mixrate = (0.41*9.8696)*u_star**2 / & + mom_mixrate = vonKar_x_pi2*u_star**2 / & (absf*h_vel**2 + 4.0*(h_vel+dz_neglect)*u_star) timescale = 0.0625 * (absf + 2.0*mom_mixrate) / (absf**2 + mom_mixrate**2) timescale = timescale * CS%ml_restrat_coef @@ -469,7 +472,7 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, Var (Rml_av_fast(i,j+1)-Rml_av_fast(i,j)) * (h_vel**2 * GV%Z_to_H) ! As above but using the slow filtered MLD h_vel = 0.5*((htot_slow(i,j) + htot_slow(i,j+1)) + h_neglect) * GV%H_to_Z - mom_mixrate = (0.41*9.8696)*u_star**2 / & + mom_mixrate = vonKar_x_pi2*u_star**2 / & (absf*h_vel**2 + 4.0*(h_vel+dz_neglect)*u_star) timescale = 0.0625 * (absf + 2.0*mom_mixrate) / (absf**2 + mom_mixrate**2) timescale = timescale * CS%ml_restrat_coef2 @@ -617,6 +620,8 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS) real :: h_vel ! htot interpolated onto velocity points [Z ~> m]. (The units are not H.) real :: absf ! absolute value of f, interpolated to velocity points [T-1 ~> s-1] real :: u_star ! surface friction velocity, interpolated to velocity points [Z T-1 ~> m s-1]. + real :: vonKar_x_pi2 ! A scaling constant that is approximately the von Karman constant times + ! pi squared [nondim] real :: mom_mixrate ! rate at which momentum is homogenized within mixed layer [T-1 ~> s-1] real :: timescale ! mixing growth timescale [T ~> s] real :: h_min ! The minimum layer thickness [H ~> m or kg m-2]. h_min could be 0. @@ -653,6 +658,7 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS) uDml(:) = 0.0 ; vDml(:) = 0.0 I4dt = 0.25 / dt g_Rho0 = GV%g_Earth / GV%Rho0 + vonKar_x_pi2 = CS%vonKar * 9.8696 use_EOS = associated(tv%eqn_of_state) h_neglect = GV%H_subroundoff dz_neglect = GV%H_subroundoff*GV%H_to_Z @@ -701,10 +707,9 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS) u_star = 0.5*(forces%ustar(i,j) + forces%ustar(i+1,j)) absf = 0.5*(abs(G%CoriolisBu(I,J-1)) + abs(G%CoriolisBu(I,J))) - ! peak ML visc: u_star * 0.41 * (h_ml*u_star)/(absf*h_ml + 4.0*u_star) + ! peak ML visc: u_star * von_Karman * (h_ml*u_star)/(absf*h_ml + 4.0*u_star) ! momentum mixing rate: pi^2*visc/h_ml^2 - ! 0.41 is the von Karmen constant, 9.8696 = pi^2. - mom_mixrate = (0.41*9.8696)*u_star**2 / & + mom_mixrate = vonKar_x_pi2*u_star**2 / & (absf*h_vel**2 + 4.0*(h_vel+dz_neglect)*u_star) timescale = 0.0625 * (absf + 2.0*mom_mixrate) / (absf**2 + mom_mixrate**2) @@ -748,10 +753,9 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS) u_star = 0.5*(forces%ustar(i,j) + forces%ustar(i,j+1)) absf = 0.5*(abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J))) - ! peak ML visc: u_star * 0.41 * (h_ml*u_star)/(absf*h_ml + 4.0*u_star) + ! peak ML visc: u_star * von_Karman * (h_ml*u_star)/(absf*h_ml + 4.0*u_star) ! momentum mixing rate: pi^2*visc/h_ml^2 - ! 0.41 is the von Karmen constant, 9.8696 = pi^2. - mom_mixrate = (0.41*9.8696)*u_star**2 / & + mom_mixrate = vonKar_x_pi2*u_star**2 / & (absf*h_vel**2 + 4.0*(h_vel+dz_neglect)*u_star) timescale = 0.0625 * (absf + 2.0*mom_mixrate) / (absf**2 + mom_mixrate**2) @@ -877,6 +881,9 @@ logical function mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, CS, call get_param(param_file, mdl, "USE_STANLEY_ML", CS%use_stanley_ml, & "If true, turn on Stanley SGS T variance parameterization "// & "in ML restrat code.", default=.false.) + call get_param(param_file, mdl, 'VON_KARMAN_CONST', CS%vonKar, & + 'The value the von Karman constant as used for mixed layer viscosity.', & + units='nondim', default=0.41) ! We use GV%nkml to distinguish between the old and new implementation of MLE. ! The old implementation only works for the layer model with nkml>0. if (GV%nkml==0) then diff --git a/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 b/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 index 7c427ab79a..9eabd8dc4d 100644 --- a/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 +++ b/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 @@ -9,7 +9,7 @@ module MOM_bkgnd_mixing use MOM_diag_mediator, only : diag_ctrl, time_type, register_diag_field use MOM_diag_mediator, only : post_data use MOM_error_handler, only : MOM_error, FATAL, WARNING, NOTE -use MOM_file_parser, only : get_param, log_version, param_file_type +use MOM_file_parser, only : get_param, log_param, log_version, param_file_type use MOM_file_parser, only : openParameterBlock, closeParameterBlock use MOM_forcing_type, only : forcing use MOM_grid, only : ocean_grid_type @@ -63,9 +63,11 @@ module MOM_bkgnd_mixing real :: Kd_tanh_lat_scale !< A nondimensional scaling for the range of !! diffusivities with Kd_tanh_lat_fn. Valid values !! are in the range of -2 to 2; 0.4 reproduces CM2M. - real :: Kdml !< mixed layer diapycnal diffusivity [Z2 T-1 ~> m2 s-1] - !! when bulkmixedlayer==.false. - real :: Hmix !< mixed layer thickness [Z ~> m] when bulkmixedlayer==.false. + real :: Kd_tot_ml !< The mixed layer diapycnal diffusivity [Z2 T-1 ~> m2 s-1] + !! when no other physically based mixed layer turbulence + !! parameterization is being used. + real :: Hmix !< mixed layer thickness [Z ~> m] when no physically based + !! ocean surface boundary layer parameterization is used. logical :: Kd_tanh_lat_fn !< If true, use the tanh dependence of Kd_sfc on !! latitude, like GFDL CM2.1/CM2M. There is no !! physical justification for this form, and it can @@ -92,7 +94,8 @@ module MOM_bkgnd_mixing !! where kd is the diapycnal diffusivity. This approach assumes that work done !! against gravity is uniformly distributed throughout the column. Whereas, kd=kd_0*e, !! as in the original version, concentrates buoyancy work in regions of strong stratification. - logical :: bulkmixedlayer !< If true, a refined bulk mixed layer scheme is used + logical :: physical_OBL_scheme !< If true, a physically-based scheme is used to determine mixing in the + !! ocean's surface boundary layer, such as ePBL, KPP, or a refined bulk mixed layer scheme. logical :: Kd_via_Kdml_bug !< If true and KDML /= KD and a number of other higher precedence !! options are not used, the background diffusivity is set incorrectly using a !! bug that was introduced in March, 2018. @@ -109,7 +112,7 @@ module MOM_bkgnd_mixing contains !> Initialize the background mixing routine. -subroutine bkgnd_mixing_init(Time, G, GV, US, param_file, diag, CS) +subroutine bkgnd_mixing_init(Time, G, GV, US, param_file, diag, CS, physical_OBL_scheme) type(time_type), intent(in) :: Time !< The current time. type(ocean_grid_type), intent(in) :: G !< Grid structure. @@ -117,15 +120,20 @@ subroutine bkgnd_mixing_init(Time, G, GV, US, param_file, diag, CS) type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(param_file_type), intent(in) :: param_file !< Run-time parameter file handle type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure. - type(bkgnd_mixing_cs), pointer :: CS !< This module's control structure. + type(bkgnd_mixing_cs), pointer :: CS !< This module's control structure. + logical, intent(in) :: physical_OBL_scheme !< If true, a physically based + !! parameterization (like KPP or ePBL or a bulk mixed + !! layer) is used outside of set_diffusivity to + !! specify the mixing that occurs in the ocean's + !! surface boundary layer. ! Local variables real :: Kv ! The interior vertical viscosity [Z2 T-1 ~> m2 s-1] - read to set Prandtl ! number unless it is provided as a parameter real :: prandtl_bkgnd_comp ! Kv/CS%Kd. Gets compared with user-specified prandtl_bkgnd. -! This include declares and sets the variable "version". -#include "version_variable.h" + ! This include declares and sets the variable "version". +# include "version_variable.h" if (associated(CS)) then call MOM_error(WARNING, "bkgnd_mixing_init called with an associated "// & @@ -154,21 +162,38 @@ subroutine bkgnd_mixing_init(Time, G, GV, US, param_file, diag, CS) ! The following is needed to set one of the choices of vertical background mixing - ! BULKMIXEDLAYER is not always defined (e.g., CM2G63L), so the following line by passes - ! the need to include BULKMIXEDLAYER in MOM_input - CS%bulkmixedlayer = (GV%nkml > 0) - if (CS%bulkmixedlayer) then + CS%physical_OBL_scheme = physical_OBL_scheme + if (CS%physical_OBL_scheme) then ! Check that Kdml is not set when using bulk mixed layer - call get_param(param_file, mdl, "KDML", CS%Kdml, default=-1.) - if (CS%Kdml>0.) call MOM_error(FATAL, & - "bkgnd_mixing_init: KDML cannot be set when using bulk mixed layer.") - CS%Kdml = CS%Kd ! This is not used with a bulk mixed layer, but also cannot be a NaN. + call get_param(param_file, mdl, "KDML", CS%Kd_tot_ml, default=-1., do_not_log=.true.) + if (CS%Kd_tot_ml>0.) call MOM_error(FATAL, & + "bkgnd_mixing_init: KDML is a depricated parameter that should not be used.") + call get_param(param_file, mdl, "KD_ML_TOT", CS%Kd_tot_ml, default=-1., do_not_log=.true.) + if (CS%Kd_tot_ml>0.) call MOM_error(FATAL, & + "bkgnd_mixing_init: KD_ML_TOT cannot be set when using a physically based ocean "//& + "boundary layer mixing parameterization.") + CS%Kd_tot_ml = CS%Kd ! This is not used with a bulk mixed layer, but also cannot be a NaN. else - call get_param(param_file, mdl, "KDML", CS%Kdml, & + call get_param(param_file, mdl, "KD_ML_TOT", CS%Kd_tot_ml, & + "The total diapcynal diffusivity in the surface mixed layer when there is "//& + "not a physically based parameterization of mixing in the mixed layer, such "//& + "as bulk mixed layer or KPP or ePBL.", & + units="m2 s-1", default=CS%Kd*US%Z2_T_to_m2_s, scale=US%m2_s_to_Z2_T, do_not_log=.true.) + if (abs(CS%Kd_tot_ml - CS%Kd) <= 1.0e-15*abs(CS%Kd)) then + call get_param(param_file, mdl, "KDML", CS%Kd_tot_ml, & "If BULKMIXEDLAYER is false, KDML is the elevated "//& "diapycnal diffusivity in the topmost HMIX of fluid. "//& "KDML is only used if BULKMIXEDLAYER is false.", & - units="m2 s-1", default=CS%Kd*US%Z2_T_to_m2_s, scale=US%m2_s_to_Z2_T) + units="m2 s-1", default=CS%Kd*US%Z2_T_to_m2_s, scale=US%m2_s_to_Z2_T, do_not_log=.true.) + if (abs(CS%Kd_tot_ml - CS%Kd) > 1.0e-15*abs(CS%Kd)) & + call MOM_error(WARNING, "KDML is a depricated parameter. Use KD_ML_TOT instead.") + endif + call log_param(param_file, mdl, "KD_ML_TOT", CS%Kd_tot_ml*US%Z2_T_to_m2_s, & + "The total diapcynal diffusivity in the surface mixed layer when there is "//& + "not a physically based parameterization of mixing in the mixed layer, such "//& + "as bulk mixed layer or KPP or ePBL.", & + units="m2 s-1", default=CS%Kd*US%Z2_T_to_m2_s) + call get_param(param_file, mdl, "HMIX_FIXED", CS%Hmix, & "The prescribed depth over which the near-surface "//& "viscosity and diffusivity are elevated when the bulk "//& @@ -290,13 +315,13 @@ subroutine bkgnd_mixing_init(Time, G, GV, US, param_file, diag, CS) "MOM_bkgnd_mixing: KD_TANH_LAT_FN can not be used with HENYEY_IGW_BACKGROUND.") CS%Kd_via_Kdml_bug = .false. - if ((CS%Kd /= CS%Kdml) .and. .not.(CS%Kd_tanh_lat_fn .or. CS%bulkmixedlayer .or. & + if ((CS%Kd /= CS%Kd_tot_ml) .and. .not.(CS%Kd_tanh_lat_fn .or. CS%physical_OBL_scheme .or. & CS%Henyey_IGW_background .or. CS%Henyey_IGW_background_new .or. & CS%horiz_varying_background .or. CS%Bryan_Lewis_diffusivity)) then call get_param(param_file, mdl, "KD_BACKGROUND_VIA_KDML_BUG", CS%Kd_via_Kdml_bug, & "If true and KDML /= KD and several other conditions apply, the background "//& "diffusivity is set incorrectly using a bug that was introduced in March, 2018.", & - default=.true.) ! The default should be changed to false and this parameter obsoleted. + default=.false.) ! The default should be changed to false and this parameter obsoleted. endif ! call closeParameterBlock(param_file) @@ -471,7 +496,7 @@ subroutine calculate_bkgnd_mixing(h, tv, N2_lay, Kd_lay, Kd_int, Kv_bkgnd, j, G, endif ! Now set background diffusivies based on these surface values, possibly with vertical structure. - if ((.not.CS%bulkmixedlayer) .and. (CS%Kd /= CS%Kdml)) then + if ((.not.CS%physical_OBL_scheme) .and. (CS%Kd /= CS%Kd_tot_ml)) then ! This is a crude way to put in a diffusive boundary layer without an explicit boundary ! layer turbulence scheme. It should not be used for any realistic ocean models. I_Hmix = 1.0 / (CS%Hmix + GV%H_subroundoff*GV%H_to_Z) @@ -481,16 +506,16 @@ subroutine calculate_bkgnd_mixing(h, tv, N2_lay, Kd_lay, Kd_int, Kv_bkgnd, j, G, if (CS%Kd_via_Kdml_bug) then ! These two lines should update Kd_lay, not Kd_int. They were correctly working on the ! same variables until MOM6 commit 7a818716 (PR#750), which was added on March 26, 2018. - if (depth_c <= CS%Hmix) then ; Kd_int(i,K) = CS%Kdml + if (depth_c <= CS%Hmix) then ; Kd_int(i,K) = CS%Kd_tot_ml elseif (depth_c >= 2.0*CS%Hmix) then ; Kd_int(i,K) = Kd_sfc(i) else - Kd_lay(i,k) = ((Kd_sfc(i) - CS%Kdml) * I_Hmix) * depth_c + (2.0*CS%Kdml - Kd_sfc(i)) + Kd_lay(i,k) = ((Kd_sfc(i) - CS%Kd_tot_ml) * I_Hmix) * depth_c + (2.0*CS%Kd_tot_ml - Kd_sfc(i)) endif else - if (depth_c <= CS%Hmix) then ; Kd_lay(i,k) = CS%Kdml + if (depth_c <= CS%Hmix) then ; Kd_lay(i,k) = CS%Kd_tot_ml elseif (depth_c >= 2.0*CS%Hmix) then ; Kd_lay(i,k) = Kd_sfc(i) else - Kd_lay(i,k) = ((Kd_sfc(i) - CS%Kdml) * I_Hmix) * depth_c + (2.0*CS%Kdml - Kd_sfc(i)) + Kd_lay(i,k) = ((Kd_sfc(i) - CS%Kd_tot_ml) * I_Hmix) * depth_c + (2.0*CS%Kd_tot_ml - Kd_sfc(i)) endif endif diff --git a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 index aa0d05ce79..49d62bbde4 100644 --- a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 +++ b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 @@ -49,6 +49,7 @@ module MOM_bulk_mixed_layer !! the mixed layer is converted to TKE [nondim]. real :: bulk_Ri_convective !< The efficiency with which convectively !! released mean kinetic energy becomes TKE [nondim]. + real :: vonKar !< The von Karman constant as used for mixed layer viscosity [nomdim] real :: Hmix_min !< The minimum mixed layer thickness [H ~> m or kg m-2]. real :: H_limit_fluxes !< When the total ocean depth is less than this !! value [H ~> m or kg m-2], scale away all surface forcing to @@ -316,7 +317,7 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, C real :: H_nbr ! A minimum thickness based on neighboring thicknesses [H ~> m or kg m-2]. real :: absf_x_H ! The absolute value of f times the mixed layer thickness [Z T-1 ~> m s-1]. - real :: kU_star ! Ustar times the Von Karmen constant [Z T-1 ~> m s-1]. + real :: kU_star ! Ustar times the Von Karman constant [Z T-1 ~> m s-1]. real :: dt__diag ! A recaled copy of dt_diag (if present) or dt [T ~> s]. logical :: write_diags ! If true, write out diagnostics with this step. logical :: reset_diags ! If true, zero out the accumulated diagnostics. @@ -618,12 +619,12 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, C ! as the third piece will then optimally describe mixed layer ! restratification. For nkml>=4 the whole strategy should be revisited. do i=is,ie - kU_star = 0.41*fluxes%ustar(i,j) ! Maybe could be replaced with u*+w*? + kU_star = CS%vonKar*fluxes%ustar(i,j) ! Maybe could be replaced with u*+w*? if (associated(fluxes%ustar_shelf) .and. & associated(fluxes%frac_shelf_h)) then if (fluxes%frac_shelf_h(i,j) > 0.0) & kU_star = (1.0 - fluxes%frac_shelf_h(i,j)) * kU_star + & - fluxes%frac_shelf_h(i,j) * (0.41*fluxes%ustar_shelf(i,j)) + fluxes%frac_shelf_h(i,j) * (CS%vonKar*fluxes%ustar_shelf(i,j)) endif absf_x_H = 0.25 * GV%H_to_Z * h(i,0) * & ((abs(G%CoriolisBu(I,J)) + abs(G%CoriolisBu(I-1,J-1))) + & @@ -1344,11 +1345,11 @@ subroutine find_starting_TKE(htot, h_CA, fluxes, Conv_En, cTKE, dKE_FC, dKE_CA, ! the equatorial areas. Although it is not cast as a parameter, it should ! be considered an empirical parameter, and it might depend strongly on the ! number of sublayers in the mixed layer and their locations. -! The 0.41 is VonKarman's constant. This equation assumes that small & large -! scales contribute to mixed layer deepening at similar rates, even though -! small scales are dissipated more rapidly (implying they are less efficient). -! Ih = 1.0/(16.0*0.41*U_star*dt) - Ih = GV%H_to_Z/(3.0*0.41*U_star*dt) +! This equation assumes that small & large scales contribute to mixed layer +! deepening at similar rates, even though small scales are dissipated more +! rapidly (implying they are less efficient). +! Ih = 1.0/(16.0*CS%vonKar*U_star*dt) + Ih = GV%H_to_Z/(3.0*CS%vonKar*U_star*dt) cMKE(1,i) = 4.0 * Ih ; cMKE(2,i) = (absf_Ustar*GV%H_to_Z) * Ih if (Idecay_len_TKE(i) > 0.0) then @@ -3387,6 +3388,9 @@ subroutine bulkmixedlayer_init(Time, G, GV, US, param_file, diag, CS) "kinetic energy is converted to turbulent kinetic "//& "energy. By default BULK_RI_CONVECTIVE=BULK_RI_ML.", & units="nondim", default=CS%bulk_Ri_ML) + call get_param(param_file, mdl, 'VON_KARMAN_CONST', CS%vonKar, & + 'The value the von Karman constant as used for mixed layer viscosity.', & + units='nondim', default=0.41) call get_param(param_file, mdl, "HMIX_MIN", CS%Hmix_min, & "The minimum mixed layer depth if the mixed layer depth "//& "is determined dynamically.", units="m", default=0.0, scale=GV%m_to_H, & diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index b4cb46830d..a3450bd6e4 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -153,8 +153,7 @@ subroutine make_frazil(h, tv, G, GV, US, CS, p_surf, halo) pressure(i,1) = ps(i) + (0.5*H_to_RL2_T2)*h(i,j,1) enddo do k=2,nz ; do i=is,ie - pressure(i,k) = pressure(i,k-1) + & - (0.5*H_to_RL2_T2) * (h(i,j,k) + h(i,j,k-1)) + pressure(i,k) = pressure(i,k-1) + (0.5*H_to_RL2_T2) * (h(i,j,k) + h(i,j,k-1)) enddo ; enddo endif diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 278fb1ddda..71e47fe792 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -2925,6 +2925,7 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di character(len=48) :: thickness_units character(len=40) :: var_name character(len=160) :: var_descript + logical :: physical_OBL_scheme integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz, nbands, m isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nz = GV%ke IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB @@ -3434,9 +3435,11 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di call internal_tides_init(Time, G, GV, US, param_file, diag, CS%int_tide) endif + physical_OBL_scheme = (CS%use_bulkmixedlayer .or. CS%use_KPP .or. CS%use_energetic_PBL) ! initialize module for setting diffusivities call set_diffusivity_init(Time, G, GV, US, param_file, diag, CS%set_diff_CSp, CS%int_tide, & - halo_TS=CS%halo_TS_diff, double_diffuse=CS%double_diffuse) + halo_TS=CS%halo_TS_diff, double_diffuse=CS%double_diffuse, & + physical_OBL_scheme=physical_OBL_scheme) if (CS%useKPP .and. (CS%double_diffuse .and. .not.CS%use_CVMix_ddiff)) & call MOM_error(FATAL, 'diabatic_driver_init: DOUBLE_DIFFUSION (old method) does not work '//& diff --git a/src/parameterizations/vertical/MOM_diapyc_energy_req.F90 b/src/parameterizations/vertical/MOM_diapyc_energy_req.F90 index 975a11d909..2ddf8b8c7a 100644 --- a/src/parameterizations/vertical/MOM_diapyc_energy_req.F90 +++ b/src/parameterizations/vertical/MOM_diapyc_energy_req.F90 @@ -27,8 +27,9 @@ module MOM_diapyc_energy_req type, public :: diapyc_energy_req_CS ; private logical :: initialized = .false. !< A variable that is here because empty !! structures are not permitted by some compilers. - real :: test_Kh_scaling !< A scaling factor for the diapycnal diffusivity. - real :: ColHt_scaling !< A scaling factor for the column height change correction term. + real :: test_Kh_scaling !< A scaling factor for the diapycnal diffusivity [nondim] + real :: ColHt_scaling !< A scaling factor for the column height change correction term [nondim] + real :: VonKar !< The von Karman coefficient as used in this module [nondim] logical :: use_test_Kh_profile !< If true, use the internal test diffusivity profile in place of !! any that might be passed in as an argument. type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to @@ -104,7 +105,7 @@ subroutine diapyc_energy_req_test(h_3d, dt, tv, G, GV, US, CS, Kd_int) do K=2,nz tmp1 = h_top(K) * h_bot(K) * GV%H_to_Z Kd(K) = CS%test_Kh_scaling * & - ustar * 0.41 * (tmp1*ustar) / (absf*tmp1 + htot*ustar) + ustar * CS%VonKar * (tmp1*ustar) / (absf*tmp1 + htot*ustar) enddo endif may_print = is_root_PE() .and. (i==ie) .and. (j==je) @@ -1292,10 +1293,12 @@ subroutine diapyc_energy_req_init(Time, G, GV, US, param_file, diag, CS) call get_param(param_file, mdl, "ENERGY_REQ_COL_HT_SCALING", CS%ColHt_scaling, & "A scaling factor for the column height change correction "//& "used in testing the energy requirements.", default=1.0, units="nondim") - call get_param(param_file, mdl, "ENERGY_REQ_USE_TEST_PROFILE", & - CS%use_test_Kh_profile, & + call get_param(param_file, mdl, "ENERGY_REQ_USE_TEST_PROFILE", CS%use_test_Kh_profile, & "If true, use the internal test diffusivity profile in "//& "place of any that might be passed in as an argument.", default=.false.) + call get_param(param_file, mdl, 'VON_KARMAN_CONST', CS%vonKar, & + 'The value the von Karman constant as used for mixed layer viscosity.', & + units='nondim', default=0.41) CS%id_ERt = register_diag_field('ocean_model', 'EnReqTest_ERt', diag%axesZi, Time, & "Diffusivity Energy Requirements, top-down", & diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index 0e090b12e3..862f775225 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -36,8 +36,7 @@ module MOM_energetic_PBL logical :: initialized = .false. !< True if this control structure has been initialized. !/ Constants - real :: VonKar = 0.41 !< The von Karman coefficient. This should be a runtime parameter, - !! but because it is set to 0.4 at runtime in KPP it might change answers. + real :: VonKar !< The von Karman coefficient as used in the ePBL module [nondim] real :: omega !< The Earth's rotation rate [T-1 ~> s-1]. real :: omega_frac !< When setting the decay scale for turbulence, use this fraction of !! the absolute rotation rate blended with the local value of f, as @@ -1982,6 +1981,9 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) "A nondimensional scaling factor controlling the inhibition "//& "of the diffusive length scale by rotation. Making this larger "//& "decreases the PBL diffusivity.", units="nondim", default=1.0) + call get_param(param_file, mdl, 'VON_KARMAN_CONST', CS%vonKar, & + 'The value the von Karman constant as used for mixed layer viscosity.', & + units='nondim', default=0.41) call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, & "This sets the default value for the various _ANSWER_DATE parameters.", & default=99991231) diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index 2e27877350..6d35616b3a 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -73,6 +73,8 @@ module MOM_set_diffusivity !! added. logical :: use_LOTW_BBL_diffusivity !< If true, use simpler/less precise, BBL diffusivity. logical :: LOTW_BBL_use_omega !< If true, use simpler/less precise, BBL diffusivity. + real :: Von_Karm !< The von Karman constant as used in the BBL diffusivity calculation + !! [nondim]. See (http://en.wikipedia.org/wiki/Von_Karman_constant) real :: BBL_effic !< efficiency with which the energy extracted !! by bottom drag drives BBL diffusion [nondim] real :: cdrag !< quadratic drag coefficient [nondim] @@ -1406,7 +1408,6 @@ subroutine add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, Kd_int logical :: Rayleigh_drag ! Set to true if there are Rayleigh drag velocities defined in visc, on ! the assumption that this extracted energy also drives diapycnal mixing. integer :: i, k, km1 - real, parameter :: von_karm = 0.41 ! Von Karman constant (http://en.wikipedia.org/wiki/Von_Karman_constant) logical :: do_diag_Kd_BBL if (.not.(CS%bottomdraglaw .and. (CS%BBL_effic > 0.0))) return @@ -1483,7 +1484,7 @@ subroutine add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, Kd_int if ( ustar_D + absf * ( z_bot * D_minus_z ) == 0.) then Kd_wall = 0. else - Kd_wall = ((von_karm * ustar2) * (z_bot * D_minus_z)) & + Kd_wall = ((CS%von_karm * ustar2) * (z_bot * D_minus_z)) & / (ustar_D + absf * (z_bot * D_minus_z)) endif @@ -1963,7 +1964,7 @@ subroutine set_density_ratios(h, tv, kb, G, GV, US, CS, j, ds_dsp1, rho_0) end subroutine set_density_ratios subroutine set_diffusivity_init(Time, G, GV, US, param_file, diag, CS, int_tide_CSp, halo_TS, & - double_diffuse) + double_diffuse, physical_OBL_scheme) type(time_type), intent(in) :: Time !< The current model time type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. @@ -1974,10 +1975,15 @@ subroutine set_diffusivity_init(Time, G, GV, US, param_file, diag, CS, int_tide_ type(set_diffusivity_CS), pointer :: CS !< pointer set to point to the module control !! structure. type(int_tide_CS), intent(in), target :: int_tide_CSp !< Internal tide control struct - integer, optional, intent(out) :: halo_TS !< The halo size of tracer points that must be + integer, intent(out) :: halo_TS !< The halo size of tracer points that must be !! valid for the calculations in set_diffusivity. logical, intent(out) :: double_diffuse !< This indicates whether some version !! of double diffusion is being used. + logical, intent(in) :: physical_OBL_scheme !< If true, a physically based + !! parameterization (like KPP or ePBL or a bulk mixed + !! layer) is used outside of set_diffusivity to + !! specify the mixing that occurs in the ocean's + !! surface boundary layer. ! Local variables real :: decay_length @@ -1990,6 +1996,7 @@ subroutine set_diffusivity_init(Time, G, GV, US, param_file, diag, CS, int_tide_ ! This include declares and sets the variable "version". # include "version_variable.h" character(len=40) :: mdl = "MOM_set_diffusivity" ! This module's name. + real :: vonKar ! The von Karman constant as used for mixed layer viscosity [nomdim] real :: omega_frac_dflt ! The default value for the fraction of the absolute rotation rate ! that is used in place of the absolute value of the local Coriolis ! parameter in the denominator of some expressions [nondim] @@ -2147,6 +2154,12 @@ subroutine set_diffusivity_init(Time, G, GV, US, param_file, diag, CS, int_tide_ call get_param(param_file, mdl, "LOTW_BBL_USE_OMEGA", CS%LOTW_BBL_use_omega, & "If true, use the maximum of Omega and N for the TKE to diffusion "//& "calculation. Otherwise, N is N.", default=.true.) + call get_param(param_file, mdl, 'VON_KARMAN_CONST', vonKar, & + 'The value the von Karman constant as used for mixed layer viscosity.', & + units='nondim', default=0.41) + call get_param(param_file, mdl, 'VON_KARMAN_BBL', CS%von_Karm, & + 'The value the von Karman constant as used in calculating the BBL diffusivity.', & + units='nondim', default=vonKar) endif else CS%use_LOTW_BBL_diffusivity = .false. ! This parameterization depends on a u* from viscous BBL @@ -2164,7 +2177,7 @@ subroutine set_diffusivity_init(Time, G, GV, US, param_file, diag, CS, int_tide_ default=.false., do_not_log=.not.TKE_to_Kd_used) ! set params related to the background mixing - call bkgnd_mixing_init(Time, G, GV, US, param_file, CS%diag, CS%bkgnd_mixing_csp) + call bkgnd_mixing_init(Time, G, GV, US, param_file, CS%diag, CS%bkgnd_mixing_csp, physical_OBL_scheme) call get_param(param_file, mdl, "KV", CS%Kv, & "The background kinematic viscosity in the interior. "//& diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index df6ed9063b..80be1ed12f 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -2050,10 +2050,9 @@ subroutine set_visc_init(Time, G, GV, US, param_file, diag, visc, CS, restart_CS do_not_log=.true.) if (adiabatic) then call log_param(param_file, mdl, "ADIABATIC",adiabatic, & - "There are no diapycnal mass fluxes if ADIABATIC is "//& - "true. This assumes that KD = KDML = 0.0 and that "//& - "there is no buoyancy forcing, but makes the model "//& - "faster by eliminating subroutine calls.", default=.false.) + "There are no diapycnal mass fluxes if ADIABATIC is true. "//& + "This assumes that KD = 0.0 and that there is no buoyancy forcing, "//& + "but makes the model faster by eliminating subroutine calls.", default=.false.) endif if (.not.adiabatic) then @@ -2111,11 +2110,11 @@ subroutine set_visc_init(Time, G, GV, US, param_file, diag, visc, CS, restart_CS endif call get_param(param_file, mdl, "HBBL", CS%Hbbl, & - "The thickness of a bottom boundary layer with a "//& - "viscosity of KVBBL if BOTTOMDRAGLAW is not defined, or "//& - "the thickness over which near-bottom velocities are "//& - "averaged for the drag law if BOTTOMDRAGLAW is defined "//& - "but LINEAR_DRAG is not.", units="m", fail_if_missing=.true.) ! Rescaled later + "The thickness of a bottom boundary layer with a viscosity increased by "//& + "KV_EXTRA_BBL if BOTTOMDRAGLAW is not defined, or the thickness over which "//& + "near-bottom velocities are averaged for the drag law if BOTTOMDRAGLAW is "//& + "defined but LINEAR_DRAG is not.", & + units="m", fail_if_missing=.true.) ! Rescaled later if (CS%bottomdraglaw) then call get_param(param_file, mdl, "CDRAG", CS%cdrag, & "CDRAG is the drag coefficient relating the magnitude of "//& diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index 21ae10fef2..f928327254 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -9,7 +9,7 @@ module MOM_vert_friction use MOM_diag_mediator, only : diag_ctrl use MOM_debugging, only : uvchksum, hchksum use MOM_error_handler, only : MOM_error, FATAL, WARNING, NOTE -use MOM_file_parser, only : get_param, log_version, param_file_type +use MOM_file_parser, only : get_param, log_param, log_version, param_file_type use MOM_forcing_type, only : mech_forcing use MOM_get_input, only : directories use MOM_grid, only : ocean_grid_type @@ -43,11 +43,15 @@ module MOM_vert_friction real :: Hmix !< The mixed layer thickness in thickness units [H ~> m or kg m-2]. real :: Hmix_stress !< The mixed layer thickness over which the wind !! stress is applied with direct_stress [H ~> m or kg m-2]. - real :: Kvml !< The mixed layer vertical viscosity [Z2 T-1 ~> m2 s-1]. + real :: Kvml_invZ2 !< The extra vertical viscosity scale in [Z2 T-1 ~> m2 s-1] in a + !! surface mixed layer with a characteristic thickness given by Hmix, + !! and scaling proportional to (Hmix/z)^2, where z is the distance + !! from the surface; this can get very large with thin layers. real :: Kv !< The interior vertical viscosity [Z2 T-1 ~> m2 s-1]. real :: Hbbl !< The static bottom boundary layer thickness [H ~> m or kg m-2]. - real :: Kvbbl !< The vertical viscosity in the bottom boundary - !! layer [Z2 T-1 ~> m2 s-1]. + real :: Kv_extra_bbl !< An extra vertical viscosity in the bottom boundary layer of thickness + !! Hbbl when there is not a bottom drag law in use [Z2 T-1 ~> m2 s-1]. + real :: vonKar !< The von Karman constant as used for mixed layer viscosity [nomdim] real :: maxvel !< Velocity components greater than maxvel are truncated [L T-1 ~> m s-1]. real :: vel_underflow !< Velocity components smaller than vel_underflow @@ -93,13 +97,23 @@ module MOM_vert_friction real :: harm_BL_val !< A scale to determine when water is in the boundary !! layers based solely on harmonic mean thicknesses !! for the purpose of determining the extent to which - !! the thicknesses used in the viscosities are upwinded. - logical :: direct_stress !< If true, the wind stress is distributed over the - !! topmost Hmix_stress of fluid and KVML may be very small. + !! the thicknesses used in the viscosities are upwinded [nondim]. + logical :: direct_stress !< If true, the wind stress is distributed over the topmost Hmix_stress + !! of fluid, and an added mixed layer viscosity or a physically based + !! boundary layer turbulence parameterization is not needed for stability. logical :: dynamic_viscous_ML !< If true, use the results from a dynamic !! calculation, perhaps based on a bulk Richardson !! number criterion, to determine the mixed layer !! thickness for viscosity. + logical :: fixed_LOTW_ML !< If true, use a Law-of-the-wall prescription for the mixed layer + !! viscosity within a boundary layer that is the lesser of Hmix and the + !! total depth of the ocean in a column. + logical :: apply_LOTW_floor !< If true, use a Law-of-the-wall prescription to set a lower bound + !! on the viscous coupling between layers within the surface boundary + !! layer, based the distance of interfaces from the surface. This only + !! acts when there are large changes in the thicknesses of successive + !! layers or when the viscosity is set externally and the wind stress + !! has subsequently increased. integer :: answer_date !< The vintage of the order of arithmetic and expressions in the viscous !! calculations. Values below 20190101 recover the answers from the end !! of 2018, while higher values use expressions that do not use an @@ -150,7 +164,7 @@ module MOM_vert_friction !! is the interfacial coupling thickness per time step, !! encompassing background viscosity as well as contributions from !! enhanced mixed and bottom layer viscosities. -!! $r_k$ is a Rayleight drag term due to channel drag. +!! $r_k$ is a Rayleigh drag term due to channel drag. !! There is an additional stress term on the right-hand side !! if DIRECT_STRESS is true, applied to the surface layer. @@ -298,7 +312,7 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & ! denote the diagonal of the system as b_k, the subdiagonal as a_k ! and the superdiagonal as c_k. The right-hand side terms are d_k. ! - ! ignoring the rayleigh drag contribution, + ! ignoring the Rayleigh drag contribution, ! we have a_k = -dt_Z_to_H * a_u(k) ! b_k = h_u(k) + dt_Z_to_H * (a_u(k) + a_u(k+1)) ! c_k = -dt_Z_to_H * a_u(k+1) @@ -713,11 +727,11 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC) real, allocatable, dimension(:,:,:) :: Kv_v !< Total vertical viscosity at v-points [Z2 T-1 ~> m2 s-1]. real :: zcol(SZI_(G)) ! The height of an interface at h-points [H ~> m or kg m-2]. real :: botfn ! A function which goes from 1 at the bottom to 0 much more - ! than Hbbl into the interior. + ! than Hbbl into the interior [nondim]. real :: topfn ! A function which goes from 1 at the top to 0 much more - ! than Htbl into the interior. + ! than Htbl into the interior [nondim]. real :: z2 ! The distance from the bottom, normalized by Hbbl [nondim] - real :: z2_wt ! A nondimensional (0-1) weight used when calculating z2. + real :: z2_wt ! A nondimensional (0-1) weight used when calculating z2 [nondim]. real :: z_clear ! The clearance of an interface above the surrounding topography [H ~> m or kg m-2]. real :: a_cpl_max ! The maximum drag coefficient across interfaces, set so that it will be ! representable as a 32-bit float in MKS units [Z T-1 ~> m s-1] @@ -726,7 +740,7 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC) real :: I_valBL ! The inverse of a scaling factor determining when water is ! still within the boundary layer, as determined by the sum - ! of the harmonic mean thicknesses. + ! of the harmonic mean thicknesses [nondim]. logical, dimension(SZIB_(G)) :: do_i, do_i_shelf logical :: do_any_shelf integer, dimension(SZIB_(G)) :: & @@ -1137,10 +1151,12 @@ subroutine find_coupling_coef(a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, intent(in) :: h_harm !< Harmonic mean of thicknesses around a velocity !! grid point [H ~> m or kg m-2] real, dimension(SZIB_(G)), intent(in) :: bbl_thick !< Bottom boundary layer thickness [H ~> m or kg m-2] - real, dimension(SZIB_(G)), intent(in) :: kv_bbl !< Bottom boundary layer viscosity [Z2 T-1 ~> m2 s-1]. + real, dimension(SZIB_(G)), intent(in) :: kv_bbl !< Bottom boundary layer viscosity, exclusive of + !! any depth-dependent contributions from + !! visc%Kv_shear [Z2 T-1 ~> m2 s-1]. real, dimension(SZIB_(G),SZK_(GV)+1), & intent(in) :: z_i !< Estimate of interface heights above the bottom, - !! normalized by the bottom boundary layer thickness + !! normalized by the bottom boundary layer thickness [nondim] real, dimension(SZIB_(G)), intent(out) :: h_ml !< Mixed layer depth [H ~> m or kg m-2] integer, intent(in) :: j !< j-index to find coupling coefficient for real, intent(in) :: dt !< Time increment [T ~> s] @@ -1159,7 +1175,6 @@ subroutine find_coupling_coef(a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, u_star, & ! ustar at a velocity point [Z T-1 ~> m s-1]. absf, & ! The average of the neighboring absolute values of f [T-1 ~> s-1]. ! h_ml, & ! The mixed layer depth [H ~> m or kg m-2]. - nk_visc, & ! The (real) interface index of the base of mixed layer. z_t, & ! The distance from the top, sometimes normalized ! by Hmix, [H ~> m or kg m-2] or [nondim]. kv_TBL, & ! The viscosity in a top boundary layer under ice [Z2 T-1 ~> m2 s-1]. @@ -1167,21 +1182,27 @@ subroutine find_coupling_coef(a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, real, dimension(SZIB_(G),SZK_(GV)+1) :: & Kv_tot, & ! The total viscosity at an interface [Z2 T-1 ~> m2 s-1]. Kv_add ! A viscosity to add [Z2 T-1 ~> m2 s-1]. + integer, dimension(SZIB_(G)) :: & + nk_in_ml ! The index of the deepest interface in the mixed layer. real :: h_shear ! The distance over which shears occur [H ~> m or kg m-2]. - real :: r ! A thickness to compare with Hbbl [H ~> m or kg m-2]. + real :: dhc ! The distance between the center of adjacent layers [H ~> m or kg m-2]. real :: visc_ml ! The mixed layer viscosity [Z2 T-1 ~> m2 s-1]. real :: I_Hmix ! The inverse of the mixed layer thickness [H-1 ~> m-1 or m2 kg-1]. real :: a_ml ! The layer coupling coefficient across an interface in ! the mixed layer [Z T-1 ~> m s-1]. + real :: a_floor ! A lower bound on the layer coupling coefficient across an interface in + ! the mixed layer [Z T-1 ~> m s-1]. real :: I_amax ! The inverse of the maximum coupling coefficient [T Z-1 ~> s m-1]. real :: temp1 ! A temporary variable [H Z ~> m2 or kg m-1] - real :: h_neglect ! A thickness that is so small it is usually lost - ! in roundoff and can be neglected [H ~> m or kg m-2]. + real :: ustar2_denom ! A temporary variable in the surface boundary layer turbulence + ! calculations [Z H-1 T-1 ~> s-1 or m3 kg-1 s-1] + real :: h_neglect ! A thickness that is so small it is usually lost + ! in roundoff and can be neglected [H ~> m or kg m-2]. real :: z2 ! A copy of z_i [nondim] real :: botfn ! A function that is 1 at the bottom and small far from it [nondim] real :: topfn ! A function that is 1 at the top and small far from it [nondim] real :: kv_top ! A viscosity associated with the top boundary layer [Z2 T-1 ~> m2 s-1] - logical :: do_shelf, do_OBCs + logical :: do_shelf, do_OBCs, can_exit integer :: i, k, is, ie, max_nk integer :: nz @@ -1207,36 +1228,33 @@ subroutine find_coupling_coef(a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, if (associated(OBC)) then ; do_OBCS = (OBC%number_of_segments > 0) ; endif h_ml(:) = 0.0 -! The following loop calculates the vertical average velocity and -! surface mixed layer contributions to the vertical viscosity. + ! This top boundary condition is appropriate when the wind stress is determined + ! externally and does not change within a timestep due to the surface velocity. do i=is,ie ; Kv_tot(i,1) = 0.0 ; enddo - if ((GV%nkml>0) .or. do_shelf) then ; do k=2,nz ; do i=is,ie - if (do_i(i)) Kv_tot(i,K) = CS%Kv - enddo ; enddo ; else + do K=2,nz+1 ; do i=is,ie + Kv_tot(i,K) = CS%Kv + enddo ; enddo + + if ((CS%Kvml_invZ2 > 0.0) .and. .not.do_shelf) then + ! This is an older (vintage ~1997) way to prevent wind stresses from driving very + ! large flows in nearly massless near-surface layers when there is not a physically- + ! based surface boundary layer parameterization. It does not have a plausible + ! physical basis, and probably should not be used. I_Hmix = 1.0 / (CS%Hmix + h_neglect) do i=is,ie ; z_t(i) = h_neglect*I_Hmix ; enddo do K=2,nz ; do i=is,ie ; if (do_i(i)) then z_t(i) = z_t(i) + h_harm(i,k-1)*I_Hmix - Kv_tot(i,K) = CS%Kv + CS%Kvml / ((z_t(i)*z_t(i)) * & + Kv_tot(i,K) = CS%Kv + CS%Kvml_invZ2 / ((z_t(i)*z_t(i)) * & (1.0 + 0.09*z_t(i)*z_t(i)*z_t(i)*z_t(i)*z_t(i)*z_t(i))) endif ; enddo ; enddo endif - do i=is,ie ; if (do_i(i)) then - if (CS%bottomdraglaw) then - r = hvel(i,nz)*0.5 - if (r < bbl_thick(i)) then - a_cpl(i,nz+1) = kv_bbl(i) / (I_amax*kv_bbl(i) + (r+h_neglect)*GV%H_to_Z) - else - a_cpl(i,nz+1) = kv_bbl(i) / (I_amax*kv_bbl(i) + (bbl_thick(i)+h_neglect)*GV%H_to_Z) - endif - else - a_cpl(i,nz+1) = CS%Kvbbl / ((0.5*hvel(i,nz)+h_neglect)*GV%H_to_Z + I_amax*CS%Kvbbl) - endif - endif ; enddo - if (associated(visc%Kv_shear)) then - ! The factor of 2 that used to be required in the viscosities is no longer needed. + ! Add in viscosities that are determined by physical processes that are handled in + ! other modules, and which do not respond immediately to the changing layer thicknesses. + ! These processes may include shear-driven mixing or contributions from some boundary + ! layer turbulence schemes. Other viscosity contributions that respond to the evolving + ! layer thicknesses or the surface wind stresses are added later. if (work_on_u) then do K=2,nz ; do i=is,ie ; if (do_i(i)) then Kv_add(i,K) = 0.5*(visc%Kv_shear(i,j,k) + visc%Kv_shear(i+1,j,k)) @@ -1273,6 +1291,9 @@ subroutine find_coupling_coef(a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, endif if (associated(visc%Kv_shear_Bu)) then + ! This is similar to what was done above, but for contributions coming from the corner + ! (vorticity) points. Because OBCs run through the faces and corners there is no need + ! to further modify these viscosities here to take OBCs into account. if (work_on_u) then do K=2,nz ; do I=Is,Ie ; If (do_i(I)) then Kv_tot(I,K) = Kv_tot(I,K) + (0.5)*(visc%Kv_shear_Bu(I,J-1,k) + visc%Kv_shear_Bu(I,J,k)) @@ -1284,29 +1305,71 @@ subroutine find_coupling_coef(a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, endif endif - do K=nz,2,-1 ; do i=is,ie ; if (do_i(i)) then - ! botfn determines when a point is within the influence of the bottom - ! boundary layer, going from 1 at the bottom to 0 in the interior. - z2 = z_i(i,k) - botfn = 1.0 / (1.0 + 0.09*z2*z2*z2*z2*z2*z2) + ! Set the viscous coupling coefficients, excluding surface mixed layer contributions + ! for now, but including viscous bottom drag, working up from the bottom. + if (CS%bottomdraglaw) then + do i=is,ie ; if (do_i(i)) then + dhc = hvel(i,nz)*0.5 + ! These expressions assume that Kv_tot(i,nz+1) = CS%Kv, consistent with + ! the suppression of turbulent mixing by the presence of a solid boundary. + if (dhc < bbl_thick(i)) then + a_cpl(i,nz+1) = kv_bbl(i) / (I_amax*kv_bbl(i) + (dhc+h_neglect)*GV%H_to_Z) + else + a_cpl(i,nz+1) = kv_bbl(i) / (I_amax*kv_bbl(i) + (bbl_thick(i)+h_neglect)*GV%H_to_Z) + endif + endif ; enddo + do K=nz,2,-1 ; do i=is,ie ; if (do_i(i)) then + ! botfn determines when a point is within the influence of the bottom + ! boundary layer, going from 1 at the bottom to 0 in the interior. + z2 = z_i(i,k) + botfn = 1.0 / (1.0 + 0.09*z2*z2*z2*z2*z2*z2) - if (CS%bottomdraglaw) then Kv_tot(i,K) = Kv_tot(i,K) + (kv_bbl(i) - CS%Kv)*botfn - r = 0.5*(hvel(i,k) + hvel(i,k-1)) - if (r > bbl_thick(i)) then - h_shear = ((1.0 - botfn) * r + botfn*bbl_thick(i)) + h_neglect + dhc = 0.5*(hvel(i,k) + hvel(i,k-1)) + if (dhc > bbl_thick(i)) then + h_shear = ((1.0 - botfn) * dhc + botfn*bbl_thick(i)) + h_neglect else - h_shear = r + h_neglect + h_shear = dhc + h_neglect endif - else - Kv_tot(i,K) = Kv_tot(i,K) + (CS%Kvbbl-CS%Kv)*botfn + + ! Calculate the coupling coefficients from the viscosities. + a_cpl(i,K) = Kv_tot(i,K) / (h_shear*GV%H_to_Z + I_amax*Kv_tot(i,K)) + endif ; enddo ; enddo ! i & k loops + elseif (abs(CS%Kv_extra_bbl) > 0.0) then + ! There is a simple enhancement of the near-bottom viscosities, but no adjustment + ! of the viscous coupling length scales to give a particular bottom stress. + do i=is,ie ; if (do_i(i)) then + a_cpl(i,nz+1) = (Kv_tot(i,nz+1) + CS%Kv_extra_bbl) / & + ((0.5*hvel(i,nz)+h_neglect)*GV%H_to_Z + I_amax*(Kv_tot(i,nz+1)+CS%Kv_extra_bbl)) + endif ; enddo + do K=nz,2,-1 ; do i=is,ie ; if (do_i(i)) then + ! botfn determines when a point is within the influence of the bottom + ! boundary layer, going from 1 at the bottom to 0 in the interior. + z2 = z_i(i,k) + botfn = 1.0 / (1.0 + 0.09*z2*z2*z2*z2*z2*z2) + + Kv_tot(i,K) = Kv_tot(i,K) + CS%Kv_extra_bbl*botfn h_shear = 0.5*(hvel(i,k) + hvel(i,k-1) + h_neglect) - endif - ! Calculate the coupling coefficients from the viscosities. - a_cpl(i,K) = Kv_tot(i,K) / (h_shear*GV%H_to_Z + I_amax*Kv_tot(i,K)) - endif ; enddo ; enddo ! i & k loops + ! Calculate the coupling coefficients from the viscosities. + a_cpl(i,K) = Kv_tot(i,K) / (h_shear*GV%H_to_Z + I_amax*Kv_tot(i,K)) + endif ; enddo ; enddo ! i & k loops + else + ! Any near-bottom viscous enhancements were already incorporated into Kv_tot, and there is + ! no adjustment of the viscous coupling length scales to give a particular bottom stress. + do i=is,ie ; if (do_i(i)) then + a_cpl(i,nz+1) = Kv_tot(i,nz+1) / ((0.5*hvel(i,nz)+h_neglect)*GV%H_to_Z + I_amax*Kv_tot(i,nz+1)) + endif ; enddo + do K=nz,2,-1 ; do i=is,ie ; if (do_i(i)) then + h_shear = 0.5*(hvel(i,k) + hvel(i,k-1) + h_neglect) + ! Calculate the coupling coefficients from the viscosities. + a_cpl(i,K) = Kv_tot(i,K) / (h_shear*GV%H_to_Z + I_amax*Kv_tot(i,K)) + endif ; enddo ; enddo ! i & k loops + endif + ! Add surface intensified viscous coupling, either as a no-slip boundary condition under a + ! rigid ice-shelf, or due to wind-stress driven surface boundary layer mixing that has not + ! already been added via visc%Kv_shear. if (do_shelf) then ! Set the coefficients to include the no-slip surface stress. do i=is,ie ; if (do_i(i)) then @@ -1331,68 +1394,165 @@ subroutine find_coupling_coef(a_cpl, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, z_t(i) = z_t(i) + hvel(i,k-1) / tbl_thick(i) topfn = 1.0 / (1.0 + 0.09 * z_t(i)**6) - r = 0.5*(hvel(i,k)+hvel(i,k-1)) - if (r > tbl_thick(i)) then - h_shear = ((1.0 - topfn) * r + topfn*tbl_thick(i)) + h_neglect + dhc = 0.5*(hvel(i,k)+hvel(i,k-1)) + if (dhc > tbl_thick(i)) then + h_shear = ((1.0 - topfn) * dhc + topfn*tbl_thick(i)) + h_neglect else - h_shear = r + h_neglect + h_shear = dhc + h_neglect endif kv_top = topfn * kv_TBL(i) a_cpl(i,K) = a_cpl(i,K) + kv_top / (h_shear*GV%H_to_Z + I_amax*kv_top) endif ; enddo ; enddo - elseif (CS%dynamic_viscous_ML .or. (GV%nkml>0)) then - max_nk = 0 - do i=is,ie ; if (do_i(i)) then - if (GV%nkml>0) nk_visc(i) = real(GV%nkml+1) - if (work_on_u) then + + elseif (CS%dynamic_viscous_ML .or. (GV%nkml>0) .or. CS%fixed_LOTW_ML .or. CS%apply_LOTW_floor) then + + ! Find the friction velocity and the absolute value of the Coriolis parameter at this point. + u_star(:) = 0.0 ! Zero out the friction velocity on land points. + if (work_on_u) then + do I=is,ie ; if (do_i(I)) then u_star(I) = 0.5*(forces%ustar(i,j) + forces%ustar(i+1,j)) absf(I) = 0.5*(abs(G%CoriolisBu(I,J-1)) + abs(G%CoriolisBu(I,J))) - if (CS%dynamic_viscous_ML) nk_visc(I) = visc%nkml_visc_u(I,j) + 1 - else - u_star(i) = 0.5*(forces%ustar(i,j) + forces%ustar(i,j+1)) - absf(i) = 0.5*(abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J))) - if (CS%dynamic_viscous_ML) nk_visc(i) = visc%nkml_visc_v(i,J) + 1 - endif - h_ml(i) = h_neglect ; z_t(i) = 0.0 - max_nk = max(max_nk,ceiling(nk_visc(i) - 1.0)) - endif ; enddo - - if (do_OBCS) then ; if (work_on_u) then - do I=is,ie ; if (do_i(I) .and. (OBC%segnum_u(I,j) /= OBC_NONE)) then + endif ; enddo + if (do_OBCs) then ; do I=is,ie ; if (do_i(I) .and. (OBC%segnum_u(I,j) /= OBC_NONE)) then if (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_E) & u_star(I) = forces%ustar(i,j) if (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_W) & u_star(I) = forces%ustar(i+1,j) - endif ; enddo + endif ; enddo ; endif else - do i=is,ie ; if (do_i(i) .and. (OBC%segnum_v(i,J) /= OBC_NONE)) then + do i=is,ie ; if (do_i(i)) then + u_star(i) = 0.5*(forces%ustar(i,j) + forces%ustar(i,j+1)) + absf(i) = 0.5*(abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J))) + endif ; enddo + if (do_OBCs) then ; do i=is,ie ; if (do_i(i) .and. (OBC%segnum_v(i,J) /= OBC_NONE)) then if (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_N) & u_star(i) = forces%ustar(i,j) if (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_S) & u_star(i) = forces%ustar(i,j+1) + endif ; enddo ; endif + endif + + ! Determine the thickness of the surface ocean boundary layer and its extent in index space. + nk_in_ml(:) = 0 + if (CS%dynamic_viscous_ML) then + ! The fractional number of layers that are within the viscous boundary layer were + ! previously stored in visc%nkml_visc_[uv]. + h_ml(:) = h_neglect + max_nk = 0 + if (work_on_u) then + do i=is,ie ; if (do_i(i)) then + nk_in_ml(I) = ceiling(visc%nkml_visc_u(I,j)) + max_nk = max(max_nk, nk_in_ml(I)) + endif ; enddo + do k=1,max_nk ; do i=is,ie ; if (do_i(i)) then + if (k <= visc%nkml_visc_u(I,j)) then ! This layer is all in the ML. + h_ml(i) = h_ml(i) + hvel(i,k) + elseif (k < visc%nkml_visc_u(I,j) + 1.0) then ! Part of this layer is in the ML. + h_ml(i) = h_ml(i) + ((visc%nkml_visc_u(I,j) + 1.0) - k) * hvel(i,k) + endif + endif ; enddo ; enddo + else + do i=is,ie ; if (do_i(i)) then + nk_in_ml(i) = ceiling(visc%nkml_visc_v(i,J)) + max_nk = max(max_nk, nk_in_ml(i)) + endif ; enddo + do k=1,max_nk ; do i=is,ie ; if (do_i(i)) then + if (k <= visc%nkml_visc_v(i,J)) then ! This layer is all in the ML. + h_ml(i) = h_ml(i) + hvel(i,k) + elseif (k < visc%nkml_visc_v(i,J) + 1.0) then ! Part of this layer is in the ML. + h_ml(i) = h_ml(i) + ((visc%nkml_visc_v(i,J) + 1.0) - k) * hvel(i,k) + endif + endif ; enddo ; enddo + endif + + elseif (GV%nkml>0) then + ! This is a simple application of a refined-bulk mixed layer with GV%nkml sublayers. + max_nk = GV%nkml + do i=is,ie ; if (do_i(i)) then + nk_in_ml(i) = GV%nkml endif ; enddo - endif ; endif - do k=1,max_nk ; do i=is,ie ; if (do_i(i)) then - if (k+1 <= nk_visc(i)) then ! This layer is all in the ML. + h_ml(:) = h_neglect + do k=1,GV%nkml ; do i=is,ie ; if (do_i(i)) then h_ml(i) = h_ml(i) + hvel(i,k) - elseif (k < nk_visc(i)) then ! Part of this layer is in the ML. - h_ml(i) = h_ml(i) + (nk_visc(i) - k) * hvel(i,k) - endif - endif ; enddo ; enddo + endif ; enddo ; enddo + elseif (CS%fixed_LOTW_ML .or. CS%apply_LOTW_floor) then + ! Determine which interfaces are within CS%Hmix of the surface, and set the viscous + ! boundary layer thickness to the the smaller of CS%Hmix and the depth of the ocean. + h_ml(:) = 0.0 + do k=1,nz + can_exit = .true. + do i=is,ie ; if (do_i(i) .and. (h_ml(i) < CS%Hmix)) then + nk_in_ml(i) = k + if (h_ml(i) + hvel(i,k) < CS%Hmix) then + h_ml(i) = h_ml(i) + hvel(i,k) + can_exit = .false. ! Part of the next deeper layer is also in the mixed layer. + else + h_ml(i) = CS%Hmix + endif + endif ; enddo + if (can_exit) exit ! All remaining layers in this row are below the mixed layer depth. + enddo + max_nk = 0 + do i=is,ie ; max_nk = max(max_nk, nk_in_ml(i)) ; enddo + endif + + ! Avoid working on land or on columns where the viscous coupling could not be increased. + do i=is,ie ; if ((u_star(i)<=0.0) .or. (.not.do_i(i))) nk_in_ml(i) = 0 ; enddo + + ! Set the viscous coupling at the interfaces as the larger of what was previously + ! set and the contributions from the surface boundary layer. + z_t(:) = 0.0 + if (CS%apply_LOTW_floor .and. & + (CS%dynamic_viscous_ML .or. (GV%nkml>0) .or. CS%fixed_LOTW_ML)) then + do K=2,max_nk ; do i=is,ie ; if (k <= nk_in_ml(i)) then + z_t(i) = z_t(i) + hvel(i,k-1) + + ! The viscosity in visc_ml is set to go to 0 at the mixed layer top and bottom + ! (in a log-layer) and be further limited by rotation to give the natural Ekman length. + temp1 = (z_t(i)*h_ml(i) - z_t(i)*z_t(i))*GV%H_to_Z + ustar2_denom = (CS%vonKar * u_star(i)**2) / (absf(i)*temp1 + (h_ml(i)+h_neglect)*u_star(i)) + visc_ml = temp1 * ustar2_denom + ! Set the viscous coupling based on the model's vertical resolution. The omission of + ! the I_amax factor here is consistent with answer dates above 20190101. + a_ml = visc_ml / (0.25*(hvel(i,k)+hvel(i,k-1) + h_neglect) * GV%H_to_Z) + + ! As a floor on the viscous coupling, assume that the length scale in the denominator can + ! not be larger than the distance from the surface, consistent with a logarithmic velocity + ! profile. This is consistent with visc_ml, but cancels out common factors of z_t. + a_floor = (h_ml(i) - z_t(i)) * ustar2_denom + + ! Choose the largest estimate of a_cpl. + a_cpl(i,K) = max(a_cpl(i,K), a_ml, a_floor) + ! An option could be added to change this to: a_cpl(i,K) = max(a_cpl(i,K) + a_ml, a_floor) + endif ; enddo ; enddo + elseif (CS%apply_LOTW_floor) then + do K=2,max_nk ; do i=is,ie ; if (k <= nk_in_ml(i)) then + z_t(i) = z_t(i) + hvel(i,k-1) + + temp1 = (z_t(i)*h_ml(i) - z_t(i)*z_t(i))*GV%H_to_Z + ustar2_denom = (CS%vonKar * u_star(i)**2) / (absf(i)*temp1 + (h_ml(i)+h_neglect)*u_star(i)) - do K=2,max_nk ; do i=is,ie ; if (do_i(i)) then ; if (k < nk_visc(i)) then - ! Set the viscosity at the interfaces. - z_t(i) = z_t(i) + hvel(i,k-1) - temp1 = (z_t(i)*h_ml(i) - z_t(i)*z_t(i))*GV%H_to_Z - ! This viscosity is set to go to 0 at the mixed layer top and bottom (in a log-layer) - ! and be further limited by rotation to give the natural Ekman length. - visc_ml = u_star(i) * 0.41 * (temp1*u_star(i)) / (absf(i)*temp1 + (h_ml(i)+h_neglect)*u_star(i)) - a_ml = visc_ml / (0.25*(hvel(i,k)+hvel(i,k-1) + h_neglect) * GV%H_to_Z + 0.5*I_amax*visc_ml) - ! Choose the largest estimate of a. - if (a_ml > a_cpl(i,K)) a_cpl(i,K) = a_ml - endif ; endif ; enddo ; enddo + ! As a floor on the viscous coupling, assume that the length scale in the denominator can not + ! be larger than the distance from the surface, consistent with a logarithmic velocity profile. + a_cpl(i,K) = max(a_cpl(i,K), (h_ml(i) - z_t(i)) * ustar2_denom) + endif ; enddo ; enddo + else + do K=2,max_nk ; do i=is,ie ; if (k <= nk_in_ml(i)) then + z_t(i) = z_t(i) + hvel(i,k-1) + + temp1 = (z_t(i)*h_ml(i) - z_t(i)*z_t(i))*GV%H_to_Z + ! This viscosity is set to go to 0 at the mixed layer top and bottom (in a log-layer) + ! and be further limited by rotation to give the natural Ekman length. + visc_ml = u_star(i) * CS%vonKar * (temp1*u_star(i)) / (absf(i)*temp1 + (h_ml(i)+h_neglect)*u_star(i)) + a_ml = visc_ml / (0.25*(hvel(i,k)+hvel(i,k-1) + h_neglect) * GV%H_to_Z + 0.5*I_amax*visc_ml) + + ! Choose the largest estimate of a_cpl, but these could be changed to be additive. + a_cpl(i,K) = max(a_cpl(i,K), a_ml) + ! An option could be added to change this to: a_cpl(i,K) = a_cpl(i,K) + a_ml + endif ; enddo ; enddo + endif endif end subroutine find_coupling_coef @@ -1626,6 +1786,7 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & ! Local variables real :: Kv_dflt ! A default viscosity [m2 s-1]. + real :: Kv_BBL ! A viscosity in the bottom boundary layer with a simple scheme [Z2 T-1 ~> m2 s-1]. real :: Hmix_m ! A boundary layer thickness [m]. integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags. logical :: default_2018_answers ! The default setting for the various 2018_ANSWERS flags. @@ -1687,23 +1848,38 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & "actual velocity in the bottommost HBBL, depending on "//& "LINEAR_DRAG.", default=.true.) call get_param(param_file, mdl, "DIRECT_STRESS", CS%direct_stress, & - "If true, the wind stress is distributed over the "//& - "topmost HMIX_STRESS of fluid (like in HYCOM), and KVML "//& - "may be set to a very small value.", default=.false.) + "If true, the wind stress is distributed over the topmost HMIX_STRESS of fluid "//& + "(like in HYCOM), and an added mixed layer viscosity or a physically based "//& + "boundary layer turbulence parameterization is not needed for stability.", & + default=.false.) call get_param(param_file, mdl, "DYNAMIC_VISCOUS_ML", CS%dynamic_viscous_ML, & "If true, use a bulk Richardson number criterion to "//& "determine the mixed layer thickness for viscosity.", & default=.false.) + call get_param(param_file, mdl, "FIXED_DEPTH_LOTW_ML", CS%fixed_LOTW_ML, & + "If true, use a Law-of-the-wall prescription for the mixed layer viscosity "//& + "within a boundary layer that is the lesser of HMIX_FIXED and the total "//& + "depth of the ocean in a column.", default=.false.) + call get_param(param_file, mdl, "LOTW_VISCOUS_ML_FLOOR", CS%apply_LOTW_floor, & + "If true, use a Law-of-the-wall prescription to set a lower bound on the "//& + "viscous coupling between layers within the surface boundary layer, based "//& + "the distance of interfaces from the surface. This only acts when there "//& + "are large changes in the thicknesses of successive layers or when the "//& + "viscosity is set externally and the wind stress has subsequently increased.", & + default=.false.) + call get_param(param_file, mdl, 'VON_KARMAN_CONST', CS%vonKar, & + 'The value the von Karman constant as used for mixed layer viscosity.', & + units='nondim', default=0.41) call get_param(param_file, mdl, "U_TRUNC_FILE", CS%u_trunc_file, & "The absolute path to a file into which the accelerations "//& "leading to zonal velocity truncations are written. "//& - "Undefine this for efficiency if this diagnostic is not "//& - "needed.", default=" ", debuggingParam=.true.) + "Undefine this for efficiency if this diagnostic is not needed.", & + default=" ", debuggingParam=.true.) call get_param(param_file, mdl, "V_TRUNC_FILE", CS%v_trunc_file, & "The absolute path to a file into which the accelerations "//& "leading to meridional velocity truncations are written. "//& - "Undefine this for efficiency if this diagnostic is not "//& - "needed.", default=" ", debuggingParam=.true.) + "Undefine this for efficiency if this diagnostic is not needed.", & + default=" ", debuggingParam=.true.) call get_param(param_file, mdl, "HARMONIC_VISC", CS%harmonic_visc, & "If true, use the harmonic mean thicknesses for "//& "calculating the vertical viscosity.", default=.false.) @@ -1724,12 +1900,12 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & if (CS%direct_stress) then if (GV%nkml < 1) then call get_param(param_file, mdl, "HMIX_STRESS", CS%Hmix_stress, & - "The depth over which the wind stress is applied if "//& - "DIRECT_STRESS is true.", units="m", default=Hmix_m, scale=GV%m_to_H) + "The depth over which the wind stress is applied if DIRECT_STRESS is true.", & + units="m", default=Hmix_m, scale=GV%m_to_H) else call get_param(param_file, mdl, "HMIX_STRESS", CS%Hmix_stress, & - "The depth over which the wind stress is applied if "//& - "DIRECT_STRESS is true.", units="m", fail_if_missing=.true., scale=GV%m_to_H) + "The depth over which the wind stress is applied if DIRECT_STRESS is true.", & + units="m", fail_if_missing=.true., scale=GV%m_to_H) endif if (CS%Hmix_stress <= 0.0) call MOM_error(FATAL, "vertvisc_init: " // & "HMIX_STRESS must be set to a positive value if DIRECT_STRESS is true.") @@ -1739,28 +1915,66 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & "The molecular value, ~1e-6 m2 s-1, may be used.", & units="m2 s-1", fail_if_missing=.true., scale=US%m2_s_to_Z2_T, unscaled=Kv_dflt) - if (GV%nkml < 1) call get_param(param_file, mdl, "KVML", CS%Kvml, & - "The kinematic viscosity in the mixed layer. A typical "//& - "value is ~1e-2 m2 s-1. KVML is not used if "//& - "BULKMIXEDLAYER is true. The default is set by KV.", & - units="m2 s-1", default=Kv_dflt, scale=US%m2_s_to_Z2_T) - if (.not.CS%bottomdraglaw) call get_param(param_file, mdl, "KVBBL", CS%Kvbbl, & - "The kinematic viscosity in the benthic boundary layer. "//& - "A typical value is ~1e-2 m2 s-1. KVBBL is not used if "//& - "BOTTOMDRAGLAW is true. The default is set by KV.", & - units="m2 s-1", default=Kv_dflt, scale=US%m2_s_to_Z2_T) + CS%Kvml_invZ2 = 0.0 + if (GV%nkml < 1) then + call get_param(param_file, mdl, "KV_ML_INVZ2", CS%Kvml_invZ2, & + "An extra kinematic viscosity in a mixed layer of thickness HMIX_FIXED, "//& + "with the actual viscosity scaling as 1/(z*HMIX_FIXED)^2, where z is the "//& + "distance from the surface, to allow for finite wind stresses to be "//& + "transmitted through infinitesimally thin surface layers. This is an "//& + "older option for numerical convenience without a strong physical basis, "//& + "and its use is now discouraged.", & + units="m2 s-1", default=-1.0, scale=US%m2_s_to_Z2_T, do_not_log=.true.) + if (CS%Kvml_invZ2 < 0.0) then + call get_param(param_file, mdl, "KVML", CS%Kvml_invZ2, & + "The scale for an extra kinematic viscosity in the mixed layer", & + units="m2 s-1", default=Kv_dflt, scale=US%m2_s_to_Z2_T, do_not_log=.true.) + if (CS%Kvml_invZ2 >= 0.0) & + call MOM_error(WARNING, "KVML is a deprecated parameter. Use KV_ML_INVZ2 instead.") + endif + if (CS%Kvml_invZ2 < 0.0) CS%Kvml_invZ2 = CS%Kv ! Change this default later to 0.0. + call log_param(param_file, mdl, "KV_ML_INVZ2", US%Z2_T_to_m2_s*CS%Kvml_invZ2, & + "An extra kinematic viscosity in a mixed layer of thickness HMIX_FIXED, "//& + "with the actual viscosity scaling as 1/(z*HMIX_FIXED)^2, where z is the "//& + "distance from the surface, to allow for finite wind stresses to be "//& + "transmitted through infinitesimally thin surface layers. This is an "//& + "older option for numerical convenience without a strong physical basis, "//& + "and its use is now discouraged.", & + units="m2 s-1", default=Kv_dflt) + endif + + if (.not.CS%bottomdraglaw) then + call get_param(param_file, mdl, "KV_EXTRA_BBL", CS%Kv_extra_bbl, & + "An extra kinematic viscosity in the benthic boundary layer. "//& + "KV_EXTRA_BBL is not used if BOTTOMDRAGLAW is true.", & + units="m2 s-1", default=0.0, scale=US%m2_s_to_Z2_T, do_not_log=.true.) + if (CS%Kv_extra_bbl == 0.0) then + call get_param(param_file, mdl, "KVBBL", Kv_BBL, & + "An extra kinematic viscosity in the benthic boundary layer. "//& + "KV_EXTRA_BBL is not used if BOTTOMDRAGLAW is true.", & + units="m2 s-1", default=US%Z2_T_to_m2_s*CS%Kv, scale=US%m2_s_to_Z2_T, do_not_log=.true.) + if (abs(Kv_BBL - CS%Kv) > 1.0e-15*abs(CS%Kv)) then + call MOM_error(WARNING, "KVBBL is a deprecated parameter. Use KV_EXTRA_BBL instead.") + CS%Kv_extra_bbl = Kv_BBL - CS%Kv + endif + endif + call log_param(param_file, mdl, "KV_EXTRA_BBL", US%Z2_T_to_m2_s*CS%Kv_extra_bbl, & + "An extra kinematic viscosity in the benthic boundary layer. "//& + "KV_EXTRA_BBL is not used if BOTTOMDRAGLAW is true.", & + units="m2 s-1", default=0.0) + endif call get_param(param_file, mdl, "HBBL", CS%Hbbl, & - "The thickness of a bottom boundary layer with a "//& - "viscosity of KVBBL if BOTTOMDRAGLAW is not defined, or "//& - "the thickness over which near-bottom velocities are "//& - "averaged for the drag law if BOTTOMDRAGLAW is defined "//& - "but LINEAR_DRAG is not.", units="m", fail_if_missing=.true., scale=GV%m_to_H) + "The thickness of a bottom boundary layer with a viscosity increased by "//& + "KV_EXTRA_BBL if BOTTOMDRAGLAW is not defined, or the thickness over which "//& + "near-bottom velocities are averaged for the drag law if BOTTOMDRAGLAW is "//& + "defined but LINEAR_DRAG is not.", & + units="m", fail_if_missing=.true., scale=GV%m_to_H) call get_param(param_file, mdl, "MAXVEL", CS%maxvel, & - "The maximum velocity allowed before the velocity "//& - "components are truncated.", units="m s-1", default=3.0e8, scale=US%m_s_to_L_T) + "The maximum velocity allowed before the velocity components are truncated.", & + units="m s-1", default=3.0e8, scale=US%m_s_to_L_T) call get_param(param_file, mdl, "CFL_BASED_TRUNCATIONS", CS%CFL_based_trunc, & - "If true, base truncations on the CFL number, and not an "//& - "absolute speed.", default=.true.) + "If true, base truncations on the CFL number, and not an absolute speed.", & + default=.true.) call get_param(param_file, mdl, "CFL_TRUNCATE", CS%CFL_trunc, & "The value of the CFL number that will cause velocity "//& "components to be truncated; instability can occur past 0.5.", & @@ -1782,7 +1996,7 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & "Flag to use Stokes drift Mixing via the Lagrangian "//& " current (Eulerian plus Stokes drift). "//& " Still needs work and testing, so not recommended for use.",& - Default=.false.) + default=.false.) !BGR 04/04/2018{ ! StokesMixing is required for MOM6 for some Langmuir mixing parameterization. ! The code used here has not been developed for vanishing layers or in