eiram.f90

Go to the documentation of this file.

!    Copyright (c) 2016  Forschungszentrum Juelich GmbH
!                        Markus Brenneis, Vladislav Kotov 
!
!    This file is part of EIRAM.
!
!    EIRAM is free software: you can redistribute it and/or modify
!    it under the terms of the GNU General Public License as published by
!    the Free Software Foundation, either version 3 of the License, or
!    (at your option) any later version.
!
!    EIRAM is distributed in the hope that it will be useful,
!    but WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!    GNU General Public License for more details.
!
!    You should have received a copy of the GNU General Public License
!    along with EIRAM.  If not, see <http://www.gnu.org/licenses/>.
!

    
 function eiram_fit(B,N,X) 
    implicit none  
    integer,parameter :: DP=selected_real_kind(13)
    real(DP) :: eiram_fit 
    real(DP),intent(in) :: B(0:n)
    integer,intent(in) :: N
    real(DP),intent(in) :: X

    integer :: i

    eiram_fit=b(n)
    do i=n-1,0,-1   
     eiram_fit=eiram_fit*x+b(i)
    end do 

 end function eiram_fit

module eiram
    implicit none

    public eiram_load        
    public eiram_deallocate  
    public eiram_get_id       
    public eiram_get_order
    public eiram_create_image
    public eiram_calc1_fast
    public eiram_calc1        
    public eiram_calc12   
    public eiram_calc2   

    integer, parameter, public :: eiram_dp = selected_real_kind(13)
    integer, public :: eiram_unit = 6
    
    integer,parameter :: units_length=32
    integer,parameter :: quantity_length=64
    integer,parameter :: reaction_length=128
    integer,parameter :: description_length=256
    integer, public, parameter :: eiram_file_name_length = 32
    integer, public, parameter :: eiram_type_length = 8
    integer, public, parameter :: eiram_index_length = 16
    integer, parameter :: max_line_length = 256

    logical, parameter :: debug = .false.

    external eiram_fit
    real(eiram_dp),public :: eiram_fit
 
    type, public:: eiram_data
        integer :: n
        integer :: m
        real(EIRAM_DP), allocatable :: creac(:,:)
        real(EIRAM_DP) :: x1min 
        real(EIRAM_DP) :: x1max 
        real(EIRAM_DP) :: x2min 
        real(EIRAM_DP) :: x2max
        character(UNITS_LENGTH) :: x1units = ""
        character(UNITS_LENGTH) ::  x2units = "" 
         character(UNITS_LENGTH) ::  yunits = ""
        character(QUANTITY_LENGTH) :: x1quantity = ""
        character(QUANTITY_LENGTH) :: x2quantity = "" 
        character(QUANTITY_LENGTH) :: yquantity = ""
        character(REACTION_LENGTH) :: reaction = ""
        character(DESCRIPTION_LENGTH) :: description = ""
        character(EIRAM_FILE_NAME_LENGTH) :: file = ""
        character(EIRAM_TYPE_LENGTH) :: TYPE = ""
        character(EIRAM_INDEX_LENGTH) :: index = ""
    end type eiram_data

    private

    integer,public :: eiram_nreact=0

    type(eiram_data), allocatable :: reactions(:)
    
    contains

!***********************************************************************************
!   INITIALIZATION/DEINITIALIZATION ROUTINES  
!***********************************************************************************
    
    subroutine eiram_load(filePath, fileName, err)

        character(*), intent(in) :: filePath
        character(*), intent(in) :: fileName
        integer, intent(out) :: err

        logical :: fileExists,isLaTeX
        integer :: nData, nOldData
        integer :: l,nDataMax
        integer, parameter :: UNIT = 20
        character(MAX_LINE_LENGTH) :: line
        character(*), parameter :: BEGIN_DATA_MARK = "##BEGIN DATA HERE##"
        !string which is used to find out that the document is in LaTeX format
        ! - w/o '\' at the beginning because for some reason this confuses g95 compiler
        character(*), parameter :: LATEX_MARK = "begin{document}"
        type(eiram_data), allocatable:: reactions_(:)

        intrinsic index,trim

         if(debug) print *, "open file ",trim(filename)

        inquire(file=trim(filepath)//trim(filename),exist=fileexists,iostat=err)
        if(.not.fileexists .or. err /= 0) then
         err=310
         write(eiram_unit, *) "ERROR in eiram_load: input file ", &
                              trim(filepath)//trim(filename)," is not found" 
         return
        end if

        if(allocated(reactions)) then
            nolddata = size(reactions)
        else
            nolddata = 0
        end if
        
        ndatamax = countlinescontaining(trim(filepath)//trim(filename), "Reaction", err)
        if(err/=0) then
            write(eiram_unit, *) "ERROR in eiram_load while counting lines, err ", err
            err = 300
            return
        end if
        
        allocate(reactions_(ndatamax+nolddata),stat=err)
        if(err/=0) then
            write(eiram_unit, *) "ERROR in eiram_load while allocating temporary array, ",&
                                 "nData",ndatamax+nolddata
            err = 200
            return
        end if

        open(unit, iostat=err, file=trim(filepath)//trim(filename), status='old', action='read')
        if(err /= 0) then
            write(eiram_unit, *) "ERROR in eiram_load while reading "//trim(filepath)//trim(filename)
            close(unit)
            err = 300
            return
        end if

        l = 0
        ndata = 0
        islatex=.false.        

        !look for beginning of the data section
        do
          read(unit, '(A)', iostat=err,end=1000) line
          l = l + 1
          if(err /= 0) goto 1300

          if(index(line,latex_mark) > 0) islatex=.true.

          if(index(line, begin_data_mark) > 0) then
            if(debug) print *, "found BEGIN_DATA_MARK in line ", l
            exit 
          end if         
        end do        

       !read the data 
       if(islatex) then 
        call eiram_read_latex(reactions_,filename,unit,l,ndata,*1300,*1303,*1304)        
       else
        call eiram_read_txt(reactions_,filename,unit,l,ndata,*1300,*1303,*1304)        
       end if
     
       if(debug) print *, "EOF in", l
       close(unit)
        
       if(allocated(reactions)) then
          !Add reactions which were already loaded in previous calls of eiram_load
           reactions_(ndata+1:ndata+nolddata) = reactions(1:nolddata)
           deallocate(reactions, stat=err)
           if(err /= 0) then
             write(eiram_unit, *) "ERROR in eiram_load: cannot re-allocate memory ", err
             err = 200
            return
           end if
       end if
        
        !(Re-)allocate REACTIONS
        allocate(reactions(ndata+nolddata),stat=err)
        if(err /= 0) then
         write(eiram_unit, *) "ERROR in eiram_load: cannot allocate the array of reactions ",&
                              "nData ",ndata+nolddata
         err = 200
         return
        end if

        reactions(1:ndata+nolddata) = reactions_(1:ndata+nolddata)
 
        call deallocate_reactions

        write(eiram_unit, *) "initalization from file "//trim(filepath)//trim(filename)//" completed"
        write(eiram_unit, *) "  new data sets loaded: ", ndata
        eiram_nreact=ndata+nolddata
        write(eiram_unit, *) "  total number of data sets: ", eiram_nreact
        return

1000    write(eiram_unit, *) "ERROR in eiram_load:  no data section is found in file "//trim(filepath)//trim(filename)
        close(unit)
        call deallocate_reactions
        err = 300
        return
1300    write(eiram_unit, *) "ERROR in eiram_load:  while reading line ", l, " of "//trim(filepath)//trim(filename)
        close(unit)
        call deallocate_reactions
        err = 300
        return
1303    write(eiram_unit, *) "ERROR in eiram_load: wrong format in line ", l, & 
                             " of "//trim(filepath)//trim(filename), err
        close(unit)
        call deallocate_reactions
        err = 300
        return
1304    write(eiram_unit, *) "ERROR in eiram_load: cannot allocate CREAC "
        write(eiram_unit, *) " line ", l, " of ",trim(filepath)//trim(filename),&
                             " reaction ",trim(reactions_(ndata)%TYPE)//" "//trim(reactions_(ndata)%INDEX)
        close(unit)
        call deallocate_reactions
        err = 200
        return       

    contains 

    !clean-up for the case if the compiler uses argument -save       
     subroutine deallocate_reactions
      integer :: i
      logical :: lerr
      lerr=.false.
      do i=1,size(reactions_)
       if(allocated(reactions_(i)%CREAC)) then
         deallocate(reactions_(i)%CREAC,stat=err)
         if(err /= 0) lerr=.true.
       end if  
      end do   
      deallocate(reactions_,stat=err)
      if(err /= 0 .or. lerr ) then
        write(eiram_unit, *) "WARNING from eiram_load: problem while deallocating temporary array "
        err = -200          
      end if
     end subroutine deallocate_reactions

    end subroutine eiram_load

    subroutine eiram_read_latex(reactions_,fileName,unit,l,nData,*,*,*)
     type(eiram_data), allocatable:: reactions_(:)
     character(*), intent(in) :: fileName
     integer,intent(in) :: unit  
     integer,intent(inout) :: l 
     integer,intent(out) :: nData    
     integer :: gotoLine,numberOfUnclosedBraces,n,err,i1,i2
     character(EIRAM_TYPE_LENGTH) :: stype
     character(MAX_LINE_LENGTH) :: line
     character(MAX_LINE_LENGTH*2) :: reactionName

     intrinsic index,count,transfer,trim,adjustl,len_trim

       ndata=0

       foreachline: do
       read(unit, '(A)', iostat=err,end=100) line
       l = l + 1
       if(err /= 0) return 1
       newsection: if(index(line, "\section") == 1) then
         if(debug) print *, "found \section in line ", l
         if(index(trim(line), " ") == 0) then
           read(unit, '(A)', iostat=err) line
           l = l + 1
           if(err /= 0) return 1
          end if
          stype= line(index(line,"{")+1:index(line," ")-1)
          if(index(stype,":")>0) stype=stype(1:index(stype,":")-1)
       end if newsection

       newreaction: if(index(line, "subsection{") > 0) then
            numberofunclosedbraces = count(transfer(line, 'a', max_line_length) == "{")
            n = count(transfer(line, 'a', max_line_length) == "}")
            numberofunclosedbraces = numberofunclosedbraces - n
            reactionname = trim(adjustl(line))
 
            !for the case when the reaction is defined like this 
            !\subsection{
            !             ............ }              
            findendofreactionname: do
                if(numberofunclosedbraces == 0) exit findendofreactionname
                read(unit, '(A)', iostat=err) line
                l = l + 1
                if(err /= 0) return 1
                n = count(transfer(line, 'a', max_line_length) == "{")
                numberofunclosedbraces = numberofunclosedbraces + n
                n = count(transfer(line, 'a', max_line_length) == "}")
                numberofunclosedbraces = numberofunclosedbraces - n
                if(len_trim(reactionname) /= index(reactionname, "{")) then
                 reactionname = trim(reactionname)//" "//trim(adjustl(line))
                else
                 reactionname = trim(reactionname)//trim(adjustl(line))
                end if
                        if(line == "}") exit findendofreactionname
             end do findendofreactionname
                          
             if(index(reactionname,'Reaction')<1) cycle

             ndata = ndata+1

             reactions_(ndata)%FILE = filename 
             reactions_(ndata)%TYPE = stype

             if(debug) print *, "found Reaction in line ", l
             !Sub-string between 1st and 2nd space after "Reaction": the index 
             i1=index(reactionname,"Reaction")+8;
             reactions_(ndata)%INDEX = reactionname(index(reactionname(i1:)," ")+i1:)
             reactions_(ndata)%INDEX = reactions_(ndata)%INDEX(:index(reactions_(ndata)%INDEX," ")-1)
             if(debug) print *, "  with index ", trim(reactions_(ndata)%INDEX)
             !Sub-string between $ and $ space: reaction equation (if available)
             i1=index(reactionname,"$")
             i2=index(reactionname,"$",.true.)
             if (i1>0 .and. i2>0 ) then
               reactions_(ndata)%REACTION = reactionname(i1:i2)
               if(debug) print *, "  with equation ", trim(reactions_(ndata)%REACTION)
             end if 
                   
             skiptillverbatim: do
                read(unit, '(A)', iostat=err) line
                l = l + 1
                if(err /= 0) return 1
                if(index(line, "verbatim") > 0) then
                   exit skiptillverbatim
                end if
             end do skiptillverbatim
                                       
             select case(reactions_(ndata)%TYPE)
                    case("H.1", "H.2", "H.8", "H.11")
                         if(debug) print *, "h1/2/8/11"                
                         call parse_single(reactions_(ndata),unit,gotoline,l)
                         if(gotoline == 1300) return 1
                         if(gotoline == 1303) return 2
                         if(gotoline == 1304) return 3
                         if(.not.allocated(reactions_(ndata)%CREAC)) ndata=ndata-1 !skip this reaction
                    case ("H.3", "H.4", "H.9", "H.10", "H.12")
                         if(debug) print *, "h3/4/9/10/12"
                         call parse_double(reactions_(ndata),unit,gotoline,l)
                         if(gotoline == 1300) return 1
                         if(gotoline == 1303) return 2
                         if(gotoline == 1304) return 3
                         if(.not.allocated(reactions_(ndata)%CREAC)) ndata=ndata-1 !skip this reaction
                    case default
       !                  write(eiram_unit, *) "WARNING from eiram_load: No handle for type "&
       !                                       //trim(reactions_(nData)%TYPE)//" in line", l
                          if(debug) print *, ">> exiting"
                          ndata = ndata-1 !skip reaction
             end select                    
         end if newreaction
            
        end do foreachline

  100   continue

    end subroutine eiram_read_latex

    subroutine eiram_read_txt(reactions_,fileName,unit,l,nData,*,*,*)
     type(eiram_data), allocatable:: reactions_(:)
     character(*), intent(in) :: fileName
     integer,intent(in) :: unit  
     integer,intent(inout) :: l 
     integer,intent(out) :: nData    
     integer :: gotoLine,err,i1,i2,itype
     character(EIRAM_TYPE_LENGTH) :: stype
     character(MAX_LINE_LENGTH) :: line, reactionName

     intrinsic index,trim,adjustl,len_trim

     ndata=0
     itype=0

     foreachline: do
       read(unit, '(A)', iostat=err,end=100) line
       l = l + 1
       if(err /= 0) return 1
       call gettype(line,stype)

      newreaction: if(index(line, "Reaction") > 0) then
         ndata = ndata+1
         reactions_(ndata)%FILE = filename 
         reactions_(ndata)%TYPE = stype
         reactionname = trim(adjustl(line))
         !Sub-string between 1st and 2nd space: the index 
         i1=index(reactionname," ")+1 
         i2=index(reactionname(i1:)," ")+i1-2
         reactions_(ndata)%INDEX = reactionname(i1:i2)          
         !Sub-string between 2nd space and end of the line: reaction equation
         reactions_(ndata)%REACTION =trim(adjustl(reactionname(i2+1:)))
         if(debug) then 
          print *, "found new reaction"
          print *, "             index ",trim(reactions_(ndata)%INDEX)
          print *, "             equation ",trim(reactions_(ndata)%REACTION)
         end if
         select case(reactions_(ndata)%TYPE)
                case("H.1", "H.2", "H.8", "H.11")
                     if(debug) print *, "h1/2/8/11"                
                     call parse_single(reactions_(ndata),unit,gotoline,l)
                     if(gotoline == 1300) return 1
                     if(gotoline == 1303) return 2
                     if(gotoline == 1304) return 3
                     if(.not.allocated(reactions_(ndata)%CREAC)) ndata=ndata-1 !skip this reaction
                case ("H.3", "H.4", "H.9", "H.10", "H.12")
                     if(debug) print *, "h3/4/9/10/12"
                     call parse_double(reactions_(ndata),unit,gotoline,l)
                     if(gotoline == 1300) return 1
                     if(gotoline == 1303) return 2
                     if(gotoline == 1304) return 3
                     if(.not.allocated(reactions_(ndata)%CREAC)) ndata=ndata-1 !skip this reaction
                case default
       !              write(eiram_unit, *) "WARNING from eiram_load: No handle for type "&
       !                                   //trim(reactions_(nData)%TYPE)//" in line", l
                if(debug) print *, ">> exiting"
                ndata = ndata-1 !skip reaction
         end select     
      end if newreaction
     end do foreachline

100  continue

    contains 

    subroutine gettype(line,stype)
     character(MAX_LINE_LENGTH),intent(in) :: line
     character(EIRAM_TYPE_LENGTH),intent(inout) :: stype

     integer,parameter :: ntp=12
     integer :: i 
     character(4),parameter :: stp(0:ntp)= &
         (/'H.0 ','H.1 ','H.2 ','H.3 ','H.4 ','H.5 ','H.6 ','H.7 ','H.8 ','H.9 ','H.10','H.11','H.12'/)
          
     if(index(line,'H.')<1) return !for speed up

     do i=ntp,itype,-1
      if(index(line,trim(stp(i)))>0) then
       stype=trim(stp(i))
       itype=i
       if(debug) print *, "found new type ", stype
       return 
      end if
     end do
 
    end subroutine gettype

    end subroutine eiram_read_txt

    
    subroutine eiram_deallocate(err)

        integer, intent(out) :: err
        
        logical :: lerr 
        integer i       

        intrinsic size

        err = 0

        eiram_nreact=0

        if(allocated(reactions)) then 

         lerr=.false.
         do i=1,size(reactions)
          if(allocated(reactions(i)%CREAC)) then
           deallocate(reactions(i)%CREAC,stat=err)
           if(err /= 0) lerr=.true.
          end if  
         end do            

         deallocate(reactions, stat=err)
         if(err /= 0) then
            write(eiram_unit, *) "WARNING from eiram_deallocate: problem while deallocating memory, err ", err
            err = -200
            return
         end if

         if(lerr) then
          write(eiram_unit, *) "WARNING from eiram_deallocate: problem while deallocating memory"
          err=-200
         end if   

       end if !if(allocated(allocated(REACTIONS)) then

    end subroutine eiram_deallocate

    subroutine parse_single(reaction,unit,gotoLine,l)
     type(eiram_data) :: reaction
     integer,intent(in) :: unit 
     integer,intent(out) :: gotoLine
     integer,intent(inout) :: l  

     character(2) :: chr
     character(80) :: zeile
     integer :: err,indff,iftflg,ind,ind2,i,j
    
     intrinsic index,mod

     gotoline = 0     
             
     reaction%M = 0
     reaction%N = 8
     if(reaction%TYPE == "H.1") then
        reaction%YUNITS = "cm^2"
        reaction%YQUANTITY = "cross section"
        reaction%X1QUANTITY = "Energy"
        chr='a0'
     elseif(reaction%TYPE == "H.2") then
        reaction%YUNITS = "cm^3/s"
        reaction%YQUANTITY = "collision rate"
        reaction%X1QUANTITY = "Temperature"
        chr='b0'
     elseif(reaction%TYPE == "H.8") then
        reaction%YUNITS = "cm^3/s\cdot eV"
        reaction%YQUANTITY = "electron cooling rate"
        reaction%X1QUANTITY = "Temperature"
        chr='h0'
     elseif(reaction%TYPE == "H.11") then
        reaction%YUNITS = ""
        reaction%YQUANTITY = "ratio"
        reaction%X1QUANTITY = "Temperature"
        chr='k0'
     end if
     reaction%X1UNITS = "eV"

     !this code mimics EIRENE, slreac.f 
3    read(unit,'(A80)',iostat=err) zeile
     l=l+1
     if(err /= 0) then
       gotoline = 1300
       return
     end if

     iftflg=0

     !Skip "An analytic formula is given in the text" in LaTeX mode
     if( index(zeile,'verbatim')>0 .or. index(zeile,'subsection')>0 ) return

     indff=index(zeile,'fit-flag')
     if (index(zeile,chr)+indff.EQ.0) GOTO 3
     if(indff > 0) THEN
       read (zeile((indff+8):80),*,iostat=err) iftflg
       if(err /= 0) then
          gotoline = 1303
          return
       end if
       goto 3
      end if

      allocate(reaction%CREAC(0:8,0:0),stat=err)
      if(err /=0) then
        gotoline=1304
        return
      end if

      if(mod(iftflg,100) == 10) then
        ind=index(zeile,chr(1:1))
        reaction%creac(:,0)=0.
        read(zeile((ind+2):80),*,iostat=err) reaction%creac(0,0)
        if(err /= 0) then
          gotoline = 1303
          return
        end if
      else
        do j=0,2
         ind=index(zeile,chr(1:1))+2                   
         do i=0,2
          ind2=index(zeile(ind:80),chr(1:1))
          if(ind2<1) then
           ind2=80
          else
           ind2=ind2+ind-2
          end if 
          read(zeile(ind:ind2),*,iostat=err) reaction%creac(j*3+i,0)          
          if(err /= 0) then
           gotoline = 1303
           return
          end if 
          ind=ind2+3
         end do
         read(unit,'(A80)',iostat=err) zeile
         l=l+1
         if(err /= 0) then
          gotoline = 1300
          return
         end if
        end do
      end if
               
      if(debug) print *, reaction%CREAC         

     end subroutine parse_single

     subroutine parse_double(reaction,unit,gotoLine,l)
      type(eiram_data) :: reaction
      integer,intent(in) :: unit 
      integer,intent(out) :: gotoLine
      integer,intent(inout) :: l  
 
      character(80) :: zeile
      integer :: err,indff,iftflg,ih,i,j,k

      intrinsic index,mod

      gotoline = 0
            
      reaction%N = 8
      reaction%M = 8
      if(reaction%TYPE == "H.3") then
         reaction%YQUANTITY = "collision rate"
         reaction%YUNITS = "cm^3/s"
         reaction%X1QUANTITY = "Energy"
         reaction%X1UNITS = "eV"
         reaction%X2QUANTITY = "Temperature"
         reaction%X2UNITS = "eV"
         reaction%X2QUANTITY = "Temperature"
      elseif(reaction%TYPE == "H.4") then
         reaction%YQUANTITY = "collision rate"
         reaction%YUNITS = "cm^3/s"
         reaction%X1QUANTITY = "Density"
         reaction%X1UNITS = "cm^{-3}"
         reaction%X2QUANTITY = "Temperature"
         reaction%X2UNITS = "eV"
      elseif(reaction%TYPE == "H.9") then
         reaction%YQUANTITY = "ion cooling rate"
         reaction%YUNITS = "cm^3/s\cdot eV"
         reaction%X1QUANTITY = "Energy"
         reaction%X1UNITS = "[eV]"
         reaction%X2QUANTITY = "Temperature"
         reaction%X2UNITS = "eV"
      elseif(reaction%TYPE == "H.10") then
         reaction%YQUANTITY = "electron cooling rate"
         reaction%YUNITS = "cm^3/s\cdot eV"
         reaction%X1QUANTITY = "Density"
         reaction%X1UNITS = "cm^{-3}"
         reaction%X2QUANTITY = "Temperature"
         reaction%X2UNITS = "eV" 
      elseif(reaction%TYPE == "H.12") then
         reaction%YQUANTITY = "ratio"
         reaction%YUNITS = ""
         reaction%X1QUANTITY = "Density"
         reaction%X1UNITS = "cm^{-3}"
         reaction%X2QUANTITY = "Temperature"
         reaction%X2UNITS = "eV" 
      end if
            
      allocate(reaction%CREAC(0:8,0:8),stat=err)
      if(err /=0) then
         gotoline=1304
         return
      end if

      iftflg=0

      !this code mimics EIRENE, slreac.f 
      do j=0,2
16      read(unit,'(A80)',iostat=err) zeile
        l=l+1
        if(err /= 0) then
           gotoline = 1300
           return
        end if

       !Skip "An analytic formula is given in the text" in LaTeX mode
        if( index(zeile,'verbatim')>0 .or. index(zeile,'subsection')>0 ) return

        indff=index(zeile,'fit-flag')
        if(index(zeile,'Index')+indff==0) goto 16
        if(indff > 0) then
          read(zeile((indff+8):80),*,iostat=err) iftflg
          if(err /= 0) then
             gotoline = 1303
             return
           end if 
          goto 16
        end if

        read(unit,'(1X)',iostat=err)
        l=l+1
        if(err /= 0) then
           gotoline = 1300
           return
        end if

        if(mod(iftflg,100) == 10) then
          read(unit,*,iostat=err) ih,reaction%CREAC(0,0)
          l=l+1
          if(err /= 0) then
            gotoline = 1300
            return
          end if
          exit
        else
         do i=0,8
          read(unit,*,iostat=err) ih,(reaction%CREAC(i,k),k=j*3,j*3+2)
          if(err /= 0) then
            gotoline = 1300
            return
          end if
         end do
        end if
      end do
           
      if(debug) print *, reaction%CREAC

      end subroutine parse_double


!***********************************************************************************
!   GET INFORMATION   
!***********************************************************************************
    
    integer function eiram_get_id(fileName, reactionType, reactionIndex, err)

        character(*), intent(in) :: fileName
        character(*), intent(in) :: reactionType
        character(*), intent(in) :: reactionIndex
        integer, intent(out) :: err
 
        integer :: i       

        eiram_get_id = 0
        if(.not.allocated(reactions)) then
            write(eiram_unit, *) "ERROR in eiram_get_id: EIRAM is not initialized"
            err = 50
            return
        end if
        
        err = 0
        
        do i = 1, size(reactions)
            if(reactions(i)%FILE == filename &
              .and. reactions(i)%TYPE == reactiontype &
              .and. reactions(i)%INDEX == reactionindex) then
                eiram_get_id=i
                return
            end if
        end do
        
        eiram_get_id = 0
        write(eiram_unit, *) "ERROR in eiram_get_id: reaction "//&
            trim(filename)//" "//trim(reactiontype)//" "//trim(reactionindex)//&
            " not found"
        err = 110
    end function eiram_get_id

    subroutine eiram_get_order(N,M,Id,err)

     integer,intent(out) :: N
     integer,intent(out) :: M
     integer,intent(in) :: Id 
     integer,intent(out) :: err

     intrinsic size

     err=0

     if(.not.allocated(reactions)) then
         write(eiram_unit, *) "ERROR in eiram_get_order: EIRAM is not initialized"
         err = 50
         return
     end if

     if(id<0 .or. id>size(reactions) ) then 
       write(eiram_unit,*) "ERROR in eiram_get_order: reaction index is out of range"
       write(eiram_unit,*) " index, min, max ",id,1,size(reactions)
       err = 100
       return
     end if 

     n=reactions(id)%N
     m=reactions(id)%M

    end subroutine eiram_get_order

    subroutine eiram_create_image(data,err)

     type(eiram_data), allocatable :: data(:)
     integer,intent(out) :: err    

     integer :: k,L,M,N,err0

     intrinsic size,trim 

     err = 0

     if(.not.allocated(reactions)) then
        write(eiram_unit, *) "ERROR in eiram_create_image: EIRAM is not initialized"
        err = 50
        return
     end if

     l=size(reactions)   

     if(allocated(data)) then 
      do k=1,l
       if(allocated(data(k)%CREAC)) then
        deallocate(data(k)%CREAC,stat=err)
        if(err /=0) then 
         write(eiram_unit, *) "ERROR in eiram_create_image: cannot deallocate memory"
         err=200
         return
        end if  
       end if
      end do
      deallocate(data,stat=err)
      if(err /=0) then
       write(eiram_unit, *) "ERROR in eiram_create_image: cannot deallocate memory"
       err=200
       return
      end if
     end if 
     
     allocate(data(l),stat=err)
     if(err /= 0) then
       write(eiram_unit, *) "ERROR in eiram_create_image: cannot allocate memory for the array of reactions"
       write(eiram_unit, *) " Total number of reactions = ",l
       err = 200
       return
     end if

     do k=1,l
      if(.not.allocated(reactions(k)%CREAC)) goto 300
      n=reactions(k)%N
      m=reactions(k)%M
      data(k)%N=n
      data(k)%M=m      
      if(m>0) then
       allocate(data(k)%creac(0:m,0:n),stat=err)
      else 
       allocate(data(k)%creac(0:n,0:0),stat=err)
      end if
      if(err /= 0) goto 200
      data(k)%creac=reactions(k)%CREAC
      data(k)%X1MIN = reactions(k)%X1MIN
      data(k)%X1MAX = reactions(k)%X1MAX
      data(k)%X2MIN = reactions(k)%X2MIN
      data(k)%X2MAX = reactions(k)%X2MAX
      data(k)%X1UNITS =reactions(k)%X1UNITS 
      data(k)%X2UNITS =reactions(k)%X2UNITS
      data(k)%YUNITS = reactions(k)%YUNITS 
      data(k)%X1QUANTITY =reactions(k)%X1QUANTITY
      data(k)%X2QUANTITY =reactions(k)%X2QUANTITY
      data(k)%YQUANTITY = reactions(k)%YQUANTITY
      data(k)%REACTION =reactions(k)%REACTION
      data(k)%DESCRIPTION =reactions(k)%DESCRIPTION
      data(k)%FILE =reactions(k)%FILE 
      data(k)%TYPE =reactions(k)%TYPE 
      data(k)%INDEX =reactions(k)%INDEX 
     end do

     return

200  continue
     write(eiram_unit, *) "ERROR in eiram_create_image: error wile allocating memory"
     write(eiram_unit, *) " Reaction ID, File, Type, Index ",&
                            k,reactions(k)%FILE, reactions(k)%TYPE, reactions(k)%INDEX 
     write(eiram_unit, *) " N, M ",n,m       
     err = 200
     do k=1,l
      if(allocated(data(k)%CREAC)) deallocate(data(k)%CREAC,stat=err0)
     end do
     deallocate(data,stat=err0)
     return
300  continue
     write(eiram_unit, *) "ERROR detected in eiram_create_image: ", &
                          "array of coefficients is not allocated"
     write(eiram_unit, *) " Reaction ID, File, Type, Index ",&
                            k,reactions(k)%FILE,reactions(k)%TYPE,reactions(k)%INDEX 
     err = 50
     do k=1,l
      if(allocated(data(k)%CREAC)) deallocate(data(k)%CREAC,stat=err0)
     end do
     deallocate(data,stat=err0)
     return
    end subroutine eiram_create_image
     
!***********************************************************************************
!   CALCULATE FITS 
!***********************************************************************************

    real(EIRAM_DP) function eiram_calc1_fast(Id,lnX)

     integer, intent(in) :: Id
     real(EIRAM_DP), intent(in) :: lnX

     intrinsic exp

     eiram_calc1_fast = eiram_fit( reactions(id)%CREAC(:,0),reactions(id)%N,lnx)
     eiram_calc1_fast = exp(eiram_calc1_fast)

     end function eiram_calc1_fast

    subroutine eiram_calc1(Y,Id,lnX,err)

        real(EIRAM_DP), intent(out) :: Y(:) 
        integer, intent(in) :: Id
        real(EIRAM_DP), intent(in) :: lnX(:)
        integer, intent(out) :: err

        integer :: k,N,M,L
 
        intrinsic size,exp

        err = 0

        l=size(y)
        if(l /= size(lnx)) then
          write(eiram_unit,*) "ERROR in eiram_calc1: arrays have different lengths"
          write(eiram_unit,*) " length(Y), length(X) ",size(y),size(lnx)
          err = 100
          return
        end if

        if(.not.allocated(reactions)) then
            write(eiram_unit, *) "ERROR in eiram_calc1: EIRAM is not initialized"
            err = 50
            return
        end if

        m=size(reactions)
        if(id<0 .or. id>m ) then 
          write(eiram_unit,*) "ERROR in eiram_calc1: reaction index is out of range"
          write(eiram_unit,*) " index, min, max ",id,1,m
          err = 100
          return
        end if 

        if(.not.allocated(reactions(id)%CREAC)) then
          write(eiram_unit,*) "ERROR in eiram_calc1: array of fit coefficients is not allocated, Id ",id
          err = 50
          return
        end if

        m=reactions(id)%M

        if(m>0) then
          write(eiram_unit,*) "ERROR in eiram_calc1: wrong type of reaction"
          write(eiram_unit,*) " Single polynomial fit is requested for reaction ",&
                              " described by a double polynomial fit" 
          err = 130
          return    
        end if 

        n=reactions(id)%N    
        do k=1,l   
         y(k)=eiram_fit( reactions(id)%CREAC(:,0),n,lnx(k))
         y(k)=exp(y(k))
        end do

    end subroutine eiram_calc1

    subroutine eiram_calc12(B,Id,lnX,err)

     real(EIRAM_DP), intent(out) :: B(:,:) 
     integer, intent(in) :: Id
     real(EIRAM_DP), intent(in) :: lnX(:)
     integer, intent(out) :: err

     integer i,k,L,M,N

     intrinsic size

     err=0  

     l=size(b,2)
     if(l /= size(lnx)) then
       write(eiram_unit,*) "ERROR in eiram_calc12: arrays have different lengths"
       write(eiram_unit,*) " length(B), length(X) ",l,size(lnx)
       err = 100
       return
     end if

     if(.not.allocated(reactions)) then
         write(eiram_unit, *) "ERROR in eiram_calc12: EIRAM is not initialized"
         err = 50
         return
     end if

     m=size(reactions)
     if(id<0 .or. id>m ) then 
       write(eiram_unit,*) "ERROR in eiram_calc12: reaction index is out of range"
       write(eiram_unit,*) " index, min, max ",id,1,m
        err = 100
        return
     end if 

     if(.not.allocated(reactions(id)%CREAC)) then
      write(eiram_unit,*) "ERROR in eiram_calc12: array of fit coefficients is not allocated, Id ",id
      err = 50
      return
     end if

     m=reactions(id)%M
     if(m<1) then
       write(eiram_unit,*) "ERROR in eiram_calc12: wrong type of reaction"
       write(eiram_unit,*) " Double polynomial fit is requested for reaction ",&
                           " described by a single polynomial fit" 
       err = 120
      return    
     end if 
              
     n=reactions(id)%N
     if( n+1 /= size(b,1) ) then
       write(eiram_unit,*) "ERROR in eiram_calc12: mismatch of array sizes"
       write(eiram_unit,*) " order of polynomial  ",n," requested order ",size(b,1)-1
       err = 100
       return
     end if

     do k=1,l 
      do i=0,n
        b(i+1,k)=eiram_fit(reactions(id)%CREAC(:,i),m,lnx(k))
      end do 
     end do 

    end subroutine eiram_calc12

    subroutine eiram_calc2(Y,Id,LnX1,LnX2,err)

        real(EIRAM_DP), intent(out) :: Y(:) 
        integer, intent(in) :: Id
        real(EIRAM_DP), intent(in) :: lnX1(:)
        real(EIRAM_DP), intent(in) :: lnX2(:)
        integer, intent(out) :: err

        integer :: k,i,N,M,L,q,qk
       
        intrinsic size,exp,min

        err = 0

        l=size(y)
        if(l /= size(lnx2)) then
          write(eiram_unit,*) "ERROR in eiram_calc2: arrays have different lengths"
          write(eiram_unit,*) " length(Y), length(X2) ",size(y),size(lnx2)
          err = 100
          return
        end if
 
        q=size(lnx1) 
        if(q /=l .and. q /=1) then
          write(eiram_unit,*) "ERROR in eiram_calc2: arrays have wrong length"
          write(eiram_unit,*) " length(X1), length(Y),  ",size(lnx1),size(y)
          write(eiram_unit,*) " Length of X1 this array must be either = length(Y), or = 1"
          err = 100
          return
        end if 

        if(.not.allocated(reactions)) then
            write(eiram_unit, *) "ERROR in eiram_calc2: EIRAM is not initialized"
            err = 50
            return
        end if

        m=size(reactions)
        if(id<0 .or. id>m) then 
          write(eiram_unit,*) "ERROR in eiram_calc2: reaction index is out of range"
          write(eiram_unit,*) " index, min, max ",id,1,m
          err = 100
          return
        end if 

        if(.not.allocated(reactions(id)%CREAC)) then
         write(eiram_unit,*) "ERROR in eiram_calc 2: array of fit coefficients is not allocated, Id ",id
         err = 50
         return
        end if

        m=reactions(id)%M
        if(m<1) then
          write(eiram_unit,*) "ERROR in eiram_calc2: wrong type of reaction"
          write(eiram_unit,*) " Double polynomial fit is requested for reaction ",&
                              " described by a single polynomial fit" 
          err = 120
          return    
        end if 
            
       n=reactions(id)%N

       do k=1,l
        y(k)=eiram_fit(reactions(id)%CREAC(:,n),m,lnx2(k))
        qk=min(q,k)
        do i=n-1,0,-1
         y(k)= y(k)*lnx1(qk)+eiram_fit(reactions(id)%CREAC(:,i),m,lnx2(k))
        end do
        y(k)=exp(y(k))            
       end do      

    end subroutine eiram_calc2

!***********************************************************************************
!   AUXILIARY
!***********************************************************************************

    integer function countlinescontaining(filePath, string, err)

        character(*), intent(in) :: filePath
        character(*), intent(in) :: string
        integer, intent(out) :: err

        integer, parameter :: UNIT = 22
        character(MAX_LINE_LENGTH) :: line

        intrinsic trim

        open(unit, iostat=err,  file=filepath, status='old', action='read')
        countlinescontaining = 0
        if(err /= 0) then
            write(eiram_unit, *) "ERROR in countLinesContaining while reading "//trim(filepath), err
            close(unit)
            err = 300
            return
        end if
        do
          read(unit, '(A)', iostat=err,end=100) line
          if(err /= 0) then
             write(eiram_unit, *) "ERROR in countLinesContaining while reading line of "//trim(filepath)
            close(unit)
             err = 300
             return
          end if
          if(index(line, string) > 0) countlinescontaining = countlinescontaining + 1
        end do 
100     err = 0
        close(unit)
    end function countlinescontaining
    
end module eiram