text:
C######################################################################|
      module calculate_hot_tail_population
C----------------------------------------------------------------------|
C     This module is used to calculate the hot-tail population.
C
C     Functions and subroutine:
C       - hot_tail_density
C       - set_Te
C       - set_ne
C       - EVDF_modified_Maxwellian
C       - tau
C       - v_c
C       - v_th
C       - print_hot_tail_info
C
C     Function arguments:
C       - time: time after onset of TQ (s)
C       - t_de: electron temperature decay time scale (s)
C       - ne:   electron density at time `time` (m**-3)
C       - ne_i: electron density at onset of TQ (m**-3)
C       - ne_f: electron density at end of TQ (m**-3)
C       - T:    electron temperature at time `time` (eV)
C       - T_i:  electron temperature at onset of TQ (eV)
C       - T_f:  electron temperature at end of TQ (eV)
C       - Epar: parallel electric field at time `time` (V/m)
C       - v:    electron velocity (m/s)
C
C----------------------------------------------------------------------|
      use double
      use physical_constants, only :
     >  e, epsilon_0, m_e, pi
      use file_io, only :
     >  store_real_arr
      use Coulomb_logarithms, only :
     >  ln_Lambda_0
      use collision_frequencies, only :
     >  nu_ee

      implicit none

      private

      public ::
     >  hot_tail_density,
     >  EVDF_modified_Maxwellian,
     >  tau,
     >  v_c,
     >  v_th

        ! ----- Parameters --------------------------------------------|
      integer, parameter ::
     >  n_eval_start= 200,              ! start # of integrand evaluatns
     >  n_eval_max  = 100000000         ! max # of integrand evaluations

      real(kind=dp) ::
     >  EVDF_min_val= 1.e-100_dp,       ! to calc. up. integration bound
     >  iter_acc    = 1.0e-8_dp         ! target accuracy of iterations

      contains

C======================================================================|
      real(kind=dp) function hot_tail_density(time, t_dec, Epar, T, T_i,
     >  T_f, ne, ne_i, ne_f, verbose, store_result, un) result(n_ht)
C----------------------------------------------------------------------|
C     Evaluates the density of hot-tail runaway electron according to
C     the model by H M Smith and E Verwichte. Phys. Plasmas 15, 072502
C     (2008), equations (9), (17), (18).
C
C     Input:
C       time:       time after onset of TQ (s)
C       t_dec:      time scale of exponential Te decay (s)
C       Epar:       parallel electric field at time `time` (V/m)
C       **T:        electron temperature at time `time` (eV)
C       **T_i:      electron temperature at onset of TQ (eV)
C       **T_f:      electron temperature at end of TQ (eV)
C       **ne:       electron density at time `time` (m**-3)
C       **ne_i:     electron density at onset of TQ (m**-3)
C       **ne_f:     electron density at end of TQ (m**-3)
C       **verbose:  write intermediate data to screen
C       **store_re: store results in text file
C       **un:       unit number of opened file to store results in
C     ** = optional
C     Note, that at least two values of T, T_i, or T_f are required!
C     Note, that at least one value of ne, ne_i or ne_f is required!
C
C     Output:
C       n_ht:      density of hot-tail runaway electrons (m**-3)
C
C----------------------------------------------------------------------|
        ! ----- Function arguments ------------------------------------|
      integer, intent(in), optional ::
     >  un

      logical, intent(in), optional ::
     >  verbose, store_result

      real(kind=dp), intent(in) ::
     >  time, t_dec, Epar
      real(kind=dp), intent(in), optional ::
     >  T, T_i, T_f, ne, ne_i, ne_f

        ! ----- Local variables ---------------------------------------|
      character(len=20) ::
     >  frmt

      integer ::
     >  i, my_un, n_eval

      logical ::
     >  verbos, store_results

      real(kind=dp) ::
     >  dv, int_bound_upp, int_bound_low, n_e, n_ef, n_ei, n_ht_old,
     >  tau_t, T_e, T_ef, T_ei, v, v_th0, weight

C----------------------------------------------------------------------|
C     Initialize hot-tail density
C----------------------------------------------------------------------|
      n_ht = 0._dp

C----------------------------------------------------------------------|
C     Check input
C----------------------------------------------------------------------|
        ! ----- Electron temperature ----------------------------------|
      T_e  = 0._dp;     if (present(T   )) T_e  = T
      T_ei = 0._dp;     if (present(T_i )) T_ei = T_i
      T_ef = 0._dp;     if (present(T_f )) T_ef = T_f

      call set_Te(time, t_dec, T_e, T_ei, T_ef)

        ! ----- Electron density --------------------------------------|
      n_e  = 0._dp;     if (present(ne  )) n_e  = ne
      n_ei = 0._dp;     if (present(ne_i)) n_ei = ne_i
      n_ef = 0._dp;     if (present(ne_f)) n_ef = ne_f

      call set_ne(time, t_dec, n_e, n_ei, n_ef)

        ! ----- Flags -------------------------------------------------|
      verbos = .false.; if (present(verbose)) verbos = verbose
      store_results = .false.
        if (present(store_result)) store_results = store_result
      my_un = 0;       if (present(un)) my_un = un

C----------------------------------------------------------------------|
C     Determine integration bounds, normalized to thermal velocity
C----------------------------------------------------------------------|
      tau_t = tau(time, t_dec, n_ei, n_ef, T_i)
      v_th0 = v_th(T_i)

      int_bound_low = v_c(n_e, T, Epar)/v_th0
      int_bound_upp =
     >  (sqrt((-log(EVDF_min_val))**3) - 3*tau_t)**(1._dp/3._dp)

      if (int_bound_upp .le. int_bound_low) return

C----------------------------------------------------------------------|
C     Print calculation parameters
C----------------------------------------------------------------------|
      call print_hot_tail_info(verbose=verbos,
     >  args= ( /time, t_dec, Epar, T_i, T, n_ei, n_e, n_ef, v_th0,
     >      v_c(n_e, T, Epar), tau_t, int_bound_low, int_bound_upp /))

C----------------------------------------------------------------------|
C     Perform numerical integration using Kepler's rule
C----------------------------------------------------------------------|
        ! Set densities such that loop is executed
      n_ht     = 1._dp
      n_ht_old = 2._dp

        ! Initial number of integrand evaluations
        ! Has to be even to have uneven number of points
      n_eval = n_eval_start

        ! ----- Iteratively integrate ---------------------------------|
      do while (abs((n_ht - n_ht_old)/n_ht_old) .gt. iter_acc
     >      .and. n_eval .le. n_eval_max)

            ! Store previous result
        n_ht_old = n_ht

            ! Integration step size
        dv = 1._dp/n_eval*(int_bound_upp - int_bound_low)

            ! Integrand is zero for v = int_bound_low
        n_ht = 0._dp
        do i=1,n_eval
                ! Weights for integrand
            if (i.lt.n_eval) then
                weight = 2._dp**(1+mod(i,2))
            else
                weight = 1._dp
            endif

            v = i*dv + int_bound_low
            n_ht = n_ht + weight*(v**2 - int_bound_low**2)
     >          * exp(-(v**3 + 3*tau_t)**(2._dp/3._dp))
        enddo
            ! Apply pre-factors
        n_ht = 4*n_ei/sqrt(pi)*dv/3*n_ht

            ! Write iteration step to screen
        call print_hot_tail_info(verbose=verbos,
     >      args = (/ 1+log(1._dp*n_eval/n_eval_start)/log(2._dp),
     >          n_ht, abs(n_ht/n_ht_old-1.), real(n_eval, kind=dp) /))

            ! Increase number of integrand evaluations
        n_eval = 2*n_eval
      enddo

C----------------------------------------------------------------------|
C     Store results if requested
C----------------------------------------------------------------------|
      if (store_results .and. my_un.gt.0) then
        call store_real_arr(my_un, (/ time, n_ht, ne, T, Epar,
     >      v_c(ne, T, Epar)/v_th0, tau_t /))
      endif

      return
      end function hot_tail_density
C======================================================================|

C======================================================================|
      subroutine set_Te(t, t_dec, y, y_i, y_f)
C----------------------------------------------------------------------|
C     Ensures the electron temperature has the correct format to
C     calculate the hot-tail density.
C----------------------------------------------------------------------|
        ! ----- Subroutine arguments ----------------------------------|
      real(kind=dp), intent(in) ::
     >  t, t_dec

      real(kind=dp), intent(inout) ::
     >  y, y_i, y_f

C----------------------------------------------------------------------|
C     Check different cases
C----------------------------------------------------------------------|
        ! Everything provided: return
      if (y.gt.0._dp .and. y_i.gt.0._dp) then
        return

        ! Initial and final value provided: calculate value at `t`
      else if (y.le.0._dp .and. y_i.gt.0._dp .and. y_f.gt.0._dp) then
        y = (y_i - y_f)*exp(-t/t_dec) + y_f

        ! Value at `t` and final value provided: calculate initial value
      else if (y.gt.0._dp .and. y_i.le.0._dp .and. y_f.gt.0._dp) then
        y_i = (y - y_f*(1._dp - exp(-t/t_dec)))*exp(t/t_dec)

        ! Otherwise, cannot set required quantities `y` and `y_i` and
        ! calculation cannot be performed
      else
        write(*,'(a)') 'ERROR: Insufficient electron temperature ' //
     >      'data provided to perform calculation of hot-tail ' //
     >      'density. STOP calculation.'
        stop

      endif

      return
      end subroutine set_Te
C======================================================================|

C======================================================================|
      subroutine set_ne(t, t_dec, y, y_i, y_f)
C----------------------------------------------------------------------|
C     Ensures the electron density has the correct format to calculate
C     the hot-tail density.
C----------------------------------------------------------------------|
        ! ----- Parameters --------------------------------------------|
      real(kind=dp), parameter ::
     >  y_def = 1.e19_dp

        ! ----- Subroutine arguments ----------------------------------|
      real(kind=dp), intent(in) ::
     >  t, t_dec

      real(kind=dp), intent(inout) ::
     >  y, y_i, y_f

C----------------------------------------------------------------------|
C     Check different cases
C----------------------------------------------------------------------|
        ! Everything provided: return
      if (y.gt.0._dp .and. y_i.gt.0._dp .and. y_f.gt.0._dp) then
        return

        ! Initial and final value provided: calculate value at `t`
      else if (y.le.0._dp .and. y_i.gt.0._dp .and. y_f.gt.0._dp) then
        y = (y_i - y_f)*exp(-t/t_dec) + y_f

        ! Initial and value at `t` provided: calculate final value
      else if (y.gt.0._dp .and. y_i.gt.0._dp .and. y_f.le.0._dp) then
        y_f = (y - y_i*exp(-t/t_dec))/(1._dp - exp(-t/t_dec))

        ! Value at `t` and final value provided: calculate initial value
      else if (y.gt.0._dp .and. y_i.le.0._dp .and. y_f.gt.0._dp) then
        y_i = (y - y_f*(1._dp - exp(-t/t_dec)))*exp(t/t_dec)

        ! Only one value provided: Set for all values
      else if (y.gt.0._dp .or. y_i.gt.0._dp .or. y_f.gt.0._dp) then
        y   = max(y, y_i, y_f)
        y_i = y
        y_f = y

        ! Nothing provided: WARNING
      else if (y.le.0._dp .and. y_i.le.0._dp .and. y_f.le.0._dp) then
        write(*,'(a,es12.4,a)') 'WARNING: No values for electron ' //
     >      'density provided for calculation of hot-tail density! ' //
     >      'Will use default density of ne=', y_def,' m**-3.'
        y   = y_def
        y_i = y_def
        y_f = y_def

      endif

      return
      end subroutine set_ne
C======================================================================|

C======================================================================|
      real(kind=dp) function EVDF_modified_Maxwellian(v, ne, v_th, tau)
     >  result(distr_func)
C----------------------------------------------------------------------|
C     Calculates the value of the Maxwellian electron velocity
C     distribution function at velocity `v` for electron density `ne`,
C     thermal velocity `v_th` and parameter `tau`, where the argument
C     of the exponential has been modified. Reduces to a standard
C     Maxwellian for tau=0. For details, see H.M. Smith and E.
C     Verwichte. Phys. Plasmas 15, 072502 (2008), especially eq. (9).
C
C     Input
C       v:          electron velocity where distribution function is
C                   to be evaluated (m/s)
C       ne:         electron density (m**-3)
C       v_th:       thermal electron velocity (m/s)
C       tau:        parameter of temporal delay
C
C     Output:
C       distr_func: magnitude of the electron distribution function
C
C----------------------------------------------------------------------|
          ! ----- Function arguments ----------------------------------|
      real(kind=dp), intent(in) ::
     >  v, ne, v_th, tau

C----------------------------------------------------------------------|
      distr_func = ne*exp(-((v/v_th)**3 + 3*tau)**(2._dp/3._dp))
     >  /(sqrt(pi)*v_th)**3

      return
      end function EVDF_modified_Maxwellian
C======================================================================|

C======================================================================|
      real(kind=dp) function tau(time, t_char, ne_i, ne_f, T_i)
C----------------------------------------------------------------------|
C     Calculates the parameter tau of temporal delay. For details, see
C     H.M. Smith and E. Verwichte. Phys. Plasmas 15, 072502 (2008).
C
C     Input:
C       time:       time since onset of TQ (s)
C       ne:         electron density at time `time` (m**-3)
C       t_char:     characteristic time scale of TQ (s)
C       ne_i:       electron density at onset of TQ (m**-3)
C       ne_f:       electron density at end of TQ (m**-3)
C       T_i:        electron temperature at onset of TQ (eV)
C
C     Output:
C       tau:        parameter of temporal delay
C
C----------------------------------------------------------------------|
          ! ----- Function arguments ----------------------------------|
      real(kind=dp), intent(in) ::
     >  time, t_char, ne_i, ne_f, T_i

C----------------------------------------------------------------------|
      if (time .lt. 2*t_char) then
        tau = time**2/4/t_char
      else if (time .ge. 2*t_char) then
        tau = time - t_char
      endif

      tau = tau*nu_ee(ne_i, T_i)*ne_f/ne_i

      return
      end function tau
C======================================================================|

C======================================================================|
      real(kind=dp) function v_c(ne, T, Epar)
C----------------------------------------------------------------------|
C     Calculates the critical velocity for electron runaway. See e.g.
C     H.M. Smith and Epar. Verwichte. Phys. Plasmas 15, 072502 (2008).
C
C     Input:
C       ne:         electron density (m**-3)
C       T:          electron temperature (eV)
C       Epar:       parallel electric field (V/m)
C
C     Output:
C       v_th:       thermal electron velocity (m/s)
C
C----------------------------------------------------------------------|
          ! ----- Function arguments ----------------------------------|
      real(kind=dp), intent(in) ::
     >  ne, T, Epar

C----------------------------------------------------------------------|
      v_c = sqrt(e**3/(4*pi*epsilon_0**2*m_e)
     >  *ne*ln_Lambda_0(ne, T)/Epar)

      return
      end function v_c
C======================================================================|

C======================================================================|
      real(kind=dp) function v_th(T)
C----------------------------------------------------------------------|
C     Calculates the thermal electron velocity v_th = sqrt(2*T/m_e).
C
C     Input:
C       T:          electron temperature (eV)
C
C     Output:
C       v_th:       thermal electron velocity (m/s)
C
C----------------------------------------------------------------------|
          ! ----- Function arguments ----------------------------------|
      real(kind=dp), intent(in) ::
     >  T
C----------------------------------------------------------------------|
      v_th = sqrt(2*T*e/m_e)

      return
      end function v_th
C======================================================================|

C======================================================================|
      subroutine print_hot_tail_info(args, verbose)
C----------------------------------------------------------------------|
C     Prints information on hot-tail calculations to the screen. If the
C     array `args` contains 13 elements, the calculation paramters are
C     printed. If `args` contains 4 elements, the density at the current
C     iteration step is printed.
C----------------------------------------------------------------------|
        ! ----- Parameters --------------------------------------------|
      character(len=10), parameter ::
     >  frmt = '(a,3e12.4)'

        ! ----- Subroutine arguments ----------------------------------|
      logical, intent(in), optional ::
     >  verbose

      real(kind=dp), dimension(:), intent(in) ::
     >  args

        ! ----- Local variables ---------------------------------------|
      logical ::
     >  verbos

C----------------------------------------------------------------------|
C     Determine print flag `verbose`
C----------------------------------------------------------------------|
      verbos=.true.; if(present(verbose)) verbos=verbose
      if (.not. verbos) return

C----------------------------------------------------------------------|
C     Print overview of calculation parameters
C----------------------------------------------------------------------|
      if ((size(args) .eq. 13) .and. verbos) then
        write(*,*) ''
        write(*,'(2a)') 'Perform numerical integration for the ',
     >  'following set of parameters:'
        write(*,frmt) 'Time (s):                  ', args(1)
        write(*,frmt) 'Decay time (s):            ', args(2)
        write(*,frmt) 'Epar (V/m):                ', args(3)
        write(*,frmt) 'Te (ini, `t`) (eV):        ', args(4), args(5)
        write(*,frmt) 'ne (ini, `t`, fin) (m**-3):', args(6), args(7),
     >      args(8)
        write(*,frmt) 'Velocities (vth, vc) (m/s):', args(9), args(10)
        write(*,frmt) 'Temporal parameter tau:    ', args(11)
        write(*,frmt) 'Integration bounds (v_th): ', args(12), args(13)
        write(*,frmt) ''
        write(*,'(a)') ' Step, Density (m**-3), Rel. Difference, # eval'
        write(*,'(a)') ' -----------------------------'

C----------------------------------------------------------------------|
C     Print current iteration step
C----------------------------------------------------------------------|
      else if ((size(args) .eq. 4) .and. verbos) then
        write(*,'(a,i0.3,2e17.8,i9)') ' #', int(args(1)), args(2),
     >      args(3), int(args(4))

      endif

      return
      end subroutine print_hot_tail_info
C======================================================================|

      end module calculate_hot_tail_population
C######################################################################|

max_stars_repo_name: o-linder/runawayelectrongeneration

max_issues_repo_name: o-linder/runawayelectrongeneration

max_forks_repo_name: o-linder/runawayelectrongeneration


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
!

! C1547 (R1526) MODULE shall appear only in the function-stmt or subroutine-stmt
! of a module subprogram or of a nonabstract interface body that is declared in
! the scoping unit of a module or submodule.

MODULE mod_test
  IMPLICIT NONE

ABSTRACT INTERFACE

    PURE INTEGER MODULE FUNCTION f1(i) !{error "PGF90-S-1055-MODULE prefix cannot be inside an abstract interface"}
      INTEGER, INTENT(IN) :: i
    END FUNCTION f1
END INTERFACE

END module mod_test

max_stars_repo_name: abrahamtovarmob/flang

max_issues_repo_name: abrahamtovarmob/flang

max_forks_repo_name: abrahamtovarmob/flang


        !COMPILER-GENERATED INTERFACE MODULE: Fri Mar  6 15:35:15 2020
        ! This source file is for reference only and may not completely
        ! represent the generated interface used by the compiler.
        MODULE LEAF_LITTER__genmod
          INTERFACE
            SUBROUTINE LEAF_LITTER(DUM_CVEG,DUM_PHOTO,DUM_RESPVEG,      &
     &DUM_EPSV,DUM_LEAF)
              REAL(KIND=8), INTENT(IN) :: DUM_CVEG
              REAL(KIND=8), INTENT(IN) :: DUM_PHOTO
              REAL(KIND=8), INTENT(IN) :: DUM_RESPVEG
              REAL(KIND=8), INTENT(OUT) :: DUM_EPSV
              REAL(KIND=8), INTENT(OUT) :: DUM_LEAF
            END SUBROUTINE LEAF_LITTER
          END INTERFACE
        END MODULE LEAF_LITTER__genmod


text:
module constants_lib
    use constants_def

    implicit none
contains

    subroutine constants_init(mesa_dir_init,ierr)
        use const_lib
        use exceptions_lib
        character(len=*), intent(in) :: mesa_dir_init
        integer, intent(out) :: ierr
        type(assertion) :: mesa_const_init = assertion(scope='constants_init')
        ierr = 0
        call const_init(mesa_dir_init,ierr)
        call mesa_const_init% assert(ierr==0)
        call initialize_constants
    end subroutine constants_init

end module constants_lib

ext: f

max_stars_repo_name: LBJ-Wade/dStar

max_issues_repo_name: LBJ-Wade/dStar

max_forks_repo_name: LBJ-Wade/dStar


text:
A political party is a group of people banded together for certain political ideals.  How they are banded together varies and what they support varies (some parties ask for strict support of their candidates, others do not).  Davis, like most towns, has a mix of people who choose to affiliate with a political party.

Davis Political Party Breakdown

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


c    Subroutines for the program sdf_lister
c
c     Contents:
c             logical function open_output_ave
c
c             subroutine assign_fluxes_ave
c
c             logical function write_output_ave
c
c
c****************************************************
c
       logical function open_output_ave(unit_ave)
c
c****************************************************
c
c     Open output data file and write the header for average matrix rates data
c     returns .true. if successful, .false. otherwise
c
c     12/4/95 by J. Dwyer
c          modified 2/28/96 by J. Dwyer added headerlines output to write_header
c                         4/8/96 by J. Dwyer changed headerlines to 40
c                         4/9/96 by J. Dwyer added year to output
c                         4/24/96 changed energy labels
c                          1/20/97 by J. Dwyer removed anisotropy data
c                          2/13/97 by J. Dwyer   added new KP list
c                          2/14/97 by J. Dwyer   removed junk, added write_kptitles
c
      include 'sdf_include.inc'   ! include type declarations
c
      integer   unit_ave   ! file unit of output file
      character*64     filename       ! complete name of output file
      character*4       file_prefix   ! prefix to be attached to output file name
      character*400   kpstring     !  string containing KP data names
      integer   headerlines
      integer   stringlength   ! length of KP string
c
      headerlines = 40
c
      file_prefix = ' AVE'
      write(filename, '(a4,a60)') file_prefix, outputfile
      open(unit=unit_ave,name=outputdir//filename,
     1     status='NEW',DISP='KEEP',RECL=32766, err=1000) ! open output data file
      print *, 'Opened output file: ', filename
c
      call write_header(unit_ave,file_prefix, headerlines)
c
      call write_kptitles(kpstring,stringlength)  !  get string containing KP data names
c
      write(unit_ave,200) kpstring(1:stringlength)
200  	format(' Time,',
     1         ' year,',
     1         ' day of the year,',
     1         ' delta_t/2,',
     1         ' saturation flag,',
     1         ' VSE, sigma,',
     1         ' VSEBAR, sigma,',
     1         ' START, sigma,',
     1         ' STOP, sigma, D, sigma,',
     1         ' STATE, sigma,',
     1         ' H E1, sigma,',
     1         ' H E2, sigma,',
     1         ' H E3, sigma,',
     1         ' H E4, sigma,',
     1         ' H E5, sigma,',
     1         ' H E6, sigma,',
     1         ' H E7, sigma,',
     1         ' H E8, sigma,',
     1         ' H E9, sigma,',
     1         ' H E10, sigma,',
     1         ' He E1, sigma,',
     1         ' He E2, sigma,',
     1         ' He E3, sigma,',
     1         ' He E4, sigma,',
     1         ' He E5, sigma,',
     1         ' He E6, sigma,',
     1         ' He E7, sigma,',
     1         ' CNO E1, sigma,',
     1         ' CNO E2, sigma,',
     1         ' CNO E3, sigma,',
     1         ' CNO E4, sigma,',
     1         ' CNO E5, sigma,',
     1         ' CNO E6, sigma,',
     1         ' CNO E7, sigma,',
     1         ' NeS E1, sigma,',
     1         ' NeS E2, sigma,',
     1         ' NeS E3, sigma,',
     1         ' NeS E4, sigma,',
     1         ' NeS E5, sigma,',
     1         ' NeS E6, sigma,',
     1         ' NeS E7, sigma,',
     1         ' Fe E1, sigma,',
     1         ' Fe E2, sigma,',
     1         ' Fe E3, sigma,',
     1         ' Fe E4, sigma,',
     1         ' Fe E5, sigma,',
     1         ' Fe E6, sigma,',
c     1         ' Junk-1, sigma, Junk-2, sigma, Junk-3, sigma,',
     1           a<stringlength>)
c
      open_output_ave = .true.
      return  !  succussful
1000   open_output_ave = .false.  ! on error go to here
      print *, 'Error opening output averege matrix rates file:',
     1     filename
      return
      end  ! end open_output_ave
c
c
c
c
c**********************************************************
c
      subroutine assign_fluxes_ave(vse_ave,
     1              sigma_vse_ave, vse_ave_flag,
     1          flux_ave, sigma_flux_ave, flux_ave_flag)
c
c**********************************************************
c
c    calculates sector averaged fluxes and sigmas for 1 or both telescopes
c
c     10/12/95 by J. Dwyer
c                 modified 4/22/96 corrected addition of two tel. when effic_flag = -1
c
      include 'sdf_include.inc'   ! include type declarations
c
      real    delta_tvse  ! delta time over which disc rate is calculated
      real    vse1, vse2
      real    vse_ave, sigma_vse_ave
      real    flux_tel(41,2), sigma_tel(41,2) ! arrays containing the fluxes and sigmas for all the elements
      real    flux_ave(41), sigma_flux_ave(41)   ! arrays containing the fluxes and sigmas for all the elements
      integer   flux_flag_tel(41,2)   ! .true. when flux is measured and is not saturated, .false. otherwise
      integer   flux_ave_flag(41)   ! .true. when flux is measured and is not saturated, .false. otherwise
      integer   vse1_flag, vse2_flag
      integer   vse_ave_flag
      integer   Ntel
c
      vse1 = 0.0
      vse2 = 0.0
      vse_ave = 0.0
      sigma_vse_ave = 0.0
      do 5 i=1,41
        flux_ave_flag(i) = 0
        flux_ave(i) = 0.0
        sigma_flux_ave(i) = 0.0
        do 5 j=1,2
           flux_flag_tel(i,j) = 0
           flux_tel(i,j) = 0.0
5         sigma_tel(i,j) = 0.0
c
      delta_tvse = time_interval   ! delta time over which disc rate is calculated
      if ((stateab_flag.eq.0).or.(stateab_flag.eq.1))   ! if selecting only state A or B
     1      delta_tvse = time_interval/2.0
c
c
      do 10 k=1,8
         if (vse_flag(1,k).eq.-1) then   ! don't include saturated counts
            vse1_flag = -1
            vse1 = 0.0
            goto 15
         else
            if (vse_flag(1,k).eq.1) then  ! include only values measured in this mf
               vse1_flag = 1
               vse1=vse1+vse(1,k)  ! do sum over sectors
            end if
         end if
10     continue
15    continue
      do 20 k=1,8
         if (vse_flag(2,k).eq.-1) then   ! don't include saturated counts
            vse2_flag = -1
            vse2 = 0.0
            goto 25
         else
            if (vse_flag(2,k).eq.1) then  ! include only values measured in this mf
               vse2_flag = 1
               vse2=vse2+vse(2,k)  ! do sum over sectors
            end if
         end if
20     continue
25     continue
c
       if (tel_flag.eq.1) then   ! telescope 1 only
         if (vse1_flag.eq.1) then
             vse_ave_flag = 1
             vse_ave = vse1/delta_tvse
             sigma_vse_ave = sqrt(1.0*vse1)/delta_tvse
         else
             vse_ave_flag = -1
         end if
      end if
c
       if (tel_flag.eq.2) then   ! telescope 2 only
        if (vse2_flag.eq.1) then
             vse_ave_flag = 1
             vse_ave = vse2/delta_tvse
             sigma_vse_ave = sqrt(1.0*vse2)/delta_tvse
         else
             vse_ave_flag = -1
         end if
      end if
c
       if (tel_flag.eq.-1) then   ! telescope 1 and 2
          if ((vse1_flag.eq.1).and.(vse2_flag.eq.1)) then
             vse_ave_flag = 1
             vse_ave = (vse1+vse2)/delta_tvse
             sigma_vse_ave = sqrt(1.0*(vse1+vse2))/delta_tvse
         else
              vse_ave_flag = -1
         end if
      end if
c
c
       do 60 i=1,41
         do 50 j=1,2
            do 40 k=1,8
                if (mdata_flag(i,j,k).eq.-1) then   ! don't include saturated counts
                   flux_flag_tel(i,j) = -1
                   flux_tel(i,j) = 0.0
                   goto 50
                end if
                if (mdata_flag(i,j,k).eq.-2) then   ! don't include saturated counts
                   flux_flag_tel(i,j) = -2
                   flux_tel(i,j) = 0.0
                end if
                if (mdata_flag(i,j,k).eq.1) then  ! include only values measured in this mf
                    flux_flag_tel(i,j) = 1
                    flux_tel(i,j)=flux_tel(i,j)+mdata(i,j,k)  ! do sum over sectors
                 end if
40        continue
45          if (flux_flag_tel(i,j).eq.1) then
              sigma_tel(i,j)=sqrt(flux_tel(i,j))/exposure(i,j)   ! calculate error in flux
              flux_tel(i,j) = flux_tel(i,j)/exposure(i,j)  ! calculate flux from counts and exposure
            end if
50        continue
60        continue
c
c           calculate total fluxes (telescope 1, 2 or 1+2)
c        **************************************
      if (tel_flag.eq.1) then   ! telescope 1 only
         do 100 i=1,41
            flux_ave_flag(i) = flux_flag_tel(i,1)
            flux_ave(i) = flux_tel(i,1)
100      sigma_flux_ave(i) = sigma_tel(i,1)
      end if
       if (tel_flag.eq.2) then   ! telescope 2 only
         do 200 i=1,40
            flux_ave_flag(i) = flux_flag_tel(i,2)
            flux_ave(i) = flux_tel(i,2)
200      sigma_flux_ave(i) = sigma_tel(i,2)
      end if
      if (tel_flag.eq.-1) then   ! both telescopes
         do 300 i=1,41
            Ntel = 0
            flux_ave_flag(i) = 0
            flux_ave(i) =0.0
            if (flux_flag_tel(i,1).eq.1) then
                Ntel = 1
                flux_ave_flag(i) = 1
                flux_ave(i) = flux_tel(i,1)
                sigma_flux_ave(i) =
     1           sigma_tel(i,1)*sigma_tel(i,1)
            end if
            if (flux_flag_tel(i,2).eq.1) then
                Ntel = Ntel+1
                flux_ave_flag(i) = 1
                flux_ave(i) = flux_ave(i)+flux_tel(i,2)
                sigma_flux_ave(i) = sigma_flux_ave(i)+
     1           sigma_tel(i,2)*sigma_tel(i,2)
            end if
            if (Ntel.gt.0) then
              flux_ave(i) = flux_ave(i)/(1.0*Ntel)
              sigma_flux_ave(i) =
     1            sqrt(sigma_flux_ave(i))/(1.0*Ntel)
            end if
            if ((flux_flag_tel(i,1).eq.-2).and.
     1        (flux_flag_tel(i,2).eq.-2))  flux_ave_flag(i) = -2
            if ((flux_flag_tel(i,1).eq.-1).or.
     1        (flux_flag_tel(i,2).eq.-1))  flux_ave_flag(i) = -1
300     continue
      end if
c
      return
      end  ! end assign_flux_ave
c
c
c
c
c****************************************************
c
       logical function write_output_ave(unit_ave, nrow_ave,
     1             vse_ave, sigma_vse_ave, vse_ave_flag,
     1             flux_ave, sigma_flux_ave,flux_ave_flag, disc,
     1             sigma_disc, disc_flag)
c
c****************************************************
c
c      writes AVE data to output file with file unit=unit_ave
c      returns .false. if number of output lines written exceeds
c      specified value, .true. otherwise.
c
c     11/27/95 by J. Dwyer
c                       modified  4/9/96 by J. Dwyer added year to output
c                                       5/6/96 changed outputstring to 4096 characters
c                          1/20/97 by J. Dwyer removed anisotropy data
c                          2/13/97 by J. Dwyer   added new KP list
c                          2/14/97 by J. Dwyer   removed junk, added write_kpdata
c
      include 'sdf_include.inc'   ! include type declarations
      include 'hex_include.inc'   ! include type declarations
c
      real   disc(4), sigma_disc(4)
      real   flux_ave(41), sigma_flux_ave(41)   ! arrays containing the fluxes and sigmas for all the elements
      real   vse_ave, sigma_vse_ave
      real  phi, theta,ani
      real  sigma_phi, sigma_theta, sigma_ani
      real  chi2
      real  doy_real
      real  half_delta_t  ! 1/2 of elapsed time interval (sec)
      integer   vse_ave_flag
      integer   flux_ave_flag(41)   ! .true. when flux is measured and is not saturated, .false. otherwise
      integer   disc_flag(4)
      integer    mid_time
      integer   unit_ave   ! output file unit
      integer   nrow_ave   !  number of rows written to output file
      integer      saturation_flag   ! 1 if any rate is saturated in mf, 0 otherwise
      integer   yr, mo, dy, hr, mn, se, doy
      integer    index     ! last index of outputstring written to
      character*20   isotime   ! sampex_timecon isotime array
      character*4096   outputstring
      logical   fit_flag
c
      mid_time = (stop_time_avg+start_time_avg)/2
      half_delta_t = (stop_time_avg-start_time_avg)/2.0
      call sampex_timcon(mid_time,yr,mo,dy,hr,mn,se,doy,isotime)
      doy_real = doy+(3600.0*hr+60.0*mn+se)/86400.0
c
      if (nrow_ave.lt.maximum_rows) then
c
         saturation_flag = 0
          if (vse_ave_flag.eq.-1) saturation_flag = 1
         do 10 i=1,4
            if (disc_flag(i).eq.-1) saturation_flag = 1
10     continue
         do 20 i=1,41
            if (flux_ave_flag(i).eq.-1) saturation_flag = 1
20     continue
c
         index = 1
c
         call write_line(index, outputstring,
     1          vse_ave, (vse_ave_flag.eq.1), .true.,
     1                dataformat_flag)
         call write_line(index, outputstring,
     1          sigma_vse_ave, (vse_ave_flag.eq.1), .true.,
     1                dataformat_flag)
c
         do 30 i=1,4
            call write_line(index, outputstring,
     1          disc(i), (disc_flag(i).eq.1), .true.,
     1                dataformat_flag)
            call write_line(index, outputstring,
     1          sigma_disc(i), (disc_flag(i).eq.1), .true.,
     1                dataformat_flag)
c
30      continue
         do 40 i=1,38
            call write_line(index, outputstring,
     1     flux_ave(i), (flux_ave_flag(i).eq.1), .true.,
     1                dataformat_flag)
            call write_line(index, outputstring,
     1     sigma_flux_ave(i), (flux_ave_flag(i).eq.1), .true.,
     1                dataformat_flag)
40     continue
c
      call write_kpdata(index, outputstring)
c
          write(unit_ave,100) mo, dy, yr, hr, mn, se,
     1   yr, doy_real,
     1   half_delta_t, saturation_flag,
     1   outputstring(1:(index-1))
c
100 	   format(1x,i2'/'i2'/'i4,2x,i2':'i2':'i2,
     1      ', 'i4', 'f12.8', 'f12.4', 'i1',',a<index-1>)
              nrow_ave = nrow_ave + 1
      else
         print *, 'Number of rows exceeds maximum specified:'
     1            , maximum_rows
         write_output_ave = .false.
         return
      end if
c
      write_output_ave = .true.
      return
      end   ! end write_output_ave
c


text:
enum E { X, Y }
struct S { x: U32, e: E }
constant s = { x = 1, e = E.X }
array A = [3] S default [ s, s, s ]

passive component C {
  enum E { X, Y }
  struct S { x: U32, e: E }
  constant s = { x = 1, e = E.X }
  array A = [3] S default [ s, s, s ]
}

module M {
  enum E { X, Y }
  struct S { x: U32, e: E }
  constant s = { x = 1, e = E.X }
  array A = [3] S default [ s, s, s ]
}

ext: fpp

max_stars_repo_name: kevin-f-ortega/fpp

max_issues_repo_name: kevin-f-ortega/fpp

max_forks_repo_name: kevin-f-ortega/fpp


text:
! KGEN-generated Fortran source file
!
! Filename    : rrsw_wvn.f90
! Generated at: 2015-07-07 00:48:24
! KGEN version: 0.4.13

    MODULE rrsw_wvn
        USE kgen_utils_mod, ONLY : kgen_dp, check_t, kgen_init_check, kgen_print_check
        USE shr_kind_mod, ONLY: r8 => shr_kind_r8
        !      use parkind, only : jpim, jprb
        USE parrrsw, ONLY: ngptsw
        USE parrrsw, ONLY: jpb1
        USE parrrsw, ONLY: jpb2
        USE parrrsw, ONLY: nbndsw
        IMPLICIT NONE
        !------------------------------------------------------------------
        ! rrtmg_sw spectral information
        ! Initial version:  JJMorcrette, ECMWF, jul1998
        ! Revised: MJIacono, AER, jul2006
        !------------------------------------------------------------------
        !  name     type     purpose
        ! -----  :  ----   : ----------------------------------------------
        ! ng     :  integer: Number of original g-intervals in each spectral band
        ! nspa   :  integer:
        ! nspb   :  integer:
        !wavenum1:  real   : Spectral band lower boundary in wavenumbers
        !wavenum2:  real   : Spectral band upper boundary in wavenumbers
        ! delwave:  real   : Spectral band width in wavenumbers
        !
        ! ngc    :  integer: The number of new g-intervals in each band
        ! ngs    :  integer: The cumulative sum of new g-intervals for each band
        ! ngm    :  integer: The index of each new g-interval relative to the
        !                    original 16 g-intervals in each band
        ! ngn    :  integer: The number of original g-intervals that are
        !                    combined to make each new g-intervals in each band
        ! ngb    :  integer: The band index for each new g-interval
        ! wt     :  real   : RRTM weights for the original 16 g-intervals
        ! rwgt   :  real   : Weights for combining original 16 g-intervals
        !                    (224 total) into reduced set of g-intervals
        !                    (112 total)
        !------------------------------------------------------------------
        INTEGER :: nspa(jpb1:jpb2)
        INTEGER :: nspb(jpb1:jpb2)
        REAL(KIND=r8) :: wavenum2(jpb1:jpb2)
        INTEGER :: ngc(nbndsw)
        INTEGER :: ngs(nbndsw)
        INTEGER :: ngb(ngptsw)
        PUBLIC kgen_read_externs_rrsw_wvn
    CONTAINS

    ! write subroutines

    ! module extern variables

    SUBROUTINE kgen_read_externs_rrsw_wvn(kgen_unit)
        INTEGER, INTENT(IN) :: kgen_unit
        READ(UNIT=kgen_unit) nspa
        READ(UNIT=kgen_unit) nspb
        READ(UNIT=kgen_unit) wavenum2
        READ(UNIT=kgen_unit) ngc
        READ(UNIT=kgen_unit) ngs
        READ(UNIT=kgen_unit) ngb
    END SUBROUTINE kgen_read_externs_rrsw_wvn

    END MODULE rrsw_wvn

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang


text:
//F# Compiler for F# 4.1

open System
printfn "Hello, World! How is it going?"


text:
      SUBROUTINE IB01PX( JOB, NOBR, N, M, L, UF, LDUF, UN, LDUN, UL,
     $                   LDUL, PGAL, LDPGAL, K, LDK, R, LDR, X, B, LDB,
     $                   D, LDD, TOL, IWORK, DWORK, LDWORK, IWARN,
     $                   INFO )
C
C     SLICOT RELEASE 5.5.
C
C     Copyright (c) 2002-2012 NICONET e.V.
C
C     PURPOSE
C
C     To build and solve the least squares problem  T*X = Kv,  and
C     estimate the matrices B and D of a linear time-invariant (LTI)
C     state space model, using the solution  X,  and the singular
C     value decomposition information and other intermediate results,
C     provided by other routines.
C
C     The matrix  T  is computed as a sum of Kronecker products,
C
C        T = T + kron(Uf(:,(i-1)*m+1:i*m),N_i),  for i = 1 : s,
C
C     (with  T  initialized by zero), where  Uf  is the triangular
C     factor of the QR factorization of the future input part (see
C     SLICOT Library routine IB01ND),  N_i  is given by the i-th block
C     row of the matrix
C
C            [ Q_11  Q_12  ...  Q_1,s-2  Q_1,s-1  Q_1s ]   [ I_L  0  ]
C            [ Q_12  Q_13  ...  Q_1,s-1    Q_1s    0   ]   [         ]
C        N = [ Q_13  Q_14  ...    Q_1s      0      0   ] * [         ],
C            [  :     :            :        :      :   ]   [         ]
C            [ Q_1s   0    ...     0        0      0   ]   [  0  GaL ]
C
C     and where
C
C               [   -L_1|1    ]          [ M_i-1 - L_1|i ]
C        Q_11 = [             ],  Q_1i = [               ],  i = 2:s,
C               [ I_L - L_2|1 ]          [     -L_2|i    ]
C
C     are  (n+L)-by-L  matrices, and  GaL  is built from the first  n
C     relevant singular vectors,  GaL = Un(1:L(s-1),1:n),  computed
C     by IB01ND.
C
C     The vector  Kv  is vec(K), with the matrix  K  defined by
C
C        K = [ K_1  K_2  K_3  ...  K_s ],
C
C     where  K_i = K(:,(i-1)*m+1:i*m),  i = 1:s,  is  (n+L)-by-m.
C     The given matrices are  Uf,  GaL,  and
C
C            [ L_1|1  ...  L_1|s ]
C        L = [                   ],   (n+L)-by-L*s,
C            [ L_2|1  ...  L_2|s ]
C
C        M = [ M_1  ...  M_s-1 ],      n-by-L*(s-1),  and
C        K,                            (n+L)-by-m*s.
C
C     Matrix  M  is the pseudoinverse of the matrix  GaL,  computed by
C     SLICOT Library routine IB01PD.
C
C     ARGUMENTS
C
C     Mode Parameters
C
C     JOB     CHARACTER*1
C             Specifies which of the matrices B and D should be
C             computed, as follows:
C             = 'B':  compute the matrix B, but not the matrix D;
C             = 'D':  compute both matrices B and D.
C
C     Input/Output Parameters
C
C     NOBR    (input) INTEGER
C             The number of block rows,  s,  in the input and output
C             Hankel matrices processed by other routines.  NOBR > 1.
C
C     N       (input) INTEGER
C             The order of the system.  NOBR > N > 0.
C
C     M       (input) INTEGER
C             The number of system inputs.  M >= 0.
C
C     L       (input) INTEGER
C             The number of system outputs.  L > 0.
C
C     UF      (input/output) DOUBLE PRECISION array, dimension
C             ( LDUF,M*NOBR )
C             On entry, the leading  M*NOBR-by-M*NOBR  upper triangular
C             part of this array must contain the upper triangular
C             factor of the QR factorization of the future input part,
C             as computed by SLICOT Library routine IB01ND.
C             The strict lower triangle need not be set to zero.
C             On exit, the leading  M*NOBR-by-M*NOBR  upper triangular
C             part of this array is unchanged, and the strict lower
C             triangle is set to zero.
C
C     LDUF    INTEGER
C             The leading dimension of the array  UF.
C             LDUF >= MAX( 1, M*NOBR ).
C
C     UN      (input) DOUBLE PRECISION array, dimension ( LDUN,N )
C             The leading  L*(NOBR-1)-by-N  part of this array must
C             contain the matrix  GaL,  i.e., the leading part of the
C             first  N  columns of the matrix  Un  of relevant singular
C             vectors.
C
C     LDUN    INTEGER
C             The leading dimension of the array  UN.
C             LDUN >= L*(NOBR-1).
C
C     UL      (input/output) DOUBLE PRECISION array, dimension
C             ( LDUL,L*NOBR )
C             On entry, the leading  (N+L)-by-L*NOBR  part of this array
C             must contain the given matrix  L.
C             On exit, if  M > 0,  the leading  (N+L)-by-L*NOBR  part of
C             this array is overwritten by the matrix
C             [ Q_11  Q_12  ...  Q_1,s-2  Q_1,s-1  Q_1s ].
C
C     LDUL    INTEGER
C             The leading dimension of the array  UL.  LDUL >= N+L.
C
C     PGAL    (input) DOUBLE PRECISION array, dimension
C             ( LDPGAL,L*(NOBR-1) )
C             The leading  N-by-L*(NOBR-1)  part of this array must
C             contain the pseudoinverse of the matrix  GaL,  computed by
C             SLICOT Library routine IB01PD.
C
C     LDPGAL  INTEGER
C             The leading dimension of the array  PGAL.  LDPGAL >= N.
C
C     K       (input) DOUBLE PRECISION array, dimension ( LDK,M*NOBR )
C             The leading  (N+L)-by-M*NOBR  part of this array must
C             contain the given matrix  K.
C
C     LDK     INTEGER
C             The leading dimension of the array  K.  LDK >= N+L.
C
C     R       (output) DOUBLE PRECISION array, dimension ( LDR,M*(N+L) )
C             The leading  (N+L)*M*NOBR-by-M*(N+L)  part of this array
C             contains details of the complete orthogonal factorization
C             of the coefficient matrix  T  of the least squares problem
C             which is solved for getting the system matrices B and D.
C
C     LDR     INTEGER
C             The leading dimension of the array  R.
C             LDR >= MAX( 1, (N+L)*M*NOBR ).
C
C     X       (output) DOUBLE PRECISION array, dimension
C             ( (N+L)*M*NOBR )
C             The leading  M*(N+L)  elements of this array contain the
C             least squares solution of the system  T*X = Kv.
C             The remaining elements are used as workspace (to store the
C             corresponding part of the vector Kv = vec(K)).
C
C     B       (output) DOUBLE PRECISION array, dimension ( LDB,M )
C             The leading N-by-M part of this array contains the system
C             input matrix.
C
C     LDB     INTEGER
C             The leading dimension of the array B.  LDB >= N.
C
C     D       (output) DOUBLE PRECISION array, dimension ( LDD,M )
C             If  JOB = 'D',  the leading L-by-M part of this array
C             contains the system input-output matrix.
C             If  JOB = 'B',  this array is not referenced.
C
C     LDD     INTEGER
C             The leading dimension of the array D.
C             LDD >= L, if  JOB = 'D';
C             LDD >= 1, if  JOB = 'B'.
C
C     Tolerances
C
C     TOL     DOUBLE PRECISION
C             The tolerance to be used for estimating the rank of
C             matrices. If the user sets  TOL > 0,  then the given value
C             of  TOL  is used as a lower bound for the reciprocal
C             condition number;  an m-by-n matrix whose estimated
C             condition number is less than  1/TOL  is considered to
C             be of full rank.  If the user sets  TOL <= 0,  then an
C             implicitly computed, default tolerance, defined by
C             TOLDEF = m*n*EPS,  is used instead, where  EPS  is the
C             relative machine precision (see LAPACK Library routine
C             DLAMCH).
C
C     Workspace
C
C     IWORK   INTEGER array, dimension ( M*(N+L) )
C
C     DWORK   DOUBLE PRECISION array, dimension ( LDWORK )
C             On exit, if  INFO = 0,  DWORK(1) returns the optimal value
C             of LDWORK,  and, if  M > 0,  DWORK(2)  contains the
C             reciprocal condition number of the triangular factor of
C             the matrix  T.
C             On exit, if  INFO = -26,  DWORK(1)  returns the minimum
C             value of LDWORK.
C
C     LDWORK  INTEGER
C             The length of the array DWORK.
C             LDWORK >= MAX( (N+L)*(N+L), 4*M*(N+L)+1 ).
C             For good performance,  LDWORK  should be larger.
C
C     Warning Indicator
C
C     IWARN   INTEGER
C             = 0:  no warning;
C             = 4:  the least squares problem to be solved has a
C                   rank-deficient coefficient matrix.
C
C     Error Indicator
C
C     INFO    INTEGER
C             = 0:  successful exit;
C             < 0:  if INFO = -i, the i-th argument had an illegal
C                   value.
C
C     METHOD
C
C     The matrix  T  is computed, evaluating the sum of Kronecker
C     products, and then the linear system  T*X = Kv  is solved in a
C     least squares sense. The matrices  B  and  D  are then directly
C     obtained from the least squares solution.
C
C     REFERENCES
C
C     [1] Verhaegen M., and Dewilde, P.
C         Subspace Model Identification. Part 1: The output-error
C         state-space model identification class of algorithms.
C         Int. J. Control, 56, pp. 1187-1210, 1992.
C
C     [2] Van Overschee, P., and De Moor, B.
C         N4SID: Two Subspace Algorithms for the Identification
C         of Combined Deterministic-Stochastic Systems.
C         Automatica, Vol.30, No.1, pp. 75-93, 1994.
C
C     [3] Van Overschee, P.
C         Subspace Identification : Theory - Implementation -
C         Applications.
C         Ph. D. Thesis, Department of Electrical Engineering,
C         Katholieke Universiteit Leuven, Belgium, Feb. 1995.
C
C     NUMERICAL ASPECTS
C
C     The implemented method is numerically stable.
C
C     CONTRIBUTOR
C
C     V. Sima, Katholieke Universiteit Leuven, Feb. 2000.
C
C     REVISIONS
C
C     V. Sima, Katholieke Universiteit Leuven, Sep. 2001.
C
C     KEYWORDS
C
C     Identification methods; least squares solutions; multivariable
C     systems; QR decomposition; singular value decomposition.
C
C     ******************************************************************
C
C     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE
      PARAMETER          ( ZERO  = 0.0D0, ONE = 1.0D0 )
C     .. Scalar Arguments ..
      DOUBLE PRECISION   TOL
      INTEGER            INFO, IWARN, L, LDB, LDD, LDK, LDPGAL, LDR,
     $                   LDUF, LDUL, LDUN, LDWORK, M, N, NOBR
      CHARACTER          JOB
C     .. Array Arguments ..
      DOUBLE PRECISION   B(LDB, *), D(LDD, *), DWORK(*), K(LDK, *),
     $                   PGAL(LDPGAL, *), R(LDR, *), UF(LDUF, *),
     $                   UL(LDUL, *), UN(LDUN, *), X(*)
      INTEGER            IWORK( * )
C     .. Local Scalars ..
      DOUBLE PRECISION   RCOND, TOLL
      INTEGER            I, IERR, J, JWORK, LDUN2, LNOBR, LP1, MAXWRK,
     $                   MINWRK, MKRON, MNOBR, NKRON, NP1, NPL, RANK
      LOGICAL            WITHB, WITHD
C     .. External Functions ..
      LOGICAL            LSAME
      DOUBLE PRECISION   DLAMCH
      EXTERNAL           DLAMCH, LSAME
C     .. External Subroutines ..
      EXTERNAL           DGELSY, DGEMM, DLACPY, DLASET, DTRCON, MB01VD,
     $                   XERBLA
C     .. Intrinsic Functions ..
      INTRINSIC          MAX
C     .. Executable Statements ..
C
C     Decode the scalar input parameters.
C
      WITHD = LSAME( JOB, 'D' )
      WITHB = LSAME( JOB, 'B' ) .OR. WITHD
      MNOBR = M*NOBR
      LNOBR = L*NOBR
      LDUN2 = LNOBR - L
      LP1   = L + 1
      NP1   = N + 1
      NPL   = N + L
      IWARN = 0
      INFO  = 0
C
C     Check the scalar input parameters.
C
      IF( .NOT.WITHB ) THEN
         INFO = -1
      ELSE IF( NOBR.LE.1 ) THEN
         INFO = -2
      ELSE IF( N.GE.NOBR .OR. N.LE.0 ) THEN
         INFO = -3
      ELSE IF( M.LT.0 ) THEN
         INFO = -4
      ELSE IF( L.LE.0 ) THEN
         INFO = -5
      ELSE IF( LDUF.LT.MAX( 1, MNOBR ) ) THEN
         INFO = -7
      ELSE IF( LDUN.LT.LDUN2 ) THEN
         INFO = -9
      ELSE IF( LDUL.LT.NPL ) THEN
         INFO = -11
      ELSE IF( LDPGAL.LT.N ) THEN
         INFO = -13
      ELSE IF( LDK.LT.NPL ) THEN
         INFO = -15
      ELSE IF( LDR.LT.MAX( 1, MNOBR*NPL ) ) THEN
         INFO = -17
      ELSE IF( LDB.LT.N ) THEN
         INFO = -20
      ELSE IF( LDD.LT.1 .OR. ( WITHD .AND. LDD.LT.L ) ) THEN
         INFO = -22
      ELSE
C
C        Compute workspace.
C        (Note: Comments in the code beginning "Workspace:" describe the
C         minimal amount of workspace needed at that point in the code,
C         as well as the preferred amount for good performance.
C         NB refers to the optimal block size for the immediately
C         following subroutine, as returned by ILAENV.)
C
         MINWRK = MAX( NPL*NPL, 4*M*NPL + 1 )
C
         IF ( LDWORK.LT.MINWRK ) THEN
            INFO = -26
            DWORK( 1 ) = MINWRK
         END IF
      END IF
C
C     Return if there are illegal arguments.
C
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'IB01PX', -INFO )
         RETURN
      END IF
C
C     Quick return if possible.
C
      IF( M.EQ.0 ) THEN
         DWORK(1) = ONE
         RETURN
      END IF
C
C     Construct the matrix  [ Q_11  Q_12  ...  Q_1,s-1  Q_1s ]  in  UL.
C
      DO 20 J = 1, L
C
         DO 10 I = 1, NPL
            UL(I,J) = -UL(I,J)
   10    CONTINUE
C
         UL(N+J,J) = ONE + UL(N+J,J)
   20 CONTINUE
C
      DO 50 J = LP1, LNOBR
C
         DO 30 I = 1, N
            UL(I,J) = PGAL(I,J-L) - UL(I,J)
   30    CONTINUE
C
         DO 40 I = NP1, NPL
            UL(I,J) = -UL(I,J)
   40    CONTINUE
C
   50 CONTINUE
C
C     Compute the coefficient matrix T using Kronecker products.
C     Workspace: (N+L)*(N+L).
C     In the same loop, vectorize K in X.
C
      CALL DLASET( 'Full', MNOBR*NPL, M*NPL, ZERO, ZERO, R, LDR )
      CALL DLASET( 'Lower', MNOBR-1, MNOBR-1, ZERO, ZERO, UF(2,1),
     $             LDUF )
      JWORK = NPL*L + 1
C
      DO 60 I = 1, NOBR
         CALL DLACPY( 'Full', NPL, L, UL(1,(I-1)*L+1), LDUL, DWORK,
     $                NPL )
         IF ( I.LT.NOBR ) THEN
            CALL DGEMM ( 'NoTranspose', 'NoTranspose', NPL, N,
     $                   L*(NOBR-I), ONE, UL(1,I*L+1), LDUL, UN, LDUN,
     $                   ZERO, DWORK(JWORK), NPL )
         ELSE
            CALL DLASET( 'Full', NPL, N, ZERO, ZERO, DWORK(JWORK), NPL )
         END IF
         CALL MB01VD( 'NoTranspose', 'NoTranspose', MNOBR, M, NPL,
     $                NPL, ONE, ONE, UF(1,(I-1)*M+1), LDUF, DWORK,
     $                NPL, R, LDR, MKRON, NKRON, IERR )
         CALL DLACPY( 'Full', NPL, M, K(1,(I-1)*M+1), LDK,
     $                X((I-1)*NKRON+1), NPL )
   60 CONTINUE
C
C     Compute the tolerance.
C
      TOLL = TOL
      IF( TOLL.LE.ZERO )
     $   TOLL = MKRON*NKRON*DLAMCH( 'Precision' )
C
C     Solve the least square problem  T*X = vec(K).
C     Workspace:  need   4*M*(N+L)+1;
C                 prefer 3*M*(N+L)+(M*(N+L)+1)*NB.
C
      DO 70 I = 1, NKRON
         IWORK(I) = 0
   70 CONTINUE
C
      CALL DGELSY( MKRON, NKRON, 1, R, LDR, X, MKRON, IWORK, TOLL, RANK,
     $             DWORK, LDWORK, IERR )
      MAXWRK = DWORK(1)
C
C     Compute the reciprocal of the condition number of the triangular
C     factor  R  of  T.
C     Workspace: need 3*M*(N+L).
C
      CALL DTRCON( '1-norm', 'Upper', 'NonUnit', NKRON, R, LDR, RCOND,
     $             DWORK, IWORK, IERR )
C
      IF ( RANK.LT.NKRON ) THEN
C
C        The least squares problem is rank-deficient.
C
         IWARN = 4
      END IF
C
C     Construct the matrix  D,  if needed.
C
      IF ( WITHD )
     $   CALL DLACPY( 'Full', L, M, X, NPL, D, LDD )
C
C     Construct the matrix  B.
C
      CALL DLACPY( 'Full', N, M, X(LP1), NPL, B, LDB )
C
C     Return optimal workspace in DWORK(1) and reciprocal condition
C     number in  DWORK(2).
C
      DWORK(1) = MAX( MINWRK, MAXWRK )
      DWORK(2) = RCOND
C
      RETURN
C
C *** Last line of IB01PX ***
      END

max_stars_repo_name: bgin/MissileSimulation

max_issues_repo_name: bgin/MissileSimulation

max_forks_repo_name: bgin/MissileSimulation


module internal_c
  implicit none

  public

  interface

    subroutine getuser_c(user, stat) bind(C)
      use, intrinsic :: ISO_C_Binding
      character(kind=C_CHAR), dimension(*), intent(out) :: user
      integer(C_INT),                       intent(out) :: stat
    end subroutine getuser_c

    subroutine gethome_c(home, stat) bind(C)
      use, intrinsic :: ISO_C_Binding
      character(kind=C_CHAR), dimension(*), intent(out) :: home
      integer(C_INT),                       intent(out) :: stat
    end subroutine gethome_c

    function PATH_MAX_c() bind(C)
      use, intrinsic :: ISO_C_Binding
      integer(C_INT) :: PATH_MAX_c
    end function PATH_MAX_c

  end interface

end module internal_c


text:
module Positions
    use Const
    implicit none

    type Pos
        real(mpc), allocatable :: mem(:)
        contains
            procedure :: x  => pos_x
            procedure :: p  => pos_p
            procedure :: q  => pos_q

            procedure :: Mx => pos_set_x
            procedure :: Mp => pos_set_p
            procedure :: Mq => pos_set_q

            procedure :: assign_pos
            generic :: assignment(=) => assign_pos
!             final :: clear_pos
    end type Pos

    interface Pos
        module procedure :: init_pos_pq, init_pos_x
    end interface Pos

    interface operator(+)
        module procedure :: add_pos_pos
    end interface
    interface operator(*)
        module procedure :: mult_real_pos, mult_pos_real
    end interface
!     interface assignment(=)
!         module procedure :: pos_assign
!     end interface

contains

    function init_pos_x(x) result(c)
        real(mpc), intent(in), dimension(:) :: x

        type(Pos) :: c
        integer :: memsize

        memsize = size(x)
        allocate(c%mem(memsize))
        c%mem = x
    end function init_pos_x

    function init_pos_pq(p,q) result(c)
        real(mpc), intent(in), dimension(:) :: p, q

        type(Pos) :: c
        integer :: memsize, dimen

        if (size(p) /= size(q)) stop "incompatible dimensions"

        memsize = size(p) + size(q)
        dimen = size(p)
        allocate(c%mem(memsize))

        c%mem(1:dimen) = p
        c%mem(dimen+1:memsize) = q
    end function init_pos_pq

!=============================
!         getters
!=============================
    pure function pos_x(self) result(v)
        class(Pos), intent(in) :: self
        real(mpc) :: v((size(self%mem)))

        v = self%mem
    end function pos_x

    pure function pos_p(self) result(v)
        class(Pos), intent(in) :: self
        real(mpc) :: v((size(self%mem))/2)
        integer :: half

        half = (size(self%mem))/2
        v = self%mem(1:half)
    end function pos_p

    pure function pos_q(self) result(v)
        class(Pos), intent(in) :: self
        real(mpc) :: v((size(self%mem))/2)
        integer :: half

        half = (size(self%mem))/2
        v = self%mem(half+1:2*half)
    end function pos_q

!=================================
!         setters
!=================================

    subroutine pos_set_x(self, v)
        class(Pos), intent(inout) :: self
        real(mpc), intent(in) :: v(:)

        self%mem(1:size(self%mem)) = v
    end subroutine pos_set_x

    subroutine pos_set_p(self, v)
        class(Pos), intent(inout) :: self
        real(mpc), intent(in) :: v(:)
        integer :: half

        half = (size(self%mem))/2
        self%mem(1:half) = v
    end subroutine pos_set_p

    subroutine pos_set_q(self, v)
        class(Pos), intent(inout) :: self
        real(mpc), intent(in) :: v(:)
        integer :: half

        half = (size(self%mem))/2
        self%mem(half+1:2*half) = v
    end subroutine pos_set_q

    subroutine clear_pos(self)
        type(Pos), intent(inout) :: self
        if (allocated(self%mem)) deallocate(self%mem)
    end subroutine clear_pos

! **************************
! overloaded operators
! **************************

    function add_pos_pos(p1, p2) result(p)
        type(Pos), intent(in) :: p1, p2
        type(Pos) :: p
        p = Pos(p1%p() + p2%p(), p1%q() + p2%q())
    end function add_pos_pos

    function mult_real_pos(r, p) result(p1)
        type(Pos), intent(in) :: p
        real(mpc), intent(in) :: r
        type(Pos) :: p1
        p1 = Pos(r * p%p(), r * p%q())
    end function mult_real_pos

    function mult_pos_real(p, r) result(p1)
        type(pos), intent(in) :: p
        real(mpc), intent(in) :: r
        type(pos) :: p1
        p1 = mult_real_pos(r, p)
    end function mult_pos_real

    subroutine assign_pos(p, p1)
        class(Pos), intent(out) :: p
        class(Pos), intent(in) :: p1
        allocate(p%mem(size(p1%mem)), source=p1%mem)
    end subroutine assign_pos

    subroutine test_pos()
        type(Pos) :: p1, p2, p
        p1 = Pos((/1.0_mpc,1.0_mpc/), (/0.1_mpc, 0.0_mpc/))
        p2 = Pos((/2.0_mpc,2.0_mpc/), (/0.2_mpc, 0.0_mpc/))
        write(*,*) "q =", p1%q()
        write(*,*) "x =", p1%x()
        call p1%Mx([1.0_mpc, 1.1_mpc, 0.0_mpc, 0.1_mpc])
        p = p1 + p2
        write(*,*) "p' =", p%p(), "q' =", p%q()
        p = p1
        write(*,*) "p' =", p%p(), "q' =", p%q()
        call p1%Mp([0.0_mpc, 0.0_mpc])
        write(*,*) "p' =", p1%p(), "p' =", p%p()
    end subroutine test_pos

end module

max_stars_repo_name: taxus-d/integrators

max_issues_repo_name: taxus-d/integrators

max_forks_repo_name: taxus-d/integrators


text:
program loop_test

      implicit none
      integer   :: i
      complex   :: C,D
      C=(0,0)
      D=(1,2)

      !$OMP TARGET PARALLEL DO REDUCTION(+:C)
      do i=1, 10
         C=C+EXP(D)
      end do
      !$OMP  END TARGET PARALLEL DO

      write(*, *) "C= ", C

end program loop_test

max_stars_repo_name: aperdeus-amd/aomp

max_issues_repo_name: aperdeus-amd/aomp

max_forks_repo_name: aperdeus-amd/aomp


text:
Jacksonfly Recycled Fashion was a thrift stores vintage/thrift clothing store located downtown next to Ace Hardware.  As they put it, Were new... were hidden... we have a lot of inventory... we are priced low... we are rad... or so we think... come and see for yourself!  It is now The Rad Vintager.

They sold new British bras, new & used designer denim, vintage tops, kids clothes, funky tshirts,dresses, and $5 Hawaiian shirts.

20101103 21:47:02 nbsp I visited this store for the first time last weekend and found it very pleasant. I like the idea of a higher quality used clothing store with some new things for sale in down town Davis. The woman behind the counter (the owner?) was friendly but not overly chatty or pushy which I appreciate when Im shopping for any clothes. I think this the kind of store that you have to check back on periodically because they get new things in often. Users/MaryCarroll

20101114 21:23:01 nbsp I really like this place. I found a fantastic selection of likenew used clothing for extremely reasonable prices. The man behind the counter was extremely nice and very helpful. They also trade and buy clothing! He mentioned that he would be more open to trade right now, which is understandable, since the store opened recently. But it might be nice to get a few bucks for the nicer things I would otherwise give away. I am definitely coming back for the great selection, value, and service.  Users/LizzyT

20110131 10:38:53 nbsp This weekend I stopped in Jacksonfly to take advantage of Davis 30$ on 1/30 day. I was pleased to find out that the entire store was 50% off, rather than just 30% off. This translated into a brand new pair of Lucky Brand jeans which were marked at $25 only costing me $12.50. I found lots of other cute things too, and even at full price they would have been a steal. There are lots of unique and vintage items as well as more popular brands like Banana Republic, Ann Taylor, Levi, Forever 21, etc.. The woman working there was very helpful in finding things my size, but was far from pushy like the previous poster said. The only thing to keep in mind when shopping here is to bring cash/checks, as they dont take debit/visa cards. Since I had to run up the street to get cash from my bank, the woman working at the store let me pick something out for free for my trouble. Worth it!! I will definitely be stopping in here frequently to check out their inventory. You should too!! Users/LaraStilling

20110131 15:32:26 nbsp This place is great! I went there on Sunday and everything was 50% off. The prices are really low and the selection is great. They dont have an overwhelmingly large amount of clothing in the store and everything there is cute. The woman working, Amanda, was really nice! She helped me find clothes in my size and gave me her opinion on the stuff I tried on (when I asked, of course). Users/RachelKerwin

20110208 09:29:09 nbsp This place is just what Davis needs. They have great customer service, even though they only take cash. However, I didnt mind because I found a silk Oscar De La Renta top for only $6!!!! Yay. Users/MelannieLavarias

20110208 09:53:12 nbsp was wondering: do they have larger sizes? Im about a 12 or sometimes 14 Users/jsbmeb

20110504 12:12:18 nbsp Just wanted to echo what is being said on their facebook the store in its previous location, by Peets Coffee and Tea, has closed, and the store has moved to the new location listed above. They should be open by this weekend (5/7 or 8ish/2011), supposedly. Users/AndrewJacobs


text:
*
*     ------------------------------------------------------------------
*	G T R A C 1
*     ------------------------------------------------------------------
*
      SUBROUTINE GTRAC1(I,L,IFIRST,NFOUND,BUF,
     &           IBUF,LBUF,JLEI,LU,NTESTG)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
*
* Obtain racah coefficient for given I and l starting
* from element IFIRST
*
* JLEI = 0 => include all terms with given I
* JLEI = 1 => include only terms with J<I
*
*
C     write(6,*) ' Old GTRACX '
C     CALL  GTRACX(I,L,IFIRST,NFOUND,BUF,
C    &                IBUF,LBUF,JLEI,LU,NTESTG)
*
C     write(6,*) ' new gtracx '
C     CALL  GTRACXN(I,L,IFIRST,NFOUND,BUF,
C    &                IBUF,LBUF,JLEI,LU,NTESTG)
* For new version where the coupling coefficients reside in core,
* the program must know which of the two lists that should be read.
* This is done by assuming that LU is as ' in the bad old days'
* i.e
*         LU = 15 => Left (initial state)
*         LU = 16 => right (final state)
      IF(LU.EQ.15) THEN
        ILIST = 1
      ELSE IF (LU.EQ.16) THEN
        ILIST = 2
      ELSE
        WRITE(6,*) ' PROBLEM in GTRAC1, Unrecognized LU = ',LU
        STOP'GTRAC1: Undefined LU'
      END IF
*
      CALL  GTRACXVN(I,L,IFIRST,NFOUND,BUF,
     &                IBUF,LBUF,JLEI,ILIST,NTESTG)
*
      RETURN
      END


text:
module example

  use, intrinsic :: iso_c_binding, only: &
       c_double, c_int, c_ptr, c_intptr_t, c_char, c_f_pointer
  implicit none
  private
  public &
       dp, foo, foo_array, foo_by_ref, make_udf, udf_ptr, just_print, &
       view_knob, turn_knob, UserDefined

  integer, parameter :: dp = kind(0.d0)
  integer(c_int) :: KNOB = 1337

  type, bind(c) :: UserDefined
     ! sequence NOT NEEDED since bind(c) is used.
     real(c_double) :: buzz
     real(c_double) :: broken
     integer(c_int) :: how_many
  end type UserDefined

  type, bind(c) :: DataContainer
     real(c_double) :: data(4, 2)
  end type DataContainer

contains

  subroutine foo(bar, baz, quux) bind(c, name='foo')
    real(c_double), intent(in), value :: bar, baz
    real(c_double), intent(out) :: quux

    quux = bar + 3.75_dp * baz

  end subroutine foo

  subroutine foo_by_ref(bar, baz, quux)
    real(dp), intent(in) :: bar, baz
    real(dp), intent(out) :: quux

    call foo(bar, baz, quux)

  end subroutine foo_by_ref

  subroutine foo_array(size_, val, two_val) bind(c, name='foo_array')
    integer(c_int), intent(in) :: size_
    real(c_double), intent(in) :: val(size_, 2)
    real(c_double), intent(out) :: two_val(size_, 2)

    two_val = 2.0_dp * val

  end subroutine foo_array

  subroutine make_udf(buzz, broken, how_many, as_bytes) bind(c, name='make_udf')
    real(c_double), intent(in) :: buzz, broken
    integer(c_int), intent(in) :: how_many
    character(c_char), intent(out) :: as_bytes(24)
    ! Outside of signature
    type(UserDefined), target :: made_it

    made_it%buzz = buzz
    made_it%broken = broken
    made_it%how_many = how_many

    ! NOTE: We need ``sizeof(as_bytes) == sizeof(made_it)``
    as_bytes = transfer(made_it, as_bytes)

  end subroutine make_udf

  subroutine udf_ptr(ptr_as_int) bind(c, name='udf_ptr')
    integer(c_intptr_t), intent(in) :: ptr_as_int
    ! Outside of signature
    type(c_ptr) :: made_it_ptr
    type(UserDefined), pointer :: made_it

    made_it_ptr = transfer(ptr_as_int, made_it_ptr)
    call c_f_pointer(made_it_ptr, made_it)

    made_it%buzz = 3.125_dp
    made_it%broken = -10.5_dp
    made_it%how_many = 101

  end subroutine udf_ptr

  subroutine make_container(contained, container) &
       bind(c, name='make_container')
    real(c_double), intent(in) :: contained(4, 2)
    type(DataContainer), intent(out) :: container

    container%data = contained

  end subroutine make_container

  subroutine just_print() bind(c, name='just_print')

    print *, "======== BEGIN FORTRAN ========"
    print *, "just_print() was called"
    print *, "========  END  FORTRAN ========"

  end subroutine just_print

  ! NOTE: This is missing ``bind(c, name='view_knob')`` because
  !       the ``f2py`` parser has issues for some reason.
  pure function view_knob() result(knob_value)
    integer(c_int) :: knob_value

    knob_value = KNOB

  end function view_knob

  subroutine turn_knob(new_value) bind(c, name='turn_knob')
    integer(c_int), intent(in) :: new_value

    KNOB = new_value

  end subroutine turn_knob

end module example

max_stars_repo_name: dhermes/foreign-fortran

max_issues_repo_name: dhermes/foreign-fortran

max_forks_repo_name: dhermes/foreign-fortran


text:
      SUBROUTINE BMG3_SymStd_UTILS_matvec( &
     &                       KG, SO, QF, Q, II, JJ, KK,&
     &                       NOG, IFD, NStncl, halof&
     &                       ) BIND(C, NAME='MPI_BMG3_SymStd_UTILS_matvec')

! ======================================================================
!  --------------------
!   DESCRIPTION:
!  --------------------
!
!     BMG3_SymStd_UTILS_matvec performs a matrix multiplication that is
!     used within the pcg routine called BMG3_SymStd_SOLVE_pcg.
!
! ======================================================================
! $license_flag$
! ======================================================================
!  --------------------
!   INPUT:
!  --------------------
!
!
! ======================================================================
!  --------------------
!   INPUT/OUTPUT:
!  --------------------
!
!
! ======================================================================
!  --------------------
!   OUTPUT:
!  --------------------
!
!
! ======================================================================
!  --------------------
!   LOCAL:
!  --------------------
!
!
! ======================================================================
      USE ModInterface
      IMPLICIT NONE

! -----------------------------
!     Includes
!
      INCLUDE 'mpif.h'
      INCLUDE 'MSG_f90.h'

      INCLUDE 'BMG_constants_f90.h'
      INCLUDE 'BMG_stencils_f90.h'
      INCLUDE 'BMG_workspace_f90.h'
      INCLUDE 'BMG_parameters_f90.h'

! ----------------------------
!     Argument Declarations
!

      integer(len_t), VALUE :: II, JJ, KK
      integer(c_int), VALUE :: IFD, KG, NOG,&
           NStncl
      real(real_t) :: Q(II,JJ,KK), QF(II,JJ,KK)
      real(real_t) :: SO(II+1,JJ+1,KK+1,NStncl)
      type(c_ptr) :: halof

! ----------------------------
!     Local Declarations
!
      integer i, i1, j, j1, k, k1

! ======================================================================

      i1=II-1
      j1=JJ-1
      k1=KK-1

      IF( KG.LT.NOG .OR. IFD.NE.BMG_STENCIL_7pt ) THEN

         DO k=2,k1
            DO j=2,j1
               DO i=2,i1
                  QF(i,j,k) =  SO( i , j , k ,kp  )*Q( i , j , k )&
     &                       - SO( i , j , k ,kpw )*Q(i-1, j , k )&
     &                       - SO( i ,j+1, k ,kpnw)*Q(i-1,j+1, k )&
     &                       - SO( i ,j+1, k ,kps )*Q( i ,j+1, k )&
     &                       - SO(i+1,j+1, k ,kpsw)*Q(i+1,j+1, k )&
     &                       - SO(i+1, j , k ,kpw )*Q(i+1, j , k )&
     &                       - SO(i+1, j , k ,kpnw)*Q(i+1,j-1, k )&
     &                       - SO( i , j , k ,kps )*Q( i ,j-1, k )&
     &                       - SO( i , j , k ,kpsw)*Q(i-1,j-1, k )&
     &                       - SO( i , j , k ,kb  )*Q( i , j ,k-1)&
     &                       - SO( i , j , k ,kbw )*Q(i-1, j ,k-1)&
     &                       - SO( i ,j+1, k ,kbnw)*Q(i-1,j+1,k-1)&
     &                       - SO( i ,j+1, k ,kbn )*Q( i ,j+1,k-1)&
     &                       - SO(i+1,j+1, k ,kbne)*Q(i+1,j+1,k-1)&
     &                       - SO(i+1, j , k ,kbe )*Q(i+1, j ,k-1)&
     &                       - SO(i+1, j , k ,kbse)*Q(i+1,j-1,k-1)&
     &                       - SO( i , j , k ,kbs )*Q( i ,j-1,k-1)&
     &                       - SO( i , j , k ,kbsw)*Q(i-1,j-1,k-1)&
     &                       - SO( i , j ,k+1,kb  )*Q( i , j ,k+1)&
     &                       - SO( i , j ,k+1,kbe )*Q(i-1, j ,k+1)&
     &                       - SO( i ,j+1,k+1,kbse)*Q(i-1,j+1,k+1)&
     &                       - SO( i ,j+1,k+1,kbs )*Q( i ,j+1,k+1)&
     &                       - SO(i+1,j+1,k+1,kbsw)*Q(i+1,j+1,k+1)&
     &                       - SO(i+1, j ,k+1,kbw )*Q(i+1, j ,k+1)&
     &                       - SO(i+1, j ,k+1,kbnw)*Q(i+1,j-1,k+1)&
     &                       - SO( i , j ,k+1,kbn )*Q( i ,j-1,k+1)&
     &                       - SO( i , j ,k+1,kbne)*Q(i-1,j-1,k+1)
               ENDDO
            ENDDO
         ENDDO

      ELSEIF( KG.EQ.NOG .AND. IFD.EQ.BMG_STENCIL_7pt ) THEN

         DO k=2,k1
            DO j=2,j1
               DO i=2,i1
                  QF(i,j,k) =  SO( i , j , k ,kp )*Q( i , j , k )&
     &                       - SO( i , j , k ,kpw)*Q(i-1, j , k )&
     &                       - SO( i ,j+1, k ,kps)*Q( i ,j+1, k )&
     &                       - SO(i+1, j , k ,kpw)*Q(i+1, j , k )&
     &                       - SO( i , j , k ,kps)*Q( i ,j-1, k )&
     &                       - SO( i , j , k ,kb )*Q( i , j ,k-1)&
     &                       - SO( i , j ,k+1,kb )*Q( i , j ,k+1)
               ENDDO
            ENDDO
         ENDDO

      ENDIF

! ======================================================================

      call halo_exchange(KG, Q, halof)

! ======================================================================

      RETURN
      END

max_stars_repo_name: cedar-framework/cedar

max_issues_repo_name: cedar-framework/cedar

max_forks_repo_name: cedar-framework/cedar


text:
! This software may be used, reproduced, and provided to others only as permitted under the terms of the agreement
! under which it was acquired from the U.S. Government. Neither title to, nor ownership of, the software is hereby
! transferred. This notice shall remain on all copies of the software.
!
! Copyright 2016 United States Government as represented by the Administrator of the National Aeronautics and
! Space Administration. No copyright is claimed in the United States under Title 17, U.S. Code. All Other Rights
! Reserved.

#include "vumatArgs.for"
#include "version.for"
#include "forlog.for"
#include "matrixUtil.for"
#include "matProp.for"
#include "stateVar.for"
#include "parameters.for"
#include "schapery.for"
#include "stress.for"
#include "strain.for"
#include "schaefer.for"
#include "plasticity.for"
#include "fiberDamage.for"
#include "friction.for"
#include "cohesive.for"
#include "vucharlength.for"
#include "DGD.for"
#include "vexternaldb.for"

Subroutine VUMAT(  &
  ! Read only (unmodifiable) variables:
  jblock,ndir,nshr,nstatev,nfieldv,nprops,lanneal,stepTime,     &
  totalTime,dt,cmname,coordMp,charLength,props,density,         &
  strainInc,relSpinInc,tempOld,stretchOld,defgradOld,fieldOld,  &
  stressOld,stateOld,enerInternOld,enerInelasOld,tempNew,       &
  stretchNew,defgradNew,fieldNew,                               &
  ! Write only (modifiable) variables:
  stressNew,stateNew,enerInternNew,enerInelasNew)

  Implicit Double Precision (a-h, o-z)
  Parameter (j_sys_Dimension = 2, n_vec_Length = 136, maxblk = n_vec_Length)

  Dimension jblock(*),props(*),density(*),coordMp(*),charLength(*),strainInc(*),relSpinInc(*),tempOld(*),tempNew(*), &
    stretchOld(*),defgradOld(*),fieldOld(*),stressOld(*),stateOld(*),enerInternOld(*),                               &
    enerInelasOld(*),stretchNew(*),defgradNew(*),fieldNew(*),stressNew(*),stateNew(*),                               &
    enerInternNew(*),enerInelasNew(*)

  Character(len=80) :: cmname

  Integer, parameter :: i_nblock=1, i_npt=2, i_layer=3, i_kspt=4, i_noel=5

  ! -------------------------------------------------------------------- !
  !    End VUMAT standard interface                                      !
  ! -------------------------------------------------------------------- !

  ! Throw error on single precision
  If (PRECISION(dt) < 15) Then
    Print *, 'ERROR: Found single precision data. CompDam must be run in double precision.'
    Call XPLB_EXIT
  Else
    Call CompDam(  &
      ! Read only (unmodifiable) variables:
      jblock(i_nblock),ndir,nshr,nstatev,nfieldv,nprops,lanneal,stepTime,     &
      totalTime,dt,cmname,coordMp,charLength,props,density,         &
      strainInc,relSpinInc,tempOld,stretchOld,defgradOld,fieldOld,  &
      stressOld,stateOld,enerInternOld,enerInelasOld,tempNew,       &
      stretchNew,defgradNew,fieldNew,                               &
      ! Write only (modifiable) variables:
      stressNew,stateNew,enerInternNew,enerInelasNew,               &
      ! Read only, extra arguments:
      jblock(i_noel), jblock(i_npt), jblock(i_layer), jblock(i_kspt))
  End If
End Subroutine VUMAT

Subroutine CompDam(  &
  ! Read only (unmodifiable) variables:
  nblock,ndir,nshr,nstatev,nfieldv,nprops,lanneal,stepTime,     &
  totalTime,dt,cmname,coordMp,charLength,props,density_abq,     &
  strainInc,relSpinInc,tempOld,stretchOld,defgradOld,fieldOld,  &
  stressOld,stateOld,enerInternOld,enerInelasOld,tempNew,       &
  stretchNew,defgradNew,fieldNew,                               &
  ! Write only (modifiable) variables:
  stressNew,stateNew,enerInternNew,enerInelasNew,               &
  ! Read only, extra arguments:
  nElement, nMatPoint, nLayer, nSecPoint)

  Use forlog_Mod
  Use matrixAlgUtil_Mod
  Use matProp_Mod
  Use stateVar_Mod
  Use parameters_Mod
  Use DGD_Mod
  Use plasticity_mod
  Use cohesive_mod
  Use friction_mod

  Implicit Double Precision (a-h, o-z)
  Parameter (j_sys_Dimension = 2, n_vec_Length = 136, maxblk = n_vec_Length)

  Double Precision :: props(nprops), density_abq(nblock), coordMp(nblock,*), charLength(nblock,*),       &
    strainInc(nblock,ndir+nshr), relSpinInc(nblock,nshr), tempOld(nblock), tempNew(nblock),              &
    stretchOld(nblock,ndir+nshr), defgradOld(nblock,ndir+nshr+nshr), fieldOld(nblock,nfieldv),           &
    stressOld(nblock,ndir+nshr), stateOld(nblock,nstatev), enerInternOld(nblock), enerInelasOld(nblock), &
    stretchNew(nblock,ndir+nshr), defgradNew(nblock,ndir+nshr+nshr), fieldNew(nblock,nfieldv),           &
    stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev), enerInternNew(nblock), enerInelasNew(nblock)

  ! Extra arguments
  Integer :: nElement(nblock), nMatPoint, nLayer, nSecPoint

  Double Precision :: stepTime, totalTime, dt

  Character(len=80) :: cmname

  ! -------------------------------------------------------------------- !
  !    End VUMAT standard interface                                      !
  ! -------------------------------------------------------------------- !
  ! Damage variables, failure indices, damage thresholds
  Double Precision :: d1C_dummy, rfC_dummy

  ! Stress
  Double Precision :: Cauchy(3,3)
  Double Precision :: CauchyABQ(3,3)
  Double Precision :: T_coh(3)  ! Cohesive tractions (shear, normal, shear)
  Double Precision :: Rot(3,3)

  ! Other
  Double Precision :: F(3,3)                                               ! Current deformation gradient tensor
  Double Precision :: F_old(3,3)                                           ! Previous deformation gradient tensor
  Double Precision :: U(3,3)
  Double Precision :: Pen(3)  ! Cohesive penalty stiffness
  Double Precision :: delta(3)  ! Cohesive displacement-jump vector
  Double Precision :: dmg_penalty
  Double Precision :: dGdGc  ! incremental change in normalized energy dissipation
  Double Precision :: density_current         ! Current density calculated from Fbulk (ignores the cohesive crack displacement)
  Logical :: Sliding

  ! Parameters
  Double Precision, parameter :: zero=0.d0, one=1.d0, two=2.d0

  ! Structure of all VUMAT arguments
  Type(vumatArg) :: args

  ! For access to the logger
  Type(forlog) :: logger

  ! Material
  Type(matProps) :: m
  Type(materialList), save :: user
  Type(materialList) :: user_temp
  Integer :: mats

  ! State variables
  Type(stateVars) :: sv

  ! Solution parameters
  Type(parameters), save :: p

  ! Fatigue
  Logical :: fatigue_step
  Integer :: fatigue_parameters(2)
  Pointer(ptr_fatigue_int, fatigue_parameters)

  Interface

      Function SMAIntArrayAccess(ID)
        ! Access an existing global integer array.
        Integer(kind=8) :: SMAIntArrayAccess  ! Returns an address that can be associated with a Fortran pointer
        Integer(kind=4) :: ID       ! Array ID
      End Function SMAIntArrayAccess

  End Interface
  ! -------------------------------------------------------------------- !

  ! Initialize structure of VUMAT args
  Call args%init(nblock, ndir, nshr, nstatev, stepTime, totalTime, dt)

  ! Only write git hash on first call
  If (stepTime <= dt .AND. p%logLevel == 0) Then

    ! Initialize the logger for use in loadParameters()
    Call log%init(level=2, vumatArgStruct=args, format=1)

    ! Load the CompDam solution parameters
    p = loadParameters()

  End If

  ! Initialize the logger
  Call log%init(level=p%logLevel, vumatArgStruct=args, format=p%logFormat)

  If (stepTime <= dt) Then

    If (Allocated(user%materials)) Then
      mats = Size(user%materials)
      readMaterial: Do I = 1,mats
        If (user%materials(I)%name == trim(cmname)) Then
          m = user%materials(I)
          Exit readMaterial
        End If

        If (I == mats) Then
          Call log%info("Saving user material #"//trim(str(I+1))//": "//trim(cmname))

          ! Store all existing material data in temp array
          Allocate (user_temp%materials(mats))
          user_temp%materials(:) = user%materials(:)

          ! Increase size of user%materials by 1
          Deallocate (user%materials)
          Allocate (user%materials(mats+1))

          ! Return old material data to resized user%materials
          user%materials(1:mats) = user_temp%materials(1:mats)
          Deallocate (user_temp%materials)

          ! Add new material data to resized user%materials
          user%materials(mats+1) = loadMatProps(cmname, nprops, props)
          m = user%materials(mats+1)  ! Abbreviate

          ! Calculate alpha0
          If (m%matrixDam) Then
            m%alpha0_deg = NINT(m%alpha0*45.d0/ATAN(one))
            m%alpha0 = alpha0_DGD(m)
          End If

          ! Run consistency checks
          If (totalTime == zero) Then
            Call consistencyChecks(m, issueWarnings=.TRUE.)
          Else
            Call consistencyChecks(m, issueWarnings=.FALSE.)
          End If

          ! Make sure updated values in m are retained
          user%materials(mats+1) = m

        End If

      End Do readMaterial

    Else

      Call log%info("Saving initial user material: "//trim(cmname))
      Allocate (user%materials(1))

      ! Load the first material
      user%materials(1) = loadMatProps(cmname, nprops, props)
      m = user%materials(1)  ! Abbreviate

      ! Calculate alpha0
      If (m%matrixDam) Then
        m%alpha0_deg = NINT(m%alpha0*45.d0/ATAN(one))
        m%alpha0 = alpha0_DGD(m)
      End If

      ! Run consistency checks
      If (totalTime == zero) Then
        Call consistencyChecks(m, issueWarnings=.TRUE.)
      Else
        Call consistencyChecks(m, issueWarnings=.FALSE.)
      End If

      ! Make sure updated values in m are retained
      user%materials(1) = m

    End If

  Else  ! After user%materials(:) has been fully populated
    Do I = 1,Size(user%materials)
      If (user%materials(I)%name == trim(cmname))  m = user%materials(I)
    End Do
  End If

  ! This forces the material to behave elastically during the packager and initial increment. This
  ! was motivated by issues with Intel MKL in the packager.
  If (totalTime == 0) Then
    m%matrixDam = .FALSE.
    m%shearNonlinearity12 = .FALSE.
    m%shearNonlinearity13 = .FALSE.
    m%schapery = .FALSE.
    m%schaefer = .FALSE.
    m%fiberTenDam = .FALSE.
    m%fiberCompDamBL = .FALSE.
    m%fiberCompDamFKT12 = .FALSE.
    m%fiberCompDamFKT13 = .FALSE.
    m%friction = .FALSE.
    fatigue_step = .FALSE.
  Else
    fatigue_step = .FALSE.
#ifndef PYEXT
    ! Is this a fatigue step?
    ptr_fatigue_int = SMAIntArrayAccess(1)
    If (fatigue_parameters(1) == 1) Then
      fatigue_step = .TRUE.
      Call log%debug("This is a fatigue step.")
    End If
#endif
  End If

  ! -------------------------------------------------------------------- !
  master: Do km = 1,nblock  ! Master Loop

  ! Update location information for debugging
  Call log%arg%update(nElement(km), nMatPoint, nLayer, nSecPoint)

  Call log%debug("Master loop. km = " // trim(str(km)))

  ! -------------------------------------------------------------------- !
  !    Recall previous elastic limits, plasticity, and damage states:    !
  ! -------------------------------------------------------------------- !
  sv = loadStateVars(nstatev, stateOld(km,:), m)

  ! -------------------------------------------------------------------- !
  !    Cohesive elements:                                                !
  ! -------------------------------------------------------------------- !
  ElementType: If (m%cohesive) Then

    ! Determine current cohesive displacement-jump
    !  Assumes a constitutive thickness equal to one.
    If (nshr == 2) Then
      delta(1) = sv%Fb1 + two*strainInc(km,3)  ! Longitudinal shear direction (1--3)
      delta(2) = sv%Fb2 + strainInc(km,1)      ! Normal direction (3--3)
      delta(3) = sv%Fb3 + two*strainInc(km,2)  ! Transverse shear direction (2--3)
    Else
      delta(1) = sv%Fb1 + two*strainInc(km,2)  ! Shear direction (1--2)
      delta(2) = sv%Fb2 + strainInc(km,1)      ! Normal normal (2--2)
      delta(3) = zero
    End If

    ! Define penalty stiffnesses
    ! Normal (Mode I) penalty stiffness. This is a user-input numerical constant.
    Pen(2) = m%E3
    ! Longitudinal (Pen(1)) and transverse (Pen(3)) shear penalty stiffnesses
    !  Defined according the material property relationships defined in Turon (2010). Also includes the
    !  normal compression load dependence of the shear strengths from LaRC04.
    Pen(1) = Pen(2)*m%GYT*(m%SL - m%etaL*MAX(-m%YC, Pen(2)*MIN(zero, delta(2))))**2/(m%GSL*m%YT**2)
    Pen(3) = Pen(2)*m%GYT*(m%ST - m%etaT*MAX(-m%YC, Pen(2)*MIN(zero, delta(2))))**2/(m%GSL*m%YT**2)

    ! Store current cohesive displacement-jump
    sv%Fb1 = delta(1)  ! Longitudinal shear displacement-jump
    sv%Fb2 = delta(2)  ! Normal displacement-jump
    sv%Fb3 = delta(3)  ! Transverse shear displacement-jump

    If (totalTime == 0) Then  ! Avoids calculating damage during the packager and the first solution increment
      Call cohesive_damage(m, p, delta, Pen, delta(2), sv%B, sv%FIm, .FALSE.)
      dmg_penalty = zero
      dGdGc = zero
    Else
      Call cohesive_damage(m, p, delta, Pen, delta(2), sv%B, sv%FIm, fatigue_step, sv%d2, dmg_penalty, dGdGc)
    End If

    If (m%friction .AND. delta(2) <= zero) Then  ! Closed cracks with friction
      Sliding = crack_is_sliding(delta, Pen, sv%slide, m%mu, m%mu)
      Call crack_traction_and_slip(delta, Pen, sv%slide, sv%slide, m%mu, m%mu, dmg_penalty, sv%d2, T_coh, Sliding)
    Else  ! Closed cracks without friction and open cracks
      T_coh = cohesive_traction(delta, Pen, dmg_penalty)
      sv%slide(1) = delta(1)
      sv%slide(2) = delta(3)
    End If

    ! Allows for cohesive elements to be deleted when their damage variable reaches one
    If (sv%d2 >= one) sv%STATUS = 0

    ! Update stress values
    stressNew(km,1) = T_coh(2)  ! Normal stress
    If (nshr == 2) Then  ! 3-D
      stressNew(km,3) = T_coh(1)  ! Longitudinal shear stress (1--3)
      stressNew(km,2) = T_coh(3)  ! Transverse shear stress (2--3)
    Else  ! 2-D
      stressNew(km,2) = T_coh(1)  ! Shear stress
    End If

    ! Update internal and inelastic energy terms
    enerInternNew(km) = DOT_PRODUCT(T_coh, delta)/two
    enerInternNew(km) = enerInternNew(km)/density_abq(km)
    enerInelasNew(km) = dGdGc*(m%GYT + (m%GSL - m%GYT)*sv%B**m%eta_BK)
    enerInelasNew(km) = enerInelasOld(km) + enerInelasNew(km)/density_abq(km)

  ! -------------------------------------------------------------------- !
  !    Solid elements:                                                   !
  ! -------------------------------------------------------------------- !
  Else ElementType

    ! -------------------------------------------------------------------- !
    !    Deformation Gradient Tensor and Right Stretch Tensor              !
    ! -------------------------------------------------------------------- !
    U = Vec2Matrix(stretchNew(km,:))
    F = Vec2Matrix(defgradNew(km,:))
    F_old = Vec2Matrix(defgradOld(km,:))

    ! -------------------------------------------------------------------- !
    ! As of Abaqus 6.16, the packager receives a defGradNew of (0.999, 0.999, 0.0, 0.001, 0.001)
    ! for S4R elements. The F(3,3) of 0.0 breaks the initial pass through the VUMAT and the model
    ! will not run. The following statement is a workaround to this problem.
    If (totalTime == 0 .AND. nshr == 1) F(3,3) = one

    ! -------------------------------------------------------------------- !
    !    Define the characteristic element lengths                         !
    ! -------------------------------------------------------------------- !
    If (sv%Lc(1) == zero) Then

      sv%Lc(1) = charLength(km, 1)
      sv%Lc(2) = charLength(km, 2)
      If (nshr == 1) Then
        sv%Lc(3) = m%thickness
      Else
        sv%Lc(3) = charLength(km, 3)
      End If

      ! Check for strange values of Lc2 (may indicate that vucharlength did not run)
      If (totalTime > 0) Then ! Don't run in packager
        If (sv%Lc(2) .EQ. zero) Then
          Call log%error("Found Lc2 = 0. Perhaps *Characteristic Length is missing from the input deck?")
        Else If (sv%Lc(3) .EQ. zero) Then
          Call log%error("Found Lc3 = 0. Perhaps *Characteristic Length is missing from the input deck?")
        ! Else If (sv%Lc(2) > sv%Lc(1)*10.d0 .OR. sv%Lc(2) < sv%Lc(1)/10.d0) Then
        !   Call log%warn("Found Lc = [" // trim(str(sv%Lc(1))) //','// trim(str(sv%Lc(2))) //','// trim(str(sv%Lc(3))) // "]; perhaps *Characteristic Length is missing from the input deck?")
        End If
      End If

      Call log%debug("Characteristic element lengths:")
      Call log%debug(trim(str(sv%Lc(1)))//' '//trim(str(sv%Lc(2)))//' '//trim(str(sv%Lc(3))))

    End If

    If (totalTime <= dt) Then
      ! Call checkForSnapBack(m, sv%Lc, nElement(km))
      sv%debugpy_count = 0  ! Initialize
    End If

    ! -------------------------------------------------------------------- !
    !    Initialize phi0                                                   !
    ! -------------------------------------------------------------------- !
    If (totalTime <= dt .AND. (m%fiberCompDamFKT12 .OR. m%fiberCompDamFKT13)) Then
      Call initializePhi0(m, sv%Lc, charLength(km, 4:6), sv%phi0_12, sv%phi0_13)
      Call log%debug("Calculated initial misalignments, phi0_12: " // trim(str(sv%phi0_12)) // " and phi0_13: " // trim(str(sv%phi0_13)))
    End If

    ! -------------------------------------------------------------------- !
    !    Initialize inelastic energy to the previous value                 !
    ! -------------------------------------------------------------------- !
    enerInelasNew(km) = enerInelasOld(km)

    ! -------------------------------------------------------------------- !
    !    Initialize fiber failure                                          !
    ! -------------------------------------------------------------------- !
    If (m%fiberCompDamFKT12 .AND. p%fkt_fiber_failure_angle > zero) Then
      sv%Inel12c = intializeFiberFailure(sv%phi0_12, p%fkt_fiber_failure_angle, m%G12, m%aPL, m%nPL)
    Else
      sv%Inel12c = Huge(zero)   ! Turn off fiber failure by setting the associate inelastic strain to a very large number
    End If
    If (m%fiberCompDamFKT13 .AND. p%fkt_fiber_failure_angle > zero) Then
      sv%Inel13c = intializeFiberFailure(sv%phi0_13, p%fkt_fiber_failure_angle, m%G13, m%aPL, m%nPL)
    Else
      sv%Inel13c = Huge(zero)   ! Turn off fiber failure by setting the associate inelastic strain to a very large number
    End If

    ! -------------------------------------------------------------------- !
    !    Damage Calculations:                                              !
    ! -------------------------------------------------------------------- !
    ! Damage initiation prediction
    If (.NOT. (m%matrixDam .AND. sv%d2 > zero) .AND. .NOT. ((m%fiberCompDamFKT12 .OR. m%fiberCompDamFKT13) .AND. sv%d1C > zero)) Then

      Call DGDInit(U,F,F_old,m,p,sv,ndir,nshr,tempNew(km),density_abq(km),Cauchy,enerInternNew(km),enerInelasNew(km),fatigue_step)

    End If

    ! Matrix crack damage evolution
    If (m%matrixDam .AND. sv%d2 > zero) Then

      Call DGDEvolve(U,F,F_old,m,p,sv,ndir,nshr,tempNew(km),density_abq(km),Cauchy,enerInternNew(km),enerInelasNew(km),fatigue_step)

    ! Fiber compression damage evolution (FKT decomposition)
    Else If ((m%fiberCompDamFKT12 .OR. m%fiberCompDamFKT13) .AND. sv%d1C > zero) Then

      Call DGDKinkband(U,F,F_old,m,p,sv,ndir,nshr,tempNew(km),density_abq(km),Cauchy,enerInternNew(km),enerInelasNew(km))

    End If

    ! -------------------------------------------------------------------- !
    !    Determine and store the stress tensor:                            !
    ! -------------------------------------------------------------------- !
    ! Rotation tensor
    Rot = MATMUL(F, MInverse(U))
    ! Cauchy stress in the current configuration
    CauchyABQ = MATMUL(TRANSPOSE(Rot), MATMUL(Cauchy, Rot))
    ! Convert to vector format
    stressNew(km,:) = Matrix2Vec(CauchyABQ, nshr)

    ! -------------------------------------------------------------------- !
    !    Excessive shear strain errors                                      !
    ! -------------------------------------------------------------------- !
    If (m%shearNonlinearity12) Then
      If (sv%Inel12 < stateOld(km,13)) Then
        Call writeDGDArgsToFile(m,p,sv,U,F,F_old,ndir,nshr,tempNew(km),density_abq(km),log%arg,'CompDam_DGD')
        Call log%terminate('Decrease in inelastic strain 12.')
      End If
    End If
    If (m%shearNonlinearity13) Then
      If (sv%Inel13 < stateOld(km,21)) Then
        Call writeDGDArgsToFile(m,p,sv,U,F,F_old,ndir,nshr,tempNew(km),density_abq(km),log%arg,'CompDam_DGD')
        Call log%terminate('Decrease in inelastic strain 13.')
      End If
    End If
    If (m%shearNonlinearity12 .OR. m%shearNonlinearity13) Then
      If (sv%Inel12 > 1.d0 .OR. sv%Inel13 > 1.d0) Then
        Call writeDGDArgsToFile(m,p,sv,U,F,F_old,ndir,nshr,tempNew(km),density_abq(km),log%arg,'CompDam_DGD')
        Call log%terminate('Excessive inelastic shear strain.')
      End If
    End If

  End If ElementType

  ! -------------------------------------------------------------------- !
  !    Store the updated state variables:                                !
  ! -------------------------------------------------------------------- !
  stateNew(km,:) = storeStateVars(sv, nstatev, m)

  ! -------------------------------------------------------------------- !
  !    Kill job for debugging purposes                                   !
  ! -------------------------------------------------------------------- !
  If (p%logLevel > 2 .AND. p%debug_kill_at_total_time > zero) Then
    If (log%arg%totalTime >= p%debug_kill_at_total_time) Then
      Call writeDGDArgsToFile(m,p,sv,U,F,F_old,ndir,nshr,tempNew(km),density_abq(km),log%arg,'CompDam_DGD')
      Call log%error('debug_kill_at_total_time condition satisfied; terminating job.')
    End If
  End If

  End Do master ! End Master Loop

  Call log%debug('End of VUMAT')

  Return
End Subroutine CompDam

ext: for


text:
program main

   use omp_lib

   integer :: id, nt

   !$OMP PARALLEL private(id, nt)
   id = omp_get_thread_num()
   nt = omp_get_num_threads()
   print *, "Hello OpenMP", id, ' of ', nt
   !$OMP END PARALLEL

   print *, omp_get_wtime()

end program

max_stars_repo_name: WVUHPC/Modern-Fortran

max_issues_repo_name: WVUHPC/Modern-Fortran

max_forks_repo_name: WVUHPC/Modern-Fortran


text:
FUNCTION FPROBDEP(TAUDDD,TAUBBB)

  !**********************************************************************C
  !
  ! **  FPROBDEP CALCULATES PROBABILITY OF DEPOSITION USING PROBABILITY
  ! **  INTEGRAL
  !
  ! **  THIS SUBROUTINE IS PART OF  EFDC-FULL VERSION 1.0a
  !
  ! **  LAST MODIFIED BY JOHN HAMRICK ON 1 NOVEMBER 2001
  !
  !----------------------------------------------------------------------C
  !
  ! CHANGE RECORD
  ! DATE MODIFIED     BY                 DATE APPROVED    BY
  !
  !----------------------------------------------------------------------C
  !
  !**********************************************************************C
  !

  IMPLICIT NONE

  INTEGER :: INEG
  REAL   :: FPROBDEP,YVAL,XVAL,POLYX,EXPY,FUNY,TMPVAL
  REAL, INTENT(IN) :: TAUDDD,TAUBBB

  ! *** EVALUATION ASSUMES TAUBBB > TAUDDD
  YVAL=2.04*LOG(0.25*((TAUBBB/TAUDDD)-1.)*EXP(1.27*TAUDDD))
  INEG=0
  IF( YVAL < 0.0 )THEN
    INEG=1
    YVAL=ABS(YVAL)
  ENDIF
  XVAL=1.0/(1.0+0.3327*YVAL)
  POLYX=XVAL*(0.4632-0.1202*XVAL+0.9373*XVAL*XVAL)
  EXPY=-0.5*YVAL*YVAL
  FUNY=0.3989*EXP(EXPY)
  TMPVAL=1.0-FUNY*POLYX
  IF( INEG == 1 ) TMPVAL=1.0-TMPVAL
  FPROBDEP=1.0-TMPVAL

  ! **  NOTES
  !     0.3989=1/SQRT(2*PI)

END FUNCTION

max_stars_repo_name: dsi-llc/EFDCPlus

max_forks_repo_name: dsi-llc/EFDCPlus


text:
      subroutine gtnods

      include "common.h"
      include "mpif.h"
      include "auxmpi.h"

!SCATTER      dimension irecvcount(numpe), numvec(numpe)

      if(numpe > 1) then
         irecvcount = 1
         numvec = nshg0
         call MPI_ALLREDUCE (numvec, nshgt, irecvcount,
     &                 MPI_INTEGER, MPI_SUM, MPI_COMM_WORLD, ierr)
c         call MPI_REDUCE_SCATTER (numvec, nshgt, irecvcount,
c     &                 MPI_INTEGER, MPI_SUM, MPI_COMM_WORLD, ierr)
      else
         nshgt = nshg0
      endif

      if (myrank .eq. master) then
         write(6,*) 'Total number of modes = ',nshgt
      endif

      return
      end