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

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

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


Revision 1.24 - (show annotations) (download)
Wed Jul 1 19:57:22 1998 UTC (25 years, 10 months ago) by adcroft
Branch: MAIN
CVS Tags: checkpoint11, checkpoint10, checkpoint12, branch-point-rdot
Branch point for: branch-rdot
Changes since 1.23: +2 -2 lines
Changed the #include at the top from CPP_EEOPTIONS.h to CPP_OPTIONS.h
This ought to be done through-out (?) but was necessary in at least
the_model_main.F, dynamics.F, ini_parms.F because of the new
macro ALLOW_DIAGNOSTICS.

1 C $Header: /u/gcmpack/models/MITgcmUV/model/src/dynamics.F,v 1.23 1998/07/01 19:49:36 adcroft Exp $
2
3 #include "CPP_OPTIONS.h"
4
5 SUBROUTINE DYNAMICS(myTime, myIter, myThid)
6 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 #include "PARAMS.h"
29 #include "DYNVARS.h"
30
31 C == Routine arguments ==
32 C myTime - Current time in simulation
33 C myIter - Current iteration number in simulation
34 C myThid - Thread number for this instance of the routine.
35 INTEGER myThid
36 _RL myTime
37 INTEGER myIter
38
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 wVel o uTrans: Zonal transport
43 C o vTrans: Meridional transport
44 C o wTrans: Vertical transport
45 C o wVel: Vertical velocity at upper and lower
46 C cell faces.
47 C maskC,maskUp o maskC: land/water mask for tracer cells
48 C o maskUp: land/water mask for W points
49 C aTerm, xTerm, cTerm - Work arrays for holding separate terms in
50 C mTerm, pTerm, tendency equations.
51 C fZon, fMer, fVer[STUV] o aTerm: Advection term
52 C o xTerm: Mixing term
53 C o cTerm: Coriolis term
54 C o mTerm: Metric term
55 C o pTerm: Pressure term
56 C o fZon: Zonal flux term
57 C o fMer: Meridional flux term
58 C o fVer: Vertical flux term - note fVer
59 C is "pipelined" in the vertical
60 C so we need an fVer for each
61 C variable.
62 C iMin, iMax - Ranges and sub-block indices on which calculations
63 C jMin, jMax are applied.
64 C bi, bj
65 C k, kUp, kDown, kM1 - Index for layer above and below. kUp and kDown
66 C are switched with layer to be the appropriate index
67 C into fVerTerm
68 _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
69 _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
70 _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
71 _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
72 _RL wTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
73 _RL wVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
74 _RS maskC (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
75 _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
76 _RL aTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
77 _RL xTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
78 _RL cTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
79 _RL mTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
80 _RL pTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
81 _RL fZon (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
82 _RL fMer (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
83 _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
84 _RL fVerS (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
85 _RL fVerU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
86 _RL fVerV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
87 _RL pH (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz)
88 _RL rhokm1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
89 _RL rhokp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
90 _RL rhok (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
91 _RL rhotmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
92 _RL pSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
93 _RL pSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
94 _RL K13 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz)
95 _RL K23 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz)
96 _RL K33 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz)
97 _RL KapGM (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98 _RL KappaZT(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nz)
99 _RL KappaZS(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nz)
100
101 INTEGER iMin, iMax
102 INTEGER jMin, jMax
103 INTEGER bi, bj
104 INTEGER i, j
105 INTEGER k, kM1, kUp, kDown
106 LOGICAL BOTTOM_LAYER
107
108 C--- The algorithm...
109 C
110 C "Correction Step"
111 C =================
112 C Here we update the horizontal velocities with the surface
113 C pressure such that the resulting flow is either consistent
114 C with the free-surface evolution or the rigid-lid:
115 C U[n] = U* + dt x d/dx P
116 C V[n] = V* + dt x d/dy P
117 C
118 C "Calculation of Gs"
119 C ===================
120 C This is where all the accelerations and tendencies (ie.
121 C physics, parameterizations etc...) are calculated
122 C w = sum_z ( div. u[n] )
123 C rho = rho ( theta[n], salt[n] )
124 C K31 = K31 ( rho )
125 C Gu[n] = Gu( u[n], v[n], w, rho, Ph, ... )
126 C Gv[n] = Gv( u[n], v[n], w, rho, Ph, ... )
127 C Gt[n] = Gt( theta[n], u[n], v[n], w, K31, ... )
128 C Gs[n] = Gs( salt[n], u[n], v[n], w, K31, ... )
129 C
130 C "Time-stepping" or "Prediction"
131 C ================================
132 C The models variables are stepped forward with the appropriate
133 C time-stepping scheme (currently we use Adams-Bashforth II)
134 C - For momentum, the result is always *only* a "prediction"
135 C in that the flow may be divergent and will be "corrected"
136 C later with a surface pressure gradient.
137 C - Normally for tracers the result is the new field at time
138 C level [n+1} *BUT* in the case of implicit diffusion the result
139 C is also *only* a prediction.
140 C - We denote "predictors" with an asterisk (*).
141 C U* = U[n] + dt x ( 3/2 Gu[n] - 1/2 Gu[n-1] )
142 C V* = V[n] + dt x ( 3/2 Gv[n] - 1/2 Gv[n-1] )
143 C theta[n+1] = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
144 C salt[n+1] = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
145 C With implicit diffusion:
146 C theta* = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
147 C salt* = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
148 C (1 + dt * K * d_zz) theta[n] = theta*
149 C (1 + dt * K * d_zz) salt[n] = salt*
150 C---
151
152 C-- Set up work arrays with valid (i.e. not NaN) values
153 C These inital values do not alter the numerical results. They
154 C just ensure that all memory references are to valid floating
155 C point numbers. This prevents spurious hardware signals due to
156 C uninitialised but inert locations.
157 DO j=1-OLy,sNy+OLy
158 DO i=1-OLx,sNx+OLx
159 xA(i,j) = 0. _d 0
160 yA(i,j) = 0. _d 0
161 uTrans(i,j) = 0. _d 0
162 vTrans(i,j) = 0. _d 0
163 aTerm(i,j) = 0. _d 0
164 xTerm(i,j) = 0. _d 0
165 cTerm(i,j) = 0. _d 0
166 mTerm(i,j) = 0. _d 0
167 pTerm(i,j) = 0. _d 0
168 fZon(i,j) = 0. _d 0
169 fMer(i,j) = 0. _d 0
170 DO K=1,nZ
171 pH (i,j,k) = 0. _d 0
172 K13(i,j,k) = 0. _d 0
173 K23(i,j,k) = 0. _d 0
174 K33(i,j,k) = 0. _d 0
175 KappaZT(i,j,k) = 0. _d 0
176 ENDDO
177 rhokm1(i,j) = 0. _d 0
178 rhok (i,j) = 0. _d 0
179 rhokp1(i,j) = 0. _d 0
180 rhotmp(i,j) = 0. _d 0
181 maskC (i,j) = 0. _d 0
182 ENDDO
183 ENDDO
184
185 DO bj=myByLo(myThid),myByHi(myThid)
186 DO bi=myBxLo(myThid),myBxHi(myThid)
187
188 C-- Set up work arrays that need valid initial values
189 DO j=1-OLy,sNy+OLy
190 DO i=1-OLx,sNx+OLx
191 wTrans(i,j) = 0. _d 0
192 wVel (i,j,1) = 0. _d 0
193 wVel (i,j,2) = 0. _d 0
194 fVerT(i,j,1) = 0. _d 0
195 fVerT(i,j,2) = 0. _d 0
196 fVerS(i,j,1) = 0. _d 0
197 fVerS(i,j,2) = 0. _d 0
198 fVerU(i,j,1) = 0. _d 0
199 fVerU(i,j,2) = 0. _d 0
200 fVerV(i,j,1) = 0. _d 0
201 fVerV(i,j,2) = 0. _d 0
202 pH(i,j,1) = 0. _d 0
203 K13(i,j,1) = 0. _d 0
204 K23(i,j,1) = 0. _d 0
205 K33(i,j,1) = 0. _d 0
206 KapGM(i,j) = GMkbackground
207 ENDDO
208 ENDDO
209
210 iMin = 1-OLx+1
211 iMax = sNx+OLx
212 jMin = 1-OLy+1
213 jMax = sNy+OLy
214
215 K = 1
216 BOTTOM_LAYER = K .EQ. Nz
217
218 C-- Calculate gradient of surface pressure
219 CALL GRAD_PSURF(
220 I bi,bj,iMin,iMax,jMin,jMax,
221 O pSurfX,pSurfY,
222 I myThid)
223
224 C-- Update fields in top level according to tendency terms
225 CALL CORRECTION_STEP(
226 I bi,bj,iMin,iMax,jMin,jMax,K,pSurfX,pSurfY,myTime,myThid)
227
228 IF ( .NOT. BOTTOM_LAYER ) THEN
229 C-- Update fields in layer below according to tendency terms
230 CALL CORRECTION_STEP(
231 I bi,bj,iMin,iMax,jMin,jMax,K+1,pSurfX,pSurfY,myTime,myThid)
232 ENDIF
233
234 C-- Density of 1st level (below W(1)) reference to level 1
235 CALL FIND_RHO(
236 I bi, bj, iMin, iMax, jMin, jMax, K, K, eosType,
237 O rhoKm1,
238 I myThid )
239
240 IF ( .NOT. BOTTOM_LAYER ) THEN
241 C-- Check static stability with layer below
242 C and mix as needed.
243 CALL FIND_RHO(
244 I bi, bj, iMin, iMax, jMin, jMax, K+1, K, eosType,
245 O rhoKp1,
246 I myThid )
247 CALL CONVECT(
248 I bi,bj,iMin,iMax,jMin,jMax,K+1,rhoKm1,rhoKp1,
249 I myTime,myIter,myThid)
250 C-- Recompute density after mixing
251 CALL FIND_RHO(
252 I bi, bj, iMin, iMax, jMin, jMax, K, K, eosType,
253 O rhoKm1,
254 I myThid )
255 ENDIF
256
257 C-- Integrate hydrostatic balance for pH with BC of pH(z=0)=0
258 CALL CALC_PH(
259 I bi,bj,iMin,iMax,jMin,jMax,K,rhoKm1,rhoKm1,
260 U pH,
261 I myThid )
262
263 DO K=2,Nz
264
265 BOTTOM_LAYER = K .EQ. Nz
266
267 IF ( .NOT. BOTTOM_LAYER ) THEN
268 C-- Update fields in layer below according to tendency terms
269 CALL CORRECTION_STEP(
270 I bi,bj,iMin,iMax,jMin,jMax,K+1,pSurfX,pSurfY,myTime,myThid)
271 ENDIF
272 C-- Update fields in layer below according to tendency terms
273 C CALL CORRECTION_STEP(
274 C I bi,bj,iMin,iMax,jMin,jMax,K,pSurfX,pSurfY,myTime,myThid)
275
276 C-- Density of K level (below W(K)) reference to K level
277 CALL FIND_RHO(
278 I bi, bj, iMin, iMax, jMin, jMax, K, K, eosType,
279 O rhoK,
280 I myThid )
281 IF ( .NOT. BOTTOM_LAYER ) THEN
282 C-- Check static stability with layer below
283 C and mix as needed.
284 C-- Density of K+1 level (below W(K+1)) reference to K level
285 CALL FIND_RHO(
286 I bi, bj, iMin, iMax, jMin, jMax, K+1, K, eosType,
287 O rhoKp1,
288 I myThid )
289 CALL CONVECT(
290 I bi,bj,iMin,iMax,jMin,jMax,K+1,rhoK,rhoKp1,
291 I myTime,myIter,myThid)
292 C-- Recompute density after mixing
293 CALL FIND_RHO(
294 I bi, bj, iMin, iMax, jMin, jMax, K, K, eosType,
295 O rhoK,
296 I myThid )
297 ENDIF
298 C-- Integrate hydrostatic balance for pH with BC of pH(z=0)=0
299 CALL CALC_PH(
300 I bi,bj,iMin,iMax,jMin,jMax,K,rhoKm1,rhoK,
301 U pH,
302 I myThid )
303 C-- Calculate iso-neutral slopes for the GM/Redi parameterisation
304 CALL FIND_RHO(
305 I bi, bj, iMin, iMax, jMin, jMax, K-1, K, eosType,
306 O rhoTmp,
307 I myThid )
308 CALL CALC_ISOSLOPES(
309 I bi, bj, iMin, iMax, jMin, jMax, K,
310 I rhoKm1, rhoK, rhotmp,
311 O K13, K23, K33, KapGM,
312 I myThid )
313 DO J=jMin,jMax
314 DO I=iMin,iMax
315 rhoKm1(I,J)=rhoK(I,J)
316 ENDDO
317 ENDDO
318
319 ENDDO ! K
320
321 DO K = Nz, 1, -1
322 kM1 =max(1,k-1) ! Points to level above k (=k-1)
323 kUp =1+MOD(k+1,2) ! Cycles through 1,2 to point to layer above
324 kDown=1+MOD(k,2) ! Cycles through 2,1 to point to current layer
325 iMin = 1-OLx+2
326 iMax = sNx+OLx-1
327 jMin = 1-OLy+2
328 jMax = sNy+OLy-1
329
330 C-- Get temporary terms used by tendency routines
331 CALL CALC_COMMON_FACTORS (
332 I bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,
333 O xA,yA,uTrans,vTrans,wTrans,wVel,maskC,maskUp,
334 I myThid)
335
336 C-- Calculate the total vertical diffusivity
337 CALL CALC_DIFFUSIVITY(
338 I bi,bj,iMin,iMax,jMin,jMax,K,
339 I maskC,maskUp,KapGM,K33,
340 O KappaZT,KappaZS,
341 I myThid)
342
343 C-- Calculate accelerations in the momentum equations
344 IF ( momStepping ) THEN
345 CALL CALC_MOM_RHS(
346 I bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,
347 I xA,yA,uTrans,vTrans,wTrans,wVel,maskC,
348 I pH,
349 U aTerm,xTerm,cTerm,mTerm,pTerm,
350 U fZon, fMer, fVerU, fVerV,
351 I myThid)
352 ENDIF
353
354 C-- Calculate active tracer tendencies
355 IF ( tempStepping ) THEN
356 CALL CALC_GT(
357 I bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown,
358 I xA,yA,uTrans,vTrans,wTrans,maskUp,maskC,
359 I K13,K23,KappaZT,KapGM,
360 U aTerm,xTerm,fZon,fMer,fVerT,
361 I myThid)
362 ENDIF
363 IF ( saltStepping ) THEN
364 CALL CALC_GS(
365 I bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown,
366 I xA,yA,uTrans,vTrans,wTrans,maskUp,maskC,
367 I K13,K23,KappaZS,KapGM,
368 U aTerm,xTerm,fZon,fMer,fVerS,
369 I myThid)
370 ENDIF
371
372 C-- Prediction step (step forward all model variables)
373 CALL TIMESTEP(
374 I bi,bj,iMin,iMax,jMin,jMax,K,
375 I myThid)
376
377 C-- Diagnose barotropic divergence of predicted fields
378 CALL DIV_G(
379 I bi,bj,iMin,iMax,jMin,jMax,K,
380 I xA,yA,
381 I myThid)
382
383 C-- Cumulative diagnostic calculations (ie. time-averaging)
384 #ifdef ALLOW_DIAGNOSTICS
385 IF (taveFreq.GT.0.) THEN
386 CALL DO_TIME_AVERAGES(
387 I myTime, myIter, bi, bj, K, kUp, kDown,
388 I K13, K23, wVel,
389 I myThid )
390 ENDIF
391 #endif
392
393 ENDDO ! K
394
395 C-- Implicit diffusion
396 IF (implicitDiffusion) THEN
397 CALL IMPLDIFF( bi, bj, iMin, iMax, jMin, jMax,
398 I KappaZT,KappaZS,
399 I myThid )
400 ENDIF
401
402 ENDDO
403 ENDDO
404
405 C write(0,*) 'dynamics: pS ',minval(cg2d_x(1:sNx,1:sNy,:,:)),
406 C & maxval(cg2d_x(1:sNx,1:sNy,:,:))
407 C write(0,*) 'dynamics: U ',minval(uVel(1:sNx,1:sNy,1,:,:),mask=uVel(1:sNx,1:sNy,1,:,:).NE.0.),
408 C & maxval(uVel(1:sNx,1:sNy,1,:,:),mask=uVel(1:sNx,1:sNy,1,:,:).NE.0.)
409 C write(0,*) 'dynamics: V ',minval(vVel(1:sNx,1:sNy,1,:,:),mask=vVel(1:sNx,1:sNy,1,:,:).NE.0.),
410 C & maxval(vVel(1:sNx,1:sNy,1,:,:),mask=vVel(1:sNx,1:sNy,1,:,:).NE.0.)
411 C write(0,*) 'dynamics: wVel(1) ',
412 C & minval(wVel(1:sNx,1:sNy,1),mask=wVel(1:sNx,1:sNy,1).NE.0.),
413 C & maxval(wVel(1:sNx,1:sNy,1),mask=wVel(1:sNx,1:sNy,1).NE.0.)
414 C write(0,*) 'dynamics: wVel(2) ',
415 C & minval(wVel(1:sNx,1:sNy,2),mask=wVel(1:sNx,1:sNy,2).NE.0.),
416 C & maxval(wVel(1:sNx,1:sNy,2),mask=wVel(1:sNx,1:sNy,2).NE.0.)
417 cblk write(0,*) 'dynamics: K13',minval(K13(1:sNx,1:sNy,:)),
418 cblk & maxval(K13(1:sNx,1:sNy,:))
419 cblk write(0,*) 'dynamics: K23',minval(K23(1:sNx,1:sNy,:)),
420 cblk & maxval(K23(1:sNx,1:sNy,:))
421 cblk write(0,*) 'dynamics: K33',minval(K33(1:sNx,1:sNy,:)),
422 cblk & maxval(K33(1:sNx,1:sNy,:))
423 C write(0,*) 'dynamics: gT ',minval(gT(1:sNx,1:sNy,:,:,:)),
424 C & maxval(gT(1:sNx,1:sNy,:,:,:))
425 C write(0,*) 'dynamics: T ',minval(Theta(1:sNx,1:sNy,:,:,:)),
426 C & maxval(Theta(1:sNx,1:sNy,:,:,:))
427 C write(0,*) 'dynamics: gS ',minval(gS(1:sNx,1:sNy,:,:,:)),
428 C & maxval(gS(1:sNx,1:sNy,:,:,:))
429 C write(0,*) 'dynamics: S ',minval(salt(1:sNx,1:sNy,:,:,:)),
430 C & maxval(salt(1:sNx,1:sNy,:,:,:))
431 C write(0,*) 'dynamics: pH ',minval(pH/(Gravity*Rhonil),mask=ph.NE.0.),
432 C & maxval(pH/(Gravity*Rhonil))
433
434 RETURN
435 END

  ViewVC Help
Powered by ViewVC 1.1.22