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