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

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.

Browse Source
This commit is contained in:
Lorenzo Figini 2015-05-15 16:02:27 +00:00
parent 045c581865
commit e31f05e9a8
7 changed files with 228 additions and 279 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
Makefile
View File

@ -4,7 +4,7 @@ EXE=gray
# Objects list # Objects list
MAINOBJ=gray.o MAINOBJ=gray.o
OTHOBJ=dqagmv.o grayl.o reflections.o green_func_p.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 # Alternative search paths
vpath %.f90 src vpath %.f90 src
@ -23,17 +23,16 @@ $(EXE): $(MAINOBJ) $(OTHOBJ)
$(FC) $(FFLAGS) -o $@ $^ $(FC) $(FFLAGS) -o $@ $^
# Dependencies on modules # 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 green_func_p.o: const_and_precisions.o
const_and_precisions.o: itm_types.o itm_constants.o reflections.o: const_and_precisions.o
itm_constants.o: itm_types.o
# General object compilation command # General object compilation command
%.o: %.f90 %.o: %.f90
$(FC) $(FFLAGS) -c $^ $(FC) $(FFLAGS) -c $<
%.o: %.f %.o: %.f
$(FC) $(FFLAGS) -c $^ $(FC) $(FFLAGS) -c $<
gray.o:gray.f gray.o:gray.f
$(FC) -cpp $(DIRECTIVES) $(FFLAGS) -c $< $(FC) -cpp $(DIRECTIVES) $(FFLAGS) -c $<

107
src/const_and_precisions.f90
View File

@ -1,17 +1,21 @@
!########################################################################! !########################################################################!
MODULE const_and_precisions 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 IMPLICIT NONE
PUBLIC PUBLIC
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! common precisions ! common precisions
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! INTEGER, PARAMETER :: sp_ = 4 ! single precision ! INTEGER, PARAMETER :: i1 = SELECTED_INT_KIND (2) ! Integer*1
! INTEGER, PARAMETER :: dp_ = 8 ! double precision ! INTEGER, PARAMETER :: i2 = SELECTED_INT_KIND (4) ! Integer*2
! INTEGER, PARAMETER :: wp_ = dp_ ! work-precision 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 :: odep_ = dp_ ! ODE-solver precision
! INTEGER, PARAMETER :: xp_ = wp_ ! for ext. modules if necessary ! INTEGER, PARAMETER :: xp_ = wp_ ! for ext. modules if necessary
!------------------------------------------------------------------------ !------------------------------------------------------------------------
@ -26,31 +30,33 @@
!!======================================================================== !!========================================================================
! Arithmetic constants ! Arithmetic constants
!======================================================================== !========================================================================
integer, parameter :: izero = 0
REAL(wp_), PARAMETER :: zero = 0.0_wp_ REAL(wp_), PARAMETER :: zero = 0.0_wp_
REAL(wp_), PARAMETER :: unit = 1.0_wp_ REAL(wp_), PARAMETER :: one = 1.0_wp_
! REAL(wp_), PARAMETER :: pi = 3.141592653589793_wp_ real(wp_), parameter :: pi = 3.141592653589793_wp_ ! 3.141592653589793238462643383280
! REAL(wp_), PARAMETER :: sqrt_pi = 1.772453850905516_wp_ REAL(wp_), PARAMETER :: sqrt_pi = 1.772453850905516_wp_
! REAL(wp_), PARAMETER :: sqrt_2 = 1.414213562373095_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 :: ex(1:3) = (/one ,zero,zero/)
! REAL(wp_), PARAMETER :: ey(1:3) = (/zero,unit,zero/) ! REAL(wp_), PARAMETER :: ey(1:3) = (/zero,one ,zero/)
! REAL(wp_), PARAMETER :: ez(1:3) = (/zero,zero,unit/) ! REAL(wp_), PARAMETER :: ez(1:3) = (/zero,zero,one /)
!--- !---
! REAL(wp_), PARAMETER :: kron(3,3) = reshape((/unit,zero,zero, & ! REAL(wp_), PARAMETER :: kron(3,3) = reshape((/one ,zero,zero, &
! zero,unit,zero, & ! zero,one ,zero, &
! zero,zero,unit/),(/3,3/)) ! zero,zero,one /),(/3,3/))
! COMPLEX(wp_), PARAMETER :: im = (0.0_wp_,1.0_wp_) COMPLEX(wp_), PARAMETER :: im = (0.0_wp_,1.0_wp_)
! COMPLEX(wp_), PARAMETER :: czero = (0.0_wp_,0.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 :: cunit = (1.0_wp_,0.0_wp_)
! COMPLEX(wp_), PARAMETER :: ctwo = (2.0_wp_,0.0_wp_) ! COMPLEX(wp_), PARAMETER :: ctwo = (2.0_wp_,0.0_wp_)
!======================================================================== !========================================================================
! Computer constants ! Computer constants
!======================================================================== !========================================================================
REAL(wp_), PARAMETER :: comp_eps = EPSILON(unit) REAL(wp_), PARAMETER :: comp_eps = EPSILON(one)
! REAL(wp_), PARAMETER :: comp_eps2 = comp_eps**2 ! REAL(wp_), PARAMETER :: comp_eps2 = comp_eps**2
! REAL(wp_), PARAMETER :: comp_tiny = TINY(unit) REAL(wp_), PARAMETER :: comp_tiny = TINY(one)
! REAL(wp_), PARAMETER :: comp_huge = HUGE(unit) REAL(wp_), PARAMETER :: comp_huge = HUGE(one)
! REAL(wp_), PARAMETER :: comp_tinylog =-200 ! LOG10(comp_tiny) ! REAL(wp_), PARAMETER :: comp_tinylog =-200 ! LOG10(comp_tiny)
! REAL(wp_), PARAMETER :: comp_hugelog =+200 ! LOG10(comp_huge) ! REAL(wp_), PARAMETER :: comp_hugelog =+200 ! LOG10(comp_huge)
! REAL(wp_), PARAMETER :: comp_tiny1 = 1d+50*comp_tiny ! REAL(wp_), PARAMETER :: comp_tiny1 = 1d+50*comp_tiny
@ -60,23 +66,39 @@
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! Conventional constants ! 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_tiny = 1.0d-66
! REAL(wp_), PARAMETER :: output_huge = 1.0d+66 ! REAL(wp_), PARAMETER :: output_huge = 1.0d+66
!======================================================================== !========================================================================
! Physical constants (SI) ! Physical constants (SI)
!======================================================================== !========================================================================
! REAL(wp_), PARAMETER :: e_ = 1.602176487d-19 ! [C] real (wp_), parameter :: e_ = 1.602176487e-19_wp_ ! elementary charge, C
! REAL(wp_), PARAMETER :: me_ = 9.10938215d-31 ! [kg] real (wp_), parameter :: me_ = 9.10938215e-31_wp_ ! electron mass, kg
! REAL(wp_), PARAMETER :: mp_ = 1.672621637d-27 ! [kg] ! real (wp_), parameter :: mp_ = 1.672621637e-27_wp_ ! proton mass, kg
! REAL(wp_), PARAMETER :: rmpe_ = mp_/me_ ! real (wp_), parameter :: md_ = 3.34358320e-27_wp_ ! deuteron mass, kg
! REAL(wp_), PARAMETER :: c_ = 2.99792458d+08 ! [m/s] ! real (wp_), parameter :: mt_ = 5.00735588e-27_wp_ ! triton mass, kg
! REAL(wp_), PARAMETER :: eps0_ = 8.854187817d-12 ! [F/m] ! 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 ! 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_SI = me_*c_**2 ! [J]
REAL(wp_), PARAMETER :: mc2_ = mc2_SI/keV_ ! [keV] REAL(wp_), PARAMETER :: mc2_ = mc2_SI/keV_ ! [keV]
REAL(wp_), PARAMETER :: mu0inv = 1._wp_/mu0_ !
! REAL(wp_), PARAMETER :: mc_ = me_*c_ ! [kg*m/s] ! REAL(wp_), PARAMETER :: mc_ = me_*c_ ! [kg*m/s]
! ! f_ce = fce1_*B (B in Tesla): ! ! ! f_ce = fce1_*B (B in Tesla): !
! REAL(wp_), PARAMETER :: wce1_ = e_/me_ ! [rad/s] ! REAL(wp_), PARAMETER :: wce1_ = e_/me_ ! [rad/s]
@ -100,6 +122,33 @@
! REAL(wp_), PARAMETER :: Npar_min = 1.0d-3 ! 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 END MODULE const_and_precisions
!########################################################################! !########################################################################!

164
src/gray.f
View File

@ -109,8 +109,8 @@ c ray integration: end
subroutine after_gray_integration subroutine after_gray_integration
use const_and_precisions, only : zero
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
parameter(zero=0.0d0)
character*255 filenmeqq,filenmprf,filenmbm character*255 filenmeqq,filenmprf,filenmbm
c c
common/ss/st common/ss/st
@ -176,11 +176,12 @@ c
c c
c c
subroutine after_onestep(i,istop) subroutine after_onestep(i,istop)
use const_and_precisions, only : pi, cvdr=>degree
use reflections, only : inside use reflections, only : inside
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
complex*16 ext,eyt,extr,eytr,exin2,eyin2 complex*16 ext,eyt,extr,eytr,exin2,eyin2
parameter(jmx=31,kmx=36,nmx=8000,nbb=5000) 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 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 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 iop(jmx,kmx),iow(jmx,kmx),ihcd(jmx,kmx),istore(jmx,kmx)
@ -471,10 +472,10 @@ c
c c
c c
subroutine print_output(i,j,k) subroutine print_output(i,j,k)
use const_and_precisions, only : pi
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
parameter(ndim=6,jmx=31,kmx=36,nmx=8000) parameter(ndim=6,jmx=31,kmx=36,nmx=8000)
parameter(taucr=12.0d0) parameter(taucr=12.0d0)
parameter(pi=3.14159265358979d0)
dimension ywrk(ndim,jmx,kmx),ypwrk(ndim,jmx,kmx) dimension ywrk(ndim,jmx,kmx),ypwrk(ndim,jmx,kmx)
dimension psjki(jmx,kmx,nmx) dimension psjki(jmx,kmx,nmx)
dimension tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx) dimension tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx)
@ -648,14 +649,13 @@ c
c c
c c
subroutine read_data subroutine read_data
use const_and_precisions, only : qe=>ecgs_,me=>mecgs_,vc=>ccgs_,
* cvdr=>degree,pi
use green_func_p, only:Setup_SpitzFunc use green_func_p, only:Setup_SpitzFunc
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
real*8 me
character*8 wdat character*8 wdat
character*10 wtim character*10 wtim
character*255 filenmeqq,filenmprf,filenmbm 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) parameter(nmx=8000,nbb=5000)
real*8 rlim(nbb),zlim(nbb) real*8 rlim(nbb),zlim(nbb)
c c
@ -1031,8 +1031,8 @@ c
c c
c c
subroutine surf_anal subroutine surf_anal
use const_and_precisions, only : pi
implicit real*8(a-h,o-z) implicit real*8(a-h,o-z)
parameter(pi=3.14159265358979d0)
common/parban/b0,rr0m,zr0m,rpam common/parban/b0,rr0m,zr0m,rpam
common/parbres/bres common/parbres/bres
c c
@ -1183,9 +1183,9 @@ c
c c
c c
subroutine equidata subroutine equidata
use const_and_precisions, only : pi
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
parameter(nnw=501,nnh=501) parameter(nnw=501,nnh=501)
parameter(pi=3.14159265358979d0)
parameter(nbb=5000) parameter(nbb=5000)
c parameter(np=100) c parameter(np=100)
character*48 stringa character*48 stringa
@ -2270,8 +2270,8 @@ c
c c
c c
subroutine contours_psi(h,np,rcn,zcn,ipr) subroutine contours_psi(h,np,rcn,zcn,ipr)
use const_and_precisions, only : pi
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
parameter(pi=3.14159265358979d0)
parameter(mest=4,kspl=3) parameter(mest=4,kspl=3)
parameter(nnw=501,nnh=501) parameter(nnw=501,nnh=501)
parameter(nrest=nnw+4,nzest=nnh+4) parameter(nrest=nnw+4,nzest=nnh+4)
@ -2335,12 +2335,11 @@ c
c c
c c
subroutine flux_average subroutine flux_average
use const_and_precisions, only : zero,one,pi,ccj=>mu0inv
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
real*8 lam real*8 lam
parameter(nnintp=101,ncnt=100,ncntt=2*ncnt+1,nlam=41) 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(ksp=3,njest=nnintp+ksp+1,nlest=nlam+ksp+1)
parameter(lwrk=4*(nnintp+nlam)+11*(njest+nlest)+ parameter(lwrk=4*(nnintp+nlam)+11*(njest+nlest)+
. njest*nnintp+nlest+54) . njest*nnintp+nlest+54)
@ -3404,8 +3403,8 @@ c
c c
c c
subroutine plas_deriv(xx,yy,zz) subroutine plas_deriv(xx,yy,zz)
use const_and_precisions, only : pi
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
parameter(pi=3.14159265358979d0)
dimension bv(3),derxg(3),deryg(3),derbv(3,3),dbtot(3) dimension bv(3),derxg(3),deryg(3),derbv(3,3),dbtot(3)
dimension bvc(3),dbvcdc(3,3),dbvdc(3,3),dbv(3,3) dimension bvc(3),dbvcdc(3,3),dbvdc(3,3),dbv(3,3)
c c
@ -3786,8 +3785,8 @@ c
c c
subroutine tor_curr(rpsim,zpsim,ajphi) subroutine tor_curr(rpsim,zpsim,ajphi)
use const_and_precisions, only : pi,ccj=>mu0inv
implicit real*8 (a-h,o-z) 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 common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv
call equinum(rpsim,zpsim) call equinum(rpsim,zpsim)
bzz= dpsidr/rpsim bzz= dpsidr/rpsim
@ -3960,10 +3959,10 @@ c
c beam tracing initial conditions igrad=1 c beam tracing initial conditions igrad=1
c c
subroutine ic_gb subroutine ic_gb
use const_and_precisions, only : izero,zero,one,pi,cvdr=>degree
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
parameter(ndim=6,ndimm=3) parameter(ndim=6,ndimm=3)
parameter(jmx=31,kmx=36,zero=0.0d0,izero=0,one=1.0d0) parameter(jmx=31,kmx=36)
parameter(pi=3.14159265358979d0,cvdr=pi/180.0d0)
dimension ywrk0(ndim,jmx,kmx),ypwrk0(ndim,jmx,kmx) dimension ywrk0(ndim,jmx,kmx),ypwrk0(ndim,jmx,kmx)
dimension ytmp(ndim),yptmp(ndim) dimension ytmp(ndim),yptmp(ndim)
dimension xc0(ndimm,jmx,kmx),du10(ndimm,jmx,kmx) dimension xc0(ndimm,jmx,kmx),du10(ndimm,jmx,kmx)
@ -4275,12 +4274,11 @@ c
c ray tracing initial conditions igrad=0 c ray tracing initial conditions igrad=0
c c
subroutine ic_rt subroutine ic_rt
use const_and_precisions, only : izero,zero,one,pi,cvdr=>degree,
* ui=>im
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
complex*16 ui
parameter(ndim=6,ndimm=3) parameter(ndim=6,ndimm=3)
parameter(jmx=31,kmx=36,zero=0.0d0,izero=0,one=1.0d0) parameter(jmx=31,kmx=36)
parameter(pi=3.14159265358979d0,cvdr=pi/180.0d0)
parameter(ui=(0.0d0,1.0d0))
dimension ywrk0(ndim,jmx,kmx),ypwrk0(ndim,jmx,kmx) dimension ywrk0(ndim,jmx,kmx),ypwrk0(ndim,jmx,kmx)
dimension ytmp(ndim),yptmp(ndim) dimension ytmp(ndim),yptmp(ndim)
dimension xc0(ndimm,jmx,kmx),du10(ndimm,jmx,kmx) dimension xc0(ndimm,jmx,kmx),du10(ndimm,jmx,kmx)
@ -4472,10 +4470,10 @@ c
subroutine ic_rt2 subroutine ic_rt2
use const_and_precisions, only : izero,zero,one,pi,cvdr=>degree
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
parameter(ndim=6,ndimm=3) parameter(ndim=6,ndimm=3)
parameter(jmx=31,kmx=36,zero=0.0d0,izero=0,one=1.0d0) parameter(jmx=31,kmx=36)
parameter(pi=3.14159265358979d0,cvdr=pi/180.0d0)
dimension ywrk0(ndim,jmx,kmx),ypwrk0(ndim,jmx,kmx) dimension ywrk0(ndim,jmx,kmx),ypwrk0(ndim,jmx,kmx)
dimension ytmp(ndim),yptmp(ndim) dimension ytmp(ndim),yptmp(ndim)
dimension yyrfl(jmx,kmx,ndim) dimension yyrfl(jmx,kmx,ndim)
@ -4685,9 +4683,9 @@ c
c c
c c
subroutine pabs_curr(i,j,k) subroutine pabs_curr(i,j,k)
use const_and_precisions, only : pi
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
parameter(jmx=31,kmx=36,nmx=8000) parameter(jmx=31,kmx=36,nmx=8000)
parameter(pi=3.14159265358979d0)
dimension psjki(jmx,kmx,nmx) dimension psjki(jmx,kmx,nmx)
dimension tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx) dimension tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx)
dimension pdjki(jmx,kmx,nmx),ppabs(jmx,kmx,nmx) dimension pdjki(jmx,kmx,nmx),ppabs(jmx,kmx,nmx)
@ -4764,12 +4762,11 @@ c
c c
c c
subroutine ecrh_cd subroutine ecrh_cd
use const_and_precisions, only : eqsi=>e_,mesi=>me_,vcsi=>c_,
* mc2=>mc2_
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
real*8 mc2,mesi
parameter(taucr=12.0d0,xxcr=16.0d0,eps=1.d-8) parameter(taucr=12.0d0,xxcr=16.0d0,eps=1.d-8)
parameter(qesi=1.602176487d-19,mesi=9.10938215d-31) complex*16 ex,ey,ez
parameter(vcsi=2.99792458d+8)
parameter(mc2 = mesi*vcsi*vcsi/qesi*1d-3)
c c
common/nharm/nharm,nhf common/nharm/nharm,nhf
common/warm/iwarm,ilarm common/warm/iwarm,ilarm
@ -4781,16 +4778,22 @@ c
common/parwv/ak0,akinv,fhz common/parwv/ak0,akinv,fhz
c c
common/nprw/anprre,anprim common/nprw/anprre,anprim
common/epolar/ex,ey,ez
c c
common/absor/alpha,effjcd,akim,tau common/absor/alpha,effjcd,akim,tau
c c
common/ierr/ierr common/ierr/ierr
common/iokh/iokhawa common/iokh/iokhawa
c c
common/psival/psinv common/psival/psinv
common/tete/tekev common/tete/tekev
common/dens/dens,ddens
common/amut/amu common/amut/amu
common/zz/Zeff common/zz/Zeff
common/btot/btot
common/bmxmn/bmax,bmin
common/fc/fc
c c
c absorption computation c absorption computation
c c
@ -4852,8 +4855,12 @@ c ratiovgr=2.0d0*anpr/derdnm*vgm
c c
ithn=1 ithn=1
if(lrm.gt.nharm) ithn=2 if(lrm.gt.nharm) ithn=2
if(ieccd.gt.0) call eccd(yg,anpl,anprre,amu,zeff, if(ieccd.gt.0) then
* ieccd,nharm,nhf,ithn,effjcd,iokhawa,ierr) 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 return
999 format(12(1x,e12.5)) 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 up to third order in Larmor radius for hermitian part
c c
subroutine diel_tens_fr(e330,epsl,lrm) 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) implicit real*8(a-h,o-z)
complex*16 ui
complex*16 e330,epsl(3,3,lrm) complex*16 e330,epsl(3,3,lrm)
complex*16 ca11,ca12,ca22,ca13,ca23,ca33 complex*16 ca11,ca12,ca22,ca13,ca23,ca33
complex*16 cq0p,cq0m,cq1p,cq1m,cq2p 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) dimension rr(-lrm:lrm,0:2,0:lrm),ri(lrm,0:2,lrm)
c c
common/xgxg/xg common/xgxg/xg
@ -5455,10 +5460,9 @@ c
c c
c c
subroutine antihermitian(ri,lrm) subroutine antihermitian(ri,lrm)
use const_and_precisions, only : rpi=>sqrt_pi
implicit none implicit none
integer lmx,nmx,lrm,n,k,m,mm integer lmx,nmx,lrm,n,k,m,mm
real*8 rpi
parameter(rpi=1.77245385090552d0)
parameter(lmx=20,nmx=lmx+2) parameter(lmx=20,nmx=lmx+2)
real*8 fsbi(nmx) real*8 fsbi(nmx)
real*8 ri(lrm,0:2,lrm) 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 Krivenki and Orefice, JPP 30,125 (1983)
c c
subroutine diel_tens_wr(ce330,cepsl,lrm) subroutine diel_tens_wr(ce330,cepsl,lrm)
use const_and_precisions, only : cui=>im
implicit real*8 (a-b,d-h,o-z) implicit real*8 (a-b,d-h,o-z)
implicit complex*16 (c) implicit complex*16 (c)
dimension cefp(0:lrm,0:2),cefm(0:lrm,0:2),cepsl(3,3,lrm) dimension cefp(0:lrm,0:2),cefm(0:lrm,0:2),cepsl(3,3,lrm)
parameter(cui=(0.0d0,1.0d0))
c c
common/xgxg/xg common/xgxg/xg
common/ygyg/yg common/ygyg/yg
@ -5651,10 +5655,10 @@ c
c c
c c
subroutine fsup(cefp,cefm,lrm) subroutine fsup(cefp,cefm,lrm)
use const_and_precisions, only : cui=>im
implicit real*8 (a-b,d-h,o-z) implicit real*8 (a-b,d-h,o-z)
implicit complex*16 (c) implicit complex*16 (c)
parameter(apsicr=0.7d0) parameter(apsicr=0.7d0)
parameter(cui=(0.0d0,1.0d0))
dimension cefp(0:lrm,0:2),cefm(0:lrm,0:2) dimension cefp(0:lrm,0:2),cefm(0:lrm,0:2)
c c
common/ygyg/yg common/ygyg/yg
@ -5796,42 +5800,32 @@ c
subroutine eccd(yg,anpl,anprre,amu,zeff, subroutine eccd(yg,anpl,anprre,amu,zeff,rbn,rbx,fc,ex,ey,ez,
* ieccd,nhmn,nhmx,ithn,effjcd,iokhawa,ierr) * dens,psinv,ieccd,nhmn,nhmx,ithn,effjcd,iokhawa,ierr)
use green_func_p, only: SpitzFuncCoeff,pi 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 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 yg,anpl,anpr,amu,anprre,anprim
real*8 mc2,anucc,canucc,ddens real*8 mc2m2,anucc,canucc,ddens
real*8 qesi,mesi,vcsi,qe,me,vc,ceff,Zeff,psinv real*8 ceff,Zeff,psinv
real*8 rbn,rbx,btot,bmin,bmax,alams,fp0s,pa,fc real*8 rbn,rbx,alams,fp0s,pa,fc
real*8 fjch,fjncl,fjch0,fconic real*8 fjch,fjncl,fjch0,fconic
real*8 cst2,eccdpar(5) real*8 cst2,eccdpar(5)
complex*16 ex,ey,ez complex*16 ex,ey,ez
integer ieccd,nhmn,nhmx,nhn,ithn,ierr,iokhawa integer ieccd,nhmn,nhmx,nhn,ithn,ierr,iokhawa
external fjch,fjncl,fjch0 external fjch,fjncl,fjch0
parameter(qesi=1.602176487d-19,mesi=9.10938215d-31) parameter(mc2m2=1.0d0/mc2**2)
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(canucc=2.0d13*pi*qe**4/(me**2*vc**3)) parameter(canucc=2.0d13*pi*qe**4/(me**2*vc**3))
parameter(ceff=qesi/(mesi*vcsi)) 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 anum=0.0d0
denom=0.0d0 denom=0.0d0
effjcd=0.0d0 effjcd=0.0d0
tekev=mc2/amu coullog=48.0d0-log(1.0d7*dens*mc2m2*amu**2)
coullog=48.0d0-log(1.0d7*dens/tekev**2)
anucc=canucc*dens*coullog anucc=canucc*dens*coullog
c nhmx=nhmn c nhmx=nhmn
@ -5848,11 +5842,9 @@ c cst2=1.0d0-B/B_max
c alams=sqrt(1-B_min/B_max) c alams=sqrt(1-B_min/B_max)
c Zeff < 31 !!! c Zeff < 31 !!!
c fp0s= P_a (alams) c fp0s= P_a (alams)
rbn=btot/bmin
rbx=btot/bmax
cst2=1.0d0-rbx cst2=1.0d0-rbx
if(cst2.lt.1d-6) cst2=0.0d0 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 pa=sqrt(32.0d0/(Zeff+1.0d0)-1.0d0)/2.0d0
fp0s=fconic(alams,pa,0) fp0s=fconic(alams,pa,0)
eccdpar(2)=rbn eccdpar(2)=rbn
@ -5879,18 +5871,15 @@ c fp0s= P_a (alams)
c neoclassical model: c neoclassical model:
c ft=1-fc trapped particle fraction c ft=1-fc trapped particle fraction
c rzfc=(1+Zeff)/fc c rzfc=(1+Zeff)/fc
rbn=btot/bmin
rbx=btot/bmax
cst2=1.0d0-rbx cst2=1.0d0-rbx
if(cst2.lt.1d-6) cst2=0.0d0 if(cst2.lt.1d-6) cst2=0.0d0
call SpitzFuncCoeff(Tekev,Zeff,fc) call SpitzFuncCoeff(amu,Zeff,fc)
eccdpar(2) = tekev eccdpar(2) = fc
eccdpar(3) = fc eccdpar(3) = rbx
eccdpar(4) = rbx eccdpar(4) = psinv
eccdpar(5) = psinv
do nhn=nhmn,nhmx do nhn=nhmn,nhmx
call curr_int(yg,anpl,anprre,amu,ex,ey,ez,nhn,ithn,cst2, 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 anum=anum+resj
denom=denom+resp denom=denom+resp
end do end do
@ -6071,9 +6060,9 @@ c extrapar(13) = uplp
c extrapar(14) = uplm c extrapar(14) = uplm
c extrapar(15) = ygn c extrapar(15) = ygn
c c
use const_and_precisions, only : ui=>im
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
complex*16 ex,ey,ez,ui,emxy,epxy complex*16 ex,ey,ez,emxy,epxy
parameter(ui=(0.0d0,1.0d0))
dimension extrapar(npar) dimension extrapar(npar)
yg=extrapar(1) yg=extrapar(1)
@ -6293,10 +6282,9 @@ c extrapar(14) = uplm
c extrapar(15) = ygn c extrapar(15) = ygn
c c
c extrapar(16) = zeff c extrapar(16) = zeff
c extrapar(17) = tekev c extrapar(17) = fc
c extrapar(18) = fc c extrapar(18) = hb
c extrapar(19) = hb c extrapar(29) = psin
c extrapar(20) = psin
c c
use green_func_p, only: GenSpitzFunc use green_func_p, only: GenSpitzFunc
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
@ -6306,10 +6294,9 @@ c
amu=extrapar(3) amu=extrapar(3)
ygn=extrapar(15) ygn=extrapar(15)
zeff=extrapar(16) zeff=extrapar(16)
tekev=extrapar(17) fc=extrapar(17)
fc=extrapar(18) hb=extrapar(18)
hb=extrapar(19) psin=extrapar(19)
psin=extrapar(20)
gam=anpl*upl+ygn gam=anpl*upl+ygn
u2=gam*gam-1.0d0 u2=gam*gam-1.0d0
@ -6317,7 +6304,7 @@ c
upr2=u2-upl*upl upr2=u2-upl*upl
bth=sqrt(2.0d0/amu) bth=sqrt(2.0d0/amu)
uth=u/bth 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) fk=fk*(4.0d0/amu**2)
dfk=dfk*(2.0d0/amu)*bth dfk=dfk*(2.0d0/amu)*bth
@ -6450,10 +6437,10 @@ c
c c
c c
subroutine pec(pabs,currt) subroutine pec(pabs,currt)
use const_and_precisions, only : pi,one
implicit real*8(a-h,o-z) implicit real*8(a-h,o-z)
parameter(nndmx=5001,jmx=31,kmx=36,nmx=8000) parameter(nndmx=5001,jmx=31,kmx=36,nmx=8000)
parameter(rtbc=1.0d0) parameter(rtbc=one)
parameter(pi=3.14159265358979d0)
dimension psjki(jmx,kmx,nmx),iiv(jmx,kmx) dimension psjki(jmx,kmx,nmx),iiv(jmx,kmx)
dimension ppabs(jmx,kmx,nmx),ccci(jmx,kmx,nmx) dimension ppabs(jmx,kmx,nmx),ccci(jmx,kmx,nmx)
dimension pdjki(jmx,kmx,nmx),currj(jmx,kmx,nmx) dimension pdjki(jmx,kmx,nmx),currj(jmx,kmx,nmx)
@ -6810,8 +6797,8 @@ c
c c
c c
subroutine profwidth(nd,xx,yy,rhotmx,rhopmx,ypk,drhot,drhop) subroutine profwidth(nd,xx,yy,rhotmx,rhopmx,ypk,drhot,drhop)
use const_and_precisions, only : emn1
implicit real*8(a-h,o-z) implicit real*8(a-h,o-z)
parameter(emn1=0.367879441171442d0)
dimension xx(nd),yy(nd) dimension xx(nd),yy(nd)
common/jmxmn/rhotp,rhotm,ypkp,ypkm common/jmxmn/rhotp,rhotm,ypkp,ypkm
common/ipec/ipec,nnd common/ipec/ipec,nnd
@ -6973,15 +6960,15 @@ c
end end
subroutine pol_limit(ext,eyt) subroutine pol_limit(ext,eyt)
use const_and_precisions, only : ui=>im,pi
implicit none implicit none
real*8 bv(3),anv(3) real*8 bv(3),anv(3)
real*8 anx,any,anz,anpl,anpr,yg,xe2om,ye2om,xe2xm,ye2xm 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 an2,an,anxy,sngam,csgam,csg2,sng2,ffo,ffx,ffo2,ffx2
real*8 deno,denx,anpl2,dnl,del0 real*8 deno,denx,anpl2,dnl,del0
real*8 pi,beta0,alpha0,gam real*8 beta0,alpha0,gam
real*8 sox real*8 sox
complex*16 ui,exom,eyom,exxm,eyxm,ext,eyt complex*16 exom,eyom,exxm,eyxm,ext,eyt
parameter(ui=(0.0d0,1.0d0),pi=3.14159265358979d0)
c c
common/anv/anv common/anv/anv
common/nplr/anpl,anpr common/nplr/anpl,anpr
@ -7043,11 +7030,10 @@ c
end subroutine stokes end subroutine stokes
subroutine polellipse(qq,uu,vv,psipol,chipol) subroutine polellipse(qq,uu,vv,psipol,chipol)
use const_and_precisions, only : pi
implicit none implicit none
real*8 qq,uu,vv,psipol,chipol real*8 qq,uu,vv,psipol,chipol
c real*8 llm,aa,bb,ell c real*8 llm,aa,bb,ell
real*8 pi
parameter(pi=3.14159265358979d0)
c llm=sqrt(qq**2+uu**2) c llm=sqrt(qq**2+uu**2)
c aa=sqrt((1+llm)/2.0d0) c aa=sqrt((1+llm)/2.0d0)
c bb=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, subroutine wall_refl(xv,anv,ext,eyt,xvrfl,anvrfl,extr,eytr,walln,
. irfl) . irfl)
use reflections, only : inters_linewall,inside use reflections, only : inters_linewall,inside
use const_and_precisions, only : ui=>im,pi
implicit none implicit none
integer*4 irfl integer*4 irfl
real*8 anv(3),anv0(3),xv(3),xvrfl(3) real*8 anv(3),anv0(3),xv(3),xvrfl(3)
real*8 walln(3),anvrfl(3),vv1(3),vv2(3),vv3(3) real*8 walln(3),anvrfl(3),vv1(3),vv2(3),vv3(3)
real*8 pi,smax,rrm,zzm real*8 smax,rrm,zzm
complex*16 ui,extr,eytr,eztr,ext,eyt complex*16 extr,eytr,eztr,ext,eyt
complex*16 evin(3),evrfl(3) complex*16 evin(3),evrfl(3)
parameter(ui=(0.0d0,1.0d0),pi=3.14159265358979d0)
integer nbb,nlim integer nbb,nlim
parameter(nbb=5000) parameter(nbb=5000)
real*8 rlim(nbb),zlim(nbb) real*8 rlim(nbb),zlim(nbb)

32
src/green_func_p.f90
View File

@ -107,8 +107,7 @@
END SUBROUTINE Setup_SpitzFunc 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 ! Author: N.B.Marushchenko
! June 2005: as start point the subroutine of Ugo Gasparino (198?) ! 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 ! dKdu = dK/du, i.e. its derivative over normalized momentum
!======================================================================= !=======================================================================
IMPLICIT NONE 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_), INTENT(out) :: K,dKdu
REAL(wp_) :: mu,gam1,gam2,gam3,w,dwdu REAL(wp_) :: gam1,gam2,gam3,w,dwdu
!======================================================================= !=======================================================================
K = 0 K = 0
dKdu = 0 dKdu = 0
IF (u < comp_eps) RETURN IF (u < comp_eps) RETURN
!--- !---
mu = mc2_/max(Te,1d-3)
SELECT CASE(adj_appr(2)) SELECT CASE(adj_appr(2))
CASE('m') !--------------- momentum conservation ------------------! CASE('m') !--------------- momentum conservation ------------------!
gam1 = gam ! gam1 = gam !
@ -210,7 +208,7 @@
!####################################################################### !#######################################################################
!####################################################################### !#######################################################################
SUBROUTINE SpitzFuncCoeff(Te,Zeff,fc) SUBROUTINE SpitzFuncCoeff(mu,Zeff,fc)
!======================================================================= !=======================================================================
! Calculates the matrix coefficients required for the subroutine ! Calculates the matrix coefficients required for the subroutine
! "GenSpitzFunc", where the Spitzer function is defined through the ! "GenSpitzFunc", where the Spitzer function is defined through the
@ -228,7 +226,7 @@
! INPUT VARIABLES: ! INPUT VARIABLES:
! rho = sqrt(SS) with SS - flux-surface label (norm. magn. flux) ! rho = sqrt(SS) with SS - flux-surface label (norm. magn. flux)
! ne - density, 1/m^3 ! ne - density, 1/m^3
! Te - temperature, keV ! mu - mc2/Te
! Zeff - effective charge ! Zeff - effective charge
! fc - fraction of circulating particles ! fc - fraction of circulating particles
! !
@ -237,9 +235,9 @@
! "Spitzer"-function (the same as in the Hirshman paper) ! "Spitzer"-function (the same as in the Hirshman paper)
!======================================================================= !=======================================================================
IMPLICIT NONE IMPLICIT NONE
REAL(wp_), INTENT(in) :: Te,Zeff,fc REAL(wp_), INTENT(in) :: mu,Zeff,fc
INTEGER :: n,i,j 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_) :: m(0:4,0:4),g(0:4)
REAL(wp_) :: om(0:4,0:4) REAL(wp_) :: om(0:4,0:4)
REAL(wp_) :: gam11,gam21,gam31,gam41,gam01, & REAL(wp_) :: gam11,gam21,gam31,gam41,gam01, &
@ -250,17 +248,17 @@
alp23,alp24,alp20, & alp23,alp24,alp20, &
alp34,alp30,alp40 alp34,alp30,alp40
REAL(wp_) :: bet0,bet1,bet2,bet3,bet4,d0 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 :: 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 rel = mu < mc2_
newTe = abs(Te -Te_old ) > delta*Te newmu = abs(mu -mu_old ) > delta*mu
newZ = abs(Zeff-Zeff_old) > delta*Zeff newZ = abs(Zeff-Zeff_old) > delta*Zeff
newfc = abs(fc -fc_old ) > delta*fc newfc = abs(fc -fc_old ) > delta*fc
SELECT CASE(adj_appr(1)) SELECT CASE(adj_appr(1))
CASE ('l','c') CASE ('l','c')
renew = (newTe .and. rel) .OR. newZ .OR. newfc renew = (newmu .and. rel) .OR. newZ .OR. newfc
END SELECT END SELECT
!--- !---
IF (.not.renew) THEN IF (.not.renew) THEN
@ -271,7 +269,7 @@
tn(:) = 0 tn(:) = 0
IF (adj_appr(4) == 'r') THEN IF (adj_appr(4) == 'r') THEN
IF (nre > 0) THEN IF (nre > 0) THEN
mu = mc2_/max(Te,1d-3) !mu = min(mu,1.e3*mc2_)
tn(1) = 1/mu tn(1) = 1/mu
DO n=2,min(2,nre) DO n=2,min(2,nre)
tn(n) = tn(n-1)/mu tn(n) = tn(n-1)/mu
@ -359,7 +357,7 @@
sfd(1) = bet1-alp10*d0-alp14*sfd(4)-alp13*sfd(3)-alp12*sfd(2) sfd(1) = bet1-alp10*d0-alp14*sfd(4)-alp13*sfd(3)-alp12*sfd(2)
!======================================================================= !=======================================================================
fc_old = fc fc_old = fc
Te_old = Te mu_old = mu
Zeff_old = Zeff Zeff_old = Zeff
!--- !---
sfdx(1:4) = sfd(1:4) sfdx(1:4) = sfd(1:4)
@ -405,7 +403,7 @@
q2 = q*q ! for the integrand, HSL_f q2 = q*q ! for the integrand, HSL_f
gp1 = gam+1 ! .. 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 gam2 = gam*gam
!--- !---

32
src/itm_constants.f90
View File

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

50
src/itm_types.f90
View File

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

111
src/reflections.f90
View File

@ -1,21 +1,20 @@
module reflections module reflections
use const_and_precisions, only : wp_, comp_tiny, comp_eps, comp_huge, zero, one
implicit none implicit none
private private
integer, parameter :: r8=selected_real_kind(15,300)
real(r8), parameter :: tinyr8=tiny(1._r8)
public :: reflect,inters_linewall,inside public :: reflect,inters_linewall,inside
public :: linecone_coord,interssegm_coord,interssegm public :: linecone_coord,interssegm_coord,interssegm
contains contains
subroutine reflect(ki,nsurf,ko) subroutine reflect(ki,nsurf,ko)
implicit none implicit none
real(r8), intent(in), dimension(3) :: ki real(wp_), intent(in), dimension(3) :: ki
real(r8), intent(in), dimension(3) :: nsurf real(wp_), intent(in), dimension(3) :: nsurf
real(r8), intent(out), dimension(3) :: ko real(wp_), intent(out), dimension(3) :: ko
real(r8) :: twokn,norm2 real(wp_) :: twokn,norm2
norm2 = dot_product(nsurf,nsurf) norm2 = dot_product(nsurf,nsurf)
if (norm2>0.0_r8) then if (norm2>zero) then
twokn = 2.0_r8*dot_product(ki,nsurf)/norm2 twokn = 2.0_wp_*dot_product(ki,nsurf)/norm2
ko=ki-twokn*nsurf ko=ki-twokn*nsurf
else else
ko=ki ko=ki
@ -24,19 +23,19 @@ end subroutine reflect
subroutine inters_linewall(xv,kv,rw,zw,nw,sint,normw) subroutine inters_linewall(xv,kv,rw,zw,nw,sint,normw)
implicit none implicit none
real(r8), intent(in), dimension(3) :: xv,kv real(wp_), intent(in), dimension(3) :: xv,kv
integer, intent(in) :: nw integer, intent(in) :: nw
real(r8), dimension(nw), intent(in) :: rw,zw real(wp_), dimension(nw), intent(in) :: rw,zw
real(r8), intent(out) :: sint real(wp_), intent(out) :: sint
real(r8), dimension(3), intent(out) :: normw real(wp_), dimension(3), intent(out) :: normw
integer :: i,j,ni,iint integer :: i,j,ni,iint
real(r8), dimension(2) :: si,ti real(wp_), dimension(2) :: si,ti
real(r8) :: drw,dzw,xint,yint,rint,l,kxy real(wp_) :: drw,dzw,xint,yint,rint,l,kxy
real(r8) :: tol real(wp_) :: tol
tol=sqrt(epsilon(1.0_r8)) tol=sqrt(comp_eps)
sint=huge(sint) sint=comp_huge
iint=0 iint=0
normw=0.0_r8 normw=zero
do i=1,nw-1 do i=1,nw-1
!search intersections with i-th wall segment !search intersections with i-th wall segment
call linecone_coord(xv,kv,rw(i:i+1),zw(i:i+1),si,ti,ni) 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) ti(1) = ti(2)
end do end do
do j=1,ni do j=1,ni
if ((si(j)<sint .or. iint==0) .and. ti(j)>=0._r8 .and. ti(j)<=1._r8) then if ((si(j)<sint .or. iint==0) .and. ti(j)>=zero .and. ti(j)<=one) then
!check intersection is in r,z range and keep the closest !check intersection is in r,z range and keep the closest
sint = si(j) sint = si(j)
iint = i iint = i
@ -64,7 +63,7 @@ subroutine inters_linewall(xv,kv,rw,zw,nw,sint,normw)
l = sqrt(drw**2+dzw**2) l = sqrt(drw**2+dzw**2)
kxy = sqrt(kv(1)**2+kv(2)**2) kxy = sqrt(kv(1)**2+kv(2)**2)
normw(3) = -drw/l normw(3) = -drw/l
if (rint>0.0_r8) then if (rint>zero) then
normw(1) = xint/rint*dzw/l normw(1) = xint/rint*dzw/l
normw(2) = yint/rint*dzw/l normw(2) = yint/rint*dzw/l
else else
@ -72,17 +71,17 @@ subroutine inters_linewall(xv,kv,rw,zw,nw,sint,normw)
normw(2) = kv(2)/kxy*dzw/l normw(2) = kv(2)/kxy*dzw/l
end if end if
!reverse normal if k.n>0 !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 end subroutine inters_linewall
subroutine linecone_coord(xv,kv,rs,zs,s,t,n) subroutine linecone_coord(xv,kv,rs,zs,s,t,n)
implicit none implicit none
real(r8), intent(in), dimension(3) :: xv,kv real(wp_), intent(in), dimension(3) :: xv,kv
real(r8), intent(in), dimension(2) :: rs,zs real(wp_), intent(in), dimension(2) :: rs,zs
real(r8), dimension(2), intent(out) :: s,t real(wp_), dimension(2), intent(out) :: s,t
integer, intent(out) :: n integer, intent(out) :: n
real(r8) :: x0,y0,z0,kx,ky,kz real(wp_) :: x0,y0,z0,kx,ky,kz
real(r8) :: dr,dz,r,a,bhalf,c,delta,tvertex,zvertex,srmin,rmin,zrmin real(wp_) :: dr,dz,r,a,bhalf,c,delta,tvertex,zvertex,srmin,rmin,zrmin
x0=xv(1) x0=xv(1)
y0=xv(2) y0=xv(2)
z0=xv(3) z0=xv(3)
@ -93,9 +92,9 @@ subroutine linecone_coord(xv,kv,rs,zs,s,t,n)
dz = zs(2)-zs(1) dz = zs(2)-zs(1)
s = 0 s = 0
t = 0 t = 0
if (abs(dz)<tinyr8) then if (abs(dz)<comp_tiny) then
!surface in horizontal plane !surface in horizontal plane
if (abs(kz)<tinyr8 .or. abs(dr)<tinyr8) then if (abs(kz)<comp_tiny .or. abs(dr)<comp_tiny) then
n = 0 n = 0
else else
s(1) = (zs(1)-z0)/kz s(1) = (zs(1)-z0)/kz
@ -107,9 +106,9 @@ subroutine linecone_coord(xv,kv,rs,zs,s,t,n)
a = (kx**2+ky**2) - (dr/dz*kz)**2 a = (kx**2+ky**2) - (dr/dz*kz)**2
bhalf = -dr/dz*kz*rs(1) + (kx*x0 + ky*y0) - (dr/dz)**2*kz*(z0-zs(1)) bhalf = -dr/dz*kz*rs(1) + (kx*x0 + ky*y0) - (dr/dz)**2*kz*(z0-zs(1))
c = (x0**2+y0**2) - (rs(1) + dr/dz*(z0-zs(1)))**2 c = (x0**2+y0**2) - (rs(1) + dr/dz*(z0-zs(1)))**2
if (abs(a)<tinyr8) then if (abs(a)<comp_tiny) then
!line parallel to cone generator !line parallel to cone generator
if (abs(dr)<tinyr8) then if (abs(dr)<comp_tiny) then
!cylinder and vertical line !cylinder and vertical line
n = 0 n = 0
else else
@ -118,14 +117,14 @@ subroutine linecone_coord(xv,kv,rs,zs,s,t,n)
srmin = -(kx*x0 + ky*y0)/(kx**2+ky**2) srmin = -(kx*x0 + ky*y0)/(kx**2+ky**2)
rmin = sqrt((x0+srmin*kx)**2+(y0+srmin*ky)**2) rmin = sqrt((x0+srmin*kx)**2+(y0+srmin*ky)**2)
zrmin = z0 + srmin*kz zrmin = z0 + srmin*kz
if (rmin<tinyr8 .and. abs(zrmin-zvertex)<tinyr8) then if (rmin<comp_tiny .and. abs(zrmin-zvertex)<comp_tiny) then
!line passing by cone vertex !line passing by cone vertex
!s(1) = srmin !s(1) = srmin
!t(1) = tvertex !t(1) = tvertex
!n = 1 !n = 1
n = 0 n = 0
else else
s(1) = -0.5_r8*c/bhalf s(1) = -0.5_wp_*c/bhalf
t(1) = (kz*s(1)+(z0-zs(1)))/dz t(1) = (kz*s(1)+(z0-zs(1)))/dz
n = 1 n = 1
end if end if
@ -147,18 +146,18 @@ end subroutine linecone_coord
subroutine interssegm_coord(xa,ya,xb,yb,s,t,ierr) subroutine interssegm_coord(xa,ya,xb,yb,s,t,ierr)
implicit none implicit none
real(r8), dimension(2), intent(in) :: xa,ya,xb,yb real(wp_), dimension(2), intent(in) :: xa,ya,xb,yb
real(r8), intent(out) :: s,t real(wp_), intent(out) :: s,t
integer, intent(out) :: ierr integer, intent(out) :: ierr
real(r8) :: crossprod,dxa,dya,dxb,dyb real(wp_) :: crossprod,dxa,dya,dxb,dyb
dxa = xa(2)-xa(1) dxa = xa(2)-xa(1)
dya = ya(2)-ya(1) dya = ya(2)-ya(1)
dxb = xb(2)-xb(1) dxb = xb(2)-xb(1)
dyb = yb(2)-yb(1) dyb = yb(2)-yb(1)
crossprod = dxb*dya - dxa*dyb crossprod = dxb*dya - dxa*dyb
if (abs(crossprod)<tiny(crossprod)) then if (abs(crossprod)<tiny(crossprod)) then
s = 0.0_r8 s = zero
t = 0.0_r8 t = zero
ierr = 1 ierr = 1
else else
s = (dyb*(xa(1)-xb(1)) - dxb*(ya(1)-yb(1)))/crossprod s = (dyb*(xa(1)-xb(1)) - dxb*(ya(1)-yb(1)))/crossprod
@ -169,25 +168,25 @@ end subroutine interssegm_coord
function interssegm(xa,ya,xb,yb) function interssegm(xa,ya,xb,yb)
implicit none implicit none
real(r8), dimension(2), intent(in) :: xa,ya,xb,yb real(wp_), dimension(2), intent(in) :: xa,ya,xb,yb
logical :: interssegm logical :: interssegm
real(r8) :: s,t real(wp_) :: s,t
integer :: ierr integer :: ierr
interssegm = .false. interssegm = .false.
call interssegm_coord(xa,ya,xb,yb,s,t,ierr) call interssegm_coord(xa,ya,xb,yb,s,t,ierr)
if (ierr==0 .and. s>=0._r8 .and. s<=1._r8 .and. & if (ierr==0 .and. s>=zero .and. s<=one .and. &
t>=0._r8 .and. t<=1._r8) interssegm = .true. t>=zero .and. t<=one) interssegm = .true.
end function interssegm end function interssegm
function inside(xc,yc,n,x,y) function inside(xc,yc,n,x,y)
implicit none implicit none
integer, intent(in) :: n integer, intent(in) :: n
real(r8), dimension(n), intent(in) :: xc,yc real(wp_), dimension(n), intent(in) :: xc,yc
real(r8), intent(in) :: x,y real(wp_), intent(in) :: x,y
logical :: inside logical :: inside
integer, dimension(n) :: jint integer, dimension(n) :: jint
real(r8), dimension(n) :: xint real(wp_), dimension(n) :: xint
real(r8), dimension(n+1) :: xclosed,yclosed real(wp_), dimension(n+1) :: xclosed,yclosed
integer :: i,nj integer :: i,nj
xclosed(1:n)=xc(1:n) xclosed(1:n)=xc(1:n)
yclosed(1:n)=yc(1:n) yclosed(1:n)=yc(1:n)
@ -208,19 +207,19 @@ function intlin(x1,y1,x2,y2,x) result(y)
!linear interpolation !linear interpolation
!must be x1 != x2 !must be x1 != x2
implicit none implicit none
real(r8),intent(in) :: x1,y1,x2,y2,x real(wp_),intent(in) :: x1,y1,x2,y2,x
real(r8) :: y real(wp_) :: y
real(r8) :: a real(wp_) :: a
a=(x2-x)/(x2-x1) a=(x2-x)/(x2-x1)
y=a*y1+(1._r8-a)*y2 y=a*y1+(one-a)*y2
end function intlin end function intlin
subroutine locate_unord(a,n,x,j,m,nj) subroutine locate_unord(a,n,x,j,m,nj)
implicit none implicit none
integer, intent(in) :: n,m integer, intent(in) :: n,m
integer, intent(out) :: nj integer, intent(out) :: nj
real(r8), dimension(n), intent(in) :: a real(wp_), dimension(n), intent(in) :: a
real(r8), intent(in) :: x real(wp_), intent(in) :: x
integer, dimension(m), intent(inout) :: j integer, dimension(m), intent(inout) :: j
integer :: i integer :: i
nj=0 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) !j=0 or j=n indicate that x is out of range (Numerical Recipes)
implicit none implicit none
integer, intent(in) :: n integer, intent(in) :: n
real(r8), dimension(n), intent(in) :: a real(wp_), dimension(n), intent(in) :: a
real(r8), intent(in) :: x real(wp_), intent(in) :: x
integer :: j integer :: j
integer :: jl,ju,jm integer :: jl,ju,jm
logical :: incr logical :: incr
@ -262,8 +261,8 @@ end function locate
subroutine order(p,q) subroutine order(p,q)
!returns p,q in ascending order !returns p,q in ascending order
implicit none implicit none
real(r8), intent(inout) :: p,q real(wp_), intent(inout) :: p,q
real(r8) :: temp real(wp_) :: temp
if (p>q) then if (p>q) then
temp=p temp=p
p=q p=q
@ -275,7 +274,7 @@ subroutine bubble(a,n)
!bubble sorting of array a !bubble sorting of array a
implicit none implicit none
integer, intent(in) :: n integer, intent(in) :: n
real(r8), dimension(n), intent(inout) :: a real(wp_), dimension(n), intent(inout) :: a
integer :: i, j integer :: i, j
do i=1,n do i=1,n
do j=n,i+1,-1 do j=n,i+1,-1