/[MITgcm]/MITgcm/model/src/thermodynamics.F
ViewVC logotype

Diff of /MITgcm/model/src/thermodynamics.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.11 by heimbach, Thu Sep 27 20:12:10 2001 UTC revision 1.39 by jmc, Thu May 1 22:30:33 2003 UTC
# Line 2  C $Header$ Line 2  C $Header$
2  C $Name$  C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
5    #ifdef ALLOW_AUTODIFF_TAMC
6    # ifdef ALLOW_GMREDI
7    #  include "GMREDI_OPTIONS.h"
8    # endif
9    # ifdef ALLOW_KPP
10    #  include "KPP_OPTIONS.h"
11    # endif
12    cswdice --- add ----
13    #ifdef ALLOW_THERM_SEAICE
14    #include "ICE.h"
15    #endif
16    cswdice ------
17    #endif /* ALLOW_AUTODIFF_TAMC */
18    
19  CBOP  CBOP
20  C     !ROUTINE: THERMODYNAMICS  C     !ROUTINE: THERMODYNAMICS
# Line 75  C     == Global variables === Line 88  C     == Global variables ===
88  # include "tamc.h"  # include "tamc.h"
89  # include "tamc_keys.h"  # include "tamc_keys.h"
90  # include "FFIELDS.h"  # include "FFIELDS.h"
91    # include "EOS.h"
92  # ifdef ALLOW_KPP  # ifdef ALLOW_KPP
93  #  include "KPP.h"  #  include "KPP.h"
94  # endif  # endif
# Line 109  C                                      i Line 123  C                                      i
123  C                                      so we need an fVer for each  C                                      so we need an fVer for each
124  C                                      variable.  C                                      variable.
125  C     rhoK, rhoKM1   - Density at current level, and level above  C     rhoK, rhoKM1   - Density at current level, and level above
 C     phiHyd         - Hydrostatic part of the potential phiHydi.  
 C                      In z coords phiHydiHyd is the hydrostatic  
 C                      Potential (=pressure/rho0) anomaly  
 C                      In p coords phiHydiHyd is the geopotential  
 C                      surface height anomaly.  
126  C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)  C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)
127  C     phiSurfY             or geopotentiel (atmos) in X and Y direction  C     phiSurfY             or geopotentiel (atmos) in X and Y direction
128  C     KappaRT,       - Total diffusion in vertical for T and S.  C     KappaRT,       - Total diffusion in vertical for T and S.
129  C     KappaRS          (background + spatially varying, isopycnal term).  C     KappaRS          (background + spatially varying, isopycnal term).
130    C     useVariableK   = T when vertical diffusion is not constant
131  C     iMin, iMax     - Ranges and sub-block indices on which calculations  C     iMin, iMax     - Ranges and sub-block indices on which calculations
132  C     jMin, jMax       are applied.  C     jMin, jMax       are applied.
133  C     bi, bj  C     bi, bj
# Line 133  C                      index into fVerTe Line 143  C                      index into fVerTe
143        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
144        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
145        _RL fVerTr1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerTr1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
       _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
146        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
147        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
148        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 145  C                      index into fVerTe Line 154  C                      index into fVerTe
154        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
155  C     This is currently used by IVDC and Diagnostics  C     This is currently used by IVDC and Diagnostics
156        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
157          LOGICAL useVariableK
158        INTEGER iMin, iMax        INTEGER iMin, iMax
159        INTEGER jMin, jMax        INTEGER jMin, jMax
160        INTEGER bi, bj        INTEGER bi, bj
161        INTEGER i, j        INTEGER i, j
162        INTEGER k, km1, kup, kDown        INTEGER k, km1, kup, kDown
163    
 Cjmc : add for phiHyd output <- but not working if multi tile per CPU  
 c     CHARACTER*(MAX_LEN_MBUF) suff  
 c     LOGICAL  DIFFERENT_MULTIPLE  
 c     EXTERNAL DIFFERENT_MULTIPLE  
 Cjmc(end)  
164  CEOP  CEOP
165    
166  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
167  C--   dummy statement to end declaration part  C--   dummy statement to end declaration part
168        ikey = 1        ikey = 1
169          itdkey = 1
170  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
171    
172  C--   Set up work arrays with valid (i.e. not NaN) values  C--   Set up work arrays with valid (i.e. not NaN) values
# Line 174  C     uninitialised but inert locations. Line 180  C     uninitialised but inert locations.
180          yA(i,j)      = 0. _d 0          yA(i,j)      = 0. _d 0
181          uTrans(i,j)  = 0. _d 0          uTrans(i,j)  = 0. _d 0
182          vTrans(i,j)  = 0. _d 0          vTrans(i,j)  = 0. _d 0
         DO k=1,Nr  
          phiHyd(i,j,k)  = 0. _d 0  
          sigmaX(i,j,k) = 0. _d 0  
          sigmaY(i,j,k) = 0. _d 0  
          sigmaR(i,j,k) = 0. _d 0  
         ENDDO  
         rhoKM1 (i,j) = 0. _d 0  
183          rhok   (i,j) = 0. _d 0          rhok   (i,j) = 0. _d 0
184          phiSurfX(i,j) = 0. _d 0          phiSurfX(i,j) = 0. _d 0
185          phiSurfY(i,j) = 0. _d 0          phiSurfY(i,j) = 0. _d 0
# Line 198  CHPF$ INDEPENDENT Line 197  CHPF$ INDEPENDENT
197  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
198  C--    HPF directive to help TAMC  C--    HPF directive to help TAMC
199  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS
200  CHPF$&                  ,phiHyd,utrans,vtrans,xA,yA  CHPF$&                  ,utrans,vtrans,xA,yA
201  CHPF$&                  ,KappaRT,KappaRS  CHPF$&                  ,KappaRT,KappaRS
202  CHPF$&                  )  CHPF$&                  )
203  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
# Line 213  CHPF$&                  ) Line 212  CHPF$&                  )
212            act3 = myThid - 1            act3 = myThid - 1
213            max3 = nTx*nTy            max3 = nTx*nTy
214            act4 = ikey_dynamics - 1            act4 = ikey_dynamics - 1
215            ikey = (act1 + 1) + act2*max1            itdkey = (act1 + 1) + act2*max1
216       &                      + act3*max1*max2       &                      + act3*max1*max2
217       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
218  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
# Line 228  C--     Set up work arrays that need val Line 227  C--     Set up work arrays that need val
227            fVerS  (i,j,2) = 0. _d 0            fVerS  (i,j,2) = 0. _d 0
228            fVerTr1(i,j,1) = 0. _d 0            fVerTr1(i,j,1) = 0. _d 0
229            fVerTr1(i,j,2) = 0. _d 0            fVerTr1(i,j,2) = 0. _d 0
230              rhoKM1 (i,j)   = 0. _d 0
231           ENDDO           ENDDO
232          ENDDO          ENDDO
233    
# Line 235  C--     Set up work arrays that need val Line 235  C--     Set up work arrays that need val
235           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
236            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
237  C This is currently also used by IVDC and Diagnostics  C This is currently also used by IVDC and Diagnostics
238               sigmaX(i,j,k) = 0. _d 0
239               sigmaY(i,j,k) = 0. _d 0
240               sigmaR(i,j,k) = 0. _d 0
241             ConvectCount(i,j,k) = 0.             ConvectCount(i,j,k) = 0.
242             KappaRT(i,j,k) = 0. _d 0             KappaRT(i,j,k)    = 0. _d 0
243             KappaRS(i,j,k) = 0. _d 0             KappaRS(i,j,k)    = 0. _d 0
244  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
245             gT(i,j,k,bi,bj) = 0. _d 0  cph all the following init. are necessary for TAF
246             gS(i,j,k,bi,bj) = 0. _d 0  cph although some of these are re-initialised later.
247  #ifdef ALLOW_PASSIVE_TRACER             gT(i,j,k,bi,bj)   = 0. _d 0
248               gS(i,j,k,bi,bj)   = 0. _d 0
249    # ifdef ALLOW_PASSIVE_TRACER
250             gTr1(i,j,k,bi,bj) = 0. _d 0             gTr1(i,j,k,bi,bj) = 0. _d 0
251  #endif  # endif
252  #endif  # ifdef ALLOW_GMREDI
253               Kwx(i,j,k,bi,bj)  = 0. _d 0
254               Kwy(i,j,k,bi,bj)  = 0. _d 0
255               Kwz(i,j,k,bi,bj)  = 0. _d 0
256    #  ifdef GM_NON_UNITY_DIAGONAL
257               Kux(i,j,k,bi,bj)  = 0. _d 0
258               Kvy(i,j,k,bi,bj)  = 0. _d 0
259    #  endif
260    #  ifdef GM_EXTRA_DIAGONAL
261               Kuz(i,j,k,bi,bj)  = 0. _d 0
262               Kvz(i,j,k,bi,bj)  = 0. _d 0
263    #  endif
264    #  ifdef GM_BOLUS_ADVEC
265               GM_PsiX(i,j,k,bi,bj)  = 0. _d 0
266               GM_PsiY(i,j,k,bi,bj)  = 0. _d 0
267    #  endif
268    #  ifdef GM_VISBECK_VARIABLE_K
269               VisbeckK(i,j,bi,bj)   = 0. _d 0
270    #  endif
271    # endif /* ALLOW_GMREDI */
272    #endif /* ALLOW_AUTODIFF_TAMC */
273            ENDDO            ENDDO
274           ENDDO           ENDDO
275          ENDDO          ENDDO
276    
277          iMin = 1-OLx+1          iMin = 1-OLx
278          iMax = sNx+OLx          iMax = sNx+OLx
279          jMin = 1-OLy+1          jMin = 1-OLy
280          jMax = sNy+OLy          jMax = sNy+OLy
281    
282    
283  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
284  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
285  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
286    CADJ STORE totphihyd
287    CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
288  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
289  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
290  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
291  #endif  #endif
292  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
293    
# Line 272  C? Patrick, is this formula correct now Line 299  C? Patrick, is this formula correct now
299  C? Do we still need this?  C? Do we still need this?
300  cph kkey formula corrected.  cph kkey formula corrected.
301  cph Needed for rhok, rhokm1, in the case useGMREDI.  cph Needed for rhok, rhokm1, in the case useGMREDI.
302           kkey = (ikey-1)*Nr + k           kkey = (itdkey-1)*Nr + k
 CADJ STORE rhokm1(:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte  
 CADJ STORE rhok  (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte  
303  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
304    
305  C--       Integrate continuity vertically for vertical velocity  C--       Integrate continuity vertically for vertical velocity
306            CALL INTEGRATE_FOR_W(  c         CALL INTEGRATE_FOR_W(
307       I                         bi, bj, k, uVel, vVel,  c    I                         bi, bj, k, uVel, vVel,
308       O                         wVel,  c    O                         wVel,
309       I                         myThid )  c    I                         myThid )
310    
311  #ifdef    ALLOW_OBCS  #ifdef    ALLOW_OBCS
312  #ifdef    ALLOW_NONHYDROSTATIC  #ifdef    ALLOW_NONHYDROSTATIC
313  C--       Apply OBC to W if in N-H mode  C--       Apply OBC to W if in N-H mode
314            IF (useOBCS.AND.nonHydrostatic) THEN  c         IF (useOBCS.AND.nonHydrostatic) THEN
315              CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )  c           CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )
316            ENDIF  c         ENDIF
317  #endif    /* ALLOW_NONHYDROSTATIC */  #endif    /* ALLOW_NONHYDROSTATIC */
318  #endif    /* ALLOW_OBCS */  #endif    /* ALLOW_OBCS */
319    
320    C--       Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
321    C--       MOST of THERMODYNAMICS will be disabled
322    #ifndef SINGLE_LAYER_MODE
323    
324  C--       Calculate gradients of potential density for isoneutral  C--       Calculate gradients of potential density for isoneutral
325  C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)  C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
326  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
# Line 301  CADJ STORE theta(:,:,k,bi,bj) = comlev1_ Line 330  CADJ STORE theta(:,:,k,bi,bj) = comlev1_
330  CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
331  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
332              CALL FIND_RHO(              CALL FIND_RHO(
333       I        bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,       I        bi, bj, iMin, iMax, jMin, jMax, k, k,
334       I        theta, salt,       I        theta, salt,
335       O        rhoK,       O        rhoK,
336       I        myThid )       I        myThid )
337    
338              IF (k.GT.1) THEN              IF (k.GT.1) THEN
339  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
340  CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
341  CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
342  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
343               CALL FIND_RHO(               CALL FIND_RHO(
344       I        bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,       I        bi, bj, iMin, iMax, jMin, jMax, k-1, k,
345       I        theta, salt,       I        theta, salt,
346       O        rhoKm1,       O        rhoKm1,
347       I        myThid )       I        myThid )
# Line 323  CADJ STORE salt (:,:,k-1,bi,bj) = comlev Line 353  CADJ STORE salt (:,:,k-1,bi,bj) = comlev
353       I             myThid )       I             myThid )
354            ENDIF            ENDIF
355    
356    #ifdef ALLOW_AUTODIFF_TAMC
357    CADJ STORE rhok   (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte
358    CADJ STORE rhokm1 (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte
359    #endif /* ALLOW_AUTODIFF_TAMC */
360  C--       Implicit Vertical Diffusion for Convection  C--       Implicit Vertical Diffusion for Convection
361  c ==> should use sigmaR !!!  c ==> should use sigmaR !!!
362            IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN            IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN
# Line 333  c ==> should use sigmaR !!! Line 367  c ==> should use sigmaR !!!
367       I        myTime, myIter, myThid)       I        myTime, myIter, myThid)
368            ENDIF            ENDIF
369    
370    #endif /* SINGLE_LAYER_MODE */
371    
372  C--     end of diagnostic k loop (Nr:1)  C--     end of diagnostic k loop (Nr:1)
373          ENDDO          ENDDO
374    
375  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
376  cph avoids recomputation of integrate_for_w  cph avoids recomputation of integrate_for_w
377  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
378  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
379    
380  #ifdef  ALLOW_OBCS  #ifdef  ALLOW_OBCS
381  C--     Calculate future values on open boundaries  C--     Calculate future values on open boundaries
382          IF (useOBCS) THEN          IF (useOBCS) THEN
383            CALL OBCS_CALC( bi, bj, myTime+deltaT,            CALL OBCS_CALC( bi, bj, myTime+deltaT, myIter+1,
384       I            uVel, vVel, wVel, theta, salt,       I            uVel, vVel, wVel, theta, salt,
385       I            myThid )       I            myThid )
386          ENDIF          ENDIF
387  #endif  /* ALLOW_OBCS */  #endif  /* ALLOW_OBCS */
388    
389    
390    c********************************************
391    cswdice --- add ---
392    #ifdef ALLOW_THERM_SEAICE
393    C--     Determines forcing terms based on external fields
394    c--     including effects from ice
395            CALL ICE_FORCING(
396         I             bi, bj, iMin, iMax, jMin, jMax,
397         I             myThid )
398    #else
399    
400    cswdice --- end add ---
401    
402  C--     Determines forcing terms based on external fields  C--     Determines forcing terms based on external fields
403  C       relaxation terms, etc.  C       relaxation terms, etc.
404          CALL EXTERNAL_FORCING_SURF(          CALL EXTERNAL_FORCING_SURF(
405       I             bi, bj, iMin, iMax, jMin, jMax,       I             bi, bj, iMin, iMax, jMin, jMax,
406       I             myThid )       I             myThid )
407    cswdice --- add ----
408    #endif
409    cswdice --- end add ---
410    c******************************************
411    
412    
413    
414    
415  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
416  cph needed for KPP  cph needed for KPP
417  CADJ STORE surfacetendencyU(:,:,bi,bj)  CADJ STORE surfacetendencyU(:,:,bi,bj)
418  CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
419  CADJ STORE surfacetendencyV(:,:,bi,bj)  CADJ STORE surfacetendencyV(:,:,bi,bj)
420  CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
421  CADJ STORE surfacetendencyS(:,:,bi,bj)  CADJ STORE surfacetendencyS(:,:,bi,bj)
422  CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
423  CADJ STORE surfacetendencyT(:,:,bi,bj)  CADJ STORE surfacetendencyT(:,:,bi,bj)
424  CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
425  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
426    
427    C--     Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
428    C--     MOST of THERMODYNAMICS will be disabled
429    #ifndef SINGLE_LAYER_MODE
430    
431  #ifdef  ALLOW_GMREDI  #ifdef  ALLOW_GMREDI
432    
433  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
434  CADJ STORE sigmaX(:,:,:) = comlev1, key=ikey, byte=isbyte  cph storing here is needed only for one GMREDI_OPTIONS:
435  CADJ STORE sigmaY(:,:,:) = comlev1, key=ikey, byte=isbyte  cph define GM_BOLUS_ADVEC
436  CADJ STORE sigmaR(:,:,:) = comlev1, key=ikey, byte=isbyte  cph but I've avoided the #ifdef for now, in case more things change
437    CADJ STORE sigmaX(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte
438    CADJ STORE sigmaY(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte
439    CADJ STORE sigmaR(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte
440  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
441    
442  C--     Calculate iso-neutral slopes for the GM/Redi parameterisation  C--     Calculate iso-neutral slopes for the GM/Redi parameterisation
443          IF (useGMRedi) THEN          IF (useGMRedi) THEN
444            DO k=1,Nr            CALL GMREDI_CALC_TENSOR(
445              CALL GMREDI_CALC_TENSOR(       I             bi, bj, iMin, iMax, jMin, jMax,
      I             bi, bj, iMin, iMax, jMin, jMax, k,  
446       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
447       I             myThid )       I             myThid )
           ENDDO  
448  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
449          ELSE          ELSE
450            DO k=1, Nr            CALL GMREDI_CALC_TENSOR_DUMMY(
451              CALL GMREDI_CALC_TENSOR_DUMMY(       I             bi, bj, iMin, iMax, jMin, jMax,
      I             bi, bj, iMin, iMax, jMin, jMax, k,  
452       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
453       I             myThid )       I             myThid )
           ENDDO  
454  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
455          ENDIF          ENDIF
456    
457  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
458  CADJ STORE Kwx(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE Kwx(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte
459  CADJ STORE Kwy(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE Kwy(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte
460  CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte
461  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
462    
463  #endif  /* ALLOW_GMREDI */  #endif  /* ALLOW_GMREDI */
# Line 418  CADJ STORE KPPghat   (:,:,:,bi,bj) Line 479  CADJ STORE KPPghat   (:,:,:,bi,bj)
479  CADJ &   , KPPdiffKzT(:,:,:,bi,bj)  CADJ &   , KPPdiffKzT(:,:,:,bi,bj)
480  CADJ &   , KPPdiffKzS(:,:,:,bi,bj)  CADJ &   , KPPdiffKzS(:,:,:,bi,bj)
481  CADJ &   , KPPfrac   (:,:  ,bi,bj)  CADJ &   , KPPfrac   (:,:  ,bi,bj)
482  CADJ &                 = comlev1_bibj, key=ikey, byte=isbyte  CADJ &                 = comlev1_bibj, key=itdkey, byte=isbyte
483  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
484    
485  #endif  /* ALLOW_KPP */  #endif  /* ALLOW_KPP */
486    
487  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
488  CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key=itdkey, byte=isbyte
489  CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key=itdkey, byte=isbyte
490  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
491  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
492  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
493  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
494  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
495  CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
496  #endif  #endif
497  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
498    
499  #ifdef ALLOW_AIM  #ifdef ALLOW_AIM
500  C       AIM - atmospheric intermediate model, physics package code.  C       AIM - atmospheric intermediate model, physics package code.
 C note(jmc) : phiHyd=0 at this point but is not really used in Molteni Physics  
501          IF ( useAIM ) THEN          IF ( useAIM ) THEN
502           CALL TIMER_START('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)           CALL TIMER_START('AIM_DO_PHYSICS   [THERMODYNAMICS]', myThid)
503           CALL AIM_DO_ATMOS_PHYSICS( phiHyd, bi, bj, myTime, myThid )           CALL AIM_DO_PHYSICS( bi, bj, myTime, myIter, myThid )
504           CALL TIMER_STOP ('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)           CALL TIMER_STOP( 'AIM_DO_PHYSICS   [THERMODYNAMICS]', myThid)
505          ENDIF          ENDIF
506  #endif /* ALLOW_AIM */  #endif /* ALLOW_AIM */
507    
508    #ifndef DISABLE_MULTIDIM_ADVECTION
509  C--     Some advection schemes are better calculated using a multi-dimensional  C--     Some advection schemes are better calculated using a multi-dimensional
510  C       method in the absence of any other terms and, if used, is done here.  C       method in the absence of any other terms and, if used, is done here.
511    C
512  #ifdef ALLOW_MULTIDIM_ADVECTION  C The CPP flag DISABLE_MULTIDIM_ADVECTION is currently unset in GAD_OPTIONS.h
513          IF (multiDimAdvection) THEN  C The default is to use multi-dimensinal advection for non-linear advection
514           IF (tempStepping .AND.  C schemes. However, for the sake of efficiency of the adjoint it is necessary
515       &       tempAdvScheme.NE.ENUM_CENTERED_2ND .AND.  C to be able to exclude this scheme to avoid excessive storage and
516       &       tempAdvScheme.NE.ENUM_UPWIND_3RD .AND.  C recomputation. It *is* differentiable, if you need it.
517       &       tempAdvScheme.NE.ENUM_CENTERED_4TH ) THEN  C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to
518    C disable this section of code.
519            IF (tempMultiDimAdvec) THEN
520            CALL GAD_ADVECTION(bi,bj,tempAdvScheme,GAD_TEMPERATURE,            CALL GAD_ADVECTION(bi,bj,tempAdvScheme,GAD_TEMPERATURE,
521       U                      theta,gT,       U                      theta,gT,
522       I                      myTime,myIter,myThid)       I                      myTime,myIter,myThid)
523           ENDIF          ENDIF
524           IF (saltStepping .AND.          IF (saltMultiDimAdvec) THEN
      &       saltAdvScheme.NE.ENUM_CENTERED_2ND .AND.  
      &       saltAdvScheme.NE.ENUM_UPWIND_3RD .AND.  
      &       saltAdvScheme.NE.ENUM_CENTERED_4TH ) THEN  
525            CALL GAD_ADVECTION(bi,bj,saltAdvScheme,GAD_SALINITY,            CALL GAD_ADVECTION(bi,bj,saltAdvScheme,GAD_SALINITY,
526       U                      salt,gS,       U                      salt,gS,
527       I                      myTime,myIter,myThid)       I                      myTime,myIter,myThid)
          ENDIF  
528          ENDIF          ENDIF
529  #endif /* ALLOW_MULTIDIM_ADVECTION */  C Since passive tracers are configurable separately from T,S we
530    C call the multi-dimensional method for PTRACERS regardless
531    C of whether multiDimAdvection is set or not.
532    #ifdef ALLOW_PTRACERS
533            IF ( usePTRACERS ) THEN
534             CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid )
535            ENDIF
536    #endif /* ALLOW_PTRACERS */
537    #endif /* DISABLE_MULTIDIM_ADVECTION */
538    
539  C--     Start of thermodynamics loop  C--     Start of thermodynamics loop
540          DO k=Nr,1,-1          DO k=Nr,1,-1
# Line 475  C--     Start of thermodynamics loop Line 542  C--     Start of thermodynamics loop
542  C? Patrick Is this formula correct?  C? Patrick Is this formula correct?
543  cph Yes, but I rewrote it.  cph Yes, but I rewrote it.
544  cph Also, the KappaR? need the index and subscript k!  cph Also, the KappaR? need the index and subscript k!
545           kkey = (ikey-1)*Nr + k           kkey = (itdkey-1)*Nr + k
546  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
547    
548  C--       km1    Points to level above k (=k-1)  C--       km1    Points to level above k (=k-1)
# Line 497  C--       Get temporary terms used by te Line 564  C--       Get temporary terms used by te
564       O         xA,yA,uTrans,vTrans,rTrans,maskUp,       O         xA,yA,uTrans,vTrans,rTrans,maskUp,
565       I         myThid)       I         myThid)
566    
567    #ifdef ALLOW_GMREDI
568    
569    C--   Residual transp = Bolus transp + Eulerian transp
570              IF (useGMRedi) THEN
571                CALL GMREDI_CALC_UVFLOW(
572         &            uTrans, vTrans, bi, bj, k, myThid)
573                IF (K.GE.2) CALL GMREDI_CALC_WFLOW(
574         &                    rTrans, bi, bj, k, myThid)
575              ENDIF
576    
577    #ifdef ALLOW_AUTODIFF_TAMC
578    #ifdef GM_BOLUS_ADVEC
579    CADJ STORE uTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
580    CADJ STORE vTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
581    CADJ STORE rTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
582    #endif
583    #endif /* ALLOW_AUTODIFF_TAMC */
584    
585    #endif /* ALLOW_GMREDI */
586    
587  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
588  CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte
589  CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte
# Line 530  C        and step forward storing result Line 617  C        and step forward storing result
617       I         theta, gT,       I         theta, gT,
618       I         myIter, myThid)       I         myIter, myThid)
619           ENDIF           ENDIF
620    cswdice ---- add ---
621    #ifdef ALLOW_THERM_SEAICE
622           if (k.eq.1) then
623            call ICE_FREEZE(bi, bj, iMin, iMax, jMin, jMax, myThid )
624           endif
625    #endif
626    cswdice -- end add ---
627           IF ( saltStepping ) THEN           IF ( saltStepping ) THEN
628             CALL CALC_GS(             CALL CALC_GS(
629       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
# Line 556  C        and step forward storing result Line 650  C        and step forward storing result
650       I         myIter,myThid)       I         myIter,myThid)
651           ENDIF           ENDIF
652  #endif  #endif
653    #ifdef ALLOW_PTRACERS
654             IF ( usePTRACERS ) THEN
655               CALL PTRACERS_INTEGERATE(
656         I         bi,bj,k,
657         I         xA,yA,uTrans,vTrans,rTrans,maskUp,
658         X         KappaRS,
659         I         myIter,myTime,myThid)
660             ENDIF
661    #endif /* ALLOW_PTRACERS */
662    
663  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
664  C--      Apply open boundary conditions  C--      Apply open boundary conditions
# Line 565  C--      Apply open boundary conditions Line 668  C--      Apply open boundary conditions
668  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
669    
670  C--      Freeze water  C--      Freeze water
671           IF (allowFreezing) THEN           IF ( allowFreezing .AND. .NOT. useSEAICE ) THEN
672  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
673  CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k  CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k
674  CADJ &   , key = kkey, byte = isbyte  CADJ &   , key = kkey, byte = isbyte
# Line 576  CADJ &   , key = kkey, byte = isbyte Line 679  CADJ &   , key = kkey, byte = isbyte
679  C--     end of thermodynamic k loop (Nr:1)  C--     end of thermodynamic k loop (Nr:1)
680          ENDDO          ENDDO
681    
682    cswdice -- add ---
683    #ifdef ALLOW_THERM_SEAICE
684    c timeaveraging for ice model values
685               CALL ICE_AVE(bi,bj,iMin,iMax,jMin,jMax,myThid )
686    #endif
687    cswdice --- end add ---
688    
689    
690    
 #ifdef ALLOW_AUTODIFF_TAMC  
 C? Patrick? What about this one?  
 cph Keys iikey and idkey dont seem to be needed  
 cph since storing occurs on different tape for each  
 cph impldiff call anyways.  
 cph Thus, common block comlev1_impl isnt needed either.  
 cph Storing below needed in the case useGMREDI.  
         iikey = (ikey-1)*maximpl  
 #endif /* ALLOW_AUTODIFF_TAMC */  
691    
692  C--     Implicit diffusion  C--     Implicit diffusion
693          IF (implicitDiffusion) THEN          IF (implicitDiffusion) THEN
694    
695           IF (tempStepping) THEN           IF (tempStepping) THEN
696  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
697              idkey = iikey + 1  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
 CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  
698  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
699              CALL IMPLDIFF(              CALL IMPLDIFF(
700       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
# Line 604  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bib Line 705  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bib
705    
706           IF (saltStepping) THEN           IF (saltStepping) THEN
707  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
708           idkey = iikey + 2  CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
 CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  
709  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
710              CALL IMPLDIFF(              CALL IMPLDIFF(
711       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
# Line 617  CADJ STORE gS(:,:,:,bi,bj) = comlev1_bib Line 717  CADJ STORE gS(:,:,:,bi,bj) = comlev1_bib
717  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
718           IF (tr1Stepping) THEN           IF (tr1Stepping) THEN
719  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
720  CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
721  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
722            CALL IMPLDIFF(            CALL IMPLDIFF(
723       I      bi, bj, iMin, iMax, jMin, jMax,       I      bi, bj, iMin, iMax, jMin, jMax,
# Line 627  CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_b Line 727  CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_b
727           ENDIF           ENDIF
728  #endif  #endif
729    
730    #ifdef ALLOW_PTRACERS
731    C Vertical diffusion (implicit) for passive tracers
732             IF ( usePTRACERS ) THEN
733               CALL PTRACERS_IMPLDIFF( bi,bj,KappaRS,myThid )
734             ENDIF
735    #endif /* ALLOW_PTRACERS */
736    
737  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
738  C--      Apply open boundary conditions  C--      Apply open boundary conditions
739           IF (useOBCS) THEN           IF (useOBCS) THEN
# Line 639  C--      Apply open boundary conditions Line 746  C--      Apply open boundary conditions
746  C--     End If implicitDiffusion  C--     End If implicitDiffusion
747          ENDIF          ENDIF
748    
749  Ccs-  #ifdef ALLOW_TIMEAVE
750            IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN
751              CALL TIMEAVE_CUMUL_1T(ConvectCountTave, ConvectCount,
752         I                           Nr, deltaTclock, bi, bj, myThid)
753            ENDIF
754            useVariableK = useKPP .OR. useGMredi .OR. ivdc_kappa.NE.0.
755            IF (taveFreq.GT.0. .AND. useVariableK ) THEN
756             IF (implicitDiffusion) THEN
757              CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, gT, kappaRT,
758         I                        Nr, 3, deltaTclock, bi, bj, myThid)
759             ELSE
760              CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, theta, kappaRT,
761         I                        Nr, 3, deltaTclock, bi, bj, myThid)
762             ENDIF
763            ENDIF
764    #endif /* ALLOW_TIMEAVE */
765    
766    #endif /* SINGLE_LAYER_MODE */
767    
768    C--   end bi,bj loops.
769         ENDDO         ENDDO
770        ENDDO        ENDDO
771    
772  #ifdef ALLOW_AIM  #ifndef DISABLE_DEBUGMODE
773        IF ( useAIM ) THEN        If (debugMode) THEN
774         CALL AIM_AIM2DYN_EXCHANGES( myTime, myThid )         CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)
775           CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)
776           CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)
777           CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)
778           CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)
779           CALL DEBUG_STATS_RL(Nr,Gt,'Gt (THERMODYNAMICS)',myThid)
780           CALL DEBUG_STATS_RL(Nr,Gs,'Gs (THERMODYNAMICS)',myThid)
781           CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (THERMODYNAMICS)',myThid)
782           CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (THERMODYNAMICS)',myThid)
783    #ifdef ALLOW_PTRACERS
784           IF ( usePTRACERS ) THEN
785             CALL PTRACERS_DEBUG(myThid)
786           ENDIF
787    #endif /* ALLOW_PTRACERS */
788        ENDIF        ENDIF
789         _EXCH_XYZ_R8(gT,myThid)  #endif
        _EXCH_XYZ_R8(gS,myThid)  
 #else  
       IF (staggerTimeStep.AND.useCubedSphereExchange) THEN  
        _EXCH_XYZ_R8(gT,myThid)  
        _EXCH_XYZ_R8(gS,myThid)  
       ENDIF  
 #endif /* ALLOW_AIM */  
790    
791        RETURN        RETURN
792        END        END
Legend:
Removed from v.1.11  
changed lines
  Added in v.1.39
  ViewVC Help
Powered by ViewVC 1.1.22