!*******************************************
!*
!* Solve the advection PDE using Explicit FTCS, 
!* Explicit Lax, Implicit FTCS, and implicit Crank-Nicolson 
!* methods for constant and varying speed.
!*
!*  Solve  dc/dt = -u dc/dx
!*  u = t/20 ft/minute
!*  and
!*  u constant
!* 
!* Compiler used: gnu 95 (g95) on Cygwin. Gcc 3.4.4
!* Date: June 20 2006
!* 
!* by Nasser Abbasi
!*******************************************

PROGRAM advection
  IMPLICIT NONE

  REAL     :: DT,DX,max_run_time,length,snapshot_delta, &
              first_limit,second_limit
  INTEGER  :: N,SNAPSHOTS 
  character(10) :: cmd_arg ! to read time step from command line

  INTEGER :: method  ! 1=FTCS, 2=LAX, 3=Implicit FTCS, 4=C-R
  INTEGER :: mode    ! 1=Fixed wind speed, 2=speed function of time

  REAL    :: t_start, t_end, cpu_time_used,end_line(1002)
  INTEGER :: ALL_DATA_FILE_ID 
  PARAMETER(ALL_DATA_FILE_ID=900)

  ! Initialize data. All methods will use the same
  ! parameters to make comparing them easier

  ! read delta t from command line.                        
  CALL getarg(1,cmd_arg)
  cmd_arg=TRIM(cmd_arg)                  
  print *,'= ', cmd_arg           
  read(cmd_arg,*)dt   !delta in time,  in minutes

  print *,'Dt=',DT                   

  N      = 1000      ! number of grid points in space           
  length = 100       ! length of space solution in feet

  first_limit  = 0.25*length
  second_limit = 0.35*length

  DX = length/N     ! delta in space, in feets

  max_run_time = 30.0  ! how long to run for in minutes
  SNAPSHOTS    = 200   ! number of snapshots per run. Used for animation

  snapshot_delta = max_run_time / SNAPSHOTS   ! time between each snap shot

  print *,'DT=',DT,' minutes, DX=',DX,' feets'
  print *,'taking snapshots every ', snapshot_delta ,' minutes'     

  DO mode=1,2
     print*,'=======> processing mode ',mode
     DO method=1,4   ! No enumeration data types in Fotran 90

        CALL CPU_TIME(t_start)  ! get current CPU time
        CALL process(mode,method,N,DT,DX,max_run_time,snapshot_delta,&
             first_limit,second_limit)
        CALL CPU_TIME(t_end)    ! get current CPU time

        cpu_time_used = t_end - t_start                                    

        WRITE(*,FMT='(A,I2,A,F12.5)') 'CPU TIME used for method', method, ' = ', cpu_time_used
        ! Now record test case parameters in last line
        end_line=0
        end_line(1)=cpu_time_used
        end_line(2)=DT
        end_line(3)=DX
        end_line(4)=mode
        end_line(5)=method
        
        WRITE(UNIT=ALL_DATA_FILE_ID,FMT=*) end_line
        CLOSE(ALL_DATA_FILE_ID)

     END DO
  END DO

END PROGRAM advection
!************************************
!*
!*
!************************************
SUBROUTINE process(mode,method,N,DT,DX,max_run_time,snapshot_delta,&
                   first_limit,second_limit)
  IMPLICIT NONE

  INTEGER, INTENT(IN) :: mode,method,N                   
  REAL,    INTENT(IN) :: DT,DX,max_run_time,snapshot_delta,&
                         first_limit,second_limit

  INTEGER :: I
  LOGICAL :: snap_shot_at_15_taken
  INTEGER :: ALL_DATA_FILE_ID 
  PARAMETER(ALL_DATA_FILE_ID=900)
  REAL    :: snap_current_time           
  REAL    :: current_time
  REAL    :: C(N)       ! current solution
  REAL    :: CNEW(N)    ! future solution
  REAL    :: CEXACT(N)  ! current exact solution
  REAL    :: current_first_limit
  REAL    :: A(N,N),aa(N),b(2:N),cc(N-1),CTEMP(N)  ! for C-R and implicit FTCS
  REAL    :: K,speed            
  REAL    :: error,RMS   ! root mean square error between current and initial sol.

  current_time      = 0.
  snap_current_time = 0.

  CALL initialize_solution(C,N,DX,first_limit,second_limit)         
  CEXACT = C
  current_first_limit = first_limit

  CALL pre_loop_initialization(mode,method,current_time,K, &
                               DT,DX,N,C,ALL_DATA_FILE_ID, &
                               A,aa,b,cc )

  snap_shot_at_15_taken=.FALSE.  

  DO WHILE(current_time < max_run_time)               

     IF( snap_current_time >= snapshot_delta ) THEN
        snap_current_time = 0.       
        WRITE( UNIT=ALL_DATA_FILE_ID, FMT=*) current_time, error, C
     END IF

     SELECT CASE(method)

     CASE( 1:2 )

        IF(method==1) THEN  ! ftcs           
           IF(mode==2)THEN
              K  = speed(mode,current_time)*DT/(2.*DX)                
           ENDIF

           DO I = 2,N-1
              CNEW(I) = C(I) - K * ( C(I+1) - C(I-1) )  
           END DO
        ELSE    !lax
           IF(mode == 2) THEN   
              K = speed(mode,current_time)*DT/(DX)    
           ENDIF

           DO I = 2,N-1
              CNEW(I) = C(I) - K/2. * ( C(I+1) - C(I-1) ) + &
                   (K**2.)/2 * ( C(I+1) +C(I-1)-2.*C(I) ) 
           END DO
        END IF

        CNEW(1) = C(1)
        CNEW(N) = C(N)  ! Boundary conditions
        C=CNEW

     CASE( 3 )  ! implicit ftcs

        IF( mode == 2) THEN ! only need to update Matrix for varying U
           K  = speed(mode,current_time)*DT/(2.*DX) 

           CALL init_A_matrix(A,K,N)
           CALL init_diagonal_vectors(N,A,cc,aa,b)
        END IF

        CALL solve_thomas_algorithm(N,aa,b,cc,C,CNEW)
        C = CNEW

     CASE( 4 )   ! C-R

        IF(mode == 2) THEN  !only need to update A if U changes                       
           K  = speed(mode,current_time)*DT/(4*DX)   ! C-R            
           CALL init_A_matrix(A,K,N)
           CALL init_diagonal_vectors(N,A,cc,aa,b)
        END IF

        CTEMP(1) = C(1)
        CTEMP(N) = C(N)

        DO I=2,N-1
           CTEMP(I)=C(I)+K*C(I-1)-K*C(I+1)
        END DO

        CALL solve_thomas_algorithm(N,aa,b,cc,CTEMP,C)

     END SELECT

     IF( current_time>=15.0 .AND. (.NOT. snap_shot_at_15_taken)) THEN
        snap_shot_at_15_taken = .TRUE.
        CALL take_one_snap_shot(mode,method,15,N,C,DX)                
     END IF

     current_time = current_time + DT
     current_first_limit = current_first_limit + speed(mode,current_time)*DT    
     CALL get_current_exact_solution(CEXACT,N,current_first_limit,DX)            
     error = RMS(CEXACT,C,N)

     snap_current_time = snap_current_time + DT

  END DO

  CALL take_one_snap_shot(mode,method,30,N,C,DX)   

END SUBROUTINE process
!************************************
!*
!*
!************************************
SUBROUTINE pre_loop_initialization(mode,method,current_time,K,&
                                   DT,DX,N,C,ALL_DATA_FILE_ID,&
                                   A,aa,b,cc)
  IMPLICIT NONE

  INTEGER, INTENT(IN)  :: mode,method,N,ALL_DATA_FILE_ID
  REAL,    INTENT(IN)  :: C(N),DT,DX,current_time
  REAL,    INTENT(OUT) :: K,A(N,N),aa(N),b(2:N),cc(N-1)
  REAL                 :: speed

  SELECT CASE(method)
  CASE( 1 )   ! FTCS

     K  = speed(mode,current_time)*DT/(2.*DX)                  

     IF(mode==1) THEN               
        OPEN(UNIT=ALL_DATA_FILE_ID, file='expAll.txt') ! all time shots     
        CALL print_to_file(C,'exp0.txt',N,DX)             
     ELSE
        OPEN(UNIT=ALL_DATA_FILE_ID, file='exp_extraAll.txt') ! all time shots     
        CALL print_to_file(C,'exp_extra0.txt',N,DX)             
     END IF

  CASE( 2 )    ! Lax

     K = speed(mode,current_time)*DT/(DX)      

     IF(mode==1) THEN
        OPEN(UNIT=ALL_DATA_FILE_ID, file='laxAll.txt') ! all time shots          
        CALL print_to_file(C,'lax0.txt',N,DX)
     ELSE
        OPEN(UNIT=ALL_DATA_FILE_ID, file='lax_extraAll.txt') ! all time shots          
        CALL print_to_file(C,'lax_extra0.txt',N,DX)        
     END IF

  CASE( 3 )    ! Implicit FTCS  

     K  = speed(mode,current_time)*DT/(2.*DX) 

     CALL init_A_matrix(A,K,N)
     CALL init_diagonal_vectors(N,A,cc,aa,b)

     IF(mode==1) THEN
        OPEN(UNIT=ALL_DATA_FILE_ID, file='impAll.txt') ! all time shots               
        CALL print_to_file(C,'imp0.txt',N,DX)
     ELSE
        OPEN(UNIT=ALL_DATA_FILE_ID, file='imp_extraAll.txt') ! all time shots               
        CALL print_to_file(C,'imp_extra0.txt',N,DX)
     END IF

  CASE( 4 )     ! C-R

     K  = speed(mode,current_time)*DT/(4*DX)   ! C-R    

     CALL init_A_matrix(A,K,N)
     CALL init_diagonal_vectors(N,A,cc,aa,b)

     IF(mode==1) THEN
        OPEN(UNIT=ALL_DATA_FILE_ID, file='crAll.txt') ! all time shots               
        CALL print_to_file(C,'cr0.txt',N,DX)
     ELSE
        OPEN(UNIT=ALL_DATA_FILE_ID, file='cr_extraAll.txt') ! all time shots               
        CALL print_to_file(C,'cr_extra0.txt',N,DX)
     END IF
  END SELECT

  WRITE( UNIT=ALL_DATA_FILE_ID, FMT=*) current_time,0, C    

END SUBROUTINE pre_loop_initialization
!************************************
!*
!*
!************************************
SUBROUTINE init_diagonal_vectors(N,A,cc,aa,b)
  IMPLICIT NONE

  INTEGER, INTENT(IN) ::N                        
  REAL, INTENT(IN)    ::A(N,N)
  REAL, INTENT(OUT)   ::aa(N),b(2:N),cc(N-1)

  INTEGER ::I,J

  J=2                         
  DO I=1,N-1
     cc(I)=A(I,J)
     J=J+1
  END DO
  cc(1)=0    

  DO I=1,N
     aa(I)=A(I,I)
  END DO

  J=1
  DO I=2,N
     b(I)=A(I,J)
     J=J+1
  END DO

END SUBROUTINE init_diagonal_vectors
!************************************
!*
!*
!************************************
SUBROUTINE initialize_solution(C,N,DX,first_limit,second_limit)
  IMPLICIT NONE

  INTEGER, INTENT(IN)    :: N
  REAL,    INTENT(IN)    :: DX,first_limit,second_limit
  REAL,    INTENT(INOUT) :: C(0:N-1)

  INTEGER :: I
  REAL    :: x, PI,av,R   

  PARAMETER( PI = ACOS(-1.) )

  x = 0
  av = (second_limit+first_limit)/2.0
  R  = av - first_limit

  C = 0.0

  DO I=0,N-1

     IF( x >= first_limit .AND. x <= second_limit ) THEN
        C(I) = 1 + COS( PI * (x-av)/R )
     END IF

     x = x + DX
  END DO

END SUBROUTINE initialize_solution
!************************************
!*
!*
!************************************
SUBROUTINE print_to_file(C,file_name,N,DX)
  IMPLICIT NONE


  REAL,    INTENT(IN) :: C(N),DX
  INTEGER, INTENT(IN) :: N

  CHARACTER* (*), INTENT(IN) :: file_name

  INTEGER :: I
  INTEGER :: FILE_ID
  PARAMETER(FILE_ID=999)
  REAL :: current_position

  OPEN(UNIT=FILE_ID, file=file_name)

  current_position = 0;
  DO I=1,N

     WRITE( UNIT=FILE_ID, FMT=* ) current_position ,'\t', C(I)
     current_position = current_position + DX

  END DO

  CLOSE(FILE_ID)

END SUBROUTINE print_to_file
!************************************
!*
!*
!************************************
SUBROUTINE init_A_matrix(A,K,N)
  IMPLICIT NONE

  INTEGER, INTENT(IN) ::N                            
  REAL, INTENT(IN)    ::K
  REAL, INTENT(OUT)   ::A(N,N)

  INTEGER ::I 

  DO I = 2,N-1
     A(I,I-1)   = -K
     A(I,I) = 1
     A(I,I+1) = K
  END DO

  A(1,1) = 1
  A(N,N) = 1

END SUBROUTINE init_A_matrix
!************************************
!*
!*
!************************************                
SUBROUTINE solve_thomas_algorithm(N,aa,b,c,old_c,new_c)
  IMPLICIT NONE

  REAL, INTENT(IN)    :: aa(N),b(2:N),c(N-1),old_c(N)
  INTEGER, INTENT(IN) :: N
  REAL, INTENT(INOUT) :: new_c(N)

  INTEGER  :: I 
  REAL :: alpha(N),beta(2:N),g(N)  

  alpha(1) = aa(1)
  DO I=2,N
     beta(I)=b(I)/alpha(I-1)
     alpha(I)=aa(I)-beta(I)*c(I-1)
  END DO

  g(1)=old_c(1)
  DO I=2,N
     g(I)=old_c(I)-beta(I)*g(I-1)
  END DO

  new_c(N)=g(N)/alpha(N)
  DO I=N-1,1,-1
     new_c(I)=(g(I)-c(I)*new_c(I+1))/alpha(I)
  END DO

END SUBROUTINE solve_thomas_algorithm
!************************************
!*
!*
!************************************                
REAL FUNCTION speed(MODE,time)
  IMPLICIT NONE

  INTEGER, INTENT(IN) :: MODE
  REAL,    INTENT(IN) :: time

  IF( MODE == 1 ) THEN
     speed=2.0
  ELSE
     speed=time/20.0
  END IF

END FUNCTION speed
!************************************
!*
!*
!************************************
SUBROUTINE take_one_snap_shot(mode,method,TIME,N,C,DX)
  IMPLICIT NONE

  INTEGER, INTENT(IN) ::TIME,mode,method,N
  REAL,    INTENT(IN) ::C(N),DX

  IF(TIME==15) THEN
     SELECT CASE(method)
     CASE(1)
        IF(mode==1) THEN
           CALL print_to_file(C,'exp15.txt',N,DX)
        ELSE
           CALL print_to_file(C,'exp_extra15.txt',N,DX)
        END IF
     CASE(2)
        IF(mode==1) THEN
           CALL print_to_file(C,'lax15.txt',N,DX)
        ELSE
           CALL print_to_file(C,'lax_extra15.txt',N,DX)
        ENDIF
     CASE(3)
        IF(mode==1) THEN
           CALL print_to_file(C,'imp15.txt',N,DX)          
        ELSE
           CALL print_to_file(C,'imp_extra15.txt',N,DX)          
        END IF
     CASE(4)
        IF(mode==1) THEN
           CALL print_to_file(C,'cr15.txt',N,DX)          
        ELSE
           CALL print_to_file(C,'cr_extra15.txt',N,DX)          
        END IF
     END SELECT
  ELSE
     SELECT CASE(method)
     CASE(1)
        IF(mode==1) THEN
           CALL print_to_file(C,'exp30.txt',N,DX)
        ELSE
           CALL print_to_file(C,'exp_extra30.txt',N,DX)
        END IF
     CASE(2)
        IF(mode==1) THEN
           CALL print_to_file(C,'lax30.txt',N,DX)
        ELSE
           CALL print_to_file(C,'lax_extra30.txt',N,DX)
        ENDIF
     CASE(3)
        IF(mode==1) THEN
           CALL print_to_file(C,'imp30.txt',N,DX)          
        ELSE
           CALL print_to_file(C,'imp_extra30.txt',N,DX)          
        END IF
     CASE(4)
        IF(mode==1) THEN
           CALL print_to_file(C,'cr30.txt',N,DX)          
        ELSE
           CALL print_to_file(C,'cr_extra30.txt',N,DX)          
        END IF
     END SELECT
  END IF
END SUBROUTINE    take_one_snap_shot
!************************************
!*
!*
!************************************
REAL FUNCTION RMS(CEXACT,C,N)
  IMPLICIT NONE

  REAL, INTENT(IN)    :: CEXACT(N),C(N)
  INTEGER, INTENT(IN) :: N

  INTEGER  :: I

  RMS=0.
  DO I=1,N
     RMS = RMS+(CEXACT(I)-C(I))**2
  END DO

  RMS = RMS/N                                      
  RMS = SQRT(RMS)
END FUNCTION RMS
!************************************
!*
!*
!************************************     
SUBROUTINE get_current_exact_solution(CEXACT,N,current_first_limit,DX)       
  IMPLICIT NONE
  REAL, INTENT(IN)    :: current_first_limit,DX
  REAL, INTENT(OUT)   :: CEXACT(0:N-1)
  INTEGER, INTENT(IN) :: N

  INTEGER :: I
  REAL :: first_limit  
  REAL :: second_limit 
  REAL :: av,R,shift,x,PI

  PARAMETER( PI = ACOS(-1.) )

  first_limit  = 25.0
  second_limit = 35.0

  shift = current_first_limit - FIRST_LIMIT
  first_limit = current_first_limit
  second_limit = second_limit + shift

  av = (second_limit+first_limit)/2.0
  R  = av - first_limit

  CEXACT = 0.
  x = 0. 
  DO I = 0,N-1

     IF( x >= first_limit .AND. x <= second_limit ) THEN
        CEXACT(I) = 1 + COS( PI * (x -av)/R )
     END IF

     x = x + DX
  END DO
END SUBROUTINE get_current_exact_solution