/[MITgcm]/MITgcm/model/src/timestep.F
ViewVC logotype

Contents of /MITgcm/model/src/timestep.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph

Revision 1.30 - (show annotations) (download)
Sun Feb 9 02:00:50 2003 UTC (21 years, 4 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint48e_post, checkpoint48d_post, checkpoint48f_post
Changes since 1.29: +7 -9 lines 
in preparation for r*:
a) use pre-computed gradient of hydrostatic potential:
   changes in timestep.F & mom_cdscheme.F affect results of ideal_2D_oce
b) move phi0surf from calc_phi_hyd to calc_grad_phi_hyd :
   => affects results of glob_oce_pressure (different truncation error)
1 C $Header: /u/gcmpack/MITgcm/model/src/timestep.F,v 1.29 2003/02/08 02:09:20 jmc Exp $
2 C $Name: $
3
4 #include "CPP_OPTIONS.h"
5
6 CBOP
7 C !ROUTINE: TIMESTEP
8 C !INTERFACE:
9 SUBROUTINE TIMESTEP( bi, bj, iMin, iMax, jMin, jMax, K,
10 I phiHyd, dPhiHydX,dPhiHydY, phiSurfX, phiSurfY,
11 I myIter, myThid )
12 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"
23 #include "DYNVARS.h"
24 #include "EEPARAMS.h"
25 #include "PARAMS.h"
26 #include "GRID.h"
27 #include "SURFACE.h"
28
29 C !INPUT/OUTPUT PARAMETERS:
30 C == Routine Arguments ==
31 C phiHyd :: Hydrostatic Potential (ocean: pressure/rho)
32 C (atmos: geopotential)
33 C dPhiHydX,Y :: Gradient (X & Y directions) of Hydrostatic Potential
34 C phiSurfX :: gradient of Surface potential (Pressure/rho, ocean)
35 C phiSurfY :: or geopotential (atmos) in X and Y direction
36 INTEGER bi,bj,iMin,iMax,jMin,jMax
37 INTEGER K
38 _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
39 _RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
40 _RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
41 _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
42 _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
43 INTEGER myIter, myThid
44
45 C !LOCAL VARIABLES:
46 C == Local variables ==
47 INTEGER i,j
48 _RL ab15,ab05
49 _RL phxFac,phyFac, psFac
50 _RL gUtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
51 _RL gVtmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
52 CEOP
53
54 C Adams-Bashforth timestepping weights
55 IF (myIter .EQ. 0) THEN
56 ab15=1.0
57 ab05=0.0
58 ELSE
59 ab15=1.5+abeps
60 ab05=-0.5-abeps
61 ENDIF
62
63 C-- stagger time step: grad Phi_Hyp is not in gU,gV => add it in this S/R
64 IF (staggerTimeStep) THEN
65 c IF (.FALSE.) THEN
66 phxFac = pfFacMom
67 phyFac = pfFacMom
68 ELSE
69 phxFac = 0.
70 phyFac = 0.
71 ENDIF
72
73 C-- explicit part of the surface potential gradient is added in this S/R
74 psFac = pfFacMom*(1. _d 0 - implicSurfPress)
75
76 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
77 C- Compute effective gU term (including Adams-Bashforth weights) :
78 DO j=jMin,jMax
79 DO i=iMin,iMax
80 gUtmp(i,j) = ab15*gU(i,j,k,bi,bj)
81 & + ab05*gUNm1(i,j,k,bi,bj)
82 #ifdef INCLUDE_CD_CODE
83 & + guCD(i,j,k,bi,bj)
84 #endif
85 ENDDO
86 ENDDO
87
88 #ifdef NONLIN_FRSURF
89 IF (.NOT. vectorInvariantMomentum
90 & .AND. nonlinFreeSurf.GT.1) THEN
91 IF (select_rStar.GT.0) THEN
92 DO j=jMin,jMax
93 DO i=iMin,iMax
94 gUtmp(i,j) = gUtmp(i,j)/rStarExpW(i,j,bi,bj)
95 ENDDO
96 ENDDO
97 ELSE
98 DO j=jMin,jMax
99 DO i=iMin,iMax
100 IF ( k.EQ.ksurfW(i,j,bi,bj) ) THEN
101 gUtmp(i,j) = gUtmp(i,j)
102 & *hFacW(i,j,k,bi,bj)/hFac_surfW(i,j,bi,bj)
103 ENDIF
104 ENDDO
105 ENDDO
106 ENDIF
107 ENDIF
108 #endif
109
110 C Step forward zonal velocity (store in Gu)
111 DO j=jMin,jMax
112 DO i=iMin,iMax
113 gUNm1(i,j,k,bi,bj) = uVel(i,j,k,bi,bj)
114 & +deltaTmom*(
115 & gUtmp(i,j)
116 & - psFac*phiSurfX(i,j)
117 & - phxFac*dPhiHydX(i,j)
118 & )*_maskW(i,j,k,bi,bj)
119 ENDDO
120 ENDDO
121
122 IF (.FALSE.) THEN
123 c IF (staggerTimeStep) THEN
124 C-- -grad Phi_Hyd has not been incorporated to gU and is added here:
125 phxFac = pfFacMom*deltaTmom
126 DO j=jMin,jMax
127 DO i=iMin,iMax
128 gUNm1(i,j,k,bi,bj)=gUNm1(i,j,k,bi,bj)
129 & - _recip_dxC(i,j,bi,bj)
130 & *(phiHyd(i,j,k)-phiHyd(i-1,j,k))*phxFac
131 & *_maskW(i,j,k,bi,bj)
132 ENDDO
133 ENDDO
134 ENDIF
135
136 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
137 C- Compute effective gV term (including Adams-Bashforth weights) :
138 DO j=jMin,jMax
139 DO i=iMin,iMax
140 gVtmp(i,j) = ab15*gV(i,j,k,bi,bj)
141 & + ab05*gVNm1(i,j,k,bi,bj)
142 #ifdef INCLUDE_CD_CODE
143 & + gvCD(i,j,k,bi,bj)
144 #endif
145 ENDDO
146 ENDDO
147
148 #ifdef NONLIN_FRSURF
149 IF (.NOT. vectorInvariantMomentum
150 & .AND. nonlinFreeSurf.GT.1) THEN
151 IF (select_rStar.GT.0) THEN
152 DO j=jMin,jMax
153 DO i=iMin,iMax
154 gVtmp(i,j) = gVtmp(i,j)/rStarExpS(i,j,bi,bj)
155 ENDDO
156 ENDDO
157 ELSE
158 DO j=jMin,jMax
159 DO i=iMin,iMax
160 IF ( k.EQ.ksurfS(i,j,bi,bj) ) THEN
161 gVtmp(i,j) = gVtmp(i,j)
162 & *hFacS(i,j,k,bi,bj)/hFac_surfS(i,j,bi,bj)
163 ENDIF
164 ENDDO
165 ENDDO
166 ENDIF
167 ENDIF
168 #endif
169
170 C Step forward meridional velocity (store in Gv)
171 DO j=jMin,jMax
172 DO i=iMin,iMax
173 gVNm1(i,j,k,bi,bj) = vVel(i,j,k,bi,bj)
174 & +deltaTmom*(
175 & gVtmp(i,j)
176 & - psFac*phiSurfY(i,j)
177 & - phyFac*dPhiHydY(i,j)
178 & )*_maskS(i,j,k,bi,bj)
179 ENDDO
180 ENDDO
181
182 IF (.FALSE.) THEN
183 c IF (staggerTimeStep) THEN
184 C-- -grad Phi_Hyd has not been incorporated to gV and is added here:
185 phyFac = pfFacMom*deltaTmom
186 DO j=jMin,jMax
187 DO i=iMin,iMax
188 gVNm1(i,j,k,bi,bj)=gVNm1(i,j,k,bi,bj)
189 & - _recip_dyC(i,j,bi,bj)
190 & *(phiHyd(i,j,k)-phiHyd(i,j-1,k))*phyFac
191 & *_maskS(i,j,k,bi,bj)
192 ENDDO
193 ENDDO
194 ENDIF
195
196 RETURN
197 END
  ViewVC Help
Powered by ViewVC 1.1.22