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)
                0418        IF ( DIAGNOSTICS_IS_ON('SItensil',myThid) ) THEN
                0419         DO bj = myByLo(myThid), myByHi(myThid)
                0420          DO bi = myBxLo(myThid), myBxHi(myThid)
                0421 C     use sig1 as a temporary field
                0422           DO j=1,sNy
                0423            DO i=1,sNx
                0424             IF ( tensileStrFac(i,j,bi,bj) .EQ. oneRL ) THEN
                0425 C     This special case of tensile strength equal to compressive strength
                0426 C     is not very physical and should actually not happen but you never know;
                0427 C     in this case, press = P-T = P*(1-k) = 0. and we have to use press0 to
                0428 C     get something
                0429              sig1(i,j) = press0(i,j,bi,bj)
                0430             ELSE
                0431 C     This is more complicated than you think because press = P-T = P*(1-k),
                0432 C     but we are looking for T = k*P = k*press/(1-k)
                0433              sig1(i,j) = tensileStrFac(i,j,bi,bj)
                0434      &            *press(i,j,bi,bj)/(1. _d 0 - tensileStrFac(i,j,bi,bj))
                0435             ENDIF
                0436            ENDDO
                0437           ENDDO
                0438           CALL DIAGNOSTICS_FILL(sig1,'SItensil',0,1,2,bi,bj,myThid)
                0439          ENDDO
                0440         ENDDO
                0441        ENDIF
                0442 C     If any of the stress or energy diagnostics are required,
                0443 C     first recompute strainrates from up-to-date velocities
                0444        diag_SIsigma_isOn  = DIAGNOSTICS_IS_ON('SIsig1  ',myThid)
                0445      &                 .OR. DIAGNOSTICS_IS_ON('SIsig2  ',myThid)
                0446        diag_SIshear_isOn  = DIAGNOSTICS_IS_ON('SIshear ',myThid)
                0447        diag_SIenpi_isOn   = DIAGNOSTICS_IS_ON('SIenpi  ',myThid)
                0448        diag_SIenpot_isOn  = DIAGNOSTICS_IS_ON('SIenpot ',myThid)
                0449        diag_SIpRfric_isOn = DIAGNOSTICS_IS_ON('SIpRfric',myThid)
                0450        diag_SIpSfric_isOn = DIAGNOSTICS_IS_ON('SIpSfric',myThid)
                0451        IF ( diag_SIsigma_isOn  .OR. diag_SIshear_isOn .OR.
                0452      &      diag_SIenpi_isOn   .OR. diag_SIenpot_isOn .OR.
                0453      &      diag_SIpRfric_isOn .OR. diag_SIpSfric_isOn ) THEN
                0454         CALL SEAICE_CALC_STRAINRATES(
                0455      I       uIce, vIce,
                0456      O       e11, e22, e12,
                0457      I       0, myTime, myIter, myThid )
                0458 C     but use old viscosities and pressure for the
                0459 C     principle stress components
                0460 CML     CALL SEAICE_CALC_VISCOSITIES(
8e32c48b8f Mart*0461 CML  I     e11, e22, e12, SEAICE_zMin, SEAICE_zMax, HEFFM, press0,
                0462 CML  I     tensileStrFac,
5fe78992ba Mart*0463 CML  O     eta, etaZ, zeta, zetaZ, press, deltaC,
                0464 CML  I     0, myTime, myIter, myThid )
                0465        ENDIF
                0466 C
                0467        DO bj = myByLo(myThid), myByHi(myThid)
                0468         DO bi = myBxLo(myThid), myBxHi(myThid)
                0469 C
                0470 C     stress diagnostics
                0471 C
                0472          IF ( diag_SIsigma_isOn ) THEN
8377b8ee87 Mart*0473 # ifdef SEAICE_ALLOW_EVP
5fe78992ba Mart*0474 C     This could go directly into EVP, but to keep a better eye on it,
                0475 C     I would like to keep it here.
                0476           IF ( SEAICEuseEVP ) THEN
                0477 C     for EVP compute principle stress components from recent
                0478 C     stress state and normalize with latest
                0479 C     PRESS = PRESS(n-1), n = number of sub-cycling steps
                0480            DO j=1,sNy
                0481             DO i=1,sNx
                0482              sigp = seaice_sigma1(i,j,bi,bj)
                0483              sigm = seaice_sigma2(i,j,bi,bj)
                0484              sig12(i,j) = 0.25 _d 0 *
                0485      &            ( seaice_sigma12(i  ,j  ,bi,bj)
                0486      &            + seaice_sigma12(i+1,j  ,bi,bj)
                0487      &            + seaice_sigma12(i+1,j+1,bi,bj)
                0488      &            + seaice_sigma12(i  ,j+1,bi,bj) )
                0489              sigTmp = SQRT( sigm*sigm + 4. _d 0*sig12(i,j)*sig12(i,j) )
                0490              recip_prs = 0. _d 0
                0491              IF ( press0(i,j,bi,bj) .GT. 1. _d -13 )
                0492      &            recip_prs = 1. _d 0 / press0(i,j,bi,bj)
                0493              sig1(i,j) = 0.5 _d 0*(sigp + sigTmp)*recip_prs
                0494              sig2(i,j) = 0.5 _d 0*(sigp - sigTmp)*recip_prs
                0495             ENDDO
                0496            ENDDO
                0497           ELSE
8377b8ee87 Mart*0498 # else
5fe78992ba Mart*0499           IF ( .TRUE. ) THEN
8377b8ee87 Mart*0500 # endif /* SEAICE_ALLOW_EVP */
5fe78992ba Mart*0501            CALL SEAICE_CALC_STRESS(
                0502      I          e11, e22, e12, press, zeta, eta, etaZ,
                0503      O          sig1, sig2, sig12,
                0504      I          bi, bj, myTime, myIter, myThid )
                0505            DO j=1,sNy
                0506             DO i=1,sNx
                0507              sigp   = sig1(i,j) + sig2(i,j)
                0508              sigm   = sig1(i,j) - sig2(i,j)
                0509 C     This should be the way of computing sig12 at C-points,
                0510 C            sigTmp = 0.25 _d 0 *
                0511 C    &            ( sig12(i  ,j  ) + sig12(i+1,j  )
                0512 C    &            + sig12(i  ,j+1) + sig12(i+1,j+1) )
                0513 C     but sig12 = 2*etaZ*e12, and because of strong gradients in eta,
                0514 C     etaZ can be very large for a cell with small eta and the straightforward
                0515 C     way of averaging mixes large etaZ with small press0, so we have to do it
                0516 C     in different way to get meaningfull sig12C (=sigTmp):
                0517              sigTmp = 2. _d 0 * eta(i,j,bi,bj) * 0.25 _d 0 *
                0518      &            (e12(i,j,bi,bj) + e12(i+1,j,bi,bj)
                0519      &            +e12(i,j+1,bi,bj)+e12(i+1,j+1,bi,bj))
                0520              sigTmp = SQRT( sigm*sigm + 4. _d 0*sigTmp*sigTmp )
                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           ENDIF
                0529           CALL DIAGNOSTICS_FILL(sig1,'SIsig1  ',0,1,2,bi,bj,myThid)
                0530           CALL DIAGNOSTICS_FILL(sig2,'SIsig2  ',0,1,2,bi,bj,myThid)
                0531          ENDIF
                0532 C
                0533          IF ( diag_SIshear_isOn  ) THEN
                0534           DO j=1,sNy
                0535            DO i=1,sNx
                0536             sigm = e11(i,j,bi,bj) - e22(i,j,bi,bj)
                0537             sigTmp =
                0538      &           ( e12(i  ,j  ,bi,bj)**2 + e12(i+1,j  ,bi,bj)**2
                0539      &           + e12(i+1,j+1,bi,bj)**2 + e12(i  ,j+1,bi,bj)**2 )
                0540 C     shear deformation as sqrt((e11-e22)**2 + 4*e12**2); the 4 pulled into
                0541 C     the average
8377b8ee87 Mart*0542             sig1(i,j) = SQRT(sigm*sigm + sigTmp)
5fe78992ba Mart*0543            ENDDO
                0544           ENDDO
                0545           CALL DIAGNOSTICS_FILL(sig1,'SIshear ',0,1,2,bi,bj,myThid)
                0546          ENDIF
                0547 C
                0548 C     most of the energy diagnostics, re-use sig1, sig2, sig12 as tmp-arrays
                0549 C
                0550          IF ( diag_SIenpi_isOn ) THEN
                0551 C     compute internal stresses with updated ice velocities
                0552 C     TAF gets confused when we use stressDivergenceX/Y as temporary arrays
                0553 C     therefore we need to use new arrays
                0554           IF ( .NOT. SEAICEuseFREEDRIFT )
                0555      &         CALL SEAICE_CALC_STRESSDIV(
                0556      I         e11, e22, e12, press, zeta, eta, etaZ,
                0557 # ifdef ALLOW_AUTODIFF
                0558      O         strDivX, strDivY,
                0559 # else /* not ALLOW_AUTODIFF */
                0560      O         stressDivergenceX, stressDivergenceY,
                0561 # endif /* ALLOW_AUTODIFF */
                0562      I         bi, bj, myTime, myIter, myThid )
                0563           DO j=1,sNy+1
                0564            DO i=1,sNx+1
                0565 # ifdef ALLOW_AUTODIFF
                0566             sig1(i,j) = uIce(i,j,bi,bj)*strDivX(i,j,bi,bj)
                0567             sig2(i,j) = vIce(i,j,bi,bj)*strDivY(i,j,bi,bj)
                0568 # else /* not ALLOW_AUTODIFF */
                0569             sig1(i,j) = uIce(i,j,bi,bj)*stressDivergenceX(i,j,bi,bj)
                0570             sig2(i,j) = vIce(i,j,bi,bj)*stressDivergenceY(i,j,bi,bj)
                0571 # endif /* ALLOW_AUTODIFF */
                0572            ENDDO
                0573           ENDDO
                0574 C     average from velocity points to pressure points
                0575           DO j=1,sNy
                0576            DO i=1,sNx
                0577             sig12(i,j)= 0.5 _d 0 * ( sig1(i,j) + sig1(i+1,j)
                0578      &                             + sig2(i,j) + sig2(i,j+1) )
                0579            ENDDO
                0580           ENDDO
                0581           CALL DIAGNOSTICS_FILL(sig12,'SIenpi  ',0,1,2,bi,bj,myThid)
                0582          ENDIF
                0583          IF ( diag_SIenpot_isOn ) THEN
                0584           DO j=1,sNy
                0585            DO i=1,sNx
                0586             sig1(i,j) = 0.5 _d 0 * press0(i,j,bi,bj) *
                0587      &           ( e11(i,j,bi,bj) + e22(i,j,bi,bj) )
                0588            ENDDO
                0589           ENDDO
                0590           CALL DIAGNOSTICS_FILL(sig1,'SIenpot ',0,1,2,bi,bj,myThid)
                0591          ENDIF
                0592          IF ( diag_SIpRfric_isOn ) THEN
                0593           DO j=1,sNy
                0594            DO i=1,sNx
                0595             sig1(i,j) = - zeta(i,j,bi,bj)
                0596      &           * ( e11(i,j,bi,bj) + e22(i,j,bi,bj) )
                0597      &           * ( e11(i,j,bi,bj) + e22(i,j,bi,bj) )
                0598            ENDDO
                0599           ENDDO
                0600           CALL DIAGNOSTICS_FILL(sig1,'SIpRfric',0,1,2,bi,bj,myThid)
                0601          ENDIF
                0602          IF ( diag_SIpSfric_isOn ) THEN
                0603           DO j=1,sNy
                0604            DO i=1,sNx
                0605             sig1(i,j) = - eta(i,j,bi,bj)
                0606      &        * ( (e11(i,j,bi,bj) - e22(i,j,bi,bj))**2 +
                0607      &            ( e12(i  ,j  ,bi,bj)**2 + e12(i+1,j  ,bi,bj)**2
                0608      &            + e12(i+1,j+1,bi,bj)**2 + e12(i  ,j+1,bi,bj)**2 )
                0609      &          )
                0610            ENDDO
                0611           ENDDO
                0612           CALL DIAGNOSTICS_FILL(sig1,'SIpSfric',0,1,2,bi,bj,myThid)
                0613          ENDIF
                0614          IF ( DIAGNOSTICS_IS_ON('SIenpa  ',myThid) ) THEN
                0615           DO j=1,sNy+1
                0616            DO i=1,sNx+1
                0617             areaW = 0.5 _d 0 * (AREA(i,j,bi,bj) + AREA(i-1,j,bi,bj))
                0618      &           * SEAICEstressFactor
                0619             areaS = 0.5 _d 0 * (AREA(i,j,bi,bj) + AREA(i,j-1,bi,bj))
                0620      &           * SEAICEstressFactor
                0621             sig1(i,j) = TAUX(i,j,bi,bj)*areaW * uIce(i,j,bi,bj)
                0622             sig2(i,j) = TAUY(i,j,bi,bj)*areaS * vIce(i,j,bi,bj)
                0623            ENDDO
                0624           ENDDO
                0625 C     average from velocity points to pressure points
                0626           DO j=1,sNy
                0627            DO i=1,sNx
                0628             sig12(i,j)= 0.5 _d 0 * ( sig1(i,j) + sig1(i+1,j)
                0629      &                             + sig2(i,j) + sig2(i,j+1) )
                0630            ENDDO
                0631           ENDDO
                0632           CALL DIAGNOSTICS_FILL(sig12,'SIenpa  ',0,1,2,bi,bj,myThid)
                0633          ENDIF
                0634          IF ( DIAGNOSTICS_IS_ON('SIenpw  ',myThid) ) THEN
                0635 C     surface level
                0636           kSrf = 1
                0637 C     introduce turning angle (default is zero)
                0638           SINWAT=SIN(SEAICE_waterTurnAngle*deg2rad)
                0639           COSWAT=COS(SEAICE_waterTurnAngle*deg2rad)
                0640           DO j=1,sNy+1
                0641            DO i=1,sNx+1
                0642             areaW = 0.5 _d 0 * (AREA(i,j,bi,bj) + AREA(i-1,j,bi,bj))
                0643      &         * SEAICEstressFactor
                0644             areaS = 0.5 _d 0 * (AREA(i,j,bi,bj) + AREA(i,j-1,bi,bj))
                0645      &         * SEAICEstressFactor
                0646             sig1(i,j) = areaW *
                0647      &           HALF*( DWATN(i,j,bi,bj)+DWATN(i-1,j,bi,bj) )*
                0648      &         COSWAT *
                0649      &         ( uIce(i,j,bi,bj)-uVel(i,j,kSrf,bi,bj) )
                0650      &         - SIGN(SINWAT, _fCori(i,j,bi,bj)) * 0.5 _d 0 *
                0651      &         ( DWATN(i  ,j,bi,bj) *
                0652      &         0.5 _d 0*(vIce(i  ,j  ,bi,bj)-vVel(i  ,j  ,kSrf,bi,bj)
                0653      &                  +vIce(i  ,j+1,bi,bj)-vVel(i  ,j+1,kSrf,bi,bj))
                0654      &         + DWATN(i-1,j,bi,bj) *
                0655      &         0.5 _d 0*(vIce(i-1,j  ,bi,bj)-vVel(i-1,j  ,kSrf,bi,bj)
                0656      &                  +vIce(i-1,j+1,bi,bj)-vVel(i-1,j+1,kSrf,bi,bj))
                0657      &         )
                0658             sig1(i,j) = sig1(i,j) * uIce(i,j,bi,bj)
                0659             sig2(i,j) = areaS *
                0660      &           HALF*( DWATN(i,j,bi,bj)+DWATN(i,j-1,bi,bj) )*
                0661      &         COSWAT *
                0662      &         ( vIce(i,j,bi,bj)-vVel(i,j,kSrf,bi,bj) )
                0663      &         + SIGN(SINWAT,  _fCori(i,j,bi,bj)) * 0.5 _d 0 *
                0664      &         ( DWATN(i,j  ,bi,bj) *
                0665      &         0.5 _d 0*(uIce(i  ,j  ,bi,bj)-uVel(i  ,j  ,kSrf,bi,bj)
                0666      &                  +uIce(i+1,j  ,bi,bj)-uVel(i+1,j  ,kSrf,bi,bj))
                0667      &         + DWATN(i,j-1,bi,bj) *
                0668      &         0.5 _d 0*(uIce(i  ,j-1,bi,bj)-uVel(i  ,j-1,kSrf,bi,bj)
                0669      &                  +uIce(i+1,j-1,bi,bj)-uVel(i+1,j-1,kSrf,bi,bj))
                0670      &         )
                0671             sig2(i,j) = sig2(i,j) * vIce(i,j,bi,bj)
                0672            ENDDO
                0673           ENDDO
                0674 C     average from velocity points to pressure points
                0675           DO j=1,sNy
                0676            DO i=1,sNx
                0677             sig12(i,j)= - 0.5 _d 0 * ( sig1(i,j) + sig1(i+1,j)
                0678      &                               + sig2(i,j) + sig2(i,j+1) )
                0679            ENDDO
                0680           ENDDO
                0681           CALL DIAGNOSTICS_FILL(sig12,'SIenpw  ',0,1,2,bi,bj,myThid)
                0682          ENDIF
                0683          IF ( SEAICEuseTilt
                0684      &        .AND. DIAGNOSTICS_IS_ON('SIenpg  ',myThid) ) THEN
                0685           DO j=1-OLy+1,sNy+OLy
                0686            DO i=1-OLx+1,sNx+OLx
                0687             sig1(i,j)=
                0688      &           - seaiceMassU(i,j,bi,bj)*_recip_dxC(i,j,bi,bj)
                0689      &           *( phiSurf(i,j)-phiSurf(i-1,j) ) * uIce(i,j,bi,bj)
                0690             sig2(i,j)=
                0691      &           - seaiceMassV(i,j,bi,bj)* _recip_dyC(i,j,bi,bj)
                0692      &           *( phiSurf(i,j)-phiSurf(i,j-1) ) * vIce(i,j,bi,bj)
                0693            ENDDO
                0694           ENDDO
                0695 C     average from velocity points to pressure points
                0696           DO j=1,sNy
                0697            DO i=1,sNx
                0698             sig12(i,j) = 0.5 _d 0 * ( sig1(i,j) + sig1(i+1,j)
                0699      &                              + sig2(i,j) + sig2(i,j+1) )
                0700            ENDDO
                0701           ENDDO
                0702           CALL DIAGNOSTICS_FILL(sig12,'SIenpg  ',0,1,2,bi,bj,myThid)
                0703          ENDIF
                0704 C     bi/bj-loop
                0705         ENDDO
                0706        ENDDO
45315406aa Mart*0707 C     useDiagnostics & SEAICEuseDYNAMICS
5fe78992ba Mart*0708       ENDIF
45315406aa Mart*0709 #endif /* ALLOW_DIAGNOSTICS and SEAICE_CGRID */
5fe78992ba Mart*0710
53092bcb42 Mart*0711       RETURN
                0712       END