!2017.6.11 !To test, modify the temporal.f90 so that the calculated HLLD flux is 0. module common_module use sktk_fortran use setup_template,only:vax,bax,sax,eq_state,eforce,eheat,Ohm,iflux_damp implicit none character(20),parameter :: outdir="/nwork/sakaue/test1/" integer n,rkstep real(DP) dx,dy,dz integer,parameter :: xbox=512 integer,parameter :: xmar(2)=(/ 10,10 /)*int(1._DP-0._DP**(xbox-1)) integer,parameter :: ybox=1 integer,parameter :: ymar(2)=(/ 10,10 /)*int(1._DP-0._DP**(ybox-1)) integer,parameter :: zbox=1 integer,parameter :: zmar(2)=(/ 0,0 /) integer,parameter :: tlen=100 integer,parameter :: xlen=xbox+xmar(1)+xmar(2) integer,parameter :: ylen=ybox+ymar(1)+ymar(2) integer,parameter :: zlen=zbox+zmar(1)+zmar(2) integer,parameter :: xr(2)=(/xmar(1)+1,xmar(1)+xbox/) integer,parameter :: yr(2)=(/ymar(1)+1,ymar(1)+ybox/) integer,parameter :: zr(2)=(/zmar(1)+1,zmar(1)+zbox/) integer,parameter :: udim(3)=1-int(0._DP**(/xbox-1,ybox-1,zbox-1/)) integer :: udir(sum(udim)) real(DP),parameter :: gamma=5._DP/3._DP,gamma1i=1.0_DP/(gamma-1.0_DP) real(DP),parameter :: pressure_inf=p_yocto !Scale of physical quantities real(DP),parameter :: scV=.1_DP**0,scL=.1_DP**0,scRho=1._DP ![cm/s], [cm], [g/cm3] real(DP),parameter :: scTmp=scV*scV/1._DP ![K] ! real(DP),parameter :: scTmp=scV*scV/c_Rg ![K] real(DP),parameter :: scPre=scRho*scV*scV real(DP),parameter :: scMag=sqrt(4._DP*pi*scRho)*scV real(DP),dimension(xlen,ylen,zlen,8)::U integer,parameter :: taxlen=10 real(DP) :: tax(taxlen+1)=0._DP !-------------------------------------------------------------! ! Coordinate System ------------------------------------------! real(DP) :: xax(xlen),yax(ylen),zax(zlen) real(DP),allocatable,dimension(:,:,:) :: cax1,cax2,cax3 real(DP),allocatable,dimension(:,:,:) ::hx,hy,hz !-------------------------------------------------------------! real(DP),parameter :: cflf=.5_DP integer :: HLLDerror=0,INTERPerror=0,pressure_blurred=-1 integer,parameter :: heat_conduction=1 integer,parameter :: miyoshi_kusano_2011=0 integer,parameter :: powell_et_al_1999=0 integer,parameter :: initial_flux_cancel=0 real(DP),allocatable,dimension(:,:,:,:) :: initial_flux integer,parameter :: torder=3 !-------------------------------! ! torder ! ! 1...Euler ! ! 2...TVD Runge-Kutta 2nd ! ! 3...TVD Runge-Kutta 3rd ! ! 4...TVD Runge-Kutta 4th ! !-------------------------------! real(DP),parameter :: pow=0._DP real(DP),parameter :: Qin=.1_DP**2 real(DP),parameter :: kappa0=.1_DP**2*scRho*scV**3*scL*scTmp**(-pow-1) real(DP),parameter :: tau_c=Qin*Qin/(kappa0/gamma1i)*100000._DP integer :: parrish_stone_2005=0 contains function outputtime(tout) implicit none real(DP) outputtime,dt_out,tout integer,parameter :: nskip=2 dt_out=tax(taxlen+1)*dble(nskip) dt_out=tau_c*.000005_DP ! dt_out=5._DP*.1_DP**2 if (parrish_stone_2005.eq.1) dt_out=2._DP outputtime=tout+dt_out*ge0(tax(mod(n,taxlen)+1)-tout) end function outputtime function heat_conductivity(temperature) implicit none real(DP) :: heat_conductivity real(DP) :: temperature,temp real(DP),parameter :: scHC=scRho*scV**3*scL/scTmp real(DP),parameter :: kappa_sp=p_micro !CGS ! for test -------------------------------------------------------------- ! heat_conductivity=10._DP**6*scTmp/(scRho*scV**3*scL) ! heat_conductivity=10._DP**3*scTmp/(scRho*scV**3*scL)*temp ! heat_conductivity=10._DP*temp**2.5_DP*scHc heat_conductivity=kappa0*temperature**pow ! heat_conductivity=.01_DP ! for test--------------------------------------------------------------- end function heat_conductivity function U2eth(U8) implicit none real(DP) :: U8(8),U2eth U2eth=U8(Sax(2))-.5_DP*(sum(U8(Vax)*U8(Vax))/U8(Sax(1))& +sum(U8(Bax)*U8(Bax))) end function U2eth function U2p(U8,gamma1i) implicit none real(DP),intent(IN) :: gamma1i,U8(8) real(DP) :: U2p U2p=U2eth(U8)/gamma1i end function U2p function CT_preprocess(U0,xi,yj,zk) implicit none real(DP),dimension(:,:,:,:),intent(IN)::U0 integer,intent(IN) :: xi,yj,zk real(DP),dimension(8) :: CT_preprocess integer s,kdlt(3) CT_preprocess=U0(xi,yj,zk,:) do s=1,3 kdlt=0**abs(s-(/1,2,3/)) CT_preprocess(Bax(s))=.5_DP*(CT_preprocess(Bax(s))& +U0(max(xi-kdlt(1),1),max(yj-kdlt(2),1),max(zk-kdlt(3),1),Bax(s))) end do end function CT_preprocess function cnsrvtv2prmtv(conservative,gamma1i) implicit none real(DP),intent(IN) :: gamma1i,conservative(8) real(DP) :: cnsrvtv2prmtv(8),Eth,Ekin,Emag Ekin=.5_DP*sum(conservative(Vax)*conservative(Vax))/conservative(Sax(1)) Emag=.5_DP*sum(conservative(Bax)*conservative(Bax)) Eth=conservative(Sax(2))-Ekin-Emag cnsrvtv2prmtv(Sax(1))=conservative(Sax(1)) cnsrvtv2prmtv(Sax(2))=max(Eth/gamma1i,pressure_inf) cnsrvtv2prmtv(Vax)=conservative(Vax)/conservative(Sax(1)) cnsrvtv2prmtv(Bax)=conservative(Bax) end function cnsrvtv2prmtv function prmtv2cnsrvtv(primitive,gamma1i) implicit none real(DP),intent(IN) :: gamma1i,primitive(8) real(DP) :: prmtv2cnsrvtv(8),Eth,Ekin,Emag Eth=max(primitive(Sax(2)),pressure_inf)*gamma1i Ekin=.5_DP*sum(primitive(Vax)*primitive(Vax))*primitive(Sax(1)) Emag=.5_DP*sum(primitive(Bax)*primitive(Bax)) prmtv2cnsrvtv(Sax(1))=primitive(Sax(1)) prmtv2cnsrvtv(Sax(2))=Eth+Ekin+Emag prmtv2cnsrvtv(Vax)=primitive(Vax)*primitive(Sax(1)) prmtv2cnsrvtv(Bax)=primitive(Bax) end function prmtv2cnsrvtv function scale_factors(xi,yj,zk) implicit none real(DP) :: scale_factors(3) integer :: xi,yj,zk,shh(3) shh=shape(hx) scale_factors& =(/hx(min(max(xi,1),shh(1)),min(max(yj,1),shh(2)),& min(max(zk,1),shh(3))),& hy(min(max(xi,1),shh(1)),min(max(yj,1),shh(2)),& min(max(zk,1),shh(3))),& hz(min(max(xi,1),shh(1)),min(max(yj,1),shh(2)),& min(max(zk,1),shh(3)))/) end function scale_factors subroutine cell_center(xi,yj,zk,pc,hc) implicit none integer,intent(IN) :: xi,yj,zk real(DP),intent(OUT) :: pc(3),hc(3) hc=scale_factors(xi,yj,zk) pc=(/xax(min(max(xi,1),xlen)),yax(min(max(yj,1),ylen)),& zax(min(max(zk,1),zlen))/) end subroutine cell_center subroutine cell_face(xi,yj,zk,axis,dir,pf,hf,sf) implicit none integer,intent(IN) :: xi,yj,zk,axis,dir real(DP),intent(OUT) :: pf(3),hf(3),sf integer :: kdlt(3),kdlt2(3),kdlt3(3),axis2,axis3,cnt real(DP) :: hc(3),hR(3),hL(3),pc(3),pR(3),pL(3),sfR,sfL axis2=mod(axis,3)+1 ; axis3=mod(axis2,3)+1 kdlt=0**abs(axis-(/1,2,3/)) call cell_center(xi,yj,zk,pc,hc) if (udim(axis).eq.0) then pf=pc; hf=hc sf=2._DP*hc(axis2)*hc(axis3)*hf(axis2)*hf(axis3)& /(hc(axis2)*hc(axis3)+hf(axis2)*hf(axis3)) return end if call cell_center(xi+kdlt(1),yj+kdlt(2),zk+kdlt(3),pR,hR) call cell_center(xi-kdlt(1),yj-kdlt(2),zk-kdlt(3),pL,hL) do cnt=1,3 pf(cnt)=Lag_interpol2(-1._DP,0._DP,1._DP,pL(cnt),pc(cnt),pR(cnt),& .5_DP*dir) end do do cnt=1,3 hf(cnt)=Lag_interpol2(pL(axis),pc(axis),pR(axis),hL(cnt),hc(cnt),& hR(cnt),pf(axis)) end do sfR=2._DP*hc(axis2)*hc(axis3)*hR(axis2)*hR(axis3)& /(hc(axis2)*hc(axis3)+hR(axis2)*hR(axis3)) sfL=2._DP*hc(axis2)*hc(axis3)*hL(axis2)*hL(axis3)& /(hc(axis2)*hc(axis3)+hL(axis2)*hL(axis3)) sf=sfR*eq0(dir-1._DP)+sfL*eq0(dir+1._DP) end subroutine cell_face subroutine cell_edge(xi,yj,zk,axis1,axis2,pe,he) implicit none integer,intent(IN) :: xi,yj,zk,axis1,axis2 real(DP),intent(OUT) :: pe(3),he(3) integer :: kdlt1(3),kdlt2(3),kdlt3(3),axis3 real(DP) :: hc(3),hp1(3),hp2(3),hp3(3) kdlt1=0**abs(axis1-(/1,2,3/)) kdlt2=0**abs(axis2-(/1,2,3/)) kdlt3=1-kdlt1-kdlt2 pe(1)=(xax(min(max(xi+1-kdlt3(1),1),xlen))+xax(min(max(xi,1),xlen)))*.5_DP pe(2)=(yax(min(max(yj+1-kdlt3(2),1),ylen))+yax(min(max(yj,1),ylen)))*.5_DP pe(3)=(zax(min(max(zk+1-kdlt3(3),1),zlen))+zax(min(max(zk,1),zlen)))*.5_DP hc=scale_factors(xi,yj,zk) hp1=scale_factors(xi+kdlt1(1),yj+kdlt1(2),zk+kdlt1(3)) hp2=scale_factors(xi+kdlt2(1),yj+kdlt2(2),zk+kdlt2(3)) hp3=scale_factors(xi+1-kdlt3(1),yj+1-kdlt3(2),zk+1-kdlt3(3)) he(1)=(hc(1)+hp1(1)+hp2(1)+hp3(1))*.25_DP he(2)=(hc(2)+hp1(2)+hp2(2)+hp3(2))*.25_DP he(3)=(hc(3)+hp1(3)+hp2(3)+hp3(3))*.25_DP end subroutine cell_edge function divBf(vec,xi,yj,zk) ! vec is defined at the cell face implicit none real(DP) :: divBf real(DP),dimension(xlen,ylen,zlen,3),intent(IN):: vec integer,intent(IN) :: xi,yj,zk integer :: s,axis,kdlt(3) real(DP) :: pfR(3),pfL(3),hfR(3),hfL(3),hc(3) real(DP) :: sfR,sfL,dBxdx,dBydy,dBzdz hc=scale_factors(xi,yj,zk) divBf=0._DP do s=1,sum(udim) axis=udir(s); kdlt=0**abs(axis-(/1,2,3/)) call cell_face(xi,yj,zk,axis, 1,pfR,hfR,sfR) call cell_face(xi,yj,zk,axis,-1,pfL,hfL,sfL) divBf=divBf+(sfR*vec(xi,yj,zk,axis)& -sfL*vec(max(xi-kdlt(1),1),max(yj-kdlt(2),1),& max(zk-kdlt(3),1),axis))& /(pfR(axis)-pfL(axis)) end do divBf=divBf/hc(1)/hc(2)/hc(3) end function divBf function rotAe(vec,xi,yj,zk,axis) !vec is defined at the cell edge (right-up) !rotAe is defined at the cell face implicit none real(DP) :: rotAe real(DP),dimension(xlen,ylen,zlen,3),intent(IN):: vec integer,intent(IN) :: xi,yj,zk,axis real(DP) :: pfR(3),hfR(3),sfR,denomi real(DP),dimension(3) :: peR,heR,peL,heL integer :: axis2,axis3,kdlt2(3),kdlt3(3) axis2=mod(axis,3)+1; axis3=mod(axis2,3)+1 kdlt2=0**abs(axis2-(/1,2,3/)); kdlt3=0**abs(axis3-(/1,2,3/)) call cell_face(xi,yj,zk,axis,1,pfR,hfR,sfR) call cell_edge(xi,yj,zk,axis,axis2,peR,heR) call cell_edge(xi-kdlt2(1),yj-kdlt2(2),zk-kdlt2(3),axis,axis2,peL,heL) denomi=peR(axis2)-peL(axis2); denomi=ne0(denomi)/(denomi+eq0(denomi)) rotAe=(heR(axis3)*vec(xi,yj,zk,axis3)& -heL(axis3)*vec(xi-kdlt2(1),yj-kdlt2(2),zk-kdlt2(3),axis3))*denomi call cell_edge(xi,yj,zk,axis,axis3,peR,heR) call cell_edge(xi-kdlt3(1),yj-kdlt3(2),zk-kdlt3(3),axis,axis3,peL,heL) denomi=peR(axis3)-peL(axis3); denomi=ne0(denomi)/(denomi+eq0(denomi)) rotAe=rotAe-(heR(axis2)*vec(xi,yj,zk,axis2)& -heL(axis2)*vec(xi-kdlt3(1),yj-kdlt3(2),zk-kdlt3(3),axis2))*denomi rotAe=rotAe/sfR end function rotAe function rotAf(vec,xi,yj,zk,axis) !vec is defined at the cell edge (right) !rotAf is defined at the cell center implicit none real(DP) :: rotAf real(DP),dimension(xlen,ylen,zlen,3),intent(IN):: vec integer,intent(IN) :: xi,yj,zk,axis real(DP) :: pc(3),hc(3),pfR(3),hfR(3),sfR,pfL(3),hfL(3),denomi integer :: axis2,axis3,kdlt2(3),kdlt3(3) axis2=mod(axis,3)+1; axis3=mod(axis2,3)+1 kdlt2=0**abs(axis2-(/1,2,3/)); kdlt3=0**abs(axis3-(/1,2,3/)) call cell_center(xi,yj,zk,pc,hc) call cell_face(xi,yj,zk,axis2, 1,pfR,hfR,sfR) call cell_face(xi,yj,zk,axis2,-1,pfL,hfL,sfR) denomi=pfR(axis2)-pfL(axis2); denomi=ne0(denomi)/(denomi+eq0(denomi)) rotAf=(hfR(axis3)*vec(xi,yj,zk,axis3)& -hfL(axis3)*vec(xi-kdlt2(1),yj-kdlt2(2),zk-kdlt2(3),axis3))*denomi call cell_face(xi,yj,zk,axis3, 1,pfR,hfR,sfR) call cell_face(xi,yj,zk,axis3,-1,pfL,hfL,sfR) denomi=pfR(axis3)-pfL(axis3); denomi=ne0(denomi)/(denomi+eq0(denomi)) rotAf=rotAf-(hfR(axis2)*vec(xi,yj,zk,axis2)& -hfL(axis2)*vec(xi-kdlt3(1),yj-kdlt3(2),zk-kdlt3(3),axis2))*denomi rotAf=rotAf/(hc(axis2)*hc(axis3)) end function rotAf end module common_module module initial_module use common_module implicit none contains subroutine ic implicit none integer :: i,j,k real(DP) :: radius,inring real(DP),dimension(xlen,ylen,zlen)::rho,pre,ene real(DP),dimension(xlen,ylen,zlen,3)::vel,mag,veca j=1 do i=1,3 if (udim(i).eq.1) then udir(j)=i; j=j+1 end if end do allocate(hx(1,1,1)) !allocate(hx(xlen,ylen,zlen)) allocate(hy(1,1,1)) !allocate(hy(xlen,ylen,zlen)) allocate(hz(1,1,1)) !allocate(hz(xlen,ylen,zlen)) allocate(cax1(1,1,1)) !allocate(cax1(xlen,ylen,zlen)) allocate(cax2(1,1,1)) !allocate(cax2(xlen,ylen,zlen)) allocate(cax3(1,1,1)) !allocate(cax3(xlen,ylen,zlen)) if (initial_flux_cancel .ge. 1) then allocate(initial_flux(xlen,ylen,zlen,8)) else allocate(initial_flux(1,1,1,1)) end if !axis coordinate +++++++++++++++++++++++++++++++++++++++++++++ do i=1,xlen; xax(i)=2._DP/dble(xbox)*dble(i-xr(1))-1._DP end do dx=minval(abs(xax(2:xlen)-xax(1:xlen-1))) do j=1,ylen; yax(j)=2._DP/dble(ybox)*dble(j-yr(1))-1._DP enddo dy=minval(abs(yax(2:ylen)-yax(1:ylen-1))) do k=1,zlen; zax(k)=1._DP/dble(zbox)*dble(k-zr(1)) enddo dz=minval(abs(zax(2:zlen)-zax(1:zlen-1))) !------------------------------------------------------------- !Orthogonal Curvilinear Coordinate System ++++++++++++++++++++ hx=1._DP hy=1._DP hz=1._DP cax1=1._DP cax2=1._DP cax3=1._DP !------------------------------------------------------------- !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! !+ +! !+ Define the initial condition with primitive variables +! !+ +! !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! vel=0._DP if (parrish_stone_2005 .eq. 1 .and. sum(udim).eq.2) then do i=1,xlen; do j=1,ylen; do k=1,zlen veca(i,j,k,:)=(/0._DP,0._DP,& sqrt((xax(i)+dx*.5_DP)**2+(yax(j)+dy*.5_DP)**2)/) end do; end do; end do do i=1,xlen; do j=1,ylen mag(i,j,:,1)=rotAe(veca,i,j,1,1) mag(i,j,:,2)=rotAe(veca,i,j,1,2) mag(i,j,:,3)=0._DP end do; end do rho=1._DP do i=1,xlen; do j=1,ylen radius=sqrt(xax(i)**2+yax(j)**2) inring=le0(radius-.7)*ge0(radius-.5)& *ge0(acos(-xax(i)/radius)-11._DP*pi/12._DP) pre(i,j,:)=eq_state(-1._DP,rho(i,j,1),12._DP)*inring+& eq_state(-1._DP,rho(i,j,1),10._DP)*(1._DP-inring) end do; end do else tax(1)=.1_DP**4 !pow=0._DP tax(1)=.1_DP**9 !pow=1._DP tax(1)=p_micro*tau_c rho=1._DP; mag=0._DP do i=1,xlen; do j=1,ylen; do k=1,zlen ! pre(i,j,k)=1._DP+10000._DP*exp(-xax(i)**2/(1._DP*dx)**2) pre(i,j,k)=pressure_analytic(tax(1),xax(i),rho(1,1,1),kappa0,pow,Qin) pre(i,j,k)=max(pre(i,j,k),pressure_inf) end do; end do; end do end if do i=1,xlen; do j=1,ylen; do k=1,zlen U(i,j,k,Bax)=mag(i,j,k,:) U(i,j,k,Vax)=vel(i,j,k,:)*rho(i,j,k) U(i,j,k,Sax(1))=rho(i,j,k) end do; end do; end do do i=1,xlen; do j=1,ylen; do k=1,zlen U(i,j,k,Sax(2))=pre(i,j,k)*gamma1i& +sum(U(i,j,k,Vax)**2)/U(i,j,k,Sax(1))/2._DP& +(((U(i,j,k,Bax(1))+U(max(i-1,1),j,k,Bax(1)))/2._DP)**2& +((U(i,j,k,Bax(2))+U(i,max(j-1,1),k,Bax(2)))/2._DP)**2& +((U(i,j,k,Bax(3))+U(i,j,max(k-1,1),Bax(3)))/2._DP)**2)/2._DP end do; end do; end do contains function pressure_analytic(t0,x0,rho0,kappa0,pow,Q0) implicit none real(DP) :: t0,x0,rho0,kappa0,pow,Q0 real(DP) :: a0,xi,xi0,f_xi,gmmfnc(2) real(DP) :: pressure_analytic a0=kappa0*(gamma-1.d0) if (pow.gt.0._DP) then xi=x0/(Q0**pow*a0*t0)**(1._DP/(pow+2._DP)) xi0=(pow+2._DP)**(pow+1._DP)*2._DP**(1._DP-pow)& *(gamma_fnc(.5_DP+1._DP/pow)/gamma_fnc(1._DP/pow))**pow& /pow/pi**(pow*.5_DP) xi0=xi0**(1._DP/(pow+2._DP)) f_xi=(.5_DP*pow*(xi0*xi0-xi*xi)/(pow+2._DP))**(1._DP/pow) else f_xi=exp(-x0*x0/4._DP/a0/t0)/2._DP/sqrt(pi) end if pressure_analytic=(Q0*Q0/a0/t0)**(1._DP/(pow+2._DP))*f_xi end function pressure_analytic function gamma_fnc(re) ! sikinote (http://slpr.sakura.ne.jp/qp) ! author : sikino ! date : 2017/01/09 ! 2017/01/10 optimize the number of a(M) ! 2017/10/27 fix for positive integer implicit none real(DP),intent(in) :: re complex(DP)::z complex(DP)::cgamma real(DP) :: gamma_fnc integer::i,n,M double precision,parameter::e1=0.3678794411714423216d0 double precision,parameter::ln2pi2=0.91893853320467274d0 double precision,parameter::sq2pi=2.5066282746310005d0 double precision::a(1:30),AB(1:14),dz,iz complex(kind(0d0))::q,s,w,r,q1,q2 z=cmplx(re,0._DP) dz=dble(z) iz=imag(z) if(dz.eq.nint(dz).and.iz.eq.0d0)then if(dz.gt.0d0)then s=dcmplx(1d0,0d0) n=nint(dz) do i=2,n-1 s=s*i enddo else s=1d+300 endif cgamma=s else q=z if(dz.lt.0d0)q=1d0-z if(abs(iz).lt.1.45d0)then ! For x=arb, |y|\lt 1.45 n=int(dble(q)) s=1d0 do i=1,n s=s*(q-i) enddo w=q-dble(n) a(1) = 1d0 a(2) = 0.57721566490153286d0 a(3) =-0.65587807152025388d0 a(4) =-0.42002635034095236d-1 a(5) = 0.16653861138229149d0 a(6) =-0.42197734555544337d-1 a(7) =-0.96219715278769736d-2 a(8) = 0.72189432466630995d-2 a(9) =-0.11651675918590651d-2 a(10)=-0.21524167411495097d-3 a(11)= 0.12805028238811619d-3 a(12)=-0.20134854780788239d-4 a(13)=-0.12504934821426707d-5 a(14)= 0.11330272319816959d-5 a(15)=-0.20563384169776071d-6 a(16)= 0.61160951044814158d-8 a(17)= 0.50020076444692229d-8 a(18)=-0.11812745704870201d-8 a(19)= 0.10434267116911005d-9 a(20)= 0.77822634399050713d-11 a(21)=-0.36968056186422057d-11 a(22)= 0.51003702874544760d-12 a(23)=-0.20583260535665068d-13 a(24)=-0.53481225394230180d-14 a(25)= 0.12267786282382608d-14 a(26)=-0.11812593016974588d-15 a(27)= 0.11866922547516003d-17 a(28)= 0.14123806553180318d-17 a(29)=-0.22987456844353702d-18 a(30)= 0.17144063219273374d-19 M=int(11.3*abs(w)+13) if(M.gt.30)M=30 r=a(M) do i=M-1,1,-1 r=r*w+a(i) enddo cgamma=s/(r*w) else ! For x=arb, |y|\gt 1.45 s=1d0 if(abs(q).lt.9d0)then do i=0,7 s=s*(q+i) enddo s=1d0/s q=q+8d0 endif AB(1) = 0.83333333333333333d-1 AB(2) =-0.27777777777777778d-2 AB(3) = 0.79365079365079365d-3 AB(4) =-0.59523809523809524d-3 AB(5) = 0.84175084175084175d-3 AB(6) =-0.19175269175269175d-2 AB(7) = 0.64102564102564103d-2 AB(8) =-0.29550653594771242d-1 q1=1d0/q q2=q1*q1 r=AB(8) do i=7,1,-1 r=r*q2+AB(i) enddo cgamma=s*exp((q-0.5d0)*log(q)-q+ln2pi2+r*q1) endif if(dz.lt.0d0)cgamma=pi/(cgamma*sin(pi*z)) endif gamma_fnc=real(cgamma) return end function gamma_fnc end subroutine ic end module initial_module