text:
This is a list of the oldest, continuously operated organizations or entities in the Davis, California area.

 1. 1864. The Tremont Mite Society, founded circa 1864, celebrated its Centennial in 1964. Contact: Ruth Bink Eason, 7076781249. The original purpose of the organization was to raise money to build a church at Tremont, a few miles south of Davis. Since that achievement in 1871, it has raised money and volunteer labor for church upkeep. The The Tremont Church and Cemetery church is still there and the task is ongoing (Larkey, Davisville ’68, p. 105).
 1. 1867. Post Office The United States Postal Service, local office in 1867. The famous William (“Solano Bill”) Dresbach, got his start as a local postmaster operating out of Solano House on the south bank of Putah Creek across from the bend of the Putah Creek Lodge Road. When he got wind a railroad was to come his way and create a grainshipping village, he illegally moved his post office from Solano to Yolo County in early 1868. The Jerome Davis farm had longsince gone bankrupt and its farm house was empty. Dresbach installed the post office in that building (which stood near the middle of the parking lot at A and First streets) and named it Yolo House (Lofland, Davisville 1868, Part I).
 1. 18681. Davis Joint Unified School District The Davis Public Schools began September, 1868 when a group of citizens petitioned the Board of Supervisors to authorize a school, which it did (Vaught, After the Gold Rush, p. 137). It was called Yolo School and joined six other alreadyexisting public schools “organized in rural districts in the 1860s” (Larkey, p. 83). (Some of these other six schools might also be counted as stillexisting and be added to this list of oldest organizations.)
 1. 18682: The Davis Democratic Club was active in 1868 as the Democratic Club of Putah Township (Vaught, p. 135, Lofland, placard 61768).
 1. 18691: The Davis Community Church “was authorized on May 8, 1869,” as The Davisville Presbyterian Church (Larkey, p. 105).
 1. 18692: A St. James Catholic Church Catholic Church group assembled in 1869 (Larkey, p. 107).
 1. 1870: Yolo Lodge, No. 169, The Independent Order of Odd Fellows was founded in ceremonies at Yolo House on April 12, 1870  (Vaught, p. 137).

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
       CHARACTER*42 FUNCTION ERRORWORDS(IT)	!Look for an explanation. One day, the system may offer coherent messages.
Curious collection of encountered codes. Will they differ on other systems?
Compaq's compiler was taken over by unintel; http://software.intel.com/sites/products/documentation/hpc/compilerpro/en-us/fortran/lin/compiler_f/bldaps_for/common/bldaps_rterrs.htm
contains a schedule of error numbers that matched those I'd found for Compaq, and so some assumptions are added.
Copying all (hundreds!) is excessive; these seem possible for the usage so far made of error diversion.
Compaq's compiler interface ("visual" blah) has a help offering, which can provide error code information.
Compaq messages also appear in http://cens.ioc.ee/local/man/CompaqCompilers/cf/dfuum028.htm#tab_runtime_errors
Combines IOSTAT codes (file open, read etc) with STAT codes (allocate/deallocate) as their numbers are distinct.
Completeness and context remains a problem. Excess brevity means cause and effect can be confused.
        INTEGER IT			!The error code in question.
        INTEGER LASTKNOWN 		!Some codes I know about.
        PARAMETER (LASTKNOWN = 26)	!But only a few, discovered by experiment and mishap.
        TYPE HINT			!For them, I can supply a table.
         INTEGER	CODE		!The code number. (But, different systems..??)
         CHARACTER*42	EXPLICATION	!An explanation. Will it be te answer?
        END TYPE HINT			!Simple enough.
        TYPE(HINT) ERROR(LASTKNOWN)	!So, let's have a collection.
        PARAMETER (ERROR = (/		!With these values.
     1   HINT(-1,"End-of-file at the start of reading!"),	!From examples supplied with the Compaq compiler involving IOSTAT.
     2   HINT( 0,"No worries."),			!Apparently the only standard value.
     3   HINT( 9,"Permissions - read only?"),
     4   HINT(10,"File already exists!"),
     5   HINT(17,"Syntax error in NameList input."),
     6   HINT(18,"Too many values for the recipient."),
     7   HINT(19,"Invalid naming of a variable."),
     8   HINT(24,"Surprise end-of-file during read!"),	!From example source.
     9   HINT(25,"Invalid record number!"),
     o   HINT(29,"File name not found."),
     1   HINT(30,"Unavailable - exclusive use?"),
     2   HINT(32,"Invalid fileunit number!"),
     3   HINT(35,"'Binary' form usage is rejected."),	!From example source.
     4   HINT(36,"Record number for a non-existing record!"),
     5   HINT(37,"No record length has been specified."),
     6   HINT(38,"I/O error during a write!"),
     7   HINT(39,"I/O error during a read!"),
     8   HINT(41,"Insufficient memory available!"),
     9   HINT(43,"Malformed file name."),
     o   HINT(47,"Attempting a write, but read-only is set."),
     1   HINT(66,"Output overflows single record size."),	!This one from experience.
     2   HINT(67,"Input demand exceeds single record size."),	!These two are for unformatted I/O.
     3   HINT(151,"Can't allocate: already allocated!"),	!These different numbers are for memory allocation failures.
     4   HINT(153,"Can't deallocate: not allocated!"),
     5   HINT(173,"The fingered item was not allocated!"),	!Such as an ordinary array that was not allocated.
     6   HINT(179,"Size exceeds addressable memory!")/))
        INTEGER I		!A stepper.
         DO I = LASTKNOWN,1,-1	!So, step through the known codes.
           IF (IT .EQ. ERROR(I).CODE) GO TO 1	!This one?
         END DO			!On to the next.
    1    IF (I.LE.0) THEN	!Fail with I = 0.
           ERRORWORDS = I8FMT(IT)//" is a novel code!"	!Reveal the mysterious number.
          ELSE			!But otherwise, it is found.
           ERRORWORDS = ERROR(I).EXPLICATION	!And these words might even apply.
         END IF			!But on all systems?
       END FUNCTION ERRORWORDS	!Hopefully, helpful.

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


max_stars_repo_name: sriyerg/systemverilog

max_issues_repo_name: sriyerg/systemverilog

max_forks_repo_name: sriyerg/systemverilog


text:
      SUBROUTINE BUNPAK (IG,I,NJ,JG)
C
C     THIS ROUTINE WORKS SIMILARLY AS BUNPK EXCEPT IT UNPACKS A WHOLE
C     (I-TH) ROW OF GRID NUMBERS (1 THRU NJ) FROM IG TABLE, AND SAVES
C     THE UNPACKED DATA IN JG ARRAY.
C     (BUNPK UNPACKS ONLY AN ELEMENT OF GRID NUMBER IN IG TABLE)
C
C     THIS ROUTINE GREATLY INCREASES BANDIT INTERNAL EFFICIENCY
C     WRITTEN BY G.CHAN/UNISYS,    MAY 1988
C
      IMPLICIT INTEGER (A-Z)
C
CDC   NEXT 2 LINES FOR CDC AND UNIVAC ONLY
C     EXTERNAL         ANDF,     RSHIFT
C     INTEGER          ANDF,     RSHIFT  ,IG(1)
C
C     NEXT LINE FOR IBM, VAX, AND MACHINES THAT HAVE INTEGER*2
      INTEGER*2        IG(1)
C
      INTEGER          JG(1),    NAM(2)
      COMMON /SYSTEM/  IBUF,     NOUT
      COMMON /BANDB /  NBIT,     DUM3B(3), IPASS,    NW
      COMMON /BANDS /  DUM4S(4), II1,      MAXDEG,   DUM4(4),  MASK
      DATA    NAM   /  4HUNPA  , 4HK       /
C
      IF (NJ .LE. MAXDEG) GO TO 20
      WRITE  (NOUT,10) NJ,MAXDEG
 10   FORMAT ('0 *** BUNPAK .GT. MAXDEG',2I7)
      CALL ERRTRC (NAM)
C
 20   IPASS = IPASS+NJ
      N1 = I
C
C     ********************************************
C     UNIVAC AND CDC MACHINES
C     ********************************************
C
C     DO 40 N=1,NJ,NW
C     N2 = IG(N1)
C     N3 = N+NW-1
C     DO 30 M=1,NW
C     JG(N3) = ANDF(N2,MASK)
C     IF (M .EQ. NW) GO TO 40
C     N2 = RSHIFT(N2,NBIT)
C  30 N3 = N3-1
C  40 N1 = N1+II1
C     RETURN
C
C     ********************************************
C     IBM AND VAX MACHINES
C     ********************************************
C
      DO 50 N=1,NJ
      JG(N) = IG(N1)
   50 N1 = N1+II1
      RETURN
      END


text:
cdis    Forecast Systems Laboratory
cdis    NOAA/OAR/ERL/FSL
cdis    325 Broadway
cdis    Boulder, CO     80303
cdis
cdis    Forecast Research Division
cdis    Local Analysis and Prediction Branch
cdis    LAPS
cdis
cdis    This software and its documentation are in the public domain and
cdis    are furnished "as is."  The United States government, its
cdis    instrumentalities, officers, employees, and agents make no
cdis    warranty, express or implied, as to the usefulness of the software
cdis    and documentation for any purpose.  They assume no responsibility
cdis    (1) for the use of the software and documentation; or (2) to provide
cdis     technical support to users.
cdis
cdis    Permission to use, copy, modify, and distribute this software is
cdis    hereby granted, provided that the entire disclaimer notice appears
cdis    in all copies.  All modifications to this software must be clearly
cdis    documented, and are solely the responsibility of the agent making
cdis    the modifications.  If significant modifications or enhancements
cdis    are made to this software, the FSL Software Policy Manager
cdis    (softwaremgr@fsl.noaa.gov) should be notified.
cdis
cdis
cdis
cdis
cdis
cdis
cdis

        subroutine get_ground_temperature(i4time,laps_cycle_time
     1                    ,imax,jmax,lat,lon,r_missing_data
     1                                      ,cvr_snow,t_sfc_k,t_gnd_k)

!       This routine returns ground temp under the assumption of clear skies
!       After this routine is called from clouds, it may be worth filtering
!       the band 8 temps to expand the cloudy areas slightly. The kernel
!       of the filter would be the coldest pixel of the nine neighbors

        character var*3,comment*125,ext*31,units*10

        integer imax,jmax         ! input
        integer i4time            ! input
        integer laps_cycle_time   ! input
        real lat(imax,jmax)       ! input
        real lon(imax,jmax)       ! input
        real r_missing_data       ! input
        real t_sfc_k(imax,jmax)   ! input
        real t_gnd_k(imax,jmax)   ! output
        real cvr_snow(imax,jmax)  ! local

        var = 'SC'
        ext = 'lm2'
        call get_laps_2d(i4time-laps_cycle_time,ext,var,units,comment
     1                            ,imax,jmax,cvr_snow,istat_cvr_snow)

        if(istat_cvr_snow .ne. 1)then    ! Try using snow cover field
            write(6,*)' Error reading snow cover'
            do i = 1,imax
            do j = 1,jmax
                cvr_snow(i,j) = r_missing_data
            enddo
            enddo
        endif

!       Calculate solar altitude and ground temperature
        do i = 1,imax
        do j = 1,jmax
            call solar_position(lat(i,j),lon(i,j),i4time,solar_alt
     1                                     ,solar_dec,solar_ha)

!           Calculate ground temperature (for now equate to sfc air temp)
!           Note that function t_ground_k needs to be moved to the lapslib
            t_gnd_k(i,j) = t_ground_k(t_sfc_k(i,j)
     1                          ,solar_alt,solar_ha
     1              ,solar_dec,lat(i,j),cvr_snow(i,j),r_missing_data
     1                          ,i,j,imax,jmax)

        enddo
        enddo

        return
        end

        function t_ground_k(t_sfc_k,solar_alt,solar_ha
     1    ,solar_dec,rlat,cvr_snow,r_missing_data,i,j,ni,nj)

        real high_alt,low_alt,corr_low,corr_high,corr_ramp

        solar_transit_alt = 90. - abs(rlat - solar_dec)

        if(solar_ha .lt. 0.)then ! morning
            high_alt = solar_transit_alt * .66
            low_alt  = solar_transit_alt * .36
!           low_alt  = high_alt - 11.
        else                     ! afternoon
            high_alt = solar_transit_alt * .79
            low_alt  = solar_transit_alt * .49
!           low_alt  = high_alt - 11.
        endif

        corr_high =+4.
        corr_low = -4.
        corr_ramp = (corr_high - corr_low) / (high_alt - low_alt)

!       Warmer at day, colder at night
        if(solar_alt .lt. low_alt)then
            corr = corr_low

        elseif(solar_alt .gt. high_alt)then
            corr = corr_high

        else ! ramp the function
            corr = corr_low + (solar_alt - low_alt) * corr_ramp

        endif

        if(i .eq. ni/2 .and. j .eq. nj/2)then
            write(6,*)' solar_ha,solar_alt,solar_transit_alt,ratio = '
     1             ,solar_ha,solar_alt,solar_transit_alt
     1             ,solar_alt/solar_transit_alt
            write(6,*)' high/low alt, corr = ',high_alt,low_alt,corr
        endif

        corr_negonly = min(corr,0.0)            ! Only colder at night

        if(cvr_snow .ne. r_missing_data)then
            t_ground_fullcorr = t_sfc_k + corr  ! Warmer at day, colder at night
            t_snow_k = t_sfc_k + corr_negonly   ! Only colder at night
            t_snow_k = min(t_snow_k,273.15)     ! Snow can't be above 0C
            t_ground_k = t_ground_fullcorr * (1.-cvr_snow) + t_snow_k *
     1cvr_snow
        else
            t_ground_k = t_sfc_k + corr_negonly ! Only colder at night
        endif

        return
        end

max_stars_repo_name: maxinye/laps-mirror

max_stars_repo_licenses: Intel Unlicense OLDAP-2.2.1 NetCDF

max_issues_repo_name: longwosion/laps-mirror

max_issues_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense

max_forks_repo_name: longwosion/laps-mirror

max_forks_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense


!
!
!                             AMG4PSBLAS version 1.0
!    Algebraic Multigrid Package
!               based on PSBLAS (Parallel Sparse BLAS version 3.7)
!
!    (C) Copyright 2021
!
!        Salvatore Filippone
!        Pasqua D'Ambra
!        Fabio Durastante
!
!    Redistribution and use in source and binary forms, with or without
!    modification, are permitted provided that the following conditions
!    are met:
!      1. Redistributions of source code must retain the above copyright
!         notice, this list of conditions and the following disclaimer.
!      2. Redistributions in binary form must reproduce the above copyright
!         notice, this list of conditions, and the following disclaimer in the
!         documentation and/or other materials provided with the distribution.
!      3. The name of the AMG4PSBLAS group or the names of its contributors may
!         not be used to endorse or promote products derived from this
!         software without specific written permission.
!
!    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!    TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!    PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS
!    BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!    CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!    SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!    INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!    POSSIBILITY OF SUCH DAMAGE.
!
!
!
!  The aggregator object hosts the aggregation method for building
!  the multilevel hierarchy.
!
module amg_z_base_aggregator_mod

  use amg_base_prec_type, only : amg_dml_parms, amg_daggr_data
  use psb_base_mod, only : psb_zspmat_type, psb_lzspmat_type, psb_z_vect_type, &
       & psb_z_base_vect_type, psb_zlinmap_type, psb_dpk_, &
       & psb_lz_csr_sparse_mat,  psb_lz_coo_sparse_mat, &
       & psb_z_csr_sparse_mat,  psb_z_coo_sparse_mat, &
       & psb_ipk_, psb_epk_, psb_lpk_, psb_desc_type, psb_i_base_vect_type, &
       & psb_erractionsave, psb_error_handler, psb_success_, psb_toupper
  !
  !
  !
  !> \class amg_z_base_aggregator_type
  !!
  !!  It is the data type containing the basic interface definition for
  !!  building a multigrid hierarchy by aggregation. The base object has no attributes,
  !!  it is intended to be essentially an abstract type.
  !!
  !!
  !!  type amg_z_base_aggregator_type
  !!  end type
  !!
  !!
  !!  Methods:
  !!
  !!  bld_tprol   -   Build a tentative prolongator
  !!
  !!  mat_bld     -   Build prolongator/restrictor and coarse matrix ac
  !!
  !!  mat_asb     -   Convert  prolongator/restrictor/coarse matrix
  !!                  and fix their descriptor(s)
  !!
  !!  update_next -   Transfer information to the next level; default is
  !!                  to do nothing, i.e. aggregators at different
  !!                  levels are independent.
  !!
  !!  default     -   Apply defaults
  !!  set_aggr_type - For aggregator that have internal options.
  !!  fmt         - Return a short string description
  !!  descr       - Print a more detailed description
  !!
  !!  cseti, csetr, csetc  - Set internal parameters, if any
  !
  type amg_z_base_aggregator_type
    ! Do we want to purge explicit zeros when aggregating?
    logical :: do_clean_zeros
  contains
    procedure, pass(ag) :: bld_tprol   => amg_z_base_aggregator_build_tprol
    procedure, pass(ag) :: mat_bld     => amg_z_base_aggregator_mat_bld
    procedure, pass(ag) :: mat_asb     => amg_z_base_aggregator_mat_asb
    procedure, pass(ag) :: bld_map     => amg_z_base_aggregator_bld_map
    procedure, pass(ag) :: update_next => amg_z_base_aggregator_update_next
    procedure, pass(ag) :: clone       => amg_z_base_aggregator_clone
    procedure, pass(ag) :: free        => amg_z_base_aggregator_free
    procedure, pass(ag) :: default     => amg_z_base_aggregator_default
    procedure, pass(ag) :: descr       => amg_z_base_aggregator_descr
    procedure, pass(ag) :: sizeof      => amg_z_base_aggregator_sizeof
    procedure, pass(ag) :: set_aggr_type => amg_z_base_aggregator_set_aggr_type
    procedure, nopass   :: fmt         => amg_z_base_aggregator_fmt
    procedure, pass(ag) :: cseti       => amg_z_base_aggregator_cseti
    procedure, pass(ag) :: csetr       => amg_z_base_aggregator_csetr
    procedure, pass(ag) :: csetc       => amg_z_base_aggregator_csetc
    generic, public     :: set         => cseti, csetr, csetc
    procedure, nopass   :: xt_desc     => amg_z_base_aggregator_xt_desc
  end type amg_z_base_aggregator_type

  abstract interface
    subroutine amg_z_soc_map_bld(iorder,theta,clean_zeros,a,desc_a,nlaggr,ilaggr,info)
      import :: psb_zspmat_type, psb_desc_type, psb_dpk_, psb_ipk_, psb_lpk_
      implicit none
      integer(psb_ipk_), intent(in)     :: iorder
      logical, intent(in)               :: clean_zeros
      type(psb_zspmat_type), intent(in) :: a
      type(psb_desc_type), intent(in)    :: desc_a
      real(psb_dpk_), intent(in)         :: theta
      integer(psb_lpk_), allocatable, intent(out)  :: ilaggr(:),nlaggr(:)
      integer(psb_ipk_), intent(out)               :: info
    end subroutine amg_z_soc_map_bld
  end interface

  interface amg_ptap_bld
    subroutine amg_z_ptap_bld(a_csr,desc_a,nlaggr,parms,ac,&
         & coo_prol,desc_cprol,coo_restr,info,desc_ax)
      import :: psb_z_csr_sparse_mat, psb_zspmat_type, psb_desc_type, &
           & psb_z_coo_sparse_mat, amg_dml_parms, psb_dpk_, psb_ipk_, psb_lpk_
      implicit none
      type(psb_z_csr_sparse_mat), intent(inout) :: a_csr
      type(psb_desc_type), intent(inout)         :: desc_a
      integer(psb_lpk_), intent(inout)           :: nlaggr(:)
      type(amg_dml_parms), intent(inout)         :: parms
      type(psb_z_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr
      type(psb_desc_type), intent(inout)         :: desc_cprol
      type(psb_zspmat_type), intent(out)        :: ac
      integer(psb_ipk_), intent(out)             :: info
      type(psb_desc_type), intent(inout), optional :: desc_ax
    end subroutine amg_z_ptap_bld
  end interface amg_ptap_bld

  interface amg_rap
    subroutine amg_z_rap(a_csr,desc_a,nlaggr,parms,ac,&
         & coo_prol,desc_cprol,coo_restr,info)
      import :: psb_z_csr_sparse_mat, psb_zspmat_type, psb_desc_type, &
           & psb_z_coo_sparse_mat, amg_dml_parms, psb_dpk_, psb_ipk_, psb_lpk_
      implicit none
      type(psb_z_csr_sparse_mat), intent(inout) :: a_csr
      type(psb_desc_type), intent(inout)         :: desc_a
      integer(psb_lpk_), intent(inout)           :: nlaggr(:)
      type(amg_dml_parms), intent(inout)         :: parms
      type(psb_z_coo_sparse_mat), intent(inout) :: coo_prol, coo_restr
      type(psb_desc_type), intent(inout)         :: desc_cprol
      type(psb_zspmat_type), intent(out)        :: ac
      integer(psb_ipk_), intent(out)             :: info
    end subroutine amg_z_rap
  end interface amg_rap

contains

  subroutine amg_z_base_aggregator_cseti(ag,what,val,info,idx)

    Implicit None

    ! Arguments
    class(amg_z_base_aggregator_type), intent(inout) :: ag
    character(len=*), intent(in)                  :: what
    integer(psb_ipk_), intent(in)                 :: val
    integer(psb_ipk_), intent(out)                :: info
    integer(psb_ipk_), intent(in), optional       :: idx
    ! Do nothing
    info = 0
  end subroutine amg_z_base_aggregator_cseti

  subroutine amg_z_base_aggregator_csetr(ag,what,val,info,idx)

    Implicit None

    ! Arguments
    class(amg_z_base_aggregator_type), intent(inout) :: ag
    character(len=*), intent(in)                  :: what
    real(psb_dpk_), intent(in)                 :: val
    integer(psb_ipk_), intent(out)                :: info
    integer(psb_ipk_), intent(in), optional       :: idx
    ! Do nothing
    info = 0
  end subroutine amg_z_base_aggregator_csetr

  subroutine amg_z_base_aggregator_csetc(ag,what,val,info,idx)

    Implicit None

    ! Arguments
    class(amg_z_base_aggregator_type), intent(inout) :: ag
    character(len=*), intent(in)                  :: what
    character(len=*), intent(in)                 :: val
    integer(psb_ipk_), intent(out)                :: info
    integer(psb_ipk_), intent(in), optional       :: idx
    ! Set clean zeros, or do nothing.
    select case (psb_toupper(trim(what)))
    case('AGGR_CLEAN_ZEROS')
      select case (psb_toupper(trim(val)))
      case('TRUE','T')
        ag%do_clean_zeros = .true.
      case('FALSE','F')
        ag%do_clean_zeros = .false.
      end select
    end select
    info = 0
  end subroutine amg_z_base_aggregator_csetc

  subroutine  amg_z_base_aggregator_update_next(ag,agnext,info)
    implicit none
    class(amg_z_base_aggregator_type), target, intent(inout) :: ag, agnext
    integer(psb_ipk_), intent(out)       :: info

    !
    ! Base version does nothing.
    !
    info = 0
  end subroutine amg_z_base_aggregator_update_next

  subroutine  amg_z_base_aggregator_clone(ag,agnext,info)
    implicit none
    class(amg_z_base_aggregator_type), intent(inout) :: ag
    class(amg_z_base_aggregator_type), allocatable, intent(inout) :: agnext
    integer(psb_ipk_), intent(out)       :: info

    info = 0
    if (allocated(agnext)) then
      call agnext%free(info)
      if (info == 0) deallocate(agnext,stat=info)
    end if
    if (info /= 0) return
    allocate(agnext,source=ag,stat=info)

  end subroutine amg_z_base_aggregator_clone

  subroutine  amg_z_base_aggregator_free(ag,info)
    implicit none
    class(amg_z_base_aggregator_type), intent(inout) :: ag
    integer(psb_ipk_), intent(out)       :: info

    info = psb_success_
    return
  end subroutine amg_z_base_aggregator_free

  subroutine  amg_z_base_aggregator_default(ag)
    implicit none
    class(amg_z_base_aggregator_type), intent(inout) :: ag
    ! Only one default setting
    ag%do_clean_zeros = .true.

    return
  end subroutine amg_z_base_aggregator_default

  function amg_z_base_aggregator_fmt() result(val)
    implicit none
    character(len=32)  :: val

    val = "Default aggregator "
  end function amg_z_base_aggregator_fmt

  function amg_z_base_aggregator_sizeof(ag) result(val)
    implicit none
    class(amg_z_base_aggregator_type), intent(in)  :: ag
    integer(psb_epk_)  :: val

    val = 1
  end function amg_z_base_aggregator_sizeof

  function amg_z_base_aggregator_xt_desc() result(val)
    implicit none
    logical  :: val

    val = .false.
  end function amg_z_base_aggregator_xt_desc

  subroutine  amg_z_base_aggregator_descr(ag,parms,iout,info)
    implicit none
    class(amg_z_base_aggregator_type), intent(in) :: ag
    type(amg_dml_parms), intent(in)   :: parms
    integer(psb_ipk_), intent(in)  :: iout
    integer(psb_ipk_), intent(out) :: info

    write(iout,*) 'Aggregator object type: ',ag%fmt()
    call parms%mldescr(iout,info)

    return
  end subroutine amg_z_base_aggregator_descr

  subroutine  amg_z_base_aggregator_set_aggr_type(ag,parms,info)
    implicit none
    class(amg_z_base_aggregator_type), intent(inout) :: ag
    type(amg_dml_parms), intent(in)   :: parms
    integer(psb_ipk_), intent(out) :: info

    ! Do nothing

    return
  end subroutine amg_z_base_aggregator_set_aggr_type

  !
  !> Function  bld_tprol:
  !! \memberof  amg_z_base_aggregator_type
  !! \brief  Build a tentative prolongator.
  !!         The routine will map the local matrix entries to aggregates.
  !!         The mapping is store in ILAGGR; for each local row index I,
  !!         ILAGGR(I) contains the index of the aggregate to which index I
  !!         will contribute, in global numbering.
  !!         Many aggregations produce a binary tentative prolongator, but some
  !!         do not, hence we also need the OP_PROL output.
  !!         AG_DATA is passed here just in case some of the
  !!         aggregators need it internally, most of them will ignore.
  !!
  !!  \param ag      The input aggregator object
  !!  \param parms   The auxiliary parameters object
  !!  \param ag_data Auxiliary global aggregation info
  !!  \param a       The local matrix part
  !!  \param desc_a  The descriptor
  !!  \param ilaggr  Output aggregation map
  !!  \param nlaggr  Sizes of ilaggr on all processes
  !!  \param op_prol The tentative prolongator operator
  !!  \param info    Return code
  !!
  !
  subroutine  amg_z_base_aggregator_build_tprol(ag,parms,ag_data,&
       & a,desc_a,ilaggr,nlaggr,t_prol,info)
    use psb_base_mod
    implicit none
    class(amg_z_base_aggregator_type), target, intent(inout) :: ag
    type(amg_dml_parms), intent(inout)  :: parms
    type(amg_daggr_data), intent(in)    :: ag_data
    type(psb_zspmat_type), intent(inout) :: a
    type(psb_desc_type), intent(inout)     :: desc_a
    integer(psb_lpk_), allocatable, intent(out) :: ilaggr(:),nlaggr(:)
    type(psb_lzspmat_type), intent(out)  :: t_prol
    integer(psb_ipk_), intent(out)      :: info

    integer(psb_ipk_) :: err_act
    character(len=20) :: name='z_base_aggregator_build_tprol'

    call psb_erractionsave(err_act)

    info = psb_err_missing_override_method_
    call psb_errpush(info,name)
    goto 9999

    call psb_erractionrestore(err_act)
    return

9999 call psb_error_handler(err_act)

    return

  end subroutine amg_z_base_aggregator_build_tprol

  !
  !> Function   mat_bld
  !! \memberof  amg_z_base_aggregator_type
  !! \brief     Build prolongator/restrictor/coarse matrix.
  !!
  !!
  !!  \param ag        The input aggregator object
  !!  \param parms     The auxiliary parameters object
  !!  \param a         The local matrix part
  !!  \param desc_a    The descriptor
  !!  \param ilaggr    Aggregation map
  !!  \param nlaggr    Sizes of ilaggr on all processes
  !!  \param ac        On output the coarse matrix
  !!  \param op_prol   On input, the  tentative prolongator operator, on output
  !!                   the final prolongator
  !!  \param op_restr  On output, the restrictor operator;
  !!                   in many cases it is the transpose of the prolongator.
  !!  \param info    Return code
  !!
  subroutine  amg_z_base_aggregator_mat_bld(ag,parms,a,desc_a,ilaggr,nlaggr,&
       & ac,desc_ac,op_prol,op_restr,t_prol,info)
    use psb_base_mod
    implicit none
    class(amg_z_base_aggregator_type), target, intent(inout) :: ag
    type(amg_dml_parms), intent(inout)   :: parms
    type(psb_zspmat_type), intent(in)    :: a
    type(psb_desc_type), intent(inout)   :: desc_a
    integer(psb_lpk_), intent(inout)     :: ilaggr(:), nlaggr(:)
    type(psb_lzspmat_type), intent(inout) :: t_prol
    type(psb_zspmat_type), intent(out)   :: op_prol, ac,op_restr
    type(psb_desc_type), intent(inout)   :: desc_ac
    integer(psb_ipk_), intent(out)       :: info
    integer(psb_ipk_) :: err_act
    character(len=20) :: name='z_base_aggregator_mat_bld'

    call psb_erractionsave(err_act)

    info = psb_err_missing_override_method_
    call psb_errpush(info,name)
    goto 9999

    call psb_erractionrestore(err_act)
    return

9999 call psb_error_handler(err_act)

    return
  end subroutine amg_z_base_aggregator_mat_bld

  !
  !> Function   mat_asb
  !! \memberof  amg_z_base_aggregator_type
  !! \brief     Build prolongator/restrictor/coarse matrix.
  !!
  !!
  !!  \param ag        The input aggregator object
  !!  \param parms     The auxiliary parameters object
  !!  \param a         The local matrix part
  !!  \param desc_a    The descriptor
  !!  \param ilaggr    Aggregation map
  !!  \param nlaggr    Sizes of ilaggr on all processes
  !!  \param ac        On output the coarse matrix
  !!  \param op_prol   On input, the  tentative prolongator operator, on output
  !!                   the final prolongator
  !!  \param op_restr  On output, the restrictor operator;
  !!                   in many cases it is the transpose of the prolongator.
  !!  \param info    Return code
  !!
  subroutine  amg_z_base_aggregator_mat_asb(ag,parms,a,desc_a,&
       & ac,desc_ac, op_prol,op_restr,info)
    use psb_base_mod
    implicit none
    class(amg_z_base_aggregator_type), target, intent(inout) :: ag
    type(amg_dml_parms), intent(inout)   :: parms
    type(psb_zspmat_type), intent(in)    :: a
    type(psb_desc_type), intent(inout)     :: desc_a
    type(psb_zspmat_type), intent(inout) :: op_prol,ac,op_restr
    type(psb_desc_type), intent(inout)      :: desc_ac
    integer(psb_ipk_), intent(out)       :: info
    integer(psb_ipk_) :: err_act
    character(len=20) :: name='z_base_aggregator_mat_asb'

    call psb_erractionsave(err_act)

    info = psb_err_missing_override_method_
    call psb_errpush(info,name)
    goto 9999

    call psb_erractionrestore(err_act)
    return

9999 call psb_error_handler(err_act)

    return
  end subroutine amg_z_base_aggregator_mat_asb

  !
  !> Function   bld_map
  !! \memberof  amg_z_base_aggregator_type
  !! \brief     Build linear map between hierarchy levels
  !!
  !!
  !!  \param ag        The input aggregator object
  !!  \param desc_a    The fine space descriptor
  !!  \param desc_ac   The coarse space descriptor
  !!  \param ilaggr    Aggregation map vector
  !!  \param nlaggr    Sizes of ilaggr on all processes
  !!  \param op_prol   The  prolongator operator
  !!  \param op_restr  The restrictor operator
  !!  \param map       The output map
  !!  \param info    Return code
  !!
  subroutine  amg_z_base_aggregator_bld_map(ag,desc_a,desc_ac,ilaggr,nlaggr,&
       & op_restr,op_prol,map,info)
    use psb_base_mod
    implicit none
    class(amg_z_base_aggregator_type), target, intent(inout) :: ag
    type(psb_desc_type), intent(in), target :: desc_a, desc_ac
    integer(psb_lpk_), intent(inout)     :: ilaggr(:), nlaggr(:)
    type(psb_zspmat_type), intent(inout)   :: op_restr, op_prol
    type(psb_zlinmap_type), intent(out)    :: map
    integer(psb_ipk_), intent(out)       :: info
    integer(psb_ipk_) :: err_act
    character(len=20) :: name='z_base_aggregator_bld_map'

    call psb_erractionsave(err_act)
    !
    ! Copy the prolongation/restriction matrices into the descriptor map.
    !  op_restr => PR^T   i.e. restriction  operator
    !  op_prol => PR     i.e. prolongation operator
    !
    !  WARNING: need to check whether the  copy into IOP_RESTR/IOP_PROL
    !  is safe or not.
    !
    !  This default implementation reuses desc_a/desc_ac through
    !  pointers in the map structure.
    !
    map = psb_linmap(psb_map_aggr_,desc_a,&
         & desc_ac,op_restr,op_prol,ilaggr,nlaggr)

    call psb_erractionrestore(err_act)
    return

9999 call psb_error_handler(err_act)

    return
  end subroutine amg_z_base_aggregator_bld_map

end module amg_z_base_aggregator_mod


text:
PROGRAM ScalarWaveTest

  USE KindModule, ONLY: &
    DP, One, Two, Zero, TwoPi
  USE ProgramHeaderModule, ONLY: &
    iX_B0, iX_B1, iX_E0, iX_E1, &
    nDimsX
  USE ReferenceElementModuleX, ONLY: &
    InitializeReferenceElementX, &
    FinalizeReferenceElementx
  USE ReferenceElementModuleX_Lagrange, ONLY: &
    InitializeReferenceElementX_Lagrange, &
    FinalizeReferenceElementX_Lagrange
  USE ScalarWave_ProgramInitializationModule, ONLY: &
    InitializeProgram, &
    FinalizeProgram
  USE InitializationModule, ONLY: &
    InitializeFields
  USE ScalarWave_InputOutputModuleHDF, ONLY: &
    WriteFieldsHDF, &
    ReadFieldsHDF
  USE TimeSteppingModule_SSPRK, ONLY: &
    InitializeScalarWave_SSPRK, &
    FinalizeScalarWave_SSPRK, &
    UpdateScalarWave_SSPRK
  USE ScalarFieldsModule, ONLY : &
    uSF, &
    ComputeTimeStep_ScalarWave
  USE ScalarWave_dgDiscretizationModule, Only : &
    ComputeIncrement_ScalarWave_DG_Explicit

  IMPLICIT NONE

  LOGICAL  :: wrt
  INTEGER  :: iCycle, nNodes, nWrites
  REAL(DP) :: t, dt, t_end, dt_wrt, t_wrt

  t = 0.0_DP
  t_end = TwoPi
  dt_wrt = 1.0d-2 * t_end

  nNodes = 2

  CALL InitializeProgram &
         ( ProgramName_Option = 'SpacetimeTest', &
           nX_Option = [ 2048, 1, 1 ], swX_Option = [ 1, 0, 0 ], &
           bcX_Option = [ 1, 0, 0 ], xL_Option = [ Zero, Zero, Zero ], &
           xR_Option = [ TwoPi, One, One ], nNodes_Option = 3 )

  CALL InitializeReferenceElementX

  CALL InitializeReferenceElementX_Lagrange

  Call InitializeFields

  CALL WriteFieldsHDF( 0.0_DP, WriteSF_Option=.true.)

  nWrites = 1

  CALL InitializeScalarWave_SSPRK( nStages = 3 )

  t_wrt = t + dt_wrt
  wrt   = .FALSE.

  iCycle = 0
  DO WHILE ( t < t_end )

    iCycle = iCycle + 1

    CALL ComputeTimeStep_ScalarWave &
           ( iX_B0, iX_E0, iX_B1, iX_E1, uSF, &
             CFL = 0.5_DP / (nDimsX * ( Two * DBLE ( nNodes ) - One ) ), &
             Timestep = dt )

    IF( t + dt > t_end )THEN

      dt = t_end - t

    END IF

    IF( t + dt > t_wrt )THEN

      dt    = t_wrt - t
      t_wrt = t_wrt + dt_wrt
      wrt   = .TRUE.

    END IF

    CALL UpdateScalarWave_SSPRK &
           ( t, dt, uSF, ComputeIncrement_ScalarWave_DG_Explicit )

    t = t + dt

    IF( wrt )THEN

      CALL WriteFieldsHDF &
             ( t , WriteSF_Option = .TRUE. )

      wrt = .FALSE.

      nWrites = nWrites + 1
    END IF

  END DO

  CALL WriteFieldsHDF &
         ( t, WriteSF_Option = .TRUE. )

  WRITE(*,*)
  WRITE(*,'(A6,A,I6.6,A)') &
    '', 'Finished ', iCycle, ' Cycles.'
  WRITE(*,*)
  WRITE(*,'(A6,A,I6.6,A)') &
    '', 'Wrote ', nWrites, ' Files.'

  CALL FinalizeScalarWave_SSPRK

  CALL FinalizeReferenceElementX_Lagrange

  CALL FinalizeReferenceElementX

  CALL FinalizeProgram

END PROGRAM ScalarWaveTest

max_stars_repo_name: srichers/thornado

max_issues_repo_name: srichers/thornado

max_forks_repo_name: srichers/thornado


text:
!-----------------------------------------------------------------------
      program test
      implicit none
      include 'ceed/fortran.h'

      integer ceed,err
      integer mult
      integer r
      integer i
      integer*8 moffset

      integer ne
      parameter(ne=3)

      integer*4 ind(4*ne)
      real*8 mm(3*ne+1)
      integer offset
      real*8 diff

      character arg*32

      call getarg(1,arg)
      call ceedinit(trim(arg)//char(0),ceed,err)

      call ceedvectorcreate(ceed,3*ne+1,mult,err)
      call ceedvectorsetvalue(mult,0.d0,err);

      do i=1,ne
        ind(4*i-3)=3*i-3
        ind(4*i-2)=3*i-2
        ind(4*i-1)=3*i-1
        ind(4*i-0)=3*i-0
      enddo
      call ceedelemrestrictioncreate(ceed,ne,4,1,1,3*ne+1,ceed_mem_host,&
     & ceed_use_pointer,ind,r,err)

      call ceedelemrestrictiongetmultiplicity(r,mult,err)

      call ceedvectorgetarrayread(mult,ceed_mem_host,mm,moffset,err)
      do i=1,3*ne+1
        if(i > 1 .and. i < 3*ne+1 .and. mod(i-1,3)==0) then
          offset = 1
        else
          offset = 0
        endif
        diff=1+offset-mm(i+moffset)
        if (abs(diff) > 1.0D-15) then
! LCOV_EXCL_START
          write(*,*) 'Error in multiplicity vector: mult(',i,')=',mm(i+moffset)
! LCOV_EXCL_STOP
        endif
      enddo
      call ceedvectorrestorearrayread(mult,mm,moffset,err)

      call ceedvectordestroy(mult,err)
      call ceedelemrestrictiondestroy(r,err)
      call ceeddestroy(ceed,err)

      end
!-----------------------------------------------------------------------

max_stars_repo_name: AdelekeBankole/libCEED

max_issues_repo_name: AdelekeBankole/libCEED

max_forks_repo_name: AdelekeBankole/libCEED


text:
!     QSATS version 1.0 (3 March 2011)

!     file name: eloc.f

! ----------------------------------------------------------------------

!     this computes the total energy and the expectation value of the
!     potential energy from the snapshots recorded by QSATS.

! ----------------------------------------------------------------------

      module domain

      integer (kind = 4), parameter :: id = 2 ! classical DoF defining the domain function
      integer (kind = 4), parameter :: nd = 8 ! number of domains

!      integer (kind = 4), parameter ::

      real (kind = 8), dimension(nd+1) ::  ydom
      real (kind = 8) ymin,ymax,dy

      save
      end module

      module cdat

      integer (kind = 4), parameter :: nq = 1 ! number of QDoF
      real (kind = 8),parameter ::  half = 0.5d0, one = 1.0d0
      real (kind = 8),parameter ::  pi = 4.0*atan(1d0)

      contains

      double precision function step(x,y,z)
      implicit double precision (a-h,o-z)
      if(z > x .and. z < y) then
         step = 1d0
      else
        step = 0d0
      endif
      return
      end function

      end module
! ----------------------------
! --- main program
! ------------------------------
      program eloc
      use cdat
      use domain

      implicit double precision (a-h, o-z)

!------------trajectory initial allocate------------------------
      real*8, dimension(:,:),allocatable :: x,p,du,ap,ar,rp,fr

      real*8, dimension(:,:),allocatable :: dv
      real*8, dimension(:),allocatable :: am,p0,w,alpha,x0,pe,ke,u,cf
      real*8 :: gasdev

      integer (kind = 4), dimension(:), allocatable :: idum

      real (kind = 8), dimension(:), allocatable :: enk,env
      real (kind = 8), dimension(:), allocatable :: q

!------------------------------------------------------
!      character*4 :: atom(NATOMS)
!      real*8 :: dvl(NATOM3)
!      integer*4  :: ntraj,ndim,kmax

c --- this common block is used to enable interpolation in the potential
c     energy lookup table in the subroutine local below.

      open(100, file = 'en.dat'  )
      open(101, file = 'xoutput' )
      open(102, file = 'traj.dat')
      open(103, file = 'den.dat' ) ! density overlap

!------read parameters from IN----------------------
      open(10,file='IN')

      read(10,*) Ntraj
      read(10,*) ndim
      read(10,*) kmax,dt
      read(10,*) a0

      close(10)

!      Ndim = NATOM3

c --- initialization.

!      call tstamp

      allocate(dv(ndim,ntraj))
!-----------------------------------------------------
      allocate(env(ntraj),enk(Ntraj),x(Ndim,Ntraj),
     +         p(Ndim,Ntraj),du(Ndim,Ntraj),w(Ntraj))

      allocate(p0(Ndim),alpha(Ndim),am(Ndim),x0(Ndim),cf(Ndim))

      allocate( idum(ndim) )
      allocate( q(ndim) )

!      binvrs=one/bin

c ------ compute the local energy and the potential energy.

      dt2 = dt/2d0
      t   = 0d0
      pow = 6d0
      am0 = 1d0

      am(1) = 1d0
      am(2) = 10d0

      alpha(1) = 1.118
      alpha(2) = 11.18

      x0(1) = 0d0
      x0(2) = -3d0

      p0(1) = 0d0
      p0(2) = 22d0

! --- domain setup

      ymin = -4.0d0
      ymax = 6.0d0
      dy = (ymax-ymin)/nd
      do i=2,nd
        ydom(i) = ymin+(i-1)*dy
      enddo
      ydom(1) = -100d0
      ydom(nd+1) = 100d0

      call seed(idum,Ndim)

      write(*,1010) Ntraj,Ndim,kmax,dt,am(1),am(2)
1010  format('Initial Conditions'//,
     +       'Ntraj   = ' , i6/,
     +       'Ndim    = ' , i6/,
     +       'kmax    = ' , i6/,
     +       'dt      = ' , f7.4/,
     +       'Mass    = ' , f7.4,f7.4/)

      write(*,1011) alpha(1),alpha(2),x0(1),x0(2),p0(1),p0(2)
1011  format('Initial Wavepacket'//,
     +       'GWP width   = ' , 2(f7.4,1x)/,
     +       'GWP center  = ' , 2(f7.4,1x)/,
     +       'Momentum    = ' , 2(f7.4,1x)/)

! --- initial Lagrangian grid points (evolve with time)
      do i=1,Ntraj
        do j=1,Ndim
1100      x(j,i)=gasdev(idum(j))
          x(j,i)=x(j,i)/dsqrt(4d0*alpha(j))+x0(j)
!          if((x(j,i)-x0(j))**2 .gt. pow/2d0/alpha(j)) goto 1100
        enddo
      enddo

! --- initial momentum for QTs and weights

      do i=1,Ntraj
        p(1,i) = p0(1)
      enddo
      do i=1,ntraj
            p(2,i) = p0(2)+alpha(2)*(x(2,i)-x0(2))
!            rp(j,i) = -2d0*alpha(j)*(x(j,i)-x0(j))
      enddo

      do i=1,ntraj
        w(i) = 1d0/dble(ntraj)
      enddo

!---- expectation value of x(1,ntraj)--------
      av = 0d0
      do i=1,ntraj
        av = av+w(i)*x(1,i)
      enddo
      write(*,6008) av
6008  format('Expectation value of x(1,ntraj) = ', f10.6)

! --- half timestep for p

      call fit(am,x,p,w,ndim,ntraj,enq,du)
      call local(ndim,ntraj,x,env,dv)
      do i=1,ntraj
          do j=1,ndim
            p(j,i) = p(j,i)-(dv(j,i)+du(j,i))*dt2
          enddo
      enddo

! --- trajectories propagation
      do 10 kt=1,kmax

        t=t+dt

! ----- half-step increments of momenta & full step increment of positions
        do i=1,Ntraj

!          call derivs(q,vloc,dv)
!          call long_force(q,vlong,dvl)
!          write(*,*) 'test2'
!          env(i) = vloc

          do j=1,Ndim
!              p(j,i)=p(j,i)+(-dv(j)-du(j,i))*dt2
              x(j,i)=x(j,i)+p(j,i)*dt/am(j)
          enddo
        enddo

!          write(*,*) 'test3',x(1,2),x(2,3)

! ----- half-step increments of momenta
        do i=1,Ntraj

          call local(ndim,ntraj,x,env,dv)
!          call long_force(q,vlong,dvl)
          call fit(am,x,p,w,ndim,ntraj,enq,du)

          do j=1,Ndim
            p(j,i)=p(j,i)+(-dv(j,i)-du(j,i))*dt
          enddo

        enddo

!-------update potential, kinetic, and total energy
          do i=1, Ntraj
            enk(i) = 0d0
            do j=1,Ndim
              enk(i) = p(j,i)**2/(2d0*am(j)) + enk(i)
            enddo
          enddo

          write(101,1000) t,(x(1,i),i=1,20)
          write(102,1000) t,(x(j,4),j=1,20)

        call aver(Ntraj,w,enk,enk_ave)
        call aver(Ntraj,w,env,env_ave)
!        call aver(Ntraj,w,u,qu)

        ent = env_ave+enk_ave+enq

!        do i=1,ntraj
!          write(105,1000) dble(i),ev(i)
!        enddo

        write(100,1000) t,enk_ave,env_ave,enq,ent
        call flush(100)

! --- density overlap
      d0 = 0d0
      d1 = 0d0
      ak1 = 5d0
      ak2 = 15d0

      alfa = dsqrt(am(1)*ak1)/2d0
      beta = dsqrt(am(1)*ak2)/2d0

      do i=1,ntraj
        if(x(2,i) > 0d0) then
        d0 = d0 + w(i)*exp(-alfa*x(1,i)**2)
        d1 = d1 + w(i)*exp(-beta*x(1,i)**2)
       endif
      enddo
      d0 = d0*dsqrt(alfa/pi)
      d1 = d1*dsqrt(beta/pi)

      write(103,1000) t,d0,d1

! --- store temporary data
!      open(106,file='temp.dat')
!      if (mod(kt,1000) .eq. 0) then
!        do i=1,Ntraj
!          do j=1,Ndim
!            write(106,1000) x(j,i),p(j,i)
!            call flush(106)
!          enddo
!        enddo
!      endif

!      close(106)

10    enddo

      write(*,1020) ent
1020  format('Total Energy =', f10.5/ ,
     +       'MISSION COMPLETE.')

!          call local(q,vloc,dv)

!       convert to eV
!          vloc = vloc*27.211d0
c ------ convert to kelvin per atom.

!         tloc=tloc/(3.1668513d-6*dble(NATOMS))
!          vloc=vloc/(3.1668513d-6*dble(NATOMS))

c ------ accumulate the results.

!         vtavg(k)=vtavg(k)+vloc
!         vtavg2(k)=vtavg2(k)+(vloc)**2
!         etavg(k)=etavg(k)+tloc+vloc
!         etavg2(k)=etavg2(k)+(tloc+vloc)**2

!350      continue

!      end do

!      goto 300

c --- account for overshooting.

!600   loop=loop-1

!      write (6, 6600) loop
!6600  format ('number of snapshots = ', i6/)

c --- compute the averages and standard deviations.

!      do k=1, NREPS, 11
!
!      vtavg(k)=vtavg(k)/dble(loop)
!      vtavg2(k)=vtavg2(k)/dble(loop)
!      etavg(k)=etavg(k)/dble(loop)
!      etavg2(k)=etavg2(k)/dble(loop)
!
!      vsd=sqrt(vtavg2(k)-vtavg(k)**2)
!      esd=sqrt(etavg2(k)-etavg(k)**2)
!
!      write (6, 6610) k, 'VAVG = ', vtavg(k)
!6610  format ('replica ', i3, 1x, a7, f10.5, ' Kelvin')
!          write (100,1000) q(1),vloc
!       end do
!1001    format('potential at equilibrim configuration [K/atom] =',
!     &         f12.7)
!
!      write (6, 6610) k, 'V SD = ', vsd
!
!      write (6, 6610) k, 'EAVG = ', etavg(k)
!
!      write (6, 6610) k, 'E SD = ', esd
!
!      end do
1000  format(20(e14.7,1x))
      end program

c ----------------------------------------------------------------------

c     quit is a subroutine used to terminate execution if there is
c     an error.

c     it is needed here because the subroutine that reads the parameters
c     (subroutine input) may call it.

c ----------------------------------------------------------------------

      subroutine quit

      write (6, *) 'termination via subroutine quit'

      stop
      end subroutine
! ---------------------------------------------------------------
!     random number seeds
! ---------------------------------------------------------------
      subroutine seed(idum,Ndim)
      implicit real*8(a-h,o-z)
      integer*4, intent(IN) :: Ndim
      integer*4, intent(OUT) :: idum(Ndim)

      do i=1,Ndim
        idum(i) = 5 + i
      enddo

      return
      end subroutine
! ---------------------------------------------------------------
!     average over trajectories
! ---------------------------------------------------------------
      subroutine aver(Ntraj,w,x,y)
      implicit real*8(a-h,o-z)
      real*8 :: x(Ntraj), w(Ntraj)

      y = 0d0

      do i=1,Ntraj
        y = y+x(i)*w(i)
      enddo

      return
      end subroutine
! ----------------------------------------------------------------


text:
module batched_tensor_permute_module
    use :: data_api, only : stream

    use :: tensor_api, only : tensor_array_element

    implicit none
    private

    public :: batched_tensor_permute

    type, abstract :: batched_tensor_permute
    contains
        procedure(btp_interface), deferred :: btp
        procedure(sync_interface), deferred :: synchronize
        procedure(empty), deferred :: cleanup
    end type batched_tensor_permute

    abstract interface
        subroutine btp_interface(this, src, dst, perm, astream)
            import :: batched_tensor_permute
            import :: tensor_array_element
            import :: stream

            class(batched_tensor_permute), target, intent(inout) :: this
            class(tensor_array_element), dimension(:), intent(in) :: src
            class(tensor_array_element), dimension(:), intent(inout) :: dst
            integer, intent(in) :: perm(:)
            type(stream), intent(in), optional :: astream
        end subroutine btp_interface

        subroutine sync_interface(this, astream)
            import :: batched_tensor_permute
            import :: stream

            class(batched_tensor_permute), intent(in) :: this
            type(stream), intent(in), optional :: astream
        end subroutine sync_interface

        subroutine empty(this)
            import :: batched_tensor_permute

            class(batched_tensor_permute), intent(inout) :: this
        end subroutine empty
    end interface
end module batched_tensor_permute_module

max_stars_repo_name: cheshyre/ntcl-algorithms

max_issues_repo_name: cheshyre/ntcl-algorithms

max_forks_repo_name: cheshyre/ntcl-algorithms


text:
!
!     file thstssp.f90
!
!     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!     *                                                               *
!     *                  copyright (c) 2005 by UCAR                   *
!     *                                                               *
!     *       University Corporation for Atmospheric Research         *
!     *                                                               *
!     *                      all rights reserved                      *
!     *                                                               *
!     *                         Fishpack                              *
!     *                                                               *
!     *                      A Package of Fortran                     *
!     *                                                               *
!     *                Subroutines and Example Programs               *
!     *                                                               *
!     *               for Modeling Geophysical Processes              *
!     *                                                               *
!     *                             by                                *
!     *                                                               *
!     *        John Adams, Paul Swarztrauber and Roland Sweet         *
!     *                                                               *
!     *                             of                                *
!     *                                                               *
!     *         the National Center for Atmospheric Research          *
!     *                                                               *
!     *                Boulder, Colorado  (80307)  U.S.A.             *
!     *                                                               *
!     *                   which is sponsored by                       *
!     *                                                               *
!     *              the National Science Foundation                  *
!     *                                                               *
!     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!
program test_hstssp

    use, intrinsic :: ISO_Fortran_env, only: &
        stdout => OUTPUT_UNIT

    use fishpack

    ! Explicit typing only
    implicit none

    !-----------------------------------------------
    ! Local variables
    !-----------------------------------------------
    integer(ip), parameter :: M = 18
    integer(ip), parameter :: N = 72
    integer(ip), parameter :: IDIMF = M
    integer(ip)            :: mbdcnd, nbdcnd, i, j, ierror
    real(wp)               :: f(IDIMF,N), sint(M)
    real(wp), dimension(1) :: bda, bdc, bdd
    real(wp), dimension(N) :: sinp, bdb
    real(wp)               :: a, b, c, d, elmbda, dtheta, dphi
    real(wp)               :: pertrb
    real(wp), parameter    :: ZERO = 0.0_wp, HALF = 0.5_wp
    real(wp), parameter    :: TWO = 2.0_wp
    !-----------------------------------------------

    ! Set domain
    a = ZERO
    b = HALF_PI
    c = ZERO
    d = TWO_PI

    ! Set boundary conditions
    nbdcnd = 0
    mbdcnd = 6

    ! Set helmholtz constant
    elmbda = ZERO

    ! Set mesh sizes
    dtheta = (b - a)/M
    dphi = (d - c)/N

    ! Generate sines for use in subsequent computations
    do i = 1, M
        sint(i) = sin((real(i, kind=wp) - HALF) * dtheta)
    end do

    do j = 1, N
        sinp(j) = sin((real(j, kind=wp) - HALF) * dphi)
    end do

    ! Compute right hand side of equation and store in f
    block
        real(wp), parameter :: SIX = 6.0_wp
        do j = 1, N
            f(:, j) = TWO - SIX * (sint * sinp(j))**2
        end do
    end block
    !
    ! Set boundary data; store derivative data at the equator
    ! In our example, bda, bdc, and bdd are 1-dimensional dummy variables.
    bdb = ZERO

    ! Solve 2D Helmholtz in spherical coordinates on staggered grid
    call hstssp(a, b, M, mbdcnd, bda, bdb, c, d, N, nbdcnd, bdc, bdd, &
        elmbda, f, IDIMF, pertrb, ierror)

    ! Compute discretization error. since problem is singular, the
    ! solution must be normalized.
    block
        real(wp), parameter :: KNOWN_PERTRB = 0.635830001454109e-3_wp
        real(wp), parameter :: KNOWN_ERROR = 0.337523232257420e-2_wp
        real(wp)            :: discretization_error
        real(wp)            :: exact_solution(M,N)

        do j = 1, N
            do i = 1, M
                exact_solution(i,j) = (sint(i) * sinp(j))**2 + f(1, 1)
            end do
        end do

        ! Set discretization error
        discretization_error = maxval(abs(exact_solution - f(:M,:N)))

        call check_output('hstssp', &
            ierror, KNOWN_ERROR, discretization_error, KNOWN_PERTRB, pertrb)
    end block

end program test_hstssp

max_stars_repo_name: jlokimlin/fishpack

max_issues_repo_name: jlokimlin/fishpack

max_forks_repo_name: jlokimlin/fishpack


text:
!!  Copyright (C)  Stichting Deltares, 2012-2016.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License version 3,
!!  as published by the Free Software Foundation.
!!
!!  This program is distributed in the hope that it will be useful,
!!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!!  GNU General Public License for more details.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program. If not, see <http://www.gnu.org/licenses/>.
!!
!!  contact: delft3d.support@deltares.nl
!!  Stichting Deltares
!!  P.O. Box 177
!!  2600 MH Delft, The Netherlands
!!
!!  All indications and logos of, and references to registered trademarks
!!  of Stichting Deltares remain the property of Stichting Deltares. All
!!  rights reserved.

!    Date:       30 Oct 1989
!    Time:       14:47
!    Program:    BVECT    FORTRAN
!    Version:    1.1
!    Programmer: Hans Los
!    (c) 1989 Deltares Sektor W&M
!    Previous versions:
!    1.0 -- 30 Oct 1989 -- 14:45
!    BVECT.FOR  -- 25 Apr 1989 -- 10:28
!
!  *********************************************************************
!  *  SUBROUTINE TO SET THE MORTALITY CONSTRAINTS INTO THE B-VECTOR    *
!  *********************************************************************
!
      SUBROUTINE BVECT(X,XDEF)
      IMPLICIT REAL*8 (A-H,O-Z)
!
      INCLUDE 'blmdim.inc'
      INCLUDE 'phyt2.inc'
      INCLUDE 'matri.inc'
      DIMENSION X(*),B2(MS),XDEF(*)
!
! To tell Bloom how much biomass of the living phytoplankton
! species is left at the end of the time-step, the 'minimum
! biomass' of each species is set in the B-vector. This value is used
! as the mortality constraint (the minimum biomass to be returned
! by simplex).
! The new, minimum biomass levels are also stored in the original
! XDEF-vector. This is to enable the program to deal with infeasible
! solutions for example due to light limitation.
!

      I1 = 0
      K1 = NUROWS
      DO 40 I=1,NUECOG
          SUMSP = 0.0
          L1 = IT2(I,1)
          L2 = IT2(I,2)
          DO 30 K=L1,L2
              I1 = I1 + 1
              K1 = K1 + 1
              SUMSP = SUMSP + X(I1)
              XDEF(K1) = X(I1)
   30     CONTINUE
          B2(I) = DMAX1(SUMSP,0.0D0)
   40 CONTINUE
!
      I1 = NUEXRO + NUECOG
      DO 53 I = 1,NUECOG
          I1 = I1 + 1
          B(I1) = B2(I)
53    CONTINUE
!
      RETURN
      END

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
program t
  implicit none
  integer::i,j
  j=-1
  i=f(j)
  print '(i0)',i
  print '(i0)',j
contains
  integer function f(i)
    integer::i
    intent(inout)::i
    print '(i0)',i
    if (i.lt.0) then
      i=999
      f=i
      return
    endif
    f=i+1
  end function f
end program t

max_stars_repo_name: maddenp/ppp

max_issues_repo_name: maddenp/ppp

max_forks_repo_name: maddenp/ppp


text:
module test_stats_random
    use stdlib_kinds, only: int8, int16, int32, int64
    use stdlib_random, only : random_seed, dist_rand
    use testdrive, only: new_unittest, unittest_type, error_type, check

    implicit none

    private
    public :: collect_stats_random

contains

    !> Collect all exported unit tests
    subroutine collect_stats_random(testsuite)
        !> Collection of tests
        type(unittest_type), allocatable, intent(out) :: testsuite(:)

        testsuite = [ &
            new_unittest("random_seed", test_random_seed), &
            new_unittest("random_rand_iint8", test_random_rand_iint8), &
            new_unittest("random_rand_iint16", test_random_rand_iint8), &
            new_unittest("random_rand_iint32", test_random_rand_iint8), &
            new_unittest("random_rand_iint64", test_random_rand_iint8) &
        ]

    end subroutine collect_stats_random

    subroutine test_random_seed(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error

        integer :: put, get, res(5)
        integer :: ans(5) = [-1859553078, -1933696596, -642834430, &
                             1711399314,  1548311463]
        integer :: i

        put = 135792468
        do i = 1, 5
            call random_seed(put, get)
            res(i) = get
            put = get
        end do
        call check(error, all(res == ans))
    end subroutine test_random_seed

    subroutine test_random_rand_iint8(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error
        integer :: put, get, i
        integer(int8) :: res(5), ans(5) = [118, -15, -72, 101, 70]

        put = 12345678
        call random_seed(put, get)
        do i = 1, 5
            res(i) = dist_rand(1_int8)
        end do
        call check(error, all(res == ans))
    end subroutine test_random_rand_iint8

    subroutine test_random_rand_iint16(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error
        integer :: put, get, i
        integer(int16) :: res(5), ans(5) = [30286, -3799, -18204, 25947, 18148]

        put = 12345678
        call random_seed(put, get)
        do i = 1, 5
            res(i) = dist_rand(1_int16)
        end do
        call check(error, all(res == ans))
    end subroutine test_random_rand_iint16

    subroutine test_random_rand_iint32(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error
        integer :: put, get, i
        integer(int32) :: res(5), ans(5)=[1984865646, -248954393, -1192993267, &
                                  1700514835, 1189401802]

        put = 12345678
        call random_seed(put, get)
        do i = 1, 5
            res(i) = dist_rand(1_int32)
        end do
        call check(error, all(res == ans))
    end subroutine test_random_rand_iint32

    subroutine test_random_rand_iint64(error)
        !> Error handling
        type(error_type), allocatable, intent(out) :: error
        integer :: put, get, i
        integer(int64) :: res(5), ans(5)=[8524933037632333570_int64, &
                                  -1069250973542918798_int64,        &
                                  -5123867065024149335_int64,        &
                                  7303655603304982073_int64,         &
                                  5108441843522503546_int64]

        put = 12345678
        call random_seed(put, get)
        do i = 1, 5
            res(i) = dist_rand(1_int64)
        end do
        call check(error, all(res == ans))
    end subroutine test_random_rand_iint64

end module test_stats_random

program tester
    use iso_fortran_env, only: error_unit
    use testdrive, only: new_testsuite, run_testsuite, testsuite_type
    use test_stats_random, only: collect_stats_random
    implicit none

    integer :: stat, is
    type(testsuite_type), allocatable :: testsuites(:)
    character(len=*), parameter :: fmt = '("#", *(1x, a))'

    stat = 0

    testsuites = [ &
        new_testsuite("stats_random", collect_stats_random) &
        ]

    do is = 1, size(testsuites)
        write(error_unit, fmt) "Testing:", testsuites(is)%name
        call run_testsuite(testsuites(is)%collect, error_unit, stat)
    end do

    if (stat > 0) then
        write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!"
        error stop
    end if

end program tester

max_stars_repo_name: freevryheid/stdlib

max_issues_repo_name: seanwallawalla-forks/stdlib

max_forks_repo_name: seanwallawalla-forks/stdlib


PROGRAM find2d
   !
   !-----------------------------------------------------------------
   ! This program finds the crystal plane with the largest distance
   ! between atoms belong to different crystal planes, and then
   ! enlarge the distance and recalculate the basis vectors and atom
   ! position.
   !-----------------------------------------------------------------
   !
   USE find2d_com
   !
   IMPLICIT NONE
   !
   ! Data dictionary-------------------------------------------------
   !
   ! counters
   INTEGER :: ii,jj,kk,nn,mm
   !
   CALL input()
   CALL prep_inf()
   CALL find_gv()
   CALL find_av()
   !
ENDPROGRAM find2d

SUBROUTINE input()
   !
   USE find2d_com
   !
   IMPLICIT NONE
   INTEGER :: ii
   !
   ! read basis vectors of lattice, regardless of units
   !
   READ(*,*) cvec(:,1)
   READ(*,*) cvec(:,2)
   READ(*,*) cvec(:,3)
   !
   ! read atom positions
   !
   READ(*,*) nat
   ALLOCATE(atom_tau(3,nat))
   ALLOCATE(atom_type(nat))
   DO ii=1,nat
      READ(*,*) atom_type(ii), atom_tau(:,ii)
   ENDDO
   !
   ! read distances
   !
   READ(*,*) dist_min
   READ(*,*) dist_a_min
   READ(*,*) dist_a_new
   READ(*,*) r_max
   !
   ! read threshold
   !
   READ(*,*) eps_scale
   !
   CALL make_xsf(cvec(:,1),cvec(:,2),cvec(:,3),nat,atom_tau(:,:),atom_type(:),"xsf.inp.xsf")
ENDSUBROUTINE input

SUBROUTINE prep_inf()
   !
   USE find2d_com
   USE my_math
   USE phy_const
   !
   IMPLICIT NONE
   INTEGER :: ii
   !
   !
   WRITE(*,*) "find2D.x, version 1.0, written by Zhida Song"
   WRITE(*,*) "Information from input file--------------------------&
                                          --------------------------"
   !
   ! lattice basis vector
   !
   WRITE(*,*) "The basis vectors of lattice:"
   DO ii=1,3
      WRITE(*,100) ii, cvec(:,ii)
      100 FORMAT( T5,'a(',I1,')=(', 3ES16.8,'   )' )
   ENDDO
   rvolume = mix_prod(cvec(:,1),cvec(:,2),cvec(:,3))
   rradius = (3*rvolume/4/const_pi)**(1.0/3.0)
   WRITE(*,102) "Volume of a primitive cell is                    ",&
                                                              rvolume
   WRITE(*,102) "Wigner-seitz radius is                           ",&
                                                              rradius
   !
   ! Reciprocal lattice basis vector
   !
   bvec_c(:,1) = 2*const_pi*x_prod(cvec(:,2),cvec(:,3))/rvolume
   bvec_c(:,2) = 2*const_pi*x_prod(cvec(:,3),cvec(:,1))/rvolume
   bvec_c(:,3) = 2*const_pi*x_prod(cvec(:,1),cvec(:,2))/rvolume
   gvolume = 8*const_pi**3/rvolume
   gradius = (3*gvolume/4/const_pi)**(1.0/3.0)
   !
   WRITE(*,*) "The basis vectors of reciprocal lattice:"
   DO ii=1,3
      WRITE(*,103) ii, bvec_c(1,ii)/2/const_pi, &
                       bvec_c(2,ii)/2/const_pi, &
                       bvec_c(3,ii)/2/const_pi
      103 FORMAT( T5,'b(',I1,')=(', 3ES16.8,'   ) * 2pi' )
   ENDDO
   WRITE(*,102) "Volume of a reciprocal primitive cell is         ",&
                                                              gvolume
   WRITE(*,102) "Reciprocal wigner-seitz radius is                ",&
                                                              gradius
   !
   ! atoms in unit cell
   !
   WRITE(*,*) "Atom positions:"
   DO ii=1,nat
      WRITE(*,101) atom_type(ii), ii, atom_tau(:,ii)
      101 FORMAT( T5, 'type=', I3,',   tau(', I3, ')=(', 3ES16.8, &
                                                            '   )' )
   ENDDO
   !
   ! computation parameter
   !
   WRITE(*,*) "Computation parameters:"
   WRITE(*,102) "Minimal distance between layers                  ",&
                                                             dist_min
   WRITE(*,102) "Minimal distance between atoms in differnt layers",&
                                                             dist_min
   WRITE(*,102) "New distance between atoms in differnt layers    ",&
                                                           dist_a_min
   WRITE(*,102) "Maximal magnitude of lattice vectors             ",&
                                                                r_max
   WRITE(*,102) "Tolerance (scale)                                ",&
                                                            eps_scale
   !
   eps_r = rradius*eps_scale
   eps_g = gradius*eps_scale
   !
   WRITE(*,102) "Tolerance (lattice)                              ",&
                                                                eps_r
   WRITE(*,102) "Tolerance (reciprocal lattice)                   ",&
                                                                eps_g
   102 FORMAT( T4,A50,':',ES16.8 )
   !
   ! reciprocal
   !
ENDSUBROUTINE

SUBROUTINE find_gv()
   !
   USE kind
   USE find2d_com
   USE my_math
   USE phy_const
   !
   IMPLICIT NONE
   !
   INTEGER :: cnt, ii, jj, kk
   REAL(KIND=DP), DIMENSION(3) :: gv_tmp1,gv_tmp2
   REAL(KIND=DP), ALLOCATABLE, DIMENSION(:,:) :: gv_vec_tmp
   !
   ALLOCATE(gv_vec_tmp(3,gv_max_num))
   !
   cnt=1
   DO ii=gv_max_n1,-gv_max_n1,-1
      DO jj=gv_max_n2,-gv_max_n2,-1
         DO kk=gv_max_n3,-gv_max_n3,-1
            gv_tmp1 = REAL(ii,KIND=DP)*bvec_c(:,1) + &
                      REAL(jj,KIND=DP)*bvec_c(:,2)
            gv_tmp2 = REAL(kk,KIND=DP)*bvec_c(:,3) + gv_tmp1
            !
            IF( vect_len(gv_tmp2) <= 2.0*const_pi/dist_min ) THEN
               gv_vec_tmp(:,cnt) = gv_tmp2(:)
               cnt = cnt + 1
            ENDIF
            !
         ENDDO ! kk
      ENDDO ! jj
   ENDDO ! ii
   !
   ! collect and output the g vectors
   !
   gv_num = cnt-1
   ALLOCATE(gv_vec(3,gv_num))
   ALLOCATE(gv_len(gv_num))
   gv_vec(:,:) = gv_vec_tmp(:,1:gv_num)
   DEALLOCATE(gv_vec_tmp)
   !
   WRITE(*,*) "All g-vectors satisfy |g| <= 2pi/dist_min------------&
                                          --------------------------"
   !
   IF (gv_num == gv_max_num) THEN
      WRITE(*,*) "    dist_min may be too small"
   ENDIF
   !
   DO cnt=1,gv_num
      gv_len(cnt) = vect_len(gv_vec(:,cnt))
      WRITE(*,104) cnt, gv_vec(1,cnt)/2/const_pi, &
                        gv_vec(2,cnt)/2/const_pi, &
                        gv_vec(3,cnt)/2/const_pi, &
                        2*const_pi/gv_len(cnt)
      104 FORMAT( T5,'g(',I4,')=(', 3ES12.4,')*2pi,  distance =',ES12.4 )
   ENDDO
   !
   ! delete parallel gvectors ---------------------------------------
   !
   CALL del_para_g()
   !
   !WRITE(*,*) "Non-parallel g-vectors:------------------------------&
   !                                       --------------------------"
   !
   !DO cnt=1,gv_num
   !   gv_len(cnt) = vect_len(gv_vec(:,cnt))
   !   WRITE(*,104) cnt, gv_vec(1,cnt)/2/const_pi, &
   !                     gv_vec(2,cnt)/2/const_pi, &
   !                     gv_vec(3,cnt)/2/const_pi, &
   !                     2*const_pi/gv_len(cnt)
   !ENDDO
   !
   ! delete equilavent gvectors--------------------------------------
   !
   CALL del_eq_g()
   !
   WRITE(*,*) "Non-equivalent g-vectors:----------------------------&
                                          --------------------------"
   !
   DO cnt=1,gv_num
      gv_len(cnt) = vect_len(gv_vec(:,cnt))
      WRITE(*,104) cnt, gv_vec(1,cnt)/2/const_pi, &
                        gv_vec(2,cnt)/2/const_pi, &
                        gv_vec(3,cnt)/2/const_pi, &
                        2*const_pi/gv_len(cnt)
   ENDDO
   !
   DEALLOCATE(gv_len)
   !
ENDSUBROUTINE

SUBROUTINE del_para_g()
   !
   USE kind
   USE find2d_com
   USE my_math
   USE phy_const
   !
   IMPLICIT NONE
   !
   INTEGER :: ii,jj,cnt
   REAL(KIND=DP), ALLOCATABLE, DIMENSION(:,:) :: gv_tmp
   REAL(KIND=DP), PARAMETER :: eps_tmp=1e-6
   LOGICAL, ALLOCATABLE, DIMENSION(:) :: neq
   !
   IF ( gv_num==1 ) THEN
      WRITE(*,*) "ERROR!!! no proper g-vector exists!!!"
   ENDIF
   !
   ALLOCATE(gv_tmp(3,gv_num))
   ALLOCATE(neq(gv_num))
   DO ii=1,gv_num
      neq(ii) = .TRUE.
   ENDDO
   !
   ! delete the parallel and equivalent g-vectors
   !
   DO ii=1,gv_num-1
      DO jj=ii+1,gv_num
         IF ( ABS(vect_len(gv_vec(:,ii))) < eps_tmp*gradius ) THEN
            neq(ii) = .FALSE.
         ELSEIF ( ABS(vect_len(gv_vec(:,jj))) < eps_tmp*gradius ) THEN
            neq(jj) = .FALSE.
         ELSE
            IF ( ABS( vect_len( x_prod(gv_vec(:,ii),gv_vec(:,jj)) ) )&
                                              < eps_tmp*gradius ) THEN
               !
               IF( ABS(vect_len(gv_vec(:,ii))) <=                   &
                                   ABS(vect_len(gv_vec(:,jj))) ) THEN
                  neq(jj) = .FALSE.
               ELSE
                  neq(ii) = .FALSE.
               ENDIF
               !
            ENDIF
         ENDIF
      ENDDO
   ENDDO
   !
   ! recollect the non-equivalent g-vectors
   !
   cnt=1
   DO ii=2,gv_num
      IF( neq(ii) ) THEN
         gv_tmp(:,cnt) = gv_vec(:,ii)
         cnt=cnt+1
      ENDIF
   ENDDO
   !
   gv_num = cnt-1
   DEALLOCATE(gv_vec)
   ALLOCATE(gv_vec(3,gv_num))
   gv_vec(:,:) = gv_tmp(:,1:gv_num)
   !
   DEALLOCATE(gv_tmp)
   DEALLOCATE(neq)
ENDSUBROUTINE del_para_g

SUBROUTINE del_eq_g()
   !
   USE kind
   USE find2d_com
   USE my_math
   USE phy_const
   !
   IMPLICIT NONE
   !
   REAL(KIND=DP), DIMENSION(3) :: rot    ! rotation direction
   REAL(KIND=DP) :: alpha                ! rotation angle
   REAL(KIND=DP), DIMENSION(3) :: r1,r2  ! tmp vector
   REAL(KIND=DP) :: n1,n2,n3
   LOGICAL :: neq
   INTEGER :: ii,jj,kk,cnt
   REAL(KIND=DP), ALLOCATABLE, DIMENSION(:,:) :: gv_tmp
   !
   ALLOCATE(gv_tmp(3,gv_num))
   !
   cnt=1
   DO ii=1,gv_num-1
      neq = .TRUE.
      !
      DO jj=ii+1, gv_num
         IF (ABS( vect_len(gv_vec(:,ii))-vect_len(gv_vec(:,jj)) )<eps_g) THEN
            !
            ! rotation
            !
            CALL rot_alpha(gv_vec(:,ii),gv_vec(:,jj),rot,alpha)
            DO kk=1,nat
               r1 = atom_tau(:,kk)
               r2 = (1-cos(alpha))*dot_prod(r1,rot)*rot &
                    + sin(alpha)*x_prod(rot,r1) &
                    + cos(alpha)*r1
               r2 = r2 - r1
               n1 = dot_prod(r2,cvec(:,1))/dot_prod(cvec(:,1),cvec(:,1))
               n2 = dot_prod(r2,cvec(:,2))/dot_prod(cvec(:,2),cvec(:,2))
               n3 = dot_prod(r2,cvec(:,3))/dot_prod(cvec(:,3),cvec(:,3))
               IF ( ABS(NINT(n1)-n1)<eps_scale .AND. &
                    ABS(NINT(n2)-n2)<eps_scale .AND. &
                    ABS(NINT(n3)-n3)<eps_scale ) THEN
                  neq = .FALSE.
               ENDIF
            ENDDO
            !
            ! rotation-reversion
            !
            !CALL rot_alpha(gv_vec(:,ii),-gv_vec(:,jj),rot,alpha)
            !DO kk=1,nat
            !   r1 = atom_tau(:,kk)
            !   r2 = -(1-cos(alpha))*dot_prod(r1,rot)*rot &
            !        - sin(alpha)*x_prod(rot,r1) &
            !        - cos(alpha)*r1
            !   r2 = r2 - r1
            !   n1 = dot_prod(r2,cvec(:,1))/dot_prod(cvec(:,1),cvec(:,1))
            !   n2 = dot_prod(r2,cvec(:,2))/dot_prod(cvec(:,2),cvec(:,2))
            !   n3 = dot_prod(r2,cvec(:,3))/dot_prod(cvec(:,3),cvec(:,3))
            !   IF ( ABS(NINT(n1)-n1)<eps_scale .AND. &
            !        ABS(NINT(n2)-n2)<eps_scale .AND. &
            !        ABS(NINT(n3)-n3)<eps_scale ) THEN
            !      neq = .FALSE.
            !   ENDIF
            !ENDDO
            !
         ENDIF
      ENDDO ! jj
      !
      IF(neq) THEN
         gv_tmp(:,cnt) = gv_vec(:,ii)
         cnt = cnt + 1
      ENDIF
   ENDDO ! ii
   !
   gv_tmp(:,cnt) = gv_vec(:,gv_num)
   gv_num = cnt
   !
   DEALLOCATE(gv_vec)
   ALLOCATE(gv_vec(3,gv_num))
   gv_vec(:,:) = gv_tmp(:,:)
   DEALLOCATE(gv_tmp)
   !
ENDSUBROUTINE del_eq_g

SUBROUTINE rot_alpha(g1,g2,rot,alpha)
   !
   USE kind
   USE my_math
   USE phy_const
   USE find2d_com, ONLY : eps_g
   !
   IMPLICIT NONE
   !
   REAL(KIND=DP), DIMENSION(3), INTENT(IN) :: g1, g2
   REAL(KIND=DP), DIMENSION(3), INTENT(OUT) :: rot
   REAL(KIND=DP), INTENT(OUT) :: alpha
   REAL(KIND=DP), DIMENSION(3) :: gtmp
   !
   alpha = ACOS( dot_prod(g1(:),g2(:))/SQRT( dot_prod(g1(:),g1(:)) * &
                                             dot_prod(g2(:),g2(:)) ) )
   !
   rot = x_prod( g1(:),g2(:) )
   rot = rot/ABS(vect_len(rot))
   !
   gtmp = SIN(alpha)*x_prod(rot(:), g1(:)) + COS(alpha)*g1(:)
   !
   IF( vect_len( g2-gtmp ) > eps_g) &
      rot = (/ -rot(1), -rot(2), -rot(3) /)
ENDSUBROUTINE rot_alpha

SUBROUTINE find_av()
   !
   USE kind
   USE my_math
   USE phy_const
   USE find2d_com, ONLY : avec, gv_num, gv_vec
   !
   IMPLICIT NONE
   !
   INTEGER :: ii, jj
   LOGICAL, ALLOCATABLE, DIMENSION(:) :: finds
   !
   ALLOCATE(avec(3,3,gv_num))
   ALLOCATE(finds(gv_num))
   !
   DO ii=1,gv_num
      CALL find_av_for_g(ii,finds(ii))
   ENDDO
   !
   !
   DO ii=1,gv_num
      !
      WRITE(*,125) ii
      125 FORMAT(//,T2,"Configuration ",I3)
      !
      WRITE(*,120) ii, gv_vec(:,ii)/2/const_pi
      120 FORMAT( T5,'g(',I4,')=(', 3ES12.4,' ) * 2pi')
      !
      IF (finds(ii)) THEN
         WRITE(*,121) 'Find inner and intra layer basis vectors:'
         121 FORMAT( T8,A41)
         !
         DO jj=1,3
            WRITE(*,122) jj, avec(:,jj,ii)
            122 FORMAT( T8,'a(',I1,')=(', 3ES12.4,' )')
         ENDDO
         !
         WRITE(*,123) abs(dot_prod(x_prod(avec(:,1,ii),avec(:,2,ii)), &
                                          avec(:,3,ii)) )
         123 FORMAT(T8,'Volume of unit cell:',ES12.4)
         CALL layer_cryst(ii)
      ELSE
         WRITE(*,124) 'Cannot find inner and intra layer basis vectors:'
         124 FORMAT( T8,A48)
      ENDIF
   ENDDO
   !
   DEALLOCATE(finds)
ENDSUBROUTINE find_av

SUBROUTINE find_av_for_g(nn,finds)
   !
   USE kind
   USE my_math
   USE phy_const, ONLY : const_pi
   USE find2d_com, ONLY : avec, cvec, gv_vec, rradius
   !
   IMPLICIT NONE
   !
   INTEGER, INTENT(IN) :: nn                  ! g-vector index
   LOGICAL, INTENT(OUT) :: finds              ! finds=.true. if sucess
   !
   REAL(KIND=DP), DIMENSION(3) :: v1, v2, v3
   REAL(KIND=DP) :: sin_tmp
   REAL(KIND=DP), PARAMETER :: eps_sin = 1.0D-6
   LOGICAL :: suces
   !
   CALL find_inner_a(gv_vec(:,nn),cvec(:,1),cvec(:,2),cvec(:,3),v1,suces)
   !
   ! if we get v1, get v2
   IF (suces) THEN
      CALL find_inner_a(gv_vec(:,nn),cvec(:,2),cvec(:,3),cvec(:,1),v2,suces)
      !
      ! if we get v2, check wether v1 // v2
      IF (suces) THEN
         sin_tmp = abs(vect_len(x_prod(v1,v2)))/(vect_len(v1)*vect_len(v2))
         !
         ! if v1 // v2, get v2 again
         IF ( abs(sin_tmp)<eps_sin ) THEN
            CALL find_inner_a(gv_vec(:,nn),cvec(:,3),cvec(:,1),cvec(:,2),v2,suces)
            !
            ! if we also cannot get v2, print an error
            IF ( .NOT. suces) THEN
               finds = .FALSE.
            ELSE
               sin_tmp = abs(vect_len(x_prod(v1,v2)))/(vect_len(v1)*vect_len(v2))
               !
               ! if v1 // v2, print an error
               IF ( abs(sin_tmp)<eps_sin ) THEN
                  finds = .FALSE.
               ELSE
                  finds = .TRUE.
               ENDIF
            ENDIF
         ELSE
            finds = .TRUE.
         ENDIF
      !
      ! if we do not get v2, get v2 again
      ELSE
         CALL find_inner_a(gv_vec(:,nn),cvec(:,3),cvec(:,1),cvec(:,2),v2,suces)
         !
         ! if we also cannot get v2, print an error
         IF ( .NOT. suces) THEN
            finds = .FALSE.
         ELSE
            sin_tmp = abs(vect_len(x_prod(v1,v2)))/(vect_len(v1)*vect_len(v2))
            !
            ! if v1 // v2, print an error
            IF ( abs(sin_tmp)<eps_sin ) THEN
               finds = .FALSE.
            ELSE
               finds = .TRUE.
            ENDIF
         ENDIF
      ENDIF
   !
   ! if we do not get v1, get v1 again
   ELSE
      CALL find_inner_a(gv_vec(:,nn),cvec(:,2),cvec(:,3),cvec(:,1),v1,suces)
      !
      ! if we also cannot get v1, print an error
      IF ( .NOT. suces) THEN
         finds = .FALSE.
      !
      ! if we get v1, get v2
      ELSE
         CALL find_inner_a(gv_vec(:,nn),cvec(:,3),cvec(:,1),cvec(:,2),v2,suces)
         !
         ! if we also cannot get v2, print an error
         IF ( .NOT. suces) THEN
            finds = .FALSE.
         ELSE
            sin_tmp = abs(vect_len(x_prod(v1,v2)))/(vect_len(v1)*vect_len(v2))
            !
            ! if v1 // v2, print an error
            IF ( abs(sin_tmp)<eps_sin ) THEN
               finds = .FALSE.
            ELSE
               finds = .TRUE.
            ENDIF
         ENDIF
      ENDIF
   ENDIF
   !
   IF ( finds ) THEN
      CALL find_intra_a(gv_vec(:,nn),cvec(:,1),cvec(:,2),cvec(:,3),v3,suces)
      IF (.NOT. suces) THEN
         CALL find_intra_a(gv_vec(:,nn),cvec(:,2),cvec(:,3),cvec(:,1),v3,suces)
         IF (.NOT. suces) THEN
            CALL find_intra_a(gv_vec(:,nn),cvec(:,3),cvec(:,1),cvec(:,2),v3,suces)
            IF (.NOT. suces) THEN
               finds = .FALSE.
            ELSE
               finds = .TRUE.
            ENDIF
         ELSE
            finds = .TRUE.
         ENDIF
      ELSE
         finds = .TRUE.
      ENDIF
   ENDIF
   avec(:,1,nn) = v1
   avec(:,2,nn) = v2
   avec(:,3,nn) = v3
   !
ENDSUBROUTINE find_av_for_g

SUBROUTINE find_inner_a(gv,ain1,ain2,ain3,aout,suces)
   !
   USE kind
   USE phy_const, ONLY : const_pi
   USE my_math
   !
   IMPLICIT NONE
   !
   REAL(KIND=DP), DIMENSION(3), INTENT(IN) :: gv, ain1, ain2, ain3
   REAL(KIND=DP), DIMENSION(3), INTENT(OUT) :: aout
   LOGICAL, INTENT(OUT) :: suces
   INTEGER :: near_to
   !
   REAL(KIND=DP), DIMENSION(3,9) :: vec_tmp
   REAL(KIND=DP), DIMENSION(9)   :: cos_tmp
   REAL(KIND=DP) :: cos_tmp2
   REAL(KIND=DP), PARAMETER :: eps_cos = 1.0D-6
   !
   INTEGER :: cnt, ii, jj
   !
   aout = ain1
   suces = .FALSE.
   DO cnt=1,50
      vec_tmp(:,1) = aout
      vec_tmp(:,2) = aout + ain2
      vec_tmp(:,3) = aout - ain2
      vec_tmp(:,4) = aout + ain3
      vec_tmp(:,5) = aout - ain3
      vec_tmp(:,6) = aout + ain2 + ain3
      vec_tmp(:,7) = aout + ain2 - ain3
      vec_tmp(:,8) = aout - ain2 + ain3
      vec_tmp(:,9) = aout - ain2 - ain3
      !
      DO ii=1,9
         cos_tmp(ii) = dot_prod(vec_tmp(:,ii),gv) / &
                       (vect_len(vec_tmp(:,ii))*vect_len(gv))
      ENDDO
      !
      jj = near_to(REAL(0.0,KIND=DP),cos_tmp,9)
      cos_tmp2 = cos_tmp(jj)
      aout = vec_tmp(:,jj)
      !
      IF (abs(cos_tmp2)<eps_cos) THEN
         suces = .TRUE.
         EXIT
      ENDIF
      !
   ENDDO
ENDSUBROUTINE find_inner_a

SUBROUTINE find_intra_a(gv,ain1,ain2,ain3,aout,suces)
   !
   USE kind
   USE phy_const, ONLY : const_pi
   USE my_math
   !
   IMPLICIT NONE
   !
   REAL(KIND=DP), DIMENSION(3), INTENT(IN) :: gv, ain1, ain2, ain3
   REAL(KIND=DP), DIMENSION(3), INTENT(OUT) :: aout
   LOGICAL, INTENT(OUT) :: suces
   INTEGER :: near_to
   !
   REAL(KIND=DP), DIMENSION(3,9) :: vec_tmp
   REAL(KIND=DP), DIMENSION(9)   :: proj_tmp   ! a.g/|g|/|distance|
   REAL(KIND=DP) :: proj_tmp2
   REAL(KIND=DP) :: dist
   REAL(KIND=DP), PARAMETER :: eps_proj = 1.0D-6
   !
   INTEGER :: cnt, ii, jj
   !
   dist = 2*const_pi/vect_len(gv)
   !
   aout = ain1
   suces = .FALSE.
   DO cnt=1,50
      vec_tmp(:,1) = aout
      vec_tmp(:,2) = aout + ain2
      vec_tmp(:,3) = aout - ain2
      vec_tmp(:,4) = aout + ain3
      vec_tmp(:,5) = aout - ain3
      vec_tmp(:,6) = aout + ain2 + ain3
      vec_tmp(:,7) = aout + ain2 - ain3
      vec_tmp(:,8) = aout - ain2 + ain3
      vec_tmp(:,9) = aout - ain2 - ain3
      !
      DO ii=1,9
         proj_tmp(ii) = dot_prod(vec_tmp(:,ii),gv)/vect_len(gv)/dist
      ENDDO
      !
      jj = near_to(REAL(1.0,KIND=DP),proj_tmp,9)
      proj_tmp2 = proj_tmp(jj)
      aout = vec_tmp(:,jj)
      !
      IF (abs(proj_tmp2-1.0)<eps_proj) THEN
         suces = .TRUE.
         EXIT
      ENDIF
      !
   ENDDO
ENDSUBROUTINE find_intra_a

FUNCTION near_to(targ,angle,num)
   !
   USE kind
   USE phy_const, ONLY : const_pi
   !
   IMPLICIT NONE
   !
   INTEGER :: near_to
   REAL(KIND=DP), INTENT(IN) :: targ
   REAL(KIND=DP), DIMENSION(5), INTENT(IN) :: angle
   INTEGER, INTENT(IN) :: num
   REAL(KIND=DP) :: angle_tmp
   !
   INTEGER :: ii
   !
   near_to = 1
   angle_tmp = angle(1)
   DO ii=2,num
      IF ( abs(targ-angle(ii))<abs(targ-angle_tmp) )  THEN
         near_to = ii
         angle_tmp = angle(ii)
      ENDIF
   ENDDO
   !
ENDFUNCTION near_to


text:
C***********************************************************************
C
C       FUNCTION RAND1(IX)
C
C       This is the heart of all parts of the program that have anything
C       to do with random numbers. It was developed by L. Schrage and
C       taken from "Simulation Modeling and Analysis" by Averill M. Law.
C
C       VARIABLES:
C
C       IX  : Originally set to the initial seed and returns the next
C           integer, ZSUBI, which is the seed for the next call to
C           function RAND.
C       RAND: Contains the next value of USUBI, which is the uniform
C           random number being requested on the interval [0,1].
C
C***********************************************************************
C
        REAL FUNCTION RAND1(IX)
        IMPLICIT NONE
        INTEGER*4 A,B15,B16,FHI,IX,K,LEFTLO,P,XALO,XHI
        DATA A/16807/, B15/32768/, B16/65536/, P/2147483647/
        XHI = IX / B16
        XALO = (IX - XHI * B16) * A
        LEFTLO = XALO / B16
        FHI = XHI * A + LEFTLO
        K = FHI / B15
        IX = (((XALO-LEFTLO*B16)-P)+(FHI-K*B15)*B16) + K
        IF(IX.LT.0) IX = IX + P
        RAND1 = FLOAT(IX) * 4.656612875E-10
        RETURN
        END

max_stars_repo_name: phaustin/empm

max_issues_repo_name: phaustin/empm

max_forks_repo_name: phaustin/empm


text:
! { dg-do compile }
! Test for fix of PR27655
!
!Contributed by Francois-Xavier Coudert  <fxcoudert@gcc.gnu.org>
  integer, pointer :: i
  print *, associated(NULL(),i) ! { dg-error "not permitted as actual argument" }
  print *, associated(i,NULL()) ! { dg-error "not permitted as actual argument" }
end

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!                                                                      %
! Copyright (C) 2015 Adrian Martinez Vargas                            %
!                                                                      %
! This software may be modified and distributed under the terms        %
! of the MIT license.  See the LICENSE.txt file for details.           %
!                                                                      %
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!                                                                      %
! Copyright (C) 1996, The Board of Trustees of the Leland Stanford     %
! Junior University.  All rights reserved.                             %
!                                                                      %
! The programs in GSLIB are distributed in the hope that they will be  %
! useful, but WITHOUT ANY WARRANTY.  No author or distributor accepts  %
! responsibility to anyone for the consequences of using them or for   %
! whether they serve any particular purpose or work at all, unless he  %
! says so in writing.  Everyone is granted permission to copy, modify  %
! and redistribute the programs in GSLIB, but only under the condition %
! that this notice and the above copyright notice remain intact.       %
!                                                                      %
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!-----------------------------------------------------------------------
!
! The functions and subroutines below were modified from the
! version 2.0 of the gslib code written in fortran 77
!
! The objective is to add functionality to GSLIB and to link
! this code with python using f2py. It uses f90 code convention.
! The arrays are dynamic and externally declared (from python)
!
! the code was converted from Fortran 77 to Fortran 90 using F2F.pl
!
! for more information please refer to:
! - gslib77 source code: http://www.statios.com/software/gslib77_ls.tar.gz
! - GSLIB: Geostatistical Software Library and User's Guide. Second edition
!   by Clayton V. Deutsch, Andre G. Journel, 1997.
! - F2PY Users Guide: http://docs.scipy.org/doc/numpy-dev/f2py/
! - F2F https://bitbucket.org/lemonlab/f2f
!-----------------------------------------------------------------------

subroutine addcoord(nx, ny, nz, xmn, ymn, zmn, xsiz, ysiz, zsiz, xx, yy, zz)

    !
    !              Calculate coordinates of an ordered grid
    !              ****************************************
    ! INPUT VARIABLES:
    !
    !   nx, ny, nz           Number nodes in the grid
    !   xmn, ymn, zmn        Coordinate of the centroid (lower left corner)
    !   xsiz, ysiz, zsiz     Sizes of the cell/block
    !
    ! OUTPUT VARIABLES:
    !
    !   xx(n), yy(n), zz(n)  coordinate arrays with size n= nx*ny*nz
    !-----------------------------------------------------------------------

    real*8,  intent(in) :: xmn, ymn, zmn, xsiz, ysiz, zsiz
    integer, intent(in) :: nx, ny, nz
    real*8,    intent(out), dimension(nx*ny*nz) :: xx, yy, zz

    real*8 :: x, y, z
    integer :: i

    !
    ! Loop over grid
    !
    i=0
    do iz=1,nz
          z = zmn + dble(iz-1)*zsiz
          do iy=1,ny
                y = ymn + dble(iy-1)*ysiz
                do ix=1,nx
                      x = xmn + dble(ix-1)*xsiz
                      i=i+1
                      xx(i)=x
                      yy(i)=y
                      zz(i)=z
                end do
          end do
    end do

    return

end subroutine addcoord

max_stars_repo_name: cmrajan/pygslib

max_issues_repo_name: kaufmanno/pygslib

max_forks_repo_name: kaufmanno/pygslib


text:
  subroutine pint (fbuObj, p, ct, st, cf, sf, t, cm)

! ======================================================================
!
!    Called by: RASTAR
!
!   Defining angles from the 3-momentum components.
!
!    Last change: 13-Aug-2003 BY NVMokhov
!    Edited by A. J. Sierk, LANL T-16, September, 2003.
!    Edited by A. J. Sierk, LANL T-2, February, 2009.
!    Edited by AJS, LANL T-2, December, 2011.
!    Edited by LMK, XCP-3, July 2013 (included error protection).
!    Edited by CMJ, XCP-3, July 2018 (creation of FermiBreakup class)
!
! ======================================================================

    use, intrinsic :: iso_fortran_env, only: real64
    use fermiBreakupParams, only : zro, one

    implicit none
    class(FermiBreakup), intent(inout) :: fbuObj
    real(real64),        intent(in   ), dimension(3) :: p
    real(real64),        intent(  out) :: ct
    real(real64),        intent(  out) :: st
    real(real64),        intent(  out) :: cf
    real(real64),        intent(  out) :: sf
    real(real64),        intent(  out) :: t
    real(real64),        intent(in   ) :: cm

    real(real64) :: ctq, pmh, pq, pz, temp

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

    pz = p(3)**2
    pq = p(1)**2 + p(2)**2 + pz
    temp = pq
    if (temp < div0Lim .and. temp > -div0Lim) then
       temp = div0Lim
       write(fbuObj%io%message,1000) "175"
       call fbuObj%io%print(4, 3, fbuObj%io%message)
    end if
    ctq = pz/temp
    temp = pq + cm**2
    if (temp < 0.0d0) then
       temp = 0.01d0
       write(fbuObj%io%message,1000) "181"
       call fbuObj%io%print(4, 3, fbuObj%io%message)
    end if
    t = sqrt(temp) - cm
    if (ctq >= one) then
       st = zro
       ct = one
       sf = zro
       cf = one
    else
       temp = ctq
       if (temp < 0.0d0) then
          temp = 0.01d0
          write(fbuObj%io%message,1050) "193"
          call fbuObj%io%print(4, 3, fbuObj%io%message)
       end if
       ct = sqrt(temp)
       if (p(3) <= zro) ct = -ct
       st = sqrt(one - ctq)
       temp = pq
       if (temp < 0.0d0) then
          temp = 0.01d0
          write(fbuObj%io%message,1100) "201"
          call fbuObj%io%print(4, 3, fbuObj%io%message)
       end if
       pmh = st*sqrt(temp)
       temp = pmh
       if (temp < div0Lim .and. temp > -div0Lim) then
          temp = div0Lim
          write(fbuObj%io%message,1000) "207, 208"
          call fbuObj%io%print(4, 3, fbuObj%io%message)
       end if
       cf = p(1)/temp
       sf = p(2)/temp
    endif

    return
! ======================================================================
1000 format("Divide by zero error prevented in 'pint.f90', line ", A)
1050 format("Square root/divide by zero error prevented in ", &
          & "'pint.f90', line ", A)
1100 format("Square root error prevented in 'pint.f90', line ", A)
! ======================================================================
  end subroutine pint

max_stars_repo_name: lanl/generalized-spallation-model

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_name: lanl/generalized-spallation-model

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_name: lanl/generalized-spallation-model

max_forks_repo_licenses: Intel Unlicense