Back to home page

MITgcm
 
 

    

File indexing completed on 2023-02-03 06:10:33 UTC

view on githubraw file Latest commit edb66560 on 2023-02-02 23:32:31 UTC
a4eca6e929 Jean*0001 #include "THSICE_OPTIONS.h"
                0002 #ifdef ALLOW_GENERIC_ADVDIFF
                0003 # include "GAD_OPTIONS.h"
                0004 #endif
6b47d550f4 Mart*0005 #ifdef ALLOW_AUTODIFF
                0006 # include "AUTODIFF_OPTIONS.h"
                0007 #endif
a4eca6e929 Jean*0008 
                0009 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0010 CBOP
                0011 C !ROUTINE: THSICE_ADVECTION
                0012 
                0013 C !INTERFACE: ==========================================================
                0014       SUBROUTINE THSICE_ADVECTION(
                0015      I     tracerIdentity,
                0016      I     advectionScheme,
204a3108f8 Jean*0017      I     useGridArea,
a4eca6e929 Jean*0018      I     uTrans, vTrans, maskOc, deltaTadvect, iceEps,
                0019      U     iceVol, iceFld,
                0020      O     afx, afy,
                0021      I     bi, bj, myTime, myIter, myThid)
                0022 
                0023 C !DESCRIPTION:
                0024 C Calculates the tendency of a sea-ice field due to advection.
                0025 C It uses the multi-dimensional method given in \ref{sect:multiDimAdvection}
                0026 C and can only be used for the non-linear advection schemes such as the
                0027 C direct-space-time method and flux-limiters.
                0028 C
                0029 C This routine is an adaption of the GAD_ADVECTION for 2D-fields.
204a3108f8 Jean*0030 C for Area, effective thickness or other sea-ice field,
a4eca6e929 Jean*0031 C  the contribution iceFld*div(u) (that is present in gad_advection)
                0032 C  is not included here.
                0033 C
                0034 C The algorithm is as follows:
                0035 C \begin{itemize}
                0036 C \item{$\theta^{(n+1/2)} = \theta^{(n)}
                0037 C      - \Delta t \partial_x (u\theta^{(n)}) + \theta^{(n)} \partial_x u$}
                0038 C \item{$\theta^{(n+2/2)} = \theta^{(n+1/2)}
                0039 C      - \Delta t \partial_y (v\theta^{(n+1/2)}) + \theta^{(n)} \partial_y v$}
                0040 C \item{$G_\theta = ( \theta^{(n+2/2)} - \theta^{(n)} )/\Delta t$}
                0041 C \end{itemize}
                0042 C
                0043 C The tendency (output) is over-written by this routine.
                0044 
                0045 C !USES: ===============================================================
                0046       IMPLICIT NONE
                0047 #include "SIZE.h"
                0048 #include "EEPARAMS.h"
                0049 #include "PARAMS.h"
                0050 #include "GRID.h"
204a3108f8 Jean*0051 #if ( defined ALLOW_DBUG_THSICE || defined ALLOW_AUTODIFF_TAMC )
                0052 # include "THSICE_SIZE.h"
                0053 #endif
a4eca6e929 Jean*0054 #ifdef ALLOW_GENERIC_ADVDIFF
                0055 # include "GAD.h"
                0056 #endif
                0057 #ifdef ALLOW_EXCH2
f9f661930b Jean*0058 #include "W2_EXCH2_SIZE.h"
a4eca6e929 Jean*0059 #include "W2_EXCH2_TOPOLOGY.h"
                0060 #endif /* ALLOW_EXCH2 */
d6f06800ae Patr*0061 #ifdef ALLOW_AUTODIFF_TAMC
                0062 # include "tamc.h"
                0063 #endif
a4eca6e929 Jean*0064 
                0065 C !INPUT PARAMETERS: ===================================================
                0066 C  tracerIdentity  :: tracer identifier (required only for OBCS)
                0067 C  advectionScheme :: advection scheme to use (Horizontal plane)
204a3108f8 Jean*0068 C  useGridArea     :: use grid-cell Area (instead of "iceVol" field)
a4eca6e929 Jean*0069 C  uTrans,vTrans   :: volume transports at U,V points
                0070 C  maskOc          :: oceanic mask
                0071 C  iceVol          :: sea-ice volume
                0072 C  iceFld          :: sea-ice field
                0073 C  deltaTadvect    :: time-step used for advection [s]
                0074 C  iceEps          :: small volume (to avoid division by zero if iceVol==0)
                0075 C  bi,bj           :: tile indices
                0076 C  myTime          :: current time in simulation [s]
                0077 C  myIter          :: current iteration number
                0078 C  myThid          :: my thread Id. number
                0079       INTEGER tracerIdentity
                0080       INTEGER advectionScheme
204a3108f8 Jean*0081       LOGICAL useGridArea
a4eca6e929 Jean*0082       _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0083       _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0084       _RS maskOc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0085       _RL iceFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0086       _RL iceVol(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0087       _RL deltaTadvect, iceEps
                0088       INTEGER bi,bj
                0089       _RL myTime
                0090       INTEGER myIter
                0091       INTEGER myThid
                0092 
                0093 C !OUTPUT PARAMETERS: ==================================================
                0094 C  iceVol (Updated):: sea-ice volume
                0095 C  iceFld (Updated):: sea-ice field
                0096 C  afx             :: horizontal advective flux, x direction
                0097 C  afy             :: horizontal advective flux, y direction
                0098       _RL afx   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0099       _RL afy   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0100 
e7f517b398 Jean*0101 #ifdef ALLOW_GENERIC_ADVDIFF
a4eca6e929 Jean*0102 C !LOCAL VARIABLES: ====================================================
204a3108f8 Jean*0103 C  maskLocC      :: 2-D array mask at grid-cell center
a4eca6e929 Jean*0104 C  maskLocW      :: 2-D array for mask at West points
                0105 C  maskLocS      :: 2-D array for mask at South points
                0106 C [iMin,iMax]Upd :: loop range to update sea-ice field
                0107 C [jMin,jMax]Upd :: loop range to update sea-ice field
                0108 C  i,j           :: loop indices
                0109 C  uCfl          :: CFL number, zonal direction
                0110 C  vCfl          :: CFL number, meridional direction
                0111 C  af            :: 2-D array for horizontal advective flux
                0112 C  calc_fluxes_X :: logical to indicate to calculate fluxes in X dir
                0113 C  calc_fluxes_Y :: logical to indicate to calculate fluxes in Y dir
                0114 C  interiorOnly  :: only update the interior of myTile, but not the edges
                0115 C  overlapOnly   :: only update the edges of myTile, but not the interior
                0116 C  nipass        :: number of passes in multi-dimensional method
                0117 C  ipass         :: number of the current pass being made
                0118 C  myTile        :: variables used to determine which cube face
                0119 C  nCFace        :: owns a tile for cube grid runs using
                0120 C                :: multi-dim advection.
                0121 C [N,S,E,W]_edge :: true if N,S,E,W edge of myTile is an Edge of the cube
204a3108f8 Jean*0122       _RS maskLocC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
a4eca6e929 Jean*0123       _RS maskLocW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0124       _RS maskLocS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0125       INTEGER iMinUpd,iMaxUpd,jMinUpd,jMaxUpd
                0126       INTEGER i,j,k
                0127       _RL uCfl    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0128       _RL vCfl    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0129       _RL af      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
204a3108f8 Jean*0130       _RL tmpVol
a4eca6e929 Jean*0131       LOGICAL calc_fluxes_X, calc_fluxes_Y, withSigns
                0132       LOGICAL interiorOnly, overlapOnly
                0133       INTEGER nipass,ipass
                0134       INTEGER nCFace
                0135       LOGICAL N_edge, S_edge, E_edge, W_edge
                0136 #ifdef ALLOW_EXCH2
                0137       INTEGER myTile
                0138 #endif
7c50f07931 Mart*0139 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0140 C     trNumber :: thsice tracer number
                0141 C     tkey     :: tape key (depends on tile indices and tracer index)
                0142 C     dkey     :: tape key (depends on direction and tkey)
                0143       INTEGER trNumber, tkey, dkey
a85293d087 Mart*0144       CHARACTER*(MAX_LEN_MBUF) msgBuf
7c50f07931 Mart*0145 #endif
204a3108f8 Jean*0146 #ifdef ALLOW_DBUG_THSICE
fef0c02f19 Jean*0147       LOGICAL dBugFlag
                0148       INTEGER idb,jdb,biDb
a4eca6e929 Jean*0149       _RL tmpFac
204a3108f8 Jean*0150 #endif /* ALLOW_DBUG_THSICE */
a4eca6e929 Jean*0151 CEOP
                0152 
d6f06800ae Patr*0153 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0154 
204a3108f8 Jean*0155       k = 1
d6f06800ae Patr*0156 #ifdef ALLOW_AUTODIFF_TAMC
a85293d087 Mart*0157 C     tracerIdentity has negative values from GAD_SI_FRAC=-5 to GAD_SI_QICE2=-9
edb6656069 Mart*0158       trNumber = GAD_SI_FRAC - tracerIdentity
                0159       tkey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
                0160       tkey = trNumber + 1 + (tkey-1)*thSIce_nAdv
                0161       IF ( trNumber .LT. 0 .OR. trNumber+1 .GT. thSIce_nAdv) THEN
a85293d087 Mart*0162        WRITE(msgBuf,'(A,2(I3,A))')
                0163      &      'THSICE_ADVECTION: tracerIdentity =', tracerIdentity,
                0164      &      ' not consistent with thSIce_nAdv =', thSIce_nAdv,
                0165      &      ', ==> check "THSICE_SIZE.h"'
                0166        CALL PRINT_ERROR( msgBuf, myThid )
                0167        STOP 'ABNORMAL END: S/R THSICE_ADVECTION'
                0168       ENDIF
d6f06800ae Patr*0169 #endif /* ALLOW_AUTODIFF_TAMC */
                0170 
204a3108f8 Jean*0171 #ifdef ALLOW_DBUG_THSICE
ae605e558b Jean*0172       dBugFlag = debugLevel.GE.debLevC
a4eca6e929 Jean*0173      &     .AND. myIter.EQ.nIter0
                0174      &     .AND. ( tracerIdentity.EQ.GAD_SI_HICE .OR.
                0175      &             tracerIdentity.EQ.GAD_SI_QICE2 )
fef0c02f19 Jean*0176 c    &     .AND. tracerIdentity.EQ.GAD_SI_FRAC
                0177       idb  = MIN(13,sNx)
                0178       jdb  = MOD(17,sNy)
                0179       biDb = nSx/2
204a3108f8 Jean*0180 #endif /* ALLOW_DBUG_THSICE */
a4eca6e929 Jean*0181 
                0182 C--   Set up work arrays with valid (i.e. not NaN) values
                0183 C     These inital values do not alter the numerical results. They
                0184 C     just ensure that all memory references are to valid floating
                0185 C     point numbers. This prevents spurious hardware signals due to
                0186 C     uninitialised but inert locations.
                0187 
                0188 C--   Set tile-specific parameters for horizontal fluxes
                0189       IF (useCubedSphereExchange) THEN
                0190        nipass=3
                0191 #ifdef ALLOW_EXCH2
c424ee7cc7 Jean*0192        myTile = W2_myTileList(bi,bj)
a4eca6e929 Jean*0193        nCFace = exch2_myFace(myTile)
                0194        N_edge = exch2_isNedge(myTile).EQ.1
                0195        S_edge = exch2_isSedge(myTile).EQ.1
                0196        E_edge = exch2_isEedge(myTile).EQ.1
                0197        W_edge = exch2_isWedge(myTile).EQ.1
                0198 #else
                0199        nCFace = bi
                0200        N_edge = .TRUE.
                0201        S_edge = .TRUE.
                0202        E_edge = .TRUE.
                0203        W_edge = .TRUE.
                0204 #endif
                0205       ELSE
                0206        nipass=2
                0207        nCFace = bi
                0208        N_edge = .FALSE.
                0209        S_edge = .FALSE.
                0210        E_edge = .FALSE.
                0211        W_edge = .FALSE.
                0212       ENDIF
                0213 
a85293d087 Mart*0214 #ifdef ALLOW_AUTODIFF_TAMC
                0215       IF (nipass .GT. maxcube) THEN
                0216        WRITE(msgBuf,'(A,2(I3,A))') 'S/R THSICE_ADVECTION: nipass =',
                0217      &      nipass, ' >', maxcube, ' = maxcube, ==> check "tamc.h"'
                0218        CALL PRINT_ERROR( msgBuf, myThid )
                0219        STOP 'ABNORMAL END: S/R THSICE_ADVECTION'
                0220       ENDIF
                0221 #endif
a4eca6e929 Jean*0222 
                0223 C--   Start horizontal fluxes
                0224 
204a3108f8 Jean*0225 C--   set mask West & South (and local Centered mask)
a4eca6e929 Jean*0226       DO j=1-OLy,sNy+OLy
6b47d550f4 Mart*0227        maskLocW(1-OLx,j) = 0.
a4eca6e929 Jean*0228        DO i=2-OLx,sNx+OLx
                0229         maskLocW(i,j) = MIN( maskOc(i-1,j), maskOc(i,j) )
204a3108f8 Jean*0230 #ifdef ALLOW_OBCS
                0231      &                * maskInW(i,j,bi,bj)
                0232 #endif /* ALLOW_OBCS */
a4eca6e929 Jean*0233        ENDDO
                0234       ENDDO
                0235       DO i=1-OLx,sNx+OLx
6b47d550f4 Mart*0236        maskLocS(i,1-OLy) = 0.
a4eca6e929 Jean*0237       ENDDO
                0238       DO j=2-OLy,sNy+OLy
                0239        DO i=1-OLx,sNx+OLx
                0240         maskLocS(i,j) = MIN( maskOc(i,j-1), maskOc(i,j) )
204a3108f8 Jean*0241 #ifdef ALLOW_OBCS
                0242      &                * maskInS(i,j,bi,bj)
                0243 #endif /* ALLOW_OBCS */
                0244        ENDDO
                0245       ENDDO
                0246 C     maskLocC is just a local copy of Ocean mask (maskOc) except if using OBCS:
                0247 C     use "maksInC" to prevent updating tracer field in OB regions
                0248       DO j=1-OLy,sNy+OLy
                0249        DO i=1-OLx,sNx+OLx
                0250         maskLocC(i,j) = maskOc(i,j)
                0251 #ifdef ALLOW_OBCS
                0252      &                * maskInC(i,j,bi,bj)
                0253 #endif /* ALLOW_OBCS */
a4eca6e929 Jean*0254        ENDDO
                0255       ENDDO
                0256 
                0257       IF (useCubedSphereExchange) THEN
                0258        withSigns = .FALSE.
                0259        CALL FILL_CS_CORNER_UV_RS(
                0260      &      withSigns, maskLocW,maskLocS, bi,bj, myThid )
                0261       ENDIF
                0262 
                0263 C--   Multiple passes for different directions on different tiles
                0264 C--   For cube need one pass for each of red, green and blue axes.
                0265       DO ipass=1,nipass
d6f06800ae Patr*0266 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0267        dkey = ipass + (tkey-1) * maxcube
d6f06800ae Patr*0268 #endif
a4eca6e929 Jean*0269 
                0270        interiorOnly = .FALSE.
                0271        overlapOnly  = .FALSE.
                0272        IF (useCubedSphereExchange) THEN
                0273 C--   CubedSphere : pass 3 times, with partial update of local seaice field
                0274         IF (ipass.EQ.1) THEN
                0275          overlapOnly  = MOD(nCFace,3).EQ.0
                0276          interiorOnly = MOD(nCFace,3).NE.0
                0277          calc_fluxes_X = nCFace.EQ.6 .OR. nCFace.EQ.1 .OR. nCFace.EQ.2
                0278          calc_fluxes_Y = nCFace.EQ.3 .OR. nCFace.EQ.4 .OR. nCFace.EQ.5
                0279         ELSEIF (ipass.EQ.2) THEN
                0280          overlapOnly  = MOD(nCFace,3).EQ.2
e3829af00c Jean*0281          interiorOnly = MOD(nCFace,3).EQ.1
a4eca6e929 Jean*0282          calc_fluxes_X = nCFace.EQ.2 .OR. nCFace.EQ.3 .OR. nCFace.EQ.4
                0283          calc_fluxes_Y = nCFace.EQ.5 .OR. nCFace.EQ.6 .OR. nCFace.EQ.1
                0284         ELSE
e3829af00c Jean*0285          interiorOnly = .TRUE.
a4eca6e929 Jean*0286          calc_fluxes_X = nCFace.EQ.5 .OR. nCFace.EQ.6
                0287          calc_fluxes_Y = nCFace.EQ.2 .OR. nCFace.EQ.3
                0288         ENDIF
                0289        ELSE
                0290 C--   not CubedSphere
                0291         calc_fluxes_X = MOD(ipass,2).EQ.1
                0292         calc_fluxes_Y = .NOT.calc_fluxes_X
                0293        ENDIF
204a3108f8 Jean*0294 #ifdef ALLOW_DBUG_THSICE
fef0c02f19 Jean*0295        IF (dBugFlag.AND.bi.EQ.biDb ) WRITE(6,'(A,3I4,2I5,4L5)')
                0296      &   'ICE_adv:', tracerIdentity, ipass, bi, idb, jdb,
                0297      &   calc_fluxes_X, calc_fluxes_Y, overlapOnly, interiorOnly
204a3108f8 Jean*0298 #endif /* ALLOW_DBUG_THSICE */
a4eca6e929 Jean*0299 
                0300 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0301 C--   X direction
                0302 
                0303        IF (calc_fluxes_X) THEN
                0304 
                0305 C-     Do not compute fluxes if
                0306 C       a) needed in overlap only
                0307 C   and b) the overlap of myTile are not cube-face Edges
d6f06800ae Patr*0308 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0309 CADJ STORE iceFld = comlev1_thsice_adv, key=dkey, byte=isbyte
                0310 CADJ STORE iceVol = comlev1_thsice_adv, key=dkey, byte=isbyte
d6f06800ae Patr*0311 #endif
a4eca6e929 Jean*0312         IF ( .NOT.overlapOnly .OR. N_edge .OR. S_edge ) THEN
                0313 
                0314 C-     Advective flux in X
6b47d550f4 Mart*0315          DO j=1-OLy,sNy+OLy
                0316           DO i=1-OLx,sNx+OLx
204a3108f8 Jean*0317             af(i,j) = 0.
a4eca6e929 Jean*0318           ENDDO
                0319          ENDDO
                0320 
                0321 C-     Internal exchange for calculations in X
e3829af00c Jean*0322          IF ( overlapOnly ) THEN
93e3461d85 Jean*0323           CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
1891130b05 Jean*0324      &                               iceFld, bi,bj, myThid )
204a3108f8 Jean*0325           IF ( .NOT.useGridArea )
93e3461d85 Jean*0326      &    CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
1891130b05 Jean*0327      &                               iceVol, bi,bj, myThid )
a4eca6e929 Jean*0328          ENDIF
                0329 
                0330 C-     Compute CFL number
204a3108f8 Jean*0331          IF ( useGridArea ) THEN
6b47d550f4 Mart*0332           DO j=1-OLy,sNy+OLy
                0333            DO i=2-OLx,sNx+OLx
a4eca6e929 Jean*0334             uCfl(i,j) = deltaTadvect*(
                0335      &          MAX( uTrans(i,j), 0. _d 0 )*recip_rA(i-1,j,bi,bj)
                0336      &         +MAX(-uTrans(i,j), 0. _d 0 )*recip_rA( i ,j,bi,bj)
                0337      &                               )
                0338            ENDDO
                0339           ENDDO
                0340          ELSE
6b47d550f4 Mart*0341           DO j=1-OLy,sNy+OLy
                0342            DO i=2-OLx,sNx+OLx
a4eca6e929 Jean*0343             uCfl(i,j) = deltaTadvect*(
                0344      &        MAX( uTrans(i,j), 0. _d 0 )/MAX( iceVol(i-1,j), iceEps )
                0345      &       +MAX(-uTrans(i,j), 0. _d 0 )/MAX( iceVol( i ,j), iceEps )
                0346      &                               )
                0347            ENDDO
                0348           ENDDO
                0349          ENDIF
                0350 
                0351          IF ( advectionScheme.EQ.ENUM_UPWIND_1RST
                0352      &        .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
                0353           CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme, .FALSE.,
                0354      I             deltaTadvect, uTrans, uCfl, iceFld,
                0355      O             af, myThid )
204a3108f8 Jean*0356 #ifdef ALLOW_DBUG_THSICE
fef0c02f19 Jean*0357           IF (dBugFlag.AND.bi.EQ.biDb ) WRITE(6,'(A,1P4E14.6)')
                0358      &      'ICE_adv: xFx=', af(idb,jdb), iceFld(idb,jdb),
                0359      &       uTrans(idb,jdb), af(idb,jdb)/uTrans(idb,jdb)
204a3108f8 Jean*0360 #endif /* ALLOW_DBUG_THSICE */
a4eca6e929 Jean*0361          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
                0362           CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, .FALSE., deltaTadvect,
                0363      I             uTrans, uCfl, maskLocW, iceFld,
                0364      O             af, myThid )
                0365          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
                0366           CALL GAD_DST3_ADV_X(      bi,bj,k, .FALSE., deltaTadvect,
                0367      I             uTrans, uCfl, maskLocW, iceFld,
                0368      O             af, myThid )
                0369          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
                0370           CALL GAD_DST3FL_ADV_X(    bi,bj,k, .FALSE., deltaTadvect,
                0371      I             uTrans, uCfl, maskLocW, iceFld,
                0372      O             af, myThid )
                0373          ELSE
                0374           STOP
                0375      & 'THSICE_ADVECTION: adv. scheme incompatibale with multi-dim'
                0376          ENDIF
                0377 
                0378 C--   Internal exchange for next calculations in Y
e3829af00c Jean*0379          IF ( overlapOnly .AND. ipass.EQ.1 ) THEN
93e3461d85 Jean*0380           CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
e3829af00c Jean*0381      &                               iceFld, bi,bj, myThid )
204a3108f8 Jean*0382           IF ( .NOT.useGridArea )
93e3461d85 Jean*0383      &    CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
e3829af00c Jean*0384      &                               iceVol, bi,bj, myThid )
                0385          ENDIF
a4eca6e929 Jean*0386 
e3829af00c Jean*0387 C--   Advective flux in X : done
                0388         ENDIF
                0389 
a4eca6e929 Jean*0390 C-     Update the local seaice field where needed:
                0391 
                0392 C     update in overlap-Only
                0393         IF ( overlapOnly ) THEN
                0394          iMinUpd = 1-OLx+1
                0395          iMaxUpd = sNx+OLx-1
                0396 C--   notes: these 2 lines below have no real effect (because recip_hFac=0
                0397 C            in corner region) but safer to keep them.
                0398          IF ( W_edge ) iMinUpd = 1
                0399          IF ( E_edge ) iMaxUpd = sNx
                0400 
d6f06800ae Patr*0401 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0402 CADJ STORE iceFld = comlev1_thsice_adv, key=dkey, byte=isbyte
                0403 CADJ STORE iceVol = comlev1_thsice_adv, key=dkey, byte=isbyte
                0404 CADJ STORE af     = comlev1_thsice_adv, key=dkey, byte=isbyte
d6f06800ae Patr*0405 #endif
204a3108f8 Jean*0406          IF ( S_edge ) THEN
                0407           IF ( useGridArea ) THEN
                0408            DO j=1-OLy,0
                0409             DO i=iMinUpd,iMaxUpd
                0410              iceFld(i,j) = iceFld(i,j)
                0411      &                    -deltaTadvect*maskLocC(i,j)
                0412      &                                 *recip_rA(i,j,bi,bj)
                0413      &                                 *( af(i+1,j)-af(i,j) )
                0414             ENDDO
a4eca6e929 Jean*0415            ENDDO
204a3108f8 Jean*0416           ELSE
                0417            DO j=1-OLy,0
                0418             DO i=iMinUpd,iMaxUpd
                0419              tmpVol = iceVol(i,j)
                0420              iceVol(i,j) = iceVol(i,j)
                0421      &                    -deltaTadvect*maskLocC(i,j)
                0422      &                           *( uTrans(i+1,j)-uTrans(i,j) )
                0423              IF ( iceVol(i,j).GT.iceEps )
                0424      &       iceFld(i,j) = ( iceFld(i,j)*tmpVol
                0425      &                      -deltaTadvect*maskLocC(i,j)
a4eca6e929 Jean*0426      &                                  *( af(i+1,j)-af(i,j) )
204a3108f8 Jean*0427      &                     )/iceVol(i,j)
                0428             ENDDO
a4eca6e929 Jean*0429            ENDDO
204a3108f8 Jean*0430           ENDIF
                0431 C-    keep advective flux (for diagnostics)
a4eca6e929 Jean*0432           DO j=1-OLy,0
                0433            DO i=1-OLx+1,sNx+OLx
204a3108f8 Jean*0434              afx(i,j) = af(i,j)
a4eca6e929 Jean*0435            ENDDO
                0436           ENDDO
204a3108f8 Jean*0437 C-    end if South Edge
a4eca6e929 Jean*0438          ENDIF
                0439          IF ( N_edge ) THEN
204a3108f8 Jean*0440           IF ( useGridArea ) THEN
                0441            DO j=sNy+1,sNy+OLy
                0442             DO i=iMinUpd,iMaxUpd
                0443              iceFld(i,j) = iceFld(i,j)
                0444      &                    -deltaTadvect*maskLocC(i,j)
                0445      &                                 *recip_rA(i,j,bi,bj)
                0446      &                                 *( af(i+1,j)-af(i,j) )
                0447             ENDDO
                0448            ENDDO
                0449           ELSE
                0450            DO j=sNy+1,sNy+OLy
                0451             DO i=iMinUpd,iMaxUpd
                0452              tmpVol = iceVol(i,j)
                0453              iceVol(i,j) = iceVol(i,j)
                0454      &                    -deltaTadvect*maskLocC(i,j)
                0455      &                           *( uTrans(i+1,j)-uTrans(i,j) )
                0456              IF ( iceVol(i,j).GT.iceEps )
                0457      &       iceFld(i,j) = ( iceFld(i,j)*tmpVol
                0458      &                      -deltaTadvect*maskLocC(i,j)
                0459      &                                   *( af(i+1,j)-af(i,j) )
                0460      &                     )/iceVol(i,j)
                0461             ENDDO
                0462            ENDDO
                0463           ENDIF
                0464 C-    keep advective flux (for diagnostics)
a4eca6e929 Jean*0465           DO j=sNy+1,sNy+OLy
                0466            DO i=1-OLx+1,sNx+OLx
204a3108f8 Jean*0467              afx(i,j) = af(i,j)
a4eca6e929 Jean*0468            ENDDO
                0469           ENDDO
204a3108f8 Jean*0470 C-    end if North Edge
a4eca6e929 Jean*0471          ENDIF
                0472 
                0473         ELSE
a85293d087 Mart*0474 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0475 CADJ STORE iceFld = comlev1_thsice_adv, key=dkey, byte=isbyte
                0476 CADJ STORE iceVol = comlev1_thsice_adv, key=dkey, byte=isbyte
                0477 CADJ STORE af     = comlev1_thsice_adv, key=dkey, byte=isbyte
a85293d087 Mart*0478 #endif
a4eca6e929 Jean*0479 C     do not only update the overlap
204a3108f8 Jean*0480           jMinUpd = 1-OLy
                0481           jMaxUpd = sNy+OLy
                0482           IF ( interiorOnly .AND. S_edge ) jMinUpd = 1
                0483           IF ( interiorOnly .AND. N_edge ) jMaxUpd = sNy
                0484           IF ( useGridArea ) THEN
                0485            DO j=jMinUpd,jMaxUpd
                0486             DO i=1-OLx+1,sNx+OLx-1
                0487              iceFld(i,j) = iceFld(i,j)
                0488      &                    -deltaTadvect*maskLocC(i,j)
                0489      &                                 *recip_rA(i,j,bi,bj)
                0490      &                                 *( af(i+1,j)-af(i,j) )
                0491             ENDDO
a4eca6e929 Jean*0492            ENDDO
204a3108f8 Jean*0493           ELSE
                0494            DO j=jMinUpd,jMaxUpd
                0495             DO i=1-OLx+1,sNx+OLx-1
                0496              tmpVol = iceVol(i,j)
                0497              iceVol(i,j) = iceVol(i,j)
                0498      &                    -deltaTadvect*maskLocC(i,j)
                0499      &                           *( uTrans(i+1,j)-uTrans(i,j) )
                0500              IF ( iceVol(i,j).GT.iceEps )
                0501      &       iceFld(i,j) = ( iceFld(i,j)*tmpVol
                0502      &                      -deltaTadvect*maskLocC(i,j)
                0503      &                                   *( af(i+1,j)-af(i,j) )
                0504      &                     )/iceVol(i,j)
                0505             ENDDO
a4eca6e929 Jean*0506            ENDDO
204a3108f8 Jean*0507           ENDIF
                0508 C-    keep advective flux (for diagnostics)
                0509           DO j=jMinUpd,jMaxUpd
a4eca6e929 Jean*0510            DO i=1-OLx+1,sNx+OLx
204a3108f8 Jean*0511              afx(i,j) = af(i,j)
a4eca6e929 Jean*0512            ENDDO
204a3108f8 Jean*0513           ENDDO
a4eca6e929 Jean*0514 
                0515 C-     end if/else update overlap-Only
                0516         ENDIF
                0517 
                0518 C--   End of X direction
                0519        ENDIF
                0520 
                0521 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0522 C--   Y direction
                0523 
                0524        IF (calc_fluxes_Y) THEN
                0525 
                0526 C-     Do not compute fluxes if
                0527 C       a) needed in overlap only
                0528 C   and b) the overlap of myTile are not cube-face edges
d6f06800ae Patr*0529 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0530 CADJ STORE iceFld = comlev1_thsice_adv, key=dkey, byte=isbyte
                0531 CADJ STORE iceVol = comlev1_thsice_adv, key=dkey, byte=isbyte
d6f06800ae Patr*0532 #endif
a4eca6e929 Jean*0533         IF ( .NOT.overlapOnly .OR. E_edge .OR. W_edge ) THEN
                0534 
                0535 C-     Advective flux in Y
                0536          DO j=1-OLy,sNy+OLy
                0537           DO i=1-OLx,sNx+OLx
204a3108f8 Jean*0538             af(i,j) = 0.
a4eca6e929 Jean*0539           ENDDO
                0540          ENDDO
                0541 
                0542 C-     Internal exchange for calculations in Y
e3829af00c Jean*0543          IF ( overlapOnly ) THEN
93e3461d85 Jean*0544           CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
1891130b05 Jean*0545      &                               iceFld, bi,bj, myThid )
204a3108f8 Jean*0546           IF ( .NOT.useGridArea )
93e3461d85 Jean*0547      &    CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
1891130b05 Jean*0548      &                               iceVol, bi,bj, myThid )
a4eca6e929 Jean*0549          ENDIF
                0550 
                0551 C-     Compute CFL number
204a3108f8 Jean*0552          IF ( useGridArea ) THEN
6b47d550f4 Mart*0553           DO j=2-OLy,sNy+OLy
                0554            DO i=1-OLx,sNx+OLx
a4eca6e929 Jean*0555             vCfl(i,j) = deltaTadvect*(
                0556      &          MAX( vTrans(i,j), 0. _d 0 )*recip_rA(i,j-1,bi,bj)
                0557      &         +MAX(-vTrans(i,j), 0. _d 0 )*recip_rA(i, j ,bi,bj)
                0558      &                               )
                0559            ENDDO
                0560           ENDDO
                0561          ELSE
6b47d550f4 Mart*0562           DO j=2-OLy,sNy+OLy
                0563            DO i=1-OLx,sNx+OLx
a4eca6e929 Jean*0564             vCfl(i,j) = deltaTadvect*(
                0565      &        MAX( vTrans(i,j), 0. _d 0 )/MAX( iceVol(i,j-1), iceEps )
                0566      &       +MAX(-vTrans(i,j), 0. _d 0 )/MAX( iceVol(i, j ), iceEps )
                0567      &                               )
                0568            ENDDO
                0569           ENDDO
                0570          ENDIF
                0571 
                0572          IF ( advectionScheme.EQ.ENUM_UPWIND_1RST
                0573      &        .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
                0574           CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme, .FALSE.,
                0575      I             deltaTadvect, vTrans, vCfl, iceFld,
                0576      O             af, myThid )
204a3108f8 Jean*0577 #ifdef ALLOW_DBUG_THSICE
fef0c02f19 Jean*0578           IF (dBugFlag.AND.bi.EQ.biDb ) WRITE(6,'(A,1P4E14.6)')
                0579      &      'ICE_adv: yFx=', af(idb,jdb), iceFld(idb,jdb),
                0580      &       vTrans(idb,jdb), af(idb,jdb)/vTrans(idb,jdb)
204a3108f8 Jean*0581 #endif /* ALLOW_DBUG_THSICE */
a4eca6e929 Jean*0582          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
                0583           CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, .FALSE., deltaTadvect,
                0584      I             vTrans, vCfl, maskLocS, iceFld,
                0585      O             af, myThid )
                0586          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
                0587           CALL GAD_DST3_ADV_Y(      bi,bj,k, .FALSE., deltaTadvect,
                0588      I             vTrans, vCfl, maskLocS, iceFld,
                0589      O             af, myThid )
                0590          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
                0591           CALL GAD_DST3FL_ADV_Y(    bi,bj,k, .FALSE., deltaTadvect,
                0592      I             vTrans, vCfl, maskLocS, iceFld,
                0593      O             af, myThid )
                0594          ELSE
                0595           STOP
                0596      &  'THSICE_ADVECTION: adv. scheme incompatibale with mutli-dim'
                0597          ENDIF
                0598 
e3829af00c Jean*0599          IF ( overlapOnly .AND. ipass.EQ.1 ) THEN
93e3461d85 Jean*0600           CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
e3829af00c Jean*0601      &                               iceFld, bi,bj, myThid )
204a3108f8 Jean*0602           IF ( .NOT.useGridArea )
93e3461d85 Jean*0603      &    CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
e3829af00c Jean*0604      &                               iceVol, bi,bj, myThid )
                0605          ENDIF
a4eca6e929 Jean*0606 
e3829af00c Jean*0607 C-     Advective flux in Y : done
                0608         ENDIF
                0609 
a4eca6e929 Jean*0610 C-     Update the local seaice field where needed:
                0611 
                0612 C      update in overlap-Only
                0613         IF ( overlapOnly ) THEN
                0614          jMinUpd = 1-OLy+1
                0615          jMaxUpd = sNy+OLy-1
                0616 C- notes: these 2 lines below have no real effect (because recip_hFac=0
                0617 C         in corner region) but safer to keep them.
                0618          IF ( S_edge ) jMinUpd = 1
                0619          IF ( N_edge ) jMaxUpd = sNy
                0620 
d6f06800ae Patr*0621 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0622 CADJ STORE iceFld = comlev1_thsice_adv, key=dkey, byte=isbyte
                0623 CADJ STORE iceVol = comlev1_thsice_adv, key=dkey, byte=isbyte
                0624 CADJ STORE af     = comlev1_thsice_adv, key=dkey, byte=isbyte
d6f06800ae Patr*0625 #endif
204a3108f8 Jean*0626          IF ( W_edge ) THEN
                0627           IF ( useGridArea ) THEN
                0628            DO j=jMinUpd,jMaxUpd
                0629             DO i=1-OLx,0
                0630              iceFld(i,j) = iceFld(i,j)
                0631      &                    -deltaTadvect*maskLocC(i,j)
                0632      &                                 *recip_rA(i,j,bi,bj)
                0633      &                                 *( af(i,j+1)-af(i,j) )
                0634             ENDDO
a4eca6e929 Jean*0635            ENDDO
204a3108f8 Jean*0636           ELSE
                0637            DO j=jMinUpd,jMaxUpd
                0638             DO i=1-OLx,0
                0639              tmpVol = iceVol(i,j)
                0640              iceVol(i,j) = iceVol(i,j)
                0641      &                    -deltaTadvect*maskLocC(i,j)
                0642      &                           *( vTrans(i,j+1)-vTrans(i,j) )
                0643              IF ( iceVol(i,j).GT.iceEps )
                0644      &       iceFld(i,j) = ( iceFld(i,j)*tmpVol
                0645      &                      -deltaTadvect*maskLocC(i,j)
                0646      &                                   *( af(i,j+1)-af(i,j) )
                0647      &                     )/iceVol(i,j)
                0648             ENDDO
a4eca6e929 Jean*0649            ENDDO
204a3108f8 Jean*0650           ENDIF
                0651 C-    keep advective flux (for diagnostics)
a4eca6e929 Jean*0652           DO j=1-OLy+1,sNy+OLy
                0653            DO i=1-OLx,0
204a3108f8 Jean*0654              afy(i,j) = af(i,j)
a4eca6e929 Jean*0655            ENDDO
                0656           ENDDO
204a3108f8 Jean*0657 C-    end if West Edge
a4eca6e929 Jean*0658          ENDIF
                0659          IF ( E_edge ) THEN
204a3108f8 Jean*0660           IF ( useGridArea ) THEN
                0661            DO j=jMinUpd,jMaxUpd
                0662             DO i=sNx+1,sNx+OLx
                0663              iceFld(i,j) = iceFld(i,j)
                0664      &                    -deltaTadvect*maskLocC(i,j)
                0665      &                                 *recip_rA(i,j,bi,bj)
                0666      &                                 *( af(i,j+1)-af(i,j) )
                0667             ENDDO
                0668            ENDDO
                0669           ELSE
                0670            DO j=jMinUpd,jMaxUpd
                0671             DO i=sNx+1,sNx+OLx
                0672              tmpVol = iceVol(i,j)
                0673              iceVol(i,j) = iceVol(i,j)
                0674      &                    -deltaTadvect*maskLocC(i,j)
                0675      &                           *( vTrans(i,j+1)-vTrans(i,j) )
                0676              IF ( iceVol(i,j).GT.iceEps )
                0677      &       iceFld(i,j) = ( iceFld(i,j)*tmpVol
                0678      &                      -deltaTadvect*maskLocC(i,j)
                0679      &                                   *( af(i,j+1)-af(i,j) )
                0680      &                     )/iceVol(i,j)
                0681             ENDDO
                0682            ENDDO
                0683           ENDIF
                0684 C-    keep advective flux (for diagnostics)
a4eca6e929 Jean*0685           DO j=1-OLy+1,sNy+OLy
                0686            DO i=sNx+1,sNx+OLx
204a3108f8 Jean*0687              afy(i,j) = af(i,j)
a4eca6e929 Jean*0688            ENDDO
                0689           ENDDO
204a3108f8 Jean*0690 C-    end if East Edge
a4eca6e929 Jean*0691          ENDIF
                0692 
                0693         ELSE
a85293d087 Mart*0694 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0695 CADJ STORE iceFld = comlev1_thsice_adv, key=dkey, byte=isbyte
                0696 CADJ STORE iceVol = comlev1_thsice_adv, key=dkey, byte=isbyte
                0697 CADJ STORE af     = comlev1_thsice_adv, key=dkey, byte=isbyte
a85293d087 Mart*0698 #endif
a4eca6e929 Jean*0699 C     do not only update the overlap
204a3108f8 Jean*0700           iMinUpd = 1-OLx
                0701           iMaxUpd = sNx+OLx
                0702           IF ( interiorOnly .AND. W_edge ) iMinUpd = 1
                0703           IF ( interiorOnly .AND. E_edge ) iMaxUpd = sNx
                0704           IF ( useGridArea ) THEN
                0705            DO j=1-OLy+1,sNy+OLy-1
                0706             DO i=iMinUpd,iMaxUpd
                0707              iceFld(i,j) = iceFld(i,j)
                0708      &                    -deltaTadvect*maskLocC(i,j)
                0709      &                                 *recip_rA(i,j,bi,bj)
                0710      &                                 *( af(i,j+1)-af(i,j) )
                0711             ENDDO
a4eca6e929 Jean*0712            ENDDO
204a3108f8 Jean*0713           ELSE
                0714            DO j=1-OLy+1,sNy+OLy-1
                0715             DO i=iMinUpd,iMaxUpd
                0716              tmpVol = iceVol(i,j)
                0717              iceVol(i,j) = iceVol(i,j)
                0718      &                    -deltaTadvect*maskLocC(i,j)
                0719      &                           *( vTrans(i,j+1)-vTrans(i,j) )
                0720              IF ( iceVol(i,j).GT.iceEps )
                0721      &       iceFld(i,j) = ( iceFld(i,j)*tmpVol
                0722      &                      -deltaTadvect*maskLocC(i,j)
                0723      &                                   *( af(i,j+1)-af(i,j) )
                0724      &                     )/iceVol(i,j)
                0725             ENDDO
a4eca6e929 Jean*0726            ENDDO
204a3108f8 Jean*0727           ENDIF
                0728 C-    keep advective flux (for diagnostics)
                0729           DO j=1-OLy+1,sNy+OLy
a4eca6e929 Jean*0730            DO i=iMinUpd,iMaxUpd
204a3108f8 Jean*0731              afy(i,j) = af(i,j)
a4eca6e929 Jean*0732            ENDDO
204a3108f8 Jean*0733           ENDDO
a4eca6e929 Jean*0734 
                0735 C      end if/else update overlap-Only
                0736         ENDIF
                0737 
                0738 C--   End of Y direction
                0739        ENDIF
                0740 
                0741 C--   End of ipass loop
                0742       ENDDO
                0743 
                0744 C-    explicit advection is done ; add some debugging print
204a3108f8 Jean*0745 #ifdef ALLOW_DBUG_THSICE
fef0c02f19 Jean*0746       IF ( dBugFlag ) THEN
a4eca6e929 Jean*0747        DO j=1-OLy,sNy+OLy
                0748         DO i=1-OLx,sNx+OLx
fef0c02f19 Jean*0749          IF ( i.EQ.idb .AND. j.EQ.jdb .AND. bi.EQ.biDb ) THEN
a4eca6e929 Jean*0750           tmpFac= deltaTadvect*recip_rA(i,j,bi,bj)
                0751           WRITE(6,'(A,1P4E14.6)') 'ICE_adv:',
                0752      &     afx(i,j)*tmpFac,afx(i+1,j)*tmpFac,
                0753      &     afy(i,j)*tmpFac,afy(i,j+1)*tmpFac
                0754          ENDIF
                0755         ENDDO
                0756        ENDDO
                0757       ENDIF
                0758 
                0759 #ifdef ALLOW_DEBUG
ae605e558b Jean*0760       IF ( debugLevel .GE. debLevC
a4eca6e929 Jean*0761      &     .AND. tracerIdentity.EQ.GAD_SI_HICE
                0762      &     .AND. k.LE.3 .AND. myIter.EQ.1+nIter0
                0763      &     .AND. nPx.EQ.1 .AND. nPy.EQ.1
                0764      &     .AND. useCubedSphereExchange ) THEN
                0765        CALL DEBUG_CS_CORNER_UV( ' afx,afy from THSICE_ADVECTION',
                0766      &      afx,afy, k, standardMessageUnit,bi,bj,myThid )
                0767       ENDIF
                0768 #endif /* ALLOW_DEBUG */
204a3108f8 Jean*0769 #endif /* ALLOW_DBUG_THSICE */
a4eca6e929 Jean*0770 
e7f517b398 Jean*0771 #endif /* ALLOW_GENERIC_ADVDIFF */
                0772 
a4eca6e929 Jean*0773       RETURN
                0774       END
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!