/[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.2 - (hide annotations) (download)
Thu Nov 7 22:53:33 2019 UTC (5 years, 8 months ago) by ou.wang
Branch: MAIN
Changes since 1.1: +5 -2 lines
Update code and namelist related to SSH and OBP

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

  ViewVC Help
Powered by ViewVC 1.1.22