snap.cpp

#include <iostream>
#include <string>
#include <vector>
#include <Kokkos_Core.hpp>

template <typename T>
Kokkos::View<T*> from_std_vector(const std::vector<T> &v)
{
    Kokkos::View<T*> r("r", v.size());
    for (size_t i=0; i < v.size(); i++) {
        r(i) = v[i];
    }
    return r;
}

int abs(const Kokkos::View<const int*> &x)
{
    int abs;
    return abs;
}

bool allocated(const Kokkos::View<const int*> &x)
{
    bool allocated;
    return allocated;
}

int floor(const Kokkos::View<const float*> &x, int kind)
{
    int floor;
    return floor;
}

int lbound(const Kokkos::View<const int*> &x, int dim)
{
    int lbound;
    return lbound;
}

int max(int a, int b)
{
    int max;
    return max;
}

int maxval(const Kokkos::View<const int*> &x)
{
    int maxval;
    return maxval;
}

int min(int a, int b)
{
    int min;
    return min;
}

int minval(const Kokkos::View<const int*> &x)
{
    int minval;
    return minval;
}

float real(const Kokkos::View<const int*> &x, int kind)
{
    float real;
    return real;
}

int sum(const Kokkos::View<const int*> &x)
{
    int sum;
    return sum;
}

float tiny(const Kokkos::View<const int*> &x)
{
    float tiny;
    return tiny;
}

int ubound(const Kokkos::View<const int*> &x, int dim)
{
    int ubound;
    return ubound;
}

float dabs(float x)
{
    float r;
    if (x >= 0) {
        r = x;
    } else {
        r = -x;
    };
    return r;
}

float daimag(std::complex x)
{
    float r;
    std::cerr << "ERROR STOP" << std::endl;
    exit(1);
    return r;
}

int dfloor(float x)
{
    int r;
    if (x >= 0) {
        r = (int)(x);
    } else {
        r = (int)(x - 1);
    };
    return r;
}

float dmodulo(float x, float y)
{
    float r;
    r = x - dfloor(x/y)*y;
    return r;
}

float dsqrt(float x)
{
    float r;
    if (x >= 0) {
        r = std::pow(x, 1.000000/2);
    } else {
        std::cerr << "ERROR STOP" << std::endl;
        exit(1);
    };
    return r;
}

int iabs(int x)
{
    int r;
    if (x >= 0) {
        r = x;
    } else {
        r = -x;
    };
    return r;
}

int imodulo(int x, int y)
{
    int r;
    r = x - sfloor(x/y)*y;
    return r;
}

float sabs(float x)
{
    float r;
    if (x >= 0) {
        r = x;
    } else {
        r = -x;
    };
    return r;
}

float saimag(std::complex x)
{
    float r;
    std::cerr << "ERROR STOP" << std::endl;
    exit(1);
    return r;
}

int sfloor(float x)
{
    int r;
    if (x >= 0) {
        r = (int)(x);
    } else {
        r = (int)(x - 1);
    };
    return r;
}

float smodulo(float x, float y)
{
    float r;
    r = x - sfloor(x/y)*y;
    return r;
}

float ssqrt(float x)
{
    float r;
    if (x >= 0) {
        r = std::pow(x, 1.000000/2);
    } else {
        std::cerr << "ERROR STOP" << std::endl;
        exit(1);
    };
    return r;
}

float cabs(std::complex x)
{
    float r;
    r = ssqrt(std::pow((double)(x), 2) + std::pow(saimag(x), 2));
    return r;
}

std::complex cacos(std::complex x)
{
    std::complex r;
    r = c_cacos(x);
    return r;
}

std::complex cacosh(std::complex x)
{
    std::complex r;
    r = c_cacosh(x);
    return r;
}

std::complex casin(std::complex x)
{
    std::complex r;
    r = c_casin(x);
    return r;
}

std::complex casinh(std::complex x)
{
    std::complex r;
    r = c_casinh(x);
    return r;
}

std::complex catan(std::complex x)
{
    std::complex r;
    r = c_catan(x);
    return r;
}

std::complex catanh(std::complex x)
{
    std::complex r;
    r = c_catanh(x);
    return r;
}

std::complex ccos(std::complex x)
{
    std::complex r;
    r = c_ccos(x);
    return r;
}

std::complex ccosh(std::complex x)
{
    std::complex r;
    r = c_ccosh(x);
    return r;
}

std::complex cexp(std::complex x)
{
    std::complex r;
    r = c_cexp(x);
    return r;
}

std::complex clog(std::complex x)
{
    std::complex r;
    r = c_clog(x);
    return r;
}

std::complex csin(std::complex x)
{
    std::complex r;
    r = c_csin(x);
    return r;
}

std::complex csinh(std::complex x)
{
    std::complex r;
    r = c_csinh(x);
    return r;
}

std::complex csqrt(std::complex x)
{
    std::complex r;
    r = c_csqrt(x);
    return r;
}

std::complex ctan(std::complex x)
{
    std::complex r;
    r = c_ctan(x);
    return r;
}

std::complex ctanh(std::complex x)
{
    std::complex r;
    r = c_ctanh(x);
    return r;
}

float dabs(float x)
{
    float r;
    if (x >= 0) {
        r = x;
    } else {
        r = -x;
    };
    return r;
}

float dacos(float x)
{
    float r;
    r = c_dacos(x);
    return r;
}

float dacosh(float x)
{
    float r;
    r = c_dacosh(x);
    return r;
}

float dasin(float x)
{
    float r;
    r = c_dasin(x);
    return r;
}

float dasinh(float x)
{
    float r;
    r = c_dasinh(x);
    return r;
}

float datan(float x)
{
    float r;
    r = c_datan(x);
    return r;
}

float datan2(float y, float x)
{
    float r;
    r = c_datan2(y, x);
    return r;
}

float datanh(float x)
{
    float r;
    r = c_datanh(x);
    return r;
}

float dcos(float x)
{
    float r;
    r = c_dcos(x);
    return r;
}

float dcosh(float x)
{
    float r;
    r = c_dcosh(x);
    return r;
}

float derf(float x)
{
    float r;
    r = c_derf(x);
    return r;
}

float dexp(float x)
{
    float r;
    r = c_dexp(x);
    return r;
}

float dlog(float x)
{
    float r;
    r = c_dlog(x);
    return r;
}

float dsin(float x)
{
    float r;
    r = c_dsin(x);
    return r;
}

float dsinh(float x)
{
    float r;
    r = c_dsinh(x);
    return r;
}

float dsqrt(float x)
{
    float r;
    if (x >= 0) {
        r = std::pow(x, 1.000000/2);
    } else {
        std::cerr << "ERROR STOP" << std::endl;
        exit(1);
    };
    return r;
}

float dtan(float x)
{
    float r;
    r = c_dtan(x);
    return r;
}

float dtanh(float x)
{
    float r;
    r = c_dtanh(x);
    return r;
}

float sabs(float x)
{
    float r;
    if (x >= 0) {
        r = x;
    } else {
        r = -x;
    };
    return r;
}

float sacos(float x)
{
    float r;
    r = c_sacos(x);
    return r;
}

float sacosh(float x)
{
    float r;
    r = c_sacosh(x);
    return r;
}

float sasin(float x)
{
    float r;
    r = c_sasin(x);
    return r;
}

float sasinh(float x)
{
    float r;
    r = c_sasinh(x);
    return r;
}

float satan(float x)
{
    float r;
    r = c_satan(x);
    return r;
}

float satan2(float y, float x)
{
    float r;
    r = c_satan2(y, x);
    return r;
}

float satanh(float x)
{
    float r;
    r = c_satanh(x);
    return r;
}

float scos(float x)
{
    float r;
    r = c_scos(x);
    return r;
}

float scosh(float x)
{
    float r;
    r = c_scosh(x);
    return r;
}

float serf(float x)
{
    float r;
    r = c_serf(x);
    return r;
}

float sexp(float x)
{
    float r;
    r = c_sexp(x);
    return r;
}

float slog(float x)
{
    float r;
    r = c_slog(x);
    return r;
}

float ssin(float x)
{
    float r;
    r = c_ssin(x);
    return r;
}

float ssinh(float x)
{
    float r;
    r = c_ssinh(x);
    return r;
}

float ssqrt(float x)
{
    float r;
    if (x >= 0) {
        r = std::pow(x, 1.000000/2);
    } else {
        std::cerr << "ERROR STOP" << std::endl;
        exit(1);
    };
    return r;
}

float stan(float x)
{
    float r;
    r = c_stan(x);
    return r;
}

float stanh(float x)
{
    float r;
    r = c_stanh(x);
    return r;
}

float zabs(std::complex x)
{
    float r;
    r = dsqrt(std::pow((double)(x), 2) + std::pow(daimag(x), 2));
    return r;
}

std::complex zacos(std::complex x)
{
    std::complex r;
    r = c_zacos(x);
    return r;
}

std::complex zacosh(std::complex x)
{
    std::complex r;
    r = c_zacosh(x);
    return r;
}

std::complex zasin(std::complex x)
{
    std::complex r;
    r = c_zasin(x);
    return r;
}

std::complex zasinh(std::complex x)
{
    std::complex r;
    r = c_zasinh(x);
    return r;
}

std::complex zatan(std::complex x)
{
    std::complex r;
    r = c_zatan(x);
    return r;
}

std::complex zatanh(std::complex x)
{
    std::complex r;
    r = c_zatanh(x);
    return r;
}

std::complex zcos(std::complex x)
{
    std::complex r;
    r = c_zcos(x);
    return r;
}

std::complex zcosh(std::complex x)
{
    std::complex r;
    r = c_zcosh(x);
    return r;
}

std::complex zexp(std::complex x)
{
    std::complex r;
    r = c_zexp(x);
    return r;
}

std::complex zlog(std::complex x)
{
    std::complex r;
    r = c_zlog(x);
    return r;
}

std::complex zsin(std::complex x)
{
    std::complex r;
    r = c_zsin(x);
    return r;
}

std::complex zsinh(std::complex x)
{
    std::complex r;
    r = c_zsinh(x);
    return r;
}

std::complex zsqrt(std::complex x)
{
    std::complex r;
    r = c_zsqrt(x);
    return r;
}

std::complex ztan(std::complex x)
{
    std::complex r;
    r = c_ztan(x);
    return r;
}

std::complex ztanh(std::complex x)
{
    std::complex r;
    r = c_ztanh(x);
    return r;
}

float abs(float x)
{
    float r;
    if (x >= 0) {
        r = x;
    } else {
        r = -x;
    };
    return r;
}

float dsin(float x)
{
    int n;
    float r;
    float y;
    y = modulo(x, 2*pi);
    y = min(y, pi - y);
    y = max(y, (-pi) - y);
    y = min(y, pi - y);
    r = kernel_dsin(y);
    return r;
}

int floor(float x)
{
    int r;
    if (x >= 0) {
        r = (int)(x);
    } else {
        r = (int)(x - 1);
    };
    return r;
}

float kernel_dsin(float x)
{
    float res;
    float s1;
    float s2;
    float s3;
    float s4;
    float s5;
    float s6;
    float s7;
    float s8;
    float z;
    z = x*x;
    res = x*(s1 + z*(s2 + z*(s3 + z*(s4 + z*(s5 + z*(s6 + z*(s7 + z*s8)))))));
    return res;
}

float max(float x, float y)
{
    float r;
    if (x > y) {
        r = x;
    } else {
        r = y;
    };
    return r;
}

float min(float x, float y)
{
    float r;
    if (x < y) {
        r = x;
    } else {
        r = y;
    };
    return r;
}

float modulo(float x, float y)
{
    float r;
    r = x - floor(x/y)*y;
    return r;
}

float ssin(float x)
{
    float r;
    float tmp;
    float x2;
    float y;
    x2 = x;
    tmp = dsin(x2);
    r = tmp;
    return r;
}

int len_trim(std::string string)
{
    int r;
    r = (string).size();
    if (r == 0) {
        return;
    };
    while (string[r-1] == " ") {
        r = r - 1;
        if (r == 0) {
            break;
        };
    };
    return r;
}

std::string trim(std::string x)
{
    int i;
    std::string r;
    for (i=1; i<=(r).size(); i++) {
        r[i-1] = x[i-1];
    };
    return r;
}

void control_allocate(int ng, int &ndpwds, int &ierr)
{
    ierr = 0;
    // FIXME: allocate(dfmxi, inrdone, );
    if (ierr != 0) {
        return;
    };
    dfmxi = -one;
    inrdone = false;
    dfmxo = -one;
    otrdone = false;
    update_ptr = false;
    ndpwds = ndpwds + dfmxi.extent(0) + inrdone.extent(0);
}

void control_deallocate()
{
    // FIXME: deallocate(dfmxi, inrdone, );
}

void data_allocate(int nx, int ny, int nz, int nmom, int nang, int noct, int timedep, int &ndpwds, int &istat)
{
    if (mat_opt > 0) {
        nmat = 2;
    };
    istat = 0;
    if (timedep == 1) {
        // FIXME: allocate(v, );
    } else {
        // FIXME: allocate(v, );
    };
    if (istat != 0) {
        return;
    };
    v = zero;
    // FIXME: allocate(mat, );
    if (istat != 0) {
        return;
    };
    mat = 1;
    if (src_opt < 3) {
        // FIXME: allocate(qi, qim, );
        if (istat != 0) {
            return;
        };
        qi = zero;
    } else {
        // FIXME: allocate(qi, qim, );
        if (istat != 0) {
            return;
        };
        qi = zero;
        qim = zero;
    };
    // FIXME: allocate(sigt, siga, sigs, slgg, );
    if (istat != 0) {
        return;
    };
    sigt = zero;
    siga = zero;
    sigs = zero;
    slgg = zero;
    // FIXME: allocate(vdelt, );
    if (istat != 0) {
        return;
    };
    vdelt = zero;
    ndpwds = ndpwds + v.extent(0) + mat.extent(0) + qi.extent(0) + qim.extent(0) + sigt.extent(0) + siga.extent(0) + sigs.extent(0) + slgg.extent(0) + vdelt.extent(0);
}

void data_deallocate()
{
    // FIXME: deallocate(v, );
    // FIXME: deallocate(mat, );
    // FIXME: deallocate(qi, qim, );
    // FIXME: deallocate(sigt, siga, sigs, );
    // FIXME: deallocate(slgg, );
    // FIXME: deallocate(vdelt, );
}

void sn_allocate(int ndimen, int &ndpwds, int &istat)
{
    cmom = nmom;
    noct = 2;
    // FIXME: allocate(mu, w, wmu, eta, weta, xi, wxi, );
    if (istat > 0) {
        return;
    };
    w = zero;
    mu = zero;
    wmu = zero;
    eta = zero;
    weta = zero;
    xi = zero;
    wxi = zero;
    if (ndimen > 1) {
        cmom = nmom*(nmom + 1)/2;
        noct = 4;
    };
    if (ndimen > 2) {
        cmom = std::pow(nmom, 2);
        noct = 8;
    };
    // FIXME: allocate(ec, lma, );
    if (istat > 0) {
        return;
    };
    ec = zero;
    lma = 0;
    ndpwds = ndpwds + mu.extent(0) + w.extent(0) + wmu.extent(0) + eta.extent(0) + weta.extent(0) + xi.extent(0) + wxi.extent(0) + ec.extent(0) + lma.extent(0);
}

void sn_deallocate()
{
    // FIXME: deallocate(mu, w, eta, xi, wmu, weta, wxi, ec, lma, );
}

void sn_expcoeff(int ndimen)
{
    int id;
    int is;
    int jd;
    int js;
    int kd;
    int ks;
    int l;
    int m;
    int mom;
    int oct;
    // FIXME: select case()
}

void geom_allocate(int nang, int ng, int swp_typ, int ichunk, int &ndpwds, int &ierr)
{
    int i;
    Kokkos::View<int*> indx("indx");
    int ing;
    int j;
    int k;
    int nn;
    ierr = 0;
    // FIXME: allocate(dinv, );
    if (ierr != 0) {
        return;
    };
    hi = zero;
    hj = zero;
    hk = zero;
    dinv = zero;
    ndpwds = ndpwds + dinv.extent(0);
    if (swp_typ == 1) {
        ndiag = ichunk + ny + nz - 2;
        if (ierr != 0) {
            return;
        };
        indx = 0;
        for (k=1; k<=nz; k++) {
            for (j=1; j<=ny; j++) {
                for (i=1; i<=ichunk; i++) {
                    nn = i + j + k - 2;
                };
            };
        };
        for (nn=1; nn<=ndiag; nn++) {
            if (ierr != 0) {
                return;
            };
        };
        for (k=1; k<=nz; k++) {
            for (j=1; j<=ny; j++) {
                for (i=1; i<=ichunk; i++) {
                    nn = i + j + k - 2;
                    indx[nn-1] = indx[nn-1] + 1;
                    ing = indx[nn-1];
                };
            };
        };
        // FIXME: deallocate(indx, );
    };
    // FIXME: implicit deallocate(indx, );
}

void geom_deallocate(int swp_typ)
{
    int i;
    // FIXME: deallocate(dinv, );
}

void geom_param_calc(int nang, int ichunk, const Kokkos::View<const float[ /* FIXME symbolic dimensions */ ]> &mu, const Kokkos::View<const float[ /* FIXME symbolic dimensions */ ]> &eta, const Kokkos::View<const float[ /* FIXME symbolic dimensions */ ]> &xi, const Kokkos::View<const float[ /* FIXME symbolic dimensions */ ][ /* FIXME symbolic dimensions */ ][ /* FIXME symbolic dimensions */ ]> &cs, float vd, const Kokkos::View<float[ /* FIXME symbolic dimensions */ ][ /* FIXME symbolic dimensions */ ][ /* FIXME symbolic dimensions */ ][ /* FIXME symbolic dimensions */ ][ /* FIXME symbolic dimensions */ ]> &d)
{
    int i;
    int ich;
    int ii;
    int j;
    int k;
    int m;
    hi = two/dx;
    if (ndimen > 1) {
        hj = two/dy;
        if (ndimen > 2) {
            hk = two/dz;
        };
    };
    for (k=1; k<=nz; k++) {
        for (j=1; j<=ny; j++) {
            ii = 0;
            ich = 1;
            for (i=1; i<=nx; i++) {
                ii = ii + 1;
                for (m=1; m<=nang; m++) {
                    d[m,ii,j,k,ich-1] = one/(cs[i,j,k-1] + vd + hi*mu[m-1] + hj*eta[m-1] + hk*xi[m-1]);
                };
                if (ii == ichunk) {
                    ii = 0;
                    ich = ich + 1;
                };
            };
        };
    };
}

void time_summ()
{
    int i;
    std::string star="*";
    tinrmisc = tinners - tinrsrc + tsweeps;
    tslvmisc = tslv - tparam + totrsrc + tinners;
    std::cout << std::endl;
    std::cout << tparset << std::endl;
    std::cout << tinp << std::endl;
    std::cout << tset << std::endl;
    std::cout << tslv << std::endl;
    std::cout << tparam << std::endl;
    std::cout << totrsrc << std::endl;
    std::cout << tinners << std::endl;
    std::cout << tinrsrc << std::endl;
    std::cout << tsweeps << std::endl;
    std::cout << tinrmisc << std::endl;
    std::cout << tslvmisc << std::endl;
    std::cout << tout << std::endl;
}

void wtime(float &time)
{
}

void barrier(int comm)
{
}

void bcast_d_1d(const Kokkos::View<const float[ /* FIXME symbolic dimensions */ ]> &value, int dlen, int comm, int bproc)
{
}

void bcast_d_scalar(float value, int comm, int bproc)
{
}

void bcast_i_1d(const Kokkos::View<const int[ /* FIXME symbolic dimensions */ ]> &value, int ilen, int comm, int bproc)
{
}

void bcast_i_scalar(int value, int comm, int bproc)
{
}

void cartrank(const Kokkos::View<const int[2]> &coord, int &rank, int comm)
{
    rank = 0;
}

void glmax_d(float value, int comm)
{
}

void glmax_d_1d(const Kokkos::View<const float[ /* FIXME symbolic dimensions */ ]> &value, int dlen, int comm)
{
}

void glmax_i(int value, int comm)
{
}

void glmin_d(float value, int comm)
{
}

void glmin_i(int value, int comm)
{
}

void glsum_d(float value, int comm)
{
}

void irecv_d_3d(int proc, int myproc, int d1, int d2, int d3, const Kokkos::View<const float[ /* FIXME symbolic dimensions */ ][ /* FIXME symbolic dimensions */ ][ /* FIXME symbolic dimensions */ ]> &value, int comm, int mtag, int req)
{
}

void isend_d_3d(int proc, int myproc, int d1, int d2, int d3, const Kokkos::View<const float[ /* FIXME symbolic dimensions */ ][ /* FIXME symbolic dimensions */ ][ /* FIXME symbolic dimensions */ ]> &value, int comm, int mtag, int req)
{
}

void pcomm_set()
{
    comm_space = 0;
    ycomm = 0;
    zcomm = 0;
    sproc = 0;
    yproc = 0;
    zproc = 0;
    firsty = true;
    ylop = 0;
    lasty = true;
    yhip = 0;
    firstz = true;
    zlop = 0;
    lastz = true;
    zhip = 0;
}

void pend()
{
}

void pinit(float &t1)
{
    int ierr;
    wtime(t1);
    thread_single = 0;
    thread_funneled = 0;
    thread_serialized = 0;
    thread_multiple = 0;
    comm_snap = 0;
    nproc = 1;
    iproc = 0;
}

void pinit_omp(int &ierr, std::string &error)
{
    int max_threads;
    ierr = 0;
    error = " ";
    max_threads = 1;
    nthreads = 1;
    do_nested = false;
}

void plock_omp(std::string dowhat, int nlock)
{
}

void precv_d_2d(int proc, int myproc, int d1, int d2, const Kokkos::View<const float[ /* FIXME symbolic dimensions */ ][ /* FIXME symbolic dimensions */ ]> &value, int comm, int mtag)
{
}

void precv_d_3d(int proc, int myproc, int d1, int d2, int d3, const Kokkos::View<const float[ /* FIXME symbolic dimensions */ ][ /* FIXME symbolic dimensions */ ][ /* FIXME symbolic dimensions */ ]> &value, int comm, int mtag)
{
}

void psend_d_2d(int proc, int myproc, int d1, int d2, const Kokkos::View<const float[ /* FIXME symbolic dimensions */ ][ /* FIXME symbolic dimensions */ ]> &value, int comm, int mtag)
{
}

void psend_d_3d(int proc, int myproc, int d1, int d2, int d3, const Kokkos::View<const float[ /* FIXME symbolic dimensions */ ][ /* FIXME symbolic dimensions */ ][ /* FIXME symbolic dimensions */ ]> &value, int comm, int mtag)
{
}

void rtsum_d_1d(const Kokkos::View<const float[ /* FIXME symbolic dimensions */ ]> &value, int dlen, int comm, int rtproc)
{
}

void tests(int req)
{
}

void testsome(const Kokkos::View<const int[ /* FIXME symbolic dimensions */ ]> &req, const Kokkos::View<const int[ /* FIXME symbolic dimensions */ ]> &indx, int d1, int outcount)
{
}

int thread_num()
{
    int thread_num;
    thread_num = 0;
    return thread_num;
}

void waitall(const Kokkos::View<const int[ /* FIXME symbolic dimensions */ ]> &req, int d1)
{
}

void waitinit(const Kokkos::View<const int[ /* FIXME symbolic dimensions */ ]> &req, int d1)
{
}

void waits(int req)
{
}

void waitsome(const Kokkos::View<const int[ /* FIXME symbolic dimensions */ ]> &req, const Kokkos::View<const int[ /* FIXME symbolic dimensions */ ]> &indx, int d1, int outcount)
{
}

void mms_allocate(int &ierr, std::string &error)
{
    ierr = 0;
    error = " ";
    // FIXME: allocate(ref_flux, ref_fluxm, );
    if (ierr != 0) {
        return;
    };
    ref_flux = zero;
    ref_fluxm = zero;
    // FIXME: allocate(ib, jb, kb, );
    if (ierr != 0) {
        return;
    };
    ib = zero;
    jb = zero;
    kb = zero;
}

void mms_cells()
{
    int i;
    int j;
    int k;
    for (i=1; i<=nx; i++) {
        ib[i + 1-1] = ib[i-1] + dx;
    };
    if (ndimen > 1) {
        for (j=1; j<=ny; j++) {
            jb[j + 1-1] = jb[j-1] + dy;
        };
        if (ndimen > 2) {
            for (k=1; k<=nz; k++) {
                kb[k + 1-1] = kb[k-1] + dz;
            };
        };
    };
    if (src_opt == 3) {
        a = pi/lx;
        if (ndimen > 1) {
            b = pi/ly;
        };
        if (ndimen > 2) {
            c = pi/lz;
        };
    };
}

void mms_deallocate()
{
    if (allocated(ref_flux)) {
        // FIXME: deallocate(ref_flux, );
    };
    if (allocated(ref_fluxm)) {
        // FIXME: deallocate(ref_fluxm, );
    };
    if (allocated(ib)) {
        // FIXME: deallocate(ib, );
    };
    if (allocated(jb)) {
        // FIXME: deallocate(jb, );
    };
    if (allocated(kb)) {
        // FIXME: deallocate(kb, );
    };
}

void mms_flux_1()
{
    int g;
    int i;
    int j;
    int k;
    int l;
    int n;
    Kokkos::View<float[ /* FIXME symbolic dimensions */ ]> p("p");
    Kokkos::View<float[ /* FIXME symbolic dimensions */ ]> tx("tx");
    Kokkos::View<float[ /* FIXME symbolic dimensions */ ]> ty("ty");
    Kokkos::View<float[ /* FIXME symbolic dimensions */ ]> tz("tz");
    mms_trigint("COS", nx, a, dx, ib, tx);
    if (ndimen > 1) {
        mms_trigint("COS", ny, b, dy, jb, ty);
        if (ndimen > 2) {
            mms_trigint("COS", nz, c, dz, kb, tz);
        } else {
            tz = one;
        };
    } else {
        ty = one;
        tz = one;
    };
    for (g=1; g<=ng; g++) {
        for (k=1; k<=nz; k++) {
            for (j=1; j<=ny; j++) {
                for (i=1; i<=nx; i++) {
                    ref_flux[i,j,k,g-1] = g*tx[i-1]*ty[j-1]*tz[k-1];
                };
            };
        };
    };
    p = zero;
    for (n=1; n<=noct; n++) {
        for (l=1; l<=cmom - 1; l++) {
            p[l-1] = p[l-1] + sum(w*ec[/* FIXME right index */,l + 1,n-1]);
        };
    };
    for (l=1; l<=cmom - 1; l++) {
        ref_fluxm[l,/* FIXME right index */,/* FIXME right index */,/* FIXME right index */,/* FIXME right index */-1] = p[l-1]*ref_flux;
    };
}

void mms_flux_1_2()
{
    float t;
    t = tf - half*dt;
    ref_flux = t*ref_flux;
}

void mms_setup(int &ierr, std::string &error)
{
    ierr = 0;
    error = " ";
    mms_allocate(ierr, error);
    glmax_i(ierr, comm_snap);
    if (ierr != 0) {
        error = "***ERROR: MMS_ALLOCATE: Error allocating MMS arrays";
        return;
    };
    ib[1-1] = zero;
    jb[1-1] = yproc*ny*dy;
    kb[1-1] = zproc*nz*dz;
    a = zero;
    b = zero;
    c = zero;
    mms_cells();
    if (src_opt == 3) {
        mms_flux_1();
        mms_src_1();
        if (timedep == 1) {
            mms_flux_1_2();
        };
    };
}

void mms_src_1()
{
    Kokkos::View<float[ /* FIXME symbolic dimensions */ ]> cx("cx");
    Kokkos::View<float[ /* FIXME symbolic dimensions */ ]> cy("cy");
    Kokkos::View<float[ /* FIXME symbolic dimensions */ ]> cz("cz");
    int g;
    int gp;
    int i;
    int id;
    int is;
    int j;
    int jd;
    int js;
    int k;
    int kd;
    int ks;
    int l;
    int ll;
    int lm;
    int m;
    int n;
    Kokkos::View<float[ /* FIXME symbolic dimensions */ ]> sx("sx");
    Kokkos::View<float[ /* FIXME symbolic dimensions */ ]> sy("sy");
    Kokkos::View<float[ /* FIXME symbolic dimensions */ ]> sz("sz");
    mms_trigint("COS", nx, a, dx, ib, cx);
    mms_trigint("SIN", nx, a, dx, ib, sx);
    if (ndimen > 1) {
        mms_trigint("COS", ny, b, dy, jb, cy);
        mms_trigint("SIN", ny, b, dy, jb, sy);
        if (ndimen > 2) {
            mms_trigint("COS", nz, c, dz, kb, cz);
            mms_trigint("SIN", nz, c, dz, kb, sz);
        } else {
            cz = one;
            sz = zero;
        };
    } else {
        cy = one;
        sy = zero;
        cz = one;
        sz = zero;
    };
    for (g=1; g<=ng; g++) {
        for (kd=1; kd<=max(ndimen - 1, 1); kd++) {
            ks = (int)(-one);
            if (kd == 2) {
                ks = (int)(one);
            };
            for (jd=1; jd<=min(ndimen, 2); jd++) {
                js = (int)(-one);
                if (jd == 2) {
                    js = (int)(one);
                };
                for (id=1; id<=2; id++) {
                    is = (int)(-one);
                    if (id == 2) {
                        is = (int)(one);
                    };
                    n = 4*(kd - 1) + 2*(jd - 1) + id;
                    for (k=1; k<=nz; k++) {
                        for (j=1; j<=ny; j++) {
                            for (i=1; i<=nx; i++) {
                                for (m=1; m<=nang; m++) {
                                    qim[m,i,j,k,n,g-1] = qim[m,i,j,k,n,g-1] + g*is*mu[m-1]*sx[i-1]*cy[j-1]*cz[k-1] + sigt[mat[i,j,k-1],g-1]*ref_flux[i,j,k,g-1];
                                    if (ndimen > 1) {
                                        qim[m,i,j,k,n,g-1] = qim[m,i,j,k,n,g-1] + g*js*eta[m-1]*cx[i-1]*sy[j-1]*cz[k-1];
                                    };
                                    if (ndimen > 2) {
                                        qim[m,i,j,k,n,g-1] = qim[m,i,j,k,n,g-1] + g*ks*xi[m-1]*cx[i-1]*cy[j-1]*sz[k-1];
                                    };
                                    for (gp=1; gp<=ng; gp++) {
                                        qim[m,i,j,k,n,g-1] = qim[m,i,j,k,n,g-1] - slgg[mat[i,j,k-1],1,gp,g-1]*ref_flux[i,j,k,gp-1];
                                        lm = 2;
                                        for (l=2; l<=nmom; l++) {
                                            for (ll=1; ll<=lma[l-1]; ll++) {
                                                qim[m,i,j,k,n,g-1] = qim[m,i,j,k,n,g-1] - ec[m,lm,n-1]*slgg[mat[i,j,k-1],l,gp,g-1]*ref_fluxm[lm - 1,i,j,k,g-1];
                                                lm = lm + 1;
                                            };
                                        };
                                    };
                                };
                            };
                        };
                    };
                };
            };
        };
    };
    if (timedep == 1) {
        for (g=1; g<=ng; g++) {
            qi[/* FIXME right index */,/* FIXME right index */,/* FIXME right index */,g-1] = ref_flux[/* FIXME right index */,/* FIXME right index */,/* FIXME right index */,g-1]/v[g-1];
        };
    };
}

void mms_trigint(std::string trig, int lc, float d, float del, const Kokkos::View<const float[ /* FIXME symbolic dimensions */ ]> &cb, const Kokkos::View<float[ /* FIXME symbolic dimensions */ ]> &fn)
{
    int i;
    fn = zero;
    if (true) {
        for (i=1; i<=lc; i++) {
            fn[i-1] = dcos(d*cb[i-1]) - dcos(d*cb[i + 1-1]);
        };
        fn = fn/d*del;
    } else {
        if (false) {
            for (i=1; i<=lc; i++) {
                fn[i-1] = dsin(d*cb[i + 1-1]) - dsin(d*cb[i-1]);
            };
            fn = fn/del;
        } else {
            return;
        };
    };
}

void mms_verify_1(const Kokkos::View<const float[ /* FIXME symbolic dimensions */ ][ /* FIXME symbolic dimensions */ ][ /* FIXME symbolic dimensions */ ][ /* FIXME symbolic dimensions */ ]> &flux)
{
    Kokkos::View<float[ /* FIXME symbolic dimensions */ ][ /* FIXME symbolic dimensions */ ][ /* FIXME symbolic dimensions */ ][ /* FIXME symbolic dimensions */ ]> df("df");
    float dfmn;
    float dfmx;
    float dfsm;
    int i;
    std::string star="*";
    df = one;
    df = dabs(flux/ref_flux - df);
    dfmx = maxval(df);
    glmax_d(dfmx, comm_snap);
    dfmn = minval(df);
    glmin_d(dfmn, comm_snap);
    dfsm = sum(df);
    glsum_d(dfsm, comm_snap);
    dfsm = dfsm/nx*ny_gl*nz_gl*ng;
    if (iproc == root) {
        std::cout << dfmx << std::endl;
        std::cout << dfmn << std::endl;
        std::cout << dfsm << std::endl;
    };
}

void solvar_allocate(int &ndpwds, int &ierr)
{
    int szcor;
    ierr = 0;
    if (timedep == 1) {
        if (angcpy == 1) {
            // FIXME: allocate(ptr_in, );
        } else {
            // FIXME: allocate(ptr_in, ptr_out, );
        };
    } else {
        // FIXME: allocate(ptr_in, ptr_out, );
    };
    if (ierr != 0) {
        return;
    };
    if (timedep == 1) {
        ptr_in = zero;
        if (angcpy == 1) {
        } else {
            ptr_out = zero;
        };
    };
    // FIXME: allocate(flux0, flux0po, flux0pi, );
    if (ierr != 0) {
        return;
    };
    if (cmom > 1) {
        // FIXME: allocate(fluxm, );
        if (ierr != 0) {
            return;
        };
    } else {
        // FIXME: allocate(fluxm, );
        if (ierr != 0) {
            return;
        };
    };
    flux0 = zero;
    flux0po = zero;
    flux0pi = zero;
    fluxm = zero;
    // FIXME: allocate(q2grp0, qtot, );
    if (ierr != 0) {
        return;
    };
    if (cmom > 1) {
        // FIXME: allocate(q2grpm, );
        if (ierr != 0) {
            return;
        };
    } else {
        // FIXME: allocate(q2grpm, );
        if (ierr != 0) {
            return;
        };
    };
    q2grp0 = zero;
    q2grpm = zero;
    qtot = zero;
    // FIXME: allocate(t_xs, a_xs, s_xs, );
    if (ierr != 0) {
        return;
    };
    t_xs = zero;
    a_xs = zero;
    s_xs = zero;
    if (multiswp == 0) {
        szcor = 1;
    } else {
        szcor = ncor;
    };
    // FIXME: allocate(psii, psij, psik, );
    if (ierr != 0) {
        return;
    };
    psii = zero;
    psij = zero;
    psik = zero;
    // FIXME: allocate(jb_in, jb_out, kb_in, kb_out, );
    if (ierr != 0) {
        return;
    };
    jb_in = zero;
    jb_out = zero;
    kb_in = zero;
    kb_out = zero;
    // FIXME: allocate(fixup_counter, );
    if (ierr != 0) {
        return;
    };
    fixup_counter = zero;
    // FIXME: allocate(flkx, flky, flkz, );
    if (ierr != 0) {
        return;
    };
    flkx = zero;
    flky = zero;
    flkz = zero;
    // FIXME: allocate(fmin, fmax, pop, );
    if (ierr != 0) {
        return;
    };
    fmin = zero;
    fmax = zero;
    pop = zero;
    ndpwds = ndpwds + ptr_in.extent(0) + flux0.extent(0) + flux0po.extent(0) + flux0pi.extent(0) + fluxm.extent(0) + q2grp0.extent(0) + q2grpm.extent(0) + qtot.extent(0) + t_xs.extent(0) + a_xs.extent(0) + s_xs.extent(0) + psii.extent(0) + psij.extent(0) + psik.extent(0) + jb_in.extent(0) + jb_out.extent(0) + kb_in.extent(0) + kb_out.extent(0) + fixup_counter.extent(0) + flkx.extent(0) + flky.extent(0) + flkz.extent(0) + fmin.extent(0) + fmax.extent(0) + pop.extent(0);
    if ((angcpy == 2) || (timedep == 0)) {
        ndpwds = ndpwds + ptr_out.extent(0);
    };
}

void solvar_deallocate()
{
    // FIXME: deallocate(ptr_in, );
    if ((angcpy == 2) || (timedep == 0)) {
        // FIXME: deallocate(ptr_out, );
    };
    // FIXME: deallocate(flux0, flux0po, flux0pi, fluxm, );
    // FIXME: deallocate(q2grp0, q2grpm, qtot, );
    // FIXME: deallocate(t_xs, a_xs, s_xs, );
    // FIXME: deallocate(psii, psij, psik, );
    // FIXME: deallocate(jb_in, jb_out, kb_in, kb_out, );
    // FIXME: deallocate(fixup_counter, );
    // FIXME: deallocate(flkx, flky, flkz, );
    // FIXME: deallocate(fmin, fmax, pop, );
}

void dealloc_input(int flg)
{
    sn_deallocate();
    if (flg > 1) {
        data_deallocate();
    };
    if (flg > 2) {
        control_deallocate();
    };
    if (flg > 3) {
        mms_deallocate();
    };
}

void dealloc_solve(int swp_typ, int flg)
{
    geom_deallocate(swp_typ);
    if (flg > 1) {
        solvar_deallocate();
    };
}

void close_file(int funit, int &ierr, std::string &error)
{
    ierr = 0;
    error = "";
    if (iproc != root) {
        return;
    };
    if (ierr != 0) {
        error = "***ERROR: CLOSE_FILE: Unable to close file, unit:";
        std::cout << trim(error) << " " << funit << std::endl;
    };
}

void cmdarg(int &ierr, std::string &error)
{
    std::string arg;
    int n;
    int narg;
    if (iproc != root) {
        return;
    };
    ierr = 0;
    error = "";
    narg = 0;
    if (narg != 2) {
        ierr = 1;
        error = "***ERROR: CMDARG: Missing command line entry";
        return;
    };
    for (n=1; n<=2; n++) {
        if (((arg[1-1] == "-") || (arg[1-1] == "<")) || (arg[1-1] == ">")) {
            ierr = 1;
            error = "***ERROR: CMDARG: Bad command line entry, arg:";
            std::cout << trim(error) << " " << n << std::endl;
        } else {
            if (n == 1) {
                ifile = arg;
            } else {
                if (n == 2) {
                    ofile = arg;
                };
            };
        };
    };
}

void open_file(int funit, std::string fname, std::string fstat, std::string faction, int &ierr, std::string &error)
{
    std::string tname;
    ierr = 0;
    error = "";
    if (iproc != root) {
        return;
    };
    tname = trim(fname);
    if (ierr != 0) {
        error = "***ERROR: OPEN_FILE: Unable to open file, unit:";
        std::cout << trim(error) << " " << funit << std::endl;
    };
}

void print_error(int funit, std::string error)
{
    if (iproc != root) {
        return;
    };
    if (funit > 0) {
        std::cout << error << std::endl;
    };
}

void stop_run(int flg1, int flg2, int flg3, int flg4)
{
    if (flg1 > 0) {
        plock_omp("destroy", nthreads);
    };
    if (flg2 > 0) {
        dealloc_input(flg2);
    };
    if (flg3 > 0) {
        dealloc_solve(swp_typ, flg3);
    };
    if (iproc == root) {
        if (flg4 == 0) {
            std::cout << "Aww SNAP. Program failed. Try again." << std::endl;
        } else {
            if (flg4 == 1) {
                std::cout << "Success! Done in a SNAP!" << std::endl;
            } else {
                if (flg4 == 2) {
                    std::cout << ("Oh SNAP. That did not converge. But ") + ("take a look at the Timing Summary anyway!") << std::endl;
                };
            };
        };
    };
    pend();
    exit(0);
}

void input_check(int &ierr)
{
    std::string error;
    ierr = 0;
    if (npey < 1) {
        ierr = ierr + 1;
        error = "***ERROR: INPUT_CHECK: NPEY must be positive";
        print_error(ounit, error);
    };
    if (npez < 1) {
        ierr = ierr + 1;
        error = "***ERROR: INPUT_CHECK: NPEZ must be positive";
        print_error(ounit, error);
    };
    if ((ndimen < 2) && (npey != 1)) {
        ierr = ierr + 1;
        error = "***ERROR: INPUT_CHECK: NPEY must be 1 if not 2-D problem";
        print_error(ounit, error);
    };
    if ((ndimen < 3) && (npez != 1)) {
        ierr = ierr + 1;
        error = "***ERROR: INPUT_CHECK: NPEZ must be 1 if not 3-D problem";
        print_error(ounit, error);
    };
    if (npey*npez != nproc) {
        ierr = ierr + 1;
        error = "***ERROR: INPUT_CHECK: NPEY*NPEZ must equal MPI NPROC";
        print_error(ounit, error);
    };
    if (ichunk > nx) {
        ichunk = nx;
        error = ("*WARNING: INPUT_CHECK: ICHUNK cannot exceed NX; ") + ("reset to NX");
        print_error(ounit, error);
    };
    if ((ndimen == 1) && (ichunk != nx)) {
        ichunk = nx;
    };
    if (nthreads < 1) {
        nthreads = 1;
        error = ("*WARNING: INPUT_CHECK: NTHREADS must be positive; ") + ("reset to 1");
        print_error(ounit, error);
    };
    if (nnested < 1) {
        nnested = 1;
        error = ("*WARNING: INPUT_CHECK: NNESTED must be positive; ") + ("reset to 1");
        print_error(ounit, error);
    };
    if ((ndimen < 1) || (ndimen > 3)) {
        ierr = ierr + 1;
        error = "***ERROR: INPUT_CHECK: NDIMEN must be 1, 2, or 3";
        print_error(ounit, error);
    };
    if (lx <= zero) {
        ierr = ierr + 1;
        error = "***ERROR: INPUT_CHECK: LX must be positive";
        print_error(ounit, error);
    };
    if (ly <= zero) {
        ierr = ierr + 1;
        error = "***ERROR: INPUT_CHECK: LY must be positive";
        print_error(ounit, error);
    };
    if (lz <= zero) {
        ierr = ierr + 1;
        error = "***ERROR: INPUT_CHECK: LZ must be positive";
        print_error(ounit, error);
    };
    if ((ndimen < 2) && ((ny > 1) || (ly != one))) {
        ny = 1;
        ly = zero;
        error = "*WARNING: INPUT_CHECK: NY and LY reset for 1-D problem";
        print_error(ounit, error);
    };
    if ((ndimen < 3) && ((nz > 1) || (lz != one))) {
        nz = 1;
        lz = zero;
        error = "*WARNING: INPUT_CHECK: NZ and LZ reset for 1/2-D problem";
        print_error(ounit, error);
    };
    if ((nmom < 1) || (nmom > 4)) {
        ierr = ierr + 1;
        error = ("***ERROR: INPUT_CHECK: NMOM must be positive and ") + ("less than 5");
        print_error(ounit, error);
    };
    if (nang < 1) {
        ierr = ierr + 1;
        error = "***ERROR: INPUT_CHECK: NANG must be positive";
        print_error(ounit, error);
    };
    if (ng < 1) {
        ierr = ierr + 1;
        error = "***ERROR: INPUT_CHECK: NG must be positive";
        print_error(ounit, error);
    };
    if (nthreads > ng) {
        nthreads = ng;
        error = ("*WARNING: INPUT_CHECK: NTHREADS should be <= NG; ") + ("reset to NG");
        print_error(ounit, error);
    };
    if ((mat_opt < 0) || (mat_opt > 2)) {
        ierr = ierr + 1;
        error = "***ERROR: INPUT_CHECK: MAT_OPT must be 0/1/2";
        print_error(ounit, error);
    };
    if ((src_opt < 0) || (src_opt > 3)) {
        ierr = ierr + 1;
        error = "***ERROR: INPUT_CHECK: SRC_OPT must be 0/1/2/3";
        print_error(ounit, error);
    };
    if ((scatp != 0) && (scatp != 1)) {
        scatp = 0;
        error = "*WARNING: INPUT_CHECK: SCATP must be 0/1; reset to 0";
        print_error(ounit, error);
    };
    if ((epsi <= zero) || (epsi >= 0.010000)) {
        ierr = ierr + 1;
        error = ("***ERROR: INPUT_CHECK: EPSI must be positive, less ") + ("than 1.0E-2");
        print_error(ounit, error);
    };
    if (iitm < 1) {
        ierr = ierr + 1;
        error = "***ERROR: INPUT_CHECK: IITM must be positive";
        print_error(ounit, error);
    };
    if (oitm < 1) {
        ierr = ierr + 1;
        error = "***ERROR: INPUT_CHECK: OITM must be positive";
        print_error(ounit, error);
    };
    if ((timedep != 0) && (timedep != 1)) {
        timedep = 0;
        error = "*WARNING: INPUT_CHECK: TIMEDEP must be 0/1; reset to 0";
        print_error(ounit, error);
    };
    if ((timedep == 1) && (tf <= zero)) {
        ierr = ierr + 1;
        error = "***ERROR: INPUT_CHECK: TF must be positive";
        print_error(ounit, error);
    };
    if ((timedep == 1) && (nsteps < 1)) {
        ierr = ierr + 1;
        error = "***ERROR: INPUT_CHECK: NSTEPS must be positive";
        print_error(ounit, error);
    };
    if ((timedep == 0) && (nsteps != 1)) {
        nsteps = 1;
        error = ("*WARNING: INPUT_CHECK: NSTEPS reset to 1 for ") + ("static calc");
        print_error(ounit, error);
    };
    if ((swp_typ != 0) && (swp_typ != 1)) {
        swp_typ = 0;
        error = ("*WARNING: INPUT_CHECK: SWP_TYP must equal 0/1; ") + ("reset to 0");
        print_error(ounit, error);
    };
    if ((swp_typ != 0) && (ndimen == 1)) {
        swp_typ = 0;
        error = ("*WARNING: INPUT_CHECK: SWP_TYP must be 0 for 1-D; ") + ("reset to 0");
        print_error(ounit, error);
    };
    if ((nnested > 1) && (swp_typ == 0)) {
        ierr = ierr + 1;
        error = ("***ERROR: INPUT_CHECK: SWP_TYP must be 1 for nested") + (" threading");
        print_error(ounit, error);
    };
    if ((multiswp != 0) && (multiswp != 1)) {
        multiswp = 0;
        error = ("*WARNING: INPUT_CHECK: MULTISWP must equal 0/1; ") + ("reset to 0");
        print_error(ounit, error);
    };
    if ((ndimen == 1) && (multiswp == 1)) {
        multiswp = 0;
        error = ("*WARNING: INPUT_CHECK: MULTISWP must equal 0 in 1D;") + (" reset to 0");
        print_error(ounit, error);
    };
    if ((multiswp == 1) && (swp_typ == 1)) {
        ierr = ierr + 1;
        error = ("***ERROR: INPUT_CHECK: MULTISWP must equal 0 if ") + ("SWP_TYP is 1");
        print_error(ounit, error);
    };
    if ((angcpy != 1) && (angcpy != 2)) {
        angcpy = 1;
        error = ("*WARNING: INPUT_CHECK: ANGCPY must be 1 or 2; ") + ("reset to 1");
        print_error(ounit, error);
    };
    if ((it_det != 0) && (it_det != 1)) {
        it_det = 0;
        error = ("*WARNING: INPUT_CHECK: IT_DET must equal 0/1; ") + ("reset to 0");
        print_error(ounit, error);
    };
    if ((fluxp < 0) || (fluxp > 2)) {
        fluxp = 0;
        error = ("*WARNING: INPUT_CHECK: FLUXP must equal 0/1/2; ") + ("reset to 0");
        print_error(ounit, error);
    };
    if ((fixup != 0) && (fixup != 1)) {
        fixup = 0;
        error = ("*WARNING: INPUT_CHECK: FIXUP must equal 0/1; ") + ("reset to 0");
        print_error(ounit, error);
    };
    if ((soloutp != 0) && (soloutp != 1)) {
        soloutp = 0;
        error = ("*WARNING: INPUT_CHECK: SOLOUTP must equal 0/1; ") + ("reset to 0");
        print_error(ounit, error);
    };
    if ((kplane < 0) || (kplane > nz)) {
        kplane = 0;
        error = ("*WARNING: INPUT_CHECK: KPLANE must be in the range ") + ("0 to NZ; reset to 0");
        print_error(ounit, error);
    };
    if ((popout < 0) || (popout > 2)) {
        popout = 0;
        error = ("*WARNING: INPUT_CHECK: POPOUT must equal 0/1/2; ") + ("reset to 0");
        print_error(ounit, error);
    };
}

void input_echo()
{
    int i;
    std::string star="*";
    std::cout <<      /* FIXME: implied do loop */  << std::endl;
    std::cout <<      /* FIXME: implied do loop */  << std::endl;
    std::cout << std::endl;
    std::cout << npey << npez << ichunk << nthreads << nnested << std::endl;
    std::cout << ndimen << nx << ny << nz << std::endl;
    std::cout << lx << ly << lz << std::endl;
    std::cout << nmom << nang << std::endl;
    std::cout << ng << mat_opt << src_opt << scatp << std::endl;
    std::cout << epsi << iitm << oitm << timedep << tf << nsteps << swp_typ << multiswp << angcpy << it_det << soloutp << kplane << popout << fluxp << fixup << std::endl;
    std::cout <<      /* FIXME: implied do loop */  << std::endl;
}

void input_read()
{
    std::string error;
    int ierr;
    float t1;
    float t2;
}

void input_var_bcast()
{
    int dlen;
    Kokkos::View<float[15]> dpak("dpak");
    int ilen;
    Kokkos::View<int[30]> ipak("ipak");
    ipak = 0;
    dpak = zero;
    if (iproc == root) {
        ipak[1-1] = npey;
        ipak[2-1] = npez;
        ipak[3-1] = ichunk;
        ipak[4-1] = nthreads;
        ipak[5-1] = ndimen;
        ipak[6-1] = nx;
        ipak[7-1] = ny;
        ipak[8-1] = nz;
        ipak[9-1] = nmom;
        ipak[10-1] = nang;
        ipak[11-1] = ng;
        ipak[12-1] = timedep;
        ipak[13-1] = nsteps;
        ipak[14-1] = iitm;
        ipak[15-1] = oitm;
        ipak[16-1] = mat_opt;
        ipak[17-1] = src_opt;
        ipak[18-1] = scatp;
        ipak[19-1] = swp_typ;
        ipak[20-1] = multiswp;
        ipak[21-1] = angcpy;
        ipak[22-1] = it_det;
        ipak[23-1] = fluxp;
        ipak[24-1] = fixup;
        ipak[25-1] = nnested;
        ipak[26-1] = soloutp;
        ipak[27-1] = kplane;
        ipak[28-1] = popout;
        dpak[1-1] = lx;
        dpak[2-1] = ly;
        dpak[3-1] = lz;
        dpak[4-1] = tf;
        dpak[5-1] = epsi;
    };
    ilen = ipak.extent(0);
    dlen = dpak.extent(0);
    bcast_i_1d(ipak, ilen, comm_snap, root);
    bcast_d_1d(dpak, dlen, comm_snap, root);
    if (iproc != root) {
        npey = ipak[1-1];
        npez = ipak[2-1];
        ichunk = ipak[3-1];
        nthreads = ipak[4-1];
        ndimen = ipak[5-1];
        nx = ipak[6-1];
        ny = ipak[7-1];
        nz = ipak[8-1];
        nmom = ipak[9-1];
        nang = ipak[10-1];
        ng = ipak[11-1];
        timedep = ipak[12-1];
        nsteps = ipak[13-1];
        iitm = ipak[14-1];
        oitm = ipak[15-1];
        mat_opt = ipak[16-1];
        src_opt = ipak[17-1];
        scatp = ipak[18-1];
        swp_typ = ipak[19-1];
        multiswp = ipak[20-1];
        angcpy = ipak[21-1];
        it_det = ipak[22-1];
        fluxp = ipak[23-1];
        fixup = ipak[24-1];
        nnested = ipak[25-1];
        soloutp = ipak[26-1];
        kplane = ipak[27-1];
        popout = ipak[28-1];
        lx = dpak[1-1];
        ly = dpak[2-1];
        lz = dpak[3-1];
        tf = dpak[4-1];
        epsi = dpak[5-1];
    };
    barrier(comm_snap);
}

void output()
{
    Kokkos::View<int[2]> co("co");
    std::string error;
    Kokkos::View<float**> fprnt("fprnt");
    int g;
    int i;
    int ierr;
    int ii;
    int is;
    int j;
    int jlb;
    int jp;
    int jub;
    int k;
    int klb;
    int kloc;
    int kub;
    int mtag;
    int rank;
    std::string star="*";
    float t1;
    float t2;
    wtime(t1);
    ierr = 0;
    error = " ";
    klb = zproc*nz + 1;
    kub = (zproc + 1)*nz;
    if (soloutp == 1) {
        if (iproc == root) {
            std::cout <<              /* FIXME: implied do loop */  << std::endl;
            // FIXME: allocate(fprnt, );
            fprnt = zero;
        };
        bcast_i_scalar(ierr, comm_snap, root);
        if (ierr != 0) {
            error = "***ERROR: OUTPUT: Allocation error";
            print_error(ounit, error);
            stop_run(1, 4, 2, 0);
        };
        k = nz_gl/2 + 1;
        if (kplane != 0) {
            k = kplane;
        };
        for (g=1; g<=ng; g++) {
            if ((klb <= k) && (k <= kub)) {
                mtag = g;
                output_send(mtag, flux0[/* FIXME right index */,/* FIXME right index */,kloc,g-1]);
            };
            if (iproc == root) {
                co[1-1] = (k - 1)/nz;
                fprnt[/* FIXME right index */,ny-1] = flux0[/* FIXME right index */,/* FIXME right index */,kloc,g-1];
                std::cout <<                  /* FIXME: implied do loop */  << g << k <<                  /* FIXME: implied do loop */  << std::endl;
                for (jp=0; jp<=npey - 1; jp++) {
                    jlb = jp*ny + 1;
                    jub = (jp + 1)*ny;
                    co[2-1] = jp;
                    cartrank(co, rank, comm_space);
                    mtag = g;
                    output_recv(mtag, rank, fprnt[/* FIXME right index */,jub-1]);
                };
                for (i=1; i<=nx; i+=6) {
                    is = i + 6 - 1;
                    if (is > nx) {
                        is = nx;
                    };
                    for (ii=i; ii<=is; ii++) {
                    };
                    for (j=ny_gl; j>=1; j--) {
                    };
                };
                std::cout <<                  /* FIXME: implied do loop */  << std::endl;
            };
        };
        if (iproc == root) {
            // FIXME: deallocate(fprnt, );
        };
    };
    barrier(comm_snap);
    if (fluxp > 0) {
        output_flux_file(klb, kub);
    };
    if (src_opt == 3) {
        mms_verify_1(flux0);
    };
    wtime(t2);
    tout = t2 - t1;
    // FIXME: implicit deallocate(fprnt, );
}

void output_flux_file(int klb, int kub)
{
    Kokkos::View<int[2]> co("co");
    std::string error;
    std::string fn="flux";
    Kokkos::View<float**> fprnt("fprnt");
    int fu;
    int g;
    int i;
    int ierr;
    int j;
    int jlb;
    int jp;
    int jub;
    int k;
    int kloc;
    int l;
    int mtag;
    int rank;
    ierr = 0;
    error = " ";
    open_file(fu, fn, "REPLACE", "WRITE", ierr, error);
    bcast_i_scalar(ierr, comm_snap, root);
    if (ierr != 0) {
        print_error(ounit, error);
        stop_run(1, 4, 2, 0);
    };
    if (iproc == root) {
        // FIXME: allocate(fprnt, );
        fprnt = zero;
    };
    bcast_i_scalar(ierr, comm_snap, root);
    if (ierr != 0) {
        error = "***ERROR: OUTPUT_FLUX_FILE: Allocation error";
        print_error(ounit, error);
        stop_run(1, 4, 2, 0);
    };
    if (iproc == root) {
        if (fluxp == 1) {
            std::cout << std::endl;
        } else {
            std::cout << std::endl;
        };
    };
    for (l=1; l<=max(1, (fluxp - 1)*cmom); l++) {
        if (iproc == root) {
            std::cout << l << std::endl;
        };
        for (g=1; g<=ng; g++) {
            for (k=1; k<=nz_gl; k++) {
                if ((klb <= k) && (k <= kub)) {
                    mtag = l*ng*nz + (g - 1)*nz + kloc;
                    if (l == 1) {
                        output_send(mtag, flux0[/* FIXME right index */,/* FIXME right index */,kloc,g-1]);
                    } else {
                        output_send(mtag, fluxm[l - 1,/* FIXME right index */,/* FIXME right index */,kloc,g-1]);
                    };
                };
                if (iproc == root) {
                    co[1-1] = (k - 1)/nz;
                    fprnt[/* FIXME right index */,ny-1] = flux0[/* FIXME right index */,/* FIXME right index */,kloc,g-1];
                    for (jp=0; jp<=npey - 1; jp++) {
                        jlb = jp*ny + 1;
                        jub = (jp + 1)*ny;
                        co[2-1] = jp;
                        cartrank(co, rank, comm_space);
                        mtag = l*ng*nz + (g - 1)*nz + kloc;
                        output_recv(mtag, rank, fprnt[/* FIXME right index */,jub-1]);
                    };
                    std::cout <<                      /* FIXME: implied do loop */  << std::endl;
                };
            };
        };
    };
    if (iproc == root) {
        // FIXME: deallocate(fprnt, );
    };
    close_file(fu, ierr, error);
    bcast_i_scalar(ierr, comm_snap, root);
    if (ierr != 0) {
        print_error(ounit, error);
        stop_run(1, 4, 2, 0);
    };
    // FIXME: implicit deallocate(fprnt, );
}

void output_recv(int mtag, int proc, const Kokkos::View<float[ /* FIXME symbolic dimensions */ ][ /* FIXME symbolic dimensions */ ]> &fprnt)
{
    precv_d_2d(proc, sproc, nx, ny, fprnt, comm_space, mtag);
}

void output_send(int mtag, const Kokkos::View<const float[ /* FIXME symbolic dimensions */ ][ /* FIXME symbolic dimensions */ ]> &fprnt)
{
    psend_d_2d(root, sproc, nx, ny, fprnt, comm_space, mtag);
}

void setup(int &ndpwds)
{
    std::string error;
    int flg;
    int hnpy;
    int hnpz;
    int i;
    int idle;
    int ierr;
    int mgpt;
    int mie;
    int mis;
    int mje;
    int mjs;
    int mke;
    int mks;
    int qie;
    int qis;
    int qje;
    int qjs;
    int qke;
    int qks;
    int stgs;
    float t1;
    float t2;
    wtime(t1);
    ny_gl = ny;
    nz_gl = nz;
    ny = ny_gl/npey;
    nz = nz_gl/npez;
    jlb = yproc*ny + 1;
    jub = (yproc + 1)*ny;
    klb = zproc*nz + 1;
    kub = (zproc + 1)*nz;
    nc = nx/ichunk;
    jdim = min(ndimen, 2);
    kdim = max(ndimen - 1, 1);
    ncor = jdim*kdim;
    setup_alloc(ndpwds, flg, ierr, error);
    if (ierr != 0) {
        print_error(ounit, error);
        stop_run(1, flg, 0, 0);
    };
    corner_sch = 0;
    if (ndimen == 3) {
        if (npey <= npez) {
            last_oct = 6;
            corner_sch[/* FIXME right index */,1-1] = from_std_vector<float>({1, 1});
            corner_sch[/* FIXME right index */,2-1] = from_std_vector<float>({2, 1});
            corner_sch[/* FIXME right index */,3-1] = from_std_vector<float>({2, 2});
            corner_sch[/* FIXME right index */,4-1] = from_std_vector<float>({1, 2});
        } else {
            last_oct = 4;
            corner_sch[/* FIXME right index */,1-1] = from_std_vector<float>({1, 1});
            corner_sch[/* FIXME right index */,2-1] = from_std_vector<float>({1, 2});
            corner_sch[/* FIXME right index */,3-1] = from_std_vector<float>({2, 2});
            corner_sch[/* FIXME right index */,4-1] = from_std_vector<float>({2, 1});
        };
    } else {
        last_oct = std::pow(2, ndimen);
        corner_sch[/* FIXME right index */,1-1] = from_std_vector<float>({1, 1});
        corner_sch[/* FIXME right index */,2-1] = from_std_vector<float>({2, 1});
    };
    yzstg = 0;
    corner_loop_order = 0;
    if (multiswp == 1) {
        nops = ng*ncor*2*nc;
        swp_typ = 0;
        yzstg[1-1] = npey - yproc + 1 + npez - zproc + 1;
        yzstg[2-1] = yproc + npez - zproc + 1;
        yzstg[3-1] = npey - yproc + 1 + zproc;
        yzstg[4-1] = yproc + zproc;
        if (ndimen == 1) {
            i = 1;
        } else {
            if (ndimen == 2) {
            };
        };
        // FIXME: select case()
    };
    if (multiswp == 0) {
        if (npey <= npez) {
            idle = 3*(npey - 1) + 2*(npez - 1);
        } else {
            idle = 2*(npey - 1) + 3*(npez - 1);
        };
    } else {
        idle = 2*(hnpy + hnpz - 2);
    };
    stgs = nc*noct*mgpt + idle;
    pce = nc*noct*ng/nthreads*stgs;
    setup_delta();
    setup_vel();
    setup_angle();
    setup_mat(mis, mie, mjs, mje, mks, mke);
    setup_data();
    sn_expcoeff(ndimen);
    setup_src(qis, qie, qjs, qje, qks, qke, ierr, error);
    if (ierr != 0) {
        print_error(ounit, error);
        stop_run(1, 3, 0, 0);
    };
    if (iproc == root) {
        setup_echo(mis, mie, mjs, mje, mks, mke, qis, qie, qjs, qje, qks, qke);
        if (scatp == 1) {
            setup_scatp(ierr, error);
        };
    };
    glmax_i(ierr, comm_snap);
    if (ierr != 0) {
        print_error(ounit, error);
        stop_run(1, 4, 0, 0);
    };
    wtime(t2);
    tset = t2 - t1;
}

void setup_alloc(int &ndpwds, int &flg, int &ierr, std::string &error)
{
    flg = 0;
    ierr = 0;
    error = " ";
    sn_allocate(ndimen, ndpwds, ierr);
    glmax_i(ierr, comm_snap);
    if (ierr != 0) {
        error = "***ERROR: SETUP_ALLOC: Allocation error in SN_ALLOCATE";
        return;
    };
    data_allocate(nx, ny, nz, nmom, nang, noct, timedep, ndpwds, ierr);
    glmax_i(ierr, comm_snap);
    if (ierr != 0) {
        flg = 1;
        error = "***ERROR: SETUP_ALLOC: Allocation error in DATA_ALLOCATE";
        return;
    };
    control_allocate(ng, ndpwds, ierr);
    glmax_i(ierr, comm_snap);
    if (ierr != 0) {
        flg = 2;
        error = ("***ERROR: CONTROL_ALLOC: Allocation error of control ") + ("arrays");
        return;
    };
}

void setup_angle()
{
    float dm;
    int m;
    Kokkos::View<float[ /* FIXME symbolic dimensions */ ]> t("t");
    dm = one/nang;
    t = zero;
    mu[1-1] = half*dm;
    for (m=2; m<=nang; m++) {
        mu[m-1] = mu[m - 1-1] + dm;
    };
    if (ndimen > 1) {
        eta[1-1] = one - half*dm;
        for (m=2; m<=nang; m++) {
            eta[m-1] = eta[m - 1-1] - dm;
        };
        if (ndimen > 2) {
            t = std::pow(mu, 2) + std::pow(eta, 2);
            for (m=1; m<=nang; m++) {
                xi[m-1] = dsqrt(one - t[m-1]);
            };
        };
    };
    if (ndimen == 1) {
        w = half/nang;
    } else {
        if (ndimen == 2) {
            w = 0.250000/nang;
        } else {
            w = 0.125000/nang;
        };
    };
}

void setup_data()
{
    int g;
    int n;
    float t;
    sigt[1,1-1] = one;
    siga[1,1-1] = half;
    sigs[1,1-1] = half;
    if (nmat == 2) {
        sigt[2,1-1] = two;
        siga[2,1-1] = 0.800000;
        sigs[2,1-1] = 1.200000;
    };
    for (g=2; g<=ng; g++) {
        sigt[/* FIXME right index */,g-1] = sigt[/* FIXME right index */,g - 1-1] + 0.010000;
        siga[/* FIXME right index */,g-1] = siga[/* FIXME right index */,g - 1-1] + 0.005000;
        sigs[/* FIXME right index */,g-1] = sigs[/* FIXME right index */,g - 1-1] + 0.005000;
    };
    for (g=1; g<=ng; g++) {
        if (ng == 1) {
            slgg[1,1,1,1-1] = sigs[1,g-1];
            break;
        };
        slgg[1,1,g,g-1] = 0.200000*sigs[1,g-1];
        if (g > 1) {
            t = one/(g - 1);
        } else {
            slgg[1,1,1,1-1] = 0.300000*sigs[1,1-1];
        };
        if (g < ng) {
            t = one/(ng - g);
        } else {
            slgg[1,1,ng,ng-1] = 0.900000*sigs[1,ng-1];
        };
    };
    if (nmat == 2) {
        for (g=1; g<=ng; g++) {
            if (ng == 1) {
                slgg[2,1,1,1-1] = sigs[2,g-1];
                break;
            };
            slgg[2,1,g,g-1] = 0.500000*sigs[2,g-1];
            if (g > 1) {
                t = one/(g - 1);
            } else {
                slgg[2,1,1,1-1] = 0.600000*sigs[2,1-1];
            };
            if (g < ng) {
                t = one/(ng - g);
            } else {
                slgg[2,1,ng,ng-1] = 0.900000*sigs[2,ng-1];
            };
        };
    };
    if (nmom > 1) {
        slgg[1,2,/* FIXME right index */,/* FIXME right index */-1] = 0.100000*slgg[1,1,/* FIXME right index */,/* FIXME right index */-1];
        for (n=3; n<=nmom; n++) {
            slgg[1,n,/* FIXME right index */,/* FIXME right index */-1] = half*slgg[1,n - 1,/* FIXME right index */,/* FIXME right index */-1];
        };
    };
    if ((nmat == 2) && (nmom > 1)) {
        slgg[2,2,/* FIXME right index */,/* FIXME right index */-1] = 0.800000*slgg[2,1,/* FIXME right index */,/* FIXME right index */-1];
        for (n=3; n<=nmom; n++) {
            slgg[2,n,/* FIXME right index */,/* FIXME right index */-1] = 0.600000*slgg[2,n - 1,/* FIXME right index */,/* FIXME right index */-1];
        };
    };
}

void setup_delta()
{
    dx = lx/nx;
    if (ndimen > 1) {
        dy = ly/ny_gl;
    };
    if (ndimen > 2) {
        dz = lz/nz_gl;
    };
    if (timedep == 1) {
        dt = tf/nsteps;
    };
}

void setup_echo(int mis, int mie, int mjs, int mje, int mks, int mke, int qis, int qie, int qjs, int qje, int qks, int qke)
{
    int i;
    int j;
    std::string star="*";
    std::cout <<      /* FIXME: implied do loop */  << std::endl;
    std::cout << std::endl;
    std::cout << ndimen << nx << ny_gl << nz_gl << lx << ly << lz << dx << dy << dz << std::endl;
    std::cout << std::endl;
    std::cout << nmom << nang << noct << w[1-1] << std::endl;
    std::cout << std::endl;
    // FIXME: select case()
    std::cout << std::endl;
    std::cout << mat_opt << nmat << std::endl;
    std::cout << std::endl;
    if (mat_opt > 0) {
        std::cout << mis << mjs << mks << mie << mje << mke << std::endl;
    };
    std::cout << std::endl;
    std::cout << src_opt << std::endl;
    if (src_opt < 3) {
        std::cout << std::endl;
        std::cout << qis << qjs << qks << qie << qje << qke << std::endl;
    } else {
        std::cout << std::endl;
    };
    std::cout << std::endl;
    std::cout << ng << std::endl;
    for (j=1; j<=nmat; j++) {
        std::cout << j << std::endl;
        std::cout << std::endl;
        std::cout <<          /* FIXME: implied do loop */  << std::endl;
    };
    if (timedep == 1) {
        std::cout << std::endl;
        std::cout << tf << nsteps << dt << std::endl;
        std::cout << std::endl;
        std::cout <<          /* FIXME: implied do loop */  << std::endl;
    };
    std::cout << std::endl;
    std::cout << epsi << iitm << oitm << timedep << swp_typ << multiswp << angcpy << it_det << soloutp << kplane << popout << fluxp << fixup << std::endl;
    std::cout << std::endl;
    std::cout << npey << npez << nthreads << std::endl;
    std::cout << thread_single << thread_funneled << thread_serialized << thread_multiple << std::endl;
    std::cout << thread_level << std::endl;
    if (do_nested) {
        std::cout << nnested << std::endl;
    } else {
        std::cout << nnested << std::endl;
    };
    std::cout << pce << std::endl;
    std::cout <<      /* FIXME: implied do loop */  << std::endl;
}

void setup_mat(int &i1, int &i2, int &j1, int &j2, int &k1, int &k2)
{
    int i;
    int j;
    int jj;
    int k;
    int kk;
    mat = 1;
    i1 = 1;
    i2 = 1;
    j1 = 1;
    j2 = 1;
    k1 = 1;
    k2 = 1;
    if (mat_opt == 1) {
        i1 = nx/4 + 1;
        i2 = 3*nx/4;
        if (ndimen > 1) {
            j1 = ny_gl/4 + 1;
            j2 = 3*ny_gl/4;
            if (ndimen > 2) {
                k1 = nz_gl/4 + 1;
                k2 = 3*nz_gl/4;
            };
        };
    } else {
        if (mat_opt == 2) {
            i2 = nx/2;
            if (ndimen > 1) {
                j2 = ny_gl/2;
                if (ndimen > 2) {
                    k2 = nz_gl/2;
                };
            };
        };
    };
    if (mat_opt > 0) {
        for (k=k1; k<=k2; k++) {
            if ((klb <= k) && (k <= kub)) {
                for (j=j1; j<=j2; j++) {
                    if ((jlb <= j) && (j <= jub)) {
                        for (i=i1; i<=i2; i++) {
                            mat[i,jj,kk-1] = 2;
                        };
                    };
                };
            };
        };
    };
}

void setup_scatp(int &ierr, std::string &error)
{
    std::string fn="slgg";
    int fu;
    int g1;
    int g2;
    int l;
    int n;
    ierr = 0;
    error = " ";
    open_file(fu, fn, "REPLACE", "WRITE", ierr, error);
    if (ierr != 0) {
        return;
    };
    std::cout << std::endl;
    std::cout <<      /* FIXME: implied do loop */  << std::endl;
    close_file(fu, ierr, error);
}

void setup_src(int &i1, int &i2, int &j1, int &j2, int &k1, int &k2, int &ierr, std::string &error)
{
    int i;
    int j;
    int jj;
    int k;
    int kk;
    if (src_opt == 3) {
        mms_setup(ierr, error);
        return;
    };
    i1 = 1;
    i2 = nx;
    j1 = 1;
    j2 = ny_gl;
    k1 = 1;
    k2 = nz_gl;
    if (src_opt == 1) {
        i1 = nx/4 + 1;
        i2 = 3*nx/4;
        if (ndimen > 1) {
            j1 = ny_gl/4 + 1;
            j2 = 3*ny_gl/4;
            if (ndimen > 2) {
                k1 = nz_gl/4 + 1;
                k2 = 3*nz_gl/4;
            };
        };
    } else {
        if (src_opt == 2) {
            i2 = nx/2;
            if (ndimen > 1) {
                j2 = ny_gl/2;
                if (ndimen > 2) {
                    k2 = nz_gl/2;
                };
            };
        };
    };
    for (k=k1; k<=k2; k++) {
        if ((klb <= k) && (k <= kub)) {
            for (j=j1; j<=j2; j++) {
                if ((jlb <= j) && (j <= jub)) {
                    for (i=i1; i<=i2; i++) {
                        qi[i,jj,kk,/* FIXME right index */-1] = one;
                    };
                };
            };
        };
    };
}

void setup_vel()
{
    int g;
    int t;
    if (timedep == 0) {
        return;
    };
    for (g=1; g<=ng; g++) {
        t = ng - g + 1;
        v[g-1] = t;
    };
}

void version_print()
{
    Kokkos::View<int[8]> dttm("dttm");
    Kokkos::View<std::string[3]> rc("rc");
}

namespace {

void main2() {
        std::string error;
        int i;
        int ierr;
        int ndpwds;
        std::string star="*";
        float t1;
        float t2;
        float t3;
        float t4;
        float t5;
        ierr = 0;
        ndpwds = 0;
        error = " ";
        pinit(t1);
        cmdarg(ierr, error);
        bcast_i_scalar(ierr, comm_snap, root);
        if (ierr != 0) {
            print_error(0, error);
            stop_run(0, 0, 0, 0);
        };
        open_file(iunit, ifile, "OLD", "READ", ierr, error);
        bcast_i_scalar(ierr, comm_snap, root);
        if (ierr != 0) {
            print_error(0, error);
            stop_run(0, 0, 0, 0);
        };
        open_file(ounit, ofile, "REPLACE", "WRITE", ierr, error);
        bcast_i_scalar(ierr, comm_snap, root);
        if (ierr != 0) {
            print_error(0, error);
            stop_run(0, 0, 0, 0);
        };
        if (iproc == root) {
            version_print();
        };
        input_read();
        close_file(iunit, ierr, error);
        bcast_i_scalar(ierr, comm_snap, root);
        if (ierr != 0) {
            print_error(ounit, error);
            stop_run(0, 0, 0, 0);
        };
        wtime(t3);
        pinit_omp(ierr, error);
        if (ierr != 0) {
            print_error(0, error);
        };
        pcomm_set();
        wtime(t4);
        tparset = tparset + t4 - t3;
        setup(ndpwds);
        output();
        if (iproc == root) {
            time_summ();
        };
        dealloc_input(4);
        dealloc_solve(swp_typ, 2);
        wtime(t5);
        tsnap = t5 - t1;
        if (iproc == root) {
            std::cout << tsnap << std::endl;
            std::cout << tgrind <<              /* FIXME: implied do loop */  << std::endl;
            std::cout << ndpwds <<              /* FIXME: implied do loop */  << std::endl;
        };
        close_file(ounit, ierr, error);
        bcast_i_scalar(ierr, comm_snap, root);
        if (ierr != 0) {
            print_error(0, error);
            stop_run(1, 0, 0, 0);
        };
        stop_run(1, 0, 0, 1);
}

}
int main(int argc, char* argv[])
{
    Kokkos::initialize(argc, argv);
    main2();
    Kokkos::finalize();
    return 0;
}