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

Contents of /MITgcm/model/src/find_alpha.F

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


Revision 1.9 - (show annotations) (download)
Thu Sep 5 20:49:33 2002 UTC (21 years, 9 months ago) by mlosch
Branch: MAIN
CVS Tags: checkpoint46g_pre, checkpoint46f_post, checkpoint46e_post
Changes since 1.8: +156 -67 lines
o Added new equation of state -> MDJWF
  - EOS of McDougall et al., 2002, JAOT, submitted
  - caveat: the equation of state is only valid for a smaller (more
    realistic?) range of values than JMD95P/Z and UNESCO
  - added masks to the calculation of pressure in store_pressure
  - added more check values for density in check_eos (ini_eos.F), some of
    the old check values are out of the range of the MDJWF-eos, so don't
    expect perfect matches for those

1 C $Header: /u/gcmpack/MITgcm/model/src/find_alpha.F,v 1.8 2002/08/19 14:21:30 mlosch Exp $
2 C $Name: $
3
4 #include "CPP_OPTIONS.h"
5 #define USE_FACTORIZED_POLY
6
7 CBOP
8 C !ROUTINE: FIND_ALPHA
9 C !INTERFACE:
10 SUBROUTINE FIND_ALPHA (
11 I bi, bj, iMin, iMax, jMin, jMax, k, kRef, eqn,
12 O alphaloc )
13
14 C !DESCRIPTION: \bv
15 C *==========================================================*
16 C | o SUBROUTINE FIND_ALPHA
17 C | Calculates [drho(S,T,z) / dT] of a horizontal slice
18 C *==========================================================*
19 C |
20 C | k - is the Theta/Salt level
21 C | kRef - determines pressure reference level
22 C | (not used in 'LINEAR' mode)
23 C | eqn - determines the eqn. of state: 'LINEAR' or 'POLY3'
24 C |
25 C | alphaloc - drho / dT (kg/m^3/C)
26 C |
27 C *==========================================================*
28 C \ev
29
30 C !USES:
31 IMPLICIT NONE
32 c Common
33 #include "SIZE.h"
34 #include "DYNVARS.h"
35 #include "EEPARAMS.h"
36 #include "PARAMS.h"
37 #include "EOS.h"
38 #include "GRID.h"
39
40 C !INPUT/OUTPUT PARAMETERS:
41 c Arguments
42 integer bi,bj,iMin,iMax,jMin,jMax
43 integer k ! Level of Theta/Salt slice
44 integer kRef ! Pressure reference level
45 character*(*) eqn
46 _RL alphaloc(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
47
48 C !LOCAL VARIABLES:
49 c Local
50 integer i,j
51 _RL refTemp,refSalt,tP,sP
52 _RL t1, t2, t3, s1, s3o2, p1, p2, sp5, p1t1
53 _RL drhoP0dtheta, drhoP0dthetaFresh, drhoP0dthetaSalt
54 _RL dKdtheta, dKdthetaFresh, dKdthetaSalt, dKdthetaPres
55 _RL rhoP0(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
56 _RL bulkMod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
57 _RL dnum_dtheta, dden_dtheta
58 _RL rhoDen(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
59 _RL rhoLoc(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
60 integer myThid
61 CEOP
62
63 Cml stop 'myThid is not properly defined in this subroutine'
64
65 if (equationOfState.eq.'LINEAR') then
66
67 do j=jMin,jMax
68 do i=iMin,iMax
69 alphaloc(i,j) = -rhonil * tAlpha
70 enddo
71 enddo
72
73 elseif (equationOfState.eq.'POLY3') then
74
75 refTemp=eosRefT(kRef)
76 refSalt=eosRefS(kRef)
77
78 do j=jMin,jMax
79 do i=iMin,iMax
80 tP=theta(i,j,k,bi,bj)-refTemp
81 sP=salt(i,j,k,bi,bj)-refSalt
82 #ifdef USE_FACTORIZED_POLY
83 alphaloc(i,j) =
84 & ( eosC(6,kRef)
85 & *tP*3.
86 & +(eosC(7,kRef)*sP + eosC(3,kRef))*2.
87 & )*tP
88 & +(eosC(8,kRef)*sP + eosC(4,kRef) )*sP + eosC(1,kRef)
89 &
90 #else
91 alphaloc(i,j) =
92 & eosC(1,kRef) +
93 & eosC(3,kRef)*tP*2. +
94 & eosC(4,kRef) *sP +
95 & eosC(6,kRef)*tP*tP*3. +
96 & eosC(7,kRef)*tP*2. *sP +
97 & eosC(8,kRef) *sP*sP
98 #endif
99 enddo
100 enddo
101
102 elseif ( equationOfState(1:5).eq.'JMD95'
103 & .or. equationOfState.eq.'UNESCO' ) then
104 C nonlinear equation of state in pressure coordinates
105
106 CALL FIND_RHOP0(
107 I bi, bj, iMin, iMax, jMin, jMax, k,
108 I theta, salt,
109 O rhoP0,
110 I myThid )
111
112 CALL FIND_BULKMOD(
113 I bi, bj, iMin, iMax, jMin, jMax, k, kRef,
114 I theta, salt,
115 O bulkMod,
116 I myThid )
117
118 do j=jMin,jMax
119 do i=iMin,iMax
120
121 C abbreviations
122 t1 = theta(i,j,k,bi,bj)
123 t2 = t1*t1
124 t3 = t2*t1
125
126 s1 = salt(i,j,k,bi,bj)
127 if ( s1 .lt. 0. _d 0 ) then
128 C issue a warning
129 write(*,'(a,i3,a,i3,a,i3,a,e13.5)')
130 & ' FIND_ALPHA: WARNING, salinity(',
131 & i,',',j,',',k,',',bi,',',bj,') = ', s1
132 s1 = 0. _d 0
133 end if
134 s3o2 = sqrt(s1*s1*s1)
135
136 p1 = pressure(i,j,kRef,bi,bj)*SItoBar
137 p2 = p1*p1
138
139 C d(rho)/d(theta)
140 C of fresh water at p = 0
141 drhoP0dthetaFresh =
142 & eosJMDCFw(2)
143 & + 2.*eosJMDCFw(3)*t1
144 & + 3.*eosJMDCFw(4)*t2
145 & + 4.*eosJMDCFw(5)*t3
146 & + 5.*eosJMDCFw(6)*t3*t1
147 C of salt water at p = 0
148 drhoP0dthetaSalt =
149 & s1*(
150 & eosJMDCSw(2)
151 & + 2.*eosJMDCSw(3)*t1
152 & + 3.*eosJMDCSw(4)*t2
153 & + 4.*eosJMDCSw(5)*t3
154 & )
155 & + s3o2*(
156 & + eosJMDCSw(7)
157 & + 2.*eosJMDCSw(8)*t1
158 & )
159 C d(bulk modulus)/d(theta)
160 C of fresh water at p = 0
161 dKdthetaFresh =
162 & eosJMDCKFw(2)
163 & + 2.*eosJMDCKFw(3)*t1
164 & + 3.*eosJMDCKFw(4)*t2
165 & + 4.*eosJMDCKFw(5)*t3
166 C of sea water at p = 0
167 dKdthetaSalt =
168 & s1*( eosJMDCKSw(2)
169 & + 2.*eosJMDCKSw(3)*t1
170 & + 3.*eosJMDCKSw(4)*t2
171 & )
172 & + s3o2*( eosJMDCKSw(6)
173 & + 2.*eosJMDCKSw(7)*t1
174 & )
175 C of sea water at p
176 dKdthetaPres =
177 & p1*( eosJMDCKP(2)
178 & + 2.*eosJMDCKP(3)*t1
179 & + 3.*eosJMDCKP(4)*t2
180 & )
181 & + p1*s1*( eosJMDCKP(6)
182 & + 2.*eosJMDCKP(7)*t1
183 & )
184 & + p2*( eosJMDCKP(10)
185 & + 2.*eosJMDCKP(11)*t1
186 & )
187 & + p2*s1*( eosJMDCKP(13)
188 & + 2.*eosJMDCKP(14)*t1
189 & )
190
191 drhoP0dtheta = drhoP0dthetaFresh
192 & + drhoP0dthetaSalt
193 dKdtheta = dKdthetaFresh
194 & + dKdthetaSalt
195 & + dKdthetaPres
196 alphaloc(i,j) =
197 & ( bulkmod(i,j)**2*drhoP0dtheta
198 & - bulkmod(i,j)*p1*drhoP0dtheta
199 & - rhoP0(i,j)*p1*dKdtheta )
200 & /( bulkmod(i,j) - p1 )**2
201
202
203 end do
204 end do
205 elseif ( equationOfState.eq.'MDJWF' ) then
206
207 CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, kRef,
208 & eqn, theta, salt, rhoLoc, myThid )
209 CALL FIND_RHODEN( bi, bj, iMin, iMax, jMin, jMax, k, kRef,
210 & theta, salt, rhoDen, myThid )
211
212 do j=jMin,jMax
213 do i=iMin,iMax
214 t1 = theta(i,j,k,bi,bj)
215 t2 = t1*t1
216 s1 = salt(i,j,k,bi,bj)
217 if ( s1 .lt. 0. _d 0 ) then
218 C issue a warning
219 write(*,'(a,i3,a,i3,a,i3,a,e13.5)')
220 & ' FIND_ALPHA: WARNING, salinity(',
221 & i,',',j,',',k,',',bi,',',bj,') = ', s1
222 s1 = 0. _d 0
223 end if
224 sp5 = sqrt(s1)
225
226 p1 = pressure(i,j,kRef,bi,bj)*SItodBar
227 p1t1 = p1*t1
228
229 dnum_dtheta = eosMDJWFnum(1)
230 & + t1*(2.*eosMDJWFnum(2) + 3.*eosMDJWFnum(3)*t1)
231 & + eosMDJWFnum(5)*s1
232 & + p1t1*(2.*eosMDJWFnum(8) + 2.*eosMDJWFnum(11)*p1)
233
234 dden_dtheta = eosMDJWFden(1)
235 & + t1*(2.*eosMDJWFden(2)
236 & + t1*(3.*eosMDJWFden(3)
237 & + 4.*eosMDJWFden(4)*t1 ) )
238 & + s1*(eosMDJWFden(6)
239 & + t1*(3.*eosMDJWFden(7)*t1
240 & + 2.*eosMDJWFden(9)*sp5 ) )
241 & + p1*p1*(3.*eosMDJWFden(11)*t2 + eosMDJWFden(12)*p1)
242
243 alphaLoc(i,j) = rhoDen(i,j)*(dnum_dtheta
244 & - rhoLoc(i,j)*dden_dtheta)
245
246 end do
247 end do
248
249 else
250 write(*,*) 'FIND_ALPHA: equationOfState = ',equationOfState
251 stop 'FIND_ALPHA: "equationOfState" has illegal value'
252 endif
253
254 return
255 end
256
257 subroutine FIND_BETA (
258 I bi, bj, iMin, iMax, jMin, jMax, k, kRef, eqn,
259 O betaloc )
260 C /==========================================================\
261 C | o SUBROUTINE FIND_BETA |
262 C | Calculates [drho(S,T,z) / dS] of a horizontal slice |
263 C |==========================================================|
264 C | |
265 C | k - is the Theta/Salt level |
266 C | kRef - determines pressure reference level |
267 C | (not used in 'LINEAR' mode) |
268 C | eqn - determines the eqn. of state: 'LINEAR' or 'POLY3' |
269 C | |
270 C | betaloc - drho / dS (kg/m^3/PSU) |
271 C | |
272 C \==========================================================/
273 implicit none
274
275 c Common
276 #include "SIZE.h"
277 #include "DYNVARS.h"
278 #include "EEPARAMS.h"
279 #include "PARAMS.h"
280 #include "EOS.h"
281 #include "GRID.h"
282
283 c Arguments
284 integer bi,bj,iMin,iMax,jMin,jMax
285 integer k ! Level of Theta/Salt slice
286 integer kRef ! Pressure reference level
287 character*(*) eqn
288 _RL betaloc(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
289
290 c Local
291 integer i,j
292 _RL refTemp,refSalt,tP,sP
293 _RL t1, t2, t3, s1, s3o2, p1, sp5, p1t1
294 _RL drhoP0dS
295 _RL dKdS, dKdSSalt, dKdSPres
296 _RL rhoP0(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
297 _RL bulkMod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
298 _RL dnum_dsalt, dden_dsalt
299 _RL rhoDen(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
300 _RL rhoLoc(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
301 integer myThid
302 CEOP
303
304 Cml stop 'myThid is not properly defined in this subroutine'
305
306 if (equationOfState.eq.'LINEAR') then
307
308 do j=jMin,jMax
309 do i=iMin,iMax
310 betaloc(i,j) = rhonil * sBeta
311 enddo
312 enddo
313
314 elseif (equationOfState.eq.'POLY3') then
315
316 refTemp=eosRefT(kRef)
317 refSalt=eosRefS(kRef)
318
319 do j=jMin,jMax
320 do i=iMin,iMax
321 tP=theta(i,j,k,bi,bj)-refTemp
322 sP=salt(i,j,k,bi,bj)-refSalt
323 #ifdef USE_FACTORIZED_POLY
324 betaloc(i,j) =
325 & ( eosC(9,kRef)*sP*3. + eosC(5,kRef)*2. )*sP + eosC(2,kRef)
326 & + ( eosC(7,kRef)*tP
327 & +eosC(8,kRef)*sP*2. + eosC(4,kRef)
328 & )*tP
329 #else
330 betaloc(i,j) =
331 & eosC(2,kRef) +
332 & eosC(4,kRef)*tP +
333 & eosC(5,kRef) *sP*2. +
334 & eosC(7,kRef)*tP*tP +
335 & eosC(8,kRef)*tP *sP*2. +
336 & eosC(9,kRef) *sP*sP*3.
337 #endif
338 enddo
339 enddo
340
341 elseif ( equationOfState(1:5).eq.'JMD95'
342 & .or. equationOfState.eq.'UNESCO' ) then
343 C nonlinear equation of state in pressure coordinates
344
345 CALL FIND_RHOP0(
346 I bi, bj, iMin, iMax, jMin, jMax, k,
347 I theta, salt,
348 O rhoP0,
349 I myThid )
350
351 CALL FIND_BULKMOD(
352 I bi, bj, iMin, iMax, jMin, jMax, k, kRef,
353 I theta, salt,
354 O bulkMod,
355 I myThid )
356
357 do j=jMin,jMax
358 do i=iMin,iMax
359
360 C abbreviations
361 t1 = theta(i,j,k,bi,bj)
362 t2 = t1*t1
363 t3 = t2*t1
364
365 s1 = salt(i,j,k,bi,bj)
366 if ( s1 .lt. 0. _d 0 ) then
367 C issue a warning
368 write(*,'(a,i3,a,i3,a,i3,a,e13.5)')
369 & ' FIND_BETA: WARNING, salinity(',
370 & i,',',j,',',k,',',bi,',',bj,') = ', s1
371 s1 = 0. _d 0
372 end if
373 s3o2 = 1.5*sqrt(s1)
374
375 p1 = pressure(i,j,kRef,bi,bj)*SItoBar
376
377 C d(rho)/d(S)
378 C of fresh water at p = 0
379 drhoP0dS = 0. _d 0
380 C of salt water at p = 0
381 drhoP0dS = drhoP0dS
382 & + eosJMDCSw(1)
383 & + eosJMDCSw(2)*t1
384 & + eosJMDCSw(3)*t2
385 & + eosJMDCSw(4)*t3
386 & + eosJMDCSw(5)*t3*t1
387 & + s3o2*(
388 & eosJMDCSw(6)
389 & + eosJMDCSw(7)*t1
390 & + eosJMDCSw(8)*t2
391 & )
392 & + 2*eosJMDCSw(9)*s1
393 C d(bulk modulus)/d(S)
394 C of fresh water at p = 0
395 dKdS = 0. _d 0
396 C of sea water at p = 0
397 dKdSSalt =
398 & eosJMDCKSw(1)
399 & + eosJMDCKSw(2)*t1
400 & + eosJMDCKSw(3)*t2
401 & + eosJMDCKSw(4)*t3
402 & + s3o2*( eosJMDCKSw(5)
403 & + eosJMDCKSw(6)*t1
404 & + eosJMDCKSw(7)*t2
405 & )
406
407 C of sea water at p
408 dKdSPres =
409 & p1*( eosJMDCKP(5)
410 & + eosJMDCKP(6)*t1
411 & + eosJMDCKP(7)*t2
412 & )
413 & + s3o2*p1*eosJMDCKP(8)
414 & + p1*p1*( eosJMDCKP(12)
415 & + eosJMDCKP(13)*t1
416 & + eosJMDCKP(14)*t2
417 & )
418
419 dKdS = dKdSSalt + dKdSPres
420
421 betaloc(i,j) =
422 & ( bulkmod(i,j)**2*drhoP0dS
423 & - bulkmod(i,j)*p1*drhoP0dS
424 & - rhoP0(i,j)*p1*dKdS )
425 & /( bulkmod(i,j) - p1 )**2
426
427
428 end do
429 end do
430 elseif ( equationOfState.eq.'MDJWF' ) then
431
432 CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, kRef,
433 & eqn, theta, salt, rhoLoc, myThid )
434 CALL FIND_RHODEN( bi, bj, iMin, iMax, jMin, jMax, k, kRef,
435 & theta, salt, rhoDen, myThid )
436
437 do j=jMin,jMax
438 do i=iMin,iMax
439 t1 = theta(i,j,k,bi,bj)
440 t2 = t1*t1
441 s1 = salt(i,j,k,bi,bj)
442 if ( s1 .lt. 0. _d 0 ) then
443 C issue a warning
444 write(*,'(a,i3,a,i3,a,i3,a,e13.5)')
445 & ' FIND_BETA: WARNING, salinity(',
446 & i,',',j,',',k,',',bi,',',bj,') = ', s1
447 s1 = 0. _d 0
448 end if
449 sp5 = sqrt(s1)
450
451 p1 = pressure(i,j,k,bi,bj)*SItodBar
452 p1t1 = p1*t1
453
454 dnum_dsalt = eosMDJWFnum(4)
455 & + eosMDJWFnum(5)*t1
456 & + 2.*eosMDJWFnum(6)*s1 + eosMDJWFnum(9)*p1
457 dden_dsalt = eosMDJWFden(5)
458 & + t1*( eosMDJWFden(6) + eosMDJWFden(7)*t2 )
459 & + 1.5*sp5*(eosMDJWFden(8) + eosMDJWFden(9)*t2)
460
461 betaLoc(i,j) = rhoDen(i,j)*( dnum_dsalt
462 & - rhoLoc(i,j)*dden_dsalt )
463
464 end do
465 end do
466
467 else
468 write(*,*) 'FIND_BETA: equationOfState = ',equationOfState
469 stop 'FIND_BETA: "equationOfState" has illegal value'
470 endif
471
472 return
473 end

  ViewVC Help
Powered by ViewVC 1.1.22