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

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

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


Revision 1.20 - (show annotations) (download)
Sat Sep 6 17:42:27 2008 UTC (15 years, 7 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint62v, checkpoint62u, checkpoint62t, checkpoint62c, checkpoint62s, checkpoint62r, checkpoint62q, checkpoint62p, checkpoint62a, checkpoint62g, checkpoint62f, checkpoint62e, checkpoint62d, checkpoint62k, checkpoint62j, checkpoint62i, checkpoint62h, checkpoint62o, checkpoint62n, checkpoint62m, checkpoint62l, checkpoint62w, checkpoint62z, checkpoint62y, checkpoint62x, checkpoint62, checkpoint63, checkpoint62b, checkpoint61f, checkpoint61n, checkpoint61q, checkpoint61e, checkpoint61g, checkpoint61d, checkpoint61o, checkpoint61l, checkpoint61m, checkpoint61j, checkpoint61k, checkpoint61h, checkpoint61i, checkpoint61v, checkpoint61w, checkpoint61t, checkpoint61u, checkpoint61r, checkpoint61s, checkpoint61p, checkpoint61z, checkpoint61x, checkpoint61y
Changes since 1.19: +2 -3 lines
change FIND_RHO_SCALAR : return rho (instead of rho - rhoConst)

1 C $Header: /u/gcmpack/MITgcm/model/src/ini_eos.F,v 1.19 2008/08/21 14:25:34 mlosch Exp $
2 C $Name: $
3
4 #include "CPP_OPTIONS.h"
5 #undef INCLUDE_EOS_CHECK
6
7 C-- File ini_eos.F: Routines to initialise Equation of State
8 C-- Contents
9 C-- o INI_EOS
10 C-- o EOS_CHECK
11
12 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
13 CBOP
14 C !ROUTINE: INI_EOS
15 C !INTERFACE:
16 SUBROUTINE INI_EOS( myThid )
17 C !DESCRIPTION: \bv
18 C *==========================================================*
19 C | SUBROUTINE INI_EOS
20 C | o Initialise coefficients of equation of state.
21 C *==========================================================*
22 C \ev
23
24 C !USES:
25
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 C myThid - Number of this instance of INI_CORI
36 INTEGER myThid
37
38 C !LOCAL VARIABLES:
39 C == Local variables ==
40 C i,k :: Loop counters
41 INTEGER i, k
42 CHARACTER*(MAX_LEN_MBUF) msgBuf
43
44 IF ( .NOT.fluidIsWater ) RETURN
45
46 _BARRIER
47 _BEGIN_MASTER(myThid)
48
49 equationOfState = eosType
50
51 DO k = 1,6
52 eosJMDCFw(k) = 0. _d 0
53 ENDDO
54 DO k = 1,9
55 eosJMDCSw(k) = 0. _d 0
56 ENDDO
57 DO k = 1,5
58 eosJMDCKFw(k) = 0. _d 0
59 ENDDO
60 DO k = 1,7
61 eosJMDCKSw(k) = 0. _d 0
62 ENDDO
63 DO k = 1,14
64 eosJMDCKP(k) = 0. _d 0
65 ENDDO
66 DO k = 0,11
67 eosMDJWFnum(k) = 0. _d 0
68 ENDDO
69 DO k = 0,12
70 eosMDJWFden(k) = 0. _d 0
71 ENDDO
72
73 IF ( equationOfState .EQ. 'LINEAR' ) THEN
74 IF ( tAlpha .EQ. UNSET_RL ) tAlpha = 2. _d -4
75 IF ( sBeta .EQ. UNSET_RL ) sBeta = 7.4 _d -4
76 ELSEIF ( equationOfState .EQ. 'POLY3' ) THEN
77 OPEN(37,FILE='POLY3.COEFFS',STATUS='OLD',FORM='FORMATTED')
78 READ(37,*) k
79 IF (k.NE.Nr) THEN
80 WRITE(msgBuf,'(A)')
81 & 'ini_eos: attempt to read POLY3.COEFFS failed'
82 CALL PRINT_ERROR( msgBuf , myThid )
83 WRITE(msgBuf,'(A)')
84 & ' because bad # of levels in data'
85 CALL PRINT_ERROR( msgBuf , myThid )
86 STOP 'Bad data in POLY3.COEFFS'
87 ENDIF
88 READ(37,*) (eosRefT(k),eosRefS(k),eosSig0(k),k=1,Nr)
89 DO k=1,Nr
90 READ(37,*) (eosC(i,k),i=1,9)
91 ENDDO
92 CLOSE(37)
93
94 ELSEIF ( equationOfState .EQ. 'JMD95Z'
95 & .OR. equationOfState .EQ. 'JMD95P'
96 & .OR. equationOfState .EQ. 'UNESCO' ) THEN
97 C
98 C Jackett & McDougall (1995, JPO) equation of state
99 C rho = R(salinity, potential temperature, pressure)
100 C pressure needs to be available (from the previous
101 C time step to linearize the problem)
102 C
103 IF ( equationOfState .EQ. 'JMD95Z' .AND. usingPCoords ) THEN
104 WRITE(msgBuf,'(A)')
105 & 'ini_eos: equation of state ''JMD95Z'' should not'
106 CALL PRINT_ERROR( msgBuf , myThid )
107 WRITE(msgBuf,'(A)')
108 & ' be used together with pressure coordinates.'
109 CALL PRINT_ERROR( msgBuf , myThid )
110 WRITE(msgBuf,'(A)')
111 & ' Use only ''JMD95P'' with ''OCEANICP''.'
112 CALL PRINT_ERROR( msgBuf , myThid )
113 STOP 'ABNORMAL END: S/R INI_EOS'
114 ENDIF
115
116 C coefficients nonlinear equation of state in pressure coordinates for
117 C 1. density of fresh water at p = 0
118 eosJMDCFw(1) = 999.842594 _d +00
119 eosJMDCFw(2) = 6.793952 _d -02
120 eosJMDCFw(3) = - 9.095290 _d -03
121 eosJMDCFw(4) = 1.001685 _d -04
122 eosJMDCFw(5) = - 1.120083 _d -06
123 eosJMDCFw(6) = 6.536332 _d -09
124 C 2. density of sea water at p = 0
125 eosJMDCSw(1) = 8.24493 _d -01
126 eosJMDCSw(2) = - 4.0899 _d -03
127 eosJMDCSw(3) = 7.6438 _d -05
128 eosJMDCSw(4) = - 8.2467 _d -07
129 eosJMDCSw(5) = 5.3875 _d -09
130 eosJMDCSw(6) = - 5.72466 _d -03
131 eosJMDCSw(7) = 1.0227 _d -04
132 eosJMDCSw(8) = - 1.6546 _d -06
133 eosJMDCSw(9) = 4.8314 _d -04
134 IF ( equationOfState(1:5) .EQ. 'JMD95' ) THEN
135 C 3. secant bulk modulus K of fresh water at p = 0
136 eosJMDCKFw(1) = 1.965933 _d +04
137 eosJMDCKFw(2) = 1.444304 _d +02
138 eosJMDCKFw(3) = - 1.706103 _d +00
139 eosJMDCKFw(4) = 9.648704 _d -03
140 eosJMDCKFw(5) = - 4.190253 _d -05
141 C 4. secant bulk modulus K of sea water at p = 0
142 eosJMDCKSw(1) = 5.284855 _d +01
143 eosJMDCKSw(2) = - 3.101089 _d -01
144 eosJMDCKSw(3) = 6.283263 _d -03
145 eosJMDCKSw(4) = - 5.084188 _d -05
146 eosJMDCKSw(5) = 3.886640 _d -01
147 eosJMDCKSw(6) = 9.085835 _d -03
148 eosJMDCKSw(7) = - 4.619924 _d -04
149 C 5. secant bulk modulus K of sea water at p
150 eosJMDCKP( 1) = 3.186519 _d +00
151 eosJMDCKP( 2) = 2.212276 _d -02
152 eosJMDCKP( 3) = - 2.984642 _d -04
153 eosJMDCKP( 4) = 1.956415 _d -06
154 eosJMDCKP( 5) = 6.704388 _d -03
155 eosJMDCKP( 6) = - 1.847318 _d -04
156 eosJMDCKP( 7) = 2.059331 _d -07
157 eosJMDCKP( 8) = 1.480266 _d -04
158 eosJMDCKP( 9) = 2.102898 _d -04
159 eosJMDCKP(10) = - 1.202016 _d -05
160 eosJMDCKP(11) = 1.394680 _d -07
161 eosJMDCKP(12) = - 2.040237 _d -06
162 eosJMDCKP(13) = 6.128773 _d -08
163 eosJMDCKP(14) = 6.207323 _d -10
164
165 ELSEIF ( equationOfState .EQ. 'UNESCO' ) THEN
166
167 WRITE(msgBuf,'(a)')
168 & 'WARNING WARNING WARNING WARNING WARNING WARNING '
169 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
170 & SQUEEZE_RIGHT , myThid )
171 WRITE(msgBuf,'(a,a)')
172 & 'WARNING: using the UNESCO formula with potential ',
173 & 'temperature'
174 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
175 & SQUEEZE_RIGHT , myThid )
176 WRITE(msgBuf,'(a)')
177 & 'WARNING: can result in density errors of up to 5%'
178 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
179 & SQUEEZE_RIGHT , myThid )
180 WRITE(msgBuf,'(a)')
181 & 'WARNING: (see Jackett and McDougall 1995, JAOT)'
182 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
183 & SQUEEZE_RIGHT , myThid )
184 WRITE(msgBuf,'(a)')
185 & 'WARNING WARNING WARNING WARNING WARNING WARNING '
186 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
187 & SQUEEZE_RIGHT , myThid )
188
189 C 3. secant bulk modulus K of fresh water at p = 0
190 eosJMDCKFw(1) = 1.965221 _d +04
191 eosJMDCKFw(2) = 1.484206 _d +02
192 eosJMDCKFw(3) = - 2.327105 _d +00
193 eosJMDCKFw(4) = 1.360477 _d -02
194 eosJMDCKFw(5) = - 5.155288 _d -05
195 C 4. secant bulk modulus K of sea water at p = 0
196 eosJMDCKSw(1) = 5.46746 _d +01
197 eosJMDCKSw(2) = - 0.603459 _d +00
198 eosJMDCKSw(3) = 1.09987 _d -02
199 eosJMDCKSw(4) = - 6.1670 _d -05
200 eosJMDCKSw(5) = 7.944 _d -02
201 eosJMDCKSw(6) = 1.6483 _d -02
202 eosJMDCKSw(7) = - 5.3009 _d -04
203 C 5. secant bulk modulus K of sea water at p
204 eosJMDCKP( 1) = 3.239908 _d +00
205 eosJMDCKP( 2) = 1.43713 _d -03
206 eosJMDCKP( 3) = 1.16092 _d -04
207 eosJMDCKP( 4) = - 5.77905 _d -07
208 eosJMDCKP( 5) = 2.2838 _d -03
209 eosJMDCKP( 6) = - 1.0981 _d -05
210 eosJMDCKP( 7) = - 1.6078 _d -06
211 eosJMDCKP( 8) = 1.91075 _d -04
212 eosJMDCKP( 9) = 8.50935 _d -05
213 eosJMDCKP(10) = - 6.12293 _d -06
214 eosJMDCKP(11) = 5.2787 _d -08
215 eosJMDCKP(12) = - 9.9348 _d -07
216 eosJMDCKP(13) = 2.0816 _d -08
217 eosJMDCKP(14) = 9.1697 _d -10
218 ELSE
219 STOP 'INI_EOS: We should never reach this point!'
220 ENDIF
221
222 ELSEIF ( equationOfState .EQ. 'MDJWF' ) THEN
223
224 eosMDJWFnum( 0) = 9.99843699 _d +02
225 eosMDJWFnum( 1) = 7.35212840 _d +00
226 eosMDJWFnum( 2) = -5.45928211 _d -02
227 eosMDJWFnum( 3) = 3.98476704 _d -04
228 eosMDJWFnum( 4) = 2.96938239 _d +00
229 eosMDJWFnum( 5) = -7.23268813 _d -03
230 eosMDJWFnum( 6) = 2.12382341 _d -03
231 eosMDJWFnum( 7) = 1.04004591 _d -02
232 eosMDJWFnum( 8) = 1.03970529 _d -07
233 eosMDJWFnum( 9) = 5.18761880 _d -06
234 eosMDJWFnum(10) = -3.24041825 _d -08
235 eosMDJWFnum(11) = -1.23869360 _d -11
236
237
238 eosMDJWFden( 0) = 1.00000000 _d +00
239 eosMDJWFden( 1) = 7.28606739 _d -03
240 eosMDJWFden( 2) = -4.60835542 _d -05
241 eosMDJWFden( 3) = 3.68390573 _d -07
242 eosMDJWFden( 4) = 1.80809186 _d -10
243 eosMDJWFden( 5) = 2.14691708 _d -03
244 eosMDJWFden( 6) = -9.27062484 _d -06
245 eosMDJWFden( 7) = -1.78343643 _d -10
246 eosMDJWFden( 8) = 4.76534122 _d -06
247 eosMDJWFden( 9) = 1.63410736 _d -09
248 eosMDJWFden(10) = 5.30848875 _d -06
249 eosMDJWFden(11) = -3.03175128 _d -16
250 eosMDJWFden(12) = -1.27934137 _d -17
251
252 ELSEIF( equationOfState .EQ. 'IDEALG' ) THEN
253
254 ELSE
255
256 WRITE(msgbuf,'(3A)') ' INI_EOS: equationOfState = "',
257 & equationOfState,'"'
258 CALL PRINT_ERROR( msgbuf, myThid )
259 STOP 'ABNORMAL END: S/R INI_EOS'
260
261 ENDIF
262
263 C-- Check EOS initialisation:
264
265 CALL EOS_CHECK( myThid )
266
267 _END_MASTER( myThid )
268 _BARRIER
269
270 RETURN
271 END
272
273 CBOP
274 C !ROUTINE: EOS_CHECK
275 C !INTERFACE:
276 SUBROUTINE EOS_CHECK( myThid )
277 C !DESCRIPTION: \bv
278 C *==========================================================*
279 C | SUBROUTINE EOS_CHECK
280 C | o check the equation of state.
281 C *==========================================================*
282 C \ev
283
284 C !USES:
285
286 IMPLICIT NONE
287 #include "SIZE.h"
288 #include "EEPARAMS.h"
289 #include "PARAMS.h"
290 #include "EOS.h"
291 #include "GRID.h"
292 #include "DYNVARS.h"
293
294 C !INPUT/OUTPUT PARAMETERS:
295 C == Routine arguments ==
296 C myThid :: Number of this instance of EOS_CHECK
297 INTEGER myThid
298
299 #ifdef INCLUDE_EOS_CHECK
300 C !LOCAL VARIABLES:
301 C == Local variables ==
302 C bi,bj - Loop counters
303 C i,j,k
304 INTEGER bi, bj
305 INTEGER iMin, iMax, jMin, jMax
306 INTEGER i, j, k
307 _RL tFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
308 _RL sFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
309 _RL pFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
310 _RL rhoLoc (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
311 _RL bulkMod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
312
313 INTEGER ncheck, kcheck
314 PARAMETER ( ncheck = 13 )
315 _RL tLoc(ncheck), ptLoc(ncheck), sLoc(ncheck), pLoc(ncheck)
316 _RL rLoc(ncheck), bLoc(ncheck)
317 _RS mskSave
318 _RS rC_Save
319
320 CHARACTER*(MAX_LEN_MBUF) msgBuf
321
322 DATA tLoc
323 & /3.25905152915860 _d 0, 20.38687090048638 _d 0,
324 & 25.44820830309568 _d 0, 20.17368557065936 _d 0,
325 & 13.43397459640398 _d 0,
326 & 5. _d 0, 25. _d 0,
327 & 5. _d 0, 25. _d 0,
328 & 5. _d 0, 25. _d 0,
329 & 5. _d 0, 25. _d 0/,
330 & ptLoc
331 & /3. _d 0, 20. _d 0,
332 & 25. _d 0, 20. _d 0,
333 & 12. _d 0,
334 & 5. _d 0, 25. _d 0,
335 & 5. _d 0, 25. _d 0,
336 & 4.03692566635316 _d 0, 22.84661726775120 _d 0,
337 & 3.62720389416752 _d 0, 22.62420229124846 _d 0/
338 & sLoc
339 & /35.5 _d 0, 35. _d 0,
340 & 35.0 _d 0, 20. _d 0,
341 & 40.0 _d 0,
342 & 0. _d 0, 0. _d 0,
343 & 35. _d 0, 35. _d 0,
344 & 0. _d 0, 0. _d 0,
345 & 35. _d 0, 35. _d 0/
346 & pLoc
347 & /300. _d 5, 200. _d 5,
348 & 200. _d 5, 100. _d 5,
349 & 800. _d 5,
350 & 0. _d 0, 0. _d 0,
351 & 0. _d 0, 0. _d 0,
352 & 1000. _d 5, 1000. _d 5,
353 & 1000. _d 5, 1000. _d 5/
354 DATA rLoc
355 & /1041.83267 _d 0, 1033.213387 _d 0,
356 & 1031.654229 _d 0, 1017.726743 _d 0,
357 & 1062.928258 _d 0,
358 & 999.96675 _d 0, 997.04796 _d 0,
359 & 1027.67547 _d 0, 1023.34306 _d 0,
360 & 1044.12802 _d 0, 1037.90204 _d 0,
361 & 1069.48914 _d 0, 1062.53817 _d 0/
362 & bLoc
363 & / -1.00000 _d 0, -1.00000 _d 0,
364 & -1.00000 _d 0, -1.00000 _d 0,
365 & -1.00000 _d 0,
366 & 20337.80375 _d 0, 22100.72106 _d 0,
367 & 22185.93358 _d 0, 23726.34949 _d 0,
368 & 23643.52599 _d 0, 25405.09717 _d 0,
369 & 25577.49819 _d 0, 27108.94504 _d 0/
370
371
372 bi = 1
373 bj = 1
374 k = 1
375 iMin = 1
376 iMax = 1
377 jMin = 1
378 jMax = 1
379 i = 1
380 j = 1
381 IF ( equationOfState.NE.'LINEAR'
382 & .AND. equationOfState.NE.'POLY3'
383 & .AND. equationOfState.NE.'IDEALG' ) THEN
384 C check nonlinear EOS
385 WRITE(msgBuf,'(A,A)')
386 & 'EOS_CHECK: Check the equation of state: Type ',
387 & equationOfState
388 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
389 & SQUEEZE_RIGHT , myThid )
390 C- to get the right pressure out of S/R PRESSURE_FOR_EOS:
391 mskSave = maskC(i,j,k,bi,bj)
392 rC_Save = rC(k)
393 maskC(i,j,k,bi,bj) = 1.
394 totPhiHyd(i,j,k,bi,bj) = 0.
395 C- + set rC accordingly
396 DO kcheck = 1,ncheck
397 IF ( usingZCoords ) THEN
398 rC(k) = -pLoc(kcheck)*recip_rhoConst*recip_gravity
399 ELSE
400 rC(k) = pLoc(kcheck)
401 ENDIF
402 IF ( equationOfState.NE.'UNESCO' ) THEN
403 tFld(i,j) = ptLoc(kcheck)
404 ELSE
405 tFld(i,j) = tLoc(kcheck)
406 ENDIF
407 sFld(i,j) = sLoc(kcheck)
408 pFld(i,j) = pLoc(kcheck)
409 rhoLoc(i,j) = 0. _d 0
410
411 CALL FIND_RHO_2D(
412 I iMin, iMax, jMin, jMax, k,
413 I tFld, sFld,
414 O rhoLoc,
415 I k, bi, bj, myThid )
416
417 IF ( equationOfState.EQ.'MDJWF' ) THEN
418 bulkMod(i,j) = -1. _d 0
419 ELSE
420 CALL FIND_BULKMOD(
421 I iMin, iMax, jMin, jMax,
422 I pFld, tFld, sFld,
423 O bulkMod,
424 I myThid )
425 ENDIF
426
427 IF ( kcheck .EQ. 1 ) THEN
428 WRITE(msgBuf,'(A)')
429 & 'EOS_CHECK: check values for eosType=JMD95:'
430 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
431 & SQUEEZE_RIGHT , myThid )
432 ELSEIF ( kcheck .EQ. 2 ) THEN
433 WRITE(msgBuf,'(A)')
434 & 'EOS_CHECK: check values for eosType=MDJWF:'
435 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
436 & SQUEEZE_RIGHT , myThid )
437 ELSEIF ( kcheck .EQ. 6 ) THEN
438 WRITE(msgBuf,'(A)')
439 & 'EOS_CHECK: check values for eosType=UNESCO:'
440 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
441 & SQUEEZE_RIGHT , myThid )
442 ENDIF
443 WRITE(msgBuf,'(2(A,F4.1),A,F5.0,2A,F10.5,1X,F11.5)')
444 & 'rho(', sFld(i,j), ' PSU,',
445 & tFld(i,j), ' degC,',
446 & pLoc(kcheck)*SItoBar, ' bar)',
447 & ' = ', rLoc(kcheck), bLoc(kcheck)
448 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
449 & SQUEEZE_RIGHT , myThid )
450 WRITE(msgBuf,'(A,F10.5,1X,F11.5)')
451 & ' rho(find_rho_2d) = ',
452 & rhoLoc(i,j)+rhoConst, bulkMod(i,j)
453 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
454 & SQUEEZE_RIGHT , myThid )
455
456 CALL FIND_RHO_SCALAR( tFld(i,j), sLoc(kcheck),
457 & pLoc(kcheck), rhoLoc(i,j), myThid )
458 WRITE(msgBuf,'(A,F10.5,1X,F11.5)')
459 & ' rho(find_rho_scalar) = ',
460 & rhoLoc(i,j)
461 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
462 & SQUEEZE_RIGHT , myThid )
463
464 ENDDO
465 C end check nonlinear EOS
466 maskC(i,j,k,bi,bj) = mskSave
467 rC(k) = rC_Save
468
469 WRITE(msgBuf,'(A)') 'EOS_CHECK: Done'
470 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
471 & SQUEEZE_RIGHT , myThid )
472
473 ENDIF
474 #endif /* INCLUDE_EOS_CHECK */
475
476 RETURN
477 END

  ViewVC Help
Powered by ViewVC 1.1.22