EIRAM
atomic and molecular data in form of polynomial fits
eiram.f90
Go to the documentation of this file.
1 
6 
7 ! Copyright (c) 2016 Forschungszentrum Juelich GmbH
8 ! Markus Brenneis, Vladislav Kotov
9 !
10 ! This file is part of EIRAM.
11 !
12 ! EIRAM is free software: you can redistribute it and/or modify
13 ! it under the terms of the GNU General Public License as published by
14 ! the Free Software Foundation, either version 3 of the License, or
15 ! (at your option) any later version.
16 !
17 ! EIRAM is distributed in the hope that it will be useful,
18 ! but WITHOUT ANY WARRANTY; without even the implied warranty of
19 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 ! GNU General Public License for more details.
21 !
22 ! You should have received a copy of the GNU General Public License
23 ! along with EIRAM. If not, see <http://www.gnu.org/licenses/>.
24 !
25 
62 
76  function eiram_fit(B,N,X)
77  implicit none
78  integer,parameter :: DP=selected_real_kind(13)
79  real(DP) :: eiram_fit
81  real(DP),intent(in) :: B(0:n)
83  integer,intent(in) :: N
85  real(DP),intent(in) :: X
86 
87  integer :: i
88 
89  eiram_fit=b(n)
90  do i=n-1,0,-1
91  eiram_fit=eiram_fit*x+b(i)
92  end do
93 
94  end function eiram_fit
95 
96 module eiram
97  implicit none
98 
99  public eiram_load
100  public eiram_deallocate
101  public eiram_get_id
102  public eiram_get_order
103  public eiram_create_image
104  public eiram_calc1_fast
105  public eiram_calc1
106  public eiram_calc12
107  public eiram_calc2
108 
110  integer, parameter, public :: eiram_dp = selected_real_kind(13)
112  integer, public :: eiram_unit = 6
113 
115  integer,parameter :: units_length=32
117  integer,parameter :: quantity_length=64
119  integer,parameter :: reaction_length=128
121  integer,parameter :: description_length=256
123  integer, public, parameter :: eiram_file_name_length = 32
125  integer, public, parameter :: eiram_type_length = 8
127  integer, public, parameter :: eiram_index_length = 16
129  integer, parameter :: max_line_length = 256
130 
131  logical, parameter :: debug = .false.
132 
133  external eiram_fit
134  real(eiram_dp),public :: eiram_fit
135 
142  type, public:: eiram_data
144  integer :: n
146  integer :: m
151  real(EIRAM_DP), allocatable :: creac(:,:)
153  real(EIRAM_DP) :: x1min
155  real(EIRAM_DP) :: x1max
157  real(EIRAM_DP) :: x2min
159  real(EIRAM_DP) :: x2max
161  character(UNITS_LENGTH) :: x1units = ""
163  character(UNITS_LENGTH) :: x2units = ""
165  character(UNITS_LENGTH) :: yunits = ""
167  character(QUANTITY_LENGTH) :: x1quantity = ""
169  character(QUANTITY_LENGTH) :: x2quantity = ""
171  character(QUANTITY_LENGTH) :: yquantity = ""
173  character(REACTION_LENGTH) :: reaction = ""
175  character(DESCRIPTION_LENGTH) :: description = ""
177  character(EIRAM_FILE_NAME_LENGTH) :: file = ""
179  character(EIRAM_TYPE_LENGTH) :: TYPE = ""
181  character(EIRAM_INDEX_LENGTH) :: index = ""
182  end type eiram_data
183 
184  private
185 
187  integer,public :: eiram_nreact=0
188 
190  type(eiram_data), allocatable :: reactions(:)
191 
192  contains
193 
194 !***********************************************************************************
195 ! INITIALIZATION/DEINITIALIZATION ROUTINES
196 !***********************************************************************************
197 
205  subroutine eiram_load(filePath, fileName, err)
207  character(*), intent(in) :: filePath
209  character(*), intent(in) :: fileName
211  integer, intent(out) :: err
212 
213  logical :: fileExists,isLaTeX
215  integer :: nData, nOldData
217  integer :: l,nDataMax
218  integer, parameter :: UNIT = 20
219  character(MAX_LINE_LENGTH) :: line
220  character(*), parameter :: BEGIN_DATA_MARK = "##BEGIN DATA HERE##"
221  !string which is used to find out that the document is in LaTeX format
222  ! - w/o '\' at the beginning because for some reason this confuses g95 compiler
223  character(*), parameter :: LATEX_MARK = "begin{document}"
224  type(eiram_data), allocatable:: reactions_(:)
225 
226  intrinsic index,trim
227 
228  if(debug) print *, "open file ",trim(filename)
229 
230  inquire(file=trim(filepath)//trim(filename),exist=fileexists,iostat=err)
231  if(.not.fileexists .or. err /= 0) then
232  err=310
233  write(eiram_unit, *) "ERROR in eiram_load: input file ", &
234  trim(filepath)//trim(filename)," is not found"
235  return
236  end if
237 
238  if(allocated(reactions)) then
239  nolddata = size(reactions)
240  else
241  nolddata = 0
242  end if
243 
244  ndatamax = countlinescontaining(trim(filepath)//trim(filename), "Reaction", err)
245  if(err/=0) then
246  write(eiram_unit, *) "ERROR in eiram_load while counting lines, err ", err
247  err = 300
248  return
249  end if
250 
251  allocate(reactions_(ndatamax+nolddata),stat=err)
252  if(err/=0) then
253  write(eiram_unit, *) "ERROR in eiram_load while allocating temporary array, ",&
254  "nData",ndatamax+nolddata
255  err = 200
256  return
257  end if
258 
259  open(unit, iostat=err, file=trim(filepath)//trim(filename), status='old', action='read')
260  if(err /= 0) then
261  write(eiram_unit, *) "ERROR in eiram_load while reading "//trim(filepath)//trim(filename)
262  close(unit)
263  err = 300
264  return
265  end if
266 
267  l = 0
268  ndata = 0
269  islatex=.false.
270 
271  !look for beginning of the data section
272  do
273  read(unit, '(A)', iostat=err,end=1000) line
274  l = l + 1
275  if(err /= 0) goto 1300
276 
277  if(index(line,latex_mark) > 0) islatex=.true.
278 
279  if(index(line, begin_data_mark) > 0) then
280  if(debug) print *, "found BEGIN_DATA_MARK in line ", l
281  exit
282  end if
283  end do
284 
285  !read the data
286  if(islatex) then
287  call eiram_read_latex(reactions_,filename,unit,l,ndata,*1300,*1303,*1304)
288  else
289  call eiram_read_txt(reactions_,filename,unit,l,ndata,*1300,*1303,*1304)
290  end if
291 
292  if(debug) print *, "EOF in", l
293  close(unit)
294 
295  if(allocated(reactions)) then
296  !Add reactions which were already loaded in previous calls of eiram_load
297  reactions_(ndata+1:ndata+nolddata) = reactions(1:nolddata)
298  deallocate(reactions, stat=err)
299  if(err /= 0) then
300  write(eiram_unit, *) "ERROR in eiram_load: cannot re-allocate memory ", err
301  err = 200
302  return
303  end if
304  end if
305 
306  !(Re-)allocate REACTIONS
307  allocate(reactions(ndata+nolddata),stat=err)
308  if(err /= 0) then
309  write(eiram_unit, *) "ERROR in eiram_load: cannot allocate the array of reactions ",&
310  "nData ",ndata+nolddata
311  err = 200
312  return
313  end if
314 
315  reactions(1:ndata+nolddata) = reactions_(1:ndata+nolddata)
316 
317  call deallocate_reactions
318 
319  write(eiram_unit, *) "initalization from file "//trim(filepath)//trim(filename)//" completed"
320  write(eiram_unit, *) " new data sets loaded: ", ndata
321  eiram_nreact=ndata+nolddata
322  write(eiram_unit, *) " total number of data sets: ", eiram_nreact
323  return
324 
325 1000 write(eiram_unit, *) "ERROR in eiram_load: no data section is found in file "//trim(filepath)//trim(filename)
326  close(unit)
327  call deallocate_reactions
328  err = 300
329  return
330 1300 write(eiram_unit, *) "ERROR in eiram_load: while reading line ", l, " of "//trim(filepath)//trim(filename)
331  close(unit)
332  call deallocate_reactions
333  err = 300
334  return
335 1303 write(eiram_unit, *) "ERROR in eiram_load: wrong format in line ", l, &
336  " of "//trim(filepath)//trim(filename), err
337  close(unit)
338  call deallocate_reactions
339  err = 300
340  return
341 1304 write(eiram_unit, *) "ERROR in eiram_load: cannot allocate CREAC "
342  write(eiram_unit, *) " line ", l, " of ",trim(filepath)//trim(filename),&
343  " reaction ",trim(reactions_(ndata)%TYPE)//" "//trim(reactions_(ndata)%INDEX)
344  close(unit)
345  call deallocate_reactions
346  err = 200
347  return
348 
349  contains
350 
351  !clean-up for the case if the compiler uses argument -save
352  subroutine deallocate_reactions
353  integer :: i
354  logical :: lerr
355  lerr=.false.
356  do i=1,size(reactions_)
357  if(allocated(reactions_(i)%CREAC)) then
358  deallocate(reactions_(i)%CREAC,stat=err)
359  if(err /= 0) lerr=.true.
360  end if
361  end do
362  deallocate(reactions_,stat=err)
363  if(err /= 0 .or. lerr ) then
364  write(eiram_unit, *) "WARNING from eiram_load: problem while deallocating temporary array "
365  err = -200
366  end if
367  end subroutine deallocate_reactions
368 
369  end subroutine eiram_load
370 
374  subroutine eiram_read_latex(reactions_,fileName,unit,l,nData,*,*,*)
375  type(eiram_data), allocatable:: reactions_(:)
376  character(*), intent(in) :: fileName
377  integer,intent(in) :: unit
378  integer,intent(inout) :: l
379  integer,intent(out) :: nData
380  integer :: gotoLine,numberOfUnclosedBraces,n,err,i1,i2
381  character(EIRAM_TYPE_LENGTH) :: stype
382  character(MAX_LINE_LENGTH) :: line
383  character(MAX_LINE_LENGTH*2) :: reactionName
384 
385  intrinsic index,count,transfer,trim,adjustl,len_trim
386 
387  ndata=0
388 
389  foreachline: do
390  read(unit, '(A)', iostat=err,end=100) line
391  l = l + 1
392  if(err /= 0) return 1
393  newsection: if(index(line, "\section") == 1) then
394  if(debug) print *, "found \section in line ", l
395  if(index(trim(line), " ") == 0) then
396  read(unit, '(A)', iostat=err) line
397  l = l + 1
398  if(err /= 0) return 1
399  end if
400  stype= line(index(line,"{")+1:index(line," ")-1)
401  if(index(stype,":")>0) stype=stype(1:index(stype,":")-1)
402  end if newsection
403 
404  newreaction: if(index(line, "subsection{") > 0) then
405  numberofunclosedbraces = count(transfer(line, 'a', max_line_length) == "{")
406  n = count(transfer(line, 'a', max_line_length) == "}")
407  numberofunclosedbraces = numberofunclosedbraces - n
408  reactionname = trim(adjustl(line))
409 
410  !for the case when the reaction is defined like this
411  !\subsection{
412  ! ............ }
413  findendofreactionname: do
414  if(numberofunclosedbraces == 0) exit findendofreactionname
415  read(unit, '(A)', iostat=err) line
416  l = l + 1
417  if(err /= 0) return 1
418  n = count(transfer(line, 'a', max_line_length) == "{")
419  numberofunclosedbraces = numberofunclosedbraces + n
420  n = count(transfer(line, 'a', max_line_length) == "}")
421  numberofunclosedbraces = numberofunclosedbraces - n
422  if(len_trim(reactionname) /= index(reactionname, "{")) then
423  reactionname = trim(reactionname)//" "//trim(adjustl(line))
424  else
425  reactionname = trim(reactionname)//trim(adjustl(line))
426  end if
427  if(line == "}") exit findendofreactionname
428  end do findendofreactionname
429 
430  if(index(reactionname,'Reaction')<1) cycle
431 
432  ndata = ndata+1
433 
434  reactions_(ndata)%FILE = filename
435  reactions_(ndata)%TYPE = stype
436 
437  if(debug) print *, "found Reaction in line ", l
438  !Sub-string between 1st and 2nd space after "Reaction": the index
439  i1=index(reactionname,"Reaction")+8;
440  reactions_(ndata)%INDEX = reactionname(index(reactionname(i1:)," ")+i1:)
441  reactions_(ndata)%INDEX = reactions_(ndata)%INDEX(:index(reactions_(ndata)%INDEX," ")-1)
442  if(debug) print *, " with index ", trim(reactions_(ndata)%INDEX)
443  !Sub-string between $ and $ space: reaction equation (if available)
444  i1=index(reactionname,"$")
445  i2=index(reactionname,"$",.true.)
446  if (i1>0 .and. i2>0 ) then
447  reactions_(ndata)%REACTION = reactionname(i1:i2)
448  if(debug) print *, " with equation ", trim(reactions_(ndata)%REACTION)
449  end if
450 
451  skiptillverbatim: do
452  read(unit, '(A)', iostat=err) line
453  l = l + 1
454  if(err /= 0) return 1
455  if(index(line, "verbatim") > 0) then
456  exit skiptillverbatim
457  end if
458  end do skiptillverbatim
459 
460  select case(reactions_(ndata)%TYPE)
461  case("H.1", "H.2", "H.8", "H.11")
462  if(debug) print *, "h1/2/8/11"
463  call parse_single(reactions_(ndata),unit,gotoline,l)
464  if(gotoline == 1300) return 1
465  if(gotoline == 1303) return 2
466  if(gotoline == 1304) return 3
467  if(.not.allocated(reactions_(ndata)%CREAC)) ndata=ndata-1 !skip this reaction
468  case ("H.3", "H.4", "H.9", "H.10", "H.12")
469  if(debug) print *, "h3/4/9/10/12"
470  call parse_double(reactions_(ndata),unit,gotoline,l)
471  if(gotoline == 1300) return 1
472  if(gotoline == 1303) return 2
473  if(gotoline == 1304) return 3
474  if(.not.allocated(reactions_(ndata)%CREAC)) ndata=ndata-1 !skip this reaction
475  case default
476  ! write(eiram_unit, *) "WARNING from eiram_load: No handle for type "&
477  ! //trim(reactions_(nData)%TYPE)//" in line", l
478  if(debug) print *, ">> exiting"
479  ndata = ndata-1 !skip reaction
480  end select
481  end if newreaction
482 
483  end do foreachline
484 
485  100 continue
486 
487  end subroutine eiram_read_latex
488 
492  subroutine eiram_read_txt(reactions_,fileName,unit,l,nData,*,*,*)
493  type(eiram_data), allocatable:: reactions_(:)
494  character(*), intent(in) :: fileName
495  integer,intent(in) :: unit
496  integer,intent(inout) :: l
497  integer,intent(out) :: nData
498  integer :: gotoLine,err,i1,i2,itype
499  character(EIRAM_TYPE_LENGTH) :: stype
500  character(MAX_LINE_LENGTH) :: line, reactionName
501 
502  intrinsic index,trim,adjustl,len_trim
503 
504  ndata=0
505  itype=0
506 
507  foreachline: do
508  read(unit, '(A)', iostat=err,end=100) line
509  l = l + 1
510  if(err /= 0) return 1
511  call gettype(line,stype)
512 
513  newreaction: if(index(line, "Reaction") > 0) then
514  ndata = ndata+1
515  reactions_(ndata)%FILE = filename
516  reactions_(ndata)%TYPE = stype
517  reactionname = trim(adjustl(line))
518  !Sub-string between 1st and 2nd space: the index
519  i1=index(reactionname," ")+1
520  i2=index(reactionname(i1:)," ")+i1-2
521  reactions_(ndata)%INDEX = reactionname(i1:i2)
522  !Sub-string between 2nd space and end of the line: reaction equation
523  reactions_(ndata)%REACTION =trim(adjustl(reactionname(i2+1:)))
524  if(debug) then
525  print *, "found new reaction"
526  print *, " index ",trim(reactions_(ndata)%INDEX)
527  print *, " equation ",trim(reactions_(ndata)%REACTION)
528  end if
529  select case(reactions_(ndata)%TYPE)
530  case("H.1", "H.2", "H.8", "H.11")
531  if(debug) print *, "h1/2/8/11"
532  call parse_single(reactions_(ndata),unit,gotoline,l)
533  if(gotoline == 1300) return 1
534  if(gotoline == 1303) return 2
535  if(gotoline == 1304) return 3
536  if(.not.allocated(reactions_(ndata)%CREAC)) ndata=ndata-1 !skip this reaction
537  case ("H.3", "H.4", "H.9", "H.10", "H.12")
538  if(debug) print *, "h3/4/9/10/12"
539  call parse_double(reactions_(ndata),unit,gotoline,l)
540  if(gotoline == 1300) return 1
541  if(gotoline == 1303) return 2
542  if(gotoline == 1304) return 3
543  if(.not.allocated(reactions_(ndata)%CREAC)) ndata=ndata-1 !skip this reaction
544  case default
545  ! write(eiram_unit, *) "WARNING from eiram_load: No handle for type "&
546  ! //trim(reactions_(nData)%TYPE)//" in line", l
547  if(debug) print *, ">> exiting"
548  ndata = ndata-1 !skip reaction
549  end select
550  end if newreaction
551  end do foreachline
552 
553 100 continue
554 
555  contains
556 
558  subroutine gettype(line,stype)
559  character(MAX_LINE_LENGTH),intent(in) :: line
560  character(EIRAM_TYPE_LENGTH),intent(inout) :: stype
561 
562  integer,parameter :: ntp=12
563  integer :: i
564  character(4),parameter :: stp(0:ntp)= &
565  (/'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'/)
566 
567  if(index(line,'H.')<1) return !for speed up
568 
569  do i=ntp,itype,-1
570  if(index(line,trim(stp(i)))>0) then
571  stype=trim(stp(i))
572  itype=i
573  if(debug) print *, "found new type ", stype
574  return
575  end if
576  end do
577 
578  end subroutine gettype
579 
580  end subroutine eiram_read_txt
581 
582 
584  subroutine eiram_deallocate(err)
586  integer, intent(out) :: err
587 
588  logical :: lerr
589  integer i
590 
591  intrinsic size
592 
593  err = 0
594 
595  eiram_nreact=0
596 
597  if(allocated(reactions)) then
598 
599  lerr=.false.
600  do i=1,size(reactions)
601  if(allocated(reactions(i)%CREAC)) then
602  deallocate(reactions(i)%CREAC,stat=err)
603  if(err /= 0) lerr=.true.
604  end if
605  end do
606 
607  deallocate(reactions, stat=err)
608  if(err /= 0) then
609  write(eiram_unit, *) "WARNING from eiram_deallocate: problem while deallocating memory, err ", err
610  err = -200
611  return
612  end if
613 
614  if(lerr) then
615  write(eiram_unit, *) "WARNING from eiram_deallocate: problem while deallocating memory"
616  err=-200
617  end if
618 
619  end if !if(allocated(allocated(REACTIONS)) then
620 
621  end subroutine eiram_deallocate
622 
624  subroutine parse_single(reaction,unit,gotoLine,l)
625  type(eiram_data) :: reaction
626  integer,intent(in) :: unit
627  integer,intent(out) :: gotoLine
628  integer,intent(inout) :: l
629 
630  character(2) :: chr
631  character(80) :: zeile
632  integer :: err,indff,iftflg,ind,ind2,i,j
633 
634  intrinsic index,mod
635 
636  gotoline = 0
637 
638  reaction%M = 0
639  reaction%N = 8
640  if(reaction%TYPE == "H.1") then
641  reaction%YUNITS = "cm^2"
642  reaction%YQUANTITY = "cross section"
643  reaction%X1QUANTITY = "Energy"
644  chr='a0'
645  elseif(reaction%TYPE == "H.2") then
646  reaction%YUNITS = "cm^3/s"
647  reaction%YQUANTITY = "collision rate"
648  reaction%X1QUANTITY = "Temperature"
649  chr='b0'
650  elseif(reaction%TYPE == "H.8") then
651  reaction%YUNITS = "cm^3/s\cdot eV"
652  reaction%YQUANTITY = "electron cooling rate"
653  reaction%X1QUANTITY = "Temperature"
654  chr='h0'
655  elseif(reaction%TYPE == "H.11") then
656  reaction%YUNITS = ""
657  reaction%YQUANTITY = "ratio"
658  reaction%X1QUANTITY = "Temperature"
659  chr='k0'
660  end if
661  reaction%X1UNITS = "eV"
662 
663  !this code mimics EIRENE, slreac.f
664 3 read(unit,'(A80)',iostat=err) zeile
665  l=l+1
666  if(err /= 0) then
667  gotoline = 1300
668  return
669  end if
670 
671  iftflg=0
672 
673  !Skip "An analytic formula is given in the text" in LaTeX mode
674  if( index(zeile,'verbatim')>0 .or. index(zeile,'subsection')>0 ) return
675 
676  indff=index(zeile,'fit-flag')
677  if (index(zeile,chr)+indff.EQ.0) GOTO 3
678  if(indff > 0) THEN
679  read (zeile((indff+8):80),*,iostat=err) iftflg
680  if(err /= 0) then
681  gotoline = 1303
682  return
683  end if
684  goto 3
685  end if
686 
687  allocate(reaction%CREAC(0:8,0:0),stat=err)
688  if(err /=0) then
689  gotoline=1304
690  return
691  end if
692 
693  if(mod(iftflg,100) == 10) then
694  ind=index(zeile,chr(1:1))
695  reaction%creac(:,0)=0.
696  read(zeile((ind+2):80),*,iostat=err) reaction%creac(0,0)
697  if(err /= 0) then
698  gotoline = 1303
699  return
700  end if
701  else
702  do j=0,2
703  ind=index(zeile,chr(1:1))+2
704  do i=0,2
705  ind2=index(zeile(ind:80),chr(1:1))
706  if(ind2<1) then
707  ind2=80
708  else
709  ind2=ind2+ind-2
710  end if
711  read(zeile(ind:ind2),*,iostat=err) reaction%creac(j*3+i,0)
712  if(err /= 0) then
713  gotoline = 1303
714  return
715  end if
716  ind=ind2+3
717  end do
718  read(unit,'(A80)',iostat=err) zeile
719  l=l+1
720  if(err /= 0) then
721  gotoline = 1300
722  return
723  end if
724  end do
725  end if
726 
727  if(debug) print *, reaction%CREAC
728 
729  end subroutine parse_single
730 
734  subroutine parse_double(reaction,unit,gotoLine,l)
735  type(eiram_data) :: reaction
736  integer,intent(in) :: unit
737  integer,intent(out) :: gotoLine
738  integer,intent(inout) :: l
739 
740  character(80) :: zeile
741  integer :: err,indff,iftflg,ih,i,j,k
742 
743  intrinsic index,mod
744 
745  gotoline = 0
746 
747  reaction%N = 8
748  reaction%M = 8
749  if(reaction%TYPE == "H.3") then
750  reaction%YQUANTITY = "collision rate"
751  reaction%YUNITS = "cm^3/s"
752  reaction%X1QUANTITY = "Energy"
753  reaction%X1UNITS = "eV"
754  reaction%X2QUANTITY = "Temperature"
755  reaction%X2UNITS = "eV"
756  reaction%X2QUANTITY = "Temperature"
757  elseif(reaction%TYPE == "H.4") then
758  reaction%YQUANTITY = "collision rate"
759  reaction%YUNITS = "cm^3/s"
760  reaction%X1QUANTITY = "Density"
761  reaction%X1UNITS = "cm^{-3}"
762  reaction%X2QUANTITY = "Temperature"
763  reaction%X2UNITS = "eV"
764  elseif(reaction%TYPE == "H.9") then
765  reaction%YQUANTITY = "ion cooling rate"
766  reaction%YUNITS = "cm^3/s\cdot eV"
767  reaction%X1QUANTITY = "Energy"
768  reaction%X1UNITS = "[eV]"
769  reaction%X2QUANTITY = "Temperature"
770  reaction%X2UNITS = "eV"
771  elseif(reaction%TYPE == "H.10") then
772  reaction%YQUANTITY = "electron cooling rate"
773  reaction%YUNITS = "cm^3/s\cdot eV"
774  reaction%X1QUANTITY = "Density"
775  reaction%X1UNITS = "cm^{-3}"
776  reaction%X2QUANTITY = "Temperature"
777  reaction%X2UNITS = "eV"
778  elseif(reaction%TYPE == "H.12") then
779  reaction%YQUANTITY = "ratio"
780  reaction%YUNITS = ""
781  reaction%X1QUANTITY = "Density"
782  reaction%X1UNITS = "cm^{-3}"
783  reaction%X2QUANTITY = "Temperature"
784  reaction%X2UNITS = "eV"
785  end if
786 
787  allocate(reaction%CREAC(0:8,0:8),stat=err)
788  if(err /=0) then
789  gotoline=1304
790  return
791  end if
792 
793  iftflg=0
794 
795  !this code mimics EIRENE, slreac.f
796  do j=0,2
797 16 read(unit,'(A80)',iostat=err) zeile
798  l=l+1
799  if(err /= 0) then
800  gotoline = 1300
801  return
802  end if
803 
804  !Skip "An analytic formula is given in the text" in LaTeX mode
805  if( index(zeile,'verbatim')>0 .or. index(zeile,'subsection')>0 ) return
806 
807  indff=index(zeile,'fit-flag')
808  if(index(zeile,'Index')+indff==0) goto 16
809  if(indff > 0) then
810  read(zeile((indff+8):80),*,iostat=err) iftflg
811  if(err /= 0) then
812  gotoline = 1303
813  return
814  end if
815  goto 16
816  end if
817 
818  read(unit,'(1X)',iostat=err)
819  l=l+1
820  if(err /= 0) then
821  gotoline = 1300
822  return
823  end if
824 
825  if(mod(iftflg,100) == 10) then
826  read(unit,*,iostat=err) ih,reaction%CREAC(0,0)
827  l=l+1
828  if(err /= 0) then
829  gotoline = 1300
830  return
831  end if
832  exit
833  else
834  do i=0,8
835  read(unit,*,iostat=err) ih,(reaction%CREAC(i,k),k=j*3,j*3+2)
836  if(err /= 0) then
837  gotoline = 1300
838  return
839  end if
840  end do
841  end if
842  end do
843 
844  if(debug) print *, reaction%CREAC
845 
846  end subroutine parse_double
847 
848 
849 !***********************************************************************************
850 ! GET INFORMATION
851 !***********************************************************************************
852 
856  integer function eiram_get_id(fileName, reactionType, reactionIndex, err)
858  character(*), intent(in) :: fileName
860  character(*), intent(in) :: reactionType
862  character(*), intent(in) :: reactionIndex
864  integer, intent(out) :: err
865 
866  integer :: i
867 
868  eiram_get_id = 0
869  if(.not.allocated(reactions)) then
870  write(eiram_unit, *) "ERROR in eiram_get_id: EIRAM is not initialized"
871  err = 50
872  return
873  end if
874 
875  err = 0
876 
877  do i = 1, size(reactions)
878  if(reactions(i)%FILE == filename &
879  .and. reactions(i)%TYPE == reactiontype &
880  .and. reactions(i)%INDEX == reactionindex) then
881  eiram_get_id=i
882  return
883  end if
884  end do
885 
886  eiram_get_id = 0
887  write(eiram_unit, *) "ERROR in eiram_get_id: reaction "//&
888  trim(filename)//" "//trim(reactiontype)//" "//trim(reactionindex)//&
889  " not found"
890  err = 110
891  end function eiram_get_id
892 
894  subroutine eiram_get_order(N,M,Id,err)
896  integer,intent(out) :: N
898  integer,intent(out) :: M
900  integer,intent(in) :: Id
902  integer,intent(out) :: err
903 
904  intrinsic size
905 
906  err=0
907 
908  if(.not.allocated(reactions)) then
909  write(eiram_unit, *) "ERROR in eiram_get_order: EIRAM is not initialized"
910  err = 50
911  return
912  end if
913 
914  if(id<0 .or. id>size(reactions) ) then
915  write(eiram_unit,*) "ERROR in eiram_get_order: reaction index is out of range"
916  write(eiram_unit,*) " index, min, max ",id,1,size(reactions)
917  err = 100
918  return
919  end if
920 
921  n=reactions(id)%N
922  m=reactions(id)%M
923 
924  end subroutine eiram_get_order
925 
930  subroutine eiram_create_image(data,err)
932  type(eiram_data), allocatable :: data(:)
934  integer,intent(out) :: err
935 
936  integer :: k,L,M,N,err0
937 
938  intrinsic size,trim
939 
940  err = 0
941 
942  if(.not.allocated(reactions)) then
943  write(eiram_unit, *) "ERROR in eiram_create_image: EIRAM is not initialized"
944  err = 50
945  return
946  end if
947 
948  l=size(reactions)
949 
950  if(allocated(data)) then
951  do k=1,l
952  if(allocated(data(k)%CREAC)) then
953  deallocate(data(k)%CREAC,stat=err)
954  if(err /=0) then
955  write(eiram_unit, *) "ERROR in eiram_create_image: cannot deallocate memory"
956  err=200
957  return
958  end if
959  end if
960  end do
961  deallocate(data,stat=err)
962  if(err /=0) then
963  write(eiram_unit, *) "ERROR in eiram_create_image: cannot deallocate memory"
964  err=200
965  return
966  end if
967  end if
968 
969  allocate(data(l),stat=err)
970  if(err /= 0) then
971  write(eiram_unit, *) "ERROR in eiram_create_image: cannot allocate memory for the array of reactions"
972  write(eiram_unit, *) " Total number of reactions = ",l
973  err = 200
974  return
975  end if
976 
977  do k=1,l
978  if(.not.allocated(reactions(k)%CREAC)) goto 300
979  n=reactions(k)%N
980  m=reactions(k)%M
981  data(k)%N=n
982  data(k)%M=m
983  if(m>0) then
984  allocate(data(k)%creac(0:m,0:n),stat=err)
985  else
986  allocate(data(k)%creac(0:n,0:0),stat=err)
987  end if
988  if(err /= 0) goto 200
989  data(k)%creac=reactions(k)%CREAC
990  data(k)%X1MIN = reactions(k)%X1MIN
991  data(k)%X1MAX = reactions(k)%X1MAX
992  data(k)%X2MIN = reactions(k)%X2MIN
993  data(k)%X2MAX = reactions(k)%X2MAX
994  data(k)%X1UNITS =reactions(k)%X1UNITS
995  data(k)%X2UNITS =reactions(k)%X2UNITS
996  data(k)%YUNITS = reactions(k)%YUNITS
997  data(k)%X1QUANTITY =reactions(k)%X1QUANTITY
998  data(k)%X2QUANTITY =reactions(k)%X2QUANTITY
999  data(k)%YQUANTITY = reactions(k)%YQUANTITY
1000  data(k)%REACTION =reactions(k)%REACTION
1001  data(k)%DESCRIPTION =reactions(k)%DESCRIPTION
1002  data(k)%FILE =reactions(k)%FILE
1003  data(k)%TYPE =reactions(k)%TYPE
1004  data(k)%INDEX =reactions(k)%INDEX
1005  end do
1006 
1007  return
1008 
1009 200 continue
1010  write(eiram_unit, *) "ERROR in eiram_create_image: error wile allocating memory"
1011  write(eiram_unit, *) " Reaction ID, File, Type, Index ",&
1012  k,reactions(k)%FILE, reactions(k)%TYPE, reactions(k)%INDEX
1013  write(eiram_unit, *) " N, M ",n,m
1014  err = 200
1015  do k=1,l
1016  if(allocated(data(k)%CREAC)) deallocate(data(k)%CREAC,stat=err0)
1017  end do
1018  deallocate(data,stat=err0)
1019  return
1020 300 continue
1021  write(eiram_unit, *) "ERROR detected in eiram_create_image: ", &
1022  "array of coefficients is not allocated"
1023  write(eiram_unit, *) " Reaction ID, File, Type, Index ",&
1024  k,reactions(k)%FILE,reactions(k)%TYPE,reactions(k)%INDEX
1025  err = 50
1026  do k=1,l
1027  if(allocated(data(k)%CREAC)) deallocate(data(k)%CREAC,stat=err0)
1028  end do
1029  deallocate(data,stat=err0)
1030  return
1031  end subroutine eiram_create_image
1032 
1033 !***********************************************************************************
1034 ! CALCULATE FITS
1035 !***********************************************************************************
1036 
1041  real(EIRAM_DP) function eiram_calc1_fast(Id,lnX)
1043  integer, intent(in) :: Id
1045  real(EIRAM_DP), intent(in) :: lnX
1046 
1047  intrinsic exp
1048 
1049  eiram_calc1_fast = eiram_fit( reactions(id)%CREAC(:,0),reactions(id)%N,lnx)
1051 
1052  end function eiram_calc1_fast
1053 
1059  subroutine eiram_calc1(Y,Id,lnX,err)
1061  real(EIRAM_DP), intent(out) :: Y(:)
1063  integer, intent(in) :: Id
1065  real(EIRAM_DP), intent(in) :: lnX(:)
1067  integer, intent(out) :: err
1068 
1069  integer :: k,N,M,L
1070 
1071  intrinsic size,exp
1072 
1073  err = 0
1074 
1075  l=size(y)
1076  if(l /= size(lnx)) then
1077  write(eiram_unit,*) "ERROR in eiram_calc1: arrays have different lengths"
1078  write(eiram_unit,*) " length(Y), length(X) ",size(y),size(lnx)
1079  err = 100
1080  return
1081  end if
1082 
1083  if(.not.allocated(reactions)) then
1084  write(eiram_unit, *) "ERROR in eiram_calc1: EIRAM is not initialized"
1085  err = 50
1086  return
1087  end if
1088 
1089  m=size(reactions)
1090  if(id<0 .or. id>m ) then
1091  write(eiram_unit,*) "ERROR in eiram_calc1: reaction index is out of range"
1092  write(eiram_unit,*) " index, min, max ",id,1,m
1093  err = 100
1094  return
1095  end if
1096 
1097  if(.not.allocated(reactions(id)%CREAC)) then
1098  write(eiram_unit,*) "ERROR in eiram_calc1: array of fit coefficients is not allocated, Id ",id
1099  err = 50
1100  return
1101  end if
1102 
1103  m=reactions(id)%M
1104 
1105  if(m>0) then
1106  write(eiram_unit,*) "ERROR in eiram_calc1: wrong type of reaction"
1107  write(eiram_unit,*) " Single polynomial fit is requested for reaction ",&
1108  " described by a double polynomial fit"
1109  err = 130
1110  return
1111  end if
1112 
1113  n=reactions(id)%N
1114  do k=1,l
1115  y(k)=eiram_fit( reactions(id)%CREAC(:,0),n,lnx(k))
1116  y(k)=exp(y(k))
1117  end do
1118 
1119  end subroutine eiram_calc1
1120 
1141  subroutine eiram_calc12(B,Id,lnX,err)
1145  real(EIRAM_DP), intent(out) :: B(:,:)
1147  integer, intent(in) :: Id
1149  real(EIRAM_DP), intent(in) :: lnX(:)
1151  integer, intent(out) :: err
1152 
1153  integer i,k,L,M,N
1154 
1155  intrinsic size
1156 
1157  err=0
1158 
1159  l=size(b,2)
1160  if(l /= size(lnx)) then
1161  write(eiram_unit,*) "ERROR in eiram_calc12: arrays have different lengths"
1162  write(eiram_unit,*) " length(B), length(X) ",l,size(lnx)
1163  err = 100
1164  return
1165  end if
1166 
1167  if(.not.allocated(reactions)) then
1168  write(eiram_unit, *) "ERROR in eiram_calc12: EIRAM is not initialized"
1169  err = 50
1170  return
1171  end if
1172 
1173  m=size(reactions)
1174  if(id<0 .or. id>m ) then
1175  write(eiram_unit,*) "ERROR in eiram_calc12: reaction index is out of range"
1176  write(eiram_unit,*) " index, min, max ",id,1,m
1177  err = 100
1178  return
1179  end if
1180 
1181  if(.not.allocated(reactions(id)%CREAC)) then
1182  write(eiram_unit,*) "ERROR in eiram_calc12: array of fit coefficients is not allocated, Id ",id
1183  err = 50
1184  return
1185  end if
1186 
1187  m=reactions(id)%M
1188  if(m<1) then
1189  write(eiram_unit,*) "ERROR in eiram_calc12: wrong type of reaction"
1190  write(eiram_unit,*) " Double polynomial fit is requested for reaction ",&
1191  " described by a single polynomial fit"
1192  err = 120
1193  return
1194  end if
1195 
1196  n=reactions(id)%N
1197  if( n+1 /= size(b,1) ) then
1198  write(eiram_unit,*) "ERROR in eiram_calc12: mismatch of array sizes"
1199  write(eiram_unit,*) " order of polynomial ",n," requested order ",size(b,1)-1
1200  err = 100
1201  return
1202  end if
1203 
1204  do k=1,l
1205  do i=0,n
1206  b(i+1,k)=eiram_fit(reactions(id)%CREAC(:,i),m,lnx(k))
1207  end do
1208  end do
1209 
1210  end subroutine eiram_calc12
1211 
1218  subroutine eiram_calc2(Y,Id,LnX1,LnX2,err)
1220  real(EIRAM_DP), intent(out) :: Y(:)
1222  integer, intent(in) :: Id
1224  real(EIRAM_DP), intent(in) :: lnX1(:)
1226  real(EIRAM_DP), intent(in) :: lnX2(:)
1228  integer, intent(out) :: err
1229 
1230  integer :: k,i,N,M,L,q,qk
1231 
1232  intrinsic size,exp,min
1233 
1234  err = 0
1235 
1236  l=size(y)
1237  if(l /= size(lnx2)) then
1238  write(eiram_unit,*) "ERROR in eiram_calc2: arrays have different lengths"
1239  write(eiram_unit,*) " length(Y), length(X2) ",size(y),size(lnx2)
1240  err = 100
1241  return
1242  end if
1243 
1244  q=size(lnx1)
1245  if(q /=l .and. q /=1) then
1246  write(eiram_unit,*) "ERROR in eiram_calc2: arrays have wrong length"
1247  write(eiram_unit,*) " length(X1), length(Y), ",size(lnx1),size(y)
1248  write(eiram_unit,*) " Length of X1 this array must be either = length(Y), or = 1"
1249  err = 100
1250  return
1251  end if
1252 
1253  if(.not.allocated(reactions)) then
1254  write(eiram_unit, *) "ERROR in eiram_calc2: EIRAM is not initialized"
1255  err = 50
1256  return
1257  end if
1258 
1259  m=size(reactions)
1260  if(id<0 .or. id>m) then
1261  write(eiram_unit,*) "ERROR in eiram_calc2: reaction index is out of range"
1262  write(eiram_unit,*) " index, min, max ",id,1,m
1263  err = 100
1264  return
1265  end if
1266 
1267  if(.not.allocated(reactions(id)%CREAC)) then
1268  write(eiram_unit,*) "ERROR in eiram_calc 2: array of fit coefficients is not allocated, Id ",id
1269  err = 50
1270  return
1271  end if
1272 
1273  m=reactions(id)%M
1274  if(m<1) then
1275  write(eiram_unit,*) "ERROR in eiram_calc2: wrong type of reaction"
1276  write(eiram_unit,*) " Double polynomial fit is requested for reaction ",&
1277  " described by a single polynomial fit"
1278  err = 120
1279  return
1280  end if
1281 
1282  n=reactions(id)%N
1283 
1284  do k=1,l
1285  y(k)=eiram_fit(reactions(id)%CREAC(:,n),m,lnx2(k))
1286  qk=min(q,k)
1287  do i=n-1,0,-1
1288  y(k)= y(k)*lnx1(qk)+eiram_fit(reactions(id)%CREAC(:,i),m,lnx2(k))
1289  end do
1290  y(k)=exp(y(k))
1291  end do
1292 
1293  end subroutine eiram_calc2
1294 
1295 !***********************************************************************************
1296 ! AUXILIARY
1297 !***********************************************************************************
1298 
1300  integer function countlinescontaining(filePath, string, err)
1301 
1302  character(*), intent(in) :: filePath
1304  character(*), intent(in) :: string
1306  integer, intent(out) :: err
1307 
1308  integer, parameter :: UNIT = 22
1309  character(MAX_LINE_LENGTH) :: line
1310 
1311  intrinsic trim
1312 
1313  open(unit, iostat=err, file=filepath, status='old', action='read')
1314  countlinescontaining = 0
1315  if(err /= 0) then
1316  write(eiram_unit, *) "ERROR in countLinesContaining while reading "//trim(filepath), err
1317  close(unit)
1318  err = 300
1319  return
1320  end if
1321  do
1322  read(unit, '(A)', iostat=err,end=100) line
1323  if(err /= 0) then
1324  write(eiram_unit, *) "ERROR in countLinesContaining while reading line of "//trim(filepath)
1325  close(unit)
1326  err = 300
1327  return
1328  end if
1329  if(index(line, string) > 0) countlinescontaining = countlinescontaining + 1
1330  end do
1331 100 err = 0
1332  close(unit)
1333  end function countlinescontaining
1334 
1335 end module eiram
1336 
logical, parameter debug
Definition: eiram.f90:131
integer, parameter reaction_length
Length of the string with the reaction equation.
Definition: eiram.f90:119
integer, parameter units_length
Length of the string with units.
Definition: eiram.f90:115
integer, parameter, public eiram_dp
Kind parameter for real numbers.
Definition: eiram.f90:110
integer, parameter max_line_length
Maximum length of a line read from an file.
Definition: eiram.f90:129
integer, public eiram_unit
Unit to which messages are written.
Definition: eiram.f90:112
integer, parameter description_length
Length of the string with description of the process.
Definition: eiram.f90:121
integer, parameter, public eiram_index_length
Length of the string with reaction index.
Definition: eiram.f90:127
subroutine, public eiram_load(filePath, fileName, err)
Initialization of the module from input files (data sets)
Definition: eiram.f90:206
subroutine, public eiram_deallocate(err)
Deallocate dynamic arrays used by this module.
Definition: eiram.f90:585
subroutine, public eiram_get_order(N, M, Id, err)
Return order of the polynomial for both variables.
Definition: eiram.f90:895
Data for one reaction.
Definition: eiram.f90:142
type(eiram_data), dimension(:), allocatable reactions
Array of reactions.
Definition: eiram.f90:190
external, public eiram_fit
Definition: eiram.f90:133
integer, parameter, public eiram_file_name_length
Length of the string with a name of an input file.
Definition: eiram.f90:123
integer function, public eiram_get_id(fileName, reactionType, reactionIndex, err)
Return the ID-index of the given reaction in the arrays of the module.
Definition: eiram.f90:857
integer, parameter, public eiram_type_length
Length of the string with reactiontype.
Definition: eiram.f90:125
integer, public eiram_nreact
current number of loaded reactions
Definition: eiram.f90:187
subroutine, public eiram_calc12(B, Id, lnX, err)
Reduce double polynomial fit to single polynomials for given values of the second variable...
Definition: eiram.f90:1142
real(dp) function eiram_fit(B, N, X)
Elementary function which calculates a single polynomial fit.
Definition: eiram.f90:77
integer, parameter quantity_length
Length of the string with description of a quantity.
Definition: eiram.f90:117
subroutine, public eiram_calc2(Y, Id, LnX1, LnX2, err)
Calculate double polynomial fit.
Definition: eiram.f90:1219
real(eiram_dp) function, public eiram_calc1_fast(Id, lnX)
Calculate single polynomial fit: scalar version W/O CHECKS.
Definition: eiram.f90:1042
Definition: eiram.f90:96
subroutine, public eiram_create_image(data, err)
Create a copy (image) of the reactions array stored in the module.
Definition: eiram.f90:931
subroutine, public eiram_calc1(Y, Id, lnX, err)
Calculate single polynomial fit.
Definition: eiram.f90:1060