/[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.54 - (hide annotations) (download)
Mon Nov 13 16:32:57 2000 UTC (23 years, 6 months ago) by heimbach
Branch: MAIN
CVS Tags: branch-atmos-merge-start, checkpoint33, checkpoint32
Branch point for: branch-atmos-merge
Changes since 1.53: +27 -13 lines
Rescaling of forcing fields done immediately after reading fields.

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

  ViewVC Help
Powered by ViewVC 1.1.22