/[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.31 by jmc, Tue Feb 18 15:27:25 2003 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     /==========================================================\  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, myThid )       I                     dPhiHydX,dPhiHydY, phiSurfX, phiSurfY,
11        implicit none       I                     myIter, myThid )
12  ! Common  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 ==
22  #include "SIZE.h"  #include "SIZE.h"
23  #include "DYNVARS.h"  #include "DYNVARS.h"
24    #include "EEPARAMS.h"
25  #include "PARAMS.h"  #include "PARAMS.h"
26  #include "GRID.h"  #include "GRID.h"
27  #include "EEPARAMS.h"  #include "SURFACE.h"
28  #include "CG2D.h"  
29    C     !INPUT/OUTPUT PARAMETERS:
30  C     == Routine Arguments ==  C     == Routine Arguments ==
31    C     dPhiHydX,Y :: Gradient (X & Y directions) of Hydrostatic Potential
32    C     phiSurfX :: gradient of Surface potential (Pressure/rho, ocean)
33    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 myThid        INTEGER K
36          _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)
39          _RL     phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
40          INTEGER myIter, myThid
41    
42    C     !LOCAL VARIABLES:
43  C     == Local variables ==  C     == Local variables ==
44  C     pg       - Pressure gradient terms. Note cg2d_x        INTEGER i,j
45  C                holds term in units so that lateral        _RL ab15,ab05
46  C                gradient is all that is needed.        _RL phxFac,phyFac, psFac
47        INTEGER i,j,k        _RL     gUtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
48        REAL ab15,ab05        _RL     gVtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
49        _RL  pg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  CEOP
50    
51  C     Adams-Bashforth timestepping weights  C     Adams-Bashforth timestepping weights
52        ab15=1.5+abeps        IF (myIter .EQ. 0) THEN
53        ab05=-0.5-abeps         ab15=1.0
54           ab05=0.0
55          ELSE
56           ab15=1.5+abeps
57           ab05=-0.5-abeps
58          ENDIF
59    
60    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)
71    
72  C     Zonal pressure term  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          pg(i,j)=rDxC(i,j,bi,bj)*          gUtmp(i,j) = ab15*gU(i,j,k,bi,bj)
77       &   (cg2d_x(i,j,bi,bj)-cg2d_x(i-1,j,bi,bj))       &             + ab05*gUNm1(i,j,k,bi,bj)  
78       &   *gravity*rhonil  #ifdef INCLUDE_CD_CODE
79         &             + guCD(i,j,k,bi,bj)
80    #endif
81         ENDDO         ENDDO
82        ENDDO        ENDDO
83  C     Step forward zonal velocity        
84        DO k=1,Nz  #ifdef NONLIN_FRSURF
85         DO j=jMin,jMax        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          DO i=iMin,iMax
109           uVel(i,j,k,bi,bj)=uVel(i,j,k,bi,bj)            gUNm1(i,j,k,bi,bj) = uVel(i,j,k,bi,bj)
110       &    +deltaTmom*(ab15*gU(i,j,k,bi,bj)+ab05*gUNm1(i,j,k,bi,bj)       &     +deltaTmom*(
111       &                -pg(i,j)/rhonil       &         gUtmp(i,j)
112       &    )*maskW(i,j,k,bi,bj)       &       - psFac*phiSurfX(i,j)
113           gUNm1(i,j,k,bi,bj)=gU(i,j,k,bi,bj)       &       - phxFac*dPhiHydX(i,j)
114         &        )*_maskW(i,j,k,bi,bj)
115          ENDDO          ENDDO
        ENDDO  
116        ENDDO        ENDDO
117  C     Meridional pressure term  
118    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
119    C-    Compute effective gV term (including Adams-Bashforth weights) :
120        DO j=jMin,jMax        DO j=jMin,jMax
121         DO i=iMin,iMax         DO i=iMin,iMax
122          pg(i,j)=rDyC(i,j,bi,bj)*          gVtmp(i,j) = ab15*gV(i,j,k,bi,bj)
123       &   (cg2d_x(i,j,bi,bj)-cg2d_x(i,j-1,bi,bj))       &             + ab05*gVNm1(i,j,k,bi,bj)  
124       &   *gravity*rhonil  #ifdef INCLUDE_CD_CODE
125         &             + gvCD(i,j,k,bi,bj)
126    #endif
127         ENDDO         ENDDO
128        ENDDO        ENDDO
129  C     Step forward meridional velocity        
130        DO k=1,Nz  #ifdef NONLIN_FRSURF
131         DO j=jMin,jMax        IF (.NOT. vectorInvariantMomentum
132          DO i=iMin,iMax       &    .AND. nonlinFreeSurf.GT.1) THEN
133           vVel(i,j,k,bi,bj)=vVel(i,j,k,bi,bj)         IF (select_rStar.GT.0) THEN
134       &    +deltaTmom*(ab15*gV(i,j,k,bi,bj)+ab05*gVNm1(i,j,k,bi,bj)          DO j=jMin,jMax
135       &                -pg(i,j)/rhonil           DO i=iMin,iMax
136       &               )*maskS(i,j,k,bi,bj)             gVtmp(i,j) = gVtmp(i,j)/rStarExpS(i,j,bi,bj)
137           gVNm1(i,j,k,bi,bj)=gV(i,j,k,bi,bj)           ENDDO
138          ENDDO          ENDDO
139         ENDDO         ELSE
140        ENDDO          DO j=jMin,jMax
141  C     Step forward temperature           DO i=iMin,iMax
142        DO k=1,Nz            IF ( k.EQ.ksurfS(i,j,bi,bj) ) THEN
143         DO j=jMin,jMax             gVtmp(i,j) = gVtmp(i,j)
144         &         *hFacS(i,j,k,bi,bj)/hFac_surfS(i,j,bi,bj)
145              ENDIF
146             ENDDO
147            ENDDO
148           ENDIF
149          ENDIF
150    #endif
151    
152    C     Step forward meridional velocity (store in Gv)
153          DO j=jMin,jMax
154          DO i=iMin,iMax          DO i=iMin,iMax
155           theta(i,j,k,bi,bj)=theta(i,j,k,bi,bj)            gVNm1(i,j,k,bi,bj) = vVel(i,j,k,bi,bj)
156       &    +deltaTtracer*(ab15*gT(i,j,k,bi,bj)+ab05*gTNm1(i,j,k,bi,bj))       &     +deltaTmom*(
157           gTNm1(i,j,k,bi,bj)=gT(i,j,k,bi,bj)       &         gVtmp(i,j)
158         &       - psFac*phiSurfY(i,j)
159         &       - phyFac*dPhiHydY(i,j)
160         &        )*_maskS(i,j,k,bi,bj)
161          ENDDO          ENDDO
        ENDDO  
162        ENDDO        ENDDO
163    
       _BARRIER  
 C     CALL PLOT_FIELD_XYZR8( uVel, 'TIEMSTEP.1 uVel',Nz,1,myThid)  
 C     CALL PLOT_FIELD_XYZR8( vVel, 'TIEMSTEP.1 vVel',Nz,1,myThid)  
 C     CALL PLOT_FIELD_XYZR8( theta, 'TIEMSTEP.1 theta',Nz,1,myThid)  
   
   
164        RETURN        RETURN
165        END        END
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.31
  ViewVC Help
Powered by ViewVC 1.1.22