do_write_correlation_extend Subroutine

  subroutine do_write_correlation_extend(file_id, time, time_old)
    use theta_grid, only: ntgrid
    use kt_grids, only: jtwist_out, ntheta0, naky
    use mp, only: proc0
    use gs2_io, only: nc_loop_corr_extend
    implicit none
    !> NetCDF ID of the file to write to
    integer, intent(in) :: file_id
    !> Current and previous simulation times
    real, intent(in) :: time, time_old
    complex, dimension (:,:,:), allocatable, save:: phicorr_sum
    real, dimension (:,:,:), allocatable, save :: phiextend_sum
    complex, dimension (:,:,:), allocatable :: phi_corr
    real, dimension (:,:,:), allocatable :: phi2_extend
    real, save :: tcorr0 = 0.0

    if (.not. allocated(phicorr_sum)) then
       ntg_extend = (2*ntgrid+1)*((ntheta0-1)/jtwist_out+1)
       nth0_extend = min(ntheta0,jtwist_out*(naky-1))
       ! Warning: For ntheta=ntheta0=naky=128 and jtwist = 6 these two arrays
       ! will require 750 MB and 375 MB respectively and they are not freed
       ! after leaving the routine.
       allocate (phicorr_sum(ntg_extend,ntheta0,naky)) ; phicorr_sum = 0.0
       allocate (phiextend_sum(ntg_extend,ntheta0,naky)) ; phiextend_sum = 0.0
       tcorr0 = time_old
    end if

    allocate (phi_corr(ntg_extend,ntheta0,naky))
    allocate (phi2_extend(ntg_extend,ntheta0,naky))
    call correlation_extend (phi_corr, phi2_extend)
    phicorr_sum = phicorr_sum + phi_corr*(time-time_old)
    phiextend_sum = phiextend_sum + phi2_extend*(time-time_old)
    if (proc0) then
       call nc_loop_corr_extend (file_id, nout_big, time, phicorr_sum/(time-tcorr0), phiextend_sum/(time-tcorr0))
       nout_big = nout_big + 1
    end if
    deallocate (phi_corr, phi2_extend)
  end subroutine do_write_correlation_extend