/[MITgcm]/MITgcm_contrib/ecco_utils/ecco_v4_release4_devel/code/dynamics.F
ViewVC logotype

Annotation of /MITgcm_contrib/ecco_utils/ecco_v4_release4_devel/code/dynamics.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.1 - (hide annotations) (download)
Mon Sep 30 20:28:05 2019 UTC (5 years, 9 months ago) by ou.wang
Branch: MAIN
Create v4r3 repository

1 ou.wang 1.1 C $Header: /u/gcmpack/MITgcm/model/src/dynamics.F,v 1.178 2016/11/28 23:05:05 jmc Exp $
2     C $Name: $
3    
4     #include "PACKAGES_CONFIG.h"
5     #include "CPP_OPTIONS.h"
6     #ifdef ALLOW_AUTODIFF
7     # include "AUTODIFF_OPTIONS.h"
8     #endif
9     #ifdef ALLOW_MOM_COMMON
10     # include "MOM_COMMON_OPTIONS.h"
11     #endif
12     #ifdef ALLOW_OBCS
13     # include "OBCS_OPTIONS.h"
14     #endif
15    
16     #undef DYNAMICS_GUGV_EXCH_CHECK
17    
18     CBOP
19     C !ROUTINE: DYNAMICS
20     C !INTERFACE:
21     SUBROUTINE DYNAMICS(myTime, myIter, myThid)
22     C !DESCRIPTION: \bv
23     C *==========================================================*
24     C | SUBROUTINE DYNAMICS
25     C | o Controlling routine for the explicit part of the model
26     C | dynamics.
27     C *==========================================================*
28     C \ev
29     C !USES:
30     IMPLICIT NONE
31     C == Global variables ===
32     #include "SIZE.h"
33     #include "EEPARAMS.h"
34     #include "PARAMS.h"
35     #include "GRID.h"
36     #include "DYNVARS.h"
37     #ifdef ALLOW_MOM_COMMON
38     # include "MOM_VISC.h"
39     #endif
40     #ifdef ALLOW_CD_CODE
41     # include "CD_CODE_VARS.h"
42     #endif
43     #ifdef ALLOW_AUTODIFF
44     # include "tamc.h"
45     # include "tamc_keys.h"
46     # include "FFIELDS.h"
47     # include "EOS.h"
48     # ifdef ALLOW_KPP
49     # include "KPP.h"
50     # endif
51     # ifdef ALLOW_PTRACERS
52     # include "PTRACERS_SIZE.h"
53     # include "PTRACERS_FIELDS.h"
54     # endif
55     # ifdef ALLOW_OBCS
56     # include "OBCS_PARAMS.h"
57     # include "OBCS_FIELDS.h"
58     # ifdef ALLOW_PTRACERS
59     # include "OBCS_PTRACERS.h"
60     # endif
61     # endif
62     # ifdef ALLOW_MOM_FLUXFORM
63     # include "MOM_FLUXFORM.h"
64     # endif
65     #endif /* ALLOW_AUTODIFF */
66     #ifdef ALLOW_ECCO
67     # include "ecco.h"
68     #endif
69    
70     C !CALLING SEQUENCE:
71     C DYNAMICS()
72     C |
73     C |-- CALC_EP_FORCING
74     C |
75     C |-- CALC_GRAD_PHI_SURF
76     C |
77     C |-- CALC_VISCOSITY
78     C |
79     C |-- MOM_CALC_3D_STRAIN
80     C |
81     C |-- CALC_EDDY_STRESS
82     C |
83     C |-- CALC_PHI_HYD
84     C |
85     C |-- MOM_FLUXFORM
86     C |
87     C |-- MOM_VECINV
88     C |
89     C |-- MOM_CALC_SMAG_3D
90     C |-- MOM_UV_SMAG_3D
91     C |
92     C |-- TIMESTEP
93     C |
94     C |-- MOM_U_IMPLICIT_R
95     C |-- MOM_V_IMPLICIT_R
96     C |
97     C |-- IMPLDIFF
98     C |
99     C |-- OBCS_APPLY_UV
100     C |
101     C |-- CALC_GW
102     C |
103     C |-- DIAGNOSTICS_FILL
104     C |-- DEBUG_STATS_RL
105    
106     C !INPUT/OUTPUT PARAMETERS:
107     C == Routine arguments ==
108     C myTime :: Current time in simulation
109     C myIter :: Current iteration number in simulation
110     C myThid :: Thread number for this instance of the routine.
111     _RL myTime
112     INTEGER myIter
113     INTEGER myThid
114    
115     C !FUNCTIONS:
116     #ifdef ALLOW_DIAGNOSTICS
117     LOGICAL DIAGNOSTICS_IS_ON
118     EXTERNAL DIAGNOSTICS_IS_ON
119     #endif
120    
121     C !LOCAL VARIABLES:
122     C == Local variables
123     C fVer[UV] o fVer: Vertical flux term - note fVer
124     C is "pipelined" in the vertical
125     C so we need an fVer for each
126     C variable.
127     C phiHydC :: hydrostatic potential anomaly at cell center
128     C In z coords phiHyd is the hydrostatic potential
129     C (=pressure/rho0) anomaly
130     C In p coords phiHyd is the geopotential height anomaly.
131     C phiHydF :: hydrostatic potential anomaly at middle between 2 centers
132     C dPhiHydX,Y :: Gradient (X & Y directions) of hydrostatic potential anom.
133     C phiSurfX, :: gradient of Surface potential (Pressure/rho, ocean)
134     C phiSurfY or geopotential (atmos) in X and Y direction
135     C guDissip :: dissipation tendency (all explicit terms), u component
136     C gvDissip :: dissipation tendency (all explicit terms), v component
137     C kappaRU :: vertical viscosity for velocity U-component
138     C kappaRV :: vertical viscosity for velocity V-component
139     C iMin, iMax :: Ranges and sub-block indices on which calculations
140     C jMin, jMax are applied.
141     C bi, bj :: tile indices
142     C k :: current level index
143     C km1, kp1 :: index of level above (k-1) and below (k+1)
144     C kUp, kDown :: Index for interface above and below. kUp and kDown are
145     C are switched with k to be the appropriate index into fVerU,V
146     _RL fVerU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
147     _RL fVerV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
148     _RL phiHydF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
149     _RL phiHydC (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
150     _RL dPhiHydX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
151     _RL dPhiHydY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
152     _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
153     _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
154     _RL guDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
155     _RL gvDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
156     _RL kappaRU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
157     _RL kappaRV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
158     #ifdef ALLOW_SMAG_3D
159     C str11 :: strain component Vxx @ grid-cell center
160     C str22 :: strain component Vyy @ grid-cell center
161     C str33 :: strain component Vzz @ grid-cell center
162     C str12 :: strain component Vxy @ grid-cell corner
163     C str13 :: strain component Vxz @ above uVel
164     C str23 :: strain component Vyz @ above vVel
165     C viscAh3d_00 :: Smagorinsky viscosity @ grid-cell center
166     C viscAh3d_12 :: Smagorinsky viscosity @ grid-cell corner
167     C viscAh3d_13 :: Smagorinsky viscosity @ above uVel
168     C viscAh3d_23 :: Smagorinsky viscosity @ above vVel
169     C addDissU :: zonal momentum tendency from 3-D Smag. viscosity
170     C addDissV :: merid momentum tendency from 3-D Smag. viscosity
171     _RL str11(1-OLx:sNx+OLx,1-OLy:sNy+OLy, Nr )
172     _RL str22(1-OLx:sNx+OLx,1-OLy:sNy+OLy, Nr )
173     _RL str33(1-OLx:sNx+OLx,1-OLy:sNy+OLy, Nr )
174     _RL str12(1-OLx:sNx+OLx,1-OLy:sNy+OLy, Nr )
175     _RL str13(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
176     _RL str23(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
177     _RL viscAh3d_00(1-OLx:sNx+OLx,1-OLy:sNy+OLy, Nr )
178     _RL viscAh3d_12(1-OLx:sNx+OLx,1-OLy:sNy+OLy, Nr )
179     _RL viscAh3d_13(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
180     _RL viscAh3d_23(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
181     _RL addDissU(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
182     _RL addDissV(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
183     #elif ( defined ALLOW_NONHYDROSTATIC )
184     _RL str13(1), str23(1), str33(1)
185     _RL viscAh3d_00(1), viscAh3d_13(1), viscAh3d_23(1)
186     #endif
187    
188     INTEGER bi, bj
189     INTEGER i, j
190     INTEGER k, km1, kp1, kUp, kDown
191     INTEGER iMin, iMax
192     INTEGER jMin, jMax
193     PARAMETER( iMin = 0 , iMax = sNx+1 )
194     PARAMETER( jMin = 0 , jMax = sNy+1 )
195    
196     #ifdef ALLOW_DIAGNOSTICS
197     LOGICAL dPhiHydDiagIsOn
198     _RL tmpFac
199     #endif /* ALLOW_DIAGNOSTICS */
200    
201     C--- The algorithm...
202     C
203     C "Correction Step"
204     C =================
205     C Here we update the horizontal velocities with the surface
206     C pressure such that the resulting flow is either consistent
207     C with the free-surface evolution or the rigid-lid:
208     C U[n] = U* + dt x d/dx P
209     C V[n] = V* + dt x d/dy P
210     C
211     C "Calculation of Gs"
212     C ===================
213     C This is where all the accelerations and tendencies (ie.
214     C physics, parameterizations etc...) are calculated
215     C rho = rho ( theta[n], salt[n] )
216     C b = b(rho, theta)
217     C K31 = K31 ( rho )
218     C Gu[n] = Gu( u[n], v[n], wVel, b, ... )
219     C Gv[n] = Gv( u[n], v[n], wVel, b, ... )
220     C Gt[n] = Gt( theta[n], u[n], v[n], wVel, K31, ... )
221     C Gs[n] = Gs( salt[n], u[n], v[n], wVel, K31, ... )
222     C
223     C "Time-stepping" or "Prediction"
224     C ================================
225     C The models variables are stepped forward with the appropriate
226     C time-stepping scheme (currently we use Adams-Bashforth II)
227     C - For momentum, the result is always *only* a "prediction"
228     C in that the flow may be divergent and will be "corrected"
229     C later with a surface pressure gradient.
230     C - Normally for tracers the result is the new field at time
231     C level [n+1} *BUT* in the case of implicit diffusion the result
232     C is also *only* a prediction.
233     C - We denote "predictors" with an asterisk (*).
234     C U* = U[n] + dt x ( 3/2 Gu[n] - 1/2 Gu[n-1] )
235     C V* = V[n] + dt x ( 3/2 Gv[n] - 1/2 Gv[n-1] )
236     C theta[n+1] = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
237     C salt[n+1] = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
238     C With implicit diffusion:
239     C theta* = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
240     C salt* = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
241     C (1 + dt * K * d_zz) theta[n] = theta*
242     C (1 + dt * K * d_zz) salt[n] = salt*
243     C---
244     CEOP
245    
246     #ifdef ALLOW_DEBUG
247     IF (debugMode) CALL DEBUG_ENTER( 'DYNAMICS', myThid )
248     #endif
249    
250     #ifdef ALLOW_DIAGNOSTICS
251     dPhiHydDiagIsOn = .FALSE.
252     IF ( useDiagnostics )
253     & dPhiHydDiagIsOn = DIAGNOSTICS_IS_ON( 'Um_dPHdx', myThid )
254     & .OR. DIAGNOSTICS_IS_ON( 'Vm_dPHdy', myThid )
255     #endif
256    
257     C-- Call to routine for calculation of Eliassen-Palm-flux-forced
258     C U-tendency, if desired:
259     #ifdef INCLUDE_EP_FORCING_CODE
260     CALL CALC_EP_FORCING(myThid)
261     #endif
262    
263     #ifdef ALLOW_AUTODIFF_MONITOR_DIAG
264     CALL DUMMY_IN_DYNAMICS( myTime, myIter, myThid )
265     #endif
266    
267     #ifdef ALLOW_AUTODIFF_TAMC
268     C-- HPF directive to help TAMC
269     CHPF$ INDEPENDENT
270     #endif /* ALLOW_AUTODIFF_TAMC */
271    
272     DO bj=myByLo(myThid),myByHi(myThid)
273    
274     #ifdef ALLOW_AUTODIFF_TAMC
275     C-- HPF directive to help TAMC
276     CHPF$ INDEPENDENT, NEW (fVerU,fVerV
277     CHPF$& ,phiHydF
278     CHPF$& ,kappaRU,kappaRV
279     CHPF$& )
280     #endif /* ALLOW_AUTODIFF_TAMC */
281    
282     DO bi=myBxLo(myThid),myBxHi(myThid)
283    
284     #ifdef ALLOW_AUTODIFF_TAMC
285     act1 = bi - myBxLo(myThid)
286     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
287     act2 = bj - myByLo(myThid)
288     max2 = myByHi(myThid) - myByLo(myThid) + 1
289     act3 = myThid - 1
290     max3 = nTx*nTy
291     act4 = ikey_dynamics - 1
292     idynkey = (act1 + 1) + act2*max1
293     & + act3*max1*max2
294     & + act4*max1*max2*max3
295     #endif /* ALLOW_AUTODIFF_TAMC */
296    
297     C-- Set up work arrays with valid (i.e. not NaN) values
298     C These initial values do not alter the numerical results. They
299     C just ensure that all memory references are to valid floating
300     C point numbers. This prevents spurious hardware signals due to
301     C uninitialised but inert locations.
302    
303     #ifdef ALLOW_AUTODIFF
304     DO k=1,Nr
305     DO j=1-OLy,sNy+OLy
306     DO i=1-OLx,sNx+OLx
307     c-- need some re-initialisation here to break dependencies
308     gU(i,j,k,bi,bj) = 0. _d 0
309     gV(i,j,k,bi,bj) = 0. _d 0
310     ENDDO
311     ENDDO
312     ENDDO
313     #endif /* ALLOW_AUTODIFF */
314     DO j=1-OLy,sNy+OLy
315     DO i=1-OLx,sNx+OLx
316     fVerU (i,j,1) = 0. _d 0
317     fVerU (i,j,2) = 0. _d 0
318     fVerV (i,j,1) = 0. _d 0
319     fVerV (i,j,2) = 0. _d 0
320     phiHydF (i,j) = 0. _d 0
321     phiHydC (i,j) = 0. _d 0
322     #ifndef INCLUDE_PHIHYD_CALCULATION_CODE
323     dPhiHydX(i,j) = 0. _d 0
324     dPhiHydY(i,j) = 0. _d 0
325     #endif
326     phiSurfX(i,j) = 0. _d 0
327     phiSurfY(i,j) = 0. _d 0
328     guDissip(i,j) = 0. _d 0
329     gvDissip(i,j) = 0. _d 0
330     #ifdef ALLOW_AUTODIFF
331     phiHydLow(i,j,bi,bj) = 0. _d 0
332     # if (defined NONLIN_FRSURF) && (defined ALLOW_MOM_FLUXFORM)
333     # ifndef DISABLE_RSTAR_CODE
334     dWtransC(i,j,bi,bj) = 0. _d 0
335     dWtransU(i,j,bi,bj) = 0. _d 0
336     dWtransV(i,j,bi,bj) = 0. _d 0
337     # endif
338     # endif
339     #endif /* ALLOW_AUTODIFF */
340     ENDDO
341     ENDDO
342    
343     C-- Start computation of dynamics
344    
345     #ifdef ALLOW_AUTODIFF_TAMC
346     CADJ STORE wVel (:,:,:,bi,bj) =
347     CADJ & comlev1_bibj, key=idynkey, byte=isbyte
348     #endif /* ALLOW_AUTODIFF_TAMC */
349    
350     C-- Explicit part of the Surface Potential Gradient (add in TIMESTEP)
351     C (note: this loop will be replaced by CALL CALC_GRAD_ETA)
352     IF (implicSurfPress.NE.1.) THEN
353     CALL CALC_GRAD_PHI_SURF(
354     I bi,bj,iMin,iMax,jMin,jMax,
355     I etaN,
356     O phiSurfX,phiSurfY,
357     I myThid )
358     ENDIF
359    
360     #ifdef ALLOW_AUTODIFF_TAMC
361     CADJ STORE uVel (:,:,:,bi,bj) = comlev1_bibj, key=idynkey, byte=isbyte
362     CADJ STORE vVel (:,:,:,bi,bj) = comlev1_bibj, key=idynkey, byte=isbyte
363     #ifdef ALLOW_KPP
364     CADJ STORE KPPviscAz (:,:,:,bi,bj)
365     CADJ & = comlev1_bibj, key=idynkey, byte=isbyte
366     #endif /* ALLOW_KPP */
367     #endif /* ALLOW_AUTODIFF_TAMC */
368    
369     #ifndef ALLOW_AUTODIFF
370     IF ( .NOT.momViscosity ) THEN
371     #endif
372     DO k=1,Nr+1
373     DO j=1-OLy,sNy+OLy
374     DO i=1-OLx,sNx+OLx
375     kappaRU(i,j,k) = 0. _d 0
376     kappaRV(i,j,k) = 0. _d 0
377     ENDDO
378     ENDDO
379     ENDDO
380     #ifndef ALLOW_AUTODIFF
381     ENDIF
382     #endif
383     #ifdef INCLUDE_CALC_DIFFUSIVITY_CALL
384     C-- Calculate the total vertical viscosity
385     IF ( momViscosity ) THEN
386     CALL CALC_VISCOSITY(
387     I bi,bj, iMin,iMax,jMin,jMax,
388     O kappaRU, kappaRV,
389     I myThid )
390     ENDIF
391     #endif /* INCLUDE_CALC_DIFFUSIVITY_CALL */
392    
393     #ifdef ALLOW_SMAG_3D
394     IF ( useSmag3D ) THEN
395     CALL MOM_CALC_3D_STRAIN(
396     O str11, str22, str33, str12, str13, str23,
397     I bi, bj, myThid )
398     ENDIF
399     #endif /* ALLOW_SMAG_3D */
400    
401     #ifdef ALLOW_AUTODIFF_TAMC
402     CADJ STORE kappaRU(:,:,:)
403     CADJ & = comlev1_bibj, key=idynkey, byte=isbyte
404     CADJ STORE kappaRV(:,:,:)
405     CADJ & = comlev1_bibj, key=idynkey, byte=isbyte
406     #endif /* ALLOW_AUTODIFF_TAMC */
407    
408     #ifdef ALLOW_OBCS
409     C-- For Stevens boundary conditions velocities need to be extrapolated
410     C (copied) to a narrow strip outside the domain
411     IF ( useOBCS ) THEN
412     CALL OBCS_COPY_UV_N(
413     U uVel(1-OLx,1-OLy,1,bi,bj),
414     U vVel(1-OLx,1-OLy,1,bi,bj),
415     I Nr, bi, bj, myThid )
416     ENDIF
417     #endif /* ALLOW_OBCS */
418    
419     #ifdef ALLOW_EDDYPSI
420     CALL CALC_EDDY_STRESS(bi,bj,myThid)
421     #endif
422    
423     C-- Start of dynamics loop
424     DO k=1,Nr
425    
426     C-- km1 Points to level above k (=k-1)
427     C-- kup Cycles through 1,2 to point to layer above
428     C-- kDown Cycles through 2,1 to point to current layer
429    
430     km1 = MAX(1,k-1)
431     kp1 = MIN(k+1,Nr)
432     kup = 1+MOD(k+1,2)
433     kDown= 1+MOD(k,2)
434    
435     #ifdef ALLOW_AUTODIFF_TAMC
436     kkey = (idynkey-1)*Nr + k
437     CADJ STORE totPhiHyd (:,:,k,bi,bj)
438     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
439     CADJ STORE phiHydLow (:,:,bi,bj)
440     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
441     CADJ STORE theta (:,:,k,bi,bj)
442     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
443     CADJ STORE salt (:,:,k,bi,bj)
444     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
445     # ifdef NONLIN_FRSURF
446     cph-test
447     CADJ STORE phiHydC (:,:)
448     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
449     CADJ STORE phiHydF (:,:)
450     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
451     CADJ STORE gU(:,:,k,bi,bj)
452     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
453     CADJ STORE gV(:,:,k,bi,bj)
454     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
455     # ifndef ALLOW_ADAMSBASHFORTH_3
456     CADJ STORE guNm1(:,:,k,bi,bj)
457     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
458     CADJ STORE gvNm1(:,:,k,bi,bj)
459     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
460     # else
461     CADJ STORE guNm(:,:,k,bi,bj,1)
462     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
463     CADJ STORE guNm(:,:,k,bi,bj,2)
464     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
465     CADJ STORE gvNm(:,:,k,bi,bj,1)
466     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
467     CADJ STORE gvNm(:,:,k,bi,bj,2)
468     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
469     # endif
470     # ifdef ALLOW_CD_CODE
471     CADJ STORE uNM1(:,:,k,bi,bj)
472     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
473     CADJ STORE vNM1(:,:,k,bi,bj)
474     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
475     CADJ STORE uVelD(:,:,k,bi,bj)
476     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
477     CADJ STORE vVelD(:,:,k,bi,bj)
478     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
479     # endif
480     # endif /* NONLIN_FRSURF */
481     #endif /* ALLOW_AUTODIFF_TAMC */
482    
483     C-- Integrate hydrostatic balance for phiHyd with BC of phiHyd(z=0)=0
484     CALL CALC_PHI_HYD(
485     I bi,bj,iMin,iMax,jMin,jMax,k,
486     I theta, salt,
487     U phiHydF,
488     O phiHydC, dPhiHydX, dPhiHydY,
489     I myTime, myIter, myThid )
490     #ifdef ALLOW_DIAGNOSTICS
491     IF ( dPhiHydDiagIsOn ) THEN
492     tmpFac = -1. _d 0
493     CALL DIAGNOSTICS_SCALE_FILL( dPhiHydX, tmpFac, 1,
494     & 'Um_dPHdx', k, 1, 2, bi, bj, myThid )
495     CALL DIAGNOSTICS_SCALE_FILL( dPhiHydY, tmpFac, 1,
496     & 'Vm_dPHdy', k, 1, 2, bi, bj, myThid )
497     ENDIF
498     #endif /* ALLOW_DIAGNOSTICS */
499    
500     C-- Calculate accelerations in the momentum equations (gU, gV, ...)
501     C and step forward storing the result in gU, gV, etc...
502     IF ( momStepping ) THEN
503     #ifdef ALLOW_AUTODIFF
504     DO j=1-OLy,sNy+OLy
505     DO i=1-OLx,sNx+OLx
506     guDissip(i,j) = 0. _d 0
507     gvDissip(i,j) = 0. _d 0
508     ENDDO
509     ENDDO
510     #endif /* ALLOW_AUTODIFF */
511     #ifdef ALLOW_AUTODIFF_TAMC
512     # if (defined NONLIN_FRSURF) && (defined ALLOW_MOM_FLUXFORM)
513     # ifndef DISABLE_RSTAR_CODE
514     CADJ STORE dWtransC(:,:,bi,bj)
515     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
516     CADJ STORE dWtransU(:,:,bi,bj)
517     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
518     CADJ STORE dWtransV(:,:,bi,bj)
519     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
520     # endif
521     # endif /* NONLIN_FRSURF and ALLOW_MOM_FLUXFORM */
522     # if (defined NONLIN_FRSURF) || (defined ALLOW_DEPTH_CONTROL)
523     CADJ STORE fVerU(:,:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
524     CADJ STORE fVerV(:,:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
525     # endif
526     #endif /* ALLOW_AUTODIFF_TAMC */
527     IF (.NOT. vectorInvariantMomentum) THEN
528     #ifdef ALLOW_MOM_FLUXFORM
529     CALL MOM_FLUXFORM(
530     I bi,bj,k,iMin,iMax,jMin,jMax,
531     I kappaRU, kappaRV,
532     U fVerU(1-OLx,1-OLy,kUp), fVerV(1-OLx,1-OLy,kUp),
533     O fVerU(1-OLx,1-OLy,kDown), fVerV(1-OLx,1-OLy,kDown),
534     O guDissip, gvDissip,
535     I myTime, myIter, myThid)
536     #endif
537     ELSE
538     #ifdef ALLOW_MOM_VECINV
539     CALL MOM_VECINV(
540     I bi,bj,k,iMin,iMax,jMin,jMax,
541     I kappaRU, kappaRV,
542     I fVerU(1-OLx,1-OLy,kUp), fVerV(1-OLx,1-OLy,kUp),
543     O fVerU(1-OLx,1-OLy,kDown), fVerV(1-OLx,1-OLy,kDown),
544     O guDissip, gvDissip,
545     I myTime, myIter, myThid)
546     #endif
547     ENDIF
548    
549     #ifdef ALLOW_SMAG_3D
550     IF ( useSmag3D ) THEN
551     CALL MOM_CALC_SMAG_3D(
552     I str11, str22, str33, str12, str13, str23,
553     O viscAh3d_00, viscAh3d_12, viscAh3d_13, viscAh3d_23,
554     I smag3D_hLsC, smag3D_hLsW, smag3D_hLsS, smag3D_hLsZ,
555     I k, bi, bj, myThid )
556     CALL MOM_UV_SMAG_3D(
557     I str11, str22, str12, str13, str23,
558     I viscAh3d_00, viscAh3d_12, viscAh3d_13, viscAh3d_23,
559     O addDissU, addDissV,
560     I iMin,iMax,jMin,jMax, k, bi, bj, myThid )
561     DO j= jMin,jMax
562     DO i= iMin,iMax
563     guDissip(i,j) = guDissip(i,j) + addDissU(i,j)
564     gvDissip(i,j) = gvDissip(i,j) + addDissV(i,j)
565     ENDDO
566     ENDDO
567     ENDIF
568     #endif /* ALLOW_SMAG_3D */
569    
570     CALL TIMESTEP(
571     I bi,bj,iMin,iMax,jMin,jMax,k,
572     I dPhiHydX,dPhiHydY, phiSurfX, phiSurfY,
573     I guDissip, gvDissip,
574     I myTime, myIter, myThid)
575    
576     ENDIF
577    
578     C-- end of dynamics k loop (1:Nr)
579     ENDDO
580    
581     C-- Implicit Vertical advection & viscosity
582     #if (defined (INCLUDE_IMPLVERTADV_CODE) && \
583     defined (ALLOW_MOM_COMMON) && !(defined ALLOW_AUTODIFF))
584     IF ( momImplVertAdv .OR. implicitViscosity
585     & .OR. selectImplicitDrag.GE.1 ) THEN
586     C to recover older (prior to 2016-10-05) results:
587     c IF ( momImplVertAdv ) THEN
588     CALL MOM_U_IMPLICIT_R( kappaRU,
589     I bi, bj, myTime, myIter, myThid )
590     CALL MOM_V_IMPLICIT_R( kappaRV,
591     I bi, bj, myTime, myIter, myThid )
592     ELSEIF ( implicitViscosity ) THEN
593     #else /* INCLUDE_IMPLVERTADV_CODE */
594     IF ( implicitViscosity ) THEN
595     #endif /* INCLUDE_IMPLVERTADV_CODE */
596     #ifdef ALLOW_AUTODIFF_TAMC
597     CADJ STORE gU(:,:,:,bi,bj) = comlev1_bibj , key=idynkey, byte=isbyte
598     #endif /* ALLOW_AUTODIFF_TAMC */
599     CALL IMPLDIFF(
600     I bi, bj, iMin, iMax, jMin, jMax,
601     I -1, kappaRU, recip_hFacW(1-OLx,1-OLy,1,bi,bj),
602     U gU(1-OLx,1-OLy,1,bi,bj),
603     I myThid )
604     #ifdef ALLOW_AUTODIFF_TAMC
605     CADJ STORE gV(:,:,:,bi,bj) = comlev1_bibj , key=idynkey, byte=isbyte
606     #endif /* ALLOW_AUTODIFF_TAMC */
607     CALL IMPLDIFF(
608     I bi, bj, iMin, iMax, jMin, jMax,
609     I -2, kappaRV, recip_hFacS(1-OLx,1-OLy,1,bi,bj),
610     U gV(1-OLx,1-OLy,1,bi,bj),
611     I myThid )
612     ENDIF
613    
614     #ifdef ALLOW_OBCS
615     C-- Apply open boundary conditions
616     IF ( useOBCS ) THEN
617     C-- but first save intermediate velocities to be used in the
618     C next time step for the Stevens boundary conditions
619     CALL OBCS_SAVE_UV_N(
620     I bi, bj, iMin, iMax, jMin, jMax, 0,
621     I gU, gV, myThid )
622     CALL OBCS_APPLY_UV( bi, bj, 0, gU, gV, myThid )
623     ENDIF
624     #endif /* ALLOW_OBCS */
625    
626     #ifdef ALLOW_CD_CODE
627     IF (implicitViscosity.AND.useCDscheme) THEN
628     #ifdef ALLOW_AUTODIFF_TAMC
629     CADJ STORE vVelD(:,:,:,bi,bj) = comlev1_bibj , key=idynkey, byte=isbyte
630     #endif /* ALLOW_AUTODIFF_TAMC */
631     CALL IMPLDIFF(
632     I bi, bj, iMin, iMax, jMin, jMax,
633     I 0, kappaRU, recip_hFacW(1-OLx,1-OLy,1,bi,bj),
634     U vVelD(1-OLx,1-OLy,1,bi,bj),
635     I myThid )
636     #ifdef ALLOW_AUTODIFF_TAMC
637     CADJ STORE uVelD(:,:,:,bi,bj) = comlev1_bibj , key=idynkey, byte=isbyte
638     #endif /* ALLOW_AUTODIFF_TAMC */
639     CALL IMPLDIFF(
640     I bi, bj, iMin, iMax, jMin, jMax,
641     I 0, kappaRV, recip_hFacS(1-OLx,1-OLy,1,bi,bj),
642     U uVelD(1-OLx,1-OLy,1,bi,bj),
643     I myThid )
644     ENDIF
645     #endif /* ALLOW_CD_CODE */
646     C-- End implicit Vertical advection & viscosity
647    
648     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
649    
650     #ifdef ALLOW_NONHYDROSTATIC
651     C-- Step forward W field in N-H algorithm
652     IF ( nonHydrostatic ) THEN
653     #ifdef ALLOW_DEBUG
654     IF (debugMode) CALL DEBUG_CALL('CALC_GW', myThid )
655     #endif
656     CALL TIMER_START('CALC_GW [DYNAMICS]',myThid)
657     CALL CALC_GW(
658     I bi,bj, kappaRU, kappaRV,
659     I str13, str23, str33,
660     I viscAh3d_00, viscAh3d_13, viscAh3d_23,
661     I myTime, myIter, myThid )
662     ENDIF
663     IF ( nonHydrostatic.OR.implicitIntGravWave )
664     & CALL TIMESTEP_WVEL( bi,bj, myTime, myIter, myThid )
665     IF ( nonHydrostatic )
666     & CALL TIMER_STOP ('CALC_GW [DYNAMICS]',myThid)
667     #endif
668    
669     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
670    
671     C- end of bi,bj loops
672     ENDDO
673     ENDDO
674    
675     #ifdef ALLOW_OBCS
676     IF (useOBCS) THEN
677     CALL OBCS_EXCHANGES( myThid )
678     ENDIF
679     #endif
680    
681     Cml(
682     C In order to compare the variance of phiHydLow of a p/z-coordinate
683     C run with etaH of a z/p-coordinate run the drift of phiHydLow
684     C has to be removed by something like the following subroutine:
685     C CALL REMOVE_MEAN_RL( 1, phiHydLow, maskInC, maskInC, rA, drF,
686     C & 'phiHydLow', myTime, myThid )
687     Cml)
688    
689     #ifdef ALLOW_DIAGNOSTICS
690     IF ( useDiagnostics ) THEN
691    
692     CALL DIAGNOSTICS_FILL(totPhihyd,'PHIHYD ',0,Nr,0,1,1,myThid)
693     CALL DIAGNOSTICS_FILL(phiHydLow,'PHIBOT ',0, 1,0,1,1,myThid)
694     CALL DIAGNOSTICS_FILL(m_bp,'OBP ',0, 1,0,1,1,myThid)
695     CALL DIAGNOSTICS_FILL(m_bp_nopabar,'OBPNOPAB',0, 1,0,1,1,myThid)
696    
697     tmpFac = 1. _d 0
698     CALL DIAGNOSTICS_SCALE_FILL(totPhihyd,tmpFac,2,
699     & 'PHIHYDSQ',0,Nr,0,1,1,myThid)
700    
701     CALL DIAGNOSTICS_SCALE_FILL(phiHydLow,tmpFac,2,
702     & 'PHIBOTSQ',0, 1,0,1,1,myThid)
703    
704     ENDIF
705     #endif /* ALLOW_DIAGNOSTICS */
706    
707     #ifdef ALLOW_DEBUG
708     IF ( debugLevel .GE. debLevD ) THEN
709     CALL DEBUG_STATS_RL(1,EtaN,'EtaN (DYNAMICS)',myThid)
710     CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (DYNAMICS)',myThid)
711     CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (DYNAMICS)',myThid)
712     CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (DYNAMICS)',myThid)
713     CALL DEBUG_STATS_RL(Nr,theta,'Theta (DYNAMICS)',myThid)
714     CALL DEBUG_STATS_RL(Nr,salt,'Salt (DYNAMICS)',myThid)
715     CALL DEBUG_STATS_RL(Nr,gU,'Gu (DYNAMICS)',myThid)
716     CALL DEBUG_STATS_RL(Nr,gV,'Gv (DYNAMICS)',myThid)
717     #ifndef ALLOW_ADAMSBASHFORTH_3
718     CALL DEBUG_STATS_RL(Nr,guNm1,'GuNm1 (DYNAMICS)',myThid)
719     CALL DEBUG_STATS_RL(Nr,gvNm1,'GvNm1 (DYNAMICS)',myThid)
720     CALL DEBUG_STATS_RL(Nr,gtNm1,'GtNm1 (DYNAMICS)',myThid)
721     CALL DEBUG_STATS_RL(Nr,gsNm1,'GsNm1 (DYNAMICS)',myThid)
722     #endif
723     ENDIF
724     #endif
725    
726     #ifdef DYNAMICS_GUGV_EXCH_CHECK
727     C- jmc: For safety checking only: This Exchange here should not change
728     C the solution. If solution changes, it means something is wrong,
729     C but it does not mean that it is less wrong with this exchange.
730     IF ( debugLevel .GE. debLevE ) THEN
731     CALL EXCH_UV_XYZ_RL(gU,gV,.TRUE.,myThid)
732     ENDIF
733     #endif
734    
735     #ifdef ALLOW_DEBUG
736     IF (debugMode) CALL DEBUG_LEAVE( 'DYNAMICS', myThid )
737     #endif
738    
739     RETURN
740     END

  ViewVC Help
Powered by ViewVC 1.1.22