From e31f05e9a80d1b1f93a272f3732a9fa93597433f Mon Sep 17 00:00:00 2001 From: Lorenzo Figini Date: Fri, 15 May 2015 16:02:27 +0000 Subject: [PATCH] moved all the definitions of arithmetical, physical and type kind constants in const_and_precisions. Changed argument from tekev to mu in spitzer function subroutines. --- Makefile | 11 ++- src/const_and_precisions.f90 | 107 ++++++++++++++++------- src/gray.f | 164 ++++++++++++++++------------------- src/green_func_p.f90 | 32 ++++--- src/itm_constants.f90 | 32 ------- src/itm_types.f90 | 50 ----------- src/reflections.f90 | 111 ++++++++++++------------ 7 files changed, 228 insertions(+), 279 deletions(-) delete mode 100644 src/itm_constants.f90 delete mode 100644 src/itm_types.f90 diff --git a/Makefile b/Makefile index 43cc315..a4080fe 100644 --- a/Makefile +++ b/Makefile @@ -4,7 +4,7 @@ EXE=gray # Objects list MAINOBJ=gray.o OTHOBJ=dqagmv.o grayl.o reflections.o green_func_p.o \ - const_and_precisions.o itm_constants.o itm_types.o + const_and_precisions.o # Alternative search paths vpath %.f90 src @@ -23,17 +23,16 @@ $(EXE): $(MAINOBJ) $(OTHOBJ) $(FC) $(FFLAGS) -o $@ $^ # Dependencies on modules -gray.o: dqagmv.o green_func_p.o reflections.o +gray.o: dqagmv.o green_func_p.o reflections.o const_and_precisions.o green_func_p.o: const_and_precisions.o -const_and_precisions.o: itm_types.o itm_constants.o -itm_constants.o: itm_types.o +reflections.o: const_and_precisions.o # General object compilation command %.o: %.f90 - $(FC) $(FFLAGS) -c $^ + $(FC) $(FFLAGS) -c $< %.o: %.f - $(FC) $(FFLAGS) -c $^ + $(FC) $(FFLAGS) -c $< gray.o:gray.f $(FC) -cpp $(DIRECTIVES) $(FFLAGS) -c $< diff --git a/src/const_and_precisions.f90 b/src/const_and_precisions.f90 index 26ce054..048ae6c 100755 --- a/src/const_and_precisions.f90 +++ b/src/const_and_precisions.f90 @@ -1,17 +1,21 @@ !########################################################################! MODULE const_and_precisions - use itm_types, only : wp_ => r8 - use itm_constants, only : pi => itm_pi, e_ => itm_qe, me_ => itm_me, c_ => itm_c !########################################################################! IMPLICIT NONE PUBLIC !------------------------------------------------------------------------ ! common precisions !------------------------------------------------------------------------ -! INTEGER, PARAMETER :: sp_ = 4 ! single precision -! INTEGER, PARAMETER :: dp_ = 8 ! double precision -! INTEGER, PARAMETER :: wp_ = dp_ ! work-precision +! INTEGER, PARAMETER :: i1 = SELECTED_INT_KIND (2) ! Integer*1 +! INTEGER, PARAMETER :: i2 = SELECTED_INT_KIND (4) ! Integer*2 + INTEGER, PARAMETER :: i4 = SELECTED_INT_KIND (9) ! Integer*4 + INTEGER, PARAMETER :: i8 = SELECTED_INT_KIND (18) ! Integer*8 + INTEGER, PARAMETER :: r4 = SELECTED_REAL_KIND (6, 37) ! Real*4 + INTEGER, PARAMETER :: r8 = SELECTED_REAL_KIND (15, 300) ! Real*8 +! INTEGER, PARAMETER :: sp_ = r4 ! single precision +! INTEGER, PARAMETER :: dp_ = r8 ! double precision + INTEGER, PARAMETER :: wp_ = r8 ! work-precision ! INTEGER, PARAMETER :: odep_ = dp_ ! ODE-solver precision ! INTEGER, PARAMETER :: xp_ = wp_ ! for ext. modules if necessary !------------------------------------------------------------------------ @@ -26,31 +30,33 @@ !!======================================================================== ! Arithmetic constants !======================================================================== + integer, parameter :: izero = 0 REAL(wp_), PARAMETER :: zero = 0.0_wp_ - REAL(wp_), PARAMETER :: unit = 1.0_wp_ -! REAL(wp_), PARAMETER :: pi = 3.141592653589793_wp_ -! REAL(wp_), PARAMETER :: sqrt_pi = 1.772453850905516_wp_ + REAL(wp_), PARAMETER :: one = 1.0_wp_ + real(wp_), parameter :: pi = 3.141592653589793_wp_ ! 3.141592653589793238462643383280 + REAL(wp_), PARAMETER :: sqrt_pi = 1.772453850905516_wp_ ! REAL(wp_), PARAMETER :: sqrt_2 = 1.414213562373095_wp_ -! REAL(wp_), PARAMETER :: rad = pi/180.0_wp_ + REAL(wp_), PARAMETER :: degree = pi/180.0_wp_ + REAL(wp_), PARAMETER :: emn1 = 0.367879441171442_wp_ ! exp(-1) !--- -! REAL(wp_), PARAMETER :: ex(1:3) = (/unit,zero,zero/) -! REAL(wp_), PARAMETER :: ey(1:3) = (/zero,unit,zero/) -! REAL(wp_), PARAMETER :: ez(1:3) = (/zero,zero,unit/) +! REAL(wp_), PARAMETER :: ex(1:3) = (/one ,zero,zero/) +! REAL(wp_), PARAMETER :: ey(1:3) = (/zero,one ,zero/) +! REAL(wp_), PARAMETER :: ez(1:3) = (/zero,zero,one /) !--- -! REAL(wp_), PARAMETER :: kron(3,3) = reshape((/unit,zero,zero, & -! zero,unit,zero, & -! zero,zero,unit/),(/3,3/)) -! COMPLEX(wp_), PARAMETER :: im = (0.0_wp_,1.0_wp_) -! COMPLEX(wp_), PARAMETER :: czero = (0.0_wp_,0.0_wp_) -! COMPLEX(wp_), PARAMETER :: cunit = (1.0_wp_,0.0_wp_) +! REAL(wp_), PARAMETER :: kron(3,3) = reshape((/one ,zero,zero, & +! zero,one ,zero, & +! zero,zero,one /),(/3,3/)) + COMPLEX(wp_), PARAMETER :: im = (0.0_wp_,1.0_wp_) + COMPLEX(wp_), PARAMETER :: czero = (0.0_wp_,0.0_wp_) + COMPLEX(wp_), PARAMETER :: cunit = (1.0_wp_,0.0_wp_) ! COMPLEX(wp_), PARAMETER :: ctwo = (2.0_wp_,0.0_wp_) !======================================================================== ! Computer constants !======================================================================== - REAL(wp_), PARAMETER :: comp_eps = EPSILON(unit) -! REAL(wp_), PARAMETER :: comp_eps2 = comp_eps**2 -! REAL(wp_), PARAMETER :: comp_tiny = TINY(unit) -! REAL(wp_), PARAMETER :: comp_huge = HUGE(unit) + REAL(wp_), PARAMETER :: comp_eps = EPSILON(one) +! REAL(wp_), PARAMETER :: comp_eps2 = comp_eps**2 + REAL(wp_), PARAMETER :: comp_tiny = TINY(one) + REAL(wp_), PARAMETER :: comp_huge = HUGE(one) ! REAL(wp_), PARAMETER :: comp_tinylog =-200 ! LOG10(comp_tiny) ! REAL(wp_), PARAMETER :: comp_hugelog =+200 ! LOG10(comp_huge) ! REAL(wp_), PARAMETER :: comp_tiny1 = 1d+50*comp_tiny @@ -60,23 +66,39 @@ !------------------------------------------------------------------------ ! Conventional constants !------------------------------------------------------------------------ + INTEGER, PARAMETER :: int_invalid = -999999999 + REAL(R8), PARAMETER :: r8_invalid = -9.0e40_r8 ! REAL(wp_), PARAMETER :: output_tiny = 1.0d-66 ! REAL(wp_), PARAMETER :: output_huge = 1.0d+66 !======================================================================== ! Physical constants (SI) !======================================================================== -! REAL(wp_), PARAMETER :: e_ = 1.602176487d-19 ! [C] -! REAL(wp_), PARAMETER :: me_ = 9.10938215d-31 ! [kg] -! REAL(wp_), PARAMETER :: mp_ = 1.672621637d-27 ! [kg] -! REAL(wp_), PARAMETER :: rmpe_ = mp_/me_ -! REAL(wp_), PARAMETER :: c_ = 2.99792458d+08 ! [m/s] -! REAL(wp_), PARAMETER :: eps0_ = 8.854187817d-12 ! [F/m] + real (wp_), parameter :: e_ = 1.602176487e-19_wp_ ! elementary charge, C + real (wp_), parameter :: me_ = 9.10938215e-31_wp_ ! electron mass, kg +! real (wp_), parameter :: mp_ = 1.672621637e-27_wp_ ! proton mass, kg +! real (wp_), parameter :: md_ = 3.34358320e-27_wp_ ! deuteron mass, kg +! real (wp_), parameter :: mt_ = 5.00735588e-27_wp_ ! triton mass, kg +! real (wp_), parameter :: ma_ = 6.64465620e-27_wp_ ! alpha mass, kg +! real (wp_), parameter :: amu_ = 1.660538782e-27_wp_ ! amu, kg +! REAL (wp_), PARAMETER :: rmpe_ = mp_/me_ ! proton-electron mass ratio + real (wp_), parameter :: c_ = 2.99792458e8_wp_ ! speed of light, m/s + real (wp_), parameter :: mu0_ = 4.0e-7_wp_ * pi ! magnetic permeability of vacuum + real (wp_), parameter :: eps0_ = 1.0_wp_ / (mu0_ * c_**2) ! dielectric constant of vacuum, F/m +! real (wp_), parameter :: avogr = 6.02214179e23_wp_ +! real (wp_), parameter :: KBolt = 1.3806504e-23_wp_ +!======================================================================== +! Physical constants (cgs) +!======================================================================== + real (wp_), parameter :: ccgs_ = c_*1.e2_wp_ ! speed of light, cm/s + real (wp_), parameter :: mecgs_ = me_*1.e3_wp_ ! electron mass, g + real (wp_), parameter :: ecgs_ = e_*c_*10._wp_ ! elementary charge, statcoul !------------------------------------------------------------------------ ! Useful definitions !------------------------------------------------------------------------ - REAL(wp_), PARAMETER :: keV_ = 1000*e_ ! [J] + REAL(wp_), PARAMETER :: keV_ = 1.e3_wp_*e_ ! [J] REAL(wp_), PARAMETER :: mc2_SI = me_*c_**2 ! [J] REAL(wp_), PARAMETER :: mc2_ = mc2_SI/keV_ ! [keV] + REAL(wp_), PARAMETER :: mu0inv = 1._wp_/mu0_ ! ! REAL(wp_), PARAMETER :: mc_ = me_*c_ ! [kg*m/s] ! ! f_ce = fce1_*B (B in Tesla): ! ! REAL(wp_), PARAMETER :: wce1_ = e_/me_ ! [rad/s] @@ -100,6 +122,33 @@ ! REAL(wp_), PARAMETER :: Npar_min = 1.0d-3 !########################################################################! + interface is_valid + module procedure is_valid_int4, is_valid_int8, is_valid_real8 + end interface + +contains + + logical function is_valid_int4(in_int) + implicit none + integer(i4), intent(in) :: in_int + is_valid_int4 = in_int /= int_invalid + return + end function is_valid_int4 + + logical function is_valid_int8(in_int) + implicit none + integer(i8), intent(in) :: in_int + is_valid_int8 = in_int /= int_invalid + return + end function is_valid_int8 + + logical function is_valid_real8(in_real) + implicit none + real(r8), intent(in) :: in_real + is_valid_real8 = abs(in_real - r8_invalid) > abs(r8_invalid) * 1.0e-15_r8 + return + end function is_valid_real8 + END MODULE const_and_precisions !########################################################################! diff --git a/src/gray.f b/src/gray.f index 86c18cc..1257aaf 100644 --- a/src/gray.f +++ b/src/gray.f @@ -109,8 +109,8 @@ c ray integration: end subroutine after_gray_integration + use const_and_precisions, only : zero implicit real*8 (a-h,o-z) - parameter(zero=0.0d0) character*255 filenmeqq,filenmprf,filenmbm c common/ss/st @@ -176,11 +176,12 @@ c c c subroutine after_onestep(i,istop) + use const_and_precisions, only : pi, cvdr=>degree use reflections, only : inside implicit real*8 (a-h,o-z) complex*16 ext,eyt,extr,eytr,exin2,eyin2 parameter(jmx=31,kmx=36,nmx=8000,nbb=5000) - parameter(taucr=12.0d0,pi=3.14159265358979d0,cvdr=pi/180.0d0) + 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) @@ -471,10 +472,10 @@ c c c subroutine print_output(i,j,k) + use const_and_precisions, only : pi implicit real*8 (a-h,o-z) parameter(ndim=6,jmx=31,kmx=36,nmx=8000) parameter(taucr=12.0d0) - parameter(pi=3.14159265358979d0) dimension ywrk(ndim,jmx,kmx),ypwrk(ndim,jmx,kmx) dimension psjki(jmx,kmx,nmx) dimension tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx) @@ -648,14 +649,13 @@ c c c subroutine read_data + use const_and_precisions, only : qe=>ecgs_,me=>mecgs_,vc=>ccgs_, + * cvdr=>degree,pi use green_func_p, only:Setup_SpitzFunc implicit real*8 (a-h,o-z) - real*8 me character*8 wdat character*10 wtim character*255 filenmeqq,filenmprf,filenmbm - parameter(qe=4.8032d-10,me=9.1095d-28,vc=2.9979d+10) - parameter(pi=3.14159265358979d0,cvdr=pi/180.0d0) parameter(nmx=8000,nbb=5000) real*8 rlim(nbb),zlim(nbb) c @@ -1031,8 +1031,8 @@ c c c subroutine surf_anal + use const_and_precisions, only : pi implicit real*8(a-h,o-z) - parameter(pi=3.14159265358979d0) common/parban/b0,rr0m,zr0m,rpam common/parbres/bres c @@ -1183,9 +1183,9 @@ c c c subroutine equidata + use const_and_precisions, only : pi implicit real*8 (a-h,o-z) parameter(nnw=501,nnh=501) - parameter(pi=3.14159265358979d0) parameter(nbb=5000) c parameter(np=100) character*48 stringa @@ -2270,8 +2270,8 @@ c c c subroutine contours_psi(h,np,rcn,zcn,ipr) + use const_and_precisions, only : pi implicit real*8 (a-h,o-z) - parameter(pi=3.14159265358979d0) parameter(mest=4,kspl=3) parameter(nnw=501,nnh=501) parameter(nrest=nnw+4,nzest=nnh+4) @@ -2335,12 +2335,11 @@ c c c subroutine flux_average + use const_and_precisions, only : zero,one,pi,ccj=>mu0inv implicit real*8 (a-h,o-z) real*8 lam parameter(nnintp=101,ncnt=100,ncntt=2*ncnt+1,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) @@ -3404,8 +3403,8 @@ c c c subroutine plas_deriv(xx,yy,zz) + use const_and_precisions, only : pi implicit real*8 (a-h,o-z) - parameter(pi=3.14159265358979d0) dimension bv(3),derxg(3),deryg(3),derbv(3,3),dbtot(3) dimension bvc(3),dbvcdc(3,3),dbvdc(3,3),dbv(3,3) c @@ -3786,8 +3785,8 @@ c c subroutine tor_curr(rpsim,zpsim,ajphi) + use const_and_precisions, only : pi,ccj=>mu0inv implicit real*8 (a-h,o-z) - parameter(pi=3.14159265358979d0,ccj=1.0d+7/(4.0d0*pi)) common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv call equinum(rpsim,zpsim) bzz= dpsidr/rpsim @@ -3960,10 +3959,10 @@ c c beam tracing initial conditions igrad=1 c subroutine ic_gb + use const_and_precisions, only : izero,zero,one,pi,cvdr=>degree implicit real*8 (a-h,o-z) parameter(ndim=6,ndimm=3) - parameter(jmx=31,kmx=36,zero=0.0d0,izero=0,one=1.0d0) - parameter(pi=3.14159265358979d0,cvdr=pi/180.0d0) + 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) @@ -4275,12 +4274,11 @@ c c ray tracing initial conditions igrad=0 c subroutine ic_rt + use const_and_precisions, only : izero,zero,one,pi,cvdr=>degree, + * ui=>im implicit real*8 (a-h,o-z) - complex*16 ui parameter(ndim=6,ndimm=3) - parameter(jmx=31,kmx=36,zero=0.0d0,izero=0,one=1.0d0) - parameter(pi=3.14159265358979d0,cvdr=pi/180.0d0) - parameter(ui=(0.0d0,1.0d0)) + 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) @@ -4472,10 +4470,10 @@ c subroutine ic_rt2 + use const_and_precisions, only : izero,zero,one,pi,cvdr=>degree implicit real*8 (a-h,o-z) parameter(ndim=6,ndimm=3) - parameter(jmx=31,kmx=36,zero=0.0d0,izero=0,one=1.0d0) - parameter(pi=3.14159265358979d0,cvdr=pi/180.0d0) + parameter(jmx=31,kmx=36) dimension ywrk0(ndim,jmx,kmx),ypwrk0(ndim,jmx,kmx) dimension ytmp(ndim),yptmp(ndim) dimension yyrfl(jmx,kmx,ndim) @@ -4685,9 +4683,9 @@ c c c subroutine pabs_curr(i,j,k) + use const_and_precisions, only : pi implicit real*8 (a-h,o-z) parameter(jmx=31,kmx=36,nmx=8000) - parameter(pi=3.14159265358979d0) dimension psjki(jmx,kmx,nmx) dimension tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx) dimension pdjki(jmx,kmx,nmx),ppabs(jmx,kmx,nmx) @@ -4764,12 +4762,11 @@ c c c subroutine ecrh_cd + use const_and_precisions, only : eqsi=>e_,mesi=>me_,vcsi=>c_, + * mc2=>mc2_ implicit real*8 (a-h,o-z) - real*8 mc2,mesi parameter(taucr=12.0d0,xxcr=16.0d0,eps=1.d-8) - parameter(qesi=1.602176487d-19,mesi=9.10938215d-31) - parameter(vcsi=2.99792458d+8) - parameter(mc2 = mesi*vcsi*vcsi/qesi*1d-3) + complex*16 ex,ey,ez c common/nharm/nharm,nhf common/warm/iwarm,ilarm @@ -4781,16 +4778,22 @@ c common/parwv/ak0,akinv,fhz 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 common/amut/amu common/zz/Zeff + common/btot/btot + common/bmxmn/bmax,bmin + common/fc/fc c c absorption computation c @@ -4852,8 +4855,12 @@ c ratiovgr=2.0d0*anpr/derdnm*vgm c ithn=1 if(lrm.gt.nharm) ithn=2 - if(ieccd.gt.0) call eccd(yg,anpl,anprre,amu,zeff, - * ieccd,nharm,nhf,ithn,effjcd,iokhawa,ierr) + if(ieccd.gt.0) then + rbn=btot/bmin + rbx=btot/bmax + call eccd(yg,anpl,anprre,amu,zeff,rbn,rbx,fc,ex,ey,ez, + * dens,psinv,ieccd,nharm,nhf,ithn,effjcd,iokhawa,ierr) + end if return 999 format(12(1x,e12.5)) @@ -5011,13 +5018,11 @@ c computation of dielectric tensor elements c up to third order in Larmor radius for hermitian part c subroutine diel_tens_fr(e330,epsl,lrm) + use const_and_precisions, only : pi,rpi=>sqrt_pi,ui=>im implicit real*8(a-h,o-z) - complex*16 ui complex*16 e330,epsl(3,3,lrm) complex*16 ca11,ca12,ca22,ca13,ca23,ca33 complex*16 cq0p,cq0m,cq1p,cq1m,cq2p - parameter(pi=3.14159265358979d0,rpi=1.77245385090552d0) - parameter(ui=(0.0d0,1.0d0)) dimension rr(-lrm:lrm,0:2,0:lrm),ri(lrm,0:2,lrm) c common/xgxg/xg @@ -5455,10 +5460,9 @@ c c c subroutine antihermitian(ri,lrm) + use const_and_precisions, only : rpi=>sqrt_pi implicit none integer lmx,nmx,lrm,n,k,m,mm - real*8 rpi - parameter(rpi=1.77245385090552d0) parameter(lmx=20,nmx=lmx+2) real*8 fsbi(nmx) real*8 ri(lrm,0:2,lrm) @@ -5581,10 +5585,10 @@ c computation of dielectric tensor elements c Krivenki and Orefice, JPP 30,125 (1983) c subroutine diel_tens_wr(ce330,cepsl,lrm) + use const_and_precisions, only : cui=>im implicit real*8 (a-b,d-h,o-z) implicit complex*16 (c) dimension cefp(0:lrm,0:2),cefm(0:lrm,0:2),cepsl(3,3,lrm) - parameter(cui=(0.0d0,1.0d0)) c common/xgxg/xg common/ygyg/yg @@ -5651,10 +5655,10 @@ c c c subroutine fsup(cefp,cefm,lrm) + use const_and_precisions, only : cui=>im implicit real*8 (a-b,d-h,o-z) implicit complex*16 (c) parameter(apsicr=0.7d0) - parameter(cui=(0.0d0,1.0d0)) dimension cefp(0:lrm,0:2),cefm(0:lrm,0:2) c common/ygyg/yg @@ -5796,42 +5800,32 @@ c - subroutine eccd(yg,anpl,anprre,amu,zeff, - * ieccd,nhmn,nhmx,ithn,effjcd,iokhawa,ierr) - use green_func_p, only: SpitzFuncCoeff,pi + subroutine eccd(yg,anpl,anprre,amu,zeff,rbn,rbx,fc,ex,ey,ez, + * dens,psinv,ieccd,nhmn,nhmx,ithn,effjcd,iokhawa,ierr) + 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_ implicit none - real*8 anum,denom,resp,resj,effjcd,coullog,dens,tekev + real*8 anum,denom,resp,resj,effjcd,coullog,dens real*8 yg,anpl,anpr,amu,anprre,anprim - real*8 mc2,anucc,canucc,ddens - real*8 qesi,mesi,vcsi,qe,me,vc,ceff,Zeff,psinv - real*8 rbn,rbx,btot,bmin,bmax,alams,fp0s,pa,fc + real*8 mc2m2,anucc,canucc,ddens + real*8 ceff,Zeff,psinv + real*8 rbn,rbx,alams,fp0s,pa,fc real*8 fjch,fjncl,fjch0,fconic real*8 cst2,eccdpar(5) complex*16 ex,ey,ez integer ieccd,nhmn,nhmx,nhn,ithn,ierr,iokhawa external fjch,fjncl,fjch0 - parameter(qesi=1.602176487d-19,mesi=9.10938215d-31) - parameter(vcsi=2.99792458d+8) - parameter(qe=4.8032068d-10,me=mesi*1.d3,vc=vcsi*1.d2) - parameter(mc2=mesi*vcsi*vcsi/qesi*1d-3) + parameter(mc2m2=1.0d0/mc2**2) parameter(canucc=2.0d13*pi*qe**4/(me**2*vc**3)) parameter(ceff=qesi/(mesi*vcsi)) - common/psival/psinv - - common/dens/dens,ddens - common/epolar/ex,ey,ez - common/btot/btot - common/bmxmn/bmax,bmin - common/fc/fc - anum=0.0d0 denom=0.0d0 effjcd=0.0d0 - tekev=mc2/amu - coullog=48.0d0-log(1.0d7*dens/tekev**2) + coullog=48.0d0-log(1.0d7*dens*mc2m2*amu**2) anucc=canucc*dens*coullog c nhmx=nhmn @@ -5848,11 +5842,9 @@ c cst2=1.0d0-B/B_max c alams=sqrt(1-B_min/B_max) c Zeff < 31 !!! c fp0s= P_a (alams) - rbn=btot/bmin - rbx=btot/bmax cst2=1.0d0-rbx if(cst2.lt.1d-6) cst2=0.0d0 - alams=sqrt(1.0d0-bmin/bmax) + alams=sqrt(1.0d0-rbx/rbn) pa=sqrt(32.0d0/(Zeff+1.0d0)-1.0d0)/2.0d0 fp0s=fconic(alams,pa,0) eccdpar(2)=rbn @@ -5879,18 +5871,15 @@ c fp0s= P_a (alams) c neoclassical model: c ft=1-fc trapped particle fraction c rzfc=(1+Zeff)/fc - rbn=btot/bmin - rbx=btot/bmax cst2=1.0d0-rbx if(cst2.lt.1d-6) cst2=0.0d0 - call SpitzFuncCoeff(Tekev,Zeff,fc) - eccdpar(2) = tekev - eccdpar(3) = fc - eccdpar(4) = rbx - eccdpar(5) = psinv + 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,5,resj,resp,iokhawa,ierr) + * fjncl,eccdpar,4,resj,resp,iokhawa,ierr) anum=anum+resj denom=denom+resp end do @@ -6071,9 +6060,9 @@ c extrapar(13) = uplp c extrapar(14) = uplm c extrapar(15) = ygn c + use const_and_precisions, only : ui=>im implicit real*8 (a-h,o-z) - complex*16 ex,ey,ez,ui,emxy,epxy - parameter(ui=(0.0d0,1.0d0)) + complex*16 ex,ey,ez,emxy,epxy dimension extrapar(npar) yg=extrapar(1) @@ -6293,10 +6282,9 @@ c extrapar(14) = uplm c extrapar(15) = ygn c c extrapar(16) = zeff -c extrapar(17) = tekev -c extrapar(18) = fc -c extrapar(19) = hb -c extrapar(20) = psin +c extrapar(17) = fc +c extrapar(18) = hb +c extrapar(29) = psin c use green_func_p, only: GenSpitzFunc implicit real*8 (a-h,o-z) @@ -6306,10 +6294,9 @@ c amu=extrapar(3) ygn=extrapar(15) zeff=extrapar(16) - tekev=extrapar(17) - fc=extrapar(18) - hb=extrapar(19) - psin=extrapar(20) + fc=extrapar(17) + hb=extrapar(18) + psin=extrapar(19) gam=anpl*upl+ygn u2=gam*gam-1.0d0 @@ -6317,7 +6304,7 @@ c upr2=u2-upl*upl bth=sqrt(2.0d0/amu) uth=u/bth - call GenSpitzFunc(Tekev,Zeff,fc,uth,u,gam,fk,dfk) + call GenSpitzFunc(amu,Zeff,fc,uth,u,gam,fk,dfk) fk=fk*(4.0d0/amu**2) dfk=dfk*(2.0d0/amu)*bth @@ -6450,10 +6437,10 @@ c c c subroutine pec(pabs,currt) + use const_and_precisions, only : pi,one implicit real*8(a-h,o-z) parameter(nndmx=5001,jmx=31,kmx=36,nmx=8000) - parameter(rtbc=1.0d0) - parameter(pi=3.14159265358979d0) + 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) @@ -6810,8 +6797,8 @@ c c c subroutine profwidth(nd,xx,yy,rhotmx,rhopmx,ypk,drhot,drhop) + use const_and_precisions, only : emn1 implicit real*8(a-h,o-z) - parameter(emn1=0.367879441171442d0) dimension xx(nd),yy(nd) common/jmxmn/rhotp,rhotm,ypkp,ypkm common/ipec/ipec,nnd @@ -6973,15 +6960,15 @@ c end subroutine pol_limit(ext,eyt) + use const_and_precisions, only : 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 pi,beta0,alpha0,gam + real*8 beta0,alpha0,gam real*8 sox - complex*16 ui,exom,eyom,exxm,eyxm,ext,eyt - parameter(ui=(0.0d0,1.0d0),pi=3.14159265358979d0) + complex*16 exom,eyom,exxm,eyxm,ext,eyt c common/anv/anv common/nplr/anpl,anpr @@ -7043,11 +7030,10 @@ c end subroutine stokes subroutine polellipse(qq,uu,vv,psipol,chipol) + use const_and_precisions, only : pi implicit none real*8 qq,uu,vv,psipol,chipol c real*8 llm,aa,bb,ell - real*8 pi - parameter(pi=3.14159265358979d0) c llm=sqrt(qq**2+uu**2) c aa=sqrt((1+llm)/2.0d0) c bb=sqrt((1-llm)/2.0d0) @@ -7089,14 +7075,14 @@ c ell=bb/aa subroutine wall_refl(xv,anv,ext,eyt,xvrfl,anvrfl,extr,eytr,walln, . irfl) use reflections, only : inters_linewall,inside + use const_and_precisions, only : ui=>im,pi 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 pi,smax,rrm,zzm - complex*16 ui,extr,eytr,eztr,ext,eyt + real*8 smax,rrm,zzm + complex*16 extr,eytr,eztr,ext,eyt complex*16 evin(3),evrfl(3) - parameter(ui=(0.0d0,1.0d0),pi=3.14159265358979d0) integer nbb,nlim parameter(nbb=5000) real*8 rlim(nbb),zlim(nbb) diff --git a/src/green_func_p.f90 b/src/green_func_p.f90 index 2fb61ee..e222928 100644 --- a/src/green_func_p.f90 +++ b/src/green_func_p.f90 @@ -107,8 +107,7 @@ END SUBROUTINE Setup_SpitzFunc - SUBROUTINE GenSpitzFunc(Te,Zeff,fc,u,q,gam, K,dKdu) - + SUBROUTINE GenSpitzFunc(mu,Zeff,fc,u,q,gam, K,dKdu) !======================================================================= ! Author: N.B.Marushchenko ! June 2005: as start point the subroutine of Ugo Gasparino (198?) @@ -172,15 +171,14 @@ ! dKdu = dK/du, i.e. its derivative over normalized momentum !======================================================================= IMPLICIT NONE - REAL(wp_), INTENT(in) :: Te,Zeff,fc,u,q,gam + REAL(wp_), INTENT(in) :: mu,Zeff,fc,u,q,gam REAL(wp_), INTENT(out) :: K,dKdu - REAL(wp_) :: mu,gam1,gam2,gam3,w,dwdu + REAL(wp_) :: gam1,gam2,gam3,w,dwdu !======================================================================= K = 0 dKdu = 0 IF (u < comp_eps) RETURN !--- - mu = mc2_/max(Te,1d-3) SELECT CASE(adj_appr(2)) CASE('m') !--------------- momentum conservation ------------------! gam1 = gam ! @@ -210,7 +208,7 @@ !####################################################################### !####################################################################### - SUBROUTINE SpitzFuncCoeff(Te,Zeff,fc) + SUBROUTINE SpitzFuncCoeff(mu,Zeff,fc) !======================================================================= ! Calculates the matrix coefficients required for the subroutine ! "GenSpitzFunc", where the Spitzer function is defined through the @@ -228,7 +226,7 @@ ! INPUT VARIABLES: ! rho = sqrt(SS) with SS - flux-surface label (norm. magn. flux) ! ne - density, 1/m^3 -! Te - temperature, keV +! mu - mc2/Te ! Zeff - effective charge ! fc - fraction of circulating particles ! @@ -237,9 +235,9 @@ ! "Spitzer"-function (the same as in the Hirshman paper) !======================================================================= IMPLICIT NONE - REAL(wp_), INTENT(in) :: Te,Zeff,fc + REAL(wp_), INTENT(in) :: mu,Zeff,fc INTEGER :: n,i,j - REAL(wp_) :: rtc,rtc1,mu,y,tn(1:nre) + REAL(wp_) :: rtc,rtc1,y,tn(1:nre) REAL(wp_) :: m(0:4,0:4),g(0:4) REAL(wp_) :: om(0:4,0:4) REAL(wp_) :: gam11,gam21,gam31,gam41,gam01, & @@ -250,17 +248,17 @@ alp23,alp24,alp20, & alp34,alp30,alp40 REAL(wp_) :: bet0,bet1,bet2,bet3,bet4,d0 - LOGICAL :: renew,rel,newTe,newne,newZ,newfc + LOGICAL :: renew,rel,newmu,newne,newZ,newfc REAL(wp_), SAVE :: sfdx(1:4) = 0 - REAL(wp_), SAVE :: ne_old =-1, Te_old =-1, Zeff_old =-1, fc_old =-1 + REAL(wp_), SAVE :: ne_old =-1, mu_old =-1, Zeff_old =-1, fc_old =-1 !======================================================================= - rel = Te > 1 - newTe = abs(Te -Te_old ) > delta*Te + rel = mu < mc2_ + newmu = abs(mu -mu_old ) > delta*mu newZ = abs(Zeff-Zeff_old) > delta*Zeff newfc = abs(fc -fc_old ) > delta*fc SELECT CASE(adj_appr(1)) CASE ('l','c') - renew = (newTe .and. rel) .OR. newZ .OR. newfc + renew = (newmu .and. rel) .OR. newZ .OR. newfc END SELECT !--- IF (.not.renew) THEN @@ -271,7 +269,7 @@ tn(:) = 0 IF (adj_appr(4) == 'r') THEN IF (nre > 0) THEN - mu = mc2_/max(Te,1d-3) + !mu = min(mu,1.e3*mc2_) tn(1) = 1/mu DO n=2,min(2,nre) tn(n) = tn(n-1)/mu @@ -359,7 +357,7 @@ sfd(1) = bet1-alp10*d0-alp14*sfd(4)-alp13*sfd(3)-alp12*sfd(2) !======================================================================= fc_old = fc - Te_old = Te + mu_old = mu Zeff_old = Zeff !--- sfdx(1:4) = sfd(1:4) @@ -405,7 +403,7 @@ q2 = q*q ! for the integrand, HSL_f gp1 = gam+1 ! .. !--- - CALL quanc8(HSL_f,zero,unit,atol,rtol,Integr,err,nfun,flag) + CALL quanc8(HSL_f,zero,one,atol,rtol,Integr,err,nfun,flag) !======================================================================= gam2 = gam*gam !--- diff --git a/src/itm_constants.f90 b/src/itm_constants.f90 deleted file mode 100644 index a617670..0000000 --- a/src/itm_constants.f90 +++ /dev/null @@ -1,32 +0,0 @@ -!> Module implementing the ITM physics constants -!> -!> Source: -!> based on SOLPS b2mod_constants.F -!> '09/12/07 xpb : source CODATA 2006 (http://www.nist.gov/)' -!> pulled from ets r100 -!> -!> \author David Coster -!> -!> \version "$Id: itm_constants.f90 37 2009-08-17 17:15:00Z coster $" - -module itm_constants - - use itm_types - - real (kind = R8), parameter :: itm_pi = 3.141592653589793238462643383280_R8 - real (kind = R8), parameter :: itm_c = 2.99792458e8_R8 ! speed of light, m/s - real (kind = R8), parameter :: itm_me = 9.10938215e-31_R8 ! electron mass, kg - real (kind = R8), parameter :: itm_mp = 1.672621637e-27_R8 ! proton mass, kg - real (kind = R8), parameter :: itm_md = 3.34358320e-27_R8 ! deuteron mass, kg - real (kind = R8), parameter :: itm_mt = 5.00735588e-27_R8 ! triton mass, kg - real (kind = R8), parameter :: itm_ma = 6.64465620e-27_R8 ! alpha mass, kg - real (kind = R8), parameter :: itm_amu = 1.660538782e-27_R8 ! amu, kg - real (kind = R8), parameter :: itm_ev = 1.602176487e-19_R8 - real (kind = R8), parameter :: itm_qe = itm_ev - real (kind = R8), parameter :: itm_mu0 = 4.0e-7_R8 * itm_pi - real (kind = R8), parameter :: itm_eps0 = 1.0_R8 / (itm_mu0 * itm_c * itm_c) - real (kind = R8), parameter :: itm_avogr = 6.02214179e23_R8 - real (kind = R8), parameter :: itm_KBolt = 1.3806504e-23_R8 - character (len=64), parameter :: itm_constants_version = '$Id: itm_constants.f90 37 2009-08-17 17:15:00Z coster $' - -end module itm_constants diff --git a/src/itm_types.f90 b/src/itm_types.f90 deleted file mode 100644 index 8a16580..0000000 --- a/src/itm_types.f90 +++ /dev/null @@ -1,50 +0,0 @@ -!> Module implementing the ITM basic types -!> -!> Source: -!> based on SOLPS b2mod_types.F -!> pulled from ets r100 and extended with input from C. Konz, T. Ribeiro & B. Scott -!> -!> \author David Coster -!> -!> \version "$Id: itm_types.f90 144 2010-10-07 09:26:24Z konz $" - -module itm_types - - INTEGER, PARAMETER :: ITM_I1 = SELECTED_INT_KIND (2) ! Integer*1 - INTEGER, PARAMETER :: ITM_I2 = SELECTED_INT_KIND (4) ! Integer*2 - INTEGER, PARAMETER :: ITM_I4 = SELECTED_INT_KIND (9) ! Integer*4 - INTEGER, PARAMETER :: ITM_I8 = SELECTED_INT_KIND (18) ! Integer*8 - INTEGER, PARAMETER :: R4 = SELECTED_REAL_KIND (6, 37) ! Real*4 - INTEGER, PARAMETER :: R8 = SELECTED_REAL_KIND (15, 300) ! Real*8 - - INTEGER, PARAMETER :: itm_int_invalid = -999999999 - REAL(R8), PARAMETER :: itm_r8_invalid = -9.0D40 - - interface itm_is_valid - module procedure itm_is_valid_int4, itm_is_valid_int8, itm_is_valid_real8 - end interface - -contains - - logical function itm_is_valid_int4(in_int) - implicit none - integer(ITM_I4) in_int - itm_is_valid_int4 = in_int .ne. itm_int_invalid - return - end function itm_is_valid_int4 - - logical function itm_is_valid_int8(in_int) - implicit none - integer(ITM_I8) in_int - itm_is_valid_int8 = in_int .ne. itm_int_invalid - return - end function itm_is_valid_int8 - - logical function itm_is_valid_real8(in_real) - implicit none - real(R8) in_real - itm_is_valid_real8 = abs(in_real - itm_r8_invalid) .gt. abs(itm_r8_invalid) * 1.0e-15_R8 - return - end function itm_is_valid_real8 - -end module itm_types diff --git a/src/reflections.f90 b/src/reflections.f90 index 6e04214..3de82da 100644 --- a/src/reflections.f90 +++ b/src/reflections.f90 @@ -1,21 +1,20 @@ module reflections + use const_and_precisions, only : wp_, comp_tiny, comp_eps, comp_huge, zero, one implicit none private - integer, parameter :: r8=selected_real_kind(15,300) - real(r8), parameter :: tinyr8=tiny(1._r8) public :: reflect,inters_linewall,inside public :: linecone_coord,interssegm_coord,interssegm contains subroutine reflect(ki,nsurf,ko) implicit none - real(r8), intent(in), dimension(3) :: ki - real(r8), intent(in), dimension(3) :: nsurf - real(r8), intent(out), dimension(3) :: ko - real(r8) :: twokn,norm2 + real(wp_), intent(in), dimension(3) :: ki + real(wp_), intent(in), dimension(3) :: nsurf + real(wp_), intent(out), dimension(3) :: ko + real(wp_) :: twokn,norm2 norm2 = dot_product(nsurf,nsurf) - if (norm2>0.0_r8) then - twokn = 2.0_r8*dot_product(ki,nsurf)/norm2 + if (norm2>zero) then + twokn = 2.0_wp_*dot_product(ki,nsurf)/norm2 ko=ki-twokn*nsurf else ko=ki @@ -24,19 +23,19 @@ end subroutine reflect subroutine inters_linewall(xv,kv,rw,zw,nw,sint,normw) implicit none - real(r8), intent(in), dimension(3) :: xv,kv + real(wp_), intent(in), dimension(3) :: xv,kv integer, intent(in) :: nw - real(r8), dimension(nw), intent(in) :: rw,zw - real(r8), intent(out) :: sint - real(r8), dimension(3), intent(out) :: normw + real(wp_), dimension(nw), intent(in) :: rw,zw + real(wp_), intent(out) :: sint + real(wp_), dimension(3), intent(out) :: normw integer :: i,j,ni,iint - real(r8), dimension(2) :: si,ti - real(r8) :: drw,dzw,xint,yint,rint,l,kxy - real(r8) :: tol - tol=sqrt(epsilon(1.0_r8)) - sint=huge(sint) + real(wp_), dimension(2) :: si,ti + real(wp_) :: drw,dzw,xint,yint,rint,l,kxy + real(wp_) :: tol + tol=sqrt(comp_eps) + sint=comp_huge iint=0 - normw=0.0_r8 + normw=zero do i=1,nw-1 !search intersections with i-th wall segment call linecone_coord(xv,kv,rw(i:i+1),zw(i:i+1),si,ti,ni) @@ -47,7 +46,7 @@ subroutine inters_linewall(xv,kv,rw,zw,nw,sint,normw) ti(1) = ti(2) end do do j=1,ni - if ((si(j)=0._r8 .and. ti(j)<=1._r8) then + if ((si(j)=zero .and. ti(j)<=one) then !check intersection is in r,z range and keep the closest sint = si(j) iint = i @@ -64,7 +63,7 @@ subroutine inters_linewall(xv,kv,rw,zw,nw,sint,normw) l = sqrt(drw**2+dzw**2) kxy = sqrt(kv(1)**2+kv(2)**2) normw(3) = -drw/l - if (rint>0.0_r8) then + if (rint>zero) then normw(1) = xint/rint*dzw/l normw(2) = yint/rint*dzw/l else @@ -72,17 +71,17 @@ subroutine inters_linewall(xv,kv,rw,zw,nw,sint,normw) normw(2) = kv(2)/kxy*dzw/l end if !reverse normal if k.n>0 - if (dot_product(normw,kv)>0.0_r8) normw=-normw + if (dot_product(normw,kv)>zero) normw=-normw end subroutine inters_linewall subroutine linecone_coord(xv,kv,rs,zs,s,t,n) implicit none - real(r8), intent(in), dimension(3) :: xv,kv - real(r8), intent(in), dimension(2) :: rs,zs - real(r8), dimension(2), intent(out) :: s,t + real(wp_), intent(in), dimension(3) :: xv,kv + real(wp_), intent(in), dimension(2) :: rs,zs + real(wp_), dimension(2), intent(out) :: s,t integer, intent(out) :: n - real(r8) :: x0,y0,z0,kx,ky,kz - real(r8) :: dr,dz,r,a,bhalf,c,delta,tvertex,zvertex,srmin,rmin,zrmin + real(wp_) :: x0,y0,z0,kx,ky,kz + real(wp_) :: dr,dz,r,a,bhalf,c,delta,tvertex,zvertex,srmin,rmin,zrmin x0=xv(1) y0=xv(2) z0=xv(3) @@ -93,9 +92,9 @@ subroutine linecone_coord(xv,kv,rs,zs,s,t,n) dz = zs(2)-zs(1) s = 0 t = 0 - if (abs(dz)=0._r8 .and. s<=1._r8 .and. & - t>=0._r8 .and. t<=1._r8) interssegm = .true. + if (ierr==0 .and. s>=zero .and. s<=one .and. & + t>=zero .and. t<=one) interssegm = .true. end function interssegm function inside(xc,yc,n,x,y) implicit none integer, intent(in) :: n - real(r8), dimension(n), intent(in) :: xc,yc - real(r8), intent(in) :: x,y + real(wp_), dimension(n), intent(in) :: xc,yc + real(wp_), intent(in) :: x,y logical :: inside integer, dimension(n) :: jint - real(r8), dimension(n) :: xint - real(r8), dimension(n+1) :: xclosed,yclosed + real(wp_), dimension(n) :: xint + real(wp_), dimension(n+1) :: xclosed,yclosed integer :: i,nj xclosed(1:n)=xc(1:n) yclosed(1:n)=yc(1:n) @@ -208,19 +207,19 @@ function intlin(x1,y1,x2,y2,x) result(y) !linear interpolation !must be x1 != x2 implicit none - real(r8),intent(in) :: x1,y1,x2,y2,x - real(r8) :: y - real(r8) :: a + real(wp_),intent(in) :: x1,y1,x2,y2,x + real(wp_) :: y + real(wp_) :: a a=(x2-x)/(x2-x1) - y=a*y1+(1._r8-a)*y2 + y=a*y1+(one-a)*y2 end function intlin subroutine locate_unord(a,n,x,j,m,nj) implicit none integer, intent(in) :: n,m integer, intent(out) :: nj - real(r8), dimension(n), intent(in) :: a - real(r8), intent(in) :: x + real(wp_), dimension(n), intent(in) :: a + real(wp_), intent(in) :: x integer, dimension(m), intent(inout) :: j integer :: i nj=0 @@ -240,8 +239,8 @@ function locate(a,n,x) result(j) !j=0 or j=n indicate that x is out of range (Numerical Recipes) implicit none integer, intent(in) :: n - real(r8), dimension(n), intent(in) :: a - real(r8), intent(in) :: x + real(wp_), dimension(n), intent(in) :: a + real(wp_), intent(in) :: x integer :: j integer :: jl,ju,jm logical :: incr @@ -262,8 +261,8 @@ end function locate subroutine order(p,q) !returns p,q in ascending order implicit none - real(r8), intent(inout) :: p,q - real(r8) :: temp + real(wp_), intent(inout) :: p,q + real(wp_) :: temp if (p>q) then temp=p p=q @@ -275,7 +274,7 @@ subroutine bubble(a,n) !bubble sorting of array a implicit none integer, intent(in) :: n - real(r8), dimension(n), intent(inout) :: a + real(wp_), dimension(n), intent(inout) :: a integer :: i, j do i=1,n do j=n,i+1,-1