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/core/MOM_PressureForce_FV.F90 b/src/core/MOM_PressureForce_FV.F90 index de268a27a8..ea2e400f7e 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,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, & 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', &