/[MITgcm]/MITgcm/model/src/seawater.F
ViewVC logotype

Diff of /MITgcm/model/src/seawater.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.4 by ce107, Tue Jun 30 20:47:32 2009 UTC revision 1.5 by jmc, Thu Oct 22 04:33:59 2009 UTC
# Line 5  C $Name$ Line 5  C $Name$
5    
6  C--  File seawater.F: routines that compute quantities related to seawater.  C--  File seawater.F: routines that compute quantities related to seawater.
7  C--   Contents  C--   Contents
 C--   o FIND_RHO_SCALAR: in-situ density for individual points  
8  C--   o SW_PTMP: function to compute potential temperature  C--   o SW_PTMP: function to compute potential temperature
9  C--   o SW_TEMP: function to compute potential temperature  C--   o SW_TEMP: function to compute in-situ temperature from pot. temp.
10  C--   o SW_ADTG: function to compute adiabatic tmperature gradient  C--   o SW_ADTG: function to compute adiabatic temperature gradient
11  C--              used by sw_ptmp  C--              (used by both SW_PTMP & SW_TEMP)
12    
13        SUBROUTINE FIND_RHO_SCALAR(  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
14       I     tLoc, sLoc, pLoc,  
15       O     rhoLoc,  CBOP
16       I     myThid )  C     !ROUTINE: SW_PTMP
17    C     !INTERFACE:
18          _RL FUNCTION SW_PTMP  (S,T,P,PR)
19    
20  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
21  C     *==========================================================*  C     *=============================================================*
22  C     | o SUBROUTINE FIND_RHO_SCALAR  C     | S/R  SW_PTMP
23  C     |   Calculates rho(S,T,p)  C     | o compute potential temperature as per UNESCO 1983 report.
24  C     *==========================================================*  C     *=============================================================*
25  C     \ev  C
26    C     started:
27    C              Armin Koehl akoehl@ucsd.edu
28    C
29    C     ==================================================================
30    C     SUBROUTINE SW_PTMP
31    C     ==================================================================
32    C     S  :: salinity    [psu      (PSS-78) ]
33    C     T  :: temperature [degree C (IPTS-68)]
34    C     P  :: pressure    [db]
35    C     PR :: Reference pressure  [db]
36    
37  C     !USES:  C     !USES:
38        IMPLICIT NONE        IMPLICIT NONE
 C     == Global variables ==  
 #include "SIZE.h"  
 #include "EEPARAMS.h"  
 #include "PARAMS.h"  
 #include "EOS.h"  
39    
40  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
 C     == Routine arguments ==  
       _RL sLoc, tLoc, pLoc  
       _RL rhoLoc  
       INTEGER myThid  
   
 C     !LOCAL VARIABLES:  
 C     == Local variables ==  
   
       _RL t1, t2, t3, t4, s1, s3o2, p1, p2, sp5, p1t1  
       _RL rfresh, rsalt, rhoP0  
       _RL bMfresh, bMsalt, bMpres, BulkMod  
       _RL rhoNum, rhoDen, den, epsln  
       parameter ( epsln = 0. _d 0 )  
   
       CHARACTER*(MAX_LEN_MBUF) msgBuf  
 CEOP  
   
       rhoLoc  = 0. _d 0  
       rhoP0   = 0. _d 0  
       bulkMod = 0. _d 0  
       rfresh  = 0. _d 0  
       rsalt   = 0. _d 0  
       bMfresh = 0. _d 0  
       bMsalt  = 0. _d 0  
       bMpres  = 0. _d 0  
       rhoNum  = 0. _d 0  
       rhoDen  = 0. _d 0  
       den     = 0. _d 0  
   
       t1 = tLoc  
       t2 = t1*t1  
       t3 = t2*t1  
       t4 = t3*t1  
   
       s1  = sLoc  
       IF ( s1 .LT. 0. _d 0 ) THEN  
 C     issue a warning  
          WRITE(msgBuf,'(A,E13.5)')  
      &        ' FIND_RHO_SCALAR:   WARNING, salinity = ', s1  
          CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,  
      &                       SQUEEZE_RIGHT , myThid )  
          s1 = 0. _d 0  
       ENDIF  
   
       IF (equationOfState.EQ.'LINEAR') THEN  
   
          rholoc = rhoNil*(  
      &                      sBeta *(sLoc-sRef(1))  
      &                     -tAlpha*(tLoc-tRef(1))  
      &                   ) + rhoNil  
 c        rhoLoc = 0. _d  0  
   
       ELSEIF (equationOfState.EQ.'POLY3') THEN  
   
 C     this is not correct, there is a field eosSig0 which should be use here  
 C     but I DO not intent to include the reference level in this routine  
          WRITE(msgBuf,'(A)')  
      &        ' FIND_RHO_SCALAR: for POLY3, the density is not'  
          CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,  
      &                       SQUEEZE_RIGHT , myThid )  
          WRITE(msgBuf,'(A)')  
      &         '                 computed correctly in this routine'  
          CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,  
      &                       SQUEEZE_RIGHT , myThid )  
          rhoLoc = 0. _d 0  
   
       ELSEIF ( equationOfState(1:5).EQ.'JMD95'  
      &      .OR. equationOfState.EQ.'UNESCO' ) THEN  
 C     nonlinear equation of state in pressure coordinates  
   
          s3o2 = s1*SQRT(s1)  
   
          p1 = pLoc*SItoBar  
          p2 = p1*p1  
   
 C     density of freshwater at the surface  
          rfresh =  
      &          eosJMDCFw(1)  
      &        + eosJMDCFw(2)*t1  
      &        + eosJMDCFw(3)*t2  
      &        + eosJMDCFw(4)*t3  
      &        + eosJMDCFw(5)*t4  
      &        + eosJMDCFw(6)*t4*t1  
 C     density of sea water at the surface  
          rsalt =  
      &        s1*(  
      &             eosJMDCSw(1)  
      &           + eosJMDCSw(2)*t1  
      &           + eosJMDCSw(3)*t2  
      &           + eosJMDCSw(4)*t3  
      &           + eosJMDCSw(5)*t4  
      &           )  
      &        + s3o2*(  
      &             eosJMDCSw(6)  
      &           + eosJMDCSw(7)*t1  
      &           + eosJMDCSw(8)*t2  
      &           )  
      &           + eosJMDCSw(9)*s1*s1  
   
          rhoP0 = rfresh + rsalt  
   
 C     secant bulk modulus of fresh water at the surface  
          bMfresh =  
      &             eosJMDCKFw(1)  
      &           + eosJMDCKFw(2)*t1  
      &           + eosJMDCKFw(3)*t2  
      &           + eosJMDCKFw(4)*t3  
      &           + eosJMDCKFw(5)*t4  
 C     secant bulk modulus of sea water at the surface  
          bMsalt =  
      &        s1*( eosJMDCKSw(1)  
      &           + eosJMDCKSw(2)*t1  
      &           + eosJMDCKSw(3)*t2  
      &           + eosJMDCKSw(4)*t3  
      &           )  
      &    + s3o2*( eosJMDCKSw(5)  
      &           + eosJMDCKSw(6)*t1  
      &           + eosJMDCKSw(7)*t2  
      &           )  
 C     secant bulk modulus of sea water at pressure p  
          bMpres =  
      &        p1*( eosJMDCKP(1)  
      &           + eosJMDCKP(2)*t1  
      &           + eosJMDCKP(3)*t2  
      &           + eosJMDCKP(4)*t3  
      &           )  
      &   + p1*s1*( eosJMDCKP(5)  
      &           + eosJMDCKP(6)*t1  
      &           + eosJMDCKP(7)*t2  
      &           )  
      &      + p1*s3o2*eosJMDCKP(8)  
      &      + p2*( eosJMDCKP(9)  
      &           + eosJMDCKP(10)*t1  
      &           + eosJMDCKP(11)*t2  
      &           )  
      &    + p2*s1*( eosJMDCKP(12)  
      &           + eosJMDCKP(13)*t1  
      &           + eosJMDCKP(14)*t2  
      &           )  
   
          bulkMod = bMfresh + bMsalt + bMpres  
   
 C     density of sea water at pressure p  
          rhoLoc = rhoP0/(1. _d 0 - p1/bulkMod)  
   
       ELSEIF ( equationOfState.EQ.'MDJWF' ) THEN  
   
          sp5 = SQRT(s1)  
   
          p1   = pLoc*SItodBar  
          p1t1 = p1*t1  
   
          rhoNum = eosMDJWFnum(0)  
      &        + t1*(eosMDJWFnum(1)  
      &        +     t1*(eosMDJWFnum(2) + eosMDJWFnum(3)*t1) )  
      &        + s1*(eosMDJWFnum(4)  
      &        +     eosMDJWFnum(5)*t1  + eosMDJWFnum(6)*s1)  
      &        + p1*(eosMDJWFnum(7) + eosMDJWFnum(8)*t2  
      &        +     eosMDJWFnum(9)*s1  
      &        +     p1*(eosMDJWFnum(10) + eosMDJWFnum(11)*t2) )  
   
   
          den = eosMDJWFden(0)  
      &        + t1*(eosMDJWFden(1)  
      &        +     t1*(eosMDJWFden(2)  
      &        +         t1*(eosMDJWFden(3) + t1*eosMDJWFden(4) ) ) )  
      &        + s1*(eosMDJWFden(5)  
      &        +     t1*(eosMDJWFden(6)  
      &        +         eosMDJWFden(7)*t2)  
      &        +     sp5*(eosMDJWFden(8) + eosMDJWFden(9)*t2) )  
      &        + p1*(eosMDJWFden(10)  
      &        +     p1t1*(eosMDJWFden(11)*t2 + eosMDJWFden(12)*p1) )  
   
          rhoDen = 1.0/(epsln+den)  
   
          rhoLoc = rhoNum*rhoDen  
   
       ELSEIF( equationOfState .EQ. 'IDEALG' ) THEN  
 C  
       ELSE  
        WRITE(msgBuf,'(3A)')  
      &        ' FIND_RHO_SCALAR : equationOfState = "',  
      &        equationOfState,'"'  
        CALL PRINT_ERROR( msgBuf, myThid )  
        STOP 'ABNORMAL END: S/R FIND_RHO_SCALAR'  
       ENDIF  
   
       RETURN  
       END  
   
 C=================================================================  
   
       _RL FUNCTION SW_PTMP  (S,T,P,PR)  
   
 c     ==================================================================  
 c     SUBROUTINE SW_PTMP  
 c     ==================================================================  
 c  
 c     o   Calculates potential temperature as per UNESCO 1983 report.  
 c  
 c     started:  
 c  
 c              Armin Koehl akoehl@ucsd.edu  
 c  
 c     ==================================================================  
 c     SUBROUTINE SW_PTMP  
 c     ==================================================================  
 C     S  = salinity    [psu      (PSS-78) ]  
 C     T  = temperature [degree C (IPTS-68)]  
 C     P  = pressure    [db]  
 C     PR = Reference pressure  [db]  
   
       implicit none  
   
 c     routine arguments  
41        _RL S,T,P,PR        _RL S,T,P,PR
42    
43  c     local arguments  C     !FUNCTIONS:
44          _RL sw_adtg
45          EXTERNAL sw_adtg
46    
47    C     !LOCAL VARIABLES
48        _RL del_P ,del_th, th, q        _RL del_P ,del_th, th, q
49        _RL onehalf, two, three        _RL onehalf, two, three
50        parameter ( onehalf = 0.5 _d 0, two = 2. _d 0, three = 3. _d 0 )        PARAMETER ( onehalf = 0.5 _d 0, two = 2. _d 0, three = 3. _d 0 )
51    CEOP
52    
53  c     externals  C theta1
       _RL sw_adtg  
       external sw_adtg  
 c theta1  
54        del_P  = PR - P        del_P  = PR - P
55        del_th = del_P*sw_adtg(S,T,P)        del_th = del_P*sw_adtg(S,T,P)
56        th     = T + onehalf*del_th        th     = T + onehalf*del_th
57        q      = del_th        q      = del_th
58  c theta2  C theta2
59        del_th = del_P*sw_adtg(S,th,P+onehalf*del_P)        del_th = del_P*sw_adtg(S,th,P+onehalf*del_P)
60    
61        th     = th + (1 - 1/sqrt(two))*(del_th - q)        th     = th + (1 - 1/sqrt(two))*(del_th - q)
62        q      = (two-sqrt(two))*del_th + (-two+three/sqrt(two))*q        q      = (two-sqrt(two))*del_th + (-two+three/sqrt(two))*q
63    
64  c theta3  C theta3
65        del_th = del_P*sw_adtg(S,th,P+onehalf*del_P)        del_th = del_P*sw_adtg(S,th,P+onehalf*del_P)
66        th     = th + (1 + 1/sqrt(two))*(del_th - q)        th     = th + (1 + 1/sqrt(two))*(del_th - q)
67        q      = (two + sqrt(two))*del_th + (-two-three/sqrt(two))*q        q      = (two + sqrt(two))*del_th + (-two-three/sqrt(two))*q
68    
69  c theta4  C theta4
70        del_th = del_P*sw_adtg(S,th,P+del_P)        del_th = del_P*sw_adtg(S,th,P+del_P)
71        SW_PTMP     = th + (del_th - two*q)/(two*three)        SW_PTMP     = th + (del_th - two*q)/(two*three)
       return  
       end  
72    
73  C======================================================================        RETURN
74          END
75    
76    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
77    
78  CBOP  CBOP
79  C     !ROUTINE: SW_TEMP  C     !ROUTINE: SW_TEMP
80  C     !INTERFACE:  C     !INTERFACE:
81        _RL FUNCTION SW_TEMP( s, t, p, pr )        _RL FUNCTION SW_TEMP( S, T, P, PR )
82  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
83  C     *=============================================================*  C     *=============================================================*
84  C     | S/R  SW_TEMP  C     | S/R  SW_TEMP
# Line 301  C     Bryden, H. 1973. Line 95  C     Bryden, H. 1973.
95  C     "New Polynomials for thermal expansion, adiabatic temperature gradient  C     "New Polynomials for thermal expansion, adiabatic temperature gradient
96  C     and potential temperature of sea water."  C     and potential temperature of sea water."
97  C     DEEP-SEA RES., 1973, Vol20,401-408.  C     DEEP-SEA RES., 1973, Vol20,401-408.
 C  
98    
99  C     !USES:  C     !USES:
100        IMPLICIT NONE        IMPLICIT NONE
   
101  C     === Global variables ===  C     === Global variables ===
 CML#include "SIZE.h"  
 CML#include "EEPARAMS.h"  
 CML#include "PARAMS.h"  
 CML#include "GRID.h"  
 CML#include "DYNVARS.h"  
 CML#include "FFIELDS.h"  
 CML#include "SHELFICE.h"  
102    
103  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
104  C     === Routine arguments ===  C     === Routine arguments ===
105  C     s      :: salinity  C     S      :: salinity
106  C     t      :: potential temperature  C     T      :: potential temperature
107  C     p      :: pressure  C     P      :: pressure
108  c     pr     :: reference pressure  C     PR     :: reference pressure
109  C     myIter :: iteration counter for this thread        _RL  S, T, P, PR
 C     myTime :: time counter for this thread  
 C     myThid :: thread number for this instance of the routine.  
       _RL  s, t, p, pr  
       _RL  myTime  
       INTEGER myIter  
       INTEGER myThid  
110  CEOP  CEOP
111    
112  C     !LOCAL VARIABLES  C     !FUNCTIONS:
113  C     === Local variables ===        _RL sw_adtg
114          EXTERNAL sw_adtg
115    
116    C     !LOCAL VARIABLES:
117        _RL del_P ,del_th, th, q        _RL del_P ,del_th, th, q
118        _RL onehalf, two, three        _RL onehalf, two, three
119        PARAMETER ( onehalf = 0.5 _d 0, two = 2. _d 0, three = 3. _d 0 )        PARAMETER ( onehalf = 0.5 _d 0, two = 2. _d 0, three = 3. _d 0 )
120    
121  c     externals  C theta1
       _RL sw_adtg  
       EXTERNAL sw_adtg  
 c theta1  
122  C--   here we swap P and PR in order to get in-situ temperature  C--   here we swap P and PR in order to get in-situ temperature
123  C     del_P  = PR - P ! to get potential from in-situ temperature  C     del_P  = PR - P ! to get potential from in-situ temperature
124        del_P  = P - PR ! to get in-situ from potential temperature        del_P  = P - PR ! to get in-situ from potential temperature
125        del_th = del_P*sw_adtg(S,T,P)        del_th = del_P*sw_adtg(S,T,P)
126        th     = T + onehalf*del_th        th     = T + onehalf*del_th
127        q      = del_th        q      = del_th
128  c theta2  C theta2
129        del_th = del_P*sw_adtg(S,th,P+onehalf*del_P)        del_th = del_P*sw_adtg(S,th,P+onehalf*del_P)
130    
131        th     = th + (1 - 1/sqrt(two))*(del_th - q)        th     = th + (1 - 1/sqrt(two))*(del_th - q)
132        q      = (two-sqrt(two))*del_th + (-two+three/sqrt(two))*q        q      = (two-sqrt(two))*del_th + (-two+three/sqrt(two))*q
133    
134  c theta3  C theta3
135        del_th = del_P*sw_adtg(S,th,P+onehalf*del_P)        del_th = del_P*sw_adtg(S,th,P+onehalf*del_P)
136        th     = th + (1 + 1/sqrt(two))*(del_th - q)        th     = th + (1 + 1/sqrt(two))*(del_th - q)
137        q      = (two + sqrt(two))*del_th + (-two-three/sqrt(two))*q        q      = (two + sqrt(two))*del_th + (-two-three/sqrt(two))*q
138    
139  c theta4  C theta4
140        del_th = del_P*sw_adtg(S,th,P+del_P)        del_th = del_P*sw_adtg(S,th,P+del_P)
141        SW_temp= th + (del_th - two*q)/(two*three)        SW_temp= th + (del_th - two*q)/(two*three)
142    
143        RETURN        RETURN
144        END        END
145    
146  C======================================================================  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
147    
148    CBOP
149    C     !ROUTINE: SW_ADTG
150    C     !INTERFACE:
151        _RL FUNCTION SW_ADTG  (S,T,P)        _RL FUNCTION SW_ADTG  (S,T,P)
152    
153  c     ==================================================================  C     !DESCRIPTION: \bv
154  c     SUBROUTINE SW_ADTG  C     *=============================================================*
155  c     ==================================================================  C     | S/R  SW_ADTG
156  c  C     | o compute adiabatic temperature gradient as per UNESCO 1983 routines.
157  c     o Calculates adiabatic temperature gradient as per UNESCO 1983 routines.  C     *=============================================================*
158  c  C
159  c     started:  C     started:
160  c  C              Armin Koehl akoehl@ucsd.edu
 c              Armin Koehl akoehl@ucsd.edu  
 c  
 c     ==================================================================  
 c     SUBROUTINE SW_ADTG  
 c     ==================================================================  
161    
162        implicit none  C     !USES:
163        _RL a0,a1,a2,a3,b0,b1,c0,c1,c2,c3,d0,d1,e0,e1,e2        IMPLICIT NONE
164    
165    C     !INPUT/OUTPUT PARAMETERS:
166        _RL S,T,P        _RL S,T,P
167    
168    C     !LOCAL VARIABLES:
169          _RL a0,a1,a2,a3,b0,b1,c0,c1,c2,c3,d0,d1,e0,e1,e2
170        _RL sref        _RL sref
171    CEOP
172    
173        sref = 35. _d 0        sref = 35. _d 0
174        a0 =  3.5803 _d -5        a0 =  3.5803 _d -5
# Line 412  c     ================================== Line 195  c     ==================================
195       &     + (b0 + b1*T)*(S-sref)       &     + (b0 + b1*T)*(S-sref)
196       &     + ( (c0 + (c1 + (c2 + c3*T)*T)*T) + (d0 + d1*T)*(S-sref) )*P       &     + ( (c0 + (c1 + (c2 + c3*T)*T)*T) + (d0 + d1*T)*(S-sref) )*P
197       &     + (  e0 + (e1 + e2*T)*T )*P*P       &     + (  e0 + (e1 + e2*T)*T )*P*P
198        return  
199        end        RETURN
200          END
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.5
  ViewVC Help
Powered by ViewVC 1.1.22