module sktk_fortran implicit none integer,parameter ::SP=kind(1.0) integer,parameter ::DP=selected_real_kind(2*precision(1.0_SP)) real(DP),parameter :: pi=acos(-1._DP) real(DP),parameter :: e1=exp(1._DP) real(DP),parameter :: p_yotta=1000000000000000000000000._DP !10^24 real(DP),parameter :: p_zetta=1000000000000000000000._DP !10^21 real(DP),parameter :: p_exa =1000000000000000000._DP !10^18 real(DP),parameter :: p_peta =1000000000000000._DP !10^15 real(DP),parameter :: p_tera =1000000000000._DP !10^12 real(DP),parameter :: p_giga =1000000000._DP !10^9 real(DP),parameter :: p_mega =1000000._DP !10^6 real(DP),parameter :: p_kilo =1000._DP !10^3 real(DP),parameter :: p_mili =.001_DP !10^-3 real(DP),parameter :: p_micro=.000001_DP !10^-6 real(DP),parameter :: p_nano =.000000001_DP !10^-9 real(DP),parameter :: p_pico =.000000000001_DP !10^-12 real(DP),parameter :: p_femto=.000000000000001_DP !10^-15 real(DP),parameter :: p_atto =.000000000000000001_DP !10^-18 real(DP),parameter :: p_zepto=.000000000000000000001_DP !10^-21 real(DP),parameter :: p_yocto=.000000000000000000000001_DP !10^-24 real(DP),parameter :: c_mp=1.672621637_DP*(.1_DP**24) real(DP),parameter :: c_mn=1.674927211_DP*(.1_DP**24) real(DP),parameter :: c_me=9.10938215_DP*(.1_DP**27) real(DP),parameter :: c_Rg=8.314472_DP*(10._DP**7) real(DP),parameter :: c_kB=1.3806504_DP*(.1_DP**16) real(DP),parameter :: c_ele=4.803204272_DP*(.1_DP**10) real(DP),parameter :: c_c0=29979245800._DP real(DP),parameter :: c_G=6.67259_DP*.1_DP**8 real(DP),parameter :: c_gsun=27478._DP real(DP),parameter :: c_rsun=695700._DP*10._DP**5 contains !------------------------------------------------------------------------ ! List of contained subroutine ! report, time_print, loopcounter, rk_butcher ! ! List of contained function ! v_mv, m_mm, eq0, gt0, ge0, lt0, le0, dfdx_c2, dfdx_c4 !------------------------------------------------------------------------ subroutine report(a,b,c) implicit none integer a,b,c,d real :: percent percent=c/100. if (mod(a/(b+0.),percent) .lt. mod((a-1)/(b+0.),percent)) then d=int(int(a*10000./(b*c))*percent) if (d .ne. 0) then print'(i10,"% have done.")',d end if end if end subroutine report subroutine time_print implicit none integer :: time(8) character(10) :: dtz(3) call date_and_time(dtz(1),dtz(2),dtz(3),time) print'(i4,".",i2,".",i2," ",i2.2,":",i2.2,":",i2.2,".",i3.3)',& time(1),time(2),time(3),time(5),time(6),time(7),time(8) end subroutine time_print subroutine stopwatch(itime,ftime) implicit none integer :: itime,t_rate,t_max,hr,min real(DP) :: sec,diff integer,optional :: ftime if (present(ftime)) then call system_clock(ftime,t_rate,t_max) if (ftime .lt. itime) then diff=((t_max-itime)+ftime+1)/dble(t_rate) else diff=(ftime-itime)/dble(t_rate) endif hr=int(diff/3600._DP) min=int((diff-3600._DP*hr)/60._DP) sec=diff-3600._DP*hr-60._DP*min print'(i5,":",i2.2,":",f6.3)',hr,min,sec else call system_clock(itime) end if end subroutine stopwatch subroutine loopcounter(count,sizes,org,subs) implicit none integer,intent(IN) :: count,sizes(:),org(:) integer,intent(OUT) :: subs(:) integer :: dim,s,quo(size(subs)) dim=size(sizes) quo=count do s=2,dim; quo(s)=(quo(s-1)-1)/sizes(s-1)+1 end do do s=1,dim; subs(s)=mod(quo(s)-1,sizes(s))+org(s) end do end subroutine loopcounter subroutine rk_butcher(rk_A,rk_b,num) !----------------------------------------------------! ! prepare the Butcher matrix for Runge-Kutta method ! !----------------------------------------------------! implicit none real(DP),allocatable,dimension(:) :: rk_A,rk_b integer :: num select case (num) case(1) allocate(rk_A(1)); allocate(rk_b(2)) rk_A=(/1._DP/) rk_B=(/.5_DP,.5_DP/) case(2) allocate(rk_A(6)); allocate(rk_b(4)) rk_A=(/.5_DP,0._DP,.5_DP,0._DP,0._DP,1._DP/) rk_b=(/1._DP,2._DP,2._DP,1._DP/)/6._DP case(3) allocate(rk_A(6)); allocate(rk_b(4)) rk_A=(/1._DP/3._DP,-1._DP/3._DP,1._DP,1._DP,-1._DP,1._DP/) rk_b=(/1._DP,3._DP,3._DP,1._DP/)/8._DP case(4) allocate(rk_A(15)); allocate(rk_b(6)) rk_A=(/1._DP/4._DP,3._DP/32._DP,9._DP/32._DP,& 1932_DP/2197._DP,-7200._DP/2197._DP,7296._DP/2197._DP,& 439._DP/216._DP,-8._DP,3680._DP/513._DP,-845._DP/4104._DP,& -8._DP/27._DP,2._DP,-3544._DP/2565._DP,1859._DP/4104._DP,& -11._DP/40._DP/) rk_b=(/25._DP/216._DP,0._DP,1408._DP/2565._DP,& 2197._DP/4104._DP,-1._DP/5._DP,0._DP/) end select end subroutine rk_butcher function eq0(a) implicit none real(DP),intent(IN) :: a real(DP) :: eq0 eq0=0._DP**abs(a) ! eq0=1._DP-(1._DP+sign(1._DP,max(exp(a)-1._DP,0._DP)))*.5_DP end function eq0 function ne0(a) implicit none real(DP),intent(IN) :: a real(DP) :: ne0 ne0=1._DP-0._DP**abs(a) end function ne0 function gt0(a) implicit none real(DP),intent(IN) :: a real(DP) :: gt0 gt0=min(1._DP-0._DP**abs(a),(1._DP+sign(1._DP,a))*.5_DP) end function gt0 function ge0(a) implicit none real(DP),intent(IN) :: a real(DP) :: ge0 ge0=max(0._DP**abs(a),(1._DP+sign(1._DP,a))*.5_DP) end function ge0 function lt0(a) implicit none real(DP),intent(IN) :: a real(DP) :: lt0 lt0=min(1._DP-0._DP**abs(a),(1._DP-sign(1._DP,a))*.5_DP) end function lt0 function le0(a) implicit none real(DP),intent(IN) :: a real(DP) :: le0 le0=max(0._DP**abs(a),(1._DP-sign(1._DP,a))*.5_DP) end function le0 function Lag_interpol2(xm,xc,xp,fm,fc,fp,x) !2nd-order Lagrange interpolation implicit none real(DP),intent(IN) :: xm,xc,xp,fm,fc,fp,x real(DP) :: Lag_interpol2 Lag_interpol2=fm*(x-xc)*(x-xp)/(xm-xc)/(xm-xp)+& fc*(x-xm)*(x-xp)/(xc-xm)/(xc-xp)+& fp*(x-xm)*(x-xc)/(xp-xm)/(xp-xc) end function Lag_interpol2 function dfdx_c2(xm1,xc,xp1,fm1,fc,fp1) !2nd-order central difference implicit none real(DP),intent(IN) :: xm1,xc,xp1,fm1,fc,fp1 real(DP) :: alpha real(DP) :: hm1,hp1 real(DP) :: dfdx_c2 hm1=xc-xm1; hp1=xp1-xc dfdx_c2=((fp1-fc)*(hm1/hp1)+(fc-fm1)*(hp1/hm1))/(hm1+hp1) end function dfdx_c2 function dfdx_c4(xm2,xm1,xc,xp1,xp2,fm2,fm1,fc,fp1,fp2) !4th-order central difference implicit none real(DP),intent(IN) :: xm2,xm1,xc,xp1,xp2 real(DP),intent(IN) :: fm2,fm1,fc,fp1,fp2 real(DP) :: hm2,hm1,hp1,hp2 real(DP) :: dfm2,dfm1,dfp1,dfp2 real(DP) :: d2fm1,d2fc,d2fp1 real(DP) :: d3fm1,d3fp1,d4fc real(DP) :: dfdx_c4 hm2=xm1-xm2; hm1=xc-xm1; hp1=xp1-xc; hp2=xp2-xp1 hm2=hm2/hp1; hm1=hm1/hp1; hp2=hp2/hp1 dfm2=(fm1-fm2)/hm2; dfm1=(fc-fm1)/hm1 dfp1=fp1-fc; dfp2=(fp2-fp1)*hp2 d2fm1=2._DP/(hm2/hp1+hm1/hp1)*(dfm1-dfm2) d2fc=2._DP/(hm1/hp1+1._DP)*(dfp1-dfm1) d2fp1=2._DP/(1._DP+hp2/hp1)*(dfp2-dfp1) d3fm1=3._DP*(d2fc-d2fm1)/(hm2+hm1+1._DP) d3fp1=3._DP*(d2fp1-d2fc)/(hm1+1._DP+hp2) d4fc=4._DP*(d3fp1-d3fm1)/(hm2+hm1+1._DP+hp2) dfdx_c4=(dfp1-.5_DP*(d2fc-(d3fp1+.25_DP*(1+hp1)*d4fc)*hm1/3._DP))/hp1 end function dfdx_c4 end module sktk_fortran