/[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.72 by heimbach, Fri Jul 13 14:26:57 2001 UTC revision 1.83.4.3 by heimbach, Wed Apr 17 14:05:34 2002 UTC
# Line 1  Line 1 
 C $Header$  
 C $Name$  
1    
2  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
3    
4    CBOP
5    C     !ROUTINE: DYNAMICS
6    C     !INTERFACE:
7        SUBROUTINE DYNAMICS(myTime, myIter, myThid)        SUBROUTINE DYNAMICS(myTime, myIter, myThid)
8  C     /==========================================================\  C     !DESCRIPTION: \bv
9  C     | SUBROUTINE DYNAMICS                                      |  C     *==========================================================*
10  C     | o Controlling routine for the explicit part of the model |  C     | SUBROUTINE DYNAMICS                                      
11  C     |   dynamics.                                              |  C     | o Controlling routine for the explicit part of the model  
12  C     |==========================================================|  C     |   dynamics.                                              
13  C     | This routine evaluates the "dynamics" terms for each     |  C     *==========================================================*
14  C     | block of ocean in turn. Because the blocks of ocean have |  C     | This routine evaluates the "dynamics" terms for each      
15  C     | overlap regions they are independent of one another.     |  C     | block of ocean in turn. Because the blocks of ocean have  
16  C     | If terms involving lateral integrals are needed in this  |  C     | overlap regions they are independent of one another.      
17  C     | routine care will be needed. Similarly finite-difference |  C     | If terms involving lateral integrals are needed in this  
18  C     | operations with stencils wider than the overlap region   |  C     | routine care will be needed. Similarly finite-difference  
19  C     | require special consideration.                           |  C     | operations with stencils wider than the overlap region    
20  C     | Notes                                                    |  C     | require special consideration.                            
21  C     | =====                                                    |  C     | The algorithm...
22  C     | C*P* comments indicating place holders for which code is |  C     |
23  C     |      presently being developed.                          |  C     | "Correction Step"
24  C     \==========================================================/  C     | =================
25    C     | Here we update the horizontal velocities with the surface
26    C     | pressure such that the resulting flow is either consistent
27    C     | with the free-surface evolution or the rigid-lid:
28    C     |   U[n] = U* + dt x d/dx P
29    C     |   V[n] = V* + dt x d/dy P
30    C     |
31    C     | "Calculation of Gs"
32    C     | ===================
33    C     | This is where all the accelerations and tendencies (ie.
34    C     | physics, parameterizations etc...) are calculated
35    C     |   rho = rho ( theta[n], salt[n] )
36    C     |   b   = b(rho, theta)
37    C     |   K31 = K31 ( rho )
38    C     |   Gu[n] = Gu( u[n], v[n], wVel, b, ... )
39    C     |   Gv[n] = Gv( u[n], v[n], wVel, b, ... )
40    C     |   Gt[n] = Gt( theta[n], u[n], v[n], wVel, K31, ... )
41    C     |   Gs[n] = Gs( salt[n], u[n], v[n], wVel, K31, ... )
42    C     |
43    C     | "Time-stepping" or "Prediction"
44    C     | ================================
45    C     | The models variables are stepped forward with the appropriate
46    C     | time-stepping scheme (currently we use Adams-Bashforth II)
47    C     | - For momentum, the result is always *only* a "prediction"
48    C     | in that the flow may be divergent and will be "corrected"
49    C     | later with a surface pressure gradient.
50    C     | - Normally for tracers the result is the new field at time
51    C     | level [n+1} *BUT* in the case of implicit diffusion the result
52    C     | is also *only* a prediction.
53    C     | - We denote "predictors" with an asterisk (*).
54    C     |   U* = U[n] + dt x ( 3/2 Gu[n] - 1/2 Gu[n-1] )
55    C     |   V* = V[n] + dt x ( 3/2 Gv[n] - 1/2 Gv[n-1] )
56    C     |   theta[n+1] = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
57    C     |   salt[n+1] = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
58    C     | With implicit diffusion:
59    C     |   theta* = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
60    C     |   salt* = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
61    C     |   (1 + dt * K * d_zz) theta[n] = theta*
62    C     |   (1 + dt * K * d_zz) salt[n] = salt*
63    C     |
64    C     *==========================================================*
65    C     \ev
66    C     !USES:
67        IMPLICIT NONE        IMPLICIT NONE
   
68  C     == Global variables ===  C     == Global variables ===
69  #include "SIZE.h"  #include "SIZE.h"
70  #include "EEPARAMS.h"  #include "EEPARAMS.h"
71  #include "PARAMS.h"  #include "PARAMS.h"
72  #include "DYNVARS.h"  #include "DYNVARS.h"
73  #include "GRID.h"  #include "GRID.h"
74    #ifdef ALLOW_PASSIVE_TRACER
75  #include "TR1.h"  #include "TR1.h"
76    #endif
77  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
78  # include "tamc.h"  # include "tamc.h"
79  # include "tamc_keys.h"  # include "tamc_keys.h"
# Line 38  C     == Global variables === Line 81  C     == Global variables ===
81  # ifdef ALLOW_KPP  # ifdef ALLOW_KPP
82  #  include "KPP.h"  #  include "KPP.h"
83  # endif  # endif
 # ifdef ALLOW_GMREDI  
 #  include "GMREDI.h"  
 # endif  
84  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
   
85  #ifdef ALLOW_TIMEAVE  #ifdef ALLOW_TIMEAVE
86  #include "TIMEAVE_STATV.h"  #include "TIMEAVE_STATV.h"
87  #endif  #endif
88    
89    C     !CALLING SEQUENCE:
90    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 55  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     maskUp                   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
# Line 75  C                      In p coords phiHy Line 133  C                      In p coords phiHy
133  C                      surface height anomaly.  C                      surface height anomaly.
134  C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)  C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)
135  C     phiSurfY             or geopotentiel (atmos) in X and Y direction  C     phiSurfY             or geopotentiel (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.
 C     tauAB - Adams-Bashforth timestepping weight: 0=forward ; 1/2=Adams-Bashf.  
       _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 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)  
       _RL fVerTr1 (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 phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
# Line 100  C     tauAB - Adams-Bashforth timesteppi Line 146  C     tauAB - Adams-Bashforth timesteppi
146        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
147        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
148        _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)  
149        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
150        _RL KappaRV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
151        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
152        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
153        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
       _RL tauAB  
   
 C This is currently used by IVDC and Diagnostics  
       _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
154    
155        INTEGER iMin, iMax        INTEGER iMin, iMax
156        INTEGER jMin, jMax        INTEGER jMin, jMax
157        INTEGER bi, bj        INTEGER bi, bj
158        INTEGER i, j        INTEGER i, j
159        INTEGER k, km1, kup, kDown        INTEGER k, km1, kp1, kup, kDown
160    
161  Cjmc : add for phiHyd output <- but not working if multi tile per CPU  Cjmc : add for phiHyd output <- but not working if multi tile per CPU
162  c     CHARACTER*(MAX_LEN_MBUF) suff  c     CHARACTER*(MAX_LEN_MBUF) suff
# Line 167  C         salt* = salt[n] + dt x ( 3/2 G Line 207  C         salt* = salt[n] + dt x ( 3/2 G
207  C         (1 + dt * K * d_zz) theta[n] = theta*  C         (1 + dt * K * d_zz) theta[n] = theta*
208  C         (1 + dt * K * d_zz) salt[n] = salt*  C         (1 + dt * K * d_zz) salt[n] = salt*
209  C---  C---
210    CEOP
 #ifdef ALLOW_AUTODIFF_TAMC  
 C--   dummy statement to end declaration part  
       ikey = 1  
 #endif /* ALLOW_AUTODIFF_TAMC */  
211    
212  C--   Set up work arrays with valid (i.e. not NaN) values  C--   Set up work arrays with valid (i.e. not NaN) values
213  C     These inital values do not alter the numerical results. They  C     These inital values do not alter the numerical results. They
# Line 180  C     point numbers. This prevents spuri Line 216  C     point numbers. This prevents spuri
216  C     uninitialised but inert locations.  C     uninitialised but inert locations.
217        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
218         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  
219          rhoKM1 (i,j) = 0. _d 0          rhoKM1 (i,j) = 0. _d 0
220          rhok   (i,j) = 0. _d 0          rhok   (i,j) = 0. _d 0
221          phiSurfX(i,j) = 0. _d 0          phiSurfX(i,j) = 0. _d 0
# Line 199  C     uninitialised but inert locations. Line 223  C     uninitialised but inert locations.
223         ENDDO         ENDDO
224        ENDDO        ENDDO
225    
   
226  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
227  C--   HPF directive to help TAMC  C--   HPF directive to help TAMC
228  CHPF$ INDEPENDENT  CHPF$ INDEPENDENT
# Line 209  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,xA,yA  CHPF$&                  ,phiHyd
237  CHPF$&                  ,KappaRT,KappaRS,KappaRU,KappaRV  CHPF$&                  ,KappaRU,KappaRV
238  CHPF$&                  )  CHPF$&                  )
239  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
240    
# Line 220  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            ikey = (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
# Line 237  CHPF$&                  ) Line 256  CHPF$&                  )
256  C--     Set up work arrays that need valid initial values  C--     Set up work arrays that need valid initial values
257          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
258           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
259            rTrans (i,j)   = 0. _d 0            DO k=1,Nr
260            fVerT  (i,j,1) = 0. _d 0             phiHyd(i,j,k)  = 0. _d 0
261            fVerT  (i,j,2) = 0. _d 0             KappaRU(i,j,k) = 0. _d 0
262            fVerS  (i,j,1) = 0. _d 0             KappaRV(i,j,k) = 0. _d 0
263            fVerS  (i,j,2) = 0. _d 0            ENDDO
           fVerTr1(i,j,1) = 0. _d 0  
           fVerTr1(i,j,2) = 0. _d 0  
264            fVerU  (i,j,1) = 0. _d 0            fVerU  (i,j,1) = 0. _d 0
265            fVerU  (i,j,2) = 0. _d 0            fVerU  (i,j,2) = 0. _d 0
266            fVerV  (i,j,1) = 0. _d 0            fVerV  (i,j,1) = 0. _d 0
# Line 251  C--     Set up work arrays that need val Line 268  C--     Set up work arrays that need val
268           ENDDO           ENDDO
269          ENDDO          ENDDO
270    
271          DO k=1,Nr  C--     Start computation of dynamics
272           DO j=1-OLy,sNy+OLy          iMin = 1-OLx+2
273            DO i=1-OLx,sNx+OLx          iMax = sNx+OLx-1
274  C This is currently also used by IVDC and Diagnostics          jMin = 1-OLy+2
275             ConvectCount(i,j,k) = 0.          jMax = sNy+OLy-1
            KappaRT(i,j,k) = 0. _d 0  
            KappaRS(i,j,k) = 0. _d 0  
           ENDDO  
          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  
 CADJ STORE tr1  (:,:,:,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  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
             CALL FIND_RHO(  
      I        bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,  
      I        theta, salt,  
      O        rhoK,  
      I        myThid )  
             IF (k.GT.1) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
              CALL FIND_RHO(  
      I        bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,  
      I        theta, salt,  
      O        rhoKm1,  
      I        myThid )  
             ENDIF  
             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)  
         ENDDO  
276    
277  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
 cph avoids recomputation of integrate_for_w  
278  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
279  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
280    
281  #ifdef  ALLOW_OBCS  C--     Explicit part of the Surface Potentiel Gradient (add in TIMESTEP)
282  C--     Calculate future values on open boundaries  C       (note: this loop will be replaced by CALL CALC_GRAD_ETA)
283          IF (useOBCS) THEN          IF (implicSurfPress.NE.1.) THEN
284            CALL OBCS_CALC( bi, bj, myTime+deltaT,            CALL CALC_GRAD_PHI_SURF(
285       I            uVel, vVel, wVel, theta, salt,       I         bi,bj,iMin,iMax,jMin,jMax,
286       I            myThid )       I         etaN,
287          ENDIF       O         phiSurfX,phiSurfY,
288  #endif  /* ALLOW_OBCS */       I         myThid )                        
   
 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_AUTODIFF_TAMC  
 cph needed for KPP  
 CADJ STORE surfacetendencyU(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE surfacetendencyV(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE surfacetendencyS(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE surfacetendencyT(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 #ifdef  ALLOW_GMREDI  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE sigmaX(:,:,:) = comlev1, key=ikey, byte=isbyte  
 CADJ STORE sigmaY(:,:,:) = comlev1, key=ikey, byte=isbyte  
 CADJ STORE sigmaR(:,:,:) = comlev1, key=ikey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
 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  
 #ifdef ALLOW_AUTODIFF_TAMC  
         ELSE  
           DO k=1, Nr  
             CALL GMREDI_CALC_TENSOR_DUMMY(  
      I             bi, bj, iMin, iMax, jMin, jMax, k,  
      I             sigmaX, sigmaY, sigmaR,  
      I             myThid )  
           ENDDO  
 #endif /* ALLOW_AUTODIFF_TAMC */  
         ENDIF  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE Kwx(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE Kwy(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 #endif  /* ALLOW_GMREDI */  
   
 #ifdef  ALLOW_KPP  
 C--     Compute KPP mixing coefficients  
         IF (useKPP) THEN  
           CALL KPP_CALC(  
      I                  bi, bj, myTime, myThid )  
 #ifdef ALLOW_AUTODIFF_TAMC  
         ELSE  
           CALL KPP_CALC_DUMMY(  
      I                  bi, bj, myTime, myThid )  
 #endif /* ALLOW_AUTODIFF_TAMC */  
289          ENDIF          ENDIF
290    
291  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
292  CADJ STORE KPPghat   (:,:,:,bi,bj)  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
293  CADJ &   , KPPviscAz (:,:,:,bi,bj)  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
294  CADJ &   , KPPdiffKzT(:,:,:,bi,bj)  #ifdef ALLOW_KPP
295  CADJ &   , KPPdiffKzS(:,:,:,bi,bj)  CADJ STORE KPPviscAz (:,:,:,bi,bj)
 CADJ &   , KPPfrac   (:,:  ,bi,bj)  
296  CADJ &                 = comlev1_bibj, key=ikey, byte=isbyte  CADJ &                 = comlev1_bibj, key=ikey, byte=isbyte
297  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_KPP */
   
 #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  
 CADJ STORE tr1  (:,:,:,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, bi, bj, 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  
 #ifdef ALLOW_AUTODIFF_TAMC  
 C? Patrick Is this formula correct?  
 cph Yes, but I rewrote it.  
 cph Also, the KappaR? need the index and subscript k!  
          kkey = (ikey-1)*Nr + k  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 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  
   
 C--      Get temporary terms used by tendency routines  
          CALL CALC_COMMON_FACTORS (  
      I        bi,bj,iMin,iMax,jMin,jMax,k,  
      O        xA,yA,uTrans,vTrans,rTrans,maskUp,  
      I        myThid)  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  
298  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
299    
300  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL
301  C--      Calculate the total vertical diffusivity  C--      Calculate the total vertical diffusivity
302           CALL CALC_DIFFUSIVITY(          DO k=1,Nr
303             CALL CALC_VISCOSITY(
304       I        bi,bj,iMin,iMax,jMin,jMax,k,       I        bi,bj,iMin,iMax,jMin,jMax,k,
305       I        maskUp,       O        KappaRU,KappaRV,
      O        KappaRT,KappaRS,KappaRU,KappaRV,  
306       I        myThid)       I        myThid)
307           ENDDO
308  #endif  #endif
309    
 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,  
      I         KappaRT,  
      U         fVerT,  
      I         myTime, myThid)  
            tauAB = 0.5d0 + abEps  
            CALL TIMESTEP_TRACER(  
      I         bi,bj,iMin,iMax,jMin,jMax,k,tauAB,  
      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,  
      I         KappaRS,  
      U         fVerS,  
      I         myTime, myThid)  
            tauAB = 0.5d0 + abEps  
            CALL TIMESTEP_TRACER(  
      I         bi,bj,iMin,iMax,jMin,jMax,k,tauAB,  
      I         salt, gS,  
      U         gSnm1,  
      I         myIter, myThid)  
          ENDIF  
          IF ( tr1Stepping ) THEN  
            CALL CALC_GTR1(  
      I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,  
      I         xA,yA,uTrans,vTrans,rTrans,maskUp,  
      I         KappaRT,  
      U         fVerTr1,  
      I         myTime, myThid)  
            tauAB = 0.5d0 + abEps  
            CALL TIMESTEP_TRACER(  
      I         bi,bj,iMin,iMax,jMin,jMax,k,tauAB,  
      I         Tr1, gTr1,  
      U         gTr1NM1,  
      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  
   
          IF (tr1Stepping) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE gTr1Nm1(:,:,:,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      gTr1Nm1,  
      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  
   
310  C--     Start of dynamics loop  C--     Start of dynamics loop
311          DO k=1,Nr          DO k=1,Nr
312    
# Line 648  C--       kup    Cycles through 1,2 to p Line 315  C--       kup    Cycles through 1,2 to p
315  C--       kDown  Cycles through 2,1 to point to current layer  C--       kDown  Cycles through 2,1 to point to current layer
316    
317            km1  = MAX(1,k-1)            km1  = MAX(1,k-1)
318              kp1  = MIN(k+1,Nr)
319            kup  = 1+MOD(k+1,2)            kup  = 1+MOD(k+1,2)
320            kDown= 1+MOD(k,2)            kDown= 1+MOD(k,2)
321    
322    #ifdef ALLOW_AUTODIFF_TAMC
323             kkey = (ikey-1)*Nr + k
324    #endif /* ALLOW_AUTODIFF_TAMC */
325    
326  C--      Integrate hydrostatic balance for phiHyd with BC of  C--      Integrate hydrostatic balance for phiHyd with BC of
327  C        phiHyd(z=0)=0  C        phiHyd(z=0)=0
328  C        distinguishe between Stagger and Non Stagger time stepping  C        distinguishe between Stagger and Non Stagger time stepping
329           IF (staggerTimeStep) THEN           IF (staggerTimeStep) THEN
330             CALL CALC_PHI_HYD(             CALL CALC_PHI_HYD(
331       I        bi,bj,iMin,iMax,jMin,jMax,k,       I        bi,bj,iMin,iMax,jMin,jMax,k,
332       I        gTnm1, gSnm1,       I        gT, gS,
333       U        phiHyd,       U        phiHyd,
334       I        myThid )       I        myThid )
335           ELSE           ELSE
# Line 671  C        distinguishe between Stagger an Line 343  C        distinguishe between Stagger an
343  C--      Calculate accelerations in the momentum equations (gU, gV, ...)  C--      Calculate accelerations in the momentum equations (gU, gV, ...)
344  C        and step forward storing the result in gUnm1, gVnm1, etc...  C        and step forward storing the result in gUnm1, gVnm1, etc...
345           IF ( momStepping ) THEN           IF ( momStepping ) THEN
346             CALL CALC_MOM_RHS(  #ifndef DISABLE_MOM_FLUXFORM
347               IF (.NOT. vectorInvariantMomentum) CALL MOM_FLUXFORM(
348       I         bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,
349       I         phiHyd,KappaRU,KappaRV,       I         phiHyd,KappaRU,KappaRV,
350       U         fVerU, fVerV,       U         fVerU, fVerV,
351       I         myTime, myThid)       I         myTime, myIter, myThid)
352    #endif
353    #ifndef DISABLE_MOM_VECINV
354               IF (vectorInvariantMomentum) CALL MOM_VECINV(
355         I         bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,
356         I         phiHyd,KappaRU,KappaRV,
357         U         fVerU, fVerV,
358         I         myTime, myIter, myThid)
359    #endif
360             CALL TIMESTEP(             CALL TIMESTEP(
361       I         bi,bj,iMin,iMax,jMin,jMax,k,       I         bi,bj,iMin,iMax,jMin,jMax,k,
362       I         phiHyd, phiSurfX, phiSurfY,       I         phiHyd, phiSurfX, phiSurfY,
# Line 772  Cjmc(end) Line 453  Cjmc(end)
453          IF (taveFreq.GT.0.) THEN          IF (taveFreq.GT.0.) THEN
454            CALL TIMEAVE_CUMUL_1T(phiHydtave, phiHyd, Nr,            CALL TIMEAVE_CUMUL_1T(phiHydtave, phiHyd, Nr,
455       I                              deltaTclock, bi, bj, myThid)       I                              deltaTclock, bi, bj, myThid)
           IF (ivdc_kappa.NE.0.) THEN  
             CALL TIMEAVE_CUMULATE(ConvectCountTave, ConvectCount, Nr,  
      I                              deltaTclock, bi, bj, myThid)  
           ENDIF  
456          ENDIF          ENDIF
457  #endif /* ALLOW_TIMEAVE */  #endif /* ALLOW_TIMEAVE */
458    
459         ENDDO         ENDDO
460        ENDDO        ENDDO
461    
462  #ifndef EXCLUDE_DEBUGMODE  #ifndef DISABLE_DEBUGMODE
463        If (debugMode) THEN        If (debugMode) THEN
464         CALL DEBUG_STATS_RL(1,EtaN,'EtaN (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(1,EtaN,'EtaN (DYNAMICS)',myThid)
465           CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (DYNAMICS)',myThid)
466         CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (DYNAMICS)',myThid)
467         CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (DYNAMICS)',myThid)
468         CALL DEBUG_STATS_RL(Nr,theta,'Theta (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,theta,'Theta (DYNAMICS)',myThid)
Legend:
Removed from v.1.72  
changed lines
  Added in v.1.83.4.3
  ViewVC Help
Powered by ViewVC 1.1.22