/[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.9 by jmc, Mon Feb 3 21:58:09 2003 UTC revision 1.14 by cnh, Mon Nov 7 18:26:02 2005 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    #define EXCLUDE_EOS_CHECK
6    
7  CBOP  CBOP
8  C     !ROUTINE: INI_EOS  C     !ROUTINE: INI_EOS
# Line 32  C     myThid -  Number of this instance Line 33  C     myThid -  Number of this instance
33    
34  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
35  C     == Local variables ==  C     == Local variables ==
36  C     bi,bj  - Loop counters  C     I,K    - Loop counters
37  C     I,J,K        INTEGER  I,  K
       INTEGER bi, bj  
       INTEGER  I,  J,  K  
38        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
39                
40          IF ( .NOT.fluidIsWater ) RETURN
41    
42        equationOfState = eosType        equationOfState = eosType
43    
# Line 63  C     I,J,K Line 63  C     I,J,K
63           eosMDJWFden(k) = 0. _d 0           eosMDJWFden(k) = 0. _d 0
64        end do        end do
65    
 C     initialise pressure to zero (this is not really the right place to  
 C     do it, but let's do it anyway; pressure will be initialised to  
 C     sensible values in ini_pressure)  
       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) = 0. _d 0  
           end do  
          end do  
         end do  
        end do  
       end do  
 C  
 C     initialise pressure to something sensible (will be overwritten)  
 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) = rhoConst * (  
      &           - gravity*rC(k)  
      &           )  
            end do  
           end do  
          end do  
         end do  
        end do  
       elseif ( buoyancyRelation .eq. 'OCEANICP' ) then  
 C      
 C     in pressure coordinates the pressure is just the coordinate of  
 C     the tracer point, that is, its constant in time  
 C  
        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  
   
66        if ( equationOfState .eq. 'LINEAR' ) then        if ( equationOfState .eq. 'LINEAR' ) then
67           if ( tAlpha .eq. UNSET_RL ) tAlpha = 2.D-4           if ( tAlpha .eq. UNSET_RL ) tAlpha = 2.  _d -4
68           if ( sBeta  .eq. UNSET_RL ) sBeta  = 7.4D-4           if ( sBeta  .eq. UNSET_RL ) sBeta  = 7.4 _d -4
69        elseif ( equationOfState .eq. 'POLY3' ) then        elseif ( equationOfState .eq. 'POLY3' ) then
70             _BEGIN_MASTER( myThid )
71           OPEN(37,FILE='POLY3.COEFFS',STATUS='OLD',FORM='FORMATTED')           OPEN(37,FILE='POLY3.COEFFS',STATUS='OLD',FORM='FORMATTED')
72           READ(37,*) I           READ(37,*) I
73           IF (I.NE.Nr) THEN           IF (I.NE.Nr) THEN
# Line 132  C Line 84  C
84              READ(37,*) (eosC(I,K),I=1,9)              READ(37,*) (eosC(I,K),I=1,9)
85           ENDDO           ENDDO
86           CLOSE(37)           CLOSE(37)
87             _END_MASTER( myThid )
88             _BARRIER
89    
90        elseif ( equationOfState(1:5) .eq. 'JMD95'        elseif ( equationOfState(1:5) .eq. 'JMD95'
91       &         .or. equationOfState .eq. 'UNESCO' ) then       &         .or. equationOfState .eq. 'UNESCO' ) then
# Line 141  C     rho = R(salinity, potential temper Line 95  C     rho = R(salinity, potential temper
95  C     pressure needs to be available (from the previous  C     pressure needs to be available (from the previous
96  C     time step to linearize the problem)  C     time step to linearize the problem)
97  C  C
98           if ( equationOfState .eq. 'JMD95Z'           if ( equationOfState .eq. 'JMD95Z' .and. usingPCoords ) then
      &        .and. buoyancyRelation .eq. 'OCEANICP' ) then  
99              write(msgBuf,'(A)')              write(msgBuf,'(A)')
100       &      'ini_eos: equation of state ''JMD95Z'' should not'       &      'ini_eos: equation of state ''JMD95Z'' should not'
101              CALL PRINT_ERROR( msgBuf , 1)              CALL PRINT_ERROR( msgBuf , 1)
# Line 157  C Line 110  C
110    
111  C     coefficients nonlinear equation of state in pressure coordinates for  C     coefficients nonlinear equation of state in pressure coordinates for
112  C     1. density of fresh water at p = 0  C     1. density of fresh water at p = 0
113           eosJMDCFw(1) =  999.842594           eosJMDCFw(1) =  999.842594 _d +00
114           eosJMDCFw(2) =    6.793952 _d -02           eosJMDCFw(2) =    6.793952 _d -02
115           eosJMDCFw(3) = -  9.095290 _d -03           eosJMDCFw(3) = -  9.095290 _d -03
116           eosJMDCFw(4) =    1.001685 _d -04           eosJMDCFw(4) =    1.001685 _d -04
# Line 166  C     1. density of fresh water at p = 0 Line 119  C     1. density of fresh water at p = 0
119  C     2. density of sea water at p = 0  C     2. density of sea water at p = 0
120           eosJMDCSw(1) =    8.24493  _d -01           eosJMDCSw(1) =    8.24493  _d -01
121           eosJMDCSw(2) = -  4.0899   _d -03           eosJMDCSw(2) = -  4.0899   _d -03
122           eosJMDCSw(3) =    7.6438   _d -05           eosJMDCSw(3) =    7.6438   _d -05
123           eosJMDCSw(4) = -  8.2467   _d -07           eosJMDCSw(4) = -  8.2467   _d -07
124           eosJMDCSw(5) =    5.3875   _d -09           eosJMDCSw(5) =    5.3875   _d -09
125           eosJMDCSw(6) = -  5.72466  _d -03           eosJMDCSw(6) = -  5.72466  _d -03
126           eosJMDCSw(7) =    1.0227   _d -04           eosJMDCSw(7) =    1.0227   _d -04
127           eosJMDCSw(8) = -  1.6546   _d -06           eosJMDCSw(8) = -  1.6546   _d -06
128           eosJMDCSw(9) =    4.8314   _d -04           eosJMDCSw(9) =    4.8314   _d -04
129           if ( equationOfState(1:5) .eq. 'JMD95' ) then           if ( equationOfState(1:5) .eq. 'JMD95' ) then
# Line 191  C     4. secant bulk modulus K of sea wa Line 144  C     4. secant bulk modulus K of sea wa
144  C     5. secant bulk modulus K of sea water at p  C     5. secant bulk modulus K of sea water at p
145              eosJMDCKP( 1) =   3.186519 _d +00              eosJMDCKP( 1) =   3.186519 _d +00
146              eosJMDCKP( 2) =   2.212276 _d -02              eosJMDCKP( 2) =   2.212276 _d -02
147              eosJMDCKP( 3) = - 2.984642 _d -04              eosJMDCKP( 3) = - 2.984642 _d -04
148              eosJMDCKP( 4) =   1.956415 _d -06              eosJMDCKP( 4) =   1.956415 _d -06
149              eosJMDCKP( 5) =   6.704388 _d -03              eosJMDCKP( 5) =   6.704388 _d -03
150              eosJMDCKP( 6) = - 1.847318 _d -04              eosJMDCKP( 6) = - 1.847318 _d -04
# Line 236  C     3. secant bulk modulus K of fresh Line 189  C     3. secant bulk modulus K of fresh
189              eosJMDCKFw(5) = - 5.155288 _d -05              eosJMDCKFw(5) = - 5.155288 _d -05
190  C     4. secant bulk modulus K of sea water at p = 0  C     4. secant bulk modulus K of sea water at p = 0
191              eosJMDCKSw(1) =   5.46746  _d +01              eosJMDCKSw(1) =   5.46746  _d +01
192              eosJMDCKSw(2) = - 0.603459 _d +00              eosJMDCKSw(2) = - 0.603459 _d +00
193              eosJMDCKSw(3) =   1.09987  _d -02              eosJMDCKSw(3) =   1.09987  _d -02
194              eosJMDCKSw(4) = - 6.1670   _d -05              eosJMDCKSw(4) = - 6.1670   _d -05
195              eosJMDCKSw(5) =   7.944    _d -02              eosJMDCKSw(5) =   7.944    _d -02
# Line 254  C     5. secant bulk modulus K of sea wa Line 207  C     5. secant bulk modulus K of sea wa
207              eosJMDCKP( 9) =   8.50935  _d -05              eosJMDCKP( 9) =   8.50935  _d -05
208              eosJMDCKP(10) = - 6.12293  _d -06              eosJMDCKP(10) = - 6.12293  _d -06
209              eosJMDCKP(11) =   5.2787   _d -08              eosJMDCKP(11) =   5.2787   _d -08
210              eosJMDCKP(12) = - 9.9348   _d -07              eosJMDCKP(12) = - 9.9348   _d -07
211              eosJMDCKP(13) =   2.0816   _d -08              eosJMDCKP(13) =   2.0816   _d -08
212              eosJMDCKP(14) =   9.1697   _d -10              eosJMDCKP(14) =   9.1697   _d -10
213           else           else
# Line 277  C     5. secant bulk modulus K of sea wa Line 230  C     5. secant bulk modulus K of sea wa
230           eosMDJWFnum(11) = -1.23869360 _d -11           eosMDJWFnum(11) = -1.23869360 _d -11
231                    
232                    
233           eosMDJWFden( 0) =  1.00000000 _d +00           eosMDJWFden( 0) =  1.00000000 _d +00
234           eosMDJWFden( 1) =  7.28606739 _d -03           eosMDJWFden( 1) =  7.28606739 _d -03
235           eosMDJWFden( 2) = -4.60835542 _d -05           eosMDJWFden( 2) = -4.60835542 _d -05
236           eosMDJWFden( 3) =  3.68390573 _d -07           eosMDJWFden( 3) =  3.68390573 _d -07
237           eosMDJWFden( 4) =  1.80809186 _d -10           eosMDJWFden( 4) =  1.80809186 _d -10
238           eosMDJWFden( 5) =  2.14691708 _d -03           eosMDJWFden( 5) =  2.14691708 _d -03
# Line 307  C--   Check EOS initialisation: Line 260  C--   Check EOS initialisation:
260    
261        call check_eos( myThid )        call check_eos( myThid )
262    
       IF ( buoyancyRelation .EQ. 'OCEANIC' .AND.  
      &     ( equationOfState .EQ. 'JMD95P'.OR.  
      &       equationOfState .EQ. 'MDJWF' .OR.  
      &       equationOfState .EQ. 'UNESCO'   ) ) THEN  
         WRITE(msgBuf,'(A)')  
      &   'S/R INI_EOS: WARNING >>> Wrong Restart !!!'  
         CALL PRINT_MESSAGE(msgBuf, errorMessageUnit,  
      &                     SQUEEZE_RIGHT,myThid)  
         WRITE(msgBuf,'(2A)') 'S/R INI_EOS: ',  
      &   'not correct using Z-coord with EOS function of P'  
         CALL PRINT_MESSAGE(msgBuf, errorMessageUnit,  
      &                     SQUEEZE_RIGHT,myThid)  
         WRITE(msgBuf,'(A)')  
      &   'S/R INI_EOS: WARNING <<< Wrong Restart !!!'  
         CALL PRINT_MESSAGE(msgBuf, errorMessageUnit,  
      &                     SQUEEZE_RIGHT,myThid)  
       ENDIF                                          
   
263        _END_MASTER( myThid )        _END_MASTER( myThid )
264    
265        return        RETURN
266        end        END
267    
268  CBOP  CBOP
269  C     !ROUTINE: CHECK_EOS  C     !ROUTINE: CHECK_EOS
# Line 354  C     == Routine arguments == Line 289  C     == Routine arguments ==
289  C     myThid -  Number of this instance of CHECK_EOS  C     myThid -  Number of this instance of CHECK_EOS
290        INTEGER myThid        INTEGER myThid
291    
292    #ifndef EXCLUDE_EOS_CHECK
293  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
294  C     == Local variables ==  C     == Local variables ==
295  C     bi,bj  - Loop counters  C     bi,bj  - Loop counters
# Line 465  C     check nonlinear EOS Line 401  C     check nonlinear EOS
401       &           '(a4,f4.1,a5,f4.1,a6,f5.0,a5,a3,f10.5,1x,f11.5)')       &           '(a4,f4.1,a5,f4.1,a6,f5.0,a5,a3,f10.5,1x,f11.5)')
402       &           'rho(', sFld(i,j,k,bi,bj), ' PSU,',       &           'rho(', sFld(i,j,k,bi,bj), ' PSU,',
403       &           tFld(i,j,k,bi,bj), ' degC,',       &           tFld(i,j,k,bi,bj), ' degC,',
404       &           pressure(i,j,k,bi,bj)*SItoBar, ' bar)',' = ',  c    &           pressure(i,j,k,bi,bj)*SItoBar, ' bar)',' = ',
405       &           rloc(kcheck), bloc(kcheck)       &           rloc(kcheck), bloc(kcheck)
406              CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,              CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
407       &           SQUEEZE_RIGHT , 1)       &           SQUEEZE_RIGHT , 1)
# Line 486  C     check nonlinear EOS Line 422  C     check nonlinear EOS
422                            
423           enddo           enddo
424  C     end check nonlinear EOS  C     end check nonlinear EOS
425           pressure(i,j,k,bi,bj) = psave  c        pressure(i,j,k,bi,bj) = psave
426    
427           write(msgBuf,'(A)') 'end check the equation of state'           write(msgBuf,'(A)') 'end check the equation of state'
428           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
429       &        SQUEEZE_RIGHT , 1)       &        SQUEEZE_RIGHT , 1)
430    
431        endif        endif
432    #endif /* EXCLUDE_EOS_CHECK */
433    
434        return        RETURN
435        end        END
436    
Legend:
Removed from v.1.9  
changed lines
  Added in v.1.14
  ViewVC Help
Powered by ViewVC 1.1.22