new branch for JINTRAC integration

Makefile

@@ -1,8 +1,9 @@
 # Executable name
 EXE=gray
+LIBFILE=lib$(EXE).a
 
 # Objects list
-OBJ=gray.o grayl.o reflections.o green_func_p.o \
+OBJ=gray_main.o gray-externals.o grayl.o reflections.o green_func_p.o \
     const_and_precisions.o itm_constants.o itm_types.o
 
 # Alternative search paths
@@ -11,18 +12,25 @@ vpath %.f src
 
 # Fortran compiler name and flags
 FC=gfortran
-FFLAGS=-O3
+FFLAGS=-O0 -g -fbounds-check -Wall
+#FFLAGS=-O3
 
-DIRECTIVES = -DREVISION="'$(shell svnversion src)'"
+DIRECTIVES = -DREVISION="'$(shell svnversion)'"
 
-all: $(EXE)
+all: library
 
 # Build executable from object files
-$(EXE): $(OBJ)
+$(EXE): main.o $(OBJ)
 	$(FC) $(FFLAGS) -o $@ $^
 
+$(LIBFILE): $(OBJ)
+	rm -f $@
+	ar -rv $@ $^
+
 # Dependencies on modules
-gray.o: green_func_p.o reflections.o
+main.o: gray_main.o
+gray_main.o: gray-externals.o itm_types.o
+gray-externals.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
@@ -31,13 +39,16 @@ itm_constants.o: itm_types.o
 %.o: %.f90
 	$(FC) $(FFLAGS) -c $<
 
-gray.o:gray.f green_func_p.o
+gray-externals.o:gray-externals.f
 	$(FC) -cpp $(DIRECTIVES) $(FFLAGS) -c $< 
 
 grayl.o:grayl.f
-	$(FC) $(FFLAGS) -c $^ 
+	$(FC) $(FFLAGS) -c $<
+
+.PHONY: library clean install
+
+library: $(LIBFILE)
 
-.PHONY: clean install
 # Remove output files
 clean:
 	rm -rf *.o *.mod $(EXE)

src/gray_main.f90

@@ -0,0 +1,95 @@
+subroutine gray_main(ijetto, mr, mz, r, z, psin, psiax, psibnd, &
+    rax, zax, nbnd, rbnd, zbnd, nrho, psijet, f, te, dne, zeff, qsf, &
+    powin, dpdv, jcd, pec, icd, ierr)
+
+  use itm_types, only : r8
+  implicit none
+
+  integer, intent(in) :: ijetto, mr, mz, nrho, nbnd
+  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
+  real(r8), intent(out), dimension(nrho) :: dpdv, jcd
+  real(r8), intent(out) :: pec, icd
+  integer, intent(out) :: ierr
+
+  integer :: istop, iercom, igrad, iopmin, iowmin, index_rt, ipass
+  real(r8) :: sox, p0mw, powrfl, taumn, taumx, pabstot, currtot
+  real(r8) :: p0mw1, powtr, pabstott, currtott
+
+  common/istop/istop
+  common/ierr/iercom
+  common/igrad/igrad
+  common/iovmin/iopmin,iowmin
+  common/mode/sox
+  common/p0/p0mw
+  common/powrfl/powrfl
+  common/index_rt/index_rt
+  common/taumnx/taumn,taumx,pabstot,currtot
+  common/ipass/ipass
+      
+! 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)
+  call vectinit
+  if(igrad.eq.0) call ic_rt
+  if(igrad.gt.0) call ic_gb
+  if(iercom.gt.0) then
+    ierr=iercom
+    write(*,*) '    IERR = ', ierr
+    return
+  end if
+
+! beam/ray propagation
+  call gray_integration
+
+! postprocessing
+  call after_gray_integration
+
+  pabstott=pabstot
+  currtott=currtot
+  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
+    pabstott=pabstott+pabstot
+    currtott=currtott+currtot
+  
+    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
+    pabstott=pabstott+pabstot
+    currtott=currtott+currtot
+  end if
+  write(*,*)
+  write(*,*) 'Pabs (MW), Icd (kA) = ', pabstott,currtott*1.e3_r8
+  pec=pabstott*1.e6_r8
+  icd=currtott*1.e6_r8
+! NOT YET COMPLETE! 
+  dpdv=0._r8
+  jcd=0._r8
+! =================
+  ierr=iercom
+end subroutine gray_main