From af5fb208b247e53807202e8034356fb071679206 Mon Sep 17 00:00:00 2001 From: Lorenzo Figini Date: Mon, 12 Oct 2015 17:09:28 +0000 Subject: [PATCH] inclusion of updated routines for analytical equilibrium/profiles; further cleaning of equidata --- Makefile | 6 +- src/coreprofiles.f90 | 81 ++- src/{interp_eqprof.f90 => equilibrium.f90} | 590 +++++++++++++++-- src/gray.f | 701 +++------------------ src/minpack.f90 | 596 +++++++++++++++++- 5 files changed, 1290 insertions(+), 684 deletions(-) rename src/{interp_eqprof.f90 => equilibrium.f90} (51%) diff --git a/Makefile b/Makefile index 3ff231f..ca70079 100644 --- a/Makefile +++ b/Makefile @@ -5,7 +5,7 @@ EXE=gray MAINOBJ=gray.o OTHOBJ=conical.o const_and_precisions.o coreprofiles.o dierckx.o dispersion.o \ eccd.o eierf.o graydata_anequil.o graydata_flags.o graydata_par.o \ - interp_eqprof.o magsurf_data.o math.o minpack.o numint.o quadpack.o \ + equilibrium.o magsurf_data.o math.o minpack.o numint.o quadpack.o \ reflections.o simplespline.o utils.o beamdata.o # Alternative search paths @@ -27,7 +27,7 @@ $(EXE): $(MAINOBJ) $(OTHOBJ) # Dependencies on modules gray.o: const_and_precisions.o coreprofiles.o dierckx.o dispersion.o eccd.o \ graydata_anequil.o graydata_flags.o graydata_par.o \ - interp_eqprof.o magsurf_data.o math.o minpack.o numint.o quadpack.o \ + equilibrium.o magsurf_data.o math.o numint.o quadpack.o \ reflections.o simplespline.o utils.o beamdata.o conical.o: const_and_precisions.o coreprofiles.o: const_and_precisions.o dierckx.o graydata_anequil.o \ @@ -39,7 +39,7 @@ eierf.o: const_and_precisions.o graydata_anequil.o: const_and_precisions.o graydata_flags.o: const_and_precisions.o graydata_par.o: const_and_precisions.o -interp_eqprof.o: const_and_precisions.o +equilibrium.o: const_and_precisions.o dierckx.o minpack.o simplespline.o utils.o magsurf_data.o: const_and_precisions.o math.o: const_and_precisions.o minpack.o: const_and_precisions.o diff --git a/src/coreprofiles.f90 b/src/coreprofiles.f90 index 9de7ef7..00c271d 100644 --- a/src/coreprofiles.f90 +++ b/src/coreprofiles.f90 @@ -6,12 +6,13 @@ module coreprofiles 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 subroutine density(psin,dens,ddens) use graydata_flags, only : iprof - use graydata_anequil, only : dens0,aln1,aln2 +! use graydata_anequil, only : dens0,aln1,aln2 use dierckx, only : splev,splder implicit none ! arguments @@ -29,18 +30,18 @@ contains ! dens=zero ddens=zero - if(psin.ge.psdbnd.or.psin.lt.zero) return + if((psin >= psdbnd).or.(psin < zero)) return ! - if(iprof.eq.0) then - if(psin.gt.one) return + if(iprof == 0) then + if(psin > one) return profd=(one-psin**aln1)**aln2 dens=dens0*profd dprofd=-aln1*aln2*psin**(aln1-one) & *(one-psin**aln1)**(aln2-one) ddens=dens0*dprofd else - if(psin.gt.psnpp) then -! + if(psin > psnpp) then + ! smooth interpolation for psnpp < psi < psdbnd ! dens = fp * fh ! fp: parabola matched at psi=psnpp with given profile density @@ -64,11 +65,11 @@ contains ier=0 call splder(tfn,nsfd,cfn,3,nu,xxs,ffs,1,wrkfd,ier) ddens=ffs(1) - if(ier.gt.0) print*,ier - if(abs(dens).lt.1.0e-10_wp_) dens=zero + if(ier > 0) print*,ier + if(abs(dens) < 1.0e-10_wp_) dens=zero end if -! if(dens.lt.zero) print*,' DENSITY NEGATIVE',dens - if(dens.lt.zero) then +! if(dens < zero) print*,' DENSITY NEGATIVE',dens + if(dens < zero) then dens=zero ddens=zero end if @@ -78,7 +79,7 @@ contains function temp(psin) use const_and_precisions, only : wp_,zero,one use graydata_flags, only : iprof - use graydata_anequil, only : te0,dte0,alt1,alt2 +! use graydata_anequil, only : te0,dte0,alt1,alt2 use utils, only : locate use simplespline, only :spli implicit none @@ -90,8 +91,8 @@ contains real(wp_) :: proft,dps temp=zero - if(psin.ge.one.or.psin.lt.zero) return - if(iprof.eq.0) then + if((psin >= one).or.(psin < zero)) return + if(iprof == 0) then proft=(1.0_wp_-psin**alt1)**alt2 temp=(te0-dte0)*proft+dte0 else @@ -105,7 +106,7 @@ contains function fzeff(psin) use const_and_precisions, only : wp_,zero,one use graydata_flags, only : iprof - use graydata_anequil, only : zeffan +! use graydata_anequil, only : zeffan use utils, only : locate use simplespline, only :spli implicit none @@ -117,8 +118,8 @@ contains real(wp_) :: dps fzeff=one - if(psin.gt.one.or.psin.lt.zero) return - if(iprof.eq.0) then + if((psin >= one).or.(psin < zero)) return + if(iprof == 0) then fzeff=zeffan else call locate(psrad,npp,psin,k) @@ -157,6 +158,34 @@ contains close(u) end subroutine read_profiles + subroutine read_profiles_an(filenm,te,ne,zeff,unit) + use utils, only : get_free_unit + 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 + + if (present(unit)) then + u=unit + else + u=get_free_unit() + end if + + if(allocated(te)) deallocate(te) + if(allocated(ne)) deallocate(ne) + if(allocated(zeff)) deallocate(zeff) + allocate(te(4),ne(3),zeff(1)) + + open(file=trim(filenm),status='old',unit=u) + read(u,*) ne(1:3) ! dens0,aln1,aln2 + read(u,*) te(1:4) ! te0,dte0,alt1,alt2 + read(u,*) zeff(1) ! zeffan + close(u) + end subroutine read_profiles_an + subroutine tene_scal(te,ne,tfact,nfact,bfact,iscal) implicit none ! arguments @@ -246,10 +275,10 @@ contains ddenpp=ddedge(1) d2denpp=d2dedge(1) psdbnd=psnpp - if(ddenpp.lt.zero) then + if(ddenpp < zero) then xnv=psnpp-ddenpp/d2denpp ynv=denpp-0.5_wp_*ddenpp**2/d2denpp - if(d2denpp.gt.zero.and.ynv.ge.zero) then + if((d2denpp > zero).and.(ynv >= zero)) then psdbnd=xnv else psdbnd=xnv+sqrt((ddenpp/d2denpp)**2-2.0_wp_*denpp/d2denpp) @@ -271,4 +300,20 @@ contains if(allocated(cfn)) deallocate(cfn) end subroutine unset_prfspl + subroutine set_prfan(te,ne,zeff) + implicit none + REAL(wp_), dimension(:), intent(in) :: te,ne,zeff + + te0=te(1) + dte0=te(2) + alt1=te(3) + alt2=te(4) + dens0=ne(1) + aln1=ne(2) + aln2=ne(3) + zeffan=zeff(1) + + psdbnd=1.0_wp_ + end subroutine set_prfan + end module coreprofiles diff --git a/src/interp_eqprof.f90 b/src/equilibrium.f90 similarity index 51% rename from src/interp_eqprof.f90 rename to src/equilibrium.f90 index 01733c9..fab68ef 100644 --- a/src/interp_eqprof.f90 +++ b/src/equilibrium.f90 @@ -1,22 +1,13 @@ -module interp_eqprof +module equilibrium use const_and_precisions, only : wp_ implicit none -! equidata - ! === 1D array Fpol(psi), q(psi), rhotor(psi) ==================== - !INTEGER, SAVE :: npsieq - !REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: rhopnr - - ! === 1D array plasma boundary Rbnd_i, Zbnd_i ==================== - INTEGER, SAVE :: nbbbs - REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: rbbbs,zbbbs - REAL(wp_), SAVE :: btaxis,rmaxis,zmaxis - REAL(wp_), SAVE :: rmnm,rmxm,zmnm,zmxm,dr,dz - REAL(wp_), SAVE :: zbmin,zbmax - REAL(wp_), SAVE :: phitedge + REAL(wp_), SAVE :: btrcen ! used only for Jcd_ASTRA def. + REAL(wp_), SAVE :: rcen ! computed as fpol(a)/btrcen + REAL(wp_), SAVE :: rmnm,rmxm,zmnm,zmxm + REAL(wp_), SAVE :: zbinf,zbsup REAL(wp_), SAVE :: rup,zup,rlw,zlw - REAL(wp_), SAVE :: rcen,btrcen ! rcen used to compute btrcen from fpol, btrcen used only for Jcd_ASTRA def. INTEGER, PARAMETER :: kspl=3,ksplp=kspl+1 ! === 2D spline psi(R,z), normalization and derivatives ========== @@ -32,9 +23,11 @@ module interp_eqprof REAL(wp_), SAVE :: fpolas ! === 1D spline rhot(rhop), rhop(rhot), q(psi) =================== ! computed on psinr,rhopnr [,rhotnr] arrays - INTEGER, SAVE :: nq - REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: rhopq,psiq,q - REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: crhot,cq + INTEGER, SAVE :: nq,nrho + REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: psinr,rhopr,rhotr + REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: cq,crhop,crhot + REAL(wp_), SAVE :: phitedge,aminor + REAL(wp_), SAVE :: q0,qa,alq contains @@ -152,6 +145,52 @@ contains psia=psiedge-psiaxis if(in==0) psin=(psin-psiaxis)/psia end subroutine read_eqdsk + + subroutine read_equil_an(filenm,rv,zv,fpol,q,unit) + use utils, only : get_free_unit + implicit none +! arguments + character(len=*), intent(in) :: filenm + integer, optional, intent(in) :: unit +! integer, intent(in) :: isgnbphi +! real(wp_), intent(in) :: factb +! real(wp_), intent(out) :: rvac,rax,zax + real(wp_), dimension(:), allocatable, intent(out) :: rv,zv,fpol,q +! local variables + integer :: u + real(wp_) :: rr0m,zr0m,rpam,b0,q0,qa,alq !,rcen,btrcen + + if (present(unit)) then + u=unit + else + u=get_free_unit() + end if + open(file=trim(filenm),status='old',unit=u) + read(u,*) rr0m,zr0m,rpam + read(u,*) b0 + read(u,*) q0,qa,alq +! rcen=rr0m +! btrcen=b0 +! b0=isgnbphi*b0*factb +! rvac=rr0m +! rax=rr0m +! zax=zr0m +! zbmin=(zr0-rpa)/1.0e2_wp_ +! zbmax=(zr0+rpa)/1.0e2_wp_ + if(allocated(rv)) deallocate(rv) + if(allocated(zv)) deallocate(zv) + if(allocated(fpol)) deallocate(fpol) + if(allocated(q)) deallocate(q) + allocate(rv(2),zv(1),fpol(1),q(3)) + rv(1)=rr0m + rv(2)=rpam + zv(1)=zr0m + fpol(1)=b0*rr0m + q(1)=q0 + q(2)=qa + q(3)=alq + close(u) + end subroutine read_equil_an subroutine change_cocos(psia,fpol,q,cocosin,cocosout,ierr) use const_and_precisions, only : zero,one,pi @@ -237,24 +276,30 @@ contains bsign=int(sign(one,fpol(size(fpol)))) end subroutine eq_scal - subroutine set_eqspl(rv,zv,psin,psinr,fpol,sspl,ssfp) + subroutine set_eqspl(rv,zv,psin,psiwbrad,psinr,fpol,sspl,ssfp, & + r0,rax,zax,rbnd,zbnd,ixp) use const_and_precisions, only : zero,one use dierckx, only : regrid,coeff_parder,curfit,splev + use utils, only : vmaxmin,vmaxmini implicit none ! arguments real(wp_), dimension(:), intent(in) :: rv,zv,psinr,fpol real(wp_), dimension(:,:), intent(in) :: psin + real(wp_), intent(in) :: psiwbrad real(wp_), intent(inout) :: sspl,ssfp + real(wp_), intent(in), optional :: r0,rax,zax + real(wp_), dimension(:), intent(in), optional :: rbnd,zbnd + integer, intent(in), optional :: ixp ! local constants integer, parameter :: iopt=0 +! local variables integer :: liwrk,lwrk,lw10,lw01,lw20,lw02,lw11,lwrkf - integer :: nr,nz,nrest,nzest,npsest,nrz,npsi - real(wp_) :: fp + integer :: nr,nz,nrest,nzest,npsest,nrz,npsi,nbnd,ibinf,ibsup + real(wp_) :: fp,rax0,zax0,psinoptmp,psinxptmp,rbmin,rbmax,rbinf,rbsup,r1,z1 real(wp_), dimension(1) :: fpoli real(wp_), dimension(:), allocatable :: fvpsi,wf,wrk integer, dimension(:), allocatable :: iwrk -! local variables - integer :: ier + integer :: ier,ixploc,info ! ! compute array sizes and prepare working space arrays nr=size(rv) @@ -328,10 +373,118 @@ contains call curfit(iopt,npsi,psinr,fpol,wf,zero,one,kspl,ssfp,npsest,nsf, & tfp,cfp,fp,wrk(1:lwrkf),lwrkf,iwrk(1:npsest),ier) call splev(tfp,nsf,cfp,kspl,psinr(npsi:npsi),fpoli,1,ier) +! set vacuum value used outside 0<=psin<=1 range fpolas=fpoli(1) - btrcen=fpolas/rcen ! free temporary arrays deallocate(wrk,iwrk,wf) +! +! re-normalize psi after spline computation +! +! start with un-corrected psi +! + psia=psiwbrad + psinop=0.0_wp_ + psiant=1.0_wp_ +! +! use provided boundary to set an initial guess for the search of O/X points +! + nbnd=0 + if(present(rbnd).and.present(zbnd)) then + nbnd=min(size(rbnd),size(zbnd)) + end if + if (nbnd>0) then + call vmaxmini(zbnd,nbnd,zbinf,zbsup,ibinf,ibsup) + rbinf=rbnd(ibinf) + rbsup=rbnd(ibsup) + call vmaxmin(rbnd,nbnd,rbmin,rbmax) + else + zbinf=zv(2) + zbsup=zv(nz-1) + rbinf=rv((nr+1)/2) + rbsup=rbinf + rbmin=rv(2) + rbmax=rv(nr-1) + end if +! +! search for exact location of the magnetic axis +! + if(present(rax)) then + rax0=rax + else + rax0=0.5_wp_*(rbmin+rbmax) + end if + if(present(zax)) then + zax0=zax + else + zax0=0.5_wp_*(zbinf+zbsup) + end if + call points_ox(rax0,zax0,rmaxis,zmaxis,psinoptmp,info) + print'(a,2f8.4,es12.5)','O-point',rmaxis,zmaxis,psinoptmp +! +! search for X-point if ixp not = 0 +! + if(present(ixp)) then + ixploc=ixp + else + ixploc=0 + end if + if(ixploc/=0) then + 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 + zbinf=z1 + psinop=psinoptmp + psiant=psinxptmp-psinop + call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbsup),r1,z1,one,info) + zbsup=z1 + else + ixploc=0 + end if + 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 + zbsup=z1 + psinop=psinoptmp + psiant=psinxptmp-psinop + call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbinf),r1,z1,one,info) + zbinf=z1 + else + ixploc=0 + end if + end if + end if + + if (ixploc==0) then + psinop=psinoptmp + psiant=one-psinop +! find upper horizontal tangent point + call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbsup),r1,z1,one,info) + zbsup=z1 + rbsup=r1 +! find lower horizontal tangent point + 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 + end if + print*,' ' +! +! save Bt value on axis (required in flux_average and used in Jcd def) +! and vacuum value B0 at ref. radius R0 (used in Jcd_astra def) +! + call equinum_fpol(0.0_wp_,btaxis) + btaxis=btaxis/rmaxis + if(present(r0)) then + btrcen=fpolas/r0 + rcen=r0 + else + btrcen=fpolas/rmaxis + rcen=rmaxis + end if + print'(a,f8.4)','BT_centr= ',btrcen + print'(a,f8.4)','BT_axis = ',btaxis end subroutine set_eqspl subroutine unset_eqspl @@ -346,14 +499,11 @@ contains if(allocated(cceq02)) deallocate(cceq02) if(allocated(cceq20)) deallocate(cceq20) if(allocated(cceq11)) deallocate(cceq11) + nsr=0 + nsz=0 + nsf=0 end subroutine unset_eqspl - subroutine dealloc_bnd - implicit none - if(allocated(rbbbs)) deallocate(rbbbs) - if(allocated(zbbbs)) deallocate(zbbbs) - end subroutine dealloc_bnd - subroutine equinum_psi(rpsim,zpsim,psinv,dpsidr,dpsidz, & ddpsidrr,ddpsidzz,ddpsidrz) use dierckx, only : fpbisp @@ -454,32 +604,6 @@ contains end if end subroutine equinum_fpol -! subroutine btor(rpsim,zpsim,bphi) -! implicit none -!! arguments -! real(wp_), intent(in) :: rpsim,zpsim -! real(wp_), intent(out) :: bphi -!! local variables -! real(wp_) :: psinv,fpolv -! -! call equinum_psi(rpsim,zpsim,psinv) -! call equinum_fpol(psinv,fpolv) -! bphi=fpolv/rpsim -! end subroutine btor - -! subroutine bpol(rpsim,zpsim,brr,bzz) -! implicit none -!! arguments -! real(wp_), intent(in) :: rpsim,zpsim -! real(wp_), intent(out) :: brr,bzz -!! local variables -! real(wp_) :: dpsidr,dpsidz -! -! call equinum_psi(rpsim,zpsim,dpsidr=dpsidr,dpsidz=dpsidz) -! brr=-dpsidz/rpsim -! bzz= dpsidr/rpsim -! end subroutine bpol - subroutine bfield(rpsim,zpsim,bphi,br,bz) implicit none ! arguments @@ -498,4 +622,358 @@ contains if (present(bz)) bz= bz/rpsim end subroutine bfield -end module interp_eqprof + subroutine setqphi_num(psinq,q,psia,rhotn) + use const_and_precisions, only : pi + use simplespline, only : difcs + implicit none +! arguments + real(wp_), dimension(:), intent(in) :: psinq,q + real(wp_), intent(in) :: psia + real(wp_), dimension(:), intent(out), optional :: rhotn +! local variables + real(wp_), dimension(size(q)) :: phit + real(wp_) :: dx + integer, parameter :: iopt=0 + integer :: k,ier + + nq=size(q) + if(allocated(psinr)) deallocate(psinr) + if(allocated(cq)) deallocate(cq) + allocate(psinr(nq),cq(nq,4)) + + psinr=psinq + call difcs(psinr,q,nq,iopt,cq,ier) +! +! Toroidal flux phi = 2*pi*Integral q dpsi +! + phit(1)=0.0_wp_ + do k=1,nq-1 + dx=psinr(k+1)-psinr(k) + phit(k+1)=phit(k) + dx*(cq(k,1) + dx*(cq(k,2)/2.0_wp_ + & + dx*(cq(k,3)/3.0_wp_ + dx* cq(k,4)/4.0_wp_) ) ) + end do + phitedge=phit(nq) + if(present(rhotn)) rhotn(1:nq)=sqrt(phit/phitedge) + phitedge=2*pi*psia*phitedge + end subroutine setqphi_num + + subroutine unset_q + implicit none + + if(allocated(psinr)) deallocate(psinr) + if(allocated(cq)) deallocate(cq) + nq=0 + end subroutine unset_q + + function fq(psin) + use const_and_precisions, only : wp_ + use simplespline, only :spli + use utils, only : locate + implicit none +! arguments + real(wp_), intent(in) :: psin + real(wp_) :: fq +! local variables + integer :: i + real(wp_) :: dps + + call locate(psinr,nq,psin,i) + i=min(max(1,i),nq-1) + dps=psin-psinr(i) + fq=spli(cq,nq,i,dps) + end function fq + + subroutine set_rhospl(rhop,rhot) + use simplespline, only : difcs + implicit none +! arguments + real(wp_), dimension(:), intent(in) :: rhop, rhot +! local variables + integer, parameter :: iopt=0 + integer :: ier + + nrho=size(rhop) + + if(allocated(rhopr)) deallocate(rhopr) + if(allocated(rhotr)) deallocate(rhotr) + if(allocated(crhop)) deallocate(crhop) + if(allocated(crhot)) deallocate(crhot) + allocate(rhopr(nrho),rhotr(nrho),crhop(nrho,4),crhot(nrho,4)) + + rhopr=rhop + rhotr=rhot + call difcs(rhotr,rhopr,nrho,iopt,crhop,ier) + call difcs(rhopr,rhotr,nrho,iopt,crhot,ier) + end subroutine set_rhospl + + subroutine unset_rhospl + implicit none + + if(allocated(rhopr)) deallocate(rhopr) + if(allocated(rhotr)) deallocate(rhotr) + if(allocated(crhop)) deallocate(crhop) + if(allocated(crhot)) deallocate(crhot) + nrho=0 + end subroutine unset_rhospl + + function frhopol(rhot) + use utils, only : locate + use simplespline, only : spli + implicit none +! arguments + real(wp_), intent(in) :: rhot + real(wp_) :: frhopol +! local variables + integer :: i + real(wp_) :: dr + + call locate(rhotr,nrho,rhot,i) + i=min(max(1,i),nrho-1) + dr=rhot-rhotr(i) + frhopol=spli(crhop,nrho,i,dr) + end function frhopol + + function frhotor(rhop) + use utils, only : locate + use simplespline, only : spli + implicit none +! arguments + real(wp_), intent(in) :: rhop + real(wp_) :: frhotor +! local variables + integer :: i + real(wp_) :: dr + + call locate(rhopr,nrho,rhop,i) + i=min(max(1,i),nrho-1) + dr=rhop-rhopr(i) + frhotor=spli(crhot,nrho,i,dr) + end function frhotor + + subroutine points_ox(rz,zz,rf,zf,psinvf,info) + use const_and_precisions, only : comp_eps + use minpack, only : hybrj1 + implicit none +! local constants + integer, parameter :: n=2,ldfjac=n,lwa=(n*(n+13))/2 +! arguments + real(wp_), intent(in) :: rz,zz + real(wp_), intent(out) :: rf,zf,psinvf + integer, intent(out) :: info +! local variables + real(wp_) :: tol + real(wp_), dimension(n) :: xvec,fvec + real(wp_), dimension(lwa) :: wa + real(wp_), dimension(ldfjac,n) :: fjac + + 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) + end subroutine points_ox + + subroutine fcnox(n,x,fvec,fjac,ldfjac,iflag) + implicit none +! 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 + real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz +! + select case(iflag) + case(1) + call equinum_psi(x(1),x(2),dpsidr=dpsidr,dpsidz=dpsidz) + fvec(1) = dpsidr/psia + fvec(2) = dpsidz/psia + case(2) + call equinum_psi(x(1),x(2),ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz, & + ddpsidrz=ddpsidrz) + fjac(1,1) = ddpsidrr/psia + fjac(1,2) = ddpsidrz/psia + fjac(2,1) = ddpsidrz/psia + fjac(2,2) = ddpsidzz/psia + case default + print*,'iflag undefined' + end select + end subroutine fcnox + + subroutine points_tgo(rz,zz,rf,zf,psin0,info) + use const_and_precisions, only : comp_eps + use minpack, only : hybrj1mv + implicit none +! local constants + integer, parameter :: n=2,ldfjac=n,lwa=(n*(n+13))/2 +! arguments + real(wp_), intent(in) :: rz,zz,psin0 + real(wp_), intent(out) :: rf,zf + integer, intent(out) :: info +! local variables + real(wp_) :: tol + real(wp_), dimension(n) :: xvec,fvec,f0 + real(wp_), dimension(lwa) :: wa + real(wp_), dimension(ldfjac,n) :: fjac +! common/external functions/variables + + xvec(1)=rz + xvec(2)=zz + f0(1)=psin0 + 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,2f8.4)',' info subr points_tgo =',info, & + ' R,z coord.',xvec + end if + rf=xvec(1) + zf=xvec(2) + end + + subroutine fcntgo(n,x,f0,fvec,fjac,ldfjac,iflag) + use const_and_precisions, only : wp_ + implicit none +! 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 + real(wp_) :: psinv,dpsidr,dpsidz,ddpsidrr,ddpsidrz + + select case(iflag) + case(1) + call equinum_psi(x(1),x(2),psinv,dpsidr) + fvec(1) = psinv-f0(1) + fvec(2) = dpsidr/psia-f0(2) + case(2) + call equinum_psi(x(1),x(2),dpsidr=dpsidr,dpsidz=dpsidz, & + ddpsidrr=ddpsidrr,ddpsidrz=ddpsidrz) + fjac(1,1) = dpsidr/psia + fjac(1,2) = dpsidz/psia + fjac(2,1) = ddpsidrr/psia + fjac(2,2) = ddpsidrz/psia + case default + print*,'iflag undefined' + end select + end subroutine fcntgo + + subroutine set_equian(rax,zax,a,bax,qax,q1,qexp,n) + use const_and_precisions, only : pi,zero,one + implicit none +! arguments + real(wp_), intent(in) :: rax,zax,a,bax,qax,q1,qexp + integer, intent(in), optional :: n +! local variables + integer, parameter :: nqdef=101 + integer :: i + real(wp_) :: dr,fq0,fq1,qq,res,rn + real(wp_), dimension(:), allocatable :: rhotn,rhopn + + btaxis=bax + rmaxis=rax + zmaxis=zax + btrcen=bax + rcen=rax + aminor=a + zbinf=zmaxis-a + zbsup=zmaxis+a + q0=qax + qa=q1 + alq=qexp + + if (present(n)) then + nq=n + else + nq=nqdef + end if + if (allocated(psinr)) deallocate(psinr) + allocate(psinr(nq),rhotn(nq),rhopn(nq)) + dr=one/(nq-1) + rhotn(1)=zero + psinr(1)=zero + res=zero + fq0=zero + do i=2,n + rn=(i-1)*dr + qq=q0+(q1-q0)*rn**qexp + fq1=rn/qq + res=res+0.5_wp_*(fq1+fq0)*dr + fq0=fq1 + rhotn(i)=rn + psinr(i)=res + end do + phitedge=btaxis*aminor**2 ! temporary + psia=res*phitedge + phitedge=pi*phitedge ! final + psinr=psinr/res + rhopn=sqrt(psinr) + call set_rhospl(rhopn,rhotn) + end subroutine set_equian + + subroutine equian(rrm,zzm,psinv,fpolv,dfpv,dpsidr,dpsidz, & + ddpsidrr,ddpsidzz,ddpsidrz) + use const_and_precisions, only : wp_ + implicit none +! arguments + real(wp_), intent(in) :: rrm,zzm + real(wp_), intent(out), optional :: psinv,fpolv,dfpv,dpsidr,dpsidz, & + ddpsidrr,ddpsidzz,ddpsidrz +! local variables + integer :: sgn + real(wp_) :: cst,dpsidrp,d2psidrp,dqq,qq,rn,rpm,snt,rhop,rhot +! real(wp_) :: frhopol + +! simple model for equilibrium: large aspect ratio +! outside plasma: analytical continuation, not solution Maxwell eqs + + rpm=sqrt((rrm-rmaxis)**2+(zzm-zmaxis)**2) !!! rpm==rho_tor[m], rn=rho_tor_norm + rn=rpm/aminor + + snt=0.0_wp_ + cst=1.0_wp_ + if (rpm > 0.0_wp_) then + snt=(zzm-zmaxis)/rpm + cst=(rrm-rmaxis)/rpm + end if + + if (present(psinv)) then + rhot=rn + if(rn <= 1.0_wp_) then + rhop=frhopol(rhot) + psinv=rhop*rhop + else + psinv=1.0_wp_+btaxis*aminor**2/2.0_wp_/psia/qa*(rn*rn-1.0_wp_) + rhop=sqrt(psinv) + end if + end if + + if(rn <= 1.0_wp_) then + qq=q0+(qa-q0)*rn**alq + dpsidrp=-btaxis*aminor*rn/qq + dqq=alq*(qa-q0)*rn**(alq-1.0_wp_) + d2psidrp=-btaxis*(1.0_wp_-rn*dqq/qq)/qq + else + dpsidrp=-btaxis*aminor/qa*rn + d2psidrp=-btaxis/qa + end if + + if(present(fpolv)) fpolv=btaxis*rmaxis + if(present(dfpv)) dfpv=0.0_wp_ + + if(present(dpsidr)) dpsidr=dpsidrp*cst + if(present(dpsidz)) dpsidz=dpsidrp*snt + if(present(ddpsidrr)) ddpsidrr=dpsidrp*snt**2/rpm+cst**2*d2psidrp + if(present(ddpsidrz)) ddpsidrz=cst*snt*(d2psidrp-dpsidrp/rpm) + if(present(ddpsidzz)) ddpsidzz=dpsidrp*cst**2/rpm+snt**2*d2psidrp + + end subroutine equian + +end module equilibrium diff --git a/src/gray.f b/src/gray.f index 4cf3aee..fc7361c 100644 --- a/src/gray.f +++ b/src/gray.f @@ -206,7 +206,7 @@ c use reflections, only : inside,rlim,zlim,nlim,wall_refl use graydata_flags, only : iwarm,istpr0,istpl0,dst,ipass,igrad use graydata_par, only : psipol0,chipol0 - use interp_eqprof, only : zbmin,zbmax + use equilibrium, only : zbinf,zbsup use beamdata, only : nrayr,nrayth,psjki,ppabs,ccci,iiv,tauv,ihcd, . istore,iop,iow,tau1v,yyrfl,ywrk,ext,eyt implicit none @@ -283,7 +283,7 @@ c exit if (iwarm.gt.0.and.i.gt.1.and.ihcd(j,k).gt.0) then if(psinv.ge.0.and.psinv.le.1.0_wp_.and. - . zzm.ge.zbmin.and.zzm.le.zbmax) then + . zzm.ge.zbinf.and.zzm.le.zbsup) then call pabs_curr(i,j,k) iiv(j,k)=i else @@ -490,6 +490,7 @@ c use coreprofiles, only : psdbnd use beamdata, only : ywrk,psjki,tauv,alphav,pdjki, . currj,didst,q,tau1v,nrayr !,nrayth + use equilibrium, only : frhotor implicit none c local constants integer, parameter :: ndim=6 @@ -504,7 +505,7 @@ c common/external functions/variables integer :: nharm,nhf,iohkawa,index_rt,istpr,istpl real(wp_) :: p0mw,st,brr,bphi,bzz,ajphi,btot,xg,yg,dens,ddens, . tekev,alpha,effjcd,akim,tau0,anpl,anpr,ddr,an2s,an2,fdia,bdotgr, - . ddi,ddr11,anprr,anpri,frhotor + . ddi,ddr11,anprr,anpri complex(wp_) :: ex,ey,ez c common/nharm/nharm,nhf @@ -656,20 +657,20 @@ c . vc=>ccgs_,cvdr=>degree,pi use graydata_flags use graydata_par - use graydata_anequil - use coreprofiles, only : psdbnd,read_profiles,tene_scal,set_prfspl - use interp_eqprof, only : rmxm,zbmin,zbmax,rcen,psia,rmaxis, - . zmaxis,rbbbs,zbbbs,read_eqdsk,change_cocos,eq_scal,set_eqspl, - . btrcen,psinr=>psiq,qpsi=>q,equinum_fpol + use coreprofiles, only : psdbnd,read_profiles,tene_scal, + . set_prfspl,set_prfan + use equilibrium, only : read_eqdsk,change_cocos,eq_scal,set_eqspl, + . rmxm,set_rhospl,setqphi_num,set_equian use reflections, only : rlim,zlim,nlim,alloc_lim use beamdata, only : nrayr,nrayth,nstep implicit none c local variables - integer :: nfil,iox,ierr,isgniphi,isgnbphi,ifreefmt - real(wp_) :: dummy,bresg,r00,anr0c,anphi0c,fghz, - . w0csi,w0eta,d0csi,d0eta,phi0,zrcsi,zreta,ssplf + integer :: nfil,iox,ierr,isgniphi,isgnbphi,ifreefmt,idesc + real(wp_) :: dummy,bresg,r00,anr0c,anphi0c,fghz,w0csi,w0eta,d0csi, + . d0eta,phi0,zrcsi,zreta,ssplf,rax,zax,rvac,psia,rr0m,zr0m,rpam,b0, + . q0,qa,alq,te0,dte0,alt1,alt2,dens0,aln1,aln2,zeffan real(wp_), dimension(:), allocatable :: psrad,terad,derad,zfc, - . rv,zv,fpol + . rv,zv,fpol,qpsi,psinr,rhotn,rbbbs,zbbbs real(wp_), dimension(:,:), allocatable :: psin character(len=8) :: wdat character(len=10) :: wtim @@ -776,7 +777,7 @@ c ipsinorm 0 standard EQDSK format, 1 format Portone summer 2004 c sspl spline parameter for psi interpolation c read(2,*) filenmeqq - read(2,*) ipsinorm,sspl,ifreefmt + read(2,*) ipsinorm,sspl,ifreefmt,idesc c c factb factor for magnetic field (only for numerical equil) c scaling adopted: beta=const, qpsi=const, nustar=const @@ -800,9 +801,26 @@ c read parameters for analytical density and temperature c profiles if iprof=0 c if(iprof.eq.0) then - read(2,*) dens0,aln1,aln2 - read(2,*) te0,dte0,alt1,alt2 - read(2,*) zeffan + ! === to be replaced with call read_profiles_an === + if (allocated(terad)) deallocate(terad) + if (allocated(derad)) deallocate(derad) + if (allocated(zfc)) deallocate(zfc) + allocate(terad(4),derad(3),zfc(1)) + read(2,*) derad(1:3) ! dens0,aln1,aln2 + read(2,*) terad(1:4) ! te0,dte0,alt1,alt2 + read(2,*) zfc(1) ! zeffan + ! === end === + call tene_scal(terad(1:2),derad(1:1),factt,factn,factb,iscal) + call set_prfan(terad,derad,zfc) + te0=terad(1) ! only for print in headers.txt + dte0=terad(2) ! only for print in headers.txt + alt1=terad(3) ! only for print in headers.txt + alt2=terad(4) ! only for print in headers.txt + dens0=derad(1) ! only for print in headers.txt + aln1=derad(2) ! only for print in headers.txt + aln2=derad(3) ! only for print in headers.txt + zeffan=zfc(1) ! only for print in headers.txt + deallocate(terad,derad,zfc) else read(2,*) dummy,dummy,dummy read(2,*) dummy,dummy,dummy,dummy @@ -813,17 +831,11 @@ c read parameters for analytical simple cylindical equilibrium c if iequil=0,1 if(iequil.lt.2) then - read(2,*) rr0,zr0,rpa + read(2,*) rr0m,zr0m,rpam read(2,*) b0 read(2,*) q0,qa,alq - rr0m=rr0/1.0e2_wp_ - zr0m=zr0/1.0e2_wp_ - rpam=rpa/1.0e2_wp_ - rcen=rr0m - btrcen=b0 - zbmin=(zr0-rpa)/100.0_wp_ - zbmax=(zr0+rpa)/100.0_wp_ b0=b0*factb + call set_equian(rr0m,zr0m,rpam,b0,q0,qa,alq) call flux_average_an c call print_prof_an else @@ -924,16 +936,24 @@ c c read equilibrium data from file if iequil=2 c if (iequil.eq.2) then +c neqdsk=99 call read_eqdsk(trim(filenmeqq)//'.eqdsk',rv,zv,psin,psia, - . psinr,fpol,qpsi,rcen,rmaxis,zmaxis,rbbbs,zbbbs,rlim,zlim, - . ipsinorm,1,ifreefmt,neqdsk) !,idesc,ifreefmt,neqdsk) + . psinr,fpol,qpsi,rvac,rax,zax,rbbbs,zbbbs,rlim,zlim, + . ipsinorm,idesc,ifreefmt,neqdsk) !,idesc,ifreefmt,neqdsk) call change_cocos(psia,fpol,qpsi,icocos,3) call eq_scal(psia,fpol,isgniphi,isgnbphi,factb) +c ssplf=0.01_wp_ - call set_eqspl(rv,zv,psin,psinr,fpol,sspl,ssplf) + call set_eqspl(rv,zv,psin,psia,psinr,fpol,sspl,ssplf,rvac, + . rax,zax,rbbbs,zbbbs,ixp) + allocate(rhotn(size(qpsi))) + call setqphi_num(psinr,abs(qpsi),abs(psia),rhotn) + call set_rhospl(sqrt(psinr),rhotn) + deallocate(rhotn) +c nlim=size(rlim) - call equidata + call equidata(psinr,qpsi,size(qpsi)) c c locate btot=bres c @@ -977,7 +997,7 @@ c set simple limiter as two cylindrical walls at rwallm and r00 if (iequil.eq.2) then write(nfil,905) trim(filenmeqq) else - if (iequil.eq.1) write(nfil,912) rr0,zr0,rpa,B0,q0,qa,alq + if (iequil.eq.1) write(nfil,912) rr0m,zr0m,rpam,B0,q0,qa,alq if (iequil.eq.0) write(nfil,'("# VACUUM CASE")') end if if (iprof.eq.1) then @@ -1224,149 +1244,28 @@ c c c c - subroutine equidata + subroutine equidata(psinq,qpsi,nq) use const_and_precisions, only : wp_,pi use utils, only : vmaxmini - use graydata_flags, only : ixp use graydata_par, only : factb - use interp_eqprof, only : - . psia,psiant,psinop,btrcen,btaxis,rmaxis,zmaxis,rmnm, - . rmxm,zmnm,zmxm,dr,dz,zbmin,zbmax,phitedge,rup,zup, - . rlw,zlw,rbbbs,zbbbs,nbbbs,qpsi=>q,nq,equinum_fpol + use equilibrium, only : rup,zup,rlw,zlw,rmaxis,zmaxis, + . zbinf,zbsup,frhotor implicit none +c arguments + integer, intent(in) :: nq + real(wp_), dimension(nq), intent(in) :: psinq,qpsi c local variables integer :: i,info,iqmn,iqmx - real(wp_) :: psinxp,rmop,zmop,psinoptmp,r10,z10, - . rxp,zxp,psinxptmp,r1,z1,psin1,rbmin,rbmax,q2,q15,qmin,qmax, + real(wp_) :: rbmin,rbmax,q2,q15,qmin,qmax, . rhot15,psi15,rhot2,psi2,rhotsx -c common/external functions/variables - real(wp_) :: frhotor -c - nbbbs=size(rbbbs) - nq=size(qpsi) -c -c compute max and min z of last closed surface -c - rbmin=rmaxis - rbmax=rmaxis - if (nbbbs.gt.1) then - zbmin=1.0e+30_wp_ - zbmax=-1.0e+30_wp_ - do i=1,nbbbs - if(zbbbs(i).le.zbmin) then - zbmin=zbbbs(i) - rbmin=rbbbs(i) - end if - if(zbbbs(i).ge.zbmax) then - zbmax=zbbbs(i) - rbmax=rbbbs(i) - end if - end do - else - zbmin=-1.0e+30_wp_ - zbmax=1.0e+30_wp_ - end if - if(zbmin.le.zmnm) zbmin=zbmin+dz - if(rbmin.le.rmnm) rbmin=rbmin+dr - if(zbmax.ge.zmxm) zbmax=zbmax-dz - if(rbmax.ge.rmxm) rbmax=rbmax-dr -c -c start with uncorrected normalized psi -c - psinop=0.0_wp_ - psinxp=1.0_wp_ - psiant=1.0_wp_ -c -c search for O-point -c - call points_ox(rmaxis,zmaxis,rmop,zmop,psinoptmp,info) - rmaxis=rmop - zmaxis=zmop - print'(a,2f8.4,es12.5)','O-point',rmop,zmop,psinoptmp -c -c search for X-point if ixp not = 0 -c - if(ixp.ne.0) then - if(ixp.lt.0) then - r10=rbmin - z10=zbmin - call points_ox(r10,z10,rxp,zxp,psinxptmp,info) - if(psinxp.ne.-1.0_wp_) then - print'(a,2f8.4,es12.5)','X-point',rxp,zxp,psinxptmp - rbmin=rxp - zbmin=zxp - psinop=psinoptmp - psinxp=psinxptmp - psiant=psinxp-psinop - psin1=1.0_wp_ - r10=rmaxis - z10=(zbmax+zmaxis)/2.0_wp_ - call points_tgo(r10,z10,r1,z1,psin1,info) - rbmax=r1 - zbmax=z1 - else - ixp=0 -c print'(a)','no X-point' - end if - else - r10=rmop - z10=zbmax - call points_ox(r10,z10,rxp,zxp,psinxptmp,info) - if(psinxp.ne.-1.0_wp_) then - print'(a,2f8.4,e16.8)','X-point',rxp,zxp,psinxptmp - zbmax=zxp - psinop=psinoptmp - psinxp=psinxptmp - psiant=psinxp-psinop - psin1=1.0_wp_ - z10=(zbmin+zmaxis)/2.0_wp_ - call points_tgo(r10,z10,r1,z1,psin1,info) - zbmin=z1 - else - ixp=0 -c print'(a)','no X-point' - end if - end if - end if -c - if (ixp.eq.0) then - psin1=1.0_wp_ - psinop=psinoptmp - psiant=psin1-psinop - r10=rmaxis - z10=(zbmax+zmaxis)/2.0_wp_ - call points_tgo(r10,z10,r1,z1,psin1,info) - zbmax=z1 - rbmax=r1 - z10=(zbmin+zmaxis)/2.0_wp_ - call points_tgo(r10,z10,r1,z1,psin1,info) - zbmin=z1 - rbmin=r1 - print'(a,4f8.4)','no X-point ',rbmin,zbmin,rbmax,zbmax - end if - print*,' ' -c compute B_toroidal on axis - - call equinum_fpol(0.0_wp_,btaxis) - btaxis=btaxis/rmaxis - print'(a,f8.4)','factb = ',factb - print'(a,f8.4)','BT_centr= ',btrcen - print'(a,f8.4)','BT_axis = ',btaxis - -c compute normalized rho_tor from eqdsk q profile - call rhotor(rhotsx) - phitedge=abs(psia)*rhotsx*2*pi -c rrtor=sqrt(phitedge/abs(btrcen)/pi) - call rhopol -c print*,rhotsx,phitedge,rrtor,abs(psia) c compute flux surface averaged quantities rup=rmaxis rlw=rmaxis - zup=zmaxis+(zbmax-zmaxis)/10.0_wp_ - zlw=zmaxis-(zmaxis-zbmin)/10.0_wp_ + zup=zmaxis+(zbsup-zmaxis)/10.0_wp_ + zlw=zmaxis-(zmaxis-zbinf)/10.0_wp_ call flux_average c ipr=1 c call contours_psi(1.0_wp_,rcon,zcon,ipr) @@ -1378,19 +1277,19 @@ c locate psi surface for q=1.5 and q=2 rup=rmaxis rlw=rmaxis - zup=(zbmax+zmaxis)/2.0_wp_ - zlw=(zmaxis+zbmin)/2.0_wp_ + zup=(zbsup+zmaxis)/2.0_wp_ + zlw=(zmaxis+zbinf)/2.0_wp_ q2=2.0_wp_ q15=1.5_wp_ call vmaxmini(abs(qpsi),nq,qmin,qmax,iqmn,iqmx) if (q15.gt.qmin.and.q15.lt.qmax) then - call surfq(q15,psi15) + call surfq(psinq,qpsi,nq,q15,psi15) rhot15=frhotor(sqrt(psi15)) print'(3(a,f8.5))','psi_1.5 = ',psi15,' rhop_1.5 = ' . ,sqrt(psi15),' rhot_1.5 = ',rhot15 end if if (q2.gt.qmin.and.q2.lt.qmax) then - call surfq(q2,psi2) + call surfq(psinq,qpsi,nq,q2,psi2) rhot2=frhotor(sqrt(psi2)) print'(3(a,f8.5))','psi_2 = ',psi2,' rhop_2 = ' . ,sqrt(psi2),' rhot_2 = ',rhot2 @@ -1400,195 +1299,33 @@ c end c c -c - subroutine points_ox(rz,zz,rf,zf,psinvf,info) - use const_and_precisions, only : wp_,comp_eps - use interp_eqprof, only : equinum_psi - use minpack, only : hybrj1 - implicit none -c local constants - integer, parameter :: n=2,ldfjac=n,lwa=(n*(n+13))/2 -c arguments - integer :: info - real(wp_) :: rz,zz,rf,zf,psinvf -c local variables - real(wp_) :: tol - real(wp_), dimension(n) :: xvec,fvec - real(wp_), dimension(lwa) :: wa - real(wp_), dimension(ldfjac,n) :: fjac - real(wp_) :: psinv -c - external fcnox -c - 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,psinv) - psinvf=psinv - return -c - end -c -c -c - subroutine fcnox(n,x,fvec,fjac,ldfjac,iflag) - use const_and_precisions, only : wp_ - use interp_eqprof, only: psia,equinum_psi - implicit none -c arguments - integer :: n,iflag,ldfjac - real(wp_), dimension(n) :: x,fvec - real(wp_), dimension(ldfjac,n) :: fjac -c local variables - real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz -c - select case(iflag) - case(1) - call equinum_psi(x(1),x(2),dpsidr=dpsidr,dpsidz=dpsidz) - fvec(1) = dpsidr/psia - fvec(2) = dpsidz/psia - case(2) - call equinum_psi(x(1),x(2),ddpsidrr=ddpsidrr, - . ddpsidzz=ddpsidzz,ddpsidrz=ddpsidrz) - fjac(1,1) = ddpsidrr/psia - fjac(1,2) = ddpsidrz/psia - fjac(2,1) = ddpsidrz/psia - fjac(2,2) = ddpsidzz/psia - case default - print*,'iflag undefined' - end select - - return - end -c -c -c - subroutine points_tgo(rz,zz,rf,zf,psin,info) - use const_and_precisions, only : wp_,comp_eps - use minpack, only : hybrj1 - implicit none -c local constants - integer, parameter :: n=2,ldfjac=n,lwa=(n*(n+13))/2 -c arguments - integer :: info - real(wp_) :: rz,zz,rf,zf,psin -c local variables - real(wp_) :: tol - real(wp_), dimension(n) :: xvec,fvec - real(wp_), dimension(lwa) :: wa - real(wp_), dimension(ldfjac,n) :: fjac -c common/external functions/variables - real(wp_) :: h -c - external fcntgo -c - common/cnpsi/h -c - h=psin - xvec(1)=rz - xvec(2)=zz - tol = sqrt(comp_eps) - call hybrj1(fcntgo,n,xvec,fvec,fjac,ldfjac,tol,info,wa,lwa) - if(info.gt.1) then - end if - rf=xvec(1) - zf=xvec(2) - return - end -c -c -c - subroutine fcntgo(n,x,fvec,fjac,ldfjac,iflag) - use const_and_precisions, only : wp_ - use interp_eqprof, only: psia, equinum_psi - implicit none -c arguments - integer :: n,ldfjac,iflag - real(wp_), dimension(n) :: x,fvec - real(wp_), dimension(ldfjac,n) :: fjac -c internal variables - real(wp_) :: psinv,dpsidr,dpsidz,ddpsidrr,ddpsidrz -c common/external functions/variables - real(wp_) :: h -c - common/cnpsi/h -c - select case(iflag) - case(1) - call equinum_psi(x(1),x(2),psinv,dpsidr) - fvec(1) = psinv-h - fvec(2) = dpsidr/psia - case(2) - call equinum_psi(x(1),x(2),dpsidr=dpsidr,dpsidz=dpsidz, - . ddpsidrr=ddpsidrr,ddpsidrz=ddpsidrz) - fjac(1,1) = dpsidr/psia - fjac(1,2) = dpsidz/psia - fjac(2,1) = ddpsidrr/psia - fjac(2,2) = ddpsidrz/psia - case default - print*,'iflag undefined' - end select -c - return - end -c -c c subroutine print_prof use const_and_precisions, only : wp_ - use interp_eqprof, only : psinr=>psiq,qpsi=>q,nq + use equilibrium, only : psinr,nq,fq,frhotor use coreprofiles, only : density, temp - use utils, only : intlin implicit none c local constants real(wp_), parameter :: eps=1.e-4_wp_ c local variables - integer :: i,ips,nst + integer :: i real(wp_) :: psin,rhop,rhot,ajphi,te,qq real(wp_) :: dens,ddens -c common/external functions/variables - real(wp_) :: frhotor c write(55,*) ' #psi rhot ne Te q Jphi' - psin=0.0_wp_ - rhop=0.0_wp_ - rhot=0.0_wp_ - call density(psin,dens,ddens) - call tor_curr_psi(eps,ajphi) - te=temp(psin) - qq=qpsi(1) -c - write(55,111) psin,rhot,dens,te,qq,ajphi*1.0e-6_wp_ -c - nst=nq - do i=2,nst - psin=dble(i-1)/dble(nst-1) + do i=1,nq + psin=psinr(i) rhop=sqrt(psin) c call density(psin,dens,ddens) te=temp(psin) -c - ips=int((nq-1)*psin+1) - if(i.lt.nst) then - call intlin(psinr(ips),qpsi(ips),psinr(ips+1),qpsi(ips+1), - . psin,qq) - else - qq=qpsi(nq) - end if + qq=fq(psin) rhot=frhotor(rhop) - call tor_curr_psi(psin,ajphi) - write(55,111) psin,rhot,dens,te,qq,ajphi*1.0e-6_wp_ + call tor_curr_psi(max(eps,psin),ajphi) + write(55,"(12(1x,e12.5))") psin,rhot,dens,te,qq,ajphi*1.e-6_wp_ end do c return - 111 format(12(1x,e12.5)) end c c @@ -1596,6 +1333,7 @@ c subroutine print_prof_an use const_and_precisions, only : wp_ use coreprofiles, only : density, temp + use equilibrium, only : frhotor implicit none c local constants integer, parameter :: nst=51 @@ -1603,8 +1341,6 @@ c local variables integer :: i real(wp_) :: psin,rhop,rhot,te real(wp_) :: dens,ddens -c common/external functions/variables - real(wp_) :: frhotor c write(55,*) ' #psi rhot ne Te' do i=1,nst @@ -1622,13 +1358,14 @@ c c c c - subroutine surfq(qval,psival) + subroutine surfq(psinq,qpsi,nq,qval,psival) use const_and_precisions, only : wp_ use magsurf_data, only : npoints - use interp_eqprof, only : psinr=>psiq,qpsi=>q,nq use utils, only : locate, intlin implicit none c arguments + integer, intent(in) :: nq + real(wp_), dimension(nq), intent(in) :: psinq,qpsi real(wp_) :: qval,psival c local variables integer :: ncnt,i1,ipr @@ -1639,7 +1376,7 @@ c c locate psi surface for q=qval c call locate(abs(qpsi),nq,qval,i1) - call intlin(abs(qpsi(i1)),psinr(i1),abs(qpsi(i1+1)),psinr(i1+1), + call intlin(abs(qpsi(i1)),psinq(i1),abs(qpsi(i1+1)),psinq(i1+1), . qval,psival) ipr=1 call contours_psi(psival,rcn,zcn,ipr) @@ -1651,7 +1388,7 @@ c c subroutine bfield_res(rv,zv,nr,nz) use const_and_precisions, only : wp_ - use interp_eqprof, only : bfield + use equilibrium, only : bfield implicit none c arguments integer, intent(in) :: nr, nz @@ -1705,215 +1442,6 @@ c 113 format(i6,12(1x,e12.5)) end c -c - subroutine rhotor(rhotsx) - use const_and_precisions, only : wp_ - use interp_eqprof, only : nq,psinr=>psiq,qpsi=>q,crhot,cq,rhopq - use simplespline, only : difcs - implicit none -c arguments - real(wp_), intent(out) :: rhotsx -c local variables - integer :: iopt,ier,k - real(wp_) :: dx,drhot - real(wp_), dimension(nq) :: rhotnr -c - if(allocated(cq)) deallocate(cq) - if(allocated(crhot)) deallocate(crhot) - if(allocated(rhopq)) deallocate(rhopq) - allocate(cq(nq,4),crhot(nq,4),rhopq(nq)) -c -c normalized toroidal rho : ~ Integral q(psi) dpsi -c - iopt=0 - call difcs(psinr,abs(qpsi),nq,iopt,cq,ier) -c - rhotnr(1)=0.0_wp_ - do k=1,nq-1 - dx=psinr(k+1)-psinr(k) - drhot=dx*(cq(k,1)+dx*(cq(k,2)/2.0_wp_+dx*(cq(k,3)/3.0_wp_ - . +dx*cq(k,4)/4.0_wp_))) - rhotnr(k+1)=rhotnr(k)+drhot - end do - rhotsx=rhotnr(nq) - rhotnr=sqrt(rhotnr/rhotsx) - rhotsx=sign(rhotsx,qpsi(nq)) -c -c spline interpolation of rhotor -c - iopt=0 - rhopq=sqrt(psinr) - call difcs(rhopq,rhotnr,nq,iopt,crhot,ier) - return - end -c -c -c - function fq_eq(psi) - use const_and_precisions, only : wp_ - use interp_eqprof, only : psinr=>psiq,nq,cq - use simplespline, only :spli - implicit none -c arguments - real(wp_) :: psi,fq_eq -c local variables - integer :: irt - real(wp_) :: dps -c - irt=int((nq-1)*psi+1) - if(irt.eq.0) irt=1 - if(irt.eq.nq) irt=nq-1 - dps=psi-psinr(irt) - fq_eq=spli(cq,nq,irt,dps) -c - return - end -c -c -c - function frhotor_eq(rhop) - use const_and_precisions, only : wp_ - use interp_eqprof, only : rhopq,nq,crhot - use simplespline, only : spli - use utils, only : locate - implicit none -c arguments - real(wp_) :: rhop,frhotor_eq -c local variables - integer :: irt - real(wp_) :: drh -c - call locate(rhopq,nq,rhop,irt) - irt=min(max(1,irt),nq-1) - drh=rhop-rhopq(irt) - frhotor_eq=spli(crhot,nq,irt,drh) -c - return - end -c -c -c - function frhotor(rhop) - use const_and_precisions, only : wp_ - use graydata_flags, only : iequil - implicit none -c arguments - real(wp_) :: rhop,frhotor -c common/external functions/variables - real(wp_) :: frhotor_eq,frhotor_an -c - if(iequil.eq.1) then - frhotor=frhotor_an(rhop) - else - frhotor=frhotor_eq(rhop) - end if -c - return - end -c -c -c - function frhotor_av(rhop) - use const_and_precisions, only : wp_ - use magsurf_data, only : rpstab, crhotq, npsi - use simplespline, only : spli - use utils, only : locate - implicit none -c arguments - real(wp_) :: rhop,frhotor_av -c local variables - integer :: ip - real(wp_) :: drh -c - call locate(rpstab,npsi,rhop,ip) - ip=min(max(1,ip),npsi-1) - drh=rhop-rpstab(ip) - frhotor_av=spli(crhotq,npsi,ip,drh) -c - return - end -c -c -c - subroutine rhopol - use const_and_precisions, only : wp_ - use dierckx, only : curfit,splev - implicit none -c local constants - integer, parameter :: nnr=101,nrest=nnr+4, - . lwrkp=nnr*4+nrest*16 -c local variables - integer :: i,iopt,ier,kspl - integer, dimension(nrest) :: iwrkp - real(wp_) :: dr,psin,xb,xe,ss,rp - real(wp_), dimension(nnr) :: rhop,rhot,rhopi - real(wp_), dimension(lwrkp) :: wrkp - real(wp_), dimension(nrest) :: wp -c common/external functions/variables - integer :: nsrp - real(wp_) :: frhotor - real(wp_), dimension(nrest) :: trp,crp -c - common/coffrtp/trp - common/coffrn/nsrp - common/coffrp/crp -c - dr=1.0_wp_/dble(nnr-1) - do i=1,nnr - rhop(i)=(i-1)*dr - psin=rhop(i)*rhop(i) - rhot(i)=frhotor(rhop(i)) - wp(i)=1.0_wp_ - end do - wp(1)=1.0e3_wp_ - wp(nnr)=1.0e3_wp_ - -c spline interpolation of rhopol versus rhotor - iopt=0 - xb=0.0_wp_ - xe=1.0_wp_ - ss=0.00001_wp_ - kspl=3 - call curfit(iopt,nnr,rhot,rhop,wp,xb,xe,kspl,ss,nrest,nsrp, - . trp,crp,rp,wrkp,lwrkp,iwrkp,ier) -c print*,ier - call splev(trp,nsrp,crp,3,rhot,rhopi,nnr,ier) -c do i=1,nnr -c write(54,*) rhop(i),rhot(i),rhopi(i) -c end do -c - return - end -c -c -c - function frhopol(rhot) - use dierckx, only : splev - use const_and_precisions, only : wp_ - implicit none -c local constants - integer, parameter :: nnr=101,nrest=nnr+4 -c arguments - real(wp_) :: rhot,frhopol -c local variables - integer :: ier - real(wp_), dimension(1) :: rrs,ffspl -c common/external functions/variables - integer :: nsrp - real(wp_), dimension(nrest) :: trp,crp -c - common/coffrtp/trp - common/coffrn/nsrp - common/coffrp/crp -c - rrs(1)=rhot - call splev(trp,nsrp,crp,3,rrs,ffspl,1,ier) - frhopol=ffspl(1) - return -c - end -c -c c subroutine cniteq(rqgrid,zqgrid,btotal,nreq,nzeq,h, . ncon,npts,icount,rcon,zcon) @@ -2102,8 +1630,8 @@ c use const_and_precisions, only : wp_,pi use graydata_par, only : rwallm use magsurf_data, only : npoints - use interp_eqprof, only : psiant,psinop,nsr,nsz, - . cc=>cceq,tr,tz,rup,zup,rlw,zlw,kspl + use equilibrium, only : psiant,psinop,nsr,nsz, + . cc=>cceq,tr,tz,rup,zup,rlw,zlw,kspl,points_tgo use dierckx, only : profil,sproota implicit none c local constants @@ -2170,8 +1698,8 @@ c subroutine flux_average use const_and_precisions, only : wp_,zero,one,pi,ccj=>mu0inv use magsurf_data - use interp_eqprof, only : btrcen,btaxis,rmaxis,zmaxis,phitedge, - . equinum_fpol,bfield + use equilibrium, only : btrcen,btaxis,rmaxis,zmaxis,phitedge, + . equinum_fpol,bfield,fq,frhotor use simplespline, only : difcs use dierckx, only : regrid,coeff_parder implicit none @@ -2187,7 +1715,7 @@ c local variables . ratio_cdbtor,ratio_pltor,qq,fc,rhot2q,height,height2,r2iav,currp, . area,volume,ajphiav,bbav,bmmx,bmmn,btot0,bpoloid0,rpsim0,dla,dlb, . dlp,drc,ph,area2,rzp,rz,rpsim,zpsim,btot,bpoloid,dlph,ajphi0, - . shlam,srl,rl2,rl0,rl,dhlam,dhlam0,ccfh,s,frhotor_eq,fq_eq,ajphi, + . shlam,srl,rl2,rl0,rl,dhlam,dhlam0,ccfh,s,ajphi, . bphi,brr,bzz,riav,fp real(wp_), dimension(nnintp) :: dadrhotv,dvdrhotv,vratjpl real(wp_), dimension(2*ncnt) :: dlpv @@ -2400,9 +1928,9 @@ c ratio_pltor = Jcd_||/J_phi Jcd_|| = qq=abs(dvdpsi*fpolv*r2iav/(4*pi*pi)) qqv(jp)=qq - dadrhotv(jp)=phitedge*frhotor_eq(height)/fq_eq(height2) + dadrhotv(jp)=phitedge*frhotor(height)/fq(height2) . *dadpsi/pi - dvdrhotv(jp)=phitedge*frhotor_eq(height)/fq_eq(height2) + dvdrhotv(jp)=phitedge*frhotor(height)/fq(height2) . *dvdpsi/pi c computation of rhot from calculated q profile @@ -2467,7 +1995,7 @@ c ffhlam = Bmn/Bmx/fc integral_lambda^1 dlam/ rpstab(jp)=1.0_wp_ pstab(jp)=1.0_wp_ end if - rhot_eq(jp)=frhotor_eq(sqrt(pstab(jp))) + rhot_eq(jp)=frhotor(sqrt(pstab(jp))) write(56,99) pstab(jp),rhot_eq(jp),rhotqv(jp), . bav(jp),bmxpsi(jp),bmnpsi(jp),varea(jp),vvol(jp), . vcurrp(jp),vajphiav(jp),qqv(jp),ffc(jp), @@ -2570,7 +2098,7 @@ c use const_and_precisions, only : wp_,zero,one,pi,ccj=>mu0inv use graydata_anequil, only : q0,qa,alq,b0,rr0m,zr0m,rpam use magsurf_data - use interp_eqprof, only : btrcen + use equilibrium, only : btrcen, frhotor use simplespline, only : difcs use dierckx, only : regrid,coeff_parder implicit none @@ -2601,7 +2129,6 @@ c local variables !to accept optional arguments c common/external functions/variables real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,dfpv - real(wp_) :: frhotor_an c common/fpol/fpolv,dfpv common/derip1/dpsidr,dpsidz @@ -2818,7 +2345,7 @@ c ratio_pltor = Jcd_||/J_phi Jcd_|| = qq=abs(dvdpsi*fpolv*r2iav/(4.0_wp_*pi*pi)) qqv(jp)=qq - rn=frhotor_an(sqrt(pstab(jp))) + rn=frhotor(sqrt(pstab(jp))) qqan=q0+(qa-q0)*rn**alq dadr=2.0_wp_*pi*rn*rpam*rpam @@ -2895,7 +2422,7 @@ c ffhlam = Bmn/Bmx/fc integral_lambda^1 dlam/ rpstab(jp)=1.0_wp_ pstab(jp)=1.0_wp_ end if - rhot_eq(jp)=frhotor_an(sqrt(pstab(jp))) + rhot_eq(jp)=frhotor(sqrt(pstab(jp))) write(56,99) pstab(jp),rhot_eq(jp),rhotqv(jp), . bav(jp),bmxpsi(jp),bmnpsi(jp),varea(jp),vvol(jp), . vcurrp(jp),vajphiav(jp),qqv(jp),ffc(jp), @@ -2956,7 +2483,7 @@ c use dierckx, only : curfit,splev use graydata_par, only : sgniphi use graydata_anequil, only : q0,qa,alq,b0,rr0m,rpam - use interp_eqprof, only : psia + use equilibrium, only : psia implicit none c local constants integer, parameter :: nnr=101,nrest=nnr+4,lwrk=nnr*4+nrest*16 @@ -3029,38 +2556,12 @@ c spline interpolation of rhotor versus rhopol end c c -c - function frhotor_an(rhop) - use const_and_precisions, only : wp_ - use dierckx, only : splev - use utils, only : locate - implicit none -c local contants - integer, parameter :: nnr=101,nrest=nnr+4 -c arguments - real(wp_) :: rhop,frhotor_an -c local variables - integer :: nsrot,ier - real(wp_), dimension(1) :: rrs(1),ffspl(1) - real(wp_), dimension(nrest) :: trot,crot -c - common/coffrptt/trot - common/coffrnt/nsrot - common/coffrpt/crot -c - rrs(1)=rhop - call splev(trot,nsrot,crot,3,rrs,ffspl,1,ier) - frhotor_an=ffspl(1) -c - return - end -c -c c subroutine contours_psi_an(h,rcn,zcn,ipr) use const_and_precisions, only : wp_,pi use graydata_anequil, only : rr0m,zr0m,rpam use magsurf_data, only : npoints + use equilibrium, only : frhotor implicit none c arguments integer :: ipr @@ -3069,12 +2570,10 @@ c arguments c local variables integer :: np,ic real(wp_) :: rn,th -c common/external functions/variables - real(wp_) :: frhotor_an c np=(npoints-1)/2 th=pi/dble(np) - rn=frhotor_an(sqrt(h)) + rn=frhotor(sqrt(h)) do ic=1,npoints zcn(ic)=zr0m+rpam*rn*sin(th*(ic-1)) rcn(ic)=rr0m+rpam*rn*cos(th*(ic-1)) @@ -3774,7 +3273,7 @@ c use const_and_precisions, only : wp_,pi,ccj=>mu0inv use graydata_flags, only : iequil use graydata_par, only : sgnbphi - use interp_eqprof, only : equinum_fpol, equinum_psi + use equilibrium, only : equinum_fpol, equinum_psi implicit none c arguments real(wp_) :: xx,yy,zz @@ -3951,7 +3450,7 @@ c use const_and_precisions, only : wp_ use graydata_par, only : sgnbphi,sgniphi use graydata_anequil, only : b0,rr0m,zr0m,rpam,q0,qa,alq - use interp_eqprof, only : psia + use equilibrium, only : psia,frhopol implicit none c arguments real(wp_) :: rrm,zzm @@ -3961,7 +3460,6 @@ c local variables c common/external functions/variables real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv, . dfpv - real(wp_) :: frhopol c common/derip1/dpsidr,dpsidz common/derip2/ddpsidrr,ddpsidzz,ddpsidrz @@ -4019,7 +3517,7 @@ c c subroutine tor_curr_psi(h,ajphi) use const_and_precisions, only : wp_ - use interp_eqprof, only : zmaxis + use equilibrium, only : zmaxis implicit none c arguments real(wp_) :: h,ajphi @@ -4036,7 +3534,7 @@ c c subroutine tor_curr(rpsim,zpsim,ajphi) use const_and_precisions, only : wp_,ccj=>mu0inv - use interp_eqprof, only : equinum_psi + use equilibrium, only : equinum_psi implicit none c arguments real(wp_) :: rpsim,zpsim,ajphi @@ -4058,7 +3556,7 @@ c c subroutine psi_raxis(h,r1,r2) use const_and_precisions, only : wp_ - use interp_eqprof, only : psiant,psinop,zmaxis,nsr, + use equilibrium, only : psiant,psinop,zmaxis,nsr, . nsz,cc=>cceq,tr,tz,kspl use dierckx, only : profil,sproota implicit none @@ -4088,7 +3586,7 @@ c c subroutine sub_xg_derxg(psinv) use const_and_precisions, only : wp_ - use interp_eqprof, only : psia + use equilibrium, only : psia use coreprofiles, only : density implicit none c arguments @@ -5089,6 +4587,7 @@ c . dvol,darea,ratjav,ratjbv,ratjplv) use const_and_precisions, only : wp_,zero,one,izero use graydata_flags, only : nnd + use equilibrium, only : frhotor,frhopol implicit none integer :: it,ipec real(wp_) :: drt,rt,rt1,rhop1 @@ -5099,7 +4598,6 @@ c real(wp_), dimension(0:nnd), intent(out) :: rtabpsi1 real(wp_), dimension(nnd), intent(out) :: dvol,darea real(wp_), dimension(nnd), intent(out) :: ratjav,ratjbv,ratjplv - real(wp_), external :: frhotor,frhopol c ipec positive build equidistant grid dimension nnd c ipec negative read input grid @@ -5330,6 +4828,7 @@ c local variables c radial average values over power and current density profile use const_and_precisions, only : wp_,one,zero,izero,pi use graydata_flags, only : nnd,ieccd + use equilibrium, only : frhopol implicit none real(wp_), intent(in) :: pabs,currt real(wp_), intent(out) :: rhotpav,rhotjava @@ -5340,7 +4839,7 @@ c radial average values over power and current density profile real(wp_) :: ratjamx,ratjbmx,ratjplmx real(wp_), dimension(nnd), intent(in) :: rhot_tab,dvol,darea real(wp_), dimension(nnd), intent(in) :: dpdv,ajphiv - real(wp_), external :: frhopol,fdadrhot,fdvdrhot + real(wp_), external :: fdadrhot,fdvdrhot real(wp_) :: sccsa real(wp_) :: rhotjav,rhot2pav,rhot2java @@ -5639,7 +5138,7 @@ c use const_and_precisions, only : wp_ use graydata_flags, only : iequil use coreprofiles, only : psdbnd - use interp_eqprof, only : zbmin,zbmax,equinum_psi + use equilibrium, only : zbinf,zbsup,equinum_psi implicit none c arguments real(wp_) :: rrm,zzm @@ -5653,7 +5152,7 @@ c end if c if (psinv.ge.0.0_wp_.and.psinv.lt.psdbnd) then - if (psinv.lt.1.0_wp_.and.zzm.lt.zbmin.or.zzm.gt.zbmax) then + if (psinv.lt.1.0_wp_.and.zzm.lt.zbinf.or.zzm.gt.zbsup) then inside_plasma=.false. else inside_plasma=.true. diff --git a/src/minpack.f90 b/src/minpack.f90 index e661639..6229b34 100644 --- a/src/minpack.f90 +++ b/src/minpack.f90 @@ -12,8 +12,8 @@ contains integer, intent(in) :: n, ldfjac, lwa integer, intent(out) :: info real(wp_), intent(in) :: tol - real(wp_), intent(out) :: fvec(n), fjac(ldfjac,n), wa(lwa) - real(wp_), intent(inout) :: x(n) + real(wp_), intent(out) :: wa(lwa) + real(wp_), intent(inout) :: fvec(n), fjac(ldfjac,n), x(n) ! ********** ! ! subroutine hybrj1 @@ -116,8 +116,9 @@ contains subroutine fcn(n,x,fvec,fjac,ldfjac,iflag) use const_and_precisions, only : wp_ implicit none - integer :: n,ldfjac,iflag - real(wp_) :: x(n),fvec(n),fjac(ldfjac,n) + integer, intent(in) :: n,ldfjac,iflag + real(wp_), intent(in) :: x(n) + real(wp_), intent(inout) :: fvec(n),fjac(ldfjac,n) end subroutine fcn end interface @@ -313,8 +314,9 @@ contains subroutine fcn(n,x,fvec,fjac,ldfjac,iflag) use const_and_precisions, only : wp_ implicit none - integer :: n,ldfjac,iflag - real(wp_) :: x(n),fvec(n),fjac(ldfjac,n) + integer, intent(in) :: n,ldfjac,iflag + real(wp_), intent(in) :: x(n) + real(wp_), intent(inout) :: fvec(n),fjac(ldfjac,n) end subroutine fcn end interface ! @@ -585,6 +587,588 @@ contains if (nprint > 0) call fcn(n,x,fvec,fjac,ldfjac,iflag) end subroutine hybrj + subroutine hybrj1mv(fcn,n,x,f0,fvec,fjac,ldfjac,tol,info,wa,lwa) + use const_and_precisions, only : zero, one + implicit none +! arguments + integer, intent(in) :: n, ldfjac, lwa + integer, intent(out) :: info + real(wp_), intent(in) :: tol,f0(n) + real(wp_), intent(out) :: wa(lwa) + real(wp_), intent(inout) :: fvec(n), fjac(ldfjac,n), x(n) +! ********** +! +! subroutine hybrj1mv +! +! the purpose of hybrj1mv is to find a zero of a system of +! n nonlinear functions in n variables by a modification +! of the powell hybrid method. this is done by using the +! more general nonlinear equation solver hybrjmv. the user +! must provide a subroutine which calculates the functions +! and the jacobian. +! +! the subroutine statement is +! +! subroutine hybrj1mv(fcn,n,x,f0,fvec,fjac,ldfjac,tol,info,wa,lwa) +! +! where +! +! fcn is the name of the user-supplied subroutine which +! calculates the functions and the jacobian. fcn must +! be declared in an external statement in the user +! calling program, and should be written as follows. +! +! subroutine fcn(n,x,f0,fvec,fjac,ldfjac,iflag) +! integer n,ldfjac,iflag +! real(8) x(n),fvec(n),fjac(ldfjac,n) +! ---------- +! if iflag = 1 calculate the functions at x and +! return this vector in fvec. do not alter fjac. +! if iflag = 2 calculate the jacobian at x and +! return this matrix in fjac. do not alter fvec. +! --------- +! return +! end +! +! the value of iflag should not be changed by fcn unless +! the user wants to terminate execution of hybrj1mv. +! in this case set iflag to a negative integer. +! +! n is a positive integer input variable set to the number +! of functions and variables. +! +! x is an array of length n. on input x must contain +! an initial estimate of the solution vector. on output x +! contains the final estimate of the solution vector. +! +! fvec is an output array of length n which contains +! the functions evaluated at the output x. +! +! fjac is an output n by n array which contains the +! orthogonal matrix q produced by the qr factorization +! of the final approximate jacobian. +! +! ldfjac is a positive integer input variable not less than n +! which specifies the leading dimension of the array fjac. +! +! tol is a nonnegative input variable. termination occurs +! when the algorithm estimates that the relative error +! between x and the solution is at most tol. +! +! info is an integer output variable. if the user has +! terminated execution, info is set to the (negative) +! value of iflag. see description of fcn. otherwise, +! info is set as follows. +! +! info = 0 improper input parameters. +! +! info = 1 algorithm estimates that the relative error +! between x and the solution is at most tol. +! +! info = 2 number of calls to fcn with iflag = 1 has +! reached 100*(n+1). +! +! info = 3 tol is too small. no further improvement in +! the approximate solution x is possible. +! +! info = 4 iteration is not making good progress. +! +! wa is a work array of length lwa. +! +! lwa is a positive integer input variable not less than +! (n*(n+13))/2. +! +! subprograms called +! +! user-supplied ...... fcn +! +! minpack-supplied ... hybrjmv +! +! argonne national laboratory. minpack project. march 1980. +! burton s. garbow, kenneth e. hillstrom, jorge j. more +! +! ********** +! local variables + integer :: j, lr, maxfev, mode, nfev, njev, nprint + real(wp_) :: xtol +! parameters + real(wp_), parameter :: factor=1.0e2_wp_ + + interface + subroutine fcn(n,x,f0,fvec,fjac,ldfjac,iflag) + use const_and_precisions, only : wp_ + implicit none + integer, intent(in) :: n,ldfjac,iflag + real(wp_), intent(in) :: x(n),f0(n) + real(wp_), intent(inout) :: fvec(n),fjac(ldfjac,n) + end subroutine fcn + end interface + + info = 0 +! +! check the input parameters for errors. +! + if (n <= 0 .or. ldfjac < n .or. tol < zero & + .or. lwa < (n*(n + 13))/2) return +! +! call hybrjmv. +! + maxfev = 100*(n + 1) + xtol = tol + mode = 2 + do j = 1, n + wa(j) = one + end do + nprint = 0 + lr = (n*(n + 1))/2 + call hybrjmv(fcn,n,x,f0,fvec,fjac,ldfjac,xtol,maxfev,wa(1),mode, & + factor,nprint,info,nfev,njev,wa(6*n+1),lr,wa(n+1), & + wa(2*n+1),wa(3*n+1),wa(4*n+1),wa(5*n+1)) + if (info == 5) info = 4 + end subroutine hybrj1mv + + subroutine hybrjmv(fcn,n,x,f0,fvec,fjac,ldfjac,xtol,maxfev,diag,mode, & + factor,nprint,info,nfev,njev,r,lr,qtf,wa1,wa2, & + wa3,wa4) + use const_and_precisions, only : zero, one, epsmch=>comp_eps + implicit none +! arguments + integer, intent(in) :: n, ldfjac, maxfev, mode, nprint, lr + integer, intent(out) :: info, nfev, njev + real(wp_), intent(in) :: xtol, factor, f0(n) + real(wp_), intent(out) :: fvec(n), fjac(ldfjac,n), r(lr), qtf(n), & + wa1(n), wa2(n), wa3(n), wa4(n) + real(wp_), intent(inout) :: x(n), diag(n) +! ********** +! +! subroutine hybrj +! +! the purpose of hybrj is to find a zero of a system of +! n nonlinear functions in n variables by a modification +! of the powell hybrid method. the user must provide a +! subroutine which calculates the functions and the jacobian. +! +! the subroutine statement is +! +! subroutine hybrj(fcn,n,x,f0,fvec,fjac,ldfjac,xtol,maxfev,diag, +! mode,factor,nprint,info,nfev,njev,r,lr,qtf, +! wa1,wa2,wa3,wa4) +! +! where +! +! fcn is the name of the user-supplied subroutine which +! calculates the functions and the jacobian. fcn must +! be declared in an external statement in the user +! calling program, and should be written as follows. +! +! subroutine fcn(n,x,f0,fvec,fjac,ldfjac,iflag) +! integer n,ldfjac,iflag +! real(8) x(n),f0(n),fvec(n),fjac(ldfjac,n) +! ---------- +! if iflag = 1 calculate the functions at x and +! return this vector in fvec. do not alter fjac. +! if iflag = 2 calculate the jacobian at x and +! return this matrix in fjac. do not alter fvec. +! --------- +! return +! end +! +! the value of iflag should not be changed by fcn unless +! the user wants to terminate execution of hybrj. +! in this case set iflag to a negative integer. +! +! n is a positive integer input variable set to the number +! of functions and variables. +! +! x is an array of length n. on input x must contain +! an initial estimate of the solution vector. on output x +! contains the final estimate of the solution vector. +! +! fvec is an output array of length n which contains +! the functions evaluated at the output x. +! +! fjac is an output n by n array which contains the +! orthogonal matrix q produced by the qr factorization +! of the final approximate jacobian. +! +! ldfjac is a positive integer input variable not less than n +! which specifies the leading dimension of the array fjac. +! +! xtol is a nonnegative input variable. termination +! occurs when the relative error between two consecutive +! iterates is at most xtol. +! +! maxfev is a positive integer input variable. termination +! occurs when the number of calls to fcn with iflag = 1 +! has reached maxfev. +! +! diag is an array of length n. if mode = 1 (see +! below), diag is internally set. if mode = 2, diag +! must contain positive entries that serve as +! multiplicative scale factors for the variables. +! +! mode is an integer input variable. if mode = 1, the +! variables will be scaled internally. if mode = 2, +! the scaling is specified by the input diag. other +! values of mode are equivalent to mode = 1. +! +! factor is a positive input variable used in determining the +! initial step bound. this bound is set to the product of +! factor and the euclidean norm of diag*x if nonzero, or else +! to factor itself. in most cases factor should lie in the +! interval (.1,100.). 100. is a generally recommended value. +! +! nprint is an integer input variable that enables controlled +! printing of iterates if it is positive. in this case, +! fcn is called with iflag = 0 at the beginning of the first +! iteration and every nprint iterations thereafter and +! immediately prior to return, with x and fvec available +! for printing. fvec and fjac should not be altered. +! if nprint is not positive, no special calls of fcn +! with iflag = 0 are made. +! +! info is an integer output variable. if the user has +! terminated execution, info is set to the (negative) +! value of iflag. see description of fcn. otherwise, +! info is set as follows. +! +! info = 0 improper input parameters. +! +! info = 1 relative error between two consecutive iterates +! is at most xtol. +! +! info = 2 number of calls to fcn with iflag = 1 has +! reached maxfev. +! +! info = 3 xtol is too small. no further improvement in +! the approximate solution x is possible. +! +! info = 4 iteration is not making good progress, as +! measured by the improvement from the last +! five jacobian evaluations. +! +! info = 5 iteration is not making good progress, as +! measured by the improvement from the last +! ten iterations. +! +! nfev is an integer output variable set to the number of +! calls to fcn with iflag = 1. +! +! njev is an integer output variable set to the number of +! calls to fcn with iflag = 2. +! +! r is an output array of length lr which contains the +! upper triangular matrix produced by the qr factorization +! of the final approximate jacobian, stored rowwise. +! +! lr is a positive integer input variable not less than +! (n*(n+1))/2. +! +! qtf is an output array of length n which contains +! the vector (q transpose)*fvec. +! +! wa1, wa2, wa3, and wa4 are work arrays of length n. +! +! subprograms called +! +! user-supplied ...... fcn +! +! minpack-supplied ... dogleg,enorm, +! qform,qrfac,r1mpyq,r1updt +! +! fortran-supplied ... abs,dmax1,dmin1,mod +! +! argonne national laboratory. minpack project. march 1980. +! burton s. garbow, kenneth e. hillstrom, jorge j. more +! +! ********** +! local variables + integer :: i, iflag, iter, j, jm1, l, ncfail, ncsuc, nslow1, nslow2 + integer, dimension(1) :: iwa + logical :: jeval, sing + real(wp_) :: actred, delta, fnorm, fnorm1, pnorm, prered, & + ratio, summ, temp, xnorm +! parameters + real(wp_), parameter :: p1 = 1.0e-1_wp_, p5 = 5.0e-1_wp_, & + p001 = 1.0e-3_wp_, p0001 = 1.0e-4_wp_ + + interface + subroutine fcn(n,x,f0,fvec,fjac,ldfjac,iflag) + use const_and_precisions, only : wp_ + implicit none + integer, intent(in) :: n,ldfjac,iflag + real(wp_), intent(in) :: x(n),f0(n) + real(wp_), intent(inout) :: fvec(n),fjac(ldfjac,n) + end subroutine fcn + end interface +! + info = 0 + iflag = 0 + nfev = 0 + njev = 0 +! +! check the input parameters for errors. +! + if (n <= 0 .or. ldfjac < n .or. xtol < zero & + .or. maxfev <= 0 .or. factor <= zero & + .or. lr < (n*(n + 1))/2) go to 300 + if (mode == 2) then + do j = 1, n + if (diag(j) <= zero) go to 300 + end do + end if +! +! evaluate the function at the starting point +! and calculate its norm. +! + iflag = 1 + call fcn(n,x,f0,fvec,fjac,ldfjac,iflag) + nfev = 1 + if (iflag < 0) go to 300 + fnorm = enorm(n,fvec) +! +! initialize iteration counter and monitors. +! + iter = 1 + ncsuc = 0 + ncfail = 0 + nslow1 = 0 + nslow2 = 0 +! +! beginning of the outer loop. +! + do + jeval = .true. +! +! calculate the jacobian matrix. +! + iflag = 2 + call fcn(n,x,f0,fvec,fjac,ldfjac,iflag) + njev = njev + 1 + if (iflag < 0) go to 300 +! +! compute the qr factorization of the jacobian. +! + call qrfac(n,n,fjac,ldfjac,.false.,iwa,1,wa1,wa2,wa3) +! +! on the first iteration and if mode is 1, scale according +! to the norms of the columns of the initial jacobian. +! + if (iter == 1) then + if (mode /= 2) then + do j = 1, n + diag(j) = wa2(j) + if (wa2(j) == zero) diag(j) = one + end do + end if +! +! on the first iteration, calculate the norm of the scaled x +! and initialize the step bound delta. +! + do j = 1, n + wa3(j) = diag(j)*x(j) + end do + xnorm = enorm(n,wa3) + delta = factor*xnorm + if (delta == zero) delta = factor + end if +! +! form (q transpose)*fvec and store in qtf. +! + do i = 1, n + qtf(i) = fvec(i) + end do + do j = 1, n + if (fjac(j,j) /= zero) then + summ = zero + do i = j, n + summ = summ + fjac(i,j)*qtf(i) + end do + temp = -summ/fjac(j,j) + do i = j, n + qtf(i) = qtf(i) + fjac(i,j)*temp + end do + end if + end do +! +! copy the triangular factor of the qr factorization into r. +! + sing = .false. + do j = 1, n + l = j + jm1 = j - 1 + do i = 1, jm1 + r(l) = fjac(i,j) + l = l + n - i + end do + r(l) = wa1(j) + if (wa1(j) == zero) sing = .true. + end do +! +! accumulate the orthogonal factor in fjac. +! + call qform(n,n,fjac,ldfjac,wa1) +! +! rescale if necessary. +! + if (mode /= 2) then + do j = 1, n + diag(j) = dmax1(diag(j),wa2(j)) + end do + end if +! +! beginning of the inner loop. +! + do +! +! if requested, call fcn to enable printing of iterates. +! + if (nprint > 0) then + iflag = 0 + if (mod(iter-1,nprint) == 0) call fcn(n,x,f0,fvec,fjac,ldfjac,iflag) + if (iflag < 0) go to 300 + end if +! +! determine the direction p. +! + call dogleg(n,r,lr,diag,qtf,delta,wa1,wa2,wa3) +! +! store the direction p and x + p. calculate the norm of p. +! + do j = 1, n + wa1(j) = -wa1(j) + wa2(j) = x(j) + wa1(j) + wa3(j) = diag(j)*wa1(j) + end do + pnorm = enorm(n,wa3) +! +! on the first iteration, adjust the initial step bound. +! + if (iter == 1) delta = dmin1(delta,pnorm) +! +! evaluate the function at x + p and calculate its norm. +! + iflag = 1 + call fcn(n,wa2,f0,wa4,fjac,ldfjac,iflag) + nfev = nfev + 1 + if (iflag < 0) go to 300 + fnorm1 = enorm(n,wa4) +! +! compute the scaled actual reduction. +! + actred = -one + if (fnorm1 < fnorm) actred = one - (fnorm1/fnorm)**2 +! +! compute the scaled predicted reduction. +! + l = 1 + do i = 1, n + summ = zero + do j = i, n + summ = summ + r(l)*wa1(j) + l = l + 1 + end do + wa3(i) = qtf(i) + summ + end do + temp = enorm(n,wa3) + prered = zero + if (temp < fnorm) prered = one - (temp/fnorm)**2 +! +! compute the ratio of the actual to the predicted +! reduction. +! + ratio = zero + if (prered > zero) ratio = actred/prered +! +! update the step bound. +! + if (ratio < p1) then + ncsuc = 0 + ncfail = ncfail + 1 + delta = p5*delta + else + ncfail = 0 + ncsuc = ncsuc + 1 + if (ratio >= p5 .or. ncsuc > 1) delta = dmax1(delta,pnorm/p5) + if (abs(ratio-one) <= p1) delta = pnorm/p5 + end if +! +! test for successful iteration. +! + if (ratio >= p0001) then +! +! successful iteration. update x, fvec, and their norms. +! + do j = 1, n + x(j) = wa2(j) + wa2(j) = diag(j)*x(j) + fvec(j) = wa4(j) + end do + xnorm = enorm(n,wa2) + fnorm = fnorm1 + iter = iter + 1 + end if +! +! determine the progress of the iteration. +! + nslow1 = nslow1 + 1 + if (actred >= p001) nslow1 = 0 + if (jeval) nslow2 = nslow2 + 1 + if (actred >= p1) nslow2 = 0 +! +! test for convergence. +! + if (delta <= xtol*xnorm .or. fnorm == zero) info = 1 + if (info /= 0) go to 300 +! +! tests for termination and stringent tolerances. +! + if (nfev >= maxfev) info = 2 + if (p1*dmax1(p1*delta,pnorm) <= epsmch*xnorm) info = 3 + if (nslow2 == 5) info = 4 + if (nslow1 == 10) info = 5 + if (info /= 0) go to 300 +! +! criterion for recalculating jacobian. +! + if (ncfail == 2) exit +! +! calculate the rank one modification to the jacobian +! and update qtf if necessary. +! + do j = 1, n + summ = zero + do i = 1, n + summ = summ + fjac(i,j)*wa4(i) + end do + wa2(j) = (summ - wa3(j))/pnorm + wa1(j) = diag(j)*((diag(j)*wa1(j))/pnorm) + if (ratio >= p0001) qtf(j) = summ + end do +! +! compute the qr factorization of the updated jacobian. +! + call r1updt(n,n,r,lr,wa1,wa2,wa3,sing) + call r1mpyq(n,n,fjac,ldfjac,wa2,wa3) + call r1mpyq(1,n,qtf,1,wa2,wa3) +! +! end of the inner loop. +! + jeval = .false. + end do +! +! end of the outer loop. +! + end do + 300 continue +! +! termination, either normal or user imposed. +! + if (iflag < 0) info = iflag + iflag = 0 + if (nprint > 0) call fcn(n,x,f0,fvec,fjac,ldfjac,iflag) + end subroutine hybrjmv + subroutine dogleg(n,r,lr,diag,qtb,delta,x,wa1,wa2) use const_and_precisions, only : zero, one, epsmch=>comp_eps implicit none