/[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.23 by jmc, Thu Mar 8 20:27:33 2001 UTC revision 1.33 by edhill, Thu Oct 9 04:19:18 2003 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
5    
6  C     /==========================================================\  CBOP
7  C     | S/R TIMESTEP                                             |  C     !ROUTINE: TIMESTEP
8  C     | o Step model fields forward in time                      |  C     !INTERFACE:
9  C     \==========================================================/        SUBROUTINE TIMESTEP( bi, bj, iMin, iMax, jMin, jMax, k,
10        SUBROUTINE TIMESTEP( bi, bj, iMin, iMax, jMin, jMax, K,       I                     dPhiHydX,dPhiHydY, phiSurfX, phiSurfY,
11       I                     phiHyd, phiSurfX, phiSurfY,       I                     myTime, myIter, myThid )
12       I                     myIter, myThid )  C     !DESCRIPTION: \bv
13        IMPLICIT NONE  C     *==========================================================*
14    C     | S/R TIMESTEP                                              
15    C     | o Step model fields forward in time                      
16    C     *==========================================================*
17    C     \ev
18    
19    C     !USES:
20          IMPLICIT NONE
21  C     == Global variables ==  C     == Global variables ==
22  #include "SIZE.h"  #include "SIZE.h"
23  #include "DYNVARS.h"  #include "DYNVARS.h"
24  #include "EEPARAMS.h"  #include "EEPARAMS.h"
25  #include "PARAMS.h"  #include "PARAMS.h"
26  #include "GRID.h"  #include "GRID.h"
27    #include "SURFACE.h"
28    
29    C     !INPUT/OUTPUT PARAMETERS:
30  C     == Routine Arguments ==  C     == Routine Arguments ==
31  C     phiHyd    - Hydrostatic Potential (ocean: pressure/rho)  C     dPhiHydX,Y :: Gradient (X & Y directions) of Hydrostatic Potential
32  C                                       (atmos: geopotentiel)  C     phiSurfX :: gradient of Surface potential (Pressure/rho, ocean)
33  C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)  C     phiSurfY ::          or geopotential (atmos) in X and Y direction
 C     phiSurfY             or geopotentiel (atmos) in X and Y direction  
34        INTEGER bi,bj,iMin,iMax,jMin,jMax        INTEGER bi,bj,iMin,iMax,jMin,jMax
35        INTEGER K        INTEGER k
36        _RL     phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
37          _RL     dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
38        _RL     phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
39        _RL     phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
40          _RL     myTime
41        INTEGER myIter, myThid        INTEGER myIter, myThid
42    
43    C     !LOCAL VARIABLES:
44  C     == Local variables ==  C     == Local variables ==
45          LOGICAL momForcing_In_AB
46        INTEGER i,j        INTEGER i,j
47        _RL ab15,ab05        _RL ab15,ab05
48        _RL phxFac,phyFac, psFac        _RL phxFac,phyFac, psFac
49          _RL     gUtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
50          _RL     gVtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
51    #ifdef INCLUDE_CD_CODE
52          _RL     guCor(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
53          _RL     gvCor(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
54    #endif
55    CEOP
56    
57  C     Adams-Bashforth timestepping weights  C     Adams-Bashforth timestepping weights
58  Caja  IF (myIter .EQ. 0) THEN        IF (myIter .EQ. 0) THEN
59  Caja   ab15=1.0         ab15=1.0
60  Caja   ab05=0.0         ab05=0.0
61  Caja  ELSE        ELSE
62         ab15=1.5+abeps         ab15=1.5+abeps
63         ab05=-0.5-abeps         ab05=-0.5-abeps
64  Caja  ENDIF        ENDIF
65    
66  C     Step forward zonal velocity (store in Gu)  C-- stagger time step: grad Phi_Hyp is not in gU,gV => add it in this S/R
67          IF (staggerTimeStep) THEN
68            phxFac = pfFacMom
69            phyFac = pfFacMom
70          ELSE
71            phxFac = 0.
72            phyFac = 0.
73          ENDIF
74    
75    C-- explicit part of the surface potential gradient is added in this S/R
76        psFac = pfFacMom*(1. _d 0 - implicSurfPress)        psFac = pfFacMom*(1. _d 0 - implicSurfPress)
77        DO j=jMin,jMax  
78          DO i=iMin,iMax  C-- including or excluding momentum forcing from Adams-Bashforth:
79            gUNm1(i,j,k,bi,bj)=uVel(i,j,k,bi,bj) + deltaTmom * (        momForcing_In_AB = forcing_In_AB
80       &       ab15*gU(i,j,k,bi,bj)+ab05*gUNm1(i,j,k,bi,bj)        momForcing_In_AB = .TRUE.
81       &     - psFac*phiSurfX(i,j)  
82    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
83    
84    C-    Initialize local arrays (not really necessary but safer)
85          DO j=1-Oly,sNy+Oly
86           DO i=1-Olx,sNx+Olx
87            gUtmp(i,j) = 0. _d 0
88            gVtmp(i,j) = 0. _d 0
89  #ifdef INCLUDE_CD_CODE  #ifdef INCLUDE_CD_CODE
90       &     + guCD(i,j,k,bi,bj)          guCor(i,j) = 0. _d 0
91            gvCor(i,j) = 0. _d 0
92  #endif  #endif
93       &        )*_maskW(i,j,k,bi,bj)         ENDDO
         ENDDO  
94        ENDDO        ENDDO
95    
96        IF (staggerTimeStep) THEN  C--   Forcing term inside the Adams-Bashforth:
97  C--   -grad Phi_Hyd has not been incorporated to gU and is added here:        IF (momForcing .AND. momForcing_In_AB) THEN
98          phxFac = pfFacMom*deltaTmom          CALL EXTERNAL_FORCING_U(
99         I     iMin,iMax,jMin,jMax,bi,bj,k,
100         I     myTime,myThid)
101            CALL EXTERNAL_FORCING_V(
102         I     iMin,iMax,jMin,jMax,bi,bj,k,
103         I     myTime,myThid)
104          ENDIF
105    
106    C-    Compute effective gU,gV_[n+1/2] terms (including Adams-Bashforth weights)
107    C     and save gU,gV_[n] into guNm1,gvNm1 for the next time step.
108          DO j=jMin,jMax
109           DO i=iMin,iMax
110            gUtmp(i,j) = ab15*gU(i,j,k,bi,bj)
111         &             + ab05*guNm1(i,j,k,bi,bj)  
112            gVtmp(i,j) = ab15*gV(i,j,k,bi,bj)
113         &             + ab05*gvNm1(i,j,k,bi,bj)  
114            guNm1(i,j,k,bi,bj)= gU(i,j,k,bi,bj)
115            gvNm1(i,j,k,bi,bj)= gV(i,j,k,bi,bj)
116            gU(i,j,k,bi,bj) = gUtmp(i,j)
117            gV(i,j,k,bi,bj) = gVtmp(i,j)
118           ENDDO
119          ENDDO
120          
121    C--   Forcing term outside the Adams-Bashforth:
122    C     (not recommanded with CD-scheme ON)
123          IF (momForcing .AND. .NOT.momForcing_In_AB) THEN
124            CALL EXTERNAL_FORCING_U(
125         I     iMin,iMax,jMin,jMax,bi,bj,k,
126         I     myTime,myThid)
127            CALL EXTERNAL_FORCING_V(
128         I     iMin,iMax,jMin,jMax,bi,bj,k,
129         I     myTime,myThid)
130           IF (useCDscheme) THEN
131    C-     for CD-scheme, compute gU,Vtmp = gU,V^n + forcing
132            DO j=jMin,jMax
133             DO i=iMin,iMax
134              gUtmp(i,j) = gU(i,j,k,bi,bj)-gUtmp(i,j)
135         &               + guNm1(i,j,k,bi,bj)
136              gVtmp(i,j) = gV(i,j,k,bi,bj)-gVtmp(i,j)
137         &               + gvNm1(i,j,k,bi,bj)
138             ENDDO
139            ENDDO
140           ELSE
141            DO j=jMin,jMax
142             DO i=iMin,iMax
143              gUtmp(i,j) = gU(i,j,k,bi,bj)
144              gVtmp(i,j) = gV(i,j,k,bi,bj)
145             ENDDO
146            ENDDO
147           ENDIF
148          ELSEIF ( useCDscheme) THEN
149          DO j=jMin,jMax          DO j=jMin,jMax
150            DO i=iMin,iMax           DO i=iMin,iMax
151              gUNm1(i,j,k,bi,bj)=gUNm1(i,j,k,bi,bj)            gUtmp(i,j) = guNm1(i,j,k,bi,bj)
152       &       - _recip_dxC(i,j,bi,bj)            gVtmp(i,j) = gvNm1(i,j,k,bi,bj)
153       &         *(phiHyd(i,j,k)-phiHyd(i-1,j,k))*phxFac           ENDDO
      &         *_maskW(i,j,k,bi,bj)  
           ENDDO  
154          ENDDO          ENDDO
155        ENDIF        ENDIF
156    
 C     Step forward meridional velocity (store in Gv)  
       psFac = pfFacMom*(1. _d 0 - implicSurfPress)  
       DO j=jMin,jMax  
         DO i=iMin,iMax  
           gVNm1(i,j,k,bi,bj)=vVel(i,j,k,bi,bj) + deltaTmom * (  
      &       ab15*gV(i,j,k,bi,bj)+ab05*gVNm1(i,j,k,bi,bj)  
      &     - psFac*phiSurfY(i,j)  
157  #ifdef INCLUDE_CD_CODE  #ifdef INCLUDE_CD_CODE
158       &     + gvCD(i,j,k,bi,bj)  #ifdef ALLOW_MOM_FLUXFORM
159  #endif        IF (useCDscheme) THEN
160       &        )*_maskS(i,j,k,bi,bj)  C-     Step forward D-grid velocity using C-grid gU,Vtmp = gU,V^n + forcing
161    C      and return coriolis terms on C-grid (guCor,gvCor)
162            CALL MOM_CDSCHEME(
163         I                    bi,bj,k, dPhiHydX,dPhiHydY, gUtmp,gVtmp,
164         O                    guCor,gvCor,
165         I                    myTime, myIter, myThid)
166            DO j=jMin,jMax
167             DO i=iMin,iMax
168              gUtmp(i,j) = gU(i,j,k,bi,bj)
169         &               + guCor(i,j)
170              gVtmp(i,j) = gV(i,j,k,bi,bj)
171         &               + gvCor(i,j)
172             ENDDO
173          ENDDO          ENDDO
174        ENDDO        ENDIF
175    #endif /* ALLOW_MOM_FLUXFORM */
176    #endif /* INCLUDE_CD_CODE */
177    
178        IF (staggerTimeStep) THEN  #ifdef NONLIN_FRSURF
179  C--   -grad Phi_Hyd has not been incorporated to gV and is added here:        IF (.NOT. vectorInvariantMomentum
180          phyFac = pfFacMom*deltaTmom       &    .AND. nonlinFreeSurf.GT.1) THEN
181           IF (select_rStar.GT.0) THEN
182            DO j=jMin,jMax
183             DO i=iMin,iMax
184               gUtmp(i,j) = gUtmp(i,j)/rStarExpW(i,j,bi,bj)
185               gVtmp(i,j) = gVtmp(i,j)/rStarExpS(i,j,bi,bj)
186             ENDDO
187            ENDDO
188           ELSE
189          DO j=jMin,jMax          DO j=jMin,jMax
190            DO i=iMin,iMax           DO i=iMin,iMax
191              gVNm1(i,j,k,bi,bj)=gVNm1(i,j,k,bi,bj)            IF ( k.EQ.ksurfW(i,j,bi,bj) ) THEN
192       &       - _recip_dyC(i,j,bi,bj)             gUtmp(i,j) = gUtmp(i,j)
193       &         *(phiHyd(i,j,k)-phiHyd(i,j-1,k))*phyFac                             &         *hFacW(i,j,k,bi,bj)/hFac_surfW(i,j,bi,bj)
194       &         *_maskS(i,j,k,bi,bj)            ENDIF
195            ENDDO            IF ( k.EQ.ksurfS(i,j,bi,bj) ) THEN
196               gVtmp(i,j) = gVtmp(i,j)
197         &         *hFacS(i,j,k,bi,bj)/hFac_surfS(i,j,bi,bj)
198              ENDIF
199             ENDDO
200          ENDDO          ENDDO
201           ENDIF
202        ENDIF        ENDIF
203    #endif /* NONLIN_FRSURF */
204    
205    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
206    
207    C     Step forward zonal velocity (store in Gu)
208          DO j=jMin,jMax
209            DO i=iMin,iMax
210              gU(i,j,k,bi,bj) = uVel(i,j,k,bi,bj)
211         &     +deltaTmom*(
212         &         gUtmp(i,j)
213         &       - psFac*phiSurfX(i,j)
214         &       - phxFac*dPhiHydX(i,j)
215         &        )*_maskW(i,j,k,bi,bj)
216            ENDDO
217          ENDDO
218    
219    C     Step forward meridional velocity (store in Gv)
220          DO j=jMin,jMax
221            DO i=iMin,iMax
222              gV(i,j,k,bi,bj) = vVel(i,j,k,bi,bj)
223         &     +deltaTmom*(
224         &         gVtmp(i,j)
225         &       - psFac*phiSurfY(i,j)
226         &       - phyFac*dPhiHydY(i,j)
227         &        )*_maskS(i,j,k,bi,bj)
228            ENDDO
229          ENDDO
230    
231        RETURN        RETURN
232        END        END
Legend:
Removed from v.1.23  
changed lines
  Added in v.1.33
  ViewVC Help
Powered by ViewVC 1.1.22