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

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

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

revision 1.8 by cnh, Fri Jun 12 19:33:33 1998 UTC revision 1.23 by heimbach, Fri Nov 15 03:01:21 2002 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "CPP_EEOPTIONS.h"  #include "CPP_OPTIONS.h"
5  #define USE_FACTORIZED_POLY  #define USE_FACTORIZED_POLY
6    
7  ! ==============================================================================  CBOP
8    C     !ROUTINE: FIND_RHO
9    C     !INTERFACE:
10        subroutine FIND_RHO(        subroutine FIND_RHO(
11       I      bi, bj, iMin, iMax, jMin, jMax,  k, kRef, eqn,       I      bi, bj, iMin, iMax, jMin, jMax,  k, kRef,
12         I      tFld, sFld,
13       O      rholoc,       O      rholoc,
14       I      myThid )       I      myThid )
15  C     /==========================================================\  
16  C     | o SUBROUTINE FIND_RHO                                    |  C     !DESCRIPTION: \bv
17  C     |   Calculates [rho(S,T,z)-Rhonil] of a slice              |  C     *==========================================================*
18  C     |==========================================================|  C     | o SUBROUTINE FIND_RHO                                    
19  C     |                                                          |  C     |   Calculates [rho(S,T,z)-rhoConst] of a slice              
20  C     | k - is the Theta/Salt level                              |  C     *==========================================================*
21  C     | kRef - determines pressure reference level               |  C     |                                                          
22  C     |        (not used in 'LINEAR' mode)                       |  C     | k - is the Theta/Salt level                              
23  C     | eqn - determines the eqn. of state: 'LINEAR' or 'POLY3'  |  C     | kRef - determines pressure reference level                
24  C     |                                                          |  C     |        (not used in 'LINEAR' mode)                        
25  C     \==========================================================/  C     |                                                          
26    C     *==========================================================*
27    C     \ev
28    
29    C     !USES:
30        implicit none        implicit none
31  ! Common  C     == Global variables ==
32  #include "SIZE.h"  #include "SIZE.h"
 #include "DYNVARS.h"  
33  #include "EEPARAMS.h"  #include "EEPARAMS.h"
34  #include "PARAMS.h"  #include "PARAMS.h"
35  ! Arguments  #include "EOS.h"
36    #include "GRID.h"
37    
38    C     !INPUT/OUTPUT PARAMETERS:
39    C     == Routine arguments ==
40        integer bi,bj,iMin,iMax,jMin,jMax        integer bi,bj,iMin,iMax,jMin,jMax
41        integer k                 ! Level of Theta/Salt slice        integer k                 ! Level of Theta/Salt slice
42        integer kRef              ! Pressure reference level        integer kRef              ! Pressure reference level
43        character*(*) eqn        _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
44          _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
45        _RL rholoc(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL rholoc(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
46        integer myThid        integer myThid
47  ! Local  
48    C     !LOCAL VARIABLES:
49    C     == Local variables ==
50        integer i,j        integer i,j
51        _RL refTemp,refSalt,sigRef,tP,sP,deltaSig        _RL refTemp,refSalt,sigRef,tP,sP,deltaSig,dRho
52  ! ------------------------------------------------------------------------------        _RL rhoP0(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
53          _RL bulkMod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
54          _RL rhoNum(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
55          _RL rhoDen(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
56          character*(max_len_mbuf) msgbuf
57    CEOP
58    
59        if (eqn.eq.'LINEAR') then  #ifdef ALLOW_AUTODIFF_TAMC
60          DO j=1-OLy,sNy+OLy
61           DO i=1-OLx,sNx+OLx
62            rholoc(i,j)  = 0. _d 0
63            rhoP0(i,j)   = 0. _d 0
64            bulkMod(i,j) = 0. _d 0
65           ENDDO
66          ENDDO
67    #endif
68    
69    #ifdef CHECK_SALINITY_FOR_NEGATIVE_VALUES
70          CALL LOOK_FOR_NEG_SALINITY( bi, bj, iMin, iMax, jMin, jMax,  k,
71         &     sFld, myThid )
72    #endif
73    
74          if (equationOfState.eq.'LINEAR') then
75    
76  C ***NOTE***  C ***NOTE***
77  C In the linear EOS, to make the static stability calculation meaningful  C In the linear EOS, to make the static stability calculation meaningful
# Line 46  C ********** Line 80  C **********
80         refTemp=tRef(kRef)         refTemp=tRef(kRef)
81         refSalt=sRef(kRef)         refSalt=sRef(kRef)
82    
83           dRho = rhoNil-rhoConst
84    
85         do j=jMin,jMax         do j=jMin,jMax
86          do i=iMin,iMax          do i=iMin,iMax
87           rholoc(i,j)=rhonil*(           rholoc(i,j)=rhoNil*(
88       &     sBeta*( salt(i,j,k,bi,bj)-refSalt)       &     sBeta*(sFld(i,j,k,bi,bj)-refSalt)
89       &   -tAlpha*(theta(i,j,k,bi,bj)-refTemp) )       &   -tAlpha*(tFld(i,j,k,bi,bj)-refTemp) )
90         &        + dRho
91          enddo          enddo
92         enddo         enddo
93                
94        elseif (eqn.eq.'POLY3') then        elseif (equationOfState.eq.'POLY3') then
95    
96         refTemp=eosRefT(kRef)         refTemp=eosRefT(kRef)
97         refSalt=eosRefS(kRef)         refSalt=eosRefS(kRef)
98         sigRef=eosSig0(kRef) + (1000.-Rhonil)         sigRef=eosSig0(kRef) + (1000.-rhoConst)
99    
100         do j=jMin,jMax         do j=jMin,jMax
101          do i=iMin,iMax          do i=iMin,iMax
102           tP=theta(i,j,k,bi,bj)-refTemp           tP=tFld(i,j,k,bi,bj)-refTemp
103           sP=salt(i,j,k,bi,bj)-refSalt           sP=sFld(i,j,k,bi,bj)-refSalt
104  #ifdef USE_FACTORIZED_POLY  #ifdef USE_FACTORIZED_POLY
105           deltaSig=           deltaSig=
106       &    (( eosC(9,kRef)*sP + eosC(5,kRef) )*sP + eosC(2,kRef) )*sP       &    (( eosC(9,kRef)*sP + eosC(5,kRef) )*sP + eosC(2,kRef) )*sP
# Line 89  C ********** Line 126  C **********
126          enddo          enddo
127         enddo         enddo
128    
129          elseif ( equationOfState(1:5).eq.'JMD95'
130         &      .or. equationOfState.eq.'UNESCO' ) then
131    C     nonlinear equation of state in pressure coordinates
132    
133             CALL FIND_RHOP0(
134         I        bi, bj, iMin, iMax, jMin, jMax, k,
135         I        tFld, sFld,
136         O        rhoP0,
137         I        myThid )
138            
139             CALL FIND_BULKMOD(
140         I        bi, bj, iMin, iMax, jMin, jMax,  k, kRef,
141         I        tFld, sFld,
142         O        bulkMod,
143         I        myThid )
144            
145    #ifdef ALLOW_AUTODIFF_TAMC
146    cph can't do storing here since find_rho is called multiple times;
147    cph additional recomp. should be acceptable
148    cphCADJ STORE rhoP0(:,:)   = comlev1_bibj_k ,  key=kkey , byte=isbyte
149    cphCADJ STORE bulkMod(:,:) = comlev1_bibj_k ,  key=kkey , byte=isbyte
150    #endif /* ALLOW_AUTODIFF_TAMC */
151             do j=jMin,jMax
152                do i=iMin,iMax
153    
154    C     density of sea water at pressure p
155                   rholoc(i,j) = rhoP0(i,j)
156         &              /(1. _d 0 -
157         &              pressure(i,j,kRef,bi,bj)*SItoBar/bulkMod(i,j))
158         &              - rhoConst
159    
160                end do
161             end do
162    
163          elseif ( equationOfState.eq.'MDJWF' ) then
164    
165             CALL FIND_RHONUM( bi, bj, iMin, iMax, jMin, jMax, k, kRef,
166         &      tFld, sFld, rhoNum, myThid )
167    
168             CALL FIND_RHODEN( bi, bj, iMin, iMax, jMin, jMax, k, kRef,
169         &      tFld, sFld, rhoDen, myThid )
170    #ifdef ALLOW_AUTODIFF_TAMC
171    cph can't do storing here since find_rho is called multiple times;
172    cph additional recomp. should be acceptable
173    cphCADJ STORE rhoNum(:,:) = comlev1_bibj_k ,  key=kkey , byte=isbyte
174    cphCADJ STORE rhoDen(:,:) = comlev1_bibj_k ,  key=kkey , byte=isbyte
175    #endif /* ALLOW_AUTODIFF_TAMC */
176             do j=jMin,jMax
177                do i=iMin,iMax
178    
179                   rholoc(i,j) = rhoNum(i,j)*rhoDen(i,j) - rhoConst
180    
181                end do
182             end do
183    
184          elseif( equationOfState .eq. 'IDEALG' ) then
185    C    
186          else
187           write(msgbuf,'(3a)')
188         &        ' FIND_RHO: equationOfState = "',equationOfState,'"'
189           call print_error( msgbuf, mythid )
190           stop 'ABNORMAL END: S/R FIND_RHO'
191          endif
192    
193          return
194          end
195    
196    CBOP
197    C     !ROUTINE: FIND_RHOP0
198    C     !INTERFACE:
199          subroutine FIND_RHOP0(
200         I      bi, bj, iMin, iMax, jMin, jMax, k,
201         I      tFld, sFld,
202         O      rhoP0,
203         I      myThid )
204    
205    C     !DESCRIPTION: \bv
206    C     *==========================================================*
207    C     | o SUBROUTINE FIND_RHOP0
208    C     |   Calculates rho(S,T,0) of a slice              
209    C     *==========================================================*
210    C     |                                                          
211    C     | k - is the surface level                              
212    C     |                                                          
213    C     *==========================================================*
214    C     \ev
215    
216    C     !USES:
217          implicit none
218    C     == Global variables ==
219    #include "SIZE.h"
220    #include "EEPARAMS.h"
221    #include "PARAMS.h"
222    #include "EOS.h"
223    
224    C     !INPUT/OUTPUT PARAMETERS:
225    C     == Routine arguments ==
226          integer bi,bj,iMin,iMax,jMin,jMax
227          integer k         ! Level of surface slice
228          _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
229          _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
230          _RL rhoP0(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
231          _RL rfresh, rsalt
232          integer myThid
233    
234    C     !LOCAL VARIABLES:
235    C     == Local variables ==
236          integer i,j
237          _RL t, t2, t3, t4, s, s3o2
238    CEOP
239    
240          do j=jMin,jMax
241             do i=iMin,iMax
242    C     abbreviations
243                t  = tFld(i,j,k,bi,bj)
244                t2 = t*t
245                t3 = t2*t
246                t4 = t3*t
247    
248                s  = sFld(i,j,k,bi,bj)
249                s3o2 = s*sqrt(s)
250                
251    C     density of freshwater at the surface
252                rfresh =
253         &             eosJMDCFw(1)
254         &           + eosJMDCFw(2)*t
255         &           + eosJMDCFw(3)*t2
256         &           + eosJMDCFw(4)*t3
257         &           + eosJMDCFw(5)*t4
258         &           + eosJMDCFw(6)*t4*t
259    C     density of sea water at the surface
260                rsalt =
261         &         s*(
262         &             eosJMDCSw(1)
263         &           + eosJMDCSw(2)*t
264         &           + eosJMDCSw(3)*t2
265         &           + eosJMDCSw(4)*t3
266         &           + eosJMDCSw(5)*t4
267         &           )
268         &       + s3o2*(
269         &             eosJMDCSw(6)
270         &           + eosJMDCSw(7)*t
271         &           + eosJMDCSw(8)*t2
272         &           )
273         &           + eosJMDCSw(9)*s*s
274    
275                rhoP0(i,j) = rfresh + rsalt
276    
277             end do
278          end do
279    
280          return
281          end
282    
283    C     !ROUTINE: FIND_BULKMOD
284    C     !INTERFACE:
285          subroutine FIND_BULKMOD(
286         I      bi, bj, iMin, iMax, jMin, jMax,  k, kRef,
287         I      tFld, sFld,
288         O      bulkMod,
289         I      myThid )
290    
291    C     !DESCRIPTION: \bv
292    C     *==========================================================*
293    C     | o SUBROUTINE FIND_BULKMOD                                    
294    C     |   Calculates the secant bulk modulus K(S,T,p) of a slice
295    C     *==========================================================*
296    C     |                                                          
297    C     | k    - is the surface level                              
298    C     | kRef - is the level to use to determine the pressure
299    C     |                                                          
300    C     *==========================================================*
301    C     \ev
302    
303    C     !USES:
304          implicit none
305    C     == Global variables ==
306    #include "SIZE.h"
307    #include "EEPARAMS.h"
308    #include "PARAMS.h"
309    #include "EOS.h"
310    
311    C     !INPUT/OUTPUT PARAMETERS:
312    C     == Routine arguments ==
313          integer bi,bj,iMin,iMax,jMin,jMax
314          integer k                 ! Level of surface slice
315          integer kRef              ! Reference pressure level
316          _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
317          _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
318          _RL bulkMod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
319          _RL bMfresh, bMsalt, bMpres
320          integer myThid
321    
322    C     !LOCAL VARIABLES:
323    C     == Local variables ==
324          integer i,j
325          _RL t, t2, t3, t4, s, s3o2, p, p2
326    CEOP
327    
328          do j=jMin,jMax
329             do i=iMin,iMax
330    C     abbreviations
331                t  = tFld(i,j,k,bi,bj)
332                t2 = t*t
333                t3 = t2*t
334                t4 = t3*t
335    
336                s  = sFld(i,j,k,bi,bj)
337                s3o2 = s*sqrt(s)
338    C
339                p = pressure(i,j,kRef,bi,bj)*SItoBar
340                p2 = p*p
341    C     secant bulk modulus of fresh water at the surface
342                bMfresh =
343         &             eosJMDCKFw(1)
344         &           + eosJMDCKFw(2)*t
345         &           + eosJMDCKFw(3)*t2
346         &           + eosJMDCKFw(4)*t3
347         &           + eosJMDCKFw(5)*t4
348    C     secant bulk modulus of sea water at the surface
349                bMsalt =
350         &         s*( eosJMDCKSw(1)
351         &           + eosJMDCKSw(2)*t
352         &           + eosJMDCKSw(3)*t2
353         &           + eosJMDCKSw(4)*t3
354         &           )
355         &    + s3o2*( eosJMDCKSw(5)
356         &           + eosJMDCKSw(6)*t
357         &           + eosJMDCKSw(7)*t2
358         &           )
359    C     secant bulk modulus of sea water at pressure p
360                bMpres =
361         &         p*( eosJMDCKP(1)
362         &           + eosJMDCKP(2)*t
363         &           + eosJMDCKP(3)*t2
364         &           + eosJMDCKP(4)*t3
365         &           )
366         &     + p*s*( eosJMDCKP(5)
367         &           + eosJMDCKP(6)*t
368         &           + eosJMDCKP(7)*t2
369         &           )
370         &      + p*s3o2*eosJMDCKP(8)
371         &      + p2*( eosJMDCKP(9)
372         &           + eosJMDCKP(10)*t
373         &           + eosJMDCKP(11)*t2
374         &           )
375         &    + p2*s*( eosJMDCKP(12)
376         &           + eosJMDCKP(13)*t
377         &           + eosJMDCKP(14)*t2
378         &           )
379    
380                bulkMod(i,j) = bMfresh + bMsalt + bMpres
381    
382             end do
383          end do
384    
385          return
386          end
387    
388    CBOP
389    C     !ROUTINE: FIND_RHONUM
390    C     !INTERFACE:
391          subroutine FIND_RHONUM(
392         I      bi, bj, iMin, iMax, jMin, jMax, k, kRef,
393         I      tFld, sFld,
394         O      rhoNum,
395         I      myThid )
396    
397    C     !DESCRIPTION: \bv
398    C     *==========================================================*
399    C     | o SUBROUTINE FIND_RHONUM
400    C     |   Calculates the numerator of the McDougall et al.
401    C     |   equation of state
402    C     |   - the code is more or less a copy of MOM4
403    C     *==========================================================*
404    C     |                                                          
405    C     | k - is the surface level                              
406    C     | kRef - is the level to use to determine the pressure
407    C     |                                                          
408    C     *==========================================================*
409    C     \ev
410    
411    C     !USES:
412          implicit none
413    C     == Global variables ==
414    #include "SIZE.h"
415    #include "EEPARAMS.h"
416    #include "PARAMS.h"
417    #include "EOS.h"
418    
419    C     !INPUT/OUTPUT PARAMETERS:
420    C     == Routine arguments ==
421          integer bi,bj,iMin,iMax,jMin,jMax
422          integer k         ! Level of surface slice
423          integer kRef              ! Reference pressure level
424          _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
425          _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
426          _RL rhoNum(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
427          integer myThid
428    
429    C     !LOCAL VARIABLES:
430    C     == Local variables ==
431          integer i,j
432          _RL t1, t2, s1, p1
433    CEOP
434          do j=jMin,jMax
435             do i=iMin,iMax
436    C     abbreviations
437                t1  = tFld(i,j,k,bi,bj)
438                t2 = t1*t1
439                s1  = sFld(i,j,k,bi,bj)
440                
441                p1   = pressure(i,j,kRef,bi,bj)*SItodBar
442                
443                rhoNum(i,j) = eosMDJWFnum(0)
444         &           + t1*(eosMDJWFnum(1)
445         &           +     t1*(eosMDJWFnum(2) + eosMDJWFnum(3)*t1) )  
446         &           + s1*(eosMDJWFnum(4)
447         &           +     eosMDJWFnum(5)*t1  + eosMDJWFnum(6)*s1)      
448         &           + p1*(eosMDJWFnum(7) + eosMDJWFnum(8)*t2
449         &           +     eosMDJWFnum(9)*s1
450         &           +     p1*(eosMDJWFnum(10) + eosMDJWFnum(11)*t2) )
451    
452             end do
453          end do
454    
455          return
456          end
457    
458    CBOP
459    C     !ROUTINE: FIND_RHODEN
460    C     !INTERFACE:
461          subroutine FIND_RHODEN(
462         I      bi, bj, iMin, iMax, jMin, jMax, k, kRef,
463         I      tFld, sFld,
464         O      rhoDen,
465         I      myThid )
466    
467    C     !DESCRIPTION: \bv
468    C     *==========================================================*
469    C     | o SUBROUTINE FIND_RHODEN
470    C     |   Calculates the denominator of the McDougall et al.
471    C     |   equation of state
472    C     |   - the code is more or less a copy of MOM4
473    C     *==========================================================*
474    C     |                                                          
475    C     | k - is the surface level                              
476    C     | kRef - is the level to use to determine the pressure
477    C     |                                                          
478    C     *==========================================================*
479    C     \ev
480    
481    C     !USES:
482          implicit none
483    C     == Global variables ==
484    #include "SIZE.h"
485    #include "EEPARAMS.h"
486    #include "PARAMS.h"
487    #include "EOS.h"
488    
489    C     !INPUT/OUTPUT PARAMETERS:
490    C     == Routine arguments ==
491          integer bi,bj,iMin,iMax,jMin,jMax
492          integer k         ! Level of surface slice
493          integer kRef              ! Reference pressure level
494          _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
495          _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
496          _RL rhoDen(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
497          integer myThid
498    
499    C     !LOCAL VARIABLES:
500    C     == Local variables ==
501          integer i,j
502          _RL t1, t2, s1, sp5, p1, p1t1
503          _RL den, epsln
504          parameter ( epsln = 0. _d 0 )
505    CEOP
506          do j=jMin,jMax
507             do i=iMin,iMax
508    C     abbreviations
509                t1  = tFld(i,j,k,bi,bj)
510                t2 = t1*t1
511                s1  = sFld(i,j,k,bi,bj)
512                sp5 = sqrt(s1)
513                
514                p1   = pressure(i,j,kRef,bi,bj)*SItodBar
515                p1t1 = p1*t1
516                
517                den = eosMDJWFden(0)
518         &           + t1*(eosMDJWFden(1)
519         &           +     t1*(eosMDJWFden(2)
520         &           +         t1*(eosMDJWFden(3) + t1*eosMDJWFden(4) ) ) )
521         &           + s1*(eosMDJWFden(5)
522         &           +     t1*(eosMDJWFden(6)
523         &           +         eosMDJWFden(7)*t2)
524         &           +     sp5*(eosMDJWFden(8) + eosMDJWFden(9)*t2) )
525         &           + p1*(eosMDJWFden(10)
526         &           +     p1t1*(eosMDJWFden(11)*t2 + eosMDJWFden(12)*p1) )
527                
528                rhoDen(i,j) = 1.0/(epsln+den)
529    
530             end do
531          end do
532    
533          return
534          end
535    
536          subroutine find_rho_scalar(
537         I     tLoc, sLoc, pLoc,
538         O     rhoLoc,
539         I     myThid )
540    
541    C     !DESCRIPTION: \bv
542    C     *==========================================================*
543    C     | o SUBROUTINE FIND_RHO_SCALAR                                  
544    C     |   Calculates [rho(S,T,p)-rhoConst]
545    C     *==========================================================*
546    C     \ev
547    
548    C     !USES:
549          implicit none
550    C     == Global variables ==
551    #include "SIZE.h"
552    #include "EEPARAMS.h"
553    #include "PARAMS.h"
554    #include "EOS.h"
555    
556    C     !INPUT/OUTPUT PARAMETERS:
557    C     == Routine arguments ==
558          _RL sLoc, tLoc, pLoc
559          _RL rhoLoc
560          integer myThid
561    
562    C     !LOCAL VARIABLES:
563    C     == Local variables ==
564    
565          _RL t1, t2, t3, t4, s1, s3o2, p1, p2, sp5, p1t1
566          _RL rfresh, rsalt, rhoP0
567          _RL bMfresh, bMsalt, bMpres, BulkMod
568          _RL rhoNum, rhoDen, den, epsln
569          parameter ( epsln = 0. _d 0 )
570    
571          character*(max_len_mbuf) msgbuf
572    CEOP
573    
574          rhoLoc  = 0. _d 0
575          rhoP0   = 0. _d 0
576          bulkMod = 0. _d 0
577          rfresh  = 0. _d 0
578          rsalt   = 0. _d 0
579          bMfresh = 0. _d 0
580          bMsalt  = 0. _d 0
581          bMpres  = 0. _d 0
582          rhoNum  = 0. _d 0
583          rhoDen  = 0. _d 0
584          den     = 0. _d 0
585    
586          t1 = tLoc
587          t2 = t1*t1
588          t3 = t2*t1
589          t4 = t3*t1
590          
591          s1  = sLoc
592          if ( s1 .lt. 0. _d 0 ) then
593    C     issue a warning
594             write(*,'(a,i3,a,i3,a,i3,a,e13.5)')
595         &        ' FIND_RHO_SCALAR:   WARNING, salinity = ', s1
596             s1 = 0. _d 0
597          end if
598    
599          if (equationOfState.eq.'LINEAR') then
600    
601             rhoLoc = 0. _d  0
602    
603          elseif (equationOfState.eq.'POLY3') then
604    
605    C     this is not correct, there is a field eosSig0 which should be use here
606    C     but I do not intent to include the reference level in this routine
607             write(*,'(a)')
608         &        ' FIND_RHO_SCALAR: for POLY3, the density is not'
609             write(*,'(a)')
610         &         '                 computed correctly in this routine'
611             rhoLoc = 0. _d 0
612    
613          elseif ( equationOfState(1:5).eq.'JMD95'
614         &      .or. equationOfState.eq.'UNESCO' ) then
615    C     nonlinear equation of state in pressure coordinates
616    
617             s3o2 = s1*sqrt(s1)
618            
619             p1 = pLoc*SItoBar
620             p2 = p1*p1
621    
622    C     density of freshwater at the surface
623             rfresh =
624         &          eosJMDCFw(1)
625         &        + eosJMDCFw(2)*t1
626         &        + eosJMDCFw(3)*t2
627         &        + eosJMDCFw(4)*t3
628         &        + eosJMDCFw(5)*t4
629         &        + eosJMDCFw(6)*t4*t1
630    C     density of sea water at the surface
631             rsalt =
632         &        s1*(
633         &             eosJMDCSw(1)
634         &           + eosJMDCSw(2)*t1
635         &           + eosJMDCSw(3)*t2
636         &           + eosJMDCSw(4)*t3
637         &           + eosJMDCSw(5)*t4
638         &           )
639         &        + s3o2*(
640         &             eosJMDCSw(6)
641         &           + eosJMDCSw(7)*t1
642         &           + eosJMDCSw(8)*t2
643         &           )
644         &           + eosJMDCSw(9)*s1*s1
645    
646             rhoP0 = rfresh + rsalt
647    
648    C     secant bulk modulus of fresh water at the surface
649             bMfresh =
650         &             eosJMDCKFw(1)
651         &           + eosJMDCKFw(2)*t1
652         &           + eosJMDCKFw(3)*t2
653         &           + eosJMDCKFw(4)*t3
654         &           + eosJMDCKFw(5)*t4
655    C     secant bulk modulus of sea water at the surface
656             bMsalt =
657         &        s1*( eosJMDCKSw(1)
658         &           + eosJMDCKSw(2)*t1
659         &           + eosJMDCKSw(3)*t2
660         &           + eosJMDCKSw(4)*t3
661         &           )
662         &    + s3o2*( eosJMDCKSw(5)
663         &           + eosJMDCKSw(6)*t1
664         &           + eosJMDCKSw(7)*t2
665         &           )
666    C     secant bulk modulus of sea water at pressure p
667             bMpres =
668         &        p1*( eosJMDCKP(1)
669         &           + eosJMDCKP(2)*t1
670         &           + eosJMDCKP(3)*t2
671         &           + eosJMDCKP(4)*t3
672         &           )
673         &   + p1*s1*( eosJMDCKP(5)
674         &           + eosJMDCKP(6)*t1
675         &           + eosJMDCKP(7)*t2
676         &           )
677         &      + p1*s3o2*eosJMDCKP(8)
678         &      + p2*( eosJMDCKP(9)
679         &           + eosJMDCKP(10)*t1
680         &           + eosJMDCKP(11)*t2
681         &           )
682         &    + p2*s1*( eosJMDCKP(12)
683         &           + eosJMDCKP(13)*t1
684         &           + eosJMDCKP(14)*t2
685         &           )
686    
687             bulkMod = bMfresh + bMsalt + bMpres
688            
689    C     density of sea water at pressure p
690             rhoLoc = rhoP0/(1. _d 0 - p1/bulkMod) - rhoConst
691    
692          elseif ( equationOfState.eq.'MDJWF' ) then
693    
694             sp5 = sqrt(s1)
695                
696             p1   = pLoc*SItodBar
697             p1t1 = p1*t1
698    
699             rhoNum = eosMDJWFnum(0)
700         &        + t1*(eosMDJWFnum(1)
701         &        +     t1*(eosMDJWFnum(2) + eosMDJWFnum(3)*t1) )  
702         &        + s1*(eosMDJWFnum(4)
703         &        +     eosMDJWFnum(5)*t1  + eosMDJWFnum(6)*s1)      
704         &        + p1*(eosMDJWFnum(7) + eosMDJWFnum(8)*t2
705         &        +     eosMDJWFnum(9)*s1
706         &        +     p1*(eosMDJWFnum(10) + eosMDJWFnum(11)*t2) )
707    
708                
709             den = eosMDJWFden(0)
710         &        + t1*(eosMDJWFden(1)
711         &        +     t1*(eosMDJWFden(2)
712         &        +         t1*(eosMDJWFden(3) + t1*eosMDJWFden(4) ) ) )
713         &        + s1*(eosMDJWFden(5)
714         &        +     t1*(eosMDJWFden(6)
715         &        +         eosMDJWFden(7)*t2)
716         &        +     sp5*(eosMDJWFden(8) + eosMDJWFden(9)*t2) )
717         &        + p1*(eosMDJWFden(10)
718         &        +     p1t1*(eosMDJWFden(11)*t2 + eosMDJWFden(12)*p1) )
719                
720             rhoDen = 1.0/(epsln+den)
721    
722             rhoLoc = rhoNum*rhoDen - rhoConst
723    
724          elseif( equationOfState .eq. 'IDEALG' ) then
725    C    
726        else        else
727         write(0,*) 'FIND_RHO: eqn = ',eqn         write(msgbuf,'(3a)')
728         stop 'FIND_RHO: argument "eqn" has illegal value'       &        ' FIND_RHO_SCALAR : equationOfState = "',
729         &        equationOfState,'"'
730           call print_error( msgbuf, mythid )
731           stop 'ABNORMAL END: S/R FIND_RHO_SCALAR'
732        endif        endif
733    
 ! ------------------------------------------------------------------------------  
734        return        return
735        end        end
736  ! ==============================================================================  
737    CBOP
738    C     !ROUTINE: LOOK_FOR_NEG_SALINITY
739    C     !INTERFACE:
740          subroutine LOOK_FOR_NEG_SALINITY(
741         I     bi, bj, iMin, iMax, jMin, jMax,  k,
742         I     sFld,
743         I     myThid )
744    
745    C     !DESCRIPTION: \bv
746    C     *==========================================================*
747    C     | o SUBROUTINE LOOK_FOR_NEG_SALINITY
748    C     |   looks for and fixes negative salinity values
749    C     |   this is necessary if the equation of state uses
750    C     |   the square root of salinity
751    C     *==========================================================*
752    C     |                                                          
753    C     | k - is the Salt level                              
754    C     |                                                          
755    C     *==========================================================*
756    C     \ev
757    
758    C     !USES:
759          implicit none
760    C     == Global variables ==
761    #include "SIZE.h"
762    #include "EEPARAMS.h"
763    #include "PARAMS.h"
764    #include "GRID.h"
765    
766    C     !INPUT/OUTPUT PARAMETERS:
767    C     == Routine arguments ==
768          integer bi,bj,iMin,iMax,jMin,jMax
769          integer k                 ! Level of Salt slice
770          _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
771          integer myThid
772    
773    C     !LOCAL VARIABLES:
774    C     == Local variables ==
775          integer i,j, localWarning
776          character*(max_len_mbuf) msgbuf
777    CEOP
778    
779          localWarning = 0
780          do j=jMin,jMax
781           do i=iMin,iMax
782    C     abbreviations
783            if ( sFld(i,j,k,bi,bj) .lt. 0. _d 0 ) then
784             localWarning = localWarning + 1
785             sFld(i,j,k,bi,bj) = 0. _d 0
786            end if
787           end do
788          end do
789    C     issue a warning
790          if ( localWarning .gt. 0 ) then
791           write(*,'(a,a)')
792         &      'S/R LOOK_FOR_NEG_SALINITY: found negative salinity',
793         &      'values and reset them to zero.'
794           write(*,'(a,I3)')
795         &      'S/R LOOK_FOR_NEG_SALINITY: current level is k = ', k
796          end if
797    
798          return
799          end
Legend:
Removed from v.1.8  
changed lines
  Added in v.1.23
  ViewVC Help
Powered by ViewVC 1.1.22