/[MITgcm]/MITgcm/model/src/set_ref_state.F
ViewVC logotype

Diff of /MITgcm/model/src/set_ref_state.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.7 by jmc, Thu Feb 4 01:16:03 2016 UTC revision 1.8 by jmc, Thu Mar 10 20:54:57 2016 UTC
# Line 39  C     == Local variables == Line 39  C     == Local variables ==
39  C     msgBuf :: Informational/error message buffer  C     msgBuf :: Informational/error message buffer
40        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
41        INTEGER k, ks, stdUnit        INTEGER k, ks, stdUnit
42        _RL rHalf(2*Nr+1)        _RL rHalf(2*Nr+1), pRefLocF(Nr+1)
43        _RL rhoRef(Nr), tLoc(Nr)        _RL rhoRef(Nr), tLoc(Nr)
44        _RL pLoc, rhoUp, rhoDw, rhoLoc        _RL pLoc, rhoUp, rhoDw, rhoLoc
45        _RL ddPI, conv_theta2T, thetaLoc        _RL ddPI, conv_theta2T, thetaLoc
46    C--
47          _RL     maxResid
48          INTEGER maxIterNb, belowCritNb
49  CEOP  CEOP
50    
51        _BEGIN_MASTER( myThid )        _BEGIN_MASTER( myThid )
# Line 56  C--   Initialise: Line 59  C--   Initialise:
59        DO k=1,Nr        DO k=1,Nr
60          rhoRef(k)  = 0.          rhoRef(k)  = 0.
61          dBdrRef(k) = 0.          dBdrRef(k) = 0.
62            pRef4EOS(k)  = 0.
63          rHalf(2*k-1) = rF(k)          rHalf(2*k-1) = rF(k)
64          rHalf(2*k)   = rC(k)          rHalf(2*k)   = rC(k)
65        ENDDO        ENDDO
# Line 66  C--   Initialise: Line 70  C--   Initialise:
70          wUnit2rVel(k) = 1. _d 0          wUnit2rVel(k) = 1. _d 0
71        ENDDO        ENDDO
72    
73  C--   Initialise density factor for anelastic formulation:  C-    Initialise density factor for anelastic formulation:
74        DO k=1,Nr        DO k=1,Nr
75          rhoFacC(k) = 1. _d 0          rhoFacC(k) = 1. _d 0
76          recip_rhoFacC(k) = 1. _d 0          recip_rhoFacC(k) = 1. _d 0
# Line 77  C--   Initialise density factor for anel Line 81  C--   Initialise density factor for anel
81        ENDDO        ENDDO
82    
83  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
84    C--   Oceanic: define reference pressure/geo-potential vertical scale:
85          IF ( buoyancyRelation.EQ.'OCEANIC' ) THEN
86            phiRef(1)   = phi0Ref
87            pRefLocF(1) = rhoConst*phi0Ref
88            IF ( gravityFile.EQ.' ' ) THEN
89             DO k=1,Nr
90              phiRef(2*k)   = phi0Ref + (rC(k) - rF(1))*gravity*gravitySign
91              phiRef(2*k+1) = phi0Ref + (rF(k+1)-rF(1))*gravity*gravitySign
92    c         pRef4EOS(k)   = rhoConst*phiRef(2*k)
93    c         pRefLocF(k+1) = rhoConst*phiRef(2*k+1)
94    C note: just to get the same previous machine truncation:
95              pRef4EOS(k)   = rhoConst*phi0Ref
96         &                  + rhoConst*(rC(k) - rF(1))*gravity*gravitySign
97              pRefLocF(k+1) = rhoConst*phi0Ref
98         &                  + rhoConst*(rF(k+1)-rF(1))*gravity*gravitySign
99             ENDDO
100            ELSEIF ( integr_GeoPot.EQ.1 ) THEN
101             DO k=1,Nr
102              phiRef(2*k)   = phiRef(2*k-1)
103         &                           + halfRL*drF(k)*gravity*gravFacC(k)
104              phiRef(2*k+1) = phiRef(2*k-1) + drF(k)*gravity*gravFacC(k)
105              pRef4EOS(k)   = rhoConst*phiRef(2*k)
106              pRefLocF(k+1) = rhoConst*phiRef(2*k+1)
107             ENDDO
108            ELSE
109             phiRef(2)   = phi0Ref + drC(1)*gravity*gravFacF(1)
110             DO k=2,Nr
111              phiRef(2*k-1) = phiRef(2*k-2)
112         &                           + halfRL*drC(k)*gravity*gravFacF(k)
113              phiRef(2*k)   = phiRef(2*k-2) + drC(k)*gravity*gravFacF(k)
114             ENDDO
115             k = Nr
116             phiRef(2*k+1) = phiRef(2*k) + drC(k+1)*gravity*gravFacF(k+1)
117             DO k=1,Nr
118              pRef4EOS(k)   = rhoConst*phiRef(2*k)
119              pRefLocF(k+1) = rhoConst*phiRef(2*k+1)
120             ENDDO
121            ENDIF
122          ELSEIF ( buoyancyRelation.EQ.'OCEANICP' ) THEN
123            phiRef(1) = phi0Ref
124            DO k=1,Nr
125              phiRef(2*k)   = phi0Ref - recip_rhoConst*(rC(k) - rF(1))
126              phiRef(2*k+1) = phi0Ref - recip_rhoConst*(rF(k+1)-rF(1))
127            ENDDO
128    c     ELSEIF (buoyancyRelation .EQ. 'ATMOSPHERIC') THEN
129    C     phiRef is computed later (see below)
130          ENDIF
131    
132    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
133        IF ( eosType.EQ.'POLY3' ) THEN        IF ( eosType.EQ.'POLY3' ) THEN
134         IF ( implicitIntGravWave ) THEN         IF ( implicitIntGravWave ) THEN
135           WRITE(msgBuf,'(2A)') 'SET_REF_STATE:',           WRITE(msgBuf,'(2A)') 'SET_REF_STATE:',
# Line 87  C---+----1----+----2----+----3----+----4 Line 140  C---+----1----+----2----+----3----+----4
140           CALL PRINT_ERROR( msgBuf , myThid )           CALL PRINT_ERROR( msgBuf , myThid )
141           STOP 'ABNORMAL END: S/R SET_REF_STATE'           STOP 'ABNORMAL END: S/R SET_REF_STATE'
142         ELSEIF ( nonHydrostatic .AND.         ELSEIF ( nonHydrostatic .AND.
143       &          buoyancyRelation .EQ. 'OCEANICP' ) THEN       &          buoyancyRelation.EQ.'OCEANICP' ) THEN
144           WRITE(msgBuf,'(2A)') 'SET_REF_STATE:',           WRITE(msgBuf,'(2A)') 'SET_REF_STATE:',
145       &    ' need to compute reference density for Non-Hyd'       &    ' need to compute reference density for Non-Hyd'
146           CALL PRINT_ERROR( msgBuf , myThid )           CALL PRINT_ERROR( msgBuf , myThid )
# Line 106  C---+----1----+----2----+----3----+----4 Line 159  C---+----1----+----2----+----3----+----4
159       &                       SQUEEZE_RIGHT , myThid)       &                       SQUEEZE_RIGHT , myThid)
160         ENDIF         ENDIF
161  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
162        ELSEIF (buoyancyRelation .EQ. 'OCEANIC') THEN        ELSEIF ( buoyancyRelation.EQ.'OCEANIC' ) THEN
163    
164  C--   Compute reference density profile and reference stratification  C--   Compute reference density profile
165          DO k=1,Nr          DO k=1,Nr
166            pLoc = -rhoConst*rC(k)*gravity            pLoc = pRef4EOS(k)
167            CALL FIND_RHO_SCALAR(            CALL FIND_RHO_SCALAR(
168       I                          tRef(k), sRef(k), pLoc,       I                          tRef(k), sRef(k), pLoc,
169       O                          rhoRef(k), myThid )       O                          rhoRef(k), myThid )
170          ENDDO          ENDDO
171    
172            IF ( selectP_inEOS_Zc.EQ.1 ) THEN
173    C--   Find reference pressure in hydrostatic balance with updated ref density
174              maxResid = rhoConst* 1. _d -14
175              belowCritNb = 5
176              maxIterNb = 10*Nr
177              CALL FIND_HYD_PRESS_1D(
178         O                            pRef4EOS, pRefLocF,
179         U                            rhoRef,
180         I                            tRef, sRef, maxResid,
181         I                            belowCritNb, maxIterNb, myThid )
182            ENDIF
183    
184  C--   Compute reference stratification: N^2 = -(g/rho_c) * d.rho/dz @ const. p  C--   Compute reference stratification: N^2 = -(g/rho_c) * d.rho/dz @ const. p
185          dBdrRef(1) = 0. _d 0          dBdrRef(1) = 0. _d 0
186          DO k=2,Nr          DO k=2,Nr
187            pLoc = -rhoConst*rF(k)*gravity            pLoc = pRefLocF(k)
188            CALL FIND_RHO_SCALAR(            CALL FIND_RHO_SCALAR(
189       I                          tRef(k-1), sRef(k-1), pLoc,       I                          tRef(k-1), sRef(k-1), pLoc,
190       O                          rhoUp, myThid )       O                          rhoUp, myThid )
# Line 127  C--   Compute reference stratification: Line 192  C--   Compute reference stratification:
192       I                          tRef(k), sRef(k), pLoc,       I                          tRef(k), sRef(k), pLoc,
193       O                          rhoDw, myThid )       O                          rhoDw, myThid )
194            dBdrRef(k) = (rhoDw - rhoUp)*recip_drC(k)            dBdrRef(k) = (rhoDw - rhoUp)*recip_drC(k)
195       &               *recip_rhoConst*gravity       &               *recip_rhoConst*gravity*gravFacF(k)
196            IF (eosType .EQ. 'LINEAR') THEN            IF ( eosType.EQ.'LINEAR' ) THEN
197  C- get more precise values (differences from above are due to machine round-off)  C- get more precise values (differences from above are due to machine round-off)
198              dBdrRef(k) = ( sBeta *(sRef(k)-sRef(k-1))              dBdrRef(k) = ( sBeta *(sRef(k)-sRef(k-1))
199       &                    -tAlpha*(tRef(k)-tRef(k-1))       &                    -tAlpha*(tRef(k)-tRef(k-1))
200       &                   )*recip_drC(k)       &                   )*recip_drC(k)
201       &                 *rhoNil*recip_rhoConst*gravity       &                 *rhoNil*recip_rhoConst*gravity*gravFacF(k)
202            ENDIF            ENDIF
203          ENDDO          ENDDO
204    
205  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
206        ELSEIF (buoyancyRelation .EQ. 'OCEANICP') THEN        ELSEIF ( buoyancyRelation.EQ.'OCEANICP' ) THEN
207    
208  C--   Compute reference density profile  C--   Compute reference density profile
209          DO k=1,Nr          DO k=1,Nr
# Line 176  C     note: wUnit2rVel & rVel2wUnit repl Line 241  C     note: wUnit2rVel & rVel2wUnit repl
241          ENDDO          ENDDO
242    
243  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
244        ELSEIF (buoyancyRelation .EQ. 'ATMOSPHERIC') THEN        ELSEIF ( buoyancyRelation.EQ.'ATMOSPHERIC' ) THEN
245    
246  C--   Compute reference stratification: -d.alpha/dp @ constant p  C--   Compute reference stratification: -d.alpha/dp @ constant p
247          dBdrRef(1) = 0. _d 0          dBdrRef(1) = 0. _d 0
# Line 214  C     note: wUnit2rVel & rVel2wUnit repl Line 279  C     note: wUnit2rVel & rVel2wUnit repl
279  C-    Compute Reference Geopotential at Half levels :  C-    Compute Reference Geopotential at Half levels :
280  C      Tracer level k: phiRef(2k)  ;  Interface_W level k: phiRef(2k-1)  C      Tracer level k: phiRef(2k)  ;  Interface_W level k: phiRef(2k-1)
281    
282         phiRef(1) = 0. _d 0         phiRef(1) = phi0Ref
283         IF ( select_rStar.GE.1 .OR. selectSigmaCoord.GE.1 ) THEN         IF ( select_rStar.GE.1 .OR. selectSigmaCoord.GE.1 ) THEN
284  C-      isothermal (theta=const) reference state  C-      isothermal (theta=const) reference state
285          DO k=1,Nr          DO k=1,Nr
# Line 263  C---+----1----+----2----+----3----+----4 Line 328  C---+----1----+----2----+----3----+----4
328  C--   endif buoyancyRelation  C--   endif buoyancyRelation
329        ENDIF        ENDIF
330    
 C--   fill-in phiRef array (presently not used)  
       IF ( buoyancyRelation.EQ.'OCEANIC' ) THEN  
 c       phiRef(1) = gravitySign*gravity*(rF(1)-Ro_SeaLevel)  
         phiRef(1) = 0. _d 0  
         DO k=1,Nr  
           phiRef(2*k)   = gravitySign*gravity*(rC(k) - Ro_SeaLevel)  
           phiRef(2*k+1) = gravitySign*gravity*(rF(k+1)-Ro_SeaLevel)  
         ENDDO  
       ELSEIF ( buoyancyRelation.EQ.'OCEANICP' ) THEN  
 C     should be : phiRef = phi_Origin - (rC - Ro_SeaLevel)/rhoConst  
 C-    but since the reference geopotential "phi_Origin" @ p = rF(k=1)  
 C     is not currently stored, we only get a relative geopotential:  
         phiRef(1) = -recip_rhoConst*rF(1)  
         DO k=1,Nr  
           phiRef(2*k)   = -recip_rhoConst*rC(k)  
           phiRef(2*k+1) = -recip_rhoConst*rF(k+1)  
         ENDDO  
       ENDIF  
   
331  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
332    
333        IF ( usingZCoords .AND. rhoRefFile .NE. ' ' ) THEN        IF ( usingZCoords .AND. rhoRefFile .NE. ' ' ) THEN
# Line 291  C       surface-interface rho-factor has Line 337  C       surface-interface rho-factor has
337  C       rhoFac(k) = density ratio between layer k and top interface  C       rhoFac(k) = density ratio between layer k and top interface
338          DO k=1,Nr          DO k=1,Nr
339            rhoFacC(k) = rho1Ref(k)/rhoConst            rhoFacC(k) = rho1Ref(k)/rhoConst
340    c         rhoFacC(k) = rho1Ref(k)*recip_rhoConst
341          ENDDO          ENDDO
342          DO k=2,Nr          DO k=2,Nr
343  C       since rkSign=-1, recip_drC(k) = 1./(rC(k-1)-rC(k))  C       since rkSign=-1, recip_drC(k) = 1./(rC(k-1)-rC(k))
# Line 313  C-      set reciprocal rho-factor: Line 360  C-      set reciprocal rho-factor:
360  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
361    
362  C--   Write to check :  C--   Write to check :
363        IF (buoyancyRelation .EQ. 'ATMOSPHERIC') THEN        IF ( gravityFile .NE. ' ' ) THEN
364    C-    write gravity vertical profile factor to binary file:
365            CALL WRITE_GLVEC_RL( 'GravFacC',' ',gravFacC, Nr, -1, myThid )
366            CALL WRITE_GLVEC_RL( 'GravFacF',' ',gravFacF,Nr+1,-1, myThid )
367          ENDIF
368          IF ( selectP_inEOS_Zc.EQ.1 ) THEN
369    C-    write (to binary file) Hyd-Pressure used in EOS:
370            CALL WRITE_GLVEC_RL( 'PRef4EOS',' ',pRef4EOS, Nr, -1, myThid )
371    c       CALL WRITE_GLVEC_RL( 'PRefLocF',' ',pRefLocF,Nr+1,-1, myThid )
372          ENDIF
373          IF ( buoyancyRelation.EQ.'ATMOSPHERIC' ) THEN
374         WRITE(msgBuf,'(A)') ' '         WRITE(msgBuf,'(A)') ' '
375         CALL PRINT_MESSAGE( msgBuf, stdUnit, SQUEEZE_RIGHT, myThid )         CALL PRINT_MESSAGE( msgBuf, stdUnit, SQUEEZE_RIGHT, myThid )
376         WRITE(msgBuf,'(A)')         WRITE(msgBuf,'(A)')
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.8
  ViewVC Help
Powered by ViewVC 1.1.22