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

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

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


Revision 1.32 - (hide annotations) (download)
Fri Nov 4 01:19:24 2005 UTC (18 years, 7 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint57y_post, checkpoint58, checkpoint57z_post, checkpoint57y_pre, checkpoint57x_post
Changes since 1.31: +2 -2 lines
remove unused variables (reduces number of compiler warning)

1 jmc 1.32 C $Header: /u/gcmpack/MITgcm/model/src/find_rho.F,v 1.31 2005/10/11 00:09:02 jmc Exp $
2 cnh 1.15 C $Name: $
3 cnh 1.1
4 cnh 1.9 #include "CPP_OPTIONS.h"
5 adcroft 1.5 #define USE_FACTORIZED_POLY
6 cnh 1.1
7 cnh 1.15 CBOP
8     C !ROUTINE: FIND_RHO
9     C !INTERFACE:
10 jmc 1.24 SUBROUTINE FIND_RHO(
11 mlosch 1.20 I bi, bj, iMin, iMax, jMin, jMax, k, kRef,
12 adcroft 1.13 I tFld, sFld,
13 adcroft 1.4 O rholoc,
14     I myThid )
15 cnh 1.15
16     C !DESCRIPTION: \bv
17     C *==========================================================*
18     C | o SUBROUTINE FIND_RHO
19 mlosch 1.21 C | Calculates [rho(S,T,z)-rhoConst] of a slice
20 cnh 1.15 C *==========================================================*
21     C |
22     C | k - is the Theta/Salt level
23     C | kRef - determines pressure reference level
24     C | (not used in 'LINEAR' mode)
25     C |
26     C *==========================================================*
27     C \ev
28    
29     C !USES:
30 jmc 1.24 IMPLICIT NONE
31 heimbach 1.12 C == Global variables ==
32 cnh 1.1 #include "SIZE.h"
33 cnh 1.7 #include "EEPARAMS.h"
34 cnh 1.1 #include "PARAMS.h"
35 mlosch 1.16 #include "EOS.h"
36     #include "GRID.h"
37 heimbach 1.12
38 cnh 1.15 C !INPUT/OUTPUT PARAMETERS:
39 heimbach 1.12 C == Routine arguments ==
40 jmc 1.24 C k :: Level of Theta/Salt slice
41     C kRef :: Pressure reference level
42     INTEGER bi,bj,iMin,iMax,jMin,jMax
43     INTEGER k
44     INTEGER kRef
45 adcroft 1.13 _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
46     _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
47 adcroft 1.4 _RL rholoc(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
48 jmc 1.24 INTEGER myThid
49 heimbach 1.12
50 cnh 1.15 C !LOCAL VARIABLES:
51 heimbach 1.12 C == Local variables ==
52 jmc 1.31 INTEGER i,j
53 mlosch 1.21 _RL refTemp,refSalt,sigRef,tP,sP,deltaSig,dRho
54 jmc 1.24 _RL locPres(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
55     _RL rhoP0 (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
56 mlosch 1.16 _RL bulkMod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
57 jmc 1.24 _RL rhoNum (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
58     _RL rhoDen (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
59     CHARACTER*(MAX_LEN_MBUF) msgbuf
60 cnh 1.15 CEOP
61 heimbach 1.11
62     #ifdef ALLOW_AUTODIFF_TAMC
63     DO j=1-OLy,sNy+OLy
64     DO i=1-OLx,sNx+OLx
65 mlosch 1.16 rholoc(i,j) = 0. _d 0
66     rhoP0(i,j) = 0. _d 0
67     bulkMod(i,j) = 0. _d 0
68 heimbach 1.11 ENDDO
69     ENDDO
70     #endif
71 cnh 1.1
72 mlosch 1.22 #ifdef CHECK_SALINITY_FOR_NEGATIVE_VALUES
73     CALL LOOK_FOR_NEG_SALINITY( bi, bj, iMin, iMax, jMin, jMax, k,
74     & sFld, myThid )
75     #endif
76    
77 jmc 1.24 IF (equationOfState.EQ.'LINEAR') THEN
78 adcroft 1.4
79 jmc 1.31 C ***NOTE***
80     C In the linear EOS, to make the static stability calculation meaningful
81     C we alway calculate the perturbation with respect to the surface level.
82     C **********
83     refTemp=tRef(kRef)
84     refSalt=sRef(kRef)
85    
86 mlosch 1.21 dRho = rhoNil-rhoConst
87    
88 jmc 1.24 DO j=jMin,jMax
89     DO i=iMin,iMax
90 mlosch 1.21 rholoc(i,j)=rhoNil*(
91 adcroft 1.13 & sBeta*(sFld(i,j,k,bi,bj)-refSalt)
92     & -tAlpha*(tFld(i,j,k,bi,bj)-refTemp) )
93 mlosch 1.21 & + dRho
94 jmc 1.24 ENDDO
95     ENDDO
96 jmc 1.31
97 jmc 1.24 ELSEIF (equationOfState.EQ.'POLY3') THEN
98 adcroft 1.4
99     refTemp=eosRefT(kRef)
100     refSalt=eosRefS(kRef)
101 mlosch 1.21 sigRef=eosSig0(kRef) + (1000.-rhoConst)
102 adcroft 1.4
103 jmc 1.24 DO j=jMin,jMax
104     DO i=iMin,iMax
105 adcroft 1.13 tP=tFld(i,j,k,bi,bj)-refTemp
106     sP=sFld(i,j,k,bi,bj)-refSalt
107 adcroft 1.5 #ifdef USE_FACTORIZED_POLY
108 adcroft 1.4 deltaSig=
109 adcroft 1.5 & (( eosC(9,kRef)*sP + eosC(5,kRef) )*sP + eosC(2,kRef) )*sP
110     & + ( ( eosC(6,kRef)
111     & *tP
112     & +eosC(7,kRef)*sP + eosC(3,kRef)
113     & )*tP
114     & +(eosC(8,kRef)*sP + eosC(4,kRef) )*sP + eosC(1,kRef)
115     & )*tP
116     #else
117     deltaSig=
118     & eosC(1,kRef)*tP
119     & +eosC(2,kRef) *sP
120     & +eosC(3,kRef)*tP*tP
121     & +eosC(4,kRef)*tP *sP
122     & +eosC(5,kRef) *sP*sP
123     & +eosC(6,kRef)*tP*tP*tP
124     & +eosC(7,kRef)*tP*tP *sP
125     & +eosC(8,kRef)*tP *sP*sP
126     & +eosC(9,kRef) *sP*sP*sP
127     #endif
128     rholoc(i,j)=sigRef+deltaSig
129 jmc 1.24 ENDDO
130     ENDDO
131 adcroft 1.4
132 jmc 1.24 ELSEIF ( equationOfState(1:5).EQ.'JMD95'
133     & .OR. equationOfState.EQ.'UNESCO' ) THEN
134 mlosch 1.16 C nonlinear equation of state in pressure coordinates
135    
136 jmc 1.24 CALL PRESSURE_FOR_EOS(
137     I bi, bj, iMin, iMax, jMin, jMax, kRef,
138     O locPres,
139     I myThid )
140    
141 mlosch 1.16 CALL FIND_RHOP0(
142     I bi, bj, iMin, iMax, jMin, jMax, k,
143     I tFld, sFld,
144     O rhoP0,
145     I myThid )
146    
147     CALL FIND_BULKMOD(
148 jmc 1.24 I bi, bj, iMin, iMax, jMin, jMax, k,
149     I locPres, tFld, sFld,
150 mlosch 1.16 O bulkMod,
151     I myThid )
152    
153 heimbach 1.23 #ifdef ALLOW_AUTODIFF_TAMC
154 edhill 1.26 cph can not DO storing here since find_rho is called multiple times;
155 heimbach 1.23 cph additional recomp. should be acceptable
156     cphCADJ STORE rhoP0(:,:) = comlev1_bibj_k , key=kkey , byte=isbyte
157     cphCADJ STORE bulkMod(:,:) = comlev1_bibj_k , key=kkey , byte=isbyte
158     #endif /* ALLOW_AUTODIFF_TAMC */
159 jmc 1.24 DO j=jMin,jMax
160     DO i=iMin,iMax
161 mlosch 1.16
162     C density of sea water at pressure p
163     rholoc(i,j) = rhoP0(i,j)
164     & /(1. _d 0 -
165 jmc 1.24 & locPres(i,j)*SItoBar/bulkMod(i,j))
166 mlosch 1.21 & - rhoConst
167 mlosch 1.16
168 jmc 1.24 ENDDO
169     ENDDO
170 mlosch 1.16
171 jmc 1.24 ELSEIF ( equationOfState.EQ.'MDJWF' ) THEN
172    
173     CALL PRESSURE_FOR_EOS(
174     I bi, bj, iMin, iMax, jMin, jMax, kRef,
175     O locPres,
176     I myThid )
177 mlosch 1.19
178 jmc 1.24 CALL FIND_RHONUM( bi, bj, iMin, iMax, jMin, jMax, k,
179     & locPres, tFld, sFld, rhoNum, myThid )
180 mlosch 1.19
181 jmc 1.24 CALL FIND_RHODEN( bi, bj, iMin, iMax, jMin, jMax, k,
182     & locPres, tFld, sFld, rhoDen, myThid )
183 heimbach 1.23 #ifdef ALLOW_AUTODIFF_TAMC
184 edhill 1.26 cph can not DO storing here since find_rho is called multiple times;
185 heimbach 1.23 cph additional recomp. should be acceptable
186     cphCADJ STORE rhoNum(:,:) = comlev1_bibj_k , key=kkey , byte=isbyte
187     cphCADJ STORE rhoDen(:,:) = comlev1_bibj_k , key=kkey , byte=isbyte
188     #endif /* ALLOW_AUTODIFF_TAMC */
189 jmc 1.24 DO j=jMin,jMax
190     DO i=iMin,iMax
191 mlosch 1.21 rholoc(i,j) = rhoNum(i,j)*rhoDen(i,j) - rhoConst
192 jmc 1.24 ENDDO
193     ENDDO
194 mlosch 1.19
195 jmc 1.24 ELSEIF( equationOfState .EQ. 'IDEALG' ) THEN
196 mlosch 1.19 C
197 jmc 1.24 ELSE
198     WRITE(msgbuf,'(3a)')
199 mlosch 1.19 & ' FIND_RHO: equationOfState = "',equationOfState,'"'
200 jmc 1.24 CALL print_error( msgbuf, mythid )
201     STOP 'ABNORMAL END: S/R FIND_RHO'
202     ENDIF
203 cnh 1.1
204 jmc 1.24 RETURN
205     END
206 mlosch 1.16
207     CBOP
208     C !ROUTINE: FIND_RHOP0
209     C !INTERFACE:
210 jmc 1.24 SUBROUTINE FIND_RHOP0(
211 mlosch 1.16 I bi, bj, iMin, iMax, jMin, jMax, k,
212     I tFld, sFld,
213     O rhoP0,
214     I myThid )
215    
216     C !DESCRIPTION: \bv
217     C *==========================================================*
218     C | o SUBROUTINE FIND_RHOP0
219     C | Calculates rho(S,T,0) of a slice
220     C *==========================================================*
221     C |
222     C | k - is the surface level
223     C |
224     C *==========================================================*
225     C \ev
226    
227     C !USES:
228 jmc 1.24 IMPLICIT NONE
229 mlosch 1.16 C == Global variables ==
230     #include "SIZE.h"
231     #include "EEPARAMS.h"
232     #include "PARAMS.h"
233     #include "EOS.h"
234    
235     C !INPUT/OUTPUT PARAMETERS:
236     C == Routine arguments ==
237 jmc 1.24 C k :: Level of Theta/Salt slice
238     INTEGER bi,bj,iMin,iMax,jMin,jMax
239     INTEGER k
240 mlosch 1.16 _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
241     _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
242     _RL rhoP0(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
243 jmc 1.24 INTEGER myThid
244 mlosch 1.16
245     C !LOCAL VARIABLES:
246     C == Local variables ==
247 jmc 1.24 INTEGER i,j
248 jmc 1.27 _RL rfresh, rsalt
249 mlosch 1.16 _RL t, t2, t3, t4, s, s3o2
250     CEOP
251    
252 jmc 1.24 DO j=jMin,jMax
253     DO i=iMin,iMax
254 mlosch 1.16 C abbreviations
255     t = tFld(i,j,k,bi,bj)
256     t2 = t*t
257     t3 = t2*t
258     t4 = t3*t
259    
260 adcroft 1.17 s = sFld(i,j,k,bi,bj)
261 jmc 1.25 IF ( s .GT. 0. _d 0 ) THEN
262 jmc 1.24 s3o2 = s*SQRT(s)
263 jmc 1.25 ELSE
264     s = 0. _d 0
265     s3o2 = 0. _d 0
266     ENDIF
267 mlosch 1.16
268     C density of freshwater at the surface
269     rfresh =
270     & eosJMDCFw(1)
271     & + eosJMDCFw(2)*t
272     & + eosJMDCFw(3)*t2
273     & + eosJMDCFw(4)*t3
274     & + eosJMDCFw(5)*t4
275     & + eosJMDCFw(6)*t4*t
276     C density of sea water at the surface
277     rsalt =
278     & s*(
279     & eosJMDCSw(1)
280     & + eosJMDCSw(2)*t
281     & + eosJMDCSw(3)*t2
282     & + eosJMDCSw(4)*t3
283     & + eosJMDCSw(5)*t4
284     & )
285     & + s3o2*(
286     & eosJMDCSw(6)
287     & + eosJMDCSw(7)*t
288     & + eosJMDCSw(8)*t2
289     & )
290     & + eosJMDCSw(9)*s*s
291    
292     rhoP0(i,j) = rfresh + rsalt
293    
294 jmc 1.24 ENDDO
295     ENDDO
296 mlosch 1.16
297 jmc 1.24 RETURN
298     END
299 mlosch 1.16
300     C !ROUTINE: FIND_BULKMOD
301     C !INTERFACE:
302 jmc 1.24 SUBROUTINE FIND_BULKMOD(
303     I bi, bj, iMin, iMax, jMin, jMax, k,
304     I locPres, tFld, sFld,
305 mlosch 1.16 O bulkMod,
306     I myThid )
307    
308     C !DESCRIPTION: \bv
309     C *==========================================================*
310     C | o SUBROUTINE FIND_BULKMOD
311     C | Calculates the secant bulk modulus K(S,T,p) of a slice
312     C *==========================================================*
313     C |
314 jmc 1.24 C | k - is the level of Theta/Salt slice
315 mlosch 1.16 C |
316     C *==========================================================*
317     C \ev
318    
319     C !USES:
320 jmc 1.24 IMPLICIT NONE
321 mlosch 1.16 C == Global variables ==
322     #include "SIZE.h"
323     #include "EEPARAMS.h"
324     #include "PARAMS.h"
325     #include "EOS.h"
326    
327     C !INPUT/OUTPUT PARAMETERS:
328     C == Routine arguments ==
329 jmc 1.24 C k :: Level of Theta/Salt slice
330     INTEGER bi,bj,iMin,iMax,jMin,jMax
331     INTEGER k
332     _RL locPres(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
333 mlosch 1.16 _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
334     _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
335     _RL bulkMod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
336 jmc 1.24 INTEGER myThid
337 mlosch 1.16
338     C !LOCAL VARIABLES:
339     C == Local variables ==
340 jmc 1.24 INTEGER i,j
341 jmc 1.27 _RL bMfresh, bMsalt, bMpres
342 mlosch 1.16 _RL t, t2, t3, t4, s, s3o2, p, p2
343     CEOP
344    
345 jmc 1.24 DO j=jMin,jMax
346     DO i=iMin,iMax
347 mlosch 1.16 C abbreviations
348     t = tFld(i,j,k,bi,bj)
349     t2 = t*t
350     t3 = t2*t
351     t4 = t3*t
352    
353 adcroft 1.17 s = sFld(i,j,k,bi,bj)
354 jmc 1.25 IF ( s .GT. 0. _d 0 ) THEN
355 jmc 1.24 s3o2 = s*SQRT(s)
356 jmc 1.25 ELSE
357     s = 0. _d 0
358     s3o2 = 0. _d 0
359     ENDIF
360 mlosch 1.16 C
361 jmc 1.24 p = locPres(i,j)*SItoBar
362 mlosch 1.16 p2 = p*p
363     C secant bulk modulus of fresh water at the surface
364     bMfresh =
365     & eosJMDCKFw(1)
366     & + eosJMDCKFw(2)*t
367     & + eosJMDCKFw(3)*t2
368     & + eosJMDCKFw(4)*t3
369     & + eosJMDCKFw(5)*t4
370     C secant bulk modulus of sea water at the surface
371     bMsalt =
372     & s*( eosJMDCKSw(1)
373     & + eosJMDCKSw(2)*t
374     & + eosJMDCKSw(3)*t2
375     & + eosJMDCKSw(4)*t3
376     & )
377     & + s3o2*( eosJMDCKSw(5)
378     & + eosJMDCKSw(6)*t
379     & + eosJMDCKSw(7)*t2
380     & )
381     C secant bulk modulus of sea water at pressure p
382     bMpres =
383     & p*( eosJMDCKP(1)
384     & + eosJMDCKP(2)*t
385     & + eosJMDCKP(3)*t2
386     & + eosJMDCKP(4)*t3
387     & )
388     & + p*s*( eosJMDCKP(5)
389     & + eosJMDCKP(6)*t
390     & + eosJMDCKP(7)*t2
391     & )
392     & + p*s3o2*eosJMDCKP(8)
393     & + p2*( eosJMDCKP(9)
394     & + eosJMDCKP(10)*t
395     & + eosJMDCKP(11)*t2
396     & )
397     & + p2*s*( eosJMDCKP(12)
398     & + eosJMDCKP(13)*t
399     & + eosJMDCKP(14)*t2
400     & )
401    
402     bulkMod(i,j) = bMfresh + bMsalt + bMpres
403    
404 jmc 1.24 ENDDO
405     ENDDO
406 mlosch 1.16
407 jmc 1.24 RETURN
408     END
409 mlosch 1.16
410 mlosch 1.19 CBOP
411     C !ROUTINE: FIND_RHONUM
412     C !INTERFACE:
413 jmc 1.24 SUBROUTINE FIND_RHONUM(
414     I bi, bj, iMin, iMax, jMin, jMax, k,
415     I locPres, tFld, sFld,
416 mlosch 1.19 O rhoNum,
417     I myThid )
418    
419     C !DESCRIPTION: \bv
420     C *==========================================================*
421     C | o SUBROUTINE FIND_RHONUM
422     C | Calculates the numerator of the McDougall et al.
423     C | equation of state
424     C | - the code is more or less a copy of MOM4
425     C *==========================================================*
426     C |
427 jmc 1.24 C | k - is the level of Theta/Salt slice
428 mlosch 1.19 C |
429     C *==========================================================*
430     C \ev
431    
432     C !USES:
433 jmc 1.24 IMPLICIT NONE
434 mlosch 1.19 C == Global variables ==
435     #include "SIZE.h"
436     #include "EEPARAMS.h"
437     #include "PARAMS.h"
438     #include "EOS.h"
439    
440     C !INPUT/OUTPUT PARAMETERS:
441     C == Routine arguments ==
442 jmc 1.24 C k :: Level of Theta/Salt slice
443     INTEGER bi,bj,iMin,iMax,jMin,jMax
444     INTEGER k
445     _RL locPres(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
446 mlosch 1.19 _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
447     _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
448     _RL rhoNum(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
449 jmc 1.24 INTEGER myThid
450 mlosch 1.19
451     C !LOCAL VARIABLES:
452     C == Local variables ==
453 jmc 1.24 INTEGER i,j
454 mlosch 1.19 _RL t1, t2, s1, p1
455     CEOP
456 jmc 1.24 DO j=jMin,jMax
457     DO i=iMin,iMax
458 mlosch 1.19 C abbreviations
459     t1 = tFld(i,j,k,bi,bj)
460     t2 = t1*t1
461     s1 = sFld(i,j,k,bi,bj)
462    
463 jmc 1.24 p1 = locPres(i,j)*SItodBar
464 mlosch 1.19
465     rhoNum(i,j) = eosMDJWFnum(0)
466     & + t1*(eosMDJWFnum(1)
467     & + t1*(eosMDJWFnum(2) + eosMDJWFnum(3)*t1) )
468     & + s1*(eosMDJWFnum(4)
469     & + eosMDJWFnum(5)*t1 + eosMDJWFnum(6)*s1)
470     & + p1*(eosMDJWFnum(7) + eosMDJWFnum(8)*t2
471     & + eosMDJWFnum(9)*s1
472     & + p1*(eosMDJWFnum(10) + eosMDJWFnum(11)*t2) )
473    
474 jmc 1.24 ENDDO
475     ENDDO
476 mlosch 1.19
477 jmc 1.24 RETURN
478 mlosch 1.19 end
479    
480     CBOP
481     C !ROUTINE: FIND_RHODEN
482     C !INTERFACE:
483 jmc 1.24 SUBROUTINE FIND_RHODEN(
484     I bi, bj, iMin, iMax, jMin, jMax, k,
485     I locPres, tFld, sFld,
486 mlosch 1.19 O rhoDen,
487     I myThid )
488    
489     C !DESCRIPTION: \bv
490     C *==========================================================*
491     C | o SUBROUTINE FIND_RHODEN
492     C | Calculates the denominator of the McDougall et al.
493     C | equation of state
494     C | - the code is more or less a copy of MOM4
495     C *==========================================================*
496     C |
497 jmc 1.24 C | k - is the level of Theta/Salt slice
498 mlosch 1.19 C |
499     C *==========================================================*
500     C \ev
501    
502     C !USES:
503 jmc 1.24 IMPLICIT NONE
504 mlosch 1.19 C == Global variables ==
505     #include "SIZE.h"
506     #include "EEPARAMS.h"
507     #include "PARAMS.h"
508     #include "EOS.h"
509    
510     C !INPUT/OUTPUT PARAMETERS:
511     C == Routine arguments ==
512 jmc 1.24 C k :: Level of Theta/Salt slice
513     INTEGER bi,bj,iMin,iMax,jMin,jMax
514     INTEGER k
515     _RL locPres(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
516 mlosch 1.19 _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
517     _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
518     _RL rhoDen(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
519 jmc 1.24 INTEGER myThid
520 mlosch 1.19
521     C !LOCAL VARIABLES:
522     C == Local variables ==
523 jmc 1.24 INTEGER i,j
524 mlosch 1.19 _RL t1, t2, s1, sp5, p1, p1t1
525     _RL den, epsln
526     parameter ( epsln = 0. _d 0 )
527     CEOP
528 jmc 1.24 DO j=jMin,jMax
529     DO i=iMin,iMax
530 mlosch 1.19 C abbreviations
531     t1 = tFld(i,j,k,bi,bj)
532     t2 = t1*t1
533     s1 = sFld(i,j,k,bi,bj)
534 jmc 1.25 IF ( s1 .GT. 0. _d 0 ) THEN
535 jmc 1.24 sp5 = SQRT(s1)
536 jmc 1.25 ELSE
537     s1 = 0. _d 0
538     sp5 = 0. _d 0
539     ENDIF
540 mlosch 1.19
541 jmc 1.24 p1 = locPres(i,j)*SItodBar
542 mlosch 1.19 p1t1 = p1*t1
543    
544     den = eosMDJWFden(0)
545     & + t1*(eosMDJWFden(1)
546     & + t1*(eosMDJWFden(2)
547     & + t1*(eosMDJWFden(3) + t1*eosMDJWFden(4) ) ) )
548     & + s1*(eosMDJWFden(5)
549     & + t1*(eosMDJWFden(6)
550     & + eosMDJWFden(7)*t2)
551     & + sp5*(eosMDJWFden(8) + eosMDJWFden(9)*t2) )
552     & + p1*(eosMDJWFden(10)
553     & + p1t1*(eosMDJWFden(11)*t2 + eosMDJWFden(12)*p1) )
554    
555     rhoDen(i,j) = 1.0/(epsln+den)
556    
557 jmc 1.24 ENDDO
558     ENDDO
559 mlosch 1.19
560 jmc 1.24 RETURN
561 mlosch 1.19 end
562 mlosch 1.20
563 jmc 1.24 SUBROUTINE find_rho_scalar(
564 mlosch 1.20 I tLoc, sLoc, pLoc,
565     O rhoLoc,
566     I myThid )
567    
568     C !DESCRIPTION: \bv
569     C *==========================================================*
570     C | o SUBROUTINE FIND_RHO_SCALAR
571 mlosch 1.21 C | Calculates [rho(S,T,p)-rhoConst]
572 mlosch 1.20 C *==========================================================*
573     C \ev
574    
575     C !USES:
576 jmc 1.24 IMPLICIT NONE
577 mlosch 1.20 C == Global variables ==
578     #include "SIZE.h"
579     #include "EEPARAMS.h"
580     #include "PARAMS.h"
581     #include "EOS.h"
582    
583     C !INPUT/OUTPUT PARAMETERS:
584     C == Routine arguments ==
585     _RL sLoc, tLoc, pLoc
586     _RL rhoLoc
587 jmc 1.24 INTEGER myThid
588 mlosch 1.20
589     C !LOCAL VARIABLES:
590     C == Local variables ==
591    
592     _RL t1, t2, t3, t4, s1, s3o2, p1, p2, sp5, p1t1
593     _RL rfresh, rsalt, rhoP0
594     _RL bMfresh, bMsalt, bMpres, BulkMod
595     _RL rhoNum, rhoDen, den, epsln
596     parameter ( epsln = 0. _d 0 )
597    
598     character*(max_len_mbuf) msgbuf
599     CEOP
600    
601     rhoLoc = 0. _d 0
602     rhoP0 = 0. _d 0
603     bulkMod = 0. _d 0
604     rfresh = 0. _d 0
605     rsalt = 0. _d 0
606     bMfresh = 0. _d 0
607     bMsalt = 0. _d 0
608     bMpres = 0. _d 0
609     rhoNum = 0. _d 0
610     rhoDen = 0. _d 0
611     den = 0. _d 0
612    
613     t1 = tLoc
614     t2 = t1*t1
615     t3 = t2*t1
616     t4 = t3*t1
617    
618     s1 = sLoc
619 jmc 1.24 IF ( s1 .LT. 0. _d 0 ) THEN
620 mlosch 1.20 C issue a warning
621 jmc 1.24 WRITE(*,'(a,i3,a,i3,a,i3,a,e13.5)')
622 mlosch 1.20 & ' FIND_RHO_SCALAR: WARNING, salinity = ', s1
623     s1 = 0. _d 0
624 jmc 1.24 ENDIF
625 mlosch 1.20
626 jmc 1.24 IF (equationOfState.EQ.'LINEAR') THEN
627 mlosch 1.20
628     rhoLoc = 0. _d 0
629    
630 jmc 1.24 ELSEIF (equationOfState.EQ.'POLY3') THEN
631 mlosch 1.20
632     C this is not correct, there is a field eosSig0 which should be use here
633 jmc 1.24 C but I DO not intent to include the reference level in this routine
634     WRITE(*,'(a)')
635 mlosch 1.20 & ' FIND_RHO_SCALAR: for POLY3, the density is not'
636 jmc 1.24 WRITE(*,'(a)')
637 mlosch 1.20 & ' computed correctly in this routine'
638     rhoLoc = 0. _d 0
639    
640 jmc 1.24 ELSEIF ( equationOfState(1:5).EQ.'JMD95'
641     & .OR. equationOfState.EQ.'UNESCO' ) THEN
642 mlosch 1.20 C nonlinear equation of state in pressure coordinates
643    
644 jmc 1.24 s3o2 = s1*SQRT(s1)
645 mlosch 1.20
646     p1 = pLoc*SItoBar
647     p2 = p1*p1
648    
649     C density of freshwater at the surface
650     rfresh =
651     & eosJMDCFw(1)
652     & + eosJMDCFw(2)*t1
653     & + eosJMDCFw(3)*t2
654     & + eosJMDCFw(4)*t3
655     & + eosJMDCFw(5)*t4
656     & + eosJMDCFw(6)*t4*t1
657     C density of sea water at the surface
658     rsalt =
659     & s1*(
660     & eosJMDCSw(1)
661     & + eosJMDCSw(2)*t1
662     & + eosJMDCSw(3)*t2
663     & + eosJMDCSw(4)*t3
664     & + eosJMDCSw(5)*t4
665     & )
666     & + s3o2*(
667     & eosJMDCSw(6)
668     & + eosJMDCSw(7)*t1
669     & + eosJMDCSw(8)*t2
670     & )
671     & + eosJMDCSw(9)*s1*s1
672    
673     rhoP0 = rfresh + rsalt
674    
675     C secant bulk modulus of fresh water at the surface
676     bMfresh =
677     & eosJMDCKFw(1)
678     & + eosJMDCKFw(2)*t1
679     & + eosJMDCKFw(3)*t2
680     & + eosJMDCKFw(4)*t3
681     & + eosJMDCKFw(5)*t4
682     C secant bulk modulus of sea water at the surface
683     bMsalt =
684     & s1*( eosJMDCKSw(1)
685     & + eosJMDCKSw(2)*t1
686     & + eosJMDCKSw(3)*t2
687     & + eosJMDCKSw(4)*t3
688     & )
689     & + s3o2*( eosJMDCKSw(5)
690     & + eosJMDCKSw(6)*t1
691     & + eosJMDCKSw(7)*t2
692     & )
693     C secant bulk modulus of sea water at pressure p
694     bMpres =
695     & p1*( eosJMDCKP(1)
696     & + eosJMDCKP(2)*t1
697     & + eosJMDCKP(3)*t2
698     & + eosJMDCKP(4)*t3
699     & )
700     & + p1*s1*( eosJMDCKP(5)
701     & + eosJMDCKP(6)*t1
702     & + eosJMDCKP(7)*t2
703     & )
704     & + p1*s3o2*eosJMDCKP(8)
705     & + p2*( eosJMDCKP(9)
706     & + eosJMDCKP(10)*t1
707     & + eosJMDCKP(11)*t2
708     & )
709     & + p2*s1*( eosJMDCKP(12)
710     & + eosJMDCKP(13)*t1
711     & + eosJMDCKP(14)*t2
712     & )
713    
714     bulkMod = bMfresh + bMsalt + bMpres
715    
716     C density of sea water at pressure p
717 mlosch 1.21 rhoLoc = rhoP0/(1. _d 0 - p1/bulkMod) - rhoConst
718 mlosch 1.20
719 jmc 1.24 ELSEIF ( equationOfState.EQ.'MDJWF' ) THEN
720 mlosch 1.20
721 jmc 1.24 sp5 = SQRT(s1)
722 mlosch 1.20
723     p1 = pLoc*SItodBar
724     p1t1 = p1*t1
725    
726     rhoNum = eosMDJWFnum(0)
727     & + t1*(eosMDJWFnum(1)
728     & + t1*(eosMDJWFnum(2) + eosMDJWFnum(3)*t1) )
729     & + s1*(eosMDJWFnum(4)
730     & + eosMDJWFnum(5)*t1 + eosMDJWFnum(6)*s1)
731     & + p1*(eosMDJWFnum(7) + eosMDJWFnum(8)*t2
732     & + eosMDJWFnum(9)*s1
733     & + p1*(eosMDJWFnum(10) + eosMDJWFnum(11)*t2) )
734    
735    
736     den = eosMDJWFden(0)
737     & + t1*(eosMDJWFden(1)
738     & + t1*(eosMDJWFden(2)
739     & + t1*(eosMDJWFden(3) + t1*eosMDJWFden(4) ) ) )
740     & + s1*(eosMDJWFden(5)
741     & + t1*(eosMDJWFden(6)
742     & + eosMDJWFden(7)*t2)
743     & + sp5*(eosMDJWFden(8) + eosMDJWFden(9)*t2) )
744     & + p1*(eosMDJWFden(10)
745     & + p1t1*(eosMDJWFden(11)*t2 + eosMDJWFden(12)*p1) )
746    
747     rhoDen = 1.0/(epsln+den)
748    
749 mlosch 1.21 rhoLoc = rhoNum*rhoDen - rhoConst
750 mlosch 1.20
751 jmc 1.24 ELSEIF( equationOfState .EQ. 'IDEALG' ) THEN
752 mlosch 1.20 C
753 jmc 1.24 ELSE
754     WRITE(msgbuf,'(3A)')
755 mlosch 1.20 & ' FIND_RHO_SCALAR : equationOfState = "',
756     & equationOfState,'"'
757 jmc 1.24 CALL PRINT_ERROR( msgbuf, mythid )
758     STOP 'ABNORMAL END: S/R FIND_RHO_SCALAR'
759     ENDIF
760 mlosch 1.22
761 jmc 1.24 RETURN
762     END
763 mlosch 1.22
764     CBOP
765     C !ROUTINE: LOOK_FOR_NEG_SALINITY
766     C !INTERFACE:
767 jmc 1.24 SUBROUTINE LOOK_FOR_NEG_SALINITY(
768 mlosch 1.22 I bi, bj, iMin, iMax, jMin, jMax, k,
769     I sFld,
770     I myThid )
771    
772     C !DESCRIPTION: \bv
773     C *==========================================================*
774     C | o SUBROUTINE LOOK_FOR_NEG_SALINITY
775     C | looks for and fixes negative salinity values
776 jmc 1.24 C | this is necessary IF the equation of state uses
777 mlosch 1.22 C | the square root of salinity
778     C *==========================================================*
779     C |
780     C | k - is the Salt level
781     C |
782     C *==========================================================*
783     C \ev
784    
785     C !USES:
786 jmc 1.24 IMPLICIT NONE
787 mlosch 1.22 C == Global variables ==
788     #include "SIZE.h"
789     #include "EEPARAMS.h"
790     #include "PARAMS.h"
791     #include "GRID.h"
792    
793     C !INPUT/OUTPUT PARAMETERS:
794     C == Routine arguments ==
795 jmc 1.24 C k :: Level of Theta/Salt slice
796     INTEGER bi,bj,iMin,iMax,jMin,jMax
797     INTEGER k
798 mlosch 1.22 _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
799 jmc 1.24 INTEGER myThid
800 mlosch 1.22
801     C !LOCAL VARIABLES:
802     C == Local variables ==
803 jmc 1.24 INTEGER i,j, localWarning
804 jmc 1.32 c character*(max_len_mbuf) msgbuf
805 mlosch 1.22 CEOP
806    
807     localWarning = 0
808 jmc 1.24 DO j=jMin,jMax
809     DO i=iMin,iMax
810 mlosch 1.22 C abbreviations
811 jmc 1.24 IF ( sFld(i,j,k,bi,bj) .LT. 0. _d 0 ) THEN
812 mlosch 1.22 localWarning = localWarning + 1
813     sFld(i,j,k,bi,bj) = 0. _d 0
814 jmc 1.24 ENDIF
815     ENDDO
816     ENDDO
817 mlosch 1.22 C issue a warning
818 jmc 1.24 IF ( localWarning .GT. 0 ) THEN
819     WRITE(standardMessageUnit,'(A,A)')
820 mlosch 1.22 & 'S/R LOOK_FOR_NEG_SALINITY: found negative salinity',
821     & 'values and reset them to zero.'
822 jmc 1.24 WRITE(standardMessageUnit,'(A,I3)')
823 mlosch 1.22 & 'S/R LOOK_FOR_NEG_SALINITY: current level is k = ', k
824 jmc 1.24 ENDIF
825 mlosch 1.20
826 jmc 1.24 RETURN
827     END

  ViewVC Help
Powered by ViewVC 1.1.22