!--------------------------------------------------------------------------------------------------
! MODULE: overlap_input_mod
!> @author Marin Sapunar, Ruđer Bošković Institute
!> @date August, 2018
!
!> @brief Hold subroutines for overlap program input handling.
!--------------------------------------------------------------------------------------------------
module overlap_input_mod
    use global_defs
    use overlap_variables
    implicit none

    private
    public :: print_help
    public :: set_defaults
    public :: command_line_interface

    character(len=4), parameter :: default_file_ao = 's_ao'
    character(len=4), parameter :: default_file_mo = 's_mo'
    character(len=4), parameter :: default_file_wf = 's_wf'

contains

    !----------------------------------------------------------------------------------------------
    ! SUBROUTINE: print_help
    !> @brief Print help message and exit program.
    !----------------------------------------------------------------------------------------------
    subroutine print_help()
        write(stdout, '(a)') 'usage: cis_overlap.exe [optional arguments] path1 path2'
        write(stdout, '(a)')
        write(stdout, '(a)') 'Calculate overlaps between two sets CIS type wave functions.'
        write(stdout, '(a)')
        write(stdout, '(a)') 'positional arguments:'
        write(stdout, '(a)') '  path1                       directory containing calculation for bra states      '
        write(stdout, '(a)') '  path2                       directory containing calculation for ket states      '
        write(stdout, '(a)')
        write(stdout, '(a)') 'optional arguments:'
        write(stdout, '(a)') '  -h, --help                  show this help message and exit                      '
        write(stdout, '(a)') '  -in, --input-format         format of input files                                '
        write(stdout, '(32x,a,a)') 'default: ', input_format_1
        write(stdout, '(a)') '  -in1, --input-format-1      format of input files for bra orbitals/states        '
        write(stdout, '(32x,a,a)') 'default: ', input_format_1
        write(stdout, '(a)') '  -in2, --input-format-2      format of input files for ket orbitals/states        '
        write(stdout, '(32x,a,a)') 'default: ', input_format_2
        write(stdout, '(a)') '  -ao, --ao-stop              stop program after AO overlap calculation            '
        write(stdout, '(a)') '  -mo, --mo-stop              stop program after MO overlap calculation            '
        write(stdout, '(a)') '  -alg, --algorithm ALG       algorithm to use for the overlap calculation         '
        write(stdout, '(a)') '                              Available options:                                   '
        write(stdout, '(a)') '                                CIS (very slow, not recommended)                   '
        write(stdout, '(a)') '                                L2M                                                '
        write(stdout, '(a)') '                                NTO                                                '
        write(stdout, '(32x,a,a)') 'default: ', cis_algorithm
        write(stdout, '(a)') '  -dys, --dyson               calculate Dyson orbitals instead of overlaps         '
        write(stdout, '(32x,a,l1)') 'default: ', dyson_c
        write(stdout, '(a)') '  -t, --wf-threshold t        truncate wave functions using given threshold        '
        write(stdout, '(a)') '  -tnex, --truncate-nex n     truncate wave functions to n dominant excitations    '
        write(stdout, '(a)') '  -ns, --(no-)norm-states     renormalize input states before calculation          '
        write(stdout, '(32x,a,l1)') 'default: ', norm_states
        write(stdout, '(a)') '  -os, --(no-)orth-states     reorthogonalize input states before calculation      '
        write(stdout, '(32x,a,l1)') 'default: ', orth_states
        write(stdout, '(a)') '  -oo, --(no-)orth-overlap    orthogonalize overlap matrix                         '
        write(stdout, '(32x,a,l1)') 'default: ', orth_overlap
        write(stdout, '(a)') '  -mp, --(no-)match-phase     match phase between assigned bra/ket states          '
        write(stdout, '(32x,a,l1)') 'default: ', match_phase
        write(stdout, '(a)') '  -op, --operator op          operator whose expectation values are cacluclated.   '
        write(stdout, '(a)') '                              Currently, only overlaps can be calculated for wave  '
        write(stdout, '(a)') '                              functions.                                           '
        write(stdout, '(a)') '                              Available options:                                   '
        write(stdout, '(a)') '                                analytic expression (eg. "x * y^2 - z^2")          '
        write(stdout, '(a)') '                                  note: use quotation marks                        '
        write(stdout, '(a)') '                                  note: see read_operator_string subroutine docs   '
        write(stdout, '(a)') '                                overlap                                            '
        write(stdout, '(a)') '                                r^2 (equivalent to "x^2 + y^2 + z^2")              '
        write(stdout, '(32x,a,a)') 'default: ', operator_string
        write(stdout, '(a)') '  -cp, --(no-)center-pairs    attempt to remove effect of basis set translation by '
        write(stdout, '(a)') '                              recentering pairs of AOs in AO overlap calculation   '
        write(stdout, '(a)') '                              (untested, not recommended)                          '
        write(stdout, '(32x,a,l1)') 'default: ', center_pairs
        write(stdout, '(a)') '  -fmn, --freeze-mo-norm t    freeze occupied ket MOs when their norm in bra basis '
        write(stdout, '(a)') '                              is smaller than given threshold. Same number of bra  '
        write(stdout, '(a)') '                              MOs with smallest norms in ket basis is also frozen. '
        write(stdout, '(a)') '                              Used when geometry of a small part of a system is    '
        write(stdout, '(a)') '                              significantly different between bra and ket states.  '
        write(stdout, '(a)') '                              (untested, not recommended)                          '
        write(stdout, '(32x,a,l1)') 'default: ', freeze_mo_norm
        write(stdout, '(a)') '  -sao, --sao-output file     output final ao overlap matrix to file               '
        write(stdout, '(32x,a,a)') 'default: ',  outfile_ao
        write(stdout, '(a)') '  -smo, --smo-output file     output final mo overlap matrix to file               '
        write(stdout, '(32x,a,a)') 'default: ',  outfile_mo
        write(stdout, '(a)') '  -swf, --swf-output file     output final wf overlap matrix to file               '
        write(stdout, '(32x,a,a)') 'default: ',  outfile_wf
        write(stdout, '(a)') '  -p, --print-level p         control output level of program (0 = quiet)          '
        write(stdout, '(32x,a,i0)') 'default: ', print_level
        stop
    end subroutine print_help

    !----------------------------------------------------------------------------------------------
    ! SUBROUTINE: set_defaults
    !> @brief Set default values for program options.
    !----------------------------------------------------------------------------------------------
    subroutine set_defaults()
        use ang_mom_defs

        call amp%init(7)
        ! CLI options.
        input_format_1 = 'turbomole'
        input_format_2 = 'turbomole'
        ao_stop = .false.
        mo_stop = .false.
        cis_algorithm = 'NTO'
        dyson_c = .false.
        wf_threshold = 1.0_dp
        truncate_nex = 0
        norm_states = .true.
        orth_states = .false.
        orth_overlap = .false.
        match_phase = .false.
        center_atoms = .false. !unused
        center_pairs = .false.
        freeze_mo_norm = .false.
        freeze_mo_norm_t = -1.0_dp
        operator_string = 'overlap'
        outfile_ao = 'None'
        outfile_mo = 'None'
        outfile_wf = default_file_wf
        prefix_dyson = 'dys'
        print_level = 2
        ! Set timings to 0
        time_in = 0.0_dp
        time_ao = 0.0_dp
        time_mo = 0.0_dp
        time_wf = 0.0_dp
        time_tot = 0.0_dp
    end subroutine set_defaults

    !----------------------------------------------------------------------------------------------
    ! SUBROUTINE: command_line_interface
    !> @brief Read command line options.
    !----------------------------------------------------------------------------------------------
    subroutine command_line_interface()
        integer :: narg, i
        character(len=1000) :: temp
        narg = command_argument_count()
        i = 0
        if (narg == 0) call print_help()
        do while (i < narg)
            i = i + 1
            if (i == narg) call print_help()
            call get_command_argument(i, temp)
            if ((i == narg - 1) .and. (temp /= '--help') .and. (temp /= '-h')) exit
            select case(temp)
            case('--help', '-h')
                call print_help()
            case('--input-format', '-in')
                i = i + 1
                call get_command_argument(i, temp)
                input_format_1 = trim(adjustl(temp))
                input_format_2 = trim(adjustl(temp))
            case('--input-format-1', '-in1')
                i = i + 1
                call get_command_argument(i, temp)
                input_format_1 = trim(adjustl(temp))
            case('--input-format-2', '-in2')
                i = i + 1
                call get_command_argument(i, temp)
                input_format_2 = trim(adjustl(temp))
            case('--ao-stop', '-ao')
                ao_stop = .true.
                if (outfile_ao == 'None') outfile_ao = default_file_ao
            case('--mo-stop', '-mo')
                mo_stop = .true.
                if (outfile_mo == 'None') outfile_mo = default_file_mo
            case('--algorithm', '-alg')
                i = i + 1
                call get_command_argument(i, temp)
                cis_algorithm = temp(1:4)
            case('--dyson', '-dys')
                dyson_c = .true.
            case('--no-dyson', '-ndys')
                dyson_c = .false.
            case('--wf-threshold', '-t')
                i = i + 1
                call get_command_argument(i, temp)
                read(temp, *) wf_threshold
            case('--truncate-nex', '-tnex')
                i = i + 1
                call get_command_argument(i, temp)
                read(temp, *) truncate_nex
            case('--norm-states', '-ns')
                norm_states = .true.
            case('--no-norm-states', '-nns')
                norm_states = .false.
            case('--orth-states', '-os')
                orth_states = .true.
            case('--no-orth-states', '-nos')
                orth_states = .false.
            case('--orth-overlap', '-oo')
                orth_overlap = .true.
            case('--no-orth-overlap', '-noo')
                orth_overlap = .false.
            case('--match-phase', '-mp')
                match_phase = .true.
            case('--no-match-phase', '-nmp')
                match_phase = .false.
            case('--operator', '-op')
                i = i + 1
                call get_command_argument(i, temp)
                operator_string = trim(adjustl(temp))
            case('--center-pairs', '-cp')
                center_pairs = .true.
            case('--no-center-pairs', '-ncp')
                center_pairs = .false.
            case('--freeze-mo-norm', '-fmn')
                i = i + 1
                call get_command_argument(i, temp)
                read(temp, *) freeze_mo_norm_t
                if (freeze_mo_norm_t > 0.0_dp) freeze_mo_norm = .true.
            case('--sao-output', '-sao')
                i = i + 1
                call get_command_argument(i, temp)
                outfile_ao = trim(adjustl(temp))
            case('--smo-output', '-smo')
                i = i + 1
                call get_command_argument(i, temp)
                outfile_mo = trim(adjustl(temp))
            case('--swf-output', '-swf')
                i = i + 1
                call get_command_argument(i, temp)
                outfile_wf = trim(adjustl(temp))
            case('--print-level', '-p')
                i = i + 1
                call get_command_argument(i, temp)
                read(temp, *) print_level
            end select
        end do
        call get_command_argument(narg-1, temp)
        path1 = trim(adjustl(temp))
        call get_command_argument(narg, temp)
        path2 = trim(adjustl(temp))

        if (operator_string /= 'overlap') then
            if (.not. mo_stop) then
                write(stderr, *) 'Warning. Only overlap operator is implemented for wave functions.'
                write(stderr, *) '  Setting mo_stop option to .true.'
                mo_stop = .true.
                if (outfile_mo == 'None') outfile_mo = default_file_mo
            end if
        end if
    end subroutine command_line_interface

end module overlap_input_mod


text:
      program prog
      integer size, dim
      parameter(size=3, dim=2)
      character strings*(size)(dim)
      strings(1) = 'AbC'
      strings(2) = 'dEf'
      write(6, '(i1,x,i1,x,a,x,a)') size, dim, strings
      end

ext: f


text:
module T { }
type T

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:
      FUNCTION K32R(J,M,K,L,R, JA)
C BY CHARLES KLUEPFEL
C EDITED FROM BASIC BY AX, 1.19.79
      k32r=ja    ! Added 7/27/93 to report own old value
      IF(M.GT.1) GO TO 520
C # OF HIS COOPS AFTER MY DEF.
      C1=0
C # OF HIS DEFECTIONS AFTER MY DEF.
      C2=0
C # OF HIS COOPS AFTER MY COOPERATION
      C3=0
C # OF HIS DEF. AFTER MY COOPERATION
      C4=0
C HIS 3RD PREV. CHOICE
      J2=0
C HIS 2ND PREV. CHOICE
      J1=0
C MY 2ND PREV. CHOICE
      I2=0
C MY PREV. CHOICE
      I1=0
C PROB. OF MY RESPONDING IN KIND
      P=0
      GO TO 550
  520 IF(M.EQ.2) GO TO 550
C RECORD HIS RESPONSES
      IF(I2.EQ.0) GO TO 530
      IF(J.EQ.0) GO TO 528
      C2=C2+1
      GO TO 540
  528 C1=C1+1
      GO TO 540
  530 IF(J.EQ.0) GO TO 537
      C4=C4+1
      GO TO 540
  537 C3=C3+1
C AFTER M. 26 TRY DETECTING RNDOM
  540 IF(M.LT.27) GO TO 550
      IF (C1.LT.((C1+C2) - 1.5*SQRT(C1+C2)) / 2 ) GO TO 550
      IF (C4.LT.((C3+C4) - 1.5*SQRT(C3+C4)) / 2 ) GO TO 550
      K32R=1
      GO TO 590
  550 K32R=0
      IF(J1.NE.J) GO TO 570
      IF(J2.NE.J1) GO TO 580
C RESPOND IN KIND TO 3 IN ROW.
      K32R = J
      GO TO 590
C PROB .6 OF GIVING DEF. AFTER SINGLE DEF.
  570 P=.6
      IF(J.EQ.1) GO TO 585
C PROB .7 GIVING COOP AFTER SINGLE COOP.
      P=.7
      GO TO 585
C PROB .9 RESP. IN KIND TO 2 IN ROW
  580 P=.9
  585 K32R=J
      IF(R.LT.P) GO TO 590
      K32R=1-J
C PUSH BACK CHOICES
  590 J2=J1
      J1=J
      I2=I1
      I1=K32R
      RETURN
      END

max_stars_repo_name: Axelrod-Python/TourExec

max_issues_repo_name: Axelrod-Python/TourExec

max_forks_repo_name: Axelrod-Python/TourExec


text:
module Mod_plcd_DamageMaterial
   use typre
   use Mod_plcd_Material
   use Mod_Listen
   use Mod_Element
   use Mod_plcd_ComparisonCriteria
   use Mod_plcd_ComparisonCriteriaFactory
   implicit none
   private
   public DamageMaterial

   type, extends(PLCDMaterial) :: DamageMaterial

      real(rp) :: DamageLimit = 0.9999
      real(rp) :: Tau0, FractureEnergy,  UniaxialStressLimit
      character(5) :: CompareCriteria = 'SIMOJ'
      logical :: DeviatoricOnlyDamage = .false.

      type(r1p), allocatable :: DPostprocess(:),TauPostprocess(:),CompPostprocess(:)

      logical :: isReadComparisonCriteria = .false.
      class(ComparisonCriteria), pointer :: ComparisonCriteria
      procedure(ComputeDamage), pointer :: ComputeDamage => NULL()
      procedure(ComputeCT), pointer :: ComputeCT => NULL()

contains
      procedure :: SpecificCreateElementData => MD_SpecificCreateElementData
      procedure :: SpecificReadData => MD_SpecificReadData
      procedure :: SpecificScatterData => MD_SpecificScatterData
      procedure :: SpecificSetup       => MD_SpecificSetup

      !Postprocess
      procedure :: SetupPostprocess    => MD_SetupPostprocess
      procedure :: DoPostprocess       => MD_DoPostprocess
      procedure :: IsElementSizeRequiredByEMDs => MD_IsElementSizeRequiredByEMDs

   end type

   type :: DamageMaterial_GPD
      logical :: IsDamaged = .false.
      real(rp) :: ReferenceTau, CurrentTau
      real(rp) :: ReferenceD = 0.0_rp, CurrentD = 0.0_rp
      real(rp), pointer :: C(:,:) => NULL()
      real(rp), pointer :: S(:) => NULL()
      real(rp) :: ComparisonRatio = 0.0_rp
   end type

   type, extends(ElementMaterialData) :: DamageMaterial_ED

      logical :: JustCreated = .false. !Required for adaptivity
      class(DamageMaterial), pointer :: Material => NULL()

      real(rp) :: ElementSize = 1.0_rp
      type(DamageMaterial_GPD), allocatable :: GP(:)

      !For UP elements
      real(rp) :: GpPressure

contains

      procedure :: ComputeHistoryAndConstitutiveTensor => MD_ComputeHistoryAndConstitutiveTensor
      procedure :: ComputeStress => MD_ComputeStress
      procedure :: GetConstitutiveTensorPointer => MD_GetConstitutiveTensorPointer
      procedure :: GetStressTensorPointer => MD_GetStressTensorPointer
      procedure :: SetElementSize => MD_SetElementSize
      procedure :: MoveHistoryVariablesToConverged => MD_MoveHistoryVariablesToConverged

      !Postprocess
      procedure :: ContributeToPostprocess => MD_ContributeToPostprocess
      procedure :: GetMaterialPointer
      procedure :: Finalize

      !Adaptivity
      procedure :: CopyTo
      procedure :: SpecificHistoryAddFromGaussArray
      procedure :: SpecificHistoryFinalRefinementUpdate
      procedure :: SpecificHistoryGetToGaussArray
      procedure :: GetHistoryVariablesNdofn
      procedure :: GetComparisonCriteriaRatio

      procedure :: GetBufferSize
      procedure :: ToBuffer
      procedure :: InitializeFromBuffer

      !For u-p elements
      procedure :: GetInverseVolumetricDeformationModulus
      procedure :: GetSecantCNorm
      procedure :: SetPressureForComputeHistoryAndConstitutiveTensor

   end type

   interface
      Subroutine ComputeDamage(a,ElementSize,GP)
         import
         implicit none
         class(DamageMaterial) :: a
         real(rp), intent(in) :: ElementSize
         class(DamageMaterial_GPD) :: GP
      end subroutine

      Subroutine ComputeCT(a,ElementSize,ElasticPredictor,GP)
         import
         implicit none
         class(DamageMaterial) :: a
         real(rp), intent(in) :: ElementSize
         real(rp), Intent(in) :: ElasticPredictor(:)
         class(DamageMaterial_GPD) :: GP
      end subroutine

   end interface

contains

   subroutine MD_SpecificCreateElementData(a,pgaus,ElementData)
      implicit none
      class(DamageMaterial), target :: a
      integer(ip) :: pgaus
      class(ElementMaterialData), pointer :: ElementData

      type(DamageMaterial_ED), pointer :: MD_ElementData
      type(DamageMaterial_GPD), pointer :: GP(:)
      integer(ip) :: igaus

      allocate(MD_ElementData)
      MD_ElementData%JustCreated = .true.
      MD_ElementData%Material => a
      allocate(MD_ElementData%GP(pgaus))

      GP => MD_ElementData%GP

      do igaus = 1,pgaus
         GP(igaus)%ReferenceTau = a%Tau0
         GP(igaus)%CurrentTau = a%Tau0
         GP(igaus)%C => a%C
         allocate(GP(igaus)%S(a%CT%VoigtSize))
      enddo

      ElementData => MD_ElementData
   end subroutine

   subroutine MD_SpecificReadData(a,Listener)
      class(DamageMaterial) :: a
      type(ListenFile) :: Listener

      if (Listener%words(1) == 'UNIAX') then
         a%UniaxialStressLimit = Listener%param(1)
      elseif (Listener%words(1) == 'FRACT') then
         a%FractureEnergy = Listener%param(1)
      elseif (Listener%words(1) == 'MAXDA' ) then
         a%DamageLimit = Listener%param(1)
      elseif (Listener%words(1) == 'COMPA' ) then
         !ComparisonCriteria
         a%CompareCriteria = Listener%words(2)
         call ComparisonCriteriaFactory(a%CompareCriteria,a%ComparisonCriteria)
         call a%ComparisonCriteria%ReadData(Listener)
         call AllocateDamageAndTypeCT(a)
         a%isReadComparisonCriteria = .true.
      elseif (Listener%words(1) == 'TYPEO') then
         if(Listener%exists('DEVIA')) then
            a%DeviatoricOnlyDamage = .true.
         elseif(Listener%exists('FULL ')) then
            a%DeviatoricOnlyDamage = .false.
         endif
      endif

   end subroutine

   subroutine MD_SpecificScatterData(a)
      use MPI
      class(DamageMaterial) :: a

      integer(ip) :: ierr

      call MPI_BCAST(a%UniaxialStressLimit,1,MPI_REAL8,a%MPIroot,a%MPIcomm,ierr)
      call MPI_BCAST(a%FractureEnergy,1,MPI_REAL8,a%MPIroot,a%MPIcomm,ierr)
      call MPI_BCAST(a%DamageLimit,1,MPI_REAL8,a%MPIroot,a%MPIcomm,ierr)
      call MPI_BCAST(a%CompareCriteria,5,MPI_CHARACTER,a%MPIroot,a%MPIcomm,ierr)
      call MPI_BCAST(a%DeviatoricOnlyDamage,1,MPI_LOGICAL,a%MPIroot,a%MPIcomm,ierr)

      !Create the comparison Criteria
      if (a%MPIrank /= a%MPIroot .or. (a%IsReadComparisonCriteria .eqv. .false.)) then
         call ComparisonCriteriaFactory(a%CompareCriteria,a%ComparisonCriteria)
         call AllocateDamageAndTypeCT(a)
      endif
      call a%ComparisonCriteria%ScatterData

      !call a%ComparisonCriteria%ComputeTauValueFromUniaxialStressTest(a%UniaxialStressLimit,a%CT,a%Tau0)

   end subroutine

   subroutine MD_SpecificSetup(a)
      class(DamageMaterial) :: a

      call a%ComparisonCriteria%ComputeTauValueFromUniaxialStressTest(a%UniaxialStressLimit,a%CT,a%Tau0)

   end subroutine

   subroutine AllocateDamageAndTypeCT(a)
      class(DamageMaterial) :: a

      select case(a%TangentTensorType)
         case(0)
            a%ComputeCT => ComputeCT_Elastic

         case(1)
            a%ComputeCT => ComputeCT_Secant

         case(2)
            select case(a%CompareCriteria)
               case('SIMOJ')
                  a%ComputeCT => ComputeCT_Secant

               case('VONMI')
                  a%ComputeCT => ComputeCT_Secant !ComputeCT_Tangent_VonMises_Exp

            end select

         case default
            a%ComputeCT => ComputeCT_Secant

      end select

      select case(a%CompareCriteria)
         case('SIMOJ')
            a%ComputeDamage => ComputeDamage_Dev_Exp !ComputeDamage_Exp !ComputeDamage_SimoJu_Exp

         case('VONMI')
            a%ComputeDamage => ComputeDamage_Dev_Exp !ComputeDamage_Exp !ComputeDamage_VonMises_Exp

         case default
!            a%ComputeDamage => ComputeDamage_SimoJu_Exp

      end select

   end subroutine

   subroutine NULLSUB(a)
      class(DamageMaterial):: a
   end subroutine

   subroutine MD_ComputeHistoryAndConstitutiveTensor(a,igaus,gradDisp)
      use Mod_plcd_ComparisonCriteria

      class(DamageMaterial_ED), target :: a
      integer(ip) :: igaus
      real(rp) :: gradDisp(:,:)

      real(rp) :: strain(a%Material%CT%VoigtSize)
      real(rp) :: ElasticPredictor(a%Material%CT%VoigtSize)
      type(DamageMaterial_GPD), pointer :: GP => NULL()
      type(DamageMaterial), pointer :: Mat => NULL()

      GP => a%GP(igaus)
      Mat => a%Material

      call Mat%CT%GetStrain(gradDisp,strain)

      !Elastic predictor
      ElasticPredictor = matmul(Mat%C,strain)

      if (Mat%UseUPFormulation) call Mat%CT%AddSphericalComponent(-a%gpPressure,ElasticPredictor)

      !Not Converged Tau
      call Mat%ComparisonCriteria%Evaluate(ElasticPredictor,strain,Mat%CT,GP%CurrentTau)
      GP%ComparisonRatio = GP%CurrentTau/Mat%Tau0

      !Update damage parameter if necessary
      if (GP%CurrentTau > GP%ReferenceTau) then

         if (GP%isDamaged .eqv. .false.) then
            allocate(GP%C(Mat%CT%VoigtSize,Mat%CT%VoigtSize))
            GP%isDamaged = .true.
         endif

         call Mat%ComputeDamage(a%ElementSize,GP)
         call Mat%ComputeCT(a%ElementSize,ElasticPredictor,GP)

      else
         GP%CurrentTau = GP%ReferenceTau
         GP%CurrentD = GP%ReferenceD
      endif

  end subroutine

   subroutine MD_ComputeStress(a,igaus,gradDisp)
      class(DamageMaterial_ED), target :: a
      integer(ip) :: igaus
      real(rp) :: gradDisp(:,:)

      real(rp) :: strain(a%Material%CT%VoigtSize)
      type(DamageMaterial_GPD), pointer :: GP => NULL()
      real(rp) :: ElasticPredictor(a%Material%CT%VoigtSize)

      real(rp), pointer :: CDev(:,:),CSph(:,:)

      call a%Material%CT%GetStrain(gradDisp,strain)

      GP => a%GP(igaus)

      !Normal damage
      if (.not. a%Material%DeviatoricOnlyDamage) then
         !Elastic predictor
         ElasticPredictor = matmul(a%Material%C,strain)

         !Final Stress
         GP%S = ElasticPredictor*(1-GP%CurrentD)

      !DeviatoricOnlyDamage
      else
         call a%Material%CT%GetPointerToCDeviatoric(CDev)
         call a%Material%CT%GetPointerToCSpherical(CSph)

         GP%S = 0.0_rp

         !Spherical part, not damaged
         GP%S = matmul(CSph,strain) + (1-GP%CurrentD)*matmul(CDev,strain)
      endif

   end subroutine

   subroutine ComputeCTensorTangent(a,ElementSize,ElasticPredictor,GP)
      class(DamageMaterial), target :: a
      real(rp), intent(in) :: ElementSize
      real(rp), Intent(in) :: ElasticPredictor(:)
      class(DamageMaterial_GPD) :: GP

      real(rp), pointer :: CDev(:,:),CSph(:,:)

      !Normal damage
      if (.not. a%DeviatoricOnlyDamage) then
         !C Tensor

         GP%C = a%C*(1-GP%CurrentD)

      !Deviatoric only damage
      else
         call a%CT%GetPointerToCDeviatoric(CDev)
         call a%CT%GetPointerToCSpherical(CSph)

         GP%C = CSph + Cdev*(1-GP%CurrentD)
      endif

   end subroutine

   subroutine MD_GetStressTensorPointer(a,igaus,S)
      class(DamageMaterial_ED), target:: a
      integer(ip) :: igaus
      real(rp), pointer :: S(:)

      S => a%GP(igaus)%S(1:a%Material%CT%VoigtSize)
   end subroutine

   subroutine MD_GetConstitutiveTensorPointer(a,igaus,C)
      class(DamageMaterial_ED), target:: a
      integer(ip) :: igaus
      real(rp), pointer :: C(:,:)

      C => a%GP(igaus)%C
   end subroutine

   subroutine MD_SetElementSize(a,h)
      class(DamageMaterial_ED) :: a
      real(rp) :: h

      a%ElementSize = h
   end subroutine

   subroutine MD_MoveHistoryVariablesToConverged(a)
      class(DamageMaterial_ED) :: a

      integer(ip) :: igaus

      do igaus = 1,size(a%GP,1)
         if (a%GP(igaus)%CurrentTau > a%GP(igaus)%ReferenceTau) then
            a%GP(igaus)%ReferenceTau = a%GP(igaus)%CurrentTau
            a%GP(igaus)%ReferenceD = a%GP(igaus)%CurrentD
         endif
      enddo
   end subroutine

   subroutine MD_SetupPostprocess(a,Mesh,Memor)
      use Mod_Mesh
      use Mod_Memor
      use Mod_RpMeshAllocator
      implicit none
      class(DamageMaterial) :: a
      class(FemMesh) :: Mesh
      type(MemoryMan) :: Memor

      !Allocate array for gauss point info
      call AllocR1p(Memor,Mesh,a%DPostprocess)
      call AllocR1p(Memor,Mesh,a%TauPostprocess)
      call AllocR1p(Memor,Mesh,a%CompPostprocess)

   end subroutine

   subroutine MD_ContributeToPostprocess(a,ielem)
      class(DamageMaterial_ED) :: a
      integer(ip) :: ielem

      integer(ip) :: igaus

      do igaus = 1,size(a%GP,1)
         a%Material%DPostprocess(ielem)%a(igaus) = a%GP(igaus)%CurrentD
         a%Material%CompPostprocess(ielem)%a(igaus) = a%GP(igaus)%ComparisonRatio
         a%Material%TauPostprocess(ielem)%a(igaus) = a%GP(igaus)%CurrentTau
      enddo
   end subroutine

   subroutine MD_DoPostprocess(a,istep,ctime,FilePostpr,Mesh,Memor,iiter_char)
      use Mod_Mesh
      use Mod_Postpr
      use Mod_Memor
      use Mod_RpMeshAllocator
      use Mod_int2str
      implicit none
      class(DamageMaterial) :: a
      class(PostprFile) :: FilePostpr
      class(FemMesh) :: Mesh
      type(MemoryMan) :: Memor
      integer(ip) :: istep
      real(rp) :: ctime
      character(5) :: iiter_char

      call FilePostpr%postgp(a%DPostprocess,'Damage'//iiter_char,istep,ctime,Mesh)
      call FilePostpr%postgp(a%TauPostprocess,'Tau'//iiter_char,istep,ctime,Mesh)
      call FilePostpr%postgp(a%CompPostprocess,'ComparisonRatio'//iiter_char,istep,ctime,Mesh)

      !Deallocate array for GaussPoint info
      call DeallocR1p(Memor,Mesh,a%DPostprocess)
      call DeallocR1p(Memor,Mesh,a%TauPostprocess)
      call DeallocR1p(Memor,Mesh,a%CompPostprocess)
   end subroutine

   subroutine GetMaterialPointer(a,Material)
      class(DamageMaterial_ED), target :: a
      class(PLCDMaterial), pointer :: Material

      Material => a%Material

   end subroutine

   subroutine GetHistoryVariablesNdofn(a,ndofn)
      class(DamageMaterial_ED) :: a
      integer(ip) :: ndofn

      !In this subroutine the dimension of the history array to be communicated
      !in the adaptive process needs to be set
      !This subroutine is to be overwritten by each material

      ndofn = 1
   end subroutine

   subroutine SpecificHistoryGetToGaussArray(a,GaussArray)
      class(DamageMaterial_ED) :: a
      real(rp) :: GaussArray(:,:)

      integer(ip) :: igaus

      do igaus = 1,size(a%GP,1)
         GaussArray(1,igaus) = a%GP(igaus)%ReferenceTau
      enddo

   end subroutine

   subroutine SpecificHistoryAddFromGaussArray(a,GaussArray)
      class(DamageMaterial_ED) :: a
      real(rp) :: GaussArray(:,:)

      integer(ip) :: igaus

      !If this is being added from multiple elements, we need to set it to zero the first time
      if (a%JustCreated) then
         a%JustCreated = .false.
         a%GP(:)%ReferenceTau = 0.0_rp
      endif

      do igaus = 1,size(a%GP,1)
         !Reference Tau
         a%GP(igaus)%ReferenceTau = a%GP(igaus)%ReferenceTau + GaussArray(1,igaus)
      enddo
   end subroutine

   subroutine SpecificHistoryFinalRefinementUpdate(a)
      class(DamageMaterial_ED) :: a
      real(rp) :: ElasPred(a%Material%CT%VoigtSize)

      integer(ip) :: igaus

      do igaus = 1,size(a%GP,1)
         !Reference Tau
         a%GP(igaus)%ReferenceTau = max(a%GP(igaus)%ReferenceTau,a%Material%Tau0)
         !Current Tau
         a%GP(igaus)%CurrentTau = a%GP(igaus)%ReferenceTau

         !Damage parameter and C tensor
         if (a%GP(igaus)%ReferenceTau > a%Material%Tau0) then
            call a%Material%ComputeDamage(a%ElementSize,a%GP(igaus))
            call ComputeCT_Secant(a%Material,a%ElementSize,ElasPred,a%GP(igaus))
            a%GP(igaus)%ReferenceD = a%GP(igaus)%CurrentD
         endif
      enddo
   end subroutine

   subroutine CopyTo(a,EMD)
      use Mod_plcd_Material
      class(DamageMaterial_ED) :: a
      class(ElementMaterialData), pointer :: EMD

      integer(ip) :: igaus

      select type(EMD)
      type is (DamageMaterial_ED)
         do igaus = 1,size(a%GP,1)
            EMD%GP(igaus)%ReferenceTau = a%GP(igaus)%ReferenceTau
            EMD%GP(igaus)%ReferenceD   = a%GP(igaus)%ReferenceD
            EMD%GP(igaus)%CurrentTau   = a%GP(igaus)%CurrentTau
            EMD%GP(igaus)%CurrentD     = a%GP(igaus)%CurrentD

            if (EMD%GP(igaus)%CurrentD /= 0.0_rp) then
               EMD%GP(igaus)%IsDamaged = .true.
               allocate(EMD%GP(igaus)%C(a%Material%CT%VoigtSize,a%Material%CT%VoigtSize))
               EMD%GP(igaus)%C = a%GP(igaus)%C
            endif
         enddo
      end select

   end subroutine

   subroutine Finalize(a)
      class(DamageMaterial_ED) :: a

      integer(ip) :: igaus

      do igaus = 1,size(a%GP,1)
         deallocate(a%GP(igaus)%S)
         if (a%GP(igaus)%IsDamaged) then
            deallocate(a%GP(igaus)%C)
         endif
      enddo
      deallocate(a%GP)

   end subroutine

   subroutine MD_IsElementSizeRequiredByEMDs(a,isRequired)
      implicit none
      class(DamageMaterial), target :: a
      logical :: IsRequired

      isRequired = .true.
   end subroutine

   subroutine GetBufferSize(a,BufferSize)
      class(DamageMaterial_ED) :: a
      integer(ip) :: BufferSize

      BufferSize = size(a%GP,1)
   end subroutine

   subroutine ToBuffer(a,Buffer)
      class(DamageMaterial_ED) :: a
      real(rp) :: Buffer(*)

      integer(ip) :: igaus

      do igaus = 1,size(a%GP,1)
         Buffer(igaus) = a%GP(igaus)%ReferenceTau
      enddo
   end subroutine

   subroutine InitializeFromBuffer(a,Buffer)
      class(DamageMaterial_ED) :: a
      real(rp) :: Buffer(*)
      real(rp) :: ElasPred(a%Material%CT%VoigtSize)

      integer(ip) :: igaus

      do igaus = 1,size(a%GP,1)
         a%GP(igaus)%ReferenceTau = Buffer(igaus)
         a%GP(igaus)%CurrentTau = a%GP(igaus)%ReferenceTau

         !Damage parameter and C tensor
         if (a%GP(igaus)%ReferenceTau > a%Material%Tau0) then
            call a%Material%ComputeDamage(a%ElementSize,a%GP(igaus))
            call ComputeCT_Secant(a%Material,a%ElementSize,ElasPred,a%GP(igaus))
            a%GP(igaus)%ReferenceD = a%GP(igaus)%CurrentD
         endif
      enddo
   end subroutine

   subroutine GetComparisonCriteriaRatio(a,ComparisonRatio)
      class(DamageMaterial_ED) :: a
      real(rp) :: ComparisonRatio

      ComparisonRatio = maxval(a%GP(:)%ComparisonRatio)

   end subroutine

    !For u-p elements
   subroutine GetInverseVolumetricDeformationModulus(a,igaus,invK)
      class(DamageMaterial_ED) :: a
      integer(ip) :: igaus
      real(rp) :: invK

      call a%Material%CT%GetInverseVolumetricDeformationModulus(invK)

      if (.not. a%Material%DeviatoricOnlyDamage) invk = invk/(1.0_rp-a%GP(igaus)%CurrentD) !It must be divided, since it is the inverse!!!!
   end subroutine

   subroutine GetSecantCNorm(a,igaus,SecantCNorm)
      class(DamageMaterial_ED) :: a
      integer(ip) :: igaus
      real(rp) :: SecantCNorm
      real(rp), pointer :: CDev(:,:), CSph(:,:)

      real(rp) :: EigenValues(size(a%Material%C,1))

      !Normal Damage
      if (.not. a%Material%DeviatoricOnlyDamage) then
         call a%Material%CT%GetCNorm(SecantCNorm)
         SecantCNorm = SecantCNorm*(1.0_rp-a%GP(igaus)%CurrentD)

      !Deviatoric Only damage
      else
         call a%Material%CT%GetPointerToCDeviatoric(CDev)
         call a%Material%CT%GetPointerToCSpherical(CSph)

         EigenValues = 0.0_rp
         call MatrixEigenValues(size(CDev,1),CDev*(1-a%GP(igaus)%CurrentD)+CSph,Eigenvalues)
         SecantCNorm = maxval(abs(Eigenvalues))

      endif

   end subroutine

   subroutine SetPressureForComputeHistoryAndConstitutiveTensor(a,Pressure)
      class(DamageMaterial_ED) :: a
      real(rp) :: Pressure

      a%gpPressure = Pressure

   end subroutine

! Constitutive tensor options

   subroutine ComputeCT_Elastic(a,ElementSize,ElasticPredictor,GP)
      class(DamageMaterial) :: a
      real(rp), intent(in) :: ElementSize
      real(rp), intent(in) :: ElasticPredictor(:)
      class(DamageMaterial_GPD) :: GP

      GP%C = a%C

   end subroutine

   subroutine ComputeCT_Secant(a,ElementSize,ElasticPredictor,GP)
      class(DamageMaterial) :: a
      real(rp), intent(in) :: ElementSize
      real(rp), intent(in) :: ElasticPredictor(:)
      class(DamageMaterial_GPD) :: GP

      real(rp), pointer :: CDev(:,:), CSph(:,:)

      !Normal damage
      if (.not. a%DeviatoricOnlyDamage) then

         GP%C = a%C*(1-GP%CurrentD)

      !Deviatoric only damage
      else
         call a%CT%GetPointerToCDeviatoric(CDev)
         call a%CT%GetPointerToCSpherical(CSph)

         GP%C = CSph + Cdev*(1-GP%CurrentD)
      endif

   end subroutine

   subroutine ComputeCT_Tangent_VonMises_Exp(a,ElementSize,ElasticPredictor,GP)
      class(DamageMaterial) :: a
      real(rp), intent(in) :: ElementSize
      real(rp), intent(in) :: ElasticPredictor(:)
      class(DamageMaterial_GPD) :: GP

      real(rp), pointer :: CDev(:,:), CSph(:,:)
      real(rp) :: Aparam, gf, raux1, facTenTang, E0
      real(rp) :: dtau_de(a%CT%VoigtSize)
      integer(ip) :: is, js

      !Normal damage
      if (.not. a%DeviatoricOnlyDamage) then

         gf = a%FractureEnergy/ElementSize

         call a%CT%GetYoungModulusEquiv(E0)

         raux1 = gf*E0/(a%Tau0**2)-0.5_rp
         Aparam = 1/raux1

         dtau_de = ElasticPredictor/GP%CurrentTau
         facTenTang = ((a%Tau0 + Aparam*GP%CurrentTau)/GP%CurrentTau**2)*exp(Aparam*(1-GP%CurrentTau/a%Tau0))
         GP%C = a%C*(1-GP%CurrentD)

         Do is = 1, a%CT%VoigtSize
            Do js = 1, a%CT%VoigtSize
               GP%C(js,is) = GP%C(js,is) - facTenTang*ElasticPredictor(js)*dtau_de(is)
            End do
         End do

      !Deviatoric only damage
      else
         call a%CT%GetPointerToCDeviatoric(CDev)
         call a%CT%GetPointerToCSpherical(CSph)

         GP%C = CSph + Cdev*(1-GP%CurrentD)
      endif

   end subroutine

! Damage types

   subroutine ComputeDamage_Exp(a,ElementSize,GP)
      class(DamageMaterial) :: a
      real(rp), intent(in) :: ElementSize
      class(DamageMaterial_GPD) :: GP
      real(rp) :: E0, Aparam, gf, raux1

      !Fracture energy
      gf = a%FractureEnergy/ElementSize

      call a%CT%GetYoungModulusEquiv(E0)

      raux1 = E0*gf/(a%Tau0**2) - 0.5_rp
      Aparam = 1/raux1

      GP%CurrentD = 1.0_rp - a%Tau0/GP%CurrentTau*exp(Aparam*(1.0_rp-GP%CurrentTau/a%Tau0))
      GP%CurrentD = min(GP%CurrentD,a%DamageLimit)
      if (GP%CurrentD < 0.0_rp) call runend('Negative damage, fracture energy too low for current mesh')

   end subroutine

   subroutine ComputeDamage_Dev_Exp(a,ElementSize,GP)
      class(DamageMaterial) :: a
      real(rp), intent(in) :: ElementSize
      class(DamageMaterial_GPD) :: GP
      real(rp) :: Aparam, gf, raux1
      real(rp) :: Cnorm

      !Fracture energy
      gf = a%FractureEnergy/ElementSize

      !Cnorm is usually 2G
      call a%CT%GetCNorm(Cnorm)
      raux1 = gf*1.5_rp*Cnorm/(a%Tau0**2)-0.5_rp
      Aparam = 1/raux1

      GP%CurrentD = 1.0_rp - a%Tau0/GP%CurrentTau*exp(Aparam*(1.0_rp-GP%CurrentTau/a%Tau0))
      GP%CurrentD = min(GP%CurrentD,a%DamageLimit)
      if (GP%CurrentD < 0.0_rp) call runend('Negative damage, fracture energy too low for current mesh')

   end subroutine

end module

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


program Pointer_array
    use mod_pt_array
    use mod_pt_array_utils
    implicit none

    integer                :: i, posi
    real                   :: x
    type(element), pointer :: pt_elem
    type(obj_elem)         :: elem_obj

    nbelements = 20

    call random_seed()

    allocate(elements(1:nbelements))
    do i = 1, nbelements
        allocate(elements(i)%p)
        call random_number(x)
        elements(i)%p  = element_constructor(i,x)
    end do

    write(*,*) 'Ids and values of elements:'
    do i = 1, nbelements
        elem_obj%p => elements(i)%p
        call elem_obj%print_elem
    end do

    do while (posi /= 0)
        write(*,*) 'Please enter the position to delete an element in pointer array:'
        write(*,*) '(Program will be terminated when typing "0").'
        read(*,*) posi
        if (posi == 0) then
            stop 'Good bye.'
        end if

        call delete_elements_at_position(posi)

        write(*,*) 'Ids and values of elements after deleting:'
        do i = 1, nbelements
            pt_elem => elements(i)%p
            elem_obj%p => pt_elem
            call elem_obj%print_elem
        end do
    end do

end program Pointer_array


text:
c###serprb.for
      SUBROUTINE SERPRB(IF)
C--------------------------------
C
C     THIS ROUTINE FINDS THE SERVICE PROBABILITY
C
C     THIS ROUTINE HAS NOT BEEN REVISED TO USE SUMMED SIGNAL
C     TLLOW(IM) = DBLOSL(IF)
C     SN(IM) = SNDB(IF)
C     TLHGH(IM) = DBLOSU(IF)
C     IM = IF
C     DO ONLY ONCE FOR IF
C
      COMMON/ANOIS/ATNU,ATNY,CC,TM,RCNSE,DU,DL,SIGM,SYGU,SYGL,KJ,JK
      COMMON / CON / D2R, DCL, GAMA, PI, PI2, PIO2, R2D, RZ, VOFL
      COMMON /DON /ALATD, AMIN, AMIND, BTR, BTRD, DLONG, DMP, ERTR, GCD,
     1 GCDKM, PMP, PWR, TLAT, TLATD, TLONG, TLONGD, RSN, SIGTR, RLAT,
     2 RLATD,RLONG,RLONGD,BRTD,FLUX,ULAT,ULATD,ULONG,ULONGD,SSN,D90R,
     3 D50R,D10R,D90S,D50S,D10S
      COMMON / ION / IHRE, IHRO, IHRS, LUFP, METHOD, NOISE, NPSL
      COMMON / SON / ANGLE(13), ANGLER(13), CPROB(13), DBLOS(13),
     A DBLOSL(13), DBLOSU(13), DBU(13), DELAY(13), DBW(13), NHP(13),
     B XNYNOIS(13), PROBMP(13), RELIAB(13), SNDB(13), SNPR(13),
     C SNRLW(13), SNRUP(13), SPROB(13), VHIGH(13), RNEFF(13),MDL(13)
       CHARACTER MDL*1
      COMMON/TON/ADJ,ADS,GNOS,GOT,REL,SL,SLS,SPR,SU,SUS
     A ,XNOISE,ZNOISE,NF
      COMMON / allMODE /ABPS(20),CREL(20),FLDST(20),HN(20),HP(20),
     1PROB(20),RELY(20),RGAIN(20),SIGPOW(20),SN(20),
     2SPRO(20),TGAIN(20),TIMED(20),TLOSS(20),B(20),FSLOS(20),
     3grlos(20),adv(20),obf(20),
     CNMODE(20),TLLOW(20),TLHGH(20),EFF(20),NREL,NMMOD
      DIMENSION TME(10),d10sa(20),d50sa(20)
C.....NORMAL DISTRIBUTION
      DATA TME/0.0,0.1257,0.2533,0.3853,0.5244,0.6745,0.8416,1.0364
     A , 1.2815, 1.6449 /
C.....DR IS THE PREDICTION ERROR IN RSN, REQUIRED SNR
      DATA DR/2./
      SPR = 0.0
C.....LUFP IS THE REQUIRED RELIABILITY
      ITM = IABS( LUFP - 50) /5 +1
      ITM = MIN0(ITM,10)
      TMX = TME(ITM)
C.....BEGINNING OF FINDING THE SERVICE PROBABILITY FOR EACH MODE
      DO 145 IM = 1,NMMOD
      IF(HP(IM).le.70.) then     !  set default for no mode
         SPRO(IM) = .001
         d10sa(im)=0.
         d50sa(im)=0.
         GO TO 145
      end if
      IF(LUFP.ge.50) then
         DN  = TMX *DU/DCL                !  DN is the noise
         DS  = TMX *TLLOW(IM)/DCL         !  DS is the signal
         XLH = -1.
C.....SIGNAL PREDICTION ERROR
         DSO = TMX *SUS
         DNO = TMX * SYGU/DCL
      else
         DN  = TMX *DL/DCL                !  DN is the noise
         DS  = TMX *TLHGH(IM)/DCL         !  DS is the signal
         XLH = 1.
         DSO = TMX * SLS
         DNO = TMX *SYGL/DCL
      end if
C.....D10S IS THE LOWER DECILE
C.....D50S IS THE MEDIAN
C.....D90S IS THE UPPER DECILE
      D50Sa(im)= SQRT(DN*DN + DS*DS)
      D10Sa(im)= D50Sa(im) +
     +        SQRT(SIGM*SIGM + ADS*ADS + DNO*DNO + DSO*DSO + DR*DR)
      D50Sa(im) =  SN(IM)  + XLH*D50Sa(im)
      Z = (RSN - D50Sa(im))/d10sa(im)
C.....SERVICE PROBABILITY FOR THIS MODE
      SPRO(IM) =  1. - FNORML(Z)
  145 CONTINUE
C.....END OF FINDING THE SERVICE PROBABILITY FOR EACH MODE
C
C.....USE THE MAXIMUM
      AMXX=-1000.
      imax=1
      DO 150 I=1,NMMOD
      IF(SPRO(I).GT.SPRO(imax)) imax=i
  150 CONTINUE
      SPROB(IF)=SPRO(imax)
      D10S=d10sa(imax)
      D50S=d50sa(imax)
      D90S=D10S
      RETURN
      END
C--------------------------------

ext: for

max_stars_repo_name: jawatson/voacapl

max_issues_repo_name: jawatson/voacapl

max_forks_repo_name: jawatson/voacapl


text:
!! R527 allocatable-decl
!    is object-name [ ( array-spec ) ] [ lbracket co-array-spec rbracket ]

   ALLOCATABLE coscalar[*]
   allocatable :: coarray[5,3:7,*]
   allocatable :: combo(:,:)[*]

end

max_stars_repo_name: OpenFortranProject/ofp-sdf

max_issues_repo_name: OpenFortranProject/ofp-sdf

max_forks_repo_name: OpenFortranProject/ofp-sdf


text:
      subroutine daxpy(nmax,scfact,array1,incx,array2,incy)
c
c     This subroutine multiplies the real*8 array "array1" by the
c     real*8 constant scfact and adds that to the real*8 array "array2".
c     The result is returned in "array2". The relevant equation is
c     (y) = a*(x) + (y). This subroutine is an adaptation of the 1979
c     Linpack subroutine of the same name. The size and order in the
c     original calling sequence has been preserved:
c     (n,da,dx,incx,dy,incy) == (nmax,scfact,array1,incx,array2,incy).
c     The increments incx and incy are not used here (each is taken
c     as having a value of 1; input of any other value constitutes
c     an error). This is a pseudo-Linpack BLAS (Basic Linear Algebra
c     Subsystem) subroutine.
c
c     This subroutine is called by:
c
c       EQLIBU/dgefa.f
c       EQLIBU/dgesl.f
c
c-----------------------------------------------------------------------
c
c     Input:
c
c       nmax   = dimension or pseudo-dimension of array
c       scfact = the scale factor
c       array1 = the first array (x)
c       array2 = the second array (y)
c       incx   = the x increment (must have a value of 1)
c       incy   = the y increment (must have a value of 1)
c
c     Output:
c
c       array2 = the original array multiplied by the scale factor.
c
c
c-----------------------------------------------------------------------
c
      implicit none
c
c-----------------------------------------------------------------------
c
c     Calling sequence variable declarations.
c
      integer nmax
      integer incx,incy
c
      real*8 array1(nmax),array2(nmax)
c
      real*8 scfact
c
c-----------------------------------------------------------------------
c
c     Local variable declarations.
c
      integer i,ileft
c
      real*8 scx
c
c-----------------------------------------------------------------------
c
c     Trap illegal values of incx and incy.
c
      if (incx .ne. 1) then
        write (6,1000) incx
 1000   format(/' * Error - (EQLIBU/daxpy) The argument incx has a',
     $  /7x,'value of ',i5,'. This version of daxpy only allows this',
     $  /7x,'argument to have a value of 1. This is a programming',
     $  /7x,'error or a linking error.')
        stop
      endif
c
      if (incy .ne. 1) then
        write (6,1010) incy
 1010   format(/' * Error - (EQLIBU/daxpy) The argument incy has a',
     $  /7x,'value of ',i5,'. This version of daxpy only allows this',
     $  /7x,'argument to have a value of 1. This is a programming',
     $  /7x,'error or a linking error.')
        stop
      endif
c
c     Note that the loop is unrolled.
c
      scx = scfact
      ileft = (nmax/8)*8
c
      do i = 1,ileft,8
        array2(i) = array2(i) + scx*array1(i)
        array2(i + 1) = array2(i + 1) + scx*array1(i + 1)
        array2(i + 2) = array2(i + 2) + scx*array1(i + 2)
        array2(i + 3) = array2(i + 3) + scx*array1(i + 3)
        array2(i + 4) = array2(i + 4) + scx*array1(i + 4)
        array2(i + 5) = array2(i + 5) + scx*array1(i + 5)
        array2(i + 6) = array2(i + 6) + scx*array1(i + 6)
        array2(i + 7) = array2(i + 7) + scx*array1(i + 7)
      enddo
c
      do i = ileft + 1,nmax
        array2(i) = array2(i) + scx*array1(i)
      enddo
c
      end


text:
! Auto-generated. DO NOT MODIFY!
module raw_tensor_permute_module
    use :: tensor_api, only : &
            tensor, &
            matrix, &
            tensor_rank3, &
            tensor_rank4, &
            tensor_rank5

    use :: tensor_permute_module, only : tensor_permute

    implicit none
    private

    public :: raw_tensor_permute

    type, abstract, extends(tensor_permute) :: raw_tensor_permute
    contains
        procedure :: is_compatible => is_compatible
    end type raw_tensor_permute
contains
    logical function is_compatible(this, src, dst, perm)
        class(raw_tensor_permute), intent(in) :: this
        class(tensor), intent(in) :: src, dst
        integer, dimension(:), intent(in) :: perm

        is_compatible = same_type_as(src, dst)

        select type (src)
        type is ( matrix )
            is_compatible = is_compatible .and. size(perm) == 2
        type is ( tensor_rank3 )
            is_compatible = is_compatible .and. size(perm) == 3
        type is ( tensor_rank4 )
            is_compatible = is_compatible .and. size(perm) == 4
        type is ( tensor_rank5 )
            is_compatible = is_compatible .and. size(perm) == 5
        class default
            is_compatible = .false.
        end select

        if ( is_compatible ) is_compatible = all(dst%dims == src%dims(perm))
    end function is_compatible
end module raw_tensor_permute_module

max_stars_repo_name: cheshyre/ntcl-algorithms

max_issues_repo_name: cheshyre/ntcl-algorithms

max_forks_repo_name: cheshyre/ntcl-algorithms


!##############################################################################
! PROGRAM matrices
!
! ## Operations on matrices
!
! This code is published under the GNU General Public License v3
!                         (https://www.gnu.org/licenses/gpl-3.0.en.html)
!
! Authors: Hans Fehr and Fabian Kindermann
!          contact@ce-fortran.com
!
! #VC# VERSION: 1.0  (23 January 2018)
!
!##############################################################################
program matrices

    implicit none
    integer :: i, j
    real*8 :: a(4), b(4)
    real*8 :: x(2, 4), y(4, 2), z(2, 2)

    ! initialize vectors and matrices
    a = (/(dble(5-i), i=1, 4)/)
    b = a+4d0
    x(1, :) = (/1d0, 2d0, 3d0, 4d0/)
    x(2, :) = (/5d0, 6d0, 7d0, 8d0/)
    y = transpose(x)
    z = matmul(x, y)

    ! show results of different functions
    write(*,'(a,4f7.1/)')'         vector a = ',(a(i), i=1,4)
    write(*,'(a,f7.1)')  '           sum(a) = ',sum(a)
    write(*,'(a,f7.1/)') '       product(a) = ',product(a)
    write(*,'(a,f7.1)')  '        maxval(a) = ',maxval(a)
    write(*,'(a,i7)')    '        maxloc(a) = ',maxloc(a)
    write(*,'(a,f7.1)')  '        minval(a) = ',minval(a)
    write(*,'(a,i7/)')   '        minloc(a) = ',minloc(a)
    write(*,'(a,4f7.1)') '    cshift(a, -1) = ',cshift(a, -1)
    write(*,'(a,4f7.1/)')'   eoshift(a, -1) = ',eoshift(a, -1)
    write(*,'(a,l7)')    '       all(a<3d0) = ',all(a < 3d0)
    write(*,'(a,l7)')    '       any(a<3d0) = ',any(a < 3d0)
    write(*,'(a,i7/)')   '     count(a<3d0) = ',count(a < 3d0)
    write(*,'(a,4f7.1/)')'         vector b = ',(b(i), i=1,4)
    write(*,'(a,f7.1/)') ' dot_product(a,b) = ',dot_product(a, b)
    write(*,'(a,4f7.1/,20x,4f7.1/)') &
                         '         matrix x = ',((x(i,j), j=1,4), i=1,2)
    write(*,'(a,2f7.1,3(/20x,2f7.1)/)')&
                         '     transpose(x) = ',((y(i,j), j=1,2), i=1,4)
    write(*,'(a,2f7.1/,20x,2f7.1/)') &
                         '      matmul(x,y) = ',((z(i,j), j=1,2), i=1,2)

end program


text:
7       	 <--SHAPES
7       	 <--LINES
id1
2       	 <--TYPE
275       	 <--LEFT
74       	 <--TOP
70       	 <--WIDTH
30       	 <--HEIGHT
16777215       	 <--BACKCOLOR
0       	 <--BORDERCOLOR
0       	 <--BORDERCOLOR
-reserved 1-
-reserved 2-
START

id2
2       	 <--TYPE
303       	 <--LEFT
650       	 <--TOP
70       	 <--WIDTH
30       	 <--HEIGHT
16777215       	 <--BACKCOLOR
0       	 <--BORDERCOLOR
0       	 <--BORDERCOLOR
-reserved 1-
-reserved 2-
STOP

id3
91       	 <--TYPE
245       	 <--LEFT
137       	 <--TOP
332       	 <--WIDTH
40       	 <--HEIGHT
16777215       	 <--BACKCOLOR
0       	 <--BORDERCOLOR
0       	 <--BORDERCOLOR
-reserved 1-
-reserved 2-
INPUT
x degeri girin x y'den buyuk olsun
x

id4
91       	 <--TYPE
245       	 <--LEFT
203       	 <--TOP
172       	 <--WIDTH
40       	 <--HEIGHT
16777215       	 <--BACKCOLOR
0       	 <--BORDERCOLOR
0       	 <--BORDERCOLOR
-reserved 1-
-reserved 2-
INPUT
y degeri girin
y

id5
0       	 <--TYPE
562       	 <--LEFT
413       	 <--TOP
108       	 <--WIDTH
30       	 <--HEIGHT
16777215       	 <--BACKCOLOR
0       	 <--BORDERCOLOR
0       	 <--BORDERCOLOR
-reserved 1-
-reserved 2-
DIVIDE
x
y
son
id6
91       	 <--TYPE
266       	 <--LEFT
408       	 <--TOP
124       	 <--WIDTH
40       	 <--HEIGHT
16777215       	 <--BACKCOLOR
0       	 <--BORDERCOLOR
0       	 <--BORDERCOLOR
-reserved 1-
-reserved 2-
OUTPUT
sonuc:
son

id7
92       	 <--TYPE
559       	 <--LEFT
235       	 <--TOP
112       	 <--WIDTH
50       	 <--HEIGHT
16777215       	 <--BACKCOLOR
0       	 <--BORDERCOLOR
0       	 <--BORDERCOLOR
-reserved 1-
-reserved 2-
IF_GREATER
x
y

---- LINES ---- from,to ----
id1,id3
reserved 1

id3,id4
reserved 1

id5,id6
reserved 1

id6,id2
reserved 1

id4,id7
reserved 1

id7,id5
reserved 1
EVET
id7,id3
reserved 1
HAYIR


text:
Outdoor Adventures has moved from their old location next to Bike Barn by the Silo. They are now located at the former Rec Pool Lodge.

Outdoor Adventures is one of the largest university outdoor programs in the United States. They run trips and classes in a variety of disciplines year round. Part of Campus Recreation, students and ARC members get a special discount, although trips are open to the general public. Trips are somewhat Cooperatives cooperative in nature, meaning that participants will often carpool to locations and coordinate food on multiday trips. They offer a wide variety trips every weekend as well as certain types of night classes and longer trips during vacation periods. While most trips require participants to be 18 years or older Outdoor Adventures offers fun for all ages, including http://campusrecreation.ucdavis.edu/cms/internal.aspx?uidcaee1f14811a42e5b3addbd8ee9f386b youth programs and trips that are held during the summer.

They also have a great Rental Center where you can rent most types of equipment that youd need to do one of the trips on your own, for a reasonable price. They dont rent downhill skis or snowboards, but you can find telemark skis, bear canisters, sleeping bags, tents, snowshoes, backpacks, and a host of other adventurerelated equipment. A price list can be found on their http://cru.ucdavis.edu/content.cfm?contentID148 website, or you can pick up a form in person. It is also a great place to hangout  good music, fun people, a variety of books, picnic tables outside and comfy couches inside.

Twice a year on the third Thursday of October and April, Outdoor Adventures sponsors an annual Equipment Swap at which you can get some fantastic deals on new and used equipment  just be sure to have cash!

Trip and class information and registration is available online at http://cru.ucdavis.edu, by phone (530) 7521995, or in person.

 Water Program
 Rafting
OA is one of the largest commercial rafting outfitters in the state. They hold several commercial permits that allow for nearly 100 people to go rafting a day and have a fleet of over 15 whitewater rafts. Outdoor Adventures utilizes expeditionsized oar boats to haul gear on longer trips. OAs rafting program is mainly focused on class III rivers, though they offer a few class IV and V trips during the Spring and Summer. Rafting Guide School, aka BGT, takes place over Spring break.
  South Fork of the American River
 Only an hour from Davis, this Class III damfed river is one of the most popular rivers in the world. Its a wonderful introduction to the sport and a great way for anyone to get out and have fun without having to drive far. The South Fork is a great river for beginner, intermediate, and advanced rafters alike. There are two different runs on the South Fork: the upper Chili Bar run and the lower Gorge run. OA rafts both: the Gorge is run in Spring & Summer Chili Bar in the Fall). When the water is cold like (April, May & October) OA provides wetsuits for rafters. Check in with OA for all the details of things you should and should not bring on a rafting trip. South Fork trips are available every Saturday and Sunday, April through October. Private trips may also be arranged, either on the weekend or during the week.
  Klamath River
 Klamath trips are three days long and thus allow you to experience a good mix of camping and rafting. The http://en.wikipedia.org/wiki/Klamath_River Klamath features rapids up to class III and is good for rafters of all levels. This scenic river is located in the beautiful and rugged woods of northern California. Bald eagle and other wildlife spottings are not uncommon on Klamath trips. Excellent food is provided as well. Klamath trips take place over the threeday weekends of Memorial Day, Fourth of July, and Labor Day.
  North Fork of the American River
 The bigger brother of the South Fork, this section is a step up in difficulty and is encouraged for intermediate to advanced rafters. Exhilarating Class IV rapids greet rafters during much of the journey and OA encourages people to have decent paddling experience before they sign up. Since trips are usually in spring and early summer, the water is quite cold and OA provides all rafters with wetsuits.  The North Fork is strictly dependent on water flows from snow melt since the river flows arent regulated by a dam. Hence, it is usually only runs in spring to early summer. After that, the flow tends to become too low for rafting. Exact trip dates vary from April to June.
  CalSalmon
 This twoday trip is open to advanced rafters only as it is a Class V river. It is a 2day trip; oar boats may be the only type of raft used for this river. The http://en.wikipedia.org/wiki/Salmon_River_(California) CalSalmon is given that moniker as there is another Salmon River in Idaho that is also popular with rafters, though in a far more remote location. The CalSalmon is run twice each May, waterflow and weather dependent.

 Whitewater Kayaking
 Brought back to the program by Jordy Margid in 2005, whitewater kayaking is an absolute blast. Classes are taught on the South Fork of the American River in the spring. This sport takes a fair amount of practice so its highly reccomended to go to a rolling clinic offered by Outdoor Adventures beforehand.  OA has a max student to teacher ratio of 3:1. Equipment that is provided on all trips: kayak, paddle, spray skirt, PFD, helmet, wetsuit and splash jacket.
  Open Rolling
 OA provides access to the Rec Pool so you can practice your http://en.wikipedia.org/wiki/Kayak_roll kayak rolling. No instruction is provided. Open rolling is every Tuesday night, April through September.
  Kayak Rolling Clinics
 OAs Rolling clinic is a twosession course and open to whitewater and sea kayakers alike. OA provides individual attention and the equipment. Clinics are offered frequently on concurrent Tuesdays, April through September.
  White Water Kayaking, Level I
 Open to all, this oneday course takes place on the Class II section of the South Fork of the American River. Classes are held April through September.
  White Water Kayaking, Level II
 This intermediate course is also one day, but meant for kayakers who have taken the Level I course and have taken the rolling clinic. Being able to roll is highly recommended. Classes are held April through September.

 Sailing
Sailing is no longer offered at OA. OA used to partner with a campus department, to use their 40 yacht and offer educational trips to students. Trips took place in and around the San Francisco Bay.

 Sea Kayaking
 OA has a fleet of 11 single and 5 double http://www.necky.com/ Necky sea kayaks. Transportation is provided on all trips, with special situation for Baja. Sea Kayaking Guide School is in May.
  Tomales Bay Day Tour
 The focus of this trip is just to get out and have fun! Some instruction is provided (depending on what people want) but really youre just out there to explore. Youll likely be able to see jellyfish, sea stars, seals, birds, elk and lots of other wildlife. Trips are held a couple of times per quarter, year round.
  Introduction to Sea Kayaking
 This class is designed for those who wish to learn the basic skills required to ocean kayak safely. The pretrip meeting covers a lot of material such as basics of navigation, planning and equipment. The day spent on http://en.wikipedia.org/wiki/Tomales_Bay Tomales Bay covers a large variety of skills including strokes, bracing and rescue. Trips are held about every other weekend, year round.
  Expedition Sea Kayaking
 Similar kayaking skills are taught on this two day trip as on the one day Intro class. The major difference is that you will be camping out on a wonderfully secluded beach, so you will learn how to camp out of your kayak. This also gives you much more time to practice skills, explore the area and just relax. Trips are once or twice per quarter, year round. There is also a threeday option that takes place during winter  this simply you gives you more time to explore the coastline!
  San Francisco Bay Intermediate Tour
 This tour is great for people who have some kayaking experience and want to do something a little more challenging. Most of the time, the trip launches from Sausilito and heads out towards Belvedere and http://en.wikipedia.org/wiki/Angel_Island,_California Angel Island. http://en.wikipedia.org/wiki/San_Francisco_Bay San Francisco Bay is a stunningly beautiful place to paddle. It can be a little challenging however. Strong currents, strong winds, boat traffic, fog and choppy water are all serious concerns. If conditions are right, you may have a chance to do a complete circumnavigation of Angel Island. Youll also get to paddle past some of the most valuable residential property in the world. This trip is offered once or twice per quarter, year round.
  Lake Tahoe
 Spend a bluebird day exploring the dramatic shores of http://en.wikipedia.org/wiki/Emerald_Bay_(Lake_Tahoe) Emerald Bay. This trip is open to all levels of kayakers and takes place during the summer.
 Ahjumawi Lava Springs
 This threeday adventure takes you up north of http://en.wikipedia.org/wiki/Mt._Lassen Mt. Lassen to http://en.wikipedia.org/wiki/Ahjumawi_Lava_Springs_State_Park Ahjumawi State Park.  A great place to paddle and explore the springs, cinder cones, and view some wildlife. Enjoyable for all levels of kayakers. This trip is offered over the threeday weekends of Memorial Day, Fourth of July, and Labor Day.
  Baja Sea Kayaking
 Baja Sea Kayaking trips through OA are no longer offered due to travel warnings. This is a sixday trip launching from http://en.wikipedia.org/wiki/Bahía_de_los_Ángeles Bahía de los Ángeles.  The drive is a little less than 2 days each way.  The itinerary varies from trip to trip, but there are a dozen or so campsites that you could visit including some beautiful white sand beaches, a dormant volcano, uninhabited islands and protected coves filled with clams.  The trip is not exclusively kayaking, as there are many opportunities to go Hiking and Backpacking hiking, snorkeling or just bumming around on the Beaches beach.  You are likely to see sea lions, migrating birds, fish, sting rays and possibly some dolphins or whales.  There are typically two trips over the winter break and two around spring break.
 Image(bahia.jpg, thumbnail, 600, noborder, View of Bahia de los Angeles from half way up a volcano)

 Wilderness Program
Offering a wide range of activities, the wilderness program offers trips all over the state of California. Activities vary seasonally due to snow. Wilderness Guide School is in April.
 Backpacking
Backpacking is an excellent yearround activity, in light of our prime location. Within a few hours, we have Pt. Reyes, Big Sur, Cache Creek, Mendocino, Mt. Lassen, Lake Tahoe/Sierras, and Yosemite.
  Backpacking Adventures
 Explore the back country and let the impressive scenery revitalize you. The course covers topics relating to safety and http://lnt.org/index.php LeaveNoTrace camping techniques. Trips are designed for the basiclevel backpackers, but most levels can be accommodated.  Group camping gear is included and personal gear is available from OA. Trips are offered nearly every weekend, year round.
  Introduction to Backpacking
 Geared toward getting you on your way toward adventures of your own, this trip is more instructional in structure (i.e. how to properly pack your backpack, ten set up, etc.).  A great way to improve your skills and get away for the weekend. Trips are offered several times per quarter during the Spring, Summer and Fall Quarters.
  Wilderness Navigation
 Learn the basics of finding your way through the back country with topographical maps, compasses, and the most important tool, common sense. Topics include map reading, on and offtrail route finding, advanced compass techniques, Leave No Trace Principles, and what to do when you get lost. This course is offered once per quarter, year round.
  Custom Backpacking Trips
 If you are looking for a trip that we do not offer, then find your friends and book a custom trip through OA.  You can choose the time and place that interests you most and OA will provide the group gear and our experiences guides to help lead you along your adventures.  These trips are a great way for any level of backpacker to gain the knowledge and skills for back country travel.
  Yosemite Photography Workshop
 Spend the weekend honing your photography skills at one of the most majestic parks in the country. OA guides have created an intense itinerary based on years of trips, taking into account seasonal highlights  Fall foliage, Winter snow, Spring flows, and Summer flowers. Relax at the end of the day by camping in the Valley, cooking a heaping group dinner, and roasting up http://en.wikipedia.org/wiki/Smores smores. This trip occurs only once per quarter, year round.
  Fly Fishing
 There are two courses on fly fishing, offered sporadically depending on demand. There is a 4session class offered on Flytying. There is also a weekend carcamping trip to the Sierras where you get handson experience. Check for OA if youre interested.

 Rock Climbing
Rock climbing is a great mental and physical challenge. It requires thoughtful approach, balance, flexibility and the ability to think several steps ahead. Davis is close to some very good climbing within 12 hours drive. Rock Climbing Guide School is in October.
  Beginning Rock Climbing
 This class is an introduction to the challenges of rock climbing.  The twopart class includes an evening session at the ARC Climbing Wall and a full day outdoors at http://en.wikipedia.org/wiki/Cosumnes_River Cosumnes River Gorge. We cover all the basic skills needed for technical climbing, including belaying, basic knots, and techniques for face and crack climbing on top rope. Shoes and all climbing gear are included. Trips are offered several times per quarter, year round.
  Introduction to Top Rope
 Take the next step and learn how to set up your own http://en.wikipedia.org/wiki/Top_rope Top Rope. Guides will go over techniques for setting up a safe Top Rope then go have some fun and climb! Trips are offered a few times per quarter, during the Spring and Fall Quarters.
  MultiPitch Climbing
MuliPitch Climbing classes are no longer offered through OA. Learn the basics of http://en.wikipedia.org/wiki/Lead_climbing lead climbing and http://en.wikipedia.org/wiki/Multipitch_climbing multipitch technique.  It is strongly recommended that you have a solid foundation in climbing and/or have taken Beginning Rock Climbing or Top Rope before enrolling. Classes are held at Lovers Leap and offered once per quarter during the Spring and Fall.
  Owens Valley Climbing
 Spend you Memorial Day Weekend doing something, er, memorial. Spend three days climbing in the famed http://en.wikipedia.org/wiki/Bishop,_California Eastern Sierras. Climbing options are limitless and can be tailored to your abilities. This trip is only offered over Memorial Day Weekend in May.

 Winter Sports
Davis is known for being fogged in during the winter. The solution for those of us wanting sun? Head to the coast (backpacking or sea kayaking) or climb out of the hills and head to the mountains. Reach for the sun. The Sierras get around 30 feet of snow every year yet there are plenty of bluebird days. Transportation is provided on all winter trips. Winter Guide School is in January.
  Wilderness Snowshoe Tour
 These great oneday trips focus primarily on going out and having fun. Light instruction in technique and winter travel is discussed. Trips often go to http://en.wikipedia.org/wiki/Castle_Peak_(California) Castle Peak or http://en.wikipedia.org/wiki/Echo_Lake_(California) Echo Lake. The OA armada of snowshoes are high quality shoes made by Atlas and MSR.
  Snow Cave and Igloo Building
 Fun yet practical, this course offers skills for any back country winter traveler. These naturalmaterial shelters can be made for emergency situations when snowboarding, crosscountry skiing, and snow shoeing. They can also be used as a lightweight solution to wintermountaineering tents. And, how many people can say theyve built and slept in a snow cave or igloo? This trip is offered everyother weekend in the Winter Quarter.
  Yosemite Snowshoe Weekend
 Leaving from http://en.wikipedia.org/wiki/Badger_Pass_Ski_Area#Current_day_Badger_Pass Glacier Point Road, these overnight trips are a great introduction to being comfortable while camping in the snow. You travel on snowshoes with heavy packs for a couple of miles to make camp in what is often a stunningly beautiful area  Dewey Point. Its definitely a trip you wont forget! This trip is offered everyother weekend during the Winter Quarter
  Yosemite Winter Exploration
 For those of you seeking a less intense experience then camping in the snow at Yosemite, you can settle for the creature comforts of car camping in http://en.wikipedia.org/wiki/Yosemite_Valley Yosemite Valley. These threeday trips take place over Presidents Weekend and MLK Weekend. Youll still be snowshoeing during the day.
  Intro to Cross Country Skiing
 Travel effortlessly through the flat meadows and rolling hills of the Sierras. Instruction in http://en.wikipedia.org/wiki/Cross_country_skiing cross country skiing is a great way to get in to the sport if you have limited skiing experience. OA uses very friendly (metal edged, waxless, slightly parabolic) skis. Trips often take place around Castle Peak but other locales are also used. Basic instruction in ski technique including kicking and gliding, going uphill, stopping and other fun stuff is practiced. This trip is offered everyother weekend in the Winter Quarter.
 Mt. Lassen CrossCountry Ski
 Probably the most intense beginner/intermediate winter trip that OA offers, http://en.wikipedia.org/wiki/Lassen_Volcanic_National_Park Lassen National Park is a great place to test yourself. You camp in the snow each night, so gear and food choice is quite important. A typical trip has the group set up base camp about half way up the mountain. Day 2 is spent skiing up to the snowcovered http://en.wikipedia.org/wiki/Lake_Helen_(Lassen_Peak) Lake Helen, at the base of http://en.wikipedia.org/wiki/Lassen_Peak Mt. Lassen. The peak is too steep to summit with crosscountry skis. Day 3 is an easy downhill descent to the trail head. This trip takes place over Presidents Day Weekend in February.
  Ice Ax and Crampon Clinic
 Safely learn the basics of using an http://en.wikipedia.org/wiki/Ice_ax ice ax and http://en.wikipedia.org/wiki/Crampons crampons in this oneday class. Before you try to climb a big mountain, like http://en.wikipedia.org/wiki/Mount_Shasta Shasta, get some formal education through a course like this. Then go out and practice, practice, practice. Ice ax and crampons are provided. This course is only offered once in the Winter Quarter, so make sure to register early!
  Intro to Avalanche Safety
 This oneday class is offered for the benefit and safety of the community. Take this basic intro class before you head out snowshoeing, snowmobiling or skiing (cross country, telemark, or Alpine Touring). A lot of information is covered in the classroom and weekend day session. If you plan on doing extensive winter backcountry travel, it is still recommended that you continue your avalanche education and get http://avtraining.org/AvalancheTrainingCourses/AIARELevel1.html AIARE Level I certified, for your own good. This course is only offered once in the Winter Quarter.
  Intro to Ski Mountaineering
 Ski untouched pow. Escape the crowds. Save $79 a day by not having to buy a lift ticket. Earn your turns. See nature, not people. Push your limits. Whatever your reason, getting out of the resorts and into the backcountry has never been easier. Learn the tips and techniques on how to explore the backcountry on http://en.wikipedia.org/wiki/Alpine_touring Alpine Touring, http://en.wikipedia.org/wiki/Telemark_skiing Telemark or http://en.wikipedia.org/wiki/Splitboard Splitboards. OA provides the ski gear, you provide the beacon, legs muscles, and desire to learn. You must be a comfortable intermediate skier to take this course, as well as in good physical shape. OA doesnt teach you how to ski, they teach you how to take your skiing skills to the backcountry. The Into to Avalanche Safety course is recommended.

 Healthcare Courses
Healthcare Instructor School is in October.
 CPR
Outdoor Adventures offers two levels of CPR classes, taught to http://www.americanheart.org/presenter.jhtml?identifier1200000 American Heart Association curriculum.  Students should decide which class to take based on their needs (layperson vs. professional).
  Healthcare Provider (HCP)
 A comprehensive course which teaches CPR, http://en.wikipedia.org/wiki/Aed AED, rescue breathing, and foreign body airway obstruction for adult, child and infant patients. Class is 4.5 hours long. Renewal course also available. Certification is for two years. This course is offered about once a week year round.
  HeartSaver AED
 Focuses on CPR and AED for adult, child and infants. Slightly easier and quicker than the HCP course; it only takes 4 hours.  Certification is for two years. This course is offered about once a week year round.
 First Aid
  HeartSaver First Aid
 Does not include CPR. Includes basic First Aid used by lay rescuers in a professional setting, including teachers, camp counselors, and others who do not need a higher level of first aid training.  Certification is two years. This course is offered several times each quarter.
  Wilderness First Aid
 A 16hour course, this is the minimum medical requirement for OA guides. This course is great for anybody heading into the wild, whether its kayaking, rafting, climbing, hiking or skiing. WFA alternates each quarter, between a full weekend or four, 4hour weekday evenings.Certification is good for 3 years. This is offered once a quarter.
 Wilderness First Responder
 Taught by Julie Munger of http://sierrarescue.com/index.php?optioncom_courses&catid4&Itemid58 Sierra Rescue, this is an intensive 9day, 80hour course for the true backcountry adventurer.  Certification is good for 3 years.  Recertification is also available at the same time. This takes place every September before school starts.
 EMT  Emergency Medical Technician
 This is a medically intense course taught by a San Francisco Firefighter/Paramedic and TAd by current UCD students. This course is taught on a semester schedule, twice a year.  Lectures are 4.5 hours in the evening once a week, labs are all day, about one Saturday per month.  This is a very popular course due to its high quality and very reasonable price.  Registration is by lottery system, so make sure to arrive promptly at 9:30 am on the first day of registration. Signups are usually late September and late January. The first semester class runs from October through January, while the second semester course is from February to June.

 Guide Schools
A unique aspect to Outdoor Adventures is that they offer guide schools for each activity. Most university programs either contract out the guiding or have a generic guide school for all of their guides, who then must apprentice each activity. OA is large enough to offer sportspecific training. It can therefore be more thorough and comprehensive. Application will be available at OAs Rental Center in the weeks prior the guide school. Applicants may be called back to do an interview. Finally, the selected few will be allowed to register for this training. After guide school, TAs will need to get CPR and WFA (or WFR/EMT) certified, in addition to working as a TA on a few trips.
 Rock Climbing
 Applications are available at the beginning of the Spring Quarter. This course is two weekends and three classroom sessions.
 Health Care
 Applications are available at the beginning of the Fall Quarter. This course is two classroom sessions and individual testing.
 Baja Sea Kayak
 This is an internal guide school, where OA guides improve their skills by learning how to lead a 10day commercial sea kayaking trip to Baja Mexico. This course takes place over Thanksgiving break and is one classroom session and 11 days in the field.
 Winter
 This is an internal guide school, which means that you must already be a guide for OA. Exceptions are sometimes granted, so definitely talk to someone at OA if youre interested. This course is one weekend and one classroom session.
 Rafting
 Also known as Beginner Guide Training (BGT), this is OAs largest, most intensive guide school. Applications are available in January and due the first week of February. Participants stay at Camp Lotus during spring break, learning to guide the Upper and Lower sections of the South Fork. http://sierrarescue.com/index.php?optioncom_courses&catid1&Itemid58 White Water Rescue certification is included; instruction provided by http://sierrarescue.com/ Sierra Rescue. This school has four classroom sessions and eight river days.
 Wilderness
 Learn the fine art of backpacking, navigation, and going light. Applications are available in mid September and classes start midOctober. This school is composed of two weekends and three classroom sessions.
 Advanced Rafting
 Another internal guide school, Advanced Guide Training (AGT) takes the step up to Class IV rivers. This used to be held on the North Fork of the American and the CalSalmon. In the future, it is likely to be held in April on the Tuolomne River, under an educational permit.
 Sea Kayaking
 The final guide school of the year, this course takes you to Tomales Bay, Bolinas, and SF Bay. It is taught in late May, but applications are available as early as April. This course is two weekends and three classroom sessions.

20120605 19:16:15 nbsp Make sure you doublecheck fees on the website before believing employees on the phone. Users/karinp

20130213 00:43:54 nbsp Does OA offer any backpacking trips designed for more experienced backpackers? Ive been looking for a better way to find fellow backpackers to go on trips with at davis, but Im not sure if this is really the right place to do that Users/ZacMillersmith

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
c Integrator Package for POWHEG
c      subroutine mint(fun,ndim,ncalls0,nitmax,imode,
c ndim=number of dimensions
c ncalls0=# of calls per iteration
c nitmax =# of iterations
c fun(xx,www,ifirst): returns the function to be integrated multiplied by www;
c                     xx(1:ndim) are the variables of integration
c                     ifirst=0: normal behaviour
c imode: integer flag
c
c imode=0:
c When called with imode=0 the routine integrates the absolute value of the function
c and sets up a grid xgrid(0:50,ndim) such that in each ndim-1 dimensional slice
c (i.e. xgrid(m-1,n)<xx(n)<xgrid(m,n)) the contribution of the integral is the same
c the array xgrid is setup at this stage; ans and err are the integral and its error
c
c imode=1 (in fact #0)
c When called with imode=1, the routine performs the integral of the function fun
c using the grid xgrid. If some number in the array ifold, (say, ifold(n))
c is different from 1, it must be a divisor of 50, and the 50 intervals xgrid(0:50,n)
c are grouped into ifold(n) groups, each group containing 50/ifold(n) nearby
c intervals. For example, if ifold(1)=5, the 50 intervals for the first dimension
c are divided in 5 groups of 10. The integral is then performed by folding on top
c of each other these 5 groups. Suppose, for example, that we choose a random point
c in xx(1) = xgrid(2,1)+x*(xgrid(3,1)-xgrid(2,1)), in the group of the first 5 interval.
c we sum the contribution of this point to the contributions of points
c xgrid(2+m*10,1)+x*(xgrid(3+m*10,1)-xgrid(2+m*10,1)), with m=1,...,4.
c In the sequence of calls to the
c function fun, the call for the first point is performed with ifirst=0, and that for
c all subsequent points with ifirst=1, so that the function can avoid to compute
c quantities that only depend upon dimensions that have ifold=1, and do not change
c in each group of folded call. The values returned by fun in a sequence of folded
c calls with ifirst=0 and ifirst=1 are not used. The function itself must accumulate
c the values, and must return them when called with ifirst=2.
c

      subroutine mint(fun,ndim,ncalls0,nitmax,imode,
     #     xgrid,xint,ymax,ans,err,old_grid,accuracy)
c imode=0: integrate and adapt the grid
c imode=1: frozen grid, compute the integral and the upper bounds
c others: same as 1 (for now)
      implicit none
      logical old_grid
      double precision accuracy
      integer nintervals,ndimmax
      parameter (nintervals=50,ndimmax=20)
      integer ncalls0,ndim,nitmax,imode
      real * 8 fun,xgrid(0:nintervals,ndim),xint,ymax(nintervals,ndim),
     #  ans,err
      real * 8 x(ndimmax),vol
      real * 8 xacc(0:nintervals,ndimmax)
      integer icell(ndimmax),ncell(ndimmax)
      integer ifold(ndimmax),kfold(ndimmax)
      common/cifold/ifold
      integer nhits(1:nintervals,ndimmax)
      real * 8 rand(ndimmax)
      real * 8 dx(ndimmax),f,vtot,etot,prod
      integer kdim,kint,kpoint,nit,ncalls,ibin,iret,nintcurr,ifirst
      real * 8 random_mint
      external random_mint,fun
      if(imode.eq.0) then
         do kdim=1,ndim
            ifold(kdim)=1
            if (.not. old_grid) then
              do kint=0,nintervals
                 xgrid(kint,kdim)=dble(kint)/nintervals
              enddo
            endif
         enddo
      elseif(imode.eq.1) then
         do kdim=1,ndim
            nintcurr=nintervals/ifold(kdim)
            if(nintcurr*ifold(kdim).ne.nintervals) then
               write(*,*)
     # 'mint: the values in the ifold array shoud be divisors of',
     #  nintervals
               stop
            endif
            do kint=1,nintcurr
               ymax(kint,kdim)=
     #              xint**(1d0/ndim)
            enddo
         enddo
      endif
      ncalls=ncalls0
      nit=0
      ans=0
      err=0
 10   continue
      nit=nit+1
      if (nit.ge.2.and.ans.eq.0d0) then
         xint=0d0
         return
      endif
      if(nit.gt.nitmax.or.(err/ans).lt.accuracy) then
         if(imode.eq.0) xint=ans
         return
      endif
      if(imode.eq.0) then
         do kdim=1,ndim
            do kint=0,nintervals
               xacc(kint,kdim)=0
               if(kint.gt.0) then
                  nhits(kint,kdim)=0
               endif
            enddo
         enddo
      endif
      vtot=0
      etot=0
      do kpoint=1,ncalls
c find random x, and its random cell
         do kdim=1,ndim
            kfold(kdim)=1
            ncell(kdim)=nintervals/ifold(kdim)*random_mint(25-kdim)+1
            rand(kdim)=random_mint(kdim)
         enddo
         f=0
         ifirst=0
 1       continue
         vol=1
         do kdim=1,ndim
            nintcurr=nintervals/ifold(kdim)
            icell(kdim)=ncell(kdim)+(kfold(kdim)-1)*nintcurr
            ibin=icell(kdim)
            dx(kdim)=xgrid(icell(kdim),kdim)-xgrid(icell(kdim)-1,kdim)
            vol=vol*dx(kdim)*nintcurr
            x(kdim)=xgrid(icell(kdim)-1,kdim)+rand(kdim)*dx(kdim)
            if(imode.eq.0) nhits(ibin,kdim)=nhits(ibin,kdim)+1
         enddo
c contribution to integral
         if(imode.eq.0) then
            f=abs(fun(x,vol,ifirst))+f
         else
c this accumulated value will not be used
            f=fun(x,vol,ifirst)+f
            ifirst=1
            call nextlexi(ndim,ifold,kfold,iret)
            if(iret.eq.0) goto 1
c closing call: accumulated value with correct sign
            f=fun(x,vol,2)
         endif
c
         if(imode.eq.0) then
c accumulate the function in xacc(icell(kdim),kdim) to adjust the grid later
            do kdim=1,ndim
               xacc(icell(kdim),kdim)=xacc(icell(kdim),kdim)+f
            enddo
         else
c update the upper bounding envelope
            prod=1
            do kdim=1,ndim
               prod=prod*ymax(ncell(kdim),kdim)
            enddo
            prod=(f/prod)
            if(prod.gt.1) then
c This guarantees a 10% increase of the upper bound in this cell
               prod=1+0.1d0/ndim
               do kdim=1,ndim
                  ymax(ncell(kdim),kdim)=ymax(ncell(kdim),kdim)
     #          * prod
               enddo
            endif
         endif
         vtot=vtot+f/ncalls
         etot=etot+f**2/ncalls
      enddo
      if(imode.eq.0) then
c iteration is finished; now rearrange the grid
         do kdim=1,ndim
            call regrid(xacc(0,kdim),xgrid(0,kdim),
     #           nhits(1,kdim),nintervals,nit)
         enddo
      endif
c the abs is to avoid tiny negative values
      etot=sqrt(abs(etot-vtot**2)/ncalls)
      write(*,*) vtot,etot
      if(nit.eq.1) then
         ans=vtot
         err=etot
      else
c prevent annoying division by zero for nearly zero
c integrands
         if(etot.eq.0.and.err.eq.0) then
            if(ans.eq.vtot) then
               goto 10
            else
               err=abs(vtot-ans)
               etot=abs(vtot-ans)
            endif
         elseif(etot.eq.0) then
            etot=err
         elseif(err.eq.0) then
            err=etot
         endif
         ans=(ans/err+vtot/etot)/(1/err+1/etot)
         err=1/sqrt(1/err**2+1/etot**2)
      endif
      goto 10
      end

      subroutine regrid(xacc,xgrid,nhits,nint,nit)
      implicit none
      integer  nint,nhits(nint),nit
      real * 8 xacc(0:nint),xgrid(0:nint)
      real * 8 xn(100),r
      integer kint,jint
      do kint=1,nint
c xacc (xerr) already containe a factor equal to the interval size
c Thus the integral of rho is performed by summing up
         if(nhits(kint).ne.0) then
            xacc(kint)= xacc(kint-1)
     #           + abs(xacc(kint))/nhits(kint)
         else
            xacc(kint)=xacc(kint-1)
         endif
      enddo
      do kint=1,nint
         xacc(kint)=xacc(kint)/xacc(nint)
      enddo
c      write(11,*) 'set limits x 0 1 y 0 1'
c      write(11,*) 0, 0
c      do kint=1,nint
c         write(11,*) xgrid(kint),xacc(kint)
c      enddo
c      write(11,*) 'join 0'

      do kint=1,nint
         r=dble(kint)/nint

c         write(11,*) 0, r
c         write(11,*) 1, r
c         write(11,*) ' join'

         do jint=1,nint
            if(r.lt.xacc(jint)) then
               xn(kint)=xgrid(jint-1)+(r-xacc(jint-1))
     #        /(xacc(jint)-xacc(jint-1))*(xgrid(jint)-xgrid(jint-1))
               goto 11
            endif
         enddo
         if(jint.ne.nint+1.and.kint.ne.nint) then
            write(*,*) ' error',jint,nint
            stop
         endif
         xn(nint)=1
 11      continue
      enddo
      do kint=1,nint
         xgrid(kint)=xn(kint)
c         xgrid(kint)=(xn(kint)+2*xgrid(kint))/3
c         xgrid(kint)=(xn(kint)+xgrid(kint)*log(dble(nit)))
c     #        /(log(dble(nit))+1)
c         write(11,*) xgrid(kint), 0
c         write(11,*) xgrid(kint), 1
c         write(11,*) ' join'
      enddo
c      write(11,*) ' newplot'
      end

      subroutine nextlexi(ndim,iii,kkk,iret)
c kkk: array of integers 1 <= kkk(j) <= iii(j), j=1,ndim
c at each call iii is increased lexicographycally.
c for example, starting from ndim=3, kkk=(1,1,1), iii=(2,3,2)
c subsequent calls to nextlexi return
c         kkk(1)      kkk(2)      kkk(3)    iret
c 0 calls   1           1           1       0
c 1         1           1           2       0
c 2         1           2           1       0
c 3         1           2           2       0
c 4         1           3           1       0
c 5         1           3           2       0
c 6         2           1           1       0
c 7         2           1           2       0
c 8         2           2           1       0
c 9         2           2           2       0
c 10        2           3           1       0
c 11        2           3           2       0
c 12        2           3           2       1
      implicit none
      integer ndim,iret,kkk(ndim),iii(ndim)
      integer k
      k=ndim
 1    continue
      if(kkk(k).lt.iii(k)) then
         kkk(k)=kkk(k)+1
         iret=0
         return
      else
         kkk(k)=1
         k=k-1
         if(k.eq.0) then
            iret=1
            return
         endif
         goto 1
      endif
      end

      subroutine gen(fun,ndim,xgrid,ymax,imode,x)
c imode=0 to initialize
c imode=1 to generate
c imode=3 store generation efficiency in x(1)
      implicit none
      integer ndim,imode
      integer nintervals,ndimmax
      parameter (nintervals=50,ndimmax=20)
      real * 8 fun,xgrid(0:nintervals,ndim),
     #         ymax(nintervals,ndim),x(ndim)
      real * 8 dx(ndimmax)
      integer icell(ndimmax),ncell(ndimmax)
      integer ifold(ndimmax),kfold(ndimmax)
      common/cifold/ifold
      real * 8 r,f,ubound,vol,random_mint,xmmm(nintervals,ndimmax)
      real * 8 rand(ndimmax)
      external fun,random_mint
      integer icalls,mcalls,kdim,kint,nintcurr,iret,ifirst
      save icalls,mcalls,xmmm
      if(imode.eq.0) then
         do kdim=1,ndim
            nintcurr=nintervals/ifold(kdim)
            xmmm(1,kdim)=ymax(1,kdim)
            do kint=2,nintcurr
               xmmm(kint,kdim)=xmmm(kint-1,kdim)+
     #              ymax(kint,kdim)
            enddo
            do kint=1,nintcurr
               xmmm(kint,kdim)=xmmm(kint,kdim)/xmmm(nintcurr,kdim)
            enddo
         enddo
         icalls=0
         mcalls=0
         return
      elseif(imode.eq.3) then
         if(icalls.gt.0) then
            x(1)=dble(mcalls)/icalls
         else
            x(1)=-1
         endif
         return
      endif
      mcalls=mcalls+1
 10   continue
      do kdim=1,ndim
         nintcurr=nintervals/ifold(kdim)
         r=random_mint(kdim)
         do kint=1,nintcurr
            if(r.lt.xmmm(kint,kdim)) then
               ncell(kdim)=kint
               goto 1
            endif
         enddo
 1       continue
         rand(kdim)=random_mint(kdim)
      enddo
      ubound=1
      do kdim=1,ndim
         ubound=ubound*ymax(ncell(kdim),kdim)
      enddo
      do kdim=1,ndim
         kfold(kdim)=1
      enddo
      f=0
      ifirst=0
 5    continue
      vol=1
      do kdim=1,ndim
         nintcurr=nintervals/ifold(kdim)
         icell(kdim)=ncell(kdim)+(kfold(kdim)-1)*nintcurr
         dx(kdim)=xgrid(icell(kdim),kdim)-xgrid(icell(kdim)-1,kdim)
         vol=vol*dx(kdim)*nintervals/ifold(kdim)
         x(kdim)=xgrid(icell(kdim)-1,kdim)+rand(kdim)*dx(kdim)
      enddo
      f=f+fun(x,vol,ifirst)
      ifirst=1
      call nextlexi(ndim,ifold,kfold,iret)
      if(iret.eq.0) goto 5
c get final value (x and vol not used in this call)
      f=fun(x,vol,2)
      if(f.lt.0) then
         write(*,*) 'gen: non positive function'
         stop
      endif
      if(f.gt.ubound) then
c         call increasecnt
c     #       ('upper bound failure in inclusive cross section')
       write(*,*) "problem"
      endif
      ubound=ubound*random_mint(1)
      icalls=icalls+1
      if(ubound.gt.f) then
       write(*,*) "problem"
c         call increasecnt
c     #       ('vetoed calls in inclusive cross section')
         goto 10
      endif
      end

      double precision function random_mint(j)
      implicit none
      integer j
      double precision x

      call ntuple(x,0.d0,1.d0,j)
      random_mint=x

      return
      end


text:
module GCL_bc_operator
    use mod_kinds,          only: ik, rk
    use type_operator,      only: operator_t
    use type_chidg_worker,  only: chidg_worker_t
    use type_properties,    only: properties_t
    use DNAD_D
    implicit none

    !>
    !!
    !!  @author Nathan A. Wukie (AFRL)
    !!  @date   9/14/2016
    !!
    !!
    !------------------------------------------------------------------------------------------
    type, public, extends(operator_t) :: GCL_bc_operator_t

    contains

        procedure   :: init
        procedure   :: compute

    end type GCL_bc_operator_t
    !*******************************************************************************************

contains

    !>
    !!
    !!  @author Nathan A. Wukie (AFRL)
    !!  @date   8/10/2017
    !!
    !--------------------------------------------------------------------------------
    subroutine init(self)
        class(GCL_bc_operator_t),   intent(inout) :: self

        !
        ! Set operator name
        !
        call self%set_name('Geometric Conservation BC Operator')

        !
        ! Set operator type
        !
        call self%set_operator_type('BC Diffusive Operator')

        !
        ! Set operator equations
        !
        call self%add_primary_field('g_bar')

    end subroutine init
    !********************************************************************************

    !> Specialized compute routine for Extrapolation Boundary Condition
    !!
    !!  @author Nathan A. Wukie (AFRL)
    !!  @date   8/10/2017
    !!
    !!  @param[in]      mesh    Mesh data containing elements and faces for the domain
    !!  @param[inout]   sdata   Solver data containing solution vector, rhs, linearization, etc.
    !!  @param[inout]   prop    properties_t object containing equations and material_t objects
    !!  @param[in]      face    face_info_t containing indices on location and misc information
    !!  @param[in]      flux    function_into_t containing info on the function being computed
    !!
    !-------------------------------------------------------------------------------------------
    subroutine compute(self,worker,prop)
        class(GCL_bc_operator_t),   intent(inout)   :: self
        type(chidg_worker_t),       intent(inout)   :: worker
        class(properties_t),        intent(inout)   :: prop

        type(AD_D), allocatable, dimension(:)   ::  &
            g_bar, flux_1, flux_2, flux_3, integrand, &
            norm_1, norm_2, norm_3

        real(rk),   allocatable, dimension(:)   ::  &
            det_jacobian_grid

        real(rk),   allocatable, dimension(:,:) :: grid_velocity

        real(rk),   allocatable, dimension(:,:,:) :: &
            inv_jacobian_grid

        ! Just to set up AD variables
        g_bar = worker%get_field('g_bar', 'value', 'boundary')

        !
        ! Get model coefficients
        !
        grid_velocity     = worker%get_grid_velocity_face('face interior')
        det_jacobian_grid = worker%get_det_jacobian_grid_face('value','face interior')
        inv_jacobian_grid = worker%get_inv_jacobian_grid_face('face interior')

        !
        ! Get normal vector
        !
        norm_1  = worker%normal(1)
        norm_2  = worker%normal(2)
        norm_3  = worker%normal(3)

        !
        ! Compute volume flux at quadrature nodes
        !
        flux_1 = g_bar  !just to initialize AD allocation
        flux_2 = g_bar  !just to initialize AD allocation
        flux_3 = g_bar  !just to initialize AD allocation

        flux_1 = (inv_jacobian_grid(:,1,1)*grid_velocity(:,1) + inv_jacobian_grid(:,1,2)*grid_velocity(:,2) + inv_jacobian_grid(:,1,3)*grid_velocity(:,3))*det_jacobian_grid
        flux_2 = (inv_jacobian_grid(:,2,1)*grid_velocity(:,1) + inv_jacobian_grid(:,2,2)*grid_velocity(:,2) + inv_jacobian_grid(:,2,3)*grid_velocity(:,3))*det_jacobian_grid
        flux_3 = (inv_jacobian_grid(:,3,1)*grid_velocity(:,1) + inv_jacobian_grid(:,3,2)*grid_velocity(:,2) + inv_jacobian_grid(:,3,3)*grid_velocity(:,3))*det_jacobian_grid

        !
        ! Integrate volume flux
        !
        call worker%integrate_boundary_condition('g_bar','Advection',flux_1,flux_2,flux_3)

    end subroutine compute
    !**********************************************************************************************

end module GCL_bc_operator

max_stars_repo_name: nwukie/ChiDG

max_issues_repo_name: nwukie/ChiDG

max_forks_repo_name: nwukie/ChiDG