model/src/integr_continuity.F

C $Header: /u/gcmpack/MITgcm/model/src/integr_continuity.F,v 1.2 2002/12/10 03:00:59 jmc Exp $
C $Name:  $

#include "CPP_OPTIONS.h"

CBOP
C     !ROUTINE: INTEGR_CONTINUITY
C     !INTERFACE:
      SUBROUTINE INTEGR_CONTINUITY(
     I                             bi, bj, uFld, vFld,
     I                             myTime, myIter, myThid )
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE INTEGR_CONTINUITY                                 
C     | o Integrate the continuity Eq : compute vertical velocity 
C     |   and free surface "r-anomaly" (etaN) to satisfy  
C     |   exactly the convervation of the Total Volume            
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
C     == Global variables
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "DYNVARS.h"
#include "GRID.h"
#include "SURFACE.h"
#include "FFIELDS.h"

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     uFld    :: Zonal velocity ( m/s )
C     vFld    :: Meridional velocity ( m/s )
C     bi,bj   :: tile index
C     myTime  :: Current time in simulation
C     myIter  :: Current iteration number in simulation
C     myThid  :: Thread number for this instance of the routine.
      _RL myTime
      INTEGER myIter
      INTEGER myThid
      INTEGER bi,bj
      LOGICAL UpdateEtaN_EtaH
      _RL uFld(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
      _RL vFld(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy) 

C     !LOCAL VARIABLES:
C     Local variables in common block

C     Local variables
C     i,j,k  :: Loop counters
C     uTrans :: Volume transports ( uVel.xA )
C     vTrans :: Volume transports ( vVel.yA )
      INTEGER i,j,k
      _RL uTrans(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL vTrans(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL wSurf, facEmP
CEOP

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

#ifdef EXACT_CONSERV
      IF (exactConserv) THEN

C--   Compute the Divergence of The Barotropic Flow :

C-    Initialise 
      DO j=1-Oly,sNy+Oly
       DO i=1-Olx,sNx+Olx
         hDivFlow(i,j,bi,bj) = 0. _d 0
         utrans(i,j)         = 0. _d 0
         vtrans(i,j)         = 0. _d 0
       ENDDO
      ENDDO

      DO k=1,Nr
       
C-    Calculate velocity field "volume transports" through tracer cell faces
        DO j=1,sNy+1
         DO i=1,sNx+1
          uTrans(i,j) = uFld(i,j,k,bi,bj)*_dyG(i,j,bi,bj)
     &                 *drF(k)*_hFacW(i,j,k,bi,bj)
          vTrans(i,j) = vFld(i,j,k,bi,bj)*_dxG(i,j,bi,bj)
     &                 *drF(k)*_hFacS(i,j,k,bi,bj)
         ENDDO
        ENDDO

C-    Integrate vertically the Horizontal Divergence 
        DO j=1,sNy
         DO i=1,sNx
           hDivFlow(i,j,bi,bj) = hDivFlow(i,j,bi,bj)
     &       +maskC(i,j,k,bi,bj)*( uTrans(i+1,j)-uTrans(i,j)
     &                            +vTrans(i,j+1)-vTrans(i,j) )
         ENDDO
        ENDDO

C-    End DO k=1,Nr
      ENDDO

      ENDIF

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

      IF (exactConserv .AND. myTime.NE.startTime) THEN
C--   Update etaN at the end of the time step : 
C     Incorporate the Implicit part of -Divergence(Barotropic_Flow)

       IF (implicDiv2Dflow.EQ. 0. _d 0) THEN
        DO j=1-Oly,sNy+Oly
         DO i=1-Olx,sNx+Olx
           etaN(i,j,bi,bj) = etaH(i,j,bi,bj) 
         ENDDO
        ENDDO
       ELSEIF (useRealFreshWaterFlux) THEN
        DO j=1,sNy
         DO i=1,sNx
           etaN(i,j,bi,bj) = etaH(i,j,bi,bj)
     &      - implicDiv2Dflow*( convertEmP2rUnit*EmPmR(i,j,bi,bj)
     &                    +hDivFlow(i,j,bi,bj)*recip_rA(i,j,bi,bj)
     &                        )*deltaTfreesurf
         ENDDO
        ENDDO
       ELSE
        DO j=1,sNy
         DO i=1,sNx
           etaN(i,j,bi,bj) = etaH(i,j,bi,bj)
     &      - implicDiv2Dflow*hDivFlow(i,j,bi,bj)
     &                       *recip_rA(i,j,bi,bj)*deltaTfreesurf
         ENDDO
        ENDDO
       ENDIF

      ENDIF
#endif /* EXACT_CONSERV */

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

#ifdef NONLIN_FRSURF
      IF (select_rStar .NE. 0) THEN
        facEmP = 0.
        IF (useRealFreshWaterFlux) facEmP = convertEmP2rUnit
       
        IF ( useRealFreshWaterFlux .AND. myTime.EQ.startTime ) THEN
C     needs previous time-step value of E-P-R, that has not been loaded
C     and was not in pickup-file ; try to use etaN & etaH instead.
          DO j=1,sNy
           DO i=1,sNx
            rStarDhCDt(i,j,bi,bj) = 
     &         ( etaN(i,j,bi,bj)-etaH(i,j,bi,bj) )
     &           /(implicDiv2Dflow*deltaTfreesurf)
     &           * recip_Rcol(i,j,bi,bj)
           ENDDO
          ENDDO      
        ELSE
         DO j=1,sNy
          DO i=1,sNx
            rStarDhCDt(i,j,bi,bj) = 
     &         (-hDivFlow(i,j,bi,bj)*recip_rA(i,j,bi,bj)
     &          -facEmP*EmPmR(i,j,bi,bj)
     &         ) * recip_Rcol(i,j,bi,bj)
          ENDDO
         ENDDO
        ENDIF

      ENDIF
#endif /* NONLIN_FRSURF */

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

      DO k=Nr,1,-1
C--    Integrate continuity vertically for vertical velocity

       CALL INTEGRATE_FOR_W(
     I                       bi, bj, k, uFld, vFld,
     O                       wVel,
     I                       myThid )
 
#ifdef NONLIN_FRSURF
       IF ( k.EQ.Nr .AND. myTime.NE. 0. _d 0 .AND.
     &      useRealFreshWaterFlux .AND. 
     &      buoyancyRelation .EQ. 'OCEANICP' ) THEN

        IF ( myTime.NE.startTime ) THEN
          DO j=1,sNy
           DO i=1,sNx
             wVel(i,j,k,bi,bj) = wVel(i,j,k,bi,bj) 
     &         -convertEmP2rUnit*EmPmR(i,j,bi,bj)*maskC(i,j,k,bi,bj)
           ENDDO
          ENDDO
        ELSEIF (nonlinFreeSurf.GE.0) THEN
C     needs previous time-step value of E-P-R, that has not been loaded
C     and was not in pickup-file ; try to use etaN & etaH instead.
          DO j=1,sNy
           DO i=1,sNx
             wSurf = hDivFlow(i,j,bi,bj)*recip_rA(i,j,bi,bj)
     &             + (etaN(i,j,bi,bj)-etaH(i,j,bi,bj))
     &              /(implicDiv2Dflow*deltaTfreesurf)
             wVel(i,j,k,bi,bj) = wVel(i,j,k,bi,bj) 
     &                  + wSurf*maskC(i,j,k,bi,bj) 
           ENDDO
          ENDDO      
        ENDIF

       ENDIF
#endif /* NONLIN_FRSURF */

#ifdef ALLOW_OBCS
#ifdef ALLOW_NONHYDROSTATIC
C--    Apply OBC to W if in N-H mode
       IF (useOBCS.AND.nonHydrostatic) THEN
          CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )
        ENDIF
#endif /* ALLOW_NONHYDROSTATIC */
#endif /* ALLOW_OBCS */            

C-    End DO k=Nr,1,-1
      ENDDO

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

      RETURN
      END
1	jmc	1.3	C $Header: /u/gcmpack/MITgcm/model/src/integr_continuity.F,v 1.2 2002/12/10 03:00:59 jmc Exp $
2	jmc	1.1	C $Name: $
3
4			#include "CPP_OPTIONS.h"
5
6			CBOP
7			C !ROUTINE: INTEGR_CONTINUITY
8			C !INTERFACE:
9			SUBROUTINE INTEGR_CONTINUITY(
10			I bi, bj, uFld, vFld,
11			I myTime, myIter, myThid )
12			C !DESCRIPTION: \bv
13			C ==========================================================
14			C \| SUBROUTINE INTEGR_CONTINUITY
15			C \| o Integrate the continuity Eq : compute vertical velocity
16			C \| and free surface "r-anomaly" (etaN) to satisfy
17			C \| exactly the convervation of the Total Volume
18			C ==========================================================
19			C \ev
20
21			C !USES:
22			IMPLICIT NONE
23			C == Global variables
24			#include "SIZE.h"
25			#include "EEPARAMS.h"
26			#include "PARAMS.h"
27			#include "DYNVARS.h"
28			#include "GRID.h"
29			#include "SURFACE.h"
30			#include "FFIELDS.h"
31
32			C !INPUT/OUTPUT PARAMETERS:
33			C == Routine arguments ==
34			C uFld :: Zonal velocity ( m/s )
35			C vFld :: Meridional velocity ( m/s )
36			C bi,bj :: tile index
37			C myTime :: Current time in simulation
38			C myIter :: Current iteration number in simulation
39			C myThid :: Thread number for this instance of the routine.
40			_RL myTime
41			INTEGER myIter
42			INTEGER myThid
43			INTEGER bi,bj
44			LOGICAL UpdateEtaN_EtaH
45			_RL uFld(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
46			_RL vFld(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
47
48			C !LOCAL VARIABLES:
49			C Local variables in common block
50
51			C Local variables
52			C i,j,k :: Loop counters
53			C uTrans :: Volume transports ( uVel.xA )
54			C vTrans :: Volume transports ( vVel.yA )
55			INTEGER i,j,k
56			_RL uTrans(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
57			_RL vTrans(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
58	jmc	1.3	_RL wSurf, facEmP
59	jmc	1.1	CEOP
60
61			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
62
63			#ifdef EXACT_CONSERV
64			IF (exactConserv) THEN
65
66			C-- Compute the Divergence of The Barotropic Flow :
67
68			C- Initialise
69			DO j=1-Oly,sNy+Oly
70			DO i=1-Olx,sNx+Olx
71			hDivFlow(i,j,bi,bj) = 0. _d 0
72			utrans(i,j) = 0. _d 0
73			vtrans(i,j) = 0. _d 0
74			ENDDO
75			ENDDO
76
77			DO k=1,Nr
78
79			C- Calculate velocity field "volume transports" through tracer cell faces
80			DO j=1,sNy+1
81			DO i=1,sNx+1
82			uTrans(i,j) = uFld(i,j,k,bi,bj)*_dyG(i,j,bi,bj)
83			& drF(k)_hFacW(i,j,k,bi,bj)
84			vTrans(i,j) = vFld(i,j,k,bi,bj)*_dxG(i,j,bi,bj)
85			& drF(k)_hFacS(i,j,k,bi,bj)
86			ENDDO
87			ENDDO
88
89			C- Integrate vertically the Horizontal Divergence
90			DO j=1,sNy
91			DO i=1,sNx
92			hDivFlow(i,j,bi,bj) = hDivFlow(i,j,bi,bj)
93			& +maskC(i,j,k,bi,bj)*( uTrans(i+1,j)-uTrans(i,j)
94			& +vTrans(i,j+1)-vTrans(i,j) )
95			ENDDO
96			ENDDO
97
98			C- End DO k=1,Nr
99			ENDDO
100
101			ENDIF
102
103			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
104
105			IF (exactConserv .AND. myTime.NE.startTime) THEN
106			C-- Update etaN at the end of the time step :
107			C Incorporate the Implicit part of -Divergence(Barotropic_Flow)
108
109			IF (implicDiv2Dflow.EQ. 0. _d 0) THEN
110			DO j=1-Oly,sNy+Oly
111			DO i=1-Olx,sNx+Olx
112			etaN(i,j,bi,bj) = etaH(i,j,bi,bj)
113			ENDDO
114			ENDDO
115			ELSEIF (useRealFreshWaterFlux) THEN
116			DO j=1,sNy
117			DO i=1,sNx
118			etaN(i,j,bi,bj) = etaH(i,j,bi,bj)
119			& - implicDiv2Dflow( convertEmP2rUnitEmPmR(i,j,bi,bj)
120			& +hDivFlow(i,j,bi,bj)*recip_rA(i,j,bi,bj)
121			& )*deltaTfreesurf
122			ENDDO
123			ENDDO
124			ELSE
125			DO j=1,sNy
126			DO i=1,sNx
127			etaN(i,j,bi,bj) = etaH(i,j,bi,bj)
128			& - implicDiv2Dflow*hDivFlow(i,j,bi,bj)
129			& recip_rA(i,j,bi,bj)deltaTfreesurf
130			ENDDO
131			ENDDO
132			ENDIF
133
134			ENDIF
135			#endif /* EXACT_CONSERV */
136	jmc	1.3
137			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
138
139			#ifdef NONLIN_FRSURF
140			IF (select_rStar .NE. 0) THEN
141			facEmP = 0.
142			IF (useRealFreshWaterFlux) facEmP = convertEmP2rUnit
143
144			IF ( useRealFreshWaterFlux .AND. myTime.EQ.startTime ) THEN
145			C needs previous time-step value of E-P-R, that has not been loaded
146			C and was not in pickup-file ; try to use etaN & etaH instead.
147			DO j=1,sNy
148			DO i=1,sNx
149			rStarDhCDt(i,j,bi,bj) =
150			& ( etaN(i,j,bi,bj)-etaH(i,j,bi,bj) )
151			& /(implicDiv2Dflow*deltaTfreesurf)
152			& * recip_Rcol(i,j,bi,bj)
153			ENDDO
154			ENDDO
155			ELSE
156			DO j=1,sNy
157			DO i=1,sNx
158			rStarDhCDt(i,j,bi,bj) =
159			& (-hDivFlow(i,j,bi,bj)*recip_rA(i,j,bi,bj)
160			& -facEmP*EmPmR(i,j,bi,bj)
161			& ) * recip_Rcol(i,j,bi,bj)
162			ENDDO
163			ENDDO
164			ENDIF
165
166			ENDIF
167			#endif /* NONLIN_FRSURF */
168	jmc	1.1
169			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
170
171			DO k=Nr,1,-1
172			C-- Integrate continuity vertically for vertical velocity
173
174			CALL INTEGRATE_FOR_W(
175			I bi, bj, k, uFld, vFld,
176			O wVel,
177			I myThid )
178
179	jmc	1.2	#ifdef NONLIN_FRSURF
180			IF ( k.EQ.Nr .AND. myTime.NE. 0. _d 0 .AND.
181			& useRealFreshWaterFlux .AND.
182			& buoyancyRelation .EQ. 'OCEANICP' ) THEN
183
184			IF ( myTime.NE.startTime ) THEN
185			DO j=1,sNy
186			DO i=1,sNx
187			wVel(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)
188			& -convertEmP2rUnitEmPmR(i,j,bi,bj)maskC(i,j,k,bi,bj)
189			ENDDO
190			ENDDO
191			ELSEIF (nonlinFreeSurf.GE.0) THEN
192			C needs previous time-step value of E-P-R, that has not been loaded
193			C and was not in pickup-file ; try to use etaN & etaH instead.
194			DO j=1,sNy
195			DO i=1,sNx
196			wSurf = hDivFlow(i,j,bi,bj)*recip_rA(i,j,bi,bj)
197			& + (etaN(i,j,bi,bj)-etaH(i,j,bi,bj))
198			& /(implicDiv2Dflow*deltaTfreesurf)
199			wVel(i,j,k,bi,bj) = wVel(i,j,k,bi,bj)
200			& + wSurf*maskC(i,j,k,bi,bj)
201			ENDDO
202			ENDDO
203			ENDIF
204
205			ENDIF
206			#endif /* NONLIN_FRSURF */
207
208	jmc	1.1	#ifdef ALLOW_OBCS
209			#ifdef ALLOW_NONHYDROSTATIC
210			C-- Apply OBC to W if in N-H mode
211			IF (useOBCS.AND.nonHydrostatic) THEN
212			CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )
213			ENDIF
214			#endif /* ALLOW_NONHYDROSTATIC */
215			#endif /* ALLOW_OBCS */
216
217			C- End DO k=Nr,1,-1
218			ENDDO
219
220			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
221
222			RETURN
223			END