/[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.5.6.1 - (show annotations) (download)
Fri Mar 7 23:10:20 2003 UTC (21 years, 3 months ago) by heimbach
Branch: ecco-branch
CVS Tags: ecco_c50_e32, ecco_c50_e33, ecco_c50_e30, ecco_c50_e31, ecco_c50_e29, ecco_c50_e28, ecco_c50_e33a
Changes since 1.5: +357 -55 lines
merging c49 and e27

1 C $Header: /u/gcmpack/MITgcm/model/src/find_alpha.F,v 1.12 2003/02/18 15:25:09 jmc 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 & kRef, theta, salt, rhoLoc, myThid )
216
217 CALL FIND_RHODEN( bi, bj, iMin, iMax, jMin, jMax, k,
218 & locPres, theta, salt, rhoDen, myThid )
219
220 DO j=jMin,jMax
221 DO i=iMin,iMax
222 t1 = theta(i,j,k,bi,bj)
223 t2 = t1*t1
224 s1 = salt(i,j,k,bi,bj)
225 sp5 = SQRT(s1)
226
227 p1 = locPres(i,j)*SItodBar
228 p1t1 = p1*t1
229
230 dnum_dtheta = eosMDJWFnum(1)
231 & + t1*(2.*eosMDJWFnum(2) + 3.*eosMDJWFnum(3)*t1)
232 & + eosMDJWFnum(5)*s1
233 & + p1t1*(2.*eosMDJWFnum(8) + 2.*eosMDJWFnum(11)*p1)
234
235 dden_dtheta = eosMDJWFden(1)
236 & + t1*(2.*eosMDJWFden(2)
237 & + t1*(3.*eosMDJWFden(3)
238 & + 4.*eosMDJWFden(4)*t1 ) )
239 & + s1*(eosMDJWFden(6)
240 & + t1*(3.*eosMDJWFden(7)*t1
241 & + 2.*eosMDJWFden(9)*sp5 ) )
242 & + p1*p1*(3.*eosMDJWFden(11)*t2 + eosMDJWFden(12)*p1)
243
244 alphaLoc(i,j) = rhoDen(i,j)*(dnum_dtheta
245 & - rhoLoc(i,j)*dden_dtheta)
246
247 ENDDO
248 ENDDO
249
250 ELSE
251 WRITE(*,*) 'FIND_ALPHA: equationOfState = ',equationOfState
252 STOP 'FIND_ALPHA: "equationOfState" has illegal value'
253 ENDIF
254
255 RETURN
256 END
257
258 SUBROUTINE FIND_BETA (
259 I bi, bj, iMin, iMax, jMin, jMax, k, kRef,
260 O betaloc )
261 C /==========================================================\
262 C | o SUBROUTINE FIND_BETA |
263 C | Calculates [drho(S,T,z) / dS] of a horizontal slice |
264 C |==========================================================|
265 C | |
266 C | k - is the Theta/Salt level |
267 C | kRef - determines pressure reference level |
268 C | (not used in 'LINEAR' mode) |
269 C | |
270 C | betaloc - drho / dS (kg/m^3/PSU) |
271 C | |
272 C \==========================================================/
273 IMPLICIT NONE
274
275 C === Global variables ===
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 == Routine arguments ==
284 C k :: Level of Theta/Salt slice
285 C kRef :: Pressure reference level
286 INTEGER bi,bj,iMin,iMax,jMin,jMax
287 INTEGER k
288 INTEGER kRef
289 _RL betaloc(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
290
291 C == Local variables ==
292 INTEGER i,j
293 _RL refTemp,refSalt,tP,sP
294 _RL t1, t2, t3, s1, s3o2, p1, sp5, p1t1
295 _RL drhoP0dS
296 _RL dKdS, dKdSSalt, dKdSPres
297 _RL locPres(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
298 _RL rhoP0 (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
299 _RL bulkMod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
300 _RL dnum_dsalt, dden_dsalt
301 _RL rhoDen (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
302 _RL rhoLoc (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
303 INTEGER myThid
304 CEOP
305
306 #ifdef CHECK_SALINITY_FOR_NEGATIVE_VALUES
307 CALL LOOK_FOR_NEG_SALINITY( bi, bj, iMin, iMax, jMin, jMax, k,
308 & sFld, myThid )
309 #endif
310
311 IF (equationOfState.EQ.'LINEAR') THEN
312
313 DO j=jMin,jMax
314 DO i=iMin,iMax
315 betaloc(i,j) = rhonil * sBeta
316 ENDDO
317 ENDDO
318
319 ELSEIF (equationOfState.EQ.'POLY3') THEN
320
321 refTemp=eosRefT(kRef)
322 refSalt=eosRefS(kRef)
323
324 DO j=jMin,jMax
325 DO i=iMin,iMax
326 tP=theta(i,j,k,bi,bj)-refTemp
327 sP=salt(i,j,k,bi,bj)-refSalt
328 #ifdef USE_FACTORIZED_POLY
329 betaloc(i,j) =
330 & ( eosC(9,kRef)*sP*3. + eosC(5,kRef)*2. )*sP + eosC(2,kRef)
331 & + ( eosC(7,kRef)*tP
332 & +eosC(8,kRef)*sP*2. + eosC(4,kRef)
333 & )*tP
334 #else
335 betaloc(i,j) =
336 & eosC(2,kRef) +
337 & eosC(4,kRef)*tP +
338 & eosC(5,kRef) *sP*2. +
339 & eosC(7,kRef)*tP*tP +
340 & eosC(8,kRef)*tP *sP*2. +
341 & eosC(9,kRef) *sP*sP*3.
342 #endif
343 ENDDO
344 ENDDO
345
346 ELSEIF ( equationOfState(1:5).EQ.'JMD95'
347 & .OR. equationOfState.EQ.'UNESCO' ) THEN
348 C nonlinear equation of state in pressure coordinates
349
350 CALL PRESSURE_FOR_EOS(
351 I bi, bj, iMin, iMax, jMin, jMax, kRef,
352 O locPres,
353 I myThid )
354
355 CALL FIND_RHOP0(
356 I bi, bj, iMin, iMax, jMin, jMax, k,
357 I theta, salt,
358 O rhoP0,
359 I myThid )
360
361 CALL FIND_BULKMOD(
362 I bi, bj, iMin, iMax, jMin, jMax, k,
363 I locPres, theta, salt,
364 O bulkMod,
365 I myThid )
366
367 DO j=jMin,jMax
368 DO i=iMin,iMax
369
370 C abbreviations
371 t1 = theta(i,j,k,bi,bj)
372 t2 = t1*t1
373 t3 = t2*t1
374
375 s1 = salt(i,j,k,bi,bj)
376 s3o2 = 1.5*SQRT(s1)
377
378 p1 = locPres(i,j)*SItoBar
379
380 C d(rho)/d(S)
381 C of fresh water at p = 0
382 drhoP0dS = 0. _d 0
383 C of salt water at p = 0
384 drhoP0dS = drhoP0dS
385 & + eosJMDCSw(1)
386 & + eosJMDCSw(2)*t1
387 & + eosJMDCSw(3)*t2
388 & + eosJMDCSw(4)*t3
389 & + eosJMDCSw(5)*t3*t1
390 & + s3o2*(
391 & eosJMDCSw(6)
392 & + eosJMDCSw(7)*t1
393 & + eosJMDCSw(8)*t2
394 & )
395 & + 2*eosJMDCSw(9)*s1
396 C d(bulk modulus)/d(S)
397 C of fresh water at p = 0
398 dKdS = 0. _d 0
399 C of sea water at p = 0
400 dKdSSalt =
401 & eosJMDCKSw(1)
402 & + eosJMDCKSw(2)*t1
403 & + eosJMDCKSw(3)*t2
404 & + eosJMDCKSw(4)*t3
405 & + s3o2*( eosJMDCKSw(5)
406 & + eosJMDCKSw(6)*t1
407 & + eosJMDCKSw(7)*t2
408 & )
409
410 C of sea water at p
411 dKdSPres =
412 & p1*( eosJMDCKP(5)
413 & + eosJMDCKP(6)*t1
414 & + eosJMDCKP(7)*t2
415 & )
416 & + s3o2*p1*eosJMDCKP(8)
417 & + p1*p1*( eosJMDCKP(12)
418 & + eosJMDCKP(13)*t1
419 & + eosJMDCKP(14)*t2
420 & )
421
422 dKdS = dKdSSalt + dKdSPres
423
424 betaloc(i,j) =
425 & ( bulkmod(i,j)**2*drhoP0dS
426 & - bulkmod(i,j)*p1*drhoP0dS
427 & - rhoP0(i,j)*p1*dKdS )
428 & /( bulkmod(i,j) - p1 )**2
429
430
431 ENDDO
432 ENDDO
433 ELSEIF ( equationOfState.EQ.'MDJWF' ) THEN
434
435 CALL PRESSURE_FOR_EOS(
436 I bi, bj, iMin, iMax, jMin, jMax, kRef,
437 O locPres,
438 I myThid )
439
440 CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k,
441 & kRef, theta, salt, rhoLoc, myThid )
442
443 CALL FIND_RHODEN( bi, bj, iMin, iMax, jMin, jMax, k,
444 & locPres, theta, salt, rhoDen, myThid )
445
446 DO j=jMin,jMax
447 DO i=iMin,iMax
448 t1 = theta(i,j,k,bi,bj)
449 t2 = t1*t1
450 s1 = salt(i,j,k,bi,bj)
451 sp5 = SQRT(s1)
452
453 p1 = locPres(i,j)*SItodBar
454 p1t1 = p1*t1
455
456 dnum_dsalt = eosMDJWFnum(4)
457 & + eosMDJWFnum(5)*t1
458 & + 2.*eosMDJWFnum(6)*s1 + eosMDJWFnum(9)*p1
459 dden_dsalt = eosMDJWFden(5)
460 & + t1*( eosMDJWFden(6) + eosMDJWFden(7)*t2 )
461 & + 1.5*sp5*(eosMDJWFden(8) + eosMDJWFden(9)*t2)
462
463 betaLoc(i,j) = rhoDen(i,j)*( dnum_dsalt
464 & - rhoLoc(i,j)*dden_dsalt )
465
466 ENDDO
467 ENDDO
468
469 ELSE
470 WRITE(*,*) 'FIND_BETA: equationOfState = ',equationOfState
471 STOP 'FIND_BETA: "equationOfState" has illegal value'
472 ENDIF
473
474 RETURN
475 END

  ViewVC Help
Powered by ViewVC 1.1.22