Back to home page

MITgcm
 
 

    

File indexing completed on 2023-08-04 05:10:41 UTC

view on githubraw file Latest commit 45315406 on 2023-08-03 16:50:12 UTC
809c36b928 Patr*0001 #include "SEAICE_OPTIONS.h"
772b2ed80e Gael*0002 #ifdef ALLOW_AUTODIFF
                0003 # include "AUTODIFF_OPTIONS.h"
                0004 #endif
cda1c18f72 Jean*0005 #ifndef SEAICE_LSRBNEW
                0006 C     Define this option below for an alternative discretization of
                0007 C     d/dx[ (zeta-eta) dV/dy] and d/dy[ (zeta-eta) dU/dx]
                0008 # undef SEAICE_TEST
                0009 #endif
809c36b928 Patr*0010 
                0011 CStartOfInterface
7109a141b2 Patr*0012       SUBROUTINE lsr( ilcall, myThid )
03c669d1ab Jean*0013 C     *==========================================================*
809c36b928 Patr*0014 C     | SUBROUTINE  lsr                                          |
09510da3bb Dimi*0015 C     | o Solve ice momentum equation with an LSR dynamics solver|
                0016 C     |   (see Zhang and Hibler,   JGR, 102, 8691-8702, 1997     |
                0017 C     |    and Zhang and Rothrock, MWR, 131,  845- 861, 2003)    |
                0018 C     |   Written by Jinlun Zhang, PSC/UW, Feb-2001              |
                0019 C     |                     zhang@apl.washington.edu             |
03c669d1ab Jean*0020 C     *==========================================================*
809c36b928 Patr*0021       IMPLICIT NONE
                0022 
                0023 C     === Global variables ===
                0024 #include "SIZE.h"
                0025 #include "EEPARAMS.h"
                0026 #include "PARAMS.h"
4366d31d92 Mart*0027 #include "GRID.h"
03c669d1ab Jean*0028 #include "SEAICE_SIZE.h"
809c36b928 Patr*0029 #include "SEAICE_PARAMS.h"
03c669d1ab Jean*0030 #include "SEAICE.h"
4366d31d92 Mart*0031 C#include "SEAICE_GRID.h"
809c36b928 Patr*0032 
7109a141b2 Patr*0033 #ifdef ALLOW_AUTODIFF_TAMC
                0034 # include "tamc.h"
                0035 #endif
809c36b928 Patr*0036 
                0037 C     === Routine arguments ===
                0038 C     myThid - Thread no. that called this routine.
7109a141b2 Patr*0039       INTEGER ilcall
809c36b928 Patr*0040       INTEGER myThid
                0041 CEndOfInterface
                0042 
45315406aa Mart*0043 #ifdef SEAICE_BGRID_DYNAMICS
809c36b928 Patr*0044 
                0045 C     === Local variables ===
cee16b76ae Dimi*0046 C     i,j,bi,bj - Loop counters
809c36b928 Patr*0047 
09510da3bb Dimi*0048       INTEGER i, j, m, bi, bj, j1, j2, im, jm
07b4a12853 Mart*0049       INTEGER ICOUNT1, ICOUNT2
e1d2983c5b Mart*0050 
                0051       _RL WFAU, WFAV, WFAU1, WFAV1, WFAU2, WFAV2
                0052       _RL AA3, S1, S2, S1A, S2A
                0053 
                0054       _RL e11loc, e22loc, e12loc
                0055 
                0056       _RL COSWAT
5acccad966 Jean*0057       _RS SINWAT
e1d2983c5b Mart*0058       _RL ECM2, DELT1, DELT2
                0059       _RL TEMPVAR
7109a141b2 Patr*0060 
e1d2983c5b Mart*0061 C     diagonals of coefficient matrices
                0062       _RL AU   (1:sNx,1:sNy,nSx,nSy)
                0063       _RL BU   (1:sNx,1:sNy,nSx,nSy)
                0064       _RL CU   (1:sNx,1:sNy,nSx,nSy)
                0065       _RL AV   (1:sNx,1:sNy,nSx,nSy)
                0066       _RL BV   (1:sNx,1:sNy,nSx,nSy)
                0067       _RL CV   (1:sNx,1:sNy,nSx,nSy)
809c36b928 Patr*0068 
e1d2983c5b Mart*0069 C     coefficients for lateral points, u(j+/-1)
                0070       _RL uRt1(1:sNx,1:sNy,nSx,nSy)
                0071       _RL uRt2(1:sNx,1:sNy,nSx,nSy)
                0072 C     coefficients for lateral points, v(i+/-1)
                0073       _RL vRt1(1:sNx,1:sNy,nSx,nSy)
                0074       _RL vRt2(1:sNx,1:sNy,nSx,nSy)
                0075 C     RHS
                0076       _RL rhsU (1:sNx,1:sNy,nSx,nSy)
                0077       _RL rhsV (1:sNx,1:sNy,nSx,nSy)
3f31c7d5de Mart*0078 C     symmetric and anti-symmetric drag coefficients
                0079       _RL DRAGS      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0080       _RL DRAGA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
e1d2983c5b Mart*0081 #ifdef SEAICE_LSRBNEW
d24acf3cfc Jean*0082       LOGICAL doIterate4u, doIterate4v
                0083       CHARACTER*(MAX_LEN_MBUF) msgBuf
e1d2983c5b Mart*0084 C     coefficients of ice velocities in coefficient matrix
                0085 C     for both U and V-equation
                0086 C     XX: double derivative in X
                0087 C     YY: double derivative in Y
                0088 C     XM: metric term with derivative in X
                0089 C     YM: metric term with derivative in Y
d24acf3cfc Jean*0090       _RL UXX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0091       _RL UYY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0092       _RL UXM  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0093       _RL UYM  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0094       _RL VXX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0095       _RL VYY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0096       _RL VXM  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0097       _RL VYM  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
e1d2983c5b Mart*0098 C     abbreviations
d24acf3cfc Jean*0099       _RL etaU (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0100       _RL zetaU(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0101       _RL etaV (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0102       _RL zetaV(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
e1d2983c5b Mart*0103 C     contribution of sigma on righ hand side
d24acf3cfc Jean*0104       _RL sig11(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0105       _RL sig22(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0106       _RL sig12(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0107       _RL sig21(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
4366d31d92 Mart*0108 
e1d2983c5b Mart*0109 C     auxillary fields
d24acf3cfc Jean*0110       _RL CUU  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0111       _RL CVV  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0112       _RL URT  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0113       _RL VRT  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
e1d2983c5b Mart*0114 #else
d24acf3cfc Jean*0115       _RL URT(1-OLx:sNx+OLx), CUU(1-OLx:sNx+OLx)
                0116       _RL VRT(1-OLy:sNy+OLy), CVV(1-OLy:sNy+OLy)
                0117 
                0118       _RL etaPlusZeta (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0119       _RL zetaMinusEta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0120       _RL ETAMEAN  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0121       _RL ZETAMEAN (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0122 
                0123       _RL dVdx     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0124       _RL dVdy     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0125       _RL dUdx     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0126       _RL dUdy     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
bb4175c1b5 Mart*0127 #ifdef SEAICE_TEST
d24acf3cfc Jean*0128       _RL uz     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0129       _RL vz     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
bb4175c1b5 Mart*0130 #endif
09510da3bb Dimi*0131 
e1d2983c5b Mart*0132       INTEGER phexit
                0133 
                0134       _RL  AA1, AA2, AA4, AA5, AA6
                0135 
                0136 #endif /* SEAICE_LSRBNEW */
                0137 
                0138       _RL UERR
                0139 
                0140 C     abbreviations
d24acf3cfc Jean*0141       _RL ucLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0142       _RL vcLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0143       _RL uTmp (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0144       _RL vTmp (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
e1d2983c5b Mart*0145 
809c36b928 Patr*0146 C SET SOME VALUES
09510da3bb Dimi*0147       WFAU1=0.95 _d 0
                0148       WFAV1=0.95 _d 0
                0149       WFAU2=ZERO
                0150       WFAV2=ZERO
                0151 
                0152       S1A=0.80 _d 0
                0153       S2A=0.80 _d 0
                0154       WFAU=WFAU1
                0155       WFAV=WFAV1
                0156 
79df32c3f1 Mart*0157       ICOUNT1=SEAICElinearIterMax
                0158       ICOUNT2=SEAICElinearIterMax
809c36b928 Patr*0159 
e1d2983c5b Mart*0160       ECM2= 1. _d 0/(SEAICE_eccen**2)
                0161 
                0162       SINWAT=SIN(SEAICE_waterTurnAngle*deg2rad)
                0163       COSWAT=COS(SEAICE_waterTurnAngle*deg2rad)
                0164 
                0165       DO bj=myByLo(myThid),myByHi(myThid)
                0166        DO bi=myBxLo(myThid),myBxHi(myThid)
d24acf3cfc Jean*0167         DO j=1-OLy,sNy+OLy-1
                0168          DO i=1-OLx,sNx+OLx-1
e1d2983c5b Mart*0169           ucLoc(I,J,bi,bj) = 0.25 _d 0 * (
45315406aa Mart*0170      &         +  uIceB(I+1,J,bi,bj) +  uIceB(I+1,J+1,bi,bj)
                0171      &         +  uIceB(I  ,J,bi,bj) +  uIceB(I,  J+1,bi,bj) )
e1d2983c5b Mart*0172           vcLoc(I,J,bi,bj) = 0.25 _d 0 * (
45315406aa Mart*0173      &         +  vIceB(I+1,J,bi,bj) +  vIceB(I+1,J+1,bi,bj)
                0174      &         +  vIceB(I  ,J,bi,bj) +  vIceB(I,  J+1,bi,bj) )
e1d2983c5b Mart*0175          ENDDO
                0176         ENDDO
d24acf3cfc Jean*0177         DO J=1-OLy,sNy+OLy-1
                0178          DO I=1-OLy,sNx+OLy-1
e1d2983c5b Mart*0179           e11loc = 0.5 _d 0 * _recip_dxF(I,J,bi,bj) *
45315406aa Mart*0180      &                (uIceB(I+1,J+1,bi,bj)+uIceB(I+1,J,bi,bj)
                0181      &                -uIceB(I,  J+1,bi,bj)-uIceB(I,  J,bi,bj))
e1d2983c5b Mart*0182      &                + vcLoc(I,J,bi,bj) * k2AtC(I,J,bi,bj)
                0183           e22loc = 0.5 _d 0 * _recip_dyF(I,J,bi,bj) *
45315406aa Mart*0184      &                (vIceB(I+1,J+1,bi,bj)+vIceB(I,J+1,bi,bj)
                0185      &                -vIceB(I+1,J,  bi,bj)-vIceB(I,J,  bi,bj))
e1d2983c5b Mart*0186      &                + ucLoc(I,J,bi,bj) * k1AtC(I,J,bi,bj)
d24acf3cfc Jean*0187           e12loc = 0.5 _d 0*(
e1d2983c5b Mart*0188      &           0.5 _d 0 * _recip_dyF(I,J,bi,bj) *
45315406aa Mart*0189      &         (uIceB(I+1,J+1,bi,bj)+uIceB(I,J+1,bi,bj)
                0190      &         -uIceB(I+1,J,  bi,bj)-uIceB(I,J,  bi,bj))
e1d2983c5b Mart*0191      &         + 0.5 _d 0 * _recip_dxF(I,J,bi,bj) *
45315406aa Mart*0192      &         (vIceB(I+1,J+1,bi,bj)+vIceB(I+1,J,bi,bj)
                0193      &         -vIceB(I,  J+1,bi,bj)-vIceB(I,  J,bi,bj))
e1d2983c5b Mart*0194      &         - vcLoc(I,J,bi,bj)*k1AtC(I,J,bi,bj)
                0195      &         - ucLoc(I,J,bi,bj)*k2AtC(I,J,bi,bj) )
                0196 C  NOW EVALUATE VISCOSITIES
                0197           DELT1=(e11loc**2+e22loc**2)*(1. _d 0+ECM2)
                0198      &         +4.0 _d 0*ECM2*e12loc**2
                0199      &         +2.0 _d 0*e11loc*e22loc*(1. _d 0-ECM2)
                0200           IF ( DELT1 .LE. SEAICE_EPS_SQ ) THEN
                0201              DELT2=SEAICE_EPS
                0202           ELSE
                0203              DELT2=SQRT(DELT1)
                0204           ENDIF
                0205           ZETA(I,J,bi,bj)  = 0.5 _d 0*PRESS0(I,J,bi,bj)/DELT2
                0206 C NOW PUT MIN AND MAX VISCOSITIES IN
8e32c48b8f Mart*0207           ZETA(I,J,bi,bj)  = MIN(SEAICE_zMax(I,J,bi,bj),ZETA(I,J,bi,bj))
                0208           ZETA(I,J,bi,bj)  = MAX(SEAICE_zMin(I,J,bi,bj),ZETA(I,J,bi,bj))
e1d2983c5b Mart*0209 C NOW SET VISCOSITIES TO ZERO AT HEFFMFLOW PTS
                0210           ZETA(I,J,bi,bj)  = ZETA(I,J,bi,bj)*HEFFM(I,J,bi,bj)
                0211           ETA(I,J,bi,bj)   = ECM2*ZETA(I,J,bi,bj)
                0212           PRESS(I,J,bi,bj) = 2.0 _d 0*ZETA(I,J,bi,bj)*DELT2
                0213          ENDDO
                0214         ENDDO
                0215         DO j=1,sNy
                0216          DO i=1,sNx
                0217 C NOW SET UP NON-LINEAR WATER DRAG, FORCEX, FORCEY
                0218           TEMPVAR=(uIce(I,J,bi,bj)-GWATX(I,J,bi,bj))**2
                0219      &           +(vIce(I,J,bi,bj)-GWATY(I,J,bi,bj))**2
                0220           IF ( YC(I,J,bi,bj) .LT. ZERO ) THEN
b8665dacca Mart*0221            DWATN(I,J,bi,bj) = SEAICE_waterDrag_south*SQRT(TEMPVAR)
e1d2983c5b Mart*0222           ELSE
b8665dacca Mart*0223            DWATN(I,J,bi,bj) = SEAICE_waterDrag*SQRT(TEMPVAR)
e1d2983c5b Mart*0224           ENDIF
b8665dacca Mart*0225           DWATN(I,J,bi,bj) = MAX(DWATN(I,J,bi,bj)*rhoConst,QUART)
e1d2983c5b Mart*0226 C NOW SET UP SYMMETTRIC DRAG
                0227           DRAGS(I,J,bi,bj)=DWATN(I,J,bi,bj)*COSWAT
                0228 C NOW SET UP ANTI SYMMETTRIC DRAG PLUS CORIOLIS
                0229           DRAGA(I,J,bi,bj)=DWATN(I,J,bi,bj)
                0230      &         *SIGN(SINWAT, _fCori(I,J,bi,bj))
                0231      &         + AMASS(I,J,bi,bj) * _fCoriG(I,J,bi,bj)
                0232 C NOW ADD IN CURRENT FORCE
                0233           FORCEX(I,J,bi,bj)=FORCEX0(I,J,bi,bj)+DWATN(I,J,bi,bj)
                0234      &         *(COSWAT*GWATX(I,J,bi,bj)
                0235      &         -SIGN(SINWAT, _fCori(I,J,bi,bj))*GWATY(I,J,bi,bj))
                0236           FORCEY(I,J,bi,bj)=FORCEY0(I,J,bi,bj)+DWATN(I,J,bi,bj)
                0237      &         *(SIGN(SINWAT, _fCori(I,J,bi,bj))*GWATX(I,J,bi,bj)
                0238      &         +COSWAT*GWATY(I,J,bi,bj))
                0239 #ifndef SEAICE_LSRBNEW
                0240 C NOW CALCULATE PRESSURE FORCE AND ADD TO EXTERNAL FORCE
                0241 C     only for old code, in the new code the pressure force
                0242 C     is added to the rhs stress tensor terms
                0243           FORCEX(I,J,bi,bj)=FORCEX(I,J,bi,bj)
                0244      &          -QUART * _recip_dxV(I,J,bi,bj)
                0245      &          *(PRESS(I,  J,bi,bj) + PRESS(I,  J-1,bi,bj)
                0246      &          - PRESS(I-1,J,bi,bj) - PRESS(I-1,J-1,bi,bj))
                0247           FORCEY(I,J,bi,bj)=FORCEY(I,J,bi,bj)
                0248      &         -QUART * _recip_dyU(I,J,bi,bj)
                0249      &          *(PRESS(I,J,  bi,bj) + PRESS(I-1,J,  bi,bj)
                0250      &          - PRESS(I,J-1,bi,bj) - PRESS(I-1,J-1,bi,bj))
                0251 #endif /* not SEAICE_LSRBNEW */
                0252           FORCEX(I,J,bi,bj)=FORCEX(I,J,bi,bj)
                0253      &         +AMASS(I,J,bi,bj)/SEAICE_deltaTdyn*uIceNm1(I,J,bi,bj)
                0254           FORCEY(I,J,bi,bj)=FORCEY(I,J,bi,bj)
                0255      &         +AMASS(I,J,bi,bj)/SEAICE_deltaTdyn*vIceNm1(I,J,bi,bj)
                0256           FORCEX(I,J,bi,bj)=FORCEX(I,J,bi,bj)*UVM(I,J,bi,bj)
                0257           FORCEY(I,J,bi,bj)=FORCEY(I,J,bi,bj)*UVM(I,J,bi,bj)
                0258          ENDDO
                0259         ENDDO
                0260        ENDDO
                0261       ENDDO
                0262 
                0263 #ifdef SEAICE_LSRBNEW
                0264       DO bj=myByLo(myThid),myByHi(myThid)
                0265        DO bi=myBxLo(myThid),myBxHi(myThid)
                0266 C     coefficients for matrices
d24acf3cfc Jean*0267         DO j=1-OLy,sNy+OLy-1
                0268          DO i=1-OLx+1,sNx+OLx-1
e1d2983c5b Mart*0269           etaU(I,J) = 0.5 _d 0 * (
                0270      &         +  eta(I,J,bi,bj) +  eta(I-1,J,bi,bj) )
                0271           zetaU(I,J)= 0.5 _d 0 *(
                0272      &         + zeta(I,J,bi,bj) + zeta(I-1,J,bi,bj) )
                0273          ENDDO
                0274         ENDDO
d24acf3cfc Jean*0275         DO j=1-OLy+1,sNy+OLy-1
                0276          DO i=1-OLx,sNx+OLx-1
e1d2983c5b Mart*0277           etaV(I,J) = 0.5 _d 0 * (
                0278      &         +  eta(I,J,bi,bj) +  eta(I,J-1,bi,bj) )
                0279           zetaV(I,J)= 0.5 _d 0 *(
                0280      &         + zeta(I,J,bi,bj) + zeta(I,J-1,bi,bj) )
                0281          ENDDO
                0282         ENDDO
                0283 C     coefficients of uIce(I,J) and vIce(I,J) belonging to ...
                0284         DO J=1,sNy
                0285          DO I=0,sNx
                0286 C     ... d/dx (eta+zeta) d/dx u
                0287           UXX(I,J) = _dyC(I,J,bi,bj) * (zetaV(I,J)+etaV(I,J))
                0288      &         * _recip_dxG(I,J,bi,bj)
                0289 C     ... d/dx (zeta-eta) k1 u
                0290           UXM(I,J) = _dyC(I,J,bi,bj) * (zetaV(I,J)-etaV(I,J))
                0291      &         * k1AtV(I,J,bi,bj) * 0.5 _d 0
                0292          ENDDO
                0293         ENDDO
                0294         DO J=0,sNy
                0295          DO I=1,sNx
                0296 C     ... d/dy eta d/dy u
                0297           UYY(I,J) = _dxC(I,J,bi,bj) * etaU(I,J)
                0298      &         * _recip_dyG(I,J,bi,bj)
                0299 C     ... d/dy eta k2 u
                0300           UYM(I,J) = _dxC(I,J,bi,bj) * etaU(I,J)
                0301      &         * k2AtU(I,J,bi,bj) * 0.5 _d 0
                0302          ENDDO
                0303         ENDDO
                0304         DO J=1,sNy
                0305          DO I=0,sNx
                0306 C     ... d/dx eta dv/dx
                0307           VXX(I,J) = _dyC(I,J,bi,bj) * etaV(I,J)
                0308      &         * _recip_dxG(I,J,bi,bj)
                0309 C     ... d/dx eta k1 v
                0310           VXM(I,J) = _dyC(I,J,bi,bj) * etaV(I,J)
                0311      &         * k1AtV(I,J,bi,bj) * 0.5 _d 0
                0312          ENDDO
                0313         ENDDO
                0314         DO J=0,sNy
                0315          DO I=1,sNx
                0316 C     ... d/dy eta+zeta dv/dy
                0317           VYY(I,J) = _dxC(I,J,bi,bj) * (zetaU(I,J)+etaU(I,J))
                0318      &         * _recip_dyG(I,J,bi,bj)
                0319 C     ... d/dy (zeta-eta) k2 v
                0320           VYM(I,J) = _dxC(I,J,bi,bj) * (zetaU(I,J)-etaU(I,J))
                0321      &         * k2AtU(I,J,bi,bj) * 0.5 _d 0
                0322          ENDDO
                0323         ENDDO
                0324 
                0325 C     assemble coefficient matrix of uIce
                0326 C     beware of sign convention: because this
                0327 C     is the left hand side we calculate -grad(sigma), but the coefficients
                0328 C     of U(I,J+/-1) are counted on the right hand side
                0329         DO J=1,sNy
                0330          DO I=1,sNx
                0331 C     coefficients for uIce(I-1,J)
                0332           AU(I,J,bi,bj)= ( - UXX(I-1,J) + UXM(I-1,J) )
                0333      &         * UVM(I,J,bi,bj)
                0334 C     coefficients for uIce(I+1,J)
                0335           CU(I,J,bi,bj)= ( - UXX(I  ,J) - UXM(I  ,J) )
                0336      &         * UVM(I,J,bi,bj)
                0337 C     coefficients for uIce(I,J)
                0338           BU(I,J,bi,bj)=(ONE - UVM(I,J,bi,bj)) +
                0339      &         ( UXX(I-1,J) + UXX(I,J) + UYY(I,J-1) + UYY(I,J)
                0340      &         + UXM(I-1,J) - UXM(I,J) - UYM(I,J-1) + UYM(I,J)
                0341      &         ) * UVM(I,J,bi,bj)
                0342 C     coefficients of uIce(I,J-1)
                0343           uRt1(I,J,bi,bj)= UYY(I,J-1) + UYM(I,J-1)
                0344 C     coefficients of uIce(I,J+1)
                0345           uRt2(I,J,bi,bj)= UYY(I,J  ) - UYM(I,J  )
                0346          ENDDO
                0347         ENDDO
                0348 
                0349 C     now we need to normalize everything by the grid cell area
                0350         DO J=1,sNy
                0351          DO I=1,sNx
                0352           AU(I,J,bi,bj)    = AU(I,J,bi,bj)    * recip_rAz(I,J,bi,bj)
                0353           CU(I,J,bi,bj)    = CU(I,J,bi,bj)    * recip_rAz(I,J,bi,bj)
                0354 C     here we need add in the contribution from the time derivative
                0355 C     and the symmetric drag term; must be done after normalizing
                0356           BU(I,J,bi,bj)    = BU(I,J,bi,bj)    * recip_rAz(I,J,bi,bj)
                0357      &         + UVM(I,J,bi,bj) *
                0358      &         ( AMASS(I,J,bi,bj)/SEAICE_deltaTdyn
                0359      &         + DRAGS(I,J,bi,bj) )
                0360           uRt1(I,J,bi,bj) = uRt1(I,J,bi,bj) * recip_rAz(I,J,bi,bj)
                0361           uRt2(I,J,bi,bj) = uRt2(I,J,bi,bj) * recip_rAz(I,J,bi,bj)
                0362          ENDDO
                0363         ENDDO
                0364 
                0365 C     assemble coefficient matrix of vIce
                0366 C     beware of sign convention: because this
                0367 C     is the left hand side we calculate -grad(sigma), but the coefficients
                0368 C     of V(I,J+/-1) are counted on the right hand side
                0369         DO J=1,sNy
                0370          DO I=1,sNx
                0371 C     coefficients for vIce(I,J-1)
                0372           AV(I,J,bi,bj)=( - VYY(I,J-1) + VYM(I,J-1)
                0373      &         ) * UVM(I,J,bi,bj)
                0374 C     coefficients for vIce(I,J+1)
                0375           CV(I,J,bi,bj)=( - VYY(I,J  ) - VYM(I,J  )
                0376      &         ) * UVM(I,J,bi,bj)
                0377 C     coefficients for vIce(I,J)
                0378           BV(I,J,bi,bj)= (ONE - UVM(I,J,bi,bj)) +
                0379      &         ( VXX(I,J) + VXX(I-1,J) + VYY(I,J) + VYY(I,J-1)
                0380      &         + VXM(I,J) - VXM(I-1,J) - VYM(I,J) + VYM(I,J-1)
                0381      &         ) * UVM(I,J,bi,bj)
                0382 C     coefficients for V(I-1,J)
                0383           vRt1(I,J,bi,bj) = VXX(I-1,J) + VXM(I-1,J)
                0384 C     coefficients for V(I+1,J)
                0385           vRt2(I,J,bi,bj) = VXX(I  ,J) - VXM(I  ,J)
                0386          ENDDO
                0387         ENDDO
                0388 
                0389 C     now we need to normalize everything by the grid cell area
                0390         DO J=1,sNy
                0391          DO I=1,sNx
                0392           AV(I,J,bi,bj)    = AV(I,J,bi,bj)    * recip_rAz(I,J,bi,bj)
                0393           CV(I,J,bi,bj)    = CV(I,J,bi,bj)    * recip_rAz(I,J,bi,bj)
                0394 C     here we need add in the contribution from the time derivative
                0395 C     and the symmetric drag term; must be done after normalizing
                0396           BV(I,J,bi,bj)    = BV(I,J,bi,bj)    * recip_rAz(I,J,bi,bj)
                0397      &         + UVM(I,J,bi,bj) *
                0398      &         ( AMASS(I,J,bi,bj)/SEAICE_deltaTdyn
                0399      &         + DRAGS(I,J,bi,bj) )
                0400           vRt1(I,J,bi,bj) = vRt1(I,J,bi,bj) * recip_rAz(I,J,bi,bj)
                0401           vRt2(I,J,bi,bj) = vRt2(I,J,bi,bj) * recip_rAz(I,J,bi,bj)
                0402          ENDDO
                0403         ENDDO
                0404 
                0405 C     set up right-hand sides
                0406         DO j=1,sNy
                0407          DO i=0,sNx
                0408 C     at V-point
d24acf3cfc Jean*0409           sig11(I,J) =
                0410      &           (zetaV(I,J)-etaV(I,J))
                0411      &         * (vcLoc(I,J,bi,bj)-vcLoc(I,J-1,bi,bj))
e1d2983c5b Mart*0412      &         * _recip_dyC(I,J,bi,bj)
                0413      &         + (zetaV(I,J)+etaV(I,J))*k2atV(I,J,bi,bj)
45315406aa Mart*0414      &         * 0.5 _d 0 * (vIceB(I,J,bi,bj)+vIceB(I+1,J,bi,bj))
e1d2983c5b Mart*0415      &         - 0.5 _d 0 * (PRESS(I,J,bi,bj)+PRESS(I,J-1,bi,bj))
d24acf3cfc Jean*0416           sig12(I,J) = etaV(I,J)
e1d2983c5b Mart*0417      &         * (ucLoc(I,J,bi,bj)-ucLoc(I,J-1,bi,bj))
                0418      &         * _recip_dyC(I,J,bi,bj)
                0419      &         - etaV(I,J) * k2AtV(I,J,bi,bj)
45315406aa Mart*0420      &         * 0.5 _d 0 * (uIceB(I,J,bi,bj)+uIceB(I+1,J,bi,bj))
e1d2983c5b Mart*0421          ENDDO
                0422         ENDDO
                0423         DO j=0,sNy
                0424          DO i=1,sNx
                0425 C     at U-point
d24acf3cfc Jean*0426           sig22(I,J) =
                0427      &           (zetaU(I,J)-etaU(I,J))
                0428      &         * (ucLoc(I,J,bi,bj)-ucLoc(I-1,J,bi,bj))
e1d2983c5b Mart*0429      &         * _recip_dxC(I,J,bi,bj)
                0430      &         + (zetaU(I,J)+etaU(I,J))*k2atU(I,J,bi,bj)
45315406aa Mart*0431      &         * 0.5 _d 0 * (uIceB(I,J,bi,bj)+uIceB(I,J+1,bi,bj))
e1d2983c5b Mart*0432      &         - 0.5 _d 0 * (PRESS(I,J,bi,bj)+PRESS(I-1,J,bi,bj))
d24acf3cfc Jean*0433           sig21(I,J) = etaU(I,J)
e1d2983c5b Mart*0434      &         * (vcLoc(I,J,bi,bj)-vcLoc(I-1,J,bi,bj))
                0435      &         * _recip_dxC(I,J,bi,bj)
                0436      &         - etaU(I,J) * k1AtU(I,J,bi,bj)
45315406aa Mart*0437      &         * 0.5 _d 0 * (vIceB(I,J,bi,bj)+vIceB(I,J+1,bi,bj))
e1d2983c5b Mart*0438          ENDDO
                0439         ENDDO
                0440         DO j=1,sNy
                0441          DO i=1,sNx
45315406aa Mart*0442           rhsU(I,J,bi,bj) =   DRAGA(I,J,bi,bj)*vIceB(I,J,bi,bj)
e1d2983c5b Mart*0443      &         +FORCEX(I,J,bi,bj)
                0444      &         + ( _dyC(I,  J,  bi,bj) * sig11(I,  J  )
                0445      &           - _dyC(I-1,J,  bi,bj) * sig11(I-1,J  )
                0446      &           + _dxC(I,  J,  bi,bj) * sig21(I,  J  )
                0447      &           - _dxC(I,  J-1,bi,bj) * sig21(I,  J-1) )
                0448      &         * recip_rAz(I,J,bi,bj) * UVM(I,J,bi,bj)
45315406aa Mart*0449           rhsV(I,J,bi,bj) = - DRAGA(I,J,bi,bj)*uIceB(I,J,bi,bj)
e1d2983c5b Mart*0450      &         +FORCEY(I,J,bi,bj)
d24acf3cfc Jean*0451      &         + ( _dyC(I,  J,  bi,bj) * sig12(I,  J  )
                0452      &           - _dyC(I-1,J,  bi,bj) * sig12(I-1,J  )
                0453      &           + _dxC(I,  J,  bi,bj) * sig22(I,  J  )
e1d2983c5b Mart*0454      &           - _dxC(I,  J-1,bi,bj) * sig22(I,  J-1) )
                0455      &         * recip_rAz(I,J,bi,bj) * UVM(I,J,bi,bj)
                0456          ENDDO
                0457         ENDDO
                0458 C     bi/bj-loops
                0459        ENDDO
                0460       ENDDO
                0461 
                0462 C NOW DO ITERATION
                0463 
                0464       doIterate4u = .TRUE.
                0465       doIterate4v = .TRUE.
                0466 
                0467 C ITERATION START -----------------------------------------------------
                0468 
79df32c3f1 Mart*0469       DO m = 1, SEAICElinearIterMax
e1d2983c5b Mart*0470        IF ( doIterate4u .OR. doIterate4v ) THEN
                0471 
                0472         IF ( useCubedSphereExchange ) THEN
                0473           doIterate4u = .TRUE.
                0474           doIterate4v = .TRUE.
                0475         ENDIF
                0476 
                0477         DO bj=myByLo(myThid),myByHi(myThid)
                0478          DO bi=myBxLo(myThid),myBxHi(myThid)
                0479 
                0480 C-jmc: get less TAF warnings when always (no if doIterate) saving uIce,vIce:
                0481 C     save uIce prior to iteration, NOW SET U(3)=U(1)
                0482           DO j=1-OLy,sNy+OLy
                0483            DO i=1-OLx,sNx+OLx
                0484             uTmp(I,J,bi,bj)=uIce(I,J,bi,bj)
                0485            ENDDO
                0486           ENDDO
                0487 C     save vIce prior to iteration, NOW SET V(3)=V(1)
                0488           DO j=1-OLy,sNy+OLy
                0489            DO i=1-OLx,sNx+OLx
                0490             vTmp(I,J,bi,bj)=vIce(I,J,bi,bj)
                0491            ENDDO
                0492           ENDDO
                0493 
                0494           IF ( doIterate4u ) THEN
                0495 C Solve for uIce :
                0496            DO J=1,sNy
                0497             DO I=1,sNx
                0498              AA3 = ZERO
                0499              IF (I.EQ.1)   AA3 = AA3 - AU(I,J,bi,bj)*uIce(I-1,J,bi,bj)
                0500              IF (I.EQ.sNx) AA3 = AA3 - CU(I,J,bi,bj)*uIce(I+1,J,bi,bj)
                0501 
                0502              URT(I,J)=rhsU(I,J,bi,bj)
                0503      &            + AA3
                0504 #ifdef SEAICE_VECTORIZE_LSR
                0505      &            + uRt1(I,J,bi,bj)*uTmp(I,J-1,bi,bj)
                0506      &            + uRt2(I,J,bi,bj)*uTmp(I,J+1,bi,bj)
                0507 #else
                0508      &            + uRt1(I,J,bi,bj)*uIce(I,J-1,bi,bj)
                0509      &            + uRt2(I,J,bi,bj)*uIce(I,J+1,bi,bj)
                0510 #endif /* SEAICE_VECTORIZE_LSR */
                0511              URT(I,J)=URT(I,J)* UVM(I,J,bi,bj)
                0512             ENDDO
                0513 
                0514             DO I=1,sNx
                0515              CUU(I,J)=CU(I,J,bi,bj)
                0516             ENDDO
                0517             CUU(1,J)=CUU(1,J)/BU(1,J,bi,bj)
                0518             URT(1,J)=URT(1,J)/BU(1,J,bi,bj)
                0519 #ifdef SEAICE_VECTORIZE_LSR
                0520            ENDDO
                0521 C     start a new loop with reversed order to support automatic vectorization
                0522            DO I=2,sNx
                0523             IM=I-1
                0524             DO J=1,sNy
                0525 #else /* do not SEAICE_VECTORIZE_LSR */
                0526             DO I=2,sNx
                0527              IM=I-1
                0528 #endif /* SEAICE_VECTORIZE_LSR */
                0529              CUU(I,J)=CUU(I,J)/(BU(I,J,bi,bj)-AU(I,J,bi,bj)*CUU(IM,J))
                0530              URT(I,J)=(URT(I,J)-AU(I,J,bi,bj)*URT(IM,J))
                0531      &           /(BU(I,J,bi,bj)-AU(I,J,bi,bj)*CUU(IM,J))
                0532             ENDDO
                0533 #ifdef SEAICE_VECTORIZE_LSR
                0534            ENDDO
                0535 C     go back to original order
                0536            DO J=1,sNy
                0537 #endif /* SEAICE_VECTORIZE_LSR */
                0538             DO I=1,sNx-1
                0539              J1=sNx-I
                0540              J2=J1+1
                0541              URT(J1,J)=URT(J1,J)-CUU(J1,J)*URT(J2,J)
                0542             ENDDO
                0543             DO I=1,sNx
                0544              uIce(I,J,bi,bj)=uTmp(I,J,bi,bj)
                0545      &           +WFAU*(URT(I,J)-uTmp(I,J,bi,bj))
                0546             ENDDO
                0547            ENDDO
                0548 C--   end doIterate4u
                0549           ENDIF
                0550 
                0551           IF ( doIterate4v ) THEN
                0552 C Solve for vIce
                0553            DO I=1,sNx
                0554             DO J=1,sNy
                0555              AA3 = ZERO
                0556              IF (J.EQ.1)   AA3 = AA3 - AV(I,J,bi,bj)*vIce(I,J-1,bi,bj)
                0557              IF (J.EQ.sNy) AA3 = AA3 - CV(I,J,bi,bj)*vIce(I,J+1,bi,bj)
                0558 
                0559              VRT(I,J)=rhsV(I,J,bi,bj)
                0560      &            + AA3
                0561 #ifdef SEAICE_VECTORIZE_LSR
                0562      &            + vRt1(I,J,bi,bj)*vTmp(I-1,J,bi,bj)
                0563      &            + vRt2(I,J,bi,bj)*vTmp(I+1,J,bi,bj)
                0564 #else
                0565      &            + vRt1(I,J,bi,bj)*vIce(I-1,J,bi,bj)
                0566      &            + vRt2(I,J,bi,bj)*vIce(I+1,J,bi,bj)
                0567 #endif /* SEAICE_VECTORIZE_LSR */
                0568              VRT(I,J)=VRT(I,J)* UVM(I,J,bi,bj)
                0569             ENDDO
                0570 
                0571             DO J=1,sNy
                0572              CVV(I,J)=CV(I,J,bi,bj)
                0573             ENDDO
                0574             CVV(I,1)=CVV(I,1)/BV(I,1,bi,bj)
                0575             VRT(I,1)=VRT(I,1)/BV(I,1,bi,bj)
                0576             DO J=2,sNy
                0577              JM=J-1
                0578              CVV(I,J)=CVV(I,J)/(BV(I,J,bi,bj)-AV(I,J,bi,bj)*CVV(I,JM))
                0579              VRT(I,J)=(VRT(I,J)-AV(I,J,bi,bj)*VRT(I,JM))
                0580      &            /(BV(I,J,bi,bj)-AV(I,J,bi,bj)*CVV(I,JM))
                0581             ENDDO
                0582             DO J=1,sNy-1
                0583              J1=sNy-J
                0584              J2=J1+1
                0585              VRT(I,J1)=VRT(I,J1)-CVV(I,J1)*VRT(I,J2)
                0586             ENDDO
                0587             DO J=1,sNy
                0588              vIce(I,J,bi,bj)=vTmp(I,J,bi,bj)
                0589      &           +WFAV*(VRT(I,J)-vTmp(I,J,bi,bj))
                0590             ENDDO
                0591            ENDDO
                0592 C--   end doIterate4v
                0593           ENDIF
                0594 
                0595 C     end bi,bj-loops
                0596          ENDDO
                0597         ENDDO
                0598 
                0599         IF ( doIterate4u.AND.MOD(m,SOLV_NCHECK).EQ.0) THEN
                0600          S1=ZERO
                0601          DO bj=myByLo(myThid),myByHi(myThid)
                0602           DO bi=myBxLo(myThid),myBxHi(myThid)
                0603            DO J=1,sNy
                0604             DO I=1,sNx
                0605              UERR=(uIce(I,J,bi,bj)-uTmp(I,J,bi,bj))
                0606      &               * UVM(I,J,bi,bj)
                0607              S1=MAX(ABS(UERR),S1)
                0608             ENDDO
                0609            ENDDO
                0610           ENDDO
                0611          ENDDO
                0612          _GLOBAL_MAX_RL( S1, myThid )
                0613 c        WRITE(standardMessageUnit,'(A,2I6,1P4E16.9)')
                0614 c    &   ' U iters,error,WF = ',ilcall,M,S1,S1A,WFAU
                0615 C SAFEGUARD AGAINST BAD FORCING ETC
                0616          IF(m.GT.1.AND.S1.GT.S1A) WFAU=WFAU2
                0617          S1A=S1
                0618          IF(S1.LT.LSR_ERROR) THEN
                0619           ICOUNT1=m
                0620           doIterate4u = .FALSE.
                0621          ENDIF
                0622         ENDIF
                0623 
                0624         IF ( doIterate4v.AND.MOD(m,SOLV_NCHECK).EQ.0) THEN
                0625          S2=ZERO
                0626          DO bj=myByLo(myThid),myByHi(myThid)
                0627           DO bi=myBxLo(myThid),myBxHi(myThid)
                0628            DO J=1,sNy
                0629             DO I=1,sNx
                0630              UERR=(vIce(I,J,bi,bj)-vTmp(I,J,bi,bj))
                0631      &               * UVM(I,J,bi,bj)
                0632              S2=MAX(ABS(UERR),S2)
                0633             ENDDO
                0634            ENDDO
                0635           ENDDO
                0636          ENDDO
                0637          _GLOBAL_MAX_RL( S2, myThid )
                0638 C SAFEGUARD AGAINST BAD FORCING ETC
                0639          IF(m.GT.1.AND.S2.GT.S2A) WFAV=WFAV2
                0640          S2A=S2
                0641          IF(S2.LT.LSR_ERROR) THEN
                0642           ICOUNT2=m
                0643           doIterate4v = .FALSE.
                0644          ENDIF
                0645         ENDIF
                0646 
                0647         CALL EXCH_UV_XY_RL( uIce, vIce,.TRUE.,myThid)
                0648 
                0649 C--    end doIterate4u or doIterate4v
                0650        ENDIF
                0651       ENDDO
                0652 C ITERATION END -----------------------------------------------------
                0653 #ifdef ALLOW_DEBUG
                0654       IF ( debugLevel .GE. debLevD ) THEN
                0655         WRITE(msgBuf,'(A,I3,A)')
                0656      &        'Uice post iter (SEAICE_LSR', MOD(ilcall,1000), ')'
                0657         CALL DEBUG_STATS_RL( 1, uIce, msgBuf, myThid )
                0658         WRITE(msgBuf,'(A,I3,A)')
                0659      &        'Vice post iter (SEAICE_LSR', MOD(ilcall,1000), ')'
                0660         CALL DEBUG_STATS_RL( 1, vIce, msgBuf, myThid )
                0661       ENDIF
                0662 #endif /* ALLOW_DEBUG */
                0663       IF ( doIterate4u .OR. doIterate4v ) THEN
                0664         WRITE(msgBuf,'(2A,I10,A,I4,A)') '** WARNING ** SEAICE_LSR ',
                0665      &    '(it=', -9999, ',', ilcall,') did not converge :'
                0666         CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
                0667      &                      SQUEEZE_RIGHT, myThid )
                0668         WRITE(msgBuf,'(2(A,I6,0PF6.3,1PE14.6))')
                0669      &      ' nIt,wFU,dU=', ICOUNT1, WFAU, S1,
                0670      &    ' ; nIt,wFV,dV=', ICOUNT2, WFAV, S2
                0671         CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
                0672      &                      SQUEEZE_RIGHT, myThid )
                0673       ENDIF
                0674 
                0675 #else /* SEAICE_LSRBNEW */
809c36b928 Patr*0676 C SOLVE FOR UICE
                0677 
7109a141b2 Patr*0678 #ifdef ALLOW_AUTODIFF_TAMC
88f72205aa Jean*0679 cph That is an important one! Note, that
7109a141b2 Patr*0680 cph * lsr is called twice, thus the icall index
                0681 cph * this storing is still outside the iteration loop
8a37a22def Patr*0682 CADJ STORE uice = comlev1_dynsol,
7109a141b2 Patr*0683 CADJ &            key = ikey_dynamics + (ilcall-1)*nchklev_1
8a37a22def Patr*0684 CADJ STORE vice = comlev1_dynsol,
7109a141b2 Patr*0685 CADJ &            key = ikey_dynamics + (ilcall-1)*nchklev_1
                0686 #endif /* ALLOW_AUTODIFF_TAMC */
                0687 
809c36b928 Patr*0688       DO bj=myByLo(myThid),myByHi(myThid)
                0689        DO bi=myBxLo(myThid),myBxHi(myThid)
d24acf3cfc Jean*0690         DO j=1-OLy,sNy+OLy
                0691          DO i=1-OLx,sNx+OLx
4366d31d92 Mart*0692           etaPlusZeta(I,J,bi,bj) = ETA(I,J,bi,bj)+ZETA(I,J,bi,bj)
                0693           zetaMinusEta(I,J,bi,bj) = ZETA(I,J,bi,bj)-ETA(I,J,bi,bj)
                0694          ENDDO
                0695         ENDDO
d24acf3cfc Jean*0696         DO j=1-OLy+1,sNy+OLy
                0697          DO i=1-OLx+1,sNx+OLx
4366d31d92 Mart*0698           ETAMEAN(I,J,bi,bj) =QUART*(
                0699      &          ETA(I,J-1,bi,bj) + ETA(I-1,J-1,bi,bj)
                0700      &         +ETA(I,J  ,bi,bj) + ETA(I-1,J  ,bi,bj))
                0701           ZETAMEAN(I,J,bi,bj)=QUART*(
                0702      &          ZETA(I,J-1,bi,bj) + ZETA(I-1,J-1,bi,bj)
                0703      &         +ZETA(I,J  ,bi,bj) + ZETA(I-1,J  ,bi,bj))
809c36b928 Patr*0704          ENDDO
                0705         ENDDO
                0706        ENDDO
                0707       ENDDO
                0708 
                0709       DO bj=myByLo(myThid),myByHi(myThid)
                0710        DO bi=myBxLo(myThid),myBxHi(myThid)
                0711 
                0712         DO J=1,sNy
                0713          DO I=1,sNx
4366d31d92 Mart*0714           AA1=( etaPlusZeta(I  ,J-1,bi,bj) * _recip_dxF(I  ,J-1,bi,bj)
                0715      &         +etaPlusZeta(I  ,J  ,bi,bj) * _recip_dxF(I  ,J  ,bi,bj)
                0716      &         )*0.5 _d 0 * _recip_dxV(I,J,bi,bj) * UVM(I,J,bi,bj)
                0717           AA2=( etaPlusZeta(I-1,J-1,bi,bj) * _recip_dxF(I-1,J-1,bi,bj)
                0718      &         +etaPlusZeta(I-1,J  ,bi,bj) * _recip_dxF(I-1,J  ,bi,bj)
                0719      &         )*0.5 _d 0 * _recip_dxV(I,J,bi,bj) * UVM(I,J,bi,bj)
                0720           AA3=  0.5 _d 0 *(ETA(I-1,J  ,bi,bj)+ETA(I,J  ,bi,bj))
                0721           AA4=  0.5 _d 0 *(ETA(I-1,J-1,bi,bj)+ETA(I,J-1,bi,bj))
df7ba7118e Jean*0722           AA5= -(AA3-AA4) * _tanPhiAtV(I,J,bi,bj)
1b4a08ade9 Mart*0723      &         * _recip_dyU(I,J,bi,bj)*recip_rSphere
                0724           AA6=TWO*ETAMEAN(I,J,bi,bj) *recip_rSphere*recip_rSphere
4366d31d92 Mart*0725      &          * _tanPhiAtV(I,J,bi,bj)  * _tanPhiAtV(I,J,bi,bj)
                0726           AU(I,J,bi,bj)=-AA2
                0727           CU(I,J,bi,bj)=-AA1
                0728           BU(I,J,bi,bj)=(ONE-UVM(I,J,bi,bj))
                0729      &         - AU(I,J,bi,bj) - CU(I,J,bi,bj)
                0730      &         + ((AA3+AA4)*_recip_dyU(I,J,bi,bj)*_recip_dyU(I,J,bi,bj)
                0731      &           + AA5 + AA6
                0732      &           + AMASS(I,J,bi,bj)/SEAICE_deltaTdyn
                0733      &           + DRAGS(I,J,bi,bj)
                0734      &           )*UVM(I,J,bi,bj)
809c36b928 Patr*0735          END DO
                0736         END DO
                0737 
                0738         DO J=1,sNy
09510da3bb Dimi*0739          AU(1,J,bi,bj)=ZERO
                0740          CU(sNx,J,bi,bj)=ZERO
                0741          CU(1,J,bi,bj)=CU(1,J,bi,bj)/BU(1,J,bi,bj)
809c36b928 Patr*0742         END DO
                0743 
4366d31d92 Mart*0744 C     now set up right-hand side
d24acf3cfc Jean*0745         DO J=1-OLy,sNy+OLy-1
                0746          DO I=1-OLx,sNx+OLx-1
bb4175c1b5 Mart*0747           dVdy(I,J) = 0.5 _d 0 * (
45315406aa Mart*0748      &         ( vIceB(I+1,J+1,bi,bj) - vIceB(I+1,J  ,bi,bj) )
bb4175c1b5 Mart*0749      &         * _recip_dyG(I+1,J,bi,bj)
45315406aa Mart*0750      &         +(vIceB(I  ,J+1,bi,bj) - vIceB(I  ,J  ,bi,bj) )
bb4175c1b5 Mart*0751      &         * _recip_dyG(I,  J,bi,bj) )
                0752           dVdx(I,J) = 0.5 _d 0 * (
45315406aa Mart*0753      &         ( vIceB(I+1,J+1,bi,bj) - vIceB(I  ,J+1,bi,bj) )
bb4175c1b5 Mart*0754      &         * _recip_dxG(I,J+1,bi,bj)
45315406aa Mart*0755      &         +(vIceB(I+1,J  ,bi,bj) - vIceB(I  ,J  ,bi,bj) )
bb4175c1b5 Mart*0756      &         * _recip_dxG(I,J,  bi,bj) )
                0757          ENDDO
                0758         ENDDO
                0759 #ifdef SEAICE_TEST
d24acf3cfc Jean*0760         DO j=1-OLy,sNy+OLy-1
                0761          DO i=1-OLx,sNx+OLx-1
04016f2c47 Mart*0762           vz(i,j) = quart * (
45315406aa Mart*0763      &         vIceB(i,j,bi,bj) + vIceB(i+1,j,bi,bj) )
04016f2c47 Mart*0764           vz(i,j)= vz(i,j) + quart * (
45315406aa Mart*0765      &         vIceB(i,j+1,bi,bj) + vIceB(i+1,j+1,bi,bj) )
4366d31d92 Mart*0766          ENDDO
                0767         ENDDO
bb4175c1b5 Mart*0768 #endif
809c36b928 Patr*0769         DO J=1,sNy
                0770          DO I=1,sNx
45315406aa Mart*0771           rhsU(I,J,bi,bj)=DRAGA(I,J,bi,bj)*vIceB(I,J,bi,bj)
4366d31d92 Mart*0772      &         +FORCEX(I,J,bi,bj)
bb4175c1b5 Mart*0773 #ifdef SEAICE_TEST
df7ba7118e Jean*0774      &        + ( 0.5 _d 0 *
bb4175c1b5 Mart*0775      &         (zetaMinusEta(i,j,bi,bj)+zetaMinusEta(i,j-1,bi,bj))
04016f2c47 Mart*0776      &         *(vz(i,j)-vz(i,j-1)) * _recip_dyC(i,j,bi,bj)
df7ba7118e Jean*0777      &          - 0.5 _d 0 *
bb4175c1b5 Mart*0778      &         (zetaMinusEta(i-1,j,bi,bj)+zetaMinusEta(i-1,j-1,bi,bj))
04016f2c47 Mart*0779      &         *(vz(i-1,j)-vz(i-1,j-1)) * _recip_dyC(i-1,j,bi,bj)
bb4175c1b5 Mart*0780      &         ) * _recip_dxV(i,j,bi,bj)
                0781 #else
                0782      &         + ( zetaMinusEta(I  ,J  ,bi,bj) * dVdy(I  ,J  )
                0783      &           + zetaMinusEta(I  ,J-1,bi,bj) * dVdy(I  ,J-1)
                0784      &           - zetaMinusEta(I-1,J  ,bi,bj) * dVdy(I-1,J  )
                0785      &           - zetaMinusEta(I-1,J-1,bi,bj) * dVdy(I-1,J-1)
df7ba7118e Jean*0786      &         )* 0.5 _d 0 * _recip_dxV(I,J,bi,bj)
bb4175c1b5 Mart*0787 #endif
4366d31d92 Mart*0788      &
bb4175c1b5 Mart*0789      &         + ( ETA         (I  ,J  ,bi,bj) * dVdx(I  ,J  )
                0790      &           + ETA         (I-1,J  ,bi,bj) * dVdx(I-1,J  )
                0791      &           - ETA         (I  ,J-1,bi,bj) * dVdx(I  ,J-1)
                0792      &           - ETA         (I-1,J-1,bi,bj) * dVdx(I-1,J-1)
4366d31d92 Mart*0793      &         ) * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)
                0794      &
                0795      &        -(etaPlusZeta(I  ,J  ,bi,bj)+etaPlusZeta(I  ,J-1,bi,bj)
                0796      &         -etaPlusZeta(I-1,J-1,bi,bj)-etaPlusZeta(I-1,J  ,bi,bj))
45315406aa Mart*0797      &         * vIceB(I,J,bi,bj)
4366d31d92 Mart*0798      &            * _tanPhiAtV(I,J,bi,bj)
1b4a08ade9 Mart*0799      &         * 0.5 _d 0 * _recip_dxV(I,J,bi,bj)*recip_rSphere
4366d31d92 Mart*0800      &
                0801      &         -(ETAMEAN(I,J,bi,bj)+ZETAMEAN(I,J,bi,bj))
45315406aa Mart*0802      &         *(vIceB(I+1,J,bi,bj) - vIceB(I-1,J,bi,bj))
4366d31d92 Mart*0803      &            * _tanPhiAtV(I,J,bi,bj)
                0804      &         * 1.0 _d 0 /( _dxG(I,J,bi,bj) + _dxG(I-1,J,bi,bj) )
1b4a08ade9 Mart*0805      &         *recip_rSphere
4366d31d92 Mart*0806      &
                0807      &         -ETAMEAN(I,J,bi,bj)
45315406aa Mart*0808      &         *(vIceB(I+1,J,bi,bj) - vIceB(I-1,J,bi,bj))
4366d31d92 Mart*0809      &            *TWO* _tanPhiAtV(I,J,bi,bj)
                0810      &         * 1.0 _d 0 /( _dxG(I,J,bi,bj) + _dxG(I-1,J,bi,bj) )
1b4a08ade9 Mart*0811      &         *recip_rSphere
4366d31d92 Mart*0812 
e1d2983c5b Mart*0813           URT1(I,J,bi,bj)=
4366d31d92 Mart*0814      &         0.5 _d 0 * (ETA(I-1,J-1,bi,bj)+ETA(I,J-1,bi,bj))
                0815      &         * _recip_dyU(I,J,bi,bj) * _recip_dyU(I,J,bi,bj)
                0816      &         -     ETAMEAN(I,J,bi,bj) * _tanPhiAtV(I,J-1,bi,bj)
1b4a08ade9 Mart*0817      &          * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere
4366d31d92 Mart*0818      &         + TWO*ETAMEAN(I,J,bi,bj) * _tanPhiAtV(I,J,bi,bj)
1b4a08ade9 Mart*0819      &          * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere
e1d2983c5b Mart*0820           URT2(I,J,bi,bj)=
4366d31d92 Mart*0821      &         0.5 _d 0 * (ETA(I-1,J,bi,bj)+ETA(I,J,bi,bj))
                0822      &         * _recip_dyU(I,J,bi,bj) * _recip_dyU(I,J,bi,bj)
                0823      &         +     ETAMEAN(I,J,bi,bj) * _tanPhiAtV(I,J+1,bi,bj)
1b4a08ade9 Mart*0824      &          * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere
4366d31d92 Mart*0825      &         - TWO*ETAMEAN(I,J,bi,bj) * _tanPhiAtV(I,J,bi,bj)
1b4a08ade9 Mart*0826      &          * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere
809c36b928 Patr*0827          END DO
                0828         END DO
                0829 
                0830        ENDDO
                0831       ENDDO
                0832 
                0833 C NOW DO ITERATION
09510da3bb Dimi*0834 
7109a141b2 Patr*0835 cph--- iteration starts here
                0836 cph--- need to kick out goto
                0837       phexit = -1
                0838 
                0839 C ITERATION START -----------------------------------------------------
                0840 #ifdef ALLOW_AUTODIFF_TAMC
                0841 CADJ LOOP = iteration uice
                0842 #endif /* ALLOW_AUTODIFF_TAMC */
79df32c3f1 Mart*0843       DO M=1, SEAICElinearIterMax
7109a141b2 Patr*0844       IF ( phexit .EQ. -1 ) THEN
df7ba7118e Jean*0845 
809c36b928 Patr*0846       DO bj=myByLo(myThid),myByHi(myThid)
                0847        DO bi=myBxLo(myThid),myBxHi(myThid)
                0848 C NOW SET U(3)=U(1)
                0849         DO J=1,sNy
                0850          DO I=1,sNx
772590b63c Mart*0851           UTMP(I,J,bi,bj)=UICE(I,J,bi,bj)
809c36b928 Patr*0852          END DO
                0853         END DO
                0854 
df7ba7118e Jean*0855         DO J=1,sNy
809c36b928 Patr*0856          DO I=1,sNx
09510da3bb Dimi*0857           IF(I.EQ.1) THEN
4366d31d92 Mart*0858            AA2=(etaPlusZeta(I-1,J-1,bi,bj) * _recip_dxF(I-1,J-1,bi,bj)
                0859      &         +etaPlusZeta(I-1,J  ,bi,bj) * _recip_dxF(I-1,J  ,bi,bj)
                0860      &          )*0.5 _d 0 * _recip_dxV(I,J,bi,bj)
772590b63c Mart*0861            AA3=AA2*UICE(I-1,J,bi,bj)*UVM(I,J,bi,bj)
09510da3bb Dimi*0862           ELSE IF(I.EQ.sNx) THEN
4366d31d92 Mart*0863            AA1=(etaPlusZeta(I  ,J-1,bi,bj) * _recip_dxF(I  ,J-1,bi,bj)
                0864      &         +etaPlusZeta(I  ,J  ,bi,bj) * _recip_dxF(I  ,J  ,bi,bj)
                0865      &          )*0.5 _d 0 * _recip_dxV(I,J,bi,bj)
772590b63c Mart*0866            AA3=AA1*UICE(I+1,J,bi,bj)*UVM(I,J,bi,bj)
09510da3bb Dimi*0867           ELSE
                0868            AA3=ZERO
                0869           END IF
e1d2983c5b Mart*0870           URT(I)=rhsU(I,J,bi,bj)+AA3
                0871      &          +URT1(I,J,bi,bj)*UICE(I,J-1,bi,bj)
                0872      &          +URT2(I,J,bi,bj)*UICE(I,J+1,bi,bj)
809c36b928 Patr*0873           URT(I)=URT(I)*UVM(I,J,bi,bj)
                0874          END DO
                0875 
                0876          DO I=1,sNx
09510da3bb Dimi*0877           CUU(I)=CU(I,J,bi,bj)
809c36b928 Patr*0878          END DO
09510da3bb Dimi*0879          URT(1)=URT(1)/BU(1,J,bi,bj)
809c36b928 Patr*0880          DO I=2,sNx
                0881           IM=I-1
09510da3bb Dimi*0882           CUU(I)=CUU(I)/(BU(I,J,bi,bj)-AU(I,J,bi,bj)*CUU(IM))
                0883           URT(I)=(URT(I)-AU(I,J,bi,bj)*URT(IM))
                0884      &          /(BU(I,J,bi,bj)-AU(I,J,bi,bj)*CUU(IM))
809c36b928 Patr*0885          END DO
                0886          DO I=1,sNx-1
                0887           J1=sNx-I
                0888           J2=J1+1
                0889           URT(J1)=URT(J1)-CUU(J1)*URT(J2)
                0890          END DO
                0891          DO I=1,sNx
772590b63c Mart*0892           UICE(I,J,bi,bj)=UTMP(I,J,bi,bj)
04016f2c47 Mart*0893      &        +WFAU*(URT(I)-UTMP(I,J,bi,bj))
809c36b928 Patr*0894          END DO
                0895 
df7ba7118e Jean*0896         END DO
809c36b928 Patr*0897 
                0898        ENDDO
                0899       ENDDO
                0900 
09510da3bb Dimi*0901       IF(MOD(M,SOLV_NCHECK).EQ.0) THEN
df7ba7118e Jean*0902        S1=ZERO
09510da3bb Dimi*0903        DO bj=myByLo(myThid),myByHi(myThid)
                0904         DO bi=myBxLo(myThid),myBxHi(myThid)
                0905          DO J=1,sNy
                0906           DO I=1,sNx
e1d2983c5b Mart*0907            UERR=(UICE(I,J,bi,bj)-UTMP(I,J,bi,bj))
809c36b928 Patr*0908      &             *UVM(I,J,bi,bj)
e1d2983c5b Mart*0909            S1=MAX(ABS(UERR),S1)
09510da3bb Dimi*0910           END DO
809c36b928 Patr*0911          END DO
09510da3bb Dimi*0912         ENDDO
809c36b928 Patr*0913        ENDDO
7163a40534 Jean*0914        _GLOBAL_MAX_RL( S1, myThid )
09510da3bb Dimi*0915 C SAFEGUARD AGAINST BAD FORCING ETC
                0916        IF(M.GT.1.AND.S1.GT.S1A) WFAU=WFAU2
                0917        S1A=S1
                0918        IF(S1.LT.LSR_ERROR) THEN
                0919         ICOUNT1=M
7109a141b2 Patr*0920         phexit = 1
09510da3bb Dimi*0921        END IF
                0922       END IF
772590b63c Mart*0923       _EXCH_XY_RL( UICE, myThid )
09510da3bb Dimi*0924 
df7ba7118e Jean*0925       ENDIF
                0926       ENDDO
7109a141b2 Patr*0927 C ITERATION END -----------------------------------------------------
                0928 
23142459d0 Jean*0929       IF ( debugLevel .GE. debLevC ) THEN
0276a7ecde Dimi*0930        _BEGIN_MASTER( myThid )
                0931         write(*,'(A,I6,1P2E22.14)')' U lsr iters, error = ',ICOUNT1,S1
ea6e02f692 Ed H*0932        _END_MASTER( myThid )
0276a7ecde Dimi*0933       ENDIF
809c36b928 Patr*0934 
                0935 C NOW FOR VICE
                0936       DO bj=myByLo(myThid),myByHi(myThid)
                0937        DO bi=myBxLo(myThid),myBxHi(myThid)
                0938 
                0939         DO J=1,sNy
                0940          DO I=1,sNx
4366d31d92 Mart*0941           AA1=0.5 _d 0 * _recip_dyU(I,J,bi,bj) * _recip_dyU(I,J,bi,bj)
                0942      &         * (etaPlusZeta(I-1,J  ,bi,bj) + etaPlusZeta(I,J  ,bi,bj))
                0943           AA2=0.5 _d 0 * _recip_dyU(I,J,bi,bj) * _recip_dyU(I,J,bi,bj)
                0944      &         * (etaPlusZeta(I-1,J-1,bi,bj) + etaPlusZeta(I,J-1,bi,bj))
                0945           AA3= (ETA(I  ,J-1,bi,bj) * _recip_dxV(I,J,bi,bj)
                0946      &         +ETA(I  ,J  ,bi,bj) * _recip_dxV(I,J,bi,bj)
                0947      &         )* 0.5 _d 0 * _recip_dxV(I,J,bi,bj)
                0948           AA4= (ETA(I-1,J-1,bi,bj)+ETA(I-1,J,bi,bj))*0.5 _d 0
                0949      &          *_recip_dxV(I,J,bi,bj) * _recip_dxV(I,J,bi,bj)
                0950           AA5=(zetaMinusEta(I-1,J  ,bi,bj) + zetaMinusEta(I,J  ,bi,bj)
                0951      &        -zetaMinusEta(I-1,J-1,bi,bj) - zetaMinusEta(I,J-1,bi,bj)
df7ba7118e Jean*0952      &          )* _tanPhiAtV(I,J,bi,bj)
1b4a08ade9 Mart*0953      &         * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere
4366d31d92 Mart*0954 
1b4a08ade9 Mart*0955           AA6=TWO*ETAMEAN(I,J,bi,bj) * recip_rSphere*recip_rSphere
4366d31d92 Mart*0956      &         * _tanPhiAtV(I,J,bi,bj) * _tanPhiAtV(I,J,bi,bj)
                0957 
                0958          AV(I,J,bi,bj)=(
                0959      &         - AA2
                0960      &         - (ZETAMEAN(I,J,bi,bj)-ETAMEAN(I,J,bi,bj))
df7ba7118e Jean*0961      &          * _tanPhiAtV(I,J-1,bi,bj)
1b4a08ade9 Mart*0962      &          * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere
4366d31d92 Mart*0963      &         -ETAMEAN(I,J,bi,bj)*TWO* _tanPhiAtV(I,J,bi,bj)
1b4a08ade9 Mart*0964      &          * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere
4366d31d92 Mart*0965      &         )*UVM(I,J,bi,bj)
                0966          CV(I,J,bi,bj)=(
df7ba7118e Jean*0967      &         -AA1
4366d31d92 Mart*0968      &        +(ZETAMEAN(I,J,bi,bj)-ETAMEAN(I,J,bi,bj))
df7ba7118e Jean*0969      &         * _tanPhiAtV(I,J+1,bi,bj)
1b4a08ade9 Mart*0970      &         * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere
4366d31d92 Mart*0971      &        +ETAMEAN(I,J,bi,bj)*TWO* _tanPhiAtV(I,J,bi,bj)
1b4a08ade9 Mart*0972      &         * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere
4366d31d92 Mart*0973      &        )*UVM(I,J,bi,bj)
                0974           BV(I,J,bi,bj)= (ONE-UVM(I,J,bi,bj))
                0975      &        +( (AA1+AA2) + (AA3+AA4)  + AA5 + AA6
                0976      &        +AMASS(I,J,bi,bj)/SEAICE_deltaTdyn+DRAGS(I,J,bi,bj))
                0977      &        *UVM(I,J,bi,bj)
809c36b928 Patr*0978          END DO
                0979         END DO
                0980 
                0981         DO I=1,sNx
09510da3bb Dimi*0982          AV(I,1,bi,bj)=ZERO
                0983          CV(I,sNy,bi,bj)=ZERO
                0984          CV(I,1,bi,bj)=CV(I,1,bi,bj)/BV(I,1,bi,bj)
809c36b928 Patr*0985         END DO
                0986 
4366d31d92 Mart*0987 C     now set up right-hand-side
d24acf3cfc Jean*0988         DO J=1-OLy,sNy+OLy-1
                0989          DO I=1-OLx,sNx+OLx-1
bb4175c1b5 Mart*0990           dUdx(I,J) = 0.5 _d 0 * (
45315406aa Mart*0991      &         ( uIceB(I+1,J+1,bi,bj) - uIceB(I ,J+1,bi,bj) )
bb4175c1b5 Mart*0992      &         * _recip_dxG(I,J+1,bi,bj)
45315406aa Mart*0993      &         +(uIceB(I+1,J  ,bi,bj) - uIceB(I ,J  ,bi,bj) )
bb4175c1b5 Mart*0994      &         * _recip_dxG(I,J  ,bi,bj) )
                0995           dUdy(I,J) = 0.5 _d 0 * (
45315406aa Mart*0996      &          ( uIceB(I+1,J+1,bi,bj) - uIceB(I+1,J  ,bi,bj) )
bb4175c1b5 Mart*0997      &         * _recip_dyG(I+1,J,bi,bj)
45315406aa Mart*0998      &          +(uIceB(I  ,J+1,bi,bj) - uIceB(I  ,J  ,bi,bj) )
bb4175c1b5 Mart*0999      &         * _recip_dyG(I,  J,bi,bj) )
4366d31d92 Mart*1000          ENDDO
                1001         ENDDO
bb4175c1b5 Mart*1002 #ifdef SEAICE_TEST
d24acf3cfc Jean*1003         DO j=1-OLy,sNy+OLy-1
                1004          DO i=1-OLx,sNx+OLx-1
04016f2c47 Mart*1005           uz(i,j) = quart * (
45315406aa Mart*1006      &         uIceB(i,j,bi,bj) + uIceB(i+1,j,bi,bj) )
04016f2c47 Mart*1007           uz(i,j)= uz(i,j) + quart * (
45315406aa Mart*1008      &         uIceB(i,j+1,bi,bj) + uIceB(i+1,j+1,bi,bj) )
2d0930fb3f Mart*1009          ENDDO
                1010         ENDDO
bb4175c1b5 Mart*1011 #endif
809c36b928 Patr*1012         DO J=1,sNy
                1013          DO I=1,sNx
45315406aa Mart*1014            rhsV(I,J,bi,bj)=-DRAGA(I,J,bi,bj)*uIceB(I,J,bi,bj)
4366d31d92 Mart*1015      &        +FORCEY(I,J,bi,bj)
                1016      &
bb4175c1b5 Mart*1017 #ifdef SEAICE_TEST
df7ba7118e Jean*1018      &        + ( 0.5 _d 0 *
2d0930fb3f Mart*1019      &         (zetaMinusEta(i,j,bi,bj)+zetaMinusEta(i-1,j,bi,bj))
04016f2c47 Mart*1020      &         *(uz(i,j)-uz(i-1,j)) * _recip_dxC(i,j,bi,bj)
df7ba7118e Jean*1021      &          - 0.5 _d 0 *
2d0930fb3f Mart*1022      &         (zetaMinusEta(i,j-1,bi,bj)+zetaMinusEta(i-1,j-1,bi,bj))
04016f2c47 Mart*1023      &         *(uz(i,j-1)-uz(i-1,j-1)) * _recip_dxC(i,j-1,bi,bj)
2d0930fb3f Mart*1024      &         ) * _recip_dyU(i,j,bi,bj)
bb4175c1b5 Mart*1025 #else
                1026      &         + ( zetaMinusEta(I  ,J  ,bi,bj) * dUdx(I  ,J  )
                1027      &           + zetaMinusEta(I-1,J  ,bi,bj) * dUdx(I-1,J  )
                1028      &           - zetaMinusEta(I  ,J-1,bi,bj) * dUdx(I  ,J-1)
                1029      &           - zetaMinusEta(I-1,J-1,bi,bj) * dUdx(I-1,J-1)
                1030      &         )* 0.5 _d 0 * _recip_dyU(I,J,bi,bj)
                1031 #endif
df7ba7118e Jean*1032      &
bb4175c1b5 Mart*1033      &        + ( ETA          (I  ,J  ,bi,bj) * dUdy(I  ,J  )
                1034      &         +  ETA          (I  ,J-1,bi,bj) * dUdy(I  ,J-1)
                1035      &         -  ETA          (I-1,J  ,bi,bj) * dUdy(I-1,J  )
                1036      &         -  ETA          (I-1,J-1,bi,bj) * dUdy(I-1,J-1)
                1037      &         )*0.5 _d 0* _recip_dxV(I,J,bi,bj)
4366d31d92 Mart*1038      &
                1039      &        +(ETA(I  ,J  ,bi,bj) + ETA(I  ,J-1,bi,bj)
                1040      &         -ETA(I-1,J-1,bi,bj) - ETA(I-1,J  ,bi,bj))
45315406aa Mart*1041      &         * uIceB(I,J,bi,bj)
4366d31d92 Mart*1042      &         * _tanPhiAtV(I,J,bi,bj)
1b4a08ade9 Mart*1043      &         * 0.5 _d 0 * _recip_dxV(I,J,bi,bj)*recip_rSphere
4366d31d92 Mart*1044      &        +ETAMEAN(I,J,bi,bj) * _tanPhiAtV(I,J,bi,bj)
45315406aa Mart*1045      &         *(uIceB(I+1,J,bi,bj)-uIceB(I-1,J,bi,bj))
1b4a08ade9 Mart*1046      &         * 0.5 _d 0 * _recip_dxV(I,J,bi,bj)*recip_rSphere
4366d31d92 Mart*1047      &
                1048      &        +ETAMEAN(I,J,bi,bj)*TWO  * _tanPhiAtV(I,J,bi,bj)
45315406aa Mart*1049      &        *(uIceB(I+1,J,bi,bj)-uIceB(I-1,J,bi,bj))
4366d31d92 Mart*1050      &         * 1. _d 0 /( _dxG(I,J,bi,bj) + _dxG(I-1,J,bi,bj))
1b4a08ade9 Mart*1051      &         *recip_rSphere
e1d2983c5b Mart*1052           VRT1(I,J,bi,bj)= 0.5 _d 0 * (
4366d31d92 Mart*1053      &           ETA(I-1,J-1,bi,bj) * _recip_dxV(I,J,bi,bj)
                1054      &          +ETA(I-1,J  ,bi,bj) * _recip_dxV(I,J,bi,bj)
                1055      &          ) * _recip_dxV(I,J,bi,bj)
e1d2983c5b Mart*1056           VRT2(I,J,bi,bj)= 0.5 _d 0 * (
4366d31d92 Mart*1057      &          ETA(I  ,J-1,bi,bj) * _recip_dxV(I,J,bi,bj)
                1058      &         +ETA(I  ,J  ,bi,bj) * _recip_dxV(I,J,bi,bj)
                1059      &         ) * _recip_dxV(I,J,bi,bj)
809c36b928 Patr*1060 
                1061          END DO
                1062         END DO
                1063 
                1064        ENDDO
                1065       ENDDO
                1066 
                1067 C NOW DO ITERATION
09510da3bb Dimi*1068 
7109a141b2 Patr*1069 cph--- iteration starts here
                1070 cph--- need to kick out goto
                1071       phexit = -1
                1072 
                1073 C ITERATION START -----------------------------------------------------
                1074 #ifdef ALLOW_AUTODIFF_TAMC
                1075 CADJ LOOP = iteration vice
                1076 #endif /* ALLOW_AUTODIFF_TAMC */
79df32c3f1 Mart*1077       DO M=1, SEAICElinearIterMax
7109a141b2 Patr*1078       IF ( phexit .EQ. -1 ) THEN
df7ba7118e Jean*1079 
809c36b928 Patr*1080 C NOW SET U(3)=U(1)
                1081       DO bj=myByLo(myThid),myByHi(myThid)
                1082        DO bi=myBxLo(myThid),myBxHi(myThid)
                1083 
                1084         DO J=1,sNy
                1085          DO I=1,sNx
772590b63c Mart*1086           VTMP(I,J,bi,bj)=VICE(I,J,bi,bj)
809c36b928 Patr*1087          END DO
                1088         END DO
                1089 
7109a141b2 Patr*1090         DO I=1,sNx
809c36b928 Patr*1091          DO J=1,sNy
09510da3bb Dimi*1092            IF(J.EQ.1) THEN
4366d31d92 Mart*1093             AA2= _recip_dyU(I,J,bi,bj) * _recip_dyU(I,J,bi,bj)
df7ba7118e Jean*1094      &           * 0.5 _d 0 *(
4366d31d92 Mart*1095      &           etaPlusZeta(I-1,J-1,bi,bj) + etaPlusZeta(I,J-1,bi,bj)
                1096      &           )
                1097               AA3=( AA2
                1098      &         +( ZETAMEAN(I,J,bi,bj)-ETAMEAN(I,J,bi,bj) )
                1099      &         * _tanPhiAtV(I,J-1,bi,bj)
df7ba7118e Jean*1100      &         * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere
4366d31d92 Mart*1101      &         + ETAMEAN(I,J,bi,bj)*TWO* _tanPhiAtV(I,J,bi,bj)
1b4a08ade9 Mart*1102      &         * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere )
772590b63c Mart*1103      &         *VICE(I,J-1,bi,bj)*UVM(I,J,bi,bj)
09510da3bb Dimi*1104            ELSE IF(J.EQ.sNy) THEN
4366d31d92 Mart*1105             AA1= _recip_dyU(I,J,bi,bj) * _recip_dyU(I,J,bi,bj)
                1106      &           * 0.5 _d 0 * (
                1107      &           etaPlusZeta(I-1,J,bi,bj) + etaPlusZeta(I,J,bi,bj)
                1108      &            )
                1109             AA3=( AA1
                1110      &         -( ZETAMEAN(I,J,bi,bj)-ETAMEAN(I,J,bi,bj))
                1111      &         * _tanPhiAtV(I,J+1,bi,bj)
df7ba7118e Jean*1112      &         * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere
4366d31d92 Mart*1113      &         - ETAMEAN(I,J,bi,bj)*TWO* _tanPhiAtV(I,J,bi,bj)
1b4a08ade9 Mart*1114      &         * 0.5 _d 0 * _recip_dyU(I,J,bi,bj)*recip_rSphere )
772590b63c Mart*1115      &         *VICE(I,J+1,bi,bj)*UVM(I,J,bi,bj)
09510da3bb Dimi*1116            ELSE
                1117             AA3=ZERO
                1118            END IF
                1119 
e1d2983c5b Mart*1120           VRT(J)=rhsV(I,J,bi,bj)+AA3+VRT1(I,J,bi,bj)*VICE(I-1,J,bi,bj)
                1121      &                       +VRT2(I,J,bi,bj)*VICE(I+1,J,bi,bj)
809c36b928 Patr*1122           VRT(J)=VRT(J)*UVM(I,J,bi,bj)
                1123          END DO
                1124 
                1125          DO J=1,sNy
09510da3bb Dimi*1126           CVV(J)=CV(I,J,bi,bj)
809c36b928 Patr*1127          END DO
09510da3bb Dimi*1128          VRT(1)=VRT(1)/BV(I,1,bi,bj)
809c36b928 Patr*1129          DO J=2,sNy
                1130           JM=J-1
09510da3bb Dimi*1131           CVV(J)=CVV(J)/(BV(I,J,bi,bj)-AV(I,J,bi,bj)*CVV(JM))
                1132           VRT(J)=(VRT(J)-AV(I,J,bi,bj)*VRT(JM))
                1133      &          /(BV(I,J,bi,bj)-AV(I,J,bi,bj)*CVV(JM))
809c36b928 Patr*1134          END DO
                1135          DO J=1,sNy-1
                1136           J1=sNy-J
                1137           J2=J1+1
                1138           VRT(J1)=VRT(J1)-CVV(J1)*VRT(J2)
                1139          END DO
                1140          DO J=1,sNy
772590b63c Mart*1141           VICE(I,J,bi,bj)=VTMP(I,J,bi,bj)
04016f2c47 Mart*1142      &        +WFAV*(VRT(J)-VTMP(I,J,bi,bj))
809c36b928 Patr*1143          END DO
7109a141b2 Patr*1144         ENDDO
809c36b928 Patr*1145 
                1146        ENDDO
                1147       ENDDO
                1148 
09510da3bb Dimi*1149       IF(MOD(M,SOLV_NCHECK).EQ.0) THEN
df7ba7118e Jean*1150        S2=ZERO
09510da3bb Dimi*1151        DO bj=myByLo(myThid),myByHi(myThid)
                1152         DO bi=myBxLo(myThid),myBxHi(myThid)
                1153          DO J=1,sNy
                1154           DO I=1,sNx
e1d2983c5b Mart*1155            UERR=(VICE(I,J,bi,bj)-VTMP(I,J,bi,bj))
809c36b928 Patr*1156      &             *UVM(I,J,bi,bj)
e1d2983c5b Mart*1157            S2=MAX(ABS(UERR),S2)
09510da3bb Dimi*1158           END DO
809c36b928 Patr*1159          END DO
09510da3bb Dimi*1160         ENDDO
809c36b928 Patr*1161        ENDDO
7163a40534 Jean*1162        _GLOBAL_MAX_RL( S2, myThid )
09510da3bb Dimi*1163 C SAFEGUARD AGAINST BAD FORCING ETC
                1164        IF(M.GT.1.AND.S2.GT.S2A) WFAV=WFAV2
                1165        S2A=S2
                1166        IF(S2.LT.LSR_ERROR) THEN
                1167         ICOUNT2=M
7109a141b2 Patr*1168         phexit = 1
09510da3bb Dimi*1169        END IF
                1170       END IF
                1171 
772590b63c Mart*1172       _EXCH_XY_RL( VICE, myThid )
09510da3bb Dimi*1173 
df7ba7118e Jean*1174       ENDIF
                1175       ENDDO
7109a141b2 Patr*1176 C ITERATION END -----------------------------------------------------
                1177 
23142459d0 Jean*1178       IF ( debugLevel .GE. debLevC ) THEN
0276a7ecde Dimi*1179        _BEGIN_MASTER( myThid )
                1180         write(*,'(A,I6,1P2E22.14)')' V lsr iters, error = ',ICOUNT2,S2
ea6e02f692 Ed H*1181        _END_MASTER( myThid )
0276a7ecde Dimi*1182       ENDIF
809c36b928 Patr*1183 
e1d2983c5b Mart*1184 #endif /* SEAICE_LSRBNEW */
                1185 
809c36b928 Patr*1186 C NOW END
                1187 C NOW MAKE COROLIS TERM IMPLICIT
                1188       DO bj=myByLo(myThid),myByHi(myThid)
                1189        DO bi=myBxLo(myThid),myBxHi(myThid)
d24acf3cfc Jean*1190         DO J=1-OLy,sNy+OLy
                1191          DO I=1-OLx,sNx+OLx
772590b63c Mart*1192           UICE(I,J,bi,bj)=UICE(I,J,bi,bj)*UVM(I,J,bi,bj)
                1193           VICE(I,J,bi,bj)=VICE(I,J,bi,bj)*UVM(I,J,bi,bj)
809c36b928 Patr*1194          END DO
                1195         END DO
                1196        ENDDO
                1197       ENDDO
                1198 
45315406aa Mart*1199 #endif /* SEAICE_BGRID_DYNAMICS */
809c36b928 Patr*1200 
                1201       RETURN
                1202       END
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!