pkg/seaice/seaice_model.F

809c36b928 Patr*0001 #include "SEAICE_OPTIONS.h"
8377b8ee87 Mart*0002 #ifdef ALLOW_EXF
                0003 # include "EXF_OPTIONS.h"
                0004 #endif
772b2ed80e Gael*0005 #ifdef ALLOW_AUTODIFF
                0006 # include "AUTODIFF_OPTIONS.h"
                0007 #endif
7e1af563d8 Jean*0008
ef080e1d37 Dimi*0009 CBOP
                0010 C !ROUTINE: SEAICE_MODEL
                0011
                0012 C !INTERFACE: ==========================================================
b92ebad606 Jean*0013       SUBROUTINE SEAICE_MODEL( myTime, myIter, myThid )
ef080e1d37 Dimi*0014
                0015 C !DESCRIPTION: \bv
ff8e41478f Jean*0016 C     *===========================================================*
09510da3bb Dimi*0017 C     | SUBROUTINE SEAICE_MODEL                                   |
                0018 C     | o Time stepping of a dynamic/thermodynamic sea ice model. |
                0019 C     |  Dynamics solver: Zhang/Hibler, JGR, 102, 8691-8702, 1997 |
                0020 C     |  Thermodynamics:        Hibler, MWR, 108, 1943-1973, 1980 |
                0021 C     |  Rheology:              Hibler, JPO,   9,  815- 846, 1979 |
                0022 C     |  Snow:          Zhang et al.  , JPO,  28,  191- 217, 1998 |
                0023 C     |  Parallel forward ice model written by Jinlun Zhang PSC/UW|
                0024 C     |  & coupled into MITgcm by Dimitris Menemenlis (JPL) 2/2001|
                0025 C     |  zhang@apl.washington.edu / menemenlis@jpl.nasa.gov       |
45315406aa Mart*0026 C     | o The code has been rewritten substantially to use the    |
                0027 C     |   MITgcm C-grid, see Losch et al. OM, 33,  129- 144, 2010 |
ff8e41478f Jean*0028 C     *===========================================================*
809c36b928 Patr*0029       IMPLICIT NONE
7e1af563d8 Jean*0030 C \ev
                0031
ef080e1d37 Dimi*0032 C !USES: ===============================================================
809c36b928 Patr*0033 #include "SIZE.h"
                0034 #include "EEPARAMS.h"
                0035 #include "DYNVARS.h"
                0036 #include "PARAMS.h"
4634ed4578 Dimi*0037 #include "GRID.h"
809c36b928 Patr*0038 #include "FFIELDS.h"
ccaa3c61f4 Patr*0039 #include "SEAICE_SIZE.h"
809c36b928 Patr*0040 #include "SEAICE_PARAMS.h"
ccaa3c61f4 Patr*0041 #include "SEAICE.h"
                0042 #include "SEAICE_TRACER.h"
ae1fb66b64 Dimi*0043 #ifdef ALLOW_EXF
                0044 # include "EXF_FIELDS.h"
                0045 #endif
7109a141b2 Patr*0046 #ifdef ALLOW_AUTODIFF_TAMC
                0047 # include "tamc.h"
6060ec2938 Dimi*0048 #endif
baa476eeba Dimi*0049
ef080e1d37 Dimi*0050 C !INPUT PARAMETERS: ===================================================
f003822524 Jean*0051 C     myTime :: Current time in simulation
                0052 C     myIter :: Current iteration number in simulation
                0053 C     myThid :: my Thread Id number
809c36b928 Patr*0054       _RL     myTime
                0055       INTEGER myIter
                0056       INTEGER myThid
                0057
ef080e1d37 Dimi*0058 C !LOCAL VARIABLES: ====================================================
a1ab12d5e7 Dimi*0059 C     i,j,bi,bj :: Loop counters
e8a5a3baa4 Jean*0060       INTEGER i, j
                0061       INTEGER bi, bj
40d541aac0 Jean*0062 #ifdef ALLOW_EXF
                0063       INTEGER grpDiag
                0064 #endif
f50f58ec54 Gael*0065 #ifdef ALLOW_SITRACER
                0066       INTEGER iTr
                0067 #endif
45315406aa Mart*0068 #ifdef SEAICE_BGRID_DYNAMICS
f003822524 Jean*0069       _RL uLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0070       _RL vLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0071 #endif
ef080e1d37 Dimi*0072 CEOP
5fe78992ba Mart*0073 #ifdef ALLOW_DIAGNOSTICS
                0074       LOGICAL  DIAGNOSTICS_IS_ON
                0075       EXTERNAL DIAGNOSTICS_IS_ON
                0076       LOGICAL  diag_SIenph_isOn
ee480e7782 amc2*0077       _RL mTmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
5fe78992ba Mart*0078       _RL recip_deltaTtherm
ee480e7782 amc2*0079       _RL addSnow
5fe78992ba Mart*0080 #endif
b92ebad606 Jean*0081 #ifdef ALLOW_DEBUG
a8b75793c4 Jean*0082       IF (debugMode) CALL DEBUG_ENTER( 'SEAICE_MODEL', myThid )
b92ebad606 Jean*0083 #endif
                0084
0320e25227 Mart*0085       IF ( usingPCoords ) THEN
                0086 C     In z-coordinates, phiHydLow is just a diagnostics and hence its
                0087 C     overlaps are never filled properly
                0088 C     in p-coordinates, phiHydLow is the sea level elevation that is required
                0089 C     for the sea level tilt term in the momentum equation, so we need
                0090 C     to fill the overlaps properly here
                0091 C     (or elsewhere, not sure where it would be best)
                0092        _EXCH_XY_RL( phiHydLow, myThid )
                0093       ENDIF
47852c9c0c Mart*0094 #ifdef ALLOW_THSICE
                0095       IF ( useThSice ) THEN
20e951135e Mart*0096 C--   Map thSice-variables to HEFF, HSNOW, and AREA, because they are
                0097 C     needed in S/R SEAICE_DYNSOLVER
47852c9c0c Mart*0098        CALL SEAICE_MAP_THSICE( myTime, myIter, myThid )
                0099       ENDIF
                0100 #endif /* ALLOW_THSICE */
d5254b4e3d Mart*0101 #ifdef ALLOW_DIAGNOSTICS
                0102       IF ( useDiagnostics ) THEN
                0103        diag_SIenph_isOn = DIAGNOSTICS_IS_ON('SIenph  ',myThid)
                0104        IF ( diag_SIenph_isOn ) THEN
ee480e7782 amc2*0105         IF ( SEAICEaddSnowMass ) THEN
                0106          addSnow=1. _d 0
                0107         ELSE
                0108          addSnow=0. _d 0
                0109         ENDIF
d5254b4e3d Mart*0110         DO bj=myByLo(myThid),myByHi(myThid)
                0111          DO bi=myBxLo(myThid),myBxHi(myThid)
                0112           DO j=1-OLy,sNy+OLy
                0113            DO i=1-OLx,sNx+OLx
ee480e7782 amc2*0114             mTmp(i,j,bi,bj)=SEAICE_rhoIce*HEFF(i,j,bi,bj)
                0115      &                     +addSnow*SEAICE_rhoSnow*HSNOW(i,j,bi,bj)
d5254b4e3d Mart*0116            ENDDO
                0117           ENDDO
                0118          ENDDO
                0119         ENDDO
                0120        ENDIF
                0121       ENDIF
                0122 #endif /* ALLOW_DIAGNOSTICS */
47852c9c0c Mart*0123
40d541aac0 Jean*0124 #ifdef ALLOW_EXF
                0125       IF ( useEXF ) THEN
                0126 C--   Winds are from pkg/exf, which does not update edges.
                0127        CALL EXCH_UV_AGRID_3D_RL( uwind, vwind, .TRUE., 1, myThid )
                0128        IF ( useDiagnostics ) THEN
                0129 C-    Fill-in EXF wind-stess diags, weighted by open-ocean fraction
                0130         grpDiag = -1
                0131         IF ( SEAICEuseDYNAMICS ) grpDiag = 1
                0132         CALL EXF_WEIGHT_SFX_DIAGS(
                0133      I                  AREA, grpDiag, myTime, myIter, myThid )
                0134        ENDIF
                0135       ENDIF
                0136 #endif /* ALLOW_EXF */
                0137
8377b8ee87 Mart*0138 #ifdef ALLOW_AUTODIFF
c20cddf271 Patr*0139       DO bj=myByLo(myThid),myByHi(myThid)
                0140        DO bi=myBxLo(myThid),myBxHi(myThid)
d24acf3cfc Jean*0141         DO j=1-OLy,sNy+OLy
                0142          DO i=1-OLx,sNx+OLx
c20cddf271 Patr*0143           uIceNm1(i,j,bi,bj) = 0. _d 0
                0144           vIceNm1(i,j,bi,bj) = 0. _d 0
aef1621d33 Patr*0145 # ifdef ALLOW_SITRACER
                0146           DO iTr = 1, SItrMaxNum
                0147            SItrBucket(i,j,bi,bj,iTr) = 0. _d 0
                0148           ENDDO
                0149 # endif
c20cddf271 Patr*0150          ENDDO
                0151         ENDDO
                0152        ENDDO
                0153       ENDDO
8377b8ee87 Mart*0154 #endif /* ALLOW_AUTODIFF */
                0155 #ifdef ALLOW_AUTODIFF_TAMC
95c72ef3a1 Patr*0156 CADJ STORE heff  = comlev1, key=ikey_dynamics, kind=isbyte
                0157 CADJ STORE area  = comlev1, key=ikey_dynamics, kind=isbyte
3cb8290a3c Mart*0158 CADJ STORE hsnow = comlev1, key=ikey_dynamics, kind=isbyte
45315406aa Mart*0159 # ifdef SEAICE_CGRID
95c72ef3a1 Patr*0160 CADJ STORE fu    = comlev1, key=ikey_dynamics, kind=isbyte
                0161 CADJ STORE fv    = comlev1, key=ikey_dynamics, kind=isbyte
d5254b4e3d Mart*0162 #  ifdef SEAICE_ALLOW_EVP
b81c343615 Patr*0163 CADJ STORE seaice_sigma1  = comlev1, key=ikey_dynamics, kind=isbyte
                0164 CADJ STORE seaice_sigma2  = comlev1, key=ikey_dynamics, kind=isbyte
                0165 CADJ STORE seaice_sigma12 = comlev1, key=ikey_dynamics, kind=isbyte
d5254b4e3d Mart*0166 #  endif /* SEAICE_ALLOW_EVP */
45315406aa Mart*0167 # endif /* SEAICE_CGRID */
aef1621d33 Patr*0168 # ifdef ALLOW_SITRACER
d5254b4e3d Mart*0169 CADJ STORE siceload = comlev1, key=ikey_dynamics, kind=isbyte
                0170 CADJ STORE sitracer = comlev1, key=ikey_dynamics, kind=isbyte
aef1621d33 Patr*0171 # endif
7109a141b2 Patr*0172 #endif /* ALLOW_AUTODIFF_TAMC */
                0173
45315406aa Mart*0174 C--   Solve ice momentum equations and calculate ocean surface stress.
                0175 C     The surface stress always needs to be updated, even if neither B-
                0176 C     or C-grid dynamics are compiled, and SEAICEuseDYNAMICS = F.
b92ebad606 Jean*0177 #ifdef ALLOW_DEBUG
a8b75793c4 Jean*0178       IF (debugMode) CALL DEBUG_CALL( 'SEAICE_DYNSOLVER', myThid )
b92ebad606 Jean*0179 #endif
45315406aa Mart*0180 #ifdef SEAICE_BGRID_DYNAMICS
02aabb8fe5 Dimi*0181       CALL TIMER_START('DYNSOLVER          [SEAICE_MODEL]',myThid)
                0182       CALL DYNSOLVER ( myTime, myIter, myThid )
                0183       CALL TIMER_STOP ('DYNSOLVER          [SEAICE_MODEL]',myThid)
45315406aa Mart*0184 #else /* use default C-grid solver */
                0185       CALL TIMER_START('SEAICE_DYNSOLVER   [SEAICE_MODEL]',myThid)
                0186       CALL SEAICE_DYNSOLVER ( myTime, myIter, myThid )
                0187       CALL TIMER_STOP ('SEAICE_DYNSOLVER   [SEAICE_MODEL]',myThid)
                0188 #endif
809c36b928 Patr*0189
0f64a7d6fb Dimi*0190 C--   Apply ice velocity open boundary conditions
5459643feb Dimi*0191 #ifdef ALLOW_OBCS
0f64a7d6fb Dimi*0192 # ifndef DISABLE_SEAICE_OBCS
d5254b4e3d Mart*0193       IF ( useOBCS ) CALL OBCS_ADJUST_UVICE( uice, vice, myThid )
0f64a7d6fb Dimi*0194 # endif /* DISABLE_SEAICE_OBCS */
5459643feb Dimi*0195 #endif /* ALLOW_OBCS */
                0196
45315406aa Mart*0197 #if ( defined ALLOW_AUTODIFF_TAMC &&  defined SEAICE_CGRID )
d5254b4e3d Mart*0198 CADJ STORE uice  = comlev1, key=ikey_dynamics, kind=isbyte
                0199 CADJ STORE vice  = comlev1, key=ikey_dynamics, kind=isbyte
                0200 C     Note: Storing u/vice **after** seaice_dynsolver (and obcs_adjust_uvice)
                0201 C     has the effect that seaice_dynsolver is not called from seaice_model_ad
                0202 C     anymore. This is important because with the numerous tricks to avoid
                0203 C     complicated code in the backward integration (see pkg/autodiff),
                0204 C     the extra call would update the ice velocities with the wrong
                0205 C     set of flag values.
45315406aa Mart*0206 #endif
d5254b4e3d Mart*0207
4a98994c14 Jean*0208 #ifdef ALLOW_THSICE
de80e32bba Jean*0209       IF ( useThSice ) THEN
                0210 #ifdef ALLOW_DEBUG
d5254b4e3d Mart*0211        IF (debugMode) CALL DEBUG_CALL( 'THSICE_DO_ADVECT', myThid )
de80e32bba Jean*0212 #endif
d5254b4e3d Mart*0213        CALL THSICE_DO_ADVECT( 0, 0, myTime, myIter, myThid )
de80e32bba Jean*0214       ELSE
bfae2ee370 jm-c 0215 #endif /* ALLOW_THSICE */
ff8e41478f Jean*0216 C--   Only call advection of heff, area, snow, and salt and
09e3b53265 Mart*0217 C--   growth for the generic 0-layer thermodynamics of seaice
                0218 C--   if useThSice=.false., otherwise the 3-layer Winton thermodynamics
                0219 C--   (called from DO_OCEANIC_PHYSICS) take care of this
                0220 C NOW DO ADVECTION and DIFFUSION
d5254b4e3d Mart*0221        IF ( SEAICEadvHeff .OR. SEAICEadvArea .OR. SEAICEadvSnow
e54fe3e1f9 Gael*0222      &        .OR. SEAICEadvSalt ) THEN
b92ebad606 Jean*0223 #ifdef ALLOW_DEBUG
d5254b4e3d Mart*0224         IF (debugMode) CALL DEBUG_CALL( 'SEAICE_ADVDIFF', myThid )
b92ebad606 Jean*0225 #endif
45315406aa Mart*0226 C--   There always needs to be advection, even if neither B- or C-grid
                0227 C     dynamics are compiled.
                0228 #ifdef SEAICE_BGRID_DYNAMICS
d5254b4e3d Mart*0229         CALL SEAICE_ADVDIFF( uLoc, vLoc, myTime, myIter, myThid )
45315406aa Mart*0230 #else /* default C-grid advection */
                0231         CALL SEAICE_ADVDIFF( uIce, vIce, myTime, myIter, myThid )
                0232 #endif
d5254b4e3d Mart*0233        ENDIF
                0234 #ifdef ALLOW_AUTODIFF_TAMC
                0235 CADJ STORE heff  = comlev1, key=ikey_dynamics, kind=isbyte
                0236 CADJ STORE area  = comlev1, key=ikey_dynamics, kind=isbyte
                0237 CADJ STORE hsnow = comlev1, key=ikey_dynamics, kind=isbyte
                0238 # ifdef SEAICE_VARIABLE_SALINITY
                0239 CADJ STORE hsalt = comlev1, key=ikey_dynamics, kind=isbyte
                0240 # endif
                0241 C     Note: This store has the effect that seaice_advdiff is not called
                0242 C     from seaice_model_ad anymore. Instead, the stored values are used.
                0243 #endif /* ALLOW_AUTODIFF_TAMC */
809c36b928 Patr*0244
1cf549c217 Mart*0245 C     After advection, the sea ice variables may have unphysical values
20970ee22f Jean*0246 C     e.g., < 0, that are regularized here. Concentration as a special case
1cf549c217 Mart*0247 C     may be > 1 in convergent motion and a ridging algorithm redistributes
                0248 C     the ice to limit the concentration to 1.
d5254b4e3d Mart*0249        CALL SEAICE_REG_RIDGE( myTime, myIter, myThid )
1cf549c217 Mart*0250
40d541aac0 Jean*0251 #ifdef ALLOW_EXF
d5254b4e3d Mart*0252        IF ( useEXF .AND. useDiagnostics ) THEN
40d541aac0 Jean*0253 C-    Fill-in EXF surface flux diags, weighted by open-ocean fraction
                0254         grpDiag = -2
                0255         IF ( usePW79thermodynamics ) grpDiag = 2
                0256         CALL EXF_WEIGHT_SFX_DIAGS(
                0257      I                  AREA, grpDiag, myTime, myIter, myThid )
d5254b4e3d Mart*0258        ENDIF
40d541aac0 Jean*0259 #endif /* ALLOW_EXF */
                0260
bfae2ee370 jm-c 0261 #ifdef DISABLE_SEAICE_GROWTH
                0262        IF ( .TRUE. ) THEN
                0263 #else /* DISABLE_SEAICE_GROWTH */
09e3b53265 Mart*0264 C     thermodynamics growth
ff8e41478f Jean*0265 C     must call growth after calling advection
09e3b53265 Mart*0266 C     because of ugly time level business
                0267        IF ( usePW79thermodynamics ) THEN
45315406aa Mart*0268 # ifdef SEAICE_USE_GROWTH_ADX
                0269 #  ifdef ALLOW_DEBUG
dc54d31829 Ian *0270         IF (debugMode) CALL DEBUG_CALL( 'SEAICE_GROWTH_ADX', myThid )
45315406aa Mart*0271 #  endif
dc54d31829 Ian *0272         CALL SEAICE_GROWTH_ADX( myTime, myIter, myThid )
382462ccb5 Mart*0273 # else /* SEAICE_USE_GROWTH_ADX */
45315406aa Mart*0274 #  ifdef ALLOW_DEBUG
a8b75793c4 Jean*0275         IF (debugMode) CALL DEBUG_CALL( 'SEAICE_GROWTH', myThid )
45315406aa Mart*0276 #  endif
09e3b53265 Mart*0277         CALL SEAICE_GROWTH( myTime, myIter, myThid )
382462ccb5 Mart*0278 # endif /* SEAICE_USE_GROWTH_ADX */
bfae2ee370 jm-c 0279        ELSE
4ce2104710 Patr*0280 #endif /* DISABLE_SEAICE_GROWTH */
bfae2ee370 jm-c 0281         DO bj=myByLo(myThid),myByHi(myThid)
                0282          DO bi=myBxLo(myThid),myBxHi(myThid)
                0283           DO j=1,sNy
                0284            DO i=1,sNx
                0285             sIceLoad(i,j,bi,bj) = HEFF(i,j,bi,bj)*SEAICE_rhoIce
                0286      &                         + HSNOW(i,j,bi,bj)*SEAICE_rhoSnow
                0287            ENDDO
                0288           ENDDO
                0289 c#ifdef SEAICE_CAP_ICELOAD
                0290 c         sIceTooHeavy = rhoConst*drF(1) / 5. _d 0
                0291 c         DO j=1,sNy
                0292 c          DO i=1,sNx
                0293 c           sIceLoad(i,j,bi,bj) = MIN( sIceLoad(i,j,bi,bj),
                0294 c    &                                 sIceTooHeavy )
                0295 c          ENDDO
                0296 c         ENDDO
                0297 c#endif
                0298          ENDDO
                0299         ENDDO
                0300        ENDIF
809c36b928 Patr*0301
f50f58ec54 Gael*0302 #ifdef ALLOW_SITRACER
aef1621d33 Patr*0303 # ifdef ALLOW_AUTODIFF_TAMC
                0304 CADJ STORE sitracer  = comlev1, key=ikey_dynamics, kind=isbyte
                0305 # endif
a746f35887 Jean*0306        CALL SEAICE_TRACER_PHYS ( myTime, myIter, myThid )
f50f58ec54 Gael*0307 #endif
                0308
0f64a7d6fb Dimi*0309 C--   Apply ice tracer open boundary conditions
a1ab12d5e7 Dimi*0310 #ifdef ALLOW_OBCS
0f64a7d6fb Dimi*0311 # ifndef DISABLE_SEAICE_OBCS
                0312        IF ( useOBCS ) CALL OBCS_APPLY_SEAICE( myThid )
                0313 # endif /* DISABLE_SEAICE_OBCS */
a1ab12d5e7 Dimi*0314 #endif /* ALLOW_OBCS */
                0315
809c36b928 Patr*0316 C--   Update overlap regions for a bunch of stuff
772590b63c Mart*0317        _EXCH_XY_RL( HEFF,  myThid )
                0318        _EXCH_XY_RL( AREA,  myThid )
                0319        _EXCH_XY_RL( HSNOW, myThid )
963c7c7f95 Torg*0320 #ifdef SEAICE_ITD
                0321        CALL EXCH_3D_RL( HEFFITD,  nITD, myThid )
                0322        CALL EXCH_3D_RL( AREAITD,  nITD, myThid )
                0323        CALL EXCH_3D_RL( HSNOWITD, nITD, myThid )
                0324 #endif
a98c4b8072 Ian *0325 #ifdef SEAICE_VARIABLE_SALINITY
772590b63c Mart*0326        _EXCH_XY_RL( HSALT, myThid )
fdfa8e151f Dimi*0327 #endif
f50f58ec54 Gael*0328 #ifdef ALLOW_SITRACER
be02c52974 Gael*0329        DO iTr = 1, SItrNumInUse
d24acf3cfc Jean*0330         _EXCH_XY_RL( SItracer(1-OLx,1-OLy,1,1,iTr),myThid )
f50f58ec54 Gael*0331        ENDDO
                0332 #endif
14cee76775 Jean*0333        _EXCH_XY_RS(EmPmR, myThid )
                0334        _EXCH_XY_RS(saltFlux, myThid )
                0335        _EXCH_XY_RS(Qnet , myThid )
809c36b928 Patr*0336 #ifdef SHORTWAVE_HEATING
14cee76775 Jean*0337        _EXCH_XY_RS(Qsw  , myThid )
e4775240e5 Dimi*0338 #endif /* SHORTWAVE_HEATING */
f3ce416a61 Jean*0339 #ifdef ATMOSPHERIC_LOADING
7e1af563d8 Jean*0340        IF ( useRealFreshWaterFlux )
09e3b53265 Mart*0341      &      _EXCH_XY_RS( sIceLoad, myThid )
f3ce416a61 Jean*0342 #endif
809c36b928 Patr*0343
4f6e464f38 Jean*0344 #ifdef ALLOW_OBCS
3352db77f4 Jean*0345 C--   In case we use scheme with a large stencil that extends into overlap:
                0346 C     no longer needed with the right masking in advection & diffusion S/R.
                0347 c      IF ( useOBCS ) THEN
                0348 c       DO bj=myByLo(myThid),myByHi(myThid)
                0349 c        DO bi=myBxLo(myThid),myBxHi(myThid)
d24acf3cfc Jean*0350 c          CALL OBCS_COPY_TRACER( HEFF(1-OLx,1-OLy,bi,bj),
3352db77f4 Jean*0351 c    I                            1, bi, bj, myThid )
d24acf3cfc Jean*0352 c          CALL OBCS_COPY_TRACER( AREA(1-OLx,1-OLy,bi,bj),
3352db77f4 Jean*0353 c    I                            1, bi, bj, myThid )
d24acf3cfc Jean*0354 c          CALL OBCS_COPY_TRACER( HSNOW(1-OLx,1-OLy,bi,bj),
3352db77f4 Jean*0355 c    I                            1, bi, bj, myThid )
a98c4b8072 Ian *0356 #ifdef SEAICE_VARIABLE_SALINITY
d24acf3cfc Jean*0357 c          CALL OBCS_COPY_TRACER( HSALT(1-OLx,1-OLy,bi,bj),
3352db77f4 Jean*0358 c    I                            1, bi, bj, myThid )
4f6e464f38 Jean*0359 #endif
3352db77f4 Jean*0360 c        ENDDO
                0361 c       ENDDO
                0362 c      ENDIF
4f6e464f38 Jean*0363 #endif /* ALLOW_OBCS */
                0364
d74b8b60ab Patr*0365 #ifdef ALLOW_DIAGNOSTICS
                0366        IF ( useDiagnostics ) THEN
9ef71a93d2 Mart*0367 C     diagnostics for "non-state variables" that are modified by
d5254b4e3d Mart*0368 C     the seaice model ...
ead4e7560e Jean*0369         CALL DIAGNOSTICS_FILL_RS(EmPmR,'SIempmr ',0,1 ,0,1,1,myThid)
                0370         CALL DIAGNOSTICS_FILL_RS(Qnet ,'SIqnet  ',0,1 ,0,1,1,myThid)
                0371         CALL DIAGNOSTICS_FILL_RS(Qsw  ,'SIqsw   ',0,1 ,0,1,1,myThid)
d5254b4e3d Mart*0372 C     ... and energy diagnostic
5fe78992ba Mart*0373         IF ( diag_SIenph_isOn ) THEN
                0374          recip_deltaTtherm = 1. _d 0 / SEAICE_deltaTtherm
                0375          DO bj=myByLo(myThid),myByHi(myThid)
                0376           DO bi=myBxLo(myThid),myBxHi(myThid)
                0377            DO j=1,sNy
                0378             DO i=1,sNx
ee480e7782 amc2*0379              mTmp(i,j,bi,bj) = 0.25 _d 0 * (
5fe78992ba Mart*0380      &              uIce (i,j,bi,bj)**2 + uIce (i+1,j,  bi,bj)**2
                0381      &            + vIce (i,j,bi,bj)**2 + vIce (i,  j+1,bi,bj)**2 )
ee480e7782 amc2*0382      &            * (SEAICE_rhoIce*HEFF(i,j,bi,bj)
                0383      &            + addSnow*SEAICE_rhoSnow*HSNOW(i,j,bi,bj)
                0384      &            - mTmp(i,j,bi,bj))
5fe78992ba Mart*0385      &             *recip_deltaTtherm
                0386             ENDDO
                0387            ENDDO
                0388           ENDDO
                0389          ENDDO
ee480e7782 amc2*0390          CALL DIAGNOSTICS_FILL(mTmp,'SIenph  ',0,1,0,1,1,myThid)
5fe78992ba Mart*0391         ENDIF
d74b8b60ab Patr*0392        ENDIF
                0393 #endif /* ALLOW_DIAGNOSTICS */
                0394
09e3b53265 Mart*0395 #ifdef ALLOW_THSICE
                0396 C     endif .not.useThSice
                0397       ENDIF
                0398 #endif /* ALLOW_THSICE */
                0399 CML   This has already been done in seaice_ocean_stress/ostres, so why repeat?
                0400 CML   CALL EXCH_UV_XY_RS(fu,fv,.TRUE.,myThid)
                0401
fb1912d055 Patr*0402 #ifdef ALLOW_EXF
5d96c49326 Jean*0403 # ifdef ALLOW_AUTODIFF
                0404 #  ifdef ALLOW_AUTODIFF_MONITOR
d5254b4e3d Mart*0405       CALL EXF_ADJOINT_SNAPSHOTS( 3, myTime, myIter, myThid )
fb1912d055 Patr*0406 #  endif
                0407 # endif
                0408 #endif
                0409
b92ebad606 Jean*0410 #ifdef ALLOW_DEBUG
a8b75793c4 Jean*0411       IF (debugMode) CALL DEBUG_LEAVE( 'SEAICE_MODEL', myThid )
b92ebad606 Jean*0412 #endif
                0413
809c36b928 Patr*0414       RETURN
                0415       END