text:
cdis    Forecast Systems Laboratory
cdis    NOAA/OAR/ERL/FSL
cdis    325 Broadway
cdis    Boulder, CO     80303
cdis
cdis    Forecast Research Division
cdis    Local Analysis and Prediction Branch
cdis    LAPS
cdis
cdis    This software and its documentation are in the public domain and
cdis    are furnished "as is."  The United States government, its
cdis    instrumentalities, officers, employees, and agents make no
cdis    warranty, express or implied, as to the usefulness of the software
cdis    and documentation for any purpose.  They assume no responsibility
cdis    (1) for the use of the software and documentation; or (2) to provide
cdis     technical support to users.
cdis
cdis    Permission to use, copy, modify, and distribute this software is
cdis    hereby granted, provided that the entire disclaimer notice appears
cdis    in all copies.  All modifications to this software must be clearly
cdis    documented, and are solely the responsibility of the agent making
cdis    the modifications.  If significant modifications or enhancements
cdis    are made to this software, the FSL Software Policy Manager
cdis    (softwaremgr@fsl.noaa.gov) should be notified.
cdis
cdis
cdis
cdis
cdis
cdis
cdis
C
        SUBROUTINE i4time_fname_lp (FNAME_IN, I4TIME, ISTATUS)

        include 'lapsparms.for'
C
cdoc    This routine converts several different file name types (i.e. yydddhhmm)
cdoc    into the corresponding I4 time.
C
C       ON INPUT
C          FILE_NAME - The file name including the directory path if wanted.
C
C       ON OUTPUT
C          I4TIME - The corresponding I4 time of the file name.
C          ISTATUS -  The return status.
C
C================================================================
C
        CHARACTER*9 wfo_fname13_to_fname9, rsa13_to_a9, a8_to_a9
        CHARACTER*9 a7_to_a9_time,yr_a10_to_a9
        CHARACTER*9 FILE_NAME
        CHARACTER*(*) FNAME_IN
        CHARACTER*256 FNAME_BUF
        CHARACTER*20 c20_type
        INTEGER I4TIME, ISTATUS
C
        INTEGER INT_FILE(9), I, NYEAR, JDAY, NHOUR, MIN, MONTH, NDAY
        INTEGER I4TIME_INT_LP
C
C================================================================
C
C       Check the length of FILE_NAME to be sure that it is valid.
C
        FNAME_BUF = FNAME_IN
        call get_filetime_type(fname_in,c20_type,leni,lent)

        if(c20_type .eq. 'yyjjjhhmm')then                   ! NIMBUS/LAPS type
           file_name = fname_in(leni+1:leni+lent)

        elseif(c20_type .eq. 'yyyymmdd_hhmm')then           ! WFO type
           file_name = wfo_fname13_to_fname9(fname_in(leni+1:leni+lent))

        elseif(c20_type .eq. 'yyyyjjjhhmmss')then           ! RSA type
           file_name = rsa13_to_a9(fname_in(leni+1:leni+lent))

        elseif(c20_type .eq. 'yymmddhh')then                ! AFWA RAOB type
           file_name = a8_to_a9(fname_in(leni+1:leni+lent))

        elseif(c20_type .eq. 'ymmddhh')then               ! Taiwan FA Model
           file_name = a7_to_a9_time(fname_in(leni+1:leni+lent))

        elseif(c20_type .eq. 'yyyymmddhh')then
           file_name = yr_a10_to_a9(fname_in(leni+1:leni+lent))

        else                                                ! Unrecognized type
c          write(6,*)'i4time_fname_lp: unable to convert to i4time',
c    1'    type = ',c20_type
           istatus = 0
           FNAME_IN = FNAME_BUF
           return
        endif
C
C       Split the file name into individual integers, checking for invalid
C       characters in the file name.
C
        DO I = 1,9
           INT_FILE(I) = ICHAR(FILE_NAME(I:I))-48
           IF (INT_FILE(I) .LT. 0 .OR. INT_FILE(I) .GT.9) THEN
              GO TO 1500
           END IF
        END DO
C
C       Convert these numbers in year, day of year (julian day), hours,
C       and minutes.
        NYEAR = 10*INT_FILE(1) + INT_FILE(2)
        JDAY = 100*INT_FILE(3) + 10*INT_FILE(4) + INT_FILE(5)
        NHOUR = 10*INT_FILE(6) + INT_FILE(7)
        MIN = 10*INT_FILE(8) + INT_FILE(9)

C       Get the actual year, here is where we assume what century it is...

        iyear_earliest_century = (iyear_earliest/100) * 100
        iyy_cutoff = iyear_earliest - iyear_earliest_century

        if(NYEAR .lt. iyy_cutoff)then
            NYEAR = NYEAR + iyear_earliest_century + 100
        else
            NYEAR = NYEAR + iyear_earliest_century
        endif
C
C       Convert the day of year (julian day) into month and day.
C
        CALL CV_JUL_MMDD_LP (JDAY, NYEAR, MONTH, NDAY, ISTATUS)
        IF (0 .EQ. ISTATUS)
     1    GO TO 2000
C
C       Convert these integers into I4 time.
C
        I4TIME = I4TIME_INT_LP (NYEAR, MONTH, NDAY, NHOUR, MIN, 0, ISTAT
     1US)
        IF (0 .EQ. ISTATUS)
     1    GO TO 2000
C
C       Successful return
C
        ISTATUS = 1
        RETURN
C
C       Error detected in this routine.
C
 1500   CONTINUE
        ISTATUS = 0
        WRITE( 6,* ) 'Error in i4time_fname_lp: bad digit in file name.'
        WRITE( 6,* ) FILE_NAME
        RETURN
C
C       Error detected in called routine
C
 2000   CONTINUE
        ISTATUS = 0
        WRITE( 6,* ) 'Error in i4time_fname_lp: error in subroutine.'
        WRITE( 6,* ) FILE_NAME
        RETURN

 1000   CONTINUE
        ISTATUS = 0
        WRITE( 6,* ) 'Error in i4time_fname_lp: wrong length for file na
     1me.'
        WRITE( 6,* ) FILE_NAME

        RETURN
C
        END
C
        SUBROUTINE CV_JUL_MMDD_LP (JULIAN_DAY, YEAR, MONTH, DAY, ISTATUS
     1)
C
cdoc    This routine converts from day of year (Julian days) to month and
cdoc    day in integer format.
C
C       ON INPUT
C          JULIAN_DAY - The Day of Year (Julian date) to be converted.
C          YEAR - The year of the Day of Year (Julian date), if .lt. 100,
C                 it is assumed to be the last two digits of 19xx.
C
C       ON OUTPUT
C          MONTH - The integer representation for the month (1-12)
C          DAY - The integer value of the day of the month.
C          ISTATUS - The return status.
C
C================================================================
C
        INTEGER JULIAN_DAY, YEAR, DAY, ISTATUS
C
        INTEGER MNTH(12), TEMP_YEAR, MONTH, MAX_DAY
        LOGICAL*1 LEAP
        DATA MNTH/31,0,31,30,31,30,31,31,30,31,30,31/
C
C================================================================
C
        IF (YEAR.LT.100) THEN
           TEMP_YEAR = YEAR + 1900
        ELSE
           TEMP_YEAR = YEAR
        ENDIF
C
C       Check to see if the year is a leap year
C
        LEAP = ( MOD(TEMP_YEAR,4)  .EQ. 0 )
C
        IF ( LEAP ) THEN
           MAX_DAY = 366
           MNTH(2) = 29
        ELSE
           MAX_DAY = 365
           MNTH(2) = 28
        END IF
C
C       Check day of year (JULIAN_DAY) for being too large or too small
C
        IF (JULIAN_DAY .LE. 0 .OR. JULIAN_DAY .GT. MAX_DAY) THEN
           ISTATUS = 0
           WRITE( 6,* ) 'Error in cv_jul_mmdd_LP.f: bad Julian day', JUL
     1IAN_DAY
           RETURN
        END IF
C
C       Convert to month and day
C
        MONTH = 1
        DAY = JULIAN_DAY
        DO WHILE (DAY .GT. MNTH(MONTH))
           DAY = DAY - MNTH(MONTH)
           MONTH = MONTH + 1
        END DO
        ISTATUS = 1
C
        RETURN
        END

C
        FUNCTION I4TIME_INT_LP (NYEAR,NMONTH,NDAY,NHOUR,NMIN,NSEC
     1                         ,ISTATUS)
C
cdoc    I4TIME_INT_LP RETURNS I4 TIME (# OF SECONDS SINCE 00:00 01-JAN-60)
cdoc    GIVEN A 6 INTEGERS CONTAINING YEAR, MONTH, DAY, HOUR,
cdoc    MINUTE, SECOND
C
C================================================================
C
        INTEGER I4TIME_INT_LP
        INTEGER NYEAR, NMONTH, NDAY, NHOUR, NMIN, NSEC, ISTATUS
C
        INTEGER NSECYR, NSECDA, NSECHR, NSECMN
        INTEGER NYR, NYRS, NLEAP
        INTEGER NDAYS(12), NSECMO(12)
        INTEGER IBASE, ISUM, I
        PARAMETER (IBASE=60)

        DATA NSECYR/31536000/NSECDA/86400/NSECHR/3600/NSECMN/60/
        DATA NSECMO/2678400,2419200,2678400,2592000,2678400,
     1  2592000,2678400,2678400,2592000,2678400,2592000,2678400/

        DATA NDAYS/31,29,31,30,31,30,31,31,30,31,30,31/
C
C================================================================
C
        ISTATUS = 1
C
C Sum the number of years.
C
        NYR=NYEAR
        IF (NYR.GT.1900) NYR=NYR-1900
        NYRS=NYR-IBASE
        IF (NYRS.LT.0.OR.NYRS.GT.67) GO TO 1000
        ISUM=NYRS*NSECYR
        NLEAP=(NYRS/4)+1                        ! Account for leap years.
        ISUM=ISUM+(NLEAP*NSECDA)
C
C Sum in the number of months.
C
        IF (NMONTH.LT.1.OR.NMONTH.GT.12) GO TO 1000
C
        IF (NMONTH.NE.1) THEN
                DO I=1,NMONTH-1
                        ISUM=ISUM + NSECMO(I)
                END DO
        END IF
C
C Correct for Jan or Feb of a leap yer.
C
!       IF (MOD(NYRS,4).EQ.0) THEN
        IF (MOD(NYR,4).EQ.0) THEN       ! Steve Albers, Linda Wharton 1993
                IF (NMONTH.LE.2) ISUM=ISUM-NSECDA
        END IF
C
C Sum in the number of days.
C
        IF (NDAY.LT.1.OR.NDAY.GT.NDAYS(NMONTH))
     1  GO TO 1000
        ISUM=ISUM+((NDAY-1)*NSECDA)
C
C Sum in the number of hours.
C
        IF (NHOUR.LT.0.OR.NHOUR.GT.23) GO TO 1000
        ISUM=ISUM+(NHOUR*NSECHR)
C
C Sum in the number of minutes.
C
        IF (NMIN.LT.0.OR.NMIN.GT.59) GO TO 1000
        ISUM=ISUM+(NMIN*NSECMN)
C
C Sum in the number of seconds.
C
        IF (NSEC.LT.0.OR.NSEC.GT.59) GO TO 1000
        ISUM=ISUM+NSEC
C
        I4TIME_INT_LP=ISUM
        RETURN
C
C Return error code.
C
1000    ISTATUS = 0
        I4TIME_INT_LP=0
        WRITE(6,*) 'Error in I4TIME_INT_LP'
        RETURN
C
        END

ext: f

max_stars_repo_name: maxinye/laps-mirror

max_stars_repo_licenses: Intel Unlicense OLDAP-2.2.1 NetCDF

max_issues_repo_name: longwosion/laps-mirror

max_issues_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense

max_forks_repo_name: longwosion/laps-mirror

max_forks_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense


text:
      subroutine evolve_dph(ctfmp,n_K,n_alpha,n_beta,n_divbeta,n_phi,np1
     &_K,np1_alpha,np1_beta,np1_divbeta,np1_phi,x,Nx,ht,hx,myzero,zepsdi
     &s,phys_bdy,res)
      implicit none
      integer i
      integer Nx
      real*8 ht
      real*8 hx
      real*8 myzero
      real*8 zepsdis
      real*8 ctfmp(Nx)
      real*8 n_K(Nx)
      real*8 n_alpha(Nx)
      real*8 n_beta(Nx)
      real*8 n_divbeta(Nx)
      real*8 n_phi(Nx)
      real*8 np1_K(Nx)
      real*8 np1_alpha(Nx)
      real*8 np1_beta(Nx)
      real*8 np1_divbeta(Nx)
      real*8 np1_phi(Nx)
      real*8 x(Nx)
      integer phys_bdy(2)
      real*8 res(Nx)
      real*8 qb
      if (phys_bdy(1) .eq. 1) then
      do i=1, 1, 1
      qb = 0.1333333333333333D1 * np1_phi(i + 1) - 0.3333333333333333D0
     #* np1_phi(i + 2)
      res(i)=qb
      end do
      endif
      if (phys_bdy(1) .eq. 1) then
      do i=2, 2, 1
      qb = np1_phi(i) - 0.32D2 / (0.32D2 + 0.7D1 * zepsdis) * ht * ((-0.
     #1D1 * n_phi(i) + np1_phi(i)) / ht + 0.2500000000000000D0 / ctfmp(i
     #) * (np1_phi(i - 1) - 0.1D1 * np1_phi(i + 1)) / hx * np1_beta(i) +
     # 0.8333333333333333D-1 * np1_alpha(i) * np1_K(i) - 0.8333333333333
     #333D-1 * np1_divbeta(i) + 0.2500000000000000D0 / ctfmp(i) * (n_phi
     #(i - 1) - 0.1D1 * n_phi(i + 1)) / hx * n_beta(i) + 0.8333333333333
     #333D-1 * n_alpha(i) * n_K(i) - 0.8333333333333333D-1 * n_divbeta(i
     #) + 0.3125000000000000D-1 * zepsdis / ht * (0.7D1 * np1_phi(i) + n
     #p1_phi(i + 2) - 0.4D1 * np1_phi(i + 1) - 0.4D1 * np1_phi(i - 1)) +
     # 0.3125000000000000D-1 * zepsdis / ht * (0.7D1 * n_phi(i) + n_phi(
     #i + 2) - 0.4D1 * n_phi(i + 1) - 0.4D1 * n_phi(i - 1)))
      res(i)=qb
      end do
      endif
      do i=3, Nx-2, 1
      qb = np1_phi(i) - 0.16D2 / (0.16D2 + 0.3D1 * zepsdis) * ht * ((-0.
     #1D1 * n_phi(i) + np1_phi(i)) / ht + 0.2500000000000000D0 / ctfmp(i
     #) * (np1_phi(i - 1) - 0.1D1 * np1_phi(i + 1)) / hx * np1_beta(i) +
     # 0.8333333333333333D-1 * np1_alpha(i) * np1_K(i) - 0.8333333333333
     #333D-1 * np1_divbeta(i) + 0.2500000000000000D0 / ctfmp(i) * (n_phi
     #(i - 1) - 0.1D1 * n_phi(i + 1)) / hx * n_beta(i) + 0.8333333333333
     #333D-1 * n_alpha(i) * n_K(i) - 0.8333333333333333D-1 * n_divbeta(i
     #) + 0.3125000000000000D-1 * zepsdis / ht * (0.6D1 * np1_phi(i) + n
     #p1_phi(i + 2) + np1_phi(i - 2) - 0.4D1 * np1_phi(i + 1) - 0.4D1 *
     #np1_phi(i - 1)) + 0.3125000000000000D-1 * zepsdis / ht * (0.6D1 *
     #n_phi(i) + n_phi(i + 2) + n_phi(i - 2) - 0.4D1 * n_phi(i + 1) - 0.
     #4D1 * n_phi(i - 1)))
      res(i)=qb
      end do
      do i=Nx-1, Nx-1, 1
      qb = np1_phi(i) - 0.1D1 * ht * ((-0.1D1 * n_phi(i) + np1_phi(i)) /
     # ht + 0.2500000000000000D0 / ctfmp(i) * (np1_phi(i - 1) - 0.1D1 *
     #np1_phi(i + 1)) / hx * np1_beta(i) + 0.8333333333333333D-1 * np1_a
     #lpha(i) * np1_K(i) - 0.8333333333333333D-1 * np1_divbeta(i) + 0.25
     #00000000000000D0 / ctfmp(i) * (n_phi(i - 1) - 0.1D1 * n_phi(i + 1)
     #) / hx * n_beta(i) + 0.8333333333333333D-1 * n_alpha(i) * n_K(i) -
     # 0.8333333333333333D-1 * n_divbeta(i))
      res(i)=qb
      end do
      if (phys_bdy(2) .eq. 1) then
      do i=Nx, Nx, 1
      qb = np1_phi(i) - 0.1D1 * ht * ((-0.1D1 * n_phi(i) + np1_phi(i)) /
     # ht - 0.1D1 * myzero * x(i))
      res(i)=qb
      end do
      endif
      END


text:
! Base class for Gaussian process types

module m_gp
  use m_util
  use m_cov
  use m_noise

  private
  public BaseGP, nlog_lik, set_hyperparams

  type, abstract :: BaseGP
     ! The noise hyperparameter(s). This is passed to the noise model
     ! (`noise_model'), which determines its precise meaning.
     real(dp), dimension(:), allocatable :: nu
     ! covariance hyperparameters
     ! meaning depends on covariance function `covariance'
     real(dp), dimension(:), allocatable :: theta
     ! inputs
     real(dp), dimension(:,:), allocatable :: x
     ! types of the observations
     integer, dimension(:), allocatable :: obs_type
     ! observations
     real(dp), dimension(:), allocatable :: t
     ! the covariance function
     class(cov_fn), allocatable :: covariance
     ! the noise model
     class(noise_model), allocatable :: noise_model
   contains
     procedure(log_lik), deferred :: log_lik
     procedure(update_matrices), deferred :: update_matrices
     procedure(predict), deferred :: predict
     procedure(write_out), deferred :: write_out
  end type BaseGP

  abstract interface
     ! The log likelihood of the hyperparameters
     function log_lik(this)
       use m_util, only: dp
       import BaseGP
       class(BaseGP), intent(in) :: this
       real(dp) log_lik
     end function log_lik

     ! Helper routine to update the internal state.  Called when an
     ! observation or hyperparameter changes and the covariance matrix
     ! must be recomputed.
     subroutine update_matrices(this)
       import BaseGP
       class(BaseGP), intent(inout) :: this
     end subroutine update_matrices

     ! Make a prediction of the underlying function value at
     ! coordinate `xnew'.
     function predict(this, xnew, obs_type_new)
       use m_util, only: dp
       import BaseGP
       class(BaseGP), intent(in) :: this
       real(dp) predict
       real(dp), dimension(:), intent(in) :: xnew
       integer, optional, intent(in) :: obs_type_new
     end function predict

     ! Serialize to a file
     subroutine write_out(this, filename)
       import BaseGP
       class(BaseGP), intent(in) :: this
       character(len=*), intent(in) :: filename
     end subroutine write_out
  end interface

contains

  ! Set the hyperparameters of `gp' to `hypers', assuming they are
  ! ordered [nu(1:nnu), theta(1:ntheta)], where nu are the noise
  ! hyperparameters and theta are the covariance hyperparameters.
  ! Calls `update_matrices'.
  subroutine set_hyperparams(gp, hypers)
    class(BaseGP), intent(inout) :: gp
    real(dp), dimension(:) :: hypers
    integer nnu
    integer ntheta
    nnu = size(gp%nu)
    ntheta = size(gp%theta)

    if ((any(hypers(1:nnu).ne.gp%nu)) &
         & .or.any(hypers(nnu+1:nnu+ntheta).ne.gp%theta)) then
       gp%nu = hypers(1:nnu)
       gp%theta = hypers(nnu+1:nnu+ntheta)
       call gp%update_matrices
    end if
  end subroutine set_hyperparams

  ! A routine with the required interface for the NLopt library, for
  ! maximizing the log-likelihood.
  subroutine nlog_lik(val, n, hypers, grad, need_gradient, gp)
    class(BaseGP) :: gp
    integer :: n, need_gradient
    real(dp) :: val, hypers(n)
    real(dp), intent(inout) :: grad(n)

    call set_hyperparams(gp, hypers)
    val = gp%log_lik()

    call output_params(gp,val)
  end subroutine nlog_lik

  ! Helper routine for nlog_lik: print out the current hyperparameters
  ! (associated with `gp') and their corresponding log-likelihood
  ! (`val').
  subroutine output_params(gp,val)
    class(BaseGP) :: gp
    real(dp) :: val
    integer, save :: u
    logical, save :: assigned = .false.
    if (.not.assigned) open(newunit=u, file="LOG_LIK_OPTIM"); assigned=.true.
    write (u,*) "log likelihood: noise = ", gp%nu, " theta = ", &
         gp%theta, "log(likelihood) = ", val
  end subroutine output_params

end module m_gp

max_stars_repo_name: JaakkoAhola/GPEmulator

max_issues_repo_name: JaakkoAhola/GPEmulator

max_forks_repo_name: JaakkoAhola/GPEmulator


module WelModule
  !
  use KindModule, only: DP, I4B
  use ConstantsModule, only: DZERO, DEM1, DONE, LENFTYPE
  use BndModule, only: BndType
  use ObsModule, only: DefaultObsIdProcessor
  use SmoothingModule,  only: sQSaturation, sQSaturationDerivative
  use TimeSeriesLinkModule, only: TimeSeriesLinkType, &
                                  GetTimeSeriesLinkFromList
  use BlockParserModule, only: BlockParserType
  !
  implicit none
  !
  private
  public :: wel_create
  !
  character(len=LENFTYPE) :: ftype = 'WEL'
  character(len=16)       :: text  = '             WEL'
  !
  type, extends(BndType) :: WelType
    integer(I4B), pointer :: iflowred => null()
    real(DP), pointer :: flowred => null()
  contains
    procedure :: allocate_scalars => wel_allocate_scalars
    procedure :: bnd_options => wel_options
    procedure :: bnd_cf => wel_cf
    procedure :: bnd_fc => wel_fc
    procedure :: bnd_fn => wel_fn
    procedure :: bnd_da => wel_da
    procedure :: define_listlabel
    ! -- methods for observations
    procedure, public :: bnd_obs_supported => wel_obs_supported
    procedure, public :: bnd_df_obs => wel_df_obs
    ! -- methods for time series
    procedure, public :: bnd_rp_ts => wel_rp_ts
  end type weltype

contains

  subroutine wel_create(packobj, id, ibcnum, inunit, iout, namemodel, pakname)
! ******************************************************************************
! wel_create -- Create a New Well Package
! Subroutine: (1) create new-style package
!             (2) point bndobj to the new package
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- dummy
    class(BndType), pointer :: packobj
    integer(I4B),intent(in) :: id
    integer(I4B),intent(in) :: ibcnum
    integer(I4B),intent(in) :: inunit
    integer(I4B),intent(in) :: iout
    character(len=*), intent(in) :: namemodel
    character(len=*), intent(in) :: pakname
    ! -- local
    type(WelType), pointer :: welobj
! ------------------------------------------------------------------------------
    !
    ! -- allocate the object and assign values to object variables
    allocate(welobj)
    packobj => welobj
    !
    ! -- create name and origin
    call packobj%set_names(ibcnum, namemodel, pakname, ftype)
    packobj%text = text
    !
    ! -- allocate scalars
    call welobj%allocate_scalars()
    !
    ! -- initialize package
    call packobj%pack_initialize()

    packobj%inunit=inunit
    packobj%iout=iout
    packobj%id=id
    packobj%ibcnum = ibcnum
    packobj%ncolbnd=1
    packobj%iscloc=1
    !
    ! -- return
    return
  end subroutine wel_create

  subroutine wel_da(this)
! ******************************************************************************
! wel_da -- deallocate
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- modules
    use MemoryManagerModule, only: mem_deallocate
    ! -- dummy
    class(WelType) :: this
! ------------------------------------------------------------------------------
    !
    ! -- Deallocate parent package
    call this%BndType%bnd_da()
    !
    ! -- scalars
    call mem_deallocate(this%iflowred)
    call mem_deallocate(this%flowred)
    !
    ! -- return
    return
  end subroutine wel_da

  subroutine wel_allocate_scalars(this)
! ******************************************************************************
! allocate_scalars -- allocate scalar members
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    use MemoryManagerModule, only: mem_allocate
    ! -- dummy
    class(WelType) :: this
! ------------------------------------------------------------------------------
    !
    ! -- call standard BndType allocate scalars
    call this%BndType%allocate_scalars()
    !
    ! -- allocate the object and assign values to object variables
    call mem_allocate(this%iflowred, 'IFLOWRED', this%origin)
    call mem_allocate(this%flowred, 'FLOWRED', this%origin)
    !
    ! -- Set values
    this%iflowred = 0
    this%flowred = DZERO
    !
    ! -- return
    return
  end subroutine wel_allocate_scalars

  subroutine wel_options(this, option, found)
! ******************************************************************************
! wel_options -- set options specific to WelType
!
! rch_options overrides BndType%bnd_options
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    use InputOutputModule, only: urword
    ! -- dummy
    class(WelType),   intent(inout) :: this
    character(len=*), intent(inout) :: option
    logical,          intent(inout) :: found
    ! -- local
    real(DP) :: r
    ! -- formats
    character(len=*),parameter :: fmtflowred = &
      "(4x, 'AUTOMATIC FLOW REDUCTION OF WELLS IMPLEMENTED.')"
    character(len=*),parameter :: fmtflowredv = &
      "(4x, 'AUTOMATIC FLOW REDUCTION FRACTION (',g15.7,').')"
! ------------------------------------------------------------------------------
    !
    ! -- Check for 'AUTO_FLOW_REDUCE' and set this%iflowred
    select case (option)
      case('AUTO_FLOW_REDUCE')
        this%iflowred = 1
        r = this%parser%GetDouble()
        if (r <= DZERO) then
          r = DEM1
        else if (r > DONE) then
          r = DONE
        end if
        this%flowred = r
        !
        ! -- Write option and return with found set to true
        if(this%iflowred > 0) &
          write(this%iout, fmtflowred)
          write(this%iout, fmtflowredv) this%flowred
        found = .true.
      case('MOVER')
        this%imover = 1
        write(this%iout, '(4x,A)') 'MOVER OPTION ENABLED'
        found = .true.
      case default
        !
        ! -- No options found
        found = .false.
    end select
    !
    ! -- return
    return
  end subroutine wel_options

  subroutine wel_cf(this)
! ******************************************************************************
! wel_cf -- Formulate the HCOF and RHS terms
! Subroutine: (1) skip in no wells
!             (2) calculate hcof and rhs
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- dummy
    class(WelType) :: this
    ! -- local
    integer(I4B) :: i, node, ict
    real(DP) :: qmult
    real(DP) :: q
    real(DP) :: tp
    real(DP) :: bt
    real(DP) :: thick
! ------------------------------------------------------------------------------
    !
    ! -- Return if no wells
    if(this%nbound == 0) return
    !
    ! -- pakmvrobj cf
    if(this%imover == 1) then
      call this%pakmvrobj%cf()
    endif
    !
    ! -- Calculate hcof and rhs for each well entry
    do i = 1, this%nbound
      node = this%nodelist(i)
      this%hcof(i) = DZERO
      if(this%ibound(node) <= 0) then
        this%rhs(i) = DZERO
        cycle
      end if
      q = this%bound(1,i)
      if (this%iflowred /= 0 .and. q < DZERO) then
        ict = this%icelltype(node)
        if (ict /= 0) then
          tp = this%dis%top(node)
          bt = this%dis%bot(node)
          thick = tp - bt
          tp = bt + this%flowred * thick
          qmult = sQSaturation(tp, bt, this%xnew(node))
          q = q * qmult
        endif
      end if
      this%rhs(i) = -q
    enddo
    !
    return
  end subroutine wel_cf

  subroutine wel_fc(this, rhs, ia, idxglo, amatsln)
! **************************************************************************
! wel_fc -- Copy rhs and hcof into solution rhs and amat
! **************************************************************************
!
!    SPECIFICATIONS:
! --------------------------------------------------------------------------
    ! -- dummy
    class(WelType) :: this
    real(DP), dimension(:), intent(inout) :: rhs
    integer(I4B), dimension(:), intent(in) :: ia
    integer(I4B), dimension(:), intent(in) :: idxglo
    real(DP), dimension(:), intent(inout) :: amatsln
    ! -- local
    integer(I4B) :: i, n, ipos
! --------------------------------------------------------------------------
    !
    ! -- pakmvrobj fc
    if(this%imover == 1) then
      call this%pakmvrobj%fc()
    endif
    !
    ! -- Copy package rhs and hcof into solution rhs and amat
    do i = 1, this%nbound
      n = this%nodelist(i)
      rhs(n) = rhs(n) + this%rhs(i)
      ipos = ia(n)
      amatsln(idxglo(ipos)) = amatsln(idxglo(ipos)) + this%hcof(i)
      !
      ! -- If mover is active and this well is discharging,
      !    store available water (as positive value).
      if(this%imover == 1 .and. this%rhs(i) > DZERO) then
        call this%pakmvrobj%accumulate_qformvr(i, this%rhs(i))
      endif
    enddo
    !
    ! -- return
    return
  end subroutine wel_fc

  subroutine wel_fn(this, rhs, ia, idxglo, amatsln)
! **************************************************************************
! wel_fn -- Fill newton terms
! **************************************************************************
!
!    SPECIFICATIONS:
! --------------------------------------------------------------------------
    implicit none
    ! -- dummy
    class(WelType) :: this
    real(DP), dimension(:), intent(inout) :: rhs
    integer(I4B), dimension(:), intent(in) :: ia
    integer(I4B), dimension(:), intent(in) :: idxglo
    real(DP), dimension(:), intent(inout) :: amatsln
    ! -- local
    integer(I4B) :: i, node, ipos, ict
    real(DP) :: drterm
    real(DP) :: q
    real(DP) :: tp
    real(DP) :: bt
    real(DP) :: thick
! --------------------------------------------------------------------------

    !
    ! -- Copy package rhs and hcof into solution rhs and amat
    do i = 1, this%nbound
      node = this%nodelist(i)
      !
      ! -- test if node is constant or inactive
      if(this%ibound(node) <= 0) then
        cycle
      end if
      !
      ! -- well rate is possibly head dependent
      ict = this%icelltype(node)
      if (this%iflowred /= 0 .and. ict /= 0) then
        ipos = ia(node)
        q = -this%rhs(i)
        if (q < DZERO) then
          ! -- calculate derivative for well
          tp = this%dis%top(node)
          bt = this%dis%bot(node)
          thick = tp - bt
          tp = bt + this%flowred * thick
          drterm = sQSaturationDerivative(tp, bt, this%xnew(node))
          drterm = drterm * this%bound(1,i)
          !--fill amat and rhs with newton-raphson terms
          amatsln(idxglo(ipos)) = amatsln(idxglo(ipos)) + drterm
          rhs(node) = rhs(node) + drterm * this%xnew(node)
        end if
      end if
    end do
    !
    ! -- return
    return
  end subroutine wel_fn

  subroutine define_listlabel(this)
! ******************************************************************************
! define_listlabel -- Define the list heading that is written to iout when
!   PRINT_INPUT option is used.
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    class(WelType), intent(inout) :: this
! ------------------------------------------------------------------------------
    !
    ! -- create the header list label
    this%listlabel = trim(this%filtyp) // ' NO.'
    if(this%dis%ndim == 3) then
      write(this%listlabel, '(a, a7)') trim(this%listlabel), 'LAYER'
      write(this%listlabel, '(a, a7)') trim(this%listlabel), 'ROW'
      write(this%listlabel, '(a, a7)') trim(this%listlabel), 'COL'
    elseif(this%dis%ndim == 2) then
      write(this%listlabel, '(a, a7)') trim(this%listlabel), 'LAYER'
      write(this%listlabel, '(a, a7)') trim(this%listlabel), 'CELL2D'
    else
      write(this%listlabel, '(a, a7)') trim(this%listlabel), 'NODE'
    endif
    write(this%listlabel, '(a, a16)') trim(this%listlabel), 'STRESS RATE'
    if(this%inamedbound == 1) then
      write(this%listlabel, '(a, a16)') trim(this%listlabel), 'BOUNDARY NAME'
    endif
    !
    ! -- return
    return
  end subroutine define_listlabel

  ! -- Procedures related to observations
  logical function wel_obs_supported(this)
  ! ******************************************************************************
  ! wel_obs_supported
  !   -- Return true because WEL package supports observations.
  !   -- Overrides BndType%bnd_obs_supported()
  ! ******************************************************************************
  !
  !    SPECIFICATIONS:
  ! ------------------------------------------------------------------------------
    implicit none
    class(WelType) :: this
  ! ------------------------------------------------------------------------------
    wel_obs_supported = .true.
    return
  end function wel_obs_supported

  subroutine wel_df_obs(this)
  ! ******************************************************************************
  ! wel_df_obs (implements bnd_df_obs)
  !   -- Store observation type supported by WEL package.
  !   -- Overrides BndType%bnd_df_obs
  ! ******************************************************************************
  !
  !    SPECIFICATIONS:
  ! ------------------------------------------------------------------------------
    implicit none
    ! -- dummy
    class(WelType) :: this
    ! -- local
    integer(I4B) :: indx
  ! ------------------------------------------------------------------------------
    call this%obs%StoreObsType('wel', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => DefaultObsIdProcessor
    !
    ! -- Store obs type and assign procedure pointer
    !    for to-mvr observation type.
    call this%obs%StoreObsType('to-mvr', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => DefaultObsIdProcessor
    !
    ! -- return
    return
  end subroutine wel_df_obs

  ! -- Procedure related to time series

  subroutine wel_rp_ts(this)
    ! -- Assign tsLink%Text appropriately for
    !    all time series in use by package.
    !    In the WEL package only the Q variable
    !    can be controlled by time series.
    ! -- dummy
    class(WelType), intent(inout) :: this
    ! -- local
    integer(I4B) :: i, nlinks
    type(TimeSeriesLinkType), pointer :: tslink => null()
    !
    nlinks = this%TsManager%boundtslinks%Count()
    do i=1,nlinks
      tslink => GetTimeSeriesLinkFromList(this%TsManager%boundtslinks, i)
      if (associated(tslink)) then
        if (tslink%JCol==1) then
          tslink%Text = 'Q'
        endif
      endif
    enddo
    !
    return
  end subroutine wel_rp_ts

end module WelModule


text:
!==============================================================================!
  subroutine Comm_Mod_Open_File_Write(fh, file_name)
!------------------------------------------------------------------------------!
!   Open file for writing for sequential runs.                                 !
!------------------------------------------------------------------------------!
  implicit none
!---------------------------------[Arguments]----------------------------------!
  integer   :: fh             ! file handle
  character :: file_name*(*)  ! file name
!==============================================================================!

  fh = 9

  open(fh, FILE=file_name, access='stream', form='unformatted')

  end subroutine

max_issues_repo_name: EdinSmartLab/T-Flows

max_forks_repo_name: EdinSmartLab/T-Flows


text:
A page to document all the fun, interesting, silly, unusual, shocking and downright weird things that make Davis what it is.

Click here if youre actually looking for the comic book store Bizarro World.

 People
See Town Characters
  Preachers with Signs
  Garage People
  Hoopies

 Groups
  Beard Liberation Front
  California Aggie Marching BandUh
  Campus Crusade for Chaos and Confusion
  Cults
  Flashmob Flashmobs a.k.a. Sheeple
  Friends Urging Campus Kindness
  Joke Fraternities
  Milk and Cookies
  Nontheist Church of St. Christopher Discordian
  Students for an Orwellian Society

Signs and Art
  Arroyo Park UFO
  Baby on a Spike Bike
  Bum, Bum, Youve Been Here Before
  Butterfly Pillar
  Cactuar Sign
  Cheshire Cat Discharge
  Concrete Canoe
  Covell Greenbelt Dog Sculptures
  Emma Goldman Institute for Anarchist Studies
  Frog Signs

Stories and News
  $5000 Pita
  Americas Weirdest City
  Bette Davis Incident
  The Blasters The Blasters Band Banned, Pardoned Decades Later

Places
Restaurants
  Delta of Venus
  The Hotdogger

Residences
  Baggins End
  The Bomb Shelter
  Dam Haus
  Gilmore 5
  Human Powered Housing Project

Shopping
  Monitor Graveyard

Buildings and Structures
  The Flying Carousel of the Delta Breeze
  Social Sciences and Humanities Building, aka The Death Star
  Stevenson Bridge
  Toad Tunnel

Streets
  Beau Vine Lane
  The Hobbit Hobbitthemed Streets

Etc.
  1st Street Topiary
  3 Ts
  Big Kids Park
  California National Primate Research Center
  California Raptor Center
  Spooky Places and ensuing Paranormal Activity
  New Zealand Consulate
  Places to Hide

Mysterious or Mystical
  Baháí Grove
  Bike Church
  Buried Treasure
  Causeway Message Rocks
  Tabatha the Titan aka the corpse flower
  Davis Cemetery Labyrinth
  Elmendorf Beast
  The Orgasmatron

Everyday Phenomena
  People being surprised by Bats
  Lots and lots of bicycling
  Seeing a Double Decker Buses doubledecker bus, the Red Wagon Car or any of Davis other unique vehicles drive past
  In winter, the presence of the Freshman Stripe
  Essence of Davis
  Perpetual Puddles
  Collecting Police Department Trading Cards

BicycleRelated
  Bike Jousting
  inevitably getting stuck in a Bike Circle
  People getting Bicycling Under The Influence BUIs (as opposed to DUIs)
  Burning Bike Festival
  Bike Polo
  Bike Sculptures
  The Bike Collective!

Media
  KDVS
  Davis Community Television
  Newsweak

Education
  Experimental College
  Davis Peoples Free School

Events
  April Fools Day, 2006
  Bank Robberies
  Break the Record Night
  Whole Earth Festival
  Picnic Day (especially the Doxie Derby and the Fistulated Cow)
  The annual Butterflies for Beer contest
  Conservative Coming Out Day
  Critical Mass
  Hamster Ball Races
  Pig Day
  International Dadaism Month
  Epic Quad Battle, Epic Quad Battle 2 2, Epic Quad Battle 3 3, and Epic Quad Battle 4 4

UCD Classes
  Classic UC Davis Courses

UCD Traditions
  40 ounce Meal Plan
  Aggie Pack antics at UCD Athletics UCD sporting events
  Cow Tipping Not really a tradition, but more of a rumor of one
  Pajamarino

Aggie BandUh Traditions
  BandUh Sushi Challenge
  The Bossy Cow Cow chant

Bizarro Davis Of The Past
  Dixie the Dinosaur
  Lamargate and Pizzagate
  Paul Moller Paul Mollers Flying Machines
  Tank Rush

Other
  Conehead Dance Extravaganza
  Power Kites
  UFOs UFO citings

Wacky Wiki Pages
  4 or more at UCD
  Buying Souls
  Conversations on bathroom walls
  District de Crozet
  Extreme Food Eating Challenges
  Hats
  Hot Sauce
  Magic Items of Davis

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      do ep=-2,0
      S000000(ep)=2d0*Cv(cc0000+C234,ep)+2d0*m1*Dv(dd0000+N0,ep)
      enddo

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
program phaml_master
use phaml
use phaml_user_mod
implicit none
type(phaml_solution_type) :: soln
probno=1
call phaml_create(soln,nproc=4,update_umod=.true.,triangle_files="period.1")
call phaml_solve_pde(soln,print_grid_when=PHASES,print_grid_who=MASTER, &
                     print_error_when=PHASES,print_error_who=MASTER, &
                     print_error_what=ENERGY_LINF_L2_ERR,max_vert=500, &
                     print_warnings=.false., &
                     refterm=DOUBLE_NVERT,mg_cycles=5)
call phaml_destroy(soln)
end program phaml_master

max_stars_repo_name: qsnake/phaml

max_issues_repo_name: qsnake/phaml

max_forks_repo_name: qsnake/phaml


text:
      REAL FUNCTION SDSDOT(N,SB,SX,INCX,SY,INCY)
c Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
c ALL RIGHTS RESERVED.
c Based on Government Sponsored Research NAS7-03001.
c>> 1996-03-30 SDSDOT Krogh   SNGL => REAL.
C>> 1994-11-11 SDSDOT Krogh   Declared all vars.
C>> 1985-08-02 SDSDOT Lawson  Initial code.
C
C     RETURNS S.P. RESULT WITH DOT PRODUCT ACCUMULATED IN D.P.
C     SDSDOT = SB + SUM FOR I = 0 TO N-1 OF SX(LX+I*INCX)*SY(LY+I*INCY),
C     WHERE LX = 1 IF INCX .GE. 0, ELSE LX = (-INCX)*N, AND LY IS
C     DEFINED IN A SIMILAR WAY USING INCY.
C
      INTEGER N, INCX, INCY, KX, KY, I, NS
      REAL              SX(*),SY(*),SB
      DOUBLE PRECISION DSDOT
C
      DSDOT = DBLE(SB)
      IF(N .LE. 0) GO TO 30
      IF(INCX.EQ.INCY.AND.INCX.GT.0) GO TO 40
      KX = 1
      KY = 1
      IF(INCX.LT.0) KX = 1+(1-N)*INCX
      IF(INCY.LT.0) KY = 1+(1-N)*INCY
          DO 10 I = 1,N
          DSDOT = DSDOT + DBLE(SX(KX))*DBLE(SY(KY))
          KX = KX + INCX
          KY = KY + INCY
   10     CONTINUE
   30 SDSDOT = REAL(DSDOT)
      RETURN
   40 CONTINUE
      NS = N*INCX
          DO 50 I=1,NS,INCX
          DSDOT = DSDOT + DBLE(SX(I))*DBLE(SY(I))
   50     CONTINUE
      SDSDOT = REAL(DSDOT)
      RETURN
      END


!        Generated by TAPENADE     (INRIA, Ecuador team)
!  tapenade 3.x
!
!  Differentiation of vz in reverse (adjoint) mode (with options noISIZE i8):
!   gradient     of useful results: *temp *propunit *pres vz
!   with respect to varying inputs: *temp *propunit *pres
!   Plus diff mem management of: temp:in propunit:in pres:in
!>    @brief calculate velocities at cell faces
!>    @param[in] i grid indices
!>    @param[in] j grid indices
!>    @param[in] k grid indices
!>    @param[in] ismpl local sample index
!>    @return z velocity (m/(Pa s))
SUBROUTINE VZ_AD(i, j, k, ismpl, vz_adv)
  use arrays

  USE ARRAYS_AD

  USE MOD_GENRL
  USE MOD_FLOW
  IMPLICIT NONE
  double precision :: vz_adv
  INTEGER :: ismpl
  INTEGER :: i, j, k
  EXTERNAL GK, RHOF, POR
  EXTERNAL GK_AD0, RHOF_AD
  DOUBLE PRECISION :: f1, f2, dif, GK, vbuoy, POR, &
& RHOF, rhav
  DOUBLE PRECISION :: dif_ad, rhav_ad
  DOUBLE PRECISION :: result1
  DOUBLE PRECISION :: result1_ad
  INTEGER :: arg1
  INTEGER :: branch
  DOUBLE PRECISION :: vz
  IF (k0 .GT. 1 .AND. k .LT. k0) THEN
!          rhav = 0.5D0*(rhof(i,j,k+1,ismpl)+rhof(i,j,k,ismpl))&
!                *(delza(k+1 )-delza(k))
    dif = pres(i, j, k+1, ismpl) - pres(i, j, k, ismpl)
    result1 = POR(i, j, k, ismpl)
    IF (result1 .GT. 1.e-19) THEN
      rhav = RHOF(i, j, k, ismpl)
      CALL PUSHCONTROL1B(0)
    ELSE
      CALL PUSHCONTROL1B(1)
      rhav = 1.29e0
    END IF
    arg1 = k + 1
    result1 = POR(i, j, arg1, ismpl)
    IF (result1 .GT. 1.e-19) THEN
      arg1 = k + 1
      result1 = RHOF(i, j, arg1, ismpl)
      rhav = rhav + result1
      CALL PUSHCONTROL1B(0)
    ELSE
      CALL PUSHCONTROL1B(1)
      rhav = rhav + 1.29e0
    END IF
    rhav = 0.5*rhav*(delza(k+1)-delza(k))
    dif = dif + rhav*grav
! - vbuoy(i,j,k,ismpl)
    result1 = GK(i, j, k, ismpl)
    result1_ad = -(dif*vz_adv)
    dif_ad = -(result1*vz_adv)
    CALL GK_AD0(i, j, k, ismpl, result1_ad)
    rhav_ad = grav*dif_ad
    rhav_ad = 0.5*(delza(k+1)-delza(k))*rhav_ad
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 0) THEN
      result1_ad = rhav_ad
      CALL RHOF_AD(i, j, arg1, ismpl, result1_ad)
    END IF
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 0) CALL RHOF_AD(i, j, k, ismpl, rhav_ad)
    pres_ad(i, j, k+1, ismpl) = pres_ad(i, j, k+1, ismpl) + dif_ad
    pres_ad(i, j, k, ismpl) = pres_ad(i, j, k, ismpl) - dif_ad
  END IF
END SUBROUTINE VZ_AD


text:
c
c
c     =====================================================
      double precision function philim(a,b,meth)
c     =====================================================
      implicit double precision (a-h,o-z)
c
c     # Compute a limiter based on wave strengths a and b.
c     # meth determines what limiter is used.
c     # a is assumed to be nonzero.
c
c     # NOTE: This routine is obsolete.  Instead of using limiter.f,
c     # which calls philim.f for every wave, it is more efficient to
c     # use inlinelimiter.f, which eliminates all these function calls
c     # to philim.  If you wish to change the limiter function and are
c     # using inlinelimiter.f, the formulas must be changed in that routine.
c
      r = b/a
      go to (10,20,30,40,50) meth

c
   10 continue
c     --------
c     # minmod
c     --------
      philim = dmax1(0.d0, dmin1(1.d0, r))
      return
c
   20 continue
c     ----------
c     # superbee
c     ----------
      philim = dmax1(0.d0, dmin1(1.d0, 2.d0*r), dmin1(2.d0, r))
      return
c
   30 continue
c     ----------
c     # van Leer
c     ----------
      philim = (r + dabs(r)) / (1.d0 + dabs(r))
      return
c
   40 continue
c     ------------------------------
c     # monotinized centered
c     ------------------------------
      c = (1.d0 + r)/2.d0
      philim = dmax1(0.d0, dmin1(c, 2.d0, 2.d0*r))
      return
c
   50 continue
c     ------------------------------
c     # Beam-Warming
c     ------------------------------
      philim = r

      return
      end

max_stars_repo_name: ECLAIRWaveS/ForestClaw

max_issues_repo_name: ECLAIRWaveS/ForestClaw

max_forks_repo_name: ECLAIRWaveS/ForestClaw


text:
        subroutine gspiv(b,n,m,eps,rnorms,ipivots,ncols)
        implicit double precision (a-h,o-z)
        save
        dimension rnorms(1),ipivots(1)
        double precision b(n,m),cd
c
c       Set rank to zero just in case.
c
        ncols=0
c
c        . . . initialize the array of pivots
c
        do 1100 i=1,m
        ipivots(i)=i
 1100 continue
c
c       . . . prepare the array of values of norms
c             of columns
c
        done=1
        dtot=0
        do 1400 i=1,m
c
        d=0
        do 1200 j=1,n
        d=d+b(j,i)*(b(j,i))
 1200 continue
        rnorms(i)=sqrt(d)
        dtot=dtot+d
 1400 continue
c
        thresh = eps

        if (dtot .lt. eps**2) then
        ncols = 0
        return
        endif

c$$$        thresh=(dtot)*eps**2
c$$$        thresh=sqrt(thresh)

cccccccccccccccccccccccccccccccccccccccccccc

c
c

c
c       . . . conduct gram-schmidt iterations
c
        do 4000 i=1,min(m,n)
c
c       find the pivot
c
        ipivot=i
        rn=rnorms(i)
c
        do 2200 j=i+1,m
        if(rnorms(j) .le. rn) goto 2200
        rn=rnorms(j)
        ipivot=j
 2200 continue
 2400 continue
c
c       put the column number ipivot in the i-th place
c
        do 2600 j=1,n
        cd=b(j,i)
        b(j,i)=b(j,ipivot)
        b(j,ipivot)=cd
 2600 continue
c
        iijj=ipivots(i)
        ipivots(i)=ipivots(ipivot)
        ipivots(ipivot)=iijj
c
        d=rnorms(ipivot)
        rnorms(ipivot)=rnorms(i)
        rnorms(i)=d
c
c       orthogonalize the i-th column to all preceding ones
c
        if(i .eq. 1) goto 2790
        do 2780 j=1,i-1
c
        call gsrleascap(b(1,i),b(1,j),n,cd)
c
        do 2770 l=1,n
        b(l,i)=b(l,i)-b(l,j)*cd
 2770 continue
 2780 continue
 2790 continue
c
c       normalize the i-th column
c
        call gsrleascap(b(1,i),b(1,i),n,cd)
c
        d=cd
        if(d .lt. thresh**2 ) return
c
        ncols=i
c
        d=done/sqrt(d)
        do 2800 j=1,n
        b(j,i)=b(j,i)*d
 2800 continue
c
        if(i .eq. m) goto 3400
c
c        orthogonalize everything else to it
c
        do 3200 j=i+1,m
c
        if(rnorms(j) .lt. thresh) goto 3200
c
        call gsrleascap(b(1,i),b(1,j),n,cd)
c
        cd=(cd)
c
        rrn=0
        do 3000 l=1,n
        b(l,j)=b(l,j)-b(l,i)*cd
        rrn=rrn+b(l,j)*(b(l,j))
 3000 continue
        rnorms(j)=sqrt(rrn)
 3200 continue
 3400 continue
c
 4000 continue
c
        return
        end
c
c
c
        subroutine gsrleascap(x,y,n,prod)
        implicit double precision (a-h,o-z)
        save
        dimension x(1),y(1)
c
        prod=0
        do 1200 i=1,n
        prod=prod+x(i)*(y(i))
 1200 continue
        return
        end
c
c
c
        subroutine gspiv_move(n,a,b)
        implicit double precision (a-h,o-z)
        dimension a(1),b(1)
        do 1000 j=1,n
        b(j)=a(j)
 1000 continue
        end

max_stars_repo_name: ButterflyLab/ButterflyLab

max_issues_repo_name: ButterflyLab/ButterflyLab

max_forks_repo_name: ButterflyLab/ButterflyLab


PROGRAM square_and_cube_roots
!
!  Purpose:
!    To calculate a table of numbers, square roots, and cube roots
!    using an implied DO loop to output the table.
!
!  Record of revisions:
!      Date       Programmer          Description of change
!      ====       ==========          =====================
!    11/15/06    S. J. Chapman        Original code
!
IMPLICIT NONE

! Data dictionary: declare constants
INTEGER, PARAMETER :: MAX_SIZE = 10      ! Max values in array

! Data dictionary: declare variable types, definitions, & units
INTEGER :: j                             ! Loop index
REAL, DIMENSION(MAX_SIZE) :: value       ! Array of numbers
REAL, DIMENSION(MAX_SIZE) :: square_root ! Array of square roots
REAL, DIMENSION(MAX_SIZE) :: cube_root   ! Array of cube roots

! Calculate the square roots & cube roots of the numbers.
DO j = 1, MAX_SIZE
   value(j) = real(j)
   square_root(j) = sqrt(value(j))
   cube_root(j) = value(j)**(1.0/3.0)
END DO

! Write out each number, its square root, and its cube root.
WRITE (*,100)
100 FORMAT ('0',20X,'Table of Square and Cube Roots',/, &
             4X,'  Number    Square Root  Cube Root', &
             3X,'  Number    Square Root  Cube Root',/, &
             4X,'  ======    ===========  =========', &
             3X,'  ======    ===========  =========')
WRITE (*,110) (value(j), square_root(j), cube_root(j), j = 1, MAX_SIZE)
110 FORMAT (2(4X,F6.0,9X,F6.4,6X,F6.4))

END PROGRAM square_and_cube_roots


text:
 Candidate Statement
> Michael A. Molnars Statement: Ignite Slate <

My name is Michael A. Molnar, but most know me as Molnar. I am a third year Economics, International Relations double major with above a 3.00 GPA. I own my own business doing landscaping and gardening for professors, I am a brother of Sigma Phi Epsilon and I am part of University Covenant Church of Davis. My hometown is Oceanside, CA.

I am going to be a real leader who takes risks. That might mean that I have to go all the way to the chancellor or the UC Regents Regents, but I am willing to go to them if that is what it takes to fight for the students. I will be different from the rest. I am not going to sit around once the campaign is over. I continue to see that many of our senators do not have a passion for student issues, but instead have a desire to have a resume booster. That is why I created the slate Ignite, because I feel that the two main parties are failing UC Davis.

  I am tired of the short Christmas and spring breaks. All the other UC quarter schools have longer breaks. Why do we not start on Mondays at the beginning of the winter and spring quarters? Why do we only get five days off for spring break? I am annoyed with this system. I will go all the way to the chancellor and even the Regents to change this short break policy at UC Davis.
  The Unitrans buses stop running at 7PM on Fridays. Therefore, parking should be free on campus after 7PM on Fridays. This would increase opportunities for clubs to use our campus, which would benefit our community.
  Friends Urging Campus Kindness Platform Other slates promised wireless printing last year, but I still cannot print from my laptop. I want to make this a reality for students. This would reduce lines in the computer labs.
  I would like to see another 24 Hour Reading Room 24hour reading room so students can have seats during finals.
  I want to employ my economic knowledge to the spending of campus funds. Currently the money allocated to ASUCD is being managed unwisely. I desire to change this.

Take it from a leader; I know how to be a senator and I am passionate about initiating change that will better UC Davis for ages to come. On Election Day, vote for a student who cares, VOTE MOLNAR, VOTE IGNITE!

Visit his home page on: Users/MichaelMolnar (page no longer available)

Visit Ignite

Molnar a ASUCD Senator (placed 1st), he works in construction/landscaping around Davis to provide for school. Also, in the future, he may have some interships with State agencies, including the Governator. He is also applying to Officer Navy Supply School now, which is very competitive.

 Please dont write nameruining comments (unless I say something during a senate meeting which is recorded and you disagree with or I do something illegal, which I will not ) because in this day of age, people can lose a job offer from a simple google search (it has happened to UC Davis Students already). I dont think that is what wiki is for, for people to lose future opportunities from google searches online. That is totally unfair, especially when information is not from a valid source. I appreciate the respect.
Senator Molnar

If you feel otherwise, tell me and we can discuss. I am a Senator, I know how to get along with people, especially if they disagree with me because I usually disagree with most of the Senate. However, at the end of the day, we are still friends. That is how politics should be, friends who disagree and can still get along with each other, but are not nasty during elections b/c that sucks.

20060224 16:24:04 nbsp You have to work with the Registrar to change the academic calendar. There is a committee set up to decide calendars three years in advance and there should be better student representation on it. Good luck with everything! Users/AdamBarr

20060225 18:12:51 nbsp Congrats on the win, I am tired of the 2 other parties and I am glad a new voice can be heard. You are now the Dark Horse in student politics at Davis. But great job on the flyers and your campaign. Users/JoshDean

20060426 00:35:43 nbsp Hey, did you find anything out about the junction? Users/ArlenAbraham

20061202 19:30:08 nbsp Which is exactly what I commented on the Ignite page.  The school calender is different each year and there was absolutely no difference made by you or ignite. What Im asking is whether you tried to do anything about it for the future, like you said you would. Users/EdWins ES


text:
module ref_type
  use iso_fortran_env

  implicit none

  type nil
  end type

  type refptr
    integer :: ref_count = 0
    class(*), allocatable :: data
  contains
    final :: refptr_finalize
  end type

  type :: ref
!     private
    type(refptr),pointer :: ptr => null()
    integer :: weak = 0 !1.. this is a weak reference, -1.. make weak references to this
  contains
    procedure :: assign_star
    generic :: assignment(=) => assign_star
    procedure :: set_res
    procedure :: val_real32
    procedure :: val_real64
    procedure :: val_int32
    procedure :: val_int64
    procedure :: val_char
    procedure :: val_logical
    generic :: value => val_real32, val_real64, &
                        val_int32, val_int64, &
                        val_char, val_logical
    procedure :: pointer => ref_pointer
    procedure :: write_formatted => ref_write_formatted
    generic :: write(formatted) => write_formatted
    final :: ref_finalize
  end type

  interface ref
    module procedure ref_init
    module procedure ref_init_1
  end interface

contains

  function ref_init(obj) result(res)
    class(*),intent(in) :: obj
    type(ref) :: res

    allocate(res%ptr)

    select type (obj)
      type is (ref)
        stop "Passing ref not allowed."
    end select
    allocate(res%ptr%data, source = obj)
  end function
  function ref_init_1(obj) result(res)
    class(*),intent(in) :: obj(:)
    type(ref) :: res

    allocate(res%ptr)

    select type (obj)
      type is (ref)
        stop "Passing ref to ref not allowed."
    end select
    allocate(res%ptr%data, source = obj(1))
  end function

  subroutine assign_star(out,in)
    class(*),intent(in) :: in
    class(ref),intent(inout) :: out

    if (associated(out%ptr).and.out%weak<=0) then
      out%ptr%ref_count = out%ptr%ref_count - 1
      if (out%ptr%ref_count<=0) then
        deallocate(out%ptr)
      end if
    end if

    select type (in)
      class is (ref)
        if (associated(in%ptr)) then
          out%ptr => in%ptr
          if (in%weak>=0) then
            out%ptr%ref_count = out%ptr%ref_count + 1
            out%weak = 0
          else
            out%weak = 1
          end if
        else
          out%ptr => null()
        end if
      class default
        if (associated(out%ptr)) nullify(out%ptr)
        allocate(out%ptr)
        allocate(out%ptr%data, source=in)
        out%ptr%ref_count = out%ptr%ref_count + 1
    end select

  end subroutine

  subroutine set_res(self)
    class(ref), intent(inout) :: self
    if (associated(self%ptr).and.self%weak<=0) &
      self%ptr%ref_count = self%ptr%ref_count + 1
  end subroutine

  recursive subroutine ref_finalize(self)
    type(ref),intent(inout) :: self

    if (associated(self%ptr).and.self%weak<=0) then
      self%ptr%ref_count = self%ptr%ref_count - 1
      if (self%ptr%ref_count<=0) then
        deallocate(self%ptr)
      end if
    end if
  end subroutine

  recursive subroutine refptr_finalize(self)
    type(refptr),intent(inout) :: self

    if (allocated(self%data)) deallocate(self%data)
  end subroutine

  subroutine ref_write_formatted(self, unit, iotype, v_list, iostat, iomsg)
    use iso_c_binding
    class(ref), intent(in) :: self
    integer, intent(in) :: unit
    character(*), intent(in) :: iotype
    integer, intent(in) :: v_list(:)
    integer, intent(out) :: iostat
    character(*), intent(inout) :: iomsg

    if (associated(self%ptr)) then
      select type (x => self%ptr%data)
        type is (character(*))
          write(unit,*, iostat=iostat, iomsg=iomsg) x
        type is (integer(int32))
          write(unit,*, iostat=iostat, iomsg=iomsg) x
        type is (integer(int64))
          write(unit,*, iostat=iostat, iomsg=iomsg) x
        type is (real(real32))
          write(unit,*, iostat=iostat, iomsg=iomsg) x
        type is (real(real64))
          write(unit,*, iostat=iostat, iomsg=iomsg) x
        type is (complex(real32))
          write(unit,*, iostat=iostat, iomsg=iomsg) x
        type is (complex(real64))
          write(unit,*, iostat=iostat, iomsg=iomsg) x
        type is (logical)
          write(unit,*, iostat=iostat, iomsg=iomsg) x
        class default
          write(unit,'(a,z0)', iostat=iostat, iomsg=iomsg) "#ref:0x",transfer(c_loc(self%ptr), 1_c_intptr_t)
      end select
    else
      write(unit,*, iostat=iostat, iomsg=iomsg) "null"
    end if
  end subroutine

#define FULL_T SHORT_T
#define RES_T SHORT_T

#define SHORT_T real(real32)
#define NAME_T val_real32
#include "value-inc.f90"
#undef SHORT_T
#undef NAME_T

#define SHORT_T real(real64)
#define NAME_T val_real64
#include "value-inc.f90"
#undef SHORT_T
#undef NAME_T

#define SHORT_T integer(int32)
#define NAME_T val_int32
#include "value-inc.f90"
#undef SHORT_T
#undef NAME_T

#define SHORT_T integer(int64)
#define NAME_T val_int64
#include "value-inc.f90"
#undef SHORT_T
#undef NAME_T

#define SHORT_T logical
#define NAME_T val_logical
#include "value-inc.f90"
#undef SHORT_T
#undef NAME_T

#define SHORT_T character(*)
#undef RES_T
#define RES_T character(:), allocatable
#define NAME_T val_char
#include "value-inc.f90"
#undef SHORT_T
#undef FULL_T
#undef NAME_T

!
  function ref_pointer(self) result(res)
    class(*), pointer             :: res
    class(ref), target, intent(in) :: self

    if (associated(self%ptr)) then
      res => self%ptr%data
    else
      res => null()
    end if
  end function

end module ref_type

#ifdef TEST
module cons_type

  use ref_type

  type :: cons
    type(ref) :: car, cdr
  end type

  interface cons
    module procedure cons_init
  end interface

  contains

  function cons_init(a, b) result(res)
    class(*),intent(in) :: a, b
    type(cons) :: res

    res%car = a
    res%cdr = b
    !res will be finalized, so we must increase the ref counts manually
    call res%car%set_res
    call res%cdr%set_res

#ifdef __GFORTRAN__
    !HACK gfortran still does not finalize function results :(
    res%car%ptr%ref_count = res%car%ptr%ref_count - 1
    res%cdr%ptr%ref_count = res%cdr%ptr%ref_count - 1
#endif
  end function

  recursive subroutine set_car(r, a)
    type(ref),intent(inout) :: r
    class(*),intent(in) :: a
    select type(x=>r%ptr%data)
      type is (cons)
        x%car = a
      class default
        stop "set-car! expects pair"
    end select
  end subroutine

  recursive subroutine set_cdr(r, a)
    type(ref),intent(inout) :: r
    class(*),intent(in) :: a
    select type(x=>r%ptr%data)
      type is (cons)
        x%cdr = a
      class default
        stop "set-cdr! expects pair"
    end select
  end subroutine

end module

program test
  use ref_type
  use cons_type

  implicit none

  call main

contains

  subroutine main
    type(ref) :: a,b,c,d, e

    class(*), pointer :: p

    a = "abc"
    print *, "a:", a

    a = 42
    b = f(f(f(f(a))))
    a = f(f(b))
    b = a
    b = f(a)

    print *, "a:", a, "b:", b

    c = cons(a, b)
    d = c

    p => d%pointer()
    select type(p)
      type is (cons)
        print *, "d: (",p%car,",",p%cdr,")"
    end select

    print *, "d:",d

    e = cons(c, d)

    print *, "reference count of a:", a%ptr%ref_count

  end subroutine

  recursive function f(a) result(res)
    type(ref) a,res
    res = a
#ifdef __GFORTRAN__
    !HACK gfortran still does not finalize function results :(
    res%ptr%ref_count = res%ptr%ref_count - 1
#endif
  end function

end program
#endif

max_stars_repo_name: LadaF/Fortran-RefCount

max_issues_repo_name: LadaF/Fortran-RefCount

max_forks_repo_name: LadaF/Fortran-RefCount


text:
      subroutine intrp(xraw,traw,praw,nraw,xgrid,tgrid,pgrid,ngrid,
     $                            zhypo,dx,branch)
c--------------------------------------------------------------------
c        takes raw input distance,time and ray parameter tabular
c  data that can be very irregularly spaced and grids this data
c  in a safe way.  that is the raw data is both decimated and
c  suplimented by extra rays to produce a table that is relatively
c  evenly spaced.  this table is then interpolated onto a regular
c  grid.  this is safer than interpolating the raw table as experience
c  has proven this can occasionally produce interpolation problems
c  caused by the erratic spacing of the usual raw input.  the
c  routine was designed to also mark the crossover between direct
c  upward ray paths and refracted arrivals as closely as possible.
c  routine is designed to only work on forward branches.  it will
c  surely fail if used on retrograde branches.
c
c  arguments-
c    xraw  - vector of length nraw of raw table of distances.
c             (xraw must be in ascending sequence.
c              i.e. forward branches)
c    traw  - vector of length nraw of raw table of travel times.
c    praw  - vector of length nraw of raw table of ray parameters.
c    nraw  - length of raw arrays.
c    xgrid - array of length ngrid of grid values to be interpolated.
c    tgrid - array of length ngrid to contain interpolated travel times
c    pgrid - array of length ngrid to contain interpolated
c             ray parameter table.
c    ngrid - length of grid arrays.  with refracted branches ngrid is
c             not altered.  with direct (upward) rays ngrid may be
c             increased to better mark the crossover between direct
c             and refracted arrivals.
c   zhypo  - source depth raw table values correspond to.
c   dx     - distance spacing of xgrid values.
c   branch - integer variable containing characters identifying
c             branch type.  (see data statements below)
c
c  external variables-
c    velocity model variable are passed through common/surface/
c         v,d,ns
c
c  written - july 1981
c  language- 1977 ansi standard fortran
c  author  - gary l. pavlis
c--------------------------------------------------------------------
	include 'model.common'
	include 'first.common'
      integer nraw,ngrid
      character branch
      real dx,zhypo
      real xraw(nraw),traw(nraw),praw(nraw)
      real xgrid(ngrid),tgrid(ngrid),pgrid(ngrid)
c--the following parameters determine amount of storage required
c--for work spaces used in this routine.
c--maxdx should be the same as in subroutine gridtp
c--nwsize should be set to some number .gt. maxdx/xminsc  (see below)
c--the actual storage allocated is nwmax = nwsize + 1
c--nwsize is used for loop termination to assure space for at least
c--one extra ray used to mark end of upward branches.
c     parameter(maxdx=50,nwsize=120)
c--rather than hardwired, we now make maxdx = first arrival table size
      parameter(maxdx=nxtab)
      parameter(nwsize=120)
      parameter(nwmax = nwsize + 1)
c--these numbers are used to determine the size of dxmin and dxmax
c--compared to dx (see below)  entries in work array will all
c--normally be .gt. dxmin and .lt.dxmax.
c--larger spacings will occasionally occur if large gaps have to be
c--filled.
      parameter(xminsc=0.5,xmaxsc=0.75)
      real xwork(nwmax),pwork(nwmax),twork(nwmax)
      character uptest
      parameter(uptest='u')
      dxmin = dx*xminsc
      dxmax = dx*xmaxsc
c--the curves are truncated beyond distance xmax
      xmax = dx*float(maxdx)
      nwork = 1
      xmark = xraw(1)
      pmark = praw(1)
      xwork(1) = xraw(1)
      twork(1) = traw(1)
      pwork(1) = praw(1)
      do 200 i=2,nraw
            deltax = abs(xraw(i) - xmark)
            if(deltax.gt.dxmin) then
                  if(deltax.le.dxmax) then
c--block for regions to be decimated or left alone
                        if(nwork.ge.nwsize) go to 250
                        nwork = nwork + 1
                        xmark = xraw(i)
                        pmark = praw(i)
                        xwork(nwork) = xraw(i)
                        twork(nwork) = traw(i)
                        pwork(nwork) = praw(i)
                  else
c--block for jumps that need to be filled.
                        nsplit = nint(deltax/dxmax) + 1
                        if((nsplit+nwork-1).gt.nwsize)
     $                          nsplit = nwsize - nwork + 1
                        if(nsplit.le.2) then
                              if(nwork.ge.nwsize) go to 250
                              nwork = nwork + 1
                              xmark = xraw(i)
                              pmark = praw(i)
                              xwork(nwork) = xraw(i)
                              twork(nwork) = traw(i)
                              pwork(nwork) = praw(i)
                        else
                              call split(pwork(nwork),xwork(nwork),
     $                            pmark,xmark,praw(i),xraw(i),nsplit)
                              if(branch.eq.uptest) then
                                    do 100 j=nwork+1,nwork+nsplit-2
                                          call txup(pwork(j),zhypo,vel,
     $                                  dep,nvel,twork(j),xwork(j),ibot)
  100                               continue
                              else
                                    do 150 j=nwork+1,nwork+nsplit-2
                                          call txref(pwork(j),zhypo,vel,
     $                                  dep,nvel,twork(j),xwork(j),ibot)
  150                               continue
                              endif
                              nwork = nwork + nsplit - 1
                              xmark = xraw(i)
                              pmark = praw(i)
                              xwork(nwork) = xraw(i)
                              twork(nwork) = traw(i)
                              pwork(nwork) = praw(i)
                        endif
                  endif
                  if(xmark.gt.xmax) go to 250
            endif
  200 continue
c--alternate loop exit
  250 continue
c--this section makes sure the final raw grid point is kept verbatim
c--to avoid disasterous truncation of regular grid
      if(xmark.le.xmax) then
            if(xmark.ne.xraw(nraw)) then
                  nwork = nwork + 1
                  xwork(nwork) = xraw(nraw)
                  twork(nwork) = traw(nraw)
                  pwork(nwork) = praw(nraw)
                  deltax = xwork(nwork) - xwork(nwork-1)
                  if(deltax.lt.0.0) nwork = nwork - 1
                  if(deltax.gt.dxmax) then
                        print *,'//////warning//////'
                        print *,'possible error in curves at depth=',
     $                           zhypo
                        print *,'gridded points between ',
     $                       xwork(nwork-1),' and ',xwork(nwork),
     $                       ' are subject to interpolation error'
                 endif
            endif
      endif
c--the following actually does the interpolation.
      call linear(xwork,twork,nwork,xgrid,tgrid,ngrid)
      call linear(xwork,pwork,nwork,xgrid,pgrid,ngrid)
      return
      end

c $Id$


SUBROUTINE p104(input_file,output_file)
!--------------------------------------------------------------------------
! Program 10.4: Eigenvalue analysis of an elastic solid in plane strain
!               using 4-node rectangular quadrilaterals.
!               Arnoldi's Method. Lumped mass. Element by element.
!               Mesh numbered in x- or y-direction.
!--------------------------------------------------------------------------
 USE main
 USE geom
 IMPLICIT NONE
 CHARACTER(len=60),INTENT(IN) :: input_file
 CHARACTER(len=60),INTENT(OUT) :: output_file
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 INTEGER::i,ido,iel,ierr,info,iparam(11),ipntr(11),ishfts,iters,j,k,      &
   logfil,lworkl,maxitr,model,msaitr,msapps,msaupd,msaup2,mseigt,mseupd,  &
   msgets,ndigit,nconv,ncv,ndim=2,nev,ndof,nels,neq,nip=4,nlen,nn,nod,    &
   nodof=2,nprops=3,np_types,nr,nst=3,nxe,nye
 REAL(iwp)::det,one=1.0_iwp,sigma,tol,zero=0.0_iwp
 CHARACTER(len=15)::argv,bmat,element='quadrilateral',which
 LOGICAL::rvec
!--------------------------- dynamic arrays--------------------------------
 REAL(iwp),ALLOCATABLE::bee(:,:),coord(:,:),d(:,:),dee(:,:),der(:,:),     &
   deriv(:,:),diag(:),ecm(:,:),fun(:),g_coord(:,:),jac(:,:),km(:,:),      &
   mm(:,:),pmul(:),points(:,:),prop(:,:),resid(:),udiag(:),utemp(:),      &
   v(:,:),vdiag(:),weights(:),workd(:),workl(:),x_coords(:),y_coords(:)
 INTEGER,ALLOCATABLE::etype(:),g(:),g_g (:,:),g_num(:,:),nf(:,:),num(:)
 LOGICAL,ALLOCATABLE::select(:)
! include 'debug.h'
!----------------------input and initialisation----------------------------

  OPEN(10,FILE=input_file)
  OPEN(11,FILE=output_file)
  READ(10,*)nxe,nye,nod,np_types
  CALL mesh_size(element,nod,nels,nn,nxe,nye);ndof=nod*nodof
  ALLOCATE(nf(nodof,nn),points(nip,ndim),dee(nst,nst),g_coord(ndim,nn),   &
    coord(nod,ndim),fun(nod),jac(ndim,ndim), weights(nip),utemp(ndof),    &
    g_num(nod,nels),der(ndim,nod),deriv(ndim,nod),bee(nst,ndof),num(nod), &
    km(ndof,ndof),g(ndof),g_g(ndof,nels),mm(ndof,ndof),ecm(ndof,ndof),    &
    pmul(ndof),x_coords(nxe+1),y_coords(nye+1),prop(nprops,np_types),     &
    etype(nels))
  READ(10,*)prop
  etype=1
  IF(np_types>1)READ(10,*)etype
  READ(10,*)x_coords,y_coords
  nf=1
  READ(10,*)nr,(k,nf(:,k),i=1,nr)
  CALL formnf(nf)
  neq=MAXVAL(nf)
  READ(10,*)nev,ncv,bmat,which,tol,maxitr
  lworkl=ncv*(ncv+8)
  ndigit=-3
  logfil=6
  msgets=0
  msaitr=0
  msapps=0
  msaupd=2
  msaup2=0
  mseigt=0
  mseupd=0
! msaupd=1
  msaup2=0
  mseigt=0
  mseupd=0
  info=0
  ido=0
  ishfts=1
  model=1
  ipntr=0
  iparam=0
  iparam(1)=ishfts
  iparam(3)=maxitr
  iparam(7)=model
!-----------------------loop the elements to find global array sizes-------
  elements_1: DO iel=1,nels
    CALL geom_rect(element,iel,x_coords,y_coords,coord,num,'y')
    CALL num_to_g (num,nf,g)
    g_num(:,iel)=num
    g_coord(:,num)=transpose(coord)
    g_g(:,iel)=g
  END DO elements_1
  CALL mesh(g_coord,g_num,argv,nlen,12)
  WRITE(11,'(A,I5,A/)')" There are",neq," equations"
  ALLOCATE(v(neq,ncv),diag(neq),udiag(neq),vdiag(neq),workd(3*neq),       &
    resid(neq),workl(lworkl),d(ncv,2),select(ncv))
  diag=zero
  v=zero
  udiag=zero
  vdiag=zero
  workd=zero
  resid=zero
  CALL sample(element,points,weights)
!--------------- element stiffness integration and assembly----------------
 elements_2: DO iel=1,nels
   num=g_num(:,iel)
   coord=transpose(g_coord(:,num ))
   g=g_g(:,iel)
   km=zero
   mm=zero
   CALL deemat(dee,prop(1,etype(iel)),prop(2,etype(iel)))
   integrating_pts_1: DO i=1,nip
     CALL shape_fun(fun,points,i)
     CALL shape_der(der,points,i)
     jac=MATMUL(der,coord)
     det=determinant(jac)
     CALL invert(jac)
     deriv=MATMUL(jac,der)
     CALL beemat(bee,deriv)
     km=km+MATMUL(MATMUL(TRANSPOSE(bee),dee),bee)*det*weights(i)
     CALL ecmat(ecm,fun,ndof,nodof)
     mm=mm+ecm*det*weights(i)**prop(3,etype(iel))
   END DO integrating_pts_1
   DO i=1,ndof
     IF(g(i)/=0)diag(g(i))=diag(g(i))+sum(mm(i,:))
   END DO
 END DO elements_2
!------------------------------find eigenvalues----------------------------
 iters=0
 diag=one/SQRT(diag)    ! diag holds l**(-1/2)
 DO
   iters=iters+1
   CALL dsaupd(ido,bmat,neq,which,nev,tol,resid,ncv,v,neq,iparam,ipntr,   &
     workd,workl,lworkl,info)
   IF(ido/=-1.AND.ido/=1)EXIT
   udiag=zero
   vdiag=workd(ipntr(1):ipntr(1)+neq-1)*diag
   elements_3: DO iel=1,nels
   g=g_g(:,iel)
     DO i=1,ndof
       IF(g(i)==0)pmul(i)=zero
       IF(g(i)/=0)pmul(i)=vdiag(g(i))
     END DO
     utemp=matmul(km,pmul)
     DO i=1,ndof
       IF(g(i)/=0)udiag(g(i))=udiag(g(i))+utemp(i)
     END DO
   END DO elements_3
   udiag=udiag*diag
   workd(ipntr(2):ipntr(2)+neq-1)=udiag
  END DO
  WRITE(11,'(A,I5,A/)')" Convergence after ",iters," iterations"
!-----Either we have convergence or there is an error----------------------
  IF(info<0)THEN
    WRITE(11,*)"Fatal error in ssaupd "
  ELSE
    rvec=.TRUE.
    CALL dseupd(rvec,'All',select,d,v,neq,sigma,bmat,neq,which,nev,tol,   &
      resid,ncv,v,neq,iparam,ipntr,workd,workl,lworkl,ierr)
!------------------------- write out the spectrum -------------------------
    IF(ierr/=0 )THEN
      WRITE(11,*)"Fatal error in sseupd"
      WRITE(11,*)"Should compute residuals"
    END IF
    WRITE(11,'(A)')"The eigenvalues are:"
    WRITE(11,'(6E12.4)')d(1:nev,1)
    DO i=1,nev
      udiag(:)=v(1:neq,i)
      udiag=udiag*diag
      WRITE(11,'(A,I3,A)')" Eigenvector number",i," is:"
      WRITE(11,'(6E12.4)')udiag(1:)/MAXVAL(ABS(udiag(1:)))
      IF(i==1)CALL dismsh((/zero,udiag/),nf,0.1_iwp,g_coord,g_num,argv,nlen,13)
    END DO
  END IF
END SUBROUTINE p104


text:
!
! Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
! See https://llvm.org/LICENSE.txt for license information.
! SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
!
      module trailz_elemental
         implicit none
         contains
            elemental integer  function do_trailz(a) result(b)
               integer, intent (in) ::  a
               b=trailz(a)
            end function
      end module

      program test
      use trailz_elemental
       implicit none
       integer, parameter :: num = 5
       integer results(num)
       integer , parameter :: expect(num) =(/2,6,6,0,0/)
       integer , parameter :: arr(num)=(/-108,-64,64,-1,1/)

       results=do_trailz (arr)

       call check(results, expect, num)
      end program

max_stars_repo_name: abrahamtovarmob/flang

max_issues_repo_name: abrahamtovarmob/flang

max_forks_repo_name: abrahamtovarmob/flang


text:
        subroutine dgemm(a,b,c,n)
          double precision a
          dimension a(*)
          double precision b
          dimension b(*)
          double precision c
          dimension c(*)
          integer n

          integer i,j,k
          do k=0,n-1,1
             do j=0,n-1,1
                do i=0,n-1,1
                   c(j*n+i) = c(j*n+i)+a(k*n+i)*b(j*n+k)
                enddo
             enddo
          enddo
        end

max_stars_repo_name: ouankou/rose

max_issues_repo_name: WildeGeist/rose

max_forks_repo_name: WildeGeist/rose


text:
!!  Copyright (C)  Stichting Deltares, 2012-2016.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License version 3,
!!  as published by the Free Software Foundation.
!!
!!  This program is distributed in the hope that it will be useful,
!!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!!  GNU General Public License for more details.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program. If not, see <http://www.gnu.org/licenses/>.
!!
!!  contact: delft3d.support@deltares.nl
!!  Stichting Deltares
!!  P.O. Box 177
!!  2600 MH Delft, The Netherlands
!!
!!  All indications and logos of, and references to registered trademarks
!!  of Stichting Deltares remain the property of Stichting Deltares. All
!!  rights reserved.

      subroutine barx(w, wbar)
!           4-69
!           for each compartment k compute wbar(k) as the sum of w(j),
!        where j is each substance in the compartment
      include  'char1.inc'
      include  'char5.inc'
      dimension w(1), wbar(1)
      do 20 k=1,ncomp
         kta = kl(k)
         ktb = kl(k+1) - 1
         wbar(k) = slacks
         if(bslack) wbar(k) = slack(k)
         do 10 j=kta,ktb
            wbar(k) = wbar(k) + w(j)
   10    continue
   20 continue
      return
      end

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
c=======================================================================
c
c     simulateur et fonction de trace de courbe pour RNLS
c
c     polynomial function f = an*x**n + .......... + a1*x + a0
c
c=======================================================================
c
c     SIMEXT
c
c     polynome: a(n)*x**n + .......... + a(1)*x + a(0)
c
c     Subroutine provided by the user to compute the value of
c     the functions vector(nfct) and the jacobian matrix(nfct*mx)
c     in function of xsol vector(mx)
c
c=======================================================================
      SUBROUTINE SIMEXT ( nfct, mxopt, xopt, spadim, grid,
     &                    liudat, iusdat, ldudat, dusdat,
     &                    funct, jacob, info )
c=======================================================================
c
c     INPUT :
c            nfct   : number of functions                        integer
c            mxopt  : size of the vector solution                integer
c            xopt   : solution  vector(mxopt)                     double
c                     in the order a(0),a(1),........,a(n-1),a(n)
c            spadim : space dimension                            integer
c            grid   : grid     vector(nfct*spadim)                double
c            liudat : size of integers data provided by the user integer
c            iusdat : integers user data provided by the user
c                     vector(liudat)                             integer
c            ldudat : size of doubles data provided by the user  integer
c            dusdat : doubles user data provided by the user
c                     vector(ldudat)                              double
c
c     OUTPUT :
c            funct  : functions   vector(nfct)                    double
c            jacob  : jacobian    matrix(nfct*mx)                 double
c            info   : = 0 successful exit                        integer
c                     = else : it's your problem
c
c-----------------------------------------------------------------------
c
      implicit none
c
      integer nfct,mxopt,spadim,liudat,ldudat,info
      integer iusdat(*)
      double precision xopt(*),grid(*),dusdat(*),funct(*),jacob(nfct,*)
c
      integer i,j
      double precision z,r
c
c-----------------------------------------------------------------------
c
      info = 0
c
       DO i=1,nfct
         z = grid(i)
         CALL POLYN ( mxopt-1, xopt,  z, funct(i) )
         r = 1.
         DO j=1,mxopt
            jacob(i,j) = r
            r = r * z
         ENDDO
      ENDDO
      RETURN
      END
c
c=======================================================================
c
c     SIMTRA
c
c     Computing RNLS files for graphic output for polynomial function
c
c=======================================================================
c
c     INPUT :
c            npc    : number of functions                        integer
c            mparam : number of parameters                       integer
c            pref   : reference parameters  vector(mparam)        double
c            prnls  : RNLS parameters  vector(mparam)             double
c            pnls   : NLS parameters  matrix(mparam,nmes)         double
c            npobs  : number of observations                     integer
c            pobs   : points of observation  vector(npobs)        double
c            obs    : mean observations  vector(npobs)            double
c            nmes   : number of measures                         integer
c            mesure : measures matrix(nmes,npobs)                 double
c
c-----------------------------------------------------------------------
c
      subroutine SIMTRA ( npc, mparam, pref, prnls, pnls, npobs, pobs,
     &                    obs, nmes, mesure )
c
      implicit none
c
      integer npc, mparam, npobs, nmes
      double precision pnls(mparam,*), mesure(nmes,*)
      double precision pref(*), prnls(*), pobs(*), obs(*)
c
c.......................................................................
c
      integer i, j, imes
      double precision pasc, z, debord, debut, fin
      double precision zc(npc), fc(npc,nmes), fref(npc), frnls(npc)
c
c-----------------------------------------------------------------------
c
      open(20,status='unknown',file='dim')
      open(21,status='unknown',file='obs')
      open(22,status='unknown',file='mod')
c
      write(20,*) npobs, nmes, npc
c
      debord = ( pobs(npobs) - pobs(1) ) / 100.
      debut =  pobs(1) - debord
      fin =  pobs(npobs) + debord
      pasc = ( fin - debut ) / ( npc - 1 )
c
      z = debut
      do i=1,npc
         zc(i) = z
         z = z + pasc
      end do
c
c     ecriture des points de mesure
c
      do i=1,npobs
         write(21,*)pobs(i),obs(i),(mesure(j,i),j=1,nmes)
      end do
c
c     calcul des nmes courbes
c
      do i=1,npc
         do imes=1,nmes
            call POLYN ( mparam-1, pnls(1,imes), zc(i), fc(i,imes) )
         end do
         call POLYN ( mparam-1, pref, zc(i), fref(i) )
         call POLYN ( mparam-1, prnls, zc(i), frnls(i) )
      end do
c
c     ecriture des courbes
c
      do i=1,npc
         write(22,*)zc(i),(fc(i,imes),imes=1,nmes),frnls(i),fref(i)
      end do
c
      close(20)
      close(21)
      close(22)
c
      return
      end

max_stars_repo_name: alpgodev/numx

max_issues_repo_name: alpgodev/numx

max_forks_repo_name: alpgodev/numx


!*==dlarzt.f90  processed by SPAG 7.51RB at 20:08 on  3 Mar 2022
!> \brief \b DLARZT forms the triangular factor T of a block reflector H = I - vtvH.
!
!  =========== DOCUMENTATION ===========
!
! Online html documentation available at
!            http://www.netlib.org/lapack/explore-html/
!
!> \htmlonly
!> Download DLARZT + dependencies
!>
!> [TGZ]
!>
!> [ZIP]
!>
!> [TXT]
!> \endhtmlonly
!
!  Definition:
!  ===========
!
!       SUBROUTINE DLARZT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
!
!       .. Scalar Arguments ..
!       CHARACTER          DIRECT, STOREV
!       INTEGER            K, LDT, LDV, N
!       ..
!       .. Array Arguments ..
!       DOUBLE PRECISION   T( LDT, * ), TAU( * ), V( LDV, * )
!       ..
!
!
!> \par Purpose:
!  =============
!>
!> \verbatim
!>
!> DLARZT forms the triangular factor T of a real block reflector
!> H of order > n, which is defined as a product of k elementary
!> reflectors.
!>
!> If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular;
!>
!> If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular.
!>
!> If STOREV = 'C', the vector which defines the elementary reflector
!> H(i) is stored in the i-th column of the array V, and
!>
!>    H  =  I - V * T * V**T
!>
!> If STOREV = 'R', the vector which defines the elementary reflector
!> H(i) is stored in the i-th row of the array V, and
!>
!>    H  =  I - V**T * T * V
!>
!> Currently, only STOREV = 'R' and DIRECT = 'B' are supported.
!> \endverbatim
!
!  Arguments:
!  ==========
!
!> \param[in] DIRECT
!> \verbatim
!>          DIRECT is CHARACTER*1
!>          Specifies the order in which the elementary reflectors are
!>          multiplied to form the block reflector:
!>          = 'F': H = H(1) H(2) . . . H(k) (Forward, not supported yet)
!>          = 'B': H = H(k) . . . H(2) H(1) (Backward)
!> \endverbatim
!>
!> \param[in] STOREV
!> \verbatim
!>          STOREV is CHARACTER*1
!>          Specifies how the vectors which define the elementary
!>          reflectors are stored (see also Further Details):
!>          = 'C': columnwise                        (not supported yet)
!>          = 'R': rowwise
!> \endverbatim
!>
!> \param[in] N
!> \verbatim
!>          N is INTEGER
!>          The order of the block reflector H. N >= 0.
!> \endverbatim
!>
!> \param[in] K
!> \verbatim
!>          K is INTEGER
!>          The order of the triangular factor T (= the number of
!>          elementary reflectors). K >= 1.
!> \endverbatim
!>
!> \param[in,out] V
!> \verbatim
!>          V is DOUBLE PRECISION array, dimension
!>                               (LDV,K) if STOREV = 'C'
!>                               (LDV,N) if STOREV = 'R'
!>          The matrix V. See further details.
!> \endverbatim
!>
!> \param[in] LDV
!> \verbatim
!>          LDV is INTEGER
!>          The leading dimension of the array V.
!>          If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K.
!> \endverbatim
!>
!> \param[in] TAU
!> \verbatim
!>          TAU is DOUBLE PRECISION array, dimension (K)
!>          TAU(i) must contain the scalar factor of the elementary
!>          reflector H(i).
!> \endverbatim
!>
!> \param[out] T
!> \verbatim
!>          T is DOUBLE PRECISION array, dimension (LDT,K)
!>          The k by k triangular factor T of the block reflector.
!>          If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is
!>          lower triangular. The rest of the array is not used.
!> \endverbatim
!>
!> \param[in] LDT
!> \verbatim
!>          LDT is INTEGER
!>          The leading dimension of the array T. LDT >= K.
!> \endverbatim
!
!  Authors:
!  ========
!
!> \author Univ. of Tennessee
!> \author Univ. of California Berkeley
!> \author Univ. of Colorado Denver
!> \author NAG Ltd.
!
!> \date December 2016
!
!> \ingroup doubleOTHERcomputational
!
!> \par Contributors:
!  ==================
!>
!>    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
!
!> \par Further Details:
!  =====================
!>
!> \verbatim
!>
!>  The shape of the matrix V and the storage of the vectors which define
!>  the H(i) is best illustrated by the following example with n = 5 and
!>  k = 3. The elements equal to 1 are not stored; the corresponding
!>  array elements are modified but restored on exit. The rest of the
!>  array is not used.
!>
!>  DIRECT = 'F' and STOREV = 'C':         DIRECT = 'F' and STOREV = 'R':
!>
!>                                              ______V_____
!>         ( v1 v2 v3 )                        /            \
!>         ( v1 v2 v3 )                      ( v1 v1 v1 v1 v1 . . . . 1 )
!>     V = ( v1 v2 v3 )                      ( v2 v2 v2 v2 v2 . . . 1   )
!>         ( v1 v2 v3 )                      ( v3 v3 v3 v3 v3 . . 1     )
!>         ( v1 v2 v3 )
!>            .  .  .
!>            .  .  .
!>            1  .  .
!>               1  .
!>                  1
!>
!>  DIRECT = 'B' and STOREV = 'C':         DIRECT = 'B' and STOREV = 'R':
!>
!>                                                        ______V_____
!>            1                                          /            \
!>            .  1                           ( 1 . . . . v1 v1 v1 v1 v1 )
!>            .  .  1                        ( . 1 . . . v2 v2 v2 v2 v2 )
!>            .  .  .                        ( . . 1 . . v3 v3 v3 v3 v3 )
!>            .  .  .
!>         ( v1 v2 v3 )
!>         ( v1 v2 v3 )
!>     V = ( v1 v2 v3 )
!>         ( v1 v2 v3 )
!>         ( v1 v2 v3 )
!> \endverbatim
!>
!  =====================================================================
      SUBROUTINE DLARZT(Direct,Storev,N,K,V,Ldv,Tau,T,Ldt)
      IMPLICIT NONE
!*--DLARZT189
!
!  -- LAPACK computational routine (version 3.7.0) --
!  -- LAPACK is a software package provided by Univ. of Tennessee,    --
!  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
!     December 2016
!
!     .. Scalar Arguments ..
      CHARACTER Direct , Storev
      INTEGER K , Ldt , Ldv , N
!     ..
!     .. Array Arguments ..
      DOUBLE PRECISION T(Ldt,*) , Tau(*) , V(Ldv,*)
!     ..
!
!  =====================================================================
!
!     .. Parameters ..
      DOUBLE PRECISION ZERO
      PARAMETER (ZERO=0.0D+0)
!     ..
!     .. Local Scalars ..
      INTEGER i , info , j
!     ..
!     .. External Subroutines ..
      EXTERNAL DGEMV , DTRMV , XERBLA
!     ..
!     .. External Functions ..
      LOGICAL LSAME
      EXTERNAL LSAME
!     ..
!     .. Executable Statements ..
!
!     Check for currently supported options
!
      info = 0
      IF ( .NOT.LSAME(Direct,'B') ) THEN
         info = -1
      ELSEIF ( .NOT.LSAME(Storev,'R') ) THEN
         info = -2
      ENDIF
      IF ( info/=0 ) THEN
         CALL XERBLA('DLARZT',-info)
         RETURN
      ENDIF
!
      DO i = K , 1 , -1
         IF ( Tau(i)==ZERO ) THEN
!
!           H(i)  =  I
!
            DO j = i , K
               T(j,i) = ZERO
            ENDDO
         ELSE
!
!           general case
!
            IF ( i<K ) THEN
!
!              T(i+1:k,i) = - tau(i) * V(i+1:k,1:n) * V(i,1:n)**T
!
               CALL DGEMV('No transpose',K-i,N,-Tau(i),V(i+1,1),Ldv,    &
     &                    V(i,1),Ldv,ZERO,T(i+1,i),1)
!
!              T(i+1:k,i) = T(i+1:k,i+1:k) * T(i+1:k,i)
!
               CALL DTRMV('Lower','No transpose','Non-unit',K-i,        &
     &                    T(i+1,i+1),Ldt,T(i+1,i),1)
            ENDIF
            T(i,i) = Tau(i)
         ENDIF
      ENDDO
!
!     End of DLARZT
!
      END SUBROUTINE DLARZT


text:
       SUBROUTINE Pl
C
C  CALCULATES PERCEIVED LEVEL OF SPECTRA BASED ON
C  1/3 OCTAVE BANDS FROM 1 Hz TO 10 kHz.
C
c                         Dave McCurdy; modified for use with a
c                         single spectrum by BMS, Feb 12 1988
c                         Extended the frequency range, Feb 19 1988
C      (Can't be directly extended up in frequency - 12.5 KHz doesn't work)
c ref1 = S.S.Stevens, `Perceived Level of Noise by Mark VII and Decibels (E)'
c                  JASA 51,2(2),575-593, Feb 1972
c ref2 = K.S.Pearsons and R.L.Bennett, `Handbook of Noises Ratings'
c                  NASA CR-2376

c This routine passes the 1/3 octave band data in a common block, which could
c be changed so the data are passed in as parameters to fit in a more modular
c program.

c L=level (43 1/3 octave bands)
c SPL=calculated loudness level (calculated from L for bands 1 to 40)
      REAL*4 L(43),SPL
      COMMON /CRUX/ SPL,L

      REAL*4 Noy(43),Noys

      REAL*4 Sm(93),Fs(93)
      REAL*4 Smax

C  READ IN DATA FOR TABLE Sm

      DATA Sm /.181,.196,.212,.230,.248,.269,.290,.314,.339,.367,
     1    .396,.428,.463,.500,.540,.583,.630,.680,.735,.794,
     1    .857,.926,1.00,1.08,1.17,1.26,1.36,1.47,1.59,1.72,
     1    1.85,2.00,2.16,2.33,2.52,2.72,2.74,3.18,3.43,3.70,
     1    4.00,4.32,4.67,5.04,5.44,5.88,6.35,6.86,7.41,8.00,
     1    8.64,9.33,10.1,10.9,11.8,12.7,13.7,14.8,16.0,17.3,
     1    18.7,20.2,21.8,23.5,25.4,27.4,29.6,32.0,34.6,37.3,
     1    40.3,43.5,47.0,50.8,54.9,59.3,64.0,69.1,74.7,80.6,
     1    87.1,94.1,102.,110.,119.,128.,138.,149.,161.,174.,
     1   188.,203.,219./

C READ IN DATA FOR TABLE Fs
      DATA Fs /.100,.122,.140,.158,.174,.187,.200,.212,.222,.232,
     1    .241,.250,.259,.267,.274,.281,.287,.293,.298,.303,
     1    .308,.312,.316,.319,.320,.322,.322,.320,.319,.317,
     1    .314,.311,.308,.304,.300,.296,.292,.288,.284,.279,
     1    .275,.270,.266,.262,.258,.253,.248,.244,.240,.235,
     1    .230,.226,.222,.217,.212,.208,.204,.200,.197,.195,
     1    .194,.193,.192,.191,.190,.190,.190,.190,.190,.190,
     1    .191,.191,.192,.193,.194,.195,.197,.199,.201,.203,
     1    .205,.208,.210,.212,.215,.217,.219,.221,.223,.224,
     1    .225,.226,.227/
C
c Initialize variables
         DO 1411 J=1,43
            Noy(J) = 0.0
 1411    CONTINUE
         Noys=0
         Smax=0

c Look at every 1/3 octave band (Nu)
         DO 1458 Nu=1,40
c Find the 1/3 o.b. level in dB
            P=L(Nu)
            IF (Nu-19 .LT. 0) THEN
c If the band center frequency is less than 80 Hz
               Pp=160.-19./Nu*(160.-P)
               Ppp=115.-26./19.*(115.-Pp)
               IF (Ppp .GT. 76.) Ppp=Pp-10.5
               IF (Ppp .GT. 121.) Ppp=160.-26./Nu*(160.-P)
            ELSE IF (Nu-26 .LT. 0) THEN
c If the band c.f. is greater than or equal to 80 and less than 400 Hz
               Pp=P
               Ppp=115.-26./Nu*(115.-Pp)
               IF (Ppp .GT. 76.) Ppp=Pp-(26.-Nu)*1.5
               IF (Ppp .GT. 121.) Ppp=160.-26./Nu*(160.-P)
            ELSE
c If the band c.f. is greater than or equal to 400 Hz
               Ppp=P
            ENDIF
C
            IF (Nu-31 .LE. 0) THEN
c If the band c.f. is less than or equal to 1250 Hz
               Pppp=40.
            ELSE IF (Nu-32 .LE. 0) THEN
c If the band c.f. is greater than 1250 and less than or equal to 1600 Hz
c (i.e. for 1600 Hz band)
               Pppp=38.
            ELSE IF (Nu-33 .LE. 0) THEN
c If the band c.f. is greater than 1600 and less than or equal to 2000 Hz
c (i.e. for 2000 Hz band)
               Pppp=36.
            ELSE IF (Nu-34 .LE. 0) THEN
c If the band c.f. is greater than 2000 and less than or equal to 2500 Hz
c (i.e. for 2500 Hz band)
               Pppp=34.
            ELSE IF (Nu-39 .LE. 0) THEN
c If the band c.f. is greater than 2500 and less than or equal to 8000 Hz
               Pppp=32.
            ELSE
c If the band c.f. is greater than 8000 Hz
               Pppp=36.
            ENDIF
c
            Noy(Nu)=2.**((Ppp-Pppp)/9.)
            Noys=Noys+Noy(Nu)
            Smax=AMAX1(Smax,Noy(Nu))
c	    print *,Nu,Noy(Nu)
 1458    CONTINUE

         Fofsm=0.
         S=Smax
         IF (S .GE. Sm(1)) THEN
            DO 1471 I=2,93
               IF (S-Sm(I) .LT. 0.) THEN
                  Fofsm=Fs(I-1)+(Fs(I)-Fs(I-1))*(S-Sm(I-1))
     1                        /(Sm(I)-Sm(I-1))
                  GOTO 1473
               ELSE IF (S-Sm(I) .EQ. 0.) THEN
                  Fofsm=Fs(I)
                  GOTO 1473
               ENDIF
 1471       CONTINUE
            Fofsm=Fs(93)
         ENDIF

 1473    St=Smax+Fofsm*(Noys-Smax)
         Spl=32+29.89735285*ALOG10(St)

C	 print *,Smax,Noys,Fofsm,Spl

         RETURN
         END

max_stars_repo_name: mandalin/DissertationPostProcessing

max_issues_repo_name: mandalin/DissertationPostProcessing

max_forks_repo_name: mandalin/DissertationPostProcessing