/[MITgcm]/MITgcm_contrib/bbl/code/mypackage_calc_rhs.F
ViewVC logotype

Diff of /MITgcm_contrib/bbl/code/mypackage_calc_rhs.F

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

revision 1.3 by dimitri, Sun Mar 6 02:11:05 2011 UTC revision 1.4 by dimitri, Sat Aug 6 02:10:27 2011 UTC
# Line 8  C     !ROUTINE: BBL_CALC_RHS Line 8  C     !ROUTINE: BBL_CALC_RHS
8    
9  C     !INTERFACE:  C     !INTERFACE:
10        SUBROUTINE MYPACKAGE_CALC_RHS(        SUBROUTINE MYPACKAGE_CALC_RHS(
11       I        bi, bj, myTime, myIter, myThid )       I        myTime, myIter, myThid )
12    
13  C     !DESCRIPTION:  C     !DESCRIPTION:
14  C     Calculate tendency of tracers due to bottom boundary layer.  C     Calculate tendency of tracers due to bottom boundary layer.
# Line 24  C     !USES: Line 24  C     !USES:
24  #include "BBL.h"  #include "BBL.h"
25    
26  C     !INPUT PARAMETERS:  C     !INPUT PARAMETERS:
 C     bi,bj     :: Tile indices  
27  C     myTime    :: Current time in simulation  C     myTime    :: Current time in simulation
28  C     myIter    :: Current time-step number  C     myIter    :: Current time-step number
29  C     myThid    :: my Thread Id number  C     myThid    :: my Thread Id number
       INTEGER bi, bj  
30        _RL     myTime        _RL     myTime
31        INTEGER myIter, myThid        INTEGER myIter, myThid
32    
33  C     !OUTPUT PARAMETERS:  C     !OUTPUT PARAMETERS:
34    
35  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
36    C     bi,bj     :: Tile indices
37  C     i,j       :: Loop indices  C     i,j       :: Loop indices
38  C     kBot      :: k index of bottommost wet grid box  C     kBot      :: k index of bottommost wet grid box
39  C     kLowC1    :: k index of bottommost (i,j) cell  C     kLowC1    :: k index of bottommost (i,j) cell
# Line 54  C     fZon      :: Zonal flux Line 53  C     fZon      :: Zonal flux
53  C     fMer      :: Meridional flux  C     fMer      :: Meridional flux
54  C     bbl_rho1  :: local (i,j) density  C     bbl_rho1  :: local (i,j) density
55  C     bbl_rho2  :: local (i+1, j) or (i,j+1) density  C     bbl_rho2  :: local (i+1, j) or (i,j+1) density
56          INTEGER bi, bj
57        INTEGER i, j, kBot, kLowC1, kLowC2, kl        INTEGER i, j, kBot, kLowC1, kLowC2, kl
58        _RL     resThk, resTheta, resSalt        _RL     resThk, resTheta, resSalt
59        _RL     deltaRho, deltaDpt, bbl_tend        _RL     deltaRho, deltaDpt, bbl_tend
# Line 67  C     bbl_rho2  :: local (i+1, j) or (i, Line 67  C     bbl_rho2  :: local (i+1, j) or (i,
67        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
68  CEOP  CEOP
69    
70    C--   Loops on tile indices bi,bj
71          DO bj=myByLo(myThid),myByHi(myThid)
72           DO bi=myBxLo(myThid),myBxHi(myThid)
73    
74  C     Initialize tendency terms and make local copy of  C     Initialize tendency terms and make local copy of
75  C     bottomost temperature, salinity, and in-situ desnity  C     bottomost temperature, salinity, and in-situ desnity
76  C     and of in-situ BBL density  C     and of in-situ BBL density
77        DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
78         DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
79          bbl_TendTheta(i,j,bi,bj) = 0. _d 0            bbl_TendTheta(i,j,bi,bj) = 0. _d 0
80          bbl_TendSalt (i,j,bi,bj) = 0. _d 0            bbl_TendSalt (i,j,bi,bj) = 0. _d 0
81          kBot = max(1,kLowC(i,j,bi,bj))            kBot = max(1,kLowC(i,j,bi,bj))
82          tLoc(i,j)   = theta(i,j,kBot,bi,bj)            tLoc(i,j)   = theta(i,j,kBot,bi,bj)
83          sLoc(i,j)   = salt (i,j,kBot,bi,bj)            sLoc(i,j)   = salt (i,j,kBot,bi,bj)
84          rholoc(i,j) = rhoInSitu(i,j,kBot,bi,bj)            rholoc(i,j) = rhoInSitu(i,j,kBot,bi,bj)
85          rhoBBL(i,j) = rhoInSitu(i,j,kBot+1, bi,bj)            rhoBBL(i,j) = rhoInSitu(i,j,kBot+1, bi,bj)
86         ENDDO           ENDDO
87        ENDDO          ENDDO
88    
89  C==== Compute and apply vertical exchange between BBL and  C==== Compute and apply vertical exchange between BBL and
90  C     residual volume of botommost wet grid box.  C     residual volume of botommost wet grid box.
91  C     This operation does not change total tracer quantity  C     This operation does not change total tracer quantity
92  C     in botommost wet grid box.  C     in botommost wet grid box.
93    
94        DO j=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
95         DO i=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
96          kBot = kLowC(i,j,bi,bj)            kBot = kLowC(i,j,bi,bj)
97          IF ( kBot .GT. 0 ) THEN            IF ( kBot .GT. 0 ) THEN
98           resThk = hFacC(i,j,kBot,bi,bj)*drF(kBot) - bbl_eta(i,j,bi,bj)             resThk = hFacC(i,j,kBot,bi,bj)*drF(kBot) - bbl_eta(i,j,bi,bj)
99    
100  C     If bbl occupies the complete bottom model grid box or  C     If bbl occupies the complete bottom model grid box or
101  C     if model density is higher than BBL then mix instantly.  C     if model density is higher than BBL then mix instantly.
102           IF ( (resThk.LE.0) .OR. (rhoLoc(i,j).GE.rhoBBL(i,j)) ) THEN             IF ( (resThk.LE.0) .OR. (rhoLoc(i,j).GE.rhoBBL(i,j)) ) THEN
103            bbl_theta(i,j,bi,bj) = tLoc(i,j)              bbl_theta(i,j,bi,bj) = tLoc(i,j)
104            bbl_salt (i,j,bi,bj) = sLoc(i,j)              bbl_salt (i,j,bi,bj) = sLoc(i,j)
105    
106  C     If model density is lower than BBL, slowly diffuse upward.  C     If model density is lower than BBL, slowly diffuse upward.
107           ELSE             ELSE
108            resTheta = ( tLoc(i,j) * (resThk+bbl_eta(i,j,bi,bj)) -              resTheta = ( tLoc(i,j) * (resThk+bbl_eta(i,j,bi,bj)) -
109       &         (bbl_theta(i,j,bi,bj)*bbl_eta(i,j,bi,bj)) ) / resThk       &           (bbl_theta(i,j,bi,bj)*bbl_eta(i,j,bi,bj)) ) / resThk
110            resSalt  = ( sLoc(i,j) * (resThk+bbl_eta(i,j,bi,bj)) -              resSalt  = ( sLoc(i,j) * (resThk+bbl_eta(i,j,bi,bj)) -
111       &         (bbl_salt (i,j,bi,bj)*bbl_eta(i,j,bi,bj)) ) / resThk       &           (bbl_salt (i,j,bi,bj)*bbl_eta(i,j,bi,bj)) ) / resThk
112            bbl_theta(i,j,bi,bj) = bbl_theta(i,j,bi,bj) +              bbl_theta(i,j,bi,bj) = bbl_theta(i,j,bi,bj) +
113       &         deltaT * (resTheta-bbl_theta(i,j,bi,bj)) / bbl_RelaxR       &           deltaT * (resTheta-bbl_theta(i,j,bi,bj)) / bbl_RelaxR
114            bbl_salt (i,j,bi,bj) = bbl_salt (i,j,bi,bj) +              bbl_salt (i,j,bi,bj) = bbl_salt (i,j,bi,bj) +
115       &         deltaT * (resSalt -bbl_salt (i,j,bi,bj)) / bbl_RelaxR       &           deltaT * (resSalt -bbl_salt (i,j,bi,bj)) / bbl_RelaxR
116           ENDIF             ENDIF
117          ENDIF            ENDIF
118         ENDDO           ENDDO
119        ENDDO          ENDDO
120    
121  C==== Compute zonal bbl exchange.  C==== Compute zonal bbl exchange.
122        DO j=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
123         DO i=1-Olx,sNx+Olx-1           DO i=1-Olx,sNx+Olx-1
124          kLowC1 = kLowC(i,j,bi,bj)            kLowC1 = kLowC(i,j,bi,bj)
125          kLowC2 = kLowC(i+1,j,bi,bj)            kLowC2 = kLowC(i+1,j,bi,bj)
126          IF ((kLowC1.GT.0).AND.(kLowC2.GT.0)) THEN            IF ((kLowC1.GT.0).AND.(kLowC2.GT.0)) THEN
127  C compare the bbl densities at the higher pressure (highest possible density of given t,s)  C     Compare the bbl densities at the higher pressure
128  C bbl in situ density is stored in kLowC + 1 index    C     (highest possible density of given t,s)
129           kl = MAX(kLowC1, kLowC2) + 1  C     bbl in situ density is stored in kLowC + 1 index  
130           bbl_rho1 = rhoInSitu(i,j,kl,bi,bj)             kl = MAX(kLowC1, kLowC2) + 1
131           bbl_rho2 = rhoInSitu(i+1,j,kl,bi,bj)                     bbl_rho1 = rhoInSitu(i,j,kl,bi,bj)
132           deltaRho = bbl_rho2 - bbl_rho1             bbl_rho2 = rhoInSitu(i+1,j,kl,bi,bj)        
133           deltaDpt = R_low(i  ,j,bi,bj) + bbl_eta(i  ,j,bi,bj) -             deltaRho = bbl_rho2 - bbl_rho1
134       &              R_low(i+1,j,bi,bj) - bbl_eta(i+1,j,bi,bj)             deltaDpt = R_low(i  ,j,bi,bj) + bbl_eta(i  ,j,bi,bj) -
135         &          R_low(i+1,j,bi,bj) - bbl_eta(i+1,j,bi,bj)
136    
137  C     If heavy BBL water is higher than light BBL water,  C     If heavy BBL water is higher than light BBL water,
138  C     exchange properties laterally.  C     exchange properties laterally.
139           IF ( (deltaRho*deltaDpt) .LE. 0. ) THEN             IF ( (deltaRho*deltaDpt) .LE. 0. ) THEN
140            bbl_TendTheta(i,j,bi,bj) = bbl_TendTheta(i,j,bi,bj) +              bbl_TendTheta(i,j,bi,bj) = bbl_TendTheta(i,j,bi,bj) +
141       &         ( bbl_theta(i+1,j,bi,bj) - bbl_theta(i,j,bi,bj) ) /       &           ( bbl_theta(i+1,j,bi,bj) - bbl_theta(i,j,bi,bj) ) /
142       &         bbl_RelaxH       &           bbl_RelaxH
143            bbl_TendTheta(i+1,j,bi,bj) = bbl_TendTheta(i+1,j,bi,bj) -              bbl_TendTheta(i+1,j,bi,bj) = bbl_TendTheta(i+1,j,bi,bj) -
144       &         ( bbl_theta(i+1,j,bi,bj) - bbl_theta(i,j,bi,bj) ) *       &           ( bbl_theta(i+1,j,bi,bj) - bbl_theta(i,j,bi,bj) ) *
145       &         ( rA(i  ,j,bi,bj) * bbl_eta(i  ,j,bi,bj) ) /       &           ( rA(i  ,j,bi,bj) * bbl_eta(i  ,j,bi,bj) ) /
146       &         ( rA(i+1,j,bi,bj) * bbl_eta(i+1,j,bi,bj) * bbl_RelaxH )       &           ( rA(i+1,j,bi,bj) * bbl_eta(i+1,j,bi,bj) * bbl_RelaxH )
147            bbl_TendSalt(i,j,bi,bj) = bbl_TendSalt(i,j,bi,bj) +              bbl_TendSalt(i,j,bi,bj) = bbl_TendSalt(i,j,bi,bj) +
148       &         ( bbl_salt(i+1,j,bi,bj) - bbl_salt(i,j,bi,bj) ) /       &           ( bbl_salt(i+1,j,bi,bj) - bbl_salt(i,j,bi,bj) ) /
149       &         bbl_RelaxH       &           bbl_RelaxH
150            bbl_TendSalt(i+1,j,bi,bj) = bbl_TendSalt(i+1,j,bi,bj) -              bbl_TendSalt(i+1,j,bi,bj) = bbl_TendSalt(i+1,j,bi,bj) -
151       &         ( bbl_salt(i+1,j,bi,bj) - bbl_salt(i,j,bi,bj) ) *       &           ( bbl_salt(i+1,j,bi,bj) - bbl_salt(i,j,bi,bj) ) *
152       &         ( rA(i  ,j,bi,bj) * bbl_eta(i  ,j,bi,bj) ) /       &           ( rA(i  ,j,bi,bj) * bbl_eta(i  ,j,bi,bj) ) /
153       &         ( rA(i+1,j,bi,bj) * bbl_eta(i+1,j,bi,bj) * bbl_RelaxH )       &           ( rA(i+1,j,bi,bj) * bbl_eta(i+1,j,bi,bj) * bbl_RelaxH )
154           ENDIF             ENDIF
155          ENDIF            ENDIF
156         ENDDO           ENDDO
157        ENDDO          ENDDO
158    
159  C==== Compute meridional bbl exchange.  C==== Compute meridional bbl exchange.
160        DO j=1-Oly,sNy+Oly-1          DO j=1-Oly,sNy+Oly-1
161         DO i=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
162          kLowC1 = kLowC(i,j,bi,bj)            kLowC1 = kLowC(i,j,bi,bj)
163          kLowC2 = kLowC(i,j+1, bi,bj)            kLowC2 = kLowC(i,j+1, bi,bj)
164          IF ((kLowC1.GT.0).AND.(kLowC2.GT.0)) THEN            IF ((kLowC1.GT.0).AND.(kLowC2.GT.0)) THEN
165  C compare the bbl densities at the higher pressure (highest possible density of given t,s)  C     compare the bbl densities at the higher pressure
166  C bbl in situ density is stored in kLowC + 1 index    C     (highest possible density of given t,s)
167           kl = MAX(kLowC1, kLowC2) + 1  C     bbl in situ density is stored in kLowC + 1 index  
168           bbl_rho1 = rhoInSitu(i,j,kl,bi,bj)             kl = MAX(kLowC1, kLowC2) + 1
169           bbl_rho2 = rhoInSitu(i,j+1,kl,bi,bj)                     bbl_rho1 = rhoInSitu(i,j,kl,bi,bj)
170           deltaRho = bbl_rho2 - bbl_rho1             bbl_rho2 = rhoInSitu(i,j+1,kl,bi,bj)        
171           deltaDpt = R_low(i,j  ,bi,bj) + bbl_eta(i,j  ,bi,bj) -             deltaRho = bbl_rho2 - bbl_rho1
172       &              R_low(i,j+1,bi,bj) - bbl_eta(i,j+1,bi,bj)             deltaDpt = R_low(i,j  ,bi,bj) + bbl_eta(i,j  ,bi,bj) -
173         &          R_low(i,j+1,bi,bj) - bbl_eta(i,j+1,bi,bj)
174    
175  C     If heavy BBL water is higher than light BBL water,  C     If heavy BBL water is higher than light BBL water,
176  C     exchange properties laterally.  C     exchange properties laterally.
177           IF ( (deltaRho*deltaDpt) .LE. 0. ) THEN             IF ( (deltaRho*deltaDpt) .LE. 0. ) THEN
178            bbl_TendTheta(i,j,bi,bj) = bbl_TendTheta(i,j,bi,bj) +              bbl_TendTheta(i,j,bi,bj) = bbl_TendTheta(i,j,bi,bj) +
179       &         ( bbl_theta(i,j+1,bi,bj) - bbl_theta(i,j,bi,bj) ) /       &           ( bbl_theta(i,j+1,bi,bj) - bbl_theta(i,j,bi,bj) ) /
180       &         bbl_RelaxH       &           bbl_RelaxH
181            bbl_TendTheta(i,j+1,bi,bj) = bbl_TendTheta(i,j+1,bi,bj) -              bbl_TendTheta(i,j+1,bi,bj) = bbl_TendTheta(i,j+1,bi,bj) -
182       &         ( bbl_theta(i,j+1,bi,bj) - bbl_theta(i,j,bi,bj) ) *       &           ( bbl_theta(i,j+1,bi,bj) - bbl_theta(i,j,bi,bj) ) *
183       &         ( rA(i  ,j,bi,bj) * bbl_eta(i  ,j,bi,bj) ) /       &           ( rA(i  ,j,bi,bj) * bbl_eta(i  ,j,bi,bj) ) /
184       &         ( rA(i,j+1,bi,bj) * bbl_eta(i,j+1,bi,bj) ) /       &           ( rA(i,j+1,bi,bj) * bbl_eta(i,j+1,bi,bj) ) /
185       &         bbl_RelaxH       &           bbl_RelaxH
186            bbl_TendSalt(i,j,bi,bj) = bbl_TendSalt(i,j,bi,bj) +              bbl_TendSalt(i,j,bi,bj) = bbl_TendSalt(i,j,bi,bj) +
187       &         ( bbl_salt(i,j+1,bi,bj) - bbl_salt(i,j,bi,bj) ) /       &           ( bbl_salt(i,j+1,bi,bj) - bbl_salt(i,j,bi,bj) ) /
188       &         bbl_RelaxH       &           bbl_RelaxH
189            bbl_TendSalt(i,j+1,bi,bj) = bbl_TendSalt(i,j+1,bi,bj) -              bbl_TendSalt(i,j+1,bi,bj) = bbl_TendSalt(i,j+1,bi,bj) -
190       &         ( bbl_salt(i,j+1,bi,bj)-bbl_salt(i,j,bi,bj)) *       &           ( bbl_salt(i,j+1,bi,bj)-bbl_salt(i,j,bi,bj)) *
191       &         ( rA(i  ,j,bi,bj) * bbl_eta(i  ,j,bi,bj) ) /       &           ( rA(i  ,j,bi,bj) * bbl_eta(i  ,j,bi,bj) ) /
192       &         ( rA(i,j+1,bi,bj) * bbl_eta(i,j+1,bi,bj) * bbl_RelaxH )       &           ( rA(i,j+1,bi,bj) * bbl_eta(i,j+1,bi,bj) * bbl_RelaxH )
193           ENDIF             ENDIF
194          ENDIF            ENDIF
195         ENDDO           ENDDO
196        ENDDO          ENDDO
197    
198  C==== Apply lateral BBL exchange then scale tendency term  C==== Apply lateral BBL exchange then scale tendency term
199  C     for botommost wet grid box.  C     for botommost wet grid box.
200        DO j=1-Oly,sNy+Oly-1          DO j=1-Oly,sNy+Oly-1
201         DO i=1-Olx,sNx+Olx-1           DO i=1-Olx,sNx+Olx-1
202          kBot = kLowC(i,j,bi,bj)            kBot = kLowC(i,j,bi,bj)
203          IF ( kBot .GT. 0 ) THEN            IF ( kBot .GT. 0 ) THEN
204           bbl_theta(i,j,bi,bj) = bbl_theta(i,j,bi,bj) +             bbl_theta(i,j,bi,bj) = bbl_theta(i,j,bi,bj) +
205       &        deltaT * bbl_TendTheta(i,j,bi,bj)       &          deltaT * bbl_TendTheta(i,j,bi,bj)
206           bbl_salt (i,j,bi,bj) = bbl_salt (i,j,bi,bj) +             bbl_salt (i,j,bi,bj) = bbl_salt (i,j,bi,bj) +
207       &        deltaT * bbl_TendSalt (i,j,bi,bj)       &          deltaT * bbl_TendSalt (i,j,bi,bj)
208           bbl_TendTheta(i,j,bi,bj) = bbl_TendTheta(i,j,bi,bj) *             bbl_TendTheta(i,j,bi,bj) = bbl_TendTheta(i,j,bi,bj) *
209       &        bbl_eta(i,j,bi,bj) / (hFacC(i,j,kBot,bi,bj)*drF(kBot))       &          bbl_eta(i,j,bi,bj) / (hFacC(i,j,kBot,bi,bj)*drF(kBot))
210           bbl_TendSalt (i,j,bi,bj) = bbl_TendSalt (i,j,bi,bj) *             bbl_TendSalt (i,j,bi,bj) = bbl_TendSalt (i,j,bi,bj) *
211       &        bbl_eta(i,j,bi,bj) / (hFacC(i,j,kBot,bi,bj)*drF(kBot))       &          bbl_eta(i,j,bi,bj) / (hFacC(i,j,kBot,bi,bj)*drF(kBot))
212          ENDIF            ENDIF
213         ENDDO           ENDDO
214        ENDDO          ENDDO
215    
216  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
217        IF ( debugLevel .GE. debLevB ) THEN          IF ( debugLevel .GE. debLevB ) THEN
218  C     Check salinity conservation  C     Check salinity conservation
219         bbl_tend=0           bbl_tend=0
220         DO j=1,sNy           DO j=1,sNy
221          DO i=1,sNx            DO i=1,sNx
222           kBot = kLowC(i,j,bi,bj)             kBot = kLowC(i,j,bi,bj)
223           IF ( kBot .GT. 0 ) THEN             IF ( kBot .GT. 0 ) THEN
224            bbl_tend = bbl_tend + bbl_TendSalt(i,j,bi,bj) *              bbl_tend = bbl_tend + bbl_TendSalt(i,j,bi,bj) *
225       &         hFacC(i,j,kBot,bi,bj) * drF(kBot) *rA(i,j,bi,bj)       &           hFacC(i,j,kBot,bi,bj) * drF(kBot) *rA(i,j,bi,bj)
226           ENDIF             ENDIF
227          ENDDO            ENDDO
228         ENDDO           ENDDO
229         _GLOBAL_SUM_RL( bbl_tend, myThid )           _GLOBAL_SUM_RL( bbl_tend, myThid )
230         WRITE(msgBuf,'(A,E10.2)') 'total salt tendency = ', bbl_tend           WRITE(msgBuf,'(A,E10.2)') 'total salt tendency = ', bbl_tend
231         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
232       &      SQUEEZE_RIGHT, myThid )       &        SQUEEZE_RIGHT, myThid )
233        ENDIF          ENDIF
234  #endif /* ALLOW_DEBUG */  #endif /* ALLOW_DEBUG */
235    
236            CALL EXCH_XY_RL( bbl_theta, myThid )
237            CALL EXCH_XY_RL( bbl_salt , myThid )
238    
239        CALL EXCH_XY_RL( bbl_theta, myThid )  C--   end bi,bj loops.
240        CALL EXCH_XY_RL( bbl_salt , myThid )         ENDDO
241          ENDDO
242    
243        RETURN        RETURN
244        END        END
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.22