model/src/the_correction_step.F

C $Header: /u/gcmpack/MITgcm/model/src/the_correction_step.F,v 1.20 2002/11/01 22:00:33 mlosch 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)
C     | Momentum:
C     |           V(n) = Gv(n) - dt * grad Eta
C     |
C     |part1: update U,V,T,S                                      
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     |  T* (contained in gT) is copied to T (theta)           
C     |  S* (contained in gS) is copied to S (salt)            
C     |                                                           
C     |part2: Adjustments & Diagnostics                                        
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 Compute vertical velocity
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_SHAP_FILT
#include "SHAP_FILT.h"
#endif
#ifdef ALLOW_ZONAL_FILT   
#include "ZONAL_FILT.h"
#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 )

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

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_PTRACERS
C-        Update passive tracer fields:  T(n) = T**, Gt(n-1) = Gt(n)
          IF (usePTRACERS)
     &        CALL PTRACERS_CYCLE(bi,bj,k,myIter,myTime,myThid)
#endif /* ALLOW_PTRACERS */


#ifdef    ALLOW_OBCS
          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.

C--   Filter (and exchange)
#ifdef ALLOW_SHAP_FILT
      IF (useSHAP_FILT) THEN
       IF ( .NOT.shap_filt_uvStar )
     &  CALL SHAP_FILT_APPLY_UV( uVel, vVel, myTime, myIter, myThid )

       IF ( .NOT.(staggerTimeStep .AND. shap_filt_TrStagg) )
     &  CALL SHAP_FILT_APPLY_TS( theta,salt, myTime, myIter, myThid )
      ENDIF
#endif 
#ifdef ALLOW_ZONAL_FILT
      IF (useZONAL_FILT) THEN
       IF ( .NOT.zonal_filt_uvStar )
     &  CALL ZONAL_FILT_APPLY_UV( uVel, vVel, myThid )

       IF ( .NOT.(staggerTimeStep .AND. zonal_filt_TrStagg) )
     &  CALL ZONAL_FILT_APPLY_TS( theta, salt, 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
        iMax = sNx
        jMin = 1
        jMax = sNy
        CALL CONVECTIVE_ADJUSTMENT(
     I       bi, bj, iMin, iMax, jMin, jMax,
     I       myTime, myIter, myThid )

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

#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

      RETURN
      END
1	C $Header: /u/gcmpack/MITgcm/model/src/the_correction_step.F,v 1.20 2002/11/01 22:00:33 mlosch Exp $
2	C Tag $Name: $
3
4	#include "CPP_OPTIONS.h"
5
6	CBOP
7	C !ROUTINE: THE_CORRECTION_STEP
8	C !INTERFACE:
9	SUBROUTINE THE_CORRECTION_STEP(myTime, myIter, myThid)
10	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)
19	C \| Momentum:
20	C \| V(n) = Gv(n) - dt * grad Eta
21	C \|
22	C \|part1: update U,V,T,S
23	C \| U,V (contained in gU,gV) have the surface
24	C \| pressure gradient term added and the result stored
25	C \| in U,V (contained in uVel, vVel)
26	C \| T* (contained in gT) is copied to T (theta)
27	C \| S* (contained in gS) is copied to S (salt)
28	C \|
29	C \|part2: Adjustments & Diagnostics
30	C \| o Filter U,V,T,S (Shapiro Filter, Zonal_Filter)
31	C \| o Convective Adjustment
32	C \| o Compute again Eta (exact volume conservation)
33	C \| o Compute vertical velocity
34	C \| o Diagmnostic of state variables (Time average)
35	C ==========================================================
36	C \ev
37
38	C !USES:
39	IMPLICIT NONE
40	C == Global variables ===
41	#include "SIZE.h"
42	#include "EEPARAMS.h"
43	#include "PARAMS.h"
44	#include "DYNVARS.h"
45	#ifdef ALLOW_PASSIVE_TRACER
46	#include "TR1.h"
47	#endif
48
49	#ifdef ALLOW_SHAP_FILT
50	#include "SHAP_FILT.h"
51	#endif
52	#ifdef ALLOW_ZONAL_FILT
53	#include "ZONAL_FILT.h"
54	#endif
55
56	C !INPUT/OUTPUT PARAMETERS:
57	C == Routine arguments ==
58	C myTime - Current time in simulation
59	C myIter - Current iteration number in simulation
60	C myThid - Thread number for this instance of the routine.
61	_RL myTime
62	INTEGER myIter
63	INTEGER myThid
64
65	C !LOCAL VARIABLES:
66	C == Local variables
67	_RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
68	_RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
69	INTEGER iMin,iMax
70	INTEGER jMin,jMax
71	INTEGER bi,bj
72	INTEGER k,i,j
73
74	CEOP
75
76	DO bj=myByLo(myThid),myByHi(myThid)
77	DO bi=myBxLo(myThid),myBxHi(myThid)
78
79	C-- Set up work arrays that need valid initial values
80	DO j=1-OLy,sNy+OLy
81	DO i=1-OLx,sNx+OLx
82	phiSurfX(i,j)=0.
83	phiSurfY(i,j)=0.
84	ENDDO
85	ENDDO
86
87	C Loop range: Gradients of Eta are evaluated so valid
88	C range is all but first row and column in overlaps.
89	iMin = 1-OLx+1
90	iMax = sNx+OLx
91	jMin = 1-OLy+1
92	jMax = sNy+OLy
93
94	C- Calculate gradient of surface Potentiel
95	CALL CALC_GRAD_PHI_SURF(
96	I bi,bj,iMin,iMax,jMin,jMax,
97	I etaN,
98	O phiSurfX,phiSurfY,
99	I myThid )
100
101	C-- Loop over all layers, top to bottom
102	DO K=1,Nr
103
104	C- Update velocity fields: V(n) = V** - dt * grad Eta
105	IF (momStepping)
106	& CALL CORRECTION_STEP(
107	I bi,bj,iMin,iMax,jMin,jMax,K,
108	I phiSurfX,phiSurfY,myTime,myThid )
109
110	C- Update tracer fields: T(n) = T**, Gt(n-1) = Gt(n)
111	IF (tempStepping)
112	& CALL CYCLE_TRACER(
113	I bi,bj,iMin,iMax,jMin,jMax,K,
114	U theta,gT,gTNm1,
115	I myTime,myThid )
116	IF (saltStepping)
117	& CALL CYCLE_TRACER(
118	I bi,bj,iMin,iMax,jMin,jMax,K,
119	U salt,gS,gSNm1,
120	I myTime,myThid )
121	#ifdef ALLOW_PASSIVE_TRACER
122	IF (tr1Stepping)
123	& CALL CYCLE_TRACER(
124	I bi,bj,iMin,iMax,jMin,jMax,K,
125	U Tr1,gTr1,gTr1Nm1,
126	I myTime,myThid )
127	#endif
128	#ifdef ALLOW_PTRACERS
129	C- Update passive tracer fields: T(n) = T**, Gt(n-1) = Gt(n)
130	IF (usePTRACERS)
131	& CALL PTRACERS_CYCLE(bi,bj,k,myIter,myTime,myThid)
132	#endif /* ALLOW_PTRACERS */
133
134
135	#ifdef ALLOW_OBCS
136	IF (useOBCS) THEN
137	CALL OBCS_APPLY_UV(bi,bj,K,uVel,vVel,myThid)
138	ENDIF
139	#endif /* ALLOW_OBCS */
140
141	C-- End DO K=1,Nr
142	ENDDO
143
144	C-- End of 1rst bi,bj loop
145	ENDDO
146	ENDDO
147
148	C--- 2nd Part : Adjustment.
149	C
150	C Static stability is calculated and the tracers are
151	C convective adjusted where statically unstable.
152
153	C-- Filter (and exchange)
154	#ifdef ALLOW_SHAP_FILT
155	IF (useSHAP_FILT) THEN
156	IF ( .NOT.shap_filt_uvStar )
157	& CALL SHAP_FILT_APPLY_UV( uVel, vVel, myTime, myIter, myThid )
158
159	IF ( .NOT.(staggerTimeStep .AND. shap_filt_TrStagg) )
160	& CALL SHAP_FILT_APPLY_TS( theta,salt, myTime, myIter, myThid )
161	ENDIF
162	#endif
163	#ifdef ALLOW_ZONAL_FILT
164	IF (useZONAL_FILT) THEN
165	IF ( .NOT.zonal_filt_uvStar )
166	& CALL ZONAL_FILT_APPLY_UV( uVel, vVel, myThid )
167
168	IF ( .NOT.(staggerTimeStep .AND. zonal_filt_TrStagg) )
169	& CALL ZONAL_FILT_APPLY_TS( theta, salt, myThid )
170
171	ENDIF
172	#endif
173
174	DO bj=myByLo(myThid),myByHi(myThid)
175	DO bi=myBxLo(myThid),myBxHi(myThid)
176
177	C-- Convectively adjust new fields to be statically stable
178	iMin = 1
179	iMax = sNx
180	jMin = 1
181	jMax = sNy
182	CALL CONVECTIVE_ADJUSTMENT(
183	I bi, bj, iMin, iMax, jMin, jMax,
184	I myTime, myIter, myThid )
185
186	C-- Integrate continuity vertically
187	C-- for vertical velocity and "etaN" (exact volume conservation) :
188	CALL INTEGR_CONTINUITY( bi, bj, uVel, vVel,
189	I myTime, myIter, myThid )
190
191	#ifdef ALLOW_TIMEAVE
192	IF (taveFreq.GT.0.) THEN
193	CALL TIMEAVE_STATVARS(myTime, myIter, bi, bj, myThid)
194	ENDIF
195	#endif /* ALLOW_TIMEAVE */
196
197	C-- End of 2nd bi,bj loop
198	ENDDO
199	ENDDO
200
201	RETURN
202	END