NTCC module -- kdsaw -- Kadomtsev Sawtooth Model with Porcelli Upgrade.
=======================================================================

D. McCune (dmc) 29 Jan 2002
Update: dmc 15 Nov 2005 -- Controls added for Porcelli "partial reconnection"
mixing model; documentation updated.

Abstract:

This module implements the Kadomtsev Sawtooth Mixing Model[1] and the 
Porcelli Sawtooth Mixing Model[2] for axisymmetric tokamaks in general 
geometry.  The module is initialized by calls to "kdsaw_pre_init" and 
to "kdsaw_init", the former specifying the mixing model selection, and
the latter providing a grid, the iota(bar) (1/q) profile, and zone volumes;
kdsaw_init computes the Kadomtsev or Porcelli mapping and stores it in the 
module's internal memory.  Subsequent calls are used to extract specific 
information about the sawtooth mixing.  The grid passed to kdsaw_init is 
specified in normalized square-root-toroidal-flux, a coordinate which 
functions in a way strictly analogous to the minor radial coordinate "r"
in Kadomtsev's original paper, but is suitable for general geometry.  The
grid is given explicitly and need not be evenly spaced.  The current profile
is expressed as (1/q), and the zone volumes (in any suitable units) are
used to adjust the mapping such that a change in a particle or energy 
density profile according to the ideal Kadomtsev mapping volume integrates
to precisely zero:  typically, a peaked density profile is flattened, with
the density reduced in the center and raised at the edge of the mixing
region.  A user-adjustable current profile remapping is also provided.

Both REAL and REAL*8 interfaces are provided to all routines in the module.

[1] Sov. J. Plasma Phys. Vol. 1 No. 5 Sept. Oct. 1975, p. 389.
[2] Porcelli et al., Plasma Physics & Controlled Fusion 38 (1996) p. 2163,
    (See section 4).

------------------------------------------------
Note: the Porcelli sawtooth trigger model is NOT included in this software.
A separate NTCC module is expected to be supplied for the sawtooth trigger
calculation.  KDSAW calculates mixing only.

------------------------------------------------
Kadomtsev mixing vs. Porcelli mixing:  the details are provided in the
references, but summary of the mixing characteristics are given here.
The coordinate x refers to normalized sqrt(toroidal flux); x=0 means the
magnetic axis, x=1 would refer to the plasma boundary.  For both models,
there is no "instantaneous" (i.e. MHD timescale) change in the plasma
current profile or plasma parameters except inside the sawtooth mixing
region.  Let x_q1 be the x location of the q=1 surface which for purposes
of this discussion is assumed to be unique; let x_mix > x_q1 be the outer
limit of the sawtooth mixing region.

----------

Kadomtsev current mixing:  within the mixing region, the q profile value is
increased in the region 0 < x < x_q1, and reduced in the region 
x_q1 < x < x_mix, and somewhat flattened.  The post sawtooth q profile 
has q(0)=1 and q'(x)>0, q(x)>1, for x>0, out to x_mix, where a 
"current sheet" (discontinuity in the q profile) forms, since the pre- 
and post-sawtooth q profile values are assumed unchanged outside the 
mixing region.  Some experimental evidence suggests q<1 persisting 
after sawtooth crashes; this was part of the motivation for the later 
development of the Porcelli model.

----------

Kadomtsev plasma mixing:  in the standard Kadomtsev model, the plasma
at the pre-sawtooth q=1 surface is transported to the axis x=0.  The 
pre-sawtooth q=1 surface represents the maximum value of a pre-sawtooth
helical flux function

  Psi_1(x) = [poloidal magnetic flux from 0 to x] -
             [poloidal magnetic flux from 0 to x if q=1 everywhere].

Away from the pre-sawtooth q=1 surface, the plasma is zoned by matching
value ranges of Psi_1, giving two plasma regions one inside and one outside
the q=1 surface for each pair of Psi_1 values.  The plasma from these two
regions is mixed and transported to a region subtending equal toroidal flux,
with zone pairs corresponding to higher pre-sawtooth Psi_1 mapped nearer to
the axis, and pairs with lower Psi_1 values mapped towards the edge of the 
mixing region.

The net effect of this is generally to flatten the plasma temperature
and density profiles but there can be residual structure depending on
the details of the pre-sawtooth profiles.

KDSAW allows the option of "ergodic" mixing instead, in which all plasma
is redistributed evenly throughout the mixing region, resulting in totally
flat post-sawtooth particle and energy density profiles.

----------

The Porcelli model depends on an additional parameter Fp, which is the
"fractional q=1 island width" prior to the sawtooth crash.  An inner
boundary of the q=1 island is defined:

  x_min = (1-Fp)*x_q1

where Fp=1 implies x_min=0, which would result in mixing similar to the
Kadomtsev model.

In the Porcelli model, the pre-sawtooth Psi_1(x_min) is found; the outer
boundary of the q=1 island and the mixing region, x_mix, is found by the
condition Psi_1(x_mix)=Psi_1(x_min).  In the Kadomtsev model, by 
contrast, Psi_1(x_mix)=Psi_1(0)=0.  The Porcelli mixing region is smaller:
the x_mix of the Porcelli model is less than the x_mix of Kadomtsev by an 
amount depending on the value of Fp.

The Porcelli model considers the region x_min <= x <= x_mix to be the q=1
island which gets "reconnected" during the sawtooth crash.  The inside
axial region 0 <= x < x_min is not included in this reconnection, but is
considered to go through a period of ergodic field lines and rapid plasma
mixing nonetheless.

Porcelli current mixing:

  for x_min <= x <= x_mix,   q = 1   is imposed.
  for 0 <= x < x_min,        q = q[pre-sawtooth](x_min)   is imposed.

So the q profile becomes a two-step function within the mixing region,
with discontinuities or "current sheets" at x_min and x_mix.  (Subsequent
application of a poloidal field diffusion equation will result in rapid
decay or spreading of the current sheets).

Porcelli plasma mixing:

The plasma is redistributed evenly but separately within the two 
subregions, resulting in two-step particle and energy density profiles
after the sawtooth crash.

KDSAW allows the option of "ergodic" mixing over the ENTIRE mixing region
instead, in which all plasma is redistributed evenly throughout, resulting
in flat post-sawtooth temperature and density profiles (with no step at
x=x_min).  This choice (ion_mix=4) does not affect current mixing.

One of the four KDSAW mixing options are chosen by means of a call to
KDSAW_PRE_INIT:

   ion_mix    zf_reconnect

      1        (ignored)        Standard Kadomtsev mixing (default).
      2        (ignored)        Kadomtsev current mixing; temperature
                                and density profiles flat post-sawtooth.

      3           =Fp           Standard Porcelli mixing.
      4           =Fp           Porcelli current mixing; temperature
                                and density profiles flat post-sawtooth.

------------------------------------------------
Current sheet smoothing:

The Kadomtsev mixing model produces 1 current sheet at x=x_mix; the
Porcelli model produces 2 current sheets at x=x_min and x=x_mix.

For finite difference numerical codes, the effective thickness of the
sheet will be determined by the grid size.  This means that the dynamical
details of the response of poloidal field diffusion post-sawtooth will
be affected somewhat by grid size.

To reduce this numerical effect, it is possible to specify a finite
width (in sqrt toroidal flux space) over which the current sheet(s) is
(are) to be smoothed out.  See the description of subroutine KDSAW_DXI_QJUMP
below.  The default (if KDSAW_DXI_QJUMP is never called) is for no such 
smoothing.

================================================
================================================

Test program:

The program "kdsaw_test" demonstrates the use of kdsaw subroutines.
It requires no input file.  To test the module, build the program and

  unix> kdsaw_test > test.output
  unix> diff test.output kdsaw_test.output


Differences, reflecting slight variations of floating point arithmetic
from one computing platform to another, should be minimal.

The program "kdsaw_test" exercises (and demonstrates) all mixing options;
current sheets are not smoothed.

================================================
================================================

Subroutine call interface:

A fortran-77 style interface is provided.  The REAL*8 interface is 
described.  Append "_r4" to get the REAL interface (e.g. use 
"kdsaw_init_r4" instead of "kdsaw_init" if REAL arguments are to be
provided.

Notes: 

  (1) If anything other than pure Kadomtsev mixing is desired, execute
      "call kdsaw_pre_init" first, to set the Porcelli mixing controls.
  (2) Next, execute "call kdsaw_init"!  All other calls refer to data 
      that was input in a preceding kdsaw_init call.  Each "kdsaw_init" 
      call updates the Kadomtsev/Porcelli sawtooth map.

-----------------------------
**** (optional) set Porcelli mixing controls ****

Default values (used if KDSAW_PRE_INIT is never called) correspond to:

   zf_reconnect = 1.0  ! real*8
   ion_mix = 1

once the defaults are modified, the new values stay in effect until the 
next kdsaw_pre_init call.  However, it is recommended to *always* call
kdsaw_pre_init prior to kdsaw_init, in codes that support non-Kadomtsev
sawtooth mixing as an option.

SUBROUTINE KDSAW_PRE_INIT(zf_reconnect,ion_mix)

  real*8, intent(in) :: zf_reconnect  ! reconnection fraction (Fp)

  integer, intent(in) :: ion_mix      ! plasma species mixing option

The meanings are as described above in the section describing the Kadomtsev
and Porcelli mixing options.

KDSAW_PRE_INIT reports no errors or messages.  The following value
ranges are silently enforced with min & max function calls:

    1  <= ion_mix      <= 4   ..and..
   0.0 <= zf_reconnect <= 1.0

NOTE: if the Porcelli model is used, current sheets, i.e. discontinuous
jumps in the q profile, will be produced.  In practice the current sheet
width will be spread out over a grid spacing in real codes.  This might
cause problems when comparing sensitivity of details of code dynamics to
grid size; therefore, a non-zero current sheet width may be specified via
a separate call to "kdsaw_dxi_qjump"...

The routine "kdsaw_pre_init" must be called *prior* to "kdsaw_init" to have
an effect.

-----------------------------
**** (optional) set current sheet width ****

Default value: dxi_qjump = 0.0  ! real*8

SUBROUTINE KDSAW_DXI_QJUMP(dxi_qjump)
  real*8, intent(in) :: dxi_qjump  ! current sheet width in xi space
                                   ! (i.e. sqrt toroidal flux normalized).

Apply a smoothing operation of half-width   dxi_qjump   at each current
sheet (i.e. the discontinuous changes in the post-sawtooth q profile at
locations corresponding to the inner and outer boundaries of the q=1
magnetic island.  For Kadomtsev mixing, only the outer boundary exists.

The maximum allowed value if dxi_qjump is 0.1 -- this is silently enforced.

Once the default is modified, the new value stays in effect until the 
next kdsaw_dxi_qjump call.  However, it is recommended to *always* call
kdsaw_dxi_qjump prior to kdsaw_init, in codes that support current sheet
smoothing as an option.

This routine must be called *prior* to "kdsaw_init" to have an effect.

-----------------------------
-----------------------------
-----------------------------

**** compute Kadomtsev/Porcelli sawtooth map ****  

**** subroutine kdsaw_init ****

SUBROUTINE KDSAW_INIT(ilun_msg,izones,zxi,ziotb,zdvol,iimix,ierr)

  !  input data:

  integer, intent(in) :: ilun_msg
                         ! fortran i/o channel for error messages

  integer, intent(in) :: izones
                         ! number of plasma zones
  real*8, dimension(izones+1) :: zxi
                         ! zone *bdy* locations in normalized sqrt(tor.flux)
                         ! in strict ascending order with...
                         ! zxi(1)=0, zxi(izones+1)=1 expected
                         ! zxi(j) is inner bdy, zxi(j+1) outer bdy of j'th zone

  real*8, dimension(izones+1) :: ziotb
                         ! iota(bar) (1/q) profile at zone bdys
                         ! These are pre-sawtooth values.

  real*8, dimension(izones) :: zdvol
                         ! "volume" of the zones 
                         ! in arbitrary units, only the relative values
                         ! amongst zones are important

  ! output data:

  integer, intent(out) :: iimix
                         ! iimix gives index to last zone involved in
                         ! Kadomtsev mixing, or 0 if none (not an error).

  integer, intent(out) :: ierr
                         ! error flag-- error in input data,
                         !   such as, izones.lt.5, or, zxi out of order.

Note:  zxi need not be evenly spaced, but it is crucial that its values
correspond to sqrt(Phi/Philim) where Phi = enclosed toroidal flux.

-----------------------------
**** retrieve mixing region information ****

**** subroutine kdsaw_xmix ****

SUBROUTINE KDSAW_XMIX(zxmin,zxq1,zxmix)

  real*8, intent(out) :: zxmin   ! x_min -- inner q=1 island boundary
  real*8, intent(out) :: zxq1    ! x_q1 -- pre-sawtooth q=1 surface location
  real*8, intent(out) :: zxmix   ! x_mix -- outer q=1 island boundary; outer
                                 !   limit of sawtooth mixing.

If Porcelli mixing is used with Fp=0 (not recommended), zxmin=zxq1=zxmix;
otherwise zxmin < zxq1 < zxmix is to be expected.  For Kadomtsev options,
zxmin=0.0 will be returned.  The higher the value of Fp, the larger zxmix
and the small zxmin will be (i.e. the larger the sawtooth mixing region).

-----------------------------
**** compute post-sawtooth iota(bar) profile with ad hoc adjustment ****

**** subroutine kdsaw_iotb ****           

SUBROUTINE KDSAW_IOTB(ziotb,zmix,ierr)

  ! return the post-sawtooth iota(bar) profile
  ! call this only after kdsaw_init has been called.

  ! user control:
  ! mixing parameter: zmix=1 ==> return fully mixed post-sawtooth profile
  !                   zmix=0 ==> return pre-sawtooth profile unchanged
  !                   0<zmix<1 ==> return
  !                       (1-zmix)*[pre-sawtooth profile] +
  !                         (zmix)*[post-sawtooth profile]
  !
  ! CAUTION: zmix is NOT the control for Porcelli partial reconnection; it
  ! is a means for further reduction of the current mixing if desired.
  !
  ! warning issued if zmix out of range; zmixl = max(0,min(1,zmix)) is used.
  !
  ! note: dimension of ziotb array is equal to "izones+1" where "izones"
  ! was the dimensioning argument provided in an earlier call to kdsaw_init.

  real*8, intent(out) :: ziotb(izones+1)
                                   ! (mixed) post-sawtooth iota(bar) profile

  real*8, intent(in) :: zmix       ! mixing parameter in range [0,1]

  integer, intent(out) :: ierr     ! completion code, ierr=0 means OK

-----------------------------
**** compute change in a density profile ****

**** subroutine kdsaw_deln ****

      SUBROUTINE KDSAW_DELN(zpre,zdel,ierr)
C
C  USING THE MIXING FACTORS CALCULATED IN THE SAWMIX ROUTINE,
C  EVALUATE THE EXPECTED CHANGE ZDEL IN PRE-SAWTOOTH DENSITY
C  PROFILE ZPRE DUE TO THE SAWTOOTH (ideal Kadomtsev/Porcelli mixing).
C
C  THE VOLUME INTEGRAL ACROSS THE WHOLE PLASMA OF ZDEL WILL COME
C  OUT TO ZERO
C
C  Note: the dimensioning of the arrays is determined by the argument
C  "izones" in a prior call to kdsaw_init.
C
      real*8, intent(in) :: zpre(izones) ! input density profile
      real*8, intent(out) :: zdel(izones) ! ideal Kadomtsev change profile
C
      integer, intent(out) :: ierr      ! completion code, 0=OK

-----------------------------
**** given a random number, choose post-sawtooth zone of particle
     starting in a particular pre-sawtooth zone.

**** subroutine kdsaw_find ****

SUBROUTINE KDSAW_FIND(zselect,izone_in,izone_out,zfrac)

  ! given an input zone and a (random) number between 0 and 1, select
  ! an output zone, according to the Kadomtsev redistribution mapping

  real*8,intent(in) :: zselect  ! a number btw 0 and 1
  integer,intent(in) :: izone_in  ! pre-sawtooth originating zone
  integer, intent(out) :: izone_out  ! post-sawtooth destination zone selected
  real*8,intent(out) :: zfrac        ! location within post-sawtooth zone
     ! zfrac=0 corresponds to inner boundary; 1 corresponds to outer
     ! boundary of mapped-to zone.

Note: if the input zone index indicates the last zone which has non-null
intersection with the mixing region (i.e. a zone which is partly inside
and partly outside the mixing region), the source point being mapped from
is assumed to be inside the mixing region.  For any call, if izone_out is
the last zone in the mixing region, the maximum value of zfrac will reflect
the fraction of the last zone that is in the mixing region.

-----------------------------
**** given a post-sawtooth zone index show what fraction of material
     from which pre-sawtooth zones were mixed to form the post-sawtooth
     zone's contents. ****

**** subroutine kdsaw_from ****

SUBROUTINE KDSAW_FROM(izone,inzmax,inzact,ifrom,zfrom,ierr)

  ! give statistics on what fraction of material from which pre-sawtooth
  ! zones constitute the material in a selected post-sawtooth zone.

  integer, intent(in) :: izone          ! selected post-sawtooth zone 
  integer, intent(in) :: inzmax         ! max number of pre-sawtooth zones

  integer, intent(out) :: inzact        ! actual number of pre-sawtooth zones
  integer, intent(out) :: ifrom(inzmax) ! pre-sawtooth zone indices
  real*8, intent(out) :: zfrom(inzmax)  ! pre-sawtooth material fractions
  integer, intent(out) :: ierr          ! error code (0=OK)

-----------------------------
**** given a pre-sawtooth zone index, show to which post-sawtooth zones
     the pre-sawtooth zone's material is distributed.

**** subroutine kdsaw_to ****

SUBROUTINE KDSAW_TO(izone,inzmax,inzact,ito,zto,ierr)

  ! give distribution of material of pre-sawtooth zone "izone" to
  ! affected post-sawtooth zones.

  integer, intent(in) :: izone          ! selected pre-sawtooth zone 
  integer, intent(in) :: inzmax         ! max number of post-sawtooth zones

  integer, intent(out) :: inzact        ! actual number of post-sawtooth zones
  integer, intent(out) :: ito(inzmax)   ! post-sawtooth zone indices
  real*8, intent(out) :: zto(inzmax)    ! post-sawtooth material fractions
  integer, intent(out) :: ierr          ! error code (0=OK)
