//
//  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"
#define PI 4.0*atan(1.0)

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;
    const int DDtt = 5000;
    const int Tmax = DDtt * 2500;
    
    double* phi = new double[n_comp];
    double* phiv = new double[n_comp];
    double* phiu = new double[n_comp];
    double* phiz = 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, cy;
    double radius = 5;
    cx = 800; cy = 800;
    
    int Poli_x, Poli_y;
    Poli_x = 1; Poli_y = 0;
    
    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[L_comp*iiy + iix] = 1.0;
                phiu[L_comp*iiy + iix] = 0.1;
                phiv[L_comp*iiy + iix] = 0.1;
                if ((iix - cx_grid)*Poli_x + (iiy - cy_grid)*Poli_y > 0) {
                    phiz[L_comp*iiy + iix] = 1.0;
                }else{
                    phiz[L_comp*iiy + iix] = 0.01;
                }
            }
            else{
                phi[L_comp*iiy + iix] = 0.0;
                phiu[L_comp*iiy + iix] = 0.0;
                phiv[L_comp*iiy + iix] = 0.0;
                phiz[L_comp*iiy + iix] = 0.0;
            }
        }
    }
    
    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)));
        }
    }
    
    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.5;
    theta = 1.4; d = 0.8; sigma = 0.075;
    double k_vz;
    k_vz = 200.0;
    double Dz, gw3, kw2, gw1, Ntot;
    Dz = 3.0; gw3 = 500.0; kw2 = 500.0; gw1 = 110.0; Ntot = 90.0;
    double tau, eta, epsilon, M, aw, A_0;
    tau = 0.83; eta = 1.0; epsilon = 1.0; M = 0.5; aw = 4.0; A_0 = 78.83;
    double lambda; lambda = 10.0;
    
    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_comp];
    double* Lap_v = new double[n_comp];
    double* Lap_z = new double[n_comp];
    double* Lap_phi = new double[n_comp];
    double* Nabla_phi_square = new double[n_comp];
    double* Nabla_phi_abs = new double[n_comp];
    double* diff_phiu = new double[n_comp];
    double* diff_phiv = new double[n_comp];
    double* diff_phiz = new double[n_comp];
    double* dif_phi0 = new double[n_comp];
    
    double* NoiseVadd = new double[n_comp];
    for (int iiy = 0; iiy < L_comp; iiy++){
        for (int iix = 0; iix < L_comp; iix++){
            NoiseVadd[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, sumz, sumv2, sumu, sumNablaphi2;
    
    double lambda_phi = 1e-4;
    
    double avephixp1, avephixm1, avephiyp1, avephiym1;
    
    //
    double p_one = theta*dt*dx*dx;
    int n_max = 10;
    int n_tmp;
    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_comp],phiv[0:n_comp],phiu[0:n_comp],phiz[0:n_comp]) create(Lap_phi[0:n_comp], Lap_u[0:n_comp], Lap_v[0:n_comp],Lap_z[0:n_comp], Nabla_phi_square[0:n_comp], Nabla_phi_abs[0:n_comp],diff_phiu[0:n_comp],diff_phiv[0:n_comp],diff_phiz[0:n_comp],dif_phi0[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++){
            double p_tmp = genrand_real3();
            for (int i = 0; i <= n_max; i++){
                n_tmp = i;
                if (p_tmp < prob_s[i]) {
                    break;
                }
            }
            if (n_tmp > 0) {
                for (int iiy = 0; iiy < L_comp; iiy++){
                    for (int iix = 0; iix < L_comp; iix++){
                        NoiseVadd[L_comp*iiy + iix] = 0.0;
                    }
                }
                
                for (int i = 0; i < n_tmp; i++){
                    iiy_noise_tmp = genrand_int32()%L_comp;
                    iix_noise_tmp = genrand_int32()%L_comp;
                    double ftmp = genrand_real3();
                    st_noise_tmp = -log(ftmp);
                    
                    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[L_comp*iiy + iix] += st_noise_tmp; }
                                if (xx*xx+yy*yy==1) { NoiseVadd[L_comp*iiy + iix] += fac11*st_noise_tmp; }
                                if (xx*xx+yy*yy==2) { NoiseVadd[L_comp*iiy + iix] += fac12*st_noise_tmp; }
                                if (xx*xx+yy*yy==4) { NoiseVadd[L_comp*iiy + iix] += fac21*st_noise_tmp; }
                                if (xx*xx+yy*yy==5) { NoiseVadd[L_comp*iiy + iix] += fac22*st_noise_tmp; }
                                if (xx*xx+yy*yy==8) { NoiseVadd[L_comp*iiy + iix] += fac23*st_noise_tmp; }
                            }
                        }
                    }
                    
                }
                
            }
            
            int Rx = (int)(cx - (L_comp / 2 - 1));
            int Ry = (int)(cy - (L_comp / 2 - 1));
            
            sum = 0.0;
#pragma acc parallel loop reduction(+:sum)
            for (int rr = 0; rr < n_comp; rr++){
                sum += phi[rr];
            }
            sumphi = sum;
            
            sum = 0.0;
#pragma acc parallel loop reduction(+:sum)
            for (int rr = 0; rr < n_comp; rr++){
                sum += phi[rr]*phiz[rr];
            }
            sumz = sum;
            
            sum = 0.0;
#pragma acc parallel loop reduction(+:sum)
            for (int rr = 0; rr < n_comp; rr++){
                sum += phi[rr]*phiu[rr];
            }
            sumu = sum;
            
            sum = 0.0;
#pragma acc parallel loop reduction(+:sum)
            for (int rr = 0; rr < n_comp; rr++){
                sum += phi[rr]*phiz[rr]*phiv[rr];
            }
            sumv2 = sum;
            
            
#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[L_comp*iiy + iix + 1] + phi[L_comp*iiy + iix]);
                    avephixm1 = 0.5*(phi[L_comp*iiy + iix - 1] + phi[L_comp*iiy + iix]);
                    avephiyp1 = 0.5*(phi[L_comp*(iiy + 1) + iix] + phi[L_comp*iiy + iix]);
                    avephiym1 = 0.5*(phi[L_comp*(iiy - 1) + iix] + phi[L_comp*iiy + iix]);
                    Lap_u[L_comp*iiy + iix] = (avephixp1*(phiu[L_comp*iiy + iix + 1] - phiu[L_comp*iiy + iix]) - avephixm1*(phiu[L_comp*iiy + iix] - phiu[L_comp*iiy + iix - 1])) + (avephiyp1*(phiu[L_comp*(iiy + 1) + iix] - phiu[L_comp*iiy + iix]) - avephiym1*(phiu[L_comp*iiy + iix] - phiu[L_comp*(iiy - 1) + iix]));
                    Lap_v[L_comp*iiy + iix] = (avephixp1*(phiv[L_comp*iiy + iix + 1] - phiv[L_comp*iiy + iix]) - avephixm1*(phiv[L_comp*iiy + iix] - phiv[L_comp*iiy + iix - 1])) + (avephiyp1*(phiv[L_comp*(iiy + 1) + iix] - phiv[L_comp*iiy + iix]) - avephiym1*(phiv[L_comp*iiy + iix] - phiv[L_comp*(iiy - 1) + iix]));
                    Lap_z[L_comp*iiy + iix] = (avephixp1*(phiz[L_comp*iiy + iix + 1] - phiz[L_comp*iiy + iix]) - avephixm1*(phiz[L_comp*iiy + iix] - phiz[L_comp*iiy + iix - 1])) + (avephiyp1*(phiz[L_comp*(iiy + 1) + iix] - phiz[L_comp*iiy + iix]) - avephiym1*(phiz[L_comp*iiy + iix] - phiz[L_comp*(iiy - 1) + iix]));
                    Lap_phi[L_comp*iiy + iix] = -4.0 * phi[L_comp*iiy + iix] + phi[L_comp*(iiy + 1) + iix] + phi[L_comp*(iiy - 1) + iix] + phi[L_comp*iiy + iix + 1] + phi[L_comp*iiy + iix - 1];
                    Nabla_phi_square[L_comp*iiy + iix] = 0.5*(phi[L_comp*(iiy + 1) + iix] - phi[L_comp*(iiy - 1) + iix])*0.5*(phi[L_comp*(iiy + 1) + iix] - phi[L_comp*(iiy - 1) + iix]) + 0.5*(phi[L_comp*iiy + iix + 1] - phi[L_comp*iiy + iix - 1])*0.5*(phi[L_comp*iiy + iix + 1] - phi[L_comp*iiy + iix - 1]);
                }
                if (iiy == 0 && iix>0 && iix < L_comp - 1) {
                    avephixp1 = 0.5*(phi[L_comp*iiy + iix + 1] + phi[L_comp*iiy + iix]);
                    avephixm1 = 0.5*(phi[L_comp*iiy + iix - 1] + phi[L_comp*iiy + iix]);
                    avephiyp1 = 0.5*(phi[L_comp*(iiy + 1) + iix] + phi[L_comp*iiy + iix]);
                    Lap_u[L_comp*iiy + iix] = (avephixp1*(phiu[L_comp*iiy + iix + 1] - phiu[L_comp*iiy + iix]) - avephixm1*(phiu[L_comp*iiy + iix] - phiu[L_comp*iiy + iix - 1])) + (avephiyp1*(phiu[L_comp*(iiy + 1) + iix] - phiu[L_comp*iiy + iix]));
                    Lap_v[L_comp*iiy + iix] = (avephixp1*(phiv[L_comp*iiy + iix + 1] - phiv[L_comp*iiy + iix]) - avephixm1*(phiv[L_comp*iiy + iix] - phiv[L_comp*iiy + iix - 1])) + (avephiyp1*(phiv[L_comp*(iiy + 1) + iix] - phiv[L_comp*iiy + iix]));
                    Lap_z[L_comp*iiy + iix] = (avephixp1*(phiz[L_comp*iiy + iix + 1] - phiz[L_comp*iiy + iix]) - avephixm1*(phiz[L_comp*iiy + iix] - phiz[L_comp*iiy + iix - 1])) + (avephiyp1*(phiz[L_comp*(iiy + 1) + iix] - phiz[L_comp*iiy + iix]));
                    Lap_phi[L_comp*iiy + iix] = -4.0 * phi[L_comp*iiy + iix] + phi[L_comp*(iiy + 1) + iix] + phi[L_comp*iiy + iix] + phi[L_comp*iiy + iix + 1] + phi[L_comp*iiy + iix - 1];
                    Nabla_phi_square[L_comp*iiy + iix] = 0.5*(phi[L_comp*(iiy + 1) + iix] - phi[L_comp*iiy + iix])*0.5*(phi[L_comp*(iiy + 1) + iix] - phi[L_comp*iiy + iix]) + 0.5*(phi[L_comp*iiy + iix + 1] - phi[L_comp*iiy + iix - 1])*0.5*(phi[L_comp*iiy + iix + 1] - phi[L_comp*iiy + iix - 1]);
                }
                if (iiy == L_comp - 1 && iix>0 && iix < L_comp - 1) {
                    avephixp1 = 0.5*(phi[L_comp*iiy + iix + 1] + phi[L_comp*iiy + iix]);
                    avephixm1 = 0.5*(phi[L_comp*iiy + iix - 1] + phi[L_comp*iiy + iix]);
                    avephiym1 = 0.5*(phi[L_comp*(iiy - 1) + iix] + phi[L_comp*iiy + iix]);
                    Lap_u[L_comp*iiy + iix] = (avephixp1*(phiu[L_comp*iiy + iix + 1] - phiu[L_comp*iiy + iix]) - avephixm1*(phiu[L_comp*iiy + iix] - phiu[L_comp*iiy + iix - 1])) - avephiym1*(phiu[L_comp*iiy + iix] - phiu[L_comp*(iiy - 1) + iix]);
                    Lap_v[L_comp*iiy + iix] = (avephixp1*(phiv[L_comp*iiy + iix + 1] - phiv[L_comp*iiy + iix]) - avephixm1*(phiv[L_comp*iiy + iix] - phiv[L_comp*iiy + iix - 1])) - avephiym1*(phiv[L_comp*iiy + iix] - phiv[L_comp*(iiy - 1) + iix]);
                    Lap_z[L_comp*iiy + iix] = (avephixp1*(phiz[L_comp*iiy + iix + 1] - phiz[L_comp*iiy + iix]) - avephixm1*(phiz[L_comp*iiy + iix] - phiz[L_comp*iiy + iix - 1])) - avephiym1*(phiz[L_comp*iiy + iix] - phiz[L_comp*(iiy - 1) + iix]);
                    Lap_phi[L_comp*iiy + iix] = -4.0 * phi[L_comp*iiy + iix] + phi[L_comp*iiy + iix] + phi[L_comp*(iiy - 1) + iix] + phi[L_comp*iiy + iix + 1] + phi[L_comp*iiy + iix - 1];
                    Nabla_phi_square[L_comp*iiy + iix] = 0.5*(phi[L_comp*iiy + iix] - phi[L_comp*(iiy - 1) + iix])*0.5*(phi[L_comp*iiy + iix] - phi[L_comp*(iiy - 1) + iix]) + 0.5*(phi[L_comp*iiy + iix + 1] - phi[L_comp*iiy + iix - 1])*0.5*(phi[L_comp*iiy + iix + 1] - phi[L_comp*iiy + iix - 1]);
                }
                if (iiy>0 && iiy < L_comp - 1 && iix == 0) {
                    avephixp1 = 0.5*(phi[L_comp*iiy + iix + 1] + phi[L_comp*iiy + iix]);
                    avephiyp1 = 0.5*(phi[L_comp*(iiy + 1) + iix] + phi[L_comp*iiy + iix]);
                    avephiym1 = 0.5*(phi[L_comp*(iiy - 1) + iix] + phi[L_comp*iiy + iix]);
                    Lap_u[L_comp*iiy + iix] = (avephixp1*(phiu[L_comp*iiy + iix + 1] - phiu[L_comp*iiy + iix])) + (avephiyp1*(phiu[L_comp*(iiy + 1) + iix] - phiu[L_comp*iiy + iix]) - avephiym1*(phiu[L_comp*iiy + iix] - phiu[L_comp*(iiy - 1) + iix]));
                    Lap_v[L_comp*iiy + iix] = (avephixp1*(phiv[L_comp*iiy + iix + 1] - phiv[L_comp*iiy + iix])) + (avephiyp1*(phiv[L_comp*(iiy + 1) + iix] - phiv[L_comp*iiy + iix]) - avephiym1*(phiv[L_comp*iiy + iix] - phiv[L_comp*(iiy - 1) + iix]));
                    Lap_z[L_comp*iiy + iix] = (avephixp1*(phiz[L_comp*iiy + iix + 1] - phiz[L_comp*iiy + iix])) + (avephiyp1*(phiz[L_comp*(iiy + 1) + iix] - phiz[L_comp*iiy + iix]) - avephiym1*(phiz[L_comp*iiy + iix] - phiz[L_comp*(iiy - 1) + iix]));
                    Lap_phi[L_comp*iiy + iix] = -4.0 * phi[L_comp*iiy + iix] + phi[L_comp*(iiy + 1) + iix] + phi[L_comp*(iiy - 1) + iix] + phi[L_comp*iiy + iix + 1] + phi[L_comp*iiy + iix];
                    Nabla_phi_square[L_comp*iiy + iix] = 0.5*(phi[L_comp*(iiy + 1) + iix] - phi[L_comp*(iiy - 1) + iix])*0.5*(phi[L_comp*(iiy + 1) + iix] - phi[L_comp*(iiy - 1) + iix]) + 0.5*(phi[L_comp*iiy + iix + 1] - phi[L_comp*iiy + iix])*0.5*(phi[L_comp*iiy + iix + 1] - phi[L_comp*iiy + iix]);
                }
                if (iiy>0 && iiy < L_comp - 1 && iix == L_comp - 1) {
                    avephixm1 = 0.5*(phi[L_comp*iiy + iix - 1] + phi[L_comp*iiy + iix]);
                    avephiyp1 = 0.5*(phi[L_comp*(iiy + 1) + iix] + phi[L_comp*iiy + iix]);
                    avephiym1 = 0.5*(phi[L_comp*(iiy - 1) + iix] + phi[L_comp*iiy + iix]);
                    Lap_u[L_comp*iiy + iix] = (-avephixm1*(phiu[L_comp*iiy + iix] - phiu[L_comp*iiy + iix - 1])) + (avephiyp1*(phiu[L_comp*(iiy + 1) + iix] - phiu[L_comp*iiy + iix]) - avephiym1*(phiu[L_comp*iiy + iix] - phiu[L_comp*(iiy - 1) + iix]));
                    Lap_v[L_comp*iiy + iix] = (-avephixm1*(phiv[L_comp*iiy + iix] - phiv[L_comp*iiy + iix - 1])) + (avephiyp1*(phiv[L_comp*(iiy + 1) + iix] - phiv[L_comp*iiy + iix]) - avephiym1*(phiv[L_comp*iiy + iix] - phiv[L_comp*(iiy - 1) + iix]));
                    Lap_z[L_comp*iiy + iix] = (-avephixm1*(phiz[L_comp*iiy + iix] - phiz[L_comp*iiy + iix - 1])) + (avephiyp1*(phiz[L_comp*(iiy + 1) + iix] - phiz[L_comp*iiy + iix]) - avephiym1*(phiz[L_comp*iiy + iix] - phiz[L_comp*(iiy - 1) + iix]));
                    Lap_phi[L_comp*iiy + iix] = -4.0 * phi[L_comp*iiy + iix] + phi[L_comp*(iiy + 1) + iix] + phi[L_comp*(iiy - 1) + iix] + phi[L_comp*iiy + iix] + phi[L_comp*iiy + iix - 1];
                    Nabla_phi_square[L_comp*iiy + iix] = 0.5*(phi[L_comp*(iiy + 1) + iix] - phi[L_comp*(iiy - 1) + iix])*0.5*(phi[L_comp*(iiy + 1) + iix] - phi[L_comp*(iiy - 1) + iix]) + 0.5*(phi[L_comp*iiy + iix] - phi[L_comp*iiy + iix - 1])*0.5*(phi[L_comp*iiy + iix] - phi[L_comp*iiy + iix - 1]);
                }
                if ((iiy == 0 || iiy == L_comp - 1) && (iix == 0 || iix == L_comp - 1)) {
                    Lap_u[L_comp*iiy + iix] = 0.0;
                    Lap_v[L_comp*iiy + iix] = 0.0;
                    Lap_z[L_comp*iiy + iix] = 0.0;
                    Lap_phi[L_comp*iiy + iix] = 0.0;
                    Nabla_phi_square[L_comp*iiy + iix] = 0.0;
                }
                
                if (Nabla_phi_square[L_comp*iiy + iix] < 0.0){
                    Nabla_phi_square[L_comp*iiy + iix] = 0.0;
                }
                Nabla_phi_abs[L_comp*iiy + iix] = sqrt(Nabla_phi_square[L_comp*iiy + iix]);
            }
            
            sum = 0.0;
#pragma acc parallel loop reduction(+:sum)
            for (int rr = 0; rr < n_comp; rr++){
                sum += Nabla_phi_square[rr];
            }
            sumNablaphi2 = sum;
            
#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[rr] = (Du*Lap_u[rr] / (dx*dx) + phi[rr] * (-alpha*phiz[rr]*phiu[rr] * phiv[rr] / (Kk + sumv2 / sumphi) + beta*phiu[rr] * phiv[rr] / (Kp + sumu / sumphi) + s - gamma*phiu[rr]) - kaiU*phiu[rr]*Nabla_phi_square[rr]/(sumNablaphi2*dx*dx));
                    diff_phiv[rr] = (Dv*Lap_v[rr] / (dx*dx) + phi[rr] * (alpha*phiz[rr]*phiu[rr] * phiv[rr] / (Kk + sumv2 / sumphi) - beta*phiu[rr] * phiv[rr] / (Kp + sumu / sumphi) - mu*phiv[rr]) - kaiV*phiv[rr]*Nabla_phi_square[rr]/(sumNablaphi2*dx*dx));
                    diff_phiz[rr] = (Dz*Lap_z[rr] / (dx*dx) + phi[rr] * (- gw3*phiz[rr]*phiz[rr]*phiz[rr] + kw2*(Ntot - sumz*dx*dx)/(sumphi*dx*dx)*phiz[rr]*phiz[rr] - gw1*phiz[rr] + k_vz*phiv[rr]));
                    dif_phi0[rr] = (1.0 / tau)*(eta*(Lap_phi[rr] / (dx*dx) - (1.0 / (epsilon*epsilon)) * 72.0 * phi[rr] * (1.0 - phi[rr])*(0.5 - phi[rr])) + M*(A_0 - sumphi * dx*dx)*Nabla_phi_abs[rr] / dx + (aw*phiz[rr])*Nabla_phi_abs[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[rr] += -lambda*phi[rr]*(1.0 - kai_grid[Ry + iiy][Rx + iix])*(1.0 - kai_grid[Ry + iiy][Rx + iix]);
                }
            }
            
#pragma acc data copyin(NoiseVadd[0:n_comp])
            {
                
#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)) {
                        if (phi[rr] > lambda_phi) {
                            if (n_tmp > 0) {
                                if (sigma*NoiseVadd[rr]>phiu[rr]){
                                    phiu[rr] = (phiu[rr]*phi[rr] + dt*diff_phiu[rr])/phi[rr] - phiu[rr];
                                    phiv[rr] = (phiv[rr]*phi[rr] + dt*diff_phiv[rr])/phi[rr] + phiu[rr];
                                }else{
                                    phiu[rr] = (phiu[rr]*phi[rr] + dt*diff_phiu[rr])/phi[rr] - sigma*NoiseVadd[rr];
                                    phiv[rr] = (phiv[rr]*phi[rr] + dt*diff_phiv[rr])/phi[rr] + sigma*NoiseVadd[rr];
                                }
                            }else{
                                phiu[rr] = (phiu[rr]*phi[rr] + dt*diff_phiu[rr])/phi[rr];
                                phiv[rr] = (phiv[rr]*phi[rr] + dt*diff_phiv[rr])/phi[rr];
                            }
                            phiz[rr] = (phiz[rr]*phi[rr] + dt*diff_phiz[rr])/phi[rr];
                        }else{
                            phiu[rr] = (phiu[rr]*phi[rr] + dt*diff_phiu[rr]);
                            phiv[rr] = (phiv[rr]*phi[rr] + dt*diff_phiv[rr]);
                            phiz[rr] = (phiz[rr]*phi[rr] + dt*diff_phiz[rr]);
                        }
                        phi[rr] += (dt*dif_phi0[rr]);
                    }
                }
            }
            
            
            if (t2 == DDtt - 1) {
                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[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[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[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);
                
#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[L_comp*iiy + iix];
                        phiz_t[L_comp*iiy + iix] = phiz[L_comp*iiy + iix];
                        phiv_t[L_comp*iiy + iix] = phiv[L_comp*iiy + iix];
                        phiu_t[L_comp*iiy + iix] = phiu[L_comp*iiy + iix];
                        phi[L_comp*iiy + iix] = 0.0;
                        phiz[L_comp*iiy + iix] = 0.0;
                        phiv[L_comp*iiy + iix] = 0.0;
                        phiu[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[L_comp*iiy + iix] = phi_t[L_comp*(iiy + DifY) + iix + DifX];
                                phiv[L_comp*iiy + iix] = phiv_t[L_comp*(iiy + DifY) + iix + DifX];
                                phiu[L_comp*iiy + iix] = phiu_t[L_comp*(iiy + DifY) + iix + DifX];
                                phiz[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[L_comp*iiy + iix] = phi_t[L_comp*(iiy + DifY) + iix + DifX];
                                phiv[L_comp*iiy + iix] = phiv_t[L_comp*(iiy + DifY) + iix + DifX];
                                phiu[L_comp*iiy + iix] = phiu_t[L_comp*(iiy + DifY) + iix + DifX];
                                phiz[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[L_comp*iiy + iix] = phi_t[L_comp*(iiy + DifY) + iix + DifX];
                                phiv[L_comp*iiy + iix] = phiv_t[L_comp*(iiy + DifY) + iix + DifX];
                                phiu[L_comp*iiy + iix] = phiu_t[L_comp*(iiy + DifY) + iix + DifX];
                                phiz[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[L_comp*iiy + iix] = phi_t[L_comp*(iiy + DifY) + iix + DifX];
                                phiv[L_comp*iiy + iix] = phiv_t[L_comp*(iiy + DifY) + iix + DifX];
                                phiu[L_comp*iiy + iix] = phiu_t[L_comp*(iiy + DifY) + iix + DifX];
                                phiz[L_comp*iiy + iix] =phiz_t[L_comp*(iiy + DifY) + iix + DifX];
                            }
                        }
                    }
                }
                
                cx = cx + DifX;
                cy = cy + 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;
        
        sprintf(filename, "data_tmp/result_tmp.txt");
        fp1 = fopen(filename, "w");
        fprintf(fp1, "%d\n", tt + 1);
        fclose(fp1);
        
        sprintf(filename, "data_tmp/result_cell%d_phi.txt", 0);
        fp1 = fopen(filename, "w");
        fprintf(fp1, "%d,%d\n", cx, cy);
        for (int yy = 0; yy < L_comp; yy++){
            for (int xx = 0; xx < L_comp; xx++){
                fprintf(fp1, "%f,", phi[L_comp*yy + xx]);
            }
            fprintf(fp1, "\n");
        }
        fclose(fp1);
        
        sprintf(filename, "data_tmp/result_cell%d_phiu.txt", 0);
        fp1 = fopen(filename, "w");
        fprintf(fp1, "%d,%d\n", cx, cy);
        for (int yy = 0; yy < L_comp; yy++){
            for (int xx = 0; xx < L_comp; xx++){
                fprintf(fp1, "%f,", phiu[L_comp*yy + xx]);
            }
            fprintf(fp1, "\n");
        }
        fclose(fp1);
        
        sprintf(filename, "data_tmp/result_cell%d_phiv.txt", 0);
        fp1 = fopen(filename, "w");
        fprintf(fp1, "%d,%d\n", cx, cy);
        for (int yy = 0; yy < L_comp; yy++){
            for (int xx = 0; xx < L_comp; xx++){
                fprintf(fp1, "%f,", phiv[L_comp*yy + xx]);
            }
            fprintf(fp1, "\n");
        }
        fclose(fp1);
        
        sprintf(filename, "data_tmp/result_cell%d_phiz.txt", 0);
        fp1 = fopen(filename, "w");
        fprintf(fp1, "%d,%d\n", cx, cy);
        for (int yy = 0; yy < L_comp; yy++){
            for (int xx = 0; xx < L_comp; xx++){
                fprintf(fp1, "%f,", phiz[L_comp*yy + xx]);
            }
            fprintf(fp1, "\n");
        }
        fclose(fp1);
        
        sprintf(filename, "data/result_t%d.txt", tt + 1);
        fp1 = fopen(filename, "w");
        fprintf(fp1, "%d,%d,%d,0,0,0,0\n", 0, cx, cy);
        for (int yy = 0; yy < L_comp; yy++){
            for (int xx = 0; xx < L_comp; xx++){
                if (phi[L_comp*yy + xx]>lambda_phi) {
                    fprintf(fp1, "%d,%d,%d,%f,%f,%f,%f\n", 0, yy, xx, phi[L_comp*yy + xx], phiv[L_comp*yy + xx], phiu[L_comp*yy + xx], phiz[L_comp*yy + xx]);
                }
            }
        }
        fclose(fp1);
        
    }
    
}

