/[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.54.2.8 by jmc, Fri Jan 12 14:39:53 2001 UTC revision 1.68 by adcroft, Tue May 29 14:01:37 2001 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
5    
# Line 25  C     \================================= Line 26  C     \=================================
26  C     == Global variables ===  C     == Global variables ===
27  #include "SIZE.h"  #include "SIZE.h"
28  #include "EEPARAMS.h"  #include "EEPARAMS.h"
 #include "CG2D.h"  
29  #include "PARAMS.h"  #include "PARAMS.h"
30  #include "DYNVARS.h"  #include "DYNVARS.h"
31  #include "GRID.h"  #include "GRID.h"
# Line 33  C     == Global variables === Line 33  C     == Global variables ===
33  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
34  # include "tamc.h"  # include "tamc.h"
35  # include "tamc_keys.h"  # include "tamc_keys.h"
36    # include "FFIELDS.h"
37    # ifdef ALLOW_KPP
38    #  include "KPP.h"
39    # endif
40    # ifdef ALLOW_GMREDI
41    #  include "GMREDI.h"
42    # endif
43  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
44    
45  #ifdef ALLOW_KPP  #ifdef ALLOW_TIMEAVE
46  # include "KPP.h"  #include "TIMEAVE_STATV.h"
47  #endif  #endif
48    
49  C     == Routine arguments ==  C     == Routine arguments ==
# Line 51  C     == Local variables Line 58  C     == Local variables
58  C     xA, yA                 - Per block temporaries holding face areas  C     xA, yA                 - Per block temporaries holding face areas
59  C     uTrans, vTrans, rTrans - Per block temporaries holding flow  C     uTrans, vTrans, rTrans - Per block temporaries holding flow
60  C                              transport  C                              transport
61  C     rVel                     o uTrans: Zonal transport  C                              o uTrans: Zonal transport
62  C                              o vTrans: Meridional transport  C                              o vTrans: Meridional transport
63  C                              o rTrans: Vertical transport  C                              o rTrans: Vertical transport
64  C                              o rVel:   Vertical velocity at upper and  C     maskUp                   o maskUp: land/water mask for W points
 C                                        lower cell faces.  
 C     maskC,maskUp             o maskC: land/water mask for tracer cells  
 C                              o maskUp: land/water mask for W points  
65  C     fVer[STUV]               o fVer: Vertical flux term - note fVer  C     fVer[STUV]               o fVer: Vertical flux term - note fVer
66  C                                      is "pipelined" in the vertical  C                                      is "pipelined" in the vertical
67  C                                      so we need an fVer for each  C                                      so we need an fVer for each
68  C                                      variable.  C                                      variable.
69  C     rhoK, rhoKM1   - Density at current level, and level above  C     rhoK, rhoKM1   - Density at current level, and level above
 C     buoyK, buoyKM1 - Buoyancy at current level and level above.  
70  C     phiHyd         - Hydrostatic part of the potential phiHydi.  C     phiHyd         - Hydrostatic part of the potential phiHydi.
71  C                      In z coords phiHydiHyd is the hydrostatic  C                      In z coords phiHydiHyd is the hydrostatic
72  C                      pressure anomaly  C                      Potential (=pressure/rho0) anomaly
73  C                      In p coords phiHydiHyd is the geopotential  C                      In p coords phiHydiHyd is the geopotential
74  C                      surface height  C                      surface height anomaly.
75  C                      anomaly.  C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)
76  C     etaSurfX,      - Holds surface elevation gradient in X and Y.  C     phiSurfY             or geopotentiel (atmos) in X and Y direction
 C     etaSurfY  
77  C     KappaRT,       - Total diffusion in vertical for T and S.  C     KappaRT,       - Total diffusion in vertical for T and S.
78  C     KappaRS          (background + spatially varying, isopycnal term).  C     KappaRS          (background + spatially varying, isopycnal term).
79  C     iMin, iMax     - Ranges and sub-block indices on which calculations  C     iMin, iMax     - Ranges and sub-block indices on which calculations
# Line 80  C     bi, bj Line 82  C     bi, bj
82  C     k, kup,        - Index for layer above and below. kup and kDown  C     k, kup,        - Index for layer above and below. kup and kDown
83  C     kDown, km1       are switched with layer to be the appropriate  C     kDown, km1       are switched with layer to be the appropriate
84  C                      index into fVerTerm.  C                      index into fVerTerm.
85    C     tauAB - Adams-Bashforth timestepping weight: 0=forward ; 1/2=Adams-Bashf.
86        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
87        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
88        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
89        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
90        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RL rVel    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)  
       _RS maskC   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
91        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
92        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
93        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
# Line 95  C                      index into fVerTe Line 96  C                      index into fVerTe
96        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
97        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
99        _RL buoyKM1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
100        _RL buoyK   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
101        _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
102        _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
103        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
# Line 104  C                      index into fVerTe Line 105  C                      index into fVerTe
105        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
106        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
107        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
108          _RL tauAB
109    
110  C This is currently also used by IVDC and Diagnostics  C This is currently used by IVDC and Diagnostics
 C #ifdef INCLUDE_CONVECT_CALL  
111        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
 C #endif  
112    
113        INTEGER iMin, iMax        INTEGER iMin, iMax
114        INTEGER jMin, jMax        INTEGER jMin, jMax
# Line 116  C #endif Line 116  C #endif
116        INTEGER i, j        INTEGER i, j
117        INTEGER k, km1, kup, kDown        INTEGER k, km1, kup, kDown
118    
119  #ifdef ALLOW_AUTODIFF_TAMC  Cjmc : add for phiHyd output <- but not working if multi tile per CPU
120        INTEGER    isbyte  c     CHARACTER*(MAX_LEN_MBUF) suff
121        PARAMETER( isbyte = 4 )  c     LOGICAL  DIFFERENT_MULTIPLE
122    c     EXTERNAL DIFFERENT_MULTIPLE
123        INTEGER act1, act2, act3, act4  Cjmc(end)
124        INTEGER max1, max2, max3  
       INTEGER iikey, kkey  
       INTEGER maximpl  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
125  C---    The algorithm...  C---    The algorithm...
126  C  C
127  C       "Correction Step"  C       "Correction Step"
# Line 140  C       "Calculation of Gs" Line 136  C       "Calculation of Gs"
136  C       ===================  C       ===================
137  C       This is where all the accelerations and tendencies (ie.  C       This is where all the accelerations and tendencies (ie.
138  C       physics, parameterizations etc...) are calculated  C       physics, parameterizations etc...) are calculated
 C         rVel = sum_r ( div. u[n] )  
139  C         rho = rho ( theta[n], salt[n] )  C         rho = rho ( theta[n], salt[n] )
140  C         b   = b(rho, theta)  C         b   = b(rho, theta)
141  C         K31 = K31 ( rho )  C         K31 = K31 ( rho )
142  C         Gu[n] = Gu( u[n], v[n], rVel, b, ... )  C         Gu[n] = Gu( u[n], v[n], wVel, b, ... )
143  C         Gv[n] = Gv( u[n], v[n], rVel, b, ... )  C         Gv[n] = Gv( u[n], v[n], wVel, b, ... )
144  C         Gt[n] = Gt( theta[n], u[n], v[n], rVel, K31, ... )  C         Gt[n] = Gt( theta[n], u[n], v[n], wVel, K31, ... )
145  C         Gs[n] = Gs( salt[n], u[n], v[n], rVel, K31, ... )  C         Gs[n] = Gs( salt[n], u[n], v[n], wVel, K31, ... )
146  C  C
147  C       "Time-stepping" or "Prediction"  C       "Time-stepping" or "Prediction"
148  C       ================================  C       ================================
# Line 188  C     uninitialised but inert locations. Line 183  C     uninitialised but inert locations.
183          uTrans(i,j)  = 0. _d 0          uTrans(i,j)  = 0. _d 0
184          vTrans(i,j)  = 0. _d 0          vTrans(i,j)  = 0. _d 0
185          DO k=1,Nr          DO k=1,Nr
186           phiHyd (i,j,k)  = 0. _d 0           phiHyd(i,j,k)  = 0. _d 0
187           KappaRU(i,j,k) = 0. _d 0           KappaRU(i,j,k) = 0. _d 0
188           KappaRV(i,j,k) = 0. _d 0           KappaRV(i,j,k) = 0. _d 0
189           sigmaX(i,j,k) = 0. _d 0           sigmaX(i,j,k) = 0. _d 0
# Line 197  C     uninitialised but inert locations. Line 192  C     uninitialised but inert locations.
192          ENDDO          ENDDO
193          rhoKM1 (i,j) = 0. _d 0          rhoKM1 (i,j) = 0. _d 0
194          rhok   (i,j) = 0. _d 0          rhok   (i,j) = 0. _d 0
195          buoyKM1(i,j) = 0. _d 0          phiSurfX(i,j) = 0. _d 0
196          buoyK  (i,j) = 0. _d 0          phiSurfY(i,j) = 0. _d 0
         maskC  (i,j) = 0. _d 0  
197         ENDDO         ENDDO
198        ENDDO        ENDDO
199    
# Line 213  CHPF$ INDEPENDENT Line 207  CHPF$ INDEPENDENT
207    
208  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
209  C--    HPF directive to help TAMC  C--    HPF directive to help TAMC
210  CHPF$  INDEPENDENT, NEW (rTrans,rVel,fVerT,fVerS,fVerU,fVerV  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS,fVerU,fVerV
211  CHPF$&                  ,phiHyd,utrans,vtrans,maskc,xA,yA  CHPF$&                  ,phiHyd,utrans,vtrans,xA,yA
212  CHPF$&                  ,KappaRT,KappaRS,KappaRU,KappaRV  CHPF$&                  ,KappaRT,KappaRS,KappaRU,KappaRV
213  CHPF$&                  )  CHPF$&                  )
214  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
# Line 242  C--     Set up work arrays that need val Line 236  C--     Set up work arrays that need val
236          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
237           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
238            rTrans(i,j)   = 0. _d 0            rTrans(i,j)   = 0. _d 0
           rVel  (i,j,1) = 0. _d 0  
           rVel  (i,j,2) = 0. _d 0  
239            fVerT (i,j,1) = 0. _d 0            fVerT (i,j,1) = 0. _d 0
240            fVerT (i,j,2) = 0. _d 0            fVerT (i,j,2) = 0. _d 0
241            fVerS (i,j,1) = 0. _d 0            fVerS (i,j,1) = 0. _d 0
# Line 252  C--     Set up work arrays that need val Line 244  C--     Set up work arrays that need val
244            fVerU (i,j,2) = 0. _d 0            fVerU (i,j,2) = 0. _d 0
245            fVerV (i,j,1) = 0. _d 0            fVerV (i,j,1) = 0. _d 0
246            fVerV (i,j,2) = 0. _d 0            fVerV (i,j,2) = 0. _d 0
           phiHyd(i,j,1) = 0. _d 0  
247           ENDDO           ENDDO
248          ENDDO          ENDDO
249    
250          DO k=1,Nr          DO k=1,Nr
251           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
252            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
253  #ifdef INCLUDE_CONVECT_CALL  C This is currently also used by IVDC and Diagnostics
254             ConvectCount(i,j,k) = 0.             ConvectCount(i,j,k) = 0.
 #endif  
255             KappaRT(i,j,k) = 0. _d 0             KappaRT(i,j,k) = 0. _d 0
256             KappaRS(i,j,k) = 0. _d 0             KappaRS(i,j,k) = 0. _d 0
257            ENDDO            ENDDO
# Line 274  C--     Set up work arrays that need val Line 264  C--     Set up work arrays that need val
264          jMax = sNy+OLy          jMax = sNy+OLy
265    
266    
267    #ifdef ALLOW_AUTODIFF_TAMC
268    CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
269    CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
270    CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
271    CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
272    #endif /* ALLOW_AUTODIFF_TAMC */
273    
274  C--     Start of diagnostic loop  C--     Start of diagnostic loop
275          DO k=Nr,1,-1          DO k=Nr,1,-1
276    
277  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
278  C? Patrick, is this formula correct now that we change the loop range?  C? Patrick, is this formula correct now that we change the loop range?
279  C? Do we still need this?  C? Do we still need this?
280           kkey = (ikey-1)*(Nr-2+1) + (k-2) + 1  cph kkey formula corrected.
281    cph Needed for rhok, rhokm1, in the case useGMREDI.
282             kkey = (ikey-1)*Nr + k
283    CADJ STORE rhokm1(:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte
284    CADJ STORE rhok  (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte
285  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
286    
287  C--       Integrate continuity vertically for vertical velocity  C--       Integrate continuity vertically for vertical velocity
# Line 290  C--       Integrate continuity verticall Line 291  C--       Integrate continuity verticall
291       I                         myThid )       I                         myThid )
292    
293  #ifdef    ALLOW_OBCS  #ifdef    ALLOW_OBCS
294  C--       Calculate future values on open boundaries  #ifdef    ALLOW_NONHYDROSTATIC
295            IF (openBoundaries) THEN  C--       Apply OBC to W if in N-H mode
296  #ifdef      ALLOW_NONHYDROSTATIC            IF (useOBCS.AND.nonHydrostatic) THEN
297              IF (nonHydrostatic) THEN              CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )
               CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )  
             ENDIF  
 #endif      /* ALLOW_NONHYDROSTATIC */  
             CALL OBCS_CALC( bi, bj, k, myTime+deltaT, myThid )  
298            ENDIF            ENDIF
299    #endif    /* ALLOW_NONHYDROSTATIC */
300  #endif    /* ALLOW_OBCS */  #endif    /* ALLOW_OBCS */
301    
302  C--       Calculate gradients of potential density for isoneutral  C--       Calculate gradients of potential density for isoneutral
303  C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)  C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
304  c         IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN  c         IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN
305            IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN            IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN
306    #ifdef ALLOW_AUTODIFF_TAMC
307    CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
308    CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
309    #endif /* ALLOW_AUTODIFF_TAMC */
310              CALL FIND_RHO(              CALL FIND_RHO(
311       I        bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,       I        bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,
312         I        theta, salt,
313       O        rhoK,       O        rhoK,
314       I        myThid )       I        myThid )
315              IF (k.GT.1) CALL FIND_RHO(              IF (k.GT.1) THEN
316    #ifdef ALLOW_AUTODIFF_TAMC
317    CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
318    CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
319    #endif /* ALLOW_AUTODIFF_TAMC */
320                 CALL FIND_RHO(
321       I        bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,       I        bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,
322         I        theta, salt,
323       O        rhoKm1,       O        rhoKm1,
324       I        myThid )       I        myThid )
325                ENDIF
326              CALL GRAD_SIGMA(              CALL GRAD_SIGMA(
327       I             bi, bj, iMin, iMax, jMin, jMax, k,       I             bi, bj, iMin, iMax, jMin, jMax, k,
328       I             rhoK, rhoKm1, rhoK,       I             rhoK, rhoKm1, rhoK,
# Line 328  c ==> should use sigmaR !!! Line 338  c ==> should use sigmaR !!!
338       I        rhoKm1, rhoK,       I        rhoKm1, rhoK,
339       U        ConvectCount, KappaRT, KappaRS,       U        ConvectCount, KappaRT, KappaRS,
340       I        myTime, myIter, myThid)       I        myTime, myIter, myThid)
341            END IF            ENDIF
342    
343  C--     end of diagnostic k loop (Nr:1)  C--     end of diagnostic k loop (Nr:1)
344          ENDDO          ENDDO
345    
346    #ifdef ALLOW_AUTODIFF_TAMC
347    cph avoids recomputation of integrate_for_w
348    CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
349    #endif /* ALLOW_AUTODIFF_TAMC */
350    
351    #ifdef  ALLOW_OBCS
352    C--     Calculate future values on open boundaries
353            IF (useOBCS) THEN
354              CALL OBCS_CALC( bi, bj, myTime+deltaT,
355         I            uVel, vVel, wVel, theta, salt,
356         I            myThid )
357            ENDIF
358    #endif  /* ALLOW_OBCS */
359    
360  C--     Determines forcing terms based on external fields  C--     Determines forcing terms based on external fields
361  C       relaxation terms, etc.  C       relaxation terms, etc.
362          CALL EXTERNAL_FORCING_SURF(          CALL EXTERNAL_FORCING_SURF(
363       I             bi, bj, iMin, iMax, jMin, jMax,       I             bi, bj, iMin, iMax, jMin, jMax,
364       I             myThid )       I             myThid )
365    #ifdef ALLOW_AUTODIFF_TAMC
366    cph needed for KPP
367    CADJ STORE surfacetendencyU(:,:,bi,bj)
368    CADJ &     = comlev1_bibj, key=ikey, byte=isbyte
369    CADJ STORE surfacetendencyV(:,:,bi,bj)
370    CADJ &     = comlev1_bibj, key=ikey, byte=isbyte
371    CADJ STORE surfacetendencyS(:,:,bi,bj)
372    CADJ &     = comlev1_bibj, key=ikey, byte=isbyte
373    CADJ STORE surfacetendencyT(:,:,bi,bj)
374    CADJ &     = comlev1_bibj, key=ikey, byte=isbyte
375    #endif /* ALLOW_AUTODIFF_TAMC */
376    
377  #ifdef  ALLOW_GMREDI  #ifdef  ALLOW_GMREDI
378    
379    #ifdef ALLOW_AUTODIFF_TAMC
380    CADJ STORE sigmaX(:,:,:) = comlev1, key=ikey, byte=isbyte
381    CADJ STORE sigmaY(:,:,:) = comlev1, key=ikey, byte=isbyte
382    CADJ STORE sigmaR(:,:,:) = comlev1, key=ikey, byte=isbyte
383    #endif /* ALLOW_AUTODIFF_TAMC */
384  C--     Calculate iso-neutral slopes for the GM/Redi parameterisation  C--     Calculate iso-neutral slopes for the GM/Redi parameterisation
385          IF (useGMRedi) THEN          IF (useGMRedi) THEN
386            DO k=1,Nr            DO k=1,Nr
# Line 348  C--     Calculate iso-neutral slopes for Line 389  C--     Calculate iso-neutral slopes for
389       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
390       I             myThid )       I             myThid )
391            ENDDO            ENDDO
392    #ifdef ALLOW_AUTODIFF_TAMC
393            ELSE
394              DO k=1, Nr
395                CALL GMREDI_CALC_TENSOR_DUMMY(
396         I             bi, bj, iMin, iMax, jMin, jMax, k,
397         I             sigmaX, sigmaY, sigmaR,
398         I             myThid )
399              ENDDO
400    #endif /* ALLOW_AUTODIFF_TAMC */
401          ENDIF          ENDIF
402    
403    #ifdef ALLOW_AUTODIFF_TAMC
404    CADJ STORE Kwx(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte
405    CADJ STORE Kwy(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte
406    CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte
407    #endif /* ALLOW_AUTODIFF_TAMC */
408    
409  #endif  /* ALLOW_GMREDI */  #endif  /* ALLOW_GMREDI */
410    
411  #ifdef  ALLOW_KPP  #ifdef  ALLOW_KPP
# Line 356  C--     Compute KPP mixing coefficients Line 413  C--     Compute KPP mixing coefficients
413          IF (useKPP) THEN          IF (useKPP) THEN
414            CALL KPP_CALC(            CALL KPP_CALC(
415       I                  bi, bj, myTime, myThid )       I                  bi, bj, myTime, myThid )
416    #ifdef ALLOW_AUTODIFF_TAMC
417            ELSE
418              CALL KPP_CALC_DUMMY(
419         I                  bi, bj, myTime, myThid )
420    #endif /* ALLOW_AUTODIFF_TAMC */
421          ENDIF          ENDIF
422    
423    #ifdef ALLOW_AUTODIFF_TAMC
424    CADJ STORE KPPghat   (:,:,:,bi,bj)
425    CADJ &   , KPPviscAz (:,:,:,bi,bj)
426    CADJ &   , KPPdiffKzT(:,:,:,bi,bj)
427    CADJ &   , KPPdiffKzS(:,:,:,bi,bj)
428    CADJ &   , KPPfrac   (:,:  ,bi,bj)
429    CADJ &                 = comlev1_bibj, key=ikey, byte=isbyte
430    #endif /* ALLOW_AUTODIFF_TAMC */
431    
432  #endif  /* ALLOW_KPP */  #endif  /* ALLOW_KPP */
433    
434  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 368  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_ Line 440  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_
440  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
441  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
442    
443    #ifdef ALLOW_AIM
444    C       AIM - atmospheric intermediate model, physics package code.
445    C note(jmc) : phiHyd=0 at this point but is not really used in Molteni Physics
446            IF ( useAIM ) THEN
447             CALL TIMER_START('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)
448             CALL AIM_DO_ATMOS_PHYSICS( phiHyd, myTime, myThid )
449             CALL TIMER_STOP ('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)
450            ENDIF
451    #endif /* ALLOW_AIM */
452    
453    
454  C--     Start of thermodynamics loop  C--     Start of thermodynamics loop
455          DO k=Nr,1,-1          DO k=Nr,1,-1
456    #ifdef ALLOW_AUTODIFF_TAMC
457    C? Patrick Is this formula correct?
458    cph Yes, but I rewrote it.
459    cph Also, the KappaR? need the index and subscript k!
460             kkey = (ikey-1)*Nr + k
461    #endif /* ALLOW_AUTODIFF_TAMC */
462    
463  C--       km1    Points to level above k (=k-1)  C--       km1    Points to level above k (=k-1)
464  C--       kup    Cycles through 1,2 to point to layer above  C--       kup    Cycles through 1,2 to point to layer above
# Line 386  C--       kDown  Cycles through 2,1 to p Line 473  C--       kDown  Cycles through 2,1 to p
473            jMin = 1-OLy+2            jMin = 1-OLy+2
474            jMax = sNy+OLy-1            jMax = sNy+OLy-1
475    
 #ifdef ALLOW_AUTODIFF_TAMC  
 CPatrick Is this formula correct?  
          kkey = (ikey-1)*(Nr-1+1) + (k-1) + 1  
 CADJ STORE rvel  (:,:,kDown) = comlev1_bibj_k, key = kkey, byte = isbyte  
 CADJ STORE rTrans(:,:)       = comlev1_bibj_k, key = kkey, byte = isbyte  
 CADJ STORE KappaRT(:,:,:)    = comlev1_bibj_k, key = kkey, byte = isbyte  
 CADJ STORE KappaRS(:,:,:)    = comlev1_bibj_k, key = kkey, byte = isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
476  C--      Get temporary terms used by tendency routines  C--      Get temporary terms used by tendency routines
477           CALL CALC_COMMON_FACTORS (           CALL CALC_COMMON_FACTORS (
478       I        bi,bj,iMin,iMax,jMin,jMax,k,km1,kup,kDown,       I        bi,bj,iMin,iMax,jMin,jMax,k,
479       O        xA,yA,uTrans,vTrans,rTrans,rVel,maskC,maskUp,       O        xA,yA,uTrans,vTrans,rTrans,maskUp,
480       I        myThid)       I        myThid)
481    
482    #ifdef ALLOW_AUTODIFF_TAMC
483    CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte
484    CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte
485    #endif /* ALLOW_AUTODIFF_TAMC */
486    
487  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL
488  C--      Calculate the total vertical diffusivity  C--      Calculate the total vertical diffusivity
489           CALL CALC_DIFFUSIVITY(           CALL CALC_DIFFUSIVITY(
490       I        bi,bj,iMin,iMax,jMin,jMax,k,       I        bi,bj,iMin,iMax,jMin,jMax,k,
491       I        maskC,maskup,       I        maskUp,
492       O        KappaRT,KappaRS,KappaRU,KappaRV,       O        KappaRT,KappaRS,KappaRU,KappaRV,
493       I        myThid)       I        myThid)
494  #endif  #endif
# Line 415  C        and step forward storing result Line 498  C        and step forward storing result
498           IF ( tempStepping ) THEN           IF ( tempStepping ) THEN
499             CALL CALC_GT(             CALL CALC_GT(
500       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
501       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,
502       I         KappaRT,       I         KappaRT,
503       U         fVerT,       U         fVerT,
504       I         myTime, myThid)       I         myTime, myThid)
505               tauAB = 0.5d0 + abEps
506             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
507       I         bi,bj,iMin,iMax,jMin,jMax,k,       I         bi,bj,iMin,iMax,jMin,jMax,k,tauAB,
508       I         theta, gT,       I         theta, gT,
509       U         gTnm1,       U         gTnm1,
510       I         myIter, myThid)       I         myIter, myThid)
# Line 428  C        and step forward storing result Line 512  C        and step forward storing result
512           IF ( saltStepping ) THEN           IF ( saltStepping ) THEN
513             CALL CALC_GS(             CALL CALC_GS(
514       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
515       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,
516       I         KappaRS,       I         KappaRS,
517       U         fVerS,       U         fVerS,
518       I         myTime, myThid)       I         myTime, myThid)
519               tauAB = 0.5d0 + abEps
520             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
521       I         bi,bj,iMin,iMax,jMin,jMax,k,       I         bi,bj,iMin,iMax,jMin,jMax,k,tauAB,
522       I         salt, gS,       I         salt, gS,
523       U         gSnm1,       U         gSnm1,
524       I         myIter, myThid)       I         myIter, myThid)
# Line 441  C        and step forward storing result Line 526  C        and step forward storing result
526    
527  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
528  C--      Apply open boundary conditions  C--      Apply open boundary conditions
529           IF (openBoundaries) THEN           IF (useOBCS) THEN
530             CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )             CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )
531           END IF           END IF
532  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
# Line 460  C--     end of thermodynamic k loop (Nr: Line 545  C--     end of thermodynamic k loop (Nr:
545    
546    
547  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
548  CPatrick? What about this one?  C? Patrick? What about this one?
549             maximpl = 6  cph Keys iikey and idkey don't seem to be needed
550             iikey = (ikey-1)*maximpl  cph since storing occurs on different tape for each
551    cph impldiff call anyways.
552    cph Thus, common block comlev1_impl isn't needed either.
553    cph Storing below needed in the case useGMREDI.
554            iikey = (ikey-1)*maximpl
555  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
556    
557  C--     Implicit diffusion  C--     Implicit diffusion
558          IF (implicitDiffusion) THEN          IF (implicitDiffusion) THEN
559    
560            IF (tempStepping) THEN           IF (tempStepping) THEN
561  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
562              idkey = iikey + 1              idkey = iikey + 1
563    CADJ STORE gTNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
564  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
565              CALL IMPLDIFF(              CALL IMPLDIFF(
566       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
# Line 482  C--     Implicit diffusion Line 572  C--     Implicit diffusion
572           IF (saltStepping) THEN           IF (saltStepping) THEN
573  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
574           idkey = iikey + 2           idkey = iikey + 2
575    CADJ STORE gSNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
576  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
577              CALL IMPLDIFF(              CALL IMPLDIFF(
578       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
# Line 492  C--     Implicit diffusion Line 583  C--     Implicit diffusion
583    
584  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
585  C--      Apply open boundary conditions  C--      Apply open boundary conditions
586           IF (openBoundaries) THEN           IF (useOBCS) THEN
587             DO K=1,Nr             DO K=1,Nr
588               CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )               CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )
589             ENDDO             ENDDO
# Line 502  C--      Apply open boundary conditions Line 593  C--      Apply open boundary conditions
593  C--     End If implicitDiffusion  C--     End If implicitDiffusion
594          ENDIF          ENDIF
595    
596    C--     Start computation of dynamics
597            iMin = 1-OLx+2
598            iMax = sNx+OLx-1
599            jMin = 1-OLy+2
600            jMax = sNy+OLy-1
601    
602    C--     Explicit part of the Surface Potentiel Gradient (add in TIMESTEP)
603    C       (note: this loop will be replaced by CALL CALC_GRAD_ETA)
604            IF (implicSurfPress.NE.1.) THEN
605              CALL CALC_GRAD_PHI_SURF(
606         I         bi,bj,iMin,iMax,jMin,jMax,
607         I         etaN,
608         O         phiSurfX,phiSurfY,
609         I         myThid )                        
610            ENDIF
611    
612  C--     Start of dynamics loop  C--     Start of dynamics loop
613          DO k=1,Nr          DO k=1,Nr
# Line 515  C--       kDown  Cycles through 2,1 to p Line 620  C--       kDown  Cycles through 2,1 to p
620            kup  = 1+MOD(k+1,2)            kup  = 1+MOD(k+1,2)
621            kDown= 1+MOD(k,2)            kDown= 1+MOD(k,2)
622    
           iMin = 1-OLx+2  
           iMax = sNx+OLx-1  
           jMin = 1-OLy+2  
           jMax = sNy+OLy-1  
   
 C--      Calculate buoyancy  
          CALL FIND_RHO(  
      I        bi, bj, iMin, iMax, jMin, jMax, km1, km1, eosType,  
      O        rhoKm1,  
      I        myThid )  
          CALL CALC_BUOYANCY(  
      I        bi,bj,iMin,iMax,jMin,jMax,k,rhoKm1,  
      O        buoyKm1,  
      I        myThid )  
          CALL FIND_RHO(  
      I        bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,  
      O        rhoK,  
      I        myThid )  
          CALL CALC_BUOYANCY(  
      I        bi,bj,iMin,iMax,jMin,jMax,k,rhoK,  
      O        buoyK,  
      I        myThid )  
   
623  C--      Integrate hydrostatic balance for phiHyd with BC of  C--      Integrate hydrostatic balance for phiHyd with BC of
624  C--      phiHyd(z=0)=0  C        phiHyd(z=0)=0
625           CALL CALC_PHI_HYD(  C        distinguishe between Stagger and Non Stagger time stepping
626       I        bi,bj,iMin,iMax,jMin,jMax,k,buoyKm1,buoyK,           IF (staggerTimeStep) THEN
627               CALL CALC_PHI_HYD(
628         I        bi,bj,iMin,iMax,jMin,jMax,k,
629         I        gTnm1, gSnm1,
630         U        phiHyd,
631         I        myThid )
632             ELSE
633               CALL CALC_PHI_HYD(
634         I        bi,bj,iMin,iMax,jMin,jMax,k,
635         I        theta, salt,
636       U        phiHyd,       U        phiHyd,
637       I        myThid )       I        myThid )
638             ENDIF
639    
640  C--      Calculate accelerations in the momentum equations (gU, gV, ...)  C--      Calculate accelerations in the momentum equations (gU, gV, ...)
641  C        and step forward storing the result in gUnm1, gVnm1, etc...  C        and step forward storing the result in gUnm1, gVnm1, etc...
# Line 555  C        and step forward storing the re Line 647  C        and step forward storing the re
647       I         myTime, myThid)       I         myTime, myThid)
648             CALL TIMESTEP(             CALL TIMESTEP(
649       I         bi,bj,iMin,iMax,jMin,jMax,k,       I         bi,bj,iMin,iMax,jMin,jMax,k,
650         I         phiHyd, phiSurfX, phiSurfY,
651       I         myIter, myThid)       I         myIter, myThid)
652    
653  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
654  C--      Apply open boundary conditions  C--      Apply open boundary conditions
655           IF (openBoundaries) THEN           IF (useOBCS) THEN
656             CALL OBCS_APPLY_UV( bi, bj, k, gUnm1, gVnm1, myThid )             CALL OBCS_APPLY_UV( bi, bj, k, gUnm1, gVnm1, myThid )
657           END IF           END IF
658  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
# Line 587  C--     Implicit viscosity Line 680  C--     Implicit viscosity
680          IF (implicitViscosity.AND.momStepping) THEN          IF (implicitViscosity.AND.momStepping) THEN
681  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
682            idkey = iikey + 3            idkey = iikey + 3
683    CADJ STORE gUNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
684  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
685            CALL IMPLDIFF(            CALL IMPLDIFF(
686       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
# Line 595  C--     Implicit viscosity Line 689  C--     Implicit viscosity
689       I         myThid )       I         myThid )
690  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
691            idkey = iikey + 4            idkey = iikey + 4
692    CADJ STORE gVNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
693  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
694            CALL IMPLDIFF(            CALL IMPLDIFF(
695       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
# Line 604  C--     Implicit viscosity Line 699  C--     Implicit viscosity
699    
700  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
701  C--      Apply open boundary conditions  C--      Apply open boundary conditions
702           IF (openBoundaries) THEN           IF (useOBCS) THEN
703             DO K=1,Nr             DO K=1,Nr
704               CALL OBCS_APPLY_UV( bi, bj, k, gUnm1, gVnm1, myThid )               CALL OBCS_APPLY_UV( bi, bj, k, gUnm1, gVnm1, myThid )
705             ENDDO             ENDDO
# Line 614  C--      Apply open boundary conditions Line 709  C--      Apply open boundary conditions
709  #ifdef    INCLUDE_CD_CODE  #ifdef    INCLUDE_CD_CODE
710  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
711            idkey = iikey + 5            idkey = iikey + 5
712    CADJ STORE vVelD(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
713  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
714            CALL IMPLDIFF(            CALL IMPLDIFF(
715       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
# Line 622  C--      Apply open boundary conditions Line 718  C--      Apply open boundary conditions
718       I         myThid )       I         myThid )
719  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
720            idkey = iikey + 6            idkey = iikey + 6
721    CADJ STORE uVelD(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
722  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
723            CALL IMPLDIFF(            CALL IMPLDIFF(
724       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
# Line 632  C--      Apply open boundary conditions Line 729  C--      Apply open boundary conditions
729  C--     End If implicitViscosity.AND.momStepping  C--     End If implicitViscosity.AND.momStepping
730          ENDIF          ENDIF
731    
732    Cjmc : add for phiHyd output <- but not working if multi tile per CPU
733    c       IF ( DIFFERENT_MULTIPLE(dumpFreq,myTime+deltaTClock,myTime)
734    c    &  .AND. buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN
735    c         WRITE(suff,'(I10.10)') myIter+1
736    c         CALL WRITE_FLD_XYZ_RL('PH.',suff,phiHyd,myIter+1,myThid)
737    c       ENDIF
738    Cjmc(end)
739    
740    #ifdef ALLOW_TIMEAVE
741            IF (taveFreq.GT.0.) THEN
742              CALL TIMEAVE_CUMUL_1T(phiHydtave, phiHyd, Nr,
743         I                              deltaTclock, bi, bj, myThid)
744              IF (ivdc_kappa.NE.0.) THEN
745                CALL TIMEAVE_CUMULATE(ConvectCountTave, ConvectCount, Nr,
746         I                              deltaTclock, bi, bj, myThid)
747              ENDIF
748            ENDIF
749    #endif /* ALLOW_TIMEAVE */
750    
751         ENDDO         ENDDO
752        ENDDO        ENDDO
753    
754        RETURN        RETURN
755        END        END
   
   
 C--      Cumulative diagnostic calculations (ie. time-averaging)  
 #ifdef INCLUDE_DIAGNOSTICS_INTERFACE_CODE  
 c        IF (taveFreq.GT.0.) THEN  
 c         CALL DO_TIME_AVERAGES(  
 c    I                           myTime, myIter, bi, bj, k, kup, kDown,  
 c    I                           ConvectCount,  
 c    I                           myThid )  
 c        ENDIF  
 #endif  
   
Legend:
Removed from v.1.54.2.8  
changed lines
  Added in v.1.68
  ViewVC Help
Powered by ViewVC 1.1.22