text:
! Set up the perpendicular wavenumbers by calling the appropriate sub-modules.
module kt_grids

  implicit none

  public :: init_kt_grids, finish_kt_grids
  public :: read_kt_grids_parameters, box
  public :: aky, theta0, akx, zed0
  public :: naky, nakx, nx, ny, reality
  public :: dx,dy,dkx, dky, dx_d
  public :: jtwist, jtwistfac, ikx_twist_shift, x0, y0
  public :: x, x_d
  public :: rho, rho_d, rho_clamped, rho_d_clamped
  public :: nalpha
  public :: ikx_max, naky_all
  public :: zonal_mode
  public :: swap_kxky, swap_kxky_back
  public :: swap_kxky_ordered, swap_kxky_back_ordered
  public :: multiply_by_rho, centered_in_rho
  public :: periodic_variation

  private

  interface swap_kxky
     module procedure swap_kxky_real
     module procedure swap_kxky_complex
  end interface

  real, dimension (:,:), allocatable :: theta0, zed0
  real, dimension (:), allocatable :: aky, akx
  real, dimension (:), allocatable :: x, x_d
  real, dimension (:), allocatable :: rho, rho_d, rho_clamped, rho_d_clamped
  complex, dimension (:,:), allocatable:: g0x
  real :: dx, dy, dkx, dky, dx_d
  real :: jtwistfac
  integer :: naky, nakx, nx, ny, nalpha
  integer :: jtwist, ikx_twist_shift
  integer :: ikx_max, naky_all
  logical :: reality = .false.
  logical :: centered_in_rho, periodic_variation
  character(20) :: grid_option
  logical, dimension (:), allocatable :: zonal_mode

  namelist /kt_grids_knobs/ grid_option

  ! internal variables
  integer :: gridopt_switch
  integer, parameter :: gridopt_range = 1, gridopt_box = 2

  real :: aky_min, aky_max
  real :: akx_min, akx_max
  real :: theta0_min, theta0_max
  real :: x0, y0
  logical :: read_kt_grids_initialized = .false.
  logical :: init_kt_grids_initialized = .false.
  logical :: box

contains

  subroutine read_kt_grids_parameters

    use mp, only: proc0

    implicit none

    if (read_kt_grids_initialized) return
    read_kt_grids_initialized = .true.

    if (proc0) then
       call read_grid_option
       select case (gridopt_switch)
       case (gridopt_range)
          call read_kt_grids_range
       case (gridopt_box)
          call read_kt_grids_box
       end select
    end if

    call broadcast_input

    call allocate_arrays

  end subroutine read_kt_grids_parameters

  subroutine read_grid_option

    use file_utils, only: input_unit, error_unit, input_unit_exist
    use text_options, only: text_option, get_option_value

    implicit none

    type (text_option), dimension (5), parameter :: gridopts = &
         (/ text_option('default', gridopt_range), &
         text_option('range', gridopt_range), &
         text_option('box', gridopt_box), &
         text_option('annulus', gridopt_box), &
         text_option('nonlinear', gridopt_box) /)

    integer :: ierr, in_file
    logical :: nml_exist

    grid_option = 'default'

    in_file = input_unit_exist ("kt_grids_knobs", nml_exist)
    if (nml_exist) read (unit=in_file, nml=kt_grids_knobs)

    ierr = error_unit()
    call get_option_value (grid_option, gridopts, gridopt_switch, &
         ierr, "grid_option in kt_grids_knobs")

  end subroutine read_grid_option

  subroutine read_kt_grids_box

    use file_utils, only: input_unit_exist
    use physics_flags, only: full_flux_surface

    implicit none

    integer :: in_file
    logical :: exist

    namelist /kt_grids_box_parameters/ nx, ny, jtwist, jtwistfac, y0, &
                                       centered_in_rho, periodic_variation

    ! note that jtwist and y0 will possibly be modified
    ! later in init_kt_grids_box if they make it out
    ! of this subroutine with negative values
    ! it is necessary to wait until later to do this check
    ! because the values to which they may be set will
    ! depend on information from the geometry module,
    ! which itself may rely on ny from here (number of alphas)
    nx = 1
    ny = 1
    jtwist = -1
    jtwistfac = 1.0
    y0 = -1.0
    nalpha = 1
    centered_in_rho = .true.
    periodic_variation = .false.

    in_file = input_unit_exist("kt_grids_box_parameters", exist)
    if (exist) read (in_file, nml=kt_grids_box_parameters)

    ! get the number of de-aliased modes in y and x
    naky = (ny-1)/3 + 1
    nakx = 2*((nx-1)/3) +  1

    reality = .true.

    if (full_flux_surface) nalpha = ny

  end subroutine read_kt_grids_box

  subroutine read_kt_grids_range

    use file_utils, only: input_unit, input_unit_exist

    implicit none

    integer :: in_file
    logical :: exist

    namelist /kt_grids_range_parameters/ naky, nakx,  &
         aky_min, aky_max, theta0_min, theta0_max, akx_min, akx_max

    nalpha = 1
    naky = 1
    nakx = 1
    aky_min = 0.0
    aky_max = 0.0
    ! set these to be nonsense values
    ! so can check later if they've been set
    akx_min = 0.0
    akx_max = -1.0
    theta0_min = 0.0
    theta0_max = -1.0

    in_file = input_unit_exist ("kt_grids_range_parameters", exist)
    if (exist) read (in_file, nml=kt_grids_range_parameters)

  end subroutine read_kt_grids_range

  subroutine init_kt_grids

    use common_types, only: flux_surface_type
    use zgrid, only: init_zgrid

    implicit none

    if (init_kt_grids_initialized) return
    init_kt_grids_initialized = .true.

    call init_zgrid

    select case (gridopt_switch)
    case (gridopt_range)
       call init_kt_grids_range
    case (gridopt_box)
       call init_kt_grids_box
    end select

    ! determine if iky corresponds to zonal mode
    if (.not.allocated(zonal_mode)) allocate (zonal_mode(naky))
    zonal_mode = .false.
    if (abs(aky(1)) < epsilon(0.)) zonal_mode(1) = .true.

  end subroutine init_kt_grids

  subroutine init_kt_grids_box

    use mp, only: mp_abort
    use common_types, only: flux_surface_type
    use constants, only: pi
    use stella_geometry, only: geo_surf, twist_and_shift_geo_fac
    use stella_geometry, only: q_as_x, get_x_to_rho, dxdXcoord, drhodpsi
    use physics_parameters, only: rhostar
    use physics_flags, only: full_flux_surface, radial_variation
    use file_utils, only: runtype_option_switch, runtype_multibox
    use zgrid, only: shat_zero

    implicit none

    integer :: ikx, iky
    real :: x_shift, dqdrho, pfac

    box = .true.

    ! set jtwist and y0 for cases where they have not been specified
    ! and for which it makes sense to set them automatically
    if (jtwist < 1) then
      jtwist = max(1,int(abs(twist_and_shift_geo_fac)+0.5))
      jtwist = max(1,int(jtwistfac*jtwist+0.5))
    endif
    ! signed version of jtwist, with sign determined by, e.g., magnetic shear
    ikx_twist_shift = -jtwist*int(sign(1.0,twist_and_shift_geo_fac))

    if (y0 < 0.) then
       if (full_flux_surface) then
          ! if simulating a flux annulus, then
          ! y0 determined by the physical
          ! extent of the device
          if (rhostar > 0.) then
             y0 = 1./(rhostar*geo_surf%rhotor)
          else
             call mp_abort ('must set rhostar if simulating a full flux surface. aborting.')
          end if
       else
          ! if simulating a flux tube
          ! makes no sense to have y0 < 0.0
          ! so abort
          call mp_abort ('y0 negative only makes sense when simulating a flux annulus.  aborting.')
       end if
    end if

    ! get the grid spacing in ky and then in kx using twist-and-shift BC
    dky = 1./y0
    ! non-quantized b/c assumed to be periodic instead
    ! of linked boundary conditions if zero magnetic shear
    if (abs(geo_surf%shat) <= shat_zero) then
       dkx = dky / real(jtwist)
    else
       dkx = dky * abs(twist_and_shift_geo_fac) / real(jtwist)
    end if

    x0 = 1./dkx

    ! ky goes from zero to ky_max
    do iky = 1, naky
       aky(iky) = real(iky-1)*dky
    end do

    ! get the ikx index corresponding to kx_max
    ikx_max = nakx/2+1

    ! get the total number of ky values, including negative ky
    naky_all = 2*naky-1

    ! kx goes from zero to kx_max down to zero...
    do ikx = 1, ikx_max
       akx(ikx) = real(ikx-1)*dkx
    end do
    ! and then from -kx_max to -|kx_min|
    do ikx = ikx_max+1, nakx
       akx(ikx) = real(ikx-nakx-1)*dkx
    end do

    ! set theta0=0 for ky=0
    theta0(1,:) = 0.0
    if (q_as_x) then
       do ikx = 1, nakx
          ! theta0 = kx/ky
          theta0(2:,ikx) = akx(ikx)/aky(2:)
       end do
    else if (abs(geo_surf%shat) > shat_zero) then
       do ikx = 1, nakx
          ! theta0 = kx/ky/shat
          theta0(2:,ikx) = akx(ikx)/(aky(2:)*geo_surf%shat)
       end do
    else
       do ikx = 1, nakx
          ! if shat=0, theta0 is meaningless, so be careful
          theta0(2:,ikx) = - akx(ikx)/aky(2:)
       end do
    end if

    ! for radial variation
    if(.not.allocated(x)) allocate (x(nx))
    if(.not.allocated(x_d)) allocate (x_d(nakx))
    if(.not.allocated(rho)) allocate (rho(nx))
    if(.not.allocated(rho_d)) allocate (rho_d(nakx))
    if(.not.allocated(rho_clamped)) allocate (rho_clamped(nx))
    if(.not.allocated(rho_d_clamped)) allocate (rho_d_clamped(nakx))

    dx = (2*pi*x0)/nx
    dy = (2*pi*y0)/ny

    x_shift = pi*x0
    pfac = 1.0
    if (periodic_variation) pfac = 0.5
    if(centered_in_rho) then
      if(q_as_x) then
        dqdrho = geo_surf%shat*geo_surf%qinp/geo_surf%rhoc
        x_shift = pi*x0*(1.0 &
                - 0.5*pfac*rhostar*pi*x0*geo_surf%d2qdr2/(dqdrho**2*dxdXcoord))
      else
        x_shift = pi*x0*(1.0 &
                - 0.5*pfac*rhostar*pi*x0*geo_surf%d2psidr2*drhodpsi**2/dxdXcoord)
      endif
    endif

    do ikx = 1, nx
      if (radial_variation.or.runtype_option_switch.eq.runtype_multibox) then
        if(periodic_variation) then
          if(ikx.le.nx/2) then
            x(ikx) = (ikx-1)*dx - 0.5*x_shift
          else
            x(ikx) = x(nx-ikx+1)
          endif
        else
          x(ikx) = (ikx-0.5)*dx - x_shift
        endif
      else
        x(ikx) = (ikx-1)*dx
      endif
    enddo

    dx_d = (2*pi*x0)/nakx
    do ikx = 1, nakx
      if (radial_variation.or.runtype_option_switch.eq.runtype_multibox) then
        if(periodic_variation) then
          if(ikx.le.(nakx+1)/2) then
            x_d(ikx) = (ikx-1)*dx_d - 0.5*x_shift
          else
            x_d(ikx) = x_d(nakx-ikx+1)
          endif
        else
          x_d(ikx) = (ikx-0.5)*dx_d - x_shift
        endif
      else
        x_d(ikx) = (ikx-1)*dx_d
      endif
    enddo

    call get_x_to_rho(1,x,rho)
    call get_x_to_rho(1,x_d,rho_d)

    if(.not.allocated(rho_clamped)) allocate(rho_clamped(nx)); rho_clamped = rho
    if(.not.allocated(rho_d_clamped)) allocate(rho_d_clamped(nakx)); rho_d_clamped = rho_d

    zed0 = theta0*geo_surf%zed0_fac

    if(radial_variation) call dump_radial_grid

    if(radial_variation.and.(any((rho+geo_surf%rhoc).lt.0.0) &
                             .or.any((rho+geo_surf%rhoc).gt.1.0))) then
      call mp_abort ('rho(x) is beyond range [0,1]. Try changing rhostar or q/psi profiles')
    endif

  end subroutine init_kt_grids_box

  subroutine init_kt_grids_range

    use mp, only: mp_abort
    use common_types, only: flux_surface_type
    use stella_geometry, only: geo_surf, q_as_x
    use zgrid, only: shat_zero

    implicit none

    integer :: i, j
    real :: dkx, dky, dtheta0, tfac
    real :: zero

    box = .false.

    ! NB: we are assuming here that all ky are positive
    ! when running in range mode
    dky = 0.0
    if (naky > 1) dky = (aky_max - aky_min)/real(naky - 1)
    aky = (/ (aky_min + dky*real(i), i = 0,naky-1) /)

    ! set default akx and theta0 to 0
    akx = 0.0 ; theta0 = 0.0

    if (q_as_x) then
      tfac = 1.0
    else
      tfac = geo_surf%shat
    endif

    zero = 100.*epsilon(0.)

    ! if theta0_min and theta0_max have been specified,
    ! use them to determine akx_min and akx_max
    if (theta0_max > theta0_min-zero) then
       if (geo_surf%shat > epsilon(0.)) then
          akx_min = theta0_min * tfac * aky(1)
          akx_max = theta0_max * tfac * aky(1)
       else
          akx_min = theta0_max * tfac * aky(1)
          akx_max = theta0_min * tfac * aky(1)
       end if
    end if

    ! shat_zero is minimum shat value below which periodic BC is enforced
    if (abs(geo_surf%shat) > shat_zero) then  ! ie assumes boundary_option .eq. 'linked'
       ! if akx_min and akx_max specified in input
       ! instead of theta0_min and theta0_max,
       ! use them to get theta0_min and theta0_max
       if (theta0_min > theta0_max+zero .and. abs(aky(1)) > zero) then
          theta0_min = akx_min/(tfac*aky(1))
          theta0_max = akx_max/(tfac*aky(1))
          dtheta0 = 0.0
          if (nakx > 1) dtheta0 = (theta0_max - theta0_min)/real(nakx - 1)

          do j = 1, naky
             theta0(j,:) &
                  = (/ (theta0_min + dtheta0*real(i), i=0,nakx-1) /)
          end do
          akx = theta0(1,:) * tfac * aky(1)
       else if (akx_max > akx_min-zero) then
          dkx = 0.0
          if (nakx > 1) dkx = (akx_max - akx_min)/real(nakx - 1)
          akx = (/ (akx_min + dkx*real(i), i = 0,nakx-1) /)

          dtheta0 = 0.0
          if (nakx > 1) dtheta0 = (theta0_max - theta0_min)/real(nakx - 1)

          if (geo_surf%shat > epsilon(0.)) then
             do j = 1, naky
                theta0(j,:) &
                     = (/ (theta0_min + dtheta0*real(i), i=0,nakx-1) /)
             end do
          else
             do j = 1, naky
                theta0(j,:) &
                     = (/ (theta0_min + dtheta0*real(i), i=nakx-1,0,-1) /)
             end do
          end if
       else
          call mp_abort ('ky=0 is inconsistent with kx_min different from kx_max. aborting.')
       end if

    else
       ! here assume boundary_option .eq. 'periodic'
       ! used for periodic finite kx ballooning space runs with shat=0
       dkx = 0.0
       if (nakx > 1) dkx = (akx_max - akx_min)/real(nakx - 1)
       akx = (/ (akx_min + dkx*real(i), i = 0,nakx-1) /)
    endif

    ikx_max = nakx
    naky_all = naky

  end subroutine init_kt_grids_range

  subroutine broadcast_input

    use mp, only: broadcast

    implicit none

    call broadcast (gridopt_switch)
    call broadcast (centered_in_rho)
    call broadcast (periodic_variation)
    call broadcast (naky)
    call broadcast (nakx)
    call broadcast (ny)
    call broadcast (nx)
    call broadcast (nalpha)
    call broadcast (reality)
    call broadcast (jtwist)
    call broadcast (jtwistfac)
    call broadcast (y0)
    call broadcast (aky_min)
    call broadcast (aky_max)
    call broadcast (akx_min)
    call broadcast (akx_max)
    call broadcast (theta0_min)
    call broadcast (theta0_max)

  end subroutine broadcast_input

  subroutine dump_radial_grid

    use file_utils, only: run_name
    use physics_parameters, only: rhostar
    use stella_geometry, only: q_as_x, geo_surf, dxdXcoord, drhodpsi

    implicit none

    integer :: ix
    character (300) :: filename

    filename = trim(trim(run_name)//'.radial_grid')
    open (1047,file=filename,status='unknown')
    if(q_as_x) then
      write (1047,'(1a12,1e12.4,1a12,1e12.4,1a12,1e12.4,1a12,1e12.4)') &
            '#dxdXcoord = ', dxdXcoord, &
            ' q    = '  , geo_surf%qinp, &
            ' dqdr = '  , geo_surf%shat*geo_surf%qinp/geo_surf%rhoc, &
            ' d2qdr2 = ', geo_surf%d2qdr2
      write (1047,'(3a12)') '#1.x', '2.q', '3.rho'
      do ix = 1, nx
        write (1047,'(3e12.4,i9)') &
          x(ix), &
          rhostar*x(ix)/dxdXcoord + geo_surf%qinp, &
          rho(ix) + geo_surf%rhoc
      enddo
    else
      write (1047,'(1a12,1e12.4,1a12,1e12.4,1a12,1e12.4,1a12,1e12.4)') &
            '#dxdXcoord = ', dxdXcoord, &
            ' dpsidr    = ' , 1.0/drhodpsi, &
            ' d2psidr2 = '  , geo_surf%d2psidr2
      write (1047,'(3a12,a9)') '#1.x', '2.psi-psi0', '3.rho'
      do ix = 1, nx
        write (1047,'(3e12.4,i9)') &
            x(ix), &
            rhostar*x(ix)/dxdXcoord, &
            rho(ix) + geo_surf%rhoc
      enddo
    endif

    close (1047)

  end subroutine dump_radial_grid

  subroutine allocate_arrays

    implicit none

    allocate (akx(nakx))
    allocate (aky(naky))
    allocate (theta0(naky,nakx))
    allocate (zed0(naky,nakx))

  end subroutine allocate_arrays

  ! take an array with ky >= 0 and all kx (ordered like 0, ..., kxmax, -kxmax, ..., -dkx)
  ! and uses reality condition to return array
  ! with kx >= 0 and all ky (ordered like 0, ..., kymax, -kymax, ..., -dky)
  subroutine swap_kxky_complex (gin, gout)

    implicit none

    complex, dimension (:,:), intent (in) :: gin
    complex, dimension (:,:), intent (out) :: gout

    integer :: ikx, ikxneg
    integer :: iky, ikyneg

    ! first set arrays equal for ky >= 0 and kx >= 0
    gout(:naky,:) = gin(:,:ikx_max)
    ! next fill in ky < 0, kx >= 0 elements of array using reality
    ikx = 1
    ikxneg = ikx
    do iky = naky+1, naky_all
       ! this is the ky index corresponding to +ky in original array
       ikyneg = naky_all-iky+2
       gout(iky,ikx) = conjg(gin(ikyneg,ikxneg))
    end do
    do ikx = 2, ikx_max
       ikxneg = nakx-ikx+2
       do iky = naky+1, naky_all
          ! this is the ky index corresponding to +ky in original array
          ikyneg = naky_all-iky+2
          gout(iky,ikx) = conjg(gin(ikyneg,ikxneg))
       end do
    end do

  end subroutine swap_kxky_complex

  ! take an array with ky >= 0 and all kx (ordered like 0, ..., kxmax, -kxmax, ..., -dkx)
  ! and uses reality condition to return array
  ! with kx >= 0 and all ky (ordered like 0, ..., kymax, -kymax, ..., -dky)
  subroutine swap_kxky_real (gin, gout)

    implicit none

    real, dimension (:,:), intent (in) :: gin
    real, dimension (:,:), intent (out) :: gout

    integer :: ikx, ikxneg
    integer :: iky, ikyneg

    ! first set arrays equal for ky >= 0 and kx >= 0
    gout(:naky,:) = gin(:,:ikx_max)
    ! next fill in ky < 0, kx >= 0 elements of array using reality
    ikx = 1
    ikxneg = ikx
    do iky = naky+1, naky_all
       ! this is the ky index corresponding to +ky in original array
       ikyneg = naky_all-iky+2
       gout(iky,ikx) = gin(ikyneg,ikxneg)
    end do
    do ikx = 2, ikx_max
       ikxneg = nakx-ikx+2
       do iky = naky+1, naky_all
          ! this is the ky index corresponding to +ky in original array
          ikyneg = naky_all-iky+2
          gout(iky,ikx) = gin(ikyneg,ikxneg)
       end do
    end do

  end subroutine swap_kxky_real

  ! take an array with ky >= 0 and all kx (ordered like 0, ..., kxmax, -kxmax, ..., -dkx)
  ! and uses reality condition to return array
  ! with kx >= 0 and all ky (ordered like -kymax, ..., 0, ..., kymax)
  subroutine swap_kxky_ordered (gin, gout)

    implicit none

    complex, dimension (:,:), intent (in) :: gin
    complex, dimension (:,:), intent (out) :: gout

    integer :: ikx, ikxneg
    integer :: iky, ikyneg

    ! first set arrays equal for ky >= 0 and kx >= 0
    gout(naky:,:) = gin(:,:ikx_max)
    ! next fill in ky < 0, kx >= 0 elements of array using reality
    ikx = 1
    ikxneg = ikx
    do iky = 1, naky-1
       ! this is the ky index corresponding to +ky in original array
       ikyneg = naky-iky+1
       gout(iky,ikx) = conjg(gin(ikyneg,ikxneg))
    end do
    do ikx = 2, ikx_max
       ikxneg = nakx-ikx+2
       do iky = 1, naky-1
          ! this is the ky index corresponding to +ky in original array
          ikyneg = naky-iky+1
          gout(iky,ikx) = conjg(gin(ikyneg,ikxneg))
       end do
    end do

  end subroutine swap_kxky_ordered

  ! take an array with kx >= 0 and all ky (ordered like 0, ..., kymax, -kymax, ..., -dky)
  ! and uses reality condition to return array
  ! with ky >= 0 and all kx (ordered like 0, ..., kxmax, -kxmax, ..., -dkx)
  subroutine swap_kxky_back (gin, gout)

    implicit none

    complex, dimension (:,:), intent (in) :: gin
    complex, dimension (:,:), intent (out) :: gout

    integer :: ikx, ikxneg
    integer :: iky, ikyneg

    ! first set arrays equal for ky >= 0 and kx >= 0
    gout(:,:ikx_max) = gin(:naky,:)
    ! next fill in kx < 0, ky >= 0 elements of array using reality
    do ikx = ikx_max+1, nakx
       ikxneg = nakx-ikx+2
       iky = 1
       ikyneg = iky
       gout(iky,ikx) = conjg(gin(ikyneg,ikxneg))
       do iky = 2, naky
          ikyneg = naky_all-iky+2
          gout(iky,ikx) = conjg(gin(ikyneg,ikxneg))
       end do
    end do

  end subroutine swap_kxky_back

  ! take an array with kx >= 0 and all ky (ordered like -kymax, ..., 0, ..., kymax)
  ! and uses reality condition to return array
  ! with ky >= 0 and all kx (ordered like 0, ..., kxmax, -kxmax, ..., -dkx)
  subroutine swap_kxky_back_ordered (gin, gout)

    implicit none

    complex, dimension (:,:), intent (in) :: gin
    complex, dimension (:,:), intent (out) :: gout

    integer :: ikx, ikxneg
    integer :: iky, ikyneg

    ! first set arrays equal for ky >= 0 and kx >= 0
    gout(:,:ikx_max) = gin(naky:,:)
    ! next fill in kx < 0, ky >= 0 elements of array using reality
    do ikx = ikx_max+1, nakx
       ikxneg = nakx-ikx+2
       do iky = 1, naky
          ikyneg = naky-iky+1
          gout(iky,ikx) = conjg(gin(ikyneg,ikxneg))
       end do
    end do

  end subroutine swap_kxky_back_ordered

  subroutine finish_kt_grids

    implicit none

    if (allocated(aky)) deallocate (aky)
    if (allocated(akx)) deallocate (akx)
    if (allocated(theta0)) deallocate (theta0)
    if (allocated(zed0)) deallocate (zed0)

    if (allocated(x)) deallocate (x)
    if (allocated(x_d)) deallocate (x_d)
    if (allocated(rho)) deallocate (rho)
    if (allocated(rho_d)) deallocate (rho_d)
    if (allocated(rho_clamped)) deallocate (rho_clamped)
    if (allocated(rho_d_clamped)) deallocate (rho_d_clamped)

    if (allocated(g0x)) deallocate (g0x)

    reality = .false.
    read_kt_grids_initialized = .false.
    init_kt_grids_initialized = .false.

  end subroutine finish_kt_grids

  subroutine multiply_by_rho (gin)

    use stella_transforms, only: transform_kx2x_unpadded, transform_x2kx_unpadded
!   use stella_transforms, only: transform_kx2x_xfirst, transform_x2kx_xfirst

    implicit none

    complex, dimension (:,:), intent (inout) :: gin

    if(.not.allocated(g0x)) allocate(g0x(naky,nakx))

    call transform_kx2x_unpadded(gin,g0x)
    g0x = spread(rho_d_clamped,1,naky)*g0x
    if(zonal_mode(1)) g0x(1,:) = real(g0x(1,:))
    call transform_x2kx_unpadded(g0x,gin)

!   if(.not.allocated(g0x)) allocate(g0x(naky,nx))

!   call transform_kx2x_xfirst(gin,g0x)
!   g0x = spread(rho_clamped,1,naky)*g0x
!   if(zonal_mode(1)) g0x(1,:) = real(g0x(1,:))
!   call transform_x2kx_xfirst(g0x,gin)

  end subroutine multiply_by_rho

end module kt_grids

max_stars_repo_name: ZedThree/stella

max_issues_repo_name: ZedThree/stella

max_forks_repo_name: ZedThree/stella


text:
c
c
c
c
c
        subroutine testit(n)
        implicit real*8 (a-h,o-z)

        call test_integration()
ccc        call utest()
ccc        call utest2()

        return
        end
c
c
c
c
c
        subroutine test_integration()
        implicit real *8 (a-h,o-z)
        real*8 a(20 000 000), y(500 000)
        real *8, allocatable :: dsums(:), dexps(:),
     1     dsums2(:), dds(:), stds(:), stds2(:)
        character(100) csv_file, filename

c
c        parameters
c
        csv_file = 'test_dense_params.dat'
        csv_file = 'abort_params.dat'
        csv_file = 'params.dat'
        call read_params(csv_file, n, k1, k2, a, y)
ccc        call prin2('a = *',a,n*k)
ccc        call prin2('y = *',y,n)
        k = k1+k2
        call prinf('n*', n, 1)
        call prinf('k1*', k1, 1)
        call prinf('k2*', k2, 1)

        k3 = k+3
        allocate(dsums(k3))
        allocate(dsums2(k3))
        allocate(dexps(k3))
        allocate(dds(k3))
        allocate(stds(k3))
        allocate(stds2(k3))

c
c        integrate
c
        nn = 60
        nn_theta = 20
        nn_theta = 4
        call dense_eval(nn_theta, nn, n, k1, k2, a, y, dsums, dsum,
     1     stds)
ccc        call prin2('dsums*', dsums, k+3)
ccc        call prin2('stds*', stds, k+3)
c
c        double nodes in all directions
c
        nn2 = 2*nn
        nn_theta2 = 2*nn_theta
        call dense_eval(nn_theta2, nn2, n, k1, k2, a, y, dsums2, dsum2,
     1     stds2)

c
c        check error
c
        call dd_abs_max(dsums2, dsums, k+3, dd_max)
        call prin2('max posterior mean error*', dd_max, 1)

        call dd_abs_max(stds2, stds, k+3, dd_max)
        call prin2('max posterior std error*', dd_max, 1)

        filename = 'exps.dat'
        call write_exps_stds(filename, k, dsums, stds)
ccc        call prin2('dsums*', dsums,k+3)

c
c        stan comparison
c
        csv_file = 'test_means_dense_stan2.dat'
        csv_file = 'test_means_dense_stan.dat'
        csv_file = 'means.dat'
        call read_means(csv_file, k, dexps)
ccc        call prin2('stan expectations*', dexps, k+3)

        call dd_abs_max(dexps, dsums, k, dd_max)
        call prin2('stan max error*', dd_max, 1)

        return
        end
c
c
c
c
c
        subroutine utest2()
        implicit real *8 (a-h,o-z)
        real*8 a(10 000 000), x(100 000), y(1000 000), dexps(10 000),
     1     dsums2(10 000), dsums(100 000), dds(10 000), stds(10 000),
     1     stds2(10 000)
        character(100) csv_file, filename

c
c        compare the results of the algorithm of this file with
c        the hardcoded results from a high accuracy run from this
c        same algorithm
c

c
c        read in unit test parameters
c
        csv_file = 'test_dense_params2.dat'
        call read_params(csv_file, n, k1, k2, a, y)
        k = k1+k2
        call prinf('n*', n, 1)
        call prinf('k1*', k1, 1)
        call prinf('k2*', k2, 1)

c
c        integrate
c
        nn = 80
        nn_theta = 60
        call dense_eval(nn_theta, nn, n, k1, k2, a, y, dsums, dsum,
     1     stds)

c
c        compare the result computed here with hardcoded results
c
        filename = 'test_means_dense_fast2.dat'
        call read_exps_stds(filename, k, dsums2, stds2)

        call dd_abs_max(dsums2, dsums, k+3, dd_max)
        call prin2('max posterior mean error*', dd_max, 1)

        call dd_abs_max(stds2, stds, k+3, dd_max)
        call prin2('max posterior std error*', dd_max, 1)

c
c        stan comparison
c
        filename = 'test_means_dense_stan2.dat'
        call read_means(filename, k, dexps)
ccc        call prin2('stan expectations*', dexps, k+3)

        call dd_abs_max(dexps, dsums, k, dd_max)
        call prin2('stan max error*', dd_max, 1)

        return
        end
c
c
c
c
c
        subroutine utest()
        implicit real *8 (a-h,o-z)
        real*8 a(10 000 000), x(100 000), y(1000 000), dexps(10 000),
     1     dsums2(10 000), dsums(100 000), dds(10 000), stds(10 000),
     1     stds2(10 000)
        character(100) csv_file, filename

c
c        compare the results of the algorithm of this file with
c        the hardcoded results from a high accuracy run from this
c        same algorithm
c

c
c        read in parameters from unit test
c
        csv_file = 'test_dense_params.dat'
        call read_params(csv_file, n, k1, k2, a, y)
        k = k1+k2
        call prinf('n*', n, 1)
        call prinf('k1*', k1, 1)
        call prinf('k2*', k2, 1)

c
c        integrate
c
        nn = 80
        nn_theta = 60
        call dense_eval(nn_theta, nn, n, k1, k2, a, y, dsums, dsum,
     1     stds)

c
c        compare to hardcoded results
c
        filename = 'test_means_dense_fast.dat'
        call read_exps_stds(filename, k, dsums2, stds2)

        call dd_abs_max(dsums2, dsums, k+3, dd_max)
        call prin2('max posterior mean error*', dd_max, 1)

        call dd_abs_max(stds2, stds, k+3, dd_max)
        call prin2('max posterior std error*', dd_max, 1)

c
c        and compare to stan
c
        filename = 'test_means_dense_stan.dat'
        call read_means(filename, k, dexps)
ccc        call prin2('stan expectations*', dexps, k+3)

        call dd_abs_max(dexps, dsums, k, dd_max)
        call prin2('stan max error*', dd_max, 1)

        return
        end
c
c
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c       this is the end of the debugging code and the beginning of the
c       least squares subroutines proper
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c
c         This file contains 8 user-callable subroutines: nrleastsq,
c         nrleastsq2, nrleamatlr, rleamatr, nrleamatrr, nrleamatll,
c         nrleamatrl, nrleas_proj. Following is a brief description
c         of these subroutines.
c
c   dense_eval - evaluates two-group normal normal model with
c         normal+(0, 1) hyperpriors on the scale parameters
c
c
c
c
c
c
        subroutine dense_eval(nn_theta, nn, n, k1, k2, a, y, dsums,
     1     dsum, stds)
        implicit real *8 (a-h,o-z)
        real*8 a(*), y(*), dsums(*), stds(*),
     1     thetas(nn_theta),rhos(nn),phis(nn), s(k1+k2),s2(k1+k2),
     1     xs(k1+k2),
     1     ysmall(k1+k2),ys2(k1+k2),vmoms(100 000), dsumsi(k1+k2+3+10),
     1     whts_thetas(nn_theta+10),whts_phis(nn+10),whts_rhos(nn+10),
     1     dsum_xsi(k1+k2+10), dsum_vars(k1+k2+10)
        real *8, allocatable :: u(:,:),vt(:,:),asca(:,:),ut(:,:),
     1     b(:,:),ynew(:),v(:,:),bb(:,:), fmaxs(:), fs(:,:),
     1     dsums_cov(:,:), dsums_covi(:,:), wwti(:,:), xxti(:,:),
     1     gs(:,:), c(:,:), ct(:,:), tmp_mat(:,:)

c
c        this subroutine computes posterior expectations of the two
c        group bayesian linear regression model
c
c          y ~ normal(A1*x1 + A2*x2, sig1)
c          x1 ~ normal(0, sig2)
c          x2 ~ normal(0, sig3)
c          sig{1,2,3} ~ normal+(0, 1)
c
c        this subroutine computes E[x_i], var(x_i), E[sig{1,2,3}].
c        after a change of variables, from sig1, sig2, sig3 to
c        spherical coordinates, we integrate over a 3-dimensional
c        grid using a tensor product of gaussian nodes.
c        for each value of the outer integral
c        we do a diagonalization where the conditional density is
c        a gaussian with a diagonal covariance.
c

        done = 1.0d0
        pi = 4*atan(done)

        k = k1+k2

        allocate (u(n,k))
        allocate (ut(k,n))
        allocate (vt(k,k))
        allocate (v(k,k))
        allocate (asca(n,k))
        allocate (b(k,k))
        allocate (bb(k,k))
        allocate (ynew(k))
        allocate (fmaxs(nn))
        allocate (fs(nn_theta,nn))
        allocate (gs(nn_theta,nn))
        allocate (dsums_cov(k,k))
        allocate (dsums_covi(k,k))
        allocate (wwti(k,k))
        allocate (xxti(k,k))
        allocate (c(k,k))
        allocate (ct(k,k))
        allocate (tmp_mat(k,k))

c
c        before computing integral find least squares solution to
c        ax = y as well as projection of y onto the
c        columns of a
c
        call dsvd(n,k, a, u, s, vt)
        call mat_trans(u, n, k, ut)
        call mat_mult(ut, k, n, a, k, b)
        call y_proj(y,u,s,n,k,ysmall,resid)
        call mat_vec_mult(ut, k, n, y, ynew)

c
c        construct grid for outer parameter of integration, theta
c
        theta0 = 0.0
        theta1 = pi/2.0d0
        call theta_lege_nodes_whts(nn_theta, theta0, theta1, thetas,
     1     whts_thetas)

c
c        initialize sums (integrals) to be computed
c          dsums - expectations
c          dsums_cov - covariance
c          ss1 - E[sig1**2]
c          ss2 - E[sig2**2]
c          ss3 - E[sig3**2]
c          fm - scaling constant
c
        ktmp = k+3
        call set_zero(ktmp, dsums)
        dsum = 0
        call set_zero(k*k, dsums_cov)
        ss1 = 0
        ss2 = 0
        ss3 = 0
        fm = -1.0d250

c
c        theta integral
c
        do 130 i=1,nn_theta
        t = thetas(i)
        wht_theta = whts_thetas(i)
c
c        compute the diagonal form of a^t*a for each theta which
c        will be used to convert the integrand to a diagonal
c        gaussian for each (theta, phi, rho)
c
        call rescale_a(t,b,asca,k,k,k1,k2)
        call dsvd(k,k, asca, u, s, vt)
        call y_proj(ynew,u,s,k,k,ysmall,res2)
        call entry_sqr(ysmall,ys2,k)
        call entry_sqr(s,s2,k)

c
c        initialize sums taken in inner loop
c          wwti is E[w*w^t] where w is in the coordinate
c          system that depends on theta
c
        dsumi = 0
        ss1i = 0
        ss2i = 0
        ss3i = 0
        ktmp = k+3
        call set_zero(ktmp, dsumsi)
        call set_zero(k**2, dsums_covi)
        call set_zero(k**2, wwti)
        call set_zero(k, dsum_vars)
        fmi = -1.0d250

c
c        for each theta, get upper integration bound for phi
c
        call get_phi1_fmax(nn, n, k, t, ysmall, ys2, s, s2, resid,
     1     phi1i, fmax_theta)
        fmaxs(i) = fmax_theta
ccc        call prin2('phi1i*', phi1i, 1)
ccc        call prin2('fmax*', fmax, 1)

c
c        ... and lay down nodes
c
        phi0 = 0.0
        call lege_nodes_whts(nn, phi0, phi1i, phis, whts_phis)
ccc        call prin2('phis new*', phis, nn)

c
c          phi integral
c
        do 120 j=1,nn
        phi = phis(j)
        call compute_beta_alpha_prefact(phi,ysmall,ys2,s,s2,
     1      n,k,resid,alpha,beta,prefact,exp_fact)
        call compute_mjs(k,s2,ysmall,phi,vmoms)

        wht_phi = whts_phis(j)

c
c          rho integral
c
        call eval_rho_int(nn, n, k, exp_fact,
     1     dsum_rho, dsum_rho1, dsum_rho2, fmax)
c
c        for fixed theta, compute integral over phi by a sum of
c        the form \sum_i exp(fi)*gi so that
c        we have an expression exp(fmi)*dsum, the reason we do
c        this is that exp(fi) is often smaller than 10^{-250}
c
        fi = prefact + fmax - fmaxs(1)
        gi = dsum_rho
        gi1 = dsum_rho1
        gi2 = dsum_rho2
        wt = wht_phi
        fs(i, j) = fi
        gs(i, j) = gi
        if (fi .gt. fmi) then
          dsumi = dsumi * exp(fmi-fi) + gi*wt

          do ijk=1,k
          dsumsi(ijk) = dsumsi(ijk) * exp(fmi - fi) + gi*wt*vmoms(ijk)
          enddo

          dsumsi(k+1) = dsumsi(k+1)*exp(fmi-fi) + gi1*cos(phi)*wt
          dsumsi(k+2) = dsumsi(k+2)*exp(fmi-fi) + gi1*sin(phi)*cos(t)*wt
          dsumsi(k+3) = dsumsi(k+3)*exp(fmi-fi) + gi1*sin(phi)*sin(t)*wt

          ss1i = ss1i*exp(fmi-fi) + gi2*cos(phi)**2*wt
          ss2i = ss2i*exp(fmi-fi) + gi2*sin(phi)**2*cos(t)**2*wt
          ss3i = ss3i*exp(fmi-fi) + gi2*sin(phi)**2*sin(t)**2*wt

          do ijk=1,k
          coef = 1/(s2(ijk)/cos(phi)**2 + 1/sin(phi)**2)
          dsum_vars(ijk) = dsum_vars(ijk) * exp(fmi-fi) + coef*gi2*wt
          enddo

          do i1=1,k
          do i2=1,k
          wwti(i1, i2)=wwti(i1,i2)*exp(fmi-fi)+gi*vmoms(i1)*vmoms(i2)*wt
          enddo
          enddo

         fmi = fi
        else
          dsumi = dsumi + exp(fi-fmi) * gi*wt

          do ijk=1,k
          dsumsi(ijk) = dsumsi(ijk) + exp(fi - fmi) * gi*wt*vmoms(ijk)
          enddo

          dsumsi(k+1) = dsumsi(k+1) + exp(fi-fmi)*gi1*cos(phi)*wt
          dsumsi(k+2) = dsumsi(k+2) + exp(fi-fmi)*gi1*sin(phi)*cos(t)*wt
          dsumsi(k+3) = dsumsi(k+3) + exp(fi-fmi)*gi1*sin(phi)*sin(t)*wt

          ss1i = ss1i + exp(fi-fmi)*gi2*cos(phi)**2*wt
          ss2i = ss2i + exp(fi-fmi)*gi2*sin(phi)**2*cos(t)**2*wt
          ss3i = ss3i + exp(fi-fmi)*gi2*sin(phi)**2*sin(t)**2*wt

          do ijk=1,k
          coef = 1/(s2(ijk)/cos(phi)**2 + 1/sin(phi)**2)
          dsum_vars(ijk) = dsum_vars(ijk) + exp(fi-fmi)*coef*gi2*wt
          enddo

          do i1=1,k
          do i2=1,k
          wwti(i1, i2)=wwti(i1,i2)+exp(fi-fmi)*gi*vmoms(i1)*vmoms(i2)*wt
          enddo
          enddo

        endif

 120    continue

c
c        for covariance we need to convert back to original coordinate
c        system. we have a different change of variables for each value
c        of theta, so do this for each theta
c
        call mat_trans(vt, k, k, v)
        call get_xs_to_ws_matrix(k,k1,k2,v,t,c,ct)
        call mat_mult_dense_diag(c, k, k, dsum_vars, tmp_mat)
        call mat_mult(tmp_mat, k, k, ct, k, dsums_covi)

ccc        call prin2('vars*', dsum_vars, k)
ccc        call prin2('vars*', dsums_covi, 1)

c
c        recover matrix E[x*x^t]
c
        call mat_mult(c, k, k, wwti, k, tmp_mat)
        call mat_mult(tmp_mat, k, k, ct, k, xxti)

        call get_xs_from_ws(dsumsi,k,k1,k2,vt,t, dsum_xsi)
ccc        call prin2('xs 2*', dsums_xsi, k)
ccc        call prin2('dsumi*', dsumi, 1)

c
c        due to underflow issues, integrate over theta by computing
c        sum of the form \sum_i exp(fi)*gi such that at the end
c        we have an expression exp(fm)*dsum
c
        fi = fmi
        wt = wht_theta
        if (fi .gt. fm) then
          dsum = dsum * exp(fm-fi) + dsumi*wt

          do ijk=1,k
          dsums(ijk) = dsums(ijk) * exp(fm-fi) + dsum_xsi(ijk)*wt
          enddo

          dsums(k+1) = dsums(k+1) * exp(fm-fi) + dsumsi(k+1)*wt
          dsums(k+2) = dsums(k+2) * exp(fm-fi) + dsumsi(k+2)*wt
          dsums(k+3) = dsums(k+3) * exp(fm-fi) + dsumsi(k+3)*wt

          ss1 = ss1 * exp(fm-fi) + ss1i*wt
          ss2 = ss2 * exp(fm-fi) + ss2i*wt
          ss3 = ss3 * exp(fm-fi) + ss3i*wt

          do i1=1,k
          do i2=1,k
          tmp = dsums_covi(i1,i2)+xxti(i1,i2)
          dsums_cov(i1,i2)=dsums_cov(i1,i2)*exp(fm-fi) + tmp*wt
          enddo
          enddo

          fm = fi
        else
          dsum = dsum + exp(fi-fm) * dsumi*wt

          do ijk=1,k
          dsums(ijk) = dsums(ijk) + exp(fi-fm)*dsum_xsi(ijk)*wt
          enddo

          dsums(k+1) = dsums(k+1) + exp(fi-fm)*dsumsi(k+1)*wt
          dsums(k+2) = dsums(k+2) + exp(fi-fm)*dsumsi(k+2)*wt
          dsums(k+3) = dsums(k+3) + exp(fi-fm)*dsumsi(k+3)*wt

          ss1 = ss1 + exp(fi-fm)*ss1i*wt
          ss2 = ss2 + exp(fi-fm)*ss2i*wt
          ss3 = ss3 + exp(fi-fm)*ss3i*wt

          do i1=1,k
          do i2=1,k
          tmp = dsums_covi(i1,i2)+xxti(i1,i2)
          dsums_cov(i1,i2)=dsums_cov(i1,i2) + exp(fi-fm)*tmp*wt
          enddo
          enddo

        endif

 130    continue

c
c        esimate error by looking at the coefficients of our
c        function in a 2d legendre expansion
ccc        call get_err(nn_theta, nn, fs, gs)

c
c        scale first and second moments by normalizing constant
c
        do i=1,k+3
        dsums(i)=dsums(i)/dsum
        enddo

        do i=1,k
        do j=1,k
        dsums_cov(i, j)=dsums_cov(i, j)/dsum
        enddo
        enddo

c
c        dsums_cov now contains E[xx^t], adjust to get covariance
c
        do i=1,k
        do j=1,k
        dsums_cov(i, j)=dsums_cov(i, j) - dsums(i)*dsums(j)
        enddo
        enddo

c
c        get stds
c
        do i=1,k
        stds(i) = sqrt(dsums_cov(i, i))
        enddo

c
c        get variances of sig1, sig2, sig3
c
        ss1 = ss1/dsum
        ss2 = ss2/dsum
        ss3 = ss3/dsum
        stds(k+1) = sqrt(ss1 - dsums(k+1)**2)
        stds(k+2) = sqrt(ss2 - dsums(k+2)**2)
        stds(k+3) = sqrt(ss3 - dsums(k+3)**2)

c
c        plot
c
ccc        call prin2('fs*', fs, nn_theta*nn)
ccc        call plot_heatmap(nn_theta, nn, fs)

        return
        end
c
c
c
c
c
        subroutine get_err(nnt, nn, fs, gs)
        implicit real *8 (a-h,o-z)
        real*8 fs(nnt, nn), gs(nnt, nn), hs(nnt, nn), u(100 000),
     1     v(100 000), x(10 000), u_theta(100 000), v_theta(100 000),
     1     coefs(nnt, nn), coefs_log(nnt, nn), whts(100 000),
     1     ut(100 000), tmp_mat(100 000)

c
c        this subroutine is not finished! it is in the middle of
c        development
c

        call max_vec(nnt*nn, fs, fmax, ind)

        do i=1,nnt
        do j=1,nn
        hs(i, j) = fs(i, j) - fmax + log10(gs(i, j))
        hs(i, j) = exp(fs(i, j) - fmax)*gs(i, j)
        enddo
        enddo
ccc        call plot_heatmap(nnt, nn, hs)

c
c        legendre expansion
c
        itype=2
        call legeexps(itype,nnt,x,u_theta,v_theta,whts)
        call legeexps(itype,nn,x,u,v,whts)

        call mat_trans(u, nn, nn, ut)

        call mat_mult(u_theta, nnt, nnt, hs, nn, tmp_mat)
        call mat_mult(tmp_mat, nnt, nn, ut, nn, coefs)

ccc        call prin2('coefs*', coefs, nnt*nn)

        dsum = 0.0d0
        ndsum = 0
        tot_coef = 0.0d0
        do i=1,nnt
        do j=1,nn
        coefs_log(i, j) = log10(abs(coefs(i,j)))
        tot_coef = tot_coef + abs(coefs(i,j))

        if ((i .gt. nnt-3) .or. (j .gt. nn-5)) then
          dsum = dsum + abs(coefs(i,j))
          ndsum = ndsum + 1
        endif

        enddo

        if (i .eq. nnt) then
          call prin2('coefs(i, j)*', coefs(i, j), 1)
        endif

        enddo

ccc        call prin2('coefs_log*', coefs_log, nnt*nn)
        call prin2('coefs*', coefs, nnt*nn)
        call plot_heatmap(nnt, nn, coefs_log)

c
c        error estimate
c
ccc        call prinf('ndsum*', ndsum, 1)
        avg_tail_coef = dsum / (nn*3 + nnt*5 - 3*5)
        call prin2('avg_tail_coef*', avg_tail_coef, 1)
        err = avg_tail_coef / tot_coef

        call prin2('err*', err, 1)

        fmax = -1.0d250
        do i=1,nn
        if (fmax .lt. coefs(nnt, i)) fmax = coefs(nnt, i)
        enddo

        do i=1,nnt
        if (fmax .lt. coefs(i, nn)) fmax = coefs(i, nn)
        enddo

ccc        call prin2('fmax*', fmax, 1)
ccc        call prin2('fmax2*', coefs(nnt, nn), 1)

        return
        end
c
c
c
c
c
        subroutine get_xs_to_ws_matrix(k,k1,k2,v,t,a,at)
        implicit real *8 (a-h,o-z)
        real*8 s1(k),v(k,k),a(*),at(*)

c
c        construct matrix that goes from diagonal coordinate
c        system (depending on t) to original coordinate system
c

        do i=1,k1
        s1(i) = cos(t)
        enddo

        do i=1,k2
        s1(i+k1) = sin(t)
        enddo

        call mat_mult_diag_dense(s1, k, v, a)
        call mat_trans(a, k, k, at)

        return
        end
c
c
c
c
c
        subroutine get_phi1_fmax(nn, n, k, t, ysmall, ys2, s, s2, resid,
     1     phi1i, fmax)
        implicit real *8 (a-h,o-z)
        real*8 rhos(nn+10),whts_rhos(nn+10),phis(nn+10),
     1     whts_phis(nn+10), ysmall(*),ys2(*),s(*),s2(*),
     1     fmaxs(nn+10)

c
c        find the upper integration bound of the integral with
c        respect to phi. do this by evaluating the integral on
c        a sparse grid and checking when the value decreases below
c        10^{-18} of its maximum. also return the maximum value of
c        the integrand
c

        done = 1.0d0
        pi = 4*atan(done)

c
c        lay down nodes
c
        phi0 = 0.0d0
        phi1 = pi/2.0d0
        call lege_nodes_whts(nn, phi0, phi1, phis, whts_phis)

c
c        tabulate function
c
        fmax = 0.0d0
        do 119 j=1,nn
        phi = phis(j)
        call compute_beta_alpha_prefact(phi,ysmall,ys2,s,s2,
     1      n,k,resid,alpha,beta,prefact,exp_fact)

c
c        get maximum of integrand and use that maximum for the
c        first value of phi for scaling
c
        dn = n - 2
        rho_max = sqrt(1.0d0/2.0d0*(sqrt(4.0d0*exp_fact + dn*dn) - dn))
        call eval_logdens_rho(rho_max,n,exp_fact,fmaxmax)
        fmaxs(j) = fmaxmax

c
c        compute integral in rho
c
        call eval_rho_int(nn, n, k, exp_fact,
     1     dsum_rho, dsum_rho1, dsum_rho2, fmaxi)
        fi = log(dsum_rho) + prefact + fmaxi - fmaxs(1)

        if (fi .gt. fmax) fmax = fi

        thresh = log(1.0d-22)
        if (fi - fmax .lt. thresh) then
          phi1i = phi
          goto 121
        endif
 119    continue

c
c        must use full interval
c
        phi1i = pi/2.0d0

 121    continue

c
c        we're going to use the log maximum of the integrand later,
c        so return that value
c
        fmax = fmax + fmaxs(1)

        return
        end
c
c
c
c
c
        subroutine eval_rho_int(nn, n, k, exp_fact,
     1     dsum_rho, dsum_rho1, dsum_rho2, fmax)
        implicit real *8 (a-h,o-z)
        real*8 whts_rhos(nn),rhos(nn)

c
c        compute the integral from 0 to infinity of
c         f(rho), rho*f(rho), and rho**2*f(rho)
c
c        where
c
c         f(rho) = exp(-rho**2/2 - exp_fact/(2*rho**2))/rho**(n-2)
c

c
c        first find integration bounds
c
        call get_rho0(exp_fact, n, rho0)
        call get_rho1(exp_fact, n, rho1)
        call lege_nodes_whts(nn, rho0, rho1, rhos, whts_rhos)

c
c        get maximum of integrand
c
        dn = n - 2
        rho_max = sqrt(1.0d0/2.0d0*(sqrt(4.0d0*exp_fact + dn*dn) - dn))
        call eval_logdens_rho(rho_max,n,exp_fact,fmax)

c
c        compute integral
c
        dsum_rho = 0
        dsum_rho1 = 0
        dsum_rho2 = 0
        do 110 i=1,nn
        rho = rhos(i)
        call eval_logdens_rho(rho,n,exp_fact,f)
        wht_rho = whts_rhos(i)

        dsum_rho = dsum_rho + exp(f-fmax)*wht_rho
        dsum_rho1 = dsum_rho1 + rho*exp(f-fmax)*wht_rho
        dsum_rho2 = dsum_rho2 + rho**2*exp(f-fmax)*wht_rho

 110    continue

        return
        end
c
c
c
c
c
        subroutine compute_mjs(k,s2,ys,phi,vmoms)
        implicit real *8 (a-h,o-z)
        dimension vmoms(*),s2(*),ys(*)

c
c        compute conditional expecations of gaussian
c

        sp = sin(phi)
        sp2= sp*sp
        cp = cos(phi)**2

        do i=1,k
        vmoms(i) = ys(i)*sqrt(s2(i))*sp2/
     1          (s2(i)*sp2+cp)
        enddo

        return
        end
c
c
c
c
c
        subroutine eval_logdens_rho(rho,n,exp_fact,val)
        implicit real *8 (a-h,o-z)

c
c        evaluate the log of the integrand of the inner integral
c

        val = -exp_fact/2.0d0/(rho*rho)
        val = val -rho*rho/2.0d0
        val = val - (n-2)*log(rho)

        return
        end
c
c
c
c
c
        subroutine compute_beta_alpha_prefact(phi,ys,ys2,s,s2,
     1      n,k,resid,alpha,beta,prefact,exp_fact)
        implicit real *8 (a-h,o-z)
        dimension ys(*),ys2(*),s(*),s2(*)

c
c        for convenience, compute quantities that appear in
c        several locations in the integrand. prefact is the
c        part of the integrand that doesn't depend on rho
c

        done = 1.0d0
        pi   = atan(done)*4

        sp = sin(phi)**2
        cp = cos(phi)**2

c
c        alpha
c
        alpha = 0

        do i=1,k
        alpha = alpha - log(cp+sp*s2(i))
        enddo

        alpha = alpha/2 + k/2.0d0*log(2*pi)

c
c        beta
c
        beta = 0

        do i=1,k
        beta = beta + ys2(i)/(s2(i)*sp+cp)
        enddo

c
c        prefact
c
        prefact = -(n-k)/2.0d0*log(cp)+alpha + log(sin(phi))
        exp_fact= resid/cp+beta

        return
        end
c
c
c
c
c
        subroutine get_xs_from_ws(ws,k,k1,k2,vt,t, xs)
        implicit real *8 (a-h,o-z)
        dimension xs(*),ws(*),vt(k,k)

c
c        convert expectations back from the diagonal coordinate (ws)
c        system to the original coordinate system (xs)
c
        do i=1,k
        z = 0
        do j=1,k
        z = z + vt(j,i)*ws(j)
        enddo
        xs(i) = z
        enddo

        do i=1,k1
        xs(i) = xs(i)*cos(t)
        enddo

        do i=1,k2
        xs(i+k1) = xs(i+k1)*sin(t)
        enddo

        return
        end
c
c
c
c
c
        subroutine entry_sqr(xin,xout,k)
        implicit real *8 (a-h,o-z)
        dimension xin(*),xout(*)
c
c        square the entries of an array
c
        do i=1,k
        xout(i) = xin(i)**2
        enddo

        return
        end
c
c
c
c
c
        subroutine y_proj(y,us,ss,n,k,ysmall,resid)
        implicit real *8 (a-h,o-z)
        dimension y(n),us(n,k),ss(*),ysmall(k)

c
c        compute the residual of the least squares solution
c

        rr = 0

        do i=1,k
        ys = 0
        do j=1,n
        ys = ys + y(j)*us(j,i)
        enddo
        ysmall(i) = ys
        rr = rr + ys*ys
        enddo

        yy = 0
        do i=1,n
        yy = yy + (y(i))**2
        enddo

        resid = yy-rr

        return
        end
c
c
c
c
c
        subroutine rescale_a(t,a,asca,n,k,k1,k2)
        implicit real *8 (a-h,o-z)
        dimension a(n,k),asca(n,k)

c
c        multiply a by diagonal matrix
c

        tc = cos(t)
        ts = sin(t)

        do i=1,n

        do j=1,k1
        asca(i,j) = a(i,j)*tc
        enddo

        do j=(k1+1),k
        asca(i,j) = a(i,j)*ts
        enddo

        enddo

        return
        end
c
c
c
c
c
        subroutine get_rho0(ci, n, rho0)
        implicit real *8 (a-h,o-z)
c
c        find the lower bound of integration of the inner integral
c        using bisection and the formula for the maximum of the
c        integrand
c

c
c        compute maximum of integrand
c
        dn = n - 2
        rho_max = sqrt(1/2.0*(sqrt(4*ci + dn**2) - dn))
        call eval_logdens_rho(rho_max,n,ci,fmax)
        fmax_log10 = fmax / log(10.0d0)

c
c        bisection
c
        xmin = 1.0d-16
        xmax = rho_max
        do 110 i=1,10
        rho0 = (xmin + xmax)/2.0d0
        call eval_logdens_rho(rho0,n,ci,f)
        f_log10 = f / log(10.0d0)
        if (f_log10 .lt. fmax_log10-18) xmin = rho0
        if (f_log10 .ge. fmax_log10-18) xmax = rho0
 110    continue

ccc        call eval_logdens_rho(rho0, n, ci, f)
ccc        call prin2('f at left endpoint*', exp(f), 1)

        return
        end
c
c
c
c
c
        subroutine get_rho1(ci, n, rho1)
        implicit real *8 (a-h,o-z)

c
c        find the upper bound of integration of the inner integral
c        using bisection and the formula for the maximum of the
c        integrand and the second derivative at the maximum
c

c
c        compute maximum of integrand and second derivative at maximum
c

        dn = n - 2
        rho_max = sqrt(1/2.0*(sqrt(4*ci + dn**2) - dn))
        call eval_logdens_rho(rho_max,n,ci,fmax)
        fmax_log10 = fmax / log(10.0d0)

        dder2 = -1-3*ci/rho_max**4 + (n-2)/rho_max**2
        sd = 1/sqrt(-dder2)

c
c        bisection
c
        xmin = rho_max
        xmax = rho_max + 40*sd
        do 120 i=1,10
        rho1 = (xmin + xmax)/2.0d0
        call eval_logdens_rho(rho1,n,ci,f)

        f_log10 = f / log(10.0d0)
        if (f_log10 .lt. fmax_log10-18) xmax = rho1
        if (f_log10 .ge. fmax_log10-18) xmin = rho1
 120    continue

ccc        call prin2('rho1*', rho1, 1)
ccc        call eval_logdens_rho(rho1, n, ci, f)
ccc        if (f .gt. -20) then
ccc          call prin2('bad!*', 1.0d0, 1)
ccc          call prin2('f*', f, 1)
ccc        endif

        return
        end
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
        subroutine plot_heatmap(nn1, nn2, fs)
        implicit real *8 (a-h,o-z)
        real*8 fs(1)
        data in/11/

        tol = -20.0d0

        call max_vec(nn1*nn2, fs, f1, ind)
ccc        call prin2('max*', f1, 1)

        do 320 i=1,nn1*nn2
        fs(i) = fs(i) - f1
        if (fs(i) .lt. tol) fs(i) = tol
ccc        call prin2('fs*', fs, nnodes**2)
 320    continue

c
c        if this function gets called twice, nothing gets
c        plotted, so change iw
c
        iw = in
        call pyimage(iw,nn1,nn2,fs,'title* ')

        in = in + 1

        return
        end
c
c
c
c
c
        subroutine cc_arr(n, x, y)
        implicit real*8 (a-h,o-z)
        real*8 x(1),y(1)

        do 100 i=1,n
        y(i) = x(i)
 100    continue

        return
        end
c
c
c
c
c
        subroutine mat_vec_mult(a, n, m, x, y)
        implicit real*8 (a-h,o-z)
        real*8 a(n, 1), x(1), y(1)

        do 450 i = 1, n
        y(i) = 0
        do 425 j = 1, m
        y(i) = y(i) + x(j) * a(i, j)
 425    continue
 450    continue

        return
        end
c
c
c
c
c
        subroutine mat_mult(a, n, m, b, k, c)
        implicit real*8 (a-h,o-z)
        real*8 a(n, 1), x(1), y(1), b(m, 1), c(n, 1)

        do 1075 i=1,n
        do 1050 j=1,k
        c(i,j) = 0
        do 1025 kk = 1,m
        c(i,j) = c(i,j) + a(i, kk) * b(kk, j)
 1025   continue
 1050   continue
 1075   continue

        return
        end
c
c
c
c
c
        subroutine inner_prod(u, v, n, y)
        implicit real *8 (a-h,o-z)
        real*8 u(1), v(1)

        y = 0
        do 1050 i=1,n
        y = y + u(i) * v(i)
 1050   continue

        return
        end
c
c
c
c
c
        subroutine mat_trans(u, n, m, ut)
        implicit real *8 (a-h,o-z)
        real*8 u(n, 1), ut(m, 1)

        do 1075 i=1,n
        do 1070 j=1,m
        ut(j, i) = u(i, j)
 1070   continue
 1075   continue

        return
        end
c
c
c
c
c
        subroutine theta_lege_nodes_whts(nn, t0, t1, ts, whts)
        implicit real *8 (a-h,o-z)
        real*8 ts(*), whts(*), u(1), v(1), x(nn+10),
     1      xleg(1000),wleg(1000)
        data nloc/0/
        save

        if (nn .ne. nloc) then
        nloc   = nn
        ifsave = 0

        itype=1
        call legeexps(itype,nn,xleg,u,v,wleg)
        endif

        do 60 i=1,nn
        tmp = (xleg(i) + 1)/2.0d0
        ts(i) = t0 + (t1-t0)*tmp
        whts(i) = wleg(i)*(t1-t0)/2.0d0
 60     continue

ccc        call prin2('wleg*', wleg, nn)
ccc        call prin2('xleg*', xleg, nn)

        return
        end
c
c
c
c
c
        subroutine lege_nodes_whts(nn, t0, t1, ts, whts)
        implicit real *8 (a-h,o-z)
        real*8 ts(*), whts(*), u(1), v(1), x(nn+10),
     1      xleg(1000),wleg(1000)
        data nloc/0/
        save

        if (nn .ne. nloc) then
        nloc   = nn
        ifsave = 0

        itype=1
        call legeexps(itype,nn,xleg,u,v,wleg)
        endif

        do 60 i=1,nn
        tmp = (xleg(i) + 1)/2.0d0
        ts(i) = t0 + (t1-t0)*tmp
        whts(i) = wleg(i)*(t1-t0)/2.0d0
 60     continue

ccc        call prin2('wleg*', wleg, nn)
ccc        call prin2('xleg*', xleg, nn)

        return
        end
c
c
c
c
c
        subroutine max_vec(n, v, f1, ind)
        implicit real *8 (a-h,o-z)
        real*8 v(1)

        f1 = v(1)
        ind = 1
        do 100 i=1,n-1
        if (v(i+1) .gt. f1) then
          f1 = v(i+1)
          ind = i+1
        endif
 100    continue

        return
        end
c
c
c
c
c
        subroutine pad_mat(b,n,m,nn,c)
        implicit real *8 (a-h,o-z)
        dimension b(n,m),c(nn,m)

        do i=1,nn
        do j=1,m
        c(i,j) = 0
        enddo
        enddo

        do i=1,n
        do j=1,m
        c(i,j) = b(i,j)
        enddo
        enddo

        return
        end
c
c
c
c
c
        subroutine read_params(csv_file, n, k1, k2, a, y)
        implicit real *8 (a-h,o-z)
        real*8 a(*), y(*)
        character(*) csv_file
c
c        read the csv file written in rstudio by test43.R
c
ccc        csv_file = 'dense_params.csv'
        open (2, file=csv_file)
c
c        read matrix size
c
 100    format(i12)
        read(2, 100) n
ccc        call prinf('n*', n, 1)
        read(2, 100) k1
ccc        call prinf('k*', k, 1)
        read(2, 100) k2
ccc        call prinf('k*', k, 1)
c
c        read a
c
        k = k1+k2
ccc        call prinf('nk*', n*k, 1)
 200    format(f12.6)
        do 250 i=1,n*k
        read(2, 200) a(i)
ccc        call prinf('i*', i, 1)
ccc        call prin2('a(i)*', a(i), 1)
 250    continue

ccc        call prin2('a*', a, n*k)
c
c        read y
c
        do 350 i=1,n
        read(2, 200) y(i)
 350    continue

ccc        call prin2('y*', y, n)
        return
        end
c
c
c
c
c
        subroutine max_vec_abs(n, v, f1, ind)
        implicit real *8 (a-h,o-z)
        real*8 v(1)

        f1 = abs(v(1))
        ind = 1
        do 100 i=1,n-1
        if (abs(v(i+1)) .gt. f1) then
          f1 = abs(v(i+1))
          ind = i+1
        endif
 100    continue

        return
        end
c
c
c
c
c
        subroutine mat_mult_dense_diag(a, n, k, s, c)
        implicit real*8 (a-h,o-z)
        real*8 a(n, *), s(*), c(n, *)

        do 1075 i=1,n
        do 1050 j=1,k
        c(i,j) = s(j)*a(i,j)
 1050   continue
 1075   continue

        return
        end
c
c
c
c
c
        subroutine mat_mult_diag_dense(s, n, a, c)
        implicit real*8 (a-h,o-z)
        real*8 a(n, 1), s(1), c(n, 1)

        do 1075 i=1,n
        do 1050 j=1,n
        c(i,j) = s(i)*a(i,j)
 1050   continue
 1075   continue

        return
        end
c
c
c
c
c
        subroutine compute_ata(a,k,aout)
        implicit real *8 (a-h,o-z)
        dimension a(k,k),aout(k,k)

        do i=1,k
        do j=1,k
        z = 0
        do ijk=1,k
        z = z + a(ijk,i)*a(ijk,j)
        enddo
        aout(i,j) = z
        enddo
        enddo

        return
        end
c
c
c
c
c
        subroutine sym_rescale_a(t,a,asca,k,k1,k2)
        implicit real *8 (a-h,o-z)
        dimension a(k,k),asca(k,k)

        tc = cos(t)
        ts = sin(t)

        tcc= tc*tc
        tss= ts*ts
        tcs= tc*ts

        do i=1,k1
        do j=1,k1
        asca(i,j) = a(i,j)*tcc
        enddo
        enddo

        do i=k1+1,k
        do j=k1+1,k
        asca(i,j) = a(i,j)*tss
        enddo
        enddo

        do i=1,k1
        do j=k1+1,k
        tt = a(i,j)*tcs
        asca(i,j) = tt
        asca(j,i) = tt
        enddo
        enddo

        return
        end
c
c
c
c
c
        subroutine read_means(csv_file, m, dexps)
        implicit real *8 (a-h,o-z)
        real*8 dexps(*)
        character(100) csv_file

c
c        read the posterior means computed by stan from
c        text file
c

        open (2, file = csv_file)

c
c        read posterior expectations
c

 200    format(f12.6)

        do 250 i=1,m+3
        read(2, 200) dexps(i)
 250    continue

        return
        end
c
c
c
c
c
        subroutine dd_abs_max(v1, v2, n, dd_max)
        implicit real*8 (a-h,o-z)
        real*8 v1(*), v2(*)

        dd_max = 0.0d0
        do i=1,n
        dd = abs(v1(i) - v2(i))
        if (dd .gt. dd_max) dd_max = dd
        enddo

        return
        end
c
c
c
c
c
        subroutine dd_abs_max_rel(v1, v2, n, dd_max)
        implicit real*8 (a-h,o-z)
        real*8 v1(*), v2(*)

        dd_max = 0.0d0
        do i=1,n
        dd = abs((v1(i) - v2(i))/v1(i))
        if (dd .gt. dd_max) dd_max = dd
        enddo

        return
        end
c
c
c
c
c
        subroutine write_exps_stds(file_out, k, dsums, stds)
        implicit real*8 (a-h,o-z)
        real*8 dsums(*), stds(*)
        character(*) file_out

        open (2, file=file_out)

 210    format(f22.16,','f22.16)
        do i=1,k+3
        write(2, 210) dsums(i), stds(i)
        enddo

        return
        end
c
c
c
c
c
        subroutine read_exps_stds(file_in, k, dsums, stds)
        implicit real*8 (a-h,o-z)
        real*8 dsums(*), stds(*)
        character(*) file_in

        open (2, file=file_in)

 210    format(f22.16,X,f22.16)
        do i=1,k+3
        read(2, 210) dsums(i), stds(i)
        enddo

        return
        end
c
c
c
c
c
        subroutine set_zero(n, v)
        implicit real *8 (a-h,o-z)
        real*8 v(*)

        do i=1,n
        v(i) = 0.0d0
        enddo

        return
        end
c
c
c
c
c
        subroutine div_arr(n, alph, v)
        implicit real *8 (a-h,o-z)
        real*8 v(*)

        do i=1,n
        v(i) = v(i) / alph
        enddo

        return
        end
c
c
c
c
c
        subroutine write_dd_exps_stds(file_out, k, dd, dsums, stds)
        implicit real*8 (a-h,o-z)
        real*8 dsums(*), stds(*)
        character(*) file_out

ccc        csv_file = 'expectations.dat'
        open (2, file=file_out)

 210    format(f22.16,','f22.16)

        write(2, 210) dd

        do i=1,k+3
        write(2, 210) dsums(i), stds(i)
        enddo

        return
        end

ext: f

max_stars_repo_name: pgree/fastNoNo

max_issues_repo_name: pgree/fastNoNo

max_forks_repo_name: pgree/fastNoNo


text:
program InterpDataFromVTKtoNCLNearestNeighbor

use NumberKindsModule
use SphereGeomModule
use STRIPACKInterfaceModule
use SSRFPACKInterfaceModule

implicit none
!
! file IO
!
character(len=MAX_STRING_LENGTH) :: namelistfile = 'InterpVTKtoNCLNearestNeighbor.namelist'
character(len=MAX_STRING_LENGTH) :: inputfile
character(len=MAX_STRING_LENGTH) :: outputPrefix, outputFile
real(kreal) :: programStart, programEnd
integer(kint) :: readStat, writeStat
!
! string processing
!
character(len=MAX_STRING_LENGTH) :: lineIn, scratch1, scratch2, foundName
integer(kint) :: lineNumber, scalarCount, polygonHeaderLine, pointDataLine
integer(kint) :: i, j, k, i1, i2, i3, ndim, ii, jj, kk
logical(klog) :: keepGoing, neighborsFound
!
! data conversion
!
real(kreal), allocatable :: xyz(:,:), lon(:), lat(:), scalars(:), outputData(:,:)
integer(kint), allocatable :: nearest(:,:)
integer(kint) :: nTotal, nTri, nLat, nLon
real(kreal) :: degreeSpacing, radSpacing, currentMin, testDist, xTest(3)

namelist /filenames/ inputfile, outputPrefix, degreeSpacing

open(unit=READ_UNIT, file=namelistfile, status='OLD',action='READ',iostat = readstat)
	if ( readstat /= 0 ) stop "ERROR : cannot open namelist file"
	read(READ_UNIT,nml=filenames)
close(READ_UNIT)

! read vtk metadata
open(unit=READ_UNIT, file=inputfile, status='old', action='READ',iostat=readStat)
	if ( readStat /= 0 ) stop "ERROR opening VTK input file."
keepGoing = .TRUE.
lineNumber = 0
scalarCount = 0
polygonHeaderLine = -1
pointDataLine = -1
print '("Getting VTK File MetaData ...")'
print '("VTK FILE INFO :")'
do while ( keepGoing )
	! read next line into lineIn
	read(READ_UNIT,'(A)',iostat = readStat) lineIn
	lineNumber = lineNumber + 1
	if ( readStat > 0 ) then
		print '("ERROR : problem encountered at line ",I8," .")' , lineNumber
		keepGoing = .FALSE.
		stop
	elseif ( readStat < 0 ) then
		print '("End of file found at line ",I8,".")', lineNumber
		keepGoing = .FALSE.
	else
		if ( lineNumber == 5) then
			! get number of particles
			lineIn = adjustl(lineIn)
			lineIn = trim(lineIn)
			i1 = scan(lineIn,' ')
			i2 = scan(lineIn,' ',back=.TRUE.)
			scratch1 = lineIn(i1+1:i2)
			read(scratch1,'(I8)') nTotal
			print '("... found ",I8," total particles ...")', nTotal
			polygonHeaderLine = 5 + nTotal + 1
		elseif ( lineNumber == polygonHeaderLine ) then
			! count the number of triangles
			lineIn = adjustl(lineIn)
			lineIn = trim(lineIn)
			i1 = scan(lineIn,' ')
			i2 = scan(trim(lineIn),' ',back=.TRUE.)
			scratch1 = lineIn(i1:i2)
			scratch1 = trim(adjustl(scratch1))
			read(scratch1,'(I8)') nTri
			print '("... found ",I8," triangles ...")', nTri
			print '("... so there are = ",I8," active panels ...")', nTri/4
			pointDataLine = polygonHeaderLine + nTri + 1
		elseif ( lineNumber > pointDataLine .AND. lineNumber > 5 ) then
			i1 = 0
			lineIn = adjustl(lineIn)
			i1 = scan(lineIn,'S')
			if ( i1 == 1 ) scalarCount = scalarCount + 1
		endif
	endif
enddo
print '("... found ", I8, " data fields.")', scalarCount

print '("Allocating memory ...")'

allocate(xyz(3,nTotal))
xyz = 0.0_kreal
allocate(scalars(nTotal))
scalars = 0.0_kreal
nLon = floor(360.0_kreal / degreeSpacing)
nLat = nLon/2 + 1
allocate(lat(nLat))
allocate(lon(nLon))
allocate(outputData(nLat,nLon))
allocate(nearest(nLat,nLon))
neighborsFound = .FALSE.
outputData = 0.0_kreal
radSpacing = degreeSpacing * PI / 180.0_kreal
do j = 1, nLon
	lon(j) = (j-1)*radSpacing
enddo
do j = 1, nLat
	lat(j) = -PI/2.0_kreal + (j-1)*radSpacing
enddo

print '("Reading particle and field data ...")'

rewind(READ_UNIT)
keepGoing = .TRUE.
lineNumber = 0
j = 1
k = 1
nDim = 1
do while (keepGoing)
	read(READ_UNIT,'(A)',iostat = readStat) lineIN
	lineNumber = lineNumber + 1
	if ( readStat > 0 ) then
		print '("ERROR : problem encountered at line ",I8," .")' , lineNumber
		keepGoing = .FALSE.
		stop
	elseif ( readStat < 0 ) then
		print '("End of file found at line ",I8,".")', lineNumber
		keepGoing = .FALSE.
	else
		if ( lineNumber > 5 .AND. lineNumber < polygonHeaderLine ) then
			! read vertex coordinate data
			read ( lineIn, '(3D)') xyz(:,j)
			j = j + 1
		endif
		if ( lineNumber == polygonHeaderLine ) then
			j = 1 ! reset counter
		! DEBUG	print '("Max lat = ",F5.3)', maxval(lat)
		endif
		if ( lineNumber > polygonHeaderLine .AND. lineNumber > pointDataLine ) then
			lineIn = adjustl(lineIn)
			i1 = scan(lineIn,'S')
			if ( i1 == 1 ) then
				print '("... processing data field ",I2," of ", I2,".")', k, scalarCount
				k = k + 1
				j = 1
				i2 = scan(lineIn,' ')
				scratch1 = lineIn(i2+1:len_trim(lineIn))
				scratch1 = adjustl(scratch1)
				i3 = scan(scratch1,' ')
				foundName = scratch1(1:i3)
				scratch2 = lineIn(len_trim(lineIN)-1:len_trim(lineIn))
				read(scratch2,'(I)') nDim
				write(6,'(A,A,A,I1)') "... data field ", trim(foundName), " has dimension ", ndim
			endif
			i1 = 0
			i1 = scan(lineIn,'L')
			if ( i1 == 0 ) then
				! read scalar data
				lineIn = trim(lineIn)
				if ( nDim == 1 ) then
					read(lineIn,'(D)') scalars(j)
					j = j + 1
				endif
			endif
			if ( j == nTotal + 1 .AND. ndim == 1) then

				if ( .NOT. neighborsFound ) then
					write(6,*) "... finding closest particles to output points."
					currentMin = PI * EARTH_RADIUS
					do jj = 1, nLon
						do ii = 1, nLat
							currentMin = PI * EARTH_RADIUS
							xTest = EARTH_RADIUS * [ cos(lat(ii)) * cos(lon(jj)), cos(lat(ii))*sin(lon(jj)), sin(lat(ii))]
							do kk = 1, nTotal
								testDist = SphereDistance( xTest, xyz(:,kk) )
								if ( testDist < currentMin ) then
									nearest(ii,jj) = kk
									currentMin = testDist
								endif
							enddo
						enddo
					enddo
					write(6,*) "... nearest neighbors found."
					neighborsFound = .TRUE.
				endif

				do jj = 1, nLon
					do ii = 1, nLat
						outputData(ii,jj) = scalars(nearest(ii,jj))
					enddo
				enddo

				write(outputFile,'(A,A,A,A)') trim(outputPrefix),'_Interp_',trim(foundName),'.txt'
				open(unit=WRITE_UNIT_1,file=outputfile,status='REPLACE',action='WRITE',iostat=writestat)
					if ( writeStat /= 0 ) stop "ERROR writing data file."
				do i1 = 1, nLat
					do i2 = 1, nLon -1
						write(WRITE_UNIT_1,'(F24.8)',advance='NO') outputData(i1,i2)
					enddo
					write(WRITE_UNIT_1,'(F24.8)',advance='YES') outputData(i1,nLon)
				enddo
				close(WRITE_UNIT_1)
			endif
		endif
	endif
enddo
close(READ_UNIT)
call cpu_time(programEnd)
print *, " PROGRAM COMPLETE "
print '("Elapsed time : ", f7.3, " minutes.")', (programEnd - programStart)/60.0

deallocate(nearest)
deallocate(outputData)
deallocate(scalars)
deallocate(lat)
deallocate(lon)
deallocate(xyz)

end program

max_stars_repo_name: pbosler/LPPM

max_issues_repo_name: pbosler/LPPM

max_forks_repo_name: pbosler/LPPM


text:
! Calculating Pi with using method:
!     integral 4.0/(1+x^2) dx = p
! 	   OpenMP Version 1
program pi_omp_v1
implicit none
integer*8 :: num_steps = 1000000000
integer i,tid,nthreads,omp_get_num_threads,omp_get_thread_num
real*8 :: PI_ref=3.1415926535897932,step,x,pi
real*8 :: summation=0.0
real*8   :: start,finish,omp_get_wtime
!call cpu_time(start)
step = 1/(1.0 * num_steps)
start=omp_get_wtime()
!$omp parallel private(x,tid) firstprivate(summation) shared(pi,nthreads,step)
	tid=omp_get_thread_num();
!$omp single
	nthreads=omp_get_num_threads();
!$omp end single
	do i =tid,num_steps-1,nthreads
		x=((i-1)+0.5) * step;
		summation=summation + (4.0/(1.0+(x*x)))
	end do
!$omp atomic update
	pi=pi+(step*summation)

!$omp end parallel
!call cpu_time(finish)
finish=omp_get_wtime()
	write(*,*) "Number of threads= ",nthreads
	write(*,*) "Estimated value of pi =",pi
	write(*,*) "Error =", abs(PI_ref - pi)
	write(*,*) "Compute time=", (finish-start)," seconds"
end

max_stars_repo_name: pelahi/Develop-with-OpenMP

max_issues_repo_name: pelahi/Develop-with-OpenMP

max_forks_repo_name: pelahi/Develop-with-OpenMP


text:
C***********************************************************************
C***********************************************************************
C
C   Version:        0.3
C   Last modified:  December 27, 1994
C   Authors:        Esmond G. Ng and Barry W. Peyton
C
C   Mathematical Sciences Section, Oak Ridge National Laboratory
C
C***********************************************************************
C***********************************************************************
C******     BFINIT ..... INITIALIZATION FOR BLOCK FACTORIZATION   ******
C***********************************************************************
C***********************************************************************
C
C   PURPOSE:
C       THIS SUBROUTINE COMPUTES ITEMS NEEDED BY THE LEFT-LOOKING
C       BLOCK-TO-BLOCK CHOLESKY FACTORITZATION ROUTINE BLKFCT.
C
C   INPUT PARAMETERS:
C       NEQNS           -   NUMBER OF EQUATIONS.
C       NSUPER          -   NUMBER OF SUPERNODES.
C       XSUPER          -   INTEGER ARRAY OF SIZE (NSUPER+1) CONTAINING
C                           THE SUPERNODE PARTITIONING.
C       SNODE           -   SUPERNODE MEMBERSHIP.
C       (XLINDX,LINDX)  -   ARRAYS DESCRIBING THE SUPERNODAL STRUCTURE.
C       CACHSZ          -   CACHE SIZE (IN KBYTES).
C
C   OUTPUT PARAMETERS:
C       TMPSIZ          -   SIZE OF WORKING STORAGE REQUIRED BY BLKFCT.
C       SPLIT           -   SPLITTING OF SUPERNODES SO THAT THEY FIT
C                           INTO CACHE.
C
C***********************************************************************
C
      SUBROUTINE  BFINIT  ( NEQNS , NSUPER, XSUPER, SNODE , XLINDX,
     &                      LINDX , CACHSZ, TMPSIZ, SPLIT           )
C
C***********************************************************************
C
        INTEGER     CACHSZ, NEQNS , NSUPER, TMPSIZ
        INTEGER     XLINDX(*)       , XSUPER(*)
        INTEGER     LINDX (*)       , SNODE (*)   ,
     &              SPLIT(*)
C
C***********************************************************************
C
C       ---------------------------------------------------
C       DETERMINE FLOATING POINT WORKING SPACE REQUIREMENT.
C       ---------------------------------------------------
        CALL  FNTSIZ (  NSUPER, XSUPER, SNODE , XLINDX, LINDX ,
     &                  TMPSIZ                                  )
C
C       -------------------------------
C       PARTITION SUPERNODES FOR CACHE.
C       -------------------------------
        CALL  FNSPLT (  NEQNS , NSUPER, XSUPER, XLINDX, CACHSZ,
     &                  SPLIT                                   )
C
        RETURN
      END

max_stars_repo_name: Pangqiyuangh/SeIInv

max_issues_repo_name: Pangqiyuangh/SeIInv

max_forks_repo_name: Pangqiyuangh/SeIInv


text:
subroutine lev_SetExmod(a)
   use typre
   use Mod_LevelSet
   implicit none

   class(LevelSetProblem) :: a

   a%exmod = 'lev'
   a%namod = 'LevSet'
end subroutine

subroutine lev_SetNdofn(a)
   use typre
   use Mod_LevelSet
   implicit none

   class(LevelSetProblem) :: a

   a%ndofn = 1
end subroutine

subroutine lev_SetNdofbc(a)
   use typre
   use Mod_LevelSet
   implicit none

   class(LevelSetProblem) :: a
   integer(ip) :: ndime

   a%ndofbc = 1
   a%ndofbcstart = 0
end subroutine

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


text:
test4/test.cpp
test4/tb.cpp
test4/main.cpp


text:
PROGRAM snap_main

!-----------------------------------------------------------------------
!
!    SNAP - SN Application Proxy
!
!    Parallel programming model based on PARTISN
!
!
!    SNAP: SN (Discrete Ordinates) Application Proxy
!    Version 1.x (C13087)
!    LA-CC-13-016
!
!    This code is Unclassified, and contains no Unclassified Controlled
!    Nuclear Information
!
!    Copyright (c) 2013, Los Alamos National Security, LLC
!    All rights reserved.
!
!    Copyright 2013. Los Alamos National Security, LLC. This software
!    was produced under U.S. Government contract DE-AC52-06NA25396 for
!    Los Alamos National Laboratory (LANL), which is operated by Los
!    Alamos National Security, LLC for the U.S. Department of Energy.
!    The U.S. Government has rights to use, reproduce, and distribute
!    this software. NEITHER THE GOVERNMENT NOR LOS ALAMOS NATIONAL
!    SECURITY, LLC MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR ASSUMES
!    ANY LIABILITY FOR THE USE OF THIS SOFTWARE. If software is
!    modified to produce derivative works, such modified software should
!    be clearly marked, so as not to confuse it with the version
!    available from LANL.
!
!    Additionally, redistribution and use in source and binary forms,
!    with or without modification, are permitted provided that the
!    following conditions are met:
!    --Redistributions of source code must retain the above copyright
!      notice, this list of conditions and the following disclaimer.
!    --Redistributions in binary form must reproduce the above copyright
!      notice, this list of conditions and the following disclaimer in
!      the documentation and/or other materials provided with the
!      distribution.
!    --Neither the name of Los Alamos National Security, LLC, Los Alamos
!      National Laboratory, LANL, the U.S. Government, nor the names of
!      its contributors may be used to endorse or promote products
!      derived from this software without specific prior written
!      permission.
!
!    THIS SOFTWARE IS PROVIDED BY LOS ALAMOS NATIONAL SECURITY, LLC AND
!    CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
!    INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
!    MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
!    DISCLAIMED. IN NO EVENT SHALL LOS ALAMOS NATIONAL SECURITY, LLC OR
!    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
!    SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
!    LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
!    USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
!    ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
!    OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
!    OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
!    SUCH DAMAGE.
!
!-----------------------------------------------------------------------

  USE global_module, ONLY: i_knd, r_knd, ifile, ofile, iunit, ounit, allocator

  USE utils_module

  USE version_module, ONLY: version_print

  USE input_module, ONLY: input_read

  USE setup_module, ONLY: setup

  USE output_module, ONLY: output

  USE time_module, ONLY: tsnap, wtime, tparset, time_summ, tgrind

  USE dealloc_module

  USE plib_module, ONLY: pinit, iproc, root, comm_snap, bcast,         &
    pcomm_set, pinit_omp

  USE control_module, ONLY: otrdone, swp_typ

  use geom_module, only: nx, ny, nz
  use data_module, only: ng
  use sn_module, only: nang

  use iso_c_binding
  use sicm_f90

  IMPLICIT NONE
!_______________________________________________________________________
!
! Local variables
!_______________________________________________________________________

  CHARACTER(LEN=1) :: star='*'

  CHARACTER(LEN=64) :: error

  INTEGER(i_knd) :: ierr, i

  REAL(r_knd) :: t1, t2, t3, t4, t5
  integer(8) :: start, finish, clock_rate
  integer(c_int) :: pin_res

  type(sf_device_list) :: devices

  call system_clock(start, clock_rate)

  ! Our SICM initialization is done in three steps:
  ! sf_init looks for the memory devices that are available on the
  ! machine, and populates the devices list
  ! sf_get_device let us actually fetch an allocator, which is stored
  ! globally
  ! since this version was set up primarily for sanity-checking and
  ! timings, the device is selected literally
  ! device 2 will typically be a huge page allocator on numa node 0,
  ! but see examples/basic.c for an example of figuring out what all the
  ! memory devices actually are
  ! (note, some of the functions have not yet been wrapped for FORTRAN,
  ! so replicating that here may require writing some more wrappers)
  ! finally, sf_pin ensure the process only runs local to the memory
  call sf_init(devices)
  allocator = sf_get_device(devices, 1)
  pin_res = sf_pin(allocator)

!_______________________________________________________________________
!
! Perform calls that set up the parallel environment in MPI and
! OpenMP. Also starts the timer. Update parallel setup time.
!_______________________________________________________________________

  ierr = 0
  error = ' '

  CALL pinit ( t1 )

  CALL wtime ( t2 )
  tparset = tparset + t2 - t1
!_______________________________________________________________________
!
! Read the command line arguments to get i/o file names. Open the two
! files.
!_______________________________________________________________________

  CALL cmdarg ( ierr, error )
  CALL bcast ( ierr, comm_snap, root )
  IF ( ierr /= 0 ) THEN
    CALL print_error ( 0, error )
    CALL stop_run ( 0, 0, 0, 0 )
  END IF

  CALL open_file ( iunit, ifile, 'OLD', 'READ', ierr, error )
  CALL bcast ( ierr, comm_snap, root )
  IF ( ierr /= 0 ) THEN
    CALL print_error ( 0, error )
    CALL stop_run ( 0, 0, 0, 0 )
  END IF

  CALL open_file ( ounit, ofile, 'REPLACE', 'WRITE', ierr, error )
  CALL bcast ( ierr, comm_snap, root )
  IF ( ierr /= 0 ) THEN
    CALL print_error ( 0, error )
    CALL stop_run ( 0, 0, 0, 0 )
  END IF
!_______________________________________________________________________
!
! Write code version and execution time to output.
!_______________________________________________________________________

  IF ( iproc == root ) CALL version_print
!_______________________________________________________________________
!
! Read input
!_______________________________________________________________________

  CALL input_read

  CALL close_file ( iunit, ierr, error )
  CALL bcast ( ierr, comm_snap, root )
  IF ( ierr /= 0 ) THEN
    CALL print_error ( ounit, error )
    CALL stop_run ( 0, 0, 0, 0 )
  END IF
!_______________________________________________________________________
!
! Get nthreads for each proc. Print the warning about resetting nthreads
! if necessary. Don't stop run. Set up the SDD MPI topology.
!_______________________________________________________________________

  CALL wtime ( t3 )

  CALL pinit_omp ( ierr, error )
  IF ( ierr /= 0 ) CALL print_error ( 0, error )

  CALL pcomm_set

  CALL wtime ( t4 )
  tparset = tparset + t4 - t3
!_______________________________________________________________________
!
! Setup problem
!_______________________________________________________________________

  CALL setup
!_______________________________________________________________________
!
! Call for the problem solution
!_______________________________________________________________________

  CALL translv
!_______________________________________________________________________
!
! Output the results. Print the timing summary.
!_______________________________________________________________________

  CALL output
  IF ( iproc == root ) CALL time_summ
!_______________________________________________________________________
!
! Final cleanup: deallocate, close output file, end the program
!_______________________________________________________________________

  CALL dealloc_input ( 4 )
  CALL dealloc_solve ( swp_typ, 2 )

  CALL wtime ( t5 )
  tsnap = t5 - t1

  !IF ( iproc == root ) THEN
    !WRITE( ounit, 501 ) tsnap
    !WRITE( ounit, 502 ) tgrind, ( star, i = 1, 80 )
  !END IF

  CALL close_file ( ounit, ierr, error )
  CALL bcast ( ierr, comm_snap, root )
  call system_clock(finish)
  write (6,*) "time:", 1000 * (finish - start) / clock_rate, "ms"
  write (6,*) "memory:", 16_8 * int(nx, c_size_t) * int(ny, c_size_t) * &
    int(nz, c_size_t) * int(nang, c_size_t) * int(ng, c_size_t) * int(r_knd, c_size_t), "B"
  IF ( ierr /= 0 ) THEN
    CALL print_error ( 0, error )
    CALL stop_run ( 1, 0, 0, 0 )
  END IF

  IF ( otrdone ) THEN
    CALL stop_run ( 1, 0, 0, 1 )
  ELSE
    CALL stop_run ( 1, 0, 0, 2 )
  END IF
!_______________________________________________________________________

  501 FORMAT( 2X, 'Total Execution time', T41, ES11.4, / )
  502 FORMAT( 2X, 'Grind Time (nanoseconds)', 8X, ES11.4, /, /, 80A )
!_______________________________________________________________________
!_______________________________________________________________________

END PROGRAM snap_main

max_stars_repo_name: gvallee/SICM

max_issues_repo_name: gvallee/SICM

max_forks_repo_name: gvallee/SICM


text:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!                                                                   !!
!!                   GNU General Public License                      !!
!!                                                                   !!
!! This file is part of the Flexible Modeling System (FMS).          !!
!!                                                                   !!
!! FMS is free software; you can redistribute it and/or modify       !!
!! it and are expected to follow the terms of the GNU General Public !!
!! License as published by the Free Software Foundation.             !!
!!                                                                   !!
!! FMS is distributed in the hope that it will be useful,            !!
!! but WITHOUT ANY WARRANTY; without even the implied warranty of    !!
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the     !!
!! GNU General Public License for more details.                      !!
!!                                                                   !!
!! You should have received a copy of the GNU General Public License !!
!! along with FMS; if not, write to:                                 !!
!!          Free Software Foundation, Inc.                           !!
!!          59 Temple Place, Suite 330                               !!
!!          Boston, MA  02111-1307  USA                              !!
!! or see:                                                           !!
!!          http://www.gnu.org/licenses/gpl.txt                      !!
!!                                                                   !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

module bgrid_sponge_mod

!-----------------------------------------------------------------------
!
!   (1) eddy damping of prognostic fields (t,q,u,v) at the top
!       level of the model (i.e., fields are damped with their
!       zonal mean values). THIS OPTION WILL NOT WORK CORRECTLY
!       WITH DOMAIN DECOMPOSITION ALONG THE X-AXIS.
!
!-----------------------------------------------------------------------

use types_mod, only : r8
 use bgrid_horiz_mod,       only: horiz_grid_type
 use bgrid_masks_mod,       only: grid_mask_type
 use bgrid_prog_var_mod,    only: prog_var_type
 use bgrid_change_grid_mod, only: mass_to_vel
 use bgrid_halo_mod,        only: update_halo, vel_flux_boundary, &
                                  TEMP, UWND, VWND, &
                                  NORTH, EAST, NOPOLE

 use         fms_mod, only: error_mesg, FATAL, write_version_number

! use mpp_domains_mod, only: domain2d, mpp_redistribute, &
!                            mpp_define_domains,         &
!                            mpp_get_global_domain,      &
!                            mpp_get_data_domain,        &
!                            CYCLIC_GLOBAL_DOMAIN

 implicit none
 private

 public   sponge_driver, sponge_init

!-----------------------------------------------------------------------

 real(r8), parameter ::  daypsec=1./86400.
 logical :: do_init=.true.

 logical :: do_topsponge
 integer :: nlev_sponge
 real(r8)    :: dfactr   ! coeff for damping eddies at the top level

 character(len=128) :: version='$Revision$'
 character(len=128) :: tagname='$Id$'
 logical :: do_log = .true.

 real(r8) :: small = .000001

 real(r8) :: coeff_vel, coeff_tmp, coeff_trs
 real(r8) :: xaxis_wt
 integer :: numlev

!--- module storage for computing exact/reproducible zonal means ---
  type zsum_type
     !!!type(domain2d) :: Domain
     integer        :: isize, jsize
  end type
  type(zsum_type) :: Zdomain_tmp, Zdomain_vel
  logical :: first = .true.
!-----------------------------------------------------------------------

contains

!#######################################################################

 subroutine sponge_driver ( Hgrid, nplev, dt, dpde, Var, Var_dt )

!-----------------------------------------------------------------------
 type (horiz_grid_type), intent(inout) :: Hgrid
 integer,                intent(in)    :: nplev
 real(r8),                   intent(in)    :: dt
 real(r8),                   intent(in)    :: dpde(Hgrid%ilb:,Hgrid%jlb:,:)
 type  (prog_var_type),  intent(in)    :: Var
 type  (prog_var_type),  intent(inout) :: Var_dt
!-----------------------------------------------------------------------

 real(r8), dimension(Hgrid%ilb:Hgrid%iub,Hgrid%jlb:Hgrid%jub,numlev) :: dpxy
 integer :: n, np, nt

!-----------------------------------------------------------------------

   if (do_init) call error_mesg ('bgrid_sponge_mod',  &
                                 'initialization not called', FATAL)

   if ( numlev == 0 ) return

!  do not allow decomposition along x-axis
!  if ( Hgrid%decompx ) call error_mesg ('bgrid_sponge_mod', &
!    'this version of sponge does not allow x-axis decomposition.', FATAL)

!-----------------------------------------------------------------------

   n  = numlev
   np = nplev

!---- temperature and tracers -----

   if ( coeff_tmp > small .or. coeff_trs > small ) then
       if (first) then
          Zdomain_tmp%isize = Hgrid%Tmp%ie - Hgrid%Tmp%is + 1
          Zdomain_tmp%jsize = Hgrid%Tmp%je - Hgrid%Tmp%js + 1
       endif
       !!!if (first) call zsum_init ( Hgrid%Tmp%Domain, Zdomain_tmp )
       nt = Var_dt%ntrace
       call local_filter_mass ( Hgrid, np, coeff_tmp, coeff_trs, dt,       &
                dpde(:,:,1:n), Var   %t(:,:,1:n), Var   %r(:,:,1:n,1:nt),  &
                               Var_dt%t(:,:,1:n), Var_dt%r(:,:,1:n,1:nt) )
   endif

!---- momentum components -----

   if ( coeff_vel > small ) then
        if (first) then
          Zdomain_vel%isize = Hgrid%Vel%ie - Hgrid%Vel%is + 1
          Zdomain_vel%jsize = Hgrid%Vel%je - Hgrid%Vel%js + 1
        endif
        !!!if (first) call zsum_init ( Hgrid%Vel%Domain, Zdomain_vel )
      ! compute pressure weights at velocity points
        dpxy(:,:,:) = dpde(:,:,1:n)
        if (np < n) then
          call mass_to_vel (Hgrid, dpxy(:,:,np+1:n), dpxy(:,:,np+1:n))
          call update_halo (Hgrid, UWND, dpxy(:,:,np+1:n), EAST+NORTH+NOPOLE)
        endif

        call local_filter_vel ( Hgrid, np, coeff_vel, dt, dpxy,       &
                                Var   %u(:,:,1:n), Var   %v(:,:,1:n), &
                                Var_dt%u(:,:,1:n), Var_dt%v(:,:,1:n)  )
   endif

   if (first) first = .false.

!-----------------------------------------------------------------------

 end subroutine sponge_driver

!#######################################################################

 subroutine sponge_init ( damp_vel, damp_tmp, damp_trs, nlev, xwt )

   real(r8),    intent(in), optional :: damp_vel, damp_tmp, damp_trs, xwt
   integer, intent(in), optional :: nlev

!
!   damp_vel, damp_tmp, damp_trs =
!         normalized (0.,1.) damping coefficients for
!         momentum, temperature, and tracers
!         [real(r8), default = 0.]
!
!   nlev   = number of levels at the top of the model where
!            damping is performed
!               [integer, default = 0]
!
!   xwt    = weight (0.,1.) applied to the x-axis
!            the y-axis will have a weight = 1.-xwt
!               [real(r8), default = .50]
!

!-----------------------------------------------------------------------

   if (do_log) then
      call write_version_number (version,tagname)
      do_log = .false.
   endif

! set values for optional arguments

   coeff_vel = 0.; if (present(damp_vel)) coeff_vel = min(max(damp_vel,0.0_r8),1.0_r8)
   coeff_tmp = 0.; if (present(damp_tmp)) coeff_tmp = min(max(damp_tmp,0.0_r8),1.0_r8)
   coeff_trs = 0.; if (present(damp_trs)) coeff_trs = min(max(damp_trs,0.0_r8),1.0_r8)

   numlev   = 0;   if (present(nlev)) numlev   = max(nlev,0)
   xaxis_wt = .50; if (present(xwt))  xaxis_wt = min(max(xwt,0.0_r8),1.0_r8)

!  do not allow more than one sponge layer
   if (numlev > 1) call error_mesg ('bgrid_sponge_mod',  &
                                    'numlev > 1 ', FATAL)

   do_init=.false.

!-----------------------------------------------------------------------

 end subroutine sponge_init

!#######################################################################

 subroutine local_filter_mass ( Hgrid, nplev, coeff_tmp, coeff_trs, dt, &
                                dp, t, tr, tdt, trdt )

   type (horiz_grid_type), intent(inout) :: Hgrid
   integer, intent(in)                   :: nplev
   real(r8),    intent(in)                   :: coeff_tmp, coeff_trs, dt
   real(r8),    intent(in),    dimension(Hgrid%ilb:,Hgrid%jlb:,:) :: dp, t
   real(r8),    intent(inout), dimension(Hgrid%ilb:,Hgrid%jlb:,:) :: tdt
   real(r8),    intent(in),    dimension(Hgrid%ilb:,Hgrid%jlb:,:,:) :: tr
   real(r8),    intent(inout), dimension(Hgrid%ilb:,Hgrid%jlb:,:,:) :: trdt

   real(r8), dimension(Hgrid%ilb:Hgrid%iub,Hgrid%jlb:Hgrid%jub,size(t,3)) :: tmp, akew, akns, akdp, ztmp
   real(r8), dimension(Hgrid%ilb:Hgrid%iub,Hgrid%jlb:Hgrid%jub) :: akew2, akns2, tew, tns
   integer :: i, j, k, is, ie, js, je, jsd, jed, n
   real(r8)    :: xwt, ywt

   is  = Hgrid%Tmp%is;   ie  = Hgrid%Tmp%ie
   js  = Hgrid%Tmp%js;   je  = Hgrid%Tmp%je
   jsd = Hgrid%Tmp%jsd;  jed = Hgrid%Tmp%jed

   xwt = 0.125 * xaxis_wt
   ywt = 0.125 * (1.-xaxis_wt)

 ! 2D geometric constants
   do j = js-1, je
   do i = is-1, ie
      akew2(i,j) = xwt * (Hgrid%Tmp%area(i,j)+Hgrid%Tmp%area(i+1,j))
      akns2(i,j) = ywt * (Hgrid%Tmp%area(i,j)+Hgrid%Tmp%area(i,j+1))
   enddo
   enddo

 ! 3D mass weighted constants
   do k = 1, size(t,3)
      akew(is-1:ie,js-1:je,k) = akew2(is-1:ie,js-1:je)
      akns(is-1:ie,js-1:je,k) = akns2(is-1:ie,js-1:je)
      akdp(:,:,k) = Hgrid%Tmp%rarea(:,:)/dt
   enddo
   do k = nplev+1, size(t,3)
      do j = js-1, je
      do i = is-1, ie
         akew(i,j,k) = akew(i,j,k) * (dp(i,j,k)+dp(i+1,j,k))
         akns(i,j,k) = akns(i,j,k) * (dp(i,j,k)+dp(i,j+1,k))
         akdp(i,j,k) = akdp(i,j,k) / (2.*dp(i,j,k))
      enddo
      enddo
   enddo

 ! temperature

   if ( coeff_tmp > small ) then
      tmp = t + dt*tdt
    !!!!! remove zonal mean !!!!!
      call zsum (Zdomain_tmp, tmp(:,jsd:jed,:), ztmp(:,jsd:jed,:))
      tmp(is:ie,js:je,:) = tmp(is:ie,js:je,:) -ztmp(is:ie,js:je,:)/(1.0 * Hgrid%nlon)
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      call update_halo ( Hgrid, TEMP, tmp )
      do k = 1, size(t,3)
         do j = js-1, je
         do i = is-1, ie
            tew(i,j) = akew(i,j,k) * (tmp(i+1,j,k)-tmp(i,j,k))
            tns(i,j) = akns(i,j,k) * (tmp(i,j+1,k)-tmp(i,j,k))
         enddo
         enddo
         do j = js, je
         do i = is, ie
            tdt(i,j,k) = tdt(i,j,k) + &
                  coeff_tmp*(tew(i,j)-tew(i-1,j)+tns(i,j)-tns(i,j-1))*akdp(i,j,k)
         enddo
         enddo
      enddo
   endif

 ! tracers

   if ( coeff_trs > small ) then
      do n = 1, size(tr,4)
         tmp = tr(:,:,:,n) + dt*trdt(:,:,:,n)
       !!!!! remove zonal mean !!!!!
         call zsum (Zdomain_tmp, tmp(:,jsd:jed,:), ztmp(:,jsd:jed,:))
         tmp(is:ie,js:je,:) = tmp(is:ie,js:je,:) -ztmp(is:ie,js:je,:)/(1.0 * Hgrid%nlon)
       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         call update_halo ( Hgrid, TEMP, tmp )
         do k = 1, size(tr,3)
            do j = js-1, je
            do i = is-1, ie
               tew(i,j) = akew(i,j,k) * (tmp(i+1,j,k)-tmp(i,j,k))
               tns(i,j) = akns(i,j,k) * (tmp(i,j+1,k)-tmp(i,j,k))
            enddo
            enddo
            do j = js, je
            do i = is, ie
               trdt(i,j,k,n) = trdt(i,j,k,n) +  &
                    coeff_trs*(tew(i,j)-tew(i-1,j)+tns(i,j)-tns(i,j-1))*akdp(i,j,k)
            enddo
            enddo
         enddo
      enddo
   endif

 end subroutine local_filter_mass

!#######################################################################

 subroutine local_filter_vel ( Hgrid, nplev, coeff, dt, dp, u, v, udt, vdt )

   type (horiz_grid_type), intent(inout) :: Hgrid
   integer, intent(in)                   :: nplev
   real(r8),    intent(in)                   :: coeff, dt
   real(r8),    intent(in),    dimension(Hgrid%ilb:,Hgrid%jlb:,:) :: dp, u, v
   real(r8),    intent(inout), dimension(Hgrid%ilb:,Hgrid%jlb:,:) :: udt, vdt

   real(r8), dimension(Hgrid%ilb:Hgrid%iub,Hgrid%jlb:Hgrid%jub,size(u,3)) :: akew, akns, akdp, uu, vv, ztmp
   real(r8), dimension(Hgrid%ilb:Hgrid%iub,Hgrid%jlb:Hgrid%jub) ::  akew2, akns2, uew, vew, uns, vns
   integer :: i, j, k, is, ie, js, je, jsd, jed
   real(r8)    :: xwt, ywt

   is  = Hgrid%Tmp%is;   ie  = Hgrid%Tmp%ie
   js  = Hgrid%Vel%js;   je  = Hgrid%Vel%je
   jsd = Hgrid%Vel%jsd;  jed = Hgrid%Vel%jed

   xwt = 0.125 * coeff * xaxis_wt
   ywt = 0.125 * coeff * (1.-xaxis_wt)

 ! 2D geometric constants
   do j = js, je+1
   do i = is, ie+1
      akew2(i,j) = xwt * (Hgrid%Vel%area(i-1,j)+Hgrid%Vel%area(i,j))
      akns2(i,j) = ywt * (Hgrid%Vel%area(i,j-1)+Hgrid%Vel%area(i,j))
   enddo
   enddo

 ! 3D mass weighted constants
   do k = 1, size(u,3)
      akew(is:ie+1,js:je+1,k) = akew2(is:ie+1,js:je+1)
      akns(is:ie+1,js:je+1,k) = akns2(is:ie+1,js:je+1)
      akdp(:,:,k) = Hgrid%Vel%rarea(:,:)/dt
   enddo
   do k = nplev+1, size(u,3)
      do j = js, je+1
      do i = is, ie+1
         akew(i,j,k) = akew(i,j,k) * (dp(i,j,k)+dp(i-1,j,k))
         akns(i,j,k) = akns(i,j,k) * (dp(i,j,k)+dp(i,j-1,k))
         akdp(i,j,k) = akdp(i,j,k) / (2.*dp(i,j,k))
      enddo
      enddo
   enddo

   uu = u + dt*udt
   vv = v + dt*vdt
 !!!!! remove zonal mean from u comp !!!!!
   call zsum (Zdomain_vel, uu(:,jsd:jed,:), ztmp(:,jsd:jed,:))
   uu(is:ie,js:je,:) = uu(is:ie,js:je,:) -ztmp(is:ie,js:je,:)/(1.0 * Hgrid%nlon)
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   call update_halo ( Hgrid, UWND, uu )
   call update_halo ( Hgrid, VWND, vv )

   do k = 1, size(u,3)
      do j = js, je+1
      do i = is, ie+1
         uew(i,j) = akew(i,j,k) * (uu(i,j,k)-uu(i-1,j,k))
         vew(i,j) = akew(i,j,k) * (vv(i,j,k)-vv(i-1,j,k))
         uns(i,j) = akns(i,j,k) * (uu(i,j,k)-uu(i,j-1,k))
         vns(i,j) = akns(i,j,k) * (vv(i,j,k)-vv(i,j-1,k))
      enddo
      enddo
     !---- remove meridional averaging at poles ----
      call vel_flux_boundary (Hgrid, uns)
      call vel_flux_boundary (Hgrid, vns)

      do j = js, je
      do i = is, ie
         udt(i,j,k) = udt(i,j,k) + (uew(i+1,j)-uew(i,j)+uns(i,j+1)-uns(i,j))*akdp(i,j,k)
         vdt(i,j,k) = vdt(i,j,k) + (vew(i+1,j)-vew(i,j)+vns(i,j+1)-vns(i,j))*akdp(i,j,k)
      enddo
      enddo
   enddo

 end subroutine local_filter_vel

!#######################################################################
! initializes domain2d type for summation in zonal direction

! subroutine zsum_init ( Domain, Zonal )
! type(domain2d),  intent(in)   :: Domain
! type(zsum_type), intent(out)  :: Zonal
! integer :: isg, ieg, jsg, jeg, layout(2)

! create new domian2d type with 1d zonal decomposition
!  call mpp_get_global_domain  ( Domain, isg, ieg, jsg, jeg )
!  layout = (/1,mpp_npes()/)

! error check
!  if (jeg-jsg+1 < mpp_npes()) call error_mesg ('bgrid_sponge_mod', 'number &
!     &of global latitude rows is less than the number of processors', FATAL)
!
!  call mpp_define_domains ( (/isg,ieg,jsg,jeg/), layout, Zonal%Domain, &
!                            xflags=CYCLIC_GLOBAL_DOMAIN )
!  call mpp_get_data_domain ( Zonal%Domain, xsize=Zonal%isize, ysize=Zonal%jsize )

! end subroutine zsum_init

!#######################################################################

 subroutine zsum (Zonal, array, zarray )
 type(zsum_type), intent(in)   :: Zonal
 real(r8),            intent(in)  ::   array(:,:,:)
 real(r8),            intent(out) ::  zarray(:,:,:)
 real(r8), dimension(Zonal%isize,Zonal%jsize,size(array,3)) :: zdat
 real(r8)    :: c
 integer :: j, k

!!!   call mpp_redistribute ( Domain, array, Zonal%Domain, zdat )
zdat(:, :, :) = array(2:Zonal%isize + 1, 2:Zonal%jsize + 1, :)

   do k = 1, size(array,3)
   do j = 1, Zonal%jsize
     c = sum(zdat(:,j,k))
     zdat(:,j,k) = c
   enddo
   enddo
!!!   call mpp_redistribute ( Zonal%Domain, zdat, Domain, zarray )
zarray(2:Zonal%isize +1, 2:Zonal%jsize +1, :) = zdat(:, :, :)

 end subroutine zsum

!#######################################################################

end module bgrid_sponge_mod

max_stars_repo_name: hkershaw-brown/feature-preprocess

max_issues_repo_name: hkershaw-brown/feature-preprocess

max_forks_repo_name: hkershaw-brown/feature-preprocess


text:
! Program :
! Author  : wansooha@gmail.com
! Date    :

        program test_cfg_read
        use cfgio_mod
        implicit none
        type(cfg_t):: cfg
        integer n1
        integer,allocatable:: n1a(:)
        real o1
        real,allocatable:: o1a(:)
        real(kind=8) d1
        real(kind=8),allocatable:: d1a(:)
        character(len=16):: label2
        character(len=16),allocatable:: label2a(:)
        logical use_swd
        logical,allocatable:: la(:)
        complex c1
        complex,allocatable:: c1a(:)
        complex(kind=8):: c2
        complex(kind=8),allocatable:: c2a(:)
        integer i,npar
        print*,'parse cfg'
        cfg=parse_cfg("test.cfg")
        print*,'n sect=',cfg%nsect()
        print*,'read i'
        call cfg%get("DEFAULTS","n1",n1)
        call cfg%get("DEFAULTS","n1a",n1a,npar)
        print*,'npar=',npar
        print*,'read f'
        call cfg%get("sec 01",'o1',o1)
        call cfg%get("sec 01",'o1a',o1a,npar)
        print*,'npar=',npar
        call cfg%get("sec 01",'d1',d1)
        call cfg%get("sec 01",'d1a',d1a,npar)
        print*,'npar=',npar
        print*,'read s'
        call cfg%get("sec 01",'label2',label2)
        call cfg%get("sec 01",'label2a',label2a,npar)
        print*,'npar=',npar
        print*,'read b'
        call cfg%get("sec 02",'use swd',use_swd)
        call cfg%get("sec 02",'la',la,npar)
        print*,'npar=',npar
        print*,'read c'
        call cfg%get("sec 03",'c1',c1)
        call cfg%get("sec 03",'c1a',c1a,npar)
        print*,'npar=',npar
        print*,'read z'
        call cfg%get("sec 03",'c2',c2)
        call cfg%get("sec 03",'c2a',c2a,npar)
        print*,'npar=',npar
        print*,'read end'

        print*,'n1=',n1
        print*,'n1a=',n1a
        print*,'o1=',o1
        print*,'o1a=',o1a
        print*,'d1=',d1
        print*,'d1a=',d1a
        print*,'label2=',trim(label2)
        do i=1,size(label2a)
            print*,'label2a=',i,label2a(i)
        enddo
        print*,'use swd=',use_swd
        print*,'la=',la
        print*,'c1',c1
        print*,'c1a',c1a
        print*,'c2',c2
        print*,'c2a',c2a
        call cfg%write()
        call cfg%finalize()

        end

max_stars_repo_name: baklanovp/cfgio

max_issues_repo_name: baklanovp/cfgio

max_forks_repo_name: baklanovp/cfgio


text:
program prgm
  implicit none
  call main
contains
  subroutine l(s)
    implicit none
    character (len = *), intent (in) :: s
    print*, s
    return
  end subroutine l

  subroutine main()
    implicit none
    call l("Fortan")
  end subroutine main
end program prgm

max_stars_repo_name: justletterh/el

max_issues_repo_name: justletterh/el

max_forks_repo_name: justletterh/el