diff --git a/src/beams.f90 b/src/beams.f90 index 34acf2a..290ac0e 100644 --- a/src/beams.f90 +++ b/src/beams.f90 @@ -5,7 +5,7 @@ module beams contains - subroutine read_beam0(params, err, unit) + subroutine read_beam0(params, err) ! Reads the wave launcher parameters for the simple case ! where w(z) and 1/R(z) are fixed. ! @@ -30,21 +30,17 @@ 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 ! subroutine arguments - type(antenna_parameters), intent(inout) :: params - integer, intent(out) :: err - integer, intent(in), optional :: unit + type(antenna_parameters), intent(inout) :: params + integer, intent(out) :: err ! local variables integer :: u real(wp_) :: k0, w0(2), d0(2), z_R(2), phi - u = get_free_unit(unit) - - open(unit=u, file=params%filenm, status='old', action='read', iostat=err) + open(newunit=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") @@ -79,7 +75,7 @@ contains end subroutine read_beam0 - subroutine read_beam1(params, err, unit) + subroutine read_beam1(params, err) ! Reads the wave launcher parameters for the case ! where w(z, α) and 1/R(z, α) depend on the launcher angle α. ! @@ -100,13 +96,12 @@ contains use gray_params, only : antenna_parameters use splines, only : spline_simple - use utils, only : get_free_unit,locate + use utils, only : locate use logger, only : log_error ! subroutine arguments type(antenna_parameters), intent(inout) :: params integer, intent(out) :: err - integer, optional, intent(in) :: unit ! local variables integer :: u, nisteer, i, k, ii @@ -118,9 +113,7 @@ contains rci1, rci2, phi1, phi2, & x0, y0, z0 - u = get_free_unit(unit) - - open(unit=u, file=params%filenm, status='old', action='read', iostat=err) + open(newunit=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") @@ -164,7 +157,7 @@ contains call z0%init(alphastv, z00v, nisteer) if((params%alpha > alphastv(1)) .and. (params%alpha < alphastv(nisteer))) then - call locate(alphastv, nisteer, params%alpha , k) + call locate(alphastv, params%alpha, k) dal = params%alpha - alphastv(k) params%beta = beta%raw_eval(k, dal) params%pos(1) = x0%raw_eval(k, dal) @@ -211,7 +204,7 @@ contains end subroutine read_beam1 - subroutine read_beam2(params, beamid, err, unit) + subroutine read_beam2(params, beamid, err) ! Reads the wave launcher parameters for the general case ! where w(z, α, β) and 1/R(z, α, β) depend on the launcher angles α, β. ! @@ -232,7 +225,8 @@ contains ! 5. Each line stores one record with the same fields as in the 1D format. use gray_params, only : antenna_parameters - use utils, only : get_free_unit, intlin, locate, inside + use utils, only : linear_interp, locate + use types, only : contour use dierckx, only : curfit, splev, surfit, bispev use logger, only : log_error @@ -240,7 +234,6 @@ contains type(antenna_parameters), intent(inout) :: params integer, intent(in) :: beamid integer, intent(out), optional :: err - integer, intent(in), optional :: unit ! local variables character(len=20) :: beamname @@ -259,17 +252,14 @@ contains ycoord, wrk, txwaist1, tywaist1, txwaist2, tywaist2, & txrci1, tyrci1, txrci2, tyrci2, txphi1, typhi1, txphi2, typhi2, & txx0, tyx0, txy0, tyy0, txz0, tyz0, txycoord, cycoord, cwaist1, & - cwaist2, crci1, crci2, cphi1,cphi2, cx0, cy0, cz0, w, wrk2, & - xpolyg, ypolyg, xpolygA, ypolygA, xpolygB, ypolygB, xpolygC, & - ypolygC, xpolygD, ypolygD, xoutA, youtA, xoutB, youtB, xoutC, youtC + cwaist2, crci1, crci2, cphi1,cphi2, cx0, cy0, cz0, w, wrk2 + type(contour) :: polyg, polygA, polygB, polygC, polygD, outA, outB, outC real(wp_), DIMENSION(4) :: xvert, yvert real(wp_), dimension(1) :: fi integer, parameter :: kspl=1 real(wp_), parameter :: sspl=0.01_wp_ - u = get_free_unit(unit) - - open(unit=u, file=params%filenm, status='old', action='read', iostat=err) + open(newunit=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") @@ -470,7 +460,7 @@ contains cy0(nxest*nyest), txy0(nmax), tyy0(nmax), & cz0(nxest*nyest), txz0(nmax), tyz0(nmax), & wrk(lwrk), wrk2(lwrk2), iwrk(kwrk), & - xpolyg(npolyg), ypolyg(npolyg), w(nisteer)) + polyg%R(npolyg), polyg%z(npolyg), w(nisteer)) ! w = 1.d0 ! @@ -511,26 +501,26 @@ contains minx,maxx,miny,maxy,kx,ky,sspl,nxest,nyest,nmax,eps, & nxz0,txz0,nyz0,tyz0,cz0,fp,wrk,lwrk,wrk2,lwrk2,iwrk,kwrk,ier) ! data range polygon - xpolyg(1:nxcoord) = xcoord(1:nxcoord) - ypolyg(1:nxcoord) = ycoord(1:nxcoord) + polyg%R(1:nxcoord) = xcoord(1:nxcoord) + polyg%z(1:nxcoord) = ycoord(1:nxcoord) ! ! c==================================================================================== do i=1,nycoord-2 - xpolyg(nxcoord+i) = xcoord((i+1)*nxcoord) - xpolyg(2*nxcoord+nycoord-2+i) = xcoord((nycoord-i-1)*nxcoord+1) - ypolyg(nxcoord+i) = ycoord((i+1)*nxcoord) - ypolyg(2*nxcoord+nycoord-2+i) = ycoord((nycoord-i-1)*nxcoord+1) + polyg%R(nxcoord+i) = xcoord((i+1)*nxcoord) + polyg%R(2*nxcoord+nycoord-2+i) = xcoord((nycoord-i-1)*nxcoord+1) + polyg%z(nxcoord+i) = ycoord((i+1)*nxcoord) + polyg%z(2*nxcoord+nycoord-2+i) = ycoord((nycoord-i-1)*nxcoord+1) end do ! c==================================================================================== do i=1,nxcoord - xpolyg(nxcoord+nycoord-2+i) = xcoord(nxcoord*nycoord-i+1) - ypolyg(nxcoord+nycoord-2+i) = ycoord(nxcoord*nycoord-i+1) + polyg%R(nxcoord+nycoord-2+i) = xcoord(nxcoord*nycoord-i+1) + polyg%z(nxcoord+nycoord-2+i) = ycoord(nxcoord*nycoord-i+1) end do ! c==================================================================================== ! ! check if (xcoord0, ycoord0) is out of table range ! incheck = 1 inside / 0 outside - if(inside(xpolyg, ypolyg, xcoord0, ycoord0)) incheck = 1 + if (polyg%contains(xcoord0, ycoord0)) incheck = 1 end if deallocate(wrk,iwrk) !####################################################################################### @@ -582,50 +572,52 @@ contains else ! c==================================================================================== if(incheck.eq.0) then - allocate(xpolygA(nxcoord), ypolygA(nxcoord), xpolygC(nxcoord), & - ypolygC(nxcoord), xpolygB(nycoord), ypolygB(nycoord), & - xpolygD(nycoord), ypolygD(nycoord), & - xoutA(nxcoord+3), youtA(nxcoord+3), xoutB(nycoord+3), & - youtB(nycoord+3), xoutC(nxcoord+3), youtC(nxcoord+3)) + allocate(polygA%R(nxcoord), polygA%z(nxcoord)) + allocate(polygC%R(nxcoord), polygC%z(nxcoord)) + allocate(polygB%R(nycoord), polygB%z(nycoord)) + allocate(polygD%R(nycoord), polygD%z(nycoord)) + allocate(outA%R(nxcoord+3), outA%z(nxcoord+3)) + allocate(outB%R(nycoord+3), outB%z(nycoord+3)) + allocate(outC%R(nxcoord+3), outC%z(nxcoord+3)) ! coordinates of vertices v1,v2,v3,v4 - xvert(1) = xpolyg(1) - xvert(2) = xpolyg(nxcoord) - xvert(3) = xpolyg(nxcoord+nycoord-1) - xvert(4) = xpolyg(2*nxcoord+nycoord-2) - yvert(1) = ypolyg(1) - yvert(2) = ypolyg(nxcoord) - yvert(3) = ypolyg(nxcoord+nycoord-1) - yvert(4) = ypolyg(2*nxcoord+nycoord-2) + xvert(1) = polyg%R(1) + xvert(2) = polyg%R(nxcoord) + xvert(3) = polyg%R(nxcoord+nycoord-1) + xvert(4) = polyg%R(2*nxcoord+nycoord-2) + yvert(1) = polyg%z(1) + yvert(2) = polyg%z(nxcoord) + yvert(3) = polyg%z(nxcoord+nycoord-1) + yvert(4) = polyg%z(2*nxcoord+nycoord-2) ! coordinates of side A,B,C,D - xpolygA = xpolyg(1:nxcoord) - ypolygA = ypolyg(1:nxcoord) - xpolygB = xpolyg(nxcoord:nxcoord+nycoord-1) - ypolygB = ypolyg(nxcoord:nxcoord+nycoord-1) - xpolygC = xpolyg(nxcoord+nycoord-1:2*nxcoord+nycoord-2) - ypolygC = ypolyg(nxcoord+nycoord-1:2*nxcoord+nycoord-2) - xpolygD(1:nycoord-1) = xpolyg(2*nxcoord+nycoord-2:npolyg) - xpolygD(nycoord) = xpolyg(1) - ypolygD(1:nycoord-1) = ypolyg(2*nxcoord+nycoord-2:npolyg) - ypolygD(nycoord) = ypolyg(1) + polygA%R = polyg%R(1:nxcoord) + polygA%z = polyg%z(1:nxcoord) + polygB%R = polyg%R(nxcoord:nxcoord+nycoord-1) + polygB%z = polyg%z(nxcoord:nxcoord+nycoord-1) + polygC%R = polyg%R(nxcoord+nycoord-1:2*nxcoord+nycoord-2) + polygC%z = polyg%z(nxcoord+nycoord-1:2*nxcoord+nycoord-2) + polygD%R(1:nycoord-1) = polyg%R(2*nxcoord+nycoord-2:npolyg) + polygD%R(nycoord) = polyg%R(1) + polygD%z(1:nycoord-1) = polyg%z(2*nxcoord+nycoord-2:npolyg) + polygD%z(nycoord) = polyg%z(1) ! contour of outside regions A (1,2), B(3,4), C(5,6) - xoutA = huge(one) - xoutA(1:nxcoord) = xpolygA - xoutA(nxcoord+3) = xvert(1) - youtA = -huge(one) - youtA(1:nxcoord) = ypolygA - youtA(nxcoord+1) = yvert(2) - xoutB = huge(one) - xoutB(1:nycoord) = xpolygB - xoutB(nycoord+1) = xvert(3) - youtB = huge(one) - youtB(1:nycoord) = ypolygB - youtB(nycoord+3) = yvert(2) - xoutC = -huge(one) - xoutC(1:nxcoord) = xpolygC - xoutC(nxcoord+3) = xvert(3) - youtC = huge(one) - youtC(1:nxcoord) = ypolygC - youtC(nxcoord+1) = yvert(4) + outA%R = huge(one) + outA%R(1:nxcoord) = polygA%R + outA%R(nxcoord+3) = xvert(1) + outA%z = -huge(one) + outA%z(1:nxcoord) = polygA%z + outA%z(nxcoord+1) = yvert(2) + outB%R = huge(one) + outB%R(1:nycoord) = polygB%R + outB%R(nycoord+1) = xvert(3) + outB%z = huge(one) + outB%z(1:nycoord) = polygB%z + outB%z(nycoord+3) = yvert(2) + outC%R = -huge(one) + outC%R(1:nxcoord) = polygC%R + outC%R(nxcoord+3) = xvert(3) + outC%z = huge(one) + outC%z(1:nxcoord) = polygC%z + outC%z(nxcoord+1) = yvert(4) ! c---------------------------------------------------------------------------------- ! search for position of xcoord0, ycoord0 with respect to (alpha,beta) data grid ! @@ -639,17 +631,17 @@ contains ! v1 A v2 ! (8) | (1) | (2) ! - if(inside(xoutA, youtA, xcoord0, ycoord0)) then + if (outA%contains(xcoord0, ycoord0)) then in = 1 if(xcoord0.GT.xvert(2)) then in = 2 end if - else if(inside(xoutB, youtB, xcoord0, ycoord0)) then + else if (outB%contains(xcoord0, ycoord0)) then in = 3 if(ycoord0.GT.yvert(3)) then in = 4 end if - else if(inside(xoutC, youtC, xcoord0, ycoord0)) then + else if (outC%contains(xcoord0, ycoord0)) then in = 5 if(xcoord0.LT.xvert(4)) then in = 6 @@ -678,9 +670,9 @@ contains CASE (1) ! params%beta outside table range ! locate position of xcoord0 with respect to x coordinates of side A - call locate(xpolygA,nxcoord,xcoord0,ii) + call locate(polygA%R, xcoord0, ii) ! find corresponding y value on side A for xcoord position - call intlin(xpolygA(ii),ypolygA(ii),xpolygA(ii+1),ypolygA(ii+1),xcoord0,ycoord0) + ycoord0 = linear_interp(polygA%R(ii), polygA%z(ii), polygA%R(ii+1), polygA%z(ii+1), xcoord0) incheck = 1 CASE (2) ! params%alpha and params%beta outside table range @@ -690,8 +682,8 @@ contains ii = nxcoord !indice per assegnare valori waist, rci, phi CASE (3) ! params%alpha outside table range - call locate(ypolygB,nycoord,ycoord0,ii) - call intlin(ypolygB(ii),xpolygB(ii),ypolygB(ii+1),xpolygB(ii+1),ycoord0,xcoord0) + call locate(polygB%z, ycoord0, ii) + xcoord0 = linear_interp(polygB%z(ii), polygB%R(ii), polygB%z(ii+1), polygB%R(ii+1), ycoord0) incheck = 1 CASE (4) ! params%alpha and params%beta outside table range @@ -700,8 +692,8 @@ contains ii = nxcoord+nycoord-1 CASE (5) ! params%beta outside table range - call locate(xpolygC,nxcoord,xcoord0,ii) - call intlin(xpolygC(ii+1),ypolygC(ii+1),xpolygC(ii),ypolygC(ii),xcoord0,ycoord0) + call locate(polygC%R, xcoord0, ii) + ycoord0 = linear_interp(polygC%R(ii+1), polygC%z(ii+1), polygC%R(ii), polygC%z(ii), xcoord0) incheck = 1 CASE (6) ! params%alpha and params%beta outside table range @@ -710,8 +702,8 @@ contains ii = 2*nxcoord+nycoord-2 CASE (7) ! params%alpha outside table range - call locate(ypolygD,nycoord,ycoord0,ii) - call intlin(ypolygD(ii),xpolygD(ii),ypolygD(ii+1),xpolygD(ii+1),ycoord0,xcoord0) + call locate(polygD%z, ycoord0, ii) + xcoord0 = linear_interp(polygD%z(ii), polygD%R(ii), polygD%z(ii+1), polygD%R(ii+1), ycoord0) incheck = 1 CASE (8) ! params%alpha and params%beta outside table range @@ -721,8 +713,8 @@ contains END SELECT ! c---------------------------------------------------------------------------------- ! - deallocate(xpolygA, ypolygA, xpolygC, ypolygC, xpolygB, ypolygB, & - xpolygD, ypolygD, xoutA, youtA, xoutB, youtB, xoutC, youtC) + deallocate(polygA%R, polygA%z, polygC%R, polygC%z, polygB%R, polygB%z, & + polygD%R, polygD%z, outA%R, outA%z, outB%R, outB%z, outC%R, outC%z) end if ! c==================================================================================== ! @@ -785,7 +777,7 @@ contains crci1, txrci1, tyrci1, crci2, txrci2, tyrci2, & cphi1, txphi1, typhi1, cphi2, txphi2, typhi2, & cx0, txx0, tyx0, cy0, txy0, tyy0, cz0, txz0, tyz0, & - wrk2, xpolyg, ypolyg, w) + wrk2, polyg%R, polyg%z, w) end if ! !####################################################################################### diff --git a/src/gray_core.f90 b/src/gray_core.f90 index b12d988..a13c574 100644 --- a/src/gray_core.f90 +++ b/src/gray_core.f90 @@ -23,7 +23,7 @@ contains use magsurf_data, only : compute_flux_averages, dealloc_surfvec use pec, only : pec_init, spec, postproc_profiles, dealloc_pec, & rhop_tab, rhot_tab - use utils, only : vmaxmin, inside + use utils, only : vmaxmin use multipass, only : initbeam, initmultipass, turnoffray, & plasma_in, plasma_out, wall_out use logger, only : log_info, log_debug @@ -508,7 +508,7 @@ contains end if ! test whether further trajectory integration is unnecessary - call vmaxmin(tau1+tau0+lgcpl1,params%raytracing%nray,taumn,taumx) ! test on tau + coupling + call vmaxmin(tau1+tau0+lgcpl1, taumn, taumx) ! test on tau + coupling if(is_critical(error)) then ! stop propagation for current beam istop_pass = istop_pass +1 ! * +1 non propagating beam @@ -531,7 +531,7 @@ contains if(i>params%raytracing%nstep) i=params%raytracing%nstep pabs_beam = sum(ppabs(:,i)) icd_beam = sum(ccci(:,i)) - call vmaxmin(tau0,params%raytracing%nray,taumn,taumx) ! taumn,taumx for print + call vmaxmin(tau0, taumn, taumx) ! taumn,taumx for print if (params%equilibrium%iequil /= EQ_VACUUM) then ! compute power and current density profiles for all rays diff --git a/src/gray_equil.f90 b/src/gray_equil.f90 index 39e7101..b772e12 100644 --- a/src/gray_equil.f90 +++ b/src/gray_equil.f90 @@ -683,7 +683,7 @@ contains use gray_params, only : equilibrium_parameters use gray_params, only : EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL use gray_params, only : X_AT_TOP, X_AT_BOTTOM, X_IS_MISSING - use utils, only : vmaxmin, vmaxmini + use utils, only : vmaxmin use logger, only : log_debug, log_info ! subroutine arguments @@ -844,7 +844,7 @@ contains ! for the search of O/X points nbnd = size(data%boundary%R) if (nbnd > 0) then - call vmaxmini(data%boundary%z, nbnd, zbinf, zbsup, ibinf, ibsup) + call vmaxmin(data%boundary%z, zbinf, zbsup, ibinf, ibsup) rbinf = data%boundary%R(ibinf) rbsup = data%boundary%R(ibsup) else @@ -1387,7 +1387,6 @@ contains ! For a description of the G-EQDSK format, see the GRAY user manual. use const_and_precisions, only : one use gray_params, only : equilibrium_parameters - use utils, only : get_free_unit use logger, only : log_error ! subroutine arguments diff --git a/src/gray_params.f90 b/src/gray_params.f90 index e08543a..574ab4a 100644 --- a/src/gray_params.f90 +++ b/src/gray_params.f90 @@ -447,7 +447,6 @@ contains subroutine read_gray_params(filename, params, err) ! Reads the GRAY parameters from the legacy gray_params.data file - use utils, only : get_free_unit use logger, only : log_error ! subrouting arguments @@ -458,9 +457,7 @@ contains ! local variables integer :: u, temp(3) - u = get_free_unit() - - open(u, file=filename, status='old', action='read', iostat=err) + open(newunit=u, file=filename, status='old', action='read', iostat=err) if (err /= 0) then call log_error('opening gray_params file ('//filename//') failed!', & mod='gray_params', proc='read_gray_params') diff --git a/src/ini_parser.f90 b/src/ini_parser.f90 index 7461a6f..549b32e 100644 --- a/src/ini_parser.f90 +++ b/src/ini_parser.f90 @@ -60,7 +60,6 @@ contains ! ERR_VALUE on invalid values for this property; ! ERR_UNKNOWN on unknown property. ! - use utils, only : get_free_unit ! function argument character(*), intent(in) :: filepath @@ -74,8 +73,7 @@ contains character(len=:), allocatable :: section, name, value ! open the INI file - ini = get_free_unit() - open(unit=ini, file=filepath, action='read', iostat=error) + open(newunit=ini, file=filepath, action='read', iostat=error) if (error /= 0) then write (msg, '("failed to open INI file: ", a)') filepath call log_error(msg, proc='parse_ini', mod='ini_parser') diff --git a/src/pec.f90 b/src/pec.f90 index 8bfc2d4..b5f2080 100644 --- a/src/pec.f90 +++ b/src/pec.f90 @@ -154,7 +154,7 @@ contains subroutine pec_tab(xxi, ypt, yamp, ii, xtab1, wdpdv, wajphiv) ! Power and current projected on ψ grid - mid points - use utils, only : locatex, intlin + use utils, only : locate, linear_interp ! subroutine arguments real(wp_), dimension(:), intent(in) :: xxi, ypt, yamp @@ -222,14 +222,15 @@ contains indi = ind if (idecr == -1) indi = ind - 1 rt1 = xtab1(indi) - call locatex(xxi,iise,iise0,iise,rt1,itb1) - if (itb1 >= iise0 .and. itb1 <= iise) then + call locate(xxi(iise0:iise), rt1, itb1) + if (itb1 >= 1 .and. itb1 <= iise) then + itb1 = itb1 + (iise0 - 1) if(itb1 == iise) then ppa2=ypt(itb1) cci2=yamp(itb1) else - call intlin(xxi(itb1), ypt(itb1),xxi(itb1+1), ypt(itb1+1),rt1,ppa2) - call intlin(xxi(itb1),yamp(itb1),xxi(itb1+1),yamp(itb1+1),rt1,cci2) + ppa2 = linear_interp(xxi(itb1), ypt(itb1), xxi(itb1+1), ypt(itb1+1), rt1) + cci2 = linear_interp(xxi(itb1), yamp(itb1), xxi(itb1+1), yamp(itb1+1), rt1) end if dppa = ppa2 - ppa1 didst = cci2 - cci1 @@ -321,7 +322,7 @@ contains subroutine profwidth(xx,yy,xpk,ypk,dxxe) use const_and_precisions, only : emn1 - use utils, only : locatex, locate, intlin, vmaxmini + use utils, only : locate, linear_interp, vmaxmin ! subroutine arguments real(wp_), intent(in) :: xx(:), yy(:) @@ -332,7 +333,7 @@ contains real(wp_) :: xmn,xmx,ymn,ymx,xpkp,xpkm,yye,rte1,rte2 real(wp_) :: ypkp,ypkm - call vmaxmini(yy,size(yy),ymn,ymx,imn,imx) + call vmaxmin(yy, ymn, ymx, imn, imx) ypk = 0 xmx = xx(imx) xmn = xx(imn) @@ -355,15 +356,16 @@ contains xpk = xpkp ypk = ypkp yye = ypk*emn1 - call locatex(yy,size(yy),1,ipk,yye,ie) + call locate(yy(1:ipk), yye, ie) if(ie > 0 .and. ie < size(yy)) then - call intlin(yy(ie),xx(ie),yy(ie+1),xx(ie+1),yye,rte1) + rte1 = linear_interp(yy(ie), xx(ie), yy(ie+1), xx(ie+1), yye) else rte1 = 0 end if - call locatex(yy,size(yy),ipk,size(yy),yye,ie) + call locate(yy(ipk:size(yy)), yye, ie) if(ie > 0 .and. ie < size(yy)) then - call intlin(yy(ie),xx(ie),yy(ie+1),xx(ie+1),yye,rte2) + ie = ie + (ipk - 1) + rte2 = linear_interp(yy(ie), xx(ie), yy(ie+1), xx(ie+1), yye) else rte2 = 0 end if @@ -373,9 +375,9 @@ contains ypk=yy(2) rte1=0.0_wp_ yye=ypk*emn1 - call locate(yy,size(yy),yye,ie) + call locate(yy, yye, ie) if(ie > 0 .and. ie < size(yy)) then - call intlin(yy(ie),xx(ie),yy(ie+1),xx(ie+1),yye,rte2) + rte2 = linear_interp(yy(ie), xx(ie), yy(ie+1), xx(ie+1), yye) else rte2 = 0 end if diff --git a/src/reflections.f90 b/src/reflections.f90 index 2c63f8f..153f896 100644 --- a/src/reflections.f90 +++ b/src/reflections.f90 @@ -86,7 +86,7 @@ end subroutine inters_linewall subroutine linecone_coord(xv,kv,rs,zs,s,t,n) - use utils, only : bubble + use utils, only : sort_pair real(wp_), intent(in), dimension(3) :: xv,kv real(wp_), intent(in), dimension(2) :: rs,zs real(wp_), dimension(2), intent(out) :: s,t @@ -147,7 +147,7 @@ subroutine linecone_coord(xv,kv,rs,zs,s,t,n) else s(1) = (-bhalf+sqrt(delta))/a s(2) = (-bhalf-sqrt(delta))/a - call bubble(s,2) + call sort_pair(s) t(:) = (kz*s(:)+(z0-zs(1)))/dz n = 2 end if diff --git a/src/splines.f90 b/src/splines.f90 index 308c7a1..2e8ed83 100644 --- a/src/splines.f90 +++ b/src/splines.f90 @@ -121,7 +121,7 @@ contains integer :: i real(wp_) :: dx - call locate(self%data, self%ndata, x, i) + call locate(self%data, x, i) i = min(max(1, i), self%ndata - 1) dx = x - self%data(i) y = self%raw_eval(i, dx) @@ -156,7 +156,7 @@ contains integer :: i real(wp_) :: dx - call locate(self%data, self%ndata, x, i) + call locate(self%data, x, i) i = min(max(1, i), self%ndata - 1) dx = x - self%data(i) y = self%coeffs(i,2) + dx*(2*self%coeffs(i,3) + 3*dx*self%coeffs(i,4)) @@ -691,7 +691,7 @@ contains ! local variables integer :: i - call locate(self%xdata, self%ndata, x, i) + call locate(self%xdata, x, i) i = min(max(1, i), self%ndata - 1) y = self%raw_eval(i, x) end function linear_1d_eval diff --git a/src/types.f90 b/src/types.f90 index 7acf18c..b4f7fc1 100644 --- a/src/types.f90 +++ b/src/types.f90 @@ -190,7 +190,7 @@ contains ! Save the table to a file ! ! Note: this operation consumes the table - use utils, only : get_free_unit, digits + use utils, only : digits ! suborutine arguments class(table), intent(inout) :: self @@ -206,8 +206,6 @@ contains if (self%empty()) return - file = get_free_unit() - ! set default file name if (present(filepath)) then filepath_ = filepath @@ -216,7 +214,7 @@ contains end if ! open output file - open(unit=file, file=filepath_, action='write', iostat=error) + open(newunit=file, file=filepath_, action='write', iostat=error) if (error /= 0) return ! write a custom file header @@ -299,7 +297,7 @@ contains pure function contour_contains(self, R0, z0) result(inside) ! Tests whether the point (`R`, `z`) lies inside the 2D contour - use utils, only : intlinf, locate_unord + use utils, only : linear_interp, locate_unord ! subroutine arguments class(contour), intent(in) :: self @@ -324,8 +322,8 @@ contains ! means that the point is outside the polygon. do i = 1, nsegs ! coordinate of the intersection between segment and z=z0 - R = intlinf(self%z(seg(i)), self%R(seg(i)), & - self%z(seg(i)+1), self%R(seg(i)+1), z0) + R = linear_interp(self%z(seg(i)), self%R(seg(i)), & + self%z(seg(i)+1), self%R(seg(i)+1), z0) if (R < R0) inside = .not. inside end do end function contour_contains diff --git a/src/utils.f90 b/src/utils.f90 index d6dfd8e..7e38e08 100644 --- a/src/utils.f90 +++ b/src/utils.f90 @@ -6,83 +6,38 @@ module utils contains - pure function locatef(a,n,x) result(j) -! Given an array a(n), and a value x, with a(n) monotonic, either -! increasing or decreasing, returns a value j such that -! a(j) < x <= a(j+1) for a increasing, and such that -! a(j+1) < x <= a(j) for a decreasing. -! j=0 or j=n indicate that x is out of range (Numerical Recipes) - integer, intent(in) :: n - real(wp_), dimension(n), intent(in) :: a - real(wp_), intent(in) :: x - integer :: j - integer :: jl,ju,jm - logical :: incr - jl=0 - ju=n+1 - incr=a(n)>a(1) - do while ((ju-jl)>1) - jm=(ju+jl)/2 - if(incr.eqv.(x>a(jm))) then - jl=jm - else - ju=jm - endif - end do - j=jl - end function locatef + pure subroutine locate(array, value, index) + ! Given an `array`, and a `value` returns the `index` + ! such that `array(index) < value < array(index+1)` + ! + ! Notes: + ! 1. `array` must be monotonic, either increasing or decreasing. + ! 2. an index equal to 0 or `size(array)` indicate the value was not found + ! Source: Numerical Recipes - pure subroutine locate(xx,n,x,j) - integer, intent(in) :: n - real(wp_), intent(in) :: xx(n), x - integer, intent(out) :: j - integer :: jl,ju,jm + ! subroutine arguments + real(wp_), intent(in) :: array(:), value + integer, intent(out) :: index + + ! local variables + integer :: jl, ju, jm logical :: incr -! -! Given an array xx(n), and a value x -! returns a value j such that xx(j) < x < xx(j+1) -! xx(n) must be monotonic, either increasing or decreasing. -! j=0 or j=n indicate that x is out of range (Numerical Recipes) -! - jl=0 - ju=n+1 - incr=xx(n)>xx(1) - do while ((ju-jl)>1) - jm=(ju+jl)/2 - if(incr .eqv. (x>xx(jm))) then - jl=jm + + jl = 0 + ju = size(array) + 1 + incr = array(size(array)) > array(1) + + do while ((ju - jl) > 1) + jm = (ju + jl)/2 + if (incr .eqv. (value > array(jm))) then + jl = jm else - ju=jm + ju = jm endif end do - j=jl + index = jl end subroutine locate - pure subroutine locatex(xx,n,n1,n2,x,j) - integer, intent(in) :: n,n1,n2 - real(wp_), intent(in) :: xx(n), x - integer, intent(out) :: j - integer :: jl,ju,jm -! -! Given an array xx(n), and a value x -! returns a value j such that xx(j) < x < xx(j+1) -! xx(n) must be monotonic, either increasing or decreasing. -! j=n1-1or j=n2+1 indicate that x is out of range -! modified from subr. locate (Numerical Recipes) -! - jl=n1-1 - ju=n2+1 - do while ((ju-jl)>1) - jm=(ju+jl)/2 - if((xx(n2)>xx(n1)) .eqv. (x>xx(jm))) then - jl=jm - else - ju=jm - endif - end do - j=jl - end subroutine locatex - pure subroutine locate_unord(array, value, locs, n, nlocs) ! Given an `array` of size `n` and a `value`, finds at most @@ -117,233 +72,74 @@ contains end subroutine locate_unord - pure function intlinf(x1,y1,x2,y2,x) result(y) - !linear interpolation - !must be x1 != x2 - use const_and_precisions, only : one - real(wp_),intent(in) :: x1,y1,x2,y2,x + pure function linear_interp(x1, y1, x2, y2, x) result(y) + ! Computes the linear interpolation between + ! the points `(x1, y1)` and `(x2, y2)` at `x`. + ! + ! Note: x1 != x2 is assumed + + ! function arguments + real(wp_), intent(in) :: x1, y1, x2, y2, x real(wp_) :: y + + ! local variables real(wp_) :: a - a=(x2-x)/(x2-x1) - y=a*y1+(one-a)*y2 - end function intlinf - pure subroutine intlin(x1,y1,x2,y2,x,y) - real(wp_), intent(in) :: x1,y1,x2,y2,x - real(wp_), intent(out) :: y - real(wp_) :: dx,aa,bb -! -! linear interpolation -! (x1,y1) < (x,y) < (x2,y2) -! - dx=x2-x1 - aa=(x2-x)/dx - bb=1.0_wp_-aa - y=aa*y1+bb*y2 - end subroutine intlin + a = (x2 - x)/(x2 - x1) + y = a*y1 + (1 - a)*y2 + end function linear_interp - pure subroutine vmax(x,n,xmax,imx) - integer, intent(in) :: n - real(wp_), intent(in) :: x(n) - real(wp_), intent(out) :: xmax - integer, intent(out) :: imx + + pure subroutine vmaxmin(x, xmin, xmax, imin, imax) + ! Computes the maximum and minimum of the array `x` + ! and optionally their indices + + ! subroutine arguments + real(wp_), intent(in) :: x(:) + real(wp_), intent(out) :: xmin, xmax + integer, intent(out), optional :: imin, imax + + ! local variables integer :: i - if (n<1) then - imx=0 + if (size(x) < 1) then + if (present(imin)) imin = 0 + if (present(imax)) imax = 0 return end if - imx=1 - xmax=x(1) - do i=2,n - if(x(i)>xmax) then - xmax=x(i) - imx=i - end if - end do - end subroutine vmax - pure subroutine vmin(x,n,xmin,imn) - integer, intent(in) :: n - real(wp_), intent(in) :: x(n) - real(wp_), intent(out) :: xmin - integer, intent(out) :: imn - integer :: i + if (present(imin)) imin = 1 + if (present(imax)) imax = 1 + xmin = x(1) + xmax = x(1) - if (n<1) then - imn=0 - return - end if - imn=1 - xmin=x(1) - do i=2,n - if(x(i)xmax) then - xmax=x(i) - imx=i - end if - end do - end subroutine vmaxmini - - pure subroutine vmaxmin(x,n,xmin,xmax) - integer, intent(in) :: n - real(wp_), intent(in) :: x(n) - real(wp_), intent(out) :: xmin, xmax - integer :: i - - if (n<1) then - return - end if - xmin=x(1) - xmax=x(1) - do i=2,n - if(x(i)xmax) then - xmax=x(i) + do i = 2, size(x) + if(x(i) < xmin) then + xmin = x(i) + if (present(imin)) imin = i + else if(x(i) > xmax) then + xmax = x(i) + if (present(imax)) imax = i end if end do end subroutine vmaxmin - pure subroutine order(p,q) -! returns p,q in ascending order - real(wp_), intent(inout) :: p,q + + pure subroutine sort_pair(p) + ! Sorts the pair `p` in ascending order + + ! subroutine inputs + real(wp_), intent(inout) :: p(2) + + ! local variables real(wp_) :: temp - if (p>q) then - temp=p - p=q - q=temp + + if (p(1) > p(2)) then + temp = p(1) + p(1) = p(2) + p(2) = temp end if - end subroutine order - - pure subroutine bubble(a,n) -! bubble sorting of array a - integer, intent(in) :: n - real(wp_), dimension(n), intent(inout) :: a - integer :: i, j - do i=1,n - do j=n,i+1,-1 - call order(a(j-1), a(j)) - end do - end do - end subroutine bubble - - - pure subroutine range2rect(xmin, xmax, ymin, ymax, x, y) - ! Given two ranges [xmin, xmax], [ymin, ymax] builds - ! the x and y vertices of the following rectangle: - ! - ! (xmin, ymax)╔═════╗(xmax, ymax) - ! ║4 3║ - ! ║ ║ - ! ║1 2║ - ! (xmin, ymin)╚═════╝(xmin, ymax) - ! - - ! subroutine arguments - real(wp_), intent(in) :: xmin, xmax, ymin, ymax - real(wp_), intent(out), dimension(5) :: x, y - - x = [xmin, xmax, xmax, xmin, xmin] - y = [ymin, ymin, ymax, ymax, ymin] - end subroutine range2rect - - - function inside(vertx, verty, x0, y0) - ! Tests whether the point (`x0`, `y0`) lies inside the - ! simple polygon of vertices `vertx`, `verty`. - - - ! subroutine arguments - real(wp_), dimension(:), intent(in) :: vertx, verty - real(wp_), intent(in) :: x0, y0 - logical :: inside - - ! local variables - integer :: seg(size(vertx)) - real(wp_) :: x, vertx_(size(vertx)+1), verty_(size(vertx)+1) - integer :: i, nsegs, n - - ! Ensure the first and last point are the same - n = size(vertx) - vertx_(1:n) = vertx(1:n) - verty_(1:n) = verty(1:n) - vertx_(n+1) = vertx(1) - verty_(n+1) = verty(1) - - inside = .false. - - ! Find the `nsegs` segments that intersect the horizontal - ! line y=y0, i.e. `verty(seg(i)) < y < verty(seg(i)+1)` - call locate_unord(verty_, y0, seg, n, nsegs) - - ! No intersections, it must be outside (above or below) - if (nsegs == 0) return - - ! Count the number of intersections that lie to the left - ! (equivalently, to the right) of the point. An even number - ! means that the point is outside the polygon. - do i = 1, nsegs - ! coordinate of the intersection between segment and y=y0 - x = intlinf(verty_(seg(i)), vertx_(seg(i)), & - verty_(seg(i)+1), vertx_(seg(i)+1), y0) - if (x < x0) inside = .not. inside - end do - end function inside - - - function get_free_unit(unit) result(i) - ! Returns `unit` back or the first free unit - ! number `i` if `unit` is absent. - ! When no unit is available, returns -1. - - ! function arguments - integer :: i - integer, intent(in), optional :: unit - - ! local variables - integer, parameter :: max_allowed = 999 - integer :: error - logical :: ex, op - - if (present(unit)) then - i = unit - return - end if - - do i=0,max_allowed - inquire(unit=i, exist=ex, opened=op, iostat=error) - ! if unit i exists and is free - if (error == 0 .and. ex .and. .not. op) return - end do - i = -1 - - end function get_free_unit + end subroutine sort_pair function dirname(filepath) result(directory)