!
!                                      ***************************************
!                                      *      module cugen                   *
!                                      *      R. J. Purser                   *
!                                      *      NOAA/NCEP/EMC April 2014       *
!                                      *      jim.purser@noaa.gov            *
!                                      ***************************************
!
! A suite of routines for generating the coefficients needed to define the
! equal-area cubic grid by the method of "map migration" defined in the
! NOAA/NCEP office note 467 (Purser and Rancic) and its supplement (a future
! office note).
!
!
! DIRECT DEPENDENCIES
! Libraries[their Modules]: pbint[pbint], pcomb[pcomb], pbend[pbend],
!                           pmat[pmat, peuc], pcoco[pcoco]
! Additional Modules      : cutable, pietc, pkind
!
!=============================================================================
module cugen
!=============================================================================
use pkind, only: dp
implicit none
interface cucogen;        module procedure cucogen;              end interface
interface migrate;        module procedure migrate;              end interface
interface outerits;       module procedure outerits;             end interface
interface innerits;       module procedure innerits;             end interface
interface curead;         module procedure curead;               end interface
interface cuwrite;        module procedure cuwrite;              end interface
interface init_map_table; module procedure init_map_table;       end interface
interface fini_map_table; module procedure fini_map_table;       end interface
interface bigaints;       module procedure bigaints;             end interface
interface senslap;        module procedure senslap;              end interface
interface plapchi;        module procedure plapchi;              end interface
interface fgradchi;       module procedure fgradchi;             end interface
interface pgradchi;       module procedure pgradchi;             end interface
interface fmetric;        module procedure fmetric;              end interface
interface pmetric;        module procedure pmetric;              end interface
interface xmtoxc;         module procedure xmtoxc,xmtoxc_s;      end interface
interface xctoxm;         module procedure xctoxm,xctoxm_s;      end interface
interface image;          module procedure image;                end interface
interface timage;         module procedure timage;               end interface
interface himage;         module procedure himage;               end interface
interface imagev;         module procedure imagev;               end interface
interface rtos;           module procedure rtos;                 end interface
interface stor;           module procedure stor;                 end interface
interface replicate_cu;   module procedure replicate_cu;         end interface
interface replicate_oh;   module procedure replicate_oh;         end interface
interface setgroup_cu;    module procedure setgroup_cu;          end interface

contains

!=============================================================================
subroutine cucogen(norh,n,ntau,mtau,nouterit,ninnerit)!              [cucogen]
!=============================================================================
! Cubic mapping coefficients generation method employing "migration".
! norh: half the nominal order of numerics used in both generation and
!       subsequent interpolations. Suggested value =3, (must be at least 1)
! n:    Size of the construction grid spanning the fundamental triangle
!       Suggested value =40 until multigrid allows this to be increased.
! ntau: Even number of "time" intervals in the full migration, which is done 2
!       steps at a time using the midpoint method.
!       Suggested value = 10. Better results might be obtained with 20.
! mtau: Actual even number of steps carried out. Obviously == ntau unless
!       during code development and improvement, one wants to look at the
!       results at an earlier stage of the migration. This parameter is
!       eventually destined for redundancy and removal!
! nouterit,ninnerit: number of allowed outer and inner iterations in the
!       modified SOR method used to solve the Poisson problem. 
!       Suggested values about 100 and 400? (More are needed when n is large.)
! 
! When successful, the constructed coefficients, both of the asymptotic 
! approximate solution near the vertex, and the gridded values, are put
! into the data module, "cutable". When done with them, please deallocate
! the space via a call to subroutine fini_map_table.
!=============================================================================
use pkind,only: dp,dpc
use pietc,only: zero,half,one,pi
use pbend,only: dif,ddif,quad
use pcoco,only: set_cuco, cuztoc0
use peuc, only: triple_product
implicit none
integer,intent(in ):: norh,n,ntau,mtau,nouterit,ninnerit
!-----------------------------------------------------------------------------
integer,parameter                :: m=55,nq=2,ns=20,              &
                                    fa=4, &! Face symmetry (=4 for cube)
                                    v=3,  &! Vertex symmetry (=3 for cu)
                                    v0=(2*fa)/(fa-2),&!Symm for planar tiling=4
                                    ngr=(4*v*v0)/(v0-v),& ! Order of group=48
                                    nfa=ngr/(2*fa),&! number of faces =6
                                    nv=ngr/(2*v) ! number of vertices =8
real(dp),parameter               :: pi4=pi*4,pio4=pi/4,               &
                                    gammac=(one*v)/v0,    &! 3/4 for cube
                                    gamma=(one*(v0-v))/v0,&! 1/4 for cube
                                    ov=one/v,             &! 1/3 for cube
                                    piov0=pi/v0,&! 45-degrees for cube's face
                                    rb=.15,ds=one/ns,                  &
                                    aahelm=.1,&! Helmholtz coefficient
                                    urc=.7     ! Relaxation coefficient
real(dp),dimension(ns)           :: rats,arint,ajac0
real(dp),dimension(3,0:n,0:n)    :: xc0s,xca
real(dp),dimension(  0:n,0:n)    :: w,wwb,jac0,ljac0,dacda0s,div
real(dp),dimension(3,2,0:n,0:n)  :: dxcdz0s
real(dp),dimension(2,2,0:n,0:n)  :: gco0s,gcontra0s
real(dp),dimension(2,0:n,0:n)    :: ddacdaz0s,hcontra0s
real(dp),dimension(m)            :: co1,co2
real(dp),dimension(3)            :: v1,v2,v3
real(dp)                         :: dx,dxx,x,y,area,r,s,drds,wa,wb,kay0,kays0,&
                                    dtau,dtau2,tau,p,dlpdt,dpdt,dbart,xran
real(dp),dimension(3,3)          :: mat
real(dp),dimension(  norh)       :: d
real(dp),dimension(0:norh)       :: dd
real(dp),dimension(0:ntau )      :: taus,bigA,dbar,kay,kays,ps
complex(dpc)                     :: z
complex(dpc),dimension(3,0:n,0:n):: dxc0s
integer                          :: i,itau,itau2,j,k,is,ix,iy,jx,jy,jpan
integer,dimension(1)             :: ii
!=============================================================================
if(norh<1)stop 'In cucogen; norh must be positive; suggest =3'
if(n<10)stop 'In cucogen; n should not be less than 10, suggest =40'
if(mod(ntau,2)/=0)stop 'In cucogen; ntau must be even'
if(ntau<=0)stop 'In cucogen; ntau must be positive; suggest =10 or 20'
if(mod(mtau,2)/=0)stop 'In cucogen; mtau must be even'
if(mtau<0)stop 'In cucogen; mtau cannot be negative'
if(mtau>ntau)stop 'In cucogen; mtau must not exceed ntau'

dx=half/n; dxx=dx*dx
call  dif(norh, d)    ! Create 1st derivative stencil, d
call ddif(norh, dd)   ! Create 2nd derivative stencil, dd
! Find the two arrays, co1 and co2, of conformal cubic expansion coefficients 
! by the "circles" method. Each array is of size m. co1 expands about the
! vertex; co2, although formally redundant for the cube, expands about the
! face-center and therefore provides very accurate results near the middle.
call set_cuco(m,co1,co2); kay0=2*co1(1)**ov/gammac; kays0=kay0**2

dtau=one/ntau;dtau2=dtau*2

! To compensate for repetition of representation of points along edges and
! at corners of the fundamental triangle when summing the contributions of
! its ngr replicates, down-weight these points, using w(:,:), in the
! appropriate proportion:
w=one
do i=1,n-1
   w(i,0)=half ! Base edge, minus end points
   w(n,i)=half ! median edge, minus end points
   w(i,i)=half ! diagonal edge, minus end points
enddo

w(0,0)=half/v ! 60-degree corner (the vertex)
w(n,0)=one/4 ! 90-degree corner (the face-edge midpoint)
w(n,n)=half/fa ! 45-degree corner (the face center)

! Zero out unused portion:
do i=0,n
   w(i,i+1:n)=zero
enddo

! Make a copy of w, except modulated by the blending weight, wb, associated
! with the grid-based contribution (as opposed to the asymptotic approximate
! contribution near the vertex, whose weight is the complement, wa=1-wb).
wwb=w
do iy=0,n
   y=iy*dx
   if(y>=rb)exit
   do ix=iy,n
      x=ix*dx
      r=sqrt(x**2+y**2)
      if(r>=rb)exit
      call rtos(nq,rb,r,wb,s)!<- get annular blending weight wb at this r
      wwb(ix,iy)=w(ix,iy)*wb ! Modulate w with blending weight wb near vertex
   enddo
enddo

! Set weighted integration weights for the radial asymptotic part of any
! integral. S is a reparametrization of r defined by r=rb*s**nq, which
! allows a uniform grid in s in (0,1] to resolve the disc r<rb especially
! well near the vertex, r=0.
arint(ns)=0
rats(ns)=rb
do is=1,ns ! ns is the number of radial intervals in s resolving disc r<rb
   s=is*ds ! effective radial location of annulus in the s coordinate...
   call stor(nq,rb,s,wb,r,drds)! <- ...and its equivalent map-radius, r.
   wa=one-wb ! Complement of annular blending weight wb
   arint(is)=piov0*r*wa*drds*ds ! <- weighted area of annular sector
   rats(is)=r                   ! <- "r at s"
   ajac0(is)=kays0*r**(gamma*2/gammac) ! <- asymptotic jacobian, conformal.
enddo

do iy=0,n
   y=iy*dx
   do ix=iy,n
      x=ix*dx
      z=cmplx(x,y,dpc)
      call cuztoc0(m,co1,co2,z,xc0s(:,ix,iy))
   enddo
enddo

call fmetric(norh,n,d,dd, xc0s, dxcdz0s,dacda0s,gco0s,gcontra0s, &
      ddacdaz0s,hcontra0s)
jac0=dacda0s; jac0(0,0)=1; ljac0=log(jac0); jac0(0,0)=0


jac0=0
do iy=0,n
   y=iy*dx
   do ix=iy,n
      x=ix*dx
      z=cmplx(x,y,dpc)
      if(ix==0 .and. iy==0)then
         call cuztoc0(m,co1,co2,z,xc0s(:,ix,iy))
      else
         call cuztoc0(m,co1,co2,z,xc0s(:,ix,iy),dxc0s(:,ix,iy))
         v1= real(dxc0s(:,ix,iy))
         v2=aimag(dxc0s(:,ix,iy))
         v3=xc0s(:,ix,iy)
         jac0(ix,iy)=triple_product(v1,v3,v2)
      endif
   enddo
enddo

do itau=0,ntau
   tau=itau*dtau; taus(itau)=tau
   ps(itau)=(one-tau*gamma)/gammac
enddo
taus(ntau)=one
ps  (ntau)=one
p=ps(0)
dlpdt=dpdt/p

do itau=0,ntau
   tau=taus(itau)
   call bigaints(norh,n,ns,kay0,gamma,tau,rats,arint,ajac0, wwb,jac0,     &
        biga(itau),dbar(itau),kay(itau),kays(itau))
enddo
xca=xc0s

call migrate(norh,n,ns,ntau,mtau,nouterit,ninnerit, &
      rb,gamma,kay0,aahelm,urc, d,dd, &
      dbar,biga,kay,ps,ljac0, xca)

call init_map_table(norh,n,rb,gamma,kay(mtau),ps(mtau),xca)
62 format(6(1x,e14.7))
64 format(i5,3(1x,e21.14))
end subroutine cucogen

!=============================================================================
subroutine migrate(norh,n,ns,ntau,mtau,nouterit,ninnerit, &
      rb,gamma,kay0,aahelm,urc, d,dd, &
      dbar,biga,kay,ps,ljac0, xca)
!=============================================================================
! Migrate the position vectors, xca, of a fine regular grid covering the
! fundamental region, for ntau/2 double time steps using the midpoint
! scheme. Denote the three time levels of each double interval by A, B, C.
! Step 1 computes the 3-vector tendendies at A and advances xc by a forward
! step to B. Step 2 computes the 3-vector tendencies at B and advances
! again, but this time from A to C. Rename C the new A.
! The vector tendencies come from the gradient of a velocity potential chi
! and chi is solved from a Poisson equation forced by divergence, D. However,
! rather than attempting to solve the Poisson directly, we use iterated
! Helmholtz equations with very small Helmholtz coefficient aahelm, and
! use a different forcing than just divergence, to ensure that the sequence
! of Helmholtz solutions converges rapidly to that of the intended Poisson.
! The details and motivation, described fully in the office note, relate to
! the fact that Poisson equations on the sphere are formally ill-posed (a
! solution is only possible when the forcing D integrates to zero; a solution
! obtained is even then not unique because an additive constant does not
! change the Laplacian. For a robust solution we need (1) to allow for an
! additive constant to D to allow any solution to be obtained; (2) a way to
! force the solution chi=0 at the vertex. This is the reason for the slightly
! complicated inner and outer iterations nested here at each step.
! Also, we recognize that the grid-based numerical solution is inadequate
! near the vertex, where an asymptotic estimate of each important field
! is available. Therefore, with the ring of radii rb/2,rb about the vertex,
! we blend the asymptotic estimates of chi, and V=grad(chi) with the 
! grid-derived versions, using the smooth blending function obtained from
! subr. wbend.
!==============================================================================
use pkind, only: dp
use pietc, only: half,one
use pbend, only: wbend
implicit none
integer,                      intent(in   ):: norh,n,ns,ntau,mtau, &
                                              nouterit,ninnerit
real(dp),                     intent(in   ):: rb,gamma,kay0,aahelm,urc
real(dp),dimension(  norh),   intent(in   ):: d
real(dp),dimension(0:norh),   intent(in   ):: dd
real(dp),dimension(0:ntau),   intent(in   ):: dbar,biga,kay,ps
real(dp),dimension(  0:n,0:n),intent(in   ):: ljac0
real(dp),dimension(3,0:n,0:n),intent(inout):: xca
!-----------------------------------------------------------------------------
integer                              :: ix,iy,iatau,ibtau,ictau,nb
real(dp)                             :: dtau,dtau2,atau,btau,ctau,         &
                                        dx,dxx,x,y,xx,yy,r,lkay0,gammac,   &
                                        dbart,bigat,kayt,p,bigr,vr,bigc,s, &
                                        wat,wbt,logr,theta,ctheta,stheta
real(dp),dimension(3)                :: xc
real(dp),dimension(    0:n,0:n)      :: dacda,div,chi,sens
real(dp),dimension(2,  0:n,0:n)      :: hcontra,ddacdaz
real(dp),dimension(3,2,0:n,0:n)      :: dxcdz
real(dp),dimension(2,2,0:n,0:n)      :: gco,gcontra
real(dp),dimension(3,0:n,0:n)        :: dchidxc,xcb
real(dp),allocatable,dimension(:,:)  :: wb,achi
real(dp),allocatable,dimension(:,:,:):: adchidxc
!=============================================================================
if(mod(mtau,2)/=0)stop'mtau must be even for a whole number of double-interval'
print'('' Just inside migrate:'')'

dx=half/n; dxx=dx**2
gammac=one-gamma
lkay0=log(kay0)
nb=floor(rb/dx)    ! <- Size of grid square just big enough to contain the
                   ! grid points where the solution is blended with the
                   ! estimates approximated asymptotically.
allocate( wb(0:nb,0:nb), achi(0:nb,0:nb) )
allocate( adchidxc(3,0:nb,0:nb) )
! Set up the mini-array of the value of the blending weight wb as it applies to
! the gridded part of the solution, near the vertex:
do iy=0,nb
   y=iy*dx; yy=y**2
   do ix=0,nb
      x=ix*dx; xx=x**2
      r=sqrt(xx + yy)
      call wbend((2*r-rb)/rb,wb(ix,iy))
   enddo
enddo

dtau=one/ntau; dtau2=dtau*2

bigc=0
chi=0 ! Consider replacing this initial estimate by the asymptotic formula?
do ictau=2,mtau,2 ! Loop 2 steps at a time; in each double interval, A, B, C
!                   refer to initial, midpoint, and final timelevel of interval
   ibtau=ictau-1
   iatau=ictau-2
   print'('' Migrate, iatau,ibtau,ictau='',5i3)',iatau,ibtau,ictau
!   read(*,*)
   atau=iatau*dtau ! "A" === Start of double-interval
   btau=ibtau*dtau ! "B" === midpoint of double-interval
   ctau=ictau*dtau ! "C" === end of double-interval (but becomes new "A")

! Apply step 1 of the midpoint scheme, advancing xc from A to B.
! Coefficients and arrays to do with asymptotics near the vertex:
   dbart=dbar(iatau)
   bigat=biga(iatau)
   kayt =kay (iatau)
   p    =ps  (iatau)
   achi(0,0)=0
   adchidxc(:,0,0)=0
   do iy=0,nb
      y=iy*dx
      yy=y**2
      do ix=max(1,iy),nb
         x=ix*dx
         xx=x**2
         r=sqrt(xx+yy)
         logr=log(r)
         bigr=gammac*kayt*r**p
         theta=atan2(y,x)/gammac; ctheta=cos(theta); stheta=sin(theta)
         achi(ix,iy)=(half*dbart-lkay0-(gamma/gammac)*(logr-one/p)) &
              *bigr**2/2
         vr=(half*dbart-lkay0-(gamma/gammac)*(logr-one/(p*2)))*bigr
         adchidxc(1,ix,iy)=vr*ctheta
         adchidxc(2,ix,iy)=vr*stheta
         adchidxc(3,ix,iy)=0
      enddo
   enddo
   call fmetric(norh,n,d,dd, xca, dxcdz,dacda,gco,gcontra, ddacdaz,hcontra)
   do iy=0,n
      do ix=max(1,iy),n
         div(ix,iy)=dbart-ljac0(ix,iy)
      enddo
   enddo
   call senslap(norh,n,d,dd,hcontra,gcontra,sens)
   call outerits(norh,n,nb,nouterit,ninnerit, rb,aahelm,urc,d,dd, &
        wb,achi,sens,div,hcontra,dxcdz,gcontra, bigc,chi)
   call fgradchi(norh,n,d,dxcdz,gcontra,chi,dchidxc)
   do iy=0,nb
      do ix=max(1,iy),nb
         wbt=wb(ix,iy);         wat=one-wbt
         dchidxc(:,ix,iy)=wat*adchidxc(:,ix,iy)+wbt*dchidxc(:,ix,iy)
      enddo
   enddo

! Increment position vectors, xc, stepping dtau from A to B
   do iy=0,n
      do ix=iy,n
         xc=xca(:,ix,iy)+dtau*dchidxc(:,ix,iy)
         s=sqrt(dot_product(xc,xc))
         xcb(:,ix,iy)=xc/s ! <- midpoint values of migrating xc
      enddo
   enddo

! Apply step 2 of the midpoint scheme, advancing xc from B to C (=== new A)

! Coefficients and arrays to do with asymptotics near the vertex:
   dbart=dbar(ibtau)
   bigat=biga(ibtau)
   kayt =kay (ibtau)
   p    =ps  (ibtau)
   achi(0,0)=0
   adchidxc(:,0,0)=0
   do iy=0,nb
      y=iy*dx
      yy=y**2
      do ix=max(1,iy),nb
         x=ix*dx
         xx=x**2
         r=sqrt(xx+yy)
         logr=log(r)
         bigr=gammac*kayt*r**p
         theta=atan2(y,x)/gammac; ctheta=cos(theta); stheta=sin(theta)
         achi(ix,iy)=(half*dbart-lkay0-(gamma/gammac)*(logr-one/p)) &
              *bigr**2/2
         vr=(half*dbart-lkay0-(gamma/gammac)*(logr-one/(p*2)))*bigr
         adchidxc(1,ix,iy)=vr*ctheta
         adchidxc(2,ix,iy)=vr*stheta
         adchidxc(3,ix,iy)=0
      enddo
   enddo
   call fmetric(norh,n,d,dd, xcb, dxcdz,dacda,gco,gcontra, &
      ddacdaz,hcontra)
   do iy=0,n
      do ix=max(1,iy),n
         div(ix,iy)=dbart-ljac0(ix,iy)
      enddo
   enddo
   call senslap(norh,n,d,dd,hcontra,gcontra,sens)
   call outerits(norh,n,nb,nouterit,ninnerit, rb,aahelm,urc,d,dd, &
        wb,achi,sens,div,hcontra,dxcdz,gcontra, bigc,chi)
   call fgradchi(norh,n,d,dxcdz,gcontra,chi,dchidxc)
! Blend velocity field with asymptotic values near vertex:
   do iy=0,nb
      do ix=max(1,iy),nb
         wbt=wb(ix,iy);         wat=one-wbt
         dchidxc(:,ix,iy)=wat*adchidxc(:,ix,iy)+wbt*dchidxc(:,ix,iy)
      enddo
   enddo

! Increment position vectors, xc, stepping dtau2 forward from A to C (new A)
   do iy=0,n
      do ix=iy,n
         xc=xca(:,ix,iy)+dtau2*dchidxc(:,ix,iy)
         s=sqrt(dot_product(xc,xc))
         xca(:,ix,iy)=xc/s ! <- normalized vectors at end of double interval
      enddo
   enddo

enddo
deallocate( wb,achi,adchidxc)

64 format(i3,3(1x,e21.14))
end subroutine migrate

!============================================================================
subroutine outerits(norh,n,nb,nouterit,ninnerit, rb,aahelm,urc,d,dd, &
      wb,achi,sens,div,hcontra,dxcdz,gcontra, bigc,chi)
!============================================================================
use pkind, only: dp
implicit none
integer,                        intent(in   ):: norh,n,nb,nouterit,ninnerit
real(dp),                       intent(in   ):: rb,aahelm,urc
real(dp),dimension(  norh),     intent(in   ):: d
real(dp),dimension(0:norh),     intent(in   ):: dd
real(dp),dimension(0:nb,0:nb),  intent(in   ):: wb,achi
real(dp),dimension(    0:n,0:n),intent(in   ):: sens,div
real(dp),dimension(2,  0:n,0:n),intent(in   ):: hcontra
real(dp),dimension(3,2,0:n,0:n),intent(in   ):: dxcdz
real(dp),dimension(2,2,0:n,0:n),intent(in   ):: gcontra
real(dp),                       intent(inout):: bigc
real(dp),dimension(    0:n,0:n),intent(inout):: chi
!-----------------------------------------------------------------------------
integer                        :: itm,ix,iy
real(dp)                       :: chi0,chi0old
real(dp),dimension(0:n,0:n)    :: ediv
real(dp),dimension(0:nouterit) :: bigcs
!=============================================================================
chi0old=0
bigcs(0)=bigc
do itm=1,nouterit
   do iy=0,n
      do ix=iy,n
         ediv(ix,iy)=div(ix,iy)-aahelm*(chi(ix,iy)-bigc)
      enddo
   enddo
   call innerits(norh,n,nb,itm,ninnerit, aahelm,urc,d,dd, &
        wb,achi,hcontra,dxcdz,gcontra,sens,ediv,chi0,chi)
   chi=chi-chi0
   bigc=bigc+chi0-chi0old
   chi0old=chi0
   bigcs(itm)=bigc
enddo
print'('' list the successive bigc:'')'
do itm=0,nouterit
   write(6,64)itm,bigcs(itm)
enddo
64 format(i5,3(1x,e13.6))

end subroutine outerits

!=============================================================================
subroutine innerits(norh,n,nb,itm,ninnerit, aahelm,urc,d,dd, &
      wb,achi,hcontra,dxcdz,gcontra,sens,div,chi0,chi)
!=============================================================================
! Perform the inner iterations for the helmholtz problem with helmholtz
! coefficient aahelm. The asymptotic part of the solution, only used where
! r<rb, is represented by chi0+achi (achi being the analytic part for the 
! POISSON equation's solution, standardized to achi(0,0)=0). The returned
! offset, chosen to optimize the fit between asymptotic and gridded solutions,
! is used by the calling routine, outerits, to nudge bigC, in order that the
! Helmholtz solution sequence can converge to the proper Poisson solution,
! even when the forcing, D, does not exactly integrate to zero over the
! domain numerically, as it is analytically required to do.
! Urc is the (under-)relaxation factor for the SOR iterations
!=============================================================================
use pkind, only: dp
use pietc, only: one
implicit none
integer,                        intent(in   ):: norh,n,nb,itm,ninnerit
real(dp),                       intent(in   ):: aahelm,urc
real(dp),dimension(  norh),     intent(in   ):: d
real(dp),dimension(0:norh),     intent(in   ):: dd
real(dp),dimension(0:nb,0:nb),  intent(in   ):: wb,achi
real(dp),dimension(2,  0:n,0:n),intent(in   ):: hcontra
real(dp),dimension(3,2,0:n,0:n),intent(in   ):: dxcdz
real(dp),dimension(2,2,0:n,0:n),intent(in   ):: gcontra
real(dp),dimension(    0:n,0:n),intent(in   ):: sens,div
real(dp),                       intent(  out):: chi0
real(dp),dimension(    0:n,0:n),intent(inout):: chi
!-----------------------------------------------------------------------------
real(dp)                     :: wwwci,wdchi,res,tchi,dchi0,lapchi,norm
real(dp),dimension(0:nb,0:nb):: wa,wawb
integer                      :: it,ix,iy
!=============================================================================
wwwci=0
do iy=0,nb
   do ix=iy,nb
      wa(ix,iy)=one-wb(ix,iy)
      wawb(ix,iy)=wa(ix,iy)*wb(ix,iy)
      wwwci=wwwci+wawb(ix,iy)*wb(ix,iy)
   enddo
enddo
chi0=0
chi(0,0)=chi0+achi(0,0)

do it=1,ninnerit
   norm=0
   do iy=0,n
      do ix=max(1,iy),n
         call plapchi(norh,n,ix,iy,d,dd, hcontra(:,ix,iy),gcontra(:,:,ix,iy),&
              chi,lapchi)
         res=-lapchi+aahelm*chi(ix,iy)+div(ix,iy)
         norm=norm+abs(res)
         tchi=chi(ix,iy)-urc*res/(sens(ix,iy)+aahelm)
         if(ix<=nb)then
! Blend with asymptotic solution, chi0+achi, with this relacation update, tchi:
            chi(ix,iy)=wa(ix,iy)*(chi0+achi(ix,iy))+wb(ix,iy)*tchi
         else
! .. or, if far from the vertex, just replace chi(ix,iy) with tchi:
            chi(ix,iy)=tchi
         endif
      enddo
   enddo
   norm=(norm*2)/((n+1)*(n+2)) ! <- norm per active grid point
   
!   print*, 'TESTE',norm,n
   if(mod(it,10)==0)print'("Inner iteration, it,norm=",i5,2x,e15.8)',it,norm
! Adjust asymptotic solution's offset, chi0, to better fit the gridded 
! solution in a weighted sense:
   wdchi=0
   do iy=1,nb
      do ix=iy,nb
         wdchi=wdchi+wawb(ix,iy)*(chi(ix,iy)-achi(ix,iy)-chi0)
      enddo
   enddo
   dchi0=wdchi/wwwci
   chi0=chi0+dchi0
! Use the adjustment of chi0 to smoothly adjust the gridded chi near the
! vertex
   do iy=0,nb
      do ix=iy,nb
         chi(ix,iy)=chi(ix,iy)+dchi0*wa(ix,iy)
      enddo
   enddo
enddo

64 format(i3,3(1x,e13.6))
end subroutine innerits

!=============================================================================
subroutine curead(filestem)
!=============================================================================
use pkind, only: dp
use pietc, only: half
use cutable, only: norh,n, rb,gamma,kay,p, xcs
implicit none
character(len=8),intent(in ):: filestem
!-----------------------------------------------------------------------------
character(len=12)                    :: filename
real(dp)                             :: rb_s,gamma_s,kay_s,p_s
real(dp),allocatable,dimension(:,:,:):: xcs_s
integer,parameter                    :: lu=10,mu=50
integer                              :: u, norh_s,n_s, iy
logical                              :: open
!=============================================================================
filename(1:8)=filestem; filename(9:12)='.dat'
print'('' full filename is '',a12)',filename

! Find the next available unit number within [lu:mu]:
do u=lu,mu
   inquire(u,opened=open)
   if(.not.open)exit
enddo
if(u>mu)stop 'In curead; unopened unit number within range is not found'

! Open the file, and write the critical data to it which defines the mapping:
open(u,file='cube_gen.dat',status='unknown',form='unformatted')
read(u)norh_s,n_s
read(u)rb_s,gamma_s,kay_s,p_s
if(norh_s<1 .or.norh_s>10)stop &
      'In curead; imported norh out of reasonable range'
if(n_s<10 .or. n_s>1000)stop &
      'In curead; imported n out of reasonable range'
if(rb_s<=0 .or. rb_s>=half)stop &
      'In curead; imported rb out of reasonable range'
if(gamma_s<=0 .or. gamma_s>half)stop &
      'In curead; imported gamma out of reasonable range'
allocate(xcs_s(3, 0:n_s, 0:n_s)); xcs_s=0
do iy=0,n_s
   read(u)xcs_s(:,iy:n_s,iy)
enddo
close(u)
call init_map_table(norh_s,n_s,rb_s,gamma_s,kay_s,p_s,xcs_s)
deallocate(xcs_s)

end subroutine curead

!=============================================================================
subroutine cuwrite(filestem)
!=============================================================================
use cutable, only: norh,n, rb,gamma,kay,p, xcs
implicit none
character(len=8),intent(in ):: filestem
!-----------------------------------------------------------------------------
character(len=12):: filename
integer,parameter:: lu=10,mu=50
integer          :: u, iy
logical          :: open
!=============================================================================
filename(1:8)=filestem; filename(9:12)='.dat'
print'('' full filename is '',a12)',filename
if(.not.allocated(xcs))stop 'In cuwrite; data xcs in cutable not allocated'

! Find the next available unit number within [lu:mu]:
do u=lu,mu
   inquire(u,opened=open)
   if(.not.open)exit
enddo
if(u>mu)stop 'In cuwrite; unopened unit number within range is not found'

! Open the file, and write the critical data to it which defines the mapping:
open(u,file='cube_gen.dat',status='unknown',form='unformatted')
write(u)norh,n
write(u)rb,gamma,kay,p
do iy=0,n
   write(u)xcs(:,iy:n,iy)
enddo
close(u)
end subroutine cuwrite

!=============================================================================
subroutine init_map_table(norh_s,n_s,rb_s,gamma_s,kay_s,p_s,xcs_s)
!=============================================================================
! Initial the map table according to the parameters and arrays provided
! in the argument list, and initialize the rotation matrices according
! to the mapping convention. That is: panels 1--4 wrap around the equator,
! panel 1 containing the 0 meridian, down its center, panel 2 the 90-degree (E)
! meridian and so on; panel 5 oriented like panel 1, except north of it, has
! the north pole at its center; likewise, panel 6 is oriented like panel 1,
! except south of it, and has the south pole at its center. Thus, unfolded,
! The panels form the lazy-T shape:
!
!  [5]
!  [1][2][3][4]
!  [6]
! 
! The standard user map coordinates, xm, in each panel run between -1 and +1.
! These xm coordinates, together with the panel index, ipan, are used to
! define a mapping to a smooth vector function, xc, and its derivative wrt xm,
! dxcdxm. 
!     The generic mapping, away from the vertex singularities, is found
! by smooth interpolation of the finely-gridded values stored in table xcs.
! For economy, this table is only defined in one fundamental triangle (1/48
! of the area) of the surface of the cube, and for a cube oriented with
! the vertex of the fundamental triangle at (0,0,1), the table's x-axis 
! pointing towards (1,0,0) and the table terminating in this direction at
! the edge midpoint. The tabulation map coordinates (x,y) differ from user 
! coordinates (for sound technical reasons) with (0,0) corresponding to
! the vertex and (1/2,0) corresponding to the square's edge midpoint. The
! third corner of the fundamental triangle, also the square's center, is
! representated by tabulation coordinates (1/2,1/2).
! The tabulation consists of the y<=x portion of a square array of the
! unit 3-vectors, xcs, array dimensions [0:n,0:n]. Interpolation is a 
! cartesian product of the smooth (norh*2+2)-point centered stencils
! obtained using subroutine getsplint.
!     The case where (x,y) are both small enough to be within a critical
! radius, rb, of the vertex (0,0) leads to a different evaluation of xc
! and its derivatives; in this case we blend the grid-interpolation result
! described above with the appropriate asymptotic approximation valid near
! vertex. The blending occurs over the annular region, rb/2<r<rb, such that
! the weight wb given to the grid=derived components rises in this radial
! interval smoothly from 0 to 1 in accordance with the "Whittaker" function
! blending defined in subr. wbend. The various coefficients by which the 
! asymptotic approximate solution is expanded are initialized in this
! subroutine from the values passed in through the argument list.
!     The actual evaluation of xc and its jacobian matrix wrt xm is done in
! subroutine xmtoxc; and inverse mapping routine (the inverse jacobian
! being only the pseudo-inverse) is carried out in xctoxm.
!     Since the tabulation array involves an allocatable array that we
! allocate and initialize in this this subroutine, for cleanliness, we need
! also to provide a finalization routine; this is to be found as subr.
! fini_map_table, whering the arrays are deallocated once more.
!=============================================================================
use pkind, only: dp
use pietc, only: zero,half,one,s60,r2,r3
use cutable
implicit none
integer,                          intent(in ):: norh_s,n_s
real(dp),                         intent(in ):: rb_s,gamma_s,kay_s,p_s
real(dp),dimension(3,0:n_s,0:n_s),intent(in ):: xcs_s
!-----------------------------------------------------------------------------
real(dp),parameter:: r6=r2*r3,r6o6=r6/6,r18o6=3*r2/6,r12o6=2*r3/6,r24o6=2*r6/6
!=============================================================================
call fini_map_table ! <- start with a clean slate
norh=norh_s; norhp=norh+1; norpp=norhp*2; n=n_s
allocate( wx(norpp),wy(norpp), wxd(norpp),wyd(norpp) )
allocate( linex(3,norpp),liney(3,norpp),lineg(3,norpp) )
allocate( xcs(3,0:n,0:n) )
rb=rb_s
gamma=gamma_s; gammac=one-gamma
kay=kay_s
p=p_s
xcs=xcs_s
dx=half/n

! Initialize the rotation matrices:
rotd(:,1)=(/-half,s60,zero/)
rotd(:,2)=(/s60,half,zero/)
rotd(:,3)=(/zero,zero,one/)

rotx(:,1)=(/-one,zero,zero/)
rotx(:,2)=(/zero,one,zero/)
rotx(:,3)=(/zero,zero,one/)

roty(:,1)=(/one,zero,zero/)
roty(:,2)=(/zero,-one,zero/)
roty(:,3)=(/zero,zero,one/)

rote(:,1)=(/r24o6,-r6o6,r6o6/)
rote(:,2)=(/zero,r18o6,r18o6/)
rote(:,3)=(/-r12o6,-r12o6,r12o6/)

rotpan(:,1,1)=(/zero,one,zero/)
rotpan(:,2,1)=(/zero,zero,one/)
rotpan(:,3,1)=(/one,zero,zero/)
!--
rotpan(:,1,2)=(/-one,zero,zero/)
rotpan(:,2,2)=(/zero,zero,one/)
rotpan(:,3,2)=(/zero,one,zero/)
!--
rotpan(:,1,3)=(/zero,-one,zero/)
rotpan(:,2,3)=(/zero,zero,one/)
rotpan(:,3,3)=(/-one,zero,zero/)
!--
rotpan(:,1,4)=(/one,zero,zero/)
rotpan(:,2,4)=(/zero,zero,one/)
rotpan(:,3,4)=(/zero,-one,zero/)
!--
rotpan(:,1,5)=(/zero,one,zero/)
rotpan(:,2,5)=(/-one,zero,zero/)
rotpan(:,3,5)=(/zero,zero,one/)
!--
rotpan(:,1,6)=(/zero,one,zero/)
rotpan(:,2,6)=(/one,zero,zero/)
rotpan(:,3,6)=(/zero,zero,-one/)

end subroutine init_map_table

!=============================================================================
subroutine fini_map_table
!=============================================================================
! "Finalize" the map table in module cutable by deallocating all allocatable 
!  arrays and resetting n and norh to 0.
!=============================================================================
use cutable
implicit none
if(allocated(wx))   deallocate( wx,wy,wxd,wyd )
if(allocated(linex))deallocate( linex,liney,lineg )
if(allocated(xcs))  deallocate( xcs )
n=0
norh=0
end subroutine fini_map_table

!=============================================================================
subroutine bigaints(norh,n,ns,kay0,gamma,tau, &
      rats,arint,ajac0, wwb,jac0, &
      biga,dbar,kay,kays)
!=============================================================================
! Carry out the integrals over the fundamental region of biga, as well as its 
! tau-derivative, and hence estimate dbar, kay and kay**2.
!=============================================================================
use pkind,only: dp
use pietc, only: zero,half,one,pi
implicit none
integer,                    intent(in ):: norh,n,ns     
real(dp),                   intent(in ):: kay0,gamma,tau
real(dp),dimension(ns),     intent(in ):: rats,arint,ajac0
real(dp),dimension(0:n,0:n),intent(in ):: wwb,jac0
real(dp),                   intent(out):: biga,dbar,kay,kays

!-----------------------------------------------------------------------------
integer, parameter      :: ngr=48
real(dp),parameter      :: pi4=pi*4
real(dp)                :: dx,dxx,tauc,gammac,p,r,dbigadt,jac0tauc,ljac0,kay0s
real(dp),dimension(ns)  :: ajac0tauc,aljac0
integer                 :: is,ix,iy
!=============================================================================
dx=half/n; dxx=dx**2; tauc=one-tau; gammac=one-gamma; p=(one-gamma*tau)/gammac
kay0s=kay0**2
do is=1,ns
   r=rats(is)
   ajac0tauc(is)=ajac0(is)**tauc ! Integrand for biga
   aljac0(is)=log(ajac0(is))     ! integrand for -dbigadt
enddo

biga   = dot_product(arint,ajac0tauc)
dbigadt=-dot_product(arint,ajac0tauc*aljac0)

do iy=0,n
   do ix=max(1,iy),n
      jac0tauc=jac0(ix,iy)**tauc
      ljac0=log(jac0(ix,iy))
      biga   =biga   +wwb(ix,iy)*jac0tauc*dxx
      dbigadt=dbigadt-wwb(ix,iy)*jac0tauc*dxx*ljac0
   enddo
enddo
biga   =biga*ngr/pi4
dbigadt=dbigadt*ngr/pi4
dbar   =-dbigadt/biga
kays   =kay0s**tauc/(biga*gammac*p)
kay    =sqrt(kays)
end subroutine bigaints

!=============================================================================
subroutine senslap(norh,n,d,dd,hcontra,gcontra,sens)
!=============================================================================
! Get the array, sens, of sensitivities of the diagnosed laplacian of chi at 
! each (ix,iy) to a change in the value of the field chi at the same point.
!=============================================================================
use pkind, only: dp
implicit none
integer,                        intent(in ):: norh,n
real(dp),dimension(  norh),     intent(in ):: d
real(dp),dimension(0:norh),     intent(in ):: dd
real(dp),dimension(2,0:n,0:n),  intent(in ):: hcontra
real(dp),dimension(2,2,0:n,0:n),intent(in ):: gcontra
real(dp),dimension(0:n,0:n),    intent(out):: sens
!-----------------------------------------------------------------------------
real(dp),dimension(0:n,0:n):: chi
integer                    :: ix,iy
!=============================================================================
sens=0
chi=0
do iy=0,n
   do ix=max(1,iy),n
      chi(ix,iy)=1
      call plapchi(norh,n,ix,iy,d,dd,hcontra(:,ix,iy),gcontra(:,:,ix,iy), &
           chi,sens(ix,iy))
      chi(ix,iy)=0
   enddo
enddo
! Change the sign of sens to make the coefficients positive:
sens=-sens
end subroutine senslap

!=============================================================================
subroutine plapchi(norh,n,ix,iy,d,dd,hcontra,gcontra,chi,lapchi)
!=============================================================================
! Get the laplacian of chi at the grid point, (ix,iy).
! Derivs. in map space are obtained by differencing along lines as follows:
! gx, hx: differencing along lines in the direction (1,0)
! gy, hy: differencing along lines in the direction (0,1)
! gu, hu: differencing along lines in the direction (1,1)
! gv, hv: differencing along lines in the direction (-1,1)
!=============================================================================
use pkind, only: dp
use pietc, only: half
implicit none
integer,                    intent(in ):: norh,n,ix,iy
real(dp),dimension(  norh), intent(in ):: d
real(dp),dimension(0:norh), intent(in ):: dd
real(dp),dimension(2),      intent(in ):: hcontra
real(dp),dimension(2,2),    intent(in ):: gcontra
real(dp),dimension(0:n,0:n),intent(in ):: chi
real(dp),                   intent(out):: lapchi
!-----------------------------------------------------------------------------
real(dp)                  :: dx,dxx,chit,chip,chim,gx,gy,gu,gv,hx,hy,hu,hv
real(dp),dimension(  norh):: line1
real(dp),dimension(0:norh):: line0
real(dp),dimension(2)     :: dchidz
real(dp),dimension(2,2)   :: ddchidzz
integer                   :: i,jx,jy
!=============================================================================
if(iy<0) stop 'In plapchi; iy < 0 is out of bounds'
if(ix<iy)stop 'In plapchi; ix < iy is out of bounds'
if(ix>n)stop  'In plapchi; ix > n is out of bounds'
dx=half/n; dxx=dx*dx
chit=chi(ix,iy)
if(iy==0.and.ix==0)then
   lapchi=0
   return
endif
line0(0)=chit
do i=1,norh
   jx=ix+i; jy=iy; call image(n,jx,jy); chip=chi(jx,jy)
   jx=ix-i; jy=iy; call image(n,jx,jy); chim=chi(jx,jy)
   line1(i)=chip-chim
   line0(i)=chip+chim
enddo
gx=dot_product(line1,d )/dx
hx=dot_product(line0,dd)/dxx
do i=1,norh
   jx=ix; jy=iy+i; call image(n,jx,jy); chip=chi(jx,jy)
   jx=ix; jy=iy-i; call image(n,jx,jy); chim=chi(jx,jy)
   line1(i)=chip-chim
   line0(i)=chip+chim
enddo
gy=dot_product(line1,d )/dx
hy=dot_product(line0,dd)/dxx
do i=1,norh
   jx=ix+i; jy=iy+i; call image(n,jx,jy); chip=chi(jx,jy)
   jx=ix-i; jy=iy-i; call image(n,jx,jy); chim=chi(jx,jy)
   line1(i)=chip-chim
   line0(i)=chip+chim
enddo
gu=dot_product(line1,d )/dx
hu=dot_product(line0,dd)/dxx
do i=1,norh
   jx=ix-i; jy=iy+i; call image(n,jx,jy); chip=chi(jx,jy)
   jx=ix+i; jy=iy-i; call image(n,jx,jy); chim=chi(jx,jy)
   line1(i)=chip-chim
   line0(i)=chip+chim
enddo
gv=dot_product(line1,d )/dx
hv=dot_product(line0,dd)/dxx

! Special cases exploit approximations, gx ~ (gu-gv)/2 and 
!    hx+hy ~ (hu+hv)/2 (Laplacians)
if(ix<=norh.and.iy==0)then ! Close in, on the base line
   dchidz(1)=(gu-gv)/2
   dchidz(2)=gy
   ddchidzz(1,1)=(hu+hv)/2-hy
   ddchidzz(1,2)=(hu-hv)/4
   ddchidzz(2,1)=(hu-hv)/4
   ddchidzz(2,2)=hy
elseif(ix<=norh.and.iy==ix)then ! Close in, on the diagonal:
   dchidz(1)=gx
   dchidz(2)=gy
   ddchidzz(1,1)=hx
   ddchidzz(1,2)=(hx+hy-hv)/2
   ddchidzz(2,1)=(hx+hy-hv)/2
   ddchidzz(2,2)=hy
else                            ! Generic case:
   dchidz(1)=gx
   dchidz(2)=gy
   ddchidzz(1,1)=hx
   ddchidzz(1,2)=(hu-hv)/4
   ddchidzz(2,1)=(hu-hv)/4
   ddchidzz(2,2)=hy
endif
lapchi=dot_product(hcontra,dchidz) +gcontra(1,1)*ddchidzz(1,1) &
      +2*gcontra(1,2)*ddchidzz(1,2)+gcontra(2,2)*ddchidzz(2,2)
!print*,'lapchi=',lapchi
end subroutine plapchi

!=============================================================================
subroutine fgradchi(norh,n,d,dxcdz,gcontra,chi,dchidxc)
!=============================================================================
! Evaluate the field of gradients of chi in the gridded fundamental triangle
use pkind, only: dp
implicit none
integer,                        intent(in ):: norh,n
real(dp),dimension(norh),       intent(in ):: d
real(dp),dimension(3,2,0:n,0:n),intent(in ):: dxcdz
real(dp),dimension(2,2,0:n,0:n),intent(in ):: gcontra
real(dp),dimension(    0:n,0:n),intent(in ):: chi
real(dp),dimension(3,  0:n,0:n),intent(out):: dchidxc
!-----------------------------------------------------------------------------
integer:: ix,iy
!=============================================================================
do iy=0,n
   do ix=iy,n
      call pgradchi(norh,n,ix,iy,d, &
           dxcdz(:,:,ix,iy),gcontra(:,:,ix,iy), chi,dchidxc(:,ix,iy))
   enddo
enddo
end subroutine fgradchi

!=============================================================================
subroutine pgradchi(norh,n,ix,iy,d,dxcdz,gcontra,chi,dchidxc)
!=============================================================================
! Get the gradient of chi as a cartesian 3-vector tangential to the sphere
! at the grid point, (ix,iy) of the fundamental triangle.
! Components in map space are obtained by differencing along lines as follows:
! gx: differencing along lines in the direction (1,0)
! gy: differencing along lines in the direction (0,1)
! gu: differencing along lines in the direction (1,1)
! gv: differencing along lines in the direction (-1,1)
! From suitable combinations that avoid using involving the vertex,
! the two map-space derivatives are first obtained (covariant vector, dchidz),
! converted to contravariant (by operating with metric tensor gcontra),
! and finally converted to cartesian 3-vector, tangential to the sphere,
! by operating with dxcdz.
!=============================================================================
use pkind, only: dp
use pietc, only: half
implicit none
integer,                    intent(in ):: norh,n,ix,iy
real(dp),dimension(norh),   intent(in ):: d
real(dp),dimension(3,2),    intent(in ):: dxcdz
real(dp),dimension(2,2),    intent(in ):: gcontra
real(dp),dimension(0:n,0:n),intent(in ):: chi
real(dp),dimension(3),      intent(out):: dchidxc
!------------------------------------------------------------------------------
real(dp),dimension(  norh):: line1
real(dp),dimension(2)     :: grads,dchidz
real(dp)                  :: dx,gx,gy,gu,gv,chim,chip
integer                   :: i,jx,jy
!==============================================================================
if(iy<0 )stop 'In pgradchi; iy < 0 is out of bounds'
if(ix<iy)stop 'In pgradchi; ix < iy is out of bounds'
if(ix>n )stop 'In pgradchi; ix > n is out of bounds'
dx=half/n
if(iy==0 .and. ix==0)then
   dchidxc=0
   return
endif
do i=1,norh
   jx=ix+i; jy=iy; call image(n,jx,jy); chip=chi(jx,jy)
   jx=ix-i; jy=iy; call image(n,jx,jy); chim=chi(jx,jy)
   line1(i)=chip-chim
enddo
gx=dot_product(line1,d )/dx
do i=1,norh
   jx=ix; jy=iy+i; call image(n,jx,jy); chip=chi(jx,jy)
   jx=ix; jy=iy-i; call image(n,jx,jy); chim=chi(jx,jy)
   line1(i)=chip-chim
enddo
gy=dot_product(line1,d )/dx
if(ix<=norh.and.iy==0)then ! Close in, on the base line
   do i=1,norh
      jx=ix+i; jy=iy+i; call image(n,jx,jy); chip=chi(jx,jy)
      jx=ix-i; jy=iy-i; call image(n,jx,jy); chim=chi(jx,jy)
      line1(i)=chip-chim
   enddo
   gu=dot_product(line1,d )/dx
   do i=1,norh
      jx=ix-i; jy=iy+i; call image(n,jx,jy); chip=chi(jx,jy)
      jx=ix+i; jy=iy-i; call image(n,jx,jy); chim=chi(jx,jy)
      line1(i)=chip-chim
   enddo
   gv=dot_product(line1,d )/dx
   
   dchidz(1)=(gu-gv)/2
   dchidz(2)=gy
else                            ! Generic case:
   dchidz(1)=gx
   dchidz(2)=gy
endif
grads=matmul(gcontra,dchidz)
dchidxc=matmul(dxcdz,grads)
end subroutine  pgradchi

!=============================================================================
subroutine fmetric(norh,n,d,dd, xc, dxcdz,dacda,gco,gcontra, &
      ddacdaz,hcontra)
!=============================================================================
! Construct the field of metric tensors and related quantities at all
! grid points of the fundamental triangle
!=============================================================================
use pkind, only: dp
implicit none
integer,                        intent(in ):: norh,n
real(dp),dimension(  norh),     intent(in ):: d
real(dp),dimension(0:norh),     intent(in ):: dd
real(dp),dimension(3,  0:n,0:n),intent(in ):: xc
real(dp),dimension(3,2,0:n,0:n),intent(out):: dxcdz
real(dp),dimension(    0:n,0:n),intent(out):: dacda
real(dp),dimension(2,2,0:n,0:n),intent(out):: gco,gcontra
real(dp),dimension(2,  0:n,0:n),intent(out):: ddacdaz,hcontra
!-----------------------------------------------------------------------------
integer:: ix,iy
!=============================================================================
do iy=0,n
   do ix=iy,n
      call pmetric(norh,n,ix,iy,d,dd, xc, &
            dxcdz(:,:,ix,iy),dacda(ix,iy),gco(:,:,ix,iy),gcontra(:,:,ix,iy),&
            ddacdaz(:,ix,iy),hcontra(:,ix,iy))
   enddo
enddo
end subroutine fmetric

!=============================================================================
subroutine pmetric(norh,n,ix,iy,d,dd, xc, dxcdz,dacda,gco,gcontra, &
      ddacdaz,hcontra)
!=============================================================================
! Get the metrical properties at a point, which derive from 1st and 2nd
! derivatives wrt map coordinates of the cartesian 3-vector positions, xc.
! Derivs. in map space are obtained by differencing along lines as follows:
! gx, hx: differencing along lines in the direction (1,0)
! gy, hy: differencing along lines in the direction (0,1)
! gu, hu: differencing along lines in the direction (1,1)
! gv, hv: differencing along lines in the direction (-1,1)
!=============================================================================
use pkind, only: dp
use pietc, only: zero,half,one
use peuc,  only: triple_product
use pmat,  only: inv
implicit none
integer,                      intent(in ):: norh,n,ix,iy
real(dp),dimension(  norh),   intent(in ):: d
real(dp),dimension(0:norh),   intent(in ):: dd
real(dp),dimension(3,0:n,0:n),intent(in ):: xc
real(dp),dimension(3,2),      intent(out):: dxcdz
real(dp),                     intent(out):: dacda
real(dp),dimension(2,2),      intent(out):: gco,gcontra
real(dp),dimension(2  ),      intent(out):: ddacdaz,hcontra
!-----------------------------------------------------------------------------
real(dp)                    :: dx,dxx
real(dp),dimension(3)       :: xct,xcp,xcm,gx,gy,gu,gv,hx,hy,hu,hv
real(dp),dimension(3,  norh):: line1
real(dp),dimension(3,0:norh):: line0
real(dp),dimension(3,3)     :: mat
real(dp),dimension(3,2,2)   :: ddxcdzz
real(dp),dimension(2,2,2)   :: dgcodz,dgcontradz
integer                     :: i,jx,jy
!=============================================================================
if(iy<0) stop 'In pmetric; iy < 0 is out of bounds'
if(ix<iy)stop 'In pmetric; ix < iy is out of bounds'
if(ix>n)stop  'In pmetric; ix > n is out of bounds'
dx=half/n; dxx=dx*dx
xct=xc(:,ix,iy)
if(iy==0.and.ix==0)then
   dxcdz=0
   dacda=0
   gco=0
   gcontra=0
   ddacdaz=0
   hcontra=0
   return
endif
line0(:,0)=xct
do i=1,norh
   jx=ix+i; jy=iy; call imagev(n,jx,jy,mat); xcp=matmul(mat,xc(:,jx,jy))
   jx=ix-i; jy=iy; call imagev(n,jx,jy,mat); xcm=matmul(mat,xc(:,jx,jy))
   line1(:,i)=xcp-xcm
   line0(:,i)=xcp+xcm
enddo
gx=matmul(line1,d )/dx
hx=matmul(line0,dd)/dxx
do i=1,norh
   jx=ix; jy=iy+i; call imagev(n,jx,jy,mat); xcp=matmul(mat,xc(:,jx,jy))
   jx=ix; jy=iy-i; call imagev(n,jx,jy,mat); xcm=matmul(mat,xc(:,jx,jy))
   line1(:,i)=xcp-xcm
   line0(:,i)=xcp+xcm
enddo
gy=matmul(line1,d )/dx
hy=matmul(line0,dd)/dxx
do i=1,norh
   jx=ix+i; jy=iy+i; call imagev(n,jx,jy,mat); xcp=matmul(mat,xc(:,jx,jy))
   jx=ix-i; jy=iy-i; call imagev(n,jx,jy,mat); xcm=matmul(mat,xc(:,jx,jy))
   line1(:,i)=xcp-xcm
   line0(:,i)=xcp+xcm
enddo
gu=matmul(line1,d )/dx
hu=matmul(line0,dd)/dxx
do i=1,norh
   jx=ix-i; jy=iy+i; call imagev(n,jx,jy,mat); xcp=matmul(mat,xc(:,jx,jy))
   jx=ix+i; jy=iy-i; call imagev(n,jx,jy,mat); xcm=matmul(mat,xc(:,jx,jy))
   line1(:,i)=xcp-xcm
   line0(:,i)=xcp+xcm
enddo
gv=matmul(line1,d )/dx
hv=matmul(line0,dd)/dxx

if(ix<=norh.and.iy==0)then ! Close in, on the base line
   dxcdz(:,1)=(gu-gv)/2
   dxcdz(:,2)=gy
   ddxcdzz(:,1,1)=(hu+hv)/2-hy
   ddxcdzz(:,1,2)=(hu-hv)/4
   ddxcdzz(:,2,1)=(hu-hv)/4
   ddxcdzz(:,2,2)=hy
elseif(ix<=norh.and.iy==ix)then ! Close in, on the diagonal:
   dxcdz(:,1)=gx
   dxcdz(:,2)=gy
   ddxcdzz(:,1,1)=hx
   ddxcdzz(:,1,2)=(hx+hy-hv)/2
   ddxcdzz(:,2,1)=(hx+hy-hv)/2
   ddxcdzz(:,2,2)=hy
else                            ! Generic case:
   dxcdz(:,1)=gx
   dxcdz(:,2)=gy
   ddxcdzz(:,1,1)=hx
   ddxcdzz(:,1,2)=(hu-hv)/4
   ddxcdzz(:,2,1)=(hu-hv)/4
   ddxcdzz(:,2,2)=hy
endif
gco=matmul(transpose(dxcdz),dxcdz) ! Covariant metric tensor
gcontra=gco; call inv(gcontra) ! Contravariant metric tensor
dacda=triple_product(xct,dxcdz(:,1),dxcdz(:,2))
ddacdaz(1)=triple_product(xct,ddxcdzz(:,1,1),dxcdz(:,2)) &
          +triple_product(xct,dxcdz(:,1),ddxcdzz(:,2,1))
ddacdaz(2)=triple_product(xct,ddxcdzz(:,1,2),dxcdz(:,2)) &
          +triple_product(xct,dxcdz(:,1),ddxcdzz(:,2,2))
dgcodz(:,:,1)=matmul(transpose(ddxcdzz(:,:,1)),dxcdz) &
             +matmul(transpose(dxcdz),ddxcdzz(:,:,1))
dgcodz(:,:,2)=matmul(transpose(ddxcdzz(:,:,2)),dxcdz) &
             +matmul(transpose(dxcdz),ddxcdzz(:,:,2))
dgcontradz(:,:,1)=-matmul(gcontra,matmul(dgcodz(:,:,1),gcontra))
dgcontradz(:,:,2)=-matmul(gcontra,matmul(dgcodz(:,:,2),gcontra))
hcontra=matmul(gcontra,ddacdaz)/dacda + dgcontradz(1,:,1)+dgcontradz(2,:,2)

end subroutine pmetric

!=============================================================================
subroutine xmtoxc_s(xm_s,xc_s,xcd_s,ipan)!                            [xmtoxc]
!=============================================================================
use pkind,  only: dp
implicit none
real,dimension(2),intent(in ):: xm_s
real,dimension(3),intent(out):: xc_s
real,dimension(3,2),intent(out):: xcd_s
integer,            intent(in ):: ipan
!-----------------------------------------------------------------------------
real(dp),dimension(2):: xm
real(dp),dimension(3):: xc
real(dp),dimension(3,2):: xcd
!=============================================================================
xm=xm_s
call xmtoxc(xm,xc,xcd,ipan)
xc_s =xc
xcd_s=xcd
end subroutine xmtoxc_s

!=============================================================================
subroutine xmtoxc(xm,xc,xcd,ipan)!                                   [xmtoxc]
!=============================================================================
use pkind,  only: dp
use pietc,  only: zero,half,one
use pbint,  only: getsplint
use pbend,  only: wbend
use cutable,only: norhp,norpp,n,rb,gammac,kay,p,dx,wx,wy,wxd,wyd, &
                  linex,liney,lineg,xcs,rotx,roty,rotd,rote,rotpan
implicit none
real(dp),dimension(2),  intent(in ):: xm
real(dp),dimension(3),  intent(out):: xc
real(dp),dimension(3,2),intent(out):: xcd
integer,                intent(in ):: ipan
!-----------------------------------------------------------------------------
logical                :: kx,ky,kd
integer                :: i,ix,iy,j,jx,jy,jys
real(dp)               :: x,y,r,rx,ry,rr,theta,ctheta,stheta,bigr,dbigrdx,&
                          dbigrdy,s,wa,wb,wbd,dwbdx,dwbdy,dthetadx,dthetady
real(dp),dimension(2)  :: xmt
real(dp),dimension(3)  :: xca
real(dp),dimension(3,2):: xcad
real(dp),dimension(3,3):: rot,mat
!=============================================================================
xmt=xm
kx=xmt(1)>zero;   if(kx)xmt(1)=-xmt(1)
ky=xmt(2)>zero;   if(ky)xmt(2)=-xmt(2)
kd=xmt(2)>xmt(1); if(kd)xmt    =xmt(2:1:-1)
x=half*(one+xmt(1))
y=half*(one+xmt(2))
rr=x**2+y**2
r=sqrt(rr)
rx=x/dx; ry=y/dx
ix=floor(rx); iy=floor(ry)
rx=rx-ix;     ry=ry-iy
call getsplint(norhp,rx,wx,wxd); call getsplint(norhp,ry,wy,wyd)
do j=1,norpp
   jys=iy-norhp+j
   do i=1,norpp
      jx=ix-norhp+i
      jy=jys
      call imagev(n,jx,jy,mat)
      linex(:,i)=matmul(mat,xcs(:,jx,jy))
   enddo
   liney (:,j)=matmul(linex,wx)
   lineg (:,j)=matmul(linex,wxd)/dx
enddo
xc=matmul(liney,wy)
xcd(:,1)=matmul(lineg,wy)
xcd(:,2)=matmul(liney,wyd)/dx
if(r<rb)then
   if(r==0)then
      xca=(/zero,zero,one/)
      xcad=0
   else ! Blend asymptotic solutions (xca, xcad), with grid-based solution:
      theta=atan2(y,x)/gammac; ctheta=cos(theta); stheta=sin(theta)
      dthetadx=-y/(gammac*rr); dthetady=x/(gammac*rr)
      bigr=gammac*kay*r**p
      dbigrdx=p*x*bigr/rr;    dbigrdy=p*y*bigr/rr
      xca(1)=bigr*ctheta
      xca(2)=bigr*stheta
      xca(3)=one
      xcad(1,1)=dbigrdx*ctheta-bigr*stheta*dthetadx
      xcad(1,2)=dbigrdy*ctheta-bigr*stheta*dthetady
      xcad(2,1)=dbigrdx*stheta+bigr*ctheta*dthetadx
      xcad(2,2)=dbigrdy*stheta+bigr*ctheta*dthetady
      xcad(3,:)=0
      call wbend((r*2-rb)/rb,wb,wbd); wa=one-wb
      dwbdx=wbd*x*2/(rb*r); dwbdy=wbd*y*2/(rb*r)
      xcd=wa*xcad+wb*xcd
      xcd(:,1)=xcd(:,1)+dwbdx*(xc-xca)
      xcd(:,2)=xcd(:,2)+dwbdy*(xc-xca)
      xc=wa*xca+wb*xc
   endif
endif
s=sqrt(dot_product(xc,xc)); xc=xc/s                ! Normalize xc
s=dot_product(xc,xcd(:,1)); xcd(:,1)=xcd(:,1)-xc*s ! Orthogonalize, wrt xc,
s=dot_product(xc,xcd(:,2)); xcd(:,2)=xcd(:,2)-xc*s ! both derivative components
xcd=xcd/2 ! Derivatives wrt final map coordinates are smaller by this factor
          ! because, although here the square is spanned by x and y in [0,1],
          ! the final coordinates will span each square in range [-1:1].
if(kd)then
   xc =matmul(rotd,xc)
   xcd=matmul(rotd,xcd); xcd=xcd(:,2:1:-1)
endif
xc =matmul(rote,xc)
xcd=matmul(rote,xcd)
if(ky)then
   xc =matmul(roty,xc)
   xcd=matmul(roty,xcd); xcd(:,2)=-xcd(:,2)
endif
if(kx)then
   xc =matmul(rotx,xc)
   xcd=matmul(rotx,xcd); xcd(:,1)=-xcd(:,1)
endif
rot=rotpan(:,:,ipan)
xc =matmul(rot,xc)
xcd=matmul(rot,xcd)

end subroutine xmtoxc

!=============================================================================
subroutine xctoxm_s(xc_s,xm_s,xmd_s,ipan)!                            [xctoxm]
!=============================================================================
use pkind,   only: dp
implicit none
real,dimension(3),  intent(in ):: xc_s
real,dimension(2),  intent(out):: xm_s
real,dimension(2,3),intent(out):: xmd_s
integer,            intent(out):: ipan
!-----------------------------------------------------------------------------
real(dp),dimension(3)  :: xc
real(dp),dimension(2)  :: xm
real(dp),dimension(2,3):: xmd
!==============================================================================
xc=xc_s
call xctoxm(xc,xm,xmd,ipan)
xm_s=xm
xmd_s=xmd
end subroutine xctoxm_s


!=============================================================================
subroutine xctoxm(xc,xm,xmd,ipan)!                                    [xctoxm]
!=============================================================================
use pkind,   only: dp
use pietc,   only: zero,one
use pmat,    only: inv
use cutable, only: rotpan
implicit none
real(dp),dimension(3),  intent(in ):: xc
real(dp),dimension(2),  intent(out):: xm
real(dp),dimension(2,3),intent(out):: xmd
integer,                intent(out):: ipan
!-----------------------------------------------------------------------------
integer,parameter      :: nit=10,nitm3=nit-3 !<- allowance of newton iterations
real(dp),parameter     :: cnorm=1.e-8
real(dp)               :: s,norm
real(dp),dimension(nit):: urc    ! <- underrelaxation coefficients
real(dp),dimension(3)  :: xct,xca,res
real(dp),dimension(6)  :: ax
real(dp),dimension(3,2):: xcd
real(dp),dimension(2,2):: gco
real(dp),dimension(3,3):: rot
integer                :: it,jpan
integer,dimension(1)   :: ii
logical                :: kx,ky
data urc/.1_dp,.2_dp,.4_dp,nitm3*one/
!=============================================================================
ax=(/xc(1),xc(2),-xc(1),-xc(2),xc(3),-xc(3)/); ii=maxloc(ax); ipan=ii(1)
s=sqrt(dot_product(xc,xc))
xct=xc/s

! Apply the rotation to xct that maps it into the polar panel, 5, temporarily
rot=matmul(rotpan(:,:,5),transpose(rotpan(:,:,ipan)))
xct=matmul(rot,xct)

! Apply reflections to bring xct into the quadrant of the square face
! next to the vertex at on the ray, (-1,-1,+1), whose xm=(-1,-1).
! Note the 90 degree rotation, owing to the orientation of panel 5:
kx=xct(2)>zero; if(kx)xct(2)=-xct(2)
ky=xct(1)<zero; if(ky)xct(1)=-xct(1)
xm(1)=xct(2)/xct(3)
xm(2)=-xct(1)/xct(3)

! Trim to eliminate any round-off error:
xm(1)=max(xm(1),-one)
xm(2)=max(xm(2),-one)

! Special treatment when xm IS the vertex:
if(xm(1)==-one .and. xm(2)==-one)then
   xmd=0
   if(kx)xm(1)=-xm(1)
   if(ky)xm(2)=-xm(2)
   return
endif

! Perform newton iterations to solve for xm, terminating when convergence
! criterion, norm<=cnorm, is met, or until all iterations are exhausted.
! Use underrelaxation for the first few steps as a precaution.
do it=1,nit
   call xmtoxc(xm,xca,xcd,5)
   gco=matmul(transpose(xcd),xcd)
   call inv(gco)
   xmd=matmul(gco,transpose(xcd))! xmd is the generalized inverse of xcd
   res=xca-xct
   norm=sqrt(dot_product(res,res))
   xm=xm-urc(it)*matmul(xmd,res)
   if(norm<=cnorm)exit
enddo
if(it>nit)print'("In xctoxm; possible convergence problem in Newton its.")'
if(kx)then
   xm(1)=-xm(1)
   xmd(1,:)=-xmd(1,:)
endif
if(ky)then
   xm(2)=-xm(2)
   xmd(2,:)=-xmd(2,:)
endif
if(kx)xmd(:,2)=-xmd(:,2)
if(ky)xmd(:,1)=-xmd(:,1)
xmd=matmul(xmd,rot)
end subroutine xctoxm

!=============================================================================
subroutine image(n,jx,jy)
!=============================================================================
! Find and return the image indices (jx,jy) inside the fundamental triangle
! of size n when given indices (jx,jy) that are initially not necessarily
! inside that region. Images are formed by kaleidoscopic reflections about
! the sides of the fundamental triangle, 0 <= jy <= jx <= n
!=============================================================================
implicit none
integer,intent(in   ):: n
integer,intent(inout):: jx,jy
!-----------------------------------------------------------------------------
integer,parameter:: nit=20
integer          :: it,count,j
!=============================================================================
count=0
do it=1,nit ! Image must pass the three tests in a row to be valid
   if(jy<0)then;  count=0; jy=-jy;            endif
   count=count+1; if(count>=3)return
   if(jy>jx)then; count=0; j=jx; jx=jy; jy=j; endif
   count=count+1; if(count>=3)return
   if(jx>n)then;  count=0; jx=n*2-jx;         endif
   count=count+1; if(count>=3)return
enddo
stop 'In image; no solution found'
end subroutine image

!=============================================================================
subroutine timage(n3,jx,jy)
!=============================================================================
! For a polyhedral face of triangular shape. Face center a gridpoint when
! n3 is divisible by 3.
! Find and return the image indices (jx,jy) inside the fundamental triangle
! of size n when given indices (jx,jy) that are initially not necessarily
! inside that region. Images are formed by kaleidoscopic reflections about
! the sides of the fundamental triangle, 0 <= jy <= jx <= (n3*2-jy)/2
!=============================================================================
implicit none
integer,intent(in   ):: n3
integer,intent(inout):: jx,jy
!-----------------------------------------------------------------------------
integer,parameter:: nit=20
integer          :: it,count,j,n6
!=============================================================================
n6=n3*2
count=0
do it=1,nit ! Image must pass the three tests in a row to be valid
   if(jy<0)then;  count=0; jx=jx+jy; jy=-jy;  endif
   count=count+1; if(count>=3)return
   if(jy>jx)then; count=0; j=jx; jx=jy; jy=j; endif
   count=count+1; if(count>=3)return
   if(2*jx>n6-jy)then;  count=0; jx=n6-jx-jy; endif
   count=count+1; if(count>=3)return
enddo
stop 'In timage; no solution found'
end subroutine timage

!=============================================================================
subroutine himage(n2,jx,jy)
!=============================================================================
! For a polyhedral face of hexagonal shape. Face center a gridpoint when
! n2 is divisible by 2.
! Find and return the image indices (jx,jy) inside the fundamental triangle
! of base size n2 when given indices (jx,jy) that are initially not necessarily
! inside that region. Images are formed by kaleidoscopic reflections about
! the sides of the fundamental triangle, 0 <= jy; 0<= jx <= (n2*2-jy)/2
!=============================================================================
implicit none
integer,intent(in   ):: n2
integer,intent(inout):: jx,jy
!-----------------------------------------------------------------------------
integer,parameter:: nit=20
integer          :: it,count,j,n4
!=============================================================================
n4=n2*2
count=0
do it=1,nit ! Image must pass the three tests in a row to be valid
   if(jy<0)then; count=0; jx=jx+jy; jy=-jy;  endif
   count=count+1; if(count>=3)return
   if(jx<0)then; count=0; jy=jx+jy; jx=-jx;  endif
   count=count+1; if(count>=3)return
   if(2*jx>n4-jy)then; count=0; jx=n4-jx-jy; endif
   count=count+1; if(count>=3)return
enddo
stop 'In himage; no solution found'
end subroutine himage

!=============================================================================
subroutine imagev(n,jx,jy,mat)
!=============================================================================
! Find and return the image indices (jx,jy) inside the fundamental triangle
! of size n when given indices (jx,jy) that are initially not necessarily
! inside that region. Images are formed by kaleidoscopic reflections about
! the sides of the fundamental triangle, 0 <= jy <= jx <= n.
! Also, return an orthogonal  matrix representing the trasformation needed 
! to project the three vector with the original indices into the 
! vector with the image indices
! mat1: Reflect about base edge, y=0
! mat2: reflect about diagonal edge x=y (at 60 degrees in cartesian projection)
! mat3: reflect about edge x=half.
!=============================================================================
use pkind, only: dp
use pietc, only: zero,third,half,one,mone,s60,r2
implicit none
integer,                intent(in   ):: n
integer,                intent(inout):: jx,jy
real(dp),dimension(3,3),intent(  out):: mat
!-----------------------------------------------------------------------------
integer,parameter      :: nit=4
real(dp),parameter     :: mhalf=-half,ca=third,mca=-ca,sa=2*r2*third 
integer                :: i,it,count,j
real(dp),dimension(9)  :: lmat!,lmat1,lmat2,lmat3
real(dp),dimension(3,3):: mat1,mat2,mat3
integer,dimension(2)   :: ii
data ii/3,3/
data mat1/one,zero,zero,  zero,mone,zero, zero,zero,one/
data mat2/mhalf,s60,zero, s60,half,zero,  zero,zero,one/
data mat3/mca,zero,sa,    zero,one,zero,    sa,zero,ca /
!=============================================================================
lmat= (/one,zero,zero,  zero,one,zero,  zero,zero,one/); mat =reshape(lmat ,ii)
!lmat1=(/one,zero,zero,  zero,-one,zero, zero,zero,one/);mat1=reshape(lmat1,ii)
!lmat2=(/-half,s60,zero, s60,half,zero,  zero,zero,one/);mat2=reshape(lmat2,ii)
!lmat3=(/-ca,zero,sa,    zero,one,zero,    sa,zero,ca /);mat3=reshape(lmat3,ii)
count=0
do it=1,nit ! Image must pass the three tests in a row to be valid
   if(jy<0)then;  count=0; mat=matmul(mat,mat1); jy=-jy;            endif
   count=count+1; if(count>=3)return
   if(jy>jx)then; count=0; mat=matmul(mat,mat2); j=jx; jx=jy; jy=j; endif
   count=count+1; if(count>=3)return
   if(jx>n)then;  count=0; mat=matmul(mat,mat3); jx=n*2-jx;         endif
   count=count+1; if(count>=3)return
enddo
stop 'In imagev; no solution found'
end subroutine imagev

!=============================================================================
subroutine rtos(nq,rb,r, ww,s)
!=============================================================================
! For points sufficiently close (in map space) to the vertex for asymptotic
! approximations to deserve consideration, determine, for a map-space 
! radial distance, r, what should be the gridded field's weighting, ww, 
! in comparison to the asymptotic solution's weighting, 1-ww. The weight
! function of r comes from applying the Whittaker blending function over
! the r interval [rb/2,rb]. In addition, return the transformed radial
! coordinate, s, such that r=rb*s**nq. Note that s=1 maps to r=rb. For nq>1
! such a transformation means that a uniform gridding in s provides a much
! higher effective resolution near the singularity than if r were used.
!=============================================================================
use pkind, only: dp
use pietc, only: one
use pbend, only: wbend
implicit none
integer, intent(in ):: nq
real(dp),intent(in ):: rb,r
real(dp),intent(out):: ww,s
!------------------------------------------------------------------------------
real(dp):: onq
!==============================================================================
onq=one/nq
s=(r/rb)**onq
call wbend((r*2-rb)/rb,ww)
end subroutine rtos

!==============================================================================
subroutine stor(nq,rb,s, ww,r,drds)
!==============================================================================
! Invert the relationship of subr. rtos between r and s, but also supply
! the function's derivative, dr/ds, for use in radial integration near the
! vertex.
!==============================================================================
use pkind, only: dp
use pietc, only: one
use pbend, only: wbend
implicit none
integer, intent(in ):: nq
real(dp),intent(in ):: rb,s
real(dp),intent(out):: ww,r,drds
!------------------------------------------------------------------------------
r   =rb*s**nq
drds=nq*rb*s**(nq-1)
call wbend((r*2-rb)/rb,ww)
end subroutine stor

!===========================================================================
subroutine replicate_cu(xc1,xc48)
!===========================================================================
! Same effect as replicate_oh except the assumed orientation of the rotation
! group is assumed to be that which is most suited to the cube in its
! position with vertex at pole (0,0,1) and an edge running towards (1,0,0)
!===========================================================================
use pkind, only: dp
use pietc, only: zero,r2,r3
implicit none
real(dp),dimension(3),   intent(in ):: xc1
real(dp),dimension(3,48),intent(out):: xc48
!----------------------------------------------------------------------------
real(dp),parameter:: r6=r2*r3,r12=2*r3,r18=3*r2,r24=2*r6, &
                     r6o=r6/6,r12o=r12/6,r18o=r18/6,r24o=r24/6, &
                     mr18o=-r18o,mr24o=-r24o
real(dp),dimension(3,3):: rot
real(dp),dimension(3)  :: v
data rot/r6o,mr18o,r12o,  r6o,r18o,r12o,  mr24o,zero,r12o/
!=============================================================================
v=matmul(transpose(rot),xc1)
call replicate_oh(v,xc48)
xc48=matmul(rot,xc48)
end subroutine replicate_cu

!============================================================================
subroutine replicate_oh(xc1,xc48)
!============================================================================
! Given a single cartesian 3-vector, xc1, replicate this vector with respect
! to the octahedral symmetry group to form a set of 48 3 vectors (including
! xc1 itslef, of course).
!============================================================================
use pkind, only: dp
implicit none
real(dp),dimension(3),   intent(in ):: xc1
real(dp),dimension(3,48),intent(out):: xc48
!----------------------------------------------------------------------------
integer:: i,im,ip,j,jp
!============================================================================
! rotate x,y,z coordinates cyclically for 1st 3 replications:
do i=1,3
   ip=i+1; if(ip>3)ip=ip-3
   im=i-1; if(im<1)im=im+3
   xc48(i ,1)=xc1(i)
   xc48(ip,2)=xc1(i)
   xc48(im,3)=xc1(i)
enddo
! Reflect about plane x-y=0 for 2nd 3 replications:
do j=1,3
   jp=j+3
   xc48(1,jp)=xc48(2,j)
   xc48(2,jp)=xc48(1,j)
   xc48(3,jp)=xc48(3,j)
enddo
! Reflect about plane x=0 for 2nd 6 replications:
do j=1,6
   jp=j+6
   xc48(1,jp)=-xc48(1,j)
   xc48(2,jp)= xc48(2,j)
   xc48(3,jp)= xc48(3,j)
enddo
! Reflect about plane y=0 for 2nd 12 replications:
do j=1,12
   jp=j+12
   xc48(1,jp)= xc48(1,j)
   xc48(2,jp)=-xc48(2,j)
   xc48(3,jp)= xc48(3,j)
enddo
! Reflect about plane z=0 for 2nd 24 replications:
do j=1,24
   jp=j+24
   xc48(1,jp)= xc48(1,j)
   xc48(2,jp)= xc48(2,j)
   xc48(3,jp)=-xc48(3,j)
enddo
end subroutine replicate_oh

!=============================================================================
subroutine setgroup_cu(ng,grouprot)
!=============================================================================
use pkind, only: dp
implicit none
integer,                   intent(in ):: ng
real(dp),dimension(3,3,ng),intent(out):: grouprot
!-----------------------------------------------------------------------------
integer                :: i
real(dp),dimension(3,3):: unit
!=============================================================================
if(ng/=48)stop 'In setgroup_cu; ng must be 48 please'
unit=0; do i=1,3; unit(i,i)=1; enddo
do i=1,3
   call replicate_cu(unit(:,i),grouprot(:,i,:))
enddo
end subroutine setgroup_cu

end module cugen