/[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.4 by mlosch, Thu Sep 5 20:49:33 2002 UTC revision 1.8 by jmc, Wed Nov 6 03:45:46 2002 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. 'ATMOSPHERIC'
98         &      .or. buoyancyRelation .eq. 'OCEANICP' ) then
99    C    
100    C     in pressure coordinates the pressure is just the coordinate of
101    C     the tracer point, that is, its constant in time
102    C
103           do bj = myByLo(myThid), myByHi(myThid)
104            do bi = myBxLo(myThid), myBxHi(myThid)
105             do K=1,Nr
106              do J=1-Oly,sNy+Oly
107               do I=1-Olx,sNx+Olx
108                pressure(i,j,k,bi,bj) = rC(k)
109               end do
110              end do
111             end do
112            end do
113           end do
114          else
115           STOP 'INI_EOS: We should never reach this point!'
116          endif
117    
118        if ( equationOfState .eq. 'LINEAR' ) then        if ( equationOfState .eq. 'LINEAR' ) then
119           if ( tAlpha .eq. UNSET_RL ) tAlpha = 2.D-4           if ( tAlpha .eq. UNSET_RL ) tAlpha = 2.D-4
120           if ( sBeta  .eq. UNSET_RL ) sBeta  = 7.4D-4           if ( sBeta  .eq. UNSET_RL ) sBeta  = 7.4D-4
# Line 105  C Line 157  C
157              CALL PRINT_ERROR( msgBuf , 1)              CALL PRINT_ERROR( msgBuf , 1)
158              STOP 'ABNORMAL END: S/R INI_EOS'              STOP 'ABNORMAL END: S/R INI_EOS'
159           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  
160    
161  C     coefficients nonlinear equation of state in pressure coordinates for  C     coefficients nonlinear equation of state in pressure coordinates for
162  C     1. density of fresh water at p = 0  C     1. density of fresh water at p = 0
# Line 202  C     5. secant bulk modulus K of sea wa Line 223  C     5. secant bulk modulus K of sea wa
223              CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,              CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
224       &           SQUEEZE_RIGHT , 1)       &           SQUEEZE_RIGHT , 1)
225              write(msgBuf,'(a)')              write(msgBuf,'(a)')
226       &           'WARNING: (see Jackett and McDougall 1995, JPO)'       &           'WARNING: (see Jackett and McDougall 1995, JAOT)'
227              CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,              CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
228       &           SQUEEZE_RIGHT , 1)       &           SQUEEZE_RIGHT , 1)
229              write(msgBuf,'(a)')              write(msgBuf,'(a)')
# Line 245  C     5. secant bulk modulus K of sea wa Line 266  C     5. secant bulk modulus K of sea wa
266    
267        elseif ( equationOfState .eq. 'MDJWF' ) then        elseif ( equationOfState .eq. 'MDJWF' ) then
268    
269           eosMDJWFnum( 0) =  9.99843699e+02           eosMDJWFnum( 0) =  9.99843699 _d +02
270           eosMDJWFnum( 1) =  7.35212840e+00           eosMDJWFnum( 1) =  7.35212840 _d +00
271           eosMDJWFnum( 2) = -5.45928211e-02           eosMDJWFnum( 2) = -5.45928211 _d -02
272           eosMDJWFnum( 3) =  3.98476704e-04           eosMDJWFnum( 3) =  3.98476704 _d -04
273           eosMDJWFnum( 4) =  2.96938239e+00           eosMDJWFnum( 4) =  2.96938239 _d +00
274           eosMDJWFnum( 5) = -7.23268813e-03           eosMDJWFnum( 5) = -7.23268813 _d -03
275           eosMDJWFnum( 6) =  2.12382341e-03           eosMDJWFnum( 6) =  2.12382341 _d -03
276           eosMDJWFnum( 7) =  1.04004591e-02           eosMDJWFnum( 7) =  1.04004591 _d -02
277           eosMDJWFnum( 8) =  1.03970529e-07           eosMDJWFnum( 8) =  1.03970529 _d -07
278           eosMDJWFnum( 9) =  5.18761880e-06           eosMDJWFnum( 9) =  5.18761880 _d -06
279           eosMDJWFnum(10) = -3.24041825e-08           eosMDJWFnum(10) = -3.24041825 _d -08
280           eosMDJWFnum(11) = -1.23869360e-11           eosMDJWFnum(11) = -1.23869360 _d -11
281                    
282                    
283           eosMDJWFden( 0) =  1.00000000e+00           eosMDJWFden( 0) =  1.00000000 _d +00
284           eosMDJWFden( 1) =  7.28606739e-03           eosMDJWFden( 1) =  7.28606739 _d -03
285           eosMDJWFden( 2) = -4.60835542e-05           eosMDJWFden( 2) = -4.60835542 _d -05
286           eosMDJWFden( 3) =  3.68390573e-07           eosMDJWFden( 3) =  3.68390573 _d -07
287           eosMDJWFden( 4) =  1.80809186e-10           eosMDJWFden( 4) =  1.80809186 _d -10
288           eosMDJWFden( 5) =  2.14691708e-03           eosMDJWFden( 5) =  2.14691708 _d -03
289           eosMDJWFden( 6) = -9.27062484e-06           eosMDJWFden( 6) = -9.27062484 _d -06
290           eosMDJWFden( 7) = -1.78343643e-10           eosMDJWFden( 7) = -1.78343643 _d -10
291           eosMDJWFden( 8) =  4.76534122e-06           eosMDJWFden( 8) =  4.76534122 _d -06
292           eosMDJWFden( 9) =  1.63410736e-09           eosMDJWFden( 9) =  1.63410736 _d -09
293           eosMDJWFden(10) =  5.30848875e-06           eosMDJWFden(10) =  5.30848875 _d -06
294           eosMDJWFden(11) = -3.03175128e-16           eosMDJWFden(11) = -3.03175128 _d -16
295           eosMDJWFden(12) = -1.27934137e-17           eosMDJWFden(12) = -1.27934137 _d -17
296                    
297        elseif( equationOfState .eq. 'IDEALG' ) then        elseif( equationOfState .eq. 'IDEALG' ) then
298  C      C    
# Line 284  C Line 305  C
305                    
306        end if        end if
307    
308          _BEGIN_MASTER( myThid )
309        call check_eos( myThid )        call check_eos( myThid )
310          _END_MASTER( myThid )
311    
312        return        return
313        end        end
# Line 414  C     check nonlinear EOS Line 437  C     check nonlinear EOS
437                    
438              call find_rho(              call find_rho(
439       &           bi, bj, iMin, iMax, jMin, jMax,  k, k,       &           bi, bj, iMin, iMax, jMin, jMax,  k, k,
      &           equationOfState,  
440       &           tFld, sFld, rholoc, myThid )       &           tFld, sFld, rholoc, myThid )
441    
442              call find_bulkmod(              call find_bulkmod(
443       &           bi, bj, imin, imax, jmin, jmax, k, k,       &           bi, bj, imin, imax, jmin, jmax, k, k,
444       &           tFld, sFld, bulkMod, myThid )       &           tFld, sFld, bulkMod, myThid )
445                
         
446              write(msgBuf,              write(msgBuf,
447       &           '(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)')
448       &           'rho(', sFld(i,j,k,bi,bj), ' PSU,',       &           'rho(', sFld(i,j,k,bi,bj), ' PSU,',
# Line 430  C     check nonlinear EOS Line 451  C     check nonlinear EOS
451       &           rloc(kcheck), bloc(kcheck)       &           rloc(kcheck), bloc(kcheck)
452              CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,              CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
453       &           SQUEEZE_RIGHT , 1)       &           SQUEEZE_RIGHT , 1)
454              write(msgBuf,'(a4,a32,f10.5,1x,f11.5)')              write(msgBuf,'(a14,a22,f10.5,1x,f11.5)')
455       &           'rho ', ' = ', rholoc(i,j)+rhoNil, bulkMod(i,j)       &           'rho(find_rho) ',
456         &           ' = ', rholoc(i,j)+rhoConst, bulkMod(i,j)
457                CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
458         &           SQUEEZE_RIGHT , 1)
459    
460                call find_rho_scalar( tFld(i,j,k,bi,bj), sLoc(kcheck),
461         &           pLoc(kcheck), rhoLoc(i,j), myThid )
462                bulkMod(i,j) = 0. _d 0
463                write(msgBuf,'(a21,a15,f10.5,1x,f11.5)')
464         &           'rho(find_rho_scalar) ',
465         &           ' = ', rholoc(i,j)+rhoConst, bulkMod(i,j)
466              CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,              CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
467       &           SQUEEZE_RIGHT , 1)       &           SQUEEZE_RIGHT , 1)
468                            
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.8
  ViewVC Help
Powered by ViewVC 1.1.22