C*************************************************************
C PLEASE KEEP THIS NOTE OF MODEL-DEVELOPMENT HISTORY
C Matrix solve uses Thomas algorithm, 10/1991, Jinlun Zhang
C Spherical coordinate system, 10/27/93, Jinlun Zhang
C Latest finite differencing scheme for treatment of NP,
C 9/9/1996,Jinlun Zhang
C Alternating direction implicit (ADI) method is used, 10/1998,
C Jinlun Zhang
C For details about ADI dynamics model, see Zhang and Rothrock,
C "Modeling
C Arctic Sea ice with an efficient plastic solution",
C submitted to JGR, 1999
C Adapted for GISS coupled model Jiping Liu/Gavin Schmidt 2000
C - Further modularised May 2001
C*************************************************************


      MODULE ICEDYN 11,4
!@sum  ICEDYN holds local variables for dynamic sea ice
!@auth Gavin Schmidt (based on code from Jinlun Zhang)
      USE CONSTANT, only : radian,radius
      USE MODEL_COM, only : im,jm
      USE DOMAIN_DECOMP, only : DYN_GRID
      USE SEAICE, only : osurf_tilt
      IMPLICIT NONE
      SAVE

C**** Definition for ice advection grid (EDIT FOR ADVSI GRID CHANGE)
      INTEGER, PARAMETER :: IMIC=IM, JMIC=JM

C**** local grid variables for ice rheology scheme
C**** Edit the definition of nx1,ny1 to change the grid for the
C**** rheology calculations without changing ADVSI grid.

!@var nx1 number of grid points in the longitudinal direction
!@+  (calculated points from 2 through nx1-1. End points are boundaries)
!@var ny1 number of grid points in the latitudinal direction
!@+  (calculated points from 2 through ny1-1. End points are boundaries)
      integer, parameter :: nx1=imic+2, ny1=jmic
      INTEGER, parameter :: NYPOLE=NY1-1,NXLCYC=NX1-1
      integer :: NPOL=1,LCYC=1
       TYPE(DYN_GRID) :: grid_MIC
       TYPE(DYN_GRID) :: grid_NXY

!@var FOCEAN land/ocean mask on ice dynamic grid
      REAL*8, DIMENSION(:,:), ALLOCATABLE :: FOCEAN

C**** input
!@var HEFFM ice mass mask (1/0)
!@var UVM ice velocity mask (1/0)
!@var COR coriolis term for ice dynamic equation
!@var GAIRX,GAIRY atmosphere-ice stress (B grid)
!@var GWATX,GWATY ocean velocities (B grid) (m/s)
!@var PGFUB,PGFVB pressure accelaration force
!@var AMASS ice mass (kg/m^2)
!@var UICEC,VICEC velocity arrays  (m/s)
!@var UICE,VICE velocity arrays  (m/s)
!@var HEFF ice thickness (mean over box) (m)
!@var AREA ice area (frac)
!@var UIB,VIB velocity arrays (m/s)  (????)
C**** internal variables
!@var PRESS ice internal pressure (Pa)
!@var FORCEX,FORCEY external force 
!@var DRAGS,DRAGA symmetric/anti-symmetric drag terms
!@var ZMAX,ZMIN max,min values of ZETA
!@var ETA,ZETA viscosities
C**** output
!@var DWATN non-linear water drag term
!@var DMU,DMV ice-ocean stress
      REAL*8, DIMENSION(:,:), ALLOCATABLE :: PRESS,HEFFM,UVM,DWATN,COR
     *     ,ZMAX,ZMIN,ETA,ZETA,DRAGS,DRAGA,GAIRX,GAIRY
     *     ,GWATX,GWATY,PGFUB,PGFVB,FORCEX,FORCEY,AMASS,UICEC,VICEC,UIB
     *     ,VIB,DMU,DMV,HEFF,AREA
      REAL*8, DIMENSION(:,:,:), ALLOCATABLE :: UICE,VICE
C**** Geometry 
!@var SINEN sin(phi)
!@var BYDXDY 
!@var DXT,DXU x-direction distances on tracer and velocity grid
!@var DYT,DYU y-direction distances on tracer and velocity grid
      REAL*8, DIMENSION(:,:), ALLOCATABLE :: SINEN,BYDXDY
      REAL*8, DIMENSION(NX1) :: DXT,DXU,BYDX2,BYDXR
      REAL*8, DIMENSION(:), ALLOCATABLE :: DYT,DYU,BYDY2,BYDYR,CST,
     &                                     CSU,TNGT,TNG,BYCSU

!@var OIPHI ice-ocean turning angle (25 degrees)
!@var ECCEN value of eccentricity for yield curve ellipse
      REAL*8, PARAMETER :: ECCEN=2.0, OIPHI=25d0*radian
!@var SINWAT,COSWAT sin and cos of ice-ocean turning angle
      REAL*8 SINWAT,COSWAT

!@var PSTAR maximum sea ice pressure (Pa)
      REAL*8, PARAMETER :: PSTAR=2.75d4

!@var BYDTS reciprocal of timestep in ice dynamics code
      REAL*8 :: BYDTS

      INTEGER :: CHECKSUM_UNIT

      CONTAINS 


      SUBROUTINE FORM 2,4
!@sum  FORM calculates ice dynamics input parameters for relaxation
!@auth Jiping Liu/Gavin Schmidt (based on code from J. Zhang)
!@ver  1.0
      USE DOMAIN_DECOMP, only : grid, GET, NORTH,SOUTH
      USE DOMAIN_DECOMP, ONLY : HALO_UPDATE, CHECKSUM

      IMPLICIT NONE
      INTEGER I,J
      REAL*8 AAA
      INTEGER :: J_0,J_1,J_0S,J_1S,J_0STG,J_1STG

C****
C**** Extract useful local domain parameters from "grid"
C****
      CALL GET(grid, J_STRT     =J_0,    J_STOP     =J_1,
     &               J_STRT_SKP =J_0S,   J_STOP_SKP =J_1S )

C****
C**** Set up non linear water drag
C****
      DO J=J_0,J_1
      DO I=1,NX1-1
        DWATN(I,J)=5.5*SQRT((UICE(I,J,1)-GWATX(I,J))**2
     1       +(VICE(I,J,1)-GWATY(I,J))**2)
      END DO
      END DO
C NOW SET UP SYMMETRIC DRAG
      DO J=J_0,J_1S
      DO I=1,NX1-1
        DRAGS(I,J)=DWATN(I,J)*COSWAT
      END DO
      END DO
C NOW SET UP ANTI SYMMETRIC DRAG PLUS CORIOLIS
      DO J=J_0,J_1S
      DO I=1,NX1-1
        IF(J.GT.NY1/2) THEN
          DRAGA(I,J)=DWATN(I,J)*SINWAT+COR(I,J)
        ELSE
          DRAGA(I,J)=DWATN(I,J)*(-SINWAT)+COR(I,J)
        END IF
      END DO
      END DO
C NOW SET UP FORCING FIELD
      DO J=J_0,J_1S
      DO I=1,NX1-1

C FIRST DO WIND
        FORCEX(I,J)=GAIRX(i,j)
        FORCEY(I,J)=GAIRY(i,j)

C NOW ADD IN CURRENT FORCE
        IF(J.GT.NY1/2) THEN
          FORCEX(I,J)=FORCEX(I,J)+DWATN(I,J)*(COSWAT*GWATX(I,J)
     1         -SINWAT*GWATY(I,J))
          FORCEY(I,J)=FORCEY(I,J)+DWATN(I,J)*(SINWAT*GWATX(I,J)
     1         +COSWAT*GWATY(I,J))
        ELSE
          FORCEX(I,J)=FORCEX(I,J)+DWATN(I,J)*(COSWAT*GWATX(I,J)
     1         +SINWAT*GWATY(I,J))
          FORCEY(I,J)=FORCEY(I,J)+DWATN(I,J)*(-SINWAT*GWATX(I,J)
     1         +COSWAT*GWATY(I,J))
        END IF
        
C     NOW ADD IN TILT 
        if (osurf_tilt.eq.1) then
C**** This assumes explicit knowledge of sea surface tilt
          FORCEX(I,J)=FORCEX(I,J)+AMASS(I,J)*PGFUB(I,J)
          FORCEY(I,J)=FORCEY(I,J)+AMASS(I,J)*PGFVB(I,J)
        else
C**** Otherwise estimate tilt using geostrophy
          FORCEX(I,J)=FORCEX(I,J)-COR(I,J)*GWATY(I,J)
          FORCEY(I,J)=FORCEY(I,J)+COR(I,J)*GWATX(I,J)
        end if
      END DO
      END DO

C NOW SET UP ICE PRESSURE AND VISCOSITIES
      DO J=J_0,J_1
      DO I=1,NX1
        PRESS(I,J)=PSTAR*HEFF(I,J)*EXP(-20.0*(1.0-AREA(I,J)))
        ZMAX(I,J)=(5d12/2d4)*PRESS(I,J)
c       ZMIN(I,J)=0.0D+00
        ZMIN(I,J)=4d8
      END DO
      END DO

      CALL PLAST

      if (grid%HAVE_NORTH_POLE) then
       AAA=0.0
       DO I=2,NX1-1
         AAA=AAA+PRESS(I,NY1-1)
       END DO
       AAA=AAA/FLOAT(NX1-2)
       DO I=1,NX1
         PRESS(I,NY1)=AAA
       END DO
      end if

 8481 CONTINUE

      DO J=J_0,J_1
        PRESS(1,J)=PRESS(NX1-1,J)
        PRESS(NX1,J)=PRESS(2,J)
      END DO

C NOW SET VISCOSITIES AND PRESSURE EQUAL TO ZERO AT OUTFLOW PTS

      DO J=J_0,J_1
      DO I=1,NX1
        PRESS(I,J)=PRESS(I,J)*HEFFM(I,J)
        ETA(I,J)=ETA(I,J)*HEFFM(I,J)
        ZETA(I,J)=ZETA(I,J)*HEFFM(I,J)
      END DO
      END DO

C NOW CALCULATE PRESSURE FORCE AND ADD TO EXTERNAL FORCE
C**** Update halo of PRESS for distributed memory implementation
      CALL CHECKSUM(grid, PRESS,  219, "ICEDYN.f")
      CALL HALO_UPDATE(grid, PRESS, FROM=NORTH)
      DO J=J_0,J_1S
        DO I=1,NX1-1
          FORCEX(I,J)=FORCEX(I,J)-(0.25/(DXU(I)*CSU(J)))
     1         *(PRESS(I+1,J)+PRESS(I+1,J+1)-PRESS(I,J)-PRESS(I,J+1))
          FORCEY(I,J)=FORCEY(I,J)-0.25/DYU(J)
     1         *(PRESS(I,J+1)+PRESS(I+1,J+1)-PRESS(I,J)-PRESS(I+1,J))
C NOW PUT IN MINIMAL MASS FOR TIME STEPPING CALCULATIONS
        END DO
      END DO

      DO J=J_0,J_1
        FORCEX(1,J)=FORCEX(NX1-1,J)
        FORCEY(1,J)=FORCEY(NX1-1,J)
        FORCEX(NX1,J)=FORCEX(2,J)
        FORCEY(NX1,J)=FORCEY(2,J)
        DWATN(1,J)=DWATN(NX1-1,J)
        DWATN(NX1,J)=DWATN(2,J)
      END DO

      RETURN
      END SUBROUTINE FORM


      SUBROUTINE PLAST 1,3
!@sum  PLAST Calculates strain rates and viscosity for dynamic ice
!@auth Jiping Liu/Gavin Schmidt (based on code from J. Zhang)
!@ver  1.0
      USE DOMAIN_DECOMP, only : grid, GET, NORTH,SOUTH
      USE DOMAIN_DECOMP, ONLY : HALO_UPDATE, CHECKSUM
      IMPLICIT NONE
      REAL*8, DIMENSION(NX1,grid%j_strt_halo:grid%j_stop_halo) 
     &        :: E11,E22,E12
c      REAL*8 :: SS11
      REAL*8, PARAMETER :: ECM2 = 1.0/(ECCEN**2),GMIN=1d-20
      REAL*8 DELT,DELT1,AAA
      INTEGER I,J

      INTEGER :: J_0,J_1,J_0S,J_1S

C****
C**** Extract useful local domain parameters from "grid"
C****
      CALL GET(grid, J_STRT     =J_0,    J_STOP     =J_1,
     &               J_STRT_SKP =J_0S,   J_STOP_SKP =J_1S )


C EVALUATE STRAIN RATES
      CALL CHECKSUM(grid, UICE, 267, "ICEDYN.f")
      CALL HALO_UPDATE(grid, UICE, FROM=SOUTH)
      CALL CHECKSUM(grid, VICE, 269, "ICEDYN.f")
      CALL HALO_UPDATE(grid, VICE, FROM=SOUTH)
      DO J=J_0S,J_1S
        DO I=2,NX1-1
          E11(I,J)=0.5/(DXT(I)*CST(J))*(UICE(I,J,1)+UICE(I,J-1,1)
     *         -UICE(I-1,J,1)-UICE(I-1,J-1,1))-0.25*(VICE(I,J,1)+
     *         VICE(I-1,J,1)+VICE(I-1,J-1,1)+VICE(I,J-1,1))*TNGT(J)
     *         /RADIUS
          E22(I,J)=0.5/DYT(J)*(VICE(I,J,1)+VICE(I-1,J,1)
     *         -VICE(I,J-1,1)-VICE(I-1,J-1,1))
          E12(I,J)=0.5*(0.5/DYT(J)*(UICE(I,J,1)+UICE(I-1,J,1)-
     *         UICE(I,J-1,1)-UICE(I-1,J-1,1))+0.5/(DXT(I)*CST(J))*
     *         (VICE(I,J,1)+VICE(I,J-1,1)-VICE(I-1,J,1)-VICE(I-1,J-1,1))
     *         +0.25*(UICE(I,J,1)+UICE(I-1,J,1)+UICE(I-1,J-1,1)+UICE(I,J
     *         -1,1))*TNGT(J)/RADIUS)
C NOW EVALUATE VISCOSITIES
          DELT=(E11(I,J)**2+E22(I,J)**2)*(1.0+ECM2)+4.0*ECM2*E12(I,J)**2
     *         +2.0*E11(I,J)*E22(I,J)*(1.0-ECM2)
          DELT1=SQRT(DELT)
          DELT1=MAX(GMIN,DELT1)
          ZETA(I,J)=0.5*PRESS(I,J)/DELT1
        END DO
      END DO

C NOW PUT MIN AND MAX VISCOSITIES IN
      DO J=J_0S,J_1S
        DO I=2,NX1-1
          ZETA(I,J)=MIN(ZMAX(I,J),ZETA(I,J))
          ZETA(I,J)=MAX(ZMIN(I,J),ZETA(I,J))
        END DO
      END DO

      if (grid%HAVE_NORTH_POLE) then
        AAA=0.0
        DO I=2,NX1-1
          AAA=AAA+ZETA(I,NY1-1)
        END DO
        AAA=AAA/FLOAT(NX1-2)
        DO I=1,NX1
          ZETA(I,NY1)=AAA
        END DO
      end if

      if (grid%HAVE_SOUTH_POLE) then
        AAA=0.0
        DO I=2,NX1-1
          AAA=AAA+ZETA(I,2)
        END DO
        AAA=AAA/FLOAT(NX1-2)
        DO I=1,NX1
          ZETA(I,1)=AAA
        END DO
      end if

      DO J=J_0,J_1
        ZETA(1,J)=ZETA(NX1-1,J)
        ZETA(NX1,J)=ZETA(2,J)
      END DO

      DO J=J_0,J_1
        DO I=1,NX1
          ETA(I,J)=ECM2*ZETA(I,J)
c         E11(I,J)=E11(I,J)*HEFFM(I,J)
c         E22(I,J)=E22(I,J)*HEFFM(I,J)
c         E12(I,J)=E12(I,J)*HEFFM(I,J)
c         SS11=(ZETA(I,J)-ETA(I,J))*(E11(I,J)+E22(I,J))-PRESS(I,J)*0.5
        END DO
      END DO
      RETURN
      END SUBROUTINE PLAST


      SUBROUTINE RELAX 2,7
!@sum  RELAX calculates ice dynamics relaxation method
!@auth Jiping Liu/Gavin Schmidt (based on code from J. Zhang)
!@ver  1.0
      USE DOMAIN_DECOMP, only : grid, GET, NORTH,SOUTH
      USE DOMAIN_DECOMP, ONLY : HALO_UPDATE, CHECKSUM
      IMPLICIT NONE

      REAL*8, DIMENSION(NX1,grid%J_STRT_HALO:grid%J_STOP_HALO) :: 
     &         AU,BU,CU,FXY,FXY1
      REAL*8, DIMENSION(grid%J_STRT_HALO:grid%J_STOP_HALO,NX1) :: 
     &         AV,BV,CV,FXYa,FXY1a
      REAL*8, DIMENSION(NX1) :: CUU,URT         !CUU,
      REAL*8, DIMENSION(grid%J_STRT_HALO:grid%J_STOP_HALO) :: 
     &         CVV,VRT,U_tmp                     !CVV,
      REAL*8, PARAMETER :: BYRAD2 = 1./(RADIUS*RADIUS)
      INTEGER I,J,J1,J2,IMD,JMD
      REAL*8 DELXY,DELXR,DELX2,DELY2,DELYR,ETAMEAN,ZETAMEAN,AA1,AA2
     *     ,AA3,AA4,AA5,AA6,AA9

      INTEGER :: J_0,J_1,J_0S,J_1S,J_0STG,J_1STG

C**** Replaces NYPOLE in loops. 
      INTEGER :: J_NYP 

C****
C**** Extract useful local domain parameters from "grid"
C****
      CALL GET(grid, J_STRT     =J_0,    J_STOP     =J_1,
     &               J_STRT_SKP =J_0S,   J_STOP_SKP =J_1S )

C****
C**** Modify if NYPOLE definition is modified.
C****
      J_NYP=J_1S

       DO J=J_0,J_1S
        DO I=1,NX1
          FORCEX(I,J)=FORCEX(I,J)*UVM(I,J)
          FORCEY(I,J)=FORCEY(I,J)*UVM(I,J)
        END DO
      END DO
C MUST UPDATE HEFF BEFORE CALLING RELAX
C FIRST SET U(2)=U(1)
       DO J=J_0,J_1S
        DO I=1,NX1
C NOW MAKE SURE BDRY PTS ARE EQUAL TO ZERO
          UICE(I,J,2)=UICE(I,J,1)
          VICE(I,J,2)=VICE(I,J,1)
          UICE(I,J,1)=UICE(I,J,3)*UVM(I,J)
          VICE(I,J,1)=VICE(I,J,3)*UVM(I,J)
        END DO
      END DO

      if (grid%HAVE_NORTH_POLE) then
        DO I=1,NX1/2
         UICE(I,NY1,1)=-UICEC(I+(NX1-2)/2,NY1-1)
         VICE(I,NY1,1)=-VICEC(I+(NX1-2)/2,NY1-1)
         UICE(I,NY1,3)=-UICEC(I+(NX1-2)/2,NY1-1)
         VICE(I,NY1,3)=-VICEC(I+(NX1-2)/2,NY1-1)
  
         UICEC(I,NY1)=-UICEC(I+(NX1-2)/2,NY1-1)
         VICEC(I,NY1)=-VICEC(I+(NX1-2)/2,NY1-1)
        END DO
  
        DO I=NX1/2+1,NX1-1
         UICE(I,NY1,1)=-UICEC(I-(NX1-2)/2,NY1-1)
         VICE(I,NY1,1)=-VICEC(I-(NX1-2)/2,NY1-1)
         UICE(I,NY1,3)=-UICEC(I-(NX1-2)/2,NY1-1)
         VICE(I,NY1,3)=-VICEC(I-(NX1-2)/2,NY1-1)
  
         UICEC(I,NY1)=-UICEC(I-(NX1-2)/2,NY1-1)
         VICEC(I,NY1)=-VICEC(I-(NX1-2)/2,NY1-1)
        END DO
      end if

      DO J=J_0,J_1S
        UICE(1,J,1)=UICEC(NX1-1,J)
        VICE(1,J,1)=VICEC(NX1-1,J)
        UICE(NX1,J,1)=UICEC(2,J)
        VICE(NX1,J,1)=VICEC(2,J)
        UICE(1,J,3)=UICEC(NX1-1,J)
        VICE(1,J,3)=VICEC(NX1-1,J)
        UICE(NX1,J,3)=UICEC(2,J)
        VICE(NX1,J,3)=VICEC(2,J)
  
        UICEC(1,J)=UICEC(NX1-1,J)
        VICEC(1,J)=VICEC(NX1-1,J)
        UICEC(NX1,J)=UICEC(2,J)
        VICEC(NX1,J)=VICEC(2,J)
      END DO

C FIRST DO UICE
C THE FIRST HALF

C**Update halos for arrays eta,zeta,vicec,bycsu as needed in the next loop
      CALL CHECKSUM(grid, ETA, 436, "ICEDYN.f")
        CALL HALO_UPDATE(grid, ETA, FROM=NORTH)
      CALL CHECKSUM(grid, ZETA, 438, "ICEDYN.f")
        CALL HALO_UPDATE(grid, ZETA, FROM=NORTH)
      CALL CHECKSUM(grid, VICEC, 440, "ICEDYN.f")
        CALL HALO_UPDATE(grid, VICEC, FROM=NORTH)
        CALL HALO_UPDATE(grid, VICEC, FROM=SOUTH)
      CALL CHECKSUM(grid, BYCSU, 443, "ICEDYN.f")
        CALL HALO_UPDATE(grid, BYCSU, FROM=NORTH)
        CALL HALO_UPDATE(grid, BYCSU, FROM=SOUTH)

      DO J=J_0S,J_NYP
      DO I=2,NXLCYC
      DELXY=BYDXDY(I,J)    ! 0.5/(DXU(I)*DYU(J))
      DELXR=BYDXR(I)       ! 0.5/(DXU(I)*RADIUS)
      DELX2=BYDX2(I)       ! 0.5/(DXU(I)*DXU(I))
      ETAMEAN=0.25*(ETA(I,J+1)+ETA(I+1,J+1)+ETA(I,J)+ETA(I+1,J))
      ZETAMEAN=0.25*(ZETA(I,J+1)+ZETA(I+1,J+1)+ZETA(I,J)+ZETA(I+1,J))

      FXY(I,J)=DRAGA(I,J)*VICEC(I,J)+FORCEX(I,J)
     3+0.5*(ZETA(I+1,J+1)*(VICEC(I+1,J+1)+VICEC(I,J+1)
     3-VICEC(I+1,J)-VICEC(I,J))+ZETA(I+1,J)*(VICEC(I+1,J)
     3+VICEC(I,J)-VICEC(I+1,J-1)-VICEC(I,J-1))+ZETA(I,J+1)
     3*(VICEC(I,J)+VICEC(I-1,J)-VICEC(I,J+1)-VICEC(I-1,J+1))
     3+ZETA(I,J)*(VICEC(I,J-1)+VICEC(I-1,J-1)-VICEC(I,J)
     3-VICEC(I-1,J)))*DELXY*BYCSU(J)
     3
     4-0.5*(ETA(I+1,J+1)*(VICEC(I+1,J+1)+VICEC(I,J+1)
     4-VICEC(I+1,J)-VICEC(I,J))+ETA(I+1,J)*(VICEC(I+1,J)
     4+VICEC(I,J)-VICEC(I+1,J-1)-VICEC(I,J-1))+ETA(I,J+1)
     4*(VICEC(I,J)+VICEC(I-1,J)-VICEC(I,J+1)-VICEC(I-1,J+1))
     4+ETA(I,J)*(VICEC(I,J-1)+VICEC(I-1,J-1)-VICEC(I,J)
     4-VICEC(I-1,J)))*DELXY*BYCSU(J)
     4
     5+0.5*(VICEC(I+1,J)-VICEC(I-1,J))*(ETA(I,J+1)+ETA(I+1,J+1)
     5-ETA(I,J)-ETA(I+1,J))*DELXY*BYCSU(J)+0.5*ETAMEAN*((VICEC(I+1,J+1)
     5-VICEC(I-1,J+1))*BYCSU(J+1)-(VICEC(I+1,J-1)-VICEC(I-1,J-1))
     5*BYCSU(J-1))*DELXY
     5
     6-((ZETA(I+1,J+1)+ZETA(I+1,J)-ZETA(I,J)-ZETA(I,J+1))
     6+(ETA(I+1,J+1)+ETA(I+1,J)-ETA(I,J)-ETA(I,J+1)))
     6*TNG(J)*VICEC(I,J)*DELXR*BYCSU(J)
     6-(ETAMEAN+ZETAMEAN)*TNG(J)*(VICEC(I+1,J)-VICEC(I-1,J))
     6*DELXR*BYCSU(J)
     6
     7-ETAMEAN*2.0*TNG(J)*(VICEC(I+1,J)-VICEC(I-1,J))*DELXR*BYCSU(J)

      END DO
      END DO

      DO J=J_0S,J_NYP
      DO I=2,NXLCYC
      DELX2=BYDX2(I)       ! 0.5/(DXU(I)*DXU(I))
      DELY2=BYDY2(J)       ! 0.5/(DYU(J)*DYU(J))
      DELYR=BYDYR(J)       ! 0.5/(DYU(J)*RADIUS)
      ETAMEAN=0.25*(ETA(I,J+1)+ETA(I+1,J+1)+ETA(I,J)+ETA(I+1,J))
      AA1=((ETA(I+1,J)  +ZETA(I+1,J)  )*BYCSU(J)+
     *     (ETA(I+1,J+1)+ZETA(I+1,J+1))*BYCSU(J))*BYCSU(J)
      AA2=((ETA(I,J)+ZETA(I,J))*BYCSU(J)+(ETA(I,J+1)+ZETA(I,J+1))
     &     *BYCSU(J))*BYCSU(J)
      AA3=ETA(I,J+1)+ETA(I+1,J+1)
      AA4=ETA(I,J)+ETA(I+1,J)
      AA5=-(ETA(I,J+1)+ETA(I+1,J+1)-ETA(I,J)-ETA(I+1,J))*TNG(J)
      AA6=2.0*ETAMEAN*TNG(J)*TNG(J)
      AU(I,J)=-AA2*DELX2*UVM(I,J)
      BU(I,J)=((AA1+AA2)*DELX2+AA6*BYRAD2
     &+AMASS(I,J)*BYDTS*2.0+DRAGS(I,J))*UVM(I,J)+(1.0-UVM(I,J))
      CU(I,J)=-AA1*DELX2*UVM(I,J)
      END DO
      END DO
      DO J=J_0S,J_NYP
      AU(2,J)=0.0
      CU(NXLCYC,J)=0.0
c      CU(2,J)=CU(2,J)/BU(2,J)  ! absorbed into TRIDIAG
      END DO

C**Update halos for UICE and TNG as needed for loop 1200
C**(ETA and ZETA were updted above)
      CALL CHECKSUM(grid, UICE, 514, "ICEDYN.f")
        CALL HALO_UPDATE(grid, UICE, FROM=SOUTH)
        CALL HALO_UPDATE(grid, UICE, FROM=NORTH)
      CALL CHECKSUM(grid, TNG, 517, "ICEDYN.f")
        CALL HALO_UPDATE(grid, TNG, FROM=SOUTH)
        CALL HALO_UPDATE(grid, TNG, FROM=NORTH)

      DO 1200 J=J_0S,J_NYP
      DO I=2,NXLCYC

      DELX2=BYDX2(I)       ! 0.5/(DXU(I)*DXU(I))
      DELY2=BYDY2(J)       ! 0.5/(DYU(J)*DYU(J))
      DELYR=BYDYR(J)       ! 0.5/(DYU(J)*RADIUS)
      ETAMEAN=0.25*(ETA(I,J+1)+ETA(I+1,J+1)+ETA(I,J)+ETA(I+1,J))

      AA1=((ETA(I+1,J)  +ZETA(I+1,J)  )*BYCSU(J)+
     *     (ETA(I+1,J+1)+ZETA(I+1,J+1))*BYCSU(J))*BYCSU(J)
      AA2=((ETA(I,J)+ZETA(I,J))*BYCSU(J)+(ETA(I,J+1)+ZETA(I,J+1))
     &     *BYCSU(J))*BYCSU(J)
      AA3=ETA(I,J+1)+ETA(I+1,J+1)
      AA4=ETA(I,J)+ETA(I+1,J)
      AA5=-(ETA(I,J+1)+ETA(I+1,J+1)-ETA(I,J)-ETA(I+1,J))*TNG(J)
      AA6=2.0*ETAMEAN*TNG(J)*TNG(J)

      IF(I.EQ.2) THEN
      AA9=AA2*DELX2*UICEC(I-1,J)*UVM(I,J)*FLOAT(LCYC-0)
      ELSE IF(I.EQ.NXLCYC) THEN
      AA9=AA1*DELX2*UICEC(I+1,J)*UVM(I,J)*FLOAT(LCYC-0)
      ELSE
      AA9=0.0
      END IF

      URT(I)=AA9+FXY(I,J)-AA5*DELYR*UICE(I,J,2)
     1-(AA3+AA4)*DELY2*UICE(I,J,2)
     1+(ETA(I,J+1)+ETA(I+1,J+1))*UICE(I,J+1,2)*DELY2
     2+(ETA(I,J)+ETA(I+1,J))*UICE(I,J-1,2)*DELY2
     3+ETAMEAN*DELYR*(UICE(I,J+1,2)*TNG(J+1)-UICE(I,J-1,2)*TNG(J-1))
     4-ETAMEAN*DELYR*2.0*TNG(J)*(UICE(I,J+1,2)-UICE(I,J-1,2))
      URT(I)=(URT(I)+AMASS(I,J)*BYDTS*UICE(I,J,2)*2.0)*UVM(I,J)
      END DO

      CALL TRIDIAG(AU(2,J),BU(2,J),CU(2,J),URT(2),UICE(2,J,1),NXLCYC-1)

c      DO I=2,NXLCYC
c       CUU(I)=CU(I,J)
c      END DO
c      URT(2)=URT(2)/BU(2,J)
c      DO I=3,NXLCYC
c       IMD=I-1
c       CUU(I)=CUU(I)/(BU(I,J)-AU(I,J)*CUU(IMD))
c       URT(I)=(URT(I)-AU(I,J)*URT(IMD))/(BU(I,J)-AU(I,J)*CUU(IMD))
c      END DO
c      DO I=1,NXLCYC-2
c       J1=NXLCYC-I
c       J2=J1+1
c       URT(J1)=URT(J1)-CUU(J1)*URT(J2)
c      END DO
c      DO I=2,NXLCYC
c       UICE(I,J,1)=URT(I)
c      END DO

 1200 CONTINUE

      DO I=2,NXLCYC
      DO J=J_0S,J_NYP
      UICE(I,J,3)=UICE(I,J,1)
      END DO
      END DO

C**ETA, ZETA, and TNG halos already updated above**
C NOW THE SECOND HALF
      DO I=2,NXLCYC
      DO J=J_0S,J_NYP
       DELX2=BYDX2(I)       ! 0.5/(DXU(I)*DXU(I))
       DELY2=BYDY2(J)       ! 0.5/(DYU(J)*DYU(J))
       DELYR=BYDYR(J)       ! 0.5/(DYU(J)*RADIUS)
       ETAMEAN=0.25*(ETA(I,J+1)+ETA(I+1,J+1)+ETA(I,J)+ETA(I+1,J))
       ZETAMEAN=0.25*(ZETA(I,J+1)+ZETA(I+1,J+1)+ZETA(I,J)+ZETA(I+1,J))

       AA1=ETA(I,J+1)+ETA(I+1,J+1)
       AA2=ETA(I,J)+ETA(I+1,J)
       AA5=-(ETA(I,J+1)+ETA(I+1,J+1)-ETA(I,J)-ETA(I+1,J))*TNG(J)
       AA6=2.0*ETAMEAN*TNG(J)*TNG(J)

       AV(J,I)=(-AA2*DELY2+ETAMEAN*DELYR*(TNG(J-1)-2.0*TNG(J)))*UVM(I,J)
       BV(J,I)=((AA1+AA2)*DELY2+AA5*DELYR+AA6*BYRAD2
     & +AMASS(I,J)*BYDTS*2.0+DRAGS(I,J))*UVM(I,J)+(1.0-UVM(I,J))
       CV(J,I)=(-AA1*DELY2-ETAMEAN*DELYR*(TNG(J+1)-2.0*TNG(J)))*UVM(I,J)
      END DO
      END DO

      if (grid%HAVE_SOUTH_POLE) then
        DO I=2,NXLCYC
          AV(2,I)=0.0
c         CV(2,I)=CV(2,I)/BV(2,I)  ! absorbed into TRIDIAG
        END DO
      end if

      if (grid%HAVE_NORTH_POLE) then
        DO I=2,NXLCYC
          CV(NYPOLE,I)=0.0
        END DO
      end if

      DO I=2,NXLCYC
      DO J=J_0S,J_NYP
        DELX2=BYDX2(I)       ! 0.5/(DXU(I)*DXU(I))
        DELY2=BYDY2(J)       ! 0.5/(DYU(J)*DYU(J))
        DELYR=BYDYR(J)       ! 0.5/(DYU(J)*RADIUS)
        ETAMEAN=0.25*(ETA(I,J+1)+ETA(I+1,J+1)+ETA(I,J)+ETA(I+1,J))
  
        AA1=((ETA(I+1,J)  +ZETA(I+1,J)  )*BYCSU(J)+
     *      (ETA(I+1,J+1)+ZETA(I+1,J+1))*BYCSU(J))*BYCSU(J)
        AA2=((ETA(I,J)+ZETA(I,J))*BYCSU(J)+(ETA(I,J+1)+ZETA(I,J+1))
     &      *BYCSU(J))*BYCSU(J)

        IF(J.EQ.NYPOLE) THEN
          AA9=(     &    +ETAMEAN*DELYR*(TNG(J+1)-2.0*TNG(J))*UICEC(I,J+1) )*UVM(I,J)
     &       *FLOAT(NPOL-0)
        ELSE
          AA9=0.0
        END IF

        FXY1a(J,I)=AA9+AMASS(I,J)*BYDTS*UICE(I,J,1)*2.0
     5  -(AA1+AA2)*DELX2*UICE(I,J,1)
     6  +((ETA(I+1,J)+ZETA(I+1,J)+ETA(I+1,J+1)+ZETA(I+1,J+1))
     6  *UICE(I+1,J,1)
     6  +(ETA(I,J)+ZETA(I,J)+ETA(I,J+1)+ZETA(I,J+1))*UICE(I-1,J,1))
     6  *DELX2*BYCSU(J)*BYCSU(J)

      END DO
      END DO

      DO 1300 I=2,NXLCYC
      DO J=J_0S,J_NYP
      VRT(J)=FXY(I,J)+FXY1a(J,I)
      VRT(J)=VRT(J)*UVM(I,J)
      END DO

      CALL TRIDIAG(AV(2,I),BV(2,I),CV(2,I),VRT(2),U_tmp(2),NYPOLE-1)

c      DO J=J_0S,J_NYP
c      CVV(J)=CV(J,I)
c      END DO
c      VRT(2)=VRT(2)/BV(2,I)
c      DO J=3,NYPOLE
c      JMD=J-1
c      CVV(J)=CVV(J)/(BV(J,I)-AV(J,I)*CVV(JMD))
c      VRT(J)=(VRT(J)-AV(J,I)*VRT(JMD))/(BV(J,I)-AV(J,I)*CVV(JMD))
c      END DO
c      DO J=1,NYPOLE-2
c      J1=NYPOLE-J
c      J2=J1+1
c      VRT(J1)=VRT(J1)-CVV(J1)*VRT(J2)
c      END DO
      DO J=J_0S,J_NYP
        UICE(I,J,1)=U_tmp(J)     ! VRT(J)   !
      END DO
 1300 CONTINUE

C NOW DO VICE
C THE FIRST HALF

      CALL CHECKSUM(grid, UICEC, 678, "ICEDYN.f")
      CALL HALO_UPDATE(grid, UICEC, FROM=NORTH)

      DO I=2,NXLCYC
      DO J=J_0S,J_NYP
      DELXY=BYDXDY(I,J)    ! 0.5/(DXU(I)*DYU(J))
      DELXR=BYDXR(I)       ! 0.5/(DXU(I)*RADIUS)
      DELY2=BYDY2(J)       ! 0.5/(DYU(J)*DYU(J))
      DELYR=BYDYR(J)       ! 0.5/(DYU(J)*RADIUS)
      ETAMEAN=0.25*(ETA(I,J+1)+ETA(I+1,J+1)+ETA(I,J)+ETA(I+1,J))
      ZETAMEAN=0.25*(ZETA(I,J+1)+ZETA(I+1,J+1)+ZETA(I,J)+ZETA(I+1,J))

      FXYa(J,I)=-DRAGA(I,J)*UICEC(I,J)+FORCEY(I,J)
     3+(0.5*(UICEC(I+1,J)-UICEC(I-1,J))*(ZETA(I,J+1)+ZETA(I+1,J+1)
     3-ZETA(I,J)-ZETA(I+1,J))*DELXY*BYCSU(J)+0.5*ZETAMEAN*
     3((UICEC(I+1,J+1)
     3-UICEC(I-1,J+1))*BYCSU(J+1)-(UICEC(I+1,J-1)-UICEC(I-1,J-1))
     3*BYCSU(J-1))*DELXY)
     3
     4-(0.5*(UICEC(I+1,J)-UICEC(I-1,J))*(ETA(I,J+1)+ETA(I+1,J+1)
     4-ETA(I,J)-ETA(I+1,J))*DELXY*BYCSU(J)+0.5*ETAMEAN*((UICEC(I+1,J+1)
     4-UICEC(I-1,J+1))*BYCSU(J+1)-(UICEC(I+1,J-1)-UICEC(I-1,J-1))
     4*BYCSU(J-1))*DELXY)
     4
     5+0.5*(ETA(I+1,J+1)*(UICEC(I+1,J+1)+UICEC(I,J+1)
     5-UICEC(I+1,J)-UICEC(I,J))+ETA(I+1,J)*(UICEC(I+1,J)
     5+UICEC(I,J)-UICEC(I+1,J-1)-UICEC(I,J-1))+ETA(I,J+1)
     5*(UICEC(I,J)+UICEC(I-1,J)-UICEC(I,J+1)-UICEC(I-1,J+1))
     5+ETA(I,J)*(UICEC(I,J-1)+UICEC(I-1,J-1)-UICEC(I,J)
     5-UICEC(I-1,J)))*DELXY*BYCSU(J)
     5
     6+(ETA(I+1,J+1)+ETA(I+1,J)-ETA(I,J)-ETA(I,J+1))
     6*TNG(J)*UICEC(I,J)*DELXR*BYCSU(J)
     6+ETAMEAN*TNG(J)*(UICEC(I+1,J)-UICEC(I-1,J))*DELXR*BYCSU(J)
     6
     7+ETAMEAN*2.0*TNG(J)*(UICEC(I+1,J)-UICEC(I-1,J))*DELXR*BYCSU(J)

      END DO
      END DO

      DO I=2,NXLCYC
      DO J=J_0S,J_NYP
      DELX2=BYDX2(I)       ! 0.5/(DXU(I)*DXU(I))
      DELY2=BYDY2(J)       ! 0.5/(DYU(J)*DYU(J))
      DELYR=BYDYR(J)       ! 0.5/(DYU(J)*RADIUS)
      ETAMEAN=0.25*(ETA(I,J+1)+ETA(I+1,J+1)+ETA(I,J)+ETA(I+1,J))
      ZETAMEAN=0.25*(ZETA(I,J+1)+ZETA(I+1,J+1)+ZETA(I,J)+ZETA(I+1,J))
      AA1=ETA(I,J+1)+ZETA(I,J+1)+ETA(I+1,J+1)+ZETA(I+1,J+1)
      AA2=ETA(I,J)+ZETA(I,J)+ETA(I+1,J)+ZETA(I+1,J)
      AA3=(ETA(I+1,J)*BYCSU(J)+ETA(I+1,J+1)*BYCSU(J))*BYCSU(J)
      AA4=(ETA(I,J)*BYCSU(J)+ETA(I,J+1)*BYCSU(J))*BYCSU(J)
      AA5=((ZETA(I,J+1)-ETA(I,J+1))+(ZETA(I+1,J+1)-ETA(I+1,J+1))
     &-(ZETA(I,J)-ETA(I,J))-(ZETA(I+1,J)-ETA(I+1,J)))*TNG(J)
      AA6=2.0*ETAMEAN*TNG(J)*TNG(J)

      AV(J,I)=(-AA2*DELY2-(ZETAMEAN-ETAMEAN)*TNG(J-1)*DELYR
     &-ETAMEAN*2.0*TNG(J)*DELYR)*UVM(I,J)
      BV(J,I)=((AA1+AA2)*DELY2+AA5*DELYR+AA6*BYRAD2
     &+AMASS(I,J)*BYDTS*2.0+DRAGS(I,J))*UVM(I,J)+(1.0-UVM(I,J))
      CV(J,I)=(-AA1*DELY2+(ZETAMEAN-ETAMEAN)*TNG(J+1)*DELYR
     &+ETAMEAN*2.0*TNG(J)*DELYR)*UVM(I,J)
      END DO
      END DO

      if (grid%HAVE_SOUTH_POLE) then
        DO I=2,NXLCYC
          AV(2,I)=0.0
c         CV(2,I)=CV(2,I)/BV(2,I)  ! absorbed into TRIDIAG
        END DO
      end if

      if (grid%HAVE_NORTH_POLE) then
        DO I=2,NXLCYC
          CV(NYPOLE,I)=0.0
        END DO
      end if

      DO 1301 I=2,NXLCYC
      DO J=J_0S,J_NYP
        DELX2=BYDX2(I)       ! 0.5/(DXU(I)*DXU(I))
        DELY2=BYDY2(J)       ! 0.5/(DYU(J)*DYU(J))
        DELYR=BYDYR(J)       ! 0.5/(DYU(J)*RADIUS)

        ETAMEAN=0.25*(ETA(I,J+1)+ETA(I+1,J+1)+ETA(I,J)+ETA(I+1,J))
        ZETAMEAN=0.25*(ZETA(I,J+1)+ZETA(I+1,J+1)+ZETA(I,J)+ZETA(I+1,J))

        AA1=ETA(I,J+1)+ZETA(I,J+1)+ETA(I+1,J+1)+ZETA(I+1,J+1)
        AA2=ETA(I,J)+ZETA(I,J)+ETA(I+1,J)+ZETA(I+1,J)
        AA3=(ETA(I+1,J)*BYCSU(J)+ETA(I+1,J+1)*BYCSU(J))*BYCSU(J)
        AA4=(ETA(I,J)*BYCSU(J)+ETA(I,J+1)*BYCSU(J))*BYCSU(J)
        AA5=((ZETA(I,J+1)-ETA(I,J+1))+(ZETA(I+1,J+1)-ETA(I+1,J+1))
     &  -(ZETA(I,J)-ETA(I,J))-(ZETA(I+1,J)-ETA(I+1,J)))*TNG(J)
        AA6=2.0*ETAMEAN*TNG(J)*TNG(J)

        IF(J.EQ.NYPOLE) THEN
          AA9=(AA1*DELY2-(ZETAMEAN-ETAMEAN)*TNG(J+1)*DELYR
     &    -ETAMEAN*2.0*TNG(J)*DELYR)*VICEC(I,J+1)*UVM(I,J)*FLOAT(NPOL-0)
        ELSE
          AA9=0.0
        END IF

        VRT(J)=AA9+FXYa(J,I)-(AA3+AA4)*DELX2*VICE(I,J,2)
     6 +((ETA(I+1,J)*BYCSU(J)+ETA(I+1,J+1)*BYCSU(J))*VICE(I+1,J,2)*DELX2
     7    +(ETA(I,J)*BYCSU(J)+ETA(I,J+1)*BYCSU(J))*VICE(I-1,J,2)*DELX2)
     *       *BYCSU(J)
        VRT(J)=(VRT(J)+AMASS(I,J)*BYDTS*VICE(I,J,2)*2.0)*UVM(I,J)
      END DO

      CALL TRIDIAG(AV(2,I),BV(2,I),CV(2,I),VRT(2),U_tmp(2),NYPOLE-1)

c      DO J=J_0S,J_NYP
c      CVV(J)=CV(J,I)
c      END DO
c      VRT(2)=VRT(2)/BV(2,I)
c      DO J=3,NYPOLE
c      JMD=J-1
c      CVV(J)=CVV(J)/(BV(J,I)-AV(J,I)*CVV(JMD))
c      VRT(J)=(VRT(J)-AV(J,I)*VRT(JMD))/(BV(J,I)-AV(J,I)*CVV(JMD))
c      END DO
c      DO J=1,NYPOLE-2
c      J1=NYPOLE-J
c      J2=J1+1
c      VRT(J1)=VRT(J1)-CVV(J1)*VRT(J2)
c      END DO
      DO J=J_0S,J_NYP
        VICE(I,J,1)=U_tmp(J)   ! VRT(J)   !
      END DO
 1301 CONTINUE

      DO I=2,NXLCYC
       DO J=J_0S,J_NYP
        VICE(I,J,3)=VICE(I,J,1)
       END DO
      END DO

C NOW THE SECOND HALF

      DO J=J_0S,J_NYP
       DO I=2,NXLCYC
        DELX2=BYDX2(I)       ! 0.5/(DXU(I)*DXU(I))
        DELY2=BYDY2(J)       ! 0.5/(DYU(J)*DYU(J))
        DELYR=BYDYR(J)       ! 0.5/(DYU(J)*RADIUS)
        ETAMEAN=0.25*(ETA(I,J+1)+ETA(I+1,J+1)+ETA(I,J)+ETA(I+1,J))
        ZETAMEAN=0.25*(ZETA(I,J+1)+ZETA(I+1,J+1)+ZETA(I,J)+ZETA(I+1,J))

        AA1=ETA(I,J+1)+ZETA(I,J+1)+ETA(I+1,J+1)+ZETA(I+1,J+1)
        AA2=ETA(I,J)+ZETA(I,J)+ETA(I+1,J)+ZETA(I+1,J)
        AA3=(ETA(I+1,J)*BYCSU(J)+ETA(I+1,J+1)*BYCSU(J))*BYCSU(J)
        AA4=(ETA(I,J)*BYCSU(J)+ETA(I,J+1)*BYCSU(J))*BYCSU(J)
        AA6=2.0*ETAMEAN*TNG(J)*TNG(J)

        AU(I,J)=-AA4*DELX2*UVM(I,J)
        BU(I,J)=((AA3+AA4)*DELX2+AA6*BYRAD2
     &  +AMASS(I,J)*BYDTS*2.0+DRAGS(I,J))*UVM(I,J)+(1.0-UVM(I,J))
        CU(I,J)=-AA3*DELX2*UVM(I,J)
       END DO
      END DO

      DO J=J_0S,J_NYP
       AU(2,J)=0.0
       CU(NXLCYC,J)=0.0
c       CU(2,J)=CU(2,J)/BU(2,J)   ! absorbed into TRIDIAG
      END DO

      CALL CHECKSUM(grid, VICE, 842, "ICEDYN.f")
      CALL HALO_UPDATE(grid, VICE, FROM=SOUTH)
      CALL HALO_UPDATE(grid, VICE, FROM=NORTH)

      DO J=J_0S,J_NYP
      DO I=2,NXLCYC

        DELX2=BYDX2(I)       ! 0.5/(DXU(I)*DXU(I))
        DELY2=BYDY2(J)       ! 0.5/(DYU(J)*DYU(J))
        DELYR=BYDYR(J)       ! 0.5/(DYU(J)*RADIUS)
        ETAMEAN=0.25*(ETA(I,J+1)+ETA(I+1,J+1)+ETA(I,J)+ETA(I+1,J))
        ZETAMEAN=0.25*(ZETA(I,J+1)+ZETA(I+1,J+1)+ZETA(I,J)+ZETA(I+1,J))

        AA1=ETA(I,J+1)+ZETA(I,J+1)+ETA(I+1,J+1)+ZETA(I+1,J+1)
        AA2=ETA(I,J)+ZETA(I,J)+ETA(I+1,J)+ZETA(I+1,J)
        AA3=(ETA(I+1,J)*BYCSU(J)+ETA(I+1,J+1)*BYCSU(J))*BYCSU(J)
        AA4=(ETA(I,J)*BYCSU(J)+ETA(I,J+1)*BYCSU(J))*BYCSU(J)
        AA5=((ZETA(I,J+1)-ETA(I,J+1))+(ZETA(I+1,J+1)-ETA(I+1,J+1))
     &  -(ZETA(I,J)-ETA(I,J))-(ZETA(I+1,J)-ETA(I+1,J)))*TNG(J)
        AA6=2.0*ETAMEAN*TNG(J)*TNG(J)

        IF(I.EQ.2) THEN
          AA9=AA4*DELX2*VICEC(I-1,J)*UVM(I,J)*FLOAT(LCYC-0)
        ELSE IF(I.EQ.NXLCYC) THEN
          AA9=AA3*DELX2*VICEC(I+1,J)*UVM(I,J)*FLOAT(LCYC-0)
        ELSE
          AA9=0.0
        END IF
        FXY1(I,J)=AA9+AMASS(I,J)*BYDTS*VICE(I,J,1)*2.0
     1  -AA5*DELYR*VICE(I,J,1)
     1  -(AA1+AA2)*DELY2*VICE(I,J,1)
     1  +AA1*DELY2*VICE(I,J+1,1)-((ZETAMEAN-ETAMEAN)*TNG(J+1)*DELYR
     1  +ETAMEAN*2.0*TNG(J)*DELYR)*VICE(I,J+1,1)
     2  +AA2*DELY2*VICE(I,J-1,1)+((ZETAMEAN-ETAMEAN)*TNG(J-1)*DELYR
     2  +ETAMEAN*2.0*TNG(J)*DELYR)*VICE(I,J-1,1)

      END DO
      END DO

      DO 1201 J=J_0S,J_NYP
        DO I=2,NXLCYC
          URT(I)=FXYa(J,I)+FXY1(I,J)
          URT(I)=URT(I)*UVM(I,J)
        END DO

        CALL TRIDIAG(AU(2,J),BU(2,J),CU(2,J),URT(2),VICE(2,J,1),
     &               NXLCYC-1)

c        DO I=2,NXLCYC
c          CUU(I)=CU(I,J)
c        END DO
c        URT(2)=URT(2)/BU(2,J)
c        DO I=3,NXLCYC
c          IMD=I-1
c          CUU(I)=CUU(I)/(BU(I,J)-AU(I,J)*CUU(IMD))
c          URT(I)=(URT(I)-AU(I,J)*URT(IMD))/(BU(I,J)-AU(I,J)*CUU(IMD))
c        END DO
c        DO I=1,NXLCYC-2
c          J1=NXLCYC-I
c          J2=J1+1
c          URT(J1)=URT(J1)-CUU(J1)*URT(J2)
c        END DO
c        DO I=2,NXLCYC
c          VICE(I,J,1)=URT(I)
c        END DO
 1201 CONTINUE

      DO J=J_0S,J_NYP
        DO I=2,NXLCYC
          UICE(I,J,1)=UICE(I,J,1)*UVM(I,J)
          VICE(I,J,1)=VICE(I,J,1)*UVM(I,J)
        END DO
      END DO

      RETURN
      END SUBROUTINE RELAX


      SUBROUTINE setup_icedyn_grid 1,4
      USE MODEL_COM, only : dts=>dtsrc
      USE DOMAIN_DECOMP, only : grid, GET, NORTH,SOUTH
      USE DOMAIN_DECOMP, ONLY : HALO_UPDATE, CHECKSUM
      IMPLICIT NONE
      REAL*8 :: dlat,dlon,phit,phiu,hemi,rms
      INTEGER I,J,n,k,kki,sumk,l

      INTEGER :: J_0,J_1,J_0S,J_1S

C****
C**** Extract useful local domain parameters from "grid"
C****
      CALL GET(grid_NXY, J_STRT     =J_0,    J_STOP     =J_1,
     &               J_STRT_SKP =J_0S,   J_STOP_SKP =J_1S )
C****
C**** calculate grid and initialise arrays
c****
      dlat=nint(180./(ny1-1))*radian
      dlon=nint(360./(nx1-2))*radian
      bydts = 1./dts
c****
      do j = j_0,j_1
        dyt(j) = dlat*radius
        dyu(j) = dlat*radius
      enddo
      do i=1,nx1
        dxt(i) = dlon*radius
        dxu(i) = dlon*radius
      enddo
      dxt(1) = dxt(nx1-1)
      dxt(nx1) = dxt(2)
      dxu(1) = dxu(nx1-1)
      dxu(nx1) = dxu(2)

      do i=1,nx1
        bydx2(i)=0.5/(dxu(i)*dxu(i))
        bydxr(i)=0.5/(dxu(i)*radius)
      end do
       do j=j_0,j_1
        bydy2(j)=0.5/(dyu(j)*dyu(j))
        bydyr(j)=0.5/(dyu(j)*radius)
        do i=1,nx1
          bydxdy(i,j) = 0.5/(dxu(i)*dyu(j))
        end do
      end do

      do j = j_0,j_1
       phit = (-90.+(j-1)*4.)*radian
       phiu = (-88.+(j-1)*4.)*radian
       cst(j) = cos(phit)
       csu(j) = cos(phiu)
       bycsu(j) = 1./csu(j)
       tng(j) = sin(phiu)/csu(j)
       TNGT(J)=SIN(PHIT)/CST(J)
       DO I=1,NX1
         SINEN(I,J)=SIN(PHIU)
       enddo
      enddo
      if (grid%have_north_pole) then
        TNGT(NY1)=TNGT(NY1-1)
        TNG(NY1)=TNG(NY1-1)
        CSU(NY1)=CSU(NY1-1)
        bycsu(NY1) = 1./csu(NY1)
      end if

C**** sin/cos ice-ocean turning angle
      SINWAT=SIN(OIPHI)
      COSWAT=COS(OIPHI)

C**** Set land masks for tracer and velocity points
       do j=j_0,j_1
        do i=2,nx1-1
          heffm(i,j)=nint(focean(i-1,j))
        enddo
        heffm(1,j)=heffm(nx1-1,j)
        heffm(nx1,j)=heffm(2,j)  
      enddo
C**** define velocity points (including exterior corners)
      CALL CHECKSUM(grid, HEFFM, 998, "ICEDYN.f")
      CALL HALO_UPDATE(grid, HEFFM, FROM=NORTH)
      do j=j_0,j_1s
        do i=1,nx1-1
c          sumk=heffm(i,j)+heffm(i+1,j)+heffm(i,j+1)+heffm(i+1,j+1)
c          if (sumk.ge.3) uvm(i,j)=1  ! includes exterior corners
          uvm(i,j) = nint(min(heffm(i,j), heffm(i+1,j), heffm(i,j+1),
     *         heffm(i+1,j+1)))
        end do
      end do
C**** reset tracer points to surround velocity points (except for single
      CALL CHECKSUM(grid, UVM, 1009, "ICEDYN.f")
      CALL HALO_UPDATE(grid, UVM, FROM=SOUTH)
c     CALL HALO_UPDATE(grid, UVM, FROM=NORTH)
      do j=j_0s,j_1s
        do i=2,nx1-1
          k = nint(max (uvm(i,j), uvm(i-1,j), uvm(i,j-1), uvm(i-1,j-1)))
c         sumk = nint(uvm(i,j)+uvm(i+1,j)+uvm(i,j+1)+uvm(i+1,j+1))
c set to k except if an island 
c         if (.not. (sumk.eq.4.and.focean(i-1,j).eq.0) ) then
            heffm(i,j) = k
c         end if
        enddo
      enddo
C**** final sweep to reinstate islands
c      do j=j_0s,j_1s
c        do i=2,nx1-1
c          sumk = nint(uvm(i,j)+uvm(i+1,j)+uvm(i,j+1)+uvm(i+1,j+1))
c          if (sumk.eq.4.and.heffm(i,j).eq.0.) then
c            uvm(i,j)=0 ; uvm(i+1,j)=0 ; uvm(i,j+1)=0 ; uvm(i+1,j+1)=0
c          end if
c        enddo
c      enddo
c set lateral boundary conditions
       do j=j_0,j_1
        heffm(1,j)   = heffm(nx1-1,j)
        heffm(nx1,j) = heffm(2,j)
      enddo

C**** Update halo of PHI for distributed memory implementation
      CALL CHECKSUM(grid, HEFFM, 1038, "ICEDYN.f")
      CALL HALO_UPDATE(grid, HEFFM, FROM=NORTH)
      do j=j_0,j_1s
        do i=1,nx1-1
          uvm(i,j) = nint(min(heffm(i,j), heffm(i+1,j), heffm(i,j+1),
     *         heffm(i+1,j+1)))
        end do
      end do

c set cyclic conditions on eastern and western boundary
       do j=j_0,j_1S
        uvm(1,j) = uvm(nx1-1,j)
        uvm(nx1,j) = uvm(2,j)
      enddo

      RETURN
      END SUBROUTINE setup_icedyn_grid

      END MODULE ICEDYN


      SUBROUTINE VPICEDYN 1,8
!@sum  vpicedyn is the entry point into the viscous-plastic ice 
!@+    dynamics code
!@auth Gavin Schmidt (based on code from J. Zhang)
      USE DOMAIN_DECOMP, only : grid, GET, NORTH,SOUTH
      USE DOMAIN_DECOMP, ONLY : HALO_UPDATE, CHECKSUM
      USE ICEDYN, only : nx1,ny1,form,relax,uice,vice,uicec,vicec
      IMPLICIT NONE
      REAL*8, DIMENSION(NX1,grid%J_STRT_HALO:grid%J_STOP_HALO) :: 
     &        USAVE,VSAVE
      REAL*8 rms
      INTEGER kki,i,j
      INTEGER :: J_0,J_1,J_0S,J_1S

C****
C**** Extract useful local domain parameters from "grid"
C****
      CALL GET(grid, J_STRT     =J_0,    J_STOP     =J_1,
     &               J_STRT_SKP =J_0S,   J_STOP_SKP =J_1S )

      rms=0.
C KKI LOOP IS FOR PSEUDO-TIMESTEPPING
      KKI=0.
 10   KKI=KKI+1

C FIRST DO PREDICTOR
      DO J=J_0,J_1
      DO I=1,NX1
       UICE(I,J,3)=UICE(I,J,1)
       VICE(I,J,3)=VICE(I,J,1)
       UICEC(I,J)=UICE(I,J,1)
       VICEC(I,J)=VICE(I,J,1)
      END DO
      END DO

      CALL FORM
      CALL RELAX

       DO J=J_0,J_1
       UICE(1,J,1)=UICE(NX1-1,J,1)
       VICE(1,J,1)=VICE(NX1-1,J,1)
       UICE(NX1,J,1)=UICE(2,J,1)
       VICE(NX1,J,1)=VICE(2,J,1)
      END DO

C NOW DO REGULAR TIME STEP
C NOW DO MODIFIED EULER STEP
c
      DO J=J_0,J_1
      DO I=1,NX1
       UICE(I,J,1)=0.5*(UICE(I,J,1)+UICE(I,J,2))
       VICE(I,J,1)=0.5*(VICE(I,J,1)+VICE(I,J,2))
      END DO
      END DO

      CALL FORM

C NOW SET U(1)=U(2) AND SAME FOR V
      DO J=J_0,J_1
      DO I=1,NX1
       UICE(I,J,3)=UICE(I,J,1)
       VICE(I,J,3)=VICE(I,J,1)
       UICEC(I,J)=UICE(I,J,1)
       VICEC(I,J)=VICE(I,J,1)
       UICE(I,J,1)=UICE(I,J,2)
       VICE(I,J,1)=VICE(I,J,2)
      END DO
      END DO

      CALL RELAX

       DO J=J_0,J_1
       UICE(1,J,1)=UICE(NX1-1,J,1)
       VICE(1,J,1)=VICE(NX1-1,J,1)
       UICE(NX1,J,1)=UICE(2,J,1)
       VICE(NX1,J,1)=VICE(2,J,1)
      END DO

      if (kki.gt.1) then ! test convergence
        rms=0.
        do i=1,nx1
           do j=j_0,j_1
            rms=rms+(USAVE(i,j)-UICE(i,j,1))**2+(VSAVE(i,j)-VICE(i,j,1))
     *           **2
          end do
        end do
      end if

      if (kki.eq.20) then
        write(6,*) "Too many iterations in VPICEDYN. kki:",kki,rms
      elseif (kki.eq.1 .or. rms.gt.0.01d0) then
        USAVE=UICE(:,:,1)
        VSAVE=VICE(:,:,1)
        goto 10 
      end if

      RETURN
      END SUBROUTINE VPICEDYN


      SUBROUTINE ALLOC_ICEDYN(grid) 1,10
!@sum ALLOC_ICEDYN allocates arrays defined in the ICEDYN module.
!@auth Rosalinda de Fainchtein

C**** arrays allocated in this routine were originally dimensioned
C**** (..,ny1,..). Since ny1=jm (see above), the grid structure as defined
C**** in DOMAIN_DECOMP can be used in the calling routine.
C**** In the case that ny1 is NOT equal to JM, a structure appropriately
C**** modified to reflect the differences should be created in DOMAIN_DECOMP 
C**** and used in the calling routine. No modification should be necesary
C**** to ALLOC_ICEDYN.

      USE DOMAIN_DECOMP, ONLY : DYN_GRID
      USE DOMAIN_DECOMP, ONLY : GET
      USE DOMAIN_DECOMP, ONLY : INIT_DECOMP
      USE ICEDYN, ONLY : FOCEAN
      USE ICEDYN, ONLY : PRESS,HEFFM,UVM,DWATN,COR,ZMAX,ZMIN,ETA,
     &                   ZETA,DRAGS,DRAGA,GAIRX,GAIRY,GWATX,GWATY,
     &                   PGFUB,PGFVB,FORCEX,FORCEY,AMASS,UICEC,
     &                   VICEC,UIB,VIB,DMU,DMV,HEFF,AREA,UICE,
     &                   VICE,SINEN,BYDXDY,DYT,DYU,BYDY2,BYDYR,
     &                   CST,CSU,TNGT,TNG,BYCSU
      USE ICEDYN, ONLY : IMIC, JMIC, NX1, NY1
      USE ICEDYN, ONLY : grid_MIC, grid_NXY
      IMPLICIT NONE
      LOGICAL, SAVE :: init = .false.
      TYPE (DYN_GRID), INTENT(IN) :: grid

      INTEGER :: I_0H, I_1H, J_1H, J_0H
      INTEGER :: IER

      INTEGER :: I,J,L

      If (init) Then
         Return ! Only invoke once
      End If
      init = .true.

      CALL INIT_DECOMP(grid_MIC, IMIC, JMIC)
      CALL INIT_DECOMP(grid_NXY, NX1, NY1)
      CALL GET( grid_NXY, I_STRT_HALO=I_0H, I_STOP_HALO=I_1H, 
     &                J_STRT_HALO=J_0H, J_STOP_HALO=J_1H  )

      ALLOCATE( FOCEAN(NX1-2,J_0H:J_1H),
     $   STAT = IER)

      ALLOCATE(  PRESS(NX1,J_0H:J_1H),
     &           HEFFM(NX1,J_0H:J_1H),
     &           UVM(NX1,J_0H:J_1H),
     &           DWATN(NX1,J_0H:J_1H),
     &           COR(NX1,J_0H:J_1H),
     *           ZMAX(NX1,J_0H:J_1H),
     &           ZMIN(NX1,J_0H:J_1H),
     &           ETA(NX1,J_0H:J_1H),
     &           ZETA(NX1,J_0H:J_1H),
     &           DRAGS(NX1,J_0H:J_1H),
     &           DRAGA(NX1,J_0H:J_1H),
     &           GAIRX(NX1,J_0H:J_1H),
     &           GAIRY(NX1,J_0H:J_1H),
     *           GWATX(NX1,J_0H:J_1H),
     &           GWATY(NX1,J_0H:J_1H),
     &           PGFUB(NX1,J_0H:J_1H),
     &           PGFVB(NX1,J_0H:J_1H),
     &           FORCEX(NX1,J_0H:J_1H),
     &           FORCEY(NX1,J_0H:J_1H),
     &           AMASS(NX1,J_0H:J_1H),
     &           UICEC(NX1,J_0H:J_1H),
     &           VICEC(NX1,J_0H:J_1H),
     &           UIB(NX1,J_0H:J_1H),
     *           VIB(NX1,J_0H:J_1H),
     &           DMU(NX1,J_0H:J_1H),
     &           DMV(NX1,J_0H:J_1H),
     &           HEFF(NX1,J_0H:J_1H),
     &           AREA(NX1,J_0H:J_1H),
     $   STAT = IER)
      ALLOCATE( UICE(NX1,J_0H:J_1H,3),
     &          VICE(NX1,J_0H:J_1H,3),
     $   STAT = IER)
C**** Geometry
      ALLOCATE( SINEN(NX1,J_0H:J_1H),
     &          BYDXDY(NX1,J_0H:J_1H),
     $   STAT = IER)

      ALLOCATE ( DYT(J_0H:J_1H),
     &           DYU(J_0H:J_1H),
     &           BYDY2(J_0H:J_1H),
     &           BYDYR(J_0H:J_1H),
     &           CST(J_0H:J_1H),
     &           CSU(J_0H:J_1H),
     &           TNGT(J_0H:J_1H),
     &           TNG(J_0H:J_1H),
     &           BYCSU(J_0H:J_1H),
     $   STAT = IER)

      RETURN
      END SUBROUTINE ALLOC_ICEDYN