rungridgen Program

program rungridgen
  use file_utils
  use gridgen4mod, only: gridgen4_1
  use constants, only: pi, twopi, run_name_size
  use mp, only: init_mp, finish_mp
  implicit none

  ! input parameters
  character(len = run_name_size) :: source, gingrid, gsource
  integer :: nthetaout, nlambdaout, nperiodout
  integer :: npadd, iperiod
  real :: alknob, epsknob, extrknob, bpknob, smoothknob
  integer :: nfinegrid
  real :: thetamax, deltaw, widthw, tension
  logical :: screenout
  logical :: auto_width, three_dim
  real :: cv_fraction, delth_max
  integer :: max_autoiter
  logical :: list = .false.
  logical :: no_end_point
  character(len = run_name_size) :: ncsource

  ! work variables
  integer :: ntheta, nlambda, ntgrid
  integer :: ntgridin, nperiodin, nthetain
  real, dimension (:), allocatable :: thetain, bmagin, bmagsm
  real, dimension (:), allocatable :: thetaout, bmagout, alambdaout

  real :: drhodpsi, rmaj, shat, kxfac, qval
  real, dimension (:), allocatable :: thetagrid
  real, dimension (:), allocatable :: gbdrift, gradpar, grho
  real, dimension (:), allocatable :: cvdrift, gds2, bmag
  real, dimension (:), allocatable :: gds21, gds22
  real, dimension (:), allocatable :: cvdrift0, gbdrift0
  real, dimension (:), allocatable :: Rplot, Rprime
  real, dimension (:), allocatable :: Zplot, Zprime
  real, dimension (:), allocatable :: aplot, aprime
  
  call init_mp
  call init_file_utils (list, name="grid")
  call read_parameters
  if (three_dim) then
     call allocate_arrays_3d
     call get_initial_grids_3d
  else
     call allocate_arrays
     call get_initial_grids
  end if
  call do_smoothing
  call generate_grids
  call write_output
  call finish_file_utils
  call finish_mp
  stop

contains
  !> FIXME : Add documentation  
  subroutine read_parameters
    implicit none
    namelist /testgridgen/ source, gsource, &
         nthetaout,nlambdaout,nperiodout, &
         npadd,alknob,epsknob,bpknob,extrknob,smoothknob, nfinegrid, &
         thetamax, deltaw, widthw, tension, gingrid, screenout, &
         auto_width, cv_fraction, delth_max, max_autoiter, three_dim, &
         iperiod, &
         no_end_point, ncsource
    gsource = "eik6.out"
    source = "eik.out"
    nthetaout = 32
    nlambdaout = 20
    nperiodout = 2
    npadd = 2
    alknob = 0.1
    bpknob = 1.e-8
    epsknob = 1e-9
    extrknob = 10.0
    smoothknob = 0.0
    nfinegrid=200
    thetamax = 0.0
    deltaw = 0.0
    widthw = 1.0
    tension = 1.0
    screenout = .false.
    auto_width = .false.
    cv_fraction = 0.6
    delth_max = 0.5
    max_autoiter = 3
    gingrid = "gingrid"
    three_dim = .false.
    iperiod = 1
    no_end_point = .false.
    ncsource = ""
    read (unit=input_unit("testgridgen"), nml=testgridgen)
    nthetaout = nthetaout + 1
  end subroutine read_parameters

  !> FIXME : Add documentation  
  subroutine allocate_arrays
    use gs2_io_grids, only:nc_grid_file_open
    use gs2_io_grids, only: nc_get_grid_sizes, nc_get_grid_scalars
    implicit none
    integer :: unit
    character(200) :: line
    logical :: ncsource_exist

    inquire(file=ncsource, exist=ncsource_exist)
    if (ncsource_exist) then
       call nc_grid_file_open(ncsource, "r")
       call nc_get_grid_sizes(nthetain, ntgridin, nperiodin)
       call nc_get_grid_scalars(shat, drhodpsi, kxfac, qval, rmaj)
    else
       call get_unused_unit(unit)
       open (unit=unit, file=trim(source), status="old")
       read (unit=unit, fmt="(a)") line
       read (unit=unit, fmt=*,err=10) ntgridin, nperiodin, nthetain, drhodpsi, rmaj, shat, kxfac, qval
10     continue
       close (unit=unit)
    end if

    allocate (thetain(nthetain+1),bmagin(nthetain+1),bmagsm(nthetain+1))
    allocate (thetaout(nthetaout),bmagout(nthetaout),alambdaout(nlambdaout))
    allocate (thetagrid(-ntgridin:ntgridin))
    allocate (gbdrift(-ntgridin:ntgridin))
    allocate (gradpar(-ntgridin:ntgridin))
    allocate (grho(-ntgridin:ntgridin))
    allocate (cvdrift(-ntgridin:ntgridin))
    allocate (gds2(-ntgridin:ntgridin))
    allocate (bmag(-ntgridin:ntgridin))
    allocate (gds21(-ntgridin:ntgridin))
    allocate (gds22(-ntgridin:ntgridin))
    allocate (cvdrift0(-ntgridin:ntgridin))
    allocate (gbdrift0(-ntgridin:ntgridin))
    allocate (Rplot(-ntgridin:ntgridin))
    allocate (Rprime(-ntgridin:ntgridin))
    allocate (Zplot(-ntgridin:ntgridin))
    allocate (Zprime(-ntgridin:ntgridin))
    allocate (aplot(-ntgridin:ntgridin))
    allocate (aprime(-ntgridin:ntgridin))
  end subroutine allocate_arrays

  !> FIXME : Add documentation    
  subroutine allocate_arrays_3d
    implicit none
    integer :: unit, ntor
!    character(200) :: line

    call get_unused_unit(unit)
    open (unit=unit, file=trim(source), status="old")
    read (unit=unit, fmt=*) ntgridin, nthetain, ntor
    close (unit=unit)
    
    drhodpsi = 1.0  ! assumes our radial variable matches VVBAL

    nthetain = nthetain * iperiod

    allocate (thetain(nthetain+1),bmagin(nthetain+1),bmagsm(nthetain+1))
    allocate (thetaout(nthetaout),bmagout(nthetaout),alambdaout(nlambdaout))
    allocate (thetagrid(-ntgridin:ntgridin))
    allocate (gbdrift(-ntgridin:ntgridin))
    allocate (gradpar(-ntgridin:ntgridin))
    allocate (grho(-ntgridin:ntgridin))
    allocate (cvdrift(-ntgridin:ntgridin))
    allocate (gds2(-ntgridin:ntgridin))
    allocate (bmag(-ntgridin:ntgridin))
    allocate (gds21(-ntgridin:ntgridin))
    allocate (gds22(-ntgridin:ntgridin))
    allocate (cvdrift0(-ntgridin:ntgridin))
    allocate (gbdrift0(-ntgridin:ntgridin))
    allocate (Rplot(-ntgridin:ntgridin))
    allocate (Rprime(-ntgridin:ntgridin))
    allocate (Zplot(-ntgridin:ntgridin))
    allocate (Zprime(-ntgridin:ntgridin))
    allocate (aplot(-ntgridin:ntgridin))
    allocate (aprime(-ntgridin:ntgridin))

    cvdrift0 = 0.   ! assumes theta0 = 0.
    gbdrift0 = 0.   ! assumes theta0 = 0.
    gds21 = 0.      ! assumes theta0 = 0.
    gds22 = 0.      ! assumes theta0 = 0.
    grho  = 1.      ! assumes linear calculation
    Rplot = 1.      ! pure fiction
    Rprime = 1.     ! pure fiction
    Zplot = 1.      ! pure fiction
    Zprime = 1.     ! pure fiction
    aplot = 1.      ! pure fiction
    aprime = 1.     ! pure fiction
    
  end subroutine allocate_arrays_3d

  !> FIXME : Add documentation  
  subroutine get_initial_grids
    use file_utils, only: get_unused_unit
    use gs2_io_grids, only: nc_get_grids
    use gs2_io_grids, only: nc_grid_file_close
    implicit none
    integer :: unit
    integer :: i
!    real :: discard
    character(200) :: line
    logical :: ncsource_exist

    inquire(file=ncsource, exist=ncsource_exist)
    if (ncsource_exist) then
       call nc_get_grids(ntgridin, &
            bmag, gradpar, gbdrift, gbdrift0, cvdrift, cvdrift0, &
            gds2, gds21, gds22, grho, thetagrid, &
            Rplot=Rplot, Rprime=Rprime, Zplot=Zplot, Zprime=Zprime, aplot=aplot, aprime=aprime)
       call nc_grid_file_close()
    else
       call get_unused_unit(unit)
       open (unit=unit, file=trim(source), status="old")
       read (unit=unit, fmt="(a)") line
       read (unit=unit, fmt="(a)") line

       ! gbdrift gradpar grho theta
       read (unit=unit, fmt="(a)") line
       do i = -ntgridin,ntgridin
          read (unit=unit, fmt=*) gbdrift(i),gradpar(i),grho(i),thetagrid(i)
       end do

       ! cvdrift gds2 bmag theta
       read (unit=unit, fmt="(a)") line
       do i = -ntgridin,ntgridin
          read (unit=unit, fmt=*) cvdrift(i),gds2(i),bmag(i),thetagrid(i)
       end do

       ! gds21 gds22 theta
       read (unit=unit, fmt="(a)") line
       do i = -ntgridin,ntgridin
          read (unit=unit, fmt=*) gds21(i),gds22(i),thetagrid(i)
       end do

       ! cvdrift0 gbdrift0 theta
       read (unit=unit, fmt="(a)") line
       do i = -ntgridin,ntgridin
          read (unit=unit, fmt=*) cvdrift0(i),gbdrift0(i),thetagrid(i)
       end do

       close (unit=unit)

       call get_plotdata
    end if

    thetain = thetagrid(-nthetain/2:nthetain/2)
    bmagin = bmag(-nthetain/2:nthetain/2)

  end subroutine get_initial_grids

  subroutine get_plotdata
    implicit none
    integer :: unit
    integer :: i
    character(200) :: line
    real, allocatable :: dummy(:)

    allocate (dummy(-ntgridin:ntgridin))

    Rplot = 0.; Rprime = 0.
    Zplot = 0.; Zprime = 0.
    aplot = 0.; aprime = 0.

    call get_unused_unit(unit)
    open (unit=unit, file=trim(gsource), status="old", err=100)

    read (unit=unit, fmt="(a)", err=100) line
    read (unit=unit, fmt="(a)", err=100) line
    read (unit=unit, fmt="(a)", err=100) line

    do i = -ntgridin, ntgridin
       read (unit=unit, fmt=*, err=100) dummy(i), Rplot(i), Zplot(i), aplot(i), Rprime(i), Zprime(i), aprime(i), dummy(i)
    end do

    close(unit=unit)

    return

100 continue
    close(unit=unit)
    write(6,*) 'read error: ', gsource

    return

  end subroutine get_plotdata

  !> FIXME : Add documentation  
  subroutine get_initial_grids_3d
    use file_utils, only: get_unused_unit
    implicit none
    integer :: unit
    integer :: i, ntor
    real :: dpdpsi, pres, bavg, cvavg

    call get_unused_unit(unit)
    open (unit=unit, file=trim(source), status="old", action="read")

    read (unit=unit, fmt=*) ntgridin, nthetain, ntor

    nthetain = nthetain * iperiod
    do i = -ntgridin,ntgridin
       read (unit=unit, fmt=*) thetagrid(i),bmag(i),gradpar(i),gds2(i),cvdrift(i),gbdrift(i),dpdpsi,pres
    end do

    close (unit=unit)

    bavg = sum(bmag(-nthetain/2:nthetain/2))/(nthetain+1)
    cvavg = sum(cvdrift)/size(cvdrift)
    
    write(*,*) 'bavg = ',bavg
    write(*,*) 'cvavg = ',cvavg,' and remember that cvdrift assumed beta_prime=0'

!    bmag = bmag / bavg

    gds2 = gds2 * ntor**2 
    gbdrift = gbdrift * ntor * 2. / bmag**2
!    cvdrift = cvdrift * ntor * 2. / bmag**3    ! not correct for beta_prime term, so use: 
    cvdrift = gbdrift 

!    alp = -1./pres*dpdpsi

!
! these are theta and mod(b) grids that are periodic
!
    thetain = thetagrid(-nthetain/2:nthetain/2)
    bmagin = bmag(-nthetain/2:nthetain/2)
  end subroutine get_initial_grids_3d

  !> FIXME : Add documentation  
  subroutine do_smoothing
    implicit none
    real :: var
    integer :: ifail
    if (smoothknob == 0.0) then
       bmagsm = bmagin
    else
       var = smoothknob
       ifail = 0
       call smooth (nthetain+1,thetain,bmagin,var,bmagsm,ifail)
       if (ifail /= 0) then
          print *, "smooth failed: ",ifail
          select case (ifail)
          case (129)
             print *, "IC IS LESS THAN N-1."
          case (130)
             print *, "N IS lESS THAN 3."
          case (131)
             print *, "INPUT ABSCISSAE ARE NOT ORDERED SO THAT X(I).LT.X(I+1)."
          case (132)
             print *, "DF(I) IS NOT POSITIVE FOR SOME I."
          case (133)
             print *, "JOB IS NOT 0 OR 1."
          case default
          end select
          stop
       end if
    end if
  end subroutine do_smoothing

  !> FIXME : Add documentation  
  subroutine generate_grids
    implicit none
    integer :: i, iter, nin
    real :: d, deltaw_ok

    if (.not. auto_width) then
       ntheta = nthetaout
       nlambda = nlambdaout
       call gridgen4_1 (iperiod,nthetain+1,thetain,bmagsm, npadd, &
            alknob,epsknob,bpknob,extrknob,thetamax,deltaw,widthw,tension, &
            ntheta,nlambda,thetaout,bmagout,alambdaout)
       return
    end if

    widthw = pi
    do i = 0, nthetain/2
       if (cvdrift(i) < 0.0) then
          widthw = thetagrid(i)
          exit
       end if
    end do

    print *, "widthw: ", widthw
    deltaw_ok = 0.0
    do iter = 1, max_autoiter
       ntheta = nthetaout
       nlambda = nlambdaout
       call gridgen4_1 (iperiod,nthetain+1,thetain,bmagsm, npadd, &
            alknob,epsknob,bpknob,extrknob,thetamax,deltaw,widthw,tension, &
            ntheta,nlambda,thetaout,bmagout,alambdaout)
       print *, "iter: ", iter
       d = maxval(thetaout(2:ntheta) - thetaout(1:ntheta-1))
       print *, "max deltheta: ", d
       nin = count(abs(thetaout(1:ntheta) - thetamax) < widthw &
                   .or. abs(thetaout(1:ntheta) + thetamax) < widthw)
       print *, "count in +ve cvdrift: ", nin
       print *, "fraction in +ve cvdrift: ", real(nin)/real(ntheta)
       print *, "deltaw: ", deltaw
       if (d > delth_max) then
          if (deltaw_ok /= 0.0) then
             deltaw = sqrt(deltaw*deltaw_ok)
          else
             deltaw = 0.5*deltaw
          end if
          print *, "max deltheta too large: ", d
          if (deltaw < deltaw_ok) exit
       else if (real(nin)/real(ntheta) < cv_fraction) then
          print *, "fraction in +ve cvdrift too small: ", &
               real(nin)/real(ntheta)
          deltaw_ok = deltaw
          deltaw = deltaw*(cv_fraction/(real(nin)/real(ntheta)))**2
       end if
    end do

    if (deltaw_ok /= 0.0 .and. d > delth_max) then
       deltaw = deltaw_ok
       ntheta = nthetaout
       nlambda = nlambdaout
       call gridgen4_1 (iperiod,nthetain+1,thetain,bmagsm, npadd, &
            alknob,epsknob,bpknob,extrknob,thetamax,deltaw,widthw,tension, &
            ntheta,nlambda,thetaout,bmagout,alambdaout)
       d = maxval(thetaout(2:ntheta) - thetaout(1:ntheta-1))
       print *, "max deltheta: ", d
       nin = count(abs(thetaout(1:ntheta) - thetamax) < widthw &
                   .or. abs(thetaout(1:ntheta) + thetamax) < widthw)
       print *, "count in +ve cvdrift: ", nin
       print *, "fraction in +ve cvdrift: ", real(nin)/real(ntheta)
       print *, "deltaw: ", deltaw
    end if

  end subroutine generate_grids

  !> FIXME : Add documentation    
  subroutine write_output
    use splines, only: inter_cspl
    use gs2_io_grids, only: nc_write_grid_sizes, nc_write_grid_scalars, nc_write_grids
    use gs2_io_grids, only: nc_write_lambda_grid
    use gs2_io_grids, only: nc_grid_file_open, nc_grid_file_close
    use constants, only: run_name_size
    use mp, only: proc0
    implicit none
    integer :: unit
    integer :: i
!    real, allocatable, dimension (:) :: bmaginaux, bmagsmaux, tmp
    real, allocatable, dimension (:) :: thetagridout
    real, allocatable, dimension (:) :: gbdriftout, gradparout, grhoout
    real, allocatable, dimension (:) :: cvdriftout, gds2out, bmaggridout
    real, allocatable, dimension (:) :: gds21out, gds22out
    real, allocatable, dimension (:) :: cvdrift0out, gbdrift0out
    real, allocatable, dimension (:) :: Rplotout, Rprimeout
    real, allocatable, dimension (:) :: Zplotout, Zprimeout
    real, allocatable, dimension (:) :: aplotout, aprimeout
!    real :: th, bmin
    real :: bmin
!    integer :: ierr
    integer, dimension (1) :: minloca
    character(len=run_name_size) :: ncout

    call open_output_file (unit, ".input.out")
    write (unit=unit, fmt="('#',i5)") nthetain+1
    do i = 1, nthetain+1
       write (unit=unit, fmt="(3(1x,g19.12))") thetain(i), bmagin(i), bmagsm(i)
    enddo
    call close_output_file (unit)

!    allocate (bmaginaux(nthetain+1), bmagsmaux(nthetain+1), tmp(2*nthetain))
!    call fitp_curvp1 (nthetain,thetain,bmagin,twopi,bmaginaux,tmp,1.0,ierr)
!    call fitp_curvp1 (nthetain,thetain,bmagsm,twopi,bmagsmaux,tmp,1.0,ierr)
!    call open_output_file (unit, ".input.fine")
!    do i = -nfinegrid, nfinegrid
!       th = pi*real(i)/real(nfinegrid)
!       write (unit=unit, fmt="(3(x,g19.12))") th, &
!            fitp_curvp2(th,nthetain,thetain,bmagin,twopi,bmaginaux,1.0), &
!            fitp_curvp2(th,nthetain,thetain,bmagsm,twopi,bmagsmaux,1.0)
!    end do
!    call close_output_file (unit)
!    deallocate (bmaginaux, bmagsmaux, tmp)

    call open_output_file (unit, ".bmag.out")
    write (unit=unit, fmt="('#',i5)") ntheta
    do i = 1, ntheta
       write (unit=unit, fmt="(2(1x,g19.12))") thetaout(i), bmagout(i)
    enddo
    call close_output_file (unit)

    call open_output_file (unit, ".lambda.out")
    write (unit=unit, fmt="('#',i5)") nlambda
    do i = 1, nlambda
       write (unit=unit, fmt=*) alambdaout(i)
    enddo
    call close_output_file (unit)

    ntgrid = ntheta/2 + (nperiodout-1)*ntheta
    allocate (thetagridout(-ntgrid:ntgrid))
    allocate (gbdriftout(-ntgrid:ntgrid))
    allocate (gradparout(-ntgrid:ntgrid))
    allocate (grhoout(-ntgrid:ntgrid))
    allocate (cvdriftout(-ntgrid:ntgrid))
    allocate (gds2out(-ntgrid:ntgrid))
    allocate (bmaggridout(-ntgrid:ntgrid))
    allocate (gds21out(-ntgrid:ntgrid))
    allocate (gds22out(-ntgrid:ntgrid))
    allocate (cvdrift0out(-ntgrid:ntgrid))
    allocate (gbdrift0out(-ntgrid:ntgrid))
    allocate (Rplotout(-ntgrid:ntgrid))
    allocate (Rprimeout(-ntgrid:ntgrid))
    allocate (Zplotout(-ntgrid:ntgrid))
    allocate (Zprimeout(-ntgrid:ntgrid))
    allocate (aplotout(-ntgrid:ntgrid))
    allocate (aprimeout(-ntgrid:ntgrid))

    thetagridout(-ntheta/2:ntheta/2-1) = thetaout(:ntheta)
    bmaggridout(-ntheta/2:ntheta/2-1) = bmagout(:ntheta)
    thetagridout(ntheta/2) = thetagridout(-ntheta/2) + twopi*iperiod
    bmaggridout(ntheta/2) = bmaggridout(-ntheta/2)
    do i = 1, nperiodout - 1
       thetagridout(-ntheta/2-ntheta*i:ntheta/2-ntheta*i-1) &
            = thetagridout(-ntheta/2:ntheta/2-1) - twopi*i
       thetagridout(-ntheta/2+ntheta*i+1:ntheta/2+ntheta*i) &
            = thetagridout(-ntheta/2+1:ntheta/2) + twopi*i
       bmaggridout(-ntheta/2-ntheta*i:ntheta/2-ntheta*i-1) &
            = bmaggridout(-ntheta/2:ntheta/2-1)
       bmaggridout(-ntheta/2+ntheta*i+1:ntheta/2+ntheta*i) &
            = bmaggridout(-ntheta/2+1:ntheta/2)
    end do

    call inter_cspl(thetagrid, gbdrift, thetagridout, gbdriftout)
    call inter_cspl(thetagrid, gradpar, thetagridout, gradparout)
    call inter_cspl(thetagrid, grho, thetagridout, grhoout)
    call inter_cspl(thetagrid, cvdrift, thetagridout, cvdriftout)
    call inter_cspl(thetagrid, gds2, thetagridout, gds2out)
    call inter_cspl(thetagrid, bmag, thetagridout, bmagout)
    call inter_cspl(thetagrid, gds21, thetagridout, gds21out)
    call inter_cspl(thetagrid, gds22, thetagridout, gds22out)
    call inter_cspl(thetagrid, cvdrift0, thetagridout, cvdrift0out)
    call inter_cspl(thetagrid, gbdrift0, thetagridout, gbdrift0out)
    call inter_cspl(thetagrid, Rplot, thetagridout, Rplotout)
    call inter_cspl(thetagrid, Rprime, thetagridout, Rprimeout)
    call inter_cspl(thetagrid, Zplot, thetagridout, Zplotout)
    call inter_cspl(thetagrid, Zprime, thetagridout, Zprimeout)
    call inter_cspl(thetagrid, aplot, thetagridout, aplotout)
    call inter_cspl(thetagrid, aprime, thetagridout, aprimeout)

    call open_output_file (unit, ".out")
    write (unit=unit, fmt=*) 'nlambda'
    write (unit=unit, fmt=*) nlambda
    write (unit=unit, fmt=*) 'lambda'
    do i = 1, nlambda
       write (unit=unit, fmt=*) alambdaout(i)
    enddo

    write (unit=unit, fmt=*) 'ntgrid nperiod ntheta drhodpsi rmaj shat kxfac q'
    write (unit=unit, fmt="(3i4,5(1x,g19.10))") ntgrid, nperiodout, ntheta, drhodpsi, rmaj, shat, kxfac, qval

    write (unit=unit, fmt=*) 'gbdrift gradpar grho tgrid'
    do i = -ntgrid,ntgrid
       write (unit=unit, fmt="(8(1x,g19.10))") &
            gbdriftout(i),gradparout(i),grhoout(i),thetagridout(i)
    end do

    write (unit=unit, fmt=*) 'cvdrift gds2 bmag tgrid'
    do i = -ntgrid,ntgrid
       write (unit=unit, fmt="(8(1x,g19.10))") &
            cvdriftout(i),gds2out(i),bmaggridout(i),thetagridout(i)
    end do

    write (unit=unit, fmt=*) 'gds21 gds22 tgrid'
    do i = -ntgrid,ntgrid
       write (unit=unit, fmt="(8(1x,g19.10))") &
            gds21out(i),gds22out(i),thetagridout(i)
    end do

    write (unit=unit, fmt=*) 'cvdrift0 gbdrift0 tgrid'
    do i = -ntgrid,ntgrid
       write (unit=unit, fmt="(8(1x,g19.10))") &
            cvdrift0out(i),gbdrift0out(i),thetagridout(i)
    end do

    write (unit=unit, fmt=*) 'Rplot Rprime tgrid'
    do i = -ntgrid,ntgrid
       write (unit=unit, fmt="(8(1x,g19.10))") &
            Rplotout(i),Rprimeout(i),thetagridout(i)
    end do

    write (unit=unit, fmt=*) 'Zplot Zprime tgrid'
    do i = -ntgrid,ntgrid
       write (unit=unit, fmt="(8(1x,g19.10))") &
            Zplotout(i),Zprimeout(i),thetagridout(i)
    end do

    write (unit=unit, fmt=*) 'aplot aprime tgrid'
    do i = -ntgrid,ntgrid
       write (unit=unit, fmt="(8(1x,g19.10))") &
            aplotout(i),aprimeout(i),thetagridout(i)
    end do

    call close_output_file (unit)

    if (screenout) then
       write (*, *) 'cvdrift       gds2          bmag          theta'
       do i = -ntheta/2,ntheta/2
          write (unit=*, fmt="(4(1x,g13.5))") &
               cvdriftout(i),gds2out(i),bmaggridout(i),thetagridout(i)
       end do
    end if

    minloca = minloc(bmagout(:ntheta))
    bmin = bmagout(minloca(1))
    print *, "theta,bmag minimum: ", thetaout(minloca(1)), bmin

    minloca = minloc(bmagin)
    print *, "theta,bmag input minimum: ", &
         thetain(minloca(1)), bmagin(minloca(1))

    if (bmin < bmagin(minloca(1))) print *, "WARNING: interpolated new minimum"

    if (proc0) then
       ncout = trim(run_name)//".out.nc"
       call nc_grid_file_open(ncout, "w")
       call nc_write_grid_sizes(ntheta, ntgrid, nperiodout)
       call nc_write_grid_scalars(shat, drhodpsi, kxfac, qval, rmaj)
       call nc_write_grids(ntgrid, bmaggridout, gradparout, gbdriftout, gbdrift0out, &
            cvdriftout, cvdrift0out, gds2out, gds21out, gds22out, grhoout, thetagridout, &
            Rplot=Rplotout, Rprime=Rprimeout, Zplot=Zplotout, Zprime=Zprimeout, &
            aplot=aplotout, aprime=aprimeout, no_end_point_in=no_end_point)
       call nc_write_lambda_grid(nlambda, alambdaout(:nlambda))
       call nc_grid_file_close()
    end if

  end subroutine write_output

  !> FIXME : Add documentation
  subroutine smooth (n, xin, yin, var, yout, ifail)
    implicit none

    integer, intent(in) :: n
    real, dimension (:), intent (in) :: xin, yin
    real, dimension (:), intent (out) :: yout
    real, intent(in out) :: var
    integer, intent(out) :: ifail

! these next arrays should be double precision, which we usually get with
! compiler options
    real, dimension (1) :: se
    real, dimension (:), allocatable :: x, f, y, df
    real, dimension (:,:), allocatable :: wk, c
    real :: dvar

    allocate (x(n), f(n), y(n), df(n), c(n, 3))
    allocate (wk(n+2,7))

    ifail = 0
    x = xin
    f = yin
    df = 1.
    dvar=var
    call cubgcv (x,f,df,n,y,c,n,dvar,0,se,wk,ifail)
    var=dvar
    yout = y

    deallocate (x, f, y, df, c, wk)

  end subroutine smooth

!     algorithm 642 collected algorithms from acm.
!     algorithm appeared in acm-trans. math. software, vol.12, no. 2,
!     jun., 1986, p. 150.
!   subroutine name     - cubgcv
!
!--------------------------------------------------------------------------
!
!   computer            - vax/double
!
!   author              - m.f.hutchinson
!                         csiro division of mathematics and statistics
!                         p.o. box 1965
!                         canberra, act 2601
!                         australia
!
!   latest revision     - 15 august 1985
!
!   purpose             - cubic spline data smoother
!
!   usage               - call cubgcv (x,f,df,n,y,c,ic,var,job,se,wk,ier)
!
!   arguments    x      - vector of length n containing the
!                           abscissae of the n data points
!                           (x(i),f(i)) i=1..n. (input) x
!                           must be ordered so that
!                           x(i) .lt. x(i+1).
!                f      - vector of length n containing the
!                           ordinates (or function values)
!                           of the n data points (input).
!                df     - vector of length n. (input/output)
!                           df(i) is the relative standard deviation
!                           of the error associated with data point i.
!                           each df(i) must be positive.  the values in
!                           df are scaled by the subroutine so that
!                           their mean square value is 1, and unscaled
!                           again on normal exit.
!                           the mean square value of the df(i) is returned
!                           in wk(7) on normal exit.
!                           if the absolute standard deviations are known,
!                           these should be provided in df and the error
!                           variance parameter var (see below) should then
!                           be set to 1.
!                           if the relative standard deviations are unknown,
!                           set each df(i)=1.
!                n      - number of data points (input).
!                           n must be .ge. 3.
!                y,c    - spline coefficients. (output) y
!                           is a vector of length n. c is
!                           an n-1 by 3 matrix. the value
!                           of the spline approximation at t is
!                           s(t)=((c(i,3)*d+c(i,2))*d+c(i,1))*d+y(i)
!                           where x(i).le.t.lt.x(i+1) and
!                           d = t-x(i).
!                ic     - row dimension of matrix c exactly
!                           as specified in the dimension
!                           statement in the calling program. (input)
!                var    - error variance. (input/output)
!                           if var is negative (i.e. unknown) then
!                           the smoothing parameter is determined
!                           by minimizing the generalized cross validation
!                           and an estimate of the error variance is
!                           returned in var.
!                           if var is non-negative (i.e. known) then the
!                           smoothing parameter is determined to minimize
!                           an estimate, which depends on var, of the true
!                           mean square error, and var is unchanged.
!                           in particular, if var is zero, then an
!                           interpolating natural cubic spline is calculated.
!                           var should be set to 1 if absolute standard
!                           deviations have been provided in df (see above).
!                job    - job selection parameter. (input)
!                         job = 0 should be selected if point standard error
!                           estimates are not required in se.
!                         job = 1 should be selected if point standard error
!                           estimates are required in se.
!                se     - vector of length n containing bayesian standard
!                           error estimates of the fitted spline values in y.
!                           se is not referenced if job=0. (output)
!                wk     - work vector of length 7*(n + 2). on normal exit the
!                           first 7 values of wk are assigned as follows:-
!
!                           wk(1) = smoothing parameter (= rho/(rho + 1))
!                           wk(2) = estimate of the number of degrees of
!                                   freedom of the residual sum of squares
!                           wk(3) = generalized cross validation
!                           wk(4) = mean square residual
!                           wk(5) = estimate of the true mean square error
!                                   at the data points
!                           wk(6) = estimate of the error variance
!                           wk(7) = mean square value of the df(i)
!
!                           if wk(1)=0 (rho=0) an interpolating natural cubic
!                           spline has been calculated.
!                           if wk(1)=1 (rho=infinite) a least squares
!                           regression line has been calculated.
!                           wk(2) is an estimate of the number of degrees of
!                           freedom of the residual which reduces to the
!                           usual value of n-2 when a least squares regression
!                           line is calculated.
!                           wk(3),wk(4),wk(5) are calculated with the df(i)
!                           scaled to have mean square value 1.  the
!                           unscaled values of wk(3),wk(4),wk(5) may be
!                           calculated by dividing by wk(7).
!                           wk(6) coincides with the output value of var if
!                           var is negative on input.  it is calculated with
!                           the unscaled values of the df(i) to facilitate
!                           comparisons with a priori variance estimates.
!
!                ier    - error parameter. (output)
!                         terminal error
!                           ier = 129, ic is less than n-1.
!                           ier = 130, n is less than 3.
!                           ier = 131, input abscissae are not
!                             ordered so that x(i).lt.x(i+1).
!                           ier = 132, df(i) is not positive for some i.
!                           ier = 133, job is not 0 or 1.
!
!   precision/hardware  - double
!
!   required routines   - spint1,spfit1,spcof1,sperr1
!
!   remarks      the number of arithmetic operations required by the
!                subroutine is proportional to n.  the subroutine
!                uses an algorithm developed by m.f. hutchinson and
!                f.r. de hoog, 'smoothing noisy data with spline
!                functions', numer. math. (in press)
!
!-----------------------------------------------------------------------
!
  !> FIXME : Add documentation
  subroutine cubgcv(x,f,df,n,y,c,ic,var,job,se,wk,ier)
!
!---specifications for arguments---
    integer n,ic,job,ier
!
! All real variables should be double precision:
!
    real, dimension (:) :: x, f, df, y, se
    real, dimension (:,:) :: c
    real, dimension (0:n+1,7) :: wk
    real :: var
!
!---specifications for local variables---
    real :: delta,err,gf1,gf2,gf3,gf4,r1,r2,r3,r4,avh,avdf,avar,stat(6),p,q
    real :: ratio = 2.
    real :: tau = 1.618033989
    real :: zero = 0.
    real :: one = 1.

    integer :: i
!
!---initialize---
    ier = 133
    if (job.lt.0 .or. job.gt.1) go to 140
    call spint1(x,avh,f,df,avdf,n,y,c,ic,wk,wk(0,4),ier)
    if (ier.ne.0) go to 140
    avar = var
    if (var.gt.zero) avar = var*avdf*avdf
!
!---check for zero variance---
    if (var.ne.zero) go to 10
    r1 = zero
    go to 90
!
!---find local minimum of gcv or the expected mean square error---
10  r1 = one
    r2 = ratio*r1
    call spfit1(x,avh,df,n,r2,p,q,gf2,avar,stat,y,c,ic,wk(:,1:3),wk(:,4:5),wk(:,6),wk(:,7))
20  call spfit1(x,avh,df,n,r1,p,q,gf1,avar,stat,y,c,ic,wk(:,1:3),wk(:,4:5),wk(:,6),wk(:,7))
    if (gf1.gt.gf2) go to 30
!
!---exit if p zero---
    if (p.le.zero) go to 100
    r2 = r1
    gf2 = gf1
    r1 = r1/ratio
    go to 20

30  r3 = ratio*r2
40  call spfit1(x,avh,df,n,r3,p,q,gf3,avar,stat,y,c,ic,wk(:,1:3),wk(:,4:5),wk(:,6),wk(:,7))
    if (gf3.gt.gf2) go to 50
!
!---exit if q zero---
    if (q.le.zero) go to 100
    r2 = r3
    gf2 = gf3
    r3 = ratio*r3
    go to 40

50  r2 = r3
    gf2 = gf3
    delta = (r2-r1)/tau
    r4 = r1 + delta
    r3 = r2 - delta
    call spfit1(x,avh,df,n,r3,p,q,gf3,avar,stat,y,c,ic,wk(:,1:3),wk(:,4:5),wk(:,6),wk(:,7))
    call spfit1(x,avh,df,n,r4,p,q,gf4,avar,stat,y,c,ic,wk(:,1:3),wk(:,4:5),wk(:,6),wk(:,7))
!
!---golden section search for local minimum---
60  if (gf3.gt.gf4) go to 70
    r2 = r4
    gf2 = gf4
    r4 = r3
    gf4 = gf3
    delta = delta/tau
    r3 = r2 - delta
    call spfit1(x,avh,df,n,r3,p,q,gf3,avar,stat,y,c,ic,wk(:,1:3),wk(:,4:5),wk(:,6),wk(:,7))
    go to 80

70  r1 = r3
    gf1 = gf3
    r3 = r4
    gf3 = gf4
    delta = delta/tau
    r4 = r1 + delta
    call spfit1(x,avh,df,n,r4,p,q,gf4,avar,stat,y,c,ic,wk(:,1:3),wk(:,4:5),wk(:,6),wk(:,7))
80  err = (r2-r1)/ (r1+r2)
    if (err*err+one.gt.one .and. err.gt.1.0e-6) go to 60
    r1 = (r1+r2)*0.5
!
!---calculate spline coefficients---
90  call spfit1(x,avh,df,n,r1,p,q,gf1,avar,stat,y,c,ic,wk(:,1:3),wk(:,4:5),wk(:,6),wk(:,7))
100 call spcof1(x,avh,f,df,n,p,q,y,c,ic,wk(:,6),wk(:,7))
!
!---optionally calculate standard error estimates---
    if (var.ge.zero) go to 110
    avar = stat(6)
    var = avar/ (avdf*avdf)
110 if (job.eq.1) call sperr1(x,avh,df,n,wk,p,avar,se)
!
!---unscale df---
    df = df * avdf
!
!--put statistics in wk---
    do i = 0,5
       wk(i,1) = stat(i+1)
    end do
    wk(5,1) = stat(6)/ (avdf*avdf)
    wk(6,1) = avdf*avdf
    go to 150
!
!---check for error condition---
140 continue
!     if (ier.ne.0) continue
150 return
  end subroutine cubgcv

  !> FIXME : Add documentation
  subroutine spint1(x,avh,y,dy,avdy,n,a,c,ic,r,t,ier)
!
! initializes the arrays c, r and t for one dimensional cubic
! smoothing spline fitting by subroutine spfit1.  the values
! df(i) are scaled so that the sum of their squares is n
! and the average of the differences x(i+1) - x(i) is calculated
! in avh in order to avoid underflow and overflow problems in
! spfit1.
!
! subroutine sets ier if elements of x are non-increasing,
! if n is less than 3, if ic is less than n-1 or if dy(i) is
! not positive for some i.
!
!---specifications for arguments---
    integer n,ic,ier
! All variables should be double precision
    real, dimension (:) :: x, y, dy, a
    real, dimension (:,:) :: c
    real, dimension (0:n+1,3) :: r
    real, dimension (0:n+1,2) :: t
    real :: avh, avdy

!
!---specifications for local variables---
    integer i
    real :: e,f,g,h
    real :: zero = 0.

!
!---initialization and input checking---
    ier = 0
    if (n.lt.3) go to 60
    if (ic.lt.n-1) go to 70
!
!---get average x spacing in avh---
    g = zero
    do i = 1,n - 1
       h = x(i+1) - x(i)
       if (h.le.zero) go to 80
       g = g + h
    end do
    avh = g/ (n-1)
!
!---scale relative weights---
    g = zero
    do i = 1,n
       if (dy(i).le.zero) go to 90
       g = g + dy(i)*dy(i)
    end do
    avdy = sqrt(g/n)

    dy = dy / avdy

!
!---initialize h,f---
    h = (x(2)-x(1))/avh
    f = (y(2)-y(1))/h
!
!---calculate a,t,r---
    do i = 2,n - 1
       g = h
       h = (x(i+1)-x(i))/avh
       e = f
       f = (y(i+1)-y(i))/h
       a(i) = f - e
       t(i,1) = 2.0d0* (g+h)/3.0d0
       t(i,2) = h/3.0d0
       r(i,3) = dy(i-1)/g
       r(i,1) = dy(i+1)/h
       r(i,2) = -dy(i)/g - dy(i)/h
    end do
!
!---calculate c = r'*r---
    r(n,2) = zero
    r(n,3) = zero
    r(n+1,3) = zero
    do i = 2,n - 1
       c(i,1) = r(i,1)*r(i,1) + r(i,2)*r(i,2) + r(i,3)*r(i,3)
       c(i,2) = r(i,1)*r(i+1,2) + r(i,2)*r(i+1,3)
       c(i,3) = r(i,1)*r(i+2,3)
    end do
    return
!
!---error conditions---
60  ier = 130
    return

70  ier = 129
    return

80  ier = 131
    return

   90 ier = 132
    return
  end subroutine spint1

  !> FIXME : Add documentation
  subroutine spfit1(x,avh,dy,n,rho,p,q,fun,var,stat,a,c,ic,r,t,u,v)
!
! fits a cubic smoothing spline to data with relative
! weighting dy for a given value of the smoothing parameter
! rho using an algorithm based on that of c.h. reinsch (1967),
! numer. math. 10, 177-183.
!
! the trace of the influence matrix is calculated using an
! algorithm developed by m.f.hutchinson and f.r.de hoog (numer.
! math., in press), enabling the generalized cross validation
! and related statistics to be calculated in order n operations.
!
! the arrays a, c, r and t are assumed to have been initialized
! by the subroutine spint1.  overflow and underflow problems are
! avoided by using p=rho/(1 + rho) and q=1/(1 + rho) instead of
! rho and by scaling the differences x(i+1) - x(i) by avh.
!
! the values in df are assumed to have been scaled so that the
! sum of their squared values is n.  the value in var, when it is
! non-negative, is assumed to have been scaled to compensate for
! the scaling of the values in df.
!
! the value returned in fun is an estimate of the true mean square
! when var is non-negative, and is the generalized cross validation
! when var is negative.
!
!---specifications for arguments---
    integer ic,n
! all variables should be double precision:
    real, dimension (:) :: x, dy, a
    real, dimension (6) :: stat
    real, dimension (ic,3) :: c
    real, dimension (0:n+1,3) :: r
    real, dimension (0:n+1,2) :: t
    real, dimension (0:n+1) :: u, v
    real :: fun, var, avh, p, q, rho
!
!---local variables---
    integer i
    real :: e,f,g,h,rho1
    real :: zero = 0.
    real :: one = 1.
    real :: two = 2.

!
!---use p and q instead of rho to prevent overflow or underflow---
    rho1 = one + rho
    p = rho/rho1
    q = one/rho1
    if (rho1.eq.one) p = zero
    if (rho1.eq.rho) q = zero
!
!---rational cholesky decomposition of p*c + q*t---
    f = zero
    g = zero
    h = zero
    do i = 0,1
       r(i,1) = zero
    end do
    do i = 2,n - 1
       r(i-2,3) = g*r(i-2,1)
       r(i-1,2) = f*r(i-1,1)
       r(i,1) = one/ (p*c(i,1)+q*t(i,1)-f*r(i-1,2)-g*r(i-2,3))
       f = p*c(i,2) + q*t(i,2) - h*r(i-1,2)
       g = h
       h = p*c(i,3)
    end do
!
!---solve for u---
    u(0) = zero
    u(1) = zero
    do i = 2,n - 1
       u(i) = a(i) - r(i-1,2)*u(i-1) - r(i-2,3)*u(i-2)
    end do
    u(n) = zero
    u(n+1) = zero
    do i = n - 1,2,-1
       u(i) = r(i,1)*u(i) - r(i,2)*u(i+1) - r(i,3)*u(i+2)
    end do
!
!---calculate residual vector v---
    e = zero
    h = zero
    do i = 1,n - 1
       g = h
       h = (u(i+1)-u(i))/ ((x(i+1)-x(i))/avh)
       v(i) = dy(i)* (h-g)
       e = e + v(i)*v(i)
    end do
    v(n) = dy(n)* (-h)
    e = e + v(n)*v(n)
!
!---calculate upper three bands of inverse matrix---
    r(n,1) = zero
    r(n,2) = zero
    r(n+1,1) = zero
    do i = n - 1,2,-1
       g = r(i,2)
       h = r(i,3)
       r(i,2) = -g*r(i+1,1) - h*r(i+1,2)
       r(i,3) = -g*r(i+1,2) - h*r(i+2,1)
       r(i,1) = r(i,1) - g*r(i,2) - h*r(i,3)
    end do
!
!---calculate trace---
    f = zero
    g = zero
    h = zero
    do i = 2,n - 1
       f = f + r(i,1)*c(i,1)
       g = g + r(i,2)*c(i,2)
       h = h + r(i,3)*c(i,3)
    end do
    f = f + two* (g+h)
!
!---calculate statistics---
    stat(1) = p
    stat(2) = f*p
    stat(3) = n*e/ (f*f)
    stat(4) = e*p*p/n
    stat(6) = e*p/f
    if (var.ge.zero) go to 80
    stat(5) = stat(6) - stat(4)
    fun = stat(3)
    go to 90

80  stat(5) = amax1(stat(4)-two*var*stat(2)/n+var,zero)
    fun = stat(5)
90  return
  end subroutine spfit1

  !> FIXME : Add documentation
  subroutine sperr1(x,avh,dy,n,r,p,var,se)
!
! calculates bayesian estimates of the standard errors of the fitted
! values of a cubic smoothing spline by calculating the diagonal elements
! of the influence matrix.
!
!---specifications for arguments---
    integer n
    real, dimension (:) :: x, dy, se
    real, dimension (0:n+1, 3) :: r
    real :: avh, p, var
!
!---specifications for local variables---
    integer i
    real :: f,g,h,f1,g1,h1
    real :: zero = 0.
    real :: one = 1.
!
!---initialize---
    h = avh/ (x(2)-x(1))
    se(1) = one - p*dy(1)*dy(1)*h*h*r(2,1)
    r(1,1) = zero
    r(1,2) = zero
    r(1,3) = zero
!
!---calculate diagonal elements---
    do i = 2,n - 1
       f = h
       h = avh/ (x(i+1)-x(i))
       g = -f - h
       f1 = f*r(i-1,1) + g*r(i-1,2) + h*r(i-1,3)
       g1 = f*r(i-1,2) + g*r(i,1) + h*r(i,2)
       h1 = f*r(i-1,3) + g*r(i,2) + h*r(i+1,1)
       se(i) = one - p*dy(i)*dy(i)* (f*f1+g*g1+h*h1)
    end do
    se(n) = one - p*dy(n)*dy(n)*h*h*r(n-1,1)
!
!---calculate standard error estimates---
    do i = 1,n
       se(i) = sqrt(amax1(se(i)*var,zero))*dy(i)
    end do

  end subroutine sperr1

  !> FIXME : Add documentation
  subroutine spcof1(x,avh,y,dy,n,p,q,a,c,ic,u,v)
!
! calculates coefficients of a cubic smoothing spline from
! parameters calculated by subroutine spfit1.
!
!---specifications for arguments---
    integer ic,n
    real, dimension (:) :: x, y, dy, a
    real, dimension (ic, 3) :: c
    real, dimension (0:n+1) :: u, v
    real :: p, q, avh
!
!---specifications for local variables---
    integer i
    real :: h,qh
!
!---calculate a---
    qh = q/ (avh*avh)
    do i = 1,n
       a(i) = y(i) - p*dy(i)*v(i)
       u(i) = qh*u(i)
    end do
!
!---calculate c---
    do i = 1,n - 1
       h = x(i+1) - x(i)
       c(i,3) = (u(i+1)-u(i))/ (3.0d0*h)
       c(i,1) = (a(i+1)-a(i))/h - (h*c(i,3)+u(i))*h
       c(i,2) = u(i)
    end do
  end subroutine spcof1

!---c cubgcv test driver
!---c ------------------
!---c
!---c author          - m.f.hutchinson
!---c                   csiro division of water and land resources
!---c                   gpo box 1666
!---c                   canberra act 2601
!---c                   australia
!---c
!---c latest revision - 7 august 1986
!---c
!---c computer        - vax/double
!---c
!---c usage           - main program
!---c
!---c required routines - cubgcv,spint1,spfit1,spcof1,sperr1,ggrand
!---c
!---c remarks   uses subroutine cubgcv to fit a cubic smoothing spline
!---c           to 50 data points which are generated by adding a random
!---c           variable with uniform density in the interval [-0.3,0.3]
!---c           to 50 points sampled from the curve  y=sin(3*pi*x/2).
!---c           random deviates in the interval [0,1] are generated by the
!---c           double precision function ggrand (similar to imsl function
!---c           ggubfs).  the abscissae are unequally spaced in [0,1].
!---c
!---c           point standard error estimates are returned in se by
!---c           setting job=1.  the error variance estimate is returned
!---c           in var.  it compares favourably with the true value of 0.03.
!---c           summary statistics from the array wk are written to
!---c           unit 6.  data values and fitted values with estimated
!---c           standard errors are also written to unit 6.
!---c
!---      parameter (n=50, ic=49)
!---c
!---      integer            job,ier
!---      double precision   x(n),f(n),y(n),df(n),c(ic,3),wk(7*(n+2)),
!---     *                   var,se(n)
!---      double precision   ggrand,dseed
!---c
!---c---initialize---
!---      dseed=1.2345d4
!---      job=1
!---      var=-1.0d0
!---c
!---c---calculate data points---
!---      do 10 i=1,n
!---      x(i)=(i - 0.5)/n + (2.0*ggrand(dseed) - 1.0)/(3.0*n)
!---      f(i)=dsin(4.71238*x(i)) + (2.0*ggrand(dseed) - 1.0)*0.3
!---      df(i)=1.0d0
!---  10  continue
!---c
!---c---fit cubic spline---
!---      call cubgcv(x,f,df,n,y,c,ic,var,job,se,wk,ier)
!---c
!---c---write out results---
!---      write(6,20)
!---  20  format(' cubgcv test driver results:')
!---      write(6,30) ier,var,wk(3),wk(4),wk(2)
!---  30  format(/' ier =',i4/' var =',f7.4/
!---     *        ' generalized cross validation =',f7.4/
!---     *        ' mean square residual         =',f7.4/
!---     *        ' residual degrees of freedom  =',f7.2)
!---      write(6,40)
!---  40  format(/' input data',17x,'output results'//
!---     *         '   i    x(i)    f(i)',6x,'    y(i)   se(i)',
!---     *          '      c(i,1)      c(i,2)      c(i,3)')
!---      do 60 i=1,n-1
!---      write(6,50) i,x(i),f(i),y(i),se(i),(c(i,j),j=1,3)
!---  50  format(i4,2f8.4,6x,2f8.4,3e12.4)
!---  60  continue
!---      write(6,50) n,x(n),f(n),y(n),se(n)
!---      stop
!---      end
!---      double precision function ggrand(dseed)
!---c
!---c double precision uniform random number generator
!---c
!---c constants: a = 7**5
!---c            b = 2**31 - 1
!---c            c = 2**31
!---c
!---c reference: imsl manual, chapter g - generation and testing of
!---c                                     random numbers
!---c
!---c---specifications for arguments---
!---      double precision dseed
!---c
!---c---specifications for local variables---
!---      double precision a,b,c,s
!---c
!---      data a,b,c/16807.0d0, 2147483647.0d0, 2147483648.0d0/
!---c
!---      s=dseed
!---      s=dmod(a*s, b)
!---      ggrand=s/c
!---      dseed=s
!---      return
!---      end
!---
!---cubgcv test driver results:
!---
!---ier =   0
!---var = 0.0279
!---generalized cross validation = 0.0318
!---mean square residual         = 0.0246
!---residual degrees of freedom  =  43.97
!---
!---input data                 output results
!---
!---  i    x(i)    f(i)          y(i)   se(i)      c(i,1)      c(i,2)      c(i,3)
!---  1  0.0046  0.2222        0.0342  0.1004  0.3630e+01  0.0000e+00  0.2542e+02
!---  2  0.0360 -0.1098        0.1488  0.0750  0.3705e+01  0.2391e+01 -0.9537e+01
!---  3  0.0435 -0.0658        0.1767  0.0707  0.3740e+01  0.2175e+01 -0.4233e+02
!---  4  0.0735  0.3906        0.2900  0.0594  0.3756e+01 -0.1642e+01 -0.2872e+02
!---  5  0.0955  0.6054        0.3714  0.0558  0.3642e+01 -0.3535e+01  0.2911e+01
!---  6  0.1078  0.3034        0.4155  0.0549  0.3557e+01 -0.3428e+01 -0.1225e+02
!---  7  0.1269  0.7386        0.4822  0.0544  0.3412e+01 -0.4131e+01  0.2242e+02
!---  8  0.1565  0.4616        0.5800  0.0543  0.3227e+01 -0.2143e+01  0.6415e+01
!---  9  0.1679  0.4315        0.6165  0.0543  0.3180e+01 -0.1923e+01 -0.1860e+02
!--- 10  0.1869  0.5716        0.6762  0.0544  0.3087e+01 -0.2985e+01 -0.3274e+02
!--- 11  0.2149  0.6736        0.7595  0.0542  0.2843e+01 -0.5733e+01 -0.4435e+02
!--- 12  0.2356  0.7388        0.8155  0.0539  0.2549e+01 -0.8486e+01 -0.5472e+02
!--- 13  0.2557  1.1953        0.8630  0.0537  0.2139e+01 -0.1180e+02 -0.9784e+01
!--- 14  0.2674  1.0299        0.8864  0.0536  0.1860e+01 -0.1214e+02  0.9619e+01
!--- 15  0.2902  0.7981        0.9225  0.0534  0.1322e+01 -0.1149e+02 -0.7202e+01
!--- 16  0.3155  0.8973        0.9485  0.0532  0.7269e+00 -0.1203e+02 -0.1412e+02
!--- 17  0.3364  1.2695        0.9583  0.0530  0.2040e+00 -0.1292e+02  0.2796e+02
!--- 18  0.3557  0.7253        0.9577  0.0527 -0.2638e+00 -0.1130e+02 -0.3453e+01
!--- 19  0.3756  1.2127        0.9479  0.0526 -0.7176e+00 -0.1151e+02  0.3235e+02
!--- 20  0.3881  0.7304        0.9373  0.0525 -0.9889e+00 -0.1030e+02  0.4381e+01
!--- 21  0.4126  0.9810        0.9069  0.0525 -0.1486e+01 -0.9977e+01  0.1440e+02
!--- 22  0.4266  0.7117        0.8842  0.0525 -0.1756e+01 -0.9373e+01 -0.8925e+01
!--- 23  0.4566  0.7203        0.8227  0.0524 -0.2344e+01 -0.1018e+02 -0.2278e+02
!--- 24  0.4704  0.9242        0.7884  0.0524 -0.2637e+01 -0.1112e+02 -0.4419e+01
!--- 25  0.4914  0.7345        0.7281  0.0523 -0.3110e+01 -0.1140e+02 -0.3562e+01
!--- 26  0.5084  0.7378        0.6720  0.0524 -0.3500e+01 -0.1158e+02  0.5336e+01
!--- 27  0.5277  0.7441        0.6002  0.0525 -0.3941e+01 -0.1127e+02  0.2479e+02
!--- 28  0.5450  0.5612        0.5286  0.0527 -0.4310e+01 -0.9980e+01  0.2920e+02
!--- 29  0.5641  0.5049        0.4429  0.0529 -0.4659e+01 -0.8309e+01  0.3758e+02
!--- 30  0.5857  0.4725        0.3390  0.0531 -0.4964e+01 -0.5878e+01  0.5563e+02
!--- 31  0.6159  0.1380        0.1850  0.0531 -0.5167e+01 -0.8307e+00  0.4928e+02
!--- 32  0.6317  0.1412        0.1036  0.0531 -0.5157e+01  0.1499e+01  0.5437e+02
!--- 33  0.6446 -0.1110        0.0371  0.0531 -0.5091e+01  0.3614e+01  0.3434e+02
!--- 34  0.6707 -0.2605       -0.0927  0.0532 -0.4832e+01  0.6302e+01  0.1164e+02
!--- 35  0.6853 -0.1284       -0.1619  0.0533 -0.4640e+01  0.6812e+01  0.1617e+02
!--- 36  0.7064 -0.3452       -0.2564  0.0536 -0.4332e+01  0.7834e+01  0.4164e+01
!--- 37  0.7310 -0.5527       -0.3582  0.0538 -0.3939e+01  0.8141e+01 -0.2214e+02
!--- 38  0.7531 -0.3459       -0.4415  0.0540 -0.3611e+01  0.6674e+01 -0.9205e+01
!--- 39  0.7686 -0.5902       -0.4961  0.0541 -0.3410e+01  0.6245e+01 -0.2193e+02
!--- 40  0.7952 -0.7644       -0.5828  0.0541 -0.3125e+01  0.4494e+01 -0.4649e+02
!--- 41  0.8087 -0.5392       -0.6242  0.0541 -0.3029e+01  0.2614e+01 -0.3499e+02
!--- 42  0.8352 -0.4247       -0.7031  0.0539 -0.2964e+01 -0.1603e+00  0.2646e+01
!--- 43  0.8501 -0.6327       -0.7476  0.0538 -0.2967e+01 -0.4132e-01  0.1817e+02
!--- 44  0.8726 -0.9983       -0.8139  0.0538 -0.2942e+01  0.1180e+01 -0.6774e+01
!--- 45  0.8874 -0.9082       -0.8574  0.0542 -0.2911e+01  0.8778e+00 -0.1364e+02
!--- 46  0.9139 -0.8930       -0.9340  0.0566 -0.2893e+01 -0.2044e+00 -0.8094e+01
!--- 47  0.9271 -1.0233       -0.9723  0.0593 -0.2903e+01 -0.5258e+00 -0.1498e+02
!--- 48  0.9473 -0.8839       -1.0313  0.0665 -0.2942e+01 -0.1433e+01  0.4945e+01
!--- 49  0.9652 -1.0172       -1.0843  0.0766 -0.2989e+01 -0.1168e+01  0.1401e+02
!--- 50  0.9930 -1.2715       -1.1679  0.0998
!---
!---
!---algorithm
!---
!---cubgcv calculates a natural cubic spline curve which smoothes a given set
!---of data points, using statistical considerations to determine the amount
!---of smoothing required, as described in reference 2.  if the error variance
!---is known, it should be supplied to the routine in var.  the degree of
!---smoothing is then determined by minimizing an unbiased estimate of the true
!---mean square error.  on the other hand, if the error variance is not known, var
!---should be set to -1.0.  the routine then determines the degree of smoothing
!---by minimizing the generalized cross validation.  this is asymptotically the
!---same as minimizing the true mean square error (see reference 1).  in this
!---case, an estimate of the error variance is returned in var which may be
!---compared with any a priori approximate estimates.  in either case, an
!---estimate of the true mean square error is returned in wk(5).  this estimate,
!---however, depends on the error variance estimate, and should only be accepted
!---if the error variance estimate is reckoned to be correct.
!---
!---if job is set to 1, bayesian estimates of the standard error of each
!---smoothed data value are returned in the array se.  these also depend on
!---the error variance estimate and should only be accepted if the error
!---variance estimate is reckoned to be correct.  see reference 4.
!---
!---the number of arithmetic operations and the amount of storage required by
!---the routine are both proportional to n, so that very large data sets may be
!---analysed.  the data points do not have to be equally spaced or uniformly
!---weighted.  the residual and the spline coefficients are calculated in the
!---manner described in reference 3, while the trace and various statistics,
!---including the generalized cross validation, are calculated in the manner
!---described in reference 2.
!---
!---when var is known, any value of n greater than 2 is acceptable.  it is
!---advisable, however, for n to be greater than about 20 when var is unknown.
!---if the degree of smoothing done by cubgcv when var is unknown is not
!---satisfactory, the user should try specifying the degree of smoothing by
!---setting var to a reasonable value.
!---
!---references:
!---
!---1.  craven, peter and wahba, grace, "smoothing noisy data with spline
!---    functions", numer. math. 31, 377-403 (1979).
!---2.  hutchinson, m.f. and de hoog, f.r., "smoothing noisy data with spline
!---    functions", numer. math. (in press).
!---3.  reinsch, c.h., "smoothing by spline functions", numer. math. 10,
!---    177-183 (1967).
!---4.  wahba, grace, "bayesian 'confidence intervals' for the cross-validated
!---    smoothing spline", j.r.statist. soc. b 45, 133-150 (1983).
!---
!---
!---example
!---
!---a sequence of 50 data points are generated by adding a random variable with
!---uniform density in the interval [-0.3,0.3] to the curve y=sin(3*pi*x/2).
!---the abscissae are unequally spaced in [0,1].  point standard error estimates
!---are returned in se by setting job to 1.  the error variance estimate is
!---returned in var.  it compares favourably with the true value of 0.03.
!---the imsl function ggubfs is used to generate sample values of a uniform
!---variable on [0,1].
!---
!---
!---input:
!---
!---      integer          n,ic,job,ier
!---      double precision x(50),f(50),y(50),df(50),c(49,3),wk(400),
!---     *                 var,se(50)
!---      double precision ggubfs,dseed,dn
!---      data dseed/1.2345d4/
!---c
!---      n=50
!---      ic=49
!---      job=1
!---      var=-1.0d0
!---      dn=n
!---c
!---      do 10 i=1,n
!---      x(i)=(i - 0.5)/dn + (2.0*ggubfs(dseed) - 1.0)/(3.0*dn)
!---      f(i)=dsin(4.71238*x(i)) + (2.0*ggubfs(dseed) - 1.0)*0.3
!---      df(i)=1.0d0
!---  10  continue
!---      call cubgcv(x,f,df,n,y,c,ic,var,job,se,wk,ier)
!---       .
!---       .
!---       .
!---      end
!---
!---output:
!---
!---ier = 0
!---var = 0.0279
!---generalized cross validation = 0.0318
!---mean square residual = 0.0246
!---residual degrees of freedom = 43.97
!---for checking purposes the following output is given:
!---
!---x(1)  = 0.0046    f(1)  =  0.2222     y(1)  =  0.0342     se(1)  = 0.1004
!---x(21) = 0.4126    f(21) =  0.9810     y(21) =  0.9069     se(21) = 0.0525
!---x(41) = 0.8087    f(41) = -0.5392     y(41) = -0.6242     se(41) = 0.0541
!---
end program rungridgen