Skip to content

Commit

Permalink
Corrections to cell-averaged density computation (#213)
Browse files Browse the repository at this point in the history
* initial hooks for stochastic EOS modifications

* remove debug statements

* add documentation

* Change ampltiude from 0.39 to sqrt(.39)

* remove global_indexing logic from stoch_eos_init

* switch to using MOM_random and add restart capability

* update random sequence to update each each time-step

* remove tseed0 from MOM_random (leftover from debugging)

* Added necessary submodules and S^2, T^2 diagnostics to MOM_diagnostics

* Added diagnostics for outputting variables related to the stochastic parameterization.

* Diagnostics in MOM_PressureForce_FV updated for stochastic (rather than deterministic) Stanley SGS T variance parameterization.

* Added parentheses for reproducibility.

* Changed diagnostics to account for possible absence of stoch_eos_pattern in MOM_PressureForce_FV,
when deterministic parameterization is on.

* remove mom6_da_hooks and geoKdTree from pkg

* add stochastic compoment to MOM_thickness_diffuse

* fix array size declaration and post_data

* Corrected indexing of loops in MOM_calc_varT

* Changed how parameterization of SGS T variance (deterministic and stochastic) is switched on in PGF and thickness diffusion codes

* Corrected a few typos

* Cleaned up indices, redundant diagnostic, printing

* Fixed diagnostic IDs

* Fixed diagnostics typo

* Corrected indices in calculation of tv%varT

* Minor index fix

* Corrected bug in pressure in Stanley diagnostics

* Fixed whitespace error

* Stoch eos clock (#5)

*Added a clock for the Stanley parameterization

Co-authored-by: jkenigson <jkenigso@gmail.com>

* add halo update to random pattern

* Update MOM_stoch_eos.F90

Fix bug for looping over compute domain (is -> isc etc.)

* Avoid unnessary computations on halo (MOM_stoch_eos) and code clean-up (MOM_thickness_diffuse)

* Removed halo updates before determ param calc

* Update MOM_stoch_eos.F90

Removed unnecessary code

* Bug - indices are transposed

* Changed Stanley stochastic coefficient from exp(X) to exp(aX) (#9)

* Changed Stanley stochastic coefficient from exp(X) to exp(aX)

* Extra spaces removed

* Stoch eos init fix (#10)

* Don't bother calculating tv%varT if stanley_coeff<0

* Missing then added

* Merge Ian Grooms Tvar Discretization (#11)

* Update MOM_stoch_eos.F90

In progress updating stencil for$ | dx \times \nabla T|^2$ calculation

* New discretization of |dx\circ\nablaT|^2

Co-authored-by: Ian Grooms <ian.grooms@colorado.edu>

* Multiplied tvar%SGS by grid cell thickness ratio

* Added limiter for tv%varT

* Stoch eos ncar linear disc (#12)

* Update MOM_stoch_eos.F90

In progress updating stencil for$ | dx \times \nabla T|^2$ calculation

* New discretization of |dx\circ\nablaT|^2

* AR1 timescale land mask

Adds land mask to the computation of the AR1 decorrelation time

* Update dt in call to MOM_stoch_eos_run

The call to `MOM_stoch_eos_run` (which time steps the noise) is from within `step_MOM_dynamics`. `step_MOM_dynamics` advances on time step `dt` (per line 957), but the noise is updated using `dt_thermo`. It seems more appropriate to update the noise using `dt`, since it gets called from within `step_MOM_dynamics`.

* Fixed the units for r_sm_H

* Remove vestigial declarations

The variables `hl`, `Tl`, `mn_T`, `mn_T2`, and `r_sm_H` are no longer used, so I removed their declarations and an OMP private clause

Co-authored-by: Ian Grooms <ian.grooms@colorado.edu>

* Update MOM_thickness_diffuse.F90

Changed index for soft convention

* Update CVMix-src

* Ensure use_varT, etc., initialized

* Don't register stanley diagnostics if scheme is off

* Stanley density second derivs at h pts (#15)

* Change discretization of Stanley correction (drho_dT_dT at h points)

* Limit Stanley noise, shrink limiting value

* Revert t variance discretization

* Reverted variable declarations

* Stanley scheme in mixed_layer_restrat, vert_fill in stoch_eos, code cleanup (#19)

* Test Stanley EOS param in mixed_layer_restrat

* Fix size of TS cov, S var in Stanley calculate_density calls

* Test move stanley scheme initialization

* Added missing openMP directives

* Revert Stanley tvar discretization (#18)

* Perform vertical filling in calculation of T variance

* Variable declaration syntax error, remove scaling from get_param

* Fix call to vert_fill_TS

* Code cleanup, whitespace cleanup

Co-authored-by: Jessica Kenigson <jessicak@cheyenne1.cheyenne.ucar.edu>

* Use Stanley (2020) variance; scheme off at coast

* Comment clean-up

* Remove factor of 0.5 in Tvar

* Don't calculate Stanley diagnostics on halo

* Change start indices in stanley_density_1d

* Stanley param in MOM_isopycnal_slopes (#22)

Stanley param in MOM_isopycnal_slopes and thickness diffuse index fix

* Set eady flag to true if use_stored_slopes is true

* Cleanup, docs, whitespace

* Docs and whitespace

* Docs and whitespace

* Docs and whitespace

* Whitespace cleanup

* Whitespace cleanup

* Clean up whitespace

* Docs cleanup

* use_stanley

* Update MOM_lateral_mixing_coeffs.F90

* Adds link to another TEOS10 module

* Set Stanley off for testing

* Line continuation

Co-authored-by: Phil Pegion <38869668+pjpegion@users.noreply.github.com>
Co-authored-by: Philip Pegion <Philip.Pegion@noaa.gov>
Co-authored-by: Jessica Kenigson <jessicak@cheyenne6.cheyenne.ucar.edu>
Co-authored-by: Jessica Kenigson <jessicak@cheyenne3.cheyenne.ucar.edu>
Co-authored-by: jkenigson <jkenigso@gmail.com>
Co-authored-by: jskenigson <jessica.kenigson@colorado.edu>
Co-authored-by: Jessica Kenigson <jessicak@cheyenne1.cheyenne.ucar.edu>
Co-authored-by: Jessica Kenigson <jessicak@cheyenne5.cheyenne.ucar.edu>
Co-authored-by: Philip Pegion <ppegion@Philips-MacBook-Pro.local>
Co-authored-by: Jessica Kenigson <jessicak@cheyenne4.cheyenne.ucar.edu>
  • Loading branch information
11 people authored May 7, 2022
1 parent 2cd228f commit 8f97b3d
Show file tree
Hide file tree
Showing 13 changed files with 515 additions and 182 deletions.
4 changes: 2 additions & 2 deletions .testing/tc2/MOM_input
Original file line number Diff line number Diff line change
Expand Up @@ -297,7 +297,7 @@ BOUND_CORIOLIS = True ! [Boolean] default = False
! v-points, and similarly at v-points. This option would
! have no effect on the SADOURNY Coriolis scheme if it
! were possible to use centered difference thickness fluxes.
PGF_STANLEY_T2_DET_COEFF = 0.5 ! [nondim] default = -1.0
PGF_STANLEY_T2_DET_COEFF = -1.0 ! [nondim] default = -1.0
! The coefficient correlating SGS temperature variance with the mean temperature
! gradient in the deterministic part of the Stanley form of the Brankart
! correction. Negative values disable the scheme.
Expand Down Expand Up @@ -430,7 +430,7 @@ KHTH = 1.0 ! [m2 s-1] default = 0.0
! The background horizontal thickness diffusivity.
KHTH_MAX = 900.0 ! [m2 s-1] default = 0.0
! The maximum horizontal thickness diffusivity.
STANLEY_PRM_DET_COEFF = 0.5 ! [nondim] default = -1.0
STANLEY_PRM_DET_COEFF = -1.0 ! [nondim] default = -1.0
! The coefficient correlating SGS temperature variance with the mean temperature
! gradient in the deterministic part of the Stanley parameterization. Negative
! values disable the scheme.
Expand Down
22 changes: 21 additions & 1 deletion src/core/MOM.F90
Original file line number Diff line number Diff line change
Expand Up @@ -108,6 +108,7 @@ module MOM
use MOM_shared_initialization, only : write_ocean_geometry_file
use MOM_sponge, only : init_sponge_diags, sponge_CS
use MOM_state_initialization, only : MOM_initialize_state
use MOM_stoch_eos, only : MOM_stoch_eos_init,MOM_stoch_eos_run,MOM_stoch_eos_CS,mom_calc_varT
use MOM_sum_output, only : write_energy, accumulate_net_input
use MOM_sum_output, only : MOM_sum_output_init, MOM_sum_output_end
use MOM_sum_output, only : sum_output_CS
Expand Down Expand Up @@ -242,6 +243,7 @@ module MOM
logical :: use_ALE_algorithm !< If true, use the ALE algorithm rather than layered
!! isopycnal/stacked shallow water mode. This logical is set by calling the
!! function useRegridding() from the MOM_regridding module.
type(MOM_stoch_eos_CS) :: stoch_eos_CS !< structure containing random pattern for stoch EOS
logical :: alternate_first_direction !< If true, alternate whether the x- or y-direction
!! updates occur first in directionally split parts of the calculation.
real :: first_dir_restart = -1.0 !< A real copy of G%first_direction for use in restart files
Expand Down Expand Up @@ -444,6 +446,8 @@ module MOM
integer :: id_clock_other
integer :: id_clock_offline_tracer
integer :: id_clock_unit_tests
integer :: id_clock_stoch
integer :: id_clock_varT
!>@}

contains
Expand Down Expand Up @@ -1063,6 +1067,15 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, &
showCallTree = callTree_showQuery()

call cpu_clock_begin(id_clock_dynamics)
call cpu_clock_begin(id_clock_stoch)
if (CS%stoch_eos_CS%use_stoch_eos) call MOM_stoch_eos_run(G,u,v,dt,Time_local,CS%stoch_eos_CS,CS%diag)
call cpu_clock_end(id_clock_stoch)
call cpu_clock_begin(id_clock_varT)
if (CS%stoch_eos_CS%stanley_coeff >= 0.0) then
call MOM_calc_varT(G,GV,h,CS%tv,CS%stoch_eos_CS,dt)
call pass_var(CS%tv%varT, G%Domain,clock=id_clock_pass,halo=1)
endif
call cpu_clock_end(id_clock_varT)

if ((CS%t_dyn_rel_adv == 0.0) .and. CS%thickness_diffuse .and. CS%thickness_diffuse_first) then

Expand Down Expand Up @@ -1229,6 +1242,9 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, &
if (IDs%id_u > 0) call post_data(IDs%id_u, u, CS%diag)
if (IDs%id_v > 0) call post_data(IDs%id_v, v, CS%diag)
if (IDs%id_h > 0) call post_data(IDs%id_h, h, CS%diag)
if (CS%stoch_eos_CS%id_stoch_eos > 0) call post_data(CS%stoch_eos_CS%id_stoch_eos, CS%stoch_eos_CS%pattern, CS%diag)
if (CS%stoch_eos_CS%id_stoch_phi > 0) call post_data(CS%stoch_eos_CS%id_stoch_phi, CS%stoch_eos_CS%phi, CS%diag)
if (CS%stoch_eos_CS%id_tvar_sgs > 0) call post_data(CS%stoch_eos_CS%id_tvar_sgs, CS%tv%varT, CS%diag)
call disable_averaging(CS%diag)
call cpu_clock_end(id_clock_diagnostics) ; call cpu_clock_end(id_clock_other)

Expand Down Expand Up @@ -2765,6 +2781,8 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
call set_visc_init(Time, G, GV, US, param_file, diag, CS%visc, CS%set_visc_CSp, restart_CSp, CS%OBC)
call thickness_diffuse_init(Time, G, GV, US, param_file, diag, CS%CDp, CS%thickness_diffuse_CSp)

new_sim = is_new_run(restart_CSp)
call MOM_stoch_eos_init(G,Time,param_file,CS%stoch_eos_CS,restart_CSp,diag)
if (CS%split) then
allocate(eta(SZI_(G),SZJ_(G)), source=0.0)
call initialize_dyn_split_RK2(CS%u, CS%v, CS%h, CS%uh, CS%vh, eta, Time, &
Expand Down Expand Up @@ -2854,7 +2872,6 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
endif

! This subroutine initializes any tracer packages.
new_sim = is_new_run(restart_CSp)
call tracer_flow_control_init(.not.new_sim, Time, G, GV, US, CS%h, param_file, &
CS%diag, CS%OBC, CS%tracer_flow_CSp, CS%sponge_CSp, &
CS%ALE_sponge_CSp, CS%tv)
Expand Down Expand Up @@ -3065,6 +3082,7 @@ subroutine register_diags(Time, G, GV, US, IDs, diag)
v_extensive=.true.)
IDs%id_ssh_inst = register_diag_field('ocean_model', 'SSH_inst', diag%axesT1, &
Time, 'Instantaneous Sea Surface Height', 'm', conversion=US%Z_to_m)

end subroutine register_diags

!> Set up CPU clock IDs for timing various subroutines.
Expand Down Expand Up @@ -3097,6 +3115,8 @@ subroutine MOM_timing_init(CS)
if (CS%offline_tracer_mode) then
id_clock_offline_tracer = cpu_clock_id('Ocean offline tracers', grain=CLOCK_SUBCOMPONENT)
endif
id_clock_stoch = cpu_clock_id('(Stochastic EOS)', grain=CLOCK_MODULE)
id_clock_varT = cpu_clock_id('(SGS Temperature Variance)', grain=CLOCK_MODULE)

end subroutine MOM_timing_init

Expand Down
109 changes: 44 additions & 65 deletions src/core/MOM_PressureForce_FV.F90
Original file line number Diff line number Diff line change
Expand Up @@ -58,11 +58,12 @@ module MOM_PressureForce_FV
integer :: Recon_Scheme !< Order of the polynomial of the reconstruction of T & S
!! for the finite volume pressure gradient calculation.
!! By the default (1) is for a piecewise linear method
real :: Stanley_T2_det_coeff !< The coefficient correlating SGS temperature variance with
!! the mean temperature gradient in the deterministic part of
!! the Stanley form of the Brankart correction.

logical :: use_stanley_pgf !< If true, turn on Stanley parameterization in the PGF
integer :: id_e_tidal = -1 !< Diagnostic identifier
integer :: id_tvar_sgs = -1 !< Diagnostic identifier
integer :: id_rho_pgf = -1 !< Diagnostic identifier
integer :: id_rho_stanley_pgf = -1 !< Diagnostic identifier
integer :: id_p_stanley = -1 !< Diagnostic identifier
type(tidal_forcing_CS), pointer :: tides_CSp => NULL() !< Tides control structure
end type PressureForce_FV_CS

Expand Down Expand Up @@ -167,7 +168,7 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_
if (.not.CS%initialized) call MOM_error(FATAL, &
"MOM_PressureForce_FV_nonBouss: Module must be initialized before it is used.")

if (CS%Stanley_T2_det_coeff>=0.) call MOM_error(FATAL, &
if (CS%use_stanley_pgf) call MOM_error(FATAL, &
"MOM_PressureForce_FV_nonBouss: The Stanley parameterization is not yet"//&
"implemented in non-Boussinesq mode.")

Expand Down Expand Up @@ -473,6 +474,13 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
S_t, S_b, T_t, T_b ! Top and bottom edge values for linear reconstructions
! of salinity and temperature within each layer.
real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
rho_pgf, rho_stanley_pgf ! Density [kg m-3] from EOS with and without SGS T variance
! in Stanley parameterization.
real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
p_stanley ! Pressure [Pa] estimated with Rho_0
real :: rho_stanley_scalar ! Scalar quantity to hold density [kg m-3] in Stanley diagnostics.
real :: p_stanley_scalar ! Scalar quantity to hold pressure [Pa] in Stanley diagnostics.
real :: rho_in_situ(SZI_(G)) ! The in situ density [R ~> kg m-3].
real :: p_ref(SZI_(G)) ! The pressure used to calculate the coordinate
! density, [R L2 T-2 ~> Pa] (usually 2e7 Pa = 2000 dbar).
Expand All @@ -487,11 +495,6 @@ 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 :: 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 [degC2]
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 @@ -511,49 +514,6 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
use_ALE = .false.
if (associated(ALE_CSp)) use_ALE = CS%reconstruct .and. use_EOS

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
! 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)
enddo ; enddo ; enddo
endif

h_neglect = GV%H_subroundoff
dz_neglect = GV%H_subroundoff * GV%H_to_Z
I_Rho0 = 1.0 / GV%Rho0
Expand Down Expand Up @@ -690,7 +650,6 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
do k=1,nz
! Calculate 4 integrals through the layer that are required in the
! subsequent calculation.

if (use_EOS) then
! The following routine computes the integrals that are needed to
! calculate the pressure gradient force. Linear profiles for T and S are
Expand All @@ -701,13 +660,13 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
if ( CS%Recon_Scheme == 1 ) then
call int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, &
rho_ref, CS%Rho0, GV%g_Earth, dz_neglect, G%bathyT, &
G%HI, GV, tv%eqn_of_state, US, dpa, intz_dpa, intx_dpa, inty_dpa, &
G%HI, GV, tv%eqn_of_state, US, CS%use_stanley_pgf, dpa, intz_dpa, intx_dpa, inty_dpa, &
useMassWghtInterp=CS%useMassWghtInterp, &
use_inaccurate_form=CS%use_inaccurate_pgf_rho_anom, Z_0p=G%Z_ref)
elseif ( CS%Recon_Scheme == 2 ) then
call int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, &
rho_ref, CS%Rho0, GV%g_Earth, dz_neglect, G%bathyT, &
G%HI, GV, tv%eqn_of_state, US, dpa, intz_dpa, intx_dpa, inty_dpa, &
G%HI, GV, tv%eqn_of_state, US, CS%use_stanley_pgf, dpa, intz_dpa, intx_dpa, inty_dpa, &
useMassWghtInterp=CS%useMassWghtInterp, Z_0p=G%Z_ref)
endif
else
Expand Down Expand Up @@ -797,8 +756,26 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
endif
endif

if (CS%use_stanley_pgf) then
do j=js,je ; do i=is,ie ;
p_stanley_scalar=0.0
do k=1, nz
p_stanley_scalar = p_stanley_scalar + 0.5 * h(i,j,k) * GV%H_to_Pa !Pressure at mid-point of layer
call calculate_density(tv%T(i,j,k), tv%S(i,j,k), p_stanley_scalar, 0.0, 0.0, 0.0, &
rho_stanley_scalar, tv%eqn_of_state)
rho_pgf(i,j,k) = rho_stanley_scalar
call calculate_density(tv%T(i,j,k), tv%S(i,j,k), p_stanley_scalar, tv%varT(i,j,k), 0.0, 0.0, &
rho_stanley_scalar, tv%eqn_of_state)
rho_stanley_pgf(i,j,k) = rho_stanley_scalar
p_stanley(i,j,k) = p_stanley_scalar
p_stanley_scalar = p_stanley_scalar + 0.5 * h(i,j,k) * GV%H_to_Pa !Pressure at bottom of layer
enddo; enddo; enddo
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)
if (CS%id_rho_pgf>0) call post_data(CS%id_rho_pgf, rho_pgf, CS%diag)
if (CS%id_rho_stanley_pgf>0) call post_data(CS%id_rho_stanley_pgf, rho_stanley_pgf, CS%diag)
if (CS%id_p_stanley>0) call post_data(CS%id_p_stanley, p_stanley, CS%diag)

end subroutine PressureForce_FV_Bouss

Expand Down Expand Up @@ -859,14 +836,16 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, tides_CS
"boundary cells is extrapolated, rather than using PCM "//&
"in these cells. If true, the same order polynomial is "//&
"used as is used for the interior cells.", default=.true.)
call get_param(param_file, mdl, "PGF_STANLEY_T2_DET_COEFF", CS%Stanley_T2_det_coeff, &
"The coefficient correlating SGS temperature variance with "// &
"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')
call get_param(param_file, mdl, "USE_STANLEY_PGF", CS%use_stanley_pgf, &
"If true, turn on Stanley SGS T variance parameterization "// &
"in PGF code.", default=.false.)
if (CS%use_stanley_pgf) then
CS%id_rho_pgf = register_diag_field('ocean_model', 'rho_pgf', diag%axesTL, &
Time, 'rho in PGF', 'kg m3')
CS%id_rho_stanley_pgf = register_diag_field('ocean_model', 'rho_stanley_pgf', diag%axesTL, &
Time, 'rho in PGF with Stanley correction', 'kg m3')
CS%id_p_stanley = register_diag_field('ocean_model', 'p_stanley', diag%axesTL, &
Time, 'p in PGF with Stanley correction', 'Pa')
endif
if (CS%tides) then
CS%id_e_tidal = register_diag_field('ocean_model', 'e_tidal', diag%axesT1, &
Expand Down
Loading

0 comments on commit 8f97b3d

Please sign in to comment.