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 and use find_ustar #406

Merged
merged 2 commits into from
Jul 26, 2023
Merged

+*Add and use find_ustar #406

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
245bf6f
+*Add and use find_ustar
Hallberg-NOAA Jul 18, 2023
029f799
Merge branch 'dev/gfdl' into add_find_ustar
marshallward Jul 26, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
141 changes: 140 additions & 1 deletion src/core/MOM_forcing_type.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@ module MOM_forcing_type

public extractFluxes1d, extractFluxes2d, optics_type
public MOM_forcing_chksum, MOM_mech_forcing_chksum
public calculateBuoyancyFlux1d, calculateBuoyancyFlux2d
public calculateBuoyancyFlux1d, calculateBuoyancyFlux2d, find_ustar
public forcing_accumulate, fluxes_accumulate
public forcing_SinglePointPrint, mech_forcing_diags, forcing_diagnostics
public register_forcing_type_diags, allocate_forcing_type, deallocate_forcing_type
Expand All @@ -53,6 +53,12 @@ module MOM_forcing_type
module procedure allocate_mech_forcing_from_ref
end interface allocate_mech_forcing

!> Determine the friction velocity from a forcing type or a mechanical forcing type.
interface find_ustar
module procedure find_ustar_fluxes
module procedure find_ustar_mech_forcing
end interface find_ustar

! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
Expand Down Expand Up @@ -1077,6 +1083,139 @@ subroutine calculateBuoyancyFlux2d(G, GV, US, fluxes, optics, h, Temp, Salt, tv,
end subroutine calculateBuoyancyFlux2d


!> Determine the friction velocity from the contenxts of a forcing type, perhaps
!! using the evolving surface density.
subroutine find_ustar_fluxes(fluxes, tv, U_star, G, GV, US, halo, H_T_units)
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
type(forcing), intent(in) :: fluxes !< Surface fluxes container
type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any
!! available thermodynamic fields.
real, dimension(SZI_(G),SZJ_(G)), &
intent(out) :: U_star !< The surface friction velocity [Z T-1 ~> m s-1] or
!! [H T-1 ~> m s-1 or kg m-2 s-1], depending on H_T_units.
integer, optional, intent(in) :: halo !< The extra halo size to fill in, 0 by default
logical, optional, intent(in) :: H_T_units !< If present and true, return U_star in units
!! of [H T-1 ~> m s-1 or kg m-2 s-1]

! Local variables
real :: I_rho ! The inverse of the reference density times a ratio of scaling
! factors [Z L-1 R-1 ~> m3 kg-1] or in some semi-Boussinesq cases
! the rescaled reference density [H2 Z-1 L-1 R-1 ~> m3 kg-1 or kg m-3]
logical :: Z_T_units ! If true, U_star is returned in units of [Z T-1 ~> m s-1], otherwise it is
! returned in [H T-1 ~> m s-1 or kg m-2 s-1]
integer :: i, j, k, is, ie, js, je, hs

hs = 0 ; if (present(halo)) hs = max(halo, 0)
is = G%isc - hs ; ie = G%iec + hs ; js = G%jsc - hs ; je = G%jec + hs

Z_T_units = .true. ; if (present(H_T_units)) Z_T_units = .not.H_T_units

if (.not.(associated(fluxes%ustar) .or. associated(fluxes%tau_mag))) &
call MOM_error(FATAL, "find_ustar_fluxes requires that either ustar or tau_mag be associated.")

if (associated(fluxes%ustar) .and. (GV%Boussinesq .or. .not.associated(fluxes%tau_mag))) then
if (Z_T_units) then
do j=js,je ; do i=is,ie
U_star(i,j) = fluxes%ustar(i,j)
enddo ; enddo
else
do j=js,je ; do i=is,ie
U_star(i,j) = GV%Z_to_H * fluxes%ustar(i,j)
enddo ; enddo
endif
elseif (allocated(tv%SpV_avg)) then
if (tv%valid_SpV_halo < 0) call MOM_error(FATAL, &
"find_ustar_fluxes called in non-Boussinesq mode with invalid values of SpV_avg.")
if (tv%valid_SpV_halo < hs) call MOM_error(FATAL, &
"find_ustar_fluxes called in non-Boussinesq mode with insufficient valid values of SpV_avg.")
if (Z_T_units) then
do j=js,je ; do i=is,ie
U_star(i,j) = sqrt(US%L_to_Z*fluxes%tau_mag(i,j) * tv%SpV_avg(i,j,1))
enddo ; enddo
else
do j=js,je ; do i=is,ie
U_star(i,j) = GV%RZ_to_H * sqrt(US%L_to_Z*fluxes%tau_mag(i,j) / tv%SpV_avg(i,j,1))
enddo ; enddo
endif
else
I_rho = US%L_to_Z * GV%Z_to_H * GV%RZ_to_H
if (Z_T_units) I_rho = US%L_to_Z * GV%H_to_Z * GV%RZ_to_H ! == US%L_to_Z / GV%Rho0
do j=js,je ; do i=is,ie
U_star(i,j) = sqrt(fluxes%tau_mag(i,j) * I_rho)
enddo ; enddo
endif

end subroutine find_ustar_fluxes


!> Determine the friction velocity from the contenxts of a forcing type, perhaps
!! using the evolving surface density.
subroutine find_ustar_mech_forcing(forces, tv, U_star, G, GV, US, halo, H_T_units)
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
type(mech_forcing), intent(in) :: forces !< Surface forces container
type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any
!! available thermodynamic fields.
real, dimension(SZI_(G),SZJ_(G)), &
intent(out) :: U_star !< The surface friction velocity [Z T-1 ~> m s-1]
integer, optional, intent(in) :: halo !< The extra halo size to fill in, 0 by default
logical, optional, intent(in) :: H_T_units !< If present and true, return U_star in units
!! of [H T-1 ~> m s-1 or kg m-2 s-1]

! Local variables
real :: I_rho ! The inverse of the reference density times a ratio of scaling
! factors [Z L-1 R-1 ~> m3 kg-1] or in some semi-Boussinesq cases
! the rescaled reference density [H2 Z-1 L-1 R-1 ~> m3 kg-1 or kg m-3]
logical :: Z_T_units ! If true, U_star is returned in units of [Z T-1 ~> m s-1], otherwise it is
! returned in [H T-1 ~> m s-1 or kg m-2 s-1]
integer :: i, j, k, is, ie, js, je, hs

hs = 0 ; if (present(halo)) hs = max(halo, 0)
is = G%isc - hs ; ie = G%iec + hs ; js = G%jsc - hs ; je = G%jec + hs

Z_T_units = .true. ; if (present(H_T_units)) Z_T_units = .not.H_T_units

if (.not.(associated(forces%ustar) .or. associated(forces%tau_mag))) &
call MOM_error(FATAL, "find_ustar_mech requires that either ustar or tau_mag be associated.")

if (associated(forces%ustar) .and. (GV%Boussinesq .or. .not.associated(forces%tau_mag))) then
if (Z_T_units) then
do j=js,je ; do i=is,ie
U_star(i,j) = forces%ustar(i,j)
enddo ; enddo
else
do j=js,je ; do i=is,ie
U_star(i,j) = GV%Z_to_H * forces%ustar(i,j)
enddo ; enddo
endif
elseif (allocated(tv%SpV_avg)) then
if (tv%valid_SpV_halo < 0) call MOM_error(FATAL, &
"find_ustar_mech called in non-Boussinesq mode with invalid values of SpV_avg.")
if (tv%valid_SpV_halo < hs) call MOM_error(FATAL, &
"find_ustar_mech called in non-Boussinesq mode with insufficient valid values of SpV_avg.")
if (Z_T_units) then
do j=js,je ; do i=is,ie
U_star(i,j) = sqrt(US%L_to_Z*forces%tau_mag(i,j) * tv%SpV_avg(i,j,1))
enddo ; enddo
else
do j=js,je ; do i=is,ie
U_star(i,j) = GV%RZ_to_H * sqrt(US%L_to_Z*forces%tau_mag(i,j) / tv%SpV_avg(i,j,1))
enddo ; enddo
endif
else
I_rho = US%L_to_Z * GV%Z_to_H * GV%RZ_to_H
if (Z_T_units) I_rho = US%L_to_Z * GV%H_to_Z * GV%RZ_to_H ! == US%L_to_Z / GV%Rho0
do j=js,je ; do i=is,ie
U_star(i,j) = sqrt(forces%tau_mag(i,j) * I_rho)
enddo ; enddo
endif

end subroutine find_ustar_mech_forcing


!> Write out chksums for thermodynamic fluxes.
subroutine MOM_forcing_chksum(mesg, fluxes, G, US, haloshift)
character(len=*), intent(in) :: mesg !< message
Expand Down
36 changes: 27 additions & 9 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 @@ -11,7 +11,7 @@ module MOM_mixed_layer_restrat
use MOM_error_handler, only : MOM_error, FATAL, WARNING
use MOM_file_parser, only : get_param, log_version, param_file_type
use MOM_file_parser, only : openParameterBlock, closeParameterBlock
use MOM_forcing_type, only : mech_forcing
use MOM_forcing_type, only : mech_forcing, find_ustar
use MOM_grid, only : ocean_grid_type
use MOM_hor_index, only : hor_index_type
use MOM_lateral_mixing_coeffs, only : VarMix_CS
Expand Down Expand Up @@ -184,6 +184,9 @@ subroutine mixedlayer_restrat_OM4(h, uhtr, vhtr, tv, forces, dt, MLD_in, VarMix,
h_avail ! The volume available for diffusion out of each face of each
! sublayer of the mixed layer, divided by dt [H L2 T-1 ~> m3 s-1 or kg s-1].
real, dimension(SZI_(G),SZJ_(G)) :: &
U_star_2d, & ! The wind friction velocity, calculated using
! the Boussinesq reference density or the time-evolving surface density
! in non-Boussinesq mode [Z T-1 ~> m s-1]
MLD_fast, & ! Mixed layer depth actually used in MLE restratification parameterization [H ~> m or kg m-2]
htot_fast, & ! The sum of the thicknesses of layers in the mixed layer [H ~> m or kg m-2]
Rml_av_fast, & ! g_Rho0 times the average mixed layer density [L2 Z-1 T-2 ~> m s-2]
Expand Down Expand Up @@ -254,6 +257,9 @@ subroutine mixedlayer_restrat_OM4(h, uhtr, vhtr, tv, forces, dt, MLD_in, VarMix,
call MOM_error(FATAL, "mixedlayer_restrat_OM4: "// &
"The resolution argument, Rd/dx, was not associated.")

! Extract the friction velocity from the forcing type.
call find_ustar(forces, tv, U_star_2d, G, GV, US, halo=1)

if (CS%MLE_density_diff > 0.) then ! We need to calculate a mixed layer depth, MLD.
!! TODO: use derivatives and mid-MLD pressure. Currently this is sigma-0. -AJA
pRef_MLD(:) = 0.
Expand Down Expand Up @@ -408,7 +414,7 @@ subroutine mixedlayer_restrat_OM4(h, uhtr, vhtr, tv, forces, dt, MLD_in, VarMix,

if (CS%debug) then
call hchksum(h,'mixed_layer_restrat: h', G%HI, haloshift=1, scale=GV%H_to_m)
call hchksum(forces%ustar,'mixed_layer_restrat: u*', G%HI, haloshift=1, scale=US%Z_to_m*US%s_to_T)
call hchksum(U_star_2d, 'mixed_layer_restrat: u*', G%HI, haloshift=1, scale=US%Z_to_m*US%s_to_T)
call hchksum(MLD_fast,'mixed_layer_restrat: MLD', G%HI, haloshift=1, scale=GV%H_to_m)
call hchksum(Rml_av_fast,'mixed_layer_restrat: rml', G%HI, haloshift=1, &
scale=US%m_to_Z*US%L_T_to_m_s**2)
Expand All @@ -421,7 +427,7 @@ subroutine mixedlayer_restrat_OM4(h, uhtr, vhtr, tv, forces, dt, MLD_in, VarMix,
! U - Component
!$OMP do
do j=js,je ; do I=is-1,ie
u_star = max(CS%ustar_min, 0.5*(forces%ustar(i,j) + forces%ustar(i+1,j)))
u_star = max(CS%ustar_min, 0.5*(U_star_2d(i,j) + U_star_2d(i+1,j)))

absf = 0.5*(abs(G%CoriolisBu(I,J-1)) + abs(G%CoriolisBu(I,J)))
! If needed, res_scaling_fac = min( ds, L_d ) / l_f
Expand Down Expand Up @@ -508,7 +514,7 @@ subroutine mixedlayer_restrat_OM4(h, uhtr, vhtr, tv, forces, dt, MLD_in, VarMix,
! V- component
!$OMP do
do J=js-1,je ; do i=is,ie
u_star = max(CS%ustar_min, 0.5*(forces%ustar(i,j) + forces%ustar(i,j+1)))
u_star = max(CS%ustar_min, 0.5*(U_star_2d(i,j) + U_star_2d(i,j+1)))

absf = 0.5*(abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J)))
! If needed, res_scaling_fac = min( ds, L_d ) / l_f
Expand Down Expand Up @@ -711,6 +717,9 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d
wpup ! Turbulent vertical momentum [ ????? ~> m2 s-2]
real :: uDml_diag(SZIB_(G),SZJ_(G)) ! A 2D copy of uDml for diagnostics [H L2 T-1 ~> m3 s-1 or kg s-1]
real :: vDml_diag(SZI_(G),SZJB_(G)) ! A 2D copy of vDml for diagnostics [H L2 T-1 ~> m3 s-1 or kg s-1]
real :: U_star_2d(SZI_(G),SZJ_(G)) ! The wind friction velocity, calculated using the Boussinesq
! reference density or the time-evolving surface density in non-Boussinesq
! mode [Z T-1 ~> m s-1]
real :: covTS(SZI_(G)) ! SGS TS covariance in Stanley param; currently 0 [degC ppt]
real :: varS(SZI_(G)) ! SGS S variance in Stanley param; currently 0 [ppt2]
real :: dmu(SZK_(GV)) ! Change in mu(z) across layer k [nondim]
Expand Down Expand Up @@ -762,12 +771,15 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d

call pass_var(bflux, G%domain, halo=1)

! Extract the friction velocity from the forcing type.
call find_ustar(forces, tv, U_star_2d, G, GV, US, halo=1)

if (CS%debug) then
call hchksum(h,'mixed_Bodner: h', G%HI, haloshift=1, scale=GV%H_to_m)
call hchksum(BLD, 'mle_Bodner: BLD in', G%HI, haloshift=1, scale=US%Z_to_m)
if (associated(bflux)) &
call hchksum(bflux, 'mle_Bodner: bflux', G%HI, haloshift=1, scale=US%Z_to_m**2*US%s_to_T**3)
call hchksum(forces%ustar,'mle_Bodner: u*', G%HI, haloshift=1, scale=US%Z_to_m*US%s_to_T)
call hchksum(U_star_2d, 'mle_Bodner: u*', G%HI, haloshift=1, scale=US%Z_to_m*US%s_to_T)
call hchksum(CS%MLD_filtered, 'mle_Bodner: MLD_filtered 1', &
G%HI, haloshift=1, scale=US%Z_to_m)
call hchksum(CS%MLD_filtered_slow,'mle_Bodner: MLD_filtered_slow 1', &
Expand All @@ -793,7 +805,7 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d
do j = js-1, je+1 ; do i = is-1, ie+1
w_star3 = max(0., -bflux(i,j)) * BLD(i,j) & ! (this line in Z3 T-3 ~> m3 s-3)
* ( ( US%Z_to_m * US%s_to_T )**3 ) ! m3 s-3
u_star3 = ( US%Z_to_m * US%s_to_T * forces%ustar(i,j) )**3 ! m3 s-3
u_star3 = ( US%Z_to_m * US%s_to_T * U_star_2d(i,j) )**3 ! m3 s-3
wpup(i,j) = max( CS%min_wstar2, & ! The max() avoids division by zero later
( CS%mstar * u_star3 + CS%nstar * w_star3 )**two_thirds ) & ! (this line m2 s-2)
* ( ( US%m_to_Z * US%T_to_s )**2 ) ! Z2 T-2 ~> m2 s-2
Expand Down Expand Up @@ -965,7 +977,7 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d

! Offer diagnostic fields for averaging.
if (query_averaging_enabled(CS%diag)) then
if (CS%id_ustar > 0) call post_data(CS%id_ustar, forces%ustar, CS%diag)
if (CS%id_ustar > 0) call post_data(CS%id_ustar, U_star_2d, CS%diag)
if (CS%id_bflux > 0) call post_data(CS%id_bflux, bflux, CS%diag)
if (CS%id_wpup > 0) call post_data(CS%id_wpup, wpup, CS%diag)
if (CS%id_Rml > 0) call post_data(CS%id_Rml, buoy_av, CS%diag)
Expand Down Expand Up @@ -1053,6 +1065,9 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS)
h_avail ! The volume available for diffusion out of each face of each
! sublayer of the mixed layer, divided by dt [H L2 T-1 ~> m3 s-1 or kg s-1].
real, dimension(SZI_(G),SZJ_(G)) :: &
U_star_2d, & ! The wind friction velocity, calculated using
! the Boussinesq reference density or the time-evolving surface density
! in non-Boussinesq mode [Z T-1 ~> m s-1]
htot, & ! The sum of the thicknesses of layers in the mixed layer [H ~> m or kg m-2]
Rml_av ! g_Rho0 times the average mixed layer density [L2 Z-1 T-2 ~> m s-2]
real :: g_Rho0 ! G_Earth/Rho0 [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1]
Expand Down Expand Up @@ -1110,6 +1125,9 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS)
if (CS%use_Stanley_ML) call MOM_error(FATAL, "mixedlayer_restrat_BML: "// &
"The Stanley parameterization is not available with the BML.")

! Extract the friction velocity from the forcing type.
call find_ustar(forces, tv, U_star_2d, G, GV, US, halo=1)

! Fix this later for nkml >= 3.

p0(:) = 0.0
Expand Down Expand Up @@ -1145,7 +1163,7 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS)
do j=js,je ; do I=is-1,ie
h_vel = 0.5*(htot(i,j) + htot(i+1,j)) * GV%H_to_Z

u_star = max(CS%ustar_min, 0.5*(forces%ustar(i,j) + forces%ustar(i+1,j)))
u_star = max(CS%ustar_min, 0.5*(U_star_2d(i,j) + U_star_2d(i+1,j)))

absf = 0.5*(abs(G%CoriolisBu(I,J-1)) + abs(G%CoriolisBu(I,J)))

Expand Down Expand Up @@ -1196,7 +1214,7 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS)
do J=js-1,je ; do i=is,ie
h_vel = 0.5*(htot(i,j) + htot(i,j+1)) * GV%H_to_Z

u_star = max(CS%ustar_min, 0.5*(forces%ustar(i,j) + forces%ustar(i,j+1)))
u_star = max(CS%ustar_min, 0.5*(U_star_2d(i,j) + U_star_2d(i,j+1)))

absf = 0.5*(abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J)))

Expand Down
Loading