/[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.12 - (hide annotations) (download)
Mon Jun 1 22:27:14 1998 UTC (25 years, 11 months ago) by adcroft
Branch: MAIN
Changes since 1.11: +41 -26 lines
Implemented implicit vertical diffusion (tracers only).
Involved introducing a "total" diffusivity array (local 3D)
calculated by calc_diffusivity().
Made some small changes to time-stepping algorithm.
Switched on by setting implicitZdiffusion.
(note: *Not* fully tested with topography. But when switched off
this does produce identical results)

1 adcroft 1.12 C $Header: /u/gcmpack/models/MITgcmUV/model/src/dynamics.F,v 1.11 1998/06/01 20:36:13 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 adcroft 1.12 _RL KappaZT(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nz)
95    
96 cnh 1.1 INTEGER iMin, iMax
97     INTEGER jMin, jMax
98     INTEGER bi, bj
99     INTEGER i, j
100     INTEGER k, kM1, kUp, kDown
101    
102 adcroft 1.11 C--- The algorithm...
103     C
104     C "Correction Step"
105     C =================
106     C Here we update the horizontal velocities with the surface
107     C pressure such that the resulting flow is either consistent
108     C with the free-surface evolution or the rigid-lid:
109     C U[n] = U* + dt x d/dx P
110     C V[n] = V* + dt x d/dy P
111     C
112     C "Calculation of Gs"
113     C ===================
114     C This is where all the accelerations and tendencies (ie.
115     C physics, parameterizations etc...) are calculated
116     C w = sum_z ( div. u[n] )
117     C rho = rho ( theta[n], salt[n] )
118     C K31 = K31 ( rho )
119     C Gu[n] = Gu( u[n], v[n], w, rho, Ph, ... )
120     C Gv[n] = Gv( u[n], v[n], w, rho, Ph, ... )
121     C Gt[n] = Gt( theta[n], u[n], v[n], w, K31, ... )
122     C Gs[n] = Gs( salt[n], u[n], v[n], w, K31, ... )
123     C
124 adcroft 1.12 C "Time-stepping" or "Prediction"
125 adcroft 1.11 C ================================
126     C The models variables are stepped forward with the appropriate
127     C time-stepping scheme (currently we use Adams-Bashforth II)
128     C - For momentum, the result is always *only* a "prediction"
129     C in that the flow may be divergent and will be "corrected"
130     C later with a surface pressure gradient.
131     C - Normally for tracers the result is the new field at time
132     C level [n+1} *BUT* in the case of implicit diffusion the result
133     C is also *only* a prediction.
134     C - We denote "predictors" with an asterisk (*).
135     C U* = U[n] + dt x ( 3/2 Gu[n] - 1/2 Gu[n-1] )
136     C V* = V[n] + dt x ( 3/2 Gv[n] - 1/2 Gv[n-1] )
137     C theta[n+1] = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
138     C salt[n+1] = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
139 adcroft 1.12 C With implicit diffusion:
140 adcroft 1.11 C theta* = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
141     C salt* = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
142 adcroft 1.12 C (1 + dt * K * d_zz) theta[n] = theta*
143     C (1 + dt * K * d_zz) salt[n] = salt*
144 adcroft 1.11 C---
145    
146 cnh 1.1 C-- Set up work arrays with valid (i.e. not NaN) values
147     C These inital values do not alter the numerical results. They
148     C just ensure that all memory references are to valid floating
149     C point numbers. This prevents spurious hardware signals due to
150     C uninitialised but inert locations.
151     DO j=1-OLy,sNy+OLy
152     DO i=1-OLx,sNx+OLx
153 adcroft 1.5 xA(i,j) = 0. _d 0
154     yA(i,j) = 0. _d 0
155     uTrans(i,j) = 0. _d 0
156     vTrans(i,j) = 0. _d 0
157     aTerm(i,j) = 0. _d 0
158     xTerm(i,j) = 0. _d 0
159     cTerm(i,j) = 0. _d 0
160     mTerm(i,j) = 0. _d 0
161     pTerm(i,j) = 0. _d 0
162     fZon(i,j) = 0. _d 0
163     fMer(i,j) = 0. _d 0
164 cnh 1.1 DO K=1,nZ
165 adcroft 1.5 pH (i,j,k) = 0. _d 0
166 adcroft 1.6 K13(i,j,k) = 0. _d 0
167     K23(i,j,k) = 0. _d 0
168     K33(i,j,k) = 0. _d 0
169 adcroft 1.12 KappaZT(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 adcroft 1.12 C-- Calculate the total vertical diffusivity
307     CALL CALC_DIFFUSIVITY(
308     I bi,bj,iMin,iMax,jMin,jMax,K,
309     I maskC,maskUp,KapGM,K33,
310     O KappaZT,
311     I myThid)
312    
313    
314 cnh 1.1 C-- Calculate accelerations in the momentum equations
315 cnh 1.9 IF ( momStepping ) THEN
316     CALL CALC_MOM_RHS(
317     I bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,
318     I xA,yA,uTrans,vTrans,wTrans,maskC,
319     I pH,
320     U aTerm,xTerm,cTerm,mTerm,pTerm,
321     U fZon, fMer, fVerU, fVerV,
322     I myThid)
323     ENDIF
324 cnh 1.1
325     C-- Calculate active tracer tendencies
326 cnh 1.9 IF ( tempStepping ) THEN
327     CALL CALC_GT(
328     I bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown,
329     I xA,yA,uTrans,vTrans,wTrans,maskUp,
330 adcroft 1.12 I K13,K23,KappaZT,KapGM,
331 cnh 1.9 U aTerm,xTerm,fZon,fMer,fVerT,
332     I myThid)
333     ENDIF
334 cnh 1.1 Cdbg CALL CALC_GS(
335     Cdbg I bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown,
336     Cdbg I xA,yA,uTrans,vTrans,wTrans,maskUp,
337 adcroft 1.6 Cdbg I K13,K23,K33,KapGM,
338 cnh 1.1 Cdbg U aTerm,xTerm,fZon,fMer,fVerS,
339     Cdbg I myThid)
340    
341 adcroft 1.11 C-- Prediction step (step forward all model variables)
342     CALL TIMESTEP(
343     I bi,bj,iMin,iMax,jMin,jMax,K,
344     I myThid)
345    
346     C-- Diagnose barotropic divergence of predicted fields
347     CALL DIV_G(
348     I bi,bj,iMin,iMax,jMin,jMax,K,
349     I xA,yA,
350     I myThid)
351    
352     ENDDO ! K
353 adcroft 1.12
354     C-- Implicit diffusion
355     IF (implicitDiffusion) THEN
356     CALL IMPLDIFF( bi, bj, iMin, iMax, jMin, jMax,
357     I KappaZT,
358     I myThid )
359     ENDIF
360 cnh 1.1
361     ENDDO
362     ENDDO
363 adcroft 1.6
364 adcroft 1.12 write(0,*) 'dynamics: pS',minval(cg2d_x),maxval(cg2d_x)
365     write(0,*) 'dynamics: U',minval(uVel(1:sNx,1:sNy,:,:,:)),
366     & maxval(uVel(1:sNx,1:sNy,:,:,:))
367     write(0,*) 'dynamics: V',minval(vVel(1:sNx,1:sNy,:,:,:)),
368     & maxval(vVel(1:sNx,1:sNy,:,:,:))
369     write(0,*) 'dynamics: K13',minval(K13(1:sNx,1:sNy,:)),
370     & maxval(K13(1:sNx,1:sNy,:))
371     write(0,*) 'dynamics: K23',minval(K23(1:sNx,1:sNy,:)),
372     & maxval(K23(1:sNx,1:sNy,:))
373     write(0,*) 'dynamics: K33',minval(K33(1:sNx,1:sNy,:)),
374     & maxval(K33(1:sNx,1:sNy,:))
375     write(0,*) 'dynamics: gT',minval(gT(1:sNx,1:sNy,:,:,:)),
376     & maxval(gT(1:sNx,1:sNy,:,:,:))
377     write(0,*) 'dynamics: T',minval(Theta(1:sNx,1:sNy,:,:,:)),
378     & maxval(Theta(1:sNx,1:sNy,:,:,:))
379     write(0,*) 'dynamics: pH',minval(pH/(Gravity*Rhonil)),
380     & maxval(pH/(Gravity*Rhonil))
381 cnh 1.1
382     RETURN
383     END

  ViewVC Help
Powered by ViewVC 1.1.22