From f56e1cbc054830cc31bdabb7bea2585c3f38136e Mon Sep 17 00:00:00 2001 From: Michele Guerini Rocco Date: Wed, 15 Dec 2021 02:31:10 +0100 Subject: [PATCH] src/coreprofiles.f90: use derived types --- src/coreprofiles.f90 | 214 ++++++++++++++++++++++++++----------------- src/main.f90 | 63 +++++-------- 2 files changed, 153 insertions(+), 124 deletions(-) diff --git a/src/coreprofiles.f90 b/src/coreprofiles.f90 index d882b92..7b1f4f3 100644 --- a/src/coreprofiles.f90 +++ b/src/coreprofiles.f90 @@ -2,11 +2,11 @@ module coreprofiles use const_and_precisions, only : wp_,zero,one implicit none - INTEGER, SAVE :: npp,nsfd - REAL(wp_), SAVE :: psdbnd,psnpp,denpp,ddenpp,d2denpp - REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: tfn,cfn,psrad - REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: ct,cz - REAL(wp_), SAVE :: dens0,aln1,aln2,te0,dte0,alt1,alt2,zeffan + integer, save :: npp, nsfd + real(wp_), save :: psdbnd, psnpp, denpp, ddenpp, d2denpp + real(wp_), dimension(:), allocatable, save :: tfn, cfn, psrad + real(wp_), dimension(:, :), allocatable, save :: ct, cz + real(wp_), save :: dens0, aln1, aln2, te0, dte0, alt1, alt2, zeffan contains @@ -124,36 +124,58 @@ contains fzeff=spli(cz,npp,k,dps) 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 - real(wp_), dimension(:), allocatable, intent(out) :: psin,te,ne,zeff - integer, optional, intent(in) :: unit -! local variables - integer :: u, i, n - + + ! subroutine arguments + character(len=*), intent(in) :: filenm + type(profiles_data), intent(out) :: data + integer, optional, intent(in) :: unit + + ! 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 - if(allocated(psin)) deallocate(psin) - if(allocated(te)) deallocate(te) - if(allocated(ne)) deallocate(ne) - if(allocated(zeff)) deallocate(zeff) - allocate(psin(n),te(n),ne(n),zeff(n)) - do i=1,n - read(u,*) psin(i),te(i),ne(i),zeff(i) + + ! Read number of rows and allocate the arrays + open(file=trim(filenm), status='old', action='read', unit=u) + read(u, *) nrows + allocate(data%psrad(nrows), data%terad(nrows), & + data%derad(nrows), data%zfc(nrows)) + + ! Read the table rows + 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 - aat=2.0_wp_/3.0_wp_ - aan=4.0_wp_/3.0_wp_ - else - 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 arguments + type(profiles_parameters), intent(in) :: params + real(wp_), intent(in) :: factb + type(profiles_data), intent(inout) :: data - 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 - real(wp_), intent(in) :: ssplne,psdbndmx -! local variables + + ! subroutine arguments + type(profiles_parameters), intent(in) :: params + 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 - call difcs(psin,te, n,iopt,ct,ier) - call difcs(psin,zeff,n,iopt,cz,ier) - psrad=psin + ! spline approximation of temperature and data%zfc + call difcs(data%psrad,data%terad, n,iopt,ct,ier) + call difcs(data%psrad,data%zfc,n,iopt,cz,ier) + 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 -! 2nd order derivative, extending the profile up to psi=psdbnd where -! ne=ne'=ne''=0 -! spline value and derivatives at the edge - call splev(tfn,nsfd,cfn,kspl,psin(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,2,psin(n:n),d2dedge(1:1),1,wrkf(1:nsfd),ier) -! determination of the boundary - psdbnd=psdbndmx + ! compute polinomial extrapolation matching the spline boundary up to the + ! 2nd order derivative, extending the profile up to psi=psdbnd where + ! data%derad=data%derad'=data%derad''=0 + ! spline value and derivatives at the edge + call splev(tfn,nsfd,cfn,kspl,data%psrad(n:n),dedge(1:1),1,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,data%psrad(n:n),d2dedge(1:1),1,wrkf(1:nsfd),ier) + ! determination of the boundary + 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 diff --git a/src/main.f90 b/src/main.f90 index e069374..aa9b7c5 100644 --- a/src/main.f90 +++ b/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, & - set_prfan, set_prfspl + use coreprofiles, only : read_profiles_an, read_profiles, & + scale_profiles, set_prfan, set_prfspl implicit none ! subroutine arguments - type(gray_parameters), intent(in) :: params - type(gray_data), intent(inout), target :: data + type(profiles_parameters), intent(in) :: params + real(wp_), intent(in) :: factb + type(profiles_data), intent(out) :: data - ! local variables - 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 + if(params%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))) - - select case (profp%irho) - case (0) ! xrad is rhot - profd%psrad = frhopolv(xrad)**2 - case (1) ! xrad is rhop - profd%psrad = xrad**2 - case default ! xrad is psi - profd%psrad = xrad + ! Convert psrad to ψ + select case (params%irho) + case (0) ! psrad is ρ_t + data%psrad = frhopolv(data%psrad)**2 + case (1) ! psrad is ρ_p + data%psrad = data%psrad**2 + case default ! psrad is already ψ end select - deallocate(xrad) - end if + ! Rescale input data - ! TODO: rewrite using derived type - call tene_scal(profd%terad, profd%derad, profp%factte, profp%factne, & - params%equilibrium%factb, profp%iscal, profp%iprof) + call scale_profiles(params, factb, data) ! Set global variables - ! TODO: rewrite using derived type - if(params%profiles%iprof == 0) then - call set_prfan(profd%terad, profd%derad, profd%zfc) + if(params%iprof == 0) then + ! Analytical profiles + ! 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, & - profp%sspld, profp%psnbnd) + ! Numerical profiles + call set_prfspl(params, data) end if end subroutine init_profiles