/[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.56 - (hide annotations) (download)
Mon Jan 29 20:05:46 2001 UTC (23 years, 3 months ago) by heimbach
Branch: MAIN
Changes since 1.55: +24 -10 lines
Corrected store directives; added one ifdef ALLOW_GMREDI.

1 heimbach 1.56 C $Header: /u/gcmpack/models/MITgcmUV/model/src/dynamics.F,v 1.55 2001/01/08 16:37:43 heimbach 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 heimbach 1.49
33     #ifdef ALLOW_AUTODIFF_TAMC
34 heimbach 1.53 # include "tamc.h"
35     # include "tamc_keys.h"
36 heimbach 1.56 # include "FFIELDS.h"
37     # ifdef ALLOW_KPP
38     # include "KPP.h"
39     # endif
40     # ifdef ALLOW_GMREDI
41     # include "GMREDI.h"
42     # endif
43 heimbach 1.53 #endif /* ALLOW_AUTODIFF_TAMC */
44 heimbach 1.49
45 cnh 1.1 C == Routine arguments ==
46 cnh 1.8 C myTime - Current time in simulation
47     C myIter - Current iteration number in simulation
48 cnh 1.1 C myThid - Thread number for this instance of the routine.
49 cnh 1.8 _RL myTime
50     INTEGER myIter
51 adcroft 1.47 INTEGER myThid
52 cnh 1.1
53     C == Local variables
54     C xA, yA - Per block temporaries holding face areas
55 cnh 1.38 C uTrans, vTrans, rTrans - Per block temporaries holding flow
56     C transport
57 cnh 1.30 C rVel o uTrans: Zonal transport
58 cnh 1.1 C o vTrans: Meridional transport
59 cnh 1.30 C o rTrans: Vertical transport
60 cnh 1.38 C o rVel: Vertical velocity at upper and
61     C lower cell faces.
62 cnh 1.1 C maskC,maskUp o maskC: land/water mask for tracer cells
63     C o maskUp: land/water mask for W points
64     C aTerm, xTerm, cTerm - Work arrays for holding separate terms in
65     C mTerm, pTerm, tendency equations.
66     C fZon, fMer, fVer[STUV] o aTerm: Advection term
67     C o xTerm: Mixing term
68     C o cTerm: Coriolis term
69     C o mTerm: Metric term
70     C o pTerm: Pressure term
71     C o fZon: Zonal flux term
72     C o fMer: Meridional flux term
73     C o fVer: Vertical flux term - note fVer
74     C is "pipelined" in the vertical
75     C so we need an fVer for each
76     C variable.
77 cnh 1.38 C rhoK, rhoKM1 - Density at current level, level above and level
78     C below.
79 cnh 1.26 C rhoKP1
80     C buoyK, buoyKM1 - Buoyancy at current level and level above.
81 cnh 1.31 C phiHyd - Hydrostatic part of the potential phiHydi.
82 cnh 1.38 C In z coords phiHydiHyd is the hydrostatic
83     C pressure anomaly
84     C In p coords phiHydiHyd is the geopotential
85     C surface height
86 cnh 1.30 C anomaly.
87     C etaSurfX, - Holds surface elevation gradient in X and Y.
88     C etaSurfY
89     C KappaRT, - Total diffusion in vertical for T and S.
90 cnh 1.38 C KappaRS (background + spatially varying, isopycnal term).
91 cnh 1.30 C iMin, iMax - Ranges and sub-block indices on which calculations
92     C jMin, jMax are applied.
93 cnh 1.1 C bi, bj
94 heimbach 1.53 C k, kup, - Index for layer above and below. kup and kDown
95     C kDown, km1 are switched with layer to be the appropriate
96 cnh 1.38 C index into fVerTerm.
97 cnh 1.30 _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98     _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
99     _RL uTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
100     _RL vTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
101     _RL rTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
102     _RL rVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
103     _RS maskC (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
104     _RS maskUp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105     _RL aTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
106     _RL xTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
107     _RL cTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
108     _RL mTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
109     _RL pTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
110     _RL fZon (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
111     _RL fMer (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
112     _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
113     _RL fVerS (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
114     _RL fVerU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
115     _RL fVerV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
116 cnh 1.31 _RL phiHyd (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
117 cnh 1.30 _RL rhokm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
118     _RL rhokp1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
119     _RL rhok (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
120     _RL buoyKM1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
121     _RL buoyK (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
122     _RL rhotmp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
123 cnh 1.29 _RL etaSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
124     _RL etaSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
125 cnh 1.31 _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
126     _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
127 adcroft 1.42 _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
128     _RL KappaRV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
129 adcroft 1.50 _RL sigmaX (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
130     _RL sigmaY (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
131     _RL sigmaR (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
132 adcroft 1.12
133 adcroft 1.52 C This is currently also used by IVDC and Diagnostics
134     C #ifdef INCLUDE_CONVECT_CALL
135 adcroft 1.45 _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
136 adcroft 1.52 C #endif
137 adcroft 1.45
138 cnh 1.1 INTEGER iMin, iMax
139     INTEGER jMin, jMax
140     INTEGER bi, bj
141     INTEGER i, j
142 heimbach 1.53 INTEGER k, km1, kup, kDown
143 cnh 1.19 LOGICAL BOTTOM_LAYER
144 cnh 1.1
145 heimbach 1.49 #ifdef ALLOW_AUTODIFF_TAMC
146     INTEGER isbyte
147     PARAMETER( isbyte = 4 )
148    
149     INTEGER act1, act2, act3, act4
150     INTEGER max1, max2, max3
151 heimbach 1.51 INTEGER iikey, kkey
152 heimbach 1.49 INTEGER maximpl
153 heimbach 1.53 #endif /* ALLOW_AUTODIFF_TAMC */
154 heimbach 1.49
155 adcroft 1.11 C--- The algorithm...
156     C
157     C "Correction Step"
158     C =================
159     C Here we update the horizontal velocities with the surface
160     C pressure such that the resulting flow is either consistent
161     C with the free-surface evolution or the rigid-lid:
162     C U[n] = U* + dt x d/dx P
163     C V[n] = V* + dt x d/dy P
164     C
165     C "Calculation of Gs"
166     C ===================
167     C This is where all the accelerations and tendencies (ie.
168 heimbach 1.53 C physics, parameterizations etc...) are calculated
169 cnh 1.27 C rVel = sum_r ( div. u[n] )
170 adcroft 1.11 C rho = rho ( theta[n], salt[n] )
171 cnh 1.27 C b = b(rho, theta)
172 adcroft 1.11 C K31 = K31 ( rho )
173 cnh 1.27 C Gu[n] = Gu( u[n], v[n], rVel, b, ... )
174     C Gv[n] = Gv( u[n], v[n], rVel, b, ... )
175     C Gt[n] = Gt( theta[n], u[n], v[n], rVel, K31, ... )
176     C Gs[n] = Gs( salt[n], u[n], v[n], rVel, K31, ... )
177 adcroft 1.11 C
178 adcroft 1.12 C "Time-stepping" or "Prediction"
179 adcroft 1.11 C ================================
180     C The models variables are stepped forward with the appropriate
181     C time-stepping scheme (currently we use Adams-Bashforth II)
182     C - For momentum, the result is always *only* a "prediction"
183     C in that the flow may be divergent and will be "corrected"
184     C later with a surface pressure gradient.
185     C - Normally for tracers the result is the new field at time
186     C level [n+1} *BUT* in the case of implicit diffusion the result
187     C is also *only* a prediction.
188     C - We denote "predictors" with an asterisk (*).
189     C U* = U[n] + dt x ( 3/2 Gu[n] - 1/2 Gu[n-1] )
190     C V* = V[n] + dt x ( 3/2 Gv[n] - 1/2 Gv[n-1] )
191     C theta[n+1] = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
192     C salt[n+1] = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
193 adcroft 1.12 C With implicit diffusion:
194 adcroft 1.11 C theta* = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
195     C salt* = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
196 adcroft 1.12 C (1 + dt * K * d_zz) theta[n] = theta*
197     C (1 + dt * K * d_zz) salt[n] = salt*
198 adcroft 1.11 C---
199    
200 heimbach 1.49 #ifdef ALLOW_AUTODIFF_TAMC
201     C-- dummy statement to end declaration part
202     ikey = 1
203 heimbach 1.53 #endif /* ALLOW_AUTODIFF_TAMC */
204 heimbach 1.49
205 heimbach 1.54
206 cnh 1.1 C-- Set up work arrays with valid (i.e. not NaN) values
207     C These inital values do not alter the numerical results. They
208     C just ensure that all memory references are to valid floating
209     C point numbers. This prevents spurious hardware signals due to
210     C uninitialised but inert locations.
211     DO j=1-OLy,sNy+OLy
212     DO i=1-OLx,sNx+OLx
213 adcroft 1.5 xA(i,j) = 0. _d 0
214     yA(i,j) = 0. _d 0
215     uTrans(i,j) = 0. _d 0
216     vTrans(i,j) = 0. _d 0
217     aTerm(i,j) = 0. _d 0
218     xTerm(i,j) = 0. _d 0
219     cTerm(i,j) = 0. _d 0
220     mTerm(i,j) = 0. _d 0
221     pTerm(i,j) = 0. _d 0
222     fZon(i,j) = 0. _d 0
223     fMer(i,j) = 0. _d 0
224 heimbach 1.53 DO k=1,Nr
225 cnh 1.31 phiHyd (i,j,k) = 0. _d 0
226 adcroft 1.45 KappaRU(i,j,k) = 0. _d 0
227     KappaRV(i,j,k) = 0. _d 0
228 adcroft 1.50 sigmaX(i,j,k) = 0. _d 0
229     sigmaY(i,j,k) = 0. _d 0
230     sigmaR(i,j,k) = 0. _d 0
231 cnh 1.1 ENDDO
232 cnh 1.30 rhoKM1 (i,j) = 0. _d 0
233     rhok (i,j) = 0. _d 0
234     rhoKP1 (i,j) = 0. _d 0
235     rhoTMP (i,j) = 0. _d 0
236 cnh 1.26 buoyKM1(i,j) = 0. _d 0
237     buoyK (i,j) = 0. _d 0
238 cnh 1.30 maskC (i,j) = 0. _d 0
239 cnh 1.1 ENDDO
240     ENDDO
241    
242 cnh 1.35
243 heimbach 1.49 #ifdef ALLOW_AUTODIFF_TAMC
244     C-- HPF directive to help TAMC
245 heimbach 1.53 CHPF$ INDEPENDENT
246     #endif /* ALLOW_AUTODIFF_TAMC */
247 heimbach 1.49
248 cnh 1.1 DO bj=myByLo(myThid),myByHi(myThid)
249 heimbach 1.49
250     #ifdef ALLOW_AUTODIFF_TAMC
251     C-- HPF directive to help TAMC
252 heimbach 1.53 CHPF$ INDEPENDENT, NEW (rTrans,rVel,fVerT,fVerS,fVerU,fVerV
253     CHPF$& ,phiHyd,utrans,vtrans,maskc,xA,yA
254     CHPF$& ,KappaRT,KappaRS,KappaRU,KappaRV
255     CHPF$& )
256     #endif /* ALLOW_AUTODIFF_TAMC */
257 heimbach 1.49
258 cnh 1.1 DO bi=myBxLo(myThid),myBxHi(myThid)
259    
260 heimbach 1.49 #ifdef ALLOW_AUTODIFF_TAMC
261     act1 = bi - myBxLo(myThid)
262     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
263    
264     act2 = bj - myByLo(myThid)
265     max2 = myByHi(myThid) - myByLo(myThid) + 1
266    
267     act3 = myThid - 1
268     max3 = nTx*nTy
269    
270     act4 = ikey_dynamics - 1
271    
272     ikey = (act1 + 1) + act2*max1
273     & + act3*max1*max2
274     & + act4*max1*max2*max3
275 heimbach 1.53 #endif /* ALLOW_AUTODIFF_TAMC */
276 heimbach 1.49
277 cnh 1.7 C-- Set up work arrays that need valid initial values
278     DO j=1-OLy,sNy+OLy
279     DO i=1-OLx,sNx+OLx
280 cnh 1.27 rTrans(i,j) = 0. _d 0
281     rVel (i,j,1) = 0. _d 0
282     rVel (i,j,2) = 0. _d 0
283 cnh 1.30 fVerT (i,j,1) = 0. _d 0
284     fVerT (i,j,2) = 0. _d 0
285     fVerS (i,j,1) = 0. _d 0
286     fVerS (i,j,2) = 0. _d 0
287     fVerU (i,j,1) = 0. _d 0
288     fVerU (i,j,2) = 0. _d 0
289     fVerV (i,j,1) = 0. _d 0
290     fVerV (i,j,2) = 0. _d 0
291 cnh 1.27 phiHyd(i,j,1) = 0. _d 0
292 cnh 1.7 ENDDO
293     ENDDO
294    
295 adcroft 1.45 DO k=1,Nr
296     DO j=1-OLy,sNy+OLy
297     DO i=1-OLx,sNx+OLx
298     #ifdef INCLUDE_CONVECT_CALL
299     ConvectCount(i,j,k) = 0.
300     #endif
301     KappaRT(i,j,k) = 0. _d 0
302     KappaRS(i,j,k) = 0. _d 0
303     ENDDO
304     ENDDO
305     ENDDO
306    
307 cnh 1.1 iMin = 1-OLx+1
308     iMax = sNx+OLx
309     jMin = 1-OLy+1
310     jMax = sNy+OLy
311 cnh 1.35
312 heimbach 1.53 k = 1
313     BOTTOM_LAYER = k .EQ. Nr
314 cnh 1.19
315 cnh 1.38 #ifdef DO_PIPELINED_CORRECTION_STEP
316 adcroft 1.4 C-- Calculate gradient of surface pressure
317 cnh 1.28 CALL CALC_GRAD_ETA_SURF(
318 adcroft 1.4 I bi,bj,iMin,iMax,jMin,jMax,
319 cnh 1.29 O etaSurfX,etaSurfY,
320 adcroft 1.4 I myThid)
321     C-- Update fields in top level according to tendency terms
322 adcroft 1.11 CALL CORRECTION_STEP(
323 heimbach 1.53 I bi,bj,iMin,iMax,jMin,jMax,k,
324 cnh 1.30 I etaSurfX,etaSurfY,myTime,myThid)
325 heimbach 1.49
326 adcroft 1.43 #ifdef ALLOW_OBCS
327 heimbach 1.49 IF (openBoundaries) THEN
328     #ifdef ALLOW_AUTODIFF_TAMC
329 heimbach 1.53 CADJ STORE uvel (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
330     CADJ STORE vvel (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
331     CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
332     CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
333     #endif /* ALLOW_AUTODIFF_TAMC */
334     CALL APPLY_OBCS1( bi, bj, k, myThid )
335 heimbach 1.49 END IF
336 adcroft 1.43 #endif
337 heimbach 1.49
338 adcroft 1.21 IF ( .NOT. BOTTOM_LAYER ) THEN
339     C-- Update fields in layer below according to tendency terms
340     CALL CORRECTION_STEP(
341 heimbach 1.53 I bi,bj,iMin,iMax,jMin,jMax,k+1,
342 cnh 1.30 I etaSurfX,etaSurfY,myTime,myThid)
343 adcroft 1.43 #ifdef ALLOW_OBCS
344 heimbach 1.49 IF (openBoundaries) THEN
345     #ifdef ALLOW_AUTODIFF_TAMC
346 heimbach 1.53 CADJ STORE uvel (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
347     CADJ STORE vvel (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
348     CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
349     CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
350     #endif /* ALLOW_AUTODIFF_TAMC */
351     CALL APPLY_OBCS1( bi, bj, k+1, myThid )
352 heimbach 1.49 END IF
353 adcroft 1.43 #endif
354 adcroft 1.21 ENDIF
355 cnh 1.38 #endif
356 heimbach 1.54
357 cnh 1.7 C-- Density of 1st level (below W(1)) reference to level 1
358 cnh 1.38 #ifdef INCLUDE_FIND_RHO_CALL
359 heimbach 1.49 #ifdef ALLOW_AUTODIFF_TAMC
360 heimbach 1.53 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
361     CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
362     #endif /* ALLOW_AUTODIFF_TAMC */
363 cnh 1.7 CALL FIND_RHO(
364 heimbach 1.53 I bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,
365 cnh 1.7 O rhoKm1,
366     I myThid )
367 cnh 1.38 #endif
368 cnh 1.19
369 heimbach 1.53 IF (.NOT. BOTTOM_LAYER) THEN
370    
371 cnh 1.19 C-- Check static stability with layer below
372 cnh 1.30 C-- and mix as needed.
373 cnh 1.38 #ifdef INCLUDE_FIND_RHO_CALL
374 heimbach 1.49 #ifdef ALLOW_AUTODIFF_TAMC
375 heimbach 1.53 CADJ STORE theta(:,:,k+1,bi,bj) = comlev1_bibj
376     CADJ & , key = ikey, byte = isbyte
377     CADJ STORE salt (:,:,k+1,bi,bj) = comlev1_bibj
378     CADJ & , key = ikey, byte = isbyte
379     #endif /* ALLOW_AUTODIFF_TAMC */
380 cnh 1.19 CALL FIND_RHO(
381 heimbach 1.53 I bi, bj, iMin, iMax, jMin, jMax, k+1, k, eosType,
382 cnh 1.19 O rhoKp1,
383     I myThid )
384 cnh 1.38 #endif
385 heimbach 1.49
386     #ifdef ALLOW_AUTODIFF_TAMC
387 heimbach 1.53 CADJ STORE rhoKm1(:,:) = comlev1_bibj, key = ikey, byte = isbyte
388     CADJ STORE rhoKp1(:,:) = comlev1_bibj, key = ikey, byte = isbyte
389     #endif /* ALLOW_AUTODIFF_TAMC */
390 heimbach 1.54
391     #ifdef INCLUDE_CONVECT_CALL
392    
393 cnh 1.19 CALL CONVECT(
394 heimbach 1.53 I bi,bj,iMin,iMax,jMin,jMax,k+1,rhoKm1,rhoKp1,
395 adcroft 1.45 U ConvectCount,
396 cnh 1.19 I myTime,myIter,myThid)
397 heimbach 1.53
398 heimbach 1.49 #ifdef ALLOW_AUTODIFF_TAMC
399     CADJ STORE theta(:,:,k+1,bi,bj),theta(:,:,k,bi,bj)
400 heimbach 1.53 CADJ & = comlev1_bibj, key = ikey, byte = isbyte
401 heimbach 1.49 CADJ STORE salt (:,:,k+1,bi,bj),salt (:,:,k,bi,bj)
402 heimbach 1.53 CADJ & = comlev1_bibj, key = ikey, byte = isbyte
403     #endif /* ALLOW_AUTODIFF_TAMC */
404 heimbach 1.49
405     #endif
406    
407 adcroft 1.45 C-- Implicit Vertical Diffusion for Convection
408 heimbach 1.53 IF (ivdc_kappa.NE.0.) THEN
409     CALL CALC_IVDC(
410     I bi,bj,iMin,iMax,jMin,jMax,k+1,rhoKm1,rhoKp1,
411 adcroft 1.45 U ConvectCount, KappaRT, KappaRS,
412     I myTime,myIter,myThid)
413 heimbach 1.53 ENDIF
414 heimbach 1.49
415 cnh 1.19 C-- Recompute density after mixing
416 cnh 1.38 #ifdef INCLUDE_FIND_RHO_CALL
417 cnh 1.19 CALL FIND_RHO(
418 heimbach 1.53 I bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,
419 cnh 1.19 O rhoKm1,
420     I myThid )
421 cnh 1.38 #endif
422 cnh 1.19 ENDIF
423 heimbach 1.53
424 cnh 1.26 C-- Calculate buoyancy
425 cnh 1.32 CALL CALC_BUOYANCY(
426 heimbach 1.53 I bi,bj,iMin,iMax,jMin,jMax,k,rhoKm1,
427 cnh 1.26 O buoyKm1,
428     I myThid )
429 heimbach 1.53
430 cnh 1.38 C-- Integrate hydrostatic balance for phiHyd with BC of
431     C-- phiHyd(z=0)=0
432 cnh 1.26 CALL CALC_PHI_HYD(
433 heimbach 1.53 I bi,bj,iMin,iMax,jMin,jMax,k,buoyKm1,buoyKm1,
434 cnh 1.26 U phiHyd,
435 adcroft 1.5 I myThid )
436 heimbach 1.53
437     #ifdef ALLOW_GMREDI
438 heimbach 1.54 IF ( useGMRedi ) THEN
439 adcroft 1.50 CALL GRAD_SIGMA(
440 heimbach 1.53 I bi, bj, iMin, iMax, jMin, jMax, k,
441 adcroft 1.50 I rhoKm1, rhoKm1, rhoKm1,
442     O sigmaX, sigmaY, sigmaR,
443     I myThid )
444 heimbach 1.55 ELSE
445     DO j=1-OLy,sNy+OLy
446     DO i=1-OLx,sNx+OLx
447     sigmaX(i,j,k) = 0. _d 0
448     sigmaY(i,j,k) = 0. _d 0
449     sigmaR(i,j,k) = 0. _d 0
450     ENDDO
451     ENDDO
452 heimbach 1.54 ENDIF
453 heimbach 1.53 #endif
454 adcroft 1.5
455 adcroft 1.50 C-- Start of downward loop
456 heimbach 1.53 DO k=2,Nr
457 heimbach 1.49
458     #ifdef ALLOW_AUTODIFF_TAMC
459 heimbach 1.51 kkey = (ikey-1)*(Nr-2+1) + (k-2) + 1
460 heimbach 1.53 #endif /* ALLOW_AUTODIFF_TAMC */
461 heimbach 1.49
462 heimbach 1.53 BOTTOM_LAYER = k .EQ. Nr
463 heimbach 1.49
464 cnh 1.38 #ifdef DO_PIPELINED_CORRECTION_STEP
465 adcroft 1.21 IF ( .NOT. BOTTOM_LAYER ) THEN
466     C-- Update fields in layer below according to tendency terms
467     CALL CORRECTION_STEP(
468 heimbach 1.53 I bi,bj,iMin,iMax,jMin,jMax,k+1,
469 cnh 1.30 I etaSurfX,etaSurfY,myTime,myThid)
470 adcroft 1.43 #ifdef ALLOW_OBCS
471 heimbach 1.49 IF (openBoundaries) THEN
472     #ifdef ALLOW_AUTODIFF_TAMC
473 heimbach 1.53 CADJ STORE uvel (:,:,k,bi,bj) = comlev1_bibj_k
474     CADJ & , key = kkey, byte = isbyte
475     CADJ STORE vvel (:,:,k,bi,bj) = comlev1_bibj_k
476     CADJ & , key = kkey, byte = isbyte
477     CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k
478     CADJ & , key = kkey, byte = isbyte
479     CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k
480     CADJ & , key = kkey, byte = isbyte
481     #endif /* ALLOW_AUTODIFF_TAMC */
482     CALL APPLY_OBCS1( bi, bj, k+1, myThid )
483 heimbach 1.49 END IF
484 adcroft 1.43 #endif
485 adcroft 1.21 ENDIF
486 heimbach 1.53 #endif /* DO_PIPELINED_CORRECTION_STEP */
487 heimbach 1.49
488 heimbach 1.53 C-- Density of k level (below W(k)) reference to k level
489 cnh 1.38 #ifdef INCLUDE_FIND_RHO_CALL
490 heimbach 1.49 #ifdef ALLOW_AUTODIFF_TAMC
491 heimbach 1.53 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k
492     CADJ & , key = kkey, byte = isbyte
493     CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k
494     CADJ & , key = kkey, byte = isbyte
495     #endif /* ALLOW_AUTODIFF_TAMC */
496 cnh 1.19 CALL FIND_RHO(
497 heimbach 1.53 I bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,
498 cnh 1.19 O rhoK,
499     I myThid )
500 heimbach 1.53
501     #ifdef ALLOW_AUTODIFF_TAMC
502 heimbach 1.54 cph( storing not necessary
503     cphCADJ STORE rhoK(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
504     cph)
505 heimbach 1.53 #endif /* ALLOW_AUTODIFF_TAMC */
506 cnh 1.38 #endif
507 heimbach 1.53
508     IF (.NOT. BOTTOM_LAYER) THEN
509    
510 cnh 1.27 C-- Check static stability with layer below and mix as needed.
511 heimbach 1.53 C-- Density of k+1 level (below W(k+1)) reference to k level.
512 cnh 1.38 #ifdef INCLUDE_FIND_RHO_CALL
513 heimbach 1.49 #ifdef ALLOW_AUTODIFF_TAMC
514 heimbach 1.53 CADJ STORE theta(:,:,k+1,bi,bj) = comlev1_bibj_k
515     CADJ & , key = kkey, byte = isbyte
516     CADJ STORE salt (:,:,k+1,bi,bj) = comlev1_bibj_k
517     CADJ & , key = kkey, byte = isbyte
518     #endif /* ALLOW_AUTODIFF_TAMC */
519 cnh 1.19 CALL FIND_RHO(
520 heimbach 1.53 I bi, bj, iMin, iMax, jMin, jMax, k+1, k, eosType,
521 cnh 1.19 O rhoKp1,
522     I myThid )
523 heimbach 1.49 #ifdef ALLOW_AUTODIFF_TAMC
524 heimbach 1.53 CADJ STORE rhoKp1(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
525     #endif /* ALLOW_AUTODIFF_TAMC */
526 heimbach 1.49 #endif
527    
528 cnh 1.38 #ifdef INCLUDE_CONVECT_CALL
529 cnh 1.19 CALL CONVECT(
530 heimbach 1.53 I bi,bj,iMin,iMax,jMin,jMax,k+1,rhoK,rhoKp1,
531 adcroft 1.45 U ConvectCount,
532 cnh 1.19 I myTime,myIter,myThid)
533 heimbach 1.53
534 cnh 1.38 #endif
535 heimbach 1.49
536 adcroft 1.45 C-- Implicit Vertical Diffusion for Convection
537 heimbach 1.49 IF (ivdc_kappa.NE.0.) THEN
538 heimbach 1.53 #ifdef ALLOW_AUTODIFF_TAMC
539     CADJ STORE rhoKm1(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
540     #endif /* ALLOW_AUTODIFF_TAMC */
541 heimbach 1.49 CALL CALC_IVDC(
542 heimbach 1.53 I bi,bj,iMin,iMax,jMin,jMax,k+1,rhoKm1,rhoKp1,
543 adcroft 1.45 U ConvectCount, KappaRT, KappaRS,
544     I myTime,myIter,myThid)
545 heimbach 1.49 END IF
546    
547 cnh 1.19 C-- Recompute density after mixing
548 cnh 1.38 #ifdef INCLUDE_FIND_RHO_CALL
549 heimbach 1.53 #ifdef ALLOW_AUTODIFF_TAMC
550     CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k
551     CADJ & , key = kkey, byte = isbyte
552     CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k
553     CADJ & , key = kkey, byte = isbyte
554     #endif /* ALLOW_AUTODIFF_TAMC */
555 cnh 1.19 CALL FIND_RHO(
556 heimbach 1.53 I bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,
557 cnh 1.19 O rhoK,
558     I myThid )
559 cnh 1.38 #endif
560 heimbach 1.53
561     C-- IF (.NOT. BOTTOM_LAYER) ends here
562 cnh 1.19 ENDIF
563 heimbach 1.53
564 cnh 1.26 C-- Calculate buoyancy
565 cnh 1.32 CALL CALC_BUOYANCY(
566 heimbach 1.53 I bi,bj,iMin,iMax,jMin,jMax,k,rhoK,
567 cnh 1.26 O buoyK,
568     I myThid )
569 heimbach 1.53
570 cnh 1.38 C-- Integrate hydrostatic balance for phiHyd with BC of
571     C-- phiHyd(z=0)=0
572 cnh 1.26 CALL CALC_PHI_HYD(
573 heimbach 1.53 I bi,bj,iMin,iMax,jMin,jMax,k,buoyKm1,buoyK,
574 cnh 1.30 U phiHyd,
575     I myThid )
576 heimbach 1.53
577     #ifdef INCLUDE_FIND_RHO_CALL
578 cnh 1.19 C-- Calculate iso-neutral slopes for the GM/Redi parameterisation
579 heimbach 1.53
580     #ifdef ALLOW_AUTODIFF_TAMC
581     CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k
582     CADJ & , key = kkey, byte = isbyte
583     CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k
584     CADJ & , key = kkey, byte = isbyte
585     #endif /* ALLOW_AUTODIFF_TAMC */
586    
587 cnh 1.19 CALL FIND_RHO(
588 heimbach 1.53 I bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,
589 cnh 1.30 O rhoTmp,
590     I myThid )
591 cnh 1.38 #endif
592 heimbach 1.53
593    
594     #ifdef ALLOW_GMREDI
595 heimbach 1.54 IF ( useGMRedi ) THEN
596 adcroft 1.50 CALL GRAD_SIGMA(
597 heimbach 1.53 I bi, bj, iMin, iMax, jMin, jMax, k,
598 adcroft 1.50 I rhoK, rhotmp, rhoK,
599     O sigmaX, sigmaY, sigmaR,
600     I myThid )
601 heimbach 1.55 ELSE
602     DO j=1-OLy,sNy+OLy
603     DO i=1-OLx,sNx+OLx
604     sigmaX(i,j,k) = 0. _d 0
605     sigmaY(i,j,k) = 0. _d 0
606     sigmaR(i,j,k) = 0. _d 0
607     ENDDO
608     ENDDO
609 heimbach 1.54 ENDIF
610 heimbach 1.53 #endif
611 heimbach 1.49
612 cnh 1.19 DO J=jMin,jMax
613     DO I=iMin,iMax
614 cnh 1.38 #ifdef INCLUDE_FIND_RHO_CALL
615 cnh 1.27 rhoKm1 (I,J) = rhoK(I,J)
616 cnh 1.38 #endif
617 cnh 1.27 buoyKm1(I,J) = buoyK(I,J)
618 cnh 1.19 ENDDO
619 adcroft 1.10 ENDDO
620 heimbach 1.53
621     C-- end of k loop
622 heimbach 1.49 ENDDO
623    
624 heimbach 1.54 C Determines forcing terms based on external fields
625     C relaxation terms, etc.
626     CALL EXTERNAL_FORCING_SURF(
627     I bi, bj, iMin, iMax, jMin, jMax,
628     I myThid )
629    
630 heimbach 1.55 #ifdef ALLOW_AUTODIFF_TAMC
631    
632     CADJ STORE surfacetendencyu(:,:,bi,bj)
633     CADJ & , surfacetendencyv(:,:,bi,bj)
634     CADJ & , surfacetendencys(:,:,bi,bj)
635     CADJ & , surfacetendencyt(:,:,bi,bj)
636     CADJ & = comlev1_bibj, key=ikey, byte=isbyte
637    
638     # ifdef ALLOW_GMREDI
639     CADJ STORE sigmaX(:,:,:) = comlev1, key=ikey, byte=isbyte
640     CADJ STORE sigmaY(:,:,:) = comlev1, key=ikey, byte=isbyte
641     CADJ STORE sigmaR(:,:,:) = comlev1, key=ikey, byte=isbyte
642     # endif /* ALLOW_GMREDI */
643    
644     #endif /* ALLOW_AUTODIFF_TAMC */
645    
646 adcroft 1.50 #ifdef ALLOW_GMREDI
647 heimbach 1.53 IF (useGMRedi) THEN
648     DO k=1, Nr
649     CALL GMREDI_CALC_TENSOR(
650     I bi, bj, iMin, iMax, jMin, jMax, k,
651 adcroft 1.50 I sigmaX, sigmaY, sigmaR,
652     I myThid )
653 heimbach 1.53 ENDDO
654 heimbach 1.55 #ifdef ALLOW_AUTODIFF_TAMC
655     ELSE
656     DO k=1, Nr
657     CALL GMREDI_CALC_TENSOR_DUMMY(
658     I bi, bj, iMin, iMax, jMin, jMax, k,
659     I sigmaX, sigmaY, sigmaR,
660     I myThid )
661     ENDDO
662     #endif /* ALLOW_AUTODIFF_TAMC */
663 heimbach 1.53 ENDIF
664 adcroft 1.50 #endif
665    
666 heimbach 1.49 #ifdef ALLOW_AUTODIFF_TAMC
667 heimbach 1.55 CADJ STORE KappaRT(:,:,:) = comlev1_bibj, key=ikey, byte=isbyte
668     CADJ STORE KappaRS(:,:,:) = comlev1_bibj, key=ikey, byte=isbyte
669 heimbach 1.53
670 heimbach 1.56 #ifdef ALLOW_GMREDI
671     C-- R.G. We need to define a new tape since Kw use mythid instead of bi,bj
672     CADJ STORE Kwx(:,:,:,myThid) = comlev1_bibj, key=ikey, byte=isbyte
673     CADJ STORE Kwy(:,:,:,myThid) = comlev1_bibj, key=ikey, byte=isbyte
674     CADJ STORE Kwz(:,:,:,myThid) = comlev1_bibj, key=ikey, byte=isbyte
675     #endif
676    
677     CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
678     CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
679     CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
680     CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
681     =======
682 heimbach 1.55 C-- R.G. We need to define a new tape since Kw use mythid instead of bi,bj
683     CADJ STORE Kwx(:,:,:,myThid) = comlev1_bibj, key=ikey, byte=isbyte
684     CADJ STORE Kwy(:,:,:,myThid) = comlev1_bibj, key=ikey, byte=isbyte
685     CADJ STORE Kwz(:,:,:,myThid) = comlev1_bibj, key=ikey, byte=isbyte
686    
687     CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
688     CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
689     CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
690     CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
691 heimbach 1.53
692     C-- dummy initialization to break data flow because
693     C-- calc_div_ghat has a condition for initialization
694     DO J=jMin,jMax
695     DO I=iMin,iMax
696     cg2d_b(i,j,bi,bj) = 0.0
697     ENDDO
698     ENDDO
699     #endif /* ALLOW_AUTODIFF_TAMC */
700    
701 adcroft 1.42 #ifdef ALLOW_KPP
702 heimbach 1.53 C-- Compute KPP mixing coefficients
703     IF (useKPP) THEN
704    
705     CALL TIMER_START('KPP_CALC [DYNAMICS]', myThid)
706     CALL KPP_CALC(
707 heimbach 1.54 I bi, bj, myTime, myThid )
708 heimbach 1.53 CALL TIMER_STOP ('KPP_CALC [DYNAMICS]', myThid)
709    
710     #ifdef ALLOW_AUTODIFF_TAMC
711     ELSE
712     DO j=1-OLy,sNy+OLy
713     DO i=1-OLx,sNx+OLx
714     KPPhbl (i,j,bi,bj) = 1.0
715     KPPfrac(i,j,bi,bj) = 0.0
716     DO k = 1,Nr
717     KPPghat (i,j,k,bi,bj) = 0.0
718     KPPviscAz (i,j,k,bi,bj) = viscAz
719     KPPdiffKzT(i,j,k,bi,bj) = diffKzT
720     KPPdiffKzS(i,j,k,bi,bj) = diffKzS
721     ENDDO
722     ENDDO
723     ENDDO
724     #endif /* ALLOW_AUTODIFF_TAMC */
725     ENDIF
726    
727     #ifdef ALLOW_AUTODIFF_TAMC
728     CADJ STORE KPPghat (:,:,:,bi,bj)
729     CADJ & , KPPviscAz (:,:,:,bi,bj)
730     CADJ & , KPPdiffKzT(:,:,:,bi,bj)
731     CADJ & , KPPdiffKzS(:,:,:,bi,bj)
732     CADJ & , KPPfrac (:,: ,bi,bj)
733 heimbach 1.55 CADJ & = comlev1_bibj, key=ikey, byte=isbyte
734 heimbach 1.53 #endif /* ALLOW_AUTODIFF_TAMC */
735    
736     #endif /* ALLOW_KPP */
737 adcroft 1.42
738 adcroft 1.50 C-- Start of upward loop
739 heimbach 1.53 DO k = Nr, 1, -1
740 cnh 1.30
741 heimbach 1.53 C-- km1 Points to level above k (=k-1)
742     C-- kup Cycles through 1,2 to point to layer above
743     C-- kDown Cycles through 2,1 to point to current layer
744    
745     km1 =max(1,k-1)
746     kup =1+MOD(k+1,2)
747     kDown=1+MOD(k,2)
748 heimbach 1.49
749 cnh 1.1 iMin = 1-OLx+2
750     iMax = sNx+OLx-1
751     jMin = 1-OLy+2
752     jMax = sNy+OLy-1
753    
754 heimbach 1.49 #ifdef ALLOW_AUTODIFF_TAMC
755 heimbach 1.51 kkey = (ikey-1)*(Nr-1+1) + (k-1) + 1
756 heimbach 1.56 CADJ STORE rvel (:,:,kdown) = comlev1_bibj_k, key=kkey, byte=isbyte
757     CADJ STORE rTrans(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
758     CADJ STORE KappaRT(:,:,k) = comlev1_bibj_k, key=kkey, byte=isbyte
759     CADJ STORE KappaRS(:,:,k) = comlev1_bibj_k, key=kkey, byte=isbyte
760 heimbach 1.53 #endif /* ALLOW_AUTODIFF_TAMC */
761 heimbach 1.49
762 cnh 1.1 C-- Get temporary terms used by tendency routines
763     CALL CALC_COMMON_FACTORS (
764 heimbach 1.53 I bi,bj,iMin,iMax,jMin,jMax,k,km1,kup,kDown,
765 cnh 1.30 O xA,yA,uTrans,vTrans,rTrans,rVel,maskC,maskUp,
766 cnh 1.1 I myThid)
767 heimbach 1.49
768 adcroft 1.47 #ifdef ALLOW_OBCS
769     IF (openBoundaries) THEN
770 heimbach 1.53 CALL APPLY_OBCS3( bi, bj, k, kup, rTrans, rVel, myThid )
771 adcroft 1.47 ENDIF
772     #endif
773 heimbach 1.49
774 cnh 1.38 #ifdef INCLUDE_CALC_DIFFUSIVITY_CALL
775 adcroft 1.12 C-- Calculate the total vertical diffusivity
776     CALL CALC_DIFFUSIVITY(
777 heimbach 1.53 I bi,bj,iMin,iMax,jMin,jMax,k,
778 adcroft 1.50 I maskC,maskUp,
779 adcroft 1.42 O KappaRT,KappaRS,KappaRU,KappaRV,
780 adcroft 1.12 I myThid)
781 cnh 1.38 #endif
782 cnh 1.1 C-- Calculate accelerations in the momentum equations
783 cnh 1.9 IF ( momStepping ) THEN
784     CALL CALC_MOM_RHS(
785 heimbach 1.53 I bi,bj,iMin,iMax,jMin,jMax,k,km1,kup,kDown,
786 cnh 1.30 I xA,yA,uTrans,vTrans,rTrans,rVel,maskC,
787 adcroft 1.42 I phiHyd,KappaRU,KappaRV,
788 cnh 1.9 U aTerm,xTerm,cTerm,mTerm,pTerm,
789     U fZon, fMer, fVerU, fVerV,
790 cnh 1.38 I myTime, myThid)
791 heimbach 1.49 #ifdef ALLOW_AUTODIFF_TAMC
792     #ifdef INCLUDE_CD_CODE
793     ELSE
794     DO j=1-OLy,sNy+OLy
795     DO i=1-OLx,sNx+OLx
796     guCD(i,j,k,bi,bj) = 0.0
797     gvCD(i,j,k,bi,bj) = 0.0
798     END DO
799     END DO
800     #endif
801 heimbach 1.53 #endif /* ALLOW_AUTODIFF_TAMC */
802 cnh 1.9 ENDIF
803 cnh 1.1 C-- Calculate active tracer tendencies
804 cnh 1.9 IF ( tempStepping ) THEN
805     CALL CALC_GT(
806 heimbach 1.53 I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
807 cnh 1.30 I xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,
808 adcroft 1.50 I KappaRT,
809 cnh 1.9 U aTerm,xTerm,fZon,fMer,fVerT,
810 cnh 1.37 I myTime, myThid)
811 cnh 1.9 ENDIF
812 adcroft 1.18 IF ( saltStepping ) THEN
813     CALL CALC_GS(
814 heimbach 1.53 I bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
815 cnh 1.30 I xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,
816 adcroft 1.50 I KappaRS,
817 adcroft 1.18 U aTerm,xTerm,fZon,fMer,fVerS,
818 cnh 1.37 I myTime, myThid)
819 adcroft 1.18 ENDIF
820 adcroft 1.47 #ifdef ALLOW_OBCS
821     C-- Calculate future values on open boundaries
822     IF (openBoundaries) THEN
823 heimbach 1.53 Caja CALL CYCLE_OBCS( k, bi, bj, myThid )
824     CALL SET_OBCS( k, bi, bj, myTime+deltaTclock, myThid )
825 adcroft 1.47 ENDIF
826     #endif
827 adcroft 1.11 C-- Prediction step (step forward all model variables)
828     CALL TIMESTEP(
829 heimbach 1.53 I bi,bj,iMin,iMax,jMin,jMax,k,
830 adcroft 1.46 I myIter, myThid)
831 adcroft 1.43 #ifdef ALLOW_OBCS
832 adcroft 1.41 C-- Apply open boundary conditions
833 heimbach 1.49 IF (openBoundaries) THEN
834     #ifdef ALLOW_AUTODIFF_TAMC
835 heimbach 1.55 CADJ STORE gunm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
836     CADJ STORE gvnm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
837     CADJ STORE gwnm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
838 heimbach 1.53 #endif /* ALLOW_AUTODIFF_TAMC */
839 heimbach 1.54
840 heimbach 1.53 CALL APPLY_OBCS2( bi, bj, k, myThid )
841 heimbach 1.49 END IF
842 adcroft 1.43 #endif
843 adcroft 1.41 C-- Freeze water
844 heimbach 1.49 IF (allowFreezing) THEN
845     #ifdef ALLOW_AUTODIFF_TAMC
846 heimbach 1.55 CADJ STORE gTNm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
847 heimbach 1.53 #endif /* ALLOW_AUTODIFF_TAMC */
848     CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
849 heimbach 1.49 END IF
850 adcroft 1.48
851     #ifdef DIVG_IN_DYNAMICS
852 adcroft 1.11 C-- Diagnose barotropic divergence of predicted fields
853 cnh 1.31 CALL CALC_DIV_GHAT(
854 heimbach 1.53 I bi,bj,iMin,iMax,jMin,jMax,k,
855 adcroft 1.11 I xA,yA,
856     I myThid)
857 adcroft 1.48 #endif /* DIVG_IN_DYNAMICS */
858 adcroft 1.23
859     C-- Cumulative diagnostic calculations (ie. time-averaging)
860 cnh 1.38 #ifdef INCLUDE_DIAGNOSTICS_INTERFACE_CODE
861 adcroft 1.23 IF (taveFreq.GT.0.) THEN
862     CALL DO_TIME_AVERAGES(
863 heimbach 1.53 I myTime, myIter, bi, bj, k, kup, kDown,
864 adcroft 1.50 I rVel, ConvectCount,
865 adcroft 1.23 I myThid )
866     ENDIF
867     #endif
868 adcroft 1.45
869 adcroft 1.11
870 heimbach 1.53 C-- k loop
871     ENDDO
872 adcroft 1.12
873 heimbach 1.49 #ifdef ALLOW_AUTODIFF_TAMC
874     maximpl = 6
875 heimbach 1.51 iikey = (ikey-1)*maximpl
876 heimbach 1.53 #endif /* ALLOW_AUTODIFF_TAMC */
877 heimbach 1.51
878     C-- Implicit diffusion
879     IF (implicitDiffusion) THEN
880 heimbach 1.49
881     IF (tempStepping) THEN
882     #ifdef ALLOW_AUTODIFF_TAMC
883     idkey = iikey + 1
884 heimbach 1.55 CADJ STORE gTNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
885 heimbach 1.53 #endif /* ALLOW_AUTODIFF_TAMC */
886 heimbach 1.49 CALL IMPLDIFF(
887 adcroft 1.42 I bi, bj, iMin, iMax, jMin, jMax,
888     I deltaTtracer, KappaRT,recip_HFacC,
889     U gTNm1,
890     I myThid )
891 heimbach 1.49 END IF
892    
893     IF (saltStepping) THEN
894     #ifdef ALLOW_AUTODIFF_TAMC
895     idkey = iikey + 2
896 heimbach 1.55 CADJ STORE gSNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
897 heimbach 1.53 #endif /* ALLOW_AUTODIFF_TAMC */
898 heimbach 1.49 CALL IMPLDIFF(
899 adcroft 1.42 I bi, bj, iMin, iMax, jMin, jMax,
900     I deltaTtracer, KappaRS,recip_HFacC,
901     U gSNm1,
902     I myThid )
903 heimbach 1.49 END IF
904    
905 heimbach 1.53 C-- implicitDiffusion
906     ENDIF
907 heimbach 1.49
908 adcroft 1.44 C-- Implicit viscosity
909     IF (implicitViscosity) THEN
910 heimbach 1.49
911 adcroft 1.42 IF (momStepping) THEN
912 heimbach 1.49 #ifdef ALLOW_AUTODIFF_TAMC
913     idkey = iikey + 3
914 heimbach 1.55 CADJ STORE gUNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
915 heimbach 1.53 #endif /* ALLOW_AUTODIFF_TAMC */
916 adcroft 1.42 CALL IMPLDIFF(
917     I bi, bj, iMin, iMax, jMin, jMax,
918     I deltaTmom, KappaRU,recip_HFacW,
919     U gUNm1,
920     I myThid )
921 heimbach 1.49 #ifdef ALLOW_AUTODIFF_TAMC
922     idkey = iikey + 4
923 heimbach 1.55 CADJ STORE gVNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
924 heimbach 1.53 #endif /* ALLOW_AUTODIFF_TAMC */
925 adcroft 1.42 CALL IMPLDIFF(
926     I bi, bj, iMin, iMax, jMin, jMax,
927     I deltaTmom, KappaRV,recip_HFacS,
928     U gVNm1,
929     I myThid )
930 heimbach 1.49
931 adcroft 1.42 #ifdef INCLUDE_CD_CODE
932 heimbach 1.49
933     #ifdef ALLOW_AUTODIFF_TAMC
934     idkey = iikey + 5
935 heimbach 1.55 CADJ STORE vVelD(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
936 heimbach 1.53 #endif /* ALLOW_AUTODIFF_TAMC */
937 adcroft 1.42 CALL IMPLDIFF(
938     I bi, bj, iMin, iMax, jMin, jMax,
939     I deltaTmom, KappaRU,recip_HFacW,
940     U vVelD,
941     I myThid )
942 heimbach 1.49 #ifdef ALLOW_AUTODIFF_TAMC
943     idkey = iikey + 6
944 heimbach 1.55 CADJ STORE uVelD(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
945 heimbach 1.53 #endif /* ALLOW_AUTODIFF_TAMC */
946 adcroft 1.42 CALL IMPLDIFF(
947     I bi, bj, iMin, iMax, jMin, jMax,
948     I deltaTmom, KappaRV,recip_HFacS,
949     U uVelD,
950     I myThid )
951 heimbach 1.49
952 adcroft 1.42 #endif
953 heimbach 1.49
954 heimbach 1.53 C-- momStepping
955     ENDIF
956    
957     C-- implicitViscosity
958     ENDIF
959 cnh 1.1
960     ENDDO
961     ENDDO
962    
963     RETURN
964     END

  ViewVC Help
Powered by ViewVC 1.1.22