fixed incompatibilities with allocatable arrays (pgf90 compiler, JETTO); fixed input/output scaling (gray/jetto interface); fixed error in read_beams2 for small beam tables

Makefile.jetto

@@ -16,7 +16,7 @@ include ../include.mk
 # Alternative search paths
 vpath %.f90 src
 
-DIRECTIVES = -DEXTBES -DREVISION="rev.155 JETTO" #'$(shell svnversion src)'"
+DIRECTIVES = -DEXTBES -DREVISION="rev.164 JETTO" #'$(shell svnversion src)'"
 
 # library name
 # ------------

src/beamdata.f90

@@ -13,13 +13,13 @@ contains
     use const_and_precisions, only : zero,half,two
     implicit none
     type(rtrparam_type), intent(in) :: rtrparam
-    real(wp_),  dimension(:,:), intent(out), allocatable :: ywork,ypwork, &
+    real(wp_),  dimension(:,:), intent(out), pointer :: ywork,ypwork, &
        gri,psjki,ppabs,ccci
-    real(wp_), dimension(:,:,:), intent(out), allocatable :: xc,du1,ggri
+    real(wp_), dimension(:,:,:), intent(out), pointer :: xc,du1,ggri
-    real(wp_),    dimension(:), intent(out), allocatable :: tau0,alphaabs0, &
+    real(wp_),    dimension(:), intent(out), pointer :: tau0,alphaabs0, &
        dids0,ccci0,p0jk
-    complex(wp_), dimension(:), intent(out), allocatable :: ext, eyt
+    complex(wp_), dimension(:), intent(out), pointer :: ext, eyt
-    integer,      dimension(:), intent(out), allocatable :: iiv
+    integer,      dimension(:), intent(out), pointer :: iiv
     
     integer :: jray1
 
@@ -164,13 +164,13 @@ contains
   subroutine alloc_beam(ywork,ypwork,xc,du1,gri,ggri,psjki,ppabs,ccci, &
                         tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv)
     implicit none
-    real(wp_),  dimension(:,:), intent(out), allocatable :: ywork,ypwork, &
+    real(wp_),  dimension(:,:), intent(out), pointer :: ywork,ypwork, &
        gri,psjki,ppabs,ccci
-    real(wp_), dimension(:,:,:), intent(out), allocatable :: xc,du1,ggri
+    real(wp_), dimension(:,:,:), intent(out), pointer :: xc,du1,ggri
-    real(wp_),    dimension(:), intent(out), allocatable :: tau0,alphaabs0, &
+    real(wp_),    dimension(:), intent(out), pointer :: tau0,alphaabs0, &
        dids0,ccci0,p0jk
-    complex(wp_), dimension(:), intent(out), allocatable :: ext, eyt
+    complex(wp_), dimension(:), intent(out), pointer :: ext, eyt
-    integer,      dimension(:), intent(out), allocatable :: iiv
+    integer,      dimension(:), intent(out), pointer :: iiv
     
     call dealloc_beam(ywork,ypwork,xc,du1,gri,ggri,psjki,ppabs,ccci, &
                       tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv)
@@ -187,31 +187,31 @@ contains
   subroutine dealloc_beam(ywork,ypwork,xc,du1,gri,ggri,psjki,ppabs,ccci, &
                           tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv)
     implicit none
-    real(wp_),  dimension(:,:), intent(out), allocatable :: ywork,ypwork, &
+    real(wp_),  dimension(:,:), intent(out), pointer :: ywork,ypwork, &
        gri,psjki,ppabs,ccci
-    real(wp_), dimension(:,:,:), intent(out), allocatable :: xc,du1,ggri
+    real(wp_), dimension(:,:,:), intent(out), pointer :: xc,du1,ggri
-    real(wp_),    dimension(:), intent(out), allocatable :: tau0,alphaabs0, &
+    real(wp_),    dimension(:), intent(out), pointer :: tau0,alphaabs0, &
        dids0,ccci0,p0jk
-    complex(wp_), dimension(:), intent(out), allocatable :: ext, eyt
+    complex(wp_), dimension(:), intent(out), pointer :: ext, eyt
-    integer,      dimension(:), intent(out), allocatable :: iiv
+    integer,      dimension(:), intent(out), pointer :: iiv
 
-    if (allocated(ywork))     deallocate(ywork)
+    if (associated(ywork))     deallocate(ywork)
-    if (allocated(ypwork))    deallocate(ypwork)
+    if (associated(ypwork))    deallocate(ypwork)
-    if (allocated(xc))        deallocate(xc)
+    if (associated(xc))        deallocate(xc)
-    if (allocated(du1))       deallocate(du1)
+    if (associated(du1))       deallocate(du1)
-    if (allocated(gri))       deallocate(gri)
+    if (associated(gri))       deallocate(gri)
-    if (allocated(ggri))      deallocate(ggri)
+    if (associated(ggri))      deallocate(ggri)
-    if (allocated(psjki))     deallocate(psjki)
+    if (associated(psjki))     deallocate(psjki)
-    if (allocated(ppabs))     deallocate(ppabs)
+    if (associated(ppabs))     deallocate(ppabs)
-    if (allocated(ccci))      deallocate(ccci)
+    if (associated(ccci))      deallocate(ccci)
-    if (allocated(tau0))      deallocate(tau0)
+    if (associated(tau0))      deallocate(tau0)
-    if (allocated(alphaabs0)) deallocate(alphaabs0)
+    if (associated(alphaabs0)) deallocate(alphaabs0)
-    if (allocated(dids0))     deallocate(dids0)
+    if (associated(dids0))     deallocate(dids0)
-    if (allocated(ccci0))     deallocate(ccci0)
+    if (associated(ccci0))     deallocate(ccci0)
-    if (allocated(p0jk))      deallocate(p0jk)
+    if (associated(p0jk))      deallocate(p0jk)
-    if (allocated(ext))       deallocate(ext)
+    if (associated(ext))       deallocate(ext)
-    if (allocated(eyt))       deallocate(eyt)
+    if (associated(eyt))       deallocate(eyt)
-    if (allocated(iiv))       deallocate(iiv)
+    if (associated(iiv))       deallocate(iiv)
   end subroutine dealloc_beam
 
 end module beamdata

src/beams.f90

@@ -269,6 +269,8 @@ contains
         rcieta=rci2v(1)
         phiw=phi1v(1)
         phir=phi2v(1)
+        deallocate(alphastv,betastv,waist1v,waist2v,rci1v,rci2v,phi1v, &
+        phi2v,x00v,y00v,z00v,xcoord,ycoord)
         return
       end if
 !#######################################################################################
@@ -368,9 +370,9 @@ contains
         maxx = maxval(xcoord)
         miny = minval(ycoord)
         maxy = maxval(ycoord)
-        nxest = kx + 1 + int(sqrt(nisteer/2.))
+        nxest = 2*(kx + 1)
-        nyest = ky + 1 + int(sqrt(nisteer/2.))
+        nyest = 2*(ky + 1)
-        nmax = max(nxest,nyest)
+        nmax = max(nxest,nyest)+max(kx,ky)
         eps = 10.**(-8)
         lwrk = (nmax-2)**2*(7*nmax-2)+18*nmax+8*nisteer-19
         lwrk2 = (nmax-2)**2*(4*nmax-1)+4*nmax-2
@@ -636,7 +638,8 @@ contains
           END SELECT
 !    c----------------------------------------------------------------------------------
 !     
-          deallocate(xpolygA, ypolygA, xpolygC, ypolygC, xpolygB, ypolygB, xpolygD, ypolygD)
+          deallocate(xpolygA, ypolygA, xpolygC, ypolygC, xpolygB, ypolygB, &
+                     xpolygD, ypolygD, xoutA, youtA, xoutB, youtB, xoutC, youtC)
         end if
 !  c====================================================================================  
 !

src/eccd.f90

@@ -103,7 +103,7 @@ contains
     implicit none
     real(wp_), intent(in) :: zeff
     real(wp_), intent(out) :: cst2
-    real(wp_), dimension(:), allocatable, intent(out) :: eccdpar
+    real(wp_), dimension(:), pointer, intent(out) :: eccdpar
 
     cst2=0.0_wp_
     allocate(eccdpar(1))
@@ -122,7 +122,7 @@ contains
     implicit none
     real(wp_), intent(in) :: zeff,rbn,rbx
     real(wp_), intent(out) :: cst2
-    real(wp_), dimension(:), allocatable, intent(out) :: eccdpar
+    real(wp_), dimension(:), pointer, intent(out) :: eccdpar
     real(wp_) :: alams,pa,fp0s
 
     cst2=1.0_wp_-rbx
@@ -145,7 +145,7 @@ contains
     integer, parameter :: ksp=3
     real(wp_), intent(in) :: zeff,rbx,fc,amu,rhop
     real(wp_), intent(out) :: cst2
-    real(wp_), dimension(:), allocatable, intent(out) :: eccdpar
+    real(wp_), dimension(:), pointer, intent(out) :: eccdpar
     real(wp_), dimension(nlmt) :: chlm
     integer :: nlm,ierr,npar
 
@@ -189,13 +189,15 @@ contains
       rdu,rdu2t,duu,uu1,uu2,xx1,xx2,resj,resp,eji,epp,anum,denom,        &
       cstrdut,anucc
     real(wp_), dimension(lw) :: w
-    real(wp_), dimension(nfpp+size(eccdpar)) :: apar
+    real(wp_), dimension(:), pointer :: apar=>null()
     real(wp_), dimension(0:1) :: uleft,uright
 ! common/external functions/variables
     real(wp_), external :: fcur
 !
 !     effjpl = <J_parallel>/<p_d> /(B_min/<B>)  [A m /W]
 !
+    allocate(apar(nfpp+size(eccdpar)))
+    
     apar(1) = yg
     apar(2) = anpl
     apar(3) = amu
@@ -310,6 +312,8 @@ contains
       effjcd=-ceff*anum/(anucc*denom)
     end if
 
+    deallocate(apar)
+
   end subroutine eccdeff
 
   function fpp(upl,extrapar,npar)

src/gray_jetto1beam.f90

@@ -1,6 +1,6 @@
 subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
     nbnd, rbnd, zbnd, nrho, psrad, fpol, te, dne, zeff, qpsi, ibeam,   &
-    powin, alphain, betain, dpdv, jcd, pabs, icd, ierr)
+    p0mw, alphain, betain, dpdv, jcd, pabs, icd, ierr)
   use const_and_precisions, only : wp_
   use units, only : ucenr,usumm,uprj0,uwbm,udisp,ubres,ucnt,uoutr,ueq,uprfin,  &
       uflx,upec,uprm,ubeam
@@ -9,13 +9,14 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
   use beams, only : read_beam2
   use graycore, only : gray_main
   use reflections, only : range2rect
+  use coreprofiles, only : tene_scal
   implicit none
 ! arguments
   integer, intent(in) :: ijetto, mr, mz, nbnd, nrho, ibeam
   real(wp_), dimension(mr), intent(in) :: r
   real(wp_), dimension(mz), intent(in) :: z
   real(wp_), dimension(mr,mz), intent(in) :: psin
-  real(wp_), intent(in) :: psia, rax, zax, powin, alphain, betain
+  real(wp_), intent(in) :: psia, rax, zax, p0mw, alphain, betain
   real(wp_), dimension(nbnd), intent(in) :: rbnd, zbnd
   real(wp_), dimension(nrho), intent(in) :: psrad, fpol, te, dne, zeff, qpsi
   real(wp_), dimension(nrho), intent(out) :: dpdv, jcd
@@ -24,7 +25,7 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
 ! local variables
   real(wp_), dimension(nrho) :: psinr
   integer :: iox0
-  real(wp_) :: r0m,rvac,alpha0,beta0,psipol0,chipol0,p0mw,rwallm
+  real(wp_) :: r0m,rvac,alpha0,beta0,psipol0,chipol0,rwallm
   real(wp_) :: fghz,x0,y0,z0,w1,w2,ri1,ri2,phiw,phir
   real(wp_), dimension(5) :: rlim,zlim
   logical, save :: firstcall=.true.
@@ -45,8 +46,12 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
     prfp%filenm='JETTO'
     outp%ipec=1
     firstcall=.false.
+    outp%nrho=nrho
   end if
   
+!  call tene_scal(te,dne,prfp%factte,prfp%factne,eqp%factb,&
+!                 prfp%iscal,prfp%iprof)
+  
   rvac=rax
   psinr=psrad
 
@@ -57,7 +62,6 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
   psipol0=antp%psi
   chipol0=antp%chi
   iox0=antp%iox
-  p0mw=powin*1.e-6_wp_
 
 ! set simple limiter
   r0m=sqrt(x0**2+y0**2)*0.01_wp_
@@ -84,4 +88,4 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
   close(uflx)
   close(upec)
 
-end subroutine gray_jetto1beam
+end subroutine gray_jetto1beam

src/graycore.f90

@@ -9,21 +9,21 @@ contains
               p0,fghz,alpha0,beta0,xv0,w1,w2,ri1,ri2,phiw,phir,iox0,           &
               psipol0,chipol0,  dpdv,jcd,pabs,icd,  outp,rtrp,hcdp,ierr, rhout)
     use const_and_precisions, only : zero, one, degree
-    use coreprofiles, only : set_prfan, set_prfspl, temp, fzeff
+    use coreprofiles, only : set_prfan, set_prfspl, temp, fzeff, unset_prfspl
     use dispersion, only : expinit
     use gray_params, only : eqparam_type, prfparam_type, outparam_type, &
        rtrparam_type, hcdparam_type, antctrl_type, set_codepar, print_params,  &
-       iequil, iprof, ieccd, iwarm, ipec, istpr0, igrad, headw, headl
+       iequil, iprof, iwarm, ipec, istpr0, igrad, headw, headl
     use beams, only : read_beam0, read_beam1, launchangles2n, xgygcoeff
     use beamdata, only : pweight, rayi2jk
     use equilibrium, only : set_equian, set_eqspl, setqphi_num, set_rhospl,    &
-       zbinf, zbsup
+       zbinf, zbsup, unset_eqspl, unset_rhospl, unset_q
     use errcodes, only : check_err, print_errn, print_errhcd
-    use magsurf_data, only : flux_average
+    use magsurf_data, only : flux_average, dealloc_surfvec
     use beamdata, only : init_btr, dealloc_beam, nray, nstep, dst
     use pec, only : pec_init, spec, postproc_profiles, dealloc_pec, &
        rhop_tab, rhot_tab
-    use limiter, only : set_lim
+    use limiter, only : set_lim, unset_lim
     use utils, only : vmaxmin
     implicit none
 ! arguments
@@ -51,8 +51,6 @@ contains
 ! local variables
     real(wp_), parameter :: taucr = 12._wp_
 
-    real(wp_), dimension(:), allocatable :: rhotn
-      
     real(wp_) :: sox,ak0,bres,xgcn,xg,yg,zzm,alpha,didp,anpl,anpr,anprim,anprre
     real(wp_) :: chipol,psipol,btot,psinv,dens,tekev,dersdst,derdnm,st,st0
     real(wp_) :: tau,pow,dids,ddr,ddi,taumn,taumx
@@ -60,16 +58,17 @@ contains
     real(wp_) :: rhotp,drhotp,rhotj,drhotj,dpdvmx,jphimx,ratjamx,ratjbmx
 
     real(wp_), dimension(3) :: xv,anv0,anv,bv
-    real(wp_), dimension(:,:), allocatable :: yw,ypw,gri
+    real(wp_), dimension(:,:), pointer :: yw=>null(),ypw=>null(),gri=>null()
-    real(wp_), dimension(:,:,:), allocatable :: xc,du1,ggri
+    real(wp_), dimension(:,:,:), pointer :: xc=>null(),du1=>null(),ggri=>null()
     integer :: i,jk,iox,nharm,nhf,nnd,iokhawa,istop,ierrn,ierrhcd,index_rt=1
     logical :: ins_pl, somein, allout
 
-    real(wp_),  dimension(:,:), allocatable :: psjki,ppabs,ccci
+    real(wp_),  dimension(:,:), pointer :: psjki=>null(),ppabs=>null(),ccci=>null()
-    real(wp_), dimension(:), allocatable :: tau0,alphaabs0,dids0,ccci0
+    real(wp_), dimension(:), pointer :: tau0=>null(),alphaabs0=>null(),dids0=>null(), &
-    real(wp_),    dimension(:), allocatable :: p0jk
+                                        ccci0=>null()
-    complex(wp_), dimension(:), allocatable :: ext, eyt
+    real(wp_),    dimension(:), pointer :: p0jk=>null()
-    integer,      dimension(:), allocatable :: iiv
+    complex(wp_), dimension(:), pointer :: ext=>null(), eyt=>null()
+    integer,      dimension(:), pointer :: iiv=>null()
 
     real(wp_),    dimension(:), allocatable :: jphi,pins,currins
 
@@ -274,15 +273,15 @@ contains
 ! ======= post-proc  END  ======
 
 ! ======= free memory BEGIN ======
-!    call unset_eqspl
+    call unset_eqspl
-!    call unset_q
+    call unset_q
-!    call unset_rhospl
+    call unset_rhospl
-!    call unset_prfspl
+    call unset_prfspl
-!    call unset_lim
+    call unset_lim
-!    call dealloc_surfvec
+    call dealloc_surfvec
-!    call dealloc_beam(yw,ypw,xc,du1,gri,ggri,psjki,ppabs,ccci, &
+    call dealloc_beam(yw,ypw,xc,du1,gri,ggri,psjki,ppabs,ccci, &
-!                      tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv)
+                      tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv)
-!    call dealloc_pec
+    call dealloc_pec
     deallocate(jphi,pins,currins)
 ! ======= free memory  END  ======
   end subroutine gray_main
@@ -1230,7 +1229,7 @@ contains
     integer :: lrm,ithn,ierrcd
     real(wp_) :: amu,ratiovgr,rbn,rbx
     real(wp_) :: zeff,cst2,bmxi,bmni,fci
-    real(wp_), dimension(:), allocatable :: eccdpar
+    real(wp_), dimension(:), pointer :: eccdpar=>null()
     real(wp_) :: effjcd,effjcdav,akim,btot
     complex(wp_) :: ex,ey,ez
 
@@ -1289,7 +1288,7 @@ contains
                          ithn,cst2,fjncl,eccdpar,effjcd,iokhawa,ierrcd)
           end select
           ierr=ierr+ierrcd
-          !deallocate(eccdpar)
+          if(associated(eccdpar)) deallocate(eccdpar)
           effjcdav=rbavi*effjcd
           didp=sgnbphi*effjcdav/(2.0_wp_*pi*rrii)
         end if
@@ -1722,7 +1721,6 @@ bb: do
 
   subroutine print_bres(bres)
     use const_and_precisions, only : wp_
-    use gray_params, only : iequil
     use equilibrium, only : rmnm, rmxm, zmnm, zmxm, bfield, nq
     use units, only : ubres
     implicit none
@@ -1841,7 +1839,7 @@ bb: do
 ! arguments
     real(wp_), dimension(:), intent(in) :: qval
 ! local variables
-    integer :: ncnt,i1,i
+    integer :: i1,i
     real(wp_) :: rup,zup,rlw,zlw,rhot,psival
     real(wp_), dimension(npoints) :: rcn,zcn
     real(wp_), dimension(nq) :: qpsi

src/magsurf_data.f90

@@ -115,7 +115,7 @@ contains
              lwrk=4*(nnintp+nlam)+11*(njest+nlest)+njest*nnintp+nlest+54, &
              kwrk=nnintp+nlam+njest+nlest+3,lw01=nnintp*4+nlam*3+nnintp*nlam
 ! local variables
-    integer :: ier,ierr,l,jp,ipr,inc,inc1,iopt,njp,nlm,ninpr
+    integer :: ier,ierr,l,jp,inc,inc1,iopt,njp,nlm,ninpr
     integer, dimension(kwrk) :: iwrk
     real(wp_) :: lam,dlam,anorm,b2av,dvdpsi,dadpsi,ratio_cdator,     &
                  ratio_cdbtor,ratio_pltor,fc,height,r2iav,currp,     &
@@ -478,7 +478,7 @@ contains
     real(wp_), dimension(:), intent(out) :: rcn,zcn
     real(wp_), intent(inout) :: rup,zup,rlw,zlw
 ! local variables
-    integer :: npoints,np,info,ic,ier,ii,iopt,m
+    integer :: npoints,np,info,ic,ier,iopt,m
     real(wp_) :: ra,rb,za,zb,rn,th,zc,val
     real(wp_), dimension(mest) :: zeroc
     real(wp_), dimension(nsr) :: czc

src/reflections.f90

@@ -203,10 +203,7 @@ end function interssegm
   subroutine range2rect(xmin,xmax,ymin,ymax,xv,yv)
     implicit none
     real(wp_), intent(in) :: xmin,xmax,ymin,ymax
-    real(wp_), intent(out), dimension(:), allocatable :: xv,yv
+    real(wp_), intent(out), dimension(5) :: xv,yv
-    if (allocated(xv)) deallocate(xv)
-    if (allocated(yv)) deallocate(yv)
-    allocate(xv(5),yv(5))
     xv=(/xmin,xmax,xmax,xmin,xmin/)
     yv=(/ymin,ymin,ymax,ymax,ymin/)
   end subroutine range2rect

srcjetto/gray.f

@@ -18,6 +18,8 @@
       integer ierr
 !   local variables
       integer i,j
+      real*8 tescal(nrho), dnescal(nrho)
+      real*8 p0mw
 
 ! === input arguments ==================================================
 !
@@ -87,9 +89,12 @@
       pec = 0.d0
       icd = 0.d0
       
+      dnescal = dne * 1.d-19
+      tescal = te * 1.d-3
 !     loop over beams with power>0
       do j=1,nbeam
-        if (powin(j).gt.0.0d0) cycle
+        if (powin(j).le.0.0d0) cycle
+        p0mw=powin(j)*1.d-6
 
 !       read j-th beam properties from file
 !       and adjust alpha/beta if out of the allowed range
@@ -97,17 +102,17 @@
 !       call main subroutine for the j-th beam
         call gray_jetto1beam(ijetto, mr, mz, r, z, psin(1:mr,:),
      .    psibnd-psiax, rax, zax, nbnd, rbnd, zbnd, nrho, psijet, -f,
-     .    te, dne, zeff, -qsf, j, powin(j), alphin(j), betain(j),
+     .    tescal, dnescal, zeff, -qsf, j, p0mw, alphin(j), 
-     .    dpdvloop, jcdloop, pecloop, icdloop, ierr)
+     .    betain(j), dpdvloop, jcdloop, pecloop, icdloop, ierr)
 
 !       add contribution of j-th beam to the total
 !       adapting output data to JETTO convention on toroidal angle
         do i=1,nrho
-          dpdv(i) = dpdv(i) + dpdvloop(i)
+          dpdv(i) = dpdv(i) + dpdvloop(i)*1.d6
-          jcd(i)  = jcd(i)  - jcdloop(i)
+          jcd(i)  = jcd(i)  - jcdloop(i)*1.d6
         end do
-        pec = pec + pecloop
+        pec = pec + pecloop*1.d6
-        icd = icd - icdloop
+        icd = icd - icdloop*1.d6
 
 !     end of loop over beams with power>0
       end do