read_neocoefs Subroutine

  subroutine read_neocoefs (theta, nspec_neo, nenergy_neo, nxi_neo, nrad_neo, ir_neo)

    use mp, only: proc0, broadcast
    use splines, only: lf_spline
    use file_utils, only: get_unused_unit
    use constants, only: pi,twopi

    implicit none

    real, dimension (:), intent (in) :: theta
    integer, intent (out) :: nspec_neo, nenergy_neo, nxi_neo, nrad_neo, ir_neo

    integer :: is, ik, ij, ig, ir, idx, ntheta, ntheta_neo
    integer :: ptheta_neo, ip
    integer :: ngrid_per_rad
    integer, save :: neo_unit, neof_unit, neophi_unit

    real, dimension (:), allocatable :: tmp, neo_coefs, dum, neo_phi, dneo_phi
    real, dimension (:), allocatable :: theta_neo_ext,neo_coefs_ext, neo_phi_ext, dneo_phi_ext !JPL
    ntheta = size(theta)

    ! for now, set ir_neo by hand, but best to derive it from neo output in future
    ir_neo = 2

    if (proc0) then
       call get_unused_unit (neo_unit)

       ! read in number of grid points from neo's grid.out file
       open (unit=neo_unit, file='neo_grid.out', status="old", action="read")
       read (neo_unit,*) nspec_neo
       read (neo_unit,*) nenergy_neo
       read (neo_unit,*) nxi_neo
       read (neo_unit,*) ntheta_neo
       if (.not. allocated(theta_neo)) allocate (theta_neo(ntheta_neo))
       do ig = 1, ntheta_neo
          read (neo_unit,*) theta_neo (ig)
       end do
       read (neo_unit,*) nrad_neo
       if (.not. allocated(rad_neo)) allocate (rad_neo(nrad_neo))
       if (.not. allocated(dens_neo)) allocate (dens_neo(nrad_neo,nspec_neo))
       if (.not. allocated(temp_neo)) allocate (temp_neo(nrad_neo,nspec_neo))
       do ir = 1, nrad_neo
          read (neo_unit,*) rad_neo(ir), dens_neo(ir,:), temp_neo(ir,:)
       end do
       close (neo_unit)
       dens_neo = dens_neo / spread(dens_neo(ir_neo,:),1,nrad_neo)
       temp_neo = temp_neo / spread(temp_neo(ir_neo,:),1,nrad_neo)
    end if

    call broadcast (nspec_neo)
    call broadcast (nenergy_neo)
    call broadcast (nxi_neo)
    call broadcast (ntheta_neo)
    call broadcast (nrad_neo)
    if (.not. allocated(theta_neo)) allocate (theta_neo(ntheta_neo))
    if (.not. allocated(rad_neo)) allocate (rad_neo(nrad_neo))
    if (.not. allocated(dens_neo)) allocate (dens_neo(nrad_neo,nspec_neo))
    if (.not. allocated(temp_neo)) allocate (temp_neo(nrad_neo,nspec_neo))
    call broadcast (theta_neo)
    call broadcast (rad_neo)
    do ir = 1, nrad_neo
       call broadcast (dens_neo)
       call broadcast (temp_neo)
    end do

    ! define number of grid points per radial location to make file i/o easier
    ngrid_per_rad = ntheta_neo*(nxi_neo+1)*nenergy_neo*nspec_neo

    if (.not. allocated(tmp)) allocate (tmp(ngrid_per_rad*nrad_neo))
    if (.not. allocated(neo_coefs)) allocate (neo_coefs(ntheta_neo), neo_phi(ntheta_neo), dneo_phi(ntheta_neo))
    if (.not. allocated(dum)) allocate (dum(ntheta))
    if (.not. allocated(coefs)) allocate (coefs(nrad_neo,ntheta,0:nxi_neo,0:nenergy_neo-1,nspec_neo))
    if (.not. allocated(phineo)) allocate (phineo(nrad_neo,ntheta))
    if (.not. allocated(dphidth)) allocate (dphidth(ntheta))
    if (.not. allocated(dcoefsdr)) then
       allocate (dcoefsdr(0:nxi_neo,0:nenergy_neo-1)) ; dcoefsdr = 0.
       allocate (dcoefsdth(ntheta,0:nxi_neo,0:nenergy_neo-1)) ; dcoefsdth = 0.
    end if

    if (proc0) then

       !JPL: if the range of theta grid in GS2 is beyond [-pi:pi](e.g. "nperiod > 1"),
       !     extend the theta grid in NEO ([-pi:pi]) to the periodic one in the theta range of GS2.
       ptheta_neo=ceiling((maxval(theta)+pi)/twopi)  !the period of 2pi in theta range
       if (ptheta_neo .gt. 1) then
          if (.not. allocated(theta_neo_ext)) allocate (theta_neo_ext(ntheta_neo*(2*ptheta_neo-1)))
          if (.not. allocated(neo_coefs_ext)) allocate (neo_coefs_ext(ntheta_neo*(2*ptheta_neo-1)))
          if (.not. allocated(neo_phi_ext)) allocate (neo_phi_ext(ntheta_neo*(2*ptheta_neo-1)))
          if (.not. allocated(dneo_phi_ext)) allocate (dneo_phi_ext(ntheta_neo*(2*ptheta_neo-1)))
       end if

       ! read in H1^{nc} (adiabatic piece of F1^{nc}) from neo's f.out file
       call get_unused_unit (neof_unit)
       open (unit=neof_unit, file='neo_f.out', status="old", action="read")

       do ir = 1, nrad_neo
          read (neof_unit,*) tmp(ngrid_per_rad*(ir-1)+1:ngrid_per_rad*ir)
       end do

       idx = 1
       do ir = 1, nrad_neo
          do is = 1, nspec_neo
             do ik = 0, nenergy_neo-1
                do ij = 0, nxi_neo
                   do ig = 1, ntheta_neo
                      neo_coefs(ig) = tmp(idx)
                      idx = idx+1
                   end do

                   if ( ptheta_neo .gt. 1) then
                      do ip = 1, ntheta_neo*(2*ptheta_neo-1)
                         theta_neo_ext(ip)=theta_neo(mod((ip-1),ntheta_neo)+1)+twopi*(int((ip-1)/ntheta_neo)-(ptheta_neo-1))
                         neo_coefs_ext(ip)=neo_coefs(mod((ip-1),ntheta_neo)+1)
                      end do
                      call lf_spline (theta_neo_ext, neo_coefs_ext, theta, coefs(ir,:,ij,ik,is), dum)  !JPL
                   else
                      ! need to interpolate coefficients from neo's theta grid to gs2's
                      call lf_spline (theta_neo, neo_coefs, theta, coefs(ir,:,ij,ik,is), dum)
                   end if
                end do
             end do
          end do
       end do

       close (neof_unit)
    end if

    deallocate (tmp)

    if (.not. allocated(tmp)) allocate (tmp(ntheta_neo*nrad_neo))

    if (proc0) then
       ! read in phi1^{nc} from neo's phi.out file
       call get_unused_unit (neophi_unit)
       open (unit=neophi_unit, file='neo_phi.out', status="old", action="read")

       do ir = 1, nrad_neo
          read (neophi_unit,*) tmp((ir-1)*ntheta_neo+1:ir*ntheta_neo)
       end do

       idx = 1
       do ir = 1, nrad_neo
          do ig = 1, ntheta_neo
             neo_phi(ig) = tmp(idx)
             idx = idx+1
          end do

          ! TMP FOR TESTING -- MAB
          !          neo_phi = 0.

          ! need to interpolate coefficients from neo's theta grid to gs2's
          if ( ptheta_neo .gt. 1) then
             do ip = 1, ntheta_neo*(2*ptheta_neo-1)
                neo_phi_ext(ip)=neo_phi(mod((ip-1),ntheta_neo)+1)
             end do
             call lf_spline (theta_neo_ext, neo_phi_ext, theta, phineo(ir,:), dum)  !JPL
          else
             call lf_spline (theta_neo, neo_phi, theta, phineo(ir,:), dum)
          end if

          ! at central radius, calculate dphi/dth and interpolate onto gs2 grid
          if (ir == 2) then
             if ( ptheta_neo .gt. 1) then
                do ip = 1, ntheta_neo*(2*ptheta_neo-1)
                   dneo_phi_ext(ip)=dneo_phi(mod((ip-1),ntheta_neo)+1)
                end do
                call get_thgrad (neo_phi_ext, theta_neo_ext, dneo_phi_ext)
                call lf_spline (theta_neo_ext, dneo_phi_ext, theta, dphidth, dum)
             else
                call get_thgrad (neo_phi, theta_neo, dneo_phi)
                call lf_spline (theta_neo, dneo_phi, theta, dphidth, dum)
             end if
          end if
       end do

       close (neophi_unit)

    end if

    call broadcast (dphidth)
    do ir = 1, nrad_neo
       call broadcast (phineo(ir,:))
    end do
    do ir = 1, nrad_neo
       do is = 1, nspec_neo
          do ik = 0, nenergy_neo-1
             do ij = 0, nxi_neo
                call broadcast (coefs(ir,:,ij,ik,is))
             end do
          end do
       end do
    end do

    deallocate (tmp, neo_coefs, theta_neo, neo_phi, dneo_phi, dum)
    if (proc0 .and. (ptheta_neo .gt. 1)) then
       deallocate (neo_coefs_ext, theta_neo_ext, neo_phi_ext, dneo_phi_ext)
    endif
  end subroutine read_neocoefs