/[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.40 by baylor, Thu Sep 29 16:47:16 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        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)
58          _RL     gVtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
59    C    Allow diagnosis of external forcing
60          _RL     gUext(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
61          _RL     gVext(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
62    #ifdef ALLOW_CD_CODE
63          _RL     guCor(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
64          _RL     gvCor(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
65    #endif
66    CEOP
67    
68  C     Adams-Bashforth timestepping weights  C     Adams-Bashforth timestepping weights
69        IF (myIter .EQ. 0) THEN        IF (myIter .EQ. 0) THEN
# Line 45  C     Adams-Bashforth timestepping weigh Line 74  C     Adams-Bashforth timestepping weigh
74         ab05=-0.5-abeps         ab05=-0.5-abeps
75        ENDIF        ENDIF
76    
77  C     Step forward zonal velocity (store in Gu)  C-- explicit part of the surface potential gradient is added in this S/R
78        psFac = pfFacMom*(1. _d 0 - implicSurfPress)        psFac = pfFacMom*(1. _d 0 - implicSurfPress)
79        DO j=jMin,jMax  
80          DO i=iMin,iMax  C--  factors for gradient (X & Y directions) of Hydrostatic Potential
81            gUNm1(i,j,k,bi,bj)=uVel(i,j,k,bi,bj) + deltaTmom * (        phxFac = pfFacMom
82       &       ab15*gU(i,j,k,bi,bj)+ab05*gUNm1(i,j,k,bi,bj)        phyFac = pfFacMom
83       &     - psFac*phiSurfX(i,j)  
84  #ifdef INCLUDE_CD_CODE  C-- including or excluding momentum forcing from Adams-Bashforth:
85       &     + guCD(i,j,k,bi,bj)        momForcing_In_AB = forcing_In_AB
86          momForcing_In_AB = .TRUE.
87          momDissip_In_AB  = .TRUE.
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            gUext(i,j) = 0. _d 0
97            gVext(i,j) = 0. _d 0
98    #ifdef ALLOW_CD_CODE
99            guCor(i,j) = 0. _d 0
100            gvCor(i,j) = 0. _d 0
101  #endif  #endif
102       &        )*_maskW(i,j,k,bi,bj)         ENDDO
         ENDDO  
103        ENDDO        ENDDO
104    
105        IF (staggerTimeStep) THEN        IF ( .NOT.staggerTimeStep ) THEN
106  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
107          phxFac = pfFacMom*deltaTmom          DO j=jMin,jMax
108             DO i=iMin,iMax
109              gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) - phxFac*dPhiHydX(i,j)
110              gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) - phyFac*dPhiHydY(i,j)
111             ENDDO
112            ENDDO
113            phxFac = 0.
114            phyFac = 0.
115    c     ELSE
116    C--   Stagger time step: grad Phi_Hyp will be added later
117          ENDIF
118    
119    C--   Dissipation term inside the Adams-Bashforth:
120          IF ( momViscosity .AND. momDissip_In_AB) THEN
121          DO j=jMin,jMax          DO j=jMin,jMax
122             DO i=iMin,iMax
123              gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) + guDissip(i,j)
124              gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) + gvDissip(i,j)
125             ENDDO
126            ENDDO
127          ENDIF
128    
129    C--   Forcing term inside the Adams-Bashforth:
130          IF (momForcing .AND. momForcing_In_AB) THEN
131    #ifdef ALLOW_DIAGNOSTICS
132            IF ( useDiagnostics ) THEN
133             DO j=jMin,jMax
134              DO i=iMin,iMax
135               gUext(i,j) = gU(i,j,k,bi,bj)
136               gVext(i,j) = gV(i,j,k,bi,bj)
137              ENDDO
138             ENDDO
139            ENDIF
140    #endif /* ALLOW_DIAGNOSTICS */
141    
142            CALL EXTERNAL_FORCING_U(
143         I     iMin,iMax,jMin,jMax,bi,bj,k,
144         I     myTime,myThid)
145            CALL EXTERNAL_FORCING_V(
146         I     iMin,iMax,jMin,jMax,bi,bj,k,
147         I     myTime,myThid)
148    
149    #ifdef ALLOW_DIAGNOSTICS
150            IF ( useDiagnostics ) THEN
151             DO j=jMin,jMax
152            DO i=iMin,iMax            DO i=iMin,iMax
153              gUNm1(i,j,k,bi,bj)=gUNm1(i,j,k,bi,bj)             gUext(i,j) = gU(i,j,k,bi,bj)-gUext(i,j)
154       &       - _recip_dxC(i,j,bi,bj)             gVext(i,j) = gV(i,j,k,bi,bj)-gVext(i,j)
      &         *(phiHyd(i,j,k)-phiHyd(i-1,j,k))*phxFac  
      &         *_maskW(i,j,k,bi,bj)  
155            ENDDO            ENDDO
156             ENDDO
157            ENDIF
158    #endif /* ALLOW_DIAGNOSTICS */
159          ENDIF
160    
161          IF (useCDscheme) THEN
162    C-     for CD-scheme, store gU,Vtmp = gU,V^n + forcing
163            DO j=jMin,jMax
164             DO i=iMin,iMax
165              gUtmp(i,j) = gU(i,j,k,bi,bj)
166              gVtmp(i,j) = gV(i,j,k,bi,bj)
167             ENDDO
168          ENDDO          ENDDO
169        ENDIF        ENDIF
170    
171    C-    Compute effective gU,gV_[n+1/2] terms (including Adams-Bashforth weights)
172    C     and save gU,gV_[n] into guNm1,gvNm1 for the next time step.
173    #ifdef ALLOW_ADAMSBASHFORTH_3
174          CALL ADAMS_BASHFORTH3(
175         I                        bi, bj, k,
176         U                        gU, guNm,
177         I                        myIter, myThid )
178          CALL ADAMS_BASHFORTH3(
179         I                        bi, bj, k,
180         U                        gV, gvNm,
181         I                        myIter, myThid )
182    #else /* ALLOW_ADAMSBASHFORTH_3 */
183          CALL ADAMS_BASHFORTH2(
184         I                        bi, bj, k,
185         U                        gU, guNm1,
186         I                        myIter, myThid )
187          CALL ADAMS_BASHFORTH2(
188         I                        bi, bj, k,
189         U                        gV, gvNm1,
190         I                        myIter, myThid )
191    #endif /* ALLOW_ADAMSBASHFORTH_3 */
192          
193    C--   Forcing term outside the Adams-Bashforth:
194    C     (not recommended with CD-scheme ON)
195          IF (momForcing .AND. .NOT.momForcing_In_AB) THEN
196           IF (useCDscheme) THEN
197            DO j=jMin,jMax
198             DO i=iMin,iMax
199              gUtmp(i,j) = gUtmp(i,j) - gU(i,j,k,bi,bj)
200              gVtmp(i,j) = gVtmp(i,j) - gV(i,j,k,bi,bj)
201             ENDDO
202            ENDDO
203           ENDIF
204    #ifdef ALLOW_DIAGNOSTICS
205           IF ( useDiagnostics ) THEN
206            DO j=jMin,jMax
207             DO i=iMin,iMax
208              gUext(i,j) = gU(i,j,k,bi,bj)
209              gVext(i,j) = gV(i,j,k,bi,bj)
210             ENDDO
211            ENDDO
212           ENDIF
213    #endif /* ALLOW_DIAGNOSTICS */
214    
215            CALL EXTERNAL_FORCING_U(
216         I     iMin,iMax,jMin,jMax,bi,bj,k,
217         I     myTime,myThid)
218            CALL EXTERNAL_FORCING_V(
219         I     iMin,iMax,jMin,jMax,bi,bj,k,
220         I     myTime,myThid)
221    
222    #ifdef ALLOW_DIAGNOSTICS
223           IF ( useDiagnostics ) THEN
224            DO j=jMin,jMax
225             DO i=iMin,iMax
226              gUext(i,j) = gU(i,j,k,bi,bj)-gUext(i,j)
227              gVext(i,j) = gV(i,j,k,bi,bj)-gVext(i,j)
228             ENDDO
229            ENDDO
230           ENDIF
231    #endif /* ALLOW_DIAGNOSTICS */
232    
233    #ifdef ALLOW_DIAGNOSTICS
234           IF ( useDiagnostics ) THEN
235             CALL DIAGNOSTICS_FILL(gUext,'Um_Ext  ',k,1,2,bi,bj,myThid)
236             CALL DIAGNOSTICS_FILL(gVext,'Vm_Ext  ',k,1,2,bi,bj,myThid)
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--   Dissipation term outside the Adams-Bashforth:
309          IF ( momViscosity .AND. .NOT.momDissip_In_AB ) THEN
310            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)
319          DO j=jMin,jMax
320            DO i=iMin,iMax
321              gU(i,j,k,bi,bj) = uVel(i,j,k,bi,bj)
322         &     +deltaTmom*(
323         &         gUtmp(i,j)
324         &       - psFac*phiSurfX(i,j)
325         &       - phxFac*dPhiHydX(i,j)
326         &        )*_maskW(i,j,k,bi,bj)
327            ENDDO
328          ENDDO
329    
330  C     Step forward meridional velocity (store in Gv)  C     Step forward meridional velocity (store in Gv)
       psFac = pfFacMom*(1. _d 0 - implicSurfPress)  
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) + deltaTmom * (            gV(i,j,k,bi,bj) = vVel(i,j,k,bi,bj)
334       &       ab15*gV(i,j,k,bi,bj)+ab05*gVNm1(i,j,k,bi,bj)       &     +deltaTmom*(
335       &     - psFac*phiSurfY(i,j)       &         gVtmp(i,j)
336  #ifdef INCLUDE_CD_CODE       &       - psFac*phiSurfY(i,j)
337       &     + gvCD(i,j,k,bi,bj)       &       - phyFac*dPhiHydY(i,j)
 #endif  
338       &        )*_maskS(i,j,k,bi,bj)       &        )*_maskS(i,j,k,bi,bj)
339          ENDDO          ENDDO
340        ENDDO        ENDDO
341    
342        IF (staggerTimeStep) THEN  #ifdef ALLOW_DIAGNOSTICS
343  C--   -grad Phi_Hyd has not been incorporated to gV and is added here:  #ifdef ALLOW_CD_CODE
344          phyFac = pfFacMom*deltaTmom        IF ( useCDscheme .AND. useDiagnostics ) THEN
345          DO j=jMin,jMax            CALL DIAGNOSTICS_FILL(guCor,'Um_Cori ',k,1,2,bi,bj,myThid)
346            DO i=iMin,iMax            CALL DIAGNOSTICS_FILL(gvCor,'Vm_Cori ',k,1,2,bi,bj,myThid)
             gVNm1(i,j,k,bi,bj)=gVNm1(i,j,k,bi,bj)  
      &       - _recip_dyC(i,j,bi,bj)  
      &         *(phiHyd(i,j,k)-phiHyd(i,j-1,k))*phyFac  
      &         *_maskS(i,j,k,bi,bj)  
           ENDDO  
         ENDDO  
347        ENDIF        ENDIF
348    #endif /* ALLOW_CD_CODE */
349    #endif /* ALLOW_DIAGNOSTICS */
350    
351        RETURN        RETURN
352        END        END
Legend:
Removed from v.1.25  
changed lines
  Added in v.1.40
  ViewVC Help
Powered by ViewVC 1.1.22