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

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

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


Revision 1.7 - (hide 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 mlosch 1.7 C $Header: /u/gcmpack/MITgcm/model/src/ini_eos.F,v 1.6 2002/09/25 19:36:50 mlosch Exp $
2 mlosch 1.2 C $Name: $
3 mlosch 1.1
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 mlosch 1.4 #include "DYNVARS.h"
27 mlosch 1.1
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 mlosch 1.4 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 mlosch 1.7 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 mlosch 1.1 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 mlosch 1.2
139 mlosch 1.1 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 mlosch 1.2 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 mlosch 1.1
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 mlosch 1.5 & 'WARNING: (see Jackett and McDougall 1995, JAOT)'
227 mlosch 1.1 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 mlosch 1.4 elseif ( equationOfState .eq. 'MDJWF' ) then
268    
269 mlosch 1.7 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 mlosch 1.4
282    
283 mlosch 1.7 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 mlosch 1.4
297 mlosch 1.1 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 mlosch 1.4 _RL rhoLoc (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
347 mlosch 1.1 _RL bulkMod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
348 mlosch 1.4 _RL psave
349 mlosch 1.1
350     INTEGER ncheck, kcheck
351 mlosch 1.4 PARAMETER ( ncheck = 13 )
352 mlosch 1.1 _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 mlosch 1.4 & 25.44820830309568 _d 0, 20.17368557065936 _d 0,
360     & 13.43397459640398 _d 0,
361 mlosch 1.1 & 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 mlosch 1.4 & 25. _d 0, 20. _d 0,
368     & 12. _d 0,
369 mlosch 1.1 & 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 mlosch 1.4 & 35.0 _d 0, 20. _d 0,
376     & 40.0 _d 0,
377 mlosch 1.1 & 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 mlosch 1.4 & 200. _d 5, 100. _d 5,
384     & 800. _d 5,
385 mlosch 1.1 & 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 mlosch 1.4 & /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 mlosch 1.1 & bloc
398     & / -1.00000 _d 0, -1.00000 _d 0,
399 mlosch 1.4 & -1.00000 _d 0, -1.00000 _d 0,
400     & -1.00000 _d 0,
401 mlosch 1.1 & 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 mlosch 1.4 & 'check_eos: Check the equation of state: Type ',
421 mlosch 1.1 & 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 mlosch 1.4 rholoc(i,j) = 0. _d 0
434 mlosch 1.1 bulkMod(i,j) = -1. _d 0
435    
436 mlosch 1.4 call find_rho(
437     & bi, bj, iMin, iMax, jMin, jMax, k, k,
438     & tFld, sFld, rholoc, myThid )
439    
440 mlosch 1.1 call find_bulkmod(
441 mlosch 1.4 & bi, bj, imin, imax, jmin, jmax, k, k,
442     & tFld, sFld, bulkMod, myThid )
443 mlosch 1.5
444 mlosch 1.1 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 mlosch 1.5 write(msgBuf,'(a14,a22,f10.5,1x,f11.5)')
453     & 'rho(find_rho) ',
454 mlosch 1.6 & ' = ', rholoc(i,j)+rhoConst, bulkMod(i,j)
455 mlosch 1.5 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 mlosch 1.6 & ' = ', rholoc(i,j)+rhoConst, bulkMod(i,j)
464 mlosch 1.1 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 mlosch 1.4

  ViewVC Help
Powered by ViewVC 1.1.22