text:
	SUBROUTINE GQMRKR  ( imark, imkhw, szmark, imkwid, iret )
C************************************************************************
C* GQMRKR								*
C*									*
C* This subroutine returns the current marker attributes including 	*
C* the marker number, the hardware/software flag and the marker 	*
C* size and line width.							*
C*									*
C* NOTE THAT THE CALLING SEQUENCE OF THIS SUBOUTINE HAS CHANGED FROM	*
C* PREVIOUS VERSIONS OF GEMPAK.						*
C*									*
C* GQMRKR  ( IMARK, IMKHW, SZMARK, IMKWID, IRET )			*
C*									*
C* Output parameters:							*
C*	IMARK		INTEGER		Marker number			*
C*	IMKHW		INTEGER		Hardware/software flag		*
C*					  1 = software			*
C*					  2 = hardware			*
C*	SZMARK		REAL		Marker size multiplier		*
C*	IMKWID		INTEGER		Marker width multiplier		*
C*	IRET		INTEGER		Return code			*
C**									*
C* Log:									*
C* I. Graffman/RDS	 5/85	GEMPLT Version 3.1			*
C* M. desJardins/GSFC	 6/88	Documentation				*
C* S. Schotz/GSC	 1/90	Added marker width			*
C************************************************************************
	INCLUDE		'ERROR.PRM'
	INCLUDE		'DEVSET.CMN'
	INCLUDE		'DEVCHR.CMN'
C------------------------------------------------------------------------
C* 	If device has not been set, return an error.
C
	IF  ( ddev .eq. ' ' )  THEN
	    imark  = 0
	    imkhw  = 0
	    szmark = 0
	    imkwid = 0
	    iret   = NDVICE
	  ELSE
C
C*	    Retrieve values from /DEVSET/.
C
	    imark  = lmark
	    imkhw  = lmkhw
	    szmark = smksz
	    imkwid = lmkwid
	    iret   = NORMAL
	END IF
C*
	RETURN
	END


text:
subroutine get_num_satobs(obsfile,npred,num_obs_tot,endian)
    use read_diag, only: diag_data_fix_list,diag_header_fix_list,diag_header_chan_list, &
    diag_data_chan_list,diag_data_extra_list,read_radiag_data,read_radiag_header, &
    diag_data_name_list
    implicit none
    character (len=500), intent(in) :: obsfile
    integer, intent(out) :: num_obs_tot, npred
    integer iunit, iflag, ios, n
    logical fexist,lretrieval,lverbose
    real(8) :: errorlimit,errorlimit2
    character(len=6), optional, intent(in) :: endian

    type(diag_header_fix_list )         :: header_fix0
    type(diag_header_chan_list),allocatable :: header_chan0(:)
    type(diag_data_fix_list   )         :: data_fix0
    type(diag_data_chan_list  ),allocatable :: data_chan0(:)
    type(diag_data_extra_list) ,allocatable :: data_extra0(:,:)
    type(diag_data_name_list)           :: data_name0

    character(len=6) :: convert_endian

    if (.not. present(endian)) then
      convert_endian = 'native'
    else
      convert_endian = endian
    endif

    errorlimit=1./sqrt(1.e9)
    errorlimit2=1./sqrt(1.e-6)
    iunit = 7
    lretrieval=.false.
    lverbose=.false.

    num_obs_tot = 0

    inquire(file=obsfile,exist=fexist)
    if (.not.fexist) goto 900

    if (trim(convert_endian) == 'big') then
       open(iunit,form="unformatted",file=trim(obsfile),iostat=ios,convert='big_endian')
    else if (trim(convert_endian) == 'little') then
       open(iunit,form="unformatted",file=trim(obsfile),iostat=ios,convert='little_endian')
    else if (trim(convert_endian) == 'native') then
       open(iunit,form="unformatted",file=trim(obsfile),iostat=ios)
    endif
    rewind(iunit)
    call read_radiag_header(iunit,lretrieval,header_fix0,header_chan0,data_name0,iflag,lverbose)
    npred = header_fix0%npred

    do
       call read_radiag_data(iunit,header_fix0,lretrieval,data_fix0,data_chan0,data_extra0,iflag )
       if( iflag /= 0 )exit
       chan: do n=1,header_fix0%nchan
         num_obs_tot = num_obs_tot + 1
       end do chan
    enddo
900 continue
    close(iunit)

end subroutine get_num_satobs

subroutine get_satobs_data(obsfile, nobs_max, npred, h_x, h_xnobc, x_obs, x_err, &
           x_lon, x_lat, x_time, x_channum, x_errorig, x_biaspred, x_use, x_qcmark, &
           x_water_frac, x_land_frac, x_ice_frac, x_snow_frac, endian)
  use read_diag, only: diag_data_fix_list,diag_header_fix_list,diag_header_chan_list, &
  diag_data_chan_list,diag_data_extra_list,read_radiag_data,read_radiag_header, &
  diag_data_name_list
  implicit none
  character(len=6), optional, intent(in) :: endian
  character*500, intent(in) :: obsfile
  integer, intent(in) :: nobs_max, npred
  real, dimension(nobs_max), intent(out) :: h_x,h_xnobc,x_obs,x_err,x_lon,&
                               x_lat,x_time,x_errorig
  integer, dimension(nobs_max), intent(out) :: x_use, x_qcmark
  real, dimension(nobs_max), intent(out) :: x_water_frac, x_land_frac, x_ice_frac, x_snow_frac
  real, dimension(npred+1,nobs_max), intent(out) :: x_biaspred
  integer, dimension(nobs_max), intent(out) ::  x_channum

  integer iunit, iflag, nobs, n, ios
  logical fexist,lretrieval,lverbose,adp_anglebc,emiss_bc
  real(8) :: errorlimit, errorlimit2

  type(diag_header_fix_list )         :: header_fix
  type(diag_header_chan_list),allocatable :: header_chan(:)
  type(diag_data_fix_list   )         :: data_fix
  type(diag_data_chan_list  ),allocatable :: data_chan(:)
  type(diag_data_extra_list) ,allocatable :: data_extra(:,:)
  type(diag_data_name_list)           :: data_name

  character(len=6) :: convert_endian

  if (.not. present(endian)) then
    convert_endian = 'native'
  else
    convert_endian = endian
  endif

  errorlimit=1./sqrt(1.e9)
  errorlimit2=1./sqrt(1.e-6)

  iunit = 7
  lretrieval=.false.
  lverbose=.false.
  emiss_bc = .true.
  adp_anglebc = .true.
  x_biaspred = 0

  nobs = 0

  inquire(file=obsfile,exist=fexist)
  if(.not.fexist) go to 900

  if (trim(convert_endian) == 'big') then
     open(iunit,form="unformatted",file=trim(obsfile),iostat=ios,convert='big_endian')
  else if (trim(convert_endian) == 'little') then
     open(iunit,form="unformatted",file=trim(obsfile),iostat=ios,convert='little_endian')
  else if (trim(convert_endian) == 'native') then
     open(iunit,form="unformatted",file=trim(obsfile),iostat=ios)
  endif
  rewind(iunit)
  call read_radiag_header(iunit,lretrieval,header_fix,header_chan,data_name,iflag,lverbose)

  do
   call read_radiag_data(iunit,header_fix,lretrieval,data_fix,data_chan,data_extra,iflag )
   if( iflag /= 0 ) exit
   chan:do n=1,header_fix%nchan
      nobs = nobs + 1
      if (nobs > nobs_max) then
          print *,'warning:  exceeding array bounds in get_satobs_data',&
          nobs,nobs_max
      end if
      x_channum(nobs) = n
      x_lon(nobs) = data_fix%lon
      x_lat(nobs) = data_fix%lat
      x_time(nobs) = data_fix%obstime
      x_obs(nobs) = data_chan(n)%tbobs
      x_use(nobs) = header_chan(n)%iuse
      x_qcmark(nobs) = data_chan(n)%qcmark
      ! bias corrected Hx
      h_x(nobs) = x_obs(nobs) - data_chan(n)%omgbc
      ! un-bias corrected Hx
      h_xnobc(nobs) = x_obs(nobs) - data_chan(n)%omgnbc
      ! data_chan%errinv is inverse error variance.
      x_errorig(nobs) = header_chan(n)%varch**2
      x_err(nobs) = (1./data_chan(n)%errinv)**2

      x_biaspred(1,nobs) = data_chan(n)%bifix(1) ! fixed angle dependent bias
      x_biaspred(2,nobs) = data_chan(n)%bicons ! constant bias correction
      x_biaspred(3,nobs) = data_chan(n)%biang ! scan angle bias correction
      x_biaspred(4,nobs) = data_chan(n)%biclw ! CLW bias correction
      x_biaspred(5,nobs) = data_chan(n)%bilap2 ! square lapse rate bias corr
      x_biaspred(6,nobs) = data_chan(n)%bilap ! lapse rate bias correction
      if (npred == 7) then
        x_biaspred(7,nobs) = data_chan(n)%bicos ! node*cos(lat) bias correction for SSMIS
        x_biaspred(8,nobs) = data_chan(n)%bisin ! sin(lat) bias correction for SSMIS
      endif
      if (emiss_bc) x_biaspred(9,nobs) = data_chan(n)%biemis

      if (adp_anglebc) then
         x_biaspred( 1,nobs)  = data_chan(n)%bifix(5) ! fixed angle dependent bias correction
         x_biaspred(npred-2,nobs)  = data_chan(n)%bifix(1) ! 4th order scan angle (predictor)
         x_biaspred(npred-1,nobs)  = data_chan(n)%bifix(2) ! 3rd order scan angle (predictor)
         x_biaspred(npred,nobs)  = data_chan(n)%bifix(3) ! 2nd order scan angle (predictor)
         x_biaspred(npred+1,nobs)    = data_chan(n)%bifix(4) ! 1st order scan angle (predictor)
      endif

      x_water_frac(nobs) = data_fix%water_frac ! fractional coverage by water
      x_land_frac(nobs)  = data_fix%land_frac  ! fractional coverage by land
      x_ice_frac(nobs)   = data_fix%ice_frac   ! fractional coverage by ice
      x_snow_frac(nobs)  = data_fix%snow_frac  ! fractional coverage by snow

  enddo chan
 enddo

900  continue
 close(iunit)

end subroutine get_satobs_data

max_stars_repo_name: whpldgs/py-gsidiag

max_issues_repo_name: whpldgs/py-gsidiag

max_forks_repo_name: whpldgs/py-gsidiag


text:
module m_Gutenberg
  use m_Useful
  implicit none
  private
!
  public :: PrintCurrents
!
contains
!
  subroutine PrintCurrents (cxz, unitout)
    character (len=*), parameter :: myname = "PrintCurrents"
    type (cxzType), intent (inout) :: cxz
    integer, intent (inout) :: unitout
!
    character (len=k_ml) :: fmtContainer, saux
    integer :: i
!
    write (unitout, '(a1,a16)', advance="no") "#", "Time "
    do i = 1, cxz%natoms
      write (unitout, '(a4,i8,a1,a2,a1)', advance="no") "Curr", i + 1, "-", trim (cxz%frames(1)%element(i)), ";"
    end do
    do i = 1, cxz%natoms
      write (unitout, '(a4,i8,a1,a2,a1)', advance="no") "GCur", i + 1 + cxz%natoms, "-", trim (cxz%frames(1)%element(i)), ";"
    end do
    write (unitout,*)
    write (unitout,'(a)') "#to get the atom number for Current component substract 1 from last number after Curr"
    write (unitout, '(a,i0,a)') "#to get the atom number for Gamma Current component substract ", 1 + cxz%natoms, " from last numbe&
   &r after GCur"
    write (fmtContainer, '("(f16.8,1x,",i0,"g,",i0,"g)")') cxz%natoms, cxz%natoms
    do i = 1, cxz%nframes
      write (unitout, trim(fmtContainer)) cxz%frames(i)%timestamp, cxz%frames(i)%curr1
    end do
  end subroutine PrintCurrents
!
!
end module m_Gutenberg

max_stars_repo_name: naumruso/tdtbuj

max_issues_repo_name: naumruso/tdtbuj

max_forks_repo_name: naumruso/tdtbuj


text:
C
C $Id: intrvpdp.f,v 1.3 2008-07-27 03:10:11 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
C NOTE: If you make any changes to this software, please remember to
C make the same changes to the corresponding single precision routine.
C
      FUNCTION INTRVPDP(T,X,N,P,TP)
      DOUBLE PRECISION TT
c
      INTEGER N
      DOUBLE PRECISION T,X(N),P,TP
c
c                                 coded by alan kaylor cline
c                           from fitpack -- january 26, 1987
c                        a curve and surface fitting package
c                      a product of pleasant valley software
c                  8603 altus cove, austin, texas 78759, usa
c
c this function determines the index of the interval
c (determined by a given increasing sequence) in which a
c given value lies, after translating the value to within
c the correct period.  it also returns this translated value.
c
c on input--
c
c   t is the given value.
c
c   x is a vector of strictly increasing values.
c
c   n is the length of x (n .ge. 2).
c
c and
c
c   p contains the period.
c
c on output--
c
c   tp contains a translated value of t (i. e. x(1) .le. tp,
c   tp .lt. x(1)+p, and tp = t + k*p for some integer k).
c
c   intrvldp returns an integer i such that
c
c          i = 1       if             tp .lt. x(2)  ,
c          i = n       if   x(n) .le. tp            ,
c          otherwise       x(i)  .le. tp .lt. x(i+1),
c
c none of the input parameters are altered.
c
c-----------------------------------------------------------
c
      SAVE I
      DATA I/1/
c
      NPER = (T-X(1))/P
      TP = T - DBLE(NPER)*P
      IF (TP.LT.X(1)) TP = TP + P
      TT = TP
c
c check for illegal i
c
      IF (I.GE.N) I = N/2
c
c check old interval and extremes
c
      IF (TT.LT.X(I)) THEN
          IF (TT.LE.X(2)) THEN
              I = 1
              INTRVPDP = 1
              RETURN
          ELSE
              IL = 2
              IH = I
          END IF
      ELSE IF (TT.LE.X(I+1)) THEN
          INTRVPDP = I
          RETURN
      ELSE IF (TT.GE.X(N)) THEN
          I = N
          INTRVPDP = N
          RETURN
      ELSE
          IL = I + 1
          IH = N
      END IF
c
c binary search loop
c
    1 I = (IL+IH)/2
      IF (TT.LT.X(I)) THEN
          IH = I
      ELSE IF (TT.GT.X(I+1)) THEN
          IL = I + 1
      ELSE
          INTRVPDP = I
          RETURN
      END IF
      GO TO 1
      END

ext: f

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
! Advent of Code 2021 : https://adventofcode.com/2021
! Puzzle for Day 14 part 01 :
! Build (or run) with:
! gfortran -static-libgcc -o aoc_day14_p1 aoc_day14_p1.f90
! fpm run aoc_day14_p1

program aoc_day14_p1

    use, intrinsic :: iso_fortran_env, stderr => error_unit
    implicit none

    ! Declare and set constants the AOC day and part
    integer, parameter :: day = 14, part = 1

    ! Declare file related variables
    ! - Prod uses: day14-input.txt
    ! - Test uses: day14-TEST-input.txt
    character(len = *), parameter :: filename = "data/day14-input.txt"

    ! Variables to track and manage the program
    integer :: total_lines = 0
    integer :: i,j,n,m
    character(len=100) :: poly_temp

    ! Output startup message
    call display_welcome(day, part)

    print *, ""
    print *, "INCOMPLETE - work in progress still..."
    print *, ""

    ! Check input puzzle data file and get its line count
    total_lines = get_line_count(filename)
    if (total_lines == -1) stop ": Input file not found."
    write (*, '(A,A,A,I3.3)') "Filename: '", filename, "' line count is: ", total_lines

contains

    ! Display the AoC puzzle day and part being resolved
    subroutine display_welcome(day, part)
        implicit none
        integer, intent(in) :: day, part
        print *, ""
        write (*, '(A,I2.2,A,I2.2)') "AoC 2021 - solution for Day ", day, " part ", part
        print *, ""
    end subroutine display_welcome

    ! If the file exists return its total line count or '-1' on error
    ! TODO: add optional flag to run 'stop' instead of returning '-1' if needs to be a fatal error?
    function get_line_count(filename) result (lines)
        implicit none
        character(len = *), intent(in) :: filename
        integer :: fileunit, io
        integer :: lines
        logical :: exists

        inquire(file = filename, exist = exists)
        if (.not. exists) then
            ! below should be printed to stderr!!
            write (stderr, '(A,A)') "ERROR: unable to locate file name: ", filename
            lines = -1
            return
        end if

        open(newunit = fileunit, file = filename, action = 'read', position = 'rewind')
        lines = 0
        do
            read(fileunit, *, iostat = io)
            if (io /= 0) exit
            lines = lines + 1
        end do

        close(fileunit)
    end function get_line_count

end program aoc_day14_p1

max_stars_repo_name: wiremoons/AdventOfCode

max_issues_repo_name: wiremoons/AdventOfCode


text:
      SUBROUTINE ARSCND( T )
*
      REAL       T
*
*  -- LAPACK auxiliary routine (preliminary version) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     July 26, 1991
*
*  Purpose
*  =======
*
*  SECOND returns the user time for a process in seconds.
*  This version gets the time from the system function CPU_TIME.
*

      CALL CPU_TIME(T)

      RETURN
*
*     End of SECOND
*
      END

max_stars_repo_name: sargas/scipy

max_issues_repo_name: sargas/scipy

max_forks_repo_name: sargas/scipy


text:
 About Alex
Alex is currently a second year student majoring in Environmental Toxicology with a specialization in Public Health and possibly a minor in Sexuality Studies. He can be found either in the 24 Hour Reading Room, Physical Sciences Library, or wandering aimlessly around campus. Alex is dedicated to being involved in campus community; particularly sexual assault awareness and prevention through education and outreach through Men Acting Against Rape(MAAR). He is currently the President of Men Acting Against Rape, a member of Students Against Sexual Violence, and well acquainted with the good people at Campus Violence Prevention Program.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      subroutine test_mpi3rma_get_acc_maxloc_lock( base, array_size, DIM1, win, comm )

      use utilities, only: check_fompi_status

#if INCLUDE_TYPE==1
      implicit none

      include 'fompif.h'
#endif
#if INCLUDE_TYPE==2
      use fompi

      implicit none
#endif

      integer, intent(in) :: array_size
      real, dimension(2,array_size), intent(inout) :: base
      integer, intent(in) :: DIM1
      integer, intent(in) :: win
      integer, intent(in) :: comm

      integer :: i
      integer :: rank
      integer :: commsize
      integer :: mpierr
      integer(kind=MPI_ADDRESS_KIND) :: start

      logical :: flag
      integer(kind=MPI_ADDRESS_KIND) :: create_flavor

      real, dimension(2,array_size) :: result_buf

      real, parameter :: fompi_epsilon = 0.00000000001

      character(50) :: filename

      write (filename,*) "get_acc_maxloc_lock.f90"

      call MPI_Comm_rank( comm, rank, mpierr )
      call MPI_Comm_size( comm, commsize, mpierr )

      call foMPI_Win_lock( foMPI_LOCK_EXCLUSIVE, rank, 0, win, mpierr )
      call check_fompi_status(mpierr, MPI_SUCCESS, filename, 34, rank)

      if ( ((commsize/2) .eq. rank) .or. ((commsize/2+1) .eq. rank) ) then
        do i=1,DIM1
          base(1,i) = real(commsize*2+i)
          base(2,i) = real(rank)
        enddo
      else
        do i=1,DIM1
          base(1,i) = real(commsize+rank+i)
          base(2,i) = real(rank)
        enddo
      endif

      call foMPI_Win_unlock( rank, win, mpierr )
      call check_fompi_status(mpierr, MPI_SUCCESS, filename, 49, rank)

      call foMPI_Win_get_attr( win, foMPI_WIN_CREATE_FLAVOR, create_flavor, flag, mpierr )
      call check_fompi_status(mpierr, MPI_SUCCESS, filename, 52, rank)
      if (flag .neqv. .true.) then
        write (*,*) "Didn't found foMPI_WIN_CREATE_FLAVOR keyval."
        STOP
      endif

      if ( create_flavor .eq. foMPI_WIN_FLAVOR_DYNAMIC ) then
        if (rank .eq. 0 ) then
          call foMPI_Get_address( base(1,1), start, mpierr )
        endif
        call MPI_Bcast( start, 1, MPI_AINT, 0, comm, mpierr ) ! also the signal that the initialization is finished
      else
        start = 0
        call MPI_Barrier( comm, mpierr ) ! guarantees that the initialization is finished
      endif

      if (rank .ne. 0 ) then
        call foMPI_Win_lock( foMPI_LOCK_SHARED, 0, 0, win, mpierr )
        call check_fompi_status(mpierr, MPI_SUCCESS, filename, 70, rank)

        call foMPI_Get_accumulate( base(1,1), DIM1, MPI_2REAL, result_buf(1,1), DIM1, MPI_2REAL, 0, start, DIM1, MPI_2REAL, &
          foMPI_MAXLOC, win, mpierr )
        call check_fompi_status(mpierr, MPI_SUCCESS, filename, 73, rank)

        call foMPI_Win_unlock( 0, win, mpierr )
        call check_fompi_status(mpierr, MPI_SUCCESS, filename, 76, rank)
      endif

      call MPI_Barrier( comm, mpierr ) ! wait for the end of the accesses

      if (rank .eq. 0) then
        do i=1,DIM1
          if ( (base(1,i) - real(commsize*2+i) .gt. fompi_epsilon) .or. (int(base(2,i)) .ne. (commsize/2)) ) then
            write (*,*) "Error on rank ", rank, " in ", trim(filename), &
              " line ", 83, ": max element ", i, " is ", base(1,i), &
              " instead of ", real(commsize*2+i), " or index is ", int(base(2,i)), &
              " instead of ", commsize/2, "."
            STOP
          endif
        enddo
      endif

      if ( rank .eq. 0 ) then
        write (*,*) "No Error"
      endif

      end subroutine test_mpi3rma_get_acc_maxloc_lock

max_stars_repo_name: spcl/foMPI-NA

max_issues_repo_name: spcl/foMPI-NA

max_forks_repo_name: spcl/foMPI-NA


text:
!********************************************************************************
!* This file is part of:
!* BAFFLE (https://github.com/ATHannington/BAFFLE)
!* BAFFLE: BAsic Fortran FLuids Engine
!*
!* BAFFLE:
!* Copyright (C) 2017 Andrew Hannington (ath4@st-andrews.ac.uk)
!*
!* This software is a derivative of work by Bert Vandenbroucke
!* (bert.vandenbroucke@gmail.com)
!* Find more work by Bert Vandenbroucke at: (https://github.com/bwvdnbro)
!*
!* BAFFLE is a free software: you can redistribute it and/or modify it under the
!* terms of the GNU Affero General Public License
!* as published by the Free Software Foundation, either version 3 of the License, or
!* (at your option) any later version.
!*
!* BAFFLE is distributed in the hope that it will
!* be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
!* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!* GNU Affero General Public License for more details.
!*
!* You should have received a copy of the GNU Affero General Public License
!* along with BAFFLE. If not, see
!* <http://www.gnu.org/licenses/>.
!********************************************************************************

!********************************************************************************
!* @file constants_module_1D.f90
!*
!* @Constants module for 1D version of BAFFLE. Can control much of the behaviour
!* of BAFFLE via this module: Fortran 90 version.
!*
!* @author Andrew Hannington (ath4@st-andrews.ac.uk)
!********************************************************************************

!-------------------------------------------------------------------------------
module constants								!Module declaring all global constants

!	Cells format: 	cells(i,j,k,1:29) = (/m,qx,qy,qz,E,rho,ux,uy,uz,p,V,Midx,Midy,Midz, &
!						& grad_d_x,grad_d_y,grad_d_z, &
!						& grad_ux_x,grad_uy_x,grad_uz_x, &
!						& grad_ux_y,grad_uy_y,grad_uz_y, &
!						& grad_ux_z,grad_uy_z,grad_uz_z, &
!						& grad_p_x,grad_p_y,grad_p_z /)

!	Define constants:
	Integer, save :: NcXr,NcX,NcYr,NcY,NcZr,NcZ,Ncellsg=2,Ncells		!Number of Cells: real, in x,y and z,(NcXr etc.) and their total with ghost counterparts (NcX etc.);
										!ghost (for implementation of boundary conditions); total
	Real*8, allocatable,dimension(:,:,:,:) :: cells				!Cells of simulation array over which the variables are sampled
	Real*8 :: xrangeMax,xrangeMin						!Min and max of x-axis range of test volume
	Real*8 :: yrangeMax,yrangeMin						!Min and max of y-axis range of test volume
	Real*8 :: zrangeMax,zrangeMin						!Min and max of z-axis range of test volume
	Real*8 :: gam								!Gamma, the ratio of specific heats Cp/Cv
	Real*8 :: delta_t,delt_default=1.d-4						!Time step for Time Integration, and default to prevent zeros being used.
	Real*8 :: t								!Time up to which the integrator will proceed
	Real*8 :: CellVol							!Cell Volume
	Real*8 :: C_CFL								!C_CFL: Courant-Friedrichs-Lewy number such that 0<C_CFL<1
	Real*8 :: Cs								!Sound speed
	Real*8 :: xArea, yArea, zArea						!Surface area for x,y,z direction faces of cuboid.
	Integer :: mn,qxn,qyn,qzn,En,rhon,uxn,uyn,uzn,pn,Vn,xn,yn,zn, &		!Integers to hold the position of each in the d-dimension (cells(i,j,k,d)) of cells array.
 		 & gdxn,gdyn,gdzn, &
		 & guxxn,guyxn,guzxn, &
		 & guxyn,guyyn,guzyn, &
		 & guxzn,guyzn,guzzn, &
	 	 & gpxn,gpyn,gpzn

	Real*8 :: Pi=DAcos(-1.d0)

!	Set save file output directory:
	character(len=100):: SavePath = "/home/alpha/ath4/summer2017/HydroCode/D_TEST/"
!***PLEASE NOTE***: A path string of greater than 100 characters will NOT be tolerated within the program unless you edit the hard-coded '(A100,~)' format identifiers in the write_data subroutine in hydrocode-main, to '(Ax,~)' where x>100
!	Set File format of outputs
	character(len=100) :: SaveFileFormat = ".dat"

!	Switch questions on or off. Reverts to following defaults if off.
!	Off is useful for running the program as a background job.
	Logical :: Qbool=.FALSE.						!Toggle to switch runtime questions on or off:	true => questions off	false => questions on

!	Settings for switching between accuracy orders:
	Logical :: orderBool							!Integer to hold what order of evaluation has been selected. FALSE => First Order TRUE=> Second Order
	Logical :: orderBooldefault= .true.					!Second Order: FALSE => First Order TRUE=> Second Order

!	Some system constants defaults:
	Real*8 ::NcXrdefault=100.d0,NcYrdefault=1.d0,NcZrdefault=1.d0
	Real*8 :: tdefault=0.2d0
	Real*8 :: gamDefault=5.0d0/3.d0
	Real*8 :: CFLDefault=0.4d0

!	Frame writing settings:
	Logical :: writeFramesbool						!Integer to hold whether to write frames or just write at end. true => write frames ; false => write at end only
	Logical :: wFramesDefault= .TRUE.					! true => Write Frames	 false => Don't write Frames
	Logical :: delt_from_frames_bool = .TRUE.				!TRUE => Adjust timesteps with number of frames (can dramatically slow program, but will enable more visualisation at 											small final times). This also ensures the number of frames selected
										!FALSE => Do not adjust the timestep from number of frames irrespective of whether frames are selected.
!	Set Max number of Frames to be written:
	Integer :: MaxFrames = 100

!	Isothermal Settings:
	Logical :: IsoBool							!" True=> Isothermal Setup ON	False => Adiabatic Setup ON"
	Logical :: IsoBoolDefault= .FALSE.					!" True=> Isothermal Setup ON	False => Adiabatic Setup ON"
	Real*8 :: CsIso								!Isothermal Global Sound-Speed
	Real*8 :: CsIsoDefault= DSQRT(0.1d0/0.5d0)				!Default Isothermal Global sounds speed. Edit this based on CsIso = SQRT(p_desired / RhoMaxSet)
	LOGICAL :: TwoPhaseBool = .FALSE.					!Toggles Two Phase Isothermal Setup. TRUE => Two Phase On	FALSE => Two Phase Off
	Real*8 :: IsoDeltaTScaler = 1.d0/DSQRT(5.d0/3.d0)			!Scales Isothermal timestep to be equivalent to that of Adiabtic case.

!	Random Number Generator Seed:
	Integer*4 :: RandSeed
	Integer*4 :: RandSeedDefault=1234

!	Set maximum and minimum density, and diameter of blobs:
!	Diameter:
	Real*8 :: DmaxSet=0.2d0
!	Density:
	Real*8 :: RhomaxSet=2.0d0
	Real*8 :: RhominSet=0.5d0
	Real*8 :: RhoCrit							!Value Halfway between Rho min and max that sets the minimum value of Rho that Blobs can take

!	Call to Random Blobs Bool:
	Logical :: RandBlobBool=.TRUE.						!This switch allows the random blob regions to be defined once, on the first call to init_conds, and no other time after 											that. 	***MUST ALWAYS BE TRUE***
	LOGICAL :: BlobSwitchBool = .TRUE.					!This toggles switches whether blobs are to be generated or not.	TRUE => Blobs Generated		FALSE => BLobs 											NOT Generated

	Integer :: NumBlobs = 25						!Number of blobs to be generated in the random blob generator
	Real*8, allocatable,dimension (:) :: x0_blob_array,y0_blob_array, &	!Arrays to hold the parameters of the blobs globally such that they can be used on every call of init_conds
				& z0_blob_array, d_blob_array, rho_blob_array

end module constants

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

max_stars_repo_name: ATHannington/BAFFLE

max_issues_repo_name: ATHannington/BAFFLE

max_forks_repo_name: ATHannington/BAFFLE


text:
Dave Reed is the current Noble Grand of the Davis Chapter of the International Order of Odd Fellows.

For more information and published articles about David Reed, visit:

http://lesliebrodie.blog.co.uk/tags/davidreed/fullposts/


text:
      program t_m2h
      use mod_za  ! HYCOM array I/O interface
      implicit none
c
      integer   i,idim,j,jdim,l,length,nfill,nzero
      real      dhmin,dhmax
      character preambl(5)*79
c
c --- read in a micom topography file,
c --- write out the corresponding hycom topography file.
c --- new header read from stdin (unit 5).
c
      integer,   allocatable :: ip(:,:)
      real,      allocatable :: dh(:,:),dm(:,:)
      character, allocatable :: util(:)*2
c
      call xcspmd  !input idm,jdm
      allocate( ip(idm,jdm) )
      allocate( dh(idm,jdm) )
      allocate( util(idm*jdm+14) )
c
      read(51,'(A79)') preambl
c
      write(6,*)       'MICOM header:'
      write(6,'(A79)') preambl
      call flush(6)
c
      read(51,   *    ) idim,jdim,length
      read(51,'(40a2)') (util(l),l=1,length)
c
      if     (idim.lt.jdm-1 .or. jdim.lt.idm-1) then
        write(6,*)
        write(6,*) 'error expected idim,jdim >= ',jdm-1,idm-1
        write(6,*) '         input idim,jdim  = ',idim,jdim
        write(6,*)
        call flush(6)
        stop
      elseif (idim.gt.jdm   .or. jdim.gt.idm  ) then
        write(6,*)
        write(6,*) 'error expected idim,jdim <= ',jdm,idm
        write(6,*) '         input idim,jdim  = ',idim,jdim
        write(6,*)
        call flush(6)
        stop
      endif
      allocate( dm(idim,jdim) )
      call unpakk(dm,idim,idim,jdim,util,length)
      close(unit=51, status='keep')
c
c --- new header read from stdin (unit 5).
c
      read(5,'(A79)') preambl
c
      write(6, *)
      write(6, *)       'HYCOM header:'
      write(6, '(A79)') preambl
      call flush(6)
      write(61,'(A79)') preambl
c
c --- micom is (NS,EW), and hycom is (EW,SN).
c
      do j=1,jdm
        do i=1,idm
          dh(i,j) = 0.0
        enddo
      enddo
      do j= 1,idim
        do i= 1,jdim
          dh(i,j) = dm(jdm-j,i)
        enddo
      enddo
c
c --- fill single-width inlets and 1-point seas.
c
 100  continue
      nfill=0
      do j=1,jdm-1
        do i=1,idm-1
          nzero=0
          if (dh(i,j).gt.0.0) then
            if (i.eq.    1.or.dh(i-1,j).le.0.0) nzero=nzero+1
            if (i.eq.idm-1.or.dh(i+1,j).le.0.0) nzero=nzero+1
            if (j.eq.    1.or.dh(i,j-1).le.0.0) nzero=nzero+1
            if (j.eq.jdm-1.or.dh(i,j+1).le.0.0) nzero=nzero+1
            if (nzero.ge.3) then
              write (6,'(a,i4,a,i4,a,i1,a)')
     +          ' dh(',i,',',j,') set to zero (',
     +          nzero,' land nieghbours)'
              dh(i,j)=0.0
              nfill=nfill+1
            end if
          end if
        enddo
      enddo
      if (nfill.gt.0) go to 100
c
c --- mask.
c
      do j= 1,jdm
        do i= 1,idm
          ip(i,j) = 0
        enddo
      enddo
      do j= 1,jdm
        do i= 1,idm
          if (dh(i,j).gt.0.0) then
            ip(i,j) = 1
          endif
        enddo
      enddo
c
c --- write out the hycom bathymetry.
c
      call zaiost
      call zaiopn('new', 61)
      call zaiowr(dh, ip,.true., dhmin,dhmax, 61, .false.)
      write(61,6100) dhmin,dhmax
      write(6, 6100) dhmin,dhmax
      write(6, *)
 6100 format('min,max depth = ',2f10.3)
      end

max_stars_repo_name: TillRasmussen/HYCOM-tools

max_issues_repo_name: TillRasmussen/HYCOM-tools

max_forks_repo_name: TillRasmussen/HYCOM-tools


text:
	SUBROUTINE DENDD  ( ieop, iret )
C************************************************************************
C* DENDD								*
C*									*
C* This subroutine must be the last subroutine called by any program 	*
C* that uses GEMPLT.  It will flush internal buffers and close files	*
C* if necessary in the device driver.  The DEVICE subprocess may be 	*
C* retained so that the current definitions are available in later 	*
C* programs.								*
C*									*
C* DENDD  ( IEOP, IRET )						*
C*									*
C* Input parameters:							*
C*	IEOP		INTEGER		End plotting flag		*
C*					  0 = retain subprocess		*
C*					  1 = stop subprocess		*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C**									*
C* Log:                                                                 *
C* M. Vilardo/RDS	 6/84	GEMPLT Version 3.0                      *
C* M. Goodman/RDS	 5/85	GEMPLT Version 3.1                      *
C* M. desJardins/GSFC	 5/88	Documentation				*
C************************************************************************
	INCLUDE		'FUNCCODE.PRM'
	INCLUDE		'ERROR.PRM'
	INTEGER		isend ( 3 )
C------------------------------------------------------------------------
C*	Load input parameters into buffer and write them to the mailbox.
C
	isend ( 1 ) = 3
	isend ( 2 ) = CENDD
	isend ( 3 ) = ieop
C
	CALL GPUT  ( isend, 3, iret )
	IF  ( iret .eq. NORMAL )  THEN
C
C*	    Get the return code.
C
	    CALL GGET  ( iret, 1, ier )
	    IF  ( ier .ne. NORMAL )  iret = ier
	END IF
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
http://i1291.photobucket.com/albums/b544/davisartcenter/PeteScully2012holidaysale_zps5ce2635d.jpg (Drawing by Pete Scully)

A wellestablished Davis tradition since 1985, the Davis Art Center Holiday Sale is nonetheless getting a fresh facelift in 2013.
April is the deadline to apply as an artisan vendor at our annual Holiday Sale fundraiser in December.  Usually due in August, these applications must be submitted by April 30, 2013.
Due to the high volume of talented artists who apply each year, the Art Center’s Holiday Sale Selection Committee decided that an earlier start would help the group ensure a diverse shopping experience for visitors of all interests, ages, and income. The new schedule will also benefit vendors by giving them plenty of notice to prepare for the busy holiday season.
Holiday Sale Artisan Vendor applications can be found online at http://www.davisartcenter.org/sale.htm. You will be asked to submit pictures and descriptions of your handmade artistic items or products – we are looking for a wide variety of items, from toffee to soaps, calendars to candles, and wallart to wearableart, so don’t be shy! Applications may be submitted by email (preferred) to office@davisartcenter.org, in person at 1919 F Street in Davis, or mailed to Davis Art Center Holiday Sale, P.O. Box 4340, Davis, CA 95617. All applications will be reviewed by the Selection Committee and artists notified of their acceptance status by May 15, 2013.
Some Holiday Sale traditions, however, will not change: the sale will still take place, as always, over the first weekend in December, spanning three days of holiday shopping paradise.

Proceeds from the Davis Art Center Holiday Sale will support the Davis Art Center’s MAKE ART HAPPEN campaign, raising money to ensure that every person who walks through the Art Center doors, literally and virtually, can have a creative experience.


text:
      program list2bin
c     ================
c
c     transforms an ASCII line list to a binary
c
      implicit real*8 (a-h,o-z)
c
      n=0
      read(5,*) itype
c
c     itype: type of the line list
c          = 0 - classical atomic line list
c          = 1 - atomic line list with Carlos' extended parameters
c          = 2 - molecular line list
c
      if(itype.eq.0) then
   10    continue
         read(19,*,err=10,end=20) alam,anum,gf,excl,ql,excu,qu,gr,gs,gw
         write(11) alam,anum,gf,excl,ql,excu,qu,gr,gs,gw
         n=n+1
         go to 10
   20    continue
       else if(itype.eq.1) then
   30    continue
         read(19,*,err=30,end=40) alam,anum,gf,excl,ql,excu,qu,i
     *        gr,gs,gw,i1,i2,i3,i4,i5,i6,i7
         write(11) alam,anum,gf,excl,ql,excu,qu,gr,gs,gw,
     *        i1,i2,i3,i4,i5,i6,i7
         n=n+1
         go to 30
   40    continue
c
       else
   50    continue
         read(20,*,err=50,end=60) alam,anum,gf,excl,gr,gs,gw
         write(11) alam,anum,gf,excl,gr,gs,gw
         n=n+1
         go to 50
   60    continue
      end if
c
      write(6,601) n
  601 format(i10,'  lines included in the binary file fort.11')
      end

max_stars_repo_name: dnidever/apogee

max_issues_repo_name: dnidever/apogee

max_forks_repo_name: dnidever/apogee


text:
!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
!> \brief This module provides a subroutine to find the eigenvalues and eigenvectors
!! of a symmetric tridiagonal matrix by the ql method.
module m_fstr_EIG_tridiag
  use hecmw

  implicit none

  public

  type fstr_eigen_vec
    real(kind=kreal), pointer :: q(:) => null()
  end type fstr_eigen_vec

  type fstr_tri_diag
    real(kind=kreal), allocatable :: alpha(:)
    real(kind=kreal), allocatable :: beta(:)
  end type fstr_tri_diag

contains

  subroutine tridiag(hecMESH, hecMAT, fstrEIG, Q, Tri, iter, is_converge)
    use hecmw
    use m_fstr
    use m_fstr_EIG_lanczos_util
    implicit none
    type(hecmwST_local_mesh) :: hecMESH
    type(hecmwST_matrix)     :: hecMAT
    type(fstr_eigen)         :: fstrEIG
    type(fstr_tri_diag)      :: Tri
    type(fstr_eigen_vec), pointer :: Q(:)

    integer(kind=kint), intent(in) :: iter
    integer(kind=kint) :: N, NP, NDOF, NNDOF, NPNDOF
    integer(kind=kint) :: i, j, k, in, jn, kn, nget
    integer(kind=kint) :: iter2, ierr, maxiter
    real(kind=kreal)   :: resid, chk, sigma, tolerance, max_eigval
    real(kind=kreal), allocatable :: alpha(:), beta(:), temp(:)
    real(kind=kreal), allocatable :: L(:,:)

    integer(kind=kint), allocatable :: iparm(:)
    real(kind=kreal), pointer       :: eigvec(:,:)
    real(kind=kreal), pointer       :: eigval(:)
    logical :: is_converge

    N      = hecMAT%N
    NP     = hecMAT%NP
    NDOF   = hecMESH%n_dof
    NNDOF  = N *NDOF
    NPNDOF = NP*NDOF
    eigval => fstrEIG%eigval
    eigvec => fstrEIG%eigvec
    nget      = fstrEIG%nget
    maxiter   = fstrEIG%maxiter
    tolerance = fstrEIG%tolerance

    allocate( iparm(maxiter) )
    allocate( temp(maxiter)  )
    allocate( alpha(iter)    )
    allocate( beta(iter)     )
    allocate( L(iter, iter)  )

    do j=1, iter
      do i = 1,iter
        L(i,j)  = 0.0d0
      enddo
    enddo

    do i=1, iter
      alpha(i) = Tri%alpha(i)
      L(i,i)   = 1.0d0
    enddo

    beta(1) = 0.0d0
    do i=2, iter
      beta(i) = Tri%beta(i)
    enddo

    call QL_decomposition(iter, iter, alpha, beta, L, ierr)

    is_converge = .true.
    chk = 0.0d0
    max_eigval = maxval(alpha)
    do i = 1, min(nget, iter)
      resid = dabs(Tri%beta(iter+1)*L(iter,i))/max_eigval
      chk = max(chk, resid)
      if(tolerance < resid) is_converge = .false.
    enddo
    if(myrank == 0) write(*,"(i8,1pe12.5)")iter, chk

    if(iter < nget) is_converge = .false.
    if(iter == maxiter-1) is_converge = .true.

    if(is_converge)then
      sigma = 0.0d0
      if(fstrEIG%is_free) sigma = fstrEIG%sigma
      do i = 1, iter
        if(alpha(i) /= 0.0d0)then
          eigval(i) = 1.0d0/alpha(i) + sigma
        endif
      enddo

      call evsort(eigval, iparm, iter)

      temp = eigval

      eigvec = 0.0d0
      do k=1, iter
        in = iparm(k)
        eigval(k) = temp(in)
        do j=1, iter
          do i=1, NPNDOF
            eigvec(i, k) = eigvec(i, k) + Q(j)%q(i) * L(j, in)
          enddo
        enddo
      enddo

      do j=1, iter
        chk = maxval(eigvec(:,j))
        call hecmw_allreduce_R1(hecMESH, chk, hecmw_max)
        if(chk /= 0.0d0)then
          chk = 1.0d0/chk
          do i = 1, NNDOF
            eigvec(i,j) = eigvec(i,j) * chk
          enddo
        endif
      enddo
    endif

    deallocate(iparm)
    deallocate(temp)
    deallocate(alpha)
    deallocate(beta)
    deallocate(L)
  end subroutine tridiag

  !======================================================================!
  !                       Description                                    !
  !======================================================================!
  !>This subroutine has been adapted from the eispack routine tql2, which is
  !!a translation of the algol procedure tql2, num. math. 11, 293-306(1968)
  !!by bowdler, martin, reinsch and wilkinson.
  !!Handbook for auto. cop., vol.ii-linear algebra, 227-240(1971).
  !!This subroutine finds the eigenvalues and eigenvectors
  !!of a symmetric tridiagonal matrix by the ql method.
  !!the eigenvectors of a full symmetric matrix can also
  !!be found if  tred2  has been used to reduce this
  !!full matrix to tridiagonal form.
  !
  !on input
  !
  !nm must be set to the row dimension of two-dimensional
  !array parameters as declared in the calling program
  !dimension statement.
  !
  !is the order of the matrix.
  !
  !contains the diagonal elements of the input matrix.
  !
  !contains the subdiagonal elements of the input matrix
  !in its last n-1 positions.  e(1) is arbitrary.
  !
  !contains the transformation matrix produced in the
  !reduction by  tred2, if performed.  if the eigenvectors
  !of the tridiagonal matrix are desired, z must contain
  !the identity matrix.
  !
  !on output
  !
  !d
  !contains the eigenvalues in ascending order.  if an
  !error exit is made, the eigenvalues are correct but
  !unordered for indices 1,2,...,ierror-1.
  !-----------------------------------------------------
  !GP
  !du
  !contains the unordered eigenvalues.  if an
  !error exit is made, the eigenvalues are correct and
  !unordered for indices 1,2,...,ierror-1.
  !-----------------------------------------------------
  !e
  !has been destroyed.
  !
  !z
  !contains orthonormal eigenvectors of the symmetric
  !tridiagonal (or full) matrix.  if an error exit is made,
  !z contains the eigenvectors associated with the stored
  !eigenvalues.
  !-----------------------------------------------------
  !zu
  !contains unordered eigenvectors of the symm. tridiag. matrix
  !-----------------------------------------------------
  !ierror is set to
  !  zero       for normal return,
  !  j          if the j-th eigenvalue has not been
  !             determined after 30 iterations.
  !
  !calls a2b2 for  dsqrt(a*a + b*b) .
  !=======================================================================

  !call QL_decomposition(iter, iter, alpha, beta, L, ierr)
  subroutine QL_decomposition(nm, n, d, e, z, ierror)
    use hecmw
    implicit none
    integer(kind=kint) :: i, j, k, l, m, n, ii, l1, l2, nm, mml, ierror
    real(kind=kreal)   :: d(n), e(n), z(nm, n)
    real(kind=kreal)   :: c, c2, c3, dl1, el1, f, g, h, p, r, s, s2, tst1, tst2

    ierror = 0
    if (n .eq. 1) go to 1001

    do i = 2, n
      e(i-1) = e(i)
    enddo

    f = 0.0d0
    tst1 = 0.0d0
    e(n) = 0.0d0

    do 240 l = 1, n
      j = 0
      h = dabs(d(l)) + dabs(e(l))
      if (tst1 .lt. h) tst1 = h
      !     .......... look for small sub-diagonal element ..........
      bb:do m = l, n
        tst2 = tst1 + dabs(e(m))
        if (tst2 .eq. tst1) exit bb
        !     .......... e(n) is always zero, so there is no exit
        !                through the bottom of the loop ..........
      enddo bb

      if (m .eq. l)  go to 220

      130    if (j .eq. 30) go to 1000
      j = j + 1
      !     .......... form shift ..........
      l1 = l + 1
      l2 = l1 + 1
      g = d(l)
      p = (d(l1) - g) / (2.0d0 * e(l))
      r = a2b2(p,1.0d0)
      d(l) = e(l) / (p + dsign(r,p))
      d(l1) = e(l) * (p + dsign(r,p))
      dl1 = d(l1)
      h = g - d(l)
      if (l2 .gt. n) go to 145

      do i = l2, n
        d(i) = d(i) - h
      enddo

      145    f = f + h
      !     .......... ql transformation ..........
      p = d(m)
      c = 1.0d0
      c2 = c
      el1 = e(l1)
      s = 0.0d0
      s2 = 0.0d0
      c3 = 0.0d0
      mml = m - l
      !     .......... for i=m-1 step -1 until l do -- ..........
      do ii = 1, mml
        c3 = c2
        c2 = c
        s2 = s
        i = m - ii
        g = c * e(i)
        h = c * p
        r = a2b2(p,e(i))
        e(i+1) = s * r
        s = e(i) / r
        c = p / r
        p = c * d(i) - s * g
        d(i+1) = h + s * (c * g + s * d(i))
        !     .......... form vector ..........
        do k = 1, n
          h = z(k,i+1)
          z(k,i+1) = s * z(k,i) + c * h
          z(k,i) = c * z(k,i) - s * h
        enddo
      enddo

      p = -s * s2 * c3 * el1 * e(l) / dl1
      e(l) = s * p
      d(l) = c * p
      tst2 = tst1 + dabs(e(l))
      if (tst2 .gt. tst1) go to 130
      220    d(l) = d(l) + f
      240 continue

      !     .......... order eigenvalues and eigenvectors ..........
      !GP: Get unordered eigenvalues and eigenvectors----------------

      do 300 ii = 2, n
        i = ii - 1
        k = i
        p = d(i)

        aa:do j = ii, n
          if (d(j) .ge. p) exit aa
          k = j
          p = d(j)
        enddo aa

        if (k .eq. i) go to 300
        d(k) = d(i)
        d(i) = p

        do j = 1, n
          p = z(j,i)
          z(j,i) = z(j,k)
          z(j,k) = p
        enddo

        300 continue

        go to 1001
        !     .......... set error -- no convergence to an
        !                eigenvalue after 30 iterations ..........
        1000 ierror = l
        1001 return
  end subroutine QL_decomposition

  function a2b2(a,b)
    use hecmw
    implicit none

    real(kind=kreal) :: a2b2
    real(kind=kreal) :: a, b
    real(kind=kreal) :: p, q, r, s, t, u

    p = dmax1(dabs(a), dabs(b))
    if (p /= 0.0d0) then
      r = (dmin1(dabs(a),dabs(b))/p) ** 2
      do
        t = 4.0d0 + r
        if (t == 4.0d0) exit
        s = r/t
        u = 1.0d0 + 2.0d0*s
        p = u * p
        q = s/u
        r = q * q * r
      end do
    end if
    a2b2 = p
    return

  end function a2b2

end module m_fstr_EIG_tridiag

max_stars_repo_name: FrontISTR/FrontISTR

max_issues_repo_name: FrontISTR/FrontISTR

max_forks_repo_name: FrontISTR/FrontISTR


text:
     TYPE CELL
        INTEGER :: VAL
        TYPE (CELL), POINTER :: NEXT_CELL
     END TYPE CELL

! is defined and objects such as HEAD and CURRENT are declared using

     TYPE (CELL), TARGET :: HEAD
     TYPE (CELL), POINTER :: CURRENT

! If a linked list has been created and attached to HEAD and the pointer CURRENT has been allocated
! space, statements such as

     CURRENT = HEAD
     CURRENT = CURRENT % NEXT_CELL

! cause the contents of the cells referenced on the right to be copied to the cell referred to by CURRENT.
! In particular, the right-hand side of the second statement causes the pointer component in the cell,
! CURRENT, to be selected. This pointer is dereferenced because it is in an expression context to produce
! the target's integer value and a pointer to a cell that is in the target's NEXT CELL component. The
! left-hand side causes the pointer CURRENT to be dereferenced to produce its present target, namely
! space to hold a cell (an integer and a cell pointer). The integer value on the right is copied to the integer
! space on the left and the pointer component is pointer assigned (the descriptor on the right is copied
! into the space for a descriptor on the left). When a statement such as

     CURRENT => CURRENT % NEXT_CELL

! is executed, the descriptor value in CURRENT % NEXT CELL is copied to the descriptor named
! CURRENT. In this case, CURRENT is made to point at a different target.
! In the intrinsic assignment statement, the space associated with the current pointer does not change but
! the values stored in that space do. In the pointer assignment, the current pointer is made to associate
! with different space. Using the intrinsic assignment causes a linked list of cells to be moved up through
! the current "window"; the pointer assignment causes the current pointer to be moved down through the
! list of cells.

end

max_stars_repo_name: OpenFortranProject/ofp-sdf

max_issues_repo_name: OpenFortranProject/ofp-sdf

max_forks_repo_name: OpenFortranProject/ofp-sdf


text:
This is the place to go if you want to have a BBQ with all your friends that live in the nearby apartment complexes. It is part of the West Manor neighborhood. It appears that the citys open container ordinance does not apply to this park.

  Im told that this park looks the way it does because it was designed by Landscape Architecture students.  Users/NickSchmalenberger

Also see Parks for info on other local parks


!
! Copyright (C) 2002-2009 Quantum ESPRESSO group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!---------------------------------------------------------------------
function qe_erf (x)
  !---------------------------------------------------------------------
  !
  !     Error function - computed from the rational approximations of
  !     W. J. Cody, Math. Comp. 22 (1969), pages 631-637.
  !
  !     for abs(x) le 0.47 erf is calculated directly
  !     for abs(x) gt 0.47 erf is calculated via erf(x)=1-erfc(x)
  !
  IMPLICIT NONE
  INTEGER, PARAMETER :: DP=8
  REAL(DP), INTENT(IN) :: x
  REAL(DP) :: x2, p1 (4), q1 (4)
  REAL(DP), EXTERNAL :: qe_erfc
  REAL(DP) :: qe_erf
  DATA p1 / 2.426679552305318E2_DP, 2.197926161829415E1_DP, &
            6.996383488619136_DP,  -3.560984370181538E-2_DP /
  DATA q1 / 2.150588758698612E2_DP, 9.116490540451490E1_DP, &
            1.508279763040779E1_DP, 1.000000000000000_DP /
  !
  IF (abs (x) > 6.0_DP) THEN
     !
     !  erf(6)=1-10^(-17) cannot be distinguished from 1
     !
     qe_erf = sign (1.0_DP, x)
  ELSE
     IF (abs (x)  <= 0.47_DP) THEN
        x2 = x**2
        qe_erf=x *(p1 (1) + x2 * (p1 (2) + x2 * (p1 (3) + x2 * p1 (4) ) ) ) &
                / (q1 (1) + x2 * (q1 (2) + x2 * (q1 (3) + x2 * q1 (4) ) ) )
     ELSE
        qe_erf = 1.0_DP - qe_erfc (x)
     ENDIF
  ENDIF
  !
  RETURN
END FUNCTION qe_erf
!
!---------------------------------------------------------------------
FUNCTION qe_erfc (x)
  !---------------------------------------------------------------------
  !
  !     erfc(x) = 1-erf(x)  - See comments in erf
  !
  IMPLICIT NONE
  INTEGER, PARAMETER :: DP=8
  REAL(DP), INTENT(IN) :: x
  REAL(DP) :: qe_erfc
  REAL(DP) :: ax, x2, xm2, p2 (8), q2 (8), p3 (5), q3 (5), pim1
  REAL(DP), EXTERNAL :: qe_erf
  DATA p2 / 3.004592610201616E2_DP,  4.519189537118719E2_DP, &
            3.393208167343437E2_DP,  1.529892850469404E2_DP, &
            4.316222722205674E1_DP,  7.211758250883094_DP,   &
            5.641955174789740E-1_DP,-1.368648573827167E-7_DP /
  DATA q2 / 3.004592609569833E2_DP,  7.909509253278980E2_DP, &
            9.313540948506096E2_DP,  6.389802644656312E2_DP, &
            2.775854447439876E2_DP,  7.700015293522947E1_DP, &
            1.278272731962942E1_DP,  1.000000000000000_DP /
  DATA p3 /-2.996107077035422E-3_DP,-4.947309106232507E-2_DP, &
           -2.269565935396869E-1_DP,-2.786613086096478E-1_DP, &
           -2.231924597341847E-2_DP /
  DATA q3 / 1.062092305284679E-2_DP, 1.913089261078298E-1_DP, &
            1.051675107067932_DP,    1.987332018171353_DP,    &
            1.000000000000000_DP /

  DATA pim1 / 0.56418958354775629_DP /
  !        ( pim1= sqrt(1/pi) )
  ax = abs (x)
  if (ax > 26.0_DP) then
     !
     !  erfc(26.0)=10^(-296); erfc( 9.0)=10^(-37);
     !
     qe_erfc = 0.0_DP
  elseif (ax > 4.0_DP) then
     x2 = x**2
     xm2 = (1.0_DP / ax) **2
     qe_erfc = (1.0_DP / ax) * exp ( - x2) * (pim1 + xm2 * (p3 (1) &
          + xm2 * (p3 (2) + xm2 * (p3 (3) + xm2 * (p3 (4) + xm2 * p3 (5) &
          ) ) ) ) / (q3 (1) + xm2 * (q3 (2) + xm2 * (q3 (3) + xm2 * &
          (q3 (4) + xm2 * q3 (5) ) ) ) ) )
  elseif (ax > 0.47_DP) then
     x2 = x**2
     qe_erfc = exp ( - x2) * (p2 (1) + ax * (p2 (2) + ax * (p2 (3) &
          + ax * (p2 (4) + ax * (p2 (5) + ax * (p2 (6) + ax * (p2 (7) &
          + ax * p2 (8) ) ) ) ) ) ) ) / (q2 (1) + ax * (q2 (2) + ax * &
          (q2 (3) + ax * (q2 (4) + ax * (q2 (5) + ax * (q2 (6) + ax * &
          (q2 (7) + ax * q2 (8) ) ) ) ) ) ) )
  else
     qe_erfc = 1.0_DP - qe_erf (ax)
  endif
  !
  ! erf(-x)=-erf(x)  =>  erfc(-x) = 2-erfc(x)
  !
  if (x < 0.0_DP) qe_erfc = 2.0_DP - qe_erfc
  !
  return
end function qe_erfc
!
!---------------------------------------------------------------------
FUNCTION gauss_freq (x)
  !---------------------------------------------------------------------
  !
  !     gauss_freq(x) = (1+erf(x/sqrt(2)))/2 = erfc(-x/sqrt(2))/2
  !             - See comments in erf
  !
  IMPLICIT NONE
  INTEGER, PARAMETER :: DP=8
  REAL(DP), INTENT(IN) :: x
  REAL(DP)            :: gauss_freq
  REAL(DP), PARAMETER :: c =  0.7071067811865475_DP
  !        ( c= sqrt(1/2) )
  REAL(DP), EXTERNAL :: qe_erfc
  !
  gauss_freq = 0.5_DP * qe_erfc ( - x * c)
  !
  return
end function gauss_freq


text:
! libRQZ
! author: daan. camps@cs.kuleuven.be
! Description:
!   Deflation monitoring
! ___________________________________________________________________
module d_deflations

use u_activeparts
use u_parameters
use d_ctransformations
use d_memorymgmt
use d_computepoles

implicit none
private
public d_check_interior_deflations, d_start_deflation_single, &
       d_start_deflation_double, d_stop_deflation_single, &
       d_stop_deflation_double
contains

  subroutine d_check_interior_deflations(defl,A, B, compAB, apC, apT, &
              Q, compQ, qcol, Z, compZ, zcol)
  ! DESCRIPTION
  !___________________________________________________________________________
  ! Check for interior deflations in active part defined by apC.
  ! When deflated poles are encountered, the pencil is updated,
  ! i.e., these poles are explicitly set to zero, and the active
  ! parts list apC is updated.The boolean defl returns true when
  ! at least 1 deflation is found.
  !
  ! ARGUMENTS
  !___________________________________________________________________________
  ! defl    boolean [OUT]
  !           Deflation indicator
  ! A       double array [INOUT]
  !           Upper block-Hessenberg matrix A
  ! B       double array [INOUT]
  !           Upper block-Hessenberg matrix B
  ! compAB  boolean [IN]
  !           If true, the equivalences are applied to (A,B). If not,
  !           then only to the active parts
  ! apC     type(tAp) [INOUT]
  !           Active parts linked list.
  ! apT     type(tAp) [IN]
  !           Total part of the pencil the method is acting on
  !           (used when compAB is true)
  !           Normally, the following order,
  !             tstrt <= cstrt < i < i+1 < cstp <= tstp,
  !           holds when calling d_swap. c = current, t = total
  ! Q       double array [INOUT]
  !           Unitary equivalence matrix Q
  ! compQ   boolean [IN]
  !           If true, the left Schur vectors are updated
  ! qcol    integer [IN]
  !           Start index of the columns of Q that are updated:
  !              qcol, qcol+1
  ! Z       double array [INOUT]
  !           Unitary equivalence matrix Z
  ! compZ   boolean [IN]
  !           If true, the right Schur vectors are updated
  ! zcol    integer [IN]
  !           Start index of the columns of Z that are updated:
  !              zcol, zcol+1
  !___________________________________________________________________________
  ! last edit: January 11, 2018
    real(kind=dp), intent(inout)      :: A(:,:), B(:,:), Q(:,:), Z(:,:)
    type(tAp), pointer, intent(inout) :: apC
    type(tAp), pointer, intent(in)    :: apT
    logical, intent(out)              :: defl
    logical, intent(in)               :: compAB, compQ, compZ
    integer, intent(in)               :: qcol, zcol

    integer                           :: i, j, cpy, cpystrt
    integer, pointer                  :: cstrt, cstp, tstrt, tstp
    real(kind=dp)                     :: l1, l2
    logical                           :: isreal, defli
    call tApGet(apC,cstrt,cstp)
    call tApGet(apT,tstrt,tstp)

    defl = .false.
    i = cstrt+1
    cpystrt = cstrt
    do while (i .le. cstp-2)
      if (abs(A(i+2,i)) .gt. dzero) then
        ! double pole
        call d_eigenvalues2x2(A,B,i+1,i,l1,l2,isreal)
        if (isreal) then
          ! Not C.C.
          call d_make_hess(l1,i+1,i,A,B,compAB,apC,apT,Q,compQ,qcol+i-cpystrt,&
                 Z,compZ,zcol+i-cpystrt-1)
          call d_interior_deflation_sng(i,defli,A,B)
          if (defli) then
            defl = .true.
            ! update active parts
            cpy = cstp
            call tApPut(apC,strt=cstrt, stp=i)
            call tApInsertAfter(apC,strt=i,stp=cpy)
            apC => tApNext(apC)
            call tApGet(apC,cstrt,cstp) ! this does update cstrt,cstp
                                        !(but does not affect the loop counter)
          end if
          i = i + 1
        else
          j = i
          call d_interior_deflation_dbl(j,defli,A,B)
          if (defli) then
            defl = .true.
            ! update active parts
            cpy = cstp
            call tApPut(apC,strt=cstrt, stp=i+j)
            call tApInsertAfter(apC,strt=i+j,stp=cpy)
            apC => tApNext(apC)
            call tApGet(apC,cstrt,cstp) ! this does update cstrt,cstp
                                        !(but does not affect the loop counter)
          end if
          i = i + 2
        end if
      else
        call d_interior_deflation_sng(i,defli,A,B)
        if (defli) then
          defl = .true.
          ! update active parts
          cpy = cstp
          call tApPut(apC,strt=cstrt, stp=i)
          call tApInsertAfter(apC,strt=i,stp=cpy)
          apC => tApNext(apC)
          call tApGet(apC,cstrt,cstp) ! this does update cstrt,cstp
                                      !(but does not affect the loop counter)
        end if
        i = i + 1
      end if
    end do
  end subroutine d_check_interior_deflations

  subroutine d_interior_deflation_sng(i,defl,A,B)
    ! DESCRIPTION
    !___________________________________________________________________________
    ! Test if the single interior pole at position (i+1,i) in (A,B)
    ! is deflated
    !
    ! ARGUMENTS
    !___________________________________________________________________________
    ! i       integer [IN]
    !           column index of single pole
    ! defl    boolean [INOUT]
    !           returns true of the ith pole can be deflated.
    !           If so, the respective entries in A,B are set to zero.
    ! A       double array [INOUT]
    !           block-Hessenberg matrix
    ! B       double array [INOUT]
    !           Hessenberg matrix
    !
    ! NOTA
    !___________________________________________________________________________
    ! This is a private subroutine called from  d_check_interior_deflations
    !___________________________________________________________________________
    ! last edit: January 11, 2019
    real(kind=dp), intent(inout)      :: A(:,:), B(:,:)
    integer, intent(in)               :: i
    logical, intent(inout)            :: defl

    defl = .false.
    if ((abs(A(i+1,i)) .lt. mp * (abs(A(i,i)) + abs(A(i+1,i+1)))) .and. &
      (abs(B(i+1,i)) .lt. mp * (abs(B(i,i)) + abs(B(i+1,i+1))))) then
      ! Set to zero
      defl = .true.
      A(i+1,i) = dzero
      B(i+1,i) = dzero
    end if
  end subroutine

  subroutine d_interior_deflation_dbl(i,defl,A,B)
  ! DESCRIPTION
  !___________________________________________________________________________
  ! Test if the double interior pole at position (i+1:i+2,i:i+1)
  ! in (A,B) is deflated
  ! ARGUMENTS
  !___________________________________________________________________________
  ! i       integer [IN]
  !           column index of double pole
  ! defl    boolean [INOUT]
  !           returns true of the ith pole can be deflated.
  !           If so, the respective entries in A,B are set to zero.
  ! A       double array [INOUT]
  !           block-Hessenberg matrix
  ! B       double array [INOUT]
  !           Hessenberg matrix
  !
  ! NOTA
  !___________________________________________________________________________
  ! This is a private subroutine called from  d_check_interior_deflations
  !___________________________________________________________________________
  ! last edit: February 7, 2019
    real(kind=dp), intent(inout)      :: A(:,:), B(:,:)
    logical, intent(inout)            :: defl
    integer, intent(inout)            :: i

    defl = .false.
    if (((abs(A(i+2,i)) + abs(A(i+2,i+1))) .lt. mp * &
        (abs(A(i+1,i+1)) + abs(A(i+2,i+2)))) .and. &
        (abs(B(i+2,i+1)) .lt. mp * &
        (abs(B(i+1,i+1)) + abs(B(i+2,i+2))))) then
        defl = .true.
        B(i+2,i:i+1) = dzero
        A(i+2,i:i+1) = dzero
        i = 1 ! Ensures that active parts are handled in calling function
    elseif (((abs(A(i+1,i)) + abs(A(i+2,i))) .lt. mp * &
        (abs(A(i,i)) + abs(A(i+1,i+1)))) .and. &
        (abs(B(i+1,i)) .lt. mp * &
        (abs(B(i,i)) + abs(B(i+1,i+1))))) then
        defl = .true.
        B(i+1:i+2,i) = dzero
        A(i+1:i+2,i) = dzero
        i = 0 ! Ensures that active parts are handled in calling function
      end if
  end subroutine

  subroutine d_start_deflation_single(defl, A, B, compAB, apC, apT, &
                                     Q, compQ, qcol)
  ! DESCRIPTION
  !___________________________________________________________________________
  ! Check for deflation at the start of the active area
  ! The pencil (A,B) and the Q matrix is updated as requested if a deflation
  ! is found.
  ! The current active part apC is also updated, which also updates
  ! it in the caller function.
  !
  ! ARGUMENTS
  !___________________________________________________________________________
  ! defl    boolean [OUT]
  !           Deflation indicator (true if deflated, else false)
  ! A       complex double array [INOUT]
  !           Upper Hessenberg matrix A
  ! B       complex double array [INOUT]
  !           Upper Hessenberg matrix B
  ! compAB  boolean [INOUT]
  !           If true, the equivalences are applied to (A,B) as defined by
  !           apT. If not, then only to the active part apC
  ! apC     type(tAp) [INOUT]
  !           Current active part
  ! apT     type(tAp) [IN]
  !           Total part of the pencil the method is acting on
  !           (used when compAB is true)
  ! Q       complex double array [INOUT]
  !           Unitary equivalence matrix Q
  ! compQ   boolean [IN]
  !           If true, the left Schur vectors are updated
  ! qcol    integer [IN]
  !           Start index of the columns of Q that are updated:
  !              qcol, qcol+1
  !___________________________________________________________________________
  ! last edit: September 12, 2018
    real(kind=dp), intent(inout)      ::  A(:,:), B(:,:), Q(:,:)
    logical, intent(in)               ::  compAB, compQ
    integer, intent(in)               ::  qcol
    logical, intent(out)              ::  defl
    type(tAp), pointer, intent(inout) ::  apC
    type(tAp), pointer, intent(in)    ::  apT

    real(kind=dp)                     :: c,s,r,t, tt, u(2)
    integer,pointer                   :: cstrt, cstp, tstrt, tstp
    integer                           :: LDQ

    call d_getLeadingDim(LDQ,Q)

    ! Get current and total ranges
    call tApGet(apC,cstrt,cstp)
    call tApGet(apT,tstrt,tstp)

    defl = .false.
    ! first check for direct deflation
    if ((abs(A(cstrt+2,cstrt+1)) .lt. &
         mp*(abs(A(cstrt+1,cstrt+1)) + abs(A(cstrt+2,cstrt+2)))) .and. &
        (abs(B(cstrt+2,cstrt+1)) .lt. &
         mp*(abs(B(cstrt+1,cstrt+1)) + abs(B(cstrt+2,cstrt+2))))) then
        defl = .true.
        A(cstrt+2,cstrt+1) = dzero
        B(cstrt+2,cstrt+1) = dzero
        ! This also updates it in the caller function
        call tApPut(apC,strt=cstrt+1,stp=cstp)
    elseif (abs(A(cstrt+1,cstrt+1) * B(cstrt+2,cstrt+1) - &
                A(cstrt+2,cstrt+1) * B(cstrt+1,cstrt+1)) .lt. &
           mp*(abs(A(cstrt+2,cstrt+2)) + abs(B(cstrt+2,cstrt+2)))) then
      if (abs(A(cstrt+2,cstrt+1)) .ge. abs(B(cstrt+2,cstrt+1))) then
        ! Compute rotation based on A
        call d_compute_ct(A(cstrt+1,cstrt+1),A(cstrt+2,cstrt+1),c,s,r)
        ! Apply to B
        r = B(cstrt+1,cstrt+1)
        t = B(cstrt+2,cstrt+1)
        u = B(cstrt+1:cstrt+2,cstrt+2)
      else
        ! Compute rotation based on B
        call d_compute_ct(B(cstrt+1,cstrt+1),B(cstrt+2,cstrt+1),c,s,r)
        ! Apply to A
        r = A(cstrt+1,cstrt+1)
        t = A(cstrt+2,cstrt+1)
        u = A(cstrt+1:cstrt+2,cstrt+2)
      end if

      tt = -s * r + c * t
      r = c * r + s * t
      t = tt
      u(2) = -s * u(1) + c * u(2)

      if (abs(t) .lt. mp*(abs(r) + abs(u(2)))) then
        defl = .true.
        if (compAB) then
          call d_apply_ct_l(A(cstrt+1,cstrt+1:tstp),A(cstrt+2,cstrt+1:tstp),c,s)
          call d_apply_ct_l(B(cstrt+1,cstrt+1:tstp),B(cstrt+2,cstrt+1:tstp),c,s)
        else
          call d_apply_ct_l(A(cstrt+1,cstrt+1:cstp),A(cstrt+2,cstrt+1:cstp),c,s)
          call d_apply_ct_l(B(cstrt+1,cstrt+1:cstp),B(cstrt+2,cstrt+1:cstp),c,s)
        end if
        if (compQ) then
          call d_apply_ct_r(Q(1:LDQ,qcol),Q(1:LDQ,qcol+1),c,-s)
        end if
        A(cstrt+2,cstrt+1) = dzero
        B(cstrt+2,cstrt+1) = dzero
        ! This also updates it in the caller function
        call tApPut(apC,strt=cstrt+1,stp=cstp)
      end if
    end if

    ! Double check if a 2x2 block has separated
    if (cstp-cstrt .gt. 3) then
      if (.not. (defl .or. (abs(A(cstrt+4,cstrt+2)) .gt. dzero))) then
      call d_interior_deflation_sng(cstrt+2,defl,A,B)
      if (defl) then
        ! update active parts
        call tApInsertAfter(apC,strt=cstrt+2,stp=cstp)
        call tApPut(apC,strt=cstrt, stp=cstrt+2)
      end if
     end if
    end if
  end subroutine d_start_deflation_single

  subroutine d_start_deflation_double(defl, A, B, compAB, apC, apT, &
                                     Q, compQ, qcol, Z, compZ, zcol)
  ! DESCRIPTION
  !___________________________________________________________________________
  ! Check for deflation at the start of the active area. It is assumed
  ! that there is a double pole at the start. It is also checked if the double
  ! pole actually is a complex conjugate pair. If not, the Hessenberg form
  ! is restored and a single deflation is checked.
  ! The pencil (A,B) and the Q and Z matrices are updated as requested
  ! if a deflation is found.
  ! The current active part apC is also updated, which also updates
  ! it in the caller function.
  !
  ! ARGUMENTS
  !___________________________________________________________________________
  ! defl    boolean [OUT]
  !           Deflation indicator (true if deflated, else false)
  ! A       complex double array [INOUT]
  !           Upper Hessenberg matrix A
  ! B       complex double array [INOUT]
  !           Upper Hessenberg matrix B
  ! compAB  boolean [INOUT]
  !           If true, the equivalences are applied to (A,B) as defined by
  !           apT. If not, then only to the active part apC
  ! apC     type(tAp) [INOUT]
  !           Current active part
  ! apT     type(tAp) [IN]
  !           Total part of the pencil the method is acting on
  !           (used when compAB is true)
  ! Q       complex double array [INOUT]
  !           Unitary equivalence matrix Q
  ! compQ   boolean [IN]
  !           If true, the left Schur vectors are updated
  ! qcol    integer [IN]
  !           Start index of the columns of Q that are updated:
  !              qcol, qcol+1, qcol+2
  ! Z       double array [INOUT]
  !           Orthonormal equivalence matrix Z
  ! compZ   boolean [IN]
  !           If true, the right Schur vectors are updated
  ! zcol    integer [IN]
  !           Start index of the columns of Z that are updated:
  !              zcol, zcol+1
  !___________________________________________________________________________
  ! last edit: September 12, 2018
  real(kind=dp), intent(inout)      ::  A(:,:), B(:,:), Q(:,:), Z(:,:)
  logical, intent(in)               ::  compAB, compQ, compZ
  integer, intent(in)               ::  qcol, zcol
  logical, intent(out)              ::  defl
  type(tAp), pointer, intent(inout) ::  apC
  type(tAp), pointer, intent(in)    ::  apT

  real(kind=dp)                     :: c,c2,s,s2,r,t
  logical                           :: isreal
  integer,pointer                   :: cstrt, cstp, tstrt, tstp
  integer                           :: LDAB, LDR, LDC, LDQ
  external                          :: DGEMM

  call d_getLeadingDim(LDAB,A)
  call d_getLeadingDim(LDR,Ru)
  call d_getLeadingDim(LDC,Cu)
  call d_getLeadingDim(LDQ,Q)

  ! Get current and total ranges
  call tApGet(apC,cstrt,cstp)
  call tApGet(apT,tstrt,tstp)

  defl = .false.

  ! Check first if it actually a cc pole
  call d_eigenvalues2x2(A,B,cstrt+2,cstrt+1,c,c2,isreal)
  if (isreal) then
    call d_make_hess(c,cstrt+2,cstrt+1,A,B,compAB,apC,apT,Q,compQ,qcol+1,&
                       Z,compZ,zcol)
    call d_start_deflation_single(defl,A,B,compAB,apC,apT,Q,compQ,qcol)
  else
    ! Compute the QR factorization
    ! We create a copy of the leading block for the deflation test
    Q4(1:3,1:2) = A(cstrt+1:cstrt+3,cstrt+1:cstrt+2)
    Z4(1:3,1:2) = B(cstrt+1:cstrt+3,cstrt+1:cstrt+2)
    ! Create zero in B(cstrt+2,cstrt+1)
    call d_compute_ct(Z4(1,1),Z4(2,1),c,s,r)
    ! Apply transformation
    Z4(1,1) = r
    Z4(2,1) = dzero
    t = -s * Z4(1,2) + c * Z4(2,2)
    Z4(1,2) = c * Z4(1,2) + s * Z4(2,2)
    Z4(2,2) = t
    call d_apply_ct_l(Q4(1,1:2),Q4(2,1:2),c,s)
    ! Try to zero row cstrt+3 all at once
    if ((abs(Z4(3,2)) .gt. abs(Q4(3,1))) .and. &
       (abs(Z4(3,2)) .gt. abs(Q4(3,2)))) then
      call d_compute_ct(Z4(2,2),Z4(3,2),c2,s2,r)
      Z4(2,2) = r
      Z4(3,2) = dzero
    else
      if (abs(Q4(3,1)) .gt. abs(Q4(3,2))) then
        call d_compute_ct(Q4(2,1),Q4(3,1),c2,s2,r)
      else
        call d_compute_ct(Q4(2,2),Q4(3,2),c2,s2,r)
      end if
      t = -s2 * Z4(2,2) + c2 * Z4(3,2)
      Z4(2,2) = c2 * Z4(2,2) + s2 * Z4(3,2)
      Z4(3,2) = t
    end if
    call d_apply_ct_l(Q4(2,1:2),Q4(3,1:2),c2,s2)
    ! Test if deflated
    if (((abs(Q4(3,1)) + abs(Q4(3,2))) .lt. &
          mp * (abs(Q4(1,1)) + abs(Q4(2,2)) + abs(Q4(1,2)) &
          + abs(Q4(2,1)))) .and. (abs(Z4(3,2)) .lt. &
           mp * (abs(Z4(1,1)) + abs(Z4(2,2)) + abs(Z4(1,2))))) then
        defl = .true.
        ! Construct Q
        Q4 = dzero
        Q4(1,1) = done
        Q4(2,2) = done
        Q4(3,3) = done
        call d_apply_ct_l(Q4(1,1:3),Q4(2,1:3),c,s)
        call d_apply_ct_l(Q4(2,1:3),Q4(3,1:3),c2,s2)
    end if

    if (defl) then
      if (compAB) then
        call DGEMM('N','N',3,tstp-cstrt,3,done,Q4(1,1),4,A(cstrt+1,cstrt+1),&
                       LDAB,dzero,Ru(1,1),LDR)
        A(cstrt+1:cstrt+3,cstrt+1:tstp) = Ru(1:3,1:tstp-cstrt)
        call DGEMM('N','N',3,tstp-cstrt,3,done,Q4(1,1),4,B(cstrt+1,cstrt+1),&
                       LDAB,dzero,Ru(1,1),LDR)
        B(cstrt+1:cstrt+3,cstrt+1:tstp) = Ru(1:3,1:tstp-cstrt)
      else
        call DGEMM('N','N',3,cstp-cstrt,3,done,Q4(1,1),4,A(cstrt+1,cstrt+1),&
                       LDAB,dzero,Ru(1,1),LDR)
        A(cstrt+1:cstrt+3,cstrt+1:cstp) = Ru(1:3,1:cstp-cstrt)
        call DGEMM('N','N',3,cstp-cstrt,3,done,Q4(1,1),4,B(cstrt+1,cstrt+1),&
                       LDAB,dzero,Ru(1,1),LDR)
        B(cstrt+1:cstrt+3,cstrt+1:cstp) = Ru(1:3,1:cstp-cstrt)
      end if
      ! Put elements to zero
      A(cstrt+3,cstrt+1:cstrt+2) = dzero
      B(cstrt+3,cstrt+1:cstrt+2) = dzero
      B(cstrt+2,cstrt+1) = dzero
      if (compQ) then
        call DGEMM('N','T',LDQ,3,3,cone,Q(1,qcol),LDQ,Q4(1,1),4,&
                      czero,Cu(1,1),LDC)
        Q(1:LDQ,qcol:qcol+2) = Cu(1:LDQ,1:3)
      end if
      ! This also updates it in the caller function
      call tApPut(apC,strt=cstrt+2,stp=cstp)
    end if
  end if
  end subroutine d_start_deflation_double

  subroutine d_stop_deflation_single(defl, A, B, compAB, apC, apT,&
                                    Z, compZ, zcol )
  !___________________________________________________________________________
  ! DESCRIPTION
  !___________________________________________________________________________
  ! Check for deflation at the end of the active area
  ! The pencil (A,B) and the Z matrix is updated as requested if a deflation
  ! is found.
  ! The current active part apC is also updated, which also updates
  ! it in the caller function.
  !
  ! ARGUMENTS
  !___________________________________________________________________________
  ! defl    boolean [OUT]
  !           Deflation indicator (true if deflated, else false)
  ! A       complex double array [INOUT]
  !           Upper Hessenberg matrix A
  ! B       complex double array [INOUT]
  !           Upper Hessenberg matrix B
  ! compAB  boolean [INOUT]
  !           If true, the equivalences are applied to (A,B). If not,
  !           then only to the current active part
  ! apC     type(tAp) [INOUT]
  !           Current active part
  ! apT     type(tAp) [IN]
  !           Total part of the pencil the method is acting on
  !           (used when compAB is true)
  ! Z       complex double array [INOUT]
  !           Unitary equivalence matrix Z
  ! compZ   boolean [IN]
  !           If true, the right Schur vectors are updated
  ! zcol    integer [IN]
  !           Start index of the columns of Z that are updated:
  !              zcol, zcol+1
  !___________________________________________________________________________
  ! last edit: September 12, 2018
    real(kind=dp), intent(inout)      ::  A(:,:), B(:,:), Z(:,:)
    logical, intent(in)               ::  compAB, compZ
    integer, intent(in)               ::  zcol
    logical, intent(out)              ::  defl
    type(tAp), pointer, intent(inout) ::  apC
    type(tAp), pointer, intent(in)    ::  apT

    real(kind=dp)                     :: c, s, r,t, tt, u(2)
    integer,pointer                   :: cstrt, cstp, tstrt, tstp
    integer                           :: LDZ

    call d_getLeadingDim(LDZ,Z)

    call tApGet(apC,cstrt,cstp)
    call tApGet(apT,tstrt,tstp)

    defl = .false.
    ! first check for direct deflation
    if ((abs(A(cstp,cstp-1)) .lt. &
         mp*(abs(A(cstp,cstp)) + abs(A(cstp-1,cstp-1)))) .and. &
        (abs(B(cstp,cstp-1)) .lt. &
         mp*(abs(B(cstp,cstp)) + abs(B(cstp-1,cstp-1))))) then
      defl = .true.
      A(cstp,cstp-1) = dzero
      B(cstp,cstp-1) = dzero
      ! This also updates it in the caller function
      call tApPut(apC,strt=cstrt,stp=cstp-1)
    elseif (abs(A(cstp,cstp) * B(cstp,cstp-1) - &
                A(cstp,cstp-1) * B(cstp,cstp)) .lt. &
           mp*(abs(A(cstp-1,cstp-1)) + abs(B(cstp-1,cstp-1)))) then
      if (abs(A(cstp,cstp-1)) .ge. abs(B(cstp,cstp-1))) then
        ! Compute rotation based on A
        call d_compute_ct(A(cstp,cstp),A(cstp,cstp-1),c,s,r)
        ! Apply to B
        r = B(cstp,cstp)
        t = B(cstp,cstp-1)
        u = B(cstp-1,cstp-1:cstp)
      else
        ! Compute rotation based on B
        call d_compute_ct(B(cstp,cstp),B(cstp,cstp-1),c,s,r)
        ! Apply to A
        r = A(cstp,cstp)
        t = A(cstp,cstp-1)
        u = A(cstp-1,cstp-1:cstp)
      end if

      tt = -s * r + c * t
      r = c * r + s * t
      t = tt
      u(2) =  -s * u(1) + c * u(2)

      if (abs(t) .lt. mp*(abs(r) + abs(u(2)))) then
        defl = .true.
        if (compAB) then
          call d_apply_ct_r(A(tstrt+1:cstp,cstp-1),A(tstrt+1:cstp,cstp),c,s)
          call d_apply_ct_r(B(tstrt+1:cstp,cstp-1),B(tstrt+1:cstp,cstp),c,s)
        else
          call d_apply_ct_r(A(cstrt+1:cstp,cstp-1),A(cstrt+1:cstp,cstp),c,s)
          call d_apply_ct_r(B(cstrt+1:cstp,cstp-1),B(cstrt+1:cstp,cstp),c,s)
        endif
        if (compZ) then
          call d_apply_ct_r(Z(1:LDZ,zcol),Z(1:LDZ,zcol+1),c,s)
        end if
        A(cstp,cstp-1) = dzero
        B(cstp,cstp-1) = dzero
        ! This also updates it in the caller function
        call tApPut(apC,strt=cstrt,stp=cstp-1)
      end if
    end if

    ! Double check if a 2x2 block has separated
    if (cstp-cstrt .gt. 3) then
      if (.not. (defl .or. (abs(A(cstp-1,cstp-3)) .gt. dzero))) then
      call d_interior_deflation_sng(cstp-2,defl,A,B)
      if (defl) then
        ! update active parts
        call tApInsertAfter(apC,strt=cstp-2,stp=cstp)
        call tApPut(apC,strt=cstrt, stp=cstp-2)
      end if
     end if
    end if
  end subroutine d_stop_deflation_single

  subroutine d_stop_deflation_double(defl, A, B, compAB, apC, apT,&
                                    Q, compQ, qcol, Z, compZ, zcol )
  ! DESCRIPTION
  !___________________________________________________________________________
  ! Check for deflation at the end of the active area. It is assumed
  ! that there is a double pole at the end. It is also checked if the double
  ! pole actually is a complex conjugate pair. If not, the Hessenberg form
  ! is restored and a single deflation is checked.
  ! The pencil (A,B) and the Q and Z matrices are updated as requested
  ! if a deflation is found.
  ! The current active part apC is also updated, which also updates
  ! it in the caller function.
  !
  ! ARGUMENTS
  !___________________________________________________________________________
  ! defl    boolean [OUT]
  !           Deflation indicator (true if deflated, else false)
  ! A       complex double array [INOUT]
  !           Upper Hessenberg matrix A
  ! B       complex double array [INOUT]
  !           Upper Hessenberg matrix B
  ! compAB  boolean [INOUT]
  !           If true, the equivalences are applied to (A,B) as defined by
  !           apT. If not, then only to the active part apC
  ! apC     type(tAp) [INOUT]
  !           Current active part
  ! apT     type(tAp) [IN]
  !           Total part of the pencil the method is acting on
  !           (used when compAB is true)
  ! Q       complex double array [INOUT]
  !           Unitary equivalence matrix Q
  ! compQ   boolean [IN]
  !           If true, the left Schur vectors are updated
  ! qcol    integer [IN]
  !           Start index of the columns of Q that are updated:
  !              qcol, qcol+1
  ! Z       double array [INOUT]
  !           Orthonormal equivalence matrix Z
  ! compZ   boolean [IN]
  !           If true, the right Schur vectors are updated
  ! zcol    integer [IN]
  !           Start index of the columns of Z that are updated:
  !              zcol, zcol+1, zcol+2
  !___________________________________________________________________________
  ! last edit: January 22, 2019
  real(kind=dp), intent(inout)      ::  A(:,:), B(:,:), Q(:,:),Z(:,:)
  logical, intent(in)               ::  compAB, compQ, compZ
  integer, intent(in)               ::  qcol, zcol
  logical, intent(out)              ::  defl
  type(tAp), pointer, intent(inout) ::  apC
  type(tAp), pointer, intent(in)    ::  apT

  real(kind=dp)                     :: c, c2, s, s2, r,t
  logical                           :: isreal
  integer,pointer                   :: cstrt, cstp, tstrt, tstp
  integer                           :: LDAB, LDC, LDZ

  external                          :: DGEMM

  call d_getLeadingDim(LDAB,A)
  call d_getLeadingDim(LDC,Cu)
  call d_getLeadingDim(LDZ,Z)

  call tApGet(apC,cstrt,cstp)
  call tApGet(apT,tstrt,tstp)

  defl = .false.

  ! Check first if it is actually a cc pole
  call d_eigenvalues2x2(A,B,cstp-1,cstp-2,c,c2,isreal)
  if (isreal) then
    call d_make_hess(c,cstp-1,cstp-2,A,B,compAB,apC,apT,Q,compQ,qcol,&
                       Z,compZ,zcol)
    call d_stop_deflation_single(defl,A,B,compAB,apC,apT,Z,compZ,zcol)
  else
    ! Compute the RQ factorization
    ! We create a copy of the trailing block for the deflation test
    Q4(1:2,1:3) = A(cstp-1:cstp,cstp-2:cstp)
    Z4(1:2,1:3) = B(cstp-1:cstp,cstp-2:cstp)
    ! Create zero in B(cstp,cstp-1)
    call d_compute_ct(Z4(2,3),Z4(2,2),c,s,r)
    ! Apply transformation
    Z4(2,3) = r
    Z4(2,2) = dzero
    t = s * Z4(1,2) + c * Z4(1,3)
    Z4(1,2) = c * Z4(1,2) - s * Z4(1,3)
    Z4(1,3) = t
    call d_apply_ct_r(Q4(1:2,2),Q4(1:2,3),c,s)
    ! Try to zero column cstp-2 all at once
    if ((abs(Z4(1,1)) .gt. abs(Q4(1,1))) .and. &
       (abs(Z4(1,1)) .gt. abs(Q4(2,1)))) then
      call d_compute_ct(Z4(1,2),Z4(1,1),c2,s2,r)
      Z4(1,1) = dzero
      Z4(1,2) = r
    else
      if (abs(Q4(1,1)) .gt. abs(Q4(2,1))) then
        call d_compute_ct(Q4(1,2),Q4(1,1),c2,s2,r)
      else
        call d_compute_ct(Q4(2,2),Q4(2,1),c2,s2,r)
      end if
      t = s2 * Z4(1,1) + c2 * Z4(1,2)
      Z4(1,1) = c2 * Z4(1,1) - s2 * Z4(1,2)
      Z4(1,2) = t
    end if
    call d_apply_ct_r(Q4(1:2,1),Q4(1:2,2),c2,s2)
    ! Test if deflated
    if (((abs(Q4(1,1)) + abs(Q4(2,1))) .lt. &
          mp * (abs(Q4(1,2)) + abs(Q4(2,3)) + abs(Q4(1,3)) &
          + abs(Q4(2,2)))) .and. (abs(Z4(1,1)) .lt. &
           mp * (abs(Z4(1,2)) + abs(Z4(2,3))+ abs(Z4(1,3))))) then
        defl = .true.
        ! Construct Z
        Z4 = dzero
        Z4(1,1) = done
        Z4(2,2) = done
        Z4(3,3) = done
        call d_apply_ct_r(Z4(1:3,2),Z4(1:3,3),c,s)
        call d_apply_ct_r(Z4(1:3,1),Z4(1:3,2),c2,s2)
    end if

    if (defl) then
      if (compAB) then
        call DGEMM('N','N',cstp-tstrt,3,3,done,A(tstrt+1,cstp-2),LDAB,&
                  Z4(1,1),4,dzero,Cu(1,1),LDC)
        A(tstrt+1:cstp,cstp-2:cstp) = Cu(1:cstp-tstrt,1:3)
        call DGEMM('N','N',cstp-tstrt,3,3,done,B(tstrt+1,cstp-2),LDAB,&
                  Z4(1,1),4,dzero,Cu(1,1),LDC)
        B(tstrt+1:cstp,cstp-2:cstp) = Cu(1:cstp-tstrt,1:3)
      else
        call DGEMM('N','N',cstp-cstrt,3,3,done,A(cstrt+1,cstp-2),LDAB,&
                 Z4(1,1),4,dzero,Cu(1,1),LDC)
        A(cstrt+1:cstp,cstp-2:cstp) = Cu(1:cstp-cstrt,1:3)
        call DGEMM('N','N',cstp-cstrt,3,3,done,B(cstrt+1,cstp-2),LDAB,&
                 Z4(1,1),4,dzero,Cu(1,1),LDC)
        B(cstrt+1:cstp,cstp-2:cstp) = Cu(1:cstp-cstrt,1:3)
      end if
      ! Put elements to zero
      A(cstp-1:cstp,cstp-2) = dzero
      B(cstp-1:cstp,cstp-2) = dzero
      B(cstp,cstp-1) = dzero
      if (compZ) then
        call DGEMM('N','N',LDZ,3,3,done,Z(1,zcol),LDZ,Z4(1,1),4,&
                   dzero,Cu(1,1),LDC)
        Z(1:LDZ,zcol:zcol+2) = Cu(1:LDZ,1:3)
      end if
      ! This also updates it in the caller function
      call tApPut(apC,strt=cstrt,stp=cstp-2)
    end if
  end if
  end subroutine d_stop_deflation_double
end module

max_stars_repo_name: campsd/libRQZ

max_issues_repo_name: campsd/libRQZ

max_forks_repo_name: campsd/libRQZ


text:
C
C $Id: mdfnme.f,v 1.4 2008-07-27 00:17:06 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      CHARACTER*128 FUNCTION MDFNME (IAIN,ILVL)
C
        INTEGER       IAIN,ILVL
C
        PARAMETER (MNAI=8000)
C
C The value of MDFNME(IAIN) is the full name of the area with the area
C identifier IAIN, including the names of all its parents up to the
C the level specified by ILVL.
C
C Assume that IHON is the area identifier for the area named "Honshu",
C which has area type 3, meaning that it exists at the third level in
C the dataset; then it is that case that
C
C   MDFNME (IHON,3) = 'Japan - Honshu'
C   MDFNME (IHON,2) = 'Eurasia - Japan - Honshu'
C   MDFNME (IHON,1) = 'Land - Eurasia - Japan - Honshu'
C
        COMMON /MAPCMX/  IATY(MNAI),ISCI(MNAI),IPAR(MNAI)
        INTEGER          IATY,ISCI,IPAR
        SAVE   /MAPCMX/
C
        COMMON /MAPCMY/  NAME(MNAI),FLNS
        CHARACTER*64     NAME
        CHARACTER*512    FLNS
        SAVE   /MAPCMY/
C
C Declare local variables.
C
        INTEGER          ITMP,NSTP
        CHARACTER*128    CTMP
C
        MDFNME=' '
C
        ITMP=IAIN
        NSTP=0
C
  101   IF (ITMP.GE.1.AND.ITMP.LE.MNAI) THEN
          IF (IATY(ITMP).NE.0) THEN
            IF (NAME(ITMP).NE.' ') THEN
              IF (MDFNME.EQ.' ') THEN
                MDFNME=NAME(ITMP)
              ELSE
                CTMP=MDFNME
                MDFNME=NAME(ITMP)(1:MDILNB(NAME(ITMP)))//' - '//CTMP
              END IF
            END IF
            ITMP=IPAR(ITMP)
            IF (ITMP.GE.1.AND.ITMP.LE.MNAI) THEN
              IF (IATY(ITMP).GE.ILVL) THEN
                NSTP=NSTP+1
                IF (NSTP.LT.10) GO TO 101
              END IF
            END IF
          END IF
        END IF
C
        RETURN
C
      END