MODULE post
     USE Compact
     USE m_parameters
     USE m_calculate
     USE m_arrays
     USE m_terms

     IMPLICIT NONE
     PRIVATE
     INTEGER :: countnum
     INTEGER :: num_, dummyu_   ! hybrid
     REAL :: tnow

     REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: old_scalar, new_scalar

     PUBLIC :: main
    CONTAINS
 
     SUBROUTINE main
      INTEGER :: fread,i

      CALL READ_INTRO

      CALL ALLOCATE_ARRAYS

      CALL PRINT_BANNER

      countnum=0


      firstloop: DO fread=startnum,endnum,skipnum

        IF ( to_omit(fread) ) THEN

            WRITE(*,'(a47,i7,a4,i5,a3,i5)') &
                 ' Current fullsavenum = ', fread, ' || ', (fread-startnum+1), ' / ', (endnum-startnum+1)

            WRITE(*,'(a12,i6,a20,i6)') &
                 '      Skip. ', omit_t(i,1), ' <= fullsavenum <= ', omit_t(i,2)

        ELSE
         countnum=countnum+1

         CALL READ_FILE(fread)

      CALL m_terms_calculate_pass1

        ENDIF
          
      ENDDO firstloop

      write(*,*) '1st loop finished'


      CALL m_terms_average_pass1(countnum)


      secondloop: DO fread=startnum,endnum,skipnum

       IF ( to_omit(fread) ) THEN

         WRITE(*,'(a40,i7,a4,i4,a3,i4)') &
              ' Current fullsavenum = ', fread, ' || ', (fread-startnum+1), ' / ', (endnum-startnum+1)

         WRITE(*,'(a12,i6,a20,i6)') &
              '     Skip. ', omit_t(i,1), ' <= fullsavenum <= ', omit_t(i,2)

       ELSE

         CALL READ_FILE(fread)

      CALL m_terms_calculate_pass2

       ENDIF

      ENDDO secondloop

      CALL m_terms_average_pass2(countnum)
      
      CALL m_terms_write_result

      CALL DEALLOCATES_CLOSE


      WRITE(*,*) ' Avergaing RAW data is FINISHED'
      WRITE(*,*)  'qEdge_X.dat is generated' 
     END SUBROUTINE main

!========================================================================================
! End of main routine
!========================================================================================

     SUBROUTINE PRINT_BANNER
      WRITE(*,*) ' This program, x-edge-cold-bc-5-hybrid, is written by D. Kim, 2018'
      WRITE(*,*) ' It is to study the statistics of the flame parameters at the leading edge'
      WRITE(*,*) ' in turbulent premixed flames.'
      WRITE(*,'(a40,i5,a11,i5,a1)') '  Postprocess will be done from "FORT.',startnum,'" to "FORT.',endnum,'"'
     END SUBROUTINE PRINT_BANNER

     SUBROUTINE READ_FILE(num)
      INTEGER, INTENT(IN) :: num
      REAL, DIMENSION(2) :: tmpr
      INTEGER :: nx, ny, nz
      REAL :: tmp1,tmp2
      REAL :: dt,dummyu
      INTEGER :: ncyc

      OPEN(num,FORM='unformatted',STATUS='unknown')

      READ (num) tnow,nx,ny,nz,tmp1,tmp2

      IF ((nx .ne. nxp) .or. (ny .ne. nyp) .or. (nz .ne. nzp)) THEN
        WRITE(0,*) "Array dimension mismatch", nx, ny, nz, " != ", nxp, nyp, nzp
        STOP -1
      ENDIF

      READ (num) ncyc,dt,dummyu
      READ (num) tmpr(1:2)
      READ (num) tmpr(1:2)
      READ (num) tmpr(1:2)

      WRITE(*,'(a40,f8.3,a2,i7,a2,i5,a4,i5,a3,i5)') ' Current time / NCYC / fullsavenum = ',tnow,&
                ' /',ncyc,' /',num,' || ',(num-startnum+1),' / ',(endnum-startnum+1)
      WRITE(*,*) '    Reading current data file and processing'

      num_=num

      IF(num.le.shiftnum) THEN
        WRITE(*,*) '    with an old fort data from Nueman-0X'
        READ (num) u,v,w,old_scalar
        y(:,:,:) = old_scalar(2,:,:,:)
      ELSE
        WRITE(*,*) '    with a new fort data from Comb-Cluster'
        READ (num) u,v,w,new_scalar
        u = u + dummyu
        y = new_scalar(:,:,:,2)
      ENDIF

      CLOSE (num)

     END SUBROUTINE READ_FILE


!                SUBROUTINE SAVE_AVG_RESULTS
!                 INTEGER :: i
!                 OPEN (200,FILE="qEdge_X.dat") 
!           ! IRE1        
!                 WRITE(200,*) 'VARIABLES = "X","<C>","<FSD>","<Yr>","<Wc/rho>","<u>","<v>","<w>"' ! 8
!                 !WRITE(200,*) 'VARIABLES = "X","<C>","<FSD>","<Yr>","<Wc/rho>","<u>","<v>","<w>","<rho>"' ! 9
!           !      WRITE(200,*) '"<T>","<D>","<dc/dx>","<dc/dy>","<dc/dz>","<du/dx>","<dv/dy>","<dw/dz>"'   ! 8 -> 17
!           !      WRITE(200,*) '"<Div(rho*V)>","<c>_g","<d2c/dn2>"'  ! 3 -> 20
!           !      WRITE(200,*) '"d(1-<c>)/dx","d(1-<c>)/dy","d(1-<c>)/dz"'   ! 3 -> 23
!           !      WRITE(200,*) '"<Lap.(c)>","<d2c/dx2>","<d2c/dy2>","<d2c/dz2>"' ! 4 -> 27
!           !      WRITE(200,*) '"<Nx*(d2c/dn2)>","<d2(1-c)/dx2>"'  ! 2 -> 29
!           
!           ! IRE2          
!           !      WRITE(200,*) '"<FSD2>","<Sd>f","<Sdr>f","<Sdd>f","<Sdn>f","<Sdt>f","<Sdd2>f","<Sd2>f"'  ! 8 -> 8
!           !      WRITE(200,*) '"<DivN>f","<|DivN|>f","<(dc/dx)/c>f"'  ! 3 -> 11
!           !      WRITE(200,*) '"<-(dc/dn)/c>f","<(1/FSD`)*d(FSD`)/dx>f","<(1/FSD`)*d(FSD`)/dy>f"'  ! 3 -> 14
!           !      WRITE(200,*) '"<(1/FSD`)*d(FSD`)/dz>f","<Nx>f","<Ny>f","<Nz>f","<u>f","<v>f","<w>f"' ! 7 -> 21   
!           !      WRITE(200,*) '"<vn>f","<N dot N>f","<d(vn)/dn>f","<d(Sd)/dn>f","<d(Vn+Sd)/dn>f"' ! 5 -> 26
!           !      WRITE(200,*) '"<v>f dot <N>f","<v``dot N``>f","<N>f dot <N>f","<N``dot N``>f"'    ! 4 -> 30
!           !      WRITE(200,*) '"|<N>f|","<vn+Sd>f","d<Nx>f/dx"'    ! 3 -> 33
!           !      WRITE(200,*) '"-<N dot (grad.(FSD`))/FSD`>f"'    ! 1 -> 34
!           
!           ! IRE3          
!                 WRITE(200,*) '"<u>b","<v>b","<w>b","<rho>b","<u>b_g","<v>b_g","<w>b_g","<rho>b_g"'  ! 8 -> 8  
!           ! IRE4          
!                 WRITE(200,*) '"<u>u","<v>u","<w>u","<rho>u","<u>u_g","<v>u_g","<w>u_g","<rho>u_g"'  ! 8 -> 8 
!           ! IRE5          
!                 WRITE(200,*) '"RMS(u`)","RMS(v`)","RMS(w`)","<k>"'      ! 4 -> 4  
!                 WRITE(200,*) '"RMS(u`)b","RMS(v`)b","RMS(w`)b","<k>b"'  ! 4 -> 8   
!                 WRITE(200,*) '"RMS(u`)u","RMS(v`)u","RMS(w`)u","<k>u"'   ! 4 -> 12  
!                 !WRITE(200,*) '"<u`c`>","<v`c`>","<w`c`>"'  ! 3 -> 15 
!                 WRITE(200,*) '"<u`c`>","RMS(ux`)","RMS(uy`)"'  ! 3 -> 15 
!                 WRITE(200,*) '"RMS(uz`)"'  ! 1 -> 16 
!                 WRITE(200,*) '"RMS( U`)","RMS(ux`)_b","RMS(uy`)_b","RMS(uz`)_b"'      ! 4 -> 20 
!                 WRITE(200,*) '"RMS(U`)_b","RMS(ux`)_u","RMS(uy`)_u","RMS(uz`)_u"'      ! 4 -> 24 
!                 WRITE(200,*) '"RMS(U`)_u"'      ! 1 -> 25 
!           !
!           ! IRE6       
!                 WRITE(200,*) '"RMS(u`)_g","RMS(v`)_g","RMS(w`)_g","<k>_g"'      ! 4 -> 4
!                 WRITE(200,*) '"RMS(u`)b_g","RMS(v`)b_g","RMS(w`)b_g","<k>b_g"'  ! 4 -> 8  
!                 WRITE(200,*) '"RMS(u`)u_g","RMS(v`)u_g","RMS(w`)u_g","<k>u_g"'  ! 4 -> 12 
!                 WRITE(200,*) '"<u`c`>_g","<v`c`>_g","<w`c`>_g"'   ! 3 -> 15  
!           !! IRE7
!           !      WRITE(200,*) '"<vn>k<d2c/dn2>","<sd>k<d2c/dn2>","<vn+sd>k<d2c/dn2>"'   ! 3 -> 3  
!           !      WRITE(200,*) '"<vn>k","<sd>k","<vn+sd>k"'   ! 3 -> 6  
!           !! IRE8
!           !      WRITE(200,*) '"(1/<c>)*(d<c>/dx)","(1/<FSD>)*(d<FSD>/dx)","Dt_x","Dt_x_g"' ! 4 -> 4  
!           !      WRITE(200,*) '"1/Lw_3_High_Turb"'  ! 1 -> 5
!           !      WRITE(200,*) '"ST1","ST2","ST3","ST4","ST5","Lm*_x","Lm*_n","Lw","Lw_3"' ! 9 -> 14
!           !      WRITE(200,*) '"d<FSD>/dx","Dts"' ! 2 -> 16 
!           !      WRITE(200,*) '"(1/(1-<c>))*(d(1-<c>)/dx)"'   ! 1 -> 17
!           !      WRITE(200,*) '"1/L_LE_3","1/L_LE_4","ST6","ST7","ST8"' ! 5 -> 22
!           !      WRITE(200,*) '"1/L_LE_5=1/Lm-<DivN>f","ST_9"' ! 2 -> 24 
!           !      WRITE(200,*) '"1/L_LE_6=1/(d<c>/dx)*(d2<c>/dx2)","ST_10"' ! 2 -> 26 
!           !
!           !! IRE9
!           !      WRITE(200,*) '"<v>dotGrad.<c>","Dts*Lap.<c>","<Sd>f*<FSD>"' ! 3 -> 3
!           !      WRITE(200,*) '"Dts*<FSD>*Div.(<n>f)","<V`` dot N``>f*<FSD>","cEqnBalance"' ! 3 -> 6 
!           !      WRITE(200,*) '"<Nx>K","<Nx>f/<Nx>K","d<Dt>/dx"' ! 3 -> 9
!           !      WRITE(200,*) '"d<Dm+Dt>/dx","1/(d<1-c>/dx)*(d2<1-c>/dx2)"' ! 2 -> 11
!           !      WRITE(200,*) '"1/(1-<c>)*(d(1-<c>)/dx)","d(L_LE)/dx","d(L_TE)/dx"' ! 3 -> 14
!           
!           
!                 DO i=1,nxp
!           !       WRITE(200,'(156e20.10)') REAL(i)*hxp,IRE1(1:28,i),IRE2(1:34,i),IRE3(1:8,i),&  ! 1+28+34+8 = 71
!           !                       IRE4(1:8,i),IRE5(1:16,i),IRE6(1:15,i),IRE7(1:6,i),IRE8(1:26,i),IRE9(1:14,i)
!           !                                                                                  ! 8+16+15+6+26+14 = 85 -> 156
!                  WRITE(200,'(64e20.10)') REAL(i)*hxp,IRE1(1:7,i),IRE3(1:8,i),&  ! 1+7+8 =16 
!                                  IRE4(1:8,i),IRE5(1:25,i),IRE6(1:15,i)
!                                                                                  ! 8+25+15=48 -> 64
!                 ENDDO
!           
!           
!                 CLOSE(200) 
!                 !CLOSE(210) 
!              
!                END SUBROUTINE SAVE_AVG_RESULTS

     SUBROUTINE ALLOCATE_ARRAYS

      INTEGER :: ierr

       CALL m_arrays_init
       CALL m_calculate_init
       CALL m_terms_init

       ALLOCATE(old_scalar(2,nxp,nyp,nzp),STAT=ierr)  ; old_scalar=0.          ! Main variables
       ALLOCATE(new_scalar(nxp,nyp,nzp,2),STAT=ierr)  ; new_scalar=0.

       WRITE(*,'(a6,i3,a8,i3,a8,i3)') '  NX = ',nxp,' / NY = ',nyp,' / NZ = ',nzp
       WRITE(*,*) ' Preparing memory space for COMPACT SCHEME'
       ! CALL ludcmp(nxp,nyp,nzp,1,0,0)      ! 1,1,0

       WRITE(*,'(a22,i3,a3,i3,a4,i3)') '  Grid number range : ',syp,' ~ ',eyp,' of ',nyp
       WRITE(*,*)

     END SUBROUTINE ALLOCATE_ARRAYS

     SUBROUTINE DEALLOCATES_CLOSE

      CALL m_arrays_finalize
      CALL m_calculate_finalize
      CALL m_terms_finalize

      DEALLOCATE(old_scalar)
      DEALLOCATE(new_scalar)
      
      IF(omitnum.gt.0)   DEALLOCATE(omit_t)
     END SUBROUTINE DEALLOCATES_CLOSE

    END MODULE post