gray/src/main.f90

program main_std
  use const_and_precisions, only : wp_,one
  use graycore,             only : gray_main, 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
  real(wp_) :: fghz
  real(wp_) :: x0, y0, z0, w1, w2, ri1, ri2, phiw, phir

  real(wp_) :: pec,icd

  integer :: ierr
  real(wp_), dimension(:), allocatable :: xrad, rhot, dpdv, jcd
  real(wp_) :: rwallm, rmxm, r0m, z0m, dzmx

  logical :: sum_mode = .FALSE.

  ! ------- sum mode variables -------
  real(wp_), dimension(:), allocatable :: jphi, currins, pins, rtin, rpin
  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
  ! ----------------------------------

  call read_params('gray_params.data', rtrp, hcdp, antp, eqp, rwallm, prfp, outp)

  ! ======= read input data BEGIN =======

  !------------ equilibrium ------------
  if(eqp%iequil<2) then
    call read_equil_an(eqp%filenm, rtrp%ipass, rv, zv, fpol, qpsi, rlim, zlim)
  ! 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
  ! ----------------------------------

  !--------------- wall ---------------
  ! set simple limiter if not read from EQDSK
  ! need to clean up...
    r0m=sqrt(x0**2+y0**2)*0.01_wp_
    dzmx=abs(rtrp%ipass)*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(dpdv(outp%nrho), jcd(outp%nrho))
  if (sum_mode) then
    allocate(jphi(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,                     &
                        antp%beta,fghz,antp%alpha,antp%beta,                        &
                        (/x0,y0,z0/),w1,w2,ri1,ri2,phiw,phir,antp%iox,              &
                        antp%psi, antp%chi, 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
  else
    call gray_main(rv,zv,psin,psia,psinr,fpol,qpsi,rvac,rax,zax,rbnd,zbnd,eqp, &
                   psrad,terad,derad,zfc,prfp,  rlim,zlim,                     &
                   antp%power,fghz,antp%alpha,antp%beta,                       &
                   (/x0,y0,z0/),w1,w2,ri1,ri2,phiw,phir,antp%iox,              &
                   antp%psi,antp%chi, dpdv,jcd,pec,icd, outp,rtrp,hcdp,ierr)
  end if
  ! ========================================================================

  ! ======= 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(dpdv, jcd)
  if(sum_mode)         deallocate(jphi, currins, pins, rtin, rpin)
  ! ======= free memory  END  ======
end program main_std