/[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.7 by mlosch, Tue Oct 29 20:16:28 2002 UTC revision 1.13 by jmc, Fri Nov 4 01:19:24 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. 'ATMOSPHERIC'  
      &      .or. 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  
       else  
        STOP 'INI_EOS: We should never reach this point!'  
       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           OPEN(37,FILE='POLY3.COEFFS',STATUS='OLD',FORM='FORMATTED')           OPEN(37,FILE='POLY3.COEFFS',STATUS='OLD',FORM='FORMATTED')
71           READ(37,*) I           READ(37,*) I
# Line 144  C     rho = R(salinity, potential temper Line 92  C     rho = R(salinity, potential temper
92  C     pressure needs to be available (from the previous  C     pressure needs to be available (from the previous
93  C     time step to linearize the problem)  C     time step to linearize the problem)
94  C  C
95           if ( equationOfState .eq. 'JMD95Z'           if ( equationOfState .eq. 'JMD95Z' .and. usingPCoords ) then
      &        .and. buoyancyRelation .eq. 'OCEANICP' ) then  
96              write(msgBuf,'(A)')              write(msgBuf,'(A)')
97       &      'ini_eos: equation of state ''JMD95Z'' should not'       &      'ini_eos: equation of state ''JMD95Z'' should not'
98              CALL PRINT_ERROR( msgBuf , 1)              CALL PRINT_ERROR( msgBuf , 1)
# Line 160  C Line 107  C
107    
108  C     coefficients nonlinear equation of state in pressure coordinates for  C     coefficients nonlinear equation of state in pressure coordinates for
109  C     1. density of fresh water at p = 0  C     1. density of fresh water at p = 0
110           eosJMDCFw(1) =  999.842594           eosJMDCFw(1) =  999.842594 _d +00
111           eosJMDCFw(2) =    6.793952 _d -02           eosJMDCFw(2) =    6.793952 _d -02
112           eosJMDCFw(3) = -  9.095290 _d -03           eosJMDCFw(3) = -  9.095290 _d -03
113           eosJMDCFw(4) =    1.001685 _d -04           eosJMDCFw(4) =    1.001685 _d -04
# Line 169  C     1. density of fresh water at p = 0 Line 116  C     1. density of fresh water at p = 0
116  C     2. density of sea water at p = 0  C     2. density of sea water at p = 0
117           eosJMDCSw(1) =    8.24493  _d -01           eosJMDCSw(1) =    8.24493  _d -01
118           eosJMDCSw(2) = -  4.0899   _d -03           eosJMDCSw(2) = -  4.0899   _d -03
119           eosJMDCSw(3) =    7.6438   _d -05           eosJMDCSw(3) =    7.6438   _d -05
120           eosJMDCSw(4) = -  8.2467   _d -07           eosJMDCSw(4) = -  8.2467   _d -07
121           eosJMDCSw(5) =    5.3875   _d -09           eosJMDCSw(5) =    5.3875   _d -09
122           eosJMDCSw(6) = -  5.72466  _d -03           eosJMDCSw(6) = -  5.72466  _d -03
123           eosJMDCSw(7) =    1.0227   _d -04           eosJMDCSw(7) =    1.0227   _d -04
124           eosJMDCSw(8) = -  1.6546   _d -06           eosJMDCSw(8) = -  1.6546   _d -06
125           eosJMDCSw(9) =    4.8314   _d -04           eosJMDCSw(9) =    4.8314   _d -04
126           if ( equationOfState(1:5) .eq. 'JMD95' ) then           if ( equationOfState(1:5) .eq. 'JMD95' ) then
# Line 194  C     4. secant bulk modulus K of sea wa Line 141  C     4. secant bulk modulus K of sea wa
141  C     5. secant bulk modulus K of sea water at p  C     5. secant bulk modulus K of sea water at p
142              eosJMDCKP( 1) =   3.186519 _d +00              eosJMDCKP( 1) =   3.186519 _d +00
143              eosJMDCKP( 2) =   2.212276 _d -02              eosJMDCKP( 2) =   2.212276 _d -02
144              eosJMDCKP( 3) = - 2.984642 _d -04              eosJMDCKP( 3) = - 2.984642 _d -04
145              eosJMDCKP( 4) =   1.956415 _d -06              eosJMDCKP( 4) =   1.956415 _d -06
146              eosJMDCKP( 5) =   6.704388 _d -03              eosJMDCKP( 5) =   6.704388 _d -03
147              eosJMDCKP( 6) = - 1.847318 _d -04              eosJMDCKP( 6) = - 1.847318 _d -04
# Line 239  C     3. secant bulk modulus K of fresh Line 186  C     3. secant bulk modulus K of fresh
186              eosJMDCKFw(5) = - 5.155288 _d -05              eosJMDCKFw(5) = - 5.155288 _d -05
187  C     4. secant bulk modulus K of sea water at p = 0  C     4. secant bulk modulus K of sea water at p = 0
188              eosJMDCKSw(1) =   5.46746  _d +01              eosJMDCKSw(1) =   5.46746  _d +01
189              eosJMDCKSw(2) = - 0.603459 _d +00              eosJMDCKSw(2) = - 0.603459 _d +00
190              eosJMDCKSw(3) =   1.09987  _d -02              eosJMDCKSw(3) =   1.09987  _d -02
191              eosJMDCKSw(4) = - 6.1670   _d -05              eosJMDCKSw(4) = - 6.1670   _d -05
192              eosJMDCKSw(5) =   7.944    _d -02              eosJMDCKSw(5) =   7.944    _d -02
# Line 257  C     5. secant bulk modulus K of sea wa Line 204  C     5. secant bulk modulus K of sea wa
204              eosJMDCKP( 9) =   8.50935  _d -05              eosJMDCKP( 9) =   8.50935  _d -05
205              eosJMDCKP(10) = - 6.12293  _d -06              eosJMDCKP(10) = - 6.12293  _d -06
206              eosJMDCKP(11) =   5.2787   _d -08              eosJMDCKP(11) =   5.2787   _d -08
207              eosJMDCKP(12) = - 9.9348   _d -07              eosJMDCKP(12) = - 9.9348   _d -07
208              eosJMDCKP(13) =   2.0816   _d -08              eosJMDCKP(13) =   2.0816   _d -08
209              eosJMDCKP(14) =   9.1697   _d -10              eosJMDCKP(14) =   9.1697   _d -10
210           else           else
# Line 280  C     5. secant bulk modulus K of sea wa Line 227  C     5. secant bulk modulus K of sea wa
227           eosMDJWFnum(11) = -1.23869360 _d -11           eosMDJWFnum(11) = -1.23869360 _d -11
228                    
229                    
230           eosMDJWFden( 0) =  1.00000000 _d +00           eosMDJWFden( 0) =  1.00000000 _d +00
231           eosMDJWFden( 1) =  7.28606739 _d -03           eosMDJWFden( 1) =  7.28606739 _d -03
232           eosMDJWFden( 2) = -4.60835542 _d -05           eosMDJWFden( 2) = -4.60835542 _d -05
233           eosMDJWFden( 3) =  3.68390573 _d -07           eosMDJWFden( 3) =  3.68390573 _d -07
234           eosMDJWFden( 4) =  1.80809186 _d -10           eosMDJWFden( 4) =  1.80809186 _d -10
235           eosMDJWFden( 5) =  2.14691708 _d -03           eosMDJWFden( 5) =  2.14691708 _d -03
# Line 305  C Line 252  C
252                    
253        end if        end if
254    
255          _BEGIN_MASTER( myThid )
256    C--   Check EOS initialisation:
257    
258        call check_eos( myThid )        call check_eos( myThid )
259    
260        return        _END_MASTER( myThid )
261        end  
262          RETURN
263          END
264    
265  CBOP  CBOP
266  C     !ROUTINE: CHECK_EOS  C     !ROUTINE: CHECK_EOS
# Line 334  C     == Routine arguments == Line 286  C     == Routine arguments ==
286  C     myThid -  Number of this instance of CHECK_EOS  C     myThid -  Number of this instance of CHECK_EOS
287        INTEGER myThid        INTEGER myThid
288    
289    #ifndef EXCLUDE_EOS_CHECK
290  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
291  C     == Local variables ==  C     == Local variables ==
292  C     bi,bj  - Loop counters  C     bi,bj  - Loop counters
# Line 445  C     check nonlinear EOS Line 398  C     check nonlinear EOS
398       &           '(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)')
399       &           'rho(', sFld(i,j,k,bi,bj), ' PSU,',       &           'rho(', sFld(i,j,k,bi,bj), ' PSU,',
400       &           tFld(i,j,k,bi,bj), ' degC,',       &           tFld(i,j,k,bi,bj), ' degC,',
401       &           pressure(i,j,k,bi,bj)*SItoBar, ' bar)',' = ',  c    &           pressure(i,j,k,bi,bj)*SItoBar, ' bar)',' = ',
402       &           rloc(kcheck), bloc(kcheck)       &           rloc(kcheck), bloc(kcheck)
403              CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,              CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
404       &           SQUEEZE_RIGHT , 1)       &           SQUEEZE_RIGHT , 1)
# Line 466  C     check nonlinear EOS Line 419  C     check nonlinear EOS
419                            
420           enddo           enddo
421  C     end check nonlinear EOS  C     end check nonlinear EOS
422           pressure(i,j,k,bi,bj) = psave  c        pressure(i,j,k,bi,bj) = psave
423    
424           write(msgBuf,'(A)') 'end check the equation of state'           write(msgBuf,'(A)') 'end check the equation of state'
425           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
426       &        SQUEEZE_RIGHT , 1)       &        SQUEEZE_RIGHT , 1)
427    
428        endif        endif
429    #endif /* EXCLUDE_EOS_CHECK */
430    
431        return        RETURN
432        end        END
433    
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.13
  ViewVC Help
Powered by ViewVC 1.1.22