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