/[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.11 by mlosch, Fri Nov 1 22:00:33 2002 UTC revision 1.13 by heimbach, Fri Mar 7 23:46:46 2003 UTC
# Line 28  C     \ev Line 28  C     \ev
28    
29  C     !USES:  C     !USES:
30        IMPLICIT NONE        IMPLICIT NONE
31  c Common  C     === Global variables ===
32  #include "SIZE.h"  #include "SIZE.h"
33  #include "DYNVARS.h"  #include "DYNVARS.h"
34  #include "EEPARAMS.h"  #include "EEPARAMS.h"
# Line 37  c Common Line 37  c Common
37  #include "GRID.h"  #include "GRID.h"
38    
39  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
40  c Arguments  C     == Routine arguments ==
41        integer bi,bj,iMin,iMax,jMin,jMax  C     k    :: Level of Theta/Salt slice
42        integer k                 ! Level of Theta/Salt slice  C     kRef :: Pressure reference level
43        integer kRef              ! Pressure reference level        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 VARIABLES:  C     !LOCAL VARIABLES:
49  c Local  C     == Local variables ==
50        integer i,j        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        _RL t1, t2, t3, s1, s3o2, p1, p2, sp5, p1t1
53        _RL drhoP0dtheta, drhoP0dthetaFresh, drhoP0dthetaSalt        _RL drhoP0dtheta, drhoP0dthetaFresh, drhoP0dthetaSalt
54        _RL dKdtheta, dKdthetaFresh, dKdthetaSalt, dKdthetaPres        _RL dKdtheta, dKdthetaFresh, dKdthetaSalt, dKdthetaPres
55        _RL rhoP0(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _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)        _RL bulkMod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
58        _RL dnum_dtheta, dden_dtheta        _RL dnum_dtheta, dden_dtheta
59        _RL rhoDen(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL rhoDen (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
60        _RL rhoLoc(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL rhoLoc (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
61        integer myThid        INTEGER myThid
62  CEOP  CEOP
63    
 Cml      stop 'myThid is not properly defined in this subroutine'  
   
64  #ifdef CHECK_SALINITY_FOR_NEGATIVE_VALUES  #ifdef CHECK_SALINITY_FOR_NEGATIVE_VALUES
65        CALL LOOK_FOR_NEG_SALINITY( bi, bj, iMin, iMax, jMin, jMax,  k,        CALL LOOK_FOR_NEG_SALINITY( bi, bj, iMin, iMax, jMin, jMax,  k,
66       &     sFld, myThid )       &     sFld, myThid )
67  #endif  #endif
68    
69        if (equationOfState.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 (equationOfState.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 99  Cml      stop 'myThid is not properly de Line 100  Cml      stop 'myThid is not properly de
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'        ELSEIF ( equationOfState(1:5).EQ.'JMD95'
107       &        .or. equationOfState.eq.'UNESCO' ) then       &        .OR. equationOfState.EQ.'UNESCO' ) THEN
108  C     nonlinear equation of state in pressure coordinates  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(           CALL FIND_RHOP0(
116       I        bi, bj, iMin, iMax, jMin, jMax, k,       I        bi, bj, iMin, iMax, jMin, jMax, k,
117       I        theta, salt,       I        theta, salt,
# Line 113  C     nonlinear equation of state in pre Line 119  C     nonlinear equation of state in pre
119       I        myThid )       I        myThid )
120                    
121           CALL FIND_BULKMOD(           CALL FIND_BULKMOD(
122       I        bi, bj, iMin, iMax, jMin, jMax,  k, kRef,       I        bi, bj, iMin, iMax, jMin, jMax,  k,
123       I        theta, salt,       I        locPres, theta, salt,
124       O        bulkMod,       O        bulkMod,
125       I        myThid )       I        myThid )
126    
127           do j=jMin,jMax           DO j=jMin,jMax
128              do i=iMin,iMax              DO i=iMin,iMax
129    
130  C     abbreviations  C     abbreviations
131                 t1 = theta(i,j,k,bi,bj)                 t1 = theta(i,j,k,bi,bj)
# Line 127  C     abbreviations Line 133  C     abbreviations
133                 t3 = t2*t1                 t3 = t2*t1
134                                
135                 s1  = salt(i,j,k,bi,bj)                 s1  = salt(i,j,k,bi,bj)
136                 s3o2 = sqrt(s1*s1*s1)                 s3o2 = SQRT(s1*s1*s1)
137                                
138                 p1  = pressure(i,j,kRef,bi,bj)*SItoBar                 p1  = locPres(i,j)*SItoBar
139                 p2 = p1*p1                 p2 = p1*p1
140    
141  C     d(rho)/d(theta)  C     d(rho)/d(theta)
# Line 196  C     of sea water at p Line 202  C     of sea water at p
202       &              /( bulkmod(i,j) - p1 )**2       &              /( bulkmod(i,j) - p1 )**2
203                                
204                                
205              end do              ENDDO
206           end do           ENDDO
207        elseif ( equationOfState.eq.'MDJWF' ) then        ELSEIF ( equationOfState.EQ.'MDJWF' ) THEN
208    
209           CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax,  k, kRef,           CALL PRESSURE_FOR_EOS(
210       &        theta, salt, rhoLoc, myThid )       I        bi, bj, iMin, iMax, jMin, jMax,  kRef,
211           CALL FIND_RHODEN( bi, bj, iMin, iMax, jMin, jMax, k, kRef,       O        locPres,
212       &      theta, salt, rhoDen, myThid )       I        myThid )
213    
214           do j=jMin,jMax           CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax,  k,
215              do i=iMin,iMax       &        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)                 t1  = theta(i,j,k,bi,bj)
223                 t2  = t1*t1                 t2  = t1*t1
224                 s1  = salt(i,j,k,bi,bj)                 s1  = salt(i,j,k,bi,bj)
225                 sp5 = sqrt(s1)                 sp5 = SQRT(s1)
226    
227                 p1   = pressure(i,j,kRef,bi,bj)*SItodBar                 p1   = locPres(i,j)*SItodBar
228                 p1t1 = p1*t1                 p1t1 = p1*t1
229                            
230                 dnum_dtheta = eosMDJWFnum(1)                 dnum_dtheta = eosMDJWFnum(1)
# Line 232  C     of sea water at p Line 244  C     of sea water at p
244                 alphaLoc(i,j)    = rhoDen(i,j)*(dnum_dtheta                 alphaLoc(i,j)    = rhoDen(i,j)*(dnum_dtheta
245       &              - rhoLoc(i,j)*dden_dtheta)       &              - rhoLoc(i,j)*dden_dtheta)
246                            
247           end do           ENDDO
248        end do        ENDDO
249    
250        else        ELSE
251           write(*,*) 'FIND_ALPHA: equationOfState = ',equationOfState           WRITE(*,*) 'FIND_ALPHA: equationOfState = ',equationOfState
252           stop 'FIND_ALPHA: "equationOfState" has illegal value'           STOP 'FIND_ALPHA: "equationOfState" has illegal value'
253        endif        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,       I     bi, bj, iMin, iMax, jMin, jMax,  k, kRef,
260       O     betaloc )       O     betaloc )
261  C     /==========================================================\  C     /==========================================================\
# Line 258  C     | Line 270  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"
# Line 268  c Common Line 280  c Common
280  #include "EOS.h"  #include "EOS.h"
281  #include "GRID.h"  #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          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        _RL t1, t2, t3, s1, s3o2, p1, sp5, p1t1
295        _RL drhoP0dS        _RL drhoP0dS
296        _RL dKdS, dKdSSalt, dKdSPres        _RL dKdS, dKdSSalt, dKdSPres
297        _RL rhoP0(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _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)        _RL bulkMod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
300        _RL dnum_dsalt, dden_dsalt        _RL dnum_dsalt, dden_dsalt
301        _RL rhoDen(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL rhoDen (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
302        _RL rhoLoc(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL rhoLoc (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
303        integer myThid        INTEGER myThid
304  CEOP  CEOP
305    
 Cml      stop 'myThid is not properly defined in this subroutine'  
   
306  #ifdef CHECK_SALINITY_FOR_NEGATIVE_VALUES  #ifdef CHECK_SALINITY_FOR_NEGATIVE_VALUES
307        CALL LOOK_FOR_NEG_SALINITY( bi, bj, iMin, iMax, jMin, jMax,  k,        CALL LOOK_FOR_NEG_SALINITY( bi, bj, iMin, iMax, jMin, jMax,  k,
308       &     sFld, myThid )       &     sFld, myThid )
309  #endif  #endif
310    
311        if (equationOfState.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 (equationOfState.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 327  Cml      stop 'myThid is not properly de Line 340  Cml      stop 'myThid is not properly de
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'        ELSEIF ( equationOfState(1:5).EQ.'JMD95'
347       &        .or. equationOfState.eq.'UNESCO' ) then       &        .OR. equationOfState.EQ.'UNESCO' ) THEN
348  C     nonlinear equation of state in pressure coordinates  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(           CALL FIND_RHOP0(
356       I        bi, bj, iMin, iMax, jMin, jMax, k,       I        bi, bj, iMin, iMax, jMin, jMax, k,
357       I        theta, salt,       I        theta, salt,
# Line 341  C     nonlinear equation of state in pre Line 359  C     nonlinear equation of state in pre
359       I        myThid )       I        myThid )
360                    
361           CALL FIND_BULKMOD(           CALL FIND_BULKMOD(
362       I        bi, bj, iMin, iMax, jMin, jMax,  k, kRef,       I        bi, bj, iMin, iMax, jMin, jMax,  k,
363       I        theta, salt,       I        locPres, theta, salt,
364       O        bulkMod,       O        bulkMod,
365       I        myThid )       I        myThid )
366    
367           do j=jMin,jMax           DO j=jMin,jMax
368              do i=iMin,iMax              DO i=iMin,iMax
369    
370  C     abbreviations  C     abbreviations
371                 t1 = theta(i,j,k,bi,bj)                 t1 = theta(i,j,k,bi,bj)
# Line 355  C     abbreviations Line 373  C     abbreviations
373                 t3 = t2*t1                 t3 = t2*t1
374                                
375                 s1  = salt(i,j,k,bi,bj)                 s1  = salt(i,j,k,bi,bj)
376                 s3o2 = 1.5*sqrt(s1)                 s3o2 = 1.5*SQRT(s1)
377    
378                 p1  = pressure(i,j,kRef,bi,bj)*SItoBar                 p1  = locPres(i,j)*SItoBar
379    
380  C     d(rho)/d(S)  C     d(rho)/d(S)
381  C     of fresh water at p = 0  C     of fresh water at p = 0
# Line 410  C     of sea water at p Line 428  C     of sea water at p
428       &              /( bulkmod(i,j) - p1 )**2       &              /( bulkmod(i,j) - p1 )**2
429                                
430                                
431              end do              ENDDO
432           end do           ENDDO
433        elseif ( equationOfState.eq.'MDJWF' ) then        ELSEIF ( equationOfState.EQ.'MDJWF' ) THEN
434    
435           CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax,  k, kRef,           CALL PRESSURE_FOR_EOS(
436       &        theta, salt, rhoLoc, myThid )       I        bi, bj, iMin, iMax, jMin, jMax,  kRef,
437           CALL FIND_RHODEN( bi, bj, iMin, iMax, jMin, jMax, k, kRef,       O        locPres,
438       &      theta, salt, rhoDen, myThid )       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           DO j=jMin,jMax
447              do i=iMin,iMax              DO i=iMin,iMax
448                 t1  = theta(i,j,k,bi,bj)                 t1  = theta(i,j,k,bi,bj)
449                 t2  = t1*t1                 t2  = t1*t1
450                 s1  = salt(i,j,k,bi,bj)                 s1  = salt(i,j,k,bi,bj)
451                 sp5 = sqrt(s1)                 sp5 = SQRT(s1)
452                    
453                 p1   = pressure(i,j,k,bi,bj)*SItodBar                 p1   = locPres(i,j)*SItodBar
454                 p1t1 = p1*t1                 p1t1 = p1*t1
455                                
456                 dnum_dsalt = eosMDJWFnum(4)                 dnum_dsalt = eosMDJWFnum(4)
# Line 439  C     of sea water at p Line 463  C     of sea water at p
463                 betaLoc(i,j) = rhoDen(i,j)*( dnum_dsalt                 betaLoc(i,j) = rhoDen(i,j)*( dnum_dsalt
464       &              - rhoLoc(i,j)*dden_dsalt )       &              - rhoLoc(i,j)*dden_dsalt )
465    
466              end do              ENDDO
467           end do           ENDDO
468    
469        else        ELSE
470           write(*,*) 'FIND_BETA: equationOfState = ',equationOfState           WRITE(*,*) 'FIND_BETA: equationOfState = ',equationOfState
471           stop 'FIND_BETA: "equationOfState" has illegal value'           STOP 'FIND_BETA: "equationOfState" has illegal value'
472        endif        ENDIF
473    
474        return        RETURN
475        end        END
Legend:
Removed from v.1.11  
changed lines
  Added in v.1.13
  ViewVC Help
Powered by ViewVC 1.1.22