/[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.3 by cnh, Fri Apr 24 03:07:12 1998 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$
3    
4  #include "CPP_EEOPTIONS.h"  #include "PACKAGES_CONFIG.h"
5    #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, myThid )       I                     dPhiHydX,dPhiHydY, phiSurfX, phiSurfY,
12        implicit none       I                     guDissip, gvDissip,
13  ! Common       I                     myTime, myIter, myThid )
14    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 ==
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 "EEPARAMS.h"  #include "SURFACE.h"
30  #include "CG2D.h"  
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 myThid        INTEGER k
41          _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  C     pg       - Pressure gradient terms. Note cg2d_x        LOGICAL momForcing_In_AB
53  C                holds term in units so that lateral        LOGICAL momDissip_In_AB
54  C                gradient is all that is needed.        LOGICAL momStartAB
55        INTEGER i,j,k        INTEGER i,j
56        REAL ab15,ab05        _RL ab15,ab05
57        _RL  pg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _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        ab15=1.5+abeps        momStartAB = nIter0.EQ.0
76        ab05=-0.5-abeps        IF (myIter .EQ. 0) THEN
77           ab15=1.0
78           ab05=0.0
79          ELSE
80           ab15=1.5+abeps
81           ab05=-0.5-abeps
82          ENDIF
83    
84  C     Zonal pressure term  C-- explicit part of the surface potential gradient is added in this S/R
85        DO j=jMin,jMax        psFac = pfFacMom*(1. _d 0 - implicSurfPress)
86         DO i=iMin,iMax  
87          pg(i,j)=rDxC(i,j,bi,bj)*  C--  factors for gradient (X & Y directions) of Hydrostatic Potential
88       &   (cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))        phxFac = pfFacMom
89       &   *gravity*rhonil        phyFac = pfFacMom
90    
91    C-- including or excluding momentum forcing from Adams-Bashforth:
92          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
115         ENDDO         ENDDO
116        ENDDO        ENDDO
117  C     Step forward zonal velocity  
118        DO k=1,Nz        IF ( .NOT.staggerTimeStep ) THEN
119         DO j=jMin,jMax  C--   Synchronous time step: add grad Phi_Hyp to gU,gV before doing Adams-Bashforth
120          DO i=iMin,iMax          DO j=jMin,jMax
121           uVel(i,j,k,bi,bj)=uVel(i,j,k,bi,bj)           DO i=iMin,iMax
122       &    +deltaTmom*(ab15*gU(i,j,k,bi,bj)+ab05*gUNm1(i,j,k,bi,bj)            gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj) - phxFac*dPhiHydX(i,j)
123       &                -pg(i,j)/rhonil            gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj) - phyFac*dPhiHydY(i,j)
124       &    )*maskW(i,j,k,bi,bj)           ENDDO
          gUNm1(i,j,k,bi,bj)=gU(i,j,k,bi,bj)  
125          ENDDO          ENDDO
126         ENDDO          phxFac = 0.
127        ENDDO          phyFac = 0.
128  C     Meridional pressure term  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
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
290    #ifdef ALLOW_CD_CODE
291          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        DO j=jMin,jMax
        DO i=iMin,iMax  
         pg(i,j)=rDyC(i,j,bi,bj)*  
      &   (cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))  
      &   *gravity*rhonil  
        ENDDO  
       ENDDO  
 C     Step forward meridional velocity  
       DO k=1,Nz  
        DO j=jMin,jMax  
335          DO i=iMin,iMax          DO i=iMin,iMax
336           vVel(i,j,k,bi,bj)=vVel(i,j,k,bi,bj)            gU(i,j,k,bi,bj) = uVel(i,j,k,bi,bj)
337       &    +deltaTmom*(ab15*gV(i,j,k,bi,bj)+ab05*gVNm1(i,j,k,bi,bj)       &     +deltaTmom*(
338       &                -pg(i,j)/rhonil       &         gUtmp(i,j)
339       &               )*maskS(i,j,k,bi,bj)       &       - psFac*phiSurfX(i,j)
340           gVNm1(i,j,k,bi,bj)=gV(i,j,k,bi,bj)       &       - phxFac*dPhiHydX(i,j)
341         &        )*_maskW(i,j,k,bi,bj)
342          ENDDO          ENDDO
        ENDDO  
343        ENDDO        ENDDO
344  C     Step forward temperature  
345        DO k=1,Nz  C     Step forward meridional velocity (store in Gv)
346         DO j=jMin,jMax        DO j=jMin,jMax
347          DO i=iMin,iMax          DO i=iMin,iMax
348           theta(i,j,k,bi,bj)=theta(i,j,k,bi,bj)            gV(i,j,k,bi,bj) = vVel(i,j,k,bi,bj)
349       &    +deltaTtracer*(ab15*gT(i,j,k,bi,bj)+ab05*gTNm1(i,j,k,bi,bj))       &     +deltaTmom*(
350           gTNm1(i,j,k,bi,bj)=gT(i,j,k,bi,bj)       &         gVtmp(i,j)
351         &       - psFac*phiSurfY(i,j)
352         &       - phyFac*dPhiHydY(i,j)
353         &        )*_maskS(i,j,k,bi,bj)
354          ENDDO          ENDDO
        ENDDO  
355        ENDDO        ENDDO
356    
357        _BARRIER  #ifdef ALLOW_DIAGNOSTICS
358  C     CALL PLOT_FIELD_XYZR8( uVel, 'TIEMSTEP.1 uVel',Nz,1,myThid)        IF ( useDiagnostics .AND. externForcDiagIsOn ) THEN
359  C     CALL PLOT_FIELD_XYZR8( vVel, 'TIEMSTEP.1 vVel',Nz,1,myThid)            CALL DIAGNOSTICS_FILL(gUext,'Um_Ext  ',k,1,2,bi,bj,myThid)
360  C     CALL PLOT_FIELD_XYZR8( theta, 'TIEMSTEP.1 theta',Nz,1,myThid)            CALL DIAGNOSTICS_FILL(gVext,'Vm_Ext  ',k,1,2,bi,bj,myThid)
361          ENDIF
362    #ifdef ALLOW_CD_CODE
363          IF ( useCDscheme .AND. useDiagnostics ) THEN
364              CALL DIAGNOSTICS_FILL(guCor,'Um_Cori ',k,1,2,bi,bj,myThid)
365              CALL DIAGNOSTICS_FILL(gvCor,'Vm_Cori ',k,1,2,bi,bj,myThid)
366          ENDIF
367    #endif /* ALLOW_CD_CODE */
368    #endif /* ALLOW_DIAGNOSTICS */
369    
370        RETURN        RETURN
371        END        END
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.42
  ViewVC Help
Powered by ViewVC 1.1.22