added modules not previously uploaded in nocommon and gray-jintrac branches

Makefile

@@ -1,50 +1,114 @@
-# Executable name
-EXE=gray
+# Makefile for JETTO fortran 90 library 'gray'
+# Author: Lorenzo Figini (figini@ifp.cnr.it)
+#
+# Derived from
+# Makefile for JETTO fortran 77 library 'frantic'
+# Author: Derek Harting (d.harting@fz-juelich.de)
+# and
+# Makefile for ITCequ Fortran 90 library - definitions
+# P. Strand, elfps@chalmers.se
+# G. Corrigan, gcor@jet.uk
+# D. Harting, d.harting@fz-juelich.de
+# Set the environment from the top-level Makefile of JETTO
+# --------------------------------------------------------
+include ../include.mk
 
-# Objects list
-OBJ=gray.o grayl.o reflections.o green_func_p.o \
-    const_and_precisions.o itm_constants.o itm_types.o
+# library name
+# ------------
+LIBNAME=$(JLIBDIR)/libgray.a
 
-# Alternative search paths
-vpath %.f90 src
-vpath %.f src
+# List source and object files
+# ----------------------------
+FSRC=$(wildcard *.f *.f90)
+FOBJ0=$(FSRC:.f=.o)
+FOBJ=$(FOBJ0:.f90=.o)
 
-# Fortran compiler name and flags
-FC=gfortran
-FFLAGS=-O3
+########################################
+# Set up compiler specific environment #
+########################################
+# set default compiler
+# --------------------
+FC=$(F90)
 
-DIRECTIVES = -DREVISION="'$(shell svnversion src)'"
+# Set compiler debug and optimization flags
+# -----------------------------------------
+ifeq ("$(DBG)","")
+     FFLAGS= $(AUTO) -O3
+else
+     FFLAGS= $(AUTO) -O0 -g -Minform=inform -Mbounds -Mchkptr
+endif
 
-all: $(EXE)
+# Set include directories
+# -----------------------
+#INC=-I$(JCOMMON_STD)
+INC=
 
-# Build executable from object files
-$(EXE): $(OBJ)
-	$(FC) $(FFLAGS) -o $@ $^
+#####################################
+# Set up RULES specific environment #
+#####################################
 
-# Dependencies on modules
-gray.o: green_func_p.o reflections.o
-green_func_p.o: const_and_precisions.o
-const_and_precisions.o: itm_types.o itm_constants.o
-itm_constants.o: itm_types.o
+# suffixes
+# --------
+.SUFFIXES: .f90 .f .o .mod
 
-# General object compilation command
+# PHONY targets
+# -------------
+.PHONY: all clean realclean $(MASTER_RULES) gray
+
+# default target
+# --------------
+all:
+	@($(MAKE) -C ../ all) || exit $$?
+
+# Set rules passed to top-level Makefile
+# --------------------------------------
+$(MASTER_RULES):
+	@($(MAKE) -C ../ $@) || exit $$?
+
+# rule for libgray.a
+# --------------------
+$(LIBNAME): $(FOBJ)
+	ar vr $(LIBNAME) $?
+
+# synonym for libgray.a
+# -----------------------
+gray: $(LIBNAME)
+
+# create rule for compiling f90 files
+# -----------------------------------
 %.o: %.f90
-	$(FC) $(FFLAGS) -c $<
+	$(FC) -c $(FFLAGS) $(INC)  $<
 
-gray.o:gray.f green_func_p.o
-	$(FC) -cpp $(DIRECTIVES) $(FFLAGS) -c $< 
+# create rule for compiling f files
+# ---------------------------------
+%.o: %.f
+	$(FC) -c $(FFLAGS) $(INC)  $<
 
-grayl.o:grayl.f
-	$(FC) $(FFLAGS) -c $^ 
+# make 'realclean' option
+# -----------------------
+realclean:
+	rm -f *.o *.mod $(LIBNAME)
 
-.PHONY: clean install
-# Remove output files
+# make 'clean' option
+# -------------------
 clean:
-	rm -rf *.o *.mod $(EXE)
+	rm -f *.o
 
-install:
-	@if [ -f $(EXE) ]; then \
-	  cp $(EXE) ~/bin/; \
-	else \
-	  echo File $(EXE) does not exist. Run \'make\' first; \
-	fi
+# Dependencies
+# ------------
+gray_main.o: const_and_precisions.o graydata_flags.o
+gray-externals.o: green_func_p.o reflections.o beamdata.o const_and_precisions.o dispersion.o \
+                  graydata_par.o graydata_flags.o graydata_anequil.o
+green_func_p.o: const_and_precisions.o
+scatterspl.o: const_and_precisions.o
+beamdata.o: const_and_precisions.o
+dispersion.o: calcei.o dqagmv.o
+graydata_par.o: const_and_precisions.o
+graydata_flags.o: const_and_precisions.o
+graydata_anequil.o: const_and_precisions.o
+
+# Special rule to preprocess source file and include svn revision
+# ---------------------------------------------------------------
+DIRECTIVES = -DREVISION="'$(shell svnversion)'"
+gray-externals.o:gray-externals.f
+	$(FC) -c -Mpreprocess $(DIRECTIVES) $(FFLAGS) $(INC)  $<

Makefile.single

@@ -0,0 +1,172 @@
+# Executable name
+EXE=gray-single
+
+# Objects list
+OBJ= gray-single.o gray_main.o gray-externals.o grayl.o beamdata.o  \
+     const_and_precisions.o green_func_p.o reflections.o scatterspl.o \
+     dispersion.o calcei.o dqagmv.o graydata_par.o graydata_flags.o \
+     graydata_anequil.o
+
+# Alternative search paths
+vpath %.f90 src
+vpath %.f src
+
+########################################
+# Set up compiler specific environment #
+########################################
+# set default compiler
+# --------------------
+FC=gfortran
+
+# Set compiler debug and optimization flags
+# -----------------------------------------
+FFLAGS= -O0 -finit-real=nan -fcheck=all 
+
+#####################################
+# Set up RULES specific environment #
+#####################################
+
+DIRECTIVES = -DREVISION="'$(shell svnversion src)'"
+
+# suffixes
+# --------
+.SUFFIXES: .f90 .f .o .mod
+
+# PHONY targets
+# -------------
+.PHONY: clean install
+
+all: $(EXE)
+
+# Build executable from object files
+$(EXE): $(OBJ)
+	$(FC) $(FFLAGS) -o $@ $^
+
+# create rule for compiling f90 files
+# -----------------------------------
+%.o: %.f90
+	$(FC) -c $(FFLAGS) $<
+
+# create rule for compiling f files
+# ---------------------------------
+%.o: %.f
+	$(FC) -c $(FFLAGS) $<
+
+# create rule for files requiring preprocessing
+gray-externals.o: gray-externals.f
+	$(FC) -cpp $(DIRECTIVES) $(FFLAGS) -c $< 
+
+# make 'clean' option
+# -------------------
+clean:
+	rm -rf *.o *.mod $(EXE)
+
+# Dependencies
+# ------------
+gray-single.o: gray_main.o grayl.o
+gray_main.o: const_and_precisions.o graydata_flags.o
+gray-externals.o: green_func_p.o reflections.o beamdata.o const_and_precisions.o dispersion.o \
+                  graydata_par.o graydata_flags.o graydata_anequil.o
+green_func_p.o: const_and_precisions.o
+scatterspl.o: const_and_precisions.o
+beamdata.o: const_and_precisions.o
+dispersion.o: calcei.o dqagmv.o
+graydata_par.o: const_and_precisions.o
+graydata_flags.o: const_and_precisions.o
+graydata_anequil.o: const_and_precisions.o
+
+## library name
+## ------------
+#LIBNAME=$(JLIBDIR)/libgray.a
+#
+## List source and object files
+## ----------------------------
+#FSRC=$(wildcard *.f *.f90)
+#FOBJ0=$(FSRC:.f=.o)
+#FOBJ=$(FOBJ0:.f90=.o)
+#
+#########################################
+## Set up compiler specific environment #
+#########################################
+## set default compiler
+## --------------------
+#FC=$(F90)
+#
+## Set compiler debug and optimization flags
+## -----------------------------------------
+#ifeq ("$(DBG)","")
+#     FFLAGS= $(AUTO) -O3
+#else
+#     FFLAGS= $(AUTO) -O0 -g -Minform=inform -Mbounds -Mchkptr
+#endif
+#
+## Set include directories
+## -----------------------
+##INC=-I$(JCOMMON_STD)
+#INC=
+#
+######################################
+## Set up RULES specific environment #
+######################################
+#
+## suffixes
+## --------
+#.SUFFIXES: .f90 .f .o .mod
+#
+## PHONY targets
+## -------------
+#.PHONY: all clean realclean $(MASTER_RULES) gray
+#
+## default target
+## --------------
+#all:
+#	@($(MAKE) -C ../ all) || exit $$?
+#
+## Set rules passed to top-level Makefile
+## --------------------------------------
+#$(MASTER_RULES):
+#	@($(MAKE) -C ../ $@) || exit $$?
+#
+## rule for libgray.a
+## --------------------
+#$(LIBNAME): $(FOBJ)
+#	ar vr $(LIBNAME) $?
+#
+## synonym for libgray.a
+## -----------------------
+#gray: $(LIBNAME)
+#
+## create rule for compiling f90 files
+## -----------------------------------
+#%.o: %.f90
+#	$(FC) -c $(FFLAGS) $(INC)  $<
+#
+## create rule for compiling f files
+## ---------------------------------
+#%.o: %.f
+#	$(FC) -c $(FFLAGS) $(INC)  $<
+#
+## make 'realclean' option
+## -----------------------
+#realclean:
+#	rm -f *.o *.mod $(LIBNAME)
+#
+## make 'clean' option
+## -------------------
+#clean:
+#	rm -f *.o
+#
+## Dependencies
+## ------------
+#gray_main.o: const_and_precisions.o
+#gray-externals.o: green_func_p.o reflections.o beamdata.o const_and_precisions.o
+#green_func_p.o: const_and_precisions.o
+#scatterspl.o: const_and_precisions.o
+#beamdata.o: const_and_precisions.o
+#
+## Special rule to preprocess source file and include svn revision
+## ---------------------------------------------------------------
+#DIRECTIVES = -DREVISION="'$(shell svnversion)'"
+#gray-externals.o:gray-externals.f
+#	$(FC) -c -Mpreprocess $(DIRECTIVES) $(FFLAGS) $(INC)  $<
+

Makefile.standalone

@@ -0,0 +1,68 @@
+# Executable name
+EXE=gray-jintrac
+LIBFILE=lib$(EXE).a
+
+# Objects list
+OBJ=gray_main.o gray-externals.o grayl.o reflections.o scatterspl.o \
+    beamdata.o green_func_p.o const_and_precisions.o dispersion.o \
+    calcei.o dqagmv.o graydata_par.o graydata_flags.o graydata_anequil.o
+
+# Alternative search paths
+vpath %.f90 src
+vpath %.f src
+
+# Fortran compiler name and flags
+FC=gfortran
+FFLAGS=-O0 -g -fbounds-check -Wall
+#FFLAGS=-O3
+
+DIRECTIVES = -DREVISION="'$(shell svnversion)'"
+
+all: library
+
+# Build executable from object files
+$(EXE): main.o $(OBJ)
+	$(FC) $(FFLAGS) -o $@ $^
+
+$(LIBFILE): $(OBJ)
+	rm -f $@
+	ar -rv $@ $^
+
+# Dependencies on modules
+main.o: const_and_precisions.o
+gray_main.o: const_and_precisions.o graydata_flags.o
+gray-externals.o: green_func_p.o reflections.o beamdata.o dispersion.o \
+                  graydata_par.o graydata_flags.o graydata_anequil.o \
+                  const_and_precisions.o
+green_func_p.o: const_and_precisions.o
+scatterspl.o: const_and_precisions.o
+beamdata.o: const_and_precisions.o
+dispersion.o: calcei.o dqagmv.o
+graydata_par.o: const_and_precisions.o
+graydata_flags.o: const_and_precisions.o
+graydata_anequil.o: const_and_precisions.o
+
+# General object compilation command
+%.o: %.f90
+	$(FC) $(FFLAGS) -c $<
+
+gray-externals.o:gray-externals.f
+	$(FC) -cpp $(DIRECTIVES) $(FFLAGS) -c $< 
+
+grayl.o:grayl.f
+	$(FC) $(FFLAGS) -c $<
+
+.PHONY: library clean install
+
+library: $(LIBFILE)
+
+# Remove output files
+clean:
+	rm -rf *.o *.mod $(EXE)
+
+install:
+	@if [ -f $(EXE) ]; then \
+	  cp $(EXE) ~/bin/; \
+	else \
+	  echo File $(EXE) does not exist. Run \'make\' first; \
+	fi

src/beamdata.f90

@@ -0,0 +1,72 @@
+module beamdata
+  use const_and_precisions, only : r8
+  implicit none
+
+  integer, parameter :: jmx=31,kmx=36,nmx=8000
+  integer, save :: nrayr,nrayth,nstep
+  real(r8), dimension(:,:,:), allocatable, save :: psjki,ppabs,ccci,tauv,alphav
+  real(r8), dimension(:,:,:), allocatable, save :: pdjki,currj,didst
+  integer, dimension(:,:), allocatable, save :: iiv,iop,iow
+  real(r8), dimension(:,:), allocatable, save :: tau1v
+  real(r8), dimension(:), allocatable, save :: q
+  real(r8), dimension(:,:,:), allocatable, save :: yyrfl !(:,:,6)
+  real(r8), dimension(:,:,:), allocatable, save :: ywrk,ypwrk !(6,:,:)
+  real(r8), dimension(:,:,:), allocatable, save :: xc,xco,du1,du1o !(3,:,:)
+  real(r8), dimension(:,:,:), allocatable, save :: gri,dgrad2v !(3,:,:)
+  real(r8), dimension(:,:,:,:), allocatable, save :: ggri !(3,3,:,:)
+  real(r8), dimension(:,:), allocatable, save :: grad2
+  real(r8), dimension(:), allocatable, save :: dffiu,ddffiu
+
+contains
+
+  subroutine alloc_beam(ierr)
+    implicit none
+    integer, intent(out) :: ierr
+    
+    call dealloc_beam
+    allocate(psjki(nrayr,nrayth,nstep),  ppabs(nrayr,nrayth,nstep), &
+             pdjki(nrayr,nrayth,nstep),   ccci(nrayr,nrayth,nstep), &
+             currj(nrayr,nrayth,nstep),  didst(nrayr,nrayth,nstep), &
+              tauv(nrayr,nrayth,nstep), alphav(nrayr,nrayth,nstep), &
+               iiv(nrayr,nrayth),          iop(nrayr,nrayth),       &
+               iow(nrayr,nrayth),        tau1v(nrayr,nrayth),       &
+                 q(nrayr),               yyrfl(nrayr,nrayth,6),     &
+            ywrk(6,nrayr,nrayth),      ypwrk(6,nrayr,nrayth),       &
+              xc(3,nrayr,nrayth),        xco(3,nrayr,nrayth),       &
+             du1(3,nrayr,nrayth),       du1o(3,nrayr,nrayth),       &
+             gri(3,nrayr,nrayth),    dgrad2v(3,nrayr,nrayth),       &
+          ggri(3,3,nrayr,nrayth),        grad2(nrayr,nrayth),       &
+             dffiu(nrayr),              ddffiu(nrayr),     stat=ierr)
+    if (ierr/=0) call dealloc_beam
+  end subroutine alloc_beam
+
+  subroutine dealloc_beam
+    implicit none
+    if (allocated(psjki))   deallocate(psjki)
+    if (allocated(ppabs))   deallocate(ppabs)
+    if (allocated(pdjki))   deallocate(pdjki)
+    if (allocated(ccci))    deallocate(ccci)
+    if (allocated(currj))   deallocate(currj)
+    if (allocated(didst))   deallocate(didst)
+    if (allocated(tauv))    deallocate(tauv)
+    if (allocated(alphav))  deallocate(alphav)
+    if (allocated(iiv))     deallocate(iiv)
+    if (allocated(iop))     deallocate(iop)
+    if (allocated(iow))     deallocate(iow)
+    if (allocated(tau1v))   deallocate(tau1v)
+    if (allocated(q))       deallocate(q)
+    if (allocated(yyrfl))   deallocate(yyrfl)
+    if (allocated(ywrk))    deallocate(ywrk)
+    if (allocated(ypwrk))   deallocate(ypwrk)
+    if (allocated(xc))      deallocate(xc)
+    if (allocated(xco))     deallocate(xco)
+    if (allocated(du1))     deallocate(du1)
+    if (allocated(du1o))    deallocate(du1o)
+    if (allocated(gri))     deallocate(gri)
+    if (allocated(dgrad2v)) deallocate(dgrad2v)
+    if (allocated(ggri))    deallocate(ggri)
+    if (allocated(grad2))   deallocate(grad2)
+    if (allocated(dffiu))   deallocate(dffiu)
+    if (allocated(ddffiu))  deallocate(ddffiu)
+  end subroutine dealloc_beam
+end module beamdata

src/calcei.f90

@@ -0,0 +1,608 @@
+module calcei_mod
+
+ contains
+
+! ======================================================================
+! nist guide to available math software.
+! fullsource for module ei from package specfun.
+! retrieved from netlib on fri mar 26 05:52:39 1999.
+! ======================================================================
+      subroutine calcei(arg,result,int)
+!----------------------------------------------------------------------
+!
+! this fortran 77 packet computes the exponential integrals eint(x),
+!  e1(x), and  exp(-x)*eint(x)  for real arguments  x  where
+!
+!           integral (from t=-infinity to t=x) (exp(t)/t),  x > 0,
+!  eint(x) =
+!          -integral (from t=-x to t=infinity) (exp(t)/t),  x < 0,
+!
+!  and where the first integral is a principal value integral.
+!  the packet contains three function type subprograms: ei, eone,
+!  and expei;  and one subroutine type subprogram: calcei.  the
+!  calling statements for the primary entries are
+!
+!                 y = eint(x),          where  x .ne. 0,
+!
+!                 y = eone(x),          where  x .gt. 0,
+!  and
+!                 y = expei(x),         where  x .ne. 0,
+!
+!  and where the entry points correspond to the functions eint(x),
+!  e1(x), and exp(-x)*eint(x), respectively.  the routine calcei
+!  is intended for internal packet use only, all computations within
+!  the packet being concentrated in this routine.  the function
+!  subprograms invoke calcei with the fortran statement
+!         call calcei(arg,result,int)
+!  where the parameter usage is as follows
+!
+!     function                  parameters for calcei
+!       call                 arg             result         int
+!
+!      eint(x)            x .ne. 0         eint(x)           1
+!      eone(x)            x .gt. 0        -eint(-x)          2
+!      expei(x)           x .ne. 0         exp(-x)*eint(x)   3
+!----------------------------------------------------------------------
+      integer i,int
+      double precision&
+     &       a,arg,b,c,d,exp40,e,ei,f,four,fourty,frac,half,one,p,&
+     &       plg,px,p037,p1,p2,q,qlg,qx,q1,q2,r,result,s,six,sump,&
+     &       sumq,t,three,twelve,two,two4,w,x,xbig,xinf,xmax,xmx0,&
+     &       x0,x01,x02,x11,y,ysq,zero
+      dimension  a(7),b(6),c(9),d(9),e(10),f(10),p(10),q(10),r(10),&
+     &   s(9),p1(10),q1(9),p2(10),q2(9),plg(4),qlg(4),px(10),qx(10)
+!----------------------------------------------------------------------
+!  mathematical constants
+!   exp40 = exp(40)
+!   x0 = zero of ei
+!   x01/x11 + x02 = zero of ei to extra precision
+!----------------------------------------------------------------------
+      data zero,p037,half,one,two/0.0d0,0.037d0,0.5d0,1.0d0,2.0d0/,&
+     &     three,four,six,twelve,two4/3.0d0,4.0d0,6.0d0,12.d0,24.0d0/,&
+     &     fourty,exp40/40.0d0,2.3538526683701998541d17/,&
+     &     x01,x11,x02/381.5d0,1024.0d0,-5.1182968633365538008d-5/,&
+     &     x0/3.7250741078136663466d-1/
+!----------------------------------------------------------------------
+! machine-dependent constants
+!----------------------------------------------------------------------
+      data xinf/1.79d+308/,xmax/716.351d0/,xbig/701.84d0/
+!----------------------------------------------------------------------
+! coefficients  for -1.0 <= x < 0.0
+!----------------------------------------------------------------------
+      data a/1.1669552669734461083368d2, 2.1500672908092918123209d3,&
+     &       1.5924175980637303639884d4, 8.9904972007457256553251d4,&
+     &       1.5026059476436982420737d5,-1.4815102102575750838086d5,&
+     &       5.0196785185439843791020d0/
+      data b/4.0205465640027706061433d1, 7.5043163907103936624165d2,&
+     &       8.1258035174768735759855d3, 5.2440529172056355429883d4,&
+     &       1.8434070063353677359298d5, 2.5666493484897117319268d5/
+!----------------------------------------------------------------------
+! coefficients for -4.0 <= x < -1.0
+!----------------------------------------------------------------------
+      data c/3.828573121022477169108d-1, 1.107326627786831743809d+1,&
+     &       7.246689782858597021199d+1, 1.700632978311516129328d+2,&
+     &       1.698106763764238382705d+2, 7.633628843705946890896d+1,&
+     &       1.487967702840464066613d+1, 9.999989642347613068437d-1,&
+     &       1.737331760720576030932d-8/
+      data d/8.258160008564488034698d-2, 4.344836335509282083360d+0,&
+     &       4.662179610356861756812d+1, 1.775728186717289799677d+2,&
+     &       2.953136335677908517423d+2, 2.342573504717625153053d+2,&
+     &       9.021658450529372642314d+1, 1.587964570758947927903d+1,&
+     &       1.000000000000000000000d+0/
+!----------------------------------------------------------------------
+! coefficients for x < -4.0
+!----------------------------------------------------------------------
+      data e/1.3276881505637444622987d+2,3.5846198743996904308695d+4,&
+     &       1.7283375773777593926828d+5,2.6181454937205639647381d+5,&
+     &       1.7503273087497081314708d+5,5.9346841538837119172356d+4,&
+     &       1.0816852399095915622498d+4,1.0611777263550331766871d03,&
+     &       5.2199632588522572481039d+1,9.9999999999999999087819d-1/
+      data f/3.9147856245556345627078d+4,2.5989762083608489777411d+5,&
+     &       5.5903756210022864003380d+5,5.4616842050691155735758d+5,&
+     &       2.7858134710520842139357d+5,7.9231787945279043698718d+4,&
+     &       1.2842808586627297365998d+4,1.1635769915320848035459d+3,&
+     &       5.4199632588522559414924d+1,1.0d0/
+!----------------------------------------------------------------------
+!  coefficients for rational approximation to ln(x/a), |1-x/a| < .1
+!----------------------------------------------------------------------
+      data plg/-2.4562334077563243311d+01,2.3642701335621505212d+02,&
+     &         -5.4989956895857911039d+02,3.5687548468071500413d+02/
+      data qlg/-3.5553900764052419184d+01,1.9400230218539473193d+02,&
+     &         -3.3442903192607538956d+02,1.7843774234035750207d+02/
+!----------------------------------------------------------------------
+! coefficients for  0.0 < x < 6.0,
+!  ratio of chebyshev polynomials
+!----------------------------------------------------------------------
+      data p/-1.2963702602474830028590d01,-1.2831220659262000678155d03,&
+     &       -1.4287072500197005777376d04,-1.4299841572091610380064d06,&
+     &       -3.1398660864247265862050d05,-3.5377809694431133484800d08,&
+     &        3.1984354235237738511048d08,-2.5301823984599019348858d10,&
+     &        1.2177698136199594677580d10,-2.0829040666802497120940d11/
+      data q/ 7.6886718750000000000000d01,-5.5648470543369082846819d03,&
+     &        1.9418469440759880361415d05,-4.2648434812177161405483d06,&
+     &        6.4698830956576428587653d07,-7.0108568774215954065376d08,&
+     &        5.4229617984472955011862d09,-2.8986272696554495342658d10,&
+     &        9.8900934262481749439886d10,-8.9673749185755048616855d10/
+!----------------------------------------------------------------------
+! j-fraction coefficients for 6.0 <= x < 12.0
+!----------------------------------------------------------------------
+      data r/-2.645677793077147237806d00,-2.378372882815725244124d00,&
+     &       -2.421106956980653511550d01, 1.052976392459015155422d01,&
+     &        1.945603779539281810439d01,-3.015761863840593359165d01,&
+     &        1.120011024227297451523d01,-3.988850730390541057912d00,&
+     &        9.565134591978630774217d00, 9.981193787537396413219d-1/
+      data s/ 1.598517957704779356479d-4, 4.644185932583286942650d00,&
+     &        3.697412299772985940785d02,-8.791401054875438925029d00,&
+     &        7.608194509086645763123d02, 2.852397548119248700147d01,&
+     &        4.731097187816050252967d02,-2.369210235636181001661d02,&
+     &        1.249884822712447891440d00/
+!----------------------------------------------------------------------
+! j-fraction coefficients for 12.0 <= x < 24.0
+!----------------------------------------------------------------------
+      data p1/-1.647721172463463140042d00,-1.860092121726437582253d01,&
+     &        -1.000641913989284829961d01,-2.105740799548040450394d01,&
+     &        -9.134835699998742552432d-1,-3.323612579343962284333d01,&
+     &         2.495487730402059440626d01, 2.652575818452799819855d01,&
+     &        -1.845086232391278674524d00, 9.999933106160568739091d-1/
+      data q1/ 9.792403599217290296840d01, 6.403800405352415551324d01,&
+     &         5.994932325667407355255d01, 2.538819315630708031713d02,&
+     &         4.429413178337928401161d01, 1.192832423968601006985d03,&
+     &         1.991004470817742470726d02,-1.093556195391091143924d01,&
+     &         1.001533852045342697818d00/
+!----------------------------------------------------------------------
+! j-fraction coefficients for  x .ge. 24.0
+!----------------------------------------------------------------------
+      data p2/ 1.75338801265465972390d02,-2.23127670777632409550d02,&
+     &        -1.81949664929868906455d01,-2.79798528624305389340d01,&
+     &        -7.63147701620253630855d00,-1.52856623636929636839d01,&
+     &        -7.06810977895029358836d00,-5.00006640413131002475d00,&
+     &        -3.00000000320981265753d00, 1.00000000000000485503d00/
+      data q2/ 3.97845977167414720840d04, 3.97277109100414518365d00,&
+     &         1.37790390235747998793d02, 1.17179220502086455287d02,&
+     &         7.04831847180424675988d01,-1.20187763547154743238d01,&
+     &        -7.99243595776339741065d00,-2.99999894040324959612d00,&
+     &         1.99999999999048104167d00/
+!----------------------------------------------------------------------
+      x = arg
+      if (x .eq. zero) then
+            ei = -xinf
+            if (int .eq. 2) ei = -ei
+         else if ((x .lt. zero) .or. (int .eq. 2)) then
+!----------------------------------------------------------------------
+! calculate ei for negative argument or for e1.
+!----------------------------------------------------------------------
+            y = abs(x)
+            if (y .le. one) then
+                  sump = a(7) * y + a(1)
+                  sumq = y + b(1)
+                  do 110 i = 2, 6
+                     sump = sump * y + a(i)
+                     sumq = sumq * y + b(i)
+  110             continue
+                  ei = log(y) - sump / sumq
+                  if (int .eq. 3) ei = ei * exp(y)
+               else if (y .le. four) then
+                  w = one / y
+                  sump = c(1)
+                  sumq = d(1)
+                  do 130 i = 2, 9
+                     sump = sump * w + c(i)
+                     sumq = sumq * w + d(i)
+  130             continue
+                  ei = - sump / sumq
+                  if (int .ne. 3) ei = ei * exp(-y)
+               else
+                  if ((y .gt. xbig) .and. (int .lt. 3)) then
+                        ei = zero
+                     else
+                        w = one / y
+                        sump = e(1)
+                        sumq = f(1)
+                        do 150 i = 2, 10
+                           sump = sump * w + e(i)
+                           sumq = sumq * w + f(i)
+  150                   continue
+                        ei = -w * (one - w * sump / sumq )
+                        if (int .ne. 3) ei = ei * exp(-y)
+                  end if
+            end if
+            if (int .eq. 2) ei = -ei
+         else if (x .lt. six) then
+!----------------------------------------------------------------------
+!  to improve conditioning, rational approximations are expressed
+!    in terms of chebyshev polynomials for 0 <= x < 6, and in
+!    continued fraction form for larger x.
+!----------------------------------------------------------------------
+            t = x + x
+            t = t / three - two
+            px(1) = zero
+            qx(1) = zero
+            px(2) = p(1)
+            qx(2) = q(1)
+            do 210 i = 2, 9
+               px(i+1) = t * px(i) - px(i-1) + p(i)
+               qx(i+1) = t * qx(i) - qx(i-1) + q(i)
+  210       continue
+            sump = half * t * px(10) - px(9) + p(10)
+            sumq = half * t * qx(10) - qx(9) + q(10)
+            frac = sump / sumq
+            xmx0 = (x - x01/x11) - x02
+            if (abs(xmx0) .ge. p037) then
+                  ei = log(x/x0) + xmx0 * frac
+                  if (int .eq. 3) ei = exp(-x) * ei
+               else
+!----------------------------------------------------------------------
+! special approximation to  ln(x/x0)  for x close to x0
+!----------------------------------------------------------------------
+                  y = xmx0 / (x + x0)
+                  ysq = y*y
+                  sump = plg(1)
+                  sumq = ysq + qlg(1)
+                  do 220 i = 2, 4
+                     sump = sump*ysq + plg(i)
+                     sumq = sumq*ysq + qlg(i)
+  220             continue
+                  ei = (sump / (sumq*(x+x0)) + frac) * xmx0
+                  if (int .eq. 3) ei = exp(-x) * ei
+            end if
+         else if (x .lt. twelve) then
+            frac = zero
+            do 230 i = 1, 9
+               frac = s(i) / (r(i) + x + frac)
+  230       continue
+            ei = (r(10) + frac) / x
+            if (int .ne. 3) ei = ei * exp(x)
+         else if (x .le. two4) then
+            frac = zero
+            do 240 i = 1, 9
+               frac = q1(i) / (p1(i) + x + frac)
+  240       continue
+            ei = (p1(10) + frac) / x
+            if (int .ne. 3) ei = ei * exp(x)
+         else
+            if ((x .ge. xmax) .and. (int .lt. 3)) then
+                  ei = xinf
+               else
+                  y = one / x
+                  frac = zero
+                  do 250 i = 1, 9
+                     frac = q2(i) / (p2(i) + x + frac)
+  250             continue
+                  frac = p2(10) + frac
+                  ei = y + y * y * frac
+                  if (int .ne. 3) then
+                        if (x .le. xmax-two4) then
+                              ei = ei * exp(x)
+                           else
+!----------------------------------------------------------------------
+! calculation reformulated to avoid premature overflow
+!----------------------------------------------------------------------
+                              ei = (ei * exp(x-fourty)) * exp40
+                        end if
+                  end if
+            end if
+      end if
+      result = ei
+      return
+!---------- last line of calcei ----------
+      end
+      function eint(x)
+!--------------------------------------------------------------------
+!
+! this function program computes approximate values for the
+!   exponential integral  eint(x), where  x  is real.
+!
+!  author: w. j. cody
+!
+!  latest modification: january 12, 1988
+!
+!--------------------------------------------------------------------
+      integer int
+      double precision  eint, x, result
+!--------------------------------------------------------------------
+      int = 1
+      call calcei(x,result,int)
+      eint = result
+      return
+!---------- last line of ei ----------
+      end
+      function expei(x)
+!--------------------------------------------------------------------
+!
+! this function program computes approximate values for the
+!   function  exp(-x) * eint(x), where  eint(x)  is the exponential
+!   integral, and  x  is real.
+!
+!  author: w. j. cody
+!
+!  latest modification: january 12, 1988
+!
+!--------------------------------------------------------------------
+      integer int
+      double precision  expei, x, result
+!--------------------------------------------------------------------
+      int = 3
+      call calcei(x,result,int)
+      expei = result
+      return
+!---------- last line of expei ----------
+      end
+      function eone(x)
+!--------------------------------------------------------------------
+!
+! this function program computes approximate values for the
+!   exponential integral e1(x), where  x  is real.
+!
+!  author: w. j. cody
+!
+!  latest modification: january 12, 1988
+!
+!--------------------------------------------------------------------
+      integer int
+      double precision  eone, x, result
+!--------------------------------------------------------------------
+      int = 2
+      call calcei(x,result,int)
+      eone = result
+      return
+!---------- last line of eone ----------
+      end
+!
+!     calcei3 = calcei for int=3
+!
+! ======================================================================
+      subroutine calcei3(arg,result)
+!----------------------------------------------------------------------
+!
+! this fortran 77 packet computes the exponential integrals eint(x),
+!  e1(x), and  exp(-x)*eint(x)  for real arguments  x  where
+!
+!           integral (from t=-infinity to t=x) (exp(t)/t),  x > 0,
+!  eint(x) =
+!          -integral (from t=-x to t=infinity) (exp(t)/t),  x < 0,
+!
+!  and where the first integral is a principal value integral.
+!  the packet contains three function type subprograms: ei, eone,
+!  and expei;  and one subroutine type subprogram: calcei.  the
+!  calling statements for the primary entries are
+!
+!                 y = eint(x),          where  x .ne. 0,
+!
+!                 y = eone(x),          where  x .gt. 0,
+!  and
+!                 y = expei(x),         where  x .ne. 0,
+!
+!  and where the entry points correspond to the functions eint(x),
+!  e1(x), and exp(-x)*eint(x), respectively.  the routine calcei
+!  is intended for internal packet use only, all computations within
+!  the packet being concentrated in this routine.  the function
+!  subprograms invoke calcei with the fortran statement
+!         call calcei(arg,result,int)
+!  where the parameter usage is as follows
+!
+!     function                  parameters for calcei
+!       call                 arg             result         int
+!
+!      eint(x)            x .ne. 0         eint(x)           1
+!      eone(x)            x .gt. 0        -eint(-x)          2
+!      expei(x)           x .ne. 0         exp(-x)*eint(x)   3
+!----------------------------------------------------------------------
+      integer i,int
+      double precision&
+     &       a,arg,b,c,d,exp40,e,ei,f,four,fourty,frac,half,one,p,&
+     &       plg,px,p037,p1,p2,q,qlg,qx,q1,q2,r,result,s,six,sump,&
+     &       sumq,t,three,twelve,two,two4,w,x,xbig,xinf,xmax,xmx0,&
+     &       x0,x01,x02,x11,y,ysq,zero
+      dimension  a(7),b(6),c(9),d(9),e(10),f(10),p(10),q(10),r(10),&
+     &   s(9),p1(10),q1(9),p2(10),q2(9),plg(4),qlg(4),px(10),qx(10)
+!----------------------------------------------------------------------
+!  mathematical constants
+!   exp40 = exp(40)
+!   x0 = zero of ei
+!   x01/x11 + x02 = zero of ei to extra precision
+!----------------------------------------------------------------------
+      data zero,p037,half,one,two/0.0d0,0.037d0,0.5d0,1.0d0,2.0d0/,&
+     &     three,four,six,twelve,two4/3.0d0,4.0d0,6.0d0,12.d0,24.0d0/,&
+     &     fourty,exp40/40.0d0,2.3538526683701998541d17/,&
+     &     x01,x11,x02/381.5d0,1024.0d0,-5.1182968633365538008d-5/,&
+     &     x0/3.7250741078136663466d-1/
+!----------------------------------------------------------------------
+! machine-dependent constants
+!----------------------------------------------------------------------
+      data xinf/1.79d+308/,xmax/716.351d0/,xbig/701.84d0/
+!----------------------------------------------------------------------
+! coefficients  for -1.0 <= x < 0.0
+!----------------------------------------------------------------------
+      data a/1.1669552669734461083368d2, 2.1500672908092918123209d3,&
+     &       1.5924175980637303639884d4, 8.9904972007457256553251d4,&
+     &       1.5026059476436982420737d5,-1.4815102102575750838086d5,&
+     &       5.0196785185439843791020d0/
+      data b/4.0205465640027706061433d1, 7.5043163907103936624165d2,&
+     &       8.1258035174768735759855d3, 5.2440529172056355429883d4,&
+     &       1.8434070063353677359298d5, 2.5666493484897117319268d5/
+!----------------------------------------------------------------------
+! coefficients for -4.0 <= x < -1.0
+!----------------------------------------------------------------------
+      data c/3.828573121022477169108d-1, 1.107326627786831743809d+1,&
+     &       7.246689782858597021199d+1, 1.700632978311516129328d+2,&
+     &       1.698106763764238382705d+2, 7.633628843705946890896d+1,&
+     &       1.487967702840464066613d+1, 9.999989642347613068437d-1,&
+     &       1.737331760720576030932d-8/
+      data d/8.258160008564488034698d-2, 4.344836335509282083360d+0,&
+     &       4.662179610356861756812d+1, 1.775728186717289799677d+2,&
+     &       2.953136335677908517423d+2, 2.342573504717625153053d+2,&
+     &       9.021658450529372642314d+1, 1.587964570758947927903d+1,&
+     &       1.000000000000000000000d+0/
+!----------------------------------------------------------------------
+! coefficients for x < -4.0
+!----------------------------------------------------------------------
+      data e/1.3276881505637444622987d+2,3.5846198743996904308695d+4,&
+     &       1.7283375773777593926828d+5,2.6181454937205639647381d+5,&
+     &       1.7503273087497081314708d+5,5.9346841538837119172356d+4,&
+     &       1.0816852399095915622498d+4,1.0611777263550331766871d03,&
+     &       5.2199632588522572481039d+1,9.9999999999999999087819d-1/
+      data f/3.9147856245556345627078d+4,2.5989762083608489777411d+5,&
+     &       5.5903756210022864003380d+5,5.4616842050691155735758d+5,&
+     &       2.7858134710520842139357d+5,7.9231787945279043698718d+4,&
+     &       1.2842808586627297365998d+4,1.1635769915320848035459d+3,&
+     &       5.4199632588522559414924d+1,1.0d0/
+!----------------------------------------------------------------------
+!  coefficients for rational approximation to ln(x/a), |1-x/a| < .1
+!----------------------------------------------------------------------
+      data plg/-2.4562334077563243311d+01,2.3642701335621505212d+02,&
+     &         -5.4989956895857911039d+02,3.5687548468071500413d+02/
+      data qlg/-3.5553900764052419184d+01,1.9400230218539473193d+02,&
+     &         -3.3442903192607538956d+02,1.7843774234035750207d+02/
+!----------------------------------------------------------------------
+! coefficients for  0.0 < x < 6.0,
+!  ratio of chebyshev polynomials
+!----------------------------------------------------------------------
+      data p/-1.2963702602474830028590d01,-1.2831220659262000678155d03,&
+     &       -1.4287072500197005777376d04,-1.4299841572091610380064d06,&
+     &       -3.1398660864247265862050d05,-3.5377809694431133484800d08,&
+     &        3.1984354235237738511048d08,-2.5301823984599019348858d10,&
+     &        1.2177698136199594677580d10,-2.0829040666802497120940d11/
+      data q/ 7.6886718750000000000000d01,-5.5648470543369082846819d03,&
+     &        1.9418469440759880361415d05,-4.2648434812177161405483d06,&
+     &        6.4698830956576428587653d07,-7.0108568774215954065376d08,&
+     &        5.4229617984472955011862d09,-2.8986272696554495342658d10,&
+     &        9.8900934262481749439886d10,-8.9673749185755048616855d10/
+!----------------------------------------------------------------------
+! j-fraction coefficients for 6.0 <= x < 12.0
+!----------------------------------------------------------------------
+      data r/-2.645677793077147237806d00,-2.378372882815725244124d00,&
+     &       -2.421106956980653511550d01, 1.052976392459015155422d01,&
+     &        1.945603779539281810439d01,-3.015761863840593359165d01,&
+     &        1.120011024227297451523d01,-3.988850730390541057912d00,&
+     &        9.565134591978630774217d00, 9.981193787537396413219d-1/
+      data s/ 1.598517957704779356479d-4, 4.644185932583286942650d00,&
+     &        3.697412299772985940785d02,-8.791401054875438925029d00,&
+     &        7.608194509086645763123d02, 2.852397548119248700147d01,&
+     &        4.731097187816050252967d02,-2.369210235636181001661d02,&
+     &        1.249884822712447891440d00/
+!----------------------------------------------------------------------
+! j-fraction coefficients for 12.0 <= x < 24.0
+!----------------------------------------------------------------------
+      data p1/-1.647721172463463140042d00,-1.860092121726437582253d01,&
+     &        -1.000641913989284829961d01,-2.105740799548040450394d01,&
+     &        -9.134835699998742552432d-1,-3.323612579343962284333d01,&
+     &         2.495487730402059440626d01, 2.652575818452799819855d01,&
+     &        -1.845086232391278674524d00, 9.999933106160568739091d-1/
+      data q1/ 9.792403599217290296840d01, 6.403800405352415551324d01,&
+     &         5.994932325667407355255d01, 2.538819315630708031713d02,&
+     &         4.429413178337928401161d01, 1.192832423968601006985d03,&
+     &         1.991004470817742470726d02,-1.093556195391091143924d01,&
+     &         1.001533852045342697818d00/
+!----------------------------------------------------------------------
+! j-fraction coefficients for  x .ge. 24.0
+!----------------------------------------------------------------------
+      data p2/ 1.75338801265465972390d02,-2.23127670777632409550d02,&
+     &        -1.81949664929868906455d01,-2.79798528624305389340d01,&
+     &        -7.63147701620253630855d00,-1.52856623636929636839d01,&
+     &        -7.06810977895029358836d00,-5.00006640413131002475d00,&
+     &        -3.00000000320981265753d00, 1.00000000000000485503d00/
+      data q2/ 3.97845977167414720840d04, 3.97277109100414518365d00,&
+     &         1.37790390235747998793d02, 1.17179220502086455287d02,&
+     &         7.04831847180424675988d01,-1.20187763547154743238d01,&
+     &        -7.99243595776339741065d00,-2.99999894040324959612d00,&
+     &         1.99999999999048104167d00/
+!----------------------------------------------------------------------
+      data int/ 3/
+      x = arg
+      if (x .eq. zero) then
+            ei = -xinf
+         else if ((x .lt. zero)) then
+!----------------------------------------------------------------------
+! calculate ei for negative argument or for e1.
+!----------------------------------------------------------------------
+            y = abs(x)
+            if (y .le. one) then
+                  sump = a(7) * y + a(1)
+                  sumq = y + b(1)
+                  do 110 i = 2, 6
+                     sump = sump * y + a(i)
+                     sumq = sumq * y + b(i)
+  110             continue
+                  ei = (log(y) - sump / sumq ) * exp(y)
+               else if (y .le. four) then
+                  w = one / y
+                  sump = c(1)
+                  sumq = d(1)
+                  do 130 i = 2, 9
+                     sump = sump * w + c(i)
+                     sumq = sumq * w + d(i)
+  130             continue
+                  ei = - sump / sumq
+               else
+                        w = one / y
+                        sump = e(1)
+                        sumq = f(1)
+                        do 150 i = 2, 10
+                           sump = sump * w + e(i)
+                           sumq = sumq * w + f(i)
+  150                   continue
+                        ei = -w * (one - w * sump / sumq )
+            end if
+         else if (x .lt. six) then
+!----------------------------------------------------------------------
+!  to improve conditioning, rational approximations are expressed
+!    in terms of chebyshev polynomials for 0 <= x < 6, and in
+!    continued fraction form for larger x.
+!----------------------------------------------------------------------
+            t = x + x
+            t = t / three - two
+            px(1) = zero
+            qx(1) = zero
+            px(2) = p(1)
+            qx(2) = q(1)
+            do 210 i = 2, 9
+               px(i+1) = t * px(i) - px(i-1) + p(i)
+               qx(i+1) = t * qx(i) - qx(i-1) + q(i)
+  210       continue
+            sump = half * t * px(10) - px(9) + p(10)
+            sumq = half * t * qx(10) - qx(9) + q(10)
+            frac = sump / sumq
+            xmx0 = (x - x01/x11) - x02
+            if (abs(xmx0) .ge. p037) then
+                  ei = exp(-x) * ( log(x/x0) + xmx0 * frac )
+               else
+!----------------------------------------------------------------------
+! special approximation to  ln(x/x0)  for x close to x0
+!----------------------------------------------------------------------
+                  y = xmx0 / (x + x0)
+                  ysq = y*y
+                  sump = plg(1)
+                  sumq = ysq + qlg(1)
+                  do 220 i = 2, 4
+                     sump = sump*ysq + plg(i)
+                     sumq = sumq*ysq + qlg(i)
+  220             continue
+                  ei = exp(-x) * (sump / (sumq*(x+x0)) + frac) * xmx0
+            end if
+         else if (x .lt. twelve) then
+            frac = zero
+            do 230 i = 1, 9
+               frac = s(i) / (r(i) + x + frac)
+  230       continue
+            ei = (r(10) + frac) / x
+         else if (x .le. two4) then
+            frac = zero
+            do 240 i = 1, 9
+               frac = q1(i) / (p1(i) + x + frac)
+  240       continue
+            ei = (p1(10) + frac) / x
+         else
+                  y = one / x
+                  frac = zero
+                  do 250 i = 1, 9
+                     frac = q2(i) / (p2(i) + x + frac)
+  250             continue
+                  frac = p2(10) + frac
+                  ei = y + y * y * frac
+      end if
+      result = ei
+      return
+!---------- last line of calcei ----------
+      end
+
+end module calcei_mod

src/const_and_precisions.f90

@@ -1,18 +1,22 @@
 !########################################################################!
 
  MODULE const_and_precisions
- use itm_types, only : wp_ => r8
- use itm_constants, only : pi => itm_pi, e_ => itm_qe, me_ => itm_me, c_ => itm_c
 !########################################################################!
  IMPLICIT NONE
  PUBLIC
+
 !------------------------------------------------------------------------
 ! common precisions
 !------------------------------------------------------------------------
-! INTEGER, PARAMETER :: sp_   = 4         ! single precision
-! INTEGER, PARAMETER :: dp_   = 8         ! double precision
-! INTEGER, PARAMETER :: wp_   = dp_       ! work-precision
-! INTEGER, PARAMETER :: odep_ = dp_       ! ODE-solver precision
+
+!  INTEGER, PARAMETER :: I1 = SELECTED_INT_KIND (2)        ! Integer*1
+!  INTEGER, PARAMETER :: I2 = SELECTED_INT_KIND (4)        ! Integer*2
+  INTEGER, PARAMETER :: I4 = SELECTED_INT_KIND (9)        ! Integer*4
+!  INTEGER, PARAMETER :: I8 = SELECTED_INT_KIND (18)       ! Integer*8
+!  INTEGER, PARAMETER :: R4 = SELECTED_REAL_KIND (6, 37)   ! Real*4
+  INTEGER, PARAMETER :: R8 = SELECTED_REAL_KIND (15, 300) ! Real*8
+  INTEGER, PARAMETER :: wp_   = R8        ! work-precision
+! INTEGER, PARAMETER :: odep_ = R8        ! ODE-solver precision
 ! INTEGER, PARAMETER :: xp_   = wp_       ! for ext. modules if necessary
 !------------------------------------------------------------------------
 ! precisions which are in use in CONFIG_yat
@@ -28,7 +32,7 @@
 !========================================================================
   REAL(wp_), PARAMETER :: zero    = 0.0_wp_
   REAL(wp_), PARAMETER :: unit    = 1.0_wp_
-!  REAL(wp_), PARAMETER :: pi      = 3.141592653589793_wp_
+  real (kind = R8), parameter :: pi = 3.141592653589793238462643383280_R8
 ! REAL(wp_), PARAMETER :: sqrt_pi = 1.772453850905516_wp_
 ! REAL(wp_), PARAMETER :: sqrt_2  = 1.414213562373095_wp_
 ! REAL(wp_), PARAMETER :: rad     = pi/180.0_wp_
@@ -47,7 +51,7 @@
 !========================================================================
 ! Computer constants
 !========================================================================
-  REAL(wp_), PARAMETER :: comp_eps      = EPSILON(unit)
+  REAL(kind = R8), PARAMETER :: comp_eps      = EPSILON(unit)
 !  REAL(wp_), PARAMETER :: comp_eps2     = comp_eps**2
 ! REAL(wp_), PARAMETER :: comp_tiny     = TINY(unit)
 ! REAL(wp_), PARAMETER :: comp_huge     = HUGE(unit)
@@ -65,16 +69,24 @@
 !========================================================================
 ! Physical constants (SI)
 !========================================================================
-!  REAL(wp_), PARAMETER :: e_     = 1.602176487d-19    ! [C]
-!  REAL(wp_), PARAMETER :: me_    = 9.10938215d-31    ! [kg]
-! REAL(wp_), PARAMETER :: mp_    = 1.672621637d-27    ! [kg]
+  real (kind = R8), parameter :: e_ = 1.602176487e-19_R8  ! [C]
+  real (kind = R8), parameter :: me_ = 9.10938215e-31_R8      ! electron mass [kg]
+!  real (kind = R8), parameter :: mp_ = 1.672621637e-27_R8     ! proton mass [kg]
 ! REAL(wp_), PARAMETER :: rmpe_  = mp_/me_
-!  REAL(wp_), PARAMETER :: c_     = 2.99792458d+08    ! [m/s]
-! REAL(wp_), PARAMETER :: eps0_  = 8.854187817d-12    ! [F/m]
+!  real (kind = R8), parameter :: md_ = 3.34358320e-27_R8      ! deuteron mass, kg
+!  real (kind = R8), parameter :: mt_ = 5.00735588e-27_R8      ! triton mass, kg
+!  real (kind = R8), parameter :: ma_ = 6.64465620e-27_R8      ! alpha mass, kg
+!  real (kind = R8), parameter :: amu_ = 1.660538782e-27_R8    ! amu, kg
+  real (kind = R8), parameter :: c_ = 2.99792458e8_R8  ! speed of light [m/s]
+  real (kind = R8), parameter :: mu0_ = 4.0e-7_R8 * pi
+!  real (kind = R8), parameter :: eps0_ = 1.0_R8 / (mu0_ * c_ * c_)    ! [F/m]
+!  real (kind = R8), parameter :: avogr_ = 6.02214179e23_R8
+!  real (kind = R8), parameter :: KBolt_ = 1.3806504e-23_R8
 !------------------------------------------------------------------------
 ! Useful definitions
 !------------------------------------------------------------------------
-  REAL(wp_), PARAMETER :: keV_   = 1000*e_         ! [J]
+  real (kind = R8), parameter :: ev_ = e_          ! [J]
+  REAL(wp_), PARAMETER :: keV_   = 1000*ev_        ! [J]
   REAL(wp_), PARAMETER :: mc2_SI = me_*c_**2       ! [J]
   REAL(wp_), PARAMETER :: mc2_   = mc2_SI/keV_     ! [keV]
 ! REAL(wp_), PARAMETER :: mc_    = me_*c_          ! [kg*m/s]

src/gray-single.f90

@@ -0,0 +1,306 @@
+
+!      (ijetto, mr, mz, mrd, r, z, psin, psiax, psibnd,
+!     .    rax, zax, nbnd, rbnd, zbnd, nrho, psijet, f, te, dne,
+!     .    zeff, qsf, nbeam, powin, alphin, betain, dpdv, jcd, pec, 
+!     .    icd, ierr)
+
+  implicit real*8 (a-h,o-z)
+  parameter(pi=3.14159265358979d0)
+!  input arguments
+  integer :: ijetto, mr, mz, nbnd, nrho, nbeam
+  real*8, DIMENSION(:), ALLOCATABLE :: r, z, psinr
+  real*8, DIMENSION(:,:), ALLOCATABLE :: psin, psi
+  real*8 psiax, psibnd, rax, zax
+  real*8, DIMENSION(:), ALLOCATABLE :: rbnd, zbnd
+  real*8, DIMENSION(:), ALLOCATABLE :: psijet, f, qsf, te, dne, zeff
+  real*8, DIMENSION(:), ALLOCATABLE :: teitpl, dneitpl, zeffitpl, cte, ctdn, ctz
+  real*8, DIMENSION(:), ALLOCATABLE :: powin, alphain, betain
+!  do loop arguments
+  real*8 :: powloop,alphaloop,betaloop
+!  do loop output arguments
+  real*8, DIMENSION(:), ALLOCATABLE :: dpdv, jcd
+  real*8 :: pec, icd
+!  gray_main output arguments
+  real*8, DIMENSION(:), ALLOCATABLE :: dpdvloop, jcdloop
+  real*8 :: pecloop, icdloop
+  integer :: ierr
+!  unused read variales
+  real*8, DIMENSION(:), ALLOCATABLE :: pres,ffprim,pprim,rlim,zlim
+  character(LEN=48) :: stringa
+!!!!!!!!!!  VARIABILI PER TEST  !!!!!!!!!!
+  real*8,dimension(2) :: alphalim, betalim
+  real*8 :: delta
+  integer :: alphanum, betanum, len
+  character(LEN=48) :: profile, eqfile
+!!!!!!!!!!  VARIABILI PER TEST  !!!!!!!!!!
+!  
+!  local variables
+!     real*8 fgray(nrho), qgray(nrho), jcdgry(nrho), icdgry
+!
+!  input arguments
+!
+!    ijetto  Equilibrium source (1 EFIT, 2 ESCO)
+!            If IJETTO=2, then PSIN values are valid only inside
+!            plasma boudary (PSIN=0 outside)
+!    mr      Size of flux map grid in R direction
+!    mz      Size of flux map grid in Z direction
+!    mrd     Leading dimension of the psin(:,:) array, mrd>mr
+!    r       R coordinates of flux map grid points [m]
+!    z       Z coordinates of flux map grid points [m]
+!    psin    Normalised poloidal flux psin=(psi-psiax)/(psibnd-psiax)
+!            on the (R, Z) grid.
+!    psiax   Poloidal flux on axis [Wb rad-1]
+!    psibnd  Poloidal flux on boundary [Wb rad-1]
+!    rax     R coordinate of magnetic axis [m]
+!    zax     Z coordinate of magnetic axis [m]
+!    nbnd    Number of points in plasma boundary contour
+!    rbnd    R coordinates of plasma boundary contour [m]
+!    zbnd    Z coordinates of plasma boundary contour [m]
+!
+!    nrho    Number of points in JETTO rho grid    -
+!    psijet  Normalised poloidal flux on JETTO radial grid
+!    f       Poloidal current stream function f=Bphi*R on JETTO
+!            radial grid [T m]
+!    te      Electron temperature on JETTO radial grid [eV]
+!    dne     Electron density on JETTO radial grid [m-3]
+!    zeff    Effective nuclear charge Zeff on JETTO radial grid
+!    qsf     Safety factor on JETTO radial grid
+!
+!    nbeam   number of beams injected
+!    powin   Input ECRH power array [W] (powin(i) =< 0 means i-th beam is unused)
+!    alphin  Beams poloidal injection angles array [rad]
+!    betain  Beams toroidal injection angles array [rad]
+!
+!  output arguments
+!
+!    dpdv    Absorbed EC power density on JETTO radial grid [W m-3]
+!    jcd     EC driven flux averaged current density on JETTO
+!            radial grid [A m-2]
+!    pec     Total absorbed EC power [W]
+!    icd     Total EC driven current [A]
+!    ierr    Return code. IERR>0 on error
+!            ierr = 90-93: error computing integrals for current drive
+!            ierr = 94: absorption coefficient alpha < 0
+!            ierr = 97: parallel comp. refract. idx N//>0.99 (warning)
+!            ierr = 98: parallel comp. refract. idx N//>1.05
+!
+!  JETTO coordinate system assumes toroidal angle increasing CW
+!  in GRAY toroidal angle increases CCW --> adapt signs on input data
+  
+!!!!!!!!!!  VARIABILI PER TEST  !!!!!!!!!!
+
+  allocate(powin(1),alphain(1),betain(1))
+  nbeam=1
+
+  ijetto=1
+  
+!  lettura gray.data
+  open(99,file='gray.data')
+  read(99,*) alphain(1),betain(1)
+  read(99,*)
+  read(99,*) powin(1)
+  powin(1)=powin(1)*1.d6
+  do i=1,13
+    read(99,*)
+  end do
+  read(99,*) eqfile
+  read(99,*) ipsinorm
+  read(99,*) factb
+  read(99,*) profile
+  read(99,*)
+  read(99,*) sgnbphi,sgniphi,icocos
+  close(99)
+  
+!  lettura profili di flusso poloidale, temperatura, densità e carica nucleare
+  len=length(profile)
+  open(98,file=profile(1:len)//'.prf',action='read')
+  read(98,*) nrho
+  allocate(psijet(nrho),te(nrho),dne(nrho),zeff(nrho))
+  read(98,*) psijet0,te0,dne0,zeff0
+  if(psijet0.ne.0.0d0) psijet0=0.0d0
+  psijet(1)=psijet0
+  te(1)=te0
+  dne(1)=dne0
+  zeff(1)=zeff0
+  do i=2,nrho
+    read(98,*) psijeti,tei,dnei,zeffi !,rhoi,qsfi
+    psijet(i)=psijeti
+    te(i)=tei
+    dne(i)=dnei
+    zeff(i)=zeffi
+  end do
+  close(98)
+  
+!  lettura dati campo magnetico
+  len=length(eqfile)
+  open(99,file=eqfile(1:len)//'.eqdsk',action='read')
+
+  if(icocos.gt.0) then 
+   read (99,'(a48,3i4)') stringa,idum,mr,mz
+  else
+   read (99,*) mr,mz
+  end if
+  
+  allocate(r(mr),z(mz),psinr(mr),psi(mr,mz),psin(mr,mz),f(mr),qsf(mr))
+  allocate(pres(mr),ffprim(mr),pprim(mr))
+  allocate(dpdv(mr),jcd(mr),dpdvloop(mr), jcdloop(mr))
+  
+  if(icocos.gt.0) then 
+   read (99,'(5e16.9)') drnr1,dznz1,rcen,rleft,zmid
+   read (99,'(5e16.9)') rax,zax,psiax,psibnd,btrcen
+   read (99,'(5e16.9)') current,xdum,xdum,xdum,xdum
+   read (99,'(5e16.9)') xdum,xdum,xdum,xdum,xdum
+   read (99,'(5e16.9)') (f(i),i=1,mr)
+   read (99,'(5e16.9)') (pres(i),i=1,mr)
+   read (99,'(5e16.9)') (ffprim(i),i=1,mr)
+   read (99,'(5e16.9)') (pprim(i),i=1,mr)
+   if (ipsinorm.eq.0) then
+     read (99,'(5e16.9)') ((psi(i,j),i=1,mr),j=1,mz)
+   else
+     read (99,'(5e16.9)') ((psin(i,j),i=1,mr),j=1,mz)
+   end if    
+   read (99,'(5e16.9)') (qsf(i),i=1,mr)
+   read (99,*) nbnd,nlim
+   allocate(rbnd(nbnd),zbnd(nbnd),rlim(nlim),zlim(nlim))
+   read (99,'(5e16.9)') (rbnd(i),zbnd(i),i=1,nbnd)
+   read (99,'(5e16.9)') (rlim(i),zlim(i),i=1,nlim)
+  else
+   read (99,*) drnr1,dznz1,rcen,rleft,zmid
+   read (99,*) rax,zax,psiax,psibnd,btrcen
+   read (99,*) current,xdum,xdum,xdum,xdum
+   read (99,*) xdum,xdum,xdum,xdum,xdum
+   read (99,*) (f(i),i=1,mr)
+   read (99,*) (pres(i),i=1,mr)
+   read (99,*) (ffprim(i),i=1,mr)
+   read (99,*) (pprim(i),i=1,mr)
+   if (ipsinorm.eq.0) then
+     read (99,*) ((psi(i,j),i=1,mr),j=1,mz)
+   else
+     read (99,*) ((psin(i,j),i=1,mr),j=1,mz)
+   end if    
+   read (99,*) (qsf(i),i=1,mr)
+   read (99,*) nbnd,nlim
+   allocate(rbnd(nbnd),zbnd(nbnd),rlim(nlim),zlim(nlim))
+   read (99,*) (rbnd(i),zbnd(i),i=1,nbnd)
+   read (99,*) (rlim(i),zlim(i),i=1,nlim)
+  end if
+  close(99)
+  
+!  adeguamento alle convenzioni e normalizzazioni
+  icocosmod=mod(icocos,10)
+
+!  icocos mod 2 = 0: toroidal angle phi CW (opposite to gray convention)
+  if (mod(icocos,2).eq.0.and.icocosmod.gt.0) then
+    btrcen=-btrcen
+    current=-current
+    do i=1,mr
+      f(i)=-f(i)
+    end do
+  end if
+
+  if (icocosmod.ne.0) sgnbphi=sign(1.0d0,f(mr))
+  btrcen=sign(btrcen,sgnbphi)
+
+  do i=1,mr
+    f(i)=sign(f(i),sgnbphi)*factb
+  end do
+  
+!  icocos mod 10 = 1,4,5,8: psi increasing with CCW Ip
+!  icocos mod 10 = 2,3,6,7: psi decreasing with CCW Ip
+  if (icocosmod.eq.1 .or. icocosmod.eq.4 .or.icocosmod.eq.5 .or. icocosmod.eq.8) then
+    psibnd=-psibnd
+    psiax=-psiax
+    if (ipsinorm.eq.0) then
+      do j=1,mz
+        do i=1,mr
+          psi(i,j)=-psi(i,j)
+        end do
+      end do
+    end if
+  end if
+  
+  psia0=psibnd-psiax
+  
+!  icocos=0: adapt psi to force Ip sign, otherwise maintain psi
+  if (icocosmod.ne.0) sgniphi=sign(1.0d0,-psia0)
+  current=sign(current,sgniphi)
+  psia=sign(psia0,-sgniphi)*factb
+  
+!  icocos>10: input psi is in Wb
+!  icocos<10: input psi is in Wb/rad (gray convention)
+  if (icocos.ge.10) psia=psia/(2.0d0*pi)
+  
+  if(ipsinorm.eq.0) then
+    do j=1,mz
+      do i=1,mr
+        psin(i,j)=(psi(i,j)-psiax)/(psia0)
+        psi(i,j)=psin(i,j)*psia
+      enddo
+    enddo
+  else
+    do j=1,mz
+      do i=1,mr
+        psi(i,j)=psin(i,j)*psia
+      enddo
+    enddo
+  end if
+  
+  dr=drnr1/dble(mr-1)
+  dz=dznz1/dble(mz-1)
+  r(1)=rleft
+  z(1)=zmid-dznz1/2.0d0
+  dpsinr=1.0d0/dble(mr-1)
+  
+  do i=1,mr
+    psinr(i)=dpsinr*(i-1)
+    r(i)=r(1)+(i-1)*dr
+  end do
+  
+  do j=1,mz
+    z(j)=z(1)+(j-1)*dz
+  end do
+  
+!  interpolazione Te, dne, Zeff (psijet -> psinr, nrho -> mr elementi)
+  
+  allocate(teitpl(mr),dneitpl(mr),zeffitpl(mr))
+  allocate(cte(4*nrho),ctdn(4*nrho),ctz(4*nrho))
+  
+  iopt=0
+  call difcsg(psijet,te,nrho,iopt,cte,ier)
+  call difcsg(psijet,dne,nrho,iopt,ctdn,ier)
+  call difcsg(psijet,zeff,nrho,iopt,ctz,ier)
+  
+  do i=1,mr
+    call vlocate(psijet,nrho,psinr(i),k)
+    k = max(1,min(k,nrho-1))
+    dal = psinr(i)-psijet(k)
+
+    teitpl(i) = fspli(cte,nrho,k,dal)*1000.d0
+    dneitpl(i) = fspli(ctdn,nrho,k,dal)*1.d19
+    zeffitpl(i) = fspli(ctz,nrho,k,dal)
+  end do
+  
+!  inizializzazione variabili output
+  dpdv = 0
+  jcd = 0
+  pec = 0
+  icd = 0
+  
+  do j=1,nbeam
+!  potenza ECRH e angoli per fascio j-esimo
+    powloop = powin(j)
+    alphaloop = alphain(j)
+    betaloop = betain(j)
+    
+    call gray_main(ijetto, mr, mz, r, z, psin, 0.0d0, psia, &
+     rax, zax, nbnd, rbnd, zbnd, mr, psinr, f, teitpl, dneitpl, &
+     zeffitpl, qsf, j, powloop, alphaloop, betaloop, &
+     dpdvloop, jcdloop, pecloop, icdloop, ierr)
+    
+    dpdv = dpdv + dpdvloop
+    jcd = jcd + jcdloop
+    pec = pec + pecloop
+    icd = icd + icdloop
+  end do
+  
+end

src/gray_main.f90

@@ -0,0 +1,105 @@
+subroutine gray_main(ijetto, mr, mz, r, z, psin, psiax, psibnd, &
+    rax, zax, nbnd, rbnd, zbnd, nrho, psijet, f, te, dne, zeff, qsf, &
+    beamid, powin, alphain, betain, dpdv, jcd, pec, icd, ierr)
+
+  use const_and_precisions, only : r8
+  use graydata_flags, only : igrad, ipass
+  
+  implicit none
+
+  integer, intent(in) :: ijetto, mr, mz, nrho, nbnd, beamid
+  real(r8), intent(in) :: r(mr), z(mz), psin(mr,mz)
+  real(r8), intent(in) :: psiax, psibnd, rax, zax
+  real(r8), intent(in), dimension(nbnd) :: rbnd, zbnd
+  real(r8), intent(in), dimension(nrho) :: psijet, f, qsf, te, dne, zeff
+  real(r8), intent(in) :: powin, alphain, betain
+  real(r8), intent(out), dimension(nrho) :: dpdv, jcd
+  real(r8), intent(out) :: pec, icd
+  integer, intent(out) :: ierr
+
+  integer :: istop, iercom, iopmin, iowmin, index_rt
+  real(r8) :: sox, p0mw, powrfl, taumn, taumx, pabstot, currtot
+  real(r8) :: p0mw1, powtr, pabstott, currtott
+  real(r8), dimension(nrho) :: rhopol,dpdv1pass, jcd1pass
+
+  common/istop/istop
+  common/ierr/iercom
+  common/iovmin/iopmin,iowmin
+  common/mode/sox
+  common/p0/p0mw
+  common/powrfl/powrfl
+  common/index_rt/index_rt
+  common/taumnx/taumn,taumx,pabstot,currtot
+  
+! read data plus initialization
+  index_rt=1
+  call prfile
+  call paraminit
+  call read_data(ijetto, mr, mz, r, z, psin, psiax, psibnd, rax, zax, &
+    nbnd, rbnd, zbnd, nrho, psijet, f, te, dne, zeff, qsf, powin, alphain, betain, &
+    beamid)
+  call vectinit
+  if(iercom.eq.0) then
+    if(igrad.eq.0) call ic_rt
+    if(igrad.gt.0) call ic_gb
+  end if
+  if(iercom.gt.0) then
+    ierr=iercom
+    write(*,*) '    IERR = ', ierr
+    return
+  end if
+
+! beam/ray propagation
+  call gray_integration
+
+! postprocessing
+  rhopol=sqrt(psijet)
+  call after_gray_integration(rhopol,nrho,dpdv1pass,jcd1pass)
+
+  pabstott=pabstot
+  currtott=currtot
+  dpdv=dpdv1pass
+  jcd=jcd1pass
+  powtr=p0mw-pabstot
+
+  if (iowmin==2 .and. ipass>1) then
+! second pass into plasma
+    p0mw1=p0mw
+    igrad=0
+    
+    index_rt=2
+    p0mw=p0mw1*powrfl
+    call prfile
+    call vectinit2
+    call paraminit
+    call ic_rt2
+    call gray_integration
+    call after_gray_integration(rhopol,nrho,dpdv1pass,jcd1pass)
+    pabstott=pabstott+pabstot
+    currtott=currtott+currtot
+    dpdv=dpdv+dpdv1pass
+    jcd=jcd+jcd1pass
+  
+    index_rt=3
+    sox=-sox
+    p0mw=p0mw1*(1._r8-powrfl)
+    call prfile
+    call vectinit2
+    call paraminit
+    call ic_rt2
+    call gray_integration
+    call after_gray_integration(rhopol,nrho,dpdv1pass,jcd1pass)
+    pabstott=pabstott+pabstot
+    currtott=currtott+currtot
+    dpdv=dpdv+dpdv1pass
+    jcd=jcd+jcd1pass
+  end if
+  write(*,*)
+  write(*,*) 'Pabs (MW), Icd (kA) = ', pabstott,currtott*1.e3_r8
+  call vectfree
+  pec=pabstott*1.e6_r8
+  icd=currtott*1.e6_r8
+  dpdv=dpdv*1.e6_r8
+  jcd=jcd*1.e6_r8
+  ierr=iercom
+end subroutine gray_main

src/graydata_anequil.f90

@@ -0,0 +1,11 @@
+module graydata_anequil
+  use const_and_precisions, only : r8
+  implicit none
+  
+  real(r8), save :: dens0,aln1,aln2
+  real(r8), save :: te0,dte0,alt1,alt2
+  real(r8), save :: rr0,zr0,rpa,rr0m,zr0m,rpam
+  real(r8), save :: b0,q0,qa,alq
+  real(r8), save :: zeffan
+
+end module graydata_anequil

src/graydata_flags.f90

@@ -0,0 +1,13 @@
+module graydata_flags
+  use const_and_precisions, only : r8
+  implicit none
+  
+  character*255, save :: filenmeqq,filenmprf,filenmbm
+  real(r8), save :: sspl, dst
+  integer, save :: ibeam,iequil,ixp,iprof
+  integer, save :: iwarm,ilarm,imx,ieccd,ipec,idst
+  integer, save :: igrad,ipass
+  integer, save :: ipsinorm,iscal,icocos
+  integer, save :: nnd,istpr0,istpl0
+  
+end module graydata_flags

src/graydata_par.f90

@@ -0,0 +1,10 @@
+module graydata_par
+  use const_and_precisions, only : r8
+  implicit none
+  
+  real(r8), save :: rwmax,rwallm
+  real(r8), save :: psipol0,chipol0
+  real(r8), save :: factb,factt,factn,psdbnd
+  real(r8), save :: sgnbphi,sgniphi
+
+end module graydata_par

src/itm_constants.f90

@@ -1,32 +0,0 @@
-!> Module implementing the ITM physics constants
-!>
-!> Source:
-!>  based on SOLPS b2mod_constants.F
-!>    '09/12/07 xpb : source CODATA 2006 (http://www.nist.gov/)'
-!>  pulled from ets r100
-!>
-!> \author David Coster
-!>
-!> \version "$Id: itm_constants.f90 37 2009-08-17 17:15:00Z coster $"
-
-module itm_constants
-
-  use itm_types
-
-  real (kind = R8), parameter :: itm_pi = 3.141592653589793238462643383280_R8
-  real (kind = R8), parameter :: itm_c = 2.99792458e8_R8         ! speed of light, m/s
-  real (kind = R8), parameter :: itm_me = 9.10938215e-31_R8      ! electron mass, kg
-  real (kind = R8), parameter :: itm_mp = 1.672621637e-27_R8     ! proton mass, kg
-  real (kind = R8), parameter :: itm_md = 3.34358320e-27_R8      ! deuteron mass, kg
-  real (kind = R8), parameter :: itm_mt = 5.00735588e-27_R8      ! triton mass, kg
-  real (kind = R8), parameter :: itm_ma = 6.64465620e-27_R8      ! alpha mass, kg
-  real (kind = R8), parameter :: itm_amu = 1.660538782e-27_R8    ! amu, kg
-  real (kind = R8), parameter :: itm_ev = 1.602176487e-19_R8
-  real (kind = R8), parameter :: itm_qe = itm_ev
-  real (kind = R8), parameter :: itm_mu0 = 4.0e-7_R8 * itm_pi
-  real (kind = R8), parameter :: itm_eps0 = 1.0_R8 / (itm_mu0 * itm_c * itm_c)
-  real (kind = R8), parameter :: itm_avogr = 6.02214179e23_R8
-  real (kind = R8), parameter :: itm_KBolt = 1.3806504e-23_R8
-  character (len=64), parameter :: itm_constants_version = '$Id: itm_constants.f90 37 2009-08-17 17:15:00Z coster $'
-
-end module itm_constants

src/itm_types.f90

@@ -1,50 +0,0 @@
-!> Module implementing the ITM basic types
-!>
-!> Source:
-!>  based on SOLPS b2mod_types.F
-!>  pulled from ets r100 and extended with input from C. Konz, T. Ribeiro & B. Scott
-!>
-!> \author David Coster
-!>
-!> \version "$Id: itm_types.f90 144 2010-10-07 09:26:24Z konz $"
-
-module itm_types
-
-  INTEGER,  PARAMETER :: ITM_I1 = SELECTED_INT_KIND (2)        ! Integer*1
-  INTEGER,  PARAMETER :: ITM_I2 = SELECTED_INT_KIND (4)        ! Integer*2
-  INTEGER,  PARAMETER :: ITM_I4 = SELECTED_INT_KIND (9)        ! Integer*4
-  INTEGER,  PARAMETER :: ITM_I8 = SELECTED_INT_KIND (18)       ! Integer*8
-  INTEGER,  PARAMETER :: R4 = SELECTED_REAL_KIND (6, 37)   ! Real*4
-  INTEGER,  PARAMETER :: R8 = SELECTED_REAL_KIND (15, 300) ! Real*8
-
-  INTEGER,  PARAMETER :: itm_int_invalid = -999999999
-  REAL(R8), PARAMETER :: itm_r8_invalid = -9.0D40
-
-  interface itm_is_valid
-     module procedure itm_is_valid_int4, itm_is_valid_int8, itm_is_valid_real8
-  end interface
-
-contains
-
-  logical function itm_is_valid_int4(in_int)
-    implicit none
-    integer(ITM_I4) in_int
-    itm_is_valid_int4 = in_int .ne. itm_int_invalid
-    return
-  end function itm_is_valid_int4
-
-  logical function itm_is_valid_int8(in_int)
-    implicit none
-    integer(ITM_I8) in_int
-    itm_is_valid_int8 = in_int .ne. itm_int_invalid
-    return
-  end function itm_is_valid_int8
-
-  logical function itm_is_valid_real8(in_real)
-    implicit none
-    real(R8) in_real
-    itm_is_valid_real8 = abs(in_real - itm_r8_invalid) .gt. abs(itm_r8_invalid) * 1.0e-15_R8
-    return
-  end function itm_is_valid_real8
-
-end module itm_types

src/main.f90

@@ -0,0 +1,286 @@
+program main
+  use const_and_precisions, only : r8, pi
+  implicit none
+  integer, parameter :: cocosout=3
+  integer :: u, ierr, cocos, exp2pi, exp2piout, idum, i, j, ipsinorm
+  logical :: ispsinorm
+  logical :: phiccw, psiincr, qpos, phiccwout, psiincrout, qposout
+  character(len=255) :: desc, eqdskfile, prffile
+  character(len=*), parameter :: fmt2000='(a48,3i4)',fmt2020='(5e16.9)',fmt2022='(2i5)'
+  integer :: nr, nz, nbnd, nprf
+  real(r8) :: deltar, deltaz, r0, r1, zmid, rax, zax, psiax, psib, bref, ipla, &
+      dummy, z1, dr, dz, dpsin, bref0, ipla0, qb0, psifact, psii, powin, pec, icd
+  real(r8), dimension(:,:), allocatable :: psi2d
+  real(r8), dimension(:),   allocatable :: r, z, psi1d, q, fdia, rbnd, zbnd, &
+      psiprf, te, dne, zeff, cte, cne, czeff, teeq, dneeq, zeffeq, dpdv, jcd
+  real(r8), external :: fspli
+
+! === read filenames and input power from gray.data ===
+  u = get_free_unit()
+  open(u,file='gray.data',status= 'unknown')
+  read(u,*)
+  read(u,*)
+  read(u,*) powin
+  do i=1,13
+    read(u,*)
+  end do
+  read(u,*) desc
+  eqdskfile=trim(desc)//'.eqdsk'
+  read(u,*) ipsinorm
+  ispsinorm=(ipsinorm==1)
+  read(u,*)
+  read(u,*) desc
+  prffile=trim(desc)//'.prf'
+  read(u,*)
+  read(u,*) dummy,dummy,cocos
+  close(u)
+
+! === read EQDSK ===
+  u = get_free_unit()
+  open(unit=u,file=eqdskfile,status='OLD',action='READ',iostat=ierr)
+  if(ierr/=0) then
+    write(*,*) 'Cannot open file '//trim(eqdskfile)
+    stop
+  end if
+
+! read grid size and allocate arrays
+  read (u,fmt2000) desc,idum,nr,nz
+  allocate(psi2d(nr,nz), r(nr), z(nz), psi1d(nr), q(nr), fdia(nr), stat=ierr)
+  if (ierr/=0) then
+    close(u)
+    call free_allocs
+    write(*,*) 'cannot allocate arrays for equilibrium data'
+    stop
+  end if
+! read 0D fields
+  read (u,fmt2020) deltar, deltaz, r0, r1, zmid
+  read (u,fmt2020) rax, zax, psiax, psib, bref
+  read (u,fmt2020) ipla, (dummy,i=1,4)
+  read (u,fmt2020) (dummy,i=1,5)
+! read 1D-2D fields
+  read (u,fmt2020) (fdia(i),i=1,nr)
+  read (u,fmt2020) (dummy,i=1,nr)
+  read (u,fmt2020) (dummy,i=1,nr)
+  read (u,fmt2020) (dummy,i=1,nr)
+  read (u,fmt2020) ((psi2d(i,j),i=1,nr),j=1,nz)
+  read (u,fmt2020) (q(i),i=1,nr)
+! read boundary size and allocate arrays
+  read (u,fmt2022) nbnd
+  if (nbnd>0) then
+    allocate(rbnd(nbnd), zbnd(nbnd), stat=ierr)
+    if (ierr/=0) then
+      close(u)
+      call free_allocs
+      write(*,*) 'cannot allocate arrays for boundary data'
+      stop
+    end if
+! read boundary shape
+    read (u,fmt2020) (rbnd(i), zbnd(i),i=1,nbnd)
+  end if
+  close(u)
+
+! normalize psi2d
+  if (.not.ispsinorm) psi2d = (psi2d - psiax)/(psib - psiax)
+
+! interpret cocos numbers
+  call decode_cocos(cocos,exp2pi,phiccw,psiincr,qpos)
+  call decode_cocos(cocosout,exp2piout,phiccwout,psiincrout,qposout)
+
+! check sign consistency
+  ipla0=ipla
+  bref0=bref
+  if (psiincr) then
+    ipla=sign(ipla,psib - psiax)
+  else
+    ipla=sign(ipla,psiax - psib)
+  end if
+  bref=sign(bref,fdia(nr))
+  if (ipla/ipla0<0 .or. bref/bref0<0) then
+    write(*,*) 'Warning: sign inconsistency in Ipla/psi or Bref/Fdia'
+  end if
+  qb0=q(nr)
+  if (qpos) then
+    q=sign(q,ipla*bref)
+  else
+    q=sign(q,-ipla*bref)
+  end if
+  if (q(nr)/qb0<0) then
+    write(*,*) 'Warning: sign inconsistency found among q, Ipla and Bref'
+  end if
+
+! convert cocosin to cocosout
+  if (phiccw.neqv.phiccwout) then
+!   opposite direction of toroidal angle phi in cocosin and cocosout
+    bref=-bref
+    ipla=-ipla
+    fdia=-fdia
+  end if
+  if ((phiccw.eqv.phiccwout) .neqv. (psiincr.eqv.psiincrout)) then
+!   psi and Ip signs don't change accordingly
+    psib=-psib
+    psiax=-psiax
+  end if
+  if (qpos .neqv. qposout) q=-q
+!   q has opposite sign for given sign of Bphi*Ip
+  if (exp2pi/=exp2piout) then
+!   convert Wb to Wb/rad or viceversa
+    psifact=(2._r8*pi)**(exp2piout-exp2pi)
+    psib=psib*psifact
+    psiax=psiax*psifact
+  end if
+
+! fill equi-spaced R, z, psi arrays
+  z1 = zmid - deltaz/2._r8
+  dr = deltar/(nr-1)
+  dz = deltaz/(nz-1)
+  dpsin = 1._r8/(nr-1)
+  do i=1,nr
+    r(i) = r1 + (i-1)*dr
+    psi1d(i) = (i-1)*dpsin
+  end do
+  do j=1,nz
+    z(j) = z1 + (j-1)*dz
+  end do
+   
+! === read profiles ===
+  u = get_free_unit()
+  open(unit=u,file=prffile,status='OLD',action='READ',iostat=ierr)
+  if(ierr/=0) then
+    write(*,*) 'Cannot open file '//trim(prffile)
+    stop
+  end if
+  read(u,*) nprf
+  if (nprf>0) then
+    allocate(psiprf(nprf),te(nprf),dne(nprf),zeff(nprf), stat=ierr)
+    if (ierr==0) allocate(cte(4*nprf),cne(4*nprf),czeff(4*nprf), stat=ierr)
+    if (ierr/=0) then
+      close(u)
+      call free_allocs
+      write(*,*) 'cannot allocate arrays for input 1d profiles'
+      stop
+    end if
+    do i=1,nprf
+      read(u,*) psiprf(i),te(i),dne(i),zeff(i)
+    end do
+  else
+    write(*,*) 'no data for 1d profiles'
+    stop
+  end if
+  close(u)
+
+! spline interpolation for resampling on uniform grid
+  call difcsg(psiprf,te  ,nprf,0,cte  ,ierr)
+  if (ierr==0) call difcsg(psiprf,dne ,nprf,0,cne  ,ierr)
+  if (ierr==0) call difcsg(psiprf,zeff,nprf,0,czeff,ierr)
+  if (ierr/=0) then
+    call free_allocs
+    write(*,*) 'error in input 1d profiles interpolation'
+    stop
+  end if
+  allocate(teeq(nr),dneeq(nr),zeffeq(nr),dpdv(nr),jcd(nr), stat=ierr)
+  if (ierr/=0) then
+    call free_allocs
+    write(*,*) 'cannot allocate arrays for resampled and output 1d profiles'
+    stop
+  end if
+  do i=1,nr
+    psii=psi1d(i)
+    call vlocate(psiprf,nprf,psii,j)
+    j=max(1,min(j,nprf-1))
+    dpsin=psii-psiprf(j)
+    teeq(i)  =fspli(cte,  nprf,j,dpsin)
+    dneeq(i) =fspli(cne,  nprf,j,dpsin)
+    zeffeq(i)=fspli(czeff,nprf,j,dpsin)
+  end do
+
+! convert keV to eV, 10^19 m^-3 to m^-3, and MW to W
+  teeq=teeq*1.e3_r8
+  dneeq=dneeq*1.e19_r8
+  powin=powin*1.e6_r8
+
+! === call GRAY subroutine ===
+  call gray_main(1, nr, nz, r, z, psi2d, psiax, psib, &
+    rax, zax, nbnd, rbnd, zbnd, nr, psi1d, fdia, teeq, dneeq, zeffeq, q, &
+    powin, dpdv, jcd, pec, icd, ierr)
+
+  call free_allocs
+  
+contains
+
+  function get_free_unit(umin,umax) result(i)
+    implicit none
+    integer :: i
+    integer, intent(in), optional :: umin, umax
+    integer, parameter :: max_allowed = 99
+    integer :: ierr, iend
+    logical :: ex, op
+
+    if (present(umin)) then
+      i = max(0,umin)    ! start searching from unit min
+    else
+      i = 0
+    end if
+    if (present(umax)) then
+      iend = min(max(0,umax),max_allowed)
+    else
+      iend = max_allowed
+    end if
+    do
+      if (i>iend) then
+        i=-1  ! no free units found
+        exit
+      end if
+      inquire(unit=i,exist=ex,opened=op,iostat=ierr)
+      if (ierr/=0) then
+        i=-2  ! cannot inquire unit i
+        exit
+      end if
+      if (ex.and..not.op) exit  ! unit i exists and is not open
+      i = i + 1
+    end do
+  end function get_free_unit
+
+  subroutine decode_cocos(cocos,exp2pi,phiccw,psiincr,qpos)
+    implicit none
+    integer, intent(in) :: cocos
+    integer, intent(out) :: exp2pi
+    logical, intent(out) :: phiccw,psiincr,qpos
+    integer :: cocosm10,cocosm4
+
+    cocosm10=mod(cocos,10)
+    cocosm4=mod(cocosm10,4)
+!  cocos>10 psi in Wb, cocos<10 psi in Wb/rad
+    exp2pi=(cocos-cocosm10)/10
+!  cocos mod 10 = 1,3,5,7: toroidal angle phi increasing CCW
+    phiccw=(mod(cocosm10,2)==1)
+!  cocos mod 10 = 1,2,5,6: psi increasing with positive Ip
+    psiincr=(cocosm4==1 .or. cocosm4==2)
+!  cocos mod 10 = 1,2,7,8: q positive for positive Bphi*Ip
+    qpos=(cocosm10<3 .or. cocosm10>6)
+  end subroutine decode_cocos
+
+  subroutine free_allocs
+    implicit none
+    if(allocated(psi2d))  deallocate(psi2d)
+    if(allocated(r))      deallocate(r)
+    if(allocated(z))      deallocate(z)
+    if(allocated(psi1d))  deallocate(psi1d)
+    if(allocated(q))      deallocate(q)
+    if(allocated(fdia))   deallocate(fdia)
+    if(allocated(rbnd))   deallocate(rbnd)
+    if(allocated(zbnd))   deallocate(zbnd)
+    if(allocated(psiprf)) deallocate(psiprf)
+    if(allocated(te))     deallocate(te)
+    if(allocated(dne))    deallocate(dne)
+    if(allocated(zeff))   deallocate(zeff)
+    if(allocated(cte))    deallocate(cte)
+    if(allocated(cne))    deallocate(cne)
+    if(allocated(czeff))  deallocate(czeff)
+    if(allocated(teeq))   deallocate(teeq)
+    if(allocated(dneeq))  deallocate(dneeq)
+    if(allocated(zeffeq)) deallocate(zeffeq)
+    if(allocated(dpdv))  deallocate(dpdv)
+    if(allocated(jcd)) deallocate(jcd)
+  end subroutine free_allocs
+
+end program main

src/reflections.f90

@@ -32,7 +32,8 @@ subroutine inters_linewall(xv,kv,rw,zw,nw,sint,normw)
   integer :: i,j,ni,iint
   real(r8), dimension(2) :: si,ti
   real(r8) :: drw,dzw,xint,yint,rint,l,kxy
-  real(r8) :: tol=sqrt(epsilon(1.0_r8))
+  real(r8) :: tol
+  tol=sqrt(epsilon(1.0_r8))
   sint=huge(sint)
   iint=0
   normw=0.0_r8
@@ -187,12 +188,17 @@ function inside(xc,yc,n,x,y)
   integer, dimension(n) :: jint
   real(r8), dimension(n) :: xint
   real(r8), dimension(n+1) :: xclosed,yclosed
-  integer :: i,nj
+  integer :: i,np,nj
   xclosed(1:n)=xc(1:n)
   yclosed(1:n)=yc(1:n)
-  xclosed(n+1)=xc(1)
-  yclosed(n+1)=yc(1)
-  call locate_unord(yclosed,n+1,y,jint,n,nj)
+  if (xc(1)==xc(n) .and. yc(1)==yc(n)) then
+    np=n
+  else
+    np=n+1
+    xclosed(np)=xc(1)
+    yclosed(np)=yc(1)
+  end if
+  call locate_unord(yclosed(1:np),np,y,jint,n,nj)
   inside=.false.
   if (nj==0) return
   do i=1,nj

src/scatterspl.f90

@@ -0,0 +1,49 @@
+subroutine scatterspl(x,y,z,w,n,kspl,sspl,xmin,xmax,ymin,ymax, &
+                        tx,nknt_x,ty,nknt_y,coeff,ierr)
+  use const_and_precisions, only : r8, comp_eps
+  implicit none
+  integer :: n
+  real(r8), dimension(n) :: x, y, z
+  real(r8), dimension(n) :: w
+  integer :: kspl
+  real(r8) :: sspl
+  real(r8) :: xmin, xmax, ymin, ymax
+  real(r8), dimension(nknt_x) :: tx
+  real(r8), dimension(nknt_y) :: ty
+  integer :: nknt_x, nknt_y
+  real(r8), dimension(nknt_x*nknt_y) :: coeff
+  integer :: ierr
+
+  integer :: iopt
+  real(r8) :: resid
+  integer :: u,v,km,ne,b1,b2,lwrk1,lwrk2,kwrk,nxest,nyest
+  real(r8), dimension(:), allocatable :: tx_tmp,ty_tmp
+  real(r8), dimension(:), allocatable :: wrk1, wrk2
+  integer, dimension(:), allocatable :: iwrk
+
+  nxest=nknt_x
+  nyest=nknt_y
+  ne = max(nxest,nyest)
+  allocate(tx_tmp(ne),ty_tmp(ne))
+
+  km = kspl+1
+  u = nxest-km
+  v = nyest-km
+  b1 = kspl*min(u,v)+kspl+1
+  b2 = (kspl+1)*min(u,v)+1
+  lwrk1 = u*v*(2+b1+b2)+2*(u+v+km*(n+ne)+ne-2*kspl)+b2+1
+  lwrk2 = u*v*(b2+1)+b2
+  kwrk = n+(nknt_x-2*kspl-1)*(nknt_y-2*kspl-1)
+  allocate(wrk1(lwrk1),wrk2(lwrk2),iwrk(kwrk))
+
+  iopt=0
+  call dierckx_surfit(iopt,n,x,y,z,w,xmin,xmax,ymin,ymax,kspl,kspl, &
+              sspl,nxest,nyest,ne,comp_eps,nknt_x,tx_tmp,nknt_y,ty_tmp, &
+              coeff,resid,wrk1,lwrk1,wrk2,lwrk2,iwrk,kwrk,ierr)
+  tx(1:nknt_x)=tx_tmp(1:nknt_x)
+  ty(1:nknt_y)=ty_tmp(1:nknt_y)
+
+  deallocate(wrk1,wrk2,iwrk)
+  deallocate(tx_tmp,ty_tmp)
+
+end subroutine scatterspl