/[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.5 - (show annotations) (download)
Wed Sep 18 16:38:02 2002 UTC (21 years, 8 months ago) by mlosch
Branch: MAIN
CVS Tags: checkpoint46h_pre, checkpoint46g_post
Changes since 1.4: +15 -7 lines
o Include a new diagnostic variable phiHydLow for the ocean model
  - in z-coordinates, it is the bottom pressure anomaly
  - in p-coordinates, it is the sea surface elevation
  - in both cases, these variable have global drift, reflecting the mass
    drift in z-coordinates and the volume drift in p-coordinates
  - included time averaging for phiHydLow, be aware of the drift!
o depth-dependent computation of Bo_surf for pressure coordinates
  in the ocean (buoyancyRelation='OCEANICP')
  - requires a new routine (FIND_RHO_SCALAR) to compute density with only
    Theta, Salinity, and Pressure in the parameter list. This routine is
    presently contained in find_rho.F. This routine does not give the
    correct density for 'POLY3', which would be a z-dependent reference
    density.
o cleaned up find_rho
  - removed obsolete 'eqn' from the parameter list.
o added two new verification experiments: gop and goz
  (4x4 degree global ocean, 15 layers in pressure and height coordinates)

1 C $Header: /u/gcmpack/MITgcm/model/src/ini_eos.F,v 1.4 2002/09/05 20:49:33 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 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
87 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 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 pressure(i,j,k,bi,bj) = rhonil * (
116 & - gravity*rC(k)
117 & )
118 end do
119 end do
120 end do
121 end do
122 end do
123 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
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 & 'WARNING: (see Jackett and McDougall 1995, JAOT)'
206 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 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 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 _RL rhoLoc (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
326 _RL bulkMod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
327 _RL psave
328
329 INTEGER ncheck, kcheck
330 PARAMETER ( ncheck = 13 )
331 _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 & 25.44820830309568 _d 0, 20.17368557065936 _d 0,
339 & 13.43397459640398 _d 0,
340 & 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 & 25. _d 0, 20. _d 0,
347 & 12. _d 0,
348 & 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 & 35.0 _d 0, 20. _d 0,
355 & 40.0 _d 0,
356 & 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 & 200. _d 5, 100. _d 5,
363 & 800. _d 5,
364 & 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 & /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 & bloc
377 & / -1.00000 _d 0, -1.00000 _d 0,
378 & -1.00000 _d 0, -1.00000 _d 0,
379 & -1.00000 _d 0,
380 & 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 & 'check_eos: Check the equation of state: Type ',
400 & 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 rholoc(i,j) = 0. _d 0
413 bulkMod(i,j) = -1. _d 0
414
415 call find_rho(
416 & bi, bj, iMin, iMax, jMin, jMax, k, k,
417 & tFld, sFld, rholoc, myThid )
418
419 call find_bulkmod(
420 & bi, bj, imin, imax, jmin, jmax, k, k,
421 & tFld, sFld, bulkMod, myThid )
422
423 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 write(msgBuf,'(a14,a22,f10.5,1x,f11.5)')
432 & 'rho(find_rho) ',
433 & ' = ', rholoc(i,j)+rhoNil, bulkMod(i,j)
434 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 & ' = ', rholoc(i,j)+rhoNil, bulkMod(i,j)
443 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

  ViewVC Help
Powered by ViewVC 1.1.22