updated makefile, added manpage, fixed few print statements

Makefile

@@ -14,8 +14,12 @@ vpath %.f src
 
 # Fortran compiler name and flags
 FC=gfortran
-FFLAGS=-O3
+LD=gfortran
-# FFLAGS=-O0 -Wall -g -finit-real=nan -ffpe-trap=invalid -fcheck=all -fbounds-check
+FFLAGS=-O3 -fPIC
+LDFLAGS= -static-libgfortran -static-libgcc
+# LDFLAGS= -static-libgfortran -lgfortran -lgcc -lSystem -nodefaultlibs
+LIBS= /usr/local/lib/libquadmath.a
+# FFLAGS=-Wall -g -finit-real=nan -ffpe-trap=invalid -fcheck=all -fbounds-check
 
 DIRECTIVES = -DREVISION="'$(shell svnversion src)'"
 
@@ -23,7 +27,7 @@ all: $(EXE)
 
 # Build executable from object files
 $(EXE): $(MAINOBJ) $(OTHOBJ)
-	$(FC) $(FFLAGS) -o $@ $^
+	$(LD) $(LDFLAGS) -o $@ $^ $(LIBS)
 
 # Dependencies on modules
 main.o: const_and_precisions.o beams.o coreprofiles.o equilibrium.o \

configure

@@ -0,0 +1,31 @@
+#!/bin/sh
+
+printf "Running on "
+case `uname -s` in
+  Linux*) echo "Linux" ;;
+  Darwin*) echo "Mac" ;;
+  [CIGWIN,MINGW]*) echo "Windows" ;;
+  *) echo "unknown OS" ;;
+esac
+
+printf "OS version " && uname -r
+
+arch=`uname -m`
+echo "Processor architecture $arch"
+
+check() {
+  printf "checking for %s... " $1 && command -v $1 || (echo "not found"; exit 1)
+  return $?
+}
+
+for FC in ifort gfortran f77; do
+  check $FC && break
+done
+if [ $? -ne 0 ]; then
+  echo "Fortran compiler not found. Cannot proceed"
+  exit 1
+else
+  echo "Set $FC as Fortran compiler"
+fi
+
+$FC --version

gray.1

@@ -0,0 +1,47 @@
+.TH GRAY 1
+.SH NAME
+gray \- A quasi-optical beam-tracing code
+.SH SYNOPSIS
+.B gray
+[\fB\-hnqvV\fR]
+[\fB\-p\fR \fIparam\-file\fR]
+[\fIoutput\-dir\fR]
+.SH DESCRIPTION
+GRAY is a quasi-optical beam-tracing code for Electron Cyclotron Gaussian beams in tokamak plasmas.
+It computes wave absorption with a fully-relativistic formulation and current drive including momentum conservation.
+.SH OPTIONS
+.TP
+\fB\-h, \-\-help\fR
+Display help and exit
+.TP
+\fB\-n, \-\-no-output\fR
+Suppress output files
+.TP
+\fB\-p \fIfilename\fB, \-\-param\-file\fR=\fIfilename\fR
+Set the parameters file name (default \fIfilename\fR=gray_params.data)
+.TP
+\fB\-q, \-\-quiet\fR
+Quiet mode.
+Suppress all messages on standard output
+.TP
+\fB\-v, \-\-verbose\fR
+Verbose mode.
+Additional information messages on standard output
+.TP
+\fB\-V, \-\-version\fR
+Display version information and exit
+.SH EXIT STATUS
+0 on success, 1 if an error occurred while reading the inputs or writing the outputs, 2 if an error occurred during the simulation.
+.SH AUTHOR
+Daniela Farina, Lorenzo Figini.
+Istituto per la Scienza e Tecnologia dei Plasmi - CNR Milano, Italy
+.PP
+The reference paper for the code is
+.PP
+.nf
+.RS
+D. Farina, Fusion Sci. Technol. 52 (2007) 154\-160
+.RE
+.fi
+.SH SEE ALSO
+The full documentation for the code is available at ...

src/graycore.f90

@@ -28,6 +28,7 @@ contains
     use reflections, only : inside
     use multipass, only : alloc_multipass, dealloc_multipass, initbeam,  &
        initmultipass, turnoffray, plasma_in, plasma_out, wall_out
+    use units, only : ucenr
 
     implicit none
 ! arguments
@@ -218,7 +219,7 @@ contains
         end if
         
         call vectinit(psjki,ppabs,ccci,tau0,alphaabs0,dids0,ccci0,iiv)
-        call print_headers((/' '/),index_rt)
+!        call print_headers((/' '/),index_rt)
         
         if(ip.eq.1) then                                                                   ! 1st pass
           igrad_b = igrad                                                                  ! * input value, igrad_b=0 from 2nd pass
@@ -322,7 +323,7 @@ contains
                 
                 if(jk.eq.1) then
                   write(*,*)
-                  write(*,'("1st pass coupling (central ray, ",a1,"-mode)",f9.4)'), &
+                  write(*,'("1st pass coupling (central ray, ",a1,"-mode)",f9.4)') &
                    mode(iox),cpl(iox)
                   psipv(index_rt) = psipol                                                 !   + polarization angles at plasma boundary for central ray
                   chipv(index_rt) = chipol
@@ -460,9 +461,12 @@ contains
               ccci(jk,i:nstep)  = ccci(jk,i-1)
               psjki(jk,i:nstep) = psjki(jk,i-1)
             else
-              call print_output(i,jk,stv(jk),p0ray(jk)/etau1(jk),xv,psinv,   &
+              call print_output(i,jk,stv(jk),p0ray(jk),xv,psinv,   &
                  btot,bv,ak0,anpl,anpr,anv,anprim,dens,tekev,alpha,tau,dids, &
                  nharm,nhf,iokhawa,index_rt,ddr,ddi,xg,yg,derxg)                           ! p0ray/etau1 [dids normalization] = fraction of p0 coupled to this ray (not including absorption from previous passes)
+!               call print_output(i,jk,stv(jk),p0ray(jk)/etau1(jk),xv,psinv,   &
+!                  btot,bv,ak0,anpl,anpr,anv,anprim,dens,tekev,alpha,tau,dids, &
+!                  nharm,nhf,iokhawa,index_rt,ddr,ddi,xg,yg,derxg)                           ! p0ray/etau1 [dids normalization] = fraction of p0 coupled to this ray (not including absorption from previous passes)
             end if
             
           end do
@@ -483,7 +487,7 @@ contains
           call check_err(ierr,istop)
     ! test whether further trajectory integration is unnecessary
           call vmaxmin(tau1+tau0+lgcpl1,nray,taumn,taumx)                                  ! test on tau + coupling
-          if(taumn > taucr .or. all(iroff(:,index_rt))) istop = 1                          ! (residual power~0) or (no ray active) => stop beam
+!          if(taumn > taucr .or. all(iroff(:,index_rt))) istop = 1                          ! (residual power~0) or (no ray active) => stop beam
           
           if(istop == 1) then                                                              ! stop propagation for current beam
             istop_pass = istop_pass +1                                                     ! * +1 non propagating beam
@@ -528,7 +532,7 @@ contains
         
     !  print final results for pass on screen
         write(*,*)
-        write(*,'("End of propagation for beam ",i5," (pass ",i3,", ",a1," mode)")'), &
+        write(*,'("End of propagation for beam ",i5," (pass ",i3,", ",a1," mode)")') &
                  index_rt,ip,mode(iox)
         write(*,'(a,f9.4)') 'final step (s, ct, Sr)  = ',stv(1)
         write(*,'(a,2e12.5)') 'taumn, taumx  = ', taumn,taumx
@@ -540,6 +544,8 @@ contains
                                           cpl_cbeam1,cpl_cbeam2                            ! central ray coupling for next O/X beams
         end if
         
+        write(ucenr,*) ''
+
         call print_pec(rhop_tab,rhot_tab,jphi_beam,jcd_beam,dpdv_beam,currins_beam, &
             pins_beam,ip)                                                                  !  *print power and current density profiles for current beam
             
@@ -562,7 +568,7 @@ contains
         
     !  print final results for pass on screen
       write(*,*)
-      write(*,'("# End of pass ",i3)'),ip
+      write(*,'("# End of pass ",i3)') ip
       write(*,'(a,f9.4,f9.4)') '# Pabs_tot (MW) [O,X mode] = ',pabs_pass(1),pabs_pass(2)
       write(*,'(a,f9.4,f9.4)') '# I_tot (kA) [O,X mode]    = ', &
                icd_pass(1)*1.0e3_wp_,icd_pass(2)*1.0e3_wp_
@@ -1995,8 +2001,6 @@ bb: do
       write(upec,'(1x,a)') strheader(i)
       write(usumm,'(1x,a)') strheader(i)
     end do
-
-    if(index_rt.gt.0) return
     
     write(uprj0,'(1x,a)')   '#sst j k xt yt zt rt'
     write(uprj0+1,'(1x,a)') '#sst j k xt yt zt rt'
@@ -2007,7 +2011,7 @@ bb: do
     write(uoutr,'(1x,a)')   '#i k sst x y R z psin tau Npl alpha index_rt'
     write(upec,'(1x,a)')    '#rhop rhot Jphi Jcdb dPdV Icdins Pins index_rt'
     write(usumm,'(1x,a)')   '#Icd Pa Jphip dPdVp rhotj rhotjava rhotp rhotpav ' // &
-        'drhotjava drhotpav ratjamx ratjbmx stmx psipol chipol index_rt' // &
+        'drhotjava drhotpav ratjamx ratjbmx stmx psipol chipol index_rt ' // &
         'Jphimx dPdVmx drhotj drhotp P0 cplO cplX'
   end subroutine print_headers
 
@@ -2285,7 +2289,7 @@ bb: do
         rhot=1.0_wp_
       end if
       akim=anprim*ak0*1.0e2_wp_
-      pt=exp(-tau)
+      pt=qj*exp(-tau)
       didsn=dids*1.0e2_wp_/qj
 
       write(ucenr,'(22(1x,e16.8e3),4i5,6(1x,e16.8e3))') stm,rrm,zzm,phideg,       &
@@ -2339,11 +2343,11 @@ bb: do
         stmx,psipol,chipol,p0,cpl1,cpl2
     integer, intent(in) :: index_rt
     
-    write(usumm,'(15(1x,e12.5),i5,4(1x,e12.5))') icd,pabs,jphip,dpdvp,    &
+    write(usumm,'(15(1x,e12.5),i5,7(1x,e12.5))') icd,pabs,jphip,dpdvp,    &
          rhotj,rhotjava,rhotp,rhotpav,drhotjava,drhotpav,ratjamx,ratjbmx, &
          stmx,psipol,chipol,index_rt,jphimx,dpdvmx,drhotj,drhotp,p0,      &
          cpl1,cpl2
-    write(usumm,*) ''
+!    write(usumm,*) ''
   end subroutine print_finals
 
 end module graycore

src/multipass.f90

@@ -147,10 +147,10 @@ contains
     if(all(iwait)) then                                                                    ! no rays active => stop beam
       iboff = .true.
       write(*,*)
-      write(*,'("Beam ",i5," inactive")'),i
+      write(*,'("Beam ",i5," inactive")') i
     else if(.not.all(.not.iwait)) then                                                     ! only some rays active
       write(*,*)
-      write(*,'("WARNING: not all rays in beam ",i5," are active")'),i
+      write(*,'("WARNING: not all rays in beam ",i5," are active")') i
     end if
 
     stv          = zero                                                                    ! starting step