/[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.58 by adcroft, Fri Feb 2 21:04:48 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
# Line 65  C                                      v Line 69  C                                      v
69  C     rhoK, rhoKM1   - Density at current level, and level above  C     rhoK, rhoKM1   - Density 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 79  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 94  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 phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
100          _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 101  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 113  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 137  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 194  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          maskC  (i,j) = 0. _d 0          phiSurfX(i,j) = 0. _d 0
196            phiSurfY(i,j) = 0. _d 0
197         ENDDO         ENDDO
198        ENDDO        ENDDO
199    
# Line 208  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 237  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 253  C--     Set up work arrays that need val Line 250  C--     Set up work arrays that need val
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 268  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 285  C--       Integrate continuity verticall Line 292  C--       Integrate continuity verticall
292    
293  #ifdef    ALLOW_OBCS  #ifdef    ALLOW_OBCS
294  #ifdef    ALLOW_NONHYDROSTATIC  #ifdef    ALLOW_NONHYDROSTATIC
295  C--       Calculate future values on open boundaries  C--       Apply OBC to W if in N-H mode
296            IF (useOBCS.AND.nonHydrostatic) THEN            IF (useOBCS.AND.nonHydrostatic) THEN
297              CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )              CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )
298            ENDIF            ENDIF
# Line 296  C--       Calculate gradients of potenti Line 303  C--       Calculate gradients of potenti
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,       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,       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 321  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  #ifdef  ALLOW_OBCS
352  C--     Calculate future values on open boundaries  C--     Calculate future values on open boundaries
353          IF (useOBCS) THEN          IF (useOBCS) THEN
# Line 340  C       relaxation terms, etc. Line 362  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 360  C--     Calculate iso-neutral slopes for Line 399  C--     Calculate iso-neutral slopes for
399            ENDDO            ENDDO
400  #endif /* ALLOW_AUTODIFF_TAMC */  #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 367  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 392  C note(jmc) : phiHyd=0 at this point but Line 453  C note(jmc) : phiHyd=0 at this point but
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 406  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 435  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 448  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 480  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 502  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 522  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 535  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  
   
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  C        distinguishe between Stagger and Non Stagger time stepping  C        distinguishe between Stagger and Non Stagger time stepping
# Line 566  C        and step forward storing the re Line 646  C        and step forward storing the re
646       U         fVerU, fVerV,       U         fVerU, fVerV,
647       I         myTime, myThid)       I         myTime, myThid)
648             CALL TIMESTEP(             CALL TIMESTEP(
649       I         bi,bj,iMin,iMax,jMin,jMax,k,phiHyd,       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
# Line 599  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 607  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 626  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 634  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 644  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    

Legend:
Removed from v.1.58  
changed lines
  Added in v.1.68

  ViewVC Help
Powered by ViewVC 1.1.22