!/===========================================================================/
! Copyright (c) 2007, The University of Massachusetts Dartmouth
! Produced at the School of Marine Science & Technology
! Marine Ecosystem Dynamics Modeling group
! All rights reserved.
!
! FVCOM has been developed by the joint UMASSD-WHOI research team. For
! details of authorship and attribution of credit please see the FVCOM
! technical manual or contact the MEDM group.
!
!
! This file is part of FVCOM. For details, see http://fvcom.smast.umassd.edu
! The full copyright notice is contained in the file COPYRIGHT located in the
! root directory of the FVCOM code. This original header must be maintained
! in all distributed versions.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
! THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED.
!
!/---------------------------------------------------------------------------/
! CVS VERSION INFORMATION
! $Id$
! $Name$
! $Revision$
!/===========================================================================/
MODULE EQS_OF_STATE
! THIS MODULE CONTAINS THREE EQUATIONS OF STATE WHICH CAN BE ACCESSED
! USING A NUMBER OF DIFFERENT INTERFACES. THE ORIGINAL FVCOM
! INTERFACE (DENS1,DENS2,DENS3) STILL EXISTS, WHERE BOTH RHO AND
! RHO1 ARE UPDATED, OR THE USER CAN ACCESS THE UNDERLIEING
! FUNCTIONS AND PASS AN ARRAY, 1D OR 2D TO THE SUBROUTINE.
!
! SEE DETAILS OF METHOD IN THE HEADER FOR EACH SUBROUTINE
!
! SUBROUTINES PUBLICLY AVAILABLE IN THIS MODULE:
!=================================================================================
! DENS1 - INTERFACE: CALL DENS1
! DENS2 - INTERFACE: CALL DENS2
! DENS3 - INTERFACE: CALL DENS3
! - THESE ROUTINES USE THE VALUES IN S1,T1,ZZ TO UPDATE THE DENSITY IN RHO AND RHO1
!==================================================================================
! FOFONOFF_MILLARD - INTERFACE: CALL FOFONOFF_MILLARD_2D(S,T,Z,PREF,RHO)
! REAL(SP), INTENT(IN),DIMENSION(:) :: S,T,P ! SALINITY,TEMPERATURE, PRESSURE
! REAL(SP), INTENT(OUT),DIMENSION(:) :: RHO ! RESULT - DENSITY
! REAL(SP), INTENT(IN) :: PREF ! REFERENCE PRESSURE
! or
! REAL(SP), INTENT(IN),DIMENSION(:,:) :: S,T,P ! SALINITY,TEMPERATURE, PRESSURE
! REAL(SP), INTENT(OUT),DIMENSION(:,:) :: RHO ! RESULT - DENSITY
! REAL(SP), INTENT(IN) :: PREF ! REFERENCE PRESSURE
!=================================================================================
!
!=================================================================================
! DENS2G - INTERFACE: CALL DENS2G(S,T,RHO)
! REAL(SP), INTENT(IN),DIMENSION(:) :: S,T ! SALINITY,TEMPERATURE
! REAL(SP), INTENT(OUT),DIMENSION(:) :: RHO ! DENSITY
! or
! REAL(SP), INTENT(IN),DIMENSION(:,:) :: S,T,P ! SALINITY,TEMPERATURE
! REAL(SP), INTENT(OUT),DIMENSION(:,:) :: RHO ! DENSITY
!=================================================================================
!
!=================================================================================
! JACKET_MCDOUGALL - INTERFACE: CALL JACKET_MCDOUGALL(S,T,Z,RHO)
! REAL(SP), INTENT(IN),DIMENSION(:) :: S,T,P ! SALINITY,TEMPERATURE, PRESSURE
! REAL(SP), INTENT(OUT),DIMENSION(:) :: RHO ! RESULT - DENSITY
! or
! REAL(SP), INTENT(IN),DIMENSION(:,:) :: S,T,P ! SALINITY,TEMPERATURE, PRESSURE
! REAL(SP), INTENT(OUT),DIMENSION(:,:) :: RHO ! RESULT - DENSITY
!=================================================================================
!=================================================================================
!=================================================================================
USE ALL_VARS
USE MOD_UTILS
IMPLICIT NONE
PUBLIC
INTERFACE FOFONOFF_MILLARD
MODULE PROCEDURE FOFONOFF_MILLARD_1D
MODULE PROCEDURE FOFONOFF_MILLARD_2D
END INTERFACE
INTERFACE DENS2G
MODULE PROCEDURE DENS2_1D
MODULE PROCEDURE DENS2_2D
END INTERFACE
INTERFACE JACKET_MCDOUGALL
MODULE PROCEDURE JACKET_MCDOUGALL_1D
MODULE PROCEDURE JACKET_MCDOUGALL_2D
END INTERFACE
PRIVATE JACKET_MCDOUGALL_1D
PRIVATE JACKET_MCDOUGALL_2D
PRIVATE DENS2_1D
PRIVATE DENS2_2D
PRIVATE FOFONOFF_MILLARD_1D
PRIVATE FOFONOFF_MILLARD_2D
PRIVATE SVAN
PRIVATE THETA
PRIVATE ATG
CONTAINS
!==============================================================================|
! Calculate Potential Density Based on Potential Temp and Salinity |
! Pressure effects are incorported (Can Model Fresh Water < 4 Deg C) |
! Ref: algorithms for computation of fundamental properties of |
! seawater , Fofonoff and Millard. |
! |
! calculates: rho1(nnode) density at nodes |
! calculates: rho (ncell) density at elements |
!==============================================================================|
SUBROUTINE FOFONOFF_MILLARD_2D(MYS,MYT,MYP,PREF, MYRHO)
!------------------------------------------------------------------------------|
IMPLICIT NONE
REAL(SP), INTENT(IN), DIMENSION(:,:) :: MYS,MYT,MYP
REAL(SP), INTENT(IN) :: PREF
REAL(SP), INTENT(INOUT), DIMENSION(:,:):: MYRHO
INTEGER :: I,K,ub1,lb1,ub2,lb2
REAL(SP) ::PT,SVA,SIGMA
!==============================================================================|
! SET DIMS USING MY S
lb1 = lbound(mys,1)
ub1 = ubound(mys,1)
lb2 = lbound(mys,2)
ub2 = ubound(mys,2)
! CHECK MY T
if(lb1 /= lbound(myt,1)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
if(ub1 /= ubound(myt,1)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
if(lb2 /= lbound(myt,2)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
if(ub2 /= ubound(myt,2)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
! CHECK MY P
if(lb1 /= lbound(myp,1)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
if(ub1 /= ubound(myp,1)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
if(lb2 /= lbound(myp,2)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
if(ub2 /= ubound(myp,2)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
! CHECK MY RHO
if(lb1 /= lbound(myrho,1)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
if(ub1 /= ubound(myrho,1)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
if(lb2 /= lbound(myrho,2)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
if(ub2 /= ubound(myrho,2)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
DO I=lb1,ub1
DO K=lb2,ub2
PT = THETA(MYS(I,K),MYT(I,K),MYP(I,K),PREF)
SVA = SVAN(MYS(I,K),PT,MYP(I,K),SIGMA)
MYRHO(I,K) = SIGMA*1.e-3_SP
END DO
END DO
!----------------TRANSFORM TO FACE CENTER--------------------------------------
END SUBROUTINE FOFONOFF_MILLARD_2D
!==============================================================================|
SUBROUTINE FOFONOFF_MILLARD_1D(MYS,MYT,MYP,PREF, MYRHO)
!------------------------------------------------------------------------------|
IMPLICIT NONE
REAL(SP), INTENT(IN), DIMENSION(:) :: MYS,MYT,MYP
REAL(SP), INTENT(IN) :: PREF
REAL(SP), INTENT(INOUT), DIMENSION(:):: MYRHO
INTEGER :: I,K,ub1,lb1
REAL(SP) ::PT,SVA,SIGMA
!==============================================================================|
! SET DIMS USING MY S
lb1 = lbound(mys,1)
ub1 = ubound(mys,1)
! CHECK MY T
if(lb1 /= lbound(myt,1)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
if(ub1 /= ubound(myt,1)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
! CHECK MY P
if(lb1 /= lbound(myp,1)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
if(ub1 /= ubound(myp,1)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
! CHECK MY RHO
if(lb1 /= lbound(myrho,1)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
if(ub1 /= ubound(myrho,1)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
DO K=lb1,ub1
PT = THETA(MYS(K),MYT(K),MYP(K),PREF)
SVA = SVAN(MYS(K),PT,MYP(K),SIGMA)
MYRHO(K) = SIGMA*1.e-3_SP
END DO
!----------------TRANSFORM TO FACE CENTER--------------------------------------
END SUBROUTINE FOFONOFF_MILLARD_1D
!==============================================================================|
! GENERIC CALL FOR FVCOM
!==============================================================================|
SUBROUTINE DENS1
!------------------------------------------------------------------------------|
IMPLICIT NONE
INTEGER :: K
REAL(SP), PARAMETER ::PR = 0.0_SP
REAL(SP), DIMENSION(0:MT,1:KB) :: RZU
!==============================================================================|
IF(DBG_SET(DBG_SBR)) WRITE(IPT,*) "START: DENS1"
! The thickness of the water column
! Is not the depth realtive to z=0
DO K=1,KBM1
RZU(:,K) = -GRAV_N*1.025_SP*(ZZ(:,K)*D(:))*0.1_SP
END DO
CALL FOFONOFF_MILLARD_2D(S1,T1,RZU,PR,RHO1)
RHO1(:,KB)=0.0_SP
RHO1(0,:)=0.0_SP
!----------------TRANSFORM TO FACE CENTER--------------------------------------
CALL N2E3D(RHO1,RHO)
IF(DBG_SET(DBG_SBR)) WRITE(IPT,*) "END: DENS1"
RETURN
END SUBROUTINE DENS1
!==============================================================================|
!==============================================================================|
! COMPUTE DENSITY USING SALINITY AND POTENTIAL TEMP |
! APPEARS TO BE BASED ON THE ROMS CODE? |
! CALCULATES: RHO1(M) DENSITY AT NODES |
! CALCULATES: RHO (N) DENSITY AT ELEMENTS |
!==============================================================================|
SUBROUTINE DENS2_2D(MYS,MYT,MYRHO)
!==============================================================================|
IMPLICIT NONE
REAL(SP), INTENT(IN), DIMENSION(:,:) :: MYS,MYT
REAL(SP), INTENT(INOUT), DIMENSION(:,:):: MYRHO
INTEGER :: I,K,ub1,lb1,ub2,lb2
!==============================================================================|
! SET DIMS USING MY S
lb1 = lbound(MYS,1)
ub1 = ubound(MYS,1)
lb2 = lbound(MYS,2)
ub2 = ubound(MYS,2)
! CHECK MY T
if(lb1 /= lbound(myt,1)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
if(ub1 /= ubound(myt,1)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
if(lb2 /= lbound(myt,2)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
if(ub2 /= ubound(myt,2)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
! CHECK MY RHO
if(lb1 /= lbound(myrho,1)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
if(ub1 /= ubound(myrho,1)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
if(lb2 /= lbound(myrho,2)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
if(ub2 /= ubound(myrho,2)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
!
! CALCULATE DENSITY FROM EQUATION OF STATE
!
MYRHO =MYS*MYS*MYS*6.76786136E-6_SP &
& - MYS*MYS*4.8249614E-4_SP &
& + MYS*8.14876577E-1_SP &
& - 0.22584586E0_SP
MYRHO = MYRHO * &
& ( MYT*MYT*MYT*1.667E-8_SP &
& - MYT*MYT*8.164E-7_SP &
& + MYT*1.803E-5_SP &
& )
MYRHO = MYRHO+1.0_SP &
& - MYT*MYT*MYT*1.0843E-6_SP &
& + MYT*MYT*9.8185E-5_SP &
& - MYT*4.786E-3_SP
MYRHO = MYRHO* &
& ( MYS*MYS*MYS*6.76786136E-6_SP &
& - MYS*MYS*4.8249614E-4_SP &
& + MYS*8.14876577E-1_SP &
& +3.895414E-2_SP &
& )
MYRHO = MYRHO &
& - (MYT-3.98_SP) *(MYT-3.98_SP) * (MYT+283.0_SP) / (503.57_SP*(MYT+67.26_SP))
!
! CALCULATE RHO
!
MYRHO = MYRHO*1.e-3_SP
END SUBROUTINE DENS2_2D
!==============================================================================|
SUBROUTINE DENS2_1D(MYS,MYT,MYRHO)
!==============================================================================|
IMPLICIT NONE
REAL(SP), INTENT(IN), DIMENSION(:) :: MYS,MYT
REAL(SP), INTENT(INOUT), DIMENSION(:):: MYRHO
INTEGER :: I,K,ub1,lb1,ub2,lb2
!==============================================================================|
! SET DIMS USING MY S
lb1 = lbound(MYS,1)
ub1 = ubound(MYS,1)
! CHECK MY T
if(lb1 /= lbound(myt,1)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
if(ub1 /= ubound(myt,1)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
! CHECK MY RHO
if(lb1 /= lbound(myrho,1)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
if(ub1 /= ubound(myrho,1)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
!
! CALCULATE DENSITY FROM EQUATION OF STATE
!
MYRHO =MYS*MYS*MYS*6.76786136E-6_SP &
& - MYS*MYS*4.8249614E-4_SP &
& + MYS*8.14876577E-1_SP &
& - 0.22584586E0_SP
MYRHO = MYRHO * &
& ( MYT*MYT*MYT*1.667E-8_SP &
& - MYT*MYT*8.164E-7_SP &
& + MYT*1.803E-5_SP &
& )
MYRHO = MYRHO+1.0_SP &
& - MYT*MYT*MYT*1.0843E-6_SP &
& + MYT*MYT*9.8185E-5_SP &
& - MYT*4.786E-3_SP
MYRHO = MYRHO* &
& ( MYS*MYS*MYS*6.76786136E-6_SP &
& - MYS*MYS*4.8249614E-4_SP &
& + MYS*8.14876577E-1_SP &
& +3.895414E-2_SP &
& )
MYRHO = MYRHO &
& - (MYT-3.98_SP) *(MYT-3.98_SP) * (MYT+283.0_SP) / (503.57_SP*(MYT+67.26_SP))
!
! CALCULATE RHO
!
MYRHO = MYRHO*1.e-3_SP
END SUBROUTINE DENS2_1D
!==============================================================================|
SUBROUTINE DENS2
!==============================================================================|
IMPLICIT NONE
!==============================================================================|
IF(DBG_SET(DBG_SBR)) WRITE(IPT,*) "START: DENS2"
!
! CALCULATE DENSITY FROM EQUATION OF STATE
!
CALL DENS2_2D(S1,T1,RHO1)
RHO1(:,KB)=0.0_SP
RHO1(0,:)=0.0_SP
!
! AVERAGE FROM NODES TO FACE CENTERS
!
CALL N2E3D(RHO1,RHO)
IF(DBG_SET(DBG_SBR)) WRITE(IPT,*) "END: DENS2"
RETURN
END SUBROUTINE DENS2
!==============================================================================|
!==============================================================================|
! COMPUTE IN SITU DENSITY - 1000 USING SALINITY, POTENTIAL TEMP, |
! AND PRESSURE FROM A POLYNOMIAL EXPRESSION (JACKETT & MCDOUGALL, |
! 1995). IT ASSUMES NO PRESSURE VARIATION ALONG GEOPOTENTIAL |
! SURFACES, THAT IS, DEPTH (METERS; NEGATIVE) AND PRESSURE (DBAR |
! ASSUMED NEGATIVE HERE) ARE INTERCHANGEABLE. |
! |
! check Values: (T=3 C, S=35.5 PSU, Z=-5000 m) |
! RHOF = 1050.3639165364 (kg/m3) |
! DEN1 = 1028.2845117925 (kg/m3) |
! |
! Reference: |
! |
! Jackett, D. R. and T. J. McDougall, 1995, Minimal Adjustment of |
! Hydrostatic Profiles to Achieve Static Stability, J. of Atmos. |
! and Oceanic Techn., vol. 12, pp. 381-389. |
! |
! CALCULATES: RHO1(M) DENSITY AT NODES |
! CALCULATES: RHO (N) DENSITY AT ELEMENTS |
!==============================================================================|
SUBROUTINE JACKET_MCDOUGALL_2D(MYS,MYT,MYP,MYRHO)
!==============================================================================|
IMPLICIT NONE
REAL(SP), INTENT(IN), DIMENSION(:,:) :: MYS,MYT,MYP
REAL(SP), INTENT(INOUT), DIMENSION(:,:):: MYRHO
REAL(SP) :: TF,SF,sqrtSF,PBAR,TEMP(10),BULK,BULK0,BULK1,BULK2
INTEGER :: I,K,ub1,lb1,ub2,lb2
!==============================================================================|
! Polynomial expansion coefficients for the computation of in situ |
! density via the nonlinear equation of state for seawater as a |
! function of potential temperature, salinity, and pressure (Jackett |
! and McDougall, 1995). |
REAL(SP), PARAMETER :: A00 = +1.965933e+04_SP
REAL(SP), PARAMETER :: A01 = +1.444304e+02_SP
REAL(SP), PARAMETER :: A02 = -1.706103e+00_SP
REAL(SP), PARAMETER :: A03 = +9.648704e-03_SP
REAL(SP), PARAMETER :: A04 = -4.190253e-05_SP
REAL(SP), PARAMETER :: B00 = +5.284855e+01_SP
REAL(SP), PARAMETER :: B01 = -3.101089e-01_SP
REAL(SP), PARAMETER :: B02 = +6.283263e-03_SP
REAL(SP), PARAMETER :: B03 = -5.084188e-05_SP
REAL(SP), PARAMETER :: D00 = +3.886640e-01_SP
REAL(SP), PARAMETER :: D01 = +9.085835e-03_SP
REAL(SP), PARAMETER :: D02 = -4.619924e-04_SP
REAL(SP), PARAMETER :: E00 = +3.186519e+00_SP
REAL(SP), PARAMETER :: E01 = +2.212276e-02_SP
REAL(SP), PARAMETER :: E02 = -2.984642e-04_SP
REAL(SP), PARAMETER :: E03 = +1.956415e-06_SP
REAL(SP), PARAMETER :: F00 = +6.704388e-03_SP
REAL(SP), PARAMETER :: F01 = -1.847318e-04_SP
REAL(SP), PARAMETER :: F02 = +2.059331e-07_SP
REAL(SP), PARAMETER :: G00 = +1.480266e-04_SP
REAL(SP), PARAMETER :: G01 = +2.102898e-04_SP
REAL(SP), PARAMETER :: G02 = -1.202016e-05_SP
REAL(SP), PARAMETER :: G03 = +1.394680e-07_SP
REAL(SP), PARAMETER :: H00 = -2.040237e-06_SP
REAL(SP), PARAMETER :: H01 = +6.128773e-08_SP
REAL(SP), PARAMETER :: H02 = +6.207323e-10_SP
REAL(SP), PARAMETER :: Q00 = +9.99842594e+02_SP
REAL(SP), PARAMETER :: Q01 = +6.793952e-02_SP
REAL(SP), PARAMETER :: Q02 = -9.095290e-03_SP
REAL(SP), PARAMETER :: Q03 = +1.001685e-04_SP
REAL(SP), PARAMETER :: Q04 = -1.120083e-06_SP
REAL(SP), PARAMETER :: Q05 = +6.536332e-09_SP
REAL(SP), PARAMETER :: U00 = +8.24493e-01_SP
REAL(SP), PARAMETER :: U01 = -4.08990e-03_SP
REAL(SP), PARAMETER :: U02 = +7.64380e-05_SP
REAL(SP), PARAMETER :: U03 = -8.24670e-07_SP
REAL(SP), PARAMETER :: U04 = +5.38750e-09_SP
REAL(SP), PARAMETER :: V00 = -5.72466e-03_SP
REAL(SP), PARAMETER :: V01 = +1.02270e-04_SP
REAL(SP), PARAMETER :: V02 = -1.65460e-06_SP
REAL(SP), PARAMETER :: W00 = +4.8314e-04_SP
! SET DIMS USING MY S
lb1 = lbound(mys,1)
ub1 = ubound(mys,1)
lb2 = lbound(mys,2)
ub2 = ubound(mys,2)
! CHECK MY T
if(lb1 /= lbound(myt,1)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
if(ub1 /= ubound(myt,1)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
if(lb2 /= lbound(myt,2)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
if(ub2 /= ubound(myt,2)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
! CHECK MY P
if(lb1 /= lbound(myp,1)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
if(ub1 /= ubound(myp,1)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
if(lb2 /= lbound(myp,2)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
if(ub2 /= ubound(myp,2)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
! CHECK MY RHO
if(lb1 /= lbound(myrho,1)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
if(ub1 /= ubound(myrho,1)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
if(lb2 /= lbound(myrho,2)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
if(ub2 /= ubound(myrho,2)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
!
! CALCULATE DENSITY FROM EQUATION OF STATE
!
DO I=lb1,ub1
DO K=lb2,ub2
TF = MYT(I,K)
SF = MYS(I,K)
sqrtSF = sqrt(SF)
PBAR = MYP(I,K)
! Compute density (kg/m3) at standard one atmosphere pressure
TEMP(1)=Q00+TF*(Q01+TF*(Q02+TF*(Q03+TF*(Q04+TF*Q05))))
TEMP(2)=U00+TF*(U01+TF*(U02+TF*(U03+TF*U04)))
TEMP(3)=V00+TF*(V01+TF*V02)
MYRHO(I,K)=TEMP(1)+SF*(TEMP(2)+sqrtSF*TEMP(3)+SF*W00)
! Compute secant bulk modulus (BULK = BULK0 + BULK1*PBAR + BULK2*PBAR*PBAR)
TEMP(4)=A00+TF*(A01+TF*(A02+TF*(A03+TF*A04)))
TEMP(5)=B00+TF*(B01+TF*(B02+TF*B03))
TEMP(6)=D00+TF*(D01+TF*D02)
TEMP(7)=E00+TF*(E01+TF*(E02+TF*E03))
TEMP(8)=F00+TF*(F01+TF*F02)
TEMP(9)=G01+TF*(G02+TF*G03)
TEMP(10)=H00+TF*(H01+TF*H02)
BULK0=TEMP(4)+SF*(TEMP(5)+sqrtSF*TEMP(6))
BULK1=TEMP(7)+SF*(TEMP(8)+sqrtSF*G00)
BULK2=TEMP(9)+SF*TEMP(10)
BULK = BULK0 + PBAR * (BULK1 + PBAR * BULK2)
! Compute "in situ" density anomaly (kg/m3)
MYRHO(I,K)=(MYRHO(I,K)*BULK)/(BULK-PBAR)
MYRHO(I,K)= MYRHO(I,K)-1000.0_SP
END DO
END DO
!
! CALCULATE RHO1
!
MYRHO = MYRHO*1.e-3_SP
END SUBROUTINE JACKET_MCDOUGALL_2D
!==============================================================================|
SUBROUTINE JACKET_MCDOUGALL_1D(MYS,MYT,MYP,MYRHO)
!==============================================================================|
IMPLICIT NONE
REAL(SP), INTENT(IN), DIMENSION(:) :: MYS,MYT,MYP
REAL(SP), INTENT(INOUT), DIMENSION(:):: MYRHO
REAL(SP) :: TF,SF,sqrtSF,PBAR,TEMP(10),BULK,BULK0,BULK1,BULK2
INTEGER :: K,ub1,lb1
!==============================================================================|
! Polynomial expansion coefficients for the computation of in situ |
! density via the nonlinear equation of state for seawater as a |
! function of potential temperature, salinity, and pressure (Jackett |
! and McDougall, 1995). |
REAL(SP), PARAMETER :: A00 = +1.965933e+04_SP
REAL(SP), PARAMETER :: A01 = +1.444304e+02_SP
REAL(SP), PARAMETER :: A02 = -1.706103e+00_SP
REAL(SP), PARAMETER :: A03 = +9.648704e-03_SP
REAL(SP), PARAMETER :: A04 = -4.190253e-05_SP
REAL(SP), PARAMETER :: B00 = +5.284855e+01_SP
REAL(SP), PARAMETER :: B01 = -3.101089e-01_SP
REAL(SP), PARAMETER :: B02 = +6.283263e-03_SP
REAL(SP), PARAMETER :: B03 = -5.084188e-05_SP
REAL(SP), PARAMETER :: D00 = +3.886640e-01_SP
REAL(SP), PARAMETER :: D01 = +9.085835e-03_SP
REAL(SP), PARAMETER :: D02 = -4.619924e-04_SP
REAL(SP), PARAMETER :: E00 = +3.186519e+00_SP
REAL(SP), PARAMETER :: E01 = +2.212276e-02_SP
REAL(SP), PARAMETER :: E02 = -2.984642e-04_SP
REAL(SP), PARAMETER :: E03 = +1.956415e-06_SP
REAL(SP), PARAMETER :: F00 = +6.704388e-03_SP
REAL(SP), PARAMETER :: F01 = -1.847318e-04_SP
REAL(SP), PARAMETER :: F02 = +2.059331e-07_SP
REAL(SP), PARAMETER :: G00 = +1.480266e-04_SP
REAL(SP), PARAMETER :: G01 = +2.102898e-04_SP
REAL(SP), PARAMETER :: G02 = -1.202016e-05_SP
REAL(SP), PARAMETER :: G03 = +1.394680e-07_SP
REAL(SP), PARAMETER :: H00 = -2.040237e-06_SP
REAL(SP), PARAMETER :: H01 = +6.128773e-08_SP
REAL(SP), PARAMETER :: H02 = +6.207323e-10_SP
REAL(SP), PARAMETER :: Q00 = +9.99842594e+02_SP
REAL(SP), PARAMETER :: Q01 = +6.793952e-02_SP
REAL(SP), PARAMETER :: Q02 = -9.095290e-03_SP
REAL(SP), PARAMETER :: Q03 = +1.001685e-04_SP
REAL(SP), PARAMETER :: Q04 = -1.120083e-06_SP
REAL(SP), PARAMETER :: Q05 = +6.536332e-09_SP
REAL(SP), PARAMETER :: U00 = +8.24493e-01_SP
REAL(SP), PARAMETER :: U01 = -4.08990e-03_SP
REAL(SP), PARAMETER :: U02 = +7.64380e-05_SP
REAL(SP), PARAMETER :: U03 = -8.24670e-07_SP
REAL(SP), PARAMETER :: U04 = +5.38750e-09_SP
REAL(SP), PARAMETER :: V00 = -5.72466e-03_SP
REAL(SP), PARAMETER :: V01 = +1.02270e-04_SP
REAL(SP), PARAMETER :: V02 = -1.65460e-06_SP
REAL(SP), PARAMETER :: W00 = +4.8314e-04_SP
! SET DIMS USING MY S
lb1 = lbound(mys,1)
ub1 = ubound(mys,1)
! CHECK MY T
if(lb1 /= lbound(myt,1)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
if(ub1 /= ubound(myt,1)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
! CHECK MY P
if(lb1 /= lbound(myp,1)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
if(ub1 /= ubound(myp,1)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
! CHECK MY RHO
if(lb1 /= lbound(myrho,1)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
if(ub1 /= ubound(myrho,1)) CALL FATAL_ERROR("EQS_OF_STATE:: Dimension mismatch!")
!
! CALCULATE DENSITY FROM EQUATION OF STATE
!
DO K=lb1,ub1
TF = MYT(K)
SF = MYS(K)
sqrtSF = sqrt(SF)
PBAR = MYP(K)
! Compute density (kg/m3) at standard one atmosphere pressure
TEMP(1)=Q00+TF*(Q01+TF*(Q02+TF*(Q03+TF*(Q04+TF*Q05))))
TEMP(2)=U00+TF*(U01+TF*(U02+TF*(U03+TF*U04)))
TEMP(3)=V00+TF*(V01+TF*V02)
MYRHO(K)=TEMP(1)+SF*(TEMP(2)+sqrtSF*TEMP(3)+SF*W00)
! Compute secant bulk modulus (BULK = BULK0 + BULK1*PBAR + BULK2*PBAR*PBAR)
TEMP(4)=A00+TF*(A01+TF*(A02+TF*(A03+TF*A04)))
TEMP(5)=B00+TF*(B01+TF*(B02+TF*B03))
TEMP(6)=D00+TF*(D01+TF*D02)
TEMP(7)=E00+TF*(E01+TF*(E02+TF*E03))
TEMP(8)=F00+TF*(F01+TF*F02)
TEMP(9)=G01+TF*(G02+TF*G03)
TEMP(10)=H00+TF*(H01+TF*H02)
BULK0=TEMP(4)+SF*(TEMP(5)+sqrtSF*TEMP(6))
BULK1=TEMP(7)+SF*(TEMP(8)+sqrtSF*G00)
BULK2=TEMP(9)+SF*TEMP(10)
BULK = BULK0 + PBAR * (BULK1 + PBAR * BULK2)
! Compute "in situ" density anomaly (kg/m3)
MYRHO(K)=(MYRHO(K)*BULK)/(BULK-PBAR)
MYRHO(K)= MYRHO(K)-1000.0_SP
END DO
!
! CALCULATE RHO1
!
MYRHO = MYRHO*1.e-3_SP
END SUBROUTINE JACKET_MCDOUGALL_1D
!==============================================================================|
SUBROUTINE DENS3
!==============================================================================|
IMPLICIT NONE
REAL(SP), DIMENSION(0:MT,KB) :: MYP
INTEGER :: K
IF(DBG_SET(DBG_SBR)) WRITE(IPT,*) "START: DENS3"
DO K=1,KBM1
MYP(:,K) = -GRAV_N(:)*1.025_SP*ZZ(:,K)*D(:) *0.01_SP
END DO
CALL JACKET_MCDOUGALL_2D(S1,T1,MYP,RHO1)
RHO1(:,KB)=0.0_SP
RHO1(0,:)=0.0_SP
!
! AVERAGE FROM NODES TO FACE CENTERS
!
CALL N2E3D(RHO1,RHO)
IF(DBG_SET(DBG_SBR)) WRITE(IPT,*) "END: DENS3"
RETURN
END SUBROUTINE DENS3
!==============================================================================!
FUNCTION SVAN(S4,T4,P04,SIGMA)
!==============================================================================!
! specific volume anomaly (steric anomaly) based on 1980 equation |
! of state for seawater and 1978 practerical salinity scale. |
! references: |
! millero, et al (1980) deep-sea res.,27a,255-264 |
! millero and poisson 1981,deep-sea res.,28a pp 625-629. |
! both above references are also found in unesco report 38 (1981) |
! |
! units: |
! pressure p04 decibars |
! temperature t4 deg celsius (ipts-68) |
! salinity s4 (ipss-78) |
! spec. vol. ana. svan m**3/kg *1.0e-8 |
! density ana. sigma kg/m**3 |
! |
! check value: svan=981.3021 e-8 m**3/kg. for s = 40 (ipss-78), |
! t = 40 deg c, p0= 10000 decibars. |
! check value: sigma = 59.82037 kg/m**3. for s = 40 (ipss-78) , |
! t = 40 deg c, p0= 10000 decibars. |
!==============================================================================!
USE MOD_PREC
IMPLICIT NONE
REAL(SP) :: SVAN
REAL(SP), INTENT(IN) :: S4,T4,P04
REAL(SP), INTENT(OUT) :: SIGMA
REAL(SP) P4,SIG,SR,RR1,RR2,RR3,V350P,DK
REAL(SP) A4,B4,C4,D4,E4,AA1,BB1,AW,BW,KK,K0,KW,K35,SVA
REAL(SP) GAM,PK,DVAN,DR35P
REAL(SP), PARAMETER :: R3500 = 1028.1063_SP
REAL(SP), PARAMETER :: RR4 = 4.8314E-4_SP
REAL(SP), PARAMETER :: DR350 = 28.106331_SP
! rr4 is refered to as c in millero and poisson 1981
! convert pressure to bars and take square root salinity.
P4=P04/10.0_SP
SR = SQRT(ABS(S4))
! pure water density at atmospheric pressure
! bigg p.h.,(1967) br. j. applied physics 8 pp 521-537.
!
RR1=((((6.536332E-9_SP*T4-1.120083E-6_SP)*T4+1.001685E-4_SP)*T4 &
-9.095290E-3_SP)*T4+6.793952E-2_SP)*T4-28.263737_SP
! seawater density atm press.
! coefficients involving salinity
! rr2 = a in notation of millero and poisson 1981
RR2=(((5.3875E-9_SP*T4-8.2467E-7_SP)*T4+7.6438E-5_SP)*T4-4.0899E-3_SP)*T4 &
+8.24493E-1_SP
! rr3 = b4 in notation of millero and poisson 1981
RR3=(-1.6546E-6_SP*T4+1.0227E-4_SP)*T4-5.72466E-3_SP
! international one-atmosphere equation of state of seawater
SIG=(RR4*S4+RR3*SR+RR2)*S4+RR1
! specific volume at atmospheric pressure
V350P = 1.0_SP/R3500
SVA = -SIG*V350P/(R3500+SIG)
SIGMA=SIG+DR350
! scale specific vol. anamoly to normally reported units
SVAN=SVA*1.0E+8_SP
IF(P4 == 0.0_SP) RETURN
!-------------------------------------------------------------|
! new high pressure equation of sate for seawater |
! |
! millero, el al., 1980 dsr 27a, pp 255-264 |
! constant notation follows article |
!-------------------------------------------------------------|
! compute compression terms
E4 = (9.1697E-10*T4+2.0816E-8_SP)*T4-9.9348E-7_SP
BW = (5.2787E-8_SP*T4-6.12293E-6_SP)*T4+3.47718E-5_SP
B4 = BW + E4*S4
D4 = 1.91075E-4
C4 = (-1.6078E-6_SP*T4-1.0981E-5_SP)*T4+2.2838E-3_SP
AW = ((-5.77905E-7_SP*T4+1.16092E-4_SP)*T4+1.43713E-3_SP)*T4 &
-0.1194975_SP
A4 = (D4*SR + C4)*S4 + AW
BB1 = (-5.3009E-4_SP*T4+1.6483E-2_SP)*T4+7.944E-2_SP
AA1 = ((-6.1670E-5_SP*T4+1.09987E-2_SP)*T4-0.603459_SP)*T4+54.6746
KW = (((-5.155288E-5_SP*T4+1.360477E-2_SP)*T4-2.327105_SP)*T4 &
+148.4206_SP)*T4-1930.06_SP
K0 = (BB1*SR + AA1)*S4 + KW
! evaluate pressure polynomial
!-----------------------------------------------------|
! k equals the secant bulk modulus of seawater |
! dk=k(s,t,p)-k(35,0,p) |
! k35=k(35,0,p) |
!-----------------------------------------------------|
DK = (B4*P4 + A4)*P4 + K0
K35 = (5.03217E-5_SP*P4+3.359406_SP)*P4+21582.27_SP
GAM=P4/K35
PK = 1.0_SP - GAM
SVA = SVA*PK + (V350P+SVA)*P4*DK/(K35*(K35+DK))
! scale specific vol. anamoly to normally reported units
SVAN=SVA*1.0E+8_SP
V350P = V350P*PK
!----------------------------------------------------------|
! compute density anamoly with respect to 1000.0 kg/m**3 |
! 1) dr350: density anamoly at 35 (ipss-78), |
! 0 deg. c and 0 decibars |
! 2) dr35p: density anamoly at 35 (ipss-78), |
! 0 deg. c, pres. variation |
! 3) dvan : density anamoly variations involving specific |
! volume anamoly |
! |
! check values: sigma = 59.82037 kg/m**3 |
! for s = 40 (ipss-78), t = 40 deg c, p0= 10000 decibars. |
!----------------------------------------------------------|
DR35P=GAM/V350P
DVAN=SVA/(V350P*(V350P+SVA))
SIGMA=DR350+DR35P-DVAN
RETURN
END FUNCTION SVAN
!==============================================================================!
!==============================================================================!
!==============================================================================|
FUNCTION THETA(S4,T04,P04,PR)
!==============================================================================|
! to compute local potential temperature at pr using |
! bryden 1973 polynomial for adiabatic lapse rate and |
! runge-kutta 4th order integration algorithm. |
! ref: bryden,h.,1973,deep-sea res.,20,401-408; |
! fofonoff,n.,1977,deep-sea res.,24,489-491 |
! |
! units: |
! pressure p04 decibars |
! temperature t04 deg celsius (ipts-68) |
! salinity s4 (ipss-78) |
! reference prs pr decibars |
! potential tmp. theta deg celsius |
! checkvalue: |
! theta= 36.89073 c,s=40 (ipss-78), |
! t0=40 deg c,p0=10000 decibars,pr=0 decibars |
! |
! |
! set up intermediate temperature and pressure variables. |
!==============================================================================|
USE MOD_PREC
IMPLICIT NONE
REAL(SP) :: THETA
REAL(SP), INTENT(IN) :: S4,T04,P04,PR
REAL(SP) :: P4,T4,H4,Q4,XK
!==============================================================================|
P4 = P04
T4 = T04
H4 = PR - P4
XK = H4*ATG(S4,T4,P4)
T4 = T4 + .5_SP*XK
Q4 = XK
P4 = P4 + 0.5_SP*H4
XK = H4*ATG(S4,T4,P4)
T4 = T4 + 0.29289322_SP*(XK-Q4)
Q4 = 0.58578644_SP*XK + 0.121320344_SP*Q4
XK = H4*ATG(S4,T4,P4)
T4 = T4 + 1.707106781_SP*(XK-Q4)
Q4 = 3.414213562_SP*XK - 4.121320344_SP*Q4
P4 = P4 + 0.5_SP*H4
XK = H4*ATG(S4,T4,P4)
THETA = T4 + (XK-2.0_SP*Q4)/6.0_SP
RETURN
END FUNCTION THETA
!==============================================================================|
!==============================================================================|
! adiabatic temperature gradient deg c per decibar |
! ref: bryden, h., 1973,deep-sea res.,20,401-408 |
! |
! units: |
! pressure P4 decibars |
! temperature T4 deg celsius(ipts-68) |
! salinity s4 (ipss-78) |
! adiabatic atg deg. c/decibar |
! checkvalue: atg=3.255976e-4 c/dbar for s=40 (ipss-78), |
! t=40 deg c,p0=10000 decibars |
!==============================================================================|
REAL(SP) FUNCTION ATG(S4,T4,P4)
USE MOD_PREC
!------------------------------------------------------------------------------|
IMPLICIT NONE
REAL(SP), INTENT(IN) :: S4,T4,P4
REAL(SP) :: DS
!==============================================================================|
DS = S4 - 35.0_SP
ATG = (((-2.1687e-16_SP*T4+1.8676e-14_SP)*T4-4.6206e-13_SP)*P4 &
+((2.7759e-12_SP*T4-1.1351e-10_SP)*DS+((-5.4481e-14_SP*T4 &
+8.733e-12_SP)*T4-6.7795e-10_SP)*T4+1.8741e-8_SP))*P4 &
+(-4.2393e-8_SP*T4+1.8932e-6_SP)*DS &
+((6.6228e-10_SP*T4-6.836e-8_SP)*T4+8.5258e-6_SP)*T4+3.5803e-5_SP
RETURN
END FUNCTION ATG
!==============================================================================|
END MODULE EQS_OF_STATE