SUBROUTINE SINT(XF, &
                   ims, ime, jms, jme, icmask , &
                   its, ite, jts, jte, nf, xstag, ystag )
      IMPLICIT NONE
      INTEGER ims, ime, jms, jme, &
              its, ite, jts, jte
      LOGICAL icmask( ims:ime, jms:jme )
      LOGICAL xstag, ystag
      INTEGER nf, ior
      REAL one12, one24, ep
      PARAMETER(one12=1./12.,one24=1./24.)
      PARAMETER(ior=2)
      REAL XF(ims:ime,jms:jme,NF)
      REAL Y(ims:ime,jms:jme,-IOR:IOR), &
           Z(ims:ime,jms:jme,-IOR:IOR), &
           F(ims:ime,jms:jme,0:1)
      INTEGER I,J,II,JJ,IIM
      INTEGER N2STAR, N2END, N1STAR, N1END
      DATA EP/ 1.E-10/
      REAL W(ims:ime,jms:jme),OV(ims:ime,jms:jme),UN(ims:ime,jms:jme)
      REAL MXM(ims:ime,jms:jme),MN(ims:ime,jms:jme)
      REAL FL(ims:ime,jms:jme,0:1)
      REAL XIG(NF*NF), XJG(NF*NF)
      integer rr
      REAL rioff, rjoff
      REAL donor, y1, y2, a
      DONOR(Y1,Y2,A)=(Y1*AMAX1(0.,SIGN(1.,A))-Y2*AMIN1(0.,SIGN(1.,A)))*A
      REAL tr4, ym1, y0, yp1, yp2
      TR4(YM1,Y0,YP1,YP2,A)=A*ONE12*(7.*(YP1+Y0)-(YP2+YM1)) &
       -A*A*ONE24*(15.*(YP1-Y0)-(YP2-YM1))-A*A*A*ONE12*((YP1+Y0) &
       -(YP2+YM1))+A*A*A*A*ONE24*(3.*(YP1-Y0)-(YP2-YM1))
      REAL pp, pn, x
      PP(X)=AMAX1(0.,X)
      PN(X)=AMIN1(0.,X)
      rr = nint(sqrt(float(nf)))
      rioff = 0
      rjoff = 0
      if(xstag .and. (mod(rr,2) .eq. 0)) rioff = 1.
      if(ystag .and. (mod(rr,2) .eq. 0)) rjoff = 1.
      DO I=1,rr
        DO J=1,rr
          XIG(J+(I-1)*rr)=(float(rr)-1.-rioff)/float(2*rr)-FLOAT(J-1)*1./float(rr)
          XJG(J+(I-1)*rr)=(float(rr)-1.-rjoff)/float(2*rr)-FLOAT(I-1)*1./float(rr)
        ENDDO
      ENDDO
      N2STAR = jts
      N2END = jte
      N1STAR = its
      N1END = ite
      DO 2000 IIM=1,NF
        DO 9000 JJ=N2STAR,N2END
          DO 50 J=-IOR,IOR
            DO 51 I=-IOR,IOR
              DO 511 II=N1STAR,N1END
                IF ( icmask(II,JJ) ) Y(II,JJ,I)=XF(II+I,JJ+J,IIM)
  511 CONTINUE
   51 CONTINUE
            DO 811 II=N1STAR,N1END
              IF ( icmask(II,JJ) ) THEN
                FL(II,JJ,0)=DONOR(Y(II,JJ,-1),Y(II,JJ,0),XIG(IIM))
                FL(II,JJ,1)=DONOR(Y(II,JJ,0),Y(II,JJ,1),XIG(IIM))
              ENDIF
  811 CONTINUE
            DO 812 II=N1STAR,N1END
              IF ( icmask(II,JJ) ) W(II,JJ)=Y(II,JJ,0)-(FL(II,JJ,1)-FL(II,JJ,0))
  812 CONTINUE
            DO 813 II=N1STAR,N1END
              IF ( icmask(II,JJ) ) THEN
                MXM(II,JJ)= &
                         AMAX1(Y(II,JJ,-1),Y(II,JJ,0),Y(II,JJ,1), &
                         W(II,JJ))
                MN(II,JJ)=AMIN1(Y(II,JJ,-1),Y(II,JJ,0),Y(II,JJ,1),W(II,JJ))
              ENDIF
  813 CONTINUE
            DO 312 II=N1STAR,N1END
              IF ( icmask(II,JJ) ) THEN
                F(II,JJ,0)= &
                           TR4(Y(II,JJ,-2),Y(II,JJ,-1),Y(II,JJ,0), &
                           Y(II,JJ,1),XIG(IIM))
                F(II,JJ,1)= &
                         TR4(Y(II,JJ,-1),Y(II,JJ,0),Y(II,JJ,1),Y(II,JJ,2),&
                         XIG(IIM))
                ENDIF
  312 CONTINUE
            DO 822 II=N1STAR,N1END
              IF ( icmask(II,JJ) ) THEN
                F(II,JJ,0)=F(II,JJ,0)-FL(II,JJ,0)
                F(II,JJ,1)=F(II,JJ,1)-FL(II,JJ,1)
              ENDIF
  822 CONTINUE
            DO 823 II=N1STAR,N1END
              IF ( icmask(II,JJ) ) THEN
                OV(II,JJ)=(MXM(II,JJ)-W(II,JJ))/(-PN(F(II,JJ,1))+ &
                        PP(F(II,JJ,0))+EP)
                UN(II,JJ)=(W(II,JJ)-MN(II,JJ))/(PP(F(II,JJ,1))- &
                      PN(F(II,JJ,0))+EP)
              ENDIF
  823 CONTINUE
            DO 824 II=N1STAR,N1END
              IF ( icmask(II,JJ) ) THEN
                F(II,JJ,0)=PP(F(II,JJ,0))*AMIN1(1.,OV(II,JJ))+ &
                           PN(F(II,JJ,0))*AMIN1(1.,UN(II,JJ))
                F(II,JJ,1)=PP(F(II,JJ,1))*AMIN1(1.,UN(II,JJ))+ &
                           PN(F(II,JJ,1))*AMIN1(1.,OV(II,JJ))
              ENDIF
  824 CONTINUE
            DO 825 II=N1STAR,N1END
              IF ( icmask(II,JJ) ) THEN
                Y(II,JJ,0)=W(II,JJ)-(F(II,JJ,1)-F(II,JJ,0))
              ENDIF
  825 CONTINUE
            DO 361 II=N1STAR,N1END
              IF ( icmask(II,JJ) ) Z(II,JJ,J)=Y(II,JJ,0)
  361 CONTINUE
 8000 CONTINUE
   50 CONTINUE
          DO 911 II=N1STAR,N1END
            IF ( icmask(II,JJ) ) THEN
              FL(II,JJ,0)=DONOR(Z(II,JJ,-1),Z(II,JJ,0),XJG(IIM))
              FL(II,JJ,1)=DONOR(Z(II,JJ,0),Z(II,JJ,1),XJG(IIM))
            ENDIF
  911 CONTINUE
          DO 912 II=N1STAR,N1END
            IF ( icmask(II,JJ) ) W(II,JJ)=Z(II,JJ,0)-(FL(II,JJ,1)-FL(II,JJ,0))
  912 CONTINUE
          DO 913 II=N1STAR,N1END
            IF ( icmask(II,JJ) ) THEN
              MXM(II,JJ)=AMAX1(Z(II,JJ,-1),Z(II,JJ,0),Z(II,JJ,1),W(II,JJ))
              MN(II,JJ)=AMIN1(Z(II,JJ,-1),Z(II,JJ,0),Z(II,JJ,1),W(II,JJ))
            ENDIF
  913 CONTINUE
          DO 412 II=N1STAR,N1END
            IF ( icmask(II,JJ) ) THEN
              F(II,JJ,0)= &
                         TR4(Z(II,JJ,-2),Z(II,JJ,-1),Z(II,JJ,0),Z(II,JJ,1)&
                         ,XJG(IIM))
              F(II,JJ,1)= &
                         TR4(Z(II,JJ,-1),Z(II,JJ,0),Z(II,JJ,1),Z(II,JJ,2), &
                         XJG(IIM))
            ENDIF
  412 CONTINUE
          DO 922 II=N1STAR,N1END
            IF ( icmask(II,JJ) ) THEN
              F(II,JJ,0)=F(II,JJ,0)-FL(II,JJ,0)
              F(II,JJ,1)=F(II,JJ,1)-FL(II,JJ,1)
            ENDIF
  922 CONTINUE
          DO 923 II=N1STAR,N1END
            IF ( icmask(II,JJ) ) THEN
              OV(II,JJ)=(MXM(II,JJ)-W(II,JJ))/(-PN(F(II,JJ,1))+ &
                        PP(F(II,JJ,0))+EP)
              UN(II,JJ)=(W(II,JJ)-MN(II,JJ))/(PP(F(II,JJ,1))-PN(F(II,JJ,0))+ &
                      EP)
            ENDIF
  923 CONTINUE
          DO 924 II=N1STAR,N1END
            IF ( icmask(II,JJ) ) THEN
              F(II,JJ,0)=PP(F(II,JJ,0))*AMIN1(1.,OV(II,JJ))+PN(F(II,JJ,0)) &
                         *AMIN1(1.,UN(II,JJ))
              F(II,JJ,1)=PP(F(II,JJ,1))*AMIN1(1.,UN(II,JJ))+PN(F(II,JJ,1)) &
                         *AMIN1(1.,OV(II,JJ))
            ENDIF
  924 CONTINUE
 9000 CONTINUE
        DO 925 JJ=N2STAR,N2END
          DO 925 II=N1STAR,N1END
            IF ( icmask(II,JJ) ) XF(II,JJ,IIM)=W(II,JJ)-(F(II,JJ,1)-F(II,JJ,0))
  925 CONTINUE
 2000 CONTINUE
      RETURN
      END
      SUBROUTINE SINTB(XF1, XF , &
                   ims, ime, jms, jme, icmask , &
                   its, ite, jts, jte, nf, xstag, ystag )
      IMPLICIT NONE
      INTEGER ims, ime, jms, jme, &
              its, ite, jts, jte
      LOGICAL icmask( ims:ime, jms:jme )
      LOGICAL xstag, ystag
      INTEGER nf, ior
      REAL one12, one24, ep
      PARAMETER(one12=1./12.,one24=1./24.)
      PARAMETER(ior=2)
      REAL XF(ims:ime,jms:jme,NF)
      REAL XF1(ims:ime,jms:jme,NF)
      REAL Y(ims:ime,jms:jme,-IOR:IOR), &
           Z(ims:ime,jms:jme,-IOR:IOR), &
           F(ims:ime,jms:jme,0:1)
      INTEGER I,J,II,JJ,IIM
      INTEGER N2STAR, N2END, N1STAR, N1END
      DATA EP/ 1.E-10/
      REAL W(ims:ime,jms:jme),OV(ims:ime,jms:jme),UN(ims:ime,jms:jme)
      REAL MXM(ims:ime,jms:jme),MN(ims:ime,jms:jme)
      REAL FL(ims:ime,jms:jme,0:1)
      REAL XIG(NF*NF), XJG(NF*NF)
      integer rr
      REAL rioff, rjoff
      REAL donor, y1, y2, a
      DONOR(Y1,Y2,A)=(Y1*AMAX1(0.,SIGN(1.,A))-Y2*AMIN1(0.,SIGN(1.,A)))*A
      REAL tr4, ym1, y0, yp1, yp2
      TR4(YM1,Y0,YP1,YP2,A)=A*ONE12*(7.*(YP1+Y0)-(YP2+YM1)) &
       -A*A*ONE24*(15.*(YP1-Y0)-(YP2-YM1))-A*A*A*ONE12*((YP1+Y0) &
       -(YP2+YM1))+A*A*A*A*ONE24*(3.*(YP1-Y0)-(YP2-YM1))
      REAL pp, pn, x
      PP(X)=AMAX1(0.,X)
      PN(X)=AMIN1(0.,X)
      rr = nint(sqrt(float(nf)))
      rioff = 0
      rjoff = 0
      if(xstag .and. (mod(rr,2) .eq. 0)) rioff = 1.
      if(ystag .and. (mod(rr,2) .eq. 0)) rjoff = 1.
      DO I=1,rr
        DO J=1,rr
          XIG(J+(I-1)*rr)=(float(rr)-1.-rioff)/float(2*rr)-FLOAT(J-1)*1./float(rr)
          XJG(J+(I-1)*rr)=(float(rr)-1.-rjoff)/float(2*rr)-FLOAT(I-1)*1./float(rr)
        ENDDO
      ENDDO
      N2STAR = jts
      N2END = jte
      N1STAR = its
      N1END = ite
      DO 2000 IIM=1,NF
        DO 9000 JJ=N2STAR,N2END
          DO 50 J=-IOR,IOR
            DO 51 I=-IOR,IOR
              DO 511 II=N1STAR,N1END
                Y(II,JJ,I)=XF1(II+I,JJ+J,IIM)
  511 CONTINUE
   51 CONTINUE
            DO 811 II=N1STAR,N1END
              FL(II,JJ,0)=DONOR(Y(II,JJ,-1),Y(II,JJ,0),XIG(IIM))
              FL(II,JJ,1)=DONOR(Y(II,JJ,0),Y(II,JJ,1),XIG(IIM))
  811 CONTINUE
            DO 812 II=N1STAR,N1END
              W(II,JJ)=Y(II,JJ,0)-(FL(II,JJ,1)-FL(II,JJ,0))
  812 CONTINUE
            DO 813 II=N1STAR,N1END
              MXM(II,JJ)= &
                       AMAX1(Y(II,JJ,-1),Y(II,JJ,0),Y(II,JJ,1), &
                       W(II,JJ))
              MN(II,JJ)=AMIN1(Y(II,JJ,-1),Y(II,JJ,0),Y(II,JJ,1),W(II,JJ))
  813 CONTINUE
            DO 312 II=N1STAR,N1END
              F(II,JJ,0)= &
                         TR4(Y(II,JJ,-2),Y(II,JJ,-1),Y(II,JJ,0), &
                         Y(II,JJ,1),XIG(IIM))
              F(II,JJ,1)= &
                       TR4(Y(II,JJ,-1),Y(II,JJ,0),Y(II,JJ,1),Y(II,JJ,2),&
                       XIG(IIM))
  312 CONTINUE
            DO 822 II=N1STAR,N1END
              F(II,JJ,0)=F(II,JJ,0)-FL(II,JJ,0)
              F(II,JJ,1)=F(II,JJ,1)-FL(II,JJ,1)
  822 CONTINUE
            DO 823 II=N1STAR,N1END
              OV(II,JJ)=(MXM(II,JJ)-W(II,JJ))/(-PN(F(II,JJ,1))+ &
                      PP(F(II,JJ,0))+EP)
              UN(II,JJ)=(W(II,JJ)-MN(II,JJ))/(PP(F(II,JJ,1))- &
                    PN(F(II,JJ,0))+EP)
  823 CONTINUE
            DO 824 II=N1STAR,N1END
              F(II,JJ,0)=PP(F(II,JJ,0))*AMIN1(1.,OV(II,JJ))+ &
                         PN(F(II,JJ,0))*AMIN1(1.,UN(II,JJ))
              F(II,JJ,1)=PP(F(II,JJ,1))*AMIN1(1.,UN(II,JJ))+ &
                         PN(F(II,JJ,1))*AMIN1(1.,OV(II,JJ))
  824 CONTINUE
            DO 825 II=N1STAR,N1END
              Y(II,JJ,0)=W(II,JJ)-(F(II,JJ,1)-F(II,JJ,0))
  825 CONTINUE
            DO 361 II=N1STAR,N1END
              Z(II,JJ,J)=Y(II,JJ,0)
  361 CONTINUE
 8000 CONTINUE
   50 CONTINUE
          DO 911 II=N1STAR,N1END
            FL(II,JJ,0)=DONOR(Z(II,JJ,-1),Z(II,JJ,0),XJG(IIM))
            FL(II,JJ,1)=DONOR(Z(II,JJ,0),Z(II,JJ,1),XJG(IIM))
  911 CONTINUE
          DO 912 II=N1STAR,N1END
            W(II,JJ)=Z(II,JJ,0)-(FL(II,JJ,1)-FL(II,JJ,0))
  912 CONTINUE
          DO 913 II=N1STAR,N1END
            MXM(II,JJ)=AMAX1(Z(II,JJ,-1),Z(II,JJ,0),Z(II,JJ,1),W(II,JJ))
            MN(II,JJ)=AMIN1(Z(II,JJ,-1),Z(II,JJ,0),Z(II,JJ,1),W(II,JJ))
  913 CONTINUE
          DO 412 II=N1STAR,N1END
            F(II,JJ,0)= &
                       TR4(Z(II,JJ,-2),Z(II,JJ,-1),Z(II,JJ,0),Z(II,JJ,1)&
                       ,XJG(IIM))
            F(II,JJ,1)= &
                       TR4(Z(II,JJ,-1),Z(II,JJ,0),Z(II,JJ,1),Z(II,JJ,2), &
                       XJG(IIM))
  412 CONTINUE
          DO 922 II=N1STAR,N1END
            F(II,JJ,0)=F(II,JJ,0)-FL(II,JJ,0)
            F(II,JJ,1)=F(II,JJ,1)-FL(II,JJ,1)
  922 CONTINUE
          DO 923 II=N1STAR,N1END
            OV(II,JJ)=(MXM(II,JJ)-W(II,JJ))/(-PN(F(II,JJ,1))+ &
                      PP(F(II,JJ,0))+EP)
            UN(II,JJ)=(W(II,JJ)-MN(II,JJ))/(PP(F(II,JJ,1))-PN(F(II,JJ,0))+ &
                    EP)
  923 CONTINUE
          DO 924 II=N1STAR,N1END
            F(II,JJ,0)=PP(F(II,JJ,0))*AMIN1(1.,OV(II,JJ))+PN(F(II,JJ,0)) &
                       *AMIN1(1.,UN(II,JJ))
            F(II,JJ,1)=PP(F(II,JJ,1))*AMIN1(1.,UN(II,JJ))+PN(F(II,JJ,1)) &
                       *AMIN1(1.,OV(II,JJ))
  924 CONTINUE
 9000 CONTINUE
        DO 925 JJ=N2STAR,N2END
          DO 925 II=N1STAR,N1END
            XF(II,JJ,IIM)=W(II,JJ)-(F(II,JJ,1)-F(II,JJ,0))
  925 CONTINUE
 2000 CONTINUE
      RETURN
      END