Program Init_round
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
!     ******************************************************************
!     *                                                                *
!     *  Jim Purser's Rounded Cube:                                    *
!     *                                                                *
!     *            Create Initial Table for Transformation             *
!     *                                                                *
!     ******************************************************************
      PARAMETER(N=64,NM=N-1,NP=N+1)
     
ccccc      PARAMETER(N=128,NM=N-1,NP=N+1) !!!!!!!!!!DRAGAN
     
      COMMON/QPAN/X(3,-1:NP,-1:NP)
		character*36 filename
c
C  NIT IS THE NUMBER OF RELAXATION ITERATIONS ( ABOUT 100 ?)
C  ORC IS THE OVER-RELAXATION COEFFICIENT (ABOUT 1.3 ?)
C  A   IS THE PARAMETER OF TERM-2 IN THE EULER-LAGRANGE PROBLEM ( 10 ?)

      nit=100
cccc      nit=400 !!!!!!!DRAGAN
      orc=1.3
!      a=0.               ! Conformal
      a=10.              ! Radi
!       a=20.

C  INITIALIZE THE GRID BY SOLVING THE EULER-LAGRANGE PROBLEM:

      CALL IRCUBE(NIT,ORC,A)

c
c Write down the table
c
      filename='round.dat'
      Open(unit=10,file=filename,form='unformatted')
      Write(10) X
      Close(10)
c-----------------------------------------------------------------------
                             Stop
			     End
c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
C------------------------------------------------------------------------------C
C   R.J.Purser, National Centers for Environmental Prediction, Washington D.C. C
C   wd23jp@sun1.wwb.noaa.gov					      1996     C
C		    SUBROUTINE IRCUBE					       C
C    Initialize the table for the "rounded cube" grid transformations	       C
C									       C
C --> NIT:  Number of relaxation iterations (more the better - at least 100)   C
C --> SORC: Single precision over-relaxation coefficient (sorc=1.3 works)      C
C --> A:    Relative weight given to 2nd term in Euler-Lagrange functional     C
C	    (A=0 reverts to conformal case, A=10. typical of rounded cube with C
C	    a more even density of grid points) 			       C
C------------------------------------------------------------------------------C
C									       C
C									       C
C   ORIENTATION CONVENTION WITHIN ONE PANEL:				       C
C									       C
C		  F   B   E						       C
C		      | 	     ^					       C
C		  C - O - A	   y |					       C
C		      | 	     |					       C
C		  G   D   H		---> x				       C
C									       C
C   PANEL NUMBERING CONVENTION AND RELATIVE ORIENTATIONS:		       C
C									       C
C    ---								       C
C   | 5 |	       (1,0,0) = center of panel 1			       C
C    --- --- --- ---							       C
C   | 1 | 2 | 3 | 4 |  (0,1,0) = center of panel 2			       C
C    --- --- --- ---							       C
C   | 6 |	       (0,0,1) = center of panel 5			       C
C    ---								       C
C									       C
C   The parameter N is half the number of grid-spaces along one side of a      C
C   basic panel of the "mapping table" used to construct future transformationsC
C									       C
C   Symmetry is exploited so that only one quadrant of a map panel's table     C
C   is explicitly stored (plus one-row or one-column of overlap all round      C
C   this quadrant). The data in this table "X" comprise the 3-vector positions C
C   (components ordered by first index) of each mesh-node of the representativeC
C   quadrant. Thus, the value X(I,IX,IY) is the Ith component of the 3-vector  C
C   location of the point on the unit-sphere that corresponds to the table's   C
C   grid point (IX,IY). The center of the standard map panel is denoted by     C
C   (IX,IY) = (0,0); the right-edge's midpoint by (N,0); the upper-right cornerC
C   of the representative panel by (N,N). The actual table has IX and IY       C
C   in [-1,N+1], the overlap values, set by symmetry, allow		       C
C   quadratic-interpolation within the proper quadrant (IX and IY in [0,N])    C
C   to be performed.							       C
C									       C
C    IMPORTANT: Make sure the parameter N used in IRCUBE is exactly	       C
C    consistent with same parameter used in subroutine XMTOXC		       C
C    (since the table is passed through common block /QPAN/)		       C
C									       C
C------------------------------------------------------------------------------C
      SUBROUTINE IRCUBE(NIT,SORC,A)
      PARAMETER(N=64,NM=N-1,NP=N+1)
ccccccc      PARAMETER(N=128,NM=N-1,NP=N+1) !!!!!DRAGAN
      
      REAL*8 XD,AODD,ORC
      COMMON/QPAN/X(3,-1:NP,-1:NP)
      DIMENSION XD(3,0:N,0:N)
      DIMENSION JOFIA(3),JOFIB(3),LOFIC(3),LOFID(3),JOFIF(3)
      DATA JOFIA/3,2,1/,JOFIB/1,3,2/,LOFIC/-1,1,1/,LOFID/1,-1,1/
     *,JOFIF/2,1,3/
      ORC=SORC
      D=1./N
      DO IX=0,N
      DO IY=0,N
       X(1,IX,IY)=IX*D
       X(2,IX,IY)=IY*D
       X(3,IX,IY)=1.
      S=0.
      DO I=1,3
       S=S+X(I,IX,IY)**2
      ENDDO
      SI=1./SQRT(S)
      DO I=1,3
       X(I,IX,IY)=X(I,IX,IY)*SI
      ENDDO
      ENDDO
      ENDDO
      DO IX=0,N
      DO IY=0,N
      DO I=1,3
       XD(I,IX,IY)=X(I,IX,IY)
      ENDDO
      ENDDO
      ENDDO

C  PERFORM RELAXATION ITERATIONS TO SOLVE THE EULER-LAGRANGE PROBLEM:
      AODD=A/(8*D*D)
      DO IT=1,NIT
       PRINT'('' IT='',I4)',IT
       DO IX=1,NM
	IXM=IX-1
	IXP=IX+1
	CALL RELAXD(XD(1,IX,0),XD(1,IXP,0),XD(1,IX,1),XD(1,IXM,0)
     *	,XD(1,IXP,1),XD(1,IXM,1),ORC,AODD)
       ENDDO
       DO IY=1,NM
	IYM=IY-1
	IYP=IY+1
	CALL RELAXC(XD(1,0,IY),XD(1,1,IY),XD(1,0,IYP),XD(1,0,IYM)
     *	,XD(1,1,IYP),XD(1,1,IYM),ORC,AODD)
	DO IX=1,NM
	 IXM=IX-1
	 IXP=IX+1
	 CALL RELAX(XD(1,IX,IY),XD(1,IXP,IY),XD(1,IX,IYP),XD(1,IXM,IY)
     *	 ,XD(1,IX,IYM),XD(1,IXP,IYP),XD(1,IXM,IYP),XD(1,IXM,IYM)
     *	 ,XD(1,IXP,IYM),ORC,AODD)
	ENDDO
	CALL RELAXA(XD(1,N,IY),XD(1,N,IYP),XD(1,NM,IY),XD(1,N,IYM)
     *	,XD(1,NM,IYP),XD(1,NM,IYM),ORC,AODD)
       ENDDO
       DO IX=1,NM
	IXM=IX-1
	IXP=IX+1
	CALL RELAXB(XD(1,IX,N),XD(1,IXP,N),XD(1,IXM,N)
     *	,XD(1,IX,NM),XD(1,IXM,NM),XD(1,IXP,NM),ORC,AODD)
       ENDDO
      ENDDO

C  SYMMETRIZE ACROSS DIAGONAL:
      DO I=1,3
       J=JOFIF(I)
       DO IX=0,N
       DO IY=0,N
	X(I,IX,IY)=.5*(XD(I,IX,IY)+XD(J,IY,IX))
       ENDDO
       ENDDO
      ENDDO
C  SYMMETRIZE ACROSS MID-LINES AND EDGES:
      DO I=1,3
       J=JOFIA(I)
       L=LOFIC(I)
       DO IY=0,N
	X(I,-1,IY)=X(I,1,IY)*L
	X(I,NP,IY)=X(J,NM,IY)
       ENDDO
      ENDDO
      DO I=1,3
       J=JOFIB(I)
       L=LOFID(I)
       DO IX=-1,NP
	X(I,IX,-1)=X(I,IX,1)*L
	X(I,IX,NP)=X(J,IX,NM)
       ENDDO
      ENDDO
      RETURN
      END

      SUBROUTINE RELAXA(XO,XB,XC,XD,XF,XG,ORC,AODD)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION XO(3),XA(3),XB(3),XC(3),XD(3),XE(3),XF(3),XG(3),XH(3)
     *,JOFI(3)
      DATA JOFI/3,2,1/
      DO I=1,3
       J=JOFI(I)
       XA(I)=XC(J)
       XE(I)=XF(J)
       XH(I)=XG(J)
      ENDDO
      CALL RELAX(XO,XA,XB,XC,XD,XE,XF,XG,XH,ORC,AODD)
      RETURN
      END

      SUBROUTINE RELAXB(XO,XA,XC,XD,XG,XH,ORC,AODD)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION XO(3),XA(3),XB(3),XC(3),XD(3),XE(3),XF(3),XG(3),XH(3)
     *,JOFI(3)
      DATA JOFI/1,3,2/
      DO I=1,3
       J=JOFI(I)
       XB(I)=XD(J)
       XE(I)=XH(J)
       XF(I)=XG(J)
      ENDDO
      CALL RELAX(XO,XA,XB,XC,XD,XE,XF,XG,XH,ORC,AODD)
      RETURN
      END

      SUBROUTINE RELAXC(XO,XA,XB,XD,XE,XH,ORC,AODD)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION XO(3),XA(3),XB(3),XC(3),XD(3),XE(3),XF(3),XG(3),XH(3)
     *,LOFI(3)
      DATA LOFI/-1,1,1/
      DO I=1,3
       L=LOFI(I)
       XC(I)=XA(I)*L
       XF(I)=XE(I)*L
       XG(I)=XH(I)*L
      ENDDO
      CALL RELAX(XO,XA,XB,XC,XD,XE,XF,XG,XH,ORC,AODD)
      RETURN
      END

      SUBROUTINE RELAXD(XO,XA,XB,XC,XE,XF,ORC,AODD)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION XO(3),XA(3),XB(3),XC(3),XD(3),XE(3),XF(3),XG(3),XH(3)
     *,LOFI(3)
      DATA LOFI/1,-1,1/
      DO I=1,3
       L=LOFI(I)
       XD(I)=XB(I)*L
       XG(I)=XF(I)*L
       XH(I)=XE(I)*L
      ENDDO
      CALL RELAX(XO,XA,XB,XC,XD,XE,XF,XG,XH,ORC,AODD)
      RETURN
      END

      SUBROUTINE RELAX(X,XA,XB,XC,XD,XE,XF,XG,XH,ORC,AODD)
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER (DELTA=.01D0,DELTAI=1.D0/DELTA)
      DIMENSION 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)
      SX=0.
      SY=0.
      DO I=1,3
       U(I)=(XA(I)-XC(I))*.5
       V(I)=(XB(I)-XD(I))*.5
       SX=SX+X(I)*U(I)
       SY=SY+X(I)*V(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 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*.5
      DYY=(DYY-DY)*DELTAI
      DETI=ORC/(DXX*DYY-DXY*DXY)
      SX=(DYY*DX-DXY*DY)*DETI
      SY=(DXX*DY-DXY*DX)*DETI
      S=0.
      DO I=1,3
       X(I)=X(I)-U(I)*SX-V(I)*SY
       S=S+X(I)**2
      ENDDO
      S=1./SQRT(S)
      DO I=1,3
       X(I)=X(I)*S
      ENDDO
      RETURN
      END

      SUBROUTINE REL(X,XA,XB,XC,XD,XE,XF,XG,XH,U,V,AODD,DX,DY)
C
C	    XF	XB  XE
C	      QF  QE
C	    XC	X   XA
C	      QG  QH
C	    XG	XD  XH
C
      IMPLICIT REAL*8(A-H,O-Z)
      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