/[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.76 by heimbach, Mon Aug 13 18:05:26 2001 UTC
# Line 29  C     == Global variables === Line 29  C     == Global variables ===
29  #include "PARAMS.h"  #include "PARAMS.h"
30  #include "DYNVARS.h"  #include "DYNVARS.h"
31  #include "GRID.h"  #include "GRID.h"
32    #ifdef ALLOW_PASSIVE_TRACER
33    #include "TR1.h"
34    #endif
35    
36  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
37  # include "tamc.h"  # include "tamc.h"
38  # include "tamc_keys.h"  # include "tamc_keys.h"
39    # include "FFIELDS.h"
40    # ifdef ALLOW_KPP
41    #  include "KPP.h"
42    # endif
43    # ifdef ALLOW_GMREDI
44    #  include "GMREDI.h"
45    # endif
46  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
47    
 #ifdef ALLOW_KPP  
 # include "KPP.h"  
 #endif  
   
48  #ifdef ALLOW_TIMEAVE  #ifdef ALLOW_TIMEAVE
49  #include "TIMEAVE_STATV.h"  #include "TIMEAVE_STATV.h"
50  #endif  #endif
# Line 52  C     myThid - Thread number for this in Line 58  C     myThid - Thread number for this in
58        INTEGER myThid        INTEGER myThid
59    
60  C     == Local variables  C     == Local variables
61  C     xA, yA                 - Per block temporaries holding face areas  C     maskUp                   o maskUp: land/water mask for W points
 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  
62  C     fVer[STUV]               o fVer: Vertical flux term - note fVer  C     fVer[STUV]               o fVer: Vertical flux term - note fVer
63  C                                      is "pipelined" in the vertical  C                                      is "pipelined" in the vertical
64  C                                      so we need an fVer for each  C                                      so we need an fVer for each
# Line 72  C                      In p coords phiHy Line 71  C                      In p coords phiHy
71  C                      surface height anomaly.  C                      surface height anomaly.
72  C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)  C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)
73  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).  
74  C     iMin, iMax     - Ranges and sub-block indices on which calculations  C     iMin, iMax     - Ranges and sub-block indices on which calculations
75  C     jMin, jMax       are applied.  C     jMin, jMax       are applied.
76  C     bi, bj  C     bi, bj
77  C     k, kup,        - Index for layer above and below. kup and kDown  C     k, kup,        - Index for layer above and below. kup and kDown
78  C     kDown, km1       are switched with layer to be the appropriate  C     kDown, km1       are switched with layer to be the appropriate
79  C                      index into fVerTerm.  C                      index into fVerTerm.
80        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  C     tauAB - Adams-Bashforth timestepping weight: 0=forward ; 1/2=Adams-Bashf.
       _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)  
81        _RS maskUp  (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)  
82        _RL fVerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
83        _RL fVerV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
84        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
# Line 96  C                      index into fVerTe Line 86  C                      index into fVerTe
86        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
87        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
88        _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)  
89        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
90        _RL KappaRV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
91        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
92        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
93        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
94          _RL tauAB
95    
96  C This is currently used by IVDC and Diagnostics  C This is currently used by IVDC and Diagnostics
97        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
# Line 163  C         (1 + dt * K * d_zz) theta[n] = Line 152  C         (1 + dt * K * d_zz) theta[n] =
152  C         (1 + dt * K * d_zz) salt[n] = salt*  C         (1 + dt * K * d_zz) salt[n] = salt*
153  C---  C---
154    
 #ifdef ALLOW_AUTODIFF_TAMC  
 C--   dummy statement to end declaration part  
       ikey = 1  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
155  C--   Set up work arrays with valid (i.e. not NaN) values  C--   Set up work arrays with valid (i.e. not NaN) values
156  C     These inital values do not alter the numerical results. They  C     These inital values do not alter the numerical results. They
157  C     just ensure that all memory references are to valid floating  C     just ensure that all memory references are to valid floating
# Line 175  C     point numbers. This prevents spuri Line 159  C     point numbers. This prevents spuri
159  C     uninitialised but inert locations.  C     uninitialised but inert locations.
160        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
161         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  
162          DO k=1,Nr          DO k=1,Nr
163           phiHyd(i,j,k)  = 0. _d 0           phiHyd(i,j,k)  = 0. _d 0
164           KappaRU(i,j,k) = 0. _d 0           KappaRU(i,j,k) = 0. _d 0
# Line 189  C     uninitialised but inert locations. Line 169  C     uninitialised but inert locations.
169          ENDDO          ENDDO
170          rhoKM1 (i,j) = 0. _d 0          rhoKM1 (i,j) = 0. _d 0
171          rhok   (i,j) = 0. _d 0          rhok   (i,j) = 0. _d 0
         maskC  (i,j) = 0. _d 0  
172          phiSurfX(i,j) = 0. _d 0          phiSurfX(i,j) = 0. _d 0
173          phiSurfY(i,j) = 0. _d 0          phiSurfY(i,j) = 0. _d 0
174         ENDDO         ENDDO
175        ENDDO        ENDDO
176    
   
177  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
178  C--   HPF directive to help TAMC  C--   HPF directive to help TAMC
179  CHPF$ INDEPENDENT  CHPF$ INDEPENDENT
# Line 205  CHPF$ INDEPENDENT Line 183  CHPF$ INDEPENDENT
183    
184  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
185  C--    HPF directive to help TAMC  C--    HPF directive to help TAMC
186  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS,fVerU,fVerV  CHPF$  INDEPENDENT, NEW (fVerU,fVerV
187  CHPF$&                  ,phiHyd,utrans,vtrans,maskc,xA,yA  CHPF$&                  ,phiHyd
188  CHPF$&                  ,KappaRT,KappaRS,KappaRU,KappaRV  CHPF$&                  ,KappaRU,KappaRV
189  CHPF$&                  )  CHPF$&                  )
190  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
191    
# Line 233  CHPF$&                  ) Line 211  CHPF$&                  )
211  C--     Set up work arrays that need valid initial values  C--     Set up work arrays that need valid initial values
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
214            rTrans(i,j)   = 0. _d 0            fVerU  (i,j,1) = 0. _d 0
215            fVerT (i,j,1) = 0. _d 0            fVerU  (i,j,2) = 0. _d 0
216            fVerT (i,j,2) = 0. _d 0            fVerV  (i,j,1) = 0. _d 0
217            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  
218           ENDDO           ENDDO
219          ENDDO          ENDDO
220    
         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  
 #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)  
         ENDDO  
   
 #ifdef  ALLOW_OBCS  
 C--     Calculate future values on open boundaries  
         IF (useOBCS) THEN  
           CALL OBCS_CALC( bi, bj, myTime+deltaT,  
      I            uVel, vVel, wVel, theta, salt,  
      I            myThid )  
         ENDIF  
 #endif  /* ALLOW_OBCS */  
   
 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  
 #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  
 #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  
           DO j=1-OLy,sNy+OLy  
             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 */  
         ENDIF  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE KPPghat   (:,:,:,bi,bj)  
 CADJ &   , KPPviscAz (:,:,:,bi,bj)  
 CADJ &   , KPPdiffKzT(:,:,:,bi,bj)  
 CADJ &   , KPPdiffKzS(:,:,:,bi,bj)  
 CADJ &   , KPPfrac   (:,:  ,bi,bj)  
 CADJ &                 = comlev1_bibj, key=ikey, byte=isbyte  
 #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  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 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  INCLUDE_CALC_DIFFUSIVITY_CALL  
 C--      Calculate the total vertical diffusivity  
          CALL CALC_DIFFUSIVITY(  
      I        bi,bj,iMin,iMax,jMin,jMax,k,  
      I        maskC,maskup,  
      O        KappaRT,KappaRS,KappaRU,KappaRV,  
      I        myThid)  
 #endif  
   
 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  
   
221  C--     Start computation of dynamics  C--     Start computation of dynamics
222          iMin = 1-OLx+2          iMin = 1-OLx+2
223          iMax = sNx+OLx-1          iMax = sNx+OLx-1
224          jMin = 1-OLy+2          jMin = 1-OLy+2
225          jMax = sNy+OLy-1          jMax = sNy+OLy-1
226    
227    #ifdef ALLOW_AUTODIFF_TAMC
228    CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
229    CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
230    CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
231    #endif /* ALLOW_AUTODIFF_TAMC */
232    
233  C--     Explicit part of the Surface Potentiel Gradient (add in TIMESTEP)  C--     Explicit part of the Surface Potentiel Gradient (add in TIMESTEP)
234  C       (note: this loop will be replaced by CALL CALC_GRAD_ETA)  C       (note: this loop will be replaced by CALL CALC_GRAD_ETA)
235          IF (implicSurfPress.NE.1.) THEN          IF (implicSurfPress.NE.1.) THEN
# Line 585  C--       kDown  Cycles through 2,1 to p Line 251  C--       kDown  Cycles through 2,1 to p
251            kup  = 1+MOD(k+1,2)            kup  = 1+MOD(k+1,2)
252            kDown= 1+MOD(k,2)            kDown= 1+MOD(k,2)
253    
254    #ifdef ALLOW_AUTODIFF_TAMC
255             kkey = (ikey-1)*Nr + k
256    #endif /* ALLOW_AUTODIFF_TAMC */
257    
258  C--      Integrate hydrostatic balance for phiHyd with BC of  C--      Integrate hydrostatic balance for phiHyd with BC of
259  C        phiHyd(z=0)=0  C        phiHyd(z=0)=0
260  C        distinguishe between Stagger and Non Stagger time stepping  C        distinguishe between Stagger and Non Stagger time stepping
# Line 602  C        distinguishe between Stagger an Line 272  C        distinguishe between Stagger an
272       I        myThid )       I        myThid )
273           ENDIF           ENDIF
274    
275    #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL
276    C--      Calculate the total vertical diffusivity
277             CALL CALC_VISCOSITY(
278         I        bi,bj,iMin,iMax,jMin,jMax,k,
279         I        maskUp,
280         O        KappaRU,KappaRV,
281         I        myThid)
282    #endif
283    
284  C--      Calculate accelerations in the momentum equations (gU, gV, ...)  C--      Calculate accelerations in the momentum equations (gU, gV, ...)
285  C        and step forward storing the result in gUnm1, gVnm1, etc...  C        and step forward storing the result in gUnm1, gVnm1, etc...
286           IF ( momStepping ) THEN           IF ( momStepping ) THEN
# Line 704  Cjmc(end) Line 383  Cjmc(end)
383    
384  #ifdef ALLOW_TIMEAVE  #ifdef ALLOW_TIMEAVE
385          IF (taveFreq.GT.0.) THEN          IF (taveFreq.GT.0.) THEN
386            CALL TIMEAVE_CUMULATE(phiHydtave, phiHyd, Nr,            CALL TIMEAVE_CUMUL_1T(phiHydtave, phiHyd, Nr,
387       I                              deltaTclock, bi, bj, myThid)       I                              deltaTclock, bi, bj, myThid)
388            IF (ivdc_kappa.NE.0.) THEN            IF (ivdc_kappa.NE.0.) THEN
389              CALL TIMEAVE_CUMULATE(ConvectCountTave, ConvectCount, Nr,              CALL TIMEAVE_CUMULATE(ConvectCountTave, ConvectCount, Nr,
# Line 716  Cjmc(end) Line 395  Cjmc(end)
395         ENDDO         ENDDO
396        ENDDO        ENDDO
397    
398    #ifndef EXCLUDE_DEBUGMODE
399          If (debugMode) THEN
400           CALL DEBUG_STATS_RL(1,EtaN,'EtaN (DYNAMICS)',myThid)
401           CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (DYNAMICS)',myThid)
402           CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (DYNAMICS)',myThid)
403           CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (DYNAMICS)',myThid)
404           CALL DEBUG_STATS_RL(Nr,theta,'Theta (DYNAMICS)',myThid)
405           CALL DEBUG_STATS_RL(Nr,salt,'Salt (DYNAMICS)',myThid)
406           CALL DEBUG_STATS_RL(Nr,Gu,'Gu (DYNAMICS)',myThid)
407           CALL DEBUG_STATS_RL(Nr,Gv,'Gv (DYNAMICS)',myThid)
408           CALL DEBUG_STATS_RL(Nr,Gt,'Gt (DYNAMICS)',myThid)
409           CALL DEBUG_STATS_RL(Nr,Gs,'Gs (DYNAMICS)',myThid)
410           CALL DEBUG_STATS_RL(Nr,GuNm1,'GuNm1 (DYNAMICS)',myThid)
411           CALL DEBUG_STATS_RL(Nr,GvNm1,'GvNm1 (DYNAMICS)',myThid)
412           CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (DYNAMICS)',myThid)
413           CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (DYNAMICS)',myThid)
414          ENDIF
415    #endif
416    
417        RETURN        RETURN
418        END        END
Legend:
Removed from v.1.66  
changed lines
  Added in v.1.76
  ViewVC Help
Powered by ViewVC 1.1.22