text:
SUBMODULE (Misc) Misc_implementation
    USE Precision, ONLY : i4k, i8k, r8k
    IMPLICIT NONE

    CONTAINS

        MODULE PROCEDURE interpret_string
        IMPLICIT NONE
        !! Interprets a string (typically read from an input file) into a user-defined # of character and/or integer inputs
        INTEGER(i4k) :: i
        CHARACTER(LEN=:), ALLOCATABLE :: string, sep, char
        TYPE :: counter
            INTEGER(i4k) :: t = 0, i = 0, r = 0, c = 0
        END TYPE counter
        TYPE (counter) :: cnt

        IF (PRESENT(ignore)) THEN
            string = TRIM(ADJUSTL(line(INDEX(line,ignore)+LEN(ignore):)))
        ELSE
            string = TRIM(ADJUSTL(line))
        END IF
        IF (PRESENT(separator)) THEN
            ALLOCATE(sep, source=separator)
        ELSE
            ALLOCATE(sep, source=' ')
        END IF
        IF (PRESENT(ints)) THEN
            IF (ALLOCATED(ints)) DEALLOCATE(ints)
            ALLOCATE(ints(1:SIZE(datatype)),source=0_i4k)
        END IF
        IF (PRESENT(reals)) THEN
            IF (ALLOCATED(reals)) DEALLOCATE(reals)
            ALLOCATE(reals(1:SIZE(datatype)),source=0.0_r8k)
        END IF
        IF (PRESENT(chars)) THEN
            IF (ALLOCATED(chars)) DEALLOCATE(chars)
            ALLOCATE(chars(1:SIZE(datatype)))
        END IF

        DO i = 1, SIZE(datatype)
            SELECT CASE (datatype(i))
            CASE ('I', 'i')
                !! Integer
                cnt%i = cnt%i + 1
                CALL get_string_value (string, sep, ints(cnt%i))
            CASE ('R', 'r')
                !! Real
                cnt%r = cnt%r + 1
                CALL get_string_value (string, sep, reals(cnt%r))
            CASE ('C', 'c')
                !! Character
                cnt%c = cnt%c + 1
                CALL get_string_value (string, sep, char)
                ALLOCATE(chars(cnt%c)%text, source=char)
            END SELECT
            CALL reduce_string (string, sep)
            cnt%t = cnt%t + 1
        END DO

        line = string

        END PROCEDURE interpret_string

        MODULE PROCEDURE reduce_string
        IMPLICIT NONE

        IF (INDEX(string,sep) == 0) THEN
            string = ''
        ELSE
            string = ADJUSTL(string(INDEX(string,sep)+LEN(sep):))
        END IF

        END PROCEDURE reduce_string

        MODULE PROCEDURE get_string_char
        IMPLICIT NONE

        IF (INDEX(string,sep) == 0) THEN
            name = string(1:)                    !! Read to end of string
        ELSE
            name = string(1:INDEX(string,sep)-1) !! Read until sep is found
        END IF

        END PROCEDURE get_string_char

        MODULE PROCEDURE get_string_int
        IMPLICIT NONE

        CHARACTER(LEN=:), ALLOCATABLE :: text

        IF (INDEX(string,sep) == 0) THEN
            ALLOCATE(text, source=string(1:))                    !! Read to end of string
        ELSE
            ALLOCATE(text, source=string(1:INDEX(string,sep)-1)) !! Read until sep is found
        END IF
        READ(text,'(i8)') name                                   !! Store value

        END PROCEDURE get_string_int

        MODULE PROCEDURE get_string_real
        IMPLICIT NONE
        CHARACTER(LEN=:), ALLOCATABLE :: text

        IF (INDEX(string,sep) == 0) THEN
            ALLOCATE(text, source=string(1:))                    !! Read to end of string
        ELSE
            ALLOCATE(text, source=string(1:INDEX(string,sep)-1)) !! Read until sep is found
        END IF
        READ(text,'(es13.6)') name                               !! Store value

        END PROCEDURE get_string_real

        MODULE PROCEDURE convert_real32_to_string
        IMPLICIT NONE
        !! Converts a real to a character string
        CHARACTER(LEN=20) :: tmp_string = '                    '

        WRITE(tmp_string,*) var
        ALLOCATE(string,source=TRIM(ADJUSTL(tmp_string)))

        END PROCEDURE convert_real32_to_string

        MODULE PROCEDURE convert_real64_to_string
        IMPLICIT NONE
        !! Converts a real to a character string
        CHARACTER(LEN=30) :: tmp_string = '                              '

        WRITE(tmp_string,*) var
        ALLOCATE(string,source=TRIM(ADJUSTL(tmp_string)))

        END PROCEDURE convert_real64_to_string

        MODULE PROCEDURE convert_real64_array_to_string
        IMPLICIT NONE
        !! Converts a real to a character string
        INTEGER(i4k) :: i
        CHARACTER(LEN=:), ALLOCATABLE :: tmp_string

        DO i = 1, SIZE(var)
            ALLOCATE(tmp_string, source=convert_real64_to_string(var(i)))
            IF (.NOT. ALLOCATED(string)) THEN
                ALLOCATE(string,source=tmp_string)
            ELSE
                string = string // " " // tmp_string
            END IF
            DEALLOCATE(tmp_string)
        END DO

        END PROCEDURE convert_real64_array_to_string

        MODULE PROCEDURE convert_int32_to_string
        IMPLICIT NONE
        CHARACTER(LEN=20) :: tmp_string = '                    '

        WRITE(tmp_string,*) var
        ALLOCATE(string,source=TRIM(ADJUSTL(tmp_string)))

        END PROCEDURE convert_int32_to_string

        MODULE PROCEDURE convert_int64_to_string
        IMPLICIT NONE
        CHARACTER(LEN=30) :: tmp_string = '                              '

        WRITE(tmp_string,*) var
        ALLOCATE(string,source=TRIM(ADJUSTL(tmp_string)))

        END PROCEDURE convert_int64_to_string

        MODULE PROCEDURE convert_logical_to_string
        IMPLICIT NONE

        IF (var) THEN
            ALLOCATE(string,source='True')
        ELSE
            ALLOCATE(string,source='False')
        END IF

        END PROCEDURE convert_logical_to_string

        MODULE PROCEDURE to_uppercase
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 01/23/2019
        !!
        !! This function changes lowercase text in a string to uppercase text
        !!
        INTEGER(i4k) :: i, j
        CHARACTER(LEN=26), PARAMETER    :: CAPL = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
        CHARACTER(LEN=26), PARAMETER    :: LOWL = 'abcdefghijklmnopqrstuvwxyz'

        new_string = string(1:LEN_TRIM(string))

        DO i = 1, LEN_TRIM(string)
            j = INDEX(LOWL, string(i:i))
            IF (j > 0) THEN
                new_string(i:i) = CAPL(j:j)
            ELSE
                new_string(i:i) = string(i:i)
            END IF
        END DO

        END PROCEDURE to_uppercase

        MODULE PROCEDURE to_lowercase
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 01/23/2019
        !!
        !! This function changes uppercase text in a string to lowercase text
        !!
        INTEGER(i4k) :: i, j
        CHARACTER(LEN=26), PARAMETER    :: CAPL = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
        CHARACTER(LEN=26), PARAMETER    :: LOWL = 'abcdefghijklmnopqrstuvwxyz'

        new_string = string(1:LEN_TRIM(string))

        DO i = 1, LEN_TRIM(string)
            j = INDEX(CAPL, string(i:i))
            IF (j > 0) THEN
                new_string(i:i) = LOWL(j:j)
            ELSE
                new_string(i:i) = string(i:i)
            END IF
        END DO

        END PROCEDURE to_lowercase

        MODULE PROCEDURE sleep_for
        IMPLICIT NONE
        !! author: Zaak Beekman, ParaTools
        !! date: 8/8/2018
        !!
        !! This performs a 'sleep' for a specified amount of time
        !!
        INTEGER(i4k), DIMENSION(8) :: time
        INTEGER(i8k) :: ms_t1, ms_t2, msecs_big

        CALL DATE_AND_TIME(values=time)

        ms_t1=(time(5)*3600+time(6)*60+time(7))*1000+time(8)
        msecs_big = msecs

        DO !! spin until elapsed time is greater than msecs
            CALL DATE_AND_TIME(values=time)
            ms_t2=(time(5)*3600+time(6)*60+time(7))*1000+time(8)
            IF ( ms_t2 - ms_t1 >= msecs_big ) EXIT
        END DO

        END PROCEDURE

END SUBMODULE Misc_implementation

max_stars_repo_name: zbeekman/vtkmofo

max_issues_repo_name: zbeekman/vtkmofo

max_forks_repo_name: zbeekman/vtkmofo


! This file is part of toml-f.
!
! Copyright (C) 2019-2020 Sebastian Ehlert
!
! Licensed under either of Apache License, Version 2.0 or MIT license
! at your option; you may not use this file except in compliance with
! the License.
!
! Unless required by applicable law or agreed to in writing, software
! distributed under the License is distributed on an "AS IS" BASIS,
! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
! See the License for the specific language governing permissions and
! limitations under the License.

!> Abstraction layer for the actual storage of the data structure.
!>
!> The structure implementations provide the actual storage for TOML values, with
!> a generic enough interface to make the definition of the TOML data structures
!> independent of the actual algorithm used for storing the TOML values.
!>
!> Every data structure defined here should strive to only use allocatable
!> data types and limit the use of pointer attributes as they interfer with
!> the automatic memory management of Fortran. A well defined data structure
!> in allocatables allows deep-copying of TOML values by assignment, data structures
!> requiring pointer attributes have to define an assignment(=) interface to
!> allow deep-copying of TOML values.
module tomlf_structure
   use tomlf_structure_base, only : toml_structure, toml_ordered
   use tomlf_structure_vector, only : toml_vector, new_vector
   implicit none
   private

   public :: toml_structure, toml_ordered
   public :: new_structure, new_ordered
   public :: len

   !> Overload len function
   interface len
      module procedure :: get_len
   end interface

contains

!> Constructor for the storage data structure
subroutine new_structure(self)

   !> Instance of the structure
   class(toml_structure), allocatable, intent(out) :: self

   type(toml_vector), allocatable :: vect

   allocate(vect)
   call new_vector(vect)
   call move_alloc(vect, self)

end subroutine new_structure

!> Constructor for the ordered storage data structure
subroutine new_ordered(self)

   !> Instance of the structure
   class(toml_ordered), allocatable, intent(out) :: self

   type(toml_vector), allocatable :: vect

   allocate(vect)
   call new_vector(vect)
   call move_alloc(vect, self)

end subroutine new_ordered

!> Get number of TOML values in the structure
pure function get_len(self) result(length)

   !> Instance of the structure
   class(toml_ordered), intent(in) :: self

   !> Current length of the ordered structure
   integer :: length

   length = self%get_len()

end function get_len

end module tomlf_structure


text:
PROGRAM PIB1D_Hydrogen_Radial
    IMPLICIT NONE

    REAL, PARAMETER :: PI = acos(-1.0)
    REAL, ALLOCATABLE:: potential(:), x(:)
    DOUBLE PRECISION, ALLOCATABLE:: hmat(:,:), eigenVal(:), work(:)
    DOUBLE PRECISION:: WF
    REAL:: xmax, xmin, dx, hbar, tconst, m, l
    INTEGER:: nx, lwork, i, ifail

    WRITE(*, fmt = '(/A)', ADVANCE = 'NO') "Enter the number of grid points: "
    READ(*, *) nx

    xmax = 10.0
    xmin = 0.0
    dx = (xmax - xmin)/(nx)

    ALLOCATE(potential(nx), hmat(nx,nx), eigenVal(nx), x(nx), work(64*nx))
    DO i = 1, nx
        x(i) = xmin + (i)*dx
    ENDDO

    m = 1.0
    hbar = 1.00
    l = 0
    tconst = (hbar**2) / (2.0 * m * dx**2)

    DO i = 1, nx
        potential(i) = (hbar**2 * l * (l+1))/(2 * m * x(i) * x(i)) - 1/x(i)
    ENDDO

    hmat = 0.0
    DO i = 1, nx-1
        hmat(i, i) = potential(i) + 2.0 * tconst
        hmat(i, i+1) = -tconst
        hmat(i+1, i) = -tconst
    ENDDO
    hmat(nx, nx) = potential(nx) + 2.0 * tconst

    lwork = 64*nx
    CALL dsyev('V', 'U', nx, hmat, nx, eigenVal, work, lwork, ifail)

    OPEN(UNIT = 20, FILE = 'hydrogen.txt')
        DO i = 1, nx
            WF = hmat(i, 3)/sqrt(dx)
            WRITE(20, 10) x(i), potential(i), eigenVal(i), hmat(i, 1)/sqrt(dx)
            10 FORMAT(f0.6, 3x, f0.6, 3x, f0.8, 3x, f0.8, 3x)
        ENDDO
    CLOSE(20)

END PROGRAM PIB1D_Hydrogen_Radial

max_stars_repo_path: Particle In a Box (PIB)/4_hydrogen.f90

max_issues_repo_path: Particle In a Box (PIB)/4_hydrogen.f90

max_forks_repo_path: Particle In a Box (PIB)/4_hydrogen.f90


text:
!=======================================================================
!  COPYRIGHT 1998-2020 DSSAT Foundation
!                      University of Florida, Gainesville, Florida
!                      Inernational Fertilizer Development Center
!
!  ALL RIGHTS RESERVED
!=======================================================================
!=======================================================================
!  MAKEFILEW, Subroutine, Fabio Olivera, Willingthon Pavan, Gerrit Hoogenboom
!
! This subroutine open the FILEX and read the weather station characters (WSTA).
! Then create the weather filename and check if FILEW exists.
!-----------------------------------------------------------------------
!  Revision history
!
!  01/23/2020 FO  Written
!  11/10/2020 CHP Need to make revisions here like in IPEXP for forecast mode
!       This mode needs to keep FILEW (WTH file) and potentially also
!       FILEWC (CLI) file and FILEWG (WTG file).
!  10/21/2021 FO  Updated GETLUN for weather file unit
!-----------------------------------------------------------------------
!  INPUT  : DSSATP, PATHEX, FILEX
!
!  OUTPUT : YEAR
!-----------------------------------------------------------------------
!=======================================================================

      SUBROUTINE MAKEFILEW(LUNEXP,DSSATP,PATHEX,FILEX,        &
                         SimLevel,LNSIM,LNPLT,LNFLD)               !Input - Factor Levels

      USE ModuleDefs

      IMPLICIT      NONE

      CHARACTER*1   BLANK, WTHER
      CHARACTER*3   PROCOD
      CHARACTER*6   ERRKEY,SECTION
      CHARACTER*5   FINDH
      CHARACTER*8   WSTA,FILEW4,FILEX
      CHARACTER*12  FILEW,NAMEF, LastFileW
      CHARACTER*78  MSG(4)
      CHARACTER*80  PATHWT,PATHEX,CHARTEST
      CHARACTER*92  FILEWW,FILETMP
      CHARACTER*102 DSSATP
      CHARACTER*255 LINE

      LOGICAL       FEXIST,SimLevel

      INTEGER       LUNEXP,LINEXP,LNUM,ISECT,ERRNUM,PATHL,IFIND,FOUND
      INTEGER       ICASAF,LUNWTH,LINWTH,UNIT
      INTEGER       SDATE,PDATE,YRDOY
      INTEGER       LNFLD,TLNFLD
      INTEGER       LNSIM,TLNSIM
      INTEGER       LNPLT,TLNPLT
      INTEGER       YEARDOY,YEAR,YR,DOY

      DATA LastFileW /" "/

      PARAMETER (ERRKEY = 'MAKEFW')
      PARAMETER (BLANK = ' ')

      SDATE = -99
      PDATE = -99
      WTHER = 'M'

!-----------------------------------------------------------------------
!    Open FILEX and read Weather station characters (WSTA)
!-----------------------------------------------------------------------

      LNUM   = 0
      LINEXP = 0
      !FO - Rewind because ipexp.for call FileX LUNEXP as a constant
      REWIND(LUNEXP)

      IF(LNFLD .GT. 0) THEN
        !FO - Find FIELD section to search for weather station
        SECTION ='*FIELD'
        CALL FIND (LUNEXP,SECTION,LINEXP,IFIND)
        IF(IFIND .NE. 0) THEN

          !FO - Find the specific section header
          FINDH ='L ID_'
          CALL HFIND(LUNEXP,FINDH,LINEXP,IFIND)
          IF(IFIND .EQ. 1) THEN
            !FO - Loop through good lines to find and read the weather station
            DO
              CALL IGNORE2(LUNEXP,LINEXP,ISECT,CHARTEST)

              IF(ISECT .EQ. 1) THEN
                READ (CHARTEST,'(I3,9X,A8)', IOSTAT=ERRNUM) TLNFLD,WSTA
                IF (ERRNUM .NE. 0) CALL ERROR (ERRKEY,4,FILEX,LINEXP)
                IF (WSTA(1:1) .EQ. BLANK) CALL ERROR (ERRKEY,3,FILEX,LINEXP)
              ELSE
                CALL ERROR (ERRKEY,3,FILEX,LINEXP)
              ENDIF

              !Until FL factor lvl equals to the current section.
              IF (TLNFLD .EQ. LNFLD .AND. LNFLD .GT. 0) THEN
                EXIT
              ENDIF
            ENDDO
          ELSE
            !Error - If FINDH ='L ID_' not found
            CALL ERROR (ERRKEY,2,FILEX,LINEXP)
          ENDIF
        ELSE
          !Error - If SECTION ='*FIELD' not found
          CALL ERROR (ERRKEY,2,FILEX,LINEXP)
        ENDIF
      ELSE
        !Error - If LNFLD  = 0
        CALL ERROR (ERRKEY,1,FILEX,LINEXP)
      ENDIF

!-----------------------------------------------------------------------
!    Read FILEX YRSIM and weather method from *SIMUALTION CONTROLS
!    Otherwise, read planting date for experiments (LNSIM = 0)
!-----------------------------------------------------------------------

      LNUM   = 0
      LINEXP = 0
      !FO - Rewind FileX, because find always set LINEXP to 0
      REWIND(LUNEXP)

      !FO - Search for *SIMUL CONTROLS if SM factor level is gt 0
      IF(LNSIM .GT. 0 .AND. SimLevel) THEN

        !FO - Find *SIMULATION CONTROLS section to search for start simulation date
        !FO - DO/UNTIL loop is used to iterate until find the right section
        SECTION = '*SIMUL'
        CALL FIND (LUNEXP,SECTION,LINEXP,IFIND)
        LNUM = LNUM + LINEXP

        IF(IFIND .NE. 0) THEN
          DO
            !FO - Find the specific header to search for SDATE
            FINDH ='N GEN'
            CALL HFIND(LUNEXP,FINDH,LINEXP,IFIND)
            IF(IFIND .EQ. 1) THEN
              CALL IGNORE2(LUNEXP,LINEXP,ISECT,CHARTEST)
              IF(ISECT .EQ. 1) THEN
                READ (CHARTEST,'(I2,30X,I6)', IOSTAT=ERRNUM) TLNSIM,SDATE
                IF (ERRNUM .NE. 0) CALL ERROR (ERRKEY,8,FILEX,LNUM)
              ELSE
                CALL ERROR (ERRKEY,7,FILEX,LNUM)
              ENDIF
            ELSE
              !Error - If !FINDH = 'N GEN'
              CALL ERROR (ERRKEY,6,FILEX,LNUM)
            ENDIF

            !FO - Find the specific header to search for WTHER
            FINDH ='N MET'
            CALL HFIND(LUNEXP,FINDH,LINEXP,IFIND)
            IF(IFIND .EQ. 1) THEN

              CALL IGNORE2(LUNEXP,LINEXP,ISECT,CHARTEST)
              IF(ISECT .EQ. 1) THEN
                READ (CHARTEST,'(I2,17X,A1)', IOSTAT=ERRNUM) TLNSIM,WTHER
                IF (ERRNUM .NE. 0) CALL ERROR (ERRKEY,11,FILEX,LNUM)
                IF (WTHER .EQ. BLANK) CALL ERROR (ERRKEY,10,FILEX,LINEXP)
              ELSE
                CALL ERROR (ERRKEY,10,FILEX,LNUM)
              ENDIF
            ELSE
              !Error - If !FINDH ='N MET'
              CALL ERROR (ERRKEY,9,FILEX,LNUM)
            ENDIF

            !Until SM factor lvl equals to the current section.
            IF (TLNSIM .EQ. LNSIM .AND. LNSIM .GT. 0) THEN
              EXIT
            ENDIF
          ENDDO
        ELSE
          !Error - If LNSIM = 0
          CALL ERROR (ERRKEY,6,FILEX,LNUM)
        ENDIF

      !FO - Search for *PLANTING DETAILS if MP factor level is gt 0
      ELSE IF(LNPLT .GT. 0) THEN

        SECTION = '*PLANT'
        CALL FIND (LUNEXP,SECTION,LINEXP,IFIND)
        LNUM = LNUM + LINEXP

        IF(IFIND .NE. 0) THEN

          !FO - Find the specific header to search for PDATE
          FINDH ='P PDA'
          CALL HFIND(LUNEXP,FINDH,LINEXP,IFIND)
          IF(IFIND .EQ. 1) THEN
            DO
              CALL IGNORE2(LUNEXP,LINEXP,ISECT,CHARTEST)
              IF(ISECT .EQ. 1) THEN

                READ (CHARTEST,'(I2,I6)', IOSTAT=ERRNUM) TLNPLT,PDATE
                IF (ERRNUM .NE. 0) CALL ERROR (ERRKEY,14,FILEX,LNUM)

              ELSE
                !Error - If IGNORE2
                CALL ERROR (ERRKEY,13,FILEX,LNUM)
              ENDIF

              !Until SM factor lvl equals to the current section.
              IF (TLNPLT .EQ. LNPLT .AND. LNPLT .GT. 0) THEN
                EXIT
              ENDIF
            ENDDO
          ELSE
            !Error - If !FINDH ='P PDA'
            CALL ERROR (ERRKEY,12,FILEX,LNUM)
          ENDIF
        ELSE
          !Error - If !SECTION ='*PLANT'
          CALL ERROR (ERRKEY,12,FILEX,LNUM)
        ENDIF

      ELSE
        !Error - If both (SM, PLT) are not found
        CALL ERROR (ERRKEY,5,FILEX,LNUM)
      ENDIF

!-----------------------------------------------------------------------
!    Converts SDATE or PDATE to YR and DOY
!-----------------------------------------------------------------------
      IF(SDATE .GT. 0 ) THEN
        YRDOY = SDATE
      ELSEIF(PDATE .GT. 0) THEN
        YRDOY = PDATE
      ELSE
        CALL ERROR (ERRKEY,15,FILEX,LNUM)
      ENDIF

      !FO - Converts YRDOY to YR and DOY
      YR  = INT(YRDOY / 1000)
      DOY = YRDOY - YR * 1000

!-----------------------------------------------------------------------
!    Create WTH file name and path
!-----------------------------------------------------------------------
!    Establish the weather file FILEW as WSTA + .WT?  where ? :
!          M = observed data
!          G = generated data
!          S = interactively generated
!-----------------------------------------------------------------------
      IF(WTHER .EQ. 'M') THEN
        IF (WSTA(5:5) .EQ. BLANK) THEN
          WRITE (FILEW(1:12),'(A4,I2.2,A6)') WSTA,YR,'01.WTH'
        ELSE
          WRITE (FILEW(1:12),'(A8,A4)') WSTA,'.WTH'
        ENDIF
        PROCOD = 'WED'
      ELSEIF(WTHER .EQ. 'G') THEN

        IF (WSTA(5:5) .EQ. BLANK) THEN
          WRITE (FILEW(1:12),'(A4,I2.2,A6)') WSTA,YR,'01.WTG'
        ELSE
          WRITE (FILEW(1:12),'(A8,A4)') WSTA,'.WTG'
        ENDIF
        PROCOD = 'WGD'
      ELSEIF (WTHER .EQ. 'S' .OR. WTHER .EQ. 'W') THEN
         WRITE (FILEW(1:12),'(A4,A8)') WSTA,'.CLI    '
         PROCOD = 'CLD'
      ELSE
         CALL ERROR (ERRKEY,16,FILEX,LINEXP)
      ENDIF

!     Check weather filename in current directory
      INQUIRE (FILE = FILEW,EXIST = FEXIST)
      IF (FEXIST) THEN
        PATHWT = BLANK
!     Check weather filename in data directory
      ELSE
        FILETMP = TRIM(PATHEX)//FILEW
        INQUIRE (FILE = FILETMP,EXIST = FEXIST)
        IF (FEXIST) THEN
          PATHWT = TRIM(PATHEX)
!       Check weather filename in default DSSAT directory
        ELSE
          CALL PATH(PROCOD,DSSATP,PATHWT,1,NAMEF)
          FILETMP = TRIM(PATHWT) // FILEW
          INQUIRE (FILE=FILETMP, EXIST = FEXIST)
          IF (FEXIST) THEN
            PATHWT = PATHWT
!         Check 4-character file name in data directory
          ELSE
            FILEW4 = FILEW(1:4) // ".WTH"
            FILETMP = TRIM(PATHEX) // FILEW4
            INQUIRE (FILE=FILETMP, EXIST = FEXIST)
            IF (FEXIST) THEN
              PATHWT = TRIM(PATHEX)
              FILEW = FILEW4
!           Check 4-character filename in default DSSAT directory
            ELSE
              FILEW4 = FILEW(1:4) // ".WTH"
              FILETMP = TRIM(PATHWT) // FILEW4
              INQUIRE (FILE=FILETMP, EXIST = FEXIST)
              IF (FEXIST) THEN
                PATHWT = TRIM(PATHWT)
                FILEW = FILEW4
              ELSE
                MSG(1) = "Weather file not found."
                MSG(2) = "  Neither " // FILEW // " nor " // FILEW4
                MSG(3) = "  were found in weather or experiment directories."
                MSG(4) = "Simulation will end."
                CALL WARNING(4,ERRKEY,MSG)
                CALL ERROR(ERRKEY,18,FILEW,0)
              ENDIF
            ENDIF
          ENDIF
        ENDIF
      ENDIF

!-----------------------------------------------------------------------
!    Open and Read Weather file data
!-----------------------------------------------------------------------
      PATHL  = INDEX(PATHWT,BLANK)
      IF (PATHL .LE. 1) THEN
         FILEWW = FILEW
       ELSE
         FILEWW = PATHWT(1:(PATHL-1)) // FILEW
      ENDIF

      !FO - Process FILEW if last FILEW is different.
      IF(LastFileW /= FILEW) THEN

        LastFileW = FILEW

        ICASAF  = 0
        YEAR    = -99
        DOY     = -99
        YEARDOY = -99

        CALL GETLUN('LUNMFW', UNIT)

        INQUIRE(FILE=FILEWW, EXIST = FEXIST)
        IF(FEXIST) THEN
          OPEN (UNIT, FILE = FILEWW,STATUS = 'OLD',IOSTAT=ERRNUM)
          IF (ERRNUM .NE. 0) CALL ERROR (ERRKEY,18,FILEW,0)
        ELSE
          CALL ERROR (ERRKEY,18,FILEW,0)
        ENDIF

        DO WHILE (ERRNUM .EQ. 0)
          READ (UNIT,'(A)', IOSTAT=ERRNUM) LINE

          !FO - Look for 1st header line beginning with '$'
          IF(INDEX(LINE, '$WEATHER') .GT. 0 .AND. ICASAF .EQ. 0) THEN
            ICASAF = 1

          !FO - New Weather file format true read 7-digit. Otherwise keep the default.
          ELSEIF(ICASAF .EQ. 1) THEN

            !FO - Look for header lines beginning with '@'
            IF(LINE(1:1) .EQ. '@' .AND. INDEX(LINE, 'DATE') .GT. 0) THEN
                !FO - Read 7-digit First Weather YEAR
                READ (UNIT,'(I7)', IOSTAT=ERRNUM) YEARDOY
                IF (ERRNUM .NE. 0) CALL ERROR (ERRKEY,17,FILEW,0)

                EXIT
            ELSE IF(LINE(1:1) .EQ. '@' .AND. INDEX(LINE, 'YEAR') .GT. 0) THEN
                !FO - Read ICASA format
                READ (UNIT,'(I4,1X,I3)', IOSTAT=ERRNUM) YEAR, DOY
                IF (ERRNUM .NE. 0) CALL ERROR (ERRKEY,17,FILEW,0)

                YEARDOY = YEAR * 1000 + DOY

                EXIT
              ENDIF
            ENDIF

        ENDDO

        CLOSE(UNIT)

!-----------------------------------------------------------------------
!    Convert SDATE to 7-digits and check if is in the date range
!-----------------------------------------------------------------------

        IF (YEARDOY .GT. 0 .AND. ICASAF .EQ. 1) THEN
          FirstWeatherDate = YEARDOY

          YEARDOY = INT(FirstWeatherDate/100000) * 100000 + YRDOY

          IF(YEARDOY .GE. FirstWeatherDate) THEN
            NEWSDATE = YEARDOY
          ELSE
            NEWSDATE = INT((FirstWeatherDate+99000)/100000)*100000+YRDOY
          ENDIF

          ! Error Checking
          IF(NEWSDATE .LT. FirstWeatherDate .OR. NEWSDATE .GT. FirstWeatherDate+99000) THEN
            CALL ERROR (ERRKEY,19,FILEX,0)
          ENDIF

        ELSE
          MSG(1) = 'Y4KDOY - Weather file'
          MSG(2) = 'Not able to find 7-digits first weather date'
          MSG(3) = 'Set simulation to 5-digits first weather date'
          FirstWeatherDate = -99

          CALL INFO(3,ERRKEY,MSG)
        ENDIF

      ENDIF

    END SUBROUTINE MAKEFILEW

max_stars_repo_name: xiangyi-wang/dssat-csm-os

max_issues_repo_name: xiangyi-wang/dssat-csm-os

max_forks_repo_name: xiangyi-wang/dssat-csm-os


text:
MODULE common_mpi_speedy
!=======================================================================
!
! [PURPOSE:] MPI procedures
!
! [ATTENTION:]
!   DO NOT COMPILE WITH BOTH INLINE EXPANSION AND OMP OPTIONS TOGETHER
!    (Use ONE if you want, but DON'T USE BOTH AT THE SAME TIME)
!
! [HISTORY:]
!   01/23/2009 Takemasa Miyoshi  created
!
!=======================================================================
!$USE OMP_LIB
  USE common
  USE common_mpi
  USE common_speedy
  IMPLICIT NONE
  PUBLIC

  INTEGER,PARAMETER :: mpibufsize=1000
  INTEGER,SAVE :: nij1
  INTEGER,SAVE :: nij1max
  INTEGER,ALLOCATABLE,SAVE :: nij1node(:)
  REAL(r_size),ALLOCATABLE,SAVE :: phi1(:)
  REAL(r_size),ALLOCATABLE,SAVE :: dx1(:),dy1(:)
  REAL(r_size),ALLOCATABLE,SAVE :: lon1(:),lat1(:)

CONTAINS
SUBROUTINE set_common_mpi_speedy
  REAL(r_sngl) :: v3dg(nlon,nlat,nlev,nv3d)
  REAL(r_sngl) :: v2dg(nlon,nlat,nv2d)
  REAL(r_size),ALLOCATABLE :: v3d(:,:,:)
  REAL(r_size),ALLOCATABLE :: v2d(:,:)
  INTEGER :: i,n

  WRITE(6,'(A)') 'Hello from set_common_mpi_speedy'
  i = MOD(nlon*nlat,nprocs)
  nij1max = (nlon*nlat - i)/nprocs + 1
  IF(myrank < i) THEN
    nij1 = nij1max
  ELSE
    nij1 = nij1max - 1
  END IF
  WRITE(6,'(A,I3.3,A,I6)') 'MYRANK ',myrank,' number of grid points: nij1= ',nij1
  ALLOCATE(nij1node(nprocs))
  DO n=1,nprocs
    IF(n-1 < i) THEN
      nij1node(n) = nij1max
    ELSE
      nij1node(n) = nij1max - 1
    END IF
  END DO

  ALLOCATE(phi1(nij1))
  ALLOCATE(dx1(nij1))
  ALLOCATE(dy1(nij1))
  ALLOCATE(lon1(nij1))
  ALLOCATE(lat1(nij1))

  ALLOCATE(v3d(nij1,nlev,nv3d))
  ALLOCATE(v2d(nij1,nv2d))
  DO i=1,nlat
    v3dg(:,i,1,1) = SNGL(dx(i))
    v3dg(:,i,1,2) = SNGL(dy2(i))
    v3dg(:,i,1,3) = SNGL(lon(:))
    v3dg(:,i,1,4) = SNGL(lat(i))
  END DO
  v2dg(:,:,1) = SNGL(phi0(:,:))
  CALL scatter_grd_mpi(0,v3dg,v2dg,v3d,v2d)
  dx1(:) = v3d(:,1,1)
  dy1(:) = v3d(:,1,2)
  lon1(:) = v3d(:,1,3)
  lat1(:) = v3d(:,1,4)
  phi1 = v2d(:,1)

  RETURN
END SUBROUTINE set_common_mpi_speedy
!-----------------------------------------------------------------------
! Scatter gridded data to processes (nrank -> all)
!-----------------------------------------------------------------------
SUBROUTINE scatter_grd_mpi(nrank,v3dg,v2dg,v3d,v2d)
  INTEGER,INTENT(IN) :: nrank
  REAL(r_sngl),INTENT(IN) :: v3dg(nlon,nlat,nlev,nv3d)
  REAL(r_sngl),INTENT(IN) :: v2dg(nlon,nlat,nv2d)
  REAL(r_size),INTENT(OUT) :: v3d(nij1,nlev,nv3d)
  REAL(r_size),INTENT(OUT) :: v2d(nij1,nv2d)

  IF(mpibufsize > nij1max) THEN
    CALL scatter_grd_mpi_fast(nrank,v3dg,v2dg,v3d,v2d)
  ELSE
    CALL scatter_grd_mpi_safe(nrank,v3dg,v2dg,v3d,v2d)
  END IF

  RETURN
END SUBROUTINE scatter_grd_mpi

SUBROUTINE scatter_grd_mpi_1d(nrank,v1dg,v1d)
  INTEGER,INTENT(IN) :: nrank
  REAL(r_sngl),INTENT(IN) :: v1dg(nlon,nlat)
  REAL(r_size),INTENT(OUT) :: v1d(nij1)
  REAL(r_sngl) :: tmp(nij1max,nprocs)
  REAL(r_sngl) :: bufs(mpibufsize,nprocs)
  REAL(r_sngl) :: bufr(mpibufsize)
  INTEGER :: i,j,ierr,ns,nr
  INTEGER :: iter,niter

  ns = mpibufsize
  nr = ns
  niter = CEILING(REAL(nij1max)/REAL(mpibufsize))

  IF(myrank == nrank) CALL grd_to_buf(v1dg(:,:),tmp)
  DO iter=1,niter
    IF(myrank == nrank) THEN
      i = mpibufsize * (iter-1)
      DO j=1,mpibufsize
        i=i+1
        IF(i > nij1max) EXIT
        bufs(j,:) = tmp(i,:)
      END DO
    END IF
    CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
    CALL MPI_SCATTER(bufs,ns,MPI_REAL,&
                   & bufr,nr,MPI_REAL,nrank,MPI_COMM_WORLD,ierr)
    i = mpibufsize * (iter-1)
    DO j=1,mpibufsize
       i=i+1
      IF(i > nij1) EXIT
      v1d(i) = REAL(bufr(j),r_size)
    END DO
  END DO
  CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)

  RETURN
END SUBROUTINE scatter_grd_mpi_1d

SUBROUTINE scatter_grd_mpi_safe(nrank,v3dg,v2dg,v3d,v2d)
  INTEGER,INTENT(IN) :: nrank
  REAL(r_sngl),INTENT(IN) :: v3dg(nlon,nlat,nlev,nv3d)
  REAL(r_sngl),INTENT(IN) :: v2dg(nlon,nlat,nv2d)
  REAL(r_size),INTENT(OUT) :: v3d(nij1,nlev,nv3d)
  REAL(r_size),INTENT(OUT) :: v2d(nij1,nv2d)
  REAL(r_sngl) :: tmp(nij1max,nprocs)
  REAL(r_sngl) :: bufs(mpibufsize,nprocs)
  REAL(r_sngl) :: bufr(mpibufsize)
  INTEGER :: i,j,k,n,ierr,ns,nr
  INTEGER :: iter,niter

  ns = mpibufsize
  nr = ns
  niter = CEILING(REAL(nij1max)/REAL(mpibufsize))

  DO n=1,nv3d
    DO k=1,nlev
      IF(myrank == nrank) CALL grd_to_buf(v3dg(:,:,k,n),tmp)
      DO iter=1,niter
        IF(myrank == nrank) THEN
          i = mpibufsize * (iter-1)
          DO j=1,mpibufsize
            i=i+1
            IF(i > nij1max) EXIT
            bufs(j,:) = tmp(i,:)
          END DO
        END IF
        CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
        CALL MPI_SCATTER(bufs,ns,MPI_REAL,&
                       & bufr,nr,MPI_REAL,nrank,MPI_COMM_WORLD,ierr)
        i = mpibufsize * (iter-1)
        DO j=1,mpibufsize
          i=i+1
          IF(i > nij1) EXIT
          v3d(i,k,n) = REAL(bufr(j),r_size)
        END DO
      END DO
    END DO
  END DO

  DO n=1,nv2d
    IF(myrank == nrank) CALL grd_to_buf(v2dg(:,:,n),tmp)
    DO iter=1,niter
      IF(myrank == nrank) THEN
        i = mpibufsize * (iter-1)
        DO j=1,mpibufsize
          i=i+1
          IF(i > nij1max) EXIT
          bufs(j,:) = tmp(i,:)
        END DO
      END IF
      CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
      CALL MPI_SCATTER(bufs,ns,MPI_REAL,&
                     & bufr,nr,MPI_REAL,nrank,MPI_COMM_WORLD,ierr)
      i = mpibufsize * (iter-1)
      DO j=1,mpibufsize
        i=i+1
        IF(i > nij1) EXIT
        v2d(i,n) = REAL(bufr(j),r_size)
      END DO
    END DO
  END DO

  CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)

  RETURN
END SUBROUTINE scatter_grd_mpi_safe

SUBROUTINE scatter_grd_mpi_fast(nrank,v3dg,v2dg,v3d,v2d)
  INTEGER,INTENT(IN) :: nrank
  REAL(r_sngl),INTENT(IN) :: v3dg(nlon,nlat,nlev,nv3d)
  REAL(r_sngl),INTENT(IN) :: v2dg(nlon,nlat,nv2d)
  REAL(r_size),INTENT(OUT) :: v3d(nij1,nlev,nv3d)
  REAL(r_size),INTENT(OUT) :: v2d(nij1,nv2d)
  REAL(r_sngl) :: bufs(nij1max,nlevall,nprocs)
  REAL(r_sngl) :: bufr(nij1max,nlevall)
  INTEGER :: j,k,n,ierr,ns,nr

  ns = nij1max * nlevall
  nr = ns
  IF(myrank == nrank) THEN
    j=0
    DO n=1,nv3d
      DO k=1,nlev
        j = j+1
        CALL grd_to_buf(v3dg(:,:,k,n),bufs(:,j,:))
      END DO
    END DO

    DO n=1,nv2d
      j = j+1
      CALL grd_to_buf(v2dg(:,:,n),bufs(:,j,:))
    END DO
  END IF

  CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
  CALL MPI_SCATTER(bufs,ns,MPI_REAL,&
                 & bufr,nr,MPI_REAL,nrank,MPI_COMM_WORLD,ierr)

  j=0
  DO n=1,nv3d
    DO k=1,nlev
      j = j+1
      v3d(:,k,n) = REAL(bufr(1:nij1,j),r_size)
    END DO
  END DO

  DO n=1,nv2d
    j = j+1
    v2d(:,n) = REAL(bufr(1:nij1,j),r_size)
  END DO

  CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)

  RETURN
END SUBROUTINE scatter_grd_mpi_fast
!-----------------------------------------------------------------------
! Gather gridded data (all -> nrank)
!-----------------------------------------------------------------------
SUBROUTINE gather_grd_mpi(nrank,v3d,v2d,v3dg,v2dg)
  INTEGER,INTENT(IN) :: nrank
  REAL(r_size),INTENT(IN) :: v3d(nij1,nlev,nv3d)
  REAL(r_size),INTENT(IN) :: v2d(nij1,nv2d)
  REAL(r_sngl),INTENT(OUT) :: v3dg(nlon,nlat,nlev,nv3d)
  REAL(r_sngl),INTENT(OUT) :: v2dg(nlon,nlat,nv2d)

  IF(mpibufsize > nij1max) THEN
    CALL gather_grd_mpi_fast(nrank,v3d,v2d,v3dg,v2dg)
  ELSE
    CALL gather_grd_mpi_safe(nrank,v3d,v2d,v3dg,v2dg)
  END IF

  RETURN
END SUBROUTINE gather_grd_mpi

SUBROUTINE gather_grd_mpi_safe(nrank,v3d,v2d,v3dg,v2dg)
  INTEGER,INTENT(IN) :: nrank
  REAL(r_size),INTENT(IN) :: v3d(nij1,nlev,nv3d)
  REAL(r_size),INTENT(IN) :: v2d(nij1,nv2d)
  REAL(r_sngl),INTENT(OUT) :: v3dg(nlon,nlat,nlev,nv3d)
  REAL(r_sngl),INTENT(OUT) :: v2dg(nlon,nlat,nv2d)
  REAL(r_sngl) :: tmp(nij1max,nprocs)
  REAL(r_sngl) :: bufs(mpibufsize)
  REAL(r_sngl) :: bufr(mpibufsize,nprocs)
  INTEGER :: i,j,k,n,ierr,ns,nr
  INTEGER :: iter,niter

  ns = mpibufsize
  nr = ns
  niter = CEILING(REAL(nij1max)/REAL(mpibufsize))

  DO n=1,nv3d
    DO k=1,nlev
      DO iter=1,niter
        i = mpibufsize * (iter-1)
        DO j=1,mpibufsize
          i=i+1
          IF(i > nij1) EXIT
          bufs(j) = REAL(v3d(i,k,n),r_sngl)
        END DO
        CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
        CALL MPI_GATHER(bufs,ns,MPI_REAL,&
                      & bufr,nr,MPI_REAL,nrank,MPI_COMM_WORLD,ierr)
        IF(myrank == nrank) THEN
          i = mpibufsize * (iter-1)
          DO j=1,mpibufsize
            i=i+1
            IF(i > nij1max) EXIT
            tmp(i,:) = bufr(j,:)
          END DO
        END IF
      END DO
      IF(myrank == nrank) CALL buf_to_grd(tmp,v3dg(:,:,k,n))
    END DO
  END DO

  DO n=1,nv2d
    DO iter=1,niter
      i = mpibufsize * (iter-1)
      DO j=1,mpibufsize
        i=i+1
        IF(i > nij1) EXIT
        bufs(j) = REAL(v2d(i,n),r_sngl)
      END DO
      CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
      CALL MPI_GATHER(bufs,ns,MPI_REAL,&
                    & bufr,nr,MPI_REAL,nrank,MPI_COMM_WORLD,ierr)
      IF(myrank == nrank) THEN
        i = mpibufsize * (iter-1)
        DO j=1,mpibufsize
          i=i+1
          IF(i > nij1max) EXIT
          tmp(i,:) = bufr(j,:)
        END DO
      END IF
    END DO
    IF(myrank == nrank) CALL buf_to_grd(tmp,v2dg(:,:,n))
  END DO

  CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)

  RETURN
END SUBROUTINE gather_grd_mpi_safe

SUBROUTINE gather_grd_mpi_1d_dp(nrank,v1d,v1dg)
  INTEGER,INTENT(IN) :: nrank
  REAL(r_size),INTENT(IN)  :: v1d(nij1)
  REAL(r_size),INTENT(OUT) :: v1dg(nlon,nlat)
  REAL(r_size) :: tmp(nij1max,nprocs)
  REAL(r_size) :: bufs(mpibufsize)
  REAL(r_size) :: bufr(mpibufsize,nprocs)
  INTEGER :: i,j,ierr,ns,nr
  INTEGER :: iter,niter

  ns = mpibufsize
  nr = ns
  niter = CEILING(REAL(nij1max)/REAL(mpibufsize))

  DO iter=1,niter
    i = mpibufsize * (iter-1)
    DO j=1,mpibufsize
      i=i+1
      IF(i > nij1) EXIT
      bufs(j) = v1d(i)
    END DO
    CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
    CALL MPI_GATHER(bufs,ns,MPI_DOUBLE_PRECISION,&
    & bufr,nr,MPI_DOUBLE_PRECISION,nrank,MPI_COMM_WORLD,ierr)
    IF(myrank == nrank) THEN
      i = mpibufsize * (iter-1)
      DO j=1,mpibufsize
        i=i+1
        IF(i > nij1max) EXIT
        tmp(i,:) = bufr(j,:)
      END DO
    END IF
  END DO

  ! check
  IF(myrank == nrank) CALL buf_to_grd_dp(tmp,v1dg(:,:))
  CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)

  RETURN
END SUBROUTINE gather_grd_mpi_1d_dp

SUBROUTINE gather_grd_mpi_fast(nrank,v3d,v2d,v3dg,v2dg)
  INTEGER,INTENT(IN) :: nrank
  REAL(r_size),INTENT(IN) :: v3d(nij1,nlev,nv3d)
  REAL(r_size),INTENT(IN) :: v2d(nij1,nv2d)
  REAL(r_sngl),INTENT(OUT) :: v3dg(nlon,nlat,nlev,nv3d)
  REAL(r_sngl),INTENT(OUT) :: v2dg(nlon,nlat,nv2d)
  REAL(r_sngl) :: bufs(nij1max,nlevall)
  REAL(r_sngl) :: bufr(nij1max,nlevall,nprocs)
  INTEGER :: j,k,n,ierr,ns,nr

  ns = nij1max * nlevall
  nr = ns
  j=0
  DO n=1,nv3d
    DO k=1,nlev
      j = j+1
      bufs(1:nij1,j) = REAL(v3d(:,k,n),r_sngl)
    END DO
  END DO

  DO n=1,nv2d
    j = j+1
    bufs(1:nij1,j) = REAL(v2d(:,n),r_sngl)
  END DO

  CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
  CALL MPI_GATHER(bufs,ns,MPI_REAL,&
                & bufr,nr,MPI_REAL,nrank,MPI_COMM_WORLD,ierr)

  IF(myrank == nrank) THEN
    j=0
    DO n=1,nv3d
      DO k=1,nlev
        j = j+1
        CALL buf_to_grd(bufr(:,j,:),v3dg(:,:,k,n))
      END DO
    END DO

    DO n=1,nv2d
      j = j+1
      CALL buf_to_grd(bufr(:,j,:),v2dg(:,:,n))
    END DO
  END IF

  CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)

  RETURN
END SUBROUTINE gather_grd_mpi_fast
!-----------------------------------------------------------------------
! Read ensemble data and distribute to processes
!-----------------------------------------------------------------------
SUBROUTINE read_ens_mpi(file,member,v3d,v2d)
  CHARACTER(4),INTENT(IN) :: file
  INTEGER,INTENT(IN) :: member
  REAL(r_size),INTENT(OUT) :: v3d(nij1,nlev,member,nv3d)
  REAL(r_size),INTENT(OUT) :: v2d(nij1,member,nv2d)
  REAL(r_sngl) :: v3dg(nlon,nlat,nlev,nv3d)
  REAL(r_sngl) :: v2dg(nlon,nlat,nv2d)
  INTEGER :: l,n,ll,im
  !!CHARACTER(11) :: filename='file000.grd'
  CHARACTER(26) :: filename='gues/MMMMMM/yyyymmddhh.grd'

  ll = CEILING(REAL(member)/REAL(nprocs))
  DO l=1,ll
    im = myrank+1 + (l-1)*nprocs
    IF(im <= member) THEN
      !!WRITE(filename(1:7),'(A4,I3.3)') file,im
      WRITE(filename(6:11), '(I6.6)') im
      WRITE(filename(13:22), '(I10.10)') ymdh
      WRITE(6,'(A,I3.3,2A)') 'MYRANK ',myrank,' is reading a file ',filename
      CALL read_grd4(filename,v3dg,v2dg)
    END IF

    DO n=0,nprocs-1
      im = n+1 + (l-1)*nprocs
      IF(im <= member) THEN
        CALL scatter_grd_mpi(n,v3dg,v2dg,v3d(:,:,im,:),v2d(:,im,:))
      END IF
    END DO
  END DO

  RETURN
END SUBROUTINE read_ens_mpi
!-----------------------------------------------------------------------
! Write ensemble data after collecting data from processes
!-----------------------------------------------------------------------
SUBROUTINE write_ens_mpi(file,member,v3d,v2d)
  CHARACTER(4),INTENT(IN) :: file
  INTEGER,INTENT(IN) :: member
  REAL(r_size),INTENT(IN) :: v3d(nij1,nlev,member,nv3d)
  REAL(r_size),INTENT(IN) :: v2d(nij1,member,nv2d)
  REAL(r_sngl) :: v3dg(nlon,nlat,nlev,nv3d)
  REAL(r_sngl) :: v2dg(nlon,nlat,nv2d)
  INTEGER :: l,n,ll,im
  !CHARACTER(11) :: filename='file000.grd'
  CHARACTER(26) :: filename='anal/MMMMMM/yyyymmddhh.grd'

  ll = CEILING(REAL(member)/REAL(nprocs))
  DO l=1,ll
    DO n=0,nprocs-1
      im = n+1 + (l-1)*nprocs
      IF(im <= member) THEN
        CALL gather_grd_mpi(n,v3d(:,:,im,:),v2d(:,im,:),v3dg,v2dg)
      END IF
    END DO

    im = myrank+1 + (l-1)*nprocs
    IF(im <= member) THEN
      !WRITE(filename(1:7),'(A4,I3.3)') file,im
      WRITE(filename(6:11), '(I6.6)') im
      WRITE(filename(13:22), '(I10.10)') ymdh
      WRITE(6,'(A,I3.3,2A)') 'MYRANK ',myrank,' is writing a file ',filename
      CALL write_grd4(filename,v3dg,v2dg)
    END IF
  END DO

  RETURN
END SUBROUTINE write_ens_mpi
!-----------------------------------------------------------------------
! gridded data -> buffer
!-----------------------------------------------------------------------
SUBROUTINE grd_to_buf(grd,buf)
  REAL(r_sngl),INTENT(IN) :: grd(nlon,nlat)
  REAL(r_sngl),INTENT(OUT) :: buf(nij1max,nprocs)
  INTEGER :: i,j,m,ilon,ilat

  DO m=1,nprocs
    DO i=1,nij1node(m)
      j = m-1 + nprocs * (i-1)
      ilon = MOD(j,nlon) + 1
      ilat = (j-ilon+1) / nlon + 1
      buf(i,m) = grd(ilon,ilat)
    END DO
  END DO

  DO m=1,nprocs
    IF(nij1node(m) < nij1max) buf(nij1max,m) = undef
  END DO

  RETURN
END SUBROUTINE grd_to_buf
!-----------------------------------------------------------------------
! buffer -> gridded data
!-----------------------------------------------------------------------
SUBROUTINE buf_to_grd(buf,grd)
  REAL(r_sngl),INTENT(IN) :: buf(nij1max,nprocs)
  REAL(r_sngl),INTENT(OUT) :: grd(nlon,nlat)
  INTEGER :: i,j,m,ilon,ilat

  DO m=1,nprocs
    DO i=1,nij1node(m)
      j = m-1 + nprocs * (i-1)
      ilon = MOD(j,nlon) + 1
      ilat = (j-ilon+1) / nlon + 1
      grd(ilon,ilat) = buf(i,m)
    END DO
  END DO

  RETURN
END SUBROUTINE buf_to_grd
!-----------------------------------------------------------------------
SUBROUTINE buf_to_grd_dp(buf,grd)
  REAL(r_size),INTENT(IN) :: buf(nij1max,nprocs)
  REAL(r_size),INTENT(OUT) :: grd(nlon,nlat)
  INTEGER :: i,j,m,ilon,ilat

  DO m=1,nprocs
    DO i=1,nij1node(m)
      j = m-1 + nprocs * (i-1)
      ilon = MOD(j,nlon) + 1
      ilat = (j-ilon+1) / nlon + 1
      grd(ilon,ilat) = buf(i,m)
    END DO
  END DO

  RETURN
END SUBROUTINE buf_to_grd_dp
!-----------------------------------------------------------------------
! STORING DATA (ensemble mean and spread)
!-----------------------------------------------------------------------
SUBROUTINE write_ensmspr_mpi(file,member,v3d,v2d)
  CHARACTER(4),INTENT(IN) :: file
  INTEGER,INTENT(IN) :: member
  REAL(r_size),INTENT(IN) :: v3d(nij1,nlev,member,nv3d)
  REAL(r_size),INTENT(IN) :: v2d(nij1,member,nv2d)
  REAL(r_size) :: v3dm(nij1,nlev,nv3d)
  REAL(r_size) :: v2dm(nij1,nv2d)
  REAL(r_size) :: v3ds(nij1,nlev,nv3d)
  REAL(r_size) :: v2ds(nij1,nv2d)
  REAL(r_size) :: v3dk(nij1,nlev,nv3d)
  REAL(r_size) :: v2dk(nij1,nv2d)
  REAL(r_sngl) :: v3dg(nlon,nlat,nlev,nv3d)
  REAL(r_sngl) :: v2dg(nlon,nlat,nv2d)
  INTEGER :: i,k,m,n
  !CHARACTER(11) :: filename='file000.grd'
  CHARACTER(24) :: filename='file/mesp/yyyymmddhh.grd'

  CALL ensmean_grd(member,nij1,v3d,v2d,v3dm,v2dm)

  CALL gather_grd_mpi(0,v3dm,v2dm,v3dg,v2dg)
  IF(myrank == 0) THEN
    !WRITE(filename(1:7),'(A4,A3)') file,'_me'
    WRITE(filename(1:4),  '(A4)') file
    WRITE(filename(6:9),  '(A4)') 'mean'
    WRITE(filename(11:20),'(I10)') ymdh

    WRITE(6,'(A,I3.3,2A)') 'MYRANK ',myrank,' is writing a file ',filename
    CALL write_grd4(filename,v3dg,v2dg)
  END IF

  !==> KLD based on KK
  CALL enskldv_grd(member,nij1,v3d,v2d,v3dk,v2dk)
  CALL gather_grd_mpi(0,v3dk,v2dk,v3dg,v2dg)
  IF(myrank == 0) THEN
    !WRITE(filename(1:7),'(A4,A3)') file,'_me'
    WRITE(filename(1:4),  '(A4)') file
    WRITE(filename(6:9),  '(A4)') 'kldv'
    WRITE(filename(11:20),'(I10)') ymdh

    WRITE(6,'(A,I3.3,2A)') 'MYRANK ',myrank,' is writing a file ',filename
    CALL write_grd4(filename,v3dg,v2dg)
  END IF

  DO n=1,nv3d
!$OMP PARALLEL DO PRIVATE(i,k,m)
    DO k=1,nlev
      DO i=1,nij1
        v3ds(i,k,n) = (v3d(i,k,1,n)-v3dm(i,k,n))**2
        DO m=2,member
          v3ds(i,k,n) = v3ds(i,k,n) + (v3d(i,k,m,n)-v3dm(i,k,n))**2
        END DO
        v3ds(i,k,n) = SQRT(v3ds(i,k,n) / REAL(member-1,r_size))
      END DO
    END DO
!$OMP END PARALLEL DO
  END DO

  DO n=1,nv2d
!$OMP PARALLEL DO PRIVATE(i,k,m)
    DO i=1,nij1
      v2ds(i,n) = (v2d(i,1,n)-v2dm(i,n))**2
      DO m=2,member
        v2ds(i,n) = v2ds(i,n) + (v2d(i,m,n)-v2dm(i,n))**2
      END DO
      v2ds(i,n) = SQRT(v2ds(i,n) / REAL(member-1,r_size))
    END DO
!$OMP END PARALLEL DO
  END DO

  CALL gather_grd_mpi(0,v3ds,v2ds,v3dg,v2dg)
  IF(myrank == 0) THEN
    !WRITE(filename(1:7),'(A4,A3)') file,'_sp'
    WRITE(filename(1:4),  '(A4)') file
    WRITE(filename(6:9),  '(A4)') 'sprd'
    WRITE(filename(11:20),'(I10)') ymdh
    WRITE(6,'(A,I3.3,2A)') 'MYRANK ',myrank,' is writing a file ',filename
    CALL write_grd4(filename,v3dg,v2dg)
  END IF

  RETURN
END SUBROUTINE write_ensmspr_mpi
!-----------------------------------------------------------------------
! Make iseed for com_randns based on grid number ! from kk
!-----------------------------------------------------------------------
SUBROUTINE mk_iseed(ymdh, ij,ilev,n,iseed)
  INTEGER, INTENT(IN) :: ymdh, ij, ilev, n
  INTEGER, INTENT(OUT) :: iseed
  INTEGER :: index_grd
  INTEGER :: index_date
  INTEGER :: dd, hh
  integer :: ymdhtmp=1982123112

! hh = ymdhtmp-nint(ymdhtmp*0.01d0)*100
! dd = nint( (ymdhtmp-nint(ymdhtmp*0.0001d0)*10000-hh) * 0.01d0)
  hh = ymdh-nint(ymdh*0.01d0)*100
  dd = nint( (ymdh-nint(ymdh*0.0001d0)*10000-hh) * 0.01d0)
  index_date = (dd-1)*24 + hh
  index_grd = (ij-1)*nprocs + myrank+1
  iseed = index_date*1000 + index_grd*100 + n*10 + ilev

  RETURN
END SUBROUTINE mk_iseed
END MODULE common_mpi_speedy

max_stars_repo_name: skotsuki/speedy-lpf

max_issues_repo_name: skotsuki/speedy-lpf

max_forks_repo_name: skotsuki/speedy-lpf


text:
program euler_056
   use EulerCommon
   integer(int64), dimension(200) :: ab
   integer(int64) :: a,b, s, max_s ,max_a, max_b

   max_s = 0
   do a=90,99
      do b=90,99
         call powers(a,b,ab)
         s = sum(ab)
         if (s > max_s) then
            max_s = s
            max_a = a
            max_b = b
         end if
      end do
   end do
   print *,max_a, max_b, max_s

contains
   subroutine powers(a,b,ab)
      integer(int64), intent(in) :: a,b
      integer(int64), intent(out) :: ab(200)
      integer(int64) :: i, max_n

      ab(:) = 0
      ab(1) = a
      max_n = max(200_int64,b*int(log10(real(a)),int64)+3_int64)
      call checkTensInArray(max_n,ab)
      do i=2,b
         ab(:) = ab(:)*a
         call checkTensInArray(max_n,ab)
      end do !- i
   end subroutine powers

end program euler_056

max_stars_repo_name: kylekanos/ProjectEuler

max_issues_repo_name: kylekanos/ProjectEuler

max_forks_repo_name: kylekanos/ProjectEuler


text:
Salsa Adicción, is a social dancing club at UC Davis that focuses on teaching componets of dance that are crucial for getting out onto the dance floor!  Guest instructors often are a part of the lessons, teaching beginning and beginningintermediate dance steps.

The focus of this group is mainly on two styles of Latin Social Dancing.
  Bachata,  a style of dancing that comes from the Dominican Republic and is known for being romantic and intimate.
  Salsa, one of the more traditional and well know of the Latin dances taught in standard salsa On1.

This Club is Currently Active: Winter Quarter 2013, Meetings at 78:30pm Thursdays, 1001 Giedt Hall.

UPCOMING EVENTS

Thursday, February 7th, 2013: 78:30pm, Salsa Adiccion Meeting in 1001 Giedt Hall

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!   EVB-QMDFF - RPMD molecular dynamics and rate constant calculations on
!               black-box generated potential energy surfaces
!
!   Copyright (c) 2021 by Julien Steffen (steffen@pctc.uni-kiel.de)
!                         Stefan Grimme (grimme@thch.uni-bonn.de) (QMDFF code)
!
!   Permission is hereby granted, free of charge, to any person obtaining a
!   copy of this software and associated documentation files (the "Software"),
!   to deal in the Software without restriction, including without limitation
!   the rights to use, copy, modify, merge, publish, distribute, sublicense,
!   and/or sell copies of the Software, and to permit persons to whom the
!   Software is furnished to do so, subject to the following conditions:
!
!   The above copyright notice and this permission notice shall be included in
!   all copies or substantial portions of the Software.
!
!   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
!   IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
!   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
!   THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
!   LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
!   FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
!   DEALINGS IN THE SOFTWARE.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!
!     subroutine ff_set: define a QMDFF amd initialize the forceconstants
!
subroutine ff_set(n,xyz,xyzvdw,at,nb,vdw,q,chb,cxb,pair, &
   &           r0ab,zab,r0094,sr42,cc6ab,echo,wbo,hyb, &
   &           cring,ringsize,mlist,mtors)
use qmdff
implicit none
integer::at(n),n,nb(20,n),pair(n*(n+1)/2),hyb(n)
integer::cring(8,n),ringsize(n),mlist(n),mtors
real(kind=8)::xyz(3,n),q(n),chb(94),cxb(94),wbo(n,n)
real(kind=8)::xyzvdw(3,n)
logical::vdw,ok,echo,pr,ring
real(kind=8)::r0ab(94,94),zab(94,94),r0094(94,94),sr42(94,94),cc6ab(n,n)

integer::iat,i,j,k,l,neig,ii,jj,kk,ll,m,iperm,iipp(3,3),ni,nj
integer::mol1,mm,iz1,iz2,nh,oo,nn,tag2(ndim),jk,jl,nxb
integer::lina,ij,tag(n*(n+1)/2),ii1,ii2,i1,i2,i3,pqn,irow
integer::iwt(2,10,n*(n+1)/2),iiwt(2,10),mem1(n),mem2(n),jrow
real(kind=8)::dum,dum1,dum2,dum3,rco,ri,rj,barr,xn,tpot,error,bav,aaa
real(kind=8)::s(ntterm),p(ntterm),o(ntterm),pp(ntterm,n*(n+1)/2),rav,bbb
real(kind=8)::angle,aa1,aa2,phi,rmax1,rmax2,rmin2,c,bb,r,fitmad,valijk,bthr
real(kind=8)::athr2,sthr1,sthr2,valijkl,wt(10,n*(n+1)/2),wwt(10),r2,omega
logical::checkpair,checktrip,checkfour,l1,l2,l3,cdbi,linear
logical::checkinvstr,doeht,samering,metaltorsion,cdbj,skiplin13

!
!     symmetrize to linear. if angle is within this value linear
!     this was 1.0 deg in the published version but has to be
!     increased for polyethiine units, SG Wed Nov 26 13:49:50 CET 2014
!
bthr=5.0d0
!
!     inlcude linear 1,3-stretch terms (if .F.)
!
skiplin13=.false.

write(10,*)
write(10,*)
write(10,*)'===================='
write(10,*)'FF setup: define all'
write(10,*)'QMDFF energy terms'
write(10,*)'===================='
write(10,*)
write(10,*) 'General parameters:'
write(10,*)
write(10,*)' parameter bthr: below which degree an angle is taken as liner?'
write(10,*) bthr
write(10,*)
write(10,*)' parameter skiplin13: if set F, all linear 1,3-strech terms are included!'
write(10,*) skiplin13
write(10,*)
iipp(1,1)=1
iipp(1,2)=2
iipp(1,3)=3
iipp(2,1)=2
iipp(2,2)=1
iipp(2,3)=3
iipp(3,1)=3
iipp(3,2)=2
iipp(3,3)=1
vdw=.false.
!
!      additional torsions around metal centers?
!
metaltorsion=.false.
if (mtors.ne.0) metaltorsion=.true.

write(10,*)' Are additional torsions around metal centers included? (True/False:)'
write(10,*) metaltorsion
write(10,*)
nangl=0
ntors=0
nhb  =0
!
!     do we have a vdW complex?
!
call splitmol(n,at,xyz,wbo,mlist,mol1,0,0)
if(mol1.ne.n) then
   write(10,*) ' found vdW complex: more than one molecule!'
   write(10,*)
!
!     make a geometry with separated monomers!
!     shift all coordinates around 1000!!
!
   vdw=.true.
   do i=1,n
      if(mlist(i).eq.1)then
         xyzvdw(1:3,i)=xyz(1:3,i)
      else
         xyzvdw(1:3,i)=xyz(1:3,i)+10000.
      endif
   enddo
else
   mlist=0
endif
!
!     stretch. set in ff_bond.f!!!!!
!
do k=1,nbond
   i=bond(1,k)
   j=bond(2,k)
   r=sqrt((xyz(1,i)-xyz(1,j))**2 &
     &   +(xyz(2,i)-xyz(2,j))**2 &
     &   +(xyz(3,i)-xyz(3,j))**2)
   vbond(1,k)=r
   vbond(2,k)=1.0/r/vbond(3,k)**2
end do
nbond12=nbond
!
!     add Hbonds
!
nhb=0
k  =0
do i=1,n
   do j=1,i
      k=k+1
      if(i.eq.j)cycle
      if(pair(k).lt.3.and.pair(k).ne.0) cycle
      iz1=at(i)
      iz2=at(j)
      c=chb(iz1)*chb(iz2)
      if (c.gt.1.d-6) then
         do nh=1,n
            if(at(nh).ne.1) cycle
            ri=sqrt((xyz(1,i)-xyz(1,nh))**2 &
              &    +(xyz(2,i)-xyz(2,nh))**2 &
              &    +(xyz(3,i)-xyz(3,nh))**2)
            rj=sqrt((xyz(1,j)-xyz(1,nh))**2 &
              &    +(xyz(2,j)-xyz(2,nh))**2 &
              &    +(xyz(3,j)-xyz(3,nh))**2)
            dum1=1.3*(rad(at(i))+rad(1))/0.52917726
            dum2=1.3*(rad(at(j))+rad(1))/0.52917726
            if (ri.lt.dum1.or.rj.lt.dum2)then
               nhb=nhb+1
               hb(1,nhb)=i
               hb(2,nhb)=j
               hb(3,nhb)=nh
               call hbpara(10.0d0,5.0d0,q(i),c)
               vhb(1,nhb)=c*chb(iz1)
               call hbpara(10.0d0,5.0d0,q(j),c)
               vhb(2,nhb)=c*chb(iz2)
            end if
         end do
         100 continue
      end if
   end do
end do
write(10,'(I3,A)') nhb ,' HB terms (hydrogen bonds)'
!
!     add halogen bonds
!
nxb=0
k  =0
do i=1,n
   do j=1,i
      k=k+1
      if (i.eq.j)cycle
      if (pair(k).lt.3.and.pair(k).ne.0) cycle
      iz1=at(i)
      iz2=at(j)
      cdbi=iz1.eq.7.or.iz1.eq.8
      cdbj=iz2.eq.7.or.iz2.eq.8
      if (iz1.eq.1.or.iz2.eq.1) cycle
      do nh=1,n
         if (nh.eq.i) cycle
         if (nh.eq.j) cycle
         if (at(nh).eq.17.or.at(nh).eq.35.or. &
           &  at(nh).eq.53.or.at(nh).eq.85) then
            ri=sqrt((xyz(1,i)-xyz(1,nh))**2 &
                &  +(xyz(2,i)-xyz(2,nh))**2 &
                &  +(xyz(3,i)-xyz(3,nh))**2)
            rj=sqrt((xyz(1,j)-xyz(1,nh))**2 &
                &  +(xyz(2,j)-xyz(2,nh))**2 &
                &  +(xyz(3,j)-xyz(3,nh))**2)
            dum1=1.2*(rad(at(i))+rad(at(nh)))/0.52917726
            dum2=1.2*(rad(at(j))+rad(at(nh)))/0.52917726
!
!     Y(I=1)-X..D(J=2)
!     the charges on X and Y determine the strength
!
            if (ri.lt.dum1.and.cdbj) then
               nxb=nxb+1
               hb(1,nxb+nhb)=i
               hb(2,nxb+nhb)=j
               hb(3,nxb+nhb)=nh
               call hbpara(-6.5d0,1.0d0,q(nh),c)
               vhb(1,nxb+nhb)=c*cxb(at(nh))
            end if
!
!     Y(j=1)-X..D(I=2)
!     another configuration
!
            if (rj.lt.dum2.and.cdbi) then
               nxb=nxb+1
               hb(1,nxb+nhb)=j
               hb(2,nxb+nhb)=i
               hb(3,nxb+nhb)=nh
               call hbpara(-6.5d0,1.0d0,q(nh),c)
               vhb(1,nxb+nhb)=c*cxb(at(nh))
            end if
         end if
      end do
   end do
end do
write(10,'(I3,A)') nxb ,' XB terms (halogen bonds)'
write(10,*)
nhb=nhb+nxb

!
!     bend (bond angle between A-B-C)
!
do i=1,n
   neig=nb(20,i)
   ii=i
!
!     no bending for highly coordinated (metal) atoms
!     without this skip, cp rings do not rotate
!     highest possible coordination is 6!
!
   if(neig.gt.6) cycle
   do j=1,neig
      do k=1,j-1
         jj=nb(j,i)
         kk=nb(k,i)
         if (j.eq.k) cycle
         call bangle (XYZ,JJ,II,KK,ANGLE)
         r=angle*180.0d0/pi
!
!     small angles are not significant and redundant
!
         athr2=80.0d0
         l1=metal(at(ii)).eq.1.or. &
            &  metal(at(jj)).eq.1.or. &
            &  metal(at(kk)).eq.1
         if (l1) athr2=60.0d0
         if (r.lt.athr2) cycle
!
!     avoid close to linear situations in highly
!     coordinated systems
!
         l1= nb(20,jj).gt.3
         l2= nb(20,kk).gt.3
         if (r.gt.160.0d0.and.l1.and.l2) cycle
         do m=1,nangl
            if(checktrip(ii,jj,kk,angl(1,m),angl(2,m),angl(3,m))) &
                 & goto 199
         end do
         nangl=nangl+1
         if (nangl.gt.ndim) stop 'too many angles'
         angl(1,nangl)=ii
         angl(2,nangl)=jj
         angl(3,nangl)=kk
         vangl(1,nangl)=r*pi/180.
!
!     symmetrize to linear. without this, close to
!     linear situations become unstable
!
         if(180.0d0-r.lt.bthr) vangl(1,nangl)=pi
         vangl(2,nangl)=0.02
         199  continue
      end do
   end do
end do

!
!     compute for each bond the
!     TB torsion potential
!
write(10,*) "Check for rings in reference system..."
do i=1,n
   call getring(n,nb,i,cring(1,i),ringsize(i))
   if (ringsize(i).gt.0) then
      write(10,'('' atom'',i4,'' is in ring '',8i4)') &
         &   i,cring(1:ringsize(i),i)
   end if
end do
write(10,*)
write(10,*) 'computing torsion potentials by Extended Hückel theory'
write(10,*) '   for rotatable bonds ...'
!
!     define the torsion potential.
!
do i=1,ntterm
   o(i)=i
   s(i)=pi
end do

!     else this part is still parity violating, see ffpot_both*
!        do i=1,ntterm/2
!        o(i)=i
!        s(i)=pi
!        enddo
!        do i=ntterm/2+1,ntterm
!        o(i)=i-ntterm/2
!        s(i)=pi/2
!        enddo
!    endif

k=0
do i=1,n
   do j=1,i
      k=k+1
      pp(1:ntterm,k)=0
      ok=.false.
!
!     check cases for which we don't need to compute the potential
!
      if (i.eq.j.or.pair(k).ne.1) cycle
      if (nb(20,j).eq.1.or.nb(20,j).gt.4) cycle
      if (nb(20,i).eq.1.or.nb(20,i).gt.4) cycle
      r=sqrt((xyz(1,i)-xyz(1,j))**2 &
         &  +(xyz(2,i)-xyz(2,j))**2 &
         &  +(xyz(3,i)-xyz(3,j))**2)
!
!     check if TB/EHT should be done
!
      l1=doeht(n,at,i,j,r,hyb,cring,ringsize,wbo)
      if(l1) then
         p(1:ntterm)=0.001
!
!     do the TB/EHT pot curve and the fit
!
         call compmodel(echo,.false.,i,j,n,xyz,at,q,nb, &
           &  r0ab,zab,r0094,sr42,c6ab,wbo,ntterm,ok,barr,fitmad, &
           &  s,o,p,iiwt,wwt,cring,ringsize)
!
!     this avoids the TB/EHT potential for conjugated systems like C60
!
         if (barr.gt.25) ok=.false.
!
!     this avoids artificial model systems
!
         if (barr.lt.0.05) ok=.false.
         write(10,'('' pair '',2i4,2x &
              &  ''barrier (kcal) :'',f6.2,2x, &
              &  '' %MAD(fit):'',f5.1, &
              &  '' fit pot used?'',l)') &
              &  i,j,barr,fitmad*100.,ok
      end if
      if (ok) then
!
!     normalize by the # of torsions for ij
!
         mm=0
         do m=1,ntterm
            pp(m,k)=p(m)/dble((nb(20,i)-1)*(nb(20,j)-1))
         end do
         iwt(1:2,1:10,k)=iiwt(1:2,1:10)
         wt(1:10,k)    =     wwt(1:10)
      end if
   end do
end do
!
!     torsional terms
!
tag=0
do i=1,n
   ii=i
!
!     no torsions if not coordinated or highly coordinated
!
   if (nb(20,ii).eq.1.or.nb(20,ii).gt.5) cycle
   do j=1,nb(20,ii)
      jj=nb(j,ii)
      ij=lina(jj,ii)
!
!     already done this ij?
!
      if (tag(ij).eq.1) cycle
      if (nb(20,jj).eq.1.or.nb(20,jj).gt.5) cycle
      do k=1,nb(20,jj)
         kk=nb(k,jj)
         if (kk.eq.ii) cycle
         do l=1,nb(20,ii)
            ll=nb(l,ii)
            if (ll.eq.jj.or.ll.eq.kk) cycle

            r=valijkl(n,xyz,ll,ii,jj,kk)*180./pi
!
!     check for linear arrangements
!
            call bangle(XYZ,kk,ii,jj,aa1)
            call bangle(XYZ,ll,ii,jj,aa2)
            aa1=aa1*180./pi
            aa2=aa2*180./pi
            if (aa1.lt.athr.or.180-aa1.lt.athr) cycle
            if (aa2.lt.athr.or.180-aa2.lt.athr) cycle
!
!     exit if rotational number is 0 (eg 5-membered rings)
!
            if (pp(1,ij).lt.1.d-6) then
               call getrot(n,at,ii,jj,hyb,cring,ringsize,xn)
               if(xn.lt.1.d-6) cycle
            end if

            ntors=ntors+1
            if(ntors.gt.ndim) stop 'too many torsions'
            tors(1,ntors)=ll
            tors(2,ntors)=ii
            tors(3,ntors)=jj
            tors(4,ntors)=kk
            vtors(1,ntors)=r*pi/180.0d0
!
!     fitted pot case (done with EHT above)
!
            if (pp(1,ij).ne.0) then
               tors(5,ntors)=ntterm
!
!     no fit
!
               tors(6,ntors)=0
               vtors(2,ntors)=1
!
!     use computed weights
!
               do m=1,iwt(1,10,ij)
                  if (ll.eq.iwt(1,m,ij).and.kk.eq.iwt(2,m,ij)) then
                     vtors(2,ntors)=wt(m,ij)
                  end if
                  if (ll.eq.iwt(2,m,ij).and.kk.eq.iwt(1,m,ij)) then
                     vtors(2,ntors)=wt(m,ij)
                  end if
               end do
               mm=3
               do m=1,ntterm
                  vtors(mm  ,ntors)=o (m)
                  vtors(mm+1,ntors)=s (m)
                  vtors(mm+2,ntors)=pp(m,ij)
                  mm=mm+3
               end do
            else
!
!     fitted to Hessian with one-term (sine: shift=pi) potential
!
               call getrot(n,at,ii,jj,hyb,cring,ringsize,xn)
               tors(5,ntors)=1
               tors(6,ntors)=1
               vtors(2,ntors)=0.001
               vtors(3,ntors)=xn
               vtors(4,ntors)=pi
               vtors(5,ntors)=1.0
            end if
!
!     mark torsion as done for this ij pair (avoids doubles)
!
            tag(ij)=1
!
!     next lijk
!
         end do
      end do
   end do
end do

!
!     special case of metal atoms
!

if(metaltorsion)then
   write(10,*)'torsion terms at metal centers with rot=',mtors,&
       &  'switched on'
   do i=1,n
      m=0
      ni=nb(20,i)
      if (ni.gt.4) then
         do i1=1,ni
            do i2=1,i1-1
               kk=nb(i1,i)
               ll=nb(i2,i)
               if (nb(20,kk).eq.1) cycle
               if (nb(20,ll).eq.1) cycle
               call bangle(xyz,kk,i,ll,angle)
               r=angle*180/pi
               if (r.gt.100.and.r.lt.160.and.m.lt.6) then
                  jk=nb(1,kk)
                  jl=nb(1,ll)
                  if (.not.checkfour(kk,i,ll,jl)) cycle
                  if (.not.checkfour(ll,i,kk,jk)) cycle
                  m=m+1
                  r=valijkl(n,xyz,kk,i,ll,jl)
                  ntors=ntors+1
                  tors(1,ntors)=kk
                  tors(2,ntors)=i
                  tors(3,ntors)=ll
                  tors(4,ntors)=jl
                  tors(5,ntors)=1
                  tors(6,ntors)=1
                  vtors(1,ntors)=r
                  vtors(2,ntors)=0.001
!   Grubbs        vtors(3,ntors)=4
                  vtors(3,ntors)=mtors
                  vtors(4,ntors)=pi
                  vtors(5,ntors)=1
                  r=valijkl(n,xyz,ll,i,kk,jk)
                  ntors=ntors+1
                  tors(1,ntors)=ll
                  tors(2,ntors)=i
                  tors(3,ntors)=kk
                  tors(4,ntors)=jk
                  tors(5,ntors)=1
                  tors(6,ntors)=1
                  vtors(1,ntors)=r
                  vtors(2,ntors)=0.001
                  vtors(3,ntors)=mtors
                  vtors(4,ntors)=pi
                  vtors(5,ntors)=1
               end if
            end do
         end do
      end if
   end do
end if

!
!     inversion type torsion (central atom has 3 bonds)
!
do i=1,n
   ii=i
   l1=nb(20,ii).eq.3
   l2=nb(20,ii).eq.4.and.nb(18,ii).eq.1
   l3=nb(20,ii).eq.4.and.metal(at(ii)).eq.1
   if (l1.or.l2.or.l3) then
      jj=nb(1,ii)
      kk=nb(2,ii)
      ll=nb(3,ii)
      r=omega (n,xyz,ii,jj,kk,ll)
      ntors=ntors+1
      if (ntors.gt.ndim) stop 'too many torsions'
      tors(1,ntors)=ii
      tors(2,ntors)=jj
      tors(3,ntors)=kk
      tors(4,ntors)=ll
      tors(5,ntors)=1
      tors(6,ntors)=2
      vtors(1,ntors)=r
      vtors(2,ntors)=0.1
!      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!      the 'if' takes the chiraliy breaking term for planar cases
!      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      if (abs(r*180.0d0/pi).le.athr) then
         vtors(3,ntors)=1
!
!     symmetrize it to planarity
!
         vtors(1,ntors)=0
      else
!
!     this is the double minimum around 0
!
         vtors(3,ntors)=0
      end if
      vtors(4,ntors)=pi
      vtors(5,ntors)=1
   end if
end do

!cccccccccccccccccccccccccccc
!     torsion terms for linear R-A-X-B-R
!     where the torsion is then taken for
!     X-A..B-Y (eg allene)
!cccccccccccccccccccccccccccc
tag2=0
do i=1,nangl
   if (abs(vangl(1,i)*180./pi-180.).le.athr.and.tag2(i).eq.0) then
       ii=angl(1,i)
       jj=angl(2,i)
       kk=angl(3,i)
       do j=1,nb(20,jj)
          if (nb(j,jj).eq.ii) cycle
          do k=1,nb(20,kk)
             if (nb(k,kk).eq.ii) cycle

             r=valijkl(n,XYZ,nb(j,jj),jj,kk,nb(k,kk))
!
!     check for linear arrangements
!
             call bangle(xyz,nb(j,jj),kk,jj,aa1)
             call bangle(xyz,nb(k,kk),kk,jj,aa2)
             aa1=aa1*180./pi
             aa2=aa2*180./pi
             if (aa1.lt.athr.or.180-aa1.lt.athr) cycle
             if (aa2.lt.athr.or.180-aa2.lt.athr) cycle

             ntors=ntors+1
             if (ntors.gt.ndim) stop 'too many torsions'
             tors(1,ntors)=nb(j,jj)
             tors(2,ntors)=jj
             tors(3,ntors)=kk
             tors(4,ntors)=nb(k,kk)
             tors(5,ntors)=1
             tors(6,ntors)=1
             vtors(1,ntors)=r
             vtors(2,ntors)=0.001
             vtors(3,ntors)=2
             vtors(4,ntors)=pi
             vtors(5,ntors)=1
          end do
       end do
   end if
end do

!cccccccccccccccccccccccccccc
!     torsion terms for linear A-X-B-R
!     where the torsion is then taken for
!     A...R-Y
!cccccccccccccccccccccccccccc
do i=1,nangl
   if (abs(vangl(1,i)*180./pi-180.).le.athr.and.tag2(i).eq.0) then
      ii=angl(1,i)
      jj=angl(2,i)
      kk=angl(3,i)
      ll=0
      if (nb(20,jj).eq.1) then
         ll=kk
         mm=jj
      end if
      if (nb(20,kk).eq.1) then
         ll=jj
         mm=kk
      end if
      if(ll.eq.0) goto 200
      do j=1,nb(20,ll)
         nn=nb(j,ll)
         if (nn.eq.ll) cycle
         if (nn.eq.ii) cycle
         do k=1,nb(20,nn)
            oo=nb(k,nn)
            if (oo.eq.ll) cycle
            if (oo.eq.nn) cycle
            r=valijkl(n,xyz,mm,ii,nn,oo)
!
!     check for linear arrangements
!
            call bangle(xyz,mm,ii,nn,aa1)
            call bangle(xyz,oo,nn,ii,aa2)
            aa1=aa1*180./pi
            aa2=aa2*180./pi
            if (aa1.lt.athr.or.180-aa1.lt.athr) cycle
            if (aa2.lt.athr.or.180-aa2.lt.athr) cycle

            tag2(i)=1
            ntors=ntors+1
            if (ntors.gt.ndim) stop 'too many torsions'
            tors(1,ntors)=mm
            tors(2,ntors)=ii
            tors(3,ntors)=nn
            tors(4,ntors)=oo
            tors(5,ntors)=1
            tors(6,ntors)=1
            vtors(1,ntors)=r
            vtors(2,ntors)=0.001
            vtors(3,ntors)=4
            vtors(4,ntors)=pi
            vtors(5,ntors)=1
         end do
      end do
   end if
   200  continue
end do

!ccccccccccccccccccccccccccccccccccccc
!     stretch for bend: only define parameters!
!ccccccccccccccccccccccccccccccccccccc
88  continue
rmax1=-1

! for fit details see ff_bond.f
!     open(unit=66,file='~/.param')
!     read(66,*) aaa,bbb
!     read(66,*) ade13(1),ade13(5:9)
!     close(66)

aaa=2.81
bbb=0.53

l=nbond
do i=1,nangl
   ii=angl(1,i)
   jj=angl(2,i)
   kk=angl(3,i)
   linear=abs(vangl(1,i)*180./pi-180.).le.athr
   if (linear) then
!
!     its not entirely clear if 1,3-str should be kept for
!     linear cases. In the published version they were
!     discarded and if included e.g. for c2h2 they are
!     automatically removed (FC<0). for other cases this is,
!     however, different
!
      if (skiplin13) cycle
   end if
   l1=metal(at(ii)).eq.1.or. &
     &  metal(at(jj)).eq.1.or. &
     &  metal(at(kk)).eq.1
   r=sqrt((xyz(1,kk)-xyz(1,jj))**2 &
     &  +(xyz(2,kk)-xyz(2,jj))**2 &
     &  +(xyz(3,kk)-xyz(3,jj))**2)*0.52917726
   if (mlist(jj).eq.mlist(ii).and.r.gt.rmax1) rmax1=r
   rco=rad(at(kk))+rad(at(jj))
!
!     already there
!
   do m=1,l
      if(jj.eq.bond(1,m).and.kk.eq.bond(2,m)) goto 299
      if(jj.eq.bond(2,m).and.kk.eq.bond(1,m)) goto 299
   end do
   l=l+1
   bond(1,l)=jj
   bond(2,l)=kk
!
!     both neighbors are in a ring
!
   l1=ringsize(kk).gt.0.and.ringsize(jj).gt.0
!
!     De not sensitive to this if be values (ff_bond) are adjusted
!
   if (l1) then
      vbond(3,l)=aaa+bbb*ade13(at(jj))*ade13(at(kk))+0.7
   else
      vbond(3,l)=aaa+bbb*ade13(at(jj))*ade13(at(kk))
   end if
!
!     adjusted to c4h2
!
   if (linear) vbond(3,l)=vbond(3,l)+2.5
   vbond(1,l)=r/0.52917726
   vbond(2,l)=1.0/vbond(1,l)/vbond(3,l)**2
   299  continue
end do
nbond=l

write(10,*)
write(10,'('' largest distance considered as covalent stretch (r(bohr), r(Ang)):'')')
write(10,'(2F8.2)')  rmax1/0.52917726,rmax1
write(10,*)

!ccccccccccccccccccccccccccccccccccccc
!     symmetrize (find out which bond is rotable)
!ccccccccccccccccccccccccccccccccccccc

!
!     larger values (higher chance for symmetrization) is better
!
sthr1=1.2
sthr2=0.2
!     first get atoms which become equivalent upon single-bond rotation
iwt=0
k=0
do i=1,n
   iz1=at(i)
   do j=1,i
      iz2=at(j)
      k=k+1
!
!     exit cases
!
      if (i.eq.j.or.pair(k).ne.1)         cycle
      if (wbo(j,i).gt.sthr1     )         cycle
      if (nb(20,j).eq.1.or.nb(20,j).gt.4) cycle
      if (nb(20,i).eq.1.or.nb(20,i).gt.4) cycle
      if (nb(18,i).eq.1.or.nb(18,j).eq.1) cycle
      if (metal(iz1).eq.1.or.metal(iz2).eq.1) cycle
      if (samering(n,i,j,cring,ringsize)) cycle
!     on i atom
      mem1=0
      do ni=1,nb(20,i)
         ii1=nb(ni,i)
         do nj=1,nb(20,i)
            ii2=nb(nj,i)
            if (ii1.eq.ii2) cycle
            if (at(ii1).ne.at(ii2)) cycle
            if (abs(wbo(i,ii1)-wbo(i,ii2)).gt.sthr2) cycle
            if (ii1.eq.i.or.ii2.eq.i) cycle
            if (ii1.eq.j.or.ii2.eq.j) cycle
            mem1(ii1)=mem1(ii1)+1
            mem1(ii2)=mem1(ii2)+1
         end do
      end do
      mm=0
      do kk=1,n
         if (mem1(kk).gt.0) then
            mm=mm+1
            iwt(1,mm,k)=kk
            iwt(1,10,k)=mm
         end if
      end do
!     on j atom
      mem2=0
      do ni=1,nb(20,j)
         ii1=nb(ni,j)
         do nj=1,nb(20,j)
            ii2=nb(nj,j)
            if (ii1.eq.ii2) cycle
            if (at(ii1).ne.at(ii2)) cycle
            if (abs(wbo(j,ii1)-wbo(j,ii2)).gt.sthr2) cycle
            if (ii1.eq.j.or.ii2.eq.j) cycle
            if (ii1.eq.i.or.ii2.eq.i) cycle
            mem2(ii1)=mem2(ii1)+1
            mem2(ii2)=mem2(ii2)+1
         end do
      end do
      mm=0
      do kk=1,n
         if (mem2(kk).gt.0) then
            mm=mm+1
            iwt(2,mm,k)=kk
            iwt(2,10,k)=mm
         end if
      end do
   end do
end do

write(10,*) 'local symmetry analysis (rotable bonds) ...'
k=0
do i=1,n
   do j=1,i
      k=k+1
      kk=iwt(1,10,k)
      if (kk.gt.0)then
         write(10,'('' rotation on atom'',i4,'' symmetrizes:'',9i4)') &
           &  i,iwt(1,1:kk,k)
!
!     bonds and 1,3 stretches
!
         rav=0
         bav=0
         do mm=1,kk
            jj=iwt(1,mm,k)
            r=sqrt((xyz(1,i)-xyz(1,jj))**2 &
              &   +(xyz(2,i)-xyz(2,jj))**2 &
              &   +(xyz(3,i)-xyz(3,jj))**2)
            rav=rav+r
            r=sqrt((xyz(1,j)-xyz(1,jj))**2 &
              &   +(xyz(2,j)-xyz(2,jj))**2 &
              &   +(xyz(3,j)-xyz(3,jj))**2)
            bav=bav+r
         end do
         rav=rav/(mm-1)
         bav=bav/(mm-1)
         do mm=1,kk
            jj=iwt(1,mm,k)
            do l=1,nbond
               i1=bond(1,l)
               i2=bond(2,l)
               if (i1.eq.i.and.i2.eq.jj) vbond(1,l)=rav
               if (i2.eq.i.and.i1.eq.jj) vbond(1,l)=rav
               if (i1.eq.j.and.i2.eq.jj) vbond(1,l)=bav
               if (i2.eq.j.and.i1.eq.jj) vbond(1,l)=bav
            end do
         end do
!
!      angles
!
         rav=0
         do mm=1,kk
            jj=iwt(1,mm,k)
            call bangle(xyz,jj,i,j,angle)
            rav=rav+angle
         end do
         rav=rav/(mm-1)
         do mm=1,kk
            jj=iwt(1,mm,k)
            do l=1,nangl
               i1=angl(1,l)
               i2=angl(2,l)
               i3=angl(3,l)
               if (i1.eq.i.and.i2.eq.j.and.i3.eq.jj)vangl(1,l)=rav
               if (i1.eq.i.and.i3.eq.j.and.i2.eq.jj)vangl(1,l)=rav
            end do
         end do
      end if
      kk=iwt(2,10,k)
      if (kk.gt.0)then
          write(10,'('' rotation on atom'',i4,'' symmetrizes:'',9i4)') &
            &  j,iwt(2,1:kk,k)
         rav=0
         bav=0
         do mm=1,kk
            jj=iwt(2,mm,k)
            r=sqrt((xyz(1,j)-xyz(1,jj))**2 &
              &   +(xyz(2,j)-xyz(2,jj))**2 &
              &   +(xyz(3,j)-xyz(3,jj))**2)
            rav=rav+r
            r=sqrt((xyz(1,i)-xyz(1,jj))**2 &
              &   +(xyz(2,i)-xyz(2,jj))**2 &
              &   +(xyz(3,i)-xyz(3,jj))**2)
            bav=bav+r
         end do
         rav=rav/(mm-1)
         bav=bav/(mm-1)
         do mm=1,kk
            jj=iwt(2,mm,k)
            do l=1,nbond
               i1=bond(1,l)
               i2=bond(2,l)
               if (i1.eq.j.and.i2.eq.jj) vbond(1,l)=rav
               if (i2.eq.j.and.i1.eq.jj) vbond(1,l)=rav
               if (i1.eq.i.and.i2.eq.jj) vbond(1,l)=bav
               if (i2.eq.i.and.i1.eq.jj) vbond(1,l)=bav
            end do
         end do
!
!      angles
!
         rav=0
         do mm=1,kk
            jj=iwt(2,mm,k)
            call bangle(xyz,jj,j,i,angle)
            rav=rav+angle
         end do
         rav=rav/(mm-1)
         do mm=1,kk
            jj=iwt(2,mm,k)
            do l=1,nangl
               i1=angl(1,l)
               i2=angl(2,l)
               i3=angl(3,l)
               if (i1.eq.j.and.i2.eq.i.and.i3.eq.jj) vangl(1,l)=rav
               if (i1.eq.j.and.i3.eq.i.and.i2.eq.jj) vangl(1,l)=rav
            end do
         end do
      end if
   end do
end do

write(10,*)
write(10,*) 'force field setup finished!.'
write(10,*)

return
end subroutine ff_set

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: chrinide/EVB-QMDFF

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: chrinide/EVB-QMDFF

max_forks_repo_licenses: MIT Unlicense


text:
      PROGRAM CTIMAA
*
*  -- LAPACK timing routine (version 3.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     June 30, 1999
*
*  Purpose
*  =======
*
*  CTIMAA is the timing program for the COMPLEX LAPACK
*  routines.  This program collects performance data for the factor,
*  solve, and inverse routines used in solving systems of linear
*  equations, and also for the orthogonal factorization and reduction
*  routines used in solving least squares problems and matrix eigenvalue
*  problems.
*
*  The subprograms call a REAL function SECOND with no
*  arguments which is assumed to return the central-processor time in
*  seconds from some fixed starting time.
*
*  The program is driven by a short data file, which specifies values
*  for the matrix dimensions M, N and K, for the blocking parameters
*  NB and NX, and for the leading array dimension LDA.  A minimum time
*  for each subroutine is included for timing small problems or for
*  obtaining results on a machine with an inaccurate SECOND function.
*
*  The matrix dimensions M, N, and K correspond to the three dimensions
*  m, n, and k in the Level 3 BLAS.  When timing the LAPACK routines for
*  square matrices, M and N correspond to the matrix dimensions m and n,
*  and K is the number of right-hand sides (nrhs) for the solves.  When
*  timing the LAPACK routines for band matrices, M is the matrix order
*  m, N is the half-bandwidth (kl, ku, or kd in the LAPACK notation),
*  and K is again the number of right-hand sides.
*
*  The first 13 records of the data file are read using list-directed
*  input.  The first line of input is printed as the first line of
*  output and can be used to identify different sets of results.  To
*  assist with debugging an input file, the values are printed out as
*  they are read in.
*
*  The following records are read using the format (A).  For these
*  records, the first 6 characters are reserved for the path or
*  subroutine name.  If a path name is used, the characters after the
*  path name indicate the routines in the path to be timed, where
*  'T' or 't' means 'Time this routine'.  If the line is blank after the
*  path name, all routines in the path are timed.  If fewer characters
*  appear than routines in a path, the remaining characters are assumed
*  to be 'F'.  For example, the following 3 lines are equivalent ways of
*  requesting timing of CGETRF:
*  CGE    T F F
*  CGE    T
*  CGETRF
*
*  An annotated example of a data file can be obtained by deleting the
*  first 3 characters from the following 32 lines:
*  LAPACK timing, COMPLEX square matrices
*  5                                Number of values of M
*  100 200 300 400 500              Values of M (row dimension)
*  5                                Number of values of N
*  100 200 300 400 500              Values of N (column dimension)
*  2                                Number of values of K
*  100 400                          Values of K
*  5                                Number of values of NB
*  1 16  32  48  64                 Values of NB (blocksize)
*  0 48 128 128 128                 Values of NX (crossover point)
*  2                                Number of values of LDA
*  512 513                          Values of LDA (leading dimension)
*  0.0                              Minimum time in seconds
*  CGE    T T T
*  CPO    T T T
*  CPP    T T T
*  CHE    T T T
*  CHP    T T T
*  CSY    T T T
*  CSP    T T T
*  CTR    T T
*  CTP    T T
*  CQR    T T F
*  CLQ    T T F
*  CQL    T T F
*  CRQ    T T F
*  CQP    T
*  CHR    T T F F
*  CTD    T T F F
*  CBR    T F F
*  CLS    T T T T T T
*
*  The routines are timed for all combinations of applicable values of
*  M, N, K, NB, NX, and LDA, and for all combinations of options such as
*  UPLO and TRANS.  For Level 2 BLAS timings, values of NB are used for
*  INCX.  Certain subroutines, such as the QR factorization, treat the
*  values of M and N as ordered pairs and operate on M x N matrices.
*
*  Internal Parameters
*  ===================
*
*  NMAX    INTEGER
*          The maximum value of M or N for square matrices.
*
*  LDAMAX  INTEGER
*          The maximum value of LDA.
*
*  NMAXB   INTEGER
*          The maximum value of N for band matrices.
*
*  MAXVAL  INTEGER
*          The maximum number of values that can be read in for M, N,
*          K, NB, or NX.
*
*  MXNLDA  INTEGER
*          The maximum number of values that can be read in for LDA.
*
*  NIN     INTEGER
*          The unit number for input.  Currently set to 5 (std input).
*
*  NOUT    INTEGER
*          The unit number for output.  Currently set to 6 (std output).
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            NMAX, LDAMAX, NMAXB
      PARAMETER          ( NMAX = 512, LDAMAX = NMAX+20, NMAXB = 5000 )
      INTEGER            LA
      PARAMETER          ( LA = NMAX*LDAMAX )
      INTEGER            MAXVAL, MXNLDA
      PARAMETER          ( MAXVAL = 12, MXNLDA = 4 )
      INTEGER            MAXPRM
      PARAMETER          ( MAXPRM = MXNLDA*(MAXVAL+1) )
      INTEGER            MAXSZS
      PARAMETER          ( MAXSZS = MAXVAL*MAXVAL*MAXVAL )
      INTEGER            NIN, NOUT
      PARAMETER          ( NIN = 5, NOUT = 6 )
*     ..
*     .. Local Scalars ..
      LOGICAL            BLAS, LDAMOK, LDANOK, LDAOK, MOK, NOK, NXNBOK
      CHARACTER          C1
      CHARACTER*2        C2
      CHARACTER*3        C3
      CHARACTER*80       LINE
      INTEGER            I, I2, J2, L, LDR1, LDR2, LDR3, MAXK, MAXLDA,
     $                   MAXM, MAXN, MAXNB, MKMAX, NEED, NK, NLDA, NM,
     $                   NN, NNB
      REAL               S1, S2, TIMMIN
*     ..
*     .. Local Arrays ..
      INTEGER            IWORK( 2*NMAXB ), KVAL( MAXVAL ),
     $                   LDAVAL( MXNLDA ), MVAL( MAXVAL ),
     $                   NBVAL( MAXVAL ), NVAL( MAXVAL ),
     $                   NXVAL( MAXVAL )
      REAL               D( 2*NMAX ),
     $                   FLPTBL( 6*6*MAXSZS*MAXPRM*5 ),
     $                   OPCTBL( 6*6*MAXSZS*MAXPRM*5 ),
     $                   RESLTS( MAXVAL, MAXVAL, 2*MXNLDA, 4*MAXVAL ),
     $                   RWORK( 150*NMAX+2*MAXVAL ), S( NMAX*2 ),
     $                   TIMTBL( 6*6*MAXSZS*MAXPRM*5 )
      COMPLEX            A( LA, 3 ), B( LA, 3 ), E( 2*NMAX ),
     $                   WORK( NMAX, NMAX+MAXVAL+10 )
*     ..
*     .. External Functions ..
      LOGICAL            LSAME, LSAMEN
      REAL               SECOND
      EXTERNAL           LSAME, LSAMEN, SECOND
*     ..
*     .. External Subroutines ..
      EXTERNAL           CTIMB2, CTIMB3, CTIMBR, CTIMGB, CTIMGE, CTIMGT,
     $                   CTIMHE, CTIMHP, CTIMHR, CTIMLQ, CTIMLS, CTIMMM,
     $                   CTIMMV, CTIMPB, CTIMPO, CTIMPP, CTIMPT, CTIMQ3,
     $                   CTIMQL, CTIMQP, CTIMQR, CTIMRQ, CTIMSP, CTIMSY,
     $                   CTIMTB, CTIMTD, CTIMTP, CTIMTR
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MAX
*     ..
*     .. Arrays in Common ..
      INTEGER            IPARMS( 100 )
*     ..
*     .. Common blocks ..
      COMMON             / CLAENV / IPARMS
*     ..
*     .. Executable Statements ..
*
      S1 = SECOND( )
      LDR1 = MAXVAL
      LDR2 = MAXVAL
      LDR3 = 2*MXNLDA
      WRITE( NOUT, FMT = 9983 )
*
*     Read the first line.  The first four characters must be 'BLAS'
*     for the BLAS data file format to be used.  Otherwise, the LAPACK
*     data file format is assumed.
*
      READ( NIN, FMT = '( A80 )' )LINE
      BLAS = LSAMEN( 4, LINE, 'BLAS' )
*
*     Find the last non-blank and print the first line of input as the
*     first line of output.
*
      DO 10 L = 80, 1, -1
         IF( LINE( L: L ).NE.' ' )
     $      GO TO 20
   10 CONTINUE
      L = 1
   20 CONTINUE
      WRITE( NOUT, FMT = '( 1X, A, / )' )LINE( 1: L )
      WRITE( NOUT, FMT = 9992 )
*
*     Read in NM and the values for M.
*
      READ( NIN, FMT = * )NM
      IF( NM.GT.MAXVAL ) THEN
         WRITE( NOUT, FMT = 9999 )'M', 'NM', MAXVAL
         NM = MAXVAL
      END IF
      READ( NIN, FMT = * )( MVAL( I ), I = 1, NM )
      WRITE( NOUT, FMT = 9991 )'M:     ', ( MVAL( I ), I = 1, NM )
*
*     Check that  M <= NMAXB for all values of M.
*
      MOK = .TRUE.
      MAXM = 0
      DO 30 I = 1, NM
         MAXM = MAX( MVAL( I ), MAXM )
         IF( MVAL( I ).GT.NMAXB ) THEN
            WRITE( NOUT, FMT = 9997 )'M', MVAL( I ), NMAXB
            MOK = .FALSE.
         END IF
   30 CONTINUE
      IF( .NOT.MOK )
     $   WRITE( NOUT, FMT = * )
*
*     Read in NN and the values for N.
*
      READ( NIN, FMT = * )NN
      IF( NN.GT.MAXVAL ) THEN
         WRITE( NOUT, FMT = 9999 )'N', 'NN', MAXVAL
         NN = MAXVAL
      END IF
      READ( NIN, FMT = * )( NVAL( I ), I = 1, NN )
      WRITE( NOUT, FMT = 9991 )'N:     ', ( NVAL( I ), I = 1, NN )
*
*     Check that  N <= NMAXB for all values of N.
*
      NOK = .TRUE.
      MAXN = 0
      DO 40 I = 1, NN
         MAXN = MAX( NVAL( I ), MAXN )
         IF( NVAL( I ).GT.NMAXB ) THEN
            WRITE( NOUT, FMT = 9997 )'N', NVAL( I ), NMAXB
            NOK = .FALSE.
         END IF
   40 CONTINUE
      IF( .NOT.NOK )
     $   WRITE( NOUT, FMT = * )
*
*     Read in NK and the values for K.
*
      READ( NIN, FMT = * )NK
      IF( NK.GT.MAXVAL ) THEN
         WRITE( NOUT, FMT = 9999 )'K', 'NK', MAXVAL
         NK = MAXVAL
      END IF
      READ( NIN, FMT = * )( KVAL( I ), I = 1, NK )
      WRITE( NOUT, FMT = 9991 )'K:     ', ( KVAL( I ), I = 1, NK )
*
*     Find the maximum value of K (= NRHS).
*
      MAXK = 0
      DO 50 I = 1, NK
         MAXK = MAX( KVAL( I ), MAXK )
   50 CONTINUE
      MKMAX = MAXM*MAX( 2, MAXK )
*
*     Read in NNB and the values for NB.  For the BLAS input files,
*     NBVAL is used to store values for INCX and INCY.
*
      READ( NIN, FMT = * )NNB
      IF( NNB.GT.MAXVAL ) THEN
         WRITE( NOUT, FMT = 9999 )'NB', 'NNB', MAXVAL
         NNB = MAXVAL
      END IF
      READ( NIN, FMT = * )( NBVAL( I ), I = 1, NNB )
*
*     Find the maximum value of NB.
*
      MAXNB = 0
      DO 60 I = 1, NNB
         MAXNB = MAX( NBVAL( I ), MAXNB )
   60 CONTINUE
*
      IF( BLAS ) THEN
         WRITE( NOUT, FMT = 9991 )'INCX:  ', ( NBVAL( I ), I = 1, NNB )
         DO 70 I = 1, NNB
            NXVAL( I ) = 0
   70    CONTINUE
      ELSE
*
*        LAPACK data files:  Read in the values for NX.
*
         READ( NIN, FMT = * )( NXVAL( I ), I = 1, NNB )
*
         WRITE( NOUT, FMT = 9991 )'NB:    ', ( NBVAL( I ), I = 1, NNB )
         WRITE( NOUT, FMT = 9991 )'NX:    ', ( NXVAL( I ), I = 1, NNB )
      END IF
*
*     Read in NLDA and the values for LDA.
*
      READ( NIN, FMT = * )NLDA
      IF( NLDA.GT.MXNLDA ) THEN
         WRITE( NOUT, FMT = 9999 )'LDA', 'NLDA', MXNLDA
         NLDA = MXNLDA
      END IF
      READ( NIN, FMT = * )( LDAVAL( I ), I = 1, NLDA )
      WRITE( NOUT, FMT = 9991 )'LDA:   ', ( LDAVAL( I ), I = 1, NLDA )
*
*     Check that LDA >= 1 for all values of LDA.
*
      LDAOK = .TRUE.
      MAXLDA = 0
      DO 80 I = 1, NLDA
         MAXLDA = MAX( LDAVAL( I ), MAXLDA )
         IF( LDAVAL( I ).LE.0 ) THEN
            WRITE( NOUT, FMT = 9998 )LDAVAL( I )
            LDAOK = .FALSE.
         END IF
   80 CONTINUE
      IF( .NOT.LDAOK )
     $   WRITE( NOUT, FMT = * )
*
*     Check that MAXLDA*MAXN <= LA (for the dense routines).
*
      LDANOK = .TRUE.
      NEED = MAXLDA*MAXN
      IF( NEED.GT.LA ) THEN
         WRITE( NOUT, FMT = 9995 )MAXLDA, MAXN, NEED
         LDANOK = .FALSE.
      END IF
*
*     Check that MAXLDA*MAXM + MAXM*MAXK <= 3*LA (for band routines).
*
      LDAMOK = .TRUE.
      NEED = MAXLDA*MAXM + MAXM*MAXK
      IF( NEED.GT.3*LA ) THEN
         NEED = ( NEED+2 ) / 3
         WRITE( NOUT, FMT = 9994 )MAXLDA, MAXM, MAXK, NEED
         LDAMOK = .FALSE.
      END IF
*
*     Check that MAXN*MAXNB (or MAXN*INCX) <= LA.
*
      NXNBOK = .TRUE.
      NEED = MAXN*MAXNB
      IF( NEED.GT.LA ) THEN
         WRITE( NOUT, FMT = 9996 )MAXN, MAXNB, NEED
         NXNBOK = .FALSE.
      END IF
*
      IF( .NOT.( MOK .AND. NOK .AND. LDAOK .AND. LDANOK .AND. NXNBOK ) )
     $     THEN
         WRITE( NOUT, FMT = 9984 )
         GO TO 110
      END IF
      IF( .NOT.LDAMOK )
     $   WRITE( NOUT, FMT = * )
*
*     Read the minimum time to time a subroutine.
*
      WRITE( NOUT, FMT = * )
      READ( NIN, FMT = * )TIMMIN
      WRITE( NOUT, FMT = 9993 )TIMMIN
      WRITE( NOUT, FMT = * )
*
*     Read the first input line.
*
      READ( NIN, FMT = '(A)', END = 100 )LINE
*
*     If the first record is the special signal 'NONE', then get the
*     next line but don't time CGEMV and CGEMM.
*
      IF( LSAMEN( 4, LINE, 'NONE' ) ) THEN
         READ( NIN, FMT = '(A)', END = 100 )LINE
      ELSE
         WRITE( NOUT, FMT = 9990 )
*
*        If the first record is the special signal 'BAND', then time
*        the band routine CGBMV and CGEMM with N = K.
*
         IF( LSAMEN( 4, LINE, 'BAND' ) ) THEN
            IF( LDAMOK ) THEN
               IF( MKMAX.GT.LA ) THEN
                  I2 = 2*LA - MKMAX + 1
                  J2 = 2
               ELSE
                  I2 = LA - MKMAX + 1
                  J2 = 3
               END IF
               CALL CTIMMV( 'CGBMV ', NM, MVAL, NN, NVAL, NLDA, LDAVAL,
     $                      TIMMIN, A( 1, 1 ), MKMAX / 2, A( I2, J2 ),
     $                      A( LA-MKMAX / 2+1, 3 ), RESLTS, LDR1, LDR2,
     $                      NOUT )
            ELSE
               WRITE( NOUT, FMT = 9989 )'CGBMV '
            END IF
            CALL CTIMMM( 'CGEMM ', 'K', NN, NVAL, NLDA, LDAVAL, TIMMIN,
     $                   A( 1, 1 ), A( 1, 2 ), A( 1, 3 ), RESLTS, LDR1,
     $                   LDR2, NOUT )
            READ( NIN, FMT = '(A)', END = 100 )LINE
*
         ELSE
*
*           Otherwise time CGEMV and CGEMM.
*
            CALL CTIMMV( 'CGEMV ', NN, NVAL, NNB, NBVAL, NLDA, LDAVAL,
     $                   TIMMIN, A( 1, 1 ), LA, A( 1, 2 ), A( 1, 3 ),
     $                   RESLTS, LDR1, LDR2, NOUT )
            CALL CTIMMM( 'CGEMM ', 'N', NN, NVAL, NLDA, LDAVAL, TIMMIN,
     $                   A( 1, 1 ), A( 1, 2 ), A( 1, 3 ), RESLTS, LDR1,
     $                   LDR2, NOUT )
         END IF
      END IF
*
*     Call the appropriate timing routine for each input line.
*
      WRITE( NOUT, FMT = 9988 )
   90 CONTINUE
      C1 = LINE( 1: 1 )
      C2 = LINE( 2: 3 )
      C3 = LINE( 4: 6 )
*
*     Check first character for correct precision.
*
      IF( .NOT.LSAME( C1, 'Complex precision' ) ) THEN
         WRITE( NOUT, FMT = 9987 )LINE( 1: 6 )
*
      ELSE IF( LSAMEN( 2, C2, 'B2' ) .OR. LSAMEN( 3, C3, 'MV ' ) .OR.
     $         LSAMEN( 3, C3, 'SV ' ) .OR. LSAMEN( 3, C3, 'R  ' ) .OR.
     $         LSAMEN( 3, C3, 'RC ' ) .OR. LSAMEN( 3, C3, 'RU ' ) .OR.
     $         LSAMEN( 3, C3, 'R2 ' ) ) THEN
*
*        Level 2 BLAS
*
         CALL CTIMB2( LINE, NM, MVAL, NN, NVAL, NK, KVAL, NNB, NBVAL,
     $                NLDA, LDAVAL, LA, TIMMIN, A( 1, 1 ), A( 1, 2 ),
     $                A( 1, 3 ), RESLTS, LDR1, LDR2, NOUT )
*
      ELSE IF( LSAMEN( 2, C2, 'B3' ) .OR. LSAMEN( 3, C3, 'MM ' ) .OR.
     $         LSAMEN( 3, C3, 'SM ' ) .OR. LSAMEN( 3, C3, 'RK ' ) .OR.
     $         LSAMEN( 3, C3, 'R2K' ) ) THEN
*
*        Level 3 BLAS
*
         CALL CTIMB3( LINE, NM, MVAL, NN, NVAL, NK, KVAL, NLDA, LDAVAL,
     $                TIMMIN, A( 1, 1 ), A( 1, 2 ), A( 1, 3 ), RESLTS,
     $                LDR1, LDR2, NOUT )
*
      ELSE IF( LSAMEN( 2, C2, 'QR' ) .OR. LSAMEN( 2, C3, 'QR' ) .OR.
     $         LSAMEN( 2, C3( 2: 3 ), 'QR' ) ) THEN
*
*        QR routines
*
         CALL CTIMQR( LINE, NN, MVAL, NVAL, NK, KVAL, NNB, NBVAL, NXVAL,
     $                NLDA, LDAVAL, TIMMIN, A( 1, 1 ), E, A( 1, 2 ),
     $                A( 1, 3 ), D, RESLTS, LDR1, LDR2, LDR3, NOUT )
*
      ELSE IF( LSAMEN( 2, C2, 'LQ' ) .OR. LSAMEN( 2, C3, 'LQ' ) .OR.
     $         LSAMEN( 2, C3( 2: 3 ), 'LQ' ) ) THEN
*
*        LQ routines
*
         CALL CTIMLQ( LINE, NN, MVAL, NVAL, NK, KVAL, NNB, NBVAL, NXVAL,
     $                NLDA, LDAVAL, TIMMIN, A( 1, 1 ), E, A( 1, 2 ),
     $                A( 1, 3 ), D, RESLTS, LDR1, LDR2, LDR3, NOUT )
*
      ELSE IF( LSAMEN( 2, C2, 'QL' ) .OR. LSAMEN( 2, C3, 'QL' ) .OR.
     $         LSAMEN( 2, C3( 2: 3 ), 'QL' ) ) THEN
*
*        QL routines
*
         CALL CTIMQL( LINE, NN, MVAL, NVAL, NK, KVAL, NNB, NBVAL, NXVAL,
     $                NLDA, LDAVAL, TIMMIN, A( 1, 1 ), E, A( 1, 2 ),
     $                A( 1, 3 ), D, RESLTS, LDR1, LDR2, LDR3, NOUT )
*
      ELSE IF( LSAMEN( 2, C2, 'RQ' ) .OR. LSAMEN( 2, C3, 'RQ' ) .OR.
     $         LSAMEN( 2, C3( 2: 3 ), 'RQ' ) ) THEN
*
*        RQ routines
*
         CALL CTIMRQ( LINE, NN, MVAL, NVAL, NK, KVAL, NNB, NBVAL, NXVAL,
     $                NLDA, LDAVAL, TIMMIN, A( 1, 1 ), E, A( 1, 2 ),
     $                A( 1, 3 ), D, RESLTS, LDR1, LDR2, LDR3, NOUT )
*
      ELSE IF( LSAMEN( 2, C2, 'QP' ) .OR. LSAMEN( 3, C3, 'QPF' ) ) THEN
*
*        QR with column pivoting
*
         CALL CTIMQP( LINE, NM, MVAL, NVAL, NLDA, LDAVAL, TIMMIN,
     $                A( 1, 1 ), A( 1, 2 ), E, A( 1, 3 ), D, IWORK,
     $                RESLTS, LDR1, LDR2, NOUT )
*
*        Rank-Revealing QR factorization
*
         CALL CTIMQ3( LINE, NM, MVAL, NVAL, NNB, NBVAL, NXVAL, NLDA,
     $                LDAVAL, TIMMIN, A( 1, 1 ), A( 1, 2 ), E,
     $                A( 1, 3 ), D, IWORK, RESLTS, LDR1, LDR2, NOUT )
*
      ELSE IF( LSAMEN( 2, C2, 'HR' ) .OR. LSAMEN( 3, C3, 'HRD' ) .OR.
     $         LSAMEN( 2, C3( 2: 3 ), 'HR' ) ) THEN
*
*        Reduction to Hessenberg form
*
         CALL CTIMHR( LINE, NN, NVAL, NK, KVAL, NNB, NBVAL, NXVAL, NLDA,
     $                LDAVAL, TIMMIN, A( 1, 1 ), E, A( 1, 2 ),
     $                A( 1, 3 ), D, RESLTS, LDR1, LDR2, LDR3, NOUT )
*
      ELSE IF( LSAMEN( 2, C2, 'TD' ) .OR. LSAMEN( 3, C3, 'TRD' ) .OR.
     $         LSAMEN( 2, C3( 2: 3 ), 'TR' ) ) THEN
*
*        Reduction to tridiagonal form
*
         CALL CTIMTD( LINE, NN, NVAL, NK, KVAL, NNB, NBVAL, NXVAL, NLDA,
     $                LDAVAL, TIMMIN, A( 1, 1 ), A( 1, 2 ), D, E,
     $                A( 1, 3 ), RESLTS, LDR1, LDR2, LDR3, NOUT )
*
      ELSE IF( LSAMEN( 2, C2, 'BR' ) .OR. LSAMEN( 3, C3, 'BRD' ) .OR.
     $         LSAMEN( 2, C3( 2: 3 ), 'BR' ) ) THEN
*
*        Reduction to bidiagonal form
*
         CALL CTIMBR( LINE, NN, MVAL, NVAL, NK, KVAL, NNB, NBVAL, NXVAL,
     $                NLDA, LDAVAL, TIMMIN, A( 1, 1 ), A( 1, 2 ), D, E,
     $                A( 1, 3 ), RESLTS, LDR1, LDR2, LDR3, NOUT )
*
      ELSE IF( LSAMEN( 2, C2, 'GE' ) ) THEN
*
*        Routines for general matrices
*
         CALL CTIMGE( LINE, NN, NVAL, NK, KVAL, NNB, NBVAL, NLDA,
     $                LDAVAL, TIMMIN, A( 1, 1 ), A( 1, 2 ), A( 1, 3 ),
     $                IWORK, RESLTS, LDR1, LDR2, LDR3, NOUT )
*
      ELSE IF( LSAMEN( 2, C2, 'GB' ) ) THEN
*
*        General band matrices
*
         IF( LDAMOK ) THEN
            CALL CTIMGB( LINE, NM, MVAL, NN, NVAL, NK, KVAL, NNB, NBVAL,
     $                   NLDA, LDAVAL, TIMMIN, A( 1, 1 ),
     $                   A( LA-MKMAX+1, 3 ), IWORK, RESLTS, LDR1, LDR2,
     $                   LDR3, NOUT )
         ELSE
            WRITE( NOUT, FMT = 9989 )LINE( 1: 6 )
         END IF
*
      ELSE IF( LSAMEN( 2, C2, 'GT' ) ) THEN
*
*        Routines for general tridiagonal matrices
*
         CALL CTIMGT( LINE, NN, NVAL, NK, KVAL, NLDA, LDAVAL, TIMMIN,
     $                A( 1, 1 ), A( 1, 2 ), IWORK, RESLTS, LDR1, LDR2,
     $                LDR3, NOUT )
*
      ELSE IF( LSAMEN( 2, C2, 'PO' ) ) THEN
*
*        Positive definite matrices
*
         CALL CTIMPO( LINE, NN, NVAL, NK, KVAL, NNB, NBVAL, NLDA,
     $                LDAVAL, TIMMIN, A( 1, 1 ), A( 1, 2 ), IWORK,
     $                RESLTS, LDR1, LDR2, LDR3, NOUT )
*
      ELSE IF( LSAMEN( 2, C2, 'PP' ) ) THEN
*
*        Positive definite packed matrices
*
         CALL CTIMPP( LINE, NN, NVAL, NK, KVAL, LA, TIMMIN, A( 1, 1 ),
     $                A( 1, 2 ), IWORK, RESLTS, LDR1, LDR2, LDR3, NOUT )
*
      ELSE IF( LSAMEN( 2, C2, 'PB' ) ) THEN
*
*        Positive definite banded matrices
*
         IF( LDAMOK ) THEN
            IF( MKMAX.GT.LA ) THEN
               J2 = 2
               I2 = 2*LA - MKMAX + 1
            ELSE
               J2 = 3
               I2 = LA - MKMAX + 1
            END IF
            CALL CTIMPB( LINE, NM, MVAL, NN, NVAL, NK, KVAL, NNB, NBVAL,
     $                   NLDA, LDAVAL, TIMMIN, A( 1, 1 ), A( I2, J2 ),
     $                   IWORK, RESLTS, LDR1, LDR2, LDR3, NOUT )
         ELSE
            WRITE( NOUT, FMT = 9989 )LINE( 1: 6 )
         END IF
*
      ELSE IF( LSAMEN( 2, C2, 'PT' ) ) THEN
*
*        Routines for positive definite tridiagonal matrices
*
         CALL CTIMPT( LINE, NN, NVAL, NK, KVAL, NLDA, LDAVAL, TIMMIN, D,
     $                A( 1, 1 ), A( 1, 2 ), RESLTS, LDR1, LDR2, LDR3,
     $                NOUT )
*
      ELSE IF( LSAMEN( 2, C2, 'HE' ) ) THEN
*
*        Hermitian indefinite matrices
*
         CALL CTIMHE( LINE, NN, NVAL, NK, KVAL, NNB, NBVAL, NLDA,
     $                LDAVAL, TIMMIN, A( 1, 1 ), A( 1, 2 ), A( 1, 3 ),
     $                IWORK, RESLTS, LDR1, LDR2, LDR3, NOUT )
*
      ELSE IF( LSAMEN( 2, C2, 'HP' ) ) THEN
*
*        Hermitian indefinite packed matrices
*
         CALL CTIMHP( LINE, NN, NVAL, NK, KVAL, LA, TIMMIN, A( 1, 1 ),
     $                A( 1, 2 ), A( 1, 3 ), IWORK, RESLTS, LDR1, LDR2,
     $                LDR3, NOUT )
*
      ELSE IF( LSAMEN( 2, C2, 'SY' ) ) THEN
*
*        Symmetric indefinite matrices
*
         CALL CTIMSY( LINE, NN, NVAL, NK, KVAL, NNB, NBVAL, NLDA,
     $                LDAVAL, TIMMIN, A( 1, 1 ), A( 1, 2 ), A( 1, 3 ),
     $                IWORK, RESLTS, LDR1, LDR2, LDR3, NOUT )
*
      ELSE IF( LSAMEN( 2, C2, 'SP' ) ) THEN
*
*        Symmetric indefinite packed matrices
*
         CALL CTIMSP( LINE, NN, NVAL, NK, KVAL, LA, TIMMIN, A( 1, 1 ),
     $                A( 1, 2 ), A( 1, 3 ), IWORK, RESLTS, LDR1, LDR2,
     $                LDR3, NOUT )
*
      ELSE IF( LSAMEN( 2, C2, 'TR' ) ) THEN
*
*        Triangular matrices
*
         CALL CTIMTR( LINE, NN, NVAL, NK, KVAL, NNB, NBVAL, NLDA,
     $                LDAVAL, TIMMIN, A( 1, 1 ), A( 1, 2 ), RESLTS,
     $                LDR1, LDR2, LDR3, NOUT )
*
      ELSE IF( LSAMEN( 2, C2, 'TP' ) ) THEN
*
*        Triangular packed matrices
*
         CALL CTIMTP( LINE, NN, NVAL, NK, KVAL, LA, TIMMIN, A( 1, 1 ),
     $                A( 1, 2 ), RESLTS, LDR1, LDR2, LDR3, NOUT )
*
      ELSE IF( LSAMEN( 2, C2, 'TB' ) ) THEN
*
*        Triangular band matrices
*
         IF( LDAMOK ) THEN
            IF( MKMAX.GT.LA ) THEN
               J2 = 2
               I2 = 2*LA - MKMAX + 1
            ELSE
               J2 = 3
               I2 = LA - MKMAX + 1
            END IF
            CALL CTIMTB( LINE, NM, MVAL, NN, NVAL, NK, KVAL, NLDA,
     $                   LDAVAL, TIMMIN, A( 1, 1 ), A( I2, J2 ), RESLTS,
     $                   LDR1, LDR2, LDR3, NOUT )
         ELSE
            WRITE( NOUT, FMT = 9989 )LINE( 1: 6 )
         END IF
*
      ELSE IF( LSAMEN( 2, C2, 'LS' ) ) THEN
*
*        Least squares drivers
*
         CALL CTIMLS( LINE, NM, MVAL, NN, NVAL, NK, KVAL, NNB, NBVAL,
     $                NXVAL, NLDA, LDAVAL, TIMMIN, A( 1, 1 ), A( 1, 2 ),
     $                B( 1, 1 ), B( 1, 2 ), S, S( NMAX+1 ), OPCTBL,
     $                TIMTBL, FLPTBL, WORK, RWORK, IWORK, NOUT )
*
      ELSE
*
         WRITE( NOUT, FMT = 9987 )LINE( 1: 6 )
      END IF
*
*     Read the next line of the input file.
*
      READ( NIN, FMT = '(A)', END = 100 )LINE
      GO TO 90
*
*     Branch to this line when the last record is read.
*
  100 CONTINUE
      S2 = SECOND( )
      WRITE( NOUT, FMT = 9986 )
      WRITE( NOUT, FMT = 9985 )S2 - S1
  110 CONTINUE
*
 9999 FORMAT( ' Too many values of ', A, ' using ', A, ' = ', I2 )
 9998 FORMAT( ' *** LDA = ', I7, ' is too small, must have ',
     $      'LDA > 0.' )
 9997 FORMAT( ' *** ', A1, ' = ', I7, ' is too big:  ',
     $      'maximum allowed is', I7 )
 9996 FORMAT( ' *** N*NB is too big for N =', I6, ', NB =', I6,
     $      / ' --> Increase LA to at least ', I8 )
 9995 FORMAT( ' *** LDA*N is too big for the dense routines ', '(LDA =',
     $      I6, ', N =', I6, ')', / ' --> Increase LA to at least ',
     $      I8 )
 9994 FORMAT( ' *** (LDA+K)*M is too big for the band routines ',
     $      '(LDA=', I6, ', M=', I6, ', K=', I6, ')',
     $      / ' --> Increase LA to at least ', I8 )
 9993 FORMAT( ' The minimum time a subroutine will be timed = ', F6.3,
     $      ' seconds' )
 9992 FORMAT( ' The following parameter values will be used:' )
 9991 FORMAT( 4X, A7, 1X, 10I6, / 12X, 10I6 )
 9990 FORMAT( / ' ------------------------------',
     $      / ' >>>>>    Sample BLAS     <<<<<',
     $      / ' ------------------------------' )
 9989 FORMAT( 1X, A6, ' not timed due to input errors', / )
 9988 FORMAT( / ' ------------------------------',
     $      / ' >>>>>    Timing data     <<<<<',
     $      / ' ------------------------------' )
 9987 FORMAT( 1X, A6, ':  Unrecognized path or subroutine name', / )
 9986 FORMAT( ' End of tests' )
 9985 FORMAT( ' Total time used = ', F12.2, ' seconds' )
 9984 FORMAT( / ' Tests not done due to input errors' )
 9983 FORMAT( ' LAPACK VERSION 3.0, released June 30, 1999 ', / )
*
*     End of CTIMAA
*
      END

max_stars_repo_name: haampie/libflame

max_issues_repo_name: haampie/libflame

max_forks_repo_name: haampie/libflame