/[MITgcm]/MITgcm/model/src/timestep.F
ViewVC logotype

Diff of /MITgcm/model/src/timestep.F

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

revision 1.31 by jmc, Tue Feb 18 15:27:25 2003 UTC revision 1.41 by jmc, Fri Sep 30 00:22:37 2005 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4    #include "PACKAGES_CONFIG.h"
5  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
6    
7  CBOP  CBOP
8  C     !ROUTINE: TIMESTEP  C     !ROUTINE: TIMESTEP
9  C     !INTERFACE:  C     !INTERFACE:
10        SUBROUTINE TIMESTEP( bi, bj, iMin, iMax, jMin, jMax, K,        SUBROUTINE TIMESTEP( bi, bj, iMin, iMax, jMin, jMax, k,
11       I                     dPhiHydX,dPhiHydY, phiSurfX, phiSurfY,       I                     dPhiHydX,dPhiHydY, phiSurfX, phiSurfY,
12       I                     myIter, myThid )       I                     guDissip, gvDissip,
13         I                     myTime, myIter, myThid )
14  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
15  C     *==========================================================*  C     *==========================================================*
16  C     | S/R TIMESTEP                                                C     | S/R TIMESTEP                                              
# Line 31  C     == Routine Arguments == Line 33  C     == Routine Arguments ==
33  C     dPhiHydX,Y :: Gradient (X & Y directions) of Hydrostatic Potential  C     dPhiHydX,Y :: Gradient (X & Y directions) of Hydrostatic Potential
34  C     phiSurfX :: gradient of Surface potential (Pressure/rho, ocean)  C     phiSurfX :: gradient of Surface potential (Pressure/rho, ocean)
35  C     phiSurfY ::          or geopotential (atmos) in X and Y direction  C     phiSurfY ::          or geopotential (atmos) in X and Y direction
36    C     guDissip :: dissipation tendency (all explicit terms), u component
37    C     gvDissip :: dissipation tendency (all explicit terms), v component
38    
39        INTEGER bi,bj,iMin,iMax,jMin,jMax        INTEGER bi,bj,iMin,iMax,jMin,jMax
40        INTEGER K        INTEGER k
41        _RL     dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL     dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
42        _RL     dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL     dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
43        _RL     phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
44        _RL     phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
45          _RL     guDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
46          _RL     gvDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
47          _RL     myTime
48        INTEGER myIter, myThid        INTEGER myIter, myThid
49    
50  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
51  C     == Local variables ==  C     == Local variables ==
52          LOGICAL momForcing_In_AB
53          LOGICAL momDissip_In_AB
54        INTEGER i,j        INTEGER i,j
55        _RL ab15,ab05        _RL ab15,ab05
56        _RL phxFac,phyFac, psFac        _RL phxFac,phyFac, psFac
57        _RL     gUtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     gUtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
58        _RL     gVtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     gVtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
59    #ifdef ALLOW_DIAGNOSTICS
60    C    Allow diagnosis of external forcing
61          LOGICAL externForcDiagIsOn
62          LOGICAL  DIAGNOSTICS_IS_ON
63          EXTERNAL DIAGNOSTICS_IS_ON
64          _RL     gUext(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
65          _RL     gVext(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
66    #endif
67    #ifdef ALLOW_CD_CODE
68          _RL     guCor(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
69          _RL     gvCor(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
70    #endif
71  CEOP  CEOP
72    
73  C     Adams-Bashforth timestepping weights  C     Adams-Bashforth timestepping weights
# Line 57  C     Adams-Bashforth timestepping weigh Line 79  C     Adams-Bashforth timestepping weigh
79         ab05=-0.5-abeps         ab05=-0.5-abeps
80        ENDIF        ENDIF
81    
 C-- stagger time step: grad Phi_Hyp is not in gU,gV => add it in this S/R  
       IF (staggerTimeStep) THEN  
         phxFac = pfFacMom  
         phyFac = pfFacMom  
       ELSE  
         phxFac = 0.  
         phyFac = 0.  
       ENDIF  
   
82  C-- explicit part of the surface potential gradient is added in this S/R  C-- explicit part of the surface potential gradient is added in this S/R
83        psFac = pfFacMom*(1. _d 0 - implicSurfPress)        psFac = pfFacMom*(1. _d 0 - implicSurfPress)
84    
85    C--  factors for gradient (X & Y directions) of Hydrostatic Potential
86          phxFac = pfFacMom
87          phyFac = pfFacMom
88    
89    C-- including or excluding momentum forcing from Adams-Bashforth:
90          momForcing_In_AB = forcing_In_AB
91          momForcing_In_AB = .TRUE.
92          momDissip_In_AB  = .TRUE.
93    
94    #ifdef ALLOW_DIAGNOSTICS
95          externForcDiagIsOn = useDiagnostics .AND. momForcing
96          IF ( externForcDiagIsOn ) THEN
97            externForcDiagIsOn = DIAGNOSTICS_IS_ON('Um_Ext  ',myThid)
98         &                  .OR. DIAGNOSTICS_IS_ON('Vm_Ext  ',myThid)
99          ENDIF
100    #endif /* ALLOW_DIAGNOSTICS */
101    
102  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
103  C-    Compute effective gU term (including Adams-Bashforth weights) :  
104        DO j=jMin,jMax  C-    Initialize local arrays (not really necessary but safer)
105         DO i=iMin,iMax        DO j=1-Oly,sNy+Oly
106          gUtmp(i,j) = ab15*gU(i,j,k,bi,bj)         DO i=1-Olx,sNx+Olx
107       &             + ab05*gUNm1(i,j,k,bi,bj)            gUtmp(i,j) = 0. _d 0
108  #ifdef INCLUDE_CD_CODE          gVtmp(i,j) = 0. _d 0
109       &             + guCD(i,j,k,bi,bj)  #ifdef ALLOW_CD_CODE
110            guCor(i,j) = 0. _d 0
111            gvCor(i,j) = 0. _d 0
112  #endif  #endif
113         ENDDO         ENDDO
114        ENDDO        ENDDO
115    
116          IF ( .NOT.staggerTimeStep ) THEN
117    C--   Synchronous time step: add grad Phi_Hyp to gU,gV before doing Adams-Bashforth
118            DO j=jMin,jMax
119             DO i=iMin,iMax
120              gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) - phxFac*dPhiHydX(i,j)
121              gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) - phyFac*dPhiHydY(i,j)
122             ENDDO
123            ENDDO
124            phxFac = 0.
125            phyFac = 0.
126    c     ELSE
127    C--   Stagger time step: grad Phi_Hyp will be added later
128          ENDIF
129    
130    C--   Dissipation term inside the Adams-Bashforth:
131          IF ( momViscosity .AND. momDissip_In_AB) THEN
132            DO j=jMin,jMax
133             DO i=iMin,iMax
134              gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) + guDissip(i,j)
135              gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) + gvDissip(i,j)
136             ENDDO
137            ENDDO
138          ENDIF
139    
140    C--   Forcing term inside the Adams-Bashforth:
141          IF (momForcing .AND. momForcing_In_AB) THEN
142    #ifdef ALLOW_DIAGNOSTICS
143            IF ( useDiagnostics .AND. externForcDiagIsOn ) THEN
144             DO j=1,sNy+1
145              DO i=1,sNx+1
146               gUext(i,j) = gU(i,j,k,bi,bj)
147               gVext(i,j) = gV(i,j,k,bi,bj)
148              ENDDO
149             ENDDO
150            ENDIF
151    #endif /* ALLOW_DIAGNOSTICS */
152    
153            CALL EXTERNAL_FORCING_U(
154         I     iMin,iMax,jMin,jMax,bi,bj,k,
155         I     myTime,myThid)
156            CALL EXTERNAL_FORCING_V(
157         I     iMin,iMax,jMin,jMax,bi,bj,k,
158         I     myTime,myThid)
159    
160    #ifdef ALLOW_DIAGNOSTICS
161            IF ( useDiagnostics .AND. externForcDiagIsOn ) THEN
162             DO j=1,sNy+1
163              DO i=1,sNx+1
164               gUext(i,j) = gU(i,j,k,bi,bj)-gUext(i,j)
165               gVext(i,j) = gV(i,j,k,bi,bj)-gVext(i,j)
166              ENDDO
167             ENDDO
168            ENDIF
169    #endif /* ALLOW_DIAGNOSTICS */
170          ENDIF
171    
172          IF (useCDscheme) THEN
173    C-     for CD-scheme, store gU,Vtmp = gU,V^n + forcing
174            DO j=jMin,jMax
175             DO i=iMin,iMax
176              gUtmp(i,j) = gU(i,j,k,bi,bj)
177              gVtmp(i,j) = gV(i,j,k,bi,bj)
178             ENDDO
179            ENDDO
180          ENDIF
181    
182    C-    Compute effective gU,gV_[n+1/2] terms (including Adams-Bashforth weights)
183    C     and save gU,gV_[n] into guNm1,gvNm1 for the next time step.
184    #ifdef ALLOW_ADAMSBASHFORTH_3
185          CALL ADAMS_BASHFORTH3(
186         I                        bi, bj, k,
187         U                        gU, guNm,
188         I                        myIter, myThid )
189          CALL ADAMS_BASHFORTH3(
190         I                        bi, bj, k,
191         U                        gV, gvNm,
192         I                        myIter, myThid )
193    #else /* ALLOW_ADAMSBASHFORTH_3 */
194          CALL ADAMS_BASHFORTH2(
195         I                        bi, bj, k,
196         U                        gU, guNm1,
197         I                        myIter, myThid )
198          CALL ADAMS_BASHFORTH2(
199         I                        bi, bj, k,
200         U                        gV, gvNm1,
201         I                        myIter, myThid )
202    #endif /* ALLOW_ADAMSBASHFORTH_3 */
203                
204  #ifdef NONLIN_FRSURF  C--   Forcing term outside the Adams-Bashforth:
205        IF (.NOT. vectorInvariantMomentum  C     (not recommended with CD-scheme ON)
206       &    .AND. nonlinFreeSurf.GT.1) THEN        IF (momForcing .AND. .NOT.momForcing_In_AB) THEN
207         IF (select_rStar.GT.0) THEN         IF (useCDscheme) THEN
208          DO j=jMin,jMax          DO j=jMin,jMax
209           DO i=iMin,iMax           DO i=iMin,iMax
210             gUtmp(i,j) = gUtmp(i,j)/rStarExpW(i,j,bi,bj)            gUtmp(i,j) = gUtmp(i,j) - gU(i,j,k,bi,bj)
211              gVtmp(i,j) = gVtmp(i,j) - gV(i,j,k,bi,bj)
212           ENDDO           ENDDO
213          ENDDO          ENDDO
214         ELSE         ENDIF
215    #ifdef ALLOW_DIAGNOSTICS
216           IF ( useDiagnostics .AND. externForcDiagIsOn ) THEN
217            DO j=1,sNy+1
218             DO i=1,sNx+1
219              gUext(i,j) = gU(i,j,k,bi,bj)
220              gVext(i,j) = gV(i,j,k,bi,bj)
221             ENDDO
222            ENDDO
223           ENDIF
224    #endif /* ALLOW_DIAGNOSTICS */
225    
226            CALL EXTERNAL_FORCING_U(
227         I     iMin,iMax,jMin,jMax,bi,bj,k,
228         I     myTime,myThid)
229            CALL EXTERNAL_FORCING_V(
230         I     iMin,iMax,jMin,jMax,bi,bj,k,
231         I     myTime,myThid)
232    
233    #ifdef ALLOW_DIAGNOSTICS
234           IF ( useDiagnostics .AND. externForcDiagIsOn ) THEN
235            DO j=1,sNy+1
236             DO i=1,sNx+1
237              gUext(i,j) = gU(i,j,k,bi,bj)-gUext(i,j)
238              gVext(i,j) = gV(i,j,k,bi,bj)-gVext(i,j)
239             ENDDO
240            ENDDO
241           ENDIF
242    #endif /* ALLOW_DIAGNOSTICS */
243    
244    C-     for CD-scheme, compute gU,Vtmp = gU,V^n + forcing
245           IF (useCDscheme) THEN
246          DO j=jMin,jMax          DO j=jMin,jMax
247           DO i=iMin,iMax           DO i=iMin,iMax
248            IF ( k.EQ.ksurfW(i,j,bi,bj) ) THEN            gUtmp(i,j) = gUtmp(i,j) + gU(i,j,k,bi,bj)
249             gUtmp(i,j) = gUtmp(i,j)            gVtmp(i,j) = gVtmp(i,j) + gV(i,j,k,bi,bj)
      &         *hFacW(i,j,k,bi,bj)/hFac_surfW(i,j,bi,bj)  
           ENDIF  
250           ENDDO           ENDDO
251          ENDDO          ENDDO
252         ENDIF         ENDIF
253        ENDIF        ENDIF
 #endif  
254    
255  C     Step forward zonal velocity (store in Gu)  #ifdef ALLOW_CD_CODE
256        DO j=jMin,jMax        IF (useCDscheme) THEN
257          DO i=iMin,iMax  C-     Step forward D-grid velocity using C-grid gU,Vtmp = gU,V^n + forcing
258            gUNm1(i,j,k,bi,bj) = uVel(i,j,k,bi,bj)  C      and return coriolis terms on C-grid (guCor,gvCor)
259       &     +deltaTmom*(          CALL CD_CODE_SCHEME(
260       &         gUtmp(i,j)       I                  bi,bj,k, dPhiHydX,dPhiHydY, gUtmp,gVtmp,
261       &       - psFac*phiSurfX(i,j)       O                  guCor,gvCor,
262       &       - phxFac*dPhiHydX(i,j)       I                  myTime, myIter, myThid)
263       &        )*_maskW(i,j,k,bi,bj)          DO j=jMin,jMax
264             DO i=iMin,iMax
265              gUtmp(i,j) = gU(i,j,k,bi,bj)
266         &               + guCor(i,j)
267              gVtmp(i,j) = gV(i,j,k,bi,bj)
268         &               + gvCor(i,j)
269             ENDDO
270          ENDDO          ENDDO
271        ENDDO        ELSE
272    #endif /* ALLOW_CD_CODE */
273  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|          DO j=jMin,jMax
274  C-    Compute effective gV term (including Adams-Bashforth weights) :           DO i=iMin,iMax
275        DO j=jMin,jMax            gUtmp(i,j) = gU(i,j,k,bi,bj)
276         DO i=iMin,iMax            gVtmp(i,j) = gV(i,j,k,bi,bj)
277          gVtmp(i,j) = ab15*gV(i,j,k,bi,bj)           ENDDO
278       &             + ab05*gVNm1(i,j,k,bi,bj)            ENDDO
279  #ifdef INCLUDE_CD_CODE  #ifdef ALLOW_CD_CODE
280       &             + gvCD(i,j,k,bi,bj)        ENDIF
281  #endif  #endif
282         ENDDO  
       ENDDO  
         
283  #ifdef NONLIN_FRSURF  #ifdef NONLIN_FRSURF
284        IF (.NOT. vectorInvariantMomentum        IF (.NOT. vectorInvariantMomentum
285       &    .AND. nonlinFreeSurf.GT.1) THEN       &    .AND. nonlinFreeSurf.GT.1) THEN
286         IF (select_rStar.GT.0) THEN         IF (select_rStar.GT.0) THEN
287          DO j=jMin,jMax          DO j=jMin,jMax
288           DO i=iMin,iMax           DO i=iMin,iMax
289               gUtmp(i,j) = gUtmp(i,j)/rStarExpW(i,j,bi,bj)
290             gVtmp(i,j) = gVtmp(i,j)/rStarExpS(i,j,bi,bj)             gVtmp(i,j) = gVtmp(i,j)/rStarExpS(i,j,bi,bj)
291           ENDDO           ENDDO
292          ENDDO          ENDDO
293         ELSE         ELSE
294          DO j=jMin,jMax          DO j=jMin,jMax
295           DO i=iMin,iMax           DO i=iMin,iMax
296              IF ( k.EQ.ksurfW(i,j,bi,bj) ) THEN
297               gUtmp(i,j) = gUtmp(i,j)
298         &         *hFacW(i,j,k,bi,bj)/hFac_surfW(i,j,bi,bj)
299              ENDIF
300            IF ( k.EQ.ksurfS(i,j,bi,bj) ) THEN            IF ( k.EQ.ksurfS(i,j,bi,bj) ) THEN
301             gVtmp(i,j) = gVtmp(i,j)             gVtmp(i,j) = gVtmp(i,j)
302       &         *hFacS(i,j,k,bi,bj)/hFac_surfS(i,j,bi,bj)       &         *hFacS(i,j,k,bi,bj)/hFac_surfS(i,j,bi,bj)
# Line 147  C-    Compute effective gV term (includi Line 305  C-    Compute effective gV term (includi
305          ENDDO          ENDDO
306         ENDIF         ENDIF
307        ENDIF        ENDIF
308  #endif  #endif /* NONLIN_FRSURF */
309    
310    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
311    
312    C--   Dissipation term outside the Adams-Bashforth:
313          IF ( momViscosity .AND. .NOT.momDissip_In_AB ) THEN
314            DO j=jMin,jMax
315             DO i=iMin,iMax
316              gUtmp(i,j) = gUtmp(i,j) + guDissip(i,j)
317              gVtmp(i,j) = gVtmp(i,j) + gvDissip(i,j)
318             ENDDO
319            ENDDO
320          ENDIF
321    
322    C     Step forward zonal velocity (store in Gu)
323          DO j=jMin,jMax
324            DO i=iMin,iMax
325              gU(i,j,k,bi,bj) = uVel(i,j,k,bi,bj)
326         &     +deltaTmom*(
327         &         gUtmp(i,j)
328         &       - psFac*phiSurfX(i,j)
329         &       - phxFac*dPhiHydX(i,j)
330         &        )*_maskW(i,j,k,bi,bj)
331            ENDDO
332          ENDDO
333    
334  C     Step forward meridional velocity (store in Gv)  C     Step forward meridional velocity (store in Gv)
335        DO j=jMin,jMax        DO j=jMin,jMax
336          DO i=iMin,iMax          DO i=iMin,iMax
337            gVNm1(i,j,k,bi,bj) = vVel(i,j,k,bi,bj)            gV(i,j,k,bi,bj) = vVel(i,j,k,bi,bj)
338       &     +deltaTmom*(       &     +deltaTmom*(
339       &         gVtmp(i,j)       &         gVtmp(i,j)
340       &       - psFac*phiSurfY(i,j)       &       - psFac*phiSurfY(i,j)
# Line 161  C     Step forward meridional velocity ( Line 343  C     Step forward meridional velocity (
343          ENDDO          ENDDO
344        ENDDO        ENDDO
345    
346    #ifdef ALLOW_DIAGNOSTICS
347          IF ( useDiagnostics .AND. externForcDiagIsOn ) THEN
348              CALL DIAGNOSTICS_FILL(gUext,'Um_Ext  ',k,1,2,bi,bj,myThid)
349              CALL DIAGNOSTICS_FILL(gVext,'Vm_Ext  ',k,1,2,bi,bj,myThid)
350          ENDIF
351    #ifdef ALLOW_CD_CODE
352          IF ( useCDscheme .AND. useDiagnostics ) THEN
353              CALL DIAGNOSTICS_FILL(guCor,'Um_Cori ',k,1,2,bi,bj,myThid)
354              CALL DIAGNOSTICS_FILL(gvCor,'Vm_Cori ',k,1,2,bi,bj,myThid)
355          ENDIF
356    #endif /* ALLOW_CD_CODE */
357    #endif /* ALLOW_DIAGNOSTICS */
358    
359        RETURN        RETURN
360        END        END
Legend:
Removed from v.1.31  
changed lines
  Added in v.1.41
  ViewVC Help
Powered by ViewVC 1.1.22