
!                         WD23JP
!                                                         *****************
!                                                         *   SOCT1.FOR   *
!                                                         *  PURSER 1996  *
!                                                         *****************
!
!------------------------------------------------------------------------------C
!   R.J.Purser, NCEP, Washington D.C.   1996                                   C
!                   SUBROUTINE INSOCT1                                         C
!  Get conformal solution, then refine it to get "smoothed" grid solution      C
!  If using multigrid method, set parameter NGEN to positive number of         C
!  "generations" of grids. Otherwise, revert to SOR and ensure that enough     C
!  iterations (NITSOR), independent of the m.g. algorithm, are given.          C
!                                                                              C
! --> HOM       homogeneity parameter                                          C
! --- RA,QA     work arrays of size NFT                                        C
! --> NFT       suitable period (eg 128) for FFT                               C
! --- JUMBLE    integer work arrays of size NFT                                C
! --- WORK      real work array of size NFT                                    C
! --> TFILE     name of file into which table is put (CHARACTER*39)            C
! --> NCY       no. of multigrid (m.g.) cycles (eg 15 for m.g., 0 for S.O.R.)  C
! --> NIT       number of SOR iterations at each stage of m.g. (eg 3)          C
! --> NITSOR    number of final SOR iterations (10 for m.g., 1000 for SOR)     C
!------------------------------------------------------------------------------C
      SUBROUTINE INSOCT1(HOM,RA,QA,NFT,JUMBLE,WORK,TFILE
     *,NCY,NIT,NITSOR)
      PARAMETER(NGEN=6,NRAT1=1,NRAT2=3
     *,NXD=9*(NRAT2**2*(1+4*(4**NGEN-1)/3)+6*NRAT2*(1+2*(2**NGEN-1))
     * +9*(NGEN+1)))
      PARAMETER(N=NRAT2*2**NGEN,NP=N+1)
      COMPLEX CI,CIR,CIQ
      CHARACTER*39 TFILE
      COMMON/QPAN/X(3,-1:NP,-1:NP)
      COMMON/CSTOCT1/A1(40),A2(40),B1(40),B2(40),A,B,CX,CY
     *,SSR,N1,ST,STI,CI,CIR,CIQ
      DIMENSION RA(0:*),QA(0:*),JUMBLE(NFT),WORK(NFT)
!jim      DIMENSION RA(0:*),QA(0:*),JUMBLE(N),WORK(N)
      DIMENSION XD(NXD)
      A=FLOAT(NRAT1)/NRAT2
      N1=20
      CALL       INDOCT(RA,QA,NFT,JUMBLE,WORK)
      CALL IROCTA(XD,NIT,NITSOR,NCY,HOM,TFILE)
      RETURN
      END

!------------------------------------------------------------------------------C
!   R.J.Purser, NCEP, Washington D.C.   1996                                   C
!                   SUBROUTINE INDOCT                                          C
!  Initialize the Taylor series for the conformal octagon solution             C
!                                                                              C
! --- RA,QA     arrays of size N used for real & imag components of FFT        C
! --> N         suitable size (eg N=128) for one cycle of data in the FFT      C
! --- JUMBLE    integer work array of size N for FFT                           C
! --- WORK      real work array of size N for FFT                              C
!------------------------------------------------------------------------------C
      SUBROUTINE INDOCT(RA,QA,N,JUMBLE,WORK)
      COMPLEX CI,CIR,CIQ,CIA,CIAP,CIRS2,Z,W,CANG,CDANG1,CDANG2
      COMMON/CSTOCT1/A1(40),A2(40),B1(40),B2(40),A,B,CX,CY
     *,SSR,N1,ST,STI,CI,CIR,CIQ
      DIMENSION RA(0:*),QA(0:*),JUMBLE(N),WORK(N)
      ST=4./3.
      STI=3./4.
      CI=(0.,1.)
      CIR=(-CI)**.75
      CIQ=CI**.25
      NH=N/2
      R2=SQRT(2.)
      ORC=.5
      SW1=SQRT(1.+A*A)
      SW2=AMIN1(A*2,R2*(1.-A))
      SSR=(SW1/SW2)**2
      RS1=.95*SW1
      RS2=.95*SW2
      RS3=RS1**4
      RS4=RS2**ST

      DO I=1,N1
       A1(I)=0.
      ENDDO
      A1(1)= 1.05060
      A1(2)= -.172629
      A1(3)=  .125716
      A1(4)=  .00015657
      A1(5)= -.0017866
      A1(6)= -.0031168

      DO I=1,N1
       A2(I)=0.
      ENDDO
      A2(1)= .688411
      A2(2)=-.180924
      A2(3)=-.186593
      A2(4)= .024811
      A2(5)=-.044888
      A2(6)=-.049190
      B=-.147179

      JUMBLE(1)=0
      PIO4=ATAN(1.)
      CDANG1=CI*PIO4/NH
      CDANG2=-3.*CDANG1
      CIA=CI*A
      CIAP=CIA+1.
      CIRS2=CI*RS2

      DO KIT=1,50
       DO I=0,NH
        CANG=CDANG1*I
        Z=RS1*CEXP(CANG)
!mish        Z=RS1*EXP(CANG)
        CALL TOCT(Z,W)
        W=W**4
        U=REAL(W)
        V=AIMAG(W)
        RA(I)=U
        QA(I)=-V
        RA(N-I)=U
        QA(N-I)=V
       ENDDO
       QA(0)=0.
       QA(NH)=0.
       CALL CFFT(RA,QA,N,1.,WORK,JUMBLE)
       RS3P=1.
       DO I=1,N1
        RS3P=RS3P*RS3
        A1(I)=A1(I) +ORC*(RA(I)/RS3P-A1(I))
       ENDDO

      DO I=0,NH
       CANG=CDANG2*I
       Z=CIAP-CIRS2*CEXP(CANG)
       CALL TOCT(Z,W)
       W=CI*(W-1.)/(W+1.)
       U=REAL(W)
       V=AIMAG(W)
       RA(I)=U
       QA(I)=V
       RA(N-I)=U
       QA(N-I)=-V
      ENDDO
      QA(0)=0.
      QA(NH)=0.
      CALL CFFT(RA,QA,N,1.,WORK,JUMBLE)
      B=B+ORC*(RA(0)-B)
      RS4P=1.
      DO I=1,N1
       RS4P=RS4P*RS4
       A2(I)=A2(I) +ORC*(RA(I)/RS4P-A2(I))
      ENDDO

      ENDDO
      W=B
      W=(W+CI)/(CI-W)
      CX=REAL(W)
      CY=AIMAG(W)
      CALL CINVRT(A1,B1,N1,RA)
      CALL CINVRT(A2,B2,N1,RA)

      PRINT'('' IMAGE OF OCTAGON-VERTEX ON UNIT-CIRCLE:'')'
      WRITE(6,62)CX,CY
      PRINT'('' TAYLOR COEFFICIENTS, SERIES A1,A2,B1,B2:'')'
      WRITE(6,63)B
      DO I=1,N1
       WRITE(6,64)I,A1(I),A2(I),B1(I),B2(I)
      ENDDO
62    FORMAT(2(1X,E12.6))
63    FORMAT(4X,'0',13X,4(1X,E12.6))
64    FORMAT(1X,I4,4(1X,E12.6))
      RETURN
      END

!------------------------------------------------------------------------------C
!   R.J.Purser, NCEP, Washington D.C.   1996                                   C
!                   SUBROUTINE IROCTA                                          C
!  Solve the variational problem for positions of reference grid of principal  C
!  quadrant of standard smoothed octagon for the given homogeneity parameter.  C
!  Solution method is either nonlinear multigrid (NGEN.GT.0, NCY.GT.0) or else C
!  reduces to SOR when NGEN=0.                                                 C
!                                                                              C
! --- XD     double precision work array                                       C
! --> NIT    number of SOR iterations used at each stage of multigrid scheme   C
! --> NITSOR number of SOR iterations used at the end to polish final solution C
! --> NCY    number of multigrid "V-cycles" (ignored when NGEN=0)              C
! --> A      single precision homogeneity parameter in the variational problem C
! --> TFILE  character*39 name of file containing final single precision table C
!------------------------------------------------------------------------------C
      SUBROUTINE IROCTA(XD,NIT,NITSOR,NCY,A,TFILE)
      CHARACTER*39 TFILE
      REAL A,X
      PARAMETER(NGEN=6,NRAT1=1,NRAT2=3,N=NRAT2*2**NGEN,NP=N+1)
      COMMON/QPAN/X(3,-1:NP,-1:NP)
      DIMENSION XD(0:*),NG(0:10),NXG(0:10),AODDG(0:10)
     *,DG(0:10),LX(0:10+1),LF(0:10),LD(0:10)
      ORC=1.2
      NG(0)=NRAT2
      DO IGEN=1,NGEN
       IGENM=IGEN-1
       NG(IGEN)=2*NG(IGENM)
      ENDDO
      LX(0)=0
      DO IGEN=0,NGEN
       DG(IGEN)=1./NG(IGEN)
       AODDG(IGEN)=A/(8.*DG(IGEN)**2)
       NXG(IGEN)=3*(NG(IGEN)+3)**2
       LF(IGEN)=LX(IGEN)+NXG(IGEN)
       LD(IGEN)=LF(IGEN)+NXG(IGEN)
       LX(IGEN+1)=LD(IGEN)+NXG(IGEN)
      ENDDO
      LXI=LX(NGEN)
      NGI=NG(NGEN)
      NGA=(NRAT1*NGI)/NRAT2
      LFI=LF(NGEN)
      LDI=LD(NGEN)
      AODDGI=AODDG(NGEN)
!  INITIALIZE FINE GRID ESTIMATE USING CONFORMAL SOLUTION:
      CALL INSUBP(XD(LXI),NGA,NGI)
!misha
!      goto 123
!misha

!  PERFORM NIT PRELIMINARY ITERATIONS OF RELAXATION (SMOOTHING) ON FINE GRID:
      DO IT=1,NIT
       CALL RELSUBP(XD(LXI),XD(LFI),XD(LDI),NGA,NGI,ORC,AODDGI,ACC,0)
      ENDDO
      PRINT'('' LEVEL '',I1,'' ACC= '',E12.6)',NGEN,ACC

!  PERFORM NCY CYCLES OF NONLINEAR MULTIGRID WITH NGEN REFINEMENTS:
      DO ICY=1,NCY

!  GO FROM FINER TO COARSER GRID REPRESENTATIONS:
       DO IGEN=NGEN-1,0,-1
        IGENP=IGEN+1
        LXIP=LX(IGENP)
        LFIP=LF(IGENP)
        LDIP=LD(IGENP)
        NGIP=NG(IGENP)
        NGAP=(NRAT1*NGIP)/NRAT2
        AODDGIP=AODDG(IGENP)
        LXI=LX(IGEN)
        LFI=LF(IGEN)
        LDI=LD(IGEN)
        NGI=NG(IGEN)
        NGA=(NRAT1*NGI)/NRAT2
        NXGI=NXG(IGEN)
        AODDGI=AODDG(IGEN)

!  EVALUATE FIELD OF RESIDUALS (D) AT FINER GRID:
        CALL RELSUBP(XD(LXIP),XD(LFIP),XD(LDIP),NGAP,NGIP,1.0
     *  ,AODDGIP,ACC,1)
        CALL ZERV(XD(LFI),NXGI)

!  COMBINE FINER GRID RESIDUALS (D) INTO A COARSER GRID REPRESENTATION (F):
        CALL ZERV(XD(LFI),NXGI)
        CALL ADDX2X(XD(LDIP),XD(LFI),NGA,NGI)

!  SET COARSER GRID ESTIMATE (X) TO BE CONSISTENT WITH FINER GRID ESTIMATE (X):
        CALL SETX2X(XD(LXIP),XD(LXI),NGA,NGI)

!  SET COARSER GRID FORCING (F) TO MISMATCH BETWEEN FINER AND COARSER RESIDUALS:
        CALL RELSUBP(XD(LXI),XD(LFI),XD(LDI),NGA,NGI,1.0,AODDGI,ACC,2)

!  PERFORM NIT RELAXATION ITERATIONS WITH THIS FORCING (F) AT COARSER GRID:
        DO IT=1,NIT
         CALL RELSUBP(XD(LXI),XD(LFI),XD(LDI),NGA,NGI,ORC,AODDGI,ACC,0)
        ENDDO
        PRINT'('' LEVEL '',I1,'' ACC= '',E12.6)',IGEN,ACC
       ENDDO

!  GO FROM COARSER TO FINER GRID REPRESENTATIONS:
       DO IGEN=0,NGEN-1
        IGENP=IGEN+1
        LXIP=LX(IGENP)
        LFIP=LF(IGENP)
        LDIP=LD(IGENP)
        NGIP=NG(IGENP)
        NGAP=(NRAT1*NGIP)/NRAT2
        AODDGIP=AODDG(IGENP)
        LXI=LX(IGEN)
        LFI=LF(IGEN)
        LDI=LD(IGEN)
        NGI=NG(IGEN)
        NGA=(NRAT1*NGI)/NRAT2
        NXGI=NXG(IGEN)
        AODDGI=AODDG(IGEN)

!  FIND DIFFERENCE (ON COARSER GRID) BETWEEN COARSER AND FINER GRID ESTIMATES X:
        CALL DIFXX2(XD(LXI),XD(LXIP),NGA,NGI)

!  INTERPOLATE AND ADD THIS DIFFERENCE TO FINER GRID ESTIMATE (X):
        CALL ADDXX2(XD(LXI),XD(LXIP),NGA,NGI)

!  PERFORM NIT RELAXATION ITERATIONS AT FINER GRID TO SMOOTH INTERPOLANTS:
        DO IT=1,NIT
         CALL RELSUBP(XD(LXIP),XD(LFIP),XD(LDIP),NGAP,NGIP,ORC
     *   ,AODDGIP,ACC,0)
        ENDDO
        PRINT'('' LEVEL '',I1,'' ACC= '',E12.6)',IGENP,ACC

       ENDDO
       PRINT'('' CYCLE '',I4,''  COMPLETE'')',ICY
      ENDDO

!  MULTIGRID CYCLES COMPLETE. NOW POLISH RESULT WITH EXTRA S.O.R. ITERATIONS:
      LXI=LX(NGEN)
      LFI=LF(NGEN)
      LDI=LD(NGEN)
      NGI=NG(NGEN)
      NGA=(NRAT1*NGI)/NRAT2
      NXGI=NXG(NGEN)
      AODDGI=AODDG(NGEN)
      DO IT=1,NITSOR
       CALL RELSUBP(XD(LXI),XD(LFI),XD(LDI),NGA,NGI,ORC
     * ,AODDGI,ACC,0)
       PRINT'('' LEVEL '',I1,'' ACC= '',E12.6)',NGEN,ACC
      ENDDO

!  COPY FINE GRID SOLUTION TO SINGLE PRECISION TABLE AND STORE FOR POSTERITY:
!misha
 123  CALL QTABLE(XD(LXI),X,NGA,NGI)
!     CALL QTABLE(XD(LXI),X,NGA,NGI)
!misha
!jim      OPEN(UNIT=9,FILE=TFILE,STATUS='UNKNOWN',FORM='BINARY')
!misha
      OPEN(UNIT=9,FILE=TFILE,STATUS='UNKNOWN',FORM='unformatted')
!misha
      WRITE(9)X
      CLOSE(UNIT=9)
!test
!      do i=-1,np
!      do j=-1,np
!        print *,x(1,i,j), x(2,i,j), x(3,i,j)
!      end do
!      end do
!test

!  TASK COMPLETE; NOW GRAB ANOTHER CUP OF COFFEE.
      RETURN
      END

!------------------------------------------------------------------------------C
!   R.J.Purser, NCEP, Washington D.C.   1996                                   C
!                   SUBROUTINE QTABLE                                          C
!  Symmetrize and copy results at finest grid of multigrid solution to a       C
!  single precision table of reference grid positions. The reference grid      C
!  covers one quadrant of the square containing the standard octagon,          C
!  plus an extra line of data points all round. Location (0,0) is the center   C
!  of the square containing the octogon, (N,N) is the corner of the quadrant   C
!  of the square, (N,NA) and (NA,N) are the two corners of the octagon in the  C
!  quadrant.                                                                   C
!                                                                              C
! <-> XD    double precision array, symmetrized on output, of table data       C
! <-- X     single precision array of symmetrized table of cartesianlocations  C
! --> NA    number of grid steps between octagon corner and center line        C
! --> N     number of grid steps from center to side of square                 C
!------------------------------------------------------------------------------C
      SUBROUTINE QTABLE(XD,X,NA,N)
      DIMENSION XD(3,-1:N+1,-1:N+1),X(3,-1:N+1,-1:N+1),LOFIE(3)
     *,JOFIF(3)
      DATA LOFIE/1,1,-1/,JOFIF/2,1,3/
!  SYMMETRIZE ACROSS E-EDGE
      NP=N+1
      NA4=NA+N
      DO I=1,3
       L=LOFIE(I)
       DO IY=NA+1,NP
        JX=NA4-IY
        DO IX=JX+2,NP
         JY=NA4-IX
         XD(I,IX,IY)=XD(I,JX,JY)*L
        ENDDO
       ENDDO
      ENDDO
!  COPY TO X WITH SYMMETRIZATION ACROSS THE DIAGONAL, X=Y:
      DO I=1,3
       J=JOFIF(I)
       DO IY=-1,NP
        DO IX=-1,IY
         X(I,IX,IY)=XD(J,IY,IX)
         X(J,IY,IX)=XD(J,IY,IX)
        ENDDO
       ENDDO
      ENDDO
      RETURN
      END

!------------------------------------------------------------------------------C
!   R.J.Purser, NCEP, Washington D.C.   1996                                   C
!                   SUBROUTINE INSUBP                                          C
!  Initialize octagon-quadrant cartesian locations to those given by the       C
!  conformal mapping solution                                                  C
!                                                                              C
! <-- X     double precision array of locations of ocatgon quadrant & margin   C
! --> MA    steps between octagon corner and center line                       C
! --> M1    steps between edge and center line                                 C
!------------------------------------------------------------------------------C
      SUBROUTINE INSUBP(X,MA,M1)
!   INITIALIZATION OF GRID VALUES USING CONFORMAL SOLUTION:
      COMPLEX Z,W
      DIMENSION X(3,-1:M1+1,-1:M1+1)
      D=1./M1
      DO IY=0,M1
       YA=IY*D
       DO IX=0,MIN(M1,M1+MA-IY)
        XA=IX*D
        Z=XA+YA*(0,1)
!mish
!      Print *, "BEFORE   ix=",ix
!      Print *, "BEFORE   z=",z
!      Print *, "BEFORE   xa=",xa
!      Print *, "BEFORE   ya=",ya
!mish
        CALL TOCT(Z,W)
!mish
!      Print *, "AFTER   ix =",ix
!mish
        XW=REAL(W)
        YW=AIMAG(W)
        XX=XW**2+YW**2
        S=1./(1.+XX)
        X(1,IX,IY)=2.*XW*S
        X(2,IX,IY)=2.*YW*S
        X(3,IX,IY)=(1.-XX)*S
       ENDDO
      ENDDO
      CALL SYMMET(X(1,-1,-1),MA,M1)
!mish
!      Print *, "AFTER"
!      stop
!mish
      RETURN
      END

!------------------------------------------------------------------------------C
!   R.J.Purser, National Meteorological Center, Washington D.C.  1995          C
!                   SUBROUTINE TOCT                                            C
!   Transform from complez-Z in the standard unit-octagon to complex-W in the  C
!   unit-circle                                                                C
!------------------------------------------------------------------------------C
      SUBROUTINE TOCT(Z,W)
      COMPLEX CI,CIR,CIQ,Z,W,ZT
      LOGICAL KX,KY,KXY,KX1,KXY1,KS
      COMMON/CSTOCT1/A1(40),A2(40),B1(40),B2(40),A,B,CX,CY
     *,SSR,N1,ST,STI,CI,CIR,CIQ
      X=REAL(Z)
      Y=AIMAG(Z)
      if(z.eq.0.) then
         W=X+Y*(0,1)
         return
      end if
      if(x.eq.1.) x=0.9999
      if(y.eq.1.) y=0.9999
      KX=X.LT.0.
      IF(KX)X=-X
      KY=Y.LT.0.
      IF(KY)Y=-Y
      KXY=Y.GT.X
      IF(KXY)THEN
       T=X
       X=Y
       Y=T
      ENDIF
      KX1=X.GT.1.
      KXY1=X+Y.GT.1.+A
      IF(KX1)THEN
       X=2.-X
      ELSEIF(KXY1)THEN
       T=X
       X=1.+A-Y
       Y=1.+A-T
      ENDIF
      DD1=X**2+Y**2
      DD2=(1.-X)**2+(A-Y)**2

      ZT=X+Y*(0,1)
      KS=DD1.LT.SSR*DD2
!misha
!        print *, "come to break1,  ks=", ks
!        print *, "              ,  x= ",x
!        print *, "              ,  y= ",y
!misha
      IF(KS)THEN
       ZT=ZT**4
       CALL TAY(ZT,A1,N1,W)
       W=CIQ*(W/CI)**.25
      ELSE
       ZT=A+CI*ZT-CI
       ZT=ZT**ST
       CALL TAY(ZT,A2,N1,W)
       W=W+B
       W=(W+CI)/(CI-W)
      ENDIF
      IF(KX1.OR.KXY1)W=W/CABS(W)**2
      X=REAL(W)
      Y=AIMAG(W)
      IF(KXY)THEN
       T=X
       X=Y
       Y=T
      ENDIF
      IF(KY)Y=-Y
      IF(KX)X=-X
      W=X+Y*(0,1)
!misha
!        print *, "              ,  x= ",x
!        print *, "              ,  y= ",y
!misha
      RETURN
      END

!------------------------------------------------------------------------------C
!   R.J.Purser, National Meteorological Center, Washington D.C.  1994          C
!                   SUBROUTINE  TAY                                            C
!  Evaluate the complex function W of Z whose real                             C
!  Taylor series coefficients are RA.                                          C
!                                                                              C
!  --> Z    function argument (complex)                                        C
!  --> RA   Taylor coefficients (real)                                         C
!  --> N    number of coeffients (starting with the linear term)               C
!  <-- W    Taylor-series approximation of the function (complex)              C
!------------------------------------------------------------------------------C
      SUBROUTINE TAY(Z,RA,N,W)
      COMPLEX Z,W
      DIMENSION RA(*)
      W=0.
      DO I=N,1,-1
       W=(W+RA(I))*Z
      ENDDO
      RETURN
      END

!------------------------------------------------------------------------------C
!   R.J.Purser, National Meteorological Center, Washington D.C.  1994          C
!                   SUBROUTINE CINVRT                                          C
!  Compute the Taylor series coefficients Z for the functional-inverse of      C
!  the function whose Taylor series coefficients are W, for the case where     C
!  the constant coefficients of both Z and W are 0                             C
!                                                                              C
!  --> W    Taylor coefficients of original function (starting with linear)    C
!  <-- Z    Taylor coefficients of inverse function                            C
!  --> M    number of Taylor series coefficients computed                      C
!  --- WORK workspace array consisting of at least 3*M elements                C
!------------------------------------------------------------------------------C
      SUBROUTINE CINVRT(W,Z,M,WORK)
      DIMENSION W(*),Z(*),WORK(M,3)
      W1=W(1)
      DO I=1,M
       WORK(I,1)=W(I)
       WORK(I,2)=0.
      ENDDO
      WORK(1,2)=1.
      DO J=1,M
       ZJ=WORK(J,2)/WORK(J,1)
       Z(J)=ZJ
       DO I=J,M
        WORK(I,2)=WORK(I,2)-ZJ*WORK(I,1)
       ENDDO
       CALL CONV(WORK(1,1),W,WORK(1,3),M)
       DO I=1,M
        WORK(I,1)=WORK(I,3)
       ENDDO
      ENDDO
      RETURN
      END

!------------------------------------------------------------------------------C
!   R.J.Purser, National Meteorological Center, Washington D.C.  1994          C
!                   SUBROUTINE CONV                                            C
!   Convolve double-precision series A with B to form C, up to M terms         C
!   starting with element 1                                                    C
!                                                                              C
! A,B   --> inputs (convolution factors)                                       C
! C     <-- output (convolution product)                                       C
! M     --> number of elements of A, B, C                                      C
!------------------------------------------------------------------------------C
      SUBROUTINE CONV(A,B,C,M)
      DIMENSION A(*),B(*),C(*)
      DO K=1,M
       C(K)=0.
      ENDDO
      DO I=1,M
       DO J=1,M-I
        K=I+J
        C(K)=C(K)+A(I)*B(J)
       ENDDO
      ENDDO
      RETURN
      END


