/[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.95 by heimbach, Fri Feb 28 02:20:52 2003 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
5    
6    CBOP
7    C     !ROUTINE: DYNAMICS
8    C     !INTERFACE:
9        SUBROUTINE DYNAMICS(myTime, myIter, myThid)        SUBROUTINE DYNAMICS(myTime, myIter, myThid)
10  C     /==========================================================\  C     !DESCRIPTION: \bv
11  C     | SUBROUTINE DYNAMICS                                      |  C     *==========================================================*
12  C     | o Controlling routine for the explicit part of the model |  C     | SUBROUTINE DYNAMICS                                      
13  C     |   dynamics.                                              |  C     | o Controlling routine for the explicit part of the model  
14  C     |==========================================================|  C     |   dynamics.                                              
15  C     | This routine evaluates the "dynamics" terms for each     |  C     *==========================================================*
16  C     | block of ocean in turn. Because the blocks of ocean have |  C     | This routine evaluates the "dynamics" terms for each      
17  C     | overlap regions they are independent of one another.     |  C     | block of ocean in turn. Because the blocks of ocean have  
18  C     | If terms involving lateral integrals are needed in this  |  C     | overlap regions they are independent of one another.      
19  C     | routine care will be needed. Similarly finite-difference |  C     | If terms involving lateral integrals are needed in this  
20  C     | operations with stencils wider than the overlap region   |  C     | routine care will be needed. Similarly finite-difference  
21  C     | require special consideration.                           |  C     | operations with stencils wider than the overlap region    
22  C     | Notes                                                    |  C     | require special consideration.                            
23  C     | =====                                                    |  C     | The algorithm...
24  C     | C*P* comments indicating place holders for which code is |  C     |
25  C     |      presently being developed.                          |  C     | "Correction Step"
26  C     \==========================================================/  C     | =================
27    C     | Here we update the horizontal velocities with the surface
28    C     | pressure such that the resulting flow is either consistent
29    C     | with the free-surface evolution or the rigid-lid:
30    C     |   U[n] = U* + dt x d/dx P
31    C     |   V[n] = V* + dt x d/dy P
32    C     |
33    C     | "Calculation of Gs"
34    C     | ===================
35    C     | This is where all the accelerations and tendencies (ie.
36    C     | physics, parameterizations etc...) are calculated
37    C     |   rho = rho ( theta[n], salt[n] )
38    C     |   b   = b(rho, theta)
39    C     |   K31 = K31 ( rho )
40    C     |   Gu[n] = Gu( u[n], v[n], wVel, b, ... )
41    C     |   Gv[n] = Gv( u[n], v[n], wVel, b, ... )
42    C     |   Gt[n] = Gt( theta[n], u[n], v[n], wVel, K31, ... )
43    C     |   Gs[n] = Gs( salt[n], u[n], v[n], wVel, K31, ... )
44    C     |
45    C     | "Time-stepping" or "Prediction"
46    C     | ================================
47    C     | The models variables are stepped forward with the appropriate
48    C     | time-stepping scheme (currently we use Adams-Bashforth II)
49    C     | - For momentum, the result is always *only* a "prediction"
50    C     | in that the flow may be divergent and will be "corrected"
51    C     | later with a surface pressure gradient.
52    C     | - Normally for tracers the result is the new field at time
53    C     | level [n+1} *BUT* in the case of implicit diffusion the result
54    C     | is also *only* a prediction.
55    C     | - We denote "predictors" with an asterisk (*).
56    C     |   U* = U[n] + dt x ( 3/2 Gu[n] - 1/2 Gu[n-1] )
57    C     |   V* = V[n] + dt x ( 3/2 Gv[n] - 1/2 Gv[n-1] )
58    C     |   theta[n+1] = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
59    C     |   salt[n+1] = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
60    C     | With implicit diffusion:
61    C     |   theta* = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
62    C     |   salt* = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
63    C     |   (1 + dt * K * d_zz) theta[n] = theta*
64    C     |   (1 + dt * K * d_zz) salt[n] = salt*
65    C     |
66    C     *==========================================================*
67    C     \ev
68    C     !USES:
69        IMPLICIT NONE        IMPLICIT NONE
   
70  C     == Global variables ===  C     == Global variables ===
71  #include "SIZE.h"  #include "SIZE.h"
72  #include "EEPARAMS.h"  #include "EEPARAMS.h"
73  #include "PARAMS.h"  #include "PARAMS.h"
74  #include "DYNVARS.h"  #include "DYNVARS.h"
75  #include "GRID.h"  #include "GRID.h"
76    #ifdef ALLOW_PASSIVE_TRACER
77    #include "TR1.h"
78    #endif
79  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
80  # include "tamc.h"  # include "tamc.h"
81  # include "tamc_keys.h"  # include "tamc_keys.h"
82  #endif /* ALLOW_AUTODIFF_TAMC */  # include "FFIELDS.h"
83    # include "EOS.h"
84  #ifdef ALLOW_KPP  # ifdef ALLOW_KPP
85  # include "KPP.h"  #  include "KPP.h"
86  #endif  # endif
87    #endif /* ALLOW_AUTODIFF_TAMC */
88  #ifdef ALLOW_TIMEAVE  
89  #include "TIMEAVE_STATV.h"  C     !CALLING SEQUENCE:
90  #endif  C     DYNAMICS()
91    C      |
92    C      |-- CALC_GRAD_PHI_SURF
93    C      |
94    C      |-- CALC_VISCOSITY
95    C      |
96    C      |-- CALC_PHI_HYD  
97    C      |
98    C      |-- MOM_FLUXFORM  
99    C      |
100    C      |-- MOM_VECINV    
101    C      |
102    C      |-- TIMESTEP      
103    C      |
104    C      |-- OBCS_APPLY_UV
105    C      |
106    C      |-- IMPLDIFF      
107    C      |
108    C      |-- OBCS_APPLY_UV
109    C      |
110    C      |-- CALL TIMEAVE_CUMUL_1T
111    C      |-- CALL DEBUG_STATS_RL
112    
113    C     !INPUT/OUTPUT PARAMETERS:
114  C     == Routine arguments ==  C     == Routine arguments ==
115  C     myTime - Current time in simulation  C     myTime - Current time in simulation
116  C     myIter - Current iteration number in simulation  C     myIter - Current iteration number in simulation
# Line 51  C     myThid - Thread number for this in Line 119  C     myThid - Thread number for this in
119        INTEGER myIter        INTEGER myIter
120        INTEGER myThid        INTEGER myThid
121    
122    C     !LOCAL VARIABLES:
123  C     == Local variables  C     == Local variables
 C     xA, yA                 - Per block temporaries holding face areas  
 C     uTrans, vTrans, rTrans - Per block temporaries holding flow  
 C                              transport  
 C                              o uTrans: Zonal transport  
 C                              o vTrans: Meridional transport  
 C                              o rTrans: Vertical transport  
 C     maskC,maskUp             o maskC: land/water mask for tracer cells  
 C                              o maskUp: land/water mask for W points  
124  C     fVer[STUV]               o fVer: Vertical flux term - note fVer  C     fVer[STUV]               o fVer: Vertical flux term - note fVer
125  C                                      is "pipelined" in the vertical  C                                      is "pipelined" in the vertical
126  C                                      so we need an fVer for each  C                                      so we need an fVer for each
127  C                                      variable.  C                                      variable.
128  C     rhoK, rhoKM1   - Density at current level, and level above  C     phiHydC    :: hydrostatic potential anomaly at cell center
129  C     phiHyd         - Hydrostatic part of the potential phiHydi.  C                   In z coords phiHyd is the hydrostatic potential
130  C                      In z coords phiHydiHyd is the hydrostatic  C                      (=pressure/rho0) anomaly
131  C                      Potential (=pressure/rho0) anomaly  C                   In p coords phiHyd is the geopotential height anomaly.
132  C                      In p coords phiHydiHyd is the geopotential  C     phiHydF    :: hydrostatic potential anomaly at middle between 2 centers
133  C                      surface height anomaly.  C     dPhiHydX,Y :: Gradient (X & Y directions) of hydrostatic potential anom.
134  C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)  C     phiSurfX,  ::  gradient of Surface potential (Pressure/rho, ocean)
135  C     phiSurfY             or geopotentiel (atmos) in X and Y direction  C     phiSurfY             or geopotential (atmos) in X and Y direction
 C     KappaRT,       - Total diffusion in vertical for T and S.  
 C     KappaRS          (background + spatially varying, isopycnal term).  
136  C     iMin, iMax     - Ranges and sub-block indices on which calculations  C     iMin, iMax     - Ranges and sub-block indices on which calculations
137  C     jMin, jMax       are applied.  C     jMin, jMax       are applied.
138  C     bi, bj  C     bi, bj
139  C     k, kup,        - Index for layer above and below. kup and kDown  C     k, kup,        - Index for layer above and below. kup and kDown
140  C     kDown, km1       are switched with layer to be the appropriate  C     kDown, km1       are switched with layer to be the appropriate
141  C                      index into fVerTerm.  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)  
142        _RL fVerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
143        _RL fVerV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
144        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL phiHydF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
145        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiHydC (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
146        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
147          _RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
148        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
149        _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)  
       _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)  
150        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
151        _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)  
152    
153        INTEGER iMin, iMax        INTEGER iMin, iMax
154        INTEGER jMin, jMax        INTEGER jMin, jMax
155        INTEGER bi, bj        INTEGER bi, bj
156        INTEGER i, j        INTEGER i, j
157        INTEGER k, km1, kup, kDown        INTEGER k, km1, kp1, kup, kDown
158    
159  Cjmc : add for phiHyd output <- but not working if multi tile per CPU        LOGICAL  DIFFERENT_MULTIPLE
160  c     CHARACTER*(MAX_LEN_MBUF) suff        EXTERNAL DIFFERENT_MULTIPLE
 c     LOGICAL  DIFFERENT_MULTIPLE  
 c     EXTERNAL DIFFERENT_MULTIPLE  
 Cjmc(end)  
161    
162  C---    The algorithm...  C---    The algorithm...
163  C  C
# Line 162  C         salt* = salt[n] + dt x ( 3/2 G Line 202  C         salt* = salt[n] + dt x ( 3/2 G
202  C         (1 + dt * K * d_zz) theta[n] = theta*  C         (1 + dt * K * d_zz) theta[n] = theta*
203  C         (1 + dt * K * d_zz) salt[n] = salt*  C         (1 + dt * K * d_zz) salt[n] = salt*
204  C---  C---
205    CEOP
 #ifdef ALLOW_AUTODIFF_TAMC  
 C--   dummy statement to end declaration part  
       ikey = 1  
 #endif /* ALLOW_AUTODIFF_TAMC */  
206    
207  C--   Set up work arrays with valid (i.e. not NaN) values  C--   Set up work arrays with valid (i.e. not NaN) values
208  C     These inital values do not alter the numerical results. They  C     These inital values do not alter the numerical results. They
# Line 175  C     point numbers. This prevents spuri Line 211  C     point numbers. This prevents spuri
211  C     uninitialised but inert locations.  C     uninitialised but inert locations.
212        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
213         DO i=1-OLx,sNx+OLx         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  
214          phiSurfX(i,j) = 0. _d 0          phiSurfX(i,j) = 0. _d 0
215          phiSurfY(i,j) = 0. _d 0          phiSurfY(i,j) = 0. _d 0
216         ENDDO         ENDDO
217        ENDDO        ENDDO
218    
219    C-- Call to routine for calculation of
220    C   Eliassen-Palm-flux-forced U-tendency,
221    C   if desired:
222    #ifdef INCLUDE_EP_FORCING_CODE
223          CALL CALC_EP_FORCING(myThid)
224    #endif
225    
226  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
227  C--   HPF directive to help TAMC  C--   HPF directive to help TAMC
# Line 205  CHPF$ INDEPENDENT Line 232  CHPF$ INDEPENDENT
232    
233  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
234  C--    HPF directive to help TAMC  C--    HPF directive to help TAMC
235  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS,fVerU,fVerV  CHPF$  INDEPENDENT, NEW (fVerU,fVerV
236  CHPF$&                  ,phiHyd,utrans,vtrans,maskc,xA,yA  CHPF$&                  ,phiHydF
237  CHPF$&                  ,KappaRT,KappaRS,KappaRU,KappaRV  CHPF$&                  ,KappaRU,KappaRV
238  CHPF$&                  )  CHPF$&                  )
239  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
240    
# Line 216  CHPF$&                  ) Line 243  CHPF$&                  )
243  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
244            act1 = bi - myBxLo(myThid)            act1 = bi - myBxLo(myThid)
245            max1 = myBxHi(myThid) - myBxLo(myThid) + 1            max1 = myBxHi(myThid) - myBxLo(myThid) + 1
   
246            act2 = bj - myByLo(myThid)            act2 = bj - myByLo(myThid)
247            max2 = myByHi(myThid) - myByLo(myThid) + 1            max2 = myByHi(myThid) - myByLo(myThid) + 1
   
248            act3 = myThid - 1            act3 = myThid - 1
249            max3 = nTx*nTy            max3 = nTx*nTy
   
250            act4 = ikey_dynamics - 1            act4 = ikey_dynamics - 1
251              idynkey = (act1 + 1) + act2*max1
           ikey = (act1 + 1) + act2*max1  
252       &                      + act3*max1*max2       &                      + act3*max1*max2
253       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
254  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
255    
256  C--     Set up work arrays that need valid initial values  C--     Set up work arrays that need valid initial values
         DO j=1-OLy,sNy+OLy  
          DO i=1-OLx,sNx+OLx  
           rTrans(i,j)   = 0. _d 0  
           fVerT (i,j,1) = 0. _d 0  
           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  
   
257          DO k=1,Nr          DO k=1,Nr
258           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
259            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
260  C This is currently also used by IVDC and Diagnostics             KappaRU(i,j,k) = 0. _d 0
261             ConvectCount(i,j,k) = 0.             KappaRV(i,j,k) = 0. _d 0
            KappaRT(i,j,k) = 0. _d 0  
            KappaRS(i,j,k) = 0. _d 0  
262            ENDDO            ENDDO
263           ENDDO           ENDDO
264          ENDDO          ENDDO
265            DO j=1-OLy,sNy+OLy
266          iMin = 1-OLx+1           DO i=1-OLx,sNx+OLx
267          iMax = sNx+OLx            fVerU  (i,j,1) = 0. _d 0
268          jMin = 1-OLy+1            fVerU  (i,j,2) = 0. _d 0
269          jMax = sNy+OLy            fVerV  (i,j,1) = 0. _d 0
270              fVerV  (i,j,2) = 0. _d 0
271              phiHydF (i,j)  = 0. _d 0
272  #ifdef ALLOW_AUTODIFF_TAMC            phiHydC (i,j)  = 0. _d 0
273  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte            dPhiHydX(i,j)  = 0. _d 0
274  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte            dPhiHydY(i,j)  = 0. _d 0
275  CADJ STORE uvel(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte           ENDDO
 CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 C--     Start of diagnostic loop  
         DO k=Nr,1,-1  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 C? Patrick, is this formula correct now that we change the loop range?  
 C? Do we still need this?  
 cph kkey formula corrected.  
 cph Needed for rhok, rhokm1, in the case useGMREDI.  
          kkey = (ikey-1)*Nr + k  
 CADJ STORE rhokm1(:,:) = comlev1_bibj_k , key = kkey, byte = isbyte  
 CADJ STORE rhok  (:,:) = comlev1_bibj_k , key = kkey, byte = isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 C--       Integrate continuity vertically for vertical velocity  
           CALL INTEGRATE_FOR_W(  
      I                         bi, bj, k, uVel, vVel,  
      O                         wVel,  
      I                         myThid )  
   
 #ifdef    ALLOW_OBCS  
 #ifdef    ALLOW_NONHYDROSTATIC  
 C--       Apply OBC to W if in N-H mode  
           IF (useOBCS.AND.nonHydrostatic) THEN  
             CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )  
           ENDIF  
 #endif    /* ALLOW_NONHYDROSTATIC */  
 #endif    /* ALLOW_OBCS */  
   
 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)  
276          ENDDO          ENDDO
277    
278  #ifdef  ALLOW_OBCS  C--     Start computation of dynamics
279  C--     Calculate future values on open boundaries          iMin = 0
280          IF (useOBCS) THEN          iMax = sNx+1
281            CALL OBCS_CALC( bi, bj, myTime+deltaT,          jMin = 0
282       I            uVel, vVel, wVel, theta, salt,          jMax = sNy+1
      I            myThid )  
         ENDIF  
 #endif  /* ALLOW_OBCS */  
283    
 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  
284  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
285          ELSE  CADJ STORE wvel (:,:,:,bi,bj) =
286            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  
287  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
         ENDIF  
 #endif  /* ALLOW_GMREDI */  
288    
289  #ifdef  ALLOW_KPP  C--     Explicit part of the Surface Potentiel Gradient (add in TIMESTEP)
290  C--     Compute KPP mixing coefficients  C       (note: this loop will be replaced by CALL CALC_GRAD_ETA)
291          IF (useKPP) THEN          IF (implicSurfPress.NE.1.) THEN
292            CALL KPP_CALC(            CALL CALC_GRAD_PHI_SURF(
293       I                  bi, bj, myTime, myThid )       I         bi,bj,iMin,iMax,jMin,jMax,
294  #ifdef ALLOW_AUTODIFF_TAMC       I         etaN,
295          ELSE       O         phiSurfX,phiSurfY,
296            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 */  
297          ENDIF          ENDIF
298    
299  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
300  CADJ STORE KPPghat   (:,:,:,bi,bj)  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=idynkey, byte=isbyte
301  CADJ &   , KPPviscAz (:,:,:,bi,bj)  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=idynkey, byte=isbyte
302  CADJ &   , KPPdiffKzT(:,:,:,bi,bj)  #ifdef ALLOW_KPP
303  CADJ &   , KPPdiffKzS(:,:,:,bi,bj)  CADJ STORE KPPviscAz (:,:,:,bi,bj)
304  CADJ &   , KPPfrac   (:,:  ,bi,bj)  CADJ &                 = comlev1_bibj, key=idynkey, byte=isbyte
305  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  
306  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
307    
 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)  
   
308  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL
309  C--      Calculate the total vertical diffusivity  C--      Calculate the total vertical diffusivity
310           CALL CALC_DIFFUSIVITY(          DO k=1,Nr
311             CALL CALC_VISCOSITY(
312       I        bi,bj,iMin,iMax,jMin,jMax,k,       I        bi,bj,iMin,iMax,jMin,jMax,k,
313       I        maskC,maskup,       O        KappaRU,KappaRV,
      O        KappaRT,KappaRS,KappaRU,KappaRV,  
314       I        myThid)       I        myThid)
315           ENDDO
316  #endif  #endif
317    
 C--      Calculate active tracer tendencies (gT,gS,...)  
 C        and step forward storing result in gTnm1, gSnm1, etc.  
          IF ( tempStepping ) THEN  
            CALL CALC_GT(  
      I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,  
      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)  
         ENDDO  
   
   
 #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  
 #ifdef ALLOW_AUTODIFF_TAMC  
             idkey = iikey + 1  
 CADJ STORE gTNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
             CALL IMPLDIFF(  
      I         bi, bj, iMin, iMax, jMin, jMax,  
      I         deltaTtracer, KappaRT, recip_HFacC,  
      U         gTNm1,  
      I         myThid )  
          ENDIF  
   
          IF (saltStepping) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
          idkey = iikey + 2  
 CADJ STORE gSNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
             CALL IMPLDIFF(  
      I         bi, bj, iMin, iMax, jMin, jMax,  
      I         deltaTtracer, KappaRS, recip_HFacC,  
      U         gSNm1,  
      I         myThid )  
          ENDIF  
   
 #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  
   
318  C--     Start of dynamics loop  C--     Start of dynamics loop
319          DO k=1,Nr          DO k=1,Nr
320    
# Line 582  C--       kup    Cycles through 1,2 to p Line 323  C--       kup    Cycles through 1,2 to p
323  C--       kDown  Cycles through 2,1 to point to current layer  C--       kDown  Cycles through 2,1 to point to current layer
324    
325            km1  = MAX(1,k-1)            km1  = MAX(1,k-1)
326              kp1  = MIN(k+1,Nr)
327            kup  = 1+MOD(k+1,2)            kup  = 1+MOD(k+1,2)
328            kDown= 1+MOD(k,2)            kDown= 1+MOD(k,2)
329    
330    #ifdef ALLOW_AUTODIFF_TAMC
331             kkey = (idynkey-1)*Nr + k
332    CADJ STORE totphihyd (:,:,k,bi,bj)
333    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
334    #endif /* ALLOW_AUTODIFF_TAMC */
335    
336  C--      Integrate hydrostatic balance for phiHyd with BC of  C--      Integrate hydrostatic balance for phiHyd with BC of
337  C        phiHyd(z=0)=0  C        phiHyd(z=0)=0
338  C        distinguishe between Stagger and Non Stagger time stepping  C        distinguishe between Stagger and Non Stagger time stepping
339           IF (staggerTimeStep) THEN           IF (staggerTimeStep) THEN
340             CALL CALC_PHI_HYD(             CALL CALC_PHI_HYD(
341       I        bi,bj,iMin,iMax,jMin,jMax,k,       I        bi,bj,iMin,iMax,jMin,jMax,k,
342       I        gTnm1, gSnm1,       I        gT, gS,
343       U        phiHyd,       U        phiHydF,
344       I        myThid )       O        phiHydC, dPhiHydX, dPhiHydY,
345         I        myTime, myIter, myThid )
346           ELSE           ELSE
347             CALL CALC_PHI_HYD(             CALL CALC_PHI_HYD(
348       I        bi,bj,iMin,iMax,jMin,jMax,k,       I        bi,bj,iMin,iMax,jMin,jMax,k,
349       I        theta, salt,       I        theta, salt,
350       U        phiHyd,       U        phiHydF,
351       I        myThid )       O        phiHydC, dPhiHydX, dPhiHydY,
352         I        myTime, myIter, myThid )
353           ENDIF           ENDIF
354    
355  C--      Calculate accelerations in the momentum equations (gU, gV, ...)  C--      Calculate accelerations in the momentum equations (gU, gV, ...)
356  C        and step forward storing the result in gUnm1, gVnm1, etc...  C        and step forward storing the result in gUnm1, gVnm1, etc...
357           IF ( momStepping ) THEN           IF ( momStepping ) THEN
358             CALL CALC_MOM_RHS(  #ifndef DISABLE_MOM_FLUXFORM
359               IF (.NOT. vectorInvariantMomentum) CALL MOM_FLUXFORM(
360       I         bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,
361       I         phiHyd,KappaRU,KappaRV,       I         dPhiHydX,dPhiHydY,KappaRU,KappaRV,
362       U         fVerU, fVerV,       U         fVerU, fVerV,
363       I         myTime, myThid)       I         myTime, myIter, myThid)
364    #endif
365    #ifndef DISABLE_MOM_VECINV
366               IF (vectorInvariantMomentum) CALL MOM_VECINV(
367         I         bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,
368         I         dPhiHydX,dPhiHydY,KappaRU,KappaRV,
369         U         fVerU, fVerV,
370         I         myTime, myIter, myThid)
371    #endif
372             CALL TIMESTEP(             CALL TIMESTEP(
373       I         bi,bj,iMin,iMax,jMin,jMax,k,       I         bi,bj,iMin,iMax,jMin,jMax,k,
374       I         phiHyd, phiSurfX, phiSurfY,       I         dPhiHydX,dPhiHydY, phiSurfX, phiSurfY,
375       I         myIter, myThid)       I         myIter, myThid)
376    
377  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
# Line 639  C--      Apply open boundary conditions Line 398  C--      Apply open boundary conditions
398  C--     end of dynamics k loop (1:Nr)  C--     end of dynamics k loop (1:Nr)
399          ENDDO          ENDDO
400    
   
   
401  C--     Implicit viscosity  C--     Implicit viscosity
402          IF (implicitViscosity.AND.momStepping) THEN          IF (implicitViscosity.AND.momStepping) THEN
403  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
404            idkey = iikey + 3  CADJ STORE gUNm1(:,:,:,bi,bj) = comlev1_bibj , key=idynkey, byte=isbyte
 CADJ STORE gUNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  
405  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
406            CALL IMPLDIFF(            CALL IMPLDIFF(
407       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
# Line 653  CADJ STORE gUNm1(:,:,:,bi,bj) = comlev1_ Line 409  CADJ STORE gUNm1(:,:,:,bi,bj) = comlev1_
409       U         gUNm1,       U         gUNm1,
410       I         myThid )       I         myThid )
411  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
412            idkey = iikey + 4  CADJ STORE gVNm1(:,:,:,bi,bj) = comlev1_bibj , key=idynkey, byte=isbyte
 CADJ STORE gVNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  
413  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
414            CALL IMPLDIFF(            CALL IMPLDIFF(
415       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
# Line 673  C--      Apply open boundary conditions Line 428  C--      Apply open boundary conditions
428    
429  #ifdef    INCLUDE_CD_CODE  #ifdef    INCLUDE_CD_CODE
430  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
431            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  
432  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
433            CALL IMPLDIFF(            CALL IMPLDIFF(
434       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
# Line 682  CADJ STORE vVelD(:,:,:,bi,bj) = comlev1_ Line 436  CADJ STORE vVelD(:,:,:,bi,bj) = comlev1_
436       U         vVelD,       U         vVelD,
437       I         myThid )       I         myThid )
438  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
439            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  
440  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
441            CALL IMPLDIFF(            CALL IMPLDIFF(
442       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
# Line 694  CADJ STORE uVelD(:,:,:,bi,bj) = comlev1_ Line 447  CADJ STORE uVelD(:,:,:,bi,bj) = comlev1_
447  C--     End If implicitViscosity.AND.momStepping  C--     End If implicitViscosity.AND.momStepping
448          ENDIF          ENDIF
449    
 Cjmc : add for phiHyd output <- but not working if multi tile per CPU  
 c       IF ( DIFFERENT_MULTIPLE(dumpFreq,myTime+deltaTClock,myTime)  
 c    &  .AND. buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN  
 c         WRITE(suff,'(I10.10)') myIter+1  
 c         CALL WRITE_FLD_XYZ_RL('PH.',suff,phiHyd,myIter+1,myThid)  
 c       ENDIF  
 Cjmc(end)  
   
 #ifdef ALLOW_TIMEAVE  
         IF (taveFreq.GT.0.) THEN  
           CALL TIMEAVE_CUMULATE(phiHydtave, phiHyd, Nr,  
      I                              deltaTclock, bi, bj, myThid)  
           IF (ivdc_kappa.NE.0.) THEN  
             CALL TIMEAVE_CUMULATE(ConvectCountTave, ConvectCount, Nr,  
      I                              deltaTclock, bi, bj, myThid)  
           ENDIF  
         ENDIF  
 #endif /* ALLOW_TIMEAVE */  
   
450         ENDDO         ENDDO
451        ENDDO        ENDDO
452    
453    Cml(
454    C     In order to compare the variance of phiHydLow of a p/z-coordinate
455    C     run with etaH of a z/p-coordinate run the drift of phiHydLow
456    C     has to be removed by something like the following subroutine:
457    C      CALL REMOVE_MEAN_RL( 1, phiHydLow, maskH, maskH, rA, drF,
458    C     &                'phiHydLow', myThid )
459    Cml)
460    
461    #ifndef DISABLE_DEBUGMODE
462          If (debugMode) THEN
463           CALL DEBUG_STATS_RL(1,EtaN,'EtaN (DYNAMICS)',myThid)
464           CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (DYNAMICS)',myThid)
465           CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (DYNAMICS)',myThid)
466           CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (DYNAMICS)',myThid)
467           CALL DEBUG_STATS_RL(Nr,theta,'Theta (DYNAMICS)',myThid)
468           CALL DEBUG_STATS_RL(Nr,salt,'Salt (DYNAMICS)',myThid)
469           CALL DEBUG_STATS_RL(Nr,Gu,'Gu (DYNAMICS)',myThid)
470           CALL DEBUG_STATS_RL(Nr,Gv,'Gv (DYNAMICS)',myThid)
471           CALL DEBUG_STATS_RL(Nr,Gt,'Gt (DYNAMICS)',myThid)
472           CALL DEBUG_STATS_RL(Nr,Gs,'Gs (DYNAMICS)',myThid)
473           CALL DEBUG_STATS_RL(Nr,GuNm1,'GuNm1 (DYNAMICS)',myThid)
474           CALL DEBUG_STATS_RL(Nr,GvNm1,'GvNm1 (DYNAMICS)',myThid)
475           CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (DYNAMICS)',myThid)
476           CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (DYNAMICS)',myThid)
477          ENDIF
478    #endif
479    
480        RETURN        RETURN
481        END        END
Legend:
Removed from v.1.66  
changed lines
  Added in v.1.95
  ViewVC Help
Powered by ViewVC 1.1.22