/[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.23 by heimbach, Fri Nov 15 03:01:21 2002 UTC revision 1.24 by jmc, Tue Feb 18 15:25:09 2003 UTC
# Line 7  C $Name$ Line 7  C $Name$
7  CBOP  CBOP
8  C     !ROUTINE: FIND_RHO  C     !ROUTINE: FIND_RHO
9  C     !INTERFACE:  C     !INTERFACE:
10        subroutine FIND_RHO(        SUBROUTINE FIND_RHO(
11       I      bi, bj, iMin, iMax, jMin, jMax,  k, kRef,       I      bi, bj, iMin, iMax, jMin, jMax,  k, kRef,
12       I      tFld, sFld,       I      tFld, sFld,
13       O      rholoc,       O      rholoc,
# Line 27  C     *================================= Line 27  C     *=================================
27  C     \ev  C     \ev
28    
29  C     !USES:  C     !USES:
30        implicit none        IMPLICIT NONE
31  C     == Global variables ==  C     == Global variables ==
32  #include "SIZE.h"  #include "SIZE.h"
33  #include "EEPARAMS.h"  #include "EEPARAMS.h"
# Line 37  C     == Global variables == Line 37  C     == Global variables ==
37    
38  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
39  C     == Routine arguments ==  C     == Routine arguments ==
40        integer bi,bj,iMin,iMax,jMin,jMax  C     k    :: Level of Theta/Salt slice
41        integer k                 ! Level of Theta/Salt slice  C     kRef :: Pressure reference level
42        integer kRef              ! Pressure reference level        INTEGER bi,bj,iMin,iMax,jMin,jMax
43          INTEGER k
44          INTEGER kRef
45        _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
46        _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
47        _RL rholoc(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL rholoc(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
48        integer myThid        INTEGER myThid
49    
50  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
51  C     == Local variables ==  C     == Local variables ==
52        integer i,j        INTEGER i,j
53        _RL refTemp,refSalt,sigRef,tP,sP,deltaSig,dRho        _RL refTemp,refSalt,sigRef,tP,sP,deltaSig,dRho
54        _RL rhoP0(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL locPres(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
55          _RL rhoP0  (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
56        _RL bulkMod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL bulkMod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
57        _RL rhoNum(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL rhoNum (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
58        _RL rhoDen(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL rhoDen (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
59        character*(max_len_mbuf) msgbuf        CHARACTER*(MAX_LEN_MBUF) msgbuf
60  CEOP  CEOP
61    
62  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 71  CEOP Line 74  CEOP
74       &     sFld, myThid )       &     sFld, myThid )
75  #endif  #endif
76    
77        if (equationOfState.eq.'LINEAR') then        IF (equationOfState.EQ.'LINEAR') THEN
78    
79  C ***NOTE***  C ***NOTE***
80  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 82  C ********** Line 85  C **********
85    
86         dRho = rhoNil-rhoConst         dRho = rhoNil-rhoConst
87    
88         do j=jMin,jMax         DO j=jMin,jMax
89          do i=iMin,iMax          DO i=iMin,iMax
90           rholoc(i,j)=rhoNil*(           rholoc(i,j)=rhoNil*(
91       &     sBeta*(sFld(i,j,k,bi,bj)-refSalt)       &     sBeta*(sFld(i,j,k,bi,bj)-refSalt)
92       &   -tAlpha*(tFld(i,j,k,bi,bj)-refTemp) )       &   -tAlpha*(tFld(i,j,k,bi,bj)-refTemp) )
93       &        + dRho       &        + dRho
94          enddo          ENDDO
95         enddo         ENDDO
96                
97        elseif (equationOfState.eq.'POLY3') then        ELSEIF (equationOfState.EQ.'POLY3') THEN
98    
99         refTemp=eosRefT(kRef)         refTemp=eosRefT(kRef)
100         refSalt=eosRefS(kRef)         refSalt=eosRefS(kRef)
101         sigRef=eosSig0(kRef) + (1000.-rhoConst)         sigRef=eosSig0(kRef) + (1000.-rhoConst)
102    
103         do j=jMin,jMax         DO j=jMin,jMax
104          do i=iMin,iMax          DO i=iMin,iMax
105           tP=tFld(i,j,k,bi,bj)-refTemp           tP=tFld(i,j,k,bi,bj)-refTemp
106           sP=sFld(i,j,k,bi,bj)-refSalt           sP=sFld(i,j,k,bi,bj)-refSalt
107  #ifdef USE_FACTORIZED_POLY  #ifdef USE_FACTORIZED_POLY
# Line 123  C ********** Line 126  C **********
126       &    +eosC(9,kRef)         *sP*sP*sP       &    +eosC(9,kRef)         *sP*sP*sP
127  #endif  #endif
128           rholoc(i,j)=sigRef+deltaSig           rholoc(i,j)=sigRef+deltaSig
129          enddo          ENDDO
130         enddo         ENDDO
131    
132        elseif ( equationOfState(1:5).eq.'JMD95'        ELSEIF ( equationOfState(1:5).EQ.'JMD95'
133       &      .or. equationOfState.eq.'UNESCO' ) then       &      .OR. equationOfState.EQ.'UNESCO' ) THEN
134  C     nonlinear equation of state in pressure coordinates  C     nonlinear equation of state in pressure coordinates
135    
136             CALL PRESSURE_FOR_EOS(
137         I        bi, bj, iMin, iMax, jMin, jMax,  kRef,
138         O        locPres,
139         I        myThid )
140    
141           CALL FIND_RHOP0(           CALL FIND_RHOP0(
142       I        bi, bj, iMin, iMax, jMin, jMax, k,       I        bi, bj, iMin, iMax, jMin, jMax, k,
143       I        tFld, sFld,       I        tFld, sFld,
# Line 137  C     nonlinear equation of state in pre Line 145  C     nonlinear equation of state in pre
145       I        myThid )       I        myThid )
146                    
147           CALL FIND_BULKMOD(           CALL FIND_BULKMOD(
148       I        bi, bj, iMin, iMax, jMin, jMax,  k, kRef,       I        bi, bj, iMin, iMax, jMin, jMax,  k,
149       I        tFld, sFld,       I        locPres, tFld, sFld,
150       O        bulkMod,       O        bulkMod,
151       I        myThid )       I        myThid )
152                    
153  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
154  cph can't do storing here since find_rho is called multiple times;  cph can't DO storing here since find_rho is called multiple times;
155  cph additional recomp. should be acceptable  cph additional recomp. should be acceptable
156  cphCADJ STORE rhoP0(:,:)   = comlev1_bibj_k ,  key=kkey , byte=isbyte  cphCADJ STORE rhoP0(:,:)   = comlev1_bibj_k ,  key=kkey , byte=isbyte
157  cphCADJ STORE bulkMod(:,:) = comlev1_bibj_k ,  key=kkey , byte=isbyte  cphCADJ STORE bulkMod(:,:) = comlev1_bibj_k ,  key=kkey , byte=isbyte
158  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
159           do j=jMin,jMax           DO j=jMin,jMax
160              do i=iMin,iMax              DO i=iMin,iMax
161    
162  C     density of sea water at pressure p  C     density of sea water at pressure p
163                 rholoc(i,j) = rhoP0(i,j)                 rholoc(i,j) = rhoP0(i,j)
164       &              /(1. _d 0 -       &              /(1. _d 0 -
165       &              pressure(i,j,kRef,bi,bj)*SItoBar/bulkMod(i,j))       &              locPres(i,j)*SItoBar/bulkMod(i,j))
166       &              - rhoConst       &              - rhoConst
167    
168              end do              ENDDO
169           end do           ENDDO
170    
171        elseif ( equationOfState.eq.'MDJWF' ) then        ELSEIF ( equationOfState.EQ.'MDJWF' ) THEN
172    
173             CALL PRESSURE_FOR_EOS(
174         I        bi, bj, iMin, iMax, jMin, jMax,  kRef,
175         O        locPres,
176         I        myThid )
177    
178           CALL FIND_RHONUM( bi, bj, iMin, iMax, jMin, jMax, k, kRef,           CALL FIND_RHONUM( bi, bj, iMin, iMax, jMin, jMax, k,
179       &      tFld, sFld, rhoNum, myThid )       &      locPres, tFld, sFld, rhoNum, myThid )
180    
181           CALL FIND_RHODEN( bi, bj, iMin, iMax, jMin, jMax, k, kRef,           CALL FIND_RHODEN( bi, bj, iMin, iMax, jMin, jMax, k,
182       &      tFld, sFld, rhoDen, myThid )       &      locPres, tFld, sFld, rhoDen, myThid )
183  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
184  cph can't do storing here since find_rho is called multiple times;  cph can't DO storing here since find_rho is called multiple times;
185  cph additional recomp. should be acceptable  cph additional recomp. should be acceptable
186  cphCADJ STORE rhoNum(:,:) = comlev1_bibj_k ,  key=kkey , byte=isbyte  cphCADJ STORE rhoNum(:,:) = comlev1_bibj_k ,  key=kkey , byte=isbyte
187  cphCADJ STORE rhoDen(:,:) = comlev1_bibj_k ,  key=kkey , byte=isbyte  cphCADJ STORE rhoDen(:,:) = comlev1_bibj_k ,  key=kkey , byte=isbyte
188  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
189           do j=jMin,jMax           DO j=jMin,jMax
190              do i=iMin,iMax              DO i=iMin,iMax
   
191                 rholoc(i,j) = rhoNum(i,j)*rhoDen(i,j) - rhoConst                 rholoc(i,j) = rhoNum(i,j)*rhoDen(i,j) - rhoConst
192                ENDDO
193             ENDDO
194    
195              end do        ELSEIF( equationOfState .EQ. 'IDEALG' ) THEN
          end do  
   
       elseif( equationOfState .eq. 'IDEALG' ) then  
196  C      C    
197        else        ELSE
198         write(msgbuf,'(3a)')         WRITE(msgbuf,'(3a)')
199       &        ' FIND_RHO: equationOfState = "',equationOfState,'"'       &        ' FIND_RHO: equationOfState = "',equationOfState,'"'
200         call print_error( msgbuf, mythid )         CALL print_error( msgbuf, mythid )
201         stop 'ABNORMAL END: S/R FIND_RHO'         STOP 'ABNORMAL END: S/R FIND_RHO'
202        endif        ENDIF
203    
204        return        RETURN
205        end        END
206    
207  CBOP  CBOP
208  C     !ROUTINE: FIND_RHOP0  C     !ROUTINE: FIND_RHOP0
209  C     !INTERFACE:  C     !INTERFACE:
210        subroutine FIND_RHOP0(        SUBROUTINE FIND_RHOP0(
211       I      bi, bj, iMin, iMax, jMin, jMax, k,       I      bi, bj, iMin, iMax, jMin, jMax, k,
212       I      tFld, sFld,       I      tFld, sFld,
213       O      rhoP0,       O      rhoP0,
# Line 214  C     *================================= Line 225  C     *=================================
225  C     \ev  C     \ev
226    
227  C     !USES:  C     !USES:
228        implicit none        IMPLICIT NONE
229  C     == Global variables ==  C     == Global variables ==
230  #include "SIZE.h"  #include "SIZE.h"
231  #include "EEPARAMS.h"  #include "EEPARAMS.h"
# Line 223  C     == Global variables == Line 234  C     == Global variables ==
234    
235  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
236  C     == Routine arguments ==  C     == Routine arguments ==
237        integer bi,bj,iMin,iMax,jMin,jMax  C     k    :: Level of Theta/Salt slice
238        integer k         ! Level of surface slice        INTEGER bi,bj,iMin,iMax,jMin,jMax
239          INTEGER k
240        _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
241        _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
242        _RL rhoP0(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL rhoP0(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
243        _RL rfresh, rsalt        _RL rfresh, rsalt
244        integer myThid        INTEGER myThid
245    
246  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
247  C     == Local variables ==  C     == Local variables ==
248        integer i,j        INTEGER i,j
249        _RL t, t2, t3, t4, s, s3o2        _RL t, t2, t3, t4, s, s3o2
250  CEOP  CEOP
251    
252        do j=jMin,jMax        DO j=jMin,jMax
253           do i=iMin,iMax           DO i=iMin,iMax
254  C     abbreviations  C     abbreviations
255              t  = tFld(i,j,k,bi,bj)              t  = tFld(i,j,k,bi,bj)
256              t2 = t*t              t2 = t*t
# Line 246  C     abbreviations Line 258  C     abbreviations
258              t4 = t3*t              t4 = t3*t
259    
260              s  = sFld(i,j,k,bi,bj)              s  = sFld(i,j,k,bi,bj)
261              s3o2 = s*sqrt(s)              s3o2 = s*SQRT(s)
262                            
263  C     density of freshwater at the surface  C     density of freshwater at the surface
264              rfresh =              rfresh =
# Line 274  C     density of sea water at the surfac Line 286  C     density of sea water at the surfac
286    
287              rhoP0(i,j) = rfresh + rsalt              rhoP0(i,j) = rfresh + rsalt
288    
289           end do           ENDDO
290        end do        ENDDO
291    
292        return        RETURN
293        end        END
294    
295  C     !ROUTINE: FIND_BULKMOD  C     !ROUTINE: FIND_BULKMOD
296  C     !INTERFACE:  C     !INTERFACE:
297        subroutine FIND_BULKMOD(        SUBROUTINE FIND_BULKMOD(
298       I      bi, bj, iMin, iMax, jMin, jMax,  k, kRef,       I      bi, bj, iMin, iMax, jMin, jMax,  k,
299       I      tFld, sFld,       I      locPres, tFld, sFld,
300       O      bulkMod,       O      bulkMod,
301       I      myThid )       I      myThid )
302    
# Line 294  C     | o SUBROUTINE FIND_BULKMOD Line 306  C     | o SUBROUTINE FIND_BULKMOD
306  C     |   Calculates the secant bulk modulus K(S,T,p) of a slice  C     |   Calculates the secant bulk modulus K(S,T,p) of a slice
307  C     *==========================================================*  C     *==========================================================*
308  C     |                                                            C     |                                                          
309  C     | k    - is the surface level                                C     | k    - is the level of Theta/Salt slice
 C     | kRef - is the level to use to determine the pressure  
310  C     |                                                            C     |                                                          
311  C     *==========================================================*  C     *==========================================================*
312  C     \ev  C     \ev
313    
314  C     !USES:  C     !USES:
315        implicit none        IMPLICIT NONE
316  C     == Global variables ==  C     == Global variables ==
317  #include "SIZE.h"  #include "SIZE.h"
318  #include "EEPARAMS.h"  #include "EEPARAMS.h"
# Line 310  C     == Global variables == Line 321  C     == Global variables ==
321    
322  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
323  C     == Routine arguments ==  C     == Routine arguments ==
324        integer bi,bj,iMin,iMax,jMin,jMax  C     k    :: Level of Theta/Salt slice
325        integer k                 ! Level of surface slice        INTEGER bi,bj,iMin,iMax,jMin,jMax
326        integer kRef              ! Reference pressure level        INTEGER k
327          INTEGER kRef
328          _RL locPres(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
329        _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
330        _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
331        _RL bulkMod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL bulkMod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
332        _RL bMfresh, bMsalt, bMpres        _RL bMfresh, bMsalt, bMpres
333        integer myThid        INTEGER myThid
334    
335  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
336  C     == Local variables ==  C     == Local variables ==
337        integer i,j        INTEGER i,j
338        _RL t, t2, t3, t4, s, s3o2, p, p2        _RL t, t2, t3, t4, s, s3o2, p, p2
339  CEOP  CEOP
340    
341        do j=jMin,jMax        DO j=jMin,jMax
342           do i=iMin,iMax           DO i=iMin,iMax
343  C     abbreviations  C     abbreviations
344              t  = tFld(i,j,k,bi,bj)              t  = tFld(i,j,k,bi,bj)
345              t2 = t*t              t2 = t*t
# Line 334  C     abbreviations Line 347  C     abbreviations
347              t4 = t3*t              t4 = t3*t
348    
349              s  = sFld(i,j,k,bi,bj)              s  = sFld(i,j,k,bi,bj)
350              s3o2 = s*sqrt(s)              s3o2 = s*SQRT(s)
351  C  C
352              p = pressure(i,j,kRef,bi,bj)*SItoBar              p = locPres(i,j)*SItoBar
353              p2 = p*p              p2 = p*p
354  C     secant bulk modulus of fresh water at the surface  C     secant bulk modulus of fresh water at the surface
355              bMfresh =              bMfresh =
# Line 379  C     secant bulk modulus of sea water a Line 392  C     secant bulk modulus of sea water a
392    
393              bulkMod(i,j) = bMfresh + bMsalt + bMpres              bulkMod(i,j) = bMfresh + bMsalt + bMpres
394    
395           end do           ENDDO
396        end do        ENDDO
397    
398        return        RETURN
399        end        END
400    
401  CBOP  CBOP
402  C     !ROUTINE: FIND_RHONUM  C     !ROUTINE: FIND_RHONUM
403  C     !INTERFACE:  C     !INTERFACE:
404        subroutine FIND_RHONUM(        SUBROUTINE FIND_RHONUM(
405       I      bi, bj, iMin, iMax, jMin, jMax, k, kRef,       I      bi, bj, iMin, iMax, jMin, jMax, k,
406       I      tFld, sFld,       I      locPres, tFld, sFld,
407       O      rhoNum,       O      rhoNum,
408       I      myThid )       I      myThid )
409    
# Line 402  C     |   equation of state Line 415  C     |   equation of state
415  C     |   - the code is more or less a copy of MOM4  C     |   - the code is more or less a copy of MOM4
416  C     *==========================================================*  C     *==========================================================*
417  C     |                                                            C     |                                                          
418  C     | k - is the surface level                                C     | k    - is the level of Theta/Salt slice
 C     | kRef - is the level to use to determine the pressure  
419  C     |                                                            C     |                                                          
420  C     *==========================================================*  C     *==========================================================*
421  C     \ev  C     \ev
422    
423  C     !USES:  C     !USES:
424        implicit none        IMPLICIT NONE
425  C     == Global variables ==  C     == Global variables ==
426  #include "SIZE.h"  #include "SIZE.h"
427  #include "EEPARAMS.h"  #include "EEPARAMS.h"
# Line 418  C     == Global variables == Line 430  C     == Global variables ==
430    
431  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
432  C     == Routine arguments ==  C     == Routine arguments ==
433        integer bi,bj,iMin,iMax,jMin,jMax  C     k    :: Level of Theta/Salt slice
434        integer k         ! Level of surface slice        INTEGER bi,bj,iMin,iMax,jMin,jMax
435        integer kRef              ! Reference pressure level        INTEGER k
436          _RL locPres(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
437        _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
438        _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
439        _RL rhoNum(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL rhoNum(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
440        integer myThid        INTEGER myThid
441    
442  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
443  C     == Local variables ==  C     == Local variables ==
444        integer i,j        INTEGER i,j
445        _RL t1, t2, s1, p1        _RL t1, t2, s1, p1
446  CEOP  CEOP
447        do j=jMin,jMax        DO j=jMin,jMax
448           do i=iMin,iMax           DO i=iMin,iMax
449  C     abbreviations  C     abbreviations
450              t1  = tFld(i,j,k,bi,bj)              t1  = tFld(i,j,k,bi,bj)
451              t2 = t1*t1              t2 = t1*t1
452              s1  = sFld(i,j,k,bi,bj)              s1  = sFld(i,j,k,bi,bj)
453                            
454              p1   = pressure(i,j,kRef,bi,bj)*SItodBar              p1   = locPres(i,j)*SItodBar
455                            
456              rhoNum(i,j) = eosMDJWFnum(0)              rhoNum(i,j) = eosMDJWFnum(0)
457       &           + t1*(eosMDJWFnum(1)       &           + t1*(eosMDJWFnum(1)
# Line 449  C     abbreviations Line 462  C     abbreviations
462       &           +     eosMDJWFnum(9)*s1       &           +     eosMDJWFnum(9)*s1
463       &           +     p1*(eosMDJWFnum(10) + eosMDJWFnum(11)*t2) )       &           +     p1*(eosMDJWFnum(10) + eosMDJWFnum(11)*t2) )
464    
465           end do           ENDDO
466        end do        ENDDO
467    
468        return        RETURN
469        end        end
470    
471  CBOP  CBOP
472  C     !ROUTINE: FIND_RHODEN  C     !ROUTINE: FIND_RHODEN
473  C     !INTERFACE:  C     !INTERFACE:
474        subroutine FIND_RHODEN(        SUBROUTINE FIND_RHODEN(
475       I      bi, bj, iMin, iMax, jMin, jMax, k, kRef,       I      bi, bj, iMin, iMax, jMin, jMax, k,
476       I      tFld, sFld,       I      locPres, tFld, sFld,
477       O      rhoDen,       O      rhoDen,
478       I      myThid )       I      myThid )
479    
# Line 472  C     |   equation of state Line 485  C     |   equation of state
485  C     |   - the code is more or less a copy of MOM4  C     |   - the code is more or less a copy of MOM4
486  C     *==========================================================*  C     *==========================================================*
487  C     |                                                            C     |                                                          
488  C     | k - is the surface level                                C     | k    - is the level of Theta/Salt slice
 C     | kRef - is the level to use to determine the pressure  
489  C     |                                                            C     |                                                          
490  C     *==========================================================*  C     *==========================================================*
491  C     \ev  C     \ev
492    
493  C     !USES:  C     !USES:
494        implicit none        IMPLICIT NONE
495  C     == Global variables ==  C     == Global variables ==
496  #include "SIZE.h"  #include "SIZE.h"
497  #include "EEPARAMS.h"  #include "EEPARAMS.h"
# Line 488  C     == Global variables == Line 500  C     == Global variables ==
500    
501  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
502  C     == Routine arguments ==  C     == Routine arguments ==
503        integer bi,bj,iMin,iMax,jMin,jMax  C     k    :: Level of Theta/Salt slice
504        integer k         ! Level of surface slice        INTEGER bi,bj,iMin,iMax,jMin,jMax
505        integer kRef              ! Reference pressure level        INTEGER k
506          _RL locPres(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
507        _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
508        _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
509        _RL rhoDen(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL rhoDen(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
510        integer myThid        INTEGER myThid
511    
512  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
513  C     == Local variables ==  C     == Local variables ==
514        integer i,j        INTEGER i,j
515        _RL t1, t2, s1, sp5, p1, p1t1        _RL t1, t2, s1, sp5, p1, p1t1
516        _RL den, epsln        _RL den, epsln
517        parameter ( epsln = 0. _d 0 )        parameter ( epsln = 0. _d 0 )
518  CEOP  CEOP
519        do j=jMin,jMax        DO j=jMin,jMax
520           do i=iMin,iMax           DO i=iMin,iMax
521  C     abbreviations  C     abbreviations
522              t1  = tFld(i,j,k,bi,bj)              t1  = tFld(i,j,k,bi,bj)
523              t2 = t1*t1              t2 = t1*t1
524              s1  = sFld(i,j,k,bi,bj)              s1  = sFld(i,j,k,bi,bj)
525              sp5 = sqrt(s1)              sp5 = SQRT(s1)
526                            
527              p1   = pressure(i,j,kRef,bi,bj)*SItodBar              p1   = locPres(i,j)*SItodBar
528              p1t1 = p1*t1              p1t1 = p1*t1
529                            
530              den = eosMDJWFden(0)              den = eosMDJWFden(0)
# Line 527  C     abbreviations Line 540  C     abbreviations
540                            
541              rhoDen(i,j) = 1.0/(epsln+den)              rhoDen(i,j) = 1.0/(epsln+den)
542    
543           end do           ENDDO
544        end do        ENDDO
545    
546        return        RETURN
547        end        end
548    
549        subroutine find_rho_scalar(        SUBROUTINE find_rho_scalar(
550       I     tLoc, sLoc, pLoc,       I     tLoc, sLoc, pLoc,
551       O     rhoLoc,       O     rhoLoc,
552       I     myThid )       I     myThid )
# Line 546  C     *================================= Line 559  C     *=================================
559  C     \ev  C     \ev
560    
561  C     !USES:  C     !USES:
562        implicit none        IMPLICIT NONE
563  C     == Global variables ==  C     == Global variables ==
564  #include "SIZE.h"  #include "SIZE.h"
565  #include "EEPARAMS.h"  #include "EEPARAMS.h"
# Line 557  C     !INPUT/OUTPUT PARAMETERS: Line 570  C     !INPUT/OUTPUT PARAMETERS:
570  C     == Routine arguments ==  C     == Routine arguments ==
571        _RL sLoc, tLoc, pLoc        _RL sLoc, tLoc, pLoc
572        _RL rhoLoc        _RL rhoLoc
573        integer myThid        INTEGER myThid
574    
575  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
576  C     == Local variables ==  C     == Local variables ==
# Line 589  CEOP Line 602  CEOP
602        t4 = t3*t1        t4 = t3*t1
603                
604        s1  = sLoc        s1  = sLoc
605        if ( s1 .lt. 0. _d 0 ) then        IF ( s1 .LT. 0. _d 0 ) THEN
606  C     issue a warning  C     issue a warning
607           write(*,'(a,i3,a,i3,a,i3,a,e13.5)')           WRITE(*,'(a,i3,a,i3,a,i3,a,e13.5)')
608       &        ' FIND_RHO_SCALAR:   WARNING, salinity = ', s1       &        ' FIND_RHO_SCALAR:   WARNING, salinity = ', s1
609           s1 = 0. _d 0           s1 = 0. _d 0
610        end if        ENDIF
611    
612        if (equationOfState.eq.'LINEAR') then        IF (equationOfState.EQ.'LINEAR') THEN
613    
614           rhoLoc = 0. _d  0           rhoLoc = 0. _d  0
615    
616        elseif (equationOfState.eq.'POLY3') then        ELSEIF (equationOfState.EQ.'POLY3') THEN
617    
618  C     this is not correct, there is a field eosSig0 which should be use here  C     this is not correct, there is a field eosSig0 which should be use here
619  C     but I do not intent to include the reference level in this routine  C     but I DO not intent to include the reference level in this routine
620           write(*,'(a)')           WRITE(*,'(a)')
621       &        ' FIND_RHO_SCALAR: for POLY3, the density is not'       &        ' FIND_RHO_SCALAR: for POLY3, the density is not'
622           write(*,'(a)')           WRITE(*,'(a)')
623       &         '                 computed correctly in this routine'       &         '                 computed correctly in this routine'
624           rhoLoc = 0. _d 0           rhoLoc = 0. _d 0
625    
626        elseif ( equationOfState(1:5).eq.'JMD95'        ELSEIF ( equationOfState(1:5).EQ.'JMD95'
627       &      .or. equationOfState.eq.'UNESCO' ) then       &      .OR. equationOfState.EQ.'UNESCO' ) THEN
628  C     nonlinear equation of state in pressure coordinates  C     nonlinear equation of state in pressure coordinates
629    
630           s3o2 = s1*sqrt(s1)           s3o2 = s1*SQRT(s1)
631                    
632           p1 = pLoc*SItoBar           p1 = pLoc*SItoBar
633           p2 = p1*p1           p2 = p1*p1
# Line 689  C     secant bulk modulus of sea water a Line 702  C     secant bulk modulus of sea water a
702  C     density of sea water at pressure p  C     density of sea water at pressure p
703           rhoLoc = rhoP0/(1. _d 0 - p1/bulkMod) - rhoConst           rhoLoc = rhoP0/(1. _d 0 - p1/bulkMod) - rhoConst
704    
705        elseif ( equationOfState.eq.'MDJWF' ) then        ELSEIF ( equationOfState.EQ.'MDJWF' ) THEN
706    
707           sp5 = sqrt(s1)           sp5 = SQRT(s1)
708                            
709           p1   = pLoc*SItodBar           p1   = pLoc*SItodBar
710           p1t1 = p1*t1           p1t1 = p1*t1
# Line 721  C     density of sea water at pressure p Line 734  C     density of sea water at pressure p
734    
735           rhoLoc = rhoNum*rhoDen - rhoConst           rhoLoc = rhoNum*rhoDen - rhoConst
736    
737        elseif( equationOfState .eq. 'IDEALG' ) then        ELSEIF( equationOfState .EQ. 'IDEALG' ) THEN
738  C      C    
739        else        ELSE
740         write(msgbuf,'(3a)')         WRITE(msgbuf,'(3A)')
741       &        ' FIND_RHO_SCALAR : equationOfState = "',       &        ' FIND_RHO_SCALAR : equationOfState = "',
742       &        equationOfState,'"'       &        equationOfState,'"'
743         call print_error( msgbuf, mythid )         CALL PRINT_ERROR( msgbuf, mythid )
744         stop 'ABNORMAL END: S/R FIND_RHO_SCALAR'         STOP 'ABNORMAL END: S/R FIND_RHO_SCALAR'
745        endif        ENDIF
746    
747        return        RETURN
748        end        END
749    
750  CBOP  CBOP
751  C     !ROUTINE: LOOK_FOR_NEG_SALINITY  C     !ROUTINE: LOOK_FOR_NEG_SALINITY
752  C     !INTERFACE:  C     !INTERFACE:
753        subroutine LOOK_FOR_NEG_SALINITY(        SUBROUTINE LOOK_FOR_NEG_SALINITY(
754       I     bi, bj, iMin, iMax, jMin, jMax,  k,       I     bi, bj, iMin, iMax, jMin, jMax,  k,
755       I     sFld,       I     sFld,
756       I     myThid )       I     myThid )
# Line 746  C     !DESCRIPTION: \bv Line 759  C     !DESCRIPTION: \bv
759  C     *==========================================================*  C     *==========================================================*
760  C     | o SUBROUTINE LOOK_FOR_NEG_SALINITY  C     | o SUBROUTINE LOOK_FOR_NEG_SALINITY
761  C     |   looks for and fixes negative salinity values  C     |   looks for and fixes negative salinity values
762  C     |   this is necessary if the equation of state uses  C     |   this is necessary IF the equation of state uses
763  C     |   the square root of salinity  C     |   the square root of salinity
764  C     *==========================================================*  C     *==========================================================*
765  C     |                                                            C     |                                                          
# Line 756  C     *================================= Line 769  C     *=================================
769  C     \ev  C     \ev
770    
771  C     !USES:  C     !USES:
772        implicit none        IMPLICIT NONE
773  C     == Global variables ==  C     == Global variables ==
774  #include "SIZE.h"  #include "SIZE.h"
775  #include "EEPARAMS.h"  #include "EEPARAMS.h"
# Line 765  C     == Global variables == Line 778  C     == Global variables ==
778    
779  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
780  C     == Routine arguments ==  C     == Routine arguments ==
781        integer bi,bj,iMin,iMax,jMin,jMax  C     k    :: Level of Theta/Salt slice
782        integer k                 ! Level of Salt slice        INTEGER bi,bj,iMin,iMax,jMin,jMax
783          INTEGER k
784        _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
785        integer myThid        INTEGER myThid
786    
787  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
788  C     == Local variables ==  C     == Local variables ==
789        integer i,j, localWarning        INTEGER i,j, localWarning
790        character*(max_len_mbuf) msgbuf        character*(max_len_mbuf) msgbuf
791  CEOP  CEOP
792    
793        localWarning = 0        localWarning = 0
794        do j=jMin,jMax        DO j=jMin,jMax
795         do i=iMin,iMax         DO i=iMin,iMax
796  C     abbreviations  C     abbreviations
797          if ( sFld(i,j,k,bi,bj) .lt. 0. _d 0 ) then          IF ( sFld(i,j,k,bi,bj) .LT. 0. _d 0 ) THEN
798           localWarning = localWarning + 1           localWarning = localWarning + 1
799           sFld(i,j,k,bi,bj) = 0. _d 0           sFld(i,j,k,bi,bj) = 0. _d 0
800          end if          ENDIF
801         end do         ENDDO
802        end do        ENDDO
803  C     issue a warning  C     issue a warning
804        if ( localWarning .gt. 0 ) then        IF ( localWarning .GT. 0 ) THEN
805         write(*,'(a,a)')         WRITE(standardMessageUnit,'(A,A)')
806       &      'S/R LOOK_FOR_NEG_SALINITY: found negative salinity',       &      'S/R LOOK_FOR_NEG_SALINITY: found negative salinity',
807       &      'values and reset them to zero.'       &      'values and reset them to zero.'
808         write(*,'(a,I3)')         WRITE(standardMessageUnit,'(A,I3)')
809       &      'S/R LOOK_FOR_NEG_SALINITY: current level is k = ', k       &      'S/R LOOK_FOR_NEG_SALINITY: current level is k = ', k
810        end if        ENDIF
811    
812        return        RETURN
813        end        END
Legend:
Removed from v.1.23  
changed lines
  Added in v.1.24
  ViewVC Help
Powered by ViewVC 1.1.22