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

Annotation of /MITgcm/model/src/dynamics.F

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


Revision 1.11 - (hide annotations) (download)
Mon Jun 1 20:36:13 1998 UTC (25 years, 11 months ago) by adcroft
Branch: MAIN
Changes since 1.10: +74 -15 lines
Swapped the time-stepping algorithm around (just a little bit).
We now officially use the predictor-corrector terminology.
We make the prediction step at the end of the dynamics() section
and store the result in the gUNm1, gVNm1, gTNm1 arrays.
The "tricky" part is that at the beginning of the dynamics section,
where the "correction" is made, theses arrays must be initialised
at the beginning of any run. A new routine init_predictor() does this.
This is "all" in preparation for implicit diffusion. Let's hope
it's enough...

1 adcroft 1.11 C $Header: /u/gcmpack/models/MITgcmUV/model/src/dynamics.F,v 1.10 1998/05/28 16:19:50 adcroft Exp $
2 cnh 1.1
3     #include "CPP_EEOPTIONS.h"
4    
5 cnh 1.8 SUBROUTINE DYNAMICS(myTime, myIter, myThid)
6 cnh 1.1 C /==========================================================\
7     C | SUBROUTINE DYNAMICS |
8     C | o Controlling routine for the explicit part of the model |
9     C | dynamics. |
10     C |==========================================================|
11     C | This routine evaluates the "dynamics" terms for each |
12     C | block of ocean in turn. Because the blocks of ocean have |
13     C | overlap regions they are independent of one another. |
14     C | If terms involving lateral integrals are needed in this |
15     C | routine care will be needed. Similarly finite-difference |
16     C | operations with stencils wider than the overlap region |
17     C | require special consideration. |
18     C | Notes |
19     C | ===== |
20     C | C*P* comments indicating place holders for which code is |
21     C | presently being developed. |
22     C \==========================================================/
23    
24     C == Global variables ===
25     #include "SIZE.h"
26     #include "EEPARAMS.h"
27     #include "CG2D.h"
28 adcroft 1.6 #include "PARAMS.h"
29 adcroft 1.3 #include "DYNVARS.h"
30 cnh 1.1
31     C == Routine arguments ==
32 cnh 1.8 C myTime - Current time in simulation
33     C myIter - Current iteration number in simulation
34 cnh 1.1 C myThid - Thread number for this instance of the routine.
35     INTEGER myThid
36 cnh 1.8 _RL myTime
37     INTEGER myIter
38 cnh 1.1
39     C == Local variables
40     C xA, yA - Per block temporaries holding face areas
41     C uTrans, vTrans, wTrans - Per block temporaries holding flow transport
42     C o uTrans: Zonal transport
43     C o vTrans: Meridional transport
44     C o wTrans: Vertical transport
45     C maskC,maskUp o maskC: land/water mask for tracer cells
46     C o maskUp: land/water mask for W points
47     C aTerm, xTerm, cTerm - Work arrays for holding separate terms in
48     C mTerm, pTerm, tendency equations.
49     C fZon, fMer, fVer[STUV] o aTerm: Advection term
50     C o xTerm: Mixing term
51     C o cTerm: Coriolis term
52     C o mTerm: Metric term
53     C o pTerm: Pressure term
54     C o fZon: Zonal flux term
55     C o fMer: Meridional flux term
56     C o fVer: Vertical flux term - note fVer
57     C is "pipelined" in the vertical
58     C so we need an fVer for each
59     C variable.
60     C iMin, iMax - Ranges and sub-block indices on which calculations
61     C jMin, jMax are applied.
62     C bi, bj
63     C k, kUp, kDown, kM1 - Index for layer above and below. kUp and kDown
64     C are switched with layer to be the appropriate index
65     C into fVerTerm
66     _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
67     _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
68     _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
69     _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
70     _RL wTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
71     _RS maskC (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
72     _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
73     _RL aTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
74     _RL xTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
75     _RL cTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
76     _RL mTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
77     _RL pTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
78     _RL fZon (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
79     _RL fMer (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
80     _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
81     _RL fVerS (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
82     _RL fVerU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
83     _RL fVerV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
84     _RL pH (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz)
85 adcroft 1.3 _RL rhokm1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
86     _RL rhokp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
87 adcroft 1.10 _RL rhotmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
88 adcroft 1.4 _RL pSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
89     _RL pSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
90 adcroft 1.6 _RL K13 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz)
91     _RL K23 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz)
92     _RL K33 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz)
93     _RL KapGM (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
94 cnh 1.1 INTEGER iMin, iMax
95     INTEGER jMin, jMax
96     INTEGER bi, bj
97     INTEGER i, j
98     INTEGER k, kM1, kUp, kDown
99    
100 adcroft 1.11 C--- The algorithm...
101     C
102     C "Correction Step"
103     C =================
104     C Here we update the horizontal velocities with the surface
105     C pressure such that the resulting flow is either consistent
106     C with the free-surface evolution or the rigid-lid:
107     C U[n] = U* + dt x d/dx P
108     C V[n] = V* + dt x d/dy P
109     C With implicit diffusion, the tracers must also be "finalized"
110     C (1 + dt * K * d_zz) theta[n] = theta*
111     C (1 + dt * K * d_zz) salt[n] = salt*
112     C
113     C "Calculation of Gs"
114     C ===================
115     C This is where all the accelerations and tendencies (ie.
116     C physics, parameterizations etc...) are calculated
117     C w = sum_z ( div. u[n] )
118     C rho = rho ( theta[n], salt[n] )
119     C K31 = K31 ( rho )
120     C Gu[n] = Gu( u[n], v[n], w, rho, Ph, ... )
121     C Gv[n] = Gv( u[n], v[n], w, rho, Ph, ... )
122     C Gt[n] = Gt( theta[n], u[n], v[n], w, K31, ... )
123     C Gs[n] = Gs( salt[n], u[n], v[n], w, K31, ... )
124     C
125     C "Time-stepping" or "Predicition"
126     C ================================
127     C The models variables are stepped forward with the appropriate
128     C time-stepping scheme (currently we use Adams-Bashforth II)
129     C - For momentum, the result is always *only* a "prediction"
130     C in that the flow may be divergent and will be "corrected"
131     C later with a surface pressure gradient.
132     C - Normally for tracers the result is the new field at time
133     C level [n+1} *BUT* in the case of implicit diffusion the result
134     C is also *only* a prediction.
135     C - We denote "predictors" with an asterisk (*).
136     C U* = U[n] + dt x ( 3/2 Gu[n] - 1/2 Gu[n-1] )
137     C V* = V[n] + dt x ( 3/2 Gv[n] - 1/2 Gv[n-1] )
138     C theta[n+1] = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
139     C salt[n+1] = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
140     C or with implicit diffusion
141     C theta* = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
142     C
143     C salt* = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
144     C---
145    
146    
147 cnh 1.1 C-- Set up work arrays with valid (i.e. not NaN) values
148     C These inital values do not alter the numerical results. They
149     C just ensure that all memory references are to valid floating
150     C point numbers. This prevents spurious hardware signals due to
151     C uninitialised but inert locations.
152     DO j=1-OLy,sNy+OLy
153     DO i=1-OLx,sNx+OLx
154 adcroft 1.5 xA(i,j) = 0. _d 0
155     yA(i,j) = 0. _d 0
156     uTrans(i,j) = 0. _d 0
157     vTrans(i,j) = 0. _d 0
158     aTerm(i,j) = 0. _d 0
159     xTerm(i,j) = 0. _d 0
160     cTerm(i,j) = 0. _d 0
161     mTerm(i,j) = 0. _d 0
162     pTerm(i,j) = 0. _d 0
163     fZon(i,j) = 0. _d 0
164     fMer(i,j) = 0. _d 0
165 cnh 1.1 DO K=1,nZ
166 adcroft 1.5 pH (i,j,k) = 0. _d 0
167 adcroft 1.6 K13(i,j,k) = 0. _d 0
168     K23(i,j,k) = 0. _d 0
169     K33(i,j,k) = 0. _d 0
170 cnh 1.1 ENDDO
171 adcroft 1.5 rhokm1(i,j) = 0. _d 0
172     rhokp1(i,j) = 0. _d 0
173 adcroft 1.10 rhotmp(i,j) = 0. _d 0
174 cnh 1.1 ENDDO
175     ENDDO
176    
177     DO bj=myByLo(myThid),myByHi(myThid)
178     DO bi=myBxLo(myThid),myBxHi(myThid)
179    
180 cnh 1.7 C-- Set up work arrays that need valid initial values
181     DO j=1-OLy,sNy+OLy
182     DO i=1-OLx,sNx+OLx
183     wTrans(i,j) = 0. _d 0
184     fVerT(i,j,1) = 0. _d 0
185     fVerT(i,j,2) = 0. _d 0
186     fVerS(i,j,1) = 0. _d 0
187     fVerS(i,j,2) = 0. _d 0
188     fVerU(i,j,1) = 0. _d 0
189     fVerU(i,j,2) = 0. _d 0
190     fVerV(i,j,1) = 0. _d 0
191     fVerV(i,j,2) = 0. _d 0
192 adcroft 1.11 pH(i,j,1) = 0. _d 0
193     K13(i,j,1) = 0. _d 0
194     K23(i,j,1) = 0. _d 0
195     K33(i,j,1) = 0. _d 0
196     KapGM(i,j) = 0. _d 0
197 cnh 1.7 ENDDO
198     ENDDO
199    
200 cnh 1.1 iMin = 1-OLx+1
201     iMax = sNx+OLx
202     jMin = 1-OLy+1
203     jMax = sNy+OLy
204    
205 adcroft 1.4 C-- Calculate gradient of surface pressure
206     CALL GRAD_PSURF(
207     I bi,bj,iMin,iMax,jMin,jMax,
208     O pSurfX,pSurfY,
209     I myThid)
210    
211     C-- Update fields in top level according to tendency terms
212 adcroft 1.11 CALL CORRECTION_STEP(
213 adcroft 1.4 I bi,bj,iMin,iMax,jMin,jMax,1,pSurfX,pSurfY,myThid)
214 cnh 1.1
215 cnh 1.7 C-- Density of 1st level (below W(1)) reference to level 1
216     CALL FIND_RHO(
217 adcroft 1.10 I bi, bj, iMin, iMax, jMin, jMax, 1, 1, eosType,
218 cnh 1.7 O rhoKm1,
219     I myThid )
220     C-- Integrate hydrostatic balance for pH with BC of pH(z=0)=0
221     CALL CALC_PH(
222     I bi,bj,iMin,iMax,jMin,jMax,1,rhoKm1,rhoKm1,
223     U pH,
224 adcroft 1.5 I myThid )
225 adcroft 1.10 DO J=1-Oly,sNy+Oly
226     DO I=1-Olx,sNx+Olx
227     rhoKp1(I,J)=rhoKm1(I,J)
228     ENDDO
229     ENDDO
230 adcroft 1.5
231 adcroft 1.3 DO K=2,Nz
232 adcroft 1.4 C-- Update fields in Kth level according to tendency terms
233 adcroft 1.11 CALL CORRECTION_STEP(
234 adcroft 1.4 I bi,bj,iMin,iMax,jMin,jMax,K,pSurfX,pSurfY,myThid)
235 adcroft 1.10 C-- Density of K-1 level (above W(K)) reference to K-1 level
236     copt CALL FIND_RHO(
237     copt I bi, bj, iMin, iMax, jMin, jMax, K-1, K-1, eosType,
238     copt O rhoKm1,
239     copt I myThid )
240     C rhoKm1=rhoKp1
241     DO J=1-Oly,sNy+Oly
242     DO I=1-Olx,sNx+Olx
243     rhoKm1(I,J)=rhoKp1(I,J)
244     ENDDO
245     ENDDO
246     C-- Density of K level (below W(K)) reference to K level
247 cnh 1.7 CALL FIND_RHO(
248 adcroft 1.10 I bi, bj, iMin, iMax, jMin, jMax, K, K, eosType,
249     O rhoKp1,
250 cnh 1.7 I myThid )
251 adcroft 1.10 C-- Density of K-1 level (above W(K)) reference to K level
252 cnh 1.7 CALL FIND_RHO(
253 adcroft 1.10 I bi, bj, iMin, iMax, jMin, jMax, K-1, K, eosType,
254     O rhotmp,
255 cnh 1.7 I myThid )
256 adcroft 1.6 C-- Calculate iso-neutral slopes for the GM/Redi parameterisation
257 cnh 1.7 CALL CALC_ISOSLOPES(
258     I bi, bj, iMin, iMax, jMin, jMax, K,
259 adcroft 1.10 I rhoKm1, rhoKp1, rhotmp,
260 cnh 1.7 O K13, K23, K33, KapGM,
261     I myThid )
262 cnh 1.1 C-- Calculate static stability and mix where convectively unstable
263 cnh 1.7 CALL CONVECT(
264 cnh 1.8 I bi,bj,iMin,iMax,jMin,jMax,K,rhoKm1,rhoKp1,
265     I myTime,myIter,myThid)
266 cnh 1.7 C-- Density of K-1 level (above W(K)) reference to K-1 level
267     CALL FIND_RHO(
268 adcroft 1.10 I bi, bj, iMin, iMax, jMin, jMax, K-1, K-1, eosType,
269 cnh 1.7 O rhoKm1,
270     I myThid )
271     C-- Density of K level (below W(K)) referenced to K level
272     CALL FIND_RHO(
273 adcroft 1.10 I bi, bj, iMin, iMax, jMin, jMax, K, K, eosType,
274 cnh 1.7 O rhoKp1,
275     I myThid )
276     C-- Integrate hydrostatic balance for pH with BC of pH(z=0)=0
277     CALL CALC_PH(
278     I bi,bj,iMin,iMax,jMin,jMax,K,rhoKm1,rhoKp1,
279     U pH,
280 adcroft 1.3 I myThid )
281 cnh 1.1
282 adcroft 1.11 ENDDO ! K
283    
284     C-- Initial boundary condition on barotropic divergence integral
285     DO j=1-OLy,sNy+OLy
286     DO i=1-OLx,sNx+OLx
287     cg2d_b(i,j,bi,bj) = 0. _d 0
288     ENDDO
289 cnh 1.7 ENDDO
290 adcroft 1.5
291 cnh 1.1 DO K = Nz, 1, -1
292     kM1 =max(1,k-1) ! Points to level above k (=k-1)
293     kUp =1+MOD(k+1,2) ! Cycles through 1,2 to point to layer above
294     kDown=1+MOD(k,2) ! Cycles through 2,1 to point to current layer
295     iMin = 1-OLx+2
296     iMax = sNx+OLx-1
297     jMin = 1-OLy+2
298     jMax = sNy+OLy-1
299    
300     C-- Get temporary terms used by tendency routines
301     CALL CALC_COMMON_FACTORS (
302     I bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,
303     O xA,yA,uTrans,vTrans,wTrans,maskC,maskUp,
304     I myThid)
305    
306     C-- Calculate accelerations in the momentum equations
307 cnh 1.9 IF ( momStepping ) THEN
308     CALL CALC_MOM_RHS(
309     I bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,
310     I xA,yA,uTrans,vTrans,wTrans,maskC,
311     I pH,
312     U aTerm,xTerm,cTerm,mTerm,pTerm,
313     U fZon, fMer, fVerU, fVerV,
314     I myThid)
315     ENDIF
316 cnh 1.1
317     C-- Calculate active tracer tendencies
318 cnh 1.9 IF ( tempStepping ) THEN
319     CALL CALC_GT(
320     I bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown,
321     I xA,yA,uTrans,vTrans,wTrans,maskUp,
322     I K13,K23,K33,KapGM,
323     U aTerm,xTerm,fZon,fMer,fVerT,
324     I myThid)
325     ENDIF
326 cnh 1.1 Cdbg CALL CALC_GS(
327     Cdbg I bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown,
328     Cdbg I xA,yA,uTrans,vTrans,wTrans,maskUp,
329 adcroft 1.6 Cdbg I K13,K23,K33,KapGM,
330 cnh 1.1 Cdbg U aTerm,xTerm,fZon,fMer,fVerS,
331     Cdbg I myThid)
332    
333 adcroft 1.11 C-- Prediction step (step forward all model variables)
334     CALL TIMESTEP(
335     I bi,bj,iMin,iMax,jMin,jMax,K,
336     I myThid)
337    
338     C-- Diagnose barotropic divergence of predicted fields
339     CALL DIV_G(
340     I bi,bj,iMin,iMax,jMin,jMax,K,
341     I xA,yA,
342     I myThid)
343    
344     ENDDO ! K
345 cnh 1.1
346     ENDDO
347     ENDDO
348 adcroft 1.6
349     !dbg write(0,*) 'dynamics: pS',minval(cg2d_x),maxval(cg2d_x)
350     !dbg write(0,*) 'dynamics: U',minval(uVel(1:sNx,1:sNy,:,:,:)),
351     !dbg & maxval(uVel(1:sNx,1:sNy,:,:,:))
352     !dbg write(0,*) 'dynamics: V',minval(vVel(1:sNx,1:sNy,:,:,:)),
353     !dbg & maxval(vVel(1:sNx,1:sNy,:,:,:))
354 adcroft 1.10 !dbg write(0,*) 'dynamics: K13',minval(K13(1:sNx,1:sNy,:)),
355     !dbg & maxval(K13(1:sNx,1:sNy,:))
356     !dbg write(0,*) 'dynamics: K23',minval(K23(1:sNx,1:sNy,:)),
357     !dbg & maxval(K23(1:sNx,1:sNy,:))
358     !dbg write(0,*) 'dynamics: K33',minval(K33(1:sNx,1:sNy,:)),
359     !dbg & maxval(K33(1:sNx,1:sNy,:))
360 adcroft 1.6 !dbg write(0,*) 'dynamics: gT',minval(gT(1:sNx,1:sNy,:,:,:)),
361     !dbg & maxval(gT(1:sNx,1:sNy,:,:,:))
362     !dbg write(0,*) 'dynamics: T',minval(Theta(1:sNx,1:sNy,:,:,:)),
363     !dbg & maxval(Theta(1:sNx,1:sNy,:,:,:))
364     !dbg write(0,*) 'dynamics: pH',minval(pH/(Gravity*Rhonil)),
365     !dbg & maxval(pH/(Gravity*Rhonil))
366 cnh 1.1
367     RETURN
368     END

  ViewVC Help
Powered by ViewVC 1.1.22