pkg/seaice/seaice_advdiff.F

03105a7583 Mart*0001 #include "SEAICE_OPTIONS.h"
772b2ed80e Gael*0002 #ifdef ALLOW_AUTODIFF
                0003 # include "AUTODIFF_OPTIONS.h"
                0004 #endif
f12f84b0ce Jean*0005
03105a7583 Mart*0006 CBOP
                0007 C !ROUTINE: SEAICE_ADVDIFF
                0008
                0009 C !INTERFACE: ==========================================================
f12f84b0ce Jean*0010       SUBROUTINE SEAICE_ADVDIFF(
f003822524 Jean*0011      U                  uc, vc,
                0012      I                  myTime, myIter, myThid )
03105a7583 Mart*0013
                0014 C !DESCRIPTION: \bv
f12f84b0ce Jean*0015 C     *===========================================================*
                0016 C     | SUBROUTINE SEAICE_ADVDIFF
                0017 C     | o driver for different advection routines
                0018 C     |   calls an adaption of gad_advection to call different
                0019 C     |   advection routines of pkg/generic_advdiff
                0020 C     *===========================================================*
                0021 C \ev
                0022
03105a7583 Mart*0023 C !USES: ===============================================================
f12f84b0ce Jean*0024       IMPLICIT NONE
                0025
03105a7583 Mart*0026 C     === Global variables ===
                0027 #include "SIZE.h"
                0028 #include "EEPARAMS.h"
                0029 #include "PARAMS.h"
                0030 #include "GRID.h"
ccaa3c61f4 Patr*0031 #include "SEAICE_SIZE.h"
03105a7583 Mart*0032 #include "SEAICE_PARAMS.h"
                0033 #include "SEAICE.h"
ccaa3c61f4 Patr*0034 #include "SEAICE_TRACER.h"
03105a7583 Mart*0035 #ifdef ALLOW_AUTODIFF_TAMC
                0036 # include "tamc.h"
                0037 #endif
                0038
f003822524 Jean*0039 C !INPUT/OUTPUT PARAMETERS: ===================================================
03105a7583 Mart*0040 C     === Routine arguments ===
f003822524 Jean*0041 C     uc/vc     :: current ice velocity on C-grid;
                0042 C               :: C-Grid : Input only ; B-grid : Output only
                0043 C     myTime    :: current time in simulation
                0044 C     myIter    :: current iteration number in simulation
                0045 C     myThid    :: my Thread Id number
                0046       _RL uc   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0047       _RL vc   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
03105a7583 Mart*0048       _RL myTime
                0049       INTEGER myIter
                0050       INTEGER myThid
                0051
                0052 C !LOCAL VARIABLES: ====================================================
                0053 C     === Local variables ===
f12f84b0ce Jean*0054 C     i,j,bi,bj :: Loop counters
1cf13b495f Torg*0055 C     it        :: Loop counter for ice thickness categories
f12f84b0ce Jean*0056 C     uTrans    :: volume transport, x direction
                0057 C     vTrans    :: volume transport, y direction
                0058 C     afx       :: horizontal advective flux, x direction
                0059 C     afy       :: horizontal advective flux, y direction
                0060 C     gFld      :: tendency of seaice field
                0061 C     xA,yA     :: "areas" of X and Y face of tracer cells
edfdf5fa1d Jean*0062       INTEGER i, j, bi, bj
86b84a92fc Patr*0063 #ifdef SEAICE_ITD
1cf13b495f Torg*0064       INTEGER it
b1ac83383d Mart*0065 #endif /* SEAICE_ITD */
d24acf3cfc Jean*0066 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0067 C     tkey :: tape key (depends on tiles)
                0068       INTEGER tkey
d24acf3cfc Jean*0069 #endif /* ALLOW_AUTODIFF_TAMC */
f50f58ec54 Gael*0070 #ifdef ALLOW_SITRACER
3dda23e098 Jean*0071       _RL hEffNm1    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0072       _RL areaNm1    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
ee0b114ca4 Jean*0073       INTEGER iTr, SEAICEadvSchSItr
                0074       _RL SEAICEdiffKhSItr
bb24b8a3e6 Gael*0075       _RL SItrExt   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
f50f58ec54 Gael*0076       _RL tmpscal1, tmpscal2
3dda23e098 Jean*0077 # ifdef ALLOW_SITRACER_ADVCAP
f50f58ec54 Gael*0078       _RL SItrPrev  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
3dda23e098 Jean*0079 # endif
                0080 # ifdef ALLOW_SITRACER_DEBUG_DIAG
f50f58ec54 Gael*0081       _RL DIAGarray (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
3dda23e098 Jean*0082 # endif
d24acf3cfc Jean*0083 #endif /* ALLOW_SITRACER */
3025a9909c Mart*0084       _RL fldNm1    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
f12f84b0ce Jean*0085       _RL uTrans    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0086       _RL vTrans    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0087       _RL afx       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0088       _RL afy       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
d24acf3cfc Jean*0089       _RL gFld      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
554418e294 Jean*0090       _RS xA        (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0091       _RS yA        (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
8852cab81a Mart*0092       _RL recip_heff(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
45315406aa Mart*0093 #if ( defined SEAICE_ITD || !defined ALLOW_GENERIC_ADVDIFF )
e0fa1cecbf Mart*0094       CHARACTER*(MAX_LEN_MBUF) msgBuf
45315406aa Mart*0095 #endif
f12f84b0ce Jean*0096 CEOP
                0097
                0098 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0099
45315406aa Mart*0100 C--   compute cell areas used by all tracers
47852c9c0c Mart*0101       DO bj=myByLo(myThid),myByHi(myThid)
                0102        DO bi=myBxLo(myThid),myBxHi(myThid)
45315406aa Mart*0103         DO j=1-OLy,sNy+OLy
                0104          DO i=1-OLx,sNx+OLx
                0105           xA(i,j,bi,bj) = _dyG(i,j,bi,bj)*SIMaskU(i,j,bi,bj)
                0106           yA(i,j,bi,bj) = _dxG(i,j,bi,bj)*SIMaskV(i,j,bi,bj)
                0107          ENDDO
                0108         ENDDO
                0109        ENDDO
                0110       ENDDO
                0111
                0112 #ifdef SEAICE_BGRID_DYNAMICS
f003822524 Jean*0113 C--   hack to ensure backward compatibility:
639dce29f3 Mart*0114 C     average B-grid seaice velocity to C-grid
45315406aa Mart*0115       DO bj=myByLo(myThid),myByHi(myThid)
                0116        DO bi=myBxLo(myThid),myBxHi(myThid)
d24acf3cfc Jean*0117         DO j=1-OLy,sNy+OLy-1
                0118          DO i=1-OLx,sNx+OLx-1
3025a9909c Mart*0119           uc(i,j,bi,bj)=.5 _d 0*(UICE(i,j,bi,bj)+UICE(i,j+1,bi,bj))
554418e294 Jean*0120           vc(i,j,bi,bj)=.5 _d 0*(VICE(i,j,bi,bj)+VICE(i+1,j,bi,bj))
47852c9c0c Mart*0121          ENDDO
                0122         ENDDO
                0123        ENDDO
                0124       ENDDO
                0125 C     Do we need this? I am afraid so.
                0126       CALL EXCH_UV_XY_RL(uc,vc,.TRUE.,myThid)
45315406aa Mart*0127 #endif
47852c9c0c Mart*0128
3dda23e098 Jean*0129 #ifdef ALLOW_AUTODIFF_TAMC
                0130 CADJ STORE uc   = comlev1, key = ikey_dynamics, kind=isbyte
                0131 CADJ STORE vc   = comlev1, key = ikey_dynamics, kind=isbyte
                0132 #endif /* ALLOW_AUTODIFF_TAMC */
                0133
03105a7583 Mart*0134       IF ( SEAICEmultiDimAdvection ) THEN
e0fa1cecbf Mart*0135 #ifdef ALLOW_GENERIC_ADVDIFF
1cb34809ba Patr*0136
03105a7583 Mart*0137       DO bj=myByLo(myThid),myByHi(myThid)
                0138        DO bi=myBxLo(myThid),myBxHi(myThid)
f12f84b0ce Jean*0139 C---   loops on tile indices bi,bj
                0140
8377b8ee87 Mart*0141 #ifdef ALLOW_AUTODIFF
                0142 C     Initialise to help AD-tools
d24acf3cfc Jean*0143         DO j=1-OLy,sNy+OLy
                0144          DO i=1-OLx,sNx+OLx
8377b8ee87 Mart*0145           gFld(i,j) = 0. _d 0
f12f84b0ce Jean*0146          ENDDO
                0147         ENDDO
8377b8ee87 Mart*0148 #endif
a68e054417 Patr*0149 C
8377b8ee87 Mart*0150 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0151         tkey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
                0152 CADJ STORE area(:,:,bi,bj)  = comlev1_bibj, key=tkey, kind=isbyte
                0153 CADJ STORE heff(:,:,bi,bj)  = comlev1_bibj, key=tkey, kind=isbyte
                0154 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
7246948615 Patr*0155 # ifdef SEAICE_VARIABLE_SALINITY
edb6656069 Mart*0156 CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
7246948615 Patr*0157 # endif
f12f84b0ce Jean*0158 #endif /* ALLOW_AUTODIFF_TAMC */
                0159
03105a7583 Mart*0160         DO j=1-OLy,sNy+OLy
                0161          DO i=1-OLx,sNx+OLx
050eb90cc6 Gael*0162 #ifdef ALLOW_SITRACER
3dda23e098 Jean*0163           hEffNm1(i,j,bi,bj) = HEFF(i,j,bi,bj)
                0164           areaNm1(i,j,bi,bj) = AREA(i,j,bi,bj)
                0165 #endif
04016f2c47 Mart*0166           recip_heff(i,j)    = 1. _d 0
03105a7583 Mart*0167          ENDDO
                0168         ENDDO
                0169
f12f84b0ce Jean*0170 C-    Calculate "volume transports" through tracer cell faces.
d24acf3cfc Jean*0171         DO j=1-OLy,sNy+OLy
                0172          DO i=1-OLx,sNx+OLx
554418e294 Jean*0173           uTrans(i,j) = uc(i,j,bi,bj)*xA(i,j,bi,bj)
                0174           vTrans(i,j) = vc(i,j,bi,bj)*yA(i,j,bi,bj)
03105a7583 Mart*0175          ENDDO
                0176         ENDDO
f12f84b0ce Jean*0177
e0fa1cecbf Mart*0178 #ifdef SEAICE_ITD
f48ea684a7 Jean*0179 C--   Effective Thickness (Volume)
09e3b53265 Mart*0180         IF ( SEAICEadvHeff ) THEN
e0fa1cecbf Mart*0181          DO it=1,SEAICE_multDim
                0182           CALL SEAICE_ADVECTION(
                0183      I         GAD_HEFF, SEAICEadvSchHeff,
                0184      I         uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj),
                0185      I         uTrans, vTrans, HEFFITD(1-OLx,1-OLy,it,bi,bj),
                0186      I         recip_heff,
                0187      O         gFld, afx, afy,
                0188      I         bi, bj, myTime, myIter, myThid )
                0189 C-    Add tendency due to diffusion
                0190           IF ( SEAICEdiffKhHeff .GT. 0. _d 0 )
                0191      &         CALL SEAICE_DIFFUSION(
                0192      I         GAD_HEFF, SEAICEdiffKhHeff, ONE,
                0193      I         HEFFITD(1-OLx,1-OLy,it,bi,bj), HEFFM,
                0194      I         xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
                0195      U         gFld,
                0196      I         bi, bj, myTime, myIter, myThid )
                0197 C     now do the "explicit" time step
                0198           DO j=1,sNy
                0199            DO i=1,sNx
                0200             HEFFITD(i,j,it,bi,bj) = HEFFM(i,j,bi,bj) * (
                0201      &           HEFFITD(i,j,it,bi,bj) + SEAICE_deltaTtherm * gFld(i,j)
                0202      &           )
                0203            ENDDO
86b84a92fc Patr*0204           ENDDO
                0205          ENDDO
e0fa1cecbf Mart*0206         ENDIF
                0207
                0208 C--   Fractional area
                0209         IF ( SEAICEadvArea ) THEN
                0210          DO it=1,SEAICE_multDim
                0211           CALL SEAICE_ADVECTION(
                0212      I         GAD_AREA, SEAICEadvSchArea,
                0213      I         uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj),
                0214      I         uTrans, vTrans, AREAITD(1-OLx,1-OLy,it,bi,bj),
                0215      I         recip_heff,
                0216      O         gFld, afx, afy,
                0217      I         bi, bj, myTime, myIter, myThid )
                0218 C-    Add tendency due to diffusion
                0219           IF ( SEAICEdiffKhArea .GT. 0. _d 0 )
                0220      &         CALL SEAICE_DIFFUSION(
                0221      I         GAD_AREA, SEAICEdiffKhArea, ONE,
                0222      I         AREAITD(1-OLx,1-OLy,it,bi,bj), HEFFM,
                0223      I         xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
                0224      U         gFld,
                0225      I         bi, bj, myTime, myIter, myThid )
                0226 C     now do the "explicit" time step
                0227           DO j=1,sNy
                0228            DO i=1,sNx
                0229             AREAITD(i,j,it,bi,bj) = HEFFM(i,j,bi,bj) * (
                0230      &           AREAITD(i,j,it,bi,bj) + SEAICE_deltaTtherm * gFld(i,j)
                0231      &           )
                0232            ENDDO
                0233           ENDDO
                0234          ENDDO
                0235 C     open water fraction needs to be advected for the ridging scheme
                0236          CALL SEAICE_ADVECTION(
                0237      I        GAD_AREA, SEAICEadvSchArea,
                0238      I        uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj),
                0239      I        uTrans, vTrans, opnWtrFrac(1-OLx,1-OLy,bi,bj), recip_heff,
                0240      O        gFld, afx, afy,
                0241      I        bi, bj, myTime, myIter, myThid )
                0242 C--   Add tendency due to diffusion
                0243          IF ( SEAICEdiffKhArea .GT. 0. _d 0 )
                0244      &        CALL SEAICE_DIFFUSION(
                0245      I        GAD_AREA, SEAICEdiffKhArea, ONE,
                0246      I        opnWtrFrac(1-OLx,1-OLy,bi,bj), HEFFM,
                0247      I        xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
                0248      U        gFld,
                0249      I        bi, bj, myTime, myIter, myThid )
                0250 C     now do the "explicit" time step
                0251          DO j=1,sNy
                0252           DO i=1,sNx
                0253            opnWtrFrac(i,j,bi,bj) = HEFFM(i,j,bi,bj) * (
                0254      &          opnWtrFrac(i,j,bi,bj) + SEAICE_deltaTtherm * gFld(i,j)
                0255      &          )
                0256           ENDDO
                0257          ENDDO
                0258         ENDIF
                0259
                0260 C--   Effective Snow Thickness (Volume)
                0261         IF ( SEAICEadvSnow ) THEN
                0262          DO it=1,SEAICE_multDim
                0263           CALL SEAICE_ADVECTION(
                0264      I         GAD_SNOW, SEAICEadvSchSnow,
                0265      I         uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj),
                0266      I         uTrans, vTrans, HSNOWITD(1-OLx,1-OLy,it,bi,bj),
                0267      I         recip_heff,
                0268      O         gFld, afx, afy,
                0269      I         bi, bj, myTime, myIter, myThid )
                0270 C--   Add tendency due to diffusion
                0271           IF ( SEAICEdiffKhSnow .GT. 0. _d 0 )
                0272      &         CALL SEAICE_DIFFUSION(
                0273      I         GAD_SNOW, SEAICEdiffKhSnow, ONE,
                0274      I         HSNOWITD(1-OLx,1-OLy,it,bi,bj), HEFFM,
                0275      I         xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
                0276      U         gFld,
                0277      I         bi, bj, myTime, myIter, myThid )
                0278 C     now do the "explicit" time step
                0279           DO j=1,sNy
                0280            DO i=1,sNx
                0281             HSNOWITD(i,j,it,bi,bj) = HEFFM(i,j,bi,bj) * (
                0282      &           HSNOWITD(i,j,it,bi,bj) + SEAICE_deltaTtherm*gFld(i,j)
                0283      &           )
                0284            ENDDO
                0285           ENDDO
                0286          ENDDO
                0287         ENDIF
                0288
                0289 C     update mean ice thickness HEFF and total ice concentration AREA
                0290 C     to match single category values
                0291 C     (necessary here because updated HEFF is used below for SItracer)
                0292         CALL SEAICE_ITD_SUM(bi, bj, myTime, myIter, myThid)
                0293
                0294 #else /* not SEAICE_ITD */
                0295 C--   Effective Thickness (Volume)
                0296         IF ( SEAICEadvHeff ) THEN
554418e294 Jean*0297           CALL SEAICE_ADVECTION(
                0298      I         GAD_HEFF, SEAICEadvSchHeff,
                0299      I         uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj),
                0300      I         uTrans, vTrans, HEFF(1-OLx,1-OLy,bi,bj), recip_heff,
                0301      O         gFld, afx, afy,
                0302      I         bi, bj, myTime, myIter, myThid )
6d78fc5463 Gael*0303          IF ( SEAICEdiffKhHeff .GT. 0. _d 0 ) THEN
f12f84b0ce Jean*0304 C-    Add tendency due to diffusion
                0305           CALL SEAICE_DIFFUSION(
b4c7dc19fe Jean*0306      I         GAD_HEFF, SEAICEdiffKhHeff, ONE,
554418e294 Jean*0307      I         HEFF(1-OLx,1-OLy,bi,bj), HEFFM,
                0308      I         xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
09e3b53265 Mart*0309      U         gFld,
                0310      I         bi, bj, myTime, myIter, myThid )
                0311          ENDIF
03105a7583 Mart*0312 C     now do the "explicit" time step
09e3b53265 Mart*0313          DO j=1,sNy
                0314           DO i=1,sNx
3025a9909c Mart*0315            HEFF(i,j,bi,bj) = HEFFM(i,j,bi,bj) * (
                0316      &          HEFF(i,j,bi,bj) + SEAICE_deltaTtherm * gFld(i,j)
09e3b53265 Mart*0317      &          )
                0318           ENDDO
03105a7583 Mart*0319          ENDDO
09e3b53265 Mart*0320         ENDIF
f12f84b0ce Jean*0321
                0322 C--   Fractional area
09e3b53265 Mart*0323         IF ( SEAICEadvArea ) THEN
554418e294 Jean*0324           CALL SEAICE_ADVECTION(
                0325      I         GAD_AREA, SEAICEadvSchArea,
                0326      I         uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj),
                0327      I         uTrans, vTrans, AREA(1-OLx,1-OLy,bi,bj), recip_heff,
                0328      O         gFld, afx, afy,
                0329      I         bi, bj, myTime, myIter, myThid )
6d78fc5463 Gael*0330          IF ( SEAICEdiffKhArea .GT. 0. _d 0 ) THEN
f12f84b0ce Jean*0331 C-    Add tendency due to diffusion
                0332           CALL SEAICE_DIFFUSION(
b4c7dc19fe Jean*0333      I         GAD_AREA, SEAICEdiffKhArea, ONE,
554418e294 Jean*0334      I         AREA(1-OLx,1-OLy,bi,bj), HEFFM,
                0335      I         xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
09e3b53265 Mart*0336      U         gFld,
                0337      I         bi, bj, myTime, myIter, myThid )
                0338          ENDIF
f12f84b0ce Jean*0339 C     now do the "explicit" time step
09e3b53265 Mart*0340          DO j=1,sNy
                0341           DO i=1,sNx
3025a9909c Mart*0342            AREA(i,j,bi,bj) = HEFFM(i,j,bi,bj) * (
                0343      &          AREA(i,j,bi,bj) + SEAICE_deltaTtherm * gFld(i,j)
09e3b53265 Mart*0344      &          )
                0345           ENDDO
f12f84b0ce Jean*0346          ENDDO
09e3b53265 Mart*0347         ENDIF
f12f84b0ce Jean*0348
eefc792916 Mart*0349 C--   Effective Snow Thickness (Volume)
09e3b53265 Mart*0350         IF ( SEAICEadvSnow ) THEN
554418e294 Jean*0351           CALL SEAICE_ADVECTION(
                0352      I         GAD_SNOW, SEAICEadvSchSnow,
                0353      I         uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj),
                0354      I         uTrans, vTrans, HSNOW(1-OLx,1-OLy,bi,bj), recip_heff,
                0355      O         gFld, afx, afy,
                0356      I         bi, bj, myTime, myIter, myThid )
6d78fc5463 Gael*0357          IF ( SEAICEdiffKhSnow .GT. 0. _d 0 ) THEN
09e3b53265 Mart*0358 C--   Add tendency due to diffusion
eefc792916 Mart*0359           CALL SEAICE_DIFFUSION(
b4c7dc19fe Jean*0360      I         GAD_SNOW, SEAICEdiffKhSnow, ONE,
554418e294 Jean*0361      I         HSNOW(1-OLx,1-OLy,bi,bj), HEFFM,
                0362      I         xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
09e3b53265 Mart*0363      U         gFld,
                0364      I         bi, bj, myTime, myIter, myThid )
                0365          ENDIF
eefc792916 Mart*0366 C     now do the "explicit" time step
09e3b53265 Mart*0367          DO j=1,sNy
                0368           DO i=1,sNx
                0369            HSNOW(i,j,bi,bj) = HEFFM(i,j,bi,bj) * (
                0370      &          HSNOW(i,j,bi,bj) + SEAICE_deltaTtherm * gFld(i,j)
                0371      &          )
                0372           ENDDO
eefc792916 Mart*0373          ENDDO
                0374         ENDIF
b1ac83383d Mart*0375 #endif /* SEAICE_ITD */
5df73465ef Torg*0376
a98c4b8072 Ian *0377 #ifdef SEAICE_VARIABLE_SALINITY
e06023d8a7 Dimi*0378 C--   Effective Sea Ice Salinity (Mass of salt)
fdfa8e151f Dimi*0379         IF ( SEAICEadvSalt ) THEN
554418e294 Jean*0380           CALL SEAICE_ADVECTION(
                0381      I         GAD_SALT, SEAICEadvSchSalt,
                0382      I         uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj),
                0383      I         uTrans, vTrans, HSALT(1-OLx,1-OLy,bi,bj), recip_heff,
                0384      O         gFld, afx, afy,
                0385      I         bi, bj, myTime, myIter, myThid )
6d78fc5463 Gael*0386          IF ( SEAICEdiffKhSalt .GT. 0. _d 0 ) THEN
fdfa8e151f Dimi*0387 C--   Add tendency due to diffusion
                0388           CALL SEAICE_DIFFUSION(
b4c7dc19fe Jean*0389      I         GAD_SALT, SEAICEdiffKhSalt, ONE,
554418e294 Jean*0390      I         HSALT(1-OLx,1-OLy,bi,bj), HEFFM,
                0391      I         xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
fdfa8e151f Dimi*0392      U         gFld,
                0393      I         bi, bj, myTime, myIter, myThid )
                0394          ENDIF
                0395 C     now do the "explicit" time step
                0396          DO j=1,sNy
                0397           DO i=1,sNx
                0398            HSALT(i,j,bi,bj) = HEFFM(i,j,bi,bj) * (
                0399      &          HSALT(i,j,bi,bj) + SEAICE_deltaTtherm * gFld(i,j)
                0400      &          )
                0401           ENDDO
                0402          ENDDO
                0403         ENDIF
a98c4b8072 Ian *0404 #endif /* SEAICE_VARIABLE_SALINITY */
fdfa8e151f Dimi*0405
f50f58ec54 Gael*0406 #ifdef ALLOW_SITRACER
                0407 C--   Sea Ice Tracers
be02c52974 Gael*0408         DO iTr = 1, SItrNumInUse
bb24b8a3e6 Gael*0409         IF ( (SEAICEadvHEFF.AND.(SItrMate(iTr).EQ.'HEFF')).OR.
                0410      &       (SEAICEadvAREA.AND.(SItrMate(iTr).EQ.'AREA')) ) THEN
                0411 C--   scale to effective value
                0412          IF (SItrMate(iTr).EQ.'HEFF') THEN
                0413           SEAICEadvSchSItr=SEAICEadvSchHEFF
                0414           SEAICEdiffKhSItr=SEAICEdiffKhHEFF
d24acf3cfc Jean*0415           DO j=1-OLy,sNy+OLy
                0416            DO i=1-OLx,sNx+OLx
bb24b8a3e6 Gael*0417             SItrExt(i,j,bi,bj) = HEFFM(i,j,bi,bj) *
3dda23e098 Jean*0418      &          SItracer(i,j,bi,bj,iTr) * hEffNm1(i,j,bi,bj)
bb24b8a3e6 Gael*0419            ENDDO
                0420           ENDDO
                0421 c TAF?   ELSEIF (SItrMate(iTr).EQ.'AREA') THEN
                0422          ELSE
                0423           SEAICEadvSchSItr=SEAICEadvSchAREA
                0424           SEAICEdiffKhSItr=SEAICEdiffKhAREA
d24acf3cfc Jean*0425           DO j=1-OLy,sNy+OLy
                0426            DO i=1-OLx,sNx+OLx
bb24b8a3e6 Gael*0427             SItrExt(i,j,bi,bj) = HEFFM(i,j,bi,bj) *
3dda23e098 Jean*0428      &          SItracer(i,j,bi,bj,iTr) * areaNm1(i,j,bi,bj)
bb24b8a3e6 Gael*0429            ENDDO
                0430           ENDDO
d24acf3cfc Jean*0431          ENDIF
bb24b8a3e6 Gael*0432 C--   store a couple things
d24acf3cfc Jean*0433           DO j=1-OLy,sNy+OLy
                0434            DO i=1-OLx,sNx+OLx
f50f58ec54 Gael*0435 #ifdef ALLOW_SITRACER_ADVCAP
                0436 C--   store previous value for spurious maxima treament
                0437             SItrPrev(i,j,bi,bj)=SItracer(i,j,bi,bj,iTr)
                0438 #endif
3721cfe5e4 Gael*0439 #ifdef ALLOW_SITRACER_DEBUG_DIAG
8377b8ee87 Mart*0440             diagArray(i,j,2+(iTr-1)*5) = SItrExt(i,j,bi,bj)
f50f58ec54 Gael*0441 #endif
                0442            ENDDO
                0443           ENDDO
d24acf3cfc Jean*0444 C--   compute advective tendency
f50f58ec54 Gael*0445           CALL SEAICE_ADVECTION(
bb24b8a3e6 Gael*0446      I         GAD_SITR+iTr-1, SEAICEadvSchSItr,
f50f58ec54 Gael*0447      I         uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj),
bb24b8a3e6 Gael*0448      I         uTrans, vTrans, SItrExt(1-OLx,1-OLy,bi,bj),
f50f58ec54 Gael*0449      I         recip_heff,
                0450      O         gFld, afx, afy,
                0451      I         bi, bj, myTime, myIter, myThid )
                0452           IF ( SEAICEdiffKhHeff .GT. 0. _d 0 ) THEN
                0453 C--   add diffusive tendency
                0454           CALL SEAICE_DIFFUSION(
bb24b8a3e6 Gael*0455      I         GAD_SITR+iTr-1, SEAICEdiffKhSItr, ONE,
                0456      I         SItrExt(1-OLx,1-OLy,bi,bj), HEFFM,
f50f58ec54 Gael*0457      I         xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
                0458      U         gFld,
                0459      I         bi, bj, myTime, myIter, myThid )
                0460           ENDIF
                0461 C--   apply tendency
                0462           DO j=1,sNy
                0463            DO i=1,sNx
bb24b8a3e6 Gael*0464             SItrExt(i,j,bi,bj) = HEFFM(i,j,bi,bj) * (
                0465      &        SItrExt(i,j,bi,bj) + SEAICE_deltaTtherm * gFld(i,j) )
f50f58ec54 Gael*0466            ENDDO
                0467           ENDDO
                0468 C--   scale back to actual value, or move effective value to ocean bucket
bb24b8a3e6 Gael*0469          IF (SItrMate(iTr).EQ.'HEFF') THEN
f50f58ec54 Gael*0470           DO j=1,sNy
                0471            DO i=1,sNx
8377b8ee87 Mart*0472             IF (HEFF(i,j,bi,bj).GE.siEps) THEN
                0473              SItracer(i,j,bi,bj,iTr)=SItrExt(i,j,bi,bj)/HEFF(i,j,bi,bj)
                0474              SItrBucket(i,j,bi,bj,iTr)=0. _d 0
                0475             ELSE
                0476              SItracer(i,j,bi,bj,iTr)=0. _d 0
                0477              SItrBucket(i,j,bi,bj,iTr)=SItrExt(i,j,bi,bj)
                0478             ENDIF
f50f58ec54 Gael*0479 #ifdef ALLOW_SITRACER_ADVCAP
3dda23e098 Jean*0480 C hack to try avoid 'spontaneous generation' of maxima, which supposedly would
                0481 C occur less frequently if we advected SItr with uXheff instead SItrXheff with u
8377b8ee87 Mart*0482             tmpscal1=max(SItrPrev(i,j,bi,bj),
f50f58ec54 Gael*0483      &       SItrPrev(i+1,j,bi,bj),SItrPrev(i-1,j,bi,bj),
                0484      &       SItrPrev(i,j+1,bi,bj),SItrPrev(i,j-1,bi,bj))
8377b8ee87 Mart*0485             tmpscal2=MAX(ZERO,SItracer(i,j,bi,bj,iTr)-tmpscal1)
                0486             SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)-tmpscal2
                0487             SItrBucket(i,j,bi,bj,iTr)=SItrBucket(i,j,bi,bj,iTr)
                0488      &           +tmpscal2*HEFF(i,j,bi,bj)
f50f58ec54 Gael*0489 #endif
3dda23e098 Jean*0490 C           treat case of potential negative value
8377b8ee87 Mart*0491             IF (HEFF(i,j,bi,bj).GE.siEps) THEN
                0492              tmpscal1=MIN(0. _d 0,SItracer(i,j,bi,bj,iTr))
                0493              SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)-tmpscal1
                0494              SItrBucket(i,j,bi,bj,iTr)=SItrBucket(i,j,bi,bj,iTr)
                0495      &                                 +HEFF(i,j,bi,bj)*tmpscal1
                0496             ENDIF
3721cfe5e4 Gael*0497 #ifdef ALLOW_SITRACER_DEBUG_DIAG
8377b8ee87 Mart*0498             diagArray(i,j,1+(iTr-1)*5)= -SItrBucket(i,j,bi,bj,iTr)
                0499      &      * HEFFM(i,j,bi,bj)/SEAICE_deltaTtherm*SEAICE_rhoIce
                0500             tmpscal1= ( HEFF(i,j,bi,bj)*SItracer(i,j,bi,bj,iTr)
                0501      &      + SItrBucket(i,j,bi,bj,iTr) )*HEFFM(i,j,bi,bj)
                0502             diagArray(i,j,2+(iTr-1)*5)= tmpscal1
                0503      &                                - diagArray(i,j,2+(iTr-1)*5)
                0504             diagArray(i,j,3+(iTr-1)*5)= HEFFM(i,j,bi,bj)
                0505      &      * SEAICE_deltaTtherm * gFld(i,j)
f50f58ec54 Gael*0506 #endif
                0507            ENDDO
                0508           ENDDO
bb24b8a3e6 Gael*0509 c TAF?   ELSEIF (SItrMate(iTr).EQ.'AREA') THEN
                0510          ELSE
                0511           DO j=1,sNy
                0512            DO i=1,sNx
8377b8ee87 Mart*0513             IF (AREA(i,j,bi,bj).GE.SEAICE_area_floor) THEN
                0514              SItracer(i,j,bi,bj,iTr)=SItrExt(i,j,bi,bj)/AREA(i,j,bi,bj)
                0515             ELSE
                0516              SItracer(i,j,bi,bj,iTr)=0. _d 0
                0517             ENDIF
bb24b8a3e6 Gael*0518             SItrBucket(i,j,bi,bj,iTr)=0. _d 0
                0519 #ifdef ALLOW_SITRACER_ADVCAP
8377b8ee87 Mart*0520             tmpscal1=max(SItrPrev(i,j,bi,bj),
bb24b8a3e6 Gael*0521      &       SItrPrev(i+1,j,bi,bj),SItrPrev(i-1,j,bi,bj),
                0522      &       SItrPrev(i,j+1,bi,bj),SItrPrev(i,j-1,bi,bj))
8377b8ee87 Mart*0523             tmpscal2=MAX(ZERO,SItracer(i,j,bi,bj,iTr)-tmpscal1)
                0524             SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)-tmpscal2
bb24b8a3e6 Gael*0525 #endif
3dda23e098 Jean*0526 C           treat case of potential negative value
8377b8ee87 Mart*0527             IF (AREA(i,j,bi,bj).GE.SEAICE_area_floor) THEN
bb24b8a3e6 Gael*0528               tmpscal1=MIN(0. _d 0,SItracer(i,j,bi,bj,iTr))
                0529               SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)-tmpscal1
8377b8ee87 Mart*0530             ENDIF
3721cfe5e4 Gael*0531 #ifdef ALLOW_SITRACER_DEBUG_DIAG
8377b8ee87 Mart*0532             diagArray(i,j,1+(iTr-1)*5)= 0. _d 0
                0533             diagArray(i,j,2+(iTr-1)*5)= - diagArray(i,j,2+(iTr-1)*5)
                0534      &      + AREA(i,j,bi,bj)*SItracer(i,j,bi,bj,iTr)*HEFFM(i,j,bi,bj)
                0535             diagArray(i,j,3+(iTr-1)*5)=HEFFM(i,j,bi,bj)
                0536      &      * SEAICE_deltaTtherm * gFld(i,j)
bb24b8a3e6 Gael*0537 #endif
                0538            ENDDO
                0539           ENDDO
                0540          ENDIF
f50f58ec54 Gael*0541 C--
bb24b8a3e6 Gael*0542          ENDIF
                0543         ENDDO
3721cfe5e4 Gael*0544 #ifdef ALLOW_SITRACER_DEBUG_DIAG
fcdcf3fa1d Jean*0545 c     CALL DIAGNOSTICS_FILL(DIAGarray,'UDIAG2  ',0,Nr,2,bi,bj,myThid)
f50f58ec54 Gael*0546 #endif
                0547 #endif /* ALLOW_SITRACER */
                0548
f12f84b0ce Jean*0549 C---   end bi,bj loops
03105a7583 Mart*0550        ENDDO
                0551       ENDDO
                0552
e0fa1cecbf Mart*0553 #else /* not ALLOW_GENERIC_ADVDIFF */
8377b8ee87 Mart*0554        WRITE(msgBuf,'(2A)')
e0fa1cecbf Mart*0555      &  'SEAICE_ADVDIFF: cannot use SEAICEmultiDimAdvection',
                0556      &  ' without pkg/generic_advdiff'
8377b8ee87 Mart*0557        CALL PRINT_ERROR( msgBuf , myThid )
                0558        WRITE(msgBuf,'(2A)') 'SEAICE_ADVDIFF: ',
e0fa1cecbf Mart*0559      &  'Re-compile with pkg "generic_advdiff" in packages.conf'
8377b8ee87 Mart*0560        CALL PRINT_ERROR( msgBuf , myThid )
                0561        CALL ALL_PROC_DIE( myThid )
                0562        STOP 'ABNORMAL END: S/R SEAICE_ADVDIFF'
e0fa1cecbf Mart*0563 #endif /* ALLOW_GENERIC_ADVDIFF */
03105a7583 Mart*0564       ELSE
47852c9c0c Mart*0565 C--   if not multiDimAdvection
e0fa1cecbf Mart*0566 #ifdef SEAICE_ITD
                0567 C     just for safety
8377b8ee87 Mart*0568        WRITE(msgBuf,'(2A)') 'SEAICE_ADVDIFF: ',
e0fa1cecbf Mart*0569      &  'ITD with SEAICEmultiDimAdvection=.False. is not allowed,'
8377b8ee87 Mart*0570        CALL PRINT_ERROR( msgBuf , myThid )
                0571        WRITE(msgBuf,'(2A)') 'SEAICE_ADVDIFF: ',
e0fa1cecbf Mart*0572      &  'use a multidimensional advection scheme (in data.seaice)'
8377b8ee87 Mart*0573        CALL PRINT_ERROR( msgBuf , myThid )
                0574        CALL ALL_PROC_DIE( myThid )
                0575        STOP 'ABNORMAL END: S/R SEAICE_ADVDIFF'
e0fa1cecbf Mart*0576 #endif /* SEAICE_ITD */
1cb34809ba Patr*0577
04016f2c47 Mart*0578        IF ( SEAICEadvHEff ) THEN
a68e054417 Patr*0579 #ifdef ALLOW_AUTODIFF_TAMC
                0580 CADJ STORE heff = comlev1, key = ikey_dynamics, kind=isbyte
                0581 #endif
3dda23e098 Jean*0582         CALL ADVECT( uc, vc, hEff, fldNm1, HEFFM, myThid )
6d78fc5463 Gael*0583         IF ( SEAICEdiffKhHeff .GT. 0. _d 0 ) THEN
                0584 C-    Add tendency due to diffusion
554418e294 Jean*0585          DO bj=myByLo(myThid),myByHi(myThid)
                0586           DO bi=myBxLo(myThid),myBxHi(myThid)
                0587             CALL SEAICE_DIFFUSION(
                0588      I           GAD_HEFF, SEAICEdiffKhHeff, SEAICE_deltaTtherm,
3dda23e098 Jean*0589      I           fldNm1(1-OLx,1-OLy,bi,bj), HEFFM,
554418e294 Jean*0590      I           xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
                0591      U           HEFF(1-OLx,1-OLy,bi,bj),
                0592      I           bi, bj, myTime, myIter, myThid )
                0593           ENDDO
                0594          ENDDO
6d78fc5463 Gael*0595         ENDIF
eefc792916 Mart*0596        ENDIF
04016f2c47 Mart*0597        IF ( SEAICEadvArea ) THEN
a68e054417 Patr*0598 #ifdef ALLOW_AUTODIFF_TAMC
                0599 CADJ STORE area = comlev1, key = ikey_dynamics, kind=isbyte
                0600 #endif
3dda23e098 Jean*0601         CALL ADVECT( uc, vc, area, fldNm1, HEFFM, myThid )
6d78fc5463 Gael*0602         IF ( SEAICEdiffKhArea .GT. 0. _d 0 ) THEN
                0603 C-    Add tendency due to diffusion
554418e294 Jean*0604          DO bj=myByLo(myThid),myByHi(myThid)
                0605           DO bi=myBxLo(myThid),myBxHi(myThid)
                0606             CALL SEAICE_DIFFUSION(
                0607      I           GAD_AREA, SEAICEdiffKhArea, SEAICE_deltaTtherm,
3dda23e098 Jean*0608      I           fldNm1(1-OLx,1-OLy,bi,bj), HEFFM,
554418e294 Jean*0609      I           xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
                0610      U           Area(1-OLx,1-OLy,bi,bj),
                0611      I           bi, bj, myTime, myIter, myThid )
                0612           ENDDO
                0613          ENDDO
6d78fc5463 Gael*0614         ENDIF
fdfa8e151f Dimi*0615        ENDIF
04016f2c47 Mart*0616        IF ( SEAICEadvSnow ) THEN
a68e054417 Patr*0617 #ifdef ALLOW_AUTODIFF_TAMC
                0618 CADJ STORE hsnow = comlev1, key = ikey_dynamics, kind=isbyte
                0619 #endif
3025a9909c Mart*0620         CALL ADVECT( uc, vc, HSNOW, fldNm1, HEFFM, myThid )
6d78fc5463 Gael*0621         IF ( SEAICEdiffKhSnow .GT. 0. _d 0 ) THEN
                0622 C-    Add tendency due to diffusion
554418e294 Jean*0623          DO bj=myByLo(myThid),myByHi(myThid)
                0624           DO bi=myBxLo(myThid),myBxHi(myThid)
                0625             CALL SEAICE_DIFFUSION(
                0626      I           GAD_SNOW, SEAICEdiffKhSnow, SEAICE_deltaTtherm,
                0627      I           fldNm1(1-OLx,1-OLy,bi,bj), HEFFM,
                0628      I           xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
                0629      U           HSNOW(1-OLx,1-OLy,bi,bj),
                0630      I           bi, bj, myTime, myIter, myThid )
                0631           ENDDO
                0632          ENDDO
6d78fc5463 Gael*0633         ENDIF
04016f2c47 Mart*0634        ENDIF
                0635
a98c4b8072 Ian *0636 #ifdef SEAICE_VARIABLE_SALINITY
04016f2c47 Mart*0637        IF ( SEAICEadvSalt ) THEN
a68e054417 Patr*0638 #ifdef ALLOW_AUTODIFF_TAMC
                0639 CADJ STORE hsalt = comlev1, key = ikey_dynamics, kind=isbyte
                0640 #endif
3025a9909c Mart*0641         CALL ADVECT( uc, vc, HSALT, fldNm1, HEFFM, myThid )
6d78fc5463 Gael*0642         IF ( SEAICEdiffKhSalt .GT. 0. _d 0 ) THEN
                0643 C-    Add tendency due to diffusion
554418e294 Jean*0644          DO bj=myByLo(myThid),myByHi(myThid)
                0645           DO bi=myBxLo(myThid),myBxHi(myThid)
                0646             CALL SEAICE_DIFFUSION(
                0647      I           GAD_SALT, SEAICEdiffKhSalt, SEAICE_deltaTtherm,
                0648      I           fldNm1(1-OLx,1-OLy,bi,bj), HEFFM,
                0649      I           xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
                0650      U           HSALT(1-OLx,1-OLy,bi,bj),
                0651      I           bi, bj, myTime, myIter, myThid )
                0652           ENDDO
                0653          ENDDO
6d78fc5463 Gael*0654         ENDIF
04016f2c47 Mart*0655        ENDIF
a98c4b8072 Ian *0656 #endif /* SEAICE_VARIABLE_SALINITY */
04016f2c47 Mart*0657
47852c9c0c Mart*0658 C--   end if multiDimAdvection
03105a7583 Mart*0659       ENDIF
0f5999dbd4 Mart*0660
03105a7583 Mart*0661       RETURN
                0662       END