/[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.5 by cnh, Wed Sep 26 18:09:14 2001 UTC revision 1.6 by mlosch, Wed Aug 7 16:55:52 2002 UTC
# Line 34  c Common Line 34  c Common
34  #include "DYNVARS.h"  #include "DYNVARS.h"
35  #include "EEPARAMS.h"  #include "EEPARAMS.h"
36  #include "PARAMS.h"  #include "PARAMS.h"
37    #include "EOS.h"
38    #include "GRID.h"
39    
40  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
41  c Arguments  c Arguments
# Line 47  C     !LOCAL VARIABLES: Line 49  C     !LOCAL VARIABLES:
49  c Local  c Local
50        integer i,j        integer i,j
51        _RL refTemp,refSalt,tP,sP        _RL refTemp,refSalt,tP,sP
52          _RL t, t2, t3, s, s3o2, p, p2
53          _RL drhoP0dtheta, drhoP0dthetaFresh, drhoP0dthetaSalt
54          _RL dKdtheta, dKdthetaFresh, dKdthetaSalt, dKdthetaPres
55          _RL rhoP0(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
56          _RL bulkMod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
57          integer myThid
58  CEOP  CEOP
59    
60    Cml      stop 'myThid is not properly defined in this subroutine'
61    
62        if (eqn.eq.'LINEAR') then        if (eqn.eq.'LINEAR') then
63    
64           do j=jMin,jMax           do j=jMin,jMax
# Line 86  CEOP Line 96  CEOP
96              enddo              enddo
97           enddo           enddo
98                    
99          elseif ( eqn(1:5).eq.'JMD95' ) then
100    C     nonlinear equation of state in pressure coordinates
101    
102             CALL FIND_RHOP0(
103         I        bi, bj, iMin, iMax, jMin, jMax, k,
104         I        theta, salt,
105         O        rhoP0,
106         I        myThid )
107            
108             CALL FIND_BULKMOD(
109         I        bi, bj, iMin, iMax, jMin, jMax,  k,
110         I        theta, salt,
111         O        bulkMod,
112         I        myThid )
113    
114             do j=jMin,jMax
115                do i=iMin,iMax
116    
117    C     abbreviations
118                   t  = theta(i,j,k,bi,bj)
119                   t2 = t*t
120                   t3 = t2*t
121                  
122                   if ( salt(i,j,k,bi,bj) .lt. 0. _d 0 ) then
123    C     issue a warning
124                      write(*,'(a,i3,a,i3,a,i3,a,e13.5)')
125         &             ' FIND_BETA: WARNING, s(',i,',',j,',',k,') = ', s
126                      s = 0. _d 0
127                   else
128                      s  = salt(i,j,k,bi,bj)
129                   end if
130                   s3o2 = sqrt(s*s*s)
131                  
132                   p  = pressure(i,j,k,bi,bj)*SItoBar
133                   p2 = p*p
134    
135    C     d(rho)/d(theta)
136    C     of fresh water at p = 0
137                   drhoP0dthetaFresh =
138         &                eosJMDCFw(2)
139         &           + 2.*eosJMDCFw(3)*t
140         &           + 3.*eosJMDCFw(4)*t2
141         &           + 4.*eosJMDCFw(5)*t3
142         &           + 5.*eosJMDCFw(6)*t3*t
143    C     of salt water at p = 0
144                   drhoP0dthetaSalt =
145         &         s*(
146         &                eosJMDCSw(2)
147         &           + 2.*eosJMDCSw(3)*t
148         &           + 3.*eosJMDCSw(4)*t2
149         &           + 4.*eosJMDCSw(5)*t3
150         &           )
151         &       + s3o2*(
152         &           +    eosJMDCSw(7)
153         &           + 2.*eosJMDCSw(8)*t
154         &           )
155    C     d(bulk modulus)/d(theta)
156    C     of fresh water at p = 0
157                   dKdthetaFresh =
158         &                eosJMDCKFw(2)
159         &           + 2.*eosJMDCKFw(3)*t
160         &           + 3.*eosJMDCKFw(4)*t2
161         &           + 4.*eosJMDCKFw(5)*t3
162    C     of sea water at p = 0
163                   dKdthetaSalt =
164         &         s*(    eosJMDCKSw(2)
165         &           + 2.*eosJMDCKSw(3)*t
166         &           + 3.*eosJMDCKSw(4)*t2
167         &           )
168         &    + s3o2*(    eosJMDCKSw(6)
169         &           + 2.*eosJMDCKSw(7)*t
170         &           )
171    C     of sea water at p
172                   dKdthetaPres =
173         &         p*(    eosJMDCKP(2)
174         &           + 2.*eosJMDCKP(3)*t
175         &           + 3.*eosJMDCKP(4)*t2
176         &           )
177         &     + p*s*(    eosJMDCKP(6)
178         &           + 2.*eosJMDCKP(7)*t
179         &           )
180         &      + p2*(    eosJMDCKP(10)
181         &           + 2.*eosJMDCKP(11)*t
182         &           )
183         &    + p2*s*(    eosJMDCKP(13)
184         &           + 2.*eosJMDCKP(14)*t
185         &           )
186    
187                   drhoP0dtheta  = drhoP0dthetaFresh
188         &                       + drhoP0dthetaSalt
189                   dKdtheta      = dKdthetaFresh
190         &                       + dKdthetaSalt
191         &                       + dKdthetaPres
192                   alphaloc(i,j) =
193         &              ( bulkmod(i,j)**2*drhoP0dtheta
194         &              - bulkmod(i,j)*p*drhoP0dtheta
195         &              - rhoP0(i,j)*p*dKdtheta )
196         &              /( bulkmod(i,j) - p )**2
197                  
198                  
199                end do
200             end do
201        else        else
202           write(*,*) 'FIND_ALPHA: eqn = ',eqn           write(*,*) 'FIND_ALPHA: eqn = ',eqn
203           stop 'FIND_ALPHA: argument "eqn" has illegal value'           stop 'FIND_ALPHA: argument "eqn" has illegal value'
# Line 117  c Common Line 229  c Common
229  #include "DYNVARS.h"  #include "DYNVARS.h"
230  #include "EEPARAMS.h"  #include "EEPARAMS.h"
231  #include "PARAMS.h"  #include "PARAMS.h"
232    #include "EOS.h"
233    #include "GRID.h"
234    
235  c Arguments  c Arguments
236        integer bi,bj,iMin,iMax,jMin,jMax        integer bi,bj,iMin,iMax,jMin,jMax
# Line 128  c Arguments Line 242  c Arguments
242  c Local  c Local
243        integer i,j        integer i,j
244        _RL refTemp,refSalt,tP,sP        _RL refTemp,refSalt,tP,sP
245          _RL t, t2, t3, s, s3o2, p
246          _RL drhoP0dS
247          _RL dKdS, dKdSSalt, dKdSPres
248          _RL rhoP0(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
249          _RL bulkMod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
250          integer myThid
251    CEOP
252    
253    Cml      stop 'myThid is not properly defined in this subroutine'
254    
255        if (eqn.eq.'LINEAR') then        if (eqn.eq.'LINEAR') then
256    
# Line 164  c Local Line 287  c Local
287              enddo              enddo
288           enddo           enddo
289    
290          elseif ( eqn(1:5).eq.'JMD95' ) then
291    C     nonlinear equation of state in pressure coordinates
292    
293             CALL FIND_RHOP0(
294         I        bi, bj, iMin, iMax, jMin, jMax, k,
295         I        theta, salt,
296         O        rhoP0,
297         I        myThid )
298            
299             CALL FIND_BULKMOD(
300         I        bi, bj, iMin, iMax, jMin, jMax,  k,
301         I        theta, salt,
302         O        bulkMod,
303         I        myThid )
304    
305             do j=jMin,jMax
306                do i=iMin,iMax
307    
308    C     abbreviations
309                   t  = theta(i,j,k,bi,bj)
310                   t2 = t*t
311                   t3 = t2*t
312                  
313                   if ( salt(i,j,k,bi,bj) .lt. 0. _d 0 ) then
314    C     issue a warning
315                      write(*,'(a,i3,a,i3,a,i3,a,e13.5)')
316         &             ' FIND_BETA: WARNING, s(',i,',',j,',',k,') = ', s
317                      s = 0. _d 0
318                   else
319                      s  = salt(i,j,k,bi,bj)
320                   end if
321                   s3o2 = 1.5*sqrt(s)
322    
323                   p  = pressure(i,j,k,bi,bj)*SItoBar
324    
325    C     d(rho)/d(S)
326    C     of fresh water at p = 0
327                   drhoP0dS = 0. _d 0
328    C     of salt water at p = 0
329                   drhoP0dS = drhoP0dS
330         &              + eosJMDCSw(1)
331         &              + eosJMDCSw(2)*t
332         &              + eosJMDCSw(3)*t2
333         &              + eosJMDCSw(4)*t3
334         &              + eosJMDCSw(5)*t3*t
335         &       + s3o2*(
336         &                eosJMDCSw(6)
337         &              + eosJMDCSw(7)*t
338         &              + eosJMDCSw(8)*t2
339         &              )
340         &              + 2*eosJMDCSw(9)*s
341    C     d(bulk modulus)/d(S)
342    C     of fresh water at p = 0
343                   dKdS = 0. _d 0
344    C     of sea water at p = 0
345                   dKdSSalt =
346         &                eosJMDCKSw(1)
347         &              + eosJMDCKSw(2)*t
348         &              + eosJMDCKSw(3)*t2
349         &              + eosJMDCKSw(4)*t3
350         &       + s3o2*( eosJMDCKSw(5)
351         &              + eosJMDCKSw(6)*t
352         &              + eosJMDCKSw(7)*t2
353         &              )
354    
355    C     of sea water at p
356                   dKdSPres =
357         &            p*( eosJMDCKP(5)
358         &              + eosJMDCKP(6)*t
359         &              + eosJMDCKP(7)*t2
360         &              )
361         &        + s3o2*p*eosJMDCKP(8)
362         &        + p*p*( eosJMDCKP(12)
363         &              + eosJMDCKP(13)*t
364         &              + eosJMDCKP(14)*t2
365         &              )
366    
367                   dKdS = dKdSSalt + dKdSPres
368    
369                   betaloc(i,j) =
370         &              ( bulkmod(i,j)**2*drhoP0dS
371         &              - bulkmod(i,j)*p*drhoP0dS
372         &              - rhoP0(i,j)*p*dKdS )
373         &              /( bulkmod(i,j) - p )**2
374                  
375                  
376                end do
377             end do
378        else        else
379           write(*,*) 'FIND_BETA: eqn = ',eqn           write(*,*) 'FIND_BETA: eqn = ',eqn
380           stop 'FIND_BETA: argument "eqn" has illegal value'           stop 'FIND_BETA: argument "eqn" has illegal value'
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.6
  ViewVC Help
Powered by ViewVC 1.1.22