/[MITgcm]/MITgcm/pkg/thsice/thsice_advection.F
ViewVC logotype

Diff of /MITgcm/pkg/thsice/thsice_advection.F

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

revision 1.14 by heimbach, Wed Jun 1 22:38:25 2011 UTC revision 1.15 by jmc, Thu Jun 2 18:34:34 2011 UTC
# Line 14  C !INTERFACE: ========================== Line 14  C !INTERFACE: ==========================
14        SUBROUTINE THSICE_ADVECTION(        SUBROUTINE THSICE_ADVECTION(
15       I     tracerIdentity,       I     tracerIdentity,
16       I     advectionScheme,       I     advectionScheme,
17       I     useGridVolume,       I     useGridArea,
18       I     uTrans, vTrans, maskOc, deltaTadvect, iceEps,       I     uTrans, vTrans, maskOc, deltaTadvect, iceEps,
19       U     iceVol, iceFld,       U     iceVol, iceFld,
20       O     afx, afy,       O     afx, afy,
# Line 27  C and can only be used for the non-linea Line 27  C and can only be used for the non-linea
27  C direct-space-time method and flux-limiters.  C direct-space-time method and flux-limiters.
28  C  C
29  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.
30  C for Area, effective thickness or other "extensive" sea-ice field,  C for Area, effective thickness or other sea-ice field,
31  C  the contribution iceFld*div(u) (that is present in gad_advection)  C  the contribution iceFld*div(u) (that is present in gad_advection)
32  C  is not included here.  C  is not included here.
33  C  C
# Line 48  C !USES: =============================== Line 48  C !USES: ===============================
48  #include "EEPARAMS.h"  #include "EEPARAMS.h"
49  #include "PARAMS.h"  #include "PARAMS.h"
50  #include "GRID.h"  #include "GRID.h"
51  #include "THSICE_SIZE.h"  #if ( defined ALLOW_DBUG_THSICE || defined ALLOW_AUTODIFF_TAMC )
52    # include "THSICE_SIZE.h"
53    #endif
54  #ifdef ALLOW_GENERIC_ADVDIFF  #ifdef ALLOW_GENERIC_ADVDIFF
55  # include "GAD.h"  # include "GAD.h"
56  #endif  #endif
# Line 57  C !USES: =============================== Line 59  C !USES: ===============================
59  #include "W2_EXCH2_TOPOLOGY.h"  #include "W2_EXCH2_TOPOLOGY.h"
60  #endif /* ALLOW_EXCH2 */  #endif /* ALLOW_EXCH2 */
61  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
 # include "THSICE_PARAMS.h"  
62  # include "tamc.h"  # include "tamc.h"
63  # include "tamc_keys.h"  # include "tamc_keys.h"
64  #endif  #endif
# Line 65  C !USES: =============================== Line 66  C !USES: ===============================
66  C !INPUT PARAMETERS: ===================================================  C !INPUT PARAMETERS: ===================================================
67  C  tracerIdentity  :: tracer identifier (required only for OBCS)  C  tracerIdentity  :: tracer identifier (required only for OBCS)
68  C  advectionScheme :: advection scheme to use (Horizontal plane)  C  advectionScheme :: advection scheme to use (Horizontal plane)
69  C  useGridVolume   :: use grid-cell Area & Volume (instead of "iceVol" field)  C  useGridArea     :: use grid-cell Area (instead of "iceVol" field)
70  C  uTrans,vTrans   :: volume transports at U,V points  C  uTrans,vTrans   :: volume transports at U,V points
71  C  maskOc          :: oceanic mask  C  maskOc          :: oceanic mask
72  C  iceVol          :: sea-ice volume  C  iceVol          :: sea-ice volume
# Line 78  C  myIter          :: current iteration Line 79  C  myIter          :: current iteration
79  C  myThid          :: my thread Id. number  C  myThid          :: my thread Id. number
80        INTEGER tracerIdentity        INTEGER tracerIdentity
81        INTEGER advectionScheme        INTEGER advectionScheme
82        LOGICAL useGridVolume        LOGICAL useGridArea
83        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
84        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
85        _RS maskOc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskOc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 100  C  afy             :: horizontal advecti Line 101  C  afy             :: horizontal advecti
101    
102  #ifdef ALLOW_GENERIC_ADVDIFF  #ifdef ALLOW_GENERIC_ADVDIFF
103  C !LOCAL VARIABLES: ====================================================  C !LOCAL VARIABLES: ====================================================
104    C  maskLocC      :: 2-D array mask at grid-cell center
105  C  maskLocW      :: 2-D array for mask at West points  C  maskLocW      :: 2-D array for mask at West points
106  C  maskLocS      :: 2-D array for mask at South points  C  maskLocS      :: 2-D array for mask at South points
107  C  iMin,iMax,    :: loop range for called routines  C  iMin,iMax,    :: loop range for called routines
# Line 120  C  myTile        :: variables used to de Line 122  C  myTile        :: variables used to de
122  C  nCFace        :: owns a tile for cube grid runs using  C  nCFace        :: owns a tile for cube grid runs using
123  C                :: multi-dim advection.  C                :: multi-dim advection.
124  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
125          _RS maskLocC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
126        _RS maskLocW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskLocW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
127        _RS maskLocS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskLocS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
128        INTEGER iMin,iMax,jMin,jMax        INTEGER iMin,iMax,jMin,jMax
# Line 128  C [N,S,E,W]_edge :: true if N,S,E,W edge Line 131  C [N,S,E,W]_edge :: true if N,S,E,W edge
131        _RL uCfl    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uCfl    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
132        _RL vCfl    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vCfl    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
133        _RL af      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL af      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
134          _RL tmpVol
135        LOGICAL calc_fluxes_X, calc_fluxes_Y, withSigns        LOGICAL calc_fluxes_X, calc_fluxes_Y, withSigns
136        LOGICAL interiorOnly, overlapOnly        LOGICAL interiorOnly, overlapOnly
137        INTEGER nipass,ipass        INTEGER nipass,ipass
# Line 136  C [N,S,E,W]_edge :: true if N,S,E,W edge Line 140  C [N,S,E,W]_edge :: true if N,S,E,W edge
140  #ifdef ALLOW_EXCH2  #ifdef ALLOW_EXCH2
141        INTEGER myTile        INTEGER myTile
142  #endif  #endif
143    #ifdef ALLOW_DBUG_THSICE
144        LOGICAL dBugFlag        LOGICAL dBugFlag
145        INTEGER idb,jdb,biDb        INTEGER idb,jdb,biDb
146        _RL tmpFac        _RL tmpFac
147        _RL tmpVol  #endif /* ALLOW_DBUG_THSICE */
148  CEOP  CEOP
149    
150  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
151    
152          k = 1
153  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
154        act1 = bi - myBxLo(myThid)        act1 = bi - myBxLo(myThid)
155        max1 = myBxHi(myThid) - myBxLo(myThid) + 1        max1 = myBxHi(myThid) - myBxLo(myThid) + 1
# Line 157  C---+----1----+----2----+----3----+----4 Line 163  C---+----1----+----2----+----3----+----4
163       &                     + act4*max1*max2*max3       &                     + act4*max1*max2*max3
164  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
165    
166        k = 1  #ifdef ALLOW_DBUG_THSICE
167        dBugFlag = debugLevel.GE.debLevB        dBugFlag = debugLevel.GE.debLevB
168       &     .AND. myIter.EQ.nIter0       &     .AND. myIter.EQ.nIter0
169       &     .AND. ( tracerIdentity.EQ.GAD_SI_HICE .OR.       &     .AND. ( tracerIdentity.EQ.GAD_SI_HICE .OR.
# Line 166  c    &     .AND. tracerIdentity.EQ.GAD_S Line 172  c    &     .AND. tracerIdentity.EQ.GAD_S
172        idb  = MIN(13,sNx)        idb  = MIN(13,sNx)
173        jdb  = MOD(17,sNy)        jdb  = MOD(17,sNy)
174        biDb = nSx/2        biDb = nSx/2
175    #endif /* ALLOW_DBUG_THSICE */
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
# Line 206  C--   Set tile-specific parameters for h Line 213  C--   Set tile-specific parameters for h
213    
214  C--   Start horizontal fluxes  C--   Start horizontal fluxes
215    
216  C--   set mask West & South  C--   set mask West & South (and local Centered mask)
217        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
218         maskLocW(1-Olx,j) = 0.         maskLocW(1-Olx,j) = 0.
219         DO i=2-OLx,sNx+OLx         DO i=2-OLx,sNx+OLx
220          maskLocW(i,j) = MIN( maskOc(i-1,j), maskOc(i,j) )          maskLocW(i,j) = MIN( maskOc(i-1,j), maskOc(i,j) )
221    #ifdef ALLOW_OBCS
222         &                * maskInW(i,j,bi,bj)
223    #endif /* ALLOW_OBCS */
224         ENDDO         ENDDO
225        ENDDO        ENDDO
226        DO i=1-OLx,sNx+OLx        DO i=1-OLx,sNx+OLx
# Line 219  C--   set mask West & South Line 229  C--   set mask West & South
229        DO j=2-OLy,sNy+OLy        DO j=2-OLy,sNy+OLy
230         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
231          maskLocS(i,j) = MIN( maskOc(i,j-1), maskOc(i,j) )          maskLocS(i,j) = MIN( maskOc(i,j-1), maskOc(i,j) )
232    #ifdef ALLOW_OBCS
233         &                * maskInS(i,j,bi,bj)
234    #endif /* ALLOW_OBCS */
235           ENDDO
236          ENDDO
237    C     maskLocC is just a local copy of Ocean mask (maskOc) except if using OBCS:
238    C     use "maksInC" to prevent updating tracer field in OB regions
239          DO j=1-OLy,sNy+OLy
240           DO i=1-OLx,sNx+OLx
241            maskLocC(i,j) = maskOc(i,j)
242    #ifdef ALLOW_OBCS
243         &                * maskInC(i,j,bi,bj)
244    #endif /* ALLOW_OBCS */
245         ENDDO         ENDDO
246        ENDDO        ENDDO
247    
# Line 271  C--   not CubedSphere Line 294  C--   not CubedSphere
294          calc_fluxes_X = MOD(ipass,2).EQ.1          calc_fluxes_X = MOD(ipass,2).EQ.1
295          calc_fluxes_Y = .NOT.calc_fluxes_X          calc_fluxes_Y = .NOT.calc_fluxes_X
296         ENDIF         ENDIF
297    #ifdef ALLOW_DBUG_THSICE
298         IF (dBugFlag.AND.bi.EQ.biDb ) WRITE(6,'(A,3I4,2I5,4L5)')         IF (dBugFlag.AND.bi.EQ.biDb ) WRITE(6,'(A,3I4,2I5,4L5)')
299       &   'ICE_adv:', tracerIdentity, ipass, bi, idb, jdb,       &   'ICE_adv:', tracerIdentity, ipass, bi, idb, jdb,
300       &   calc_fluxes_X, calc_fluxes_Y, overlapOnly, interiorOnly       &   calc_fluxes_X, calc_fluxes_Y, overlapOnly, interiorOnly
301    #endif /* ALLOW_DBUG_THSICE */
302    
303  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
304  C--   X direction  C--   X direction
# Line 293  CADJ STORE af(:,:)     = comlev1_thsice_ Line 318  CADJ STORE af(:,:)     = comlev1_thsice_
318  C-     Advective flux in X  C-     Advective flux in X
319           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
320            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
321             af(i,j) = 0.              af(i,j) = 0.
322            ENDDO            ENDDO
323           ENDDO           ENDDO
324    
# Line 302  C-     Internal exchange for calculation Line 327  C-     Internal exchange for calculation
327           IF ( overlapOnly ) THEN           IF ( overlapOnly ) THEN
328            CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,            CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
329       &                               iceFld, bi,bj, myThid )       &                               iceFld, bi,bj, myThid )
330            IF ( .NOT.useGridVolume )            IF ( .NOT.useGridArea )
331       &    CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,       &    CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
332       &                               iceVol, bi,bj, myThid )       &                               iceVol, bi,bj, myThid )
333           ENDIF           ENDIF
334  cph-exch2#endif  cph-exch2#endif
335    
336  C-     Compute CFL number  C-     Compute CFL number
337           IF ( useGridVolume ) THEN           IF ( useGridArea ) THEN
338            DO j=1-Oly,sNy+Oly            DO j=1-Oly,sNy+Oly
339             DO i=2-Olx,sNx+Olx             DO i=2-Olx,sNx+Olx
340              uCfl(i,j) = deltaTadvect*(              uCfl(i,j) = deltaTadvect*(
# Line 334  C-     Compute CFL number Line 359  C-     Compute CFL number
359            CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme, .FALSE.,            CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme, .FALSE.,
360       I             deltaTadvect, uTrans, uCfl, iceFld,       I             deltaTadvect, uTrans, uCfl, iceFld,
361       O             af, myThid )       O             af, myThid )
362    #ifdef ALLOW_DBUG_THSICE
363            IF (dBugFlag.AND.bi.EQ.biDb ) WRITE(6,'(A,1P4E14.6)')            IF (dBugFlag.AND.bi.EQ.biDb ) WRITE(6,'(A,1P4E14.6)')
364       &      'ICE_adv: xFx=', af(idb,jdb), iceFld(idb,jdb),       &      'ICE_adv: xFx=', af(idb,jdb), iceFld(idb,jdb),
365       &       uTrans(idb,jdb), af(idb,jdb)/uTrans(idb,jdb)       &       uTrans(idb,jdb), af(idb,jdb)/uTrans(idb,jdb)
366    #endif /* ALLOW_DBUG_THSICE */
367           ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN           ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
368            CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, .FALSE., deltaTadvect,            CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, .FALSE., deltaTadvect,
369       I             uTrans, uCfl, maskLocW, iceFld,       I             uTrans, uCfl, maskLocW, iceFld,
# Line 359  C--   Internal exchange for next calcula Line 386  C--   Internal exchange for next calcula
386           IF ( overlapOnly .AND. ipass.EQ.1 ) THEN           IF ( overlapOnly .AND. ipass.EQ.1 ) THEN
387            CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,            CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
388       &                               iceFld, bi,bj, myThid )       &                               iceFld, bi,bj, myThid )
389            IF ( .NOT.useGridVolume )            IF ( .NOT.useGridArea )
390       &    CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,       &    CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
391       &                               iceVol, bi,bj, myThid )       &                               iceVol, bi,bj, myThid )
392           ENDIF           ENDIF
# Line 383  C            in corner region) but safer Line 410  C            in corner region) but safer
410  CADJ STORE iceFld(:,:) = comlev1_thsice_4, key=ikey_4, byte=isbyte  CADJ STORE iceFld(:,:) = comlev1_thsice_4, key=ikey_4, byte=isbyte
411  CADJ STORE iceVol(:,:) = comlev1_thsice_4, key=ikey_4, byte=isbyte  CADJ STORE iceVol(:,:) = comlev1_thsice_4, key=ikey_4, byte=isbyte
412  #endif  #endif
413           IF ( S_edge .AND. useGridVolume ) THEN           IF ( S_edge ) THEN
414            DO j=1-OLy,0            IF ( useGridArea ) THEN
415             DO i=iMinUpd,iMaxUpd             DO j=1-OLy,0
416              iceFld(i,j) = iceFld(i,j)              DO i=iMinUpd,iMaxUpd
417       &                   -deltaTadvect*maskOc(i,j)               iceFld(i,j) = iceFld(i,j)
418       &                                *recip_rA(i,j,bi,bj)       &                    -deltaTadvect*maskLocC(i,j)
419       &                                *( af(i+1,j)-af(i,j) )       &                                 *recip_rA(i,j,bi,bj)
420             ENDDO       &                                 *( af(i+1,j)-af(i,j) )
421            ENDDO              ENDDO
422           ELSEIF ( S_edge ) THEN             ENDDO
423            DO j=1-OLy,0            ELSE
424             DO i=iMinUpd,iMaxUpd             DO j=1-OLy,0
425              tmpVol = iceVol(i,j)              DO i=iMinUpd,iMaxUpd
426              iceVol(i,j) = iceVol(i,j)               tmpVol = iceVol(i,j)
427       &                   -deltaTadvect*maskOc(i,j)               iceVol(i,j) = iceVol(i,j)
428       &                          *( uTrans(i+1,j)-uTrans(i,j) )       &                    -deltaTadvect*maskLocC(i,j)
429              IF ( iceVol(i,j).GT.iceEps )       &                           *( uTrans(i+1,j)-uTrans(i,j) )
430       &      iceFld(i,j) = ( iceFld(i,j)*tmpVol               IF ( iceVol(i,j).GT.iceEps )
431       &                     -deltaTadvect*maskOc(i,j)       &       iceFld(i,j) = ( iceFld(i,j)*tmpVol
432       &                                  *( af(i+1,j)-af(i,j) )       &                      -deltaTadvect*maskLocC(i,j)
      &                    )/iceVol(i,j)  
            ENDDO  
           ENDDO  
          ENDIF  
          IF ( N_edge .AND. useGridVolume ) THEN  
           DO j=sNy+1,sNy+OLy  
            DO i=iMinUpd,iMaxUpd  
             iceFld(i,j) = iceFld(i,j)  
      &                   -deltaTadvect*maskOc(i,j)  
      &                                *recip_rA(i,j,bi,bj)  
      &                                *( af(i+1,j)-af(i,j) )  
            ENDDO  
           ENDDO  
          ELSEIF ( N_edge ) THEN  
           DO j=sNy+1,sNy+OLy  
            DO i=iMinUpd,iMaxUpd  
             tmpVol = iceVol(i,j)  
             iceVol(i,j) = iceVol(i,j)  
      &                   -deltaTadvect*maskOc(i,j)  
      &                          *( uTrans(i+1,j)-uTrans(i,j) )  
             IF ( iceVol(i,j).GT.iceEps )  
      &      iceFld(i,j) = ( iceFld(i,j)*tmpVol  
      &                     -deltaTadvect*maskOc(i,j)  
433       &                                  *( af(i+1,j)-af(i,j) )       &                                  *( af(i+1,j)-af(i,j) )
434       &                    )/iceVol(i,j)       &                     )/iceVol(i,j)
435                ENDDO
436             ENDDO             ENDDO
437            ENDDO            ENDIF
438           ENDIF  C-    keep advective flux (for diagnostics)
 C--   keep advective flux (for diagnostics)  
          IF ( S_edge ) THEN  
439            DO j=1-OLy,0            DO j=1-OLy,0
440             DO i=1-OLx+1,sNx+OLx             DO i=1-OLx+1,sNx+OLx
441              afx(i,j) = af(i,j)               afx(i,j) = af(i,j)
442             ENDDO             ENDDO
443            ENDDO            ENDDO
444    C-    end if South Edge
445           ENDIF           ENDIF
446           IF ( N_edge ) THEN           IF ( N_edge ) THEN
447              IF ( useGridArea ) THEN
448               DO j=sNy+1,sNy+OLy
449                DO i=iMinUpd,iMaxUpd
450                 iceFld(i,j) = iceFld(i,j)
451         &                    -deltaTadvect*maskLocC(i,j)
452         &                                 *recip_rA(i,j,bi,bj)
453         &                                 *( af(i+1,j)-af(i,j) )
454                ENDDO
455               ENDDO
456              ELSE
457               DO j=sNy+1,sNy+OLy
458                DO i=iMinUpd,iMaxUpd
459                 tmpVol = iceVol(i,j)
460                 iceVol(i,j) = iceVol(i,j)
461         &                    -deltaTadvect*maskLocC(i,j)
462         &                           *( uTrans(i+1,j)-uTrans(i,j) )
463                 IF ( iceVol(i,j).GT.iceEps )
464         &       iceFld(i,j) = ( iceFld(i,j)*tmpVol
465         &                      -deltaTadvect*maskLocC(i,j)
466         &                                   *( af(i+1,j)-af(i,j) )
467         &                     )/iceVol(i,j)
468                ENDDO
469               ENDDO
470              ENDIF
471    C-    keep advective flux (for diagnostics)
472            DO j=sNy+1,sNy+OLy            DO j=sNy+1,sNy+OLy
473             DO i=1-OLx+1,sNx+OLx             DO i=1-OLx+1,sNx+OLx
474              afx(i,j) = af(i,j)               afx(i,j) = af(i,j)
475             ENDDO             ENDDO
476            ENDDO            ENDDO
477    C-    end if North Edge
478           ENDIF           ENDIF
479    
480          ELSE          ELSE
481  C     do not only update the overlap  C     do not only update the overlap
482           jMinUpd = 1-OLy            jMinUpd = 1-OLy
483           jMaxUpd = sNy+OLy            jMaxUpd = sNy+OLy
484           IF ( interiorOnly .AND. S_edge ) jMinUpd = 1            IF ( interiorOnly .AND. S_edge ) jMinUpd = 1
485           IF ( interiorOnly .AND. N_edge ) jMaxUpd = sNy            IF ( interiorOnly .AND. N_edge ) jMaxUpd = sNy
486           IF ( useGridVolume ) THEN            IF ( useGridArea ) THEN
487            DO j=jMinUpd,jMaxUpd             DO j=jMinUpd,jMaxUpd
488             DO i=1-OLx+1,sNx+OLx-1              DO i=1-OLx+1,sNx+OLx-1
489              iceFld(i,j) = iceFld(i,j)               iceFld(i,j) = iceFld(i,j)
490       &                   -deltaTadvect*maskOc(i,j)       &                    -deltaTadvect*maskLocC(i,j)
491       &                                *recip_rA(i,j,bi,bj)       &                                 *recip_rA(i,j,bi,bj)
492       &                                *( af(i+1,j)-af(i,j) )       &                                 *( af(i+1,j)-af(i,j) )
493                ENDDO
494               ENDDO
495              ELSE
496               DO j=jMinUpd,jMaxUpd
497                DO i=1-OLx+1,sNx+OLx-1
498                 tmpVol = iceVol(i,j)
499                 iceVol(i,j) = iceVol(i,j)
500         &                    -deltaTadvect*maskLocC(i,j)
501         &                           *( uTrans(i+1,j)-uTrans(i,j) )
502                 IF ( iceVol(i,j).GT.iceEps )
503         &       iceFld(i,j) = ( iceFld(i,j)*tmpVol
504         &                      -deltaTadvect*maskLocC(i,j)
505         &                                   *( af(i+1,j)-af(i,j) )
506         &                     )/iceVol(i,j)
507                ENDDO
508             ENDDO             ENDDO
509            ENDDO            ENDIF
510           ELSE  C-    keep advective flux (for diagnostics)
511            DO j=jMinUpd,jMaxUpd            DO j=jMinUpd,jMaxUpd
            DO i=1-OLx+1,sNx+OLx-1  
             tmpVol = iceVol(i,j)  
             iceVol(i,j) = iceVol(i,j)  
      &                   -deltaTadvect*maskOc(i,j)  
      &                          *( uTrans(i+1,j)-uTrans(i,j) )  
             IF ( iceVol(i,j).GT.iceEps )  
      &      iceFld(i,j) = ( iceFld(i,j)*tmpVol  
      &                     -deltaTadvect*maskOc(i,j)  
      &                                  *( af(i+1,j)-af(i,j) )  
      &                    )/iceVol(i,j)  
            ENDDO  
           ENDDO  
          ENDIF  
 C--   keep advective flux (for diagnostics)  
          DO j=jMinUpd,jMaxUpd  
512             DO i=1-OLx+1,sNx+OLx             DO i=1-OLx+1,sNx+OLx
513              afx(i,j) = af(i,j)               afx(i,j) = af(i,j)
514             ENDDO             ENDDO
515           ENDDO            ENDDO
516    
517  C-     end if/else update overlap-Only  C-     end if/else update overlap-Only
518          ENDIF          ENDIF
# Line 508  CADJ STORE af(:,:)     = comlev1_thsice_ Line 538  CADJ STORE af(:,:)     = comlev1_thsice_
538  C-     Advective flux in Y  C-     Advective flux in Y
539           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
540            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
541             af(i,j) = 0.              af(i,j) = 0.
542            ENDDO            ENDDO
543           ENDDO           ENDDO
544    
# Line 517  C-     Internal exchange for calculation Line 547  C-     Internal exchange for calculation
547           IF ( overlapOnly ) THEN           IF ( overlapOnly ) THEN
548            CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,            CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
549       &                               iceFld, bi,bj, myThid )       &                               iceFld, bi,bj, myThid )
550            IF ( .NOT.useGridVolume )            IF ( .NOT.useGridArea )
551       &    CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,       &    CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
552       &                               iceVol, bi,bj, myThid )       &                               iceVol, bi,bj, myThid )
553           ENDIF           ENDIF
554  cph-exch2#endif  cph-exch2#endif
555    
556  C-     Compute CFL number  C-     Compute CFL number
557           IF ( useGridVolume ) THEN           IF ( useGridArea ) THEN
558            DO j=2-Oly,sNy+Oly            DO j=2-Oly,sNy+Oly
559             DO i=1-Olx,sNx+Olx             DO i=1-Olx,sNx+Olx
560              vCfl(i,j) = deltaTadvect*(              vCfl(i,j) = deltaTadvect*(
# Line 549  C-     Compute CFL number Line 579  C-     Compute CFL number
579            CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme, .FALSE.,            CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme, .FALSE.,
580       I             deltaTadvect, vTrans, vCfl, iceFld,       I             deltaTadvect, vTrans, vCfl, iceFld,
581       O             af, myThid )       O             af, myThid )
582    #ifdef ALLOW_DBUG_THSICE
583            IF (dBugFlag.AND.bi.EQ.biDb ) WRITE(6,'(A,1P4E14.6)')            IF (dBugFlag.AND.bi.EQ.biDb ) WRITE(6,'(A,1P4E14.6)')
584       &      'ICE_adv: yFx=', af(idb,jdb), iceFld(idb,jdb),       &      'ICE_adv: yFx=', af(idb,jdb), iceFld(idb,jdb),
585       &       vTrans(idb,jdb), af(idb,jdb)/vTrans(idb,jdb)       &       vTrans(idb,jdb), af(idb,jdb)/vTrans(idb,jdb)
586    #endif /* ALLOW_DBUG_THSICE */
587           ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN           ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
588            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, .FALSE., deltaTadvect,            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, .FALSE., deltaTadvect,
589       I             vTrans, vCfl, maskLocS, iceFld,       I             vTrans, vCfl, maskLocS, iceFld,
# Line 573  cph-exch2#ifndef ALLOW_AUTODIFF_TAMC Line 605  cph-exch2#ifndef ALLOW_AUTODIFF_TAMC
605           IF ( overlapOnly .AND. ipass.EQ.1 ) THEN           IF ( overlapOnly .AND. ipass.EQ.1 ) THEN
606            CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,            CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
607       &                               iceFld, bi,bj, myThid )       &                               iceFld, bi,bj, myThid )
608            IF ( .NOT.useGridVolume )            IF ( .NOT.useGridArea )
609       &    CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,       &    CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
610       &                               iceVol, bi,bj, myThid )       &                               iceVol, bi,bj, myThid )
611           ENDIF           ENDIF
# Line 597  C         in corner region) but safer to Line 629  C         in corner region) but safer to
629  CADJ STORE iceFld(:,:) = comlev1_thsice_4, key=ikey_4, byte=isbyte  CADJ STORE iceFld(:,:) = comlev1_thsice_4, key=ikey_4, byte=isbyte
630  CADJ STORE iceVol(:,:) = comlev1_thsice_4, key=ikey_4, byte=isbyte  CADJ STORE iceVol(:,:) = comlev1_thsice_4, key=ikey_4, byte=isbyte
631  #endif  #endif
          IF ( W_edge .AND. useGridVolume ) THEN  
           DO j=jMinUpd,jMaxUpd  
            DO i=1-OLx,0  
             iceFld(i,j) = iceFld(i,j)  
      &                   -deltaTadvect*maskOc(i,j)  
      &                                *recip_rA(i,j,bi,bj)  
      &                                *( af(i,j+1)-af(i,j) )  
            ENDDO  
           ENDDO  
          ELSEIF ( W_edge ) THEN  
           DO j=jMinUpd,jMaxUpd  
            DO i=1-OLx,0  
             tmpVol = iceVol(i,j)  
             iceVol(i,j) = iceVol(i,j)  
      &                   -deltaTadvect*maskOc(i,j)  
      &                          *( vTrans(i,j+1)-vTrans(i,j) )  
             IF ( iceVol(i,j).GT.iceEps )  
      &      iceFld(i,j) = ( iceFld(i,j)*tmpVol  
      &                     -deltaTadvect*maskOc(i,j)  
      &                                  *( af(i,j+1)-af(i,j) )  
      &                    )/iceVol(i,j)  
            ENDDO  
           ENDDO  
          ENDIF  
          IF ( E_edge .AND. useGridVolume ) THEN  
           DO j=jMinUpd,jMaxUpd  
            DO i=sNx+1,sNx+OLx  
             iceFld(i,j) = iceFld(i,j)  
      &                   -deltaTadvect*maskOc(i,j)  
      &                                *recip_rA(i,j,bi,bj)  
      &                                *( af(i,j+1)-af(i,j) )  
            ENDDO  
           ENDDO  
          ELSEIF ( E_edge ) THEN  
           DO j=jMinUpd,jMaxUpd  
            DO i=sNx+1,sNx+OLx  
             tmpVol = iceVol(i,j)  
             iceVol(i,j) = iceVol(i,j)  
      &                   -deltaTadvect*maskOc(i,j)  
      &                          *( vTrans(i,j+1)-vTrans(i,j) )  
             IF ( iceVol(i,j).GT.iceEps )  
      &      iceFld(i,j) = ( iceFld(i,j)*tmpVol  
      &                     -deltaTadvect*maskOc(i,j)  
      &                                  *( af(i,j+1)-af(i,j) )  
      &                    )/iceVol(i,j)  
            ENDDO  
           ENDDO  
          ENDIF  
 C--   keep advective flux (for diagnostics)  
632           IF ( W_edge ) THEN           IF ( W_edge ) THEN
633              IF ( useGridArea ) THEN
634               DO j=jMinUpd,jMaxUpd
635                DO i=1-OLx,0
636                 iceFld(i,j) = iceFld(i,j)
637         &                    -deltaTadvect*maskLocC(i,j)
638         &                                 *recip_rA(i,j,bi,bj)
639         &                                 *( af(i,j+1)-af(i,j) )
640                ENDDO
641               ENDDO
642              ELSE
643               DO j=jMinUpd,jMaxUpd
644                DO i=1-OLx,0
645                 tmpVol = iceVol(i,j)
646                 iceVol(i,j) = iceVol(i,j)
647         &                    -deltaTadvect*maskLocC(i,j)
648         &                           *( vTrans(i,j+1)-vTrans(i,j) )
649                 IF ( iceVol(i,j).GT.iceEps )
650         &       iceFld(i,j) = ( iceFld(i,j)*tmpVol
651         &                      -deltaTadvect*maskLocC(i,j)
652         &                                   *( af(i,j+1)-af(i,j) )
653         &                     )/iceVol(i,j)
654                ENDDO
655               ENDDO
656              ENDIF
657    C-    keep advective flux (for diagnostics)
658            DO j=1-OLy+1,sNy+OLy            DO j=1-OLy+1,sNy+OLy
659             DO i=1-OLx,0             DO i=1-OLx,0
660              afy(i,j) = af(i,j)               afy(i,j) = af(i,j)
661             ENDDO             ENDDO
662            ENDDO            ENDDO
663    C-    end if West Edge
664           ENDIF           ENDIF
665           IF ( E_edge ) THEN           IF ( E_edge ) THEN
666              IF ( useGridArea ) THEN
667               DO j=jMinUpd,jMaxUpd
668                DO i=sNx+1,sNx+OLx
669                 iceFld(i,j) = iceFld(i,j)
670         &                    -deltaTadvect*maskLocC(i,j)
671         &                                 *recip_rA(i,j,bi,bj)
672         &                                 *( af(i,j+1)-af(i,j) )
673                ENDDO
674               ENDDO
675              ELSE
676               DO j=jMinUpd,jMaxUpd
677                DO i=sNx+1,sNx+OLx
678                 tmpVol = iceVol(i,j)
679                 iceVol(i,j) = iceVol(i,j)
680         &                    -deltaTadvect*maskLocC(i,j)
681         &                           *( vTrans(i,j+1)-vTrans(i,j) )
682                 IF ( iceVol(i,j).GT.iceEps )
683         &       iceFld(i,j) = ( iceFld(i,j)*tmpVol
684         &                      -deltaTadvect*maskLocC(i,j)
685         &                                   *( af(i,j+1)-af(i,j) )
686         &                     )/iceVol(i,j)
687                ENDDO
688               ENDDO
689              ENDIF
690    C-    keep advective flux (for diagnostics)
691            DO j=1-OLy+1,sNy+OLy            DO j=1-OLy+1,sNy+OLy
692             DO i=sNx+1,sNx+OLx             DO i=sNx+1,sNx+OLx
693              afy(i,j) = af(i,j)               afy(i,j) = af(i,j)
694             ENDDO             ENDDO
695            ENDDO            ENDDO
696    C-    end if East Edge
697           ENDIF           ENDIF
698    
699          ELSE          ELSE
700  C     do not only update the overlap  C     do not only update the overlap
701           iMinUpd = 1-OLx            iMinUpd = 1-OLx
702           iMaxUpd = sNx+OLx            iMaxUpd = sNx+OLx
703           IF ( interiorOnly .AND. W_edge ) iMinUpd = 1            IF ( interiorOnly .AND. W_edge ) iMinUpd = 1
704           IF ( interiorOnly .AND. E_edge ) iMaxUpd = sNx            IF ( interiorOnly .AND. E_edge ) iMaxUpd = sNx
705           IF ( useGridVolume ) THEN            IF ( useGridArea ) THEN
706            DO j=1-OLy+1,sNy+OLy-1             DO j=1-OLy+1,sNy+OLy-1
707             DO i=iMinUpd,iMaxUpd              DO i=iMinUpd,iMaxUpd
708              iceFld(i,j) = iceFld(i,j)               iceFld(i,j) = iceFld(i,j)
709       &                   -deltaTadvect*maskOc(i,j)       &                    -deltaTadvect*maskLocC(i,j)
710       &                                *recip_rA(i,j,bi,bj)       &                                 *recip_rA(i,j,bi,bj)
711       &                                *( af(i,j+1)-af(i,j) )       &                                 *( af(i,j+1)-af(i,j) )
712                ENDDO
713               ENDDO
714              ELSE
715               DO j=1-OLy+1,sNy+OLy-1
716                DO i=iMinUpd,iMaxUpd
717                 tmpVol = iceVol(i,j)
718                 iceVol(i,j) = iceVol(i,j)
719         &                    -deltaTadvect*maskLocC(i,j)
720         &                           *( vTrans(i,j+1)-vTrans(i,j) )
721                 IF ( iceVol(i,j).GT.iceEps )
722         &       iceFld(i,j) = ( iceFld(i,j)*tmpVol
723         &                      -deltaTadvect*maskLocC(i,j)
724         &                                   *( af(i,j+1)-af(i,j) )
725         &                     )/iceVol(i,j)
726                ENDDO
727             ENDDO             ENDDO
728            ENDDO            ENDIF
729           ELSE  C-    keep advective flux (for diagnostics)
730            DO j=1-OLy+1,sNy+OLy-1            DO j=1-OLy+1,sNy+OLy
731             DO i=iMinUpd,iMaxUpd             DO i=iMinUpd,iMaxUpd
732              tmpVol = iceVol(i,j)               afy(i,j) = af(i,j)
             iceVol(i,j) = iceVol(i,j)  
      &                   -deltaTadvect*maskOc(i,j)  
      &                          *( vTrans(i,j+1)-vTrans(i,j) )  
             IF ( iceVol(i,j).GT.iceEps )  
      &      iceFld(i,j) = ( iceFld(i,j)*tmpVol  
      &                     -deltaTadvect*maskOc(i,j)  
      &                                  *( af(i,j+1)-af(i,j) )  
      &                    )/iceVol(i,j)  
733             ENDDO             ENDDO
734            ENDDO            ENDDO
          ENDIF  
 C--   keep advective flux (for diagnostics)  
          DO j=1-OLy+1,sNy+OLy  
            DO i=iMinUpd,iMaxUpd  
             afy(i,j) = af(i,j)  
            ENDDO  
          ENDDO  
735    
736  C      end if/else update overlap-Only  C      end if/else update overlap-Only
737          ENDIF          ENDIF
# Line 708  C--   End of ipass loop Line 743  C--   End of ipass loop
743        ENDDO        ENDDO
744    
745  C-    explicit advection is done ; add some debugging print  C-    explicit advection is done ; add some debugging print
746    #ifdef ALLOW_DBUG_THSICE
747        IF ( dBugFlag ) THEN        IF ( dBugFlag ) THEN
748         DO j=1-OLy,sNy+OLy         DO j=1-OLy,sNy+OLy
749          DO i=1-OLx,sNx+OLx          DO i=1-OLx,sNx+OLx
# Line 731  C-    explicit advection is done ; add s Line 767  C-    explicit advection is done ; add s
767       &      afx,afy, k, standardMessageUnit,bi,bj,myThid )       &      afx,afy, k, standardMessageUnit,bi,bj,myThid )
768        ENDIF        ENDIF
769  #endif /* ALLOW_DEBUG */  #endif /* ALLOW_DEBUG */
770    #endif /* ALLOW_DBUG_THSICE */
771    
772  #endif /* ALLOW_GENERIC_ADVDIFF */  #endif /* ALLOW_GENERIC_ADVDIFF */
773    
Legend:
Removed from v.1.14  
changed lines
  Added in v.1.15
  ViewVC Help
Powered by ViewVC 1.1.22