text:
!> 装载或产生初始数据的子程序。
!> subroutine for loading or generating initial particle information
subroutine input(x, vx, mass, rho, p, u, itype, hsml, ntotal)

    use sph_kind, only: rk
    use parameter
    implicit none

    !> 粒子的位置
    !> coordinates of particles
    real(rk), intent(out) :: x(dim, maxn)
    !> 粒子的速度
    !> velocities of particles
    real(rk), intent(out) :: vx(dim, maxn)
    !> 粒子的质量
    !> mass of particles
    real(rk), intent(out) :: mass(maxn)
    !> 粒子的密度
    !> dnesities of particles
    real(rk), intent(out) :: rho(maxn)
    !> 粒子的压力
    !> pressure  of particles
    real(rk), intent(out) :: p(maxn)
    !> 粒子的内部能量
    !> internal energy of particles
    real(rk), intent(out) :: u(maxn)
    !> 粒子的类型
    !> types of particles
    integer, intent(out) :: itype(maxn)
    !> 粒子的光滑长度
    !> smoothing lengths of particles
    real(rk), intent(out) :: hsml(maxn)
    !> 在模拟中所使用的粒子总数
    !> number of particles in simulation
    integer, intent(out) :: ntotal

    integer :: i, d, im

    !> load initial particle information from external disk file

    if (config_input) then

        open (1, file='./data/f_xv.dat')
        open (2, file='./data/f_state.dat')
        open (3, file='./data/f_other.dat')

        write (*, *) '  **************************************************'
        write (*, *) '      loading initial particle configuration...   '
        read (1, *) ntotal
        write (*, *) '      total number of particles   ', ntotal
        write (*, *) '  **************************************************'
        do i = 1, ntotal
            read (1, *) im, (x(d, i), d=1, dim), (vx(d, i), d=1, dim)
            read (2, *) im, mass(i), rho(i), p(i), u(i)
            read (3, *) im, itype(i), hsml(i)
        end do

    else

        open (1, file='./data/ini_xv.dat')
        open (2, file='./data/ini_state.dat')
        open (3, file='./data/ini_other.dat')

        if (shocktube) call shock_tube(x, vx, mass, rho, p, u, itype, hsml, ntotal)

        if (shearcavity) call shear_cavity(x, vx, mass, rho, p, u, itype, hsml, ntotal)
        do i = 1, ntotal
            write (1, 1001) i, (x(d, i), d=1, dim), (vx(d, i), d=1, dim)
            write (2, 1002) i, mass(i), rho(i), p(i), u(i)
            write (3, 1003) i, itype(i), hsml(i)
        end do
        write (*, *) '  **************************************************'
        write (*, *) '      initial particle configuration generated   '
        write (*, *) '      total number of particles   ', ntotal
        write (*, *) '  **************************************************'

    end if

    close (1)
    close (2)
    close (3)

1001 format(1x, i5, 6(2x, e14.8))
1002 format(1x, i5, 7(2x, e14.8))
1003 format(1x, i5, 2x, i2, 2x, e14.8)

end subroutine input

!> 一维振荡管的初始数据。
!> this subroutine is used to generate initial data for the
!> 1 d noh shock tube problem
subroutine shock_tube(x, vx, mass, rho, p, u, itype, hsml, ntotal)

    use sph_kind, only: rk
    use parameter
    implicit none

    !> 粒子的位置
    !> coordinates of particles
    real(rk), intent(out) :: x(dim, maxn)
    !> 粒子的速度
    !> velocities of particles
    real(rk), intent(out) :: vx(dim, maxn)
    !> 粒子的质量
    !> mass of particles
    real(rk), intent(out) :: mass(maxn)
    !> 粒子的密度
    !> dnesities of particles
    real(rk), intent(out) :: rho(maxn)
    !> 粒子的压力
    !> pressure  of particles
    real(rk), intent(out) :: p(maxn)
    !> 粒子的内部能量
    !> internal energy of particles
    real(rk), intent(out) :: u(maxn)
    !> 粒子的类型(1: ideal gas; 2: water)
    !> types of particles
    integer, intent(out) :: itype(maxn)
    !> 粒子的光滑长度
    !> smoothing lengths of particles
    real(rk), intent(out) :: hsml(maxn)
    !> 在模拟中所使用的粒子总数
    !> number of particles in simulation
    integer, intent(out) :: ntotal

    integer :: i, d
    real(rk) :: space_x

    ntotal = 400
    space_x = 0.6_rk/80._rk

    do i = 1, ntotal
        mass(i) = 0.75_rk/400._rk
        hsml(i) = 0.015_rk
        itype(i) = 1
        do d = 1, dim
            x(d, i) = 0._rk
            vx(d, i) = 0._rk
        end do
    end do

    do i = 1, 320
        x(1, i) = -0.6_rk + space_x/4._rk*(i - 1)
    end do

    do i = 320 + 1, ntotal
        x(1, i) = 0._rk + space_x*(i - 320)
    end do

    do i = 1, ntotal
        if (x(1, i) <= 1.e-8_rk) then
            u(i) = 2.5_rk
            rho(i) = 1._rk
            p(i) = 1._rk
        end if
        if (x(1, i) > 1.e-8_rk) then
            u(i) = 1.795_rk
            rho(i) = 0.25_rk
            p(i) = 0.1795_rk
        end if
    end do

end subroutine shock_tube

!> 二维剪切腔的初始数据。
!> this subroutine is used to generate initial data for the
!> 2 d shear driven cavity probem with re = 1
subroutine shear_cavity(x, vx, mass, rho, p, u, itype, hsml, ntotal)

    use sph_kind, only: rk
    use parameter
    implicit none

    !> 粒子的位置
    !> coordinates of particles
    real(rk), intent(out) :: x(dim, maxn)
    !> 粒子的速度
    !> velocities of particles
    real(rk), intent(out) :: vx(dim, maxn)
    !> 粒子的质量
    !> mass of particles
    real(rk), intent(out) :: mass(maxn)
    !> 粒子的密度
    !> dnesities of particles
    real(rk), intent(out) :: rho(maxn)
    !> 粒子的压力
    !> pressure  of particles
    real(rk), intent(out) :: p(maxn)
    !> 粒子的内部能量
    !> internal energy of particles
    real(rk), intent(out) :: u(maxn)
    !> 粒子的类型(1: ideal gas; 2: water)
    !> types of particles
    integer, intent(out) :: itype(maxn)
    !> 粒子的光滑长度
    !> smoothing lengths of particles
    real(rk), intent(out) :: hsml(maxn)
    !> 在模拟中所使用的粒子总数
    !> number of particles in simulation
    integer, intent(out) :: ntotal

    integer :: i, j, d, m, n, mp, np, k
    real(rk) :: xl, yl, dx, dy

    !> giving mass and smoothing length as well as other data.

    m = 41
    n = 41
    mp = m - 1
    np = n - 1
    ntotal = mp*np
    xl = 1.e-3_rk
    yl = 1.e-3_rk
    dx = xl/mp
    dy = yl/np

    do i = 1, mp
        do j = 1, np
            k = j + (i - 1)*np
            x(1, k) = (i - 1)*dx + dx/2._rk
            x(2, k) = (j - 1)*dy + dy/2._rk
        end do
    end do

    do i = 1, mp*np
        vx(1, i) = 0._rk
        vx(2, i) = 0._rk
        rho(i) = 1000._rk
        mass(i) = dx*dy*rho(i)
        p(i) = 0._rk
        u(i) = 357.1_rk
        itype(i) = 2
        hsml(i) = dx
    end do

end subroutine shear_cavity

max_stars_repo_name: zoziha/SPH

max_issues_repo_name: zoziha/SPH

max_forks_repo_name: zoziha/SPH


PROGRAM ufld2gp_weight_sum

  implicit none
  integer,parameter :: u1=1, u2=10, u3=11
  integer,parameter :: maxloop=1000000
  integer :: loop,ik,k,n,i,j,s,ne,nk,nb,ns,kmax,nmax,nmin
  integer :: nkp(0:9),nbp(0:9),nsp(0:9),kmaxp(0:9),nmaxp(0:9),nminp(0:9)
  integer :: kminp(0:9),k0
  integer :: loop_read, n_loop_read
  integer :: unit_conv=1
  character(1) :: c1
  character(5) :: cbuf
  character(10) :: file_name
  character(128) :: cbuf_long
  real(8),parameter :: aB_=0.529177d0, HT_=27.2116d0
  real(8) :: aB, HT
  real(8),allocatable :: kxyz_pc(:,:),kxyz_sc(:,:)
  real(8),allocatable :: esp(:,:,:), weight(:,:,:)
  real(8) :: emax,emin,e_mergin,eta,de,pi,f(2),e,x,dummy(4),d
  real(8) :: weight_tmp(2)
  logical :: flag

! ---

  ne = 500

  eta = 0.10d0/27.2116d0

  pi = acos(-1.0d0)

! ---

  write(*,*) "eta(eV), ne="
  read(*,*) eta, ne
  eta = eta/HT_

! ---

  nsp(:)=0
  nbp(:)=0
  nkp(:)=0
  kminp(:)=1000000000
  kmaxp(:)=0
  nmaxp(:)=0
  nminp(:)=1000000000

  do loop_read=0,9

     if ( loop_read == 0 ) then
        file_name="band_ufld"
     else
        write(c1,'(i1)') loop_read
        file_name = "band_ufld"//c1
     end if

     inquire( FILE=file_name, EXIST=flag )
     if ( .not.flag ) exit

     write(*,*) loop_read,file_name

     open(u1,file=file_name,status="old")

     do loop=1,maxloop
        read(u1,*,END=1) cbuf,k
        if ( cbuf == "iktrj" ) then
           kminp(loop_read)=min(k,kminp(loop_read))
           kmaxp(loop_read)=max(k,kmaxp(loop_read))
           nkp(loop_read)=nkp(loop_read)+1
           if ( nbp(loop_read) /= 0 ) then
              nminp(loop_read)=min(nminp(loop_read),nbp(loop_read))
           end if
           nbp(loop_read)=0
        else
           nbp(loop_read)=nbp(loop_read)+1
           nmaxp(loop_read)=max(nmaxp(loop_read),nbp(loop_read))
        end if
     end do ! loop
1    continue

     rewind u1
     read(u1,*)
     read(u1,'(a)') cbuf_long
     read(cbuf_long,*,END=2) i,j,dummy(1:4)
     nsp(loop_read)=2
2    continue
     if ( nsp(loop_read) == 0 ) nsp(loop_read)=1

!     if ( loop_read > 0 .and. kminp(loop_read) == 1 ) then
!        kminp(loop_read)=kminp(loop_read)+kmaxp(loop_read-1)
!        kmaxp(loop_read)=kmaxp(loop_read)+kmaxp(loop_read-1)
!     end if

     write(*,'(a50,"loop_read=",i1)') repeat("-",50),loop_read
     write(*,*) file_name, loop_read
     write(*,*) "ns, nk, nb=",nsp(loop_read),nkp(loop_read),nbp(loop_read)
     write(*,*) "kmin,max=",kminp(loop_read),kmaxp(loop_read)
     write(*,*) "nmin,nmax=",nminp(loop_read),nmaxp(loop_read)

     close(u1)

  end do ! loop_read

  n_loop_read=loop_read-1

!  nk=sum(nkp)
  nk=maxval(nkp) ; if ( any(nkp(1:n_loop_read)/=nk) ) stop "stop(1)"
  ns=maxval(nsp)
  nmin=minval(nmaxp(0:n_loop_read))
  nmax=maxval(nmaxp)

  write(*,*) "ns, nk=",ns,nk
  write(*,*) "nmin,nmax=",nmin,nmax

! ---

  allocate( kxyz_pc(3,0:nk)    ) ; kxyz_pc=0.0d0
  allocate( kxyz_sc(3,nk)      ) ; kxyz_sc=0.0d0
  allocate( esp(nmax,nk,ns)    ) ; esp=0.0d0
  allocate( weight(nmax,nk,ns) ) ; weight=0.0d0

! ---

  do loop_read=0,n_loop_read

     if ( loop_read == 0 ) then
        file_name="band_ufld"
     else
        write(c1,'(i1)') loop_read
        file_name = "band_ufld"//c1
     end if

     k0=0
!     if ( loop_read > 0 ) k0=kmaxp(loop_read-1)

     open(u1,file=file_name,status="old")

     do ik=1,nkp(loop_read)

        read(u1,*) cbuf,k,kxyz_pc(1:3,k+k0),kxyz_sc(1:3,k+k0)

        k=k+k0

        if ( k == kminp(loop_read) .or. k == kmaxp(loop_read) ) then
           write(*,'(1x,2i4,2(3f15.8,2x))') &
                loop_read, k, kxyz_pc(:,k), kxyz_sc(:,k)
        end if

        do n=1,nmaxp(loop_read)
           read(u1,*) i,j,(esp(n,k,s),weight_tmp(s),s=1,nsp(loop_read))
           do s=1,nsp(loop_read)
              weight(n,k,s)=weight(n,k,s)+weight_tmp(s)
           end do
           read(u1,*,END=9) cbuf
           backspace(u1)
           if ( cbuf == "iktrj" ) exit
        end do

     end do
9    continue

!     if ( loop_read > 0 ) then
!        d=sum((kxyz_pc(:,kminp(loop_read))-kxyz_pc(:,kmaxp(loop_read-1)))**2)
!        if ( d < 1.d-8 ) then
!           write(*,*) "omit the first k points"
!           do k=kminp(loop_read),kmaxp(loop_read)-1
!              kxyz_pc(:,k)=kxyz_pc(:,k+1)
!              esp(:,k,:)=esp(:,k+1,:)
!              weight(:,k,:)=weight(:,k+1,:)
!           end do
!           kminp(loop_read)=kminp(loop_read)+1
!           nkp(loop_read)=nkp(loop_read)-1
!        end if
!     end if

     close(u1)

  end do ! loop_read

!  nk=sum(nkp)
!  write(*,*) "nk(new)=",nk

! ---

  emin = minval( esp(1:nmin,:,:) )
  emax = maxval( esp(1:nmin,:,:) )

  write(*,*) "emin, emax=",emin,emax

  e_mergin = (emax-emin)*0.1d0

  emax = emax + e_mergin
  emin = emin - e_mergin

  write(*,*) "emin, emax, e_mergin=",emin,emax,e_mergin

  de = (emax-emin)/ne

  write(*,*) "ne, de=",ne,de

  write(*,*) "eta=",eta

! ---

  write(*,*) "Chose the unit [ 0: Hartree a.u.,  1: eV & angstrome^-1 ]"
  read(*,*) unit_conv

  if ( unit_conv == 0 ) then
     aB=1.0d0
     HT=1.0d0
  else
     aB=aB_
     HT=HT_
  end if

! ---

  rewind u2
  rewind u3

  kxyz_pc(:,0) = kxyz_pc(:,1)
  x=0.0d0

  do k=1,nk

     x = x + sqrt( sum( (kxyz_pc(:,k)-kxyz_pc(:,k-1))**2 ) )/aB

     do i=1,ne

        e = emin + (i-1)*de

        do s=1,ns
           f(s)=0.0d0
           do n=1,nmin
              f(s) = f(s) + eta/( (e-esp(n,k,s))**2 + eta**2 )*weight(n,k,s)
           end do
           f(s)=f(s)/pi
        end do

        write(u2,'(1x,4f25.15)') x,e*HT,(f(s),s=1,ns)
        write(u3,'(1x,4f25.15)') x,e*HT,(log(f(s)),s=1,ns)

     end do ! i

     write(u2,*)
     write(u3,*)

  end do ! k

! ---

  deallocate( weight  )
  deallocate( esp     )
  deallocate( kxyz_sc )
  deallocate( kxyz_pc )

END PROGRAM ufld2gp_weight_sum


text:
C PR rtl-optimization/42294
C { dg-do compile { target powerpc*-*-* ia64-*-* i?86-*-* x86_64-*-* } }
C { dg-options "-O2 -fselective-scheduling2 -fsel-sched-pipelining -funroll-all-loops" }

      SUBROUTINE ORIEN(IW,NATOT,NTOTORB,NATORB,P,T)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION NATORB(NATOT),P(NTOTORB*(NTOTORB+1)/2)
      DIMENSION T(NTOTORB,NTOTORB)
      DO 9000 IATOM=1,NATOT
         ILAST = NTOTORB
         IF (IATOM.NE.NATOT) ILAST=NATORB(IATOM+1)-1
         DO 8000 IAOI=NATORB(IATOM),ILAST
            DO 7000 IAOJ = IAOI+1,ILAST
               R2 = 0.0D+00
               R3 = 0.0D+00
               DO 6000 INOTA=1,NATOT
                  DO 5000 IK=NATORB(INOTA),NTOTORB
                     IMAI=MAX(IK,IAOI)
                     IMII=MIN(IK,IAOI)
                     IMAJ=MAX(IK,IAOJ)
                     IMIJ=MIN(IK,IAOJ)
                     IKI=(IMAI*(IMAI-1))/2 + IMII
                     IKJ=(IMAJ*(IMAJ-1))/2 + IMIJ
                     PIKI=P(IKI)
                     PIKJ=P(IKJ)
                     R2 = R2 + (PIKI**4)-6*(PIKI*PIKI*PIKJ*PIKJ)+(PIKJ)
 5000             CONTINUE
 6000          CONTINUE
               R2 = (R2/4.0D+00)
               Q = SQRT(R2*R2 + R3*R3)
               IF (Q.LT.1.0D-08) GO TO 7000
               A = COS(THETA)
               B = -SIN(THETA)
               CALL ROT1INT(NTOTORB,IAOI,IAOJ,A,B,P)
 7000       CONTINUE
 8000    CONTINUE
 9000 CONTINUE
      RETURN
      END

ext: f


program Test1075_2         !同为试验1075.f90中的功能而编写的代码
real :: a=1,b=2
print*, 1.0/2.0
end


text:
*DECK DBESK1
      DOUBLE PRECISION FUNCTION DBESK1 (X)
C***BEGIN PROLOGUE  DBESK1
C***PURPOSE  Compute the modified (hyperbolic) Bessel function of the
C            third kind of order one.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C10B1
C***TYPE      DOUBLE PRECISION (BESK1-S, DBESK1-D)
C***KEYWORDS  FNLIB, HYPERBOLIC BESSEL FUNCTION,
C             MODIFIED BESSEL FUNCTION, ORDER ONE, SPECIAL FUNCTIONS,
C             THIRD KIND
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C DBESK1(X) calculates the double precision modified (hyperbolic)
C Bessel function of the third kind of order one for double precision
C argument X.  The argument must be large enough that the result does
C not overflow and small enough that the result does not underflow.
C
C Series for BK1        on the interval  0.          to  4.00000E+00
C                                        with weighted error   9.16E-32
C                                         log weighted error  31.04
C                               significant figures required  30.61
C                                    decimal places required  31.64
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  D1MACH, DBESI1, DBSK1E, DCSEVL, INITDS, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770701  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C***END PROLOGUE  DBESK1
      DOUBLE PRECISION X, BK1CS(16), XMAX, XMAXT, XMIN, XSML, Y,
     1  D1MACH, DCSEVL, DBESI1, DBSK1E
      LOGICAL FIRST
      SAVE BK1CS, NTK1, XMIN, XSML, XMAX, FIRST
      DATA BK1CS(  1) / +.2530022733 8947770532 5311208685 33 D-1     /
      DATA BK1CS(  2) / -.3531559607 7654487566 7238316918 01 D+0     /
      DATA BK1CS(  3) / -.1226111808 2265714823 4790679300 42 D+0     /
      DATA BK1CS(  4) / -.6975723859 6398643501 8129202960 83 D-2     /
      DATA BK1CS(  5) / -.1730288957 5130520630 1765073689 79 D-3     /
      DATA BK1CS(  6) / -.2433406141 5659682349 6007350301 64 D-5     /
      DATA BK1CS(  7) / -.2213387630 7347258558 3152525451 26 D-7     /
      DATA BK1CS(  8) / -.1411488392 6335277610 9583302126 08 D-9     /
      DATA BK1CS(  9) / -.6666901694 1993290060 8537512643 73 D-12    /
      DATA BK1CS( 10) / -.2427449850 5193659339 2631968648 53 D-14    /
      DATA BK1CS( 11) / -.7023863479 3862875971 7837971200 00 D-17    /
      DATA BK1CS( 12) / -.1654327515 5100994675 4910293333 33 D-19    /
      DATA BK1CS( 13) / -.3233834745 9944491991 8933333333 33 D-22    /
      DATA BK1CS( 14) / -.5331275052 9265274999 4666666666 66 D-25    /
      DATA BK1CS( 15) / -.7513040716 2157226666 6666666666 66 D-28    /
      DATA BK1CS( 16) / -.9155085717 6541866666 6666666666 66 D-31    /
      DATA FIRST /.TRUE./
C***FIRST EXECUTABLE STATEMENT  DBESK1
      IF (FIRST) THEN
         NTK1 = INITDS (BK1CS, 16, 0.1*REAL(D1MACH(3)))
         XMIN = EXP(MAX(LOG(D1MACH(1)), -LOG(D1MACH(2))) + 0.01D0)
         XSML = SQRT(4.0D0*D1MACH(3))
         XMAXT = -LOG(D1MACH(1))
         XMAX = XMAXT - 0.5D0*XMAXT*LOG(XMAXT)/(XMAXT+0.5D0)
      ENDIF
      FIRST = .FALSE.
C
      IF (X .LE. 0.D0) CALL XERMSG ('SLATEC', 'DBESK1',
     +   'X IS ZERO OR NEGATIVE', 2, 2)
      IF (X.GT.2.0D0) GO TO 20
C
      IF (X .LT. XMIN) CALL XERMSG ('SLATEC', 'DBESK1',
     +   'X SO SMALL K1 OVERFLOWS', 3, 2)
      Y = 0.D0
      IF (X.GT.XSML) Y = X*X
      DBESK1 = LOG(0.5D0*X)*DBESI1(X) + (0.75D0 + DCSEVL (.5D0*Y-1.D0,
     1  BK1CS, NTK1))/X
      RETURN
C
 20   DBESK1 = 0.D0
      IF (X .GT. XMAX) CALL XERMSG ('SLATEC', 'DBESK1',
     +   'X SO BIG K1 UNDERFLOWS', 1, 1)
      IF (X.GT.XMAX) RETURN
C
      DBESK1 = EXP(-X) * DBSK1E(X)
C
      RETURN
      END


text: #redirect Davis Auto Works

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


! READ function

! read (unit=1,rec=1) a => READ (UNIT=1,FMT=*,REC=1) A
!
! Compilation error: FMT=* is not allowed with a REC= specifier at (1)

! Bug: the function return type replaced by "OVERLOADED"

program main

    implicit none
    integer :: a

    a = 0
    open(unit=1,file='main08.f90',form='unformatted', access='direct',status='old',recl=1000)

    read(1,rec=1) a

end program main


text:
module constants_mod
    implicit none
    integer, parameter :: str_len = 20
end module constants_mod

max_stars_repo_name: bblay/fab

max_issues_repo_name: bblay/fab

max_forks_repo_name: bblay/fab


text: guero is now Users/ashiggins.


text:
      FUNCTION RAN4(IDUM)
      PARAMETER (IM=11979,IA=430,IC=2531,NACC=24)
      DIMENSION INP(64),JOT(64),KEY(64),POW(65)
      DATA IFF/0/
      IF(IDUM.LT.0.OR.IFF.EQ.0)THEN
        IFF=1
        IDUM=MOD(IDUM,IM)
        POW(1)=0.5
        DO 11 J=1,64
          IDUM=MOD(IDUM*IA+IC,IM)
          KEY(J)=(2*IDUM)/IM
          INP(J)=MOD((4*IDUM)/IM,2)
          POW(J+1)=0.5*POW(J)
11      CONTINUE
        NEWKEY=1
      ENDIF
      ISAV=INP(64)
      IF(ISAV.NE.0)THEN
        INP(4)=1-INP(4)
        INP(3)=1-INP(3)
        INP(1)=1-INP(1)
      ENDIF
      DO 12 J=64,2,-1
        INP(J)=INP(J-1)
12    CONTINUE
      INP(1)=ISAV
      CALL DES(INP,KEY,NEWKEY,0,JOT)
      RAN4=0.0
      DO 13 J=1,NACC
        IF(JOT(J).NE.0)RAN4=RAN4+POW(J)
13    CONTINUE
      RETURN
      END

ext: for


text:
      SUBROUTINE IB01OY( NS, NMAX, N, SV, INFO )
C
C     SLICOT RELEASE 5.5.
C
C     Copyright (c) 2002-2012 NICONET e.V.
C
C     PURPOSE
C
C     To ask for user's confirmation of the system order found by
C     SLICOT Library routine IB01OD. This routine may be modified,
C     but its interface must be preserved.
C
C     ARGUMENTS
C
C     Input/Output Parameters
C
C     NS      (input) INTEGER
C             The number of singular values.  NS > 0.
C
C     NMAX    (input) INTEGER
C             The maximum value of the system order.  0 <= NMAX <= NS.
C
C     N       (input/output) INTEGER
C             On entry, the estimate of the system order computed by
C             IB01OD routine.  0 <= N <= NS.
C             On exit, the user's estimate of the system order, which
C             could be identical with the input value of  N.
C             Note that the output value of  N  should be less than
C             or equal to  NMAX.
C
C     SV      (input) DOUBLE PRECISION array, dimension ( NS )
C             The singular values, in descending order, used for
C             determining the system order.
C
C     Error Indicator
C
C     INFO    INTEGER
C             = 0:  successful exit;
C             < 0:  if INFO = -i, the i-th argument had an illegal
C                   value.
C
C     CONTRIBUTORS
C
C     V. Sima, Research Institute for Informatics, Bucharest, Aug. 1999.
C
C     REVISIONS
C
C     -
C
C     KEYWORDS
C
C     Identification, parameter estimation, singular values, structure
C     identification.
C
C  *********************************************************************
C
C     .. Parameters ..
      INTEGER            INTRMN, OUTRMN
      PARAMETER          ( INTRMN = 5, OUTRMN = 6 )
C        INTRMN is the unit number for the (terminal) input device.
C        OUTRMN is the unit number for the (terminal) output device.
C     ..
C     .. Scalar Arguments ..
      INTEGER            INFO, N, NMAX, NS
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION   SV( * )
C     ..
C     .. Local Scalars ..
      LOGICAL            YES
      INTEGER            I
      CHARACTER          ANS
C     ..
C     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           LSAME
C     ..
C     .. External Subroutines ..
      EXTERNAL           XERBLA
C
C     .. Executable Statements ..
C
C     Check the scalar input parameters.
C
      INFO = 0
      IF( NS.LE.0 ) THEN
         INFO = -1
      ELSE IF( NMAX.LT.0 .OR. NMAX.GT.NS ) THEN
         INFO = -2
      ELSE IF( N.LT.0 .OR. N.GT.NS ) THEN
         INFO = -3
      END IF
C
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'IB01OY', -INFO )
         RETURN
      END IF
C
      WRITE( OUTRMN, '(/'' Singular values (in descending order) used'',
     $                  '' to estimate the system order:'', //
     $                     (5D15.8) )' ) ( SV(I), I = 1, NS )
      WRITE( OUTRMN, '(/'' Estimated order of the system,  n = '', I5 )'
     $     )               N
      WRITE( OUTRMN, '(/'' Do you want this value of  n  to be used'',
     $                  '' to determine the system matrices?'' )' )
C
   10 CONTINUE
         WRITE( OUTRMN, '(/''  Type "yes" or "no":  '' )' )
         READ ( INTRMN,  '( A )' ) ANS
         YES = LSAME( ANS, 'Y' )
         IF( YES ) THEN
            IF( N.LE.NMAX ) THEN
C
C              The value of n is adequate and has been confirmed.
C
               RETURN
            ELSE
C
C              The estimated value of n is not acceptable.
C
               WRITE( OUTRMN, '(/'' n  should be less than or equal'',
     $                           '' to '', I5 )' ) NMAX
               WRITE( OUTRMN, '( '' (It may be useful to restart'',
     $                           '' with a larger tolerance.)'' )' )
               GO TO 20
            END IF
C
         ELSE IF( LSAME( ANS, 'N' ) ) THEN
            GO TO 20
         ELSE
C
C           Wrong answer should be re-entered.
C
            GO TO 10
         END IF
C
C     Enter the desired value of n.
C
   20 CONTINUE
         WRITE( OUTRMN,'(/'' Enter the desired value of n (n <= '', I5,
     $                    '');  n = '' )' ) NMAX
         READ ( INTRMN, * ) N
         IF ( N.LT.0 ) THEN
C
C           The specified value of n is not acceptable.
C
            WRITE( OUTRMN, '(/'' n  should be larger than zero.'' )' )
            GO TO 20
         ELSE IF ( N.GT.NMAX ) THEN
C
C           The specified value of n is not acceptable.
C
            WRITE( OUTRMN, '(/'' n  should be less than or equal to '',
     $                   I5 )' ) NMAX
            GO TO 20
         END IF
C
      RETURN
C
C *** Last line of IB01OY ***
      END

max_stars_repo_name: bgin/MissileSimulation

max_issues_repo_name: bgin/MissileSimulation

max_forks_repo_name: bgin/MissileSimulation


text:
C    Copyright (c) 2014, Sandia Corporation.
C    Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
C    the U.S. Government retains certain rights in this software.
C
C    Redistribution and use in source and binary forms, with or without
C    modification, are permitted provided that the following conditions are
C    met:
C
C        * Redistributions of source code must retain the above copyright
C          notice, this list of conditions and the following disclaimer.
C
C        * Redistributions in binary form must reproduce the above
C          copyright notice, this list of conditions and the following
C          disclaimer in the documentation and/or other materials provided
C          with the distribution.
C
C        * Neither the name of Sandia Corporation nor the names of its
C          contributors may be used to endorse or promote products derived
C          from this software without specific prior written permission.
C
C    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C    A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C    OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C    SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C    LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C    DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C    THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
C

      program testcpnl

c
c This is a test program for the Fortran binding of the EXODUS II
c database copy function (excopy).
c
      implicit none

      include 'exodusII.inc'

      integer iin, iout, exoid, exoid1, ierr, cpu_ws, io_ws, mod_sz

      real vers

      data iin /5/, iout /6/

c
c open EXODUS II input file
c

c the setting of cpu_ws isn't used for copying but will test the
c conversion routines

      cpu_ws = 8
      io_ws = 4

      exoid = exopen ("test.exo", EXREAD, cpu_ws, io_ws, vers, ierr)
      write (iout, '(/"after exopen, error = ",i3)')
     1			ierr

      write (iout, '("test.exo is an EXODUSII file; version ",
     1                f4.2)') vers
      write (iout, '(" I/O word size: ",i4)') io_ws
      mod_sz = exlgmd(exoid)
      write (iout, '(" Model Size",i2)') mod_sz

c
c  create EXODUS II output file with default size reals
c
c the setting of cpu_ws isn't used for copying but will test the
c conversion routines

      cpu_ws = 8
      io_ws = 0

      exoid1 = excre ("testcp.exo",
     1               EXCLOB+EXNORM, cpu_ws, io_ws, ierr)
      write (iout,'("after excre, id = ", i6, ", error = ",i4)')
     1               exoid1, ierr
      write (iout,'(" I/O word size: ",i4)') io_ws

      mod_sz = exlgmd(exoid1)
      write (iout, '(" Model Size",i2)') mod_sz

      write (iout,'("after excre, error = ", i4)') ierr

      call excopy (exoid, exoid1, ierr)
      write (iout, '(/"after excopy, error = ", i3)' ) ierr

      call exclos (exoid, ierr)
      write (iout, '(/"after exclos, error = ", i3)' ) ierr

      call exclos (exoid1, ierr)
      write (iout, '(/"after exclos, error = ", i3)' ) ierr

      stop
      end

max_stars_repo_name: milljm/seacas

max_issues_repo_name: milljm/seacas

max_forks_repo_name: milljm/seacas


text:
! This source file is part of the Limited-area GAME version (L-GAME), which is released under the MIT license.
! Github repository: https://github.com/OpenNWP/L-GAME

module column_solvers

  ! This module contains the implicit vertical routines (implicit part of the HEVI scheme).

  use run_nml,        only: nlins,ncols,wp,nlays,dtime,p_0,toa,impl_weight,partial_impl_weight
  use definitions,    only: t_grid,t_state,t_tend
  use thermodynamics, only: spec_heat_cap_diagnostics_v,spec_heat_cap_diagnostics_p,gas_constant_diagnostics
  use diff_nml,       only: lklemp,klemp_damp_max,klemp_begin_rel

  implicit none

  private

  public :: three_band_solver_ver

  contains

  subroutine three_band_solver_ver(state_old,state_new,tend,grid,rk_step)

    type(t_state), intent(in)    :: state_old ! state at the old timestep
    type(t_state), intent(inout) :: state_new ! state at the new timestep
    type(t_tend),  intent(in)    :: tend      ! explicit tendencies
    type(t_grid),  intent(in)    :: grid      ! model grid
    integer,       intent(in)    :: rk_step   ! Runge Kutta substep

    ! local variables
    real(wp)                 :: c_vector(nlays-2)       ! needed for the vertical solver
    real(wp)                 :: d_vector(nlays-1)       ! needed for the vertical solver
    real(wp)                 :: e_vector(nlays-2)       ! needed for the vertical solver
    real(wp)                 :: r_vector(nlays-1)       ! needed for the vertical solver
    real(wp)                 :: solution(nlays-1)       ! covariant mass flux density at the interfaces (solution)
    real(wp)                 :: rho_expl(nlays)         ! explicit mass density
    real(wp)                 :: rhotheta_expl(nlays)    ! explicit potential temperature density
    real(wp)                 :: exner_pert_expl(nlays)  ! explicit Exner pressure perturbation
    real(wp)                 :: theta_pert_expl(nlays)  ! explicit potential temperature perturbation
    real(wp)                 :: rho_int_old(nlays-1)    ! old interface mass density
    real(wp)                 :: rho_int_expl(nlays-1)   ! explicit interface mass density
    real(wp)                 :: theta_int_new(nlays-1)  ! preliminary new potential temperature interface values
    integer                  :: ji,jk,jl                ! loop variables
    real(wp)                 :: rho_int_new             ! new density interface value
    real(wp)                 :: alpha_old(nlays)        ! alpha at the old time step
    real(wp)                 :: beta_old(nlays)         ! beta at the old time step
    real(wp)                 :: gamma_old(nlays)        ! gamma at the old time step
    real(wp)                 :: alpha_new(nlays)        ! alpha at the new time step
    real(wp)                 :: beta_new(nlays)         ! beta at the new time step
    real(wp)                 :: gamma_new(nlays)        ! gamma at the new time step
    real(wp)                 :: alpha(nlays)            ! alpha
    real(wp)                 :: beta(nlays)             ! beta
    real(wp)                 :: gammaa(nlays)           ! gamma
    real(wp)                 :: c_v                     ! specific heat capacity at constant volume
    real(wp)                 :: c_p                     ! specific heat capacity at constant pressure
    real(wp)                 :: r_d                     ! individual gas constant of dry air
    real(wp)                 :: damping_start_height    ! lower boundary height of the Klemp layer
    real(wp)                 :: damping_coeff           ! damping coefficient of the Klemp layer
    real(wp)                 :: z_above_damping         ! height above the lower boundary of the damping height

    c_v = spec_heat_cap_diagnostics_v(1)
    c_p = spec_heat_cap_diagnostics_p(1)
    r_d = gas_constant_diagnostics(1)
    damping_start_height = klemp_begin_rel*toa

    !$OMP PARALLEL
    !$OMP DO PRIVATE(ji,jk,jl,c_vector,d_vector,e_vector,r_vector,solution, &
    !$OMP rho_expl,rhotheta_expl,exner_pert_expl,theta_pert_expl,rho_int_old, &
    !$OMP rho_int_expl,theta_int_new,rho_int_new,alpha_old,beta_old,gamma_old,alpha_new, &
    !$OMP beta_new,gamma_new,alpha,beta,gammaa,damping_coeff,z_above_damping)
    do ji=1,nlins
      do jk=1,ncols

        ! explicit quantities
        do jl=1,nlays
          ! explicit density
          rho_expl(jl) = state_old%rho(ji+1,jk+1,jl) + dtime*tend%rho(ji,jk,jl)
          ! explicit potential temperature density
          rhotheta_expl(jl) = state_old%rhotheta(ji+1,jk+1,jl) + dtime*tend%rhotheta(ji,jk,jl)
          if (rk_Step == 1) then
            ! old time step partial derivatives of theta and Pi (divided by the volume)
            alpha(jl) = -state_old%rhotheta(ji+1,jk+1,jl)/state_old%rho(ji+1,jk+1,jl)**2/grid%volume(ji,jk,jl)
            beta(jl)  = 1._wp/state_old%rho(ji+1,jk+1,jl)/grid%volume(ji,jk,jl)
            gammaa(jl) = r_d/(c_v*state_old%rhotheta(ji+1,jk+1,jl))* &
            (grid%exner_bg(ji+1,jk+1,jl)+state_old%exner_pert(ji+1,jk+1,jl))/grid%volume(ji,jk,jl)
          else
            ! old time step partial derivatives of theta and Pi
            alpha_old(jl) = -state_old%rhotheta(ji+1,jk+1,jl)/state_old%rho(ji+1,jk+1,jl)**2
            beta_old(jl)  = 1._wp/state_old%rho(ji+1,jk+1,jl)
            gamma_old(jl) = r_d/(c_v*state_old%rhotheta(ji+1,jk+1,jl))* &
            (grid%exner_bg(ji+1,jk+1,jl)+state_old%exner_pert(ji+1,jk+1,jl))
            ! new time step partial derivatives of theta and Pi
            alpha_new(jl) = -state_new%rhotheta(ji+1,jk+1,jl)/state_new%rho(ji+1,jk+1,jl)**2
            beta_new(jl)  = 1._wp/state_new%rho(ji+1,jk+1,jl)
            gamma_new(jl) = r_d/(c_v*state_new%rhotheta(ji+1,jk+1,jl)) &
            *(grid%exner_bg(ji+1,jk+1,jl)+state_new%exner_pert(ji+1,jk+1,jl))
            ! interpolation of partial derivatives of theta and Pi (divided by the volume)
            alpha (jl) = ((1._wp - partial_impl_weight)*alpha_old(jl) + partial_impl_weight*alpha_new(jl))/grid%volume(ji,jk,jl)
            beta  (jl) = ((1._wp - partial_impl_weight)*beta_old (jl) + partial_impl_weight*beta_new (jl))/grid%volume(ji,jk,jl)
            gammaa(jl) = ((1._wp - partial_impl_weight)*gamma_old(jl) + partial_impl_weight*gamma_new(jl))/grid%volume(ji,jk,jl)
          endif
          ! explicit potential temperature perturbation
          theta_pert_expl(jl) = state_old%theta_pert(ji+1,jk+1,jl) + dtime*grid%volume(ji,jk,jl)*(alpha(jl)*tend%rho(ji,jk,jl) &
          + beta(jl)*tend%rhotheta(ji,jk,jl))
          ! explicit Exner pressure perturbation
          exner_pert_expl(jl) = state_old%exner_pert(ji+1,jk+1,jl) + dtime*grid%volume(ji,jk,jl)*gammaa(jl)*tend%rhotheta(ji,jk,jl)
        enddo

        ! interface values
        do jl=1,nlays-1
          rho_int_old(jl) = 0.5_wp*(state_old%rho(ji+1,jk+1,jl)+state_old%rho(ji+1,jk+1,jl+1))
          rho_int_expl(jl) = 0.5_wp*(rho_expl(jl)+rho_expl(jl+1))
          theta_int_new(jl) = 0.5_wp*(state_new%rhotheta(ji+1,jk+1,jl)/state_new%rho(ji+1,jk+1,jl) &
          + state_new%rhotheta(ji+1,jk+1,jl+1)/state_new%rho(ji+1,jk+1,jl+1))
        enddo

        ! filling up the coefficient vectors
        do jl=1,nlays-1
          ! main diagonal
          d_vector(jl) = -theta_int_new(jl)**2*(gammaa(jl)+gammaa(jl+1)) &
          + 0.5_wp*(grid%exner_bg(ji+1,jk+1,jl)-grid%exner_bg(ji+1,jk+1,jl+1)) &
          *(alpha(jl+1)-alpha(jl)+theta_int_new(jl)*(beta(jl+1)-beta(jl))) &
          - (grid%z_geo_scal(ji+1,jk+1,jl)-grid%z_geo_scal(ji+1,jk+1,jl+1))/(impl_weight*dtime**2*c_p*rho_int_old(jl)) &
          *(2._wp/grid%area_z(ji,jk,jl+1)+dtime*state_old%wind_w(ji+1,jk+1,jl+1)*0.5_wp &
          *(-1._wp/grid%volume(ji,jk,jl)+1._wp/grid%volume(ji,jk,jl+1)))
          ! right hand side
          r_vector(jl) = -(state_old%wind_w(ji+1,jk+1,jl+1)+dtime*tend%wind_w(ji,jk,jl+1))* &
          (grid%z_geo_scal(ji+1,jk+1,jl)-grid%z_geo_scal(ji+1,jk+1,jl+1)) &
          /(impl_weight*dtime**2*c_p) &
          + theta_int_new(jl)*(exner_pert_expl(jl)-exner_pert_expl(jl+1))/dtime &
          + 0.5_wp/dtime*(theta_pert_expl(jl)+theta_pert_expl(jl+1))*(grid%exner_bg(ji+1,jk+1,jl)-grid%exner_bg(ji+1,jk+1,jl+1)) &
          - (grid%z_geo_scal(ji+1,jk+1,jl)-grid%z_geo_scal(ji+1,jk+1,jl+1))/(impl_weight*dtime**2*c_p) &
          *state_old%wind_w(ji+1,jk+1,jl+1)*rho_int_expl(jl)/rho_int_old(jl)
        enddo

        do jl=1,nlays-2
          ! lower diagonal
          c_vector(jl) = theta_int_new(jl+1)*gammaa(jl+1)*theta_int_new(jl) &
          + 0.5_wp*(grid%exner_bg(ji+1,jk+1,jl+1)-grid%exner_bg(ji+1,jk+1,jl+2)) &
          *(alpha(jl+1)+beta(jl+1)*theta_int_new(jl)) &
          - (grid%z_geo_scal(ji+1,jk+1,jl+1)-grid%z_geo_scal(ji+1,jk+1,jl+2))/(impl_weight*dtime*c_p)*0.5_wp &
          *state_old%wind_w(ji+1,jk+1,jl+2)/(grid%volume(ji,jk,jl+1)*rho_int_old(jl+1))
          ! upper diagonal
          e_vector(jl) = theta_int_new(jl)*gammaa(jl+1)*theta_int_new(jl+1) &
          - 0.5_wp*(grid%exner_bg(ji+1,jk+1,jl)-grid%exner_bg(ji+1,jk+1,jl+1)) &
          *(alpha(jl+1)+beta(jl+1)*theta_int_new(jl+1)) &
          + (grid%z_geo_scal(ji+1,jk+1,jl)-grid%z_geo_scal(ji+1,jk+1,jl+1))/(impl_weight*dtime*c_p)*0.5_wp &
          *state_old%wind_w(ji+1,jk+1,jl+1)/(grid%volume(ji,jk,jl+1)*rho_int_old(jl))
        enddo

        call thomas_algorithm(c_vector,d_vector,e_vector,r_vector,solution,nlays-1)

        ! Klemp swamp layer
        do jl=1,nlays-1
          z_above_damping = grid%z_geo_w(ji+1,jk+1,jl+1)-damping_start_height
          if (z_above_damping < 0._wp .or. .not. lklemp) then
            damping_coeff = 0._wp
          else
            damping_coeff = klemp_damp_max*sin(0.5_wp*4*atan(1.d0)*z_above_damping/(toa-damping_start_height))**2
          endif
          solution(jl) = solution(jl)/(1._wp + damping_coeff*dtime)
        enddo

        ! results
        ! density, potential temperature density
        do jl=2,nlays-1
          state_new%rho(ji+1,jk+1,jl) = rho_expl(jl) + dtime*(-solution(jl-1)+solution(jl))/grid%volume(ji,jk,jl)
          state_new%rhotheta(ji+1,jk+1,jl) = rhotheta_expl(jl) &
          + dtime*(-theta_int_new(jl-1)*solution(jl-1)+theta_int_new(jl)*solution(jl))/grid%volume(ji,jk,jl)
        enddo
        ! uppermost layer
        state_new%rho(ji+1,jk+1,1) = rho_expl(1) + dtime*solution(1)/grid%volume(ji,jk,1)
        state_new%rhotheta(ji+1,jk+1,1) = rhotheta_expl(1)+dtime*theta_int_new(1)*solution(1)/grid%volume(ji,jk,1)
        ! lowest layer
        state_new%rho(ji+1,jk+1,nlays) = rho_expl(nlays) - dtime*solution(nlays-1)/grid%volume(ji,jk,nlays)
        state_new%rhotheta(ji+1,jk+1,nlays) = rhotheta_expl(nlays) &
        -dtime*theta_int_new(nlays-1)*solution(nlays-1)/grid%volume(ji,jk,nlays)
        ! vertical velocity
        do jl=2,nlays
          rho_int_new = 0.5_wp*(state_new%rho(ji+1,jk+1,jl-1)+state_new%rho(ji+1,jk+1,jl))
          state_new%wind_w(ji+1,jk+1,jl) = (2._wp*solution(jl-1)/grid%area_z(ji,jk,jl) &
          - rho_int_new*state_old%wind_w(ji+1,jk+1,jl))/rho_int_old(jl-1)
        enddo
        ! Exner pressure
        do jl=1,nlays
          state_new%exner_pert(ji+1,jk+1,jl) = state_old%exner_pert(ji+1,jk+1,jl) &
          + grid%volume(ji,jk,jl)*gammaa(jl)*(state_new%rhotheta(ji+1,jk+1,jl)-state_old%rhotheta(ji+1,jk+1,jl))
        enddo

      enddo
    enddo
    !$OMP END DO
    !$OMP END PARALLEL

    ! potential temperature perturbation at the new time step
    state_new%theta_pert(:,:,:) = state_new%rhotheta(:,:,:)/state_new%rho(:,:,:) - grid%theta_bg(:,:,:)

  end subroutine three_band_solver_ver

  subroutine thomas_algorithm(c_vector,d_vector,e_vector,r_vector,solution_vector,solution_length)

    ! This subroutine solves a system of linear equations with a three-band matrix.

    real(wp), intent(in)    :: c_vector(:)
    real(wp), intent(in)    :: d_vector(:)
    real(wp), intent(in)    :: e_vector(:)
    real(wp), intent(in)    :: r_vector(:)
    real(wp), intent(inout) :: solution_vector(:) ! vector containing the solution
    integer,  intent(in)    :: solution_length    ! length of the solution vector

    ! local variables
    real(wp) :: e_prime_vector(solution_length-1) ! help vector for solving the matrix equation
    real(wp) :: r_prime_vector(solution_length)   ! help vector for solving the matrix equation
    integer  :: jl                                ! loop index

    ! downward sweep (matrix)
    e_prime_vector(1) = e_vector(1)/d_vector(1)
    do jl=2,solution_length-1
      e_prime_vector(jl) = e_vector(jl)/(d_vector(jl) - e_prime_vector(jl-1)*c_vector(jl-1))
    enddo
    ! downward sweep (right-hand side)
    r_prime_vector(1) = r_vector(1)/d_vector(1)
    do jl=2,solution_length
      r_prime_vector(jl) = (r_vector(jl) - r_prime_vector(jl-1)*c_vector(jl-1)) &
      /(d_vector(jl) - e_prime_vector(jl-1)*c_vector(jl-1))
    enddo

    ! upward sweep (final solution)
    solution_vector(solution_length) = r_prime_vector(solution_length)
    do jl=solution_length-1,1,-1
      solution_vector(jl) = r_prime_vector(jl) - e_prime_vector(jl)*solution_vector(jl+1)
    enddo

  end subroutine thomas_algorithm

end module column_solvers

max_stars_repo_name: OpenNWP/L-GAME

max_issues_repo_name: OpenNWP/L-GAME

max_forks_repo_name: OpenNWP/L-GAME