/[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.32 by jmc, Thu Apr 17 13:41:34 2003 UTC
# Line 6  C $Name$ Line 6  C $Name$
6  CBOP  CBOP
7  C     !ROUTINE: TIMESTEP  C     !ROUTINE: TIMESTEP
8  C     !INTERFACE:  C     !INTERFACE:
9        SUBROUTINE TIMESTEP( bi, bj, iMin, iMax, jMin, jMax, K,        SUBROUTINE TIMESTEP( bi, bj, iMin, iMax, jMin, jMax, k,
10       I                     dPhiHydX,dPhiHydY, phiSurfX, phiSurfY,       I                     dPhiHydX,dPhiHydY, phiSurfX, phiSurfY,
11       I                     myIter, myThid )       I                     myTime, myIter, myThid )
12  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
13  C     *==========================================================*  C     *==========================================================*
14  C     | S/R TIMESTEP                                                C     | S/R TIMESTEP                                              
# Line 32  C     dPhiHydX,Y :: Gradient (X & Y dire Line 32  C     dPhiHydX,Y :: Gradient (X & Y dire
32  C     phiSurfX :: gradient of Surface potential (Pressure/rho, ocean)  C     phiSurfX :: gradient of Surface potential (Pressure/rho, ocean)
33  C     phiSurfY ::          or geopotential (atmos) in X and Y direction  C     phiSurfY ::          or geopotential (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     dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL     dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
37        _RL     dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _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:  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)        _RL     gUtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
50        _RL     gVtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _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  CEOP
56    
57  C     Adams-Bashforth timestepping weights  C     Adams-Bashforth timestepping weights
# Line 69  C-- stagger time step: grad Phi_Hyp is n Line 75  C-- stagger time step: grad Phi_Hyp is n
75  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
76        psFac = pfFacMom*(1. _d 0 - implicSurfPress)        psFac = pfFacMom*(1. _d 0 - implicSurfPress)
77    
78    C-- including or excluding momentum forcing from Adams-Bashforth:
79          momForcing_In_AB = forcing_In_AB
80          momForcing_In_AB = .TRUE.
81    
82  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
83  C-    Compute effective gU term (including Adams-Bashforth weights) :  
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
90            guCor(i,j) = 0. _d 0
91            gvCor(i,j) = 0. _d 0
92    #endif
93           ENDDO
94          ENDDO
95    
96    C--   Forcing term inside the Adams-Bashforth:
97          IF (momForcing .AND. momForcing_In_AB) THEN
98            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        DO j=jMin,jMax
109         DO i=iMin,iMax         DO i=iMin,iMax
110          gUtmp(i,j) = ab15*gU(i,j,k,bi,bj)          gUtmp(i,j) = ab15*gU(i,j,k,bi,bj)
111       &             + ab05*gUNm1(i,j,k,bi,bj)         &             + ab05*guNm1(i,j,k,bi,bj)  
112  #ifdef INCLUDE_CD_CODE          gVtmp(i,j) = ab15*gV(i,j,k,bi,bj)
113       &             + guCD(i,j,k,bi,bj)       &             + ab05*gvNm1(i,j,k,bi,bj)  
114  #endif          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         ENDDO
119        ENDDO        ENDDO
120                
121  #ifdef NONLIN_FRSURF  C--   Forcing term outside the Adams-Bashforth:
122        IF (.NOT. vectorInvariantMomentum  C     (not recommanded with CD-scheme ON)
123       &    .AND. nonlinFreeSurf.GT.1) THEN        IF (momForcing .AND. .NOT.momForcing_In_AB) THEN
124         IF (select_rStar.GT.0) THEN          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          DO j=jMin,jMax
133           DO i=iMin,iMax           DO i=iMin,iMax
134             gUtmp(i,j) = gUtmp(i,j)/rStarExpW(i,j,bi,bj)            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           ENDDO
139          ENDDO          ENDDO
140         ELSE         ELSE
141          DO j=jMin,jMax          DO j=jMin,jMax
142           DO i=iMin,iMax           DO i=iMin,iMax
143            IF ( k.EQ.ksurfW(i,j,bi,bj) ) THEN            gUtmp(i,j) = gU(i,j,k,bi,bj)
144             gUtmp(i,j) = gUtmp(i,j)            gVtmp(i,j) = gV(i,j,k,bi,bj)
      &         *hFacW(i,j,k,bi,bj)/hFac_surfW(i,j,bi,bj)  
           ENDIF  
145           ENDDO           ENDDO
146          ENDDO          ENDDO
147         ENDIF         ENDIF
148          ELSEIF ( useCDscheme) THEN
149            DO j=jMin,jMax
150             DO i=iMin,iMax
151              gUtmp(i,j) = guNm1(i,j,k,bi,bj)
152              gVtmp(i,j) = gvNm1(i,j,k,bi,bj)
153             ENDDO
154            ENDDO
155        ENDIF        ENDIF
 #endif  
156    
157  C     Step forward zonal velocity (store in Gu)  #ifdef INCLUDE_CD_CODE
158        DO j=jMin,jMax        IF (useCDscheme) THEN
159          DO i=iMin,iMax  C-     Step forward D-grid velocity using C-grid gU,Vtmp = gU,V^n + forcing
160            gUNm1(i,j,k,bi,bj) = uVel(i,j,k,bi,bj)  C      and return coriolis terms on C-grid (guCor,gvCor)
161       &     +deltaTmom*(          CALL MOM_CDSCHEME(
162       &         gUtmp(i,j)       I                    bi,bj,k, dPhiHydX,dPhiHydY, gUtmp,gVtmp,
163       &       - psFac*phiSurfX(i,j)       O                    guCor,gvCor,
164       &       - phxFac*dPhiHydX(i,j)       I                    myTime, myIter, myThid)
165       &        )*_maskW(i,j,k,bi,bj)          DO j=jMin,jMax
166             DO i=iMin,iMax
167              gUtmp(i,j) = gU(i,j,k,bi,bj)
168         &               + guCor(i,j)
169              gVtmp(i,j) = gV(i,j,k,bi,bj)
170         &               + gvCor(i,j)
171             ENDDO
172          ENDDO          ENDDO
173        ENDDO        ENDIF
174    #endif /* INCLUDE_CD_CODE */
175    
 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  
 C-    Compute effective gV term (including Adams-Bashforth weights) :  
       DO j=jMin,jMax  
        DO i=iMin,iMax  
         gVtmp(i,j) = ab15*gV(i,j,k,bi,bj)  
      &             + ab05*gVNm1(i,j,k,bi,bj)    
 #ifdef INCLUDE_CD_CODE  
      &             + gvCD(i,j,k,bi,bj)  
 #endif  
        ENDDO  
       ENDDO  
         
176  #ifdef NONLIN_FRSURF  #ifdef NONLIN_FRSURF
177        IF (.NOT. vectorInvariantMomentum        IF (.NOT. vectorInvariantMomentum
178       &    .AND. nonlinFreeSurf.GT.1) THEN       &    .AND. nonlinFreeSurf.GT.1) THEN
179         IF (select_rStar.GT.0) THEN         IF (select_rStar.GT.0) THEN
180          DO j=jMin,jMax          DO j=jMin,jMax
181           DO i=iMin,iMax           DO i=iMin,iMax
182               gUtmp(i,j) = gUtmp(i,j)/rStarExpW(i,j,bi,bj)
183             gVtmp(i,j) = gVtmp(i,j)/rStarExpS(i,j,bi,bj)             gVtmp(i,j) = gVtmp(i,j)/rStarExpS(i,j,bi,bj)
184           ENDDO           ENDDO
185          ENDDO          ENDDO
186         ELSE         ELSE
187          DO j=jMin,jMax          DO j=jMin,jMax
188           DO i=iMin,iMax           DO i=iMin,iMax
189              IF ( k.EQ.ksurfW(i,j,bi,bj) ) THEN
190               gUtmp(i,j) = gUtmp(i,j)
191         &         *hFacW(i,j,k,bi,bj)/hFac_surfW(i,j,bi,bj)
192              ENDIF
193            IF ( k.EQ.ksurfS(i,j,bi,bj) ) THEN            IF ( k.EQ.ksurfS(i,j,bi,bj) ) THEN
194             gVtmp(i,j) = gVtmp(i,j)             gVtmp(i,j) = gVtmp(i,j)
195       &         *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 198  C-    Compute effective gV term (includi
198          ENDDO          ENDDO
199         ENDIF         ENDIF
200        ENDIF        ENDIF
201  #endif  #endif /* NONLIN_FRSURF */
202    
203    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
204    
205    C     Step forward zonal velocity (store in Gu)
206          DO j=jMin,jMax
207            DO i=iMin,iMax
208              gU(i,j,k,bi,bj) = uVel(i,j,k,bi,bj)
209         &     +deltaTmom*(
210         &         gUtmp(i,j)
211         &       - psFac*phiSurfX(i,j)
212         &       - phxFac*dPhiHydX(i,j)
213         &        )*_maskW(i,j,k,bi,bj)
214            ENDDO
215          ENDDO
216    
217  C     Step forward meridional velocity (store in Gv)  C     Step forward meridional velocity (store in Gv)
218        DO j=jMin,jMax        DO j=jMin,jMax
219          DO i=iMin,iMax          DO i=iMin,iMax
220            gVNm1(i,j,k,bi,bj) = vVel(i,j,k,bi,bj)            gV(i,j,k,bi,bj) = vVel(i,j,k,bi,bj)
221       &     +deltaTmom*(       &     +deltaTmom*(
222       &         gVtmp(i,j)       &         gVtmp(i,j)
223       &       - psFac*phiSurfY(i,j)       &       - psFac*phiSurfY(i,j)
Legend:
Removed from v.1.31  
changed lines
  Added in v.1.32
  ViewVC Help
Powered by ViewVC 1.1.22