/[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.20 - (hide annotations) (download)
Mon Jun 15 05:17:42 1998 UTC (25 years, 11 months ago) by cnh
Branch: MAIN
Changes since 1.19: +13 -13 lines
Commented out debugging at end

1 cnh 1.20 C $Header: /u/gcmpack/models/MITgcmUV/model/src/dynamics.F,v 1.19 1998/06/15 05:13:56 cnh 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 cnh 1.14 C wVel o uTrans: Zonal transport
43 cnh 1.1 C o vTrans: Meridional transport
44     C o wTrans: Vertical transport
45 cnh 1.14 C o wVel: Vertical velocity at upper and lower
46     C cell faces.
47 cnh 1.1 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 cnh 1.14 _RL wVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
74 cnh 1.1 _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 adcroft 1.3 _RL rhokm1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
89     _RL rhokp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
90 cnh 1.19 _RL rhok (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
91 adcroft 1.10 _RL rhotmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
92 adcroft 1.4 _RL pSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
93     _RL pSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
94 adcroft 1.6 _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 adcroft 1.12 _RL KappaZT(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nz)
99 adcroft 1.18 _RL KappaZS(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nz)
100 adcroft 1.12
101 cnh 1.1 INTEGER iMin, iMax
102     INTEGER jMin, jMax
103     INTEGER bi, bj
104     INTEGER i, j
105     INTEGER k, kM1, kUp, kDown
106 cnh 1.19 LOGICAL BOTTOM_LAYER
107 cnh 1.1
108 adcroft 1.11 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 adcroft 1.12 C "Time-stepping" or "Prediction"
131 adcroft 1.11 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 adcroft 1.12 C With implicit diffusion:
146 adcroft 1.11 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 adcroft 1.12 C (1 + dt * K * d_zz) theta[n] = theta*
149     C (1 + dt * K * d_zz) salt[n] = salt*
150 adcroft 1.11 C---
151    
152 cnh 1.1 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 adcroft 1.5 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 cnh 1.1 DO K=1,nZ
171 adcroft 1.5 pH (i,j,k) = 0. _d 0
172 adcroft 1.6 K13(i,j,k) = 0. _d 0
173     K23(i,j,k) = 0. _d 0
174     K33(i,j,k) = 0. _d 0
175 adcroft 1.12 KappaZT(i,j,k) = 0. _d 0
176 cnh 1.1 ENDDO
177 adcroft 1.5 rhokm1(i,j) = 0. _d 0
178 cnh 1.19 rhok (i,j) = 0. _d 0
179 adcroft 1.5 rhokp1(i,j) = 0. _d 0
180 adcroft 1.10 rhotmp(i,j) = 0. _d 0
181 cnh 1.16 maskC (i,j) = 0. _d 0
182 cnh 1.1 ENDDO
183     ENDDO
184    
185     DO bj=myByLo(myThid),myByHi(myThid)
186     DO bi=myBxLo(myThid),myBxHi(myThid)
187    
188 cnh 1.7 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 cnh 1.14 wVel (i,j,1) = 0. _d 0
193     wVel (i,j,2) = 0. _d 0
194 cnh 1.7 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 adcroft 1.11 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) = 0. _d 0
207 cnh 1.7 ENDDO
208     ENDDO
209    
210 cnh 1.1 iMin = 1-OLx+1
211     iMax = sNx+OLx
212     jMin = 1-OLy+1
213     jMax = sNy+OLy
214    
215 cnh 1.19 K = 1
216     BOTTOM_LAYER = K .EQ. Nz
217    
218 adcroft 1.4 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 adcroft 1.11 CALL CORRECTION_STEP(
226 cnh 1.19 I bi,bj,iMin,iMax,jMin,jMax,K,pSurfX,pSurfY,myThid)
227 cnh 1.1
228 cnh 1.7 C-- Density of 1st level (below W(1)) reference to level 1
229     CALL FIND_RHO(
230 cnh 1.19 I bi, bj, iMin, iMax, jMin, jMax, K, K, eosType,
231 cnh 1.7 O rhoKm1,
232     I myThid )
233 cnh 1.19
234     IF ( .NOT. BOTTOM_LAYER ) THEN
235     C-- Check static stability with layer below
236     C and mix as needed.
237     CALL FIND_RHO(
238     I bi, bj, iMin, iMax, jMin, jMax, K+1, K, eosType,
239     O rhoKp1,
240     I myThid )
241     CALL CONVECT(
242     I bi,bj,iMin,iMax,jMin,jMax,K+1,rhoKm1,rhoKp1,
243     I myTime,myIter,myThid)
244     C-- Recompute density after mixing
245     CALL FIND_RHO(
246     I bi, bj, iMin, iMax, jMin, jMax, K, K, eosType,
247     O rhoKm1,
248     I myThid )
249     ENDIF
250    
251 cnh 1.7 C-- Integrate hydrostatic balance for pH with BC of pH(z=0)=0
252     CALL CALC_PH(
253 cnh 1.19 I bi,bj,iMin,iMax,jMin,jMax,K,rhoKm1,rhoKm1,
254 cnh 1.7 U pH,
255 adcroft 1.5 I myThid )
256    
257 adcroft 1.3 DO K=2,Nz
258 cnh 1.19
259     BOTTOM_LAYER = K .EQ. Nz
260    
261     C-- Update fields in Kth level according to tendency terms
262     CALL CORRECTION_STEP(
263     I bi,bj,iMin,iMax,jMin,jMax,K,pSurfX,pSurfY,myThid)
264     C-- Density of K level (below W(K)) reference to K level
265     CALL FIND_RHO(
266     I bi, bj, iMin, iMax, jMin, jMax, K, K, eosType,
267     O rhoK,
268     I myThid )
269     IF ( .NOT. BOTTOM_LAYER ) THEN
270     C-- Check static stability with layer below
271     C and mix as needed.
272     C-- Density of K+1 level (below W(K+1)) reference to K level
273     CALL FIND_RHO(
274     I bi, bj, iMin, iMax, jMin, jMax, K+1, K, eosType,
275     O rhoKp1,
276     I myThid )
277     CALL CONVECT(
278     I bi,bj,iMin,iMax,jMin,jMax,K+1,rhoK,rhoKp1,
279     I myTime,myIter,myThid)
280     C-- Recompute density after mixing
281     CALL FIND_RHO(
282     I bi, bj, iMin, iMax, jMin, jMax, K, K, eosType,
283     O rhoK,
284     I myThid )
285     ENDIF
286     C-- Integrate hydrostatic balance for pH with BC of pH(z=0)=0
287     CALL CALC_PH(
288     I bi,bj,iMin,iMax,jMin,jMax,K,rhoKm1,rhoK,
289     U pH,
290     I myThid )
291     C-- Calculate iso-neutral slopes for the GM/Redi parameterisation
292     CALL FIND_RHO(
293     I bi, bj, iMin, iMax, jMin, jMax, K-1, K, eosType,
294     O rhoTmp,
295     I myThid )
296     CALL CALC_ISOSLOPES(
297     I bi, bj, iMin, iMax, jMin, jMax, K,
298     I rhoKm1, rhoK, rhotmp,
299     O K13, K23, K33, KapGM,
300     I myThid )
301     DO J=jMin,jMax
302     DO I=iMin,iMax
303     rhoKm1(I,J)=rhoK(I,J)
304     ENDDO
305 adcroft 1.10 ENDDO
306 cnh 1.1
307 adcroft 1.11 ENDDO ! K
308    
309 cnh 1.1 DO K = Nz, 1, -1
310     kM1 =max(1,k-1) ! Points to level above k (=k-1)
311     kUp =1+MOD(k+1,2) ! Cycles through 1,2 to point to layer above
312     kDown=1+MOD(k,2) ! Cycles through 2,1 to point to current layer
313     iMin = 1-OLx+2
314     iMax = sNx+OLx-1
315     jMin = 1-OLy+2
316     jMax = sNy+OLy-1
317    
318     C-- Get temporary terms used by tendency routines
319     CALL CALC_COMMON_FACTORS (
320     I bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,
321 cnh 1.14 O xA,yA,uTrans,vTrans,wTrans,wVel,maskC,maskUp,
322 cnh 1.1 I myThid)
323    
324 adcroft 1.12 C-- Calculate the total vertical diffusivity
325     CALL CALC_DIFFUSIVITY(
326     I bi,bj,iMin,iMax,jMin,jMax,K,
327     I maskC,maskUp,KapGM,K33,
328 adcroft 1.18 O KappaZT,KappaZS,
329 adcroft 1.12 I myThid)
330    
331 cnh 1.1 C-- Calculate accelerations in the momentum equations
332 cnh 1.9 IF ( momStepping ) THEN
333     CALL CALC_MOM_RHS(
334     I bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,
335 cnh 1.14 I xA,yA,uTrans,vTrans,wTrans,wVel,maskC,
336 cnh 1.9 I pH,
337     U aTerm,xTerm,cTerm,mTerm,pTerm,
338     U fZon, fMer, fVerU, fVerV,
339     I myThid)
340     ENDIF
341 cnh 1.1
342     C-- Calculate active tracer tendencies
343 cnh 1.9 IF ( tempStepping ) THEN
344     CALL CALC_GT(
345     I bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown,
346     I xA,yA,uTrans,vTrans,wTrans,maskUp,
347 adcroft 1.12 I K13,K23,KappaZT,KapGM,
348 cnh 1.9 U aTerm,xTerm,fZon,fMer,fVerT,
349     I myThid)
350     ENDIF
351 adcroft 1.18 IF ( saltStepping ) THEN
352     CALL CALC_GS(
353     I bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown,
354     I xA,yA,uTrans,vTrans,wTrans,maskUp,
355     I K13,K23,KappaZS,KapGM,
356     U aTerm,xTerm,fZon,fMer,fVerS,
357     I myThid)
358     ENDIF
359 cnh 1.1
360 adcroft 1.11 C-- Prediction step (step forward all model variables)
361     CALL TIMESTEP(
362     I bi,bj,iMin,iMax,jMin,jMax,K,
363     I myThid)
364    
365     C-- Diagnose barotropic divergence of predicted fields
366     CALL DIV_G(
367     I bi,bj,iMin,iMax,jMin,jMax,K,
368     I xA,yA,
369     I myThid)
370    
371     ENDDO ! K
372 adcroft 1.12
373     C-- Implicit diffusion
374     IF (implicitDiffusion) THEN
375     CALL IMPLDIFF( bi, bj, iMin, iMax, jMin, jMax,
376 adcroft 1.18 I KappaZT,KappaZS,
377 adcroft 1.12 I myThid )
378     ENDIF
379 cnh 1.1
380     ENDDO
381     ENDDO
382 adcroft 1.6
383 cnh 1.19 C write(0,*) 'dynamics: pS ',minval(cg2d_x(1:sNx,1:sNy,:,:)),
384     C & maxval(cg2d_x(1:sNx,1:sNy,:,:))
385 cnh 1.20 C write(0,*) 'dynamics: U ',minval(uVel(1:sNx,1:sNy,1,:,:),mask=uVel(1:sNx,1:sNy,1,:,:).NE.0.),
386     C & maxval(uVel(1:sNx,1:sNy,1,:,:))
387     C write(0,*) 'dynamics: V ',minval(vVel(1:sNx,1:sNy,1,:,:),mask=vVel(1:sNx,1:sNy,1,:,:).NE.0.),
388     C & maxval(vVel(1:sNx,1:sNy,1,:,:))
389     C write(0,*) 'dynamics: wVel(1) ',
390     C & minval(wVel(1:sNx,1:sNy,1),mask=wVel(1:sNx,1:sNy,1).NE.0.),
391     C & maxval(wVel(1:sNx,1:sNy,1))
392     C write(0,*) 'dynamics: wVel(2) ',
393     C & minval(wVel(1:sNx,1:sNy,2),mask=wVel(1:sNx,1:sNy,2).NE.0.),
394     C & maxval(wVel(1:sNx,1:sNy,2))
395 adcroft 1.15 cblk write(0,*) 'dynamics: K13',minval(K13(1:sNx,1:sNy,:)),
396     cblk & maxval(K13(1:sNx,1:sNy,:))
397     cblk write(0,*) 'dynamics: K23',minval(K23(1:sNx,1:sNy,:)),
398     cblk & maxval(K23(1:sNx,1:sNy,:))
399     cblk write(0,*) 'dynamics: K33',minval(K33(1:sNx,1:sNy,:)),
400     cblk & maxval(K33(1:sNx,1:sNy,:))
401 cnh 1.19 C write(0,*) 'dynamics: gT ',minval(gT(1:sNx,1:sNy,:,:,:)),
402     C & maxval(gT(1:sNx,1:sNy,:,:,:))
403     C write(0,*) 'dynamics: T ',minval(Theta(1:sNx,1:sNy,:,:,:)),
404     C & maxval(Theta(1:sNx,1:sNy,:,:,:))
405     C write(0,*) 'dynamics: gS ',minval(gS(1:sNx,1:sNy,:,:,:)),
406     C & maxval(gS(1:sNx,1:sNy,:,:,:))
407     C write(0,*) 'dynamics: S ',minval(salt(1:sNx,1:sNy,:,:,:)),
408     C & maxval(salt(1:sNx,1:sNy,:,:,:))
409 cnh 1.20 C write(0,*) 'dynamics: pH ',minval(pH/(Gravity*Rhonil),mask=ph.NE.0.),
410     C & maxval(pH/(Gravity*Rhonil))
411 cnh 1.1
412     RETURN
413     END

  ViewVC Help
Powered by ViewVC 1.1.22