/[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.15 - (show annotations) (download)
Tue Oct 17 18:52:34 2006 UTC (17 years, 7 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint59, checkpoint58y_post, checkpoint58t_post, checkpoint58w_post, checkpoint58q_post, mitgcm_mapl_00, checkpoint58r_post, checkpoint59e, checkpoint59d, checkpoint59g, checkpoint59f, checkpoint59a, checkpoint59c, checkpoint59b, checkpoint59m, checkpoint59l, checkpoint59i, checkpoint59h, checkpoint59k, checkpoint58v_post, checkpoint58x_post, checkpoint59j, checkpoint58u_post, checkpoint58s_post
Changes since 1.14: +118 -118 lines
clean-up multi-threaded problems (reported by debugger tcheck on ACES).

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

  ViewVC Help
Powered by ViewVC 1.1.22