/[MITgcm]/MITgcm/pkg/generic_advdiff/gad_advection.F
ViewVC logotype

Diff of /MITgcm/pkg/generic_advdiff/gad_advection.F

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

revision 1.71 by jmc, Fri Mar 9 20:13:44 2012 UTC revision 1.72 by jmc, Mon Mar 4 18:20:46 2013 UTC
# Line 83  C !OUTPUT PARAMETERS: ================== Line 83  C !OUTPUT PARAMETERS: ==================
83  C  gTracer           :: tendency array  C  gTracer           :: tendency array
84        _RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
85    
86    C !FUNCTIONS: ==========================================================
87    #ifdef ALLOW_DIAGNOSTICS
88          CHARACTER*4 GAD_DIAG_SUFX
89          EXTERNAL    GAD_DIAG_SUFX
90          LOGICAL  DIAGNOSTICS_IS_ON
91          EXTERNAL DIAGNOSTICS_IS_ON
92    #endif
93    
94  C !LOCAL VARIABLES: ====================================================  C !LOCAL VARIABLES: ====================================================
95  C  maskUp        :: 2-D array for mask at W points  C  maskUp        :: 2-D array for mask at W points
96  C  maskLocW      :: 2-D array for mask at West points  C  maskLocW      :: 2-D array for mask at West points
# Line 98  C  uFld,vFld     :: 2-D local copy of ho Line 106  C  uFld,vFld     :: 2-D local copy of ho
106  C  wFld          :: 2-D local copy of vertical velocity  C  wFld          :: 2-D local copy of vertical velocity
107  C  uTrans,vTrans :: 2-D arrays of volume transports at U,V points  C  uTrans,vTrans :: 2-D arrays of volume transports at U,V points
108  C  rTrans        :: 2-D arrays of volume transports at W points  C  rTrans        :: 2-D arrays of volume transports at W points
109  C  rTransKp1     :: vertical volume transport at interface k+1  C  rTransKp      :: vertical volume transport at interface k+1
110  C  af            :: 2-D array for horizontal advective flux  C  af            :: 2-D array for horizontal advective flux
111  C  afx           :: 2-D array for horizontal advective flux, x direction  C  afx           :: 2-D array for horizontal advective flux, x direction
112  C  afy           :: 2-D array for horizontal advective flux, y direction  C  afy           :: 2-D array for horizontal advective flux, y direction
113  C  fVerT         :: 2 1/2D arrays for vertical advective flux  C  fVerT         :: 2 1/2D arrays for vertical advective flux
114  C  localTij      :: 2-D array, temporary local copy of tracer fld  C  localTij      :: 2-D array, temporary local copy of tracer field
115  C  localTijk     :: 3-D array, temporary local copy of tracer fld  C  localT3d      :: 3-D array, temporary local copy of tracer field
116  C  kp1Msk        :: flag (0,1) for over-riding mask for W levels  C  kp1Msk        :: flag (0,1) for over-riding mask for W levels
117  C  calc_fluxes_X :: logical to indicate to calculate fluxes in X dir  C  calc_fluxes_X :: logical to indicate to calculate fluxes in X dir
118  C  calc_fluxes_Y :: logical to indicate to calculate fluxes in Y dir  C  calc_fluxes_Y :: logical to indicate to calculate fluxes in Y dir
# Line 129  c     _RS maskUp  (1-OLx:sNx+OLx,1-OLy:s Line 137  c     _RS maskUp  (1-OLx:sNx+OLx,1-OLy:s
137        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
138        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
139        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
140        _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTransKp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
141        _RL af      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL af      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
142        _RL afx     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL afx     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
143        _RL afy     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL afy     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
144        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
145        _RL localTij(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL localTij(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
146        _RL localTijk(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL localT3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
147    #ifdef GAD_MULTIDIM_COMPRESSIBLE
148          _RL tmpTrac
149          _RL localVol(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
150          _RL locVol3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
151    #endif
152        _RL kp1Msk        _RL kp1Msk
153        LOGICAL calc_fluxes_X, calc_fluxes_Y, withSigns        LOGICAL calc_fluxes_X, calc_fluxes_Y, withSigns
154        LOGICAL interiorOnly, overlapOnly        LOGICAL interiorOnly, overlapOnly
# Line 153  C     msgBuf     :: Informational/error Line 166  C     msgBuf     :: Informational/error
166        CHARACTER*8 diagName        CHARACTER*8 diagName
167        CHARACTER*4 diagSufx        CHARACTER*4 diagSufx
168        LOGICAL     doDiagAdvX, doDiagAdvY, doDiagAdvR        LOGICAL     doDiagAdvX, doDiagAdvY, doDiagAdvR
169  C-    Functions:  #endif /* ALLOW_DIAGNOSTICS */
       CHARACTER*4 GAD_DIAG_SUFX  
       EXTERNAL    GAD_DIAG_SUFX  
       LOGICAL  DIAGNOSTICS_IS_ON  
       EXTERNAL DIAGNOSTICS_IS_ON  
 #endif  
170  CEOP  CEOP
171    
172  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 199  C--   Set diagnostics flags and suffix f Line 207  C--   Set diagnostics flags and suffix f
207          diagName = 'ADVr'//diagSufx          diagName = 'ADVr'//diagSufx
208          doDiagAdvR = DIAGNOSTICS_IS_ON( diagName, myThid )          doDiagAdvR = DIAGNOSTICS_IS_ON( diagName, myThid )
209        ENDIF        ENDIF
210  #endif  #endif /* ALLOW_DIAGNOSTICS */
211    
212  C--   Set up work arrays with valid (i.e. not NaN) values  C--   Set up work arrays with valid (i.e. not NaN) values
213  C     These inital values do not alter the numerical results. They  C     These inital values do not alter the numerical results. They
# Line 215  C     uninitialised but inert locations. Line 223  C     uninitialised but inert locations.
223          rTrans(i,j)  = 0. _d 0          rTrans(i,j)  = 0. _d 0
224          fVerT(i,j,1) = 0. _d 0          fVerT(i,j,1) = 0. _d 0
225          fVerT(i,j,2) = 0. _d 0          fVerT(i,j,2) = 0. _d 0
226          rTransKp1(i,j)= 0. _d 0          rTransKp(i,j)= 0. _d 0
227  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
228    # ifdef GAD_MULTIDIM_COMPRESSIBLE
229            localVol(i,j) = 0. _d 0
230    # endif
231          localTij(i,j) = 0. _d 0          localTij(i,j) = 0. _d 0
232          wfld(i,j)    = 0. _d 0          wFld(i,j)    = 0. _d 0
233  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
234         ENDDO         ENDDO
235        ENDDO        ENDDO
236    
# Line 278  C--   Make local copy of tracer array an Line 289  C--   Make local copy of tracer array an
289        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
290         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
291           localTij(i,j) = tracer(i,j,k,bi,bj)           localTij(i,j) = tracer(i,j,k,bi,bj)
292    #ifdef GAD_MULTIDIM_COMPRESSIBLE
293             localVol(i,j) = rA(i,j,bi,bj)*deepFac2C(k)
294         &                  *rhoFacC(k)*drF(k)*hFacC(i,j,k,bi,bj)
295         &                 + ( oneRS - maskC(i,j,k,bi,bj) )
296    #endif
297  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
298           maskLocW(i,j) = _maskW(i,j,k,bi,bj)*maskInW(i,j,bi,bj)           maskLocW(i,j) = _maskW(i,j,k,bi,bj)*maskInW(i,j,bi,bj)
299           maskLocS(i,j) = _maskS(i,j,k,bi,bj)*maskInS(i,j,bi,bj)           maskLocS(i,j) = _maskS(i,j,k,bi,bj)*maskInS(i,j,bi,bj)
# Line 288  C--   Make local copy of tracer array an Line 304  C--   Make local copy of tracer array an
304         ENDDO         ENDDO
305        ENDDO        ENDDO
306    
 cph-exch2#ifndef ALLOW_AUTODIFF_TAMC  
307        IF (useCubedSphereExchange) THEN        IF (useCubedSphereExchange) THEN
308          withSigns = .FALSE.          withSigns = .FALSE.
309          CALL FILL_CS_CORNER_UV_RS(          CALL FILL_CS_CORNER_UV_RS(
310       &            withSigns, maskLocW,maskLocS, bi,bj, myThid )       &            withSigns, maskLocW,maskLocS, bi,bj, myThid )
311        ENDIF        ENDIF
 cph-exch2#endif  
312    
313  C--   Multiple passes for different directions on different tiles  C--   Multiple passes for different directions on different tiles
314  C--   For cube need one pass for each of red, green and blue axes.  C--   For cube need one pass for each of red, green and blue axes.
# Line 335  C-    not CubedSphere Line 349  C-    not CubedSphere
349    
350  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
351  C--   X direction  C--   X direction
352  C-     Advective flux in X  
353    #ifdef ALLOW_AUTODIFF_TAMC
354    C-     Always reset advective flux in X
355          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
356           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
357            af(i,j) = 0.            af(i,j) = 0.
358           ENDDO           ENDDO
359          ENDDO          ENDDO
 C  
 #ifdef ALLOW_AUTODIFF_TAMC  
360  # ifndef DISABLE_MULTIDIM_ADVECTION  # ifndef DISABLE_MULTIDIM_ADVECTION
361  CADJ STORE localTij(:,:)  =  CADJ STORE localTij(:,:)  =
362  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
# Line 350  CADJ STORE af(:,:)  = Line 364  CADJ STORE af(:,:)  =
364  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
365  # endif  # endif
366  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
367  C  
368        IF (calc_fluxes_X) THEN        IF (calc_fluxes_X) THEN
369    
370  C-     Do not compute fluxes if  C-     Do not compute fluxes if
# Line 364  C-     Internal exchange for calculation Line 378  C-     Internal exchange for calculation
378       &                              localTij, bi,bj, myThid )       &                              localTij, bi,bj, myThid )
379          ENDIF          ENDIF
380    
381  #ifdef ALLOW_AUTODIFF_TAMC  C-     Advective flux in X
382    #ifndef ALLOW_AUTODIFF_TAMC
383            DO j=1-OLy,sNy+OLy
384             DO i=1-OLx,sNx+OLx
385              af(i,j) = 0.
386             ENDDO
387            ENDDO
388    #else /* ALLOW_AUTODIFF_TAMC */
389  # ifndef DISABLE_MULTIDIM_ADVECTION  # ifndef DISABLE_MULTIDIM_ADVECTION
390  CADJ STORE localTij(:,:)  =  CADJ STORE localTij(:,:)  =
391  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
# Line 418  C-     Advective flux in X : done Line 439  C-     Advective flux in X : done
439    
440  C-     Update the local tracer field where needed:  C-     Update the local tracer field where needed:
441  C      use "maksInC" to prevent updating tracer field in OB regions  C      use "maksInC" to prevent updating tracer field in OB regions
442    #ifdef ALLOW_AUTODIFF_TAMC
443    # ifdef GAD_MULTIDIM_COMPRESSIBLE
444    CADJ STORE localVol(:,:)  =
445    CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
446    CADJ STORE localTij(:,:)  =
447    CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
448    # endif
449    #endif /* ALLOW_AUTODIFF_TAMC */
450    
451  C      update in overlap-Only  C      update in overlap-Only
452         IF ( overlapOnly ) THEN         IF ( overlapOnly ) THEN
# Line 431  C         in corner region) but safer to Line 460  C         in corner region) but safer to
460          IF ( S_edge ) THEN          IF ( S_edge ) THEN
461           DO j=1-OLy,0           DO j=1-OLy,0
462            DO i=iMinUpd,iMaxUpd            DO i=iMinUpd,iMaxUpd
463    #ifdef GAD_MULTIDIM_COMPRESSIBLE
464               tmpTrac = localTij(i,j)*localVol(i,j)
465         &      -deltaTLev(k)*( af(i+1,j) - af(i,j) )
466         &                   *maskInC(i,j,bi,bj)
467               localVol(i,j) = localVol(i,j)
468         &      -deltaTLev(k)*( uTrans(i+1,j) - uTrans(i,j) )
469         &                   *maskInC(i,j,bi,bj)
470               localTij(i,j) = tmpTrac/localVol(i,j)
471    #else /* GAD_MULTIDIM_COMPRESSIBLE */
472             localTij(i,j) = localTij(i,j)             localTij(i,j) = localTij(i,j)
473       &      -deltaTLev(k)*recip_rhoFacC(k)       &      -deltaTLev(k)*recip_rhoFacC(k)
474       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
# Line 438  C         in corner region) but safer to Line 476  C         in corner region) but safer to
476       &       *( af(i+1,j)-af(i,j)       &       *( af(i+1,j)-af(i,j)
477       &         -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))       &         -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))
478       &        )*maskInC(i,j,bi,bj)       &        )*maskInC(i,j,bi,bj)
479    #endif /* GAD_MULTIDIM_COMPRESSIBLE */
480            ENDDO            ENDDO
481           ENDDO           ENDDO
482          ENDIF          ENDIF
483          IF ( N_edge ) THEN          IF ( N_edge ) THEN
484           DO j=sNy+1,sNy+OLy           DO j=sNy+1,sNy+OLy
485            DO i=iMinUpd,iMaxUpd            DO i=iMinUpd,iMaxUpd
486    #ifdef GAD_MULTIDIM_COMPRESSIBLE
487               tmpTrac = localTij(i,j)*localVol(i,j)
488         &      -deltaTLev(k)*( af(i+1,j) - af(i,j) )
489         &                   *maskInC(i,j,bi,bj)
490               localVol(i,j) = localVol(i,j)
491         &      -deltaTLev(k)*( uTrans(i+1,j) - uTrans(i,j) )
492         &                   *maskInC(i,j,bi,bj)
493               localTij(i,j) = tmpTrac/localVol(i,j)
494    #else /* GAD_MULTIDIM_COMPRESSIBLE */
495             localTij(i,j) = localTij(i,j)             localTij(i,j) = localTij(i,j)
496       &      -deltaTLev(k)*recip_rhoFacC(k)       &      -deltaTLev(k)*recip_rhoFacC(k)
497       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
# Line 451  C         in corner region) but safer to Line 499  C         in corner region) but safer to
499       &       *( af(i+1,j)-af(i,j)       &       *( af(i+1,j)-af(i,j)
500       &         -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))       &         -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))
501       &        )*maskInC(i,j,bi,bj)       &        )*maskInC(i,j,bi,bj)
502    #endif /* GAD_MULTIDIM_COMPRESSIBLE */
503            ENDDO            ENDDO
504           ENDDO           ENDDO
505          ENDIF          ENDIF
# Line 463  C      do not only update the overlap Line 512  C      do not only update the overlap
512          IF ( interiorOnly .AND. N_edge ) jMaxUpd = sNy          IF ( interiorOnly .AND. N_edge ) jMaxUpd = sNy
513          DO j=jMinUpd,jMaxUpd          DO j=jMinUpd,jMaxUpd
514           DO i=1-OLx+1,sNx+OLx-1           DO i=1-OLx+1,sNx+OLx-1
515    #ifdef GAD_MULTIDIM_COMPRESSIBLE
516               tmpTrac = localTij(i,j)*localVol(i,j)
517         &      -deltaTLev(k)*( af(i+1,j) - af(i,j) )
518         &                   *maskInC(i,j,bi,bj)
519               localVol(i,j) = localVol(i,j)
520         &      -deltaTLev(k)*( uTrans(i+1,j) - uTrans(i,j) )
521         &                   *maskInC(i,j,bi,bj)
522               localTij(i,j) = tmpTrac/localVol(i,j)
523    #else /* GAD_MULTIDIM_COMPRESSIBLE */
524             localTij(i,j) = localTij(i,j)             localTij(i,j) = localTij(i,j)
525       &      -deltaTLev(k)*recip_rhoFacC(k)       &      -deltaTLev(k)*recip_rhoFacC(k)
526       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
# Line 470  C      do not only update the overlap Line 528  C      do not only update the overlap
528       &       *( af(i+1,j)-af(i,j)       &       *( af(i+1,j)-af(i,j)
529       &         -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))       &         -tracer(i,j,k,bi,bj)*(uTrans(i+1,j)-uTrans(i,j))
530       &        )*maskInC(i,j,bi,bj)       &        )*maskInC(i,j,bi,bj)
531    #endif /* GAD_MULTIDIM_COMPRESSIBLE */
532           ENDDO           ENDDO
533          ENDDO          ENDDO
534  C-      keep advective flux (for diagnostics)  C-      keep advective flux (for diagnostics)
# Line 487  C--   End of X direction Line 546  C--   End of X direction
546    
547  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
548  C--   Y direction  C--   Y direction
549  cph-test  
550  C-     Advective flux in Y  #ifdef ALLOW_AUTODIFF_TAMC
551    C-     Always reset advective flux in Y
552          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
553           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
554            af(i,j) = 0.            af(i,j) = 0.
555           ENDDO           ENDDO
556          ENDDO          ENDDO
 C  
 #ifdef ALLOW_AUTODIFF_TAMC  
557  # ifndef DISABLE_MULTIDIM_ADVECTION  # ifndef DISABLE_MULTIDIM_ADVECTION
558  CADJ STORE localTij(:,:)  =  CADJ STORE localTij(:,:)  =
559  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
# Line 503  CADJ STORE af(:,:)  = Line 561  CADJ STORE af(:,:)  =
561  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
562  # endif  # endif
563  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
564  C  
565        IF (calc_fluxes_Y) THEN        IF (calc_fluxes_Y) THEN
566    
567  C-     Do not compute fluxes if  C-     Do not compute fluxes if
# Line 518  C-     Internal exchange for calculation Line 576  C-     Internal exchange for calculation
576          ENDIF          ENDIF
577    
578  C-     Advective flux in Y  C-     Advective flux in Y
579    #ifndef ALLOW_AUTODIFF_TAMC
580          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
581           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
582            af(i,j) = 0.            af(i,j) = 0.
583           ENDDO           ENDDO
584          ENDDO          ENDDO
585    #else /* ALLOW_AUTODIFF_TAMC */
 #ifdef ALLOW_AUTODIFF_TAMC  
586  #ifndef DISABLE_MULTIDIM_ADVECTION  #ifndef DISABLE_MULTIDIM_ADVECTION
587  CADJ STORE localTij(:,:)  =  CADJ STORE localTij(:,:)  =
588  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
# Line 578  C-     Advective flux in Y : done Line 636  C-     Advective flux in Y : done
636    
637  C-     Update the local tracer field where needed:  C-     Update the local tracer field where needed:
638  C      use "maksInC" to prevent updating tracer field in OB regions  C      use "maksInC" to prevent updating tracer field in OB regions
639    #ifdef ALLOW_AUTODIFF_TAMC
640    # ifdef GAD_MULTIDIM_COMPRESSIBLE
641    CADJ STORE localVol(:,:)  =
642    CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
643    CADJ STORE localTij(:,:)  =
644    CADJ &     comlev1_bibj_k_gad_pass, key=passkey, kind=isbyte
645    # endif
646    #endif /* ALLOW_AUTODIFF_TAMC */
647    
648  C      update in overlap-Only  C      update in overlap-Only
649         IF ( overlapOnly ) THEN         IF ( overlapOnly ) THEN
# Line 591  C         in corner region) but safer to Line 657  C         in corner region) but safer to
657          IF ( W_edge ) THEN          IF ( W_edge ) THEN
658           DO j=jMinUpd,jMaxUpd           DO j=jMinUpd,jMaxUpd
659            DO i=1-OLx,0            DO i=1-OLx,0
660    #ifdef GAD_MULTIDIM_COMPRESSIBLE
661               tmpTrac = localTij(i,j)*localVol(i,j)
662         &      -deltaTLev(k)*( af(i,j+1) - af(i,j) )
663         &                   *maskInC(i,j,bi,bj)
664               localVol(i,j) = localVol(i,j)
665         &      -deltaTLev(k)*( vTrans(i,j+1) - vTrans(i,j) )
666         &                   *maskInC(i,j,bi,bj)
667               localTij(i,j) = tmpTrac/localVol(i,j)
668    #else /* GAD_MULTIDIM_COMPRESSIBLE */
669             localTij(i,j) = localTij(i,j)             localTij(i,j) = localTij(i,j)
670       &      -deltaTLev(k)*recip_rhoFacC(k)       &      -deltaTLev(k)*recip_rhoFacC(k)
671       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
# Line 598  C         in corner region) but safer to Line 673  C         in corner region) but safer to
673       &       *( af(i,j+1)-af(i,j)       &       *( af(i,j+1)-af(i,j)
674       &         -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))       &         -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))
675       &        )*maskInC(i,j,bi,bj)       &        )*maskInC(i,j,bi,bj)
676    #endif /* GAD_MULTIDIM_COMPRESSIBLE */
677            ENDDO            ENDDO
678           ENDDO           ENDDO
679          ENDIF          ENDIF
680          IF ( E_edge ) THEN          IF ( E_edge ) THEN
681           DO j=jMinUpd,jMaxUpd           DO j=jMinUpd,jMaxUpd
682            DO i=sNx+1,sNx+OLx            DO i=sNx+1,sNx+OLx
683    #ifdef GAD_MULTIDIM_COMPRESSIBLE
684               tmpTrac = localTij(i,j)*localVol(i,j)
685         &      -deltaTLev(k)*( af(i,j+1) - af(i,j) )
686         &                   *maskInC(i,j,bi,bj)
687               localVol(i,j) = localVol(i,j)
688         &      -deltaTLev(k)*( vTrans(i,j+1) - vTrans(i,j) )
689         &                   *maskInC(i,j,bi,bj)
690               localTij(i,j) = tmpTrac/localVol(i,j)
691    #else /* GAD_MULTIDIM_COMPRESSIBLE */
692             localTij(i,j) = localTij(i,j)             localTij(i,j) = localTij(i,j)
693       &      -deltaTLev(k)*recip_rhoFacC(k)       &      -deltaTLev(k)*recip_rhoFacC(k)
694       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
# Line 611  C         in corner region) but safer to Line 696  C         in corner region) but safer to
696       &       *( af(i,j+1)-af(i,j)       &       *( af(i,j+1)-af(i,j)
697       &         -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))       &         -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))
698       &        )*maskInC(i,j,bi,bj)       &        )*maskInC(i,j,bi,bj)
699    #endif /* GAD_MULTIDIM_COMPRESSIBLE */
700            ENDDO            ENDDO
701           ENDDO           ENDDO
702          ENDIF          ENDIF
# Line 623  C      do not only update the overlap Line 709  C      do not only update the overlap
709          IF ( interiorOnly .AND. E_edge ) iMaxUpd = sNx          IF ( interiorOnly .AND. E_edge ) iMaxUpd = sNx
710          DO j=1-OLy+1,sNy+OLy-1          DO j=1-OLy+1,sNy+OLy-1
711           DO i=iMinUpd,iMaxUpd           DO i=iMinUpd,iMaxUpd
712    #ifdef GAD_MULTIDIM_COMPRESSIBLE
713               tmpTrac = localTij(i,j)*localVol(i,j)
714         &      -deltaTLev(k)*( af(i,j+1) - af(i,j) )
715         &                   *maskInC(i,j,bi,bj)
716               localVol(i,j) = localVol(i,j)
717         &      -deltaTLev(k)*( vTrans(i,j+1) - vTrans(i,j) )
718         &                   *maskInC(i,j,bi,bj)
719               localTij(i,j) = tmpTrac/localVol(i,j)
720    #else /* GAD_MULTIDIM_COMPRESSIBLE */
721             localTij(i,j) = localTij(i,j)             localTij(i,j) = localTij(i,j)
722       &      -deltaTLev(k)*recip_rhoFacC(k)       &      -deltaTLev(k)*recip_rhoFacC(k)
723       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
# Line 630  C      do not only update the overlap Line 725  C      do not only update the overlap
725       &       *( af(i,j+1)-af(i,j)       &       *( af(i,j+1)-af(i,j)
726       &         -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))       &         -tracer(i,j,k,bi,bj)*(vTrans(i,j+1)-vTrans(i,j))
727       &        )*maskInC(i,j,bi,bj)       &        )*maskInC(i,j,bi,bj)
728    #endif /* GAD_MULTIDIM_COMPRESSIBLE */
729           ENDDO           ENDDO
730          ENDDO          ENDDO
731  C-      keep advective flux (for diagnostics)  C-      keep advective flux (for diagnostics)
# Line 650  C--   End of ipass loop Line 746  C--   End of ipass loop
746    
747        IF ( implicitAdvection ) THEN        IF ( implicitAdvection ) THEN
748  C-    explicit advection is done ; store tendency in gTracer:  C-    explicit advection is done ; store tendency in gTracer:
749    #ifdef GAD_MULTIDIM_COMPRESSIBLE
750            STOP 'GAD_ADVECTION: missing code for implicitAdvection'
751    #endif /* GAD_MULTIDIM_COMPRESSIBLE */
752          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
753           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
754            gTracer(i,j,k,bi,bj)=            gTracer(i,j,k,bi,bj)=
# Line 660  C-    explicit advection is done ; store Line 759  C-    explicit advection is done ; store
759  C-    horizontal advection done; store intermediate result in 3D array:  C-    horizontal advection done; store intermediate result in 3D array:
760          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
761           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
762            localTijk(i,j,k) = localTij(i,j)  #ifdef GAD_MULTIDIM_COMPRESSIBLE
763              locVol3d(i,j,k) = localVol(i,j)
764    #endif /* GAD_MULTIDIM_COMPRESSIBLE */
765              localT3d(i,j,k) = localTij(i,j)
766           ENDDO           ENDDO
767          ENDDO          ENDDO
768        ENDIF        ENDIF
# Line 674  C-    horizontal advection done; store i Line 776  C-    horizontal advection done; store i
776            diagName = 'ADVy'//diagSufx            diagName = 'ADVy'//diagSufx
777            CALL DIAGNOSTICS_FILL( afy, diagName, k,1, 2,bi,bj, myThid )            CALL DIAGNOSTICS_FILL( afy, diagName, k,1, 2,bi,bj, myThid )
778          ENDIF          ENDIF
779  #endif  #endif /* ALLOW_DIAGNOSTICS */
780    
781  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
782        IF ( debugLevel .GE. debLevC        IF ( debugLevel .GE. debLevC
# Line 709  c       kp1=min(Nr,k+1) Line 811  c       kp1=min(Nr,k+1)
811  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
812  CADJ STORE rtrans(:,:)  =  CADJ STORE rtrans(:,:)  =
813  CADJ &     comlev1_bibj_k_gad, key=kkey, kind=isbyte  CADJ &     comlev1_bibj_k_gad, key=kkey, kind=isbyte
814  cphCADJ STORE wfld(:,:)  =  cphCADJ STORE wFld(:,:)  =
815  cphCADJ &     comlev1_bibj_k_gad, key=kkey, kind=isbyte  cphCADJ &     comlev1_bibj_k_gad, key=kkey, kind=isbyte
816  #endif  #endif
817    
# Line 724  C- a hack to prevent Water-Vapor vert.tr Line 826  C- a hack to prevent Water-Vapor vert.tr
826  #endif  #endif
827    
828  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
829  cphmultiCADJ STORE wfld(:,:)  =  cphmultiCADJ STORE wFld(:,:)  =
830  cphmultiCADJ &     comlev1_bibj_k_gad, key=kkey, kind=isbyte  cphmultiCADJ &     comlev1_bibj_k_gad, key=kkey, kind=isbyte
831  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
832    
833  C- Surface interface :  C- Surface interface :
834           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
835            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
836             rTransKp1(i,j) = kp1Msk*rTrans(i,j)             rTransKp(i,j) = kp1Msk*rTrans(i,j)
837             wFld(i,j)   = 0.             wFld(i,j)   = 0.
838             rTrans(i,j) = 0.             rTrans(i,j) = 0.
839             fVerT(i,j,kUp) = 0.             fVerT(i,j,kUp) = 0.
# Line 741  C- Surface interface : Line 843  C- Surface interface :
843          ELSE          ELSE
844    
845  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
846  cphmultiCADJ STORE wfld(:,:)  =  cphmultiCADJ STORE wFld(:,:)  =
847  cphmultiCADJ &     comlev1_bibj_k_gad, key=kkey, kind=isbyte  cphmultiCADJ &     comlev1_bibj_k_gad, key=kkey, kind=isbyte
848  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
849    
850  C- Interior interface :  C- Interior interface :
851           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
852            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
853             rTransKp1(i,j) = kp1Msk*rTrans(i,j)             rTransKp(i,j) = kp1Msk*rTrans(i,j)
854             wFld(i,j)   = wVel(i,j,k,bi,bj)             wFld(i,j)   = wVel(i,j,k,bi,bj)
855             rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)             rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)
856       &                 *deepFac2F(k)*rhoFacF(k)       &                 *deepFac2F(k)*rhoFacF(k)
# Line 766  C--   Residual transp = Bolus transp + E Line 868  C--   Residual transp = Bolus transp + E
868  #endif /* ALLOW_GMREDI */  #endif /* ALLOW_GMREDI */
869    
870  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
871  cphmultiCADJ STORE localTijk(:,:,k)  cphmultiCADJ STORE localT3d(:,:,k)
872  cphmultiCADJ &     = comlev1_bibj_k_gad, key=kkey, kind=isbyte  cphmultiCADJ &     = comlev1_bibj_k_gad, key=kkey, kind=isbyte
873  cphmultiCADJ STORE rTrans(:,:)  cphmultiCADJ STORE rTrans(:,:)
874  cphmultiCADJ &     = comlev1_bibj_k_gad, key=kkey, kind=isbyte  cphmultiCADJ &     = comlev1_bibj_k_gad, key=kkey, kind=isbyte
# Line 776  C-    Compute vertical advective flux in Line 878  C-    Compute vertical advective flux in
878           IF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST           IF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST
879       &      .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN       &      .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN
880             CALL GAD_DST2U1_ADV_R(    bi,bj,k, advectionScheme,             CALL GAD_DST2U1_ADV_R(    bi,bj,k, advectionScheme,
881       I                          deltaTLev(k),rTrans,wFld,localTijk,       I                          deltaTLev(k),rTrans,wFld,localT3d,
882       O                               fVerT(1-OLx,1-OLy,kUp), myThid )       O                               fVerT(1-OLx,1-OLy,kUp), myThid )
883           ELSEIF( vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN           ELSEIF( vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
884             CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k, deltaTLev(k),             CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k, deltaTLev(k),
885       I                               rTrans, wFld, localTijk,       I                               rTrans, wFld, localT3d,
886       O                               fVerT(1-OLx,1-OLy,kUp), myThid )       O                               fVerT(1-OLx,1-OLy,kUp), myThid )
887           ELSEIF( vertAdvecScheme.EQ.ENUM_DST3 ) THEN           ELSEIF( vertAdvecScheme.EQ.ENUM_DST3 ) THEN
888             CALL GAD_DST3_ADV_R(      bi,bj,k, deltaTLev(k),             CALL GAD_DST3_ADV_R(      bi,bj,k, deltaTLev(k),
889       I                               rTrans, wFld, localTijk,       I                               rTrans, wFld, localT3d,
890       O                               fVerT(1-OLx,1-OLy,kUp), myThid )       O                               fVerT(1-OLx,1-OLy,kUp), myThid )
891           ELSEIF( vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN           ELSEIF( vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
892             CALL GAD_DST3FL_ADV_R(    bi,bj,k, deltaTLev(k),             CALL GAD_DST3FL_ADV_R(    bi,bj,k, deltaTLev(k),
893       I                               rTrans, wFld, localTijk,       I                               rTrans, wFld, localT3d,
894       O                               fVerT(1-OLx,1-OLy,kUp), myThid )       O                               fVerT(1-OLx,1-OLy,kUp), myThid )
895  #ifndef ALLOW_AUTODIFF_TAMC  #ifndef ALLOW_AUTODIFF_TAMC
896           ELSEIF (vertAdvecScheme.EQ.ENUM_OS7MP ) THEN           ELSEIF (vertAdvecScheme.EQ.ENUM_OS7MP ) THEN
897             CALL GAD_OS7MP_ADV_R(     bi,bj,k, deltaTLev(k),             CALL GAD_OS7MP_ADV_R(     bi,bj,k, deltaTLev(k),
898       I                               rTrans, wFld, localTijk,       I                               rTrans, wFld, localT3d,
899       O                               fVerT(1-OLx,1-OLy,kUp), myThid )       O                               fVerT(1-OLx,1-OLy,kUp), myThid )
900  #endif  #endif
901           ELSE           ELSE
# Line 806  C- end Surface/Interior if bloc Line 908  C- end Surface/Interior if bloc
908  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
909  cphmultiCADJ STORE rTrans(:,:)  cphmultiCADJ STORE rTrans(:,:)
910  cphmultiCADJ &     = comlev1_bibj_k_gad, key=kkey, kind=isbyte  cphmultiCADJ &     = comlev1_bibj_k_gad, key=kkey, kind=isbyte
911  cphmultiCADJ STORE rTranskp1(:,:)  cphmultiCADJ STORE rTranskp(:,:)
912  cphmultiCADJ &     = comlev1_bibj_k_gad, key=kkey, kind=isbyte  cphmultiCADJ &     = comlev1_bibj_k_gad, key=kkey, kind=isbyte
913  cph --- following storing of fVerT is critical for correct  cph --- following storing of fVerT is critical for correct
914  cph --- gradient with multiDimAdvection  cph --- gradient with multiDimAdvection
# Line 817  CADJ &     = comlev1_bibj_k_gad, key=kke Line 919  CADJ &     = comlev1_bibj_k_gad, key=kke
919  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
920    
921  C--   Divergence of vertical fluxes  C--   Divergence of vertical fluxes
922    #ifdef GAD_MULTIDIM_COMPRESSIBLE
923          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
924           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
925            localTij(i,j) = localTijk(i,j,k)            tmpTrac = localT3d(i,j,k)*locVol3d(i,j,k)
926         &      -deltaTLev(k)*( fVerT(i,j,kDown)-fVerT(i,j,kUp) )
927         &                   *rkSign*maskInC(i,j,bi,bj)
928              localVol(i,j) = locVol3d(i,j,k)
929         &      -deltaTLev(k)*( rTransKp(i,j) - rTrans(i,j) )
930         &                   *rkSign*maskInC(i,j,bi,bj)
931    C- localTij only needed for Variance Bugget: can be move there
932              localTij(i,j) = tmpTrac/localVol(i,j)
933    C--  without rescaling of tendencies:
934    c         gTracer(i,j,k,bi,bj)=
935    c    &     (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTLev(k)
936    C--  Non-Lin Free-Surf: consistent with rescaling of tendencies
937    C     (in FREESURF_RESCALE_G) and RealFreshFlux/addMass.
938    C    Also valid for linear Free-Surf (r & r* coords) w/wout RealFreshFlux
939    C     and consistent with linFSConserveTr and "surfExpan_" monitor.
940              gTracer(i,j,k,bi,bj) =
941         &          ( tmpTrac - tracer(i,j,k,bi,bj)*localVol(i,j) )
942         &            *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
943         &            *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
944         &            *recip_rhoFacC(k)
945         &            /deltaTLev(k)
946             ENDDO
947            ENDDO
948    #else /* GAD_MULTIDIM_COMPRESSIBLE */
949            DO j=1-OLy,sNy+OLy
950             DO i=1-OLx,sNx+OLx
951              localTij(i,j) = localT3d(i,j,k)
952       &      -deltaTLev(k)*recip_rhoFacC(k)       &      -deltaTLev(k)*recip_rhoFacC(k)
953       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
954       &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)       &       *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)
955       &       *( fVerT(i,j,kDown)-fVerT(i,j,kUp)       &       *( fVerT(i,j,kDown)-fVerT(i,j,kUp)
956       &         -tracer(i,j,k,bi,bj)*(rTransKp1(i,j)-rTrans(i,j))       &         -tracer(i,j,k,bi,bj)*(rTransKp(i,j)-rTrans(i,j))
957       &        )*rkSign*maskInC(i,j,bi,bj)       &        )*rkSign*maskInC(i,j,bi,bj)
958            gTracer(i,j,k,bi,bj)=            gTracer(i,j,k,bi,bj)=
959       &     (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTLev(k)       &     (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTLev(k)
960           ENDDO           ENDDO
961          ENDDO          ENDDO
962    #endif /* GAD_MULTIDIM_COMPRESSIBLE */
963    
964  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
965          IF ( doDiagAdvR ) THEN          IF ( doDiagAdvR ) THEN
# Line 837  C--   Divergence of vertical fluxes Line 967  C--   Divergence of vertical fluxes
967            CALL DIAGNOSTICS_FILL( fVerT(1-OLx,1-OLy,kUp),            CALL DIAGNOSTICS_FILL( fVerT(1-OLx,1-OLy,kUp),
968       &                           diagName, k,1, 2,bi,bj, myThid )       &                           diagName, k,1, 2,bi,bj, myThid )
969          ENDIF          ENDIF
970  #endif  #endif /* ALLOW_DIAGNOSTICS */
971    
972  C--   End of K loop for vertical flux  C--   End of K loop for vertical flux
973         ENDDO         ENDDO
Legend:
Removed from v.1.71  
changed lines
  Added in v.1.72
  ViewVC Help
Powered by ViewVC 1.1.22