PROGRAM MAIN

          use spud
          use input, only : get_filename, &
                            read_in_options_from_xml!, & 
                            !read_in_numerical_options
          use output, only : write_results
          use boundary_intial_conditions, only : get_initial_conditions,  &
                                                 get_material_coefficients 

	  IMPLICIT NONE
         
          INTEGER I1,I2, N_1,N_2,NTIME,ITIME,II,ITS,NITS,NSUPCV,NSUBCV,NW,I,J,dump_frequency,ierror 

          REAL, ALLOCATABLE :: RHO(:,:),U1(:,:),U2(:,:),KAPPA(:,:),SIGMA(:,:),S(:,:),W(:,:)
          REAL, ALLOCATABLE :: TNP1(:,:),TN(:,:)
          ReaL, ALLOCATABLE :: T_SPATIAL(:,:)

          CHARACTER(LEN=256) filename, output_style


          REAL DT,DX_I1,DX_I2,XCORD,A,b,DX 
!          PARAMETER(N_1=12,N_2=12,NITS=50,NSUPCV=12,NW=7,NSUBCV=NW+1,dump_frequency=500)
!          REAL RHO(N_1,N_2),U1(N_1,N_2),U2(N_1,N_2),KAPPA(N_1,N_2),SIGMA(N_1,N_2),S(N_1,N_2),W(NSUPCV,NW)
!          REAL TNP1(N_1,N_2),TN(N_1,N_2)
!          ReaL T_SPATIAL(N_1-1,N_2-1)


       !---------------------------------------------------------------
       ! read in options from the input file    
       call get_filename(filename)
   
       call read_in_options_from_xml(N_1,N_2,NSUPCV,NW,NSUBCV,debug,output_style,&
                                     dump_frequency, NITS,NTIME, filename)

       !---------------------------------------------------------------
       ! allocate and initialise arrays
       ierror = 0
       allocate(TNP1(N_1,N_2),TN(N_1,N_2,RHO(N_1,N_2),U1(N_1,N_2),U2(N_1,N_2),KAPPA(N_1,N_2),SIGMA(N_1,N_2),S(N_1,N_2),W(NSUPCV,NW)), STAT = ierror)

       if (ierror/=0) then
            print *, "allocation error ... aborting"
            stop
       endif
       !---------------------------------------------------------------
       ! read in initial conditions
       call get_initial_conditions(N_1,N_2,TN,TNP1,filename)
       !---------------------------------------------------------------
       ! read in material parameters (source, absorption coeff, diffusion coeff)
       !allocate(S_C(NSUPCV),STAT=ierror) ! NOT USED
       allocate(S(N_1,N_2),STAT=ierror)
       allocate(SIGMA(N_1,N_2),STAT=ierror)
       allocate(KAPPA(N_1,N_2),STAT=ierror)
       allocate(RHO(N_1,N_2),STAT=ierror)

       call get_material_coefficients(S,SIGMA,KAPPA,N_1,N_2,filename)




! Space and time step sizes
         DT=3./120.0
  
         DX_I1 = 1./REAL(N_1)
         DX_I2 = 1./REAL(N_2)
         NTIME=INT(5.0/DT)
         PRINT *,'DT,NTIME',DT,NTIME,DX_I1

! Inititional conditions
         TN=0.0; TNP1=0.0
         TN(5:6,5:6)=1.0
         TNP1(5:6,5:6)=1.0
         PRINT *,'TNP1=',TNP1(:,5:6)
! Parameter values
         RHO=1.0
         U1=0.0
         U2=0.0
         KAPPA=1.0
         SIGMA=0.0
         S=0.0

!***** material properties
! Diffusion Coefficients
!!H2O
!         DO I1=1,FLOOR(N_1/5.0)
!            DO I2=1,FLOOR(N_2/5.0)
!                KAPPA(I1,I2)=0.1587
!            END DO 
!         END DO
!! FUEL
!         DO I1=FLOOR(N_1/5.0)+1,2*FLOOR(N_1/5.0)
!             DO I2=2*FLOOR(N_2/5.0)+1,N_2
!                KAPPA(I1,I2)=0.1205
!             END DO 
!         END DO
!! H2O
!         DO I1=2*FLOOR(N_1/5.0)+1,3*FLOOR(N_1/5.0)
!             DO I2=2*FLOOR(N_2/5.0)+1,N_2
!                KAPPA(I1,I2)=0.1587
!             END DO 
!         END DO
!! FUEL
!         DO I1=3*FLOOR(N_1/5.0)+1,4*FLOOR(N_1/5.0)
!            DO I2=FLOOR(N_2/5.0)+1,FLOOR(N_2/5.0)
!               KAPPA(I1,I2)=0.1205
!            END DO
!         END DO
!! H2O
!         DO I1=4*FLOOR(N_1/5.0)+1,N_1
!             DO I2=2*FLOOR(N_2/5.0)+1,N_2
!                KAPPA(I1,I2)=0.1587
!             END DO 
!         END DO
!!Source
!    
!!H2O
!         DO I1=1,FLOOR(N_1/5.0)
!            DO I2=1,FLOOR(N_2/5.0)
!                S(I1,I2)=0.0
!            END DO 
!         END DO
!! FUEL
!         DO I1=FLOOR(N_1/5.0)+1,2*FLOOR(N_1/5.0)
!             DO I2=2*FLOOR(N_2/5.0)+1,N_2
!                S(I1,I2)=1.0
!             END DO 
!         END DO
!! H2O
!         DO I1=2*FLOOR(N_1/5.0)+1,3*FLOOR(N_1/5.0)
!             DO I2=2*FLOOR(N_2/5.0)+1,N_2
!                S(I1,I2)=0.0
!             END DO 
!         END DO
!! FUEL
!         DO I1=3*FLOOR(N_1/5.0)+1,4*FLOOR(N_1/5.0)
!            DO I2=FLOOR(N_2/5.0)+1,FLOOR(N_2/5.0)
!               S(I1,I2)=1.0
!            END DO
!         END DO
!! H2O
!         DO I1=4*FLOOR(N_1/5.0)+1,N_1
!             DO I2=2*FLOOR(N_2/5.0)+1,N_2
!                S(I1,I2)=0.0
!             END DO 
!         END DO
!!Cross sections
!!H2O
!        SIGMA=0.0
!         DO I1=1,FLOOR(N_1/5.0)
!            DO I2=1,FLOOR(N_2/5.0)
!                SIGMA(I1,I2)=0.06284
!            END DO 
!         END DO
!! FUEL
!         DO I1=FLOOR(N_1/5.0)+1,2*FLOOR(N_1/5.0)
!             DO I2=2*FLOOR(N_2/5.0)+1,N_2
!                SIGMA(I1,I2)=0.05822
!             END DO 
!         END DO
!! H2O
!         DO I1=2*FLOOR(N_1/5.0)+1,3*FLOOR(N_1/5.0)
!             DO I2=2*FLOOR(N_2/5.0)+1,N_2
!                SIGMA(I1,I2)=0.06284
!             END DO 
!         END DO
!! FUEL
!         DO I1=3*FLOOR(N_1/5.0)+1,4*FLOOR(N_1/5.0)
!            DO I2=FLOOR(N_2/5.0)+1,FLOOR(N_2/5.0)
!               SIGMA(I1,I2)=0.05822
!            END DO
!         END DO
!! H2O
!         DO I1=4*FLOOR(N_1/5.0)+1,N_1
!             DO I2=2*FLOOR(N_2/5.0)+1,N_2
!                SIGMA(I1,I2)=0.06284
!             END DO 
!         END DO
! ****end of material properties






! Time loop
         DO ITIME=1, NTIME !10
         DO ITS=1, NITS
! Space loop
            DO I1=2, N_1-1
            DO I2=2,N_2-1 

!            PRINT *,U(I1,I2)
            A = ( (RHO(I1,I2)/DT) + (0.5*(U1(I1-1,I2)+U1(I1,I2)))/DX_I1 + (0.5*(U2(I1,I2-1)+U2(I1,I2)))/DX_I2 )  + SIGMA(I1,I2)  &
             + ( (0.5*(KAPPA(I1+1,I2)+KAPPA(I1,I2)))/DX_I1 + (0.5*(KAPPA(I1-1,I2)+KAPPA(I1,I2)))/DX_I1 ) /DX_I1 &
              + ( (0.5*(KAPPA(I1,I2+1)+KAPPA(I1,I2)))/DX_I2 + (0.5*(KAPPA(I1,I2-1)+KAPPA(I1,I2)))/DX_I2 ) /DX_I2 

            b = TN(I1,I2)*(RHO(I1,I2)/DT) + (0.5*(U1(I1-1,I2)+U1(I1,I2))/DX_I1) *TNP1(I1-1,I2) +(0.5*(U2(I1,I2-1)+U2(I1,I2))/DX_I2)*TNP1(I1,I2-1) + S(I1,I2) &
              + ( (0.5*(KAPPA(I1+1,I2)+KAPPA(I1,I2)))*TNP1(I1+1,I2)/DX_I1 + (0.5*(KAPPA(I1-1,I2)+KAPPA(I1,I2)))*TNP1(I1-1,I2)/DX_I1 ) /DX_I1 &
              + ( (0.5*(KAPPA(I1,I2+1)+KAPPA(I1,I2)))*TNP1(I1,I2+1)/DX_I2 + (0.5*(KAPPA(I1,I2-1)+KAPPA(I1,I2)))*TNP1(I1,I2-1)/DX_I2 ) /DX_I2  

!            PRINT *, I1,I2,'A',A, 'b',b
            TNP1(I1,I2)=b/A
!           PRINT *,I1,I2,TN(I1,I2),TNP1(I1,I2)

! no Diff

!            A = ( (RHO(I1,I2)/DT) + (0.5*(U(I1+1,I2)-U(I1,I2)))/DX_I1 + (0.5*(U(I1,I2-1)-U(I1,I2)))/DX_I2 - SIGMA(I1,I2) ) 

!            b = TN(I1,I2)*(RHO(I1,I2)/DT) + (0.5*(U(I1+1,I2)-U(I1,I2))/DX_I1) *TNP1(I1,I2) +(0.5*(U(I1,I2-1)-U(I1,I2))/DX_I2)*TNP1(I1,I2) + S(I1,I2)

!            TNP1(I1,I2)=b/A
                              
            END DO ! End of I2=1,N_2
            END DO ! End of I1=1,N_1
!            stop 234
!       IF (ITS==49)  THEN
!          DO I1=2, N_1
!          DO I2=5,5
!          PRINT *,I1,I2,TN(I1,I2),TNP1(I1,I2)
!          END DO
!          END DO
!       END IF
         END DO ! END of Its
            TN = TNP1

             ! prepare for output
             T_SPATIAL = 0.0   
             DO I1=1,N_1
             DO I2=1,N_2
                   T_SPATIAL(I1,I2) = T_SPATIAL(I1,I2) + TNP1(I1,I2)
             END DO
             END DO
!            PRINT *,T_SPATIAL


             if (MOD(ITIME,int(NTIME/dump_frequency))==0) then
                 print *, "writing results at ITIME=", ITIME
                 call write_results(DX_I1,N_1,N_2,TNP1,T_SPATIAL,ITIME,output_style) 

             endif 



         END DO ! End of time loop


! Print results to terminal
!           PRINT *,I1,I2,TN,TNP1
             XCORD = 0.5*DX_I1
             DO I1=1,N_1
             DO I2=2,2 !NSUBCV_2
                II=II+1
!                PRINT *,II,TCV_NEW(I,K) 
                PRINT *,XCORD, TNP1(I1,I2),T_SPATIAL 
                XCORD = XCORD + DX_I1
             END DO
             END DO


!! Time loop
!         DO ITIME=1,NTIME
!! Space loop
!            DO I1=2, N_1-1
!            DO I2=2,N_2-1      
!               TNP1(I1,I2) = (DT/RHO(I1,I2)) * (TN(I1,I2) &
!                                 + (0.5*(U(I1+1,I2)-U(I1,I2)))*(TNP1(I1,I2)-TNP1(I1-1,I2))/DX_I1 + (0.5*(U(I1,I2-1)-U(I1,I2)))*(TNP1(I1,I2)-TNP1(I1,I2-1))/DX_I2  &
!                                - SIGMA(I1,I2) * TN(I1,I2) - S(I1,I2) &
!                                 +((0.5*(KAPPA(I1+1,I2)-KAPPA(I1,I2))) * (TN(I1+1,I2)-TN(I1,I2))/DX_I1  - (0.5*(KAPPA(I1-1,I2)-KAPPA(I1,I2)))  * (TN(I1,I2)-TN(I1-1,I1))/DX_I1  ) /DX_I1 &  
!                               +((0.5*(KAPPA(I1,I2+1)-KAPPA(I1,I2))) * (TN(I1,I2+1)-TN(I1,I2))/DX_I2  - (0.5*(KAPPA(I1,I2-1)-KAPPA(I1,I2)))  * (TN(I1,I2)-TN(I1,I1-1))/DX_I2  ) /DX_I2 )  
!                              
!                               ! + (0.5*(KAPPA(I1,I2+1)-KAPPA(I1,I2)))  * (TN(I1,I2-1)-TN(I1,I2)/DX_I2)) - (0.5*(KAPPA(I1,I2)-KAPPA(I1,I2-1)))  * (TN(I1,I2)-TN(I1,I2-1)/DX_I2)) /DX_I2 ) 
!            END DO ! End of I2=1,N_2
!            END DO ! End of I1=1,N_1
!            TN = TNP1
!         END DO ! End of time loop

    DO I=2,NSUPCV-1
    DO J=1,NW
          W(I,J)=DX(I,J,DX_I1)
    END DO
    END DO
!    PRINT *,'w= ', W

         END PROGRAM MAIN



        
          
          REAL FUNCTION DX(I,J,DX_I)
	  IMPLICIT NONE
          INTEGER I,J
          REAL DX_I
! Local variables...
          INTEGER HW ! the function.
          INTEGER NO_WAVELETS_ACROSS

          NO_WAVELETS_ACROSS = 2**(HW(J)-1)
          DX = DX_I /REAL( NO_WAVELETS_ACROSS )
          
          RETURN
          END FUNCTION DX




          INTEGER FUNCTION JW(K,L,NO_W_LEVELS)
	  IMPLICIT NONE
          INTEGER K,L,NO_W_LEVELS
! local variables...
          INTEGER PRODUCT_NO ! the function

          JW = INT( (K-1)/ 2**(NO_W_LEVELS-L+1)  +1 )    + PRODUCT_NO(L-1)
          RETURN
          END FUNCTION JW

       


          INTEGER FUNCTION PRODUCT_NO(II)
! returns the triangular number of II
	  IMPLICIT NONE
          INTEGER II
! local variables...
          INTEGER K

          PRODUCT_NO = 0
          DO K=1,II
             PRODUCT_NO = PRODUCT_NO + 2**(K-1) 
          END DO
          RETURN
          END FUNCTION PRODUCT_NO
 




          REAL FUNCTION SW(K,L,NO_W_LEVELS)
	  IMPLICIT NONE
          INTEGER K,L,NO_W_LEVELS

!          SW = 2*MOD( INT(K/2**(NO_W_LEVELS-1)), 2) -1
          SW = 2.0*(  2*INT((K-1)/2**(NO_W_LEVELS-L+1) +1) - INT((K-1)/2**(NO_W_LEVELS-L) +1) )  -1.0
          RETURN
          END FUNCTION SW




          INTEGER FUNCTION HW(J) 
	  IMPLICIT NONE
          INTEGER J
! local variables...
          INTEGER L,LL,L_KEEP, JJ
          INTEGER PRODUCT_NO ! the function
!          JJ=J
!          JJ=0
          DO L=10,1,-1
             L_KEEP=L
             LL = J -PRODUCT_NO(L-1)
!             PRINT *,'J,L-1,PRODUCT_NO(L-1),ll:',J,L-1,PRODUCT_NO(L-1),ll
!             IF(LL.GE.0) CYCLE
             IF(LL.GT.0) EXIT
          END DO
          IF(J==0) L_KEEP =0
!          print *,'L_KEEP:',L_KEEP
!          STOP 45
          HW = L_KEEP
          RETURN
          END FUNCTION HW