/[MITgcm]/MITgcm_contrib/heimbach/OpenAD/code_ad_moc/seawater.F
ViewVC logotype

Annotation of /MITgcm_contrib/heimbach/OpenAD/code_ad_moc/seawater.F

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


Revision 1.2 - (hide annotations) (download)
Tue Nov 18 17:33:11 2008 UTC (16 years, 8 months ago) by utke
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +1 -1 lines
FILE REMOVED
no longer needed here

1 utke 1.2 C$Header: /u/gcmpack/MITgcm_contrib/heimbach/OpenAD/code_ad_moc/seawater.F,v 1.1 2008/09/15 22:16:01 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 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     I tLoc, sLoc, pLoc,
16     O rhoLoc,
17     I myThid )
18    
19     C !DESCRIPTION: \bv
20     C *==========================================================*
21     C | o SUBROUTINE FIND_RHO_SCALAR
22     C | Calculates rho(S,T,p)
23     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    
69     s1 = sLoc
70     IF ( s1 .LT. 0. _d 0 ) THEN
71     C issue a warning
72     WRITE(msgBuf,'(A,E13.5)')
73     & ' 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
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     WRITE(msgBuf,'(A)')
92     & ' FIND_RHO_SCALAR: for POLY3, the density is not'
93     CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
94     & SQUEEZE_RIGHT , myThid )
95     WRITE(msgBuf,'(A)')
96     & ' 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    
107     p1 = pLoc*SItoBar
108     p2 = p1*p1
109    
110     C density of freshwater at the surface
111     rfresh =
112     & eosJMDCFw(1)
113     & + 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     rsalt =
120     & 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     bMfresh =
138     & 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     bMpres =
156     & 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    
177     C density of sea water at pressure p
178     rhoLoc = rhoP0/(1. _d 0 - p1/bulkMod)
179    
180     ELSEIF ( equationOfState.EQ.'MDJWF' ) THEN
181    
182     sp5 = SQRT(s1)
183    
184     p1 = pLoc*SItodBar
185     p1t1 = p1*t1
186    
187     rhoNum = eosMDJWFnum(0)
188     & + t1*(eosMDJWFnum(1)
189     & + t1*(eosMDJWFnum(2) + eosMDJWFnum(3)*t1) )
190     & + s1*(eosMDJWFnum(4)
191     & + eosMDJWFnum(5)*t1 + eosMDJWFnum(6)*s1)
192     & + p1*(eosMDJWFnum(7) + eosMDJWFnum(8)*t2
193     & + eosMDJWFnum(9)*s1
194     & + p1*(eosMDJWFnum(10) + eosMDJWFnum(11)*t2) )
195    
196    
197     den = eosMDJWFden(0)
198     & + t1*(eosMDJWFden(1)
199     & + t1*(eosMDJWFden(2)
200     & + t1*(eosMDJWFden(3) + t1*eosMDJWFden(4) ) ) )
201     & + s1*(eosMDJWFden(5)
202     & + t1*(eosMDJWFden(6)
203     & + eosMDJWFden(7)*t2)
204     & + sp5*(eosMDJWFden(8) + eosMDJWFden(9)*t2) )
205     & + p1*(eosMDJWFden(10)
206     & + p1t1*(eosMDJWFden(11)*t2 + eosMDJWFden(12)*p1) )
207    
208     rhoDen = 1.0/(epsln+den)
209    
210     rhoLoc = rhoNum*rhoDen
211    
212     ELSEIF( equationOfState .EQ. 'IDEALG' ) THEN
213     C
214     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     RETURN
223     END
224    
225     subroutine SW_PTMP (S,T,P,PR, rv)
226    
227     c ==================================================================
228     c SUBROUTINE SW_PTMP
229     c ==================================================================
230     c
231     c o Calculates potential temperature as per UNESCO 1983 report.
232     c
233     c started:
234     c
235     c Armin Koehl akoehl@ucsd.edu
236     c
237     c ==================================================================
238     c SUBROUTINE SW_PTMP
239     c ==================================================================
240     C S = salinity [psu (PSS-78) ]
241     C T = temperature [degree C (IPTS-68)]
242     C P = pressure [db]
243     C PR = Reference pressure [db]
244    
245     implicit none
246    
247     c routine arguments
248     _RL S,T,P,PR
249    
250     _RL rv
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 adtg_val
259     c theta1
260     del_P = PR - P
261     call sw_adtg(S,T,P, adtg_val)
262     del_th = del_P*adtg_val
263     th = T + onehalf*del_th
264     q = del_th
265     c theta2
266     call sw_adtg(S,th,P+onehalf*del_P, adtg_val)
267     del_th = del_P*adtg_val
268    
269     th = th + (1 - 1/sqrt(two))*(del_th - q)
270     q = (two-sqrt(two))*del_th + (-two+three/sqrt(two))*q
271    
272     c theta3
273     call sw_adtg(S,th,P+onehalf*del_P, adtg_val)
274     del_th = del_P*adtg_val
275     th = th + (1 + 1/sqrt(two))*(del_th - q)
276     q = (two + sqrt(two))*del_th + (-two-three/sqrt(two))*q
277    
278     c theta4
279     call sw_adtg(S,th,P+del_P, adtg_val)
280     del_th = del_P*adtg_val
281     rv = th + (del_th - two*q)/(two*three)
282     return
283     end
284    
285     C======================================================================
286    
287     CBOP
288     C !ROUTINE: SW_TEMP
289     C !INTERFACE:
290     SUBROUTINE SW_TEMP( s, t, p, pr, rv)
291     C !DESCRIPTION: \bv
292     C *=============================================================*
293     C | S/R SW_TEMP
294     C | o compute in-situ temperature from potential temperature
295     C *=============================================================*
296     C
297     C REFERENCES:
298     C Fofonoff, P. and Millard, R.C. Jr
299     C Unesco 1983. Algorithms for computation of fundamental properties of
300     C seawater, 1983. _Unesco Tech. Pap. in Mar. Sci._, No. 44, 53 pp.
301     C Eqn.(31) p.39
302     C
303     C Bryden, H. 1973.
304     C "New Polynomials for thermal expansion, adiabatic temperature gradient
305     C and potential temperature of sea water."
306     C DEEP-SEA RES., 1973, Vol20,401-408.
307     C
308    
309     C !USES:
310     IMPLICIT NONE
311    
312     C === Global variables ===
313     CML#include "SIZE.h"
314     CML#include "EEPARAMS.h"
315     CML#include "PARAMS.h"
316     CML#include "GRID.h"
317     CML#include "DYNVARS.h"
318     CML#include "FFIELDS.h"
319     CML#include "SHELFICE.h"
320    
321     C !INPUT/OUTPUT PARAMETERS:
322     C === Routine arguments ===
323     C s :: salinity
324     C t :: potential temperature
325     C p :: pressure
326     c pr :: reference pressure
327     C myIter :: iteration counter for this thread
328     C myTime :: time counter for this thread
329     C myThid :: thread number for this instance of the routine.
330     _RL s, t, p, pr
331     _RL myTime
332     INTEGER myIter
333     INTEGER myThid
334     _RL rv
335     CEOP
336    
337     C !LOCAL VARIABLES
338     C === Local variables ===
339     _RL del_P ,del_th, th, q
340     _RL onehalf, two, three
341     PARAMETER ( onehalf = 0.5 _d 0, two = 2. _d 0, three = 3. _d 0 )
342    
343     c externals
344     _RL adtg_val
345     c theta1
346     C-- here we swap P and PR in order to get in-situ temperature
347     C del_P = PR - P ! to get potential from in-situ temperature
348     del_P = P - PR ! to get in-situ from potential temperature
349     call sw_adtg(S,T,P, adtg_val)
350     del_th = del_P*adtg_val
351     th = T + onehalf*del_th
352     q = del_th
353     c theta2
354     call sw_adtg(S,th,P+onehalf*del_P, adtg_val)
355     del_th = del_P*adtg_val
356    
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 theta3
361     call sw_adtg(S,th,P+onehalf*del_P, adtg_val)
362     del_th = del_P*adtg_val
363     th = th + (1 + 1/sqrt(two))*(del_th - q)
364     q = (two + sqrt(two))*del_th + (-two-three/sqrt(two))*q
365    
366     c theta4
367     call sw_adtg(S,th,P+del_P, adtg_val)
368     del_th = del_P*adtg_val
369     rv = th + (del_th - two*q)/(two*three)
370    
371     RETURN
372     END
373    
374     C======================================================================
375    
376     SUBROUTINE SW_ADTG (S,T,P, rv)
377    
378     c ==================================================================
379     c SUBROUTINE SW_ADTG
380     c ==================================================================
381     c
382     c o Calculates adiabatic temperature gradient as per UNESCO 1983 routines.
383     c
384     c started:
385     c
386     c Armin Koehl akoehl@ucsd.edu
387     c
388     c ==================================================================
389     c SUBROUTINE SW_ADTG
390     c ==================================================================
391    
392     implicit none
393     _RL a0,a1,a2,a3,b0,b1,c0,c1,c2,c3,d0,d1,e0,e1,e2
394     _RL S,T,P
395     _RL sref
396     _RL rv
397    
398     sref = 35. _d 0
399     a0 = 3.5803 _d -5
400     a1 = +8.5258 _d -6
401     a2 = -6.836 _d -8
402     a3 = 6.6228 _d -10
403    
404     b0 = +1.8932 _d -6
405     b1 = -4.2393 _d -8
406    
407     c0 = +1.8741 _d -8
408     c1 = -6.7795 _d -10
409     c2 = +8.733 _d -12
410     c3 = -5.4481 _d -14
411    
412     d0 = -1.1351 _d -10
413     d1 = 2.7759 _d -12
414    
415     e0 = -4.6206 _d -13
416     e1 = +1.8676 _d -14
417     e2 = -2.1687 _d -16
418    
419     rv = a0 + (a1 + (a2 + a3*T)*T)*T
420     & + (b0 + b1*T)*(S-sref)
421     & + ( (c0 + (c1 + (c2 + c3*T)*T)*T) + (d0 + d1*T)*(S-sref) )*P
422     & + ( e0 + (e1 + e2*T)*T )*P*P
423     end

  ViewVC Help
Powered by ViewVC 1.1.22