Skip to content

Commit

Permalink
+Eliminated the calculate_2_densities routines
Browse files Browse the repository at this point in the history 
  Eliminated the routine calculate_2_densities and each of the EOS-specific
variants  that MOM6 used to calculate densities at two different reference
pressures at the same time.  These routines were no longer being used and this
functionality can easily be recovered with two calls to calculate_density.
All answers are bitwise identical.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed Jan 12, 2017
1 parent 4952f3d commit 5147edc
Show file tree
Hide file tree
Showing 11 changed files with 13 additions and 273 deletions.
47 changes: 6 additions & 41 deletions src/equation_of_state/MOM_EOS.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,19 +5,19 @@ module MOM_EOS

use MOM_EOS_linear, only : calculate_density_scalar_linear, calculate_density_array_linear
use MOM_EOS_linear, only : calculate_density_derivs_linear, calculate_specvol_derivs_linear, int_density_dz_linear
use MOM_EOS_linear, only : calculate_compress_linear, calculate_2_densities_linear, int_spec_vol_dp_linear
use MOM_EOS_linear, only : calculate_compress_linear, int_spec_vol_dp_linear
use MOM_EOS_Wright, only : calculate_density_scalar_wright, calculate_density_array_wright
use MOM_EOS_Wright, only : calculate_density_derivs_wright, calculate_specvol_derivs_wright, int_density_dz_wright
use MOM_EOS_Wright, only : calculate_compress_wright, calculate_2_densities_wright, int_spec_vol_dp_wright
use MOM_EOS_Wright, only : calculate_compress_wright, int_spec_vol_dp_wright
use MOM_EOS_UNESCO, only : calculate_density_scalar_unesco, calculate_density_array_unesco
use MOM_EOS_UNESCO, only : calculate_density_derivs_unesco, calculate_density_unesco
use MOM_EOS_UNESCO, only : calculate_compress_unesco, calculate_2_densities_unesco
use MOM_EOS_UNESCO, only : calculate_compress_unesco
use MOM_EOS_NEMO, only : calculate_density_scalar_nemo, calculate_density_array_nemo
use MOM_EOS_NEMO, only : calculate_density_derivs_nemo, calculate_density_nemo
use MOM_EOS_NEMO, only : calculate_compress_nemo, calculate_2_densities_nemo
use MOM_EOS_NEMO, only : calculate_compress_nemo
use MOM_EOS_TEOS10, only : calculate_density_scalar_teos10, calculate_density_array_teos10
use MOM_EOS_TEOS10, only : calculate_density_derivs_teos10, calculate_specvol_derivs_teos10
use MOM_EOS_TEOS10, only : calculate_compress_teos10, calculate_2_densities_teos10
use MOM_EOS_TEOS10, only : calculate_compress_teos10
use MOM_EOS_TEOS10, only : gsw_sp_from_sr, gsw_pt_from_ct
use MOM_TFreeze, only : calculate_TFreeze_linear, calculate_TFreeze_Millero
use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg
Expand All @@ -30,8 +30,7 @@ module MOM_EOS
#include <MOM_memory.h>

public calculate_compress, calculate_density, query_compressible
public calculate_density_derivs, calculate_2_densities
public calculate_specific_vol_derivs
public calculate_density_derivs, calculate_specific_vol_derivs
public EOS_init, EOS_end, EOS_allocate
public EOS_use_linear
public int_density_dz, int_specific_vol_dp
Expand Down Expand Up @@ -317,40 +316,6 @@ subroutine calculate_compress(T, S, pressure, rho, drho_dp, start, npts, EOS)

end subroutine calculate_compress

!> Calls the appropriate subroutine to calculate density for two arrays.
subroutine calculate_2_densities( T, S, pressure1, pressure2, rho1, rho2, start, npts, EOS)
real, dimension(:), intent(in) :: T !< Potential temperature referenced to the surface (degC)
real, dimension(:), intent(in) :: S !< Salinity (PSU)
real, intent(in) :: pressure1 !< Pressure (Pa)
real, intent(in) :: pressure2 !< A second pressure (Pa)
real, dimension(:), intent(out) :: rho1 !< Density at pressure1, in kg m-3.
real, dimension(:), intent(out) :: rho2 !< Density at pressure2, in kg m-3.
integer, intent(in) :: start !< Starting index within the array
integer, intent(in) :: npts !< The number of values to calculate
type(EOS_type), pointer :: EOS !< Equation of state structure

if (.not.associated(EOS)) call MOM_error(FATAL, &
"calculate_2_densities called with an unassociated EOS_type EOS.")

select case (EOS%form_of_EOS)
case (EOS_LINEAR)
call calculate_2_densities_linear(T, S, pressure1, pressure2, rho1, rho2, start, &
npts, EOS%Rho_T0_S0, EOS%dRho_dT, EOS%dRho_dS)
case (EOS_UNESCO)
call calculate_2_densities_unesco(T, S, pressure1, pressure2, rho1, rho2, start, npts)
case (EOS_WRIGHT)
call calculate_2_densities_wright(T, S, pressure1, pressure2, rho1, rho2, start, npts)
case (EOS_TEOS10)
call calculate_2_densities_teos10(T, S, pressure1, pressure2, rho1, rho2, start, npts)
case (EOS_NEMO)
call calculate_2_densities_nemo(T, S, pressure1, pressure2, rho1, rho2, start, npts)
case default
call MOM_error(FATAL, &
"calculate_2_densities: EOS%form_of_EOS is not valid.")
end select

end subroutine calculate_2_densities

!> Calls the appropriate subroutine to alculate analytical and nearly-analytical
!! integrals in pressure across layers of geopotential anomalies, which are
!! required for calculating the finite-volume form pressure accelerations in a
Expand Down
71 changes: 1 addition & 70 deletions src/equation_of_state/MOM_EOS_NEMO.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,11 +32,10 @@ module MOM_EOS_NEMO
!use gsw_mod_toolbox, only : gsw_sr_from_sp, gsw_ct_from_pt
use gsw_mod_toolbox, only : gsw_rho_first_derivatives


implicit none ; private

public calculate_compress_nemo, calculate_density_nemo
public calculate_density_derivs_nemo, calculate_2_densities_nemo
public calculate_density_derivs_nemo
public calculate_density_scalar_nemo, calculate_density_array_nemo

interface calculate_density_nemo
Expand Down Expand Up @@ -374,72 +373,4 @@ subroutine calculate_compress_nemo(T, S, pressure, rho, drho_dp, start, npts)
enddo
end subroutine calculate_compress_nemo

subroutine calculate_2_densities_nemo( T, S, pressure1, pressure2, rho1, rho2, start, npts)
real, intent(in), dimension(:) :: T, S
real, intent(in) :: pressure1, pressure2
real, intent(out), dimension(:) :: rho1, rho2
integer, intent(in) :: start, npts
! * Arguments: T - conservative temperature in C. *
! * (in) S - absolute salinity in g/kg. *
! * (in) pressure1 - the first pressure in Pa. *
! * (in) pressure2 - the second pressure in Pa. *
! * (out) rho1 - density at pressure1 in kg m-3. *
! * (out) rho2 - density at pressure2 in kg m-3. *
! * (in) start - the starting point in the arrays. *
! * (in) npts - the number of values to calculate. *

real :: zp1, zp2, zt , zh , zs , zr0, zn , zn0, zn1, zn2, zn3
integer :: j

zp1 = pressure1 * Pa2db !Convert pressure from Pascal to decibar
zp2 = pressure2 * Pa2db !Convert pressure from Pascal to decibar

do j=start,start+npts-1
!Conversions
zs = S(j) !gsw_sr_from_sp(S(j)) !Convert practical salinity to absolute salinity
zt = T(j) !gsw_ct_from_pt(S(j),T(j)) !Convert potantial temp to conservative temp

!The following algorithm was provided by Roquet in a private communication.
!It is not necessarily the algorithm used in NEMO ocean!
zp1 = pressure1 * r1_P0 ! pressure (first converted to decibar)
zp2 = pressure2 * r1_P0 ! pressure (first converted to decibar)
zt = zt * r1_T0 ! temperature
zs = SQRT( ABS( zs + rdeltaS ) * r1_S0 ) ! square root salinity
!
zn3 = EOS013*zt &
& + EOS103*zs+EOS003
!
zn2 = (EOS022*zt &
& + EOS112*zs+EOS012)*zt &
& + (EOS202*zs+EOS102)*zs+EOS002
!
zn1 = (((EOS041*zt &
& + EOS131*zs+EOS031)*zt &
& + (EOS221*zs+EOS121)*zs+EOS021)*zt &
& + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt &
& + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001
!
zn0 = (((((EOS060*zt &
& + EOS150*zs+EOS050)*zt &
& + (EOS240*zs+EOS140)*zs+EOS040)*zt &
& + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt &
& + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt &
& + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt &
& + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
!
zn = ( ( zn3 * zp1 + zn2 ) * zp1 + zn1 ) * zp1 + zn0
!
zr0 = (((((R05 * zp1+R04) * zp1+R03 ) * zp1+R02 ) * zp1+R01) * zp1+R00) * zp1
!
rho1(j) = ( zn + zr0 ) ! density
!
zn = ( ( zn3 * zp2 + zn2 ) * zp2 + zn1 ) * zp2 + zn0
!
zr0 = (((((R05 * zp2+R04) * zp2+R03 ) * zp2+R02 ) * zp2+R01) * zp2+R00) * zp2
!
rho2(j) = ( zn + zr0 ) ! density
enddo
end subroutine calculate_2_densities_nemo


end module MOM_EOS_NEMO
36 changes: 2 additions & 34 deletions src/equation_of_state/MOM_EOS_TEOS10.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,16 +32,15 @@ module MOM_EOS_TEOS10
implicit none ; private

public calculate_compress_teos10, calculate_density_teos10
public calculate_density_derivs_teos10, calculate_2_densities_teos10
public calculate_specvol_derivs_teos10
public calculate_density_derivs_teos10, calculate_specvol_derivs_teos10
public calculate_density_scalar_teos10, calculate_density_array_teos10
public gsw_sp_from_sr, gsw_pt_from_ct

interface calculate_density_teos10
module procedure calculate_density_scalar_teos10, calculate_density_array_teos10
end interface calculate_density_teos10

real, parameter :: Pa2db = 1.e-4
real, parameter :: Pa2db = 1.e-4 ! The conversion factor from Pa to dbar.

contains

Expand Down Expand Up @@ -195,35 +194,4 @@ subroutine calculate_compress_teos10(T, S, pressure, rho, drho_dp, start, npts)
enddo
end subroutine calculate_compress_teos10

subroutine calculate_2_densities_teos10( T, S, pressure1, pressure2, rho1, rho2, start, npts)
real, intent(in), dimension(:) :: T, S
real, intent(in) :: pressure1, pressure2
real, intent(out), dimension(:) :: rho1, rho2
integer, intent(in) :: start, npts
! * Arguments: T - conservative temperature in C. *
! * (in) S - absolute salinity in g/kg. *
! * (in) pressure1 - the first pressure in Pa. *
! * (in) pressure2 - the second pressure in Pa. *
! * (out) rho1 - density at pressure1 in kg m-3. *
! * (out) rho2 - density at pressure2 in kg m-3. *
! * (in) start - the starting point in the arrays. *
! * (in) npts - the number of values to calculate. *

real :: zp1, zp2, zt, zs
integer :: j

zp1 = pressure1 * Pa2db !Convert pressure from Pascal to decibar
zp2 = pressure2 * Pa2db !Convert pressure from Pascal to decibar

do j=start,start+npts-1
!Conversions
zs = S(j) !gsw_sr_from_sp(S(j)) !Convert practical salinity to absolute salinity
zt = T(j) !gsw_ct_from_pt(S(j),T(j)) !Convert potantial temp to conservative temp
if(S(j).lt.-1.0e-10) cycle !Can we assume safely that this is a missing value?
rho1(j) = gsw_rho(zs,zt,zp1)
rho2(j) = gsw_rho(zs,zt,zp2)
enddo
end subroutine calculate_2_densities_teos10


end module MOM_EOS_TEOS10
56 changes: 1 addition & 55 deletions src/equation_of_state/MOM_EOS_UNESCO.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@ module MOM_EOS_UNESCO
implicit none ; private

public calculate_compress_UNESCO, calculate_density_UNESCO
public calculate_density_derivs_UNESCO, calculate_2_densities_UNESCO
public calculate_density_derivs_UNESCO
public calculate_density_scalar_UNESCO, calculate_density_array_UNESCO

interface calculate_density_UNESCO
Expand Down Expand Up @@ -255,59 +255,5 @@ subroutine calculate_compress_UNESCO(T, S, pressure, rho, drho_dp, start, npts)
enddo
end subroutine calculate_compress_UNESCO

subroutine calculate_2_densities_UNESCO( T, S, pressure1, pressure2, rho1, rho2, start, npts)
real, intent(in), dimension(:) :: T, S
real, intent(in) :: pressure1, pressure2
real, intent(out), dimension(:) :: rho1, rho2
integer, intent(in) :: start, npts
! * This subroutine computes the densities of sea water (rho1 and *
! * rho2) at two reference pressures (pressure1 and pressure2) from *
! * salinity and potential temperature. *
! * *
! * Arguments: T - potential temperature relative to the surface in C. *
! * (in) S - salinity in PSU. *
! * (in) pressure1 - the first pressure in Pa. *
! * (in) pressure2 - the second pressure in Pa. *
! * (out) rho1 - density at pressure1 in kg m-3. *
! * (out) rho2 - density at pressure2 in kg m-3. *
! * (in) start - the starting point in the arrays. *
! * (in) npts - the number of values to calculate. *
real :: t_local, t2, t3, t4, t5; ! Temperature to the 1st - 5th power.
real :: s_local, s32, s2; ! Salinity to the 1st, 3/2, & 2nd power.
real :: p1a, p2a; ! Pressure1 (in bars) to the 1st and 2nd power.
real :: p1b, p2b; ! Pressure2 (in bars) to the 1st and 2nd power.
real :: rho0; ! Density at 1 bar pressure, in kg m-3.
real :: ksa, ksb; ! The secant bulk modulus in bar.
real :: ks_0, ks_1, ks_2;
integer :: j

p1a = pressure1*1.0e-5; p2a = p1a*p1a;
p1b = pressure2*1.0e-5; p2b = p1b*p1b;

do j=start, start+npts-1
t_local = T(j); t2 = t_local*t_local; t3 = t_local*t2; t4 = t2*t2; t5 = t3*t2;
s_local = S(j); s2 = s_local*s_local; s32 = s_local*sqrt(s_local);

! Compute rho(s,theta,p=0) - (same as rho(s,t_insitu,p=0) ).

rho0 = R00 + R10*t_local + R20*t2 + R30*t3 + R40*t4 + R50*t5 + &
s_local*(R01 + R11*t_local + R21*t2 + R31*t3 + R41*t4) + &
s32*(R032 + R132*t_local + R232*t2) + R02*s2;

! Compute rho(s,theta,p), first calculating the secant bulk modulus.

ks_0 = S00 + S10*t_local + S20*t2 + S30*t3 + S40*t4 + &
s_local*(S01 + S11*t_local + S21*t2 + S31*t3) + s32*(S032 + S132*t_local + S232*t2);
ks_1 = Sp00 + Sp10*t_local + Sp20*t2 + Sp30*t3 + &
s_local*(Sp01 + Sp11*t_local + Sp21*t2) + Sp032*s32;
ks_2 = SP000 + SP010*t_local + SP020*t2 + s_local*(SP001 + SP011*t_local + SP021*t2);

ksa = ks_0 + p1a*ks_1 + p2a*ks_2;
ksb = ks_0 + p1b*ks_1 + p2b*ks_2;

rho1(j) = rho0*ksa / (ksa - p1a);
rho2(j) = rho0*ksb / (ksb - p1b);
enddo
end subroutine calculate_2_densities_UNESCO

end module MOM_EOS_UNESCO
37 changes: 1 addition & 36 deletions src/equation_of_state/MOM_EOS_Wright.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,8 +32,7 @@ module MOM_EOS_Wright
#include <MOM_memory.h>

public calculate_compress_wright, calculate_density_wright
public calculate_density_derivs_wright, calculate_2_densities_wright
public calculate_specvol_derivs_wright
public calculate_density_derivs_wright, calculate_specvol_derivs_wright
public calculate_density_scalar_wright, calculate_density_array_wright
public int_density_dz_wright, int_spec_vol_dp_wright

Expand Down Expand Up @@ -223,40 +222,6 @@ subroutine calculate_compress_wright(T, S, pressure, rho, drho_dp, start, npts)
enddo
end subroutine calculate_compress_wright

subroutine calculate_2_densities_wright( T, S, pressure1, pressure2, rho1, rho2, start, npts)
real, intent(in), dimension(:) :: T, S
real, intent(in) :: pressure1, pressure2
real, intent(out), dimension(:) :: rho1, rho2
integer, intent(in) :: start, npts
! * Arguments: T - potential temperature relative to the surface in C. *
! * (in) S - salinity in PSU. *
! * (in) pressure1 - the first pressure in Pa. *
! * (in) pressure2 - the second pressure in Pa. *
! * (out) rho1 - density at pressure1 in kg m-3. *
! * (out) rho2 - density at pressure2 in kg m-3. *
! * (in) start - the starting point in the arrays. *
! * (in) npts - the number of values to calculate. *

! *====================================================================*
! * This subroutine computes the in situ density of sea water (rho in *
! * units of kg/m^3) from salinity (sal in psu), potential temperature*
! * (th in deg C), and pressure in Pa. It uses the expression from *
! * Wright, 1997, J. Atmos. Ocean. Tech., 14, 735-740. *
! * Coded by R. Hallberg, 7/00 *
! *====================================================================*
real :: al0, p0, lambda
integer :: j

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))

rho1(j) = (pressure1 + p0) / (lambda + al0*(pressure1 + p0))
rho2(j) = (pressure2 + p0) / (lambda + al0*(pressure2 + p0))
enddo
end subroutine calculate_2_densities_wright

subroutine 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)
type(hor_index_type), intent(in) :: HII, HIO
Expand Down
33 changes: 1 addition & 32 deletions src/equation_of_state/MOM_EOS_linear.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,8 +32,7 @@ module MOM_EOS_linear
#include <MOM_memory.h>

public calculate_compress_linear, calculate_density_linear
public calculate_density_derivs_linear, calculate_2_densities_linear
public calculate_specvol_derivs_linear
public calculate_density_derivs_linear, calculate_specvol_derivs_linear
public calculate_density_scalar_linear, calculate_density_array_linear
public int_density_dz_linear, int_spec_vol_dp_linear

Expand Down Expand Up @@ -181,36 +180,6 @@ subroutine calculate_compress_linear(T, S, pressure, rho, drho_dp, start, npts,&
enddo
end subroutine calculate_compress_linear

subroutine calculate_2_densities_linear(T, S, pressure1, pressure2, rho1, rho2,&
start, npts, Rho_T0_S0, dRho_dT, dRho_dS)
real, intent(in), dimension(:) :: T, S
real, intent(in) :: pressure1, pressure2
real, intent(out), dimension(:) :: rho1, rho2
integer, intent(in) :: start, npts
real, intent(in) :: Rho_T0_S0, dRho_dT, dRho_dS
! * This subroutine computes the densities of sea water (rho1 and *
! * rho2) at two reference pressures (pressure1 and pressure2) from *
! * salinity and potential temperature. *
! * *
! * Arguments: T - potential temperature relative to the surface in C. *
! * (in) S - salinity in PSU. *
! * (in) pressure1 - the first pressure in Pa. *
! * (in) pressure2 - the second pressure in Pa. *
! * (out) rho1 - density at pressure1 in kg m-3. *
! * (out) rho2 - density at pressure2 in kg m-3. *
! * (in) start - the starting point in the arrays. *
! * (in) npts - the number of values to calculate. *
! * (in) Rho_T0_S0 - The density at T=0, S=0, in kg m-3. *
! * (in) dRho_dT - The derivatives of density with temperature *
! * (in) dRho_dS - and salinity, in kg m-3 C-1 and kg m-3 psu-1. *
integer :: j

do j=start, start+npts-1
rho1(j) = Rho_T0_S0 + dRho_dT*T(j) + dRho_dS*S(j);
rho2(j) = Rho_T0_S0 + dRho_dT*T(j) + dRho_dS*S(j);
enddo
end subroutine calculate_2_densities_linear

subroutine int_density_dz_linear(T, S, z_t, z_b, rho_ref, rho_0_pres, G_e, HII, HIO, &
Rho_T0_S0, dRho_dT, dRho_dS, dpa, intz_dpa, intx_dpa, inty_dpa)
type(hor_index_type), intent(in) :: HII, HIO
Expand Down
1 change: 0 additions & 1 deletion src/parameterizations/vertical/MOM_bulk_mixed_layer.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,7 +65,6 @@ module MOM_bulk_mixed_layer
use MOM_variables, only : thermo_var_ptrs
use MOM_verticalGrid, only : verticalGrid_type
use MOM_EOS, only : calculate_density, calculate_density_derivs
use MOM_EOS, only : calculate_2_densities

implicit none ; private

Expand Down
Loading

0 comments on commit 5147edc

Please sign in to comment.