module RWE2D_polar_wemig
  !! Migration driver for the RWE2D migration program
  !! contains
  !!		1) wemig_init		- wemig initialization routine
  !!		2) wemig		- main migration driver routine
  !!		3) wemig_ssfboth	- split-step fourier propagation
  !!		4) wemig_phsboth	- phase-shift propagation
  !!		5) wemig_refboth	- interpolation for various phase-shifts
  !!		6) rwetaper_init	- a cosine-taper initialization routine
  !!		7) rwetaper		- apply cosine-taper to reduce boundary artifacts
  use sep !! basic SEPlib routine I/O module
  use util_ !! handle axes module
  use fft !! FFT module - could use others though
  use lloyd_mod !! Lloyd's code for velocity picking
 !! use omppar_
  !! use sep_timers_mod

  implicit none
  integer, private :: nr,nt,nw,nloop,nref,ntap
  real,    private :: dr,dt,dw,or,ot,ow,okr,dkr,forward,eps,pi,okt,dkt,r0
  logical, private :: kin,verbose

  !! Geometry
  real,    allocatable,private :: rline(:),oline(:)


  !! Global references
  integer, allocatable,private :: rwf_allmask(:,:  ),swf_allmask(:,:  )
  real,    allocatable,private :: rwf_allrefs(:,:,:),rwf_allflds(:,:,:)
  real,    allocatable,private :: swf_allrefs(:,:,:),swf_allflds(:,:,:)

  !! Temp in/out to set up global refs
  integer, allocatable,private :: rwf_mask(:  ),swf_mask(:  )
  real,    allocatable,private :: rwf_refs(:,:),rwf_flds(:,:)
  real,    allocatable,private :: swf_refs(:,:),swf_flds(:,:)
  integer, allocatable,private :: rmask(:  ),smask(:  )
  real,    allocatable,private :: rrefs(:,:),rflds(:,:)
  real,    allocatable,private :: srefs(:,:),sflds(:,:)

  !! Inside lloyds
  integer, allocatable,private :: iregion_rwf(:),iregion_swf(:),nstart(:)
  real,    allocatable,private :: slice_rwf(:,:),slice_swf(:,:),scaly(:)

  !! Propagation
  real,    allocatable,private :: rmsk(:),smsk(:),rmtt(:),smtt(:),rtap(:),stap(:)
  complex, allocatable,private :: swfl(:,:),swtt(:,:),rwfl(:,:),rwtt(:,:),rdax(:,:),sdax(:,:)
  complex, allocatable,private :: swf_save(:,:,:),rwf_save(:,:,:)

  !! Temp. image
  real,    allocatable,private :: Timg(:,:,:)

contains

  !----------------------------------------------------------------
  !! Initialization routine - parameter passing and array allocation
  subroutine wemig_init(ar,at,aw,forward_in,kin_in,nref_in,verbose_in,r0_in)
    type(myaxis)     :: ar,at,aw
    integer          ::                            nref_in
    real             ::          forward_in,                          r0_in
    logical          ::                     kin_in,        verbose_in
    real             :: pi

  !! For OMP
    integer, external :: omp_get_num_threads
    integer, external :: omp_get_thread_num
    integer  :: nth,ith,i1

    pi=acos(-1.)
    forward=forward_in
    nr=ar%n;    nt=at%n;
    dr=ar%d;    dt=at%d; 
    or=ar%o;    ot=at%o; 
    nw=aw%n;

    dw=aw%d * 2.*pi; ow=aw%o * 2.*pi
    verbose=verbose_in;    
    r0=r0_in
    nref=nref_in

    call from_param("ntap",ntap,5) !! Width of cosine taper
    call from_param("nloop",nloop,3) !! Number of smoothing iterations 
    okr= -pi/dr;  dkr=2.*pi/(nr*dr)
    okt= -pi/dt;  dkt=2.*pi/(nt*dt)
    kin=kin_in
    eps=0.0000001

    !$OMP PARALLEL
    nth = omp_get_num_threads()
    !$OMP END PARALLEL

    !! Global Reference allocation
    allocate(rwf_allmask(nr,    nt),swf_allmask(nt    ,nr))
    allocate(rwf_allflds(nr,2  ,nt),swf_allflds(nt,2  ,nr)) 
    allocate(rwf_allrefs(nref,2,nt),swf_allrefs(nref,2,nr))

    !! Temp in/out to set up global refs
    allocate(rwf_mask(nr),swf_mask(nt),rwf_flds(nr,2),swf_flds(nt,2) ) 
    allocate(rwf_refs(nref,2),swf_refs(nref,2),rwf_save(nr,nt,nth),swf_save(nr,nt,nth) )
    allocate(rmask(nr),smask(nt),rflds(nr,2),sflds(nt,2),rrefs(nref,2),srefs(nref,2))

    !! Inside lloyds
    allocate( scaly(2),iregion_rwf(nr),iregion_swf(nt))
    allocate(nstart(2),slice_rwf(nr,2),slice_swf(nt,2))

    !! Propgation
    allocate( rdax(nr,nth), sdax(nt,nth), rwfl(nr,nth), swfl(nt,nth), rwtt(nr,nth), swtt(nt,nth))
    allocate( rmsk(nr), smsk(nt), rmtt(nr), smtt(nt), rtap(nr), stap(nt))

    !! Geometry
    allocate(rline(nr),oline(nt))
    rline=(/ (i1,i1=0,nr-1) /)*dr+r0
    oline=(/ (i1,i1=0,nt-1) /)*dt+ot
    rtap=1.;    stap=1.
    rmsk=0.;    smsk=0.;
    call rwetaper_rwf_init() !! Initialize the taper program
    call rwetaper_swf_init() !! Initialize the taper program

    !! Temp Image
    allocate( Timg(nr,nt,nth) )
    Timg=0.

    if (verbose) then
       call swf_print()
       call rwf_print()
       call img_print()
       call ref1_print()
       call ref2_print()
       call ref3_print()
       call ref4_print()
       call a_print()
    end if


  end subroutine wemig_init

  !----------------------------------------------------------------
  !! WE migration driver
  !!
  !! Basic routine is:
  !!
  !!	for all frequencies
  !!		for all depths
  !!			Split-step Fourier Propagation
  !!			FFT wavefields
  !!			for all References
  !!				Phase-shift propagation
  !!				IFFT wavefields
  !!				Wavefield Interpolation
  !!			end References
  !!			Imaging Condition
  !!		end depths
  !!	end frequencies
  !!
  subroutine wemig(rwf,swf,Pimg,ishot,vel)
    integer            :: iw,ir,it,inn,is,id,ip,ishot
    real               :: rarg,saxmax,sarg,ss,w,Pimg(:,:),vel(:,:)
    complex            :: rwf(:,:),swf(:,:)!!,rax(:),sax(:)
    complex,allocatable:: rax(:,:),sax(:,:) !! rwf(:,:),swf(:,:),
    !! For OMP
    integer, external :: omp_get_num_threads
    integer, external :: omp_get_thread_num
    integer  :: nth,ith

    !$OMP PARALLEL
    nth = omp_get_num_threads()
    !$OMP END PARALLEL


    allocate( sax(nt,nth),rax(nr,nth)) 

    !! set and freeze the # of threads
    call omp_set_dynamic(.false.)
    call omp_set_num_threads(nth)

    !! . . Forward Scattering Option
    sarg=-1.; !! if ( forward .ne. 0. ) sarg=1.
    rarg= 1.

    !! Calculate Geometry for swf and rwf
    call init_lloyd()

    write(0,*) 'CALCULATING SOURCE REFERENCES'
    !! Source Geometry

    !$OMP PARALLEL PRIVATE(id)
    id = omp_get_thread_num()
    write(0,*) 'HELLO WORLD from thread ', id
    !$OMP barrier
    if (id .eq. 0) then
       nth = omp_get_num_threads()
       write(0,*) 'THERE ARE ',nth,' threads'
    end if
    !$OMP END PARALLEL


    do ir=1,nr
       !where (vel(ir,:)==0 ) 
       !   vel(ir,:)=maxval(vel(ir,:))
       !end where
       do it=1,nt
          swf_flds(it,1)=1./vel(ir,it)
          swf_flds(it,2)=1./rline(ir)
       end do

       if (modulo(8*ir,nr) .eq. 0) write(0,*) "S Ref",ir,nr
       call get_references_swf(swf_flds,swf_refs,swf_mask)

       do it=1,nt
          swf_allflds(it,1,ir)=swf_flds(it,1)
          swf_allflds(it,2,ir)=swf_flds(it,2)
          swf_allmask(it,  ir)=swf_mask(it  )
       end do

       do is=1,nref
          swf_allrefs(is,1,ir)=swf_refs(is,1)
          swf_allrefs(is,2,ir)=swf_refs(is,2)
       end do
       
       swf_refs=0.; swf_flds=0.; swf_mask=1
       
    end do

    do it=1,nt

       do ir=1,nr
          rwf_flds(ir,1)=rline(ir)/vel(ir,it)
          rwf_flds(ir,2)=rline(ir)
       end do

       if (modulo(8*it,nt) .eq. 0) write(0,*) "R Ref",it,nt

       call get_references_rwf(rwf_flds,rwf_refs,rwf_mask)

       do ir=1,nr
          rwf_allflds(ir,1,it)=rwf_flds(ir,1)
          rwf_allflds(ir,2,it)=rwf_flds(ir,2)
          rwf_allmask(ir,  it)=rwf_mask(ir  )
       end do

       do is=1,nref
          rwf_allrefs(is,1,it) = rwf_refs(is,1)
          rwf_allrefs(is,2,it) = rwf_refs(is,2)
       end do

        rwf_refs=0.; rwf_flds=0.; rwf_mask=1
   end do


    write(0,*) 'STARTING MIGRATION'

    !$OMP PARALLEL DO SCHEDULE(DYNAMIC) PRIVATE(iw,w,id,it,ir,inn,stap,rtap) SHARED(swf_allflds,swf_allmask,swf_allrefs,rwf_allflds,rwf_allmask,rwf_allrefs,swf_save,rwf_save,Timg,nw,ow,dw,nt,nr,nref,rarg,sarg,sax,rax,swtt,rwtt,swfl,rwfl)
    do iw=1,nw  !! FREQUENCY LOOP         
       id = omp_get_thread_num()+1
       w=ow+(iw-1)*dw !! Assign frequency

       do it=1,nt
          sax(it,id) = swf(it,iw)*sin(oline(it)) 
       end do
       saxmax=maxval(cabs(sax(:,id)))

       do it=1,nt
          swf_save(1,it,id)=sax(it,id)   
       end do

       call rwetaper_swf(sax,id) 
       do ir=2,nr  !! RADIUS LOOP

          call wemig_ssf_swf(id,sax,w,swf_allmask(:,ir),swf_allflds(:,:,ir),swf_allrefs(:,:,ir),sarg)

          do it=1,nt
             swtt(it,id)=sax(it,id);
             sax(it,id)=0.
          end do

          !$OMP CRITICAL
          call fth(.false.,.false.,swtt(:,id) )
          !$OMP END CRITICAL

          do inn=1,nref
             if (swf_allrefs(inn,1,ir) .eq. -1.) cycle 
             do it=1,nt
                swfl(it,id)=swtt(it,id)
             end do
             call wemig_phs_swf(id,swfl,w,inn,swf_allrefs(:,:,ir),sarg,ir)
             !$OMP CRITICAL
             call fth( .true.,.false.,swfl(:,id) )
             !$OMP END CRITICAL
             call wemig_ref_swf(id,swfl,ir,inn,swf_allmask,sax)
          end do

          call rwetaper_swf(sax,id) 

          ss = sqrt(maxval(abs((sax*conjg(sax))))+0.01*saxmax) 
          do it=1,nt 
             swf_save(ir,it,id) = sax(it,id)!/ss
          end do

       end do
       !!write(0,*) 'FINAL SWF_SAVE ',maxval(abs(swf_save(nr,:,id))),iw,id

       !! RECEIVER WAVEFIELD
       do ir=1,nr
          rax(ir,id) = rwf(ir,iw) !! Choose receiver wavefield frequency
       end do

       do ir=1,nr
          rwf_save(ir,1,id)=rwf(ir,iw) 
       end do

       call rwetaper_rwf(rax,id) 
       do it=2,nt !! THETA LOOP

          call wemig_ssf_rwf(id,rax,w,rwf_allmask(:,it),rwf_allflds(:,:,it),rwf_allrefs(:,:,it),rarg)

          do ir=1,nr
             rwtt(ir,id)=rax(ir,id);
             rax(ir,id)=0. 
          end do

          !$OMP CRITICAL
          call fth(.false.,.false.,rwtt(:,id) )
          !$OMP END CRITICAL

          do inn=1,nref  !! REFERENCE VELOCITY LOOP
             if (rwf_allrefs(inn,2,it) .eq. -1.) cycle
             do ir=1,nr
                rwfl(ir,id)=rwtt(ir,id)
             end do
             call wemig_phs_rwf(id,rwfl,w,inn,rwf_allrefs(:,:,it),rarg,it)
             !$OMP CRITICAL
             call fth( .true.,.false.,rwfl(:,id) )
             !$OMP END CRITICAL
             call wemig_ref_rwf(id,rwfl,inn,it,rwf_allmask,rax)
          end do !! END Reference loop

          call rwetaper_rwf(rax,id) 
 
          do ir=1,nr
             rwf_save(ir,it,id)=rax(ir,id)
          end do

          !! Imaging condition
          !$OMP CRITICAL
          do ir=1,nr
             Timg(ir,it,id)=Timg(ir,it,id)+&
             real(conjg(rax(ir,id))*swf_save(ir,it,id))
          end do
          !$OMP END CRITICAL

       end do !! END Depth loop

     !!  write(0,*) 'FINAL RWF_SAVE ',maxval(abs(rwf_save(:,nt,id))), iw,id,ishot

       if (modulo(iw,nw/10) .eq. 0 ) then
          write(0,*) "SSF",iw,nw,ishot,maxval(abs(Timg(:,:,id)))
       end if

       ! if (verbose) then
       !    call srite("swf.H",swf_save(:,:,ip),size(swf_save)*8/nth)
       !    call srite("rwf.H",rwf_save(:,:,ip),size(rwf_save)*8/nth)
       !    call srite("img.H",Timg(:,:,ip),size(Timg)*4/nth)
       ! end if

      ! do it=1,nt
      !    do ir=1,nr
      !       rwf_save(ir,it,id)=0.
      !       swf_save(ir,it,id)=0.
      !    end do
      ! end do
    end do  
    !$OMP END PARALLEL DO

    !! SUMMATION OVER OUTSIDE IMAGE CONDITION
    do it=1,nt
       do ir=1,nr
          do ith=1,nth
             Pimg(ir,it)=Pimg(ir,it)+Timg(ir,it,ith)
             Timg(ir,it,ith)=0.
          end do
       end do
    end do
    
    
  end subroutine wemig

  !----------------------------------------------------------------
  !!
  !! Split-step Fourier Propagation subroutine 
  !!	
  subroutine wemig_ssf_rwf(id,rdax,w,rwf_mask,rwf_flds,rwf_refs,rarg) 
    integer :: ir,rwf_mask(:),id
    real    :: w,rarg,rwf_flds(:,:),rwf_refs(:,:),iikt2d
    complex :: rdax(:,:)

    !! Calculate w-r domain correction step
    do ir=1,nr
       iikt2d   =  -dt*w*(rwf_flds(ir,1)-rwf_refs(rwf_mask(ir),1))
       rdax(ir,id) =  rdax(ir,id) * cmplx( cos(iikt2d) ,rarg*sin(iikt2d) )
    end do
  end subroutine wemig_ssf_rwf

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

  subroutine wemig_ssf_swf(id,sdax,w,swf_mask,swf_flds,swf_refs,sarg)
    integer :: it,swf_mask(:),id
    real    :: w,sarg,swf_flds(:,:),swf_refs(:,:),iikr2d
    complex :: sdax(:,:)

    !! Calculate w-theta domain correction step
    do it=1,nt
       iikr2d = -dr*w*( swf_flds(it,1) - swf_refs(swf_mask(it),1) )
       sdax(it,id) = sdax(it,id)*cmplx( cos(iikr2d) ,sarg*sin(iikr2d) )
    end do

  end subroutine wemig_ssf_swf

  !----------------------------------------------------------------  
  !!
  !! Fourier-domain phase shift
  !!
  subroutine wemig_phs_rwf(id,rdax,w,inn,rwf_refs,rarg,it)
    integer :: ir,inn,it,id
    real    :: w,kr,rarg,carg,carg1,rwf_refs(:,:),pkr
    complex :: rdax(:,:)

    carg1= (w*rwf_refs(inn,1))**2
    do ir=1,nr 
       kr=okr+(ir-1)*dkr  
       carg = carg1 - (kr*rwf_refs(inn,2))**2  
       if ( carg .lt. 0) then 
          pkr = -dt*abs(sqrt(-carg))  
          rdax(ir,id)=rdax(ir,id)*exp(pkr) 
       else 
          pkr = -dt*sqrt(carg) 
          rdax(ir,id) = rdax(ir,id)*cmplx(cos(pkr),rarg*sin(pkr)) 
       end if
    end do  

  end subroutine wemig_phs_rwf

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

  subroutine wemig_phs_swf(id,sdax,w,inn,swf_refs,sarg,ir)
    integer :: it,inn,ir,id
    real    :: w,kt,sarg,carg,carg1,swf_refs(:,:),pkt
    complex :: sdax(:,:)

    carg1= (w*swf_refs(inn,1))**2 
    do it=1,nt  
       kt=okt+(it-1)*dkt    
       carg = carg1 - (kt*swf_refs(inn,2))**2  
       if ( carg .lt. 0) then        
          pkt = -dr*abs(sqrt(-carg)) 
          sdax(it,id)=sdax(it,id)*exp(pkt) 
       else 
          pkt =-dr*sqrt(carg) 
          sdax(it,id)=sdax(it,id)*cmplx(cos(pkt),sarg*sin(pkt))
       end if
    end do 

  end subroutine wemig_phs_swf

  !!----------------------------------------------------------------
  !!
  !! Interpolation routine for different phase-shifts
  !!
  subroutine wemig_ref_rwf(id,rwfl,inn,it,rwf_allmask,rdax)
    integer :: rwf_allmask(:,:),id
    complex :: rdax(:,:),rwfl(:,:)
    integer :: ir,inn,ii,it

    !! Find out locations where want to use wavefield for given
    !! Reference values
    rmtt=0.
    where( rwf_allmask(:,it) .eq. inn) 
       rmtt=1.
    end where

    do ii=1,nloop
       rmsk(1 ) = rmtt(1   )
       rmsk(nr) = rmtt(nr-1)
       do ir=2,nr-1
          rmsk(ir)=(rmtt(ir-1)+rmtt(ir+1) )/2.
       end do
       do ir=1,nr
          rmtt(ir)=rmsk(ir)
       end do
    end do

    do ir=1,nr
       rdax(ir,id) = rdax(ir,id) + rmsk(ir)*rwfl(ir,id)
    end do
  end subroutine wemig_ref_rwf

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

  subroutine wemig_ref_swf(id,swfl,ir,inn,swf_allmask,sdax)
    integer :: swf_allmask(:,:),id
    complex :: sdax(:,:),swfl(:,:)
    integer :: it,inn,ii,ir

    smtt=0.
    where( swf_allmask(:,ir) .eq. inn) 
       smtt=1.
    end where

    do ii=1,nloop
       smsk(1 ) = smtt(1   )
       smsk(nt) = smtt(nt-1)
       do it=2,nt-1
          smsk(it)=( smtt(it-1)+smtt(it+1) )/2.
       end do
       do it=1,nt
          smtt(it)=smsk(it)
       end do
    end do

    do it=1,nt
       sdax(it,id) = sdax(it,id) + smsk(it)*swfl(it,id)
    end do

  end subroutine wemig_ref_swf

  !----------------------------------------------------------------
  !! Taper Initialization Program
  subroutine rwetaper_rwf_init()
    integer :: itap,jtap
    real :: pi
    pi=acos(-1.)

    if (ntap .gt. 1) then

       do itap=1,ntap
          jtap = abs(ntap-itap)
          rtap(   itap) = cos(pi/2.* jtap/(ntap-1))
       end do

       do itap=1,ntap
          jtap = abs(ntap-itap)
          rtap(nr-itap+1) = cos(pi/2.* jtap/(ntap-1))
       end do

    end if

  end subroutine rwetaper_rwf_init
  !!----------------------------------------------------------------
  subroutine rwetaper_swf_init()
    integer :: itap,jtap
    real :: pi
    pi=acos(-1.)

    if (ntap .gt. 1) then

       do itap=1,ntap
          jtap = abs(ntap-itap)
          stap(   itap) = cos(pi/2.* jtap/(ntap-1))
       end do

       do itap=1,ntap
          jtap = abs(ntap-itap)
          stap(nt-itap+1) = cos(pi/2.* jtap/(ntap-1))
       end do

    end if

  end subroutine rwetaper_swf_init

  !!----------------------------------------------------------------
  !!
  !! Apply Taper to wavefield
  !!
  subroutine rwetaper_rwf(dax,id)
    complex  :: dax(:,:)
    integer :: ip,id
    do ip=1,nr
       dax(ip,id) = dax(ip,id) * rtap(ip)
    end do
  end subroutine rwetaper_rwf

  !----------------------------------------------------------------
  subroutine rwetaper_swf(dax,id)
    complex  :: dax(:,:)
    integer :: ip,id
    do ip=1,nt
       dax(ip,id) = dax(ip,id) * stap(ip)
    end do
 end subroutine rwetaper_swf

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

  subroutine get_references_rwf(rwf_flds,rwf_refs,rwf_mask)
    integer ::iref,ii,ix,ir,id,icnt
    real   :: s1,s2
    real   :: rwf_flds(:,:)
    real   :: rwf_refs(:,:)
    integer:: rwf_mask(:)

    nstart = (/ 6,6 /)

    s1 = sum(rwf_flds(:,1))/dble(nr) 
    s2 = sum(rwf_flds(:,2))/dble(nr) 

    do ir=1,nr
       slice_rwf(ir,1)=rwf_flds(ir,1)/s1
       slice_rwf(ir,2)=rwf_flds(ir,2)/s2
    end do

    iref=lloyd_go(slice_rwf,rwf_refs,nref,iregion_rwf,nstart) 

    do ir=1,nref
       rwf_refs(ir,1)=rwf_refs(ir,1)*s1
       rwf_refs(ir,2)=rwf_refs(ir,2)*s2
    end do

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

    do ir=1,nr
       rwf_mask(ir)=iregion_rwf(ir) 
    end do

    if (verbose) then
       call srite("ref1.H" ,rwf_refs(:,1),size(rwf_refs(:,1))*4)
       call srite("ref2.H" ,rwf_refs(:,2),size(rwf_refs(:,2))*4)
       call srite("a1.H" ,rwf_flds(:,1),size(rwf_flds(:,1))*4)
       call srite("a2.H" ,rwf_flds(:,2),size(rwf_flds(:,2))*4)
    end if

  end subroutine get_references_rwf
  !----------------------------------------------------------------
  subroutine get_references_swf(swf_flds,swf_refs,swf_mask)
    integer :: iref,ii,ix,it,id,icnt
    real    :: s1,s2
    real    :: swf_flds(:,:)
    real    :: swf_refs(:,:)
    integer :: swf_mask(:)

    nstart = (/ 6,6 /)

    s1 = sum(swf_flds(:,1))/dble(nt) 
    s2 = sum(swf_flds(:,2))/dble(nt) 

    do it=1,nt
       slice_swf(it,1) = swf_flds(it,1)/s1
       slice_swf(it,2) = swf_flds(it,2)/s2
    end do

    iref=lloyd_go(slice_swf,swf_refs,nref,iregion_swf,nstart) 

    do ix=1,nref
       swf_refs(ix,1)=swf_refs(ix,1)*s1
       swf_refs(ix,2)=swf_refs(ix,2)*s2
    end do

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

    do it=1,nt
       swf_mask(it)=iregion_swf(it)
    end do

    if (verbose) then
       call srite("ref3.H" ,swf_refs(:,1),size(swf_refs(:,1))*4)
       call srite("ref4.H" ,swf_refs(:,2),size(swf_refs(:,2))*4)
       call srite("a3.H"   ,swf_flds(:,1),size(swf_flds(:,1))*4)
       call srite("a4.H"   ,swf_flds(:,2),size(swf_flds(:,2))*4)
    end if

  
  end subroutine get_references_swf
  !----------------------------------------------------------------

  subroutine swf_print()
    call auxputch("n1","i",nr,"swf.H")
    call auxputch("n2","i",nt,"swf.H")
    call auxputch("d1","r",dr,"swf.H")
    call auxputch("d2","r",dt,"swf.H")
    call auxputch("o1","r",or,"swf.H")
    call auxputch("o2","r",ot,"swf.H")
    call auxputch("esize","i",8,"swf.H")
    call auxputch("n3","i",nw,"swf.H")
    call auxputch("d3","r",dw,"swf.H")
    call auxputch("o3","r",ow,"swf.H")
  end subroutine swf_print
  !----------------------------------------------------------------

  subroutine rwf_print()
    call auxputch("n2","i",nr,"rwf.H")
    call auxputch("n1","i",nt,"rwf.H")
    call auxputch("d2","r",dr,"rwf.H")
    call auxputch("d1","r",dt,"rwf.H")
    call auxputch("o2","r",or,"rwf.H")
    call auxputch("o1","r",ot,"rwf.H")
    call auxputch("esize","i",8,"rwf.H")
    call auxputch("n3","i",nw,"rwf.H")
    call auxputch("d3","r",dw,"rwf.H")
    call auxputch("o3","r",ow,"rwf.H")
  end subroutine rwf_print
  !----------------------------------------------------------------

  subroutine img_print()
    call auxputch("n1","i",nr,"img.H")
    call auxputch("n2","i",nt,"img.H")
    call auxputch("d1","r",dr,"img.H")
    call auxputch("d2","r",dt,"img.H")
    call auxputch("o1","r",or,"img.H")
    call auxputch("o2","r",ot,"img.H")
    call auxputch("esize","i",4,"img.H")
    call auxputch("n3","i",nw,"img.H")
    call auxputch("d3","r",dw,"img.H")
    call auxputch("o3","r",ow,"img.H")
  end subroutine img_print

  subroutine ref1_print()
    call auxputch("n1","i",nr  ,"ref1.H")
    call auxputch("n2","i",nref,"ref1.H")
    call auxputch("d1","r",dr  ,"ref1.H")
    call auxputch("d2","r",dt  ,"ref1.H")
    call auxputch("o1","r",or  ,"ref1.H")
    call auxputch("o2","r",ot  ,"ref1.H")
    call auxputch("esize","i",4,"ref1.H")
  end subroutine ref1_print

  subroutine ref2_print()
    call auxputch("n1","i",nr  ,"ref2.H")
    call auxputch("n2","i",nref,"ref2.H")
    call auxputch("d1","r",dr  ,"ref2.H")
    call auxputch("d2","r",dt  ,"ref2.H")
    call auxputch("o1","r",or  ,"ref2.H")
    call auxputch("o2","r",ot  ,"ref2.H")
    call auxputch("esize","i",4,"ref2.H")
  end subroutine ref2_print

  subroutine ref3_print()
    call auxputch("n1","i",nt  ,"ref3.H")
    call auxputch("n2","i",nref,"ref3.H")
    call auxputch("d1","r",dr  ,"ref3.H")
    call auxputch("d2","r",dt  ,"ref3.H")
    call auxputch("o1","r",or  ,"ref3.H")
    call auxputch("o2","r",ot  ,"ref3.H")
    call auxputch("esize","i",4,"ref3.H")
  end subroutine ref3_print


  subroutine ref4_print()
    call auxputch("n1","i",nr  ,"ref4.H")
    call auxputch("n2","i",nref,"ref4.H")
    call auxputch("d1","r",dr  ,"ref4.H")
    call auxputch("d2","r",dt  ,"ref4.H")
    call auxputch("o1","r",or  ,"ref4.H")
    call auxputch("o2","r",ot  ,"ref4.H")
    call auxputch("esize","i",4,"ref4.H")
  end subroutine ref4_print

  subroutine a_print()
    call auxputch("n1","i",nr,"a1.H" )
    call auxputch("n1","i",nr,"a2.H" )
    call auxputch("n2","i",nt-2,"a1.H" )
    call auxputch("n2","i",nt-2,"a2.H" )
    call auxputch("o1","r",0.,"a1.H" )
    call auxputch("o1","r",0.,"a2.H" )
    call auxputch("o2","r",0.,"a1.H" )
    call auxputch("o2","r",0.,"a2.H" )
    call auxputch("d1","r",1.,"a1.H" )
    call auxputch("d1","r",1.,"a2.H" )
    call auxputch("d2","r",1.,"a1.H" )
    call auxputch("d2","r",1.,"a2.H" )
    call auxputch("n1","i",nr,"a4.H" )
    call auxputch("n1","i",nr,"a3.H" )
    call auxputch("n2","i",nt-2,"a4.H" )
    call auxputch("n2","i",nt-2,"a3.H" )
    call auxputch("o1","r",0.,"a4.H" )
    call auxputch("o1","r",0.,"a3.H" )
    call auxputch("o2","r",0.,"a4.H" )
    call auxputch("o2","r",0.,"a3.H" )
    call auxputch("d1","r",1.,"a4.H" )
    call auxputch("d1","r",1.,"a3.H" )
    call auxputch("d2","r",1.,"a4.H" )
    call auxputch("d2","r",1.,"a3.H" )
  end subroutine a_print
 
  !----------------------------------------------------------------
end module  RWE2D_polar_wemig
