/[MITgcm]/MITgcm/model/src/dynamics.F
ViewVC logotype

Diff of /MITgcm/model/src/dynamics.F

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

revision 1.66 by heimbach, Sun Mar 25 22:33:52 2001 UTC revision 1.164 by jahn, Thu Mar 21 18:15:44 2013 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4    #include "PACKAGES_CONFIG.h"
5  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
6    #ifdef ALLOW_OBCS
7    # include "OBCS_OPTIONS.h"
8    #endif
9    
10    #undef DYNAMICS_GUGV_EXCH_CHECK
11    
12    CBOP
13    C     !ROUTINE: DYNAMICS
14    C     !INTERFACE:
15        SUBROUTINE DYNAMICS(myTime, myIter, myThid)        SUBROUTINE DYNAMICS(myTime, myIter, myThid)
16  C     /==========================================================\  C     !DESCRIPTION: \bv
17  C     | SUBROUTINE DYNAMICS                                      |  C     *==========================================================*
18  C     | o Controlling routine for the explicit part of the model |  C     | SUBROUTINE DYNAMICS
19  C     |   dynamics.                                              |  C     | o Controlling routine for the explicit part of the model
20  C     |==========================================================|  C     |   dynamics.
21  C     | This routine evaluates the "dynamics" terms for each     |  C     *==========================================================*
22  C     | block of ocean in turn. Because the blocks of ocean have |  C     | This routine evaluates the "dynamics" terms for each
23  C     | overlap regions they are independent of one another.     |  C     | block of ocean in turn. Because the blocks of ocean have
24  C     | If terms involving lateral integrals are needed in this  |  C     | overlap regions they are independent of one another.
25  C     | routine care will be needed. Similarly finite-difference |  C     | If terms involving lateral integrals are needed in this
26  C     | operations with stencils wider than the overlap region   |  C     | routine care will be needed. Similarly finite-difference
27  C     | require special consideration.                           |  C     | operations with stencils wider than the overlap region
28  C     | Notes                                                    |  C     | require special consideration.
29  C     | =====                                                    |  C     | The algorithm...
30  C     | C*P* comments indicating place holders for which code is |  C     |
31  C     |      presently being developed.                          |  C     | "Correction Step"
32  C     \==========================================================/  C     | =================
33    C     | Here we update the horizontal velocities with the surface
34    C     | pressure such that the resulting flow is either consistent
35    C     | with the free-surface evolution or the rigid-lid:
36    C     |   U[n] = U* + dt x d/dx P
37    C     |   V[n] = V* + dt x d/dy P
38    C     |   W[n] = W* + dt x d/dz P  (NH mode)
39    C     |
40    C     | "Calculation of Gs"
41    C     | ===================
42    C     | This is where all the accelerations and tendencies (ie.
43    C     | physics, parameterizations etc...) are calculated
44    C     |   rho = rho ( theta[n], salt[n] )
45    C     |   b   = b(rho, theta)
46    C     |   K31 = K31 ( rho )
47    C     |   Gu[n] = Gu( u[n], v[n], wVel, b, ... )
48    C     |   Gv[n] = Gv( u[n], v[n], wVel, b, ... )
49    C     |   Gt[n] = Gt( theta[n], u[n], v[n], wVel, K31, ... )
50    C     |   Gs[n] = Gs( salt[n], u[n], v[n], wVel, K31, ... )
51    C     |
52    C     | "Time-stepping" or "Prediction"
53    C     | ================================
54    C     | The models variables are stepped forward with the appropriate
55    C     | time-stepping scheme (currently we use Adams-Bashforth II)
56    C     | - For momentum, the result is always *only* a "prediction"
57    C     | in that the flow may be divergent and will be "corrected"
58    C     | later with a surface pressure gradient.
59    C     | - Normally for tracers the result is the new field at time
60    C     | level [n+1} *BUT* in the case of implicit diffusion the result
61    C     | is also *only* a prediction.
62    C     | - We denote "predictors" with an asterisk (*).
63    C     |   U* = U[n] + dt x ( 3/2 Gu[n] - 1/2 Gu[n-1] )
64    C     |   V* = V[n] + dt x ( 3/2 Gv[n] - 1/2 Gv[n-1] )
65    C     |   theta[n+1] = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
66    C     |   salt[n+1] = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
67    C     | With implicit diffusion:
68    C     |   theta* = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
69    C     |   salt* = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
70    C     |   (1 + dt * K * d_zz) theta[n] = theta*
71    C     |   (1 + dt * K * d_zz) salt[n] = salt*
72    C     |
73    C     *==========================================================*
74    C     \ev
75    C     !USES:
76        IMPLICIT NONE        IMPLICIT NONE
   
77  C     == Global variables ===  C     == Global variables ===
78  #include "SIZE.h"  #include "SIZE.h"
79  #include "EEPARAMS.h"  #include "EEPARAMS.h"
80  #include "PARAMS.h"  #include "PARAMS.h"
81  #include "DYNVARS.h"  #include "DYNVARS.h"
82    #ifdef ALLOW_CD_CODE
83    #include "CD_CODE_VARS.h"
84    #endif
85  #include "GRID.h"  #include "GRID.h"
   
86  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
87  # include "tamc.h"  # include "tamc.h"
88  # include "tamc_keys.h"  # include "tamc_keys.h"
89  #endif /* ALLOW_AUTODIFF_TAMC */  # include "FFIELDS.h"
90    # include "EOS.h"
91  #ifdef ALLOW_KPP  # ifdef ALLOW_KPP
92  # include "KPP.h"  #  include "KPP.h"
93  #endif  # endif
94    # ifdef ALLOW_PTRACERS
95  #ifdef ALLOW_TIMEAVE  #  include "PTRACERS_SIZE.h"
96  #include "TIMEAVE_STATV.h"  #  include "PTRACERS_FIELDS.h"
97  #endif  # endif
98    # ifdef ALLOW_OBCS
99    #include "OBCS_PARAMS.h"
100    #  include "OBCS_FIELDS.h"
101    #  ifdef ALLOW_PTRACERS
102    #   include "OBCS_PTRACERS.h"
103    #  endif
104    # endif
105    # ifdef ALLOW_MOM_FLUXFORM
106    #  include "MOM_FLUXFORM.h"
107    # endif
108    #endif /* ALLOW_AUTODIFF_TAMC */
109    
110    C     !CALLING SEQUENCE:
111    C     DYNAMICS()
112    C      |
113    C      |-- CALC_EP_FORCING
114    C      |
115    C      |-- CALC_GRAD_PHI_SURF
116    C      |
117    C      |-- CALC_VISCOSITY
118    C      |
119    C      |-- CALC_PHI_HYD
120    C      |
121    C      |-- MOM_FLUXFORM
122    C      |
123    C      |-- MOM_VECINV
124    C      |
125    C      |-- TIMESTEP
126    C      |
127    C      |-- MOM_U_IMPLICIT_R
128    C      |-- MOM_V_IMPLICIT_R
129    C      |
130    C      |-- IMPLDIFF
131    C      |
132    C      |-- OBCS_APPLY_UV
133    C      |
134    C      |-- CALC_GW
135    C      |
136    C      |-- DIAGNOSTICS_FILL
137    C      |-- DEBUG_STATS_RL
138    
139    C     !INPUT/OUTPUT PARAMETERS:
140  C     == Routine arguments ==  C     == Routine arguments ==
141  C     myTime - Current time in simulation  C     myTime :: Current time in simulation
142  C     myIter - Current iteration number in simulation  C     myIter :: Current iteration number in simulation
143  C     myThid - Thread number for this instance of the routine.  C     myThid :: Thread number for this instance of the routine.
144        _RL myTime        _RL myTime
145        INTEGER myIter        INTEGER myIter
146        INTEGER myThid        INTEGER myThid
147    
148    C     !FUNCTIONS:
149    #ifdef ALLOW_DIAGNOSTICS
150          LOGICAL  DIAGNOSTICS_IS_ON
151          EXTERNAL DIAGNOSTICS_IS_ON
152    #endif
153    
154    C     !LOCAL VARIABLES:
155  C     == Local variables  C     == Local variables
156  C     xA, yA                 - Per block temporaries holding face areas  C     fVer[UV]               o fVer: Vertical flux term - note fVer
157  C     uTrans, vTrans, rTrans - Per block temporaries holding flow  C                                    is "pipelined" in the vertical
158  C                              transport  C                                    so we need an fVer for each
159  C                              o uTrans: Zonal transport  C                                    variable.
160  C                              o vTrans: Meridional transport  C     phiHydC    :: hydrostatic potential anomaly at cell center
161  C                              o rTrans: Vertical transport  C                   In z coords phiHyd is the hydrostatic potential
162  C     maskC,maskUp             o maskC: land/water mask for tracer cells  C                      (=pressure/rho0) anomaly
163  C                              o maskUp: land/water mask for W points  C                   In p coords phiHyd is the geopotential height anomaly.
164  C     fVer[STUV]               o fVer: Vertical flux term - note fVer  C     phiHydF    :: hydrostatic potential anomaly at middle between 2 centers
165  C                                      is "pipelined" in the vertical  C     dPhiHydX,Y :: Gradient (X & Y directions) of hydrostatic potential anom.
166  C                                      so we need an fVer for each  C     phiSurfX,  ::  gradient of Surface potential (Pressure/rho, ocean)
167  C                                      variable.  C     phiSurfY             or geopotential (atmos) in X and Y direction
168  C     rhoK, rhoKM1   - Density at current level, and level above  C     guDissip   :: dissipation tendency (all explicit terms), u component
169  C     phiHyd         - Hydrostatic part of the potential phiHydi.  C     gvDissip   :: dissipation tendency (all explicit terms), v component
170  C                      In z coords phiHydiHyd is the hydrostatic  C     KappaRU    :: vertical viscosity
171  C                      Potential (=pressure/rho0) anomaly  C     KappaRV    :: vertical viscosity
172  C                      In p coords phiHydiHyd is the geopotential  C     iMin, iMax :: Ranges and sub-block indices on which calculations
173  C                      surface height anomaly.  C     jMin, jMax    are applied.
174  C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)  C     bi, bj     :: tile indices
175  C     phiSurfY             or geopotentiel (atmos) in X and Y direction  C     k          :: current level index
176  C     KappaRT,       - Total diffusion in vertical for T and S.  C     km1, kp1   :: index of level above (k-1) and below (k+1)
177  C     KappaRS          (background + spatially varying, isopycnal term).  C     kUp, kDown :: Index for interface above and below. kUp and kDown are
178  C     iMin, iMax     - Ranges and sub-block indices on which calculations  C                   are switched with k to be the appropriate index into fVerU,V
 C     jMin, jMax       are applied.  
 C     bi, bj  
 C     k, kup,        - Index for layer above and below. kup and kDown  
 C     kDown, km1       are switched with layer to be the appropriate  
 C                      index into fVerTerm.  
       _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RS maskC   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)  
       _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)  
179        _RL fVerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
180        _RL fVerV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
181        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL phiHydF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
182        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiHydC (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
183        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL dPhiHydX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
184          _RL dPhiHydY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
185        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
186        _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
187        _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL guDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
188        _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL gvDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
189        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
190        _RL KappaRV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
       _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
       _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
       _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
   
 C This is currently used by IVDC and Diagnostics  
       _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
191    
192        INTEGER iMin, iMax        INTEGER iMin, iMax
193        INTEGER jMin, jMax        INTEGER jMin, jMax
194        INTEGER bi, bj        INTEGER bi, bj
195        INTEGER i, j        INTEGER i, j
196        INTEGER k, km1, kup, kDown        INTEGER k, km1, kp1, kUp, kDown
197    
198    #ifdef ALLOW_DIAGNOSTICS
199          LOGICAL dPhiHydDiagIsOn
200          _RL tmpFac
201    #endif /* ALLOW_DIAGNOSTICS */
202    
203    
 Cjmc : add for phiHyd output <- but not working if multi tile per CPU  
 c     CHARACTER*(MAX_LEN_MBUF) suff  
 c     LOGICAL  DIFFERENT_MULTIPLE  
 c     EXTERNAL DIFFERENT_MULTIPLE  
 Cjmc(end)  
   
204  C---    The algorithm...  C---    The algorithm...
205  C  C
206  C       "Correction Step"  C       "Correction Step"
# Line 162  C         salt* = salt[n] + dt x ( 3/2 G Line 244  C         salt* = salt[n] + dt x ( 3/2 G
244  C         (1 + dt * K * d_zz) theta[n] = theta*  C         (1 + dt * K * d_zz) theta[n] = theta*
245  C         (1 + dt * K * d_zz) salt[n] = salt*  C         (1 + dt * K * d_zz) salt[n] = salt*
246  C---  C---
247    CEOP
248    
249  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_DEBUG
250  C--   dummy statement to end declaration part        IF (debugMode) CALL DEBUG_ENTER( 'DYNAMICS', myThid )
251        ikey = 1  #endif
 #endif /* ALLOW_AUTODIFF_TAMC */  
252    
253  C--   Set up work arrays with valid (i.e. not NaN) values  #ifdef ALLOW_DIAGNOSTICS
254  C     These inital values do not alter the numerical results. They        dPhiHydDiagIsOn = .FALSE.
255  C     just ensure that all memory references are to valid floating        IF ( useDiagnostics )
256  C     point numbers. This prevents spurious hardware signals due to       &  dPhiHydDiagIsOn = DIAGNOSTICS_IS_ON( 'Um_dPHdx', myThid )
257  C     uninitialised but inert locations.       &               .OR. DIAGNOSTICS_IS_ON( 'Vm_dPHdy', myThid )
258        DO j=1-OLy,sNy+OLy  #endif
        DO i=1-OLx,sNx+OLx  
         xA(i,j)      = 0. _d 0  
         yA(i,j)      = 0. _d 0  
         uTrans(i,j)  = 0. _d 0  
         vTrans(i,j)  = 0. _d 0  
         DO k=1,Nr  
          phiHyd(i,j,k)  = 0. _d 0  
          KappaRU(i,j,k) = 0. _d 0  
          KappaRV(i,j,k) = 0. _d 0  
          sigmaX(i,j,k) = 0. _d 0  
          sigmaY(i,j,k) = 0. _d 0  
          sigmaR(i,j,k) = 0. _d 0  
         ENDDO  
         rhoKM1 (i,j) = 0. _d 0  
         rhok   (i,j) = 0. _d 0  
         maskC  (i,j) = 0. _d 0  
         phiSurfX(i,j) = 0. _d 0  
         phiSurfY(i,j) = 0. _d 0  
        ENDDO  
       ENDDO  
259    
260    C-- Call to routine for calculation of
261    C   Eliassen-Palm-flux-forced U-tendency,
262    C   if desired:
263    #ifdef INCLUDE_EP_FORCING_CODE
264          CALL CALC_EP_FORCING(myThid)
265    #endif
266    
267    #ifdef ALLOW_AUTODIFF_MONITOR_DIAG
268          CALL DUMMY_IN_DYNAMICS( myTime, myIter, myThid )
269    #endif
270    
271  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
272  C--   HPF directive to help TAMC  C--   HPF directive to help TAMC
# Line 205  CHPF$ INDEPENDENT Line 277  CHPF$ INDEPENDENT
277    
278  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
279  C--    HPF directive to help TAMC  C--    HPF directive to help TAMC
280  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS,fVerU,fVerV  CHPF$  INDEPENDENT, NEW (fVerU,fVerV
281  CHPF$&                  ,phiHyd,utrans,vtrans,maskc,xA,yA  CHPF$&                  ,phiHydF
282  CHPF$&                  ,KappaRT,KappaRS,KappaRU,KappaRV  CHPF$&                  ,KappaRU,KappaRV
283  CHPF$&                  )  CHPF$&                  )
284  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
285    
# Line 216  CHPF$&                  ) Line 288  CHPF$&                  )
288  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
289            act1 = bi - myBxLo(myThid)            act1 = bi - myBxLo(myThid)
290            max1 = myBxHi(myThid) - myBxLo(myThid) + 1            max1 = myBxHi(myThid) - myBxLo(myThid) + 1
   
291            act2 = bj - myByLo(myThid)            act2 = bj - myByLo(myThid)
292            max2 = myByHi(myThid) - myByLo(myThid) + 1            max2 = myByHi(myThid) - myByLo(myThid) + 1
   
293            act3 = myThid - 1            act3 = myThid - 1
294            max3 = nTx*nTy            max3 = nTx*nTy
   
295            act4 = ikey_dynamics - 1            act4 = ikey_dynamics - 1
296              idynkey = (act1 + 1) + act2*max1
           ikey = (act1 + 1) + act2*max1  
297       &                      + act3*max1*max2       &                      + act3*max1*max2
298       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
299  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
300    
301  C--     Set up work arrays that need valid initial values  C--   Set up work arrays with valid (i.e. not NaN) values
302          DO j=1-OLy,sNy+OLy  C     These initial values do not alter the numerical results. They
303           DO i=1-OLx,sNx+OLx  C     just ensure that all memory references are to valid floating
304            rTrans(i,j)   = 0. _d 0  C     point numbers. This prevents spurious hardware signals due to
305            fVerT (i,j,1) = 0. _d 0  C     uninitialised but inert locations.
           fVerT (i,j,2) = 0. _d 0  
           fVerS (i,j,1) = 0. _d 0  
           fVerS (i,j,2) = 0. _d 0  
           fVerU (i,j,1) = 0. _d 0  
           fVerU (i,j,2) = 0. _d 0  
           fVerV (i,j,1) = 0. _d 0  
           fVerV (i,j,2) = 0. _d 0  
          ENDDO  
         ENDDO  
306    
307    #ifdef ALLOW_AUTODIFF_TAMC
308          DO k=1,Nr          DO k=1,Nr
309           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
310            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
311  C This is currently also used by IVDC and Diagnostics             KappaRU(i,j,k) = 0. _d 0
312             ConvectCount(i,j,k) = 0.             KappaRV(i,j,k) = 0. _d 0
313             KappaRT(i,j,k) = 0. _d 0  cph(
314             KappaRS(i,j,k) = 0. _d 0  c--   need some re-initialisation here to break dependencies
315    cph)
316               gU(i,j,k,bi,bj) = 0. _d 0
317               gV(i,j,k,bi,bj) = 0. _d 0
318            ENDDO            ENDDO
319           ENDDO           ENDDO
320          ENDDO          ENDDO
   
         iMin = 1-OLx+1  
         iMax = sNx+OLx  
         jMin = 1-OLy+1  
         jMax = sNy+OLy  
   
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 CADJ STORE uvel(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
321  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
322            DO j=1-OLy,sNy+OLy
323  C--     Start of diagnostic loop           DO i=1-OLx,sNx+OLx
324          DO k=Nr,1,-1            fVerU  (i,j,1) = 0. _d 0
325              fVerU  (i,j,2) = 0. _d 0
326  #ifdef ALLOW_AUTODIFF_TAMC            fVerV  (i,j,1) = 0. _d 0
327  C? Patrick, is this formula correct now that we change the loop range?            fVerV  (i,j,2) = 0. _d 0
328  C? Do we still need this?            phiHydF (i,j)  = 0. _d 0
329  cph kkey formula corrected.            phiHydC (i,j)  = 0. _d 0
330  cph Needed for rhok, rhokm1, in the case useGMREDI.  #ifndef INCLUDE_PHIHYD_CALCULATION_CODE
331           kkey = (ikey-1)*Nr + k            dPhiHydX(i,j)  = 0. _d 0
332  CADJ STORE rhokm1(:,:) = comlev1_bibj_k , key = kkey, byte = isbyte            dPhiHydY(i,j)  = 0. _d 0
333  CADJ STORE rhok  (:,:) = comlev1_bibj_k , key = kkey, byte = isbyte  #endif
334  #endif /* ALLOW_AUTODIFF_TAMC */            phiSurfX(i,j)  = 0. _d 0
335              phiSurfY(i,j)  = 0. _d 0
336  C--       Integrate continuity vertically for vertical velocity            guDissip(i,j)  = 0. _d 0
337            CALL INTEGRATE_FOR_W(            gvDissip(i,j)  = 0. _d 0
338       I                         bi, bj, k, uVel, vVel,  #ifdef ALLOW_AUTODIFF_TAMC
339       O                         wVel,            phiHydLow(i,j,bi,bj) = 0. _d 0
340       I                         myThid )  # if (defined NONLIN_FRSURF) && (defined ALLOW_MOM_FLUXFORM)
341    #  ifndef DISABLE_RSTAR_CODE
342  #ifdef    ALLOW_OBCS  #   ifndef ALLOW_AUTODIFF_OPENAD
343  #ifdef    ALLOW_NONHYDROSTATIC            dWtransC(i,j,bi,bj) = 0. _d 0
344  C--       Apply OBC to W if in N-H mode            dWtransU(i,j,bi,bj) = 0. _d 0
345            IF (useOBCS.AND.nonHydrostatic) THEN            dWtransV(i,j,bi,bj) = 0. _d 0
346              CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )  #   endif
347            ENDIF  #  endif
348  #endif    /* ALLOW_NONHYDROSTATIC */  # endif
349  #endif    /* ALLOW_OBCS */  #endif
350             ENDDO
 C--       Calculate gradients of potential density for isoneutral  
 C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)  
 c         IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN  
           IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN  
             CALL FIND_RHO(  
      I        bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,  
      I        theta, salt,  
      O        rhoK,  
      I        myThid )  
             IF (k.GT.1) CALL FIND_RHO(  
      I        bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,  
      I        theta, salt,  
      O        rhoKm1,  
      I        myThid )  
             CALL GRAD_SIGMA(  
      I             bi, bj, iMin, iMax, jMin, jMax, k,  
      I             rhoK, rhoKm1, rhoK,  
      O             sigmaX, sigmaY, sigmaR,  
      I             myThid )  
           ENDIF  
   
 C--       Implicit Vertical Diffusion for Convection  
 c ==> should use sigmaR !!!  
           IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN  
             CALL CALC_IVDC(  
      I        bi, bj, iMin, iMax, jMin, jMax, k,  
      I        rhoKm1, rhoK,  
      U        ConvectCount, KappaRT, KappaRS,  
      I        myTime, myIter, myThid)  
           ENDIF  
   
 C--     end of diagnostic k loop (Nr:1)  
351          ENDDO          ENDDO
352    
353  #ifdef  ALLOW_OBCS  C--     Start computation of dynamics
354  C--     Calculate future values on open boundaries          iMin = 0
355          IF (useOBCS) THEN          iMax = sNx+1
356            CALL OBCS_CALC( bi, bj, myTime+deltaT,          jMin = 0
357       I            uVel, vVel, wVel, theta, salt,          jMax = sNy+1
      I            myThid )  
         ENDIF  
 #endif  /* ALLOW_OBCS */  
358    
 C--     Determines forcing terms based on external fields  
 C       relaxation terms, etc.  
         CALL EXTERNAL_FORCING_SURF(  
      I             bi, bj, iMin, iMax, jMin, jMax,  
      I             myThid )  
   
 #ifdef  ALLOW_GMREDI  
 C--     Calculate iso-neutral slopes for the GM/Redi parameterisation  
         IF (useGMRedi) THEN  
           DO k=1,Nr  
             CALL GMREDI_CALC_TENSOR(  
      I             bi, bj, iMin, iMax, jMin, jMax, k,  
      I             sigmaX, sigmaY, sigmaR,  
      I             myThid )  
           ENDDO  
359  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
360          ELSE  CADJ STORE wVel (:,:,:,bi,bj) =
361            DO k=1, Nr  CADJ &     comlev1_bibj, key=idynkey, byte=isbyte
             CALL GMREDI_CALC_TENSOR_DUMMY(  
      I             bi, bj, iMin, iMax, jMin, jMax, k,  
      I             sigmaX, sigmaY, sigmaR,  
      I             myThid )  
           ENDDO  
362  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
         ENDIF  
 #endif  /* ALLOW_GMREDI */  
363    
364  #ifdef  ALLOW_KPP  C--     Explicit part of the Surface Potential Gradient (add in TIMESTEP)
365  C--     Compute KPP mixing coefficients  C       (note: this loop will be replaced by CALL CALC_GRAD_ETA)
366          IF (useKPP) THEN          IF (implicSurfPress.NE.1.) THEN
367            CALL KPP_CALC(            CALL CALC_GRAD_PHI_SURF(
368       I                  bi, bj, myTime, myThid )       I         bi,bj,iMin,iMax,jMin,jMax,
369  #ifdef ALLOW_AUTODIFF_TAMC       I         etaN,
370          ELSE       O         phiSurfX,phiSurfY,
371            DO j=1-OLy,sNy+OLy       I         myThid )
             DO i=1-OLx,sNx+OLx  
               KPPhbl (i,j,bi,bj) = 1.0  
               KPPfrac(i,j,bi,bj) = 0.0  
               DO k = 1,Nr  
                  KPPghat   (i,j,k,bi,bj) = 0.0  
                  KPPviscAz (i,j,k,bi,bj) = viscAz  
                  KPPdiffKzT(i,j,k,bi,bj) = diffKzT  
                  KPPdiffKzS(i,j,k,bi,bj) = diffKzS  
               ENDDO  
             ENDDO  
           ENDDO  
 #endif /* ALLOW_AUTODIFF_TAMC */  
372          ENDIF          ENDIF
373    
374  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
375  CADJ STORE KPPghat   (:,:,:,bi,bj)  CADJ STORE uVel (:,:,:,bi,bj) = comlev1_bibj, key=idynkey, byte=isbyte
376  CADJ &   , KPPviscAz (:,:,:,bi,bj)  CADJ STORE vVel (:,:,:,bi,bj) = comlev1_bibj, key=idynkey, byte=isbyte
377  CADJ &   , KPPdiffKzT(:,:,:,bi,bj)  #ifdef ALLOW_KPP
378  CADJ &   , KPPdiffKzS(:,:,:,bi,bj)  CADJ STORE KPPviscAz (:,:,:,bi,bj)
379  CADJ &   , KPPfrac   (:,:  ,bi,bj)  CADJ &                 = comlev1_bibj, key=idynkey, byte=isbyte
380  CADJ &                 = comlev1_bibj, key=ikey, byte=isbyte  #endif /* ALLOW_KPP */
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 #endif  /* ALLOW_KPP */  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte  
 CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte  
 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 #ifdef ALLOW_AIM  
 C       AIM - atmospheric intermediate model, physics package code.  
 C note(jmc) : phiHyd=0 at this point but is not really used in Molteni Physics  
         IF ( useAIM ) THEN  
          CALL TIMER_START('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)  
          CALL AIM_DO_ATMOS_PHYSICS( phiHyd, myTime, myThid )  
          CALL TIMER_STOP ('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)  
         ENDIF  
 #endif /* ALLOW_AIM */  
   
   
 C--     Start of thermodynamics loop  
         DO k=Nr,1,-1  
   
 C--       km1    Points to level above k (=k-1)  
 C--       kup    Cycles through 1,2 to point to layer above  
 C--       kDown  Cycles through 2,1 to point to current layer  
   
           km1  = MAX(1,k-1)  
           kup  = 1+MOD(k+1,2)  
           kDown= 1+MOD(k,2)  
   
           iMin = 1-OLx+2  
           iMax = sNx+OLx-1  
           jMin = 1-OLy+2  
           jMax = sNy+OLy-1  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 C? Patrick Is this formula correct?  
 cph Yes, but I rewrote it.  
 cph Also, the KappaR? need the index k!  
          kkey = (ikey-1)*Nr + k  
 CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key = kkey, byte = isbyte  
 CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key = kkey, byte = isbyte  
381  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
382    
 C--      Get temporary terms used by tendency routines  
          CALL CALC_COMMON_FACTORS (  
      I        bi,bj,iMin,iMax,jMin,jMax,k,km1,kup,kDown,  
      O        xA,yA,uTrans,vTrans,rTrans,maskC,maskUp,  
      I        myThid)  
   
383  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL
384  C--      Calculate the total vertical diffusivity  C--     Calculate the total vertical viscosity
385           CALL CALC_DIFFUSIVITY(          CALL CALC_VISCOSITY(
386       I        bi,bj,iMin,iMax,jMin,jMax,k,       I            bi,bj, iMin,iMax,jMin,jMax,
387       I        maskC,maskup,       O            KappaRU, KappaRV,
388       O        KappaRT,KappaRS,KappaRU,KappaRV,       I            myThid )
389       I        myThid)  #else
390  #endif          DO k=1,Nr
391             DO j=1-OLy,sNy+OLy
392  C--      Calculate active tracer tendencies (gT,gS,...)            DO i=1-OLx,sNx+OLx
393  C        and step forward storing result in gTnm1, gSnm1, etc.             KappaRU(i,j,k) = 0. _d 0
394           IF ( tempStepping ) THEN             KappaRV(i,j,k) = 0. _d 0
395             CALL CALC_GT(            ENDDO
396       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,           ENDDO
      I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,  
      I         KappaRT,  
      U         fVerT,  
      I         myTime, myThid)  
            CALL TIMESTEP_TRACER(  
      I         bi,bj,iMin,iMax,jMin,jMax,k,  
      I         theta, gT,  
      U         gTnm1,  
      I         myIter, myThid)  
          ENDIF  
          IF ( saltStepping ) THEN  
            CALL CALC_GS(  
      I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,  
      I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,  
      I         KappaRS,  
      U         fVerS,  
      I         myTime, myThid)  
            CALL TIMESTEP_TRACER(  
      I         bi,bj,iMin,iMax,jMin,jMax,k,  
      I         salt, gS,  
      U         gSnm1,  
      I         myIter, myThid)  
          ENDIF  
   
 #ifdef   ALLOW_OBCS  
 C--      Apply open boundary conditions  
          IF (useOBCS) THEN  
            CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )  
          END IF  
 #endif   /* ALLOW_OBCS */  
   
 C--      Freeze water  
          IF (allowFreezing) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE gTNm1(:,:,k,bi,bj) = comlev1_bibj_k  
 CADJ &   , key = kkey, byte = isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
             CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )  
          END IF  
   
 C--     end of thermodynamic k loop (Nr:1)  
397          ENDDO          ENDDO
398    #endif
399    
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 C? Patrick? What about this one?  
 cph Keys iikey and idkey don't seem to be needed  
 cph since storing occurs on different tape for each  
 cph impldiff call anyways.  
 cph Thus, common block comlev1_impl isn't needed either.  
 cph Storing below needed in the case useGMREDI.  
         iikey = (ikey-1)*maximpl  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 C--     Implicit diffusion  
         IF (implicitDiffusion) THEN  
   
          IF (tempStepping) THEN  
400  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
401              idkey = iikey + 1  CADJ STORE KappaRU(:,:,:)
402  CADJ STORE gTNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ &     = comlev1_bibj, key=idynkey, byte=isbyte
403    CADJ STORE KappaRV(:,:,:)
404    CADJ &     = comlev1_bibj, key=idynkey, byte=isbyte
405  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
             CALL IMPLDIFF(  
      I         bi, bj, iMin, iMax, jMin, jMax,  
      I         deltaTtracer, KappaRT, recip_HFacC,  
      U         gTNm1,  
      I         myThid )  
          ENDIF  
406    
407           IF (saltStepping) THEN  #ifdef ALLOW_OBCS
408  #ifdef ALLOW_AUTODIFF_TAMC  C--   For Stevens boundary conditions velocities need to be extrapolated
409           idkey = iikey + 2  C     (copied) to a narrow strip outside the domain
410  CADJ STORE gSNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte           IF ( useOBCS ) THEN
411  #endif /* ALLOW_AUTODIFF_TAMC */            CALL OBCS_COPY_UV_N(
412              CALL IMPLDIFF(       U         uVel(1-OLx,1-OLy,1,bi,bj),
413       I         bi, bj, iMin, iMax, jMin, jMax,       U         vVel(1-OLx,1-OLy,1,bi,bj),
414       I         deltaTtracer, KappaRS, recip_HFacC,       I         Nr, bi, bj, myThid )
      U         gSNm1,  
      I         myThid )  
415           ENDIF           ENDIF
416    #endif /* ALLOW_OBCS */
 #ifdef   ALLOW_OBCS  
 C--      Apply open boundary conditions  
          IF (useOBCS) THEN  
            DO K=1,Nr  
              CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )  
            ENDDO  
          END IF  
 #endif   /* ALLOW_OBCS */  
   
 C--     End If implicitDiffusion  
         ENDIF  
   
 C--     Start computation of dynamics  
         iMin = 1-OLx+2  
         iMax = sNx+OLx-1  
         jMin = 1-OLy+2  
         jMax = sNy+OLy-1  
   
 C--     Explicit part of the Surface Potentiel Gradient (add in TIMESTEP)  
 C       (note: this loop will be replaced by CALL CALC_GRAD_ETA)  
         IF (implicSurfPress.NE.1.) THEN  
           CALL CALC_GRAD_PHI_SURF(  
      I         bi,bj,iMin,iMax,jMin,jMax,  
      I         etaN,  
      O         phiSurfX,phiSurfY,  
      I         myThid )                          
         ENDIF  
417    
418  C--     Start of dynamics loop  C--     Start of dynamics loop
419          DO k=1,Nr          DO k=1,Nr
# Line 582  C--       kup    Cycles through 1,2 to p Line 423  C--       kup    Cycles through 1,2 to p
423  C--       kDown  Cycles through 2,1 to point to current layer  C--       kDown  Cycles through 2,1 to point to current layer
424    
425            km1  = MAX(1,k-1)            km1  = MAX(1,k-1)
426              kp1  = MIN(k+1,Nr)
427            kup  = 1+MOD(k+1,2)            kup  = 1+MOD(k+1,2)
428            kDown= 1+MOD(k,2)            kDown= 1+MOD(k,2)
429    
430  C--      Integrate hydrostatic balance for phiHyd with BC of  #ifdef ALLOW_AUTODIFF_TAMC
431             kkey = (idynkey-1)*Nr + k
432    c
433    CADJ STORE totPhiHyd (:,:,k,bi,bj)
434    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
435    CADJ STORE phiHydLow (:,:,bi,bj)
436    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
437    CADJ STORE theta (:,:,k,bi,bj)
438    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
439    CADJ STORE salt  (:,:,k,bi,bj)
440    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
441    CADJ STORE gT(:,:,k,bi,bj)
442    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
443    CADJ STORE gS(:,:,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  guDissip (:,:)
452    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
453    CADJ STORE  gvDissip (:,:)
454    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
455    CADJ STORE  fVerU (:,:,:)
456    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
457    CADJ STORE  fVerV (:,:,:)
458    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
459    CADJ STORE gU(:,:,k,bi,bj)
460    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
461    CADJ STORE gV(:,:,k,bi,bj)
462    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
463    #  ifndef ALLOW_ADAMSBASHFORTH_3
464    CADJ STORE guNm1(:,:,k,bi,bj)
465    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
466    CADJ STORE gvNm1(:,:,k,bi,bj)
467    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
468    #  else
469    CADJ STORE guNm(:,:,k,bi,bj,1)
470    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
471    CADJ STORE guNm(:,:,k,bi,bj,2)
472    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
473    CADJ STORE gvNm(:,:,k,bi,bj,1)
474    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
475    CADJ STORE gvNm(:,:,k,bi,bj,2)
476    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
477    #  endif
478    #  ifdef ALLOW_CD_CODE
479    CADJ STORE uNM1(:,:,k,bi,bj)
480    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
481    CADJ STORE vNM1(:,:,k,bi,bj)
482    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
483    CADJ STORE uVelD(:,:,k,bi,bj)
484    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
485    CADJ STORE vVelD(:,:,k,bi,bj)
486    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
487    #  endif
488    # endif
489    # ifdef ALLOW_DEPTH_CONTROL
490    CADJ STORE  fVerU (:,:,:)
491    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
492    CADJ STORE  fVerV (:,:,:)
493    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
494    # endif
495    #endif /* ALLOW_AUTODIFF_TAMC */
496    
497    C--      Integrate hydrostatic balance for phiHyd with BC of
498  C        phiHyd(z=0)=0  C        phiHyd(z=0)=0
499  C        distinguishe between Stagger and Non Stagger time stepping           IF ( implicitIntGravWave ) THEN
          IF (staggerTimeStep) THEN  
500             CALL CALC_PHI_HYD(             CALL CALC_PHI_HYD(
501       I        bi,bj,iMin,iMax,jMin,jMax,k,       I        bi,bj,iMin,iMax,jMin,jMax,k,
502       I        gTnm1, gSnm1,       I        gT, gS,
503       U        phiHyd,       U        phiHydF,
504       I        myThid )       O        phiHydC, dPhiHydX, dPhiHydY,
505         I        myTime, myIter, myThid )
506           ELSE           ELSE
507             CALL CALC_PHI_HYD(             CALL CALC_PHI_HYD(
508       I        bi,bj,iMin,iMax,jMin,jMax,k,       I        bi,bj,iMin,iMax,jMin,jMax,k,
509       I        theta, salt,       I        theta, salt,
510       U        phiHyd,       U        phiHydF,
511       I        myThid )       O        phiHydC, dPhiHydX, dPhiHydY,
512         I        myTime, myIter, myThid )
513             ENDIF
514    #ifdef ALLOW_DIAGNOSTICS
515             IF ( dPhiHydDiagIsOn ) THEN
516               tmpFac = -1. _d 0
517               CALL DIAGNOSTICS_SCALE_FILL( dPhiHydX, tmpFac, 1,
518         &                           'Um_dPHdx', k, 1, 2, bi, bj, myThid )
519               CALL DIAGNOSTICS_SCALE_FILL( dPhiHydY, tmpFac, 1,
520         &                           'Vm_dPHdy', k, 1, 2, bi, bj, myThid )
521           ENDIF           ENDIF
522    #endif /* ALLOW_DIAGNOSTICS */
523    
524  C--      Calculate accelerations in the momentum equations (gU, gV, ...)  C--      Calculate accelerations in the momentum equations (gU, gV, ...)
525  C        and step forward storing the result in gUnm1, gVnm1, etc...  C        and step forward storing the result in gU, gV, etc...
526           IF ( momStepping ) THEN           IF ( momStepping ) THEN
527             CALL CALC_MOM_RHS(  #ifdef ALLOW_AUTODIFF_TAMC
528       I         bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,  # ifdef NONLIN_FRSURF
529       I         phiHyd,KappaRU,KappaRV,  #  if (defined ALLOW_MOM_FLUXFORM) && !(defined DISABLE_RSTAR_CODE)
530       U         fVerU, fVerV,  CADJ STORE dWtransC(:,:,bi,bj)
531       I         myTime, myThid)  CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
532    CADJ STORE dWtransU(:,:,bi,bj)
533    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
534    CADJ STORE dWtransV(:,:,bi,bj)
535    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
536    #  endif
537    CADJ STORE fVerU(:,:,:)
538    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
539    CADJ STORE fVerV(:,:,:)
540    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
541    # endif /* NONLIN_FRSURF */
542    #endif /* ALLOW_AUTODIFF_TAMC */
543               IF (.NOT. vectorInvariantMomentum) THEN
544    #ifdef ALLOW_MOM_FLUXFORM
545                  CALL MOM_FLUXFORM(
546         I         bi,bj,k,iMin,iMax,jMin,jMax,
547         I         KappaRU, KappaRV,
548         U         fVerU(1-OLx,1-OLy,kUp),   fVerV(1-OLx,1-OLy,kUp),
549         O         fVerU(1-OLx,1-OLy,kDown), fVerV(1-OLx,1-OLy,kDown),
550         O         guDissip, gvDissip,
551         I         myTime, myIter, myThid)
552    #endif
553               ELSE
554    #ifdef ALLOW_MOM_VECINV
555                 CALL MOM_VECINV(
556         I         bi,bj,k,iMin,iMax,jMin,jMax,
557         I         KappaRU, KappaRV,
558         I         fVerU(1-OLx,1-OLy,kUp),   fVerV(1-OLx,1-OLy,kUp),
559         O         fVerU(1-OLx,1-OLy,kDown), fVerV(1-OLx,1-OLy,kDown),
560         O         guDissip, gvDissip,
561         I         myTime, myIter, myThid)
562    #endif
563               ENDIF
564    C
565             CALL TIMESTEP(             CALL TIMESTEP(
566       I         bi,bj,iMin,iMax,jMin,jMax,k,       I         bi,bj,iMin,iMax,jMin,jMax,k,
567       I         phiHyd, phiSurfX, phiSurfY,       I         dPhiHydX,dPhiHydY, phiSurfX, phiSurfY,
568       I         myIter, myThid)       I         guDissip, gvDissip,
569         I         myTime, myIter, myThid)
570    
 #ifdef   ALLOW_OBCS  
 C--      Apply open boundary conditions  
          IF (useOBCS) THEN  
            CALL OBCS_APPLY_UV( bi, bj, k, gUnm1, gVnm1, myThid )  
          END IF  
 #endif   /* ALLOW_OBCS */  
   
 #ifdef   ALLOW_AUTODIFF_TAMC  
 #ifdef   INCLUDE_CD_CODE  
          ELSE  
            DO j=1-OLy,sNy+OLy  
              DO i=1-OLx,sNx+OLx  
                guCD(i,j,k,bi,bj) = 0.0  
                gvCD(i,j,k,bi,bj) = 0.0  
              END DO  
            END DO  
 #endif   /* INCLUDE_CD_CODE */  
 #endif   /* ALLOW_AUTODIFF_TAMC */  
571           ENDIF           ENDIF
572    
   
573  C--     end of dynamics k loop (1:Nr)  C--     end of dynamics k loop (1:Nr)
574          ENDDO          ENDDO
575    
576    C--     Implicit Vertical advection & viscosity
577    #if (defined (INCLUDE_IMPLVERTADV_CODE) && \
578  C--     Implicit viscosity       defined (ALLOW_MOM_COMMON) && !(defined ALLOW_AUTODIFF_TAMC))
579          IF (implicitViscosity.AND.momStepping) THEN          IF ( momImplVertAdv ) THEN
580              CALL MOM_U_IMPLICIT_R( kappaRU,
581         I                           bi, bj, myTime, myIter, myThid )
582              CALL MOM_V_IMPLICIT_R( kappaRV,
583         I                           bi, bj, myTime, myIter, myThid )
584            ELSEIF ( implicitViscosity ) THEN
585    #else /* INCLUDE_IMPLVERTADV_CODE */
586            IF     ( implicitViscosity ) THEN
587    #endif /* INCLUDE_IMPLVERTADV_CODE */
588  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
589            idkey = iikey + 3  CADJ STORE KappaRU(:,:,:) = comlev1_bibj , key=idynkey, byte=isbyte
590  CADJ STORE gUNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gU(:,:,:,bi,bj) = comlev1_bibj , key=idynkey, byte=isbyte
591  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
592            CALL IMPLDIFF(            CALL IMPLDIFF(
593       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
594       I         deltaTmom, KappaRU,recip_HFacW,       I         -1, KappaRU, recip_hFacW(1-OLx,1-OLy,1,bi,bj),
595       U         gUNm1,       U         gU,
596       I         myThid )       I         myThid )
597  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
598            idkey = iikey + 4  CADJ STORE KappaRV(:,:,:) = comlev1_bibj , key=idynkey, byte=isbyte
599  CADJ STORE gVNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gV(:,:,:,bi,bj) = comlev1_bibj , key=idynkey, byte=isbyte
600  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
601            CALL IMPLDIFF(            CALL IMPLDIFF(
602       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
603       I         deltaTmom, KappaRV,recip_HFacS,       I         -2, KappaRV, recip_hFacS(1-OLx,1-OLy,1,bi,bj),
604       U         gVNm1,       U         gV,
605       I         myThid )       I         myThid )
606            ENDIF
607    
608  #ifdef   ALLOW_OBCS  #ifdef ALLOW_OBCS
609  C--      Apply open boundary conditions  C--      Apply open boundary conditions
610           IF (useOBCS) THEN          IF ( useOBCS ) THEN
611             DO K=1,Nr  C--      but first save intermediate velocities to be used in the
612               CALL OBCS_APPLY_UV( bi, bj, k, gUnm1, gVnm1, myThid )  C        next time step for the Stevens boundary conditions
613             ENDDO            CALL OBCS_SAVE_UV_N(
614           END IF       I        bi, bj, iMin, iMax, jMin, jMax, 0,
615  #endif   /* ALLOW_OBCS */       I        gU, gV, myThid )
616              CALL OBCS_APPLY_UV( bi, bj, 0, gU, gV, myThid )
617            ENDIF
618    #endif /* ALLOW_OBCS */
619    
620  #ifdef    INCLUDE_CD_CODE  #ifdef    ALLOW_CD_CODE
621            IF (implicitViscosity.AND.useCDscheme) THEN
622  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
623            idkey = iikey + 5  CADJ STORE vVelD(:,:,:,bi,bj) = comlev1_bibj , key=idynkey, byte=isbyte
 CADJ STORE vVelD(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  
624  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
625            CALL IMPLDIFF(            CALL IMPLDIFF(
626       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
627       I         deltaTmom, KappaRU,recip_HFacW,       I         0, KappaRU, recip_hFacW(1-OLx,1-OLy,1,bi,bj),
628       U         vVelD,       U         vVelD,
629       I         myThid )       I         myThid )
630  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
631            idkey = iikey + 6  CADJ STORE uVelD(:,:,:,bi,bj) = comlev1_bibj , key=idynkey, byte=isbyte
 CADJ STORE uVelD(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  
632  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
633            CALL IMPLDIFF(            CALL IMPLDIFF(
634       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
635       I         deltaTmom, KappaRV,recip_HFacS,       I         0, KappaRV, recip_hFacS(1-OLx,1-OLy,1,bi,bj),
636       U         uVelD,       U         uVelD,
637       I         myThid )       I         myThid )
 #endif    /* INCLUDE_CD_CODE */  
 C--     End If implicitViscosity.AND.momStepping  
638          ENDIF          ENDIF
639    #endif    /* ALLOW_CD_CODE */
640  Cjmc : add for phiHyd output <- but not working if multi tile per CPU  C--     End implicit Vertical advection & viscosity
641  c       IF ( DIFFERENT_MULTIPLE(dumpFreq,myTime+deltaTClock,myTime)  
642  c    &  .AND. buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
643  c         WRITE(suff,'(I10.10)') myIter+1  
644  c         CALL WRITE_FLD_XYZ_RL('PH.',suff,phiHyd,myIter+1,myThid)  #ifdef ALLOW_NONHYDROSTATIC
645  c       ENDIF  C--   Step forward W field in N-H algorithm
646  Cjmc(end)          IF ( nonHydrostatic ) THEN
647    #ifdef ALLOW_DEBUG
648  #ifdef ALLOW_TIMEAVE           IF (debugMode) CALL DEBUG_CALL('CALC_GW', myThid )
649          IF (taveFreq.GT.0.) THEN  #endif
650            CALL TIMEAVE_CUMULATE(phiHydtave, phiHyd, Nr,           CALL TIMER_START('CALC_GW          [DYNAMICS]',myThid)
651       I                              deltaTclock, bi, bj, myThid)           CALL CALC_GW(
652            IF (ivdc_kappa.NE.0.) THEN       I                 bi,bj, KappaRU, KappaRV,
653              CALL TIMEAVE_CUMULATE(ConvectCountTave, ConvectCount, Nr,       I                 myTime, myIter, myThid )
      I                              deltaTclock, bi, bj, myThid)  
           ENDIF  
654          ENDIF          ENDIF
655  #endif /* ALLOW_TIMEAVE */          IF ( nonHydrostatic.OR.implicitIntGravWave )
656         &   CALL TIMESTEP_WVEL( bi,bj, myTime, myIter, myThid )
657            IF ( nonHydrostatic )
658         &   CALL TIMER_STOP ('CALC_GW          [DYNAMICS]',myThid)
659    #endif
660    
661    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
662    
663    C-    end of bi,bj loops
664         ENDDO         ENDDO
665        ENDDO        ENDDO
666    
667    #ifdef ALLOW_OBCS
668          IF (useOBCS) THEN
669            CALL OBCS_EXCHANGES( myThid )
670          ENDIF
671    #endif
672    
673    Cml(
674    C     In order to compare the variance of phiHydLow of a p/z-coordinate
675    C     run with etaH of a z/p-coordinate run the drift of phiHydLow
676    C     has to be removed by something like the following subroutine:
677    C      CALL REMOVE_MEAN_RL( 1, phiHydLow, maskInC, maskInC, rA, drF,
678    C     &                     'phiHydLow', myTime, myThid )
679    Cml)
680    
681    #ifdef ALLOW_DIAGNOSTICS
682          IF ( useDiagnostics ) THEN
683    
684           CALL DIAGNOSTICS_FILL(totPhihyd,'PHIHYD  ',0,Nr,0,1,1,myThid)
685           CALL DIAGNOSTICS_FILL(phiHydLow,'PHIBOT  ',0, 1,0,1,1,myThid)
686    
687           tmpFac = 1. _d 0
688           CALL DIAGNOSTICS_SCALE_FILL(totPhihyd,tmpFac,2,
689         &                                 'PHIHYDSQ',0,Nr,0,1,1,myThid)
690    
691           CALL DIAGNOSTICS_SCALE_FILL(phiHydLow,tmpFac,2,
692         &                                 'PHIBOTSQ',0, 1,0,1,1,myThid)
693    
694          ENDIF
695    #endif /* ALLOW_DIAGNOSTICS */
696    
697    #ifdef ALLOW_DEBUG
698          IF ( debugLevel .GE. debLevD ) THEN
699           CALL DEBUG_STATS_RL(1,EtaN,'EtaN (DYNAMICS)',myThid)
700           CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (DYNAMICS)',myThid)
701           CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (DYNAMICS)',myThid)
702           CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (DYNAMICS)',myThid)
703           CALL DEBUG_STATS_RL(Nr,theta,'Theta (DYNAMICS)',myThid)
704           CALL DEBUG_STATS_RL(Nr,salt,'Salt (DYNAMICS)',myThid)
705           CALL DEBUG_STATS_RL(Nr,gU,'Gu (DYNAMICS)',myThid)
706           CALL DEBUG_STATS_RL(Nr,gV,'Gv (DYNAMICS)',myThid)
707           CALL DEBUG_STATS_RL(Nr,gT,'Gt (DYNAMICS)',myThid)
708           CALL DEBUG_STATS_RL(Nr,gS,'Gs (DYNAMICS)',myThid)
709    #ifndef ALLOW_ADAMSBASHFORTH_3
710           CALL DEBUG_STATS_RL(Nr,guNm1,'GuNm1 (DYNAMICS)',myThid)
711           CALL DEBUG_STATS_RL(Nr,gvNm1,'GvNm1 (DYNAMICS)',myThid)
712           CALL DEBUG_STATS_RL(Nr,gtNm1,'GtNm1 (DYNAMICS)',myThid)
713           CALL DEBUG_STATS_RL(Nr,gsNm1,'GsNm1 (DYNAMICS)',myThid)
714    #endif
715          ENDIF
716    #endif
717    
718    #ifdef DYNAMICS_GUGV_EXCH_CHECK
719    C- jmc: For safety checking only: This Exchange here should not change
720    C       the solution. If solution changes, it means something is wrong,
721    C       but it does not mean that it is less wrong with this exchange.
722          IF ( debugLevel .GE. debLevE ) THEN
723           CALL EXCH_UV_XYZ_RL(gU,gV,.TRUE.,myThid)
724          ENDIF
725    #endif
726    
727    #ifdef ALLOW_DEBUG
728          IF (debugMode) CALL DEBUG_LEAVE( 'DYNAMICS', myThid )
729    #endif
730    
731        RETURN        RETURN
732        END        END
Legend:
Removed from v.1.66  
changed lines
  Added in v.1.164
  ViewVC Help
Powered by ViewVC 1.1.22