PROGRAM MAIN
	  IMPLICIT NONE
	  INTEGER NSUPCV,NW,NSUBCV,NTIME,NITS,SAMPLE
          REAL THETA
          LOGICAL UPWIND
! if use_g_for_pert then use g in the calculation of purterbations.
! NSUPCV=no of course grid or super control volumes 
! NSUBCV=no of fine grid CVs per supercontrol volume. 
! NW=no of wavelets within each course grid node I (NW=NO_W_LEVELS! so not free to choose any NW)
! NITS=no of non-linear iterations used for anything other than THETA=0.0 NITS>1.
!          PARAMETER(NSUPCV=10,NW=7,NSUBCV=NW+1,NTIME=100*0.8,NITS=1,UPWIND=.false.)
!          PARAMETER(NSUPCV=10,NW=7,NSUBCV=NW+1,NTIME=100,NITS=1,UPWIND=.false.)
!          PARAMETER(NSUPCV=10,NW=3,NSUBCV=NW+1,NITS=1,UPWIND=.false.)
          PARAMETER(NSUPCV=100,NW=1,NSUBCV=NW+1,NITS=10,UPWIND=.false.)
!          PARAMETER(NSUPCV=10,NW=15,NSUBCV=NW+1,NTIME=100,NITS=1,UPWIND=.false.)
!          PARAMETER(NSUPCV=10,NW=15,NSUBCV=NW+1,NITS=100,UPWIND=.false.)
!          PARAMETER(NSUPCV=10,NW=31,NSUBCV=NW+1,NTIME=100,NITS=1,UPWIND=.false.)
!          PARAMETER(NSUPCV=10,NW=31,NSUBCV=NW+1,NTIME=100*50/15,NITS=1,UPWIND=.false.)
! THETA=0.0 is forward Euler, =0.5 CranckNickolson, =1.0 Backward Euler. 
          PARAMETER(THETA=0.0)
          PARAMETER(SAMPLE = 2) 

       REAL TCV_OLD(NSUPCV,NSUBCV), TCV_NEW(NSUPCV,NSUBCV)
       REAL TSGS_OLD(NSUPCV,NW), TSGS_NEW(NSUPCV,NW)
       REAL TCV(NSUPCV,NSUBCV)
       REAL TC_OLD(NSUPCV), TC_NEW(NSUPCV)
       REAL S_C(NSUPCV), SIGMA_C(NSUPCV)
       REAL S_SGS(NSUPCV,NSUPCV), SIGMA_SGS(NSUPCV,NSUPCV)


! Functions...
       INTEGER JW,HW
       INTEGER MPL, MPR, MMR
       REAL SW, DX, XCORD
! Local variables...
       INTEGER MPL_J, MPR_J, MMR_J ! functions
       INTEGER MML_J
       INTEGER I,J,K,L, NO_W_LEVELS, V,V_M,U_M, ITS, ITIME, II,Q,unit_n
       INTEGER MPL_JP, MPR_JP, MML_JP, MMR_JP
       INTEGER W_M_RIGHT
       REAL DT, DX_I,U,D(NSUPCV,NSUBCV),V_N,DS
       real v_diff_left, v_diff_right, dx_small
       real v_diff_left_plus, v_diff_right_plus, v_diff_left_minus, v_diff_right_minus
! 


! Source and sigma:
       S_C =0.0; SIGMA_C  =0.0
       S_SGS=0.0; SIGMA_SGS=0.0
!       S_C(8) =20.0;  S_SGS(8,8)=20.0
! Adevction Speed
         U=2.0
! Diffusion Coefficients
         D=0.01
! Thermal mean speed
!         V_N=2.2e-4
         V_N=1
                            

! Space and time step sizes
!         DT=0.0001 /REAL(NSUPCV*NSUBCV) 
         DT=0.1 /REAL(NSUPCV*NSUBCV) 
!         DT=0.5 /REAL(NSUPCV*NSUBCV) 
         DX_I = 1./REAL(NSUPCV)
         NTIME=INT(0.5/DT)
!         NTIME=INT(0.005/DT)
         PRINT *,'DT,NTIME',DT,NTIME
!         STOP 22

! Inititional conditions...
         TCV_OLD=0.0; TCV_NEW=0.0
         TSGS_OLD=0.0; TSGS_NEW=0.0
         TC_OLD=0.0; TC_NEW=0.0
!         TC_OLD(3)=1.0
!         TC_NEW(3)=1.0
!   DO I=1, NSUPCV
!      TC_OLD(I)=exp(-Real(I-10)**2/6.0) 
!      TC_NEW(I)=exp(-Real(I-10)**2/6.0) 
!   END DO

   DO I=1+1, NSUPCV-1
      TC_OLD(I)=sin(2*3.14159*Real(I)/100.0)
      TC_NEW(I)=sin(2*3.14159*Real(I)/100.0)
   END DO
!         TSGS_OLD(3,3)=1.0; TSGS_NEW(3,3)=1.0


!         PRINT *,'TC_NEW=',TC_NEW
         NO_W_LEVELS = HW(NW)
!          print *,'NO_W_LEVELS,nw:',NO_W_LEVELS,nw
!          stop 2828

! M calculate:
	unit_n = 1
	OPEN (unit=unit_n,file="./snapshots_burgerc.dat",action="write",status="replace")
         DO ITIME=1,NTIME

         DO ITS=1,NITS

! MAP to CV space: 
            DO I=1,NSUPCV
            DO K=1,NSUBCV
               TCV_OLD(I,K) = TC_OLD(I) 
               TCV_NEW(I,K) = TC_NEW(I) 
               DO L=1,NO_W_LEVELS 
                  TCV_OLD(I,K) = TCV_OLD(I,K) + SW(K,L,NO_W_LEVELS) * TSGS_OLD(I,JW(K,L,NO_W_LEVELS))
                  TCV_NEW(I,K) = TCV_NEW(I,K) + SW(K,L,NO_W_LEVELS) * TSGS_NEW(I,JW(K,L,NO_W_LEVELS))
!                  print *,'i,k,l,SW(K,L,NO_W_LEVELS):',i,k,l,SW(K,L,NO_W_LEVELS)
               END DO
            END DO
            END DO
!            PRINT *,TCV_OLD(1:10,1:3)
!             stop 991
!            PRINT *,'HERE1'
! Use theta time stepping:
            TCV = THETA * TCV_NEW + (1.0 - THETA) * TCV_OLD
!            print *,'NO_W_LEVELS,nw:',NO_W_LEVELS,nw
!            print *,'tcv_old:',tcv_old
!            print *,'tcv:',tcv
!            stop 229

! COARSE GRID...
            DO I=2, NSUPCV-1
               v_diff_left  = V_N* 0.5* (D(I,1)+D(I-1,NSUBCV))
               v_diff_right = V_N* 0.5* (D(I+1,1)+D(I,NSUBCV))
               dx_small  = DX_I/real( NSUBCV)
               TC_NEW(I) = TC_OLD(I) + V_N*U*(DT/DX_I)* (0.5*TCV(I-1,NSUBCV)**2 - 0.5*TCV(I,NSUBCV)**2) &
                                    -v_diff_left  * (DT/DX_I) *  ((TCV(I,1)-TCV(I-1,NSUBCV)) /dx_small) &
                                    +v_diff_right*(DT/DX_I)*     ((TCV(I+1,1)-TCV(I,NSUBCV)) /dx_small) &
                                    +DT*V_N*SIGMA_C(I)*TC_NEW(I) &
                                    +DT*V_N*S_C(I)
            END DO
!            PRINT *,'HERE2'

! FINE GRID (WAVELET)...
            DO I=2, NSUPCV-1
               DO J=1, NW
!                  PRINT *,'1-I,J=',I,J
                  MPL_J = MPL(J,NSUBCV) 
                  MPR_J = MPR(J,NSUBCV) 
                  MML_J = MPR_J 
                  MMR_J = MMR(J,NSUBCV)
                  IF( MPL_J==0 ) THEN
                     MPL_J= NW+1
                     V=I-1
                  ELSE
                     V=I
                  ENDIF

                  MPL_JP = MPL_J+1
                  MPR_JP = MPR_J+1
                  MML_JP = MML_J+1
                  MMR_JP = MMR_J+1

                  IF( MPL_JP==NSUBCV+1 ) THEN
                     MPL_JP= 1
                  ENDIF
                  IF( MMR_JP==NSUBCV+1 ) THEN
                     MMR_JP= 1
                     W_M_RIGHT=I+1
                  ELSE
                     W_M_RIGHT=I
                  ENDIF
!                  PRINT *,'2-I,J=',I,J
!                  print *,'i,j,v,MPL_J,MPR_J,MML_J,MMR_J,DX(I,J,DX_I):',i,j,v,MPL_J,MPR_J,MML_J,MMR_J,DX(I,J,DX_I)
               v_diff_left_plus  = V_N* 0.5* (D(I,MPL_JP)+D(V,MPL_J))
               v_diff_right_plus = V_N* 0.5* (D(I,MPR_JP)+D(I,MPR_J))
               v_diff_left_minus = V_N* 0.5* (D(I,MML_JP)+D(I,MML_J))
               v_diff_right_minus = V_N* 0.5* (D(W_M_RIGHT,MMR_JP)+D(I,MMR_J))
               dx_small  = DX_I/real( NSUBCV)

                  TSGS_NEW(I,J) = TSGS_OLD(I,J) + V_N*U*(DT/DX(I,J,DX_I))*(0.5*TCV(V,MPL_J)**2 - 0.5*TCV(I,MPR_J)**2 - 0.5*TCV(I,MML_J)**2 + 0.5*TCV(I,MMR_J)**2) &
                                    -v_diff_left_plus * (DT/DX(I,J,DX_I)) * (TCV(I,MPL_JP)-TCV(V,MPL_J))  /dx_small &
                                    +v_diff_right_plus* (DT/DX(I,J,DX_I)) * (TCV(I,MPR_JP)-TCV(I,MPR_J))/dx_small &
                                    +v_diff_left_minus* (DT/DX(I,J,DX_I)) * (TCV(I,MML_JP)-TCV(I,MML_J))/dx_small &
                                    -v_diff_right_minus*(DT/DX(I,J,DX_I)) * (TCV(W_M_RIGHT,MMR_JP)-TCV(I,MMR_J))  /dx_small &
                                    -DT*V_N*SIGMA_SGS(I,J)*TSGS_NEW(I,J) &
                                    +DT*V_N*S_SGS(I,J)
!                  print *,'NSUBCV,i,j:',NSUBCV,i,j
!                  print *,'I,MPL_JP,V,MPL_J:',I,MPL_JP,V,MPL_J
!                  print *,'I,MPR_JP,I,MPR_J:',I,MPR_JP,I,MPR_J
!                  print *,'I,MML_JP,I,MML_J:',I,MML_JP,I,MML_J
!                  print *,'W_M_RIGHT,MMR_JP,I,MMR_J:',W_M_RIGHT,MMR_JP,I,MMR_J
!                  print *,'W_M_RIGHT:',W_M_RIGHT

               END DO ! END DO J=1, NW-1
!               print *,' '
            END DO ! END DO I=2, NSUPCV-1
!       stop 22991
!            PRINT *,'HERE3'

         END DO ! DO ITS=1,NITS
!          print *,'here4'
            TSGS_OLD = TSGS_NEW
            TC_OLD = TC_NEW
!          print *,'here5'



!       Write data to file to be passed to python
!       Take sample of time steps
        
        DO Q=1,SAMPLE
             IF (MODULO(ITIME, SAMPLE) == 0) THEN
             DO I=1,NSUPCV
             DO K=1,NSUBCV 
!                PRINT *, TCV_NEW(I,K)
!                WRITE(unit_n,*) TSGS_NEW(I,K)
                 WRITE(unit_n,*) TC_NEW(I)
             END DO
             END DO
             END IF
       END DO
             


         END DO ! DO ITIME=1,NTIME
       CLOSE (unit_n)

!          print *,'dx_i,DX_I/real(NSUBCV):', dx_i,DX_I/real(NSUBCV)
!          stop 819
!             XCORD = 0.5*DX_I/real(NSUBCV)
!             DO I=1,NSUPCV
!             DO K=1,NSUBCV
!                II=II+1
!!                PRINT *,II,TCV_NEW(I,K) 
!                PRINT *,XCORD, TCV_NEW(I,K) 
!                XCORD = XCORD + DX_I/real(NSUBCV)
!             END DO
!             END DO

	  STOP
	  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




          INTEGER FUNCTION MPL(J,NSUBCV) 
	  IMPLICIT NONE
          INTEGER J,NSUBCV
! local variables...
! Funcitons...
          INTEGER HW ! the function
          INTEGER PRODUCT_NO ! the function
!            PRINT *,'INSIDE MPL J,NCV:',J,NCV
!            PRINT *,'HW(J):',HW(J)
!            PRINT *,'PRODUCT_NO(HW(J)-1):',PRODUCT_NO(HW(J)-1)
             
        
          MPL = INT(  ((2*(J-PRODUCT_NO(HW(J)-1))-2)*NSUBCV)/(2**HW(J))  )
!            PRINT *,'J,NCV,MPL:',J,NCV,MPL
          RETURN
          END FUNCTION MPL



          INTEGER FUNCTION MPR(J,NSUBCV) 
	  IMPLICIT NONE
          INTEGER J,NSUBCV
! Funcitons...
          INTEGER HW ! the function
          INTEGER PRODUCT_NO ! the function
        
          MPR = INT(  ((2*(J-PRODUCT_NO(HW(J)-1))-1)*NSUBCV)/(2**HW(J))  )
          RETURN
          END FUNCTION MPR



          INTEGER FUNCTION MMR(J,NSUBCV) 
	  IMPLICIT NONE
          INTEGER J,NSUBCV
! Funcitons...
          INTEGER HW ! the function
          INTEGER PRODUCT_NO ! the function
        
          MMR = INT(  ((2*(J-PRODUCT_NO(HW(J)-1))-0)*NSUBCV)/(2**HW(J))  )
          RETURN
          END FUNCTION MMR