diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index c61ed72e0c..af8481fd1c 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -91,7 +91,7 @@ module MOM use MOM_grid, only : set_first_direction, rescale_grid_bathymetry use MOM_hor_index, only : hor_index_type, hor_index_init use MOM_hor_index, only : rotate_hor_index -use MOM_interface_heights, only : find_eta +use MOM_interface_heights, only : find_eta, calc_derived_thermo use MOM_interface_filter, only : interface_filter, interface_filter_init, interface_filter_end use MOM_interface_filter, only : interface_filter_CS use MOM_lateral_mixing_coeffs, only : calc_slope_functions, VarMix_init, VarMix_end @@ -1400,6 +1400,12 @@ subroutine step_MOM_tracer_dyn(CS, G, GV, US, h, Time_local) call create_group_pass(pass_T_S, CS%tv%T, G%Domain, To_All+Omit_Corners, halo=1) call create_group_pass(pass_T_S, CS%tv%S, G%Domain, To_All+Omit_Corners, halo=1) call do_group_pass(pass_T_S, G%Domain, clock=id_clock_pass) + halo_sz = 1 + endif + + ! Update derived thermodynamic quantities. + if (allocated(CS%tv%SpV_avg)) then + call calc_derived_thermo(CS%tv, h, G, GV, US, halo=halo_sz) endif endif @@ -1581,6 +1587,11 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, & call create_group_pass(pass_uv_T_S_h, h, G%Domain, halo=dynamics_stencil) call do_group_pass(pass_uv_T_S_h, G%Domain, clock=id_clock_pass) + ! Update derived thermodynamic quantities. + if (allocated(tv%SpV_avg)) then + call calc_derived_thermo(tv, h, G, GV, US, halo=dynamics_stencil) + endif + if (CS%debug .and. CS%use_ALE_algorithm) then call MOM_state_chksum("Post-ALE ", u, v, h, CS%uh, CS%vh, G, GV, US) call hchksum(tv%T, "Post-ALE T", G%HI, haloshift=1, scale=US%C_to_degC) @@ -1623,13 +1634,19 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, & call cpu_clock_end(id_clock_adiabatic) if (associated(tv%T)) then - call create_group_pass(pass_T_S, tv%T, G%Domain, To_All+Omit_Corners, halo=1) - call create_group_pass(pass_T_S, tv%S, G%Domain, To_All+Omit_Corners, halo=1) + dynamics_stencil = min(3, G%Domain%nihalo, G%Domain%njhalo) + call create_group_pass(pass_T_S, tv%T, G%Domain, To_All+Omit_Corners, halo=dynamics_stencil) + call create_group_pass(pass_T_S, tv%S, G%Domain, To_All+Omit_Corners, halo=dynamics_stencil) call do_group_pass(pass_T_S, G%Domain, clock=id_clock_pass) if (CS%debug) then if (associated(tv%T)) call hchksum(tv%T, "Post-diabatic T", G%HI, haloshift=1, scale=US%C_to_degC) if (associated(tv%S)) call hchksum(tv%S, "Post-diabatic S", G%HI, haloshift=1, scale=US%S_to_ppt) endif + + ! Update derived thermodynamic quantities. + if (allocated(tv%SpV_avg)) then + call calc_derived_thermo(tv, h, G, GV, US, halo=dynamics_stencil) + endif endif endif ! endif for the block "if (.not.CS%adiabatic)" @@ -1676,6 +1693,8 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS type(time_type), pointer :: accumulated_time => NULL() type(time_type), pointer :: vertical_time => NULL() + integer :: dynamics_stencil ! The computational stencil for the calculations + ! in the dynamic core. integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz ! 3D pointers @@ -1848,6 +1867,12 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS fluxes%fluxes_used = .true. + ! Update derived thermodynamic quantities. + if (allocated(CS%tv%SpV_avg)) then + dynamics_stencil = min(3, G%Domain%nihalo, G%Domain%njhalo) + call calc_derived_thermo(CS%tv, CS%h, G, GV, US, halo=dynamics_stencil) + endif + if (last_iter) then accumulated_time = real_to_time(0.0) endif @@ -2817,6 +2842,11 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & endif endif + ! Allocate any derived equation of state fields. + if (use_temperature .and. .not.(GV%Boussinesq .or. GV%semi_Boussinesq)) then + allocate(CS%tv%SpV_avg(isd:ied,jsd:jed,nz), source=0.0) + endif + if (use_ice_shelf .and. CS%debug) then call hchksum(CS%frac_shelf_h, "MOM:frac_shelf_h", G%HI, haloshift=0) call hchksum(CS%mass_shelf, "MOM:mass_shelf", G%HI, haloshift=0,scale=US%RZ_to_kg_m2) @@ -3103,6 +3133,11 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & call do_group_pass(pass_uv_T_S_h, G%Domain) + ! Update derived thermodynamic quantities. + if (allocated(CS%tv%SpV_avg)) then + call calc_derived_thermo(CS%tv, CS%h, G, GV, US, halo=dynamics_stencil) + endif + if (associated(CS%visc%Kv_shear)) & call pass_var(CS%visc%Kv_shear, G%Domain, To_All+Omit_Corners, halo=1) @@ -3931,6 +3966,7 @@ subroutine MOM_end(CS) if (associated(CS%Hml)) deallocate(CS%Hml) if (associated(CS%tv%salt_deficit)) deallocate(CS%tv%salt_deficit) if (associated(CS%tv%frazil)) deallocate(CS%tv%frazil) + if (allocated(CS%tv%SpV_avg)) deallocate(CS%tv%SpV_avg) if (associated(CS%tv%T)) then DEALLOC_(CS%T) ; CS%tv%T => NULL() ; DEALLOC_(CS%S) ; CS%tv%S => NULL() diff --git a/src/core/MOM_density_integrals.F90 b/src/core/MOM_density_integrals.F90 index e1fb3d3278..9fed528e71 100644 --- a/src/core/MOM_density_integrals.F90 +++ b/src/core/MOM_density_integrals.F90 @@ -4,12 +4,13 @@ module MOM_density_integrals ! This file is part of MOM6. See LICENSE.md for the license. use MOM_EOS, only : EOS_type -use MOM_EOS, only : EOS_quadrature +use MOM_EOS, only : EOS_quadrature, EOS_domain use MOM_EOS, only : analytic_int_density_dz use MOM_EOS, only : analytic_int_specific_vol_dp use MOM_EOS, only : calculate_density use MOM_EOS, only : calculate_spec_vol use MOM_EOS, only : calculate_specific_vol_derivs +use MOM_EOS, only : average_specific_vol use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg use MOM_hor_index, only : hor_index_type use MOM_string_functions, only : uppercase @@ -28,6 +29,7 @@ module MOM_density_integrals public int_specific_vol_dp public int_spec_vol_dp_generic_pcm public int_spec_vol_dp_generic_plm +public avg_specific_vol public find_depth_of_pressure_in_cell contains @@ -1613,6 +1615,36 @@ subroutine find_depth_of_pressure_in_cell(T_t, T_b, S_t, S_b, z_t, z_b, P_t, P_t end subroutine find_depth_of_pressure_in_cell +!> Calculate the average in situ specific volume across layers +subroutine avg_specific_vol(T, S, p_t, dp, HI, EOS, SpV_avg, halo_size) + type(hor_index_type), intent(in) :: HI !< The horizontal index structure + real, dimension(SZI_(HI),SZJ_(HI)), & + intent(in) :: T !< Potential temperature of the layer [C ~> degC] + real, dimension(SZI_(HI),SZJ_(HI)), & + intent(in) :: S !< Salinity of the layer [S ~> ppt] + real, dimension(SZI_(HI),SZJ_(HI)), & + intent(in) :: p_t !< Pressure at the top of the layer [R L2 T-2 ~> Pa] + real, dimension(SZI_(HI),SZJ_(HI)), & + intent(in) :: dp !< Pressure change in the layer [R L2 T-2 ~> Pa] + type(EOS_type), intent(in) :: EOS !< Equation of state structure + real, dimension(SZI_(HI),SZJ_(HI)), & + intent(inout) :: SpV_avg !< The vertical average specific volume + !! in the layer [R-1 ~> m3 kg-1] + integer, optional, intent(in) :: halo_size !< The number of halo points in which to work. + + ! Local variables + integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state + integer :: jsh, jeh, j, halo + + halo = 0 ; if (present(halo_size)) halo = MAX(halo_size,0) + jsh = HI%jsc-halo ; jeh = HI%jec+halo + + EOSdom(:) = EOS_domain(HI, halo_size) + do j=jsh,jeh + call average_specific_vol(T(:,j), S(:,j), p_t(:,j), dp(:,j), SpV_avg(:,j), EOS, EOSdom) + enddo + +end subroutine avg_specific_vol !> Returns change in anomalous pressure change from top to non-dimensional !! position pos between z_t and z_b [R L2 T-2 ~> Pa] diff --git a/src/core/MOM_interface_heights.F90 b/src/core/MOM_interface_heights.F90 index 4f41cb074b..befeb1c2ad 100644 --- a/src/core/MOM_interface_heights.F90 +++ b/src/core/MOM_interface_heights.F90 @@ -3,7 +3,7 @@ module MOM_interface_heights ! This file is part of MOM6. See LICENSE.md for the license. -use MOM_density_integrals, only : int_specific_vol_dp +use MOM_density_integrals, only : int_specific_vol_dp, avg_specific_vol use MOM_error_handler, only : MOM_error, FATAL use MOM_EOS, only : calculate_density, EOS_type, EOS_domain use MOM_file_parser, only : log_version @@ -16,18 +16,26 @@ module MOM_interface_heights #include -public find_eta, dz_to_thickness, dz_to_thickness_simple +public find_eta, dz_to_thickness, thickness_to_dz, dz_to_thickness_simple +public calc_derived_thermo !> Calculates the heights of the free surface or all interfaces from layer thicknesses. interface find_eta module procedure find_eta_2d, find_eta_3d end interface find_eta -!> Calculates layer thickness in thickness units from geometric thicknesses in height units. +!> Calculates layer thickness in thickness units from geometric distance between the +!! interfaces around that layer in height units. interface dz_to_thickness module procedure dz_to_thickness_tv, dz_to_thickness_EoS end interface dz_to_thickness +!> Converts layer thickness in thickness units into the vertical distance between the +!! interfaces around a layer in height units. +interface thickness_to_dz + module procedure thickness_to_dz_3d, thickness_to_dz_jslice +end interface thickness_to_dz + contains !> Calculates the heights of all interfaces between layers, using the appropriate @@ -253,6 +261,45 @@ subroutine find_eta_2d(h, tv, G, GV, US, eta, eta_bt, halo_size, dZref) end subroutine find_eta_2d +!> Calculate derived thermodynamic quantities for re-use later. +subroutine calc_derived_thermo(tv, h, G, GV, US, halo) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(thermo_var_ptrs), intent(inout) :: tv !< A structure pointing to various + !! thermodynamic variables, some of + !! which will be set here. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]. + integer, optional, intent(in) :: halo !< Width of halo within which to + !! calculate thicknesses + ! Local variables + real, dimension(SZI_(G),SZJ_(G)) :: p_t ! Hydrostatic pressure atop a layer [R L2 T-2 ~> Pa] + real, dimension(SZI_(G),SZJ_(G)) :: dp ! Pressure change across a layer [R L2 T-2 ~> Pa] + integer :: i, j, k, is, ie, js, je, halos, nz + + halos = 0 ; if (present(halo)) halos = max(0,halo) + is = G%isc-halos ; ie = G%iec+halos ; js = G%jsc-halos ; je = G%jec+halos ; nz = GV%ke + + if (allocated(tv%Spv_avg) .and. associated(tv%eqn_of_state)) then + if (associated(tv%p_surf)) then + do j=js,je ; do i=is,ie ; p_t(i,j) = tv%p_surf(i,j) ; enddo ; enddo + else + do j=js,je ; do i=is,ie ; p_t(i,j) = 0.0 ; enddo ; enddo + endif + do k=1,nz + do j=js,je ; do i=is,ie + dp(i,j) = GV%g_Earth*GV%H_to_RZ*h(i,j,k) + enddo ; enddo + call avg_specific_vol(tv%T(:,:,k), tv%S(:,:,k), p_t, dp, G%HI, tv%eqn_of_state, tv%SpV_avg(:,:,k), halo) + if (k Converts thickness from geometric height units to thickness units, perhaps via an !! inversion of the integral of the density in pressure using variables stored in !! the thermo_var_ptrs type when in non-Boussinesq mode. @@ -428,4 +475,75 @@ subroutine dz_to_thickness_simple(dz, h, G, GV, US, halo_size, layer_mode) end subroutine dz_to_thickness_simple +!> Converts layer thicknesses in thickness units to the vertical distance between edges in height +!! units, perhaps by multiplication by the precomputed layer-mean specific volume stored in an +!! array in the thermo_var_ptrs type when in non-Boussinesq mode. +subroutine thickness_to_dz_3d(h, tv, dz, G, GV, US, halo_size) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: h !< Input thicknesses in thickness units [H ~> m or kg m-2]. + type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various + !! thermodynamic variables + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(inout) :: dz !< Geometric layer thicknesses in height units [Z ~> m] + !! This is essentially intent out, but declared as intent + !! inout to preserve any initialized values in halo points. + integer, optional, intent(in) :: halo_size !< Width of halo within which to + !! calculate thicknesses + ! Local variables + integer :: i, j, k, is, ie, js, je, halo, nz + + halo = 0 ; if (present(halo_size)) halo = max(0,halo_size) + is = G%isc-halo ; ie = G%iec+halo ; js = G%jsc-halo ; je = G%jec+halo ; nz = GV%ke + + if ((.not.GV%Boussinesq) .and. allocated(tv%SpV_avg)) then + do k=1,nz ; do j=js,je ; do i=is,ie + dz(i,j,k) = GV%H_to_RZ * h(i,j,k) * tv%SpV_avg(i,j,k) + enddo ; enddo ; enddo + else + do k=1,nz ; do j=js,je ; do i=is,ie + dz(i,j,k) = GV%H_to_Z * h(i,j,k) + enddo ; enddo ; enddo + endif + +end subroutine thickness_to_dz_3d + + +!> Converts a vertical i- / k- slice of layer thicknesses in thickness units to the vertical +!! distance between edges in height units, perhaps by multiplication by the precomputed layer-mean +!! specific volume stored in an array in the thermo_var_ptrs type when in non-Boussinesq mode. +subroutine thickness_to_dz_jslice(h, tv, dz, j, G, GV, halo_size) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: h !< Input thicknesses in thickness units [H ~> m or kg m-2]. + type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various + !! thermodynamic variables + real, dimension(SZI_(G),SZK_(GV)), & + intent(inout) :: dz !< Geometric layer thicknesses in height units [Z ~> m] + !! This is essentially intent out, but declared as intent + !! inout to preserve any initialized values in halo points. + integer, intent(in) :: j !< The second (j-) index of the input thicknesses to work with + integer, optional, intent(in) :: halo_size !< Width of halo within which to + !! calculate thicknesses + ! Local variables + integer :: i, k, is, ie, halo, nz + + halo = 0 ; if (present(halo_size)) halo = max(0,halo_size) + is = G%isc-halo ; ie = G%iec+halo ; nz = GV%ke + + if ((.not.GV%Boussinesq) .and. allocated(tv%SpV_avg)) then + do k=1,nz ; do i=is,ie + dz(i,k) = GV%H_to_RZ * h(i,j,k) * tv%SpV_avg(i,j,k) + enddo ; enddo + else + do k=1,nz ; do i=is,ie + dz(i,k) = GV%H_to_Z * h(i,j,k) + enddo ; enddo + endif + +end subroutine thickness_to_dz_jslice + end module MOM_interface_heights diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index 5efb02fe44..bec93376af 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -93,6 +93,9 @@ module MOM_variables logical :: S_is_absS = .false. !< If true, the salinity variable tv%S is !! actually the absolute salinity in units of [gSalt kg-1]. real :: min_salinity !< The minimum value of salinity when BOUND_SALINITY=True [S ~> ppt]. + real, allocatable, dimension(:,:,:) :: SpV_avg + !< The layer averaged in situ specific volume [R-1 ~> m3 kg-1]. + ! These arrays are accumulated fluxes for communication with other components. real, dimension(:,:), pointer :: frazil => NULL() !< The energy needed to heat the ocean column to the diff --git a/src/equation_of_state/MOM_EOS.F90 b/src/equation_of_state/MOM_EOS.F90 index 276c4c3019..c68dc7b661 100644 --- a/src/equation_of_state/MOM_EOS.F90 +++ b/src/equation_of_state/MOM_EOS.F90 @@ -8,24 +8,25 @@ module MOM_EOS use MOM_EOS_linear, only : calculate_specvol_derivs_linear, int_density_dz_linear use MOM_EOS_linear, only : calculate_density_second_derivs_linear, EoS_fit_range_linear use MOM_EOS_linear, only : calculate_compress_linear, int_spec_vol_dp_linear +use MOM_EOS_linear, only : avg_spec_vol_linear use MOM_EOS_Wright, only : calculate_density_wright, calculate_spec_vol_wright use MOM_EOS_Wright, only : calculate_density_derivs_wright use MOM_EOS_Wright, only : calculate_specvol_derivs_wright, int_density_dz_wright use MOM_EOS_Wright, only : calculate_compress_wright, int_spec_vol_dp_wright use MOM_EOS_Wright, only : calculate_density_second_derivs_wright, calc_density_second_derivs_wright_buggy -use MOM_EOS_Wright, only : EoS_fit_range_Wright +use MOM_EOS_Wright, only : EoS_fit_range_Wright, avg_spec_vol_Wright use MOM_EOS_Wright_full, only : calculate_density_wright_full, calculate_spec_vol_wright_full use MOM_EOS_Wright_full, only : calculate_density_derivs_wright_full use MOM_EOS_Wright_full, only : calculate_specvol_derivs_wright_full, int_density_dz_wright_full use MOM_EOS_Wright_full, only : calculate_compress_wright_full, int_spec_vol_dp_wright_full use MOM_EOS_Wright_full, only : calculate_density_second_derivs_wright_full -use MOM_EOS_Wright_full, only : EoS_fit_range_Wright_full +use MOM_EOS_Wright_full, only : EoS_fit_range_Wright_full, avg_spec_vol_Wright_full use MOM_EOS_Wright_red, only : calculate_density_wright_red, calculate_spec_vol_wright_red use MOM_EOS_Wright_red, only : calculate_density_derivs_wright_red use MOM_EOS_Wright_red, only : calculate_specvol_derivs_wright_red, int_density_dz_wright_red use MOM_EOS_Wright_red, only : calculate_compress_wright_red, int_spec_vol_dp_wright_red use MOM_EOS_Wright_red, only : calculate_density_second_derivs_wright_red -use MOM_EOS_Wright_red, only : EoS_fit_range_Wright_red +use MOM_EOS_Wright_red, only : EoS_fit_range_Wright_red, avg_spec_vol_Wright_red use MOM_EOS_Jackett06, only : calculate_density_Jackett06, calculate_spec_vol_Jackett06 use MOM_EOS_Jackett06, only : calculate_density_derivs_Jackett06, calculate_specvol_derivs_Jackett06 use MOM_EOS_Jackett06, only : calculate_compress_Jackett06, calculate_density_second_derivs_Jackett06 @@ -68,6 +69,7 @@ module MOM_EOS public EOS_unit_tests public analytic_int_density_dz public analytic_int_specific_vol_dp +public average_specific_vol public calculate_compress public calculate_density public calculate_density_derivs @@ -1324,6 +1326,97 @@ subroutine calculate_compress_scalar(T, S, pressure, rho, drho_dp, EOS) end subroutine calculate_compress_scalar +!> Calls the appropriate subroutine to calculate the layer averaged specific volume either using +!! Boole's rule quadrature or analytical and nearly-analytical averages in pressure. +subroutine average_specific_vol(T, S, p_t, dp, SpV_avg, EOS, dom, scale) + real, dimension(:), intent(in) :: T !< Potential temperature referenced to the surface [C ~> degC] + real, dimension(:), intent(in) :: S !< Salinity [S ~> ppt] + real, dimension(:), intent(in) :: p_t !< Pressure at the top of the layer [R L2 T-2 ~> Pa] + real, dimension(:), intent(in) :: dp !< Pressure change in the layer [R L2 T-2 ~> Pa] + real, dimension(:), intent(inout) :: SpV_avg !< The vertical average specific volume + !! in the layer [R-1 ~> m3 kg-1] + type(EOS_type), intent(in) :: EOS !< Equation of state structure + integer, dimension(2), optional, intent(in) :: dom !< The domain of indices to work on, taking + !! into account that arrays start at 1. + real, optional, intent(in) :: scale !< A multiplicative factor by which to scale + !! output specific volume in combination with + !! scaling stored in EOS [various] + + ! Local variables + real, dimension(size(T)) :: pres ! Layer-top pressure converted to [Pa] + real, dimension(size(T)) :: dpres ! Pressure change converted to [Pa] + real, dimension(size(T)) :: Ta ! Temperature converted to [degC] + real, dimension(size(T)) :: Sa ! Salinity converted to [ppt] + real :: T5(5) ! Temperatures at five quadrature points [C ~> degC] + real :: S5(5) ! Salinities at five quadrature points [S ~> ppt] + real :: p5(5) ! Pressures at five quadrature points [R L2 T-2 ~> Pa] + real :: a5(5) ! Specific volumes at five quadrature points [R-1 ~> m3 kg-1] + real, parameter :: C1_90 = 1.0/90.0 ! A rational constant [nondim] + real :: spv_scale ! A factor to convert specific volume from m3 kg-1 to the desired units [kg m-3 R-1 ~> 1] + integer :: i, n, is, ie, npts + + if (present(dom)) then + is = dom(1) ; ie = dom(2) ; npts = 1 + ie - is + else + is = 1 ; ie = size(T) ; npts = 1 + ie - is + endif + + if (EOS%EOS_quadrature) then + do i=is,ie + do n=1,5 + T5(n) = T(i) ; S5(n) = S(i) + p5(n) = p_t(i) + 0.25*real(5-n)*dp(i) + enddo + call calculate_spec_vol(T5, S5, p5, a5, EOS) + + ! Use Boole's rule to estimate the average specific volume. + SpV_avg(i) = C1_90*(7.0*(a5(1)+a5(5)) + 32.0*(a5(2)+a5(4)) + 12.0*a5(3)) + enddo + elseif ((EOS%RL2_T2_to_Pa == 1.0) .and. (EOS%C_to_degC == 1.0) .and. (EOS%S_to_ppt == 1.0)) then + select case (EOS%form_of_EOS) + case (EOS_LINEAR) + call avg_spec_vol_linear(T, S, p_t, dp, SpV_avg, is, npts, EOS%Rho_T0_S0, & + EOS%dRho_dT, EOS%dRho_dS) + case (EOS_WRIGHT) + call avg_spec_vol_wright(T, S, p_t, dp, SpV_avg, is, npts) + case (EOS_WRIGHT_FULL) + call avg_spec_vol_wright_full(T, S, p_t, dp, SpV_avg, is, npts) + case (EOS_WRIGHT_REDUCED) + call avg_spec_vol_wright_red(T, S, p_t, dp, SpV_avg, is, npts) + case default + call MOM_error(FATAL, "No analytic average specific volume option is available with this EOS!") + end select + else + do i=is,ie + pres(i) = EOS%RL2_T2_to_Pa * p_t(i) + dpres(i) = EOS%RL2_T2_to_Pa * dp(i) + Ta(i) = EOS%C_to_degC * T(i) + Sa(i) = EOS%S_to_ppt * S(i) + enddo + select case (EOS%form_of_EOS) + case (EOS_LINEAR) + call avg_spec_vol_linear(Ta, Sa, pres, dpres, SpV_avg, is, npts, EOS%Rho_T0_S0, & + EOS%dRho_dT, EOS%dRho_dS) + case (EOS_WRIGHT) + call avg_spec_vol_wright(Ta, Sa, pres, dpres, SpV_avg, is, npts) + case (EOS_WRIGHT_FULL) + call avg_spec_vol_wright_full(Ta, Sa, pres, dpres, SpV_avg, is, npts) + case (EOS_WRIGHT_REDUCED) + call avg_spec_vol_wright_red(Ta, Sa, pres, dpres, SpV_avg, is, npts) + case default + call MOM_error(FATAL, "No analytic average specific volume option is available with this EOS!") + end select + endif + + spv_scale = EOS%R_to_kg_m3 + if (EOS%EOS_quadrature) spv_scale = 1.0 + if (present(scale)) spv_scale = spv_scale * scale + if (spv_scale /= 1.0) then ; do i=is,ie + SpV_avg(i) = spv_scale * SpV_avg(i) + enddo ; endif + +end subroutine average_specific_vol + !> Return the range of temperatures, salinities and pressures for which the equation of state that !! is being used has been fitted to observations. Care should be taken when applying !! this equation of state outside of its fit range. @@ -2057,13 +2150,13 @@ logical function EOS_unit_tests(verbose) call EOS_manual_init(EOS_tmp, form_of_EOS=EOS_WRIGHT_FULL) fail = test_EOS_consistency(25.0, 35.0, 1.0e7, EOS_tmp, verbose, "WRIGHT_FULL", & - rho_check=1027.55177447616*EOS_tmp%kg_m3_to_R) + rho_check=1027.55177447616*EOS_tmp%kg_m3_to_R, avg_Sv_check=.true.) if (verbose .and. fail) call MOM_error(WARNING, "WRIGHT_FULL EOS has failed some self-consistency tests.") EOS_unit_tests = EOS_unit_tests .or. fail call EOS_manual_init(EOS_tmp, form_of_EOS=EOS_WRIGHT_REDUCED) fail = test_EOS_consistency(25.0, 35.0, 1.0e7, EOS_tmp, verbose, "WRIGHT_REDUCED", & - rho_check=1027.54303596346*EOS_tmp%kg_m3_to_R) + rho_check=1027.54303596346*EOS_tmp%kg_m3_to_R, avg_Sv_check=.true.) if (verbose .and. fail) call MOM_error(WARNING, "WRIGHT_REDUCED EOS has failed some self-consistency tests.") EOS_unit_tests = EOS_unit_tests .or. fail @@ -2076,7 +2169,7 @@ logical function EOS_unit_tests(verbose) call EOS_manual_init(EOS_tmp, form_of_EOS=EOS_WRIGHT) fail = test_EOS_consistency(25.0, 35.0, 1.0e7, EOS_tmp, verbose, "WRIGHT", & - rho_check=1027.54303596346*EOS_tmp%kg_m3_to_R) + rho_check=1027.54303596346*EOS_tmp%kg_m3_to_R, avg_Sv_check=.true.) if (verbose .and. fail) call MOM_error(WARNING, "WRIGHT EOS has failed some self-consistency tests.") EOS_unit_tests = EOS_unit_tests .or. fail @@ -2126,7 +2219,7 @@ logical function EOS_unit_tests(verbose) call EOS_manual_init(EOS_tmp, form_of_EOS=EOS_LINEAR, Rho_T0_S0=1000.0, drho_dT=-0.2, dRho_dS=0.8) fail = test_EOS_consistency(25.0, 35.0, 1.0e7, EOS_tmp, verbose, "LINEAR", & - rho_check=1023.0*EOS_tmp%kg_m3_to_R) + rho_check=1023.0*EOS_tmp%kg_m3_to_R, avg_Sv_check=.true.) if (verbose .and. fail) call MOM_error(WARNING, "LINEAR EOS has failed some self-consistency tests.") EOS_unit_tests = EOS_unit_tests .or. fail @@ -2293,7 +2386,7 @@ end subroutine write_check_msg !> Test an equation of state for self-consistency and consistency with check values, returning false !! if it is consistent by all tests, and true if it fails any test. logical function test_EOS_consistency(T_test, S_test, p_test, EOS, verbose, & - EOS_name, rho_check, spv_check, skip_2nd) result(inconsistent) + EOS_name, rho_check, spv_check, skip_2nd, avg_Sv_check) result(inconsistent) real, intent(in) :: T_test !< Potential temperature or conservative temperature [C ~> degC] real, intent(in) :: S_test !< Salinity or absolute salinity [S ~> ppt] real, intent(in) :: p_test !< Pressure [R L2 T-2 ~> Pa] @@ -2302,7 +2395,9 @@ logical function test_EOS_consistency(T_test, S_test, p_test, EOS, verbose, & character(len=*), intent(in) :: EOS_name !< A name used in error messages to describe the EoS real, optional, intent(in) :: rho_check !< A check value for the density [R ~> kg m-3] real, optional, intent(in) :: spv_check !< A check value for the specific volume [R-1 ~> m3 kg-1] - logical, optional, intent(in) :: skip_2nd !< If present and true, do not check the 2nd derivatives. + logical, optional, intent(in) :: skip_2nd !< If present and true, do not check the 2nd derivatives. + logical, optional, intent(in) :: avg_Sv_check !< If present and true, compare analytical and numerical + !! quadrature estimates of the layer-averaged specific volume. ! Local variables real, dimension(-3:3,-3:3,-3:3) :: T ! Temperatures at the test value and perturbed points [C ~> degC] @@ -2329,6 +2424,8 @@ logical function test_EOS_consistency(T_test, S_test, p_test, EOS, verbose, & ! temperature [R-1 C-1 ~> m3 kg-1 degC-1] real :: dSV_dS(1) ! The partial derivative of specific volume with salinity ! [R-1 S-1 ~> m3 kg-1 ppt-1] + real :: SpV_avg_a(1) ! The pressure-averaged specific volume determined analytically [R-1 ~> m3 kg-1] + real :: SpV_avg_q(1) ! The pressure-averaged specific volume determined via quadrature [R-1 ~> m3 kg-1] real :: drho_dS_dS ! Second derivative of density with respect to S [R S-2 ~> kg m-3 ppt-2] real :: drho_dS_dT ! Second derivative of density with respect to T and S [R S-1 C-1 ~> kg m-3 ppt-1 degC-1] real :: drho_dT_dT ! Second derivative of density with respect to T [R C-2 ~> kg m-3 degC-2] @@ -2370,13 +2467,17 @@ logical function test_EOS_consistency(T_test, S_test, p_test, EOS, verbose, & real :: count_fac2 ! A factor in the roundoff estimates based on the factors in the numerator and ! denominator in the finite difference second derivative expression [nondim] character(len=200) :: mesg + type(EOS_type) :: EOS_tmp logical :: test_OK ! True if a particular test is consistent. logical :: OK ! True if all checks so far are consistent. logical :: test_2nd ! If true, do tests on the 2nd derivative calculations + logical :: test_avg_Sv ! If true, compare numerical and analytical estimates of the vertically + ! averaged specific volume integer :: order ! The order of accuracy of the centered finite difference estimates (2, 4 or 6). integer :: i, j, k, n test_2nd = .true. ; if (present(skip_2nd)) test_2nd = .not.skip_2nd + test_avg_Sv = .false. ; if (present(avg_Sv_check)) test_avg_Sv = avg_Sv_check dT = 0.1*EOS%degC_to_C ! Temperature perturbations [C ~> degC] dS = 0.5*EOS%ppt_to_S ! Salinity perturbations [S ~> ppt] @@ -2442,6 +2543,14 @@ logical function test_EOS_consistency(T_test, S_test, p_test, EOS, verbose, & drho_dS_dS, drho_dS_dT, drho_dT_dT, drho_dS_dP, drho_dT_dP, EOS) call calculate_compress(T(0,0,0), S(0,0,0), p(0,0,0), rho_tmp, drho_dp, EOS) + if (test_avg_Sv) then + EOS_tmp = EOS + call EOS_manual_init(EOS_tmp, EOS_quadrature=.false.) + call average_specific_vol(T(0:0,0,0), S(0:0,0,0), p(0:0,0,0), p(0:0,0,0), SpV_avg_a, EOS_tmp) + call EOS_manual_init(EOS_tmp, EOS_quadrature=.true.) + call average_specific_vol(T(0:0,0,0), S(0:0,0,0), p(0:0,0,0), p(0:0,0,0), SpV_avg_q, EOS_tmp) + endif + OK = .true. tol = 1000.0*epsilon(tol) @@ -2532,6 +2641,23 @@ logical function test_EOS_consistency(T_test, S_test, p_test, EOS, verbose, & OK = OK .and. check_FD(drho_dS_dP, drho_dS_dP_fd, tol_here, verbose, trim(EOS_name)//" drho_dS_dP", order) endif + if (test_avg_Sv) then + tol_here = 0.5*tol*(abs(SpV_avg_a(1)) + abs(SpV_avg_q(1))) + test_OK = (abs(SpV_avg_a(1) - SpV_avg_q(1)) < tol_here) + if (verbose) then + write(mesg, '(ES24.16," and ",ES24.16," differ by ",ES16.8," (",ES10.2"), tol=",ES16.8)') & + SpV_avg_a(1), SpV_avg_q(1), SpV_avg_a(1) - SpV_avg_q(1), & + 2.0*(SpV_avg_a(1) - SpV_avg_q(1)) / (abs(SpV_avg_a(1)) + abs(SpV_avg_q(1)) + tiny(SpV_avg_a(1))), & + tol_here + if (verbose .and. .not.test_OK) then + call MOM_error(WARNING, "The values of "//trim(EOS_name)//" SpV_avg disagree. "//trim(mesg)) + elseif (verbose) then + call MOM_mesg("The values of "//trim(EOS_name)//" SpV_avg agree: "//trim(mesg)) + endif + endif + OK = OK .and. test_OK + endif + inconsistent = .not.OK contains diff --git a/src/equation_of_state/MOM_EOS_Wright.F90 b/src/equation_of_state/MOM_EOS_Wright.F90 index 25ae9219a8..d8dee28aa2 100644 --- a/src/equation_of_state/MOM_EOS_Wright.F90 +++ b/src/equation_of_state/MOM_EOS_Wright.F90 @@ -10,7 +10,7 @@ module MOM_EOS_Wright public calculate_compress_wright, calculate_density_wright, calculate_spec_vol_wright public calculate_density_derivs_wright, calculate_specvol_derivs_wright public calculate_density_second_derivs_wright, calc_density_second_derivs_wright_buggy -public EoS_fit_range_Wright +public EoS_fit_range_Wright, avg_spec_vol_Wright public int_density_dz_wright, int_spec_vol_dp_wright !> Compute the in situ density of sea water (in [kg m-3]), or its anomaly with respect to @@ -547,6 +547,42 @@ subroutine calculate_compress_wright(T, S, pressure, rho, drho_dp, start, npts) enddo end subroutine calculate_compress_wright +!> Calculates analytical and nearly-analytical integrals, in pressure across layers, to determine +!! the layer-average specific volumes. There are essentially no free assumptions, apart from a +!! truncation in the series for log(1-eps/1+eps) that assumes that |eps| < 0.34. +subroutine avg_spec_vol_Wright(T, S, p_t, dp, SpV_avg, start, npts) + real, dimension(:), intent(in) :: T !< Potential temperature relative to the surface + !! [degC]. + real, dimension(:), intent(in) :: S !< Salinity [PSU]. + real, dimension(:), intent(in) :: p_t !< Pressure at the top of the layer [Pa] + real, dimension(:), intent(in) :: dp !< Pressure change in the layer [Pa] + real, dimension(:), intent(inout) :: SpV_avg !< The vertical average specific volume + !! in the layer [m3 kg-1] + integer, intent(in) :: start !< the starting point in the arrays. + integer, intent(in) :: npts !< the number of values to calculate. + + ! Local variables + real :: al0 ! A term in the Wright EOS [m3 kg-1] + real :: p0 ! A term in the Wright EOS [Pa] + real :: lambda ! A term in the Wright EOS [m2 s-2] + real :: eps2 ! The square of a nondimensional ratio [nondim] + real :: I_pterm ! The inverse of p0 plus p_ave [Pa-1]. + real, parameter :: C1_3 = 1.0/3.0, C1_7 = 1.0/7.0, C1_9 = 1.0/9.0 ! Rational constants [nondim] + integer :: j + + ! alpha(j) = al0 + lambda / (pressure(j) + p0) + do j=start,start+npts-1 + al0 = a0 + (a1*T(j) + a2*S(j)) + p0 = b0 + ( b4*S(j) + T(j) * (b1 + (T(j)*(b2 + b3*T(j)) + b5*S(j))) ) + lambda = c0 + ( c4*S(j) + T(j) * (c1 + (T(j)*(c2 + c3*T(j)) + c5*S(j))) ) + + I_pterm = 1.0 / (p0 + (p_t(j) + 0.5*dp(j))) + eps2 = (0.5 * dp(j) * I_pterm)**2 + SpV_avg(j) = al0 + (lambda * I_pterm) * & + (1.0 + eps2*(C1_3 + eps2*(0.2 + eps2*(C1_7 + eps2*C1_9)))) + enddo +end subroutine avg_spec_vol_Wright + !> Return the range of temperatures, salinities and pressures for which the reduced-range equation !! of state from Wright (1997) has been fitted to observations. Care should be taken when applying !! this equation of state outside of its fit range. @@ -1066,6 +1102,7 @@ subroutine int_spec_vol_dp_wright(T, S, p_t, p_b, spv_ref, HI, dza, & enddo ; enddo ; endif end subroutine int_spec_vol_dp_wright + !> \namespace mom_eos_wright !! !! \section section_EOS_Wright Wright equation of state diff --git a/src/equation_of_state/MOM_EOS_Wright_full.F90 b/src/equation_of_state/MOM_EOS_Wright_full.F90 index 3f00a92cef..107ced3f5b 100644 --- a/src/equation_of_state/MOM_EOS_Wright_full.F90 +++ b/src/equation_of_state/MOM_EOS_Wright_full.F90 @@ -11,6 +11,7 @@ module MOM_EOS_Wright_full public calculate_density_derivs_wright_full, calculate_specvol_derivs_wright_full public calculate_density_second_derivs_wright_full, EoS_fit_range_Wright_full public int_density_dz_wright_full, int_spec_vol_dp_wright_full +public avg_spec_vol_Wright_full !> Compute the in situ density of sea water (in [kg m-3]), or its anomaly with respect to !! a reference density, from salinity in practical salinity units ([PSU]), potential @@ -450,6 +451,42 @@ subroutine calculate_compress_wright_full(T, S, pressure, rho, drho_dp, start, n enddo end subroutine calculate_compress_wright_full +!> Calculates analytical and nearly-analytical integrals, in pressure across layers, to determine +!! the layer-average specific volumes. There are essentially no free assumptions, apart from a +!! truncation in the series for log(1-eps/1+eps) that assumes that |eps| < 0.34. +subroutine avg_spec_vol_Wright_full(T, S, p_t, dp, SpV_avg, start, npts) + real, dimension(:), intent(in) :: T !< Potential temperature relative to the surface + !! [degC]. + real, dimension(:), intent(in) :: S !< Salinity [PSU]. + real, dimension(:), intent(in) :: p_t !< Pressure at the top of the layer [Pa] + real, dimension(:), intent(in) :: dp !< Pressure change in the layer [Pa] + real, dimension(:), intent(inout) :: SpV_avg !< The vertical average specific volume + !! in the layer [m3 kg-1] + integer, intent(in) :: start !< the starting point in the arrays. + integer, intent(in) :: npts !< the number of values to calculate. + + ! Local variables + real :: al0 ! A term in the Wright EOS [m3 kg-1] + real :: p0 ! A term in the Wright EOS [Pa] + real :: lambda ! A term in the Wright EOS [m2 s-2] + real :: eps2 ! The square of a nondimensional ratio [nondim] + real :: I_pterm ! The inverse of p0 plus p_ave [Pa-1]. + real, parameter :: C1_3 = 1.0/3.0, C1_7 = 1.0/7.0, C1_9 = 1.0/9.0 ! Rational constants [nondim] + integer :: j + + ! alpha(j) = al0 + lambda / (pressure(j) + p0) + do j=start,start+npts-1 + al0 = a0 + (a1*T(j) + a2*S(j)) + p0 = b0 + ( b4*S(j) + T(j) * (b1 + (T(j)*(b2 + b3*T(j)) + b5*S(j))) ) + lambda = c0 + ( c4*S(j) + T(j) * (c1 + (T(j)*(c2 + c3*T(j)) + c5*S(j))) ) + + I_pterm = 1.0 / (p0 + (p_t(j) + 0.5*dp(j))) + eps2 = (0.5 * dp(j) * I_pterm)**2 + SpV_avg(j) = al0 + (lambda * I_pterm) * & + (1.0 + eps2*(C1_3 + eps2*(0.2 + eps2*(C1_7 + eps2*C1_9)))) + enddo +end subroutine avg_spec_vol_Wright_full + !> Return the range of temperatures, salinities and pressures for which full-range equation !! of state from Wright (1997) has been fitted to observations. Care should be taken when applying !! this equation of state outside of its fit range. @@ -972,6 +1009,7 @@ subroutine int_spec_vol_dp_wright_full(T, S, p_t, p_b, spv_ref, HI, dza, & enddo ; enddo ; endif end subroutine int_spec_vol_dp_wright_full + !> \namespace mom_eos_wright_full !! !! \section section_EOS_Wright_full Wright equation of state diff --git a/src/equation_of_state/MOM_EOS_Wright_red.F90 b/src/equation_of_state/MOM_EOS_Wright_red.F90 index cf78ce2211..5553112274 100644 --- a/src/equation_of_state/MOM_EOS_Wright_red.F90 +++ b/src/equation_of_state/MOM_EOS_Wright_red.F90 @@ -11,6 +11,7 @@ module MOM_EOS_Wright_red public calculate_density_derivs_wright_red, calculate_specvol_derivs_wright_red public calculate_density_second_derivs_wright_red, EoS_fit_range_Wright_red public int_density_dz_wright_red, int_spec_vol_dp_wright_red +public avg_spec_vol_Wright_red !> Compute the in situ density of sea water (in [kg m-3]), or its anomaly with respect to !! a reference density, from salinity in practical salinity units ([PSU]), potential @@ -450,6 +451,42 @@ subroutine calculate_compress_wright_red(T, S, pressure, rho, drho_dp, start, np enddo end subroutine calculate_compress_wright_red +!> Calculates analytical and nearly-analytical integrals, in pressure across layers, to determine +!! the layer-average specific volumes. There are essentially no free assumptions, apart from a +!! truncation in the series for log(1-eps/1+eps) that assumes that |eps| < 0.34. +subroutine avg_spec_vol_Wright_red(T, S, p_t, dp, SpV_avg, start, npts) + real, dimension(:), intent(in) :: T !< Potential temperature relative to the surface + !! [degC]. + real, dimension(:), intent(in) :: S !< Salinity [PSU]. + real, dimension(:), intent(in) :: p_t !< Pressure at the top of the layer [Pa] + real, dimension(:), intent(in) :: dp !< Pressure change in the layer [Pa] + real, dimension(:), intent(inout) :: SpV_avg !< The vertical average specific volume + !! in the layer [m3 kg-1] + integer, intent(in) :: start !< the starting point in the arrays. + integer, intent(in) :: npts !< the number of values to calculate. + + ! Local variables + real :: al0 ! A term in the Wright EOS [m3 kg-1] + real :: p0 ! A term in the Wright EOS [Pa] + real :: lambda ! A term in the Wright EOS [m2 s-2] + real :: eps2 ! The square of a nondimensional ratio [nondim] + real :: I_pterm ! The inverse of p0 plus p_ave [Pa-1]. + real, parameter :: C1_3 = 1.0/3.0, C1_7 = 1.0/7.0, C1_9 = 1.0/9.0 ! Rational constants [nondim] + integer :: j + + ! alpha(j) = al0 + lambda / (pressure(j) + p0) + do j=start,start+npts-1 + al0 = a0 + (a1*T(j) + a2*S(j)) + p0 = b0 + ( b4*S(j) + T(j) * (b1 + (T(j)*(b2 + b3*T(j)) + b5*S(j))) ) + lambda = c0 + ( c4*S(j) + T(j) * (c1 + (T(j)*(c2 + c3*T(j)) + c5*S(j))) ) + + I_pterm = 1.0 / (p0 + (p_t(j) + 0.5*dp(j))) + eps2 = (0.5 * dp(j) * I_pterm)**2 + SpV_avg(j) = al0 + (lambda * I_pterm) * & + (1.0 + eps2*(C1_3 + eps2*(0.2 + eps2*(C1_7 + eps2*C1_9)))) + enddo +end subroutine avg_spec_vol_Wright_red + !> Return the range of temperatures, salinities and pressures for which the reduced-range equation !! of state from Wright (1997) has been fitted to observations. Care should be taken when applying !! this equation of state outside of its fit range. @@ -972,6 +1009,7 @@ subroutine int_spec_vol_dp_wright_red(T, S, p_t, p_b, spv_ref, HI, dza, & enddo ; enddo ; endif end subroutine int_spec_vol_dp_wright_red + !> \namespace mom_eos_wright_red !! !! \section section_EOS_Wright_red Wright equation of state diff --git a/src/equation_of_state/MOM_EOS_linear.F90 b/src/equation_of_state/MOM_EOS_linear.F90 index 1899103f5d..b1dacf2780 100644 --- a/src/equation_of_state/MOM_EOS_linear.F90 +++ b/src/equation_of_state/MOM_EOS_linear.F90 @@ -13,6 +13,7 @@ module MOM_EOS_linear public calculate_density_scalar_linear, calculate_density_array_linear public calculate_density_second_derivs_linear, EoS_fit_range_linear public int_density_dz_linear, int_spec_vol_dp_linear +public avg_spec_vol_linear ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with @@ -292,7 +293,7 @@ end subroutine calculate_specvol_derivs_linear !> This subroutine computes the in situ density of sea water (rho) !! and the compressibility (drho/dp == C_sound^-2) at the given !! salinity, potential temperature, and pressure. -subroutine calculate_compress_linear(T, S, pressure, rho, drho_dp, start, npts,& +subroutine calculate_compress_linear(T, S, pressure, rho, drho_dp, start, npts, & Rho_T0_S0, dRho_dT, dRho_dS) real, intent(in), dimension(:) :: T !< Potential temperature relative to the surface !! [degC]. @@ -318,6 +319,29 @@ subroutine calculate_compress_linear(T, S, pressure, rho, drho_dp, start, npts,& enddo end subroutine calculate_compress_linear +!> Calculates the layer average specific volumes. +subroutine avg_spec_vol_linear(T, S, p_t, dp, SpV_avg, start, npts, Rho_T0_S0, dRho_dT, dRho_dS) + real, dimension(:), intent(in) :: T !< Potential temperature [degC] + real, dimension(:), intent(in) :: S !< Salinity [PSU] + real, dimension(:), intent(in) :: p_t !< Pressure at the top of the layer [Pa] + real, dimension(:), intent(in) :: dp !< Pressure change in the layer [Pa] + real, dimension(:), intent(inout) :: SpV_avg !< The vertical average specific volume + !! in the layer [m3 kg-1] + integer, intent(in) :: start !< the starting point in the arrays. + integer, intent(in) :: npts !< the number of values to calculate. + real, intent(in) :: Rho_T0_S0 !< The density at T=0, S=0 [kg m-3] + real, intent(in) :: dRho_dT !< The derivative of density with temperature + !! [kg m-3 degC-1] + real, intent(in) :: dRho_dS !< The derivative of density with salinity + !! [kg m-3 ppt-1] + ! Local variables + integer :: j + + do j=start,start+npts-1 + SpV_avg(j) = 1.0 / (Rho_T0_S0 + (dRho_dT*T(j) + dRho_dS*S(j))) + enddo +end subroutine avg_spec_vol_linear + !> Return the range of temperatures, salinities and pressures for which the reduced-range equation !! of state from Wright (1997) has been fitted to observations. Care should be taken when applying !! this equation of state outside of its fit range. diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 59d5aaf60a..0fe08a06b2 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -43,7 +43,7 @@ module MOM_diabatic_driver use MOM_grid, only : ocean_grid_type use MOM_int_tide_input, only : set_int_tide_input, int_tide_input_init use MOM_int_tide_input, only : int_tide_input_end, int_tide_input_CS, int_tide_input_type -use MOM_interface_heights, only : find_eta +use MOM_interface_heights, only : find_eta, calc_derived_thermo use MOM_internal_tides, only : propagate_int_tide use MOM_internal_tides, only : internal_tides_init, internal_tides_end, int_tide_CS use MOM_kappa_shear, only : kappa_shear_is_used @@ -1844,9 +1844,15 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e ! Also changes: visc%Kd_shear and visc%Kv_shear if ((CS%halo_TS_diff > 0) .and. (CS%ML_mix_first > 0.0)) then if (associated(tv%T)) call pass_var(tv%T, G%Domain, halo=CS%halo_TS_diff, complete=.false.) - if (associated(tv%T)) call pass_var(tv%S, G%Domain, halo=CS%halo_TS_diff, complete=.false.) + if (associated(tv%S)) call pass_var(tv%S, G%Domain, halo=CS%halo_TS_diff, complete=.false.) call pass_var(h, G%domain, halo=CS%halo_TS_diff, complete=.true.) endif + + ! Update derived thermodynamic quantities. + if ((CS%ML_mix_first > 0.0) .and. allocated(tv%SpV_avg)) then + call calc_derived_thermo(tv, h, G, GV, US, halo=CS%halo_TS_diff) + endif + if (CS%debug) & call MOM_state_chksum("before set_diffusivity", u, v, h, G, GV, US, haloshift=CS%halo_TS_diff) if (CS%double_diffuse) then diff --git a/src/tracer/MOM_offline_main.F90 b/src/tracer/MOM_offline_main.F90 index ea6167a6b8..40dced9b20 100644 --- a/src/tracer/MOM_offline_main.F90 +++ b/src/tracer/MOM_offline_main.F90 @@ -22,6 +22,7 @@ module MOM_offline_main use MOM_file_parser, only : read_param, get_param, log_version, param_file_type use MOM_forcing_type, only : forcing use MOM_grid, only : ocean_grid_type +use MOM_interface_heights, only : calc_derived_thermo use MOM_io, only : MOM_read_data, MOM_read_vector, CENTER use MOM_offline_aux, only : update_offline_from_arrays, update_offline_from_files use MOM_offline_aux, only : next_modulo_time, offline_add_diurnal_sw @@ -1025,6 +1026,7 @@ subroutine update_offline_fields(CS, G, GV, US, h, fluxes, do_ale) type(forcing), intent(inout) :: fluxes !< Pointers to forcing fields logical, intent(in ) :: do_ale !< True if using ALE ! Local variables + integer :: stencil integer :: i, j, k, is, ie, js, je, nz real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_start ! Initial thicknesses [H ~> m or kg m-2] is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke @@ -1086,6 +1088,12 @@ subroutine update_offline_fields(CS, G, GV, US, h, fluxes, do_ale) call pass_var(CS%tv%T, G%Domain) call pass_var(CS%tv%S, G%Domain) + ! Update derived thermodynamic quantities. + if (allocated(CS%tv%SpV_avg)) then + stencil = min(3, G%Domain%nihalo, G%Domain%njhalo) + call calc_derived_thermo(CS%tv, CS%h_end, G, GV, US, halo=stencil) + endif + ! Update the read indices CS%ridx_snap = next_modulo_time(CS%ridx_snap,CS%numtime) CS%ridx_sum = next_modulo_time(CS%ridx_sum,CS%numtime)