/[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.12 - (show annotations) (download)
Tue Feb 18 15:25:09 2003 UTC (21 years, 3 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint48i_post, checkpoint48h_post, checkpoint49, checkpoint48g_post
Changes since 1.11: +127 -105 lines
o compute locally the pressure for use in EOS: UNESCO, JMD95P or MDJWF
o store total Potential in totPhyHyd for diagnostic & EOS funct. of P
o fix restart and overlap Pb when using Z-coord and EOS funct. of P

1 C $Header: /u/gcmpack/MITgcm/model/src/find_alpha.F,v 1.11 2002/11/01 22:00:33 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,
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 |
24 C | alphaloc - drho / dT (kg/m^3/C)
25 C |
26 C *==========================================================*
27 C \ev
28
29 C !USES:
30 IMPLICIT NONE
31 C === Global variables ===
32 #include "SIZE.h"
33 #include "DYNVARS.h"
34 #include "EEPARAMS.h"
35 #include "PARAMS.h"
36 #include "EOS.h"
37 #include "GRID.h"
38
39 C !INPUT/OUTPUT PARAMETERS:
40 C == Routine arguments ==
41 C k :: Level of Theta/Salt slice
42 C kRef :: Pressure reference level
43 INTEGER bi,bj,iMin,iMax,jMin,jMax
44 INTEGER k
45 INTEGER kRef
46 _RL alphaloc(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
47
48 C !LOCAL VARIABLES:
49 C == Local variables ==
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 locPres(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
56 _RL rhoP0 (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
57 _RL bulkMod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
58 _RL dnum_dtheta, dden_dtheta
59 _RL rhoDen (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
60 _RL rhoLoc (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
61 INTEGER myThid
62 CEOP
63
64 #ifdef CHECK_SALINITY_FOR_NEGATIVE_VALUES
65 CALL LOOK_FOR_NEG_SALINITY( bi, bj, iMin, iMax, jMin, jMax, k,
66 & sFld, myThid )
67 #endif
68
69 IF (equationOfState.EQ.'LINEAR') THEN
70
71 DO j=jMin,jMax
72 DO i=iMin,iMax
73 alphaloc(i,j) = -rhonil * tAlpha
74 ENDDO
75 ENDDO
76
77 ELSEIF (equationOfState.EQ.'POLY3') THEN
78
79 refTemp=eosRefT(kRef)
80 refSalt=eosRefS(kRef)
81
82 DO j=jMin,jMax
83 DO i=iMin,iMax
84 tP=theta(i,j,k,bi,bj)-refTemp
85 sP=salt(i,j,k,bi,bj)-refSalt
86 #ifdef USE_FACTORIZED_POLY
87 alphaloc(i,j) =
88 & ( eosC(6,kRef)
89 & *tP*3.
90 & +(eosC(7,kRef)*sP + eosC(3,kRef))*2.
91 & )*tP
92 & +(eosC(8,kRef)*sP + eosC(4,kRef) )*sP + eosC(1,kRef)
93 &
94 #else
95 alphaloc(i,j) =
96 & eosC(1,kRef) +
97 & eosC(3,kRef)*tP*2. +
98 & eosC(4,kRef) *sP +
99 & eosC(6,kRef)*tP*tP*3. +
100 & eosC(7,kRef)*tP*2. *sP +
101 & eosC(8,kRef) *sP*sP
102 #endif
103 ENDDO
104 ENDDO
105
106 ELSEIF ( equationOfState(1:5).EQ.'JMD95'
107 & .OR. equationOfState.EQ.'UNESCO' ) THEN
108 C nonlinear equation of state in pressure coordinates
109
110 CALL PRESSURE_FOR_EOS(
111 I bi, bj, iMin, iMax, jMin, jMax, kRef,
112 O locPres,
113 I myThid )
114
115 CALL FIND_RHOP0(
116 I bi, bj, iMin, iMax, jMin, jMax, k,
117 I theta, salt,
118 O rhoP0,
119 I myThid )
120
121 CALL FIND_BULKMOD(
122 I bi, bj, iMin, iMax, jMin, jMax, k,
123 I locPres, theta, salt,
124 O bulkMod,
125 I myThid )
126
127 DO j=jMin,jMax
128 DO i=iMin,iMax
129
130 C abbreviations
131 t1 = theta(i,j,k,bi,bj)
132 t2 = t1*t1
133 t3 = t2*t1
134
135 s1 = salt(i,j,k,bi,bj)
136 s3o2 = SQRT(s1*s1*s1)
137
138 p1 = locPres(i,j)*SItoBar
139 p2 = p1*p1
140
141 C d(rho)/d(theta)
142 C of fresh water at p = 0
143 drhoP0dthetaFresh =
144 & eosJMDCFw(2)
145 & + 2.*eosJMDCFw(3)*t1
146 & + 3.*eosJMDCFw(4)*t2
147 & + 4.*eosJMDCFw(5)*t3
148 & + 5.*eosJMDCFw(6)*t3*t1
149 C of salt water at p = 0
150 drhoP0dthetaSalt =
151 & s1*(
152 & eosJMDCSw(2)
153 & + 2.*eosJMDCSw(3)*t1
154 & + 3.*eosJMDCSw(4)*t2
155 & + 4.*eosJMDCSw(5)*t3
156 & )
157 & + s3o2*(
158 & + eosJMDCSw(7)
159 & + 2.*eosJMDCSw(8)*t1
160 & )
161 C d(bulk modulus)/d(theta)
162 C of fresh water at p = 0
163 dKdthetaFresh =
164 & eosJMDCKFw(2)
165 & + 2.*eosJMDCKFw(3)*t1
166 & + 3.*eosJMDCKFw(4)*t2
167 & + 4.*eosJMDCKFw(5)*t3
168 C of sea water at p = 0
169 dKdthetaSalt =
170 & s1*( eosJMDCKSw(2)
171 & + 2.*eosJMDCKSw(3)*t1
172 & + 3.*eosJMDCKSw(4)*t2
173 & )
174 & + s3o2*( eosJMDCKSw(6)
175 & + 2.*eosJMDCKSw(7)*t1
176 & )
177 C of sea water at p
178 dKdthetaPres =
179 & p1*( eosJMDCKP(2)
180 & + 2.*eosJMDCKP(3)*t1
181 & + 3.*eosJMDCKP(4)*t2
182 & )
183 & + p1*s1*( eosJMDCKP(6)
184 & + 2.*eosJMDCKP(7)*t1
185 & )
186 & + p2*( eosJMDCKP(10)
187 & + 2.*eosJMDCKP(11)*t1
188 & )
189 & + p2*s1*( eosJMDCKP(13)
190 & + 2.*eosJMDCKP(14)*t1
191 & )
192
193 drhoP0dtheta = drhoP0dthetaFresh
194 & + drhoP0dthetaSalt
195 dKdtheta = dKdthetaFresh
196 & + dKdthetaSalt
197 & + dKdthetaPres
198 alphaloc(i,j) =
199 & ( bulkmod(i,j)**2*drhoP0dtheta
200 & - bulkmod(i,j)*p1*drhoP0dtheta
201 & - rhoP0(i,j)*p1*dKdtheta )
202 & /( bulkmod(i,j) - p1 )**2
203
204
205 ENDDO
206 ENDDO
207 ELSEIF ( equationOfState.EQ.'MDJWF' ) THEN
208
209 CALL PRESSURE_FOR_EOS(
210 I bi, bj, iMin, iMax, jMin, jMax, kRef,
211 O locPres,
212 I myThid )
213
214 CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k,
215 & locPres, theta, salt, rhoLoc, myThid )
216 CALL FIND_RHODEN( bi, bj, iMin, iMax, jMin, jMax, k,
217 & locPres, theta, salt, rhoDen, myThid )
218
219 DO j=jMin,jMax
220 DO i=iMin,iMax
221 t1 = theta(i,j,k,bi,bj)
222 t2 = t1*t1
223 s1 = salt(i,j,k,bi,bj)
224 sp5 = SQRT(s1)
225
226 p1 = locPres(i,j)*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 ENDDO
247 ENDDO
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,
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 | |
269 C | betaloc - drho / dS (kg/m^3/PSU) |
270 C | |
271 C \==========================================================/
272 IMPLICIT NONE
273
274 C === Global variables ===
275 #include "SIZE.h"
276 #include "DYNVARS.h"
277 #include "EEPARAMS.h"
278 #include "PARAMS.h"
279 #include "EOS.h"
280 #include "GRID.h"
281
282 C == Routine arguments ==
283 C k :: Level of Theta/Salt slice
284 C kRef :: Pressure reference level
285 INTEGER bi,bj,iMin,iMax,jMin,jMax
286 INTEGER k
287 INTEGER kRef
288 _RL betaloc(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
289
290 C == Local variables ==
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 locPres(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
297 _RL rhoP0 (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
298 _RL bulkMod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
299 _RL dnum_dsalt, dden_dsalt
300 _RL rhoDen (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
301 _RL rhoLoc (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
302 INTEGER myThid
303 CEOP
304
305 #ifdef CHECK_SALINITY_FOR_NEGATIVE_VALUES
306 CALL LOOK_FOR_NEG_SALINITY( bi, bj, iMin, iMax, jMin, jMax, k,
307 & sFld, myThid )
308 #endif
309
310 IF (equationOfState.EQ.'LINEAR') THEN
311
312 DO j=jMin,jMax
313 DO i=iMin,iMax
314 betaloc(i,j) = rhonil * sBeta
315 ENDDO
316 ENDDO
317
318 ELSEIF (equationOfState.EQ.'POLY3') THEN
319
320 refTemp=eosRefT(kRef)
321 refSalt=eosRefS(kRef)
322
323 DO j=jMin,jMax
324 DO i=iMin,iMax
325 tP=theta(i,j,k,bi,bj)-refTemp
326 sP=salt(i,j,k,bi,bj)-refSalt
327 #ifdef USE_FACTORIZED_POLY
328 betaloc(i,j) =
329 & ( eosC(9,kRef)*sP*3. + eosC(5,kRef)*2. )*sP + eosC(2,kRef)
330 & + ( eosC(7,kRef)*tP
331 & +eosC(8,kRef)*sP*2. + eosC(4,kRef)
332 & )*tP
333 #else
334 betaloc(i,j) =
335 & eosC(2,kRef) +
336 & eosC(4,kRef)*tP +
337 & eosC(5,kRef) *sP*2. +
338 & eosC(7,kRef)*tP*tP +
339 & eosC(8,kRef)*tP *sP*2. +
340 & eosC(9,kRef) *sP*sP*3.
341 #endif
342 ENDDO
343 ENDDO
344
345 ELSEIF ( equationOfState(1:5).EQ.'JMD95'
346 & .OR. equationOfState.EQ.'UNESCO' ) THEN
347 C nonlinear equation of state in pressure coordinates
348
349 CALL PRESSURE_FOR_EOS(
350 I bi, bj, iMin, iMax, jMin, jMax, kRef,
351 O locPres,
352 I myThid )
353
354 CALL FIND_RHOP0(
355 I bi, bj, iMin, iMax, jMin, jMax, k,
356 I theta, salt,
357 O rhoP0,
358 I myThid )
359
360 CALL FIND_BULKMOD(
361 I bi, bj, iMin, iMax, jMin, jMax, k,
362 I locPres, theta, salt,
363 O bulkMod,
364 I myThid )
365
366 DO j=jMin,jMax
367 DO i=iMin,iMax
368
369 C abbreviations
370 t1 = theta(i,j,k,bi,bj)
371 t2 = t1*t1
372 t3 = t2*t1
373
374 s1 = salt(i,j,k,bi,bj)
375 s3o2 = 1.5*SQRT(s1)
376
377 p1 = locPres(i,j)*SItoBar
378
379 C d(rho)/d(S)
380 C of fresh water at p = 0
381 drhoP0dS = 0. _d 0
382 C of salt water at p = 0
383 drhoP0dS = drhoP0dS
384 & + eosJMDCSw(1)
385 & + eosJMDCSw(2)*t1
386 & + eosJMDCSw(3)*t2
387 & + eosJMDCSw(4)*t3
388 & + eosJMDCSw(5)*t3*t1
389 & + s3o2*(
390 & eosJMDCSw(6)
391 & + eosJMDCSw(7)*t1
392 & + eosJMDCSw(8)*t2
393 & )
394 & + 2*eosJMDCSw(9)*s1
395 C d(bulk modulus)/d(S)
396 C of fresh water at p = 0
397 dKdS = 0. _d 0
398 C of sea water at p = 0
399 dKdSSalt =
400 & eosJMDCKSw(1)
401 & + eosJMDCKSw(2)*t1
402 & + eosJMDCKSw(3)*t2
403 & + eosJMDCKSw(4)*t3
404 & + s3o2*( eosJMDCKSw(5)
405 & + eosJMDCKSw(6)*t1
406 & + eosJMDCKSw(7)*t2
407 & )
408
409 C of sea water at p
410 dKdSPres =
411 & p1*( eosJMDCKP(5)
412 & + eosJMDCKP(6)*t1
413 & + eosJMDCKP(7)*t2
414 & )
415 & + s3o2*p1*eosJMDCKP(8)
416 & + p1*p1*( eosJMDCKP(12)
417 & + eosJMDCKP(13)*t1
418 & + eosJMDCKP(14)*t2
419 & )
420
421 dKdS = dKdSSalt + dKdSPres
422
423 betaloc(i,j) =
424 & ( bulkmod(i,j)**2*drhoP0dS
425 & - bulkmod(i,j)*p1*drhoP0dS
426 & - rhoP0(i,j)*p1*dKdS )
427 & /( bulkmod(i,j) - p1 )**2
428
429
430 ENDDO
431 ENDDO
432 ELSEIF ( equationOfState.EQ.'MDJWF' ) THEN
433
434 CALL PRESSURE_FOR_EOS(
435 I bi, bj, iMin, iMax, jMin, jMax, kRef,
436 O locPres,
437 I myThid )
438
439 CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k,
440 & locPres, theta, salt, rhoLoc, myThid )
441 CALL FIND_RHODEN( bi, bj, iMin, iMax, jMin, jMax, k,
442 & locPres, theta, salt, rhoDen, myThid )
443
444 DO j=jMin,jMax
445 DO i=iMin,iMax
446 t1 = theta(i,j,k,bi,bj)
447 t2 = t1*t1
448 s1 = salt(i,j,k,bi,bj)
449 sp5 = SQRT(s1)
450
451 p1 = locPres(i,j)*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 ENDDO
465 ENDDO
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