/[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.5 by mlosch, Wed Sep 18 16:38:02 2002 UTC revision 1.9 by jmc, Mon Feb 3 21:58:09 2003 UTC
# 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    
66    C     initialise pressure to zero (this is not really the right place to
67    C     do it, but let's do it anyway; pressure will be initialised to
68    C     sensible values in ini_pressure)
69          do bj = myByLo(myThid), myByHi(myThid)
70           do bi = myBxLo(myThid), myBxHi(myThid)
71            do K=1,Nr
72             do J=1-Oly,sNy+Oly
73              do I=1-Olx,sNx+Olx
74               pressure(i,j,k,bi,bj) = 0. _d 0
75              end do
76             end do
77            end do
78           end do
79          end do
80    C
81    C     initialise pressure to something sensible (will be overwritten)
82    C
83          if ( buoyancyRelation .eq. 'OCEANIC' ) then
84           do bj = myByLo(myThid), myByHi(myThid)
85            do bi = myBxLo(myThid), myBxHi(myThid)
86             do K=1,Nr
87              do J=1-Oly,sNy+Oly
88               do I=1-Olx,sNx+Olx
89                pressure(i,j,k,bi,bj) = rhoConst * (
90         &           - gravity*rC(k)
91         &           )
92               end do
93              end do
94             end do
95            end do
96           end do
97          elseif ( buoyancyRelation .eq. 'OCEANICP' ) then
98    C    
99    C     in pressure coordinates the pressure is just the coordinate of
100    C     the tracer point, that is, its constant in time
101    C
102           do bj = myByLo(myThid), myByHi(myThid)
103            do bi = myBxLo(myThid), myBxHi(myThid)
104             do K=1,Nr
105              do J=1-Oly,sNy+Oly
106               do I=1-Olx,sNx+Olx
107                pressure(i,j,k,bi,bj) = rC(k)
108               end do
109              end do
110             end do
111            end do
112           end do
113          endif
114    
115        if ( equationOfState .eq. 'LINEAR' ) then        if ( equationOfState .eq. 'LINEAR' ) then
116           if ( tAlpha .eq. UNSET_RL ) tAlpha = 2.D-4           if ( tAlpha .eq. UNSET_RL ) tAlpha = 2.D-4
117           if ( sBeta  .eq. UNSET_RL ) sBeta  = 7.4D-4           if ( sBeta  .eq. UNSET_RL ) sBeta  = 7.4D-4
# Line 105  C Line 154  C
154              CALL PRINT_ERROR( msgBuf , 1)              CALL PRINT_ERROR( msgBuf , 1)
155              STOP 'ABNORMAL END: S/R INI_EOS'              STOP 'ABNORMAL END: S/R INI_EOS'
156           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  
157    
158  C     coefficients nonlinear equation of state in pressure coordinates for  C     coefficients nonlinear equation of state in pressure coordinates for
159  C     1. density of fresh water at p = 0  C     1. density of fresh water at p = 0
# Line 245  C     5. secant bulk modulus K of sea wa Line 263  C     5. secant bulk modulus K of sea wa
263    
264        elseif ( equationOfState .eq. 'MDJWF' ) then        elseif ( equationOfState .eq. 'MDJWF' ) then
265    
266           eosMDJWFnum( 0) =  9.99843699e+02           eosMDJWFnum( 0) =  9.99843699 _d +02
267           eosMDJWFnum( 1) =  7.35212840e+00           eosMDJWFnum( 1) =  7.35212840 _d +00
268           eosMDJWFnum( 2) = -5.45928211e-02           eosMDJWFnum( 2) = -5.45928211 _d -02
269           eosMDJWFnum( 3) =  3.98476704e-04           eosMDJWFnum( 3) =  3.98476704 _d -04
270           eosMDJWFnum( 4) =  2.96938239e+00           eosMDJWFnum( 4) =  2.96938239 _d +00
271           eosMDJWFnum( 5) = -7.23268813e-03           eosMDJWFnum( 5) = -7.23268813 _d -03
272           eosMDJWFnum( 6) =  2.12382341e-03           eosMDJWFnum( 6) =  2.12382341 _d -03
273           eosMDJWFnum( 7) =  1.04004591e-02           eosMDJWFnum( 7) =  1.04004591 _d -02
274           eosMDJWFnum( 8) =  1.03970529e-07           eosMDJWFnum( 8) =  1.03970529 _d -07
275           eosMDJWFnum( 9) =  5.18761880e-06           eosMDJWFnum( 9) =  5.18761880 _d -06
276           eosMDJWFnum(10) = -3.24041825e-08           eosMDJWFnum(10) = -3.24041825 _d -08
277           eosMDJWFnum(11) = -1.23869360e-11           eosMDJWFnum(11) = -1.23869360 _d -11
278                    
279                    
280           eosMDJWFden( 0) =  1.00000000e+00           eosMDJWFden( 0) =  1.00000000 _d +00
281           eosMDJWFden( 1) =  7.28606739e-03           eosMDJWFden( 1) =  7.28606739 _d -03
282           eosMDJWFden( 2) = -4.60835542e-05           eosMDJWFden( 2) = -4.60835542 _d -05
283           eosMDJWFden( 3) =  3.68390573e-07           eosMDJWFden( 3) =  3.68390573 _d -07
284           eosMDJWFden( 4) =  1.80809186e-10           eosMDJWFden( 4) =  1.80809186 _d -10
285           eosMDJWFden( 5) =  2.14691708e-03           eosMDJWFden( 5) =  2.14691708 _d -03
286           eosMDJWFden( 6) = -9.27062484e-06           eosMDJWFden( 6) = -9.27062484 _d -06
287           eosMDJWFden( 7) = -1.78343643e-10           eosMDJWFden( 7) = -1.78343643 _d -10
288           eosMDJWFden( 8) =  4.76534122e-06           eosMDJWFden( 8) =  4.76534122 _d -06
289           eosMDJWFden( 9) =  1.63410736e-09           eosMDJWFden( 9) =  1.63410736 _d -09
290           eosMDJWFden(10) =  5.30848875e-06           eosMDJWFden(10) =  5.30848875 _d -06
291           eosMDJWFden(11) = -3.03175128e-16           eosMDJWFden(11) = -3.03175128 _d -16
292           eosMDJWFden(12) = -1.27934137e-17           eosMDJWFden(12) = -1.27934137 _d -17
293                    
294        elseif( equationOfState .eq. 'IDEALG' ) then        elseif( equationOfState .eq. 'IDEALG' ) then
295  C      C    
# Line 284  C Line 302  C
302                    
303        end if        end if
304    
305          _BEGIN_MASTER( myThid )
306    C--   Check EOS initialisation:
307    
308        call check_eos( myThid )        call check_eos( myThid )
309    
310          IF ( buoyancyRelation .EQ. 'OCEANIC' .AND.
311         &     ( equationOfState .EQ. 'JMD95P'.OR.
312         &       equationOfState .EQ. 'MDJWF' .OR.
313         &       equationOfState .EQ. 'UNESCO'   ) ) THEN
314            WRITE(msgBuf,'(A)')
315         &   'S/R INI_EOS: WARNING >>> Wrong Restart !!!'
316            CALL PRINT_MESSAGE(msgBuf, errorMessageUnit,
317         &                     SQUEEZE_RIGHT,myThid)
318            WRITE(msgBuf,'(2A)') 'S/R INI_EOS: ',
319         &   'not correct using Z-coord with EOS function of P'
320            CALL PRINT_MESSAGE(msgBuf, errorMessageUnit,
321         &                     SQUEEZE_RIGHT,myThid)
322            WRITE(msgBuf,'(A)')
323         &   'S/R INI_EOS: WARNING <<< Wrong Restart !!!'
324            CALL PRINT_MESSAGE(msgBuf, errorMessageUnit,
325         &                     SQUEEZE_RIGHT,myThid)
326          ENDIF                                        
327    
328          _END_MASTER( myThid )
329    
330        return        return
331        end        end
332    
# Line 430  C     check nonlinear EOS Line 471  C     check nonlinear EOS
471       &           SQUEEZE_RIGHT , 1)       &           SQUEEZE_RIGHT , 1)
472              write(msgBuf,'(a14,a22,f10.5,1x,f11.5)')              write(msgBuf,'(a14,a22,f10.5,1x,f11.5)')
473       &           'rho(find_rho) ',       &           'rho(find_rho) ',
474       &           ' = ', rholoc(i,j)+rhoNil, bulkMod(i,j)       &           ' = ', rholoc(i,j)+rhoConst, bulkMod(i,j)
475              CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,              CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
476       &           SQUEEZE_RIGHT , 1)       &           SQUEEZE_RIGHT , 1)
477    
# Line 439  C     check nonlinear EOS Line 480  C     check nonlinear EOS
480              bulkMod(i,j) = 0. _d 0              bulkMod(i,j) = 0. _d 0
481              write(msgBuf,'(a21,a15,f10.5,1x,f11.5)')              write(msgBuf,'(a21,a15,f10.5,1x,f11.5)')
482       &           'rho(find_rho_scalar) ',       &           'rho(find_rho_scalar) ',
483       &           ' = ', rholoc(i,j)+rhoNil, bulkMod(i,j)       &           ' = ', rholoc(i,j)+rhoConst, bulkMod(i,j)
484              CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,              CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
485       &           SQUEEZE_RIGHT , 1)       &           SQUEEZE_RIGHT , 1)
486                            
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.9
  ViewVC Help
Powered by ViewVC 1.1.22