diff --git a/src/dispersion.f90 b/src/dispersion.f90 index 233fc6a..3898575 100644 --- a/src/dispersion.f90 +++ b/src/dispersion.f90 @@ -6,7 +6,7 @@ module dispersion integer, parameter :: npts=500 real(wp_), parameter :: tmax=5.0_wp_,dtex=2.0_wp_*tmax/dble(npts) ! variables - real(wp_), save :: ttv(npts+1),extv(npts+1) + real(wp_), dimension(npts+1), save :: ttv,extv ! contains ! diff --git a/src/gray.f b/src/gray.f index 8debb5e..53abe33 100644 --- a/src/gray.f +++ b/src/gray.f @@ -1,16 +1,22 @@ + use const_and_precisions, only : wp_ use graydata_flags, only : ipass,igrad - - implicit real*8 (a-h,o-z) - + implicit none +c local variables + real(wp_) :: p0mw1 +c common/external functions/variables + integer :: ierr,index_rt + real(wp_) :: sox,p0mw,powrfl,taumn,taumx,pabstot,currtot, + . pabstott,currtott +c common/ierr/ierr common/mode/sox common/p0/p0mw common/powrfl/powrfl common/index_rt/index_rt common/taumnx/taumn,taumx,pabstot,currtot - +c c read data plus initialization - +c index_rt=1 call prfile call paraminit @@ -22,22 +28,21 @@ c read data plus initialization print*,' IERR = ', ierr stop end if - +c c beam/ray propagation call gray_integration - +c c postprocessing - call after_gray_integration - +c pabstott=pabstot currtott=currtot - +c if (ipass.gt.1) then c second pass into plasma p0mw1=p0mw igrad=0 - +c index_rt=2 p0mw=p0mw1*powrfl call prfile @@ -48,10 +53,10 @@ c second pass into plasma call after_gray_integration pabstott=pabstott+pabstot currtott=currtott+currtot - +c index_rt=3 sox=-sox - p0mw=p0mw1*(1.0d0-powrfl) + p0mw=p0mw1*(1.0_wp_-powrfl) call prfile call vectinit2 call paraminit @@ -62,59 +67,65 @@ c second pass into plasma currtott=currtott+currtot end if print*,' ' - write(6,*) 'Pabs (MW), Icd (kA) = ', pabstott,currtott*1.0d+3 - + write(6,*) 'Pabs (MW), Icd (kA) = ', pabstott,currtott*1.0e3_wp_ +c end - - +c +c +c subroutine gray_integration + use const_and_precisions, only : wp_ use graydata_flags, only : dst - implicit none - integer istep,istop,index_rt - integer nstep - real*8 st,strfl11 - integer i - real*8 st0 - +c local variables + integer :: i + real(wp_) :: st0 +c common/external functions/variables + integer :: istep,nstep,istop,index_rt + real(wp_) :: st,strfl11 +c common/ss/st common/istep/istep common/nstep/nstep common/istop/istop common/strfl11/strfl11 common/index_rt/index_rt - +c c ray integration: begin - st0=0.0d0 + st0=0.0_wp_ if(index_rt.gt.1) st0=strfl11 do i=1,nstep istep=i st=i*dst+st0 - +c c advance one step - call rkint4 - +c c calculations after one step: - call after_onestep(istep,istop) if(istop.eq.1) exit c end do - +c c ray integration: end - +c return end - - +c +c +c subroutine after_gray_integration - use const_and_precisions, only : zero + use const_and_precisions, only : wp_,zero use graydata_flags, only : ibeam,iwarm,ilarm,iequil,iprof, . filenmeqq,filenmprf,filenmbm use graydata_anequil, only : dens0,te0 - - implicit real*8 (a-h,o-z) + implicit none +c local variables + integer :: iproj,nfilp + real(wp_) :: currtka,pabs,currt +c common/external functions/variables + integer :: nrayr,nrayth,index_rt + real(wp_) :: st,taumn,taumx,pabstot,currtot c common/ss/st common/nray/nrayr,nrayth @@ -140,7 +151,7 @@ c end if c write(*,'(a,f9.4)') 'Pabs_tot (MW) = ',pabstot - currtka =currtot*1.0d3 + currtka =currtot*1.0e3_wp_ write(*,'(a,f9.4)') 'I_tot (kA) = ',currtka c if (index_rt.eq.1) then @@ -164,34 +175,47 @@ c c c subroutine after_onestep(i,istop) - use const_and_precisions, only : pi + use const_and_precisions, only : wp_,pi use reflections, only : inside use graydata_flags, only : iwarm,istpr0,istpl0,dst,ipass,igrad use graydata_par, only : psipol0,chipol0 use interp_eqprof, only : zbmin,zbmax,rlim,zlim,nlim - - implicit real*8 (a-h,o-z) - complex*16 ext,eyt,extr,eytr,exin2,eyin2 - parameter(jmx=31,kmx=36,nmx=8000) - parameter(taucr=12.0d0) - dimension psjki(jmx,kmx,nmx),ppabs(jmx,kmx,nmx),ccci(jmx,kmx,nmx) - dimension iiv(jmx,kmx),tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx) - dimension iop(jmx,kmx),iow(jmx,kmx),ihcd(jmx,kmx),istore(jmx,kmx) - dimension tau1v(jmx,kmx),yyrfl(jmx,kmx,6),y(6),dery(6) - dimension ext(jmx,kmx,0:4),eyt(jmx,kmx,0:4) - dimension xv(3),anv(3),xvrfl(3),anvrfl(3),xvvac(3),anw(3),anwcl(3) - dimension xwcl(3),xvjk(3,jmx,kmx),anvjk(3,jmx,kmx) - dimension ywrk(6,jmx,kmx),ypwrk(6,jmx,kmx) - logical ins_pl,inside_plasma + implicit none +c local constants + integer, parameter :: jmx=31,kmx=36,nmx=8000 + real(wp_), parameter :: taucr=12.0_wp_ +c arguments + integer :: i, istop +c local variables + integer iproj,nfilp,irfl,j,k,kkk,iopmin,iowmax,iowmin,ivac,j1,k1, + . kkkk + real(wp_) :: qqout,uuout,vvout,qqin2,uuin2,vvin2,rrm,zzm,anvjkr, + . aknmin,akdotn,anwclr + real(wp_), dimension(6) :: y,dery + real(wp_), dimension(3) :: xvrfl,anvrfl,xvvac,anw + real(wp_), dimension(3,jmx,kmx) :: xvjk,anvjk + complex(wp_) :: extr,eytr,exin2,eyin2 + logical :: ins_pl +c common/external functions/variables + integer :: nrayr,nrayth,istpr,istpl,jclosest,ierr,index_rt + integer, dimension(jmx,kmx) :: iiv,iop,iow,ihcd,istore + real(wp_) :: psinv,psinv11,taumn,taumx,pabstot,currtot,psipol, + . chipol,powrfl,strfl11 + real(wp_), dimension(3) :: anwcl,xwcl,xv,anv + real(wp_), dimension(jmx,kmx) :: tau1v + real(wp_), dimension(jmx,kmx,6) :: yyrfl + real(wp_), dimension(6,jmx,kmx) :: ywrk,ypwrk + real(wp_), dimension(jmx,kmx,nmx) :: ppabs,ccci,tauv,alphav,psjki + complex(wp_), dimension(jmx,kmx,0:4) :: ext,eyt + logical :: inside_plasma +c external inside_plasma c common/pcjki/ppabs,ccci common/atjki/tauv,alphav common/tau1v/tau1v -c common/nray/nrayr,nrayth common/istgr/istpr,istpl -c common/iiv/iiv common/iov/iop,iow,ihcd,istore common/refln/anwcl,xwcl,jclosest @@ -202,21 +226,19 @@ c common/taumnx/taumn,taumx,pabstot,currtot common/xv/xv common/anv/anv -c common/polcof/psipol,chipol common/evt/ext,eyt common/yyrfl/yyrfl common/powrfl/powrfl common/strfl11/strfl11 common/index_rt/index_rt -c common/wrk/ywrk,ypwrk c - pabstot=0.0d0 - currtot=0.0d0 - taumn=1d+30 - taumx=-1d+30 - psinv11=1.0d0 + pabstot=0.0_wp_ + currtot=0.0_wp_ + taumn=1.0e+30_wp_ + taumx=-1.0e+30_wp_ + psinv11=1.0_wp_ iopmin=100 iowmin=100 iowmax=0 @@ -244,11 +266,11 @@ c exit end if psjki(j,k,i)=psinv - rrm=sqrt(xv(1)**2+xv(2)**2)/100.d0 - zzm=xv(3)/100.d0 + rrm=sqrt(xv(1)**2+xv(2)**2)/100.0_wp_ + zzm=xv(3)/100.0_wp_ if (iwarm.gt.0.and.i.gt.1.and.ihcd(j,k).gt.0) then - if(psinv.ge.0.and.psinv.le.1.0d0.and. + if(psinv.ge.0.and.psinv.le.1.0_wp_.and. . zzm.ge.zbmin.and.zzm.le.zbmax) then call pabs_curr(i,j,k) iiv(j,k)=i @@ -306,8 +328,8 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc end if xv=xvrfl anv=anvrfl - rrm=1.d-2*sqrt(xv(1)**2+xv(2)**2) - zzm=1.d-2*xv(3) + rrm=1.0e-2_wp_*sqrt(xv(1)**2+xv(2)**2) + zzm=1.0e-2_wp_*xv(3) ywrk(1:3,j,k)=xv ywrk(4:6,j,k)=anv igrad=0 @@ -332,21 +354,21 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc end do if(jclosest.le.nrayr) then - aknmin=1.0d0 + aknmin=1.0_wp_ do j=1,nrayr kkk=nrayth if(j.eq.1) kkk=1 do k=1,kkk anwclr=(anwcl(1)*xwcl(1)+anwcl(2)*xwcl(2)) - . /sqrt(xwcl(1)**2+xwcl(2)**2) + . /sqrt(xwcl(1)**2+xwcl(2)**2) anvjkr=(anvjk(1,j,k)*xvjk(1,j,k)+anvjk(2,j,k)*xvjk(2,j,k)) - . /sqrt(xvjk(1,j,k)**2+xvjk(2,j,k)**2) + . /sqrt(xvjk(1,j,k)**2+xvjk(2,j,k)**2) akdotn=anwclr*anvjkr+anwcl(3)*anvjk(3,j,k) if(akdotn.lt.aknmin) aknmin=akdotn end do end do else - aknmin=-1.0d0 + aknmin=-1.0_wp_ end if cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc @@ -380,18 +402,18 @@ c reflection point, or when no reflection has occurred and c central ray re-enters the plasma if((ipass.eq.1 .and. ((iopmin.gt.1) .or. - . (taumn.lt.1d+30.and.taumn.gt.taucr))) - . .or.(index_rt.gt.1 .and. (iopmin.gt.1 .or. iowmin.gt.1 .or. - . (taumn.lt.1d+30.and.taumn.gt.taucr)))) then + . (taumn.lt.1.0e+30_wp_.and.taumn.gt.taucr))) + . .or.(index_rt.gt.1 .and. (iopmin.gt.1 .or. iowmin.gt.1 .or. + . (taumn.lt.1.0e+30_wp_.and.taumn.gt.taucr)))) then istop=1 else if(ipass.gt.1 .and. index_rt.eq.1 .and. . ((iowmin.gt.1 .and. aknmin.gt.0) .or. . (iowmax.le.1 .and. iop(1,1).gt.2))) then c flag second pass mode coupling as unset - powrfl=-1.0d0 - qqout=0.0d0 - uuout=0.0d0 - vvout=0.0d0 + powrfl=-1.0_wp_ + qqout=0.0_wp_ + uuout=0.0_wp_ + vvout=0.0_wp_ do j=1,nrayr kkk=nrayth if(j.eq.1) kkk=1 @@ -404,7 +426,7 @@ c store missing initial conditions for the second pass tau1v(j,k)=tauv(j,k,iiv(j,k)) end if c determine mode coupling at the plasma boundary - if (powrfl.lt.0.0d0) then + if (powrfl.lt.0.0_wp_) then call vacuum_rt(xvjk(:,j,k),anvjk(:,j,k),xvvac,ivac) c look for first ray hitting the plasma, starting from the central c and evaluate polarization @@ -429,10 +451,10 @@ c if found, use its polarization state to compute mode coupling end do powloop c if no ray completed a first pass in the plasma, use central ray c initial polarization (possibly reflected) - if (qqout.le.0.0d0) then + if (qqout.le.0.0_wp_) then call stokes(ext(1,1,0),eyt(1,1,0),qqout,uuout,vvout) end if - powrfl=0.5d0*(1.0d0+vvout*vvin2+ + powrfl=0.5_wp_*(1.0_wp_+vvout*vvin2+ . uuout*uuin2+qqout*qqin2) end if end if @@ -451,21 +473,30 @@ c c c subroutine print_output(i,j,k) - use const_and_precisions, only : pi + use const_and_precisions, only : wp_,pi use graydata_flags, only : iequil,istpl0 use graydata_par, only : psdbnd - - implicit real*8 (a-h,o-z) - parameter(ndim=6,jmx=31,kmx=36,nmx=8000) - parameter(taucr=12.0d0) - dimension ywrk(ndim,jmx,kmx),ypwrk(ndim,jmx,kmx) - dimension psjki(jmx,kmx,nmx) - dimension tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx) - dimension pdjki(jmx,kmx,nmx),ppabs(jmx,kmx,nmx) - dimension currj(jmx,kmx,nmx),didst(jmx,kmx,nmx),ccci(jmx,kmx,nmx) - dimension q(jmx),tau1v(jmx,kmx) - complex*16 ex,ey,ez - + implicit none +c local constants + integer, parameter :: ndim=6,jmx=31,kmx=36,nmx=8000 + real(wp_), parameter :: taucr=12.0_wp_ +c arguments + integer i,j,k +c local variables + real(wp_) :: x,y,z,rr,rrm,zzm,stm,xxm,yym,anx,any,anz,anr,anphi, + . phi,phideg,psi,rhot,bbr,bbz,bpol,dens11,dids11,dpdv11,ajcd11, + . alpha11,taujk,tau11,pt11 +c common/external functions/variables + integer :: nharm,nhf,iohkawa,index_rt,nrayr,nrayth,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 + real(wp_), dimension(jmx) :: q + real(wp_), dimension(jmx,kmx) :: tau1v + real(wp_), dimension(ndim,jmx,kmx) :: ywrk,ypwrk + real(wp_), dimension(jmx,kmx,nmx) :: psjki,tauv,alphav,pdjki, + . ppabs,currj,didst,ccci + complex(wp_) :: ex,ey,ez c common/psjki/psjki common/atjki/tauv,alphav @@ -478,11 +509,9 @@ c common/tau1v/tau1v common/qw/q common/index_rt/index_rt - common/ss/st common/nray/nrayr,nrayth common/istgr/istpr,istpl -c common/parpl/brr,bphi,bzz,ajphi common/btot/btot common/xgxg/xg @@ -491,11 +520,9 @@ c common/tete/tekev common/absor/alpha,effjcd,akim,tau0 common/epolar/ex,ey,ez -c common/nplr/anpl,anpr common/ddd/ddr,an2s,an2,fdia,bdotgr,ddi,ddr11 common/nprw/anprr,anpri -c common/wrk/ywrk,ypwrk c x=ywrk(1,j,k) @@ -510,11 +537,11 @@ c anphi=(any*x-anx*y)/rr c cnphi=(any*x-anx*y) c - rrm=rr*1.0d-2 - zzm=z*1.0d-2 - stm=st*1.0d-2 - xxm=x*1.0d-2 - yym=y*1.0d-2 + rrm=rr*1.0e-2_wp_ + zzm=z*1.0e-2_wp_ + stm=st*1.0e-2_wp_ + xxm=x*1.0e-2_wp_ + yym=y*1.0e-2_wp_ c if(index_rt.gt.1) then taujk=tauv(j,k,i)+tau1v(j,k) @@ -525,31 +552,31 @@ c c central ray only begin if(j.eq.1) then phi=acos(x/sqrt(y*y+x*x)) - if(y.lt.0.0d0) phi=-phi - phideg=phi*180.0d0/pi + if(y.lt.0.0_wp_) phi=-phi + phideg=phi*180.0_wp_/pi psi=psjki(j,k,i) - rhot=1.0d0 - bbr=0.0d0 - bbz=0.0d0 - bpol=0.0d0 - rhot=1.0d0 - dens11=0.0d0 - if(psi.ge.0.0d0) then - if (psi.le.1.0d0) rhot=frhotor(psi) + rhot=1.0_wp_ + bbr=0.0_wp_ + bbz=0.0_wp_ + bpol=0.0_wp_ + rhot=1.0_wp_ + dens11=0.0_wp_ + if(psi.ge.0.0_wp_) then + if (psi.le.1.0_wp_) rhot=frhotor(psi) bbr=brr bbz=bzz bpol=sqrt(bbr**2+bbz**2) else - tekev=0.0d0 - akim=0.0d0 + tekev=0.0_wp_ + akim=0.0_wp_ end if ddr11=ddr c cutoff=xg-(1-yg)*(1-anpl**2) c print dIds in A/m/W, Jphi and Jt in MA/m2, ki in m^-1 - if(psi.le.psdbnd.and.psi.ge.0.0d0) dens11=dens - dids11=didst(j,k,i)*1.0d2/(p0mw*q(j)) + if(psi.le.psdbnd.and.psi.ge.0.0_wp_) dens11=dens + dids11=didst(j,k,i)*1.0e2_wp_/(p0mw*q(j)) dpdv11=pdjki(j,k,i)/(p0mw*q(j)) ajcd11=currj(j,k,i)/(p0mw*q(j)) alpha11=alphav(j,k,i) @@ -558,7 +585,7 @@ c print dIds in A/m/W, Jphi and Jt in MA/m2, ki in m^-1 write(4,99) stm,rrm,zzm,phideg, . psi,rhot,dens11,tekev, - . bbr,bphi,bbz,ajphi*1.0d-6,sqrt(an2),anpl,akim*100, + . bbr,bphi,bbz,ajphi*1.0e-6_wp_,sqrt(an2),anpl,akim*100, . alpha11,tau11,pt11,ajcd11,dids11, . dble(nhf),dble(iohkawa),dble(index_rt) @@ -587,71 +614,76 @@ c c subroutine prfile implicit none - integer*4 index_rt +c common/external functions/variables + integer :: index_rt +c common/index_rt/index_rt - - If(index_rt.eq.1) then - - write(4,*)' #sst R z phi psi rhot ne Te BR Bphi Bz Jphi '// - .'N Npl ki alpha tau Pt Jcd dIds nh iohkw index_rt' - write(8,*) ' #istep j k xt yt zt rt psin' - write(9,*) ' #istep j k xt yt zt rt psin' - write(17,*) ' #sst Dr_11 Dr_Nr1 Di_Nr1' - write(33,*) ' #i j k sst x y R z psi tauv Npl alpha index_rt' - write(12,*) ' #i sst psi w1 w2' - write(7,*)'#Icd Pa Jphip dPdVp '// - .'rhotj rhotjava rhotp rhotpav drhotjava drhotpav ratjamx '// - .'ratjbmx stmx psipol chipol index_rt Jphimx dPdVmx drhotj drhotp' - write(48,*) '#psi rhot Jphi Jcdb dPdV Icdins Pins P% index_rt' - write(66,*) "# psipol0 chipol0 powrfl" - +c + if(index_rt.eq.1) then + write(4,*)' #sst R z phi psi rhot ne Te BR Bphi Bz Jphi '// + . 'N Npl ki alpha tau Pt Jcd dIds nh iohkw index_rt' + write(8,*) ' #istep j k xt yt zt rt psin' + write(9,*) ' #istep j k xt yt zt rt psin' + write(17,*) ' #sst Dr_11 Dr_Nr1 Di_Nr1' + write(33,*) ' #i j k sst x y R z psi tauv Npl alpha index_rt' + write(12,*) ' #i sst psi w1 w2' + write(7,*)'#Icd Pa Jphip dPdVp '// + . 'rhotj rhotjava rhotp rhotpav drhotjava drhotpav ratjamx '// + . 'ratjbmx stmx psipol chipol index_rt Jphimx dPdVmx drhotj '// + . 'drhotp' + write(48,*) '#psi rhot Jphi Jcdb dPdV Icdins Pins P% index_rt' + write(66,*) "# psipol0 chipol0 powrfl" +c else - - -c If(index_rt.eq.3) then - write(4,*) ' ' - write(8,*) ' ' - write(9,*) ' ' - write(17,*) ' ' - write(12,*) ' ' - write(48,*) ' ' - +c +c if(index_rt.eq.3) then + write(4,*) ' ' + write(8,*) ' ' + write(9,*) ' ' + write(17,*) ' ' + write(12,*) ' ' + write(48,*) ' ' c end if end if +c return end c c c subroutine read_data - use const_and_precisions, only : qe=>ecgs_,me=>mecgs_,vc=>ccgs_, - * cvdr=>degree,pi + use const_and_precisions, only : wp_,qe=>ecgs_,me=>mecgs_, + . vc=>ccgs_,cvdr=>degree,pi use green_func_p, only:Setup_SpitzFunc use graydata_flags use graydata_par use graydata_anequil use interp_eqprof, only : rmxm,rlim,zlim,nlim,zbmin,zbmax, . btrcen,rcen,alloc_lim - - implicit real*8 (a-h,o-z) - character*8 wdat - character*10 wtim - parameter(nmx=8000) + implicit none +c local constants + integer, parameter :: nmx=8000 +c local variables + integer :: nfil,iox,ierr,leqq,lprf + real(wp_) :: dummy,bresg,r00,anr0c,anphi0c,fghz,zeff, + . w0csi,w0eta,d0csi,d0eta,phi0,zrcsi,zreta + character(len=8) :: wdat + character(len=10) :: wtim +c common/external functions/variables + integer :: nstep,nrayr,nrayth + real(wp_) :: xgcn,ak0,akinv,fhz,wcsi,weta,rcicsi,rcieta,phiw, + . phir,anx0c,any0c,anz0c,x00,y00,z00,bres,p0mw,sox,alpha0,beta0 + integer :: length c common/xgcn/xgcn - common/nstep/nstep -c common/nray/nrayr,nrayth -c common/parwv/ak0,akinv,fhz common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir common/anic/anx0c,any0c,anz0c common/mirr/x00,y00,z00 -c common/parbres/bres common/p0/p0mw -c common/mode/sox common/angles/alpha0,beta0 c @@ -752,13 +784,13 @@ c iscal ne Te scaling 0: nustar=const, 1: n_greenw=const; 2 no rescaling c factT factn factor for Te&ne scaling c read(2,*) factb, iscal,factt,factn - if (factb.eq.0.0d0) factb=1.0d0 + if (factb.eq.0.0_wp_) factb=1.0_wp_ c c psbnd value of psi ( > 1 ) of density boundary c read(2,*) filenmprf read(2,*) psdbnd - if(iprof.eq.0) psdbnd=1.0d0 + if(iprof.eq.0) psdbnd=1.0_wp_ c c signum of toroidal B and I c icocos index for equilibrium from COCOS - O. Sauter Feb 2012 @@ -783,13 +815,13 @@ c if iequil=0,1 read(2,*) rr0,zr0,rpa read(2,*) b0 read(2,*) q0,qa,alq - rr0m=rr0/1.0d2 - zr0m=zr0/1.0d2 - rpam=rpa/1.0d2 + rr0m=rr0/1.0e2_wp_ + zr0m=zr0/1.0e2_wp_ + rpam=rpa/1.0e2_wp_ rcen=rr0m btrcen=b0 - zbmin=(zr0-rpa)/100.d0 - zbmax=(zr0+rpa)/100.d0 + zbmin=(zr0-rpa)/100.0_wp_ + zbmax=(zr0+rpa)/100.0_wp_ b0=b0*factb call flux_average_an c call print_prof_an @@ -808,30 +840,30 @@ c print*,' ' end if c - fhz=fghz*1.0d9 - ak0=2.0d9*pi*fghz/vc - akinv=1.0d0/ak0 + fhz=fghz*1.0e9_wp_ + ak0=2.0e9_wp_*pi*fghz/vc + akinv=1.0_wp_/ak0 c - bresg=2.0d0*pi*fhz*me*vc/qe - bres=bresg*1.0d-4 + bresg=2.0_wp_*pi*fhz*me*vc/qe + bres=bresg*1.0e-4_wp_ c c xg=xgcn*dens19 c - xgcn=1.0d-5*qe**2/(pi*me*fghz**2) + xgcn=1.0e-5_wp_*qe**2/(pi*me*fghz**2) c - sox=-1.0d0 - if(iox.eq.2) sox=1.0d0 + sox=-1.0_wp_ + if(iox.eq.2) sox=1.0_wp_ c c read data for beam from file if ibeam>0 c if(ibeam.gt.0) then call read_beams else - zrcsi=0.5d0*ak0*w0csi**2 - zreta=0.5d0*ak0*w0eta**2 + zrcsi=0.5_wp_*ak0*w0csi**2 + zreta=0.5_wp_*ak0*w0eta**2 if(igrad.gt.0) then - wcsi=w0csi*sqrt(1.0d0+(d0csi/zrcsi)**2) - weta=w0eta*sqrt(1.0d0+(d0eta/zreta)**2) + wcsi=w0csi*sqrt(1.0_wp_+(d0csi/zrcsi)**2) + weta=w0eta*sqrt(1.0_wp_+(d0eta/zreta)**2) rcicsi=-d0csi/(d0csi**2+zrcsi**2) rcieta=-d0eta/(d0eta**2+zreta**2) phir=phiw @@ -911,13 +943,13 @@ c set simple limiter as two cylindrical walls at rwallm and r00 call alloc_lim(ierr) if (ierr.ne.0) stop rlim(1)=rwallm - rlim(2)=max(r00*1.d-2,rmxm) + rlim(2)=max(r00*1.0e-2_wp_,rmxm) rlim(3)=rlim(2) rlim(4)=rlim(1) rlim(5)=rlim(1) - zlim(1)=(z00-dst*nmx)*1.d-2 + zlim(1)=(z00-dst*nmx)*1.0e-2_wp_ zlim(2)=zlim(1) - zlim(3)=(z00+dst*nmx)*1.d-2 + zlim(3)=(z00+dst*nmx)*1.0e-2_wp_ zlim(4)=zlim(3) zlim(5)=zlim(1) ipass=abs(ipass) @@ -981,14 +1013,19 @@ c c c subroutine surf_anal - use const_and_precisions, only : pi + use const_and_precisions, only : wp_,pi use graydata_anequil, only : b0,rr0m,zr0m,rpam use magsurf_data, only : npsi,npoints,psicon,rcon,zcon, . alloc_surf_anal - - implicit real*8(a-h,o-z) + implicit none +c local variables + integer :: i,k,ierr + real(wp_) :: rni,rres,zzres,zmx,zmn,dal,drn,drrm,dzzm,rrm,zzm,zrm +c common/external functions/variables + real(wp_) :: bres +c common/parbres/bres - +c npsi=10 npoints=101 call alloc_surf_anal(ierr) @@ -997,8 +1034,8 @@ c c print circular magnetic surfaces iequil=1 c write(71,*) '#i psin R z' - dal=2.0d-2*pi - drn=0.1d0 + dal=2.0e-2_wp_*pi + drn=0.1_wp_ do i=1,npsi rni=i*drn psicon(i)=rni @@ -1022,28 +1059,35 @@ c zmx=zr0m+rpam zmn=zr0m-rpam do i=1,21 - zzres=zmn+(i-1)*(zmx-zmn)/2.0d1 + zzres=zmn+(i-1)*(zmx-zmn)/2.0e1_wp_ write(70,111) i,bres,rres,zzres end do return 111 format(i4,20(1x,e16.8e3)) end - +c +c +c subroutine bres_anal - implicit real*8(a-h,o-z) - parameter(pi=3.14159265358979d0) - common/parban/b0,rr0m,zr0m,rpam + use const_and_precisions, only : wp_,pi + use graydata_anequil, only : b0,rr0m,zr0m,rpam + implicit none +c local variables + integer :: i + real(wp_) :: rres,zmx,zmn,zzres +c common/external functions/variables + real(wp_) :: bres +c common/parbres/bres c c print resonant B iequil=1 -c write(70,*)'#i Btot R z' rres=b0*rr0m/bres zmx=zr0m+rpam zmn=zr0m-rpam do i=1,21 - zzres=zmn+(i-1)*(zmx-zmn)/2.0d1 + zzres=zmn+(i-1)*(zmx-zmn)/2.0e1_wp_ write(70,111) i,bres,rres,zzres end do @@ -1051,24 +1095,28 @@ c 111 format(i4,20(1x,e16.8e3)) end c +c c subroutine read_beams + use const_and_precisions, only : wp_ use graydata_flags, only : filenmbm, ibeam use utils, only : locate use simplespline, only : difcs,spli - - implicit real*8(a-h,o-z) - parameter(nstrmx=50) -c - dimension alphastv(nstrmx),betastv(nstrmx),cbeta(nstrmx,4) - dimension x00v(nstrmx),y00v(nstrmx),z00v(nstrmx) - dimension cx0(nstrmx,4),cy0(nstrmx,4),cz0(nstrmx,4) - dimension waist1v(nstrmx),waist2v(nstrmx) - dimension rci1v(nstrmx),rci2v(nstrmx) - dimension cwaist1(nstrmx,4),cwaist2(nstrmx,4) - dimension crci1(nstrmx,4),crci2(nstrmx,4) - dimension phi1v(nstrmx),phi2v(nstrmx) - dimension cphi1(nstrmx,4),cphi2(nstrmx,4) + implicit none +c local constants + integer, parameter :: nstrmx=50 +c local variables + integer :: i,k,ii,nfbeam,lbm,nisteer,steer,iopt,ier + real(wp_) :: alphast,betast,x00mm,y00mm,z00mm,waist1,waist2, + . dvir1,dvir2,delta,phi1,phi2,zr1,zr2,w1,w2,rci1,rci2,dal,betst + real(wp_), dimension(nstrmx) :: alphastv,betastv,x00v,y00v, + . z00v,waist1v,waist2v,rci1v,rci2v,phi1v,phi2v + real(wp_), dimension(nstrmx,4) :: cbeta,cx0,cy0,cz0,cwaist1, + . cwaist2,crci1,crci2,cphi1,cphi2 +c common/external functions/variables + real(wp_) :: wcsi,weta,rcicsi,rcieta,phiw,phir,x00,y00,z00, + . alpha0,beta0,ak0,akinv,fhz + integer :: length c common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir common/mirr/x00,y00,z00 @@ -1093,25 +1141,25 @@ c . waist1,dvir1,waist2,dvir2,delta phi1=delta phi2=delta - zr1=0.5d-1*ak0*waist1**2 - zr2=0.5d-1*ak0*waist2**2 - w1=waist1*sqrt(1.0d0+(dvir1/zr1)**2) - w2=waist2*sqrt(1.0d0+(dvir2/zr2)**2) + zr1=0.5e-1_wp_*ak0*waist1**2 + zr2=0.5e-1_wp_*ak0*waist2**2 + w1=waist1*sqrt(1.0_wp_+(dvir1/zr1)**2) + w2=waist2*sqrt(1.0_wp_+(dvir2/zr2)**2) rci1=-dvir1/(dvir1**2+zr1**2) rci2=-dvir2/(dvir2**2+zr2**2) else read(nfbeam,*) steer,alphast,betast,x00mm,y00mm,z00mm, . w1,w2,rci1,rci2,phi1,phi2 end if - x00v(i)=0.1d0*x00mm - y00v(i)=0.1d0*y00mm - z00v(i)=0.1d0*z00mm + x00v(i)=0.1_wp_*x00mm + y00v(i)=0.1_wp_*y00mm + z00v(i)=0.1_wp_*z00mm alphastv(i)=alphast betastv(i)=betast - waist1v(i)=0.1d0*w1 - rci1v(i)=1.0d1*rci1 - waist2v(i)=0.1d0*w2 - rci2v(i)=1.0d1*rci2 + waist1v(i)=0.1_wp_*w1 + rci1v(i)=1.0e1_wp_*rci1 + waist2v(i)=0.1_wp_*w2 + rci2v(i)=1.0e1_wp_*rci2 phi1v(i)=phi1 phi2v(i)=phi2 end do @@ -1166,7 +1214,7 @@ c c c subroutine equidata - use const_and_precisions, only : pi + use const_and_precisions, only : wp_,pi use utils, only : vmaxmini use dierckx, only : curfit,splev,regrid,bispev,coeff_parder use graydata_flags, only : ipsinorm,sspl,ixp,icocos,neqdsk @@ -1177,15 +1225,24 @@ c . rlw,zlw,tr,tz,tfp,cc=>cceq,cfp,cc01=>cceq01,cc10=>cceq10, . cc20=>cceq20,cc02=>cceq02,cc11=>cceq11,psinr,qpsi,rv,zv,psin, . psi,rbbbs,zbbbs,nbbbs,alloc_equilvec,alloc_bnd - - implicit real*8 (a-h,o-z) -c parameter(np=100) - parameter(kspl=3) - character*48 stringa -c dimension rcon(2*np+1),zcon(2*np+1) + implicit none +c local constants + integer, parameter :: kspl=3 +c local variables + integer :: idum,nrest,nzest,lwrk,liwrk,lw10,lw01,lw20,lw02,lw11, + . lwrkf,lwrkbsp,liwrkbsp,i,j,ij,iopt,ier,nsr,nsz,nur,nuz,info, + . iqmn,iqmx,icocosmod,ierr integer, dimension(:), allocatable :: iwrk,iwrkf,iwrkbsp - real*8, dimension(:), allocatable :: fpoli,fvpsi,ffvpsi,fpol, + real(wp_) :: xdum,drnr1,dznz1,rleft,zmid,psiaxis,psiedge,current, + . dpsinr,psia0,fp,xb,xe,ssfp,psinxp,rmop,zmop,psinoptmp,r10,z10, + . rxp,zxp,psinxptmp,r1,z1,psin1,rbmin,rbmax,q2,q15,qmin,qmax, + . rhot15,psi15,rhot2,psi2 + real(wp_), dimension(:), allocatable :: fpoli,fvpsi,ffvpsi,fpol, . wrkf,wf,wrk,wrkbsp,pres,ffprim,pprim + character(len=48) :: stringa +c common/external functions/variables + real(wp_) :: rhotsx + real(wp_) :: frhotor c common/rhotsx/rhotsx c @@ -1215,7 +1272,7 @@ c lwrkf=nr*4+nrest*16 lwrkbsp=4*(nr+nz) liwrkbsp=nr+nz - +c if(allocated(fvpsi)) deallocate(fvpsi) if(allocated(pres)) deallocate(pres) if(allocated(ffprim)) deallocate(ffprim) @@ -1230,7 +1287,7 @@ c if(allocated(iwrk)) deallocate(iwrk) if(allocated(wrkbsp)) deallocate(wrkbsp) if(allocated(iwrkbsp)) deallocate(iwrkbsp) - +c allocate(fvpsi(nr*nz),pres(nr),ffprim(nr),pprim(nr), . ffvpsi(nr*nz),fpol(nr),fpoli(nr),wrkf(lwrkf),iwrkf(nrest), . wf(nr),wrk(lwrk),iwrk(liwrk),wrkbsp(lwrkbsp),iwrkbsp(liwrkbsp)) @@ -1259,12 +1316,11 @@ c read (neqdsk,*) (qpsi(i),i=1,nr) end if 2020 format (5e16.9) - c c compensate for different reference systems c icocosmod=mod(icocos,10) - +c if (mod(icocos,2).eq.0.and.icocosmod.gt.0) then c icocos mod 2 = 0: toroidal angle phi CW (opposite to gray convention) btrcen=-btrcen @@ -1288,21 +1344,19 @@ c icocos mod 10 = 2,3,6,7: psi decreasing with CCW Ip end do end if end if - c c add check for Ip/psi and B0/Fpol sign consistency? c current=sign(current,psiaxis-psiedge) btrcen=sign(btrcen,fpol(nr)) - c c length in m !!! c dr=drnr1/dble(nr-1) dz=dznz1/dble(nz-1) rv(1)=rleft - zv(1)=zmid-dznz1/2.0d0 - dpsinr=1.0d0/dble(nr-1) + zv(1)=zmid-dznz1/2.0_wp_ + dpsinr=1.0_wp_/dble(nr-1) c do i=1,nr psinr(i)=(i-1)*dpsinr @@ -1317,18 +1371,18 @@ c rmxm=rv(nr) zmnm=zv(1) zmxm=zv(nz) - +c c psi function - +c psia0=psiedge-psiaxis c icocos=0: adapt psi to force Ip sign, otherwise maintain psi - if (icocosmod.ne.0) sgniphi=sign(1.0d0,-psia0) + if (icocosmod.ne.0) sgniphi=sign(1.0_wp_,-psia0) current=sign(current,sgniphi) - +c psia=-sgniphi*abs(psia0)*factb c icocos>10: input psi is in Wb c icocos<10: input psi is in Wb/rad (gray convention) - if (icocos.ge.10) psia=psia/(2.0d0*pi) + if (icocos.ge.10) psia=psia/(2.0_wp_*pi) c c do j=1,nz c do i=1,nr @@ -1336,7 +1390,7 @@ c write(80,2021) rv(i),zv(j),psi(i,j) c enddo c write(80,*) ' ' c enddo - +c do j=1,nz do i=1,nr if(ipsinorm.eq.0) then @@ -1358,14 +1412,14 @@ c . wrk,lwrk,iwrk,liwrk,ier) c if ier=-1 data are fitted using sspl=0 if(ier.eq.-1) then - sspl=0.0d0 + sspl=0.0_wp_ call regrid(iopt,nr,rv,nz,zv,fvpsi,rmnm,rmxm,zmnm,zmxm, . kspl,kspl,sspl,nrest,nzest,nsr,tr,nsz,tz,cc,fp, . wrk,lwrk,iwrk,liwrk,ier) end if nsrt=nsr nszt=nsz - if (sspl.gt.0.0d0) then + if (sspl.gt.0.0_wp_) then call bispev(tr,nsr,tz,nsz,cc,kspl,kspl,rv,nr,zv,nz,ffvpsi, . wrkbsp,lwrkbsp,iwrkbsp,liwrkbsp,ier) c @@ -1409,21 +1463,21 @@ c c scaling of f_poloidal c c icocos=0: adapt fpol to force Ip sign, otherwise maintain fpol - if (icocosmod.ne.0) sgnbphi=sign(1.0d0,fpol(nr)) + if (icocosmod.ne.0) sgnbphi=sign(1.0_wp_,fpol(nr)) btrcen=sign(btrcen,sgnbphi) - +c do i=1,nr fpol(i)=sgnbphi*abs(fpol(i))*factb - wf(i)=1.0d0 + wf(i)=1.0_wp_ end do - wf(nr)=1.0d2 + wf(nr)=1.0e2_wp_ c c spline interpolation of fpol(i) c iopt=0 - xb=0.0d0 - xe=1.0d0 - ssfp=0.01d0 + xb=0.0_wp_ + xe=1.0_wp_ + ssfp=0.01_wp_ call curfit(iopt,nr,psinr,fpol,wf,xb,xe,kspl,ssfp,nrest,nsft, . tfp,cfp,fp,wrkf,lwrkf,iwrkf,ier) c @@ -1433,24 +1487,24 @@ c c read plasma boundaries from eqdsk file c read (neqdsk,*) nbbbs,nlim - +c call alloc_bnd(ierr) if (ierr.ne.0) stop - +c if(nbbbs.gt.0) then if(ipsinorm.eq.1) - . read (neqdsk,*) (rbbbs(i),zbbbs(i),i=1,nbbbs) + . read(neqdsk,*) (rbbbs(i),zbbbs(i),i=1,nbbbs) if(ipsinorm.eq.0) - . read (neqdsk,2020) (rbbbs(i),zbbbs(i),i=1,nbbbs) + . read(neqdsk,2020) (rbbbs(i),zbbbs(i),i=1,nbbbs) c do i=1,nbbbs c write(51,*) rbbbs(i),zbbbs(i) c end do end if if(nlim.gt.0) then if(ipsinorm.eq.1) - . read (neqdsk,*) (rlim(i),zlim(i),i=1,nlim) - if(ipsinorm.eq.0) - . read (neqdsk,2020) (rlim(i),zlim(i),i=1,nlim) + . read(neqdsk,*) (rlim(i),zlim(i),i=1,nlim) + if(ipsinorm.eq.0) + . read(neqdsk,2020) (rlim(i),zlim(i),i=1,nlim) end if c c compute max and min z of last closed surface @@ -1458,8 +1512,8 @@ c rbmin=rmaxis rbmax=rmaxis if (nbbbs.gt.1) then - zbmin=1.0d+30 - zbmax=-1.0d+30 + zbmin=1.0e+30_wp_ + zbmax=-1.0e+30_wp_ do i=1,nbbbs if(zbbbs(i).le.zbmin) then zbmin=zbbbs(i) @@ -1471,8 +1525,8 @@ c end if end do else - zbmin=-1.0d+30 - zbmax=1.0d+30 + 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 @@ -1481,9 +1535,9 @@ c c c start with uncorrected normalized psi c - psinop=0.0d0 - psinxp=1.0d0 - psiant=1.0d0 + psinop=0.0_wp_ + psinxp=1.0_wp_ + psiant=1.0_wp_ c c search for O-point c @@ -1499,16 +1553,16 @@ c r10=rbmin z10=zbmin call points_ox(r10,z10,rxp,zxp,psinxptmp,info) - if(psinxp.ne.-1.0d0) then + 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.0d0 + psin1=1.0_wp_ r10=rmaxis - z10=(zbmax+zmaxis)/2.0d0 + z10=(zbmax+zmaxis)/2.0_wp_ call points_tgo(r10,z10,r1,z1,psin1,info) rbmax=r1 zbmax=z1 @@ -1520,14 +1574,14 @@ c print'(a)','no X-point' r10=rmop z10=zbmax call points_ox(r10,z10,rxp,zxp,psinxptmp,info) - if(psinxp.ne.-1.0d0) then + 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.0d0 - z10=(zbmin+zmaxis)/2.0d0 + psin1=1.0_wp_ + z10=(zbmin+zmaxis)/2.0_wp_ call points_tgo(r10,z10,r1,z1,psin1,info) zbmin=z1 else @@ -1538,15 +1592,15 @@ c print'(a)','no X-point' end if c if (ixp.eq.0) then - psin1=1.0d0 + psin1=1.0_wp_ psinop=psinoptmp psiant=psin1-psinop r10=rmaxis - z10=(zbmax+zmaxis)/2.0d0 + z10=(zbmax+zmaxis)/2.0_wp_ call points_tgo(r10,z10,r1,z1,psin1,info) zbmax=z1 rbmax=r1 - z10=(zbmin+zmaxis)/2.0d0 + z10=(zbmin+zmaxis)/2.0_wp_ call points_tgo(r10,z10,r1,z1,psin1,info) zbmin=z1 rbmin=r1 @@ -1573,24 +1627,23 @@ c compute flux surface averaged quantities rup=rmaxis rlw=rmaxis - zup=zmaxis+(zbmax-zmaxis)/10.0d0 - zlw=zmaxis-(zmaxis-zbmin)/10.0d0 + zup=zmaxis+(zbmax-zmaxis)/10.0_wp_ + zlw=zmaxis-(zmaxis-zbmin)/10.0_wp_ call flux_average c ipr=1 -c call contours_psi(1.0d0,rcon,zcon,ipr) +c call contours_psi(1.0_wp_,rcon,zcon,ipr) c do ii=1,2*np+1 c write(52,*) rcon(ii), zcon(ii) c end do c - c locate psi surface for q=1.5 and q=2 rup=rmaxis rlw=rmaxis - zup=(zbmax+zmaxis)/2.0d0 - zlw=(zmaxis+zbmin)/2.0d0 - q2=2.0d0 - q15=1.5d0 + zup=(zbmax+zmaxis)/2.0_wp_ + zlw=(zmaxis+zbmin)/2.0_wp_ + q2=2.0_wp_ + q15=1.5_wp_ call vmaxmini(qpsi,nr,qmin,qmax,iqmn,iqmx) if (q15.gt.qmin.and.q15.lt.qmax) then call surfq(q15,psi15) @@ -1617,13 +1670,26 @@ c c c subroutine points_ox(rz,zz,rf,zf,psinvf,info) - use const_and_precisions, only : comp_eps + use const_and_precisions, only : wp_,comp_eps use minpack, only : hybrj1 - implicit real*8 (a-h,o-z) - parameter(n=2,ldfjac=n,lwa=(n*(n+13))/2) - dimension xvec(n),fvec(n),fjac(ldfjac,n),wa(lwa) + 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 +c common/external functions/variables + real(wp_) :: psinv +c external fcnox +c common/psival/psinv +c xvec(1)=rz xvec(2)=zz tol = sqrt(comp_eps) @@ -1637,18 +1703,25 @@ c call equinum_psi(rf,zf) 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 - - implicit real*8 (a-h,o-z) - dimension x(n),fvec(n),fjac(ldfjac,n) + implicit none +c arguments + integer :: n,iflag,ldfjac + real(wp_), dimension(n) :: x,fvec + real(wp_), dimension(ldfjac,n) :: fjac +c common/external functions/variables + real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz +c common/derip1/dpsidr,dpsidz common/derip2/ddpsidrr,ddpsidzz,ddpsidrz - +c select case(iflag) case(1) call equinum_derpsi(x(1),x(2),iflag) @@ -1670,13 +1743,26 @@ c c c subroutine points_tgo(rz,zz,rf,zf,psin,info) - use const_and_precisions, only : comp_eps + use const_and_precisions, only : wp_,comp_eps use minpack, only : hybrj1 - implicit real*8 (a-h,o-z) - parameter(n=2,ldfjac=n,lwa=(n*(n+13))/2) - dimension xvec(n),fvec(n),fjac(ldfjac,n),wa(lwa) + 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 @@ -1692,15 +1778,24 @@ c c c subroutine fcntgo(n,x,fvec,fjac,ldfjac,iflag) + use const_and_precisions, only : wp_ use interp_eqprof, only: psia - - implicit real*8 (a-h,o-z) - dimension x(n),fvec(n),fjac(ldfjac,n) + implicit none +c arguments + integer :: n,ldfjac,iflag + real(wp_), dimension(n) :: x,fvec + real(wp_), dimension(ldfjac,n) :: fjac +c internal variables + integer :: ii +c common/external functions/variables + real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz, + . psinv,h +c common/derip1/dpsidr,dpsidz common/derip2/ddpsidrr,ddpsidzz,ddpsidrz common/psival/psinv common/cnpsi/h - +c select case(iflag) case(1) call equinum_psi(x(1),x(2)) @@ -1717,31 +1812,38 @@ c 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,qpsi,nr use utils, only : intlin - - implicit real*8 (a-h,o-z) - parameter(eps=1.d-4) + implicit none +c local constants + real(wp_), parameter :: eps=1.e-4_wp_ +c local variables + integer :: i,ips,nst + real(wp_) :: psin,rhop,rhot,ajphi,te,qq +c common/external functions/variables + real(wp_) :: dens,ddens + real(wp_) :: frhotor,temp c common/dens/dens,ddens c write(55,*) ' #psi rhot ne Te q Jphi' - psin=0.0d0 - rhop=0.0d0 - rhot=0.0d0 + psin=0.0_wp_ + rhop=0.0_wp_ + rhot=0.0_wp_ call density(psin) call tor_curr_psi(eps,ajphi) te=temp(psin) qq=qpsi(1) c - write(55,111) psin,rhot,dens,te,qq,ajphi*1.d-6 + write(55,111) psin,rhot,dens,te,qq,ajphi*1.0e-6_wp_ c nst=nr do i=2,nst @@ -1760,18 +1862,29 @@ c end if rhot=frhotor(psin) call tor_curr_psi(psin,ajphi) - write(55,111) psin,rhot,dens,te,qq,ajphi*1.d-6 + write(55,111) psin,rhot,dens,te,qq,ajphi*1.0e-6_wp_ end do c return 111 format(12(1x,e12.5)) end - +c +c +c subroutine print_prof_an - implicit real*8 (a-h,o-z) - parameter(nst=51) + use const_and_precisions, only : wp_ + implicit none +c local constants + integer, parameter :: nst=51 +c local variables + integer :: i + real(wp_) :: psin,rhop,rhot,te +c common/external functions/variables + real(wp_) :: dens,ddens + real(wp_) :: temp,frhotor +c common/dens/dens,ddens - +c write(55,*) ' #psi rhot ne Te' do i=1,nst psin=dble(i-1)/dble(nst-1) @@ -1781,7 +1894,7 @@ c te=temp(psin) write(55,111) psin,rhot,dens,te end do - +c return 111 format(12(1x,e12.5)) end @@ -1789,40 +1902,52 @@ c c c subroutine surfq(qval,psival) + use const_and_precisions, only : wp_ use magsurf_data, only : npoints use interp_eqprof, only : nr,psinr,qpsi use utils, only : locate, intlin - - implicit real*8 (a-h,o-z) - dimension rcn(npoints),zcn(npoints) + implicit none +c arguments + real(wp_) qval,psival +c local variables + integer ncnt,i1,ipr + real(wp_), dimension(npoints) :: rcn,zcn +c + ncnt=(npoints-1)/2 c c locate psi surface for q=qval c - ncnt=(npoints-1)/2 call locate(qpsi,nr,qval,i1) call intlin(qpsi(i1),psinr(i1),qpsi(i1+1),psinr(i1+1), . qval,psival) ipr=1 call contours_psi(psival,rcn,zcn,ipr) +c return end c c c subroutine bfield_res + use const_and_precisions, only : wp_ use interp_eqprof, only : rv,zv,nr,nz,btotal - - implicit real*8 (a-h,o-z) - parameter(nnw=501,nnh=501) - parameter(icmx=2002) - dimension rrcb(icmx),zzcb(icmx),ncpts(10) + implicit none +c local constants + integer, parameter :: icmx=2002 +c local variables + integer j,k,n,nconts,inc,nctot + integer, dimension(10) :: ncpts + real(wp_) btmx,btmn,zzk,rrj,bbphi,bbr,bbz,bbb + real(wp_), dimension(icmx) :: rrcb,zzcb +c common/external functions/variables + real(wp_) :: bres c common/parbres/bres c c Btotal on psi grid c - btmx=-1.0d30 - btmn=1.0d30 + btmx=-1.0e30_wp_ + btmn=1.0e30_wp_ do j=1,nr rrj=rv(j) do k=1,nz @@ -1856,39 +1981,42 @@ c c c subroutine profdata + use const_and_precisions, only : wp_ use graydata_flags, only : iprof,iscal,nprof use graydata_par, only : psdbnd,factb,factt,factn use interp_eqprof, only : nsfd,npp,psnpp,aa,bb,cc, . psrad,derad,terad,zfc,tfn,cfn,ct,cz,alloc_profvec use simplespline, only : difcsn use dierckx, only : curfit,splev,splder - - implicit real*8 (a-h,o-z) + implicit none +c local variables + integer ierr,i,iopt,ier,kspl,npest,lwrkf,nu,nn,nn1,nn2 integer, dimension(:), allocatable :: iwrkf - real*8, dimension(:), allocatable :: wrkf,wf,wrkfd, - . densi,ddensi,d2densi + real(wp_) aat,aan,ffact,psrad0,terad0,derad0,zfc0,psradi, + . teradi,deradi,zfci,xb,xe,sspl,dnpp,ddnpp,d2dnpp,dpsb,fp + real(wp_), dimension(:), allocatable :: wf,wrkf,wrkfd,densi, + . ddensi,d2densi c c read plasma profiles from file if iprof>0 c if (iscal.eq.0) then - aat=2.0d0/3.0d0 - aan=4.0d0/3.0d0 + aat=2.0_wp_/3.0_wp_ + aan=4.0_wp_/3.0_wp_ else - aat=1.0d0 - aan=1.0d0 + aat=1.0_wp_ + aan=1.0_wp_ end if ffact=factb - if(iscal.eq.2) ffact=1.0d0 + if(iscal.eq.2) ffact=1.0_wp_ c if (iprof.gt.0) then read(nprof,*) npp - +c call alloc_profvec(ierr) if (ierr.ne.0) stop - +c npest=npp+4 lwrkf=npp*4+npest*16 - if(allocated(wrkf)) deallocate(wrkf) if(allocated(iwrkf)) deallocate(iwrkf) if(allocated(wf)) deallocate(wf) @@ -1896,24 +2024,23 @@ c if(allocated(densi)) deallocate(densi) if(allocated(ddensi)) deallocate(ddensi) if(allocated(d2densi)) deallocate(d2densi) - allocate(wrkf(lwrkf),iwrkf(npest),wf(npp),densi(npest), . wrkfd(npest),ddensi(npest),d2densi(npest)) - +c read(nprof,*) psrad0,terad0,derad0,zfc0 - if(psrad0.ne.0.0d0) psrad0=0.0d0 + if(psrad0.ne.0.0_wp_) psrad0=0.0_wp_ psrad(1)=psrad0 terad(1)=terad0*ffact**aat*factt derad(1)=derad0*ffact**aan*factn zfc(1)=zfc0 - wf(1)=1.0d0 + wf(1)=1.0_wp_ do i=2,npp read(nprof,*) psradi,teradi,deradi,zfci psrad(i)=psradi terad(i)=teradi*ffact**aat*factt derad(i)=deradi*ffact**aan*factn zfc(i)=zfci - wf(i)=1.0d0 + wf(i)=1.0_wp_ end do c c spline approximation of temperature and Zeff @@ -1927,10 +2054,10 @@ c c spline approximation of density c iopt=0 - xb=0.0d0 + xb=0.0_wp_ xe=psrad(npp) kspl=3 - sspl=.001d0 + sspl=.001_wp_ c call curfit(iopt,npp,psrad,derad,wf,xb,xe,kspl,sspl,npest,nsfd, . tfn,cfn,fp,wrkf,lwrkf,iwrkf,ier) @@ -1945,7 +2072,7 @@ c call splder(tfn,nsfd,cfn,3,nu,psrad,d2densi,npp,wrkfd,ier) d2dnpp=d2densi(npp) - if(derad(npp).eq.0.0d0) then + if(derad(npp).eq.0.0_wp_) then psdbnd=psrad(npp) else psnpp=psrad(npp) @@ -1954,11 +2081,11 @@ c nn1=nn+1 nn2=nn+2 aa=(nn1*nn2*dnpp+2*nn1*ddnpp*dpsb+d2dnpp*dpsb**2) - aa=aa/(-dpsb)**nn/2.0d0 + aa=aa/(-dpsb)**nn/2.0_wp_ bb=-(nn*nn2*dnpp+(2*nn+1)*ddnpp*dpsb+d2dnpp*dpsb**2) bb=bb/(-dpsb)**nn1 cc=(nn1*nn*dnpp+2*nn*ddnpp*dpsb+d2dnpp*dpsb**2) - cc=cc/(-dpsb)**nn2/2.0d0 + cc=cc/(-dpsb)**nn2/2.0_wp_ end if c end if @@ -1971,11 +2098,19 @@ c c c subroutine rhotor(nnr) + use const_and_precisions, only : wp_ use interp_eqprof, only : nr,psinr,qpsi,crhot,cq use simplespline, only : difcsn - - implicit real*8(a-h,o-z) - real*8, dimension(:), allocatable :: rhotnr + implicit none +c arguments + integer nnr +c local variables + integer iopt,ier,k + real(wp_) dx,drhot + real(wp_), dimension(:), allocatable :: rhotnr +c common/external functions/variables + real(wp_) rhotsx +c common/rhotsx/rhotsx c if(allocated(cq)) deallocate(cq) @@ -1987,11 +2122,11 @@ c iopt=0 call difcsn(psinr,qpsi,nr,nnr,iopt,cq,ier) c - rhotnr(1)=0.0d0 + rhotnr(1)=0.0_wp_ do k=1,nnr-1 dx=psinr(k+1)-psinr(k) - drhot=dx*(cq(k,1)+dx*(cq(k,2)/2.0d0+dx*(cq(k,3)/3.0d0 - . +dx*cq(k,4)/4.0d0))) + 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(nnr) @@ -2003,31 +2138,46 @@ c spline interpolation of rhotor c iopt=0 call difcsn(psinr,rhotnr,nr,nnr,iopt,crhot,ier) - +c deallocate(rhotnr) - +c return end - +c +c +c function fq_eq(psi) + use const_and_precisions, only : wp_ use interp_eqprof, only : psinr,nr,cq use simplespline, only :spli - - implicit real*8(a-h,o-z) - + implicit none +c arguments + real(wp_) psi,fq_eq +c local variables + integer irt + real(wp_) dps +c irt=int((nr-1)*psi+1) if(irt.eq.0) irt=1 if(irt.eq.nr) irt=nr-1 dps=psi-psinr(irt) fq_eq=spli(cq,nr,irt,dps) +c return end - +c +c +c function frhotor_eq(psi) + use const_and_precisions, only : wp_ use interp_eqprof, only : psinr,nr,crhot use simplespline, only :spli - - implicit real*8(a-h,o-z) + implicit none +c arguments + real(wp_) :: psi,frhotor_eq +c local variables + integer :: irt + real(wp_) :: dps c irt=int((nr-1)*psi+1) c if(irt.eq.0) irt=1 @@ -2035,24 +2185,40 @@ c if(irt.eq.nr) irt=nr-1 irt=min(max(1,irt),nr-1) dps=psi-psinr(irt) frhotor_eq=spli(crhot,nr,irt,dps) +c return end - +c +c +c function frhotor(psi) + use const_and_precisions, only : wp_ use graydata_flags, only : iequil - - implicit real*8(a-h,o-z) + implicit none +c arguments + real(wp_) :: psi,frhotor +c common/external functions/variables + real(wp_) :: frhotor_eq,frhotor_an +c if(iequil.eq.2) frhotor=frhotor_eq(psi) if(iequil.eq.1) frhotor=frhotor_an(sqrt(psi)) +c return end - +c +c +c function frhotor_av(psi) + use const_and_precisions, only : wp_ use magsurf_data, only : rpstab, crhotq, npsi use simplespline, only :spli - - implicit real*8(a-h,o-z) - + implicit none +c arguments + real(wp_) :: psi,frhotor_av +c local variables + integer :: ip + real(wp_) :: rpsi,dps +c rpsi=sqrt(psi) ip=int((npsi-1)*rpsi+1) c if(ip.eq.0) ip=1 @@ -2060,75 +2226,114 @@ c if(ip.eq.npsi) ip=npsi-1 ip=min(max(1,ip),npsi-1) dps=rpsi-rpstab(ip) frhotor_av=spli(crhotq,npsi,ip,dps) +c return end - +c +c +c subroutine rhopol use dierckx, only : curfit,splev - implicit real*8(a-h,o-z) - parameter(nnr=101,nrest=nnr+4) - parameter(lwrkp=nnr*4+nrest*16) - dimension rhop(nnr),rhot(nnr),rhopi(nnr) - dimension trp(nrest),crp(nrest),wp(nrest) - dimension wrkp(lwrkp),iwrkp(nrest) + use const_and_precisions, only : wp_ + 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_), dimension(nrest) :: trp,crp + real(wp_) :: frhotor +c common/coffrtp/trp common/coffrn/nsrp common/coffrp/crp - - dr=1.0d0/dble(nnr-1) +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(psin) - wp(i)=1.0d0 + wp(i)=1.0_wp_ end do - wp(1)=1.0d3 - wp(nnr)=1.0d3 + wp(1)=1.0e3_wp_ + wp(nnr)=1.0e3_wp_ c spline interpolation of rhopol versus rhotor iopt=0 - xb=0.0d0 - xe=1.0d0 - ss=0.00001d0 + 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) + 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 - implicit real*8(a-h,o-z) - parameter(nnr=101,nrest=nnr+4) - dimension trp(nrest),crp(nrest),rrs(1),ffspl(1) + 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 - rrs(1)=rhot - call splev(trp,nsrp,crp,3,rrs,ffspl,1,ier) - frhopol=ffspl(1) - return - end - - subroutine cniteq(h, ncon, npts, icount, rcon, zcon,ichoi) - use interp_eqprof, only : dr,dz,nr,nz,rqgrid=>rv,zqgrid=>zv,psin, - . btotal - - implicit real*8 (a-h,o-z) 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(h,ncon,npts,icount,rcon,zcon,ichoi) c v2.01 12/07/95 -- written by d v bartlett, jet joint undertaking. c (based on an older code) c - parameter(icmx=2002) - dimension ja(3,2),lx(1000),npts(10) - dimension rcon(icmx),zcon(icmx) - real*8, dimension(:), allocatable :: a + use const_and_precisions, only : wp_ + use interp_eqprof, only : dr,dz,nr,nz,rqgrid=>rv,zqgrid=>zv,psin, + . btotal + implicit none +c local constants + integer, parameter :: icmx=2002 +c arguments + integer :: ncon,icount,ichoi + integer, dimension(10) :: npts + real(wp_) :: h + real(wp_), dimension(icmx) :: rcon,zcon +c local variables + integer i,j,k,l,ico,nrqmax,nreq,nzeq,iclast,mpl,ix,jx, + . mxr,n1,jm,jfor,lda,ldb,jabs,jnb,kx,ikx,itm,inext,in + integer, dimension(3,2) :: ja + integer, dimension(1000) :: lx + real(wp_) :: drgrd,dzgrd,ah,adn,px,x,y + real(wp_), dimension(:), allocatable :: a c data px/0.5d0/ c @@ -2152,8 +2357,8 @@ c endif c do ico=1,icmx - rcon(ico)=0.0d0 - zcon(ico)=0.0d0 + rcon(ico)=0.0_wp_ + zcon(ico)=0.0_wp_ enddo c nrqmax=nr @@ -2311,17 +2516,24 @@ c c c subroutine contours_psi(h,rcn,zcn,ipr) + use const_and_precisions, only : wp_,pi use graydata_par, only : rwallm use magsurf_data, only : npoints use interp_eqprof, only : psia,psiant,psinop,nsr=>nsrt,nsz=>nszt, . cc=>cceq,tr,tz,rup,zup,rlw,zlw,nr - use const_and_precisions, only : pi use dierckx, only : profil,sproota - implicit real*8 (a-h,o-z) - parameter(mest=4,kspl=3) - dimension zeroc(mest) - dimension rcn(npoints),zcn(npoints) - real*8, dimension(:), allocatable :: czc + implicit none +c local constants + integer, parameter :: mest=4,kspl=3 +c arguments + integer :: ipr + real(wp_) :: h + real(wp_), dimension(npoints) :: rcn,zcn +c local variables + integer :: np,info,ic,ier,ii,iopt,m + real(wp_) :: ra,rb,za,zb,th,zc,val + real(wp_), dimension(mest) :: zeroc + real(wp_), dimension(:), allocatable :: czc c np=(npoints-1)/2 if(allocated(czc)) deallocate(czc) @@ -2333,7 +2545,7 @@ c zb=zlw call points_tgo(ra,za,rup,zup,h,info) call points_tgo(rb,zb,rlw,zlw,h,info) - +c th=pi/dble(np) rcn(1)=rlw zcn(1)=zlw @@ -2342,7 +2554,7 @@ c rcn(np+1)=rup zcn(np+1)=zup do ic=2,np - zc=zlw+(zup-zlw)*(1.0d0-cos(th*(ic-1)))/2.0d0 + zc=zlw+(zup-zlw)*(1.0_wp_-cos(th*(ic-1)))/2.0_wp_ iopt=1 call profil(iopt,tr,nsr,tz,nsz,cc,kspl,kspl,zc,nsr,czc,ier) if(ier.gt.0) print*,' profil =',ier @@ -2367,9 +2579,9 @@ c write(71,*) write(71,*) end if - +c deallocate(czc) - +c return 111 format(i6,12(1x,e12.5)) end @@ -2377,31 +2589,36 @@ c c c subroutine flux_average - use const_and_precisions, only : zero,one,pi,ccj=>mu0inv + use const_and_precisions, only : wp_,zero,one,pi,ccj=>mu0inv use magsurf_data use interp_eqprof, only : btrcen,btaxis,rmaxis,zmaxis,phitedge use simplespline, only : difcs use dierckx, only : regrid,coeff_parder - - implicit real*8 (a-h,o-z) - real*8 lam - real*8, dimension(:), allocatable :: rctemp,zctemp - - parameter(nnintp=101,ncnt=100,nlam=41) - parameter(ksp=3,njest=nnintp+ksp+1,nlest=nlam+ksp+1) - parameter(lwrk=4*(nnintp+nlam)+11*(njest+nlest)+ - . njest*nnintp+nlest+54) - parameter(kwrk=nnintp+nlam+njest+nlest+3) - parameter(lw01=nnintp*4+nlam*3+nnintp*nlam) -c - dimension dadrhotv(nnintp) - dimension dvdrhotv(nnintp) - dimension dlpv(2*ncnt),bv(2*ncnt+1),bpv(2*ncnt+1) - dimension vratjpl(nnintp) - dimension alam(nlam),fhlam(nnintp,nlam) - dimension ffhlam(nnintp*nlam),dffhlam(nnintp*nlam) - dimension iwrk(kwrk),wrk(lwrk) - dimension weights(nlam) + implicit none +c local constants + integer, parameter :: nnintp=101,ncnt=100,nlam=41,ksp=3, + . njest=nnintp+ksp+1,nlest=nlam+ksp+1, + . lwrk=4*(nnintp+nlam)+11*(njest+nlest)+njest*nnintp+nlest+54, + . kwrk=nnintp+nlam+njest+nlest+3,lw01=nnintp*4+nlam*3+nnintp*nlam +c local variables + integer ier,ierr,l,jp,ipr,jpr,inc,inc1,iopt,njp,nlm,ninpr + integer, dimension(kwrk) :: iwrk + real(wp_) :: lam,dlam,anorm,b2av,dvdpsi,dadpsi,ratio_cdator, + . 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, + . bphi,brr,bzz,riav,fp + real(wp_), dimension(nnintp) :: dadrhotv,dvdrhotv,vratjpl + real(wp_), dimension(2*ncnt) :: dlpv + real(wp_), dimension(2*ncnt+1) :: bv,bpv + real(wp_), dimension(nlam) :: alam,weights + real(wp_), dimension(nnintp,nlam) :: fhlam + real(wp_), dimension(nnintp*nlam) :: ffhlam,dffhlam + real(wp_), dimension(lwrk) :: wrk + real(wp_), dimension(:), allocatable :: rctemp,zctemp +c common/external functions/variables + real(wp_) :: fpolv,ffpv c common/fpol/fpolv,ffpv c @@ -2422,67 +2639,67 @@ c computation of flux surface averaged quantities write(71,*)' #i psin R z' - dlam=1.0d0/dble(nlam-1) + dlam=1.0_wp_/dble(nlam-1) do l=1,nlam-1 alam(l)=dble(l-1)*dlam - fhlam(1,l)=sqrt(1.0d0-alam(l)) + fhlam(1,l)=sqrt(1.0_wp_-alam(l)) ffhlam(l)=fhlam(1,l) - dffhlam(l)=-0.5d0/sqrt(1.0d0-alam(l)) - weights(l)=1.0d0 + dffhlam(l)=-0.5_wp_/sqrt(1.0_wp_-alam(l)) + weights(l)=1.0_wp_ end do - weights(1)=0.5d0 - weights(nlam)=0.5d0 - alam(nlam)=1.0d0 - fhlam(1,nlam)=0.0d0 - ffhlam(nlam)=0.0d0 - dffhlam(nlam)=-99999.0d0 + weights(1)=0.5_wp_ + weights(nlam)=0.5_wp_ + alam(nlam)=1.0_wp_ + fhlam(1,nlam)=0.0_wp_ + ffhlam(nlam)=0.0_wp_ + dffhlam(nlam)=-99999.0_wp_ jp=1 - anorm=2.0d0*pi*rmaxis/abs(btaxis) + anorm=2.0_wp_*pi*rmaxis/abs(btaxis) b2av=btaxis**2 - dvdpsi=2.0d0*pi*anorm - dadpsi=2.0d0*pi/abs(btaxis) + dvdpsi=2.0_wp_*pi*anorm + dadpsi=2.0_wp_*pi/abs(btaxis) ratio_cdator=abs(btaxis/btrcen) - ratio_cdbtor=1.0d0 - ratio_pltor=1.0d0 - qq=1.0d0 - fc=1.0d0 + ratio_cdbtor=1.0_wp_ + ratio_pltor=1.0_wp_ + qq=1.0_wp_ + fc=1.0_wp_ - psicon(1)=0.0d0 - rcon(1,:)=0.0d0 - zcon(1,:)=0.0d0 - pstab(1)=0.0d0 - rhot_eq(1)=0.0d0 - rpstab(1)=0.0d0 - rhotqv(1)=0.0d0 - vcurrp(1)=0.0d0 - vajphiav(1)=0.0d0 + psicon(1)=0.0_wp_ + rcon(1,:)=0.0_wp_ + zcon(1,:)=0.0_wp_ + pstab(1)=0.0_wp_ + rhot_eq(1)=0.0_wp_ + rpstab(1)=0.0_wp_ + rhotqv(1)=0.0_wp_ + vcurrp(1)=0.0_wp_ + vajphiav(1)=0.0_wp_ bmxpsi(1)=abs(btaxis) bmnpsi(1)=abs(btaxis) bav(1)=abs(btaxis) - rbav(1)=1.0d0 + rbav(1)=1.0_wp_ rri(1)=rmaxis - varea(1)=0.0d0 - vvol(1)=0.0d0 + varea(1)=0.0_wp_ + vvol(1)=0.0_wp_ vratjpl(1)=ratio_pltor vratja(1)=ratio_cdator vratjb(1)=ratio_cdbtor ffc(1)=fc - rhot2q=0.0d0 - qqv(1)=1.0d0 + rhot2q=0.0_wp_ + qqv(1)=1.0_wp_ - dadrhotv(1)=0.0d0 - dvdrhotv(1)=0.0d0 + dadrhotv(1)=0.0_wp_ + dvdrhotv(1)=0.0_wp_ -C rup=rmaxis -C rlw=rmaxis -C zup=(zbmax+zmaxis)/2.0d0 -C zlw=(zmaxis+zbmin)/2.0d0 +c rup=rmaxis +c rlw=rmaxis +c zup=(zbmax+zmaxis)/2.0_wp_ +c zlw=(zmaxis+zbmin)/2.0_wp_ do jp=2,npsi height=dble(jp-1)/dble(npsi-1) psicon(jp)=height - if(jp.eq.npsi) height=0.9999d0 + if(jp.eq.npsi) height=0.9999_wp_ ipr=0 jpr=mod(jp,ninpr) if(jpr.eq.1) ipr=1 @@ -2491,17 +2708,17 @@ C zlw=(zmaxis+zbmin)/2.0d0 rcon(jp,:) = rctemp zcon(jp,:) = zctemp c - r2iav=0.0d0 - anorm=0.0d0 - dadpsi=0.0d0 - currp=0.0d0 - b2av=0.0d0 - area=0.0d0 - volume=0.0d0 - ajphiav=0.0d0 - bbav=0.0d0 - bmmx=-1.0d+30 - bmmn=1.0d+30 + r2iav=0.0_wp_ + anorm=0.0_wp_ + dadpsi=0.0_wp_ + currp=0.0_wp_ + b2av=0.0_wp_ + area=0.0_wp_ + volume=0.0_wp_ + ajphiav=0.0_wp_ + bbav=0.0_wp_ + bmmx=-1.0e+30_wp_ + bmmn=1.0e+30_wp_ call tor_curr(rctemp(1),zctemp(1),ajphi0) call bpol(rctemp(1),zctemp(1),brr,bzz) @@ -2523,7 +2740,7 @@ c c compute length, area and volume defined by psi=height^2 - ph=0.5d0*(dla+dlb+dlp) + ph=0.5_wp_*(dla+dlb+dlp) area2=ph*(ph-dla)*(ph-dlb)*(ph-dlp) area=area+sqrt(area2) rzp=rctemp(inc1)*zctemp(inc1) @@ -2544,15 +2761,15 @@ c compute line integral on the contour psi=height^2 bv(inc1)=btot bpv(inc1)=bpoloid - dlph=0.5d0*dlp - anorm=anorm+dlph*(1.0d0/bpoloid+1.0d0/bpoloid0) + dlph=0.5_wp_*dlp + anorm=anorm+dlph*(1.0_wp_/bpoloid+1.0_wp_/bpoloid0) dadpsi=dadpsi+dlph* - . (1.0d0/(bpoloid*rpsim)+1.0d0/(bpoloid0*rpsim0)) + . (1.0_wp_/(bpoloid*rpsim)+1.0_wp_/(bpoloid0*rpsim0)) currp=currp+dlph*(bpoloid+bpoloid0) b2av=b2av+dlph*(btot0**2/bpoloid0+btot**2/bpoloid) bbav=bbav+dlph*(btot/bpoloid+btot0/bpoloid0) - r2iav=r2iav+dlph* - . (1.0d0/(bpoloid*rpsim**2)+1.0d0/(bpoloid0*rpsim0**2)) + r2iav=r2iav+dlph*(1.0_wp_/(bpoloid*rpsim**2) + . +1.0_wp_/(bpoloid0*rpsim0**2)) ajphiav=ajphiav+dlph* . (ajphi0/(bpoloid0*rpsim0)+ajphi/(bpoloid*rpsim)) @@ -2576,7 +2793,7 @@ c currp = plasma current within psi=const bbav=bbav/anorm r2iav=r2iav/anorm - dvdpsi=2.0d0*pi*anorm + dvdpsi=2.0_wp_*pi*anorm riav=dadpsi/anorm b2av=b2av/anorm vcurrp(jp)=ccj*currp @@ -2621,37 +2838,37 @@ c computation of fraction of circulating/trapped fraction fc, ft c and of function H(lambda,rhop) c ffhlam = Bmn/Bmx/fc integral_lambda^1 dlam/ - fc=0.0d0 - shlam=0.0d0 + fc=0.0_wp_ + shlam=0.0_wp_ do l=nlam,1,-1 lam=alam(l) - srl=0.0d0 - rl2=1.0d0-lam*bv(1)/bmmx - rl0=0.d0 + srl=0.0_wp_ + rl2=1.0_wp_-lam*bv(1)/bmmx + rl0=0.0_wp_ if(rl2.gt.0) rl0=sqrt(rl2) do inc=1,npoints-1 - rl2=1.0d0-lam*bv(inc+1)/bmmx - rl=0.0d0 + rl2=1.0_wp_-lam*bv(inc+1)/bmmx + rl=0.0_wp_ if(rl2.gt.0) rl=sqrt(rl2) - srl=srl+0.5d0*dlpv(inc)*(rl/bpv(inc+1)+rl0/bpv(inc)) + srl=srl+0.5_wp_*dlpv(inc)*(rl/bpv(inc+1)+rl0/bpv(inc)) rl0=rl end do srl=srl/anorm - dhlam=0.5d0/srl + dhlam=0.5_wp_/srl fc=fc+lam/srl*weights(l) if(l.eq.nlam) then - fhlam(jp,l)=0.0d0 - ffhlam(nlam*(jp-1)+l)=0.0d0 + fhlam(jp,l)=0.0_wp_ + ffhlam(nlam*(jp-1)+l)=0.0_wp_ dffhlam(nlam*(jp-1)+l)=-dhlam dhlam0=dhlam else - shlam=shlam+0.5d0*(dhlam+dhlam0)*dlam + shlam=shlam+0.5_wp_*(dhlam+dhlam0)*dlam fhlam(jp,l)=shlam dffhlam(nlam*(jp-1)+l)=-dhlam dhlam0=dhlam end if end do - fc=0.75d0*b2av/bmmx**2*fc*dlam + fc=0.75_wp_*b2av/bmmx**2*fc*dlam ffc(jp)=fc ccfh=bmmn/bmmx/fc @@ -2669,9 +2886,9 @@ c ffhlam = Bmn/Bmx/fc integral_lambda^1 dlam/ do jp=1,npsi rhotqv(jp)=rhotqv(jp)/rhotqv(npsi) if(jp.eq.npsi) then - rhotqv(jp)=1.0d0 - rpstab(jp)=1.0d0 - pstab(jp)=1.0d0 + rhotqv(jp)=1.0_wp_ + rpstab(jp)=1.0_wp_ + pstab(jp)=1.0_wp_ end if rhot_eq(jp)=frhotor_eq(pstab(jp)) write(56,99) pstab(jp),rhot_eq(jp),rhotqv(jp), @@ -2716,7 +2933,7 @@ c spline coefficients of rhot c spline interpolation of H(lambda,rhop) and dH/dlambda iopt=0 - s=0.0d0 + s=0.0_wp_ call regrid(iopt,npsi,rpstab,nlam,alam,ffhlam, . zero,one,zero,one,ksp,ksp,s, . njest,nlest,njp,tjp,nlm,tlm,ch,fp,wrk,lwrk,iwrk,kwrk,ier) @@ -2725,18 +2942,24 @@ c spline interpolation of H(lambda,rhop) and dH/dlambda call coeff_parder(tjp,njp,tlm,nlm,ch,ksp,ksp,0,1, . ch01,lw01,ier) - - +c return 99 format(20(1x,e12.5)) end - +c +c +c function fdadrhot(rpsi) + use const_and_precisions, only : wp_ use magsurf_data, only : rpstab,cdadrhot,npsi use simplespline, only :spli - - implicit real*8(a-h,o-z) - + implicit none +c arguments + real(wp_) :: rpsi,fdadrhot +c local variables + integer :: ip + real(wp_) :: dps +c ip=int((npsi-1)*rpsi+1) c if(ip.eq.0) ip=1 c if(ip.eq.npsi) ip=npsi-1 @@ -2745,52 +2968,66 @@ c if(ip.eq.npsi) ip=npsi-1 fdadrhot=spli(cdadrhot,npsi,ip,dps) return end - +c +c +c function fdvdrhot(rpsi) + use const_and_precisions, only : wp_ use magsurf_data, only : rpstab,cdvdrhot,npsi use simplespline, only :spli - - implicit real*8(a-h,o-z) - + implicit none +c arguments + real(wp_) :: rpsi,fdvdrhot +c local variables + integer :: ip + real(wp_) :: dps +c ip=int((npsi-1)*rpsi+1) ip=min(max(1,ip),npsi-1) dps=rpsi-rpstab(ip) fdvdrhot=spli(cdvdrhot,npsi,ip,dps) +c return end - +c +c +c subroutine flux_average_an + 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 simplespline, only : difcs use dierckx, only : regrid,coeff_parder - - implicit real*8 (a-h,o-z) - real*8 lam - real*8, dimension(:), allocatable :: rctemp,zctemp - - parameter(nnintp=101,ncnt=100,nlam=41) - parameter(zero=0.0d0,one=1.0d0) - parameter(pi=3.14159265358979d0,ccj=1.0d+7/(4.0d0*pi)) - parameter(ksp=3,njest=nnintp+ksp+1,nlest=nlam+ksp+1) - parameter(lwrk=4*(nnintp+nlam)+11*(njest+nlest)+ - . njest*nnintp+nlest+54) - parameter(kwrk=nnintp+nlam+njest+nlest+3) - parameter(lw01=nnintp*4+nlam*3+nnintp*nlam) - - dimension dadrhotv(nnintp) - dimension dvdrhotv(nnintp) - dimension dlpv(2*ncnt),bv(2*ncnt+1),bpv(2*ncnt+1) - dimension vratjpl(nnintp) - dimension alam(nlam),fhlam(nnintp,nlam) - dimension ffhlam(nnintp*nlam),dffhlam(nnintp*nlam) - dimension iwrk(kwrk),wrk(lwrk) - dimension weights(nlam) + implicit none +c local constants + integer, parameter :: nnintp=101,ncnt=100,nlam=41,ksp=3, + . njest=nnintp+ksp+1,nlest=nlam+ksp+1, + . lwrk=4*(nnintp+nlam)+11*(njest+nlest)+njest*nnintp+nlest+54, + . kwrk=nnintp+nlam+njest+nlest+3,lw01=nnintp*4+nlam*3+nnintp*nlam +c local variables + integer ier,ierr,l,jp,ipr,jpr,inc,inc1,iopt,njp,nlm,ninpr + integer, dimension(kwrk) :: iwrk + real(wp_) :: lam,dlam,anorm,b2av,dvdpsi,dadpsi,ratio_cdator, + . ratio_cdbtor,ratio_pltor,qq,fc,rhot2q,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,ajphi,bphi,brr,bzz,riav,fp, + . btaxis,dadr,dbvcdc13,dbvcdc31,dvdr,height,phitedge,qqan,rmaxis, + . zmaxis,rn + real(wp_), dimension(nnintp) :: dadrhotv,dvdrhotv,vratjpl + real(wp_), dimension(2*ncnt) :: dlpv + real(wp_), dimension(2*ncnt+1) :: bv,bpv + real(wp_), dimension(nlam) :: alam,weights + real(wp_), dimension(nnintp,nlam) :: fhlam + real(wp_), dimension(nnintp*nlam) :: ffhlam,dffhlam + real(wp_), dimension(lwrk) :: wrk + real(wp_), dimension(:), allocatable :: rctemp,zctemp +c common/external functions/variables + real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv + real(wp_) :: frhotor_an c common/fpol/fpolv,ffpv -c common/psival/psin - common/derip1/dpsidr,dpsidz common/derip2/ddpsidrr,ddpsidzz,ddpsidrz c @@ -2818,66 +3055,66 @@ c computation of flux surface averaged quantities write(71,*)' #i psin R z' - dlam=1.0d0/dble(nlam-1) + dlam=1.0_wp_/dble(nlam-1) do l=1,nlam-1 alam(l)=dble(l-1)*dlam - fhlam(1,l)=sqrt(1.0d0-alam(l)) + fhlam(1,l)=sqrt(1.0_wp_-alam(l)) ffhlam(l)=fhlam(1,l) - dffhlam(l)=-0.5d0/sqrt(1.0d0-alam(l)) - weights(l)=1.0d0 + dffhlam(l)=-0.5_wp_/sqrt(1.0_wp_-alam(l)) + weights(l)=1.0_wp_ end do - weights(1)=0.5d0 - weights(nlam)=0.5d0 - alam(nlam)=1.0d0 - fhlam(1,nlam)=0.0d0 - ffhlam(nlam)=0.0d0 - dffhlam(nlam)=-99999.0d0 + weights(1)=0.5_wp_ + weights(nlam)=0.5_wp_ + alam(nlam)=1.0_wp_ + fhlam(1,nlam)=0.0_wp_ + ffhlam(nlam)=0.0_wp_ + dffhlam(nlam)=-99999.0_wp_ jp=1 - anorm=2.0d0*pi*rmaxis/abs(btaxis) + anorm=2.0_wp_*pi*rmaxis/abs(btaxis) b2av=btaxis**2 - dvdpsi=2.0d0*pi*anorm - dadpsi=2.0d0*pi/abs(btaxis) + dvdpsi=2.0_wp_*pi*anorm + dadpsi=2.0_wp_*pi/abs(btaxis) ratio_cdator=abs(btaxis/btrcen) - ratio_cdbtor=1.0d0 - ratio_pltor=1.0d0 - qq=1.0d0 - fc=1.0d0 + ratio_cdbtor=1.0_wp_ + ratio_pltor=1.0_wp_ + qq=1.0_wp_ + fc=1.0_wp_ - psicon(1)=0.0d0 - rcon(1,:)=0.0d0 - zcon(1,:)=0.0d0 - pstab(1)=0.0d0 - rhot_eq(1)=0.0d0 - rpstab(1)=0.0d0 - rhotqv(1)=0.0d0 - vcurrp(1)=0.0d0 - vajphiav(1)=0.0d0 + psicon(1)=0.0_wp_ + rcon(1,:)=0.0_wp_ + zcon(1,:)=0.0_wp_ + pstab(1)=0.0_wp_ + rhot_eq(1)=0.0_wp_ + rpstab(1)=0.0_wp_ + rhotqv(1)=0.0_wp_ + vcurrp(1)=0.0_wp_ + vajphiav(1)=0.0_wp_ bmxpsi(1)=abs(btaxis) bmnpsi(1)=abs(btaxis) bav(1)=abs(btaxis) - rbav(1)=1.0d0 + rbav(1)=1.0_wp_ rri(1)=rmaxis - varea(1)=0.0d0 - vvol(1)=0.0d0 + varea(1)=0.0_wp_ + vvol(1)=0.0_wp_ vratjpl(1)=ratio_pltor vratja(1)=ratio_cdator vratjb(1)=ratio_cdbtor ffc(1)=fc - rhot2q=0.0d0 - dadrhotv(1)=0.0d0 - dvdrhotv(1)=0.0d0 - qqv(1)=1.0d0 + rhot2q=0.0_wp_ + dadrhotv(1)=0.0_wp_ + dvdrhotv(1)=0.0_wp_ + qqv(1)=1.0_wp_ -C rup=rmaxis -C rlw=rmaxis -C zup=(zbmax+zmaxis)/2.0d0 -C zlw=(zmaxis+zbmin)/2.0d0 +c rup=rmaxis +c rlw=rmaxis +c zup=(zbmax+zmaxis)/2.0_wp_ +c zlw=(zmaxis+zbmin)/2.0_wp_ do jp=2,npsi height=dble(jp-1)/dble(npsi-1) psicon(jp)=height - if(jp.eq.npsi) height=0.9999d0 + if(jp.eq.npsi) height=0.9999_wp_ height2=height*height ipr=0 jpr=mod(jp,ninpr) @@ -2886,17 +3123,17 @@ C zlw=(zmaxis+zbmin)/2.0d0 rcon(jp,:) = rctemp zcon(jp,:) = zctemp c - r2iav=0.0d0 - anorm=0.0d0 - dadpsi=0.0d0 - currp=0.0d0 - b2av=0.0d0 - area=0.0d0 - volume=0.0d0 - ajphiav=0.0d0 - bbav=0.0d0 - bmmx=-1.0d+30 - bmmn=1.0d+30 + r2iav=0.0_wp_ + anorm=0.0_wp_ + dadpsi=0.0_wp_ + currp=0.0_wp_ + b2av=0.0_wp_ + area=0.0_wp_ + volume=0.0_wp_ + ajphiav=0.0_wp_ + bbav=0.0_wp_ + bmmx=-1.0e+30_wp_ + bmmn=1.0e+30_wp_ call equian(rctemp(1),zctemp(1)) dbvcdc13=-ddpsidzz/rctemp(1) @@ -2921,7 +3158,7 @@ c c compute length, area and volume defined by psi=height^2 - ph=0.5d0*(dla+dlb+dlp) + ph=0.5_wp_*(dla+dlb+dlp) area2=ph*(ph-dla)*(ph-dlb)*(ph-dlp) area=area+sqrt(area2) rzp=rctemp(inc1)*zctemp(inc1) @@ -2945,15 +3182,15 @@ c compute line integral on the contour psi=height^2 bv(inc1)=btot bpv(inc1)=bpoloid - dlph=0.5d0*dlp - anorm=anorm+dlph*(1.0d0/bpoloid+1.0d0/bpoloid0) + dlph=0.5_wp_*dlp + anorm=anorm+dlph*(1.0_wp_/bpoloid+1.0_wp_/bpoloid0) dadpsi=dadpsi+dlph* - . (1.0d0/(bpoloid*rpsim)+1.0d0/(bpoloid0*rpsim0)) + . (1.0_wp_/(bpoloid*rpsim)+1.0_wp_/(bpoloid0*rpsim0)) currp=currp+dlph*(bpoloid+bpoloid0) b2av=b2av+dlph*(btot0**2/bpoloid0+btot**2/bpoloid) bbav=bbav+dlph*(btot/bpoloid+btot0/bpoloid0) - r2iav=r2iav+dlph* - . (1.0d0/(bpoloid*rpsim**2)+1.0d0/(bpoloid0*rpsim0**2)) + r2iav=r2iav+dlph*(1.0_wp_/(bpoloid*rpsim**2) + . +1.0_wp_/(bpoloid0*rpsim0**2)) ajphiav=ajphiav+dlph* . (ajphi0/(bpoloid0*rpsim0)+ajphi/(bpoloid*rpsim)) @@ -2976,7 +3213,7 @@ c currp = plasma current within psi=const bbav=bbav/anorm r2iav=r2iav/anorm - dvdpsi=2.0d0*pi*anorm + dvdpsi=2.0_wp_*pi*anorm riav=dadpsi/anorm b2av=b2av/anorm vcurrp(jp)=ccj*currp @@ -3003,14 +3240,14 @@ c ratio_pltor = Jcd_||/J_phi Jcd_|| = vratjpl(jp)=ratio_pltor vratja(jp)=ratio_cdator vratjb(jp)=ratio_cdbtor - qq=abs(dvdpsi*fpolv*r2iav/(4.0d0*pi*pi)) + qq=abs(dvdpsi*fpolv*r2iav/(4.0_wp_*pi*pi)) qqv(jp)=qq rn=frhotor_an(sqrt(pstab(jp))) qqan=q0+(qa-q0)*rn**alq - dadr=2.0d0*pi*rn*rpam*rpam - dvdr=4.0d0*pi*pi*rn*rmaxis*rpam*rpam + dadr=2.0_wp_*pi*rn*rpam*rpam + dvdr=4.0_wp_*pi*pi*rn*rmaxis*rpam*rpam c dadrhotv(jp)=rpam*rpam*btaxis*rn/qqan*dadpsi c dvdrhotv(jp)=rpam*rpam*btaxis*rn/qqan*dvdpsi @@ -3031,37 +3268,37 @@ c computation of fraction of circulating/trapped fraction fc, ft c and of function H(lambda,rhop) c ffhlam = Bmn/Bmx/fc integral_lambda^1 dlam/ - fc=0.0d0 - shlam=0.0d0 + fc=0.0_wp_ + shlam=0.0_wp_ do l=nlam,1,-1 lam=alam(l) - srl=0.0d0 - rl2=1.0d0-lam*bv(1)/bmmx - rl0=0.d0 + srl=0.0_wp_ + rl2=1.0_wp_-lam*bv(1)/bmmx + rl0=0.0_wp_ if(rl2.gt.0) rl0=sqrt(rl2) do inc=1,npoints-1 - rl2=1.0d0-lam*bv(inc+1)/bmmx - rl=0.0d0 + rl2=1.0_wp_-lam*bv(inc+1)/bmmx + rl=0.0_wp_ if(rl2.gt.0) rl=sqrt(rl2) - srl=srl+0.5d0*dlpv(inc)*(rl/bpv(inc+1)+rl0/bpv(inc)) + srl=srl+0.5_wp_*dlpv(inc)*(rl/bpv(inc+1)+rl0/bpv(inc)) rl0=rl end do srl=srl/anorm - dhlam=0.5d0/srl + dhlam=0.5_wp_/srl fc=fc+lam/srl*weights(l) if(l.eq.nlam) then - fhlam(jp,l)=0.0d0 - ffhlam(nlam*(jp-1)+l)=0.0d0 + fhlam(jp,l)=0.0_wp_ + ffhlam(nlam*(jp-1)+l)=0.0_wp_ dffhlam(nlam*(jp-1)+l)=-dhlam dhlam0=dhlam else - shlam=shlam+0.5d0*(dhlam+dhlam0)*dlam + shlam=shlam+0.5_wp_*(dhlam+dhlam0)*dlam fhlam(jp,l)=shlam dffhlam(nlam*(jp-1)+l)=-dhlam dhlam0=dhlam end if end do - fc=0.75d0*b2av/bmmx**2*fc*dlam + fc=0.75_wp_*b2av/bmmx**2*fc*dlam ffc(jp)=fc ccfh=bmmn/bmmx/fc @@ -3079,9 +3316,9 @@ c ffhlam = Bmn/Bmx/fc integral_lambda^1 dlam/ do jp=1,npsi rhotqv(jp)=rhotqv(jp)/rhotqv(npsi) if(jp.eq.npsi) then - rhotqv(jp)=1.0d0 - rpstab(jp)=1.0d0 - pstab(jp)=1.0d0 + rhotqv(jp)=1.0_wp_ + rpstab(jp)=1.0_wp_ + pstab(jp)=1.0_wp_ end if rhot_eq(jp)=frhotor_an(sqrt(pstab(jp))) write(56,99) pstab(jp),rhot_eq(jp),rhotqv(jp), @@ -3126,7 +3363,7 @@ c spline coefficients of rhot c spline interpolation of H(lambda,rhop) and dH/dlambda iopt=0 - s=0.0d0 + s=0.0_wp_ call regrid(iopt,npsi,rpstab,nlam,alam,ffhlam, . zero,one,zero,one,ksp,ksp,s, . njest,nlest,njp,tjp,nlm,tlm,ch,fp,wrk,lwrk,iwrk,kwrk,ier) @@ -3135,63 +3372,68 @@ c spline interpolation of H(lambda,rhop) and dH/dlambda call coeff_parder(tjp,njp,tlm,nlm,ch,ksp,ksp,0,1, . ch01,lw01,ier) - - +c return 99 format(20(1x,e12.5)) end - - +c +c +c subroutine rhopol_an use dierckx, only : curfit,splev + use const_and_precisions, only : wp_ use graydata_par, only : sgniphi use graydata_anequil, only : q0,qa,alq,b0,rr0m,rpam use interp_eqprof, only : psia - - implicit real*8(a-h,o-z) - parameter(nnr=101,nrest=nnr+4) - parameter(lwrk=nnr*4+nrest*16) - dimension rhop(nnr),rhot(nnr),rhopi(nnr),rhoti(nnr) - dimension psin(nnr) - dimension trp(nrest),crp(nrest) - dimension trot(nrest),crot(nrest) - dimension wrk(lwrk),iwrk(nrest),wp(nrest) + implicit none +c local constants + integer, parameter :: nnr=101,nrest=nnr+4,lwrk=nnr*4+nrest*16 +c local variables + integer :: i,ier,iopt,kspl + integer, dimension(nrest) :: iwrk + real(wp_) :: dr,xb,xe,ss,rp,fq0,fq1,qq,res,rn + real(wp_), dimension(nnr) :: rhop,rhot,rhopi,rhoti,psin + real(wp_), dimension(lwrk) :: wrk + real(wp_), dimension(nrest) :: wp +c common/external functions/variables + integer :: nsrp,nsrot + real(wp_), dimension(nrest) :: trp,crp,trot,crot +c common/coffrtp/trp common/coffrn/nsrp common/coffrp/crp common/coffrptt/trot common/coffrnt/nsrot common/coffrpt/crot - - - dr=1.0d0/dble(nnr-1) - rhot(1)=0.0d0 - psin(1)=0.0d0 - res=0.0d0 - fq0=0.0d0 +c + dr=1.0_wp_/dble(nnr-1) + rhot(1)=0.0_wp_ + psin(1)=0.0_wp_ + res=0.0_wp_ + fq0=0.0_wp_ do i=2,nnr rhot(i)=(i-1)*dr rn=rhot(i) qq=q0+(qa-q0)*rn**alq fq1=rn/qq - res=res+0.5d0*(fq1+fq0)*dr + res=res+0.5_wp_*(fq1+fq0)*dr fq0=fq1 psin(i)=res end do - wp=1.0d0 + wp=1.0_wp_ psin=psin/res rhop=sqrt(psin) psia=-sgniphi*abs(res*rpam*rpam*b0) - print*,psia,log(8.0d0*rr0m/rpam)-2.0d0 - wp(1)=1.0d3 - wp(nnr)=1.0d3 + print*,psia,log(8.0_wp_*rr0m/rpam)-2.0_wp_ + wp(1)=1.0e3_wp_ + wp(nnr)=1.0e3_wp_ c spline interpolation of rhopol versus rhotor iopt=0 - xb=0.0d0 - xe=1.0d0 - ss=0.00001d0 + 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,wrk,lwrk,iwrk,ier) @@ -3199,9 +3441,9 @@ c spline interpolation of rhopol versus rhotor c spline interpolation of rhotor versus rhopol iopt=0 - xb=0.0d0 - xe=1.0d0 - ss=0.00001d0 + xb=0.0_wp_ + xe=1.0_wp_ + ss=0.00001_wp_ kspl=3 call curfit(iopt,nnr,rhop,rhot,wp,xb,xe,kspl,ss,nrest,nsrot, . trot,crot,rp,wrk,lwrk,iwrk,ier) @@ -3213,65 +3455,84 @@ c spline interpolation of rhotor versus rhopol return end - - +c +c +c function frhotor_an(rhop) use dierckx, only : splev - implicit real*8(a-h,o-z) - parameter(nnr=101,nrest=nnr+4) - dimension trot(nrest),crot(nrest),rrs(1),ffspl(1) + use const_and_precisions, only : wp_ + 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 - rrs(1)=rhop - call splev(trot,nsrot,crot,3,rrs,ffspl,1,ier) - frhotor_an=ffspl(1) +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 - - implicit real*8 (a-h,o-z) - parameter(pi=3.14159265358979d0) - parameter(mest=4,kspl=3) - dimension rcn(npoints),zcn(npoints) - + implicit none +c arguments + integer :: ipr + real(wp_) :: h + real(wp_), dimension(npoints) :: rcn,zcn +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)) do ic=1,npoints zcn(ic)=zr0m+rpam*rn*sin(th*(ic-1)) rcn(ic)=rr0m+rpam*rn*cos(th*(ic-1)) - - if (ipr.gt.0) then - write(71,111) ic,h,rcn(ic),zcn(ic) - end if +c + if (ipr.gt.0) write(71,111) ic,h,rcn(ic),zcn(ic) end do - write(71,*) - -111 format(i6,12(1x,e12.5)) + write(71,*) +c return +111 format(i6,12(1x,e12.5)) end - - - - +c +c +c subroutine vectinit + use const_and_precisions, only : wp_ use dispersion, only: expinit use graydata_flags, only : iwarm - - implicit real*8 (a-h,o-z) - parameter(jmx=31,kmx=36,nmx=8000) - dimension psjki(jmx,kmx,nmx) - dimension tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx) - dimension pdjki(jmx,kmx,nmx),ppabs(jmx,kmx,nmx),tau1v(jmx,kmx) - dimension currj(jmx,kmx,nmx),didst(jmx,kmx,nmx),ccci(jmx,kmx,nmx) - dimension iiv(jmx,kmx),iop(jmx,kmx),iow(jmx,kmx),ihcd(jmx,kmx) - dimension istore(jmx,kmx),anwcl(3),xwcl(3) + implicit none +c local constants + integer, parameter :: jmx=31,kmx=36,nmx=8000 +c local variables + integer :: i,j,k +c common/external functions/variables + integer nrayr,nrayth,nstep + integer, dimension(jmx,kmx) :: iiv,iop,iow,ihcd,istore + real(wp_) :: jclosest + real(wp_), dimension(3) :: anwcl,xwcl + real(wp_), dimension(jmx,kmx) :: tau1v + real(wp_), dimension(jmx,kmx,nmx) :: psjki,tauv,alphav,pdjki, + . ppabs,currj,didst,ccci c common/iiv/iiv common/iov/iop,iow,ihcd,istore @@ -3289,26 +3550,26 @@ c if(nrayr.gt.jmx) nrayr=jmx if(nrayth.gt.kmx) nrayth=kmx jclosest=nrayr+1 - anwcl(1:3)=0.0d0 - xwcl(1:3)=0.0d0 + anwcl(1:3)=0.0_wp_ + xwcl(1:3)=0.0_wp_ c do i=1,nstep do k=1,nrayth do j=1,nrayr - psjki(j,k,i)=0.0d0 - tauv(j,k,i)=0.0d0 - alphav(j,k,i)=0.0d0 - pdjki(j,k,i)=0.0d0 - ppabs(j,k,i)=0.0d0 - didst(j,k,i)=0.0d0 - ccci(j,k,i)=0.0d0 - currj(j,k,i)=0.0d0 + psjki(j,k,i)=0.0_wp_ + tauv(j,k,i)=0.0_wp_ + alphav(j,k,i)=0.0_wp_ + pdjki(j,k,i)=0.0_wp_ + ppabs(j,k,i)=0.0_wp_ + didst(j,k,i)=0.0_wp_ + ccci(j,k,i)=0.0_wp_ + currj(j,k,i)=0.0_wp_ iiv(j,k)=1 iop(j,k)=0 iow(j,k)=0 ihcd(j,k)=1 istore(j,k)=0 - tau1v(j,k)=0.0d0 + tau1v(j,k)=0.0_wp_ end do end do end do @@ -3317,19 +3578,22 @@ c c return end - - - +c +c +c subroutine vectinit2 - implicit real*8 (a-h,o-z) - parameter(jmx=31,kmx=36,nmx=8000) - dimension psjki(jmx,kmx,nmx) - dimension tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx) - dimension pdjki(jmx,kmx,nmx),ppabs(jmx,kmx,nmx) - dimension currj(jmx,kmx,nmx),didst(jmx,kmx,nmx),ccci(jmx,kmx,nmx) - dimension iiv(jmx,kmx),iop(jmx,kmx),iow(jmx,kmx),ihcd(jmx,kmx) - dimension istore(jmx,kmx) - + use const_and_precisions, only : wp_ + implicit none +c local constants + integer, parameter :: jmx=31,kmx=36,nmx=8000 +c local variables + integer :: i,j,k +c common/external functions/variables + integer nrayr,nrayth,nstep + integer, dimension(jmx,kmx) :: iiv,iop,iow,ihcd,istore + real(wp_), dimension(jmx,kmx,nmx) :: psjki,tauv,alphav,pdjki, + . ppabs,currj,didst,ccci +c common/iiv/iiv common/iov/iop,iow,ihcd,istore common/psjki/psjki @@ -3339,18 +3603,18 @@ c common/dcjki/didst common/nray/nrayr,nrayth common/nstep/nstep - +c do i=1,nstep do k=1,nrayth do j=1,nrayr - psjki(j,k,i)=0.0d0 - tauv(j,k,i)=0.0d0 - alphav(j,k,i)=0.0d0 - pdjki(j,k,i)=0.0d0 - ppabs(j,k,i)=0.0d0 - didst(j,k,i)=0.0d0 - ccci(j,k,i)=0.0d0 - currj(j,k,i)=0.0d0 + psjki(j,k,i)=0.0_wp_ + tauv(j,k,i)=0.0_wp_ + alphav(j,k,i)=0.0_wp_ + pdjki(j,k,i)=0.0_wp_ + ppabs(j,k,i)=0.0_wp_ + didst(j,k,i)=0.0_wp_ + ccci(j,k,i)=0.0_wp_ + currj(j,k,i)=0.0_wp_ iiv(j,k)=1 iop(j,k)=0 iow(j,k)=0 @@ -3365,13 +3629,16 @@ c c c subroutine paraminit - implicit real*8(a-h,o-z) + use const_and_precisions, only : wp_ + implicit none +c common/external functions/variables + integer istep,istpr,istpl,ierr,istop c common/istep/istep common/istgr/istpr,istpl common/ierr/ierr common/istop/istop - +c istpr=0 istpl=1 ierr=0 @@ -3381,13 +3648,19 @@ c return end c +c c subroutine updatepos - implicit real*8 (a-h,o-z) - parameter(jmx=31,kmx=36) - dimension xc(3,jmx,kmx),xco(3,jmx,kmx) - dimension du1(3,jmx,kmx),du1o(3,jmx,kmx) - dimension ywrk(6,jmx,kmx),ypwrk(6,jmx,kmx) + use const_and_precisions, only : wp_ + implicit none +c local constants + integer, parameter :: jmx=31,kmx=36 +c local variables + integer j,k +c common/external functions/variables + integer :: nrayr,nrayth + real(wp_), dimension(3,jmx,kmx) :: xc,xco,du1,du1o + real(wp_), dimension(6,jmx,kmx) :: ywrk,ypwrk c common/nray/nrayr,nrayth common/grco/xco,du1o @@ -3423,17 +3696,22 @@ c c c subroutine gradi - implicit real*8 (a-h,o-z) - parameter(jmx=31,kmx=36) - dimension dffiu(jmx),ddffiu(jmx) - dimension xc(3,jmx,kmx),xco(3,jmx,kmx) - dimension du1(3,jmx,kmx),du1o(3,jmx,kmx) - dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx) - dimension dgrad2v(3,jmx,kmx) - dimension grad2(jmx,kmx) - dimension dxv1(3),dxv2(3),dxv3(3),dgu(3) - dimension dgg1(3),dgg2(3),dgg3(3) - dimension df1(3),df2(3),df3(3) + use const_and_precisions, only : wp_ + implicit none +c local constants + integer, parameter :: jmx=31,kmx=36 +c local variables + integer iv,j,jm,jp,k,km,kp + real(wp_) ux,uxx,uxy,uxz,uy,uyy,uyz,uz,uzz, + . dfu,dfuu,gx,gxx,gxy,gxz,gy,gyy,gyz,gz,gzz + real(wp_), dimension(3) :: dxv1,dxv2,dxv3,dgu, + . dgg1,dgg2,dgg3,df1,df2,df3 +c common/external functions/variables + integer :: nrayr,nrayth + real(wp_), dimension(jmx) :: dffiu,ddffiu + real(wp_), dimension(jmx,kmx) :: grad2 + real(wp_), dimension(3,jmx,kmx) :: xc,xco,du1,du1o,gri,dgrad2v + real(wp_), dimension(3,3,jmx,kmx) :: ggri c common/nray/nrayr,nrayth common/fi/dffiu,ddffiu @@ -3449,9 +3727,9 @@ c do k=1,nrayth do j=1,nrayr if(j.eq.1) then - gri(1,j,k)=0.0d0 - gri(2,j,k)=0.0d0 - gri(3,j,k)=0.0d0 + gri(1,j,k)=0.0_wp_ + gri(2,j,k)=0.0_wp_ + gri(3,j,k)=0.0_wp_ jp=j+1 km=k-1 if(k.eq.1) km=nrayth @@ -3542,9 +3820,9 @@ c uxx=dgg1(1) uyy=dgg2(2) uzz=dgg3(3) - uxy=(dgg1(2)+dgg2(1))/2.0d0 - uxz=(dgg1(3)+dgg3(1))/2.0d0 - uyz=(dgg2(3)+dgg3(2))/2.0d0 + uxy=(dgg1(2)+dgg2(1))/2.0_wp_ + uxz=(dgg1(3)+dgg3(1))/2.0_wp_ + uyz=(dgg2(3)+dgg3(2))/2.0_wp_ c c derivatives of S_I and Grad(S_I) c @@ -3572,24 +3850,30 @@ c c c derivatives of |Grad(S_I)|^2 c - dgrad2v(1,j,k)=2.0d0*(gx*gxx+gy*gxy+gz*gxz) - dgrad2v(2,j,k)=2.0d0*(gx*gxy+gy*gyy+gz*gyz) - dgrad2v(3,j,k)=2.0d0*(gx*gxz+gy*gyz+gz*gzz) + dgrad2v(1,j,k)=2.0_wp_*(gx*gxx+gy*gxy+gz*gxz) + dgrad2v(2,j,k)=2.0_wp_*(gx*gxy+gy*gyy+gz*gyz) + dgrad2v(3,j,k)=2.0_wp_*(gx*gxz+gy*gyz+gz*gzz) end do end do c return end c +c +c + subroutine solg0(dxv1,dxv2,dxv3,dgg) c solution of the linear system of 3 eqs : dgg . dxv = dff c input vectors : dxv1, dxv2, dxv3, dff c output vector : dgg -c c dff=(1,0,0) c - subroutine solg0(dxv1,dxv2,dxv3,dgg) - double precision denom,aa1,aa2,aa3 - double precision dxv1(3),dxv2(3),dxv3(3),dgg(3) + use const_and_precisions, only : wp_ + implicit none +c arguments + real(wp_), dimension(3) :: dxv1,dxv2,dxv3,dgg +c local variables + real(wp_) denom,aa1,aa2,aa3 +c aa1=(dxv2(2)*dxv3(3)-dxv3(2)*dxv2(3)) aa2=(dxv1(2)*dxv3(3)-dxv1(3)*dxv3(2)) aa3=(dxv1(2)*dxv2(3)-dxv1(3)*dxv2(2)) @@ -3597,16 +3881,23 @@ c dgg(1) = aa1/denom dgg(2) = -(dxv2(1)*dxv3(3)-dxv3(1)*dxv2(3))/denom dgg(3) = (dxv2(1)*dxv3(2)-dxv3(1)*dxv2(2))/denom +c return end c -c three rhs vectors df1, df2, df3 +c c subroutine solg3(dxv1,dxv2,dxv3,df1,df2,df3,dg1,dg2,dg3) - double precision denom,a11,a21,a31,a12,a22,a32,a13,a23,a33 - double precision dxv1(3),dxv2(3),dxv3(3) - double precision df1(3),df2(3),df3(3) - double precision dg1(3),dg2(3),dg3(3) +c three rhs vectors df1, df2, df3 +c + use const_and_precisions, only : wp_ + implicit none +c arguments + real(wp_), dimension(3) :: dxv1,dxv2,dxv3,df1,df2,df3, + . dg1,dg2,dg3 +c local variables + real(wp_) denom,a11,a21,a31,a12,a22,a32,a13,a23,a33 +c a11=(dxv2(2)*dxv3(3)-dxv3(2)*dxv2(3)) a21=(dxv1(2)*dxv3(3)-dxv1(3)*dxv3(2)) a31=(dxv1(2)*dxv2(3)-dxv1(3)*dxv2(2)) @@ -3626,25 +3917,35 @@ c dg3(1) = (df3(1)*a11-df3(2)*a21+df3(3)*a31)/denom dg3(2) = -(df3(1)*a12-df3(2)*a22+df3(3)*a32)/denom dg3(3) = (df3(1)*a13-df3(2)*a23+df3(3)*a33)/denom +c return end c -c Runge-Kutta integrator +c c subroutine rkint4 +c Runge-Kutta integrator +c + use const_and_precisions, only : wp_ use graydata_flags, only : dst,igrad - - implicit real*8 (a-h,o-z) - parameter(ndim=6,jmx=31,kmx=36) - dimension y(ndim),yy(ndim) - dimension fk1(ndim),fk2(ndim),fk3(ndim),fk4(ndim) - dimension dgr2(3),dgr(3),ddgr(3,3) - dimension ywrk(ndim,jmx,kmx),ypwrk(ndim,jmx,kmx) - dimension grad2(jmx,kmx),dgrad2v(3,jmx,kmx) - dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx) + implicit none +c parameter + integer, parameter :: ndim=6,jmx=31,kmx=36 +c local variables + integer ieq,iv,j,jv,k,kkk + real(wp_) h,hh,h6 + real(wp_), dimension(ndim) :: y,yy,fk1,fk2,fk3,fk4 +c common/external functions/variables + integer :: nrayr,nrayth + real(wp_) gr2 + real(wp_), dimension(3) :: dgr2,dgr + real(wp_), dimension(3,3) :: ddgr + real(wp_), dimension(jmx,kmx) :: grad2 + real(wp_), dimension(3,jmx,kmx) :: dgrad2v,gri + real(wp_), dimension(ndim,jmx,kmx) :: ywrk,ypwrk + real(wp_), dimension(3,3,jmx,kmx) :: ggri c common/nray/nrayr,nrayth -c common/wrk/ywrk,ypwrk common/grad2jk/grad2 common/dgrad2jk/dgrad2v @@ -3654,8 +3955,8 @@ c common/dgr/dgr2,dgr,ddgr c h=dst - hh=h*0.5d0 - h6=h/6.0d0 + hh=h*0.5_wp_ + h6=h/6.0_wp_ c do j=1,nrayr kkk=nrayth @@ -3688,7 +3989,7 @@ c c do ieq=1,ndim ywrk(ieq,j,k)=y(ieq) - . +h6*(fk1(ieq)+2.0d0*fk2(ieq)+2.0d0*fk3(ieq)+fk4(ieq)) + . +h6*(fk1(ieq)+2.0_wp_*fk2(ieq)+2.0_wp_*fk3(ieq)+fk4(ieq)) end do end do end do @@ -3703,15 +4004,24 @@ c c c subroutine gwork(j,k) + use const_and_precisions, only : wp_ use graydata_flags, only : igrad - - implicit real*8 (a-h,o-z) - parameter(ndim=6,jmx=31,kmx=36) - dimension yy(ndim),yyp(ndim) - dimension dgr2(3),dgr(3),ddgr(3,3) - dimension ywrk(ndim,jmx,kmx),ypwrk(ndim,jmx,kmx) - dimension grad2(jmx,kmx),dgrad2v(3,jmx,kmx) - dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx) + implicit none +c local constants + integer, parameter :: ndim=6,jmx=31,kmx=36 +c arguments + integer j,k +c local variables + integer :: ieq,iv,jv + real(wp_), dimension(ndim) :: yy,yyp +c common/external functions/variables + real(wp_) :: gr2 + real(wp_), dimension(3) :: dgr2,dgr + real(wp_), dimension(3,3) :: ddgr + real(wp_), dimension(jmx,kmx) :: grad2 + real(wp_), dimension(3,jmx,kmx) :: dgrad2v,gri + real(wp_), dimension(ndim,jmx,kmx) :: ywrk,ypwrk + real(wp_), dimension(3,3,jmx,kmx) :: ggri c common/wrk/ywrk,ypwrk common/grad2jk/grad2 @@ -3752,14 +4062,25 @@ c c c subroutine fwork(y,dery) + use const_and_precisions, only : wp_ use graydata_flags, only : idst,igrad - - implicit real*8 (a-h,o-z) - parameter(ndim=6) - dimension y(ndim),dery(ndim) - dimension xv(3),anv(3),vgv(3),bv(3),derbv(3,3),derxg(3),deryg(3) - dimension derdxv(3),danpldxv(3),derdnv(3) - dimension dgr2(3),dgr(3),ddgr(3,3),dbgr(3) + implicit none +c local constants + integer, parameter :: ndim=6 +c arguments + real(wp_), dimension(ndim) :: y,dery +c local variables + integer iv,jv + real(wp_) :: xx,yy,zz,yg2,anpl2,anpr2,del,dnl,duh,dan2sdnpl, + . dan2sdxg,dan2sdyg,ddelnpl2,ddelnpl2x,ddelnpl2y,denom,derdel, + . derdom,dfdiadnpl,dfdiadxg,dfdiadyg + real(wp_), dimension(3) :: vgv,derdxv,danpldxv,derdnv,dbgr +c common/external functions/variables + integer ierr + real(wp_) sox,gr2,dd,an2s,an2,fdia,bdotgr,ddi,ddr11,anpl, + . anpr,xg,yg,vgm,derdnm,dersdst + real(wp_), dimension(3) :: dgr2,dgr,xv,anv,bv,derxg,deryg + real(wp_), dimension(3,3) :: ddgr,derbv c common/gr/gr2 common/dgr/dgr2,dgr,ddgr @@ -3789,19 +4110,19 @@ c anv(2) = y(5) anv(3) = y(6) c -C rr=sqrt(xx**2+yy**2) -C phi=acos(xx/rr) -C if (yy.lt.0.0d0) then -C phi=-phi -C end if +c rr=sqrt(xx**2+yy**2) +c phi=acos(xx/rr) +c if (yy.lt.0.0_wp_) then +c phi=-phi +c end if c call plas_deriv(xx,yy,zz) c an2=anv(1)*anv(1)+anv(2)*anv(2)+anv(3)*anv(3) anpl=anv(1)*bv(1)+anv(2)*bv(2)+anv(3)*bv(3) c - if(abs(anpl).gt.0.99d0) then - if(abs(anpl).le.1.05d0) then + if(abs(anpl).gt.0.99_wp_) then + if(abs(anpl).le.1.05_wp_) then ierr=97 else ierr=98 @@ -3809,64 +4130,64 @@ c end if c anpl2=anpl*anpl - dnl=1.0d0-anpl2 + dnl=1.0_wp_-anpl2 anpr2=an2-anpl2 - anpr=0.0d0 - if(anpr2.gt.0.0d0) anpr=sqrt(anpr2) + anpr=0.0_wp_ + if(anpr2.gt.0.0_wp_) anpr=sqrt(anpr2) yg2=yg**2 c - an2s=1.0d0 - dan2sdxg=0.0d0 - dan2sdyg=0.0d0 - dan2sdnpl=0.0d0 - del=0.0d0 - fdia=0.0d0 - dfdiadnpl=0.0d0 - dfdiadxg=0.0d0 - dfdiadyg=0.0d0 + an2s=1.0_wp_ + dan2sdxg=0.0_wp_ + dan2sdyg=0.0_wp_ + dan2sdnpl=0.0_wp_ + del=0.0_wp_ + fdia=0.0_wp_ + dfdiadnpl=0.0_wp_ + dfdiadxg=0.0_wp_ + dfdiadyg=0.0_wp_ c - duh=1.0d0-xg-yg2 - if(xg.gt.0.0d0) then - del=sqrt(dnl*dnl+4.0d0*anpl*anpl*(1.0d0-xg)/yg2) - an2s=1.0d0-xg -xg*yg2*(1.0d0+anpl2+sox*del)/duh/2.0d0 + duh=1.0_wp_-xg-yg2 + if(xg.gt.0.0_wp_) then + del=sqrt(dnl*dnl+4.0_wp_*anpl*anpl*(1.0_wp_-xg)/yg2) + an2s=1.0_wp_-xg -xg*yg2*(1.0_wp_+anpl2+sox*del)/duh/2.0_wp_ c - dan2sdxg=-yg2*(1.0d0-yg2)*(1.0d0+anpl2+sox*del)/duh**2/2.0d0 - . +sox*xg*anpl2/(del*duh)-1.0d0 - dan2sdyg=-xg*yg*(1.0d0-xg)*(1.0d0+anpl2+sox*del)/duh**2 - . +2.0d0*sox*xg*(1.0d0-xg)*anpl2/(yg*del*duh) + dan2sdxg=-yg2*(1.0_wp_-yg2)*(1.0_wp_+anpl2+sox*del)/ + . duh**2/2.0_wp_+sox*xg*anpl2/(del*duh)-1.0_wp_ + dan2sdyg=-xg*yg*(1.0_wp_-xg)*(1.0_wp_+anpl2+sox*del)/duh**2 + . +2.0_wp_*sox*xg*(1.0_wp_-xg)*anpl2/(yg*del*duh) dan2sdnpl=-xg*yg2*anpl/duh - . -sox*xg*anpl*(2.0d0*(1.0d0-xg)-yg2*dnl)/(del*duh) + . -sox*xg*anpl*(2.0_wp_*(1.0_wp_-xg)-yg2*dnl)/(del*duh) c if(igrad.gt.0) then - ddelnpl2=2.0d0*(2.0d0*(1.0d0-xg)*(1.0d0+3.0d0*anpl2**2) - . -yg2*dnl**3)/yg2/del**3 - fdia=-xg*yg2*(1.0d0+sox*ddelnpl2/2.0d0)/duh - derdel=2.0d0*(1.0d0-xg)*anpl2*(1.0d0+3.0d0*anpl2**2) - . -dnl**2*(1.0d0+3.0d0*anpl2)*yg2 - derdel=4.0d0*derdel/(yg**5*del**5) - ddelnpl2y=2.0d0*(1.0d0-xg)*derdel + ddelnpl2=2.0_wp_*(2.0_wp_*(1.0_wp_-xg) + . *(1.0_wp_+3.0_wp_*anpl2**2)-yg2*dnl**3)/yg2/del**3 + fdia=-xg*yg2*(1.0_wp_+sox*ddelnpl2/2.0_wp_)/duh + derdel=2.0_wp_*(1.0_wp_-xg)*anpl2*(1.0_wp_+3.0_wp_*anpl2**2) + . -dnl**2*(1.0_wp_+3.0_wp_*anpl2)*yg2 + derdel=4.0_wp_*derdel/(yg**5*del**5) + ddelnpl2y=2.0_wp_*(1.0_wp_-xg)*derdel ddelnpl2x=yg*derdel - dfdiadnpl=24.0d0*sox*xg*(1.0d0-xg)*anpl*(1.0d0-anpl**4) + dfdiadnpl=24.0_wp_*sox*xg*(1.0_wp_-xg)*anpl*(1.0_wp_-anpl**4) . /(yg2*del**5) - dfdiadxg=-yg2*(1.0d0-yg2)/duh**2-sox*yg2*((1.0d0-yg2)*ddelnpl2 - . +xg*duh*ddelnpl2x)/(2.0d0*duh**2) - dfdiadyg=-2.0d0*yg*xg*(1.0d0-xg)/duh**2 - . -sox*xg*yg*(2.0d0*(1.0d0-xg)*ddelnpl2 - . +yg*duh*ddelnpl2y)/(2.0d0*duh**2) + dfdiadxg=-yg2*(1.0_wp_-yg2)/duh**2-sox*yg2*((1.0_wp_-yg2) + . *ddelnpl2+xg*duh*ddelnpl2x)/(2.0_wp_*duh**2) + dfdiadyg=-2.0_wp_*yg*xg*(1.0_wp_-xg)/duh**2 + . -sox*xg*yg*(2.0_wp_*(1.0_wp_-xg)*ddelnpl2 + . +yg*duh*ddelnpl2y)/(2.0_wp_*duh**2) c end if end if c - bdotgr=0.0d0 + bdotgr=0.0_wp_ do iv=1,3 bdotgr=bdotgr+bv(iv)*dgr(iv) - dbgr(iv)=0.0d0 + dbgr(iv)=0.0_wp_ do jv=1,3 dbgr(iv)=dgr(jv)*derbv(jv,iv)+bv(jv)*ddgr(jv,iv) end do end do c - derdnm=0.0d0 + derdnm=0.0_wp_ c do iv=1,3 danpldxv(iv)=anv(1)*derbv(1,iv)+anv(2)*derbv(2,iv) @@ -3875,21 +4196,21 @@ c . +deryg(iv)*dan2sdyg+danpldxv(iv)*dan2sdnpl) derdxv(iv)=derdxv(iv)-igrad*dgr2(iv) c - derdxv(iv)=derdxv(iv)+fdia*bdotgr*dbgr(iv)+0.5d0*bdotgr**2* + derdxv(iv)=derdxv(iv)+fdia*bdotgr*dbgr(iv)+0.5_wp_*bdotgr**2* . (derxg(iv)*dfdiadxg+deryg(iv)*dfdiadyg+danpldxv(iv)*dfdiadnpl) c - derdnv(iv)=2.0d0*anv(iv)-bv(iv)*dan2sdnpl - . +0.5d0*bdotgr**2*bv(iv)*dfdiadnpl + derdnv(iv)=2.0_wp_*anv(iv)-bv(iv)*dan2sdnpl + . +0.5_wp_*bdotgr**2*bv(iv)*dfdiadnpl c derdnm=derdnm+derdnv(iv)**2 end do c derdnm=sqrt(derdnm) c - derdom=-2.0d0*an2+2.0d0*xg*dan2sdxg+yg*dan2sdyg+anpl*dan2sdnpl - . +2.0d0*igrad*gr2 - . -bdotgr**2*(fdia+xg*dfdiadxg+yg*dfdiadyg/2.0d0 - . +anpl*dfdiadnpl/2.0d0) + derdom=-2.0_wp_*an2+2.0_wp_*xg*dan2sdxg+yg*dan2sdyg+anpl*dan2sdnpl + . +2.0_wp_*igrad*gr2 + . -bdotgr**2*(fdia+xg*dfdiadxg+yg*dfdiadyg/2.0_wp_ + . +anpl*dfdiadnpl/2.0_wp_) c if (idst.eq.0) then c integration variable: s @@ -3925,7 +4246,7 @@ c c dd : dispersion relation (real part) c ddi : dispersion relation (imaginary part) c - dd=an2-an2s-igrad*(gr2-0.5d0*bdotgr**2*fdia) + dd=an2-an2s-igrad*(gr2-0.5_wp_*bdotgr**2*fdia) ddi=derdnv(1)*dgr(1)+derdnv(2)*dgr(2)+derdnv(3)*dgr(3) c return @@ -3934,13 +4255,23 @@ c c c subroutine plas_deriv(xx,yy,zz) - use const_and_precisions, only : pi + use const_and_precisions, only : wp_,pi,ccj=>mu0inv use graydata_flags, only : iequil use graydata_par, only : sgnbphi - - implicit real*8 (a-h,o-z) - dimension bv(3),derxg(3),deryg(3),derbv(3,3),dbtot(3) - dimension bvc(3),dbvcdc(3,3),dbvdc(3,3),dbv(3,3) + implicit none +c arguments + real(wp_) :: xx,yy,zz +c local variables + integer :: iv,jv + real(wp_) :: b2tot,csphi,dfpolv,drrdx,drrdy,dphidx,dphidy, + . rr,rr2,rrm,snphi,zzm + real(wp_), dimension(3) :: dbtot,bvc + real(wp_), dimension(3,3) :: dbvcdc,dbvdc,dbv +c common/external functions/variables + real(wp_) bres,brr,bphi,bzz,ajphi,btot,xg,yg,dxgdpsi,dpsidr, + . dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv,psinv + real(wp_), dimension(3) :: bv,derxg,deryg + real(wp_), dimension(3,3) :: derbv c common/parbres/bres common/parpl/brr,bphi,bzz,ajphi @@ -3956,20 +4287,20 @@ c common/fpol/fpolv,ffpv common/psival/psinv - xg=0.0d0 - yg=9.9d1 + xg=0.0_wp_ + yg=9.9e1_wp_ c do iv=1,3 - derxg(iv)=0.0d0 - deryg(iv)=0.0d0 - bv(iv)=0.0d0 - dbtot(iv)=0.0d0 + derxg(iv)=0.0_wp_ + deryg(iv)=0.0_wp_ + bv(iv)=0.0_wp_ + dbtot(iv)=0.0_wp_ do jv=1,3 - dbv(iv,jv)=0.0d0 - derbv(iv,jv)=0.0d0 - dbvcdc(iv,jv)=0.0d0 - dbvcdc(iv,jv)=0.0d0 - dbvdc(iv,jv)=0.0d0 + dbv(iv,jv)=0.0_wp_ + derbv(iv,jv)=0.0_wp_ + dbvcdc(iv,jv)=0.0_wp_ + dbvcdc(iv,jv)=0.0_wp_ + dbvdc(iv,jv)=0.0_wp_ end do end do c @@ -3987,8 +4318,8 @@ c c c convert from cm to meters c - zzm=1.0d-2*zz - rrm=1.0d-2*rr + zzm=1.0e-2_wp_*zz + rrm=1.0e-2_wp_*rr c if(iequil.eq.1) then call equian(rrm,zzm) @@ -4004,7 +4335,7 @@ c yg=fpolv/(rrm*bres) bphi=fpolv/rrm btot=abs(bphi) - if(psinv.lt.0.0d0) return + if(psinv.lt.0.0_wp_) return c c B = f(psi)/R e_phi+ grad psi x e_phi/R c @@ -4078,21 +4409,20 @@ c c bv(i) = B_i / B ; derbv(i,j) = d (B_i / B) /d x,y,z c do iv=1,3 - deryg(iv)=1.0d-2*dbtot(iv)/Bres + deryg(iv)=1.0e-2_wp_*dbtot(iv)/Bres bv(iv)=bv(iv)/btot do jv=1,3 - derbv(iv,jv)=1.0d-2*(dbv(iv,jv)-bv(iv)*dbtot(jv))/btot + derbv(iv,jv)=1.0e-2_wp_*(dbv(iv,jv)-bv(iv)*dbtot(jv))/btot end do end do c - derxg(1)=1.0d-2*drrdx*dpsidr*dxgdpsi - derxg(2)=1.0d-2*drrdy*dpsidr*dxgdpsi - derxg(3)=1.0d-2*dpsidz*dxgdpsi + derxg(1)=1.0e-2_wp_*drrdx*dpsidr*dxgdpsi + derxg(2)=1.0e-2_wp_*drrdy*dpsidr*dxgdpsi + derxg(3)=1.0e-2_wp_*dpsidz*dxgdpsi c c current density computation in Ampere/m^2 c ccj==1/mu_0 c - ccj=1.0d+7/(4.0d0*pi) ajphi=ccj*(dbvcdc(1,3)-dbvcdc(3,1)) c ajr=ccj*(dbvcdc(3,2)/rrm-dbvcdc(2,3)) c ajz=ccj*(bvc(2)/rrm+dbvcdc(2,1)-dbvcdc(1,2)) @@ -4103,11 +4433,20 @@ c c c subroutine equian(rrm,zzm) + 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 - - implicit real*8 (a-h,o-z) + implicit none +c arguments + real(wp_) :: rrm,zzm +c local variables + real(wp_) :: cst,dpsidrp,d2psidrp,dqq,qq,rn,rpm,sgn,snt,dfpolv, + . rhop,rhot +c common/external functions/variables + real(wp_) :: psinv,dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv, + . ffpv,xg,xgcn,dxgdpsi,dens,ddens + real(wp_) :: frhopol c common/psival/psinv common/derip1/dpsidr,dpsidz @@ -4117,9 +4456,9 @@ c common/dxgdps/dxgdpsi common/xgcn/xgcn common/dens/dens,ddens - - dpsidrp=0.0d0 - d2psidrp=0.0d0 +c + dpsidrp=0.0_wp_ + d2psidrp=0.0_wp_ c c simple model for equilibrium: large aspect ratio c outside plasma: analytical continuation, not solution Maxwell eqs @@ -4127,9 +4466,9 @@ c rpm=sqrt((rrm-rr0m)**2+(zzm-zr0m)**2) rn=rpm/rpam c - snt=0.0d0 - cst=1.0d0 - if (rpm.gt.0.0d0) then + snt=0.0_wp_ + cst=1.0_wp_ + if (rpm.gt.0.0_wp_) then snt=(zzm-zr0m)/rpm cst=(rrm-rr0m)/rpm end if @@ -4139,24 +4478,24 @@ c rhop=frhopol(rhot) psinv=rhop*rhop else - psinv=1.0d0+B0*rpam**2/2.0d0/abs(psia)/qa*(rn*rn-1.0d0) + psinv=1.0_wp_+B0*rpam**2/2.0_wp_/abs(psia)/qa*(rn*rn-1.0_wp_) rhop=sqrt(psinv) end if c sgn=sgniphi*sgnbphi - if(rn.le.1.0d0) then + if(rn.le.1.0_wp_) then qq=q0+(qa-q0)*rn**alq dpsidrp=B0*rpam*rn/qq*sgn - dqq=alq*(qa-q0)*rn**(alq-1.0d0) - d2psidrp=B0*(1.0d0-rn*dqq/qq)/qq*sgn + dqq=alq*(qa-q0)*rn**(alq-1.0_wp_) + d2psidrp=B0*(1.0_wp_-rn*dqq/qq)/qq*sgn else dpsidrp=B0*rpam/qa*rn*sgn d2psidrp=B0/qa*sgn end if c fpolv=sgnbphi*b0*rr0m - dfpolv=0.0d0 - ffpv=0.0d0 + dfpolv=0.0_wp_ + ffpv=0.0_wp_ c dpsidr=dpsidrp*cst dpsidz=dpsidrp*snt @@ -4166,22 +4505,30 @@ c return end - - +c +c +c subroutine equinum_psi(rpsim,zpsim) + use const_and_precisions, only : wp_ use interp_eqprof, only : rmnm,rmxm,zmnm,zmxm,psiant,psinop, . tr,tz,ccspl=>cceq,nsrt,nszt use dierckx, only : fpbisp - - implicit real*8 (a-h,o-z) - parameter(lwrk=8,liwrk=2) - dimension wrk(lwrk),iwrk(liwrk) - dimension rrs(1),zzs(1),ffspl(1) - parameter(nrs=1,nzs=1) + implicit none +c local constants + integer, parameter :: lwrk=8,liwrk=2,nrs=1,nzs=1 +c arguments + real(wp_) :: rpsim,zpsim +c local variables + integer :: nsr,nsz + integer, dimension(liwrk) :: iwrk + real(wp_), dimension(1) :: rrs,zzs,ffspl + real(wp_), dimension(lwrk) :: wrk +c common/external functions/variables + real(wp_) :: psinv c common/psival/psinv c - psinv=-1.0d0 + psinv=-1.0_wp_ c c here lengths are measured in meters c @@ -4198,13 +4545,22 @@ c c return end - +c +c +c subroutine equinum_derpsi(rpsim,zpsim,iderpsi) + use const_and_precisions, only : wp_ use interp_eqprof, only : nr,nz,rmnm,rmxm,zmnm,zmxm,cc01=>cceq01, . cc10=>cceq10,cc20=>cceq20,cc02=>cceq02,cc11=>cceq11 - - implicit real*8 (a-h,o-z) - integer*4 iderpsi + implicit none +c arguments + integer :: iderpsi + real(wp_) :: rpsim,zpsim +c local variables + integer :: lw10,lw01,lw20,lw02,lw11,nur,nuz + real(wp_) :: derpsi +c common/external functions/variables + real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz c common/derip1/dpsidr,dpsidz common/derip2/ddpsidrr,ddpsidzz,ddpsidrz @@ -4215,18 +4571,18 @@ c lw02=nr*4+nz*2+nr*nz lw11=nr*3+nz*3+nr*nz c - dpsidr=0.0d0 - dpsidz=0.0d0 - ddpsidrr=0.0d0 - ddpsidzz=0.0d0 - ddpsidrz=0.0d0 + dpsidr=0.0_wp_ + dpsidz=0.0_wp_ + ddpsidrr=0.0_wp_ + ddpsidzz=0.0_wp_ + ddpsidrz=0.0_wp_ c here lengths are measured in meters if (rpsim.gt.rmxm.or.rpsim.lt.rmnm) return if (zpsim.gt.zmxm.or.zpsim.lt.zmnm) return - +c select case(iderpsi) - + case(1) nur=1 nuz=0 @@ -4236,7 +4592,7 @@ c here lengths are measured in meters nuz=1 call sub_derpsi(rpsim,zpsim,nur,nuz,derpsi,cc01,lw01) dpsidz=derpsi - + case(2) nur=2 nuz=0 @@ -4272,26 +4628,33 @@ c here lengths are measured in meters nuz=1 call sub_derpsi(rpsim,zpsim,nur,nuz,derpsi,cc11,lw11) ddpsidrz=derpsi - + case default print*,'iderpsi undefined' end select - +c return end - +c +c +c subroutine sub_derpsi(rpsim,zpsim,nur,nuz,derpsi,cc,lw) + use const_and_precisions, only : wp_ use dierckx, only : fpbisp use interp_eqprof, only : nr,nz,psia,tr,tz,nsrt,nszt - - implicit real*8 (a-h,o-z) - parameter(liwrk=2) - parameter(nrs=1,nzs=1) - dimension iwrk(liwrk) - dimension rrs(1),zzs(1),ffspl(1) - dimension cc(lw) - + implicit none +c local constants + integer, parameter :: liwrk=2,nrs=1,nzs=1 +c arguments + integer :: nur,nuz,lw + real(wp_) :: rpsim,zpsim,derpsi + real(wp_), dimension(lw) :: cc +c local variables + integer :: iwr,iwz,kkr,kkz,nsr,nsz + integer, dimension(liwrk) :: iwrk + real(wp_), dimension(1) :: rrs,zzs,ffspl +c rrs(1)=rpsim zzs(1)=zpsim nsr=nsrt @@ -4303,28 +4666,36 @@ c here lengths are measured in meters call fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc,kkr,kkz . ,rrs,nrs,zzs,nzs,ffspl,cc(iwr),cc(iwz),iwrk(1),iwrk(2)) derpsi=ffspl(1)*psia - +c return end - +c +c +c subroutine equinum_fpol(iderfpol) + use const_and_precisions, only : wp_ use dierckx, only : splev,splder use interp_eqprof, only : nr,psia,tfp,nsft,cfp,fpolas - - implicit real*8 (a-h,o-z) - dimension rrs(1),ffspl(1) - dimension wrkfd(nr+4) - integer*4 iderfpol + implicit none +c arguments + integer :: iderfpol +c local variables + integer :: ier + real(wp_) :: dfpolv + real(wp_), dimension(1) :: rrs,ffspl + real(wp_), dimension(nr+4) :: wrkfd +c common/external functions/variables + real(wp_) :: psinv,fpolv,ffpv c common/psival/psinv common/fpol/fpolv,ffpv - +c fpolv=fpolas - dfpolv=0.0d0 - ffpv=0.0d0 + dfpolv=0.0_wp_ + ffpv=0.0_wp_ if(iderfpol.lt.0.or.iderfpol.gt.1) return - - if(psinv.le.1.0d0.and.psinv.gt.0.0d0) then +c + if(psinv.le.1.0_wp_.and.psinv.gt.0.0_wp_) then rrs(1)=psinv call splev(tfp,nsft,cfp,3,rrs,ffspl,1,ier) fpolv=ffspl(1) @@ -4334,67 +4705,121 @@ c ffpv=fpolv*dfpolv/psia end if end if - - return - end - - subroutine bfield(rpsim,zpsim,bphi,brr,bzz) - implicit real*8 (a-h,o-z) - call btor(rpsim,zpsim,bphi) - call bpol(rpsim,zpsim,brr,bzz) - return - end - - subroutine btor(rpsim,zpsim,bphi) - implicit real*8 (a-h,o-z) - common/psival/psinv - common/fpol/fpolv,ffpv - call equinum_psi(rpsim,zpsim) - call equinum_fpol(0) - bphi=fpolv/rpsim - return - end - - subroutine bpol(rpsim,zpsim,brr,bzz) - implicit real*8 (a-h,o-z) - common/derip1/dpsidr,dpsidz - call equinum_derpsi(rpsim,zpsim,1) - brr=-dpsidz/rpsim - bzz= dpsidr/rpsim - return - end - - subroutine tor_curr_psi(h,ajphi) - use interp_eqprof, only : zmaxis - - implicit real*8 (a-h,o-z) - - call psi_raxis(h,r1,r2) - call tor_curr(r2,zmaxis,ajphi) +c return end +c +c +c + subroutine bfield(rpsim,zpsim,bphi,brr,bzz) + use const_and_precisions, only : wp_ + implicit none +c arguments + real(wp_) :: rpsim,zpsim,bphi,brr,bzz +c + call btor(rpsim,zpsim,bphi) + call bpol(rpsim,zpsim,brr,bzz) +c + return + end +c +c +c + subroutine btor(rpsim,zpsim,bphi) + use const_and_precisions, only : wp_ + implicit none +c arguments + real(wp_) :: rpsim,zpsim,bphi +c common/external functions/variables + real(wp_) :: psinv,fpolv,ffpv +c + common/psival/psinv + common/fpol/fpolv,ffpv +c + call equinum_psi(rpsim,zpsim) + call equinum_fpol(0) + bphi=fpolv/rpsim +c + return + end +c +c +c + subroutine bpol(rpsim,zpsim,brr,bzz) + use const_and_precisions, only : wp_ + implicit none +c arguments + real(wp_) :: rpsim,zpsim,brr,bzz +c common/external functions/variables + real(wp_) :: dpsidr,dpsidz +c + common/derip1/dpsidr,dpsidz +c + call equinum_derpsi(rpsim,zpsim,1) + brr=-dpsidz/rpsim + bzz= dpsidr/rpsim +c + return + end +c +c +c + subroutine tor_curr_psi(h,ajphi) + use const_and_precisions, only : wp_ + use interp_eqprof, only : zmaxis + implicit none +c arguments + real(wp_) :: h,ajphi +c local variables + real(wp_) :: r1,r2 +c + call psi_raxis(h,r1,r2) + call tor_curr(r2,zmaxis,ajphi) +c + return + end +c +c c subroutine tor_curr(rpsim,zpsim,ajphi) - use const_and_precisions, only : pi,ccj=>mu0inv - implicit real*8 (a-h,o-z) + use const_and_precisions, only : wp_,ccj=>mu0inv + implicit none +c arguments + real(wp_) :: rpsim,zpsim,ajphi +c local variables + real(wp_) :: bzz,dbvcdc13,dbvcdc31 +c common/external functions/variables + real(wp_) dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz +c common/derip1/dpsidr,dpsidz common/derip2/ddpsidrr,ddpsidzz,ddpsidrz +c call equinum_derpsi(rpsim,zpsim,3) bzz= dpsidr/rpsim dbvcdc13=-ddpsidzz/rpsim dbvcdc31= ddpsidrr/rpsim-bzz/rpsim ajphi=ccj*(dbvcdc13-dbvcdc31) +c return end +c +c c subroutine psi_raxis(h,r1,r2) + use const_and_precisions, only : wp_ use interp_eqprof, only : psiant,psinop,zmaxis,nsr=>nsrt, . nsz=>nszt,cc=>cceq,tr,tz,nr use dierckx, only : profil,sproota - implicit real*8 (a-h,o-z) - parameter(mest=4,kspl=3) - dimension zeroc(mest) - real*8, dimension(:), allocatable :: czc + implicit none +c local constants + integer, parameter :: mest=4,kspl=3 +c arguments + real(wp_) :: h,r1,r2 +c local variables + integer :: iopt,ier,m + real(wp_) :: zc,val + real(wp_), dimension(mest) :: zeroc + real(wp_), dimension(:), allocatable :: czc c if(allocated(czc)) deallocate(czc) allocate(czc(nr+4)) @@ -4407,61 +4832,73 @@ c call sproota(val,tr,nsr,czc,zeroc,mest,m,ier) r1=zeroc(1) r2=zeroc(2) - +c deallocate(czc) - +c return end c c c subroutine sub_xg_derxg + use const_and_precisions, only : wp_ use interp_eqprof, only : psia - - implicit real*8 (a-h,o-z) + implicit none +c common/external functions/variables + real(wp_) psinv,xg,xgcn,dxgdpsi,dens,ddenspsin +c common/psival/psinv common/xgxg/xg common/dxgdps/dxgdpsi common/xgcn/xgcn common/dens/dens,ddenspsin - xg=0.0d0 - dxgdpsi=0.0d0 +c + xg=0.0_wp_ + dxgdpsi=0.0_wp_ c if(psinv.le.psdbnd.and.psinv.ge.0) then call density(psinv) xg=xgcn*dens dxgdpsi=xgcn*ddenspsin/psia c end if +c return end c c c subroutine density(arg) + use const_and_precisions, only : wp_ use graydata_flags, only : iprof use graydata_par, only : psdbnd use graydata_anequil, only : dens0,aln1,aln2 use interp_eqprof, only : psnpp,aad=>aa,bbd=>bb,ccd=>cc,tfn,nsfd, . cfn,npp use dierckx, only : splev,splder - - implicit real*8 (a-h,o-z) - dimension xxs(1),ffs(1) - dimension wrkfd(npp+4) + implicit none +c arguments + real(wp_) :: arg +c local variables + integer :: ier,nn,nn1,nn2,nu + real(wp_) :: profd,dprofd,dpsib + real(wp_), dimension(1) :: xxs,ffs + real(wp_), dimension(npp+4) :: wrkfd +c common/external functions/variables + real(wp_) :: dens,ddens c common/dens/dens,ddens c c computation of density [10^19 m^-3] and derivative wrt psi c - dens=0.0d0 - ddens=0.0d0 - if(arg.gt.psdbnd.or.arg.lt.0.0d0) return + dens=0.0_wp_ + ddens=0.0_wp_ + if(arg.gt.psdbnd.or.arg.lt.0.0_wp_) return c if(iprof.eq.0) then - if(arg.gt.1.0d0) return - profd=(1.0d0-arg**aln1)**aln2 + if(arg.gt.1.0_wp_) return + profd=(1.0_wp_-arg**aln1)**aln2 dens=dens0*profd - dprofd=-aln1*aln2*arg**(aln1-1.0d0) - . *(1.0d0-arg**aln1)**(aln2-1.0d0) + dprofd=-aln1*aln2*arg**(aln1-1.0_wp_) + . *(1.0_wp_-arg**aln1)**(aln2-1.0_wp_) ddens=dens0*dprofd else if(arg.le.psdbnd.and.arg.gt.psnpp) then @@ -4485,29 +4922,35 @@ c 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.0d-10) dens=0.0d0 + if(abs(dens).lt.1.0e-10_wp_) dens=0.0_wp_ end if - if(dens.lt.0.0d0) print*,' DENSITY NEGATIVE',dens + if(dens.lt.0.0_wp_) print*,' DENSITY NEGATIVE',dens c end if +c return end c c c function temp(arg) + use const_and_precisions, only : wp_ use graydata_flags, only : iprof use graydata_anequil, only : te0,dte0,alt1,alt2 use interp_eqprof, only : psrad,ct,npp use utils, only : locate use simplespline, only :spli - - implicit real*8 (a-h,o-z) + implicit none +c arguments + real(wp_) :: arg,temp +c local variables + integer :: k + real(wp_) :: proft,dps c - temp=0.0d0 - if(arg.ge.1.0d0.or.arg.lt.0.0d0) return + temp=0.0_wp_ + if(arg.ge.1.0_wp_.or.arg.lt.0.0_wp_) return if(iprof.eq.0) then - proft=(1.0d0-arg**alt1)**alt2 + proft=(1.0_wp_-arg**alt1)**alt2 temp=(te0-dte0)*proft+dte0 else call locate(psrad,npp,arg,k) @@ -4516,23 +4959,29 @@ c dps=arg-psrad(k) temp=spli(ct,npp,k,dps) endif +c return end c c c function fzeff(arg) + use const_and_precisions, only : wp_ use graydata_flags, only : iprof use graydata_anequil, only : zeffan use interp_eqprof, only : psrad,cz,npp use utils, only : locate use simplespline, only :spli - - implicit real*8 (a-h,o-z) + implicit none +c arguments + real(wp_) :: arg,fzeff +c local variables + integer :: k + real(wp_) :: ps,dps c fzeff=1 ps=arg - if(ps.gt.1.0d0.and.ps.lt.0.0d0) return + if(ps.gt.1.0_wp_.and.ps.lt.0.0_wp_) return if(iprof.eq.0) then fzeff=zeffan else @@ -4541,30 +4990,48 @@ c dps=ps-psrad(k) fzeff=spli(cz,npp,k,dps) endif +c return end c -c beam tracing initial conditions igrad=1 +c c subroutine ic_gb - use const_and_precisions, only : izero,zero,one,pi,cvdr=>degree +c beam tracing initial conditions igrad=1 +c + use const_and_precisions, only : wp_,izero,zero,one,pi, + . cvdr=>degree,ui=>im use math, only : catand use graydata_flags, only : ipol use graydata_par, only : rwmax,psipol0,chipol0 - - implicit real*8 (a-h,o-z) - parameter(ndim=6,ndimm=3) - parameter(jmx=31,kmx=36) - dimension ywrk0(ndim,jmx,kmx),ypwrk0(ndim,jmx,kmx) - dimension ytmp(ndim),yptmp(ndim) - dimension xc0(ndimm,jmx,kmx),du10(ndimm,jmx,kmx) - dimension dffiu(jmx),ddffiu(jmx) - dimension grad2(jmx,kmx),dgrad2v(ndimm,jmx,kmx) - dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx) - complex*16 ext(jmx,kmx,0:4),eyt(jmx,kmx,0:4) - complex*16 ui,sss,ddd,phic,qi1,qi2,tc,ts,qqxx,qqxy,qqyy - complex*16 dqi1,dqi2,dqqxx,dqqyy,dqqxy - complex*16 d2qi1,d2qi2,d2qqxx,d2qqyy,d2qqxy + implicit none +c local constants + integer, parameter :: ndim=6,ndimm=3,jmx=31,kmx=36 +c local variables + integer :: j,k,iproj,nfilp + real(wp_) :: csth,snth,csps,snps,phiwrad,phirrad,csphiw,snphiw, + . wwcsi,wweta,sk,sw,dk,dw,rci1,ww1,rci2,ww2,wwxx,wwyy,wwxy,rcixx, + . rciyy,rcixy,dwwxx,dwwyy,dwwxy,d2wwxx,d2wwyy,d2wwxy,drcixx,drciyy, + . drcixy,dr,da,ddfu,u,alfak,dcsiw,detaw,dx0t,dy0t,x0t,y0t,z0t,dx0, + . dy0,dz0,x0,y0,z0,gxt,gyt,gzt,gr2,gxxt,gyyt,gzzt,gxyt,gxzt,gyzt, + . dgr2xt,dgr2yt,dgr2zt,dgr2x,dgr2y,dgr2z,pppx,pppy,denpp,ppx,ppy, + . anzt,anxt,anyt,anx,any,anz,an20,an0,anx0,any0,anz0,du1tx,du1ty, + . du1tz,vgradi,r0,x0m,y0m,r0m,z0m,ddr110,deltapol,qq,uu,vv + real(wp_), dimension(ndim) :: ytmp,yptmp + complex(wp_) :: sss,ddd,phic,qi1,qi2,tc,ts,qqxx,qqxy,qqyy,dqi1, + . dqi2,dqqxx,dqqyy,dqqxy,d2qi1,d2qi2,d2qqxx,d2qqyy,d2qqxy +c common/external functions/variables + integer :: nrayr,nrayth + real(wp_) :: ak0,akinv,fhz,wcsi,weta,rcicsi,rcieta,phiw,phir, + . anx0c,any0c,anz0c,x00,y00,z00,dd,an2s,an2,fdia,bdotgr,ddi, + . ddr11,psinv,dens,ddens,tekev,anpl,anpr,brr,bphi,bzz,ajphi + real(wp_), dimension(jmx) :: dffiu,ddffiu + real(wp_), dimension(jmx,kmx) :: grad2 + real(wp_), dimension(3,jmx,kmx) :: gri + real(wp_), dimension(ndim,jmx,kmx) :: ywrk0,ypwrk0 + real(wp_), dimension(ndimm,jmx,kmx) :: xc0,du10,dgrad2v + real(wp_), dimension(3,3,jmx,kmx) :: ggri + complex(wp_), dimension(jmx,kmx,0:4) :: ext,eyt c common/nray/nrayr,nrayth common/parwv/ak0,akinv,fhz @@ -4585,14 +5052,12 @@ c common/parpl/brr,bphi,bzz,ajphi common/dens/dens,ddens common/tete/tekev - c - ui=(0.0d0,1.0d0) csth=anz0c - snth=sqrt(1.0d0-csth**2) - csps=1.0d0 - snps=0.0d0 - if(snth.gt.0.0d0) then + snth=sqrt(1.0_wp_-csth**2) + csps=1.0_wp_ + snps=0.0_wp_ + if(snth.gt.0.0_wp_) then csps=any0c/snth snps=anx0c/snth end if @@ -4604,21 +5069,21 @@ c c csphir=cos(phirrad) c snphir=sin(phirrad) c - wwcsi=2.0d0*akinv/wcsi**2 - wweta=2.0d0*akinv/weta**2 + wwcsi=2.0_wp_*akinv/wcsi**2 + wweta=2.0_wp_*akinv/weta**2 c if(phir.ne.phiw) then sk=(rcicsi+rcieta) sw=(wwcsi+wweta) dk=(rcicsi-rcieta) dw=(wwcsi-wweta) - ts=-(dk*sin(2.0d0*phirrad)-ui*dw*sin(2.0d0*phiwrad)) - tc=(dk*cos(2.0d0*phirrad)-ui*dw*cos(2.0d0*phiwrad)) - phic=0.5d0*catand(ts/tc) + ts=-(dk*sin(2.0_wp_*phirrad)-ui*dw*sin(2.0_wp_*phiwrad)) + tc=(dk*cos(2.0_wp_*phirrad)-ui*dw*cos(2.0_wp_*phiwrad)) + phic=0.5_wp_*catand(ts/tc) ddd=dk*cos(2*(phirrad+phic))-ui*dw*cos(2*(phiwrad+phic)) sss=sk-ui*sw - qi1=0.5d0*(sss+ddd) - qi2=0.5d0*(sss-ddd) + qi1=0.5_wp_*(sss+ddd) + qi2=0.5_wp_*(sss-ddd) rci1=dble(qi1) ww1=-dimag(qi1) rci2=dble(qi2) @@ -4633,20 +5098,20 @@ c qi2=rci2-ui*ww2 end if -c w01=sqrt(2.0d0*akinv/ww1) +c w01=sqrt(2.0_wp_*akinv/ww1) c z01=-rci1/(rci1**2+ww1**2) -c w02=sqrt(2.0d0*akinv/ww2) +c w02=sqrt(2.0_wp_*akinv/ww2) c z02=-rci2/(rci2**2+ww2**2) qqxx=qi1*cos(phic)**2+qi2*sin(phic)**2 qqyy=qi1*sin(phic)**2+qi2*cos(phic)**2 - qqxy=-(qi1-qi2)*sin(2.0d0*phic) + qqxy=-(qi1-qi2)*sin(2.0_wp_*phic) wwxx=-dimag(qqxx) wwyy=-dimag(qqyy) - wwxy=-dimag(qqxy)/2.0d0 + wwxy=-dimag(qqxy)/2.0_wp_ rcixx=dble(qqxx) rciyy=dble(qqyy) - rcixy=dble(qqxy)/2.0d0 + rcixy=dble(qqxy)/2.0_wp_ dqi1=-qi1**2 dqi2=-qi2**2 @@ -4654,28 +5119,28 @@ c z02=-rci2/(rci2**2+ww2**2) d2qi2=2*qi2**3 dqqxx=dqi1*cos(phic)**2+dqi2*sin(phic)**2 dqqyy=dqi1*sin(phic)**2+dqi2*cos(phic)**2 - dqqxy=-(dqi1-dqi2)*sin(2.0d0*phic) + dqqxy=-(dqi1-dqi2)*sin(2.0_wp_*phic) d2qqxx=d2qi1*cos(phic)**2+d2qi2*sin(phic)**2 d2qqyy=d2qi1*sin(phic)**2+d2qi2*cos(phic)**2 - d2qqxy=-(d2qi1-d2qi2)*sin(2.0d0*phic) + d2qqxy=-(d2qi1-d2qi2)*sin(2.0_wp_*phic) dwwxx=-dimag(dqqxx) dwwyy=-dimag(dqqyy) - dwwxy=-dimag(dqqxy)/2.0d0 + dwwxy=-dimag(dqqxy)/2.0_wp_ d2wwxx=-dimag(d2qqxx) d2wwyy=-dimag(d2qqyy) - d2wwxy=-dimag(d2qqxy)/2.0d0 + d2wwxy=-dimag(d2qqxy)/2.0_wp_ drcixx=dble(dqqxx) drciyy=dble(dqqyy) - drcixy=dble(dqqxy)/2.0d0 + drcixy=dble(dqqxy)/2.0_wp_ - dr=1.0d0 + dr=1.0_wp_ if(nrayr.gt.1) dr=rwmax/dble(nrayr-1) - da=2.0d0*pi/dble(nrayth) + da=2.0_wp_*pi/dble(nrayth) c - ddfu=2.0d0*dr**2*akinv + ddfu=2.0_wp_*dr**2*akinv do j=1,nrayr u=dble(j-1) -c ffi=u**2*ddfu/2.0d0 +c ffi=u**2*ddfu/2.0_wp_ dffiu(j)=u*ddfu ddffiu(j)=ddfu do k=1,nrayth @@ -4686,7 +5151,7 @@ c ffi=u**2*ddfu/2.0d0 dy0t=dcsiw*snphiw+detaw*csphiw x0t=u*dx0t y0t=u*dy0t - z0t=-0.5d0*(rcixx*x0t**2+rciyy*y0t**2+2*rcixy*x0t*y0t) + z0t=-0.5_wp_*(rcixx*x0t**2+rciyy*y0t**2+2*rcixy*x0t*y0t) c c csiw=u*dcsiw c etaw=u*detaw @@ -4702,58 +5167,61 @@ c gxt=x0t*wwxx+y0t*wwxy gyt=x0t*wwxy+y0t*wwyy - gzt=0.5d0*(x0t**2*dwwxx+y0t**2*dwwyy)+x0t*y0t*dwwxy + gzt=0.5_wp_*(x0t**2*dwwxx+y0t**2*dwwyy)+x0t*y0t*dwwxy gr2=gxt*gxt+gyt*gyt+gzt*gzt gxxt=wwxx gyyt=wwyy - gzzt=0.5d0*(x0t**2*d2wwxx+y0t**2*d2wwyy)+x0t*y0t*d2wwxy + gzzt=0.5_wp_*(x0t**2*d2wwxx+y0t**2*d2wwyy)+x0t*y0t*d2wwxy gxyt=wwxy gxzt=x0t*dwwxx+y0t*dwwxy gyzt=x0t*dwwxy+y0t*dwwyy - dgr2xt=2.0d0*(gxt*gxxt+gyt*gxyt+gzt*gxzt) - dgr2yt=2.0d0*(gxt*gxyt+gyt*gyyt+gzt*gyzt) - dgr2zt=2.0d0*(gxt*gxzt+gyt*gyzt+gzt*gzzt) + dgr2xt=2.0_wp_*(gxt*gxxt+gyt*gxyt+gzt*gxzt) + dgr2yt=2.0_wp_*(gxt*gxyt+gyt*gyyt+gzt*gyzt) + dgr2zt=2.0_wp_*(gxt*gxzt+gyt*gyzt+gzt*gzzt) dgr2x= dgr2xt*csps+snps*(dgr2yt*csth+dgr2zt*snth) dgr2y=-dgr2xt*snps+csps*(dgr2yt*csth+dgr2zt*snth) dgr2z= dgr2zt*csth-dgr2yt*snth gri(1,j,k)=gxt*csps+snps*(gyt*csth+gzt*snth) gri(2,j,k)=-gxt*snps+csps*(gyt*csth+gzt*snth) gri(3,j,k)=gzt*csth-gyt*snth - ggri(1,1,j,k)=gxxt*csps**2+ - . snps**2*(gyyt*csth**2+gzzt*snth**2+2.0d0*snth*csth*gyzt)+ - . 2.0d0*snps*csps*(gxyt*csth+gxzt*snth) - ggri(2,2,j,k)=gxxt*snps**2+ - . csps**2*(gyyt*csth**2+gzzt*snth**2+2.0d0*snth*csth*gyzt)- - . 2.0d0*snps*csps*(gxyt*csth+gxzt*snth) - ggri(3,3,j,k)=gzzt*csth**2+gyyt*snth**2-2.0d0*csth*snth*gyzt + ggri(1,1,j,k)=gxxt*csps**2+snps**2 + . *(gyyt*csth**2+gzzt*snth**2 + . +2.0_wp_*snth*csth*gyzt) + . +2.0_wp_*snps*csps*(gxyt*csth+gxzt*snth) + ggri(2,2,j,k)=gxxt*snps**2+csps**2 + . *(gyyt*csth**2+gzzt*snth**2 + . +2.0_wp_*snth*csth*gyzt) + . -2.0_wp_*snps*csps*(gxyt*csth+gxzt*snth) + ggri(3,3,j,k)=gzzt*csth**2+gyyt*snth**2 + . -2.0_wp_*csth*snth*gyzt ggri(1,2,j,k)=csps*snps*(-gxxt+csth**2*gyyt+snth**2*gzzt - . +2.0d0*csth*snth*gyzt) + . +2.0_wp_*csth*snth*gyzt) . +(csps**2-snps**2)*(snth*gxzt+csth*gxyt) - ggri(1,3,j,k)=csth*snth*snps*(gzzt-gyyt) - . +(csth**2-snth**2)*snps*gyzt+csps*(csth*gxzt-snth*gxyt) - ggri(2,3,j,k)=csth*snth*csps*(gzzt-gyyt) - . +(csth**2-snth**2)*csps*gyzt+snps*(snth*gxyt-csth*gxzt) + ggri(1,3,j,k)=csth*snth*snps*(gzzt-gyyt)+(csth**2-snth**2) + . *snps*gyzt+csps*(csth*gxzt-snth*gxyt) + ggri(2,3,j,k)=csth*snth*csps*(gzzt-gyyt)+(csth**2-snth**2) + . *csps*gyzt+snps*(snth*gxyt-csth*gxzt) ggri(2,1,j,k)=ggri(1,2,j,k) ggri(3,1,j,k)=ggri(1,3,j,k) ggri(3,2,j,k)=ggri(2,3,j,k) c du1tx=(dx0t*wwxx+dy0t*wwxy)/ddfu du1ty=(dx0t*wwxy+dy0t*wwyy)/ddfu - du1tz=0.5d0*u*(dx0t**2*dwwxx+dy0t**2*dwwyy+ - . 2.0d0*dx0t*dy0t*dwwxy)/ddfu + du1tz=0.5_wp_*u*(dx0t**2*dwwxx+dy0t**2*dwwyy + . +2.0_wp_*dx0t*dy0t*dwwxy)/ddfu c pppx=x0t*rcixx+y0t*rcixy pppy=x0t*rcixy+y0t*rciyy denpp=pppx*gxt+pppy*gyt - if (denpp.ne.0.0d0) then + if (denpp.ne.0.0_wp_) then ppx=-pppx*gzt/denpp ppy=-pppy*gzt/denpp else - ppx=0.0d0 - ppy=0.0d0 + ppx=0.0_wp_ + ppy=0.0_wp_ end if c - anzt=sqrt((1.0d0+gr2)/(1.0d0+ppx**2+ppy**2)) + anzt=sqrt((1.0_wp_+gr2)/(1.0_wp_+ppx**2+ppy**2)) anxt=ppx*anzt anyt=ppy*anzt c @@ -4761,7 +5229,7 @@ c any=-anxt*snps+csps*(anyt*csth+anzt*snth) anz= anzt*csth-anyt*snth c - an20=1.0d0+gr2 + an20=1.0_wp_+gr2 an0=sqrt(an20) anx0=anx any0=any @@ -4781,9 +5249,9 @@ c ypwrk0(1,j,k) = anx0/an0 ypwrk0(2,j,k) = any0/an0 ypwrk0(3,j,k) = anz0/an0 - ypwrk0(4,j,k) = dgr2x/an0/2.0d0 - ypwrk0(5,j,k) = dgr2y/an0/2.0d0 - ypwrk0(6,j,k) = dgr2z/an0/2.0d0 + ypwrk0(4,j,k) = dgr2x/an0/2.0_wp_ + ypwrk0(5,j,k) = dgr2y/an0/2.0_wp_ + ypwrk0(6,j,k) = dgr2z/an0/2.0_wp_ c ytmp=ywrk0(:,j,k) yptmp=ypwrk0(:,j,k) @@ -4792,20 +5260,20 @@ c if(ipol.eq.0) then call pol_limit(ext(j,k,0),eyt(j,k,0)) qq=abs(ext(j,k,0))**2-abs(eyt(j,k,0))**2 - uu=2.0d0*dble(ext(j,k,0)*dconjg(eyt(j,k,0))) - vv=2.0d0*dimag(ext(j,k,0)*dconjg(eyt(j,k,0))) + uu=2.0_wp_*dble(ext(j,k,0)*dconjg(eyt(j,k,0))) + vv=2.0_wp_*dimag(ext(j,k,0)*dconjg(eyt(j,k,0))) call polellipse(qq,uu,vv,psipol0,chipol0) else - qq=cos(2.0d0*chipol0*cvdr)*cos(2.0d0*psipol0*cvdr) - uu=cos(2.0d0*chipol0*cvdr)*sin(2.0d0*psipol0*cvdr) - vv=sin(2.0d0*chipol0*cvdr) + qq=cos(2.0_wp_*chipol0*cvdr)*cos(2.0_wp_*psipol0*cvdr) + uu=cos(2.0_wp_*chipol0*cvdr)*sin(2.0_wp_*psipol0*cvdr) + vv=sin(2.0_wp_*chipol0*cvdr) if(qq**2.lt.1) then - deltapol=asin(vv/sqrt(1.0d0-qq**2)) - ext(j,k,0)= sqrt((1.0d0+qq)/2) - eyt(j,k,0)= sqrt((1.0d0-qq)/2)*exp(-ui*deltapol) + deltapol=asin(vv/sqrt(1.0_wp_-qq**2)) + ext(j,k,0)= sqrt((1.0_wp_+qq)/2) + eyt(j,k,0)= sqrt((1.0_wp_-qq)/2)*exp(-ui*deltapol) else - ext(j,k,0)= 1.0d0 - eyt(j,k,0)= 0.0d0 + ext(j,k,0)= 1.0_wp_ + eyt(j,k,0)= 0.0_wp_ end if endif c @@ -4820,21 +5288,21 @@ c c dd=anx0**2+any0**2+anz0**2-an20 vgradi=anxt*gxt+anyt*gyt+anzt*gzt - ddi=2.0d0*vgradi + ddi=2.0_wp_*vgradi c r0=sqrt(x0**2+y0**2) - x0m=x0/1.0d2 - y0m=y0/1.0d2 - r0m=r0/1.0d2 - z0m=z0/1.0d2 + x0m=x0/1.0e2_wp_ + y0m=y0/1.0e2_wp_ + r0m=r0/1.0e2_wp_ + z0m=z0/1.0e2_wp_ if(j.eq.nrayr) then write(33,111) izero,j,k,zero,x0m,y0m,r0m,z0m, . psinv,zero,anpl,zero,one end if if(j.eq.1.and.k.eq.1) then - write(4,99) zero,r0m,z0m,atan2(y0m,x0m)*180.0d0/pi, + write(4,99) zero,r0m,z0m,atan2(y0m,x0m)*180.0_wp_/pi, . psinv,one,dens,tekev,brr,bphi,bzz, - . ajphi*1.0d-6,sqrt(anpl**2+anpr**2),anpl,zero, + . ajphi*1.0e-6_wp_,sqrt(anpl**2+anpr**2),anpl,zero, . zero,zero,zero,zero,zero,zero,zero,one ddr110=dd @@ -4858,24 +5326,37 @@ c 111 format(3i5,20(1x,e16.8e3)) end c -c ray tracing initial conditions igrad=0 +c c subroutine ic_rt - use const_and_precisions, only : izero,zero,one,pi,cvdr=>degree, - * ui=>im +c ray tracing initial conditions igrad=0 +c + use const_and_precisions, only : wp_,izero,zero,one,pi, + . cvdr=>degree,ui=>im use graydata_flags, only : ipol use graydata_par, only : rwmax,psipol0,chipol0 - - implicit real*8 (a-h,o-z) - parameter(ndim=6,ndimm=3) - parameter(jmx=31,kmx=36) - dimension ywrk0(ndim,jmx,kmx),ypwrk0(ndim,jmx,kmx) - dimension ytmp(ndim),yptmp(ndim) - dimension xc0(ndimm,jmx,kmx),du10(ndimm,jmx,kmx) - dimension dffiu(jmx),ddffiu(jmx) - dimension grad2(jmx,kmx),dgrad2v(ndimm,jmx,kmx) - dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx) - complex*16 ext(jmx,kmx,0:4),eyt(jmx,kmx,0:4) + implicit none +c local constants + integer, parameter :: ndim=6,ndimm=3,jmx=31,kmx=36 +c local variables + integer :: j,k,iv,jv,iproj,nfilp + real(wp_) :: csth,snth,csps,snps,phiwrad,csphiw,snphiw,dr,da,u, + . alfak,dcsiw,detaw,dx0t,dy0t,x0t,y0t,z0t,dx0,dy0,dz0,x0,y0,z0, + . anzt,anxt,anyt,anx,any,anz,an20,an0,anx0,any0,anz0,vgradi,r0, + . x0m,y0m,r0m,z0m,ancsi,aneta,ppcsi,ppeta,deltapol,qq,uu,vv + real(wp_), dimension(ndim) :: ytmp,yptmp +c common/external functions/variables + integer :: nrayr,nrayth + real(wp_) :: wcsi,weta,rcicsi,rcieta,phiw,phir,anx0c,any0c,anz0c, + . x00,y00,z00,dd,an2s,an2,fdia,bdotgr,ddi,ddr11,psinv,dens,ddens, + . tekev,anpl,anpr,brr,bphi,bzz,ajphi + real(wp_), dimension(jmx) :: dffiu,ddffiu + real(wp_), dimension(jmx,kmx) :: grad2 + real(wp_), dimension(3,jmx,kmx) :: gri + real(wp_), dimension(ndim,jmx,kmx) :: ywrk0,ypwrk0 + real(wp_), dimension(ndimm,jmx,kmx) :: xc0,du10,dgrad2v + real(wp_), dimension(3,3,jmx,kmx) :: ggri + complex(wp_), dimension(jmx,kmx,0:4) :: ext,eyt c common/nray/nrayr,nrayth common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir @@ -4897,10 +5378,10 @@ c common/tete/tekev c csth=anz0c - snth=sqrt(1.0d0-csth**2) - csps=1.0d0 - snps=0.0d0 - if(snth.gt.0.0d0) then + snth=sqrt(1.0_wp_-csth**2) + csps=1.0_wp_ + snps=0.0_wp_ + if(snth.gt.0.0_wp_) then csps=any0c/snth snps=anx0c/snth end if @@ -4909,15 +5390,15 @@ c csphiw=cos(phiwrad) snphiw=sin(phiwrad) c - dr=1.0d0 + dr=1.0_wp_ if(nrayr.gt.1) dr=rwmax/dble(nrayr-1) - da=2.0d0*pi/dble(nrayth) - z0t=0.0d0 + da=2.0_wp_*pi/dble(nrayth) + z0t=0.0_wp_ c do j=1,nrayr u=dble(j-1) - dffiu(j)=0.0d0 - ddffiu(j)=0.0d0 + dffiu(j)=0.0_wp_ + ddffiu(j)=0.0_wp_ do k=1,nrayth alfak=(k-1)*da dcsiw=dr*cos(alfak)*wcsi @@ -4943,7 +5424,7 @@ c ppcsi=u*dr*cos(alfak)*rcicsi ppeta=u*dr*sin(alfak)*rcieta c - anzt=1.0d0/sqrt(1.0d0+ppcsi**2+ppeta**2) + anzt=1.0_wp_/sqrt(1.0_wp_+ppcsi**2+ppeta**2) ancsi=ppcsi*anzt aneta=ppeta*anzt c @@ -4954,7 +5435,7 @@ c any=-anxt*snps+csps*(anyt*csth+anzt*snth) anz= anzt*csth-anyt*snth c - an20=1.0d0 + an20=1.0_wp_ an0=sqrt(an20) anx0=anx any0=any @@ -4974,9 +5455,9 @@ c ypwrk0(1,j,k) = anx0/an0 ypwrk0(2,j,k) = any0/an0 ypwrk0(3,j,k) = anz0/an0 - ypwrk0(4,j,k) = 0.0d0 - ypwrk0(5,j,k) = 0.0d0 - ypwrk0(6,j,k) = 0.0d0 + ypwrk0(4,j,k) = 0.0_wp_ + ypwrk0(5,j,k) = 0.0_wp_ + ypwrk0(6,j,k) = 0.0_wp_ c ytmp=ywrk0(:,j,k) yptmp=ypwrk0(:,j,k) @@ -4985,57 +5466,57 @@ c if(ipol.eq.0) then call pol_limit(ext(j,k,0),eyt(j,k,0)) qq=abs(ext(j,k,0))**2-abs(eyt(j,k,0))**2 - uu=2.0d0*dble(ext(j,k,0)*dconjg(eyt(j,k,0))) - vv=2.0d0*dimag(ext(j,k,0)*dconjg(eyt(j,k,0))) + uu=2.0_wp_*dble(ext(j,k,0)*dconjg(eyt(j,k,0))) + vv=2.0_wp_*dimag(ext(j,k,0)*dconjg(eyt(j,k,0))) call polellipse(qq,uu,vv,psipol0,chipol0) else - qq=cos(2.0d0*chipol0*cvdr)*cos(2.0d0*psipol0*cvdr) - uu=cos(2.0d0*chipol0*cvdr)*sin(2.0d0*psipol0*cvdr) - vv=sin(2.0d0*chipol0*cvdr) - if(qq**2.lt.1.0d0) then + qq=cos(2.0_wp_*chipol0*cvdr)*cos(2.0_wp_*psipol0*cvdr) + uu=cos(2.0_wp_*chipol0*cvdr)*sin(2.0_wp_*psipol0*cvdr) + vv=sin(2.0_wp_*chipol0*cvdr) + if(qq**2.lt.1.0_wp_) then c deltapol=phix-phiy, phix =0 deltapol=atan2(vv,uu) - ext(j,k,0)= sqrt((1.0d0+qq)/2) - eyt(j,k,0)= sqrt((1.0d0-qq)/2)*exp(-ui*deltapol) + ext(j,k,0)= sqrt((1.0_wp_+qq)/2) + eyt(j,k,0)= sqrt((1.0_wp_-qq)/2)*exp(-ui*deltapol) else - if(qq.gt.0.0d0) then - ext(j,k,0)= 1.0d0 - eyt(j,k,0)= 0.0d0 + if(qq.gt.0.0_wp_) then + ext(j,k,0)= 1.0_wp_ + eyt(j,k,0)= 0.0_wp_ else - eyt(j,k,0)= 1.0d0 - ext(j,k,0)= 0.0d0 + eyt(j,k,0)= 1.0_wp_ + ext(j,k,0)= 0.0_wp_ end if end if endif c do iv=1,3 - gri(iv,j,k)=0.0d0 - dgrad2v(iv,j,k)=0.0d0 - du10(iv,j,k)=0.0d0 + gri(iv,j,k)=0.0_wp_ + dgrad2v(iv,j,k)=0.0_wp_ + du10(iv,j,k)=0.0_wp_ do jv=1,3 - ggri(iv,jv,j,k)=0.0d0 + ggri(iv,jv,j,k)=0.0_wp_ end do end do - grad2(j,k)=0.0d0 + grad2(j,k)=0.0_wp_ c dd=anx0**2+any0**2+anz0**2-an20 - vgradi=0.0d0 - ddi=2.0d0*vgradi + vgradi=0.0_wp_ + ddi=2.0_wp_*vgradi c r0=sqrt(x0**2+y0**2) - x0m=x0/1.0d2 - y0m=y0/1.0d2 - r0m=r0/1.0d2 - z0m=z0/1.0d2 + x0m=x0/1.0e2_wp_ + y0m=y0/1.0e2_wp_ + r0m=r0/1.0e2_wp_ + z0m=z0/1.0e2_wp_ if(j.eq.nrayr) then write(33,111) izero,j,k,zero,x0m,y0m,r0m,z0m, . psinv,zero,anpl,zero,one end if if(j.eq.1.and.k.eq.1) then write(17,99) zero,zero,zero,zero - write(4,99) zero,r0m,z0m,atan2(y0m,x0m)*180.0d0/pi, + write(4,99) zero,r0m,z0m,atan2(y0m,x0m)*180.0_wp_/pi, . psinv,one,dens,tekev,brr,bphi,bzz, - . ajphi*1.0d-6,sqrt(anpl**2+anpr**2),anpl,zero, + . ajphi*1.0e-6_wp_,sqrt(anpl**2+anpr**2),anpl,zero, . zero,zero,zero,zero,zero,zero,zero,one end if end do @@ -5053,23 +5534,32 @@ c 99 format(24(1x,e16.8e3)) 111 format(3i5,20(1x,e16.8e3)) end - - - +c +c +c subroutine ic_rt2 - use const_and_precisions, only : izero,zero,one,pi,cvdr=>degree + use const_and_precisions, only : wp_,izero,zero,one,pi, + . cvdr=>degree use graydata_par, only : psipol0,chipol0 - - implicit real*8 (a-h,o-z) - parameter(ndim=6,ndimm=3) - parameter(jmx=31,kmx=36) - dimension ywrk0(ndim,jmx,kmx),ypwrk0(ndim,jmx,kmx) - dimension ytmp(ndim),yptmp(ndim) - dimension yyrfl(jmx,kmx,ndim) - dimension xc0(ndimm,jmx,kmx),du10(ndimm,jmx,kmx) - dimension grad2(jmx,kmx),dgrad2v(ndimm,jmx,kmx) - dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx) - complex*16 ext(jmx,kmx,0:4),eyt(jmx,kmx,0:4) + implicit none +c local constants + integer, parameter :: ndim=6,ndimm=3,jmx=31,kmx=36 +c local variables + integer :: j,k,iv,jv,iproj,nfilp + real(wp_) :: x0,y0,z0,an20,an0,anx0,any0,anz0,vgradi, + . r0,x0m,y0m,r0m,z0m,strfl11,qq,uu,vv + real(wp_), dimension(ndim) :: ytmp,yptmp +c common/external functions/variables + integer :: nrayr,nrayth + real(wp_) :: dd,an2s,an2,fdia,bdotgr,ddi,ddr11,psinv, + . dens,ddens,tekev,anpl,anpr,brr,bphi,bzz,ajphi + real(wp_), dimension(jmx,kmx) :: grad2 + real(wp_), dimension(3,jmx,kmx) :: gri + real(wp_), dimension(ndim,jmx,kmx) :: ywrk0,ypwrk0 + real(wp_), dimension(jmx,kmx,ndim) :: yyrfl + real(wp_), dimension(ndimm,jmx,kmx) :: xc0,du10,dgrad2v + real(wp_), dimension(3,3,jmx,kmx) :: ggri + complex(wp_), dimension(jmx,kmx,0:4) :: ext,eyt c common/nray/nrayr,nrayth common/wrk/ywrk0,ypwrk0 @@ -5086,7 +5576,7 @@ c common/parpl/brr,bphi,bzz,ajphi common/dens/dens,ddens common/tete/tekev - +c do j=1,nrayr do k=1,nrayth x0=yyrfl(j,k,1) @@ -5112,9 +5602,9 @@ c ypwrk0(1,j,k) = anx0/an0 ypwrk0(2,j,k) = any0/an0 ypwrk0(3,j,k) = anz0/an0 - ypwrk0(4,j,k) = 0.0d0 - ypwrk0(5,j,k) = 0.0d0 - ypwrk0(6,j,k) = 0.0d0 + ypwrk0(4,j,k) = 0.0_wp_ + ypwrk0(5,j,k) = 0.0_wp_ + ypwrk0(6,j,k) = 0.0_wp_ c ytmp=ywrk0(:,j,k) yptmp=ypwrk0(:,j,k) @@ -5122,38 +5612,38 @@ c call pol_limit(ext(j,k,0),eyt(j,k,0)) qq=abs(ext(j,k,0))**2-abs(eyt(j,k,0))**2 - uu=2.0d0*dble(ext(j,k,0)*dconjg(eyt(j,k,0))) - vv=2.0d0*dimag(ext(j,k,0)*dconjg(eyt(j,k,0))) + uu=2.0_wp_*dble(ext(j,k,0)*dconjg(eyt(j,k,0))) + vv=2.0_wp_*dimag(ext(j,k,0)*dconjg(eyt(j,k,0))) call polellipse(qq,uu,vv,psipol0,chipol0) c do iv=1,3 - gri(iv,j,k)=0.0d0 - dgrad2v(iv,j,k)=0.0d0 - du10(iv,j,k)=0.0d0 + gri(iv,j,k)=0.0_wp_ + dgrad2v(iv,j,k)=0.0_wp_ + du10(iv,j,k)=0.0_wp_ do jv=1,3 - ggri(iv,jv,j,k)=0.0d0 + ggri(iv,jv,j,k)=0.0_wp_ end do end do - grad2(j,k)=0.0d0 + grad2(j,k)=0.0_wp_ c dd=anx0**2+any0**2+anz0**2-an20 - vgradi=0.0d0 - ddi=2.0d0*vgradi + vgradi=0.0_wp_ + ddi=2.0_wp_*vgradi c r0=sqrt(x0**2+y0**2) - x0m=x0/1.0d2 - y0m=y0/1.0d2 - r0m=r0/1.0d2 - z0m=z0/1.0d2 + x0m=x0/1.0e2_wp_ + y0m=y0/1.0e2_wp_ + r0m=r0/1.0e2_wp_ + z0m=z0/1.0e2_wp_ if(j.eq.nrayr) then write(33,111) izero,j,k,zero,x0m,y0m,r0m,z0m, . psinv,zero,anpl,zero,one end if if(j.eq.1.and.k.eq.1) then write(17,99) zero,zero,zero,zero - write(4,99) strfl11,r0m,z0m,atan2(y0m,x0m)*180.0d0/pi, + write(4,99) strfl11,r0m,z0m,atan2(y0m,x0m)*180.0_wp_/pi, . psinv,one,dens,tekev,brr,bphi,bzz, - . ajphi*1.0d-6,sqrt(anpl**2+anpr**2),anpl,zero, + . ajphi*1.0e-6_wp_,sqrt(anpl**2+anpr**2),anpl,zero, . zero,zero,zero,zero,zero,zero,zero,one end if end do @@ -5172,26 +5662,35 @@ c 111 format(3i5,20(1x,e16.8e3)) end c -c ray power weigth coefficient q(j) +c c subroutine pweigth +c ray power weigth coefficient q(j) +c + use const_and_precisions, only : wp_ use graydata_par, only : rwmax - - implicit real*8(a-h,o-z) - parameter(jmx=31) - dimension q(jmx) + implicit none +c local constants + integer, parameter :: jmx=31 +c local variables + integer :: j,k + real(wp_) :: dr,r1,r2,summ,sm +c common/external functions/variables + integer :: nrayr,nrayth + real(wp_), dimension(jmx) :: q +c common/qw/q common/nray/nrayr,nrayth c - dr=1.0d0 + dr=1.0_wp_ if(nrayr.gt.1) dr=rwmax/dble(nrayr-1) - r1=0.0d0 - summ=0.0d0 - q(1)=1.0d0 + r1=0.0_wp_ + summ=0.0_wp_ + q(1)=1.0_wp_ if(nrayr.gt.1) then do j=1,nrayr - r2=(dble(j)-0.5d0)*dr - q(j)=(exp(-2.0d0*r1**2)-exp(-2.0d0*r2**2)) + r2=(dble(j)-0.5_wp_)*dr + q(j)=(exp(-2.0_wp_*r1**2)-exp(-2.0_wp_*r2**2)) r1=r2 summ=summ+q(j) end do @@ -5215,11 +5714,17 @@ c c subroutine valpsispl(rpsi,voli,dervoli,areai,rrii,rbavi, . bmxi,bmni,fci,intp) + use const_and_precisions, only : wp_ use simplespline, only :spli,splid use magsurf_data, only : rpstab,cvol,crri,crbav,cbmx,cbmn, . carea,cfc,npsi - - implicit real*8 (a-h,o-z) + implicit none +c arguments + integer :: intp + real(wp_) :: rpsi,voli,dervoli,areai,rrii,rbavi,bmxi,bmni,fci +c local variables + integer :: ip + real(wp_) :: dps c ip=int((npsi-1)*rpsi+1) c if(ip.eq.0) ip=1 @@ -5246,11 +5751,16 @@ c c c subroutine ratioj(rpsi,ratjai,ratjbi,ratjpli) + use const_and_precisions, only : wp_ use magsurf_data, only : rpstab,cratja,cratjb,cratjpl,npsi use simplespline, only :spli - - implicit real*8 (a-h,o-z) - + implicit none +c arguments + real(wp_) :: rpsi,ratjai,ratjbi,ratjpli +c local variables + integer :: ip + real(wp_) :: dps +c ip=int((npsi-1)*rpsi+1) c if(ip.eq.0) ip=1 c if(ip.eq.npsi) ip=npsi-1 @@ -5259,24 +5769,33 @@ c if(ip.eq.npsi) ip=npsi-1 ratjai=spli(cratja,npsi,ip,dps) ratjbi=spli(cratjb,npsi,ip,dps) ratjpli=spli(cratjpl,npsi,ip,dps) +c return end c c c subroutine pabs_curr(i,j,k) - use const_and_precisions, only : pi + use const_and_precisions, only : wp_,pi use graydata_flags, only : iequil,dst use graydata_par, only : sgnbphi use graydata_anequil, only : b0,rr0m,rpam - - implicit real*8 (a-h,o-z) - parameter(jmx=31,kmx=36,nmx=8000) - dimension psjki(jmx,kmx,nmx) - dimension tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx) - dimension pdjki(jmx,kmx,nmx),ppabs(jmx,kmx,nmx) - dimension currj(jmx,kmx,nmx),didst(jmx,kmx,nmx),ccci(jmx,kmx,nmx) - dimension q(jmx),tau1v(jmx,kmx) + implicit none +c local constants + integer, parameter :: jmx=31,kmx=36,nmx=8000 +c arguments + integer :: i,j,k +c local variables + integer :: intp + real(wp_) :: dvvl,rbavi,rrii,rhop0,rhop,voli,dervoli,area,adnm, + . dtau,effjcdav,dcidst,dpdst +c common/external functions/variables + real(wp_) :: p0mw,dersdst,alpha,effjcd,akim,tau0,psinv, + . bmxi,bmni,fci + real(wp_), dimension(jmx) :: q + real(wp_), dimension(jmx,kmx) :: tau1v + real(wp_), dimension(jmx,kmx,nmx) :: psjki,tauv,alphav,pdjki, + . ppabs,currj,didst,ccci c common/psjki/psjki common/atjki/tauv,alphav @@ -5284,20 +5803,16 @@ c common/pcjki/ppabs,ccci common/dcjki/didst common/tau1v/tau1v -c common/qw/q common/p0/p0mw -c common/dersdst/dersdst -c common/absor/alpha,effjcd,akim,tau0 -c common/psival/psinv common/bmxmn/bmxi,bmni common/fc/fci c - dvvl=1.0d0 - rbavi=1.0d0 + dvvl=1.0_wp_ + rbavi=1.0_wp_ rrii=rr0m rhop0=sqrt(psjki(j,k,i-1)) intp=1 @@ -5306,9 +5821,9 @@ c call valpsispl(rhop,voli,dervoli,area,rrii, . rbavi,bmxi,bmni,fci,intp) dvvl=dervoli*abs(rhop-rhop0) - if(dvvl.le.0.0d0) dvvl=1.0d0 + if(dvvl.le.0.0_wp_) dvvl=1.0_wp_ c - adnm=2.0d0*pi*rrii + adnm=2.0_wp_*pi*rrii c c calcolo della corrente cohen: currtot [MA] c calcolo della densita' corrente cohen: currj [MA/m^2] @@ -5320,17 +5835,17 @@ c call ecrh_cd c alphav(j,k,i)=alpha - dtau=(alphav(j,k,i)+alphav(j,k,i-1))*dersdst*dst/2.0d0 + dtau=(alphav(j,k,i)+alphav(j,k,i-1))*dersdst*dst/2.0_wp_ tauv(j,k,i)=tauv(j,k,i-1)+dtau dpdst=p0mw*q(j)*exp(-tauv(j,k,i)-tau1v(j,k))* . alphav(j,k,i)*dersdst pdjki(j,k,i)=dpdst*dst/dvvl ppabs(j,k,i)=p0mw*q(j)*exp(-tau1v(j,k))* - . (1.0d0-exp(-tauv(j,k,i))) + . (1.0_wp_-exp(-tauv(j,k,i))) effjcdav=rbavi*effjcd currj(j,k,i)=sgnbphi*effjcdav*pdjki(j,k,i) didst(j,k,i)=sgnbphi*effjcdav*dpdst/adnm - dcidst=(didst(j,k,i)+didst(j,k,i-1))/2.0d0 + dcidst=(didst(j,k,i)+didst(j,k,i-1))/2.0_wp_ ccci(j,k,i)=ccci(j,k,i-1)+dcidst*dst return end @@ -5338,31 +5853,35 @@ c c c subroutine ecrh_cd - use const_and_precisions, only : mc2=>mc2_ + use const_and_precisions, only : wp_,mc2=>mc2_ use dispersion, only : harmnumber, warmdisp use graydata_flags, only : iwarm,ilarm,ieccd,imx - - implicit real*8 (a-h,o-z) - parameter(taucr=12.0d0,xxcr=16.0d0,eps=1.d-8) - complex*16 ex,ey,ez + implicit none +c local constants + real(wp_), parameter :: taucr=12.0_wp_,xxcr=16.0_wp_,eps=1.e-8_wp_ +c local variables + integer :: lrm,ithn + real(wp_) :: amu,ratiovgr,rbn,rbx,zeff +c common/external functions/variables + integer :: nhmin,nhmax,iokhawa,ierr + real(wp_) :: xg,yg,anpl,anpr,vgm,derdnm,ak0,akinv,fhz,sox, + . anprre,anprim,alpha,effjcd,akim,tau,psinv,tekev,dens,ddens, + . btot,bmax,bmin,fc + complex(wp_) :: ex,ey,ez + real(wp_) :: fzeff,temp c common/nharm/nhmin,nhmax -c common/xgxg/xg common/ygyg/yg common/nplr/anpl,anpr common/vgv/vgm,derdnm common/parwv/ak0,akinv,fhz common/mode/sox -c common/nprw/anprre,anprim common/epolar/ex,ey,ez -c common/absor/alpha,effjcd,akim,tau -c common/ierr/ierr common/iokh/iokhawa -c common/psival/psinv common/tete/tekev common/dens/dens,ddens @@ -5372,13 +5891,13 @@ c c c absorption computation c - alpha=0.0d0 - akim=0.0d0 - effjcd=0.0d0 + alpha=0.0_wp_ + akim=0.0_wp_ + effjcd=0.0_wp_ c tekev=temp(psinv) c - if(tekev.le.0.0d0.or.tau.gt.taucr) return + if(tekev.le.0.0_wp_.or.tau.gt.taucr) return c amu=mc2/tekev call harmnumber(yg,amu,anpl,nhmin,nhmax,iwarm) @@ -5389,10 +5908,10 @@ c * iwarm,imx,ex,ey,ez) c akim=ak0*anprim - ratiovgr=2.0d0*anpr/derdnm -c ratiovgr=2.0d0*anpr/derdnm*vgm - alpha=2.0d0*akim*ratiovgr - if(alpha.lt.0.0d0) then + ratiovgr=2.0_wp_*anpr/derdnm +c ratiovgr=2.0_wp_*anpr/derdnm*vgm + alpha=2.0_wp_*akim*ratiovgr + if(alpha.lt.0.0_wp_) then ierr=94 print*,' IERR = ', ierr,' alpha negative' end if @@ -5413,52 +5932,53 @@ c c c subroutine eccd(yg,anpl,anprre,amu,zeff,rbn,rbx,fc,ex,ey,ez, - * dens,psinv,ieccd,nhmn,nhmx,ithn,effjcd,iokhawa,ierr) + . dens,psinv,ieccd,nhmn,nhmx,ithn,effjcd,iokhawa,ierr) + use const_and_precisions, only : wp_,pi,qesi=>e_,mesi=>me_, + . vcsi=>c_,qe=>ecgs_,me=>mecgs_,vc=>ccgs_,mc2=>mc2_ use green_func_p, only: SpitzFuncCoeff - use const_and_precisions, only : pi,qesi=>e_,mesi=>me_,vcsi=>c_, - * qe=>ecgs_,me=>mecgs_,vc=>ccgs_,mc2=>mc2_ use conical, only : fconic implicit none - real*8 anum,denom,resp,resj,effjcd,coullog,dens - real*8 yg,anpl,anpr,amu,anprre,anprim - real*8 mc2m2,anucc,canucc,ddens - real*8 ceff,Zeff,psinv - real*8 rbn,rbx,alams,fp0s,pa,fc - real*8 fjch,fjncl,fjch0 - real*8 cst2,eccdpar(5) - complex*16 ex,ey,ez - integer ieccd,nhmn,nhmx,nhn,ithn,ierr,iokhawa - external fjch,fjncl,fjch0 - - parameter(mc2m2=1.0d0/mc2**2) - parameter(canucc=2.0d13*pi*qe**4/(me**2*vc**3)) - parameter(ceff=qesi/(mesi*vcsi)) - - anum=0.0d0 - denom=0.0d0 - effjcd=0.0d0 - - coullog=48.0d0-log(1.0d7*dens*mc2m2*amu**2) +c local constants + real(wp_), parameter :: mc2m2=1.0_wp_/mc2**2, + . canucc=2.0e13_wp_*pi*qe**4/(me**2*vc**3),ceff=qesi/(mesi*vcsi) +c arguments + integer ieccd,nhmn,nhmx,ithn,iokhawa,ierr + real(wp_) :: yg,anpl,anprre,amu,Zeff,rbn,rbx, + . fc,dens,psinv,effjcd + complex(wp_) :: ex,ey,ez +c local variables + integer nhn + real(wp_) :: anum,denom,resp,resj,coullog,anucc,alams, + . fp0s,pa,cst2 + real(wp_), dimension(5) :: eccdpar +c + real(wp_), external :: fjch,fjncl,fjch0 +c + anum=0.0_wp_ + denom=0.0_wp_ + effjcd=0.0_wp_ +c + coullog=48.0_wp_-log(1.0e7_wp_*dens*mc2m2*amu**2) anucc=canucc*dens*coullog - +c c nhmx=nhmn - +c eccdpar(1)=zeff - +c select case(ieccd) - +c case(1) c cohen model c rbn=B/B_min c rbx=B/B_max -c cst2=1.0d0-B/B_max +c cst2=1.0_wp_-B/B_max c alams=sqrt(1-B_min/B_max) c Zeff < 31 !!! c fp0s= P_a (alams) - cst2=1.0d0-rbx - if(cst2.lt.1d-6) cst2=0.0d0 - alams=sqrt(1.0d0-rbx/rbn) - pa=sqrt(32.0d0/(Zeff+1.0d0)-1.0d0)/2.0d0 + cst2=1.0_wp_-rbx + if(cst2.lt.1.0e-6_wp_) cst2=0.0_wp_ + alams=sqrt(1.0_wp_-rbx/rbn) + pa=sqrt(32.0_wp_/(Zeff+1.0_wp_)-1.0_wp_)/2.0_wp_ fp0s=fconic(alams,pa,0) eccdpar(2)=rbn eccdpar(3)=alams @@ -5466,16 +5986,16 @@ c fp0s= P_a (alams) eccdpar(5)=fp0s do nhn=nhmn,nhmx call curr_int(yg,anpl,anprre,amu,ex,ey,ez,nhn,ithn,cst2, - * fjch,eccdpar,5,resj,resp,iokhawa,ierr) + . fjch,eccdpar,5,resj,resp,iokhawa,ierr) anum=anum+resj denom=denom+resp end do case(2:9) - cst2=0.0d0 + cst2=0.0_wp_ do nhn=nhmn,nhmx call curr_int(yg,anpl,anprre,amu,ex,ey,ez,nhn,ithn,cst2, - * fjch0,eccdpar,1,resj,resp,iokhawa,ierr) + . fjch0,eccdpar,1,resj,resp,iokhawa,ierr) anum=anum+resj denom=denom+resp end do @@ -5484,15 +6004,15 @@ c fp0s= P_a (alams) c neoclassical model: c ft=1-fc trapped particle fraction c rzfc=(1+Zeff)/fc - cst2=1.0d0-rbx - if(cst2.lt.1d-6) cst2=0.0d0 + cst2=1.0_wp_-rbx + if(cst2.lt.1.0e-6_wp_) cst2=0.0_wp_ call SpitzFuncCoeff(amu,Zeff,fc) eccdpar(2) = fc eccdpar(3) = rbx eccdpar(4) = psinv do nhn=nhmn,nhmx call curr_int(yg,anpl,anprre,amu,ex,ey,ez,nhn,ithn,cst2, - * fjncl,eccdpar,4,resj,resp,iokhawa,ierr) + . fjncl,eccdpar,4,resj,resp,iokhawa,ierr) anum=anum+resj denom=denom+resp end do @@ -5504,7 +6024,7 @@ c rzfc=(1+Zeff)/fc c c effjpl = / /(B_min/) [A m /W] c - if(denom.gt.0.0d0) effjcd=-ceff*anum/(anucc*denom) + if(denom.gt.0.0_wp_) effjcd=-ceff*anum/(anucc*denom) return end c @@ -5512,46 +6032,58 @@ c c subroutine curr_int(yg,anpl,anprre,amu,ex,ey,ez,nhn,ithn,cst2, * fcur,eccdpar,necp,resj,resp,iokhawa,ierr) + use const_and_precisions, only : wp_ use quadpack, only : dqagsmv - implicit real*8(a-h,o-z) - parameter(epsa=0.0d0,epsr=1.0d-2) - parameter(xxcr=16.0d0) - parameter (lw=5000,liw=lw/4) - parameter (nfpp=15) - complex*16 ex,ey,ez - dimension w(lw),iw(liw),eccdpar(necp),apar(nfpp+necp) - external fcur,fpp - + implicit none +c local constants + real(wp_), parameter :: epsa=0.0_wp_,epsr=1.0e-2_wp_,xxcr=16.0_wp_ + integer, parameter :: lw=5000,liw=lw/4,nfpp=15 +c arguments + integer :: i,nhn,ithn,necp,iokhawa,ierr + real(wp_) :: yg,anpl,anprre,amu,cst2,resj,resp + real(wp_), dimension(necp) :: eccdpar + complex(wp_) :: ex,ey,ez +c local variables + integer :: neval,ier,last,npar + integer, dimension(liw) :: iw + real(wp_) :: anpl2,dnl,ygn,ygn2,resj1,resj2,rdu2,upltp,upltm, + . rdu,rdut,rdu2t,duu,uu1,uu2,xx1,xx2,ej,ej1,ej2,epp,uplp,uplm, + . cst + real(wp_), dimension(lw) :: w + real(wp_), dimension(nfpp+necp) :: apar +c + real(wp_), external :: fcur,fpp +c c EC power and current densities - +c anpl2=anpl*anpl - dnl=1.0d0-anpl2 + dnl=1.0_wp_-anpl2 ygn=nhn*yg ygn2=ygn*ygn - resj=0.0d0 - resj1=0.0d0 - resj2=0.0d0 - resp=0.0d0 + resj=0.0_wp_ + resj1=0.0_wp_ + resj2=0.0_wp_ + resp=0.0_wp_ c - rdu2=anpl2+ygn2-1.0d0 - uplp=0.0d0 - uplm=0.0d0 - upltp=0.0d0 - upltm=0.0d0 + rdu2=anpl2+ygn2-1.0_wp_ + uplp=0.0_wp_ + uplm=0.0_wp_ + upltp=0.0_wp_ + upltm=0.0_wp_ c - if (rdu2.ge.0.0d0) then + if (rdu2.ge.0.0_wp_) then rdu=sqrt(rdu2) uplp=(anpl*ygn+rdu)/dnl uplm=(anpl*ygn-rdu)/dnl c uu1=uplm uu2=uplp - xx1=amu*(anpl*uu1+ygn-1.0d0) - xx2=amu*(anpl*uu2+ygn-1.0d0) + xx1=amu*(anpl*uu1+ygn-1.0_wp_) + xx2=amu*(anpl*uu2+ygn-1.0_wp_) c - if(xx2.gt.xxcr) uu2=(xxcr/amu-ygn+1.0d0)/anpl - if(xx1.gt.xxcr) uu1=(xxcr/amu-ygn+1.0d0)/anpl + if(xx2.gt.xxcr) uu2=(xxcr/amu-ygn+1.0_wp_)/anpl + if(xx1.gt.xxcr) uu1=(xxcr/amu-ygn+1.0_wp_)/anpl duu=abs(uu1-uu2) c apar(1) = yg @@ -5575,20 +6107,20 @@ c apar(nfpp+i) = eccdpar(i) end do c - if(duu.gt.1.d-6) then + if(duu.gt.1.0e-6_wp_) then call dqagsmv(fpp,uu1,uu2,apar,npar,epsa,epsr,resp,epp,neval, . ier,liw,lw,last,iw,w) if (ier.gt.0) ierr=90 end if c - rdu2t=cst2*anpl2+ygn2-1.0d0 + rdu2t=cst2*anpl2+ygn2-1.0_wp_ c - if (rdu2t.lt.0.0d0.or.cst2.eq.0.0d0) then + if (rdu2t.lt.0.0_wp_.or.cst2.eq.0.0_wp_) then c c resonance curve does not cross the trapping region c iokhawa=0 - if(duu.gt.1.d-4) then + if(duu.gt.1.0e-4_wp_) then call dqagsmv(fcur,uu1,uu2,apar,npar,epsa,epsr, . resj,ej,neval,ier,liw,lw,last,iw,w) if (ier.gt.0) ierr=91 @@ -5600,23 +6132,23 @@ c iokhawa=1 rdut=sqrt(rdu2t) cst=sqrt(cst2) - upltm=(cst2*anpl*ygn-cst*rdut)/(1.0d0-cst2*anpl2) - upltp=(cst2*anpl*ygn+cst*rdut)/(1.0d0-cst2*anpl2) + upltm=(cst2*anpl*ygn-cst*rdut)/(1.0_wp_-cst2*anpl2) + upltp=(cst2*anpl*ygn+cst*rdut)/(1.0_wp_-cst2*anpl2) c uu1=uplm uu2=upltm - xx1=amu*(anpl*uu1+ygn-1.0d0) - xx2=amu*(anpl*uu2+ygn-1.0d0) + xx1=amu*(anpl*uu1+ygn-1.0_wp_) + xx2=amu*(anpl*uu2+ygn-1.0_wp_) if(xx1.lt.xxcr.or.xx2.lt.xxcr) then - if(xx2.gt.xxcr) uu2=(xxcr/amu-ygn+1.0d0)/anpl - if(xx1.gt.xxcr) uu1=(xxcr/amu-ygn+1.0d0)/anpl + if(xx2.gt.xxcr) uu2=(xxcr/amu-ygn+1.0_wp_)/anpl + if(xx1.gt.xxcr) uu1=(xxcr/amu-ygn+1.0_wp_)/anpl duu=abs(uu1-uu2) - if(duu.gt.1.d-6) then + if(duu.gt.1.0e-6_wp_) then call dqagsmv(fcur,uu1,uu2,apar,npar,epsa,epsr, . resj1,ej1,neval,ier,liw,lw,last,iw,w) if (ier.gt.0) then - if (abs(resj1).lt.1.0d-10) then - resj1=0.0d0 + if (abs(resj1).lt.1.0e-10_wp_) then + resj1=0.0_wp_ else ierr=92 end if @@ -5626,13 +6158,13 @@ c c uu1=upltp uu2=uplp - xx1=amu*(anpl*uu1+ygn-1.0d0) - xx2=amu*(anpl*uu2+ygn-1.0d0) + xx1=amu*(anpl*uu1+ygn-1.0_wp_) + xx2=amu*(anpl*uu2+ygn-1.0_wp_) if(xx1.lt.xxcr.or.xx2.lt.xxcr) then - if(xx2.gt.xxcr) uu2=(xxcr/amu-ygn+1.0d0)/anpl - if(xx1.gt.xxcr) uu1=(xxcr/amu-ygn+1.0d0)/anpl + if(xx2.gt.xxcr) uu2=(xxcr/amu-ygn+1.0_wp_)/anpl + if(xx1.gt.xxcr) uu1=(xxcr/amu-ygn+1.0_wp_)/anpl duu=abs(uu1-uu2) - if(duu.gt.1.d-6) then + if(duu.gt.1.0e-6_wp_) then call dqagsmv(fcur,uu1,uu2,apar,npar,epsa,epsr, . resj2,ej2,neval,ier,liw,lw,last,iw,w) if (ier.gt.0) then @@ -5654,7 +6186,6 @@ c ith=1 : polarization term Larmor radius expansion to lowest order c ith=2 : full polarization term (J Bessel) c function fpp(upl,extrapar,npar) -c c integration variable upl passed explicitly. other variables passed c as array of extra parameters of length npar=size(extrapar) c @@ -5674,12 +6205,19 @@ c extrapar(13) = uplp c extrapar(14) = uplm c extrapar(15) = ygn c - use const_and_precisions, only : ui=>im + use const_and_precisions, only : wp_,ui=>im use math, only : fact - implicit real*8 (a-h,o-z) - complex*16 ex,ey,ez,emxy,epxy - dimension extrapar(npar) - + implicit none +c arguments + integer :: npar + real(wp_) :: upl,fpp + real(wp_), dimension(npar) :: extrapar +c local variables + integer :: ithn,nhn,nm,np + real(wp_) :: yg,anpl,amu,anprre,uplp,uplm,ygn,upr,upr2,gam,ee, + . thn2,thn2u,bb,cth,ajbnm,ajbnp,ajbn + complex(wp_) :: ex,ey,ez,emxy,epxy +c yg=extrapar(1) anpl=extrapar(2) amu=extrapar(3) @@ -5693,23 +6231,23 @@ c uplm=extrapar(14) ygn=extrapar(15) - upr2=(1.0d0-anpl**2)*(uplp-upl)*(upl-uplm) + upr2=(1.0_wp_-anpl**2)*(uplp-upl)*(upl-uplm) gam=anpl*upl+ygn ee=exp(-amu*(gam-1)) - thn2=1.0d0 + thn2=1.0_wp_ thn2u=upr2*thn2 if(ithn.gt.0) then emxy=ex-ui*ey epxy=ex+ui*ey - if(upr2.gt.0.0d0) then + if(upr2.gt.0.0_wp_) then upr=sqrt(upr2) bb=anprre*upr/yg if(ithn.eq.1) then c Larmor radius expansion polarization term at lowest order - cth=1.0d0 - if(nhn.gt.1) cth=(0.5d0*bb)**(nhn-1)*nhn/fact(nhn) - thn2=(0.5d0*cth*abs(emxy+ez*anprre*upl/ygn))**2 + cth=1.0_wp_ + if(nhn.gt.1) cth=(0.5_wp_*bb)**(nhn-1)*nhn/fact(nhn) + thn2=(0.5_wp_*cth*abs(emxy+ez*anprre*upl/ygn))**2 thn2u=upr2*thn2 else c Full polarization term @@ -5718,15 +6256,17 @@ c Full polarization term ajbnm=dbesjn(nm, bb) ajbnp=dbesjn(np, bb) ajbn=dbesjn(nhn, bb) - thn2u=(abs(ez*ajbn*upl+upr*(ajbnp*epxy+ajbnm*emxy)/2.0d0))**2 + thn2u=(abs(ez*ajbn*upl+upr*(ajbnp*epxy+ajbnm*emxy) + . /2.0_wp_))**2 end if end if end if fpp=ee*thn2u +c return end - +c c computation of integral for current density c fjch integrand for Cohen model with trapping c fjch0 integrand for Cohen model without trapping @@ -5734,7 +6274,6 @@ c fjncl integrand for momentum conserv. model K(u) from Maruschenko c gg=F(u)/u with F(u) as in Cohen paper c function fjch(upl,extrapar,npar) -c c integration variable upl passed explicitly. Other variables passed c as array of extra parameters of length npar=size(extrapar). c variables with index 1..15 must be passed to fpp @@ -5763,11 +6302,20 @@ c extrapar(18) = alams c extrapar(19) = pa c extrapar(20) = fp0s c + use const_and_precisions, only : wp_ use conical, only : fconic - implicit real*8 (a-h,o-z) - complex*16 ex,ey,ez - dimension extrapar(npar) - + implicit none +c arguments + integer :: npar + real(wp_) :: upl,fjch + real(wp_), dimension(npar) :: extrapar +c local variables + integer :: nhn + real(wp_) :: anpl,anprre,uplp,uplm,ygn,zeff,rb,alams,pa,fp0s, + . upr2,gam,u2,u,z5,xi,xib,xibi,fu2b,fu2,gu,gg,dgg,alam,fp0, + . dfp0,fh,dfhl,eta,fpp + complex(wp_) :: ex,ey,ez +c anpl=extrapar(2) anprre=extrapar(4) ex=dcmplx(extrapar(5),extrapar(6)) @@ -5783,40 +6331,40 @@ c pa=extrapar(19) fp0s=extrapar(20) - upr2=(1.0d0-anpl**2)*(uplp-upl)*(upl-uplm) + upr2=(1.0_wp_-anpl**2)*(uplp-upl)*(upl-uplm) gam=anpl*upl+ygn - u2=gam*gam-1.0d0 + u2=gam*gam-1.0_wp_ u=sqrt(u2) - z5=Zeff+5.0d0 - xi=1.0d0/z5**2 - xib=1.0d0-xi - xibi=1.0d0/xib - fu2b=1.0d0+xib*u2 - fu2=1.0d0+xi*u2 - gu=(1.0d0-1.0d0/fu2b**xibi)/sqrt(fu2) + z5=Zeff+5.0_wp_ + xi=1.0_wp_/z5**2 + xib=1.0_wp_-xi + xibi=1.0_wp_/xib + fu2b=1.0_wp_+xib*u2 + fu2=1.0_wp_+xi*u2 + gu=(1.0_wp_-1.0_wp_/fu2b**xibi)/sqrt(fu2) gg=u*gu/z5 - dgg=(gu+u2*(2.0d0/fu2b**(1.0d0+xibi)/sqrt(fu2)-xi*gu/fu2))/z5 + dgg=(gu+u2*(2.0_wp_/fu2b**(1.0_wp_+xibi)/sqrt(fu2)-xi*gu/fu2))/z5 - alam=sqrt(1.0d0-upr2/u2/rb) + alam=sqrt(1.0_wp_-upr2/u2/rb) fp0=fconic(alam,pa,0) - dfp0=-(pa*pa/2.0d0+0.125d0) - if (alam.lt.1.0d0) then - dfp0=-fconic(alam,pa,1)/sqrt(1.0d0-alam**2) + dfp0=-(pa*pa/2.0_wp_+0.125_wp_) + if (alam.lt.1.0_wp_) then + dfp0=-fconic(alam,pa,1)/sqrt(1.0_wp_-alam**2) end if - fh=alam*(1.0d0-alams*fp0/(alam*fp0s)) - dfhl=1.0d0-alams*dfp0/fp0s + fh=alam*(1.0_wp_-alams*fp0/(alam*fp0s)) + dfhl=1.0_wp_-alams*dfp0/fp0s eta=gam*fh*(gg/u+dgg)+upl*(anpl*u2-upl*gam)*gg*dfhl/(u2*u*rb*alam) - if(upl.lt.0.0d0) eta=-eta + if(upl.lt.0.0_wp_) eta=-eta fjch=eta*fpp(upl,extrapar,npar) +c return end - - - - function fjch0(upl,extrapar,npar) c +c +c + function fjch0(upl,extrapar,npar) c integration variable upl passed explicitly. Other variables passed c as array of extra parameters of length npar=size(extrapar). c variables with index 1..15 must be passed to fpp @@ -5841,10 +6389,18 @@ c extrapar(15) = ygn c c extrapar(16) = zeff c - implicit real*8 (a-h,o-z) - complex*16 ex,ey,ez - dimension extrapar(npar) - + use const_and_precisions, only : wp_ + implicit none +c arguments + real(wp_) upl,fjch0 + integer :: npar + real(wp_), dimension(npar) :: extrapar +c local variables + integer :: nhn + real(wp_) :: anpl,anprre,ygn,zeff,gam,u2,u,z5,xi,xib,xibi, + . fu2b,fu2,gu,gg,dgg,eta,fpp + complex(wp_) :: ex,ey,ez +c anpl=extrapar(2) anprre=extrapar(4) ex=dcmplx(extrapar(5),extrapar(6)) @@ -5855,26 +6411,26 @@ c zeff=extrapar(16) gam=anpl*upl+ygn - u2=gam*gam-1.0d0 + u2=gam*gam-1.0_wp_ u=sqrt(u2) - z5=Zeff+5.0d0 - xi=1.0d0/z5**2 - xib=1.0d0-xi - xibi=1.0d0/xib - fu2b=1.0d0+xib*u2 - fu2=1.0d0+xi*u2 - gu=(1.0d0-1.0d0/fu2b**xibi)/sqrt(fu2) + z5=Zeff+5.0_wp_ + xi=1.0_wp_/z5**2 + xib=1.0_wp_-xi + xibi=1.0_wp_/xib + fu2b=1.0_wp_+xib*u2 + fu2=1.0_wp_+xi*u2 + gu=(1.0_wp_-1.0_wp_/fu2b**xibi)/sqrt(fu2) gg=u*gu/z5 - dgg=(gu+u2*(2.0d0/fu2b**(1.0d0+xibi)/sqrt(fu2)-xi*gu/fu2))/z5 + dgg=(gu+u2*(2.0_wp_/fu2b**(1.0_wp_+xibi)/sqrt(fu2)-xi*gu/fu2))/z5 eta=anpl*gg+gam*upl*dgg/u fjch0=eta*fpp(upl,extrapar,npar) +c return end - - - - function fjncl(upl,extrapar,npar) c +c +c + function fjncl(upl,extrapar,npar) c integration variable upl passed explicitly. Other variables passed c as array of extra parameters of length npar=size(extrapar). c variables with index 1..15 must be passed to fpp @@ -5900,12 +6456,19 @@ c c extrapar(16) = zeff c extrapar(17) = fc c extrapar(18) = hb -c extrapar(29) = psin +c extrapar(19) = psin c + use const_and_precisions, only : wp_ use green_func_p, only: GenSpitzFunc - implicit real*8 (a-h,o-z) - dimension extrapar(npar) - + implicit none +c arguments + integer :: npar + real(wp_) :: upl,fjncl + real(wp_), dimension(npar) :: extrapar +c local variables + real(wp_) :: anpl,amu,ygn,zeff,fc,hb,psin,gam,u2,u,upr2, + . bth,uth,fk,dfk,alam,fu,dfu,eta,fpp +c anpl=extrapar(2) amu=extrapar(3) ygn=extrapar(15) @@ -5915,64 +6478,81 @@ c psin=extrapar(19) gam=anpl*upl+ygn - u2=gam*gam-1.0d0 + u2=gam*gam-1.0_wp_ u=sqrt(u2) upr2=u2-upl*upl - bth=sqrt(2.0d0/amu) + bth=sqrt(2.0_wp_/amu) uth=u/bth call GenSpitzFunc(amu,Zeff,fc,uth,u,gam,fk,dfk) - fk=fk*(4.0d0/amu**2) - dfk=dfk*(2.0d0/amu)*bth + fk=fk*(4.0_wp_/amu**2) + dfk=dfk*(2.0_wp_/amu)*bth alam=upr2/u2/hb call vlambda(alam,psin,fu,dfu) - eta=gam*fu*dfk/u-2.0d0*(anpl-gam*upl/u2)*fk*dfu*upl/u2/hb + eta=gam*fu*dfk/u-2.0_wp_*(anpl-gam*upl/u2)*fk*dfu*upl/u2/hb if(upl.lt.0) eta=-eta fjncl=eta*fpp(upl,extrapar,npar) return end - +c +c +c subroutine vlambda(alam,psi,fv,dfv) + use const_and_precisions, only : wp_ use magsurf_data, only : ch,ch01,tjp,tlm,njpt,nlmt,npsi use dierckx, only : fpbisp - - implicit real*8 (a-h,o-z) - parameter (nlam=41) - parameter(ksp=3,nlest=nlam+ksp+1) - - dimension xxs(1),yys(1),ffs(1) + implicit none +c local constants + integer, parameter :: nlam=41,ksp=3,nlest=nlam+ksp+1 +c arguments + real(wp_) alam,psi,fv,dfv +c local variables + integer :: njp,nlm,iwp,iwl,njest,lwrk,kwrk integer, dimension(:), allocatable :: iwrk - real*8, dimension(:), allocatable :: wrk - + real(wp_), dimension(1) :: xxs,yys,ffs + real(wp_), dimension(:), allocatable :: wrk +c njest=npsi+ksp+1 kwrk=npsi+nlam+njest+nlest+3 lwrk=4*(npsi+nlam)+11*(njest+nlest)+njest*npsi+nlest+54 allocate(iwrk(kwrk),wrk(lwrk)) - +c njp=njpt nlm=nlmt xxs(1)=sqrt(psi) yys(1)=alam - +c call fpbisp(tjp,njp,tlm,nlm,ch,ksp,ksp,xxs,1,yys,1,ffs, . wrk(1),wrk(5),iwrk(1),iwrk(2)) fv=ffs(1) - +c iwp=1+(njp-4)*(nlm-5) iwl=iwp+4 call fpbisp(tjp(1),njp,tlm(2),nlm-2,ch01,3,2, . xxs,1,yys,1,ffs,ch01(iwp),ch01(iwl),iwrk(1),iwrk(2)) dfv=ffs(1) - +c return end - - +c +c +c subroutine projxyzt(iproj,nfile) - implicit real*8 (a-h,o-z) - parameter(jmx=31,kmx=36) - dimension ywrk(6,jmx,kmx),ypwrk(6,jmx,kmx) + use const_and_precisions, only : wp_ + implicit none +c local constants + integer, parameter :: jmx=31,kmx=36 +c arguments + integer iproj,nfile +c local variables + integer :: jd,j,kkk,k + real(wp_) :: dir,rtimn,rtimx,dx,dy,dz,dirx,diry,dirz, + . csth1,snth1,csps1,snps1,xti,yti,zti,rti +c common/external functions/variables + integer :: nrayr,nrayth,istep + real(wp_) :: psinv11,st + real(wp_), dimension(6,jmx,kmx) :: ywrk,ypwrk c common/nray/nrayr,nrayth common/wrk/ywrk,ypwrk @@ -5980,8 +6560,8 @@ c common/istep/istep common/ss/st c - rtimn=1.d+30 - rtimx=-1.d-30 + rtimn=1.0e+30_wp_ + rtimx=-1.0e-30_wp_ c jd=1 if(iproj.eq.0) jd=nrayr-1 @@ -6001,10 +6581,10 @@ c c if(j.eq.1.and.k.eq.1) then csth1=dirz/dir - snth1=sqrt(1.0d0-csth1**2) - csps1=1.0d0 - snps1=0.0d0 - if(snth1.gt.0.0d0) then + snth1=sqrt(1.0_wp_-csth1**2) + csps1=1.0_wp_ + snps1=0.0_wp_ + if(snth1.gt.0.0_wp_) then csps1=diry/(dir*snth1) snps1=dirx/(dir*snth1) end if @@ -6014,16 +6594,16 @@ c zti=(dx*snps1+dy*csps1)*snth1+dz*csth1 rti=sqrt(xti**2+yti**2) Cc -C dirxt= (dirx*csps1-diry*snps1)/dir -C diryt=((dirx*snps1+diry*csps1)*csth1-dirz*snth1)/dir -C dirzt=((dirx*snps1+diry*csps1)*snth1+dirz*csth1)/dir +c dirxt= (dirx*csps1-diry*snps1)/dir +c diryt=((dirx*snps1+diry*csps1)*csth1-dirz*snth1)/dir +c dirzt=((dirx*snps1+diry*csps1)*snth1+dirz*csth1)/dir Cc -C if(k.eq.1) then -C xti1=xti -C yti1=yti -C zti1=zti -C rti1=rti -C end if +c if(k.eq.1) then +c xti1=xti +c yti1=yti +c zti1=zti +c rti1=rti +c end if if(.not.(iproj.eq.0.and.j.eq.1)) . write(nfile,111) istep,j,k,xti,yti,zti,rti,psinv11 @@ -6049,67 +6629,76 @@ c c c subroutine pec(pabs,currt) - use const_and_precisions, only : pi,one + use const_and_precisions, only : wp_,pi,one use numint, only : simpson use graydata_flags, only : ipec,ieccd,iequil,nnd,dst use graydata_anequil, only : rr0m,rpam use utils, only : locatex, locate, intlin - - implicit real*8(a-h,o-z) - parameter(nndmx=5001,jmx=31,kmx=36,nmx=8000) - parameter(rtbc=one) - dimension psjki(jmx,kmx,nmx),iiv(jmx,kmx) - dimension ppabs(jmx,kmx,nmx),ccci(jmx,kmx,nmx) - dimension pdjki(jmx,kmx,nmx),currj(jmx,kmx,nmx) - dimension xxi(nmx),ypt(nmx),yamp(nmx) - dimension rtab(nndmx),rhotv(nndmx) - dimension rtab1(0:nndmx) - dimension ajphiv(nndmx),dpdv(nndmx) - dimension dvolt(nndmx),darea(nndmx) - dimension ratjav(nndmx),ratjbv(nndmx),ratjplv(nndmx) - dimension ajplv(nndmx),ajcdav(nndmx),ajcdbv(nndmx) - dimension pins(nndmx),currins(nndmx),fi(nndmx) - parameter(llmx=21) - dimension isev(llmx) + implicit none +c local constants + integer, parameter :: nndmx=5001,jmx=31,kmx=36,nmx=8000,llmx=21 + real(wp_), parameter :: rtbc=one +c arguments + real(wp_) :: pabs,currt +c local variables + integer :: ll,it,nd,intp,kkk,i,j,k,ii,is,idecr,ise0,iis,iis1, + . iise0,iise,ind1,ind2,iind,indi,iif1,istmx,i1,ind,itb1 + integer, dimension(llmx) :: isev + real(wp_) :: stf,voli0,areai0,drt,rt,rt1,psit,psit1,voli1, + . dervoli,areai1,rrii,rbavi,bmxi,bmni,fci,ppa1,cci1,ppa2,cci2, + . dppa,didst,rhotpav,rhot2pav,rhotjav,rhotjava,rhot2java,spds, + . sccs,sccsa,drhot2pav,drhotpav,drhot2java,drhotjava,facpds, + . facjs,rhpp,rhpj,dpdvp,ajphip,ratjamx,ratjbmx,ratjplmx,rhopp, + . drhopp,rhotjfi,rhopfi,ajmxfi,drhotjfi,drhopfi,xps,dpdvmx, + . rhotp,drhotp,stmx,pins_02,pins_05,pins_085,xrhot,currtka,rhop, + . pinsr,ratjpli,ratjai,ratjbi,frhotor,frhopol,fdadrhot, + . fdvdrhot,h,psin + real(wp_), dimension(nmx) :: xxi,ypt,yamp + real(wp_), dimension(nndmx) :: ajphiv,dpdv,dvolt,darea,ratjav, + . ratjbv,ratjplv,ajplv,ajcdav,ajcdbv,pins,currins,fi,rtab,rhotv + real(wp_), dimension(0:nndmx) :: rtab1 +c common/external functions/variables + integer :: nrayr,nrayth,istep,index_rt + integer, dimension(jmx,kmx) :: iiv + real(wp_) :: alpha0,beta0,taumn,taumx,pabstot,currtot, + . psipol,chipol + real(wp_), dimension(jmx,kmx,nmx) :: psjki,ppabs,ccci,pdjki,currj c common/nray/nrayr,nrayth common/istep/istep common/index_rt/index_rt -c common/iiv/iiv common/psjki/psjki common/pcjki/ppabs,ccci common/dpjjki/pdjki,currj -c common/angles/alpha0,beta0 common/taumnx/taumn,taumx,pabstot,currtot -c common/polcof/psipol,chipol c stf=istep*dst nd=nnd - if(pabstot.gt.0.0d0) then + if(pabstot.gt.0.0_wp_) then do ll=1,llmx isev(ll)=0 end do intp=0 - voli0=0.0d0 - areai0=0.0d0 - rtab1(0)=0.0d0 + voli0=0.0_wp_ + areai0=0.0_wp_ + rtab1(0)=0.0_wp_ do it=1,nd - drt=1.0d0/dble(nd-1) + drt=1.0_wp_/dble(nd-1) rt=dble(it-1)*drt if(it.lt.nd) then - rt1=rt+drt/2.0d0 + rt1=rt+drt/2.0_wp_ else rt1=rt end if rtab(it)=rt rtab1(it)=rt1 - dpdv(it)=0.0d0 - ajphiv(it)=0.0d0 + dpdv(it)=0.0_wp_ + ajphiv(it)=0.0_wp_ if (ipec.eq.0) then psit=rt psit1=rt1 @@ -6137,7 +6726,7 @@ c ise0=0 ii=iiv(j,k) if (ii.lt.nmx) then - if(psjki(j,k,ii+1).ne.0.0d0) ii=ii+1 + if(psjki(j,k,ii+1).ne.0.0_wp_) ii=ii+1 end if idecr=-1 is=1 @@ -6148,8 +6737,8 @@ c ypt(i)=ppabs(j,k,i) yamp(i)=ccci(j,k,i) else - ypt(i)=0.0d0 - yamp(i)=0.0d0 + ypt(i)=0.0_wp_ + yamp(i)=0.0_wp_ end if if(ise0.eq.0) then if(xxi(i).lt.rtbc) then @@ -6167,7 +6756,7 @@ c else if(xxi(i).gt.rtbc) exit if(xxi(i).lt.xxi(i-1)) then -! isev(is)=i !!!!!!!!!! it should be isev(is)=i-1 +c isev(is)=i !!!!!!!!!! it should be isev(is)=i-1 isev(is)=i-1 is=is+1 idecr=-1 @@ -6177,8 +6766,8 @@ c end do c isev(is)=i-1 - ppa1=0.0d0 - cci1=0.0d0 + ppa1=0.0_wp_ + cci1=0.0_wp_ do iis=1,is-1 iis1=iis+1 iise0=isev(iis) @@ -6195,7 +6784,7 @@ c iind=1 end if do ind=ind1,ind2,iind !!!!!!!!!! is it safe? -! iind=ind !!!!!!!!!! iind reused in the loop! +c iind=ind !!!!!!!!!! iind reused in the loop! indi=ind !!!!!!!!!! iind reused in the loop! if (idecr.eq.-1) indi=ind-1 rt1=rtab1(indi) @@ -6217,12 +6806,12 @@ c end do end do - h=1.0d0/dble(nd-1) - rhotpav=0.0d0 - drhotpav=0.0d0 - rhotjav=0.0d0 - rhotjava=0.0d0 - rhot2java=0.0d0 + h=1.0_wp_/dble(nd-1) + rhotpav=0.0_wp_ + drhotpav=0.0_wp_ + rhotjav=0.0_wp_ + rhotjava=0.0_wp_ + rhot2java=0.0_wp_ fi=dpdv/h call simpson (nd,h,fi,spds) @@ -6252,12 +6841,12 @@ c c factor sqrt(8)=2 sqrt(2) to match with full width c of gaussian profile drhot2pav=rhot2pav-rhotpav**2 - drhotpav=sqrt(8.d0*drhot2pav) + drhotpav=sqrt(8.0_wp_*drhot2pav) drhot2java=rhot2java-rhotjava**2 - drhotjava=sqrt(8.d0*drhot2java) + drhotjava=sqrt(8.0_wp_*drhot2java) - spds=0.0d0 - sccs=0.0d0 + spds=0.0_wp_ + sccs=0.0_wp_ do i=1,nd spds=spds+dpdv(i) sccs=sccs+ajphiv(i) @@ -6267,10 +6856,10 @@ c of gaussian profile ajphiv(i)=ajphiv(i)/darea(i) end do - facpds=1.0d0 - facjs=1.0d0 - if(spds.gt.0.0d0) facpds=pabs/spds - if(sccs.ne.0.0d0) facjs=currt/sccs + facpds=1.0_wp_ + facjs=1.0_wp_ + if(spds.gt.0.0_wp_) facpds=pabs/spds + if(sccs.ne.0.0_wp_) facjs=currt/sccs do i=1,nd dpdv(i)=facpds*dpdv(i) @@ -6282,8 +6871,8 @@ c of gaussian profile rhpp=frhopol(rhotpav) rhpj=frhopol(rhotjava) - dpdvp=pabs*2.0d0/(sqrt(pi)*drhotpav*fdvdrhot(rhpp)) - ajphip=currt*2.0d0/(sqrt(pi)*drhotjava*fdadrhot(rhpj)) + dpdvp=pabs*2.0_wp_/(sqrt(pi)*drhotpav*fdvdrhot(rhpp)) + ajphip=currt*2.0_wp_/(sqrt(pi)*drhotjava*fdadrhot(rhpj)) call ratioj(rhpj,ratjamx,ratjbmx,ratjplmx) call profwidth(nd,rtab,dpdv,rhotp,rhopp,dpdvmx, @@ -6293,73 +6882,73 @@ c of gaussian profile . drhotjfi,drhopfi) xps=rhopfi else - rhotjfi=0.0d0 - rhopfi=0.0d0 - ajmxfi=0.0d0 - ajphip=0.0d0 - drhotjfi=0.0d0 - drhopfi=0.0d0 + rhotjfi=0.0_wp_ + rhopfi=0.0_wp_ + ajmxfi=0.0_wp_ + ajphip=0.0_wp_ + drhotjfi=0.0_wp_ + drhopfi=0.0_wp_ xps=rhopp - ratjamx=1.0d0 - ratjbmx=1.0d0 - ratjplmx=1.0d0 + ratjamx=1.0_wp_ + ratjbmx=1.0_wp_ + ratjplmx=1.0_wp_ end if iif1=iiv(1,1) istmx=1 do i=2,iif1 - if(psjki(1,1,i).ge.0.0d0) then + if(psjki(1,1,i).ge.0.0_wp_) then if(pdjki(1,1,i).gt.pdjki(1,1,i-1)) istmx=i end if end do stmx=istmx*dst - pins_02=0.0d0 - pins_05=0.0d0 - pins_085=0.0d0 + pins_02=0.0_wp_ + pins_05=0.0_wp_ + pins_085=0.0_wp_ - xrhot=0.2d0 + xrhot=0.2_wp_ call locate(rhotv,nd,xrhot,i1) call intlin(rhotv(i1),pins(i1),rhotv(i1+1),pins(i1+1), . xrhot,pins_02) - xrhot=0.5d0 + xrhot=0.5_wp_ call locate(rhotv,nd,xrhot,i1) call intlin(rhotv(i1),pins(i1),rhotv(i1+1),pins(i1+1), . xrhot,pins_05) - xrhot=0.85d0 + xrhot=0.85_wp_ call locate(rhotv,nd,xrhot,i1) call intlin(rhotv(i1),pins(i1),rhotv(i1+1),pins(i1+1), . xrhot,pins_085) else - ajmxfi=0.0d0 - ajphip=0.0d0 - dpdvp=0.0d0 - dpdvmx=0.0d0 - rhotjfi=1.0d0 - rhotjav=1.0d0 - rhotjava=1.0d0 - rhotp=1.0d0 - rhotpav=1.0d0 - drhotjfi=0.0d0 - drhotjava=0.0d0 - drhotp=0.0d0 - drhotpav=0.0d0 - ratjamx=1.0d0 - ratjbmx=1.0d0 + ajmxfi=0.0_wp_ + ajphip=0.0_wp_ + dpdvp=0.0_wp_ + dpdvmx=0.0_wp_ + rhotjfi=1.0_wp_ + rhotjav=1.0_wp_ + rhotjava=1.0_wp_ + rhotp=1.0_wp_ + rhotpav=1.0_wp_ + drhotjfi=0.0_wp_ + drhotjava=0.0_wp_ + drhotp=0.0_wp_ + drhotpav=0.0_wp_ + ratjamx=1.0_wp_ + ratjbmx=1.0_wp_ taumn=0 taumx=0 stmx=stf - pins_02=0.0d0 - pins_05=0.0d0 - pins_085=0.0d0 + pins_02=0.0_wp_ + pins_05=0.0_wp_ + pins_085=0.0_wp_ c end of pabstot > 0 end if c dPdV [MW/m^3], Jcd [MA/m^2] - if(ieccd.eq.0) currt=0.0d0 - currtka=currt*1.0d3 + if(ieccd.eq.0) currt=0.0_wp_ + currtka=currt*1.0e3_wp_ write(6,*)' ' write(6,*)'#beta0 alpha0 Icd Pa Jphip dPdVp '// @@ -6383,7 +6972,7 @@ c dPdV [MW/m^3], Jcd [MA/m^2] psin=rtab(i)**2 rhop=rtab(i) end if - pinsr=0.0d0 + pinsr=0.0_wp_ if(pabstot.gt.0) pinsr=pins(i)/pabstot write(48,99) psin,rhotv(i),ajphiv(i),ajcdbv(i),dpdv(i), . currins(i),pins(i),pinsr,real(index_rt) @@ -6396,34 +6985,44 @@ c c c subroutine profwidth(nd,xx,yy,rhotmx,rhopmx,ypk,drhot,drhop) - use const_and_precisions, only : emn1 + use const_and_precisions, only : wp_,emn1 use graydata_flags, only : ipec,iequil use utils, only : locatex, locate, intlin, vmaxmini - - implicit real*8(a-h,o-z) - dimension xx(nd),yy(nd) + implicit none +c arguments + integer :: nd + real(wp_) :: rhotmx,rhopmx,ypk,drhot,drhop + real(wp_), dimension(nd) :: xx,yy +c local variables + integer :: imn,imx,ipk,ie1,ie2 + real(wp_) :: xmn,xmx,ymn,ymx,xpkp,xpkm,xpk,yye,rte1,rte2, + . psie1,psie2,rhopmn,rhote2,rhote1,rhotmn +c common/external functions/variables + real(wp_) :: rhotp,rhotm,ypkp,ypkm + real(wp_) :: frhotor +c common/jmxmn/rhotp,rhotm,ypkp,ypkm c call vmaxmini(yy,nd,ymn,ymx,imn,imx) - ypk=0.0d0 + ypk=0.0_wp_ xmx=xx(imx) xmn=xx(imn) if (abs(ymx).gt.abs(ymn)) then ipk=imx ypkp=ymx xpkp=xmx - if(abs(ymn/ymx).lt.1d-2) ymn=0.0d0 + if(abs(ymn/ymx).lt.1.0e-2_wp_) ymn=0.0_wp_ ypkm=ymn xpkm=xmn else ipk=imn ypkp=ymn xpkp=xmn - if(abs(ymx/ymn).lt.1d-2) ymx=0.0d0 + if(abs(ymx/ymn).lt.1.0e-2_wp_) ymx=0.0_wp_ ypkm=ymx xpkm=xmx end if - if(xpkp.gt.0.0d0) then + if(xpkp.gt.0.0_wp_) then xpk=xpkp ypk=ypkp yye=ypk*emn1 @@ -6431,19 +7030,19 @@ c if(ie1.gt.0.and.ie1.lt.nd) then call intlin(yy(ie1),xx(ie1),yy(ie1+1),xx(ie1+1),yye,rte1) else - rte1=0.0d0 + rte1=0.0_wp_ end if call locatex(yy,nd,ipk,nd,yye,ie2) if(ie2.gt.0.and.ie2.lt.nd) then call intlin(yy(ie2),xx(ie2),yy(ie2+1),xx(ie2+1),yye,rte2) else - rte2=0.0d0 + rte2=0.0_wp_ end if else ipk=2 xpk=xx(2) ypk=yy(2) - rte1=0.0d0 + rte1=0.0_wp_ yye=ypk*emn1 call locate(yy,nd,yye,ie2) call intlin(yy(ie2),xx(ie2),yy(ie1+2),xx(ie1+2),yye,rte2) @@ -6452,12 +7051,12 @@ c ipk=1 if(ymx.ne.0.and.ymn.ne.0) ipk=2 c - drhop=0.0d0 - drhot=0.0d0 - psie1=0.0d0 - psie2=1.0d0 - rhopmx=1.0d0 - rhopmn=0.0d0 + drhop=0.0_wp_ + drhot=0.0_wp_ + psie1=0.0_wp_ + psie2=1.0_wp_ + rhopmx=1.0_wp_ + rhopmn=0.0_wp_ if (ie1.gt.0.or.ie2.gt.0) then if(ipec.eq.0) then rhopmx=sqrt(xpkp) @@ -6491,9 +7090,19 @@ c c return end +c +c c subroutine polarcold(exf,eyif,ezf,elf,etf) - implicit real*8(a-h,o-z) + use const_and_precisions, only : wp_ + implicit none +c arguments + real(wp_) :: exf,eyif,ezf,elf,etf +c local variables + real(wp_) :: anpl2,anpr2,an2,yg2,dy2,aa,e3,qq,p +c common/external functions/variables + real(wp_) anpl,anpr,xg,yg,sox +c common/nplr/anpl,anpr common/xgxg/xg common/ygyg/yg @@ -6502,49 +7111,49 @@ c c dcold dispersion c dielectric tensor (transposed) c -c exf=0.0d0 -c eyif=0.0d0 -c ezf=0.0d0 +c exf=0.0_wp_ +c eyif=0.0_wp_ +c ezf=0.0_wp_ c if(xg.le.0) return c anpl2=anpl*anpl anpr2=anpr*anpr an2=anpl2+anpr2 yg2=yg**2 - dy2=1.0d0-yg2 - aa=1.0d0-xg-yg2 - e3=1.0d0-xg + dy2=1.0_wp_-yg2 + aa=1.0_wp_-xg-yg2 + e3=1.0_wp_-xg c - if(xg.gt.0.0d0) then - if (anpl.ne.0.0d0) then + if(xg.gt.0.0_wp_) then + if (anpl.ne.0.0_wp_) then qq=xg*yg/(an2*dy2-aa) p=(anpr2-e3)/(anpl*anpr) - ezf=1.0d0/sqrt(1.0d0+p*p*(1.0d0+qq*qq)) + ezf=1.0_wp_/sqrt(1.0_wp_+p*p*(1.0_wp_+qq*qq)) exf=p*ezf eyif=qq*exf else - if(sox.lt.0.d0) then + if(sox.lt.0.0_wp_) then ezf=1 exf=0 eyif=0 else ezf=0 qq=-aa/(xg*yg) - exf=1.0d0/sqrt(1.0d0+qq*qq) + exf=1.0_wp_/sqrt(1.0_wp_+qq*qq) eyif=qq*exf end if end if elf=(anpl*ezf+anpr*exf)/sqrt(an2) - etf=sqrt(1.0d0-elf*elf) + etf=sqrt(1.0_wp_-elf*elf) else - if(sox.lt.0.0d0) then + if(sox.lt.0.0_wp_) then ezf=1 exf=0 eyif=0 else ezf=0 - exf=0.0d0 - eyif=1.0d0 + exf=0.0_wp_ + eyif=1.0_wp_ end if elf=0 etf=1 @@ -6552,17 +7161,22 @@ c c return end - +c +c +c subroutine pol_limit(ext,eyt) - use const_and_precisions, only : ui=>im,pi + use const_and_precisions, only : wp_,ui=>im,pi implicit none - real*8 bv(3),anv(3) - real*8 anx,any,anz,anpl,anpr,yg,xe2om,ye2om,xe2xm,ye2xm - real*8 an2,an,anxy,sngam,csgam,csg2,sng2,ffo,ffx,ffo2,ffx2 - real*8 deno,denx,anpl2,dnl,del0 - real*8 gam - real*8 sox - complex*16 exom,eyom,exxm,eyxm,ext,eyt +c arguments + complex(wp_) :: ext,eyt +c local variables + real(wp_) :: anx,any,anz,xe2om,ye2om,xe2xm,ye2xm,an2,an,anxy, + . sngam,csgam,csg2,sng2,ffo,ffx,ffo2,ffx2,deno,denx,anpl2,dnl, + . del0,gam + complex(wp_) :: exom,eyom,exxm,eyxm +c common/external functions/variables + real(wp_) :: anpl,anpr,yg,sox + real(wp_), dimension(3) :: bv,anv c common/anv/anv common/nplr/anpl,anpr @@ -6574,10 +7188,10 @@ c any=anv(2) anz=anv(3) anpl2=anpl*anpl - dnl=1.0d0-anpl2 - del0=sqrt(dnl**2+4.0d0*anpl2/yg**2) - ffo=0.5d0*yg*(dnl+del0) - ffx=0.5d0*yg*(dnl-del0) + dnl=1.0_wp_-anpl2 + del0=sqrt(dnl**2+4.0_wp_*anpl2/yg**2) + ffo=0.5_wp_*yg*(dnl+del0) + ffx=0.5_wp_*yg*(dnl-del0) an2=anx*anx+any*any+anz*anz an=sqrt(an2) anxy=sqrt(anx*anx+any*any) @@ -6600,7 +7214,7 @@ c exxm=(ffx*csgam-ui*anpl*sngam)/sqrt(denx) eyxm=(-ffx*sngam-ui*anpl*csgam)/sqrt(denx) c - if (sox.lt.0.0d0) then + if (sox.lt.0.0_wp_) then ext=exom eyt=eyom else @@ -6608,52 +7222,67 @@ c eyt=eyxm endif - gam=atan(sngam/csgam)*180.d0/pi + gam=atan(sngam/csgam)*180.0_wp_/pi return end - +c +c +c subroutine stokes(ext,eyt,qq,uu,vv) + use const_and_precisions, only : wp_ implicit none - complex*16 ext,eyt - real*8 qq,uu,vv +c arguments + complex(wp_) :: ext,eyt + real(wp_) :: qq,uu,vv +c qq=abs(ext)**2-abs(eyt)**2 - uu=2.0d0*dble(ext*dconjg(eyt)) - vv=2.0d0*dimag(ext*dconjg(eyt)) + uu=2.0_wp_*dble(ext*dconjg(eyt)) + vv=2.0_wp_*dimag(ext*dconjg(eyt)) +c end subroutine stokes - +c +c +c subroutine polellipse(qq,uu,vv,psipol,chipol) - use const_and_precisions, only : pi + use const_and_precisions, only : wp_,pi implicit none - real*8 qq,uu,vv,psipol,chipol +c arguments + real(wp_) :: qq,uu,vv,psipol,chipol +c c real*8 llm,aa,bb,ell c llm=sqrt(qq**2+uu**2) -c aa=sqrt((1+llm)/2.0d0) -c bb=sqrt((1-llm)/2.0d0) +c aa=sqrt((1+llm)/2.0_wp_) +c bb=sqrt((1-llm)/2.0_wp_) c ell=bb/aa - psipol=0.5d0*atan2(uu,qq)*180.d0/pi - chipol=0.5d0*asin(vv)*180.d0/pi + psipol=0.5_wp_*atan2(uu,qq)*180.0_wp_/pi + chipol=0.5_wp_*asin(vv)*180.0_wp_/pi +c end subroutine polellipse - - +c +c +c logical function inside_plasma(rrm,zzm) + use const_and_precisions, only : wp_ use graydata_flags, only : iequil use graydata_par, only : psdbnd use interp_eqprof, only : zbmin,zbmax - implicit none - real*8 rrm,zzm,psinv - +c arguments + real(wp_) :: rrm,zzm +c common/external functions/variables + real(wp_) :: psinv +c common/psival/psinv - +c if(iequil.eq.1) then call equian(rrm,zzm) else call equinum_psi(rrm,zzm) end if - - if (psinv.ge.0.0d0.and.psinv.lt.psdbnd) then - if (psinv.lt.1.0d0.and.zzm.lt.zbmin.or.zzm.gt.zbmax) then +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 inside_plasma=.false. else inside_plasma=.true. @@ -6661,31 +7290,35 @@ c ell=bb/aa else inside_plasma=.false. end if - +c end function inside_plasma - - +c +c +c subroutine wall_refl(xv,anv,ext,eyt,xvrfl,anvrfl,extr,eytr,walln, . irfl) + use const_and_precisions, only : wp_,ui=>im,pi use reflections, only : inters_linewall,inside - use const_and_precisions, only : ui=>im,pi use interp_eqprof, only : rlim,zlim,nlim implicit none - integer*4 irfl - real*8 anv(3),anv0(3),xv(3),xvrfl(3) - real*8 walln(3),anvrfl(3),vv1(3),vv2(3),vv3(3) - real*8 smax,rrm,zzm - complex*16 extr,eytr,eztr,ext,eyt - complex*16 evin(3),evrfl(3) - +c arguments + integer :: irfl + real(wp_), dimension(3) :: xv,anv,xvrfl,anvrfl,walln + complex(wp_) :: ext,eyt,extr,eytr +c local variables + real(wp_) :: smax,rrm,zzm + real(wp_), dimension(3) :: anv0,vv1,vv2,vv3 + complex(wp_) :: eztr + complex(wp_), dimension(3) :: evin,evrfl +c anv0=anv/sqrt(anv(1)**2+anv(2)**2+anv(3)**2) - rrm=1.d-2*sqrt(xv(1)**2+xv(2)**2) - zzm=1.d-2*xv(3) - + rrm=1.0e-2_wp_*sqrt(xv(1)**2+xv(2)**2) + zzm=1.0e-2_wp_*xv(3) +c c computation of reflection coordinates and normal to the wall - call inters_linewall(xv/1.d2,anv0,rlim(1:nlim),zlim(1:nlim), + call inters_linewall(xv/1.0e2_wp_,anv0,rlim(1:nlim),zlim(1:nlim), . nlim,smax,walln) - smax=smax*1.d2 + smax=smax*1.0e2_wp_ xvrfl=xv+smax*anv0 irfl=1 if (.not.inside(rlim,zlim,nlim,rrm,zzm)) then @@ -6700,82 +7333,90 @@ c computation of reflection coordinates and normal to the wall return end if ! search second wall interface (inside-outside) - call inters_linewall(xvrfl/1.d2,anv0,rlim(1:nlim),zlim(1:nlim), - . nlim,smax,walln) - smax=smax*1.d2 + call inters_linewall(xvrfl/1.0e2_wp_,anv0,rlim(1:nlim), + . zlim(1:nlim),nlim,smax,walln) + smax=smax*1.0e2_wp_ xvrfl=xvrfl+smax*anv0 irfl=2 end if - +c c rotation matrix from local to lab frame vv1(1)=anv0(2) vv1(2)=-anv0(1) - vv1(3)=0.0d0 + vv1(3)=0.0_wp_ vv2(1)=anv0(1)*anv0(3) vv2(2)=anv0(2)*anv0(3) vv2(3)=-anv0(1)*anv0(1)-anv0(2)*anv0(2) vv1=vv1/sqrt(vv1(1)**2+vv1(2)**2+vv1(3)**2) vv2=vv2/sqrt(vv2(1)**2+vv2(2)**2+vv2(3)**2) vv3=anv0 - +c evin=ext*vv1+eyt*vv2 c wave vector and electric field after reflection in lab frame - anvrfl=anv0-2.0d0* + anvrfl=anv0-2.0_wp_* . (anv0(1)*walln(1)+anv0(2)*walln(2)+anv0(3)*walln(3))*walln - evrfl=-evin+2.0d0* + evrfl=-evin+2.0_wp_* . (evin(1)*walln(1)+evin(2)*walln(2)+evin(3)*walln(3))*walln - +c vv1(1)=anvrfl(2) vv1(2)=-anvrfl(1) - vv1(3)=0.0d0 + vv1(3)=0.0_wp_ vv2(1)=anvrfl(1)*anvrfl(3) vv2(2)=anvrfl(2)*anvrfl(3) vv2(3)=-anvrfl(1)*anvrfl(1)-anvrfl(2)*anvrfl(2) vv1=vv1/sqrt(vv1(1)**2+vv1(2)**2+vv1(3)**2) vv2=vv2/sqrt(vv2(1)**2+vv2(2)**2+vv2(3)**2) vv3=anvrfl/sqrt(anvrfl(1)**2+anvrfl(2)**2+anvrfl(3)**2) - +c extr=dot_product(vv1,evrfl) eytr=dot_product(vv2,evrfl) eztr=dot_product(vv3,evrfl) - +c return end - +c +c +c subroutine vacuum_rt(xvstart,anv,xvend,ivac) + use const_and_precisions, only : wp_ use reflections, only : inters_linewall,inside use graydata_flags, only : dst use interp_eqprof, only : rlim,zlim,nlim - implicit none - integer*4 ivac - integer i - real*8 st,rrm,zzm,dstvac,smax - real*8 anv(3),xvstart(3),xvend(3),walln(3) - real*8 xv0(3),anv0(3) - logical plfound,inside_plasma +c arguments + integer :: ivac + real(wp_), dimension(3) :: xvstart,anv,xvend +c local variables + integer :: i + real(wp_) :: st,rrm,zzm,smax + real(wp_), dimension(3) :: walln,xv0,anv0 + logical :: plfound +c common/external functions/variables + real(wp_) :: dstvac + logical :: inside_plasma +c external inside_plasma - +c common/dstvac/dstvac - +c c ivac=1 plasma hit before wall reflection c ivac=2 wall hit before plasma c ivac=-1 vessel (and thus plasma) never crossed -! the real argument associated to xvstart is in a common block -! used by fwork!!! +c the real argument associated to xvstart is in a common block +c used by fwork!!! xv0=xvstart anv0=anv - call inters_linewall(xv0/1.d2,anv,rlim(1:nlim),zlim(1:nlim), + call inters_linewall(xv0/1.0e2_wp_,anv,rlim(1:nlim),zlim(1:nlim), . nlim,smax,walln) - smax=smax*1.d2 - rrm=1.d-2*sqrt(xv0(1)**2+xv0(2)**2) - zzm=1.d-2*xv0(3) + smax=smax*1.0e2_wp_ + rrm=1.0e-2_wp_*sqrt(xv0(1)**2+xv0(2)**2) + zzm=1.0e-2_wp_*xv0(3) if (.not.inside(rlim,zlim,nlim,rrm,zzm)) then ! first wall interface is outside-inside if (dot_product(walln,walln)