/[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.7 - (show annotations) (download)
Tue Oct 29 20:16:28 2002 UTC (21 years, 7 months ago) by mlosch
Branch: MAIN
CVS Tags: checkpoint46l_post, checkpoint46l_pre
Changes since 1.6: +78 -57 lines
o fixed a bug in ini_eos: changed e to _d for all MDJWF-coefficients

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

  ViewVC Help
Powered by ViewVC 1.1.22