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

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

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

revision 1.3 by mlosch, Tue Aug 20 14:45:18 2002 UTC revision 1.19 by mlosch, Thu Aug 21 14:25:34 2008 UTC
# Line 2  C $Header$ Line 2  C $Header$
2  C $Name$  C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
5    #undef INCLUDE_EOS_CHECK
6    
7    C--  File ini_eos.F: Routines to initialise Equation of State
8    C--   Contents
9    C--   o INI_EOS
10    C--   o EOS_CHECK
11    
12    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
13  CBOP  CBOP
14  C     !ROUTINE: INI_EOS  C     !ROUTINE: INI_EOS
15  C     !INTERFACE:  C     !INTERFACE:
16        subroutine ini_eos( myThid )        SUBROUTINE INI_EOS( myThid )
17  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
18  C     *==========================================================*  C     *==========================================================*
19  C     | SUBROUTINE INI_EOS                                        C     | SUBROUTINE INI_EOS
20  C     | o Initialise coefficients of equation of state.  C     | o Initialise coefficients of equation of state.
21  C     *==========================================================*  C     *==========================================================*
22  C     \ev  C     \ev
23    
24  C     !USES:  C     !USES:
25    
26        implicit none        IMPLICIT NONE
27  C     == Global variables ==  C     == Global variables ==
28  #include "SIZE.h"  #include "SIZE.h"
29  #include "EEPARAMS.h"  #include "EEPARAMS.h"
30  #include "PARAMS.h"  #include "PARAMS.h"
31  #include "EOS.h"  #include "EOS.h"
 #include "GRID.h"  
32    
33  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
34  C     == Routine arguments ==  C     == Routine arguments ==
# Line 31  C     myThid -  Number of this instance Line 37  C     myThid -  Number of this instance
37    
38  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
39  C     == Local variables ==  C     == Local variables ==
40  C     bi,bj  - Loop counters  C     i,k    :: Loop counters
41  C     I,J,K        INTEGER  i, k
       INTEGER bi, bj  
       INTEGER  I,  J,  K  
42        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
43          
44          IF ( .NOT.fluidIsWater ) RETURN
45    
46          _BARRIER
47          _BEGIN_MASTER(myThid)
48    
49        equationOfState = eosType        equationOfState = eosType
50    
51        if ( equationOfState .eq. 'LINEAR' ) then        DO k = 1,6
52           if ( tAlpha .eq. UNSET_RL ) tAlpha = 2.D-4           eosJMDCFw(k) = 0. _d 0
53           if ( sBeta  .eq. UNSET_RL ) sBeta  = 7.4D-4        ENDDO
54        elseif ( equationOfState .eq. 'POLY3' ) then        DO k = 1,9
55             eosJMDCSw(k) = 0. _d 0
56          ENDDO
57          DO k = 1,5
58             eosJMDCKFw(k) = 0. _d 0
59          ENDDO
60          DO k = 1,7
61             eosJMDCKSw(k) = 0. _d 0
62          ENDDO
63          DO k = 1,14
64             eosJMDCKP(k) = 0. _d 0
65          ENDDO
66          DO k = 0,11
67             eosMDJWFnum(k) = 0. _d 0
68          ENDDO
69          DO k = 0,12
70             eosMDJWFden(k) = 0. _d 0
71          ENDDO
72    
73          IF ( equationOfState .EQ. 'LINEAR' ) THEN
74             IF ( tAlpha .EQ. UNSET_RL ) tAlpha = 2.  _d -4
75             IF ( sBeta  .EQ. UNSET_RL ) sBeta  = 7.4 _d -4
76          ELSEIF ( equationOfState .EQ. 'POLY3' ) THEN
77           OPEN(37,FILE='POLY3.COEFFS',STATUS='OLD',FORM='FORMATTED')           OPEN(37,FILE='POLY3.COEFFS',STATUS='OLD',FORM='FORMATTED')
78           READ(37,*) I           READ(37,*) k
79           IF (I.NE.Nr) THEN           IF (k.NE.Nr) THEN
80              WRITE(msgBuf,'(A)')              WRITE(msgBuf,'(A)')
81       &           'ini_eos: attempt to read POLY3.COEFFS failed'       &           'ini_eos: attempt to read POLY3.COEFFS failed'
82              CALL PRINT_ERROR( msgBuf , 1)              CALL PRINT_ERROR( msgBuf , myThid )
83              WRITE(msgBuf,'(A)')              WRITE(msgBuf,'(A)')
84       &           '           because bad # of levels in data'       &           '           because bad # of levels in data'
85              CALL PRINT_ERROR( msgBuf , 1)              CALL PRINT_ERROR( msgBuf , myThid )
86              STOP 'Bad data in POLY3.COEFFS'              STOP 'Bad data in POLY3.COEFFS'
87           ENDIF           ENDIF
88           READ(37,*) (eosRefT(K),eosRefS(K),eosSig0(K),K=1,Nr)           READ(37,*) (eosRefT(k),eosRefS(k),eosSig0(k),k=1,Nr)
89           DO K=1,Nr           DO k=1,Nr
90              READ(37,*) (eosC(I,K),I=1,9)              READ(37,*) (eosC(i,k),i=1,9)
91           ENDDO           ENDDO
92           CLOSE(37)           CLOSE(37)
93    
94        elseif ( equationOfState(1:5) .eq. 'JMD95'        ELSEIF ( equationOfState .EQ. 'JMD95Z'
95       &         .or. equationOfState .eq. 'UNESCO' ) then       &    .OR. equationOfState .EQ. 'JMD95P'
96         &    .OR. equationOfState .EQ. 'UNESCO' ) THEN
97  C  C
98  C     Jackett & McDougall (1995, JPO) equation of state  C     Jackett & McDougall (1995, JPO) equation of state
99  C     rho = R(salinity, potential temperature, pressure)  C     rho = R(salinity, potential temperature, pressure)
100  C     pressure needs to be available (from the previous  C     pressure needs to be available (from the previous
101  C     time step to linearize the problem)  C     time step to linearize the problem)
102  C  C
103           if ( equationOfState .eq. 'JMD95Z'           IF ( equationOfState .EQ. 'JMD95Z' .AND. usingPCoords ) THEN
104       &        .and. buoyancyRelation .eq. 'OCEANICP' ) then              WRITE(msgBuf,'(A)')
             write(msgBuf,'(A)')  
105       &      'ini_eos: equation of state ''JMD95Z'' should not'       &      'ini_eos: equation of state ''JMD95Z'' should not'
106              CALL PRINT_ERROR( msgBuf , 1)              CALL PRINT_ERROR( msgBuf , myThid )
107              write(msgBuf,'(A)')              WRITE(msgBuf,'(A)')
108       &      '         be used together with pressure coordinates.'       &      '         be used together with pressure coordinates.'
109              CALL PRINT_ERROR( msgBuf , 1)              CALL PRINT_ERROR( msgBuf , myThid )
110              write(msgBuf,'(A)')              WRITE(msgBuf,'(A)')
111       &      '         Use only ''JMD95P'' with ''OCEANICP''.'       &      '         Use only ''JMD95P'' with ''OCEANICP''.'
112              CALL PRINT_ERROR( msgBuf , 1)              CALL PRINT_ERROR( msgBuf , myThid )
113              STOP 'ABNORMAL END: S/R INI_EOS'              STOP 'ABNORMAL END: S/R INI_EOS'
114           endif           ENDIF
 C      
          if ( buoyancyRelation .eq. 'OCEANIC' ) then  
             do bj = myByLo(myThid), myByHi(myThid)  
                do bi = myBxLo(myThid), myBxHi(myThid)  
                   do K=1,Nr  
                      do J=1-Oly,sNy+Oly  
                         do I=1-Olx,sNx+Olx  
                            pressure(i,j,k,bi,bj) = rhonil * (  
      &                          - gravity*rC(k)  
      &                          )  
                         end do  
                      end do  
                   end do  
                end do  
             end do  
          elseif ( buoyancyRelation .eq. 'ATMOSPHERIC'  
      &           .or. buoyancyRelation .eq. 'OCEANICP' ) then  
 C     in pressure coordinates the pressure is just the coordinate of  
 C     the tracer point  
             do bj = myByLo(myThid), myByHi(myThid)  
                do bi = myBxLo(myThid), myBxHi(myThid)  
                   do K=1,Nr  
                      do J=1-Oly,sNy+Oly  
                         do I=1-Olx,sNx+Olx  
                            pressure(i,j,k,bi,bj) = rC(k)  
                         end do  
                      end do  
                   end do  
                end do  
             end do  
          endif  
115    
116  C     coefficients nonlinear equation of state in pressure coordinates for  C     coefficients nonlinear equation of state in pressure coordinates for
117  C     1. density of fresh water at p = 0  C     1. density of fresh water at p = 0
118           eosJMDCFw(1) =  999.842594           eosJMDCFw(1) =  999.842594 _d +00
119           eosJMDCFw(2) =    6.793952 _d -02           eosJMDCFw(2) =    6.793952 _d -02
120           eosJMDCFw(3) = -  9.095290 _d -03           eosJMDCFw(3) = -  9.095290 _d -03
121           eosJMDCFw(4) =    1.001685 _d -04           eosJMDCFw(4) =    1.001685 _d -04
# Line 125  C     1. density of fresh water at p = 0 Line 124  C     1. density of fresh water at p = 0
124  C     2. density of sea water at p = 0  C     2. density of sea water at p = 0
125           eosJMDCSw(1) =    8.24493  _d -01           eosJMDCSw(1) =    8.24493  _d -01
126           eosJMDCSw(2) = -  4.0899   _d -03           eosJMDCSw(2) = -  4.0899   _d -03
127           eosJMDCSw(3) =    7.6438   _d -05           eosJMDCSw(3) =    7.6438   _d -05
128           eosJMDCSw(4) = -  8.2467   _d -07           eosJMDCSw(4) = -  8.2467   _d -07
129           eosJMDCSw(5) =    5.3875   _d -09           eosJMDCSw(5) =    5.3875   _d -09
130           eosJMDCSw(6) = -  5.72466  _d -03           eosJMDCSw(6) = -  5.72466  _d -03
131           eosJMDCSw(7) =    1.0227   _d -04           eosJMDCSw(7) =    1.0227   _d -04
132           eosJMDCSw(8) = -  1.6546   _d -06           eosJMDCSw(8) = -  1.6546   _d -06
133           eosJMDCSw(9) =    4.8314   _d -04           eosJMDCSw(9) =    4.8314   _d -04
134           if ( equationOfState(1:5) .eq. 'JMD95' ) then           IF ( equationOfState(1:5) .EQ. 'JMD95' ) THEN
135  C     3. secant bulk modulus K of fresh water at p = 0  C     3. secant bulk modulus K of fresh water at p = 0
136              eosJMDCKFw(1) =   1.965933 _d +04              eosJMDCKFw(1) =   1.965933 _d +04
137              eosJMDCKFw(2) =   1.444304 _d +02              eosJMDCKFw(2) =   1.444304 _d +02
# Line 150  C     4. secant bulk modulus K of sea wa Line 149  C     4. secant bulk modulus K of sea wa
149  C     5. secant bulk modulus K of sea water at p  C     5. secant bulk modulus K of sea water at p
150              eosJMDCKP( 1) =   3.186519 _d +00              eosJMDCKP( 1) =   3.186519 _d +00
151              eosJMDCKP( 2) =   2.212276 _d -02              eosJMDCKP( 2) =   2.212276 _d -02
152              eosJMDCKP( 3) = - 2.984642 _d -04              eosJMDCKP( 3) = - 2.984642 _d -04
153              eosJMDCKP( 4) =   1.956415 _d -06              eosJMDCKP( 4) =   1.956415 _d -06
154              eosJMDCKP( 5) =   6.704388 _d -03              eosJMDCKP( 5) =   6.704388 _d -03
155              eosJMDCKP( 6) = - 1.847318 _d -04              eosJMDCKP( 6) = - 1.847318 _d -04
# Line 163  C     5. secant bulk modulus K of sea wa Line 162  C     5. secant bulk modulus K of sea wa
162              eosJMDCKP(13) =   6.128773 _d -08              eosJMDCKP(13) =   6.128773 _d -08
163              eosJMDCKP(14) =   6.207323 _d -10              eosJMDCKP(14) =   6.207323 _d -10
164    
165           elseif ( equationOfState .eq. 'UNESCO' ) then           ELSEIF ( equationOfState .EQ. 'UNESCO' ) THEN
166    
167              write(msgBuf,'(a)')              WRITE(msgBuf,'(a)')
168       &           'WARNING WARNING WARNING WARNING WARNING WARNING '       &           'WARNING WARNING WARNING WARNING WARNING WARNING '
169              CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,              CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
170       &           SQUEEZE_RIGHT , 1)       &           SQUEEZE_RIGHT , myThid )
171              write(msgBuf,'(a,a)')              WRITE(msgBuf,'(a,a)')
172       &           'WARNING: using the UNESCO formula with potential ',       &           'WARNING: using the UNESCO formula with potential ',
173       &           'temperature'       &           'temperature'
174              CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,              CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
175       &           SQUEEZE_RIGHT , 1)       &           SQUEEZE_RIGHT , myThid )
176              write(msgBuf,'(a)')              WRITE(msgBuf,'(a)')
177       &           'WARNING: can result in density errors of up to 5%'       &           'WARNING: can result in density errors of up to 5%'
178              CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,              CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
179       &           SQUEEZE_RIGHT , 1)       &           SQUEEZE_RIGHT , myThid )
180              write(msgBuf,'(a)')              WRITE(msgBuf,'(a)')
181       &           'WARNING: (see Jackett and McDougall 1995, JPO)'       &           'WARNING: (see Jackett and McDougall 1995, JAOT)'
182              CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,              CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
183       &           SQUEEZE_RIGHT , 1)       &           SQUEEZE_RIGHT , myThid )
184              write(msgBuf,'(a)')              WRITE(msgBuf,'(a)')
185       &           'WARNING WARNING WARNING WARNING WARNING WARNING '       &           'WARNING WARNING WARNING WARNING WARNING WARNING '
186              CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,              CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
187       &           SQUEEZE_RIGHT , 1)       &           SQUEEZE_RIGHT , myThid )
188    
189  C     3. secant bulk modulus K of fresh water at p = 0  C     3. secant bulk modulus K of fresh water at p = 0
190              eosJMDCKFw(1) =   1.965221 _d +04              eosJMDCKFw(1) =   1.965221 _d +04
# Line 195  C     3. secant bulk modulus K of fresh Line 194  C     3. secant bulk modulus K of fresh
194              eosJMDCKFw(5) = - 5.155288 _d -05              eosJMDCKFw(5) = - 5.155288 _d -05
195  C     4. secant bulk modulus K of sea water at p = 0  C     4. secant bulk modulus K of sea water at p = 0
196              eosJMDCKSw(1) =   5.46746  _d +01              eosJMDCKSw(1) =   5.46746  _d +01
197              eosJMDCKSw(2) = - 0.603459 _d +00              eosJMDCKSw(2) = - 0.603459 _d +00
198              eosJMDCKSw(3) =   1.09987  _d -02              eosJMDCKSw(3) =   1.09987  _d -02
199              eosJMDCKSw(4) = - 6.1670   _d -05              eosJMDCKSw(4) = - 6.1670   _d -05
200              eosJMDCKSw(5) =   7.944    _d -02              eosJMDCKSw(5) =   7.944    _d -02
# Line 213  C     5. secant bulk modulus K of sea wa Line 212  C     5. secant bulk modulus K of sea wa
212              eosJMDCKP( 9) =   8.50935  _d -05              eosJMDCKP( 9) =   8.50935  _d -05
213              eosJMDCKP(10) = - 6.12293  _d -06              eosJMDCKP(10) = - 6.12293  _d -06
214              eosJMDCKP(11) =   5.2787   _d -08              eosJMDCKP(11) =   5.2787   _d -08
215              eosJMDCKP(12) = - 9.9348   _d -07              eosJMDCKP(12) = - 9.9348   _d -07
216              eosJMDCKP(13) =   2.0816   _d -08              eosJMDCKP(13) =   2.0816   _d -08
217              eosJMDCKP(14) =   9.1697   _d -10              eosJMDCKP(14) =   9.1697   _d -10
218           else           ELSE
219              STOP 'INI_EOS: We should never reach this point!'              STOP 'INI_EOS: We should never reach this point!'
220           endif           ENDIF
221    
222          ELSEIF ( equationOfState .EQ. 'MDJWF' ) THEN
223    
224             eosMDJWFnum( 0) =  9.99843699 _d +02
225             eosMDJWFnum( 1) =  7.35212840 _d +00
226             eosMDJWFnum( 2) = -5.45928211 _d -02
227             eosMDJWFnum( 3) =  3.98476704 _d -04
228             eosMDJWFnum( 4) =  2.96938239 _d +00
229             eosMDJWFnum( 5) = -7.23268813 _d -03
230             eosMDJWFnum( 6) =  2.12382341 _d -03
231             eosMDJWFnum( 7) =  1.04004591 _d -02
232             eosMDJWFnum( 8) =  1.03970529 _d -07
233             eosMDJWFnum( 9) =  5.18761880 _d -06
234             eosMDJWFnum(10) = -3.24041825 _d -08
235             eosMDJWFnum(11) = -1.23869360 _d -11
236    
237    
238             eosMDJWFden( 0) =  1.00000000 _d +00
239             eosMDJWFden( 1) =  7.28606739 _d -03
240             eosMDJWFden( 2) = -4.60835542 _d -05
241             eosMDJWFden( 3) =  3.68390573 _d -07
242             eosMDJWFden( 4) =  1.80809186 _d -10
243             eosMDJWFden( 5) =  2.14691708 _d -03
244             eosMDJWFden( 6) = -9.27062484 _d -06
245             eosMDJWFden( 7) = -1.78343643 _d -10
246             eosMDJWFden( 8) =  4.76534122 _d -06
247             eosMDJWFden( 9) =  1.63410736 _d -09
248             eosMDJWFden(10) =  5.30848875 _d -06
249             eosMDJWFden(11) = -3.03175128 _d -16
250             eosMDJWFden(12) = -1.27934137 _d -17
251    
252        elseif( equationOfState .eq. 'IDEALG' ) then        ELSEIF( equationOfState .EQ. 'IDEALG' ) THEN
 C      
       else  
253    
254           write(msgbuf,'(3a)') ' INI_EOS: equationOfState = "',        ELSE
255    
256             WRITE(msgbuf,'(3A)') ' INI_EOS: equationOfState = "',
257       &        equationOfState,'"'       &        equationOfState,'"'
258           call print_error( msgbuf, mythid )           CALL PRINT_ERROR( msgbuf, myThid )
259           stop 'ABNORMAL END: S/R INI_EOS'           STOP 'ABNORMAL END: S/R INI_EOS'
260            
261        end if        ENDIF
262    
263    C--   Check EOS initialisation:
264    
265        call check_eos( myThid )        CALL EOS_CHECK( myThid )
266    
267        return        _END_MASTER( myThid )
268        end        _BARRIER
269    
270          RETURN
271          END
272    
273  CBOP  CBOP
274  C     !ROUTINE: CHECK_EOS  C     !ROUTINE: EOS_CHECK
275  C     !INTERFACE:  C     !INTERFACE:
276        subroutine check_eos( myThid )        SUBROUTINE EOS_CHECK( myThid )
277  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
278  C     *==========================================================*  C     *==========================================================*
279  C     | SUBROUTINE CHECK_EOS                                        C     | SUBROUTINE EOS_CHECK
280  C     | o check the equation of state.  C     | o check the equation of state.
281  C     *==========================================================*  C     *==========================================================*
282  C     \ev  C     \ev
283    
284  C     !USES:  C     !USES:
285    
286        implicit none        IMPLICIT NONE
287  #include "SIZE.h"  #include "SIZE.h"
288  #include "EEPARAMS.h"  #include "EEPARAMS.h"
289  #include "PARAMS.h"  #include "PARAMS.h"
290  #include "EOS.h"  #include "EOS.h"
291    #include "GRID.h"
292    #include "DYNVARS.h"
293    
294  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
295  C     == Routine arguments ==  C     == Routine arguments ==
296  C     myThid -  Number of this instance of CHECK_EOS  C     myThid ::  Number of this instance of EOS_CHECK
297        INTEGER myThid        INTEGER myThid
298    
299    #ifdef INCLUDE_EOS_CHECK
300  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
301  C     == Local variables ==  C     == Local variables ==
302  C     bi,bj  - Loop counters  C     bi,bj  - Loop counters
303  C     I,J,K  C     i,j,k
304        INTEGER bi, bj        INTEGER bi, bj
305        INTEGER imin, imax, jmin, jmax        INTEGER iMin, iMax, jMin, jMax
306        INTEGER  I,  J,  K        INTEGER  i, j, k
307        _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL tFld   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
308        _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL sFld   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
309        _RL rhoP0  (1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL pFld   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
310          _RL rhoLoc (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
311        _RL bulkMod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL bulkMod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
       _RL rholoc, psave  
312    
313        INTEGER ncheck, kcheck        INTEGER ncheck, kcheck
314        PARAMETER ( ncheck = 10 )        PARAMETER ( ncheck = 13 )
315        _RL tloc(ncheck), ptloc(ncheck), sloc(ncheck), ploc(ncheck)        _RL tLoc(ncheck), ptLoc(ncheck), sLoc(ncheck), pLoc(ncheck)
316        _RL rloc(ncheck), bloc(ncheck)        _RL rLoc(ncheck), bLoc(ncheck)
317          _RS mskSave
318          _RS rC_Save
319    
320        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
321    
322        DATA tloc        DATA tLoc
323       &     /3.25905152915860 _d 0, 20.38687090048638 _d 0,       &     /3.25905152915860 _d 0, 20.38687090048638 _d 0,
324         &     25.44820830309568 _d 0, 20.17368557065936 _d 0,
325         &     13.43397459640398 _d 0,
326       &      5.               _d 0, 25.               _d 0,       &      5.               _d 0, 25.               _d 0,
327       &      5.               _d 0, 25.               _d 0,       &      5.               _d 0, 25.               _d 0,
328       &      5.               _d 0, 25.               _d 0,       &      5.               _d 0, 25.               _d 0,
329       &      5.               _d 0, 25.               _d 0/,       &      5.               _d 0, 25.               _d 0/,
330       &     ptloc       &     ptLoc
331       &     /3.               _d 0, 20.               _d 0,       &     /3.               _d 0, 20.               _d 0,
332         &     25.               _d 0, 20.               _d 0,
333         &     12.               _d 0,
334       &      5.               _d 0, 25.               _d 0,       &      5.               _d 0, 25.               _d 0,
335       &      5.               _d 0, 25.               _d 0,       &      5.               _d 0, 25.               _d 0,
336       &      4.03692566635316 _d 0, 22.84661726775120 _d 0,       &      4.03692566635316 _d 0, 22.84661726775120 _d 0,
337       &      3.62720389416752 _d 0, 22.62420229124846 _d 0/       &      3.62720389416752 _d 0, 22.62420229124846 _d 0/
338       &     sloc       &     sLoc
339       &     /35.5 _d 0, 35. _d 0,       &     /35.5 _d 0, 35. _d 0,
340         &      35.0 _d 0, 20. _d 0,
341         &      40.0 _d 0,
342       &       0.  _d 0,  0. _d 0,       &       0.  _d 0,  0. _d 0,
343       &      35.  _d 0, 35. _d 0,       &      35.  _d 0, 35. _d 0,
344       &       0.  _d 0,  0. _d 0,       &       0.  _d 0,  0. _d 0,
345       &      35.  _d 0, 35. _d 0/       &      35.  _d 0, 35. _d 0/
346       &     ploc       &     pLoc
347       &     /300. _d 5,  200. _d 5,       &     /300. _d 5,  200. _d 5,
348         &      200. _d 5,  100. _d 5,
349         &      800. _d 5,
350       &        0. _d 0,    0. _d 0,       &        0. _d 0,    0. _d 0,
351       &        0. _d 0,    0. _d 0,       &        0. _d 0,    0. _d 0,
352       &     1000. _d 5, 1000. _d 5,       &     1000. _d 5, 1000. _d 5,
353       &     1000. _d 5, 1000. _d 5/       &     1000. _d 5, 1000. _d 5/
354        DATA rloc        DATA rLoc
355       &     /1041.83267 _d 0, 1033.213387 _d 0,       &     /1041.83267  _d 0, 1033.213387 _d 0,
356       &       999.96675 _d 0,  997.04796  _d 0,       &      1031.654229 _d 0, 1017.726743 _d 0,
357       &      1027.67547 _d 0, 1023.34306  _d 0,       &      1062.928258 _d 0,
358       &      1044.12802 _d 0, 1037.90204  _d 0,       &       999.96675  _d 0,  997.04796  _d 0,
359       &      1069.48914 _d 0, 1062.53817  _d 0/       &      1027.67547  _d 0, 1023.34306  _d 0,
360       &     bloc       &      1044.12802  _d 0, 1037.90204  _d 0,
361         &      1069.48914  _d 0, 1062.53817  _d 0/
362         &     bLoc
363       &     /   -1.00000 _d 0,    -1.00000 _d 0,       &     /   -1.00000 _d 0,    -1.00000 _d 0,
364         &         -1.00000 _d 0,    -1.00000 _d 0,
365         &         -1.00000 _d 0,
366       &      20337.80375 _d 0, 22100.72106 _d 0,       &      20337.80375 _d 0, 22100.72106 _d 0,
367       &      22185.93358 _d 0, 23726.34949 _d 0,       &      22185.93358 _d 0, 23726.34949 _d 0,
368       &      23643.52599 _d 0, 25405.09717 _d 0,       &      23643.52599 _d 0, 25405.09717 _d 0,
369       &      25577.49819 _d 0, 27108.94504 _d 0/       &      25577.49819 _d 0, 27108.94504 _d 0/
370    
371          
372        bi   = 1        bi   = 1
373        bj   = 1        bj   = 1
374        k    = 1        k    = 1
375        imin = 1        iMin = 1
376        imax = 1        iMax = 1
377        jmin = 1        jMin = 1
378        jmax = 1        jMax = 1
379        i    = 1        i    = 1
380        j    = 1        j    = 1
381        if ( equationOfState.ne.'LINEAR'        IF (       equationOfState.NE.'LINEAR'
382       &     .and. equationOfState.ne.'POLY3' ) then       &     .AND. equationOfState.NE.'POLY3'
383         &     .AND. equationOfState.NE.'IDEALG' ) THEN
384  C     check nonlinear EOS  C     check nonlinear EOS
385           write(msgBuf,'(a,a)')          WRITE(msgBuf,'(A,A)')
386       &        'ml-eos: Check the equation of state: Type ',       &        'EOS_CHECK: Check the equation of state: Type ',
387       &        equationOfState       &        equationOfState
388           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
389       &        SQUEEZE_RIGHT , 1)       &       SQUEEZE_RIGHT , myThid )
390           psave = pressure(i,j,k,bi,bj)  C-    to get the right pressure out of S/R PRESSURE_FOR_EOS:
391           do kcheck = 1,ncheck          mskSave = maskC(i,j,k,bi,bj)
392              pressure(i,j,k,bi,bj) = ploc(kcheck)          rC_Save = rC(k)
393              if ( equationOfState.ne.'UNESCO' ) then          maskC(i,j,k,bi,bj) = 1.
394                 tFld(i,j,k,bi,bj) = ptloc(kcheck)          totPhiHyd(i,j,k,bi,bj) = 0.
395              else  C-    + set rC accordingly
396                 tFld(i,j,k,bi,bj) = tloc(kcheck)          DO kcheck = 1,ncheck
397              endif            IF ( usingZCoords ) THEN
398              sFld(i,j,k,bi,bj)    = sloc(kcheck)              rC(k) = -pLoc(kcheck)*recip_rhoConst*recip_gravity
399              rholoc               =  0. _d 0            ELSE
400              bulkMod(i,j)         = -1. _d 0              rC(k) =  pLoc(kcheck)
401                      ENDIF
402              call find_rhop0(            IF ( equationOfState.NE.'UNESCO' ) THEN
403       I           bi, bj, imin, imax, jmin, jmax, k,               tFld(i,j) = ptLoc(kcheck)
404       I           tFld, sFld,            ELSE
405       O           rhoP0,               tFld(i,j) = tLoc(kcheck)
406       I           myThid )            ENDIF
407                          sFld(i,j)    = sLoc(kcheck)
408              call find_bulkmod(            pFld(i,j)    = pLoc(kcheck)
409       I           bi, bj, imin, imax, jmin, jmax, k, k,            rhoLoc(i,j)  =  0. _d 0
410    
411              CALL FIND_RHO_2D(
412         I           iMin, iMax, jMin, jMax, k,
413       I           tFld, sFld,       I           tFld, sFld,
414         O           rhoLoc,
415         I           k, bi, bj, myThid )
416    
417              IF ( equationOfState.EQ.'MDJWF' ) THEN
418                bulkMod(i,j) = -1. _d 0
419              ELSE
420                CALL FIND_BULKMOD(
421         I           iMin, iMax, jMin, jMax,
422         I           pFld, tFld, sFld,
423       O           bulkMod,       O           bulkMod,
424       I           myThid )       I           myThid )
425                          ENDIF
426              rholoc = rhoP0(i,j)  
427       &           /(1. _d 0 -            IF ( kcheck .EQ. 1 ) THEN
428       &           pressure(i,j,k,bi,bj)*SItoBar/bulkMod(i,j))             WRITE(msgBuf,'(A)')
429                     &          'EOS_CHECK: check values for eosType=JMD95:'
430                     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
431              write(msgBuf,       &          SQUEEZE_RIGHT , myThid )
432       &           '(a4,f4.1,a5,f4.1,a6,f5.0,a5,a3,f10.5,1x,f11.5)')            ELSEIF ( kcheck .EQ. 2 ) THEN
433       &           'rho(', sFld(i,j,k,bi,bj), ' PSU,',             WRITE(msgBuf,'(A)')
434       &           tFld(i,j,k,bi,bj), ' degC,',       &          'EOS_CHECK: check values for eosType=MDJWF:'
435       &           pressure(i,j,k,bi,bj)*SItoBar, ' bar)',' = ',             CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
436       &           rloc(kcheck), bloc(kcheck)       &          SQUEEZE_RIGHT , myThid )
437              CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,            ELSEIF ( kcheck .EQ. 6 ) THEN
438       &           SQUEEZE_RIGHT , 1)             WRITE(msgBuf,'(A)')
439              write(msgBuf,'(a4,a32,f10.5,1x,f11.5)')       &          'EOS_CHECK: check values for eosType=UNESCO:'
440       &           'rho ', ' = ', rholoc, bulkMod(i,j)             CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
441              CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,       &          SQUEEZE_RIGHT , myThid )
442       &           SQUEEZE_RIGHT , 1)            ENDIF
443                          WRITE(msgBuf,'(2(A,F4.1),A,F5.0,2A,F10.5,1X,F11.5)')
444           enddo       &         'rho(', sFld(i,j), ' PSU,',
445         &         tFld(i,j), ' degC,',
446         &         pLoc(kcheck)*SItoBar, ' bar)',
447         &         ' = ', rLoc(kcheck), bLoc(kcheck)
448              CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
449         &         SQUEEZE_RIGHT , myThid )
450              WRITE(msgBuf,'(A,F10.5,1X,F11.5)')
451         &         ' rho(find_rho_2d)                 = ',
452         &         rhoLoc(i,j)+rhoConst, bulkMod(i,j)
453              CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
454         &         SQUEEZE_RIGHT , myThid )
455    
456              CALL FIND_RHO_SCALAR( tFld(i,j), sLoc(kcheck),
457         &         pLoc(kcheck), rhoLoc(i,j), myThid )
458              WRITE(msgBuf,'(A,F10.5,1X,F11.5)')
459         &         ' rho(find_rho_scalar)             = ',
460         &         rhoLoc(i,j)+rhoConst
461    c    &         rhoLoc(i,j)+rhoConst, bulkMod(i,j)
462              CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
463         &         SQUEEZE_RIGHT , myThid )
464    
465            ENDDO
466  C     end check nonlinear EOS  C     end check nonlinear EOS
467           pressure(i,j,k,bi,bj) = psave          maskC(i,j,k,bi,bj) = mskSave
468            rC(k) = rC_Save
469    
470           write(msgBuf,'(A)') 'end check the equation of state'          WRITE(msgBuf,'(A)') 'EOS_CHECK: Done'
471           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
472       &        SQUEEZE_RIGHT , 1)       &       SQUEEZE_RIGHT , myThid )
473    
474        endif        ENDIF
475    #endif /* INCLUDE_EOS_CHECK */
476    
477        return        RETURN
478        end        END
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.19
  ViewVC Help
Powered by ViewVC 1.1.22