/[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.44 - (hide annotations) (download)
Wed Jul 28 16:32:12 1999 UTC (24 years, 10 months ago) by adcroft
Branch: MAIN
CVS Tags: checkpoint24
Changes since 1.43: +5 -2 lines
Added a parameter "implicitViscosity" to separately control implicit
viscosity and diffusion.

1 adcroft 1.44 C $Header: /u/gcmpack/models/MITgcmUV/model/src/dynamics.F,v 1.43 1999/05/24 15:42:23 adcroft Exp $
2 cnh 1.1
3 adcroft 1.24 #include "CPP_OPTIONS.h"
4 cnh 1.1
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 adcroft 1.40 IMPLICIT NONE
24 cnh 1.1
25     C == Global variables ===
26     #include "SIZE.h"
27     #include "EEPARAMS.h"
28     #include "CG2D.h"
29 adcroft 1.6 #include "PARAMS.h"
30 adcroft 1.3 #include "DYNVARS.h"
31 adcroft 1.42 #include "GRID.h"
32     #ifdef ALLOW_KPP
33     #include "KPPMIX.h"
34     #endif
35 cnh 1.1
36     C == Routine arguments ==
37 cnh 1.8 C myTime - Current time in simulation
38     C myIter - Current iteration number in simulation
39 cnh 1.1 C myThid - Thread number for this instance of the routine.
40     INTEGER myThid
41 cnh 1.8 _RL myTime
42     INTEGER myIter
43 cnh 1.1
44     C == Local variables
45     C xA, yA - Per block temporaries holding face areas
46 cnh 1.38 C uTrans, vTrans, rTrans - Per block temporaries holding flow
47     C transport
48 cnh 1.30 C rVel o uTrans: Zonal transport
49 cnh 1.1 C o vTrans: Meridional transport
50 cnh 1.30 C o rTrans: Vertical transport
51 cnh 1.38 C o rVel: Vertical velocity at upper and
52     C lower cell faces.
53 cnh 1.1 C maskC,maskUp o maskC: land/water mask for tracer cells
54     C o maskUp: land/water mask for W points
55     C aTerm, xTerm, cTerm - Work arrays for holding separate terms in
56     C mTerm, pTerm, tendency equations.
57     C fZon, fMer, fVer[STUV] o aTerm: Advection term
58     C o xTerm: Mixing term
59     C o cTerm: Coriolis term
60     C o mTerm: Metric term
61     C o pTerm: Pressure term
62     C o fZon: Zonal flux term
63     C o fMer: Meridional flux term
64     C o fVer: Vertical flux term - note fVer
65     C is "pipelined" in the vertical
66     C so we need an fVer for each
67     C variable.
68 cnh 1.38 C rhoK, rhoKM1 - Density at current level, level above and level
69     C below.
70 cnh 1.26 C rhoKP1
71     C buoyK, buoyKM1 - Buoyancy at current level and level above.
72 cnh 1.31 C phiHyd - Hydrostatic part of the potential phiHydi.
73 cnh 1.38 C In z coords phiHydiHyd is the hydrostatic
74     C pressure anomaly
75     C In p coords phiHydiHyd is the geopotential
76     C surface height
77 cnh 1.30 C anomaly.
78     C etaSurfX, - Holds surface elevation gradient in X and Y.
79     C etaSurfY
80     C K13, K23, K33 - Non-zero elements of small-angle approximation
81     C diffusion tensor.
82     C KapGM - Spatially varying Visbeck et. al mixing coeff.
83     C KappaRT, - Total diffusion in vertical for T and S.
84 cnh 1.38 C KappaRS (background + spatially varying, isopycnal term).
85 cnh 1.30 C iMin, iMax - Ranges and sub-block indices on which calculations
86     C jMin, jMax are applied.
87 cnh 1.1 C bi, bj
88 cnh 1.30 C k, kUp, - Index for layer above and below. kUp and kDown
89 cnh 1.38 C kDown, kM1 are switched with layer to be the appropriate
90     C index into fVerTerm.
91 cnh 1.30 _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
92     _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
93     _RL uTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
94     _RL vTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
95     _RL rTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
96     _RL rVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
97     _RS maskC (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98     _RS maskUp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
99     _RL aTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
100     _RL xTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
101     _RL cTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
102     _RL mTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
103     _RL pTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
104     _RL fZon (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105     _RL fMer (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
106     _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
107     _RL fVerS (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
108     _RL fVerU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
109     _RL fVerV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
110 cnh 1.31 _RL phiHyd (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
111 cnh 1.30 _RL rhokm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
112     _RL rhokp1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
113     _RL rhok (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
114     _RL buoyKM1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
115     _RL buoyK (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
116     _RL rhotmp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
117 cnh 1.29 _RL etaSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
118     _RL etaSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
119 cnh 1.31 _RL K13 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
120     _RL K23 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
121     _RL K33 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
122 cnh 1.30 _RL KapGM (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
123 cnh 1.31 _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
124     _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
125 adcroft 1.42 _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
126     _RL KappaRV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
127 adcroft 1.12
128 cnh 1.1 INTEGER iMin, iMax
129     INTEGER jMin, jMax
130     INTEGER bi, bj
131     INTEGER i, j
132     INTEGER k, kM1, kUp, kDown
133 cnh 1.19 LOGICAL BOTTOM_LAYER
134 cnh 1.1
135 adcroft 1.11 C--- The algorithm...
136     C
137     C "Correction Step"
138     C =================
139     C Here we update the horizontal velocities with the surface
140     C pressure such that the resulting flow is either consistent
141     C with the free-surface evolution or the rigid-lid:
142     C U[n] = U* + dt x d/dx P
143     C V[n] = V* + dt x d/dy P
144     C
145     C "Calculation of Gs"
146     C ===================
147     C This is where all the accelerations and tendencies (ie.
148 cnh 1.31 C phiHydysics, parameterizations etc...) are calculated
149 cnh 1.27 C rVel = sum_r ( div. u[n] )
150 adcroft 1.11 C rho = rho ( theta[n], salt[n] )
151 cnh 1.27 C b = b(rho, theta)
152 adcroft 1.11 C K31 = K31 ( rho )
153 cnh 1.27 C Gu[n] = Gu( u[n], v[n], rVel, b, ... )
154     C Gv[n] = Gv( u[n], v[n], rVel, b, ... )
155     C Gt[n] = Gt( theta[n], u[n], v[n], rVel, K31, ... )
156     C Gs[n] = Gs( salt[n], u[n], v[n], rVel, K31, ... )
157 adcroft 1.11 C
158 adcroft 1.12 C "Time-stepping" or "Prediction"
159 adcroft 1.11 C ================================
160     C The models variables are stepped forward with the appropriate
161     C time-stepping scheme (currently we use Adams-Bashforth II)
162     C - For momentum, the result is always *only* a "prediction"
163     C in that the flow may be divergent and will be "corrected"
164     C later with a surface pressure gradient.
165     C - Normally for tracers the result is the new field at time
166     C level [n+1} *BUT* in the case of implicit diffusion the result
167     C is also *only* a prediction.
168     C - We denote "predictors" with an asterisk (*).
169     C U* = U[n] + dt x ( 3/2 Gu[n] - 1/2 Gu[n-1] )
170     C V* = V[n] + dt x ( 3/2 Gv[n] - 1/2 Gv[n-1] )
171     C theta[n+1] = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
172     C salt[n+1] = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
173 adcroft 1.12 C With implicit diffusion:
174 adcroft 1.11 C theta* = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
175     C salt* = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
176 adcroft 1.12 C (1 + dt * K * d_zz) theta[n] = theta*
177     C (1 + dt * K * d_zz) salt[n] = salt*
178 adcroft 1.11 C---
179    
180 cnh 1.1 C-- Set up work arrays with valid (i.e. not NaN) values
181     C These inital values do not alter the numerical results. They
182     C just ensure that all memory references are to valid floating
183     C point numbers. This prevents spurious hardware signals due to
184     C uninitialised but inert locations.
185     DO j=1-OLy,sNy+OLy
186     DO i=1-OLx,sNx+OLx
187 adcroft 1.5 xA(i,j) = 0. _d 0
188     yA(i,j) = 0. _d 0
189     uTrans(i,j) = 0. _d 0
190     vTrans(i,j) = 0. _d 0
191     aTerm(i,j) = 0. _d 0
192     xTerm(i,j) = 0. _d 0
193     cTerm(i,j) = 0. _d 0
194     mTerm(i,j) = 0. _d 0
195     pTerm(i,j) = 0. _d 0
196     fZon(i,j) = 0. _d 0
197     fMer(i,j) = 0. _d 0
198 cnh 1.31 DO K=1,Nr
199     phiHyd (i,j,k) = 0. _d 0
200 cnh 1.30 K13(i,j,k) = 0. _d 0
201     K23(i,j,k) = 0. _d 0
202     K33(i,j,k) = 0. _d 0
203 cnh 1.31 KappaRT(i,j,k) = 0. _d 0
204     KappaRS(i,j,k) = 0. _d 0
205 cnh 1.1 ENDDO
206 cnh 1.30 rhoKM1 (i,j) = 0. _d 0
207     rhok (i,j) = 0. _d 0
208     rhoKP1 (i,j) = 0. _d 0
209     rhoTMP (i,j) = 0. _d 0
210 cnh 1.26 buoyKM1(i,j) = 0. _d 0
211     buoyK (i,j) = 0. _d 0
212 cnh 1.30 maskC (i,j) = 0. _d 0
213 cnh 1.1 ENDDO
214     ENDDO
215    
216 cnh 1.35
217 cnh 1.1 DO bj=myByLo(myThid),myByHi(myThid)
218     DO bi=myBxLo(myThid),myBxHi(myThid)
219    
220 cnh 1.7 C-- Set up work arrays that need valid initial values
221     DO j=1-OLy,sNy+OLy
222     DO i=1-OLx,sNx+OLx
223 cnh 1.27 rTrans(i,j) = 0. _d 0
224     rVel (i,j,1) = 0. _d 0
225     rVel (i,j,2) = 0. _d 0
226 cnh 1.30 fVerT (i,j,1) = 0. _d 0
227     fVerT (i,j,2) = 0. _d 0
228     fVerS (i,j,1) = 0. _d 0
229     fVerS (i,j,2) = 0. _d 0
230     fVerU (i,j,1) = 0. _d 0
231     fVerU (i,j,2) = 0. _d 0
232     fVerV (i,j,1) = 0. _d 0
233     fVerV (i,j,2) = 0. _d 0
234 cnh 1.27 phiHyd(i,j,1) = 0. _d 0
235 cnh 1.30 K13 (i,j,1) = 0. _d 0
236     K23 (i,j,1) = 0. _d 0
237     K33 (i,j,1) = 0. _d 0
238     KapGM (i,j) = GMkbackground
239 cnh 1.7 ENDDO
240     ENDDO
241    
242 cnh 1.1 iMin = 1-OLx+1
243     iMax = sNx+OLx
244     jMin = 1-OLy+1
245     jMax = sNy+OLy
246 cnh 1.35
247 cnh 1.1
248 cnh 1.19 K = 1
249 cnh 1.31 BOTTOM_LAYER = K .EQ. Nr
250 cnh 1.19
251 cnh 1.38 #ifdef DO_PIPELINED_CORRECTION_STEP
252 adcroft 1.4 C-- Calculate gradient of surface pressure
253 cnh 1.28 CALL CALC_GRAD_ETA_SURF(
254 adcroft 1.4 I bi,bj,iMin,iMax,jMin,jMax,
255 cnh 1.29 O etaSurfX,etaSurfY,
256 adcroft 1.4 I myThid)
257     C-- Update fields in top level according to tendency terms
258 adcroft 1.11 CALL CORRECTION_STEP(
259 cnh 1.30 I bi,bj,iMin,iMax,jMin,jMax,K,
260     I etaSurfX,etaSurfY,myTime,myThid)
261 adcroft 1.43 #ifdef ALLOW_OBCS
262 adcroft 1.39 IF (openBoundaries) CALL APPLY_OBCS1( bi, bj, K, myThid )
263 adcroft 1.43 #endif
264 adcroft 1.21 IF ( .NOT. BOTTOM_LAYER ) THEN
265     C-- Update fields in layer below according to tendency terms
266     CALL CORRECTION_STEP(
267 cnh 1.30 I bi,bj,iMin,iMax,jMin,jMax,K+1,
268     I etaSurfX,etaSurfY,myTime,myThid)
269 adcroft 1.43 #ifdef ALLOW_OBCS
270 adcroft 1.39 IF (openBoundaries) CALL APPLY_OBCS1( bi, bj, K+1, myThid )
271 adcroft 1.43 #endif
272 adcroft 1.21 ENDIF
273 cnh 1.38 #endif
274 cnh 1.7 C-- Density of 1st level (below W(1)) reference to level 1
275 cnh 1.38 #ifdef INCLUDE_FIND_RHO_CALL
276 cnh 1.7 CALL FIND_RHO(
277 cnh 1.19 I bi, bj, iMin, iMax, jMin, jMax, K, K, eosType,
278 cnh 1.7 O rhoKm1,
279     I myThid )
280 cnh 1.38 #endif
281 cnh 1.19
282 adcroft 1.42 IF ( (.NOT. BOTTOM_LAYER)
283     #ifdef ALLOW_KPP
284     & .AND. (.NOT.usingKPPmixing) ! CONVECT not needed with KPP mixing
285     #endif
286     & ) THEN
287 cnh 1.19 C-- Check static stability with layer below
288 cnh 1.30 C-- and mix as needed.
289 cnh 1.38 #ifdef INCLUDE_FIND_RHO_CALL
290 cnh 1.19 CALL FIND_RHO(
291     I bi, bj, iMin, iMax, jMin, jMax, K+1, K, eosType,
292     O rhoKp1,
293     I myThid )
294 cnh 1.38 #endif
295     #ifdef INCLUDE_CONVECT_CALL
296 cnh 1.19 CALL CONVECT(
297     I bi,bj,iMin,iMax,jMin,jMax,K+1,rhoKm1,rhoKp1,
298     I myTime,myIter,myThid)
299 cnh 1.38 #endif
300 cnh 1.19 C-- Recompute density after mixing
301 cnh 1.38 #ifdef INCLUDE_FIND_RHO_CALL
302 cnh 1.19 CALL FIND_RHO(
303     I bi, bj, iMin, iMax, jMin, jMax, K, K, eosType,
304     O rhoKm1,
305     I myThid )
306 cnh 1.38 #endif
307 cnh 1.19 ENDIF
308 cnh 1.26 C-- Calculate buoyancy
309 cnh 1.32 CALL CALC_BUOYANCY(
310 cnh 1.26 I bi,bj,iMin,iMax,jMin,jMax,K,rhoKm1,
311     O buoyKm1,
312     I myThid )
313 cnh 1.38 C-- Integrate hydrostatic balance for phiHyd with BC of
314     C-- phiHyd(z=0)=0
315 cnh 1.26 CALL CALC_PHI_HYD(
316     I bi,bj,iMin,iMax,jMin,jMax,K,buoyKm1,buoyKm1,
317     U phiHyd,
318 adcroft 1.5 I myThid )
319    
320 cnh 1.31 DO K=2,Nr
321     BOTTOM_LAYER = K .EQ. Nr
322 cnh 1.38 #ifdef DO_PIPELINED_CORRECTION_STEP
323 adcroft 1.21 IF ( .NOT. BOTTOM_LAYER ) THEN
324     C-- Update fields in layer below according to tendency terms
325     CALL CORRECTION_STEP(
326 cnh 1.30 I bi,bj,iMin,iMax,jMin,jMax,K+1,
327     I etaSurfX,etaSurfY,myTime,myThid)
328 adcroft 1.43 #ifdef ALLOW_OBCS
329 adcroft 1.39 IF (openBoundaries) CALL APPLY_OBCS1( bi, bj, K+1, myThid )
330 adcroft 1.43 #endif
331 adcroft 1.21 ENDIF
332 cnh 1.38 #endif
333 cnh 1.19 C-- Density of K level (below W(K)) reference to K level
334 cnh 1.38 #ifdef INCLUDE_FIND_RHO_CALL
335 cnh 1.19 CALL FIND_RHO(
336     I bi, bj, iMin, iMax, jMin, jMax, K, K, eosType,
337     O rhoK,
338     I myThid )
339 cnh 1.38 #endif
340 adcroft 1.42 IF ( (.NOT. BOTTOM_LAYER)
341     #ifdef ALLOW_KPP
342     & .AND. (.NOT.usingKPPmixing) ! CONVECT not needed with KPP mixing
343     #endif
344     & ) THEN
345 cnh 1.27 C-- Check static stability with layer below and mix as needed.
346     C-- Density of K+1 level (below W(K+1)) reference to K level.
347 cnh 1.38 #ifdef INCLUDE_FIND_RHO_CALL
348 cnh 1.19 CALL FIND_RHO(
349     I bi, bj, iMin, iMax, jMin, jMax, K+1, K, eosType,
350     O rhoKp1,
351     I myThid )
352 cnh 1.38 #endif
353     #ifdef INCLUDE_CONVECT_CALL
354 cnh 1.19 CALL CONVECT(
355     I bi,bj,iMin,iMax,jMin,jMax,K+1,rhoK,rhoKp1,
356     I myTime,myIter,myThid)
357 cnh 1.38 #endif
358 cnh 1.19 C-- Recompute density after mixing
359 cnh 1.38 #ifdef INCLUDE_FIND_RHO_CALL
360 cnh 1.19 CALL FIND_RHO(
361     I bi, bj, iMin, iMax, jMin, jMax, K, K, eosType,
362     O rhoK,
363     I myThid )
364 cnh 1.38 #endif
365 cnh 1.19 ENDIF
366 cnh 1.26 C-- Calculate buoyancy
367 cnh 1.32 CALL CALC_BUOYANCY(
368 cnh 1.26 I bi,bj,iMin,iMax,jMin,jMax,K,rhoK,
369     O buoyK,
370     I myThid )
371 cnh 1.38 C-- Integrate hydrostatic balance for phiHyd with BC of
372     C-- phiHyd(z=0)=0
373 cnh 1.26 CALL CALC_PHI_HYD(
374 cnh 1.30 I bi,bj,iMin,iMax,jMin,jMax,K,buoyKm1,buoyK,
375     U phiHyd,
376     I myThid )
377 cnh 1.19 C-- Calculate iso-neutral slopes for the GM/Redi parameterisation
378 cnh 1.38 #ifdef INCLUDE_FIND_RHO_CALL
379 cnh 1.19 CALL FIND_RHO(
380 cnh 1.30 I bi, bj, iMin, iMax, jMin, jMax, K-1, K, eosType,
381     O rhoTmp,
382     I myThid )
383 cnh 1.38 #endif
384     #ifdef INCLUDE_CALC_ISOSLOPES_CALL
385 cnh 1.19 CALL CALC_ISOSLOPES(
386 cnh 1.30 I bi, bj, iMin, iMax, jMin, jMax, K,
387     I rhoKm1, rhoK, rhotmp,
388     O K13, K23, K33, KapGM,
389     I myThid )
390 cnh 1.38 #endif
391 cnh 1.19 DO J=jMin,jMax
392     DO I=iMin,iMax
393 cnh 1.38 #ifdef INCLUDE_FIND_RHO_CALL
394 cnh 1.27 rhoKm1 (I,J) = rhoK(I,J)
395 cnh 1.38 #endif
396 cnh 1.27 buoyKm1(I,J) = buoyK(I,J)
397 cnh 1.19 ENDDO
398 adcroft 1.10 ENDDO
399 adcroft 1.11 ENDDO ! K
400    
401 adcroft 1.42 #ifdef ALLOW_KPP
402     C-- Compute KPP mixing coefficients
403     IF (usingKPPmixing) THEN
404     CALL TIMER_START('KVMIX (FIND KPP COEFFICIENTS) [DYNAMICS]'
405     I , myThid)
406     CALL KVMIX(
407     I bi, bj, myTime, myThid )
408     CALL TIMER_STOP ('KVMIX (FIND KPP COEFFICIENTS) [DYNAMICS]'
409     I , myThid)
410     ENDIF
411     #endif
412    
413 cnh 1.31 DO K = Nr, 1, -1
414 cnh 1.30
415 cnh 1.1 kM1 =max(1,k-1) ! Points to level above k (=k-1)
416     kUp =1+MOD(k+1,2) ! Cycles through 1,2 to point to layer above
417     kDown=1+MOD(k,2) ! Cycles through 2,1 to point to current layer
418     iMin = 1-OLx+2
419     iMax = sNx+OLx-1
420     jMin = 1-OLy+2
421     jMax = sNy+OLy-1
422    
423     C-- Get temporary terms used by tendency routines
424     CALL CALC_COMMON_FACTORS (
425     I bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,
426 cnh 1.30 O xA,yA,uTrans,vTrans,rTrans,rVel,maskC,maskUp,
427 cnh 1.1 I myThid)
428 cnh 1.38 #ifdef INCLUDE_CALC_DIFFUSIVITY_CALL
429 adcroft 1.12 C-- Calculate the total vertical diffusivity
430     CALL CALC_DIFFUSIVITY(
431     I bi,bj,iMin,iMax,jMin,jMax,K,
432     I maskC,maskUp,KapGM,K33,
433 adcroft 1.42 O KappaRT,KappaRS,KappaRU,KappaRV,
434 adcroft 1.12 I myThid)
435 cnh 1.38 #endif
436 cnh 1.1 C-- Calculate accelerations in the momentum equations
437 cnh 1.9 IF ( momStepping ) THEN
438     CALL CALC_MOM_RHS(
439     I bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,
440 cnh 1.30 I xA,yA,uTrans,vTrans,rTrans,rVel,maskC,
441 adcroft 1.42 I phiHyd,KappaRU,KappaRV,
442 cnh 1.9 U aTerm,xTerm,cTerm,mTerm,pTerm,
443     U fZon, fMer, fVerU, fVerV,
444 cnh 1.38 I myTime, myThid)
445 cnh 1.9 ENDIF
446 cnh 1.1 C-- Calculate active tracer tendencies
447 cnh 1.9 IF ( tempStepping ) THEN
448     CALL CALC_GT(
449     I bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown,
450 cnh 1.30 I xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,
451     I K13,K23,KappaRT,KapGM,
452 cnh 1.9 U aTerm,xTerm,fZon,fMer,fVerT,
453 cnh 1.37 I myTime, myThid)
454 cnh 1.9 ENDIF
455 adcroft 1.18 IF ( saltStepping ) THEN
456     CALL CALC_GS(
457     I bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown,
458 cnh 1.30 I xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,
459     I K13,K23,KappaRS,KapGM,
460 adcroft 1.18 U aTerm,xTerm,fZon,fMer,fVerS,
461 cnh 1.37 I myTime, myThid)
462 adcroft 1.18 ENDIF
463 adcroft 1.11 C-- Prediction step (step forward all model variables)
464     CALL TIMESTEP(
465     I bi,bj,iMin,iMax,jMin,jMax,K,
466     I myThid)
467 adcroft 1.43 #ifdef ALLOW_OBCS
468 adcroft 1.41 C-- Apply open boundary conditions
469 adcroft 1.39 IF (openBoundaries) CALL APPLY_OBCS2( bi, bj, K, myThid )
470 adcroft 1.43 #endif
471 adcroft 1.41 C-- Freeze water
472     IF (allowFreezing)
473     & CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, K, myThid )
474 adcroft 1.11 C-- Diagnose barotropic divergence of predicted fields
475 cnh 1.31 CALL CALC_DIV_GHAT(
476 adcroft 1.11 I bi,bj,iMin,iMax,jMin,jMax,K,
477     I xA,yA,
478     I myThid)
479 adcroft 1.23
480     C-- Cumulative diagnostic calculations (ie. time-averaging)
481 cnh 1.38 #ifdef INCLUDE_DIAGNOSTICS_INTERFACE_CODE
482 adcroft 1.23 IF (taveFreq.GT.0.) THEN
483     CALL DO_TIME_AVERAGES(
484     I myTime, myIter, bi, bj, K, kUp, kDown,
485 cnh 1.30 I K13, K23, rVel, KapGM,
486 adcroft 1.23 I myThid )
487     ENDIF
488     #endif
489 adcroft 1.11
490     ENDDO ! K
491 adcroft 1.12
492     C-- Implicit diffusion
493     IF (implicitDiffusion) THEN
494 adcroft 1.42 IF (tempStepping) CALL IMPLDIFF(
495     I bi, bj, iMin, iMax, jMin, jMax,
496     I deltaTtracer, KappaRT,recip_HFacC,
497     U gTNm1,
498     I myThid )
499     IF (saltStepping) CALL IMPLDIFF(
500     I bi, bj, iMin, iMax, jMin, jMax,
501     I deltaTtracer, KappaRS,recip_HFacC,
502     U gSNm1,
503     I myThid )
504 adcroft 1.44 ENDIF ! implicitDiffusion
505     C-- Implicit viscosity
506     IF (implicitViscosity) THEN
507 adcroft 1.42 IF (momStepping) THEN
508     CALL IMPLDIFF(
509     I bi, bj, iMin, iMax, jMin, jMax,
510     I deltaTmom, KappaRU,recip_HFacW,
511     U gUNm1,
512     I myThid )
513     CALL IMPLDIFF(
514     I bi, bj, iMin, iMax, jMin, jMax,
515     I deltaTmom, KappaRV,recip_HFacS,
516     U gVNm1,
517     I myThid )
518     #ifdef INCLUDE_CD_CODE
519     CALL IMPLDIFF(
520     I bi, bj, iMin, iMax, jMin, jMax,
521     I deltaTmom, KappaRU,recip_HFacW,
522     U vVelD,
523     I myThid )
524     CALL IMPLDIFF(
525     I bi, bj, iMin, iMax, jMin, jMax,
526     I deltaTmom, KappaRV,recip_HFacS,
527     U uVelD,
528     I myThid )
529     #endif
530     ENDIF ! momStepping
531 adcroft 1.44 ENDIF ! implicitViscosity
532 cnh 1.1
533     ENDDO
534     ENDDO
535 adcroft 1.6
536 cnh 1.19 C write(0,*) 'dynamics: pS ',minval(cg2d_x(1:sNx,1:sNy,:,:)),
537     C & maxval(cg2d_x(1:sNx,1:sNy,:,:))
538 cnh 1.20 C write(0,*) 'dynamics: U ',minval(uVel(1:sNx,1:sNy,1,:,:),mask=uVel(1:sNx,1:sNy,1,:,:).NE.0.),
539 adcroft 1.21 C & maxval(uVel(1:sNx,1:sNy,1,:,:),mask=uVel(1:sNx,1:sNy,1,:,:).NE.0.)
540 cnh 1.20 C write(0,*) 'dynamics: V ',minval(vVel(1:sNx,1:sNy,1,:,:),mask=vVel(1:sNx,1:sNy,1,:,:).NE.0.),
541 adcroft 1.21 C & maxval(vVel(1:sNx,1:sNy,1,:,:),mask=vVel(1:sNx,1:sNy,1,:,:).NE.0.)
542 cnh 1.30 C write(0,*) 'dynamics: rVel(1) ',
543     C & minval(rVel(1:sNx,1:sNy,1),mask=rVel(1:sNx,1:sNy,1).NE.0.),
544     C & maxval(rVel(1:sNx,1:sNy,1),mask=rVel(1:sNx,1:sNy,1).NE.0.)
545     C write(0,*) 'dynamics: rVel(2) ',
546     C & minval(rVel(1:sNx,1:sNy,2),mask=rVel(1:sNx,1:sNy,2).NE.0.),
547     C & maxval(rVel(1:sNx,1:sNy,2),mask=rVel(1:sNx,1:sNy,2).NE.0.)
548 adcroft 1.15 cblk write(0,*) 'dynamics: K13',minval(K13(1:sNx,1:sNy,:)),
549     cblk & maxval(K13(1:sNx,1:sNy,:))
550     cblk write(0,*) 'dynamics: K23',minval(K23(1:sNx,1:sNy,:)),
551     cblk & maxval(K23(1:sNx,1:sNy,:))
552     cblk write(0,*) 'dynamics: K33',minval(K33(1:sNx,1:sNy,:)),
553     cblk & maxval(K33(1:sNx,1:sNy,:))
554 cnh 1.19 C write(0,*) 'dynamics: gT ',minval(gT(1:sNx,1:sNy,:,:,:)),
555     C & maxval(gT(1:sNx,1:sNy,:,:,:))
556     C write(0,*) 'dynamics: T ',minval(Theta(1:sNx,1:sNy,:,:,:)),
557     C & maxval(Theta(1:sNx,1:sNy,:,:,:))
558     C write(0,*) 'dynamics: gS ',minval(gS(1:sNx,1:sNy,:,:,:)),
559     C & maxval(gS(1:sNx,1:sNy,:,:,:))
560     C write(0,*) 'dynamics: S ',minval(salt(1:sNx,1:sNy,:,:,:)),
561     C & maxval(salt(1:sNx,1:sNy,:,:,:))
562 cnh 1.31 C write(0,*) 'dynamics: phiHyd ',minval(phiHyd/(Gravity*Rhonil),mask=phiHyd.NE.0.),
563     C & maxval(phiHyd/(Gravity*Rhonil))
564 cnh 1.36 C CALL PLOT_FIELD_XYZRL( gU, ' GU exiting dyanmics ' ,
565     C &Nr, 1, myThid )
566     C CALL PLOT_FIELD_XYZRL( gV, ' GV exiting dyanmics ' ,
567     C &Nr, 1, myThid )
568     C CALL PLOT_FIELD_XYZRL( gS, ' GS exiting dyanmics ' ,
569     C &Nr, 1, myThid )
570     C CALL PLOT_FIELD_XYZRL( gT, ' GT exiting dyanmics ' ,
571 cnh 1.38 C &Nr, 1, myThid )
572     C CALL PLOT_FIELD_XYZRL( phiHyd, ' phiHyd exiting dyanmics ' ,
573 cnh 1.36 C &Nr, 1, myThid )
574    
575 cnh 1.1
576     RETURN
577     END

  ViewVC Help
Powered by ViewVC 1.1.22