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

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

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


Revision 1.1 - (hide annotations) (download)
Tue Mar 18 19:45:01 2008 UTC (17 years, 4 months ago) by utke
Branch: MAIN
CVS Tags: HEAD
treat all exch routines in split mode,
retain the templates directive that already are in the aMPI stubs

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

  ViewVC Help
Powered by ViewVC 1.1.22