/[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.39 by jmc, Sun Sep 4 19:24:26 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_CD_CODE
60          _RL     guCor(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
61          _RL     gvCor(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
62    #endif
63  CEOP  CEOP
64    
65  C     Adams-Bashforth timestepping weights  C     Adams-Bashforth timestepping weights
# Line 57  C     Adams-Bashforth timestepping weigh Line 71  C     Adams-Bashforth timestepping weigh
71         ab05=-0.5-abeps         ab05=-0.5-abeps
72        ENDIF        ENDIF
73    
 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  
   
74  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
75        psFac = pfFacMom*(1. _d 0 - implicSurfPress)        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    C-- including or excluding momentum forcing from Adams-Bashforth:
82          momForcing_In_AB = forcing_In_AB
83          momForcing_In_AB = .TRUE.
84          momDissip_In_AB  = .TRUE.
85    
86  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
87  C-    Compute effective gU term (including Adams-Bashforth weights) :  
88        DO j=jMin,jMax  C-    Initialize local arrays (not really necessary but safer)
89         DO i=iMin,iMax        DO j=1-Oly,sNy+Oly
90          gUtmp(i,j) = ab15*gU(i,j,k,bi,bj)         DO i=1-Olx,sNx+Olx
91       &             + ab05*gUNm1(i,j,k,bi,bj)            gUtmp(i,j) = 0. _d 0
92  #ifdef INCLUDE_CD_CODE          gVtmp(i,j) = 0. _d 0
93       &             + guCD(i,j,k,bi,bj)  #ifdef ALLOW_CD_CODE
94            guCor(i,j) = 0. _d 0
95            gvCor(i,j) = 0. _d 0
96  #endif  #endif
97         ENDDO         ENDDO
98        ENDDO        ENDDO
99    
100          IF ( .NOT.staggerTimeStep ) THEN
101    C--   Synchronous time step: add grad Phi_Hyp to gU,gV before doing Adams-Bashforth
102            DO j=jMin,jMax
103             DO i=iMin,iMax
104              gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) - phxFac*dPhiHydX(i,j)
105              gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) - phyFac*dPhiHydY(i,j)
106             ENDDO
107            ENDDO
108            phxFac = 0.
109            phyFac = 0.
110    c     ELSE
111    C--   Stagger time step: grad Phi_Hyp will be added later
112          ENDIF
113    
114    C--   Dissipation term inside the Adams-Bashforth:
115          IF ( momViscosity .AND. momDissip_In_AB) THEN
116            DO j=jMin,jMax
117             DO i=iMin,iMax
118              gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) + guDissip(i,j)
119              gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) + gvDissip(i,j)
120             ENDDO
121            ENDDO
122          ENDIF
123    
124    C--   Forcing term inside the Adams-Bashforth:
125          IF (momForcing .AND. momForcing_In_AB) THEN
126            CALL EXTERNAL_FORCING_U(
127         I     iMin,iMax,jMin,jMax,bi,bj,k,
128         I     myTime,myThid)
129            CALL EXTERNAL_FORCING_V(
130         I     iMin,iMax,jMin,jMax,bi,bj,k,
131         I     myTime,myThid)
132          ENDIF
133    
134          IF (useCDscheme) THEN
135    C-     for CD-scheme, store gU,Vtmp = gU,V^n + forcing
136            DO j=jMin,jMax
137             DO i=iMin,iMax
138              gUtmp(i,j) = gU(i,j,k,bi,bj)
139              gVtmp(i,j) = gV(i,j,k,bi,bj)
140             ENDDO
141            ENDDO
142          ENDIF
143    
144    C-    Compute effective gU,gV_[n+1/2] terms (including Adams-Bashforth weights)
145    C     and save gU,gV_[n] into guNm1,gvNm1 for the next time step.
146    #ifdef ALLOW_ADAMSBASHFORTH_3
147          CALL ADAMS_BASHFORTH3(
148         I                        bi, bj, k,
149         U                        gU, guNm,
150         I                        myIter, myThid )
151          CALL ADAMS_BASHFORTH3(
152         I                        bi, bj, k,
153         U                        gV, gvNm,
154         I                        myIter, myThid )
155    #else /* ALLOW_ADAMSBASHFORTH_3 */
156          CALL ADAMS_BASHFORTH2(
157         I                        bi, bj, k,
158         U                        gU, guNm1,
159         I                        myIter, myThid )
160          CALL ADAMS_BASHFORTH2(
161         I                        bi, bj, k,
162         U                        gV, gvNm1,
163         I                        myIter, myThid )
164    #endif /* ALLOW_ADAMSBASHFORTH_3 */
165                
166  #ifdef NONLIN_FRSURF  C--   Forcing term outside the Adams-Bashforth:
167        IF (.NOT. vectorInvariantMomentum  C     (not recommanded with CD-scheme ON)
168       &    .AND. nonlinFreeSurf.GT.1) THEN        IF (momForcing .AND. .NOT.momForcing_In_AB) THEN
169         IF (select_rStar.GT.0) THEN         IF (useCDscheme) THEN
170          DO j=jMin,jMax          DO j=jMin,jMax
171           DO i=iMin,iMax           DO i=iMin,iMax
172             gUtmp(i,j) = gUtmp(i,j)/rStarExpW(i,j,bi,bj)            gUtmp(i,j) = gUtmp(i,j) - gU(i,j,k,bi,bj)
173              gVtmp(i,j) = gVtmp(i,j) - gV(i,j,k,bi,bj)
174           ENDDO           ENDDO
175          ENDDO          ENDDO
176         ELSE         ENDIF
177            CALL EXTERNAL_FORCING_U(
178         I     iMin,iMax,jMin,jMax,bi,bj,k,
179         I     myTime,myThid)
180            CALL EXTERNAL_FORCING_V(
181         I     iMin,iMax,jMin,jMax,bi,bj,k,
182         I     myTime,myThid)
183    
184    C-     for CD-scheme, compute gU,Vtmp = gU,V^n + forcing
185           IF (useCDscheme) THEN
186          DO j=jMin,jMax          DO j=jMin,jMax
187           DO i=iMin,iMax           DO i=iMin,iMax
188            IF ( k.EQ.ksurfW(i,j,bi,bj) ) THEN            gUtmp(i,j) = gUtmp(i,j) + gU(i,j,k,bi,bj)
189             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  
190           ENDDO           ENDDO
191          ENDDO          ENDDO
192         ENDIF         ENDIF
193        ENDIF        ENDIF
 #endif  
194    
195  C     Step forward zonal velocity (store in Gu)  #ifdef ALLOW_CD_CODE
196        DO j=jMin,jMax        IF (useCDscheme) THEN
197          DO i=iMin,iMax  C-     Step forward D-grid velocity using C-grid gU,Vtmp = gU,V^n + forcing
198            gUNm1(i,j,k,bi,bj) = uVel(i,j,k,bi,bj)  C      and return coriolis terms on C-grid (guCor,gvCor)
199       &     +deltaTmom*(          CALL CD_CODE_SCHEME(
200       &         gUtmp(i,j)       I                  bi,bj,k, dPhiHydX,dPhiHydY, gUtmp,gVtmp,
201       &       - psFac*phiSurfX(i,j)       O                  guCor,gvCor,
202       &       - phxFac*dPhiHydX(i,j)       I                  myTime, myIter, myThid)
203       &        )*_maskW(i,j,k,bi,bj)          DO j=jMin,jMax
204             DO i=iMin,iMax
205              gUtmp(i,j) = gU(i,j,k,bi,bj)
206         &               + guCor(i,j)
207              gVtmp(i,j) = gV(i,j,k,bi,bj)
208         &               + gvCor(i,j)
209             ENDDO
210          ENDDO          ENDDO
211        ENDDO        ELSE
212    #endif /* ALLOW_CD_CODE */
213  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|          DO j=jMin,jMax
214  C-    Compute effective gV term (including Adams-Bashforth weights) :           DO i=iMin,iMax
215        DO j=jMin,jMax            gUtmp(i,j) = gU(i,j,k,bi,bj)
216         DO i=iMin,iMax            gVtmp(i,j) = gV(i,j,k,bi,bj)
217          gVtmp(i,j) = ab15*gV(i,j,k,bi,bj)           ENDDO
218       &             + ab05*gVNm1(i,j,k,bi,bj)            ENDDO
219  #ifdef INCLUDE_CD_CODE  #ifdef ALLOW_CD_CODE
220       &             + gvCD(i,j,k,bi,bj)        ENDIF
221  #endif  #endif
222         ENDDO  
       ENDDO  
         
223  #ifdef NONLIN_FRSURF  #ifdef NONLIN_FRSURF
224        IF (.NOT. vectorInvariantMomentum        IF (.NOT. vectorInvariantMomentum
225       &    .AND. nonlinFreeSurf.GT.1) THEN       &    .AND. nonlinFreeSurf.GT.1) THEN
226         IF (select_rStar.GT.0) THEN         IF (select_rStar.GT.0) THEN
227          DO j=jMin,jMax          DO j=jMin,jMax
228           DO i=iMin,iMax           DO i=iMin,iMax
229               gUtmp(i,j) = gUtmp(i,j)/rStarExpW(i,j,bi,bj)
230             gVtmp(i,j) = gVtmp(i,j)/rStarExpS(i,j,bi,bj)             gVtmp(i,j) = gVtmp(i,j)/rStarExpS(i,j,bi,bj)
231           ENDDO           ENDDO
232          ENDDO          ENDDO
233         ELSE         ELSE
234          DO j=jMin,jMax          DO j=jMin,jMax
235           DO i=iMin,iMax           DO i=iMin,iMax
236              IF ( k.EQ.ksurfW(i,j,bi,bj) ) THEN
237               gUtmp(i,j) = gUtmp(i,j)
238         &         *hFacW(i,j,k,bi,bj)/hFac_surfW(i,j,bi,bj)
239              ENDIF
240            IF ( k.EQ.ksurfS(i,j,bi,bj) ) THEN            IF ( k.EQ.ksurfS(i,j,bi,bj) ) THEN
241             gVtmp(i,j) = gVtmp(i,j)             gVtmp(i,j) = gVtmp(i,j)
242       &         *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 245  C-    Compute effective gV term (includi
245          ENDDO          ENDDO
246         ENDIF         ENDIF
247        ENDIF        ENDIF
248  #endif  #endif /* NONLIN_FRSURF */
249    
250    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
251    
252    C--   Dissipation term outside the Adams-Bashforth:
253          IF ( momViscosity .AND. .NOT.momDissip_In_AB ) THEN
254            DO j=jMin,jMax
255             DO i=iMin,iMax
256              gUtmp(i,j) = gUtmp(i,j) + guDissip(i,j)
257              gVtmp(i,j) = gVtmp(i,j) + gvDissip(i,j)
258             ENDDO
259            ENDDO
260          ENDIF
261    
262    C     Step forward zonal velocity (store in Gu)
263          DO j=jMin,jMax
264            DO i=iMin,iMax
265              gU(i,j,k,bi,bj) = uVel(i,j,k,bi,bj)
266         &     +deltaTmom*(
267         &         gUtmp(i,j)
268         &       - psFac*phiSurfX(i,j)
269         &       - phxFac*dPhiHydX(i,j)
270         &        )*_maskW(i,j,k,bi,bj)
271            ENDDO
272          ENDDO
273    
274  C     Step forward meridional velocity (store in Gv)  C     Step forward meridional velocity (store in Gv)
275        DO j=jMin,jMax        DO j=jMin,jMax
276          DO i=iMin,iMax          DO i=iMin,iMax
277            gVNm1(i,j,k,bi,bj) = vVel(i,j,k,bi,bj)            gV(i,j,k,bi,bj) = vVel(i,j,k,bi,bj)
278       &     +deltaTmom*(       &     +deltaTmom*(
279       &         gVtmp(i,j)       &         gVtmp(i,j)
280       &       - psFac*phiSurfY(i,j)       &       - psFac*phiSurfY(i,j)
# Line 161  C     Step forward meridional velocity ( Line 283  C     Step forward meridional velocity (
283          ENDDO          ENDDO
284        ENDDO        ENDDO
285    
286    #ifdef ALLOW_DIAGNOSTICS
287    #ifdef ALLOW_CD_CODE
288          IF ( useCDscheme .AND. useDiagnostics ) THEN
289              CALL DIAGNOSTICS_FILL(guCor,'Um_Cori ',k,1,2,bi,bj,myThid)
290              CALL DIAGNOSTICS_FILL(gvCor,'Vm_Cori ',k,1,2,bi,bj,myThid)
291          ENDIF
292    #endif /* ALLOW_CD_CODE */
293    #endif /* ALLOW_DIAGNOSTICS */
294    
295        RETURN        RETURN
296        END        END
Legend:
Removed from v.1.31  
changed lines
  Added in v.1.39
  ViewVC Help
Powered by ViewVC 1.1.22