/[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.31 by jmc, Tue Feb 18 15:27:25 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:
 C     \==========================================================/  
9        SUBROUTINE TIMESTEP( bi, bj, iMin, iMax, jMin, jMax, K,        SUBROUTINE TIMESTEP( bi, bj, iMin, iMax, jMin, jMax, K,
10       I                     phiHyd, phiSurfX, phiSurfY,       I                     dPhiHydX,dPhiHydY, phiSurfX, phiSurfY,
11       I                     myIter, myThid )       I                     myIter, myThid )
12        IMPLICIT NONE  C     !DESCRIPTION: \bv
13    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        INTEGER myIter, myThid        INTEGER myIter, myThid
41    
42    C     !LOCAL VARIABLES:
43  C     == Local variables ==  C     == Local variables ==
44        INTEGER i,j        INTEGER i,j
45        _RL ab15,ab05        _RL ab15,ab05
46        _RL phxFac,phyFac, psFac        _RL phxFac,phyFac, psFac
47          _RL     gUtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
48          _RL     gVtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
49    CEOP
50    
51  C     Adams-Bashforth timestepping weights  C     Adams-Bashforth timestepping weights
52  Caja  IF (myIter .EQ. 0) THEN        IF (myIter .EQ. 0) THEN
53  Caja   ab15=1.0         ab15=1.0
54  Caja   ab05=0.0         ab05=0.0
55  Caja  ELSE        ELSE
56         ab15=1.5+abeps         ab15=1.5+abeps
57         ab05=-0.5-abeps         ab05=-0.5-abeps
58  Caja  ENDIF        ENDIF
59    
60  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
61          IF (staggerTimeStep) THEN
62            phxFac = pfFacMom
63            phyFac = pfFacMom
64          ELSE
65            phxFac = 0.
66            phyFac = 0.
67          ENDIF
68    
69    C-- explicit part of the surface potential gradient is added in this S/R
70        psFac = pfFacMom*(1. _d 0 - implicSurfPress)        psFac = pfFacMom*(1. _d 0 - implicSurfPress)
71    
72    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
73    C-    Compute effective gU term (including Adams-Bashforth weights) :
74        DO j=jMin,jMax        DO j=jMin,jMax
75          DO i=iMin,iMax         DO i=iMin,iMax
76            gUNm1(i,j,k,bi,bj)=uVel(i,j,k,bi,bj) + deltaTmom * (          gUtmp(i,j) = ab15*gU(i,j,k,bi,bj)
77       &       ab15*gU(i,j,k,bi,bj)+ab05*gUNm1(i,j,k,bi,bj)       &             + ab05*gUNm1(i,j,k,bi,bj)  
      &     - psFac*phiSurfX(i,j)  
78  #ifdef INCLUDE_CD_CODE  #ifdef INCLUDE_CD_CODE
79       &     + guCD(i,j,k,bi,bj)       &             + guCD(i,j,k,bi,bj)
80  #endif  #endif
81           ENDDO
82          ENDDO
83          
84    #ifdef NONLIN_FRSURF
85          IF (.NOT. vectorInvariantMomentum
86         &    .AND. nonlinFreeSurf.GT.1) THEN
87           IF (select_rStar.GT.0) THEN
88            DO j=jMin,jMax
89             DO i=iMin,iMax
90               gUtmp(i,j) = gUtmp(i,j)/rStarExpW(i,j,bi,bj)
91             ENDDO
92            ENDDO
93           ELSE
94            DO j=jMin,jMax
95             DO i=iMin,iMax
96              IF ( k.EQ.ksurfW(i,j,bi,bj) ) THEN
97               gUtmp(i,j) = gUtmp(i,j)
98         &         *hFacW(i,j,k,bi,bj)/hFac_surfW(i,j,bi,bj)
99              ENDIF
100             ENDDO
101            ENDDO
102           ENDIF
103          ENDIF
104    #endif
105    
106    C     Step forward zonal velocity (store in Gu)
107          DO j=jMin,jMax
108            DO i=iMin,iMax
109              gUNm1(i,j,k,bi,bj) = uVel(i,j,k,bi,bj)
110         &     +deltaTmom*(
111         &         gUtmp(i,j)
112         &       - psFac*phiSurfX(i,j)
113         &       - phxFac*dPhiHydX(i,j)
114       &        )*_maskW(i,j,k,bi,bj)       &        )*_maskW(i,j,k,bi,bj)
115          ENDDO          ENDDO
116        ENDDO        ENDDO
117    
118        IF (staggerTimeStep) THEN  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
119  C--   -grad Phi_Hyd has not been incorporated to gU and is added here:  C-    Compute effective gV term (including Adams-Bashforth weights) :
120          phxFac = pfFacMom*deltaTmom        DO j=jMin,jMax
121           DO i=iMin,iMax
122            gVtmp(i,j) = ab15*gV(i,j,k,bi,bj)
123         &             + ab05*gVNm1(i,j,k,bi,bj)  
124    #ifdef INCLUDE_CD_CODE
125         &             + gvCD(i,j,k,bi,bj)
126    #endif
127           ENDDO
128          ENDDO
129          
130    #ifdef NONLIN_FRSURF
131          IF (.NOT. vectorInvariantMomentum
132         &    .AND. nonlinFreeSurf.GT.1) THEN
133           IF (select_rStar.GT.0) THEN
134            DO j=jMin,jMax
135             DO i=iMin,iMax
136               gVtmp(i,j) = gVtmp(i,j)/rStarExpS(i,j,bi,bj)
137             ENDDO
138            ENDDO
139           ELSE
140          DO j=jMin,jMax          DO j=jMin,jMax
141            DO i=iMin,iMax           DO i=iMin,iMax
142              gUNm1(i,j,k,bi,bj)=gUNm1(i,j,k,bi,bj)            IF ( k.EQ.ksurfS(i,j,bi,bj) ) THEN
143       &       - _recip_dxC(i,j,bi,bj)             gVtmp(i,j) = gVtmp(i,j)
144       &         *(phiHyd(i,j,k)-phiHyd(i-1,j,k))*phxFac       &         *hFacS(i,j,k,bi,bj)/hFac_surfS(i,j,bi,bj)
145       &         *_maskW(i,j,k,bi,bj)            ENDIF
146            ENDDO           ENDDO
147          ENDDO          ENDDO
148           ENDIF
149        ENDIF        ENDIF
150    #endif
151    
152  C     Step forward meridional velocity (store in Gv)  C     Step forward meridional velocity (store in Gv)
       psFac = pfFacMom*(1. _d 0 - implicSurfPress)  
153        DO j=jMin,jMax        DO j=jMin,jMax
154          DO i=iMin,iMax          DO i=iMin,iMax
155            gVNm1(i,j,k,bi,bj)=vVel(i,j,k,bi,bj) + deltaTmom * (            gVNm1(i,j,k,bi,bj) = vVel(i,j,k,bi,bj)
156       &       ab15*gV(i,j,k,bi,bj)+ab05*gVNm1(i,j,k,bi,bj)       &     +deltaTmom*(
157       &     - psFac*phiSurfY(i,j)       &         gVtmp(i,j)
158  #ifdef INCLUDE_CD_CODE       &       - psFac*phiSurfY(i,j)
159       &     + gvCD(i,j,k,bi,bj)       &       - phyFac*dPhiHydY(i,j)
 #endif  
160       &        )*_maskS(i,j,k,bi,bj)       &        )*_maskS(i,j,k,bi,bj)
161          ENDDO          ENDDO
162        ENDDO        ENDDO
163    
       IF (staggerTimeStep) THEN  
 C--   -grad Phi_Hyd has not been incorporated to gV and is added here:  
         phyFac = pfFacMom*deltaTmom  
         DO j=jMin,jMax  
           DO i=iMin,iMax  
             gVNm1(i,j,k,bi,bj)=gVNm1(i,j,k,bi,bj)  
      &       - _recip_dyC(i,j,bi,bj)  
      &         *(phiHyd(i,j,k)-phiHyd(i,j-1,k))*phyFac                        
      &         *_maskS(i,j,k,bi,bj)  
           ENDDO  
         ENDDO  
       ENDIF  
   
164        RETURN        RETURN
165        END        END
Legend:
Removed from v.1.23  
changed lines
  Added in v.1.31
  ViewVC Help
Powered by ViewVC 1.1.22