/[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.67 by heimbach, Mon May 14 21:46:17 2001 UTC revision 1.83 by heimbach, Thu Sep 27 20:12:10 2001 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"
# Line 41  C     == Global variables === Line 87  C     == Global variables ===
87  #  include "GMREDI.h"  #  include "GMREDI.h"
88  # endif  # endif
89  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
   
90  #ifdef ALLOW_TIMEAVE  #ifdef ALLOW_TIMEAVE
91  #include "TIMEAVE_STATV.h"  #include "TIMEAVE_STATV.h"
92  #endif  #endif
93    
94    C     !CALLING SEQUENCE:
95    C     DYNAMICS()
96    C      |
97    C      |-- CALC_GRAD_PHI_SURF
98    C      |
99    C      |-- CALC_VISCOSITY
100    C      |
101    C      |-- CALC_PHI_HYD  
102    C      |
103    C      |-- MOM_FLUXFORM  
104    C      |
105    C      |-- MOM_VECINV    
106    C      |
107    C      |-- TIMESTEP      
108    C      |
109    C      |-- OBCS_APPLY_UV
110    C      |
111    C      |-- IMPLDIFF      
112    C      |
113    C      |-- OBCS_APPLY_UV
114    C      |
115    C      |-- CALL TIMEAVE_CUMUL_1T
116    C      |-- CALL TIMEAVE_CUMULATE
117    C      |-- CALL DEBUG_STATS_RL
118    
119    C     !INPUT/OUTPUT PARAMETERS:
120  C     == Routine arguments ==  C     == Routine arguments ==
121  C     myTime - Current time in simulation  C     myTime - Current time in simulation
122  C     myIter - Current iteration number in simulation  C     myIter - Current iteration number in simulation
# Line 54  C     myThid - Thread number for this in Line 125  C     myThid - Thread number for this in
125        INTEGER myIter        INTEGER myIter
126        INTEGER myThid        INTEGER myThid
127    
128    C     !LOCAL VARIABLES:
129  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  
130  C     fVer[STUV]               o fVer: Vertical flux term - note fVer  C     fVer[STUV]               o fVer: Vertical flux term - note fVer
131  C                                      is "pipelined" in the vertical  C                                      is "pipelined" in the vertical
132  C                                      so we need an fVer for each  C                                      so we need an fVer for each
# Line 75  C                      In p coords phiHy Line 139  C                      In p coords phiHy
139  C                      surface height anomaly.  C                      surface height anomaly.
140  C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)  C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)
141  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).  
142  C     iMin, iMax     - Ranges and sub-block indices on which calculations  C     iMin, iMax     - Ranges and sub-block indices on which calculations
143  C     jMin, jMax       are applied.  C     jMin, jMax       are applied.
144  C     bi, bj  C     bi, bj
145  C     k, kup,        - Index for layer above and below. kup and kDown  C     k, kup,        - Index for layer above and below. kup and kDown
146  C     kDown, km1       are switched with layer to be the appropriate  C     kDown, km1       are switched with layer to be the appropriate
147  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)  
148        _RL fVerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
149        _RL fVerV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
150        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
# Line 99  C                      index into fVerTe Line 152  C                      index into fVerTe
152        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
153        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
154        _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)  
155        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
156        _RL KappaRV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
157        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
# Line 114  C This is currently used by IVDC and Dia Line 165  C This is currently used by IVDC and Dia
165        INTEGER jMin, jMax        INTEGER jMin, jMax
166        INTEGER bi, bj        INTEGER bi, bj
167        INTEGER i, j        INTEGER i, j
168        INTEGER k, km1, kup, kDown        INTEGER k, km1, kp1, kup, kDown
169    
170  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
171  c     CHARACTER*(MAX_LEN_MBUF) suff  c     CHARACTER*(MAX_LEN_MBUF) suff
# Line 165  C         salt* = salt[n] + dt x ( 3/2 G Line 216  C         salt* = salt[n] + dt x ( 3/2 G
216  C         (1 + dt * K * d_zz) theta[n] = theta*  C         (1 + dt * K * d_zz) theta[n] = theta*
217  C         (1 + dt * K * d_zz) salt[n] = salt*  C         (1 + dt * K * d_zz) salt[n] = salt*
218  C---  C---
219    CEOP
 #ifdef ALLOW_AUTODIFF_TAMC  
 C--   dummy statement to end declaration part  
       ikey = 1  
 #endif /* ALLOW_AUTODIFF_TAMC */  
220    
221  C--   Set up work arrays with valid (i.e. not NaN) values  C--   Set up work arrays with valid (i.e. not NaN) values
222  C     These inital values do not alter the numerical results. They  C     These inital values do not alter the numerical results. They
# Line 178  C     point numbers. This prevents spuri Line 225  C     point numbers. This prevents spuri
225  C     uninitialised but inert locations.  C     uninitialised but inert locations.
226        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
227         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  
228          DO k=1,Nr          DO k=1,Nr
229           phiHyd(i,j,k)  = 0. _d 0           phiHyd(i,j,k)  = 0. _d 0
230           KappaRU(i,j,k) = 0. _d 0           KappaRU(i,j,k) = 0. _d 0
# Line 192  C     uninitialised but inert locations. Line 235  C     uninitialised but inert locations.
235          ENDDO          ENDDO
236          rhoKM1 (i,j) = 0. _d 0          rhoKM1 (i,j) = 0. _d 0
237          rhok   (i,j) = 0. _d 0          rhok   (i,j) = 0. _d 0
         maskC  (i,j) = 0. _d 0  
238          phiSurfX(i,j) = 0. _d 0          phiSurfX(i,j) = 0. _d 0
239          phiSurfY(i,j) = 0. _d 0          phiSurfY(i,j) = 0. _d 0
240         ENDDO         ENDDO
241        ENDDO        ENDDO
242    
   
243  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
244  C--   HPF directive to help TAMC  C--   HPF directive to help TAMC
245  CHPF$ INDEPENDENT  CHPF$ INDEPENDENT
# Line 208  CHPF$ INDEPENDENT Line 249  CHPF$ INDEPENDENT
249    
250  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
251  C--    HPF directive to help TAMC  C--    HPF directive to help TAMC
252  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS,fVerU,fVerV  CHPF$  INDEPENDENT, NEW (fVerU,fVerV
253  CHPF$&                  ,phiHyd,utrans,vtrans,maskc,xA,yA  CHPF$&                  ,phiHyd
254  CHPF$&                  ,KappaRT,KappaRS,KappaRU,KappaRV  CHPF$&                  ,KappaRU,KappaRV
255  CHPF$&                  )  CHPF$&                  )
256  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
257    
# Line 219  CHPF$&                  ) Line 260  CHPF$&                  )
260  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
261            act1 = bi - myBxLo(myThid)            act1 = bi - myBxLo(myThid)
262            max1 = myBxHi(myThid) - myBxLo(myThid) + 1            max1 = myBxHi(myThid) - myBxLo(myThid) + 1
   
263            act2 = bj - myByLo(myThid)            act2 = bj - myByLo(myThid)
264            max2 = myByHi(myThid) - myByLo(myThid) + 1            max2 = myByHi(myThid) - myByLo(myThid) + 1
   
265            act3 = myThid - 1            act3 = myThid - 1
266            max3 = nTx*nTy            max3 = nTx*nTy
   
267            act4 = ikey_dynamics - 1            act4 = ikey_dynamics - 1
   
268            ikey = (act1 + 1) + act2*max1            ikey = (act1 + 1) + act2*max1
269       &                      + act3*max1*max2       &                      + act3*max1*max2
270       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
# Line 236  CHPF$&                  ) Line 273  CHPF$&                  )
273  C--     Set up work arrays that need valid initial values  C--     Set up work arrays that need valid initial values
274          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
275           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
276            rTrans(i,j)   = 0. _d 0            fVerU  (i,j,1) = 0. _d 0
277            fVerT (i,j,1) = 0. _d 0            fVerU  (i,j,2) = 0. _d 0
278            fVerT (i,j,2) = 0. _d 0            fVerV  (i,j,1) = 0. _d 0
279            fVerS (i,j,1) = 0. _d 0            fVerV  (i,j,2) = 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  
   
         DO k=1,Nr  
          DO j=1-OLy,sNy+OLy  
           DO i=1-OLx,sNx+OLx  
 C This is currently also used by IVDC and Diagnostics  
            ConvectCount(i,j,k) = 0.  
            KappaRT(i,j,k) = 0. _d 0  
            KappaRS(i,j,k) = 0. _d 0  
           ENDDO  
280           ENDDO           ENDDO
281          ENDDO          ENDDO
282    
283          iMin = 1-OLx+1  C--     Start computation of dynamics
284          iMax = sNx+OLx          iMin = 1-OLx+2
285          jMin = 1-OLy+1          iMax = sNx+OLx-1
286          jMax = sNy+OLy          jMin = 1-OLy+2
287            jMax = sNy+OLy-1
   
 #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  
 #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  
288    
289  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
 cph avoids recomputation of integrate_for_w  
290  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
291  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
292    
293  #ifdef  ALLOW_OBCS  C--     Explicit part of the Surface Potentiel Gradient (add in TIMESTEP)
294  C--     Calculate future values on open boundaries  C       (note: this loop will be replaced by CALL CALC_GRAD_ETA)
295          IF (useOBCS) THEN          IF (implicSurfPress.NE.1.) THEN
296            CALL OBCS_CALC( bi, bj, myTime+deltaT,            CALL CALC_GRAD_PHI_SURF(
297       I            uVel, vVel, wVel, theta, salt,       I         bi,bj,iMin,iMax,jMin,jMax,
298       I            myThid )       I         etaN,
299          ENDIF       O         phiSurfX,phiSurfY,
300  #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 */  
301          ENDIF          ENDIF
302    
303  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
304  CADJ STORE KPPghat   (:,:,:,bi,bj)  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
305  CADJ &   , KPPviscAz (:,:,:,bi,bj)  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
306  CADJ &   , KPPdiffKzT(:,:,:,bi,bj)  #ifdef ALLOW_KPP
307  CADJ &   , KPPdiffKzS(:,:,:,bi,bj)  CADJ STORE KPPviscAz (:,:,:,bi,bj)
 CADJ &   , KPPfrac   (:,:  ,bi,bj)  
308  CADJ &                 = comlev1_bibj, key=ikey, byte=isbyte  CADJ &                 = comlev1_bibj, key=ikey, byte=isbyte
309  #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  
 #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  
 #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,km1,kup,kDown,  
      O        xA,yA,uTrans,vTrans,rTrans,maskC,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  
310  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
311    
312  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL
313  C--      Calculate the total vertical diffusivity  C--      Calculate the total vertical diffusivity
314           CALL CALC_DIFFUSIVITY(          DO k=1,Nr
315             CALL CALC_VISCOSITY(
316       I        bi,bj,iMin,iMax,jMin,jMax,k,       I        bi,bj,iMin,iMax,jMin,jMax,k,
317       I        maskC,maskup,       O        KappaRU,KappaRV,
      O        KappaRT,KappaRS,KappaRU,KappaRV,  
318       I        myThid)       I        myThid)
319           ENDDO
320  #endif  #endif
321    
 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  
   
322  C--     Start of dynamics loop  C--     Start of dynamics loop
323          DO k=1,Nr          DO k=1,Nr
324    
# Line 616  C--       kup    Cycles through 1,2 to p Line 327  C--       kup    Cycles through 1,2 to p
327  C--       kDown  Cycles through 2,1 to point to current layer  C--       kDown  Cycles through 2,1 to point to current layer
328    
329            km1  = MAX(1,k-1)            km1  = MAX(1,k-1)
330              kp1  = MIN(k+1,Nr)
331            kup  = 1+MOD(k+1,2)            kup  = 1+MOD(k+1,2)
332            kDown= 1+MOD(k,2)            kDown= 1+MOD(k,2)
333    
334    #ifdef ALLOW_AUTODIFF_TAMC
335             kkey = (ikey-1)*Nr + k
336    #endif /* ALLOW_AUTODIFF_TAMC */
337    
338  C--      Integrate hydrostatic balance for phiHyd with BC of  C--      Integrate hydrostatic balance for phiHyd with BC of
339  C        phiHyd(z=0)=0  C        phiHyd(z=0)=0
340  C        distinguishe between Stagger and Non Stagger time stepping  C        distinguishe between Stagger and Non Stagger time stepping
341           IF (staggerTimeStep) THEN           IF (staggerTimeStep) THEN
342             CALL CALC_PHI_HYD(             CALL CALC_PHI_HYD(
343       I        bi,bj,iMin,iMax,jMin,jMax,k,       I        bi,bj,iMin,iMax,jMin,jMax,k,
344       I        gTnm1, gSnm1,       I        gT, gS,
345       U        phiHyd,       U        phiHyd,
346       I        myThid )       I        myThid )
347           ELSE           ELSE
# Line 639  C        distinguishe between Stagger an Line 355  C        distinguishe between Stagger an
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,
361         I         phiHyd,KappaRU,KappaRV,
362         U         fVerU, fVerV,
363         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,       I         bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,
368       I         phiHyd,KappaRU,KappaRV,       I         phiHyd,KappaRU,KappaRV,
369       U         fVerU, fVerV,       U         fVerU, fVerV,
370       I         myTime, myThid)       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         phiHyd, phiSurfX, phiSurfY,
# Line 738  Cjmc(end) Line 463  Cjmc(end)
463    
464  #ifdef ALLOW_TIMEAVE  #ifdef ALLOW_TIMEAVE
465          IF (taveFreq.GT.0.) THEN          IF (taveFreq.GT.0.) THEN
466            CALL TIMEAVE_CUMULATE(phiHydtave, phiHyd, Nr,            CALL TIMEAVE_CUMUL_1T(phiHydtave, phiHyd, Nr,
467       I                              deltaTclock, bi, bj, myThid)       I                              deltaTclock, bi, bj, myThid)
468            IF (ivdc_kappa.NE.0.) THEN            IF (ivdc_kappa.NE.0.) THEN
469              CALL TIMEAVE_CUMULATE(ConvectCountTave, ConvectCount, Nr,              CALL TIMEAVE_CUMULATE(ConvectCountTave, ConvectCount, Nr,
# Line 750  Cjmc(end) Line 475  Cjmc(end)
475         ENDDO         ENDDO
476        ENDDO        ENDDO
477    
478    #ifndef DISABLE_DEBUGMODE
479          If (debugMode) THEN
480           CALL DEBUG_STATS_RL(1,EtaN,'EtaN (DYNAMICS)',myThid)
481           CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (DYNAMICS)',myThid)
482           CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (DYNAMICS)',myThid)
483           CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (DYNAMICS)',myThid)
484           CALL DEBUG_STATS_RL(Nr,theta,'Theta (DYNAMICS)',myThid)
485           CALL DEBUG_STATS_RL(Nr,salt,'Salt (DYNAMICS)',myThid)
486           CALL DEBUG_STATS_RL(Nr,Gu,'Gu (DYNAMICS)',myThid)
487           CALL DEBUG_STATS_RL(Nr,Gv,'Gv (DYNAMICS)',myThid)
488           CALL DEBUG_STATS_RL(Nr,Gt,'Gt (DYNAMICS)',myThid)
489           CALL DEBUG_STATS_RL(Nr,Gs,'Gs (DYNAMICS)',myThid)
490           CALL DEBUG_STATS_RL(Nr,GuNm1,'GuNm1 (DYNAMICS)',myThid)
491           CALL DEBUG_STATS_RL(Nr,GvNm1,'GvNm1 (DYNAMICS)',myThid)
492           CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (DYNAMICS)',myThid)
493           CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (DYNAMICS)',myThid)
494          ENDIF
495    #endif
496    
497        RETURN        RETURN
498        END        END
Legend:
Removed from v.1.67  
changed lines
  Added in v.1.83
  ViewVC Help
Powered by ViewVC 1.1.22