program kdsaw_test
  !
  ! test "kdsaw" -- sawtooth module
  !
  ! construct artificial iota(bar) and density profiles; examine
  ! ideal Kadomtsev sawtooth based on these profiles
  !   dmc 25 Jan 2002

  implicit NONE

  integer, parameter :: nzones = 20
  integer, parameter :: nzp1 = nzones + 1

  real xi(nzp1)                     ! grid
  real xiotb(nzp1)                  ! pre-sawtooth iota(bar) profile
  real dvol(nzones)                 ! nominal zone volumes
  real rho(nzones)                  ! pre-sawtooth density

  real xiotb_new(nzp1),xiotb_mix(nzp1)
                                    ! post-sawtooth (ideal Kadomtsev) iota(bar)
                                    ! and 25% of Kadomtsev changed iota(bar)

  real :: zmix = 0.25               ! mixing parameter for xiotb_mix

  real :: zfp(4)                    ! Porcelli fractions

  real drho(nzones)                 ! ideal Kadomtsev delta(rho)

  integer :: kzonex(nzp1),kzoney(nzp1)  ! source for post-sawtooth zone

  real :: zdxi = 0.0                ! qjump smoothing width (0 for none)

  real zzonex(nzp1),zfrac

  real :: q0 = 0.7
  real :: qa = 3.5

  real :: rho0 = 3.0e13
  real :: rhoa = 1.0e12

  real :: one = 1.0

  real zq,zrho,x,xc
  integer ipor,i,j,jp1,imix,ierr,igot

  !-----------------------------------------------
  !  initialize local data

  zfp(1) = 1.0          ! Kadomtsev (standard)
  zfp(2) = 1.0          ! Kadomtsev with flat n & T profiles post-sawtooth
  zfp(3) = 0.7          ! Porcelli, Fp=0.7 (standard)
  zfp(4) = 0.7          ! Porcelli, Fp=0.7, flat n & T profiles post-sawtooth

  do ipor = 1,4

     xi=0;  xiotb=0
     dvol=0;  rho=0;  drho=0
     xiotb_new=0;  xiotb_mix=0

     !  setup profiles

     xiotb(1)=1/q0

     do j=1,nzones

        jp1=j+1
        x=float(j)/nzones  ! outer zone bdy  (.05, .1, .15, ...)
        xi(jp1)=x          ! = surface j+1
        zq=q0+x**2*(qa-q0)
        xiotb(jp1)=one/zq  ! at outer bdy = surface j+1

        xc=(float(j-1)+0.5)/nzones  ! zone center (.025, .075, .125, ...)
        dvol(j)=xc         ! "dvol" need not have units; here ~ x*dx
        zrho=rho0+xc**2*(rhoa-rho0)
        rho(j)=zrho
 
     enddo

     !--------------------------------------------------
     !  fire up Kadomtsev model -- this computes the Kadomtsev mapping

     imix=ipor
     if(zfp(imix).gt.0.99) then
        write(6,*) ' KDSAW -- Kadomtsev calculation...'
        if(imix.eq.1) write(6,*) ' STANDARD Kadomtsev mixing'
        if(imix.eq.2) write(6,*) ' Mixing => flat n & T profiles'
     else
        write(6,*) ' ==============================================='
        write(6,'(a,1pe10.3)') ' KDSAW -- Porcelli mixing fraction: ',zfp(imix)
        if(imix.eq.3) write(6,*) ' STANDARD Porcelli mixing'
        if(imix.eq.4) write(6,*) ' Mixing => flat n & T profiles'
     endif

     call kdsaw_pre_init_r4(zfp(imix),imix)  ! default is Kadomtsev mixing.

     call kdsaw_dxi_qjump_r4(zdxi)   ! (no) smoothing -- default is none.

     call kdsaw_init_r4(6,nzones,xi,xiotb,dvol,imix,ierr)  ! calculate mixing

     if(ierr.ne.0) then
        write(6,*) ' ?kdsaw_test: after kdsaw_init_r4, ierr=',ierr
        stop
     endif

     write(6,1001) imix
1001 format(/ &
          ' kdsaw_test: last sawtooth mixing zone = ',i2)

     !--------------------------------------------------
     !  full current profile mixing:

     call kdsaw_iotb_r4(xiotb_new,one,ierr)

     if(ierr.ne.0) then
        write(6,*) ' ?kdsaw_test: after kdsaw_iotb_r4 (1), ierr=',ierr
        stop
     endif

     !--------------------------------------------------
     !  partial current mixing:

     call kdsaw_iotb_r4(xiotb_mix,zmix,ierr)

     if(ierr.ne.0) then
        write(6,*) ' ?kdsaw_test: after kdsaw_iotb_r4 (mix), ierr=',ierr
        stop
     endif

     !--------------------------------------------------
     !  ideal density profile change

     call kdsaw_deln_r4(rho,drho,ierr)

     if(ierr.ne.0) then
        write(6,*) ' ?kdsaw_test: after kdsaw_deln_r4, ierr=',ierr
        stop
     endif

     !--------------------------------------------------
     !  printout...

     write(6,1002)
1002 format(/'  j  x            iota(orig)   iota(mix)    iota(25%)')
     do j=1,nzones+1
        write(6,1003) j,xi(j),xiotb(j),xiotb_new(j),xiotb_mix(j)
     enddo
1003 format(1x,i2,4(1x,1pe12.5))
1004 format(1x,i2,5(1x,1pe12.5))

     write(6,1005)
1005 format(/'  j  x            dvol         rho(pre)    +drho        =rho')
     do j=1,nzones
        write(6,1004) j,(xi(j)+xi(j+1))/2,dvol(j),rho(j),drho(j),rho(j)+drho(j)
     enddo

     !--------------------------------------------------
     ! "find" routines used by Monte Carlo code:  given a random
     ! number choose a post-sawtooth zone for a particle from a given
     ! pre-sawtooth zone; zfrac gives displacement within zone.  Note
     ! that the last numerical zone is partly inside and partly
     ! partly outside the mixing region; the mapping is assumed to
     ! apply to particles from inside the mixing region only.

     do j=1,nzones+1
        call kdsaw_find_r4(xi(j),imix-1,kzonex(j),zfrac)
        call kdsaw_find_r4(xi(j),imix,kzoney(j),zfrac)
     enddo

     !--------------------------------------------------
     !  printout...

     write(6,1009)
1009 format(/ &
          ' the following table shows the post-sawtooth zones to which'/ &
          ' material from the last two pre-sawtooth zones in the mixing'/ &
          ' region is sent, as a function of a sampling parameter "x"'// &
          '  j  x          (from imix-1) (from imix)')
     do j=1,nzones+1
        write(6,1010) j,xi(j),kzonex(j),kzoney(j)
     enddo
1010 format(1x,i2,1x,1pe12.5,1x,i2,11x,i2)

     !--------------------------------------------------
     !  "kdsaw_from" can be used to find out where the "material" in
     !  a post-sawtooth zone comes from...
     ! 
     !  outside the mixing region the material comes from the same zone.
     !
     !  calls and printout...

     write(6,*) ' '
     do j=1,imix+1
        call kdsaw_from_r4(j,nzones,igot,kzonex,zzonex,ierr)
        if(ierr.ne.0) write(6,*) ' ?? kdsaw_from_r4 returned ierr = ',ierr
        write(6,1012) j
1012    format(' material in post-sawtooth zone ',i2,' contains:')
        do i=1,igot
           write(6,1013) zzonex(i),kzonex(i)
1013       format('   ',f10.3,' x 100% of material from pre-sawtooth zone ',i2)
        enddo
     enddo

     !--------------------------------------------------
     !  "kdsaw_to" can be used to find out where the "material" in
     !  a pre-sawtooth zone goes to...
     !
     !  calls and printout...

     write(6,*) ' '
     do j=1,imix+1
        call kdsaw_to_r4(j,nzones,igot,kzonex,zzonex,ierr)
        if(ierr.ne.0) write(6,*) ' ?? kdsaw_from_r4 returned ierr = ',ierr
        write(6,1022) j
1022    format(' material in pre-sawtooth zone ',i2,' is distributed:')
        do i=1,igot
           write(6,1023) zzonex(i),kzonex(i)
1023       format('   ',f10.3,' x 100% goes to post-sawtooth zone ',i2)
        enddo
     enddo

  enddo

  write(6,*) ' '
  write(6,*) ' kdsaw_test:  normal exit.'

  stop

end program kdsaw_test
