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