PROGRAM main

  USE var_dep
  IMPLICIT NONE
  
  COMMON/count_alloc/nalloc
  INTEGER :: nalloc
  
  nalloc=0
  
  CALL vardep_test(1000,4,1000)
   
  CONTAINS 
   
  SUBROUTINE vardep_test(nx,ng,niter)
  
   INTEGER,INTENT(in) :: nx    ! number of meshes
   INTEGER,INTENT(in) :: ng    ! number of ideal gases
   INTEGER,INTENT(in) :: niter ! number of iterations to increase the CPU
   
   TYPE(vardep),ALLOCATABLE :: PP(:,:),P(:),T(:),n(:),rho(:),v(:),q(:)
   REAL(r8)                 :: M=18.e-3_r8,Vol,R=8.314_r8,S=2_r8,fric=30_r8
   REAL(r8)    ,ALLOCATABLE :: sp(:),st(:),sn(:),srho(:),sv(:),sq(:)
   INTEGER                  :: i,k
   REAL                     :: t0,t1,t2,t3,t4,t5
   INTEGER                  :: n0,n1,n2,n3,n4,n5
   
   Vol=50._r8/nx ! volume of a mesh
   
   ALLOCATE(pp(ng,nx),p(nx),t(nx),n(nx),rho(nx),v(nx-1),q(nx-1))
   ALLOCATE(sp(nx),st(nx),sn(nx),srho(nx),sv(nx-1),sq(nx-1))
   
   ! Initializing main variables (pressures and temperatures)
   
   P=0
   DO i=1,nx
      DO k=1,ng
        CALL newMainVar(PP(k,i),1.e5_r8*(5-DBLE(i+k)/nx),(ng+1)*(i-1)+k)
        P(i)=P(i)+PP(k,i) ! total pressure
      ENDDO
      CALL newMainVar(T(i),1000._r8*(1+DBLE(i)/nx),(ng+1)*i)
      sp(i)=p(i)%value
      st(i)=t(i)%value
   ENDDO
   
   ! Allocating main results only once (for the ALLOCATABLE version)
   
   n=(P/T)*(Vol/R)                   ! number of moles
   rho=n*(M/Vol)                     ! density in assuming an average molar mass
   v=(P(1:nx-1)-P(2:nx))/fric        ! velocity
   q=v*(rho(1:nx-1)+rho(2:nx))*(S/2) ! mass flow rate
   
   ! Use of subroutines
   
   CALL cpu_time(t0) ; n0=nalloc
   DO k=1,niter
   DO i=1,nx
      CALL newVar(n(i),P(i))
      CALL divVar(n(i),T(i))
      CALL multVar(n(i),Vol/R)
      CALL newVar(rho(i),n(i))
      CALL multVar(rho(i),M/Vol)
      IF(i < nx) THEN 
        CALL newVar  (v(i),P(i))
        CALL minusVar(v(i),P(i+1))
        CALL divVar  (v(i),fric)
        CALL newVar  (q(i),rho(i))
        CALL addVar  (q(i),rho(i+1))
        CALL multVar (q(i),v(i))
        CALL multVar (q(i),S/2)
      ENDIF
   ENDDO
   ENDDO
   write(*,10) 'NR ',n(nx)%value,rho(nx)%value,v(1)%value,q(1)%value
   CALL cpu_time(t1) ; n1=nalloc
   
   ! Use of elemental subroutines
   
   DO k=1,niter
      CALL newVar(n,P)
      CALL divVar(n,T)
      CALL multVar(n,Vol/R)
      CALL newVar(rho,n)
      CALL multVar(rho,M/Vol)
      CALL newVar  (v,P(1:nx-1))
      CALL minusVar(v,P(2:nx))
      CALL divVar  (v,fric)
      CALL newVar  (q,rho(1:nx-1))
      CALL addVar  (q,rho(2:nx))
      CALL multVar (q,v)
      CALL multVar (q,S/2)
    ENDDO
    write(*,10) 'ER ',n(nx)%value,rho(nx)%value,v(1)%value,q(1)%value
   
    ! Use of functions (operators)
    
    CALL cpu_time(t2) ; n2=nalloc
    DO k=1,niter
    DO i=1,nx
       n(i)=(P(i)/T(i))*(Vol/R)
       rho(i)=n(i)*(M/Vol)
       IF(i < nx) THEN
         v(i)=(P(i)-P(i+1))/fric
         q(i)=v(i)*(rho(i)+rho(i+1))*(S/2)
       ENDIF  
    ENDDO
    ENDDO
    write(*,10) 'NF ',n(nx)%value,rho(nx)%value,v(1)%value,q(1)%value
   
    ! Use of elemental functions (vector operators)
    
    CALL cpu_time(t3) ; n3=nalloc
    DO k=1,niter
       n=(P/T)*(Vol/R)
       rho=n*(M/Vol)
       v=(P(1:nx-1)-P(2:nx))/fric
       q=v*(rho(1:nx-1)+rho(2:nx))*(S/2)
    ENDDO
    write(*,10) 'EF ',n(nx)%value,rho(nx)%value,v(1)%value,q(1)%value
   
    ! Scalar compupation
    
    CALL cpu_time(t4) ; n4=nalloc
    DO k=1,niter
       sn=(sp/st)*(Vol/R)
       srho=sn*(M/Vol)
       sv=(sp(1:nx-1)-sp(2:nx))/fric
       sq=sv*(srho(1:nx-1)+srho(2:nx))*(S/2)
    ENDDO
    write(*,10) 'SC ',sn(nx),srho(nx),sv(1),sq(1)
    CALL cpu_time(t5) ; n5=nalloc
    
    WRITE(*,*) 'NR ',t1-t0,n1-n0
    WRITE(*,*) 'ER ',t2-t1,n2-n1
    WRITE(*,*) 'NF ',t3-t2,n3-n2
    WRITE(*,*) 'EF ',t4-t3,n4-n3
    WRITE(*,*) 'SC ',t5-t4,n5-n4
    
 10 FORMAT(1x,a,10(1x,1pe11.4))
 
  END SUBROUTINE

END PROGRAM

SUBROUTINE count
  COMMON/count_alloc/nalloc
  INTEGER :: nalloc
  nalloc=nalloc+1
END SUBROUTINE
