/[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.9 by mlosch, Thu Sep 5 20:49:33 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 t1, t2, t3, s1, s3o2, p1, p2, sp5, p1t1
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          _RL dnum_dtheta, dden_dtheta
58          _RL rhoDen(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
59          _RL rhoLoc(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
60          integer myThid
61  CEOP  CEOP
62    
63        if (eqn.eq.'LINEAR') then  Cml      stop 'myThid is not properly defined in this subroutine'
64    
65          if (equationOfState.eq.'LINEAR') then
66    
67           do j=jMin,jMax           do j=jMin,jMax
68              do i=iMin,iMax              do i=iMin,iMax
# Line 57  CEOP Line 70  CEOP
70              enddo              enddo
71           enddo           enddo
72                    
73        elseif (eqn.eq.'POLY3') then        elseif (equationOfState.eq.'POLY3') then
74    
75           refTemp=eosRefT(kRef)           refTemp=eosRefT(kRef)
76           refSalt=eosRefS(kRef)           refSalt=eosRefS(kRef)
# Line 86  CEOP Line 99  CEOP
99              enddo              enddo
100           enddo           enddo
101                    
102          elseif ( equationOfState(1:5).eq.'JMD95'
103         &        .or. equationOfState.eq.'UNESCO' ) then
104    C     nonlinear equation of state in pressure coordinates
105    
106             CALL FIND_RHOP0(
107         I        bi, bj, iMin, iMax, jMin, jMax, k,
108         I        theta, salt,
109         O        rhoP0,
110         I        myThid )
111            
112             CALL FIND_BULKMOD(
113         I        bi, bj, iMin, iMax, jMin, jMax,  k, kRef,
114         I        theta, salt,
115         O        bulkMod,
116         I        myThid )
117    
118             do j=jMin,jMax
119                do i=iMin,iMax
120    
121    C     abbreviations
122                   t1 = theta(i,j,k,bi,bj)
123                   t2 = t1*t1
124                   t3 = t2*t1
125                  
126                   s1  = salt(i,j,k,bi,bj)
127                   if ( s1 .lt. 0. _d 0 ) then
128    C     issue a warning
129                      write(*,'(a,i3,a,i3,a,i3,a,e13.5)')
130         &                 ' FIND_ALPHA: WARNING, salinity(',
131         &                 i,',',j,',',k,',',bi,',',bj,') = ', s1
132                      s1 = 0. _d 0
133                   end if
134                   s3o2 = sqrt(s1*s1*s1)
135                  
136                   p1  = pressure(i,j,kRef,bi,bj)*SItoBar
137                   p2 = p1*p1
138    
139    C     d(rho)/d(theta)
140    C     of fresh water at p = 0
141                   drhoP0dthetaFresh =
142         &                eosJMDCFw(2)
143         &           + 2.*eosJMDCFw(3)*t1
144         &           + 3.*eosJMDCFw(4)*t2
145         &           + 4.*eosJMDCFw(5)*t3
146         &           + 5.*eosJMDCFw(6)*t3*t1
147    C     of salt water at p = 0
148                   drhoP0dthetaSalt =
149         &        s1*(
150         &                eosJMDCSw(2)
151         &           + 2.*eosJMDCSw(3)*t1
152         &           + 3.*eosJMDCSw(4)*t2
153         &           + 4.*eosJMDCSw(5)*t3
154         &           )
155         &       + s3o2*(
156         &           +    eosJMDCSw(7)
157         &           + 2.*eosJMDCSw(8)*t1
158         &           )
159    C     d(bulk modulus)/d(theta)
160    C     of fresh water at p = 0
161                   dKdthetaFresh =
162         &                eosJMDCKFw(2)
163         &           + 2.*eosJMDCKFw(3)*t1
164         &           + 3.*eosJMDCKFw(4)*t2
165         &           + 4.*eosJMDCKFw(5)*t3
166    C     of sea water at p = 0
167                   dKdthetaSalt =
168         &        s1*(    eosJMDCKSw(2)
169         &           + 2.*eosJMDCKSw(3)*t1
170         &           + 3.*eosJMDCKSw(4)*t2
171         &           )
172         &    + s3o2*(    eosJMDCKSw(6)
173         &           + 2.*eosJMDCKSw(7)*t1
174         &           )
175    C     of sea water at p
176                   dKdthetaPres =
177         &        p1*(    eosJMDCKP(2)
178         &           + 2.*eosJMDCKP(3)*t1
179         &           + 3.*eosJMDCKP(4)*t2
180         &           )
181         &   + p1*s1*(    eosJMDCKP(6)
182         &           + 2.*eosJMDCKP(7)*t1
183         &           )
184         &      + p2*(    eosJMDCKP(10)
185         &           + 2.*eosJMDCKP(11)*t1
186         &           )
187         &   + p2*s1*(    eosJMDCKP(13)
188         &           + 2.*eosJMDCKP(14)*t1
189         &           )
190    
191                   drhoP0dtheta  = drhoP0dthetaFresh
192         &                       + drhoP0dthetaSalt
193                   dKdtheta      = dKdthetaFresh
194         &                       + dKdthetaSalt
195         &                       + dKdthetaPres
196                   alphaloc(i,j) =
197         &              ( bulkmod(i,j)**2*drhoP0dtheta
198         &              - bulkmod(i,j)*p1*drhoP0dtheta
199         &              - rhoP0(i,j)*p1*dKdtheta )
200         &              /( bulkmod(i,j) - p1 )**2
201                  
202                  
203                end do
204             end do
205          elseif ( equationOfState.eq.'MDJWF' ) then
206    
207             CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax,  k, kRef,
208         &        eqn, theta, salt, rhoLoc, myThid )
209             CALL FIND_RHODEN( bi, bj, iMin, iMax, jMin, jMax, k, kRef,
210         &      theta, salt, rhoDen, myThid )
211    
212             do j=jMin,jMax
213                do i=iMin,iMax
214                   t1  = theta(i,j,k,bi,bj)
215                   t2  = t1*t1
216                   s1  = salt(i,j,k,bi,bj)
217                   if ( s1 .lt. 0. _d 0 ) then
218    C     issue a warning
219                      write(*,'(a,i3,a,i3,a,i3,a,e13.5)')
220         &                 ' FIND_ALPHA: WARNING, salinity(',
221         &                 i,',',j,',',k,',',bi,',',bj,') = ', s1
222                      s1 = 0. _d 0
223                   end if
224                   sp5 = sqrt(s1)
225    
226                   p1   = pressure(i,j,kRef,bi,bj)*SItodBar
227                   p1t1 = p1*t1
228                
229                   dnum_dtheta = eosMDJWFnum(1)
230         &              + t1*(2.*eosMDJWFnum(2) + 3.*eosMDJWFnum(3)*t1)
231         &              + eosMDJWFnum(5)*s1                                
232         &              + p1t1*(2.*eosMDJWFnum(8) + 2.*eosMDJWFnum(11)*p1)    
233                
234                   dden_dtheta = eosMDJWFden(1)
235         &              + t1*(2.*eosMDJWFden(2)
236         &              +     t1*(3.*eosMDJWFden(3)
237         &              +         4.*eosMDJWFden(4)*t1 ) )
238         &              + s1*(eosMDJWFden(6)
239         &              +     t1*(3.*eosMDJWFden(7)*t1
240         &              +         2.*eosMDJWFden(9)*sp5 ) )
241         &              + p1*p1*(3.*eosMDJWFden(11)*t2 + eosMDJWFden(12)*p1)
242                
243                   alphaLoc(i,j)    = rhoDen(i,j)*(dnum_dtheta
244         &              - rhoLoc(i,j)*dden_dtheta)
245                
246             end do
247          end do
248    
249        else        else
250           write(*,*) 'FIND_ALPHA: eqn = ',eqn           write(*,*) 'FIND_ALPHA: equationOfState = ',equationOfState
251           stop 'FIND_ALPHA: argument "eqn" has illegal value'           stop 'FIND_ALPHA: "equationOfState" has illegal value'
252        endif        endif
253    
254        return        return
# Line 117  c Common Line 277  c Common
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 Arguments
284        integer bi,bj,iMin,iMax,jMin,jMax        integer bi,bj,iMin,iMax,jMin,jMax
# Line 128  c Arguments Line 290  c Arguments
290  c Local  c Local
291        integer i,j        integer i,j
292        _RL refTemp,refSalt,tP,sP        _RL refTemp,refSalt,tP,sP
293          _RL t1, t2, t3, s1, s3o2, p1, sp5, p1t1
294          _RL drhoP0dS
295          _RL dKdS, dKdSSalt, dKdSPres
296          _RL rhoP0(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
297          _RL bulkMod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
298          _RL dnum_dsalt, dden_dsalt
299          _RL rhoDen(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
300          _RL rhoLoc(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
301          integer myThid
302    CEOP
303    
304    Cml      stop 'myThid is not properly defined in this subroutine'
305    
306        if (eqn.eq.'LINEAR') then        if (equationOfState.eq.'LINEAR') then
307    
308           do j=jMin,jMax           do j=jMin,jMax
309              do i=iMin,iMax              do i=iMin,iMax
# Line 137  c Local Line 311  c Local
311              enddo              enddo
312           enddo           enddo
313                    
314        elseif (eqn.eq.'POLY3') then        elseif (equationOfState.eq.'POLY3') then
315    
316           refTemp=eosRefT(kRef)           refTemp=eosRefT(kRef)
317           refSalt=eosRefS(kRef)           refSalt=eosRefS(kRef)
# Line 164  c Local Line 338  c Local
338              enddo              enddo
339           enddo           enddo
340    
341          elseif ( equationOfState(1:5).eq.'JMD95'
342         &        .or. equationOfState.eq.'UNESCO' ) then
343    C     nonlinear equation of state in pressure coordinates
344    
345             CALL FIND_RHOP0(
346         I        bi, bj, iMin, iMax, jMin, jMax, k,
347         I        theta, salt,
348         O        rhoP0,
349         I        myThid )
350            
351             CALL FIND_BULKMOD(
352         I        bi, bj, iMin, iMax, jMin, jMax,  k, kRef,
353         I        theta, salt,
354         O        bulkMod,
355         I        myThid )
356    
357             do j=jMin,jMax
358                do i=iMin,iMax
359    
360    C     abbreviations
361                   t1 = theta(i,j,k,bi,bj)
362                   t2 = t1*t1
363                   t3 = t2*t1
364                  
365                   s1  = salt(i,j,k,bi,bj)
366                   if ( s1 .lt. 0. _d 0 ) then
367    C     issue a warning
368                      write(*,'(a,i3,a,i3,a,i3,a,e13.5)')
369         &                 ' FIND_BETA: WARNING, salinity(',
370         &                 i,',',j,',',k,',',bi,',',bj,') = ', s1
371                      s1 = 0. _d 0
372                   end if
373                   s3o2 = 1.5*sqrt(s1)
374    
375                   p1  = pressure(i,j,kRef,bi,bj)*SItoBar
376    
377    C     d(rho)/d(S)
378    C     of fresh water at p = 0
379                   drhoP0dS = 0. _d 0
380    C     of salt water at p = 0
381                   drhoP0dS = drhoP0dS
382         &              + eosJMDCSw(1)
383         &              + eosJMDCSw(2)*t1
384         &              + eosJMDCSw(3)*t2
385         &              + eosJMDCSw(4)*t3
386         &              + eosJMDCSw(5)*t3*t1
387         &       + s3o2*(
388         &                eosJMDCSw(6)
389         &              + eosJMDCSw(7)*t1
390         &              + eosJMDCSw(8)*t2
391         &              )
392         &              + 2*eosJMDCSw(9)*s1
393    C     d(bulk modulus)/d(S)
394    C     of fresh water at p = 0
395                   dKdS = 0. _d 0
396    C     of sea water at p = 0
397                   dKdSSalt =
398         &                eosJMDCKSw(1)
399         &              + eosJMDCKSw(2)*t1
400         &              + eosJMDCKSw(3)*t2
401         &              + eosJMDCKSw(4)*t3
402         &       + s3o2*( eosJMDCKSw(5)
403         &              + eosJMDCKSw(6)*t1
404         &              + eosJMDCKSw(7)*t2
405         &              )
406    
407    C     of sea water at p
408                   dKdSPres =
409         &           p1*( eosJMDCKP(5)
410         &              + eosJMDCKP(6)*t1
411         &              + eosJMDCKP(7)*t2
412         &              )
413         &        + s3o2*p1*eosJMDCKP(8)
414         &      + p1*p1*( eosJMDCKP(12)
415         &              + eosJMDCKP(13)*t1
416         &              + eosJMDCKP(14)*t2
417         &              )
418    
419                   dKdS = dKdSSalt + dKdSPres
420    
421                   betaloc(i,j) =
422         &              ( bulkmod(i,j)**2*drhoP0dS
423         &              - bulkmod(i,j)*p1*drhoP0dS
424         &              - rhoP0(i,j)*p1*dKdS )
425         &              /( bulkmod(i,j) - p1 )**2
426                  
427                  
428                end do
429             end do
430          elseif ( equationOfState.eq.'MDJWF' ) then
431    
432             CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax,  k, kRef,
433         &        eqn, theta, salt, rhoLoc, myThid )
434             CALL FIND_RHODEN( bi, bj, iMin, iMax, jMin, jMax, k, kRef,
435         &      theta, salt, rhoDen, myThid )
436            
437             do j=jMin,jMax
438                do i=iMin,iMax
439                   t1  = theta(i,j,k,bi,bj)
440                   t2  = t1*t1
441                   s1  = salt(i,j,k,bi,bj)
442                   if ( s1 .lt. 0. _d 0 ) then
443    C     issue a warning
444                      write(*,'(a,i3,a,i3,a,i3,a,e13.5)')
445         &                 ' FIND_BETA: WARNING, salinity(',
446         &                 i,',',j,',',k,',',bi,',',bj,') = ', s1
447                      s1 = 0. _d 0
448                   end if
449                   sp5 = sqrt(s1)
450            
451                   p1   = pressure(i,j,k,bi,bj)*SItodBar
452                   p1t1 = p1*t1
453                  
454                   dnum_dsalt = eosMDJWFnum(4)
455         &              + eosMDJWFnum(5)*t1
456         &              + 2.*eosMDJWFnum(6)*s1 + eosMDJWFnum(9)*p1
457                   dden_dsalt = eosMDJWFden(5)
458         &              + t1*( eosMDJWFden(6) + eosMDJWFden(7)*t2 )
459         &              + 1.5*sp5*(eosMDJWFden(8) + eosMDJWFden(9)*t2)
460    
461                   betaLoc(i,j) = rhoDen(i,j)*( dnum_dsalt
462         &              - rhoLoc(i,j)*dden_dsalt )
463    
464                end do
465             end do
466    
467        else        else
468           write(*,*) 'FIND_BETA: eqn = ',eqn           write(*,*) 'FIND_BETA: equationOfState = ',equationOfState
469           stop 'FIND_BETA: argument "eqn" has illegal value'           stop 'FIND_BETA: "equationOfState" has illegal value'
470        endif        endif
471    
472        return        return
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.9
  ViewVC Help
Powered by ViewVC 1.1.22