Skip to content
New issue

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

Sign up for GitHub

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

Already on GitHub? Sign in to your account

+(*)Add new viscosity options and discourage the use of KVML #201

Merged
merged 9 commits into from
Sep 10, 2022
Merged

+(*)Add new viscosity options and discourage the use of KVML #201

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
39f8ee0
(*)+Rename KDML to KD_ML_TOT and restrict its use
Hallberg-NOAA Sep 1, 2022
107d367
+(*)Replace KVBBL with KV_EXTRA_BBL
Hallberg-NOAA Sep 1, 2022
fa6d8e3
+Spell HARMONIC right in PORBAR_ETA_INTERP docs
Hallberg-NOAA Sep 1, 2022
4ca424a
+Replace KVML with KV_ML_INVZ2
Hallberg-NOAA Sep 2, 2022
f267462
Refactor boundary layer code in find_coupling_coef
Hallberg-NOAA Sep 2, 2022
5b53ff9
+Add FIXED_DEPTH_LOTW_ML and LOTW_VISCOUS_ML_FLOOR
Hallberg-NOAA Sep 3, 2022
318303a
+(*)Set ustar with wind_forcing_2gyre
Hallberg-NOAA Sep 5, 2022
59b8b06
+Made the von Karman constant a runtime parameter
Hallberg-NOAA Sep 5, 2022
98465af
Merge branch 'dev/gfdl' into KVML_cleanup
marshallward Sep 8, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
72 changes: 41 additions & 31 deletions config_src/drivers/solo_driver/MOM_surface_forcing.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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

Expand All @@ -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

Expand Down Expand Up @@ -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")
Expand All @@ -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]
Expand Down Expand Up @@ -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

Expand All @@ -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. /)
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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)
Expand Down
7 changes: 3 additions & 4 deletions src/core/MOM.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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, "//&
Expand Down
6 changes: 3 additions & 3 deletions src/core/MOM_porous_barriers.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down Expand Up @@ -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
Expand Down
35 changes: 21 additions & 14 deletions src/parameterizations/lateral/MOM_mixed_layer_restrat.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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].
Expand Down Expand Up @@ -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
Expand All @@ -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.) &
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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
Expand Down
Loading