istp/gray
2
0
Fork 0
Code Issues Pull Requests Packages Projects Releases 1 Wiki Activity

inclusion of updated routines for analytical equilibrium/profiles; further cleaning of equidata

Browse Source
This commit is contained in:
Lorenzo Figini 2015-10-12 17:09:28 +00:00
parent 79c080b14d
commit af5fb208b2
5 changed files with 1290 additions and 684 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
Makefile
View File

@ -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

81
src/coreprofiles.f90
View File

@ -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

590
src/interp_eqprof.f90 → src/equilibrium.f90
View File

@ -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

701
src/gray.f
View File

File diff suppressed because it is too large Load Diff

596
src/minpack.f90
View File

@ -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