From ce0905e17e7fc4b0e1f4f55f6ece1d75234d8826 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 10 Apr 2020 06:03:18 -0400 Subject: [PATCH] +Add hor_index_type variants of EOS routines Added variants of calculate_density routines that use a hor_index_type to specify array extents and unit_scale_types for dimensional consistency testing, further overloading existing interfaces. Also replaced the recently added rho_scale and pres_scale arguments to int_density_dz with an optional unit_scale_type argument, and modified calls to use this new argument. All answers are bitwise identical, but there are changes to external interfaces. --- src/core/MOM_PressureForce_Montgomery.F90 | 27 +- src/core/MOM_PressureForce_analytic_FV.F90 | 39 +- src/core/MOM_PressureForce_blocked_AFV.F90 | 14 +- src/core/MOM_interface_heights.F90 | 6 +- src/diagnostics/MOM_diagnostics.F90 | 3 +- src/equation_of_state/MOM_EOS.F90 | 722 +++++++++++++++++---- 6 files changed, 627 insertions(+), 184 deletions(-) diff --git a/src/core/MOM_PressureForce_Montgomery.F90 b/src/core/MOM_PressureForce_Montgomery.F90 index 58687b874f..1665e28def 100644 --- a/src/core/MOM_PressureForce_Montgomery.F90 +++ b/src/core/MOM_PressureForce_Montgomery.F90 @@ -186,8 +186,7 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pb !$OMP parallel do default(shared) do k=1,nz call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p(:,:,k), p(:,:,k+1), & - 0.0, G%HI, tv%eqn_of_state, dz_geo(:,:,k), halo_size=1, & - SV_scale=US%R_to_kg_m3, pres_scale=US%RL2_T2_to_Pa) + 0.0, G%HI, tv%eqn_of_state, dz_geo(:,:,k), halo_size=1, US=US) enddo !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do k=1,nz ; do i=Isq,Ieq+1 @@ -659,8 +658,8 @@ subroutine Set_pbce_Bouss(e, tv, G, GV, US, Rho0, GFS_scale, pbce, rho_star) Ihtot(i) = GV%H_to_Z / ((e(i,j,1)-e(i,j,nz+1)) + z_neglect) press(i) = -Rho0xG*e(i,j,1) enddo - call calculate_density(tv%T(:,j,1), tv%S(:,j,1), press, rho_in_situ, Isq, Ieq-Isq+2, & - tv%eqn_of_state, scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa) + call calculate_density(tv%T(:,j,1), tv%S(:,j,1), press, rho_in_situ, G%HI, & + tv%eqn_of_state, US, halo=1) do i=Isq,Ieq+1 pbce(i,j,1) = G_Rho0*(GFS_scale * rho_in_situ(i)) * GV%H_to_Z enddo @@ -670,8 +669,8 @@ subroutine Set_pbce_Bouss(e, tv, G, GV, US, Rho0, GFS_scale, pbce, rho_star) T_int(i) = 0.5*(tv%T(i,j,k-1)+tv%T(i,j,k)) S_int(i) = 0.5*(tv%S(i,j,k-1)+tv%S(i,j,k)) enddo - call calculate_density_derivs(T_int, S_int, press, dR_dT, dR_dS, Isq, Ieq-Isq+2, & - tv%eqn_of_state, scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa) + call calculate_density_derivs(T_int, S_int, press, dR_dT, dR_dS, G%HI, & + tv%eqn_of_state, US, halo=1) do i=Isq,Ieq+1 pbce(i,j,k) = pbce(i,j,k-1) + G_Rho0 * & ((e(i,j,K) - e(i,j,nz+1)) * Ihtot(i)) * & @@ -755,9 +754,8 @@ subroutine Set_pbce_nonBouss(p, tv, G, GV, US, GFS_scale, pbce, alpha_star) else !$OMP parallel do default(shared) private(T_int,S_int,dR_dT,dR_dS,rho_in_situ) do j=Jsq,Jeq+1 - call calculate_density(tv%T(:,j,nz), tv%S(:,j,nz), p(:,j,nz+1), rho_in_situ, & - Isq, Ieq-Isq+2, tv%eqn_of_state, scale=US%kg_m3_to_R, & - pres_scale=US%RL2_T2_to_Pa) + call calculate_density(tv%T(:,j,nz), tv%S(:,j,nz), p(:,j,nz+1), rho_in_situ, G%HI, & + tv%eqn_of_state, US, halo=1) do i=Isq,Ieq+1 C_htot(i,j) = dP_dH / ((p(i,j,nz+1)-p(i,j,1)) + dp_neglect) pbce(i,j,nz) = dP_dH / (rho_in_situ(i)) @@ -767,11 +765,12 @@ subroutine Set_pbce_nonBouss(p, tv, G, GV, US, GFS_scale, pbce, alpha_star) T_int(i) = 0.5*(tv%T(i,j,k)+tv%T(i,j,k+1)) S_int(i) = 0.5*(tv%S(i,j,k)+tv%S(i,j,k+1)) enddo - call calculate_density(T_int, S_int, p(:,j,k+1), rho_in_situ, Isq, Ieq-Isq+2, tv%eqn_of_state, & - scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa) - call calculate_density_derivs(T_int, S_int, p(:,j,k+1), dR_dT, dR_dS, & - Isq, Ieq-Isq+2, tv%eqn_of_state, & - scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa) + call calculate_density(T_int, S_int, p(:,j,k+1), rho_in_situ, G%HI, tv%eqn_of_state, US, halo=1) +! call calculate_density_derivs(T_int, S_int, p(:,j,k+1), dR_dT, dR_dS, & +! Isq, Ieq-Isq+2, tv%eqn_of_state, & +! scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa) + call calculate_density_derivs(T_int, S_int, p(:,j,k+1), dR_dT, dR_dS, G%HI, & + tv%eqn_of_state, US, halo=1) do i=Isq,Ieq+1 pbce(i,j,k) = pbce(i,j,k+1) + ((p(i,j,K+1)-p(i,j,1))*C_htot(i,j)) * & ((dR_dT(i)*(tv%T(i,j,k+1)-tv%T(i,j,k)) + & diff --git a/src/core/MOM_PressureForce_analytic_FV.F90 b/src/core/MOM_PressureForce_analytic_FV.F90 index aaab3d822f..10842f9e7e 100644 --- a/src/core/MOM_PressureForce_analytic_FV.F90 +++ b/src/core/MOM_PressureForce_analytic_FV.F90 @@ -264,24 +264,22 @@ subroutine PressureForce_AFV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p if ( CS%Recon_Scheme == 1 ) then call int_spec_vol_dp_generic_plm( T_t(:,:,k), T_b(:,:,k), S_t(:,:,k), S_b(:,:,k), & p(:,:,K), p(:,:,K+1), alpha_ref, dp_neglect, p(:,:,nz+1), G%HI, & - tv%eqn_of_state, dza(:,:,k), intp_dza(:,:,k), & - intx_dza(:,:,k), inty_dza(:,:,k), useMassWghtInterp=CS%useMassWghtInterp, & - SV_scale=US%R_to_kg_m3, pres_scale=US%RL2_T2_to_Pa) + tv%eqn_of_state, dza(:,:,k), intp_dza(:,:,k), intx_dza(:,:,k), inty_dza(:,:,k), & + useMassWghtInterp=CS%useMassWghtInterp, US=US) elseif ( CS%Recon_Scheme == 2 ) then call MOM_error(FATAL, "PressureForce_AFV_nonBouss: "//& "int_spec_vol_dp_generic_ppm does not exist yet.") ! call int_spec_vol_dp_generic_ppm ( tv%T(:,:,k), T_t(:,:,k), T_b(:,:,k), & ! tv%S(:,:,k), S_t(:,:,k), S_b(:,:,k), p(:,:,K), p(:,:,K+1), & ! alpha_ref, G%HI, tv%eqn_of_state, dza(:,:,k), intp_dza(:,:,k), & - ! intx_dza(:,:,k), inty_dza(:,:,k)) + ! intx_dza(:,:,k), inty_dza(:,:,k), US=US) endif else call int_specific_vol_dp(tv_tmp%T(:,:,k), tv_tmp%S(:,:,k), p(:,:,K), & p(:,:,K+1), alpha_ref, G%HI, tv%eqn_of_state, & dza(:,:,k), intp_dza(:,:,k), intx_dza(:,:,k), & inty_dza(:,:,k), bathyP=p(:,:,nz+1), dP_tiny=dp_neglect, & - useMassWghtInterp=CS%useMassWghtInterp, & - SV_scale=US%R_to_kg_m3, pres_scale=US%RL2_T2_to_Pa) + useMassWghtInterp=CS%useMassWghtInterp, US=US) endif else alpha_anom = 1.0 / GV%Rlay(k) - alpha_ref @@ -335,8 +333,8 @@ subroutine PressureForce_AFV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p if (use_EOS) then !$OMP parallel do default(shared) private(rho_in_situ) do j=Jsq,Jeq+1 - call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p(:,j,1), rho_in_situ, Isq, Ieq-Isq+2, & - tv%eqn_of_state, scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa) + call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p(:,j,1), rho_in_situ, G%HI, & + tv%eqn_of_state, US, halo=1) do i=Isq,Ieq+1 dM(i,j) = (CS%GFS_scale - 1.0) * (p(i,j,1)*(1.0/rho_in_situ(i) - alpha_ref) + za(i,j)) @@ -601,11 +599,11 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_at !$OMP parallel do default(shared) do j=Jsq,Jeq+1 if (use_p_atm) then - call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p_atm(:,j), rho_in_situ, & - Isq, Ieq-Isq+2, tv%eqn_of_state, scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa) + call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p_atm(:,j), rho_in_situ, G%HI, & + tv%eqn_of_state, US, halo=1) else - call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p0, rho_in_situ, & - Isq, Ieq-Isq+2, tv%eqn_of_state, scale=US%kg_m3_to_R) + call calculate_density(tv_tmp%T(:,j,1), tv_tmp%S(:,j,1), p0, rho_in_situ, G%HI, & + tv%eqn_of_state, US, halo=1) endif do i=Isq,Ieq+1 dM(i,j) = (CS%GFS_scale - 1.0) * (G_Rho0 * rho_in_situ(i)) * e(i,j,1) @@ -668,20 +666,19 @@ subroutine PressureForce_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_at if ( use_ALE ) then if ( CS%Recon_Scheme == 1 ) then call int_density_dz_generic_plm( T_t(:,:,k), T_b(:,:,k), S_t(:,:,k), S_b(:,:,k),& - e(:,:,K), e(:,:,K+1), rho_ref, CS%Rho0, GV%g_Earth, dz_neglect, G%bathyT, G%HI, G%HI, & - tv%eqn_of_state, dpa, intz_dpa, intx_dpa, inty_dpa, useMassWghtInterp=CS%useMassWghtInterp, & - rho_scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa) + e(:,:,K), e(:,:,K+1), rho_ref, CS%Rho0, GV%g_Earth, dz_neglect, G%bathyT, & + G%HI, G%HI, tv%eqn_of_state, dpa, intz_dpa, intx_dpa, inty_dpa, & + useMassWghtInterp=CS%useMassWghtInterp, US=US) elseif ( CS%Recon_Scheme == 2 ) then call int_density_dz_generic_ppm( tv%T(:,:,k), T_t(:,:,k), T_b(:,:,k), & - tv%S(:,:,k), S_t(:,:,k), S_b(:,:,k), e(:,:,K), e(:,:,K+1), rho_ref, CS%Rho0, & - GV%g_Earth, G%HI, G%HI, tv%eqn_of_state, dpa, intz_dpa, intx_dpa, inty_dpa, & - rho_scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa) + tv%S(:,:,k), S_t(:,:,k), S_b(:,:,k), e(:,:,K), e(:,:,K+1), & + rho_ref, CS%Rho0, GV%g_Earth, G%HI, G%HI, tv%eqn_of_state, dpa, & + intz_dpa, intx_dpa, inty_dpa, US=US) endif else call int_density_dz(tv_tmp%T(:,:,k), tv_tmp%S(:,:,k), e(:,:,K), e(:,:,K+1), & - rho_ref, CS%Rho0, GV%g_Earth, G%HI, G%HI, tv%eqn_of_state, & - dpa, intz_dpa, intx_dpa, inty_dpa, G%bathyT, dz_neglect, CS%useMassWghtInterp, & - rho_scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa) + rho_ref, CS%Rho0, GV%g_Earth, G%HI, G%HI, tv%eqn_of_state, dpa, & + intz_dpa, intx_dpa, inty_dpa, G%bathyT, dz_neglect, CS%useMassWghtInterp, US=US) endif !$OMP parallel do default(shared) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 diff --git a/src/core/MOM_PressureForce_blocked_AFV.F90 b/src/core/MOM_PressureForce_blocked_AFV.F90 index ae50019987..c532e17001 100644 --- a/src/core/MOM_PressureForce_blocked_AFV.F90 +++ b/src/core/MOM_PressureForce_blocked_AFV.F90 @@ -248,8 +248,7 @@ subroutine PressureForce_blk_AFV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, p(:,:,K+1), alpha_ref, G%HI, tv%eqn_of_state, & dza(:,:,k), intp_dza(:,:,k), intx_dza(:,:,k), & inty_dza(:,:,k), bathyP=p(:,:,nz+1), dP_tiny=dp_neglect, & - useMassWghtInterp=CS%useMassWghtInterp, & - SV_scale=US%R_to_kg_m3, pres_scale=US%RL2_T2_to_Pa) + useMassWghtInterp=CS%useMassWghtInterp, US=US) else alpha_anom = 1.0 / GV%Rlay(k) - alpha_ref do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 @@ -669,21 +668,18 @@ subroutine PressureForce_blk_AFV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, call int_density_dz_generic_plm( T_t(:,:,k), T_b(:,:,k), S_t(:,:,k), S_b(:,:,k), & e(:,:,K), e(:,:,K+1), rho_ref, CS%Rho0, GV%g_Earth, dz_neglect, G%bathyT, G%HI, & G%Block(n), tv%eqn_of_state, dpa_bk, intz_dpa_bk, intx_dpa_bk, inty_dpa_bk, & - useMassWghtInterp=CS%useMassWghtInterp, & - rho_scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa) + useMassWghtInterp=CS%useMassWghtInterp, US=US) elseif ( CS%Recon_Scheme == 2 ) then call int_density_dz_generic_ppm( tv%T(:,:,k), T_t(:,:,k), T_b(:,:,k), & tv%S(:,:,k), S_t(:,:,k), S_b(:,:,k), e(:,:,K), e(:,:,K+1), rho_ref, CS%Rho0, & - GV%g_Earth, G%HI, G%Block(n), tv%eqn_of_state, dpa_bk, intz_dpa_bk, & - intx_dpa_bk, inty_dpa_bk, & - rho_scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa) + GV%g_Earth, G%HI, G%Block(n), tv%eqn_of_state, dpa_bk, intz_dpa_bk, & + intx_dpa_bk, inty_dpa_bk, US=US) endif else call int_density_dz(tv_tmp%T(:,:,k), tv_tmp%S(:,:,k), e(:,:,K), e(:,:,K+1), & rho_ref, CS%Rho0, GV%g_Earth, G%HI, G%Block(n), tv%eqn_of_state, & dpa_bk, intz_dpa_bk, intx_dpa_bk, inty_dpa_bk, & - G%bathyT, dz_neglect, CS%useMassWghtInterp, & - rho_scale=US%kg_m3_to_R, pres_scale=US%RL2_T2_to_Pa) + G%bathyT, dz_neglect, CS%useMassWghtInterp, US=US) endif do jb=Jsq_bk,Jeq_bk+1 ; do ib=Isq_bk,Ieq_bk+1 intz_dpa_bk(ib,jb) = intz_dpa_bk(ib,jb)*GV%Z_to_H diff --git a/src/core/MOM_interface_heights.F90 b/src/core/MOM_interface_heights.F90 index c7147669dd..bfb9ad2703 100644 --- a/src/core/MOM_interface_heights.F90 +++ b/src/core/MOM_interface_heights.F90 @@ -106,8 +106,7 @@ subroutine find_eta_3d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m) !$OMP do do k=1,nz call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p(:,:,K), p(:,:,K+1), & - 0.0, G%HI, tv%eqn_of_state, dz_geo(:,:,k), halo_size=halo, & - SV_scale=US%R_to_kg_m3, pres_scale=US%RL2_T2_to_Pa) + 0.0, G%HI, tv%eqn_of_state, dz_geo(:,:,k), halo_size=halo, US=US) enddo !$OMP do do j=jsv,jev @@ -209,8 +208,7 @@ subroutine find_eta_2d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m) !$OMP do do k = 1, nz call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p(:,:,k), p(:,:,k+1), 0.0, & - G%HI, tv%eqn_of_state, dz_geo(:,:,k), halo_size=halo, & - SV_scale=US%R_to_kg_m3, pres_scale=US%RL2_T2_to_Pa) + G%HI, tv%eqn_of_state, dz_geo(:,:,k), halo_size=halo, US=US) enddo !$OMP do do j=js,je ; do k=1,nz ; do i=is,ie diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index be97723ee2..8c69853f3d 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -842,8 +842,7 @@ subroutine calculate_vertical_integrals(h, tv, p_surf, G, GV, US, CS) z_bot(i,j) = z_top(i,j) - GV%H_to_Z*h(i,j,k) enddo ; enddo call int_density_dz(tv%T(:,:,k), tv%S(:,:,k), z_top, z_bot, 0.0, GV%Rho0, GV%g_Earth, & - G%HI, G%HI, tv%eqn_of_state, dpress, rho_scale=US%kg_m3_to_R, & - pres_scale=US%RL2_T2_to_Pa) + G%HI, G%HI, tv%eqn_of_state, dpress, US=US) do j=js,je ; do i=is,ie mass(i,j) = mass(i,j) + dpress(i,j) * IG_Earth enddo ; enddo diff --git a/src/equation_of_state/MOM_EOS.F90 b/src/equation_of_state/MOM_EOS.F90 index d329b718bd..cfd286450e 100644 --- a/src/equation_of_state/MOM_EOS.F90 +++ b/src/equation_of_state/MOM_EOS.F90 @@ -29,8 +29,9 @@ module MOM_EOS use MOM_TFreeze, only : calculate_TFreeze_teos10 use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg use MOM_file_parser, only : get_param, log_version, param_file_type +use MOM_hor_index, only : hor_index_type use MOM_string_functions, only : uppercase -use MOM_hor_index, only : hor_index_type +use MOM_unit_scaling, only : unit_scale_type implicit none ; private @@ -58,19 +59,24 @@ module MOM_EOS !> Calculates density of sea water from T, S and P interface calculate_density - module procedure calculate_density_scalar, calculate_density_array + module procedure calculate_density_scalar, calculate_density_array, calculate_density_HI_1d end interface calculate_density !> Calculates specific volume of sea water from T, S and P interface calculate_spec_vol - module procedure calculate_spec_vol_scalar, calculate_spec_vol_array + module procedure calculate_spec_vol_scalar , calculate_spec_vol_array, calculate_spec_vol_HI_1d end interface calculate_spec_vol !> Calculate the derivatives of density with temperature and salinity from T, S, and P interface calculate_density_derivs - module procedure calculate_density_derivs_scalar, calculate_density_derivs_array + module procedure calculate_density_derivs_scalar, calculate_density_derivs_array, & + calculate_density_derivs_HI_1d end interface calculate_density_derivs +interface calculate_specific_vol_derivs + module procedure calculate_spec_vol_derivs_array, calculate_spec_vol_derivs_HI_1d +end interface calculate_specific_vol_derivs + !> Calculates the second derivatives of density with various combinations of temperature, !! salinity, and pressure from T, S and P interface calculate_density_second_derivs @@ -180,7 +186,7 @@ end subroutine calculate_density_scalar subroutine calculate_density_array(T, S, pressure, rho, start, npts, EOS, rho_ref, scale, pres_scale) real, dimension(:), intent(in) :: T !< Potential temperature referenced to the surface [degC] real, dimension(:), intent(in) :: S !< Salinity [ppt] - real, dimension(:), intent(in) :: pressure !< Pressure [Pa] + real, dimension(:), intent(in) :: pressure !< Pressure [Pa] or [other ~> Pa] real, dimension(:), intent(inout) :: rho !< Density (in-situ if pressure is local) [kg m-3] or [R ~> kg m-3] integer, intent(in) :: start !< Start index for computation integer, intent(in) :: npts !< Number of point to compute @@ -240,6 +246,94 @@ subroutine calculate_density_array(T, S, pressure, rho, start, npts, EOS, rho_re end subroutine calculate_density_array +!> Calls the appropriate subroutine to calculate the density of sea water for 1-D array inputs +!! using array extents determined from a hor_index_type. +!! If rho_ref is present, the anomaly with respect to rho_ref is returned. +subroutine calculate_density_HI_1d(T, S, pressure, rho, HI, EOS, US, halo, rho_ref) + type(hor_index_type), intent(in) :: HI !< The horizontal index structure + real, dimension(HI%isd:HI%ied), intent(in) :: T !< Potential temperature referenced to the surface [degC] + real, dimension(HI%isd:HI%ied), intent(in) :: S !< Salinity [ppt] + real, dimension(HI%isd:HI%ied), intent(in) :: pressure !< Pressure [R L2 T-2 ~> Pa] + real, dimension(HI%isd:HI%ied), intent(inout) :: rho !< Density (in-situ if pressure is local) [R ~> kg m-3] + type(EOS_type), pointer :: EOS !< Equation of state structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + integer, optional, intent(in) :: halo !< The halo size to work on; missing is equivalent to 0. + real, optional, intent(in) :: rho_ref !< A reference density [kg m-3]. + + ! Local variables + real, dimension(HI%isd:HI%ied) :: pres ! Pressure converted to [Pa] + real :: rho_reference ! rho_ref converted to [kg m-3] + integer :: i, is, ie, start, npts, halo_sz + + if (.not.associated(EOS)) call MOM_error(FATAL, & + "calculate_density_HI_1d called with an unassociated EOS_type EOS.") + + halo_sz = 0 ; if (present(halo)) halo_sz = halo + + start = HI%isc - (HI%isd-1) - halo_sz + npts = HI%iec - HI%isc + 1 + 2*halo_sz + is = HI%isc - halo_sz ; ie = HI%iec + halo_sz + + if ((US%RL2_T2_to_Pa == 1.0) .and. (US%R_to_kg_M3 == 1.0)) then + select case (EOS%form_of_EOS) + case (EOS_LINEAR) + call calculate_density_linear(T, S, pressure, rho, start, npts, & + EOS%Rho_T0_S0, EOS%dRho_dT, EOS%dRho_dS, rho_ref) + case (EOS_UNESCO) + call calculate_density_unesco(T, S, pressure, rho, start, npts, rho_ref) + case (EOS_WRIGHT) + call calculate_density_wright(T, S, pressure, rho, start, npts, rho_ref) + case (EOS_TEOS10) + call calculate_density_teos10(T, S, pressure, rho, start, npts, rho_ref) + case (EOS_NEMO) + call calculate_density_nemo(T, S, pressure, rho, start, npts, rho_ref) + case default + call MOM_error(FATAL, "calculate_density_HI_1d: EOS%form_of_EOS is not valid.") + end select + elseif (present(rho_ref)) then ! This is the same as above, but with some extra work to rescale variables. + do i=is,ie ; pres(i) = US%RL2_T2_to_Pa * pressure(i) ; enddo + rho_reference = 0.0 ; if (present(rho_ref)) rho_reference = US%R_to_kg_m3*rho_ref + select case (EOS%form_of_EOS) + case (EOS_LINEAR) + call calculate_density_linear(T, S, pres, rho, start, npts, & + EOS%Rho_T0_S0-rho_reference, EOS%dRho_dT, EOS%dRho_dS) + case (EOS_UNESCO) + call calculate_density_unesco(T, S, pres, rho, start, npts, rho_reference) + case (EOS_WRIGHT) + call calculate_density_wright(T, S, pres, rho, start, npts, rho_reference) + case (EOS_TEOS10) + call calculate_density_teos10(T, S, pres, rho, start, npts, rho_reference) + case (EOS_NEMO) + call calculate_density_nemo(T, S, pres, rho, start, npts, rho_reference) + case default + call MOM_error(FATAL, "calculate_density_HI_1d: EOS%form_of_EOS is not valid.") + end select + else ! There is rescaling of variables, but rho_ref is not present. Passing a 0 value of rho_ref + ! changes answers at roundoff for some equations of state, like Wright and UNESCO. + do i=is,ie ; pres(i) = US%RL2_T2_to_Pa * pressure(i) ; enddo + select case (EOS%form_of_EOS) + case (EOS_LINEAR) + call calculate_density_linear(T, S, pres, rho, start, npts, & + EOS%Rho_T0_S0, EOS%dRho_dT, EOS%dRho_dS) + case (EOS_UNESCO) + call calculate_density_unesco(T, S, pres, rho, start, npts) + case (EOS_WRIGHT) + call calculate_density_wright(T, S, pres, rho, start, npts) + case (EOS_TEOS10) + call calculate_density_teos10(T, S, pres, rho, start, npts) + case (EOS_NEMO) + call calculate_density_nemo(T, S, pres, rho, start, npts) + case default + call MOM_error(FATAL, "calculate_density_HI_1d: EOS%form_of_EOS is not valid.") + end select + endif + + if (US%kg_m3_to_R /= 1.0) then ; do i=is,ie + rho(i) = US%kg_m3_to_R * rho(i) + enddo ; endif + +end subroutine calculate_density_HI_1d + !> Calls the appropriate subroutine to calculate specific volume of sea water !! for scalar inputs. subroutine calculate_spec_vol_scalar(T, S, pressure, specvol, EOS, spv_ref, scale, pres_scale) @@ -288,7 +382,6 @@ subroutine calculate_spec_vol_scalar(T, S, pressure, specvol, EOS, spv_ref, scal end subroutine calculate_spec_vol_scalar - !> Calls the appropriate subroutine to calculate the specific volume of sea water !! for 1-D array inputs. subroutine calculate_spec_vol_array(T, S, pressure, specvol, start, npts, EOS, spv_ref, scale, pres_scale) @@ -365,6 +458,100 @@ subroutine calculate_spec_vol_array(T, S, pressure, specvol, start, npts, EOS, s end subroutine calculate_spec_vol_array +!> Calls the appropriate subroutine to calculate the specific volume of sea water for 1-D array +!! inputs using array extents determined from a hor_index_type. +subroutine calculate_spec_vol_HI_1d(T, S, pressure, specvol, HI, EOS, US, halo, spv_ref) + type(hor_index_type), intent(in) :: HI !< The horizontal index structure + real, dimension(HI%isd:HI%ied), intent(in) :: T !< Potential temperature referenced to the surface [degC] + real, dimension(HI%isd:HI%ied), intent(in) :: S !< Salinity [ppt] + real, dimension(HI%isd:HI%ied), intent(in) :: pressure !< Pressure [R L2 T-2 ~> Pa] + real, dimension(HI%isd:HI%ied), intent(inout) :: specvol !< In situ specific volume [R-1 ~> m3 kg-1] + type(EOS_type), pointer :: EOS !< Equation of state structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + integer, optional, intent(in) :: halo !< The halo size to work on; missing is equivalent to 0. + real, optional, intent(in) :: spv_ref !< A reference specific volume [R-1 ~> m3 kg-1]. + + ! Local variables + real, dimension(HI%isd:HI%ied) :: pres ! Pressure converted to [Pa] + real, dimension(HI%isd:HI%ied) :: rho ! Density [kg m-3] + real :: spv_reference ! spv_ref converted to [m3 kg-1] + integer :: i, is, ie, start, npts, halo_sz + + if (.not.associated(EOS)) call MOM_error(FATAL, & + "calculate_spec_vol_HI_1d called with an unassociated EOS_type EOS.") + + halo_sz = 0 ; if (present(halo)) halo_sz = halo + + start = HI%isc - (HI%isd-1) - halo_sz + npts = HI%iec - HI%isc + 1 + 2*halo_sz + is = HI%isc - halo_sz ; ie = HI%iec + halo_sz + + if ((US%RL2_T2_to_Pa == 1.0) .and. (US%R_to_kg_m3 == 1.0)) then + select case (EOS%form_of_EOS) + case (EOS_LINEAR) + call calculate_spec_vol_linear(T, S, pressure, specvol, start, npts, & + EOS%rho_T0_S0, EOS%drho_dT, EOS%drho_dS, spv_ref) + case (EOS_UNESCO) + call calculate_spec_vol_unesco(T, S, pressure, specvol, start, npts, spv_ref) + case (EOS_WRIGHT) + call calculate_spec_vol_wright(T, S, pressure, specvol, start, npts, spv_ref) + case (EOS_TEOS10) + call calculate_spec_vol_teos10(T, S, pressure, specvol, start, npts, spv_ref) + case (EOS_NEMO) + call calculate_density_nemo(T, S, pressure, rho, start, npts) + if (present(spv_ref)) then + specvol(:) = 1.0 / rho(:) - spv_ref + else + specvol(:) = 1.0 / rho(:) + endif + case default + call MOM_error(FATAL, "calculate_spec_vol_HI_1d: EOS%form_of_EOS is not valid.") + end select + elseif (present(spv_ref)) then ! This is the same as above, but with some extra work to rescale variables. + do i=is,ie ; pres(i) = US%RL2_T2_to_Pa * pressure(i) ; enddo + spv_reference = US%kg_m3_to_R*spv_ref + select case (EOS%form_of_EOS) + case (EOS_LINEAR) + call calculate_spec_vol_linear(T, S, pres, specvol, start, npts, & + EOS%rho_T0_S0, EOS%drho_dT, EOS%drho_dS, spv_reference) + case (EOS_UNESCO) + call calculate_spec_vol_unesco(T, S, pres, specvol, start, npts, spv_reference) + case (EOS_WRIGHT) + call calculate_spec_vol_wright(T, S, pres, specvol, start, npts, spv_reference) + case (EOS_TEOS10) + call calculate_spec_vol_teos10(T, S, pres, specvol, start, npts, spv_reference) + case (EOS_NEMO) + call calculate_density_nemo(T, S, pres, rho, start, npts) + specvol = 1.0 / rho - spv_reference + case default + call MOM_error(FATAL, "calculate_spec_vol_HI_1d: EOS%form_of_EOS is not valid.") + end select + else ! There is rescaling of variables, but spv_ref is not present. Passing a 0 value of spv_ref + ! changes answers at roundoff for some equations of state, like Wright and UNESCO. + do i=is,ie ; pres(i) = US%RL2_T2_to_Pa * pressure(i) ; enddo + select case (EOS%form_of_EOS) + case (EOS_LINEAR) + call calculate_spec_vol_linear(T, S, pres, specvol, start, npts, & + EOS%rho_T0_S0, EOS%drho_dT, EOS%drho_dS) + case (EOS_UNESCO) + call calculate_spec_vol_unesco(T, S, pres, specvol, start, npts) + case (EOS_WRIGHT) + call calculate_spec_vol_wright(T, S, pres, specvol, start, npts) + case (EOS_TEOS10) + call calculate_spec_vol_teos10(T, S, pres, specvol, start, npts) + case (EOS_NEMO) + call calculate_density_nemo(T, S, pres, rho, start, npts) + do i=is,ie ; specvol(i) = 1.0 / rho(i) ; enddo + case default + call MOM_error(FATAL, "calculate_spec_vol_HI_1d: EOS%form_of_EOS is not valid.") + end select + endif + + if (US%R_to_kg_m3 /= 1.0) then ; do i=is,ie + specvol(i) = US%R_to_kg_m3 * specvol(i) + enddo ; endif + +end subroutine calculate_spec_vol_HI_1d !> Calls the appropriate subroutine to calculate the freezing point for scalar inputs. subroutine calculate_TFreeze_scalar(S, pressure, T_fr, EOS, pres_scale) @@ -515,6 +702,78 @@ subroutine calculate_density_derivs_array(T, S, pressure, drho_dT, drho_dS, star end subroutine calculate_density_derivs_array + +!> Calls the appropriate subroutine to calculate density derivatives for 1-D array inputs. +subroutine calculate_density_derivs_HI_1d(T, S, pressure, drho_dT, drho_dS, HI, EOS, US, halo) + type(hor_index_type), intent(in) :: HI !< The horizontal index structure + real, dimension(HI%isd:HI%ied), intent(in) :: T !< Potential temperature referenced to the surface [degC] + real, dimension(HI%isd:HI%ied), intent(in) :: S !< Salinity [ppt] + real, dimension(HI%isd:HI%ied), intent(in) :: pressure !< Pressure [R L2 T-2 ~> Pa] + real, dimension(HI%isd:HI%ied), intent(inout) :: drho_dT !< The partial derivative of density with potential + !! temperature [R degC-1 ~> kg m-3 degC-1] + real, dimension(HI%isd:HI%ied), intent(inout) :: drho_dS !< The partial derivative of density with salinity + !! [R degC-1 ~> kg m-3 ppt-1] + type(EOS_type), pointer :: EOS !< Equation of state structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + integer, optional, intent(in) :: halo !< The halo size to work on; missing is equivalent to 0. + + ! Local variables + real, dimension(HI%isd:HI%ied) :: pres ! Pressure converted to [Pa] + real :: rho_reference ! rho_ref converted to [kg m-3] + integer :: i, is, ie, start, npts, halo_sz + + if (.not.associated(EOS)) call MOM_error(FATAL, & + "calculate_density_derivs called with an unassociated EOS_type EOS.") + + halo_sz = 0 ; if (present(halo)) halo_sz = halo + + start = HI%isc - (HI%isd-1) - halo_sz + npts = HI%iec - HI%isc + 1 + 2*halo_sz + is = HI%isc - halo_sz ; ie = HI%iec + halo_sz + + if (US%RL2_T2_to_Pa == 1.0) then + select case (EOS%form_of_EOS) + case (EOS_LINEAR) + call calculate_density_derivs_linear(T, S, pressure, drho_dT, drho_dS, EOS%Rho_T0_S0, & + EOS%dRho_dT, EOS%dRho_dS, start, npts) + case (EOS_UNESCO) + call calculate_density_derivs_unesco(T, S, pressure, drho_dT, drho_dS, start, npts) + case (EOS_WRIGHT) + call calculate_density_derivs_wright(T, S, pressure, drho_dT, drho_dS, start, npts) + case (EOS_TEOS10) + call calculate_density_derivs_teos10(T, S, pressure, drho_dT, drho_dS, start, npts) + case (EOS_NEMO) + call calculate_density_derivs_nemo(T, S, pressure, drho_dT, drho_dS, start, npts) + case default + call MOM_error(FATAL, "calculate_density_derivs_HI_1d: EOS%form_of_EOS is not valid.") + end select + else + do i=is,ie ; pres(i) = US%RL2_T2_to_Pa * pressure(i) ; enddo + select case (EOS%form_of_EOS) + case (EOS_LINEAR) + call calculate_density_derivs_linear(T, S, pres, drho_dT, drho_dS, EOS%Rho_T0_S0, & + EOS%dRho_dT, EOS%dRho_dS, start, npts) + case (EOS_UNESCO) + call calculate_density_derivs_unesco(T, S, pres, drho_dT, drho_dS, start, npts) + case (EOS_WRIGHT) + call calculate_density_derivs_wright(T, S, pres, drho_dT, drho_dS, start, npts) + case (EOS_TEOS10) + call calculate_density_derivs_teos10(T, S, pres, drho_dT, drho_dS, start, npts) + case (EOS_NEMO) + call calculate_density_derivs_nemo(T, S, pres, drho_dT, drho_dS, start, npts) + case default + call MOM_error(FATAL, "calculate_density_derivs_HI_1d: EOS%form_of_EOS is not valid.") + end select + endif + + if (US%kg_m3_to_R /= 1.0) then ; do i=is,ie + drho_dT(i) = US%kg_m3_to_R * drho_dT(i) + drho_dS(i) = US%kg_m3_to_R * drho_dS(i) + enddo ; endif + +end subroutine calculate_density_derivs_HI_1d + + !> Calls the appropriate subroutines to calculate density derivatives by promoting a scalar !! to a one-element array subroutine calculate_density_derivs_scalar(T, S, pressure, drho_dT, drho_dS, EOS, scale, pres_scale) @@ -620,8 +879,7 @@ subroutine calculate_density_second_derivs_array(T, S, pressure, drho_dS_dS, drh end select endif - if (present(scale)) then ; if (scale /= 1.0) then - ; do j=start,start+npts-1 + if (present(scale)) then ; if (scale /= 1.0) then ; do j=start,start+npts-1 drho_dS_dS(j) = scale * drho_dS_dS(j) drho_dS_dT(j) = scale * drho_dS_dT(j) drho_dT_dT(j) = scale * drho_dT_dT(j) @@ -699,7 +957,7 @@ subroutine calculate_density_second_derivs_scalar(T, S, pressure, drho_dS_dS, dr end subroutine calculate_density_second_derivs_scalar !> Calls the appropriate subroutine to calculate specific volume derivatives for an array. -subroutine calculate_specific_vol_derivs(T, S, pressure, dSV_dT, dSV_dS, start, npts, EOS, scale, pres_scale) +subroutine calculate_spec_vol_derivs_array(T, S, pressure, dSV_dT, dSV_dS, start, npts, EOS, scale, pres_scale) real, dimension(:), intent(in) :: T !< Potential temperature referenced to the surface [degC] real, dimension(:), intent(in) :: S !< Salinity [ppt] real, dimension(:), intent(in) :: pressure !< Pressure [Pa] @@ -715,48 +973,173 @@ subroutine calculate_specific_vol_derivs(T, S, pressure, dSV_dT, dSV_dS, start, real, optional, intent(in) :: pres_scale !< A multiplicative factor to convert pressure into Pa. ! Local variables - real, dimension(size(T)) :: dRho_dT, dRho_dS, rho + real, dimension(size(T)) :: press ! Pressure converted to [Pa] + real, dimension(size(T)) :: rho ! In situ density [kg m-3] + real, dimension(size(T)) :: dRho_dT ! Derivative of density with temperature [kg m-3 degC-1] + real, dimension(size(T)) :: dRho_dS ! Derivative of density with salinity [kg m-3 ppt-1] real :: p_scale ! A factor to convert pressure to units of Pa. integer :: j if (.not.associated(EOS)) call MOM_error(FATAL, & - "calculate_density_derivs called with an unassociated EOS_type EOS.") + "calculate_spec_vol_derivs_array called with an unassociated EOS_type EOS.") p_scale = 1.0 ; if (present(pres_scale)) p_scale = pres_scale - select case (EOS%form_of_EOS) - case (EOS_LINEAR) - call calculate_specvol_derivs_linear(T, S, pressure, dSV_dT, dSV_dS, start, & - npts, EOS%Rho_T0_S0, EOS%dRho_dT, EOS%dRho_dS) - case (EOS_UNESCO) - call calculate_density_unesco(T, S, pressure, rho, start, npts) - call calculate_density_derivs_unesco(T, S, pressure, drho_dT, drho_dS, start, npts) - do j=start,start+npts-1 - dSV_dT(j) = -dRho_DT(j)/(rho(j)**2) - dSV_dS(j) = -dRho_DS(j)/(rho(j)**2) - enddo - case (EOS_WRIGHT) - call calculate_specvol_derivs_wright(T, S, pressure, dSV_dT, dSV_dS, start, npts) - case (EOS_TEOS10) - call calculate_specvol_derivs_teos10(T, S, pressure, dSV_dT, dSV_dS, start, npts) - case (EOS_NEMO) - call calculate_density_nemo(T, S, pressure, rho, start, npts) - call calculate_density_derivs_nemo(T, S, pressure, drho_dT, drho_dS, start, npts) - do j=start,start+npts-1 - dSV_dT(j) = -dRho_DT(j)/(rho(j)**2) - dSV_dS(j) = -dRho_DS(j)/(rho(j)**2) - enddo - case default - call MOM_error(FATAL, "calculate_density_derivs: EOS%form_of_EOS is not valid.") - end select + if (p_scale == 1.0) then + select case (EOS%form_of_EOS) + case (EOS_LINEAR) + call calculate_specvol_derivs_linear(T, S, pressure, dSV_dT, dSV_dS, start, & + npts, EOS%Rho_T0_S0, EOS%dRho_dT, EOS%dRho_dS) + case (EOS_UNESCO) + call calculate_density_unesco(T, S, pressure, rho, start, npts) + call calculate_density_derivs_unesco(T, S, pressure, drho_dT, drho_dS, start, npts) + do j=start,start+npts-1 + dSV_dT(j) = -dRho_DT(j)/(rho(j)**2) + dSV_dS(j) = -dRho_DS(j)/(rho(j)**2) + enddo + case (EOS_WRIGHT) + call calculate_specvol_derivs_wright(T, S, pressure, dSV_dT, dSV_dS, start, npts) + case (EOS_TEOS10) + call calculate_specvol_derivs_teos10(T, S, pressure, dSV_dT, dSV_dS, start, npts) + case (EOS_NEMO) + call calculate_density_nemo(T, S, pressure, rho, start, npts) + call calculate_density_derivs_nemo(T, S, pressure, drho_dT, drho_dS, start, npts) + do j=start,start+npts-1 + dSV_dT(j) = -dRho_DT(j)/(rho(j)**2) + dSV_dS(j) = -dRho_DS(j)/(rho(j)**2) + enddo + case default + call MOM_error(FATAL, "calculate_spec_vol_derivs_array: EOS%form_of_EOS is not valid.") + end select + else + do j=start,start+npts-1 ; press(j) = p_scale * pressure(j) ; enddo + select case (EOS%form_of_EOS) + case (EOS_LINEAR) + call calculate_specvol_derivs_linear(T, S, press, dSV_dT, dSV_dS, start, & + npts, EOS%Rho_T0_S0, EOS%dRho_dT, EOS%dRho_dS) + case (EOS_UNESCO) + call calculate_density_unesco(T, S, press, rho, start, npts) + call calculate_density_derivs_unesco(T, S, press, drho_dT, drho_dS, start, npts) + do j=start,start+npts-1 + dSV_dT(j) = -dRho_DT(j)/(rho(j)**2) + dSV_dS(j) = -dRho_DS(j)/(rho(j)**2) + enddo + case (EOS_WRIGHT) + call calculate_specvol_derivs_wright(T, S, press, dSV_dT, dSV_dS, start, npts) + case (EOS_TEOS10) + call calculate_specvol_derivs_teos10(T, S, press, dSV_dT, dSV_dS, start, npts) + case (EOS_NEMO) + call calculate_density_nemo(T, S, pressure, rho, start, npts) + call calculate_density_derivs_nemo(T, S, press, drho_dT, drho_dS, start, npts) + do j=start,start+npts-1 + dSV_dT(j) = -dRho_DT(j)/(rho(j)**2) + dSV_dS(j) = -dRho_DS(j)/(rho(j)**2) + enddo + case default + call MOM_error(FATAL, "calculate_spec_vol_derivs_array: EOS%form_of_EOS is not valid.") + end select + endif if (present(scale)) then ; if (scale /= 1.0) then ; do j=start,start+npts-1 dSV_dT(j) = scale * dSV_dT(j) dSV_dS(j) = scale * dSV_dS(j) enddo ; endif ; endif +end subroutine calculate_spec_vol_derivs_array + +!> Calls the appropriate subroutine to calculate specific volume derivatives for an array. +subroutine calculate_spec_vol_derivs_HI_1d(T, S, pressure, dSV_dT, dSV_dS, HI, EOS, US, halo) + type(hor_index_type), intent(in) :: HI !< The horizontal index structure + real, dimension(HI%isd:HI%ied), intent(in) :: T !< Potential temperature referenced to the surface [degC] + real, dimension(HI%isd:HI%ied), intent(in) :: S !< Salinity [ppt] + real, dimension(HI%isd:HI%ied), intent(in) :: pressure !< Pressure [R L2 T-2 ~> Pa] + real, dimension(HI%isd:HI%ied), intent(inout) :: dSV_dT !< The partial derivative of specific volume with potential + !! temperature [R-1 degC-1 ~> m3 kg-1 degC-1] + real, dimension(HI%isd:HI%ied), intent(inout) :: dSV_dS !< The partial derivative of specific volume with salinity + !! [R-1 ppt-1 ~> m3 kg-1 ppt-1]. + type(EOS_type), pointer :: EOS !< Equation of state structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + integer, optional, intent(in) :: halo !< The halo size to work on; missing is equivalent to 0. + + ! Local variables + real, dimension(HI%isd:HI%ied) :: rho ! In situ density [kg m-3] + real, dimension(HI%isd:HI%ied) :: dRho_dT ! Derivative of density with temperature [kg m-3 degC-1] + real, dimension(HI%isd:HI%ied) :: dRho_dS ! Derivative of density with salinity [kg m-3 ppt-1] + real, dimension(HI%isd:HI%ied) :: press ! Pressure converted to [Pa] + real :: rho_reference ! rho_ref converted to [kg m-3] + integer :: i, is, ie, start, npts, halo_sz + + if (.not.associated(EOS)) call MOM_error(FATAL, & + "calculate_spec_vol_derivs_HI_1d called with an unassociated EOS_type EOS.") + + halo_sz = 0 ; if (present(halo)) halo_sz = halo + + start = HI%isc - (HI%isd-1) - halo_sz + npts = HI%iec - HI%isc + 1 + 2*halo_sz + is = HI%isc - halo_sz ; ie = HI%iec + halo_sz + + if (US%RL2_T2_to_Pa == 1.0) then + select case (EOS%form_of_EOS) + case (EOS_LINEAR) + call calculate_specvol_derivs_linear(T, S, pressure, dSV_dT, dSV_dS, start, & + npts, EOS%Rho_T0_S0, EOS%dRho_dT, EOS%dRho_dS) + case (EOS_UNESCO) + call calculate_density_unesco(T, S, pressure, rho, start, npts) + call calculate_density_derivs_unesco(T, S, pressure, drho_dT, drho_dS, start, npts) + do i=is,ie + dSV_dT(i) = -dRho_DT(i)/(rho(i)**2) + dSV_dS(i) = -dRho_DS(i)/(rho(i)**2) + enddo + case (EOS_WRIGHT) + call calculate_specvol_derivs_wright(T, S, pressure, dSV_dT, dSV_dS, start, npts) + case (EOS_TEOS10) + call calculate_specvol_derivs_teos10(T, S, pressure, dSV_dT, dSV_dS, start, npts) + case (EOS_NEMO) + call calculate_density_nemo(T, S, pressure, rho, start, npts) + call calculate_density_derivs_nemo(T, S, pressure, drho_dT, drho_dS, start, npts) + do i=is,ie + dSV_dT(i) = -dRho_DT(i)/(rho(i)**2) + dSV_dS(i) = -dRho_DS(i)/(rho(i)**2) + enddo + case default + call MOM_error(FATAL, "calculate_spec_vol_derivs_HI_1d: EOS%form_of_EOS is not valid.") + end select + else + do i=is,ie ; press(i) = US%RL2_T2_to_Pa * pressure(i) ; enddo + select case (EOS%form_of_EOS) + case (EOS_LINEAR) + call calculate_specvol_derivs_linear(T, S, press, dSV_dT, dSV_dS, start, & + npts, EOS%Rho_T0_S0, EOS%dRho_dT, EOS%dRho_dS) + case (EOS_UNESCO) + call calculate_density_unesco(T, S, press, rho, start, npts) + call calculate_density_derivs_unesco(T, S, press, drho_dT, drho_dS, start, npts) + do i=is,ie + dSV_dT(i) = -dRho_DT(i)/(rho(i)**2) + dSV_dS(i) = -dRho_DS(i)/(rho(i)**2) + enddo + case (EOS_WRIGHT) + call calculate_specvol_derivs_wright(T, S, press, dSV_dT, dSV_dS, start, npts) + case (EOS_TEOS10) + call calculate_specvol_derivs_teos10(T, S, press, dSV_dT, dSV_dS, start, npts) + case (EOS_NEMO) + call calculate_density_nemo(T, S, pressure, rho, start, npts) + call calculate_density_derivs_nemo(T, S, press, drho_dT, drho_dS, start, npts) + do i=is,ie + dSV_dT(i) = -dRho_DT(i)/(rho(i)**2) + dSV_dS(i) = -dRho_DS(i)/(rho(i)**2) + enddo + case default + call MOM_error(FATAL, "calculate_spec_vol_derivs_HI_1d: EOS%form_of_EOS is not valid.") + end select + endif + + if (US%R_to_kg_m3 /= 1.0) then ; do i=is,ie + drho_dT(i) = US%R_to_kg_m3 * drho_dT(i) + drho_dS(i) = US%R_to_kg_m3 * drho_dS(i) + enddo ; endif + +end subroutine calculate_spec_vol_derivs_HI_1d -end subroutine calculate_specific_vol_derivs !> Calls the appropriate subroutine to calculate the density and compressibility for 1-D array inputs. subroutine calculate_compress_array(T, S, pressure, rho, drho_dp, start, npts, EOS) @@ -791,15 +1174,15 @@ subroutine calculate_compress_array(T, S, pressure, rho, drho_dp, start, npts, E end subroutine calculate_compress_array -!> Calculate density and compressibility for a scalar. This just promotes the scalar to an array with a singleton -!! dimension and calls calculate_compress_array +!> Calculate density and compressibility for a scalar. This just promotes the scalar to an array +!! with a singleton dimension and calls calculate_compress_array subroutine calculate_compress_scalar(T, S, pressure, rho, drho_dp, EOS) - real, intent(in) :: T !< Potential temperature referenced to the surface (degC) - real, intent(in) :: S !< Salinity (PSU) - real, intent(in) :: pressure !< Pressure (Pa) - real, intent(out) :: rho !< In situ density in kg m-3. + real, intent(in) :: T !< Potential temperature referenced to the surface [degC] + real, intent(in) :: S !< Salinity [ppt] + real, intent(in) :: pressure !< Pressure [Pa] + real, intent(out) :: rho !< In situ density [kg m-3] real, intent(out) :: drho_dp !< The partial derivative of density with pressure - !! (also the inverse of the square of sound speed) in s2 m-2. + !! (also the inverse of the square of sound speed) [s2 m-2]. type(EOS_type), pointer :: EOS !< Equation of state structure real, dimension(1) :: Ta, Sa, pa, rhoa, drho_dpa @@ -820,7 +1203,7 @@ end subroutine calculate_compress_scalar !! series for log(1-eps/1+eps) that assumes that |eps| < . subroutine int_specific_vol_dp(T, S, p_t, p_b, alpha_ref, HI, EOS, & dza, intp_dza, intx_dza, inty_dza, halo_size, & - bathyP, dP_tiny, useMassWghtInterp, SV_scale, pres_scale) + bathyP, dP_tiny, useMassWghtInterp, US) type(hor_index_type), intent(in) :: HI !< The horizontal index structure real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & intent(in) :: T !< Potential temperature referenced to the surface [degC] @@ -857,13 +1240,12 @@ subroutine int_specific_vol_dp(T, S, p_t, p_b, alpha_ref, HI, EOS, & !! the same units as p_t [R L2 T-2 ~> Pa] or [Pa] logical, optional, intent(in) :: useMassWghtInterp !< If true, uses mass weighting !! to interpolate T/S for top and bottom integrals. - real, optional, intent(in) :: SV_scale !< A multiplicative factor by which to scale specific - !! volume from m3 kg-1 to the desired units [kg m-3 R-1 ~> 1] - real, optional, intent(in) :: pres_scale !< A multiplicative factor to convert pressure - !! into Pa [Pa T2 R-1 L-2 ~> 1]. + type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type ! Local variables - real :: rho_scale ! A scaling factor for densities from kg m-3 to R [R m3 kg-1 ~> 1] + real :: pres_scale ! A unit conversion factor from the rescaled units of pressure to Pa [Pa T2 R-1 L-2 ~> 1] + real :: SV_scale ! A multiplicative factor by which to scale specific + ! volume from m3 kg-1 to the desired units [kg m-3 R-1 ~> 1] if (.not.associated(EOS)) call MOM_error(FATAL, & "int_specific_vol_dp called with an unassociated EOS_type EOS.") @@ -871,14 +1253,13 @@ subroutine int_specific_vol_dp(T, S, p_t, p_b, alpha_ref, HI, EOS, & if (EOS%EOS_quadrature) then call int_spec_vol_dp_generic(T, S, p_t, p_b, alpha_ref, HI, EOS, & dza, intp_dza, intx_dza, inty_dza, halo_size, & - bathyP, dP_tiny, useMassWghtInterp, SV_scale, pres_scale) + bathyP, dP_tiny, useMassWghtInterp, US) else ; select case (EOS%form_of_EOS) case (EOS_LINEAR) - if (present(SV_scale)) then - rho_scale = 1.0 / SV_scale - call int_spec_vol_dp_linear(T, S, p_t, p_b, alpha_ref, HI, rho_scale*EOS%Rho_T0_S0, & - rho_scale*EOS%dRho_dT, rho_scale*EOS%dRho_dS, dza, intp_dza, & - intx_dza, inty_dza, halo_size, & + if (present(US)) then + call int_spec_vol_dp_linear(T, S, p_t, p_b, alpha_ref, HI, US%kg_m3_to_R*EOS%Rho_T0_S0, & + US%kg_m3_to_R*EOS%dRho_dT, US%kg_m3_to_R*EOS%dRho_dS, dza, & + intp_dza, intx_dza, inty_dza, halo_size, & bathyP, dP_tiny, useMassWghtInterp) else call int_spec_vol_dp_linear(T, S, p_t, p_b, alpha_ref, HI, EOS%Rho_T0_S0, & @@ -887,13 +1268,18 @@ subroutine int_specific_vol_dp(T, S, p_t, p_b, alpha_ref, HI, EOS, & bathyP, dP_tiny, useMassWghtInterp) endif case (EOS_WRIGHT) - call int_spec_vol_dp_wright(T, S, p_t, p_b, alpha_ref, HI, dza, & - intp_dza, intx_dza, inty_dza, halo_size, & - bathyP, dP_tiny, useMassWghtInterp, SV_scale, pres_scale) + if (present(US)) then + call int_spec_vol_dp_wright(T, S, p_t, p_b, alpha_ref, HI, dza, intp_dza, intx_dza, & + inty_dza, halo_size, bathyP, dP_tiny, useMassWghtInterp, & + SV_scale=US%R_to_kg_m3, pres_scale=US%RL2_T2_to_Pa) + else + call int_spec_vol_dp_wright(T, S, p_t, p_b, alpha_ref, HI, dza, intp_dza, intx_dza, & + inty_dza, halo_size, bathyP, dP_tiny, useMassWghtInterp) + endif case default call int_spec_vol_dp_generic(T, S, p_t, p_b, alpha_ref, HI, EOS, & dza, intp_dza, intx_dza, inty_dza, halo_size, & - bathyP, dP_tiny, useMassWghtInterp, SV_scale, pres_scale) + bathyP, dP_tiny, useMassWghtInterp, US) end select ; endif end subroutine int_specific_vol_dp @@ -902,7 +1288,7 @@ end subroutine int_specific_vol_dp !! pressure anomalies across layers, which are required for calculating the !! finite-volume form pressure accelerations in a Boussinesq model. subroutine int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, EOS, dpa, & - intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp, rho_scale, pres_scale) + intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp, US) type(hor_index_type), intent(in) :: HII !< Ocean horizontal index structures for the input arrays type(hor_index_type), intent(in) :: HIO !< Ocean horizontal index structures for the output arrays real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), & @@ -942,21 +1328,23 @@ subroutine int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, EOS, dp real, optional, intent(in) :: dz_neglect !< A miniscule thickness change [Z ~> m]. logical, optional, intent(in) :: useMassWghtInterp !< If true, uses mass weighting to !! interpolate T/S for top and bottom integrals. - real, optional, intent(in) :: rho_scale !< A multiplicative factor by which to scale density - !! from kg m-3 to the desired units [R m3 kg-1 ~> 1] - real, optional, intent(in) :: pres_scale !< A multiplicative factor to convert pressure - !! into Pa [R L2 T-2 Pa-1 ~> 1]. + type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type + + ! Local variables + real :: rho_scale ! A multiplicative factor by which to scale density from kg m-3 to the + ! desired units [R m3 kg-1 ~> 1]. + real :: pres_scale ! A multiplicative factor to convert pressure into Pa [Pa T2 R-1 L-2 ~> 1]. if (.not.associated(EOS)) call MOM_error(FATAL, & "int_density_dz called with an unassociated EOS_type EOS.") if (EOS%EOS_quadrature) then call int_density_dz_generic(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, EOS, dpa, & - intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp, & - rho_scale, pres_scale) + intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp, US) else ; select case (EOS%form_of_EOS) case (EOS_LINEAR) - if (present(rho_scale)) then + rho_scale = 1.0 ; if (present(US)) rho_scale = US%kg_m3_to_R + if (rho_scale /= 1.0) then call int_density_dz_linear(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, & rho_scale*EOS%Rho_T0_S0, rho_scale*EOS%dRho_dT, rho_scale*EOS%dRho_dS, & dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp) @@ -966,13 +1354,20 @@ subroutine int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, EOS, dp dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp) endif case (EOS_WRIGHT) - call int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, & - dpa, intz_dpa, intx_dpa, inty_dpa, & - bathyT, dz_neglect, useMassWghtInterp, rho_scale, pres_scale) + rho_scale = 1.0 ; if (present(US)) rho_scale = US%kg_m3_to_R + pres_scale = 1.0 ; if (present(US)) pres_scale = US%RL2_T2_to_Pa + if ((rho_scale /= 1.0) .or. (pres_scale /= 1.0)) then + call int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, & + dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, & + dz_neglect, useMassWghtInterp, rho_scale, pres_scale) + else + call int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, & + dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, & + dz_neglect, useMassWghtInterp) + endif case default call int_density_dz_generic(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, EOS, dpa, & - intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp, & - rho_scale, pres_scale) + intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp, US) end select ; endif end subroutine int_density_dz @@ -1166,7 +1561,7 @@ end subroutine EOS_use_linear !! finite-volume form pressure accelerations in a Boussinesq model. subroutine int_density_dz_generic(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, & EOS, dpa, intz_dpa, intx_dpa, inty_dpa, & - bathyT, dz_neglect, useMassWghtInterp, rho_scale, pres_scale) + bathyT, dz_neglect, useMassWghtInterp, US) type(hor_index_type), intent(in) :: HII !< Horizontal index type for input variables. type(hor_index_type), intent(in) :: HIO !< Horizontal index type for output variables. real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), & @@ -1187,18 +1582,18 @@ subroutine int_density_dz_generic(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, !! [L2 Z-1 T-2 ~> m s-2] or [m2 Z-1 s-2 ~> m s-2]. type(EOS_type), pointer :: EOS !< Equation of state structure real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), & - intent(inout) :: dpa !< The change in the pressure anomaly + intent(inout) :: dpa !< The change in the pressure anomaly !! across the layer [R L2 T-2 ~> Pa] or [Pa]. real, dimension(HIO%isd:HIO%ied,HIO%jsd:HIO%jed), & - optional, intent(inout) :: intz_dpa !< The integral through the thickness of the + optional, intent(inout) :: intz_dpa !< The integral through the thickness of the !! layer of the pressure anomaly relative to the !! anomaly at the top of the layer [R L2 Z T-2 ~> Pa m]. real, dimension(HIO%IsdB:HIO%IedB,HIO%jsd:HIO%jed), & - optional, intent(inout) :: intx_dpa !< The integral in x of the difference between + optional, intent(inout) :: intx_dpa !< The integral in x of the difference between !! the pressure anomaly at the top and bottom of the !! layer divided by the x grid spacing [R L2 T-2 ~> Pa]. real, dimension(HIO%isd:HIO%ied,HIO%JsdB:HIO%JedB), & - optional, intent(inout) :: inty_dpa !< The integral in y of the difference between + optional, intent(inout) :: inty_dpa !< The integral in y of the difference between !! the pressure anomaly at the top and bottom of the !! layer divided by the y grid spacing [R L2 T-2 ~> Pa]. real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), & @@ -1206,10 +1601,8 @@ subroutine int_density_dz_generic(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, real, optional, intent(in) :: dz_neglect !< A miniscule thickness change [Z ~> m]. logical, optional, intent(in) :: useMassWghtInterp !< If true, uses mass weighting to !! interpolate T/S for top and bottom integrals. - real, optional, intent(in) :: rho_scale !< A multiplicative factor by which to scale density - !! from kg m-3 to the desired units [R m3 kg-1 ~> 1] - real, optional, intent(in) :: pres_scale !< A multiplicative factor to convert pressure - !! into Pa [Pa T2 R-1 L-2 ~> 1]. + type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type + ! Local variables real :: T5(5), S5(5) ! Temperatures and salinities at five quadrature points [degC] and [ppt] real :: p5(5) ! Pressures at five quadrature points, never rescaled from Pa [Pa] @@ -1219,6 +1612,7 @@ subroutine int_density_dz_generic(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, real, parameter :: C1_90 = 1.0/90.0 ! Rational constants. real :: GxRho ! The gravitational acceleration times density and unit conversion factors [Pa Z-1 ~> kg m-2 s-2] real :: I_Rho ! The inverse of the Boussinesq density [R-1 ~> m3 kg-1] or [m3 kg-1] + real :: rho_scale ! A scaling factor for densities from kg m-3 to R [R m3 kg-1 ~> 1] real :: rho_ref_mks ! The reference density in MKS units, never rescaled from kg m-3 [kg m-3] real :: dz ! The layer thickness [Z ~> m]. real :: hWght ! A pressure-thickness below topography [Z ~> m]. @@ -1243,8 +1637,9 @@ subroutine int_density_dz_generic(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, is = HIO%isc + ioff ; ie = HIO%iec + ioff js = HIO%jsc + joff ; je = HIO%jec + joff - GxRho = G_e * rho_0 ; if (present(pres_scale)) GxRho = pres_scale * G_e * rho_0 - rho_ref_mks = rho_ref ; if (present(rho_scale)) rho_ref_mks = rho_ref / rho_scale + rho_scale = 1.0 ; if (present(US)) rho_scale = US%kg_m3_to_R + GxRho = G_e * rho_0 ; if (present(US)) GxRho = US%RL2_T2_to_Pa * G_e * rho_0 + rho_ref_mks = rho_ref ; if (present(US)) rho_ref_mks = rho_ref * US%R_to_kg_m3 I_Rho = 1.0 / rho_0 do_massWeight = .false. @@ -1262,7 +1657,11 @@ subroutine int_density_dz_generic(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, T5(n) = T(i,j) ; S5(n) = S(i,j) p5(n) = -GxRho*(z_t(i,j) - 0.25*real(n-1)*dz) enddo - call calculate_density(T5, S5, p5, r5, 1, 5, EOS, rho_ref_mks, scale=rho_scale) + if (rho_scale /= 1.0) then + call calculate_density(T5, S5, p5, r5, 1, 5, EOS, rho_ref_mks, scale=rho_scale) + else + call calculate_density(T5, S5, p5, r5, 1, 5, EOS, rho_ref_mks) + endif ! Use Bode's rule to estimate the pressure anomaly change. rho_anom = C1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3)) @@ -1304,7 +1703,11 @@ subroutine int_density_dz_generic(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, do n=2,5 T5(n) = T5(1) ; S5(n) = S5(1) ; p5(n) = p5(n-1) + GxRho*0.25*dz enddo - call calculate_density(T5, S5, p5, r5, 1, 5, EOS, rho_ref_mks, scale=rho_scale) + if (rho_scale /= 1.0) then + call calculate_density(T5, S5, p5, r5, 1, 5, EOS, rho_ref_mks, scale=rho_scale) + else + call calculate_density(T5, S5, p5, r5, 1, 5, EOS, rho_ref_mks) + endif ! Use Bode's rule to estimate the pressure anomaly change. intz(m) = G_e*dz*( C1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3))) @@ -1346,7 +1749,11 @@ subroutine int_density_dz_generic(T, S, z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, T5(n) = T5(1) ; S5(n) = S5(1) p5(n) = p5(n-1) + GxRho*0.25*dz enddo - call calculate_density(T5, S5, p5, r5, 1, 5, EOS, rho_ref_mks, scale=rho_scale) + if (rho_scale /= 1.0) then + call calculate_density(T5, S5, p5, r5, 1, 5, EOS, rho_ref_mks, scale=rho_scale) + else + call calculate_density(T5, S5, p5, r5, 1, 5, EOS, rho_ref_mks) + endif ! Use Bode's rule to estimate the pressure anomaly change. intz(m) = G_e*dz*( C1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3))) @@ -1363,8 +1770,7 @@ end subroutine int_density_dz_generic !! T and S are linear profiles. subroutine int_density_dz_generic_plm (T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, & rho_0, G_e, dz_subroundoff, bathyT, HII, HIO, EOS, dpa, & - intz_dpa, intx_dpa, inty_dpa, & - useMassWghtInterp, rho_scale, pres_scale) + intz_dpa, intx_dpa, inty_dpa, useMassWghtInterp, US) type(hor_index_type), intent(in) :: HII !< Ocean horizontal index structures for the input arrays type(hor_index_type), intent(in) :: HIO !< Ocean horizontal index structures for the output arrays real, dimension(HII%isd:HII%ied,HII%jsd:HII%jed), & @@ -1404,10 +1810,8 @@ subroutine int_density_dz_generic_plm (T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, & !! divided by the y grid spacing [R L2 T-2 ~> Pa]. logical, optional, intent(in) :: useMassWghtInterp !< If true, uses mass weighting to !! interpolate T/S for top and bottom integrals. - real, optional, intent(in) :: rho_scale !< A multiplicative factor by which to scale density - !! from kg m-3 to the desired units [R m3 kg-1 ~> 1] - real, optional, intent(in) :: pres_scale !< A multiplicative factor to convert pressure - !! into Pa [Pa T2 R-1 L-2 ~> 1]. + type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type + ! This subroutine calculates (by numerical quadrature) integrals of ! pressure anomalies across layers, which are required for calculating the ! finite-volume form pressure accelerations in a Boussinesq model. The one @@ -1439,6 +1843,7 @@ subroutine int_density_dz_generic_plm (T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, & real, parameter :: C1_90 = 1.0/90.0 ! A rational constant [nondim]. real :: GxRho ! The gravitational acceleration times density and unit conversion factors [Pa Z-1 ~> kg m-2 s-2] real :: I_Rho ! The inverse of the Boussinesq density [R-1 ~> m3 kg-1] or [m3 kg-1] + real :: rho_scale ! A scaling factor for densities from kg m-3 to R [R m3 kg-1 ~> 1] real :: rho_ref_mks ! The reference density in MKS units, never rescaled from kg m-3 [kg m-3] real :: dz(HIO%iscB:HIO%iecB+1) ! Layer thicknesses at tracer points [Z ~> m]. real :: dz_x(5,HIO%iscB:HIO%iecB) ! Layer thicknesses along an x-line of subrid locations [Z ~> m]. @@ -1459,8 +1864,9 @@ subroutine int_density_dz_generic_plm (T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, & Isq = HIO%IscB ; Ieq = HIO%IecB ; Jsq = HIO%JscB ; Jeq = HIO%JecB - GxRho = G_e * rho_0 ; if (present(pres_scale)) GxRho = pres_scale * G_e * rho_0 - rho_ref_mks = rho_ref ; if (present(rho_scale)) rho_ref_mks = rho_ref / rho_scale + rho_scale = 1.0 ; if (present(US)) rho_scale = US%kg_m3_to_R + GxRho = G_e * rho_0 ; if (present(US)) GxRho = US%RL2_T2_to_Pa * G_e * rho_0 + rho_ref_mks = rho_ref ; if (present(US)) rho_ref_mks = rho_ref * US%R_to_kg_m3 I_Rho = 1.0 / rho_0 massWeightToggle = 0. if (present(useMassWghtInterp)) then @@ -1486,7 +1892,11 @@ subroutine int_density_dz_generic_plm (T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, & T5(i*5+n) = wt_t(n) * T_t(iin,jin) + wt_b(n) * T_b(iin,jin) enddo enddo - call calculate_density_array(T5, S5, p5, r5, 1, (ieq-isq+2)*5, EOS, rho_ref_mks, scale=rho_scale) + if (rho_scale /= 1.0) then + call calculate_density_array(T5, S5, p5, r5, 1, (ieq-isq+2)*5, EOS, rho_ref_mks, scale=rho_scale) + else + call calculate_density_array(T5, S5, p5, r5, 1, (ieq-isq+2)*5, EOS, rho_ref_mks) + endif do i=isq,ieq+1 ; iin = i+ioff ! Use Bode's rule to estimate the pressure anomaly change. @@ -1566,7 +1976,11 @@ subroutine int_density_dz_generic_plm (T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, & enddo enddo - call calculate_density(T15, S15, p15, r15, 1, 15*(ieq-isq+1), EOS, rho_ref_mks, scale=rho_scale) + if (rho_scale /= 1.0) then + call calculate_density(T15, S15, p15, r15, 1, 15*(ieq-isq+1), EOS, rho_ref_mks, scale=rho_scale) + else + call calculate_density(T15, S15, p15, r15, 1, 15*(ieq-isq+1), EOS, rho_ref_mks) + endif do I=Isq,Ieq ; iin = i+ioff intz(1) = dpa(i,j) ; intz(5) = dpa(i+1,j) @@ -1645,8 +2059,13 @@ subroutine int_density_dz_generic_plm (T_t, T_b, S_t, S_b, z_t, z_b, rho_ref, & enddo enddo - call calculate_density_array(T15(15*HIO%isc+1:), S15(15*HIO%isc+1:), p15(15*HIO%isc+1:), & - r15(15*HIO%isc+1:), 1, 15*(HIO%iec-HIO%isc+1), EOS, rho_ref_mks, scale=rho_scale) + if (rho_scale /= 1.0) then + call calculate_density_array(T15(15*HIO%isc+1:), S15(15*HIO%isc+1:), p15(15*HIO%isc+1:), & + r15(15*HIO%isc+1:), 1, 15*(HIO%iec-HIO%isc+1), EOS, rho_ref_mks, scale=rho_scale) + else + call calculate_density_array(T15(15*HIO%isc+1:), S15(15*HIO%isc+1:), p15(15*HIO%isc+1:), & + r15(15*HIO%isc+1:), 1, 15*(HIO%iec-HIO%isc+1), EOS, rho_ref_mks) + endif do i=HIO%isc,HIO%iec ; iin = i+ioff intz(1) = dpa(i,j) ; intz(5) = dpa(i,j+1) @@ -1784,9 +2203,9 @@ end function frac_dp_at_pos ! ========================================================================== !> Compute pressure gradient force integrals for the case where T and S !! are parabolic profiles -subroutine int_density_dz_generic_ppm (T, T_t, T_b, S, S_t, S_b, & - z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, & - EOS, dpa, intz_dpa, intx_dpa, inty_dpa, rho_scale, pres_scale) +subroutine int_density_dz_generic_ppm(T, T_t, T_b, S, S_t, S_b, & + z_t, z_b, rho_ref, rho_0, G_e, HII, HIO, & + EOS, dpa, intz_dpa, intx_dpa, inty_dpa, US) type(hor_index_type), intent(in) :: HII !< Ocean horizontal index structures for the input arrays type(hor_index_type), intent(in) :: HIO !< Ocean horizontal index structures for the output arrays @@ -1826,10 +2245,7 @@ subroutine int_density_dz_generic_ppm (T, T_t, T_b, S, S_t, S_b, & optional, intent(inout) :: inty_dpa !< The integral in y of the difference between the !! pressure anomaly at the top and bottom of the layer !! divided by the y grid spacing [Pa]. - real, optional, intent(in) :: rho_scale !< A multiplicative factor by which to scale density - !! from kg m-3 to the desired units [R m3 kg-1 ~> 1] - real, optional, intent(in) :: pres_scale !< A multiplicative factor to convert pressure - !! into Pa [Pa T2 R-1 L-2 ~> 1]. + type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type ! This subroutine calculates (by numerical quadrature) integrals of ! pressure anomalies across layers, which are required for calculating the @@ -1855,6 +2271,7 @@ subroutine int_density_dz_generic_ppm (T, T_t, T_b, S, S_t, S_b, & real, parameter :: C1_90 = 1.0/90.0 ! Rational constants. real :: GxRho ! The gravitational acceleration times density and unit conversion factors [Pa Z-1 ~> kg m-2 s-2] real :: I_Rho ! The inverse of the Boussinesq density [R-1 ~> m3 kg-1] or [m3 kg-1] + real :: rho_scale ! A scaling factor for densities from kg m-3 to R [R m3 kg-1 ~> 1] real :: rho_ref_mks ! The reference density in MKS units, never rescaled from kg m-3 [kg m-3] real :: dz real :: weight_t, weight_b @@ -1881,8 +2298,9 @@ subroutine int_density_dz_generic_ppm (T, T_t, T_b, S, S_t, S_b, & is = HIO%isc + ioff ; ie = HIO%iec + ioff js = HIO%jsc + joff ; je = HIO%jec + joff - GxRho = G_e * rho_0 ; if (present(pres_scale)) GxRho = pres_scale * G_e * rho_0 - rho_ref_mks = rho_ref ; if (present(rho_scale)) rho_ref_mks = rho_ref / rho_scale + rho_scale = 1.0 ; if (present(US)) rho_scale = US%kg_m3_to_R + GxRho = G_e * rho_0 ; if (present(US)) GxRho = US%RL2_T2_to_Pa * G_e * rho_0 + rho_ref_mks = rho_ref ; if (present(US)) rho_ref_mks = rho_ref * US%R_to_kg_m3 I_Rho = 1.0 / rho_0 ! ============================= @@ -1910,7 +2328,11 @@ subroutine int_density_dz_generic_ppm (T, T_t, T_b, S, S_t, S_b, & T5(n) = t0 + t1 * xi + t2 * xi**2 enddo - call calculate_density(T5, S5, p5, r5, 1, 5, EOS, rho_ref_mks, scale=rho_scale) + if (rho_scale /= 1.0) then + call calculate_density(T5, S5, p5, r5, 1, 5, EOS, rho_ref_mks, scale=rho_scale) + else + call calculate_density(T5, S5, p5, r5, 1, 5, EOS, rho_ref_mks) + endif ! Use Bode's rule to estimate the pressure anomaly change. rho_anom = C1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + 12.0*r5(3)) @@ -1969,7 +2391,11 @@ subroutine int_density_dz_generic_ppm (T, T_t, T_b, S, S_t, S_b, & T5(n) = t0 + t1 * xi + t2 * xi**2 enddo - call calculate_density(T5, S5, p5, r5, 1, 5, EOS, rho_ref_mks, scale=rho_scale) + if (rho_scale /= 1.0) then + call calculate_density(T5, S5, p5, r5, 1, 5, EOS, rho_ref_mks, scale=rho_scale) + else + call calculate_density(T5, S5, p5, r5, 1, 5, EOS, rho_ref_mks) + endif ! Use Bode's rule to estimate the pressure anomaly change. intz(m) = G_e*dz*( C1_90*(7.0*(r5(1)+r5(5)) + 32.0*(r5(2)+r5(4)) + & @@ -2016,7 +2442,11 @@ subroutine int_density_dz_generic_ppm (T, T_t, T_b, S, S_t, S_b, & S_node(9) = 0.5 * ( S_node(6) + S_node(8) ) S_node(7) = 0.5 * ( S_node(3) + S_node(4) ) - call calculate_density( T_node, S_node, p_node, r_node, 1, 9, EOS, rho_ref_mks, scale=rho_scale ) + if (rho_scale /= 1.0) then + call calculate_density( T_node, S_node, p_node, r_node, 1, 9, EOS, rho_ref_mks, scale=rho_scale ) + else + call calculate_density( T_node, S_node, p_node, r_node, 1, 9, EOS, rho_ref_mks) + endif r_node = r_node - rho_ref call compute_integral_quadratic( x, y, r_node, intx_dpa(i-ioff,j-joff) ) @@ -2231,9 +2661,9 @@ end subroutine evaluate_shape_quadratic !! pressure across layers, which are required for calculating the finite-volume !! form pressure accelerations in a non-Boussinesq model. There are essentially !! no free assumptions, apart from the use of Bode's rule quadrature to do the integrals. -subroutine int_spec_vol_dp_generic(T, S, p_t, p_b, alpha_ref, HI, EOS, & - dza, intp_dza, intx_dza, inty_dza, halo_size, & - bathyP, dP_neglect, useMassWghtInterp, SV_scale, pres_scale) +subroutine int_spec_vol_dp_generic(T, S, p_t, p_b, alpha_ref, HI, EOS, dza, & + intp_dza, intx_dza, inty_dza, halo_size, & + bathyP, dP_neglect, useMassWghtInterp, US) type(hor_index_type), intent(in) :: HI !< A horizontal index type structure. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & intent(in) :: T !< Potential temperature of the layer [degC]. @@ -2271,10 +2701,7 @@ subroutine int_spec_vol_dp_generic(T, S, p_t, p_b, alpha_ref, HI, EOS, & !! the same units as p_t [R L2 T-2 ~> Pa] or [Pa] logical, optional, intent(in) :: useMassWghtInterp !< If true, uses mass weighting !! to interpolate T/S for top and bottom integrals. - real, optional, intent(in) :: SV_scale !< A multiplicative factor by which to scale specific - !! volume from m3 kg-1 to the desired units [kg m-3 R-1 ~> 1] - real, optional, intent(in) :: pres_scale !< A multiplicative factor to convert pressure - !! into Pa [Pa T2 R-1 L-2 ~> 1]. + type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type ! This subroutine calculates analytical and nearly-analytical integrals in ! pressure across layers of geopotential anomalies, which are required for @@ -2301,6 +2728,8 @@ subroutine int_spec_vol_dp_generic(T, S, p_t, p_b, alpha_ref, HI, EOS, & real :: intp(5) ! The integrals of specific volume with pressure at the ! 5 sub-column locations [L2 T-2 ~> m2 s-2]. real :: RL2_T2_to_Pa ! A unit conversion factor from the rescaled units of pressure to Pa [Pa T2 R-1 L-2 ~> 1] + real :: SV_scale ! A multiplicative factor by which to scale specific + ! volume from m3 kg-1 to the desired units [kg m-3 R-1 ~> 1] logical :: do_massWeight ! Indicates whether to do mass weighting. real, parameter :: C1_90 = 1.0/90.0 ! A rational constant. integer :: Isq, Ieq, Jsq, Jeq, ish, ieh, jsh, jeh, i, j, m, n, halo @@ -2311,8 +2740,9 @@ subroutine int_spec_vol_dp_generic(T, S, p_t, p_b, alpha_ref, HI, EOS, & if (present(intx_dza)) then ; ish = MIN(Isq,ish) ; ieh = MAX(Ieq+1,ieh); endif if (present(inty_dza)) then ; jsh = MIN(Jsq,jsh) ; jeh = MAX(Jeq+1,jeh); endif - RL2_T2_to_Pa = 1.0 ; if (present(pres_scale)) RL2_T2_to_Pa = pres_scale - alpha_ref_mks = alpha_ref ; if (present(SV_scale)) alpha_ref_mks = alpha_ref / SV_scale + SV_scale = 1.0 ; if (present(US)) SV_scale = US%R_to_kg_m3 + RL2_T2_to_Pa = 1.0 ; if (present(US)) RL2_T2_to_Pa = US%RL2_T2_to_Pa + alpha_ref_mks = alpha_ref ; if (present(US)) alpha_ref_mks = alpha_ref * US%kg_m3_to_R do_massWeight = .false. if (present(useMassWghtInterp)) then ; if (useMassWghtInterp) then @@ -2329,7 +2759,12 @@ subroutine int_spec_vol_dp_generic(T, S, p_t, p_b, alpha_ref, HI, EOS, & T5(n) = T(i,j) ; S5(n) = S(i,j) p5(n) = RL2_T2_to_Pa * (p_b(i,j) - 0.25*real(n-1)*dp) enddo - call calculate_spec_vol(T5, S5, p5, a5, 1, 5, EOS, alpha_ref_mks, scale=SV_scale) + + if (SV_scale /= 1.0) then + call calculate_spec_vol(T5, S5, p5, a5, 1, 5, EOS, alpha_ref_mks, scale=SV_scale) + else + call calculate_spec_vol(T5, S5, p5, a5, 1, 5, EOS, alpha_ref_mks) + endif ! Use Bode's rule to estimate the interface height anomaly change. alpha_anom = C1_90*(7.0*(a5(1)+a5(5)) + 32.0*(a5(2)+a5(4)) + 12.0*a5(3)) @@ -2373,7 +2808,11 @@ subroutine int_spec_vol_dp_generic(T, S, p_t, p_b, alpha_ref, HI, EOS, & do n=2,5 T5(n) = T5(1) ; S5(n) = S5(1) ; p5(n) = p5(n-1) - RL2_T2_to_Pa * 0.25*dp enddo - call calculate_spec_vol(T5, S5, p5, a5, 1, 5, EOS, alpha_ref_mks, scale=SV_scale) + if (SV_scale /= 1.0) then + call calculate_spec_vol(T5, S5, p5, a5, 1, 5, EOS, alpha_ref_mks, scale=SV_scale) + else + call calculate_spec_vol(T5, S5, p5, a5, 1, 5, EOS, alpha_ref_mks) + endif ! Use Bode's rule to estimate the interface height anomaly change. intp(m) = dp*( C1_90*(7.0*(a5(1)+a5(5)) + 32.0*(a5(2)+a5(4)) + & @@ -2416,7 +2855,11 @@ subroutine int_spec_vol_dp_generic(T, S, p_t, p_b, alpha_ref, HI, EOS, & do n=2,5 T5(n) = T5(1) ; S5(n) = S5(1) ; p5(n) = RL2_T2_to_Pa * (p5(n-1) - 0.25*dp) enddo - call calculate_spec_vol(T5, S5, p5, a5, 1, 5, EOS, alpha_ref_mks, scale=SV_scale) + if (SV_scale /= 1.0) then + call calculate_spec_vol(T5, S5, p5, a5, 1, 5, EOS, alpha_ref_mks, scale=SV_scale) + else + call calculate_spec_vol(T5, S5, p5, a5, 1, 5, EOS, alpha_ref_mks) + endif ! Use Bode's rule to estimate the interface height anomaly change. intp(m) = dp*( C1_90*(7.0*(a5(1)+a5(5)) + 32.0*(a5(2)+a5(4)) + & @@ -2435,7 +2878,7 @@ end subroutine int_spec_vol_dp_generic !! no free assumptions, apart from the use of Bode's rule quadrature to do the integrals. subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, & dP_neglect, bathyP, HI, EOS, dza, & - intp_dza, intx_dza, inty_dza, useMassWghtInterp, SV_scale, pres_scale) + intp_dza, intx_dza, inty_dza, useMassWghtInterp, US) type(hor_index_type), intent(in) :: HI !< A horizontal index type structure. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & intent(in) :: T_t !< Potential temperature at the top of the layer [degC]. @@ -2476,10 +2919,7 @@ subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, !! by the y grid spacing [L2 T-2 ~> m2 s-2] or [m2 s-2]. logical, optional, intent(in) :: useMassWghtInterp !< If true, uses mass weighting !! to interpolate T/S for top and bottom integrals. - real, optional, intent(in) :: SV_scale !< A multiplicative factor by which to scale specific - !! volume from m3 kg-1 to the desired units [kg m-3 R-1 ~> 1] - real, optional, intent(in) :: pres_scale !< A multiplicative factor to convert pressure - !! into Pa [Pa T2 R-1 L-2 ~> 1]. + type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type ! This subroutine calculates analytical and nearly-analytical integrals in ! pressure across layers of geopotential anomalies, which are required for @@ -2513,6 +2953,8 @@ subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, real :: intp(5) ! The integrals of specific volume with pressure at the ! 5 sub-column locations [L2 T-2 ~> m2 s-2]. real :: RL2_T2_to_Pa ! A unit conversion factor from the rescaled units of pressure to Pa [Pa T2 R-1 L-2 ~> 1] + real :: SV_scale ! A multiplicative factor by which to scale specific + ! volume from m3 kg-1 to the desired units [kg m-3 R-1 ~> 1] real, parameter :: C1_90 = 1.0/90.0 ! A rational constant. logical :: do_massWeight ! Indicates whether to do mass weighting. integer :: Isq, Ieq, Jsq, Jeq, i, j, m, n, pos @@ -2522,9 +2964,9 @@ subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, do_massWeight = .false. if (present(useMassWghtInterp)) do_massWeight = useMassWghtInterp - RL2_T2_to_Pa = 1.0 ; if (present(pres_scale)) RL2_T2_to_Pa = pres_scale - alpha_ref_mks = alpha_ref ; if (present(SV_scale)) alpha_ref_mks = alpha_ref / SV_scale - + SV_scale = 1.0 ; if (present(US)) SV_scale = US%R_to_kg_m3 + RL2_T2_to_Pa = 1.0 ; if (present(US)) RL2_T2_to_Pa = US%RL2_T2_to_Pa + alpha_ref_mks = alpha_ref ; if (present(US)) alpha_ref_mks = alpha_ref * US%kg_m3_to_R do n = 1, 5 ! Note that these are reversed from int_density_dz. wt_t(n) = 0.25 * real(n-1) @@ -2541,7 +2983,11 @@ subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, S5(n) = wt_t(n) * S_t(i,j) + wt_b(n) * S_b(i,j) T5(n) = wt_t(n) * T_t(i,j) + wt_b(n) * T_b(i,j) enddo - call calculate_spec_vol(T5, S5, p5, a5, 1, 5, EOS, alpha_ref_mks, scale=SV_scale) + if (SV_scale /= 1.0) then + call calculate_spec_vol(T5, S5, p5, a5, 1, 5, EOS, alpha_ref_mks, scale=SV_scale) + else + call calculate_spec_vol(T5, S5, p5, a5, 1, 5, EOS, alpha_ref_mks) + endif ! Use Bode's rule to estimate the interface height anomaly change. alpha_anom = C1_90*((7.0*(a5(1)+a5(5)) + 32.0*(a5(2)+a5(4))) + 12.0*a5(3)) @@ -2598,7 +3044,11 @@ subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, enddo enddo - call calculate_spec_vol(T15, S15, p15, a15, 1, 15, EOS, alpha_ref_mks, scale=SV_scale) + if (SV_scale /= 1.0) then + call calculate_spec_vol(T15, S15, p15, a15, 1, 15, EOS, alpha_ref_mks, scale=SV_scale) + else + call calculate_spec_vol(T15, S15, p15, a15, 1, 15, EOS, alpha_ref_mks) + endif intp(1) = dza(i,j) ; intp(5) = dza(i+1,j) do m=2,4 @@ -2657,7 +3107,11 @@ subroutine int_spec_vol_dp_generic_plm(T_t, T_b, S_t, S_b, p_t, p_b, alpha_ref, enddo enddo - call calculate_spec_vol(T15, S15, p15, a15, 1, 15, EOS, alpha_ref_mks, scale=SV_scale) + if (SV_scale /= 1.0) then + call calculate_spec_vol(T15, S15, p15, a15, 1, 15, EOS, alpha_ref_mks, scale=SV_scale) + else + call calculate_spec_vol(T15, S15, p15, a15, 1, 15, EOS, alpha_ref_mks) + endif intp(1) = dza(i,j) ; intp(5) = dza(i,j+1) do m=2,4