/[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.3 - (hide annotations) (download)
Tue Aug 20 14:45:18 2002 UTC (21 years, 9 months ago) by mlosch
Branch: MAIN
Changes since 1.2: +2 -2 lines
o bug-fix in ini_eos.F: check_eos still called find_rhop0 with the wrong
  number of arguments

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

  ViewVC Help
Powered by ViewVC 1.1.22