Stoch eos

.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

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,12 @@ 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 :: 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).
@@ -796,8 +804,27 @@ 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)))
+      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
+      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
+   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 +896,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, &

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,33 +436,57 @@ 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
   if (tv%S_is_absS) then
     ! 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', &