/[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.12 by adcroft, Wed Jun 10 16:05:39 1998 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$
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"  #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 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
57          _RL     gUtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
58          _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
64    
65  C     Adams-Bashforth timestepping weights  C     Adams-Bashforth timestepping weights
66        ab15=1.5+abeps        IF (myIter .EQ. 0) THEN
67        ab05=-0.5-abeps         ab15=1.0
68           ab05=0.0
69          ELSE
70           ab15=1.5+abeps
71           ab05=-0.5-abeps
72          ENDIF
73    
74  C     Step forward zonal velocity (store in Gu)  C-- explicit part of the surface potential gradient is added in this S/R
75         DO j=jMin,jMax        psFac = pfFacMom*(1. _d 0 - implicSurfPress)
76          DO i=iMin,iMax  
77           gUNm1(i,j,k,bi,bj)=uVel(i,j,k,bi,bj)  C--  factors for gradient (X & Y directions) of Hydrostatic Potential
78       &    +deltaTmom*(ab15*gU(i,j,k,bi,bj)+ab05*gUNm1(i,j,k,bi,bj)        phxFac = pfFacMom
79  #ifdef ALLOW_CD        phyFac = pfFacMom
80       &    +guCD(i,j,k,bi,bj)  
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-|--+----|
87    
88    C-    Initialize local arrays (not really necessary but safer)
89          DO j=1-Oly,sNy+Oly
90           DO i=1-Olx,sNx+Olx
91            gUtmp(i,j) = 0. _d 0
92            gVtmp(i,j) = 0. _d 0
93    #ifdef ALLOW_CD_CODE
94            guCor(i,j) = 0. _d 0
95            gvCor(i,j) = 0. _d 0
96  #endif  #endif
      &        )*_maskW(i,j,k,bi,bj)  
         ENDDO  
97         ENDDO         ENDDO
98          ENDDO
99    
100  C     Step forward meridional velocity (store in Gv)        IF ( .NOT.staggerTimeStep ) THEN
101         DO j=jMin,jMax  C--   Synchronous time step: add grad Phi_Hyp to gU,gV before doing Adams-Bashforth
102          DO i=iMin,iMax          DO j=jMin,jMax
103           gVNm1(i,j,k,bi,bj)=vVel(i,j,k,bi,bj)           DO i=iMin,iMax
104       &    +deltaTmom*(ab15*gV(i,j,k,bi,bj)+ab05*gVNm1(i,j,k,bi,bj)            gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) - phxFac*dPhiHydX(i,j)
105  #ifdef ALLOW_CD            gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) - phyFac*dPhiHydY(i,j)
106       &    +gvCD(i,j,k,bi,bj)           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    C--   Forcing term outside the Adams-Bashforth:
167    C     (not recommanded with CD-scheme ON)
168          IF (momForcing .AND. .NOT.momForcing_In_AB) THEN
169           IF (useCDscheme) THEN
170            DO j=jMin,jMax
171             DO i=iMin,iMax
172              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
175            ENDDO
176           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
187             DO i=iMin,iMax
188              gUtmp(i,j) = gUtmp(i,j) + gU(i,j,k,bi,bj)
189              gVtmp(i,j) = gVtmp(i,j) + gV(i,j,k,bi,bj)
190             ENDDO
191            ENDDO
192           ENDIF
193          ENDIF
194    
195    #ifdef ALLOW_CD_CODE
196          IF (useCDscheme) THEN
197    C-     Step forward D-grid velocity using C-grid gU,Vtmp = gU,V^n + forcing
198    C      and return coriolis terms on C-grid (guCor,gvCor)
199            CALL CD_CODE_SCHEME(
200         I                  bi,bj,k, dPhiHydX,dPhiHydY, gUtmp,gVtmp,
201         O                  guCor,gvCor,
202         I                  myTime, myIter, myThid)
203            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
211          ELSE
212    #endif /* ALLOW_CD_CODE */
213            DO j=jMin,jMax
214             DO i=iMin,iMax
215              gUtmp(i,j) = gU(i,j,k,bi,bj)
216              gVtmp(i,j) = gV(i,j,k,bi,bj)
217             ENDDO
218            ENDDO
219    #ifdef ALLOW_CD_CODE
220          ENDIF
221  #endif  #endif
222       &        )*_maskS(i,j,k,bi,bj)  
223    #ifdef NONLIN_FRSURF
224          IF (.NOT. vectorInvariantMomentum
225         &    .AND. nonlinFreeSurf.GT.1) THEN
226           IF (select_rStar.GT.0) THEN
227            DO j=jMin,jMax
228             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)
231             ENDDO
232          ENDDO          ENDDO
233         ENDDO         ELSE
234            DO j=jMin,jMax
235             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
241               gVtmp(i,j) = gVtmp(i,j)
242         &         *hFacS(i,j,k,bi,bj)/hFac_surfS(i,j,bi,bj)
243              ENDIF
244             ENDDO
245            ENDDO
246           ENDIF
247          ENDIF
248    #endif /* NONLIN_FRSURF */
249    
250  C     Step forward temperature  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
251        IF (tempStepping) THEN  
252         DO j=jMin,jMax  C--   Dissipation term outside the Adams-Bashforth:
253          DO i=iMin,iMax        IF ( momViscosity .AND. .NOT.momDissip_In_AB ) THEN
254           theta(i,j,k,bi,bj)=theta(i,j,k,bi,bj)          DO j=jMin,jMax
255       &    +deltaTtracer*(ab15*gT(i,j,k,bi,bj)+ab05*gTNm1(i,j,k,bi,bj))           DO i=iMin,iMax
256           gTNm1(i,j,k,bi,bj)=gT(i,j,k,bi,bj)            gUtmp(i,j) = gUtmp(i,j) + guDissip(i,j)
257              gVtmp(i,j) = gVtmp(i,j) + gvDissip(i,j)
258             ENDDO
259          ENDDO          ENDDO
        ENDDO  
260        ENDIF        ENDIF
261    
262  C     Step forward salt  C     Step forward zonal velocity (store in Gu)
263        IF (saltStepping) THEN        DO j=jMin,jMax
        DO j=jMin,jMax  
264          DO i=iMin,iMax          DO i=iMin,iMax
265           salt(i,j,k,bi,bj)=salt(i,j,k,bi,bj)            gU(i,j,k,bi,bj) = uVel(i,j,k,bi,bj)
266       &    +deltaTtracer*(ab15*gS(i,j,k,bi,bj)+ab05*gSNm1(i,j,k,bi,bj))       &     +deltaTmom*(
267           gSNm1(i,j,k,bi,bj)=gS(i,j,k,bi,bj)       &         gUtmp(i,j)
268         &       - psFac*phiSurfX(i,j)
269         &       - phxFac*dPhiHydX(i,j)
270         &        )*_maskW(i,j,k,bi,bj)
271          ENDDO          ENDDO
272         ENDDO        ENDDO
273    
274    C     Step forward meridional velocity (store in Gv)
275          DO j=jMin,jMax
276            DO i=iMin,iMax
277              gV(i,j,k,bi,bj) = vVel(i,j,k,bi,bj)
278         &     +deltaTmom*(
279         &         gVtmp(i,j)
280         &       - psFac*phiSurfY(i,j)
281         &       - phyFac*dPhiHydY(i,j)
282         &        )*_maskS(i,j,k,bi,bj)
283            ENDDO
284          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        ENDIF
292    #endif /* ALLOW_CD_CODE */
293    #endif /* ALLOW_DIAGNOSTICS */
294    
295        RETURN        RETURN
296        END        END
Legend:
Removed from v.1.12  
changed lines
  Added in v.1.39
  ViewVC Help
Powered by ViewVC 1.1.22