/[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.7 by mlosch, Tue Oct 29 20:16:28 2002 UTC
# Line 23  C     == Global variables == Line 23  C     == Global variables ==
23  #include "PARAMS.h"  #include "PARAMS.h"
24  #include "EOS.h"  #include "EOS.h"
25  #include "GRID.h"  #include "GRID.h"
26    #include "DYNVARS.h"
27    
28  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
29  C     == Routine arguments ==  C     == Routine arguments ==
# Line 40  C     I,J,K Line 41  C     I,J,K
41    
42        equationOfState = eosType        equationOfState = eosType
43    
44          do k = 1,6
45             eosJMDCFw(k) = 0. _d 0
46          end do
47          do k = 1,9
48             eosJMDCSw(k) = 0. _d 0
49          end do
50          do k = 1,5
51             eosJMDCKFw(k) = 0. _d 0
52          end do
53          do k = 1,7
54             eosJMDCKSw(k) = 0. _d 0
55          end do
56          do k = 1,14
57             eosJMDCKP(k) = 0. _d 0
58          end do
59          do k = 0,11
60             eosMDJWFnum(k) = 0. _d 0
61          end do
62          do k = 0,12
63             eosMDJWFden(k) = 0. _d 0
64          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 82  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 179  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 220  C     5. secant bulk modulus K of sea wa Line 264  C     5. secant bulk modulus K of sea wa
264              STOP 'INI_EOS: We should never reach this point!'              STOP 'INI_EOS: We should never reach this point!'
265           endif           endif
266    
267          elseif ( equationOfState .eq. 'MDJWF' ) then
268    
269             eosMDJWFnum( 0) =  9.99843699 _d +02
270             eosMDJWFnum( 1) =  7.35212840 _d +00
271             eosMDJWFnum( 2) = -5.45928211 _d -02
272             eosMDJWFnum( 3) =  3.98476704 _d -04
273             eosMDJWFnum( 4) =  2.96938239 _d +00
274             eosMDJWFnum( 5) = -7.23268813 _d -03
275             eosMDJWFnum( 6) =  2.12382341 _d -03
276             eosMDJWFnum( 7) =  1.04004591 _d -02
277             eosMDJWFnum( 8) =  1.03970529 _d -07
278             eosMDJWFnum( 9) =  5.18761880 _d -06
279             eosMDJWFnum(10) = -3.24041825 _d -08
280             eosMDJWFnum(11) = -1.23869360 _d -11
281            
282            
283             eosMDJWFden( 0) =  1.00000000 _d +00
284             eosMDJWFden( 1) =  7.28606739 _d -03
285             eosMDJWFden( 2) = -4.60835542 _d -05
286             eosMDJWFden( 3) =  3.68390573 _d -07
287             eosMDJWFden( 4) =  1.80809186 _d -10
288             eosMDJWFden( 5) =  2.14691708 _d -03
289             eosMDJWFden( 6) = -9.27062484 _d -06
290             eosMDJWFden( 7) = -1.78343643 _d -10
291             eosMDJWFden( 8) =  4.76534122 _d -06
292             eosMDJWFden( 9) =  1.63410736 _d -09
293             eosMDJWFden(10) =  5.30848875 _d -06
294             eosMDJWFden(11) = -3.03175128 _d -16
295             eosMDJWFden(12) = -1.27934137 _d -17
296            
297        elseif( equationOfState .eq. 'IDEALG' ) then        elseif( equationOfState .eq. 'IDEALG' ) then
298  C      C    
299        else        else
# Line 269  C     I,J,K Line 343  C     I,J,K
343        INTEGER  I,  J,  K        INTEGER  I,  J,  K
344        _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
345        _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
346        _RL rhoP0  (1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL rhoLoc (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
347        _RL bulkMod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL bulkMod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
348        _RL rholoc, psave        _RL psave
349    
350        INTEGER ncheck, kcheck        INTEGER ncheck, kcheck
351        PARAMETER ( ncheck = 10 )        PARAMETER ( ncheck = 13 )
352        _RL tloc(ncheck), ptloc(ncheck), sloc(ncheck), ploc(ncheck)        _RL tloc(ncheck), ptloc(ncheck), sloc(ncheck), ploc(ncheck)
353        _RL rloc(ncheck), bloc(ncheck)        _RL rloc(ncheck), bloc(ncheck)
354    
# Line 282  C     I,J,K Line 356  C     I,J,K
356    
357        DATA tloc        DATA tloc
358       &     /3.25905152915860 _d 0, 20.38687090048638 _d 0,       &     /3.25905152915860 _d 0, 20.38687090048638 _d 0,
359         &     25.44820830309568 _d 0, 20.17368557065936 _d 0,
360         &     13.43397459640398 _d 0,
361       &      5.               _d 0, 25.               _d 0,       &      5.               _d 0, 25.               _d 0,
362       &      5.               _d 0, 25.               _d 0,       &      5.               _d 0, 25.               _d 0,
363       &      5.               _d 0, 25.               _d 0,       &      5.               _d 0, 25.               _d 0,
364       &      5.               _d 0, 25.               _d 0/,       &      5.               _d 0, 25.               _d 0/,
365       &     ptloc       &     ptloc
366       &     /3.               _d 0, 20.               _d 0,       &     /3.               _d 0, 20.               _d 0,
367         &     25.               _d 0, 20.               _d 0,
368         &     12.               _d 0,
369       &      5.               _d 0, 25.               _d 0,       &      5.               _d 0, 25.               _d 0,
370       &      5.               _d 0, 25.               _d 0,       &      5.               _d 0, 25.               _d 0,
371       &      4.03692566635316 _d 0, 22.84661726775120 _d 0,       &      4.03692566635316 _d 0, 22.84661726775120 _d 0,
372       &      3.62720389416752 _d 0, 22.62420229124846 _d 0/       &      3.62720389416752 _d 0, 22.62420229124846 _d 0/
373       &     sloc       &     sloc
374       &     /35.5 _d 0, 35. _d 0,       &     /35.5 _d 0, 35. _d 0,
375         &      35.0 _d 0, 20. _d 0,
376         &      40.0 _d 0,
377       &       0.  _d 0,  0. _d 0,       &       0.  _d 0,  0. _d 0,
378       &      35.  _d 0, 35. _d 0,       &      35.  _d 0, 35. _d 0,
379       &       0.  _d 0,  0. _d 0,       &       0.  _d 0,  0. _d 0,
380       &      35.  _d 0, 35. _d 0/       &      35.  _d 0, 35. _d 0/
381       &     ploc       &     ploc
382       &     /300. _d 5,  200. _d 5,       &     /300. _d 5,  200. _d 5,
383         &      200. _d 5,  100. _d 5,
384         &      800. _d 5,
385       &        0. _d 0,    0. _d 0,       &        0. _d 0,    0. _d 0,
386       &        0. _d 0,    0. _d 0,       &        0. _d 0,    0. _d 0,
387       &     1000. _d 5, 1000. _d 5,       &     1000. _d 5, 1000. _d 5,
388       &     1000. _d 5, 1000. _d 5/       &     1000. _d 5, 1000. _d 5/
389        DATA rloc        DATA rloc
390       &     /1041.83267 _d 0, 1033.213387 _d 0,       &     /1041.83267  _d 0, 1033.213387 _d 0,
391       &       999.96675 _d 0,  997.04796  _d 0,       &      1031.654229 _d 0, 1017.726743 _d 0,
392       &      1027.67547 _d 0, 1023.34306  _d 0,       &      1062.928258 _d 0,
393       &      1044.12802 _d 0, 1037.90204  _d 0,       &       999.96675  _d 0,  997.04796  _d 0,
394       &      1069.48914 _d 0, 1062.53817  _d 0/       &      1027.67547  _d 0, 1023.34306  _d 0,
395         &      1044.12802  _d 0, 1037.90204  _d 0,
396         &      1069.48914  _d 0, 1062.53817  _d 0/
397       &     bloc       &     bloc
398       &     /   -1.00000 _d 0,    -1.00000 _d 0,       &     /   -1.00000 _d 0,    -1.00000 _d 0,
399         &         -1.00000 _d 0,    -1.00000 _d 0,
400         &         -1.00000 _d 0,
401       &      20337.80375 _d 0, 22100.72106 _d 0,       &      20337.80375 _d 0, 22100.72106 _d 0,
402       &      22185.93358 _d 0, 23726.34949 _d 0,       &      22185.93358 _d 0, 23726.34949 _d 0,
403       &      23643.52599 _d 0, 25405.09717 _d 0,       &      23643.52599 _d 0, 25405.09717 _d 0,
# Line 331  C     I,J,K Line 417  C     I,J,K
417       &     .and. equationOfState.ne.'POLY3' ) then       &     .and. equationOfState.ne.'POLY3' ) then
418  C     check nonlinear EOS  C     check nonlinear EOS
419           write(msgBuf,'(a,a)')           write(msgBuf,'(a,a)')
420       &        'ml-eos: Check the equation of state: Type ',       &        'check_eos: Check the equation of state: Type ',
421       &        equationOfState       &        equationOfState
422           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
423       &        SQUEEZE_RIGHT , 1)       &        SQUEEZE_RIGHT , 1)
# Line 344  C     check nonlinear EOS Line 430  C     check nonlinear EOS
430                 tFld(i,j,k,bi,bj) = tloc(kcheck)                 tFld(i,j,k,bi,bj) = tloc(kcheck)
431              endif              endif
432              sFld(i,j,k,bi,bj)    = sloc(kcheck)              sFld(i,j,k,bi,bj)    = sloc(kcheck)
433              rholoc               =  0. _d 0              rholoc(i,j)          =  0. _d 0
434              bulkMod(i,j)         = -1. _d 0              bulkMod(i,j)         = -1. _d 0
435                    
436              call find_rhop0(              call find_rho(
437       I           bi, bj, imin, imax, jmin, jmax, k,       &           bi, bj, iMin, iMax, jMin, jMax,  k, k,
438       I           tFld, sFld,       &           tFld, sFld, rholoc, myThid )
439       O           rhoP0,  
      I           myThid )  
               
440              call find_bulkmod(              call find_bulkmod(
441       I           bi, bj, imin, imax, jmin, jmax, k, k,       &           bi, bj, imin, imax, jmin, jmax, k, k,
442       I           tFld, sFld,       &           tFld, sFld, bulkMod, myThid )
443       O           bulkMod,  
      I           myThid )  
               
             rholoc = rhoP0(i,j)  
      &           /(1. _d 0 -  
      &           pressure(i,j,k,bi,bj)*SItoBar/bulkMod(i,j))  
               
         
444              write(msgBuf,              write(msgBuf,
445       &           '(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)')
446       &           'rho(', sFld(i,j,k,bi,bj), ' PSU,',       &           'rho(', sFld(i,j,k,bi,bj), ' PSU,',
# Line 372  C     check nonlinear EOS Line 449  C     check nonlinear EOS
449       &           rloc(kcheck), bloc(kcheck)       &           rloc(kcheck), bloc(kcheck)
450              CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,              CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
451       &           SQUEEZE_RIGHT , 1)       &           SQUEEZE_RIGHT , 1)
452              write(msgBuf,'(a4,a32,f10.5,1x,f11.5)')              write(msgBuf,'(a14,a22,f10.5,1x,f11.5)')
453       &           'rho ', ' = ', rholoc, bulkMod(i,j)       &           'rho(find_rho) ',
454         &           ' = ', rholoc(i,j)+rhoConst, bulkMod(i,j)
455                CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
456         &           SQUEEZE_RIGHT , 1)
457    
458                call find_rho_scalar( tFld(i,j,k,bi,bj), sLoc(kcheck),
459         &           pLoc(kcheck), rhoLoc(i,j), myThid )
460                bulkMod(i,j) = 0. _d 0
461                write(msgBuf,'(a21,a15,f10.5,1x,f11.5)')
462         &           'rho(find_rho_scalar) ',
463         &           ' = ', rholoc(i,j)+rhoConst, bulkMod(i,j)
464              CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,              CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
465       &           SQUEEZE_RIGHT , 1)       &           SQUEEZE_RIGHT , 1)
466                            
# Line 389  C     end check nonlinear EOS Line 476  C     end check nonlinear EOS
476    
477        return        return
478        end        end
479    
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.7
  ViewVC Help
Powered by ViewVC 1.1.22