/[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.11 - (hide annotations) (download)
Tue Oct 19 02:39:58 2004 UTC (19 years, 7 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint57b_post, checkpoint56b_post, checkpoint57, checkpoint56, checkpoint55i_post, checkpoint57a_post, checkpoint55j_post, checkpoint55h_post, checkpoint56a_post, checkpoint56c_post, checkpoint57a_pre
Changes since 1.10: +3 -4 lines
use flags: fluidIsAir/Water, usingP/ZCoords instead of buoyancyRelation

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

  ViewVC Help
Powered by ViewVC 1.1.22