model/src/the_correction_step.F

C $Header: /u/gcmpack/models/MITgcmUV/model/src/the_correction_step.F,v 1.13 2001/09/19 13:58:08 jmc Exp $
C Tag $Name:  $

#include "CPP_OPTIONS.h"

CBOP
C     !ROUTINE: THE_CORRECTION_STEP
C     !INTERFACE:
      SUBROUTINE THE_CORRECTION_STEP(myTime, myIter, myThid)
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE THE_CORRECTION_STEP                            
C     *==========================================================*
C     |1rst Part : Update U,V,T,S.
C     |
C     | The arrays used for time stepping are cycled.
C     | Tracers:
C     |           T(n) = Gt(n-1)
C     |           Gt(n-1) = Gt(n)
C     | Momentum:
C     |           V(n) = Gv(n-1) - dt * grad Eta
C     |           Gv(n-1) = Gv(n)
C     |
C     |part1: update U,V,T,S                                      
C     |  U*,V* (contained in gUnm1,gVnm1) have the surface        
C     |     pressure gradient term added and the result stored    
C     |     in U,V (contained in uVel, vVel)                      
C     |  T* (contained in gTnm1) is copied to T (theta)           
C     |  S* (contained in gSnm1) is copied to S (salt)            
C     |                                                           
C     |part2: Adjustments.                                        
C     |   o Filter  U,V,T,S (Shapiro Filter, Zonal_Filter)        
C     |   o Convective Adjustment                                 
C     |   o Compute again Eta (exact volume conservation)         
C     |   o Diagmnostic of state variables (Time average)         
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
C     == Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "DYNVARS.h"
#ifdef ALLOW_PASSIVE_TRACER
#include "TR1.h"
#endif
#ifdef ALLOW_AUTODIFF_TAMC
#include "tamc.h"
#include "tamc_keys.h"
#endif /* ALLOW_AUTODIFF_TAMC */

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)

#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

          ikey = (act1 + 1) + act2*max1
     &                      + act3*max1*max2
     &                      + act4*max1*max2*max3
#endif /* ALLOW_AUTODIFF_TAMC */

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 )

C--     Loop over all layers, top to bottom
        DO K=1,Nr

#ifdef ALLOW_AUTODIFF_TAMC
          kkey = (ikey-1)*Nr + k
#endif

C-        Update velocity fields:  V(n) = V** - dt * grad Eta
          IF (momStepping)
     &      CALL CORRECTION_STEP(
     I         bi,bj,iMin,iMax,jMin,jMax,K,
     I         phiSurfX,phiSurfY,myTime,myThid )

C-        Update tracer fields:  T(n) = T**, Gt(n-1) = Gt(n)
          IF (tempStepping)
     &      CALL CYCLE_TRACER(
     I           bi,bj,iMin,iMax,jMin,jMax,K,
     U           theta,gT,gTNm1,
     I           myTime,myThid )
          IF (saltStepping)
     &      CALL CYCLE_TRACER(
     I           bi,bj,iMin,iMax,jMin,jMax,K,
     U           salt,gS,gSNm1,
     I           myTime,myThid )
#ifdef ALLOW_PASSIVE_TRACER
          IF (tr1Stepping)
     &      CALL CYCLE_TRACER(
     I           bi,bj,iMin,iMax,jMin,jMax,K,
     U           Tr1,gTr1,gTr1Nm1,
     I           myTime,myThid )
#endif

#ifdef    ALLOW_OBCS
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE uvel (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
CADJ STORE vvel (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
          IF (useOBCS) THEN
            CALL OBCS_APPLY_UV(bi,bj,K,uVel,vVel,myThid)
          ENDIF
#endif    /* ALLOW_OBCS */

C--     End DO K=1,Nr
        ENDDO

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

C--- 2nd Part : Adjustment.
C
C       Static stability is calculated and the tracers are
C       convective adjusted where statically unstable.

#ifdef ALLOW_SHAP_FILT
      IF (useSHAP_FILT) THEN
C--   Filter (and exchange).
      CALL SHAP_FILT_APPLY(
     I                     uVel, vVel, theta, salt,
     I                     myTime, myIter, myThid )
      ENDIF
#endif 

#ifdef ALLOW_ZONAL_FILT
      IF (zonal_filt_lat.LT.90.) THEN
      CALL ZONAL_FILT_APPLY(
     U                      uVel, vVel, theta, salt,
     I                      myThid )
      ENDIF
#endif 

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

C--     Convectively adjust new fields to be statically stable
        iMin = 1-OLx+1
        iMax = sNx+OLx
        jMin = 1-OLy+1
        jMax = sNy+OLy
        CALL CONVECTIVE_ADJUSTMENT(
     I       bi, bj, iMin, iMax, jMin, jMax,
     I       myTime, myIter, myThid )

#ifdef EXACT_CONSERV
        IF (exactConserv) THEN
C--     Compute again "eta" to satisfy exactly the total Volume Conservation :
          CALL CALC_EXACT_ETA( .TRUE., bi,bj, uVel,vVel,
     I                         myTime, myIter, myThid )
        ENDIF
#endif /* EXACT_CONSERV */ 

#ifdef ALLOW_TIMEAVE
        IF (taveFreq.GT.0.) THEN 
          CALL TIMEAVE_STATVARS(myTime, myIter, bi, bj, myThid)
        ENDIF
#endif /* ALLOW_TIMEAVE */

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

#ifdef EXACT_CONSERV
       IF (exactConserv .AND. implicDiv2Dflow .NE. 0. _d 0) 
     &    _EXCH_XY_R8(etaN, myThid ) 
#endif /* EXACT_CONSERV */ 

      RETURN
      END
1	cnh	1.14	C $Header: /u/gcmpack/models/MITgcmUV/model/src/the_correction_step.F,v 1.13 2001/09/19 13:58:08 jmc Exp $
2	jmc	1.13	C Tag $Name: $
3	adcroft	1.2
4			#include "CPP_OPTIONS.h"
5
6	cnh	1.14	CBOP
7			C !ROUTINE: THE_CORRECTION_STEP
8			C !INTERFACE:
9	adcroft	1.2	SUBROUTINE THE_CORRECTION_STEP(myTime, myIter, myThid)
10	cnh	1.14	C !DESCRIPTION: \bv
11			C ==========================================================
12			C \| SUBROUTINE THE_CORRECTION_STEP
13			C ==========================================================
14			C \|1rst Part : Update U,V,T,S.
15			C \|
16			C \| The arrays used for time stepping are cycled.
17			C \| Tracers:
18			C \| T(n) = Gt(n-1)
19			C \| Gt(n-1) = Gt(n)
20			C \| Momentum:
21			C \| V(n) = Gv(n-1) - dt * grad Eta
22			C \| Gv(n-1) = Gv(n)
23			C \|
24			C \|part1: update U,V,T,S
25			C \| U,V (contained in gUnm1,gVnm1) have the surface
26			C \| pressure gradient term added and the result stored
27			C \| in U,V (contained in uVel, vVel)
28			C \| T* (contained in gTnm1) is copied to T (theta)
29			C \| S* (contained in gSnm1) is copied to S (salt)
30			C \|
31			C \|part2: Adjustments.
32			C \| o Filter U,V,T,S (Shapiro Filter, Zonal_Filter)
33			C \| o Convective Adjustment
34			C \| o Compute again Eta (exact volume conservation)
35			C \| o Diagmnostic of state variables (Time average)
36			C ==========================================================
37			C \ev
38
39			C !USES:
40	adcroft	1.2	IMPLICIT NONE
41			C == Global variables ===
42			#include "SIZE.h"
43			#include "EEPARAMS.h"
44			#include "PARAMS.h"
45			#include "DYNVARS.h"
46	heimbach	1.12	#ifdef ALLOW_PASSIVE_TRACER
47	heimbach	1.11	#include "TR1.h"
48	heimbach	1.12	#endif
49	heimbach	1.7	#ifdef ALLOW_AUTODIFF_TAMC
50			#include "tamc.h"
51			#include "tamc_keys.h"
52			#endif /* ALLOW_AUTODIFF_TAMC */
53
54	cnh	1.14	C !INPUT/OUTPUT PARAMETERS:
55	adcroft	1.2	C == Routine arguments ==
56			C myTime - Current time in simulation
57			C myIter - Current iteration number in simulation
58			C myThid - Thread number for this instance of the routine.
59			_RL myTime
60			INTEGER myIter
61			INTEGER myThid
62
63	cnh	1.14	C !LOCAL VARIABLES:
64	adcroft	1.2	C == Local variables
65	jmc	1.6	_RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
66			_RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
67	adcroft	1.2	INTEGER iMin,iMax
68			INTEGER jMin,jMax
69			INTEGER bi,bj
70			INTEGER k,i,j
71
72	cnh	1.14	CEOP
73	adcroft	1.9
74	adcroft	1.2	DO bj=myByLo(myThid),myByHi(myThid)
75			DO bi=myBxLo(myThid),myBxHi(myThid)
76	heimbach	1.7
77			#ifdef ALLOW_AUTODIFF_TAMC
78			act1 = bi - myBxLo(myThid)
79			max1 = myBxHi(myThid) - myBxLo(myThid) + 1
80
81			act2 = bj - myByLo(myThid)
82			max2 = myByHi(myThid) - myByLo(myThid) + 1
83
84			act3 = myThid - 1
85			max3 = nTx*nTy
86
87			act4 = ikey_dynamics - 1
88
89			ikey = (act1 + 1) + act2*max1
90			& + act3max1max2
91			& + act4max1max2*max3
92			#endif /* ALLOW_AUTODIFF_TAMC */
93	adcroft	1.2
94			C-- Set up work arrays that need valid initial values
95			DO j=1-OLy,sNy+OLy
96			DO i=1-OLx,sNx+OLx
97	jmc	1.6	phiSurfX(i,j)=0.
98			phiSurfY(i,j)=0.
99	adcroft	1.2	ENDDO
100			ENDDO
101
102			C Loop range: Gradients of Eta are evaluated so valid
103			C range is all but first row and column in overlaps.
104			iMin = 1-OLx+1
105			iMax = sNx+OLx
106			jMin = 1-OLy+1
107			jMax = sNy+OLy
108
109	jmc	1.6	C- Calculate gradient of surface Potentiel
110			CALL CALC_GRAD_PHI_SURF(
111	adcroft	1.2	I bi,bj,iMin,iMax,jMin,jMax,
112	jmc	1.5	I etaN,
113	jmc	1.6	O phiSurfX,phiSurfY,
114	adcroft	1.2	I myThid )
115
116			C-- Loop over all layers, top to bottom
117			DO K=1,Nr
118
119	heimbach	1.8	#ifdef ALLOW_AUTODIFF_TAMC
120			kkey = (ikey-1)*Nr + k
121			#endif
122
123	adcroft	1.2	C- Update velocity fields: V(n) = V** - dt * grad Eta
124	adcroft	1.3	IF (momStepping)
125			& CALL CORRECTION_STEP(
126	adcroft	1.2	I bi,bj,iMin,iMax,jMin,jMax,K,
127	jmc	1.6	I phiSurfX,phiSurfY,myTime,myThid )
128	adcroft	1.2
129			C- Update tracer fields: T(n) = T**, Gt(n-1) = Gt(n)
130			IF (tempStepping)
131			& CALL CYCLE_TRACER(
132			I bi,bj,iMin,iMax,jMin,jMax,K,
133			U theta,gT,gTNm1,
134			I myTime,myThid )
135			IF (saltStepping)
136			& CALL CYCLE_TRACER(
137			I bi,bj,iMin,iMax,jMin,jMax,K,
138			U salt,gS,gSNm1,
139	heimbach	1.11	I myTime,myThid )
140	heimbach	1.12	#ifdef ALLOW_PASSIVE_TRACER
141	heimbach	1.11	IF (tr1Stepping)
142			& CALL CYCLE_TRACER(
143			I bi,bj,iMin,iMax,jMin,jMax,K,
144			U Tr1,gTr1,gTr1Nm1,
145	adcroft	1.2	I myTime,myThid )
146	heimbach	1.12	#endif
147	adcroft	1.2
148			#ifdef ALLOW_OBCS
149			#ifdef ALLOW_AUTODIFF_TAMC
150	heimbach	1.8	CADJ STORE uvel (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
151			CADJ STORE vvel (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
152			CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
153			CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
154	adcroft	1.2	#endif /* ALLOW_AUTODIFF_TAMC */
155			IF (useOBCS) THEN
156			CALL OBCS_APPLY_UV(bi,bj,K,uVel,vVel,myThid)
157			ENDIF
158			#endif /* ALLOW_OBCS */
159
160			C-- End DO K=1,Nr
161			ENDDO
162
163	adcroft	1.9	C-- End of 1rst bi,bj loop
164			ENDDO
165			ENDDO
166
167			C--- 2nd Part : Adjustment.
168			C
169			C Static stability is calculated and the tracers are
170			C convective adjusted where statically unstable.
171
172			#ifdef ALLOW_SHAP_FILT
173			IF (useSHAP_FILT) THEN
174			C-- Filter (and exchange).
175			CALL SHAP_FILT_APPLY(
176			I uVel, vVel, theta, salt,
177			I myTime, myIter, myThid )
178			ENDIF
179			#endif
180
181			#ifdef ALLOW_ZONAL_FILT
182			IF (zonal_filt_lat.LT.90.) THEN
183			CALL ZONAL_FILT_APPLY(
184			U uVel, vVel, theta, salt,
185			I myThid )
186			ENDIF
187			#endif
188
189			DO bj=myByLo(myThid),myByHi(myThid)
190			DO bi=myBxLo(myThid),myBxHi(myThid)
191
192	adcroft	1.2	C-- Convectively adjust new fields to be statically stable
193	adcroft	1.9	iMin = 1-OLx+1
194			iMax = sNx+OLx
195			jMin = 1-OLy+1
196			jMax = sNy+OLy
197	adcroft	1.2	CALL CONVECTIVE_ADJUSTMENT(
198			I bi, bj, iMin, iMax, jMin, jMax,
199			I myTime, myIter, myThid )
200	jmc	1.4
201	adcroft	1.9	#ifdef EXACT_CONSERV
202			IF (exactConserv) THEN
203			C-- Compute again "eta" to satisfy exactly the total Volume Conservation :
204	jmc	1.13	CALL CALC_EXACT_ETA( .TRUE., bi,bj, uVel,vVel,
205	adcroft	1.9	I myTime, myIter, myThid )
206			ENDIF
207			#endif /* EXACT_CONSERV */
208
209	jmc	1.5	#ifdef ALLOW_TIMEAVE
210	jmc	1.4	IF (taveFreq.GT.0.) THEN
211	jmc	1.5	CALL TIMEAVE_STATVARS(myTime, myIter, bi, bj, myThid)
212	jmc	1.4	ENDIF
213	jmc	1.5	#endif /* ALLOW_TIMEAVE */
214	adcroft	1.2
215	adcroft	1.9	C-- End of 2nd bi,bj loop
216	adcroft	1.2	ENDDO
217			ENDDO
218	adcroft	1.9
219			#ifdef EXACT_CONSERV
220	jmc	1.13	IF (exactConserv .AND. implicDiv2Dflow .NE. 0. _d 0)
221			& _EXCH_XY_R8(etaN, myThid )
222	adcroft	1.9	#endif /* EXACT_CONSERV */
223	adcroft	1.2
224			RETURN
225			END