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

C $Id: mmschm.f,v 1.2 1998/07/14 18:19:26 gdsjaar Exp $
C $Log: mmschm.f,v $
C Revision 1.2  1998/07/14 18:19:26  gdsjaar
C Removed unused variables, cleaned up a little.
C
C Changed BLUE labels to GREEN to help visibility on black background
C (indirectly requested by a couple users)
C
C Revision 1.1.1.1  1990/11/30 11:12:21  gdsjaar
C FASTQ Version 2.0X
C
c Revision 1.1  90/11/30  11:12:20  gdsjaar
c Initial revision
c
C
CC* FILE: [.QMESH]MMSCHM.FOR
CC* MODIFIED BY: TED BLACKER
CC* MODIFICATION DATE: 7/6/90
CC* MODIFICATION: COMPLETED HEADER INFORMATION
C
      SUBROUTINE MMSCHM (NPER, KKK, LLL, NNN, ML, MS, NSPR, ISLIST,
     &   NINT, IFLINE, NLPS, ILLIST, LINKL, LINKS, MXNPER, MAXPRM, MAX3,
     &   MXND, X, Y, NID, NNPS, ANGLE, XN, YN, NUID, LXK, KXL, NXL, LXN,
     &   IAVAIL, NAVAIL, CCW, REAL, SCHSTR, M1, ERR)
C***********************************************************************
C
C  MMSCHM - "M" MESH SCHEME; CALCULATE A REGULAR RECTANGULAR MESH
C
C***********************************************************************
C
      DIMENSION ISLIST(NSPR), NINT(ML), IFLINE(MS), NLPS(MS)
      DIMENSION ILLIST(MS*3), LINKL(2, ML), LINKS(2, MS)
      DIMENSION X(MXNPER), Y(MXNPER), NID(MXNPER*MAXPRM), NNPS(MAX3)
      DIMENSION ANGLE(MXNPER), XN(MXND), YN(MXND), NUID(MXND)
      DIMENSION LXK(4, MXND), KXL(2, 3*MXND), NXL(2, 3*MXND)
      DIMENSION LXN(4, MXND)
C
      CHARACTER*72 SCHSTR
C
      LOGICAL CCW, ERR, NORM, REAL
C
C  CALCULATE THE BASE OF THE RECTANGLE FOR THE REGION
C
      CALL GETM1 (ML, MS, MAX3, NSPR, ISLIST, NINT, IFLINE, NLPS,
     &   ILLIST, LINKL, LINKS, X, Y, NID, NNPS, ANGLE, NPER, SCHSTR, M1,
     &   CCW, NORM, REAL, ERR)
      IF (NORM) THEN
         CALL MESAGE ('FORCED RECTANGLE PRIMITIVE PROCESSING USED')
      ELSE
         CALL MESAGE ('GENERAL RECTANGLE PRIMITIVE PROCESSING USED')
      END IF
      M2 = NPER/2 - M1
C
C  CALCULATE A REGUALR MAPPED "RECTANGULAR" MESH
C
      KKKOLD = KKK
      LLLOLD = LLL
      NNNOLD = NNN
      CALL RMESH (NPER, MXND, X, Y, NID, XN, YN, NUID, LXK, KXL, NXL,
     &   LXN, M1, M2, KKK, KKKOLD, NNN, NNNOLD, LLL, LLLOLD, IAVAIL,
     &   NAVAIL, ERR)
C
      RETURN
      END

ext: f

max_stars_repo_name: milljm/seacas

max_issues_repo_name: milljm/seacas

max_forks_repo_name: milljm/seacas


text:
! -*- mode: F90; mode: font-lock -*-
! ------------------------------------------------------------------------------
! $Id$
! ------------------------------------------------------------------------------
! Module read_support_spec
! ------------------------------------------------------------------------------
! Code area 11: Basis operations
! ------------------------------------------------------------------------------

!!****h* Conquest/read_support_spec *
!!  NAME
!!   read_support_spec
!!  PURPOSE
!!   Sole purpose of the module is to contain the sbrt read_support
!!  USES
!!
!!  AUTHOR
!!   Mike Gillan
!!  CREATION DATE
!!   22/6/02
!!  MODIFICATION HISTORY
!!   2006/06/12 08:10 dave
!!    Added per-species PAO coefficient information derived type
!!   2008/02/06 08:10 dave
!!    Changed for output to file not stdout
!!  SOURCE
module read_support_spec

  use datatypes
  use global_module, ONLY: io_lun

  implicit none

  save

  type acz_spec
     integer :: which_ac
     integer :: which_zeta
     real(double) :: coeff
  end type acz_spec

  type support_spec
     integer :: no_of_acz
     type(acz_spec), pointer, dimension(:) :: acz_info
  end type support_spec

  type(support_spec), allocatable, dimension(:,:) :: support_info

  character(len=80) :: support_spec_file
  logical :: flag_read_support_spec

!!***

contains

! -----------------------------------------------------------
! Subroutine read_support
! -----------------------------------------------------------

!!****f* read_support_spec/read_support *
!!
!!  NAME
!!   read_support
!!  USAGE
!!
!!  PURPOSE
!!   Reads all data specifying the initial support functions of
!!   all atomic species in terms of pseudo-atomic orbitals.
!!
!!  INPUTS
!!   inode: id of node running the executable
!!   ionode: id of input/output node
!!   n_species: number of atomic species in the simulated system
!!
!!  USES
!!
!!  AUTHOR
!!   Mike Gillan
!!  CREATION DATE
!!  22/6/02
!!  MODIFICATION HISTORY
!!   2008/09/01 08:21 dave
!!    Added io_ routines from input_module
!!  SOURCE
  subroutine read_support(inode,ionode,n_species)

    use datatypes
    use GenComms, ONLY: cq_abort, gcopy, my_barrier
    use global_module, ONLY : iprint_basis
    use species_module, ONLY: nsf_species, npao_species
    use maxima_module, ONLY: maxnsf
    use numbers, ONLY: one
    use input_module, ONLY: io_assign, io_close

    implicit none

    integer, intent(in) :: inode, ionode, n_species
    integer :: ios, lun, n_acz, n_sp, n_s, info, l, n_coeff, i, maxacz
    integer, allocatable, dimension(:) :: ibuff1, ibuff2 ! n_species_max*NSF*max_acz_in_spec
    real(double), allocatable, dimension(:) :: rbuff

    allocate(support_info(n_species,maxnsf),STAT=info)
    if(info/=0) call cq_abort("read_support_spec: support info: ",n_species)
    if(inode == ionode) then
       if(iprint_basis >= 3) then
          write(unit=io_lun,fmt='(//1x,60("+")/10x,"no of species for reading",&
               &" support spec:",i5/1x,60("+"))') n_species
       end if
       call io_assign(lun)
       if(iprint_basis >= 3) then
          write(unit=io_lun,fmt='(/" read_support_spec: input file unit no:",&
               &i3)') lun
       end if
       if(iprint_basis >= 3) then
          write(unit=io_lun,fmt='(/" read_support_spec: name of input file:",&
               &a80)') support_spec_file
       end if
       if(flag_read_support_spec) then
          open(unit=lun,file=support_spec_file,status='old',iostat=ios)
          if(ios /= 0) then
             call cq_abort('read_support_spec: input file not found')
          end if
          do n_sp = 1, n_species
             if(iprint_basis >= 1) then
                write(unit=io_lun,fmt='(/1x,60("*")/5x,"species no",i5/1x,60("*")/)') n_sp
                write(unit=io_lun,fmt='(/" no. of support functions:",i5)') nsf_species(n_sp)
             end if

             do n_s = 1, nsf_species(n_sp)
                if(iprint_basis >= 1) then
                   write(unit=io_lun,fmt='(/1x,40("-")/8x," support function no:",i5/&
                        &1x,40("-"))') n_s
                end if

                read(unit=lun,fmt=*) support_info(n_sp,n_s)%no_of_acz
                if(iprint_basis >= 1) then
                   write(unit=io_lun,fmt='(" no. of acz to specify support&
                        & function:",i5)') support_info(n_sp,n_s)%no_of_acz
                end if
                allocate(support_info(n_sp,n_s)%acz_info(support_info(n_sp,n_s)%no_of_acz))
                do n_acz = 1, support_info(n_sp,n_s)%no_of_acz
                   read(unit=lun,fmt=*) &
                        support_info(n_sp,n_s)%acz_info(n_acz)%which_ac, &
                        support_info(n_sp,n_s)%acz_info(n_acz)%which_zeta, &
                        support_info(n_sp,n_s)%acz_info(n_acz)%coeff
                   if(iprint_basis >= 1) then
                      write(unit=io_lun,fmt='(" acz no:",i5," which ac:",i5,&
                           &" which zeta:",i5," coeff:",e15.6)') n_acz, &
                           &support_info(n_sp,n_s)%acz_info(n_acz)%which_ac, &
                           &support_info(n_sp,n_s)%acz_info(n_acz)%which_zeta, &
                           &support_info(n_sp,n_s)%acz_info(n_acz)%coeff
                   end if
                end do ! n_acz
             end do ! n_s = nsf_species
          end do ! n_sp = n_species
          call io_close(lun)
       else
          do n_sp = 1, n_species
             if(iprint_basis >= 1) then
                write(unit=io_lun,fmt='(/1x,60("*")/5x,"species no",i5/1x,60("*")/)') n_sp
                write(unit=io_lun,fmt='(/" no. of support functions:",i5)') nsf_species(n_sp)
             end if
             if(npao_species(n_sp)/=nsf_species(n_sp)) then
                write(io_lun,fmt='(10x,"If there are more PAOs than SFs you must specify how they map using support.dat")')
                call cq_abort("PAO-to-blip map not specified.")
             end if
             do n_s = 1, nsf_species(n_sp)
                if(iprint_basis >= 1) then
                   write(unit=io_lun,fmt='(/1x,40("-")/8x," support function no:",i5/&
                        &1x,40("-"))') n_s
                end if

                support_info(n_sp,n_s)%no_of_acz = 1
                if(iprint_basis >= 1) then
                   write(unit=io_lun,fmt='(" no. of acz to specify support&
                        & function:",i5)') support_info(n_sp,n_s)%no_of_acz
                end if
                allocate(support_info(n_sp,n_s)%acz_info(support_info(n_sp,n_s)%no_of_acz))
                support_info(n_sp,n_s)%acz_info(1)%which_ac = n_s
                support_info(n_sp,n_s)%acz_info(1)%which_zeta = 1
                support_info(n_sp,n_s)%acz_info(1)%coeff = one
                if(iprint_basis >= 1) then
                   write(unit=io_lun,fmt='(" acz no:",i5," which ac:",i5,&
                        &" which zeta:",i5," coeff:",e15.6)') 1, &
                        &support_info(n_sp,n_s)%acz_info(1)%which_ac, &
                        &support_info(n_sp,n_s)%acz_info(1)%which_zeta, &
                        &support_info(n_sp,n_s)%acz_info(1)%coeff
                end if
             end do ! n_s = nsf_species
          end do ! n_sp = n_species
       end if
    end if ! inode==ionode

    ! distribute to all nodes all the data read above
    allocate(ibuff1(n_species*maxnsf))
    if(inode == ionode) then
       do n_sp = 1, n_species
          do n_s = 1, nsf_species(n_sp)
             ibuff1((n_sp-1)*maxnsf + n_s) = support_info(n_sp,n_s)%no_of_acz
          end do
       end do
    end if

    call gcopy(ibuff1,n_species*maxnsf)

    maxacz = 0
    do n_sp = 1, n_species
       do n_s = 1, nsf_species(n_sp)
          support_info(n_sp,n_s)%no_of_acz = ibuff1((n_sp-1)*maxnsf + n_s)
          if(support_info(n_sp,n_s)%no_of_acz>maxacz) maxacz = support_info(n_sp,n_s)%no_of_acz
          if(inode/=ionode) allocate(support_info(n_sp,n_s)%acz_info(support_info(n_sp,n_s)%no_of_acz))
       end do
    end do
    deallocate(ibuff1)
    if((inode == ionode).and.(iprint_basis >= 2)) then
       write(unit=io_lun,fmt='(/" read_support_spec:&
            & support_info(n_sp,n_s)%no_of_acz distributed")')
    end if

    allocate(ibuff1(n_species*maxnsf*maxacz),ibuff2(n_species*maxnsf*maxacz),rbuff(n_species*maxnsf*maxacz))
    if(inode == ionode) then
       do n_sp = 1, n_species
          do n_s = 1, nsf_species(n_sp)
             do n_acz = 1, support_info(n_sp,n_s)%no_of_acz
                ibuff1((n_sp-1)*maxnsf*maxacz +&
                     &(n_s-1)*maxacz + n_acz) = &
                     &support_info(n_sp,n_s)%acz_info(n_acz)%which_ac
                ibuff2((n_sp-1)*maxnsf*maxacz +&
                     &(n_s-1)*maxacz + n_acz) = &
                     &support_info(n_sp,n_s)%acz_info(n_acz)%which_zeta
                rbuff((n_sp-1)*maxnsf*maxacz +&
                     &(n_s-1)*maxacz + n_acz) = &
                     &support_info(n_sp,n_s)%acz_info(n_acz)%coeff
             end do
          end do
       end do
    end if

    call gcopy(ibuff1,n_species*maxnsf*maxacz)
    call gcopy(ibuff2,n_species*maxnsf*maxacz)
    call gcopy(rbuff,n_species*maxnsf*maxacz)

    do n_sp = 1, n_species
       do n_s = 1, nsf_species(n_sp)
          do n_acz = 1, support_info(n_sp,n_s)%no_of_acz
             support_info(n_sp,n_s)%acz_info(n_acz)%which_ac = &
                  &ibuff1((n_sp-1)*maxnsf*maxacz +&
                  &(n_s-1)*maxacz + n_acz)
             support_info(n_sp,n_s)%acz_info(n_acz)%which_zeta = &
                  &ibuff2((n_sp-1)*maxnsf*maxacz +&
                  &(n_s-1)*maxacz + n_acz)
             support_info(n_sp,n_s)%acz_info(n_acz)%coeff = &
                  &rbuff((n_sp-1)*maxnsf*maxacz +&
                  &(n_s-1)*maxacz + n_acz)
          end do
       end do
    end do
    deallocate(rbuff,ibuff2,ibuff1)
    if((inode == ionode).and.(iprint_basis >= 2)) then
       write(unit=io_lun,fmt='(/" read_support_spec: support_info(n_sp,n_s)%&
            &acz_info(n_acz)%which_ac distributed")')
       write(unit=io_lun,fmt='(/" read_support_spec: support_info(n_sp,n_s)%&
            &acz_info(n_acz)%which_zeta distributed")')
       write(unit=io_lun,fmt='(/" read_support_spec: support_info(n_sp,n_s)%&
            &acz_info(n_acz)%coeff distributed")')
    end if
    if(inode==ionode) then
       close(unit=lun)
    end if

    !do n_sp = n_species,1,-1
    !   do n_s = nsf_species(n_sp),1,-1
    !      deallocate(support_info(n_sp,n_s)%acz_info)
    !   end do
    !end do
    !deallocate(support_info,STAT=info)
    call my_barrier()

  end subroutine read_support
!!***

end module read_support_spec

max_stars_repo_name: negf/CONQUEST-release

max_issues_repo_name: negf/CONQUEST-release

max_forks_repo_name: negf/CONQUEST-release


text:
      MODULE h1elec_I
      INTERFACE
!...Generated by Pacific-Sierra Research 77to90  4.4G  10:47:20  03/09/06
      SUBROUTINE h1elec (NI, NJ, XI, XJ, SMAT)
      USE vast_kind_param,ONLY: DOUBLE
      INTEGER, INTENT(IN) :: NI
      INTEGER, INTENT(IN) :: NJ
      REAL(DOUBLE), DIMENSION(3), INTENT(IN) :: XI
      REAL(DOUBLE), DIMENSION(3), INTENT(IN) :: XJ
      REAL(DOUBLE), DIMENSION(9,9), INTENT(OUT) :: SMAT
      END SUBROUTINE
      END INTERFACE
      END MODULE

max_stars_repo_name: openmopac/MOPAC-archive

max_issues_repo_name: openmopac/MOPAC-archive

max_forks_repo_name: openmopac/MOPAC-archive


text:
      SUBROUTINE LAMX
C
C     LAMX MAKES OR EDITS THE LAMA DATA BLOCK
C
C     LAMX  EDIT,LAMA/LAMB/C,Y,NLAM=0 $
C     IF NLAM LT 0 MAKE LAMB A MATRIX OF 5 COLUMNS
C     LAMA  OMEGA FREQ GM GS
C     UNTIL GM = 0.0
C
C
      INTEGER SYSBUF,IST(10),TRL(7),BUFA,BUFB,BUFE,EDIT
C
      DIMENSION D(3),Z(7)
C
      COMMON /BLANK / NLAM
      COMMON /ZZZZZZ/ IZ(1)
      COMMON /SYSTEM/ SYSBUF,NOUT
      COMMON /CONDAS/ PI,TWOPI
      COMMON /UNPAKX/ ITO,II,IE,INCR
      COMMON /PACKX / ITYIN,ITYOUT,III,NNN,INCR1
      COMMON /OUTPUT/ HDG(96)
C
      EQUIVALENCE     (D(1),A),(D(2),B),(D(3),C)
      EQUIVALENCE     (Z(1),IZ(1))
C
      DATA    EDIT  , LAMA,LAMB /101,102,201/
      DATA    IST   / 21,6,7*0,7/
      DATA    LMA   / 1/,   IED /1/,  IZ2 /2/
      DATA    NAM   / 4HLAMX    /
C
C     INITILIZE AND DECIDE MODE OF OPERATIONS
C
      ICORE  = KORSZ(Z)
      TRL(1) = LAMA
      CALL RDTRL (TRL)
      IF (TRL(1) .LT. 0) LMA = 0
      TRL(1) = EDIT
      CALL RDTRL (TRL)
      IF (TRL(1) .LT. 0) IED = 0
      NCOL = TRL(2)
      IF (NCOL .EQ. 0) IED = 0
      IF (LMA.EQ.0 .AND. IED.EQ.0) GO TO 1000
      ITO = 1
      II  = 1
      INCR= 1
      IE  = TRL(3)
      IF (IE .GT. 3) IE = 3
      B = 0.0
      C = 0.0
      BUFB = ICORE - SYSBUF
      CALL GOPEN (LAMB,Z(BUFB),1)
      IF (LMA .EQ. 0) GO TO 200
      BUFA = BUFB - SYSBUF
      CALL GOPEN (LAMA,Z(BUFA),0)
      IF (NLAM .LT. 0) GO TO 500
      BUFE = BUFA - SYSBUF
      IF (IED .EQ. 0) GO TO 5
C
C      EDITING LAMA FROM EDIT
C
      CALL GOPEN (EDIT,Z(BUFE),0)
C
C     WRITE HEADER
C
    5 CALL READ (*10,*10,LAMA,Z,BUFE,1,NWR)
   10 CALL WRITE (LAMB,Z,NWR,1)
      IF (IED .EQ. 0) GO TO 100
C
C     MAKE RECORDS
C
      J = 0
      DO 50 I = 1,NCOL
      CALL READ (*60,*60,LAMA,Z,7,0,NWR)
      CALL UNPACK (*40,EDIT,D)
      IF (A.EQ.0.0 .AND. B.EQ.0.0 .AND. C.EQ.0.0) GO TO 40
      IF (C .LT. 0.0) GO TO 50
      Z(5) = Z(5)*(1.0+B) + A
      Z(4) = Z(5)*TWOPI
      Z(3) = Z(4)*Z(4)
      IF (C .NE. 0.0) Z(6) = C
      Z(7) = Z(6)*Z(3)
   40 J    = J+1
      IZ(1)= J
      IF (NLAM .LE. 0) GO TO 45
      IF (J .GT. NLAM) GO TO 180
   45 CALL WRITE (LAMB,Z,7,0)
   50 CONTINUE
   60 GO TO 180
C
C     COPY LAMA TO LAMB FOR NLAM RECORDS
C
  100 IF (NLAM .EQ. 0) GO TO 190
      J = NLAM
      M = 7*NLAM
      CALL READ (*180,*110,LAMA,Z,M,0,NWR)
      CALL WRITE (LAMB,Z,7*NLAM,0)
      GO TO 180
  110 CALL WRITE (LAMB,Z,NWR,0)
  180 TRL(1) = LAMB
      TRL(2) = J
      CALL WRTTRL (TRL)
  190 CALL CLOSE (LAMA,1)
      CALL CLOSE (LAMB,1)
      CALL CLOSE (EDIT,1)
      GO TO 1000
C
C      MAKE A NEW LAMB
C
  200 BUFE = BUFB - SYSBUF
      CALL GOPEN (EDIT,Z(BUFE),0)
      IF (NLAM .GT. 0) NCOL = MIN0(NCOL,NLAM)
C
C     WRITE HEADER
C
      CALL WRITE (LAMB,IST,50,0)
      CALL WRITE (LAMB,HDG,96,1)
C
C     MAKE RECORDS
C
      DO 220 I = 1,NCOL
      CALL UNPACK (*310,EDIT,D)
      GO TO 210
  310 D(1) = 0.0
      D(2) = 0.0
      D(3) = 0.0
  210 IZ(  1) = I
      IZ(IZ2) = I
      Z(5) = A
      Z(4) = TWOPI*A
      Z(3) = Z(4)*Z(4)
      Z(6) = C
      Z(7) = C*Z(3)
      CALL WRITE (LAMB,Z,7,0)
  220 CONTINUE
      J = NCOL
      GO TO 180
C
C     BUILD LAMB AS A MATRIX
C
  500 TRL(1) = LAMB
      TRL(2) = 0
      TRL(4) = 1
      TRL(5) = 1
      TRL(6) = 0
      TRL(7) = 0
      ITYIN  = 1
      ITYOUT = 1
      III    = 1
      INCR1  = 7
      CALL FWDREC (*190,LAMA)
      CALL READ (*190,*510,LAMA,Z,BUFA,0,NWR)
      CALL MESAGE (8,0,NAM)
      GO TO 190
  510 NLOOP = 0
      DO 520 I = 1,NWR,7
      IF (Z(I+5) .EQ.0.0) GO TO 530
      NLOOP = NLOOP +1
  520 CONTINUE
  530 IF (NLOOP .EQ. 0) GO TO 190
      TRL(3) = NLOOP
      NNN = NLOOP
      L   = 3
      DO 540 I = 1,5
      CALL PACK (Z(L),LAMB,TRL)
      L = L + 1
  540 CONTINUE
      CALL WRTTRL (TRL)
      GO TO 190
 1000 RETURN
      END


text:
!===============================================================================
!> Doubly linked list (DLL), carrying anything, as simply as possible
!===============================================================================
MODULE dll_mod
  implicit none
  private
  type, public:: dll_node_t
    class(dll_node_t), pointer:: prev => null()
    class(dll_node_t), pointer:: next => null()
    class(*), pointer:: car
    integer:: id=0
  end type
  type, public:: dll_t
    class(dll_node_t), pointer:: head => null()
    class(dll_node_t), pointer:: tail => null()
    integer:: n=0
    character(len=32):: name
  contains
    procedure:: init
    procedure:: append
    procedure:: prepend
    procedure:: find
    procedure:: insert_before
    procedure:: insert_after
    procedure:: remove
    procedure:: delete
    procedure:: pop
    procedure:: test
  end type
CONTAINS

!===============================================================================
!> Initialize a doubly-linked-list (DLL), with anonymous carry
!===============================================================================
SUBROUTINE init (self, name)
  class(dll_t):: self
  character(len=*), optional:: name
  !.............................................................................
  nullify (self%head, self%tail)
  !$omp atomic write
  self%n = 0
  if (present(name)) then
    self%name = name
  else
    self%name = 'dll'
  end if
END SUBROUTINE init

!===============================================================================
!> Add an item at the tail
!===============================================================================
SUBROUTINE append (self, new)
  class(dll_t):: self
  class(dll_node_t), pointer:: new
  !.............................................................................
  new%prev => self%tail
  nullify(new%next)
  if (.not.associated(self%head)) self%head => new
  if (associated(self%tail)) self%tail%next => new
  self%tail => new
  !$omp atomic
  self%n = self%n+1
END SUBROUTINE append

!===============================================================================
!> Add an item at the head
!===============================================================================
SUBROUTINE prepend (self, new)
  class(dll_t):: self
  class(dll_node_t), pointer:: new
  !.............................................................................
  if (.not.associated(self%tail)) self%tail => new
  new%next => self%head
  nullify (new%prev)
  self%head => new
  !$omp atomic
  self%n = self%n+1
END SUBROUTINE prepend

!===============================================================================
!> Find an item, by matching with a known member
!===============================================================================
SUBROUTINE find (self, old, prev, next)
  class(dll_t):: self
  class(dll_node_t), pointer:: prev, next
  class(*), pointer:: old
  !.............................................................................
  next => self%head
  nullify (prev)
  do while (associated(next))
    if (associated(next%car, old)) then
      return
    end if
    prev => next
    next => next%next
  end do
  nullify (prev, next)
END SUBROUTINE find

!===============================================================================
!> Insert a new item before an old one
!===============================================================================
SUBROUTINE insert_before (self, old, new)
  class(dll_t):: self
  class(dll_node_t), pointer:: old, new
  !.............................................................................
  new%prev => old%prev
  new%next => old
  if (associated(new%prev)) new%prev%next => new
  if (associated(new%next)) new%next%prev => new
  !$omp atomic
  self%n = self%n+1
END SUBROUTINE insert_before

!===============================================================================
!> Insert a new item after an old one
!===============================================================================
SUBROUTINE insert_after (self, old, new)
  class(dll_t):: self
  class(dll_node_t), pointer:: old, new
  !.............................................................................
  new%next => old%next
  old%next => new
  if (associated(new%prev)) new%prev%next => new
  if (associated(new%next)) new%next%prev => new
  !$omp atomic
  self%n = self%n+1
END SUBROUTINE insert_after

!===============================================================================
!> Remove an item from the DLL; note that the item is not deallocated!
!===============================================================================
SUBROUTINE remove (self, old)
  class(dll_t):: self
  class(dll_node_t), pointer:: old
  !.............................................................................
  if (associated(old%prev)) old%prev%next => old%next
  if (associated(old%next)) old%next%prev => old%prev
  if (associated(old, self%head)) then
    self%head => old%next
  end if
  if (associated(old, self%tail)) then
    self%tail => old%prev
  end if
  !$omp atomic
  self%n = self%n-1
END SUBROUTINE remove

!===============================================================================
!> Pop the last appended item off the list.  The list node and its carry should
!> be deallocated after use.
!===============================================================================
FUNCTION pop (self)
  class(dll_t):: self
  class(dll_node_t), pointer:: pop
  !.............................................................................
  pop => self%tail
  if (associated(pop%prev)) then
    self%tail => pop%prev
    nullify(self%tail%next)
  end if
  !$omp atomic
  self%n = self%n-1
END FUNCTION pop

!===============================================================================
!> Run a self test
!===============================================================================
SUBROUTINE test (self)
  class(dll_t):: self
  !.............................................................................
  integer:: n
  class(*), pointer:: car, save
  integer, pointer:: ip
  real, pointer:: rp
  class(dll_node_t), pointer:: prev, next
  integer:: i
  !-----------------------------------------------------------------------------
  ! Append five numbers
  !-----------------------------------------------------------------------------
  print *,'===================== Doubly linked list test ======================='
  print *, 'fill:', n
  do i=1,n
    if (i == 2) then
      allocate (next, ip)
      ip = i
      next%car => ip
    else
      allocate (next, rp)
      rp = i
      next%car => rp
    end if
    if (i == 3) save => next%car
    call self%append (next)
  end do
  !-----------------------------------------------------------------------------
  ! List them
  !-----------------------------------------------------------------------------
  print *, 'list:', self%n
  next => self%head
  do while (associated(next))
    call print (next%car)
    next => next%next
  end do
  !-----------------------------------------------------------------------------
  ! Find number 3
  !-----------------------------------------------------------------------------
  print *, 'find:'
  call print (save)
  call self%find (save, prev, next)
  call print (next%car)
  print *, 'previous:'
  call print (prev%car)
contains
  subroutine print (car)
    class(*), pointer:: car
    select type (car)
    type is (integer)
      print *, car
    type is (real)
      print *, car
    end select
  end subroutine print
END SUBROUTINE test

!===============================================================================
!> Delete a double-linked list, including carried content
!===============================================================================
SUBROUTINE delete (self)
  class(dll_t):: self
  !.............................................................................
  class (dll_node_t), pointer:: item, next
  !-----------------------------------------------------------------------------
  item => self%head
  do while (associated(item))
    next => item%next
    if (associated(item%car)) deallocate (item%car)
    deallocate (item)
    item => next
  end do
  nullify (self%head)
  nullify (self%tail)
  self%n = 0
  call self%init
END SUBROUTINE delete

END MODULE dll_mod

max_stars_repo_name: applejwjcat/dispatch

max_issues_repo_name: applejwjcat/dispatch

max_forks_repo_name: applejwjcat/dispatch


text:
      SUBROUTINE BSYNBplatt(NALLIN)
*
*-----------------------------------------------------------------------
*
* BSYNB is merely a continuation of BSYN, data transfered
* via scratch-file 14
*
* B takes line-absorption coefficients generated by A and
* calculates the spectra.
*
* Export version  1988-03-24  ********* Olof Morell *** Uppsala
*
*-----------------------------------------------------------------------
*
      include 'spectrum.inc'
*
      INTEGER    MAXNL
      PARAMETER (MAXNL=200)
      INTEGER    NLBL
      PARAMETER (NLBL=100)
*
      CHARACTER*20  LELE(100)
      character*50 mcode
      DIMENSION ION(100),SRXLAL(2*MAXNL),
     &          GFELOG(100),ABUND(100),ETAK(2*MAXNL),ETAR(ndp*2*MAXNL),
     &          ETAD(ndp),PROF(100),XC(ndp),CHIE(100)
      real fluxme(nlbl),icenter(nlbl)
      real mum,prf,iprf
      doubleprecision XL1,XL2,eqwidth
      doubleprecision XLBEG,XLEND,DEL,XLB(100),XL(100),XLAL(2*MAXNL)
      COMMON/ATMOS/ T(NDP),PE(NDP),PG(NDP),XI(NDP),MUM(NDP),RO(NDP),
     &  nnNTAU
*
      COMMON /CTRAN/X(NDP),S(NDP),BPLAN(NDP),XJ(NDP),HFLUX(NDP),XK(NDP)
     & ,FJ(NDP),SOURCE(NDP),TAUS(NDP),DTAUS(NDP),JTAU0,JTAU1,ISCAT
      COMMON/CSURF/ HSURF,Y1(NRAYS)
      COMMON/CANGLE/ NMY,XMY(nrays),XMY2(nrays),WMY(nrays)
      COMMON/TAUC/ TAU(ndp),DTAULN(ndp),NTAU
      COMMON/PIECES/ XL1,XL2,DEL,EPS,NMX,NLBLDU,IINT,XMYC,IWEAK
*
* extension for large number of wavelengths and lines (monster II)
      doubleprecision xlambda
      common/large/ xlambda(lpoint),maxlam,ABSO(NDP,lpoint),
     & absos(ndp,lpoint),absocont(ndp,lpoint)
*
      dimension fcfc(lpoint),y1cy1c(lpoint),
     & xlm(lpoint),jlcont(lpoint)

      logical findtau1,hydrovelo
      real velocity
      common/velo/velocity(ndp),hydrovelo
      logical debug
      data debug /.false./
*
* NLBLDU is a dummy, real val. of NLBL is set as parameter
*
      DATA eqwidth/0./,profold/0./

      PI=3.141593
*
* Initiate angle quadrature points and number of wavelengths
* per block (NLBL)
*
      NMY=NMX
      CALL GAUSI(NMY,0.,1.,WMY,XMY)
      DO 1 I=1,NMY
        XMY2(I)=XMY(I)*XMY(I)
    1 CONTINUE
*
* Initiate mode of calculation
* IINT =1  : intensity at MY=XMYC
* IINT =0  : flux
* XL1      : wavelength where syntetic spectrum starts
* XL2      : wavelength where syntetic spectrum stops
* DEL      : wavelength step
* IWEAK =1 : weak line approximation for L/KAPPA le. EPS
* note  XL2.gt.XL1
*
      iweak=0
      IF(IINT.GT.0) THEN
        NMY=NMY+1
        XMY(NMY)=XMYC
        wmy(nmy)=0.
        WRITE(7,200) XMYC,XL1,XL2,DEL
      END IF
      IF(IINT.EQ.0) WRITE(7,201) XL1,XL2,DEL
      IF(IWEAK.GT.0) WRITE(7,202) EPS
*
* Read model atmosphere
*
      READ(14,rec=1) MCODE,nlcont,xlm(1),BPLAN,XC,S,XI
cc       print*,'bsynbplatt read 14 ',mcode,nlcont,xlm(1),xi
      WRITE(7,203) MCODE(1:lenstr(mcode))
*
*
      if (iint.gt.0) then
* header for intensity file
*
        write(66) mcode
        write(66) 1.,1.
        write(66) maxlam
      endif

*
* Continuum calculations:
*
      do 1963 jc=1,nlcont
        READ(14,rec=jc) MCODE,idum,xlm(jc),BPLAN,XC,S,XI
        if (debug) then
          print*,'bsynbplatt read 14 ',mcode,jc,idum,xlm(jc)
          print*,' XC     S     BPlan'
          do k=1,ntau
            print*,xc(k),s(k),bplan(k)
          enddo
        endif
      DO 9 K=1,NTAU
        X(K)=XC(K)
    9 CONTINUE

      call traneqplatt(0)
      Y1CY1C(jc)=Y1(NMY)
      FCFC(jc)=4.*HSURF*pi
      IF(IINT.LE.0) WRITE(7,204) fcFC(jc),xlm(jc)
      IF(IINT.GT.0) WRITE(7,205) Y1Cy1c(jc),xlm(jc)
      if (debug) then
        print*,' continuum;  lambda = ',jc,xlm(jc),fcfc(jc)
      endif
1963  continue
*
* Start loop over wavelength blocks, each containing NLBL
* wavelength points
*
        numb=0
        DO 39 jjj=0,maxlam,nlbl
         xlend=xlambda(jjj+nlbl)
         do 390 j=1,nlbl
          xl(j)=xlambda(j+jjj)
          if (xl(j).le.xl2) then
          IF(IWEAK.LE.0.OR.IINT.LE.0) THEN
            DO 30 K=1,NTAU
* the continuum opacity is already included in abso
ccc              X(K)=XC(K)+ABSO(K,J+jjj)
              X(K)=ABSO(K,J+jjj)
              s(k)=absos(k,j+jjj)
              xlsingle=xl(j)
              BPLAN(k)=BPL(T(k),xlsingle)
   30       CONTINUE
c            if (debug) then
c              print*,' line;  lambda = ',xl(j)
c            endif
cc            if (abs(xl(j)-5240.41d0).lt.1.d-4.or.
cc     &          abs(xl(j)-5242.49d0).lt.1.d-4.or.
cc     &          abs(xl(j)-5241.62d0).lt.1.d-4)  then
cc              print*,'lambda = ',xl(j),' calling traneqplatt',idebug
cc              idebug=1
cc            else
              idebug=0
cc            endif
            call traneqplatt(idebug)
* interpolate continuum flux
            do 3691 jc=2,nlcont
              if(xlm(jc)-xlsingle.gt.0.) then
                jjc=jc-1
                goto 3692
              endif
3691        continue
3692        continue
	    if (nlcont.gt.1) then
              jjc=min(jjc,nlcont-1)
              fc=(fcfc(jjc+1)-fcfc(jjc))/
     &           (xlm(jjc+1)-xlm(jjc))*(xlsingle-xlm(jjc)) + fcfc(jjc)
              y1c=(y1cy1c(jjc+1)-y1cy1c(jjc))/
     &           (xlm(jjc+1)-xlm(jjc))*(xlsingle-xlm(jjc)) + y1cy1c(jjc)
              if (debug) then
                print*,'interpolating cont: ',jjc,xlm(jjc),xlsingle,
     &                                        xlm(jjc+1)
                print*,'fluxes ' ,fcfc(jjc),fc,fcfc(jjc+1)
              endif
	    else
              fc=fcfc(1)
              y1c=y1cy1c(1)
            endif
            PRF=4.*pi*HSURF/FC
* starting with version 12.1, flux is not divided by pi anymore.
* F_lambda integrated over lambda is sigma.Teff^4
*            fluxme(j)=hsurf*4.
            fluxme(j)=hsurf*4.*pi
            IF(IINT.GT.0) then
              iPRF=Y1(NMY)/Y1C
              icenter(j)=y1(nmy)
            endif
            PROF(J)=1.-PRF
          ELSE
            DO 31 K=1,NTAU
              ETAD(K)=ABSO(K,J+jjj)
              IF(ETAD(K).GT.EPS) GOTO 32
   31       CONTINUE
            CALL TRANW(NTAU,TAU,XMYC,BPLAN,XC,ETAD,DELI)
            PROF(J)=DELI/Y1C
            GOTO 39
   32       DO 33 K=1,NTAU
              X(K)=ABSO(K,J+jjj)
              s(k)=absos(k,j+jjj)
              xlsingle=xl(j)
              BPLAN(k)=BPL(T(k),xlsingle)
   33       CONTINUE
c            if (debug) then
c              print*,' line;  lambda = ',xl(j)
c            endif
             print*,'debug ',xlsingle
            call traneqplatt(0)
* interpolate continuum flux
            do 3693 jc=2,nlcont
              if(xlm(jc)-xlsingle.gt.0.) then
                jjc=jc-1
                goto 3694
              endif
3693        continue
3694        continue
	    if (nlcont.gt.1) then
              jjc=min(jjc,nlcont-1)
              fc=(fcfc(jjc+1)-fcfc(jjc))/
     &           (xlm(jjc+1)-xlm(jjc))*(xlsingle-xlm(jjc)) + fcfc(jjc)
              y1c=(y1cy1c(jjc+1)-y1cy1c(jjc))/
     &           (xlm(jjc+1)-xlm(jjc))*(xlsingle-xlm(jjc)) + y1cy1c(jjc)
            else
              fc=fcfc(1)
              y1c=y1cy1c(1)
	    endif
            PRF=4.*pi*HSURF/FC
* starting with version 12.1, flux is not divided by pi anymore.
* F_lambda integrated over lambda is sigma.Teff^4
*            fluxme(j)=hsurf*4.
            fluxme(j)=hsurf*4.*pi
            IF(IINT.GT.0) then
              iPRF=Y1(NMY)/Y1C
              icenter(j)=y1(nmy)
            endif
            PROF(J)=1.-PRF
          END IF
ccc          eqwidth=eqwidth+(prof(j)+profold)*del/2.
ccc          profold=prof(j)
*
* find depth where tau_lambda=1
          if (hydrovelo) then
            findtau1=.true.
          else
            findtau1=.false.
          endif
          if (findtau1) then
          taulambda=tau(1)*(x(1)+s(1))
          do k=2,ntau
            taulambda=taulambda+(tau(k)-tau(k-1))*(x(k)+s(k)+
     &                   x(k-1)+s(k-1))*0.5
            if (taulambda.ge.1.) then
               print333,xl(j),k,tau(k),taulambda,t(k),ro(k)
333            format(f10.3,x,i3,x,1pe11.4,x,1pe11.4,x,0pf7.1,x,
     &                1pe11.4,0p)

               goto 1966
            endif
          enddo
1966      continue
          endif
          endif
*
390      continue
cc         PRINT 2341,NUMB
ccc        WRITE(7,206) PROF
         IF(IINT .LE. 0) DUMMY=FC
         IF(IINT .GT. 0) DUMMY=Y1C

*
* Write spectrum on file for convolution with instrument profile
*
         do iplez=1,nlbl
           if (xl(iplez).le.xl2) then
             plezflux=1.-prof(iplez)

             if (iint.eq.0) then
* fluxme is sigma teff^4
               write(46,1964) xl(iplez),plezflux,fluxme(iplez)
1964           format(f11.3,1x,f10.5,1x,1pe12.5)

C OUTPUT IN CASE OF intensity calculation : STORE LIMB FUNCTION
             else
* We add intensity at center of disk in spectrum output.
               write(46,1965) xl(iplez),plezflux,fluxme(iplez),
     &                        icenter(iplez)
1965           format(f11.3,1x,f10.5,2(1x,1pe12.5))
               write(66) xl(iplez),nmy,fluxme(iplez),
     &                (xmy(nlx),y1(nlx),nlx=1,nmy)
             endif

           endif
         enddo
*
* Check whether there are more wavelength blocks
*
         NUMB=NUMB+1

39      CONTINUE
*
* End loop over wavelength blocks
*
ccc        WRITE(7,207) NALLIN
*
*
ccc      DO 442 LIN=1,NALLIN
ccc        WRITE(7,208) LELE(LIN),ION(LIN),XLB(LIN),CHIE(LIN),GFELOG(LIN),
ccc     &               ABUND(LIN)
ccc  442 CONTINUE
ccc      print 1111,eqwidth*1000.
1111    format(' eqwidth: ',f9.1,' mA')
        call clock
      RETURN
*
  100 FORMAT(4X,I1,6X,I3)
  207 FORMAT(' SPECTRUM CALCULATED FOR THE FOLLOWING ',I3,' LINES'
     &      /' ELEMENT LAMBDA XI(EV) LOG(GF) LOG(ABUND) LINE NO')
  208 FORMAT('   ',A2,I2,F9.2,F5.2,1X,F6.2,3X,F7.2,4X,I5)
  200 FORMAT(' INTENSITY SPECTRUM AT MU=',F6.2,' BETWEEN',F10.3,
     &      ' AND',F10.3,' A WITH STEP',F6.3,' A')
  201 FORMAT(' FLUX SPECTRUM BETWEEN',F10.3,' A AND ',F10.3,' A WITH',
     &       ' STEP',F6.3,' A')
  202 FORMAT(' ***WEAK LINE APPROX FOR L/KAPPA L.T.',F7.4)
  203 FORMAT(' MODEL ATMOSPHERE:',A)
  204 FORMAT(' CONTINUUM FLUX=',E12.4, ' AT ',f10.2,'AA')
  205 FORMAT(' CONTINUUM INTENSITY=',E12.4,' at ',f10.2,'AA')
  206 FORMAT(1X,10F8.4)
 2341 FORMAT(1X,'Block number: ',I8)
*
      END

max_stars_repo_name: dnidever/apogee

max_issues_repo_name: dnidever/apogee

max_forks_repo_name: dnidever/apogee


text:
module test2008_55_temp
!   type real_num
!   real :: x
!   end type

!   interface operator (.add.)
!   module procedure real_add
!   end interface

   interface operator (.not.)
   module procedure real_not
   end interface

!   contains
!   function real_add(a,b)
!   type(real_num) :: real_add
!   type(real_num), intent(in) :: a,b
!   real_add%x = a%x+b%x
!   end function real_add

end module

program main
 ! This is a F03 feature, allows to rename operators that are not intrinsic operators.
 ! use test2008_55_temp , operator(.plus.) => operator(.add.)
 ! use test2008_55_temp , operator(.plus.) => operator(.mynot.)

 ! type(real_num) :: a,b,c

 ! This is a F03 feature
 ! c = a .plus. b
 ! c = a .gt. b

 ! d = .mynot. b
end program

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
      SUBROUTINE STABL4(ITERM)
*
*
*       Stability test of four-body system.
*       -----------------------------------
*
      IMPLICIT REAL*8  (A-H,M,O-Z)
      LOGICAL  SWITCH,GTYPE,GTYPE0,AUX
      COMMON/CREG/  M(4),X(12),XD(12),P(12),Q(12),TIME4,ENERGY,EPSR2,
     &              XR(9),W(9),R(6),TA(6),MIJ(6),CM(10),RMAX4,TMAX,
     &              DS,TSTEP,EPS,NSTEP4,NAME4(4),KZ15,KZ27,NREG,NFN
      COMMON/TPR/   SWITCH,GTYPE,GTYPE0
      COMMON/KSAVE/  K1,K2
      REAL*8  RC(3),VC(3),RC0(3),VC0(3)
*
*
*       Transform to physical variables (retain value of SWITCH).
      AUX = SWITCH
      SWITCH = .FALSE.
      CALL ENDREG
      SWITCH = AUX
*
*       Specify indices of two least dominant bodies (denoted K3 & K4).
      K3 = 0
      DO 1 L = 1,4
          IF (L.EQ.K1.OR.L.EQ.K2) GO TO 1
          IF (K3.EQ.0) THEN
              K3 = L
          ELSE
              K4 = L
          END IF
    1 CONTINUE
*
*       Initialize scalars for orbital elements.
      RB0 = 0.0D0
      RB = 0.0D0
      RB1 = 0.0D0
      RB2 = 0.0D0
      RDOT = 0.0D0
      RDOT1 = 0.0D0
      RDOT2 = 0.0D0
      VREL2 = 0.0D0
      VREL20 = 0.0D0
      VREL21 = 0.0D0
      VREL22 = 0.0D0
*
*       Define binary masses of smallest & widest pair (K1 & K2 and K3 & K4).
      MB0 = M(K1) + M(K2)
      MB = M(K3) + M(K4)
*
*       Form separations & velocities of MB0, MB and their relative orbit.
      DO 10 K = 1,3
          J1 = 3*(K1-1) + K
          J2 = 3*(K2-1) + K
          J3 = 3*(K3-1) + K
          J4 = 3*(K4-1) + K
          RC0(K) = (M(K1)*X(J1) + M(K2)*X(J2))/MB0
          RC(K) = (M(K3)*X(J3) + M(K4)*X(J4))/MB
          VC0(K) = (M(K1)*XD(J1) + M(K2)*XD(J2))/MB0
          VC(K) = (M(K3)*XD(J3) + M(K4)*XD(J4))/MB
          RB0 = RB0 + (X(J1) - X(J2))**2
          RB = RB + (X(J3) - X(J4))**2
          RB1 = RB1 + (RC(K) - RC0(K))**2
          RB2 = RB2 + (RC0(K) - X(J3))**2
          RDOT = RDOT + (X(J3) - X(J4))*(XD(J3) - XD(J4))
          RDOT1 = RDOT1 + (RC(K) - RC0(K))*(VC(K) - VC0(K))
          RDOT2 = RDOT2 + (RC0(K) - X(J3))*(VC0(K) - XD(J3))
          VREL2 = VREL2 + (XD(J3) - XD(J4))**2
          VREL20 = VREL20 + (XD(J1) - XD(J2))**2
          VREL21 = VREL21 + (VC(K) - VC0(K))**2
          VREL22 = VREL22 + (VC0(K) - XD(J3))**2
   10 CONTINUE
*
*       Determine semi-major axis of inner binary.
      RB0 = SQRT(RB0)
      SEMI0 = 2.0/RB0 - VREL20/MB0
      SEMI0 = 1.0/SEMI0
      E0 = 1.0 - RB0/SEMI0
*
*       Form semi-major axis & eccentricity of outer pair.
      RB = SQRT(RB)
      SEMI = 2.0D0/RB - VREL2/MB
      SEMI = 1.0/SEMI
*     E = SQRT((1.0D0 - RB/SEMI)**2 + RDOT**2/(SEMI*MB))
*
*       Evaluate orbital elements of relative c.m. motion.
      RB1 = SQRT(RB1)
      SEMI1 = 2.0D0/RB1 - VREL21/CM(7)
      SEMI1 = 1.0/SEMI1
      E1 = SQRT((1.0D0 - RB1/SEMI1)**2 + RDOT1**2/(SEMI1*CM(7)))
*
*       Consider the inner triple.
      MB2 = MB0 + M(K3)
      RB2 = SQRT(RB2)
      SEMI2 = 2.0D0/RB2 - VREL22/MB2
      SEMI2 = 1.0/SEMI2
      E2 = SQRT((1.0D0 - RB2/SEMI2)**2 + RDOT2**2/(SEMI2*MB2))
*
*       Obtain standard stability ratio (outer pericentre / inner apocentre).
      RATIO = SEMI1*(1.0D0 - E1)/(SEMI0*(1.0D0 + E0))
*
*       Form coefficients for stability test (Valtonen, Vistas Ast 32, 1988).
*     AM = (2.65 + E0)*(1.0 + MB0/MB)**0.3333
*     FM = (2.0*MB0 - MB)/(3.0*MB)
*
*       Expand natural logarithm for small arguments.
*     IF (ABS(FM).LT.0.67) THEN
*         BM = FM*(1.0 - (0.5 - 0.3333*FM)*FM)
*     ELSE
*         BM = LOG(1.0D0 + FM)
*     END IF
*
*       Adopt mass dependent criterion of Harrington (A.J. 80) & Bailyn.
*     PCRIT = AM*(1.0 + 0.7*BM)*SEMI0
*
*       Form hierarchical stability ratio (Kiseleva & Eggleton 1995).
*     Q0 = MB/MB0
*     Q1 = MAX(M(K2)/M(K1),M(K1)/M(K2))
*     Q3 = Q0**0.33333
*     Q13 = Q1**0.33333
*     AR = 1.0 + 3.7/Q3 - 2.2/(1.0 + Q3) + 1.4/Q13*(Q3 - 1.0)/(Q3 + 1.0)
*     PCRIT = AR*SEMI0*(1.0D0 + E0)
*
*       Check stability (AM 1997; inner triple or well separated quadruple).
      ITERM = 0
      IF (RB1.GT.5.0*RB2.AND.E2.LT.1.0) THEN
          Q1 = M(K3)/MB0
          XFAC = (1.0 + Q1)*(1.0 + E2)/SQRT(1.0 - E2)
          PCRIT = 2.8*XFAC**0.4*SEMI0
          PMIN = SEMI2*(1.0 - E2)
          IF (PCRIT.LT.PMIN) THEN
              ITERM = -1
              RATIO = SEMI2*(1.0D0 - E2)/(SEMI0*(1.0D0 + E0))
              WRITE (6,15)  SEMI0, SEMI2, E0, E2, RATIO, RB0, RB2,
     &                      PCRIT, PMIN
   15         FORMAT ('  STABT:    A0 A2 E0 E2 RATIO R0 R2 PCR PM ',
     &                             1P,2E10.2,0P,2F7.3,F6.2,1P,4E9.1)
          END IF
      ELSE IF (RB1.GT.5.0*MAX(RB0,RB).AND.E1.LT.1.0.AND.
     &         MIN(SEMI0,SEMI).GT.0.0) THEN
*       Choose smallest binary as third body and ignore fudge factor.
          IF (SEMI.GT.SEMI0) THEN
              Q1 = MB0/MB
              AIN = SEMI
          ELSE
              Q1 = MB/MB0
              AIN = SEMI0
          END IF
          XFAC = (1.0 + Q1)*(1.0 + E1)/SQRT(1.0 - E1)
          PCRIT = 2.8*XFAC**0.4*AIN
          PMIN = SEMI1*(1.0 - E1)
          IF (PCRIT.LT.PMIN) THEN
              ITERM = -1
              WRITE (6,20)  AIN, SEMI1, E0, E1, RATIO, RB0, RB1,
     &                      PCRIT, PMIN
   20         FORMAT ('  STABQ:    AIN A1 E0 E1 RATIO R0 R1 PCR PM ',
     &                             1P,2E10.2,0P,2F7.3,F6.2,1P,4E9.1)
          END IF
      END IF
*
      RETURN
*
      END


! { dg-do run }
!
! PR 40176:  Fortran 2003: Procedure pointers with array return value
!
! Original test case by Barron Bichon <barron.bichon@swri.org>
! Modified by Janus Weil <janus@gcc.gnu.org>

PROGRAM test_prog

 TYPE ProcPointerType
   PROCEDURE(triple), POINTER, NOPASS :: f
 END TYPE ProcPointerType

 TYPE (ProcPointerType) :: ppt
 PROCEDURE(triple), POINTER :: f
 REAL :: tres(2)

 ppt%f => triple
 f => ppt%f
 tres = f(2,[2.,4.])
 if (abs(tres(1)-6.)>1E-3) STOP 1
 if (abs(tres(2)-12.)>1E-3) STOP 2
 tres = ppt%f(2,[3.,5.])
 if (abs(tres(1)-9.)>1E-3) STOP 3
 if (abs(tres(2)-15.)>1E-3) STOP 4

CONTAINS

 FUNCTION triple(n,x) RESULT(tre)
   INTEGER, INTENT(in) :: n
   REAL, INTENT(in) :: x(2)
   REAL :: tre(2)
   tre = 3.*x
 END FUNCTION triple

END PROGRAM test_prog


text:
      Czero5(z5(1,1,1,1,1),ep) =  - Cv(cc22222 + C234,ep) - 5*Cv(
     & cc12222 + C234,ep) - 10*Cv(cc11222 + C234,ep) - 10*Cv(cc11122 +
     & C234,ep) - 5*Cv(cc11112 + C234,ep) - Cv(cc11111 + C234,ep) - 5*
     & Cv(cc2222 + C234,ep) - 20*Cv(cc1222 + C234,ep) - 30*Cv(cc1122 +
     & C234,ep) - 20*Cv(cc1112 + C234,ep) - 5*Cv(cc1111 + C234,ep) - 10
     & *Cv(cc222 + C234,ep) - 30*Cv(cc122 + C234,ep) - 30*Cv(cc112 +
     & C234,ep) - 10*Cv(cc111 + C234,ep) - 10*Cv(cc22 + C234,ep) - 20*
     & Cv(cc12 + C234,ep) - 10*Cv(cc11 + C234,ep) - 5*Cv(cc2 + C234,ep)
     &  - 5*Cv(cc1 + C234,ep) - Cv(cc0 + C234,ep)

      Czero5(z5(1,1,1,1,2),ep) = Cv(cc12222 + C234,ep) + 4*Cv(cc11222
     &  + C234,ep) + 6*Cv(cc11122 + C234,ep) + 4*Cv(cc11112 + C234,ep)
     &  + Cv(cc11111 + C234,ep) + 4*Cv(cc1222 + C234,ep) + 12*Cv(cc1122
     &  + C234,ep) + 12*Cv(cc1112 + C234,ep) + 4*Cv(cc1111 + C234,ep)
     &  + 6*Cv(cc122 + C234,ep) + 12*Cv(cc112 + C234,ep) + 6*Cv(cc111
     &  + C234,ep) + 4*Cv(cc12 + C234,ep) + 4*Cv(cc11 + C234,ep) + Cv(
     & cc1 + C234,ep)

      Czero5(z5(1,1,1,1,3),ep) = Cv(cc22222 + C234,ep) + 4*Cv(cc12222
     &  + C234,ep) + 6*Cv(cc11222 + C234,ep) + 4*Cv(cc11122 + C234,ep)
     &  + Cv(cc11112 + C234,ep) + 4*Cv(cc2222 + C234,ep) + 12*Cv(cc1222
     &  + C234,ep) + 12*Cv(cc1122 + C234,ep) + 4*Cv(cc1112 + C234,ep)
     &  + 6*Cv(cc222 + C234,ep) + 12*Cv(cc122 + C234,ep) + 6*Cv(cc112
     &  + C234,ep) + 4*Cv(cc22 + C234,ep) + 4*Cv(cc12 + C234,ep) + Cv(
     & cc2 + C234,ep)

      Czero5(z5(1,1,1,2,2),ep) =  - Cv(cc11222 + C234,ep) - 3*Cv(
     & cc11122 + C234,ep) - 3*Cv(cc11112 + C234,ep) - Cv(cc11111 + C234
     & ,ep) - 3*Cv(cc1122 + C234,ep) - 6*Cv(cc1112 + C234,ep) - 3*Cv(
     & cc1111 + C234,ep) - 3*Cv(cc112 + C234,ep) - 3*Cv(cc111 + C234,ep
     & ) - Cv(cc11 + C234,ep)

      Czero5(z5(1,1,1,2,3),ep) =  - Cv(cc12222 + C234,ep) - 3*Cv(
     & cc11222 + C234,ep) - 3*Cv(cc11122 + C234,ep) - Cv(cc11112 + C234
     & ,ep) - 3*Cv(cc1222 + C234,ep) - 6*Cv(cc1122 + C234,ep) - 3*Cv(
     & cc1112 + C234,ep) - 3*Cv(cc122 + C234,ep) - 3*Cv(cc112 + C234,ep
     & ) - Cv(cc12 + C234,ep)

      Czero5(z5(1,1,1,3,3),ep) =  - Cv(cc22222 + C234,ep) - 3*Cv(
     & cc12222 + C234,ep) - 3*Cv(cc11222 + C234,ep) - Cv(cc11122 + C234
     & ,ep) - 3*Cv(cc2222 + C234,ep) - 6*Cv(cc1222 + C234,ep) - 3*Cv(
     & cc1122 + C234,ep) - 3*Cv(cc222 + C234,ep) - 3*Cv(cc122 + C234,ep
     & ) - Cv(cc22 + C234,ep)

      Czero5(z5(1,1,2,2,2),ep) = Cv(cc11122 + C234,ep) + 2*Cv(cc11112
     &  + C234,ep) + Cv(cc11111 + C234,ep) + 2*Cv(cc1112 + C234,ep) + 2
     & *Cv(cc1111 + C234,ep) + Cv(cc111 + C234,ep)

      Czero5(z5(1,1,2,2,3),ep) = Cv(cc11222 + C234,ep) + 2*Cv(cc11122
     &  + C234,ep) + Cv(cc11112 + C234,ep) + 2*Cv(cc1122 + C234,ep) + 2
     & *Cv(cc1112 + C234,ep) + Cv(cc112 + C234,ep)

      Czero5(z5(1,1,2,3,3),ep) = Cv(cc12222 + C234,ep) + 2*Cv(cc11222
     &  + C234,ep) + Cv(cc11122 + C234,ep) + 2*Cv(cc1222 + C234,ep) + 2
     & *Cv(cc1122 + C234,ep) + Cv(cc122 + C234,ep)

      Czero5(z5(1,1,3,3,3),ep) = Cv(cc22222 + C234,ep) + 2*Cv(cc12222
     &  + C234,ep) + Cv(cc11222 + C234,ep) + 2*Cv(cc2222 + C234,ep) + 2
     & *Cv(cc1222 + C234,ep) + Cv(cc222 + C234,ep)

      Czero5(z5(1,2,2,2,2),ep) =  - Cv(cc11112 + C234,ep) - Cv(cc11111
     &  + C234,ep) - Cv(cc1111 + C234,ep)

      Czero5(z5(1,2,2,2,3),ep) =  - Cv(cc11122 + C234,ep) - Cv(cc11112
     &  + C234,ep) - Cv(cc1112 + C234,ep)

      Czero5(z5(1,2,2,3,3),ep) =  - Cv(cc11222 + C234,ep) - Cv(cc11122
     &  + C234,ep) - Cv(cc1122 + C234,ep)

      Czero5(z5(1,2,3,3,3),ep) =  - Cv(cc12222 + C234,ep) - Cv(cc11222
     &  + C234,ep) - Cv(cc1222 + C234,ep)

      Czero5(z5(1,3,3,3,3),ep) =  - Cv(cc22222 + C234,ep) - Cv(cc12222
     &  + C234,ep) - Cv(cc2222 + C234,ep)

      Czero5(z5(2,2,2,2,2),ep) = Cv(cc11111 + C234,ep)

      Czero5(z5(2,2,2,2,3),ep) = Cv(cc11112 + C234,ep)

      Czero5(z5(2,2,2,3,3),ep) = Cv(cc11122 + C234,ep)

      Czero5(z5(2,2,3,3,3),ep) = Cv(cc11222 + C234,ep)

      Czero5(z5(2,3,3,3,3),ep) = Cv(cc12222 + C234,ep)

      Czero5(z5(3,3,3,3,3),ep) = Cv(cc22222 + C234,ep)

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


!-----------------------------------------------------------------------
subroutine sub_rhs_fwd
!-----------------------------------------------------------------------
  use mod_variables
  implicit none
!
  integer :: ic,jc
  real(8) :: dci,dci2,r12,dc12,dc12al
  real(8) :: xibl,yibl
!
  real(8) :: uxc,uxm,uxp,uym,uyp
  real(8) :: vyc,vxm,vxp,vym,vyp
  real(8) :: uxmm,uxpp,uymm,uyc,uypp
  real(8) :: vxmm,vxc,vxpp,vymm,vypp
!
  real(8) :: dudx1,dudy1
  real(8) :: dvdx1,dvdy1
  real(8) :: dudx2,dudy2
  real(8) :: dvdx2,dvdy2
!
  real(8) :: abs_uxc,abs_uyc
  real(8) :: abs_vxc,abs_vyc
!
  real(8) :: uold(imx1,jmax)
  real(8) :: vold(imax,jmx1)
!-----------------------------------------------------------------------
!
!
!---> Coefficients
  r12    = 1.d0/12.d0
  dci    = 1.d0/dx
  dci2   = dci**2/Re
  dc12   = dci*r12
  dc12al = dc12*3.d0
!
!
  udlt(1:imx1,1:jmax) = 0.d0
  vdlt(1:imax,1:jmx1) = 0.d0
  uold(1:imx1,1:jmax) = 0.d0
  vold(1:imax,1:jmx1) = 0.d0
!
!-----------------------------------------------------------------------
!---- x-direction
!-----------------------------------------------------------------------
    do jc=3,jmax-2
    do ic=3,imx1-2
!
      uxm  = ustg(ic-1,jc)
      uxp  = ustg(ic+1,jc)
      uym  = ustg(ic,jc-1)
      uyp  = ustg(ic,jc+1)
!
      uxmm = ustg(ic-2,jc)
      uxpp = ustg(ic+2,jc)
      uymm = ustg(ic,jc-2)
      uypp = ustg(ic,jc+2)
!
      uxc  = ustg(ic,jc)
      vxc  = 0.25d0*(vstg(ic  ,jc+1)+vstg(ic  ,jc) &
                    +vstg(ic-1,jc+1)+vstg(ic-1,jc))
!
      dudx1 = dc12  *(uxmm-8.d0*uxm         +8.d0*uxp-uxpp)
      dudy1 = dc12  *(uymm-8.d0*uym         +8.d0*uyp-uypp)
      dudx2 = dc12al*(uxmm-4.d0*uxm+6.d0*uxc-4.d0*uxp+uxpp)
      dudy2 = dc12al*(uymm-4.d0*uym+6.d0*uxc-4.d0*uyp+uypp)
!
      abs_uxc = dabs(uxc)
      abs_vxc = dabs(vxc)
!
!---- Inviscid term
      udlt(ic,jc) = (uxc*dudx1+abs_uxc*dudx2) &
                   +(vxc*dudy1+abs_vxc*dudy2)
!---- Viscous term
      udlt(ic,jc) = udlt(ic,jc)-(uxp-2.d0*uxc+uxm)*dci2 &
                               -(uyp-2.d0*uxc+uym)*dci2
!---- Pressure gradient term
      udlt(ic,jc) = udlt(ic,jc)+(pcnt(ic,jc)-pcnt(ic-1,jc))*dci
!
    enddo
    enddo
!#endif
!
!
!-----------------------------------------------------------------------
!---- y-direction
!-----------------------------------------------------------------------
    do jc=3,jmx1-2
    do ic=3,imax-2
!
      vxm  = vstg(ic-1,jc)
      vxp  = vstg(ic+1,jc)
      vym  = vstg(ic,jc-1)
      vyp  = vstg(ic,jc+1)
!
      vxmm = vstg(ic-2,jc)
      vxpp = vstg(ic+2,jc)
      vymm = vstg(ic,jc-2)
      vypp = vstg(ic,jc+2)
!
      uyc  = 0.25d0*(ustg(ic+1,jc  )+ustg(ic,jc  ) &
                    +ustg(ic+1,jc-1)+ustg(ic,jc-1))
      vyc  = vstg(ic,jc)
!
      dvdx1 = dc12  *(vxmm-8.d0*vxm         +8.d0*vxp-vxpp)
      dvdy1 = dc12  *(vymm-8.d0*vym         +8.d0*vyp-vypp)
      dvdx2 = dc12al*(vxmm-4.d0*vxm+6.d0*vyc-4.d0*vxp+vxpp)
      dvdy2 = dc12al*(vymm-4.d0*vym+6.d0*vyc-4.d0*vyp+vypp)
!
      abs_uyc = dabs(uyc)
      abs_vyc = dabs(vyc)
!
!---- Inviscid term
      vdlt(ic,jc) = (uyc*dvdx1+abs_uyc*dvdx2) &
                   +(vyc*dvdy1+abs_vyc*dvdy2)
!---- Viscous term
      vdlt(ic,jc) = vdlt(ic,jc)-(vxp-2.d0*vyc+vxm)*dci2 &
                               -(vyp-2.d0*vyc+vym)*dci2
!---- Pressure gradient term
      vdlt(ic,jc) = vdlt(ic,jc)+(pcnt(ic,jc)-pcnt(ic,jc-1))*dci
!
  enddo
  enddo
!
!
  uold(1:imx1,1:jmax) = ustg(1:imx1,1:jmax)
  vold(1:imax,1:jmx1) = vstg(1:imax,1:jmx1)
  ustg(1:imx1,1:jmax) = 0.d0
  vstg(1:imax,1:jmx1) = 0.d0
!
  do jc=3,jmx1-2
  do ic=3,imx1-2
    if(jc==jmx1-2) goto 10
    xibl = dble(min0(iblk(ic-1,jc),iblk(ic,jc)))
    ustg(ic,jc) = uold(ic,jc)-dt*udlt(ic,jc)*xibl
 10 continue
    if(ic==imx1-2) goto 20
    yibl = dble(min0(iblk(ic,jc-1),iblk(ic,jc)))
    vstg(ic,jc) = vold(ic,jc)-dt*vdlt(ic,jc)*yibl
 20 continue
  enddo
  enddo
!
  return
end subroutine sub_rhs_fwd
!-----------------------------------------------------------------------