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

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

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


Revision 1.2 - (hide annotations) (download)
Sat Aug 9 00:51:34 2008 UTC (15 years, 9 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint61c
Changes since 1.1: +47 -45 lines
list content (ALL s/r & functions) at the top of the file

1 jmc 1.2 C$Header: /u/gcmpack/MITgcm/model/src/seawater.F,v 1.1 2008/02/27 19:42:02 mlosch Exp $
2     C$Name: $
3 mlosch 1.1
4     #include "CPP_OPTIONS.h"
5 jmc 1.2
6     C-- File seawater.F: routines that compute quantities related to seawater.
7     C-- Contents
8     C-- o FIND_RHO_SCALAR: in-situ density for individual points
9     C-- o SW_PTMP: function to compute potential temperature
10     C-- o SW_TEMP: function to compute potential temperature
11     C-- o SW_ADTG: function to compute adiabatic tmperature gradient
12     C-- used by sw_ptmp
13    
14     SUBROUTINE FIND_RHO_SCALAR(
15 mlosch 1.1 I tLoc, sLoc, pLoc,
16     O rhoLoc,
17     I myThid )
18    
19     C !DESCRIPTION: \bv
20     C *==========================================================*
21 jmc 1.2 C | o SUBROUTINE FIND_RHO_SCALAR
22     C | Calculates [rho(S,T,p)-rhoConst]
23 mlosch 1.1 C *==========================================================*
24     C \ev
25    
26     C !USES:
27     IMPLICIT NONE
28     C == Global variables ==
29     #include "SIZE.h"
30     #include "EEPARAMS.h"
31     #include "PARAMS.h"
32     #include "EOS.h"
33    
34     C !INPUT/OUTPUT PARAMETERS:
35     C == Routine arguments ==
36     _RL sLoc, tLoc, pLoc
37     _RL rhoLoc
38     INTEGER myThid
39    
40     C !LOCAL VARIABLES:
41     C == Local variables ==
42    
43     _RL t1, t2, t3, t4, s1, s3o2, p1, p2, sp5, p1t1
44     _RL rfresh, rsalt, rhoP0
45     _RL bMfresh, bMsalt, bMpres, BulkMod
46     _RL rhoNum, rhoDen, den, epsln
47     parameter ( epsln = 0. _d 0 )
48    
49     CHARACTER*(MAX_LEN_MBUF) msgBuf
50     CEOP
51    
52     rhoLoc = 0. _d 0
53     rhoP0 = 0. _d 0
54     bulkMod = 0. _d 0
55     rfresh = 0. _d 0
56     rsalt = 0. _d 0
57     bMfresh = 0. _d 0
58     bMsalt = 0. _d 0
59     bMpres = 0. _d 0
60     rhoNum = 0. _d 0
61     rhoDen = 0. _d 0
62     den = 0. _d 0
63    
64     t1 = tLoc
65     t2 = t1*t1
66     t3 = t2*t1
67     t4 = t3*t1
68 jmc 1.2
69 mlosch 1.1 s1 = sLoc
70     IF ( s1 .LT. 0. _d 0 ) THEN
71     C issue a warning
72 jmc 1.2 WRITE(msgBuf,'(A,E13.5)')
73 mlosch 1.1 & ' FIND_RHO_SCALAR: WARNING, salinity = ', s1
74     CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
75     & SQUEEZE_RIGHT , myThid )
76     s1 = 0. _d 0
77     ENDIF
78    
79     IF (equationOfState.EQ.'LINEAR') THEN
80    
81     rholoc = rhoNil*(
82     & sBeta *(sLoc-sRef(1))
83     & -tAlpha*(tLoc-tRef(1))
84     & ) + (rhoNil-rhoConst)
85     c rhoLoc = 0. _d 0
86    
87     ELSEIF (equationOfState.EQ.'POLY3') THEN
88    
89     C this is not correct, there is a field eosSig0 which should be use here
90     C but I DO not intent to include the reference level in this routine
91 jmc 1.2 WRITE(msgBuf,'(A)')
92 mlosch 1.1 & ' FIND_RHO_SCALAR: for POLY3, the density is not'
93     CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
94     & SQUEEZE_RIGHT , myThid )
95 jmc 1.2 WRITE(msgBuf,'(A)')
96 mlosch 1.1 & ' computed correctly in this routine'
97     CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
98     & SQUEEZE_RIGHT , myThid )
99     rhoLoc = 0. _d 0
100    
101     ELSEIF ( equationOfState(1:5).EQ.'JMD95'
102     & .OR. equationOfState.EQ.'UNESCO' ) THEN
103     C nonlinear equation of state in pressure coordinates
104    
105     s3o2 = s1*SQRT(s1)
106 jmc 1.2
107 mlosch 1.1 p1 = pLoc*SItoBar
108     p2 = p1*p1
109    
110     C density of freshwater at the surface
111 jmc 1.2 rfresh =
112     & eosJMDCFw(1)
113 mlosch 1.1 & + eosJMDCFw(2)*t1
114     & + eosJMDCFw(3)*t2
115     & + eosJMDCFw(4)*t3
116     & + eosJMDCFw(5)*t4
117     & + eosJMDCFw(6)*t4*t1
118     C density of sea water at the surface
119 jmc 1.2 rsalt =
120 mlosch 1.1 & s1*(
121     & eosJMDCSw(1)
122     & + eosJMDCSw(2)*t1
123     & + eosJMDCSw(3)*t2
124     & + eosJMDCSw(4)*t3
125     & + eosJMDCSw(5)*t4
126     & )
127     & + s3o2*(
128     & eosJMDCSw(6)
129     & + eosJMDCSw(7)*t1
130     & + eosJMDCSw(8)*t2
131     & )
132     & + eosJMDCSw(9)*s1*s1
133    
134     rhoP0 = rfresh + rsalt
135    
136     C secant bulk modulus of fresh water at the surface
137 jmc 1.2 bMfresh =
138 mlosch 1.1 & eosJMDCKFw(1)
139     & + eosJMDCKFw(2)*t1
140     & + eosJMDCKFw(3)*t2
141     & + eosJMDCKFw(4)*t3
142     & + eosJMDCKFw(5)*t4
143     C secant bulk modulus of sea water at the surface
144     bMsalt =
145     & s1*( eosJMDCKSw(1)
146     & + eosJMDCKSw(2)*t1
147     & + eosJMDCKSw(3)*t2
148     & + eosJMDCKSw(4)*t3
149     & )
150     & + s3o2*( eosJMDCKSw(5)
151     & + eosJMDCKSw(6)*t1
152     & + eosJMDCKSw(7)*t2
153     & )
154     C secant bulk modulus of sea water at pressure p
155 jmc 1.2 bMpres =
156 mlosch 1.1 & p1*( eosJMDCKP(1)
157     & + eosJMDCKP(2)*t1
158     & + eosJMDCKP(3)*t2
159     & + eosJMDCKP(4)*t3
160     & )
161     & + p1*s1*( eosJMDCKP(5)
162     & + eosJMDCKP(6)*t1
163     & + eosJMDCKP(7)*t2
164     & )
165     & + p1*s3o2*eosJMDCKP(8)
166     & + p2*( eosJMDCKP(9)
167     & + eosJMDCKP(10)*t1
168     & + eosJMDCKP(11)*t2
169     & )
170     & + p2*s1*( eosJMDCKP(12)
171     & + eosJMDCKP(13)*t1
172     & + eosJMDCKP(14)*t2
173     & )
174    
175     bulkMod = bMfresh + bMsalt + bMpres
176 jmc 1.2
177 mlosch 1.1 C density of sea water at pressure p
178     rhoLoc = rhoP0/(1. _d 0 - p1/bulkMod) - rhoConst
179    
180     ELSEIF ( equationOfState.EQ.'MDJWF' ) THEN
181    
182 jmc 1.2 sp5 = SQRT(s1)
183    
184 mlosch 1.1 p1 = pLoc*SItodBar
185     p1t1 = p1*t1
186    
187     rhoNum = eosMDJWFnum(0)
188 jmc 1.2 & + t1*(eosMDJWFnum(1)
189     & + t1*(eosMDJWFnum(2) + eosMDJWFnum(3)*t1) )
190 mlosch 1.1 & + s1*(eosMDJWFnum(4)
191 jmc 1.2 & + eosMDJWFnum(5)*t1 + eosMDJWFnum(6)*s1)
192     & + p1*(eosMDJWFnum(7) + eosMDJWFnum(8)*t2
193     & + eosMDJWFnum(9)*s1
194 mlosch 1.1 & + p1*(eosMDJWFnum(10) + eosMDJWFnum(11)*t2) )
195    
196 jmc 1.2
197 mlosch 1.1 den = eosMDJWFden(0)
198     & + t1*(eosMDJWFden(1)
199     & + t1*(eosMDJWFden(2)
200     & + t1*(eosMDJWFden(3) + t1*eosMDJWFden(4) ) ) )
201     & + s1*(eosMDJWFden(5)
202 jmc 1.2 & + t1*(eosMDJWFden(6)
203     & + eosMDJWFden(7)*t2)
204     & + sp5*(eosMDJWFden(8) + eosMDJWFden(9)*t2) )
205 mlosch 1.1 & + p1*(eosMDJWFden(10)
206     & + p1t1*(eosMDJWFden(11)*t2 + eosMDJWFden(12)*p1) )
207 jmc 1.2
208     rhoDen = 1.0/(epsln+den)
209 mlosch 1.1
210     rhoLoc = rhoNum*rhoDen - rhoConst
211    
212     ELSEIF( equationOfState .EQ. 'IDEALG' ) THEN
213 jmc 1.2 C
214 mlosch 1.1 ELSE
215     WRITE(msgBuf,'(3A)')
216     & ' FIND_RHO_SCALAR : equationOfState = "',
217     & equationOfState,'"'
218     CALL PRINT_ERROR( msgBuf, myThid )
219     STOP 'ABNORMAL END: S/R FIND_RHO_SCALAR'
220     ENDIF
221    
222 jmc 1.2 RETURN
223 mlosch 1.1 END
224    
225     C=================================================================
226    
227 jmc 1.2 _RL FUNCTION SW_PTMP (S,T,P,PR)
228 mlosch 1.1
229     c ==================================================================
230     c SUBROUTINE SW_PTMP
231     c ==================================================================
232     c
233     c o Calculates potential temperature as per UNESCO 1983 report.
234     c
235     c started:
236     c
237     c Armin Koehl akoehl@ucsd.edu
238     c
239     c ==================================================================
240     c SUBROUTINE SW_PTMP
241     c ==================================================================
242     C S = salinity [psu (PSS-78) ]
243     C T = temperature [degree C (IPTS-68)]
244     C P = pressure [db]
245     C PR = Reference pressure [db]
246    
247     implicit none
248    
249     c routine arguments
250     _RL S,T,P,PR
251    
252     c local arguments
253     _RL del_P ,del_th, th, q
254     _RL onehalf, two, three
255     parameter ( onehalf = 0.5 _d 0, two = 2. _d 0, three = 3. _d 0 )
256    
257     c externals
258     _RL sw_adtg
259     external sw_adtg
260     c theta1
261     del_P = PR - P
262     del_th = del_P*sw_adtg(S,T,P)
263     th = T + onehalf*del_th
264     q = del_th
265     c theta2
266     del_th = del_P*sw_adtg(S,th,P+onehalf*del_P)
267    
268     th = th + (1 - 1/sqrt(two))*(del_th - q)
269     q = (two-sqrt(two))*del_th + (-two+three/sqrt(two))*q
270    
271     c theta3
272     del_th = del_P*sw_adtg(S,th,P+onehalf*del_P)
273     th = th + (1 + 1/sqrt(two))*(del_th - q)
274     q = (two + sqrt(two))*del_th + (-two-three/sqrt(two))*q
275    
276     c theta4
277     del_th = del_P*sw_adtg(S,th,P+del_P)
278     SW_PTMP = th + (del_th - two*q)/(two*three)
279     return
280     end
281    
282     C======================================================================
283    
284     CBOP
285     C !ROUTINE: SW_TEMP
286     C !INTERFACE:
287     _RL FUNCTION SW_TEMP( s, t, p, pr )
288     C !DESCRIPTION: \bv
289     C *=============================================================*
290     C | S/R SW_TEMP
291 jmc 1.2 C | o compute in-situ temperature from potential temperature
292 mlosch 1.1 C *=============================================================*
293     C
294     C REFERENCES:
295     C Fofonoff, P. and Millard, R.C. Jr
296     C Unesco 1983. Algorithms for computation of fundamental properties of
297     C seawater, 1983. _Unesco Tech. Pap. in Mar. Sci._, No. 44, 53 pp.
298     C Eqn.(31) p.39
299 jmc 1.2 C
300 mlosch 1.1 C Bryden, H. 1973.
301     C "New Polynomials for thermal expansion, adiabatic temperature gradient
302     C and potential temperature of sea water."
303     C DEEP-SEA RES., 1973, Vol20,401-408.
304     C
305    
306     C !USES:
307     IMPLICIT NONE
308    
309     C === Global variables ===
310     CML#include "SIZE.h"
311     CML#include "EEPARAMS.h"
312     CML#include "PARAMS.h"
313     CML#include "GRID.h"
314     CML#include "DYNVARS.h"
315     CML#include "FFIELDS.h"
316     CML#include "SHELFICE.h"
317 jmc 1.2
318 mlosch 1.1 C !INPUT/OUTPUT PARAMETERS:
319     C === Routine arguments ===
320     C s :: salinity
321     C t :: potential temperature
322     C p :: pressure
323     c pr :: reference pressure
324     C myIter :: iteration counter for this thread
325     C myTime :: time counter for this thread
326     C myThid :: thread number for this instance of the routine.
327     _RL s, t, p, pr
328     _RL myTime
329     INTEGER myIter
330     INTEGER myThid
331     CEOP
332    
333 jmc 1.2 C !LOCAL VARIABLES
334 mlosch 1.1 C === Local variables ===
335     _RL del_P ,del_th, th, q
336     _RL onehalf, two, three
337     PARAMETER ( onehalf = 0.5 _d 0, two = 2. _d 0, three = 3. _d 0 )
338    
339     c externals
340     _RL sw_adtg
341     EXTERNAL sw_adtg
342     c theta1
343     C-- here we swap P and PR in order to get in-situ temperature
344     C del_P = PR - P ! to get potential from in-situ temperature
345     del_P = P - PR ! to get in-situ from potential temperature
346     del_th = del_P*sw_adtg(S,T,P)
347     th = T + onehalf*del_th
348     q = del_th
349     c theta2
350     del_th = del_P*sw_adtg(S,th,P+onehalf*del_P)
351    
352     th = th + (1 - 1/sqrt(two))*(del_th - q)
353     q = (two-sqrt(two))*del_th + (-two+three/sqrt(two))*q
354    
355     c theta3
356     del_th = del_P*sw_adtg(S,th,P+onehalf*del_P)
357     th = th + (1 + 1/sqrt(two))*(del_th - q)
358     q = (two + sqrt(two))*del_th + (-two-three/sqrt(two))*q
359    
360     c theta4
361     del_th = del_P*sw_adtg(S,th,P+del_P)
362     SW_temp= th + (del_th - two*q)/(two*three)
363    
364     RETURN
365     END
366    
367     C======================================================================
368    
369 jmc 1.2 _RL FUNCTION SW_ADTG (S,T,P)
370 mlosch 1.1
371     c ==================================================================
372     c SUBROUTINE SW_ADTG
373     c ==================================================================
374     c
375     c o Calculates adiabatic temperature gradient as per UNESCO 1983 routines.
376     c
377     c started:
378     c
379     c Armin Koehl akoehl@ucsd.edu
380     c
381     c ==================================================================
382     c SUBROUTINE SW_ADTG
383     c ==================================================================
384    
385     implicit none
386     _RL a0,a1,a2,a3,b0,b1,c0,c1,c2,c3,d0,d1,e0,e1,e2
387     _RL S,T,P
388     _RL sref
389    
390     sref = 35. _d 0
391     a0 = 3.5803 _d -5
392     a1 = +8.5258 _d -6
393     a2 = -6.836 _d -8
394     a3 = 6.6228 _d -10
395    
396     b0 = +1.8932 _d -6
397     b1 = -4.2393 _d -8
398    
399     c0 = +1.8741 _d -8
400     c1 = -6.7795 _d -10
401     c2 = +8.733 _d -12
402     c3 = -5.4481 _d -14
403    
404     d0 = -1.1351 _d -10
405     d1 = 2.7759 _d -12
406    
407     e0 = -4.6206 _d -13
408     e1 = +1.8676 _d -14
409     e2 = -2.1687 _d -16
410    
411     SW_ADTG = a0 + (a1 + (a2 + a3*T)*T)*T
412     & + (b0 + b1*T)*(S-sref)
413     & + ( (c0 + (c1 + (c2 + c3*T)*T)*T) + (d0 + d1*T)*(S-sref) )*P
414     & + ( e0 + (e1 + e2*T)*T )*P*P
415 jmc 1.2 return
416 mlosch 1.1 end

  ViewVC Help
Powered by ViewVC 1.1.22