From 9a6b5faef156b3c886084195b1db68b8e166d5c0 Mon Sep 17 00:00:00 2001 From: Jessica Kenigson Date: Wed, 7 Oct 2020 11:09:21 -0600 Subject: [PATCH 1/5] Added necessary submodules and S^2, T^2 diagnostics to MOM_diagnostics --- .gitmodules | 6 +++++ pkg/geoKdTree | 1 + pkg/mom6_da_hooks | 1 + src/diagnostics/MOM_diagnostics.F90 | 39 ++++++++++++++++++++++++++--- 4 files changed, 44 insertions(+), 3 deletions(-) create mode 160000 pkg/geoKdTree create mode 160000 pkg/mom6_da_hooks diff --git a/.gitmodules b/.gitmodules index 637f1188ed..1ad22f9f7f 100644 --- a/.gitmodules +++ b/.gitmodules @@ -4,3 +4,9 @@ [submodule "pkg/GSW-Fortran"] path = pkg/GSW-Fortran url = https://github.com/TEOS-10/GSW-Fortran.git +[submodule "pkg/mom6_da_hooks"] + path = pkg/mom6_da_hooks + url = https://github.com/NOAA-GFDL/MOM6_DA_hooks.git +[submodule "pkg/geoKdTree"] + path = pkg/geoKdTree + url = https://github.com/travissluka/geoKdTree.git diff --git a/pkg/geoKdTree b/pkg/geoKdTree new file mode 160000 index 0000000000..f8ac844ac5 --- /dev/null +++ b/pkg/geoKdTree @@ -0,0 +1 @@ +Subproject commit f8ac844ac558979e43697a6f5e7d9305efea088e diff --git a/pkg/mom6_da_hooks b/pkg/mom6_da_hooks new file mode 160000 index 0000000000..9c930afc5e --- /dev/null +++ b/pkg/mom6_da_hooks @@ -0,0 +1 @@ +Subproject commit 9c930afc5e2c4f86085476f524fc71dec321f68b diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 3936a788d0..cc143b299e 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -142,6 +142,8 @@ module MOM_diagnostics integer :: id_rhopot0 = -1, id_rhopot2 = -1 integer :: id_drho_dT = -1, id_drho_dS = -1 integer :: id_h_pre_sync = -1 + integer :: id_tosq = -1, id_sosq = -1 + !>@} !> The control structure for calculating wave speed. type(wave_speed_CS), pointer :: wave_speed_CSp => NULL() @@ -434,16 +436,28 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & ! Internal T&S variables are conservative temperature & absolute salinity, ! so they need to converted to potential temperature and practical salinity ! for some diagnostics using TEOS-10 function calls. - if ((CS%id_Tpot > 0) .or. (CS%id_tob > 0)) then + if ((CS%id_Tpot > 0) .or. (CS%id_tob > 0) .or. (CS%id_tosq > 0)) then do k=1,nz ; do j=js,je ; do i=is,ie work_3d(i,j,k) = gsw_pt_from_ct(tv%S(i,j,k),tv%T(i,j,k)) enddo ; enddo ; enddo if (CS%id_Tpot > 0) call post_data(CS%id_Tpot, work_3d, CS%diag) if (CS%id_tob > 0) call post_data(CS%id_tob, work_3d(:,:,nz), CS%diag, mask=G%mask2dT) + if (CS%id_tosq > 0) then + do k=1,nz ; do j=js,je ; do i=is,ie + work_3d(i,j,k) = work_3d(i,j,k)*work_3d(i,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_tosq, work_3d, CS%diag) + endif endif else ! Internal T&S variables are potential temperature & practical salinity if (CS%id_tob > 0) call post_data(CS%id_tob, tv%T(:,:,nz), CS%diag, mask=G%mask2dT) + if (CS%id_tosq > 0) then + do k=1,nz ; do j=js,je ; do i=is,ie + work_3d(i,j,k) = tv%T(i,j,k)*tv%T(i,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_tosq, work_3d, CS%diag) + endif endif ! Calculate additional, potentially derived salinity diagnostics @@ -451,16 +465,28 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & ! Internal T&S variables are conservative temperature & absolute salinity, ! so they need to converted to potential temperature and practical salinity ! for some diagnostics using TEOS-10 function calls. - if ((CS%id_Sprac > 0) .or. (CS%id_sob > 0)) then + if ((CS%id_Sprac > 0) .or. (CS%id_sob > 0) .or. (CS%id_sosq >0)) then do k=1,nz ; do j=js,je ; do i=is,ie work_3d(i,j,k) = gsw_sp_from_sr(tv%S(i,j,k)) enddo ; enddo ; enddo if (CS%id_Sprac > 0) call post_data(CS%id_Sprac, work_3d, CS%diag) - if (CS%id_sob > 0) call post_data(CS%id_sob, work_3d(:,:,nz), CS%diag, mask=G%mask2dT) + if (CS%id_sob > 0) call post_data(CS%id_sob, work_3d(:,:,nz), CS%diag, mask=G%mask2dT) + if (CS%id_sosq > 0) then + do k=1,nz ; do j=js,je ; do i=is,ie + work_3d(i,j,k) = work_3d(i,j,k)*work_3d(i,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_sosq, work_3d, CS%diag) + endif endif else ! Internal T&S variables are potential temperature & practical salinity if (CS%id_sob > 0) call post_data(CS%id_sob, tv%S(:,:,nz), CS%diag, mask=G%mask2dT) + if (CS%id_sosq > 0) then + do k=1,nz ; do j=js,je ; do i=is,ie + work_3d(i,j,k) = tv%S(i,j,k)*tv%S(i,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_sosq, work_3d, CS%diag) + endif endif ! volume mean potential temperature @@ -1611,6 +1637,13 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag long_name='Sea Water Salinity at Sea Floor', & standard_name='sea_water_salinity_at_sea_floor', units='psu') + CS%id_tosq = register_diag_field('ocean_model', 'tosq', diag%axesTL,& + Time, 'Square of Potential Temperature', 'degc2', & + standard_name='Potential Temperature Squared') + CS%id_sosq = register_diag_field('ocean_model', 'sosq', diag%axesTL,& + Time, 'Square of Salinity', 'psu2', & + standard_name='Salinity Squared') + CS%id_temp_layer_ave = register_diag_field('ocean_model', 'temp_layer_ave', & diag%axesZL, Time, 'Layer Average Ocean Temperature', 'degC') CS%id_salt_layer_ave = register_diag_field('ocean_model', 'salt_layer_ave', & From 9f649411cc55885c7a702f8cf546e789be6259f8 Mon Sep 17 00:00:00 2001 From: Jessica Kenigson Date: Wed, 21 Oct 2020 08:27:20 -0600 Subject: [PATCH 2/5] Added diagnostics for outputting variables related to the stochastic parameterization. --- src/core/MOM_PressureForce_FV.F90 | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/src/core/MOM_PressureForce_FV.F90 b/src/core/MOM_PressureForce_FV.F90 index de268a27a8..0957c50c64 100644 --- a/src/core/MOM_PressureForce_FV.F90 +++ b/src/core/MOM_PressureForce_FV.F90 @@ -59,6 +59,8 @@ module MOM_PressureForce_FV !! the Stanley form of the Brankart correction. 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 type(tidal_forcing_CS), pointer :: tides_CSp => NULL() !< Tides control structure end type PressureForce_FV_CS @@ -467,6 +469,11 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: & 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)) :: & + p_stanley, rho_pgf, rho_stanley_pgf ! Pressure [Pa] estimated with Rho_0 and + ! density [kg m-3] from EOS with and without SGS T variance + ! in Stanley parameterization. + real :: rho_stanley_scalar ! Scalar quantity to hold density [kg m-3] 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). @@ -796,8 +803,22 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm endif endif + if (CS%Stanley_T2_det_coeff>=0.) then + do k=1, nz ; do j=js-1,je+1 ; do i=is-1,ie+1 + p_stanley(i,j,k)= -1 * rho_ref * GV%g_Earth * e(i,j,k) + call calculate_density(tv%T(i,j,k), tv%S(i,j,k), p_stanley(i,j,k), 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(i,j,k), 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 + 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_tvar_sgs>0) call post_data(CS%id_rho_pgf, rho_pgf, CS%diag) + if (CS%id_tvar_sgs>0) call post_data(CS%id_rho_stanley_pgf, rho_stanley_pgf, CS%diag) end subroutine PressureForce_FV_Bouss @@ -869,6 +890,10 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, tides_CS 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') + 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') endif if (CS%tides) then CS%id_e_tidal = register_diag_field('ocean_model', 'e_tidal', diag%axesT1, & From 363871a0de142abba929c0226d7e7ad64f513ab5 Mon Sep 17 00:00:00 2001 From: Jessica Kenigson Date: Thu, 22 Oct 2020 14:34:56 -0600 Subject: [PATCH 3/5] Diagnostics in MOM_PressureForce_FV updated for stochastic (rather than deterministic) Stanley SGS T variance parameterization. --- src/core/MOM_PressureForce_FV.F90 | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/core/MOM_PressureForce_FV.F90 b/src/core/MOM_PressureForce_FV.F90 index 0957c50c64..20ea54d56c 100644 --- a/src/core/MOM_PressureForce_FV.F90 +++ b/src/core/MOM_PressureForce_FV.F90 @@ -473,7 +473,8 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm p_stanley, rho_pgf, rho_stanley_pgf ! Pressure [Pa] estimated with Rho_0 and ! density [kg m-3] from EOS with and without SGS T variance ! in Stanley parameterization. - real :: rho_stanley_scalar ! Scalar quantity to hold density [kg m-3] in Stanley diagnostics. + real :: rho_stanley_scalar ! Scalar quantity to hold density [kg m-3] in Stanley diagnostics. + real :: tv_stanley_scalar ! Scalar quantity to hold SGS T variance [degc2] 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). @@ -806,10 +807,11 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm if (CS%Stanley_T2_det_coeff>=0.) then do k=1, nz ; do j=js-1,je+1 ; do i=is-1,ie+1 p_stanley(i,j,k)= -1 * rho_ref * GV%g_Earth * e(i,j,k) + tv_stanley_scalar=exp(stoch_eos_pattern(i,j))*tv%varT(i,j,k) call calculate_density(tv%T(i,j,k), tv%S(i,j,k), p_stanley(i,j,k), 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(i,j,k), tv%varT(i,j,k), 0.0, 0.0, & + call calculate_density(tv%T(i,j,k), tv%S(i,j,k), p_stanley(i,j,k), tv_stanley_scalar, 0.0, 0.0, & rho_stanley_scalar, tv%eqn_of_state) rho_stanley_pgf(i,j,k)=rho_stanley_scalar enddo; enddo; enddo From d039cd60ff66ef55b81d92bd72fdcd5133772fef Mon Sep 17 00:00:00 2001 From: Jessica Kenigson Date: Thu, 29 Oct 2020 14:38:27 -0600 Subject: [PATCH 4/5] Added parentheses for reproducibility. --- src/core/MOM_PressureForce_FV.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/core/MOM_PressureForce_FV.F90 b/src/core/MOM_PressureForce_FV.F90 index 20ea54d56c..1c101a7e38 100644 --- a/src/core/MOM_PressureForce_FV.F90 +++ b/src/core/MOM_PressureForce_FV.F90 @@ -806,7 +806,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm if (CS%Stanley_T2_det_coeff>=0.) then do k=1, nz ; do j=js-1,je+1 ; do i=is-1,ie+1 - p_stanley(i,j,k)= -1 * rho_ref * GV%g_Earth * e(i,j,k) + p_stanley(i,j,k)= -1 * (rho_ref * (GV%g_Earth * e(i,j,k))) tv_stanley_scalar=exp(stoch_eos_pattern(i,j))*tv%varT(i,j,k) call calculate_density(tv%T(i,j,k), tv%S(i,j,k), p_stanley(i,j,k), 0.0, 0.0, 0.0, & rho_stanley_scalar, tv%eqn_of_state) From c6a12ce0b14cee634c36f19c1a5308319ab973ab Mon Sep 17 00:00:00 2001 From: Jessica Kenigson Date: Thu, 5 Nov 2020 10:16:21 -0700 Subject: [PATCH 5/5] Changed diagnostics to account for possible absence of stoch_eos_pattern in MOM_PressureForce_FV, when deterministic parameterization is on. --- src/core/MOM_PressureForce_FV.F90 | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/core/MOM_PressureForce_FV.F90 b/src/core/MOM_PressureForce_FV.F90 index 1c101a7e38..ea2e400f7e 100644 --- a/src/core/MOM_PressureForce_FV.F90 +++ b/src/core/MOM_PressureForce_FV.F90 @@ -807,7 +807,11 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm if (CS%Stanley_T2_det_coeff>=0.) then do k=1, nz ; do j=js-1,je+1 ; do i=is-1,ie+1 p_stanley(i,j,k)= -1 * (rho_ref * (GV%g_Earth * e(i,j,k))) - tv_stanley_scalar=exp(stoch_eos_pattern(i,j))*tv%varT(i,j,k) + if (present(stoch_eos_pattern)) then + tv_stanley_scalar=exp(stoch_eos_pattern(i,j))*tv%varT(i,j,k) + else + tv_stanley_scalar=tv%varT(i,j,k) + endif call calculate_density(tv%T(i,j,k), tv%S(i,j,k), p_stanley(i,j,k), 0.0, 0.0, 0.0, & rho_stanley_scalar, tv%eqn_of_state) rho_pgf(i,j,k)=rho_stanley_scalar