/[MITgcm]/MITgcm/pkg/seaice/seaice_advection.F
ViewVC logotype

Diff of /MITgcm/pkg/seaice/seaice_advection.F

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

revision 1.3 by jmc, Mon Jun 19 15:48:35 2006 UTC revision 1.4 by jmc, Wed Nov 1 01:56:23 2006 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4  #include "GAD_OPTIONS.h"  #include "SEAICE_OPTIONS.h"
5    #ifdef ALLOW_GENERIC_ADVDIFF
6    # include "GAD_OPTIONS.h"
7    #endif
8    
9  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
10  CBOP  CBOP
# Line 9  C !ROUTINE: SEAICE_ADVECTION Line 12  C !ROUTINE: SEAICE_ADVECTION
12    
13  C !INTERFACE: ==========================================================  C !INTERFACE: ==========================================================
14        SUBROUTINE SEAICE_ADVECTION(        SUBROUTINE SEAICE_ADVECTION(
      I     advectionScheme,  
15       I     tracerIdentity,       I     tracerIdentity,
16       I     uVel, vVel, tracer,       I     advectionScheme,
17       O     gTracer,       I     extensiveFld,
18         I     uFld, vFld, uTrans, vTrans, iceFld, r_hFld,
19         O     gFld, afx, afy,
20       I     bi, bj, myTime, myIter, myThid)       I     bi, bj, myTime, myIter, myThid)
21    
22  C !DESCRIPTION:  C !DESCRIPTION:
23  C Calculates the tendency of a tracer due to advection.  C Calculates the tendency of a sea-ice field due to advection.
24  C It uses the multi-dimensional method given in \ref{sect:multiDimAdvection}  C It uses the multi-dimensional method given in \ref{sect:multiDimAdvection}
25  C and can only be used for the non-linear advection schemes such as the  C and can only be used for the non-linear advection schemes such as the
26  C direct-space-time method and flux-limiters.  C direct-space-time method and flux-limiters.
27  C  C
28  C This routine is an adaption of the GAD_ADVECTION for 2D-fields.  C This routine is an adaption of the GAD_ADVECTION for 2D-fields.
29  C Seaice velocities are not divergence free; therefore the contribution  C for Area, effective thickness or other "extensive" sea-ice field,
30  C tracer*div(u) that is present in gad_advection is removed in this routine.  C  the contribution iceFld*div(u) (that is present in gad_advection)
31    C  is not included here.
32  C  C
33  C The algorithm is as follows:  C The algorithm is as follows:
34  C \begin{itemize}  C \begin{itemize}
# Line 42  C !USES: =============================== Line 47  C !USES: ===============================
47  #include "EEPARAMS.h"  #include "EEPARAMS.h"
48  #include "PARAMS.h"  #include "PARAMS.h"
49  #include "GRID.h"  #include "GRID.h"
 #include "GAD.h"  
50  #include "SEAICE_PARAMS.h"  #include "SEAICE_PARAMS.h"
51    #ifdef ALLOW_GENERIC_ADVDIFF
52    # include "GAD.h"
53    #endif
54  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
55  # include "tamc.h"  # include "tamc.h"
56  # include "tamc_keys.h"  # include "tamc_keys.h"
# Line 54  C !USES: =============================== Line 61  C !USES: ===============================
61  #endif /* ALLOW_EXCH2 */  #endif /* ALLOW_EXCH2 */
62    
63  C !INPUT PARAMETERS: ===================================================  C !INPUT PARAMETERS: ===================================================
64  C  implicitAdvection :: implicit vertical advection (later on)  C  tracerIdentity  :: tracer identifier (required only for OBCS)
65  C  advectionScheme   :: advection scheme to use (Horizontal plane)  C  advectionScheme :: advection scheme to use (Horizontal plane)
66  C  tracerIdentity    :: tracer identifier (required only for OBCS)  C  extensiveFld    :: indicates to advect an "extensive" type of ice field
67  C  uVel              :: velocity, zonal component  C  uFld            :: velocity, zonal component
68  C  vVel              :: velocity, meridional component  C  vFld            :: velocity, meridional component
69  C  tracer            :: tracer field  C  uTrans,vTrans   :: volume transports at U,V points
70  C  bi,bj             :: tile indices  C  iceFld          :: sea-ice field
71  C  myTime            :: current time  C  r_hFld          :: reciprol of ice-thickness (only used for "intensive"
72  C  myIter            :: iteration number  C                     type of sea-ice field)
73  C  myThid            :: thread number  C  bi,bj           :: tile indices
74        INTEGER advectionScheme  C  myTime          :: current time
75    C  myIter          :: iteration number
76    C  myThid          :: thread number
77        INTEGER tracerIdentity        INTEGER tracerIdentity
78        _RL uVel  (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)        INTEGER advectionScheme
79        _RL vVel  (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)        LOGICAL extensiveFld
80        _RL tracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)        _RL uFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
81          _RL vFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
82          _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
83          _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
84          _RL iceFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
85          _RL r_hFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
86        INTEGER bi,bj        INTEGER bi,bj
87        _RL myTime        _RL myTime
88        INTEGER myIter        INTEGER myIter
89        INTEGER myThid        INTEGER myThid
90    
91  C !OUTPUT PARAMETERS: ==================================================  C !OUTPUT PARAMETERS: ==================================================
92  C  gTracer           :: tendancy array  C  gFld          :: tendency array
93        _RL gTracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)  C  afx           :: horizontal advective flux, x direction
94    C  afy           :: horizontal advective flux, y direction
95          _RL gFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
96          _RL afx   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
97          _RL afy   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98    
99  C !LOCAL VARIABLES: ====================================================  C !LOCAL VARIABLES: ====================================================
100  C  maskLocW      :: 2-D array for mask at West points  C  maskLocW      :: 2-D array for mask at West points
101  C  maskLocS      :: 2-D array for mask at South points  C  maskLocS      :: 2-D array for mask at South points
102  C  iMin,iMax,    :: loop range for called routines  C  iMin,iMax,    :: loop range for called routines
103  C  jMin,jMax     :: loop range for called routines  C  jMin,jMax     :: loop range for called routines
104  C [iMin,iMax]Upd :: loop range to update tracer field  C [iMin,iMax]Upd :: loop range to update sea-ice field
105  C [jMin,jMax]Upd :: loop range to update tracer field  C [jMin,jMax]Upd :: loop range to update sea-ice field
106  C  i,j,k         :: loop indices  C  i,j,k         :: loop indices
 C  kp1           :: =k+1 for k<Nr, =Nr for k=Nr  
 C  xA,yA         :: areas of X and Y face of tracer cells  
 C  uTrans,vTrans :: 2-D arrays of volume transports at U,V points  
107  C  af            :: 2-D array for horizontal advective flux  C  af            :: 2-D array for horizontal advective flux
108  C  afx           :: 2-D array for horizontal advective flux, x direction  C  localTij      :: 2-D array, temporary local copy of sea-ice fld
 C  afy           :: 2-D array for horizontal advective flux, y direction  
 C  fVerT         :: 2 1/2D arrays for vertical advective flux  
 C  localTij      :: 2-D array, temporary local copy of tracer fld  
 C  localTijk     :: 3-D array, temporary local copy of tracer fld  
 C  kp1Msk        :: flag (0,1) for over-riding mask for W levels  
109  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
110  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
111  C  interiorOnly  :: only update the interior of myTile, but not the edges  C  interiorOnly  :: only update the interior of myTile, but not the edges
112  C  overlapOnly   :: only update the edges of myTile, but not the interior  C  overlapOnly   :: only update the edges of myTile, but not the interior
113  C  nipass        :: number of passes in multi-dimensional method  C  nipass        :: number of passes in multi-dimensional method
114  C  ipass         :: number of the current pass being made  C  ipass         :: number of the current pass being made
115  C  myTile        :: variables used to determine which cube face  C  myTile        :: variables used to determine which cube face
116  C  nCFace        :: owns a tile for cube grid runs using  C  nCFace        :: owns a tile for cube grid runs using
117  C                :: multi-dim advection.  C                :: multi-dim advection.
118  C [N,S,E,W]_edge :: true if N,S,E,W edge of myTile is an Edge of the cube  C [N,S,E,W]_edge :: true if N,S,E,W edge of myTile is an Edge of the cube
# Line 111  C [N,S,E,W]_edge :: true if N,S,E,W edge Line 121  C [N,S,E,W]_edge :: true if N,S,E,W edge
121        INTEGER iMin,iMax,jMin,jMax        INTEGER iMin,iMax,jMin,jMax
122        INTEGER iMinUpd,iMaxUpd,jMinUpd,jMaxUpd        INTEGER iMinUpd,iMaxUpd,jMinUpd,jMaxUpd
123        INTEGER i,j,k        INTEGER i,j,k
       _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL uFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL vFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
124        _RL af      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL af      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RL afx     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL afy     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)  
125        _RL localTij(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL localTij(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
126        LOGICAL calc_fluxes_X, calc_fluxes_Y, withSigns        LOGICAL calc_fluxes_X, calc_fluxes_Y, withSigns
127        LOGICAL interiorOnly, overlapOnly        LOGICAL interiorOnly, overlapOnly
# Line 135  C [N,S,E,W]_edge :: true if N,S,E,W edge Line 136  C [N,S,E,W]_edge :: true if N,S,E,W edge
136        CHARACTER*4 GAD_DIAG_SUFX, diagSufx        CHARACTER*4 GAD_DIAG_SUFX, diagSufx
137        EXTERNAL    GAD_DIAG_SUFX        EXTERNAL    GAD_DIAG_SUFX
138  #endif  #endif
139          LOGICAL dBug
140          _RL tmpFac
141  CEOP  CEOP
142    
143  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 147  CEOP Line 150  CEOP
150            act3 = myThid - 1            act3 = myThid - 1
151            max3 = nTx*nTy            max3 = nTx*nTy
152            act4 = ikey_dynamics - 1            act4 = ikey_dynamics - 1
153            igadkey = (act0 + 1)            igadkey = (act0 + 1)
154       &                      + act1*max0       &                      + act1*max0
155       &                      + act2*max0*max1       &                      + act2*max0*max1
156       &                      + act3*max0*max1*max2       &                      + act3*max0*max1*max2
# Line 159  CEOP Line 162  CEOP
162  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
163    
164  CML#ifdef ALLOW_DIAGNOSTICS  CML#ifdef ALLOW_DIAGNOSTICS
165  CMLC--   Set diagnostic suffix for the current tracer  CMLC--   Set diagnostic suffix for the current tracer
166  CML      IF ( useDiagnostics ) THEN  CML      IF ( useDiagnostics ) THEN
167  CML        diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid )  CML        diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid )
168  CML      ENDIF  CML      ENDIF
169  CML#endif  CML#endif
170    
171          dBug = debugLevel.GE.debLevB
172         &     .AND. myIter.EQ.nIter0
173         &     .AND. ( tracerIdentity.EQ.GAD_HEFF .OR.
174         &             tracerIdentity.EQ.GAD_QICE2 )
175    c    &     .AND. tracerIdentity.EQ.GAD_HEFF
176    
177  C--   Set up work arrays with valid (i.e. not NaN) values  C--   Set up work arrays with valid (i.e. not NaN) values
178  C     These inital values do not alter the numerical results. They  C     These inital values do not alter the numerical results. They
179  C     just ensure that all memory references are to valid floating  C     just ensure that all memory references are to valid floating
180  C     point numbers. This prevents spurious hardware signals due to  C     point numbers. This prevents spurious hardware signals due to
181  C     uninitialised but inert locations.  C     uninitialised but inert locations.
182    #ifdef ALLOW_AUTODIFF_TAMC
183        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
184         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
         xA(i,j)      = 0. _d 0  
         yA(i,j)      = 0. _d 0  
         uTrans(i,j)  = 0. _d 0  
         vTrans(i,j)  = 0. _d 0  
         fVerT(i,j,1) = 0. _d 0  
         fVerT(i,j,2) = 0. _d 0  
 #ifdef ALLOW_AUTODIFF_TAMC  
185          localTij(i,j) = 0. _d 0          localTij(i,j) = 0. _d 0
 #endif  
186         ENDDO         ENDDO
187        ENDDO        ENDDO
188    #endif
189    
190  C--   Set tile-specific parameters for horizontal fluxes  C--   Set tile-specific parameters for horizontal fluxes
191        IF (useCubedSphereExchange) THEN        IF (useCubedSphereExchange) THEN
# Line 218  C--   Set tile-specific parameters for h Line 221  C--   Set tile-specific parameters for h
221        jMin = 1-OLy        jMin = 1-OLy
222        jMax = sNy+OLy        jMax = sNy+OLy
223    
 C--   Start of k loop for horizontal fluxes  
224        k = 1        k = 1
225  #ifdef ALLOW_AUTODIFF_TAMC  C--   Start of k loop for horizontal fluxes
226    #ifdef ALLOW_AUTODIFF_TAMC
227        kkey = (igadkey-1)*Nr + k        kkey = (igadkey-1)*Nr + k
228  CADJ STORE tracer(:,:,bi,bj) =  CADJ STORE iceFld =
229  CADJ &     comlev1_bibj_k_gad, key=kkey, byte=isbyte  CADJ &     comlev1_bibj_k_gad, key=kkey, byte=isbyte
230  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
231    
232  C     Content of CALC_COMMON_FACTORS, adapted for 2D fields  C     Content of CALC_COMMON_FACTORS, adapted for 2D fields
233  C--   Get temporary terms used by tendency routines  C--   Get temporary terms used by tendency routines
234    
235  C--   Calculate tracer cell face open areas  C--   Make local copy of sea-ice field and mask West & South
       DO j=jMin,jMax  
        DO i=iMin,iMax  
         xA(i,j) = _dyG(i,j,bi,bj) * maskW(i,j,k,bi,bj)  
         yA(i,j) = _dxG(i,j,bi,bj) * maskS(i,j,k,bi,bj)  
        ENDDO  
       ENDDO  
   
 C--   Calculate velocity field "volume transports" through  
 C--   tracer cell faces.  
       DO j=jMin,jMax  
        DO i=iMin,iMax  
         uFld(i,j) = uVel(i,j,bi,bj)  
         vFld(i,j) = vVel(i,j,bi,bj)  
         uTrans(i,j) = uVel(i,j,bi,bj)*xA(i,j)  
         vTrans(i,j) = vVel(i,j,bi,bj)*yA(i,j)  
        ENDDO  
       ENDDO  
 C     end of CALC_COMMON_FACTORS, adapted for 2D fields  
   
 C--   Make local copy of tracer array and mask West & South  
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          localTij(i,j)=tracer(i,j,bi,bj)          localTij(i,j)=iceFld(i,j)
239          maskLocW(i,j)=maskW(i,j,k,bi,bj)          maskLocW(i,j)=maskW(i,j,k,bi,bj)
240          maskLocS(i,j)=maskS(i,j,k,bi,bj)          maskLocS(i,j)=maskS(i,j,k,bi,bj)
241         ENDDO         ENDDO
242        ENDDO        ENDDO
243          
244    #ifdef ALLOW_AUTODIFF_TAMC
245    C-     Initialise Advective flux in X & Y
246           DO j=1-OLy,sNy+OLy
247            DO i=1-OLx,sNx+OLx
248             afx(i,j) = 0.
249             afy(i,j) = 0.
250            ENDDO
251           ENDDO
252    #endif
253    
254  #ifndef ALLOW_AUTODIFF_TAMC  #ifndef ALLOW_AUTODIFF_TAMC
255        IF (useCubedSphereExchange) THEN        IF (useCubedSphereExchange) THEN
256         withSigns = .FALSE.         withSigns = .FALSE.
257         CALL FILL_CS_CORNER_UV_RS(         CALL FILL_CS_CORNER_UV_RS(
258       &      withSigns, maskLocW,maskLocS, bi,bj, myThid )       &      withSigns, maskLocW,maskLocS, bi,bj, myThid )
259        ENDIF        ENDIF
260  #endif  #endif
# Line 280  C--   For cube need one pass for each of Line 273  C--   For cube need one pass for each of
273         interiorOnly = .FALSE.         interiorOnly = .FALSE.
274         overlapOnly  = .FALSE.         overlapOnly  = .FALSE.
275         IF (useCubedSphereExchange) THEN         IF (useCubedSphereExchange) THEN
276  C--   CubedSphere : pass 3 times, with partial update of local tracer field  C--   CubedSphere : pass 3 times, with partial update of local seaice field
277          IF (ipass.EQ.1) THEN          IF (ipass.EQ.1) THEN
278           overlapOnly  = MOD(nCFace,3).EQ.0           overlapOnly  = MOD(nCFace,3).EQ.0
279           interiorOnly = MOD(nCFace,3).NE.0           interiorOnly = MOD(nCFace,3).NE.0
# Line 299  C--   not CubedSphere Line 292  C--   not CubedSphere
292          calc_fluxes_X = MOD(ipass,2).EQ.1          calc_fluxes_X = MOD(ipass,2).EQ.1
293          calc_fluxes_Y = .NOT.calc_fluxes_X          calc_fluxes_Y = .NOT.calc_fluxes_X
294         ENDIF         ENDIF
295           IF (dBug.AND.bi.EQ.3 ) WRITE(6,*) 'ICE_adv:',tracerIdentity,
296         &   ipass,calc_fluxes_X,calc_fluxes_Y,overlapOnly,interiorOnly
297    
298  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
299  C--   X direction  C--   X direction
300  C-     Advective flux in X  
        DO j=1-Oly,sNy+Oly  
         DO i=1-Olx,sNx+Olx  
          af(i,j) = 0.  
         ENDDO  
        ENDDO  
 C  
301  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
302  CADJ STORE localTij(:,:)  =  CADJ STORE localTij(:,:)  =
303  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte
304  # ifndef DISABLE_MULTIDIM_ADVECTION  # ifndef DISABLE_MULTIDIM_ADVECTION
305  CADJ STORE af(:,:)  =  CADJ STORE af(:,:)  =
306  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte
307  # endif  # endif
308  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
309  C  C
310         IF (calc_fluxes_X) THEN         IF (calc_fluxes_X) THEN
311            
312  C-     Do not compute fluxes if  C-     Do not compute fluxes if
313  C       a) needed in overlap only  C       a) needed in overlap only
314  C   and b) the overlap of myTile are not cube-face Edges  C   and b) the overlap of myTile are not cube-face Edges
315          IF ( .NOT.overlapOnly .OR. N_edge .OR. S_edge ) THEN          IF ( .NOT.overlapOnly .OR. N_edge .OR. S_edge ) THEN
316    
317    C-     Advective flux in X
318             DO j=1-Oly,sNy+Oly
319              DO i=1-Olx,sNx+Olx
320               af(i,j) = 0.
321              ENDDO
322             ENDDO
323    
324  #ifndef ALLOW_AUTODIFF_TAMC  #ifndef ALLOW_AUTODIFF_TAMC
325  C-     Internal exchange for calculations in X  C-     Internal exchange for calculations in X
326  #ifdef MULTIDIM_OLD_VERSION  #ifdef MULTIDIM_OLD_VERSION
# Line 339  C-     Internal exchange for calculation Line 335  C-     Internal exchange for calculation
335    
336  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
337  # ifndef DISABLE_MULTIDIM_ADVECTION  # ifndef DISABLE_MULTIDIM_ADVECTION
338  CADJ STORE localTij(:,:)  =  CADJ STORE localTij(:,:)  =
339  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte
340  # endif  # endif
341  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
# Line 349  CADJ &     comlev1_bibj_k_gad_pass, key= Line 345  CADJ &     comlev1_bibj_k_gad_pass, key=
345            CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme,            CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme,
346       I         SEAICE_deltaTtherm,uTrans,uFld,localTij,       I         SEAICE_deltaTtherm,uTrans,uFld,localTij,
347       O         af, myThid )       O         af, myThid )
348              IF (dBug.AND. bi.EQ.3) WRITE(6,'(A,1P4E14.6)')
349         &      'ICE_adv: xFx=',af(13,11),localTij(12,11),uTrans(13,11),
350         &       af(13,11)/uTrans(13,11)
351           ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN           ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
352            CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, SEAICE_deltaTtherm,            CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, SEAICE_deltaTtherm,
353       I         uTrans, uFld, maskLocW, localTij,       I         uTrans, uFld, maskLocW, localTij,
# Line 362  CADJ &     comlev1_bibj_k_gad_pass, key= Line 361  CADJ &     comlev1_bibj_k_gad_pass, key=
361       I         uTrans, uFld, maskLocW, localTij,       I         uTrans, uFld, maskLocW, localTij,
362       O         af, myThid )       O         af, myThid )
363           ELSE           ELSE
364            STOP            STOP
365       & 'SEAICE_ADVECTION: adv. scheme incompatibale with multi-dim'       & 'SEAICE_ADVECTION: adv. scheme incompatibale with multi-dim'
366           ENDIF           ENDIF
367            
368  C--   Advective flux in X : done  C--   Advective flux in X : done
369          ENDIF          ENDIF
370            
371  #ifndef ALLOW_AUTODIFF_TAMC  #ifndef ALLOW_AUTODIFF_TAMC
372  C--   Internal exchange for next calculations in Y  C--   Internal exchange for next calculations in Y
373          IF ( overlapOnly .AND. ipass.EQ.1 ) THEN          IF ( overlapOnly .AND. ipass.EQ.1 ) THEN
# Line 376  C--   Internal exchange for next calcula Line 375  C--   Internal exchange for next calcula
375          ENDIF          ENDIF
376  #endif  #endif
377    
378  C-     Update the local tracer field where needed:  C-     Update the local seaice field where needed:
379    
380  C     update in overlap-Only  C     update in overlap-Only
381          IF ( overlapOnly ) THEN          IF ( overlapOnly ) THEN
382           iMinUpd = 1-Olx+1           iMinUpd = 1-OLx+1
383           iMaxUpd = sNx+Olx-1           iMaxUpd = sNx+OLx-1
384  C--   notes: these 2 lines below have no real effect (because recip_hFac=0  C--   notes: these 2 lines below have no real effect (because recip_hFac=0
385  C            in corner region) but safer to keep them.  C            in corner region) but safer to keep them.
386           IF ( W_edge ) iMinUpd = 1           IF ( W_edge ) iMinUpd = 1
387           IF ( E_edge ) iMaxUpd = sNx           IF ( E_edge ) iMaxUpd = sNx
388            
389           IF ( S_edge ) THEN           IF ( S_edge .AND. extensiveFld ) THEN
390            DO j=1-Oly,0            DO j=1-OLy,0
391             DO i=iMinUpd,iMaxUpd             DO i=iMinUpd,iMaxUpd
392              localTij(i,j)=localTij(i,j)-SEAICE_deltaTtherm*              localTij(i,j)=localTij(i,j)
393       &           maskC(i,j,k,bi,bj)       &         -SEAICE_deltaTtherm*maskC(i,j,k,bi,bj)
394       &           *recip_rA(i,j,bi,bj)       &           *recip_rA(i,j,bi,bj)
395       &           *( af(i+1,j)-af(i,j)       &           *(  af(i+1,j)-af(i,j)
396       &           )       &            )
397               ENDDO
398              ENDDO
399             ELSEIF ( S_edge ) THEN
400              DO j=1-OLy,0
401               DO i=iMinUpd,iMaxUpd
402                localTij(i,j)=localTij(i,j)
403         &         -SEAICE_deltaTtherm*maskC(i,j,k,bi,bj)
404         &           *recip_rA(i,j,bi,bj)*r_hFld(i,j)
405         &           *( (af(i+1,j)-af(i,j))
406         &             -(uTrans(i+1,j)-uTrans(i,j))*iceFld(i,j)
407         &            )
408             ENDDO             ENDDO
409            ENDDO            ENDDO
410           ENDIF           ENDIF
411           IF ( N_edge ) THEN           IF ( N_edge .AND. extensiveFld ) THEN
412            DO j=sNy+1,sNy+Oly            DO j=sNy+1,sNy+OLy
413             DO i=iMinUpd,iMaxUpd             DO i=iMinUpd,iMaxUpd
414              localTij(i,j)=localTij(i,j)-SEAICE_deltaTtherm*              localTij(i,j)=localTij(i,j)
415       &           maskC(i,j,k,bi,bj)       &         -SEAICE_deltaTtherm*maskC(i,j,k,bi,bj)
416       &           *recip_rA(i,j,bi,bj)       &           *recip_rA(i,j,bi,bj)
417       &           *( af(i+1,j)-af(i,j)       &           *(  af(i+1,j)-af(i,j)
418       &           )       &            )
419               ENDDO
420              ENDDO
421             ELSEIF ( N_edge ) THEN
422              DO j=sNy+1,sNy+OLy
423               DO i=iMinUpd,iMaxUpd
424                localTij(i,j)=localTij(i,j)
425         &         -SEAICE_deltaTtherm*maskC(i,j,k,bi,bj)
426         &           *recip_rA(i,j,bi,bj)*r_hFld(i,j)
427         &           *( (af(i+1,j)-af(i,j))
428         &             -(uTrans(i+1,j)-uTrans(i,j))*iceFld(i,j)
429         &            )
430             ENDDO             ENDDO
431            ENDDO            ENDDO
432           ENDIF           ENDIF
433            C--   keep advective flux (for diagnostics)
434             IF ( S_edge ) THEN
435              DO j=1-OLy,0
436               DO i=1-OLx+1,sNx+OLx
437                afx(i,j) = af(i,j)
438               ENDDO
439              ENDDO
440             ENDIF
441             IF ( N_edge ) THEN
442              DO j=sNy+1,sNy+OLy
443               DO i=1-OLx+1,sNx+OLx
444                afx(i,j) = af(i,j)
445               ENDDO
446              ENDDO
447             ENDIF
448    
449          ELSE          ELSE
450  C     do not only update the overlap  C     do not only update the overlap
451           jMinUpd = 1-Oly           jMinUpd = 1-OLy
452           jMaxUpd = sNy+Oly           jMaxUpd = sNy+OLy
453           IF ( interiorOnly .AND. S_edge ) jMinUpd = 1           IF ( interiorOnly .AND. S_edge ) jMinUpd = 1
454           IF ( interiorOnly .AND. N_edge ) jMaxUpd = sNy           IF ( interiorOnly .AND. N_edge ) jMaxUpd = sNy
455           DO j=jMinUpd,jMaxUpd           IF ( extensiveFld ) THEN
456            DO i=1-Olx+1,sNx+Olx-1            DO j=jMinUpd,jMaxUpd
457             localTij(i,j)=localTij(i,j)-SEAICE_deltaTtherm*             DO i=1-OLx+1,sNx+OLx-1
458       &           maskC(i,j,k,bi,bj)              localTij(i,j)=localTij(i,j)
459       &          *recip_rA(i,j,bi,bj)       &         -SEAICE_deltaTtherm*maskC(i,j,k,bi,bj)
460       &          *( af(i+1,j)-af(i,j)       &           *recip_rA(i,j,bi,bj)
461       &          )       &           *(  af(i+1,j)-af(i,j)
462         &            )
463               ENDDO
464            ENDDO            ENDDO
465           ENDDO           ELSE
466  C--   keep advective flux (for diagnostics)            DO j=jMinUpd,jMaxUpd
467           DO j=1-Oly,sNy+Oly             DO i=1-OLx+1,sNx+OLx-1
468            DO i=1-Olx,sNx+Olx              localTij(i,j)=localTij(i,j)
469             afx(i,j) = af(i,j)       &         -SEAICE_deltaTtherm*maskC(i,j,k,bi,bj)
470         &           *recip_rA(i,j,bi,bj)*r_hFld(i,j)
471         &           *( (af(i+1,j)-af(i,j))
472         &             -(uTrans(i+1,j)-uTrans(i,j))*iceFld(i,j)
473         &            )
474               ENDDO
475            ENDDO            ENDDO
476             ENDIF
477    C--   keep advective flux (for diagnostics)
478             DO j=jMinUpd,jMaxUpd
479               DO i=1-OLx+1,sNx+OLx
480                afx(i,j) = af(i,j)
481               ENDDO
482           ENDDO           ENDDO
483    
484  C     This is for later  C     This is for later
# Line 446  CML#endif /* ALLOW_OBCS */ Line 495  CML#endif /* ALLOW_OBCS */
495    
496  C-     end if/else update overlap-Only  C-     end if/else update overlap-Only
497          ENDIF          ENDIF
498            
499  C--   End of X direction  C--   End of X direction
500         ENDIF         ENDIF
501            
502  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
503  C--   Y direction  C--   Y direction
504  cph-test  
 C-     Advective flux in Y  
        DO j=1-Oly,sNy+Oly  
         DO i=1-Olx,sNx+Olx  
          af(i,j) = 0.  
         ENDDO  
        ENDDO  
 C  
505  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
506  # ifndef DISABLE_MULTIDIM_ADVECTION  # ifndef DISABLE_MULTIDIM_ADVECTION
507  CADJ STORE localTij(:,:)  =  CADJ STORE localTij(:,:)  =
508  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte
509  CADJ STORE af(:,:)  =  CADJ STORE af(:,:)  =
510  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte
511  # endif  # endif
512  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
513  C  
514         IF (calc_fluxes_Y) THEN         IF (calc_fluxes_Y) THEN
515    
516  C-     Do not compute fluxes if  C-     Do not compute fluxes if
# Line 476  C       a) needed in overlap only Line 518  C       a) needed in overlap only
518  C   and b) the overlap of myTile are not cube-face edges  C   and b) the overlap of myTile are not cube-face edges
519          IF ( .NOT.overlapOnly .OR. E_edge .OR. W_edge ) THEN          IF ( .NOT.overlapOnly .OR. E_edge .OR. W_edge ) THEN
520    
521    C-     Advective flux in Y
522             DO j=1-OLy,sNy+OLy
523              DO i=1-OLx,sNx+OLx
524               af(i,j) = 0.
525              ENDDO
526             ENDDO
527    
528  #ifndef ALLOW_AUTODIFF_TAMC  #ifndef ALLOW_AUTODIFF_TAMC
529  C-     Internal exchange for calculations in Y  C-     Internal exchange for calculations in Y
530  #ifdef MULTIDIM_OLD_VERSION  #ifdef MULTIDIM_OLD_VERSION
531           IF ( useCubedSphereExchange ) THEN           IF ( useCubedSphereExchange ) THEN
532  #else  #else
533           IF ( useCubedSphereExchange .AND.           IF ( useCubedSphereExchange .AND.
534       &      ( overlapOnly .OR. ipass.EQ.1 ) ) THEN       &      ( overlapOnly .OR. ipass.EQ.1 ) ) THEN
535  #endif  #endif
# Line 488  C-     Internal exchange for calculation Line 537  C-     Internal exchange for calculation
537           ENDIF           ENDIF
538  #endif  #endif
539    
540  C-     Advective flux in Y  #ifdef ALLOW_AUTODIFF_TAMC
          DO j=1-Oly,sNy+Oly  
           DO i=1-Olx,sNx+Olx  
            af(i,j) = 0.  
           ENDDO  
          ENDDO  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
541  #ifndef DISABLE_MULTIDIM_ADVECTION  #ifndef DISABLE_MULTIDIM_ADVECTION
542  CADJ STORE localTij(:,:)  =  CADJ STORE localTij(:,:)  =
543  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte  CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte
544  #endif  #endif
545  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
# Line 507  CADJ &     comlev1_bibj_k_gad_pass, key= Line 549  CADJ &     comlev1_bibj_k_gad_pass, key=
549            CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme,            CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme,
550       I         SEAICE_deltaTtherm,vTrans,vFld,localTij,       I         SEAICE_deltaTtherm,vTrans,vFld,localTij,
551       O         af, myThid )       O         af, myThid )
552              IF (dBug.AND. bi.EQ.3) WRITE(6,'(A,1P4E14.6)')
553         &      'ICE_adv: yFx=',af(12,12),localTij(12,11),vTrans(12,12),
554         &       af(12,12)/vTrans(12,12)
555           ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN           ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
556            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, SEAICE_deltaTtherm,            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, SEAICE_deltaTtherm,
557       I         vTrans, vFld, maskLocS, localTij,       I         vTrans, vFld, maskLocS, localTij,
# Line 520  CADJ &     comlev1_bibj_k_gad_pass, key= Line 565  CADJ &     comlev1_bibj_k_gad_pass, key=
565       I         vTrans, vFld, maskLocS, localTij,       I         vTrans, vFld, maskLocS, localTij,
566       O         af, myThid )       O         af, myThid )
567           ELSE           ELSE
568            STOP            STOP
569       &  'SEAICE_ADVECTION: adv. scheme incompatibale with mutli-dim'       &  'SEAICE_ADVECTION: adv. scheme incompatibale with mutli-dim'
570           ENDIF           ENDIF
571    
# Line 534  C-     Internal exchange for next calcul Line 579  C-     Internal exchange for next calcul
579          ENDIF          ENDIF
580  #endif  #endif
581    
582  C-     Update the local tracer field where needed:  C-     Update the local seaice field where needed:
583    
584  C      update in overlap-Only  C      update in overlap-Only
585          IF ( overlapOnly ) THEN          IF ( overlapOnly ) THEN
586           jMinUpd = 1-Oly+1           jMinUpd = 1-OLy+1
587           jMaxUpd = sNy+Oly-1           jMaxUpd = sNy+OLy-1
588  C- notes: these 2 lines below have no real effect (because recip_hFac=0  C- notes: these 2 lines below have no real effect (because recip_hFac=0
589  C         in corner region) but safer to keep them.  C         in corner region) but safer to keep them.
590           IF ( S_edge ) jMinUpd = 1           IF ( S_edge ) jMinUpd = 1
591           IF ( N_edge ) jMaxUpd = sNy           IF ( N_edge ) jMaxUpd = sNy
592            
593           IF ( W_edge ) THEN           IF ( W_edge .AND. extensiveFld ) THEN
594            DO j=jMinUpd,jMaxUpd            DO j=jMinUpd,jMaxUpd
595             DO i=1-Olx,0             DO i=1-OLx,0
596              localTij(i,j)=localTij(i,j)-SEAICE_deltaTtherm*              localTij(i,j)=localTij(i,j)
597       &           maskC(i,j,k,bi,bj)       &         -SEAICE_deltaTtherm*maskC(i,j,k,bi,bj)
598       &           *recip_rA(i,j,bi,bj)       &           *recip_rA(i,j,bi,bj)
599       &           *( af(i,j+1)-af(i,j)       &           *(  af(i,j+1)-af(i,j)
600       &           )       &            )
601               ENDDO
602              ENDDO
603             ELSEIF ( W_edge ) THEN
604              DO j=jMinUpd,jMaxUpd
605               DO i=1-OLx,0
606                localTij(i,j)=localTij(i,j)
607         &         -SEAICE_deltaTtherm*maskC(i,j,k,bi,bj)
608         &           *recip_rA(i,j,bi,bj)*r_hFld(i,j)
609         &           *( (af(i,j+1)-af(i,j))
610         &             -(vTrans(i,j+1)-vTrans(i,j))*iceFld(i,j)
611         &            )
612             ENDDO             ENDDO
613            ENDDO            ENDDO
614           ENDIF           ENDIF
615           IF ( E_edge ) THEN           IF ( E_edge .AND. extensiveFld ) THEN
616            DO j=jMinUpd,jMaxUpd            DO j=jMinUpd,jMaxUpd
617             DO i=sNx+1,sNx+Olx             DO i=sNx+1,sNx+OLx
618              localTij(i,j)=localTij(i,j)-SEAICE_deltaTtherm*              localTij(i,j)=localTij(i,j)
619       &           maskC(i,j,k,bi,bj)       &         -SEAICE_deltaTtherm*maskC(i,j,k,bi,bj)
620       &           *recip_rA(i,j,bi,bj)       &           *recip_rA(i,j,bi,bj)
621       &           *( af(i,j+1)-af(i,j)       &           *(  af(i,j+1)-af(i,j)
622       &           )       &            )
623               ENDDO
624              ENDDO
625             ELSEIF ( E_edge ) THEN
626              DO j=jMinUpd,jMaxUpd
627               DO i=sNx+1,sNx+OLx
628                localTij(i,j)=localTij(i,j)
629         &         -SEAICE_deltaTtherm*maskC(i,j,k,bi,bj)
630         &           *recip_rA(i,j,bi,bj)*r_hFld(i,j)
631         &           *( (af(i,j+1)-af(i,j))
632         &             -(vTrans(i,j+1)-vTrans(i,j))*iceFld(i,j)
633         &            )
634             ENDDO             ENDDO
635            ENDDO            ENDDO
636           ENDIF           ENDIF
637            C--   keep advective flux (for diagnostics)
638             IF ( W_edge ) THEN
639              DO j=1-OLy+1,sNy+OLy
640               DO i=1-OLx,0
641                afy(i,j) = af(i,j)
642               ENDDO
643              ENDDO
644             ENDIF
645             IF ( E_edge ) THEN
646              DO j=1-OLy+1,sNy+OLy
647               DO i=sNx+1,sNx+OLx
648                afy(i,j) = af(i,j)
649               ENDDO
650              ENDDO
651             ENDIF
652    
653          ELSE          ELSE
654  C     do not only update the overlap  C     do not only update the overlap
655           iMinUpd = 1-Olx           iMinUpd = 1-OLx
656           iMaxUpd = sNx+Olx           iMaxUpd = sNx+OLx
657           IF ( interiorOnly .AND. W_edge ) iMinUpd = 1           IF ( interiorOnly .AND. W_edge ) iMinUpd = 1
658           IF ( interiorOnly .AND. E_edge ) iMaxUpd = sNx           IF ( interiorOnly .AND. E_edge ) iMaxUpd = sNx
659           DO j=1-Oly+1,sNy+Oly-1           IF ( extensiveFld ) THEN
660            DO i=iMinUpd,iMaxUpd            DO j=1-OLy+1,sNy+OLy-1
661             localTij(i,j)=localTij(i,j)-SEAICE_deltaTtherm*             DO i=iMinUpd,iMaxUpd
662       &           maskC(i,j,k,bi,bj)              localTij(i,j)=localTij(i,j)
663       &          *recip_rA(i,j,bi,bj)       &         -SEAICE_deltaTtherm*maskC(i,j,k,bi,bj)
664       &          *( af(i,j+1)-af(i,j)       &           *recip_rA(i,j,bi,bj)
665       &          )       &           *(  af(i,j+1)-af(i,j)
666         &            )
667               ENDDO
668            ENDDO            ENDDO
669           ENDDO           ELSE
670  C--   keep advective flux (for diagnostics)            DO j=1-OLy+1,sNy+OLy-1
671           DO j=1-Oly,sNy+Oly             DO i=iMinUpd,iMaxUpd
672            DO i=1-Olx,sNx+Olx              localTij(i,j)=localTij(i,j)
673             afy(i,j) = af(i,j)       &         -SEAICE_deltaTtherm*maskC(i,j,k,bi,bj)
674         &           *recip_rA(i,j,bi,bj)*r_hFld(i,j)
675         &           *( (af(i,j+1)-af(i,j))
676         &             -(vTrans(i,j+1)-vTrans(i,j))*iceFld(i,j)
677         &            )
678               ENDDO
679            ENDDO            ENDDO
680             ENDIF
681    C--   keep advective flux (for diagnostics)
682             DO j=1-OLy+1,sNy+OLy
683               DO i=iMinUpd,iMaxUpd
684                afy(i,j) = af(i,j)
685               ENDDO
686           ENDDO           ENDDO
687    
688  C--   Save this for later  C--   Save this for later
# Line 611  C--   End of Y direction Line 706  C--   End of Y direction
706  C--   End of ipass loop  C--   End of ipass loop
707        ENDDO        ENDDO
708    
709  C-    explicit advection is done ; store tendency in gTracer:  C-    explicit advection is done ; store tendency in gFld:
710        DO j=1-Oly,sNy+Oly        DO j=1-OLy,sNy+OLy
711         DO i=1-Olx,sNx+Olx         DO i=1-OLx,sNx+OLx
712          gTracer(i,j,bi,bj)=          gFld(i,j)=
713       &       (localTij(i,j)-tracer(i,j,bi,bj))/SEAICE_deltaTtherm       &       (localTij(i,j)-iceFld(i,j))/SEAICE_deltaTtherm
714            IF ( dBug
715         &     .AND. i.EQ.12 .AND. j.EQ.11 .AND. bi.EQ.3 ) THEN
716             tmpFac= SEAICE_deltaTtherm*recip_rA(i,j,bi,bj)
717             WRITE(6,'(A,1P4E14.6)') 'ICE_adv:',
718         &    afx(i,j)*tmpFac,afx(i+1,j)*tmpFac,
719         &    afy(i,j)*tmpFac,afy(i,j+1)*tmpFac
720            ENDIF
721         ENDDO         ENDDO
722        ENDDO        ENDDO
723          
724  CML#ifdef ALLOW_DIAGNOSTICS  CML#ifdef ALLOW_DIAGNOSTICS
725  CML        IF ( useDiagnostics ) THEN  CML        IF ( useDiagnostics ) THEN
726  CML         diagName = 'ADVx'//diagSufx  CML         diagName = 'ADVx'//diagSufx
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.22