MODULE var_dep

   ! management of general variables with dependencies (scalar value + derivatives)
   ! in computing formal derivatives for usual operations.
   
   ! From a technical point of view, derivative indexes are always ordered
   
   ! the derivatives can be managed as allocatable arrays but, unfortunately,
   ! the performances are then strongly degraded : it seems better to fix the maximum
   ! number of derivatives (parameter ndmax). 
   
   ! It is possible to activate the version with allocatable arrays : you just have
   ! to uncomment lines starting by !! and delete few useless lines ending by !!
   ! Two entities are concerned : the vardep derived type and the routine AllocDeriv
   
   IMPLICIT NONE

   INTEGER,PARAMETER :: ndmax=19,r8=SELECTED_REAL_KIND(15,100) ! may be changed if necessary
   
   PRIVATE
   
   PUBLIC vardep,r8 ! the derived type
   PUBLIC assignment(=)
   PUBLIC operator(+),operator(-),operator(*),operator(/),operator(**)                       ! classical operator
   PUBLIC addvar,multvar,minusvar,divvar,powvar                                              ! idem but subroutines : v1 <- v1 ? v2
   PUBLIC operator(==),operator(<=),operator(>=),operator(>),operator(/=)                    ! logical operators
   PUBLIC sin,cos,tan,tanh,atan,exp,log,sqrt                                                 ! classical functions
   PUBLIC sinus,cosinus,tangent,tangenthyperbolic,arctangent,exponent,logarithm,squaredroot  ! idem but as subroutines
   PUBLIC min,max,abs
   PUBLIC newVar,newMainVar,inverseVar,setDeriv                                              ! useful other functions
   PUBLIC numericalDeriv                                                                     ! useful other functions
   
   TYPE vardep
      REAL(r8) :: value
       !!REAL(r8) :: deriv(ndmax) !!
      !!INTEGER  :: index(ndmax) !!
      REAL(r8),ALLOCATABLE :: deriv(:)   ! may be activated instead of the two previous lines
      INTEGER ,ALLOCATABLE :: index(:)   ! the routine AllocDeriv must be modified too
      INTEGER  :: nderiv=-1
   END TYPE
   
   INTERFACE assignment(=)     ; MODULE PROCEDURE newVar1,newVarReal,newVarInt                                  ; END INTERFACE
   INTERFACE operator(+)       ; MODULE PROCEDURE VarPlusVar,VarPlusReal,VarPlusInt,RealPlusVar,IntPlusVar      ; END INTERFACE
   INTERFACE operator(-)       ; MODULE PROCEDURE VarMinusVar,VarMinusReal,VarMinusInt,RealMinusVar,IntMinusVar ; END INTERFACE
   INTERFACE operator(*)       ; MODULE PROCEDURE VarMultVar,VarMultReal,VarMultInt,RealMultVar,IntMultVar      ; END INTERFACE
   INTERFACE operator(/)       ; MODULE PROCEDURE VarDivVar,VarDivReal,VarDivInt,RealDivVar,IntDivVar           ; END INTERFACE
   INTERFACE operator(**)      ; MODULE PROCEDURE VarPowVar,VarPowReal,VarPowInt,RealPowVar,IntPowVar           ; END INTERFACE
   INTERFACE operator(==)      ; MODULE PROCEDURE VarEqVar,VarEqReal,VarEqInt,RealEqVar,IntEqVar                ; END INTERFACE
   INTERFACE operator(/=)      ; MODULE PROCEDURE VarNeVar,VarNeReal,VarNeInt,RealNeVar,IntNeVar                ; END INTERFACE
   INTERFACE operator(<)       ; MODULE PROCEDURE VarLtVar,VarLtReal,VarLtInt,RealLtVar,IntLtVar                ; END INTERFACE
   INTERFACE operator(>)       ; MODULE PROCEDURE VarGtVar,VarGtReal,VarGtInt,RealGtVar,IntGtVar                ; END INTERFACE
   INTERFACE operator(<=)      ; MODULE PROCEDURE VarLeVar,VarLeReal,VarLeInt,RealLeVar,IntLeVar                ; END INTERFACE
   INTERFACE operator(>=)      ; MODULE PROCEDURE VarGeVar,VarGeReal,VarGeInt,RealGeVar,IntGeVar                ; END INTERFACE
   INTERFACE newVar            ; MODULE PROCEDURE newVar0,newVar1,newVar2,newVar3,newVarReal,newVarInt          ; END INTERFACE
   INTERFACE setDeriv          ; MODULE PROCEDURE setDeriv0,setPartialDeriv                                     ; END INTERFACE
   INTERFACE addVar            ; MODULE PROCEDURE addVars,addVarReal,addVarInt                                  ; END INTERFACE
   INTERFACE multVar           ; MODULE PROCEDURE multVars,multVarReal,multVarInt                               ; END INTERFACE
   INTERFACE divVar            ; MODULE PROCEDURE divVars,divVarReal,divVarInt                                  ; END INTERFACE
   INTERFACE minusVar          ; MODULE PROCEDURE minusVars,minusVarReal,minusVarInt                            ; END INTERFACE
   INTERFACE powVar            ; MODULE PROCEDURE powVars,powVarReal,powVarInt                                  ; END INTERFACE
   INTERFACE min               ; MODULE PROCEDURE minVarVar,minRealVar,minIntVar,minVarReal,minVarInt           ; END INTERFACE
   INTERFACE max               ; MODULE PROCEDURE maxVarVar,maxRealVar,maxIntVar,maxVarReal,maxVarInt           ; END INTERFACE
   INTERFACE abs               ; MODULE PROCEDURE absFun                                                        ; END INTERFACE
   INTERFACE absolutevalue     ; MODULE PROCEDURE absSub                                                        ; END INTERFACE
   INTERFACE sqrt              ; MODULE PROCEDURE sqrtFun                                                       ; END INTERFACE
   INTERFACE squaredroot       ; MODULE PROCEDURE sqrtSub                                                       ; END INTERFACE
   INTERFACE sin               ; MODULE PROCEDURE sinFun                                                        ; END INTERFACE
   INTERFACE sinus             ; MODULE PROCEDURE sinSub                                                        ; END INTERFACE
   INTERFACE cos               ; MODULE PROCEDURE cosFun                                                        ; END INTERFACE
   INTERFACE cosinus           ; MODULE PROCEDURE cosSub                                                        ; END INTERFACE
   INTERFACE tan               ; MODULE PROCEDURE tanFun                                                        ; END INTERFACE
   INTERFACE tangent           ; MODULE PROCEDURE tanSub                                                        ; END INTERFACE
   INTERFACE tanh              ; MODULE PROCEDURE tanhFun                                                       ; END INTERFACE
   INTERFACE tangenthyperbolic ; MODULE PROCEDURE tanhSub                                                       ; END INTERFACE
   INTERFACE atan              ; MODULE PROCEDURE atanFun                                                       ; END INTERFACE
   INTERFACE arctangent        ; MODULE PROCEDURE atanSub                                                       ; END INTERFACE
   INTERFACE exp               ; MODULE PROCEDURE expFun                                                        ; END INTERFACE
   INTERFACE exponent          ; MODULE PROCEDURE expSub                                                        ; END INTERFACE
   INTERFACE log               ; MODULE PROCEDURE logFun                                                        ; END INTERFACE
   INTERFACE logarithm         ; MODULE PROCEDURE logSub                                                        ; END INTERFACE
   INTERFACE numericalDeriv    ; MODULE PROCEDURE numericalDeriv1,numericalDerivn                               ; END INTERFACE
   
   INTERFACE 
      PURE SUBROUTINE count
      END SUBROUTINE
   END INTERFACE
        
   CONTAINS
   
   SUBROUTINE numericalDeriv1(proc,x,dx,y)
     INTERFACE 
       SUBROUTINE proc(x,y)
         IMPORT r8
         REAL(r8),INTENT(in)  ::  x(:) ! main variables
         REAL(r8),INTENT(inout) ::  y    ! output value
       END SUBROUTINE
     END INTERFACE
     TYPE(vardep),INTENT(in)  ::  x(:) ! main variables
     REAL(r8)    ,INTENT(in)  :: dx(:) ! increment of state variables
     TYPE(vardep),INTENT(inout) ::    y  ! output value
     INTEGER  :: j,n
     REAL(r8) :: xx(SIZE(x)),yy,a(SIZE(x)),z
     n=SIZE(x)
     xx=x(:)%value     
     CALL proc(xx,yy) ! value of the function
     DO j=1,n
       z=xx(j)
       xx(j)=xx(j)+dx(j)
       CALL proc(xx,a(j))
       xx(j)=z
       a(j)=(a(j)-yy)/dx(j)
     ENDDO
     y=yy
     DO j=1,n
       CALL setPartialDeriv(y,a(j),x(j))
     ENDDO
   END SUBROUTINE
   
   SUBROUTINE numericalDerivN(proc,x,dx,y)
     INTERFACE 
       SUBROUTINE proc(x,y)
         IMPORT r8
         REAL(r8),INTENT(in)  ::  x(:) ! main variables
         REAL(r8),INTENT(inout) ::  y(:) ! output values
       END SUBROUTINE
     END INTERFACE
     TYPE(vardep),INTENT(in)  ::  x(:) ! main variables
     REAL(r8)    ,INTENT(in)  :: dx(:) ! increment of state variables
     TYPE(vardep),INTENT(inout) ::  y(:) ! output values
     INTEGER :: i,j,nx,ny
     REAL(r8) :: xx(SIZE(x)),yy(SIZE(y)),a(SIZE(y),SIZE(x)),z
     nx=SIZE(x)
     ny=SIZE(y)
     xx=x(:)%value     
     CALL proc(xx,yy) ! value of the function
     DO j=1,nx
       z=xx(j)
       xx(j)=xx(j)+dx(j)
       CALL proc(xx,a(:,j))
       xx(j)=z
       a(:,j)=(a(:,j)-yy)/dx(j)
     ENDDO
     y=yy
     DO i=1,ny
       y(i)=yy(i)
       DO j=1,nx
         CALL setPartialDeriv(y(i),a(i,j),x(j))
       ENDDO
     ENDDO
   END SUBROUTINE
  
   PURE ELEMENTAL SUBROUTINE addVars(v1,v2)  ! v1 <- v1+v2
      TYPE(vardep),INTENT(inout) :: v1
      TYPE(vardep),INTENT(in)    :: v2
      INTEGER   :: n,i,j
      TYPE(vardep) :: r
      CALL newVar(r,v1,v2)
      n=r%nderiv
      r%value=v1%value+v2%value
      j=0
      DO i=1,v1%nderiv
         DO j=j+1,n
            IF(v1%index(i) == r%index(j)) THEN
               r%deriv(j)=v1%deriv(i)
               EXIT
            ENDIF
         ENDDO
      ENDDO
      j=0
      DO i=1,v2%nderiv
         DO j=j+1,n
            IF(v2%index(i) == r%index(j)) THEN
               r%deriv(j)=r%deriv(j)+v2%deriv(i)
               EXIT
            ENDIF
         ENDDO
      ENDDO
      CALL newVar(v1,r)
      CALL deleteVar(r)
   END SUBROUTINE
   
   PURE ELEMENTAL SUBROUTINE addVarReal(v,c)  ! v <- v+c
      TYPE(vardep) ,INTENT(inout) :: v
      REAL(r8)     ,INTENT(in)    :: c
      v%value=v%value+c
   END SUBROUTINE
   
   PURE ELEMENTAL SUBROUTINE addVarInt(v,c)  ! v <- v+c
      TYPE(vardep) ,INTENT(inout) :: v
      INTEGER      ,INTENT(in)    :: c
      v%value=v%value+c
   END SUBROUTINE
   
   PURE ELEMENTAL SUBROUTINE minusVars(v1,v2)  ! v1 <- v1-v2
      TYPE(vardep),INTENT(inout) :: v1
      TYPE(vardep),INTENT(in)    :: v2
      TYPE(vardep) :: r
      INTEGER :: n,i,j
      CALL newVar(r,v1,v2)
      n=r%nderiv
      r%deriv(1:n)=0._r8
      r%value=v1%value-v2%value
      j=0
      DO i=1,v1%nderiv
         DO j=j+1,n
            IF(v1%index(i) == r%index(j)) THEN
               r%deriv(j)=v1%deriv(i)
               EXIT
            ENDIF
         ENDDO
      ENDDO
      j=0
      DO i=1,v2%nderiv
         DO j=j+1,n
            IF(v2%index(i) == r%index(j)) THEN
               r%deriv(j)=r%deriv(j)-v2%deriv(i)
               EXIT
            ENDIF
         ENDDO
      ENDDO
      CALL newVar(v1,r)
      CALL deleteVar(r)
   END SUBROUTINE
   
   PURE ELEMENTAL SUBROUTINE minusVarReal(v,c) ! v <- v-c
      TYPE(vardep),INTENT(inout)    :: v
      REAL(r8)      ,INTENT(in)    :: c
      v%value=v%value-c
   END SUBROUTINE
   
   PURE ELEMENTAL SUBROUTINE minusVarInt(v,c) ! v <- v-c
      TYPE(vardep),INTENT(inout)    :: v
      INTEGER     ,INTENT(in)    :: c
      v%value=v%value-c
   END SUBROUTINE
   
   PURE ELEMENTAL SUBROUTINE multVars(v1,v2) ! v1 <- v1*v2
      TYPE(vardep),INTENT(inout) :: v1
      TYPE(vardep),INTENT(in)    :: v2
      INTEGER :: n,i,j
      TYPE(vardep) :: r
      CALL newVar(r,v1,v2)
      n=r%nderiv
      r%value=v1%value*v2%value
      j=0
      DO i=1,v1%nderiv
         DO j=j+1,n
            IF(v1%index(i) == r%index(j)) THEN
               r%deriv(j)=v1%deriv(i)*v2%value
               EXIT
            ENDIF
         ENDDO
      ENDDO
      j=0
      DO i=1,v2%nderiv
         DO j=j+1,n
            IF(v2%index(i) == r%index(j)) THEN
               r%deriv(j)=r%deriv(j)+v2%deriv(i)*v1%value
               EXIT
            ENDIF
         ENDDO
      ENDDO
      CALL newVar(v1,r)
      CALL deleteVar(r)
   END SUBROUTINE
   
   PURE ELEMENTAL SUBROUTINE multVarReal(v,c) ! v <- v*c
      TYPE(vardep),INTENT(inout) :: v 
      REAL(r8)      ,INTENT(in)    :: c
      INTEGER ::  n
      n=v%nderiv
      v%value     =v%value     *c
      v%deriv(1:n)=v%deriv(1:n)*c
   END SUBROUTINE
   
   PURE ELEMENTAL SUBROUTINE multVarInt(v,c) ! v <- v*c
      TYPE(vardep),INTENT(inout) :: v 
      INTEGER     ,INTENT(in)    :: c
      INTEGER ::  n
      n=v%nderiv
      v%value     =v%value     *c
      v%deriv(1:n)=v%deriv(1:n)*c
   END SUBROUTINE
   
   PURE ELEMENTAL SUBROUTINE divVars(v1,v2) ! v1 <- v1/v2
      TYPE(vardep),INTENT(inout) :: v1
      TYPE(vardep),INTENT(in)    :: v2
      TYPE(vardep)               :: r     ! the resulting variable
      INTEGER :: n,i,j
      CALL newVar(r,v1,v2)
      n=r%nderiv
      r%value=v1%value/v2%value
      j=0
      DO i=1,v1%nderiv
         DO j=j+1,n
            IF(v1%index(i) == r%index(j)) THEN
               r%deriv(j)=v1%deriv(i)/v2%value
               EXIT
            ENDIF
         ENDDO
      ENDDO
      j=0
      DO i=1,v2%nderiv
         DO j=j+1,n
            IF(v2%index(i) == r%index(j)) THEN
               r%deriv(j)=r%deriv(j)-v2%deriv(i)*v1%value/v2%value**2
               EXIT
            ENDIF
         ENDDO
      ENDDO
      CALL newVar(v1,r)
      CALL deleteVar(r)
   END SUBROUTINE
   
   PURE ELEMENTAL SUBROUTINE divVarReal(v,c) ! v <- v/c
      TYPE(vardep),INTENT(inout) :: v
      REAL(r8) ,INTENT(in)    :: c
      INTEGER :: n
      n=v%nderiv
      v%value     =v%value     / c
      v%deriv(1:n)=v%deriv(1:n)/ c
   END SUBROUTINE
   
   PURE ELEMENTAL SUBROUTINE divVarInt(v,c) ! v <- v/c
      TYPE(vardep),INTENT(inout) :: v
      INTEGER     ,INTENT(in)    :: c
      INTEGER :: n
      n=v%nderiv
      v%value     =v%value     / c
      v%deriv(1:n)=v%deriv(1:n)/ c
   END SUBROUTINE
    
   PURE ELEMENTAL SUBROUTINE powVarReal(v,c) ! v <- v**c
      TYPE(vardep),INTENT(inout) :: v
      REAL(r8) ,INTENT(in)    :: c
      INTEGER :: n
      REAL(r8) :: d
      n=v%nderiv
      d=v%value**(c-1)
      v%value     =d*v%value
      d=d*c
      v%deriv(1:n)=v%deriv(1:n)*d
   END SUBROUTINE
    
   PURE ELEMENTAL SUBROUTINE powVarInt(v,c) ! v <- v**c
      TYPE(vardep),INTENT(inout) :: v
      INTEGER  ,INTENT(in)    :: c
      INTEGER :: n
      REAL(r8) :: d
      n=v%nderiv
      d=v%value**(c-1)
      v%value     =d*v%value
      d=d*c
      v%deriv(1:n)=v%deriv(1:n)*d
   END SUBROUTINE
    
   PURE ELEMENTAL SUBROUTINE powVars(v1,v2) ! v1 <- v1**v2
      TYPE(vardep),INTENT(inout) :: v1
      TYPE(vardep),INTENT(in)    :: v2
      INTEGER   :: i,j,n
      REAL(r8)  :: d
      TYPE(vardep) :: r
      CALL newVar(r,v1,v2)
      n=r%nderiv
      r%value=v1%value**v2%value
      j=0
      d=r%value*v2%value/v1%value
      DO i=1,v1%nderiv
         DO j=j+1,n
            IF(v1%index(i) == r%index(j)) THEN
               r%deriv(j)=v1%deriv(i)*d
               EXIT
            ENDIF
         ENDDO
      ENDDO
      j=0
      d=r%value*LOG(v1%value)
      DO i=1,v2%nderiv
         DO j=j+1,n
            IF(v2%index(i) == r%index(j)) THEN
               r%deriv(j)=r%deriv(j)+v2%deriv(i)*d
               EXIT
            ENDIF
         ENDDO
      ENDDO
      CALL newVar(v1,r)
      CALL deleteVar(r)
   END SUBROUTINE
   
   PURE ELEMENTAL FUNCTION VarPowVar(v1,v2) RESULT(r)
      TYPE(vardep),INTENT(in) :: v1,v2
      TYPE(vardep) :: r
      r=v1
      CALL powVar(r,v2)
   END FUNCTION
   
   PURE ELEMENTAL FUNCTION VarPowReal(v1,v2) RESULT(r)
      TYPE(vardep),INTENT(in) :: v1
      REAL(r8) ,INTENT(in) :: v2
      TYPE(vardep) :: r
      r=v1
      CALL powVar(r,v2)
   END FUNCTION
   
   PURE ELEMENTAL FUNCTION VarPowInt(v1,v2) RESULT(r)
      TYPE(vardep),INTENT(in) :: v1
      INTEGER  ,INTENT(in) :: v2
      TYPE(vardep) :: r
      r=v1
      CALL powVar(r,v2)
   END FUNCTION
   
   PURE ELEMENTAL FUNCTION RealPowVar(v1,v2) RESULT(r)
      REAL(r8) ,INTENT(in) :: v1
      TYPE(vardep),INTENT(in) :: v2
      TYPE(vardep) :: r
      r=v1
      CALL powVar(r,v2)
   END FUNCTION
   
   PURE ELEMENTAL FUNCTION IntPowVar(v1,v2) RESULT(r)
      INTEGER     ,INTENT(in) :: v1
      TYPE(vardep),INTENT(in) :: v2
      TYPE(vardep) :: r
      r=v1
      CALL powVar(r,v2)
   END FUNCTION
        
   PURE ELEMENTAL SUBROUTINE newMainVar(r,value,index) ! defining a main variable
      TYPE(vardep),INTENT(inout) :: r
      REAL(r8)      ,INTENT(IN)  :: value
      INTEGER       ,INTENT(IN)  :: index
      CALL allocDeriv(r,1,(/index/))
      r%value=value
      r%deriv(1)=1._r8
   END SUBROUTINE
   
   PURE SUBROUTINE allocDeriv(r,n,index)
      TYPE(vardep),INTENT(INOUT) :: r
      INTEGER,INTENT(in) :: n
      INTEGER,INTENT(in) :: index(*)
      IF(r%nderiv < n) THEN                     ! uncomment these lines to activate
         IF(ALLOCATED(r%deriv)) THEN            ! allocatable arrays in the vardep type
            IF(n > SIZE(r%deriv)) THEN          ! (the vardep type must also be modified)
               DEALLOCATE(r%deriv,r%index)
               ALLOCATE(r%deriv(n),r%index(n))
               CALL count
            ENDIF
         ELSE
            ALLOCATE(r%deriv(n),r%index(n))
            CALL count
         ENDIF
      ENDIF 
      r%nderiv=n
      IF(n > 0) THEN
         r%index(1:n)=index(1:n)
         r%deriv(1:n)=0._r8
      ENDIF
   END SUBROUTINE
   
   PURE ELEMENTAL SUBROUTINE newVar0(r)
      TYPE(vardep),INTENT(INOUT) :: r
      CALL allocDeriv(r,0,(/0/))
      r%value=0._r8
   END SUBROUTINE
      
   PURE ELEMENTAL SUBROUTINE newVar1(r,v)
      TYPE(vardep),INTENT(INOUT) :: r
      TYPE(vardep),INTENT(IN)    :: v
      CALL allocDeriv(r,v%nderiv,v%index)
      r%value=v%value
      r%deriv(1:r%nderiv)=v%deriv(1:v%nderiv)
   END SUBROUTINE
      
   PURE SUBROUTINE newVar3(r,value,n,list,deriv)
      TYPE(vardep),INTENT(INOUT) :: r
      INTEGER ,INTENT(IN)    :: n,list(*)
      REAL(r8),INTENT(IN)    :: value,deriv(*)
      CALL allocDeriv(r,n,list)
      r%value=value
      r%deriv(1:n)=deriv(1:n)
   END SUBROUTINE
   
   PURE ELEMENTAL SUBROUTINE newVar2(r,v1,v2)
      TYPE(vardep),INTENT(INOUT) :: r
      TYPE(vardep),INTENT(IN)    :: v1,v2
      INTEGER :: n,i1,i2,n1,n2
      INTEGER,ALLOCATABLE :: itab(:)
      n1=v1%nderiv
      n2=v2%nderiv
      ALLOCATE(itab(n1+n2))
      n=0
      i1=1
      i2=1
      DO
         IF(i1 > n1) THEN
            itab(n+1:n+1+n2-i2)=v2%index(i2:n2)
            n=n+1+n2-i2
            EXIT
         ELSE IF(i2 > n2) THEN
            itab(n+1:n+1+n1-i1)=v1%index(i1:n1)
            n=n+1+n1-i1
            EXIT
         ELSE IF(v1%index(i1) == v2%index(i2)) THEN
            n=n+1
            itab(n)=v1%index(i1)
            i1=i1+1
            i2=i2+1
         ELSE IF(v1%index(i1) < v2%index(i2)) THEN
            n=n+1
            itab(n)=v1%index(i1)
            i1=i1+1
         ELSE
            n=n+1
            itab(n)=v2%index(i2)
            i2=i2+1
         ENDIF
      ENDDO
      CALL allocDeriv(r,n,itab)
      DEALLOCATE(itab)
   END SUBROUTINE
   
   PURE ELEMENTAL SUBROUTINE newVarReal(r,c)
      TYPE(vardep),INTENT(inout) :: r
      REAL(r8),INTENT(in) :: c
      CALL allocDeriv(r,0,(/0/))
      r%value=c
   END SUBROUTINE
   
   PURE ELEMENTAL SUBROUTINE newVarInt(r,c)
      TYPE(vardep),INTENT(inout) :: r
      INTEGER     ,INTENT(in) :: c
      CALL allocDeriv(r,0,(/0/))
      r%value=c
   END SUBROUTINE
   
   PURE ELEMENTAL SUBROUTINE inverseVar(r)
      TYPE(vardep),INTENT(inout) :: r     ! the resulting variable
      INTEGER :: n
      n=r%nderiv
      r%value=1._r8/r%value
      r%deriv(1:n)=-r%deriv(1:n)*r%value**2
   END SUBROUTINE     
   
   PURE SUBROUTINE setDeriv0(r,index,deriv)
      TYPE(vardep),INTENT(inout) :: r
      INTEGER       ,INTENT(IN)    :: index
      REAL(r8)      ,INTENT(IN)    :: deriv
      INTEGER  :: n,i,itab(r%nderiv+1)
      REAL(r8) :: rtab(r%nderiv+1),value
      n=r%nderiv
      DO i=1,n
         IF(index == r%index(i)) THEN
            r%deriv(i)=r%deriv(i)+deriv
            RETURN
         ENDIF
         IF(index <  r%index(i)) EXIT
      ENDDO
      value=r%value
      itab(1:i-1)=r%index(1:i-1)
      itab(i)=index
      itab(i+1:n+1)=r%index(i:n)
      rtab(1:i-1)=r%deriv(1:i-1)
      rtab(i)=deriv
      rtab(i+1:n+1)=r%deriv(i:n)
      n=n+1
      CALL deleteVar(r)
      CALL allocDeriv(r,n,itab)
      r%value=value
      r%deriv(1:n)=rtab(1:n)
   END SUBROUTINE
   
   PURE ELEMENTAL SUBROUTINE deleteVar(r)
      TYPE(vardep),INTENT(inout) :: r
      r%value=0._r8
      r%nderiv=0
   END SUBROUTINE
   
   PURE SUBROUTINE setPartialDeriv(v1,partial,v2)
      TYPE(vardep),INTENT(inout) :: v1
      REAL(r8)    ,INTENT(in)    :: partial ! d(v1)/d(v2)
      TYPE(vardep),INTENT(in)    :: v2
      INTEGER  :: i1,i2,n,n1,n2
      REAL(r8) :: value
      INTEGER ,ALLOCATABLE :: itab(:)
      REAL(r8),ALLOCATABLE :: rtab(:)
      n1=v1%nderiv
      n2=v2%nderiv
      ALLOCATE(itab(n1+n2),rtab(n1+n2))
      n=0
      i1=1
      i2=1
      value=v1%value
      DO
         IF(i1 > n1) THEN
            itab(n+1:n+1+n2-i2)=v2%index(i2:n2)
            rtab(n+1:n+1+n2-i2)=v2%deriv(i2:n2)*partial
            n=n+1+n2-i2
            EXIT
         ELSE IF(i2 > n2) THEN
            itab(n+1:n+1+n1-i1)=v1%index(i1:n1)
            rtab(n+1:n+1+n1-i1)=v1%deriv(i1:n1)
            n=n+1+n1-i1
            EXIT
         ELSE IF(v1%index(i1) == v2%index(i2)) THEN
            n=n+1
            itab(n)=v1%index(i1)
            rtab(n)=v1%deriv(i1)+v2%deriv(i2)*partial
            i1=i1+1
            i2=i2+1
         ELSE IF(v1%index(i1) < v2%index(i2)) THEN
            n=n+1
            itab(n)=v1%index(i1)
            rtab(n)=v1%deriv(i1)
            i1=i1+1
         ELSE
            n=n+1
            itab(n)=v2%index(i2)
            rtab(n)=v2%deriv(i2)*partial
            i2=i2+1
         ENDIF
      ENDDO
      CALL newVar(v1,value,n,itab,rtab)
      DEALLOCATE(itab,rtab)
   END SUBROUTINE
   
   PURE ELEMENTAL FUNCTION VarPlusVar (v1,v2) RESULT(r)
       TYPE(vardep),INTENT(in) :: v1,v2
       TYPE(vardep)            :: r
       CALL newVar(r,v1)
       CALL addVar(r,v2)
   END FUNCTION
   
   PURE ELEMENTAL FUNCTION VarPlusReal(v,c) RESULT(r)
       TYPE(vardep),INTENT(in) :: v
       REAL(r8) ,INTENT(in) :: c
       TYPE(vardep)            :: r
       CALL newVar(r,v)
       CALL addVar(r,c)
   END FUNCTION
   
   PURE ELEMENTAL FUNCTION VarPlusInt(v,c) RESULT(r)
       TYPE(vardep),INTENT(in) :: v
       INTEGER     ,INTENT(in) :: c
       TYPE(vardep)            :: r
       CALL newVar(r,v)
       CALL addVar(r,c)
   END FUNCTION
   
   PURE ELEMENTAL FUNCTION RealPlusVar (c,v) RESULT(r)
       TYPE(vardep),INTENT(in) :: v
       REAL(r8) ,INTENT(in) :: c
       TYPE(vardep)            :: r
       CALL newVar(r,v)
       CALL addVar(r,c)
   END FUNCTION
   
   PURE ELEMENTAL FUNCTION IntPlusVar (c,v) RESULT(r)
       TYPE(vardep),INTENT(in) :: v
       INTEGER     ,INTENT(in) :: c
       TYPE(vardep)            :: r
       CALL newVar(r,v)
       CALL addVar(r,c)
   END FUNCTION
   
   PURE ELEMENTAL FUNCTION VarMultVar (v1,v2) RESULT(r)
       TYPE(vardep),INTENT(in) :: v1,v2
       TYPE(vardep)            :: r
       CALL newVar (r,v1)
       CALL multVar(r,v2)
   END FUNCTION
   
   PURE ELEMENTAL FUNCTION VarMultReal(v,c) RESULT(r)
       TYPE(vardep),INTENT(in) :: v
       REAL(r8) ,INTENT(in) :: c
       TYPE(vardep)            :: r
       CALL newVar(r,v)
       CALL multVar(r,c)
   END FUNCTION
   
   PURE ELEMENTAL FUNCTION VarMultInt(v,c) RESULT(r)
       TYPE(vardep),INTENT(in) :: v
       INTEGER     ,INTENT(in) :: c
       TYPE(vardep)            :: r
       CALL newVar(r,v)
       CALL multVar(r,c)
   END FUNCTION
   
   PURE ELEMENTAL FUNCTION RealMultVar (c,v) RESULT(r)
       TYPE(vardep),INTENT(in) :: v
       REAL(r8) ,INTENT(in) :: c
       TYPE(vardep)            :: r
       CALL newVar(r,v)
       CALL multVar(r,c)
   END FUNCTION
   
   PURE ELEMENTAL FUNCTION IntMultVar(c,v) RESULT(r)
       TYPE(vardep),INTENT(in) :: v
       INTEGER     ,INTENT(in) :: c
       TYPE(vardep)            :: r
       CALL newVar(r,v)
       CALL multVar(r,c)
   END FUNCTION
   
   PURE ELEMENTAL FUNCTION VarDivVar(v1,v2) RESULT(r)
       TYPE(vardep),INTENT(in) :: v1,v2
       TYPE(vardep)            :: r
       CALL newVar(r,v1)
       CALL divVar(r,v2)
   END FUNCTION
   
   PURE ELEMENTAL FUNCTION VarDivReal(v,c) RESULT(r)
       TYPE(vardep),INTENT(in) :: v
       REAL(r8) ,INTENT(in) :: c
       TYPE(vardep)            :: r
       CALL newVar(r,v)
       CALL divVar(r,c)
   END FUNCTION
   
   PURE ELEMENTAL FUNCTION VarDivInt(v,c) RESULT(r)
       TYPE(vardep),INTENT(in) :: v
       INTEGER     ,INTENT(in) :: c
       TYPE(vardep)            :: r
       CALL newVar(r,v)
       CALL divVar(r,c)
   END FUNCTION
   
   PURE ELEMENTAL FUNCTION RealDivVar(c,v) RESULT(r)
       TYPE(vardep),INTENT(in) :: v
       REAL(r8) ,INTENT(in) :: c
       TYPE(vardep)            :: r
       CALL newVar(r,v)
       CALL inverseVar(r)
       CALL multVar(r,c)
   END FUNCTION
   
   PURE ELEMENTAL FUNCTION IntDivVar(c,v) RESULT(r)
       TYPE(vardep),INTENT(in) :: v
       INTEGER     ,INTENT(in) :: c
       TYPE(vardep)            :: r
       CALL newVar(r,v)
       CALL inverseVar(r)
       CALL multVar(r,c)
   END FUNCTION
   
   PURE ELEMENTAL FUNCTION VarMinusVar(v1,v2) RESULT(r)
       TYPE(vardep),INTENT(in) :: v1,v2
       TYPE(vardep)            :: r
       CALL newVar(r,v1)
       CALL minusVar(r,v2)
   END FUNCTION
   
   PURE ELEMENTAL FUNCTION VarMinusReal(v,c) RESULT(r)
       TYPE(vardep),INTENT(in) :: v
       REAL(r8) ,INTENT(in) :: c
       TYPE(vardep)            :: r
       CALL newVar(r,v)
       CALL addVar(r,-c)
   END FUNCTION
   
   PURE ELEMENTAL FUNCTION VarMinusInt(v,c) RESULT(r)
       TYPE(vardep),INTENT(in) :: v
       INTEGER     ,INTENT(in) :: c
       TYPE(vardep)            :: r
       CALL newVar(r,v)
       CALL addVar(r,-c)
   END FUNCTION
   
   PURE ELEMENTAL FUNCTION RealMinusVar(c,v) RESULT(r)
       TYPE(vardep),INTENT(in) :: v
       REAL(r8)    ,INTENT(in) :: c
       TYPE(vardep)            :: r
       CALL newVar(r,v)
       CALL multVar(r,-1._r8)
       CALL addVar(r,c)
   END FUNCTION
   
   PURE ELEMENTAL FUNCTION IntMinusVar(c,v) RESULT(r)
       TYPE(vardep),INTENT(in) :: v
       INTEGER     ,INTENT(in) :: c
       TYPE(vardep)            :: r
       CALL newVar(r,v)
       CALL multVar(r,-1._r8)
       CALL addVar(r,c)
   END FUNCTION
   
   PURE SUBROUTINE deleteIndex(var,index)
      TYPE(vardep),INTENT(inout) :: var
      INTEGER  ,INTENT(in)    :: index
      INTEGER :: i,n
      n=var%nderiv
      DO i=1,n
         IF(var % index(i) == index) THEN
            var % index(i:n-1) = var % index(i+1:n)
            var % deriv(i:n-1) = var % deriv(i+1:n)
            var % nderiv       = n-1
            EXIT
         ENDIF
      ENDDO
   END SUBROUTINE
   
   PURE SUBROUTINE replaceIndex(var,index,ref) ! replacing a particular index by a variable
      TYPE(vardep),INTENT(inout) :: var   ! the variable to modify
      INTEGER  ,INTENT(in)    :: index ! the index to replace
      TYPE(vardep),INTENT(in)    :: ref   ! new reference variable
      INTEGER  :: i,n
      REAL(r8) :: coef,value
      n=var%nderiv
      value=var%value
      DO i=1,n
         IF(var % index(i) == index) THEN
            coef=var%deriv(i)
            CALL deleteIndex(var,index)
            var = var + (coef * ref)
            EXIT
         ENDIF
      ENDDO
      var%value=value
   END SUBROUTINE
   
   PURE ELEMENTAL FUNCTION minVarVar(v1,v2) RESULT(v)   ! computing v=MIN(v1,v2)
      TYPE(vardep),INTENT(in) :: v1,v2
      TYPE(vardep)            :: v
      IF(v1%value < v2%value) THEN
         v=v1
      ELSE
         v=v2
      ENDIF
   END FUNCTION
   
   PURE ELEMENTAL FUNCTION minVarReal(v1,c2) RESULT(v)   ! computing v=MIN(v1,c2)
      TYPE(vardep),INTENT(in) :: v1
      REAL(r8) ,INTENT(in) :: c2
      TYPE(vardep)            :: v
      IF(v1%value < c2) THEN
         v=v1
      ELSE
         v=c2
      ENDIF
   END FUNCTION
   
   PURE ELEMENTAL FUNCTION minVarInt(v1,c2) RESULT(v)   ! computing v=MIN(v1,c2)
      TYPE(vardep),INTENT(in) :: v1
      INTEGER     ,INTENT(in) :: c2
      TYPE(vardep)            :: v
      IF(v1%value < c2) THEN
         v=v1
      ELSE
         v=c2
      ENDIF
   END FUNCTION
   
   PURE ELEMENTAL FUNCTION minRealVar(c1,v2) RESULT(v)   ! computing v=MIN(c1,v2)
      REAL(r8) ,INTENT(in) :: c1
      TYPE(vardep),INTENT(in) :: v2
      TYPE(vardep)            :: v
      IF(c1 < v2%value) THEN
         v=c1
      ELSE
         v=v2
      ENDIF
   END FUNCTION
   
   PURE ELEMENTAL FUNCTION minIntVar(c1,v2) RESULT(v)   ! computing v=MIN(c1,v2)
      INTEGER     ,INTENT(in) :: c1
      TYPE(vardep),INTENT(in) :: v2
      TYPE(vardep)            :: v
      IF(c1 < v2%value) THEN
         v=c1
      ELSE
         v=v2
      ENDIF
   END FUNCTION
   
   PURE ELEMENTAL FUNCTION maxVarVar(v1,v2) RESULT(v)   ! computing v=MAX(v1,v2)
      TYPE(vardep),INTENT(in) :: v1,v2
      TYPE(vardep)            :: v
      IF(v1%value > V2%value) THEN
         v=v1
      ELSE
         v=v2
      ENDIF
   END FUNCTION
   
   PURE ELEMENTAL FUNCTION maxVarReal(v1,c2) RESULT(v)   ! computing v=MAX(v1,c2)
      TYPE(vardep),INTENT(in) :: v1
      REAL(r8)    ,INTENT(in) :: c2
      TYPE(vardep) :: v
      IF(v1%value > c2) THEN
         v=v1
      ELSE
         v=c2
      ENDIF
   END FUNCTION
   
   PURE ELEMENTAL FUNCTION maxVarInt(v1,c2) RESULT(v)   ! computing v=MAX(v1,c2)
      TYPE(vardep),INTENT(in) :: v1
      INTEGER     ,INTENT(in) :: c2
      TYPE(vardep) :: v
      IF(v1%value > c2) THEN
         v=v1
      ELSE
         v=c2
      ENDIF
   END FUNCTION
   
   PURE ELEMENTAL FUNCTION maxRealVar(c1,v2) RESULT(v)   ! computing v=MAX(c1,v2)
      REAL(r8)    ,INTENT(in) :: c1
      TYPE(vardep),INTENT(in) :: v2
      TYPE(vardep) :: v
      IF(c1 > v2%value) THEN
         v=c1
      ELSE
         v=v2
      ENDIF
   END FUNCTION
   
   PURE ELEMENTAL FUNCTION maxIntVar(c1,v2) RESULT(v)   ! computing v=MAX(c1,v2)
      INTEGER     ,INTENT(in) :: c1
      TYPE(vardep),INTENT(in) :: v2
      TYPE(vardep) :: v
      IF(c1 > v2%value) THEN
         v=c1
      ELSE
         v=v2
      ENDIF
   END FUNCTION
   
   PURE ELEMENTAL FUNCTION sinFun(c) RESULT (v)
      TYPE(vardep),INTENT(in) :: c
      TYPE(vardep)            :: v
      INTEGER  :: i
      REAL(r8) :: d
      CALL newVar(v,c)
      d=cos(c%value)
      v%value=sin(c%value)
      DO i=1,c%nderiv
        v%deriv(i)=c%deriv(i)*d
      ENDDO
   END FUNCTION
   
   PURE ELEMENTAL SUBROUTINE sinSub(c)
      TYPE(vardep),INTENT(inout) :: c
      INTEGER  :: i
      REAL(r8) :: d
      d=cos(c%value)
      c%value=sin(c%value)
      DO i=1,c%nderiv
        c%deriv(i)=c%deriv(i)*d
      ENDDO
   END SUBROUTINE
   
   PURE ELEMENTAL FUNCTION cosFun(c) RESULT (v)
      TYPE(vardep),INTENT(in) :: c
      TYPE(vardep)            :: v
      INTEGER  :: i
      REAL(r8) :: d
      CALL newVar(v,c)
      d=-sin(c%value)
      v%value=cos(c%value)
      DO i=1,c%nderiv
        v%deriv(i)=c%deriv(i)*d
      ENDDO
   END FUNCTION
   
   PURE ELEMENTAL SUBROUTINE cosSub(c)
      TYPE(vardep),INTENT(inout) :: c
      INTEGER  :: i
      REAL(r8) :: d
      d=-sin(c%value)
      c%value=cos(c%value)
      DO i=1,c%nderiv
        c%deriv(i)=c%deriv(i)*d
      ENDDO
   END SUBROUTINE
   
   PURE ELEMENTAL FUNCTION tanFun(c) RESULT (v)
      TYPE(vardep),INTENT(in) :: c
      TYPE(vardep)            :: v
      INTEGER  :: i
      REAL(r8) :: d
      CALL newVar(v,c)
      v%value=tan(c%value)
      d=1+v%value**2
      DO i=1,c%nderiv
        v%deriv(i)=c%deriv(i)*d
      ENDDO
   END FUNCTION
   
   PURE ELEMENTAL SUBROUTINE tanSub(c)
      TYPE(vardep),INTENT(inout) :: c
      INTEGER  :: i
      REAL(r8) :: d
      c%value=tan(c%value)
      d=1+c%value**2
      DO i=1,c%nderiv
        c%deriv(i)=c%deriv(i)*d
      ENDDO
   END SUBROUTINE
   
   PURE ELEMENTAL FUNCTION tanhFun(c) RESULT (v)
      TYPE(vardep),INTENT(in) :: c
      TYPE(vardep)            :: v
      INTEGER  :: i
      REAL(r8) :: d
      CALL newVar(v,c)
      v%value=tanh(c%value)
      d=1+v%value**2
      DO i=1,c%nderiv
        v%deriv(i)=c%deriv(i)*d
      ENDDO
   END FUNCTION
   
   PURE ELEMENTAL SUBROUTINE tanhSub(c)
      TYPE(vardep),INTENT(inout) :: c
      INTEGER  :: i
      REAL(r8) :: d
      c%value=tanh(c%value)
      d=1+c%value**2
      DO i=1,c%nderiv
        c%deriv(i)=c%deriv(i)*d
      ENDDO
   END SUBROUTINE
   
   PURE ELEMENTAL FUNCTION sqrtFun(c) RESULT (v)
      TYPE(vardep),INTENT(in) :: c
      TYPE(vardep)            :: v
      INTEGER  :: i
      REAL(r8) :: d
      CALL newVar(v,c)
      v%value=sqrt(c%value)
      d=1/(2*v%value)
      DO i=1,c%nderiv
        v%deriv(i)=c%deriv(i)*d
      ENDDO
   END FUNCTION
   
   PURE ELEMENTAL SUBROUTINE sqrtSub(c)
      TYPE(vardep),INTENT(inout) :: c
      INTEGER  :: i
      REAL(r8) :: d
      c%value=sqrt(c%value)
      d=1/(2*c%value)
      DO i=1,c%nderiv
        c%deriv(i)=c%deriv(i)*d
      ENDDO
   END SUBROUTINE
   
   PURE ELEMENTAL FUNCTION logFun(c) RESULT (v)
      TYPE(vardep),INTENT(in) :: c
      TYPE(vardep)            :: v
      INTEGER  :: i
      REAL(r8) :: d
      CALL newVar(v,c)
      v%value=log(c%value)
      d=1/c%value
      DO i=1,c%nderiv
        v%deriv(i)=c%deriv(i)*d
      ENDDO
   END FUNCTION
   
   PURE ELEMENTAL SUBROUTINE logSub(c)
      TYPE(vardep),INTENT(inout) :: c
      INTEGER  :: i
      REAL(r8) :: d
      d=1/c%value
      c%value=log(c%value)
      DO i=1,c%nderiv
        c%deriv(i)=c%deriv(i)*d
      ENDDO
   END SUBROUTINE
   
   PURE ELEMENTAL FUNCTION expFun(c) RESULT (v)
      TYPE(vardep),INTENT(in) :: c
      TYPE(vardep)            :: v
      INTEGER  :: i
      REAL(r8) :: d
      CALL newVar(v,c)
      v%value=exp(c%value)
      d=v%value
      DO i=1,c%nderiv
        v%deriv(i)=c%deriv(i)*d
      ENDDO
   END FUNCTION
   
   PURE ELEMENTAL SUBROUTINE expSub(c)
      TYPE(vardep),INTENT(inout) :: c
      INTEGER  :: i
      REAL(r8) :: d
      c%value=exp(c%value)
      d=c%value
      DO i=1,c%nderiv
        c%deriv(i)=c%deriv(i)*d
      ENDDO
   END SUBROUTINE
   
   PURE ELEMENTAL FUNCTION atanFun(c) RESULT (v)
      TYPE(vardep),INTENT(in) :: c
      TYPE(vardep)            :: v
      INTEGER  :: i
      REAL(r8) :: d
      CALL newVar(v,c)
      v%value=atan(c%value)
      d=1/(1+v%value**2)
      DO i=1,c%nderiv
        v%deriv(i)=c%deriv(i)*d
      ENDDO
   END FUNCTION
   
   PURE ELEMENTAL SUBROUTINE atanSub(c)
      TYPE(vardep),INTENT(inout) :: c
      INTEGER  :: i
      REAL(r8) :: d
      c%value=atan(c%value)
      d=1/(1+c%value**2)
      DO i=1,c%nderiv
        c%deriv(i)=c%deriv(i)*d
      ENDDO
   END SUBROUTINE
   
   PURE ELEMENTAL FUNCTION absFun(c) RESULT (v)
      TYPE(vardep),INTENT(in) :: c
      TYPE(vardep)            :: v
      IF(c%value >= 0) THEN
        v=c
      ELSE
        v=(-1._r8)*c
      ENDIF
   END FUNCTION
   
   PURE ELEMENTAL SUBROUTINE absSub(c)
      TYPE(vardep),INTENT(inout) :: c
      IF(c%value >= 0) RETURN
      c%value=-c%value
      c%deriv(1:c%nderiv)=-c%deriv(1:c%nderiv)
   END SUBROUTINE
   
   PURE FUNCTION VarEqVar(v1,v2) RESULT(r)
      TYPE(vardep),INTENT(in) :: v1,v2
      LOGICAL                 :: r
      r=v1%value == v2%value
   END FUNCTION
   
   PURE FUNCTION VarEqReal(v1,v2) RESULT(r)
      TYPE(vardep),INTENT(in) :: v1
      REAL(r8)    ,INTENT(in) :: v2
      LOGICAL                 :: r
      r=v1%value == v2
   END FUNCTION
   
   PURE FUNCTION VarEqInt(v1,v2) RESULT(r)
      TYPE(vardep),INTENT(in) :: v1
      INTEGER     ,INTENT(in) :: v2
      LOGICAL                 :: r
      r=v1%value == v2
   END FUNCTION
   
   PURE FUNCTION RealEqVar(v1,v2) RESULT(r)
      REAL(r8)    ,INTENT(in) :: v1
      TYPE(vardep),INTENT(in) :: v2
      LOGICAL                 :: r
      r=v1 == v2%value
   END FUNCTION
   
   PURE FUNCTION IntEqVar(v1,v2) RESULT(r)
      INTEGER     ,INTENT(in) :: v1
      TYPE(vardep),INTENT(in) :: v2
      LOGICAL                 :: r
      r=v1 == v2%value
   END FUNCTION
   
   PURE FUNCTION VarNeVar(v1,v2) RESULT(r)
      TYPE(vardep),INTENT(in) :: v1,v2
      LOGICAL                 :: r
      r=v1%value /= v2%value
   END FUNCTION
   
   PURE FUNCTION VarNeReal(v1,v2) RESULT(r)
      TYPE(vardep),INTENT(in) :: v1
      REAL(r8)    ,INTENT(in) :: v2
      LOGICAL                 :: r
      r=v1%value /= v2
   END FUNCTION
   
   PURE FUNCTION VarNeInt(v1,v2) RESULT(r)
      TYPE(vardep),INTENT(in) :: v1
      INTEGER     ,INTENT(in) :: v2
      LOGICAL                 :: r
      r=v1%value /= v2
   END FUNCTION
   
   PURE FUNCTION RealNeVar(v1,v2) RESULT(r)
      REAL(r8)    ,INTENT(in) :: v1
      TYPE(vardep),INTENT(in) :: v2
      LOGICAL                 :: r
      r=v1 /= v2%value
   END FUNCTION
   
   PURE FUNCTION IntNeVar(v1,v2) RESULT(r)
      INTEGER     ,INTENT(in) :: v1
      TYPE(vardep),INTENT(in) :: v2
      LOGICAL                 :: r
      r=v1 /= v2%value
   END FUNCTION
   
   PURE FUNCTION VarLeVar(v1,v2) RESULT(r)
      TYPE(vardep),INTENT(in) :: v1,v2
      LOGICAL                 :: r
      r=v1%value <= v2%value
   END FUNCTION
   
   PURE FUNCTION VarLeReal(v1,v2) RESULT(r)
      TYPE(vardep),INTENT(in) :: v1
      REAL(r8)    ,INTENT(in) :: v2
      LOGICAL                 :: r
      r=v1%value <= v2
   END FUNCTION
   
   PURE FUNCTION VarLeInt(v1,v2) RESULT(r)
      TYPE(vardep),INTENT(in) :: v1
      INTEGER     ,INTENT(in) :: v2
      LOGICAL                 :: r
      r=v1%value <= v2
   END FUNCTION
   
   PURE FUNCTION RealLeVar(v1,v2) RESULT(r)
      REAL(r8)    ,INTENT(in) :: v1
      TYPE(vardep),INTENT(in) :: v2
      LOGICAL                 :: r
      r=v1 <= v2%value
   END FUNCTION
   
   PURE FUNCTION IntLeVar(v1,v2) RESULT(r)
      INTEGER     ,INTENT(in) :: v1
      TYPE(vardep),INTENT(in) :: v2
      LOGICAL                 :: r
      r=v1 <= v2%value
   END FUNCTION
   
   PURE FUNCTION VarLtVar(v1,v2) RESULT(r)
      TYPE(vardep),INTENT(in) :: v1,v2
      LOGICAL                 :: r
      r=v1%value < v2%value
   END FUNCTION
   
   PURE FUNCTION VarLtReal(v1,v2) RESULT(r)
      TYPE(vardep),INTENT(in) :: v1
      REAL(r8)    ,INTENT(in) :: v2
      LOGICAL                 :: r
      r=v1%value < v2
   END FUNCTION
   
   PURE FUNCTION VarLtInt(v1,v2) RESULT(r)
      TYPE(vardep),INTENT(in) :: v1
      INTEGER     ,INTENT(in) :: v2
      LOGICAL                 :: r
      r=v1%value < v2
   END FUNCTION
   
   PURE FUNCTION RealLtVar(v1,v2) RESULT(r)
      REAL(r8)    ,INTENT(in) :: v1
      TYPE(vardep),INTENT(in) :: v2
      LOGICAL                 :: r
      r=v1 < v2%value
   END FUNCTION
   
   PURE FUNCTION IntLtVar(v1,v2) RESULT(r)
      INTEGER     ,INTENT(in) :: v1
      TYPE(vardep),INTENT(in) :: v2
      LOGICAL                 :: r
      r=v1 < v2%value
   END FUNCTION
   
   PURE FUNCTION VarGtVar(v1,v2) RESULT(r)
      TYPE(vardep),INTENT(in) :: v1,v2
      LOGICAL                 :: r
      r=v1%value > v2%value
   END FUNCTION
   
   PURE FUNCTION VarGtReal(v1,v2) RESULT(r)
      TYPE(vardep),INTENT(in) :: v1
      REAL(r8)    ,INTENT(in) :: v2
      LOGICAL                 :: r
      r=v1%value > v2
   END FUNCTION
   
   PURE FUNCTION VarGtInt(v1,v2) RESULT(r)
      TYPE(vardep),INTENT(in) :: v1
      INTEGER     ,INTENT(in) :: v2
      LOGICAL                 :: r
      r=v1%value > v2
   END FUNCTION
   
   PURE FUNCTION RealGtVar(v1,v2) RESULT(r)
      REAL(r8)    ,INTENT(in) :: v1
      TYPE(vardep),INTENT(in) :: v2
      LOGICAL                 :: r
      r=v1 > v2%value
   END FUNCTION
   
   PURE FUNCTION IntGtVar(v1,v2) RESULT(r)
      INTEGER     ,INTENT(in) :: v1
      TYPE(vardep),INTENT(in) :: v2
      LOGICAL                 :: r
      r=v1 > v2%value
   END FUNCTION
   
   PURE FUNCTION VarGeVar(v1,v2) RESULT(r)
      TYPE(vardep),INTENT(in) :: v1,v2
      LOGICAL                 :: r
      r=v1%value >= v2%value
   END FUNCTION
   
   PURE FUNCTION VarGeReal(v1,v2) RESULT(r)
      TYPE(vardep),INTENT(in) :: v1
      REAL(r8)    ,INTENT(in) :: v2
      LOGICAL                 :: r
      r=v1%value >= v2
   END FUNCTION
   
   PURE FUNCTION VarGeInt(v1,v2) RESULT(r)
      TYPE(vardep),INTENT(in) :: v1
      INTEGER     ,INTENT(in) :: v2
      LOGICAL                 :: r
      r=v1%value >= v2
   END FUNCTION
   
   PURE FUNCTION RealGeVar(v1,v2) RESULT(r)
      REAL(r8)    ,INTENT(in) :: v1
      TYPE(vardep),INTENT(in) :: v2
      LOGICAL                 :: r
      r=v1 >= v2%value
   END FUNCTION
   
   PURE FUNCTION IntGeVar(v1,v2) RESULT(r)
      INTEGER     ,INTENT(in) :: v1
      TYPE(vardep),INTENT(in) :: v2
      LOGICAL                 :: r
      r=v1 >= v2%value
   END FUNCTION

END MODULE
