module geometry
  use SEP
  use util
  use lloyd_mod

  implicit none
  integer, private :: nz,nx,ny
  real,    private :: dz,dx,dy

  integer, private :: nref

  real   , private :: eps
  logical, private :: kin,verbose

  !! GEOMETRY
  real,    allocatable, private :: ray(:,:,:,:)
  real,    allocatable,dimension(:,:), private  :: g11,g12,g13,g22,g23,g33
  real,    allocatable,dimension(:,:), private  :: iG11,iG12,iG13,iG22,iG23,iG33
  real,    allocatable,dimension(:,:), private  :: m11,m12,m13,m22,m23,m33
  real,    allocatable,dimension(:,:), private  :: gdetr,n1,n2,n3

  integer, allocatable,dimension(:,:  ), private  :: allmask
  real,    allocatable,dimension(:,:,:), private  :: allfields

  !! For Reference Calculation
  integer, allocatable,private :: iregion(:),nstart(:)
  real,    allocatable,private :: slice(:,:),scales(:),have_zeros(:)
  real,    allocatable,private :: vslice(:,:),allrefs(:,:)
contains

  !----------------------------------------------------------------  
  subroutine geometry_init(nz_in,nx_in,ny_in,dz_in,dx_in,dy_in,kin_in,nref_in,verbose_in)
    integer          ::    nz_in,nx_in,ny_in,                         nref_in
    real             ::                      dz_in,dx_in,dy_in
    logical          ::                                        kin_in,        verbose_in

    nz=nz_in;  nx=nx_in;  ny=ny_in
    dz=dz_in;  dx=dx_in;  dy=dy_in
    nref=nref_in;    kin=kin_in
    verbose=verbose_in;    eps=0.0001;

    !! For Calculating Geometry
    allocate(ray(3,nx,ny,3) )
    allocate(g11(nx,ny),g12(nx,ny),g13(nx,ny),g22(nx,ny),g23(nx,ny),g33(nx,ny))
    allocate(iG11(nx,ny),iG12(nx,ny),iG13(nx,ny),iG22(nx,ny),iG23(nx,ny),iG33(nx,ny))
    allocate(m11(nx,ny),m12(nx,ny),m13(nx,ny),m22(nx,ny),m23(nx,ny),m33(nx,ny))
    allocate(gdetr(nx,ny),n1(nx,ny),n2(nx,ny),n3(nx,ny),vslice(nx,ny) )
    g11=0.;g12=0.;g13=0.;g22=0.;g23=0.; g33=0.;
    iG11=0.;iG12=0.;iG13=0.;iG22=0.;iG23=0.;iG33=0.;
    m11=0.;m12=0.;m13=0.;m22=0.;m23=0.;m33=0.;
    gdetr=0.;n1=0.;n2=0.;n3=0.; vslice=0.

    allocate( allfields(nx,ny,10), allrefs(nref,10), allmask(nx,ny) )
    allfields=0.; allrefs=0.;allmask=0;

    !! For Calculating Reference Parameters
    allocate(have_zeros(10));have_zeros=0.;
    allocate(scales(10)); scales=0.
    allocate(slice(nx*ny,10)); slice=0.
    allocate(iregion(nx*ny)); iregion=0;
    allocate(nstart(10));nstart=0;

    call init_lloyd()

  end subroutine geometry_init

  !----------------------------------------------------------------  

  subroutine geometry_calc(rays,Vel,fields,refs,mask)
    !! Main subroutine dealing with geometry
    integer :: iz,inn,mask(:,:,:)
    integer, allocatable :: masky(:,:,:)
    real    :: rays(:,:,:,:),Vel(:,:,:),fields(:,:,:,:),refs(:,:,:)
    allocate( masky(nx,ny,nz) );    masky=0

    if (ny .eq. 1) then

       do iz=2,nz-1
          ray    = rays(iz-1:iz+1,:,:,:)
          vslice = Vel(iz,:,:)
          call calc_g11_2d(ray,g11);!   where (    g11  .lt. 10e-4) g11=0
          call calc_g12_2d(ray,g12);   where (abs(g12) .lt. 10e-4) g12=0
          call calc_g13_2d(ray,g13);   where (abs(g13) .lt. 10e-4) g13=0
          call calc_g22_2d(ray,g22);!   where (    g22  .lt. 10e-4) g22=0
          call calc_g23_2d(ray,g23);   where (abs(g23) .lt. 10e-4) g23=0
          call calc_g33_2d(ray,g33);!   where (    g33  .lt. 10e-4) g33=0
          call calc_determinant(g11,g12,g13,g22,g23,g33,gdetr)
          call calc_inverse(g11,g12,g13,g22,g23,g33,iG11,iG12,iG13,iG22,iG23,iG33,gdetr)
          m11=gdetr*iG11;
          m12=gdetr*iG12;
          m13=gdetr*iG13;
          m22=gdetr*iG22;
          m23=gdetr*iG23;
          m33=gdetr*iG33;
          call calc_n_2d(m11,m12,m13,n1)
          call calc_n_2d(m12,m22,m23,n2)
          call calc_n_2d(m13,m23,m33,n3)

          !! Non-stationary coefficient fields
          allfields(:,:,1) = iG13/iG33
          allfields(:,:,2) = iG23/iG33
          allfields(:,:,3) = n3/(2.*m33)
          allfields(:,:,4) = sqrt(1./(vslice**2*iG33))
          allfields(:,:,5) = sqrt(iG11/iG33 - (iG13/iG33)**2)
          allfields(:,:,6) = sqrt(iG22/iG33 - (iG23/iG33)**2)
          allfields(:,:,7) = 2.*iG12/iG33 - 2.*iG13*iG23/(iG33**2)
          allfields(:,:,8) = n1/m33 - m13*n3/(m33**2)
          allfields(:,:,9) = n2/m33 - m23*n3/(m33**2)
          allfields(:,:,10)= n3/m33

          !! Set kinematic stuff to zero
          if (kin) then
             allfields(:,:,3)=0.
             allfields(:,:,8)=0.
             allfields(:,:,9)=0.
          end if

          !! Compute Reference Fields and Masks
          call get_references(allfields,allrefs,allmask)
          fields(:,:,:,    iz) = allfields
          masky (1:nx,1:ny,iz) = allmask(1:nx,1:ny)
          refs  (:,:,iz  )=allrefs

       end do

    else


       do iz=2,nz-1
          ray    = rays(iz-1:iz+1,:,:,:)
          vslice = Vel(iz,:,:)

          call calc_g11(ray,g11);!   where (g11 .lt. 10e-4) g11=0
          call calc_g12(ray,g12);   where (abs(g12) .lt. 10e-4) g12=0
          call calc_g13(ray,g13);   where (abs(g13) .lt. 10e-4) g13=0
          call calc_g22(ray,g22);!   where (g22 .lt. 10e-4) g22=0
          call calc_g23(ray,g23);   where (abs(g23) .lt. 10e-4) g23=0
          call calc_g33(ray,g33);!   where (g33 .lt. 10e-4) g33=0
          call calc_determinant(g11,g12,g13,g22,g23,g33,gdetr)
          call calc_inverse(g11,g12,g13,g22,g23,g33,iG11,iG12,iG13,iG22,iG23,iG33,gdetr)
          m11=gdetr*iG11;
          m12=gdetr*iG12;
          m13=gdetr*iG13;
          m22=gdetr*iG22;
          m23=gdetr*iG23;
          m33=gdetr*iG33;
          call calc_n(m11,m12,m13,n1)
          call calc_n(m12,m22,m23,n2)
          call calc_n(m13,m23,m33,n3)

          !! Non-stationary coefficient fields
          allfields(:,:,1) = iG13/iG33
          allfields(:,:,2) = iG23/iG33
          allfields(:,:,3) = n3/(2.*m33)
          allfields(:,:,4) = sqrt(1./(vslice**2*iG33))
          allfields(:,:,5) = sqrt(iG11/iG33 - (iG13/iG33)**2)
          allfields(:,:,6) = sqrt(iG22/iG33 - (iG23/iG33)**2)
          allfields(:,:,7) = 2.*iG12/iG33 - 2.*iG13*iG23/(iG33**2)
          allfields(:,:,8) = n1/m33 - m13*n3/(m33**2)
          allfields(:,:,9) = n2/m33 - m23*n3/(m33**2)
          allfields(:,:,10)= n3/m33

          !! Set kinematic stuff to zero
          if (kin) then
             allfields(:,:,3)=0.
             allfields(:,:,8)=0.
             allfields(:,:,9)=0.
          end if

          !! Compute Reference Fields and Masks
          call get_references(allfields,allrefs,allmask)
          fields(:,:,:,    iz) = allfields
          masky (1:nx,1:ny,iz) = allmask(1:nx,1:ny)
          refs  (:,:,iz  )=allrefs

       end do

    end if

    !! Treat ends
    fields(:,:,:,1 )=fields(:,:,:,2); fields(:,:,:,nz)=fields(:,:,:,nz-1)
    refs  (:,:,1   )=refs  (:,:,2  ); refs  (:,:,nz  )=refs  (:,:,nz-1  )
    masky (:,:,1   )=masky (:,:,2  ); masky (:,:,nz  )=masky (:,:,nz-1  )
    mask=masky

    write(0,*) 'Computed reference Parameters and Geometry'
    call coeff_report(fields,refs)

  end subroutine geometry_calc


  !----------------------------------------------------------------

  subroutine calc_g11(ray,g11)
    real,intent(in ) :: ray(:,:,:,:)
    real,intent(out) :: g11(:,:)



    g11(2:nx-1,2:ny-1) = &
    ( ( ray(2,3:nx,2:ny-1,1)-ray(2,1:nx-2,2:ny-1,1) ) / (2.*dx) )**2 + &
    ( ( ray(2,3:nx,2:ny-1,2)-ray(2,1:nx-2,2:ny-1,2) ) / (2.*dx) )**2 + &
    ( ( ray(2,3:nx,2:ny-1,3)-ray(2,1:nx-2,2:ny-1,3) ) / (2.*dx) )**2 
    g11(1 ,:) = g11(   2,:) 
    g11(nx,:) = g11(nx-1,:)
    g11(:,1 ) = g11(:,2   ) 
    g11(:,ny) = g11(:,ny-1)


  end subroutine calc_g11

  !----------------------------------------------------------------

  subroutine calc_g22(ray,g22)
    real,intent(in ) :: ray(:,:,:,:)
    real,intent(out) :: g22(:,:)


    g22(2:nx-1,2:ny-1) = &
    ( ( ray(2,2:nx-1,3:ny,1)-ray(2,2:nx-1,1:ny-2,1) ) / (2.*dy) )**2 + &
    ( ( ray(2,2:nx-1,3:ny,2)-ray(2,2:nx-1,1:ny-2,2) ) / (2.*dy) )**2 + &
    ( ( ray(2,2:nx-1,3:ny,3)-ray(2,2:nx-1,1:ny-2,3) ) / (2.*dy) )**2 

    g22(1 ,:) = g22(   2,:) 
    g22(nx,:) = g22(nx-1,:)
    g22(:,1 ) = g22(:,2   ) 
    g22(:,ny) = g22(:,ny-1)

  end subroutine calc_g22

  !----------------------------------------------------------------

  subroutine calc_g33(ray,g33)
    real,intent(in ) :: ray(:,:,:,:)
    real,intent(out) :: g33(:,:)



       g33(2:nx-1,2:ny-1) = &
       ( ( ray(3,2:nx-1,2:ny-1,1)-ray(1,2:nx-1,2:ny-1,1) ) / (2.*dz) )**2 + &
       ( ( ray(3,2:nx-1,2:ny-1,2)-ray(1,2:nx-1,2:ny-1,2) ) / (2.*dz) )**2 + &
       ( ( ray(3,2:nx-1,2:ny-1,3)-ray(1,2:nx-1,2:ny-1,3) ) / (2.*dz) )**2 

       g33(1 ,:) = g33(   2,:) 
       g33(nx,:) = g33(nx-1,:)
       g33(:,1 ) = g33(:,2   ) 
       g33(:,ny) = g33(:,ny-1)

  end subroutine calc_g33

  !----------------------------------------------------------------

  subroutine calc_g12(ray,g12)
    real,intent(in ) :: ray(:,:,:,:)
    real,intent(out) :: g12(:,:)
  
       g12(2:nx-1,2:ny-1) = &
       ( ray(2,3:nx  ,2:ny-1,1)-ray(2,1:nx-2,2:ny-1,1) ) * &
       ( ray(2,2:nx-1,3:ny  ,1)-ray(2,2:nx-1,1:ny-2,1) ) / (4.*dx*dy) + &
       ( ray(2,3:nx  ,2:ny-1,2)-ray(2,1:nx-2,2:ny-1,2) ) * &
       ( ray(2,2:nx-1,3:ny  ,2)-ray(2,2:nx-1,1:ny-2,2) ) / (4.*dx*dy) + &
       ( ray(2,3:nx  ,2:ny-1,3)-ray(2,1:nx-2,2:ny-1,3) ) * &
       ( ray(2,2:nx-1,3:ny  ,3)-ray(2,2:nx-1,1:ny-2,3) ) / (4.*dx*dy) 

       g12(1 ,:) = g12(   2,:) 
       g12(nx,:) = g12(nx-1,:)
       g12(:,1 ) = g12(:,2   ) 
       g12(:,ny) = g12(:,ny-1)
  end subroutine calc_g12

  !----------------------------------------------------------------

  subroutine calc_g13(ray,g13)
    real,intent(in ) :: ray(:,:,:,:)
    real,intent(out) :: g13(:,:)


       g13(2:nx-1,2:ny-1) = &
       ( ray(2,3:nx  ,2:ny-1,1)-ray(2,1:nx-2,2:ny-1,1) ) * &
       ( ray(3,2:nx-1,2:ny-1,1)-ray(1,2:nx-1,2:ny-1,1) ) / (4.*dx*dz) + &
       ( ray(2,3:nx  ,2:ny-1,2)-ray(2,1:nx-2,2:ny-1,2) ) * &
       ( ray(3,2:nx-1,2:ny-1,2)-ray(1,2:nx-1,2:ny-1,2) ) / (4.*dx*dz) + &
       ( ray(2,3:nx  ,2:ny-1,3)-ray(2,1:nx-2,2:ny-1,3) ) * &
       ( ray(3,2:nx-1,2:ny-1,3)-ray(1,2:nx-1,2:ny-1,3) ) / (4.*dx*dz) 

       g13(1 ,:) = g13(   2,:) 
       g13(nx,:) = g13(nx-1,:)
       g13(:,1 ) = g13(:,2   ) 
       g13(:,ny) = g13(:,ny-1)

  end subroutine calc_g13

  !----------------------------------------------------------------

  subroutine calc_g23(ray,g23)
    real,intent(in ) :: ray(:,:,:,:)
    real,intent(out) :: g23(:,:)

    g23(2:nx-1,2:ny-1) = &
    ( ray(2,2:nx-1,3:ny  ,1)-ray(2,2:nx-1,1:ny-2,1) ) * &
    ( ray(3,2:nx-1,2:ny-1,1)-ray(1,2:nx-1,2:ny-1,1) ) / (4.*dy*dz) + &
    ( ray(2,2:nx-1,3:ny  ,2)-ray(2,2:nx-1,1:ny-2,2) ) * &
    ( ray(3,2:nx-1,2:ny-1,2)-ray(1,2:nx-1,2:ny-1,2) ) / (4.*dy*dz) + &
    ( ray(2,2:nx-1,3:ny  ,3)-ray(2,2:nx-1,1:ny-2,3) ) * &
    ( ray(3,2:nx-1,2:ny-1,3)-ray(1,2:nx-1,2:ny-1,3) ) / (4.*dy*dz)       

    g23(1 ,:) = g23(   2,:) 
    g23(nx,:) = g23(nx-1,:)
    g23(:,1 ) = g23(:,2   ) 
    g23(:,ny) = g23(:,ny-1)

  end subroutine calc_g23

 !----------------------------------------------------------------

  subroutine calc_g11_2d(ray,g11)
    real,intent(in ) :: ray(:,:,:,:)
    real,intent(out) :: g11(:,:)


    g11(2:nx-1,1) = &
    ( ( ray(2,3:nx,1,1)-ray(2,1:nx-2,1,1) ) / (2.*dx) )**2 + &
    ( ( ray(2,3:nx,1,2)-ray(2,1:nx-2,1,2) ) / (2.*dx) )**2 + &
    ( ( ray(2,3:nx,1,3)-ray(2,1:nx-2,1,3) ) / (2.*dx) )**2 

    g11(1 ,:) = g11(   2,1) 
    g11(nx,:) = g11(nx-1,1)



  end subroutine calc_g11_2d

  !----------------------------------------------------------------

  subroutine calc_g22_2d(ray,g22)
    real,intent(in ) :: ray(:,:,:,:)
    real,intent(out) :: g22(:,:)

       g22=1.

  end subroutine calc_g22_2d

  !----------------------------------------------------------------

  subroutine calc_g33_2d(ray,g33)
    real,intent(in ) :: ray(:,:,:,:)
    real,intent(out) :: g33(:,:)

    g33(2:nx-1,1) = &
    ( ( ray(3,2:nx-1,1,1)-ray(1,2:nx-1,1,1) ) / (2.*dz) )**2 + &
    ( ( ray(3,2:nx-1,1,2)-ray(1,2:nx-1,1,2) ) / (2.*dz) )**2 + &
    ( ( ray(3,2:nx-1,1,3)-ray(1,2:nx-1,1,3) ) / (2.*dz) )**2 

    g33(1 ,1) = g33(   2,1) 
    g33(nx,1) = g33(nx-1,1)

  end subroutine calc_g33_2d

  !----------------------------------------------------------------

  subroutine calc_g12_2d(ray,g12)
    real,intent(in ) :: ray(:,:,:,:)
    real,intent(out) :: g12(:,:)
       g12 = 0.
   end subroutine calc_g12_2d

  !----------------------------------------------------------------

   subroutine calc_g13_2d(ray,g13)
     real,intent(in ) :: ray(:,:,:,:)
     real,intent(out) :: g13(:,:)

     g13(2:nx-1,1) = &
     ( ray(2,3:nx  ,1,1)-ray(2,1:nx-2,1,1) ) * &
     ( ray(3,2:nx-1,1,1)-ray(1,2:nx-1,1,1) ) / (4.*dx*dz) + &
     ( ray(2,3:nx  ,1,2)-ray(2,1:nx-2,1,2) ) * &
     ( ray(3,2:nx-1,1,2)-ray(1,2:nx-1,1,2) ) / (4.*dx*dz) + &
     ( ray(2,3:nx  ,1,3)-ray(2,1:nx-2,1,3) ) * &
     ( ray(3,2:nx-1,1,3)-ray(1,2:nx-1,1,3) ) / (4.*dx*dz) 

     g13(1 ,1) = g13(   2,1) 
     g13(nx,1) = g13(nx-1,1)


   end subroutine calc_g13_2d

  !----------------------------------------------------------------

   subroutine calc_g23_2d(ray,g23)
     real,intent(in ) :: ray(:,:,:,:)
     real,intent(out) :: g23(:,:)

     g23 = 0.

   end subroutine calc_g23_2d

  !----------------------------------------------------------------

  subroutine calc_determinant(g11,g12,g13,g22,g23,g33,gdetr)
    real, dimension(:,:),intent(in ) :: g11,g12,g13,g22,g23,g33
    real, dimension(:,:),intent(out) :: gdetr

    gdetr = sqrt( g11*g22*g33-g11*g23**2-g12**2*g33+2.*g12*g23*g13-g22*g13**2 )	

  end subroutine calc_determinant

  !----------------------------------------------------------------

  subroutine calc_inverse(g11,g12,g13,g22,g23,g33,iG11,iG12,iG13,iG22,iG23,iG33,gdetr)
    real,dimension(:,:),intent(in ) :: g11,g12,g13,g22,g23,g33,gdetr
    real,dimension(:,:),intent(out) :: iG11,iG12,iG13,iG22,iG23,iG33

    iG11 = ( g22*g33-g23*g23 ) / gdetr**2
    iG12 = ( g13*g23-g12*g33 ) / gdetr**2
    iG13 = ( g12*g23-g13*g22 ) / gdetr**2
    iG22 = ( g11*g33-g13*g13 ) / gdetr**2
    iG23 = ( g13*g12-g11*g23 ) / gdetr**2
    iG33 = ( g11*g22-g12*g12 ) / gdetr**2

  end subroutine calc_inverse

  !----------------------------------------------------------------

  subroutine calc_n(m11,m12,m13,n)
    real, intent(in) :: m11(:,:),m12(:,:),m13(:,:)
    real, intent(out):: n(:,:)


    n(2:nx-1,2:ny-1) = &
    ( m11(3:nx  ,2:ny-1) - m11(1:nx-2,2:ny-1) )/ (2.*dx) + &
    ( m12(2:nx-1,3:ny  ) - m12(2:nx-1,1:ny-2) )/ (2.*dy) 
    !! Oops ( m13(2:nx-1,2:ny-1) - m13(2:nx-1,2:ny-1) )/ (2.*dz) 

    n(1 ,:) = n(   2,:) 
    n(nx,:) = n(nx-1,:)
    n(:,1 ) = n(:,2   ) 
    n(:,ny) = n(:,ny-1)
  end subroutine calc_n
  !----------------------------------------------------------------

  subroutine calc_n_2d(m11,m12,m13,n)
    real, intent(in) :: m11(:,:),m12(:,:),m13(:,:)
    real, intent(out):: n(:,:)

       n(2:nx-1,1) = &
       ( m11(3:nx  ,1) - m11(1:nx-2,1) )/ (2.*dx)

       n(1 ,1) = n(   2,1) 
       n(nx,1) = n(nx-1,1)

    end subroutine calc_n_2d

  !----------------------------------------------------------------  

  subroutine get_references(allfields,allrefs,allmask)
    integer :: iref,ii
    real    :: allfields(:,:,:)
    real    :: allrefs(:,:)
    integer :: allmask(:,:)

    !! allfields (in)  - Input fields to be used (nx,ny,10)
    !! scales (tmp)    - Scaling function (10)
    !! slice (tmp)     - 1D array used to get refs (nx*ny,10)
    !! iregion (tmp)   - What region has the best attribute set (nx*ny)
    !! nstart (tmp)    - Number of regions (nref)
    !! have_zeros(tmp) - zero mask(10)
    !! allmask(out)    - Output mask (nx,ny)
    !! allrefs (out)   - Reference values in (nref,10)
    !! Global for reference call

    nstart = (/ 1,1,1,1,1,1,1,1,1,1 /); 

    slice=reshape(allfields , (/nx*ny,10/) )

    do ii=1,10 !! Get new scaling function
       scales(ii) = sum( abs( slice(:,ii) ) )  / real(nx*ny)
    end do

    where (scales .eq. 0.) 
       have_zeros = 0.
    elsewhere
       have_zeros = 1. 
    end where

    do ii=1,10 !! Rescale Data
       if ( have_zeros(ii) .ne. 0.) then
          slice(:,ii)=slice(:,ii)/scales(ii)
       else
          slice(:,ii)=0. 
       end if
    end do

    !! Call ND Lloyd's algorithm
    iref=lloyd_go(slice,allrefs,nref,iregion,nstart) 

    do ii=1,10  !! rescaling 
       if ( have_zeros(ii) .ne. 0.) then
          allrefs(:,ii)=allrefs(:,ii)*scales(ii)
       else
          allrefs(:,ii)=0. 
       end if
    end do

    if(iref < nref)  allrefs(iref+1:nref,:)=-1.

    allmask=reshape(iregion,(/nx,ny/) ) !! Mask
  end subroutine get_references

  !----------------------------------------------------------------

  subroutine coeff_report(fields,refs)
    real  ::  fields(:,:,:,:),refs(:,:,:)
    write(0,*) 'METRIC COEFFICIENTS'
    write(0,*) "g11 min/max", minval(g11),maxval(g11)
    write(0,*) "g12 min/max", minval(g12),maxval(g12)
    write(0,*) "g13 min/max", minval(g13),maxval(g13)
    write(0,*) "g22 min/max", minval(g22),maxval(g22)
    write(0,*) "g23 min/max", minval(g23),maxval(g23)
    write(0,*) "g33 min/max", minval(g33),maxval(g33)
    write(0,*) "gdetr min/max", minval(gdetr),maxval(gdetr)
    write(0,*) " "
    write(0,*) 'INVERSE METRIC COEFFICIENTS'
    write(0,*) "G11 min/max", minval(iG11),maxval(iG11)
    write(0,*) "G12 min/max", minval(iG12),maxval(iG12)
    write(0,*) "G13 min/max", minval(iG13),maxval(iG13)
    write(0,*) "G22 min/max", minval(iG22),maxval(iG22)
    write(0,*) "G23 min/max", minval(iG23),maxval(iG23)
    write(0,*) "G33 min/max", minval(iG33),maxval(iG33)
    write(0,*) " "
    write(0,*) "COEFFICIENT RANGES"
    write(0,*) "a1 min/max", minval(fields(:,:,1,:)),maxval(fields(:,:,1,:))
    write(0,*) "a2 min/max", minval(fields(:,:,2,:)),maxval(fields(:,:,2,:))
    write(0,*) "a3 min/max", minval(fields(:,:,3,:)),maxval(fields(:,:,3,:))
    write(0,*) "a4 min/max", minval(fields(:,:,4,:)),maxval(fields(:,:,4,:))
    write(0,*) "a5 min/max", minval(fields(:,:,5,:)),maxval(fields(:,:,5,:))
    write(0,*) "a6 min/max", minval(fields(:,:,6,:)),maxval(fields(:,:,6,:))
    write(0,*) "a7 min/max", minval(fields(:,:,7,:)),maxval(fields(:,:,7,:))
    write(0,*) "a8 min/max", minval(fields(:,:,8,:)),maxval(fields(:,:,8,:))
    write(0,*) "a9 min/max", minval(fields(:,:,9,:)),maxval(fields(:,:,9,:))
    write(0,*) "a10 min/max", minval(fields(:,:,10,:)),maxval(fields(:,:,10,:))
    write(0,*) " "
    write(0,*) "REFERENCE PARAMETERS"
    write(0,*) "b1 min/max", minval(refs(:,1,:)),maxval(refs(:,1,:))
    write(0,*) "b2 min/max", minval(refs(:,2,:)),maxval(refs(:,2,:))
    write(0,*) "b3 min/max", minval(refs(:,3,:)),maxval(refs(:,3,:))
    write(0,*) "b4 min/max", minval(refs(:,4,:)),maxval(refs(:,4,:))
    write(0,*) "b5 min/max", minval(refs(:,5,:)),maxval(refs(:,5,:))
    write(0,*) "b6 min/max", minval(refs(:,6,:)),maxval(refs(:,6,:))
    write(0,*) "b7 min/max", minval(refs(:,7,:)),maxval(refs(:,7,:))
    write(0,*) "b8 min/max", minval(refs(:,8,:)),maxval(refs(:,8,:))
    write(0,*) "b9 min/max", minval(refs(:,9,:)),maxval(refs(:,9,:))
    write(0,*) "b10 min/max", minval(refs(:,10,:)),maxval(refs(:,10,:))

  end subroutine coeff_report

  !----------------------------------------------------------------

  subroutine geometry_close()

    deallocate(g11,g12,g13,g23,g33,iG11,iG12,iG13,iG22,iG23,iG33,m11,m12,m13)
    deallocate(m22,m23,m33,gdetr,n1,n2,n3,allfields,allrefs,allmask)
    deallocate(have_zeros,scales,slice,iregion,nstart)

  end subroutine geometry_close

  !----------------------------------------------------------------  

end module geometry
