make parameter of number of scalar equations

Makefile

@@ -4,76 +4,13 @@
 
 MPIF90 = blah
 
-ifeq ($(OMPI), 1)
-        FFTW_HOME = ../libs/fftw-3.2.2/lib
-        OPENMPI_HOME = /home/comb/ignis/openmpi/
-        MPIF90 = $(OPENMPI_HOME)/bin/mpif90
-        FCFLAGS = -i8 -r8 -c -openmp $(MPI_COMPILE_FLAGS) -I$(FFTW_HOME)/include -I$(OPENMPI_HOME)/include
-        LDFLAGS = -i8 -r8 -openmp $(MPI_LD_FLAGS) -L$(FFTW_HOME)/lib -lfftw3 -L$(OPENMPI_HOME)/lib -lfftw3 -lm  -shared-intel
-        FCFLAGS_F77 = -Mextend
+        FFTW_HOME = ../libs/fftw
+        MPIF90 = mpif90
+        FCFLAGS = -fdefault-double-8  -fdefault-integer-8  -fdefault-real-8 -c $(MPI_COMPILE_FLAGS) -I$(FFTW_HOME)/include
+        LDFLAGS = -fdefault-double-8  -fdefault-integer-8  -fdefault-real-8 $(MPI_LD_FLAGS) -L$(FFTW_HOME)/lib -lfftw3 -lm
+        FCFLAGS_F77 =
 
         MODULES += m_openmpi.o
-else
-        FFTW_HOME = ../libs/fftw-3.2.2/lib
-        MPIF90 = mpif90
-        FCFLAGS = -i8 -r8 -c $(MPI_COMPILE_FLAGS) -I$(FFTW_HOME)/include
-        LDFLAGS = -i8 -r8 $(MPI_LD_FLAGS) -L$(FFTW_HOME)/lib -lfftw3 -lm
-        FCFLAGS_F77 = -Mextend
-
-        MODULES += m_mpi.o
-endif
-
-# Franklin cluster, NERSC
-ifeq ($(NERSC_HOST), franklin)
-        MPIF90 = ftn
-
-        FCFLAGS = -target=linux -i8 -r8 -O4 -c
-        LDFLAGS = -target=linux -i8 -r8 -O4 -lfftw3
-        FCFLAGS_F77 = -Mextend
-endif
-
-# Hopper cluster, NERSC
-ifeq ($(NERSC_HOST), hopper)
-        MPIF90 = ftn
-
-        FCFLAGS = -target=linux -i8 -r8 -O4 -c
-        LDFLAGS = -target=linux -i8 -r8 -O4 -lfftw3
-        FCFLAGS_F77 = -Mextend
-endif
-
-# Yellowrail cluster, LANL
-ifeq ($(HOSTNAME), yr-fe1.lanl.gov)
-	MPIF90 = mpif90
-
-	FCFLAGS = -i8 -r8 -O4 -c $(MPI_COMPILE_FLAGS) -I$(FFTW_INCLUDE)
-	LDFLAGS = -i8 -r8 -O4 -lfftw3 $(MPI_LD_FLAGS) -L$(FFTW_HOME)/lib
-endif
-
-# Coyote cluster, LANL
-ifeq ($(HOSTNAME), cy-c2.lanl.gov)
-	MPIF90 = mpif90
-	FCFLAGS = -i8 -r8 -O4 -c $(MPI_COMPILE_FLAGS) -I$(FFTW_INCLUDE)
-	LDFLAGS = -i8 -r8 -O4 -lfftw3 $(MPI_LD_FLAGS) -L$(FFTW_HOME)/lib
-endif
-
-# Linux Cluster, CMU 
-ifeq ($(HOSTNAME), karman.me.cmu.edu)
-        FFTW_HOME = ../libs/fftw-3.2.2/lib
-        MPIF90 = mpif90
-        FCFLAGS = -i8 -r8 -c $(MPI_COMPILE_FLAGS) -I$(FFTW_HOME)/include
-        LDFLAGS = -i8 -r8 $(MPI_LD_FLAGS) -L$(FFTW_HOME)/lib -lfftw3 -lm
-        FCFLAGS_F77 = -Mextend
-
-endif
-
-# sparrow.stanford.edu, Mac OS X box
-ifeq ($(HOSTNAME), sparrow.stanford.edu)
-        MPIF90 = mpif90
-        FCFLAGS = -fdefault-real-8 -fdefault-integer-8 -finit-integer=0 -finit-real=zero -c 
-        FCFLAGS_F77 = -ffixed-form -ffixed-line-length-none
-        FCFLAGS_F90 = -ffree-form -ffree-line-length-none
-        LDFLAGS = -fdefault-real-8 -fdefault-integer-8 -finit-integer=0 -finit-real=zero -lmpi -lmpi -lfftw3 -lm
-endif
 
 # Program name
 PROG    = hit3d.x

gather_4.f90

@@ -45,18 +45,15 @@ subroutine gather_4
   ! velocities
   call xFFT3d(-1,1)
   fname = 'u.'//file_ext
-  tmp4(1:nx,1:ny,1:nz) = wrk(1:nx,1:ny,1:nz,1)
-  call gather_tmp4(u_)
+  u_(1:nx,1:ny,1:nz) = wrk(1:nx,1:ny,1:nz,1)
 
   call xFFT3d(-1,2)
   fname = 'v.'//file_ext
-  tmp4(1:nx,1:ny,1:nz) = wrk(1:nx,1:ny,1:nz,2)
-  call gather_tmp4(v_)
+  v_(1:nx,1:ny,1:nz) = wrk(1:nx,1:ny,1:nz,2)
 
   call xFFT3d(-1,3)
   fname = 'w.'//file_ext
-  tmp4(1:nx,1:ny,1:nz) = wrk(1:nx,1:ny,1:nz,3)
-  call gather_tmp4(w_)
+  w_(1:nx,1:ny,1:nz) = wrk(1:nx,1:ny,1:nz,3)
 
   ! scalars
   if (int_scalars) then

init_velocity.f90

@@ -266,8 +266,10 @@ subroutine init_velocity
         do i = 1,nx+2
 
            n_shell = nint(sqrt(real(akx(i)**2 + aky(k)**2 + akz(j)**2, 4)))
-           if (n_shell .gt. 0 .and. n_shell .le. kmax .and. e_spec(n_shell) .gt. zip) then
+           if (n_shell .gt. 0 .and. n_shell .le. kmax ) then
+              if (e_spec(n_shell) .gt. zip) then
               fields(i,j,k,1:3) = fields(i,j,k,1:3) * sqrt(e_spec1(n_shell)/e_spec(n_shell))
+              endif
            else
               fields(i,j,k,1:3) = zip
            end if

io_write_4.f90

@@ -65,7 +65,6 @@ subroutine io_write_4
   wmax = maxval(w_(1:nx,1:ny,1:nz_all))
   wmin = minval(w_(1:nx,1:ny,1:nz_all))
 
-  call write_tp_field
   call write_vel_max ( umax,vmax,wmax,umin,vmin,wmin )
 
   endif

m_compact.f90

@@ -1,4 +1,7 @@
     MODULE m_compact
+
+  use m_openmpi
+
      IMPLICIT NONE
      PRIVATE
      REAL*8, DIMENSION(:), ALLOCATABLE :: lxf,lxs,wxf,wxs, & 
@@ -8,7 +11,7 @@
      INTEGER :: nxc,nyc,nzc
      REAL*8, PARAMETER :: ezero = 1.0e-14
       
-     PUBLIC :: ludcmp,dfnonp,d2fnonp,dfp,d2fp
+     PUBLIC :: ludcmp,dfnonp,d2fnonp,dfp,d2fp,par_dfp,par_d2fp
 
     CONTAINS
 
@@ -219,6 +222,85 @@
 
      END SUBROUTINE dfp
 
+     SUBROUTINE par_dfp(x,xu,xl,dx,h,nd,n,nall,dir)
+      INTEGER,INTENT(IN) :: nall,n,nd,dir
+      REAL*8,INTENT(IN) :: h
+      REAL*8,INTENT(IN),DIMENSION(nd,n) :: x
+      REAL*8,INTENT(OUT),DIMENSION(nd,n) :: dx
+      INTEGER :: i,j
+      REAL*8 :: r1,r2,h1
+
+      REAL*8,DIMENSION(nd,2) :: xu, xl
+
+      h1=1./h
+      r1=7./3.
+      r2=1./12.
+
+      !if (myid.eq.0) write(*,*) "parallel dfp commuication"
+      if (myid.eq.master) then
+
+              call MPI_ISEND (x(1,1),   2*nd, MPI_REAL8, numprocs-1, myid, &
+                              MPI_COMM_TASK, mpi_request, mpi_err)
+
+              call MPI_RECV  (xu,       2*nd, MPI_REAL8, myid+1,     myid+1, &
+                              MPI_COMM_TASK, mpi_status, mpi_err)
+
+              call MPI_ISEND (x(1,n-1), 2*nd, MPI_REAL8, myid+1,     myid, &
+                              MPI_COMM_TASK, mpi_request, mpi_err)
+
+              call MPI_RECV  (xl,       2*nd, MPI_REAL8, numprocs-1, numprocs-1, &
+                              MPI_COMM_TASK, mpi_status, mpi_err)
+
+      elseif (myid.eq.numprocs-1) then
+
+              call MPI_ISEND (x(1,1),   2*nd, MPI_REAL8, myid-1, myid, &
+                              MPI_COMM_TASK, mpi_request, mpi_err)
+
+              call MPI_RECV  (xu,       2*nd, MPI_REAL8, 0,      0, &
+                              MPI_COMM_TASK, mpi_status, mpi_err)
+
+              call MPI_ISEND (x(1,n-1), 2*nd, MPI_REAL8, 0,      myid, &
+                              MPI_COMM_TASK, mpi_request, mpi_err)
+
+              call MPI_RECV  (xl,       2*nd, MPI_REAL8, myid-1, myid-1, &
+                              MPI_COMM_TASK, mpi_status, mpi_err)
+
+      else
+
+              call MPI_ISEND (x(1,1),   2*nd, MPI_REAL8, myid-1, myid, &
+                              MPI_COMM_TASK, mpi_request, mpi_err)
+
+              call MPI_RECV  (xu,       2*nd, MPI_REAL8, myid+1, myid+1, &
+                              MPI_COMM_TASK, mpi_status, mpi_err)
+
+              call MPI_ISEND (x(1,n-1), 2*nd, MPI_REAL8, myid+1, myid, &
+                              MPI_COMM_TASK, mpi_request, mpi_err)
+
+              call MPI_RECV  (xl,       2*nd, MPI_REAL8, myid-1, myid-1, &
+                              MPI_COMM_TASK, mpi_status, mpi_err)
+
+      endif
+
+      DO j=1,nd
+        dx(j,n-1)=(r1*( x(j,n)- x(j,n-2))+r2*(xu(j,1)- x(j,n-3)))*h1
+        dx(j,n)=  (r1*(xu(j,1)- x(j,n-1))+r2*(xu(j,2)- x(j,n-2)))*h1
+        dx(j,1)=  (r1*( x(j,2)-xl(j,2))  +r2*( x(j,3)-xl(j,1)))*h1
+        dx(j,2)=  (r1*( x(j,3)- x(j,1))  +r2*( x(j,4)-xl(j,2)))*h1
+      ENDDO
+
+      DO i=3,n-2
+       DO j=1,nd
+        dx(j,i)=(r1*(x(j,i+1)-x(j,i-1))+r2*(x(j,i+2)-x(j,i-2)))*h1
+       ENDDO
+      ENDDO
+
+
+      IF (dir.eq.1) CALL par_ptdslv(dx,lxf,wxf,nd,n,nall) ! x-direction
+      IF (dir.eq.2) CALL par_ptdslv(dx,lyf,wyf,nd,n,nall) ! x-direction
+      IF (dir.eq.3) CALL par_ptdslv(dx,lzf,wzf,nd,n,nall) ! x-direction
+
+     END SUBROUTINE par_dfp
+
      SUBROUTINE ptdslv(r,n,l,w,nd)
       INTEGER,INTENT(IN) :: n,nd
       REAL*8,INTENT(INOUT),DIMENSION(nd,n) :: r
@@ -247,6 +329,183 @@
       ENDDO
      END SUBROUTINE ptdslv
 
+     SUBROUTINE par_ptdslv(r,la,wa,nd,n,nall)
+      INTEGER,PARAMETER :: nb = 128
+      INTEGER,INTENT(IN) :: n,nd,nall
+      REAL*8,INTENT(INOUT),DIMENSION(nd,n) :: r
+      REAL*8,INTENT(IN),DIMENSION(nall) :: la,wa
+      REAL*8,DIMENSION(n) :: l,w
+      INTEGER i,j
+      INTEGER ii,jj
+      REAL*8, DIMENSION(nd) :: sum
+
+      INTEGER npart, nbase, nlow, nupp
+      INTEGER :: pid, np
+      REAL*8, DIMENSION(nb) :: r0, r1, sum0
+      REAL*8, DIMENSION(nb,2) :: buf
+
+      pid = myid
+      np = numprocs
+
+      npart = nall / np
+      nbase = pid * npart
+      nlow = pid * npart + 1
+      nupp = (pid + 1) * npart
+
+      l = la(nlow:nupp)
+      w = wa(nlow:nupp)
+
+      if (npart.lt.4) then
+        ! assertion fail
+      endif
+
+      call MPI_BARRIER(MPI_COMM_WORLD, mpi_err)
+
+      ! first process
+      if (myid.eq.0) then
+
+          DO jj=1,nd,nb
+               DO j=jj,jj+nb-1
+                sum(j)=w(1)*r(j,1)
+                r(j,1)=r(j,1)*l(1)
+               ENDDO
+
+              DO i=2,n
+               DO j=jj,jj+nb-1
+                r(j,i)=r(j,i)-r(j,i-1)
+                sum(j)=sum(j)+w(i)*r(j,i)
+                r(j,i)=r(j,i)*l(i)
+               ENDDO
+              ENDDO
+
+              buf(:,1) = r(jj:jj+nb-1,n)
+              buf(:,2) = sum(jj:jj+nb-1)
+              CALL MPI_SSEND(buf, 2*nb, MPI_REAL8, pid+1, 1, MPI_COMM_TASK, mpi_err)
+
+          ENDDO
+
+          DO jj=1,nd,nb
+
+              CALL MPI_RECV(buf, 2*nb, MPI_REAL8, pid+1, 2, MPI_COMM_TASK, mpi_status, mpi_err)
+              r0 = buf(:,1)
+              r1 = buf(:,2)
+
+              i=n
+               DO j=jj,jj+nb-1
+                r(j,i)=r(j,i)-l(i)*r0(j-jj+1)-w(i)*r1(j-jj+1)
+               ENDDO
+              DO i=n-1,1,-1
+               DO j=jj,jj+nb-1
+                r(j,i)=r(j,i)-l(i)*r(j,i+1)-w(i)*r1(j-jj+1)
+               ENDDO
+              ENDDO
+
+          ENDDO
+
+      ! last process
+      elseif (pid.eq.(np-1)) then
+
+          DO jj=1,nd,nb
+
+              CALL MPI_RECV(buf, 2*nb, MPI_REAL8, pid-1, 1, MPI_COMM_TASK, mpi_status, mpi_err)
+              r0 = buf(:,1)
+              sum0 = buf(:,2)
+
+              i=1
+               DO j=jj,jj+nb-1
+                r(j,i)=r(j,i)-r0(j-jj+1)
+                sum(j)=sum0(j-jj+1)+w(i)*r(j,i)
+                r(j,i)=r(j,i)*l(i)
+               ENDDO
+
+              DO i=2,n-1
+               DO j=jj,jj+nb-1
+                r(j,i)=r(j,i)-r(j,i-1)
+                sum(j)=sum(j)+w(i)*r(j,i)
+                r(j,i)=r(j,i)*l(i)
+               ENDDO
+              ENDDO
+
+               DO j=jj,jj+nb-1
+                r(j,n)=l(n)*(r(j,n)-sum(j))
+               ENDDO
+
+          ENDDO
+
+          DO jj=1,nd,nb
+
+               DO j=jj,jj+nb-1
+                r(j,n-1)=r(j,n-1)-w(n-1)*r(j,n)
+               ENDDO
+
+              DO i=n-2,1,-1
+               DO j=jj,jj+nb-1
+                r(j,i)=r(j,i)-l(i)*r(j,i+1)-w(i)*r(j,n)
+               ENDDO
+              ENDDO
+
+              buf(:,1) = r(jj:jj+nb-1,1)
+              buf(:,2) = r(jj:jj+nb-1,n)
+              CALL MPI_SSEND(buf, 2*nb, MPI_REAL8, pid-1, 2, MPI_COMM_TASK, mpi_err)
+
+          ENDDO
+
+      ! intermediate process
+      else
+
+          DO jj=1,nd,nb
+
+              CALL MPI_RECV(buf, 2*nb, MPI_REAL8, pid-1, 1, MPI_COMM_TASK, mpi_status, mpi_err)
+              r0 = buf(:,1)
+              sum0 = buf(:,2)
+
+              i=1
+               DO j=jj,jj+nb-1
+                r(j,i)=r(j,i)-r0(j-jj+1)
+                sum(j)=sum0(j-jj+1)+w(i)*r(j,i)
+                r(j,i)=r(j,i)*l(i)
+               ENDDO
+
+              DO i=2,n
+               DO j=jj,jj+nb-1
+                r(j,i)=r(j,i)-r(j,i-1)
+                sum(j)=sum(j)+w(i)*r(j,i)
+                r(j,i)=r(j,i)*l(i)
+               ENDDO
+              ENDDO
+
+              buf(:,1) = r(jj:jj+nb-1,n)
+              buf(:,2) = sum(jj:jj+nb-1)
+              CALL MPI_SSEND(buf, 2*nb, MPI_REAL8, pid+1, 1, MPI_COMM_TASK, mpi_err)
+
+          ENDDO
+
+          DO jj=1,nd,nb
+
+              CALL MPI_RECV(buf, 2*nb, MPI_REAL8, pid+1, 2, MPI_COMM_TASK, mpi_status, mpi_err)
+              r0 = buf(:,1)
+              r1 = buf(:,2)
+
+              i=n
+               DO j=jj,jj+nb-1
+                r(j,i)=r(j,i)-l(i)*r0(j-jj+1)-w(i)*r1(j-jj+1)
+               ENDDO
+              DO i=n-1,1,-1
+               DO j=jj,jj+nb-1
+                r(j,i)=r(j,i)-l(i)*r(j,i+1)-w(i)*r1(j-jj+1)
+               ENDDO
+              ENDDO
+
+              buf(:,1) = r(jj:jj+nb-1,1)
+              !buf(:,2) = r1
+              CALL MPI_SSEND(buf, 2*nb, MPI_REAL8, pid-1, 2, MPI_COMM_TASK, mpi_err)
+
+          ENDDO
+
+      endif
+
+     END SUBROUTINE par_ptdslv
+
      SUBROUTINE d2fp(n,h,x,dx,nd,dir)
       INTEGER,INTENT(IN) :: n,nd,dir
       REAL*8,INTENT(IN) :: h
@@ -310,6 +569,119 @@
       IF (dir.eq.3) CALL ptdslv(dx,n,lzs,wzs,nd) ! z-direction
      END SUBROUTINE d2fp
 
+     SUBROUTINE par_d2fp(x,xu,xl,dx,h,nd,n,nall,dir)
+      INTEGER,INTENT(IN) :: nall,n,nd,dir
+      REAL*8,INTENT(IN) :: h
+      REAL*8,INTENT(IN),DIMENSION(nd,n) :: x
+      REAL*8,INTENT(OUT),DIMENSION(nd,n) :: dx
+      INTEGER :: i,j
+      REAL*8 :: h2,r1,r2,t1,t2
+
+      REAL*8,DIMENSION(nd,2) :: xu, xl
+
+      h2=1./(h*h)
+      r1=6.
+      r2=3./8.
+
+      !if (myid.eq.0) write(*,*) "parallel d2fp commuication"
+      if (myid.eq.master) then
+
+              call MPI_ISEND (x(1,1),   2*nd, MPI_REAL8, numprocs-1, myid, &
+                              MPI_COMM_TASK, mpi_request, mpi_err)
+
+              call MPI_RECV  (xu,       2*nd, MPI_REAL8, myid+1,     myid+1, &
+                              MPI_COMM_TASK, mpi_status, mpi_err)
+
+              call MPI_ISEND (x(1,n-1), 2*nd, MPI_REAL8, myid+1,     myid, &
+                              MPI_COMM_TASK, mpi_request, mpi_err)
+
+              call MPI_RECV  (xl,       2*nd, MPI_REAL8, numprocs-1, numprocs-1, &
+                              MPI_COMM_TASK, mpi_status, mpi_err)
+
+      elseif (myid.eq.numprocs-1) then
+
+              call MPI_ISEND (x(1,1),   2*nd, MPI_REAL8, myid-1, myid, &
+                              MPI_COMM_TASK, mpi_request, mpi_err)
+
+              call MPI_RECV  (xu,       2*nd, MPI_REAL8, 0,      0, &
+                              MPI_COMM_TASK, mpi_status, mpi_err)
+
+              call MPI_ISEND (x(1,n-1), 2*nd, MPI_REAL8, 0,      myid, &
+                              MPI_COMM_TASK, mpi_request, mpi_err)
+
+              call MPI_RECV  (xl,       2*nd, MPI_REAL8, myid-1, myid-1, &
+                              MPI_COMM_TASK, mpi_status, mpi_err)
+
+      else
+
+              call MPI_ISEND (x(1,1),   2*nd, MPI_REAL8, myid-1, myid, &
+                              MPI_COMM_TASK, mpi_request, mpi_err)
+
+              call MPI_RECV  (xu,       2*nd, MPI_REAL8, myid+1, myid+1, &
+                              MPI_COMM_TASK, mpi_status, mpi_err)
+
+              call MPI_ISEND (x(1,n-1), 2*nd, MPI_REAL8, myid+1, myid, &
+                              MPI_COMM_TASK, mpi_request, mpi_err)
+
+              call MPI_RECV  (xl,       2*nd, MPI_REAL8, myid-1, myid-1, &
+                              MPI_COMM_TASK, mpi_status, mpi_err)
+
+      endif
+
+
+      DO j=1,nd
+
+        t1 = ( x(j,n)-2.*x(j,n-1)+x(j,n-2))
+        t2 = (xu(j,1)-2.*x(j,n-1)+x(j,n-3))
+        IF ( x(j,n).eq.x(j,n-1).and.x(j,n-1).eq.x(j,n-2)) t1=0.
+        IF (xu(j,1).eq.x(j,n-1).and.x(j,n-1).eq.x(j,n-3)) t2=0.
+        dx(j,n-1)=(r1*t1+r2*t2)*h2
+
+        t1 = (xu(j,1)-2.*x(j,n)+x(j,n-1))
+        t2 = (xu(j,2)-2.*x(j,n)+x(j,n-2))
+        IF (xu(j,1).eq.x(j,n).and.x(j,n).eq.x(j,n-1)) t1=0.
+        IF (xu(j,2).eq.x(j,n).and.x(j,n).eq.x(j,n-2)) t2=0.
+!        dx(j,n)=(r1*(x(j,1)-2.*x(j,n)+x(j,n-1)) &
+!               +r2*(x(j,2)-2.*x(j,n)+x(j,n-2)))
+        dx(j,n)=(r1*t1+r2*t2)*h2
+
+        t1 = (x(j,2)-2.*x(j,1)+xl(j,2))
+        t2 = (x(j,3)-2.*x(j,1)+xl(j,1))
+        IF (x(j,2).eq.x(j,1).and.x(j,1).eq.xl(j,2)) t1=0.
+        IF (x(j,3).eq.x(j,1).and.x(j,1).eq.xl(j,1)) t2=0.
+!        dx(j,1)=(r1*(x(j,2)-2.*x(j,1)+x(j,n)) &
+!               +r2*(x(j,3)-2.*x(j,1)+x(j,n-1)))
+        dx(j,1)=(r1*t1+r2*t2)*h2
+
+        t1 = (x(j,3)-2.*x(j,2)+ x(j,1))
+        t2 = (x(j,4)-2.*x(j,2)+xl(j,2))
+        IF (x(j,3).eq.x(j,2).and.x(j,2).eq.x(j,1)) t1=0.
+        IF (x(j,4).eq.x(j,2).and.x(j,2).eq.xl(j,2)) t2=0.
+!        dx(j,2)=(r1*(x(j,3)-2.*x(j,2)+x(j,1)) &
+!               +r2*(x(j,4)-2.*x(j,2)+x(j,n)))
+        dx(j,2)=(r1*t1+r2*t2)*h2
+
+      ENDDO
+
+      DO i=3,n-2
+       DO j=1,nd
+
+        t1 = (x(j,i+1)-2.*x(j,i)+x(j,i-1))
+        t2 = (x(j,i+2)-2.*x(j,i)+x(j,i-2))
+        IF (x(j,i+1).eq.x(j,i).and.x(j,i).eq.x(j,i-1)) t1=0.
+        IF (x(j,i+2).eq.x(j,i).and.x(j,i).eq.x(j,i-2)) t2=0.
+!        dx(j,i)=(r1*(x(j,i+1)-2.*x(j,i)+x(j,i-1)) &
+!               +r2*(x(j,i+2)-2.*x(j,i)+x(j,i-2)))
+        dx(j,i)=(r1*t1+r2*t2)*h2
+
+       ENDDO
+      ENDDO
+
+      IF (dir.eq.1) CALL par_ptdslv(dx,lxs,wxs,nd,n,nall) ! x-direction
+      IF (dir.eq.2) CALL par_ptdslv(dx,lys,wys,nd,n,nall) ! y-direction
+      IF (dir.eq.3) CALL par_ptdslv(dx,lzs,wzs,nd,n,nall) ! z-direction
+     END SUBROUTINE par_d2fp
+
      SUBROUTINE tdslv(r,n,l,nd)
       INTEGER,INTENT(IN) :: n,nd
       REAL*8,INTENT(INOUT),DIMENSION(nd,n) :: r

m_fdm_calc.f90

@@ -1,24 +1,34 @@
 module m_fdm_calc
 
   use m_parameters
+  use m_io
   use m_compact
 
       implicit none
 
 !variables
 
+      integer, parameter :: nscalar = 2
+      integer, parameter :: neq = 1 + nscalar
+      integer, parameter :: nd1 = 2 + 2 * nscalar
+      integer, parameter :: nd2 = nscalar
+
       real*8, dimension(:,:,:), allocatable :: u_,v_,w_
       real*8, dimension(:,:,:,:), allocatable :: y1,y2,yf
+      real*8, dimension(:,:), allocatable :: fz, dfz, fzz, dfzz
+      real*8, dimension(:,:), allocatable :: fzl, fzu, fzzl, fzzu
  
       real*8 :: in_yr,out_yr,refwr,minf
       integer :: fullsavenum   !,svfx,svfy
       integer :: fdmcyc,fdmsavecount
       real*8 :: fdmtime,fdmdt
       real*8 :: sumc,oldsumc,time_int,sum_wrate
+      real*8 :: sumc1,sum_wrate1
       real*8 :: visdt, convdt
       integer :: fdmstep
       real*8 :: coe,tt1,tt2,tt3,tt4
       real*8 :: umax,umin,vmax,vmin,wmax,wmin  ! J. Kwon
+      real*8, dimension(3) :: velmax, velmin, velmax1, velmin1
 
 !===========================================================================
 !===========================================================================
@@ -26,134 +36,28 @@ module m_fdm_calc
   contains
 
 
-  subroutine fdm_restart_write
-      write(*,*) '======================================================='
-      write(*,*) 'Full results are being written',fullsavenum
-      OPEN (fullsavenum,form='unformatted',status='unknown')
-      write (fullsavenum) fdmtime,nx,ny,nz_all,oldsumc,time_int
-      write (fullsavenum) fdmcyc,DT,dummyu_
-      write (fullsavenum) dt_fdmsave,dt_fullsave
-      write (fullsavenum) t_fdmsave,t_fullsave
-      write (fullsavenum) in_yr,out_yr
-      write (fullsavenum) u_+dummyu_,v_,w_,y1
-      CLOSE (fullsavenum)
-      write(*,*) '======================================================='
-      fullsavenum=fullsavenum+1
-  end subroutine fdm_restart_write
+  subroutine fdm_restart_read
+    integer(kind=MPI_INTEGER_KIND) :: fh
+    integer(kind=MPI_OFFSET_KIND)  :: offset
 
-  subroutine prepare_fdm
+      integer  :: header_length
 
-      implicit none
+    write(fname,"('fort.',i4)") restartnum
 
-      integer :: i,j,ii,k
-      real*8 :: fl_location,wrate,yr,c
-      real*8 :: tar_sumc
-
-      allocate(u_(nx,ny,nz_all))
-      allocate(v_(nx,ny,nz_all))
-      allocate(w_(nx,ny,nz_all))
-
-      u_=0.0
-      v_=0.0
-      w_=0.0
-
-      if (fdm_sw .eq. 0) then
-         return
-      endif
-
-! DQ initializing
-      fdmsavecount=1        !FDM save count
-      sum_wrate=0.
-      sumc=0.       !for adjusting mean velocity
-      oldsumc=0.
-      time_int=0.
-      visdt=9999.
-      convdt=9999.
-      fdmcyc=0
-      fdmtime=0.
-      fullsavenum=1000      !full save file
-
-      allocate(y1(2,nx,ny,nz_all))
-      allocate(y2(2,nx,ny,nz_all))
-      allocate(yf(2,nx,ny,nz_all))
-
-      CALL ludcmp(nx,ny,nz_all,1,0,0)
-
-      refwr=pre*1.*exp(-ac/(1.+bc*c_ref))  ! Kwon
-
-
-
-!FDM normal start==================================
-      if (restartnum==0) then
-
-       OPEN(305,FILE='sfield.bin',form='unformatted',status='unknown')
-         DO i=1,nx
-           READ (305) y1(2,i,1,1)   ! Yr
-         ENDDO
-       CLOSE (305)
-
-       in_yr=y1(2,1,1,1)         ! inlet_Yr
-       out_yr=y1(2,nx,1,1)       ! outlet_Yr
- 
-       do i=1,ny 
-        do j=1,nz_all 
-         do ii=1,nx 
-           y1(1,ii,i,j)=1.            ! rho initializing
-           y1(2,ii,i,j)=y1(2,ii,1,1)  ! Yr  initializing
-         enddo 
-        enddo 
-       enddo 
- 
-
-! flame location setup
-         sum_wrate=0.; sumc=0.
-         DO ii=1,nz_all
-          DO j=1,ny
-           DO i=1,nx
-            yr=y1(2,i,j,ii)/y1(1,i,j,ii)
-            c=1.-yr
-            IF (c.lt.0.) c=0.
-            wrate=pre*yr*exp(-ac/(1.+bc*c))
-            ! cold boundary difficulty treatment
-            IF(c.le.c_ref) THEN
-               wrate=min_wr
-               IF(c.gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(c-c_ref))+ &
-                  min_wr-refwr*exp(prof_wr*(c_cut-c_ref)))/(1.-exp(prof_wr*(c_cut-c_ref)))
-            ENDIF
-            sum_wrate=sum_wrate+wrate*(hx*hy*hy)
-            sumc=sumc+(1.-yr)
-           ENDDO
-          ENDDO
-         ENDDO
-
-         sum_wrate=sum_wrate/(hy*hy*ny*nz_all)
-         write(*,633) sum_wrate
-
-         fl_location=(hx*(nx-1.))*(1.-(sumc/(nx*ny*nz_all)))
-         sumc=sumc*(hx*hy*hy)/(hy*(ny-1.)*hy*(ny-1.))
-         oldsumc=sumc
-         write(*,634) fl_location/(REAL(nx-1)*hx)*100.
-
-633     format (' ** Consumption Speed, Sc = ',f7.4)
-634     format (' ** Flame Location = ',f7.3, ' % point of x-domain.')
-
-! Restart setup ================================================================
-      else
+    if (myid.eq.master) then
          write (6,*) '********************************************************'
          write (6,*) ' FDM field is being initialized for restarting.'
          write (6,457) restartnum
  457     format('  Restart file number : ',i6)
 
-         OPEN (restartnum,form='unformatted',status='unknown')
+         OPEN (restartnum,file=fname,form='unformatted',status='unknown')
          read (restartnum) fdmtime,tt1,tt2,tt3,oldsumc,time_int
-         read (restartnum) fdmcyc,tt4,dummyu_
+         read (restartnum) fdmcyc,fdmdt,dummyu_
          read (restartnum) tt1,tt2    !dt_fdmsave,dt_fullsave
          read (restartnum) t_fdmsave,t_fullsave
          read (restartnum) in_yr,out_yr
-         read (restartnum) u_,v_,w_,y1
          CLOSE (restartnum)
 
-
          if(fdmtime.ge.t_fdmsave) t_fdmsave=t_fdmsave+dt_fdmsave
          if(fdmtime.ge.t_fullsave) t_fullsave=t_fullsave+dt_fullsave
 
@@ -176,6 +80,264 @@ module m_fdm_calc
  455      format('  !! Spectral Time : ',f10.5,' / FDM Time : ',f10.5)
          endif
 
+    end if
+
+    count = 1
+
+    call MPI_BCAST(fdmtime, count, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err)
+    call MPI_BCAST(tt1, count, MPI_INTEGER8, 0, MPI_COMM_TASK, mpi_err)
+    call MPI_BCAST(tt2, count, MPI_INTEGER8, 0, MPI_COMM_TASK, mpi_err)
+    call MPI_BCAST(tt3, count, MPI_INTEGER8, 0, MPI_COMM_TASK, mpi_err)
+    call MPI_BCAST(oldsumc, count, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err)
+    call MPI_BCAST(time_int, count, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err)
+
+    call MPI_BCAST(fdmcyc, count, MPI_INTEGER8, 0, MPI_COMM_TASK, mpi_err)
+    call MPI_BCAST(fdmdt, count, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err)
+    call MPI_BCAST(dummyu_, count, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err)
+
+    !call MPI_BCAST(dt_fdmsave, count, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err)
+    !call MPI_BCAST(dt_fullsave, count, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err)
+
+    !call MPI_BCAST(t_fdmsave, count, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err)
+    !call MPI_BCAST(t_fullsave, count, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err)
+
+    call MPI_BCAST(in_yr, count, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err)
+    call MPI_BCAST(out_yr, count, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err)
+
+
+  ! opening the file
+    call MPI_INFO_CREATE(mpi_info, mpi_err)
+    call MPI_FILE_OPEN(MPI_COMM_TASK,fname,MPI_MODE_RDONLY,mpi_info,fh,mpi_err)
+
+    header_length = 2*4+6*8 + 2*4+3*8 + 3*(2*4+2*8)
+
+    count = nx * ny * nz
+    offset = header_length + 4 + 0*nx*ny*nz_all*8 + myid*nx*ny*nz*8
+    call MPI_FILE_READ_AT_ALL(fh, offset, u_, count, MPI_REAL8, mpi_status, mpi_err)
+
+    offset = header_length + 4 + 1*nx*ny*nz_all*8 + myid*nx*ny*nz*8
+    call MPI_FILE_READ_AT_ALL(fh, offset, v_, count, MPI_REAL8, mpi_status, mpi_err)
+
+    offset = header_length + 4 + 2*nx*ny*nz_all*8 + myid*nx*ny*nz*8
+    call MPI_FILE_READ_AT_ALL(fh, offset, w_, count, MPI_REAL8, mpi_status, mpi_err)
+
+    count = neq * nx * ny * nz
+    offset = header_length + 4 + 3*nx*ny*nz_all*8 + myid*neq*nx*ny*nz*8
+    call MPI_FILE_READ_AT_ALL(fh, offset, y1, count, MPI_REAL8, mpi_status, mpi_err)
+
+    call MPI_FILE_CLOSE(fh, mpi_err)
+    call MPI_INFO_FREE(mpi_info, mpi_err)
+
+  end subroutine fdm_restart_read
+
+  subroutine fdm_restart_write
+      integer(kind=MPI_INTEGER_KIND) :: fh
+      integer(kind=MPI_OFFSET_KIND)  :: offset
+
+      integer  :: header_length
+  ! --------------- writing process ------------------
+
+      if (myid.eq.0) write(*,*) '======================================================='
+      if (myid.eq.0) write(*,*) 'Full results are being written',fullsavenum
+
+  !fname = 'fort.'//savenum
+      write(fname,"('fort.',i4)") fullsavenum
+
+  ! opening the file
+      call MPI_INFO_CREATE(mpi_info, mpi_err)
+      call MPI_FILE_OPEN(MPI_COMM_TASK,fname,MPI_MODE_WRONLY+MPI_MODE_CREATE,mpi_info,fh,mpi_err)
+
+  ! the master node writes the header with parameters
+      if (myid.eq.0) then
+         count = 1
+         call MPI_FILE_WRITE(fh,   8*6, count, MPI_INTEGER4, mpi_status, mpi_err)
+         call MPI_FILE_WRITE(fh,   fdmtime, count, MPI_REAL8, mpi_status, mpi_err)
+         call MPI_FILE_WRITE(fh,   nx, count, MPI_INTEGER8, mpi_status, mpi_err)
+         call MPI_FILE_WRITE(fh,   ny, count, MPI_INTEGER8, mpi_status, mpi_err)
+         call MPI_FILE_WRITE(fh,   nz_all, count, MPI_INTEGER8, mpi_status, mpi_err)
+         call MPI_FILE_WRITE(fh,   oldsumc, count, MPI_REAL8, mpi_status, mpi_err)
+         call MPI_FILE_WRITE(fh,   time_int, count, MPI_REAL8, mpi_status, mpi_err)
+         call MPI_FILE_WRITE(fh,   8*6, count, MPI_INTEGER4, mpi_status, mpi_err)
+
+         call MPI_FILE_WRITE(fh,   8*3, count, MPI_INTEGER4, mpi_status, mpi_err)
+         call MPI_FILE_WRITE(fh,   fdmcyc, count, MPI_INTEGER8, mpi_status, mpi_err)
+         call MPI_FILE_WRITE(fh,   DT, count, MPI_REAL8, mpi_status, mpi_err)
+         call MPI_FILE_WRITE(fh,   dummyu_, count, MPI_REAL8, mpi_status, mpi_err)
+         call MPI_FILE_WRITE(fh,   8*3, count, MPI_INTEGER4, mpi_status, mpi_err)
+
+         call MPI_FILE_WRITE(fh,   8*2, count, MPI_INTEGER4, mpi_status, mpi_err)
+         call MPI_FILE_WRITE(fh,   dt_fdmsave, count, MPI_REAL8, mpi_status, mpi_err)
+         call MPI_FILE_WRITE(fh,   dt_fullsave, count, MPI_REAL8, mpi_status, mpi_err)
+         call MPI_FILE_WRITE(fh,   8*2, count, MPI_INTEGER4, mpi_status, mpi_err)
+
+         call MPI_FILE_WRITE(fh,   8*2, count, MPI_INTEGER4, mpi_status, mpi_err)
+         call MPI_FILE_WRITE(fh,   t_fdmsave, count, MPI_REAL8, mpi_status, mpi_err)
+         call MPI_FILE_WRITE(fh,   t_fullsave, count, MPI_REAL8, mpi_status, mpi_err)
+         call MPI_FILE_WRITE(fh,   8*2, count, MPI_INTEGER4, mpi_status, mpi_err)
+
+         call MPI_FILE_WRITE(fh,   8*2, count, MPI_INTEGER4, mpi_status, mpi_err)
+         call MPI_FILE_WRITE(fh,   in_yr, count, MPI_REAL8, mpi_status, mpi_err)
+         call MPI_FILE_WRITE(fh,   out_yr, count, MPI_REAL8, mpi_status, mpi_err)
+         call MPI_FILE_WRITE(fh,   8*2, count, MPI_INTEGER4, mpi_status, mpi_err)
+
+         call MPI_FILE_WRITE(fh,   8*nx*ny*nz_all*5, count, MPI_INTEGER4, mpi_status, mpi_err)
+
+         offset = 2*4+6*8 + 2*4+3*8 + 3*(2*4+2*8) + 4 + 5*nx*ny*nz_all*8
+         call MPI_FILE_WRITE_AT(fh, offset, 8*nx*ny*nz_all*5, count, MPI_INTEGER4, mpi_status, mpi_err)
+
+      end if
+
+
+      header_length = 2*4+6*8 + 2*4+3*8 + 3*(2*4+2*8)
+
+      count = nx * ny * nz
+      offset = header_length + 4 + 0*nx*ny*nz_all*8 + myid*nx*ny*nz*8
+      call MPI_FILE_WRITE_AT(fh, offset, u_, count, MPI_REAL8, mpi_status, mpi_err)
+
+      offset = header_length + 4 + 1*nx*ny*nz_all*8 + myid*nx*ny*nz*8
+      call MPI_FILE_WRITE_AT(fh, offset, v_, count, MPI_REAL8, mpi_status, mpi_err)
+
+      offset = header_length + 4 + 2*nx*ny*nz_all*8 + myid*nx*ny*nz*8
+      call MPI_FILE_WRITE_AT(fh, offset, w_, count, MPI_REAL8, mpi_status, mpi_err)
+
+      count = neq * nx * ny * nz
+      offset = header_length + 4 + 3*nx*ny*nz_all*8 + myid*neq*nx*ny*nz*8
+      call MPI_FILE_WRITE_AT(fh, offset, y1, count, MPI_REAL8, mpi_status, mpi_err)
+
+
+      call MPI_FILE_CLOSE(fh, mpi_err)
+      call MPI_INFO_FREE(mpi_info, mpi_err)
+
+      if (myid.eq.0) write(*,*) '======================================================='
+      fullsavenum=fullsavenum+1
+  end subroutine fdm_restart_write
+
+  subroutine prepare_fdm
+
+      implicit none
+   
+      integer :: i,j,ii,k
+      real*8 :: fl_location,wrate,yr,c
+      real*8 :: tar_sumc
+      real*8 :: yy(nx)
+
+      if (fdm_sw .eq. 0) then
+         return
+      endif
+
+      allocate(u_(nx,ny,nz))
+      allocate(v_(nx,ny,nz))
+      allocate(w_(nx,ny,nz))
+
+      u_=0.0
+      v_=0.0
+      w_=0.0
+
+! DQ initializing
+      fdmsavecount=1        !FDM save count
+      sum_wrate=0.
+      sum_wrate1=0.
+      sumc=0.       !for adjusting mean velocity
+      sumc1=0.       !for adjusting mean velocity
+      oldsumc=0.
+      time_int=0.
+      visdt=9999.
+      convdt=9999.
+      fdmcyc=0
+      fdmtime=0.
+      fullsavenum=1000      !full save file
+
+      allocate(y1(neq,nx,ny,nz))
+      allocate(y2(neq,nx,ny,nz))
+      allocate(yf(neq,nx,ny,nz))
+
+      allocate( fz (nd1*nx*ny,nz))
+      allocate(dfz (nd1*nx*ny,nz))
+      allocate( fzz(nd2*nx*ny,nz))
+      allocate(dfzz(nd2*nx*ny,nz))
+
+      allocate(fzu(nd1*nx*ny,2))
+      allocate(fzl(nd1*nx*ny,2))
+
+      allocate(fzzu(nd2*nx*ny,2))
+      allocate(fzzl(nd2*nx*ny,2))
+
+      y1=0.0
+      y2=0.0
+      yf=0.0
+
+      CALL ludcmp(nx,ny,nz_all,1,0,0)
+
+      refwr=pre*1.*exp(-ac/(1.+bc*c_ref))  ! Kwon
+
+
+
+!FDM normal start==================================
+      if (restartnum==0) then
+
+       if (myid.eq.master) then
+       OPEN(305,FILE='sfield.bin',form='unformatted',status='unknown')
+         DO i=1,nx
+           READ (305) yy(i)   ! Yr
+         ENDDO
+       CLOSE (305)
+       end if
+
+       call MPI_BCAST(yy, nx, MPI_REAL8, 0, MPI_COMM_TASK, mpi_err)
+
+       in_yr=yy(1)         ! inlet_Yr
+       out_yr=yy(nx)       ! outlet_Yr
+ 
+       do j=1,nz
+        do i=1,ny 
+         do ii=1,nx 
+           y1(1,ii,i,j)=1.            ! rho initializing
+           y1(2,ii,i,j)=yy(ii)  ! Yr  initializing
+           y1(3,ii,i,j)=0.0     ! Yr  initializing
+         enddo 
+        enddo 
+       enddo 
+ 
+
+! flame location setup
+         sum_wrate=0.; sumc=0.
+         DO ii=1,nz
+          DO j=1,ny
+           DO i=1,nx
+            yr=y1(2,i,j,ii)/y1(1,i,j,ii)
+            c=1.-yr
+            IF (c.lt.0.) c=0.
+            wrate=pre*yr*exp(-ac/(1.+bc*c))
+            ! cold boundary difficulty treatment
+            IF(c.le.c_ref) THEN
+               wrate=min_wr
+               IF(c.gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(c-c_ref))+ &
+                  min_wr-refwr*exp(prof_wr*(c_cut-c_ref)))/(1.-exp(prof_wr*(c_cut-c_ref)))
+            ENDIF
+            sum_wrate=sum_wrate+wrate*(hx*hy*hy)
+            sumc=sumc+(1.-yr)
+           ENDDO
+          ENDDO
+         ENDDO
+
+         sum_wrate1=sum_wrate/(hy*hy*ny*nz_all)
+         CALL MPI_ALLREDUCE(sum_wrate1,sum_wrate,1,MPI_REAL8,MPI_SUM,MPI_COMM_TASK,mpi_err)
+         if (myid.eq.0) write(*,633) sum_wrate
+
+         CALL MPI_ALLREDUCE(sumc,sumc1,1,MPI_REAL8,MPI_SUM,MPI_COMM_TASK,mpi_err)
+         fl_location=(hx*(nx-1.))*(1.-(sumc1/(nx*ny*nz_all)))
+         sumc=sumc1*(hx*hy*hy)/(hy*(ny-1.)*hy*(ny-1.))
+         !CALL MPI_ALLREDUCE(sumc1,sumc,1,MPI_REAL8,MPI_SUM,MPI_COMM_TASK,mpi_err)
+         oldsumc=sumc
+         if (myid.eq.0) write(*,634) fl_location/(REAL(nx-1)*hx)*100.
+
+633     format (' ** Consumption Speed, Sc = ',f7.4)
+634     format (' ** Flame Location = ',f7.3, ' % point of x-domain.')
+
+! Restart setup ================================================================
+      else
+
+        CALL fdm_restart_read
         fullsavenum=restartnum+1
 
       endif
@@ -193,9 +355,12 @@ module m_fdm_calc
       real*8 :: fl_location,delu,wrate,yr
       real*8 :: c
       real*8 :: uvel,vvel,wvel  ! J. Kwon
+
+
+
       fdmcyc=itime
 ! Mean velocty setup  
-      do k=1,nz_all
+      do k=1,nz
        do i=1,nx
         do j=1,ny
           u_(i,j,k)=u_(i,j,k)+dummyu_
@@ -203,32 +368,34 @@ module m_fdm_calc
        enddo
       enddo
 
-       write(*,*)
-       write(*,*) '=========================================================='
+       if (myid.eq.0) write(*,*)
+       if (myid.eq.0) write(*,*) '=========================================================='
 
-      WRITE(*,932) TIME,DT
+      if (myid.eq.0) WRITE(*,932) TIME,DT
 932   format(' Spectral results at time = ',f10.5,', dT = ',f7.5)
-      write(*,933) dummyu_
+      if (myid.eq.0) write(*,933) dummyu_
 933   format(' ** Mean U = ',f7.4)
 
 !----------------------------------------------------------------
 ! FDM DNS BGN
 !----------------------------------------------------------------
-          umax=-1.0d300;umin=1.0d300;vmax=-1.0d300;vmin=1.0d300;wmax=-1.0d300;wmin=1.0d300
+          velmax1(1)=maxval(u_)
+          velmin1(1)=minval(u_)
+          velmax1(2)=maxval(v_)
+          velmin1(2)=minval(v_)
+          velmax1(3)=maxval(w_)
+          velmin1(3)=minval(w_)
 
-!$omp parallel do private(i,j,k) reduction(max:umax,vmax,wmax) reduction(min:umin,vmin,wmin)  NUM_THREADS(6) schedule(static)
-          do k=1,nz_all
-          do j=1,ny
-          do i=1,nx
-            umax=max(umax,u_(i,j,k))
-            umin=min(umin,u_(i,j,k))
-            vmax=max(vmax,v_(i,j,k))
-            vmin=min(vmin,v_(i,j,k))
-            wmax=max(wmax,w_(i,j,k))
-            wmin=min(wmin,w_(i,j,k))
-          enddo
-          enddo
-          enddo
+          ! MPI_REDUCE [uvw](max|min)
+          CALL MPI_ALLREDUCE(velmax1,velmax,3,MPI_REAL8,MPI_MAX,MPI_COMM_TASK,mpi_err)
+          CALL MPI_ALLREDUCE(velmin1,velmin,3,MPI_REAL8,MPI_MIN,MPI_COMM_TASK,mpi_err)
+
+          umax=velmax(1)
+          umin=velmin(1)
+          vmax=velmax(2)
+          vmin=velmin(2)
+          wmax=velmax(3)
+          wmin=velmin(3)
 
           uvel = max(abs(umax), abs(umin))
           vvel = max(abs(vmax), abs(vmin))
@@ -241,12 +408,14 @@ module m_fdm_calc
 
         fdmdt=min(DT,visdt,convdt)        
         DT=fdmdt
-        call solve(nx,ny,nz_all,u_,v_,w_,y1,y2,yf)
+        call solve(nx,ny,nz,u_,v_,w_,y1,y2,yf)
         fdmtime=time+fdmdt
-        write(*,'(a30,3x,4f12.7)')' ** DT, visdt, convdt, fdmdt =' , DT,visdt,convdt,fdmdt
+        if (myid.eq.0) write(*,'(a30,3x,4f12.7)')' ** DT, visdt, convdt, fdmdt =' , DT,visdt,convdt,fdmdt
 
-!$omp parallel do private(i,j,k,yr,wrate) reduction(+:sum_wrate,sumc) NUM_THREADS(6) schedule(static)
-        do k=1,nz_all
+        sumc=0.
+        sum_wrate=0.
+
+        do k=1,nz
         do j=1,ny
         do i=1,nx
            yr=y1(2,i,j,k)/y1(1,i,j,k)
@@ -263,21 +432,25 @@ module m_fdm_calc
         enddo
         enddo
 
+
 ! Real time results for Sc and flame location.
-        sum_wrate=sum_wrate/(hy*hy*REAL(ny*nz_all))
-        write(*,633) sum_wrate
+        sum_wrate1=sum_wrate/(hy*hy*REAL(ny*nz_all))
+        CALL MPI_ALLREDUCE(sum_wrate1,sum_wrate,1,MPI_REAL8,MPI_SUM,MPI_COMM_TASK,mpi_err)
+        if (myid.eq.0) write(*,633) sum_wrate
 633     format (' ** Consumption Speed, Sc = ',f7.4)
 
-        fl_location=(hx*REAL(nx-1))*(1.-(sumc/(REAL(nx*ny*nz_all))))
-        sumc=sumc*(hx*hy*hy)/(hy*(REAL(ny)-1.)*hy*(REAL(ny)-1.))
-        write(*,634) fl_location/(REAL(nx-1)*hx)*100.
+        CALL MPI_ALLREDUCE(sumc,sumc1,1,MPI_REAL8,MPI_SUM,MPI_COMM_TASK,mpi_err)
+        fl_location=(hx*REAL(nx-1))*(1.-(sumc1/(REAL(nx*ny*nz_all))))
+        sumc1=sumc*(hx*hy*hy)/(hy*(REAL(ny)-1.)*hy*(REAL(ny)-1.))
+        CALL MPI_ALLREDUCE(sumc1,sumc,1,MPI_REAL8,MPI_SUM,MPI_COMM_TASK,mpi_err)
+        if (myid.eq.0) write(*,634) fl_location/(REAL(nx-1)*hx)*100.
 634     format (' ** Flame Location = ',f7.3, ' % point of x-domain.')
 
 
-        write(St_data,'(f10.5,10e25.15)')fdmtime,dummyu_,sum_wrate,fl_location/(REAL(nx-1)*hx)*100.
+        if (myid.eq.0) write(St_data,'(f10.5,10e25.15)')fdmtime,dummyu_,sum_wrate,fl_location/(REAL(nx-1)*hx)*100.
 
 ! Control the inflow mean velocity, dummyu_
-        write(*,'(a20,2f12.8)') ' ** oldsumc, sumc = ',oldsumc,sumc
+        if (myid.eq.0) write(*,'(a20,2f12.8)') ' ** oldsumc, sumc = ',oldsumc,sumc
         if(swadtv.ne.0.and.mod((fdmcyc+1),swadtv).eq.0) then
           time_int=fdmtime-time_int
           delu=(sumc-oldsumc)/time_int
@@ -291,11 +464,11 @@ module m_fdm_calc
 ! FDM DNS END
 !----------------------------------------------------------------------------
        
-      write(*,*) 'Number of whole calculation = ',fdmcyc
-      write(*,*) '=========================================================='
-      write(*,*)
+      if (myid.eq.0) write(*,*) 'Number of whole calculation = ',fdmcyc
+      if (myid.eq.0) write(*,*) '=========================================================='
+      if (myid.eq.0) write(*,*)
 
-       call write_vel_max ( umax,vmax,wmax,umin,vmin,wmin )
+      if (myid.eq.0) call write_vel_max ( umax,vmax,wmax,umin,vmin,wmin )
 
      return 
 
@@ -307,85 +480,63 @@ module m_fdm_calc
       implicit none
 
       integer :: i,j,k,xx,yy,zz,ii
+      integer :: n
+      integer :: idx1, idx2, idx3
       real*8 :: wrate,yr,yp
-      real*8 :: r1_(2,xx,yy,zz),f_(2,xx,yy,zz)
+      real*8 :: r1_(neq,xx,yy,zz),f_(neq,xx,yy,zz)
       real*8 :: uu_(xx,yy,zz),vv_(xx,yy,zz),ww_(xx,yy,zz)
-      real*8 :: ux(4,xx),dux(4,xx),d2ux(xx)
-      real*8 :: uy(4,yy),duy(4,yy),d2uy(yy)
-      real*8 :: uz(4,yy),duz(4,yy),d2uz(yy)
-      real*8 :: Ly(yy,zz),Dy
+      real*8 :: ux(nd1,xx),dux(nd1,xx),d2ux(nd2,xx)
+      real*8 :: uy(nd1,yy),duy(nd1,yy),d2uy(nd2,yy)
+      real*8 :: uz(nd1,zz),duz(nd1,zz),d2uz(nd2,zz)
+      real*8 :: uux(nd2,xx)
+      real*8 :: uuy(nd2,yy)
+      real*8 :: uuz(nd2,zz)
+      real*8 :: y
 
-! x-direction
-!$omp parallel do private(i,j,k,ux,dux,d2ux,wrate,uy,duy,d2uy) NUM_THREADS(6) schedule(static)
+
+! reaction source term
       DO k=1,zz
        DO j=1,yy
         DO i=1,xx
-         ux(1,i)=r1_(1,i,j,k)                            ! 1:rho
-         ux(2,i)=r1_(2,i,j,k)/r1_(1,i,j,k)               ! 2:Y
-         ux(3,i)=ux(1,i)*uu_(i,j,k)                      ! 3:rho*u
-         ux(4,i)=ux(3,i)*ux(2,i)                         ! 4:rho*u*Y
-        ENDDO
 
-         CALL dfnonp(xx,hx,ux(1:4,:),dux(1:4,:),4,1)
-         CALL d2fnonp(xx,hx,ux(2,:),d2ux(:),1,1)
+!        y=r1_(2,i,j,k)/r1_(1,i,j,k)               ! 2:Y
 
-        DO i=1,xx
-         wrate=pre*ux(2,i)*exp(-ac/(1.+bc*(1.-ux(2,i))))   !wrate
-         IF ((1.-ux(2,i)).le.c_ref) THEN
-           wrate=min_wr
-          IF ((1.-ux(2,i)).gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(1.-ux(2,i)-c_ref))+ &
-          min_wr-refwr*exp(prof_wr*(c_cut-c_ref)))/(1.-exp(prof_wr*(c_cut-c_ref)))
-         ENDIF
+         wrate=0.0
+!        wrate=pre*y*exp(-ac/(1.+bc*(1.-y)))   !wrate
+!        IF ((1.-y).le.c_ref) THEN
+!          wrate=min_wr
+!         IF ((1.-y).gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(1.-y-c_ref))+ &
+!         min_wr-refwr*exp(prof_wr*(c_cut-c_ref)))/(1.-exp(prof_wr*(c_cut-c_ref)))
+!        ENDIF
 
-
-! -( d(rho*u)/dx )
-         f_(1,i,j,k)=-dux(3,i)   ! continuity
-
-! -( d(rho*u*Yr)/dx ) + d(rho*D* d(Yr)/dx)/dx
-! = -( d(rho*u*Yr)/dx )
-!    + D* (rho* d2(Yr)/dx2 + d(rho)/dx * d(Yr)/dx )
-         f_(2,i,j,k)=-dux(4,i) + diff*(ux(1,i)*d2ux(i)+dux(1,i)*dux(2,i)) - wrate  ! species conservation
+         f_(1,i,j,k) = 0.0 ! continuity
+         f_(2,i,j,k) = - wrate  ! species conservation
+         f_(3,i,j,k) = - wrate  ! species conservation
 
         ENDDO
        ENDDO
-
-!! y-direction
-       DO i=1,xx
-        DO j=1,yy
-         uy(1,j)=r1_(1,i,j,k)                            ! 1:rho
-         uy(2,j)=r1_(2,i,j,k)/r1_(1,i,j,k)               ! 2:Y
-         uy(3,j)=uy(1,j)*vv_(i,j,k)                       ! 3:rho*v
-         uy(4,j)=uy(3,j)*uy(2,j)                         ! 4:rho*v*Y
       ENDDO
 
-        CALL dfp(yy,hy,uy(1:4,:),duy(1:4,:),4,2)
-        CALL d2fp(yy,hy,uy(2,:),d2uy(:),1,2)
-
-        DO j=1,yy
-! -( d(rho*v)/dy ) 
-         f_(1,i,j,k)=f_(1,i,j,k)-duy(3,j)   ! continuity
-
-! -( d(rho*v*Yr)/dy ) + d(rho*D* d(Yr)/dy)/dy
-! = -( d(rho*v*Yr)/dy ) 
-!    + D* (rho* d2(Yr)/dyy2 + d(rho)/dy * d(Yr)/dy )
-         f_(2,i,j,k)=f_(2,i,j,k)-duy(4,j) + diff*(uy(1,j)*d2uy(j)+duy(1,j)*duy(2,j))  ! species conserv.
-        ENDDO
-       ENDDO
-      ENDDO
 
 !! z-direction
-!$omp parallel do private(i,j,k,uz,duz,d2uz) NUM_THREADS(6) schedule(static)
+    if (numprocs.eq.1) then
+
       DO j=1,yy
        DO i=1,xx
         DO k=1,zz
          uz(1,k)=r1_(1,i,j,k)                             ! 1:rho
-         uz(2,k)=r1_(2,i,j,k)/r1_(1,i,j,k)               ! 2:Y
+         uz(2,k)=r1_(2,i,j,k)/r1_(1,i,j,k)                ! 2:YF
          uz(3,k)=uz(1,k)*ww_(i,j,k)                       ! 3:rho*w
-         uz(4,k)=uz(3,k)*uz(2,k)                         ! 4:rho*w*Y
+         uz(4,k)=uz(3,k)*uz(2,k)                          ! 4:rho*w*YF
+         uz(5,k)=r1_(3,i,j,k)/r1_(1,i,j,k)                ! 5:YO
+         uz(6,k)=uz(3,k)*uz(5,k)                          ! 6:rho*w*YO
+
+         uuz (1,k)=uz(2,k)
+         uuz (2,k)=uz(5,k)
         ENDDO
 
-       CALL dfp(zz,hy,uz(1:4,:),duz(1:4,:),4,3)
-       CALL d2fp(zz,hy,uz(2,:),d2uz(:),1,3)    
+       CALL dfp (yy,hy, uz(1:nd1,:), duz(1:nd1,:),nd1,3)
+       CALL d2fp(yy,hy,uuz(1:nd2,:),d2uz(1:nd2,:),nd2,3)
 
         DO k=1,zz
 ! -( d(rho*w)/dz )
@@ -394,19 +545,131 @@ module m_fdm_calc
 ! -( d(rho*w*Yr)/dz ) + d(rho*D* d(Yr)/dz)/dz
 ! = -( d(rho*w*Yr)/dz )
 !    + D* (rho* d2(Yr)/dz2 + d(rho)/dz * d(Yr)/dz )
-         f_(2,i,j,k)=f_(2,i,j,k)-duz(4,k) + diff*(uz(1,k)*d2uz(k)+duz(1,k)*duz(2,k))   ! species conserv.
+         f_(2,i,j,k) = f_(2,i,j,k) - duz(4,k) + diff*(uz(1,k)*d2uz(1,k)+duz(1,k)*duz(2,k))   ! species conserv.
+         f_(3,i,j,k) = f_(3,i,j,k) - duz(6,k) + diff*(uz(1,k)*d2uz(2,k)+duz(1,k)*duz(5,k))   ! species conserv.
         ENDDO
        ENDDO
       ENDDO
 
+    else
+
+      DO k=1,zz
+       DO j=1,yy
+        DO i=1,xx
+         idx3 = xx*(j-1)+i
+         idx1 = (idx3-1)*nd1
+         idx2 = (idx3-1)*nd2
+         fz(idx1+1,k) = r1_(1,i,j,k)                             ! 1:rho
+         fz(idx1+2,k) = r1_(2,i,j,k)/r1_(1,i,j,k)                ! 2:YF
+         fz(idx1+3,k) = r1_(1,i,j,k)*ww_(i,j,k)                  ! 3:rho*w
+         fz(idx1+4,k) = r1_(2,i,j,k)*ww_(i,j,k)                  ! 4:rho*w*YF
+         fz(idx1+5,k) = r1_(3,i,j,k)/r1_(1,i,j,k)                ! 5:YO
+         fz(idx1+6,k) = r1_(3,i,j,k)*ww_(i,j,k)                  ! 6:rho*w*YO
+
+         fzz(idx2+1,k) = r1_(2,i,j,k)/r1_(1,i,j,k)
+         fzz(idx2+2,k) = r1_(3,i,j,k)/r1_(1,i,j,k)
+        ENDDO
+       ENDDO
+      ENDDO
+
+       CALL par_dfp (fz,  fzl,  fzu,  dfz,  hy, nd1*xx*yy, zz, yy, 3)
+       CALL par_d2fp(fzz, fzzl, fzzu, dfzz, hy, nd2*xx*yy, zz, yy, 3)
+
+      DO k=1,zz
+       DO j=1,yy
+        DO i=1,xx
+         idx3 = xx*(j-1)+i
+         idx1 = (idx3-1)*nd1
+         idx2 = (idx3-1)*nd2
+! -( d(rho*w)/dz )
+         f_(1,i,j,k) = f_(1,i,j,k) - dfz(idx1+3,k)         ! continuity
+
+! -( d(rho*w*Yr)/dz ) + d(rho*D* d(Yr)/dz)/dz
+! = -( d(rho*w*Yr)/dz )
+!    + D* (rho* d2(Yr)/dz2 + d(rho)/dz * d(Yr)/dz )
+         f_(2,i,j,k) = f_(2,i,j,k) - dfz(idx1+4,k) + diff*(fz(idx1+1,k)*dfzz(idx2+1,k)+dfz(idx1+1,k)*dfz(idx1+2,k))   ! species conserv.
+         f_(3,i,j,k) = f_(3,i,j,k) - dfz(idx1+6,k) + diff*(fz(idx1+1,k)*dfzz(idx2+2,k)+dfz(idx1+1,k)*dfz(idx1+5,k))   ! species conserv.
+
+        ENDDO
+       ENDDO
+      ENDDO
+
+    endif
+
+! x-direction
+      DO k=1,zz
+       DO j=1,yy
+        DO i=1,xx
+         ux(1,i)=r1_(1,i,j,k)                            ! 1:rho
+         ux(2,i)=r1_(2,i,j,k)/r1_(1,i,j,k)               ! 2:YF
+         ux(3,i)=ux(1,i)*uu_(i,j,k)                      ! 3:rho*u
+         ux(4,i)=ux(3,i)*ux(2,i)                         ! 4:rho*u*YF
+         ux(5,i)=r1_(3,i,j,k)/r1_(1,i,j,k)               ! 5:YO
+         ux(6,i)=ux(3,i)*ux(5,i)                         ! 6:rho*u*YO
+
+         uux (1,i)=ux(2,i)
+         uux (2,i)=ux(5,i)
+        ENDDO
+
+         CALL dfnonp (xx,hx, ux(1:nd1,:), dux(1:nd1,:),nd1,1)
+         CALL d2fnonp(xx,hx,uux(1:nd2,:),d2ux(1:nd2,:),nd2,1)
+
+        DO i=1,xx
+
+! -( d(rho*u)/dx )
+         f_(1,i,j,k) = f_(1,i,j,k) - dux(3,i)   ! continuity
+
+! -( d(rho*u*Yr)/dx ) + d(rho*D* d(Yr)/dx)/dx
+! = -( d(rho*u*Yr)/dx )
+!    + D* (rho* d2(Yr)/dx2 + d(rho)/dx * d(Yr)/dx )
+         f_(2,i,j,k) = f_(2,i,j,k) - dux(4,i) + diff*(ux(1,i)*d2ux(1,i)+dux(1,i)*dux(2,i)) ! species conservation
+         f_(3,i,j,k) = f_(3,i,j,k) - dux(6,i) + diff*(ux(1,i)*d2ux(2,i)+dux(1,i)*dux(5,i)) ! species conservation
+
+        ENDDO
+       ENDDO
+      ENDDO
+
+!! y-direction
+      DO k=1,zz
+       DO i=1,xx
+        DO j=1,yy
+         uy(1,j)=r1_(1,i,j,k)                            ! 1:rho
+         uy(2,j)=r1_(2,i,j,k)/r1_(1,i,j,k)               ! 2:YF
+         uy(3,j)=uy(1,j)*vv_(i,j,k)                      ! 3:rho*v
+         uy(4,j)=uy(3,j)*uy(2,j)                         ! 4:rho*v*YF
+         uy(5,j)=r1_(3,i,j,k)/r1_(1,i,j,k)               ! 5:YO
+         uy(6,j)=uy(3,j)*uy(5,j)                         ! 6:rho*v*YO
+
+         uuy (1,j)=uy(2,j)
+         uuy (2,j)=uy(5,j)
+        ENDDO
+
+        CALL dfp (yy,hy, uy(1:nd1,:), duy(1:nd1,:),nd1,2)
+        CALL d2fp(yy,hy,uuy(1:nd2,:),d2uy(1:nd2,:),nd2,2)
+
+        DO j=1,yy
+! -( d(rho*v)/dy ) 
+         f_(1,i,j,k)=f_(1,i,j,k)-duy(3,j)   ! continuity
+
+! -( d(rho*v*Yr)/dy ) + d(rho*D* d(Yr)/dy)/dy
+! = -( d(rho*v*Yr)/dy ) 
+!    + D* (rho* d2(Yr)/dyy2 + d(rho)/dy * d(Yr)/dy )
+         f_(2,i,j,k)=f_(2,i,j,k)-duy(4,j) + diff*(uy(1,j)*d2uy(1,j)+duy(1,j)*duy(2,j))  ! species conserv.
+         f_(3,i,j,k)=f_(3,i,j,k)-duy(6,j) + diff*(uy(1,j)*d2uy(2,j)+duy(1,j)*duy(5,j))  ! species conserv.
+        ENDDO
+       ENDDO
+      ENDDO
+
+
 !    Boundary condition
-!$omp parallel do private(k,j,i) NUM_THREADS(6) schedule(static)
       DO k=1,zz
        DO j=1,yy
         DO i=1,yrsw
          f_(2,i,j,k)=r1_(1,i,j,k)*0.+f_(1,i,j,k)*in_yr
+         f_(3,i,j,k)=r1_(1,i,j,k)*0.+f_(1,i,j,k)*in_yr
         ENDDO
         IF (uu_(xx,j,k).lt.0.) f_(2,xx,j,k)=f_(1,xx,j,k)*r1_(2,xx,j,k)/r1_(1,xx,j,k)
+        IF (uu_(xx,j,k).lt.0.) f_(3,xx,j,k)=f_(1,xx,j,k)*r1_(3,xx,j,k)/r1_(1,xx,j,k)
        ENDDO
       ENDDO
 
@@ -419,7 +682,7 @@ module m_fdm_calc
 
       integer :: istage,xx,yy,zz,i
       real*8 :: uu_(xx,yy,zz),vv_(xx,yy,zz),ww_(xx,yy,zz)
-      real*8 :: yy1(2,xx,yy,zz),yy2(2,xx,yy,zz),yyf(2,xx,yy,zz)
+      real*8 :: yy1(neq,xx,yy,zz),yy2(neq,xx,yy,zz),yyf(neq,xx,yy,zz)
 
 
        istage=1; CALL substep(yy1,yy1,yy2,yyf,xx,yy,zz,istage,uu_,vv_,ww_)
@@ -437,7 +700,7 @@ module m_fdm_calc
 
        integer :: i,j,k,xx,yy,zz
        real*8 :: uu_(xx,yy,zz),vv_(xx,yy,zz),ww_(xx,yy,zz)
-       real*8 :: yy1(2,xx,yy,zz),yy2(2,xx,yy,zz),yyf(2,xx,yy,zz)
+       real*8 :: yy1(neq,xx,yy,zz),yy2(neq,xx,yy,zz),yyf(neq,xx,yy,zz)
 
 
        CALL RK4(xx,yy,zz,uu_,vv_,ww_,yy1,yy2,yyf)
@@ -451,7 +714,7 @@ module m_fdm_calc
 
       integer :: i,j,k,xx,yy,zz,istage
       real*8 :: at,bt , wrate , yr
-      real*8 :: ri(2,xx,yy,zz),r1(2,xx,yy,zz),r2(2,xx,yy,zz),f(2,xx,yy,zz)
+      real*8 :: ri(neq,xx,yy,zz),r1(neq,xx,yy,zz),r2(neq,xx,yy,zz),f(neq,xx,yy,zz)
       real*8 :: a(5),b(5)
       real*8 :: uu_(xx,yy,zz),vv_(xx,yy,zz),ww_(xx,yy,zz)
       integer :: nfinal, iscr, mspec, mpict, msave, nmindt, nv
@@ -473,11 +736,10 @@ module m_fdm_calc
       IF(istage<5) THEN
         at=a(istage)*fdmdt
         bt=(b(istage)-a(istage))*fdmdt
-!$omp parallel do private(k,j,i,nv) NUM_THREADS(6) schedule(static)
         DO k=1,zz
          DO j=1,yy
           DO i=1,xx
-           DO nv=1,2
+           DO nv=1,neq
             r1(nv,i,j,k)=r1(nv,i,j,k)+at*f(nv,i,j,k)
             r2(nv,i,j,k)=r1(nv,i,j,k)+bt*f(nv,i,j,k)
            ENDDO
@@ -486,21 +748,20 @@ module m_fdm_calc
         ENDDO
       ELSE
         bt=b(istage)*fdmdt
-        sumc=0.
-        sum_wrate=0.
-!$omp parallel do private(k,j,i,nv) NUM_THREADS(6) schedule(static)
         DO k=1,zz
          DO j=1,yy
           DO i=1,xx
-           DO nv=1,2
+           DO nv=1,neq
             r1(nv,i,j,k)=r1(nv,i,j,k)+bt*f(nv,i,j,k)
            ENDDO
 
 !==========rho=1 treatment
            r1(2,i,j,k)=r1(2,i,j,k)/r1(1,i,j,k)
+           r1(3,i,j,k)=r1(3,i,j,k)/r1(1,i,j,k)
            r1(1,i,j,k)=1.
 !==========Max Yr=1 treatment
            r1(2,i,j,k)=MIN(in_yr,r1(2,i,j,k))
+           r1(3,i,j,k)=MIN(in_yr,r1(3,i,j,k))
 !==========Min Yr=0 treatment
 !          r1(2,i,j,k)=MAX(out_yr,r1(2,i,j,k))
 

m_mpi.f90

@@ -64,7 +64,8 @@ contains
 
 
     ! initializing MPI environment
-    call MPI_INIT_THREAD(MPI_THREAD_SERIALIZED, mpi_provide, mpi_err)
+    !call MPI_INIT_THREAD(MPI_THREAD_SERIALIZED, mpi_provide, mpi_err)
+    call MPI_INIT(mpi_err)
     call MPI_Comm_size(MPI_COMM_WORLD,numprocs_world,mpi_err)
     call MPI_Comm_rank(MPI_COMM_WORLD,myid_world,mpi_err)
 

m_openmpi.f90

@@ -64,7 +64,8 @@ contains
 
 
     ! initializing MPI environment
-    call MPI_INIT_THREAD(MPI_THREAD_SERIALIZED, mpi_provide, mpi_err)
+    !call MPI_INIT_THREAD(MPI_THREAD_SERIALIZED, mpi_provide, mpi_err)
+    call MPI_INIT(mpi_err)
     call MPI_Comm_size(MPI_COMM_WORLD,numprocs_world,mpi_err)
     call MPI_Comm_rank(MPI_COMM_WORLD,myid_world,mpi_err)
 

m_timing.f90

@@ -19,8 +19,8 @@ contains
     ! make the default job runlimit to be 6 months
     job_runlimit = 6 * 30 * 24 * 60
 
-    ! make the default job runlimit to be 12 hours
-    job_runlimit = 12 * 60
+    ! make the default job runlimit to be 12 hours -> one month
+    job_runlimit = 12 * 60 * 30 * 2
 
     write(out,'(a30,3x,i20)') 'job_runlimit (default):',job_runlimit
 

main.f90

@@ -67,10 +67,7 @@ program x_code
   ! Initializing the LES stuff
   if (les) call m_les_begin
 
-     if (task.eq.'hydro' .and. myid.eq.master) then
      call prepare_fdm
-     end if
-     call MPI_BARRIER(MPI_COMM_TASK, mpi_err)
 
      call stat_main
 
@@ -134,20 +131,7 @@ program x_code
 
         if (fdm_sw.ne.0) then  ! Reacting
             call gather_4
-            if (myid.eq.master) then
                call fdm_exe
-               do id_to=1,numprocs-1
-                 call MPI_Isend (dummy, 1, MPI_REAL, id_to, 0, MPI_COMM_TASK, mpi_request, mpi_err)
-               end do
-            else
-               call MPI_Irecv (dummy, 1, MPI_REAL, master, 0, MPI_COMM_TASK, mpi_request, mpi_err)
-               do
-                 call MPI_Test(mpi_request, flag, mpi_status, mpi_err)
-                 if (flag) exit
-                 call SLEEPQQ (200)
-               end do
-            end if
-            call MPI_BARRIER(MPI_COMM_TASK, mpi_err)
         endif
 
         ! RHS for scalars
@@ -202,18 +186,18 @@ program x_code
         ! write the restart file if it's the time
         if (mod(itime,IPRINT2).eq.0) then
            call restart_write_parallel
-           !remove hit restart file
-           call system ('rm -f '//run_name//'.64.'//ex_file_ext)
            if (fdm_sw.ne.0) then  ! Reacting
-               call gather_4
+               !call gather_4
+               if (TIME.ge.bgn_save) then
+                 call fdm_restart_write
                  if (myid.eq.master) then
-                  !write itime, savenum pair
-                  if(TIME.ge.bgn_save) call fdm_restart_write
-                  write (restart_pair, *) file_ext, fullsavenum-1
-                  call write_tp_field
+                   write (restart_pair, *) file_ext, fullsavenum-1, '\n'
+                   !call write_tp_field
                  end if
-               call MPI_BARRIER(MPI_COMM_TASK, mpi_err)
                end if
+           endif
+           !remove HIT3D restart file except the latest
+           call system ('rm -f '//run_name//'.64.'//ex_file_ext)
            ex_file_ext = file_ext
         endif
 
@@ -266,7 +250,7 @@ program x_code
            if (task_split) call fields_to_stats
            if (mod(itime,iprint1).eq.0) call stat_main
            if (fdm_sw.eq.0) then  ! non-reacting
-              if (mod(itime,iprint2).eq.0) call io_write_4
+              !if (mod(itime,iprint2).eq.0) call io_write_4
            endif
         end if
      end if stats

write_stats.f90

@@ -31,6 +31,7 @@ subroutine write_tp_field
      do i=1,nx,spx
         write(tp_field,'(6e18.9)') dx*real(i),dy*real(j),dz*real(k),u_(i,j,k),v_(i,j,k),w_(i,j,k)   ! u',v',w'
      enddo
+     call flush(tp_field)
     enddo
    enddo
 
@@ -49,13 +50,12 @@ subroutine write_tp_field
                                  (1.-y1(2,i,j,ii)/y1(1,i,j,ii)),u_(i,j,ii), &
                                   v_(i,j,ii),w_(i,j,ii)
      enddo
+     call flush(tp_field)
     enddo
    enddo
 
  endif
 
- call flush(tp_field)
-
 897      format('ZONE T="time= ',f10.5,'" I=  ',i4,' J=  ',i4,' K=  ',i4,' F=POINT')
 
 end subroutine write_tp_field