pkg/seaice/seaice_dynsolver.F

53092bcb42 Mart*0001 #include "SEAICE_OPTIONS.h"
772b2ed80e Gael*0002 #ifdef ALLOW_AUTODIFF
                0003 # include "AUTODIFF_OPTIONS.h"
                0004 #endif
53092bcb42 Mart*0005
91e72625af Jean*0006 CBOP
                0007 C     !ROUTINE: SEAICE_DYNSOLVER
                0008 C     !INTERFACE:
53092bcb42 Mart*0009       SUBROUTINE SEAICE_DYNSOLVER( myTime, myIter, myThid )
91e72625af Jean*0010
45315406aa Mart*0011 C     !DESCRIPTION:
                0012 C     *=============================================================*
                0013 C     | SUBROUTINE SEAICE_DYNSOLVER                                 |
                0014 C     | C-grid version of ice dynamics using either                 |
                0015 C     | o free drift                                                |
                0016 C     | o LSR solver, Zhang and Hibler, JGR, 102, 8691-8702, 1997   |
                0017 C     | o Krylov solver, after Lemieux and Tremblay, JGR, 114, 2009 |
                0018 C     | o JFNK solver, Losch et al., JCP, 257 901-911, 2014         |
                0019 C     | o EVP solver, Hunke and Dukowicz, JPO 27, 1849-1867 1997    |
                0020 C     *=============================================================*
91e72625af Jean*0021
                0022 C     !USES:
53092bcb42 Mart*0023       IMPLICIT NONE
                0024 C     === Global variables ===
                0025 #include "SIZE.h"
                0026 #include "EEPARAMS.h"
                0027 #include "PARAMS.h"
                0028 #include "GRID.h"
8468e0a1f9 Jean*0029 #include "SURFACE.h"
53092bcb42 Mart*0030 #include "DYNVARS.h"
                0031 #include "FFIELDS.h"
03c669d1ab Jean*0032 #include "SEAICE_SIZE.h"
53092bcb42 Mart*0033 #include "SEAICE_PARAMS.h"
03c669d1ab Jean*0034 #include "SEAICE.h"
53092bcb42 Mart*0035
                0036 #ifdef ALLOW_AUTODIFF_TAMC
                0037 # include "tamc.h"
                0038 #endif
                0039
45315406aa Mart*0040 C     !INPUT PARAMETERS:
53092bcb42 Mart*0041 C     === Routine arguments ===
91e72625af Jean*0042 C     myTime     :: Simulation time
                0043 C     myIter     :: Simulation timestep number
                0044 C     myThid     :: my Thread Id. number
53092bcb42 Mart*0045       _RL     myTime
                0046       INTEGER myIter
                0047       INTEGER myThid
                0048
91e72625af Jean*0049 C     !FUNCTIONS:
                0050       LOGICAL  DIFFERENT_MULTIPLE
                0051       EXTERNAL DIFFERENT_MULTIPLE
45315406aa Mart*0052 #if ( defined ALLOW_DIAGNOSTICS && defined SEAICE_CGRID )
5fe78992ba Mart*0053       LOGICAL  DIAGNOSTICS_IS_ON
                0054       EXTERNAL DIAGNOSTICS_IS_ON
45315406aa Mart*0055 #endif
91e72625af Jean*0056
                0057 C     !LOCAL VARIABLES:
53092bcb42 Mart*0058 C     === Local variables ===
17bd6b9422 jm-c 0059 C     i,j    :: Loop counters
                0060 C     bi,bj  :: tile counters
7abe6d1375 Mart*0061       INTEGER i, j, bi, bj
45315406aa Mart*0062 #ifdef SEAICE_CGRID
8377b8ee87 Mart*0063 # ifndef ALLOW_AUTODIFF
17bd6b9422 jm-c 0064       _RL mask_uice, mask_vice
8377b8ee87 Mart*0065 # endif
                0066 C     phiSurf :: geopotential height at sea surface (including pressure load)
                0067       _RL phiSurf(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0068 #endif
                0069 C     TAUX    :: zonal      wind stress over seaice at U point
                0070 C     TAUY    :: meridional wind stress over seaice at V point
                0071       _RL TAUX   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0072       _RL TAUY   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
45315406aa Mart*0073 #if ( defined ALLOW_DIAGNOSTICS && defined SEAICE_CGRID )
5fe78992ba Mart*0074 # ifdef ALLOW_AUTODIFF
                0075       _RL strDivX(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0076       _RL strDivY(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0077 # endif /* ALLOW_AUTODIFF */
                0078       _RL sig1   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0079       _RL sig2   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0080       _RL sig12  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0081       _RL sigp, sigm, sigTmp, recip_prs
99da6b3183 Jean*0082       _RL areaW, areaS, COSWAT
                0083       _RS SINWAT
5fe78992ba Mart*0084       INTEGER kSrf
                0085       LOGICAL diag_SIsigma_isOn, diag_SIshear_isOn
                0086       LOGICAL diag_SIenpi_isOn, diag_SIenpot_isOn
                0087       LOGICAL diag_SIpRfric_isOn, diag_SIpSfric_isOn
45315406aa Mart*0088 #endif
                0089 CEOP
17bd6b9422 jm-c 0090
fb1912d055 Patr*0091       DO bj=myByLo(myThid),myByHi(myThid)
                0092        DO bi=myBxLo(myThid),myBxHi(myThid)
5fe78992ba Mart*0093         DO j=1-OLy,sNy+OLy
                0094          DO i=1-OLx,sNx+OLx
45315406aa Mart*0095 #if ( defined ALLOW_AUTODIFF && defined SEAICE_CGRID )
8377b8ee87 Mart*0096 C Following re-initialisation breaks some "artificial" AD dependencies
                0097 C incured by IF (DIFFERENT_MULTIPLE ... statement
8e32c48b8f Mart*0098           PRESS0     (i,j,bi,bj) = SEAICE_strength*HEFF(i,j,bi,bj)
ba6cfc5714 Mart*0099      &         *EXP(-SEAICE_cStar*(ONE-AREA(i,j,bi,bj)))
8e32c48b8f Mart*0100           SEAICE_zMax(i,j,bi,bj) = SEAICE_zetaMaxFac*PRESS0(i,j,bi,bj)
                0101           SEAICE_zMin(i,j,bi,bj) = SEAICE_zetaMin
                0102           PRESS0     (i,j,bi,bj) = PRESS0(i,j,bi,bj)*HEFFM(i,j,bi,bj)
5fe78992ba Mart*0103 # ifdef SEAICE_ALLOW_FREEDRIFT
72c7f62b0e Patr*0104           uice_fd(i,j,bi,bj)= 0. _d 0
                0105           vice_fd(i,j,bi,bj)= 0. _d 0
5fe78992ba Mart*0106 # endif
45315406aa Mart*0107 #if ( defined ALLOW_DIAGNOSTICS && defined SEAICE_CGRID )
5fe78992ba Mart*0108           strDivX(i,j,bi,bj)= 0. _d 0
                0109           strDivY(i,j,bi,bj)= 0. _d 0
                0110 # endif
45315406aa Mart*0111 #endif /* ALLOW_AUTODIFF and SEAICE_CGRID */
8377b8ee87 Mart*0112 C     Always initialise these local variables, needed for TAF, but also
                0113 C     because they are not completely filled in S/R seaice_get_dynforcing
                0114           TAUX(i,j,bi,bj) = 0. _d 0
                0115           TAUY(i,j,bi,bj) = 0. _d 0
fb1912d055 Patr*0116          ENDDO
                0117         ENDDO
                0118        ENDDO
                0119       ENDDO
45315406aa Mart*0120 #ifdef ALLOW_AUTODIFF_TAMC
c20cddf271 Patr*0121 CADJ STORE uice    = comlev1, key=ikey_dynamics, kind=isbyte
                0122 CADJ STORE vice    = comlev1, key=ikey_dynamics, kind=isbyte
                0123 CADJ STORE uicenm1 = comlev1, key=ikey_dynamics, kind=isbyte
                0124 CADJ STORE vicenm1 = comlev1, key=ikey_dynamics, kind=isbyte
45315406aa Mart*0125 #endif /* ALLOW_AUTODIFF_TAMC */
8377b8ee87 Mart*0126 C--   interface of dynamics with atmopheric forcing fields (wind/stress)
                0127 C     Call this in each time step so that we can use the surface stress
                0128 C     in S/R seaice_ocean_stress
                0129       CALL SEAICE_GET_DYNFORCING (
                0130      I     uIce, vIce, AREA, SIMaskU, SIMaskV,
                0131      O     TAUX, TAUY,
                0132      I     myTime, myIter, myThid )
fb1912d055 Patr*0133
45315406aa Mart*0134 #ifdef SEAICE_CGRID
8377b8ee87 Mart*0135       IF ( SEAICEuseDYNAMICS .AND.
02aabb8fe5 Dimi*0136      &  DIFFERENT_MULTIPLE(SEAICE_deltaTdyn,myTime,SEAICE_deltaTtherm)
                0137      &   ) THEN
53092bcb42 Mart*0138
8377b8ee87 Mart*0139 #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_ALLOW_EVP)
8e32c48b8f Mart*0140 CADJ STORE press0      = comlev1, key=ikey_dynamics, kind=isbyte
                0141 CADJ STORE SEAICE_zMax = comlev1, key=ikey_dynamics, kind=isbyte
8377b8ee87 Mart*0142 #endif /* ALLOW_AUTODIFF_TAMC and SEAICE_ALLOW_EVP */
36feb88151 Patr*0143
53092bcb42 Mart*0144 C--   NOW SET UP MASS PER UNIT AREA AND CORIOLIS TERM
8377b8ee87 Mart*0145        DO bj=myByLo(myThid),myByHi(myThid)
                0146         DO bi=myBxLo(myThid),myBxHi(myThid)
79022779f5 Mart*0147          DO j=1-OLy+1,sNy+OLy
                0148           DO i=1-OLx+1,sNx+OLx
8377b8ee87 Mart*0149            seaiceMassC(i,j,bi,bj)=SEAICE_rhoIce*HEFF(i,j,bi,bj)
                0150            seaiceMassU(i,j,bi,bj)=SEAICE_rhoIce*HALF*(
                0151      &          HEFF(i,j,bi,bj) + HEFF(i-1,j  ,bi,bj) )
                0152            seaiceMassV(i,j,bi,bj)=SEAICE_rhoIce*HALF*(
                0153      &          HEFF(i,j,bi,bj) + HEFF(i  ,j-1,bi,bj) )
79022779f5 Mart*0154           ENDDO
                0155          ENDDO
8377b8ee87 Mart*0156          IF ( SEAICEaddSnowMass ) THEN
                0157           DO j=1-OLy+1,sNy+OLy
                0158            DO i=1-OLx+1,sNx+OLx
                0159             seaiceMassC(i,j,bi,bj)=seaiceMassC(i,j,bi,bj)
                0160      &           +                 SEAICE_rhoSnow*HSNOW(i,j,bi,bj)
                0161             seaiceMassU(i,j,bi,bj)=seaiceMassU(i,j,bi,bj)
                0162      &           +                  SEAICE_rhoSnow*HALF*(
                0163      &           HSNOW(i,j,bi,bj) + HSNOW(i-1,j  ,bi,bj) )
                0164
                0165             seaiceMassV(i,j,bi,bj)=seaiceMassV(i,j,bi,bj)
                0166      &           +                  SEAICE_rhoSnow*HALF*(
                0167      &           HSNOW(i,j,bi,bj) + HSNOW(i  ,j-1,bi,bj) )
                0168            ENDDO
                0169           ENDDO
                0170          ENDIF
                0171         ENDDO
53092bcb42 Mart*0172        ENDDO
                0173
45315406aa Mart*0174 # ifndef ALLOW_AUTODIFF
8377b8ee87 Mart*0175        IF ( SEAICE_maskRHS ) THEN
de5b48fa34 Mart*0176 C     dynamic masking of areas with no ice, not recommended
                0177 C     and only kept for testing purposes
8377b8ee87 Mart*0178         DO bj=myByLo(myThid),myByHi(myThid)
                0179          DO bi=myBxLo(myThid),myBxHi(myThid)
                0180           DO j=1-OLy+1,sNy+OLy
                0181            DO i=1-OLx+1,sNx+OLx
                0182             seaiceMaskU(i,j,bi,bj)=AREA(i,j,bi,bj)+AREA(i-1,j,bi,bj)
                0183             mask_uice=HEFFM(i,j,bi,bj)+HEFFM(i-1,j  ,bi,bj)
                0184             IF ( (seaiceMaskU(i,j,bi,bj) .GT. 0. _d 0) .AND.
                0185      &           (mask_uice .GT. 1.5 _d 0) ) THEN
                0186              seaiceMaskU(i,j,bi,bj) = 1. _d 0
                0187             ELSE
                0188              seaiceMaskU(i,j,bi,bj) = 0. _d 0
                0189             ENDIF
                0190             seaiceMaskV(i,j,bi,bj)=AREA(i,j,bi,bj)+AREA(i,j-1,bi,bj)
                0191             mask_vice=HEFFM(i,j,bi,bj)+HEFFM(i  ,j-1,bi,bj)
                0192             IF ( (seaiceMaskV(i,j,bi,bj) .GT. 0. _d 0) .AND.
                0193      &           (mask_vice .GT. 1.5 _d 0) ) THEN
                0194              seaiceMaskV(i,j,bi,bj) = 1. _d 0
                0195             ELSE
                0196              seaiceMaskV(i,j,bi,bj) = 0. _d 0
                0197             ENDIF
                0198            ENDDO
7abe6d1375 Mart*0199           ENDDO
                0200          ENDDO
                0201         ENDDO
8377b8ee87 Mart*0202         CALL EXCH_UV_XY_RL( seaiceMaskU, seaiceMaskV, .FALSE., myThid )
                0203        ENDIF
45315406aa Mart*0204 # endif /* ndef ALLOW_AUTODIFF */
7abe6d1375 Mart*0205
53092bcb42 Mart*0206 C--   NOW SET UP FORCING FIELDS
                0207
0dfb864288 Mart*0208 C     initialise fields
45315406aa Mart*0209 # if (defined ALLOW_AUTODIFF && defined SEAICE_ALLOW_EVP)
8377b8ee87 Mart*0210        DO bj=myByLo(myThid),myByHi(myThid)
                0211         DO bi=myBxLo(myThid),myBxHi(myThid)
                0212          DO j=1-OLy,sNy+OLy
                0213           DO i=1-OLx,sNx+OLx
                0214            stressDivergenceX(i,j,bi,bj) = 0. _d 0
                0215            stressDivergenceY(i,j,bi,bj) = 0. _d 0
                0216           ENDDO
484317af64 Mart*0217          ENDDO
                0218         ENDDO
                0219        ENDDO
45315406aa Mart*0220 # endif
53092bcb42 Mart*0221
8377b8ee87 Mart*0222        DO bj=myByLo(myThid),myByHi(myThid)
                0223         DO bi=myBxLo(myThid),myBxHi(myThid)
8468e0a1f9 Jean*0224 C--   Compute surface pressure at z==0:
                0225 C-    use actual sea surface height for tilt computations
8377b8ee87 Mart*0226          IF ( usingPCoords ) THEN
03c669d1ab Jean*0227           DO j=1-OLy,sNy+OLy
                0228            DO i=1-OLx,sNx+OLx
8377b8ee87 Mart*0229             phiSurf(i,j) = phiHydLow(i,j,bi,bj)
8468e0a1f9 Jean*0230            ENDDO
                0231           ENDDO
0320e25227 Mart*0232          ELSE
03c669d1ab Jean*0233           DO j=1-OLy,sNy+OLy
                0234            DO i=1-OLx,sNx+OLx
8377b8ee87 Mart*0235             phiSurf(i,j) = Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
8468e0a1f9 Jean*0236            ENDDO
                0237           ENDDO
0320e25227 Mart*0238          ENDIF
45315406aa Mart*0239 # ifdef ATMOSPHERIC_LOADING
8377b8ee87 Mart*0240 C--   add atmospheric loading and Sea-Ice loading as it is done for phi0surf
                0241 C--   in S/R external_forcing_surf
                0242          IF ( usingZCoords ) THEN
                0243           IF ( useRealFreshWaterFlux ) THEN
                0244            DO j=1-OLy,sNy+OLy
                0245             DO i=1-OLx,sNx+OLx
                0246              phiSurf(i,j) = phiSurf(i,j)
                0247      &                    + ( pload(i,j,bi,bj)
                0248      &                       +sIceLoad(i,j,bi,bj)*gravity*sIceLoadFac
                0249      &                      )*recip_rhoConst
                0250             ENDDO
                0251            ENDDO
                0252           ELSE
                0253            DO j=1-OLy,sNy+OLy
                0254             DO i=1-OLx,sNx+OLx
                0255              phiSurf(i,j) = phiSurf(i,j)
                0256      &                    + pload(i,j,bi,bj)*recip_rhoConst
                0257             ENDDO
                0258            ENDDO
                0259           ENDIF
                0260 C        ELSEIF ( usingPCoords ) THEN
0320e25227 Mart*0261 C     The true atmospheric P-loading is not yet implemented for P-coord
                0262 C     (requires time varying dP(Nr) like dP(k-bottom) with NonLin FS).
8377b8ee87 Mart*0263          ENDIF
45315406aa Mart*0264 # endif /* ATMOSPHERIC_LOADING */
21936d7dea Gael*0265 C--   basic forcing by wind stress
8377b8ee87 Mart*0266          IF ( SEAICEscaleSurfStress ) THEN
                0267           DO j=1-OLy+1,sNy+OLy
                0268            DO i=1-OLx+1,sNx+OLx
                0269             FORCEX0(i,j,bi,bj)=TAUX(i,j,bi,bj)
                0270      &           * 0.5 _d 0*(AREA(i,j,bi,bj)+AREA(i-1,j,bi,bj))
                0271             FORCEY0(i,j,bi,bj)=TAUY(i,j,bi,bj)
                0272      &           * 0.5 _d 0*(AREA(i,j,bi,bj)+AREA(i,j-1,bi,bj))
                0273            ENDDO
70e078b38a Mart*0274           ENDDO
8377b8ee87 Mart*0275          ELSE
                0276           DO j=1-OLy+1,sNy+OLy
                0277            DO i=1-OLx+1,sNx+OLx
                0278             FORCEX0(i,j,bi,bj)=TAUX(i,j,bi,bj)
                0279             FORCEY0(i,j,bi,bj)=TAUY(i,j,bi,bj)
                0280            ENDDO
70e078b38a Mart*0281           ENDDO
8377b8ee87 Mart*0282          ENDIF
21936d7dea Gael*0283
8377b8ee87 Mart*0284          IF ( SEAICEuseTILT ) THEN
                0285           DO j=1-OLy+1,sNy+OLy
                0286            DO i=1-OLx+1,sNx+OLx
21936d7dea Gael*0287 C--   now add in tilt
8377b8ee87 Mart*0288             FORCEX0(i,j,bi,bj)=FORCEX0(i,j,bi,bj)
                0289      &           -seaiceMassU(i,j,bi,bj)*_recip_dxC(i,j,bi,bj)
                0290      &           *( phiSurf(i,j)-phiSurf(i-1,j) )
                0291             FORCEY0(i,j,bi,bj)=FORCEY0(i,j,bi,bj)
                0292      &           -seaiceMassV(i,j,bi,bj)* _recip_dyC(i,j,bi,bj)
                0293      &           *( phiSurf(i,j)-phiSurf(i,j-1) )
                0294            ENDDO
                0295           ENDDO
                0296          ENDIF
000ae6c470 Mart*0297
8377b8ee87 Mart*0298          CALL SEAICE_CALC_ICE_STRENGTH( bi, bj, myTime, myIter, myThid )
0c32bd3cb0 Mart*0299
8377b8ee87 Mart*0300         ENDDO
53092bcb42 Mart*0301        ENDDO
                0302
45315406aa Mart*0303 # ifdef SEAICE_ALLOW_FREEDRIFT
2c255e54a7 Jean*0304        IF ( SEAICEuseFREEDRIFT .OR. SEAICEuseEVP
                0305      &                         .OR. LSR_mixIniGuess.EQ.0 ) THEN
                0306         CALL SEAICE_FREEDRIFT( myTime, myIter, myThid )
                0307        ENDIF
                0308        IF ( SEAICEuseFREEDRIFT ) THEN
                0309         DO bj=myByLo(myThid),myByHi(myThid)
                0310          DO bi=myBxLo(myThid),myBxHi(myThid)
                0311           DO j=1-OLy,sNy+OLy
                0312            DO i=1-OLx,sNx+OLx
                0313             uIce(i,j,bi,bj) = uIce_fd(i,j,bi,bj)
                0314             vIce(i,j,bi,bj) = vIce_fd(i,j,bi,bj)
                0315             stressDivergenceX(i,j,bi,bj) = 0. _d 0
                0316             stressDivergenceY(i,j,bi,bj) = 0. _d 0
                0317            ENDDO
                0318           ENDDO
                0319          ENDDO
                0320         ENDDO
                0321        ENDIF
45315406aa Mart*0322 # endif /* SEAICE_ALLOW_FREEDRIFT */
2c255e54a7 Jean*0323
45315406aa Mart*0324 # ifdef ALLOW_OBCS
91e72625af Jean*0325        IF ( useOBCS ) THEN
                0326          CALL OBCS_APPLY_UVICE( uIce, vIce, myThid )
                0327        ENDIF
45315406aa Mart*0328 # endif
91e72625af Jean*0329
45315406aa Mart*0330 # ifdef SEAICE_ALLOW_EVP
                0331 #  ifdef ALLOW_AUTODIFF_TAMC
c20cddf271 Patr*0332 CADJ STORE uice    = comlev1, key=ikey_dynamics, kind=isbyte
                0333 CADJ STORE vice    = comlev1, key=ikey_dynamics, kind=isbyte
                0334 CADJ STORE uicenm1 = comlev1, key=ikey_dynamics, kind=isbyte
                0335 CADJ STORE vicenm1 = comlev1, key=ikey_dynamics, kind=isbyte
45315406aa Mart*0336 #  endif /* ALLOW_AUTODIFF_TAMC */
6e2f4e58fa Mart*0337        IF ( SEAICEuseEVP ) THEN
b0fd37e69b Mart*0338 C     Elastic-Viscous-Plastic solver, following Hunke (2001)
6e2f4e58fa Mart*0339         CALL SEAICE_EVP( myTime, myIter, myThid )
2c255e54a7 Jean*0340        ENDIF
45315406aa Mart*0341 # endif /* SEAICE_ALLOW_EVP */
2c255e54a7 Jean*0342
c8739d4898 Mart*0343        IF ( SEAICEuseLSR ) THEN
10811105fa Mart*0344 C     Picard solver with LSR scheme (Zhang-J/Hibler 1997), ported to a C-grid
c8739d4898 Mart*0345         CALL SEAICE_LSR( myTime, myIter, myThid )
                0346        ENDIF
                0347
45315406aa Mart*0348 # ifdef SEAICE_ALLOW_KRYLOV
                0349 #  ifdef ALLOW_AUTODIFF
10811105fa Mart*0350        STOP 'Adjoint does not work with Picard-Krylov solver.'
45315406aa Mart*0351 #  else
10811105fa Mart*0352        IF ( SEAICEuseKrylov ) THEN
                0353 C     Picard solver with Matrix-free Krylov solver (Lemieux et al. 2008)
                0354         CALL SEAICE_KRYLOV( myTime, myIter, myThid )
                0355        ENDIF
45315406aa Mart*0356 #  endif /*  ALLOW_AUTODIFF */
                0357 # endif /* SEAICE_ALLOW_KRYLOV */
10811105fa Mart*0358
45315406aa Mart*0359 # ifdef SEAICE_ALLOW_JFNK
                0360 #  ifdef ALLOW_AUTODIFF
10811105fa Mart*0361        STOP 'Adjoint does not work with JFNK solver.'
45315406aa Mart*0362 #  else
10811105fa Mart*0363        IF ( SEAICEuseJFNK ) THEN
b0fd37e69b Mart*0364 C     Jacobian-free Newton Krylov solver (Lemieux et al. 2010, 2012)
e89e650bb4 Mart*0365         CALL SEAICE_JFNK( myTime, myIter, myThid )
                0366        ENDIF
45315406aa Mart*0367 #  endif /*  ALLOW_AUTODIFF */
                0368 # endif /* SEAICE_ALLOW_JFNK */
e89e650bb4 Mart*0369
8377b8ee87 Mart*0370 C End of IF (SEAICEuseDYNAMICS and DIFFERENT_MULTIPLE ...
53092bcb42 Mart*0371       ENDIF
45315406aa Mart*0372 #endif /* SEAICE_CGRID */
53092bcb42 Mart*0373
210ee8461e jm-c 0374 C Update ocean surface stress
                0375       IF ( SEAICEupdateOceanStress ) THEN
8377b8ee87 Mart*0376 #ifdef ALLOW_AUTODIFF_TAMC
                0377 CADJ STORE uice, vice, DWATN = comlev1, key=ikey_dynamics, kind=isbyte
45315406aa Mart*0378 # ifdef SEAICE_CGRID
8377b8ee87 Mart*0379 CADJ STORE stressDivergenceX = comlev1, key=ikey_dynamics, kind=isbyte
                0380 CADJ STORE stressDivergenceY = comlev1, key=ikey_dynamics, kind=isbyte
45315406aa Mart*0381 # endif
8377b8ee87 Mart*0382 #endif /* ALLOW_AUTODIFF_TAMC */
210ee8461e jm-c 0383         CALL SEAICE_OCEAN_STRESS (
                0384      I              TAUX, TAUY, myTime, myIter, myThid )
                0385       ENDIF
53092bcb42 Mart*0386
45315406aa Mart*0387 #ifdef SEAICE_ALLOW_CLIPVELS
7abe6d1375 Mart*0388       IF ( SEAICEuseDYNAMICS .AND. SEAICE_clipVelocities) THEN
8377b8ee87 Mart*0389 # ifdef ALLOW_AUTODIFF_TAMC
95c72ef3a1 Patr*0390 CADJ STORE uice = comlev1, key=ikey_dynamics, kind=isbyte
                0391 CADJ STORE vice = comlev1, key=ikey_dynamics, kind=isbyte
8377b8ee87 Mart*0392 # endif /* ALLOW_AUTODIFF_TAMC */
7abe6d1375 Mart*0393 c Put a cap on ice velocity
                0394 c limit velocity to 0.40 m s-1 to avoid potential CFL violations
                0395 c in open water areas (drift of zero thickness ice)
                0396        DO bj=myByLo(myThid),myByHi(myThid)
                0397         DO bi=myBxLo(myThid),myBxHi(myThid)
                0398          DO j=1-OLy,sNy+OLy
                0399           DO i=1-OLx,sNx+OLx
772590b63c Mart*0400            uIce(i,j,bi,bj)=
                0401      &          MAX(MIN(uIce(i,j,bi,bj),0.40 _d +00),-0.40 _d +00)
                0402            vIce(i,j,bi,bj)=
                0403      &          MAX(MIN(vIce(i,j,bi,bj),0.40 _d +00),-0.40 _d +00)
7abe6d1375 Mart*0404           ENDDO
                0405          ENDDO
                0406         ENDDO
                0407        ENDDO
                0408       ENDIF
45315406aa Mart*0409 #endif /* SEAICE_ALLOW_CLIPVELS */
f0d90fb111 Jean*0410
45315406aa Mart*0411 #if ( defined ALLOW_DIAGNOSTICS && defined SEAICE_CGRID )
5fe78992ba Mart*0412 C     diagnostics related to mechanics/dynamics/momentum equations
                0413       IF ( useDiagnostics .AND. SEAICEuseDYNAMICS ) THEN
                0414        CALL DIAGNOSTICS_FILL(zeta   ,'SIzeta  ',0,1,0,1,1,myThid)
                0415        CALL DIAGNOSTICS_FILL(eta    ,'SIeta   ',0,1,0,1,1,myThid)
                0416        CALL DIAGNOSTICS_FILL(press  ,'SIpress ',0,1,0,1,1,myThid)
                0417        CALL DIAGNOSTICS_FILL(deltaC ,'SIdelta ',0,1,0,1,1,myThid)
5bb179ddc2 Mart*0418 # ifdef SEAICE_ALLOW_SIDEDRAG
                0419 C     recompute lateral coast drag terms
                0420        IF ( DIAGNOSTICS_IS_ON('SIlatDgU',myThid) ) THEN
                0421         DO bj = myByLo(myThid), myByHi(myThid)
                0422          DO bi = myBxLo(myThid), myBxHi(myThid)
                0423 C     use sig1 as a temporary field
                0424           DO j=1,sNy
                0425            DO i=1,sNx+1
                0426             sig1(i,j) = sideDragU(i,j,bi,bj)*uIce(i,j,bi,bj)
                0427            ENDDO
                0428           ENDDO
                0429           CALL DIAGNOSTICS_FILL(sig1,'SIlatDgU',0,1,2,bi,bj,myThid)
                0430          ENDDO
                0431         ENDDO
                0432        ENDIF
                0433        IF ( DIAGNOSTICS_IS_ON('SIlatDgV',myThid) ) THEN
                0434         DO bj = myByLo(myThid), myByHi(myThid)
                0435          DO bi = myBxLo(myThid), myBxHi(myThid)
                0436 C     use sig1 as a temporary field
                0437           DO j=1,sNy+1
                0438            DO i=1,sNx
                0439             sig1(i,j) = sideDragV(i,j,bi,bj)*vIce(i,j,bi,bj)
                0440            ENDDO
                0441           ENDDO
                0442           CALL DIAGNOSTICS_FILL(sig1,'SIlatDgV',0,1,2,bi,bj,myThid)
                0443          ENDDO
                0444         ENDDO
                0445        ENDIF
                0446 # endif /* SEAICE_ALLOW_SIDEDRAG */
                0447
5fe78992ba Mart*0448        IF ( DIAGNOSTICS_IS_ON('SItensil',myThid) ) THEN
                0449         DO bj = myByLo(myThid), myByHi(myThid)
                0450          DO bi = myBxLo(myThid), myBxHi(myThid)
                0451 C     use sig1 as a temporary field
                0452           DO j=1,sNy
                0453            DO i=1,sNx
                0454             IF ( tensileStrFac(i,j,bi,bj) .EQ. oneRL ) THEN
                0455 C     This special case of tensile strength equal to compressive strength
                0456 C     is not very physical and should actually not happen but you never know;
                0457 C     in this case, press = P-T = P*(1-k) = 0. and we have to use press0 to
                0458 C     get something
                0459              sig1(i,j) = press0(i,j,bi,bj)
                0460             ELSE
                0461 C     This is more complicated than you think because press = P-T = P*(1-k),
                0462 C     but we are looking for T = k*P = k*press/(1-k)
                0463              sig1(i,j) = tensileStrFac(i,j,bi,bj)
                0464      &            *press(i,j,bi,bj)/(1. _d 0 - tensileStrFac(i,j,bi,bj))
                0465             ENDIF
                0466            ENDDO
                0467           ENDDO
                0468           CALL DIAGNOSTICS_FILL(sig1,'SItensil',0,1,2,bi,bj,myThid)
                0469          ENDDO
                0470         ENDDO
                0471        ENDIF
                0472 C     If any of the stress or energy diagnostics are required,
                0473 C     first recompute strainrates from up-to-date velocities
                0474        diag_SIsigma_isOn  = DIAGNOSTICS_IS_ON('SIsig1  ',myThid)
                0475      &                 .OR. DIAGNOSTICS_IS_ON('SIsig2  ',myThid)
                0476        diag_SIshear_isOn  = DIAGNOSTICS_IS_ON('SIshear ',myThid)
                0477        diag_SIenpi_isOn   = DIAGNOSTICS_IS_ON('SIenpi  ',myThid)
                0478        diag_SIenpot_isOn  = DIAGNOSTICS_IS_ON('SIenpot ',myThid)
                0479        diag_SIpRfric_isOn = DIAGNOSTICS_IS_ON('SIpRfric',myThid)
                0480        diag_SIpSfric_isOn = DIAGNOSTICS_IS_ON('SIpSfric',myThid)
                0481        IF ( diag_SIsigma_isOn  .OR. diag_SIshear_isOn .OR.
                0482      &      diag_SIenpi_isOn   .OR. diag_SIenpot_isOn .OR.
                0483      &      diag_SIpRfric_isOn .OR. diag_SIpSfric_isOn ) THEN
                0484         CALL SEAICE_CALC_STRAINRATES(
                0485      I       uIce, vIce,
                0486      O       e11, e22, e12,
                0487      I       0, myTime, myIter, myThid )
                0488 C     but use old viscosities and pressure for the
                0489 C     principle stress components
                0490 CML     CALL SEAICE_CALC_VISCOSITIES(
8e32c48b8f Mart*0491 CML  I     e11, e22, e12, SEAICE_zMin, SEAICE_zMax, HEFFM, press0,
                0492 CML  I     tensileStrFac,
5fe78992ba Mart*0493 CML  O     eta, etaZ, zeta, zetaZ, press, deltaC,
                0494 CML  I     0, myTime, myIter, myThid )
                0495        ENDIF
                0496 C
                0497        DO bj = myByLo(myThid), myByHi(myThid)
                0498         DO bi = myBxLo(myThid), myBxHi(myThid)
                0499 C
                0500 C     stress diagnostics
                0501 C
                0502          IF ( diag_SIsigma_isOn ) THEN
8377b8ee87 Mart*0503 # ifdef SEAICE_ALLOW_EVP
5fe78992ba Mart*0504 C     This could go directly into EVP, but to keep a better eye on it,
                0505 C     I would like to keep it here.
                0506           IF ( SEAICEuseEVP ) THEN
                0507 C     for EVP compute principle stress components from recent
                0508 C     stress state and normalize with latest
                0509 C     PRESS = PRESS(n-1), n = number of sub-cycling steps
                0510            DO j=1,sNy
                0511             DO i=1,sNx
                0512              sigp = seaice_sigma1(i,j,bi,bj)
                0513              sigm = seaice_sigma2(i,j,bi,bj)
                0514              sig12(i,j) = 0.25 _d 0 *
                0515      &            ( seaice_sigma12(i  ,j  ,bi,bj)
                0516      &            + seaice_sigma12(i+1,j  ,bi,bj)
                0517      &            + seaice_sigma12(i+1,j+1,bi,bj)
                0518      &            + seaice_sigma12(i  ,j+1,bi,bj) )
                0519              sigTmp = SQRT( sigm*sigm + 4. _d 0*sig12(i,j)*sig12(i,j) )
                0520              recip_prs = 0. _d 0
                0521              IF ( press0(i,j,bi,bj) .GT. 1. _d -13 )
                0522      &            recip_prs = 1. _d 0 / press0(i,j,bi,bj)
                0523              sig1(i,j) = 0.5 _d 0*(sigp + sigTmp)*recip_prs
                0524              sig2(i,j) = 0.5 _d 0*(sigp - sigTmp)*recip_prs
                0525             ENDDO
                0526            ENDDO
                0527           ELSE
8377b8ee87 Mart*0528 # else
5fe78992ba Mart*0529           IF ( .TRUE. ) THEN
8377b8ee87 Mart*0530 # endif /* SEAICE_ALLOW_EVP */
5fe78992ba Mart*0531            CALL SEAICE_CALC_STRESS(
                0532      I          e11, e22, e12, press, zeta, eta, etaZ,
                0533      O          sig1, sig2, sig12,
                0534      I          bi, bj, myTime, myIter, myThid )
                0535            DO j=1,sNy
                0536             DO i=1,sNx
                0537              sigp   = sig1(i,j) + sig2(i,j)
                0538              sigm   = sig1(i,j) - sig2(i,j)
                0539 C     This should be the way of computing sig12 at C-points,
                0540 C            sigTmp = 0.25 _d 0 *
                0541 C    &            ( sig12(i  ,j  ) + sig12(i+1,j  )
                0542 C    &            + sig12(i  ,j+1) + sig12(i+1,j+1) )
                0543 C     but sig12 = 2*etaZ*e12, and because of strong gradients in eta,
                0544 C     etaZ can be very large for a cell with small eta and the straightforward
                0545 C     way of averaging mixes large etaZ with small press0, so we have to do it
                0546 C     in different way to get meaningfull sig12C (=sigTmp):
                0547              sigTmp = 2. _d 0 * eta(i,j,bi,bj) * 0.25 _d 0 *
                0548      &            (e12(i,j,bi,bj) + e12(i+1,j,bi,bj)
                0549      &            +e12(i,j+1,bi,bj)+e12(i+1,j+1,bi,bj))
                0550              sigTmp = SQRT( sigm*sigm + 4. _d 0*sigTmp*sigTmp )
                0551              recip_prs = 0. _d 0
                0552              IF ( press0(i,j,bi,bj) .GT. 1. _d -13 )
                0553      &            recip_prs = 1. _d 0 / press0(i,j,bi,bj)
                0554              sig1(i,j) = 0.5 _d 0*(sigp + sigTmp)*recip_prs
                0555              sig2(i,j) = 0.5 _d 0*(sigp - sigTmp)*recip_prs
                0556             ENDDO
                0557            ENDDO
                0558           ENDIF
                0559           CALL DIAGNOSTICS_FILL(sig1,'SIsig1  ',0,1,2,bi,bj,myThid)
                0560           CALL DIAGNOSTICS_FILL(sig2,'SIsig2  ',0,1,2,bi,bj,myThid)
                0561          ENDIF
                0562 C
                0563          IF ( diag_SIshear_isOn  ) THEN
                0564           DO j=1,sNy
                0565            DO i=1,sNx
                0566             sigm = e11(i,j,bi,bj) - e22(i,j,bi,bj)
                0567             sigTmp =
                0568      &           ( e12(i  ,j  ,bi,bj)**2 + e12(i+1,j  ,bi,bj)**2
                0569      &           + e12(i+1,j+1,bi,bj)**2 + e12(i  ,j+1,bi,bj)**2 )
                0570 C     shear deformation as sqrt((e11-e22)**2 + 4*e12**2); the 4 pulled into
                0571 C     the average
8377b8ee87 Mart*0572             sig1(i,j) = SQRT(sigm*sigm + sigTmp)
5fe78992ba Mart*0573            ENDDO
                0574           ENDDO
                0575           CALL DIAGNOSTICS_FILL(sig1,'SIshear ',0,1,2,bi,bj,myThid)
                0576          ENDIF
                0577 C
                0578 C     most of the energy diagnostics, re-use sig1, sig2, sig12 as tmp-arrays
                0579 C
                0580          IF ( diag_SIenpi_isOn ) THEN
                0581 C     compute internal stresses with updated ice velocities
                0582 C     TAF gets confused when we use stressDivergenceX/Y as temporary arrays
                0583 C     therefore we need to use new arrays
                0584           IF ( .NOT. SEAICEuseFREEDRIFT )
                0585      &         CALL SEAICE_CALC_STRESSDIV(
                0586      I         e11, e22, e12, press, zeta, eta, etaZ,
                0587 # ifdef ALLOW_AUTODIFF
                0588      O         strDivX, strDivY,
                0589 # else /* not ALLOW_AUTODIFF */
                0590      O         stressDivergenceX, stressDivergenceY,
                0591 # endif /* ALLOW_AUTODIFF */
                0592      I         bi, bj, myTime, myIter, myThid )
                0593           DO j=1,sNy+1
                0594            DO i=1,sNx+1
                0595 # ifdef ALLOW_AUTODIFF
                0596             sig1(i,j) = uIce(i,j,bi,bj)*strDivX(i,j,bi,bj)
                0597             sig2(i,j) = vIce(i,j,bi,bj)*strDivY(i,j,bi,bj)
                0598 # else /* not ALLOW_AUTODIFF */
                0599             sig1(i,j) = uIce(i,j,bi,bj)*stressDivergenceX(i,j,bi,bj)
                0600             sig2(i,j) = vIce(i,j,bi,bj)*stressDivergenceY(i,j,bi,bj)
                0601 # endif /* ALLOW_AUTODIFF */
                0602            ENDDO
                0603           ENDDO
                0604 C     average from velocity points to pressure points
                0605           DO j=1,sNy
                0606            DO i=1,sNx
                0607             sig12(i,j)= 0.5 _d 0 * ( sig1(i,j) + sig1(i+1,j)
                0608      &                             + sig2(i,j) + sig2(i,j+1) )
                0609            ENDDO
                0610           ENDDO
                0611           CALL DIAGNOSTICS_FILL(sig12,'SIenpi  ',0,1,2,bi,bj,myThid)
                0612          ENDIF
                0613          IF ( diag_SIenpot_isOn ) THEN
                0614           DO j=1,sNy
                0615            DO i=1,sNx
                0616             sig1(i,j) = 0.5 _d 0 * press0(i,j,bi,bj) *
                0617      &           ( e11(i,j,bi,bj) + e22(i,j,bi,bj) )
                0618            ENDDO
                0619           ENDDO
                0620           CALL DIAGNOSTICS_FILL(sig1,'SIenpot ',0,1,2,bi,bj,myThid)
                0621          ENDIF
                0622          IF ( diag_SIpRfric_isOn ) THEN
                0623           DO j=1,sNy
                0624            DO i=1,sNx
                0625             sig1(i,j) = - zeta(i,j,bi,bj)
                0626      &           * ( e11(i,j,bi,bj) + e22(i,j,bi,bj) )
                0627      &           * ( e11(i,j,bi,bj) + e22(i,j,bi,bj) )
                0628            ENDDO
                0629           ENDDO
                0630           CALL DIAGNOSTICS_FILL(sig1,'SIpRfric',0,1,2,bi,bj,myThid)
                0631          ENDIF
                0632          IF ( diag_SIpSfric_isOn ) THEN
                0633           DO j=1,sNy
                0634            DO i=1,sNx
                0635             sig1(i,j) = - eta(i,j,bi,bj)
                0636      &        * ( (e11(i,j,bi,bj) - e22(i,j,bi,bj))**2 +
                0637      &            ( e12(i  ,j  ,bi,bj)**2 + e12(i+1,j  ,bi,bj)**2
                0638      &            + e12(i+1,j+1,bi,bj)**2 + e12(i  ,j+1,bi,bj)**2 )
                0639      &          )
                0640            ENDDO
                0641           ENDDO
                0642           CALL DIAGNOSTICS_FILL(sig1,'SIpSfric',0,1,2,bi,bj,myThid)
                0643          ENDIF
                0644          IF ( DIAGNOSTICS_IS_ON('SIenpa  ',myThid) ) THEN
                0645           DO j=1,sNy+1
                0646            DO i=1,sNx+1
                0647             areaW = 0.5 _d 0 * (AREA(i,j,bi,bj) + AREA(i-1,j,bi,bj))
                0648      &           * SEAICEstressFactor
                0649             areaS = 0.5 _d 0 * (AREA(i,j,bi,bj) + AREA(i,j-1,bi,bj))
                0650      &           * SEAICEstressFactor
                0651             sig1(i,j) = TAUX(i,j,bi,bj)*areaW * uIce(i,j,bi,bj)
                0652             sig2(i,j) = TAUY(i,j,bi,bj)*areaS * vIce(i,j,bi,bj)
                0653            ENDDO
                0654           ENDDO
                0655 C     average from velocity points to pressure points
                0656           DO j=1,sNy
                0657            DO i=1,sNx
                0658             sig12(i,j)= 0.5 _d 0 * ( sig1(i,j) + sig1(i+1,j)
                0659      &                             + sig2(i,j) + sig2(i,j+1) )
                0660            ENDDO
                0661           ENDDO
                0662           CALL DIAGNOSTICS_FILL(sig12,'SIenpa  ',0,1,2,bi,bj,myThid)
                0663          ENDIF
                0664          IF ( DIAGNOSTICS_IS_ON('SIenpw  ',myThid) ) THEN
                0665 C     surface level
                0666           kSrf = 1
                0667 C     introduce turning angle (default is zero)
                0668           SINWAT=SIN(SEAICE_waterTurnAngle*deg2rad)
                0669           COSWAT=COS(SEAICE_waterTurnAngle*deg2rad)
                0670           DO j=1,sNy+1
                0671            DO i=1,sNx+1
                0672             areaW = 0.5 _d 0 * (AREA(i,j,bi,bj) + AREA(i-1,j,bi,bj))
                0673      &         * SEAICEstressFactor
                0674             areaS = 0.5 _d 0 * (AREA(i,j,bi,bj) + AREA(i,j-1,bi,bj))
                0675      &         * SEAICEstressFactor
                0676             sig1(i,j) = areaW *
                0677      &           HALF*( DWATN(i,j,bi,bj)+DWATN(i-1,j,bi,bj) )*
                0678      &         COSWAT *
                0679      &         ( uIce(i,j,bi,bj)-uVel(i,j,kSrf,bi,bj) )
                0680      &         - SIGN(SINWAT, _fCori(i,j,bi,bj)) * 0.5 _d 0 *
                0681      &         ( DWATN(i  ,j,bi,bj) *
                0682      &         0.5 _d 0*(vIce(i  ,j  ,bi,bj)-vVel(i  ,j  ,kSrf,bi,bj)
                0683      &                  +vIce(i  ,j+1,bi,bj)-vVel(i  ,j+1,kSrf,bi,bj))
                0684      &         + DWATN(i-1,j,bi,bj) *
                0685      &         0.5 _d 0*(vIce(i-1,j  ,bi,bj)-vVel(i-1,j  ,kSrf,bi,bj)
                0686      &                  +vIce(i-1,j+1,bi,bj)-vVel(i-1,j+1,kSrf,bi,bj))
                0687      &         )
                0688             sig1(i,j) = sig1(i,j) * uIce(i,j,bi,bj)
                0689             sig2(i,j) = areaS *
                0690      &           HALF*( DWATN(i,j,bi,bj)+DWATN(i,j-1,bi,bj) )*
                0691      &         COSWAT *
                0692      &         ( vIce(i,j,bi,bj)-vVel(i,j,kSrf,bi,bj) )
                0693      &         + SIGN(SINWAT,  _fCori(i,j,bi,bj)) * 0.5 _d 0 *
                0694      &         ( DWATN(i,j  ,bi,bj) *
                0695      &         0.5 _d 0*(uIce(i  ,j  ,bi,bj)-uVel(i  ,j  ,kSrf,bi,bj)
                0696      &                  +uIce(i+1,j  ,bi,bj)-uVel(i+1,j  ,kSrf,bi,bj))
                0697      &         + DWATN(i,j-1,bi,bj) *
                0698      &         0.5 _d 0*(uIce(i  ,j-1,bi,bj)-uVel(i  ,j-1,kSrf,bi,bj)
                0699      &                  +uIce(i+1,j-1,bi,bj)-uVel(i+1,j-1,kSrf,bi,bj))
                0700      &         )
                0701             sig2(i,j) = sig2(i,j) * vIce(i,j,bi,bj)
                0702            ENDDO
                0703           ENDDO
                0704 C     average from velocity points to pressure points
                0705           DO j=1,sNy
                0706            DO i=1,sNx
                0707             sig12(i,j)= - 0.5 _d 0 * ( sig1(i,j) + sig1(i+1,j)
                0708      &                               + sig2(i,j) + sig2(i,j+1) )
                0709            ENDDO
                0710           ENDDO
                0711           CALL DIAGNOSTICS_FILL(sig12,'SIenpw  ',0,1,2,bi,bj,myThid)
                0712          ENDIF
                0713          IF ( SEAICEuseTilt
                0714      &        .AND. DIAGNOSTICS_IS_ON('SIenpg  ',myThid) ) THEN
                0715           DO j=1-OLy+1,sNy+OLy
                0716            DO i=1-OLx+1,sNx+OLx
                0717             sig1(i,j)=
                0718      &           - seaiceMassU(i,j,bi,bj)*_recip_dxC(i,j,bi,bj)
                0719      &           *( phiSurf(i,j)-phiSurf(i-1,j) ) * uIce(i,j,bi,bj)
                0720             sig2(i,j)=
                0721      &           - seaiceMassV(i,j,bi,bj)* _recip_dyC(i,j,bi,bj)
                0722      &           *( phiSurf(i,j)-phiSurf(i,j-1) ) * vIce(i,j,bi,bj)
                0723            ENDDO
                0724           ENDDO
                0725 C     average from velocity points to pressure points
                0726           DO j=1,sNy
                0727            DO i=1,sNx
                0728             sig12(i,j) = 0.5 _d 0 * ( sig1(i,j) + sig1(i+1,j)
                0729      &                              + sig2(i,j) + sig2(i,j+1) )
                0730            ENDDO
                0731           ENDDO
                0732           CALL DIAGNOSTICS_FILL(sig12,'SIenpg  ',0,1,2,bi,bj,myThid)
                0733          ENDIF
                0734 C     bi/bj-loop
                0735         ENDDO
                0736        ENDDO
45315406aa Mart*0737 C     useDiagnostics & SEAICEuseDYNAMICS
5fe78992ba Mart*0738       ENDIF
45315406aa Mart*0739 #endif /* ALLOW_DIAGNOSTICS and SEAICE_CGRID */
5fe78992ba Mart*0740
53092bcb42 Mart*0741       RETURN
                0742       END