地域密着型エリア広告配信リクルートの地域密着型広告ドコイク?アドネットワーク Ads by ドコイク?[無料でホームページを作成] [通報・削除依頼]
[無料でホームページを作成] [通報・削除依頼]

テトラヘドロン法

 ここではテトラヘドロン法について解説する。非常に重要なアルゴリズムとプログラムであるが、一般的に利用できるプログラムは容易に見つからない。筆者ももう若くないし、職も続かないので、若い人に頼みたい。
--------------------------------------------------------------------------------
■ 理論的背景の理解の仕方
  テトラヘドロン法の簡単な概略を知りたい場合は文献[1]を読むことをオススメしたい。計算式は文献[2]も参照しながら読み込んでいくとよい。
Reference
[1] H. Sato et al., Phil. Mag. 93 (2013) 3029.
[2] G.Lehmann and M.Taut, phys. stat. sol. (b) 54 (1972) 469.
--------------------------------------------------------------------------------
■ 引数
NPX,NPY, NPZ: それぞれh,k,l方向でのテトラへドロン法での分割数。
NP: 分割したメッシュに付けられた番号。
EOA: Electrons over atoms
I2: NPX, NPY, NPZで分割したメッシュの番号NPに対応するklistにある規約表現上の番号(IEX).
AMA: NPT(NPX, NPY, NPZで分割した総メッシュ数)に対応させたklistでの番号に対するエネルギー(eV) (AMA2)。
AKK: NPT(NPX, NPY, NPZで分割した総メッシュ数)に対応させたklistでの番号に対するC^2 (BFFF)。
AKKB: NPT(NPX, NPY, NPZで分割した総メッシュ数)に対応させたklistでの番号に対するK^2 (BKKK)。
S: エネルギー
F: C^2
FB: K^2
FC: (K^2)^2
CALL SORT: 4つあるSを昇順に並べた基準でデータが揃えられる。
NE1: エネルギー差分(S(1)-EM(1))/DE。もし、EBよりもS(1)が小さければNE1=-1
V: 四面体の体積
DELTA:
D1:
D2:
D4:
AM: テキスト中のMに該当する。
V1: V/DELTA
DEE: E-S(1 or 2 or 3 or 4)
CN: よくテキスト見かける式(ΔS/|∇E|)を計算する。
DF:
DF1:
DF2:
CN:
CN1:
CN2:
FN:C^2に(ΔS/|∇E|)を掛ける。
FNB:K^2に(ΔS/|∇E|)を掛ける。Kcgを計算。
FND:(K^2)^2に(ΔS/|∇E|)を掛ける。
SNT:SNの総和。CNT(つまり、DOS) をバンドの底からEまで積分した状態数。
CNT,CC:CNの総和。よくテキスト見かける式Σ(ΔS/|∇E|)を計算する。DOSが得られる。CCが0以下のときは割られる方は0にする。
FNT,C2:FNの総和
FNTB,C2B:FNBの総和
FNTC,C2C:FNCの総和
FNTD,C2D:FNDの総和
TNE: トータルの電子数
EE: エネルギー
CPP,CFC,CCC: K^2
CPP1,CFC1,CCC1: EOA
CP,CF: C^2

■ サブルーチンの解説(下記のプログラムは計算に時間が掛かるが、単純であり、どのような構造に対してもプログラムを変えずに計算が可能であるという利点がある)
CALL TGEN: outputkgenの形式にあわせてNPZ, NPY, NPXを変えていき、それに順番に番号をつける。そして、N1(NP)-N7(NP)までの七つの分類に分けて整理する。つまり、全体を六面体で分割し、さらに六面体を6つの四面体に分割したデータを作る。

CALL FDENS: TGENで分割した四面体の番号NPに対してPDNSで値を計算する。

CALL PDNS: 四面体(テトラへドロン)での計算をするルーチン。 SORTを呼び出して昇順に並べた後に処理が行われる。「■ 理論的背景の理解の仕方」で分かるように、幾つかの場合に分けて計算をしている。基本的に、[1] H. Sato et al., Phil. Mag. 93 (2013) 3029.と [2] G.Lehmann and M.Taut, phys. stat. sol. (b) 54 (1972) 469.の方法に従っているので式とプログラムの対応付けや理解がしやすいと思う。

CALL SORT: ソートのルーチン。昇順に並べる。
--------------------------------------------------------------------------------
■ プログラム

□ TGEN: outputkgenの形式にあわせてNPZ, NPY, NPXを変えていき、それに順番に番号をつける。そして、N1(NP)-N7(NP)までの七つの分類に分けて整理する。つまり、全体を六面体で分割し、さらに六面体を6つの四面体に分割したデータを作る。
!----------------------------------------------------------------------
subroutine TGEN(NPXD,NPYD,NPZD,IC,N1,N2,N3,N4,N5,N6,N7)
  use common
  implicit none
!---------------------------------------------------------------
!   *   DEFINE TETRAHEDRA ASSOCIATED WITH EACH K-POINT         *
!   *   THE K-SPACE MESH IMPLIED SHOULD BE CONSISTENT WITH THAT*
!   *   GENERATED BY SUBROUTINE <MESH> IN <STR>                *
!   *   SIMPLE CUBIC                                           *
!---------------------------------------------------------------
  integer(4),intent(in):: NPXD,NPYD,NPZD
  integer(4),intent(inout):: IC(NPXD,NPYD,NPZD)
  integer(4),intent(inout):: N1(NPT),N2(NPT),N3(NPT),N4(NPT),N5(NPT),N6(NPT),N7(NPT)
  !
  integer(4):: X,Y,Z
  integer(4):: NPXM,NPYM,NPZM
  integer(4):: NP,LP
  real(8):: DKX,DKY,DKZ
  !
  NPX=NPXD
  NPY=NPYD
  NPZ=NPZD
  !
  NPXM=NPX-1
  NPYM=NPY-1
  NPZM=NPZ-1
  !
  DKX=1.d0/(float(NPXM))
  DKY=1.d0/(float(NPYM))
  DKZ=1.d0/(float(NPZM))
  V=  DKX*DKY*DKZ / 3.0
  write(6,*) ' start!! PNT=', 'V=',V
!
! be careful! h,k,l start at 0, but fortran array start at 1.
! Please, remember X:1 = h:0, Y:1 = k:0 , Z:1 = l:0
!
  NP=0
  ! write(6,*) NPX, NPY, NPZ
  do X=1,NPX
    do Y=1,NPY
      do Z=1,NPZ
        NP=NP+1
        IC(X,Y,Z)=NP
      end do
    end do
  end do
  !
  ! write(6,*) ' total k-point',NP
  !
  NP=0
  DO X=1,NPX-1
    DO Y=1,NPY-1
      DO Z=1,NPZ-1
        NP=IC(X,Y,Z)
        N1(NP)=IC(X+1,Y  ,Z  )
        N2(NP)=IC(X+1,Y+1,Z  )
        N3(NP)=IC(X+1,Y+1,Z+1)
        N4(NP)=IC(X+1,Y  ,Z+1)
        N5(NP)=IC(X  ,Y  ,Z+1)
        N6(NP)=IC(X  ,Y+1,Z+1)
        N7(NP)=IC(X  ,Y+1,Z  )
      end do
    end do
  end do
  return
end subroutine TGEN
!----------------------------------------------------------------------

□ FDENS: TGENで分割した四面体の番号NPに対してPDNSで値を計算する。
!----------------------------------------------------------------------
subroutine FDENS(NPXD,NPYD,NPZD,IC,NPT,N1,N2,N3,N4,N5,N6,N7)
  implicit none
  integer(4),intent(in):: NPXD, NPYD, NPZD, NPT
  integer(4),intent(inout):: IC(NPXD,NPYD,NPZD)
  integer(4),intent(inout):: N1(NPT),N2(NPT),N3(NPT),N4(NPT),N5(NPT),N6(NPT),N7(NPT)
  integer(4):: I,J,K
  integer(4):: NPX,NPY,NPZ
  integer(4):: NP,LP
  integer(4):: I0,I1,I2,I3,I4,I5,I6,I7
  !
  NPX=NPXD
  NPY=NPYD
  NPZ=NPZD
  LP=0
  DO I=1,NPX-1
    DO J=1,NPY-1
      DO K=1,NPZ-1
        ! be careful! Tihs method is calculate tetrahedron including I0+1 range.
        ! NPX-1 is OK.
        ! NP=NP+1
        LP=IC(I,J,K)
        I0=LP
        I1=N1(LP)
        I2=N2(LP)
        I3=N3(LP)
        I4=N4(LP)
        I5=N5(LP)
        I6=N6(LP)
        I7=N7(LP)
        !
        CALL PDNS(I0,I1,I2,I3,LP)
        CALL PDNS(I0,I1,I3,I4,LP)
        CALL PDNS(I0,I3,I4,I5,LP)
        CALL PDNS(I0,I3,I5,I6,LP)
        CALL PDNS(I0,I3,I6,I7,LP)
        CALL PDNS(I0,I2,I3,I7,LP)
        !
      end do
    end do
  end do
  return
end subroutine FDENS
!----------------------------------------------------------------------

□ PDNS: 四面体(テトラへドロン)での計算をするルーチン。 SORTを呼び出して昇順に並べた後に処理が行われる。「■ 理論的背景の理解の仕方」で分かるように、幾つかの場合に分けて計算をしている。基本的に、[1] H. Sato et al., Phil. Mag. 93 (2013) 3029.と [2] G.Lehmann and M.Taut, phys. stat. sol. (b) 54 (1972) 469.の方法に従っているので式とプログラムの対応付けや理解がしやすいと思う。
!----------------------------------------------------------------------
subroutine PDNS(I1,I2,I3,I4,LP)
  use constants
  use common
  use allocate_data
  implicit none
!------------------------------------------------------------------
!   *   CALCULATE STATE DENSITY AND NUMBER OF STATES FOR A SINGLE *
!   *   TETRAHEDRON                                               *
!   *   CN : STATE DENSITY                                        *
!   *   SN : NUMBER OF STATES                                     *
!------------------------------------------------------------------
  !
  integer(4),intent(in):: I1,I2,I3,I4,LP
  !
  real(8):: S(4),F(4),FB(4),FC(4),G(4)
  real(8):: DELTA,AM,D1,D2,D4
  real(8):: V1,DEE,FRAC
  real(8):: CN,SN,SSS,CN1,CN2,DF,DF1,DF2
  real(8):: FN,FNB,FNC,FND
  real(8):: EE
  integer(4):: LE,NE1,NE2
  !
  G(1)=1.0
  G(2)=1.0
  G(3)=1.0
  G(4)=1.0
  !
  S(1)=E(I1)
  S(2)=E(I2)
  S(3)=E(I3)
  S(4)=E(I4)
  F(1)=FK(I1)
  F(2)=FK(I2)
  F(3)=FK(I3)
  F(4)=FK(I4)
  FB(1)=FKB(I1)
  FB(2)=FKB(I2)
  FB(3)=FKB(I3)
  FB(4)=FKB(I4)
  FC(1)=FKC(I1)
  FC(2)=FKC(I2)
  FC(3)=FKC(I3)
  FC(4)=FKC(I4)
  !
  CALL SORT(S,F,FB,FC)
  !
  if( S(1)>=(EM(NE)) ) return
  EB=EM(1)
  NE1=int ( (S(1)-EB)/DE )
  if( S(1)<EB ) NE1=-1
  NE1=NE1+2
  NE2=NEMAX
  DO 25 LE=NE1,NE2
    EE=EM(LE)
    if( EE <= S(1) ) goto 25
    if( EE >= S(4) ) goto 20
    if( EE <= S(2) ) goto 21
    if( EE >= S(3) ) goto 22
! CASE_2 E2 <E< E3
    DELTA=S(4)+S(3)-S(2)-S(1)
    AM=(S(4)*S(3)-S(2)*S(1))/DELTA
    D1=S(1)-AM
    D2=S(2)-AM
    V1=V/DELTA
    DEE=EE-AM
    FRAC=DEE*DEE/D2/D1
    CN=3.D0*V1*(1.D0-FRAC)
    SN=3.D0*V1*DEE-V1*(D1+D2+DEE*FRAC)
    if( ABS(S(2)-S(1)) <= 0.0 ) then
      SSS=3.D0*V*(EE-S(1))/(S(3)-S(1))/(S(4)-S(1))
      CN1=SSS*(S(3)-EE)/(S(3)-S(1))
      CN2=SSS*(S(4)-EE)/(S(4)-S(1))
      ! FN
      DF1=(F(2)-F(1))*( (S(3)-EE)/(S(3)-S(1))+ (S(4)-EE)/(S(4)-S(1)) ) &
        +2*(F(3)-F(1))*(EE-S(1))/(S(3)-S(1)) &
        +  (F(4)-F(1))*(EE-S(1))/(S(4)-S(1))
      DF2=(F(2)-F(1))*( (S(3)-EE)/(S(3)-S(1))+ (S(4)-EE)/(S(4)-S(1)) ) &
        +  (F(3)-F(1))*(EE-S(1))/(S(3)-S(1)) &
        +2*(F(4)-F(1))*(EE-S(1))/(S(4)-S(1))
      FN=F(1)*CN+ (DF1*CN1+DF2*CN2)/3.D0
      ! FNB
      DF1=(FB(2)-FB(1))*( (S(3)-EE)/(S(3)-S(1))+ (S(4)-EE)/(S(4)-S(1)) ) &
        +2*(FB(3)-FB(1))*(EE-S(1))/(S(3)-S(1)) &
        +  (FB(4)-FB(1))*(EE-S(1))/(S(4)-S(1))
      DF2=(FB(2)-FB(1))*( (S(3)-EE)/(S(3)-S(1))+ (S(4)-EE)/(S(4)-S(1)) ) &
        +  (FB(3)-FB(1))*(EE-S(1))/(S(3)-S(1)) &
        +2*(FB(4)-FB(1))*(EE-S(1))/(S(4)-S(1))
      FNB=FB(1)*CN+(DF1*CN1+DF2*CN2)/3.D0
      ! FNC
      DF1=(F(2)*FB(2)-F(1)*FB(1))*( (S(3)-EE)/(S(3)-S(1))+ (S(4)-EE)/(S(4)-S(1)) ) &
        +2*(F(3)*FB(3)-F(1)*FB(1))*(EE-S(1))/(S(3)-S(1)) &
        +  (F(4)*FB(4)-F(1)*FB(1))*(EE-S(1))/(S(4)-S(1))
      DF2=(F(2)*FB(2)-F(1)*FB(1))*( (S(3)-EE)/(S(3)-S(1))+ (S(4)-EE)/(S(4)-S(1)) ) &
        +  (F(3)*FB(3)-F(1)*FB(1))*(EE-S(1))/(S(3)-S(1)) &
        +2*(F(4)*FB(4)-F(1)*FB(1))*(EE-S(1))/(S(4)-S(1))
      FNC=F(1)*FB(1)*CN+(DF1*CN1+DF2*CN2)/3.D0
      ! FND
      DF1=(G(2)*FC(2)-G(1)*FC(1))*( (S(3)-EE)/(S(3)-S(1))+ (S(4)-EE)/(S(4)-S(1)) ) &
        +2*(G(3)*FC(3)-G(1)*FC(1))*(EE-S(1))/(S(3)-S(1)) &
        +  (G(4)*FC(4)-G(1)*FC(1))*(EE-S(1))/(S(4)-S(1))
      DF2=(G(2)*FC(2)-G(1)*FC(1))*( (S(3)-EE)/(S(3)-S(1))+ (S(4)-EE)/(S(4)-S(1)) ) &
        +  (G(3)*FC(3)-G(1)*FC(1))*(EE-S(1))/(S(3)-S(1)) &
        +2*(G(4)*FC(4)-G(1)*FC(1))*(EE-S(1))/(S(4)-S(1))
      FND=G(1)*FC(1)*CN+(DF1*CN1+DF2*CN2)/3.D0
      !
      goto 23
    end if
    if( ABS(S(3)-S(2)) <= 0.0) then
      ! FN
      DF=F(2)-F(1)+F(3)-F(1)+(F(4)-F(1))*(EE-S(1))/(S(4)-S(1))
      FN=F(1)*CN+DF*CN/3.D0
      ! FNB
      DF=FB(2)-FB(1)+FB(3)-FB(1)+(FB(4)-FB(1))*(EE-S(1))/(S(4)-S(1))
      FNB=FB(1)*CN+DF*CN/3.D0
      ! FNC
      DF=F(2)*FB(2)-F(1)*FB(1)+F(3)*FB(3)-F(1)*FB(1)+ &
      (F(4)*FB(4)-F(1)*FB(1))*(EE-S(1))/(S(4)-S(1))
      FNC=F(1)*FB(1)*CN+DF*CN/3.D0
      ! FND
      DF=G(2)*FC(2)-G(1)*FC(1)+G(3)*FC(3)-G(1)*FC(1)+ &
      (G(4)*FC(4)-G(1)*FC(1))*(EE-S(1))/(S(4)-S(1))
      FND=G(1)*FC(1)*CN+DF*CN/3.D0
      !
      goto 23
    end if
! CN1=small tri, CN2=big tri
!     CN2=3.D0*V1*(-DDE*DDE/D1+2.D0*DEE-D1) /(S(2)-S(1))
    CN2=3.D0*V *(EE-S(1))**2/(S(2)-S(1))/(S(3)-S(1))/(S(4)-S(1))
    CN1=CN2-CN
    ! FN
    DF2=( (F(2)-F(1))/(S(2)-S(1)) +(F(3)-F(1))/(S(3)-S(1)) &
      +(F(4)-F(1))/(S(4)-S(1)) )*(EE-S(1))
    DF1=(F(2)-F(1))*( (EE-S(1))/(S(2)-S(1))+(S(3)-EE)/(S(3)-S(2)) &
      +(S(4)-EE)/(S(4)-S(2)) ) &
      +( (F(3)-F(1))/(S(3)-S(2))+(F(4)-F(1))/(S(4)-S(2)) )*(EE-S(2))
    FN=F(1)*CN +(DF2*CN2-DF1*CN1)/3.D0
    ! FNB
    DF2=( (FB(2)-FB(1))/(S(2)-S(1)) +(FB(3)-FB(1))/(S(3)-S(1)) &
      +(FB(4)-FB(1))/(S(4)-S(1)) )*(EE-S(1))
    DF1=(FB(2)-FB(1))*( (EE-S(1))/(S(2)-S(1))+(S(3)-EE)/(S(3)-S(2)) &
      +(S(4)-EE)/(S(4)-S(2)) ) &
      +( (FB(3)-FB(1))/(S(3)-S(2))+(FB(4)-FB(1))/(S(4)-S(2)) )*(EE-S(2))
    FNB=FB(1)*CN+(DF2*CN2-DF1*CN1)/3.D0
    ! FNC
    DF2=( (F(2)*FB(2)-F(1)*FB(1))/(S(2)-S(1)) +(F(3)*FB(3)-F(1)*FB(1))/(S(3)-S(1)) &
      +(F(4)*FB(4)-F(1)*FB(1))/(S(4)-S(1)) )*(EE-S(1))
    DF1=(F(2)*FB(2)-F(1)*FB(1))*( (EE-S(1))/(S(2)-S(1))+(S(3)-EE)/(S(3)-S(2)) &
      +(S(4)-EE)/(S(4)-S(2)) ) &
      +( (F(3)*FB(3)-F(1)*FB(1))/(S(3)-S(2))+(F(4)*FB(4)-F(1)*FB(1))/ &
      (S(4)-S(2)) )*(EE-S(2))
    FNC=F(1)*FB(1)*CN+(DF2*CN2-DF1*CN1)/3.D0
    ! FND
    DF2=( (G(2)*FC(2)-G(1)*FC(1))/(S(2)-S(1)) +(G(3)*FC(3)-G(1)*FC(1))/(S(3)-S(1)) &
      +(G(4)*FC(4)-G(1)*FC(1))/(S(4)-S(1)) )*(EE-S(1))
    DF1=(G(2)*FC(2)-G(1)*FC(1))*( (EE-S(1))/(S(2)-S(1))+(S(3)-EE)/(S(3)-S(2)) &
      +(S(4)-EE)/(S(4)-S(2)) ) &
      +( (G(3)*FC(3)-G(1)*FC(1))/(S(3)-S(2))+(G(4)*FC(4)-G(1)*FC(1))/ &
      (S(4)-S(2)) )*(EE-S(2))
    FND=G(1)*FC(1)*CN+(DF2*CN2-DF1*CN1)/3.D0
    !
    goto 23
! CASE_3 E1<E2< E3<E<E4
 22 if( ABS(S(4)-S(3)) <= 0.0 ) then
      CN=0.0
      SN=0.0
      FN=0.0
      FNB=0.0
      FNC=0.0
      FND=0.0
      goto 23
    end if
    DELTA=S(4)+S(3)-S(2)-S(1)
    AM=(S(4)*S(3)-S(2)*S(1))/DELTA
    D4=S(4)-AM
    V1=V/DELTA
    DEE=EE-S(4)
    FRAC=DEE*DEE/D4/(S(4)-S(3))
    CN=3.D0*V1*FRAC
    SN=V1*(DELTA+DEE*FRAC)
    ! FN
    DF=(F(2)-F(1))*(S(4)-EE)/(S(4)-S(2)) &
        +(F(3)-F(1))*(S(4)-EE)/(S(4)-S(3)) &
        +(F(4)-F(1))*( (EE-S(1))/(S(4)-S(1)) &
        +(EE-S(2))/(S(4)-S(2))+(EE-S(3))/(S(4)-S(3)) )
    FN=F(1)*CN+DF*CN/3.D0
    ! FNB
    DF=(FB(2)-FB(1))*(S(4)-EE)/(S(4)-S(2)) &
        +(FB(3)-FB(1))*(S(4)-EE)/(S(4)-S(3)) &
        +(FB(4)-FB(1))*( (EE-S(1))/(S(4)-S(1)) &
        +(EE-S(2))/(S(4)-S(2))+(EE-S(3))/(S(4)-S(3)) )
    FNB=FB(1)*CN+DF*CN/3.D0
    ! FNC
    DF=(F(2)*FB(2)-F(1)*FB(1))*(S(4)-EE)/(S(4)-S(2)) &
        +(F(3)*FB(3)-F(1)*FB(1))*(S(4)-EE)/(S(4)-S(3)) &
        +(F(4)*FB(4)-F(1)*FB(1))*( (EE-S(1))/(S(4)-S(1)) &
        +(EE-S(2))/(S(4)-S(2))+(EE-S(3))/(S(4)-S(3)) )
    FNC=F(1)*FB(1)*CN+DF*CN/3.D0
    ! FND
    DF=(G(2)*FC(2)-G(1)*FC(1))*(S(4)-EE)/(S(4)-S(2)) &
        +(G(3)*FC(3)-G(1)*FC(1))*(S(4)-EE)/(S(4)-S(3)) &
        +(G(4)*FC(4)-G(1)*FC(1))*( (EE-S(1))/(S(4)-S(1)) &
        +(EE-S(2))/(S(4)-S(2))+(EE-S(3))/(S(4)-S(3)) )
    FND=G(1)*FC(1)*CN+DF*CN/3.D0
    !
    goto 23
! CASE_1 E1<E<E2 <E3<E4
 21 if( ABS(S(2)-S(1)) <= 0.0 ) then
      CN=0.0
      SN=0.0
      FN=0.0
      FNB=0.0
      FNC=0.0
      FND=0.0
      goto 23
    endif
    DELTA=S(4)+S(3)-S(2)-S(1)
    AM=(S(4)*S(3)-S(2)*S(1))/DELTA
    D1=S(1)-AM
    V1=V/DELTA
    DEE=EE-S(1)
    FRAC=DEE*DEE/D1/(S(1)-S(2))
    CN=3.D0*V1*FRAC
    SN=V1*FRAC*DEE
    ! FN
    DF=( (F(2)-F(1))/(S(2)-S(1))+ (F(3)-F(1))/(S(3)-S(1)) &
           +(F(4)-F(1))/(S(4)-S(1)) )*DEE
    FN=F(1)*CN+DF*CN/3.D0
    ! FNB
    DF=( (FB(2)-FB(1))/(S(2)-S(1))+ (FB(3)-FB(1))/(S(3)-S(1)) &
           +(FB(4)-FB(1))/(S(4)-S(1)) )*DEE
    FNB=FB(1)*CN+DF*CN/3.D0
    ! FNC
    DF=( (F(2)*FB(2)-F(1)*FB(1))/(S(2)-S(1))+ (F(3)*FB(3)-F(1)*FB(1))/(S(3)-S(1)) &
           +(F(4)*FB(4)-F(1)*FB(1))/(S(4)-S(1)) )*DEE
    FNC=F(1)*FB(1)*CN+DF*CN/3.D0
    ! FND
    DF=( (G(2)*FC(2)-G(1)*FC(1))/(S(2)-S(1))+ (G(3)*FC(3)-G(1)*FC(1))/(S(3)-S(1)) &
           +(G(4)*FC(4)-G(1)*FC(1))/(S(4)-S(1)) )*DEE
    FND=G(1)*FC(1)*CN+DF*CN/3.D0
    !
    goto 23
!
 20 SN =V
    CN =0.d0
    FN =0.d0
    FNB=0.d0
    FNC=0.d0
    FND=0.d0
!
 23 SNT(LE) =SNT(LE) +SN !*IWW(LP)
    CNT(LE) =CNT(LE) +CN !*IWW(LP)
    FNT(LE) =FNT(LE) +FN !*IWW(LP)
    FNTB(LE)=FNTB(LE)+FNB !*IWW(LP)
    FNTC(LE)=FNTC(LE)+FNC !*IWW(LP)
    FNTD(LE)=FNTD(LE)+FND !*IWW(LP)
 25 CONTINUE
  return
end subroutine PDNS
!----------------------------------------------------------------------

□ SORT: ソートのルーチン。昇順に並べる。
!----------------------------------------------------------------------
subroutine SORT(S,F,G,O)
  implicit none
  integer(4):: I,J,NP
  real(8),intent(inout):: S(4),F(4),G(4),O(4)
  real(8):: ST, FT, GT, OT
  DO J=1,3
    NP=4-J
    DO I=1,NP
      if(S(I+1)<S(I)) then
        ST=S(I)
        FT=F(I)
        GT=G(I)
        OT=O(I)
        !
        S(I)=S(I+1)
        F(I)=F(I+1)
        G(I)=G(I+1)
        O(I)=O(I+1)
        !
        S(I+1)=ST
        F(I+1)=FT
        G(I+1)=GT
        O(I+1)=OT
      end if
    end do
  end do
  return
end subroutine SORT
!----------------------------------------------------------------------
--------------------------------------------------------------------------------
QRコード
携帯用QRコード
アクセス数
ページビュー数