module active_module
        use w2f__types
        implicit none
        private
        public :: deriv, active, saxpy, sax, setderiv, zero_deriv, &
                & getderiv, init_adjoint, set_adjoint, get_adjoint, &
                & adjoint_interpreter, dump_tape, dump_adjoint, &
                & init_fm_rm, finish

        integer, parameter :: tape_chunk_size=100000
        integer :: tape_counter=1
        integer, parameter :: addresses_chunk_size=100
        integer :: addresses_counter=0
        integer :: nr_adjoints
! need a static here if active common blocks present like in box_timestep
        integer, parameter :: nr_tangents=6

        type deriv
          sequence
          double precision, dimension(nr_tangents) :: d
        end type deriv

        type adj
          sequence
          double precision, dimension(:), pointer :: d=>NULL()
        end type adj
        !
        ! active needs to be a sequence type
        !  with no initialization
        !
        type active
          sequence
          double precision :: v 
          type(deriv) :: d
        end type active

        !
        ! x=x+a*y
        ! 
        type adjoint_saxpy
          integer :: type ! 0=zero_deriv, 1=setderiv, 2=saxpy, 3=sax
          integer :: addr_x
          integer :: addr_y
          double precision :: a
        end type adjoint_saxpy

        !
        ! can make this allocatable etc.
        !
        type(adjoint_saxpy), dimension(:), allocatable :: tape
        integer, dimension(:), allocatable :: addresses
        type(adj), dimension(:), allocatable :: adjoint

        interface init_fm_rm
          module procedure init_fm_rm
        end interface

        interface finish
          module procedure finish
        end interface

        interface saxpy
          module procedure saxpy_a_a
        end interface
        
        interface setderiv
          module procedure setderiv_a_a, setderiv_a_c
        end interface

        interface zero_deriv
          module procedure zero_deriv_a
        end interface
        
        interface sax
          module procedure sax_d_a_a, sax_i_a_a
        end interface

        interface getderiv
          module procedure getderiv
        end interface

        interface set_adjoint
          module procedure set_adjoint_scal, set_adjoint_vec
        end interface

        interface get_adjoint
          module procedure get_adjoint
        end interface

        interface adjoint_interpreter
          module procedure adjoint_interpreter
        end interface

        interface dump_tape
          module procedure dump_tape
        end interface

        interface dump_adjoint
          module procedure dump_adjoint
        end interface

        contains

        subroutine init_fm_rm(nr_tan, nr_adj)
          integer, intent(in) :: nr_tan, nr_adj
    !      nr_tangents=nr_tan
          nr_adjoints=nr_adj
        end subroutine init_fm_rm

        subroutine finish
          deallocate(tape)
          deallocate(addresses)
          deallocate(adjoint)
        end subroutine finish

        subroutine check_addresses_size
          implicit none
          integer, dimension(:), allocatable :: addresses_tmp
          integer s
          if (.not.allocated(addresses)) then
            allocate(addresses(addresses_chunk_size))
          endif
          s=size(addresses)
          if (addresses_counter>=s) then
            allocate(addresses_tmp(s))
            addresses_tmp=addresses 
            deallocate(addresses)
            allocate(addresses(s+addresses_chunk_size))
            addresses=addresses_tmp
            deallocate(addresses_tmp)
          endif
        end subroutine check_addresses_size

        integer function hash(x)
          integer, intent(in) :: x
          integer i
          call check_addresses_size
          do i=1,addresses_counter
            if (addresses(i)==x) then
              hash=i
              return
            endif
          end do
          addresses_counter=addresses_counter+1
          addresses(addresses_counter)=x 
          hash=addresses_counter
          return
        end function hash

!
! Local routine to handle dynamic allocation of tape
!

        subroutine check_tape_size
          type(adjoint_saxpy), dimension(:), allocatable :: tape_tmp
          integer s
          if (.not.allocated(tape)) then
            allocate(tape(tape_chunk_size))
          endif
          s=size(tape)
          if (tape_counter>s) then
            allocate(tape_tmp(s))
            tape_tmp=tape 
            deallocate(tape)
            allocate(tape(s+tape_chunk_size))
            tape(1:s)=tape_tmp
            deallocate(tape_tmp)
          endif
            
        end subroutine check_tape_size
        
        !
        ! chain rule saxpy to be used in forward and reverse modes
        !
        
        subroutine saxpy_a_a(a,x,y)
          double precision, intent(in) :: a
          type(active), intent(in) :: x
          type(active), intent(inout) :: y
          type(adjoint_saxpy) :: adj_saxpy
          integer iaddr

          call check_tape_size
        
          adj_saxpy%type=2
          adj_saxpy%addr_x=hash(iaddr(x))
          adj_saxpy%addr_y=hash(iaddr(y))
          adj_saxpy%a=a
          tape(tape_counter)=adj_saxpy
          tape_counter=tape_counter+1
         
          y%d%d=y%d%d+x%d%d*a
        end subroutine saxpy_a_a
        
        !
        ! chain rule saxpy to be used in forward and reverse modes
        ! derivative component of y is equal to zero initially
        ! note: y needs to be inout as otherwise value component gets
        ! zeroed out
        !
        
        subroutine sax_d_a_a(a,x,y)
          double precision, intent(in) :: a
          type(active), intent(in) :: x
          type(active), intent(inout) :: y
          type(adjoint_saxpy) :: adj_saxpy
          integer iaddr
        
          call check_tape_size

          adj_saxpy%type=3
          adj_saxpy%addr_x=hash(iaddr(x))
          adj_saxpy%addr_y=hash(iaddr(y))
          adj_saxpy%a=a
          tape(tape_counter)=adj_saxpy
          tape_counter=tape_counter+1
       
          y%d%d=x%d%d*a
        end subroutine sax_d_a_a

        subroutine sax_i_a_a(a,x,y)
          integer(kind=w2f__i8), intent(in) :: a
          type(active), intent(in) :: x
          type(active), intent(inout) :: y
          type(adjoint_saxpy) :: adj_saxpy
          integer iaddr
        
          call check_tape_size

          adj_saxpy%type=3
          adj_saxpy%addr_x=hash(iaddr(x))
          adj_saxpy%addr_y=hash(iaddr(y))
          adj_saxpy%a=a
          tape(tape_counter)=adj_saxpy
          tape_counter=tape_counter+1

          y%d%d=x%d%d*a
        end subroutine sax_i_a_a
        
        !
        ! set derivative of y to be equal to derivative of x
        ! note: making y inout allows for already existing active
        ! variables to become the target of a derivative assignment
        !
        
        subroutine setderiv_a_a(y,x)
          type(active), intent(inout) :: y
          type(active), intent(in) :: x
          type(adjoint_saxpy) :: adj_saxpy
          integer iaddr
        
          call check_tape_size

          adj_saxpy%type=1
          adj_saxpy%addr_x=hash(iaddr(x))
          adj_saxpy%addr_y=hash(iaddr(y))
          adj_saxpy%a=1.D0
          tape(tape_counter)=adj_saxpy
          tape_counter=tape_counter+1

          y%d%d=x%d%d
        end subroutine setderiv_a_a

        !
        ! initialization of derivative component with a constant
        !

        subroutine setderiv_a_c(y,x)
          type(active), intent(inout) :: y
          double precision, dimension(:), intent(in) :: x
          y%d%d=x
        end subroutine setderiv_a_c

        !
        ! set derivative components to 0.0
        !
        subroutine zero_deriv_a(x)
          type(active) :: x
          type(adjoint_saxpy) :: adj_saxpy
          integer iaddr
        
          call check_tape_size

          adj_saxpy%type=0
          adj_saxpy%addr_x=hash(iaddr(x))
          adj_saxpy%addr_y=0
          adj_saxpy%a=0.D0
          tape(tape_counter)=adj_saxpy
          tape_counter=tape_counter+1

          x%d%d=0.0d0
        end subroutine zero_deriv_a

        !
        ! return double value of deriv
        !

        function getderiv(x) result(res)
          type(active), intent(in) :: x
          double precision, dimension(nr_tangents) :: res
          res=x%d%d
          return
        end function getderiv

        
        subroutine init_adjoint
          integer i
          allocate(adjoint(addresses_counter)) 
          do i=1,addresses_counter
            allocate(adjoint(i)%d(nr_adjoints))
            adjoint(i)%d=0.D0
          end do
        end subroutine init_adjoint

        !
        ! initialization of adjoint component with a constant
        !

        subroutine set_adjoint_vec(y,x)
          type(active), intent(inout) :: y
          double precision, dimension(nr_adjoints), &
            & intent(in) :: x
          integer iaddr
          adjoint(hash(iaddr(y)))%d=x
        end subroutine set_adjoint_vec

        subroutine set_adjoint_scal(y,x)
          type(active), intent(inout) :: y
          double precision, intent(in) :: x
          integer iaddr
          adjoint(hash(iaddr(y)))%d=x
        end subroutine set_adjoint_scal

        function get_adjoint(y) result(res)
          type(active), intent(in) :: y
          double precision, dimension(nr_adjoints) :: res
          integer iaddr
          res=adjoint(hash(iaddr(y)))%d
          return
        end function get_adjoint


        subroutine adjoint_interpreter
          integer i
          do i=tape_counter-1,1,-1
            select case (tape(i)%type) 
             case (0) 
               adjoint(tape(i)%addr_x)%d=0.D0
             case (1) 
               adjoint(tape(i)%addr_x)%d=adjoint(tape(i)%addr_x)%d &
                  & +adjoint(tape(i)%addr_y)%d
               adjoint(tape(i)%addr_y)%d=0.D0
             case (2) 
               adjoint(tape(i)%addr_x)%d= adjoint(tape(i)%addr_x)%d &
                  & +adjoint(tape(i)%addr_y)%d* &
                  & tape(i)%a
             case (3) 
               adjoint(tape(i)%addr_x)%d= adjoint(tape(i)%addr_x)%d &
                  & +adjoint(tape(i)%addr_y)%d &
                  & *tape(i)%a
               adjoint(tape(i)%addr_y)%d=0.D0
            end select
          end do
        
        end subroutine adjoint_interpreter

        subroutine dump_tape
          integer i
          print*, "tape"
          do i=1,tape_counter-1
            print*, tape(i)
          end do 
        end subroutine dump_tape

        subroutine dump_adjoint
          integer i
          print*, "adjoint"
          do i=1,addresses_counter-1
            print*, adjoint(i)%d
          end do 
        end subroutine dump_adjoint

        end module