OpenAD/code_common/seawater.F

C$Header: /u/gcmpack/MITgcm_contrib/heimbach/OpenAD/code_common/seawater.F,v 1.1 2008/10/08 20:44:39 utke Exp $
C$Name:  $

#include "CPP_OPTIONS.h"

C--  File seawater.F: routines that compute quantities related to seawater.
C--   Contents
C--   o SW_PTMP: function to compute potential temperature
C--   o SW_TEMP: function to compute in-situ temperature from pot. temp.
C--   o SW_ADTG: function to compute adiabatic temperature gradient
C--              (used by both SW_PTMP & SW_TEMP)

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

CBOP
C     !ROUTINE: SW_PTMP
C     !INTERFACE:
      subroutine SW_PTMP  (S,T,P,PR, rv)

C     !DESCRIPTION: \bv
C     *=============================================================*
C     | S/R  SW_PTMP
C     | o compute potential temperature as per UNESCO 1983 report.
C     *=============================================================*
C
C     started:
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]

C     !USES:
      IMPLICIT NONE

C     !INPUT/OUTPUT PARAMETERS:
      _RL S,T,P,PR
      _RL rv

C     !LOCAL VARIABLES
      _RL del_P ,del_th, th, q
      _RL onehalf, two, three
      PARAMETER ( onehalf = 0.5 _d 0, two = 2. _d 0, three = 3. _d 0 )
      _RL adtg_val
CEOP

C theta1
      del_P  = PR - P
      call sw_adtg(S,T,P, adtg_val)
      del_th = del_P*adtg_val
      th     = T + onehalf*del_th
      q      = del_th
C theta2
      call sw_adtg(S,th,P+onehalf*del_P, adtg_val)
      del_th = del_P*adtg_val

      th     = th + (1 - 1/sqrt(two))*(del_th - q)
      q      = (two-sqrt(two))*del_th + (-two+three/sqrt(two))*q

C theta3
      call sw_adtg(S,th,P+onehalf*del_P, adtg_val)
      del_th = del_P*adtg_val
      th     = th + (1 + 1/sqrt(two))*(del_th - q)
      q      = (two + sqrt(two))*del_th + (-two-three/sqrt(two))*q

C theta4
      call sw_adtg(S,th,P+del_P, adtg_val)
      del_th = del_P*adtg_val
      rv     = th + (del_th - two*q)/(two*three)
      RETURN
      END

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

CBOP
C     !ROUTINE: SW_TEMP
C     !INTERFACE:
      SUBROUTINE SW_TEMP( s, t, p, pr, rv)
C     !DESCRIPTION: \bv
C     *=============================================================*
C     | S/R  SW_TEMP
C     | o compute in-situ temperature from potential temperature
C     *=============================================================*
C
C     REFERENCES:
C     Fofonoff, P. and Millard, R.C. Jr
C     Unesco 1983. Algorithms for computation of fundamental properties of
C     seawater, 1983. _Unesco Tech. Pap. in Mar. Sci._, No. 44, 53 pp.
C     Eqn.(31) p.39
C
C     Bryden, H. 1973.
C     "New Polynomials for thermal expansion, adiabatic temperature gradient
C     and potential temperature of sea water."
C     DEEP-SEA RES., 1973, Vol20,401-408.

C     !USES:
      IMPLICIT NONE
C     === Global variables ===

C     !INPUT/OUTPUT PARAMETERS:
C     === Routine arguments ===
C     S      :: salinity
C     T      :: potential temperature
C     P      :: pressure
C     PR     :: reference pressure
      _RL  S, T, P, PR
      _RL rv
CEOP


C     !LOCAL VARIABLES:
      _RL del_P ,del_th, th, q
      _RL onehalf, two, three
      PARAMETER ( onehalf = 0.5 _d 0, two = 2. _d 0, three = 3. _d 0 )
      _RL adtg_val

C theta1
C--   here we swap P and PR in order to get in-situ temperature
C     del_P  = PR - P ! to get potential from in-situ temperature
      del_P  = P - PR ! to get in-situ from potential temperature
      call sw_adtg(S,T,P, adtg_val)
      del_th = del_P*adtg_val
      th     = T + onehalf*del_th
      q      = del_th
C theta2
      call sw_adtg(S,th,P+onehalf*del_P, adtg_val)
      del_th = del_P*adtg_val

      th     = th + (1 - 1/sqrt(two))*(del_th - q)
      q      = (two-sqrt(two))*del_th + (-two+three/sqrt(two))*q

C theta3
      call sw_adtg(S,th,P+onehalf*del_P, adtg_val)
      del_th = del_P*adtg_val
      th     = th + (1 + 1/sqrt(two))*(del_th - q)
      q      = (two + sqrt(two))*del_th + (-two-three/sqrt(two))*q

C theta4
      call sw_adtg(S,th,P+del_P, adtg_val)
      del_th = del_P*adtg_val
      rv = th + (del_th - two*q)/(two*three)

      RETURN
      END

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

CBOP
C     !ROUTINE: SW_ADTG
C     !INTERFACE:
      SUBROUTINE SW_ADTG  (S,T,P, rv)

C     !DESCRIPTION: \bv
C     *=============================================================*
C     | S/R  SW_ADTG
C     | o compute adiabatic temperature gradient as per UNESCO 1983 routines.
C     *=============================================================*
C
C     started:
C              Armin Koehl akoehl@ucsd.edu

C     !USES:
      IMPLICIT NONE

C     !INPUT/OUTPUT PARAMETERS:
      _RL S,T,P

C     !LOCAL VARIABLES:
      _RL a0,a1,a2,a3,b0,b1,c0,c1,c2,c3,d0,d1,e0,e1,e2
      _RL sref
      _RL rv
CEOP

      sref = 35. _d 0
      a0 =  3.5803 _d -5
      a1 = +8.5258 _d -6
      a2 = -6.836 _d -8
      a3 =  6.6228 _d -10

      b0 = +1.8932 _d -6
      b1 = -4.2393 _d -8

      c0 = +1.8741 _d -8
      c1 = -6.7795 _d -10
      c2 = +8.733 _d -12
      c3 = -5.4481 _d -14

      d0 = -1.1351 _d -10
      d1 =  2.7759 _d -12

      e0 = -4.6206 _d -13
      e1 = +1.8676 _d -14
      e2 = -2.1687 _d -16

      rv =      a0 + (a1 + (a2 + a3*T)*T)*T
     &     + (b0 + b1*T)*(S-sref)
     &     + ( (c0 + (c1 + (c2 + c3*T)*T)*T) + (d0 + d1*T)*(S-sref) )*P
     &     + (  e0 + (e1 + e2*T)*T )*P*P

      RETURN
      END
1	utke	1.2	C$Header: /u/gcmpack/MITgcm_contrib/heimbach/OpenAD/code_common/seawater.F,v 1.1 2008/10/08 20:44:39 utke Exp $
2	utke	1.1	C$Name: $
3
4			#include "CPP_OPTIONS.h"
5
6			C-- File seawater.F: routines that compute quantities related to seawater.
7			C-- Contents
8			C-- o SW_PTMP: function to compute potential temperature
9	utke	1.2	C-- o SW_TEMP: function to compute in-situ temperature from pot. temp.
10			C-- o SW_ADTG: function to compute adiabatic temperature gradient
11			C-- (used by both SW_PTMP & SW_TEMP)
12
13			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
14
15			CBOP
16			C !ROUTINE: SW_PTMP
17			C !INTERFACE:
18			subroutine SW_PTMP (S,T,P,PR, rv)
19	utke	1.1
20			C !DESCRIPTION: \bv
21	utke	1.2	C =============================================================
22			C \| S/R SW_PTMP
23			C \| o compute potential temperature as per UNESCO 1983 report.
24			C =============================================================
25			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	utke	1.1
37			C !USES:
38			IMPLICIT NONE
39
40			C !INPUT/OUTPUT PARAMETERS:
41			_RL S,T,P,PR
42			_RL rv
43
44	utke	1.2	C !LOCAL VARIABLES
45	utke	1.1	_RL del_P ,del_th, th, q
46			_RL onehalf, two, three
47	utke	1.2	PARAMETER ( onehalf = 0.5 _d 0, two = 2. _d 0, three = 3. _d 0 )
48			_RL adtg_val
49			CEOP
50	utke	1.1
51	utke	1.2	C theta1
52	utke	1.1	del_P = PR - P
53			call sw_adtg(S,T,P, adtg_val)
54			del_th = del_P*adtg_val
55			th = T + onehalf*del_th
56			q = del_th
57	utke	1.2	C theta2
58	utke	1.1	call sw_adtg(S,th,P+onehalf*del_P, adtg_val)
59			del_th = del_P*adtg_val
60
61			th = th + (1 - 1/sqrt(two))*(del_th - q)
62			q = (two-sqrt(two))del_th + (-two+three/sqrt(two))q
63
64	utke	1.2	C theta3
65	utke	1.1	call sw_adtg(S,th,P+onehalf*del_P, adtg_val)
66			del_th = del_P*adtg_val
67			th = th + (1 + 1/sqrt(two))*(del_th - q)
68			q = (two + sqrt(two))del_th + (-two-three/sqrt(two))q
69
70	utke	1.2	C theta4
71	utke	1.1	call sw_adtg(S,th,P+del_P, adtg_val)
72			del_th = del_P*adtg_val
73			rv = th + (del_th - twoq)/(twothree)
74	utke	1.2	RETURN
75			END
76	utke	1.1
77	utke	1.2	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
78	utke	1.1
79			CBOP
80			C !ROUTINE: SW_TEMP
81			C !INTERFACE:
82			SUBROUTINE SW_TEMP( s, t, p, pr, rv)
83			C !DESCRIPTION: \bv
84			C =============================================================
85			C \| S/R SW_TEMP
86			C \| o compute in-situ temperature from potential temperature
87			C =============================================================
88			C
89			C REFERENCES:
90			C Fofonoff, P. and Millard, R.C. Jr
91			C Unesco 1983. Algorithms for computation of fundamental properties of
92			C seawater, 1983. _Unesco Tech. Pap. in Mar. Sci._, No. 44, 53 pp.
93			C Eqn.(31) p.39
94			C
95			C Bryden, H. 1973.
96			C "New Polynomials for thermal expansion, adiabatic temperature gradient
97			C and potential temperature of sea water."
98			C DEEP-SEA RES., 1973, Vol20,401-408.
99
100			C !USES:
101			IMPLICIT NONE
102			C === Global variables ===
103
104			C !INPUT/OUTPUT PARAMETERS:
105			C === Routine arguments ===
106	utke	1.2	C S :: salinity
107			C T :: potential temperature
108			C P :: pressure
109			C PR :: reference pressure
110			_RL S, T, P, PR
111	utke	1.1	_RL rv
112			CEOP
113
114	utke	1.2
115			C !LOCAL VARIABLES:
116	utke	1.1	_RL del_P ,del_th, th, q
117			_RL onehalf, two, three
118			PARAMETER ( onehalf = 0.5 _d 0, two = 2. _d 0, three = 3. _d 0 )
119	utke	1.2	_RL adtg_val
120	utke	1.1
121	utke	1.2	C theta1
122	utke	1.1	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
124			del_P = P - PR ! to get in-situ from potential temperature
125			call sw_adtg(S,T,P, adtg_val)
126			del_th = del_P*adtg_val
127			th = T + onehalf*del_th
128			q = del_th
129	utke	1.2	C theta2
130	utke	1.1	call sw_adtg(S,th,P+onehalf*del_P, adtg_val)
131			del_th = del_P*adtg_val
132
133			th = th + (1 - 1/sqrt(two))*(del_th - q)
134			q = (two-sqrt(two))del_th + (-two+three/sqrt(two))q
135
136	utke	1.2	C theta3
137	utke	1.1	call sw_adtg(S,th,P+onehalf*del_P, adtg_val)
138			del_th = del_P*adtg_val
139			th = th + (1 + 1/sqrt(two))*(del_th - q)
140			q = (two + sqrt(two))del_th + (-two-three/sqrt(two))q
141
142	utke	1.2	C theta4
143	utke	1.1	call sw_adtg(S,th,P+del_P, adtg_val)
144			del_th = del_P*adtg_val
145			rv = th + (del_th - twoq)/(twothree)
146
147			RETURN
148			END
149
150	utke	1.2	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
151	utke	1.1
152	utke	1.2	CBOP
153			C !ROUTINE: SW_ADTG
154			C !INTERFACE:
155	utke	1.1	SUBROUTINE SW_ADTG (S,T,P, rv)
156
157	utke	1.2	C !DESCRIPTION: \bv
158			C =============================================================
159			C \| S/R SW_ADTG
160			C \| o compute adiabatic temperature gradient as per UNESCO 1983 routines.
161			C =============================================================
162			C
163			C started:
164			C Armin Koehl akoehl@ucsd.edu
165
166			C !USES:
167			IMPLICIT NONE
168	utke	1.1
169	utke	1.2	C !INPUT/OUTPUT PARAMETERS:
170			_RL S,T,P
171
172			C !LOCAL VARIABLES:
173	utke	1.1	_RL a0,a1,a2,a3,b0,b1,c0,c1,c2,c3,d0,d1,e0,e1,e2
174			_RL sref
175			_RL rv
176	utke	1.2	CEOP
177	utke	1.1
178			sref = 35. _d 0
179			a0 = 3.5803 _d -5
180			a1 = +8.5258 _d -6
181			a2 = -6.836 _d -8
182			a3 = 6.6228 _d -10
183
184			b0 = +1.8932 _d -6
185			b1 = -4.2393 _d -8
186
187			c0 = +1.8741 _d -8
188			c1 = -6.7795 _d -10
189			c2 = +8.733 _d -12
190			c3 = -5.4481 _d -14
191
192			d0 = -1.1351 _d -10
193			d1 = 2.7759 _d -12
194
195			e0 = -4.6206 _d -13
196			e1 = +1.8676 _d -14
197			e2 = -2.1687 _d -16
198
199			rv = a0 + (a1 + (a2 + a3T)T)*T
200			& + (b0 + b1T)(S-sref)
201			& + ( (c0 + (c1 + (c2 + c3T)T)T) + (d0 + d1T)(S-sref) )P
202			& + ( e0 + (e1 + e2T)T )PP
203	utke	1.2
204			RETURN
205			END