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

Brankart PGF update #1193

Merged
merged 8 commits into from
Sep 9, 2020
Merged

Brankart PGF update #1193

Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
859abbd
Add diagnostic of SGS T variance in Stanley PGF
adcroft Aug 4, 2020
f549dd0
Use thickness weighted variance estimate in PGF
adcroft Aug 4, 2020
995e422
Fix spelling in comments and line length
adcroft Aug 31, 2020
54cf60d
Merge branch 'dev/gfdl' into brankart
adcroft Aug 31, 2020
63a2e2a
Merge branch 'dev/gfdl' into brankart
adcroft Sep 1, 2020
4dbc126
Adds land masking of thicknesses in SGS estimate
adcroft Sep 1, 2020
edff25e
Rename intermediate variable in Stanley PGF parameterization
adcroft Sep 4, 2020
365dd7c
Merge branch 'dev/gfdl' into brankart
marshallward Sep 9, 2020
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
58 changes: 48 additions & 10 deletions src/core/MOM_PressureForce_FV.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -58,6 +58,7 @@ module MOM_PressureForce_FV
!! the mean temperature gradient in the deterministic part of
!! the Stanley form of the Brankart correction.
integer :: id_e_tidal = -1 !< Diagnostic identifier
integer :: id_tvar_sgs = -1 !< Diagnostic identifier
type(tidal_forcing_CS), pointer :: tides_CSp => NULL() !< Tides control structure
end type PressureForce_FV_CS

Expand Down Expand Up @@ -479,7 +480,11 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
logical :: use_ALE ! If true, use an ALE pressure reconstruction.
logical :: use_EOS ! If true, density is calculated from T & S using an equation of state.
type(thermo_var_ptrs) :: tv_tmp! A structure of temporary T & S.
real :: dTdi2, dTdj2 ! Differences in T variance [degC2]
real :: Tl(5) ! copy and T in local stencil [degC]
real :: mn_T ! mean of T in local stencil [degC]
real :: mn_T2 ! mean of T**2 in local stencil [degC]
real :: hl(5) ! Copy of local stencil of H [H ~> m]
real :: r_sm_H ! Reciprocal of sum of H in local stencil [H-1 ~> m-1]
real, parameter :: C1_6 = 1.0/6.0
integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb
Expand All @@ -503,15 +508,43 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
if (CS%Stanley_T2_det_coeff>=0.) then
if (.not. associated(tv%varT)) call safe_alloc_ptr(tv%varT, G%isd, G%ied, G%jsd, G%jed, GV%ke)
do k=1, nz ; do j=js-1,je+1 ; do i=is-1,ie+1
! SGS variance in i-direction [degC2]
dTdi2 = ( ( G%mask2dCu(I ,j) * G%IdxCu(I ,j) * ( tv%T(i+1,j,k) - tv%T(i,j,k) ) &
+ G%mask2dCu(I-1,j) * G%IdxCu(I-1,j) * ( tv%T(i,j,k) - tv%T(i-1,j,k) ) &
) * G%dxT(i,j) * 0.5 )**2
! SGS variance in j-direction [degC2]
dTdj2 = ( ( G%mask2dCv(i,J ) * G%IdyCv(i,J ) * ( tv%T(i,j+1,k) - tv%T(i,j,k) ) &
+ G%mask2dCv(i,J-1) * G%IdyCv(i,J-1) * ( tv%T(i,j,k) - tv%T(i,j-1,k) ) &
) * G%dyT(i,j) * 0.5 )**2
tv%varT(i,j,k) = CS%Stanley_T2_det_coeff * 0.5 * ( dTdi2 + dTdj2 )
! Strictly speaking we should estimate the *horizontal* grid-scale variance
! but neither of the following blocks make a rotation to the horizontal
! and instead work along coordinate.

! This block calculates a simple |delta T| along coordinates and does
! not allow vanishing layer thicknesses or layers tracking topography
!! SGS variance in i-direction [degC2]
!dTdi2 = ( ( G%mask2dCu(I ,j) * G%IdxCu(I ,j) * ( tv%T(i+1,j,k) - tv%T(i,j,k) ) &
! + G%mask2dCu(I-1,j) * G%IdxCu(I-1,j) * ( tv%T(i,j,k) - tv%T(i-1,j,k) ) &
! ) * G%dxT(i,j) * 0.5 )**2
!! SGS variance in j-direction [degC2]
!dTdj2 = ( ( G%mask2dCv(i,J ) * G%IdyCv(i,J ) * ( tv%T(i,j+1,k) - tv%T(i,j,k) ) &
! + G%mask2dCv(i,J-1) * G%IdyCv(i,J-1) * ( tv%T(i,j,k) - tv%T(i,j-1,k) ) &
! ) * G%dyT(i,j) * 0.5 )**2
!tv%varT(i,j,k) = CS%Stanley_T2_det_coeff * 0.5 * ( dTdi2 + dTdj2 )

! This block does a thickness weighted variance calculation and helps control for
! extreme gradients along layers which are vanished against topography. It is
! still a poor approximation in the interior when coordinates are strongly tilted.
hl(1) = h(i,j,k) * G%mask2dT(i,j)
hl(2) = h(i-1,j,k) * G%mask2dCu(I-1,j)
hl(3) = h(i+1,j,k) * G%mask2dCu(I,j)
hl(4) = h(i,j-1,k) * G%mask2dCv(i,J-1)
hl(5) = h(i,j+1,k) * G%mask2dCv(i,J)
r_sm_H = 1. / ( ( hl(1) + ( ( hl(2) + hl(3) ) + ( hl(4) + hl(5) ) ) ) + GV%H_subroundoff )
! Mean of T
Tl(1) = tv%T(i,j,k) ; Tl(2) = tv%T(i-1,j,k) ; Tl(3) = tv%T(i+1,j,k)
Tl(4) = tv%T(i,j-1,k) ; Tl(5) = tv%T(i,j+1,k)
mn_T = ( hl(1)*Tl(1) + ( ( hl(2)*Tl(2) + hl(3)*Tl(3) ) + ( hl(4)*Tl(4) + hl(5)*Tl(5) ) ) ) * r_sm_H
! Adjust T vectors to have zero mean
Tl(:) = Tl(:) - mn_T ; mn_T = 0.
! Variance of T
mn_T2 = ( hl(1)*Tl(1)*Tl(1) + ( ( hl(2)*Tl(2)*Tl(2) + hl(3)*Tl(3)*Tl(3) ) &
+ ( hl(4)*Tl(4)*Tl(4) + hl(5)*Tl(5)*Tl(5) ) ) ) * r_sm_H
! Variance should be positive but round-off can violate this. Calculating
! variance directly would fix this but requires more operations.
tv%varT(i,j,k) = CS%Stanley_T2_det_coeff * max(0., mn_T2 - mn_T*mn_T)
enddo ; enddo ; enddo
endif

Expand Down Expand Up @@ -758,6 +791,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
endif

if (CS%id_e_tidal>0) call post_data(CS%id_e_tidal, e_tidal, CS%diag)
if (CS%id_tvar_sgs>0) call post_data(CS%id_tvar_sgs, tv%varT, CS%diag)

end subroutine PressureForce_FV_Bouss

Expand Down Expand Up @@ -826,6 +860,10 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, tides_CS
"the mean temperature gradient in the deterministic part of "// &
"the Stanley form of the Brankart correction. "// &
"Negative values disable the scheme.", units="nondim", default=-1.0)
if (CS%Stanley_T2_det_coeff>=0.) then
CS%id_tvar_sgs = register_diag_field('ocean_model', 'tvar_sgs_pgf', diag%axesTL, &
Time, 'SGS temperature variance used in PGF', 'degC2')
endif
if (CS%tides) then
CS%id_e_tidal = register_diag_field('ocean_model', 'e_tidal', diag%axesT1, &
Time, 'Tidal Forcing Astronomical and SAL Height Anomaly', 'meter', conversion=US%Z_to_m)
Expand Down