write_spx1.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: WRITE_SPX1                                             C
!  Purpose: write out the time-dependent restart records (REAL)        C
!                                                                      C
!  Author: P. Nicoletti                               Date: 13-DEC-91  C
!  Reviewer: P. Nicoletti, W. Rogers, M. Syamlal      Date: 24-JAN-92  C
!                                                                      C
!  Revision Number:                                                    C
!  Purpose:                                                            C
!  Author:                                            Date: dd-mmm-yy  C
!  Reviewer:                                          Date: dd-mmm-yy  C
!                                                                      C
!  Literature/Document References:                                     C
!                                                                      C
!  Variables referenced: TIME, NSTEP, EP_g, P_g, P_star, U_g           C
!                        V_g, W_g, U_s, V_s, W_s, ROP_s, T_g, T_s      C
!                        IJKMAX2, MMAX                                 C
!  Variables modified: None                                            C
!                                                                      C
!  Local variables:  LC, N, NEXT_REC, NUM_REC                          C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
!
      SUBROUTINE write_spx1(L, unit_add)
!...Translated by Pacific-Sierra Research VAST-90 2.06G5  12:17:31  12/09/98
!...Switches: -xf
!
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE cdist
      USE compar
      USE cutcell
      USE exit, only: mfix_exit
      USE fldvar
      USE funits
      USE geometry
      USE machine
      USE mpi_utility
      USE output
      USE param
      USE param1
      USE physprop
      USE run
      USE rxns
      USE scalars
      use discretelement, only: des_continuum_coupled
      use discretelement, only: discrete_element
      use discretelement, only: print_des_data

!//       USE tmp_array
      IMPLICIT NONE
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
!
!             flag whether to write a particular SPx file
      INTEGER L

!              offset for use in post_mfix
      INTEGER  unit_add
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
!
! local variables
!
!//
      double precision, allocatable :: array1(:)     !//
      double precision, allocatable :: array2(:)     !//

!             loop counters
      INTEGER LC, NN
!
!             Pointer to the next record
      INTEGER NEXT_REC
!
!              Number of records written each time step
      INTEGER  NUM_REC

      INTEGER  uspx   ! UNIT_SPX + offset from post_mfix
      CHARACTER(LEN=50), DIMENSION(1) :: LINE   !error message
      double precision, dimension(:), allocatable :: TMP_VAR

      allocate(tmp_var(dimension_3))

!-----------------------------------------------
      uspx = unit_spx + unit_add

!
      if (mype .eq.pe_io) then
         allocate (array1(ijkmax2))   !//
         allocate (array2(ijkmax3))   !//
      else
         allocate (array1(1))   !//
         allocate (array2(1))   !//
      end if
!
! ".SP1" FILE         EP_g    [ ROP_g, RO_g  must be calculated ...
!                                        not written out ]
!
      SELECT CASE (l)
      CASE (1)
         if (mype.eq.pe_io.or.bdist_io) then
            READ (uspx + l, rec=3) next_rec, num_rec
            num_rec = next_rec
            WRITE (uspx + l, rec=next_rec) REAL(TIME), NSTEP
            next_rec = next_rec + 1
         end if

! Explicitly coupled simulations do not need to rebin particles to
! the fluid grid every time step. However, this implies that the
! fluid cell information and interpolation weights become stale.
         IF(discrete_element .AND. .NOT.des_continuum_coupled) THEN
! Bin particles to fluid grid.
            CALL particles_in_cell
! Calculate interpolation weights
            CALL calc_interp_weights
! Calculate mean fields (EPg).
            CALL comp_mean_fields
         ENDIF

         if (bdist_io) then
            IF(re_indexing) THEN
               CALL unshift_dp_array(ep_g,tmp_var)
               call out_bin_r(uspx+l,tmp_var,size(tmp_var),next_rec)
            ELSE
               call out_bin_r(uspx+l,ep_g,size(ep_g),next_rec)
            ENDIF
!           call OUT_BIN_R(uspx+L,EP_g,size(EP_g),NEXT_REC)
         else
            call gatherwritespx (ep_g,array2, array1, uspx+l, next_rec)   !//
         end if
         if (mype .eq. pe_io.or.bdist_io) then
            num_rec = next_rec - num_rec
            WRITE (uspx + l, rec=3) next_rec, num_rec
            if(unit_add == 0) CALL flush_bin (uspx + l)
         end if

! The call made in make_arrays captures the initial state of the system
! as the input and RES files for DES runs are read afte the the first
! call to this routine.
         IF(discrete_element.AND.print_des_data) THEN
            IF(time /= zero .OR. trim(run_type)=='RESTART_1') &
               CALL write_des_data
         ENDIF

!        call MPI_Barrier(MPI_COMM_WORLD,mpierr)  !//PAR_I/O enforce barrier here
!
! ".SP2" FILE         P_g , P_star
!
      CASE (2)
         if (mype.eq.pe_io.or.bdist_io) then
            READ (uspx + l, rec=3) next_rec, num_rec
            num_rec = next_rec
            WRITE (uspx + l, rec=next_rec) REAL(TIME), NSTEP
            next_rec = next_rec + 1
         end if
         if (bdist_io) then
            IF(re_indexing) THEN
               CALL unshift_dp_array(p_g,tmp_var)
               call out_bin_r(uspx+l,tmp_var,size(tmp_var),next_rec)
               CALL unshift_dp_array(p_star,tmp_var)
               call out_bin_r(uspx+l,tmp_var,size(tmp_var),next_rec)
            ELSE
              call out_bin_r(uspx+l,p_g,size(p_g),next_rec)
              call out_bin_r(uspx+l,p_star,size(p_star),next_rec)
            ENDIF
!           call OUT_BIN_R(uspx+L,P_g,size(P_g),NEXT_REC)
!           call OUT_BIN_R(uspx+L,P_star,size(P_star),NEXT_REC)
         else
           call gatherwritespx (p_g,array2, array1, uspx+l, next_rec)   !//
           call gatherwritespx (p_star,array2, array1, uspx+l, next_rec)   !//
         end if
         if (mype.eq.pe_io.or.bdist_io) then
            num_rec = next_rec - num_rec
            WRITE (uspx + l, rec=3) next_rec, num_rec
            if(unit_add == 0) CALL flush_bin (uspx + l)
         end if
!        call MPI_Barrier(MPI_COMM_WORLD,mpierr)  !//PAR_I/O enforce barrier here
!
! ".SP3" FILE         U_g , V_g , W_g
!
      CASE (3)
         if (mype.eq.pe_io.or.bdist_io) then
            READ (uspx + l, rec=3) next_rec, num_rec
            num_rec = next_rec
            WRITE (uspx + l, rec=next_rec) REAL(TIME), NSTEP
            next_rec = next_rec + 1
         end if
         if (bdist_io) then
            IF(re_indexing) THEN
               CALL unshift_dp_array(u_g,tmp_var)
               call out_bin_r(uspx+l,tmp_var,size(tmp_var),next_rec)
               CALL unshift_dp_array(v_g,tmp_var)
               call out_bin_r(uspx+l,tmp_var,size(tmp_var),next_rec)
               CALL unshift_dp_array(w_g,tmp_var)
               call out_bin_r(uspx+l,tmp_var,size(tmp_var),next_rec)
            ELSE
               call out_bin_r(uspx+l,u_g,size(u_g),next_rec)
               call out_bin_r(uspx+l,v_g,size(v_g),next_rec)
               call out_bin_r(uspx+l,w_g,size(w_g),next_rec)
            ENDIF
!           call OUT_BIN_R(uspx+L,U_g,size(U_g),NEXT_REC)
!           call OUT_BIN_R(uspx+L,V_g,size(V_g),NEXT_REC)
!           call OUT_BIN_R(uspx+L,W_g,size(W_g),NEXT_REC)
         else
           call gatherwritespx (u_g,array2, array1, uspx+l, next_rec)   !//
           call gatherwritespx (v_g,array2, array1, uspx+l, next_rec)   !//
           call gatherwritespx (w_g,array2, array1, uspx+l, next_rec)   !//
         end if
         if (mype.eq.pe_io.or.bdist_io) then
            num_rec = next_rec - num_rec
            WRITE (uspx + l, rec=3) next_rec, num_rec
            if(unit_add == 0) CALL flush_bin (uspx + l)
         end if
!        call MPI_Barrier(MPI_COMM_WORLD,mpierr)  !//PAR_I/O enforce barrier here
!
! ".SP4" FILE         U_s , V_s , W_s
!
      CASE (4)
         if (mype.eq.pe_io.or.bdist_io) then
            READ (uspx + l, rec=3) next_rec, num_rec
            num_rec = next_rec
            WRITE (uspx + l, rec=next_rec) REAL(TIME), NSTEP
            next_rec = next_rec + 1
         end if
         if (bdist_io) then
            DO lc = 1, mmax
               IF(re_indexing) THEN
                  CALL unshift_dp_array(u_s(:,lc),tmp_var)
                  call out_bin_r(uspx+l,tmp_var,size(tmp_var),next_rec)
                  CALL unshift_dp_array(v_s(:,lc),tmp_var)
                  call out_bin_r(uspx+l,tmp_var,size(tmp_var),next_rec)
                  CALL unshift_dp_array(w_s(:,lc),tmp_var)
                  call out_bin_r(uspx+l,tmp_var,size(tmp_var),next_rec)
               ELSE
                  call out_bin_r(uspx+l,u_s(:,lc),Size(u_s(:,lc)),next_rec)
                  call out_bin_r(uspx+l,v_s(:,lc),Size(v_s(:,lc)),next_rec)
                  call out_bin_r(uspx+l,w_s(:,lc),Size(w_s(:,lc)),next_rec)
               ENDIF
            ENDDO
!        DO LC = 1, MMAX
!          call OUT_BIN_R(uspx+L,U_s(:,LC),Size(U_s(:,LC)),NEXT_REC)
!          call OUT_BIN_R(uspx+L,V_s(:,LC),Size(V_s(:,LC)),NEXT_REC)
!          call OUT_BIN_R(uspx+L,W_s(:,LC),Size(W_s(:,LC)),NEXT_REC)
!        END DO
         else
         DO lc = 1, mmax
            call gatherwritespx (u_s(:,lc),array2, array1, uspx+l, next_rec)
            call gatherwritespx (v_s(:,lc),array2, array1, uspx+l, next_rec)
            call gatherwritespx (w_s(:,lc),array2, array1, uspx+l, next_rec)
         END DO
         end if
         if (mype.eq.pe_io.or.bdist_io) then
            num_rec = next_rec - num_rec
            WRITE (uspx + l, rec=3) next_rec, num_rec
            if(unit_add == 0) CALL flush_bin (uspx + l)
         end if
!        call MPI_Barrier(MPI_COMM_WORLD,mpierr)  !//PAR_I/O enforce barrier here
!
! ".SP5" FILE         ROP_s
!
      CASE (5)
         if (mype.eq.pe_io.or.bdist_io) then
            READ (uspx + l, rec=3) next_rec, num_rec
            num_rec = next_rec
            WRITE (uspx + l, rec=next_rec) REAL(TIME), NSTEP
            next_rec = next_rec + 1
         end if
         if (bdist_io) then
            DO lc = 1, mmax
               IF(re_indexing) THEN
                  CALL unshift_dp_array(rop_s(:,lc),tmp_var)
                  call out_bin_r(uspx+l,tmp_var,size(tmp_var),next_rec)
               ELSE
                  call out_bin_r(uspx+l,rop_s(:,lc),size(rop_s(:,lc)), next_rec)
               ENDIF
            ENDDO
!          DO LC = 1, MMAX
!            call OUT_BIN_R(uspx+L,ROP_s(:,LC),size(ROP_s(:,LC)), NEXT_REC)
!         END DO
         else
         DO lc = 1, mmax
            call gatherwritespx (rop_s(:,lc),array2, array1, uspx+l, next_rec)
         END DO
         end if
         if (mype.eq.pe_io.or.bdist_io) then
            num_rec = next_rec - num_rec
            WRITE (uspx + l, rec=3) next_rec, num_rec
            if(unit_add == 0) CALL flush_bin (uspx + l)
         end if
!        call MPI_Barrier(MPI_COMM_WORLD,mpierr)  !//PAR_I/O enforce barrier here
!
! ".SP6" FILE         T_g  , T_s
!
      CASE (6)
         if (mype.eq.pe_io.or.bdist_io) then
            READ (uspx + l, rec=3) next_rec, num_rec
            num_rec = next_rec
            WRITE (uspx + l, rec=next_rec) REAL(TIME), NSTEP
            next_rec = next_rec + 1
         end if
         if (bdist_io) then
            IF(re_indexing) THEN
               CALL unshift_dp_array(t_g,tmp_var)
               call out_bin_r(uspx+l,tmp_var,size(tmp_var),next_rec)
               DO lc = 1, mmax
                  CALL unshift_dp_array(t_s(:,lc),tmp_var)
                  call out_bin_r(uspx+l,tmp_var,size(tmp_var), next_rec)
               END DO
            ELSE
               call out_bin_r(uspx+l,t_g,size(t_g), next_rec)
               DO lc = 1, mmax
                  call out_bin_r(uspx+l,t_s(:,lc),size(t_s(:,lc)), next_rec)
               END DO
            ENDIF
!           call OUT_BIN_R(uspx+L,T_g,size(T_g), NEXT_REC)
!          DO LC = 1, MMAX
!            call OUT_BIN_R(uspx+L,T_s(:,LC),size(T_s(:,LC)), NEXT_REC)
!          END DO
         else
         call gatherwritespx (t_g,array2, array1, uspx+l, next_rec)   !//
          DO lc = 1, mmax
            call gatherwritespx (t_s(:,lc),array2, array1, uspx+l, next_rec)
          END DO
         end if
         if (mype.eq.pe_io.or.bdist_io) then
            num_rec = next_rec - num_rec
            WRITE (uspx + l, rec=3) next_rec, num_rec
            if(unit_add == 0) CALL flush_bin (uspx + l)
         end if
!
! ".SP7" FILE         X_g, X_s
!
      CASE (7)
         if (mype.eq.pe_io.or.bdist_io) then
            READ (uspx + l, rec=3) next_rec, num_rec
            num_rec = next_rec
            WRITE (uspx + l, rec=next_rec) REAL(TIME), NSTEP
            next_rec = next_rec + 1
         end if
         if (bdist_io) then
            IF(re_indexing) THEN
               DO nn = 1, nmax(0)
                  CALL unshift_dp_array(x_g(:,nn),tmp_var)
                  call out_bin_r(uspx+l,tmp_var,size(tmp_var),next_rec)
               END DO
               DO lc = 1, mmax
                  DO nn = 1, nmax(lc)
                     CALL unshift_dp_array(x_s(:,lc,nn),tmp_var)
                     call out_bin_r(uspx+l,tmp_var,size(tmp_var), next_rec)
                  ENDDO
               END DO
            ELSE
               DO nn = 1, nmax(0)
                  call out_bin_r(uspx+l,x_g(:,nn),size(x_g(:,nn)), next_rec)
               END DO
               DO lc = 1, mmax
                  DO nn = 1, nmax(lc)
                     call out_bin_r(uspx+l,x_s(:,lc,nn),size(x_s(:,lc,nn)), next_rec)
                  END DO
               END DO
            ENDIF

!           DO NN = 1, NMAX(0)
!             call OUT_BIN_R(uspx+L,X_G(:,nn),size(X_G(:,nn)), NEXT_REC)
!           END DO
!           DO LC = 1, MMAX
!            DO NN = 1, NMAX(LC)
!               call OUT_BIN_R(uspx+L,X_s(:,LC,nn),size(X_s(:,LC,nn)), NEXT_REC)
!            END DO
!           END DO

         else
          DO nn = 1, nmax(0)
            call gatherwritespx (x_g(:,nn),array2, array1, uspx+l, next_rec)
          END DO
          DO lc = 1, mmax
            DO nn = 1, nmax(lc)
               call gatherwritespx (x_s(:,lc,nn),array2, array1, uspx+l, next_rec)
            END DO
          END DO
         end if
         if (mype.eq.pe_io.or.bdist_io) then
            num_rec = next_rec - num_rec
            WRITE (uspx + l, rec=3) next_rec, num_rec
            if(unit_add == 0) CALL flush_bin (uspx + l)
         end if
!        call MPI_Barrier(MPI_COMM_WORLD,mpierr)  !//PAR_I/O enforce barrier here
!
! ".SP8" FILE         THETA_m
!
      CASE (8)
         if (mype.eq.pe_io.or.bdist_io) then
            READ (uspx + l, rec=3) next_rec, num_rec
            num_rec = next_rec
            WRITE (uspx + l, rec=next_rec) REAL(TIME), NSTEP
            next_rec = next_rec + 1
         end if
         if (bdist_io) then
            IF(re_indexing) THEN
               DO lc = 1, mmax
                  CALL unshift_dp_array(theta_m(:,lc),tmp_var)
                  call out_bin_r(uspx+l,tmp_var,size(tmp_var),next_rec)
               ENDDO
            ELSE
               DO lc = 1, mmax
                  call out_bin_r(uspx+l,theta_m(:,lc),size(theta_m(:,lc)), next_rec)
               END DO
            ENDIF
!          DO LC = 1, MMAX
!            call OUT_BIN_R(uspx+L,THETA_m(:,LC),size(THETA_m(:,LC)), NEXT_REC)
!         END DO
         else
          DO lc = 1, mmax
            call gatherwritespx (theta_m(:,lc),array2, array1, uspx+l, next_rec)
          END DO
         end if
         if (mype.eq.pe_io .or. bdist_io) then
            num_rec = next_rec - num_rec
            WRITE (uspx + l, rec=3) next_rec, num_rec
            if(unit_add == 0) CALL flush_bin (uspx + l)
         end if
!        call MPI_Barrier(MPI_COMM_WORLD,mpierr)  !//PAR_I/O enforce barrier here
!
! ".SP9" FILE         Scalar
!
      CASE (9)
         if (mype.eq.pe_io.or.bdist_io) then
            READ (uspx + l, rec=3) next_rec, num_rec
            num_rec = next_rec
            WRITE (uspx + l, rec=next_rec) REAL(TIME), NSTEP
            next_rec = next_rec + 1
         end if
         if (bdist_io) then
            IF(re_indexing) THEN
               DO lc = 1, nscalar
                  CALL unshift_dp_array(scalar(:,lc),tmp_var)
                  call out_bin_r(uspx+l,tmp_var,size(tmp_var),next_rec)
               ENDDO
            ELSE
               DO lc = 1, nscalar
                 call out_bin_r(uspx+l,scalar(:,lc),size(scalar(:,lc)), next_rec)
               END DO
            ENDIF
!          DO LC = 1, Nscalar
!            call OUT_BIN_R(uspx+L,Scalar(:,LC),size(Scalar(:,LC)), NEXT_REC)
!          END DO
         else
          DO lc = 1, nscalar
            call gatherwritespx (scalar(:,lc),array2, array1, uspx+l, next_rec)
          END DO
         end if
         if (mype.eq.pe_io.or.bdist_io) then
            num_rec = next_rec - num_rec
            WRITE (uspx + l, rec=3) next_rec, num_rec
            if(unit_add == 0) CALL flush_bin (uspx + l)
         end if
!        call MPI_Barrier(MPI_COMM_WORLD,mpierr)  !//PAR_I/O enforce barrier here
!
      CASE (10)  ! Reaction rates

         if (mype.eq.pe_io.or.bdist_io) then
            READ (uspx + l, rec=3) next_rec, num_rec
            num_rec = next_rec
            WRITE (uspx + l, rec=next_rec) REAL(TIME), NSTEP
            next_rec = next_rec + 1
         end if
         if (bdist_io) then
            IF(re_indexing) THEN
               DO lc = 1, nrr
                  CALL unshift_dp_array(reactionrates(:,lc),tmp_var)
                  call out_bin_r(uspx+l,tmp_var,size(tmp_var),next_rec)
               ENDDO
            ELSE
               DO lc = 1, nrr
                  call out_bin_r(uspx+l,reactionrates(:,lc),size(reactionrates(:,lc)), next_rec)
               END DO
            ENDIF
!          DO LC = 1, nRR
!            call OUT_BIN_R(uspx+L,ReactionRates(:,LC),size(ReactionRates(:,LC)), NEXT_REC)
!          END DO
         else
          DO lc = 1, nrr
            call gatherwritespx (reactionrates(:,lc),array2, array1, uspx+l, next_rec)
          END DO
         end if
         if (mype.eq.pe_io.or.bdist_io) then
            num_rec = next_rec - num_rec
            WRITE (uspx + l, rec=3) next_rec, num_rec
            if(unit_add == 0) CALL flush_bin (uspx + l)
         end if
!
! ".SP11" FILE         turbulence
!
      CASE (11)
         if (mype.eq.pe_io.or.bdist_io) then
            READ (uspx + l, rec=3) next_rec, num_rec
            num_rec = next_rec
            WRITE (uspx + l, rec=next_rec) REAL(TIME), NSTEP
            next_rec = next_rec + 1
         end if
         if (k_epsilon) then
            if (bdist_io) then
            IF(re_indexing) THEN
               CALL unshift_dp_array(k_turb_g,tmp_var)
               call out_bin_r(uspx+l,tmp_var,size(tmp_var),next_rec)
               CALL unshift_dp_array(e_turb_g,tmp_var)
               call out_bin_r(uspx+l,tmp_var,size(tmp_var),next_rec)
            ELSE
               call out_bin_r(uspx+l,k_turb_g,size(k_turb_g), next_rec)
               call out_bin_r(uspx+l,e_turb_g,size(e_turb_g), next_rec)
            ENDIF
!             call OUT_BIN_R(uspx+L,K_Turb_G,size(K_Turb_G), NEXT_REC)
!             call OUT_BIN_R(uspx+L,E_Turb_G,size(E_Turb_G), NEXT_REC)

            else
             call gatherwritespx (k_turb_g,array2, array1, uspx+l, next_rec)
             call gatherwritespx (e_turb_g,array2, array1, uspx+l, next_rec)
            end if
          end if

           if (mype.eq.pe_io.or.bdist_io) then
            num_rec = next_rec - num_rec
            WRITE (uspx + l, rec=3) next_rec, num_rec
            if(unit_add == 0) CALL flush_bin (uspx + l)
           end if
!        call MPI_Barrier(MPI_COMM_WORLD,mpierr)  !//PAR_I/O enforce barrier here
!
!
      CASE DEFAULT
            line(1) = 'Unknown SPx file index'
            CALL write_error ('WRITE_SPX1', line, 1)
            CALL mfix_exit(mype)
      END SELECT

!//      call unlock_tmp_array
!
      deallocate (array1)
      deallocate (array2)
      deallocate (tmp_var)
!
      RETURN
      END SUBROUTINE write_spx1

      subroutine gatherwritespx(VAR, array2, array1, uspxL, NEXT_REC)
        USE geometry
        USE compar           !//
        USE mpi_utility      !//d pnicol : for gatherWriteSpx
        USE sendrecv         !//d pnicol : for gatherWriteSpx
        USE cutcell
        USE in_binary_512
        USE param, only: dimension_3
        USE param1, only: undefined
        IMPLICIT NONE
        integer uspxL, NEXT_REC
        double precision, dimension(ijkmax2) :: array1
        double precision, dimension(ijkmax3) :: array2
        double precision, dimension(DIMENSION_3) :: VAR
        double precision, dimension(:), allocatable :: TMP_VAR

        allocate(tmp_var(dimension_3))

!       call MPI_Barrier(MPI_COMM_WORLD,mpierr)  !//PAR_I/O enforce barrier here
        IF(re_indexing) THEN
           tmp_var = undefined
           CALL unshift_dp_array(var,tmp_var)
           CALL gather (tmp_var,array2,root)
        ELSE
           CALL gather (var,array2,root)
        ENDIF
!        call gather (VAR,array2,root)  !//d pnicol

!       call MPI_Barrier(MPI_COMM_WORLD,mpierr)  !//PAR_I/O enforce barrier here
        if (mype.eq.pe_io) then
           call convert_to_io_dp(array2,array1,ijkmax2)
           CALL out_bin_r (uspxl, array1, ijkmax2, next_rec)
        end if

        deallocate(tmp_var)

      End subroutine gatherwritespx



      subroutine gatherwritespx_netcdf(VAR, arr1, arr2 , arr4d, ncid, varid , &
                nx,ny,nz,ijkmax2_use , ijkmax3_use)


     USE geometry
     use param, only: dimension_3
     USE compar           !//
     USE mpi_utility      !//d pnicol : for gatherWriteSpx
     USE sendrecv         !//d pnicol : for gatherWriteSpx
     USE mfix_netcdf
     USE in_binary_512

     IMPLICIT NONE

     integer          :: ncid , varid , nx,ny,nz , ijkmax2_use , ijkmax3_use
     integer          :: ii , jj , kk , ijk

     double precision ::  arr1(ijkmax2_use)
     double precision ::  arr2(ijkmax3_use)
     double precision ::  arr4d(nx,ny,nz,1)
     double precision ::  var(dimension_3)

     call gather(var,arr2,root)
     if (mype .eq. pe_io) then
        call convert_to_io_dp(arr2,arr1,ijkmax2_use)

        ijk = 0
        do kk = 1,nz
           do jj = 1,ny
              do ii = 1,nx
                 ijk = ijk + 1
                 arr4d(ii,jj,kk,1) = arr1(ijk)
              end do
           end do
        end do

        call mfix_check_netcdf( mfix_nf90_put_var(ncid,varid,arr4d) )

     end if


    End subroutine gatherwritespx_netcdf



      subroutine gatherwritespx_netcdf_int(VAR, arr1, arr2 , arr4d, ncid, &
                varid , nx,ny,nz,ijkmax2_use , ijkmax3_use)


     USE geometry
     use param, only: dimension_3
     USE compar           !//
     USE mpi_utility      !//d pnicol : for gatherWriteSpx
     USE sendrecv         !//d pnicol : for gatherWriteSpx
     USE mfix_netcdf
     USE in_binary_512i

     IMPLICIT NONE

     integer          :: ncid , varid , nx,ny,nz , ijkmax2_use , ijkmax3_use
     integer          :: ii , jj , kk , ijk

     integer ::  arr1(ijkmax2_use)
     integer ::  arr2(ijkmax3_use)
     integer :: arr4d(nx,ny,nz,1)
     integer ::   var(dimension_3)

     call gather(var,arr2,root)
     if (mype .eq. pe_io) then
        call convert_to_io_i(arr2,arr1,ijkmax2_use)

        ijk = 0
        do kk = 1,nz
           do jj = 1,ny
              do ii = 1,nx
                 ijk = ijk + 1
                 arr4d(ii,jj,kk,1) = arr1(ijk)
              end do
           end do
        end do

        call mfix_check_netcdf( mfix_nf90_put_var(ncid,varid,arr4d) )

     end if


    End subroutine gatherwritespx_netcdf_int

        subroutine copy_d_to_r(darr,rarr,nx,ny,nz)
        implicit none

        integer          :: nx , ny , nz
        double precision :: darr(*)
        real             :: rarr(nx,ny,*)
        integer          :: i , j , k , ijk


        ijk = 0

        do i = 1,nx
           do j = 1,ny
              do k = 1,nz
                 ijk = ijk + 1
                 rarr(i,j,k) = real(darr(ijk))
              end do
           end do
        end do

        return
        end subroutine copy_d_to_r

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!                          write_mesh_netcdf
!
!
        SUBROUTINE write_mesh_netcdf

        USE param
        USE param1
        USE fldvar
        USE geometry
        USE physprop
        USE run
!       USE funits
        USE scalars
!       USE output
        USE rxns
        USE cdist
        USE compar
        USE mpi_utility
        USE mfix_netcdf
!       USE tmp_array

        implicit none

        integer   :: ncid    , x_dimid , y_dimid , z_dimid  , t_dimid
        integer   :: varid_x , varid_y , varid_z , L , dimids(4)
        integer   :: varid_flag , coords_dimid , varid_coords , coords

        character(LEN=80) :: fname

        double precision, dimension(:) , allocatable :: xloc
        double precision, dimension(:) , allocatable :: yloc
        double precision, dimension(:) , allocatable :: zloc

        integer, dimension(:) , allocatable :: arr1
        integer, dimension(:) , allocatable :: arr2
        integer, dimension(:,:,:,:) , allocatable :: arr4d


        if (.not. mfix_usingnetcdf()) return

        if (.not. bglobalnetcdf) return  ! no netCDF writes asked for

        if (.not. bfirst_netcdf_write) return

        if (mype.eq.pe_io) then
                allocate ( arr1(ijkmax2))
                allocate ( arr2(ijkmax3))
                allocate ( arr4d(imax2,jmax2,kmax2,1))
                allocate ( xloc(imax2) )
                allocate ( yloc(jmax2) )
                allocate ( zloc(kmax2) )
        else
                allocate ( arr1(1))
                allocate ( arr2(1))
                allocate ( arr4d(1,1,1,1))
                allocate ( xloc(1) )
                allocate ( yloc(1) )
                allocate ( zloc(1) )
        end if

        if (mype.eq.pe_io) then
                xloc(1) = -dx(1)
                do l = 2,imax2
                        xloc(l) = xloc(l-1) + dx(l)
                end do

                yloc(1) = -dy(1)
                do l = 2,jmax2
                        yloc(l) = yloc(l-1) + dy(l)
                end do

                zloc(1) = -dz(1)
                do l = 2,kmax2
                        zloc(l) = zloc(l-1) + dz(l)
                end do

                fname = trim(run_name) // "_MESH.nc"
                call mfix_check_netcdf( mfix_nf90_create(fname, nf90_64bit_offset, ncid) )

                call mfix_check_netcdf( mfix_nf90_def_dim(ncid, "x"   , imax2   , x_dimid  ) )
                call mfix_check_netcdf( mfix_nf90_def_dim(ncid, "y"   , jmax2   , y_dimid  ) )
                call mfix_check_netcdf( mfix_nf90_def_dim(ncid, "z"   , kmax2   , z_dimid  ) )
                call mfix_check_netcdf( mfix_nf90_def_dim(ncid, "coordinates"   , 1   , coords_dimid  ) )
                call mfix_check_netcdf( mfix_nf90_def_dim(ncid, "t"   , 1       , t_dimid  ) )  ! 4

                dimids =  (/ x_dimid , y_dimid, z_dimid , t_dimid/)


                ! The dimids array is used to pass the IDs of the dimensions of
                ! the variables. Note that in fortran arrays are stored in
                ! column-major format.
                call mfix_check_netcdf( mfix_nf90_def_var(ncid, "x"  , nf90_double, x_dimid, varid_x  ) )
                call mfix_check_netcdf( mfix_nf90_def_var(ncid, "y"  , nf90_double, y_dimid, varid_y  ) )
                call mfix_check_netcdf( mfix_nf90_def_var(ncid, "z"  , nf90_double, z_dimid, varid_z  ) )
                call mfix_check_netcdf( mfix_nf90_def_var(ncid, "coordinates"  , nf90_int, coords_dimid, varid_coords  ) )
                call mfix_check_netcdf( mfix_nf90_def_var(ncid, "flag"  , nf90_int, dimids, varid_flag  ) )  ! 9


                call mfix_check_netcdf( mfix_nf90_enddef(ncid) )

                coords = 0
                if (coordinates .eq. 'CYLINDRICAL') coords = 1

                call mfix_check_netcdf( mfix_nf90_put_var(ncid,varid_coords,coords) )
                call mfix_check_netcdf( mfix_nf90_put_var(ncid,varid_x,xloc) )
                call mfix_check_netcdf( mfix_nf90_put_var(ncid,varid_y,yloc) )
                call mfix_check_netcdf( mfix_nf90_put_var(ncid,varid_z,zloc) )
        end if


        ! needs to be called by all processes

        call gatherwritespx_netcdf_int(flag, arr1, arr2 , arr4d, ncid, &
               varid_flag , imax2,jmax2,kmax2,ijkmax2 , ijkmax3)



        if (mype.eq.pe_io) then
                call mfix_check_netcdf( mfix_nf90_close(ncid) )
        end if

        deallocate ( arr1  )
        deallocate ( arr2  )
        deallocate ( arr4d )
        deallocate ( xloc  )
        deallocate ( yloc  )
        deallocate ( zloc  )

        return
        end subroutine write_mesh_netcdf

        SUBROUTINE write_netcdf(L, unit_add, the_time)

        USE param
        USE param1
        USE fldvar
        USE geometry
        USE physprop
        USE run
!       USE funits
        USE scalars
!       USE output
        USE rxns
        USE cdist
        USE compar
        USE mpi_utility
        USE mfix_netcdf
!       USE tmp_array


        implicit none

        integer :: L , unit_add , I , nn , ii

        integer   :: ncid , x_dimid , y_dimid , z_dimid
        integer   :: t_dimid
        integer   :: dimids(4) , varid_epg , varid_pg

        integer   :: varid_pstar  , varid_ug , varid_vg , varid_wg
        integer   :: varid_tg , varid_x , varid_y , varid_z , varid_t
        integer   :: varid_coords , coords_dimid , coords

        integer   :: varid_us(20) , varid_vs(20) , varid_ws(20)  !! MMAX
        integer   :: varid_rops(20)  , varid_ts(20) !! mmax
        integer   :: varid_thetam(20) !! mmax

        integer   :: varid_xg(20)  ! nmax(0)
        integer   :: varid_xs(20,20)  ! mmax , MAX(nmax(1:mmax))

        integer   :: varid_scalar(20)  ! nscalar
        integer   :: varid_rr(20)      ! nRR

        integer   :: varid_kturbg , varid_eturbg


        character(LEN=80) :: fname, var_name
        character(LEN=9) :: fname_index

        double precision, dimension(:) , allocatable :: arr1
        double precision, dimension(:) , allocatable :: arr2

        double precision, dimension(:,:,:,:) , allocatable :: arr4d


        double precision, dimension(:) , allocatable :: xloc
        double precision, dimension(:) , allocatable :: yloc
        double precision, dimension(:) , allocatable :: zloc

        double precision :: the_time
        logical          :: file_exists


! bWrite_netcdf(1)  : EP_g
! bWrite_netcdf(2)  : P_g
! bWrite_netcdf(3)  : P_star
! bWrite_netcdf(4)  : U_g / V_g / W_g
! bWrite_netcdf(5)  : U_s / V_s / W_s
! bWrite_netcdf(6)  : ROP_s
! bWrite_netcdf(7)  : T_g
! bWrite_netcdf(8)  : T_s
! bWrite_netcdf(9)  : X_g
! bWrite_netcdf(10) : X_s
! bWrite_netcdf(11) : Theta_m
! bWrite_netCDF(12) : Scalar
! bWrite_netCDF(13) : ReactionRates
! bWrite_netCDF(14) : k_turb_g , e_turb_g

        if (.not. mfix_usingnetcdf()) return
        if (.not. bglobalnetcdf) return

        call write_mesh_netcdf

        if (mype.eq.pe_io .and. .not.bdist_io) then
           allocate (arr1(ijkmax2))
           allocate (arr2(ijkmax3))
           allocate (arr4d(imax2,jmax2,kmax2,1))
           allocate ( xloc(imax2) )
           allocate ( yloc(jmax2) )
           allocate ( zloc(kmax2) )




        else
           allocate (arr1(1))
           allocate (arr2(1))
           allocate (arr4d(1,1,1,1))
           allocate ( xloc(1) )
           allocate ( yloc(1) )
           allocate ( zloc(1) )
       end if

!       CALL MPI_BARRIER(MPI_COMM_WORLD,mpierr)
        if (mype .ne. pe_io) goto 1234


        xloc(1) = -dx(1)
        do ii = 2,imax2
           xloc(ii) = xloc(ii-1) + dx(ii)
        end do


        yloc(1) = -dy(1)
        do ii = 2,jmax2
           yloc(ii) = yloc(ii-1) + dy(ii)
        end do

        zloc(1) = -dz(1)
        do ii = 2,kmax2
           zloc(ii) = zloc(ii-1) + dz(ii)
        end do




        if (bfirst_netcdf_write .and. mfix_usingnetcdf()) then
           bfirst_netcdf_write = .false.
           fname = trim(run_name) // '_netcdf_index.txt'
           inquire (file=fname,exist=file_exists)

           if (file_exists) then
               open (unit=11,file=fname,status='old')
               read (11,*) netcdf_file_index
               close (unit=11)
           else
               netcdf_file_index = 0
           end if
        end if

        fname_index = '_xxxxx.nc'
        write (fname_index(2:6),'(i5.5)') netcdf_file_index
        fname = trim(run_name)// fname_index
        call mfix_check_netcdf( mfix_nf90_create(fname, nf90_64bit_offset, ncid) )
        netcdf_file_index = netcdf_file_index + 1

        if (mfix_usingnetcdf()) then
                fname = trim(run_name) // '_netcdf_index.txt'
                open (unit=11,file=fname,status='unknown')
                write (11,*) netcdf_file_index
                close (unit=11)
        end if


        call mfix_check_netcdf( mfix_nf90_def_dim(ncid, "x"   , imax2   , x_dimid  ) )  ! 1
        call mfix_check_netcdf( mfix_nf90_def_dim(ncid, "y"   , jmax2   , y_dimid  ) )  ! 2
        call mfix_check_netcdf( mfix_nf90_def_dim(ncid, "z"   , kmax2   , z_dimid  ) )  ! 3
        call mfix_check_netcdf( mfix_nf90_def_dim(ncid, "coordinates"   , 1   , coords_dimid  ) )  ! 3
        call mfix_check_netcdf( mfix_nf90_def_dim(ncid, "t"   , 1       , t_dimid  ) )  ! 4


        call mfix_check_netcdf( mfix_nf90_def_var(ncid, "x"  , nf90_double, x_dimid, varid_x  ) )  ! 5
        call mfix_check_netcdf( mfix_nf90_def_var(ncid, "y"  , nf90_double, y_dimid, varid_y  ) )  ! 6
        call mfix_check_netcdf( mfix_nf90_def_var(ncid, "z"  , nf90_double, z_dimid, varid_z  ) )  ! 7
        call mfix_check_netcdf( mfix_nf90_def_var(ncid, "coordinates"  , nf90_int, coords_dimid, varid_coords  ) )  ! 7
        call mfix_check_netcdf( mfix_nf90_def_var(ncid, "t"  , nf90_double, t_dimid, varid_t  ) )  ! 8

        dimids =  (/ x_dimid , y_dimid, z_dimid , t_dimid/)

        if (bwrite_netcdf(1)) call mfix_check_netcdf( mfix_nf90_def_var(ncid, "EP_g"  , nf90_double, dimids, varid_epg  ) )  ! 9
        if (bwrite_netcdf(2)) call mfix_check_netcdf( mfix_nf90_def_var(ncid, "P_g"   , nf90_double, dimids, varid_pg  ) )   ! 10

        if (bwrite_netcdf(3)) call mfix_check_netcdf( mfix_nf90_def_var(ncid, "P_star", nf90_double, dimids, varid_pstar) )
        if (bwrite_netcdf(4)) call mfix_check_netcdf( mfix_nf90_def_var(ncid, "U_g"   , nf90_double, dimids, varid_ug   ) )
        if (bwrite_netcdf(4)) call mfix_check_netcdf( mfix_nf90_def_var(ncid, "V_g"   , nf90_double, dimids, varid_vg   ) )
        if (bwrite_netcdf(4)) call mfix_check_netcdf( mfix_nf90_def_var(ncid, "W_g"   , nf90_double, dimids, varid_wg   ) )
        if (bwrite_netcdf(7)) call mfix_check_netcdf( mfix_nf90_def_var(ncid, "T_g"   , nf90_double, dimids, varid_tg   ) )
        do i = 1,1   ! mmax
           var_name = 'U_s_xxx'
           write (var_name(5:7),'(i3.3)') i
           if (bwrite_netcdf(5)) call mfix_check_netcdf( mfix_nf90_def_var(ncid, var_name, nf90_double, dimids, varid_us(i)) )

           var_name = 'V_s_xxx'
           write (var_name(5:7),'(i3.3)') i
           if (bwrite_netcdf(5)) call mfix_check_netcdf( mfix_nf90_def_var(ncid, var_name, nf90_double, dimids, varid_vs(i)) )

           var_name = 'W_s_xxx'
           write (var_name(5:7),'(i3.3)') i
           if (bwrite_netcdf(5)) call mfix_check_netcdf( mfix_nf90_def_var(ncid, var_name, nf90_double, dimids, varid_ws(i)) )

           var_name = 'ROP_s_xxx'
           write (var_name(7:10),'(i3.3)') i
           if (bwrite_netcdf(6)) call mfix_check_netcdf( mfix_nf90_def_var(ncid, var_name, nf90_double, dimids, varid_rops(i)) )

           var_name = 'T_s_xxx'
           write (var_name(5:7),'(i3.3)') i
           if (bwrite_netcdf(8)) call mfix_check_netcdf( mfix_nf90_def_var(ncid, var_name, nf90_double, dimids, varid_ts(i)) )

           var_name = 'Theta_m_xxx'
           write (var_name(9:11),'(i3.3)') i
           if (bwrite_netcdf(11)) call mfix_check_netcdf( mfix_nf90_def_var(ncid, var_name, nf90_double, dimids, varid_thetam(i)) )

           DO nn = 1, nmax(i)
              var_name = 'X_s_xxx_xxx'
              write (var_name(5:7) ,'(i3.3)') i
              write (var_name(9:11),'(i3.3)') nn
              if (bwrite_netcdf(10)) call mfix_check_netcdf( mfix_nf90_def_var(ncid, var_name, nf90_double, dimids, varid_xs(i,nn)) )
           END DO


        end do

        do i = 1,nmax(0)
           var_name = 'X_g_xxx'
           write (var_name(5:7),'(i3.3)') i
           if (bwrite_netcdf(9)) call mfix_check_netcdf( mfix_nf90_def_var(ncid, var_name, nf90_double, dimids, varid_xg(i)) )
        end do

        do i = 1,nscalar
           var_name = 'Scalar_xxx'
           write (var_name(8:10),'(i3.3)') i
           if (bwrite_netcdf(12)) call mfix_check_netcdf( mfix_nf90_def_var(ncid, var_name, nf90_double, dimids, varid_scalar(i)) )
        end do

        do i = 1,nrr
           var_name = 'RRates_xxx'
           write (var_name(8:10),'(i3.3)') i
           if (bwrite_netcdf(13)) call mfix_check_netcdf( mfix_nf90_def_var(ncid, var_name, nf90_double, dimids, varid_rr(i)) )
        end do


        if (bwrite_netcdf(14) .and. k_epsilon) then
           call mfix_check_netcdf( mfix_nf90_def_var(ncid, 'k_turb_g', nf90_double, dimids, varid_kturbg) )
           call mfix_check_netcdf( mfix_nf90_def_var(ncid, 'e_turb_g', nf90_double, dimids, varid_eturbg) )
        end if



        call mfix_check_netcdf( mfix_nf90_enddef(ncid) ) ! 11

 1234   continue
           bfirst_netcdf_write = .false.

!       CALL MPI_BARRIER(MPI_COMM_WORLD,mpierr)


        if (mype .eq. pe_io) then
           call mfix_check_netcdf( mfix_nf90_put_var(ncid,varid_t,the_time) )   ! 12
           call mfix_check_netcdf( mfix_nf90_put_var(ncid,varid_x,xloc) )
           call mfix_check_netcdf( mfix_nf90_put_var(ncid,varid_y,yloc) )
           call mfix_check_netcdf( mfix_nf90_put_var(ncid,varid_z,zloc) )
           coords = 0
           if (coordinates .eq. 'CYLINDRICAL') coords = 1
           call mfix_check_netcdf( mfix_nf90_put_var(ncid,varid_coords,coords) )
        end if

        if (bwrite_netcdf(1)) then

            call gatherwritespx_netcdf(ep_g, arr1, arr2 , arr4d, ncid, varid_epg , &
                imax2,jmax2,kmax2,ijkmax2 , ijkmax3)

        end if

        if (bwrite_netcdf(2)) then

            call gatherwritespx_netcdf(p_g, arr1, arr2 , arr4d, ncid, varid_pg , &
                imax2,jmax2,kmax2,ijkmax2 , ijkmax3)

        end if




        if (bwrite_netcdf(3)) then

            call gatherwritespx_netcdf(p_star, arr1, arr2 , arr4d, ncid, varid_pstar , &
                imax2,jmax2,kmax2,ijkmax2 , ijkmax3)

        end if

        if (bwrite_netcdf(4)) then

            call gatherwritespx_netcdf(u_g, arr1, arr2 , arr4d, ncid, varid_ug , &
                imax2,jmax2,kmax2,ijkmax2 , ijkmax3)

            call gatherwritespx_netcdf(v_g, arr1, arr2 , arr4d, ncid, varid_vg , &
                imax2,jmax2,kmax2,ijkmax2 , ijkmax3)

            call gatherwritespx_netcdf(w_g, arr1, arr2 , arr4d, ncid, varid_wg , &
                imax2,jmax2,kmax2,ijkmax2 , ijkmax3)

        end if

        if (bwrite_netcdf(7)) then

            call gatherwritespx_netcdf(t_g, arr1, arr2 , arr4d, ncid, varid_tg , &
                imax2,jmax2,kmax2,ijkmax2 , ijkmax3)

        end if


        do i = 1,1   ! mmax

           if (bwrite_netcdf(5)) then

              call gatherwritespx_netcdf(u_s(:,i) , arr1, arr2 , arr4d, ncid, varid_us(i) , &
                   imax2,jmax2,kmax2,ijkmax2 , ijkmax3)

              call gatherwritespx_netcdf(v_s(:,i) , arr1, arr2 , arr4d, ncid, varid_vs(i) , &
                   imax2,jmax2,kmax2,ijkmax2 , ijkmax3)

              call gatherwritespx_netcdf(w_s(:,i) , arr1, arr2 , arr4d, ncid, varid_ws(i) , &
                   imax2,jmax2,kmax2,ijkmax2 , ijkmax3)

           end if

           if (bwrite_netcdf(6)) then

              call gatherwritespx_netcdf(rop_s(:,i) , arr1, arr2 , arr4d, ncid, varid_rops(i) , &
                   imax2,jmax2,kmax2,ijkmax2 , ijkmax3)

           end if

           if (bwrite_netcdf(8)) then

              call gatherwritespx_netcdf(t_s(:,i) , arr1, arr2 , arr4d, ncid, varid_ts(i) , &
                   imax2,jmax2,kmax2,ijkmax2 , ijkmax3)

           end if

           if (bwrite_netcdf(11)) then

              call gatherwritespx_netcdf(theta_m(:,i) , arr1, arr2 , arr4d, ncid, varid_thetam(i) , &
                   imax2,jmax2,kmax2,ijkmax2 , ijkmax3)

           end if

           if (bwrite_netcdf(10)) then
              do nn = 1,nmax(i)

                 call gatherwritespx_netcdf(x_s(:,i,nn) , arr1, arr2 , arr4d, ncid, varid_xs(i,nn) , &
                      imax2,jmax2,kmax2,ijkmax2 , ijkmax3)

              end do
           end if


        end do

        if (bwrite_netcdf(9)) then
           do i = 1,nmax(0)

               call gatherwritespx_netcdf(x_g(:,i) , arr1, arr2 , arr4d, ncid, varid_xg(i) , &
                   imax2,jmax2,kmax2,ijkmax2 , ijkmax3)

           end do
        end if

       if (bwrite_netcdf(12)) then
          do i = 1,nscalar

             call gatherwritespx_netcdf(scalar(:,i) , arr1, arr2 , arr4d, ncid, varid_scalar(i) , &
                   imax2,jmax2,kmax2,ijkmax2 , ijkmax3)

           end do
        end if


        if (bwrite_netcdf(13)) then
           do i = 1,nrr

             call gatherwritespx_netcdf(reactionrates(:,i) , arr1, arr2 , arr4d, ncid, varid_rr(i) , &
                   imax2,jmax2,kmax2,ijkmax2 , ijkmax3)

          end do
        end if

        if (bwrite_netcdf(14) .and. k_epsilon) then

           call gatherwritespx_netcdf(k_turb_g , arr1, arr2 , arr4d, ncid, varid_kturbg , &
                   imax2,jmax2,kmax2,ijkmax2 , ijkmax3)

           call gatherwritespx_netcdf(e_turb_g , arr1, arr2 , arr4d, ncid, varid_eturbg , &
                   imax2,jmax2,kmax2,ijkmax2 , ijkmax3)

        end if

        ! Close the file. This frees up any internal netCDF resources
        ! associated with the file, and flushes any buffers.
!       CALL MPI_BARRIER(MPI_COMM_WORLD,mpierr)
        if (mype .eq. pe_io) then
           call mfix_check_netcdf( mfix_nf90_close(ncid) )
        end if

        deallocate (arr1)
        deallocate (arr2)
        deallocate (arr4d)
        deallocate (xloc)
        deallocate (yloc)
        deallocate (zloc)

        return
        end subroutine write_netcdf