re-added routines to sum profiles (removed by mistake in rev 176)

Makefile.sum

@@ -0,0 +1,75 @@
+# Executable name
+EXE=sumgray
+
+# Objects list
+MAINOBJ=main-sum.o
+OTHOBJ= beamdata.o beams.o conical.o const_and_precisions.o coreprofiles.o \
+    dierckx.o dispersion.o eccd.o eierf.o errcodes.o sumcore.o \
+    gray_params.o equilibrium.o limiter.o magsurf_data.o math.o minpack.o numint.o \
+    pec.o polarization.o quadpack.o reflections.o simplespline.o units.o utils.o
+    
+# Alternative search paths
+vpath %.f90 src
+vpath %.f src
+
+# Fortran compiler name and flags
+FC=gfortran
+FFLAGS=-O3
+#FFLAGS=-Wall -g -finit-real=nan -ffpe-trap=invalid -fcheck=all -fbounds-check
+
+DIRECTIVES = -DREVISION="'$(shell svnversion src)'"
+
+all: $(EXE)
+
+# Build executable from object files
+$(EXE): $(MAINOBJ) $(OTHOBJ)
+	$(FC) $(FFLAGS) -o $@ $^
+
+# Dependencies on modules
+main-sum.o: const_and_precisions.o beams.o coreprofiles.o equilibrium.o \
+        sumcore.o gray_params.o reflections.o
+sumcore.o: const_and_precisions.o beamdata.o beams.o coreprofiles.o \
+            dispersion.o eccd.o equilibrium.o errcodes.o gray_params.o \
+            pec.o polarization.o limiter.o units.o utils.o
+beams.o: const_and_precisions.o dierckx.o reflections.o simplespline.o utils.o
+beamdata.o: const_and_precisions.o gray_params.o 
+conical.o: const_and_precisions.o
+coreprofiles.o: const_and_precisions.o dierckx.o gray_params.o simplespline.o \
+                utils.o
+dierckx.o: const_and_precisions.o
+dispersion.o: const_and_precisions.o eierf.o errcodes.o math.o quadpack.o
+eccd.o: const_and_precisions.o conical.o dierckx.o errcodes.o magsurf_data.o \
+        numint.o
+eierf.o: const_and_precisions.o
+errcodes.o: const_and_precisions.o
+gray_params.o: const_and_precisions.o utils.o
+equilibrium.o: const_and_precisions.o dierckx.o limiter.o minpack.o \
+               reflections.o simplespline.o utils.o gray_params.o
+magsurf_data.o: const_and_precisions.o gray_params.o equilibrium.o dierckx.o \
+                reflections.o simplespline.o units.o utils.o
+math.o: const_and_precisions.o
+minpack.o: const_and_precisions.o
+numint.o: const_and_precisions.o
+pec.o: const_and_precisions.o beamdata.o equilibrium.o gray_params.o \
+       magsurf_data.o utils.o
+polarization.o: const_and_precisions.o
+quadpack.o: const_and_precisions.o
+reflections.o: const_and_precisions.o limiter.o utils.o
+simplespline.o: const_and_precisions.o
+utils.o: const_and_precisions.o
+
+# General object compilation command
+%.o: %.f90
+	$(FC) -cpp $(DIRECTIVES) $(FFLAGS) -c $<
+
+.PHONY: clean install
+# 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/main-sum.f90

@@ -0,0 +1,195 @@
+program main_sum
+  use const_and_precisions, only : wp_,one
+  use sumcore, only : sum_profiles
+  use gray_params, only : read_params, antctrl_type,eqparam_type,  &
+     prfparam_type,outparam_type,rtrparam_type,hcdparam_type
+  use beams, only : read_beam0, read_beam1, read_beam2
+  use equilibrium, only : read_equil_an,read_eqdsk,change_cocos,eq_scal, &
+     set_rhospl,setqphi_num,frhopolv
+  use coreprofiles, only : read_profiles_an,read_profiles,tene_scal
+  use reflections, only : range2rect
+  implicit none
+  type(antctrl_type) :: antp
+  type(eqparam_type) :: eqp
+  type(prfparam_type) :: prfp
+  type(outparam_type) :: outp
+  type(rtrparam_type) :: rtrp
+  type(hcdparam_type) :: hcdp
+
+  real(wp_), dimension(:), allocatable :: psrad, terad, derad, zfc
+  real(wp_), dimension(:), allocatable :: rv, zv, psinr, fpol, qpsi
+  real(wp_), dimension(:), allocatable :: rbnd, zbnd, rlim, zlim
+  real(wp_), dimension(:,:), allocatable :: psin
+  real(wp_) :: psia, rvac, rax, zax
+  integer :: iox0
+  real(wp_) :: p0mw, fghz, psipol0, chipol0
+  real(wp_) :: alpha0, beta0, x0, y0, z0, w1, w2, ri1, ri2, phiw, phir
+
+  real(wp_) :: pec,icd
+
+  integer :: ierr
+  real(wp_), dimension(:), allocatable :: xrad, rhot
+  real(wp_), dimension(:), allocatable :: jphi,jcd,dpdv,currins,pins,rtin,rpin
+  real(wp_) :: rwallm, rmxm, r0m, z0m, dzmx
+  
+  integer :: i,j,k,n,ngam,irt
+  character(len=255) :: fn,dumstr
+  real(wp_), dimension(5) :: f48v
+  real(wp_) :: gam,alp,bet,        jphip,dpdvp, &
+        rhotj,rhotjava,rhotp,rhotpav,drhotjava,drhotpav,ratjamx,ratjbmx
+
+
+! ======= read parameters BEGIN =======
+  call read_params('gray_params.data',rtrp,hcdp,antp,eqp,rwallm,prfp,outp)
+! ======= read parameters  END  =======
+
+! ======= read input data BEGIN =======
+!------------ equilibrium ------------
+  if(eqp%iequil<2) then
+    call read_equil_an(eqp%filenm, rv, zv, fpol, qpsi)
+! psia sign set to give the correct sign to Iphi (COCOS=3: psia<0 for Iphi>0)
+    psia = sign(one,qpsi(2)*fpol(1))
+  else
+    call read_eqdsk(eqp%filenm, rv,zv,psin, psia, psinr,fpol,qpsi, rvac,       &
+       rax,zax, rbnd,zbnd, rlim,zlim, eqp%ipsinorm,eqp%idesc,eqp%ifreefmt)
+    call change_cocos(psia, fpol, qpsi, eqp%icocos, 3)
+  end if
+! re-scale B/I and/or force signs. If sgn=0 on input, set to fpol/-psia signs on output
+  call eq_scal(psia, fpol, eqp%sgni, eqp%sgnb, eqp%factb)
+! ??? analytical only? change for numerical!
+!  qpsi(1) = sign(qpsi(1),qpsi(1)*qpsi(2)*psia*fpol(1))
+!  qpsi(2) = sign(qpsi(2),psia*fpol(1))
+!------------- profiles -------------
+  if(prfp%iprof==0) then
+    call read_profiles_an(prfp%filenm, terad, derad, zfc)
+  else
+    call read_profiles(prfp%filenm, xrad, terad, derad, zfc)
+    allocate(psrad(size(xrad)))
+    if(prfp%irho==0) then ! xrad==rhot
+      allocate(rhot(size(psinr)))
+      call setqphi_num(psinr,qpsi,psia,rhot)
+      call set_rhospl(sqrt(psinr),rhot)
+      deallocate(rhot)
+      psrad=frhopolv(xrad)**2
+    else if(prfp%irho == 1) then ! xrad==rhop
+      psrad=xrad**2
+    else
+      psrad=xrad   
+    end if
+    deallocate(xrad)
+  end if
+! re-scale input data
+  call tene_scal(terad,derad,prfp%factte,prfp%factne,eqp%factb,prfp%iscal,     &
+                 prfp%iprof)
+!------------- antenna --------------
+! interpolate beam table if antctrl%ibeam>0
+  select case (antp%ibeam)
+  case (2)
+! to be completed: now 1st beamd always selected, iox read from table
+    call read_beam2(antp%filenm,1,antp%alpha,antp%beta,fghz,antp%iox,x0,y0,z0, &
+                        w1,w2,ri1,ri2,phiw,phir)
+  case (1)
+    call read_beam1(antp%filenm,antp%alpha,antp%beta,fghz,x0,y0,z0,            &
+                        w1,w2,ri1,ri2,phiw,phir)
+  case default
+    call read_beam0(antp%filenm,fghz,x0,y0,z0,w1,w2,ri1,ri2,phiw,phir)
+  end select
+  alpha0=antp%alpha
+  beta0=antp%beta
+  p0mw=antp%power
+  psipol0=antp%psi
+  chipol0=antp%chi
+  iox0=antp%iox
+!--------------- wall ---------------
+! set simple limiter if not read from EQDSK
+! need to clean up...
+  r0m=sqrt(x0**2+y0**2)*0.01_wp_
+  dzmx=rtrp%dst*rtrp%nstep*0.01_wp_
+  z0m=z0*0.01_wp_
+  if (.not.allocated(rlim).or.rtrp%ipass<0) then
+    if (allocated(rlim)) deallocate(rlim)
+    if (allocated(zlim)) deallocate(zlim)
+    allocate(rlim(5))
+    allocate(zlim(5))
+    if (rtrp%ipass<0) rtrp%ipass = -rtrp%ipass
+    if(eqp%iequil<2) then
+      rmxm=(rv(1)+rv(2))*0.01_wp_
+    else
+      rmxm=rv(size(rv))
+    end if
+    call range2rect(rwallm,max(r0m,rmxm),z0m-dzmx,z0m+dzmx,rlim,zlim)
+  end if
+! ======= read input data  END  =======
+
+! ========================= MAIN SUBROUTINE CALL =========================
+  allocate(jphi(outp%nrho),jcd(outp%nrho),dpdv(outp%nrho), &
+           currins(outp%nrho),pins(outp%nrho),rtin(outp%nrho),rpin(outp%nrho))
+
+  open(100,file='filelist.txt',action='read',status='old')
+  read(100,*) n,ngam
+  do i=1,n
+    read(100,*) fn
+    open(100+i,file=fn,action='read',status='old')
+    do j=1,22
+      read(100+i,*) dumstr
+    end do
+  end do
+  close(100)
+
+  open(100+n+1,file='f48sum.txt',action='write',status='unknown')
+  open(100+n+2,file='f7sum.txt',action='write',status='unknown')
+
+  do k=1,ngam
+    jphi=0.0_wp_
+    jcd=0.0_wp_
+    dpdv=0.0_wp_
+    currins=0.0_wp_
+    pins=0.0_wp_
+    do j=1,outp%nrho
+      do i=1,n
+        read(100+i,*) gam,alp,bet,rpin(j),rtin(j),f48v(1:5),irt
+        jphi(j)=jphi(j)+f48v(1)
+        jcd(j)=jcd(j)+f48v(2)
+        dpdv(j)=dpdv(j)+f48v(3)
+        currins(j)=currins(j)+f48v(4)
+        pins(j)=pins(j)+f48v(5)
+      end do
+      write(100+n+1,'(10(1x,e16.8e3),i5)') gam,alp,bet,rpin(j),rtin(j), &
+          jphi(j),jcd(j),dpdv(j),currins(j),pins(j),irt
+    end do
+    pec=pins(outp%nrho)
+    icd=currins(outp%nrho)
+    write(100+n+1,*)
+    call sum_profiles(rv,zv,psin,psia,psinr,fpol,qpsi,rvac,rax,zax,rbnd,zbnd,eqp,&
+        psrad,terad,derad,zfc,prfp,  rlim,zlim,                            &
+        p0mw,fghz,alpha0,beta0,(/x0,y0,z0/),w1,w2,ri1,ri2,phiw,phir,iox0,  &
+        psipol0,chipol0,  jphi,jcd,dpdv,currins,pins,pec,icd,              &
+        jphip,dpdvp, &
+        rhotj,rhotjava,rhotp,rhotpav,drhotjava,drhotpav,ratjamx,ratjbmx, &
+        outp,rtrp,hcdp,ierr)
+    write(100+n+2,'(15(1x,e12.5),i5,4(1x,e12.5))') gam,alp,bet,icd,pec, &
+         jphip,dpdvp, &
+         rhotj,rhotjava,rhotp,rhotpav,drhotjava,drhotpav,ratjamx,ratjbmx
+  end do
+  do i=1,n+2
+    close(100+i)
+  end do
+! ========================================================================
+
+! ======= control prints BEGIN =======
+  if(ierr/=0) print*,'    IERR = ', ierr
+  print*,' '
+  print*,'Pabs (MW) = ', pec
+  print*,'Icd (kA) = ', icd*1.0e3_wp_
+! ======= control prints  END  =======
+
+! ======= free memory BEGIN =======
+  if(allocated(psrad)) deallocate(psrad)
+  if(allocated(terad)) deallocate(terad, derad, zfc)
+  if(allocated(rv))    deallocate(rv, zv, fpol, qpsi)
+  if(allocated(psin))  deallocate(psin, psinr)
+  if(allocated(rbnd))  deallocate(rbnd,zbnd)
+  if(allocated(rlim))  deallocate(rlim,zlim)
+  if(allocated(dpdv))  deallocate(jphi,jcd,dpdv,currins,pins,rtin,rpin)
+! ======= free memory  END  ======
+end program main_sum