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

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

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


Revision 1.13 - (hide annotations) (download)
Fri Mar 7 23:46:46 2003 UTC (21 years, 3 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint50c_post, checkpoint50c_pre, checkpoint51, checkpoint50, checkpoint50d_post, checkpoint50b_pre, checkpoint51d_post, checkpoint51b_pre, checkpoint50f_post, checkpoint50a_post, checkpoint50f_pre, checkpoint51b_post, checkpoint51c_post, checkpoint50g_post, checkpoint50h_post, checkpoint50e_pre, checkpoint50i_post, checkpoint50e_post, checkpoint50d_pre, checkpoint51e_post, checkpoint51f_pre, checkpoint50b_post, checkpoint51a_post
Changes since 1.12: +5 -3 lines
bug fix for call find_rho

1 heimbach 1.13 C $Header: /u/gcmpack/MITgcm/model/src/find_alpha.F,v 1.12 2003/02/18 15:25:09 jmc Exp $
2 cnh 1.5 C $Name: $
3 adcroft 1.1
4     #include "CPP_OPTIONS.h"
5     #define USE_FACTORIZED_POLY
6    
7 cnh 1.5 CBOP
8     C !ROUTINE: FIND_ALPHA
9     C !INTERFACE:
10     SUBROUTINE FIND_ALPHA (
11 mlosch 1.10 I bi, bj, iMin, iMax, jMin, jMax, k, kRef,
12 adcroft 1.1 O alphaloc )
13    
14 cnh 1.5 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 jmc 1.12 C === Global variables ===
32 adcroft 1.1 #include "SIZE.h"
33     #include "DYNVARS.h"
34     #include "EEPARAMS.h"
35     #include "PARAMS.h"
36 mlosch 1.6 #include "EOS.h"
37     #include "GRID.h"
38 adcroft 1.1
39 cnh 1.5 C !INPUT/OUTPUT PARAMETERS:
40 jmc 1.12 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 adcroft 1.1 _RL alphaloc(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
47    
48 cnh 1.5 C !LOCAL VARIABLES:
49 jmc 1.12 C == Local variables ==
50     INTEGER i,j
51 adcroft 1.1 _RL refTemp,refSalt,tP,sP
52 mlosch 1.9 _RL t1, t2, t3, s1, s3o2, p1, p2, sp5, p1t1
53 mlosch 1.6 _RL drhoP0dtheta, drhoP0dthetaFresh, drhoP0dthetaSalt
54     _RL dKdtheta, dKdthetaFresh, dKdthetaSalt, dKdthetaPres
55 jmc 1.12 _RL locPres(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
56     _RL rhoP0 (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
57 mlosch 1.6 _RL bulkMod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
58 mlosch 1.9 _RL dnum_dtheta, dden_dtheta
59 jmc 1.12 _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 cnh 1.5 CEOP
63 adcroft 1.1
64 mlosch 1.11 #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 jmc 1.12 IF (equationOfState.EQ.'LINEAR') THEN
70 adcroft 1.1
71 jmc 1.12 DO j=jMin,jMax
72     DO i=iMin,iMax
73 adcroft 1.1 alphaloc(i,j) = -rhonil * tAlpha
74 jmc 1.12 ENDDO
75     ENDDO
76 adcroft 1.1
77 jmc 1.12 ELSEIF (equationOfState.EQ.'POLY3') THEN
78 adcroft 1.1
79     refTemp=eosRefT(kRef)
80     refSalt=eosRefS(kRef)
81    
82 jmc 1.12 DO j=jMin,jMax
83     DO i=iMin,iMax
84 adcroft 1.1 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 jmc 1.12 ENDDO
104     ENDDO
105 adcroft 1.1
106 jmc 1.12 ELSEIF ( equationOfState(1:5).EQ.'JMD95'
107     & .OR. equationOfState.EQ.'UNESCO' ) THEN
108 mlosch 1.6 C nonlinear equation of state in pressure coordinates
109    
110 jmc 1.12 CALL PRESSURE_FOR_EOS(
111     I bi, bj, iMin, iMax, jMin, jMax, kRef,
112     O locPres,
113     I myThid )
114    
115 mlosch 1.6 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 jmc 1.12 I bi, bj, iMin, iMax, jMin, jMax, k,
123     I locPres, theta, salt,
124 mlosch 1.6 O bulkMod,
125     I myThid )
126    
127 jmc 1.12 DO j=jMin,jMax
128     DO i=iMin,iMax
129 mlosch 1.6
130     C abbreviations
131 mlosch 1.9 t1 = theta(i,j,k,bi,bj)
132     t2 = t1*t1
133     t3 = t2*t1
134 mlosch 1.6
135 mlosch 1.9 s1 = salt(i,j,k,bi,bj)
136 jmc 1.12 s3o2 = SQRT(s1*s1*s1)
137 mlosch 1.6
138 jmc 1.12 p1 = locPres(i,j)*SItoBar
139 mlosch 1.9 p2 = p1*p1
140 mlosch 1.6
141     C d(rho)/d(theta)
142     C of fresh water at p = 0
143     drhoP0dthetaFresh =
144     & eosJMDCFw(2)
145 mlosch 1.9 & + 2.*eosJMDCFw(3)*t1
146 mlosch 1.6 & + 3.*eosJMDCFw(4)*t2
147     & + 4.*eosJMDCFw(5)*t3
148 mlosch 1.9 & + 5.*eosJMDCFw(6)*t3*t1
149 mlosch 1.6 C of salt water at p = 0
150     drhoP0dthetaSalt =
151 mlosch 1.9 & s1*(
152 mlosch 1.6 & eosJMDCSw(2)
153 mlosch 1.9 & + 2.*eosJMDCSw(3)*t1
154 mlosch 1.6 & + 3.*eosJMDCSw(4)*t2
155     & + 4.*eosJMDCSw(5)*t3
156     & )
157     & + s3o2*(
158     & + eosJMDCSw(7)
159 mlosch 1.9 & + 2.*eosJMDCSw(8)*t1
160 mlosch 1.6 & )
161     C d(bulk modulus)/d(theta)
162     C of fresh water at p = 0
163     dKdthetaFresh =
164     & eosJMDCKFw(2)
165 mlosch 1.9 & + 2.*eosJMDCKFw(3)*t1
166 mlosch 1.6 & + 3.*eosJMDCKFw(4)*t2
167     & + 4.*eosJMDCKFw(5)*t3
168     C of sea water at p = 0
169     dKdthetaSalt =
170 mlosch 1.9 & s1*( eosJMDCKSw(2)
171     & + 2.*eosJMDCKSw(3)*t1
172 mlosch 1.6 & + 3.*eosJMDCKSw(4)*t2
173     & )
174     & + s3o2*( eosJMDCKSw(6)
175 mlosch 1.9 & + 2.*eosJMDCKSw(7)*t1
176 mlosch 1.6 & )
177     C of sea water at p
178     dKdthetaPres =
179 mlosch 1.9 & p1*( eosJMDCKP(2)
180     & + 2.*eosJMDCKP(3)*t1
181 mlosch 1.6 & + 3.*eosJMDCKP(4)*t2
182     & )
183 mlosch 1.9 & + p1*s1*( eosJMDCKP(6)
184     & + 2.*eosJMDCKP(7)*t1
185 mlosch 1.6 & )
186     & + p2*( eosJMDCKP(10)
187 mlosch 1.9 & + 2.*eosJMDCKP(11)*t1
188 mlosch 1.6 & )
189 mlosch 1.9 & + p2*s1*( eosJMDCKP(13)
190     & + 2.*eosJMDCKP(14)*t1
191 mlosch 1.6 & )
192    
193     drhoP0dtheta = drhoP0dthetaFresh
194     & + drhoP0dthetaSalt
195     dKdtheta = dKdthetaFresh
196     & + dKdthetaSalt
197     & + dKdthetaPres
198     alphaloc(i,j) =
199     & ( bulkmod(i,j)**2*drhoP0dtheta
200 mlosch 1.9 & - bulkmod(i,j)*p1*drhoP0dtheta
201     & - rhoP0(i,j)*p1*dKdtheta )
202     & /( bulkmod(i,j) - p1 )**2
203 mlosch 1.6
204    
205 jmc 1.12 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 heimbach 1.13 & kRef, theta, salt, rhoLoc, myThid )
216    
217 jmc 1.12 CALL FIND_RHODEN( bi, bj, iMin, iMax, jMin, jMax, k,
218     & locPres, theta, salt, rhoDen, myThid )
219 mlosch 1.9
220 jmc 1.12 DO j=jMin,jMax
221     DO i=iMin,iMax
222 mlosch 1.9 t1 = theta(i,j,k,bi,bj)
223     t2 = t1*t1
224     s1 = salt(i,j,k,bi,bj)
225 jmc 1.12 sp5 = SQRT(s1)
226 mlosch 1.9
227 jmc 1.12 p1 = locPres(i,j)*SItodBar
228 mlosch 1.9 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 jmc 1.12 ENDDO
248     ENDDO
249 mlosch 1.9
250 jmc 1.12 ELSE
251     WRITE(*,*) 'FIND_ALPHA: equationOfState = ',equationOfState
252     STOP 'FIND_ALPHA: "equationOfState" has illegal value'
253     ENDIF
254 adcroft 1.1
255 jmc 1.12 RETURN
256     END
257 adcroft 1.1
258 jmc 1.12 SUBROUTINE FIND_BETA (
259 mlosch 1.10 I bi, bj, iMin, iMax, jMin, jMax, k, kRef,
260 adcroft 1.1 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 jmc 1.12 IMPLICIT NONE
274 adcroft 1.1
275 jmc 1.12 C === Global variables ===
276 adcroft 1.1 #include "SIZE.h"
277     #include "DYNVARS.h"
278     #include "EEPARAMS.h"
279     #include "PARAMS.h"
280 mlosch 1.6 #include "EOS.h"
281     #include "GRID.h"
282 adcroft 1.1
283 jmc 1.12 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 adcroft 1.1 _RL betaloc(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
290    
291 jmc 1.12 C == Local variables ==
292     INTEGER i,j
293 adcroft 1.1 _RL refTemp,refSalt,tP,sP
294 mlosch 1.9 _RL t1, t2, t3, s1, s3o2, p1, sp5, p1t1
295 mlosch 1.6 _RL drhoP0dS
296     _RL dKdS, dKdSSalt, dKdSPres
297 jmc 1.12 _RL locPres(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
298     _RL rhoP0 (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
299 mlosch 1.6 _RL bulkMod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
300 mlosch 1.9 _RL dnum_dsalt, dden_dsalt
301 jmc 1.12 _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 mlosch 1.6 CEOP
305    
306 mlosch 1.11 #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 jmc 1.12 IF (equationOfState.EQ.'LINEAR') THEN
312 adcroft 1.1
313 jmc 1.12 DO j=jMin,jMax
314     DO i=iMin,iMax
315 adcroft 1.1 betaloc(i,j) = rhonil * sBeta
316 jmc 1.12 ENDDO
317     ENDDO
318 adcroft 1.1
319 jmc 1.12 ELSEIF (equationOfState.EQ.'POLY3') THEN
320 adcroft 1.1
321     refTemp=eosRefT(kRef)
322     refSalt=eosRefS(kRef)
323    
324 jmc 1.12 DO j=jMin,jMax
325     DO i=iMin,iMax
326 adcroft 1.1 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 jmc 1.12 ENDDO
344     ENDDO
345 adcroft 1.1
346 jmc 1.12 ELSEIF ( equationOfState(1:5).EQ.'JMD95'
347     & .OR. equationOfState.EQ.'UNESCO' ) THEN
348 mlosch 1.6 C nonlinear equation of state in pressure coordinates
349    
350 jmc 1.12 CALL PRESSURE_FOR_EOS(
351     I bi, bj, iMin, iMax, jMin, jMax, kRef,
352     O locPres,
353     I myThid )
354    
355 mlosch 1.6 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 jmc 1.12 I bi, bj, iMin, iMax, jMin, jMax, k,
363     I locPres, theta, salt,
364 mlosch 1.6 O bulkMod,
365     I myThid )
366    
367 jmc 1.12 DO j=jMin,jMax
368     DO i=iMin,iMax
369 mlosch 1.6
370     C abbreviations
371 mlosch 1.9 t1 = theta(i,j,k,bi,bj)
372     t2 = t1*t1
373     t3 = t2*t1
374 mlosch 1.6
375 mlosch 1.9 s1 = salt(i,j,k,bi,bj)
376 jmc 1.12 s3o2 = 1.5*SQRT(s1)
377 mlosch 1.6
378 jmc 1.12 p1 = locPres(i,j)*SItoBar
379 mlosch 1.6
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 mlosch 1.9 & + eosJMDCSw(2)*t1
387 mlosch 1.6 & + eosJMDCSw(3)*t2
388     & + eosJMDCSw(4)*t3
389 mlosch 1.9 & + eosJMDCSw(5)*t3*t1
390 mlosch 1.6 & + s3o2*(
391     & eosJMDCSw(6)
392 mlosch 1.9 & + eosJMDCSw(7)*t1
393 mlosch 1.6 & + eosJMDCSw(8)*t2
394     & )
395 mlosch 1.9 & + 2*eosJMDCSw(9)*s1
396 mlosch 1.6 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 mlosch 1.9 & + eosJMDCKSw(2)*t1
403 mlosch 1.6 & + eosJMDCKSw(3)*t2
404     & + eosJMDCKSw(4)*t3
405     & + s3o2*( eosJMDCKSw(5)
406 mlosch 1.9 & + eosJMDCKSw(6)*t1
407 mlosch 1.6 & + eosJMDCKSw(7)*t2
408     & )
409    
410     C of sea water at p
411     dKdSPres =
412 mlosch 1.9 & p1*( eosJMDCKP(5)
413     & + eosJMDCKP(6)*t1
414 mlosch 1.6 & + eosJMDCKP(7)*t2
415     & )
416 mlosch 1.9 & + s3o2*p1*eosJMDCKP(8)
417     & + p1*p1*( eosJMDCKP(12)
418     & + eosJMDCKP(13)*t1
419 mlosch 1.6 & + eosJMDCKP(14)*t2
420     & )
421    
422     dKdS = dKdSSalt + dKdSPres
423    
424     betaloc(i,j) =
425     & ( bulkmod(i,j)**2*drhoP0dS
426 mlosch 1.9 & - bulkmod(i,j)*p1*drhoP0dS
427     & - rhoP0(i,j)*p1*dKdS )
428     & /( bulkmod(i,j) - p1 )**2
429 mlosch 1.6
430    
431 jmc 1.12 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 heimbach 1.13 & kRef, theta, salt, rhoLoc, myThid )
442    
443 jmc 1.12 CALL FIND_RHODEN( bi, bj, iMin, iMax, jMin, jMax, k,
444     & locPres, theta, salt, rhoDen, myThid )
445 mlosch 1.9
446 jmc 1.12 DO j=jMin,jMax
447     DO i=iMin,iMax
448 mlosch 1.9 t1 = theta(i,j,k,bi,bj)
449     t2 = t1*t1
450     s1 = salt(i,j,k,bi,bj)
451 jmc 1.12 sp5 = SQRT(s1)
452 mlosch 1.9
453 jmc 1.12 p1 = locPres(i,j)*SItodBar
454 mlosch 1.9 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 jmc 1.12 ENDDO
467     ENDDO
468 mlosch 1.9
469 jmc 1.12 ELSE
470     WRITE(*,*) 'FIND_BETA: equationOfState = ',equationOfState
471     STOP 'FIND_BETA: "equationOfState" has illegal value'
472     ENDIF
473 adcroft 1.1
474 jmc 1.12 RETURN
475     END

  ViewVC Help
Powered by ViewVC 1.1.22