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

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

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


Revision 1.38 - (show annotations) (download)
Tue Jul 19 12:53:24 2011 UTC (12 years, 9 months ago) by mlosch
Branch: MAIN
CVS Tags: checkpoint63g, checkpoint63h, checkpoint63d, checkpoint63e, checkpoint63f, checkpoint63a, checkpoint63b, checkpoint63c
Changes since 1.37: +150 -1 lines
add code for TEOS-10 (www.teos-10.org, McDougall et al. 2011). Use
this eos with eosType = 'TEOS10', in data (PARM01). This eos implies
that THETA and SALT are "conservative temperature" and "absolute
salinity"

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

  ViewVC Help
Powered by ViewVC 1.1.22