        program mud_test_fortran

        include 'mud.f77'

        integer*4 status
        integer*4 i
        character*32 filename
        integer*4 fileHandle
        character*20 title

        integer*4 MUD_hist_hdr

*  Hist header structure, implemented as a common block.
*  All hist header elements are named hh_... The core elements
*  are named hh_c_...

            integer*4   hh_c_pNext              !* pointer to next section *
            integer*4   hh_c_size
            integer*4   hh_c_secID              !* Ident of section type *
            integer*4   hh_c_instanceID
            integer*4   hh_c_sizeof
            integer*4   hh_c_proc
            integer*4   hh_histType
            integer*4   hh_nBytes
            integer*4   hh_nBins
            integer*4   hh_bytesPerBin
            integer*4   hh_fsPerBin
            integer*4   hh_t0_ps
            integer*4   hh_t0_bin
            integer*4   hh_goodBin1
            integer*4   hh_goodBin2
            integer*4   hh_bkgd1
            integer*4   hh_bkgd2
            integer*4   hh_nEvents
            integer*4   hh_pcsTitle

        common /cmn_hdr/
     +   hh_c_pNext, hh_c_size, hh_c_secID, hh_c_instanceID,
     +   hh_c_sizeof, hh_c_proc,
     +   hh_histType, hh_nBytes, hh_nBins, hh_bytesPerBin,
     +   hh_fsPerBin, hh_t0_ps, hh_t0_bin, hh_goodBin1, hh_goodBin2,
     +   hh_bkgd1, hh_bkgd2, hh_nEvents, hh_pcsTitle

        equivalence(hh_c_pNext,MUD_hist_hdr)

        !
        !  Open an MUD format file
        !
        filename = '001234.msr'

        fileHandle = fMUD_openInput( filename )
        if (fileHandle .eq. 0) then
           write (*,*) 'Could not open file ',filename,
     +          '  (',fileHandle,')'
           stop
        endif
        write (*,*) 'Opened file ', filename

        !
        !  Position the file before the first histogram of the
        !  TD histogram group
        !
        status = fMUD_fseek( fileHandle,
     +                       MUD_SEC_GRP_ID, MUD_GRP_TRI_TD_HIST_ID,
     +                       0 )
        if( status .eq. -1 ) then
           write (*,*) 'Failed to find histogram group!  status=',status
           goto 999
        endif
        !
        !  Read the histogram headers
        !
        do i=1,16

            status = fMUD_fseek( fileHandle,
     +                     MUD_SEC_GEN_HIST_HDR_ID, i,
     +                     0 )

            if (status .eq. -1 ) goto 999

            status = fMUD_read( fileHandle, MUD_hist_hdr )

            !
            !  Access the histogram title
            !
            if (status.eq.1) then
               call fMUD_ctofString( title, hh_pcsTitle )
               write (*,*) ' histogram title = <', title, '>'
            else
               write (*,*) ' Failed to read header for histogram',i
            endif

        end do

 999    continue
        call fMUD_close( fileHandle )
        stop
        end


text:
!
!***********************************************************************
!                                                                      *
      PROGRAM rlevelsj
!                                                                      *
!     Purpose:                                                         *
!     Extract energy levels from rscf/rci output files (.m, .cm file)  *
!     and from jj2lsj output file (*.lsj.lbl).                         *
!     And then:                                                        *
!         (1) sort                                                     *
!         (2) difference with the lowest                               *
!         (3) difference with the nearby lower                         *
!     Files names are provided by user.                                *
!     Usage:                                                           *
!           $ grlevels file1 file2 file3 ...                           *
!           or                                                         *
!           $ grlevels                                                 *
!           file1                                                      *
!           file2                                                      *
!           file3                                                      *
!           ... (return to terminate)                                  *
!                                                                      *
!     Calls: LDLBL, indexS.                                            *
!                                                                      *
!     Xinghong He  98-10-16                                            *
!                                                                      *
!     Rewritten by  G. Gaigalas                                        *
!     for LSJ calssification of levels                                 *
!     NIST                                                 May 2011    *
!                                                                      *
!***********************************************************************
      IMPLICIT NONE
      INTEGER, PARAMETER:: JMax = 22      ! max J value, see JFraction !
      INTEGER, PARAMETER:: ndim = 20000   ! max number of states
      INTEGER, PARAMETER:: maxFile = 1000 ! max number of files
      DOUBLE PRECISION, PARAMETER:: Rydberg = 109737.31568508d0
      DOUBLE PRECISION, PARAMETER:: eV = 27.211386018D0
!
      CHARACTER(LEN=80) strInFile(maxFile), strFile
      CHARACTER*1 iaspa(ndim), PlusMinus(-1:1)        ! Parity
      CHARACTER*1 Lev_par(ndim)                       ! Parity
      CHARACTER*4 iatjp(ndim), JFraction(1:2*Jmax+1)  ! J
      CHARACTER*4 Lev_J(ndim)                         ! J
      CHARACTER g92mix*6
      CHARACTER*64 string_CSF(ndim), string_PRN(ndim) ! String in LSJ
      CHARACTER util_lbl_file*256
!
      INTEGER i, j, iargc, ios, ncountState, nFile, mFile
      INTEGER nelec, ncftot, nw, nvectot, nvecsiz, nblock, jblock
      INTEGER nb, ncfblk, nevblk, iiatjp, iiaspa
      INTEGER K, Iprint, IMaxCount
      INTEGER ivec(ndim), indx(ndim), Lev_POS(ndim)
!
      DOUBLE PRECISION eav, eval(ndim), evec, RLev_ENER(ndim),ZERO
!
      COMMON/JJ2LSJ/ Lev_POS,Lev_J,Lev_Par,RLev_ENER,string_CSF,     &
                      IMaxCount
!
      DATA PlusMinus/'-', ' ', '+'/
      DATA JFraction/'  0 ', ' 1/2', '  1 ', ' 3/2', '  2 ', ' 5/2', &
                     '  3 ', ' 7/2', '  4 ', ' 9/2', '  5 ', '11/2', &
                     '  6 ', '13/2', '  7 ', '15/2', '  8 ', '17/2', &
                     '  9 ', '19/2', ' 10 ',                         &
                             '21/2', ' 11 ', '23/2', ' 12 ', '25/2', &
                     ' 13 ', '27/2', ' 14 ', '29/2', ' 15 ', '31/2', &
                     ' 16 ', '33/2', ' 17 ', '35/2', ' 18 ', '37/2', &
                     ' 19 ', '39/2', ' 20 ', '41/2', ' 21 ', '43/2', &
                     ' 22 '/
!
      mFile = iargc()
      string_PRN = ''
      Iprint = 0
      ZERO = 0.0D00
      IF (mFile .EQ. 0) THEN                 ! Get file names interactively
         WRITE (0,*)'  You can also use command-line option:'
         WRITE (0,*)'    %grlevels file1 file2 (wild cards allowed)...'
         WRITE (0,*)'  Now, carry on'
         WRITE (0,*)
         WRITE (0,*)'Type the input file name, one for each line',  &
                  ' (NULL to terminate)'
         WRITE (0,*)
         i = 0
         DO  ! Don't know the exact number of files.
            WRITE (0,'(A12)', ADVANCE='NO') 'File name ? '
            READ (5, '(A80)') strFile
            strFile = ADJUSTL (strFile)
            IF (LEN_TRIM (strFile) .GT. 0) THEN ! a valid input
               i = i + 1
               IF (i .GT. maxFile) THEN      !  impose an upper limit
                  WRITE (0,*)'Too many files opened. processing first '&
                             , i-1, ' files.'
                  EXIT
               ENDIF
               strInFile(i) = strFile
            ELSE
               EXIT
            ENDIF
         ENDDO
         mFile = i
      ELSEIF (mFile .GT. 0 .AND. mFile .LE. maxFile) THEN
         DO i = 1, mFile
            CALL getarg (i, strInFile(i))
         ENDDO
      ELSE
         WRITE (0,*) 'More than ', maxFile, ' files entered,',         &
                  ' modify parameter maxFile'
      ENDIF
!
!     Open files, read energies, concatenate to a single place
!     Open mix file, check header
!
      ncountState = 0
      IMaxCount = 0
      DO nFile = 1, mFile
         strFile = strInFile(nFile)
         OPEN (3, FILE = strFile, FORM = 'UNFORMATTED', STATUS = 'OLD' &
            , IOSTAT = IOS)
         IF (IOS .NE. 0) THEN
           WRITE (0,*) 'Failed to open file "',                   &
                strFile(1:LEN_TRIM (strFile)), '", skipping...'
            CLOSE (3)
            CYCLE
         ENDIF
         READ (3) g92mix
         IF (g92mix .NE. 'G92MIX') THEN
            WRITE (0,*) 'Not a mixing coefficient file, skipping "',&
                      strFile(1:LEN_TRIM (strFile)), '"'
            CLOSE (3)
            CYCLE
         ENDIF
         READ (3) nelec, ncftot, nw, nvectot, nvecsiz, nblock
!
         PRINT *, 'nblock = ', nblock, '  ncftot = ', ncftot,      &
                      '  nw = ', nw, '  nelec = ', nelec

         DO jblock = 1, nblock
            READ (3) nb, ncfblk, nevblk, iiatjp, iiaspa
            IF (jblock .NE. nb) THEN
!
!     This error can occur anywhere and therefore cannot
!     be simply skipped - stop instead.
!
               WRITE (0,*) 'jblock .NE. nb, stopping...'
               CLOSE (3)
               STOP
            ENDIF
            IF (nevblk .LE. 0) CYCLE
            READ (3) (ivec(i+ncountState), i = 1, nevblk)
            READ (3) eav, (eval(i+ncountState), i = 1, nevblk)
            READ (3) (evec, i = 1, ncfblk*nevblk)
!
!     Assign J and parity to every individual state
!     Also add the average energy (back) to energy
!
            DO i = 1, nevblk
               iatjp(i+ncountState) = JFraction(iiatjp)
               iaspa(i+ncountState) = PlusMinus(iiaspa)
               eval(i+ncountState) = eval(i+ncountState) + eav
            END DO
!
!     Update ncountState
!
            ncountState = ncountState + nevblk
         END DO
         CLOSE (3)
!
!     Opens the file  *.lsj.lbl
!
         K = INDEX(strFile,' ')
         if(strFile(K-2:K-1) .EQ. '.m') then
            util_lbl_file = strFile(1:K-3)//'.lsj.lbl'
         else if(strFile(K-3:K-1) .EQ. '.cm') then
            util_lbl_file = strFile(1:K-4)//'.lsj.lbl'
         end if
         OPEN (31, FILE = util_lbl_file, FORM = 'FORMATTED',  &
            STATUS = 'OLD', IOSTAT = IOS)
         IF (IOS .NE. 0) THEN
!GG            WRITE (0,*) 'Failed to open file "',
!GG     &      util_lbl_file(1:LEN_TRIM (util_lbl_file)), '", skipping...'
            CLOSE (31)
            CYCLE
         END IF
         Iprint = 1
!
!     Defines the LSJ string for the levels
!
         CALL LDLBL
         CLOSE (31)
         do i = 1,ncountState
            do j = 1,IMaxCount
               if(DABS(DABS(eval(i))-DABS(RLev_ENER(j))) &
                                 .LT. 0.00000001) then
                   string_PRN(i) =  string_CSF(j)
                   exit
               end if
            end do
         end do
      END DO
      CALL indexS (ncountState, eval, .FALSE., indx)
!
!     The output of the levels
!
      PRINT *
      WRITE (6,1)
      WRITE (6,2) Rydberg
      if(Iprint .eq. 1) then
         WRITE (6,*) 'Splitting is the energy difference ',  &
                     'with the lower neighbor'
         WRITE (6,5)
         WRITE (6,*) 'No Pos  J Parity Energy Total    Levels', &
                     '     Splitting     Configuration'
         WRITE (6,*) '                     (a.u.)         (eV)', &
                     '         (eV)'
         WRITE (6,5)
         j = 1
         i = indx(j)
         WRITE (6,3) j,ivec(i),iatjp(i),iaspa(i),eval(i), &
            ZERO,ZERO,                                    &
            Trim(string_PRN(i))
         DO j = 2, ncountState
            i = indx(j)
            WRITE (6,3) j,ivec(i),iatjp(i),iaspa(i),eval(i),  &
            (eval(i)-eval(indx(1)))*eV,                       &
            (eval(i)-eval(indx(j-1)))*eV,                     &
            Trim(string_PRN(i))
         END DO
         WRITE (6,5)
      ELSE
         WRITE (6,*) 'No - Serial number of the state; ',    &
                     'Pos - Position of the state within the '
         WRITE (6,*) 'J/P block; Splitting is the energy difference ',&
                     'with the lower neighbor'
         WRITE (6,6)
         WRITE (6,*) 'No Pos  J Parity Energy Total    Levels',  &
                     '     Splitting '
         WRITE (6,*) '                     (a.u.)         (eV)', &
                     '         (eV)'
         WRITE (6,6)
         j = 1
         i = indx(j)
         WRITE (6,3) j,ivec(i),iatjp(i),iaspa(i),eval(i)
         DO j = 2, ncountState
            i = indx(j)
            WRITE (6,3) j,ivec(i),iatjp(i),iaspa(i),eval(i),  &
            (eval(i)-eval(indx(1)))*eV,                       &
            (eval(i)-eval(indx(j-1)))*eV
         END DO
         WRITE (6,6)
      END IF
   1  FORMAT (' Energy levels for ...')
   2  FORMAT (' Rydberg constant is ',F14.5)
   3  FORMAT (2I3,1X,A4,1x,A1,2X,F14.7,F15.5,F15.5,2X,A)
   5  FORMAT ('---------------------------------------------',&
              '---------------------------------------------')
   6  FORMAT ('------------------------------------------',   &
              '-------------------------------')
      STOP

      CONTAINS
!
!***********************************************************************
!                                                                      *
      SUBROUTINE indexS(n,a,ldown,indx)
!                                                                      *
!     Sort out the order of array a and store the index in indx        *
!                                                        (a pointer)   *
!     The input array a is unchanged  written in the bases of UpDown   *
!                                                                      *
!     !$Id: rlevels.f,v 1.2 2003/10/02 07:56:22 per Exp $              *
!     $Log: rlevels.f,v $                                              *
!     Revision 1.2  2003/10/02 07:56:22  per                           *
!     *** empty log message ***                                        *
!                                                                      *
!     Revision 1.1.1.1  2003/01/04 21:45:39  georgio                   *
!                                                                      *
!***********************************************************************
      IMPLICIT NONE
      LOGICAL          ldown              ! .TRUE. then Big ---> Small
      INTEGER          n, indx(n)
      DOUBLE PRECISION a(n)
      INTEGER          i, j, ipos, jpos, jhere
      DOUBLE PRECISION aimx
!
!     Initialize the index array
      DO i = 1, n
         indx(i) = i
      ENDDO
      IF (ldown) THEN
         DO i = 1, n
            ipos = indx(i)
            aimx = a(ipos)
            jhere = i
            DO j = i+1, n
               jpos = indx(j)
               IF(a(jpos) .GT. aimx) THEN
                  aimx = a(jpos)
                  jhere = j
               ENDIF
            ENDDO
            indx(i) = indx(jhere)
            indx(jhere) = ipos
         ENDDO
      ELSE
         DO i = 1, n
            ipos = indx(i)
            aimx = a(ipos)
            jhere = i
            DO j = i+1, n
               jpos = indx(j)
               IF(a(jpos) .LT. aimx) THEN
                  aimx = a(jpos)
                  jhere = j
               ENDIF
            ENDDO
            indx(i) = indx(jhere)
            indx(jhere) = ipos
         ENDDO
      ENDIF
      RETURN
      END SUBROUTINE  IndexS
!
!***********************************************************************
!                                                                      *
      SUBROUTINE LDLBL
!                                                                      *
!     Open, check and load data from the  .lsj.lbl   file of the       *
!     inital state.                                                    *
!                                                                      *
!     Calls:                                                           *
!                                                                      *
!     Written by G. Gaigalas,                                          *
!     NIST                                                  May 2011   *
!                                                                      *
!***********************************************************************
!
      IMPLICIT NONE
      INTEGER, PARAMETER :: ndim = 20000          ! max number of states
!
      CHARACTER*1 Lev_par(ndim)
      CHARACTER*4 Lev_J(ndim)
      CHARACTER*15 RECORD
      CHARACTER*64 string_CSF(ndim)
!
      INTEGER IOS, ITEST, Lev_POS(ndim), Icount, ImaxCount
      REAL WEIGHTS
      DOUBLE PRECISION RLev_ENER(ndim)
!
      COMMON/JJ2LSJ/ Lev_POS,Lev_J,Lev_Par,RLev_ENER,string_CSF, &
                      IMaxCount
!
      READ (31,'(1A15)',IOSTAT = IOS) RECORD
      ICount = 0
      IF (IOS .NE. 0) GO TO 1
      ICount = 1
      READ (31,'(1X,I2,1X,A4,5X,A1,8X,F16.9)',IOSTAT = IOS)  &
        Lev_Pos(ICount),Lev_J(ICount),Lev_Par(ICount),       &
        RLev_ENER(ICount)
      IF (IOS .NE. 0) GO TO 1
!
      READ (31,'(7X,F12.8,17X,A)') WEIGHTS,string_CSF(ICount)
!
    2 READ (31,'(1X,I2)',IOSTAT = IOS) ITEST
      IF (IOS .NE. 0) GO TO 1
      IF (ITEST .EQ. 0) GO TO 2
      BACKSPACE 31
      ICount = ICount + 1
      READ (31,'(1X,I2,1X,A4,5X,A1,8X,F16.9)',IOSTAT = IOS) &
        Lev_Pos(ICount),Lev_J(ICount),Lev_Par(ICount),      &
        RLev_ENER(ICount)
      IF (IOS .NE. 0) GO TO 1
      READ (31,'(7X,F12.8,17X,A)') WEIGHTS,string_CSF(ICount)
      GO TO 2
    1 CONTINUE
      IMaxCount = ICount + IMaxCount
      RETURN
      END  SUBROUTINE LDLBL
      END PROGRAM

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: sylas/grasp-continuum

max_forks_repo_name: sylas/grasp-continuum


!-------------------------------------------------------------------------------
!> module User
!!
!! @par Description
!!          calc perturbation
!!
!! @author Team SCALE
!!
!! @par History
!! @li      2013-09-05 (H.Yashiro)   [new]
!!
!<
!-------------------------------------------------------------------------------
module mod_user
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_stdio
  use scale_prof
  use scale_grid_index
  use scale_tracer
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: USER_setup0
  public :: USER_setup
  public :: USER_step

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  logical,  private, save :: USER_do = .false. !< do user step?

  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  !> Setup0
  subroutine USER_setup0
  end subroutine USER_setup0

  !-----------------------------------------------------------------------------
  !> Setup
  subroutine USER_setup
    use scale_process, only: &
       PRC_MPIstop
    implicit none

    namelist / PARAM_USER / &
       USER_do

    integer :: ierr
    !---------------------------------------------------------------------------

    if( IO_L ) write(IO_FID_LOG,*)
    if( IO_L ) write(IO_FID_LOG,*) '+++ Module[USER]/Categ[MAIN]'

    !--- read namelist
    rewind(IO_FID_CONF)
    read(IO_FID_CONF,nml=PARAM_USER,iostat=ierr)

    if( ierr < 0 ) then !--- missing
       if( IO_L ) write(IO_FID_LOG,*) '*** Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
       write(*,*) 'xxx Not appropriate names in namelist PARAM_USER. Check!'
       call PRC_MPIstop
    endif
    if( IO_LNML ) write(IO_FID_LOG,nml=PARAM_USER)

    call USER_step

    return
  end subroutine USER_setup

  !-----------------------------------------------------------------------------
  !> Step
  subroutine USER_step
    use scale_process, only: &
       PRC_MPIstop
    use scale_grid_real, only : &
       CZ => REAL_CZ, &
       FZ => REAL_FZ
    use scale_time, only : &
       NOWSEC => TIME_NOWSEC
    use mod_atmos_vars, only: &
       QTRC
    implicit none

    integer  :: modsec
    real(RP) :: dist

    integer  :: k, i, j
    !---------------------------------------------------------------------------

    if ( USER_do ) then
       if( IO_L ) write(IO_FID_LOG,*) '*** Add rain.'
       USER_do = .false.

       do j = JS, JE
       do i = IS, IE
       do k = KS, KE
          dist = ( ( CZ(k,i,j) - 3000.0_RP ) / 50.0_RP )**2

          QTRC(k,i,j,I_QR) = 1.E-3 / ( 1.0_RP + dist )
!          if (       FZ(k,  i,j) >= 3000.0_RP &
!               .AND. FZ(k-1,i,j) <  3000.0_RP ) then
!
!             if( IO_L ) write(IO_FID_LOG,*) k,i,j,CZ(k,i,j),dist,1.E-3/( 1.0_RP + dist )
!             QTRC(k,i,j,I_QR) = 1.E-3
!          endif
       enddo
       enddo
       enddo
    endif

    return
  end subroutine USER_step

end module mod_user


text:
      SUBROUTINE MLOSS
*
*
*       Mass loss from evolving stars.
*       ------------------------------
*
*       Original scheme of Elena Terlevich 1983.
*       ----------------------------------------
*
      INCLUDE 'common6.h'
      REAL*8  A(9)
*
*
*       Find the heaviest body (exclude any merged binary components).
   99 IMAX = 0
    1 BMAX = 0.0D0
      DO 2 I = 1,N
          IF (BODY(I).LT.BMAX.OR.I.LE.IMAX) GO TO 2
          BMAX = BODY(I)
          IMAX = I
    2 CONTINUE
*
      ZMSTAR = ZMBAR*BMAX
*       Obtain total evolution time for ZMSTAR in solar masses.
      TMS = (2.55D+03 + 6.69D2*ZMSTAR**2.5D0 + ZMSTAR**4.5D0)/
     &      (3.27D-02*ZMSTAR**1.5D0 + 3.46D-01*ZMSTAR**4.5D0)
      TG = 0.15*TMS
      THE = TMS*1.37*ZMSTAR**(-0.881D0)
      TEV1 = TMS + TG + THE
*       Scale the evolution time to model units.
      TMDOT = TEV1/TSCALE
*
      IF (TMDOT.GT.TIME) THEN
*       Update the maximum mass.
          BODY1 = BMAX
*       Set phase indicator = -1 for new time-step list in INTGRT.
          IPHASE = -1
          RETURN
      END IF
*
      IBODY = IMAX
*       Include special treatment for regularized components.
      IF (IMAX.GT.2*NPAIRS) GO TO 4
      IPAIR = KVEC(IMAX)
      IBODY = N + IPAIR
*       Continue search if merged binary component is identified.
      IF (NAME(IBODY).LT.0) GO TO 1
*
*       Obtain current coordinates & velocities (unperturbed KS pair OK).
      CALL RESOLV(IPAIR,2)
*
*       Copy X0DOT since routine RESOLV skipped in KSTERM (IPHASE = -1).
      DO 3 K = 1,3
          X0DOT(K,2*IPAIR-1) = XDOT(K,2*IPAIR-1)
          X0DOT(K,2*IPAIR) = XDOT(K,2*IPAIR)
    3 CONTINUE
*
*       Set variable white dwarf mass (Iben & Renzini, Ann. Rev. 21, 298).
    4 ZMSTAR = 0.38 + 0.15*ZMSTAR
*       Assume neutron star instead if mass > 6 solar masses.
      IF (ZMBAR*BMAX.GT.6.0) ZMSTAR = 1.5
      DM = BODY(IMAX) - ZMSTAR/ZMBAR
      ZMASS = ZMASS - DM
      ZMDOT = ZMDOT + DM*ZMBAR
      NMDOT = NMDOT + 1
      BODY(IMAX) = ZMSTAR/ZMBAR
      IF (TIDAL(1).NE.0.0D0) RTIDE = (ZMASS/TIDAL(1))**0.3333
*
*       Save old velocity for FDOT correction.
      VI2 = 0.0
      DO 5 K = 1,3
          A(K+6) = XDOT(K,IMAX)
          VI2 = VI2 + XDOT(K,IMAX)**2
    5 CONTINUE
*
      DE = 0.5*DM*(VI2 - TIDAL(1)*X(1,IMAX)**2 - TIDAL(3)*X(3,IMAX)**2)
      A0 = 1.0
      IF (ZMBAR*BMAX.LE.6.0) GO TO 8
*
*       Assign a high velocity to neutron star (at least 4*rms velocity).
      A0 = 4.0*SQRT(0.5D0*ZMASS/(RSCALE*VI2))
      IF (IMAX.GT.2*NPAIRS) GO TO 6
      A2 = 2.0*(BODY(2*IPAIR-1) + BODY(2*IPAIR))/R(IPAIR)
*       Add escape velocity from the regularized pair.
      A0 = SQRT(A0**2 + A2/VI2)
    6 DO 7 K = 1,3
          XDOT(K,IMAX) = A0*XDOT(K,IMAX)
          X0DOT(K,IMAX) = XDOT(K,IMAX)
    7 CONTINUE
*
*       Correct total energy, forces & first derivatives.
    8 POTI = 0.0D0
      DO 20 J = 1,NTOT
          IF (J.EQ.IMAX) GO TO 20
          RIJDOT = 0.0D0
          RDVDOT = 0.0D0
*
          DO 10 K = 1,3
              A(K) = X(K,IMAX) - X(K,J)
              A(K+3) = A(K+6) - XDOT(K,J)
              RIJDOT = RIJDOT + A(K)*A(K+3)
              RDVDOT = RDVDOT + A(K)*(XDOT(K,IMAX) - A(K+6))
   10     CONTINUE
*
          RIJ2 = A(1)**2 + A(2)**2 + A(3)**2
          IF (J.LE.N) POTI = POTI + BODY(J)/SQRT(RIJ2)
          IF (J.LE.2*NPAIRS.OR.J.EQ.IBODY) GO TO 20
          A3 = 1.0/(RIJ2*SQRT(RIJ2))
          A4 = BODY(IMAX)*A3
          A5 = DM*A3
          A6 = 3.0*RIJDOT/RIJ2
          A7 = 3.0*RDVDOT/RIJ2
*
          DO 11 K = 1,3
              A(K+3) = (A(K+3) - A(K)*A6)*A5
              IF (A0.GT.1.0) THEN
*       Include FDOT corrections due to increased velocity.
                  A(K+3) = A(K+3) + (XDOT(K,IMAX) - A(K+6))*A4
                  A(K+3) = A(K+3) - A7*A(K)*A4
              END IF
   11     CONTINUE
*
*       Use neighbour list to distinguish irregular & regular terms.
          NNB = LIST(1,J) + 1
          DO 12 L = 2,NNB
              IF (LIST(L,J).EQ.IBODY)  GO TO 16
              IF (LIST(L,J).GT.IBODY)  GO TO 13
   12     CONTINUE
*
   13     DO 14 K = 1,3
              F(K,J) = F(K,J) - 0.5*A(K)*A5
              FR(K,J) = FR(K,J) - A(K)*A5
              FDOT(K,J) = FDOT(K,J) - ONE6*A(K+3)
              D1R(K,J) = D1R(K,J) - A(K+3)
   14     CONTINUE
          GO TO 20
*
   16     DO 18 K = 1,3
              F(K,J) = F(K,J) - 0.5*A(K)*A5
              FI(K,J) = FI(K,J) - A(K)*A5
              FDOT(K,J) = FDOT(K,J) - ONE6*A(K+3)
              D1(K,J) = D1(K,J) - A(K+3)
   18     CONTINUE
   20 CONTINUE
*
*       Correct total energy for mass loss effect.
      DE = DE - DM*POTI
      BE(3) = BE(3) - DE
      RI = SQRT((X(1,IMAX) - RDENS(1))**2 + (X(2,IMAX) - RDENS(2))**2 +
     &                                      (X(3,IMAX) - RDENS(3))**2)
      if(rank.eq.0)
     &WRITE (6,30)  NAME(IMAX), BMAX, BMAX*ZMBAR, ZMDOT, DE, BE(3),
     &              TMDOT, RI, LIST(1,IMAX)
   30 FORMAT (/,'   MASS LOSS   NAME =',I5,'  MI =',F8.4,'  M* =',F5.1,
     &                       '  ZM* =',F6.1,'  DE =',F9.5,'  E =',F10.6,
     &                      '  TMDOT =',F6.1,'  RI =',F5.2,'  NNB =',I3)
*
*       Set option = 2 to skip next energy check (reduced by 1 in ADJUST).
      KZ(19) = 2
*
      IF (ZMBAR*BMAX.LE.6.0) GO TO 50
*       Include correction to total energy from the increased velocity.
      DE = 0.5*BODY(IMAX)*VI2*(A0**2 - 1.0)
      BE(3) = BE(3) + DE
      if(rank.eq.0)
     &WRITE (6,45)  IMAX, DE, TEV1, A0, SQRT(VI2), BE(3)
   45 FORMAT ('   RECOIL   I =',I5,'  DE =',F9.5,'  T* =',F5.1,
     &                      '  V/V0 =',F6.2,'  V0 =',F5.2,'  E =',F10.6)
*
   50 IF (IMAX.LE.2*NPAIRS) THEN
*       Save pair index in KSPAIR and terminate regularization.
          KSPAIR = IPAIR
          IPHASE = -1
*       Indicator for skipping routine RESOLV to preserve velocity X0DOT.
          CALL KSTERM
      ELSE
*       Reduce the steps by velocity factor.
          STEP(IMAX) = MAX(STEP(IMAX)/A0,TIME - T0(IMAX))
          STEPR(IMAX) = MAX(STEPR(IMAX)/A0,STEP(IMAX))
      END IF
*
*       Set next mass loss time & current maximum mass before returning.
      GO TO 99
*
      END

ext: f

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
!-----------------------------------------------------------------------
! Interface file for push1.f
      module push1_h
      implicit none
!
      interface
         subroutine DISTR1(part,vtx,vdx,npx,idimp,nop,nx,ipbc)
         implicit none
         integer, intent(in) :: npx, idimp, nop, nx, ipbc
         real, intent(in) :: vtx, vdx
         real, dimension(idimp,nop), intent(inout) :: part
         end subroutine
      end interface
!
      interface
         subroutine GPUSH1L(part,fx,qbm,dt,ek,idimp,nop,nx,nxv,ipbc)
         implicit none
         integer, intent(in) :: idimp, nop, nx, nxv, ipbc
         real, intent(in) :: qbm, dt
         real, intent(inout) :: ek
         real, dimension(idimp,nop), intent(inout) :: part
         real, dimension(nxv), intent(in) :: fx
         end subroutine
      end interface
!
      interface
         subroutine GPOST1L(part,q,qm,nop,idimp,nxv)
         implicit none
         integer, intent(in) :: nop, idimp, nxv
         real, intent(in) :: qm
         real, dimension(idimp,nop), intent(in) :: part
         real, dimension(nxv), intent(inout) :: q
         end subroutine
      end interface
!
      interface
         subroutine DSORTP1XL(parta,partb,npic,idimp,nop,nx1)
         implicit none
         integer, intent(in) :: idimp, nop, nx1
         real, dimension(idimp,nop), intent(in) :: parta
         real, dimension(idimp,nop), intent(inout) :: partb
         integer, dimension(nx1), intent(inout) :: npic
         end subroutine
      end interface
!
      interface
         subroutine CGUARD1L(fx,nx,nxe)
         implicit none
         integer, intent(in) :: nx, nxe
         real, dimension(nxe), intent(inout) :: fx
         end subroutine
      end interface
!
      interface
         subroutine AGUARD1L(q,nx,nxe)
         implicit none
         integer, intent(in) :: nx, nxe
         real, dimension(nxe), intent(inout) :: q
         end subroutine
      end interface
!
      interface
         subroutine POIS1(q,fx,isign,ffc,ax,affp,we,nx)
         implicit none
         integer, intent(in) :: isign, nx
         real, intent(in) :: ax, affp
         real, intent(inout) :: we
         real, dimension(nx), intent(in) :: q
         real, dimension(nx), intent(inout) :: fx
         complex, dimension(nx/2), intent(inout) :: ffc
         end subroutine
      end interface
!
      interface
         subroutine WFFT1RINIT(mixup,sct,indx,nxhd)
         implicit none
         integer, intent(in) :: indx, nxhd
         integer, dimension(nxhd), intent(inout) :: mixup
         complex, dimension(nxhd), intent(inout) :: sct
         end subroutine
      end interface
!
      interface
         subroutine FFT1RXX(f,t,isign,mixup,sct,indx,nxd,nxhd)
         implicit none
         integer, intent(in) :: isign, indx, nxd, nxhd
         real, dimension(nxd), intent(inout) :: f, t
         integer, dimension(nxhd), intent(in) :: mixup
         complex, dimension(nxhd), intent(in) :: sct
         end subroutine
      end interface
!
      interface
         function ranorm()
         implicit none
         double precision :: ranorm
         end function
      end interface
!
      interface
         function randum()
         implicit none
         double precision :: randum
         end function
      end interface
!
      end module

max_stars_repo_name: gcasabona/cuda

max_issues_repo_name: gcasabona/cuda

max_forks_repo_name: gcasabona/cuda


text:
!> The FlowModule contains types which define material flows within the surface water network.
!! Separate types are provided for water, SPM, NM and dissolved species flows.
module FlowModule
    use Globals, only: dp, C
    implicit none

    !> The WaterFlows object stores information on water flows in and out of a reach.
    type, public :: WaterFlows
        real(dp) :: inflow
        real(dp) :: runoff
        real(dp) :: transfers
        real(dp) :: demands
        real(dp) :: outflow
      contains
        procedure :: init => initWaterFlows
        procedure :: empty => emptyWaterFlows
        procedure :: asArray => asArrayWaterFlows
        procedure :: assignWaterFlows
        generic :: assignment(=) => assignWaterFlows
    end type

    !> The SPMFlows object stores information on SPM flows in and out of a reach.
    type, public :: SPMFlows
        real(dp), allocatable :: inflow(:)
        real(dp), allocatable :: soilErosion(:)
        real(dp), allocatable :: bankErosion(:)
        real(dp), allocatable :: transfers(:)
        real(dp), allocatable :: demands(:)
        real(dp), allocatable :: deposition(:)
        real(dp), allocatable :: resuspension(:)
        real(dp), allocatable :: outflow(:)
      contains
        procedure :: init => initSPMFlows
        procedure :: empty => emptySPMFlows
        procedure :: asArray => asArraySPMFlows
        procedure :: assignSPMFlows
        generic :: assignment(=) => assignSPMFlows
    end type

    !> The NMFlows object stores information on NM flows in and out of a reach.
    type, public :: NMFlows
        real(dp), allocatable :: inflow(:,:,:)
        real(dp), allocatable :: soilErosion(:,:,:)
        real(dp), allocatable :: bankErosion(:,:,:)
        real(dp), allocatable :: transfers(:,:,:)
        real(dp), allocatable :: demands(:,:,:)
        real(dp), allocatable :: deposition(:,:,:)
        real(dp), allocatable :: resuspension(:,:,:)
        real(dp), allocatable :: outflow(:,:,:)
        real(dp), allocatable :: pointSources(:,:,:)
        real(dp), allocatable :: diffuseSources(:,:,:)
      contains
        procedure :: init => initNMFlows
        procedure :: empty => emptyNMFlows
        procedure :: asArray => asArrayNMFlows
        procedure :: assignNMFlows
        generic :: assignment(=) => assignNMFlows
    end type

    !> The DissolvedFlows object stores information on dissolved species flows in and out of a reach.
    type, public :: DissolvedFlows
        real(dp), allocatable :: inflow
        real(dp), allocatable :: transfers
        real(dp), allocatable :: demands
        real(dp), allocatable :: outflow
        real(dp), allocatable :: pointSources
        real(dp), allocatable :: diffuseSources
      contains
        procedure :: init => initDissolvedFlows
        procedure :: empty => emptyDissolvedFlows
        procedure :: asArray => asArrayDissolvedFlows
        procedure :: assignDissolvedFlows
        generic :: assignment(=) => assignDissolvedFlows
    end type

  contains

    subroutine initWaterFlows(me)
        class(WaterFlows) :: me
        call me%empty()
    end subroutine

    subroutine initSPMFlows(me)
        class(SPMFlows) :: me
        allocate(me%inflow(C%nSizeClassesSPM))
        allocate(me%soilErosion(C%nSizeClassesSPM))
        allocate(me%bankErosion(C%nSizeClassesSPM))
        allocate(me%transfers(C%nSizeClassesSPM))
        allocate(me%demands(C%nSizeClassesSPM))
        allocate(me%deposition(C%nSizeClassesSPM))
        allocate(me%resuspension(C%nSizeClassesSPM))
        allocate(me%outflow(C%nSizeClassesSPM))
        call me%empty()
    end subroutine

    subroutine initNMFlows(me)
        class(NMFlows) :: me
        allocate(me%inflow(C%npDim(1),C%npDim(2),C%npDim(3)))
        allocate(me%soilErosion(C%npDim(1),C%npDim(2),C%npDim(3)))
        allocate(me%bankErosion(C%npDim(1),C%npDim(2),C%npDim(3)))
        allocate(me%transfers(C%npDim(1),C%npDim(2),C%npDim(3)))
        allocate(me%demands(C%npDim(1),C%npDim(2),C%npDim(3)))
        allocate(me%deposition(C%npDim(1),C%npDim(2),C%npDim(3)))
        allocate(me%resuspension(C%npDim(1),C%npDim(2),C%npDim(3)))
        allocate(me%outflow(C%npDim(1),C%npDim(2),C%npDim(3)))
        allocate(me%pointSources(C%npDim(1),C%npDim(2),C%npDim(3)))
        allocate(me%diffuseSources(C%npDim(1),C%npDim(2),C%npDim(3)))
        call me%empty()
    end subroutine

    subroutine initDissolvedFlows(me)
        class(DissolvedFlows) :: me
        call me%empty()
    end subroutine

    subroutine emptyWaterFlows(me)
        class(WaterFlows) :: me
        me%inflow = 0.0_dp
        me%runoff = 0.0_dp
        me%transfers = 0.0_dp
        me%demands = 0.0_dp
        me%outflow = 0.0_dp
    end subroutine

    subroutine emptySPMFlows(me)
        class(SPMFlows) :: me
        me%inflow = 0.0_dp
        me%soilErosion = 0.0_dp
        me%bankErosion = 0.0_dp
        me%transfers = 0.0_dp
        me%demands = 0.0_dp
        me%deposition = 0.0_dp
        me%resuspension = 0.0_dp
        me%outflow = 0.0_dp
    end subroutine

    subroutine emptyNMFlows(me)
        class(NMFlows) :: me
        me%inflow = 0.0_dp
        me%soilErosion = 0.0_dp
        me%bankErosion = 0.0_dp
        me%transfers = 0.0_dp
        me%demands = 0.0_dp
        me%deposition = 0.0_dp
        me%resuspension = 0.0_dp
        me%outflow = 0.0_dp
        me%pointSources = 0.0_dp
        me%diffuseSources = 0.0_dp
    end subroutine

    subroutine emptyDissolvedFlows(me)
        class(DissolvedFlows) :: me
        me%inflow = 0.0_dp
        me%transfers = 0.0_dp
        me%demands = 0.0_dp
        me%outflow = 0.0_dp
        me%diffuseSources = 0.0_dp
        me%pointSources = 0.0_dp
    end subroutine

    function asArrayWaterFlows(me) result(arr)
        class(WaterFlows)   :: me
        real(dp)            :: arr(5)
        arr = [me%inflow, me%runoff, me%transfers, me%demands, me%outflow]
    end function

    function asArraySPMFlows(me) result(arr)
        class(SPMFlows) :: me
        real(dp)        :: arr(8,C%nSizeClassesSpm)
        arr(1,:) = me%inflow
        arr(2,:) = me%soilErosion
        arr(3,:) = me%bankErosion
        arr(4,:) = me%transfers
        arr(5,:) = me%demands
        arr(6,:) = me%deposition
        arr(7,:) = me%resuspension
        arr(8,:) = me%outflow
    end function

    function asArrayNMFlows(me) result(arr)
        class(NMFlows)  :: me
        real(dp)        :: arr(10,C%npDim(1),C%npDim(2),C%npDim(3))
        arr(1,:,:,:) = me%inflow
        arr(2,:,:,:) = me%soilErosion
        arr(3,:,:,:) = me%bankErosion
        arr(4,:,:,:) = me%transfers
        arr(5,:,:,:) = me%demands
        arr(6,:,:,:) = me%deposition
        arr(7,:,:,:) = me%resuspension
        arr(8,:,:,:) = me%outflow
        arr(9,:,:,:) = me%pointSources
        arr(10,:,:,:) = me%diffuseSources
    end function

    function asArrayDissolvedFlows(me) result(arr)
        class(DissolvedFlows)   :: me
        real(dp)                :: arr(6)
        arr = [me%inflow, me%transfers, me%demands, me%outflow, me%pointSources, me%diffuseSources]
    end function

    subroutine assignWaterFlows(obj, arr)
        class(WaterFlows), intent(out)  :: obj
        real(dp), intent(in)            :: arr(5)
        obj%inflow = arr(1)
        obj%runoff = arr(2)
        obj%transfers = arr(3)
        obj%demands = arr(4)
        obj%outflow = arr(5)
    end subroutine

    subroutine assignSPMFlows(obj, arr)
        class(SPMFlows), intent(out)    :: obj
        real(dp), intent(in)            :: arr(8,C%nSizeClassesSpm)
        obj%inflow = arr(1,:)
        obj%soilErosion = arr(2,:)
        obj%bankErosion = arr(3,:)
        obj%transfers = arr(4,:)
        obj%demands = arr(5,:)
        obj%deposition = arr(6,:)
        obj%resuspension = arr(7,:)
        obj%outflow = arr(8,:)
    end subroutine

    subroutine assignNMFlows(obj, arr)
        class(NMFlows), intent(out) :: obj
        real(dp), intent(in)        :: arr(10,C%npDim(1),C%npDim(2),C%npDim(3))
        obj%inflow = arr(1,:,:,:)
        obj%soilErosion = arr(2,:,:,:)
        obj%bankErosion = arr(3,:,:,:)
        obj%transfers = arr(4,:,:,:)
        obj%demands = arr(5,:,:,:)
        obj%deposition = arr(6,:,:,:)
        obj%resuspension = arr(7,:,:,:)
        obj%outflow = arr(8,:,:,:)
        obj%pointSources = arr(9,:,:,:)
        obj%diffuseSources = arr(10,:,:,:)
    end subroutine

    subroutine assignDissolvedFlows(obj, arr)
        class(DissolvedFlows), intent(out)  :: obj
        real(dp), intent(in)                :: arr(6)
        obj%inflow = arr(1)
        obj%transfers = arr(2)
        obj%demands = arr(3)
        obj%outflow = arr(4)
        obj%pointSources = arr(5)
        obj%diffuseSources = arr(6)
    end subroutine

end module

max_stars_repo_name: NERC-CEH/nanofase

max_issues_repo_name: NERC-CEH/nanofase

max_forks_repo_name: NERC-CEH/nanofase


text:
program phaml_master
use phaml
implicit none
type(phaml_solution_type) :: soln
call phaml_create(soln,nproc=4)
call phaml_solve_pde(soln,print_grid_when=FINAL,print_grid_who=MASTER, &
                     print_error_when=FINAL,print_error_who=MASTER, &
                     print_errest_what=ENERGY_ERREST, &
                     term_energy_err=4.0e-2_my_real)
call phaml_destroy(soln)
end program phaml_master

max_stars_repo_name: qsnake/phaml

max_issues_repo_name: qsnake/phaml

max_forks_repo_name: qsnake/phaml


submodule (grid:grid_read) readgrid_hdf5

use phys_consts, only: debug
use h5fortran, only: hdf5_file,h5exist

implicit none (type, external)

contains

module procedure get_grid3_coords_hdf5

  type(hdf5_file) :: hf
  character(:), allocatable :: fn

  if (index(path, 'simgrid.h5') /= 0) then
    fn = path
  else
    fn = path // '/simgrid.h5'
  endif
  if (debug) print '(A,/,A)', 'READ 3D (B-parallel, B-perp, B-perp) grid:', fn

  call hf%open(fn, action='r')

  !> reads common to 2D and 3D

  !> 1-D variables
  call hf%read('/x1', x1)
  call hf%read('/x2', x2all)
  call hf%read('/x3', x3all)

  if (.not. (hf%exist('/glonctr') .and. &
             hf%exist('/glatctr'))) error stop "please specify /glonctr and /glatctr in " // fn

  call hf%read('/glonctr',glonctr)
  call hf%read('/glatctr',glatctr)

  call hf%close()

end procedure get_grid3_coords_hdf5

!module procedure get_grid3_hdf5
!
!  type(hdf5_file) :: hout
!  character(:), allocatable :: fn
!
!  if (index(path, 'simgrid.h5') /= 0) then
!    fn = path
!  else
!    fn = path // '/simgrid.h5'
!  endif
!  if (debug) print '(A,/,A)', 'READ 3D (B-parallel, B-perp, B-perp) grid:', fn
!
!  call hout%open(fn, action='r')
!
!  !> reads common to 2D and 3D
!
!  !> 1-D variables
!  call hout%read('/x1', x%x1)
!  call hout%read('/x1i', x%x1i)
!  call hout%read('/dx1b', x%dx1)
!  call hout%read('/dx1h', x%dx1i)
!
!  if (flagswap/=1) then
!    !! normal (i.e. full 3D) grid ordering, or a 2D grid with 1 element naturally in the second dimension
!
!    call hout%read('/x2', x%x2all)
!    call hout%read('/x2i', x%x2iall)
!    call hout%read('/dx2b', x%dx2all)
!    call hout%read('/dx2h', x%dx2iall)
!
!    call hout%read('/x3', x%x3all)
!    call hout%read('/x3i', x%x3iall)
!    call hout%read('/dx3b', x%dx3all)
!    call hout%read('/dx3h', x%dx3iall)
!
!    !> 3D variables
!    call hout%read('/h1', x%h1all)
!    call hout%read('/h2', x%h2all)
!    call hout%read('/h3', x%h3all)
!
!    call hout%read('/h1x1i', x%h1x1iall)
!    call hout%read('/h2x1i', x%h2x1iall)
!    call hout%read('/h3x1i', x%h3x1iall)
!
!    call hout%read('/h1x2i', x%h1x2iall)
!    call hout%read('/h2x2i', x%h2x2iall)
!    call hout%read('/h3x2i', x%h3x2iall)
!
!    call hout%read('/h1x3i', x%h1x3iall)
!    call hout%read('/h2x3i', x%h2x3iall)
!    call hout%read('/h3x3i', x%h3x3iall)
!
!    call hout%read('/gx1', g1all)
!    call hout%read('/gx2', g2all)
!    call hout%read('/gx3', g3all)
!
!    call hout%read('/alt', x%altall)
!    call hout%read('/glat', glatall)
!    call hout%read('/glon', x%glonall)
!
!    call hout%read('/Bmag', x%Bmagall)
!    call hout%read('/I', Incall)
!    call hout%read('/nullpts', nullptsall)
!
!    call hout%read('/e1', e1all)
!    call hout%read('/e2', e2all)
!    call hout%read('/e3', e3all)
!
!    call hout%read('/er', erall)
!    call hout%read('/etheta', ethetaall)
!    call hout%read('/ephi', ephiall)
!
!    call hout%read('/r', rall)
!    call hout%read('/theta', thetaall)
!    call hout%read('/phi', phiall)
!
!  else
!    !! 2D with swapped axes
!
!    call hout%read('/x2', x%x3all)
!    call hout%read('/x2i', x%x3iall)
!    call hout%read('/dx2b', x%dx3all)
!    call hout%read('/dx2h', x%dx3iall)
!    !! for a 3D grid this is x2, but now considered x3(all)
!    call hout%read('/x3', x%x2all)
!    call hout%read('/x3i', x%x2iall)
!    call hout%read('/dx3b', x%dx2all)
!    call hout%read('/dx3h', x%dx2iall)
!    !! formerly x3, now x2
!
!    block
!    !> NOTE: workaround for Intel 2020, may not really be a bug
!    !> notice this is the only one with negative indices
!    !real(wp), dimension(-1:lx1+2,-1:lx3all+2,-1:lx2all+2) :: htmp
!    real(wp), allocatable :: htmp(:,:,:)
!    allocate(htmp(-1:lx1+2,-1:lx3all+2,-1:lx2all+2))
!    !! end workaround
!
!    call hout%read('/h1', htmp)
!    x%h1all = reshape(htmp, [lx1+4,lx2all+4,lx3all+4], order=[1,3,2])
!    call hout%read('/h2', htmp)
!    !! this would be h3, but with the input structure shape
!    x%h3all = reshape(htmp, [lx1+4,lx2all+4,lx3all+4], order=[1,3,2])
!    !! permute the dimensions of the array 3 --> 2, 2 --> 3
!    call hout%read('/h3', htmp)
!    !! this would be h3, but with the input structure shape
!    x%h2all = reshape(htmp, [lx1+4,lx2all+4,lx3all+4], order=[1,3,2])
!    end block
!
!    block
!    !! htmp includes degenerate dimension
!    real(wp), dimension(1:lx1+1,1:lx3all, 1) :: htmp
!    if(lx3all==1) error stop 'lx3 is assumed to be degenerate'
!    !! input 2 vs. 3 dimensions swapped from this program
!    call hout%read('/h1x1i', htmp)
!    x%h1x1iall = reshape(htmp,[lx1+1,lx2all,lx3all],order=[1,3,2])
!    call hout%read('/h2x1i', htmp)
!    x%h3x1iall = reshape(htmp,[lx1+1,lx2all,lx3all],order=[1,3,2])
!    call hout%read('/h3x1i', htmp)
!    x%h2x1iall = reshape(htmp,[lx1+1,lx2all,lx3all],order=[1,3,2])
!    end block
!
!    block
!    !! htmp includes degenerate dimension
!    real(wp), dimension(1:lx1,1:lx3all+1, 1) :: htmp
!    call hout%read('/h1x2i', htmp)
!    x%h1x3iall = reshape(htmp,[lx1,lx2all,lx3all+1],order=[1,3,2])
!    !! Note also that the x2 interface from teh input file is x3i in this simulation
!    call hout%read('/h2x2i', htmp)
!    x%h3x3iall = reshape(htmp,[lx1,lx2all,lx3all+1],order=[1,3,2])
!    call hout%read('/h3x2i', htmp)
!    x%h2x3iall = reshape(htmp,[lx1,lx2all,lx3all+1],order=[1,3,2])
!    end block
!
!    block
!    real(wp), dimension(1:lx1,1:lx3all,1:lx2all+1) :: htmp
!    call hout%read('/h1x3i', htmp)
!    x%h1x2iall = reshape(htmp,[lx1,lx2all+1,lx3all],order=[1,3,2])
!    call hout%read('/h2x3i', htmp)
!    x%h3x2iall = reshape(htmp,[lx1,lx2all+1,lx3all],order=[1,3,2])
!    call hout%read('/h3x3i', htmp)
!    x%h2x2iall = reshape(htmp,[lx1,lx2all+1,lx3all],order=[1,3,2])
!    end block
!
!    block
!    !! htmp includes degenerate dimension
!    real(wp), dimension(lx1, lx3all, 1) :: htmp
!    call hout%read('/gx1', htmp)
!    g1all = reshape(htmp,[lx1,lx2all,lx3all],order=[1,3,2])
!    call hout%read('/gx2', htmp)
!    g3all = reshape(htmp,[lx1,lx2all,lx3all],order=[1,3,2])
!    call hout%read('/gx3', htmp)
!    g2all = reshape(htmp,[lx1,lx2all,lx3all],order=[1,3,2])
!
!    call hout%read('/alt', htmp)
!    x%altall = reshape(htmp,[lx1,lx2all,lx3all],order=[1,3,2])
!    call hout%read('/glat', htmp)
!    glatall = reshape(htmp,[lx1,lx2all,lx3all],order=[1,3,2])
!    call hout%read('/glon', htmp)
!    x%glonall = reshape(htmp,[lx1,lx2all,lx3all],order=[1,3,2])
!
!    call hout%read('/Bmag', htmp)
!    x%Bmagall = reshape(htmp,[lx1,lx2all,lx3all],order=[1,3,2])
!    end block
!
!    block
!    !! htmp includes degenerate dimension
!    real(wp), dimension(lx3all, 1) :: htmp
!    call hout%read('/I', htmp)
!    Incall = reshape(htmp,[lx2all,lx3all],order=[2,1])
!    end block
!
!    block
!    !! htmp includes degenerate dimension
!    real(wp), dimension(lx1, lx3all, 1) :: htmp
!    call hout%read('/nullpts', htmp)
!    nullptsall = reshape(htmp,[lx1,lx2all,lx3all],order=[1,3,2])
!    end block
!
!    block
!    real(wp), dimension(lx1,lx3all,lx2all,3) :: htmp
!    call hout%read('/e1', htmp)
!    e1all = reshape(htmp,[lx1,lx2all,lx3all,3],order=[1,3,2,4])
!
!    !> swap the x2/x3 unit vectors
!    call hout%read('/e2', htmp)
!    e3all = reshape(htmp,[lx1,lx2all,lx3all,3],order=[1,3,2,4])
!    call hout%read('/e3', htmp)
!    e2all = reshape(htmp,[lx1,lx2all,lx3all,3],order=[1,3,2,4])
!    call hout%read('/er', htmp)
!    erall = reshape(htmp,[lx1,lx2all,lx3all,3],order=[1,3,2,4])
!    call hout%read('/etheta', htmp)
!    ethetaall = reshape(htmp,[lx1,lx2all,lx3all,3],order=[1,3,2,4])
!    call hout%read('/ephi', htmp)
!    ephiall = reshape(htmp,[lx1,lx2all,lx3all,3],order=[1,3,2,4])
!    end block
!
!    block
!    !! htmp includes degenerate dimension
!    real(wp), dimension(lx1, lx3all, 1) :: htmp
!    call hout%read('/r', htmp)
!    rall = reshape(htmp,[lx1,lx2all,lx3all],order=[1,3,2])
!    call hout%read('/theta', htmp)
!    thetaall = reshape(htmp,[lx1,lx2all,lx3all],order=[1,3,2])
!    call hout%read('/phi', htmp)
!    phiall = reshape(htmp,[lx1,lx2all,lx3all],order=[1,3,2])
!    end block
!  endif
!
!  call hout%close()
!
!  x%rall = rall
!  x%thetaall = thetaall
!  x%phiall = phiall
!
!end procedure get_grid3_hdf5
!
end submodule readgrid_hdf5


text:
      subroutine getorb (iz, ihole, ion, norb, norbco,
     1                  den, nqn, nk, nel)

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

c     Save internal variables in case this gets re-entered
      save

c     Gets orbital data for chosen element.  Input is iz, atomic number
c     of desired element, other arguments are output.

c     Written by Steven Zabinsky, July 1989
c
c     last modified (20 aug 1989)  table increased to at no 95

c     Table for each element has occupation of the various levels.
c     The order of the levels in each array is:

c     element  level     principal qn (nqn), kappa qn (nk)
c           1  1s        1  -1
c           2  2s        2  -1
c           3  2p1/2     2   1
c           4  2p3/2     2  -2
c           5  3s        3  -1
c           6  3p1/2     3   1
c           7  3p3/2     3  -2
c           8  3d3/2     3   2
c           9  3d5/2     3  -3
c          10  4s        4  -1
c          11  4p1/2     4   1
c          12  4p3/2     4  -2
c          13  4d3/2     4   2
c          14  4d5/2     4  -3
c          15  4f5/2     4   3
c          16  4f7/2     4  -4
c          17  5s        5  -1
c          18  5p1/2     5   1
c          19  5p3/2     5  -2
c          20  5d3/2     5   2
c          21  5d5/2     5  -3
c          22  5f5/2     5   3
c          23  5f7/2     5  -4
c          24  6s        6  -1
c          25  6p1/2     6   1
c          26  6p3/2     6  -2
c          27  6d3/2     6   2
c          28  6d5/2     6  -3
c          29  7s        7  -1

      dimension den(30), nqn(30), nk(30), nel(30)
      dimension kappa (29)
      dimension iocc (95, 29)
      dimension nnum (29)
c     dimension ncore(95)

c     kappa quantum number for each orbital
c     k = - (j + 1/2)  if l = j - 1/2
c     k = + (j + 1/2)  if l = j + 1/2
      data kappa /-1,-1, 1,-2,-1,   1,-2, 2,-3,-1,   1,-2, 2,-3, 3,
     1            -4,-1, 1,-2, 2,  -3, 3,-4,-1, 1,  -2, 2,-3,-1/

c     principal quantum number (energy eigenvalue)
      data nnum  /1,2,2,2,3,  3,3,3,3,4,  4,4,4,4,4,
     1            4,5,5,5,5,  5,5,5,6,6,  6,6,6,7/

c     number of core orbitals for z = 1 to 95
c     data ncore
c    1  /0, 0, 1, 1, 1,  1, 1, 1, 1, 1,  4, 4, 4, 4, 4,  4, 4, 4, 4, 4,
c    2   4, 4, 4, 4, 4,  4, 4, 4, 9, 9,  9, 9, 9, 9, 9,  9, 9, 9, 9, 9,
c    3   9, 9, 9, 9, 9,  9, 9, 9, 9, 9,  9, 9, 9, 9, 9,  9, 9, 9, 9, 9,
c    4   9, 9, 9, 9, 9,  9, 9, 9, 9, 9, 16,16,16,16,16, 16,16,16,16,16,
c    5  16,16,16,16,16, 16,16,16,16,16, 16,16,16,16,16/

c     occupation of each level for z = 1, 95
      data (iocc( 1,i),i=1,29)  /1,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,
     1                           0,0,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc( 2,i),i=1,29)  /2,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,
     1                           0,0,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc( 3,i),i=1,29)  /2,1,0,0,0,  0,0,0,0,0,  0,0,0,0,0,
     1                           0,0,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc( 4,i),i=1,29)  /2,2,0,0,0,  0,0,0,0,0,  0,0,0,0,0,
     1                           0,0,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc( 5,i),i=1,29)  /2,2,1,0,0,  0,0,0,0,0,  0,0,0,0,0,
     1                           0,0,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc( 6,i),i=1,29)  /2,2,2,0,0,  0,0,0,0,0,  0,0,0,0,0,
     1                           0,0,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc( 7,i),i=1,29)  /2,2,2,1,0,  0,0,0,0,0,  0,0,0,0,0,
     1                           0,0,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc( 8,i),i=1,29)  /2,2,2,2,0,  0,0,0,0,0,  0,0,0,0,0,
     1                           0,0,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc( 9,i),i=1,29)  /2,2,2,3,0,  0,0,0,0,0,  0,0,0,0,0,
     1                           0,0,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(10,i),i=1,29)  /2,2,2,4,0,  0,0,0,0,0,  0,0,0,0,0,
     1                           0,0,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(11,i),i=1,29)  /2,2,2,4,1,  0,0,0,0,0,  0,0,0,0,0,
     1                           0,0,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(12,i),i=1,29)  /2,2,2,4,2,  0,0,0,0,0,  0,0,0,0,0,
     1                           0,0,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(13,i),i=1,29)  /2,2,2,4,2,  1,0,0,0,0,  0,0,0,0,0,
     1                           0,0,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(14,i),i=1,29)  /2,2,2,4,2,  2,0,0,0,0,  0,0,0,0,0,
     1                           0,0,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(15,i),i=1,29)  /2,2,2,4,2,  2,1,0,0,0,  0,0,0,0,0,
     1                           0,0,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(16,i),i=1,29)  /2,2,2,4,2,  2,2,0,0,0,  0,0,0,0,0,
     1                           0,0,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(17,i),i=1,29)  /2,2,2,4,2,  2,3,0,0,0,  0,0,0,0,0,
     1                           0,0,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(18,i),i=1,29)  /2,2,2,4,2,  2,4,0,0,0,  0,0,0,0,0,
     1                           0,0,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(19,i),i=1,29)  /2,2,2,4,2,  2,4,0,0,1,  0,0,0,0,0,
     1                           0,0,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(20,i),i=1,29)  /2,2,2,4,2,  2,4,0,0,2,  0,0,0,0,0,
     1                           0,0,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(21,i),i=1,29)  /2,2,2,4,2,  2,4,1,0,2,  0,0,0,0,0,
     1                           0,0,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(22,i),i=1,29)  /2,2,2,4,2,  2,4,2,0,2,  0,0,0,0,0,
     1                           0,0,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(23,i),i=1,29)  /2,2,2,4,2,  2,4,3,0,2,  0,0,0,0,0,
     1                           0,0,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(24,i),i=1,29)  /2,2,2,4,2,  2,4,4,1,1,  0,0,0,0,0,
     1                           0,0,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(25,i),i=1,29)  /2,2,2,4,2,  2,4,4,1,2,  0,0,0,0,0,
     1                           0,0,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(26,i),i=1,29)  /2,2,2,4,2,  2,4,4,2,2,  0,0,0,0,0,
     1                           0,0,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(27,i),i=1,29)  /2,2,2,4,2,  2,4,4,3,2,  0,0,0,0,0,
     1                           0,0,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(28,i),i=1,29)  /2,2,2,4,2,  2,4,4,4,2,  0,0,0,0,0,
     1                           0,0,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(29,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,1,  0,0,0,0,0,
     1                           0,0,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(30,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  0,0,0,0,0,
     1                           0,0,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(31,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  1,0,0,0,0,
     1                           0,0,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(32,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,0,0,0,0,
     1                           0,0,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(33,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,1,0,0,0,
     1                           0,0,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(34,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,2,0,0,0,
     1                           0,0,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(35,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,3,0,0,0,
     1                           0,0,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(36,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,0,0,0,
     1                           0,0,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(37,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,0,0,0,
     1                           0,1,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(38,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,0,0,0,
     1                           0,2,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(39,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,1,0,0,
     1                           0,2,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(40,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,2,0,0,
     1                           0,2,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(41,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,0,0,
     1                           0,1,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(42,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,1,0,
     1                           0,1,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(43,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,1,0,
     1                           0,2,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(44,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,3,0,
     1                           0,1,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(45,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,4,0,
     1                           0,1,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(46,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,0,
     1                           0,0,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(47,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,0,
     1                           0,1,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(48,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,0,
     1                           0,2,0,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(49,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,0,
     1                           0,2,1,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(50,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,0,
     1                           0,2,2,0,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(51,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,0,
     1                           0,2,2,1,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(52,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,0,
     1                           0,2,2,2,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(53,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,0,
     1                           0,2,2,3,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(54,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,0,
     1                           0,2,2,4,0,  0,0,0,0,0,  0,0,0,0/
      data (iocc(55,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,0,
     1                           0,2,2,4,0,  0,0,0,1,0,  0,0,0,0/
      data (iocc(56,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,0,
     1                           0,2,2,4,0,  0,0,0,2,0,  0,0,0,0/
      data (iocc(57,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,0,
     1                           0,2,2,4,1,  0,0,0,2,0,  0,0,0,0/
      data (iocc(58,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,2,
     1                           0,2,2,4,0,  0,0,0,2,0,  0,0,0,0/
      data (iocc(59,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,3,
     1                           0,2,2,4,0,  0,0,0,2,0,  0,0,0,0/
      data (iocc(60,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,4,
     1                           0,2,2,4,0,  0,0,0,2,0,  0,0,0,0/
      data (iocc(61,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,5,
     1                           0,2,2,4,0,  0,0,0,2,0,  0,0,0,0/
      data (iocc(62,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,6,
     1                           0,2,2,4,0,  0,0,0,2,0,  0,0,0,0/
      data (iocc(63,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,6,
     1                           1,2,2,4,0,  0,0,0,2,0,  0,0,0,0/
      data (iocc(64,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,6,
     1                           1,2,2,4,1,  0,0,0,2,0,  0,0,0,0/
      data (iocc(65,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,6,
     1                           3,2,2,4,0,  0,0,0,2,0,  0,0,0,0/
      data (iocc(66,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,6,
     1                           4,2,2,4,0,  0,0,0,2,0,  0,0,0,0/
      data (iocc(67,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,6,
     1                           5,2,2,4,0,  0,0,0,2,0,  0,0,0,0/
      data (iocc(68,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,6,
     1                           6,2,2,4,0,  0,0,0,2,0,  0,0,0,0/
      data (iocc(69,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,6,
     1                           7,2,2,4,0,  0,0,0,2,0,  0,0,0,0/
      data (iocc(70,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,6,
     1                           8,2,2,4,0,  0,0,0,2,0,  0,0,0,0/
      data (iocc(71,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,6,
     1                           8,2,2,4,1,  0,0,0,2,0,  0,0,0,0/
      data (iocc(72,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,6,
     1                           8,2,2,4,2,  0,0,0,2,0,  0,0,0,0/
      data (iocc(73,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,6,
     1                           8,2,2,4,3,  0,0,0,2,0,  0,0,0,0/
      data (iocc(74,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,6,
     1                           8,2,2,4,4,  0,0,0,2,0,  0,0,0,0/
      data (iocc(75,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,6,
     1                           8,2,2,4,4,  1,0,0,2,0,  0,0,0,0/
      data (iocc(76,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,6,
     1                           8,2,2,4,4,  2,0,0,2,0,  0,0,0,0/
      data (iocc(77,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,6,
     1                           8,2,2,4,4,  3,0,0,2,0,  0,0,0,0/
      data (iocc(78,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,6,
     1                           8,2,2,4,4,  5,0,0,1,0,  0,0,0,0/
      data (iocc(79,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,6,
     1                           8,2,2,4,4,  6,0,0,1,0,  0,0,0,0/
      data (iocc(80,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,6,
     1                           8,2,2,4,4,  6,0,0,2,0,  0,0,0,0/
      data (iocc(81,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,6,
     1                           8,2,2,4,4,  6,0,0,2,1,  0,0,0,0/
      data (iocc(82,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,6,
     1                           8,2,2,4,4,  6,0,0,2,2,  0,0,0,0/
      data (iocc(83,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,6,
     1                           8,2,2,4,4,  6,0,0,2,2,  1,0,0,0/
      data (iocc(84,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,6,
     1                           8,2,2,4,4,  6,0,0,2,2,  2,0,0,0/
      data (iocc(85,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,6,
     1                           8,2,2,4,4,  6,0,0,2,2,  3,0,0,0/
      data (iocc(86,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,6,
     1                           8,2,2,4,4,  6,0,0,2,2,  4,0,0,0/
      data (iocc(87,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,6,
     1                           8,2,2,4,4,  6,0,0,2,2,  4,0,0,1/
      data (iocc(88,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,6,
     1                           8,2,2,4,4,  6,0,0,2,2,  4,0,0,2/
      data (iocc(89,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,6,
     1                           8,2,2,4,4,  6,0,0,2,2,  4,1,0,2/
      data (iocc(90,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,6,
     1                           8,2,2,4,4,  6,0,0,2,2,  4,2,0,2/
      data (iocc(91,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,6,
     1                           8,2,2,4,4,  6,2,0,2,2,  4,1,0,2/
      data (iocc(92,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,6,
     1                           8,2,2,4,4,  6,3,0,2,2,  4,1,0,2/
      data (iocc(93,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,6,
     1                           8,2,2,4,4,  6,4,0,2,2,  4,1,0,2/
      data (iocc(94,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,6,
     1                           8,2,2,4,4,  6,6,0,2,2,  4,0,0,2/
      data (iocc(95,i),i=1,29)  /2,2,2,4,2,  2,4,4,6,2,  2,4,4,6,6,
     1                           8,2,2,4,4,  6,6,1,2,2,  4,0,0,2/

      if (iz .lt. 1  .or.  iz .gt. 95)  then
         call echo('   Atomic number not available')
         call fstop(' at GETORB')
      endif

      index = iz - ion
      if (ihole .gt. 0)  then
         index = index + 1
c        remove an electron from the level specified by ihole
         if (iocc(index,ihole) .lt. 1)  then
            call echo(' cannot remove an electron from this level')
            call fstop(' at GETORB-1')
         endif
         iocc(index,ihole) = iocc(index,ihole) - 1
      endif

      norb = 0
      do 10  i = 1, 29
         if (iocc(index,i) .ne. 0)  then
            norb = norb + 1
            nqn(norb) = nnum(i)
            nk(norb)  = kappa(i)
            nel(norb) = iocc(index,i)
            den(norb) = 0.0
         endif
   10 continue

c     restore iocc array for neatness
      if (ihole .gt. 0)  then
         iocc(index,ihole) = iocc(index,ihole) + 1
      endif

      norbco = norb

      return
      end

max_stars_repo_name: keechul/ifeffit

max_stars_repo_licenses: Naumen Condor-1.1 MS-PL

max_issues_repo_name: bruceravel/ifeffit

max_issues_repo_licenses: Naumen Condor-1.1 MS-PL

max_forks_repo_name: bruceravel/ifeffit

max_forks_repo_licenses: Naumen Condor-1.1 MS-PL


        !COMPILER-GENERATED INTERFACE MODULE: Mon Jul 19 15:11:37 2021
        ! This source file is for reference only and may not completely
        ! represent the generated interface used by the compiler.
        MODULE CDFUNC__genmod
          INTERFACE
            RECURSIVE FUNCTION CDFUNC(DEPTH,FLOW)
              REAL(KIND=8) :: DEPTH
              REAL(KIND=8) :: FLOW
              REAL(KIND=8) :: CDFUNC
            END FUNCTION CDFUNC
          END INTERFACE
        END MODULE CDFUNC__genmod


text:
      SUBROUTINE MB03XP( JOB, COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB,
     $                   Q, LDQ, Z, LDZ, ALPHAR, ALPHAI, BETA, DWORK,
     $                   LDWORK, INFO )
C
C     SLICOT RELEASE 5.5.
C
C     Copyright (c) 2002-2012 NICONET e.V.
C
C     PURPOSE
C
C     To compute the periodic Schur decomposition and the eigenvalues of
C     a product of matrices, H = A*B, with A upper Hessenberg and B
C     upper triangular without evaluating any part of the product.
C     Specifically, the matrices Q and Z are computed, so that
C
C          Q' * A * Z = S,    Z' * B * Q = T
C
C     where S is in real Schur form, and T is upper triangular.
C
C     ARGUMENTS
C
C     Mode Parameters
C
C     JOB     CHARACTER*1
C             Indicates whether the user wishes to compute the full
C             Schur form or the eigenvalues only, as follows:
C             = 'E':  Compute the eigenvalues only;
C             = 'S':  compute the factors S and T of the full
C                     Schur form.
C
C     COMPQ   CHARACTER*1
C             Indicates whether or not the user wishes to accumulate
C             the matrix Q as follows:
C             = 'N':  The matrix Q is not required;
C             = 'I':  Q is initialized to the unit matrix and the
C                     orthogonal transformation matrix Q is returned;
C             = 'V':  Q must contain an orthogonal matrix U on entry,
C                     and the product U*Q is returned.
C
C     COMPZ   CHARACTER*1
C             Indicates whether or not the user wishes to accumulate
C             the matrix Z as follows:
C             = 'N':  The matrix Z is not required;
C             = 'I':  Z is initialized to the unit matrix and the
C                     orthogonal transformation matrix Z is returned;
C             = 'V':  Z must contain an orthogonal matrix U on entry,
C                     and the product U*Z is returned.
C
C     Input/Output Parameters
C
C     N       (input) INTEGER
C             The order of the matrices A and B. N >= 0.
C
C     ILO     (input) INTEGER
C     IHI     (input) INTEGER
C             It is assumed that the matrices A and B are already upper
C             triangular in rows and columns 1:ILO-1 and IHI+1:N.
C             The routine works primarily with the submatrices in rows
C             and columns ILO to IHI, but applies the transformations to
C             all the rows and columns of the matrices A and B, if
C             JOB = 'S'.
C             1 <= ILO <= max(1,N+1); min(ILO,N) <= IHI <= N.
C
C     A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
C             On entry, the leading N-by-N part of this array A must
C             contain the upper Hessenberg matrix A.
C             On exit, if JOB = 'S', the leading N-by-N part of this
C             array is upper quasi-triangular with any 2-by-2 diagonal
C             blocks corresponding to a pair of complex conjugated
C             eigenvalues.
C             If JOB = 'E', the diagonal elements and 2-by-2 diagonal
C             blocks of A will be correct, but the remaining parts of A
C             are unspecified on exit.
C
C     LDA     INTEGER
C             The leading dimension of the array A.  LDA >= MAX(1,N).
C
C     B       (input/output) DOUBLE PRECISION array, dimension (LDB,N)
C             On entry, the leading N-by-N part of this array B must
C             contain the upper triangular matrix B.
C             On exit, if JOB = 'S', the leading N-by-N part of this
C             array contains the transformed upper triangular matrix.
C             2-by-2 blocks in B corresponding to 2-by-2 blocks in A
C             will be reduced to positive diagonal form. (I.e., if
C             A(j+1,j) is non-zero, then B(j+1,j)=B(j,j+1)=0 and B(j,j)
C             and B(j+1,j+1) will be positive.)
C             If JOB = 'E', the elements corresponding to diagonal
C             elements and 2-by-2 diagonal blocks in A will be correct,
C             but the remaining parts of B are unspecified on exit.
C
C     LDB     INTEGER
C             The leading dimension of the array B.  LDB >= MAX(1,N).
C
C     Q       (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
C             On entry, if COMPQ = 'V', then the leading N-by-N part of
C             this array must contain a matrix Q which is assumed to be
C             equal to the unit matrix except for the submatrix
C             Q(ILO:IHI,ILO:IHI).
C             If COMPQ = 'I', Q need not be set on entry.
C             On exit, if COMPQ = 'V' or COMPQ = 'I' the leading N-by-N
C             part of this array contains the transformation matrix
C             which produced the Schur form.
C             If COMPQ = 'N', Q is not referenced.
C
C     LDQ     INTEGER
C             The leading dimension of the array Q.  LDQ >= 1.
C             If COMPQ <> 'N', LDQ >= MAX(1,N).
C
C     Z       (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
C             On entry, if COMPZ = 'V', then the leading N-by-N part of
C             this array must contain a matrix Z which is assumed to be
C             equal to the unit matrix except for the submatrix
C             Z(ILO:IHI,ILO:IHI).
C             If COMPZ = 'I', Z need not be set on entry.
C             On exit, if COMPZ = 'V' or COMPZ = 'I' the leading N-by-N
C             part of this array contains the transformation matrix
C             which produced the Schur form.
C             If COMPZ = 'N', Z is not referenced.
C
C     LDZ     INTEGER
C             The leading dimension of the array Z.  LDZ >= 1.
C             If COMPZ <> 'N', LDZ >= MAX(1,N).
C
C     ALPHAR  (output) DOUBLE PRECISION array, dimension (N)
C     ALPHAI  (output) DOUBLE PRECISION array, dimension (N)
C     BETA    (output) DOUBLE PRECISION array, dimension (N)
C             The i-th (1 <= i <= N) computed eigenvalue is given by
C             BETA(I) * ( ALPHAR(I) + sqrt(-1)*ALPHAI(I) ). If two
C             eigenvalues are computed as a complex conjugate pair,
C             they are stored in consecutive elements of ALPHAR, ALPHAI
C             and BETA. If JOB = 'S', the eigenvalues are stored in the
C             same order as on the diagonales of the Schur forms of A
C             and B.
C
C     Workspace
C
C     DWORK   DOUBLE PRECISION array, dimension (LDWORK)
C             On exit, if INFO = 0,  DWORK(1)  returns the optimal
C             value of LDWORK.
C             On exit, if  INFO = -19,  DWORK(1)  returns the minimum
C             value of LDWORK.
C
C     LDWORK  INTEGER
C             The length of the array DWORK.  LDWORK >= MAX(1,N).
C
C     Error Indicator
C
C     INFO    INTEGER
C             = 0:  successful exit;
C             < 0:  if INFO = -i, the i-th argument had an illegal
C                   value;
C             > 0:  if INFO = i, then MB03XP failed to compute the Schur
C                   form in a total of 30*(IHI-ILO+1) iterations;
C                   elements 1:ilo-1 and i+1:n of ALPHAR, ALPHAI and
C                   BETA contain successfully computed eigenvalues.
C
C     METHOD
C
C     The implemented algorithm is a multi-shift version of the periodic
C     QR algorithm described in [1,3] with some minor modifications
C     proposed in [2].
C
C     REFERENCES
C
C     [1] Bojanczyk, A.W., Golub, G.H., and Van Dooren, P.
C         The periodic Schur decomposition: Algorithms and applications.
C         Proc. of the SPIE Conference (F.T. Luk, Ed.), 1770, pp. 31-42,
C         1992.
C
C     [2] Kressner, D.
C         An efficient and reliable implementation of the periodic QZ
C         algorithm. Proc. of the IFAC Workshop on Periodic Control
C         Systems, pp. 187-192, 2001.
C
C     [3] Van Loan, C.
C         Generalized Singular Values with Algorithms and Applications.
C         Ph. D. Thesis, University of Michigan, 1973.
C
C     NUMERICAL ASPECTS
C
C     The algorithm requires O(N**3) floating point operations and is
C     backward stable.
C
C     CONTRIBUTORS
C
C     D. Kressner, Technical Univ. Berlin, Germany, and
C     P. Benner, Technical Univ. Chemnitz, Germany, December 2003.
C
C     REVISIONS
C
C     V. Sima, May 2008 (SLICOT version of the HAPACK routine DHGPQR).
C     V. Sima, Oct. 2011.
C
C     KEYWORDS
C
C     Eigenvalue, eigenvalue decomposition, Hessenberg form, orthogonal
C     transformation, (periodic) Schur form
C
C     ******************************************************************
C
C     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE
      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
      INTEGER            NSMAX, LDAS, LDBS
      PARAMETER          ( NSMAX = 15, LDAS = NSMAX, LDBS = NSMAX )
C     .. Scalar Arguments ..
      CHARACTER          COMPQ, COMPZ, JOB
      INTEGER            IHI, ILO, INFO, LDA, LDB, LDQ, LDWORK, LDZ, N
C     .. Array Arguments ..
      DOUBLE PRECISION   A(LDA,*), ALPHAI(*), ALPHAR(*), B(LDB,*),
     $                   BETA(*), DWORK(*), Q(LDQ,*), Z(LDZ,*)
C     .. Local Scalars ..
      LOGICAL            INITQ, INITZ, WANTQ, WANTT, WANTZ
      INTEGER            DUM, I, I1, I2, IERR, ITEMP, ITN, ITS, J, K,
     $                   KK, L, MAXB, NH, NR, NS, NV, PV2, PV3
      DOUBLE PRECISION   OVFL, SMLNUM, TAUV, TAUW, TEMP, TST, ULP, UNFL
C     .. Local Arrays ..
      INTEGER            ISEED(4)
      DOUBLE PRECISION   AS(LDAS,LDAS), BS(LDBS,LDBS), V(3*NSMAX+6)
C     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            IDAMAX, UE01MD
      DOUBLE PRECISION   DLAMCH, DLANHS
      EXTERNAL           DLAMCH, DLANHS, IDAMAX, LSAME, UE01MD
C     .. External Subroutines ..
      EXTERNAL           DAXPY, DCOPY, DGEMV, DLABAD, DLACPY, DLARFG,
     $                   DLARFX, DLARNV, DLASET, DSCAL, DTRMV, MB03YA,
     $                   MB03YD, XERBLA
C     .. Intrinsic Functions ..
      INTRINSIC          ABS, DBLE, MAX, MIN
C
C     .. Executable Statements ..
C
C     Decode the scalar input parameters.
C
      WANTT = LSAME( JOB, 'S' )
      INITQ = LSAME( COMPQ, 'I' )
      WANTQ = INITQ.OR.LSAME( COMPQ, 'V' )
      INITZ = LSAME( COMPZ, 'I' )
      WANTZ = INITZ.OR.LSAME( COMPZ, 'V' )
C
C     Check the scalar input parameters.
C
      INFO = 0
      IF ( .NOT.LSAME( JOB, 'E' ) .AND. .NOT.WANTT ) THEN
         INFO = -1
      ELSE IF ( .NOT.LSAME( COMPQ, 'N' ) .AND. .NOT.WANTQ ) THEN
         INFO = -2
      ELSE IF ( .NOT.LSAME( COMPZ, 'N' ) .AND. .NOT.WANTZ ) THEN
         INFO = -3
      ELSE IF ( N.LT.0 ) THEN
         INFO = -4
      ELSE IF ( ILO.LT.1 .OR. ILO.GT.MAX( 1, N+1 ) ) THEN
         INFO = -5
      ELSE IF ( IHI.LT.MIN( ILO,N ).OR.IHI.GT.N ) THEN
         INFO = -6
      ELSE IF ( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -8
      ELSE IF ( LDB.LT.MAX( 1, N ) ) THEN
         INFO = -10
      ELSE IF ( LDQ.LT.1 .OR. WANTQ .AND. LDQ.LT.N ) THEN
         INFO = -12
      ELSE IF ( LDZ.LT.1 .OR. WANTZ .AND. LDZ.LT.N ) THEN
         INFO = -14
      ELSE IF ( LDWORK.LT.MAX( 1, N ) ) THEN
         DWORK(1) = DBLE( MAX( 1, N ) )
         INFO = -19
      END IF
C
C     Return if there were illegal values.
C
      IF ( INFO.NE.0 ) THEN
         CALL XERBLA( 'MB03XP', -INFO )
         RETURN
      END IF
C
C     Initialize Q and Z, if necessary.
C
      IF ( INITQ )
     $   CALL DLASET( 'All', N, N, ZERO, ONE, Q, LDQ )
      IF ( INITZ )
     $   CALL DLASET( 'All', N, N, ZERO, ONE, Z, LDZ )
C
C     Store isolated eigenvalues and standardize B.
C
C     FOR I = [1:ILO-1, IHI+1:N]
      I = 1
   10 CONTINUE
      IF ( I.EQ.ILO ) THEN
         I = IHI+1
      END IF
      IF ( I.LE.N ) THEN
         IF ( B(I,I).LT.ZERO ) THEN
            IF ( WANTT ) THEN
               DO 20  K = ILO, I
                  B(K,I) = -B(K,I)
   20          CONTINUE
               DO 30  K = I, IHI
                  A(I,K) = -A(I,K)
   30          CONTINUE
            ELSE
               B(I,I) = -B(I,I)
               A(I,I) = -A(I,I)
            END IF
            IF ( WANTQ ) THEN
               DO 40  K = ILO, IHI
                  Q(K,I) = -Q(K,I)
   40          CONTINUE
            END IF
         END IF
         ALPHAR(I) = A(I,I)
         ALPHAI(I) = ZERO
         BETA(I) = B(I,I)
         I = I + 1
C        END FOR
         GO TO 10
      END IF
C
C     Quick return if possible.
C
      IF ( N.EQ.0 .OR. ILO.EQ.IHI+1 ) THEN
         DWORK(1) = ONE
         RETURN
      END IF
C
C     Set rows and coloms ILO to IHI of B (A) to zero below the first
C     (sub)diagonal.
C
      DO  60 J = ILO, IHI - 2
         DO 50  I = J + 2, N
            A(I,J) = ZERO
   50    CONTINUE
   60 CONTINUE
      DO  80 J = ILO, IHI - 1
         DO 70  I = J + 1, N
            B(I,J) = ZERO
   70    CONTINUE
   80 CONTINUE
      NH = IHI - ILO + 1
C
C     Suboptimal choice of the number of shifts.
C
      IF ( WANTQ ) THEN
         NS   = UE01MD( 4, 'MB03XP', JOB // COMPQ, N, ILO, IHI )
         MAXB = UE01MD( 8, 'MB03XP', JOB // COMPQ, N, ILO, IHI )
      ELSE
         NS   = UE01MD( 4, 'MB03XP', JOB // COMPZ, N, ILO, IHI )
         MAXB = UE01MD( 8, 'MB03XP', JOB // COMPZ, N, ILO, IHI )
      END IF
C
      IF ( NS.LE.2 .OR. NS.GT.NH .OR. MAXB.GE.NH ) THEN
C
C        Standard double-shift product QR.
C
         CALL MB03YD( WANTT, WANTQ, WANTZ, N, ILO, IHI, ILO, IHI, A,
     $                LDA, B, LDB, Q, LDQ, Z, LDZ, ALPHAR, ALPHAI, BETA,
     $                DWORK, LDWORK, INFO )
         RETURN
      END IF
      MAXB = MAX( 3, MAXB )
      NS = MIN( NS, MAXB, NSMAX )
C
C     Set machine-dependent constants for the stopping criterion.
C     If max(norm(A),norm(B)) <= sqrt(OVFL), then overflow should not
C     occur.
C
      UNFL = DLAMCH( 'Safe minimum' )
      OVFL = ONE / UNFL
      CALL DLABAD( UNFL, OVFL )
      ULP = DLAMCH( 'Precision' )
      SMLNUM = UNFL*( DBLE( NH ) / ULP )
C
C     I1 and I2 are the indices of the first rows and last columns of
C     A and B to which transformations must be applied.
C
      IF ( WANTT ) THEN
         I1 = 1
         I2 = N
      END IF
      ISEED(1) = 1
      ISEED(2) = 0
      ISEED(3) = 0
      ISEED(4) = 1
C
C     ITN is the maximal number of QR iterations.
C
      ITN = 30*NH
      DUM = 0
C
C     Main loop. Eigenvalues I+1:IHI have converged. Either L = ILO
C     or A(L,L-1) is negligible.
C
      I = IHI
   90 CONTINUE
      L = ILO
      IF ( I.LT.ILO )
     $   GO TO 210
C
      DO 190  ITS = 0, ITN
         DUM = DUM + (IHI-ILO)*(IHI-ILO)
C
C        Look for deflations in A.
C
         DO 100  K = I, L + 1, -1
            TST = ABS( A(K-1,K-1) ) + ABS( A(K,K) )
            IF ( TST.EQ.ZERO )
     $         TST = DLANHS( '1', I-L+1, A(L,L), LDA, DWORK )
            IF ( ABS( A(K,K-1) ).LE.MAX( ULP*TST, SMLNUM ) )
     $         GO TO 110
  100    CONTINUE
  110    CONTINUE
C
C        Look for deflation in B if problem size is greater than 1.
C
         IF ( I-K.GE.1 ) THEN
            DO 120  KK = I, K, -1
               IF ( KK.EQ.I ) THEN
                  TST = ABS( B(KK-1,KK) )
               ELSE IF ( KK.EQ.K ) THEN
                  TST = ABS( B(KK,KK+1) )
               ELSE
                  TST = ABS( B(KK-1,KK) ) + ABS( B(KK,KK+1) )
               END IF
               IF ( TST.EQ.ZERO )
     $            TST = DLANHS( '1', I-K+1, B(K,K), LDB, DWORK )
               IF ( ABS( B(KK,KK) ).LE.MAX( ULP*TST, SMLNUM ) )
     $            GO TO 130
  120       CONTINUE
         ELSE
            KK = K-1
         END IF
  130    CONTINUE
         IF ( KK.GE.K ) THEN
C
C           B has an element close to zero at position (KK,KK).
C
            B(KK,KK) = ZERO
            CALL MB03YA( WANTT, WANTQ, WANTZ, N, K, I, ILO, IHI, KK,
     $                   A, LDA, B, LDB, Q, LDQ, Z, LDZ, INFO )
            K = KK+1
         END IF
         L = K
         IF( L.GT.ILO ) THEN
C
C           A(L,L-1) is negligible.
C
            A(L,L-1) = ZERO
         END IF
C
C        Exit from loop if a submatrix of order <= MAXB has split off.
C
         IF ( L.GE.I-MAXB+1 )
     $      GO TO 200
C
C        The active submatrices are now in rows and columns L:I.
C
         IF ( .NOT.WANTT ) THEN
            I1 = L
            I2 = I
         END IF
         IF ( ITS.EQ.10.OR.ITS.EQ.20 ) THEN
C
C           Exceptional shift. The first column of the shift polynomial
C           is a pseudo-random vector.
C
            CALL DLARNV( 3, ISEED, NS+1, V )
         ELSE
C
C           Use eigenvalues of trailing submatrix as shifts.
C
            CALL DLACPY( 'Full', NS, NS, A(I-NS+1,I-NS+1), LDA, AS,
     $                   LDAS )
            CALL DLACPY( 'Full', NS, NS, B(I-NS+1,I-NS+1), LDB, BS,
     $                   LDBS )
            CALL MB03YD( .FALSE., .FALSE., .FALSE., NS, 1, NS, 1, NS,
     $                   AS, LDAS, BS, LDBS, Q, LDQ, Z, LDZ,
     $                   ALPHAR(I-NS+1), ALPHAI(I-NS+1), BETA(I-NS+1),
     $                   DWORK, LDWORK, IERR )
         END IF
C
C        Compute the nonzero elements of the first column of
C        (A*B-w(1)) (A*B-w(2)) .. (A*B-w(ns)).
C
         V(1) = ONE
         NV = 1
C        WHILE NV <= NS
  140    CONTINUE
         IF ( NV.LE.NS ) THEN
            IF ( NV.EQ.NS .OR. AS(NV+1,NV).EQ.ZERO ) THEN
C
C              Real shift.
C
               V(NV+1) = ZERO
               PV2 = NV+2
               CALL DCOPY( NV, V, 1, V(PV2), 1 )
               CALL DTRMV( 'Upper', 'No transpose', 'No unit diagonal',
     $                     NV, B(L,L), LDB, V(PV2), 1 )
               CALL DSCAL( NV, BS(NV,NV), V, 1 )
               ITEMP = IDAMAX( 2*NV+1, V, 1 )
               TEMP = ONE / MAX( ABS( V(ITEMP) ), SMLNUM )
               CALL DSCAL( 2*NV+1, TEMP, V, 1 )
               CALL DGEMV( 'No transpose', NV+1, NV, ONE, A(L,L), LDA,
     $                     V(PV2), 1, -AS(NV,NV), V, 1 )
               NV = NV + 1
            ELSE
C
C              Double shift using a product formulation of the shift
C              polynomial [2].
C
               V(NV+1) = ZERO
               V(NV+2) = ZERO
               PV2 = NV+3
               PV3 = 2*NV+5
               CALL DCOPY( NV+2, V, 1, V(PV2), 1 )
               CALL DCOPY( NV+1, V, 1, V(PV3), 1 )
               CALL DSCAL( NV, BS(NV+1,NV+1), V(PV2), 1 )
               CALL DTRMV( 'Upper', 'No transpose', 'No unit diagonal',
     $                     NV, B(L,L), LDB, V(PV3), 1 )
               ITEMP = IDAMAX( 2*NV+3, V(PV2), 1 )
               TEMP = ONE / MAX( ABS( V(PV2+ITEMP-1) ), SMLNUM )
               CALL DSCAL( 2*NV+3, TEMP, V(PV2), 1 )
C
               CALL DCOPY( NV, V(PV2), 1, V, 1 )
               CALL DGEMV( 'No transpose', NV+1, NV, -ONE, A(L,L), LDA,
     $                     V(PV3), 1, AS(NV+1,NV+1), V(PV2), 1 )
               CALL DSCAL( NV, AS(NV,NV+1), V, 1 )
               ITEMP = IDAMAX( 2*NV+3, V, 1 )
               TEMP = ONE / MAX( ABS( V(ITEMP) ), SMLNUM )
               CALL DSCAL( 2*NV+3, TEMP, V, 1 )
C
               CALL DSCAL( NV, -AS(NV+1,NV), V, 1 )
               CALL DAXPY( NV+1, AS(NV,NV), V(PV2), 1, V, 1)
               ITEMP = IDAMAX( 2*NV+3, V, 1 )
               TEMP = ONE / MAX( ABS( V(ITEMP) ), SMLNUM )
               CALL DSCAL( 2*NV+3, TEMP, V, 1 )
C
               CALL DSCAL( NV+1, BS(NV,NV), V, 1 )
               CALL DTRMV( 'Upper', 'No transpose', 'No unit diagonal',
     $                     NV+1, B(L,L), LDB, V(PV2), 1 )
               ITEMP = IDAMAX( 2*NV+3, V, 1 )
               TEMP = ONE / MAX( ABS( V(ITEMP) ), SMLNUM )
               CALL DSCAL( 2*NV+3, TEMP, V, 1 )
C
               CALL DGEMV( 'No transpose', NV+2, NV+1, -ONE, A(L,L),
     $                     LDA, V(PV2), 1, ONE, V, 1 )
               NV = NV + 2
            END IF
            ITEMP = IDAMAX( NV, V, 1 )
            TEMP = ABS( V(ITEMP) )
            IF ( TEMP.EQ.ZERO ) THEN
               V(1) = ONE
               DO 150  K = 2, NV
                  V(K) = ZERO
  150          CONTINUE
            ELSE
               TEMP = MAX( TEMP, SMLNUM )
               CALL DSCAL( NV, ONE/TEMP, V, 1 )
            END IF
            GO TO 140
C        END WHILE
         END IF
C
C        Multi-shift product QR step.
C
         PV2 = NS+2
         DO 180  K = L,I-1
            NR = MIN( NS+1,I-K+1 )
            IF ( K.GT.L )
     $         CALL DCOPY( NR, A(K,K-1), 1, V, 1 )
            CALL DLARFG( NR, V(1), V(2), 1, TAUV )
            IF ( K.GT.L ) THEN
               A(K,K-1) = V(1)
               DO 160  KK = K+1,I
                  A(KK,K-1) = ZERO
  160          CONTINUE
            END IF
C
C           Apply reflector V from the right to B in rows
C           I1:min(K+NS,I).
C
            V(1) = ONE
            CALL DLARFX( 'Right', MIN(K+NS,I)-I1+1, NR, V, TAUV,
     $                   B(I1,K), LDB, DWORK )
C
C           Annihilate the introduced nonzeros in the K-th column.
C
            CALL DCOPY( NR, B(K,K), 1, V(PV2), 1 )
            CALL DLARFG( NR, V(PV2), V(PV2+1), 1, TAUW )
            B(K,K) = V(PV2)
            DO 170  KK = K+1,I
               B(KK,K) = ZERO
  170       CONTINUE
            V(PV2) = ONE
C
C           Apply reflector W from the left to transform the rows of the
C           matrix B in columns K+1:I2.
C
            CALL DLARFX( 'Left', NR, I2-K, V(PV2), TAUW, B(K,K+1), LDB,
     $                   DWORK )
C
C           Apply reflector V from the left to transform the rows of the
C           matrix A in columns K:I2.
C
            CALL DLARFX( 'Left', NR, I2-K+1, V, TAUV, A(K,K), LDA,
     $                   DWORK )
C
C           Apply reflector W from the right to transform the columns of
C           the matrix A in rows I1:min(K+NS,I).
C
            CALL DLARFX( 'Right', MIN(K+NS+1,I)-I1+1, NR, V(PV2), TAUW,
     $                   A(I1,K), LDA, DWORK )
C
C           Accumulate transformations in the matrices Q and Z.
C
            IF ( WANTQ )
     $         CALL DLARFX( 'Right', NH, NR, V, TAUV, Q(ILO,K), LDQ,
     $                      DWORK )
            IF ( WANTZ )
     $         CALL DLARFX( 'Right', NH, NR, V(PV2), TAUW, Z(ILO,K),
     $                      LDZ, DWORK )
  180    CONTINUE
  190 CONTINUE
C
C     Failure to converge.
C
      INFO = I
      RETURN
  200 CONTINUE
C
C     Submatrix of order <= MAXB has split off. Use double-shift
C     periodic QR algorithm.
C
      CALL MB03YD( WANTT, WANTQ, WANTZ, N, L, I, ILO, IHI, A, LDA, B,
     $             LDB, Q, LDQ, Z, LDZ, ALPHAR, ALPHAI, BETA, DWORK,
     $             LDWORK, INFO )
      IF ( INFO.GT.0 )
     $   RETURN
      ITN = ITN - ITS
      I = L - 1
      GO TO 90
C
  210 CONTINUE
      DWORK(1) = DBLE( MAX( 1,N ) )
      RETURN
C *** Last line of MB03XP ***
      END

max_stars_repo_name: bgin/MissileSimulation

max_issues_repo_name: bgin/MissileSimulation

max_forks_repo_name: bgin/MissileSimulation


text:
+incdir+${ROCKET_SOC}/ve/rocket_soc_uvm/tb
+incdir+${ROCKET_SOC}/ve/rocket_soc_uvm/tests

-f ${SV_BFMS}/uart/uvm/uvm.f

-f ${VMON}/src/client/sv/vmon_sv_client_hvl.f
-f ${VMON}/src/client/sv/vmon_sv_client_uvm.f
${VMON}/src/client/sv/wb_vmon_monitor.sv

${ROCKET_SOC}/ve/rocket_soc_uvm/tb/RocketSocTBEnvBasePkg.sv
${ROCKET_SOC}/ve/rocket_soc_uvm/tb/rocket_soc_env_pkg.sv
${ROCKET_SOC}/ve/rocket_soc_uvm/tests/rocket_soc_tests_pkg.sv
${ROCKET_SOC}/ve/rocket_soc_uvm/tb/RocketSocTBHvl.sv


10 program main
      implicit none
      integer :: a,b
      a = 1
      a = b - 1
      a = b * 1
      a = b / 1
      a = b < 1_4
20 end program main