src/coreprofiles.f90: use derived types

src/coreprofiles.f90

@@ -2,11 +2,11 @@ module coreprofiles
   use const_and_precisions, only : wp_,zero,one
   implicit none
 
-  INTEGER, SAVE :: npp,nsfd
+  integer, save :: npp, nsfd
-  REAL(wp_), SAVE :: psdbnd,psnpp,denpp,ddenpp,d2denpp
+  real(wp_), save :: psdbnd, psnpp, denpp, ddenpp, d2denpp
-  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: tfn,cfn,psrad
+  real(wp_), dimension(:),  allocatable,  save :: tfn, cfn, psrad
-  REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: ct,cz
+  real(wp_), dimension(:, :),  allocatable,  save :: ct, cz
-  REAL(wp_), SAVE :: dens0,aln1,aln2,te0,dte0,alt1,alt2,zeffan
+  real(wp_), save :: dens0, aln1, aln2, te0, dte0, alt1, alt2, zeffan
 
 contains
 
@@ -125,35 +125,57 @@ contains
     endif
   end function fzeff
 
-  subroutine read_profiles(filenm,psin,te,ne,zeff,unit)
+  
-    use utils, only : get_free_unit
+  subroutine read_profiles(filenm, data, unit)
+    ! Reads the radial plasma profiles from `file` and store them
+    ! into `data`. If given, the file is opened in the `unit` number.
+    ! Format notes:
+    ! 1. The file is formatted as a table with the following columns:
+    !   radial coordinate, temperature, density, effective charge.
+    ! 2. The first line is a header specifying the number of rows.
+    use utils,       only : get_free_unit
+    use gray_params, only : profiles_data
+
     implicit none
-! arguments
+
-    character(len=*), intent(in) :: filenm
+    ! subroutine arguments
-    real(wp_), dimension(:), allocatable, intent(out) :: psin,te,ne,zeff
+    character(len=*),    intent(in)  :: filenm
-    integer, optional, intent(in) :: unit
+    type(profiles_data), intent(out) :: data
-! local variables
+    integer, optional,   intent(in)  :: unit
-    integer :: u, i, n
+
+    ! local variables
+    integer :: u, i, nrows
+
+    ! Free the arrays when already allocated
+    if(allocated(data%psrad)) deallocate(data%psrad)
+    if(allocated(data%terad)) deallocate(data%terad)
+    if(allocated(data%derad)) deallocate(data%derad)
+    if(allocated(data%zfc))   deallocate(data%zfc)
 
     if (present(unit)) then
-      u=unit
+      u = unit
     else
-      u=get_free_unit()
+      u = get_free_unit()
     end if
-    open(file=trim(filenm),status='old',action='read',unit=u)
+
-    read(u,*) n
+    ! Read number of rows and allocate the arrays
-    if(allocated(psin)) deallocate(psin)
+    open(file=trim(filenm), status='old', action='read', unit=u)
-    if(allocated(te)) deallocate(te)
+    read(u, *) nrows
-    if(allocated(ne)) deallocate(ne)
+    allocate(data%psrad(nrows), data%terad(nrows), &
-    if(allocated(zeff)) deallocate(zeff)
+             data%derad(nrows), data%zfc(nrows))
-    allocate(psin(n),te(n),ne(n),zeff(n))
+
-    do i=1,n
+    ! Read the table rows
-      read(u,*) psin(i),te(i),ne(i),zeff(i)
+    do i=1,nrows
+      read(u, *) data%psrad(i), data%terad(i), data%derad(i), data%zfc(i)
     end do
-    psin(1)=max(psin(1),zero)
     close(u)
+
+    ! ??
+    data%psrad(1) = max(data%psrad(1), zero)
+
   end subroutine read_profiles
 
+
   subroutine read_profiles_an(filenm,te,ne,zeff,unit)
     use utils, only : get_free_unit
      implicit none
@@ -182,47 +204,66 @@ contains
     close(u)
   end subroutine read_profiles_an
 
-  subroutine tene_scal(te,ne,tfact,nfact,bfact,iscal,iprof)
+
+  subroutine scale_profiles(params, factb, data)
+    ! Rescales the temperature and density profiles by `params%factte`
+    ! and `params%factne`, respectively, according to the model
+    ! specified by `params%iscal`.
+    ! See the GRAY user manual for the explanation.
+    use gray_params, only : profiles_parameters, profiles_data
+
     implicit none
-! arguments
-    real(wp_), dimension(:), intent(inout) :: te,ne
-    real(wp_), intent(in) :: tfact,nfact,bfact
-    integer, intent(in) :: iscal,iprof
-! local variables
-    real(wp_) :: aat,aan,ffact
-    integer :: lastte,lastne
 
-    if (iscal==0) then
+    ! subroutine arguments
-      aat=2.0_wp_/3.0_wp_
+    type(profiles_parameters), intent(in)    :: params
-      aan=4.0_wp_/3.0_wp_
+    real(wp_),                 intent(in)    :: factb
-    else
+    type(profiles_data),       intent(inout) :: data
-      aat=1.0_wp_
-      aan=1.0_wp_
-    end if
-    if(iscal==2) then
-      ffact=1.0_wp_
-    else
-      ffact=bfact
-    end if
-    if (iprof==0) then
-      lastte=2
-      lastne=1
-    else
-      lastte=size(te)
-      lastne=size(ne)
-    end if
-    te(1:lastte)=te(1:lastte)*ffact**aat*tfact
-    ne(1:lastne)=ne(1:lastne)*ffact**aan*nfact
-  end subroutine tene_scal
 
-  subroutine set_prfspl(psin,te,ne,zeff,ssplne,psdbndmx)
+    ! local variables
+    real(wp_) :: aat, aan, ffact
+    integer   :: last_te, last_ne
+
+    if (params%iscal==0) then
+      aat = 2.0_wp_/3.0_wp_
+      aan = 4.0_wp_/3.0_wp_
+    else
+      aat = one
+      aan = one
+    end if
+
+    if(params%iscal==2) then
+      ffact = one
+    else
+      ffact = factb
+    end if
+
+    if(params%iprof==0) then
+      last_te = 2
+      last_ne = 1
+    else
+      last_te = size(data%terad)
+      last_ne = size(data%derad)
+    end if
+
+    data%terad(1:last_te) = data%terad(1:last_te) * ffact**aat * params%factte
+    data%derad(1:last_ne) = data%derad(1:last_ne) * ffact**aan * params%factne
+  end subroutine scale_profiles
+
+
+  subroutine set_prfspl(params, data)
+    ! Computes splines for the plasma profiles data and stores them
+    ! in their respective global variables, see the top of this file.
     use simplespline, only : difcs
-    use dierckx, only : curfit, splev, splder
+    use dierckx,      only : curfit, splev, splder
+    use gray_params,  only : profiles_parameters, profiles_data
+
     implicit none
-! arguments
+
-    real(wp_), dimension(:), intent(in) :: psin,te,ne,zeff
+    ! subroutine arguments
-    real(wp_), intent(in) :: ssplne,psdbndmx
+    type(profiles_parameters), intent(in)    :: params
-! local variables
+    type(profiles_data),       intent(inout) :: data
+
+    ! local variables
     integer, parameter :: iopt=0, kspl=3
     integer :: n, npest, lwrkf, ier
     real(wp_) :: xb, xe, fp, xnv, xxp,xxm,delta2,ssplne_loc
@@ -230,14 +271,14 @@ contains
     integer, dimension(:), allocatable :: iwrkf
     real(wp_), dimension(1) :: dedge,ddedge,d2dedge
 
-    n=size(psin)
+    n=size(data%psrad)
     npest=n+4
     lwrkf=n*4+npest*16
     allocate(wrkf(lwrkf),iwrkf(npest),wf(n))
 
-    ssplne_loc=ssplne
+    ssplne_loc=params%sspld
 
-!   if necessary, reallocate spline arrays
+    ! if necessary, reallocate spline arrays
     if(.not.allocated(psrad)) then
       allocate(psrad(n),ct(n,4),cz(n,4))
     else
@@ -255,37 +296,37 @@ contains
       end if
     end if
 
-!   spline approximation of  temperature and Zeff
+    ! spline approximation of  temperature and data%zfc
-    call difcs(psin,te,  n,iopt,ct,ier)
+    call difcs(data%psrad,data%terad,  n,iopt,ct,ier)
-    call difcs(psin,zeff,n,iopt,cz,ier)
+    call difcs(data%psrad,data%zfc,n,iopt,cz,ier)
-    psrad=psin
+    psrad=data%psrad
     npp=n
 
-!   spline approximation of density
+    ! spline approximation of density
     xb=zero
-    xe=psin(n)
+    xe=data%psrad(n)
     wf(:)=one
-    call curfit(iopt,n,psin,ne,wf,xb,xe,kspl,ssplne_loc,npest, &
+    call curfit(iopt,n,data%psrad,data%derad,wf,xb,xe,kspl,ssplne_loc,npest, &
        nsfd,tfn,cfn,fp,wrkf,lwrkf,iwrkf,ier)
-! if ier=-1 data are re-fitted using sspl=0
+    ! if ier=-1 data are re-fitted using sspl=0
     if(ier==-1) then
         write(*,*) 'density curfit: ier=-1. Re-fitting with interpolating spline'
         ssplne_loc=0.0_wp_
-        call curfit(iopt,n,psin,ne,wf,xb,xe,kspl,ssplne_loc,npest, &
+        call curfit(iopt,n,data%psrad,data%derad,wf,xb,xe,kspl,ssplne_loc,npest, &
            nsfd,tfn,cfn,fp,wrkf,lwrkf,iwrkf,ier)
     end if
 
-!   compute polinomial extrapolation matching the spline boundary up to the
+    ! compute polinomial extrapolation matching the spline boundary up to the
-!   2nd order derivative, extending the profile up to psi=psdbnd where
+    ! 2nd order derivative, extending the profile up to psi=psdbnd where
-!   ne=ne'=ne''=0
+    ! data%derad=data%derad'=data%derad''=0
-!   spline value and derivatives at the edge
+    ! spline value and derivatives at the edge
-    call splev(tfn,nsfd,cfn,kspl,psin(n:n),dedge(1:1),1,ier)
+    call splev(tfn,nsfd,cfn,kspl,data%psrad(n:n),dedge(1:1),1,ier)
-    call splder(tfn,nsfd,cfn,kspl,1,psin(n:n),ddedge(1:1), 1,wrkf(1:nsfd),ier)
+    call splder(tfn,nsfd,cfn,kspl,1,data%psrad(n:n),ddedge(1:1), 1,wrkf(1:nsfd),ier)
-    call splder(tfn,nsfd,cfn,kspl,2,psin(n:n),d2dedge(1:1),1,wrkf(1:nsfd),ier)
+    call splder(tfn,nsfd,cfn,kspl,2,data%psrad(n:n),d2dedge(1:1),1,wrkf(1:nsfd),ier)
-!   determination of the boundary
+    ! determination of the boundary
-    psdbnd=psdbndmx
+    psdbnd=params%psnbnd
 
-    psnpp=psin(n)
+    psnpp=data%psrad(n)
     denpp=dedge(1)
     ddenpp=ddedge(1)
     d2denpp=d2dedge(1)
@@ -293,7 +334,7 @@ contains
     delta2=(ddenpp/d2denpp)**2-2.0_wp_*denpp/d2denpp
     xnv=psnpp-ddenpp/d2denpp
     if(delta2 < zero) then
-!       if(xnv > psnpp) psdbnd=min(psdbnd,xnv)
+    ! if(xnv > psnpp) psdbnd=min(psdbnd,xnv)
     else
       xxm=xnv-sqrt(delta2)
       xxp=xnv+sqrt(delta2)
@@ -302,12 +343,13 @@ contains
       else if (xxp > psnpp) then
         psdbnd=min(psdbnd,xxp)
       end if    
-      write(*,*) 'density psdbnd=',psdbnd
+      write(*,'(a,f5.3)') 'density psdbnd=', psdbnd
     end if
 
     deallocate(iwrkf,wrkf,wf)
   end subroutine set_prfspl
 
+
   subroutine unset_prfspl
     implicit none
 

src/main.f90

@@ -21,7 +21,7 @@ program main
   ! Read the input data into set the global variables
   ! of the respective module. Note: order matters.
   call init_equilibrium(params, data)
-  call init_profiles(params, data)
+  call init_profiles(params%profiles, params%equilibrium%factb, data%profiles)
   call init_antenna(params%antenna)
   call init_misc(params, data)
 
@@ -182,63 +182,50 @@ contains
   end subroutine deinit_equilibrium
 
 
-  subroutine init_profiles(params, data)
+  subroutine init_profiles(params, factb, data)
     ! Reads the plasma kinetic profiles file (containing the elecron
     ! temperature, density and plasma effective charge) and initialises
     ! the respective GRAY data structure.
     use gray_params,  only : profiles_parameters, profiles_data
     use equilibrium,  only : frhopolv
-    use coreprofiles, only : read_profiles_an, read_profiles, tene_scal, &
+    use coreprofiles, only : read_profiles_an, read_profiles, &
-                             set_prfan, set_prfspl
+                             scale_profiles, set_prfan, set_prfspl
     implicit none
 
     ! subroutine arguments
-    type(gray_parameters), intent(in)            :: params
+    type(profiles_parameters), intent(in)  :: params
-    type(gray_data),       intent(inout), target :: data
+    real(wp_),                 intent(in)  :: factb
+    type(profiles_data),       intent(out) :: data
 
-    ! local variables
+    if(params%iprof == 0) then
-    type(profiles_parameters) :: profp
-    type(profiles_data), pointer :: profd
-
-    ! Radial coordinate (depending on profp%irho: ρ_t, ρ_p, or ψ)
-    real(wp_), allocatable :: xrad(:)
-
-    profp = params%profiles
-    profd => data%profiles
-
-    if(params%profiles%iprof == 0) then
       ! Analytical profiles
       ! TODO: rewrite using derived type
-      call read_profiles_an(profp%filenm, profd%terad, profd%derad, profd%zfc)
+      call read_profiles_an(params%filenm, data%terad, data%derad, data%zfc)
     else
       ! Numerical profiles
-      call read_profiles(profp%filenm, xrad, profd%terad, profd%derad, profd%zfc)
+      call read_profiles(params%filenm, data)
 
-      allocate(profd%psrad(size(xrad)))
+      ! Convert psrad to ψ
-
+      select case (params%irho)
-      select case (profp%irho)
+        case (0) ! psrad is ρ_t
-        case (0) ! xrad is rhot
+          data%psrad = frhopolv(data%psrad)**2
-          profd%psrad = frhopolv(xrad)**2
+        case (1) ! psrad is ρ_p
-        case (1) ! xrad is rhop
+          data%psrad = data%psrad**2
-          profd%psrad = xrad**2
+        case default ! psrad is already ψ
-        case default ! xrad is psi
-          profd%psrad = xrad
       end select
-      deallocate(xrad)
-
     end if
+
     ! Rescale input data
-    ! TODO: rewrite using derived type
+    call scale_profiles(params, factb, data)
-    call tene_scal(profd%terad, profd%derad, profp%factte, profp%factne, &
-                   params%equilibrium%factb, profp%iscal, profp%iprof)
 
     ! Set global variables
-    ! TODO: rewrite using derived type
+    if(params%iprof == 0) then
-    if(params%profiles%iprof == 0) then
+      ! Analytical profiles
-      call set_prfan(profd%terad, profd%derad, profd%zfc)
+      ! TODO: rewrite using derived type
+      call set_prfan(data%terad, data%derad, data%zfc)
     else
-      call set_prfspl(profd%psrad, profd%terad, profd%derad, profd%zfc, &
+      ! Numerical profiles
-                      profp%sspld, profp%psnbnd)
+      call set_prfspl(params, data)
     end if
   end subroutine init_profiles