module geometry2D_WI
  !! ---------------------------------------------
  !!
  !! Geometry calculation module for RWE2D
  !! Calculates:
  !!		1) Reference Parmeters
  !!		2) Coefficient fields
  !!		3) Mask for most appropriate reference
  !!
  !! Written Jeff Shragge May 16, 2006
  !!			email: jshragge@gmail.com
  !!
  !! ---------------------------------------------
  use SEP
  use util_			  !! Regular SEP Utilities
  use lloyd_mod		  !! Lloyd's module
  !  use sep_timers_mod  !! Timer module for benchmarking

  implicit none
  integer, private :: nz,nx,nnn
  real,    private :: dz,dx
  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  :: gmask,vmask
  real,    allocatable,dimension(:,:), private  :: gfields

  !! For Reference Calculation
  integer, allocatable,private :: iregion(:),nstart(:),tstart(:)
  real,    allocatable,private :: scaly(:),have_zeros(:)
  real,    allocatable,private :: vslice(:),grefs(:,:),slice(:,:),tslice(:,:),trefs(:,:)
contains

  !----------------------------------------------------------------  
  !!
  !! Geometry initialization routine
  !!
  subroutine geometry_init(nz_in,nx_in,dz_in,dx_in,kin_in,nref_in,verbose_in)
    integer          ::    nz_in,nx_in,                   nref_in
    real             ::                dz_in,dx_in
    logical          ::                            kin_in,        verbose_in
    integer          :: nnn
    nz=nz_in;  nx=nx_in
    dz=dz_in;  dx=dx_in
    nref=nref_in;    kin=kin_in
    verbose=verbose_in;    eps=0.0001;

    !! For Calculating Geometry
    nnn=6
    allocate(ray(3,nx,2) )
    allocate(g11(nx),g12(nx),g13(nx),g22(nx),g23(nx),g33(nx))
    allocate(iG11(nx),iG12(nx),iG13(nx),iG22(nx),iG23(nx),iG33(nx))
    allocate(m11(nx),m12(nx),m13(nx),m22(nx),m23(nx),m33(nx))
    allocate(gdetr(nx),n1(nx),n2(nx),n3(nx),vslice(nx) )
    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( gfields(nx,nnn), grefs(nref,nnn), gmask(nx) )
    allocate(                  trefs(nref,2),tstart(2));
    allocate( slice(nx,nnn),tslice(nx,2) )
    gfields=0.; grefs=0.;gmask=0;

    !! For Calculating Reference Parameters
    allocate(have_zeros(nnn));have_zeros=0.;
    allocate(scaly(nnn)); scaly=0.
    allocate(iregion(nx)); iregion=0;
    allocate(nstart(nnn));nstart=0;
    
    call init_lloyd()
    call lookatfields()

  end subroutine geometry_init

  !!----------------------------------------------------------------  
  !! 
  !! Geometry calcuation routine
  !!
  !! Outputs:
  !!		1) fields - coefficient fields
  !!		2) refs	  - references
  !!		3) mask	  - mask for approprite reference
  !!		4) gnorm  - normalization array
  subroutine geometry_calc(rays,Vel,fields,refs,mask,gnorm)
    !! Main subroutine dealing with geometry
    integer :: iz,inn,mask(:,:)!,t1,t2,t3,t4,t5,t6
    !    integer, allocatable :: masky(:,:)
    real    :: rays(:,:,:),Vel(:,:),fields(:,:,:),refs(:,:,:),gnorm(:,:)
    logical          :: logic
    real, allocatable, dimension(:) :: m11old,m13old,m33old
    !    allocate( masky(nx,nz) );    masky=0
    allocate( m11old(nx),m13old(nx),m33old(nx) )
    m11old=0.;m13old=0.;m33old=0.

    !    logic=init_sep_timers()
    !    logic=setup_next_timer("ZL",t1)
    !    logic=setup_next_timer("C1",t2)
    !   logic=setup_next_timer("C2",t3)
    !   logic=setup_next_timer("C3",t4)
    !   logic=setup_next_timer("C4",t5)
    !   logic=setup_next_timer("LL",t6)

    !! 2D Constants
    m12=0.
    m23=0. 
    g12=0.; 
    g22=1.; 
    g23=0.; 
    n2=0.

    iG12 = 0. 
    iG22 = 1. 
    iG23 = 0. 

    do iz=2,nz-1 !! For all extrapolation steps

       !       call start_timer_num(t1)

       !! Metric tensor calculation
       !       call start_timer_num(t2)
       call calc_g11_2d(rays(iz-1:iz+1,:,:),g11);
       call calc_g13_2d(rays(iz-1:iz+1,:,:),g13); 
       if ( maxval(abs(g13)) .lt. 0.01) then
          g13=0.
          !      write(0,*) 'SET g13=0.'
       end if
       call calc_g33_2d(rays(iz-1:iz+1,:,:),g33);
       call calc_determinant(g11,g13,g33,gdetr)
       !       call stop_timer_num(t2)

       !! Inverse metric tensor calculation
       !       call start_timer_num(t3)
       iG11 = g33 / gdetr**2
       iG13 =-g13 / gdetr**2
       iG33 = g11 / gdetr**2

       !! Weighted metric tensor calculation
       m11=gdetr*iG11
       m13=gdetr*iG13
       m22=gdetr
       m33=gdetr*iG33


       !! Stretch/rotate parameter calculation
       call calc_n1(m11,m13,m11old,m13old,n1)
       call calc_n3(m13,m33,m13old,m33old,n3)
       m11old=m11
       m13old=m13
       m33old=m33
       !       call stop_timer_num(t3)

       !       call start_timer_num(t4)
       !! Non-stationary coefficient fields
       gfields(:,1) = 1.*iG13/iG33
       !       if (maxval(abs(gfields(:,1))) .lt. 0.01 ) gfields(:,1)=0.
       gfields(:,2) = n3/(2.*m33)
       gfields(:,3) = sqrt(1./(Vel(iz,:)**2*iG33))
       gfields(:,4) = sqrt(abs(iG11/iG33 - (iG13/iG33)**2))
       gfields(:,5) = n1/m33
       !       if (maxval(abs(gfields(:,5))) .lt. 0.01 ) gfields(:,5)=0.
       gfields(:,6) = n3/m33
       !       call stop_timer_num(t4)

       !       write(0,*) "F1: ",minval(gfields(:,1)),maxval(gfields(:,1))
       !       write(0,*) "F2: ",minval(gfields(:,2)),maxval(gfields(:,2))
       !       write(0,*) "F3: ",minval(gfields(:,3)),maxval(gfields(:,3))
       !       write(0,*) "F4: ",minval(gfields(:,4)),maxval(gfields(:,4))
       !       write(0,*) "F5: ",minval(gfields(:,5)),maxval(gfields(:,5))
       !       write(0,*) "F6: ",minval(gfields(:,6)),maxval(gfields(:,6))



       !! Compute Reference Fields and Masks using N-D Lloyd's algorithm
       !       call start_timer_num(t5)
       call get_references(gfields,grefs,gmask)
       !       call stop_timer_num(t5)

       !! Store all fields in arrays
       !      call start_timer_num(t6)
       gnorm (iz,:   ) = gdetr
       fields(: ,:,iz) = gfields
       mask (: ,  iz) = gmask
       refs  (: ,:,iz) = grefs
       !      call stop_timer_num(t6)

       !     call stop_timer_num(t1)

    end do

    !! Treat ends
    gnorm(1 ,:)=gnorm(2   ,:)
    gnorm(nz,:)=gnorm(nz-1,:)

    fields(:,2,2) =fields(:,2,3)
    fields(:,6,2) =fields(:,6,3)
    fields(:,:,1 )=fields(:,:,2); 
    fields(:,:,nz)=fields(:,:,nz-1)

    refs(:,2,2 )=refs(:,2,3   ) 
    refs(:,6,2 )=refs(:,6,3   ) 
    refs(:,:,1 )=refs(:,:,2   ) 
    refs(:,:,nz)=refs(:,:,nz-1)

    mask (:,1   )=mask (:,2   ) 
    mask (:,nz  )=mask (:,nz-1)


    !    write(0,*) 'Computed reference Parameters and Geometry'
    !    call coeff_report(fields,refs)
    !    call print_timers()
  end subroutine geometry_calc

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

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

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

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

  end subroutine calc_g11_2d

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

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

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

  end subroutine calc_g33_2d

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

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

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

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

  end subroutine calc_g13_2d

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

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

    gdetr = sqrt( abs(g11*g33-g13**2) )	

  end subroutine calc_determinant

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

  subroutine calc_inverse(g11,g13,g33,iG11,iG12,iG13,iG22,iG23,iG33,gdetr)
    real,dimension(:) :: g11,g13,g33,gdetr
    real,dimension(:) :: iG11,iG12,iG13,iG22,iG23,iG33

    iG11 = g33 / gdetr**2
    iG12 = 0. 
    iG13 =-g13 / gdetr**2
    iG22 = 1. 
    iG23 = 0. 
    iG33 = g11 / gdetr**2

  end subroutine calc_inverse

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

  subroutine calc_n3(m13,m33,m13old,m33old,n)
    real, intent(in) :: m13(:),m33(:),m13old(:),m33old(:)
    real, intent(out):: n(:)

    n  = ( m13 - m13old )/ dz + (m33-m33old)/dz

  end subroutine calc_n3
  !----------------------------------------------------------------

  subroutine calc_n1(m11,m13,m11old,m13old,n)
    real, intent(in) :: m11(:),m13(:),m11old(:),m13old(:)
    real, intent(out):: n(:)

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

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

  end subroutine calc_n1

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

  subroutine get_references(gfields,grefs,gmask)
    integer :: iref,ii,ix
    real   :: gfields(:,:)
    real ,intent(out)  :: grefs(:,:)
    integer,intent(out):: gmask(:)

    !! gfields (in)  - Input fields to be used (nx,nnn)
    !! scaly (tmp)    - Scaling function (nnn)
    !! slice (tmp)     - 1D array used to get refs (nx,nnn)
    !! iregion (tmp)   - What region has the best attribute set (nx)
    !! nstart (tmp)    - Number of regions (nref)
    !! have_zeros(tmp) - zero mask(nnn)
    !! gmask(out)    - Output mask (nx)
    !! grefs (out)   - Reference values in (nref,nnn)

    scaly=0.

    scaly(1) = 1.! sum(abs(gfields(:,1)))/dble(nx) 
    scaly(2) = 1. !maxval(gfields(:,2))
    scaly(3) = (maxval(gfields(:,3)))
    scaly(4) = (maxval(gfields(:,4))) 
    scaly(5) = 1.!maxval(gfields(:,5)) 
    scaly(6) = 1.!sum(abs(gfields(:,4)))/dble(nx) 

  !  where (scaly .eq. 0.) 
  !     have_zeros = 0.
  !  elsewhere
  !     have_zeros = 1. 
  !  end where

    slice = gfields

    !if (maxval(abs(gfields(:,1))) .eq. 0.) slice(:,1)=1.
    !if (maxval(abs(gfields(:,2))) .eq. 0.) slice(:,2)=1.
    !if (maxval(abs(gfields(:,5))) .eq. 0.) slice(:,5)=1.
    !if (maxval(abs(gfields(:,6))) .eq. 0.) slice(:,6)=1.
    slice(:,1)=0.
    slice(:,2)=0.
    slice(:,5)=0.
    slice(:,6)=0.

    do ii=3,4 !! Rescale Data
      ! if ( have_zeros(ii) .eq. 1.) then
          slice(:,ii)=slice(:,ii)/scaly(ii)
      ! else
      !    slice(:,ii)=0. 
      ! end if
    end do

    tstart = (/ 4,4 /)
    tslice=slice(:,3:4)
    iref=lloyd_go(tslice,trefs,nref,iregion,tstart) 
    grefs(1:nref,3:4)=trefs
   ! write(0,*) minval(grefs(1:nref,3)),maxval(grefs(1:nref,3))
    !write(0,*) minval(grefs(1:nref,4)),maxval(grefs(1:nref,4))

    do ii=3,4  
   !    if ( have_zeros(ii) .eq. 1.) then
          grefs(1:nref,ii)=(grefs(1:nref,ii))*scaly(ii)
    !   else
    !      grefs(1:nref,ii)=0. 
    !   end if
    end do

   ! write(0,*) minval(grefs(1:nref,3)),maxval(grefs(1:nref,3))
   ! write(0,*) minval(grefs(1:nref,4)),maxval(grefs(1:nref,4))

    !if (maxval(abs(gfields(:,1))) .eq. 0.) grefs(:,1)=0.
    !if (maxval(abs(gfields(:,2))) .eq. 0.) grefs(:,2)=0.
    !if (maxval(abs(gfields(:,5))) .eq. 0.) grefs(:,5)=0.
    !if (maxval(abs(gfields(:,6))) .eq. 0.) grefs(:,6)=0.

    call srite("ref1.H" ,grefs(:,1),size(grefs(:,1))*4)
    call srite("ref2.H" ,grefs(:,2),size(grefs(:,2))*4)
    call srite("ref3.H" ,grefs(:,3),size(grefs(:,3))*4)
    call srite("ref4.H" ,grefs(:,4),size(grefs(:,4))*4)
    call srite("ref5.H" ,grefs(:,5),size(grefs(:,5))*4)
    call srite("ref6.H" ,grefs(:,6),size(grefs(:,6))*4)

    call srite("a1.H" ,gfields(:,1),size(gfields(:,1))*4)
    call srite("a2.H" ,gfields(:,2),size(gfields(:,2))*4)
    call srite("a3.H" ,gfields(:,3),size(gfields(:,3))*4)
    call srite("a4.H" ,gfields(:,4),size(gfields(:,4))*4)
    call srite("a5.H" ,gfields(:,5),size(gfields(:,5))*4)
    call srite("a6.H" ,gfields(:,6),size(gfields(:,6))*4)

    if(iref < nref)  grefs(iref+1:nref,:)=0.

    gmask=reshape(iregion,(/nx/) ) !! Mask

  end subroutine get_references

  !----------------------------------------------------------------  
!!$  subroutine update_velocities(Vel,refs,mask,a3tmp)
!!$    integer :: iref
!!$    real   :: Vel(:,:),vfield(:)
!!$    real   :: vrefs(:),a3tmp(:,:)
!!$    real ,intent(out)  :: refs(:,:,:)
!!$    integer, intent(out):: mask(:,:,:)
!!$
!!$    vstart = (/ nref /)
!!$
!!$    do it=1,nt
!!$       vfield=Vel(:,it)*a3tmp(:,it)
!!$       iref=lloyd_go(vfield,vrefs,nref,vregion,vstart)  
!!$       refs(:,3,it)=vrefs
!!$       if (iref < nref ) refs(iref+1:nref,3,it)=-999.
!!$
!!$       mask(:,3,it)=reshape(vregion(,(/nx/) )
!!$    end do
!!$
!!$    !      refs(:,3,1 )=refs(:,3,2   )
!!$    !      refs(:,3,nt)=refs(:,3,nt-1)
!!$    !      mask(:,3,1 )=mask(:,3,2   )
!!$    !      mask(:,3,nt)=mask(:,3,nt-1)
!!$
!!$  end subroutine update_velocities
!!$  !----------------------------------------------------------------  

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

  subroutine coeff_report(fields,refs)
    real  ::  fields(:,:,:),refs(:,:,:)
    write(0,*) "COEFFICIENT RANGES"
    write(0,*) "a1  min/max", minval(fields(2:,1,:)),maxval(fields(2:,1,:))
    write(0,*) "a3  min/max", minval(fields(2:,2,:)),maxval(fields(2:,2,:))
    write(0,*) "a4  min/max", minval(fields(2:,2,:)),maxval(fields(2:,3,:))
    write(0,*) "a5  min/max", minval(fields(2:,3,:)),maxval(fields(2:,4,:))
    write(0,*) "a8  min/max", minval(fields(2:,3,:)),maxval(fields(2:,5,:))
    write(0,*) "a10 min/max", minval(fields(2:,4,:)),maxval(fields(2:,6,:))
    write(0,*) " "
    write(0,*) "REFERENCE PARAMETERS"
    write(0,*) "b1  min/max", minval(refs(:,1,:)),maxval(refs(:,1,:))
    write(0,*) "b3  min/max", minval(refs(:,1,:)),maxval(refs(:,2,:))
    write(0,*) "b4  min/max", minval(refs(:,2,:)),maxval(refs(:,3,:))
    write(0,*) "b5  min/max", minval(refs(:,3,:)),maxval(refs(:,4,:))
    write(0,*) "b8  min/max", minval(refs(:,3,:)),maxval(refs(:,5,:))
    write(0,*) "b10 min/max", minval(refs(:,4,:)),maxval(refs(:,6,:))

  end subroutine coeff_report
  !----------------------------------------------------------------

  subroutine geom_report( )
    write(0,*) 'METRIC COEFFICIENTS'
    write(0,*) "g11 min/max", minval(g11(2:)),maxval(g11(2:))
    write(0,*) "g13 min/max", minval(g13(2:)),maxval(g13(2:))
    write(0,*) "g33 min/max", minval(g33(2:)),maxval(g33(2:))
    write(0,*) "gdetr min/max", minval(gdetr(2:)),maxval(gdetr(2:))
    write(0,*) " "
    write(0,*) 'INVERSE METRIC COEFFICIENTS'
    write(0,*) "G11 min/max", minval(iG11(2:)),maxval(iG11(2:))
    write(0,*) "G13 min/max", minval(iG13(2:)),maxval(iG13(2:))
    write(0,*) "G22 min/max", minval(iG22(2:)),maxval(iG22(2:))
    write(0,*) "G33 min/max", minval(iG33(2:)),maxval(iG33(2:))
    write(0,*) " "

  end subroutine geom_report

  subroutine lookatfields()
    call auxputch("n1","i",nref,"ref1.H" )
    call auxputch("n1","i",nref,"ref2.H" )
    call auxputch("n1","i",nref,"ref3.H" )
    call auxputch("n1","i",nref,"ref4.H" )
    call auxputch("n1","i",nref,"ref5.H" )
    call auxputch("n1","i",nref,"ref6.H" )

    call auxputch("n2","i",nz-2,"ref1.H" )
    call auxputch("n2","i",nz-2,"ref2.H" )
    call auxputch("n2","i",nz-2,"ref3.H" )
    call auxputch("n2","i",nz-2,"ref4.H" )
    call auxputch("n2","i",nz-2,"ref5.H" )
    call auxputch("n2","i",nz-2,"ref6.H" )

    call auxputch("o1","r",0.,"ref1.H" )
    call auxputch("o1","r",0.,"ref2.H" )
    call auxputch("o1","r",0.,"ref3.H" )
    call auxputch("o1","r",0.,"ref4.H" )
    call auxputch("o1","r",0.,"ref5.H" )
    call auxputch("o1","r",0.,"ref6.H" )

    call auxputch("o2","r",0.,"ref1.H" )
    call auxputch("o2","r",0.,"ref2.H" )
    call auxputch("o2","r",0.,"ref3.H" )
    call auxputch("o2","r",0.,"ref4.H" )
    call auxputch("o2","r",0.,"ref5.H" )
    call auxputch("o2","r",0.,"ref6.H" )

    call auxputch("d1","r",1.,"ref1.H" )
    call auxputch("d1","r",1.,"ref2.H" )
    call auxputch("d1","r",1.,"ref3.H" )
    call auxputch("d1","r",1.,"ref4.H" )
    call auxputch("d1","r",1.,"ref5.H" )
    call auxputch("d1","r",1.,"ref6.H" )

    call auxputch("d2","r",1.,"ref1.H" )
    call auxputch("d2","r",1.,"ref2.H" )
    call auxputch("d2","r",1.,"ref3.H" )
    call auxputch("d2","r",1.,"ref4.H" )
    call auxputch("d2","r",1.,"ref5.H" )
    call auxputch("d2","r",1.,"ref6.H" )


    call auxputch("n1","i",nx,"a1.H" )
    call auxputch("n1","i",nx,"a2.H" )
    call auxputch("n1","i",nx,"a3.H" )
    call auxputch("n1","i",nx,"a4.H" )
    call auxputch("n1","i",nx,"a5.H" )
    call auxputch("n1","i",nx,"a6.H" )

    call auxputch("n2","i",nz-2,"a1.H" )
    call auxputch("n2","i",nz-2,"a2.H" )
    call auxputch("n2","i",nz-2,"a3.H" )
    call auxputch("n2","i",nz-2,"a4.H" )
    call auxputch("n2","i",nz-2,"a5.H" )
    call auxputch("n2","i",nz-2,"a6.H" )

    call auxputch("o1","r",0.,"a1.H" )
    call auxputch("o1","r",0.,"a2.H" )
    call auxputch("o1","r",0.,"a3.H" )
    call auxputch("o1","r",0.,"a4.H" )
    call auxputch("o1","r",0.,"a5.H" )
    call auxputch("o1","r",0.,"a6.H" )

    call auxputch("o2","r",0.,"a1.H" )
    call auxputch("o2","r",0.,"a2.H" )
    call auxputch("o2","r",0.,"a3.H" )
    call auxputch("o2","r",0.,"a4.H" )
    call auxputch("o2","r",0.,"a5.H" )
    call auxputch("o2","r",0.,"a6.H" )

    call auxputch("d1","r",1.,"a1.H" )
    call auxputch("d1","r",1.,"a2.H" )
    call auxputch("d1","r",1.,"a3.H" )
    call auxputch("d1","r",1.,"a4.H" )
    call auxputch("d1","r",1.,"a5.H" )
    call auxputch("d1","r",1.,"a6.H" )

    call auxputch("d2","r",1.,"a1.H" )
    call auxputch("d2","r",1.,"a2.H" )
    call auxputch("d2","r",1.,"a3.H" )
    call auxputch("d2","r",1.,"a4.H" )
    call auxputch("d2","r",1.,"a5.H" )
    call auxputch("d2","r",1.,"a6.H" )

  end subroutine lookatfields
  !----------------------------------------------------------------

  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,gfields,grefs,gmask)
    ! deallocate(have_zeros,scaly,slice,iregion,nstart,vslice)

  end subroutine geometry_close



end module geometry2D_WI
