/[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.25 by adcroft, Tue Jul 31 15:01:33 2001 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "CPP_EEOPTIONS.h"  #include "CPP_OPTIONS.h"
5    
6  C     /==========================================================\  C     /==========================================================\
7  C     | S/R TIMESTEP                                             |  C     | S/R TIMESTEP                                             |
8  C     | o Step model fields forward in time                      |  C     | o Step model fields forward in time                      |
9  C     \==========================================================/  C     \==========================================================/
10        SUBROUTINE TIMESTEP( bi, bj, iMin, iMax, jMin, jMax, myThid )        SUBROUTINE TIMESTEP( bi, bj, iMin, iMax, jMin, jMax, K,
11        implicit none       I                     phiHyd, phiSurfX, phiSurfY,
12  ! Common       I                     myIter, myThid )
13          IMPLICIT NONE
14    
15    C     == Global variables ==
16  #include "SIZE.h"  #include "SIZE.h"
17  #include "DYNVARS.h"  #include "DYNVARS.h"
18    #include "EEPARAMS.h"
19  #include "PARAMS.h"  #include "PARAMS.h"
20  #include "GRID.h"  #include "GRID.h"
21  #include "EEPARAMS.h"  
 #include "CG2D.h"  
22  C     == Routine Arguments ==  C     == Routine Arguments ==
23    C     phiHyd    - Hydrostatic Potential (ocean: pressure/rho)
24    C                                       (atmos: geopotentiel)
25    C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)
26    C     phiSurfY             or geopotentiel (atmos) in X and Y direction
27        INTEGER bi,bj,iMin,iMax,jMin,jMax        INTEGER bi,bj,iMin,iMax,jMin,jMax
28        INTEGER myThid        INTEGER K
29          _RL     phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
30          _RL     phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
31          _RL     phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
32          INTEGER myIter, myThid
33    
34  C     == Local variables ==  C     == Local variables ==
35  C     pg       - Pressure gradient terms. Note cg2d_x        INTEGER i,j
36  C                holds term in units so that lateral        _RL ab15,ab05
37  C                gradient is all that is needed.        _RL phxFac,phyFac, psFac
       INTEGER i,j,k  
       REAL ab15,ab05  
       _RL  pg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
38    
39  C     Adams-Bashforth timestepping weights  C     Adams-Bashforth timestepping weights
40        ab15=1.5+abeps        IF (myIter .EQ. 0) THEN
41        ab05=-0.5-abeps         ab15=1.0
42           ab05=0.0
43          ELSE
44           ab15=1.5+abeps
45           ab05=-0.5-abeps
46          ENDIF
47    
48  C     Zonal pressure term  C     Step forward zonal velocity (store in Gu)
49          psFac = pfFacMom*(1. _d 0 - implicSurfPress)
50        DO j=jMin,jMax        DO j=jMin,jMax
        DO i=iMin,iMax  
         pg(i,j)=rDxC(i,j,bi,bj)*  
      &   (cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))  
      &   *gravity*rhonil  
        ENDDO  
       ENDDO  
 C     Step forward zonal velocity  
       DO k=1,Nz  
        DO j=jMin,jMax  
51          DO i=iMin,iMax          DO i=iMin,iMax
52           uVel(i,j,k,bi,bj)=uVel(i,j,k,bi,bj)            gUNm1(i,j,k,bi,bj)=uVel(i,j,k,bi,bj) + deltaTmom * (
53       &    +deltaTmom*(ab15*gU(i,j,k,bi,bj)+ab05*gUNm1(i,j,k,bi,bj)       &       ab15*gU(i,j,k,bi,bj)+ab05*gUNm1(i,j,k,bi,bj)
54       &                -pg(i,j)/rhonil       &     - psFac*phiSurfX(i,j)
55       &    )*maskW(i,j,k,bi,bj)  #ifdef INCLUDE_CD_CODE
56           gUNm1(i,j,k,bi,bj)=gU(i,j,k,bi,bj)       &     + guCD(i,j,k,bi,bj)
57    #endif
58         &        )*_maskW(i,j,k,bi,bj)
59          ENDDO          ENDDO
        ENDDO  
       ENDDO  
 C     Meridional pressure term  
       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  
60        ENDDO        ENDDO
61  C     Step forward meridional velocity  
62        DO k=1,Nz        IF (staggerTimeStep) THEN
63         DO j=jMin,jMax  C--   -grad Phi_Hyd has not been incorporated to gU and is added here:
64          DO i=iMin,iMax          phxFac = pfFacMom*deltaTmom
65           vVel(i,j,k,bi,bj)=vVel(i,j,k,bi,bj)          DO j=jMin,jMax
66       &    +deltaTmom*(ab15*gV(i,j,k,bi,bj)+ab05*gVNm1(i,j,k,bi,bj)            DO i=iMin,iMax
67       &                -pg(i,j)/rhonil              gUNm1(i,j,k,bi,bj)=gUNm1(i,j,k,bi,bj)
68       &               )*maskS(i,j,k,bi,bj)       &       - _recip_dxC(i,j,bi,bj)
69           gVNm1(i,j,k,bi,bj)=gV(i,j,k,bi,bj)       &         *(phiHyd(i,j,k)-phiHyd(i-1,j,k))*phxFac
70         &         *_maskW(i,j,k,bi,bj)
71              ENDDO
72          ENDDO          ENDDO
73         ENDDO        ENDIF
74        ENDDO  
75  C     Step forward temperature  C     Step forward meridional velocity (store in Gv)
76        DO k=1,Nz        psFac = pfFacMom*(1. _d 0 - implicSurfPress)
77         DO j=jMin,jMax        DO j=jMin,jMax
78          DO i=iMin,iMax          DO i=iMin,iMax
79           theta(i,j,k,bi,bj)=theta(i,j,k,bi,bj)            gVNm1(i,j,k,bi,bj)=vVel(i,j,k,bi,bj) + deltaTmom * (
80       &    +deltaTtracer*(ab15*gT(i,j,k,bi,bj)+ab05*gTNm1(i,j,k,bi,bj))       &       ab15*gV(i,j,k,bi,bj)+ab05*gVNm1(i,j,k,bi,bj)
81           gTNm1(i,j,k,bi,bj)=gT(i,j,k,bi,bj)       &     - psFac*phiSurfY(i,j)
82    #ifdef INCLUDE_CD_CODE
83         &     + gvCD(i,j,k,bi,bj)
84    #endif
85         &        )*_maskS(i,j,k,bi,bj)
86          ENDDO          ENDDO
        ENDDO  
87        ENDDO        ENDDO
88    
89        _BARRIER        IF (staggerTimeStep) THEN
90  C     CALL PLOT_FIELD_XYZR8( uVel, 'TIEMSTEP.1 uVel',Nz,1,myThid)  C--   -grad Phi_Hyd has not been incorporated to gV and is added here:
91  C     CALL PLOT_FIELD_XYZR8( vVel, 'TIEMSTEP.1 vVel',Nz,1,myThid)          phyFac = pfFacMom*deltaTmom
92  C     CALL PLOT_FIELD_XYZR8( theta, 'TIEMSTEP.1 theta',Nz,1,myThid)          DO j=jMin,jMax
93              DO i=iMin,iMax
94                gVNm1(i,j,k,bi,bj)=gVNm1(i,j,k,bi,bj)
95         &       - _recip_dyC(i,j,bi,bj)
96         &         *(phiHyd(i,j,k)-phiHyd(i,j-1,k))*phyFac
97         &         *_maskS(i,j,k,bi,bj)
98              ENDDO
99            ENDDO
100          ENDIF
101    
102        RETURN        RETURN
103        END        END
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.25
  ViewVC Help
Powered by ViewVC 1.1.22