/[MITgcm]/MITgcm_contrib/ksnow/press_release/code_expt/dynamics.F
ViewVC logotype

Annotation of /MITgcm_contrib/ksnow/press_release/code_expt/dynamics.F

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


Revision 1.1 - (hide annotations) (download)
Sun Feb 12 13:32:38 2017 UTC (8 years, 5 months ago) by dgoldberg
Branch: MAIN
CVS Tags: HEAD
overlap update for phi0surf required. fix to TOPDR to avoid zero division

1 dgoldberg 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     #ifdef ALLOW_SHELFICE
16     # include "SHELFICE_OPTIONS.h"
17     #endif
18    
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_SHELFICE_GROUNDED_ICE
70     # include "SURFACE.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     #ifdef ALLOW_SHELFICE_GROUNDED_ICE
286     _EXCH_XY_RL (phi0surf, mythid)
287     #endif
288    
289     DO bi=myBxLo(myThid),myBxHi(myThid)
290    
291     #ifdef ALLOW_AUTODIFF_TAMC
292     act1 = bi - myBxLo(myThid)
293     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
294     act2 = bj - myByLo(myThid)
295     max2 = myByHi(myThid) - myByLo(myThid) + 1
296     act3 = myThid - 1
297     max3 = nTx*nTy
298     act4 = ikey_dynamics - 1
299     idynkey = (act1 + 1) + act2*max1
300     & + act3*max1*max2
301     & + act4*max1*max2*max3
302     #endif /* ALLOW_AUTODIFF_TAMC */
303    
304     C-- Set up work arrays with valid (i.e. not NaN) values
305     C These initial values do not alter the numerical results. They
306     C just ensure that all memory references are to valid floating
307     C point numbers. This prevents spurious hardware signals due to
308     C uninitialised but inert locations.
309    
310     #ifdef ALLOW_AUTODIFF
311     DO k=1,Nr
312     DO j=1-OLy,sNy+OLy
313     DO i=1-OLx,sNx+OLx
314     c-- need some re-initialisation here to break dependencies
315     gU(i,j,k,bi,bj) = 0. _d 0
316     gV(i,j,k,bi,bj) = 0. _d 0
317     ENDDO
318     ENDDO
319     ENDDO
320     #endif /* ALLOW_AUTODIFF */
321     DO j=1-OLy,sNy+OLy
322     DO i=1-OLx,sNx+OLx
323     fVerU (i,j,1) = 0. _d 0
324     fVerU (i,j,2) = 0. _d 0
325     fVerV (i,j,1) = 0. _d 0
326     fVerV (i,j,2) = 0. _d 0
327     phiHydF (i,j) = 0. _d 0
328     phiHydC (i,j) = 0. _d 0
329     #ifndef INCLUDE_PHIHYD_CALCULATION_CODE
330     dPhiHydX(i,j) = 0. _d 0
331     dPhiHydY(i,j) = 0. _d 0
332     #endif
333     phiSurfX(i,j) = 0. _d 0
334     phiSurfY(i,j) = 0. _d 0
335     guDissip(i,j) = 0. _d 0
336     gvDissip(i,j) = 0. _d 0
337     #ifdef ALLOW_AUTODIFF
338     phiHydLow(i,j,bi,bj) = 0. _d 0
339     # if (defined NONLIN_FRSURF) && (defined ALLOW_MOM_FLUXFORM)
340     # ifndef DISABLE_RSTAR_CODE
341     dWtransC(i,j,bi,bj) = 0. _d 0
342     dWtransU(i,j,bi,bj) = 0. _d 0
343     dWtransV(i,j,bi,bj) = 0. _d 0
344     # endif
345     # endif
346     #endif /* ALLOW_AUTODIFF */
347     ENDDO
348     ENDDO
349    
350     C-- Start computation of dynamics
351    
352     #ifdef ALLOW_AUTODIFF_TAMC
353     CADJ STORE wVel (:,:,:,bi,bj) =
354     CADJ & comlev1_bibj, key=idynkey, byte=isbyte
355     #endif /* ALLOW_AUTODIFF_TAMC */
356    
357     C-- Explicit part of the Surface Potential Gradient (add in TIMESTEP)
358     C (note: this loop will be replaced by CALL CALC_GRAD_ETA)
359     IF (implicSurfPress.NE.1.) THEN
360     CALL CALC_GRAD_PHI_SURF(
361     I bi,bj,iMin,iMax,jMin,jMax,
362     I etaN,
363     O phiSurfX,phiSurfY,
364     I myThid )
365     ENDIF
366    
367     #ifdef ALLOW_AUTODIFF_TAMC
368     CADJ STORE uVel (:,:,:,bi,bj) = comlev1_bibj, key=idynkey, byte=isbyte
369     CADJ STORE vVel (:,:,:,bi,bj) = comlev1_bibj, key=idynkey, byte=isbyte
370     #ifdef ALLOW_KPP
371     CADJ STORE KPPviscAz (:,:,:,bi,bj)
372     CADJ & = comlev1_bibj, key=idynkey, byte=isbyte
373     #endif /* ALLOW_KPP */
374     #endif /* ALLOW_AUTODIFF_TAMC */
375    
376     #ifndef ALLOW_AUTODIFF
377     IF ( .NOT.momViscosity ) THEN
378     #endif
379     DO k=1,Nr+1
380     DO j=1-OLy,sNy+OLy
381     DO i=1-OLx,sNx+OLx
382     kappaRU(i,j,k) = 0. _d 0
383     kappaRV(i,j,k) = 0. _d 0
384     ENDDO
385     ENDDO
386     ENDDO
387     #ifndef ALLOW_AUTODIFF
388     ENDIF
389     #endif
390     #ifdef INCLUDE_CALC_DIFFUSIVITY_CALL
391     C-- Calculate the total vertical viscosity
392     IF ( momViscosity ) THEN
393     CALL CALC_VISCOSITY(
394     I bi,bj, iMin,iMax,jMin,jMax,
395     O kappaRU, kappaRV,
396     I myThid )
397     ENDIF
398     #endif /* INCLUDE_CALC_DIFFUSIVITY_CALL */
399    
400     #ifdef ALLOW_SMAG_3D
401     IF ( useSmag3D ) THEN
402     CALL MOM_CALC_3D_STRAIN(
403     O str11, str22, str33, str12, str13, str23,
404     I bi, bj, myThid )
405     ENDIF
406     #endif /* ALLOW_SMAG_3D */
407    
408     #ifdef ALLOW_AUTODIFF_TAMC
409     CADJ STORE kappaRU(:,:,:)
410     CADJ & = comlev1_bibj, key=idynkey, byte=isbyte
411     CADJ STORE kappaRV(:,:,:)
412     CADJ & = comlev1_bibj, key=idynkey, byte=isbyte
413     #endif /* ALLOW_AUTODIFF_TAMC */
414    
415     #ifdef ALLOW_OBCS
416     C-- For Stevens boundary conditions velocities need to be extrapolated
417     C (copied) to a narrow strip outside the domain
418     IF ( useOBCS ) THEN
419     CALL OBCS_COPY_UV_N(
420     U uVel(1-OLx,1-OLy,1,bi,bj),
421     U vVel(1-OLx,1-OLy,1,bi,bj),
422     I Nr, bi, bj, myThid )
423     ENDIF
424     #endif /* ALLOW_OBCS */
425    
426     #ifdef ALLOW_EDDYPSI
427     CALL CALC_EDDY_STRESS(bi,bj,myThid)
428     #endif
429    
430     C-- Start of dynamics loop
431     DO k=1,Nr
432    
433     C-- km1 Points to level above k (=k-1)
434     C-- kup Cycles through 1,2 to point to layer above
435     C-- kDown Cycles through 2,1 to point to current layer
436    
437     km1 = MAX(1,k-1)
438     kp1 = MIN(k+1,Nr)
439     kup = 1+MOD(k+1,2)
440     kDown= 1+MOD(k,2)
441    
442     #ifdef ALLOW_AUTODIFF_TAMC
443     kkey = (idynkey-1)*Nr + k
444     CADJ STORE totPhiHyd (:,:,k,bi,bj)
445     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
446     CADJ STORE phiHydLow (:,:,bi,bj)
447     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
448     CADJ STORE theta (:,:,k,bi,bj)
449     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
450     CADJ STORE salt (:,:,k,bi,bj)
451     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
452     # ifdef NONLIN_FRSURF
453     cph-test
454     CADJ STORE phiHydC (:,:)
455     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
456     CADJ STORE phiHydF (:,:)
457     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
458     CADJ STORE gU(:,:,k,bi,bj)
459     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
460     CADJ STORE gV(:,:,k,bi,bj)
461     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
462     # ifndef ALLOW_ADAMSBASHFORTH_3
463     CADJ STORE guNm1(:,:,k,bi,bj)
464     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
465     CADJ STORE gvNm1(:,:,k,bi,bj)
466     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
467     # else
468     CADJ STORE guNm(:,:,k,bi,bj,1)
469     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
470     CADJ STORE guNm(:,:,k,bi,bj,2)
471     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
472     CADJ STORE gvNm(:,:,k,bi,bj,1)
473     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
474     CADJ STORE gvNm(:,:,k,bi,bj,2)
475     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
476     # endif
477     # ifdef ALLOW_CD_CODE
478     CADJ STORE uNM1(:,:,k,bi,bj)
479     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
480     CADJ STORE vNM1(:,:,k,bi,bj)
481     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
482     CADJ STORE uVelD(:,:,k,bi,bj)
483     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
484     CADJ STORE vVelD(:,:,k,bi,bj)
485     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
486     # endif
487     # endif /* NONLIN_FRSURF */
488     #endif /* ALLOW_AUTODIFF_TAMC */
489    
490     C-- Integrate hydrostatic balance for phiHyd with BC of phiHyd(z=0)=0
491     CALL CALC_PHI_HYD(
492     I bi,bj,iMin,iMax,jMin,jMax,k,
493     I theta, salt,
494     U phiHydF,
495     O phiHydC, dPhiHydX, dPhiHydY,
496     I myTime, myIter, myThid )
497     #ifdef ALLOW_DIAGNOSTICS
498     IF ( dPhiHydDiagIsOn ) THEN
499     tmpFac = -1. _d 0
500     CALL DIAGNOSTICS_SCALE_FILL( dPhiHydX, tmpFac, 1,
501     & 'Um_dPHdx', k, 1, 2, bi, bj, myThid )
502     CALL DIAGNOSTICS_SCALE_FILL( dPhiHydY, tmpFac, 1,
503     & 'Vm_dPHdy', k, 1, 2, bi, bj, myThid )
504     ENDIF
505     #endif /* ALLOW_DIAGNOSTICS */
506    
507     C-- Calculate accelerations in the momentum equations (gU, gV, ...)
508     C and step forward storing the result in gU, gV, etc...
509     IF ( momStepping ) THEN
510     #ifdef ALLOW_AUTODIFF
511     DO j=1-OLy,sNy+OLy
512     DO i=1-OLx,sNx+OLx
513     guDissip(i,j) = 0. _d 0
514     gvDissip(i,j) = 0. _d 0
515     ENDDO
516     ENDDO
517     #endif /* ALLOW_AUTODIFF */
518     #ifdef ALLOW_AUTODIFF_TAMC
519     # if (defined NONLIN_FRSURF) && (defined ALLOW_MOM_FLUXFORM)
520     # ifndef DISABLE_RSTAR_CODE
521     CADJ STORE dWtransC(:,:,bi,bj)
522     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
523     CADJ STORE dWtransU(:,:,bi,bj)
524     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
525     CADJ STORE dWtransV(:,:,bi,bj)
526     CADJ & = comlev1_bibj_k, key=kkey, byte=isbyte
527     # endif
528     # endif /* NONLIN_FRSURF and ALLOW_MOM_FLUXFORM */
529     # if (defined NONLIN_FRSURF) || (defined ALLOW_DEPTH_CONTROL)
530     CADJ STORE fVerU(:,:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
531     CADJ STORE fVerV(:,:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
532     # endif
533     #endif /* ALLOW_AUTODIFF_TAMC */
534     IF (.NOT. vectorInvariantMomentum) THEN
535     #ifdef ALLOW_MOM_FLUXFORM
536     CALL MOM_FLUXFORM(
537     I bi,bj,k,iMin,iMax,jMin,jMax,
538     I kappaRU, kappaRV,
539     U fVerU(1-OLx,1-OLy,kUp), fVerV(1-OLx,1-OLy,kUp),
540     O fVerU(1-OLx,1-OLy,kDown), fVerV(1-OLx,1-OLy,kDown),
541     O guDissip, gvDissip,
542     I myTime, myIter, myThid)
543     #endif
544     ELSE
545     #ifdef ALLOW_MOM_VECINV
546     CALL MOM_VECINV(
547     I bi,bj,k,iMin,iMax,jMin,jMax,
548     I kappaRU, kappaRV,
549     I fVerU(1-OLx,1-OLy,kUp), fVerV(1-OLx,1-OLy,kUp),
550     O fVerU(1-OLx,1-OLy,kDown), fVerV(1-OLx,1-OLy,kDown),
551     O guDissip, gvDissip,
552     I myTime, myIter, myThid)
553     #endif
554     ENDIF
555    
556     #ifdef ALLOW_SMAG_3D
557     IF ( useSmag3D ) THEN
558     CALL MOM_CALC_SMAG_3D(
559     I str11, str22, str33, str12, str13, str23,
560     O viscAh3d_00, viscAh3d_12, viscAh3d_13, viscAh3d_23,
561     I smag3D_hLsC, smag3D_hLsW, smag3D_hLsS, smag3D_hLsZ,
562     I k, bi, bj, myThid )
563     CALL MOM_UV_SMAG_3D(
564     I str11, str22, str12, str13, str23,
565     I viscAh3d_00, viscAh3d_12, viscAh3d_13, viscAh3d_23,
566     O addDissU, addDissV,
567     I iMin,iMax,jMin,jMax, k, bi, bj, myThid )
568     DO j= jMin,jMax
569     DO i= iMin,iMax
570     guDissip(i,j) = guDissip(i,j) + addDissU(i,j)
571     gvDissip(i,j) = gvDissip(i,j) + addDissV(i,j)
572     ENDDO
573     ENDDO
574     ENDIF
575     #endif /* ALLOW_SMAG_3D */
576    
577     CALL TIMESTEP(
578     I bi,bj,iMin,iMax,jMin,jMax,k,
579     I dPhiHydX,dPhiHydY, phiSurfX, phiSurfY,
580     I guDissip, gvDissip,
581     I myTime, myIter, myThid)
582    
583     ENDIF
584    
585     C-- end of dynamics k loop (1:Nr)
586     ENDDO
587    
588     C-- Implicit Vertical advection & viscosity
589     #if (defined (INCLUDE_IMPLVERTADV_CODE) && \
590     defined (ALLOW_MOM_COMMON) && !(defined ALLOW_AUTODIFF))
591     IF ( momImplVertAdv .OR. implicitViscosity
592     & .OR. selectImplicitDrag.GE.1 ) THEN
593     C to recover older (prior to 2016-10-05) results:
594     c IF ( momImplVertAdv ) THEN
595     CALL MOM_U_IMPLICIT_R( kappaRU,
596     I bi, bj, myTime, myIter, myThid )
597     CALL MOM_V_IMPLICIT_R( kappaRV,
598     I bi, bj, myTime, myIter, myThid )
599     ELSEIF ( implicitViscosity ) THEN
600     #else /* INCLUDE_IMPLVERTADV_CODE */
601     IF ( implicitViscosity ) THEN
602     #endif /* INCLUDE_IMPLVERTADV_CODE */
603     #ifdef ALLOW_AUTODIFF_TAMC
604     CADJ STORE gU(:,:,:,bi,bj) = comlev1_bibj , key=idynkey, byte=isbyte
605     #endif /* ALLOW_AUTODIFF_TAMC */
606     CALL IMPLDIFF(
607     I bi, bj, iMin, iMax, jMin, jMax,
608     I -1, kappaRU, recip_hFacW(1-OLx,1-OLy,1,bi,bj),
609     U gU(1-OLx,1-OLy,1,bi,bj),
610     I myThid )
611     #ifdef ALLOW_AUTODIFF_TAMC
612     CADJ STORE gV(:,:,:,bi,bj) = comlev1_bibj , key=idynkey, byte=isbyte
613     #endif /* ALLOW_AUTODIFF_TAMC */
614     CALL IMPLDIFF(
615     I bi, bj, iMin, iMax, jMin, jMax,
616     I -2, kappaRV, recip_hFacS(1-OLx,1-OLy,1,bi,bj),
617     U gV(1-OLx,1-OLy,1,bi,bj),
618     I myThid )
619     ENDIF
620    
621     #ifdef ALLOW_OBCS
622     C-- Apply open boundary conditions
623     IF ( useOBCS ) THEN
624     C-- but first save intermediate velocities to be used in the
625     C next time step for the Stevens boundary conditions
626     CALL OBCS_SAVE_UV_N(
627     I bi, bj, iMin, iMax, jMin, jMax, 0,
628     I gU, gV, myThid )
629     CALL OBCS_APPLY_UV( bi, bj, 0, gU, gV, myThid )
630     ENDIF
631     #endif /* ALLOW_OBCS */
632    
633     #ifdef ALLOW_CD_CODE
634     IF (implicitViscosity.AND.useCDscheme) THEN
635     #ifdef ALLOW_AUTODIFF_TAMC
636     CADJ STORE vVelD(:,:,:,bi,bj) = comlev1_bibj , key=idynkey, byte=isbyte
637     #endif /* ALLOW_AUTODIFF_TAMC */
638     CALL IMPLDIFF(
639     I bi, bj, iMin, iMax, jMin, jMax,
640     I 0, kappaRU, recip_hFacW(1-OLx,1-OLy,1,bi,bj),
641     U vVelD(1-OLx,1-OLy,1,bi,bj),
642     I myThid )
643     #ifdef ALLOW_AUTODIFF_TAMC
644     CADJ STORE uVelD(:,:,:,bi,bj) = comlev1_bibj , key=idynkey, byte=isbyte
645     #endif /* ALLOW_AUTODIFF_TAMC */
646     CALL IMPLDIFF(
647     I bi, bj, iMin, iMax, jMin, jMax,
648     I 0, kappaRV, recip_hFacS(1-OLx,1-OLy,1,bi,bj),
649     U uVelD(1-OLx,1-OLy,1,bi,bj),
650     I myThid )
651     ENDIF
652     #endif /* ALLOW_CD_CODE */
653     C-- End implicit Vertical advection & viscosity
654    
655     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
656    
657     #ifdef ALLOW_NONHYDROSTATIC
658     C-- Step forward W field in N-H algorithm
659     IF ( nonHydrostatic ) THEN
660     #ifdef ALLOW_DEBUG
661     IF (debugMode) CALL DEBUG_CALL('CALC_GW', myThid )
662     #endif
663     CALL TIMER_START('CALC_GW [DYNAMICS]',myThid)
664     CALL CALC_GW(
665     I bi,bj, kappaRU, kappaRV,
666     I str13, str23, str33,
667     I viscAh3d_00, viscAh3d_13, viscAh3d_23,
668     I myTime, myIter, myThid )
669     ENDIF
670     IF ( nonHydrostatic.OR.implicitIntGravWave )
671     & CALL TIMESTEP_WVEL( bi,bj, myTime, myIter, myThid )
672     IF ( nonHydrostatic )
673     & CALL TIMER_STOP ('CALC_GW [DYNAMICS]',myThid)
674     #endif
675    
676     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
677    
678     C- end of bi,bj loops
679     ENDDO
680     ENDDO
681    
682     #ifdef ALLOW_OBCS
683     IF (useOBCS) THEN
684     CALL OBCS_EXCHANGES( myThid )
685     ENDIF
686     #endif
687    
688     Cml(
689     C In order to compare the variance of phiHydLow of a p/z-coordinate
690     C run with etaH of a z/p-coordinate run the drift of phiHydLow
691     C has to be removed by something like the following subroutine:
692     C CALL REMOVE_MEAN_RL( 1, phiHydLow, maskInC, maskInC, rA, drF,
693     C & 'phiHydLow', myTime, myThid )
694     Cml)
695    
696     #ifdef ALLOW_DIAGNOSTICS
697     IF ( useDiagnostics ) THEN
698    
699     CALL DIAGNOSTICS_FILL(totPhihyd,'PHIHYD ',0,Nr,0,1,1,myThid)
700     CALL DIAGNOSTICS_FILL(phiHydLow,'PHIBOT ',0, 1,0,1,1,myThid)
701    
702     tmpFac = 1. _d 0
703     CALL DIAGNOSTICS_SCALE_FILL(totPhihyd,tmpFac,2,
704     & 'PHIHYDSQ',0,Nr,0,1,1,myThid)
705    
706     CALL DIAGNOSTICS_SCALE_FILL(phiHydLow,tmpFac,2,
707     & 'PHIBOTSQ',0, 1,0,1,1,myThid)
708    
709     ENDIF
710     #endif /* ALLOW_DIAGNOSTICS */
711    
712     #ifdef ALLOW_DEBUG
713     IF ( debugLevel .GE. debLevD ) THEN
714     CALL DEBUG_STATS_RL(1,EtaN,'EtaN (DYNAMICS)',myThid)
715     CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (DYNAMICS)',myThid)
716     CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (DYNAMICS)',myThid)
717     CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (DYNAMICS)',myThid)
718     CALL DEBUG_STATS_RL(Nr,theta,'Theta (DYNAMICS)',myThid)
719     CALL DEBUG_STATS_RL(Nr,salt,'Salt (DYNAMICS)',myThid)
720     CALL DEBUG_STATS_RL(Nr,gU,'Gu (DYNAMICS)',myThid)
721     CALL DEBUG_STATS_RL(Nr,gV,'Gv (DYNAMICS)',myThid)
722     #ifndef ALLOW_ADAMSBASHFORTH_3
723     CALL DEBUG_STATS_RL(Nr,guNm1,'GuNm1 (DYNAMICS)',myThid)
724     CALL DEBUG_STATS_RL(Nr,gvNm1,'GvNm1 (DYNAMICS)',myThid)
725     CALL DEBUG_STATS_RL(Nr,gtNm1,'GtNm1 (DYNAMICS)',myThid)
726     CALL DEBUG_STATS_RL(Nr,gsNm1,'GsNm1 (DYNAMICS)',myThid)
727     #endif
728     ENDIF
729     #endif
730    
731     #ifdef DYNAMICS_GUGV_EXCH_CHECK
732     C- jmc: For safety checking only: This Exchange here should not change
733     C the solution. If solution changes, it means something is wrong,
734     C but it does not mean that it is less wrong with this exchange.
735     IF ( debugLevel .GE. debLevE ) THEN
736     CALL EXCH_UV_XYZ_RL(gU,gV,.TRUE.,myThid)
737     ENDIF
738     #endif
739    
740     #ifdef ALLOW_DEBUG
741     IF (debugMode) CALL DEBUG_LEAVE( 'DYNAMICS', myThid )
742     #endif
743    
744     RETURN
745     END

  ViewVC Help
Powered by ViewVC 1.1.22