text:
! Created by EverLookNeverSee@GitHub on 6/8/20
! For more information see FCS/img/Exercise_09_b.png

program main
    implicit none
    ! declaring and initializing variables
    integer :: i, j, m, n
    real :: norm, tmp, sum_of_elements = 0.0
    real, allocatable, dimension(:,:) :: A
    ! getting shape of matrix from user input
    print *, "Enter the shape of matrix(rows, columns):"
    read *, m, n
    ! allocating mem space to matrix
    allocate(A(m, n))
    ! getting matrix elements from user input
    do i = 1, m
        do j = 1, n
            print *, "A(", i, j, "):"
            read *, A(i, j)
        end do
    end do
    ! sum of all squared elements in matrix
    do i = 1, m
        do j = 1, n
            tmp = (abs(A(i, j))) ** 2
            sum_of_elements = sum_of_elements + tmp
        end do
    end do
    ! square root of sum
    norm = sqrt(sum_of_elements)
    ! printing results
    print *, "Norm =", norm
end program main

max_stars_repo_path: src/numerical methods/Exercise_09_b.f90

max_stars_repo_name: EverLookNeverSee/FCS

max_issues_repo_path: src/numerical methods/Exercise_09_b.f90

max_issues_repo_name: EverLookNeverSee/FCS

max_forks_repo_path: src/numerical methods/Exercise_09_b.f90

max_forks_repo_name: EverLookNeverSee/FCS


text:
*> \brief \b SERRGT
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*  Definition:
*  ===========
*
*       SUBROUTINE SERRGT( PATH, NUNIT )
*
*       .. Scalar Arguments ..
*       CHARACTER*3        PATH
*       INTEGER            NUNIT
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> SERRGT tests the error exits for the REAL tridiagonal
*> routines.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] PATH
*> \verbatim
*>          PATH is CHARACTER*3
*>          The LAPACK path name for the routines to be tested.
*> \endverbatim
*>
*> \param[in] NUNIT
*> \verbatim
*>          NUNIT is INTEGER
*>          The unit number for output.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup single_lin
*
*  =====================================================================
      SUBROUTINE SERRGT( PATH, NUNIT )
*
*  -- LAPACK test routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      CHARACTER*3        PATH
      INTEGER            NUNIT
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            NMAX
      PARAMETER          ( NMAX = 2 )
*     ..
*     .. Local Scalars ..
      CHARACTER*2        C2
      INTEGER            INFO
      REAL               ANORM, RCOND
*     ..
*     .. Local Arrays ..
      INTEGER            IP( NMAX ), IW( NMAX )
      REAL               B( NMAX ), C( NMAX ), CF( NMAX ), D( NMAX ),
     $                   DF( NMAX ), E( NMAX ), EF( NMAX ), F( NMAX ),
     $                   R1( NMAX ), R2( NMAX ), W( NMAX ), X( NMAX )
*     ..
*     .. External Functions ..
      LOGICAL            LSAMEN
      EXTERNAL           LSAMEN
*     ..
*     .. External Subroutines ..
      EXTERNAL           ALAESM, CHKXER, SGTCON, SGTRFS, SGTTRF, SGTTRS,
     $                   SPTCON, SPTRFS, SPTTRF, SPTTRS
*     ..
*     .. Scalars in Common ..
      LOGICAL            LERR, OK
      CHARACTER*32       SRNAMT
      INTEGER            INFOT, NOUT
*     ..
*     .. Common blocks ..
      COMMON             / INFOC / INFOT, NOUT, OK, LERR
      COMMON             / SRNAMC / SRNAMT
*     ..
*     .. Executable Statements ..
*
      NOUT = NUNIT
      WRITE( NOUT, FMT = * )
      C2 = PATH( 2: 3 )
      D( 1 ) = 1.
      D( 2 ) = 2.
      DF( 1 ) = 1.
      DF( 2 ) = 2.
      E( 1 ) = 3.
      E( 2 ) = 4.
      EF( 1 ) = 3.
      EF( 2 ) = 4.
      ANORM = 1.0
      OK = .TRUE.
*
      IF( LSAMEN( 2, C2, 'GT' ) ) THEN
*
*        Test error exits for the general tridiagonal routines.
*
*        SGTTRF
*
         SRNAMT = 'SGTTRF'
         INFOT = 1
         CALL SGTTRF( -1, C, D, E, F, IP, INFO )
         CALL CHKXER( 'SGTTRF', INFOT, NOUT, LERR, OK )
*
*        SGTTRS
*
         SRNAMT = 'SGTTRS'
         INFOT = 1
         CALL SGTTRS( '/', 0, 0, C, D, E, F, IP, X, 1, INFO )
         CALL CHKXER( 'SGTTRS', INFOT, NOUT, LERR, OK )
         INFOT = 2
         CALL SGTTRS( 'N', -1, 0, C, D, E, F, IP, X, 1, INFO )
         CALL CHKXER( 'SGTTRS', INFOT, NOUT, LERR, OK )
         INFOT = 3
         CALL SGTTRS( 'N', 0, -1, C, D, E, F, IP, X, 1, INFO )
         CALL CHKXER( 'SGTTRS', INFOT, NOUT, LERR, OK )
         INFOT = 10
         CALL SGTTRS( 'N', 2, 1, C, D, E, F, IP, X, 1, INFO )
         CALL CHKXER( 'SGTTRS', INFOT, NOUT, LERR, OK )
*
*        SGTRFS
*
         SRNAMT = 'SGTRFS'
         INFOT = 1
         CALL SGTRFS( '/', 0, 0, C, D, E, CF, DF, EF, F, IP, B, 1, X, 1,
     $                R1, R2, W, IW, INFO )
         CALL CHKXER( 'SGTRFS', INFOT, NOUT, LERR, OK )
         INFOT = 2
         CALL SGTRFS( 'N', -1, 0, C, D, E, CF, DF, EF, F, IP, B, 1, X,
     $                1, R1, R2, W, IW, INFO )
         CALL CHKXER( 'SGTRFS', INFOT, NOUT, LERR, OK )
         INFOT = 3
         CALL SGTRFS( 'N', 0, -1, C, D, E, CF, DF, EF, F, IP, B, 1, X,
     $                1, R1, R2, W, IW, INFO )
         CALL CHKXER( 'SGTRFS', INFOT, NOUT, LERR, OK )
         INFOT = 13
         CALL SGTRFS( 'N', 2, 1, C, D, E, CF, DF, EF, F, IP, B, 1, X, 2,
     $                R1, R2, W, IW, INFO )
         CALL CHKXER( 'SGTRFS', INFOT, NOUT, LERR, OK )
         INFOT = 15
         CALL SGTRFS( 'N', 2, 1, C, D, E, CF, DF, EF, F, IP, B, 2, X, 1,
     $                R1, R2, W, IW, INFO )
         CALL CHKXER( 'SGTRFS', INFOT, NOUT, LERR, OK )
*
*        SGTCON
*
         SRNAMT = 'SGTCON'
         INFOT = 1
         CALL SGTCON( '/', 0, C, D, E, F, IP, ANORM, RCOND, W, IW,
     $                INFO )
         CALL CHKXER( 'SGTCON', INFOT, NOUT, LERR, OK )
         INFOT = 2
         CALL SGTCON( 'I', -1, C, D, E, F, IP, ANORM, RCOND, W, IW,
     $                INFO )
         CALL CHKXER( 'SGTCON', INFOT, NOUT, LERR, OK )
         INFOT = 8
         CALL SGTCON( 'I', 0, C, D, E, F, IP, -ANORM, RCOND, W, IW,
     $                INFO )
         CALL CHKXER( 'SGTCON', INFOT, NOUT, LERR, OK )
*
      ELSE IF( LSAMEN( 2, C2, 'PT' ) ) THEN
*
*        Test error exits for the positive definite tridiagonal
*        routines.
*
*        SPTTRF
*
         SRNAMT = 'SPTTRF'
         INFOT = 1
         CALL SPTTRF( -1, D, E, INFO )
         CALL CHKXER( 'SPTTRF', INFOT, NOUT, LERR, OK )
*
*        SPTTRS
*
         SRNAMT = 'SPTTRS'
         INFOT = 1
         CALL SPTTRS( -1, 0, D, E, X, 1, INFO )
         CALL CHKXER( 'SPTTRS', INFOT, NOUT, LERR, OK )
         INFOT = 2
         CALL SPTTRS( 0, -1, D, E, X, 1, INFO )
         CALL CHKXER( 'SPTTRS', INFOT, NOUT, LERR, OK )
         INFOT = 6
         CALL SPTTRS( 2, 1, D, E, X, 1, INFO )
         CALL CHKXER( 'SPTTRS', INFOT, NOUT, LERR, OK )
*
*        SPTRFS
*
         SRNAMT = 'SPTRFS'
         INFOT = 1
         CALL SPTRFS( -1, 0, D, E, DF, EF, B, 1, X, 1, R1, R2, W, INFO )
         CALL CHKXER( 'SPTRFS', INFOT, NOUT, LERR, OK )
         INFOT = 2
         CALL SPTRFS( 0, -1, D, E, DF, EF, B, 1, X, 1, R1, R2, W, INFO )
         CALL CHKXER( 'SPTRFS', INFOT, NOUT, LERR, OK )
         INFOT = 8
         CALL SPTRFS( 2, 1, D, E, DF, EF, B, 1, X, 2, R1, R2, W, INFO )
         CALL CHKXER( 'SPTRFS', INFOT, NOUT, LERR, OK )
         INFOT = 10
         CALL SPTRFS( 2, 1, D, E, DF, EF, B, 2, X, 1, R1, R2, W, INFO )
         CALL CHKXER( 'SPTRFS', INFOT, NOUT, LERR, OK )
*
*        SPTCON
*
         SRNAMT = 'SPTCON'
         INFOT = 1
         CALL SPTCON( -1, D, E, ANORM, RCOND, W, INFO )
         CALL CHKXER( 'SPTCON', INFOT, NOUT, LERR, OK )
         INFOT = 4
         CALL SPTCON( 0, D, E, -ANORM, RCOND, W, INFO )
         CALL CHKXER( 'SPTCON', INFOT, NOUT, LERR, OK )
      END IF
*
*     Print a summary line.
*
      CALL ALAESM( PATH, OK, NOUT )
*
      RETURN
*
*     End of SERRGT
*
      END

ext: f

max_stars_repo_name: sergey-v-kuznetsov/lapack

max_issues_repo_name: sergey-v-kuznetsov/lapack

max_forks_repo_name: sergey-v-kuznetsov/lapack


text:
program dim_assgn_test
! Warning:
integer, private, dimension(2,2) :: a(2,3)
! Error:
x = 5
end program

max_stars_repo_name: Thirumalai-Shaktivel/lfortran

max_issues_repo_name: Thirumalai-Shaktivel/lfortran

max_forks_repo_name: Thirumalai-Shaktivel/lfortran


text:
TYPE,BIND(C) :: GDALDatasetH
  PRIVATE
  TYPE(c_ptr) :: ptr = C_NULL_PTR
END TYPE GDALDatasetH

TYPE,BIND(C) :: GDALMajorObjectH
  PRIVATE
  TYPE(c_ptr) :: ptr = C_NULL_PTR
END TYPE GDALMajorObjectH

TYPE,BIND(C) :: GDALDriverH
  PRIVATE
  TYPE(c_ptr) :: ptr = C_NULL_PTR
END TYPE GDALDriverH

TYPE,BIND(C) :: GDALRasterAttributeTableH
  PRIVATE
  TYPE(c_ptr) :: ptr = C_NULL_PTR
END TYPE GDALRasterAttributeTableH

TYPE,BIND(C) :: GDALColorTableH
  PRIVATE
  TYPE(c_ptr) :: ptr = C_NULL_PTR
END TYPE GDALColorTableH

TYPE,BIND(C) :: GDALRasterBandH
  PRIVATE
  TYPE(c_ptr) :: ptr = C_NULL_PTR
END TYPE GDALRasterBandH

INTERFACE gdalassociated
  MODULE PROCEDURE gdalassociated_gdaldataseth
  MODULE PROCEDURE gdalassociated_gdalmajorobjecth
  MODULE PROCEDURE gdalassociated_gdaldriverh
  MODULE PROCEDURE gdalassociated_gdalrasterattributetableh
  MODULE PROCEDURE gdalassociated_gdalcolortableh
  MODULE PROCEDURE gdalassociated_gdalrasterbandh
END INTERFACE
INTERFACE gdalnullify
  MODULE PROCEDURE gdalnullify_gdaldataseth
  MODULE PROCEDURE gdalnullify_gdalmajorobjecth
  MODULE PROCEDURE gdalnullify_gdaldriverh
  MODULE PROCEDURE gdalnullify_gdalrasterattributetableh
  MODULE PROCEDURE gdalnullify_gdalcolortableh
  MODULE PROCEDURE gdalnullify_gdalrasterbandh
END INTERFACE

max_stars_repo_name: gronki/slezanbear

max_issues_repo_name: gronki/slezanbear

max_forks_repo_name: gronki/slezanbear


text:
c-----------------------------------------------------------------
      subroutine create_initial_grid
c------------------------------------------------------------------

      include 'header.h'

      integer i

      nelt=1
      ntot=nelt*lx1*lx1*lx1
      tree(1)=1
      mt_to_id(1)=1
      do i=1,7,2
        xc(i,1)=0.d0
        xc(i+1,1)=1.d0
      end do

      do i=1,2
        yc(i,1)=0.d0
        yc(2+i,1)=1.d0
        yc(4+i,1)=0.d0
        yc(6+i,1)=1.d0
      end do

      do i=1,4
        zc(i,1)=0.d0
        zc(4+i,1)=1.d0
      end do

      do i=1,6
        cbc(i,1)=0
      end do

      return

      end

c-----------------------------------------------------------------
      subroutine coef
c-----------------------------------------------------------------
c
c     generate
c
c            - collocation points
c            - weights
c            - derivative matrices
c            - projection matrices
c            - interpolation matrices
c
c     associated with the
c
c            - gauss-legendre lobatto mesh (suffix m1)
c
c----------------------------------------------------------------

      include 'header.h'

      integer i,j,k

c.....for gauss-legendre lobatto mesh (suffix m1)
c.....generate collocation points and weights

      zgm1(1)=-1.d0
      zgm1(2)=-0.6546536707079771d0
      zgm1(3)=0.d0
      zgm1(4)= 0.6546536707079771d0
      zgm1(5)=1.d0
      wxm1(1)=0.1d0
      wxm1(2)=49.d0/90.d0
      wxm1(3)=32.d0/45.d0
      wxm1(4)=wxm1(2)
      wxm1(5)=0.1d0

      do k=1,lx1
        do j=1,lx1
          do i=1,lx1
            w3m1(i,j,k)=wxm1(i)*wxm1(j)*wxm1(k)
          end do
        end do
      end do

c.....generate derivative matrices

      dxm1(1,1)=-5.0d0
      dxm1(2,1)=-1.240990253030982d0
      dxm1(3,1)= 0.375d0
      dxm1(4,1)=-0.2590097469690172d0
      dxm1(5,1)= 0.5d0
      dxm1(1,2)= 6.756502488724238d0
      dxm1(2,2)= 0.d0
      dxm1(3,2)=-1.336584577695453d0
      dxm1(4,2)= 0.7637626158259734d0
      dxm1(5,2)=-1.410164177942427d0
      dxm1(1,3)=-2.666666666666667d0
      dxm1(2,3)= 1.745743121887939d0
      dxm1(3,3)= 0.d0
      dxm1(4,3)=-dxm1(2,3)
      dxm1(5,3)=-dxm1(1,3)
      do j=4,lx1
        do i=1,lx1
          dxm1(i,j)=-dxm1(lx1+1-i,lx1+1-j)
        end do
      end do
      do j=1,lx1
        do i=1,lx1
          dxtm1(i,j)=dxm1(j,i)
        end do
      end do

c.....generate projection (mapping) matrices

      qbnew(1,1,1)=-0.1772843218615690d0
      qbnew(2,1,1)=9.375d-02
      qbnew(3,1,1)=-3.700139242414530d-02
      qbnew(1,2,1)= 0.7152146412463197d0
      qbnew(2,2,1)=-0.2285757930375471d0
      qbnew(3,2,1)= 8.333333333333333d-02
      qbnew(1,3,1)= 0.4398680650316104d0
      qbnew(2,3,1)= 0.2083333333333333d0
      qbnew(3,3,1)=-5.891568407922938d-02
      qbnew(1,4,1)= 8.333333333333333d-02
      qbnew(2,4,1)= 0.3561799597042137d0
      qbnew(3,4,1)=-4.854797457965334d-02
      qbnew(1,5,1)= 0.d0
      qbnew(2,5,1)=7.03125d-02
      qbnew(3,5,1)=0.d0

      do j=1,lx1
        do i=1,3
          qbnew(i,j,2)=qbnew(4-i,lx1+1-j,1)
        end do
      end do

c.....generate interpolation matrices for mesh refinement

      ixtmc1(1,1)=1.d0
      ixtmc1(2,1)=0.d0
      ixtmc1(3,1)=0.d0
      ixtmc1(4,1)=0.d0
      ixtmc1(5,1)=0.d0
      ixtmc1(1,2)= 0.3385078435248143d0
      ixtmc1(2,2)= 0.7898516348912331d0
      ixtmc1(3,2)=-0.1884018684471238d0
      ixtmc1(4,2)= 9.202967302175333d-02
      ixtmc1(5,2)=-3.198728299067715d-02
      ixtmc1(1,3)=-0.1171875d0
      ixtmc1(2,3)= 0.8840317166357952d0
      ixtmc1(3,3)= 0.3125d0
      ixtmc1(4,3)=-0.118406716635795d0
      ixtmc1(5,3)= 0.0390625d0
      ixtmc1(1,4)=-7.065070066767144d-02
      ixtmc1(2,4)= 0.2829703269782467d0
      ixtmc1(3,4)= 0.902687582732838d0
      ixtmc1(4,4)=-0.1648516348912333d0
      ixtmc1(5,4)= 4.984442584781999d-02
      ixtmc1(1,5)=0.d0
      ixtmc1(2,5)=0.d0
      ixtmc1(3,5)=1.d0
      ixtmc1(4,5)=0.d0
      ixtmc1(5,5)=0.d0
      do j=1,lx1
        do i=1,lx1
          ixmc1(i,j)=ixtmc1(j,i)
        end do
      end do

      do j=1,lx1
        do i=1,lx1
          ixtmc2(i,j)=ixtmc1(lx1+1-i,lx1+1-j)
        end do
      end do

      do j=1,lx1
        do i=1,lx1
          ixmc2(i,j)=ixtmc2(j,i)
        end do
      end do

c.....solution interpolation matrix for mesh coarsening

      map2(1)=-0.1179652785083428d0
      map2(2)= 0.5505046330389332d0
      map2(3)= 0.7024534364259963d0
      map2(4)=-0.1972224518285866d0
      map2(5)= 6.222966087199998d-02

      do i=1,lx1
        map4(i)=map2(lx1+1-i)
      end do

      return
      end

c-------------------------------------------------------------------
      subroutine geom1
c-------------------------------------------------------------------
c
c     routine to generate elemental geometry information on mesh m1,
c     (gauss-legendre lobatto mesh).
c
c         xrm1_s   -   dx/dr, dy/dr, dz/dr
c         rxm1_s   -   dr/dx, dr/dy, dr/dz
c         g1m1_s  geometric factors used in preconditioner computation
c         g4m1_s  g5m1_s  g6m1_s :
c         geometric factors used in lapacian opertor
c         jacm1    -   jacobian
c         bm1      -   mass matrix
c         xfrac    -   will be used in prepwork for calculating collocation
c                      coordinates
c         idel     -   collocation points index on element boundaries
c------------------------------------------------------------------

      include 'header.h'

      double precision temp,temp1,temp2,dtemp
      integer isize,i,j,k,ntemp,iel

      do i=1,lx1
        xfrac(i)=zgm1(i)*0.5d0 + 0.5d0
      end do

      do isize=1,refine_max
        temp=2.d0**(-isize-1)
        dtemp=1.d0/temp
        temp1=temp**3
        temp2=temp**2
        do k=1,lx1
          do j=1,lx1
            do i=1,lx1
              xrm1_s(i,j,k,isize)=dtemp
              jacm1_s(i,j,k,isize)=temp1
              rxm1_s(i,j,k,isize)=temp2
              g1m1_s(i,j,k,isize)=w3m1(i,j,k)*temp
              bm1_s(i,j,k,isize)=w3m1(i,j,k)*temp1
              g4m1_s(i,j,k,isize)=g1m1_s(i,j,k,isize)/wxm1(i)
              g5m1_s(i,j,k,isize)=g1m1_s(i,j,k,isize)/wxm1(j)
              g6m1_s(i,j,k,isize)=g1m1_s(i,j,k,isize)/wxm1(k)
            end do
          end do
        end do
      end do

      do iel = 1, lelt
        ntemp=lx1*lx1*lx1*(iel-1)
        do j = 1, lx1
          do i = 1, lx1
            idel(i,j,1,iel)=ntemp+(i-1)*lx1 + (j-1)*lx1*lx1+lx1
            idel(i,j,2,iel)=ntemp+(i-1)*lx1 + (j-1)*lx1*lx1+1
            idel(i,j,3,iel)=ntemp+(i-1)*1 + (j-1)*lx1*lx1+lx1*(lx1-1)+1
            idel(i,j,4,iel)=ntemp+(i-1)*1 + (j-1)*lx1*lx1+1
            idel(i,j,5,iel)=ntemp+(i-1)*1 + (j-1)*lx1+lx1*lx1*(lx1-1)+1
            idel(i,j,6,iel)=ntemp+(i-1)*1 + (j-1)*lx1+1
          end do
        end do
      end do

      return
      end

c------------------------------------------------------------------
      subroutine setdef
c------------------------------------------------------------------
c     compute the discrete laplacian operators
c------------------------------------------------------------------

      include 'header.h'

      integer i,j,ip

      call r_init(wdtdr(1,1),lx1*lx1,0.d0)

      do i=1,lx1
        do j=1,lx1
          do ip=1,lx1
            wdtdr(i,j) = wdtdr(i,j) + wxm1(ip)*dxm1(ip,i)*dxm1(ip,j)
          end do
        end do
      end do

      return
      end

c------------------------------------------------------------------
      subroutine prepwork
c------------------------------------------------------------------
c     mesh information preparations: calculate refinement levels of
c     each element, mask matrix for domain boundary and element
c     boundaries
c------------------------------------------------------------------

      include 'header.h'

      integer i, j, k, iel, iface, cb
      double precision rdlog2

      ntot = nelt*nxyz
      rdlog2 = 1.d0/dlog(2.d0)

c.....calculate the refinement levels of each element

      do iel = 1, nelt
        size_e(iel)=-dlog(xc(2,iel)-xc(1,iel))*rdlog2+1.d-8
      end do

c.....mask matrix for element boundary

      do iel = 1, nelt
        call r_init(tmult(1,1,1,iel),nxyz,1.d0)
        do iface=1,nsides
          call facev(tmult(1,1,1,iel),iface,0.0d0)
        end do
      end do

c.....masks for domain boundary at mortar

      call r_init(tmmor,nmor,1.d0)

      do iel = 1, nelt
        do iface = 1,nsides
          cb=cbc(iface,iel)
          if(cb.eq.0) then
            do j=2,lx1-1
              do i=2,lx1-1
               tmmor(idmo(i,j,1,1,iface,iel))=0.d0
              end do
            end do

            j=1
            do i = 1, lx1-1
               tmmor(idmo(i,j,1,1,iface,iel))=0.d0
            end do

            if(idmo(lx1,1,1,1,iface,iel).eq.0)then
              tmmor(idmo(lx1,1,1,2,iface,iel))=0.d0
            else
              tmmor(idmo(lx1,1,1,1,iface,iel))=0.d0
              do i=1,lx1
                tmmor(idmo(i,j,1,2,iface,iel))=0.d0
              end do
            end if

            i=lx1
            if(idmo(lx1,2,1,2,iface,iel).eq.0)then
              do j=2,lx1-1
                tmmor(idmo(i,j,1,1,iface,iel))=0.d0
              end do
              tmmor(idmo(lx1,lx1,2,2,iface,iel))=0.d0
            else
              do j=2,lx1
                tmmor(idmo(i,j,1,2,iface,iel))=0.d0
              end do
              do j=1,lx1
                tmmor(idmo(i,j,2,2,iface,iel))=0.d0
              end do
            end if

            j=lx1
            tmmor(idmo(1,lx1,2,1,iface,iel))=0.d0
            if(idmo(2,lx1,2,1,iface,iel).eq.0)then
              do i=2,lx1-1
                tmmor(idmo(i,j,1,1,iface,iel))=0.d0
              end do
            else
              do i=2,lx1
                tmmor(idmo(i,j,2,1,iface,iel))=0.d0
              end do
              do i=1,lx1-1
                tmmor(idmo(i,j,2,2,iface,iel))=0.d0
              end do
            end if

            i=1
            do j=2,lx1-1
             tmmor(idmo(i,j,1,1,iface,iel))=0.d0
            end do
            if(idmo(1,lx1,1,1,iface,iel).ne.0)then
              tmmor(idmo(i,lx1,1,1,iface,iel))=0.d0
              do j=1,lx1-1
               tmmor(idmo(i,j,2,1,iface,iel))=0.d0
              end do
            end if

          endif
        end do
       end do

      return
      end

c------------------------------------------------------------------
      block data top_constants

c------------------------------------------------------------------
c.....We store some tables of useful topological constants
c------------------------------------------------------------------
      include 'header.h'

c     f_e_ef(e,f) returns the other face sharing the e'th local edge of face f.
      data f_e_ef/6,3,5,4, 6,3,5,4, 6,1,5,2, 6,1,5,2, 4,1,3,2, 4,1,3,2/

c.....e_c(n,j) returns n'th edge sharing the vertex j of an element
      data e_c /5,8,11, 1,4,11,  5,6,9, 1,2,9,
     &          7,8,12, 3,4,12, 6,7,10, 2,3,10/

c.....local_corner(n,i) returns the local corner index of vertex n on face i
      data local_corner /0,1,0,2,0,3,0,4, 1,0,2,0,3,0,4,0,
     &                   0,0,1,2,0,0,3,4, 1,2,0,0,3,4,0,0,
     &                   0,0,0,0,1,2,3,4, 1,2,3,4,0,0,0,0/

c.....cal_nnb(n,i) returns the neighbor elements neighbored by n'th edge
c     among the three edges sharing vertex i
c     the elements are the eight children elements ordered as 1 to 8.
      data cal_nnb/5,2,3, 6,1,4, 7,4,1, 8,3,2,
     &             1,6,7, 2,5,8, 3,8,5, 4,7,6/

c.....returns the opposite local corner index: 1-4,2-3
      data oplc /4,3,2,1/

c.....cal_iijj(i,n) returns the location of local corner number n on a face
c     i =1  to get ii, i=2 to get jj
c     (ii,jj) is defined the same as in mortar location (ii,jj)
      data cal_iijj /1,1, 1,2, 2,1, 2,2/

c.....returns the adjacent(neighbored by a face) element's children,
c     assumming a vertex is shared by eight child elements 1-8.
c     index n is local corner number on the face which is being
c     assigned the mortar index number
      data cal_intempx /8,6,4,2, 7,5,3,1, 8,7,4,3,
     $                  6,5,2,1, 8,7,6,5, 4,3,2,1/

c.....c_f(i,f) returns the vertex number of i'th local corner on face f
      data c_f /2,4,6,8, 1,3,5,7, 3,4,7,8, 1,2,5,6, 5,6,7,8, 1,2,3,4/

c.....on each face of the parent element, there are four children element.
c     le_arr(i,j,n) returns the i'th elements among the four children elements
c     n refers to the direction: 1 for x, 2 for y and 3 for z direction.
c     j refers to positive(0) or negative(1) direction on x, y or z direction.
c     n=1,j=0 refers to face 1 and n=1, j=1 refers to face 2, n=2,j=0 refers to
c     face 3....
c     The current eight children are ordered as 8,1,2,3,4,5,6,7
      data    le_arr/8,2,4,6, 1,3,5,7,
     $               8,1,4,5, 2,3,6,7,
     $               8,1,2,3, 4,5,6,7/

c.....jjface(n) returns the face opposite to face n
      data jjface /2,1,4,3,6,5/

cc.....edgeface(n,f) returns OTHER face which shares local edge n on face f
c      integer edgeface(4,6)
c      data edgeface /6,3,5,4, 6,3,5,4, 6,1,5,2,
c     $               6,1,5,2, 4,1,3,2, 4,1,3,2/

c.....e_face2(n,f) returns the local edge number of edge n on the
c     other face sharing local edge n on face f
      data e_face2 /2,2,2,2, 4,4,4,4, 3,2,3,2,
     $              1,4,1,4, 3,3,3,3, 1,1,1,1/

c.....op(n) returns the local edge number of the edge which
c     is opposite to local edge n on the same face
      data op /3,4,1,2/

c.....localedgenumber(f,e) returns the local edge number for edge e
c     on face f. A zero result value signifies illegal input
      data localedgenumber /1,0,0,0,0,2, 2,0,2,0,0,0, 3,0,0,0,2,0,
     $                      4,0,0,2,0,0, 0,1,0,0,0,4, 0,2,4,0,0,0,
     $                      0,3,0,0,4,0, 0,4,0,4,0,0, 0,0,1,0,0,3,
     $                      0,0,3,0,3,0, 0,0,0,1,0,1, 0,0,0,3,1,0/

c.....edgenumber(e,f) returns the edge index of local edge e on face f
      data edgenumber / 1,2, 3,4,  5,6, 7,8,  9,2,10,6,
     $                 11,4,12,8, 12,3,10,7, 11,1, 9,5/

c.....f_c(c,n) returns the face index of i'th face sharing vertex n
      data f_c /2,4,6, 1,4,6, 2,3,6, 1,3,6,
     &          2,4,5, 1,4,5, 2,3,5, 1,3,5/

c.....if two elements are neighbor by one edge,
c     e1v1(f1,f2) returns the smaller index of the two vertices on this
c     edge on one element
c     e1v2 returns the larger index of the two vertices of this edge on
c     on element. exfor a vertex on element
c     e2v1 returns the smaller index of the two vertices on this edge on
c     another element
c     e2v2 returns the larger index of the two vertiex on this edge on
c     another element
      data e1v1/0,0,4,2,6,2, 0,0,3,1,5,1, 4,3,0,0,7,3,
     &          2,1,0,0,5,1, 6,5,7,5,0,0, 2,1,3,1,0,0/
      data e2v1/0,0,1,3,1,5, 0,0,2,4,2,6, 1,2,0,0,1,5,
     &          3,4,0,0,3,7, 1,2,1,3,0,0, 5,6,5,7,0,0/
      data e1v2/0,0,8,6,8,4, 0,0,7,5,7,3, 8,7,0,0,8,4,
     &          6,5,0,0,6,2, 8,7,8,6,0,0, 4,3,4,2,0,0/
      data e2v2/0,0,5,7,3,7, 0,0,6,8,4,8, 5,6,0,0,2,6,
     &          7,8,0,0,4,8, 3,4,2,4,0,0, 7,8,6,8,0,0/

c.....children(n1,n)returns the four elements among the eight children
c     elements to be merged on face n of the parent element
c     the IDs for the eight children are 1,2,3,4,5,6,7,8
      data children/2,4,6,8, 1,3,5,7, 3,4,7,8,
     &              1,2,5,6, 5,6,7,8, 1,2,3,4/

c.....iijj(n1,n) returns the location of n's mortar on an element face
c     n1=1 refers to x direction location and n1=2 refers to y direction
      data iijj/1,1,1,2,2,1,2,2/

c.....v_end(n) returns the index of collocation points at two ends of each
c     direction
      data v_end /1,lx1/

c.....face_l1,face_l2,face_ld return for start,end,stride for a loop over faces
c     used on subroutine  mortar_vertex
      data face_l1 /2,3,1/, face_l2 /3,1,2/, face_ld /1,-2,1/

      end

max_stars_repo_name: jonarbo/KUBE

max_issues_repo_name: jonarbo/KUBE

max_forks_repo_name: jonarbo/KUBE


text:
      SUBROUTINE BMG3_SymStd_SETUP_cg_LU(&
     &                SO, II, JJ, KK, NStncl, ABD, NABD1, NABD2,&
     &                WS, NMSGr, NProcI, NProcJ, NProcK, NProc, MyProc,&
     &                ProcGrid, ProcCoord, LocArrSize, MPICOMM&
     &                ) BIND(C, NAME='MPI_BMG3_SymStd_SETUP_cg_LU')

! ======================================================================
!  --------------------
!   DESCRIPTION:
!  --------------------
!
!     BMG3_SymStd_SETUP_cg_LU sets up the matrix on the coarsest grid,
!     and using the LAPACK routine DPBTRF, it forms the LU decomposition
!     of the matrix.
!
! ======================================================================
! $license_flag$
! ======================================================================
!  --------------------
!   INPUT:
!  --------------------
!
!     SO        Refer to BMG2_SymStd_SOLVE_boxmg
!     II        Number of grid points in x direction, including
!               two fictitious points.
!     JJ        Number of grid points in y direction, including
!               two fictitious points.
!     KK        Number of grid points in z direction, including
!               two fictitious points.
!
! ======================================================================
!  --------------------
!   INPUT/OUTPUT:
!  --------------------
!
!
! ======================================================================
!  --------------------
!   OUTPUT:
!  --------------------
!
!     ABD       Refer to BMG3_SymStd_SOLVE_boxmg
!
! ======================================================================
!  --------------------
!   LOCAL:
!  --------------------
!
!
! ======================================================================

      USE ModInterface
      IMPLICIT NONE

! ----------------------------
!     Includes
!
      INCLUDE 'BMG_stencils_f90.h'
      INCLUDE 'BMG_constants_f90.h'
      INCLUDE 'BMG_parameters_f90.h'

      INCLUDE 'mpif.h'

! ----------------------------
!     Argument Declarations
!
      integer(len_t), value :: ii, jj, kk, nabd1, nabd2, NMSGr
      integer(c_int), value :: NStncl, NProcI, NProcj, NProcK, NProc,&
           MyProc, MPICOMM
      integer(c_int) :: ProcGrid(NProcI,NProcJ,NProcK), ProcCoord(3, NProc)
      integer(len_t) :: LocArrSize(3,*)
      real(real_t) :: abd(nabd1,nabd2), so(ii+1,jj+1,kk+1,NStncl)
      real(real_t) :: WS(NMSGr)

! ----------------------------
!     Local Declarations
!
      INTEGER i, i1, i2, ibw, j, j1, k, k1, kl, KK_t, n
      INTEGER IIG, JJG, KKG, proc, KKMAX, P1, P2, P3
      INTEGER P1SUM, P2SUM, P3SUM, PXP1, PYP2, PZP3, p_WS
      INTEGER INT_TEMP, INT_TEMP1, INT_TEMP2, INT_TEMP3
      INTEGER LARGESTNODES, IERR, MPI_IERROR
      INTEGER cg_comm_type, info

! ======================================================================

      ! hardcoded for now
      cg_comm_type = BMG_CG_ALLGATHER
!
!     set up and factor matrix for direct solve on coarsest grid.
!
      info = 0

      i1=ii-1
      j1=jj-1
      k1=kk-1

      i2=i1-1
!      n=i2*(j1-1)*(kk-2)

! ------------------------------------------------
!     Calculate the global number of points
! ------------------------------------------------

      !
      ! Global number in x
      !
      IIG=2
      DO i=1,NProcI
         IIG = IIG + LocArrSize(1,ProcGrid(i,1,1)) - 2
      ENDDO

      !
      ! Global number in y
      !
      JJG=2
      DO j=1,NProcJ
         JJG = JJG + LocArrSize(2,ProcGrid(1,j,1)) - 2
      ENDDO

      !
      ! Global number in z
      !
      KKG=2
      DO k=1,NProcK
         KKG = KKG + LocArrSize(3,ProcGrid(1,1,k)) - 2
      ENDDO

! ------------------------------------------------
!     Find the largest local array
! ------------------------------------------------

      !
      ! Loop over local x dimensions
      !
      INT_TEMP = LocArrSize(1,1) - 2
      INT_TEMP1 = INT_TEMP
      DO proc=2, NProcI*NProcJ*NProcK
         INT_TEMP = LocArrSize(1,proc) - 2
         IF (INT_TEMP.gt.INT_TEMP1) THEN
            INT_TEMP1 = INT_TEMP
         END IF
      END DO

      !
      ! Loop over local y dimensions
      !
      INT_TEMP = LocArrSize(2,1) - 2
      INT_TEMP2 = INT_TEMP
      DO proc=2, NProcI*NProcJ*NProcK
         INT_TEMP = LocArrSize(2,proc) - 2
         IF (INT_TEMP.gt.INT_TEMP2) THEN
            INT_TEMP2 = INT_TEMP
         END IF
      END DO

      !
      ! Loop over local z dimensions
      !
      INT_TEMP = LocArrSize(3,1) - 2
      INT_TEMP3 = INT_TEMP
      DO proc=2, NProcI*NProcJ*NProcK
         INT_TEMP = LocArrSize(3,proc) - 2
         IF (INT_TEMP.gt.INT_TEMP3) THEN
            INT_TEMP3 = INT_TEMP
         END IF
      END DO

      !
      ! Conservative: take largest from each
      !
      LARGESTNODES = INT_TEMP1 * INT_TEMP2 * INT_TEMP3 * 14 + 4

! ------------------------------------------------
!     Copy all information into the buffer
! ------------------------------------------------

      WS(1) = MyProc
      WS(2) = II
      WS(3) = JJ
      WS(4) = KK
      INT_TEMP = 5

      IF ( NStncl.EQ.14 ) THEN
         !
         DO K = 2, KK - 1
            DO J = 2, JJ - 1
               DO I = 2, II - 1
                  WS(INT_TEMP)    =  SO(I  ,J  ,k  ,kp)
                  WS(INT_TEMP+ 1) = -SO(I  ,J  ,k  ,kpw)
                  WS(INT_TEMP+ 2) = -SO(I+1,J  ,k  ,kpnw)
                  WS(INT_TEMP+ 3) = -SO(I  ,J  ,k  ,kps)
                  WS(INT_TEMP+ 4) = -SO(I  ,J  ,k  ,kpsw)
                  WS(INT_TEMP+ 5) = -SO(I+1,J+1,k  ,kbne)
                  WS(INT_TEMP+ 6) = -SO(I  ,J+1,k  ,kbn)
                  WS(INT_TEMP+ 7) = -SO(I  ,J+1,k  ,kbnw)
                  WS(INT_TEMP+ 8) = -SO(I+1,J  ,k  ,kbe)
                  WS(INT_TEMP+ 9) = -SO(I  ,J  ,k  ,kb)
                  WS(INT_TEMP+10) = -SO(I  ,J  ,k  ,kbw)
                  WS(INT_TEMP+11) = -SO(I+1,J  ,k  ,kbse)
                  WS(INT_TEMP+12) = -SO(I  ,J  ,k  ,kbs)
                  WS(INT_TEMP+13) = -SO(I  ,J  ,k  ,kbsw)
                  INT_TEMP = INT_TEMP + 14
               END DO
            END DO
         END DO
         !
      ELSE IF ( NStncl.EQ.4 ) THEN
         !
         DO K = 2, KK - 1
            DO J = 2, JJ - 1
               DO I = 2, II - 1
                  WS(INT_TEMP)    =  SO(I  ,J  ,k  ,kp)
                  WS(INT_TEMP+ 1) = -SO(I  ,J  ,k  ,kpw)
                  WS(INT_TEMP+ 2) = rZERO
                  WS(INT_TEMP+ 3) = -SO(I  ,J  ,k  ,kps)
                  WS(INT_TEMP+ 4) = rZERO
                  WS(INT_TEMP+ 5) = rZERO
                  WS(INT_TEMP+ 6) = rZERO
                  WS(INT_TEMP+ 7) = rZERO
                  WS(INT_TEMP+ 8) = rZERO
                  WS(INT_TEMP+ 9) = -SO(I  ,J  ,k  ,kb)
                  WS(INT_TEMP+10) = rZERO
                  WS(INT_TEMP+11) = rZERO
                  WS(INT_TEMP+12) = rZERO
                  WS(INT_TEMP+13) = rZERO
                  INT_TEMP = INT_TEMP + 14
               END DO
            END DO
         END DO
         !
      ELSE
         call print_error(C_CHAR_"Cholesky failed! (incorrect NStncl)"//C_NULL_CHAR)
         RETURN
      ENDIF

! ------------------------------------------------
!     Send/Receive information to/from everybody
! ------------------------------------------------

      IF (cg_comm_type .eq. BMG_CG_ALLGATHER) THEN

         CALL MPI_ALLGATHER ( &
     &        WS(1), LARGESTNODES, MPI_DOUBLE_PRECISION,&
     &        WS(LARGESTNODES+1), LARGESTNODES, MPI_DOUBLE_PRECISION,&
     &        MPICOMM, IERR&
     &        )

      ELSE

         CALL MPI_GATHER ( &
     &        WS(1), LARGESTNODES, MPI_DOUBLE_PRECISION,&
     &        WS(LARGESTNODES+1), LARGESTNODES, MPI_DOUBLE_PRECISION,&
     &        iZERO, MPICOMM, IERR&
     &        )

      ENDIF

!     ------------------------------------------------
!     Assemble the global coarse grid matrix
!     ------------------------------------------------

      IF ( cg_comm_type .eq. BMG_CG_ALLGATHER .OR.&
     &     MyProc .eq. iONE ) THEN

         KKMAX=0

         DO proc=1, NProcI*NProcJ*NProcK

            INT_TEMP = proc * LARGESTNODES + 1
            INT_TEMP = WS(INT_TEMP)

            P1 = ProcCoord(1,INT_TEMP)
            P2 = ProcCoord(2,INT_TEMP)
            P3 = ProcCoord(3,INT_TEMP)

            P1SUM = 0
            DO i=1, P1-1
               P1SUM = P1SUM + LocArrSize(1,ProcGrid(i,1,1)) - 2
            END DO

            P2SUM = 0
            DO j=1, P2-1
               P2SUM = P2SUM + LocArrSize(2,ProcGrid(1,j,1)) - 2
            END DO

            P3SUM = 0
            DO k=1, P3-1
               P3SUM = P3SUM + LocArrSize(3,ProcGrid(1,1,k)) - 2
            END DO

            PXP1 = LocArrSize(1,INT_TEMP) - 2
            PYP2 = LocArrSize(2,INT_TEMP) - 2
            PZP3 = LocArrSize(3,INT_TEMP) - 2

            IBW = (IIG-2)*(JJG-1)+1

            DO k=1, PZP3
               DO j=1, PYP2
                  DO i=1, PXP1

                     KK_t = (P1SUM + i - 1) + (IIG-2) * (P2SUM + j - 1) &
     &                    + (IIG-2)*(JJG-2)*(P3SUM + k - 1) + 1

                     p_WS = proc*LARGESTNODES + 5 &
     &                    + 14*((i-1) + (j-1)*PXP1 + (k-1)*PXP1*PYP2)

                     ABD(IBW+1,KK_t)                 = WS( p_WS ) ! kp
                     ABD(IBW  ,KK_t)                 = WS( p_WS + 1 ) !

                     ABD(IBW-IIG+4,KK_t)             = WS( p_WS + 2 ) !
                     ABD(IBW-IIG+3,KK_t)             = WS( p_WS + 3 ) !
                     ABD(IBW-IIG+2,KK_t)             = WS( p_WS + 4 ) !

                     ABD(IBW-(JJG-3)*(IIG-2)+2,KK_t) = WS( p_WS + 5 ) !
                     ABD(IBW-(JJG-3)*(IIG-2)+1,KK_t) = WS( p_WS + 6 ) !
                     ABD(IBW-(JJG-3)*(IIG-2)  ,KK_t) = WS( p_WS + 7 ) !

                     ABD(IBW-(JJG-2)*(IIG-2)+2,kk_t) = WS( p_WS + 8 ) !
                     ABD(IBW-(JJG-2)*(IIG-2)+1,kk_t) = WS( p_WS + 9 ) !
                     ABD(IBW-(JJG-2)*(IIG-2)  ,kk_t) = WS( p_WS + 10 ) !

                     ABD(3,KK_t)                     = WS( p_WS + 11 ) !
                     ABD(2,KK_t)                     = WS( p_WS + 12 ) !
                     ABD(1,KK_t)                     = WS( p_WS + 13 ) !

                     IF (KK_t.gt.KKMAX) KKMAX = KK_t

                  END DO
               END DO
            END DO

         END DO

         CALL DPBTRF('U', KKMAX, IBW, ABD, NABD1, INFO)

         IF (INFO .NE. 0) THEN

            call print_error("Coarse grid Cholesky decomp failed!"//C_NULL_CHAR)

            if (MyProc .eq. 1) then
               WRITE(*,*) 'INFO = ', INFO
            endif

         END IF

      END IF

! ======================================================================

 150  FORMAT(I4,4X,11(E10.4,4X))

! ===========================================

      RETURN
      END

max_stars_repo_name: cedar-framework/cedar

max_issues_repo_name: cedar-framework/cedar

max_forks_repo_name: cedar-framework/cedar


text:
!>------------------------------------------------
!! Defines model constants (e.g. gravity, and MAXFILELENGTH)
!!
!!------------------------------------------------
module icar_constants

    implicit none

    character(len=5) :: kVERSION_STRING = "2.0a1"

    ! string lengths
    integer, parameter :: kMAX_FILE_LENGTH = 1024
    integer, parameter :: kMAX_DIM_LENGTH  = 1024
    integer, parameter :: kMAX_NAME_LENGTH = 1024
    integer, parameter :: kMAX_ATTR_LENGTH = 1024

    !>--------------------------------------------
    ! list of integer constants to be used when accessing various arrays that track variable allocation, usage, etc. requests
    !
    ! NOTE: IF YOU ADD TO THIS LIST BE SURE TO ADD AN INTEGER TO THE kVARS STRUCTURE CONSTRUCTOR BELOW IT!
    ! This could be transitioned to an enum... but then one can't just "use, only:kVARS"...
    ! enum, bind(C)
    !   enumerator ::  u, v, w,...
    ! end enum
    ! --------------------------------------------
    type var_constants_type
        SEQUENCE    ! technically SEQUENCE just requires the compiler leave them in order,
                    ! but it can also keep compilers (e.g. ifort) from padding for alignment,
                    ! as long as there is no padding we can test last_var = sizeof(kVARS)

        integer :: u
        integer :: v
        integer :: w
        integer :: w_real
        integer :: pressure
        integer :: pressure_interface
        integer :: potential_temperature
        integer :: temperature
        integer :: water_vapor
        integer :: cloud_water
        integer :: cloud_number_concentration
        integer :: cloud_ice
        integer :: ice_number_concentration
        integer :: rain_in_air
        integer :: rain_number_concentration
        integer :: snow_in_air
        integer :: snow_number_concentration
        integer :: graupel_in_air
        integer :: graupel_number_concentration
        integer :: precipitation
        integer :: convective_precipitation
        integer :: snowfall
        integer :: graupel
        integer :: exner
        integer :: density
        integer :: z
        integer :: z_interface
        integer :: dz
        integer :: dz_interface
        integer :: cloud_fraction
        integer :: shortwave
        integer :: longwave
        integer :: vegetation_fraction
        integer :: veg_type
        integer :: soil_type
        integer :: lai
        integer :: canopy_water
        integer :: sensible_heat
        integer :: latent_heat
        integer :: u_10m
        integer :: v_10m
        integer :: ustar
        integer :: temperature_2m
        integer :: humidity_2m
        integer :: surface_pressure
        integer :: longwave_up
        integer :: ground_heat_flux
        integer :: soil_totalmoisture
        integer :: soil_deep_temperature
        integer :: roughness_z0
        integer :: snow_water_equivalent
        integer :: soil_water_content
        integer :: soil_temperature
        integer :: skin_temperature
        integer :: sst
        integer :: land_mask
        integer :: terrain
        integer :: latitude
        integer :: longitude
        integer :: u_latitude
        integer :: u_longitude
        integer :: v_latitude
        integer :: v_longitude
        integer :: tend_qv_adv
        integer :: tend_qv_pbl
        integer :: tend_qv
        integer :: tend_th
        integer :: tend_qc
        integer :: tend_qi
        integer :: tend_qs
        integer :: tend_qr
        integer :: tend_u
        integer :: tend_v
        integer :: znu
        integer :: znw
        integer :: last_var
    end type var_constants_type

    type(var_constants_type) :: kVARS = var_constants_type( 1,  2,  3,  4,  5,  6,  7,  8,  9, 10,  &
                                                           11, 12, 13, 14, 15, 16, 17, 18, 19, 20,  &
                                                           21, 22, 23, 24, 25, 26, 27, 28, 29, 30,  &
                                                           31, 32, 33, 34, 35, 36, 37, 38, 39, 40,  &
                                                           41, 42, 43, 44, 45, 46, 47, 48, 49, 50,  &
                                                           51, 52, 53, 54, 55, 56, 57, 58, 59, 60,  &
                                                           61, 62, 63, 64, 65, 66, 67, 68, 69, 70,  &
                                                           71, 72, 73, 74, 75, 76 )

    integer, parameter :: kINTEGER_BITS     = storage_size(kINTEGER_BITS)
    integer, parameter :: kMAX_STORAGE_VARS = storage_size(kVARS) / kINTEGER_BITS

    ! Initial number of output variables for which pointers are created
    integer, parameter :: kINITIAL_VAR_SIZE= 128

    ! Maximum number of dimensions
    ! Note this is defined in NetCDF, though not enforced (a file can have more than 1024 dimensions)
    integer, parameter :: kMAX_DIMENSIONS  = 1024

!>------------------------------------------------
!! Model constants (mostly string lengths)
!! ------------------------------------------------
    integer, parameter :: MAXFILELENGTH      =   1024  ! maximum file name length
    integer, parameter :: MAXVARLENGTH       =   1024  ! maximum variable name length
    integer, parameter :: MAXLEVELS          =    500  ! maximum number of vertical layers (should typically be ~10-20)
    integer, parameter :: MAX_NUMBER_FILES   =  50000  ! maximum number of permitted input files (probably a bit extreme)
    integer, parameter :: MAXSTRINGLENGTH    =   1024  ! maximum length of other strings (e.g. netcdf attributes)
    integer, parameter :: kMAX_STRING_LENGTH =   1024  ! maximum length of other strings (e.g. netcdf attributes)

!>------------------------------------------------
!!  Default width of coarray halos, ideally might be physics dependant (e.g. based on advection spatial order)
!! ------------------------------------------------
    integer,parameter :: kDEFAULT_HALO_SIZE = 1

!>------------------------------------------------
!! Value to accept for difference between real numbers should be as a fraction but then have to test for non-zero...
!! For some variables (cloud ice) 1e-6 is not that small, for others (pressure) it might be below precision...
!! ------------------------------------------------
    real,   parameter :: kSMALL_VALUE = 1e-6

    integer, parameter :: kMAINTAIN_LON      = 0
    integer, parameter :: kPRIME_CENTERED    = 1
    integer, parameter :: kDATELINE_CENTERED = 2
    integer, parameter :: kGUESS_LON         = 3

! ------------------------------------------------
! Physics scheme selection definitions
!
! NB: BASIC typically means "use the data from the low res model"
!     SIMPLE typically means a relatively simple formulation written for ICAR
! These could all be switched to enums too, but this makes it easy to see what number each has for the options file...
! ------------------------------------------------
    integer, parameter :: kNO_STOCHASTIC = -9999
    integer, parameter :: kCU_TIEDTKE    = 1
    integer, parameter :: kCU_SIMPLE     = 2
    integer, parameter :: kCU_KAINFR     = 3

    integer, parameter :: kMP_THOMPSON   = 1
    integer, parameter :: kMP_SB04       = 2
    integer, parameter :: kMP_MORRISON   = 3
    integer, parameter :: kMP_WSM6       = 4
    integer, parameter :: kMP_THOMP_AER  = 5

    integer, parameter :: kPBL_BASIC     = 1
    integer, parameter :: kPBL_SIMPLE    = 2
    integer, parameter :: kPBL_YSU       = 3

    integer, parameter :: kWATER_BASIC   = 1
    integer, parameter :: kWATER_SIMPLE  = 2

    integer, parameter :: kLSM_BASIC     = 1
    integer, parameter :: kLSM_SIMPLE    = 2
    integer, parameter :: kLSM_NOAH      = 3

    integer, parameter :: kRA_BASIC      = 1
    integer, parameter :: kRA_SIMPLE     = 2

    integer, parameter :: kADV_UPWIND    = 1
    integer, parameter :: kADV_MPDATA    = 2

    integer, parameter :: kWIND_LINEAR   = 1
    integer, parameter :: kCONSERVE_MASS = 2
    integer, parameter :: kITERATIVE_WINDS = 3

    integer, parameter :: kLC_LAND       = 1
    integer, parameter :: kLC_WATER      = 2

    ! mm of accumulated precip before "tipping" into the bucket
    ! only performed on output operations
    integer, parameter :: kPRECIP_BUCKET_SIZE=100

! ------------------------------------------------
! Physical Constants
! ------------------------------------------------
    real, parameter :: LH_vaporization=2260000.0 ! J/kg
    ! could be calculated as 2.5E6 + (-2112.0)*temp_degC ?
    real, parameter :: Rd  = 287.058   ! J/(kg K) specific gas constant for dry air
    real, parameter :: Rw  = 461.5     ! J/(kg K) specific gas constant for moist air
    real, parameter :: cp  = 1012.0    ! J/kg/K   specific heat capacity of moist STP air?
    real, parameter :: gravity= 9.81   ! m/s^2    gravity
    real, parameter :: pi  = 3.1415927 ! pi
    real, parameter :: stefan_boltzmann = 5.67e-8 ! the Stefan-Boltzmann constant
    real, parameter :: karman = 0.41   ! the von Karman constant

    ! convenience parameters for various physics packages
    real, parameter :: rovcp = Rd/cp
    real, parameter :: rovg  = Rd/gravity

    ! from wrf module_model_constants
    ! parameters for calculating latent heat as a function of temperature for
    ! vaporization
    real, parameter ::  XLV0 = 3.15E6
    real, parameter ::  XLV1 = 2370.
    ! sublimation
    real, parameter ::  XLS0 = 2.905E6
    real, parameter ::  XLS1 = 259.532

    ! saturated vapor pressure parameters (?)
    real, parameter ::  SVP1 = 0.6112
    real, parameter ::  SVP2 = 17.67
    real, parameter ::  SVP3 = 29.65
    real, parameter ::  SVPT0= 273.15

    real, parameter ::  EP1  = Rw/Rd-1.
    real, parameter ::  EP2  = Rd/Rw

end module

max_stars_repo_name: scrasmussen/icar

max_issues_repo_name: scrasmussen/icar

max_forks_repo_name: scrasmussen/icar


text:
Payphone of Venus on the Downtown Phones node

Phone Number: (530) 7579736

Location: This phone is located on the Delta patio. Delta of Venus is on 2nd Street 2nd and B St.

Description:Youre likely to have somone pick up the phone if you call during business hours. Even in cold weather, theres likely somebody outside smoking who may pick up the phone. This phone has a reputation for ringing on a regular basis. If you pick it up, there is often no one on the other end. Unless perhaps its YOU that is calling it after finding the number through the Payphone Project.

Charge: 50 local (15 mins); 75 instate (3 mins); $1 US or Mexico (3 min)

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
	SUBROUTINE GGDRIV ( grid, grid1, kx, ky, hist, work1, work2,
     +			    work3, data, infoflg, iret )
C************************************************************************
C* GGDRIV                                                               *
C*                                                                      *
C* This subroutine takes all the original GRPHGD input parameters and	*
C* processes it into a GEMPAK grid file.				*
C*                                                                      *
C* GGDRIV ( GRID, GRID1, KX, KY, HIST, WORK1, WORK2, WORK3, DATA	*
C*          INFOFLG, IRET )						*
C*                                                                      *
C* Input parameters:                                                    *
C*                                                                      *
C* Output parameters:                                                   *
C*      GRID(*)         REAL            Data grid			*
C*      GRID1(*)        REAL            Data grid1			*
C*      KX              INTEGER         Data grid x-dimension		*
C*      KY              INTEGER         Data grid y-dimension		*
C*      HIST(*)         REAL            History data grid		*
C*      WORK1(*)        REAL            Latitude in degrees             *
C*      WORK2(*)        REAL            Longitude in degrees            *
C*      WORK3(*)        REAL            Cosine of latitude              *
C*      DATA(*)         REAL            Surface data			*
C*      INFOFLG         LOGICAL         Info flag			*
C*      IRET            INTEGER         Return code                     *
C*					= >0  - misc warnings		*
C*					= -1  - misc error 		*
C*					= -8  - unable to create gd file*
C*					= -13 - grid write problem	*
C*                                                                     	*
C**                                                                     *
C* Log:                                                                 *
C* D.W.Plummer/NCEP	10/99   From the original GRPHGD program	*
C* D.W.Plummer/NCEP	10/99   Remove user option to exit calls	*
C* S. Jacobs/NCEP	11/99	Changed GSMPRJ to GSGPRJ		*
C* D.W.Plummer/NCEP	 2/00	Added grid,kx,ky to calling sequence	*
C* D.W.Plummer/NCEP	 3/00	Added subarea option to CPYFIL		*
C* D.W.Plummer/NCEP	 5/00	Calling sequence change; use IP library	*
C* T. Lee/GSC		 6/00	Checked navigation and error message	*
C* T. Lee/GSC		10/00	Output default grid value		*
C* S. Jacobs/NCEP	 3/01	Increased file template to 48 chars	*
C* D.W.Plummer/NCEP	12/02	Changes for bounds, minmax, gen upgrade	*
C* D.W.Plummer/NCEP	 7/03	Add calls to GR_FIXA and GR_MAPS	*
C* D.W.Plummer/NCEP	 7/03	Chgs for intersection efficiency	*
C* D.W.Plummer/NCEP	 9/03	Chgs for relative mm, bounds		*
C* T. Lee/SAIC		 7/04	Added type, oaattr and oaguess to GRGINP*
C* H. Zeng/SAIC         03/05   Added catmap to GRGINP                  *
C* D.W.Plummer/NCEP     03/05   Increased grid limit to LLMXTG          *
C* R. Tian/SAIC          3/05   Modfied to adapt new file/time mgnt     *
C* H. Zeng/SAIC		04/05	Added discrete&dlines to grpinp		*
C* M. Li/SAIC		05/05	Process vector information		*
C* D.W.Plummer/NCEP	 5/05	Add DG_INTL; add ER_WMSG after DG_NWDT	*
C* D.W.Plummer/NCEP	 6/05	Bug fix to assign RMISSD properly	*
C* M. Li/SAIC		11/05	Updated error msg for ggrdinfo		*
C* m.gamazaychikov/SAIC	12/05	Add call to GG_RDVGF			*
C* T. Lee/SAIC          12/05   Initialized ighdr                       *
C* M. Li/SAIC		04/06	Added input to IN_DLINQ			*
C* S. Gilbert/NCEP	10/06	Added another ER_WMSG after GDGCFL	*
C* M. Li/SAIC		03/07	Added new input parameter EDGEOPTS	*
C* M. Li/SAIC		04/07   Added new function GGBNDCONDS		*
C* F. J. Yen/NCEP	01/08	Added parameters OLKDAY and KEYLINE	*
C************************************************************************
	INCLUDE 	'GEMPRM.PRM'
	INCLUDE		'grphgd.cmn'

	INTEGER 	iret, ier, kx, ky
C
	REAL 		grid(*), grid1(*), hist(*), work1(*), work2(*),
     +  		work3(*), data(*)
        LOGICAL		infoflg
C*
        CHARACTER       gdfile*(LLMXLN), proj*(LLMXLN), cpyfil*(LLMXLN),
     +                  gdarea*(LLMXLN), anlyss*(LLMXLN), kxky*(LLMXLN),
     +                  maxgrd*(LLMXLN), gparm*(LLMXLN),
     +                  cntrfl*(LLMXLN), gdatim*(LLMXLN),
     +                  gvcord*(LLMXLN), glevel*(LLMXLN),
     +                  ckeycl*(LLMXLN), gglimt*(LLMXLN),
     +                  hstgrd*(LLMXLN), bounds*(LLMXLN),
     +			type*(LLMXLN),   catmap*(LLMXLN),
     +			discrete*(LLMXLN), dlines*(LLMXLN),
     +                  guess*(LLMXLN),  guesfun*(LLMXLN),
     +                  cgamma*(LLMXLN), csrch*(LLMXLN),
     +                  pass*(LLMXLN),   cqcntl*(LLMXLN),
     +                  ggvgf*(LLMXLN),  edgeopts*(LLMXLN),
     +                  olkday*(LLMXLN), keylin*(LLMXLN)
C*
        CHARACTER       newfil*256, filnam*(LLMXLN)
        LOGICAL         proces, exist, gottm
C
	INTEGER		ighdr(LLGDHD), level(2), ivcord
	CHARACTER	parm*12, gdattm(2)*20, vparm*(LLMXLN),
     +			vgsub(2)*20
	LOGICAL		rewrit/.true./, hstflg
C
	REAL		grdmin, grdmax
	INTEGER		nmiss, ngpts
	CHARACTER	garout*72, prjout*72, imgfil*72/' '/,
     +	                gpack*12/' '/
C*
C-----------------------------------------------------------------------
C
	iret = 0
	proces = .true.
C
C*	Initialize DG library
C
	CALL DG_INTL ( ier )
C
C*	Get user input
C
        CALL GRGINP   ( gdfile, guess, proj, gdarea, kxky, maxgrd,
     +			cpyfil, anlyss, cntrfl, gdatim, gparm,
     +			glevel, gvcord, ckeycl, gglimt,
     +			hstgrd, bounds, type, cgamma, csrch,
     +			pass, cqcntl, guesfun, catmap, discrete,
     +			dlines, ggvgf, edgeopts, olkday, keylin, iperr )

        CALL ST_NULL ( ggvgf, ggvgf, lens, ier )
	IF  ( iperr .ne. 0 )  THEN
	    CALL ER_WMSG ( 'GRPHGD', iperr, ' ', ier )
	    iret = -1
	    RETURN
	END IF
C
        CALL ST_NULL ( catmap, catmap, lens, ier )
	CALL IN_CATMINP ( catmap, ier )
C
        CALL ST_NULL ( discrete, discrete, lens, ier )
	CALL IN_DISCRETE ( discrete, ier )
C
        CALL ST_NULL ( dlines, dlines, lens, ier )
	CALL IN_DLINES ( dlines, ier )
C
C*	Process LIMITS input.
C
	CALL GLIMIT ( gglimt, ier )
C
C*      Create the grid file
C
	CALL FL_INQR ( gdfile, exist, newfil, ier )

	IF ( .not. exist ) THEN
	    CALL GDGCFL ( gdfile, proj, gdarea, kxky, maxgrd, cpyfil,
     +                    anlyss, ier )
	    IF ( ier .ne. 0 ) THEN
	        CALL ER_WMSG  ( 'GDCFIL', ier, ' ', ier1 )
	        iret = -8
	        CALL ER_WMSG  ( 'GRPHGD', iret, gdfile, ier )
	        RETURN
	    END IF
	END IF
C
C*      Process the GDFILE input
C
	IF ( proces ) THEN
            CALL DG_NFIL ( gdfile, gdfile, ier )
            IF ( ier .ne. 0 ) THEN
                CALL ER_WMSG ( 'DG', ier, ' ', irr )
                proces = .false.
            END IF
	END IF
C
C*      Process the GDATTIM input
C
	IF ( proces ) THEN
            CALL DG_NDTM ( gdatim, ier )
            IF ( ier .ne. 0 ) THEN
                CALL ER_WMSG ( 'DG', ier, gdatim, irr )
                proces = .false.
            END IF
	END IF
C
C*      Get the next time to process.
C
	IF ( proces ) THEN
            CALL DG_NTIM ( .true., .false., gdattm, gottm, ier )
            proces = ( ier .eq. 0 .and. gottm )
	END IF
C
C*      Set the map projection and graphics area.
C
	IF ( proces )  THEN
            CALL DG_FIXA  ( 'dset', ' ', garout, prjout, ier )
 	    CALL GG_MAPS ( prjout, garout, imgfil, idrpfl, ier )
	END IF
C
	CALL DG_QKXY ( kx, ky, ier )
C
C*      Check to see if grid is too large.
C
        IF ( ( kx * ky ) .gt. LLMXTG )  THEN
            iret = +1
            CALL ER_WMSG  ( 'GRPHGD', iret, ' ', ier )
            WRITE(6,*) 'Number of grid points (', kx*ky,
     &                  ') exceeds maximum allowed (', LLMXTG, ')'
            proces = .false.
        END IF
C
C*      Check to see if grid dimensions are too big for
C*	intersection arrays.
C

        IF ( ( kx + ky ) .gt. MAXDIM )  THEN
	    iret = +2
	    CALL ER_WMSG  ( 'GRPHGD', iret, ' ', ier )
         WRITE(6,*) 'Total of grid dimensions (', kx+ky,
     & 			') exceeds maximum allowed (', MAXDIM, ')'
	    proces = .false.
        END IF
C
	DO ii = 1, LLGDHD
	    ighdr ( ii ) = 0
	END DO
C
	IF ( proces )  THEN
	    CALL ST_LCUC  ( hstgrd, hstgrd, ier )
	    hstflg = hstgrd (1:1) .eq. 'Y'
C
C*	    Initialize history and data arrays
C
	    DO  ij = 1, kx*ky
		hist(ij) = INIT
		grid(ij) = RMISSD
		grid1(ij) = RMISSD
		work1(ij) = RMISSD
		work2(ij) = RMISSD
		work3(ij) = RMISSD
	    END DO
C
C*	    Get level value, vertical coord
C
	    CALL ST_ILST ( glevel, ':', -1, 2, level, num, ier )
	    CALL ST_LCUC ( gvcord, gvcord, ier )
	    CALL LV_CORD ( gvcord, vparm, ivcord, ier )
C
            imn = 1
	    imx = kx
	    jmn = 1
	    jmx = ky
C
C*	    Get keycol
C
	    CALL ST_ILST ( ckeycl, ' ', 0, 1, keycol, knum, ier )
C
C*	    Get iolkdy from GRPHGD parameter OLKDAY
C
	    CALL ST_ILST ( olkday, ' ', 0, 1, iolkdy, knum, ier )
C
C*	    Get and store vgtype (in kvgtyp) and subtyp (in ksubtp) of
C*	    line from GRPHGD parameter KEYLINE, which defines the
C*	    line types to process.  KEYLINE can be 0 or vgtype/subtyp.
C*	    Allowable values for vgtype are 'LINE' and 'SPLN'.
C*	    Invalid or default value for KEYLINE will be 0, which
C*	    means process all line types.
C
	    CALL ST_CLST ( keylin, '/', '0', 2, vgsub, knum, ier )
	    CALL ST_NULL (vgsub(1), vgsub(1), lens, ier)
	    CALL ST_LCUC ( vgsub(1), vgsub(1), ier )
	    IF ( vgsub(1)(1:lens) .eq. 'LINE' ) THEN
	        CALL ST_ILST ( vgsub(2), ' ', 0, 1, ksubtp, knum, ier )
		IF ( ksubtp .lt. 1 .or. ksubtp .gt. 9 ) THEN
		    ksubtp = 0
		    kvgtyp = 0
		ELSE
		    kvgtyp = 1
		END IF
	      ELSE IF ( vgsub(1)(1:lens) .eq. 'SPLN' ) THEN
	        CALL ST_ILST ( vgsub(2), ' ', 0, 1, ksubtp, knum, ier )
		IF ( ksubtp .lt. 1 .or. ksubtp .gt. 26 ) THEN
		    ksubtp = 0
		    kvgtyp = 0
		ELSE
		    kvgtyp = 20
		END IF
	      ELSE
		kvgtyp = 0
		ksubtp = 0
	    END IF
C
C*	    Read in the VGF file
C
            CALL ST_NULL ( cntrfl, cntrfl, lens, ier )
            CALL GG_RDVGF ( ggvgf, cntrfl, catmap, ier )
C
	    IF  ( ier .ne. 0 )  THEN
	         CALL ER_WMSG ( 'GRPHGD', ier, ' ', ier2 )
	         iret = -1
	         RETURN
	    END IF
C
C*	    Read in contour line information into common
C
	    CALL GGRDINFO ( cntrfl, keycol, kvgtyp, ksubtp, iolkdy, ier)
C
            IF ( ier .ne. 0 ) THEN
                CALL ER_WMSG ( 'GRPHGD', ier, cntrfl, ier2 )
		iret = -1
                RETURN
            END IF
C
C*          Create a surface file with vector information
C
	    CALL ST_LCUC ( type, type, ier )
  	    filnam = 'g2g_' // gdatim(:6) // gdatim(8:9) //
     +               '.sfc'
	    IF ( type(1:1) .eq. 'A' .or. type(1:1) .eq. 'B' )
     +	        CALL GGOACR ( filnam, gdatim(:11), type, ier )
C
	    IF ( type(1:1) .ne. 'A' .and. type(1:1) .ne. 'B' )  THEN
C
		IF ( hstflg )
     &		    WRITE(6,*) 'Total number of lines = ', nlines,
     &		    ' (',nint((nlines*100.0)/MAXLIN),'% of MAX ALLOWED)'
C
	        IF ( nlines .eq. 0 )  THEN
C
C*		    Set data array to default value and write to the file.
C
		    DO  ij = 1, kx*ky
		        grid(ij) = defalt
		    END DO
	            WRITE(6,*) 'Set data array to default = ', defalt
C
		    parm = gparm(:12)

                    CALL DG_NWDT ( grid, gdattm, level, ivcord, parm,
     +                             ighdr, gpack, rewrit, ier )
C
		    proces = ier .eq. 0
		    IF ( hstflg .and. proces )  THEN
		        parm = 'hist' // gparm(:8)
                        CALL DG_NWDT ( hist, gdattm, level, ivcord,
     +                             parm, ighdr, gpack, rewrit, ier )
		    END IF
		    proces = .false.
C
		END IF
C
	    END IF
C
	END IF
C
	IF ( proces )  THEN
	    bnds = bounds
C
C*	    Perform BOUNDS SEARCH
C
            CALL GBOUND ( imn, imx, jmn, jmx, grid, hist, kx, ky, ier )
            IF ( hstflg )
     +        WRITE (6,*) 'BOUNDS SEARCH COMPLETE, ier =', ier
C
	    IF ( type(1:1) .eq. 'A' .or. type(1:1) .eq. 'B' ) THEN
C
C*	        Process vector.
C
	        CALL IP_PUTV ( 'SFFILE', filnam, ier )
	        CALL IP_PUTV ( 'GDFILE', gdfile, ier )
	        CALL IP_PUTV ( 'SFPARM', 'UWND;VWND', ier )
	        CALL IP_PUTV ( 'DATTIM', gdatim(:11), ier )
	        CALL IP_PUTV ( 'DTAAREA', '-90;-180;90;180', ier )
C
 	        CALL OABSDR ( work1, work2, work3, data, infoflg, ier )

	        hstflg = .false.
C
C*	        Delete temporary surface data file.
C
		CALL ST_NULL ( filnam, filnam, len, ier )
		CALL GGOADL ( filnam, ier )
C
 	    ELSE
C
C*	       Process scalar.
C*	       Scan contour lines for all intersections; put into common.
C
	       CALL GGEDGE ( edgeopts, ier )
	       CALL GGBNDCONDS ( kx, ky, ier )
	       CALL GGGINT ( grid, hist, kx, ky, ier )
               IF ( hstflg )  THEN
                  WRITE (6,*)'INTERSECTIONS SEARCH COMPLETE, ier =', ier
		  numints = 0
		  DO  jj = 1, 4
		    ntot = 0
		    DO  ii = 1, intdim(jj)
		    DO  nn = 1, intptrs(ii,jj,1)
		        ntot = ntot + 1
			numints = numints + 1
		    END DO
		    END DO
		    WRITE(6,*) 'Total number of ints in dir ',jj,
     &			' = ',ntot
		  END DO
		  WRITE(6,*)'TOTAL NUMBER of INTERSECTIONS = ', numints,
     &		'  (',nint((numints*100.0)/MAXINT),'% of MAX ALLOWED)'
 	       END IF
C
C*	    Check whether it is discrete or continuous
C
	    CALL IN_DISCQ ( istate, ier )
C
C
C*	    The following is for discrete case
C
	    IF ( istate .eq. 1 ) THEN
C
C*	      Perform DLINE SEARCH
C
	      CALL GDLINE ( imn, imx, jmn, jmx, grid, hist, kx, ky, ier )
	      IF ( hstflg )
     +	        WRITE (6,*) 'DLINE SEARCH COMPLETE, ier =', ier
C
C*	      Perform DLINE WEIGHTED SEARCH
C
	      CALL GDLWEI ( imn, imx, jmn, jmx, grid, hist, kx, ky, ier )
	      IF ( hstflg )
     +	        WRITE (6,*) 'DLINE WEIGHTED SEARCH COMPLETE, ier =', ier
C
	    ELSE
C
C*	      The following is for continuous case
C
C*	       Perform RADIAL SEARCH
C
	       CALL GRADLS (imn, imx, jmn, jmx, grid, hist, kx, ky, ier)
	       IF ( hstflg )
     +	         WRITE (6,*) 'RADIAL SEARCH COMPLETE, ier =', ier
C
C*	       Perform LIKE-VALUE SEARCH
C
	      CALL GLVALU ( imn, imx, jmn, jmx, grid, hist, kx,
     +							ky, ier )
	       IF ( hstflg )
     +	         WRITE (6,*) 'LIKE-VALUE SEARCH COMPLETE, ier =', ier
C
C*	      Check the state of DLINES. If yes, add calls to
C*            GDLINE and GWEIGS
C
	      CALL IN_DLINQ ( istate, istateL, epsi, ier )
C
	      IF ( istate .eq. 1 ) THEN
C
C*	        Perform DLINE SEARCH
C
	        CALL GDLINE ( imn, imx, jmn, jmx, grid, hist, kx,
     +							ky, ier )
	        IF ( hstflg )
     +	          WRITE (6,*) 'DLINE SEARCH COMPLETE, ier =', ier
C
	      END IF
C
C*	       Perform WEIGHTED SEARCH
C
 	       CALL GWEIGS (imn, imx, jmn, jmx, grid, hist, kx, ky, ier)
 	       IF ( hstflg )
     +	         WRITE (6,*) 'WEIGHTED SEARCH COMPLETE, ier =', ier
C
C*	        Perform SMOOTHING
C
 	        CALL GSMOOT (imn, imx, jmn, jmx, grid, hist, kx, ky, ier)
 	        IF ( hstflg )
     +	          WRITE (6,*) 'SMOOTHING COMPLETE, ier =', ier
C
C*	        Apply limits, if any.
C
 	        CALL GGLIMS (imn, imx, jmn, jmx, grid, hist, kx, ky, ier)
 	        IF ( hstflg )
     +	          WRITE (6,*) 'LIMITS COMPLETE, ier =', ier
C
C*	      This is the end of continous case
C
	    ENDIF
C
C*	        Write computed grid and history grid to grid file.
C
	        parm = gparm(:12)
                CALL DG_NWDT ( grid, gdattm, level, ivcord, parm,
     +                         ighdr, gpack, rewrit, ier )
	        IF ( hstflg )  WRITE (6,*) 'RETURN from DG_NWDT =', ier
		IF  ( ier .ne. 0 )  THEN
	    	    CALL ER_WMSG ( 'DG', ier, ' ', ier )
		    iret = -13
	    	    RETURN
		END IF
	        proces = ier .eq. 0
C
C*	        Figure whether history grid needs to be written.
C
                IF  ( hstflg .and. proces )  THEN
C
	          parm = 'hist' // gparm(:8)
                  CALL DG_NWDT ( hist, gdattm, level, ivcord, parm,
     +                           ighdr, gpack, rewrit, ier )
	          IF ( hstflg )  WRITE(6,*) 'RETURN from DG_NWDT =', ier
		  IF  ( ier .ne. 0 )  THEN
	    	    CALL ER_WMSG ( 'DG', ier, ' ', ier )
		    iret = -13
	    	    RETURN
		  END IF
C
	        END IF

	    END IF
C
	END IF
C
        IF  ( hstflg .and. proces )  THEN
C
	    nmiss = 0
	    grdmin =  1.0E10
	    grdmax = -1.0E10
	    ngpts = (imx-imn+1)*(jmx-jmn+1)
	    DO  ii = 1, ngpts
		IF ( grid(ii) .ne. RMISSD )  THEN
		    grdmin = min(grdmin,grid(ii))
		    grdmax = max(grdmax,grid(ii))
		ELSE
		    nmiss = nmiss + 1
		END IF
	    END DO
	    WRITE(6,*) 'Total number of grid points = ', ngpts
	    WRITE(6,*) 'Grid minimum and maximum = ', grdmin, grdmax
	    WRITE(6,*) 'Number of missing values = ', nmiss
C
	END IF
C
	CALL DG_NEND ( ier )
C
	RETURN
C*
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
program if_stop
  if (.true.) stop  ! assertion failure
  if (.true.) stop 2  ! 2 becomes "2" in translation
  if (n .ne. 2) then
    stop
  endif
  if (.true.) then
    stop 2
  endif
  if (.false.) then
    n = 1
  else
    stop
  endif
  if (n .eq. 1) then
    stop 'Done #123'
  else
    stop 2
  endif
end program

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
!---------------------------------------------------------------
!   This subroutine checks the number of external dipoles----
!---absorbing the correct number into the fragmenation functions
!   it then returns the finite (msq_qcd*dip) ---------------------
!---------------------------------------------------------------
!--- Author C. Williams April 2012
!-----------------------------------------------------------------
!---- modified March 2013 during frag re-write

      subroutine qqb_zaa_fragdips(p,p_phys,qcd_tree,msq_out)
      implicit none
      include 'constants.f'
      include 'ewcouple.f'
      include 'ewcharge.f'
      include 'frag.f'
      double precision p(mxpart,4),p_phys(mxpart,4)
      double precision msq_qcd(-nf:nf,-nf:nf),msq_out(-nf:nf,-nf:nf)
      integer j,k
      double precision virt_dips,xl,dot,fsq
      double precision aewo2pi,fi_gaq
      external qcd_tree

      aewo2pi=esq/(fourpi*twopi)

      fsq=frag_scale**2

      xl=dlog(-two*dot(p_phys,2,6)/fsq)
      virt_dips=+aewo2pi*(fi_gaq(z_frag,p_phys,xl,6,2,2))

      do j=-nf,nf
         do k=-nf,nf
            msq_qcd(j,k)=0d0
            msq_out(j,k)=0d0
         enddo
      enddo

      call qcd_tree(p,msq_qcd)

      do j=-nf,nf
         do k=-nf,nf

!   factor of two cancelled by statistical factor because two photons
            if((j.eq.0).and.(k.ne.0)) then
                  msq_out(j,k)=msq_qcd(j,k)*Q(k)**2*virt_dips
            elseif((j.ne.0).and.(k.eq.0)) then
                  msq_out(j,k)=msq_qcd(j,k)*Q(j)**2*virt_dips
            endif

         enddo
      enddo

      return
      end subroutine

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


module cvode_interface
  use fcvode_mod
  use fsunmat_dense_mod
  use fsunlinsol_dense_mod
  contains
  integer(c_int) function FCVDense(cvode_mem, N)
    use, intrinsic :: iso_c_binding
    use fnvector_serial
    implicit none
    type(c_ptr),     value :: cvode_mem
    integer(c_long), value :: N
    type(c_ptr)            :: sunmat_A
    type(c_ptr)            :: sunlinsol_LS
    type(c_ptr)            :: sunvec_y
    integer(c_int)         :: ierr

    sunvec_y = N_VNewEmpty_Serial(N)
    sunmat_A = FSUNDenseMatrix(N, N)
    sunlinsol_LS = FSUNDenseLinearSolver(sunvec_y, sunmat_A)
    ierr = FCVDlsSetLinearSolver(cvode_mem, sunlinsol_LS, sunmat_A)
  end function FCVDense
end module cvode_interface


text:
      subroutine dgthrz ( nz, y, x, indx )
C
C     ==================================================================
C     ==================================================================
C     ====  DGTHRZ -- DOUBLE PRECISION GATHER AND ZERO              ====
C     ==================================================================
C     ==================================================================
C
C     PURPOSE
C     -------
C
C         DGTHRZ GATHERS THE SPECIFIED ELEMENTS FROM
C             A DOUBLE PRECISION VECTOR  Y  IN FULL STORAGE FORM
C         INTO
C             A DOUBLE PRECISION VECTOR  X  IN COMPRESSED FORM  (X,INDX).
C         FURTHERMORE THE GATHERED ELEMENTS OF  Y  ARE SET TO ZERO.
C
C         ONLY THE ELEMENTS OF  Y  WHOSE INDICES ARE LISTED IN  INDX
C         ARE REFERENCED OR MODIFIED.
C
C     ARGUMENTS
C     ---------
C
C     INPUT ...
C
C         NZ      INTEGER     NUMBER OF ELEMENTS TO BE GATHERED INTO
C                             COMPRESSED FORM.
C         INDX    INTEGER     ARRAY CONTAINING THE INDICES OF THE
C                             VALUES TO BE GATHERED INTO COMPRESSED FORM.
C
C     UPDATED ...
C
C         Y       DOUBLE      ARRAY, ON INPUT, WHICH CONTAINS THE
C                             VECTOR  Y  IN FULL STORAGE FORM.  THE
C                             GATHERED COMPONENTS IN  Y  ARE SET TO ZERO.
C                             ONLY THE ELEMENTS CORRESPONDING TO THE
C                             INDICES IN  INDX  HAVE BEEN ACCESSED.
C
C     OUTPUT ...
C
C         X       DOUBLE      ARRAY CONTAINING THE VALUES GATHERED INTO
C                             THE COMPRESSED FORM.
C
C     SPARSE BASIC LINEAR ALGEBRA SUBPROGRAM
C
C     FORTRAN VERSION WRITTEN OCTOBER 1984
C     ROGER G GRIMES, BOEING COMPUTER SERVICES
C
C     ==================================================================
C
C     -------------
C     ... ARGUMENTS
C     -------------
C
      INTEGER             NZ, INDX (*)
C
      DOUBLE PRECISION    Y (*), X (*)
C
C     -------------------
C     ... LOCAL VARIABLES
C     -------------------
C
      INTEGER             I
C
C     ==================================================================
C
      IF  ( NZ .LE. 0 )  RETURN
C
      DO 10 I = 1, NZ
          X(I)       = Y(INDX(I))
          Y(INDX(I)) = 0.0D0
   10 CONTINUE
C
      RETURN
      END


text:
MODULE cartesian_mod
  USE task_list_mod
  USE experiment_mod
  USE patch_mod
  USE link_mod
  USE task_mod
  USE trace_mod
  USE io_mod
  USE timer_mod
  USE mpi_coords_mod
  USE mpi_mod
  USE omp_mod
  USE omp_timer_mod
  USE bits_mod
  implicit none
  private
  type, public, extends(task_t):: rank_t
  end type
  type, public:: cartesian_t
    real(8):: size(3)
    integer:: dims(3)
    type(task_list_t), pointer:: task_list
  contains
    procedure init
    procedure, nopass:: diagnostics
  end type
  logical:: omp_init=.false.
  type(rank_t):: rank
  integer, save:: np=0                  ! patches per process
  public task_list_t                    ! export, to avoid excessive USE
CONTAINS

!===============================================================================
!> Distribute patches in a Cartesian arrangement in a box
!===============================================================================
SUBROUTINE init (self, label)
  class(cartesian_t):: self
  character(len=*), optional:: label
  !.............................................................................
  type(task_list_t), pointer:: patch_list
  class(link_t), pointer:: link, next
  class(experiment_t), pointer:: patch
  class(task_t), pointer:: task
  integer:: i, j, k, id, ip
  integer:: dims(3) = [4,4,4]
  integer:: patches_per_mpi(3)
  integer:: mpi_dims(3) = [1,1,1]
  integer:: per_rank(3) = [0,0,0]
  real(8):: size(3) = [1.0,1.0,1.0]
  real(8):: origin(3) = [-0.5,-0.5,-0.5]
  real(8):: wc
  logical:: face_nbors=.false.
  logical:: periodic(3)=.true.
  integer, save:: itimer=0
  namelist /cartesian_params/ size, dims, mpi_dims, per_rank, origin, face_nbors, &
    omp_init, periodic
  character(len=120):: ids = &
    '$Id: e71d387b943008f46d0f441b4ecad422725030df $ components/cartesian_mod.f90'
  !-----------------------------------------------------------------------------
  call trace%begin ('cartesian_t%init')
  call io%header('begin cartesian_t%init: Cartesian patch arrangement')
  call trace%print_id (ids)
  !-----------------------------------------------------------------------------
  rewind (io%input); read (io%input, cartesian_params)
  write (io%output, cartesian_params)
  if (any(per_rank/=0)) then
    where (per_rank/=0) dims=per_rank*mpi_dims
  end if
  per_rank = dims/mpi_dims
  self%dims = dims
  self%size = size
  io%dims     = dims
  if (mpi%size==1) &
    mpi_dims = 1
  !-----------------------------------------------------------------------------
  ! Generate a (periodic) MPI geometry, with dims specified in the input file,
  ! or assigned automatically by mpi_coords%init.
  !-----------------------------------------------------------------------------
  if (mpi%size == 1) then
     mpi_coords%dims = 1
  else
    if (product(mpi_dims)==1) then
      call mpi_coords%init (dims=dims)
    else
      call mpi_coords%init (mpi_dims=mpi_dims)
    end if
  end if
  mpi_dims = mpi_coords%dims
  patches_per_mpi = self%dims/mpi_dims
  mpi_coords%npatch = patches_per_mpi
  io%mpi_dims = mpi_dims
  !-----------------------------------------------------------------------------
  ! Set attritbutes of the rank_t data type -- just enough to use distance().
  !-----------------------------------------------------------------------------
  rank%periodic = .true.
  rank%box = self%size
  rank%size = rank%box/mpi_dims
  rank%ds = rank%size/patches_per_mpi
  rank%position = (mpi_coords%rank_to_coords(mpi%rank) + 0.5_8)*rank%size
  !-----------------------------------------------------------------------------
  ! Begin on task list
  !-----------------------------------------------------------------------------
  allocate (self%task_list)
  call self%task_list%init ('tlist')
  self%task_list%size = self%size
  self%task_list%dims = self%dims
  io%dims = self%dims
  self%task_list%n_tasks = product(self%dims)
  self%task_list%face_nbors = face_nbors
  !-----------------------------------------------------------------------------
  ! Allocate an array of links, pointing to patches, and give the patches each
  ! id, rank, size, position, integer position, and status bits.  Discard
  ! external patches and add the rest to the task_list.
  !-----------------------------------------------------------------------------
  ip = 0
  do k=0,self%dims(3)-1
  do j=0,self%dims(2)-1
  do i=0,self%dims(1)-1
    allocate (link, patch); link%task => patch; patch%link => link
    task => link%task
    ! --- patch 3D and 1D coordinates ---
    patch%ipos = [i,j,k]
    patch%rank = mpi_coords%coords_to_rank ([i,j,k]/patches_per_mpi)
    if (patch%rank == mpi%rank) then
      ip = ip + 1
      patch%ip = ip
    end if
    !---------------------------------------------------------------------------
    ! Generate the same IDs on all ranks, pruning the external tasks below
    !---------------------------------------------------------------------------
    call patch%task_t%init_unique (same=.true.)
    patch%box = self%size
    patch%size = self%size / self%dims
    patch%position = ([i,j,k]+0.5d0)*patch%size + origin
    if (any(abs(rank%distance (patch)) > 0.5*(rank%size + 1.5*rank%ds))) then
      deallocate (link, patch)
      cycle
    end if
    call patch%set (bits%virtual)
    patch%ds = patch%size*0.001
    patch%origin = origin
    patch%llc_cart = patch%position - 0.5 * patch%size
    patch%centre_nat = patch%position
    patch%llc_nat = patch%llc_cart
    patch%periodic = periodic
    call patch%set(bits%static)
    if (all(self%dims==1)) call task%set(bits%root_grid)
    ! -- initially append all tasks to the task list, to build nbor lists ------
    call self%task_list%append_link (link)
  end do
  end do
  end do
  !-----------------------------------------------------------------------------
  ! Generate nbor lists and set status bits for all tasks
  !-----------------------------------------------------------------------------
  if (io%verbose > 0) write(stdout,*) 'cartesian_t%init: init_all_nbors'
  call self%task_list%init_all_nbors
  if (io%verbose > 0) write(stdout,*) 'cartesian_t%init: reset_status'
  call self%task_list%reset_status
  !-----------------------------------------------------------------------------
  ! Remove tasks that are external, and make sure to also remove them from the
  ! nbor lists of virtual tasks
  !-----------------------------------------------------------------------------
  link => self%task_list%head
  do while (associated(link))
    next => link%next
    if (link%task%is_set (bits%external)) then
      call link%remove_nbor_list2 (link%nbor)
      call self%task_list%remove (link)
      deallocate (link%task)
      deallocate (link)
    end if
    link => next
  end do
  !-----------------------------------------------------------------------------
  ! This call is not valid until the patch meshes have been allocated
  !-----------------------------------------------------------------------------
  if (io%verbose > 0) write(stdout,*) 'cartesian_t%init: count status'
  call self%task_list%count_status
  !-----------------------------------------------------------------------------
  ! These must be set before reading input snapshots
  !-----------------------------------------------------------------------------
  io%ntask = self%task_list%na
  io%nwrite = self%task_list%na
  io%ntotal = product(self%dims)
  write (stdout,*) 'Number of tasks in task list:', self%task_list%n
  write(stdout,*) 'ntask, nwrite, ntotal =', io%ntask, io%nwrite, io%ntotal
  !-----------------------------------------------------------------------------
  ! Initialize tasks
  !-----------------------------------------------------------------------------
  if (omp_init) then
    !$omp parallel default(shared)
    call init_exp (self, origin)
    !$omp end parallel
  else
    call init_exp (self, origin)
  end if
  !-----------------------------------------------------------------------------
  ! This call is not valid until the patch mem and meshes have been allocated
  !-----------------------------------------------------------------------------
  if (io%verbose > 0) &
    write(stdout,*) 'cartesian_t%init: init boundaries'
  call self%task_list%init_bdries
  !-----------------------------------------------------------------------------
  if (io_unit%master) then
    write (io_unit%nml, cartesian_params)
    flush (io_unit%nml)
  end if
  if (io%verbose > 0) &
    write(stdout,*) 'cartesian_t%init: diagnostics'
  call self%diagnostics()
  if (io%verbose > 0) &
    write(stdout,*) 'cartesian_t%init: print tasks'
  call self%task_list%print_tasks
  if (io%verbose > 0) &
    write(stdout,*) 'cartesian_t%init: done'
  flush (stdout)
  !-----------------------------------------------------------------------------
  call trace%end ()
END SUBROUTINE init

SUBROUTINE init_exp (self, origin)
  class(cartesian_t):: self
  real(8):: origin(3)
  !.............................................................................
  type(link_t), pointer:: link
  class(task_t), pointer:: task
  class(experiment_t), pointer:: patch
  integer:: id
  !-----------------------------------------------------------------------------
  id = 0
  link => self%task_list%head
  do while (associated(link))
    id = id+1
    if (.not.omp_init .or. mod(id,omp%nthreads)==omp%thread) then
      task => link%task
      select type (task)
      class is (experiment_t)
        patch => task
        call patch%init
        if (omp_init .and. omp%nthreads>1) then
          task%mem_thread = omp%thread
        else
          task%mem_thread = -1
        end if
        patch%mesh%origin = origin
        !$omp atomic
        np = np + 1
        ! Set physical boundary bits if required.
        if (.not.patch%periodic(1).and.patch%ipos(1)==0)              call patch%boundaries%set(bits%xl)
        if (.not.patch%periodic(2).and.patch%ipos(2)==0)              call patch%boundaries%set(bits%yl)
        if (.not.patch%periodic(3).and.patch%ipos(3)==0)              call patch%boundaries%set(bits%zl)
        if (.not.patch%periodic(1).and.patch%ipos(1)==self%dims(1)-1) call patch%boundaries%set(bits%xu)
        if (.not.patch%periodic(2).and.patch%ipos(2)==self%dims(2)-1) call patch%boundaries%set(bits%yu)
        if (.not.patch%periodic(3).and.patch%ipos(3)==self%dims(3)-1) call patch%boundaries%set(bits%zu)
        if (io%verbose>0) then
          write (io_unit%log,1) &
            'cartesian_t%init: id, rank, size, pos =', &
            patch%id, patch%rank, patch%size, patch%position, &
            patch%is_set(bits%boundary), patch%is_set(bits%virtual)
          1 format(a,i8,i6,1p,2(2x,3e14.6),2x,2l1,2x,a)
          flush (io_unit%log)
        end if
        if (patch%id == io%id_debug) then
          write (io_unit%output,1) &
            'cartesian_t%init: id, rank, size, pos =', &
            patch%id, patch%rank, patch%size, patch%position, &
            patch%is_set(bits%boundary), patch%is_set(bits%virtual), 'DBG'
          flush (io_unit%log)
        else if (io_unit%verbose>1 .and. io_unit%master) then
          write (io_unit%output,1) &
            'cartesian_t%init: id, rank, size, pos =', &
            patch%id, patch%rank, patch%size, patch%position, &
            patch%is_set(bits%boundary), patch%is_set(bits%virtual)
          flush (io_unit%log)
        end if
      end select
    end if
    link => link%next
  end do
END SUBROUTINE init_exp

!===============================================================================
!> Diagnostics: Number of patches and GB of mem
!===============================================================================
SUBROUTINE diagnostics
  if (mpi%master) then
    print '(a,i8)',     ' Patches per process:', np
  end if
END SUBROUTINE diagnostics

END MODULE cartesian_mod

max_stars_repo_name: applejwjcat/dispatch

max_issues_repo_name: applejwjcat/dispatch

max_forks_repo_name: applejwjcat/dispatch


text:
*DECK DCSEVL
      DOUBLE PRECISION FUNCTION DCSEVL (X, CS, N)
C***BEGIN PROLOGUE  DCSEVL
C***PURPOSE  Evaluate a Chebyshev series.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C3A2
C***TYPE      DOUBLE PRECISION (CSEVL-S, DCSEVL-D)
C***KEYWORDS  CHEBYSHEV SERIES, FNLIB, SPECIAL FUNCTIONS
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C  Evaluate the N-term Chebyshev series CS at X.  Adapted from
C  a method presented in the paper by Broucke referenced below.
C
C       Input Arguments --
C  X    value at which the series is to be evaluated.
C  CS   array of N terms of a Chebyshev series.  In evaluating
C       CS, only half the first coefficient is summed.
C  N    number of terms in array CS.
C
C***REFERENCES  R. Broucke, Ten subroutines for the manipulation of
C                 Chebyshev series, Algorithm 446, Communications of
C                 the A.C.M. 16, (1973) pp. 254-256.
C               L. Fox and I. B. Parker, Chebyshev Polynomials in
C                 Numerical Analysis, Oxford University Press, 1968,
C                 page 56.
C***ROUTINES CALLED  D1MACH, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770401  DATE WRITTEN
C   890831  Modified array declarations.  (WRB)
C   890831  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   900329  Prologued revised extensively and code rewritten to allow
C           X to be slightly outside interval (-1,+1).  (WRB)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  DCSEVL
      DOUBLE PRECISION B0, B1, B2, CS(*), ONEPL, TWOX, X, D1MACH
      LOGICAL FIRST
      SAVE FIRST, ONEPL
      DATA FIRST /.TRUE./
C***FIRST EXECUTABLE STATEMENT  DCSEVL
      IF (FIRST) ONEPL = 1.0D0 + D1MACH(4)
      FIRST = .FALSE.
      IF (N .LT. 1) CALL XERMSG ('SLATEC', 'DCSEVL',
     +   'NUMBER OF TERMS .LE. 0', 2, 2)
      IF (N .GT. 1000) CALL XERMSG ('SLATEC', 'DCSEVL',
     +   'NUMBER OF TERMS .GT. 1000', 3, 2)
      IF (ABS(X) .GT. ONEPL) CALL XERMSG ('SLATEC', 'DCSEVL',
     +   'X OUTSIDE THE INTERVAL (-1,+1)', 1, 1)
C
      B1 = 0.0D0
      B0 = 0.0D0
      TWOX = 2.0D0*X
      DO 10 I = 1,N
         B2 = B1
         B1 = B0
         NI = N + 1 - I
         B0 = TWOX*B1 - B2 + CS(NI)
   10 CONTINUE
C
      DCSEVL = 0.5D0*(B0-B2)
C
      RETURN
      END

max_stars_repo_name: floryst/ITK

max_issues_repo_name: floryst/ITK

max_forks_repo_name: floryst/ITK