/[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.6 - (hide annotations) (download)
Wed Sep 25 19:36:50 2002 UTC (21 years, 8 months ago) by mlosch
Branch: MAIN
CVS Tags: checkpoint46j_pre, checkpoint46j_post, checkpoint46k_post, checkpoint46i_post, checkpoint46h_post
Changes since 1.5: +4 -4 lines
o cleaned up the use of rhoNil and rhoConst.
  - rhoNil should only appear in the LINEAR equation of state, everywhere
    else rhoNil is replaced by rhoConst, e.g. find_rho computes rho-rhoConst
    and the dynamical equations are all divided by rhoConst
o introduced new parameter rhoConstFresh, a reference density of fresh
  water, to remove the fresh water flux's dependence on rhoNil. The default
  value is 999.8 kg/m^3
o cleanup up external_forcing.F and external_forcing_surf.F
  - can now be used by both OCEANIC and OCEANICP

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

  ViewVC Help
Powered by ViewVC 1.1.22