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