//
//  Created by Daisuke Imoto on 2015/07/17.
//  Copyright (c) 2015-Daisuke-Imoto All rights reserved.
//
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <sys/time.h>
#include "MT.h"

const int N_cell = 1;
const int size_x = 1600;
const int size_y = 1600;
const int L_comp = 350;
const int n_comp = L_comp*L_comp;
double Cphi_x, Cphi_y, me;

int main(){
    FILE *fp1;
    char filename[500];
    
    init_genrand((unsigned)time(NULL));
    
    double dx = 0.1;
    double dt = 0.00004;
    double dt0 = 0.00004;
    const int DDtt = 5000;
    const int Tmax = DDtt * 500;
    
    double** phi = new double*[N_cell];
    for (int ii = 0; ii < N_cell; ii++){
        phi[ii] = new double[n_comp];
    }
    double** phiv = new double*[N_cell];
    for (int ii = 0; ii < N_cell; ii++){
        phiv[ii] = new double[n_comp];
    }
    double** phiu = new double*[N_cell];
    for (int ii = 0; ii < N_cell; ii++){
        phiu[ii] = new double[n_comp];
    }
    double** phiz = new double*[N_cell];
    for (int ii = 0; ii < N_cell; ii++){
        phiz[ii] = new double[n_comp];
    }
    double** kai_grid = new double*[size_y];
    for (int iiy = 0; iiy < size_y; iiy++){
        kai_grid[iiy] = new double[size_x];
    }
    
    int cx_grid, cy_grid;
    cx_grid = L_comp / 2 - 1; cy_grid = L_comp / 2 - 1;
    
    int cx[N_cell], cy[N_cell];
    double radius = 5;
    cx[0] = 800; cy[0] = 800;
    //cx[1] = 800 + 50; cy[1] = 800;
    /*for (int i1 = 0; i1 < 4; i1++){
        for (int i2 = 0; i2 < 4; i2++){
            cx[4*i1 + i2] = 650 + 100*i1; cy[4*i1 + i2] = 650 + 100*i2;
        }
    }
    for (int i1 = 0; i1 < 7; i1++){
        for (int i2 = 0; i2 < 7; i2++){
            cx[7*i1 + i2] = 500 + 100*i1; cy[7*i1 + i2] = 500 + 100*i2;
        }
    }*/
    
    int Poli_x[N_cell], Poli_y[N_cell];
    Poli_x[0] = 1; Poli_y[0] = 0;
    
    for (int ii = 0; ii < N_cell; ii++){
        for (int iiy = 0; iiy < L_comp; iiy++){
            for (int iix = 0; iix < L_comp; iix++){
                if ((iix - cx_grid)*(iix - cx_grid) + (iiy - cy_grid)*(iiy - cy_grid) < (radius / dx)*(radius / dx)){
                    phi[ii][L_comp*iiy + iix] = 1.0;
                    phiu[ii][L_comp*iiy + iix] = 0.1;
                    phiv[ii][L_comp*iiy + iix] = 0.1;
                    if ((iix - cx_grid)*Poli_x[ii] + (iiy - cy_grid)*Poli_y[ii] > 0) {
                        phiz[ii][L_comp*iiy + iix] = 1.0;
                    }else{
                        phiz[ii][L_comp*iiy + iix] = 0.01;
                    }
                }
                else{
                    phi[ii][L_comp*iiy + iix] = 0.0;
                    phiu[ii][L_comp*iiy + iix] = 0.0;
                    phiv[ii][L_comp*iiy + iix] = 0.0;
                    phiz[ii][L_comp*iiy + iix] = 0.0;
                }
            }
        }
    }
    /*for (int ii = 0; ii < N_cell; ii++){
        sprintf(filename, "data/result_cell%d_t%d_phi.txt", ii, 0);
        fp1 = fopen(filename, "w");
        for (int yy = 0; yy < L_comp; yy++){
            for (int xx = 0; xx < L_comp; xx++){
                fprintf(fp1, "%f,", phi[ii][L_comp*yy + xx]);
            }
            fprintf(fp1, "\n");
        }
        fclose(fp1);
        
        sprintf(filename, "data/result_cell%d_t%d_phiu.txt", ii, 0);
        fp1 = fopen(filename, "w");
        for (int yy = 0; yy < L_comp; yy++){
            for (int xx = 0; xx < L_comp; xx++){
                fprintf(fp1, "%f,", phiu[ii][L_comp*yy + xx]);
            }
            fprintf(fp1, "\n");
        }
        fclose(fp1);
        
        sprintf(filename, "data/result_cell%d_t%d_phiv.txt", ii, 0);
        fp1 = fopen(filename, "w");
        for (int yy = 0; yy < L_comp; yy++){
            for (int xx = 0; xx < L_comp; xx++){
                fprintf(fp1, "%f,", phiv[ii][L_comp*yy + xx]);
            }
            fprintf(fp1, "\n");
        }
        fclose(fp1);
        
    }*/
    
    double epsilon_grid = 1.0;
    double ratio_grid = 0.4673;
    for (int iiy = 0; iiy < size_y; iiy++){
        for (int iix = 0; iix < size_x; iix++){
            kai_grid[iiy][iix] = 0.5*(1.0 + tanh((-sqrt(((double)iix - (double)size_x / 2.0)*((double)iix - (double)size_x / 2.0) + ((double)iiy - (double)size_y / 2.0)*((double)iiy - (double)size_y / 2.0)) + (double)size_x*ratio_grid) / (epsilon_grid / dx)));
        }
    }
    /*fp1 = fopen("data/kai_grid.txt", "w");
    for (int iiy = 0; iiy < size_y; iiy++){
        for (int iix = 0; iix < size_x; iix++){
            fprintf(fp1, "%f,", kai_grid[iiy][iix]);
        }
        fprintf(fp1, "\n");
    }
    fclose(fp1);*/
    
    double Du, Dv, alpha, Kk, beta, Kp, s, gamma, mu, kaiU, kaiV;
    double theta, d, sigma;
    Du = 0.05; Dv = 0.20; alpha = 2.0; beta = 3.5; Kk = 3.5; Kp = 3.2; kaiU = 50.0; kaiV = 0.0;
    s = 1.0; gamma = 0.1; mu = 0.1;
    theta = 1.4; d = 0.8; sigma = 0.075;
    double k_vz;
    k_vz = 200.0;
    double Dz, gw3, kw2, gw1, Ntot;
    Dz = 0.3*10.0; gw3 = 500.0; kw2 = 500.0; gw1 = 100.0; Ntot = 200.0*0.4;
    double tau, eta, epsilon, M, av, bu, aw, A_0;
    double lambda, kappa;
    tau = 0.83; eta = 1.0; epsilon = 1.0; M = 0.5; aw = 4.0*0.6; A_0 = 78.83;
    lambda = 10.0; kappa = 2.0/(dx*dx);
    
    double fac11, fac12, fac21, fac22, fac23;
    fac11 = expf(- 1.0 / (2.0*(d / dx)*(d / dx)));
    fac12 = expf(- 2.0 / (2.0*(d / dx)*(d / dx)));
    fac21 = expf(- 4.0 / (2.0*(d / dx)*(d / dx)));
    fac22 = expf(- 5.0 / (2.0*(d / dx)*(d / dx)));
    fac23 = expf(- 8.0 / (2.0*(d / dx)*(d / dx)));
    
    double** Lap_u = new double*[N_cell];
    for (int ii = 0; ii < N_cell; ii++){
        Lap_u[ii] = new double[n_comp];
    }
    double** Lap_v = new double*[N_cell];
    for (int ii = 0; ii < N_cell; ii++){
        Lap_v[ii] = new double[n_comp];
    }
    double** Lap_z = new double*[N_cell];
    for (int ii = 0; ii < N_cell; ii++){
        Lap_z[ii] = new double[n_comp];
    }
    double** Lap_phi = new double*[N_cell];
    for (int ii = 0; ii < N_cell; ii++){
        Lap_phi[ii] = new double[n_comp];
    }
    double** Nabla_phi_square = new double*[N_cell];
    for (int ii = 0; ii < N_cell; ii++){
        Nabla_phi_square[ii] = new double[n_comp];
    }
    double** Nabla_phi_abs = new double*[N_cell];
    for (int ii = 0; ii < N_cell; ii++){
        Nabla_phi_abs[ii] = new double[n_comp];
    }
    double** diff_phiu = new double*[N_cell];
    for (int ii = 0; ii < N_cell; ii++){
        diff_phiu[ii] = new double[n_comp];
    }
    double** diff_phiv = new double*[N_cell];
    for (int ii = 0; ii < N_cell; ii++){
        diff_phiv[ii] = new double[n_comp];
    }
    double** diff_phiz = new double*[N_cell];
    for (int ii = 0; ii < N_cell; ii++){
        diff_phiz[ii] = new double[n_comp];
    }
    double** dif_phi0 = new double*[N_cell];
    for (int ii = 0; ii < N_cell; ii++){
        dif_phi0[ii] = new double[n_comp];
    }
    double** NoiseVadd = new double*[N_cell];
    for (int ii = 0; ii < N_cell; ii++){
        NoiseVadd[ii] = new double[n_comp];
    }
    for (int ii = 0; ii < N_cell; ii++){
        for (int iiy = 0; iiy < L_comp; iiy++){
            for (int iix = 0; iix < L_comp; iix++){
                NoiseVadd[ii][L_comp*iiy + iix] = 0.0;
            }
        }
    }
    
    double* phi_t = new double[n_comp];
    double* phiv_t = new double[n_comp];
    double* phiu_t = new double[n_comp];
    double* phiz_t = new double[n_comp];
    
    double sum;
    
    double sumphi[N_cell], sumz[N_cell], sumv2[N_cell], sumu[N_cell], sumNablaphi2[N_cell];
    
    double lambda_phi = 1e-4;
    
    double avephixp1, avephixm1, avephiyp1, avephiym1;
    
    //
    double p_one = theta*dt*dx*dx;
    int n_max = 10;
    int n_tmp[N_cell];
    int iiy_noise_tmp;
    int iix_noise_tmp;
    double st_noise_tmp;
    double prob_s[n_max+1];
    prob_s[0] = pow(1-p_one,n_comp);
    for (int i = 1; i <= n_max; i++){
        double nci = 1.0;
        for (int j = 0; j < i; j++){
            nci *= (n_comp-j)/(j+1);
        }
        prob_s[i] = prob_s[i - 1] + nci*pow(p_one,i)*pow(1-p_one,n_comp-i);
    }
    
    printf( "tic\n" );
    //開始時刻
    struct timeval start_time;
    gettimeofday(&start_time, NULL);
    
    int Tmax0=(int)(Tmax/DDtt);
    //int Tmax0 = 1;
    for (int tt = 0; tt < Tmax0; tt++){
        
        // 外場のプロファイル
#pragma acc data copy(phi[0:N_cell][0:n_comp],phiv[0:N_cell][0:n_comp],phiu[0:N_cell][0:n_comp],phiz[0:N_cell][0:n_comp]) create(Lap_phi[0:N_cell][0:n_comp], Lap_u[0:N_cell][0:n_comp], Lap_v[0:N_cell][0:n_comp],Lap_z[0:N_cell][0:n_comp], Nabla_phi_square[0:N_cell][0:n_comp], Nabla_phi_abs[0:N_cell][0:n_comp],diff_phiu[0:N_cell][0:n_comp],diff_phiv[0:N_cell][0:n_comp],diff_phiz[0:N_cell][0:n_comp],dif_phi0[0:N_cell][0:n_comp],phi_t[0:n_comp],phiv_t[0:n_comp],phiu_t[0:n_comp],phiz_t[0:n_comp]) copyin(kai_grid[0:size_y][0:size_x])
        for (int t2 = 0; t2 < DDtt; t2++){
            
            for (int ii = 0; ii < N_cell; ii++){
                double p_tmp = genrand_real3();
                for (int i = 0; i <= n_max; i++){
                    n_tmp[ii] = i;
                    if (p_tmp < prob_s[i]) {
                        break;
                    }
                }
                if (n_tmp[ii] > 0) {
                    for (int iiy = 0; iiy < L_comp; iiy++){
                        for (int iix = 0; iix < L_comp; iix++){
                            NoiseVadd[ii][L_comp*iiy + iix] = 0.0;
                        }
                    }
                    
                    for (int i = 0; i < n_tmp[ii]; i++){
                        iiy_noise_tmp = genrand_int32()%L_comp;
                        iix_noise_tmp = genrand_int32()%L_comp;
                        double ftmp = genrand_real3();
                        st_noise_tmp = -log(ftmp)*sigma;
                        //printf("(t2,ii)=%d,%d\n",t2,ii);
                        //printf("(y,x)=%d,%d\n",iiy_noise_tmp,iix_noise_tmp);
                        //printf("s=%f\n",st_noise_tmp);
                        
                        if ((iiy_noise_tmp>=2&&iiy_noise_tmp<L_comp-2)&&(iix_noise_tmp>=2&&iix_noise_tmp<L_comp-2)) {
                            for (int yy = - 2; yy < 2; yy++){
                                for (int xx = - 2; xx < 2; xx++){
                                    int iiy = iiy_noise_tmp + yy;
                                    int iix = iix_noise_tmp + xx;
                                    if (xx*xx+yy*yy==0) { NoiseVadd[ii][L_comp*iiy + iix] += st_noise_tmp; }
                                    if (xx*xx+yy*yy==1) { NoiseVadd[ii][L_comp*iiy + iix] += fac11*st_noise_tmp; }
                                    if (xx*xx+yy*yy==2) { NoiseVadd[ii][L_comp*iiy + iix] += fac12*st_noise_tmp; }
                                    if (xx*xx+yy*yy==4) { NoiseVadd[ii][L_comp*iiy + iix] += fac21*st_noise_tmp; }
                                    if (xx*xx+yy*yy==5) { NoiseVadd[ii][L_comp*iiy + iix] += fac22*st_noise_tmp; }
                                    if (xx*xx+yy*yy==8) { NoiseVadd[ii][L_comp*iiy + iix] += fac23*st_noise_tmp; }
                                }
                            }
                        }
                        
                    }
                    
                }
            }
            
            for (int ii = 0; ii < N_cell; ii++){
                sum = 0.0;
#pragma acc parallel loop reduction(+:sum)
                for (int rr = 0; rr < n_comp; rr++){
                    sum += phi[ii][rr];
                }
                sumphi[ii] = sum;
                
                sum = 0.0;
#pragma acc parallel loop reduction(+:sum)
                for (int rr = 0; rr < n_comp; rr++){
                    sum += phi[ii][rr]*phiz[ii][rr];
                }
                sumz[ii] = sum;
                
                sum = 0.0;
#pragma acc parallel loop reduction(+:sum)
                for (int rr = 0; rr < n_comp; rr++){
                    sum += phi[ii][rr]*phiu[ii][rr];
                }
                sumu[ii] = sum;
                
                sum = 0.0;
#pragma acc parallel loop reduction(+:sum)
                for (int rr = 0; rr < n_comp; rr++){
                    sum += phi[ii][rr]*phiz[ii][rr]*phiv[ii][rr];
                }
                sumv2[ii] = sum;
            }
            
            for (int ii = 0; ii < N_cell; ii++){
                
#pragma acc parallel loop
                for (int rr = 0; rr < n_comp; rr++){
                    int iiy = (int)(rr/L_comp);
                    int iix = rr - L_comp*iiy;
                    
                    if (iiy > 0 && iiy<L_comp - 1 && iix>0 && iix < L_comp - 1) {
                        avephixp1 = 0.5*(phi[ii][L_comp*iiy + iix + 1] + phi[ii][L_comp*iiy + iix]);
                        avephixm1 = 0.5*(phi[ii][L_comp*iiy + iix - 1] + phi[ii][L_comp*iiy + iix]);
                        avephiyp1 = 0.5*(phi[ii][L_comp*(iiy + 1) + iix] + phi[ii][L_comp*iiy + iix]);
                        avephiym1 = 0.5*(phi[ii][L_comp*(iiy - 1) + iix] + phi[ii][L_comp*iiy + iix]);
                        Lap_u[ii][L_comp*iiy + iix] = (avephixp1*(phiu[ii][L_comp*iiy + iix + 1] - phiu[ii][L_comp*iiy + iix]) - avephixm1*(phiu[ii][L_comp*iiy + iix] - phiu[ii][L_comp*iiy + iix - 1])) + (avephiyp1*(phiu[ii][L_comp*(iiy + 1) + iix] - phiu[ii][L_comp*iiy + iix]) - avephiym1*(phiu[ii][L_comp*iiy + iix] - phiu[ii][L_comp*(iiy - 1) + iix]));
                        Lap_v[ii][L_comp*iiy + iix] = (avephixp1*(phiv[ii][L_comp*iiy + iix + 1] - phiv[ii][L_comp*iiy + iix]) - avephixm1*(phiv[ii][L_comp*iiy + iix] - phiv[ii][L_comp*iiy + iix - 1])) + (avephiyp1*(phiv[ii][L_comp*(iiy + 1) + iix] - phiv[ii][L_comp*iiy + iix]) - avephiym1*(phiv[ii][L_comp*iiy + iix] - phiv[ii][L_comp*(iiy - 1) + iix]));
                        Lap_z[ii][L_comp*iiy + iix] = (avephixp1*(phiz[ii][L_comp*iiy + iix + 1] - phiz[ii][L_comp*iiy + iix]) - avephixm1*(phiz[ii][L_comp*iiy + iix] - phiz[ii][L_comp*iiy + iix - 1])) + (avephiyp1*(phiz[ii][L_comp*(iiy + 1) + iix] - phiz[ii][L_comp*iiy + iix]) - avephiym1*(phiz[ii][L_comp*iiy + iix] - phiz[ii][L_comp*(iiy - 1) + iix]));
                        Lap_phi[ii][L_comp*iiy + iix] = -4.0 * phi[ii][L_comp*iiy + iix] + phi[ii][L_comp*(iiy + 1) + iix] + phi[ii][L_comp*(iiy - 1) + iix] + phi[ii][L_comp*iiy + iix + 1] + phi[ii][L_comp*iiy + iix - 1];
                        Nabla_phi_square[ii][L_comp*iiy + iix] = 0.5*(phi[ii][L_comp*(iiy + 1) + iix] - phi[ii][L_comp*(iiy - 1) + iix])*0.5*(phi[ii][L_comp*(iiy + 1) + iix] - phi[ii][L_comp*(iiy - 1) + iix]) + 0.5*(phi[ii][L_comp*iiy + iix + 1] - phi[ii][L_comp*iiy + iix - 1])*0.5*(phi[ii][L_comp*iiy + iix + 1] - phi[ii][L_comp*iiy + iix - 1]);
                        if (Nabla_phi_square[ii][L_comp*iiy + iix] < 0.0){
                            Nabla_phi_square[ii][L_comp*iiy + iix] = 0.0;
                        }
                        Nabla_phi_abs[ii][L_comp*iiy + iix] = sqrt(Nabla_phi_square[ii][L_comp*iiy + iix]);
                    }
                    if (iiy == 0 && iix>0 && iix < L_comp - 1) {
                        avephixp1 = 0.5*(phi[ii][L_comp*iiy + iix + 1] + phi[ii][L_comp*iiy + iix]);
                        avephixm1 = 0.5*(phi[ii][L_comp*iiy + iix - 1] + phi[ii][L_comp*iiy + iix]);
                        avephiyp1 = 0.5*(phi[ii][L_comp*(iiy + 1) + iix] + phi[ii][L_comp*iiy + iix]);
                        Lap_u[ii][L_comp*iiy + iix] = (avephixp1*(phiu[ii][L_comp*iiy + iix + 1] - phiu[ii][L_comp*iiy + iix]) - avephixm1*(phiu[ii][L_comp*iiy + iix] - phiu[ii][L_comp*iiy + iix - 1])) + (avephiyp1*(phiu[ii][L_comp*(iiy + 1) + iix] - phiu[ii][L_comp*iiy + iix]));
                        Lap_v[ii][L_comp*iiy + iix] = (avephixp1*(phiv[ii][L_comp*iiy + iix + 1] - phiv[ii][L_comp*iiy + iix]) - avephixm1*(phiv[ii][L_comp*iiy + iix] - phiv[ii][L_comp*iiy + iix - 1])) + (avephiyp1*(phiv[ii][L_comp*(iiy + 1) + iix] - phiv[ii][L_comp*iiy + iix]));
                        Lap_z[ii][L_comp*iiy + iix] = (avephixp1*(phiz[ii][L_comp*iiy + iix + 1] - phiz[ii][L_comp*iiy + iix]) - avephixm1*(phiz[ii][L_comp*iiy + iix] - phiz[ii][L_comp*iiy + iix - 1])) + (avephiyp1*(phiz[ii][L_comp*(iiy + 1) + iix] - phiz[ii][L_comp*iiy + iix]));
                        Lap_phi[ii][L_comp*iiy + iix] = -4.0 * phi[ii][L_comp*iiy + iix] + phi[ii][L_comp*(iiy + 1) + iix] + phi[ii][L_comp*iiy + iix] + phi[ii][L_comp*iiy + iix + 1] + phi[ii][L_comp*iiy + iix - 1];
                        Nabla_phi_square[ii][L_comp*iiy + iix] = 0.5*(phi[ii][L_comp*(iiy + 1) + iix] - phi[ii][L_comp*iiy + iix])*0.5*(phi[ii][L_comp*(iiy + 1) + iix] - phi[ii][L_comp*iiy + iix]) + 0.5*(phi[ii][L_comp*iiy + iix + 1] - phi[ii][L_comp*iiy + iix - 1])*0.5*(phi[ii][L_comp*iiy + iix + 1] - phi[ii][L_comp*iiy + iix - 1]);
                        if (Nabla_phi_square[ii][L_comp*iiy + iix] < 0.0){
                            Nabla_phi_square[ii][L_comp*iiy + iix] = 0.0;
                        }
                        Nabla_phi_abs[ii][L_comp*iiy + iix] = sqrt(Nabla_phi_square[ii][L_comp*iiy + iix]);
                    }
                    if (iiy == L_comp - 1 && iix>0 && iix < L_comp - 1) {
                        avephixp1 = 0.5*(phi[ii][L_comp*iiy + iix + 1] + phi[ii][L_comp*iiy + iix]);
                        avephixm1 = 0.5*(phi[ii][L_comp*iiy + iix - 1] + phi[ii][L_comp*iiy + iix]);
                        avephiym1 = 0.5*(phi[ii][L_comp*(iiy - 1) + iix] + phi[ii][L_comp*iiy + iix]);
                        Lap_u[ii][L_comp*iiy + iix] = (avephixp1*(phiu[ii][L_comp*iiy + iix + 1] - phiu[ii][L_comp*iiy + iix]) - avephixm1*(phiu[ii][L_comp*iiy + iix] - phiu[ii][L_comp*iiy + iix - 1])) - avephiym1*(phiu[ii][L_comp*iiy + iix] - phiu[ii][L_comp*(iiy - 1) + iix]);
                        Lap_v[ii][L_comp*iiy + iix] = (avephixp1*(phiv[ii][L_comp*iiy + iix + 1] - phiv[ii][L_comp*iiy + iix]) - avephixm1*(phiv[ii][L_comp*iiy + iix] - phiv[ii][L_comp*iiy + iix - 1])) - avephiym1*(phiv[ii][L_comp*iiy + iix] - phiv[ii][L_comp*(iiy - 1) + iix]);
                        Lap_z[ii][L_comp*iiy + iix] = (avephixp1*(phiz[ii][L_comp*iiy + iix + 1] - phiz[ii][L_comp*iiy + iix]) - avephixm1*(phiz[ii][L_comp*iiy + iix] - phiz[ii][L_comp*iiy + iix - 1])) - avephiym1*(phiz[ii][L_comp*iiy + iix] - phiz[ii][L_comp*(iiy - 1) + iix]);
                        Lap_phi[ii][L_comp*iiy + iix] = -4.0 * phi[ii][L_comp*iiy + iix] + phi[ii][L_comp*iiy + iix] + phi[ii][L_comp*(iiy - 1) + iix] + phi[ii][L_comp*iiy + iix + 1] + phi[ii][L_comp*iiy + iix - 1];
                        Nabla_phi_square[ii][L_comp*iiy + iix] = 0.5*(phi[ii][L_comp*iiy + iix] - phi[ii][L_comp*(iiy - 1) + iix])*0.5*(phi[ii][L_comp*iiy + iix] - phi[ii][L_comp*(iiy - 1) + iix]) + 0.5*(phi[ii][L_comp*iiy + iix + 1] - phi[ii][L_comp*iiy + iix - 1])*0.5*(phi[ii][L_comp*iiy + iix + 1] - phi[ii][L_comp*iiy + iix - 1]);
                        if (Nabla_phi_square[ii][L_comp*iiy + iix] < 0.0){
                            Nabla_phi_square[ii][L_comp*iiy + iix] = 0.0;
                        }
                        Nabla_phi_abs[ii][L_comp*iiy + iix] = sqrt(Nabla_phi_square[ii][L_comp*iiy + iix]);
                    }
                    if (iiy>0 && iiy < L_comp - 1 && iix == 0) {
                        avephixp1 = 0.5*(phi[ii][L_comp*iiy + iix + 1] + phi[ii][L_comp*iiy + iix]);
                        avephiyp1 = 0.5*(phi[ii][L_comp*(iiy + 1) + iix] + phi[ii][L_comp*iiy + iix]);
                        avephiym1 = 0.5*(phi[ii][L_comp*(iiy - 1) + iix] + phi[ii][L_comp*iiy + iix]);
                        Lap_u[ii][L_comp*iiy + iix] = (avephixp1*(phiu[ii][L_comp*iiy + iix + 1] - phiu[ii][L_comp*iiy + iix])) + (avephiyp1*(phiu[ii][L_comp*(iiy + 1) + iix] - phiu[ii][L_comp*iiy + iix]) - avephiym1*(phiu[ii][L_comp*iiy + iix] - phiu[ii][L_comp*(iiy - 1) + iix]));
                        Lap_v[ii][L_comp*iiy + iix] = (avephixp1*(phiv[ii][L_comp*iiy + iix + 1] - phiv[ii][L_comp*iiy + iix])) + (avephiyp1*(phiv[ii][L_comp*(iiy + 1) + iix] - phiv[ii][L_comp*iiy + iix]) - avephiym1*(phiv[ii][L_comp*iiy + iix] - phiv[ii][L_comp*(iiy - 1) + iix]));
                        Lap_z[ii][L_comp*iiy + iix] = (avephixp1*(phiz[ii][L_comp*iiy + iix + 1] - phiz[ii][L_comp*iiy + iix])) + (avephiyp1*(phiz[ii][L_comp*(iiy + 1) + iix] - phiz[ii][L_comp*iiy + iix]) - avephiym1*(phiz[ii][L_comp*iiy + iix] - phiz[ii][L_comp*(iiy - 1) + iix]));
                        Lap_phi[ii][L_comp*iiy + iix] = -4.0 * phi[ii][L_comp*iiy + iix] + phi[ii][L_comp*(iiy + 1) + iix] + phi[ii][L_comp*(iiy - 1) + iix] + phi[ii][L_comp*iiy + iix + 1] + phi[ii][L_comp*iiy + iix];
                        Nabla_phi_square[ii][L_comp*iiy + iix] = 0.5*(phi[ii][L_comp*(iiy + 1) + iix] - phi[ii][L_comp*(iiy - 1) + iix])*0.5*(phi[ii][L_comp*(iiy + 1) + iix] - phi[ii][L_comp*(iiy - 1) + iix]) + 0.5*(phi[ii][L_comp*iiy + iix + 1] - phi[ii][L_comp*iiy + iix])*0.5*(phi[ii][L_comp*iiy + iix + 1] - phi[ii][L_comp*iiy + iix]);
                        if (Nabla_phi_square[ii][L_comp*iiy + iix] < 0.0){
                            Nabla_phi_square[ii][L_comp*iiy + iix] = 0.0;
                        }
                        Nabla_phi_abs[ii][L_comp*iiy + iix] = sqrt(Nabla_phi_square[ii][L_comp*iiy + iix]);
                    }
                    if (iiy>0 && iiy < L_comp - 1 && iix == L_comp - 1) {
                        avephixm1 = 0.5*(phi[ii][L_comp*iiy + iix - 1] + phi[ii][L_comp*iiy + iix]);
                        avephiyp1 = 0.5*(phi[ii][L_comp*(iiy + 1) + iix] + phi[ii][L_comp*iiy + iix]);
                        avephiym1 = 0.5*(phi[ii][L_comp*(iiy - 1) + iix] + phi[ii][L_comp*iiy + iix]);
                        Lap_u[ii][L_comp*iiy + iix] = (-avephixm1*(phiu[ii][L_comp*iiy + iix] - phiu[ii][L_comp*iiy + iix - 1])) + (avephiyp1*(phiu[ii][L_comp*(iiy + 1) + iix] - phiu[ii][L_comp*iiy + iix]) - avephiym1*(phiu[ii][L_comp*iiy + iix] - phiu[ii][L_comp*(iiy - 1) + iix]));
                        Lap_v[ii][L_comp*iiy + iix] = (-avephixm1*(phiv[ii][L_comp*iiy + iix] - phiv[ii][L_comp*iiy + iix - 1])) + (avephiyp1*(phiv[ii][L_comp*(iiy + 1) + iix] - phiv[ii][L_comp*iiy + iix]) - avephiym1*(phiv[ii][L_comp*iiy + iix] - phiv[ii][L_comp*(iiy - 1) + iix]));
                        Lap_z[ii][L_comp*iiy + iix] = (-avephixm1*(phiz[ii][L_comp*iiy + iix] - phiz[ii][L_comp*iiy + iix - 1])) + (avephiyp1*(phiz[ii][L_comp*(iiy + 1) + iix] - phiz[ii][L_comp*iiy + iix]) - avephiym1*(phiz[ii][L_comp*iiy + iix] - phiz[ii][L_comp*(iiy - 1) + iix]));
                        Lap_phi[ii][L_comp*iiy + iix] = -4.0 * phi[ii][L_comp*iiy + iix] + phi[ii][L_comp*(iiy + 1) + iix] + phi[ii][L_comp*(iiy - 1) + iix] + phi[ii][L_comp*iiy + iix] + phi[ii][L_comp*iiy + iix - 1];
                        Nabla_phi_square[ii][L_comp*iiy + iix] = 0.5*(phi[ii][L_comp*(iiy + 1) + iix] - phi[ii][L_comp*(iiy - 1) + iix])*0.5*(phi[ii][L_comp*(iiy + 1) + iix] - phi[ii][L_comp*(iiy - 1) + iix]) + 0.5*(phi[ii][L_comp*iiy + iix] - phi[ii][L_comp*iiy + iix - 1])*0.5*(phi[ii][L_comp*iiy + iix] - phi[ii][L_comp*iiy + iix - 1]);
                        if (Nabla_phi_square[ii][L_comp*iiy + iix] < 0.0){
                            Nabla_phi_square[ii][L_comp*iiy + iix] = 0.0;
                        }
                        Nabla_phi_abs[ii][L_comp*iiy + iix] = sqrt(Nabla_phi_square[ii][L_comp*iiy + iix]);
                    }
                    if ((iiy == 0 || iiy == L_comp - 1) && (iix == 0 || iix == L_comp - 1)) {
                        Lap_u[ii][L_comp*iiy + iix] = 0.0;
                        Lap_v[ii][L_comp*iiy + iix] = 0.0;
                        Lap_z[ii][L_comp*iiy + iix] = 0.0;
                        Lap_phi[ii][L_comp*iiy + iix] = 0.0;
                        Nabla_phi_square[ii][L_comp*iiy + iix] = 0.0;
                        Nabla_phi_abs[ii][L_comp*iiy + iix] = 0.0;
                    }
                }
            }
            
            for (int ii = 0; ii < N_cell; ii++){
                sum = 0.0;
#pragma acc parallel loop reduction(+:sum)
                for (int rr = 0; rr < n_comp; rr++){
                    sum += Nabla_phi_square[ii][rr];
                }
                sumNablaphi2[ii] = sum;
            }
            
            for (int ii = 0; ii < N_cell; ii++){
                int Rx = (int)(cx[ii] - (L_comp / 2 - 1));
                int Ry = (int)(cy[ii] - (L_comp / 2 - 1));
                
#pragma acc parallel loop
                for (int rr = 0; rr < n_comp; rr++){
                    int iiy = (int)(rr/L_comp);
                    int iix = rr - L_comp*iiy;
                    if ((Rx + iix >= 0 && Rx + iix < size_x) && (Ry + iiy >= 0 && Ry + iiy < size_y)) {
                        diff_phiu[ii][rr] = (Du*Lap_u[ii][rr] / (dx*dx) + phi[ii][rr] * (-alpha*phiz[ii][rr]*phiu[ii][rr] * phiv[ii][rr] / (Kk + sumv2[ii] / sumphi[ii]) + beta*phiu[ii][rr] * phiv[ii][rr] / (Kp + sumu[ii] / sumphi[ii]) + s - gamma*phiu[ii][rr]) - kaiU*phiu[ii][rr]*Nabla_phi_square[ii][rr]/(sumNablaphi2[ii]*dx*dx));
                        diff_phiv[ii][rr] = (Dv*Lap_v[ii][rr] / (dx*dx) + phi[ii][rr] * (alpha*phiz[ii][rr]*phiu[ii][rr] * phiv[ii][rr] / (Kk + sumv2[ii] / sumphi[ii]) - beta*phiu[ii][rr] * phiv[ii][rr] / (Kp + sumu[ii] / sumphi[ii]) - mu*phiv[ii][rr]) - kaiV*phiv[ii][rr]*Nabla_phi_square[ii][rr]/(sumNablaphi2[ii]*dx*dx));
                        diff_phiz[ii][rr] = (Dz*Lap_z[ii][rr] / (dx*dx) + phi[ii][rr] * (- gw3*phiz[ii][rr]*phiz[ii][rr]*phiz[ii][rr] + kw2*(Ntot - sumz[ii]*dx*dx)/(sumphi[ii]*dx*dx)*phiz[ii][rr]*phiz[ii][rr] - gw1*phiz[ii][rr] + k_vz*phiv[ii][rr]));
                        dif_phi0[ii][rr] = (1.0 / tau)*(eta*(Lap_phi[ii][rr] / (dx*dx) - (1.0 / (epsilon*epsilon)) * 72.0 * phi[ii][rr] * (1.0 - phi[ii][rr])*(0.5 - phi[ii][rr])) + M*(A_0 - sumphi[ii] * dx*dx)*Nabla_phi_abs[ii][rr] / dx + (aw*phiz[ii][rr])*Nabla_phi_abs[ii][rr] / dx);
                    }
                }
                
#pragma acc parallel loop
                for (int rr = 0; rr < n_comp; rr++){
                    int iiy = (int)(rr/L_comp);
                    int iix = rr - L_comp*iiy;
                    if ((Rx + iix >= 0 && Rx + iix < size_x) && (Ry + iiy >= 0 && Ry + iiy < size_y)) {
                        dif_phi0[ii][rr] += -lambda*phi[ii][rr]*(1.0 - kai_grid[Ry + iiy][Rx + iix])*(1.0 - kai_grid[Ry + iiy][Rx + iix]);
                    }
                }
                
                for (int jj = 0; jj < N_cell; jj++){
                    if (ii != jj){
#pragma acc parallel loop
                        for (int rr = 0; rr < n_comp; rr++){
                            int iiy = (int)(rr/L_comp);
                            int iix = rr - L_comp*iiy;
                            if ((Rx + iix >= 0 && Rx + iix < size_x) && (Ry + iiy >= 0 && Ry + iiy < size_y)) {
                                int ixj = (int)(iix + cx[ii] - cx[jj]);
                                int iyj = (int)(iiy + cy[ii] - cy[jj]);
                                if (ixj >= 0 && ixj < L_comp && iyj >= 0 && iyj < L_comp){
                                    dif_phi0[ii][rr] += -lambda*phi[ii][rr]*phi[jj][L_comp*iyj + ixj]*phi[jj][L_comp*iyj + ixj];
                                    dif_phi0[ii][rr] += kappa * Nabla_phi_abs[ii][rr] * Nabla_phi_abs[jj][L_comp*iyj + ixj];
                                }
                            }
                        }
                    }
                    
                }
            }
            
#pragma acc data copyin(NoiseVadd[0:N_cell][0:n_comp])
            {
                
#pragma acc parallel loop
                for (int i = 0; i < N_cell*n_comp; i++){
                    int ii = (int)(i/n_comp);
                    int iiy = (int)((i - n_comp*ii)/L_comp);
                    int iix = i - n_comp*ii - L_comp*iiy;
                    int rr = L_comp*iiy + iix;
                    int Rx = (int)(cx[ii] - (L_comp / 2 - 1));
                    int Ry = (int)(cy[ii] - (L_comp / 2 - 1));
                    
                    if ((Rx + iix >= 0 && Rx + iix<size_x) && (Ry + iiy >= 0 && Ry + iiy<size_y)) {
                        if (phi[ii][rr] > lambda_phi) {
                            if (n_tmp[ii] > 0) {
                                if (NoiseVadd[ii][rr]>phiu[ii][rr]){
                                    phiu[ii][rr] = (phiu[ii][rr]*phi[ii][rr] + dt*diff_phiu[ii][rr])/phi[ii][rr] - phiu[ii][rr];
                                    phiv[ii][rr] = (phiv[ii][rr]*phi[ii][rr] + dt*diff_phiv[ii][rr])/phi[ii][rr] + phiu[ii][rr];
                                }else{
                                    phiu[ii][rr] = (phiu[ii][rr]*phi[ii][rr] + dt*diff_phiu[ii][rr])/phi[ii][rr] - NoiseVadd[ii][rr];
                                    phiv[ii][rr] = (phiv[ii][rr]*phi[ii][rr] + dt*diff_phiv[ii][rr])/phi[ii][rr] + NoiseVadd[ii][rr];
                                }
                            }else{
                                phiu[ii][rr] = (phiu[ii][rr]*phi[ii][rr] + dt*diff_phiu[ii][rr])/phi[ii][rr];
                                phiv[ii][rr] = (phiv[ii][rr]*phi[ii][rr] + dt*diff_phiv[ii][rr])/phi[ii][rr];
                            }
                            phiz[ii][rr] = (phiz[ii][rr]*phi[ii][rr] + dt*diff_phiz[ii][rr])/phi[ii][rr];
                        }else{
                            phiu[ii][rr] = (phiu[ii][rr]*phi[ii][rr] + dt*diff_phiu[ii][rr]);
                            phiv[ii][rr] = (phiv[ii][rr]*phi[ii][rr] + dt*diff_phiv[ii][rr]);
                            phiz[ii][rr] = (phiz[ii][rr]*phi[ii][rr] + dt*diff_phiz[ii][rr]);
                        }
                        phi[ii][rr] += (dt*dif_phi0[ii][rr]);
                    }
                }
            }
            
            
            if (t2 == DDtt - 1) {
                
                for (int ii = 0; ii < N_cell; ii++){
                    
                    Cphi_x = 0.0; Cphi_y = 0.0; me = 0.0;
                    
#pragma acc parallel loop
                    for (int iiy = 0; iiy < L_comp; iiy++){
#pragma acc loop
                        for (int iix = 0; iix < L_comp; iix++){
                            Cphi_x += (double)iix*phi[ii][L_comp*iiy + iix];
                        }
                    }
                    
#pragma acc parallel loop
                    for (int iiy = 0; iiy < L_comp; iiy++){
#pragma acc loop
                        for (int iix = 0; iix < L_comp; iix++){
                            Cphi_y += (double)iiy*phi[ii][L_comp*iiy + iix];
                        }
                    }
                    
#pragma acc parallel loop
                    for (int iiy = 0; iiy < L_comp; iiy++){
#pragma acc loop
                        for (int iix = 0; iix < L_comp; iix++){
                            me += phi[ii][L_comp*iiy + iix];
                        }
                    }
                    
                    int Cp_x, Cp_y, DifX, DifY;
                    Cp_x = (int)(Cphi_x / me);
                    Cp_y = (int)(Cphi_y / me);
                    DifX = Cp_x - (L_comp / 2 - 1);
                    DifY = Cp_y - (L_comp / 2 - 1);
                    
                    //printf("%d,%d\n",DifX,DifY);
                    
#pragma acc parallel loop
                    for (int iiy = 0; iiy < L_comp; iiy++){
#pragma acc loop
                        for (int iix = 0; iix < L_comp; iix++){
                            phi_t[L_comp*iiy + iix] = phi[ii][L_comp*iiy + iix];
                            phiz_t[L_comp*iiy + iix] = phiz[ii][L_comp*iiy + iix];
                            phiv_t[L_comp*iiy + iix] = phiv[ii][L_comp*iiy + iix];
                            phiu_t[L_comp*iiy + iix] = phiu[ii][L_comp*iiy + iix];
                            phi[ii][L_comp*iiy + iix] = 0.0;
                            phiz[ii][L_comp*iiy + iix] = 0.0;
                            phiv[ii][L_comp*iiy + iix] = 0.0;
                            phiu[ii][L_comp*iiy + iix] = 0.0;
                        }
                    }
                    
                    if (DifX > 0){
                        if (DifY > 0){
#pragma acc parallel loop
                            for (int iiy = 0; iiy < L_comp - DifY; iiy++){
#pragma acc loop
                                for (int iix = 0; iix < L_comp - DifX; iix++){
                                    phi[ii][L_comp*iiy + iix] = phi_t[L_comp*(iiy + DifY) + iix + DifX];
                                    phiv[ii][L_comp*iiy + iix] = phiv_t[L_comp*(iiy + DifY) + iix + DifX];
                                    phiu[ii][L_comp*iiy + iix] = phiu_t[L_comp*(iiy + DifY) + iix + DifX];
                                    phiz[ii][L_comp*iiy + iix] =phiz_t[L_comp*(iiy + DifY) + iix + DifX];
                                }
                            }
                        }
                        else{
#pragma acc parallel loop
                            for (int iiy = -DifY; iiy < L_comp; iiy++){
#pragma acc loop
                                for (int iix = 0; iix < L_comp - DifX; iix++){
                                    phi[ii][L_comp*iiy + iix] = phi_t[L_comp*(iiy + DifY) + iix + DifX];
                                    phiv[ii][L_comp*iiy + iix] = phiv_t[L_comp*(iiy + DifY) + iix + DifX];
                                    phiu[ii][L_comp*iiy + iix] = phiu_t[L_comp*(iiy + DifY) + iix + DifX];
                                    phiz[ii][L_comp*iiy + iix] =phiz_t[L_comp*(iiy + DifY) + iix + DifX];
                                }
                            }
                        }
                    }
                    else{
                        
                        if (DifY > 0){
#pragma acc parallel loop
                            for (int iiy = 0; iiy < L_comp - DifY; iiy++){
#pragma acc loop
                                for (int iix = -DifX; iix < L_comp; iix++){
                                    phi[ii][L_comp*iiy + iix] = phi_t[L_comp*(iiy + DifY) + iix + DifX];
                                    phiv[ii][L_comp*iiy + iix] = phiv_t[L_comp*(iiy + DifY) + iix + DifX];
                                    phiu[ii][L_comp*iiy + iix] = phiu_t[L_comp*(iiy + DifY) + iix + DifX];
                                    phiz[ii][L_comp*iiy + iix] =phiz_t[L_comp*(iiy + DifY) + iix + DifX];
                                }
                            }
                        }
                        else{
#pragma acc parallel loop
                            for (int iiy = -DifY; iiy < L_comp; iiy++){
#pragma acc loop
                                for (int iix = -DifX; iix < L_comp; iix++){
                                    phi[ii][L_comp*iiy + iix] = phi_t[L_comp*(iiy + DifY) + iix + DifX];
                                    phiv[ii][L_comp*iiy + iix] = phiv_t[L_comp*(iiy + DifY) + iix + DifX];
                                    phiu[ii][L_comp*iiy + iix] = phiu_t[L_comp*(iiy + DifY) + iix + DifX];
                                    phiz[ii][L_comp*iiy + iix] =phiz_t[L_comp*(iiy + DifY) + iix + DifX];
                                }
                            }
                        }
                    }
                    
                    cx[ii] = cx[ii] + DifX;
                    cy[ii] = cy[ii] + DifY;
                }
            }
        }
        
        //時間計測
        struct timeval now;
        gettimeofday(&now, NULL);
        double sec = (double)(now.tv_sec - start_time.tv_sec);
        double micro = (double)(now.tv_usec - start_time.tv_usec);
        double passed = sec + micro / 1000.0 / 1000.0;
        printf( "passed:%.12f\n", passed );
        printf( "toc\n" );
        
        start_time=now;
        
        for (int ii = 0; ii < N_cell; ii++){
            sprintf(filename, "data/result_cell%d_t%d_phi.txt", ii, tt + 1);
            fp1 = fopen(filename, "w");
            fprintf(fp1, "%d,%d\n", cx[ii], cy[ii]);
            for (int yy = 0; yy < L_comp; yy++){
                for (int xx = 0; xx < L_comp; xx++){
                    fprintf(fp1, "%f,", phi[ii][L_comp*yy + xx]);
                }
                fprintf(fp1, "\n");
            }
            fclose(fp1);
            
            sprintf(filename, "data/result_cell%d_t%d_phiu.txt", ii, tt + 1);
            fp1 = fopen(filename, "w");
            fprintf(fp1, "%d,%d\n", cx[ii], cy[ii]);
            for (int yy = 0; yy < L_comp; yy++){
                for (int xx = 0; xx < L_comp; xx++){
                    fprintf(fp1, "%f,", phiu[ii][L_comp*yy + xx]);
                }
                fprintf(fp1, "\n");
            }
            fclose(fp1);
            
            sprintf(filename, "data/result_cell%d_t%d_phiv.txt", ii, tt + 1);
            fp1 = fopen(filename, "w");
            fprintf(fp1, "%d,%d\n", cx[ii], cy[ii]);
            for (int yy = 0; yy < L_comp; yy++){
                for (int xx = 0; xx < L_comp; xx++){
                    fprintf(fp1, "%f,", phiv[ii][L_comp*yy + xx]);
                }
                fprintf(fp1, "\n");
            }
            fclose(fp1);
            
            sprintf(filename, "data/result_cell%d_t%d_phiz.txt", ii, tt + 1);
            fp1 = fopen(filename, "w");
            fprintf(fp1, "%d,%d\n", cx[ii], cy[ii]);
            for (int yy = 0; yy < L_comp; yy++){
                for (int xx = 0; xx < L_comp; xx++){
                    fprintf(fp1, "%f,", phiz[ii][L_comp*yy + xx]);
                }
                fprintf(fp1, "\n");
            }
            fclose(fp1);
            
        }
    }
    
}

