C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/seaice/Attic/adi.F,v 1.12 2004/05/03 06:09:39 dimitri dead $ C $Name: $ #include "SEAICE_OPTIONS.h" CStartOfInterface SUBROUTINE adi( myThid ) C /==========================================================\ C | SUBROUTINE adi | C | o Solve ice momentum equation with an ADI dynamics solver| C | (see Zhang and Rothrock, JGR, 105, 3325-3338, 2000) | C | written by Jinlun Zhang, PSC/UW, Feb-2001 | C | zhang@apl.washington.edu | C |==========================================================| C \==========================================================/ IMPLICIT NONE C === Global variables === #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "SEAICE.h" #include "SEAICE_PARAMS.h" #include "SEAICE_GRID.h" C === Routine arguments === C myThid - Thread no. that called this routine. INTEGER myThid CEndOfInterface #ifdef ALLOW_SEAICE #ifdef SEAICE_ALLOW_DYNAMICS C === Local variables === C i,j,bi,bj - Loop counters INTEGER i, j, bi, bj, j1, j2, im, jm _RL AA1, AA2, AA3, AA4, AA5, AA6, AA9, RADIUS, RADIUS2 _RL AU (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL BU (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL CU (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL AV (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL BV (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL CV (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL DELXY (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL DELXR (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL DELYR (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL DELX2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL DELY2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL ETAMEAN (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL ZETAMEAN(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL FXY (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL FXY1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL URT(1-OLx:sNx+OLx), VRT(1-OLy:sNy+OLy) _RL CUU(1-OLx:sNx+OLx), CVV(1-OLy:sNy+OLy) C SET SOME VALUES RADIUS=6370. _d 3 RADIUS2=ONE/(RADIUS*RADIUS) C SOLVE FOR UICE C FIRST HALF FIRST c DO bj=myByLo(myThid),myByHi(myThid) c DO bi=myBxLo(myThid),myBxHi(myThid) c DO j=1,sNy c DO i=1,sNx c FORCEX(I,J,bi,bj)=FORCEX(I,J,bi,bj)*UVM(I,J,bi,bj) c FORCEY(I,J,bi,bj)=FORCEY(I,J,bi,bj)*UVM(I,J,bi,bj) c ENDDO c ENDDO c ENDDO c ENDDO C-- Update overlap regions CALL EXCH_UV_XY_RL(FORCEX,FORCEY,.TRUE.,myThid) _EXCH_XY_R8(DRAGS, myThid) _EXCH_XY_R8(DRAGA, myThid) _EXCH_XY_R8(AMASS, myThid) c$taf loop = parallel DO bj=myByLo(myThid),myByHi(myThid) c$taf loop = parallel DO bi=myBxLo(myThid),myBxHi(myThid) DO J=1-OLy+1,sNy+OLy-1 DO I=1-OLx+1,sNx+OLx-1 DELXY(I,J)=HALF/(DXUICE(I,J,bi,bj)*DYUICE(I,J,bi,bj)) DELXR(I,J)=HALF/(DXUICE(I,J,bi,bj)*RADIUS) DELX2(I,J)=HALF/(DXUICE(I,J,bi,bj)*DXUICE(I,J,bi,bj)) ETAMEAN(I,J)=QUART*(ETA(I,J-1,bi,bj)+ETA(I-1,J-1,bi,bj) & +ETA(I,J,bi,bj)+ETA(I-1,J,bi,bj)) ZETAMEAN(I,J)=QUART*(ZETA(I,J-1,bi,bj)+ZETA(I-1,J-1,bi,bj) & +ZETA(I,J,bi,bj)+ZETA(I-1,J,bi,bj)) FXY(I,J)=DRAGA(I,J,bi,bj)*VICEC(I,J,bi,bj)+FORCEX(I,J,bi,bj) 3 +HALF*(ZETA(I,J,bi,bj)*(VICEC(I+1,J+1,bi,bj) 3 +VICEC(I,J+1,bi,bj) 3 -VICEC(I+1,J,bi,bj)-VICEC(I,J,bi,bj)) 3 +ZETA(I,J-1,bi,bj)*(VICEC(I+1,J,bi,bj) 3 +VICEC(I,J,bi,bj)-VICEC(I+1,J-1,bi,bj) 3 -VICEC(I,J-1,bi,bj))+ZETA(I-1,J,bi,bj) 3 *(VICEC(I,J,bi,bj)+VICEC(I-1,J,bi,bj) 3 -VICEC(I,J+1,bi,bj)-VICEC(I-1,J+1,bi,bj)) 3 +ZETA(I-1,J-1,bi,bj)*(VICEC(I,J-1,bi,bj) 3 +VICEC(I-1,J-1,bi,bj)-VICEC(I,J,bi,bj) 3 -VICEC(I-1,J,bi,bj)))*DELXY(I,J) 3 *RECIP_CSUICE(I,J,bi,bj) 3 4 -HALF*(ETA(I,J,bi,bj)*(VICEC(I+1,J+1,bi,bj) 4 +VICEC(I,J+1,bi,bj) 4 -VICEC(I+1,J,bi,bj)-VICEC(I,J,bi,bj)) 4 +ETA(I,J-1,bi,bj)*(VICEC(I+1,J,bi,bj) 4 +VICEC(I,J,bi,bj)-VICEC(I+1,J-1,bi,bj) 4 -VICEC(I,J-1,bi,bj))+ETA(I-1,J,bi,bj) 4 *(VICEC(I,J,bi,bj)+VICEC(I-1,J,bi,bj) 4 -VICEC(I,J+1,bi,bj)-VICEC(I-1,J+1,bi,bj)) 4 +ETA(I-1,J-1,bi,bj)*(VICEC(I,J-1,bi,bj) 4 +VICEC(I-1,J-1,bi,bj)-VICEC(I,J,bi,bj) 4 -VICEC(I-1,J,bi,bj)))*DELXY(I,J) 4 *RECIP_CSUICE(I,J,bi,bj) 4 5 +HALF*(VICEC(I+1,J,bi,bj)-VICEC(I-1,J,bi,bj)) 5 *(ETA(I-1,J,bi,bj)+ETA(I,J,bi,bj) 5 -ETA(I-1,J-1,bi,bj)-ETA(I,J-1,bi,bj))*DELXY(I,J) 5 *RECIP_CSUICE(I,J,bi,bj)+HALF*ETAMEAN(I,J) 5 *((VICEC(I+1,J+1,bi,bj) 5 -VICEC(I-1,J+1,bi,bj))*RECIP_CSUICE(I,J+1,bi,bj) 5 -(VICEC(I+1,J-1,bi,bj)-VICEC(I-1,J-1,bi,bj)) 5 *RECIP_CSUICE(I,J-1,bi,bj))*DELXY(I,J) 5 6 -((ZETA(I,J,bi,bj)+ZETA(I,J-1,bi,bj) 6 -ZETA(I-1,J-1,bi,bj)-ZETA(I-1,J,bi,bj)) 6 +(ETA(I,J,bi,bj)+ETA(I,J-1,bi,bj) 6 -ETA(I-1,J-1,bi,bj)-ETA(I-1,J,bi,bj))) 6 *TNGICE(I,J,bi,bj)*VICEC(I,J,bi,bj) & *DELXR(I,J)*RECIP_CSUICE(I,J,bi,bj) 6 -(ETAMEAN(I,J)+ZETAMEAN(I,J))*TNGICE(I,J,bi,bj) 6 *(VICEC(I+1,J,bi,bj)-VICEC(I-1,J,bi,bj)) 6 *DELXR(I,J)*RECIP_CSUICE(I,J,bi,bj) 6 7 -ETAMEAN(I,J)*TWO*TNGICE(I,J,bi,bj) 7 *(VICEC(I+1,J,bi,bj)-VICEC(I-1,J,bi,bj)) 7 *DELXR(I,J)*RECIP_CSUICE(I,J,bi,bj) END DO END DO DO J=1-OLy+1,sNy+OLy-1 DO I=1-OLx+1,sNx+OLx-1 DELY2(I,J)=HALF/(DYUICE(I,J,bi,bj)*DYUICE(I,J,bi,bj)) DELYR(I,J)=HALF/(DYUICE(I,J,bi,bj)*RADIUS) AA1=((ETA(I,J-1,bi,bj)+ZETA(I,J-1,bi,bj)) & *RECIP_CSUICE(I,J,bi,bj) & +(ETA(I,J,bi,bj)+ZETA(I,J,bi,bj)) & *RECIP_CSUICE(I,J,bi,bj))*RECIP_CSUICE(I,J,bi,bj) AA2=((ETA(I-1,J-1,bi,bj)+ZETA(I-1,J-1,bi,bj)) & *RECIP_CSUICE(I,J,bi,bj) & +(ETA(I-1,J,bi,bj)+ZETA(I-1,J,bi,bj)) & *RECIP_CSUICE(I,J,bi,bj))*RECIP_CSUICE(I,J,bi,bj) AA3=ETA(I-1,J,bi,bj)+ETA(I,J,bi,bj) AA4=ETA(I-1,J-1,bi,bj)+ETA(I,J-1,bi,bj) AA5=-(ETA(I-1,J,bi,bj)+ETA(I,J,bi,bj)-ETA(I-1,J-1,bi,bj) & -ETA(I,J-1,bi,bj))*TNGICE(I,J,bi,bj) AA6=TWO*ETAMEAN(I,J)*TNGICE(I,J,bi,bj)*TNGICE(I,J,bi,bj) AU(I,J)=-AA2*DELX2(I,J)*UVM(I,J,bi,bj) BU(I,J)=((AA1+AA2)*DELX2(I,J)+AA6*RADIUS2 & +AMASS(I,J,bi,bj)/SEAICE_DT*TWO & +DRAGS(I,J,bi,bj))*UVM(I,J,bi,bj) & +(ONE-UVM(I,J,bi,bj)) CU(I,J)=-AA1*DELX2(I,J)*UVM(I,J,bi,bj) END DO END DO DO J=1-OLy+1,sNy+OLy-1 AU(1-OLx+1,J)=ZERO CU(sNx+OLx-1,J)=ZERO CU(1-OLx+1,J)=CU(1-OLx+1,J)/BU(1-OLx+1,J) END DO c$taf loop = parallel DO 1200 J=1-OLy+1,sNy+OLy-1 DO I=1-OLx+1,sNx+OLx-1 AA1=((ETA(I,J-1,bi,bj)+ZETA(I,J-1,bi,bj)) & *RECIP_CSUICE(I,J,bi,bj) & +(ETA(I,J,bi,bj)+ZETA(I,J,bi,bj)) & *RECIP_CSUICE(I,J,bi,bj))*RECIP_CSUICE(I,J,bi,bj) AA2=((ETA(I-1,J-1,bi,bj)+ZETA(I-1,J-1,bi,bj)) & *RECIP_CSUICE(I,J,bi,bj) & +(ETA(I-1,J,bi,bj)+ZETA(I-1,J,bi,bj)) & *RECIP_CSUICE(I,J,bi,bj))*RECIP_CSUICE(I,J,bi,bj) AA3=ETA(I-1,J,bi,bj)+ETA(I,J,bi,bj) AA4=ETA(I-1,J-1,bi,bj)+ETA(I,J-1,bi,bj) AA5=-(ETA(I-1,J,bi,bj)+ETA(I,J,bi,bj)-ETA(I-1,J-1,bi,bj) & -ETA(I,J-1,bi,bj))*TNGICE(I,J,bi,bj) AA6=TWO*ETAMEAN(I,J)*TNGICE(I,J,bi,bj)*TNGICE(I,J,bi,bj) IF(I.EQ.1-OLx+1) THEN AA9=AA2*DELX2(I,J)*UICEC(I-1,J,bi,bj)*UVM(I,J,bi,bj) ELSE IF(I.EQ.sNx+OLx-1) THEN AA9=AA1*DELX2(I,J)*UICEC(I+1,J,bi,bj)*UVM(I,J,bi,bj) ELSE AA9=ZERO END IF URT(I)=AA9+FXY(I,J)-AA5*DELYR(I,J)*UICE(I,J,2,bi,bj) 1 -(AA3+AA4)*DELY2(I,J)*UICE(I,J,2,bi,bj) 1 +(ETA(I-1,J,bi,bj)+ETA(I,J,bi,bj)) 1 *UICE(I,J+1,2,bi,bj)*DELY2(I,J) 2 +(ETA(I-1,J-1,bi,bj)+ETA(I,J-1,bi,bj)) 2 *UICE(I,J-1,2,bi,bj)*DELY2(I,J) 3 +ETAMEAN(I,J)*DELYR(I,J)*(UICE(I,J+1,2,bi,bj) 3 *TNGICE(I,J+1,bi,bj) 3 -UICE(I,J-1,2,bi,bj)*TNGICE(I,J-1,bi,bj)) 4 -ETAMEAN(I,J)*DELYR(I,J)*TWO*TNGICE(I,J,bi,bj) 4 *(UICE(I,J+1,2,bi,bj) 4 -UICE(I,J-1,2,bi,bj)) URT(I)=(URT(I)+AMASS(I,J,bi,bj)/SEAICE_DT & *UICE(I,J,2,bi,bj)*TWO)*UVM(I,J,bi,bj) END DO DO I=1-OLx+1,sNx+OLx-1 CUU(I)=CU(I,J) END DO URT(1-OLx+1)=URT(1-OLx+1)/BU(1-OLx+1,J) DO I=1-OLx+2,sNx+OLx-1 IM=I-1 CUU(I)=CUU(I)/(BU(I,J)-AU(I,J)*CUU(IM)) URT(I)=(URT(I)-AU(I,J)*URT(IM))/(BU(I,J)-AU(I,J)*CUU(IM)) END DO DO I=1-OLx+1,sNx+OLx-1-1 J1=sNx+OLx-1-I J2=J1+1 URT(J1)=URT(J1)-CUU(J1)*URT(J2) END DO DO I=1-OLx+1,sNx+OLx-1 UICE(I,J,1,bi,bj)=URT(I) END DO 1200 CONTINUE c DO J=1,sNy c DO I=1,sNx c UICE(I,J,3,bi,bj)=UICE(I,J,1,bi,bj) c END DO c END DO ENDDO ENDDO CALL SEAICE_EXCH( UICE, myThid ) C NOW SECOND HALF c$taf loop = parallel DO bj=myByLo(myThid),myByHi(myThid) c$taf loop = parallel DO bi=myBxLo(myThid),myBxHi(myThid) DO I=1-OLx+1,sNx+OLx-1 DO J=1-OLy+1,sNy+OLy-1 AA1=ETA(I-1,J,bi,bj)+ETA(I,J,bi,bj) AA2=ETA(I-1,J-1,bi,bj)+ETA(I,J-1,bi,bj) AA5=-(ETA(I-1,J,bi,bj)+ETA(I,J,bi,bj) & -ETA(I-1,J-1,bi,bj)-ETA(I,J-1,bi,bj))*TNGICE(I,J,bi,bj) AA6=TWO*ETAMEAN(I,J)*TNGICE(I,J,bi,bj)*TNGICE(I,J,bi,bj) AV(I,J)=(-AA2*DELY2(I,J)+ETAMEAN(I,J)*DELYR(I,J) & *(TNGICE(I,J-1,bi,bj) & -TWO*TNGICE(I,J,bi,bj)))*UVM(I,J,bi,bj) BV(I,J)=((AA1+AA2)*DELY2(I,J)+AA5*DELYR(I,J)+AA6*RADIUS2 & +AMASS(I,J,bi,bj)/SEAICE_DT*TWO & +DRAGS(I,J,bi,bj))*UVM(I,J,bi,bj) & +(ONE-UVM(I,J,bi,bj)) CV(I,J)=(-AA1*DELY2(I,J)-ETAMEAN(I,J)*DELYR(I,J) & *(TNGICE(I,J+1,bi,bj) & -TWO*TNGICE(I,J,bi,bj)))*UVM(I,J,bi,bj) END DO END DO DO I=1-OLx+1,sNx+OLx-1 AV(I,1-OLy+1)=ZERO CV(I,sNy+OLy-1)=ZERO CV(I,1-OLy+1)=CV(I,1-OLy+1)/BV(I,1-OLy+1) END DO DO I=1-OLx+1,sNx+OLx-1 DO J=1-OLy+1,sNy+OLy-1 AA1=((ETA(I,J-1,bi,bj)+ZETA(I,J-1,bi,bj)) & *RECIP_CSUICE(I,J,bi,bj) & +(ETA(I,J,bi,bj)+ZETA(I,J,bi,bj)) & *RECIP_CSUICE(I,J,bi,bj))*RECIP_CSUICE(I,J,bi,bj) AA2=((ETA(I-1,J-1,bi,bj)+ZETA(I-1,J-1,bi,bj)) & *RECIP_CSUICE(I,J,bi,bj) & +(ETA(I-1,J,bi,bj)+ZETA(I-1,J,bi,bj)) & *RECIP_CSUICE(I,J,bi,bj))*RECIP_CSUICE(I,J,bi,bj) AA3=ETA(I-1,J,bi,bj)+ETA(I,J,bi,bj) AA4=ETA(I-1,J-1,bi,bj)+ETA(I,J-1,bi,bj) IF(J.EQ.1-OLy+1) THEN AA9=( AA4*DELY2(I,J)*UICEC(I,J-1,bi,bj) & -ETAMEAN(I,J)*DELYR(I,J)*(TNGICE(I,J-1,bi,bj) & -TWO*TNGICE(I,J,bi,bj)) & *UICEC(I,J-1,bi,bj) )*UVM(I,J,bi,bj) ELSE IF(J.EQ.sNy+OLy-1) THEN AA9=( AA3*DELY2(I,J)*UICEC(I,J+1,bi,bj) & +ETAMEAN(I,J)*DELYR(I,J)*(TNGICE(I,J+1,bi,bj) & -TWO*TNGICE(I,J,bi,bj)) & *UICEC(I,J+1,bi,bj) )*UVM(I,J,bi,bj) ELSE AA9=ZERO END IF FXY1(I,J)=AA9+AMASS(I,J,bi,bj)/SEAICE_DT*UICE(I,J,1,bi,bj)*TWO 5 -(AA1+AA2)*DELX2(I,J)*UICE(I,J,1,bi,bj) 6 +((ETA(I,J-1,bi,bj)+ZETA(I,J-1,bi,bj) 6 +ETA(I,J,bi,bj)+ZETA(I,J,bi,bj)) 6 *UICE(I+1,J,1,bi,bj) 6 +(ETA(I-1,J-1,bi,bj)+ZETA(I-1,J-1,bi,bj) 6 + ETA(I-1,J ,bi,bj)+ZETA(I-1,J ,bi,bj)) 6 *UICE(I-1,J,1,bi,bj))*DELX2(I,J) 6 *RECIP_CSUICE(I,J,bi,bj)*RECIP_CSUICE(I,J,bi,bj) END DO END DO DO 1300 I=1-OLx+1,sNx+OLx-1 DO J=1-OLy+1,sNy+OLy-1 VRT(J)=FXY(I,J)+FXY1(I,J) VRT(J)=VRT(J)*UVM(I,J,bi,bj) END DO DO J=1-OLy+1,sNy+OLy-1 CVV(J)=CV(I,J) END DO VRT(1-OLy+1)=VRT(1-OLy+1)/BV(I,1-OLy+1) DO J=1-OLy+2,sNy+OLy-1 JM=J-1 CVV(J)=CVV(J)/(BV(I,J)-AV(I,J)*CVV(JM)) VRT(J)=(VRT(J)-AV(I,J)*VRT(JM))/(BV(I,J)-AV(I,J)*CVV(JM)) END DO DO J=1-OLy+1,sNy+OLy-1-1 J1=sNy+OLy-1-J J2=J1+1 VRT(J1)=VRT(J1)-CVV(J1)*VRT(J2) END DO DO J=1-OLy+1,sNy+OLy-1 UICE(I,J,1,bi,bj)=VRT(J) END DO 1300 CONTINUE ENDDO ENDDO C SOLVE FOR VICE C FIRST HALF FIRST c$taf loop = parallel DO bj=myByLo(myThid),myByHi(myThid) c$taf loop = parallel DO bi=myBxLo(myThid),myBxHi(myThid) DO I=1-OLx+1,sNx+OLx-1 DO J=1-OLy+1,sNy+OLy-1 FXY(I,J)=-DRAGA(I,J,bi,bj)*UICEC(I,J,bi,bj)+FORCEY(I,J,bi,bj) 3 +(HALF*(UICEC(I+1,J,bi,bj)-UICEC(I-1,J,bi,bj)) 3 *(ZETA(I-1,J,bi,bj)+ZETA(I,J,bi,bj) 3 -ZETA(I-1,J-1,bi,bj)-ZETA(I,J-1,bi,bj))*DELXY(I,J) 3 *RECIP_CSUICE(I,J,bi,bj)+HALF*ZETAMEAN(I,J) 3 *((UICEC(I+1,J+1,bi,bj) 3 -UICEC(I-1,J+1,bi,bj))*RECIP_CSUICE(I,J+1,bi,bj) 3 -(UICEC(I+1,J-1,bi,bj)-UICEC(I-1,J-1,bi,bj)) 3 *RECIP_CSUICE(I,J-1,bi,bj))*DELXY(I,J)) 3 4 -(HALF*(UICEC(I+1,J,bi,bj)-UICEC(I-1,J,bi,bj)) 4 *(ETA(I-1,J,bi,bj)+ETA(I,J,bi,bj) 4 -ETA(I-1,J-1,bi,bj)-ETA(I,J-1,bi,bj))*DELXY(I,J) 4 *RECIP_CSUICE(I,J,bi,bj)+HALF*ETAMEAN(I,J) 4 *((UICEC(I+1,J+1,bi,bj) 4 -UICEC(I-1,J+1,bi,bj))*RECIP_CSUICE(I,J+1,bi,bj) 4 -(UICEC(I+1,J-1,bi,bj)-UICEC(I-1,J-1,bi,bj)) 4 *RECIP_CSUICE(I,J-1,bi,bj))*DELXY(I,J)) 4 5 +HALF*(ETA(I,J,bi,bj)*(UICEC(I+1,J+1,bi,bj) 5 +UICEC(I,J+1,bi,bj) 5 -UICEC(I+1,J,bi,bj)-UICEC(I,J,bi,bj))+ETA(I,J-1,bi,bj) 5 *(UICEC(I+1,J,bi,bj) 5 +UICEC(I,J,bi,bj)-UICEC(I+1,J-1,bi,bj)-UICEC(I,J-1,bi,bj)) 5 +ETA(I-1,J,bi,bj) 5 *(UICEC(I,J,bi,bj)+UICEC(I-1,J,bi,bj)-UICEC(I,J+1,bi,bj) 5 -UICEC(I-1,J+1,bi,bj)) 5 +ETA(I-1,J-1,bi,bj)*(UICEC(I,J-1,bi,bj) 5 +UICEC(I-1,J-1,bi,bj) 5 -UICEC(I,J,bi,bj) 5 -UICEC(I-1,J,bi,bj)))*DELXY(I,J)*RECIP_CSUICE(I,J,bi,bj) 5 6 +(ETA(I,J,bi,bj)+ETA(I,J-1,bi,bj) 6 -ETA(I-1,J-1,bi,bj)-ETA(I-1,J,bi,bj)) 6 *TNGICE(I,J,bi,bj)*UICEC(I,J,bi,bj) 6 *DELXR(I,J)*RECIP_CSUICE(I,J,bi,bj) 6 +ETAMEAN(I,J)*TNGICE(I,J,bi,bj)*(UICEC(I+1,J,bi,bj) 6 -UICEC(I-1,J,bi,bj))*DELXR(I,J)*RECIP_CSUICE(I,J,bi,bj) 6 7 +ETAMEAN(I,J)*TWO*TNGICE(I,J,bi,bj)*(UICEC(I+1,J,bi,bj) 7 -UICEC(I-1,J,bi,bj))*DELXR(I,J)*RECIP_CSUICE(I,J,bi,bj) END DO END DO DO I=1-OLx+1,sNx+OLx-1 DO J=1-OLy+1,sNy+OLy-1 AA1=ETA(I-1,J,bi,bj)+ZETA(I-1,J,bi,bj) & +ETA(I,J,bi,bj)+ZETA(I,J,bi,bj) AA2=ETA(I-1,J-1,bi,bj)+ZETA(I-1,J-1,bi,bj) & +ETA(I,J-1,bi,bj)+ZETA(I,J-1,bi,bj) AA3=(ETA(I,J-1,bi,bj)*RECIP_CSUICE(I,J,bi,bj)+ETA(I,J,bi,bj) & *RECIP_CSUICE(I,J,bi,bj))*RECIP_CSUICE(I,J,bi,bj) AA4=(ETA(I-1,J-1,bi,bj)+ETA(I-1,J,bi,bj)) & *RECIP_CSUICE(I,J,bi,bj)*RECIP_CSUICE(I,J,bi,bj) AA5=((ZETA(I-1,J,bi,bj)-ETA(I-1,J,bi,bj)) & +(ZETA(I,J,bi,bj)-ETA(I,J,bi,bj)) & -(ZETA(I-1,J-1,bi,bj)-ETA(I-1,J-1,bi,bj)) & -(ZETA(I,J-1,bi,bj) & -ETA(I,J-1,bi,bj)))*TNGICE(I,J,bi,bj) AA6=TWO*ETAMEAN(I,J)*TNGICE(I,J,bi,bj)*TNGICE(I,J,bi,bj) AV(I,J)=(-AA2*DELY2(I,J)-(ZETAMEAN(I,J)-ETAMEAN(I,J)) & *TNGICE(I,J-1,bi,bj) & *DELYR(I,J)-ETAMEAN(I,J)*TWO*TNGICE(I,J,bi,bj) & *DELYR(I,J))*UVM(I,J,bi,bj) BV(I,J)=((AA1+AA2)*DELY2(I,J)+AA5*DELYR(I,J)+AA6*RADIUS2 & +AMASS(I,J,bi,bj)/SEAICE_DT*TWO & +DRAGS(I,J,bi,bj))*UVM(I,J,bi,bj) & +(ONE-UVM(I,J,bi,bj)) CV(I,J)=(-AA1*DELY2(I,J)+(ZETAMEAN(I,J)-ETAMEAN(I,J)) & *TNGICE(I,J+1,bi,bj) & *DELYR(I,J)+ETAMEAN(I,J)*TWO*TNGICE(I,J,bi,bj) & *DELYR(I,J))*UVM(I,J,bi,bj) END DO END DO DO I=1-OLx+1,sNx+OLx-1 AV(I,1-OLy+1)=ZERO CV(I,sNy+OLy-1)=ZERO CV(I,1-OLy+1)=CV(I,1-OLy+1)/BV(I,1-OLy+1) END DO c$taf loop = parallel DO 1301 I=1-OLx+1,sNx+OLx-1 DO J=1-OLy+1,sNy+OLy-1 AA1=ETA(I-1,J,bi,bj)+ZETA(I-1,J,bi,bj) & +ETA(I,J,bi,bj)+ZETA(I,J,bi,bj) AA2=ETA(I-1,J-1,bi,bj)+ZETA(I-1,J-1,bi,bj) & +ETA(I,J-1,bi,bj)+ZETA(I,J-1,bi,bj) AA3=(ETA(I,J-1,bi,bj)*RECIP_CSUICE(I,J,bi,bj)+ETA(I,J,bi,bj) & *RECIP_CSUICE(I,J,bi,bj))*RECIP_CSUICE(I,J,bi,bj) AA4=(ETA(I-1,J-1,bi,bj)+ETA(I-1,J,bi,bj)) & *RECIP_CSUICE(I,J,bi,bj)*RECIP_CSUICE(I,J,bi,bj) AA5=((ZETA(I-1,J,bi,bj)-ETA(I-1,J,bi,bj)) & +(ZETA(I,J,bi,bj)-ETA(I,J,bi,bj)) & -(ZETA(I-1,J-1,bi,bj)-ETA(I-1,J-1,bi,bj)) & -(ZETA(I,J-1,bi,bj)-ETA(I,J-1,bi,bj)))*TNGICE(I,J,bi,bj) AA6=TWO*ETAMEAN(I,J)*TNGICE(I,J,bi,bj)*TNGICE(I,J,bi,bj) IF(J.EQ.1-OLy+1) THEN AA9=(AA2*DELY2(I,J)+(ZETAMEAN(I,J)-ETAMEAN(I,J)) & *TNGICE(I,J-1,bi,bj)*DELYR(I,J) & +ETAMEAN(I,J)*TWO*TNGICE(I,J,bi,bj)*DELYR(I,J)) & *VICEC(I,J-1,bi,bj)*UVM(I,J,bi,bj) ELSE IF(J.EQ.sNy+OLy-1) THEN AA9=(AA1*DELY2(I,J)-(ZETAMEAN(I,J)-ETAMEAN(I,J)) & *TNGICE(I,J+1,bi,bj)*DELYR(I,J) & -ETAMEAN(I,J)*TWO*TNGICE(I,J,bi,bj)*DELYR(I,J)) & *VICEC(I,J+1,bi,bj)*UVM(I,J,bi,bj) ELSE AA9=ZERO END IF VRT(J)=AA9+FXY(I,J)-(AA3+AA4)*DELX2(I,J)*VICE(I,J,2,bi,bj) 6 +((ETA(I,J-1,bi,bj)*RECIP_CSUICE(I,J,bi,bj)+ETA(I,J,bi,bj) 6 *RECIP_CSUICE(I,J,bi,bj))*VICE(I+1,J,2,bi,bj)*DELX2(I,J) 7 +(ETA(I-1,J-1,bi,bj)+ETA(I-1,J,bi,bj)) 7 *VICE(I-1,J,2,bi,bj)*DELX2(I,J)) 7 *RECIP_CSUICE(I,J,bi,bj)*RECIP_CSUICE(I,J,bi,bj) VRT(J)=(VRT(J)+AMASS(I,J,bi,bj)/SEAICE_DT & *VICE(I,J,2,bi,bj)*TWO) & *UVM(I,J,bi,bj) END DO DO J=1-OLy+1,sNy+OLy-1 CVV(J)=CV(I,J) END DO VRT(1-OLy+1)=VRT(1-OLy+1)/BV(I,1-OLy+1) DO J=1-OLy+2,sNy+OLy-1 JM=J-1 CVV(J)=CVV(J)/(BV(I,J)-AV(I,J)*CVV(JM)) VRT(J)=(VRT(J)-AV(I,J)*VRT(JM))/(BV(I,J)-AV(I,J)*CVV(JM)) END DO DO J=1-OLy+1,sNy+OLy-1-1 J1=sNy+OLy-1-J J2=J1+1 VRT(J1)=VRT(J1)-CVV(J1)*VRT(J2) END DO DO J=1-OLy+1,sNy+OLy-1 VICE(I,J,1,bi,bj)=VRT(J) END DO 1301 CONTINUE c DO J=1,sNy c DO I=1,sNx c VICE(I,J,3,bi,bj)=VICE(I,J,1,bi,bj) c END DO c END DO ENDDO ENDDO CALL SEAICE_EXCH( VICE, myThid ) C NOW SECOND HALF c$taf loop = parallel DO bj=myByLo(myThid),myByHi(myThid) c$taf loop = parallel DO bi=myBxLo(myThid),myBxHi(myThid) DO J=1-OLy+1,sNy+OLy-1 DO I=1-OLx+1,sNx+OLx-1 AA1=ETA(I-1,J,bi,bj)+ZETA(I-1,J,bi,bj) & +ETA(I,J,bi,bj)+ZETA(I,J,bi,bj) AA2=ETA(I-1,J-1,bi,bj)+ZETA(I-1,J-1,bi,bj) & +ETA(I,J-1,bi,bj)+ZETA(I,J-1,bi,bj) AA3=(ETA(I,J-1,bi,bj)*RECIP_CSUICE(I,J,bi,bj)+ETA(I,J,bi,bj) & *RECIP_CSUICE(I,J,bi,bj))*RECIP_CSUICE(I,J,bi,bj) AA4=(ETA(I-1,J-1,bi,bj)+ETA(I-1,J,bi,bj)) & *RECIP_CSUICE(I,J,bi,bj)*RECIP_CSUICE(I,J,bi,bj) AA6=TWO*ETAMEAN(I,J)*TNGICE(I,J,bi,bj)*TNGICE(I,J,bi,bj) AU(I,J)=-AA4*DELX2(I,J)*UVM(I,J,bi,bj) BU(I,J)=((AA3+AA4)*DELX2(I,J)+AA6*RADIUS2 & +AMASS(I,J,bi,bj)/SEAICE_DT*TWO & +DRAGS(I,J,bi,bj))*UVM(I,J,bi,bj) & +(ONE-UVM(I,J,bi,bj)) CU(I,J)=-AA3*DELX2(I,J)*UVM(I,J,bi,bj) END DO END DO DO J=1-OLy+1,sNy+OLy-1 AU(1-OLx+1,J)=ZERO CU(sNx+OLx-1,J)=ZERO CU(1-OLx+1,J)=CU(1-OLx+1,J)/BU(1-OLx+1,J) END DO DO J=1-OLy+1,sNy+OLy-1 DO I=1-OLx+1,sNx+OLx-1 AA1=ETA(I-1,J,bi,bj)+ZETA(I-1,J,bi,bj) & +ETA(I,J,bi,bj)+ZETA(I,J,bi,bj) AA2=ETA(I-1,J-1,bi,bj)+ZETA(I-1,J-1,bi,bj) & +ETA(I,J-1,bi,bj)+ZETA(I,J-1,bi,bj) AA3=(ETA(I,J-1,bi,bj)*RECIP_CSUICE(I,J,bi,bj)+ETA(I,J,bi,bj) & *RECIP_CSUICE(I,J,bi,bj))*RECIP_CSUICE(I,J,bi,bj) AA4=(ETA(I-1,J-1,bi,bj)+ETA(I-1,J,bi,bj)) & *RECIP_CSUICE(I,J,bi,bj)*RECIP_CSUICE(I,J,bi,bj) AA5=((ZETA(I-1,J,bi,bj)-ETA(I-1,J,bi,bj)) & +(ZETA(I,J,bi,bj)-ETA(I,J,bi,bj)) & -(ZETA(I-1,J-1,bi,bj)-ETA(I-1,J-1,bi,bj)) & -(ZETA(I,J-1,bi,bj) & -ETA(I,J-1,bi,bj)))*TNGICE(I,J,bi,bj) AA6=TWO*ETAMEAN(I,J)*TNGICE(I,J,bi,bj)*TNGICE(I,J,bi,bj) IF(I.EQ.1-OLx+1) THEN AA9=AA4*DELX2(I,J)*VICEC(I-1,J,bi,bj)*UVM(I,J,bi,bj) ELSE IF(I.EQ.sNx+OLx-1) THEN AA9=AA3*DELX2(I,J)*VICEC(I+1,J,bi,bj)*UVM(I,J,bi,bj) ELSE AA9=ZERO END IF FXY1(I,J)=AA9+AMASS(I,J,bi,bj)/SEAICE_DT 1 *VICE(I,J,1,bi,bj)*TWO 1 -AA5*DELYR(I,J)*VICE(I,J,1,bi,bj) 1 -(AA1+AA2)*DELY2(I,J)*VICE(I,J,1,bi,bj) 1 +AA1*DELY2(I,J)*VICE(I,J+1,1,bi,bj) 1 -((ZETAMEAN(I,J)-ETAMEAN(I,J)) 1 *TNGICE(I,J+1,bi,bj)*DELYR(I,J) 1 +ETAMEAN(I,J)*TWO*TNGICE(I,J,bi,bj)*DELYR(I,J)) 1 *VICE(I,J+1,1,bi,bj) 2 +AA2*DELY2(I,J)*VICE(I,J-1,1,bi,bj) 2 +((ZETAMEAN(I,J)-ETAMEAN(I,J)) 2 *TNGICE(I,J-1,bi,bj)*DELYR(I,J) 2 +ETAMEAN(I,J)*TWO*TNGICE(I,J,bi,bj)*DELYR(I,J)) 2 *VICE(I,J-1,1,bi,bj) END DO END DO DO 1201 J=1-OLy+1,sNy+OLy-1 DO I=1-OLx+1,sNx+OLx-1 URT(I)=FXY(I,J)+FXY1(I,J) URT(I)=URT(I)*UVM(I,J,bi,bj) END DO DO I=1-OLx+1,sNx+OLx-1 CUU(I)=CU(I,J) END DO URT(1-OLx+1)=URT(1-OLx+1)/BU(1-OLx+1,J) DO I=1-OLx+2,sNx+OLx-1 IM=I-1 CUU(I)=CUU(I)/(BU(I,J)-AU(I,J)*CUU(IM)) URT(I)=(URT(I)-AU(I,J)*URT(IM))/(BU(I,J)-AU(I,J)*CUU(IM)) END DO DO I=1-OLx+1,sNx+OLx-1-1 J1=sNx+OLx-1-I J2=J1+1 URT(J1)=URT(J1)-CUU(J1)*URT(J2) END DO DO I=1-OLx+1,sNx+OLx-1 VICE(I,J,1,bi,bj)=URT(I) END DO 1201 CONTINUE ENDDO ENDDO DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO J=1,sNy DO I=1,sNx UICEC(I,J,bi,bj)=UICE(I,J,1,bi,bj)*UVM(I,J,bi,bj) VICEC(I,J,bi,bj)=VICE(I,J,1,bi,bj)*UVM(I,J,bi,bj) END DO END DO ENDDO ENDDO C-- Update overlap regions CALL EXCH_UV_XY_RL(UICEC,VICEC,.TRUE.,myThid) DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx UICE(I,J,1,bi,bj)=UICEC(I,J,bi,bj) VICE(I,J,1,bi,bj)=VICEC(I,J,bi,bj) ENDDO ENDDO ENDDO ENDDO #endif /* SEAICE_ALLOW_DYNAMICS */ #endif /* ALLOW_SEAICE */ RETURN END