model/src/momentum_correction_step.F

C $Header: /u/gcmpack/MITgcm/model/src/momentum_correction_step.F,v 1.7 2009/04/28 18:01:14 jmc Exp $
C $Name:  $

#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"

CBOP
C     !ROUTINE: MOMENTUM_CORRECTION_STEP
C     !INTERFACE:
      SUBROUTINE MOMENTUM_CORRECTION_STEP(myTime, myIter, myThid)
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE MOMENTUM_CORRECTION_STEP
C     *==========================================================*
C     |1rst Part : Update U,V.
C     |
C     | The arrays used for time stepping are cycled.
C     | Momentum:
C     |           V(n) = Gv(n) - dt * grad Eta
C     |
C     |part1: update U,V
C     |  U*,V* (contained in gU,gV) have the surface
C     |     pressure gradient term added and the result stored
C     |     in U,V (contained in uVel, vVel)
C     |
C     |part2: Adjustments
C     |   o Filter  U,V (Shapiro Filter, Zonal_Filter)
C     |   o Compute again Eta (exact volume conservation)
C     |   o Compute vertical velocity
C     *==========================================================*
C     \ev

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

#ifdef ALLOW_SHAP_FILT
#include "SHAP_FILT.h"
#endif
#ifdef ALLOW_ZONAL_FILT
#include "ZONAL_FILT.h"
#endif
#ifdef ALLOW_AUTODIFF_TAMC
#include "tamc.h"
#include "tamc_keys.h"
# ifdef NONLIN_FRSURF
#  include "SURFACE.h"
# endif
#endif

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
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

C     !LOCAL VARIABLES:
C     == Local variables
      _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      INTEGER iMin,iMax
      INTEGER jMin,jMax
      INTEGER bi,bj
      INTEGER k,i,j
CEOP

      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)

C--     Set up work arrays that need valid initial values
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          phiSurfX(i,j)=0.
          phiSurfY(i,j)=0.
         ENDDO
        ENDDO

C       Loop range: Gradients of Eta are evaluated so valid
C       range is all but first row and column in overlaps.
        iMin = 1-OLx+1
        iMax = sNx+OLx
        jMin = 1-OLy+1
        jMax = sNy+OLy

C-      Calculate gradient of surface Potentiel
        CALL CALC_GRAD_PHI_SURF(
     I       bi,bj,iMin,iMax,jMin,jMax,
     I       etaN,
     O       phiSurfX,phiSurfY,
     I       myThid )

        IF (momStepping) THEN
C--       Loop over all layers, top to bottom
          DO k=1,Nr
C-        Update velocity fields:  V(n) = V** - dt * grad Eta
            CALL CORRECTION_STEP(
     I         bi,bj,iMin,iMax,jMin,jMax,k,
     I         phiSurfX,phiSurfY,myTime,myThid )

#ifdef    ALLOW_OBCS
c         IF (useOBCS) THEN
c           CALL OBCS_APPLY_UV(bi,bj,k,uVel,vVel,myThid)
c         ENDIF
#endif    /* ALLOW_OBCS */
C--       End DO k=1,Nr
          ENDDO
        ENDIF

#ifdef ALLOW_OBCS
        IF (useOBCS) THEN
          CALL OBCS_APPLY_UV( bi, bj, 0, uVel, vVel, myThid )
        ENDIF
#endif /* ALLOW_OBCS */

C--    End of 1rst bi,bj loop
       ENDDO
      ENDDO

C--- 2nd Part : Adjustment.

C--   Filter (and exchange)
#ifdef ALLOW_SHAP_FILT
      IF (useSHAP_FILT) THEN
       IF ( .NOT.shap_filt_uvStar ) THEN
        CALL TIMER_START('SHAP_FILT_UV       [MOM_CORR_STEP]',myThid)
        CALL SHAP_FILT_APPLY_UV( uVel, vVel, myTime, myIter, myThid )
        CALL TIMER_STOP ('SHAP_FILT_UV       [MOM_CORR_STEP]',myThid)
       ENDIF
      ENDIF
#endif
#ifdef ALLOW_ZONAL_FILT
      IF (useZONAL_FILT) THEN
       IF ( .NOT.zonal_filt_uvStar ) THEN
        CALL TIMER_START('ZONAL_FILT_UV      [MOM_CORR_STEP]',myThid)
        CALL ZONAL_FILT_APPLY_UV( uVel, vVel, myThid )
        CALL TIMER_STOP ('ZONAL_FILT_UV      [MOM_CORR_STEP]',myThid)
       ENDIF
      ENDIF
#endif

#ifdef ALLOW_AUTODIFF_TAMC
# ifdef NONLIN_FRSURF
CADJ STORE uvel, vvel  = comlev1, key = ikey_dynamics, byte = isbyte
# endif
#endif
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)

#ifdef ALLOW_AUTODIFF_TAMC
          act1 = bi - myBxLo(myThid)
          max1 = myBxHi(myThid) - myBxLo(myThid) + 1
          act2 = bj - myByLo(myThid)
          max2 = myByHi(myThid) - myByLo(myThid) + 1
          act3 = myThid - 1
          max3 = nTx*nTy
          act4 = ikey_dynamics - 1
          idynkey = (act1 + 1) + act2*max1
     &                      + act3*max1*max2
     &                      + act4*max1*max2*max3
# ifdef NONLIN_FRSURF
#  ifndef DISABLE_RSTAR_CODE
cph-test
CADJ STORE detahdt(:,:,bi,bj)    = comlev1_bibj, key = idynkey, byte = isbyte
CADJ STORE etan(:,:,bi,bj)       = comlev1_bibj, key = idynkey, byte = isbyte
CADJ STORE rstardhcdt(:,:,bi,bj) = comlev1_bibj, key = idynkey, byte = isbyte
#  endif
# endif
#endif

C--     Integrate continuity vertically
C--     for vertical velocity and "etaN" (exact volume conservation) :
        CALL INTEGR_CONTINUITY( bi, bj, uVel, vVel,
     I                          myTime, myIter, myThid )

C--    End of 2nd bi,bj loop
       ENDDO
      ENDDO

      IF ( exactConserv .AND. implicDiv2Dflow .NE. 0. _d 0)
     & _EXCH_XY_RL( etaN , myThid )
      IF ( implicitIntGravWave )
     & _EXCH_XYZ_RL( wVel , myThid )

      RETURN
      END
1	jmc	1.8	C $Header: /u/gcmpack/MITgcm/model/src/momentum_correction_step.F,v 1.7 2009/04/28 18:01:14 jmc Exp $
2	jmc	1.1	C $Name: $
3
4			#include "PACKAGES_CONFIG.h"
5			#include "CPP_OPTIONS.h"
6
7			CBOP
8			C !ROUTINE: MOMENTUM_CORRECTION_STEP
9			C !INTERFACE:
10			SUBROUTINE MOMENTUM_CORRECTION_STEP(myTime, myIter, myThid)
11			C !DESCRIPTION: \bv
12			C ==========================================================
13	jmc	1.2	C \| SUBROUTINE MOMENTUM_CORRECTION_STEP
14	jmc	1.1	C ==========================================================
15			C \|1rst Part : Update U,V.
16			C \|
17			C \| The arrays used for time stepping are cycled.
18			C \| Momentum:
19			C \| V(n) = Gv(n) - dt * grad Eta
20			C \|
21			C \|part1: update U,V
22	jmc	1.8	C \| U,V (contained in gU,gV) have the surface
23			C \| pressure gradient term added and the result stored
24			C \| in U,V (contained in uVel, vVel)
25			C \|
26	jmc	1.1	C \|part2: Adjustments
27	jmc	1.8	C \| o Filter U,V (Shapiro Filter, Zonal_Filter)
28			C \| o Compute again Eta (exact volume conservation)
29	jmc	1.1	C \| o Compute vertical velocity
30			C ==========================================================
31			C \ev
32
33			C !USES:
34			IMPLICIT NONE
35			C == Global variables ===
36			#include "SIZE.h"
37			#include "EEPARAMS.h"
38			#include "PARAMS.h"
39			#include "DYNVARS.h"
40
41			#ifdef ALLOW_SHAP_FILT
42			#include "SHAP_FILT.h"
43			#endif
44	jmc	1.2	#ifdef ALLOW_ZONAL_FILT
45	jmc	1.1	#include "ZONAL_FILT.h"
46			#endif
47	heimbach	1.3	#ifdef ALLOW_AUTODIFF_TAMC
48			#include "tamc.h"
49	heimbach	1.5	#include "tamc_keys.h"
50			# ifdef NONLIN_FRSURF
51			# include "SURFACE.h"
52			# endif
53	heimbach	1.3	#endif
54	jmc	1.1
55			C !INPUT/OUTPUT PARAMETERS:
56			C == Routine arguments ==
57	jmc	1.8	C myTime :: Current time in simulation
58			C myIter :: Current iteration number in simulation
59			C myThid :: Thread number for this instance of the routine.
60	jmc	1.1	_RL myTime
61			INTEGER myIter
62			INTEGER myThid
63
64			C !LOCAL VARIABLES:
65			C == Local variables
66			_RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
67			_RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
68			INTEGER iMin,iMax
69			INTEGER jMin,jMax
70			INTEGER bi,bj
71			INTEGER k,i,j
72			CEOP
73
74			DO bj=myByLo(myThid),myByHi(myThid)
75			DO bi=myBxLo(myThid),myBxHi(myThid)
76
77			C-- Set up work arrays that need valid initial values
78			DO j=1-OLy,sNy+OLy
79			DO i=1-OLx,sNx+OLx
80			phiSurfX(i,j)=0.
81			phiSurfY(i,j)=0.
82			ENDDO
83			ENDDO
84
85			C Loop range: Gradients of Eta are evaluated so valid
86			C range is all but first row and column in overlaps.
87			iMin = 1-OLx+1
88			iMax = sNx+OLx
89			jMin = 1-OLy+1
90			jMax = sNy+OLy
91
92			C- Calculate gradient of surface Potentiel
93			CALL CALC_GRAD_PHI_SURF(
94			I bi,bj,iMin,iMax,jMin,jMax,
95			I etaN,
96			O phiSurfX,phiSurfY,
97			I myThid )
98
99	jmc	1.8	IF (momStepping) THEN
100			C-- Loop over all layers, top to bottom
101			DO k=1,Nr
102	jmc	1.1	C- Update velocity fields: V(n) = V** - dt * grad Eta
103	jmc	1.8	CALL CORRECTION_STEP(
104			I bi,bj,iMin,iMax,jMin,jMax,k,
105	jmc	1.1	I phiSurfX,phiSurfY,myTime,myThid )
106
107			#ifdef ALLOW_OBCS
108	jmc	1.8	c IF (useOBCS) THEN
109			c CALL OBCS_APPLY_UV(bi,bj,k,uVel,vVel,myThid)
110			c ENDIF
111	jmc	1.1	#endif /* ALLOW_OBCS */
112	jmc	1.8	C-- End DO k=1,Nr
113			ENDDO
114			ENDIF
115
116			#ifdef ALLOW_OBCS
117			IF (useOBCS) THEN
118			CALL OBCS_APPLY_UV( bi, bj, 0, uVel, vVel, myThid )
119			ENDIF
120			#endif /* ALLOW_OBCS */
121	jmc	1.1
122			C-- End of 1rst bi,bj loop
123			ENDDO
124			ENDDO
125
126			C--- 2nd Part : Adjustment.
127
128			C-- Filter (and exchange)
129			#ifdef ALLOW_SHAP_FILT
130			IF (useSHAP_FILT) THEN
131	jmc	1.2	IF ( .NOT.shap_filt_uvStar ) THEN
132			CALL TIMER_START('SHAP_FILT_UV [MOM_CORR_STEP]',myThid)
133			CALL SHAP_FILT_APPLY_UV( uVel, vVel, myTime, myIter, myThid )
134			CALL TIMER_STOP ('SHAP_FILT_UV [MOM_CORR_STEP]',myThid)
135			ENDIF
136	jmc	1.1	ENDIF
137	jmc	1.2	#endif
138	jmc	1.1	#ifdef ALLOW_ZONAL_FILT
139			IF (useZONAL_FILT) THEN
140	jmc	1.2	IF ( .NOT.zonal_filt_uvStar ) THEN
141			CALL TIMER_START('ZONAL_FILT_UV [MOM_CORR_STEP]',myThid)
142			CALL ZONAL_FILT_APPLY_UV( uVel, vVel, myThid )
143			CALL TIMER_STOP ('ZONAL_FILT_UV [MOM_CORR_STEP]',myThid)
144			ENDIF
145	jmc	1.1	ENDIF
146	jmc	1.2	#endif
147	jmc	1.1
148	heimbach	1.3	#ifdef ALLOW_AUTODIFF_TAMC
149			# ifdef NONLIN_FRSURF
150	heimbach	1.6	CADJ STORE uvel, vvel = comlev1, key = ikey_dynamics, byte = isbyte
151	heimbach	1.3	# endif
152			#endif
153	jmc	1.1	DO bj=myByLo(myThid),myByHi(myThid)
154			DO bi=myBxLo(myThid),myBxHi(myThid)
155
156	heimbach	1.5	#ifdef ALLOW_AUTODIFF_TAMC
157			act1 = bi - myBxLo(myThid)
158			max1 = myBxHi(myThid) - myBxLo(myThid) + 1
159			act2 = bj - myByLo(myThid)
160			max2 = myByHi(myThid) - myByLo(myThid) + 1
161			act3 = myThid - 1
162			max3 = nTx*nTy
163			act4 = ikey_dynamics - 1
164			idynkey = (act1 + 1) + act2*max1
165			& + act3max1max2
166			& + act4max1max2*max3
167			# ifdef NONLIN_FRSURF
168			# ifndef DISABLE_RSTAR_CODE
169			cph-test
170			CADJ STORE detahdt(:,:,bi,bj) = comlev1_bibj, key = idynkey, byte = isbyte
171			CADJ STORE etan(:,:,bi,bj) = comlev1_bibj, key = idynkey, byte = isbyte
172			CADJ STORE rstardhcdt(:,:,bi,bj) = comlev1_bibj, key = idynkey, byte = isbyte
173			# endif
174			# endif
175			#endif
176
177	jmc	1.1	C-- Integrate continuity vertically
178			C-- for vertical velocity and "etaN" (exact volume conservation) :
179			CALL INTEGR_CONTINUITY( bi, bj, uVel, vVel,
180			I myTime, myIter, myThid )
181
182			C-- End of 2nd bi,bj loop
183			ENDDO
184			ENDDO
185
186			IF ( exactConserv .AND. implicDiv2Dflow .NE. 0. _d 0)
187	jmc	1.7	& _EXCH_XY_RL( etaN , myThid )
188	jmc	1.8	IF ( implicitIntGravWave )
189	jmc	1.7	& _EXCH_XYZ_RL( wVel , myThid )
190	jmc	1.1
191			RETURN
192			END