!                         WD23JP
!                                                         *****************
!                                                         *   SOCT2.FOR   *
!                                                         *  PURSER 1996  *
!                                                         *****************
!
!------------------------------------------------------------------------------C
!   R.J.Purser, NCEP, Washington D.C.   1996                                   C
!                   SUBROUTINE RELSUBP                                         C
!  Perform one sweep of the relaxations of the variational problem within the  C
!  first octant (0 to 45 degrees from x axis), or subvert the code to compute  C
!  residuals instead. Everything done in double precision.                     C
!                                                                              C
! <-> X     estimated cartesian locations for each octogon grid point          C
! <-> FORX  forcings associated with locations X                               C
! <-> DELX  residuals associated with X                                        C
! --> NA    steps between octagon corner and center line                       C
! --> N     steps between octagon center and edge                              C
! --> ORC   overrelaxation coefficient (typically 1.2 or 1.3)                  C
! --> AODD  homogeneity parameter normalized by step size                      C
! <-- ACC   accumulated rms norm                                               C
! --> MODE  mode of application: 0 => standard relaxation iteration            C
!                                1 => get residuals, put into DELX             C
!                                2 => adjust forcings, put into FORX           C
!------------------------------------------------------------------------------C
      SUBROUTINE RELSUBP(X,FORX,DELX,NA,N,ORC,AODD,ACC,MODE)
      DIMENSION X(3,-1:N+1,-1:N+1),DELX(3,-1:N+1,-1:N+1)
     *,FORX(3,-1:N+1,-1:N+1),XCOR0(3,3),XCOR1(3,3)
      NM=N-1
      NP=N+1
      NA4=NA+N
      NA2=(NA4+1)/2
      NAM=NA-1
      NAP=NA+1

      NACC=0
      ACC=0.

!  STANDARD RELAXATIONS APPLIED TO LOWER-BAND OF PRINCIPAL SUB-PANEL..
      DO IY=0,NA-1
       IYM=IY-1
       IYP=IY+1
       DO IX=IY,N
        IXM=IX-1
        IXP=IX+1
        CALL RELAX(FORX(1,IX,IY),DELX(1,IX,IY)
     *   ,X(1,IX,IY),X(1,IXP,IY),X(1,IX,IYP),X(1,IXM,IY)
     *   ,X(1,IX,IYM),X(1,IXP,IYP),X(1,IXM,IYP),X(1,IXM,IYM)
     *   ,X(1,IXP,IYM),ORC,AODD,ACC,NACC,MODE)
       ENDDO
      ENDDO
!                        ..APPLIED TO MIDLINE..
      IY=NA
       IYM=IY-1
       IYP=IY+1
       DO IX=IY,NM
        IXM=IX-1
        IXP=IX+1
        CALL RELAX(FORX(1,IX,IY),DELX(1,IX,IY)
     *   ,X(1,IX,IY),X(1,IXP,IY),X(1,IX,IYP),X(1,IXM,IY)
     *   ,X(1,IX,IYM),X(1,IXP,IYP),X(1,IXM,IYP),X(1,IXM,IYM)
     *   ,X(1,IXP,IYM),ORC,AODD,ACC,NACC,MODE)
       ENDDO
!                        ..APPLIED TO UPPER BAND:
      DO IY=NA+1,NA2
       IYM=IY-1
       IYP=IY+1
       DO IX=IY,NA4-IY
        IXM=IX-1
        IXP=IX+1
        CALL RELAX(FORX(1,IX,IY),DELX(1,IX,IY)
     *   ,X(1,IX,IY),X(1,IXP,IY),X(1,IX,IYP),X(1,IXM,IY)
     *   ,X(1,IX,IYM),X(1,IXP,IYP),X(1,IXM,IYP),X(1,IXM,IYM)
     *   ,X(1,IXP,IYM),ORC,AODD,ACC,NACC,MODE)
       ENDDO
      ENDDO
      DO I=1,3
       XCOR0(I,1)=X(I,N,NAM)
       XCOR1(I,1)=X(I,NM,NAM)
       XCOR0(I,3)=X(I,NM,NA)
       XCOR1(I,3)=X(I,NM,NAP)
       XCOR0(I,2)=X(I,NP,NA)
       XCOR1(I,2)=X(I,NP,NAM)
      ENDDO
      CALL RELV(FORX(1,N,NA),DELX(1,N,NA)
     *,X(1,N,NA),XCOR0,XCOR1,3,ORC,AODD,ACC,NACC,MODE)
      IF(MODE.EQ.0)THEN
       ACC=ACC/NACC
       ACC=SQRT(ACC)
       CALL SYMMET(X,NA,N)
      ELSEIF(MODE.EQ.1)THEN
       CALL SYMMET(DELX,NA,N)
      ELSE
       CALL SYMMET(FORX,NA,N)
      ENDIF
      RETURN
      END

!------------------------------------------------------------------------------C
!   R.J.Purser, NCEP, Washington D.C.   1996                                   C
!                   SUBROUTINE RELAX                                           C
!  Relaxation (MODE.eq.0), or residual determination (MODE.eq.1), or           C
!  adjustment of forcing (MODE.eq.2), at a generic grid point.                 C
!                                                                              C
!  <-> FORX,DELX,X:  vectors of forcing, residual and solution at this point   C
!                                                                              C
!      XF XB XE:                                                               C
!  --> XC    XA:     surrounding solution vectors in the arrangement shown     C
!      XG XD XH:                                                               C
!                                                                              C
! --> ORC:           overrelaxation cooeficient                                C
! --> AODD:          homogeneity parameter rescaled to grid units              C
! --> ACC:           accumulated squared residual                              C
! --> NACC:          counter of contributors to ACC                            C
! --> MODE:          mode of application                                       C
!------------------------------------------------------------------------------C
      SUBROUTINE RELAX(FORX,DELX,X
     *,XA,XB,XC,XD,XE,XF,XG,XH,ORC,AODD,ACC,NACC,MODE)
      PARAMETER (DELTA=.3E-4,DELTAI=1./DELTA)
      DIMENSION FORX(3),DELX(3),X(3),XA(3),XB(3),XC(3),XD(3)
     *,XE(3),XF(3),XG(3),XH(3),XU(3),XV(3),U(3),V(3),XS(3)
      SX=0.
      SY=0.
      DO I=1,3
       U(I)=(XA(I)-XC(I))*.5D0
       V(I)=(XB(I)-XD(I))*.5D0
       SX=SX+X(I)*U(I)
       SY=SY+X(I)*V(I)
       XS(I)=X(I)
      ENDDO
!misha
!       print *,'=================  XC   =',XC
!       print *,'=================  XD   =',XD
!misha
      DO I=1,3
       U(I)=U(I)-X(I)*SX
       V(I)=V(I)-X(I)*SY
      ENDDO
      SX=0.
      SY=0.
      DO I=1,3
       XU(I)=X(I)+U(I)*DELTA
       XV(I)=X(I)+V(I)*DELTA
       SX=SX+XU(I)**2
       SY=SY+XV(I)**2
      ENDDO
      SX=1./SQRT(SX)
      SY=1./SQRT(SY)
      DO I=1,3
       XU(I)=XU(I)*SX
       XV(I)=XV(I)*SY
      ENDDO
      CALL REL(X ,XA,XB,XC,XD,XE,XF,XG,XH,U,V,AODD,DX,DY)
      CALL REL(XU,XA,XB,XC,XD,XE,XF,XG,XH,U,V,AODD,DXX,DYX)
      CALL REL(XV,XA,XB,XC,XD,XE,XF,XG,XH,U,V,AODD,DXY,DYY)
      DXX=(DXX-DX)*DELTAI
      DXY=(DXY-DY+DYX-DX)*DELTAI*.5D0
      DYY=(DYY-DY)*DELTAI
      DET=DXX*DYY-DXY*DXY
!      IF(DET.EQ.0.D0)CALL WHAT(X,XA,XB,XC,XD,XE,XF,XG,XH)
      DETI=1./DET
      SX=(DYY*DX-DXY*DY)*DETI
      SY=(DXX*DY-DXY*DX)*DETI
      IF(MODE.EQ.0)THEN       ! RELAXATION WITH LAGRANGE MULTIPLIER DELX
      S=0.
      DO I=1,3
       X(I)=X(I)+(FORX(I)-U(I)*SX-V(I)*SY)*ORC
       S=S+X(I)**2
      ENDDO
      S=1./SQRT(S)
      DO I=1,3
       X(I)=X(I)*S
      ENDDO
      DO I=1,3
        ACC=ACC+(XS(I)-X(I))**2
      ENDDO
      ELSEIF(MODE.EQ.1)THEN  ! DIAGNOSE RESIDUAL WITHOUT CORRECTING
      DO I=1,3
       DELX(I)=(FORX(I)-U(I)*SX-V(I)*SY)
       ACC=ACC+DELX(I)**2
      ENDDO
      ELSEIF(MODE.EQ.2)THEN ! RESET LAGRANGE MULTIPLIER
      DO I=1,3
       ACC=ACC+FORX(I)**2
       FORX(I)=FORX(I)+U(I)*SX+V(I)*SY
      ENDDO
      ELSE
        STOP'INVALID MODE IN RELAX'
      ENDIF
      NACC=NACC+1

      RETURN
      END

!------------------------------------------------------------------------------C
!   R.J.Purser, NCEP, Washington D.C.   1996                                   C
!                   SUBROUTINE RELV                                            C
!  Relaxation (MODE.eq.0), or residual determination (MODE.eq.1), or           C
!  adjustment of forcing (MODE.eq.2), at a vertex of the octagon               C
!                                                                              C
!  <-> FORX,DELX,X:  vectors of forcing, residual and solution at this point   C
!  --> XCOR0: X at the NC nearest surrounding grid points                      C
!  --> XCOR1: X at the NC next nearest surrounding points                      C
!  --> NC:    number of panels coming together at this corner (=3 for octagon) C
!  --> ORC:   over relaxation coefficient                                      C
! --> AODD:          homogeneity parameter rescaled to grid units              C
! --> ACC:           accumulated squared residual                              C
! --> NACC:          counter of contributors to ACC                            C
! --> MODE:          mode of application                                       C
!------------------------------------------------------------------------------C
      SUBROUTINE RELV(FORX,DELX,X
     *,XCOR0,XCOR1,NC,ORC,AODD,ACC,NACC,MODE)
      PARAMETER (DELTA=.3E-4,DELTAI=1./DELTA)
      DIMENSION FORX(3),DELX(3),X(3),XCOR0(3,3),XCOR1(3,3)
     *,XU(3),XV(3),U(3),V(3),XS(3)
      SX=0.
      SY=0.
      DO I=1,3
       U(I)=XCOR0(I,1)-X(I)
       V(I)=XCOR1(I,1)-XCOR0(I,1)
       SX=SX+X(I)*U(I)
       SY=SY+X(I)*V(I)
       XS(I)=X(I)
      ENDDO
      DO I=1,3
       U(I)=U(I)-X(I)*SX
       V(I)=V(I)-X(I)*SY
      ENDDO

      SX=0.
      SY=0.
      DO I=1,3
       XU(I)=X(I)+U(I)*DELTA
       XV(I)=X(I)+V(I)*DELTA
       SX=SX+XU(I)**2
       SY=SY+XV(I)**2
      ENDDO
      SX=1./SQRT(SX)
      SY=1./SQRT(SY)
      DO I=1,3
       XU(I)=XU(I)*SX
       XV(I)=XV(I)*SY
      ENDDO
      CALL DELV(X,XCOR0,XCOR1,NC
     *,U,V,AODD,DX,DY)
      CALL DELV(XU,XCOR0,XCOR1,NC
     *,U,V,AODD,DXX,DYX)
      CALL DELV(XV,XCOR0,XCOR1,NC
     *,U,V,AODD,DXY,DYY)
      DXX=(DXX-DX)*DELTAI
      DXY=(DXY-DY+DYX-DX)*DELTAI*.5D0
      DYY=(DYY-DY)*DELTAI
      DETI=1./(DXX*DYY-DXY*DXY)
      SX=(DYY*DX-DXY*DY)*DETI
      SY=(DXX*DY-DXY*DX)*DETI
      IF(MODE.EQ.0)THEN       ! RELAXATION WITH LAGRANGE MULTIPLIER DELX
      S=0.
      DO I=1,3
       X(I)=X(I)+(FORX(I)-U(I)*SX-V(I)*SY)*ORC
       S=S+X(I)**2
      ENDDO
      S=1./SQRT(S)
      DO I=1,3
       X(I)=X(I)*S
      ENDDO
      DO I=1,3
        ACC=ACC+(XS(I)-X(I))**2
      ENDDO
      ELSEIF(MODE.EQ.1)THEN  ! DIAGNOSE RESIDUAL WITHOUT CORRECTING
      DO I=1,3
       DELX(I)=(FORX(I)-U(I)*SX-V(I)*SY)
       ACC=ACC+DELX(I)**2
      ENDDO
      ELSEIF(MODE.EQ.2)THEN ! RESET LAGRANGE MULTIPLIER
      DO I=1,3
       ACC=ACC+FORX(I)**2
       FORX(I)=FORX(I)+U(I)*SX+V(I)*SY
      ENDDO
      ELSE
        STOP'INVALID MODE IN RELAX'
      ENDIF
      NACC=NACC+1
      RETURN
      END

      SUBROUTINE DELV(X,XCOR0,XCOR1,NC,U,V,AODD,DX,DY)
      DIMENSION X(3),XCOR0(3,*),XCOR1(3,*),U(3),V(3)
      DX=0.
      DY=0.
      DO IC=1,NC
       ICP=IC+1
       IF(ICP.GT.NC)ICP=1
       UE=0.
       UAB=0.
       VAB=0.
       VE=0.
       SAB=0.
       SE=0.
       DO I=1,3
        UI=U(I)
        VI=V(I)
        XI=X(I)
        D1=XCOR1(I,IC)-XI
        D2=XCOR0(I,ICP)-XCOR1(I,IC)
        D3=XCOR1(I,IC)-XI
        DX=DX+UI*D3
        DY=DY+VI*D3
        UE=UE+UI*D1
        UAB=UAB+UI*D2
        VE=VE+VI*D1
        VAB=VAB+VI*D2
        SAB=SAB+D2*D2
        SE=SE+D1*D2
       ENDDO
       DX=DX+AODD*(UE*SAB-UAB*SE)
       DY=DY+AODD*(VE*SAB-VAB*SE)
      ENDDO
      RETURN
      END

      SUBROUTINE REL(X,XA,XB,XC,XD,XE,XF,XG,XH,U,V,AODD,DX,DY)
!
!           XF  XB  XE
!             QF  QE
!           XC  X   XA
!             QG  QH
!           XG  XD  XH
!
      DIMENSION X(3),XA(3),XB(3),XC(3),XD(3),XE(3),XF(3),XG(3),XH(3)
     *,U(3),V(3)
      DX=0.
      DY=0.
      UAB=0.
      UBC=0.
      UCD=0.
      UDA=0.
      UE=0.
      UF=0.
      UG=0.
      UH=0.
      VAB=0.
      VBC=0.
      VCD=0.
      VDA=0.
      VE=0.
      VF=0.
      VG=0.
      VH=0.
      SAB=0.
      SBC=0.
      SCD=0.
      SDA=0.
      SE=0.
      SF=0.
      SG=0.
      SH=0.
      DO I=1,3
       UI=U(I)
       VI=V(I)
       XI=X(I)
       D1=XE(I)-XI
       D2=XB(I)-XA(I)
       D3=XA(I)-XI
       UE=UE+UI*D1
       UAB=UAB+UI*D2
       DX=DX+UI*D3
       VE=VE+VI*D1
       VAB=VAB+VI*D2
       DY=DY+VI*D3
       SAB=SAB+D2*D2
       SE=SE+D1*D2

       D1=XF(I)-XI
       D2=XC(I)-XB(I)
       D3=XB(I)-XI
       UF=UF+UI*D1
       UBC=UBC+UI*D2
       DX=DX+UI*D3
       VF=VF+VI*D1
       VBC=VBC+VI*D2
       DY=DY+VI*D3
       SBC=SBC+D2*D2
       SF=SF+D1*D2

       D1=XG(I)-XI
       D2=XD(I)-XC(I)
       D3=XC(I)-XI
       UG=UG+UI*D1
       UCD=UCD+UI*D2
       DX=DX+UI*D3
       VG=VG+VI*D1
       VCD=VCD+VI*D2
       DY=DY+VI*D3
       SCD=SCD+D2*D2
       SG=SG+D1*D2

       D1=XH(I)-XI
       D2=XA(I)-XD(I)
       D3=XD(I)-XI
       UH=UH+UI*D1
       UDA=UDA+UI*D2
       DX=DX+UI*D3
       VH=VH+VI*D1
       VDA=VDA+VI*D2
       DY=DY+VI*D3
       SDA=SDA+D2*D2
       SH=SH+D1*D2
      ENDDO
      DX=DX+AODD*(UE*SAB+UF*SBC+UG*SCD+UH*SDA
     *           -UAB*SE-UBC*SF-UCD*SG-UDA*SH)
      DY=DY+AODD*(VE*SAB+VF*SBC+VG*SCD+VH*SDA
     *           -VAB*SE-VBC*SF-VCD*SG-VDA*SH)
      RETURN
      END

      SUBROUTINE SYMMET(X,NA,N)
      DIMENSION X(3,-1:N+1,-1:N+1)
      DIMENSION LOFIA(3),LOFIC(3),LOFID(3),LOFIE(3)
     *,JOFIF(3)
      DATA LOFIA/1,1,-1/,LOFIC/-1,1,1/,LOFID/1,-1,1/,LOFIE/1,1,-1/
      DATA JOFIF/2,1,3/
      NM=N-1
      NP=N+1
      NA4=NA+N
      NA2=(NA4+1)/2

!  SYMMETRIZE ON BOUNDARY OF SUB-PANEL:
      X(1,0,0)=0.
      X(2,0,0)=0.
      DO IX=1,N
       X(2,IX,0)=0.
      ENDDO
      DO IY=0,NA
       X(3,N,IY)=0.
      ENDDO
      DO IY=NA+1,NA2
       IX=NA4-IY
       X(3,IX,IY)=0.
      ENDDO
      DO IY=1,NA2
       IX=IY
       S=(X(1,IX,IY)+X(2,IX,IY))/2.D0
       X(1,IX,IY)=S
       X(2,IX,IY)=S
      ENDDO

!  SYMMETRIZE JUST OUTSIDE BOUNDARY OF SUB-PANEL:
      DO I=1,3
       L=LOFIC(I)
       DO IY=0,1
        X(I,-1,IY)=L*X(I,1,IY)
       ENDDO
       L=LOFIA(I)
       DO IY=0,NA
        X(I,NP,IY)=L*X(I,NM,IY)
       ENDDO
       J=JOFIF(I)
       DO IY=1,NA2
        IX=IY-1
        X(I,IX,IY)=X(J,IX+1,IY-1)
       ENDDO
       DO IY=2,NA2+1
        IX=IY-2
        X(I,IX,IY)=X(J,IX+2,IY-2)
       ENDDO
       L=LOFID(I)
       DO IX=-1,NP
        X(I,IX,-1)=L*X(I,IX,1)
       ENDDO
      ENDDO
      DO I=1,3
       L=LOFIE(I)
       DO IY=NA+1,NA2+1
        IX=NA4-IY+1
        X(I,IX,IY)=L*X(I,IX-1,IY-1)
       ENDDO
       DO IY=NA+1,NA2+1
        IX=NA4-IY+2
        X(I,IX,IY)=L*X(I,IX-2,IY-2)
       ENDDO
      ENDDO
      RETURN
      END

      SUBROUTINE ADDXX2(X,X2,NA,N)
      DIMENSION X(3,-1:N+1,-1:N+1),X2(3,-1:N*2+1,-1:N*2+1)
      NA4=NA+N
      NA2=(NA4+1)/2
      N2=N*2
      N2A=NA*2

      IY=0
       IY2=IY*2
       IY2M=IY2-1
       IY2P=IY2+1
       DO IX=0,N
        IX2=IX*2
        IX2M=IX2-1
        IX2P=IX2+1
        DO I=1,3
         R=X(I,IX,IY)
         X2(I,IX2,IY2)=X2(I,IX2,IY2)+R
         R=.5D0*R
         X2(I,IX2P,IY2)=X2(I,IX2P,IY2)+R
         X2(I,IX2,IY2P)=X2(I,IX2,IY2P)+R
         X2(I,IX2M,IY2)=X2(I,IX2M,IY2)+R
         X2(I,IX2,IY2M)=X2(I,IX2,IY2M)+R
         R=R*.5D0
         X2(I,IX2P,IY2P)=X2(I,IX2P,IY2P)+R
         X2(I,IX2M,IY2P)=X2(I,IX2M,IY2P)+R
         X2(I,IX2M,IY2M)=X2(I,IX2M,IY2M)+R
         X2(I,IX2P,IY2M)=X2(I,IX2P,IY2M)+R
        ENDDO
       ENDDO
      DO IY=1,NA
       IY2=IY*2
       IY2M=IY2-1
       IY2P=IY2+1
       DO IX=IY-1,N
        IX2=IX*2
        IX2M=IX2-1
        IX2P=IX2+1
        DO I=1,3
         R=X(I,IX,IY)
         X2(I,IX2,IY2)=X2(I,IX2,IY2)+R
         R=.5D0*R
         X2(I,IX2P,IY2)=X2(I,IX2P,IY2)+R
         X2(I,IX2,IY2P)=X2(I,IX2,IY2P)+R
         X2(I,IX2M,IY2)=X2(I,IX2M,IY2)+R
         X2(I,IX2,IY2M)=X2(I,IX2,IY2M)+R
         R=R*.5D0
         X2(I,IX2P,IY2P)=X2(I,IX2P,IY2P)+R
         X2(I,IX2M,IY2P)=X2(I,IX2M,IY2P)+R
         X2(I,IX2M,IY2M)=X2(I,IX2M,IY2M)+R
         X2(I,IX2P,IY2M)=X2(I,IX2P,IY2M)+R
        ENDDO
       ENDDO
      ENDDO
      DO IY=NA+1,NA2+1
       IY2=IY*2
       IY2M=IY2-1
       IY2P=IY2+1
       DO IX=IY-1,NA4-IY+1
        IX2=IX*2
        IX2M=IX2-1
        IX2P=IX2+1
        DO I=1,3
         R=X(I,IX,IY)
         X2(I,IX2,IY2)=X2(I,IX2,IY2)+R
         R=.5D0*R
         X2(I,IX2P,IY2)=X2(I,IX2P,IY2)+R
         X2(I,IX2,IY2P)=X2(I,IX2,IY2P)+R
         X2(I,IX2M,IY2)=X2(I,IX2M,IY2)+R
         X2(I,IX2,IY2M)=X2(I,IX2,IY2M)+R
         R=R*.5D0
         X2(I,IX2P,IY2P)=X2(I,IX2P,IY2P)+R
         X2(I,IX2M,IY2P)=X2(I,IX2M,IY2P)+R
         X2(I,IX2M,IY2M)=X2(I,IX2M,IY2M)+R
         X2(I,IX2P,IY2M)=X2(I,IX2P,IY2M)+R
        ENDDO
       ENDDO
      ENDDO
      CALL FNORM(X2,N2A,N2)
      CALL SYMMET(X2,N2A,N2)
      RETURN
      END

      SUBROUTINE ADDX2X(X2,X,NA,N)
      DIMENSION X(3,-1:N+1,-1:N+1),X2(3,-1:N*2+1,-1:N*2+1)
      NA4=NA+N
      NA2=(NA4+1)/2
      N2=N*2
      N2A=NA*2
      N2M=N2-1
      N2P=N2+1
      N2AM=N2A-1
      N2AP=N2A+1
      DO IY=0,NA-1
       IY2=IY*2
       IY2M=IY2-1
       IY2P=IY2+1
       DO IX=IY,N
        IX2=IX*2
        IX2M=IX2-1
        IX2P=IX2+1
        DO I=1,3
         X(I,IX,IY)=X(I,IX,IY)+X2(I,IX2,IY2)+.5D0*(
     * X2(I,IX2P,IY2)+X2(I,IX2,IY2P)+X2(I,IX2M,IY2)+X2(I,IX2,IY2M)+.5*(
     * X2(I,IX2P,IY2P)+X2(I,IX2M,IY2P)+X2(I,IX2M,IY2M)+X2(I,IX2P,IY2M)
     *  ))
        ENDDO
       ENDDO
      ENDDO
      IY=NA
       IY2=IY*2
       IY2M=IY2-1
       IY2P=IY2+1
       DO IX=IY,N-1
        IX2=IX*2
        IX2M=IX2-1
        IX2P=IX2+1
        DO I=1,3
         X(I,IX,IY)=X(I,IX,IY)+X2(I,IX2,IY2)   +.5D0*(
     * X2(I,IX2P,IY2)+X2(I,IX2,IY2P)+X2(I,IX2M,IY2)+X2(I,IX2,IY2M)+.5*(
     *X2(I,IX2P,IY2P)+X2(I,IX2M,IY2P)+X2(I,IX2M,IY2M)+X2(I,IX2P,IY2M)))
        ENDDO
       ENDDO
      DO IY=NA+1,NA2
       IY2=IY*2
       IY2M=IY2-1
       IY2P=IY2+1
       DO IX=IY,NA4-IY
        IX2=IX*2
        IX2M=IX2-1
        IX2P=IX2+1
        DO I=1,3
         X(I,IX,IY)=X(I,IX,IY)+X2(I,IX2,IY2)   +.5D0*(
     * X2(I,IX2P,IY2)+X2(I,IX2,IY2P)+X2(I,IX2M,IY2)+X2(I,IX2,IY2M)+.5*(
     *X2(I,IX2P,IY2P)+X2(I,IX2M,IY2P)+X2(I,IX2M,IY2M)+X2(I,IX2P,IY2M)))
        ENDDO
       ENDDO
      ENDDO
      DO I=1,3
       X(I,N,NA)=X(I,N,NA)+X2(I,N2,N2A)
     * +.5D0*(X2(I,N2,N2AM)+X2(I,N2P,N2A)+X2(I,N2M,N2A)
     * +.5D0*(X2(I,N2M,N2AM)+X2(I,N2P,N2AM)+X2(I,N2M,N2AP)))
      ENDDO
      RETURN
      END

      SUBROUTINE SETX2X(X2,X,NA,N)
      DIMENSION X(3,-1:N+1,-1:N+1),X2(3,-1:N*2+1,-1:N*2+1)
      NA4=NA+N
      NA2=(NA4+1)/2
      DO IY=0,NA
       IY2=IY*2
       DO IX=IY,N
        IX2=IX*2
        DO I=1,3
         X(I,IX,IY)=X2(I,IX2,IY2)
        ENDDO
       ENDDO
      ENDDO
      DO IY=NA+1,NA2
       IY2=IY*2
       DO IX=IY,NA4-IY
        IX2=IX*2
        DO I=1,3
         X(I,IX,IY)=X2(I,IX2,IY2)
        ENDDO
       ENDDO
      ENDDO
      CALL SYMMET(X,NA,N)
      RETURN
      END

      SUBROUTINE DIFXX2(X,X2,NA,N)
      DIMENSION X(3,-1:N+1,-1:N+1),X2(3,-1:N*2+1,-1:N*2+1)
      NA4=NA+N
      NA2=(NA4+1)/2
      DO IY=0,NA
       IYM=IY-1
       IYP=IY+1
       IY2=IY*2
       DO IX=0,N
        IXM=IX-1
        IXP=IX+1
        IX2=IX*2
        DO I=1,3
         X(I,IX,IY)=X(I,IX,IY)-X2(I,IX2,IY2)
        ENDDO
       ENDDO
      ENDDO
      DO IY=NA+1,NA2
       IYM=IY-1
       IYP=IY+1
       IY2=IY*2
       DO IX=IY,NA4-IY
        IXM=IX-1
        IXP=IX+1
        IX2=IX*2
        DO I=1,3
         X(I,IX,IY)=X(I,IX,IY)-X2(I,IX2,IY2)
        ENDDO
       ENDDO
      ENDDO
      CALL SYMMET(X,NA,N)
      RETURN
      END

      SUBROUTINE ZERV(A,N)
      DIMENSION A(*)
      DO I=1,N
       A(I)=0.
      ENDDO
      RETURN
      END

      SUBROUTINE FNORM(X,NA,N)
      DIMENSION X(3,-1:N+1,-1:N+1)
      NA4=NA+N
      NA2=(NA4+1)/2
      DO IY=0,NA
       DO IX=IY,N
        CALL NORM3(X(1,IX,IY))
       ENDDO
      ENDDO
      DO IY=NA+1,NA2
       DO IX=IY,NA4-IY
        CALL NORM3(X(1,IX,IY))
       ENDDO
      ENDDO
      RETURN
      END

      SUBROUTINE NORM3(X)
!  NORMALIZE THE 3-VECTOR:
      DIMENSION X(3)
      S=0.
      DO I=1,3
       S=S+X(I)*X(I)
      ENDDO
      S=1./SQRT(S)
      DO I=1,3
       X(I)=X(I)*S
      ENDDO
      RETURN
      END


