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

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

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

revision 1.1 by adcroft, Tue May 18 18:01:13 1999 UTC revision 1.13 by heimbach, Fri Mar 7 23:46:46 2003 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
5  #define USE_FACTORIZED_POLY  #define USE_FACTORIZED_POLY
6    
7        subroutine FIND_ALPHA (  CBOP
8       I     bi, bj, iMin, iMax, jMin, jMax,  k, kRef, eqn,  C     !ROUTINE: FIND_ALPHA
9    C     !INTERFACE:
10          SUBROUTINE FIND_ALPHA (
11         I     bi, bj, iMin, iMax, jMin, jMax,  k, kRef,
12       O     alphaloc )       O     alphaloc )
 C     /==========================================================\  
 C     | o SUBROUTINE FIND_ALPHA                                  |  
 C     |   Calculates [drho(S,T,z) / dT] of a horizontal slice    |  
 C     |==========================================================|  
 C     |                                                          |  
 C     | k - is the Theta/Salt level                              |  
 C     | kRef - determines pressure reference level               |  
 C     |        (not used in 'LINEAR' mode)                       |  
 C     | eqn - determines the eqn. of state: 'LINEAR' or 'POLY3'  |  
 C     |                                                          |  
 C     | alphaloc - drho / dT (kg/m^3/C)                          |  
 C     |                                                          |  
 C     \==========================================================/  
       implicit none  
13    
14  c Common  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"  #include "SIZE.h"
33  #include "DYNVARS.h"  #include "DYNVARS.h"
34  #include "EEPARAMS.h"  #include "EEPARAMS.h"
35  #include "PARAMS.h"  #include "PARAMS.h"
36    #include "EOS.h"
37    #include "GRID.h"
38    
39  c Arguments  C     !INPUT/OUTPUT PARAMETERS:
40        integer bi,bj,iMin,iMax,jMin,jMax  C     == Routine arguments ==
41        integer k                 ! Level of Theta/Salt slice  C     k    :: Level of Theta/Salt slice
42        integer kRef              ! Pressure reference level  C     kRef :: Pressure reference level
43        character*(*) eqn        INTEGER bi,bj,iMin,iMax,jMin,jMax
44          INTEGER k
45          INTEGER kRef
46        _RL alphaloc(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL alphaloc(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
47    
48  c Local  C     !LOCAL VARIABLES:
49        integer i,j  C     == Local variables ==
50          INTEGER i,j
51        _RL refTemp,refSalt,tP,sP        _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 (eqn.eq.'LINEAR') then        IF (equationOfState.EQ.'LINEAR') THEN
70    
71           do j=jMin,jMax           DO j=jMin,jMax
72              do i=iMin,iMax              DO i=iMin,iMax
73                 alphaloc(i,j) = -rhonil * tAlpha                 alphaloc(i,j) = -rhonil * tAlpha
74              enddo              ENDDO
75           enddo           ENDDO
76                    
77        elseif (eqn.eq.'POLY3') then        ELSEIF (equationOfState.EQ.'POLY3') THEN
78    
79           refTemp=eosRefT(kRef)           refTemp=eosRefT(kRef)
80           refSalt=eosRefS(kRef)           refSalt=eosRefS(kRef)
81    
82           do j=jMin,jMax           DO j=jMin,jMax
83              do i=iMin,iMax              DO i=iMin,iMax
84                 tP=theta(i,j,k,bi,bj)-refTemp                 tP=theta(i,j,k,bi,bj)-refTemp
85                 sP=salt(i,j,k,bi,bj)-refSalt                 sP=salt(i,j,k,bi,bj)-refSalt
86  #ifdef USE_FACTORIZED_POLY  #ifdef USE_FACTORIZED_POLY
# Line 72  c Local Line 100  c Local
100       &              eosC(7,kRef)*tP*2.   *sP    +       &              eosC(7,kRef)*tP*2.   *sP    +
101       &              eosC(8,kRef)         *sP*sP       &              eosC(8,kRef)         *sP*sP
102  #endif  #endif
103              enddo              ENDDO
104           enddo           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        else           CALL FIND_BULKMOD(
122           write(0,*) 'FIND_ALPHA: eqn = ',eqn       I        bi, bj, iMin, iMax, jMin, jMax,  k,
123           stop 'FIND_ALPHA: argument "eqn" has illegal value'       I        locPres, theta, salt,
124        endif       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        RETURN
256        end        END
257    
258        subroutine FIND_BETA (        SUBROUTINE FIND_BETA (
259       I     bi, bj, iMin, iMax, jMin, jMax,  k, kRef, eqn,       I     bi, bj, iMin, iMax, jMin, jMax,  k, kRef,
260       O     betaloc )       O     betaloc )
261  C     /==========================================================\  C     /==========================================================\
262  C     | o SUBROUTINE FIND_BETA                                   |  C     | o SUBROUTINE FIND_BETA                                   |
# Line 94  C     | Line 266  C     |
266  C     | k - is the Theta/Salt level                              |  C     | k - is the Theta/Salt level                              |
267  C     | kRef - determines pressure reference level               |  C     | kRef - determines pressure reference level               |
268  C     |        (not used in 'LINEAR' mode)                       |  C     |        (not used in 'LINEAR' mode)                       |
 C     | eqn - determines the eqn. of state: 'LINEAR' or 'POLY3'  |  
269  C     |                                                          |  C     |                                                          |
270  C     | betaloc - drho / dS (kg/m^3/PSU)                         |  C     | betaloc - drho / dS (kg/m^3/PSU)                         |
271  C     |                                                          |  C     |                                                          |
272  C     \==========================================================/  C     \==========================================================/
273        implicit none        IMPLICIT NONE
274    
275  c Common  C     === Global variables ===
276  #include "SIZE.h"  #include "SIZE.h"
277  #include "DYNVARS.h"  #include "DYNVARS.h"
278  #include "EEPARAMS.h"  #include "EEPARAMS.h"
279  #include "PARAMS.h"  #include "PARAMS.h"
280    #include "EOS.h"
281    #include "GRID.h"
282    
283  c Arguments  C     == Routine arguments ==
284        integer bi,bj,iMin,iMax,jMin,jMax  C     k    :: Level of Theta/Salt slice
285        integer k                 ! Level of Theta/Salt slice  C     kRef :: Pressure reference level
286        integer kRef              ! Pressure reference level        INTEGER bi,bj,iMin,iMax,jMin,jMax
287        character*(*) eqn        INTEGER k
288          INTEGER kRef
289        _RL betaloc(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL betaloc(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
290    
291  c Local  C     == Local variables ==
292        integer i,j        INTEGER i,j
293        _RL refTemp,refSalt,tP,sP        _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 (eqn.eq.'LINEAR') then        IF (equationOfState.EQ.'LINEAR') THEN
312    
313           do j=jMin,jMax           DO j=jMin,jMax
314              do i=iMin,iMax              DO i=iMin,iMax
315                 betaloc(i,j) = rhonil * sBeta                 betaloc(i,j) = rhonil * sBeta
316              enddo              ENDDO
317           enddo           ENDDO
318                    
319        elseif (eqn.eq.'POLY3') then        ELSEIF (equationOfState.EQ.'POLY3') THEN
320    
321           refTemp=eosRefT(kRef)           refTemp=eosRefT(kRef)
322           refSalt=eosRefS(kRef)           refSalt=eosRefS(kRef)
323    
324           do j=jMin,jMax           DO j=jMin,jMax
325              do i=iMin,iMax              DO i=iMin,iMax
326                 tP=theta(i,j,k,bi,bj)-refTemp                 tP=theta(i,j,k,bi,bj)-refTemp
327                 sP=salt(i,j,k,bi,bj)-refSalt                 sP=salt(i,j,k,bi,bj)-refSalt
328  #ifdef USE_FACTORIZED_POLY  #ifdef USE_FACTORIZED_POLY
# Line 150  c Local Line 340  c Local
340       &              eosC(8,kRef)*tP      *sP*2. +       &              eosC(8,kRef)*tP      *sP*2. +
341       &              eosC(9,kRef)         *sP*sP*3.       &              eosC(9,kRef)         *sP*sP*3.
342  #endif  #endif
343              enddo              ENDDO
344           enddo           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        else           CALL FIND_RHODEN( bi, bj, iMin, iMax, jMin, jMax, k,
444           write(0,*) 'FIND_BETA: eqn = ',eqn       &        locPres, theta, salt, rhoDen, myThid )
445           stop 'FIND_BETA: argument "eqn" has illegal value'          
446        endif           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        RETURN
475        end        END
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.13
  ViewVC Help
Powered by ViewVC 1.1.22