text:
!-----------------------------------------------------------------------
      subroutine ppiclf_user_SetYdot
!
      implicit none
!
      include "PPICLF"
!
! Internal:
!
      real*8 rmass,fbx,fby,fbz,fcx,fcy,fcz
      integer*4 i
!
! evaluate ydot
      do i=1,ppiclf_npart
         ! Particle mass
         rmass = ppiclf_rprop(PPICLF_R_JVOLP,i)
     >          *ppiclf_rprop(PPICLF_R_JRHOP,i)

         ! Gravity
         fbx  = 0.0d0
         fby  = -9.8d0*rmass
         fbz  = 0.0d0

         ! Collision search
         call ppiclf_solve_NearestNeighbor(i)

         ! User implemented collision force from EvalNearestNeighbor()
         fcx  = ppiclf_ydotc(PPICLF_JVX,i)
         fcy  = ppiclf_ydotc(PPICLF_JVY,i)
         fcz  = ppiclf_ydotc(PPICLF_JVZ,i)

         ! set ydot for all PPICLF_LRS number of equations
         ppiclf_ydot(PPICLF_JX ,i) = ppiclf_y(PPICLF_JVX,i)
         ppiclf_ydot(PPICLF_JY ,i) = ppiclf_y(PPICLF_JVY,i)
         ppiclf_ydot(PPICLF_JZ ,i) = ppiclf_y(PPICLF_JVZ,i)
         ppiclf_ydot(PPICLF_JVX,i) = (fbx+fcx)/rmass
         ppiclf_ydot(PPICLF_JVY,i) = (fby+fcy)/rmass
         ppiclf_ydot(PPICLF_JVZ,i) = (fbz+fcz)/rmass
      enddo
! evaluate ydot

      return
      end
!-----------------------------------------------------------------------
      subroutine ppiclf_user_MapProjPart(map,y,ydot,ydotc,rprop)
!
      implicit none
!
! Input:
!
      real*8 y    (PPICLF_LRS)
      real*8 ydot (PPICLF_LRS)
      real*8 ydotc(PPICLF_LRS)
      real*8 rprop(PPICLF_LRP)
!
! Output:
!
      real*8 map  (PPICLF_LRP_PRO)
!

      return
      end
!-----------------------------------------------------------------------
      subroutine ppiclf_user_EvalNearestNeighbor
     >                                        (i,j,yi,rpropi,yj,rpropj)
!
      implicit none
!
      include "PPICLF"
!
! Input:
!
      integer*4 i
      integer*4 j
      real*8 yi    (PPICLF_LRS)
      real*8 rpropi(PPICLF_LRP)
      real*8 yj    (PPICLF_LRS)
      real*8 rpropj(PPICLF_LRP)
!
! Internal:
!
      real*8 ksp,erest
      common /ucollision/ ksp,erest
#ifdef PPICLC
      BIND(C, name="ucollision") :: /ucollision/ ! c binding
#endif

      real*8 rpi2, rthresh, rxdiff, rydiff, rzdiff, rdiff, rm1, rm2,
     >       rmult, eta, rbot, rn_12x, rn_12y, rn_12z, rdelta12,
     >       rv12_mag, rv12_mage, rksp_max, rnmag, rksp_wall, rextra
!
      rpi2  =  9.869604401089358d0

      ! other particles
      if (j .ne. 0) then
         rthresh  = 0.5d0*(rpropi(PPICLF_R_JDP) + rpropj(PPICLF_R_JDP))

         rxdiff = yj(PPICLF_JX) - yi(PPICLF_JX)
         rydiff = yj(PPICLF_JY) - yi(PPICLF_JY)
         rzdiff = yj(PPICLF_JZ) - yi(PPICLF_JZ)

         rdiff = sqrt(rxdiff**2 + rydiff**2 + rzdiff**2)

         if (rdiff .gt. rthresh) return

         rm1 = rpropi(PPICLF_R_JRHOP)*rpropi(PPICLF_R_JVOLP)
         rm2 = rpropj(PPICLF_R_JRHOP)*rpropj(PPICLF_R_JVOLP)

         rmult = 1.0d0/sqrt(1.0d0/rm1+1.0d0/rm2)
         eta   = 2.0d0*sqrt(ksp)*log(erest)/sqrt(log(erest)**2+rpi2)
     >           *rmult

         rbot = 1.0d0/rdiff
         rn_12x = rxdiff*rbot
         rn_12y = rydiff*rbot
         rn_12z = rzdiff*rbot

         rdelta12 = rthresh - rdiff

         rv12_mag = (yj(PPICLF_JVX)-yi(PPICLF_JVX))*rn_12x +
     >              (yj(PPICLF_JVY)-yi(PPICLF_JVY))*rn_12y +
     >              (yj(PPICLF_JVZ)-yi(PPICLF_JVZ))*rn_12z

         rv12_mage = rv12_mag*eta
         rksp_max  = ksp*rdelta12
         rnmag     = -rksp_max - rv12_mage

         ppiclf_ydotc(PPICLF_JVX,i) = ppiclf_ydotc(PPICLF_JVX,i)
     >                              + rnmag*rn_12x
         ppiclf_ydotc(PPICLF_JVY,i) = ppiclf_ydotc(PPICLF_JVY,i)
     >                              + rnmag*rn_12y
         ppiclf_ydotc(PPICLF_JVZ,i) = ppiclf_ydotc(PPICLF_JVZ,i)
     >                              + rnmag*rn_12z

      ! boundaries
      elseif (j .eq. 0) then

         rksp_wall = ksp

         ! give a bit larger collision threshold for walls
         rextra   = 0.5d0
         rthresh  = (0.5d0+rextra)*rpropi(PPICLF_R_JDP)

         rxdiff = yj(PPICLF_JX) - yi(PPICLF_JX)
         rydiff = yj(PPICLF_JY) - yi(PPICLF_JY)
         rzdiff = yj(PPICLF_JZ) - yi(PPICLF_JZ)

         rdiff = sqrt(rxdiff**2 + rydiff**2 + rzdiff**2)

         if (rdiff .gt. rthresh) return

         rm1 = rpropi(PPICLF_R_JRHOP)*rpropi(PPICLF_R_JVOLP)

         rmult = sqrt(rm1)
         eta   = 2.0d0*sqrt(rksp_wall)*log(erest)
     >           /sqrt(log(erest)**2+rpi2)*rmult

         rbot = 1.0d0/rdiff
         rn_12x = rxdiff*rbot
         rn_12y = rydiff*rbot
         rn_12z = rzdiff*rbot

         rdelta12 = rthresh - rdiff

         rv12_mag = -1.0d0*(yi(PPICLF_JVX)*rn_12x +
     >                      yi(PPICLF_JVY)*rn_12y +
     >                      yi(PPICLF_JVZ)*rn_12z)

         rv12_mage = rv12_mag*eta
         rksp_max  = rksp_wall*rdelta12
         rnmag     = -rksp_max - rv12_mage

         ppiclf_ydotc(PPICLF_JVX,i) = ppiclf_ydotc(PPICLF_JVX,i)
     >                              + rnmag*rn_12x
         ppiclf_ydotc(PPICLF_JVY,i) = ppiclf_ydotc(PPICLF_JVY,i)
     >                              + rnmag*rn_12y
         ppiclf_ydotc(PPICLF_JVZ,i) = ppiclf_ydotc(PPICLF_JVZ,i)
     >                              + rnmag*rn_12z

      endif

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

max_stars_repo_name: dpzwick/ppiclF

max_issues_repo_name: dpzwick/ppiclF

max_forks_repo_name: dpzwick/ppiclF


text:
C Copyright (C) 2002, Carnegie Mellon University and others.
C All Rights Reserved.
C This code is published under the Common Public License.
C*******************************************************************************
C
      logical function CHECK_BASIS(ITER, M, C, PY, CONDC)
C
C*******************************************************************************
C
C    $Id: check_basis.f,v 1.1 2002/05/02 18:52:17 andreasw Exp $
C
C-------------------------------------------------------------------------------
C                                 Title
C-------------------------------------------------------------------------------
C
CT    Test for conditioning of basis matrix
C
C-------------------------------------------------------------------------------
C                          Programm description
C-------------------------------------------------------------------------------
C
CB
C
C-------------------------------------------------------------------------------
C                             Author, date
C-------------------------------------------------------------------------------
C
CA    Andreas Waechter      05/01/02  Release as version IPOPT 2.0
C
C-------------------------------------------------------------------------------
C                             Documentation
C-------------------------------------------------------------------------------
C
CD
C
C-------------------------------------------------------------------------------
C                             Parameter list
C-------------------------------------------------------------------------------
C
C    Name     I/O   Type   Meaning
C
CP   ITER      I    INT    iteration counter
CP   M         I    INT    number of constraints
CP   C         I    DP     values of constraints
CP   PY        I    DP     range space step
CP   CONDC     O    DP     estimated condition number of C
CP                            (only for output; is stored internally)
CP   CHECK_BASIS O  LOG    new basis required if true
C
C-------------------------------------------------------------------------------
C                             local variables
C-------------------------------------------------------------------------------
C
CL
C
C-------------------------------------------------------------------------------
C                             used subroutines
C-------------------------------------------------------------------------------
C
CS    DNRM2
CS    C_OUT
C
C*******************************************************************************
C
C                              Declarations
C
C*******************************************************************************
C
      IMPLICIT NONE
C
C*******************************************************************************
C
C                              Include files
C
C*******************************************************************************
C
      include 'IPOPT.INC'
C
C-------------------------------------------------------------------------------
C                             Parameter list
C-------------------------------------------------------------------------------
C
      integer ITER
      integer M
      double precision C(M)
      double precision PY(M)
      double precision CONDC
C
C-------------------------------------------------------------------------------
C                            Local varibales
C-------------------------------------------------------------------------------
C
      double precision RLARGE, RSMALL
      save             RLARGE, RSMALL

      double precision DNRM2

      double precision pynrm, ratio, cnrm, quot
      character*80 line(2)
C
C*******************************************************************************
C
C                           Executable Statements
C
C*******************************************************************************
C
      if( ITER .eq. 0 ) then
         RLARGE = 1.d-20
         RSMALL = 1.d20

         CHECK_BASIS = .false.
         CONDC = 0.d0
      endif

      pynrm = DNRM2(M, PY, 1)
      cnrm  = DNRM2(M, C , 1)

      if( pynrm .gt. 0 ) then
         ratio = cnrm/pynrm
      else
         ratio = RLARGE
      endif

      if( ITER.gt.0 .and. pynrm.gt.0 ) then
         RLARGE = dmax1(ratio, RLARGE)
         RSMALL = dmin1(ratio, RSMALL)

         quot = RLARGE/RSMALL
         write(line,700) ITER, quot
 700     format(/,' In Iter ',i5,
     1        ' Check Basis with RLARGE/RSMALL = ',
     1        d11.3)
         call C_OUT(1,2,2,line)

         CONDC = RLARGE/RSMALL
         CHECK_BASIS = (CONDC .gt. QMAXCOND)

         if( CHECK_BASIS ) then
            RLARGE = 1.d-20
            RSMALL = 1.d20
         endif

      else

         CHECK_BASIS = .false.

      endif

      return
      end

ext: f

max_stars_repo_name: kant/CMMPPT

max_issues_repo_name: kant/CMMPPT

max_forks_repo_name: kant/CMMPPT


text:
This editor can edit this entry and tell us a bit about themselves by clicking the Edit icon.

20090814 21:45:58 nbsp  Welcome to the Wiki. You can establish identity here in a variety of ways. Doing so will give more weight to your comment on the Alvarado Sunset Apartments page. Youve already chosen your Importance of using your RealName username, so the next step would be to introduce yourself and make some positive contributions like adding business reviews, Photo Requests uploading photos, or adding to Seed Pages List pages that need more information. Users/JasonAller

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
C$Procedure      SPHLAT ( Spherical to latitudinal coordinates )

      SUBROUTINE SPHLAT ( R, COLAT, LONGS, RADIUS, LONG, LAT )

C$ Abstract
C
C     Convert from spherical coordinates to latitudinal coordinates.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C     None.
C
C$ Keywords
C
C     CONVERSION, COORDINATES
C
C$ Declarations

      DOUBLE PRECISION R
      DOUBLE PRECISION COLAT
      DOUBLE PRECISION LONGS
      DOUBLE PRECISION RADIUS
      DOUBLE PRECISION LONG
      DOUBLE PRECISION LAT

C$ Brief_I/O
C
C     VARIABLE  I/O  DESCRIPTION
C     --------  ---  --------------------------------------------------
C     R          I   Distance of the point from the origin.
C     COLAT      I   Angle of the point from positive Z axis (radians).
C     LONGS      I   Angle of the point from the XZ plane (radians).
C     RADIUS     O   Distance of a point from the origin
C     LONG       O   Angle of the point from the XZ plane in radians
C     LAT        O   Angle of the point from the XY plane in radians
C
C$ Detailed_Input
C
C     R          Distance of the point from the origin.
C
C     COLAT      Angle between the vector from the origin to the point
C                and the positive Z axis in radians.
C
C     LONGS      Angle of the point from the XZ plane (radians).
C
C$ Detailed_Output
C
C     RADIUS     Distance of a point from the origin
C
C     LONG       Angle of the point from the XZ plane in radians. LONG
C                is set equal to LONGS.
C
C     LAT        Angle of the point from the XY plane in radians. LAT
C                is computed as pi/2 - COLAT.
C
C$ Parameters
C
C     None.
C
C$ Exceptions
C
C     Error free.
C
C$ Files
C
C     None.
C
C$ Particulars
C
C     This routine returns the latitudinal coordinates of a point
C     whose position is input in spherical coordinates.
C
C     Latitudinal coordinates are defined by a distance from a central
C     reference point, an angle from a reference meridian, and an angle
C     above the equator of a sphere centered at the central reference
C     point.
C
C     Spherical coordinates are defined by a distance from a central
C     reference point, an angle from a reference meridian, and an angle
C     from the z-axis.
C
C$ Examples
C
C
C     Latitude is obtained by subtracting co-latitude from HALFPI()
C     Radius and longitude mean the same thing in both latitudinal
C     and spherical coordinates.  The table below lists LAT and
C     corresponding COLAT in terms of degrees.
C
C             LAT            COLAT
C            ------         ------
C              0             90
C             20             70
C             45             45
C            -30            120
C             90              0
C            -45            135
C
C$ Restrictions
C
C     None.
C
C$ Literature_References
C
C     None.
C
C$ Author_and_Institution
C
C     W.L. Taber      (JPL)
C
C$ Version
C
C-    SPICELIB Version 1.0.2, 26-JUL-2016 (BVS)
C
C        Minor headers edits.
C
C-    SPICELIB Version 1.0.1, 10-MAR-1992 (WLT)
C
C        Comment section for permuted index source lines was added
C        following the header.
C
C-    SPICELIB Version 1.0.0, 31-JAN-1990 (WLT)
C
C-&

C$ Index_Entries
C
C     spherical to latitudinal coordinates
C
C-&

C$ Revisions
C
C-    Beta Version 1.0.1, 1-Feb-1989 (WLT)
C
C        Example section of header upgraded.
C
C-&

C
C     SPICELIB functions
C
      DOUBLE PRECISION HALFPI

C
C     Local Variables
C
      DOUBLE PRECISION RR
      DOUBLE PRECISION LATTUD

C
C     Convert to latitudinal coordinates, storing the results in
C     temporary variables
C
      RR     = R
      LATTUD = HALFPI() - COLAT

C
C     Move the results to the output variables.
C
      LONG   = LONGS
      RADIUS = RR
      LAT    = LATTUD

      RETURN
      END

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
!
! MONORTM_r3p1_Module
!
! Module containing procedures for MONORTM Record r3p1.
!
!
! CREATION HISTORY:
!       Written by:   Paul van Delst, 24-Dec-2012
!                     paul.vandelst@noaa.gov
!

MODULE MONORTM_r3p1_Module

  ! -----------------
  ! Environment setup
  ! -----------------
  ! Module usage
  USE Type_Kinds     , ONLY: fp
  USE File_Utility   , ONLY: File_Open
  USE Message_Handler, ONLY: SUCCESS, FAILURE, INFORMATION, Display_Message
  ! Line-by-line model parameters
  USE LBL_Parameters
  ! Disable implicit typing
  IMPLICIT NONE

  ! ----------
  ! Visibility
  ! ----------
  ! Everything private by default
  PRIVATE
  ! Datatypes
  PUBLIC :: MONORTM_r3p1_type
  ! Procedures
  PUBLIC :: MONORTM_r3p1_Write

  ! -----------------
  ! Module parameters
  ! -----------------
  CHARACTER(*), PARAMETER :: MODULE_VERSION_ID = &
  ! Message string length
  INTEGER, PARAMETER :: ML = 256
  ! The record I/O format
  CHARACTER(*), PARAMETER :: MONORTM_R3P1_FMT = '(7i5,i2,1x,i2,3f10.3,10x,f10.3)'

  ! -------------
  ! Derived types
  ! -------------
  TYPE :: MONORTM_r3p1_type
    INTEGER  :: model   = 0       !  Selects atmospheric profile
    INTEGER  :: itype   = 0       !  Selects type of path
    INTEGER  :: ibmax   = 0       !  Selects layering
    INTEGER  :: nozero  = 1       !  Controls zeroing of small absorber amounts.
    INTEGER  :: noprnt  = 0       !  Controls printout length
    INTEGER  :: nmol    = 0       !  Number of molecular species
    INTEGER  :: ipunch  = 1       !  Controls output of layer data to TAPE7
    INTEGER  :: ifxtyp  = 0       !  Controls cross-layer DV field output
    INTEGER  :: munits  = 0       !  Controls molecular amounts units output
    REAL(fp) :: re      = 0.0_fp  !  Radius of earth (km). 0.0 == model defaults
    REAL(fp) :: hspace  = 0.0_fp  !  Altitude defintion for space. 100km == default
    REAL(fp) :: vbar    = 0.0_fp  !  Frequency for refractive geometry. Average == default
    REAL(fp) :: ref_lat = 0.0_fp  !  LAtitude of location of calculation. 0.0 == model defaults
  END TYPE MONORTM_r3p1_type

CONTAINS

  FUNCTION MONORTM_r3p1_Write(r3p1,fid) RESULT(err_stat)

    ! Arguments
    TYPE(MONORTM_r3p1_type), INTENT(IN) :: r3p1
    INTEGER                , INTENT(IN) :: fid
    ! Function result
    INTEGER :: err_stat
    ! Function parameters
    CHARACTER(*), PARAMETER :: ROUTINE_NAME = 'MONORTM_r3p1_Write'
    ! Function variables
    CHARACTER(ML) :: msg
    CHARACTER(ML) :: io_msg
    INTEGER :: io_stat

    ! Setup
    err_stat = SUCCESS
    ! ...Check unit is open
    IF ( .NOT. File_Open(fid) ) THEN
      msg = 'File unit is not connected'
      CALL Cleanup(); RETURN
    END IF

    ! Write the record
    WRITE( fid,FMT=MONORTM_R3P1_FMT,IOSTAT=io_stat,IOMSG=io_msg) r3p1
    IF ( io_stat /= 0 ) THEN
      msg = 'Error writing record - '//TRIM(io_msg)
      CALL Cleanup(); RETURN
    END IF

  CONTAINS

    SUBROUTINE CleanUp()
      err_stat = FAILURE
      CALL Display_Message( ROUTINE_NAME,msg,err_stat )
    END SUBROUTINE CleanUp

  END FUNCTION MONORTM_r3p1_Write

END MODULE MONORTM_r3p1_Module

max_stars_repo_name: hsbadr/crtm

max_issues_repo_name: hsbadr/crtm

max_forks_repo_name: hsbadr/crtm


!!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
! TYPE  : MODULE
! NAME  : ED_Subspace
! OBJECT: TYPE(EDSubSpace)
! USED  : fermion_table,FermionHamiltonian,MatrixOperations
! DATE  : 2017-11-28
! AUTHOR: hengyueli@gmail.com
!--------------
! Open-Source : No
!------------------
! DESCRIPTION:
!            consider a subspace in ED. This subspace is identified by a subid which is correpsonding to that in table.
!
! STANDARD:
!            *CALL Initialization( )
!             call diagonalization()
!
! USING LIST:
!            :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!
!
! avalable sets:
!                  [sub] Initialization(Ta,IsReal,H,subid,*PRINT)
!                        type(table)::Ta  (will be pointed and be read only)
!                        logical    ::IsReal ! if H is a real problem.
!                        type(Ham)  ::H , will be  pointed.
!                        integer    ::subid
!                        Ta shoule be initiated before this subroutine since some information would be used.
!                  [sub] NormolizedEnergy(Eg)
!                        real(8)::Eg
!                        E -> E -Eg
!                  [sub] diagonalization()
!
! avalable gets:
!
!                  [fun] get_E(i)
!                        integer::i   ;   real(8)::get_E
!                  [fun] get_Vec(i)
!                        integer(8)::i     ;complex*16::get_Vec(d)  !where d is the dimension of the subspace.
!                  [fun] get_act(opt,i,d)
!                        type(FermOper)::opt
!                        integer(8)::i
!                        integer(8)::d
!                        complex*16::get_act(d)
!                        for a input operator opt%act
!                        opt%act |i> = S(d)    ,  thus the dimension of outputstate should be input.
!                  [fun] get_product(S,i)
!                        complex*16::S(self%d),get_product
!                        integer::i to identify the i-th eignstate in this subspace.
!                        returnV = <i|S> .      The input array S should be the same dimension of this subspace.
!
!
! avalable is :
!                  [fun] i
! others      :
!                  [sub] printE()
!                        print out eigenenergy
!
!
!
!!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

module ED_Subspace
   use fermion_table
   use FermionHamiltonian
   use MatrixOperations
   implicit none

   TYPE::EDSubSpace
       private
       logical                ::initiated  = .false.
       logical                ::Normolized = .false.
       logical                ::IsReal     = .false.
       integer                ::subid      = -1_8
       integer                ::d          = -1_8
       class(table),pointer   ::T          => null()
       class(Ham),pointer     ::Hi         => NULL()
       complex*16,allocatable ::Hs(:,:)
       real*8    ,allocatable ::E(:)
       real*8::Eg ! after normolized, this will save the ground state energy.
       integer::print = 6
   contains
       procedure,pass::Initialization
       procedure,pass::UnInitialization
       final::Finalization
       !
       procedure,pass::diagonalization
       procedure,pass::NormolizedEnergy
       procedure,pass::get_E
       procedure,pass::get_act
       procedure,pass::get_product
       procedure,pass::get_Vec
       !
       !
       procedure,pass::printE
   end type
   !
   private::Initialization,UnInitialization,Finalization
   !
   private::set_Hamiltonian_matrix_element,diagonalization,NormolizedEnergy,get_E
   private::get_act,get_product,printE,get_Vec
   !

 contains

  subroutine Initialization(self,Ta,IsReal,H,subid,print_)
    implicit none
    class(EDSubSpace),intent(inout)::self
    class(table),target,intent(inout)::Ta
    logical,intent(in)::IsReal
    class(Ham),target,intent(inout)::H
    integer,intent(in)::subid
    integer,intent(in),optional::print_
    !---------------------------------------
    call UnInitialization(self)
    self%initiated = .true.
    !----- check if T is initiated
    if (.not. ( Ta%Is_initiated()   ) )then
       write(self%print,*)"When initialize EDSubSpace, the input Table should be Initiated."
       stop
    endif
    self%T     => Ta
    self%IsReal= IsReal
    self%Hi    => H
    self%subid = subid
    self%d     = self%T%get_sub_d(self%subid)

    if (present(print_)) self%print = print_
    allocate(  self%Hs(self%d,self%d)   )
    allocate(  self%E(self%d)           )

  endsubroutine

  subroutine UnInitialization(self)
    implicit none
    class(EDSubSpace),intent(inout)::self
    !---------------------------------------
    if (self%initiated)then
       self%Initiated = .false.
       self%subid     = -1_8
       self%d         = -1_8
       self%T         => null()
       self%Hi        => NULL()
       deallocate( self%Hs  )
       deallocate( self%E   )
    endif
  endsubroutine

  subroutine Finalization(self)
    implicit none
    type(EDSubSpace),intent(inout)::self
    !---------------------------------------
    call UnInitialization(self)
  endsubroutine

  subroutine set_Hamiltonian_matrix_element(self)  ! only upper part of H will be setted.
    implicit none
    class(EDSubSpace),intent(inout)::self
    !---------------------------------------
    integer*8::basisIn,basisOut
    integer::jcopra,jci,jcj,Sign_,subid
    complex*16::v
    class(FermOper),pointer::opr
    self%Hs = (0._8,0._8)
    subid   = self%subid
    do jcopra = 1_8 , self%hi%GetOptN()
       opr    => self%hi%GetOptact(jcopra)
       v      =  self%Hi%GetOptV(jcopra)

                                        ! write(*,"(  (I2),',',(I2) )")jcopra,self%hi%GetOptN()
                                        ! write(*,"(I2)")opr%get_optid()
                                        ! write(*,*)opr%get_para()
                                        ! write(*,"(f6.3)")real(v)
                                        ! write(*,*)"-------------------"

       !------------------------------
       !    for ( opr  ,   v   )
       do jcj = 1 , self%d
          !-------- for subbasi  |jcs>---
          basisIn = self%T%get_subindex_to_basis(subid,jcj)
          call opr%act(basisIn,basisOut,Sign_)
          if (basisOut.ne.-1_8)then
             jci = self%T%get_basis_to_sub_index(basisOut)
             if ( jci.le.jcj  )then
               self%Hs(jci,jcj) = self%Hs(jci,jcj) + V * Sign_
             endif
          endif
       enddo
       !---------------------------------
   enddo
  endsubroutine

  subroutine diagonalization(self)
    implicit none
    class(EDSubSpace),intent(inout)::self
    !---------------------------------------
    TYPE(nummethod)::num
    real*8,allocatable::Hr(:,:)

    call set_Hamiltonian_matrix_element(self)

    !--------------
    if (self%IsReal)then
      allocate(Hr(self%d,self%d))
      Hr = Real( self%Hs ,8 )
      call num%ED_Hermitian_matrix("U",self%d,Hr     ,self%E)
      self%Hs = cmplx(Hr,0._8,8)
      deallocate( Hr  )
    else
      call num%ED_Hermitian_matrix("U",self%d,self%Hs,self%E)
    endif
    !--------------
    self%Normolized = .false.

  endsubroutine

  subroutine NormolizedEnergy(self,Eg)
            implicit none
            class(EDSubSpace),intent(inout)::self
            real(8),intent(in)::Eg
            !---------------------------------------
            integer(8)::jc
            if (.Not.self%normolized)then
               do jc = 1_8 , self%d
                  self%E(jc) = self%E(jc) -Eg
               enddo
               self%Eg = Eg
               self%normolized = .true.
            endif
  endsubroutine

  function get_E(self,i) result(E)
            implicit none
            class(EDSubSpace),intent(inout)::self
            integer,intent(in)::i
            real*8::E
            !---------------------------------------
            E = self%E(i)
  endfunction

  function get_Vec(self,i) result(v)
    implicit none
    class(EDSubSpace),intent(inout)::self
    integer,intent(in)::i
    complex*16::V(self%d)
    !---------------------------------------
    V = self%Hs(:,i)
endfunction

  function get_act(self,opt,i,d)  result(s)
            implicit none
            class(EDSubSpace),intent(inout)::self
            class(FermOper),intent(inout)::opt
            integer,intent(in)::i,d
            complex*16::s(d)
            !---------------------------------------
            integer::jc,sign_,id,subid
            integer*8::bi,bo
            subid = self%subid
            s = (0._8,0._8)
            do jc = 1 ,self%d
               bi = self%T%get_subindex_to_basis(subid,jc)
               call opt%act(bi,bo,sign_ )
               if ( bo.ne.-1_8 )then
                  id =  self%T%get_basis_to_sub_index(bo)
                  s(id) = s(id) + sign_ * self%Hs(jc,i)
               endif
            enddo
  endfunction

  function get_product(self,S,i) result(returnV)
            implicit none
            class(EDSubSpace),intent(inout)::self
            complex*16,intent(in)::s(self%d)
            integer,intent(in)::i
            complex*16::returnV
            !---------------------------------------
            returnV = sum( conjg(self%Hs(:,i))*S )
  endfunction

  subroutine printE(self)
            implicit none
            class(EDSubSpace),intent(inout)::self
            !---------------------------------------
            integer::jc
            write(self%print,*) self%T%get_subspace_marker(self%subid)
            do jc = 1 , self%d
              write(self%print,*)jc-1,self%e(jc)
            enddo
            ! write(self%print,*) self%E
  endsubroutine

end module


module euler_multiprecision_util_m

    implicit none
    private
    public :: add, sub, compare, mul, pow2

contains

    !> The core function to cut leading zeros of a contiguous array.
    !>
    !> As the name of the function suggests, the function
    !> deletes all the leading zeros of a contiguous array.
    !> This procedure uses the intrinsic function `findloc`
    !> that is only available with gfortran > 9.0.0
    pure function cut_leading_zeros(arr) result(ret)
        integer, contiguous, intent(in) :: arr(:)
        integer, allocatable :: ret(:)

        ret = arr(findloc(arr == 0, value=.false., dim=1):)
    end function cut_leading_zeros

    !> The core function to carry digits during add and multiplication.
    !>
    !> Each element of the input array is an 32bit integer greater than or
    !> equal to zero and less than or equal to nine.
    pure function carry(arr) result(ret)
        integer, contiguous, intent(in) :: arr(:)
        integer, dimension(size(arr) + 2) :: tmp1, tmp2, ret

        tmp1 = 0; tmp2 = 0; ret(1:2) = 0; ret(3:) = arr(:)
        do
            if (all(ret <= 9)) exit
            tmp1(:) = ret(:) - ret(:)/10*10
            tmp2(1:size(tmp2) - 1) = ret(2:)/10
            ret(:) = tmp1(:) + tmp2(:)
        end do
    end function carry

    !> The core function to provide the ability of addition.
    !>
    !> Each element of the input array is an 32bit integer greater than or
    !> equal to zero and less than or equal to nine. When two such arrays
    !> are "added", first of all, the shorter array will be extended to the
    !> length of the longer array and two arrays with the same dimension
    !> are added. Then the function `carry`
    !> "[carries](https://en.wikipedia.org/wiki/Carry_(arithmetic))" the
    !> added array iteratively until all elements of the array are greater
    !> than or equal to zero and less than or equal to nine. Finally, the
    !> leading zeros are deleted by `cut_leading_zeros`.
    pure function add(arr1, arr2) result(ret)
        integer, contiguous, intent(in) :: arr1(:), arr2(:)
        integer, allocatable, dimension(:) :: ret, tmp1, tmp2

        associate (x => max(size(arr1), size(arr2)) + 1)
            allocate (tmp1(x), tmp2(x))
            tmp1 = 0; tmp1(x - size(arr1) + 1:x) = arr1(:)
            tmp2 = 0; tmp2(x - size(arr2) + 1:x) = arr2(:)
        end associate
        ret = cut_leading_zeros(carry(tmp1 + tmp2))
    end function add

    !> The core function to provide the ability of subtraction.
    !>
    !> Each element of the input array is an 32bit integer greater than or
    !> equal to zero and less than or equal to nine.
    pure function sub(arr1, arr2) result(ret)
        integer, contiguous, intent(in) :: arr1(:), arr2(:)
        integer, allocatable, dimension(:) :: ret, tmp1, tmp2, tmp
        integer :: i

        associate (x => max(size(arr1), size(arr2)) + 1)
            allocate (tmp1(x), tmp2(x), tmp(x))
            tmp1 = 0; tmp2 = 0; tmp = 0
            tmp1(x - size(arr1) + 1:x) = arr1(:)
            tmp2(x - size(arr2) + 1:x) = arr2(:)
        end associate

        do i = size(tmp1), 2, -1
            if (tmp1(i) >= tmp2(i)) then
                tmp(i) = tmp1(i) - tmp2(i)
            else
                tmp(i) = tmp1(i) - tmp2(i) + 10
                tmp1(i - 1) = tmp1(i - 1) - 1
            end if
        end do
        ret = cut_leading_zeros(tmp)
    end function sub

    !> The core function to provide the ability of comparison.
    !>
    !>### Result
    !> | return value | result |
    !> |:------------:|:------:|
    !> |`ret = 1`  |arr1 > arr2|
    !> |`ret = 0`  |arr1 = arr2|
    !> |`ret = -1` |arr1 < arr2|
    pure function compare(arr1, arr2) result(ret)
        integer, contiguous, intent(in) :: arr1(:), arr2(:)
        integer :: ret, i

        if (size(arr1) > size(arr2)) then
            ret = 1; return
        else if (size(arr1) < size(arr2)) then
            ret = -1; return
        end if

        do i = size(arr1), 1, -1
            if (arr1(i) > arr2(i)) then
                ret = 1; return
            else if (arr1(i) < arr2(i)) then
                ret = -1; return
            else
                ret = 0
            end if
        end do
        if (ret /= 0) error stop 'compare: invalid output.'
    end function compare

    !> The core function to provide the ability of multiplication.
    !>
    !> Each element of the input array is an 32bit integer greater than or
    !> equal to zero and less than or equal to nine.
    pure function mul(arr1, arr2) result(ret)
        integer, contiguous, intent(in) :: arr1(:), arr2(:)
        integer, allocatable :: tmp(:), tmp_row(:), ret(:)
        integer :: i

        associate (s1 => size(arr1), s2 => size(arr2))
            allocate (tmp_row(s1 + s2), tmp(s1 + s2 + 2))
            tmp = 0

            do i = s2, 1, -1
                tmp_row = 0; tmp_row(i + 1:i + s1) = arr1
                tmp = tmp + carry(tmp_row*arr2(i))
            end do
        end associate
        ret = cut_leading_zeros(carry(tmp))
    end function mul

    !> The core function to provide the ability of power.
    !>
    !> Each element of the input array is an 32bit integer greater than or
    !> equal to zero and less than or equal to nine. The function uses
    !> a very simple implementation of exponentiation
    !> ([wiki links](https://en.wikipedia.org/wiki/
    !>Exponentiation_by_squaring)). The function currently doesn't
    !> support negative powers.
    pure recursive function pow2(arr, n) result(ret)
        integer, contiguous, intent(in) :: arr(:)
        integer, intent(in) :: n
        integer, allocatable :: ret(:)

        if (n < 0) error stop 'pow2: n is nonnegative.'

        if (n == 0) then
            ret = [1]
        else if (n == 1) then
            ret = arr
        else if (mod(n, 2) == 0) then
            ret = pow2(mul(arr, arr), n/2)
        else if (mod(n, 2) /= 0) then
            ret = pow2(mul(arr, arr), (n - 1)/2)
            ret = mul(arr, ret)
        end if
    end function pow2

end module euler_multiprecision_util_m


text:
module corr_arrays

use corr_parms

real :: lst(ilg,jlg)
real :: tsfc(ci,cj)
real :: qsfc(ci,cj)
real :: psfc(ci,cj)
real :: pres(kz)
real :: theta(ci,cj,kz)
real :: spfh(ci,cj,kz)
real :: prsmd(ilg,jlg)
real :: tc15(ilg,jlg)
real :: tc55(ilg,jlg)
real :: radini(ilg,jlg)
real :: radfin(ilg,jlg)
real :: satang(ilg,jlg)
integer :: lookup_i(ilg,jlg)
integer :: lookup_j(ilg,jlg)

end module

max_stars_repo_name: bucricket/projectMASviirsMAC

max_issues_repo_name: bucricket/projectMASviirsMAC

max_forks_repo_name: bucricket/projectMASviirsMAC


! PR fortran/63938
! { dg-do run }

program pr63938_1
  integer :: i, x(1)
  x(1) = 0
!$omp parallel do
  do i = 1, 1000
    !$omp atomic
    x(1) = x(1) + 1
  end do
!$omp end parallel do
  if (x(1) .ne. 1000) call abort
end program pr63938_1


text:
      SUBROUTINE INVARS
C--
C--   SUBROUTINE INVARS (ISTART)
C--
C--   Purpose : initialize all spectral variables starting from
C--             either a reference atmosphere or a restart file
C--   Input :   ISTART = 0 : reference atmosphere (at rest)
C--                    = 1 : restart file
C--   Initialized common blocks : DATE1, DYNSP1, DYNSP2 (PHIS only),
C--                               SFCANOM, SFCFLUX
C--
      include "atparam.h"
      include "atparam1.h"

      include "com_date.h"
      include "com_tsteps.h"

      include "com_surfcon.h"

      include "com_dyncon0.h"
      include "com_dyncon1.h"

      include "com_dynvar.h"

      COMPLEX ZERO, CCON, SURFS(MX,NX)
      REAL  SURFG(IX,IL)

      GAM1 = GAMMA/(1000.*GRAV)
      ZERO = (0.,0.)
      CCON = (1.,0.)*SQRT(2.)

C--   1. Compute spectral surface geopotential

      CALL SPEC (PHI0,PHIS)
      IF (IX.EQ.IY*4) CALL TRUNCT (PHIS)

      CALL GRID (PHIS,PHIS0,1)

      IF (ISTART.EQ.0) THEN

C--   2. Start from reference atmosphere (at rest)

        print*, ' starting from rest'

        IYEAR  = IYEAR0
        IMONTH = IMONT0

C       2.1 Set vorticity, divergence and tracers to zero

        DO K=1,KX
         DO N=1,NX
          DO M=1,MX
            VOR(M,N,K,1)=ZERO
            DIV(M,N,K,1)=ZERO
          ENDDO
         ENDDO
        ENDDO

        DO ITR=1,NTR
         DO K=1,KX
          DO N=1,NX
           DO M=1,MX
             TR(M,N,K,1,ITR)=ZERO
           ENDDO
          ENDDO
         ENDDO
        ENDDO

C       2.2 Set reference temperature :
C           tropos:  T = 288 degK at z = 0, constant lapse rate
C           stratos: T = 216 degK, lapse rate = 0

        TREF  = 288.
        TTOP  = 216.
        GAM2  = GAM1/TREF
        RGAM  = RGAS*GAM1
        RGAMR = 1./RGAM

C       Surface and stratospheric air temperature

        DO N=1,NX
         DO M=1,MX
           T(M,N,1,1)=ZERO
	   T(M,N,2,1)=ZERO
           SURFS(M,N)=-GAM1*PHIS(M,N)
         ENDDO
        ENDDO

        T(1,1,1,1)=CCON*TTOP
        T(1,1,2,1)=CCON*TTOP
        SURFS(1,1)=CCON*TREF-GAM1*PHIS(1,1)

C       Temperature at tropospheric levels
        DO K=3,KX
          FACTK=FSG(K)**RGAM
          DO N=1,NX
           DO M=1,MX
             T(M,N,K,1)=SURFS(M,N)*FACTK
           ENDDO
          ENDDO
        ENDDO

C       2.3 Set log(ps) consistent with temperature profile
C           p_ref = 1013 hPa at z = 0

        RLOG0=LOG(1.013)
        DO J=1,IL
         DO I=1,IX
           SURFG(I,J)=RLOG0+RGAMR*LOG(1.-GAM2*PHIS0(I,J))
         ENDDO
        ENDDO

        CALL SPEC (SURFG,PS)
        IF (IX.EQ.IY*4) CALL TRUNCT (PS)

C       2.4 Set tropospheric spec. humidity in g/kg
C           Qref = RHref * Qsat(288K, 1013hPa)

        ESREF=17.
        QREF=REFRH1*0.622*ESREF
        QEXP=HSCALE/HSHUM

C       Spec. humidity at the surface

        DO J=1,IL
         DO I=1,IX
           SURFG(I,J)=QREF*EXP(QEXP*SURFG(I,J))
         ENDDO
         print *, ' Q0 jlat = ', j, surfg(1,j)
        ENDDO

        CALL SPEC (SURFG,SURFS)
        IF (IX.EQ.IY*4) CALL TRUNCT (SURFS)

C       Spec. humidity at tropospheric levels

        DO K=3,KX
          FACTK=FSG(K)**QEXP
          print *, 'vertical scale factor at level ', k, factk
          DO N=1,NX
           DO M=1,MX
             TR(M,N,K,1,1)=SURFS(M,N)*FACTK
           ENDDO
          ENDDO
        ENDDO

C       Print diagnostics from initial conditions

        CALL DIAGNS (1,0)

      ELSE

C--   3. Start from restart file

        print*,' reading a restart file'

        CALL RESTART (0)

C       Print diagnostics from initial conditions

        CALL DIAGNS (2,0)

      ENDIF

C--
      RETURN
      END


! C function declarations

type, bind(C) :: float3
    real(kind = 4) :: x, y, z
end type

interface
    function create_npcf_c(timesRans, numShells, volBox, rMin, rMax) bind(C, name="create_npcf")
        use iso_c_binding
        implicit none
        type(c_ptr) :: create_npcf_c
        integer(c_int), value :: timesRans
        integer(c_int), value :: numShells
        real(c_double), value :: volBox
        real(c_double), value :: rMax
        real(c_double), value :: rMin
    end function

    subroutine delete_npcf_c(obj) bind(C, name="delete_npcf")
        use iso_c_binding
        implicit none
        type(c_ptr), value :: obj
    end subroutine

    function get_shells_c(obj, shells) bind(C, name="get_shells")
        use iso_c_binding
        implicit none
        integer(c_int) :: get_shells_c
        type(c_ptr), value :: obj
        real(c_double), dimension(:) :: shells
    end function

    function get_num_triangles_c(obj) bind(C, name="get_num_triangles")
        use iso_c_binding
        implicit none
        integer(c_int) :: get_num_triangles_c
        type(c_ptr), value :: obj
    end function

    function get_triangles_c(obj, tris) bind(C, name="get_triangles")
        use iso_c_binding
        import :: float3
        integer(c_int) :: get_triangles_c
        type(c_ptr), value :: obj
        type(float3), dimension(:) :: tris
    end function

    function set_num_particles_c(obj, numParts) bind(C, name="set_num_particles")
        use iso_c_binding
        implicit none
        integer(c_int) :: set_num_particles_c
        type(c_ptr), value :: obj
        integer(c_int), value :: numParts
    end function

    function get_2pt_size_c(obj) bind(C, name="get_2pt_size")
        use iso_c_binding
        implicit none
        integer(c_int) :: get_2pt_size_c
        type(c_ptr), value :: obj
    end function

    function get_3pt_size_c(obj) bind(C, name="get_3pt_size")
        use iso_c_binding
        implicit none
        integer(c_int) :: get_3pt_size_c
        type(c_ptr), value :: obj
    end function

    function calculate_correlations_c(obj, galaxies) bind(C, name="calculate_correlations")
        use iso_c_binding
        import :: float3
        integer(c_int) :: calculate_correlations_c
        type(c_ptr), value :: obj
        type(float3), dimension(:) :: galaxies
    end function

    function calculate_2pt_c(obj, galaxies) bind(C, name="calculate_2pt")
        use iso_c_binding
        import :: float3
        integer(c_int) :: calculate_2pt_c
        type(c_ptr), value :: obj
        type(float3), dimension(:) :: galaxies
    end function

    function get_2pt_c(obj, twoPoint) bind(C, name="get_2pt")
        use iso_c_binding
        implicit none
        integer(c_int) :: get_2pt_c
        type(c_ptr), value :: obj
        real(c_double), dimension(:) :: twoPoint
    end function

    function get_3pt_c(obj, threePoint) bind(C, name="get_3pt")
        use iso_c_binding
        implicit none
        integer(c_int) :: get_3pt_c
        type(c_ptr), value :: obj
        real(c_double), dimension(:) :: threePoint
    end function
end interface


text:
! H0 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
! H0 X
! H0 X   libAtoms+QUIP: atomistic simulation library
! H0 X
! H0 X   Portions of this code were written by
! H0 X     Albert Bartok-Partay, Silvia Cereda, Gabor Csanyi, James Kermode,
! H0 X     Ivan Solt, Wojciech Szlachta, Csilla Varnai, Steven Winfield.
! H0 X
! H0 X   Copyright 2006-2010.
! H0 X
! H0 X   These portions of the source code are released under the GNU General
! H0 X   Public License, version 2, http://www.gnu.org/copyleft/gpl.html
! H0 X
! H0 X   If you would like to license the source code under different terms,
! H0 X   please contact Gabor Csanyi, gabor@csanyi.net
! H0 X
! H0 X   Portions of this code were written by Noam Bernstein as part of
! H0 X   his employment for the U.S. Government, and are not subject
! H0 X   to copyright in the USA.
! H0 X
! H0 X
! H0 X   When using this software, please cite the following reference:
! H0 X
! H0 X   http://www.libatoms.org
! H0 X
! H0 X  Additional contributions by
! H0 X    Alessio Comisso, Chiara Gattinoni, and Gianpietro Moras
! H0 X
! H0 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

! print 1 out of every n frames
#include "error.inc"
program decimate
use libatoms_module
implicit none
  type(Atoms) :: at
  integer :: i, n, error
  type(dictionary) :: cli_params
  character(len=STRING_LENGTH) :: infilename, outfilename
  type(cinoutput) :: infile_io, outfile_io

  call system_initialise(verbosity=PRINT_SILENT)
  call verbosity_push(PRINT_NORMAL)

  call initialise(cli_params)
  call param_register(cli_params,"n",PARAM_MANDATORY, n, help_string="Number of frames to skip between printing")
  call param_register(cli_params,"infile","stdin", infilename, help_string="input filename")
  call param_register(cli_params,"outfile","stdout", outfilename, help_string="output filename")
  if (.not. param_read_args(cli_params)) then
    call system_abort("Usage: decimate n=(1) [infile=(stdin)] [outfile=(stdout)]")
  endif
  call finalise(cli_params)

  call initialise(infile_io, infilename, INPUT)
  call initialise(outfile_io, outfilename, OUTPUT)

  call read(at, infile_io, error=error)
  HANDLE_ERROR(error)
  call write(at, outfile_io)
  i = 0
  do while (error == 0)
    i = i + n
    call read(at, infile_io, frame=i, error=error)
    if (error /= 0) then
       if (error == ERROR_IO_EOF) exit
       HANDLE_ERROR(error)
    endif
    call write(at, outfile_io)
  end do

  call finalise(outfile_io)
  call finalise(infile_io)

  call verbosity_pop()
  call system_finalise()
end program decimate

max_stars_repo_name: albapa/QUIP

max_issues_repo_name: albapa/QUIP

max_forks_repo_name: albapa/QUIP


text:
      program srtest

      include 'p4f.h'

      call p4init()
      call p4crpg()
      if (p4myid() .eq. 0) then
	  call fmaster()
      else
	  call fworker()
      endif
      call p4cleanup()
      print *,'mainline exiting normally'
      end

      subroutine fmaster()

      include 'p4f.h'

      integer*4 i,slaves,toid,type,from,retcde,recvlen,buflen,length
      character*40 buffer
      integer*4 TAGCNT, TAGDAT, TAGEND
      parameter (TAGCNT = 10)
      parameter (TAGDAT = 20)
      parameter (TAGEND = 30)

      print 11,'Entering fmaster'
11    format(a)
      slaves = p4ntotids() - 1
      length = 0
      buflen = 40

      do 10 i = 1,slaves
         call p4sendr(TAGCNT,i,buffer,length,retcde)
10    continue

20    print *,'Enter a string: '
      read (*,99,end=50) buffer
99    format(a40)

      do 30 length=40,1,-1
         if(buffer(length:length) .ne. ' ') goto 40
30    continue
      length = 0
40    continue

      toid = 1
      print *,'master sending msg of length ',length
      call p4send(TAGDAT,toid,buffer,length,retcde)
      buffer = ' '
      type = TAGDAT
      from = -1
      call p4recv(type,from,buffer,buflen,recvlen,retcde)

      print *,'MASTER receives from=',from,' buffer=',buffer
      length = 0
      goto 20
50    continue

      do 60 i = 1,slaves
         call p4sendr(TAGEND,i,buffer,buflen,retcde)
60    continue

      print *,'Master exiting normally'
      end

      subroutine fworker()

      include 'p4f.h'

      character*40 buffer
      integer*4 type, from, next, done, procid, length, buflen
      integer*4 numsl, retcde, recvlen
      integer*4 TAGCNT, TAGDAT, TAGEND
      parameter (TAGCNT = 10)
      parameter (TAGDAT = 20)
      parameter (TAGEND = 30)

      numsl = p4ntotids() - 1
      procid = p4myid()
      buflen = 40

      print 200,'slave ',procid,' has started'
 200  format(a,i2,a)
      call p4flush

      if (procid .eq. numsl) then
         next = 0
      else
         next = procid + 1
      endif

      print 201,'slave ',procid,' next = ',next
 201  format(a,i2,a,i2)
      call p4flush

      length = 40
      from = -1
      type = TAGCNT
      call p4recv(type,from,buffer,length,recvlen,retcde)
      call p4flush
      done = 0

50    if (done .ne. 0) goto 100

         buffer = ' '
         length = 40
         from = -1
         type = -1
         call p4recv(type,from,buffer,length,recvlen,retcde)
         call p4flush
         if (type .eq. TAGEND) then
            done = 1
         else
            call p4send(TAGDAT,next,buffer,recvlen,retcde)
         endif
         goto 50

100   continue

      end

max_stars_repo_name: RWTH-OS/MP-MPICH

max_issues_repo_name: RWTH-OS/MP-MPICH

max_forks_repo_name: RWTH-OS/MP-MPICH


text:
      SUBROUTINE TRAIL
C
C     MODULE TO INTERROGATE OR ALTER ANY VALUE OF A 6 WORD MATRIX
C     OR TABLE TRAILER
C
C     DMAP CALL
C
C        TRAILER  DB / /*OPT*/WORD/S,N,VALUE $
C
C     INPUT DATA BLOCKS
C
C        DB - DATA BLOCK FOR WHICH TRAILER IS TO BE ALTERED OR READ
C
C     PARAMETERS
C
C        OPT   - BCD,INPUT.
C                RETURN - VALUE OF SPECIFIED TRAILER WORD IS TO
C                         BE RETURNED
C                STORE  - VALUE OF SPECIFIED TRAILER WORD IS TO
C                         CHANGED
C        WORD  - INTEGER,INPUT. DESIRED WORD OF TRAILER
C        VALUE - INTEGER,INPUT OR OUTPUT. LOCATION WHERE VALUED WILL
C                RETURNED OR FROM WHICH REPLACEMENT VALUE WILL BE
C                TAKEN.  RETURNED NEGATIVE IF DB IS PURGED.
C
C     FOR MATRIX DATA BLOCKS, THE TRAILER POSITIONS ARE AS FOLLOWS
C
C        WORD 1 - NUMBER OF COLUMNS
C        WORD 2 - MUNBER OF ROWS
C        WORD 3 - MATRIX FORM
C        WORD 4 - TYPE OF ELEMENTS
C        WORD 5 - MAXIMUM NUMBER OF NON-ZERO WORDS IN ANY ONE COLUMN
C        WORD 6 - MATRIX DENSITY * 100
C
      EXTERNAL        LSHIFT  ,ANDF    ,ORF
      INTEGER         DB      ,OPT     ,WORD     ,VALUE    ,STORE(2) ,
     1                MCB(7)  ,FIAT    ,FIST     ,RETURN(2),MODNAM(2),
     2                ORF     ,ANDF
      CHARACTER       UFM*23
      COMMON /XMSSG / UFM
      COMMON /XFIAT / FIAT(2)
      COMMON /XFIST / FIST(2)
      COMMON /BLANK / OPT(2)  ,WORD    ,VALUE
      COMMON /SYSTEM/ IBUF    ,NOUT    ,DUM21(21),ICFIAT
      DATA    STORE / 4HSTOR  ,4HE    /
      DATA    RETURN/ 4HRETU  ,4HRN   /
      DATA    MODNAM/ 4HTRAI  ,4HLER  /
C
C     GET TRAILER
C
      DB = 101
      MCB(1) = DB
      CALL RDTRL (MCB)
      IF (MCB(1) .LE. 0) GO TO 70
C
C     TEST ILLEGAL PARAMETER VALUES AND BRANCH ON OPT
C
      IF (WORD.LT.1 .OR. WORD.GT.6) GO TO 100
      IF (OPT(1).EQ.RETURN(1) .AND. OPT(2).EQ.RETURN(2)) GO TO 10
      IF (OPT(1).EQ.STORE(1)  .AND. OPT(2).EQ.STORE(2) ) GO TO 20
      GO TO 300
C
C     RETURN OPTION
C
   10 VALUE = MCB(WORD+1)
      RETURN
C
C     STORE OPTION
C
C     SEARCH FIST FOR THE FILE
C
   20 N = FIST(2)*2 + 1
      DO 30 I = 3,N,2
      IF (FIST(I) .NE. DB) GO TO 30
      INDEX = FIST(I+1) + 1
      GO TO 40
   30 CONTINUE
      GO TO 70
C
C     PACK THE TRAILER INFORMATION INTO THE REQUESTED WORD.
C     MAKE SURE THE NUMBER IS POSITIVE AND .LE. 16 BITS IF ICFIAT=8
C
   40 IF (VALUE  .LT.  0) GO TO 200
      IF (ICFIAT .EQ. 11) GO TO 60
      IF (VALUE .GT. 65535) GO TO 200
      IW = (WORD+1)/2 + 2
      IF (WORD .EQ. (WORD/2*2)) GO TO 50
C
C     WORD IS ODD
C
      MASK = 65535
      FIAT(INDEX+IW) = ORF(ANDF(FIAT(INDEX+IW),MASK),LSHIFT(VALUE,16))
      RETURN
C
C     WORD IS EVEN
C
   50 MASK = LSHIFT(65535,16)
      FIAT(INDEX+IW) = ORF(ANDF(FIAT(INDEX+IW),MASK),VALUE)
      RETURN
C
C     ICFIAT = 11, TRAILER WORDS ARE NOT PACKED
C
   60 IW = 2
      IF (WORD .GE. 4) IW = 4
      FIAT(INDEX+IW+WORD) = VALUE
      RETURN
C
C     PURGED DATA BLOCK
C
   70 VALUE = -1
      RETURN
C
C     ERROR CONDITIONS
C
  100 WRITE  (NOUT,110) UFM,WORD
  110 FORMAT (A23,' 2202.  PARAMETER, WORD, HAS ILLEGAL VALUE OF',I9)
      GO TO 500
C
  200 WRITE  (6,210) UFM,VALUE
  210 FORMAT (A23,' 2202.  PARAMETER, VALUE, HAS ILLEGAL VALUE OF',I9)
      GO TO 500
C
  300 WRITE  (NOUT,310) UFM,OPT
  310 FORMAT (A23,' 2202.  PARAMETER, OPT, HAS ILLEGAL VALUE OF ',2A4)
C
  500 CALL MESAGE (-37,0,MODNAM)
      RETURN
      END