From add59dbdda15a747ba5e8162362b54ede6d8af2d Mon Sep 17 00:00:00 2001 From: Michele Guerini Rocco Date: Sat, 18 Dec 2021 18:57:38 +0100 Subject: [PATCH] src: use the logging system everywhere --- src/beams.f90 | 29 ++++++++- src/conical.f90 | 24 ++++--- src/coreprofiles.f90 | 81 +++++++++++++++-------- src/equilibrium.f90 | 149 +++++++++++++++++++++++++++++++------------ src/gray_cli.f90 | 18 +++--- src/gray_core.f90 | 144 +++++++++++++++++++++++++---------------- src/gray_params.f90 | 21 ++++-- src/main.f90 | 40 +++++++++--- src/multipass.f90 | 9 +-- 9 files changed, 352 insertions(+), 163 deletions(-) diff --git a/src/beams.f90 b/src/beams.f90 index c6241ee..f3b6d42 100644 --- a/src/beams.f90 +++ b/src/beams.f90 @@ -11,6 +11,7 @@ contains use const_and_precisions, only : pi, vc=>ccgs_ use gray_params, only : antenna_parameters use utils, only : get_free_unit + use logger, only : log_error implicit none @@ -20,11 +21,18 @@ contains ! local variables integer :: u + integer :: err real(wp_) :: ak0,zrcsi,zreta u = get_free_unit(unit) - open(unit=u, file=trim(params%filenm), status='OLD', action='READ') + open(unit=u, file=params%filenm, status='old', action='read', iostat=err) + if (err /= 0) then + call log_error('opening beams file ('//trim(params%filenm)//') failed!', & + mod='beams', proc="read_beam0") + call exit(1) + end if + read(u, *) params%fghz read(u, *) params%pos read(u, *) params%w, params%ri, params%phi(1) @@ -49,6 +57,7 @@ contains use gray_params, only : antenna_parameters use simplespline, only : spli, difcs use utils, only : get_free_unit,locate + use logger, only : log_error implicit none @@ -64,10 +73,16 @@ contains z00v, waist1v, waist2v, rci1v, rci2v, phi1v, phi2v, & cbeta, cx0, cy0, cz0, cwaist1, cwaist2, & crci1, crci2, cphi1, cphi2 + integer :: err u = get_free_unit(unit) - open(unit=u, file=params%filenm, status='OLD', action='READ') + open(unit=u, file=params%filenm, status='old', action='read', iostat=err) + if (err /= 0) then + call log_error('opening beams file ('//trim(params%filenm)//') failed!', & + mod='beams', proc="read_beam1") + call exit(1) + end if read(u,*) params%fghz read(u,*) nisteer @@ -155,6 +170,7 @@ contains use utils, only : get_free_unit, intlin, locate use reflections, only : inside use dierckx, only : curfit, splev, surfit, bispev + use logger, only : log_error implicit none @@ -187,10 +203,17 @@ contains real(wp_), dimension(1) :: fi integer, parameter :: kspl=1 real(wp_), parameter :: sspl=0.01_wp_ + integer :: err u = get_free_unit(unit) - open(unit=u, file=params%filenm, status='OLD', action='READ') + open(unit=u, file=params%filenm, status='old', action='read', iostat=err) + if (err /= 0) then + call log_error('opening beams file ('//trim(params%filenm)//') failed!', & + mod='beams', proc="read_beam1") + call exit(1) + end if + !======================================================================================= ! # of beams read(u,*) nbeam diff --git a/src/conical.f90 b/src/conical.f90 index fb117d2..596e52b 100644 --- a/src/conical.f90 +++ b/src/conical.f90 @@ -29,12 +29,14 @@ contains complex(wp_) v0,v1,v2,vv,w(19) logical lm0,lm1,lta + character(256) :: msg fconic=0.0_wp_ lm0=m == 0 lm1=m == 1 if(.not.(lm0 .or. lm1)) then - write(*,"(1x,'fconic ... illegal value for m = ',i4)") m + write (msg, '("invalid value for m: ",g0)') m + call log_error(msg, mod='conical', proc='fconic') return end if fm=m @@ -204,7 +206,9 @@ contains do n=n+1 if(n > nmax) then - write(*,200) x,tau,m + write (msg, '(a," x=",g0.4," tau=",g0.4," m=",g0)') & + "convergence difficulties for c function:", x, tau, m + call log_error(msg, mod='conical', proc='fconical') return end if rr=r @@ -258,7 +262,9 @@ contains if(abs(r-rr) < eps) exit end do if (n > nmax) then - write(*,200) x,tau,m + write (msg, '(a," x=",g0.4," tau=",g0.4," m=",g0)') & + "convergence difficulties for c function:", x, tau, m + call log_error(msg, mod='conical', proc='fconical') return end if end if @@ -287,10 +293,6 @@ contains return end if end if -! - 200 format(1x,'fconic ... convergence difficulties for c function, x = ', & - e12.4,5x,'tau = ',e12.4,5x,'m = ',i5) -! end function fconic ! @@ -312,11 +314,10 @@ contains +8.4175084175084e-4_wp_, -1.9175269175269e-3_wp_, & +6.4102564102564e-3_wp_, -2.9550653594771e-2_wp_, & +1.7964437236883e-1_wp_, -1.3924322169059e+0_wp_/) -! + x=real(z) t=aimag(z) if(-abs(x) == aint(x) .and. t == 0.0_wp_) then - write(*,'(1x,f20.2)') x clogam=(0.0_wp_,0.0_wp_) return end if @@ -360,7 +361,12 @@ contains end function clogam function ellick(xk) + ! Computes the complete elliptic integrals K(x), E(x): + ! entry ellick(X)= E(x) + ! entry ellice(X)= E(x) + implicit none + real(wp_), intent(in) :: xk real(wp_) :: ellick, ellice integer :: i diff --git a/src/coreprofiles.f90 b/src/coreprofiles.f90 index a0a7cc6..2bad8d9 100644 --- a/src/coreprofiles.f90 +++ b/src/coreprofiles.f90 @@ -12,24 +12,29 @@ contains subroutine density(psin,dens,ddens) use gray_params, only : iprof - use dierckx, only : splev,splder + use dierckx, only : splev,splder + use logger, only : log_error + implicit none -! arguments + + ! subroutine arguments real(wp_), intent(in) :: psin real(wp_), intent(out) :: dens,ddens -! local variables + + ! local variables integer :: ier,nu real(wp_) :: profd,dprofd,dpsib,tt,fp,dfp,fh,dfh real(wp_), dimension(1) :: xxs,ffs real(wp_), dimension(npp+4) :: wrkfd + character(256) :: msg -! -! computation of density [10^19 m^-3] and derivative wrt psi -! + ! + ! Computation of density [10¹⁹ m⁻³] and derivative wrt ψ + ! dens=zero ddens=zero - if((psin >= psdbnd).or.(psin < zero)) return -! + if((psin >= psdbnd) .or. (psin < zero)) return + if(iprof == 0) then if(psin > one) return profd=(one-psin**aln1)**aln2 @@ -40,12 +45,12 @@ contains else if(psin > psnpp) then -! smooth interpolation for psnpp < psi < psdbnd -! dens = fp * fh -! fp: parabola matched at psi=psnpp with given profile density -! fh=(1-t)^3(1+3t+6t^2) is a smoothing function: -! fh(0)=1, fh(1)=0 and zero first and second deriv at t=0,1 -! + ! Smooth interpolation for psnpp < psi < psdbnd + ! dens = fp * fh + ! fp: parabola matched at psi=psnpp with given profile density + ! fh=(1-t)^3(1+3t+6t^2) is a smoothing function: + ! fh(0)=1, fh(1)=0 and zero first and second deriv at t=0,1 + ! dpsib=psin-psnpp fp=denpp+dpsib*ddenpp+0.5_wp_*dpsib**2*d2denpp dfp=ddenpp+dpsib*d2denpp @@ -65,14 +70,15 @@ contains ddens=ffs(1) if(abs(dens) < 1.0e-10_wp_) dens=zero end if - if(dens < zero) print*,'psin = ',psin,': DENSITY NEGATIVE ne=',dens -! if(dens < zero) then -! dens=zero -! ddens=zero -! end if + if(dens < zero) then + write (msg, '("negative density:", 2(x,a,"=",g0.3))') & + 'ne', dens, 'ψ', psin + call log_error(msg, mod='coreprofiles', proc='density') + end if end if end subroutine density + function temp(psin) use const_and_precisions, only : wp_,zero,one use gray_params, only : iprof @@ -134,6 +140,7 @@ contains ! 2. The first line is a header specifying the number of rows. use utils, only : get_free_unit use gray_params, only : profiles_data + use logger, only : log_error implicit none @@ -144,6 +151,7 @@ contains ! local variables integer :: u, i, nrows + integer :: err ! Free the arrays when already allocated if(allocated(data%psrad)) deallocate(data%psrad) @@ -154,7 +162,13 @@ contains u = get_free_unit(unit) ! Read number of rows and allocate the arrays - open(file=trim(filenm), status='old', action='read', unit=u) + open(file=filenm, status='old', action='read', unit=u, iostat=err) + if (err /= 0) then + call log_error('opening profiles file ('//trim(filenm)//') failed!', & + mod='coreprofiles', proc="read_profiles") + call exit(1) + end if + read(u, *) nrows allocate(data%psrad(nrows), data%terad(nrows), & data%derad(nrows), data%zfc(nrows)) @@ -172,14 +186,17 @@ contains subroutine read_profiles_an(filenm,te,ne,zeff,unit) - use utils, only : get_free_unit - implicit none + use utils, only : get_free_unit + use logger, only : log_error + + implicit none ! arguments character(len=*), intent(in) :: filenm real(wp_), dimension(:), allocatable, intent(out) :: te,ne,zeff integer, optional, intent(in) :: unit ! local variables integer :: u + integer :: err u = get_free_unit(unit) @@ -188,7 +205,13 @@ contains if(allocated(zeff)) deallocate(zeff) allocate(te(4),ne(3),zeff(1)) - open(file=trim(filenm),status='old',action='read',unit=u) + open(file=filenm, status='old', action='read', unit=u, iostat=err) + if (err /= 0) then + call log_error('opening profiles file ('//trim(filenm)//') failed!', & + mod='coreprofiles', proc='read_profiles_an') + call exit(1) + end if + read(u,*) ne(1:3) ! dens0,aln1,aln2 read(u,*) te(1:4) ! te0,dte0,alt1,alt2 read(u,*) zeff(1) ! zeffan @@ -247,6 +270,7 @@ contains use simplespline, only : difcs use dierckx, only : curfit, splev, splder use gray_params, only : profiles_parameters, profiles_data + use logger, only : log_info, log_warning implicit none @@ -261,6 +285,7 @@ contains real(wp_), dimension(:), allocatable :: wf, wrkf integer, dimension(:), allocatable :: iwrkf real(wp_), dimension(1) :: dedge,ddedge,d2dedge + character(256) :: msg ! for log messages formatting n=size(data%psrad) npest=n+4 @@ -300,8 +325,9 @@ contains 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) then - write(*,*) 'density curfit: ier=-1. Re-fitting with interpolating spline' + if(ier == -1) then + call log_warning('curfit failed with error -1: re-fitting with '// & + 's=0', mod='coreprofiles', proc='density') ssplne_loc=0.0_wp_ call curfit(iopt,n,data%psrad,data%derad,wf,xb,xe,kspl,ssplne_loc,npest, & nsfd,tfn,cfn,fp,wrkf,lwrkf,iwrkf,ier) @@ -333,8 +359,9 @@ contains psdbnd=min(psdbnd,xxm) else if (xxp > psnpp) then psdbnd=min(psdbnd,xxp) - end if - write(*,'(a,f5.3)') 'density psdbnd=', psdbnd + end if + write (msg, '(a,g0.3)') 'density boundary: ψ=', psdbnd + call log_info(msg, mod="coreprofiles", proc="set_prfspl") end if deallocate(iwrkf,wrkf,wf) diff --git a/src/equilibrium.f90 b/src/equilibrium.f90 index 2c2be51..21fbcb1 100644 --- a/src/equilibrium.f90 +++ b/src/equilibrium.f90 @@ -41,6 +41,7 @@ contains use const_and_precisions, only : one use gray_params, only : equilibrium_parameters, equilibrium_data use utils, only : get_free_unit + use logger, only : log_error implicit none @@ -54,11 +55,17 @@ contains character(len=48) :: string real(wp_) :: dr, dz, dps, rleft, zmid, zleft, psiedge, psiaxis real(wp_) :: xdum ! dummy variable, used to discard data + integer :: err u = get_free_unit(unit) ! Open the G-EQDSK file - open(file=trim(params%filenm), status='old', action='read', unit=u) + open(file=params%filenm, status='old', action='read', unit=u, iostat=err) + if (err /= 0) then + call log_error('opening eqdsk file ('//trim(params%filenm)//') failed!', & + mod='equilibrium', proc='read_eqdsk') + call exit(1) + end if ! get size of main arrays and allocate them if (params%idesc == 1) then @@ -154,7 +161,9 @@ contains subroutine read_equil_an(filenm,ipass,rv,zv,fpol,q,rlim,zlim,unit) - use utils, only : get_free_unit + use utils, only : get_free_unit + use logger, only : log_error + implicit none ! arguments character(len=*), intent(in) :: filenm @@ -163,11 +172,18 @@ contains real(wp_), dimension(:), allocatable, intent(out) :: rv,zv,fpol,q,rlim,zlim ! local variables integer :: i, u, nlim + integer :: err real(wp_) :: rr0m,zr0m,rpam,b0,q0,qa,alq !,rcen,btrcen u = get_free_unit(unit) - open(file=trim(filenm),status='old',action='read',unit=u) + open(file=filenm, status='old', action='read', unit=u, iostat=err) + if (err /= 0) then + call log_error('opening equilibrium file ('//trim(filenm)//') failed!', & + mod='equilibrium', proc='read_equil_an') + call exit(1) + end if + read(u,*) rr0m,zr0m,rpam read(u,*) b0 read(u,*) q0,qa,alq @@ -338,6 +354,7 @@ contains use gray_params, only : iequil use reflections, only : inside use utils, only : vmaxmin, vmaxmini + use logger, only : log_info implicit none @@ -358,6 +375,7 @@ contains real(wp_), dimension(:), allocatable :: rv1d,zv1d,fvpsi,wf,wrk integer, dimension(:), allocatable :: iwrk integer :: ier,ixploc,info,i,j,ij + character(256) :: msg ! for log messages formatting ! compute array sizes and prepare working space arrays nr=size(data%rv) @@ -526,7 +544,10 @@ contains rax0=data%rax zax0=data%zax call points_ox(rax0,zax0,rmaxis,zmaxis,psinoptmp,info) - print'(a,2f8.4,es12.5)','O-point',rmaxis,zmaxis,psinoptmp + + write (msg, '("O-point found:", 3(x,a,"=",g0.3))') & + 'r', rmaxis, 'z', zmaxis, 'ψ', psinoptmp + call log_info(msg, mod='equilibrium', proc='set_eqspl') ! search for X-point if params%ixp /= 0 @@ -535,7 +556,10 @@ contains if(ixploc<0) then call points_ox(rbinf,zbinf,r1,z1,psinxptmp,info) if(psinxptmp/=-1.0_wp_) then - print'(a,2f8.4,es12.5)','X-point',r1,z1,psinxptmp + write (msg, '("X-point found:", 3(x,a,"=",g0.3))') & + 'r', r1, 'z', z1, 'ψ', psinxptmp + call log_info(msg, mod='equilibrium', proc='set_eqspl') + zbinf=z1 psinop=psinoptmp psiant=psinxptmp-psinop @@ -547,7 +571,10 @@ contains else call points_ox(rbsup,zbsup,r1,z1,psinxptmp,info) if(psinxptmp.ne.-1.0_wp_) then - print'(a,2f8.4,e16.8)','X-point',r1,z1,psinxptmp + write (msg, '("X-point found:", 3(x,a,"=",g0.3))') & + 'r', r1, 'z', z1, 'ψ', psinxptmp + call log_info(msg, mod='equilibrium', proc='set_eqspl') + zbsup=z1 psinop=psinoptmp psiant=psinxptmp-psinop @@ -572,9 +599,10 @@ contains call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbinf),r1,z1,one,info) zbinf=z1 rbinf=r1 - print'(a,4f8.4)','no X-point ',rbinf,zbinf,rbsup,zbsup + write (msg, '("X-point not found in", 2(x,a,"∈[",g0.3,",",g0.3,"]"))') & + 'r', rbinf, rbsup, 'z', zbinf, zbsup + call log_info(msg, mod='equilibrium', proc='set_eqspl') end if - print*,' ' ! Save Bt value on axis (required in flux_average and used in Jcd def) ! and vacuum value B0 at ref. radius data%rvac (used in Jcd_astra def) @@ -583,8 +611,9 @@ contains btaxis = btaxis/rmaxis btrcen = fpolas/data%rvac rcen = data%rvac - print '(a,f8.4)', 'BT_centr=', btrcen - print '(a,f8.4)', 'BT_axis =', btaxis + + write (msg, '(2(a,g0.3))') 'Bt_center=', btrcen, ' Bt_axis=', btaxis + call log_info(msg, mod='equilibrium', proc='set_eqspl') ! Compute rho_pol/rho_tor mapping based on input q profile call setqphi_num(data%psinr,abs(data%qpsi),abs(psia),rhotn) @@ -1072,18 +1101,24 @@ contains subroutine psi_raxis(psin,r1,r2) use const_and_precisions, only : wp_ - use gray_params, only : iequil - use dierckx, only : profil,sproota + use gray_params, only : iequil + use dierckx, only : profil, sproota + use logger, only : log_error + implicit none -! local constants - integer, parameter :: mest=4 -! arguments + + ! subroutine arguments real(wp_) :: psin,r1,r2 -! local variables + + ! local constants + integer, parameter :: mest=4 + + ! local variables integer :: iopt,ier,m real(wp_) :: zc,val real(wp_), dimension(mest) :: zeroc real(wp_), dimension(nsr) :: czc + character(64) :: msg if (iequil < 2) then val=frhotor(sqrt(psin)) @@ -1093,7 +1128,11 @@ contains iopt=1 zc=zmaxis call profil(iopt,tr,nsr,tz,nsz,cceq,kspl,kspl,zc,nsr,czc,ier) - if(ier.gt.0) print*,' profil =',ier + if (ier > 0) then + write (msg, '("profil failed with error ",g0)') ier + call log_error(msg, mod='equilibrium', proc='psi_raxis') + end if + val=psin*psiant+psinop call sproota(val,tr,nsr,czc,zeroc,mest,m,ier) r1=zeroc(1) @@ -1118,44 +1157,57 @@ contains subroutine points_ox(rz,zz,rf,zf,psinvf,info) use const_and_precisions, only : comp_eps - use minpack, only : hybrj1 + use minpack, only : hybrj1 + use logger, only : log_error, log_debug + implicit none -! local constants + + ! local constants integer, parameter :: n=2,ldfjac=n,lwa=(n*(n+13))/2 -! arguments + + ! arguments real(wp_), intent(in) :: rz,zz real(wp_), intent(out) :: rf,zf,psinvf integer, intent(out) :: info -! local variables + + ! local variables real(wp_) :: tol real(wp_), dimension(n) :: xvec,fvec real(wp_), dimension(lwa) :: wa real(wp_), dimension(ldfjac,n) :: fjac + character(256) :: msg xvec(1)=rz xvec(2)=zz tol = sqrt(comp_eps) call hybrj1(fcnox,n,xvec,fvec,fjac,ldfjac,tol,info,wa,lwa) - if(info.gt.1) then - print'(a,i2,a,2f8.4)',' info subr points_ox =',info, & - ' O/X coord.',xvec - end if - rf=xvec(1) - zf=xvec(2) - call equinum_psi(rf,zf,psinvf) + if(info /= 1) then + write (msg, '("O,X coordinates:",2(x,", ",g0.3))') xvec + call log_debug(msg, mod='equilibrium', proc='points_ox') + write (msg, '("hybrj1 failed with error ",g0)') info + call log_error(msg, mod='equilibrium', proc='points_ox') + end if + rf=xvec(1) + zf=xvec(2) + call equinum_psi(rf,zf,psinvf) end subroutine points_ox subroutine fcnox(n,x,fvec,fjac,ldfjac,iflag) + use logger, only : log_error + implicit none -! arguments + + ! subroutine arguments integer, intent(in) :: n,iflag,ldfjac real(wp_), dimension(n), intent(in) :: x real(wp_), dimension(n), intent(inout) :: fvec real(wp_), dimension(ldfjac,n), intent(inout) :: fjac -! local variables + + ! local variables real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz + character(64) :: msg select case(iflag) case(1) @@ -1170,7 +1222,8 @@ contains fjac(2,1) = ddpsidrz/psia fjac(2,2) = ddpsidzz/psia case default - print*,'iflag undefined' + write (msg, '("invalid iflag: ",g0)') + call log_error(msg, mod='equilibrium', proc='fcnox') end select end subroutine fcnox @@ -1178,15 +1231,21 @@ contains subroutine points_tgo(rz,zz,rf,zf,psin0,info) use const_and_precisions, only : comp_eps - use minpack, only : hybrj1mv + use minpack, only : hybrj1mv + use logger, only : log_error, log_debug + implicit none -! local constants + + ! local constants integer, parameter :: n=2,ldfjac=n,lwa=(n*(n+13))/2 -! arguments + + ! arguments real(wp_), intent(in) :: rz,zz,psin0 real(wp_), intent(out) :: rf,zf integer, intent(out) :: info -! local variables + character(256) :: msg + + ! local variables real(wp_) :: tol real(wp_), dimension(n) :: xvec,fvec,f0 real(wp_), dimension(lwa) :: wa @@ -1198,9 +1257,11 @@ contains f0(2)=0.0_wp_ tol = sqrt(comp_eps) call hybrj1mv(fcntgo,n,xvec,f0,fvec,fjac,ldfjac,tol,info,wa,lwa) - if(info.gt.1) then - print'(a,i2,a,5f8.4)',' info subr points_tgo =',info, & - ' R,z coord.',xvec,rz,zz,psin0 + if(info /= 1) then + write (msg, '("R,z coordinates:",5(x,g0.3))') xvec, rz, zz, psin0 + call log_debug(msg, mod='equilibrium', proc='points_tgo') + write (msg, '("hybrj1mv failed with error ",g0)') info + call log_error(msg, mod='equilibrium', proc='points_tgo') end if rf=xvec(1) zf=xvec(2) @@ -1210,14 +1271,19 @@ contains subroutine fcntgo(n,x,f0,fvec,fjac,ldfjac,iflag) use const_and_precisions, only : wp_ + use logger, only : log_error + implicit none -! arguments + + ! subroutine arguments integer, intent(in) :: n,ldfjac,iflag real(wp_), dimension(n), intent(in) :: x,f0 real(wp_), dimension(n), intent(inout) :: fvec real(wp_), dimension(ldfjac,n), intent(inout) :: fjac -! internal variables + + ! local variables real(wp_) :: psinv,dpsidr,dpsidz,ddpsidrr,ddpsidrz + character(64) :: msg select case(iflag) case(1) @@ -1232,7 +1298,8 @@ contains fjac(2,1) = ddpsidrr/psia fjac(2,2) = ddpsidrz/psia case default - print*,'iflag undefined' + write (msg, '("invalid iflag: ",g0)') + call log_error(msg, mod='equilibrium', proc='fcntgo') end select end subroutine fcntgo diff --git a/src/gray_cli.f90 b/src/gray_cli.f90 index 93daa65..7a5304d 100644 --- a/src/gray_cli.f90 +++ b/src/gray_cli.f90 @@ -10,13 +10,13 @@ module gray_cli ! 2. the print_cli_options() subroutine type cli_options ! Switches - logical :: verbose logical :: quiet ! Files character(len=:), allocatable :: output_dir character(len=:), allocatable :: params_file character(len=:), allocatable :: sum_filelist ! others + integer :: verbose integer, allocatable :: units(:) end type @@ -40,8 +40,9 @@ contains print '(a)', 'Options:' print '(a)', ' -h, --help display this help and exit' print '(a)', ' -V, --version display version information and exit' - print '(a)', ' -v, --verbose print additional information messages' - print '(a)', ' -q, --quiet suppress all messages on standard output' + print '(a)', ' -v, --verbose print more information messages;' + print '(a)', ' repeating -v increase the log verbosity' + print '(a)', ' -q, --quiet suppress all non-critical messages' print '(a)', ' -o, --output-dir DIR specify where to write the output files' print '(a)', ' (default: current directory)' print '(a)', ' -p, --params-file FILE set the parameters file' @@ -77,21 +78,22 @@ contains type(cli_options), intent(in) :: opts print '(a)' , 'switches:' - print '(a,l)' , ' - verbose: ' , opts%verbose print '(a,l)' , ' - quiet: ' , opts%quiet print '(a)' , 'files:' print '(a,a)' , ' output-dir: ' , opts%output_dir print '(a,a)' , ' param-file: ' , opts%params_file print '(a,a)' , ' sum: ' , opts%sum_filelist print '(a)' , 'others:' - print '(a,20i3)' , ' - units: ' , opts%units + print '(a,20i3)' , ' - units: ' , opts%units + print '(a,l)' , ' - verbose: ' , opts%verbose end subroutine subroutine parse_cli_options(opts) ! Parse the CLI arguments and initialise the options - use units, only : ucenr, usumm + use units, only : ucenr, usumm + use logger, only : WARNING implicit none @@ -105,7 +107,7 @@ contains integer :: error, commas ! Default option values - opts%verbose = .false. + opts%verbose = WARNING opts%quiet = .false. opts%params_file = 'gray_params.data' opts%units = [ucenr, usumm] @@ -131,7 +133,7 @@ contains call exit(0) case ('-v', '--verbose') - opts%verbose = .true. + opts%verbose = opts%verbose + 1 case ('-q', '--quiet') opts%quiet = .true. diff --git a/src/gray_core.f90 b/src/gray_core.f90 index 81c1dcf..12ec932 100644 --- a/src/gray_core.f90 +++ b/src/gray_core.f90 @@ -24,6 +24,7 @@ contains use multipass, only : alloc_multipass, dealloc_multipass, initbeam, & initmultipass, turnoffray, plasma_in, plasma_out, wall_out use units, only : ucenr + use logger, only : log_info, log_debug implicit none @@ -41,7 +42,7 @@ contains ! local variables real(wp_), parameter :: taucr = 12._wp_, etaucr = exp(-taucr) character, dimension(2), parameter :: mode=(/'O','X'/) - + real(wp_) :: sox,ak0,bres,xgcn,xg,yg,rrm,zzm,alpha,didp,anpl,anpr,anprim,anprre real(wp_) :: chipol,psipol,btot,psinv,dens,tekev,dersdst,derdnm real(wp_) :: tau,pow,dids,ddr,ddi,taumn,taumx @@ -81,6 +82,8 @@ contains ! parameters log in file headers character(len=headw), dimension(headl) :: strheader + ! buffer for formatting log messages + character(256) :: msg ! ======== set environment BEGIN ======== ! Number of limiter contourn points @@ -129,13 +132,16 @@ contains call print_parameters(params, strheader) call print_headers(strheader) - ! print ψ surface for q=1.5 and q=2 on file and psi,rhot,rhop on stdout + ! print ψ surface for q=1.5 and q=2 on file and log psi,rhot,rhop call print_surfq([1.5_wp_, 2.0_wp_]) ! print initial position - print *, '' - print '(a,2f8.3)', 'alpha0, beta0 = ', params%antenna%alpha, params%antenna%beta - print '(a,4f8.3)', 'x00, y00, z00 = ', params%antenna%pos + write (msg, '("initial position:",3(x,g0.3))') params%antenna%pos + call log_info(msg, mod='gray_core', proc='gray_main') + + write (msg, '("initial direction:",2(x,a,"=",g0.2))') & + 'α', params%antenna%alpha, 'β', params%antenna%beta + call log_info(msg, mod='gray_core', proc='gray_main') ! print Btot=Bres ! print ne, Te, q, Jphi versus psi, rhop, rhot @@ -169,8 +175,11 @@ contains nbeam_pass=1 ! max n of beam per pass index_rt=0 ! global beam index: 1,O 2,X 1st pass ! | | | | - do ip=1,ipass ! 3,O 4,X 5,O 6,X 2nd pass - + call log_debug('pass loop start', mod='gray_core', proc='gray_main') ! 3,O 4,X 5,O 6,X 2nd pass + do ip=1,ipass + write (msg, '("pass: ",g0)') ip + call log_info(msg, mod='gray_core', proc='gray_main') + pabs_pass = zero icd_pass = zero istop_pass = 0 ! stop flag for current pass @@ -184,6 +193,7 @@ contains end if ! =========== beam loop BEGIN =========== + call log_debug('beam loop start', mod='gray_core', proc='gray_main') do ib=1,nbeam_pass sox = -sox ! invert mode @@ -194,9 +204,13 @@ contains call initbeam(index_rt,iroff,iboff,iwait,stv,jphi_beam, & pins_beam,currins_beam,dpdv_beam,jcd_beam) - + + write(msg, '(" beam: ",g0," (",a1," mode)")') index_rt, mode(iox) + call log_info(msg, mod='gray_core', proc='gray_main') + if(iboff) then ! no propagation for current beam istop_pass = istop_pass +1 ! * +1 non propagating beam + call log_info(" beam is off", mod='gray_core', proc='gray_main') cycle end if @@ -248,8 +262,9 @@ contains if(nray>1 .and. all(.not.iwait)) call print_projxyzt(stv,yw,0) ! iproj=0 ==> nfilp=8 ! ======= propagation loop BEGIN ======= + call log_debug(' propagation loop start', mod='gray_core', proc='gray_main') + do i=1,nstep - ! advance one step with "frozen" grad(S_I) do jk=1,nray if(iwait(jk)) cycle ! jk ray is waiting for next pass @@ -309,9 +324,9 @@ contains cpls(jk,index_rt) = cpl(iox) if(jk.eq.1) then - write(*,*) - write(*,'("1st pass coupling (central ray, ",a1,"-mode)",f9.4)') & - mode(iox),cpl(iox) + write (msg,'(" 1st pass - central ray (",a1,"-mode) c=",g0.4)') & + mode(iox), cpl(iox) + call log_info(msg, mod='gray_core', proc='gray_main') psipv(index_rt) = psipol ! + polarization angles at plasma boundary for central ray chipv(index_rt) = chipol end if @@ -481,6 +496,7 @@ contains exit end if end do + call log_debug(' propagation loop end', mod='gray_core', proc='gray_main') ! ======== propagation loop END ======== ! print all ray positions in local reference system @@ -513,21 +529,31 @@ contains cpl_beam1 = zero cpl_beam2 = zero end if - - ! print final results for pass on screen - write(*,*) - write(*,'("End of propagation for beam ",i5," (pass ",i3,", ",a1," mode)")') & - index_rt,ip,mode(iox) - write(*,'(a,f9.4)') 'final step (s, ct, Sr) = ',stv(1) - write(*,'(a,2e12.5)') 'taumn, taumx = ', taumn,taumx - write(*,'(a,f9.4)') 'Pabs_tot (MW) = ',pabs_beam - write(*,'(a,f9.4)') 'I_tot (kA) = ',icd_beam*1.0e3_wp_ - if(ip.lt.ipass) then - write(*,'(a,2(f9.4,1x))') 'Coupling (average, O/X):',cpl_beam1,cpl_beam2 ! average coupling for next O/X beams (=0 if no ray re-entered plasma) - if(iop(1).gt.2) write(*,'(a,2(f9.4,1x))') 'Coupling (ctr ray, O/X):', & - cpl_cbeam1,cpl_cbeam2 ! central ray coupling for next O/X beams + + ! print final results for pass on screen + call log_info(' partial results:', mod='gray_core', proc='gray_main') + write(msg, '(3x,a,g0.4)') 'final step: (s, ct, Sr)=' ,stv(1) + call log_info(msg, mod='gray_core', proc='gray_main') + + write(msg, '(3x,a,2(x,a,"=",g0.4))') 'optical depth:', 'τ_min', taumn, 'τ_max', taumx + call log_info(msg, mod='gray_core', proc='gray_main') + + write(msg, '(3x,a,g0.3," MW")') 'absoption: P=', pabs_beam + call log_info(msg, mod='gray_core', proc='gray_main') + + write(msg, '(3x,a,g0.3," MW")') 'current drive: I=', icd_beam*1.0e3_wp_ + call log_info(msg, mod='gray_core', proc='gray_main') + + if(ip < ipass) then + write (msg,'(3x,a,(g0.4,", ",g0.4))') & ! average coupling for next O/X beams (=0 if no ray re-entered plasma) + 'next couplings [O,X mode]: c=', cpl_beam1, cpl_beam2 + call log_info(msg, mod='gray_core', proc='gray_main') + if(iop(1) > 2) then + write(msg, '(3x,a,(g0.4,", ",g0.4))') & + 'coupling [ctr ray, O/X]:', cpl_cbeam1, cpl_cbeam2 ! central ray coupling for next O/X beams + end if end if - + write(ucenr,*) '' call print_pec(rhop_tab,rhot_tab,jphi_beam,jcd_beam,dpdv_beam,currins_beam, & @@ -544,31 +570,30 @@ contains ! ============ post-proc END ============ end do + call log_debug('beam loop end', mod='gray_core', proc='gray_main') ! ============ beam loop END ============ ! ======= cumulative prints BEGIN ======= results%pabs = results%pabs + sum(pabs_pass) ! *final results (O+X) [gray_main output] results%icd = results%icd + sum(icd_pass) - - ! print final results for pass on screen - write(*,*) - write(*,'("# End of pass ",i3)') ip - write(*,'(a,f9.4,f9.4)') '# Pabs_tot (MW) [O,X mode] = ',pabs_pass(1),pabs_pass(2) - write(*,'(a,f9.4,f9.4)') '# I_tot (kA) [O,X mode] = ', & - icd_pass(1)*1.0e3_wp_,icd_pass(2)*1.0e3_wp_ + + ! print final results for pass on screen + call log_info(' comulative results:', mod='gray_core', proc='gray_main') + + write(msg, '(" absoption [O,X mode] P=",g0.4,", ",g0.4," MW")') & + pabs_pass(1), pabs_pass(2) + call log_info(msg, mod='gray_core', proc='gray_main') + + write(msg, '(" current drive [O,X mode] I=",g0.4,", ",g0.4," kA")') & + icd_pass(1)*1.0e3_wp_, icd_pass(2)*1.0e3_wp_ + call log_info(msg, mod='gray_core', proc='gray_main') ! ======== cumulative prints END ======== if(istop_pass == nbeam_pass) exit ! no active beams - end do + call log_debug('pass loop end', mod='gray_core', proc='gray_main') ! ============ main loop END ============ -! print final results on screen - write(*,*) - write(*,'(a)') '## Final results:' - write(*,'(a,f9.4)') '## Pabs_tot (MW) = ', results%pabs - write(*,'(a,f9.4)') '## I_tot (kA) = ', results%icd*1.0e3_wp_ - ! ========== free memory BEGIN ========== call dealloc_surfvec call dealloc_beam(yw,ypw,xc,du1,gri,ggri,psjki,ppabs,ccci,tau0, & @@ -2039,30 +2064,36 @@ bb: do end subroutine print_maps - - subroutine print_surfq(qval) - use equilibrium, only : psinr,nq,fq,frhotor,rmaxis,zmaxis, & - zbsup,zbinf - use magsurf_data, only : contours_psi,npoints,print_contour - use utils, only : locate, intlin + use equilibrium, only : psinr, nq, fq, frhotor, & + rmaxis, zmaxis, zbsup, zbinf + use magsurf_data, only : contours_psi, npoints, print_contour + use utils, only : locate, intlin + use logger, only : log_info + implicit none -! arguments + + ! subroutine arguments real(wp_), dimension(:), intent(in) :: qval -! local variables + + ! local variables integer :: i1,i real(wp_) :: rup,zup,rlw,zlw,rhot,psival real(wp_), dimension(npoints) :: rcn,zcn real(wp_), dimension(nq) :: qpsi + character(256) :: msg ! for log messages formatting -! build q profile on psin grid + ! build q profile on psin grid do i=1,nq - qpsi(i) = fq(psinr(i)) + qpsi(i) = fq(psinr(i)) end do -! locate psi surface for q=qval - print* + + ! locate ψ surface for q=qval + call log_info('constant ψ surfaces for:', & + mod='gray_core', proc='print_surfq') do i=1,size(qval) - call locate(abs(qpsi),nq,qval(i),i1) !!!! check for non monotonous q profile + ! FIXME: check for non monotonous q profile + call locate(abs(qpsi),nq,qval(i),i1) if (i1>0.and.i1 stop beam iboff = .true. - write(*,*) - write(*,'("Beam ",i5," inactive")') i else if(.not.all(.not.iwait)) then ! only some rays active - write(*,*) - write(*,'("WARNING: not all rays in beam ",i5," are active")') i + write (msg,'(" beam ",g0,": some rays are active!")') i + call log_warning(msg, mod='multipass', proc='initbeam') end if stv = zero ! starting step