text:
*----------------------------------------------------------------------*
      subroutine set_orbinf(orb_info,env_type,hole_rv)
*----------------------------------------------------------------------*
*     set up orbital info arrays
*
*     symmetry ordering:
*        orbitals are ordered according to symmetry only
*
*     type ordering:
*        orbitals are ordered first according to the shell they
*        belong to (and within shell according to their symmetry)
*
*----------------------------------------------------------------------*

      implicit none

      include 'stdunit.h'
      include 'opdim.h'
      include 'def_orbinf.h'

      integer, parameter ::
     &     ntest = 00

      type(orbinf), intent(inout) ::
     &     orb_info
      character, intent(in) ::
     &     env_type*(*)
      logical, intent(in) ::
     &     hole_rv

      integer ::
     &     ngas, nsym, ntoob, nspin,
     &     iprint, idx, jdx, kdx, isym, igas, igastp, ispin, ipass,
     &     idxst, idxnd, ist, ind, inc, igasr, j, caborb, iloop, loop
      integer ::
     &     icount(ngastp)
      integer, allocatable ::
     &     iadscr(:),koffs(:)
      logical, pointer ::
     &     assigned(:)

      iprint = max(ntest,iprlvl)

      if (iprint.ge.100) then
        write(lulog,*) '************'
        write(lulog,*) ' set_orbinf'
        write(lulog,*) '************'
        write(lulog,*) ' hole_rv: ',hole_rv
        write(lulog,*) ' nspin  = ',orb_info%nspin
        write(lulog,*) ' ngas   = ',orb_info%ngas
        write(lulog,*) ' nsym   = ',orb_info%nsym
        write(lulog,*) ' ntoob  = ',orb_info%ntoob
        write(lulog,*) ' caborb = ',orb_info%caborb
      end if

      ! for convenience
      nspin = orb_info%nspin
      ngas = orb_info%ngas
      nsym = orb_info%nsym
      ntoob = orb_info%ntoob
      caborb=orb_info%caborb

      if (hole_rv.and.nspin.gt.1) then
        call quit(1,'set_orbinf',
     &       'hole_rv is not compatible with nspin.gt.1')
      end if

      ! allocate some arrays on orb_info structure
      allocate(orb_info%ireots(ntoob+caborb),
     &     orb_info%ireost(ntoob+caborb),orb_info%igamorb(ntoob+caborb),
     &     orb_info%igasorb(ntoob+caborb),orb_info%mostnd(2,nsym,ngas),
     &     orb_info%ngas_hpv(ngastp),orb_info%nactt_hpv(ngastp),
     &     orb_info%idx_gas(ngastp),orb_info%ioff_gas(ngastp),
     &     orb_info%gas_reo(ngas),orb_info%norb_hpv(ngastp,nspin))

      ! set the info arrays:

      ! set ngas_hpv (number of spaces per H/P/V space)
      ! set nactt_hpv (number of active orbitals per H/P/V)
      ! open shell with p/h splitting: we double-count!
      orb_info%ngas_hpv(1:ngastp) = 0
      orb_info%nactt_hpv(1:ngastp) = 0
      orb_info%norb_hpv(1:ngastp,1:nspin) = 0
      do ispin = 1, nspin
        icount(1:ngastp) = 0
        do igas = 1, ngas
          icount(orb_info%ihpvgas(igas,ispin)) =
     &         icount(orb_info%ihpvgas(igas,ispin))+1
        end do
        do igastp = 1, ngastp
          orb_info%ngas_hpv(igastp) =
     &         max(orb_info%ngas_hpv(igastp),icount(igastp))
        end do

        icount(1:ngastp) = 0
        do igas = 1, ngas
          if (orb_info%iad_gas(igas).ne.2) cycle
          do isym = 1, nsym
            icount(orb_info%ihpvgas(igas,ispin)) =
     &           icount(orb_info%ihpvgas(igas,ispin))+
     &           orb_info%igassh(isym,igas)
          end do
        end do
        do igastp = 1, ngastp
          orb_info%nactt_hpv(igastp) =
     &         max(orb_info%nactt_hpv(igastp),icount(igastp))
        end do

        do igas = 1, ngas
          do isym = 1, nsym
            orb_info%norb_hpv(orb_info%ihpvgas(igas,ispin),ispin) =
     &           orb_info%norb_hpv(orb_info%ihpvgas(igas,ispin),ispin)+
     &           orb_info%igassh(isym,igas)
          end do
        end do

      end do

      ! index and offset array from ngas_hpv:
c      idx = 1
cmh adapted for use of valence space
      orb_info%ioff_gas(1:ngas) = 0
      orb_info%idx_gas(1:ngas) = 0
      do igastp = 1, ngastp
        do idx = ngas,1,-1
          orb_info%ioff_gas(orb_info%ihpvgas(idx,1)) = idx-1
          orb_info%idx_gas(orb_info%ihpvgas(idx,1)) = idx
        end do
c        orb_info%ioff_gas(igastp) = idx-1
c        orb_info%idx_gas(igastp) = idx
c        idx = idx+orb_info%ngas_hpv(igastp)
      end do

      if (iprint.ge.100) then
        write(lulog,*) 'ihpvgas:   ',orb_info%ihpvgas(1:ngas,1)
        if (orb_info%nspin.gt.1)
     &       write(lulog,*) '           ',orb_info%ihpvgas(1:ngas,2)
        write(lulog,*) 'ngas_hpv:  ',orb_info%ngas_hpv(1:ngastp)
        write(lulog,*) 'nactt_hpv: ',orb_info%nactt_hpv(1:ngastp)
        write(lulog,*) 'ioff_gas:  ',orb_info%ioff_gas(1:ngastp)
        write(lulog,*) 'idx_gas:   ',orb_info%idx_gas(1:ngastp)
        write(lulog,*) 'igassh:'
        do igas = 1, ngas
          write(lulog,'(3x,8i4)') orb_info%igassh(1:nsym,igas)
        end do
      end if

      do igas = 1, ngas
        orb_info%gas_reo(igas) = igas
      end do
      if (hole_rv) then
        igasr = 0
        do igas = orb_info%ngas_hpv(ihole), 1, -1
          igasr = igasr + 1
          orb_info%gas_reo(igas) = igasr
        end do
      end if

      ! set up igamorb (IRREP per orbital in type ordering) and
      ! set up igasorb (shell per orbital in type ordering)
      ! set up mostnd (start and end indices)
      ! start with hole spaces
      if (hole_rv) then
        ! reverse counting
        ist = orb_info%ngas_hpv(ihole)
        ind = 1
        inc = -1
        call quit(1,'set_orbinf','who needs hole_rv?? (1)')
      else
        ist = 1
        ind = orb_info%ngas_hpv(ihole)
        inc = +1
      end if
      idxst = 1
      igasr = 0
      do igas = ist, ind, inc
        igasr = igasr+1
        do isym = 1, nsym
          idxnd = idxst+orb_info%igassh(isym,igas)-1
          orb_info%mostnd(1:2,isym,igasr) = (/idxst,idxnd/)
          if (idxst.le.idxnd) then
            orb_info%igamorb(idxst:idxnd) = isym
            orb_info%igasorb(idxst:idxnd) = igasr
          end if
          idxst = idxnd+1
        end do
      end do
      ! remaining spaces
      ist = ind+1
      ind = ngas
      do igas = ist, ind
        do isym = 1, nsym
          idxnd = idxst+orb_info%igassh(isym,igas)-1
          orb_info%mostnd(1:2,isym,igas) = (/idxst,idxnd/)
          if (idxst.le.idxnd) then
            orb_info%igamorb(idxst:idxnd) = isym
            orb_info%igasorb(idxst:idxnd) = igas
          end if
          idxst = idxnd+1
        end do
      end do

      if (iprint.ge.100) then
        write(lulog,*) 'gas_reo:',orb_info%gas_reo(1:ngas)
        write(lulog,*) 'igamorb:'
        call iwrtma(orb_info%igamorb,1,ntoob+caborb,1,ntoob+caborb)
        write(lulog,*) 'igasorb:'
        call iwrtma(orb_info%igasorb,1,ntoob+caborb,1,ntoob+caborb)
        write(lulog,*) 'mostnd:'
        do igas = 1, ngas
          write(lulog,'(2x,i4,2x,8(x,2i4))')
     &         igas,orb_info%mostnd(1:2,1:nsym,igas)
        end do
      end if

      ! generate symmetry ordering -> type ordering mapping
      select case(trim(env_type))
      case('dalton','DALTON','dalton_special','DALTON_SPECIAL',
     &     'dalton64','DALTON64','molpro_dump','MOLPRO_DUMP')
        jdx = 0
        do isym = 1, nsym
          do igas = 1, ngas
            ! ignore cabs orbitals
            if(caborb.gt.0.and.igas.eq.ngas)cycle
            if (orb_info%ihpvgas(igas,1).eq.ihole.and.hole_rv) then
              igasr = orb_info%ngas_hpv(ihole)-igas+1
              ist = orb_info%mostnd(2,isym,igasr)
              ind = orb_info%mostnd(1,isym,igasr)
              inc = -1
            else
              ist = orb_info%mostnd(1,isym,igas)
              ind = orb_info%mostnd(2,isym,igas)
              inc = +1
            end if
            do idx = ist, ind, inc
              jdx = jdx+1
              orb_info%ireost(jdx) = idx
            end do
          end do
        end do

        ! append CABS orbitals (no reordering)
        do idx = ntoob+1, ntoob+caborb
          orb_info%ireost(idx) = idx
        end do
      case ('molpro_ifc','MOLPRO_IFC','virtual','VIRTUAL')
        ! the interface removes the core orbitals and shifts the numbering
        ! we just pretend that the symmetry ordering is done for core and other
        ! orbitals separately
        ! the "virtual molecule" is just handled analogous to that
        ! (although we will never do actual calculations with it)
        if (hole_rv) call quit(1,'set_orbinf','who needs hole_rv??')
        jdx = 0
        do ipass = 1,2 ! two passes: 1-deleted core, 2-others
          do isym = 1, nsym
            do igas = 1, ngas
              if (ipass.eq.1.and.orb_info%iad_gas(igas).ne.1) cycle
              if (ipass.eq.2.and.orb_info%iad_gas(igas).eq.1) cycle
              ! ignore cabs orbitals
              if(caborb.gt.0.and.igas.eq.ngas)cycle
              ist = orb_info%mostnd(1,isym,igas)
              ind = orb_info%mostnd(2,isym,igas)
              inc = +1
              do idx = ist, ind, inc
                jdx = jdx+1
                orb_info%ireost(jdx) = idx
              end do
            end do
          end do
        end do

        ! append CABS orbitals (no reordering)
        do idx = ntoob+1, ntoob+caborb
          orb_info%ireost(idx) = idx
        end do
      case ('cfour','CFOUR')
        ! actually: symmetry ordering means "the way the external program
        !   has ordered the orbitals"; cfour orders kind-of-type-ordering like
        !   we try to get this info by interpreting the ext_gamorb array
        call make_ext2typ_reo(orb_info)
      case ('gamess','GAMESS')
        ! (We loop over all orbitals to avoid errors when the active
        !  space is modified such that orbitals must be reordered.
        !  The DALTON inferface above also needs to be changed
        !  if we want to use this nasty kind of tricks there as well)
        if (orb_info%n_bound_orbs.ne.ntoob.or.caborb.ne.0)
     &        call quit(1,'set_orbinf','need full orb-sym list!')
        allocate(assigned(ntoob))
        assigned = .false.
        do igas = 1, ngas
          do isym = 1, nsym
            if (orb_info%igassh(isym,igas).eq.0) cycle
            jdx = orb_info%mostnd(1,isym,igas)
            do idx = 1, ntoob
              if (assigned(idx)) cycle
              if (orb_info%isym_bound_orbs(idx).eq.isym) then
                orb_info%ireost(idx) = jdx
                jdx = jdx + 1
                assigned(idx) = .true.
                if (jdx.gt.orb_info%mostnd(2,isym,igas)) exit
              end if
            end do
          end do
        end do
        if (.not.all(assigned(1:ntoob)))
     &    call quit(1,'set_orbinf','unable to assign all orbitals')
        deallocate(assigned)
      case default
        call quit(1,'set_orbinf','unknown type '//trim(env_type))
      end select

        ! generate reverse mapping
        do idx = 1, ntoob+caborb
          orb_info%ireots(orb_info%ireost(idx)) = idx
        end do

      if (iprint.ge.100) then
        write(lulog,*) 'ireost:'
        call iwrtma(orb_info%ireost,1,ntoob,1,ntoob+caborb)
        write(lulog,*) 'ireots:'
        call iwrtma(orb_info%ireots,1,ntoob,1,ntoob+caborb)
      end if

      ! reverse iad_array
      if (hole_rv) then
        allocate(iadscr(ngas))
        iadscr(1:ngas) = orb_info%iad_gas(1:ngas)
        do igas = 1, orb_info%ngas_hpv(ihole)
          orb_info%iad_gas(igas) =
     &         iadscr(orb_info%ngas_hpv(ihole)-igas+1)
        end do
        deallocate(iadscr)
      end if

      if (iprint.ge.100) then
        write(lulog,*) 'iad_gas:'
        call iwrtma(orb_info%iad_gas,ngas,1,ngas,1)
      end if

      return
      end

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public


text:
!>----------------------------------------
!! Simple program to test and time sorting routines
!!
!! Generates n random numbers, uses heap_sort to sort them,
!!  Then tests that the sort succeeded and prints timing info
!!
!! Next this is repeated for quick_sort.
!!
!! Initial data show that this implementation of quick_sort is ~2x faster than
!!  the (somewhat naive) implementation of heap_sort
!!
!!----------------------------------------
program test_sort
    use sort_mod

    implicit none

    real, parameter :: valid_error = 1e-5
    integer*8, parameter :: n = 1000000

    real, allocatable, dimension(:) :: to_sort, sorted
    integer :: i, start_time, end_time, COUNT_RATE, COUNT_MAX

    allocate(to_sort(n))
    allocate(sorted(n))

    !-------------------------------
    ! SETUP sorting
    !-------------------------------
    call random_number(to_sort)
    sorted=0
    ! test heap_sort (system_clock is used to time how long it takes)

    call system_clock(start_time)
    !-------------------------------
    ! HEAP SORT
    !-------------------------------
    call heap_sort(to_sort, sorted)

    call system_clock(end_time, COUNT_RATE, COUNT_MAX)
    if (start_time>end_time) end_time=end_time+COUNT_MAX

    call show_results(sorted, "HEAP_SORT", (end_time-start_time) / real(COUNT_RATE))

    !-------------------------------
    ! SETUP sorting
    !-------------------------------
    ! now test quick_sort
    call random_number(to_sort)
    sorted = 0

    call system_clock(start_time)
    !-------------------------------
    ! QUICK SORT
    !-------------------------------
    call quick_sort(to_sort, sorted)

    call system_clock(end_time, COUNT_RATE, COUNT_MAX)
    if (start_time>end_time) end_time=end_time+COUNT_MAX

    call show_results(sorted, "QUICK_SORT", (end_time-start_time) / real(COUNT_RATE))

    !-------------------------------
    ! SETUP sorting
    !-------------------------------
    ! now test quick_sort
    call random_number(to_sort)

    call system_clock(start_time)
    !-------------------------------
    ! QUICK SORT in place
    !-------------------------------
    call sort(to_sort)

    call system_clock(end_time, COUNT_RATE, COUNT_MAX)
    if (start_time>end_time) end_time=end_time+COUNT_MAX

    call show_results(to_sort, "QUICK_SORT (in place)", (end_time-start_time) / real(COUNT_RATE))

contains
    !>--------------------------
    !! Print out results from a given sort test
    !!
    !!--------------------------
    subroutine show_results(data, name, time)
        implicit none
        real, dimension(:), intent(in) :: data
        character(len=*),   intent(in) :: name
        real,               intent(in) :: time

        integer :: i, n, err

        n = size(data)

        print*, "---------------------------"
        err = 0
        do i=1,n-1
            if (sorted(i+1)<sorted(i)) then
                err = 1
            endif
        end do
        print*, trim(name)
        if ( err < valid_error ) then
            print*, "  PASSED"
        else
            print*, "  FAILED:", err
        endif
        print*, "---------------------------"
        print*, "   Elapsed time = ", time

    end subroutine show_results

end program test_sort

max_stars_repo_name: gutmann/downscale

max_issues_repo_name: gutmann/downscale

max_forks_repo_name: gutmann/downscale


! { dg-do compile }

!
! PR fortran/64726
!
subroutine oacc1()
  implicit none
  integer :: i
  integer  :: a
  !$acc parallel loop reduction(+:a)
  do i = 1,5
  enddo
  !$acc end parallel loop
  !$acc kernels loop collapse(2)
  do i = 2,6
     do a = 3,5
     enddo
  enddo
  !$acc end kernels loop
end subroutine oacc1


text:
      function diskdenspsi(r,z,psi)
      parameter (pi=3.1415926535)
      common /gparameters/  a, b, c, v0, q, psi00,
     +                      psiout, rho1, sigbulge2,
     +                      rmdisk, rdisk, zdisk, outdisk, drtrunc,
     +                      potcor1
      common /moreconstants/ v02, v03, rdisk2, diskconst, bulgea

c the truncation factor is derived assuming that the azimuthal velocity dist
c is gaussian with dispersion sigma_R * 2 omega/kappa, as per epicycle theory.
      if( abs(z/zdisk) .gt. 10.) goto 99
      erfarg=(r-outdisk)/1.4142136/drtrunc
      if (erfarg.gt.4.) goto 99
      if (erfarg.lt.-4) then
         trunfac=1
      else
         trunfac=0.5*erfc(erfarg)
      endif
      if (z.eq.0.) then
         con=1
      else
        psizh= pot(r,zdisk)
        psi0 = pot(r,0.)
        con = 0.419974**((psi-psi0)/(psizh-psi0))
      endif
      diskdenspsi = diskconst*exp(-r/rdisk)*con*trunfac
      return
 99   diskdenspsi = 0.0
      return
      end

ext: f

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
!==============================================================================!
  subroutine Comm_Mod_Create_New_Types(grid)
!------------------------------------------------------------------------------!
!   Creates new data type for MPI I/O.                                         !
!------------------------------------------------------------------------------!
!----------------------------------[Modules]-----------------------------------!
  use Grid_Mod
!------------------------------------------------------------------------------!
  implicit none
!---------------------------------[Arguments]----------------------------------!
  type(Grid_Type) :: grid
!-----------------------------------[Locals]-----------------------------------!
  integer :: error
!==============================================================================!
!   There is an issue with this procedure, but it's more related to MPI/IO     !
!   functions than T-Flows.  In cases a subdomain has no physical boundary     !
!   cells, variable "nb_s" turns out to be zero.  This, per se, should not     !
!   be an issue if MPI/IO functions could handle calls to:                     !
!   "Mpi_Type_Create_Indexed_Block(...)" and "Mpi_File_Write(...)" with zero   !
!   length.  But they don't.  Therefore, I avoid allocation with zero size     !
!   (max(nb_s,1)) here and creation of new types with zero size in             !
!   "Load_Maps".  It is a bit of a dirty trick :-(                             !
!------------------------------------------------------------------------------!

  ! Create new type for cells
  call Mpi_Type_Create_Indexed_Block(nc_s,                   & ! length of map
                                     1,                      & ! block size
                                     grid % comm % cell_map, & ! displacements
                                     MPI_DOUBLE,             & ! old data type
                                     cell_map_type,          & ! new data type
                                     error)                    ! integer error
  call Mpi_Type_Commit(cell_map_type, error)

  ! Create new type for boundary cells
  call Mpi_Type_Create_Indexed_Block(max(nb_s,1),                & ! map length
                                     1,                          & ! block size
                                     grid % comm % bnd_cell_map, & ! displacem.
                                     MPI_DOUBLE,                 & ! old type
                                     bnd_cell_map_type,          & ! new type
                                     error)                        ! int. error
  call Mpi_Type_Commit(bnd_cell_map_type, error)

  end subroutine

max_stars_repo_name: MassimoBorrelliPhysicist/T-Flows

max_issues_repo_name: MassimoBorrelliPhysicist/T-Flows

max_forks_repo_name: MassimoBorrelliPhysicist/T-Flows


text:
namespace AulaSoftwareLibre\Gata\Domain\Idea\Model {
    data IdeaId = IdeaId deriving (Uuid);
    data IdeaTitle = String deriving (ToString, FromString, Equals);
    data IdeaDescription = String deriving (ToString, FromString, Equals);
    data IdeaStatus = Pending | Accepted | Rejected deriving (Enum);
}

namespace AulaSoftwareLibre\Gata\Domain\Idea\Event {
    data IdeaAttendeeWasRegistered = IdeaAttendeeWasRegistered {
        \AulaSoftwareLibre\Gata\Domain\Idea\Model\IdeaId $ideaId,
        \AulaSoftwareLibre\Gata\Domain\User\Model\UserId $userId,
    } deriving (AggregateChanged);

    data IdeaAttendeeWasUnregistered = IdeaAttendeeWasUnregistered {
        \AulaSoftwareLibre\Gata\Domain\Idea\Model\IdeaId $ideaId,
        \AulaSoftwareLibre\Gata\Domain\User\Model\UserId $userId,
    } deriving (AggregateChanged);

    data IdeaCapacityWasLimited = IdeaCapacityWasLimited {
        \AulaSoftwareLibre\Gata\Domain\Idea\Model\IdeaId $ideaId,
        int $limit
    } deriving (AggregateChanged);

    data IdeaCapacityWasUnLimited = IdeaCapacityWasUnlimited {
        \AulaSoftwareLibre\Gata\Domain\Idea\Model\IdeaId $ideaId
    } deriving (AggregateChanged);

    data IdeaWasAccepted = IdeaWasAccepted {
        \AulaSoftwareLibre\Gata\Domain\Idea\Model\IdeaId $ideaId
    } deriving (AggregateChanged);

    data IdeaWasAdded = IdeaWasAdded {
        \AulaSoftwareLibre\Gata\Domain\Idea\Model\IdeaId $ideaId,
        \AulaSoftwareLibre\Gata\Domain\Group\Model\GroupId $groupId,
        \AulaSoftwareLibre\Gata\Domain\Idea\Model\IdeaTitle $title,
        \AulaSoftwareLibre\Gata\Domain\Idea\Model\IdeaDescription $description
    } deriving (AggregateChanged);

    data IdeaWasRedescribed = IdeaWasRedescribed {
        \AulaSoftwareLibre\Gata\Domain\Idea\Model\IdeaId $ideaId,
        \AulaSoftwareLibre\Gata\Domain\Idea\Model\IdeaDescription $description
    } deriving (AggregateChanged);

    data IdeaWasRejected = IdeaWasRejected {
        \AulaSoftwareLibre\Gata\Domain\Idea\Model\IdeaId $ideaId
    } deriving (AggregateChanged);

    data IdeaWasRetitled = IdeaWasRetitled {
        \AulaSoftwareLibre\Gata\Domain\Idea\Model\IdeaId $ideaId,
        \AulaSoftwareLibre\Gata\Domain\Idea\Model\IdeaTitle $title
    } deriving (AggregateChanged);
}

namespace AulaSoftwareLibre\Gata\Application\Idea\Command {
    data AcceptIdea = AcceptIdea {
        \AulaSoftwareLibre\Gata\Domain\Idea\Model\IdeaId $ideaId,
    } deriving (Command);

    data AddIdea = AddIdea {
        \AulaSoftwareLibre\Gata\Domain\Idea\Model\IdeaId $ideaId,
        \AulaSoftwareLibre\Gata\Domain\Group\Model\GroupId $groupId,
        \AulaSoftwareLibre\Gata\Domain\Idea\Model\IdeaTitle $title,
        \AulaSoftwareLibre\Gata\Domain\Idea\Model\IdeaDescription $description
    } deriving (Command);

    data RedescribeIdea = RedescribeIdea {
        \AulaSoftwareLibre\Gata\Domain\Idea\Model\IdeaId $ideaId,
        \AulaSoftwareLibre\Gata\Domain\Idea\Model\IdeaDescription $description
    } deriving (Command);

    data RegisterIdeaAttendee = RegisterIdeaAttendee {
        \AulaSoftwareLibre\Gata\Domain\Idea\Model\IdeaId $ideaId,
        \AulaSoftwareLibre\Gata\Domain\User\Model\UserId $userId,
    } deriving (Command);

    data RejectIdea = RejectIdea {
        \AulaSoftwareLibre\Gata\Domain\Idea\Model\IdeaId $ideaId,
    } deriving (Command);

    data RetitleIdea = RetitleIdea {
        \AulaSoftwareLibre\Gata\Domain\Idea\Model\IdeaId $ideaId,
        \AulaSoftwareLibre\Gata\Domain\Idea\Model\IdeaTitle $title,
    } deriving (Command);

    data UnregisterIdeaAttendee = UnregisterIdeaAttendee {
        \AulaSoftwareLibre\Gata\Domain\Idea\Model\IdeaId $ideaId,
        \AulaSoftwareLibre\Gata\Domain\User\Model\UserId $userId,
    } deriving (Command);
}

ext: fpp

max_stars_repo_name: aulasoftwarelibre/gata-legacy

max_issues_repo_name: aulasoftwarelibre/gata

max_forks_repo_name: aulasoftwarelibre/gata-legacy


text:
*----------------------------------------------------------------------*
      subroutine update_graphs(str_info,mel,orb_info)
*----------------------------------------------------------------------*
*
*     successor of unique_graph:
*
*     scan through occupation classes of operator mel%op assiciated
*     with ME-list mel and look for new graphs (i.e. occupation class
*     for H/P/V alone); set up the list on string info structure str_info
*
*     the idx_graph array of mel%op is set up accordingly
*
*     andreas, autumn 2006
*              december 2007
*
*----------------------------------------------------------------------*

      implicit none

      include 'stdunit.h'
      include 'opdim.h'
      include 'def_filinf.h'
      include 'def_operator.h'
      include 'def_me_list.h'
      include 'def_graph.h'
      include 'def_strinf.h'
      include 'def_orbinf.h'

      integer, parameter ::
     &     ntest = 00

      type(strinf), intent(inout) ::
     &     str_info
      type(orbinf), intent(in), target ::
     &     orb_info
      type(me_list), intent(inout) ::
     &     mel

      type(operator), pointer ::
     &     op
      logical ::
     &     unique, same
      integer ::
     &     ngas, nspin,
     &     iocc_cls, ihpv, ica, igas, jgas, idx, idxgr, len(ngastp),
     &     nocc, njoined, ijoin, ioff, igraph
      integer, pointer ::
     &     ihpvgas(:,:)

      logical, external ::
     &     restr_cmp

      if (ntest.ge.100) then
        call write_title(lulog,wst_dbg_subr,'update_graphs')
        write(lulog,*) ' unique graphs on input: ',str_info%ngraph
      end if

      nspin = orb_info%nspin
      ngas = orb_info%ngas
      ihpvgas => orb_info%ihpvgas

      op => mel%op

      njoined = op%njoined

      ! loop over occupation classes of op
      do iocc_cls = 1, op%n_occ_cls
       if(op%formal_blk(iocc_cls))cycle

       ioff = (iocc_cls-1)*njoined
       do ijoin = 1, njoined
        idx = ioff+ijoin

        ! do not forget to init
        mel%idx_graph(1:ngastp,1:2,idx) = 0

        ! loop over h/p/v/x and c/a
        do ihpv = 1,ngastp
          ica_loop: do ica = 1,2

            nocc = op%ihpvca_occ(ihpv,ica,idx)
            if (nocc.eq.0) cycle ica_loop

            unique = .true.
            ! loop over all graphs that are already defined
            igraph_loop: do igraph = 1, str_info%ngraph
c dbg
c            print *,'look at graph # ',igraph,
c     &           str_info%ispc_typ(igraph),
c     &           str_info%ispc_occ(igraph)
c dbg

              ! we might compare number of orbitals instead of ihpv
              if (ihpv.eq.str_info%ispc_typ(igraph).and.
     &             nocc.eq.str_info%ispc_occ(igraph)) then

c dbg
c                print *,'comparing: ',ihpv, nocc
c dbg
                ! ... but are restrictions (=shape of graph) the same?
                same = restr_cmp(op%igasca_restr(1,1,1,1,1,idx),
     &                       str_info%igas_restr(1,1,1,1,igraph),
     &                       ica,ihpv,
     &                       ihpvgas,ngas,nspin)

                unique = .not.same

              end if
c dbg
c              print *,'still unique? ',unique
c dbg

              ! .not.unique == we have this graph already
              ! let op%idx_graph point to this graph ...
              if (.not.unique) mel%idx_graph(ihpv,ica,idx) = igraph
              ! ... and go to next space
              if (.not.unique) cycle ica_loop

            end do igraph_loop
c dbg
c              print *,'final: unique = ',unique
c dbg

            ! a new unique graph was found:
            call add_graph(ihpv,nocc,ica,
     &           op%igasca_restr(1,1,1,1,1,idx),
     &           str_info,orb_info)

            ! set idxgr to that graph
            idxgr = str_info%ngraph
            mel%idx_graph(ihpv,ica,idx) = idxgr

          end do ica_loop
        end do

       end do ! ijoin
      end do ! iocc_cls

      if (ntest.eq.100) then
        write(lulog,*) 'unique graphs on output: ', str_info%ngraph
        len(1:ngastp) = 0
        do igas = 1, ngas
          len(ihpvgas(igas,1)) = len(ihpvgas(igas,1))+1
        end do
        write(lulog,*) '    #    typ occ'
        write(lulog,*) '                 restrictions'
        do igraph = 1, str_info%ngraph
          write(lulog,'(2x,i3,2x,2i4)') igraph,
     &         str_info%ispc_typ(igraph),
     &         str_info%ispc_occ(igraph)
          write(lulog,'(15x,10(2i3,x))')
     &         str_info%igas_restr(1:2,
     &                  1:len(str_info%ispc_typ(igraph)),1,1,igraph)
          if (nspin.eq.2)
     &         write(lulog,'(15x,10(2i3,x))')
     &         str_info%igas_restr(1:2,
     &                  1:len(str_info%ispc_typ(igraph)),1,2,igraph)
        end do
        write(lulog,*) 'operator->graph assignments:'
        do iocc_cls = 1, op%n_occ_cls
          if(op%formal_blk(iocc_cls))cycle
          ioff = (iocc_cls-1)*njoined
          write(lulog,'(2x,i3,4x,5(4i3,2x))') iocc_cls,
     &         mel%idx_graph(1:ngastp,1,ioff+1:ioff+njoined)
          write(lulog,'(5x,4x,5(4i3))')
     &         mel%idx_graph(1:ngastp,2,ioff+1:ioff+njoined)
        end do
      end if

      return
      end


text:
!**************************************************************
!* AceGen    6.702 Windows (4 May 16)                         *
!*           Co. J. Korelc  2013           19 Mar 17 15:23:33 *
!**************************************************************
! User     : Full professional version
! Notebook : ResidualFunction
! Evaluation time                 : 15 s    Mode  : Optimal
! Number of formulae              : 240     Method: Automatic
! Subroutine                      : RF4 size: 5790
! Total size of Mathematica  code : 5790 subexpressions
! Total size of Fortran code      : 12570 bytes

!******************* S U B R O U T I N E **********************
SUBROUTINE RF4(v,x,props,statev,Fnew,R)
USE SMSUtility
IMPLICIT NONE
LOGICAL b267
DOUBLE PRECISION v(489),x(46),props(18),statev(46),Fnew(9),R(46)
v(439)=props(4)*x(1)
v(428)=x(39)*x(40)-x(42)*x(46)
v(427)=-(x(40)*x(41))+x(44)*x(46)
v(426)=x(41)*x(43)-x(38)*x(46)
v(425)=-(x(39)*x(43))+x(45)*x(46)
v(424)=x(42)*x(43)-x(40)*x(45)
v(423)=-(x(38)*x(42))+x(44)*x(45)
v(422)=x(38)*x(39)-x(41)*x(45)
v(421)=x(38)*x(40)-x(43)*x(44)
v(420)=x(41)*x(42)-x(39)*x(44)
v(419)=x(30)*x(31)-x(33)*x(37)
v(418)=-(x(31)*x(32))+x(35)*x(37)
v(417)=x(32)*x(34)-x(29)*x(37)
v(416)=-(x(30)*x(34))+x(36)*x(37)
v(415)=x(33)*x(34)-x(31)*x(36)
v(414)=-(x(29)*x(33))+x(35)*x(36)
v(413)=x(29)*x(30)-x(32)*x(36)
v(412)=x(29)*x(31)-x(34)*x(35)
v(411)=x(32)*x(33)-x(30)*x(35)
v(410)=x(21)*x(22)-x(24)*x(28)
v(409)=-(x(22)*x(23))+x(26)*x(28)
v(408)=x(23)*x(25)-x(20)*x(28)
v(407)=-(x(21)*x(25))+x(27)*x(28)
v(406)=x(24)*x(25)-x(22)*x(27)
v(405)=-(x(20)*x(24))+x(26)*x(27)
v(404)=x(20)*x(21)-x(23)*x(27)
v(403)=x(20)*x(22)-x(25)*x(26)
v(402)=x(23)*x(24)-x(21)*x(26)
v(401)=x(12)*x(13)-x(15)*x(19)
v(400)=-(x(13)*x(14))+x(17)*x(19)
v(399)=x(14)*x(16)-x(11)*x(19)
v(398)=-(x(12)*x(16))+x(18)*x(19)
v(397)=x(15)*x(16)-x(13)*x(18)
v(396)=-(x(11)*x(15))+x(17)*x(18)
v(395)=x(11)*x(12)-x(14)*x(18)
v(394)=x(11)*x(13)-x(16)*x(17)
v(393)=x(14)*x(15)-x(12)*x(17)
v(392)=x(3)*x(4)-x(10)*x(6)
v(391)=-(x(3)*x(7))+x(10)*x(9)
v(390)=-(x(10)*x(2))+x(5)*x(7)
v(389)=-(x(4)*x(5))+x(10)*x(8)
v(388)=x(6)*x(7)-x(4)*x(9)
v(387)=x(2)*x(3)-x(5)*x(9)
v(386)=-(x(2)*x(6))+x(8)*x(9)
v(385)=x(2)*x(4)-x(7)*x(8)
v(384)=x(5)*x(6)-x(3)*x(8)
v(383)=-statev(1)+x(1)
v(382)=0.15d1*props(6)
v(381)=props(8)*v(382)
v(380)=1d0-props(6)
v(329)=props(11)*v(382)
v(347)=props(14)*v(382)
v(365)=props(17)*v(382)
v(121)=1d0/(v(386)*x(10)+v(387)*x(4)+v(384)*x(7))
v(306)=v(395)*x(13)+v(393)*x(16)+v(396)*x(19)
v(457)=-(v(396)/v(306))
v(456)=-(v(395)/v(306))
v(455)=-(v(393)/v(306))
v(454)=-(v(400)/v(306))
v(453)=-(v(399)/v(306))
v(452)=-(v(394)/v(306))
v(451)=-(v(401)/v(306))
v(450)=-(v(398)/v(306))
v(449)=-(v(397)/v(306))
v(324)=v(404)*x(22)+v(402)*x(25)+v(405)*x(28)
v(466)=-(v(405)/v(324))
v(465)=-(v(404)/v(324))
v(464)=-(v(402)/v(324))
v(463)=-(v(409)/v(324))
v(462)=-(v(408)/v(324))
v(461)=-(v(403)/v(324))
v(460)=-(v(410)/v(324))
v(459)=-(v(407)/v(324))
v(458)=-(v(406)/v(324))
v(343)=v(413)*x(31)+v(411)*x(34)+v(414)*x(37)
v(475)=-(v(414)/v(343))
v(474)=-(v(413)/v(343))
v(473)=-(v(411)/v(343))
v(472)=-(v(418)/v(343))
v(471)=-(v(417)/v(343))
v(470)=-(v(412)/v(343))
v(469)=-(v(419)/v(343))
v(468)=-(v(416)/v(343))
v(467)=-(v(415)/v(343))
v(361)=v(422)*x(40)+v(420)*x(43)+v(423)*x(46)
v(484)=-(v(423)/v(361))
v(483)=-(v(422)/v(361))
v(482)=-(v(420)/v(361))
v(481)=-(v(427)/v(361))
v(480)=-(v(426)/v(361))
v(479)=-(v(421)/v(361))
v(478)=-(v(428)/v(361))
v(477)=-(v(425)/v(361))
v(476)=-(v(424)/v(361))
v(120)=v(121)*v(392)
v(122)=v(121)*v(385)
v(123)=v(121)*v(387)
v(124)=v(121)*v(389)
v(125)=v(121)*v(386)
v(126)=v(121)*v(391)
v(127)=v(121)*v(384)
v(128)=v(121)*v(388)
v(129)=v(121)*v(390)
v(130)=Fnew(1)*v(120)+Fnew(7)*v(126)+Fnew(4)*v(128)
v(131)=Fnew(2)*v(122)+Fnew(8)*v(124)+Fnew(5)*v(129)
v(132)=Fnew(3)*v(123)+Fnew(9)*v(125)+Fnew(6)*v(127)
v(133)=Fnew(4)*v(122)+Fnew(1)*v(124)+Fnew(7)*v(129)
v(134)=Fnew(5)*v(123)+Fnew(2)*v(125)+Fnew(8)*v(127)
v(135)=Fnew(6)*v(120)+Fnew(3)*v(126)+Fnew(9)*v(128)
v(136)=Fnew(7)*v(123)+Fnew(4)*v(125)+Fnew(1)*v(127)
v(137)=Fnew(8)*v(120)+Fnew(5)*v(126)+Fnew(2)*v(128)
v(138)=Fnew(9)*v(122)+Fnew(6)*v(124)+Fnew(3)*v(129)
v(139)=(v(130)*v(130))+(v(135)*v(135))+(v(137)*v(137))
v(151)=-v(139)/3d0
v(140)=(v(131)*v(131))+(v(133)*v(133))+(v(138)*v(138))
v(152)=-v(140)/3d0
v(141)=(v(132)*v(132))+(v(134)*v(134))+(v(136)*v(136))
v(147)=-v(141)/3d0
v(142)=v(130)*v(133)+v(131)*v(137)+v(135)*v(138)
v(143)=v(131)*v(134)+v(133)*v(136)+v(132)*v(138)
v(144)=v(132)*v(135)+v(130)*v(136)+v(134)*v(137)
v(145)=v(141)*(v(139)*v(140)-(v(142)*v(142)))-v(139)*(v(143)*v(143))+2d0*v(142)*v(143)*v(144)-v(140)*(v(144)*v(144))
v(149)=props(2)*(v(145)-sqrt(v(145)))
v(148)=1d0/v(145)**0.3333333333333333d0
v(437)=props(1)*v(148)
v(157)=1d0/v(306)**2
v(163)=1d0/v(157)**0.3333333333333333d0
v(429)=props(7)*v(163)
v(430)=v(157)*v(429)
v(158)=v(157)*((v(393)*v(393))+(v(395)*v(395))+(v(396)*v(396)))
v(162)=-v(158)/3d0
v(159)=v(157)*((v(394)*v(394))+(v(399)*v(399))+(v(400)*v(400)))
v(165)=-v(159)/3d0
v(160)=v(157)*((v(397)*v(397))+(v(398)*v(398))+(v(401)*v(401)))
v(166)=-v(160)/3d0
v(161)=((2d0/3d0)*v(160)+v(162)+v(165))*v(429)
v(164)=((2d0/3d0)*v(159)+v(162)+v(166))*v(429)
v(167)=((2d0/3d0)*v(158)+v(165)+v(166))*v(429)
v(175)=(v(394)*v(397)+v(398)*v(399)+v(400)*v(401))*v(430)
v(180)=(v(394)*v(396)+v(395)*v(399)+v(393)*v(400))*v(430)
v(182)=(v(396)*v(397)+v(395)*v(398)+v(393)*v(401))*v(430)
v(183)=1d0/v(324)**2
v(189)=1d0/v(183)**0.3333333333333333d0
v(431)=props(10)*v(189)
v(432)=v(183)*v(431)
v(184)=v(183)*((v(402)*v(402))+(v(404)*v(404))+(v(405)*v(405)))
v(188)=-v(184)/3d0
v(185)=v(183)*((v(403)*v(403))+(v(408)*v(408))+(v(409)*v(409)))
v(191)=-v(185)/3d0
v(186)=v(183)*((v(406)*v(406))+(v(407)*v(407))+(v(410)*v(410)))
v(192)=-v(186)/3d0
v(187)=((2d0/3d0)*v(186)+v(188)+v(191))*v(431)
v(190)=((2d0/3d0)*v(185)+v(188)+v(192))*v(431)
v(193)=((2d0/3d0)*v(184)+v(191)+v(192))*v(431)
v(201)=(v(403)*v(406)+v(407)*v(408)+v(409)*v(410))*v(432)
v(206)=(v(403)*v(405)+v(404)*v(408)+v(402)*v(409))*v(432)
v(208)=(v(405)*v(406)+v(404)*v(407)+v(402)*v(410))*v(432)
v(209)=1d0/v(343)**2
v(215)=1d0/v(209)**0.3333333333333333d0
v(433)=props(13)*v(215)
v(434)=v(209)*v(433)
v(210)=v(209)*((v(411)*v(411))+(v(413)*v(413))+(v(414)*v(414)))
v(214)=-v(210)/3d0
v(211)=v(209)*((v(412)*v(412))+(v(417)*v(417))+(v(418)*v(418)))
v(217)=-v(211)/3d0
v(212)=v(209)*((v(415)*v(415))+(v(416)*v(416))+(v(419)*v(419)))
v(218)=-v(212)/3d0
v(213)=((2d0/3d0)*v(212)+v(214)+v(217))*v(433)
v(216)=((2d0/3d0)*v(211)+v(214)+v(218))*v(433)
v(219)=((2d0/3d0)*v(210)+v(217)+v(218))*v(433)
v(227)=(v(412)*v(415)+v(416)*v(417)+v(418)*v(419))*v(434)
v(232)=(v(412)*v(414)+v(413)*v(417)+v(411)*v(418))*v(434)
v(234)=(v(414)*v(415)+v(413)*v(416)+v(411)*v(419))*v(434)
v(235)=1d0/v(361)**2
v(241)=1d0/v(235)**0.3333333333333333d0
v(435)=props(16)*v(241)
v(436)=v(235)*v(435)
v(236)=v(235)*((v(420)*v(420))+(v(422)*v(422))+(v(423)*v(423)))
v(240)=-v(236)/3d0
v(237)=v(235)*((v(421)*v(421))+(v(426)*v(426))+(v(427)*v(427)))
v(243)=-v(237)/3d0
v(238)=v(235)*((v(424)*v(424))+(v(425)*v(425))+(v(428)*v(428)))
v(244)=-v(238)/3d0
v(239)=((2d0/3d0)*v(238)+v(240)+v(243))*v(435)
v(242)=((2d0/3d0)*v(237)+v(240)+v(244))*v(435)
v(245)=((2d0/3d0)*v(236)+v(243)+v(244))*v(435)
v(253)=(v(421)*v(424)+v(425)*v(426)+v(427)*v(428))*v(436)
v(258)=(v(421)*v(423)+v(422)*v(426)+v(420)*v(427))*v(436)
v(260)=(v(423)*v(424)+v(422)*v(425)+v(420)*v(428))*v(436)
v(261)=v(149)-v(161)-v(187)-v(213)-v(239)+((2d0/3d0)*v(139)+v(147)+v(152))*v(437)
v(262)=v(149)-v(164)-v(190)-v(216)-v(242)+((2d0/3d0)*v(140)+v(147)+v(151))*v(437)
v(263)=v(149)-v(167)-v(193)-v(219)-v(245)+((2d0/3d0)*v(141)+v(151)+v(152))*v(437)
v(264)=-v(175)-v(201)-v(227)-v(253)+v(142)*v(437)
v(265)=-v(180)-v(206)-v(232)-v(258)+v(143)*v(437)
v(266)=-v(182)-v(208)-v(234)-v(260)+v(144)*v(437)
IF(dabs(props(5)).lt.0.1d-11) THEN
 v(438)=v(439)
 v(268)=v(438)
ELSE
 v(440)=1d0/props(5)
 v(268)=v(440)*(1d0-dexp(-(props(5)*v(439))))
ENDIF
v(276)=1d0/(props(3)+v(268))
v(441)=0.15d1*v(276)
v(269)=-v(262)/3d0
v(270)=-v(263)/3d0
v(273)=(2d0/3d0)*v(261)+v(269)+v(270)
v(271)=-v(261)/3d0
v(278)=(2d0/3d0)*v(263)+v(269)+v(271)
v(275)=(2d0/3d0)*v(262)+v(270)+v(271)
v(274)=v(273)*v(441)
v(277)=v(275)*v(441)
v(279)=v(278)*v(441)
v(280)=v(264)*v(441)
v(443)=2d0*v(280)
v(302)=-(v(280)*v(383))
v(281)=v(265)*v(441)
v(444)=2d0*v(281)
v(304)=-(v(281)*v(383))
v(282)=v(266)*v(441)
v(445)=2d0*v(282)
v(366)=props(17)*v(380)*(v(239)*v(274)+v(242)*v(277)+v(245)*v(279)+v(253)*v(443)+v(258)*v(444)+v(260)*v(445))
v(442)=(-1d0)+v(366)
v(377)=-(v(383)*(v(258)*v(365)+v(281)*v(442)))
v(375)=-(v(383)*(v(253)*v(365)+v(280)*v(442)))
v(373)=-(v(383)*(v(260)*v(365)+v(282)*v(442)))
v(348)=props(14)*v(380)*(v(213)*v(274)+v(216)*v(277)+v(219)*v(279)+v(227)*v(443)+v(232)*v(444)+v(234)*v(445))
v(446)=(-1d0)+v(348)
v(359)=-(v(383)*(v(232)*v(347)+v(281)*v(446)))
v(357)=-(v(383)*(v(227)*v(347)+v(280)*v(446)))
v(355)=-(v(383)*(v(234)*v(347)+v(282)*v(446)))
v(330)=props(11)*v(380)*(v(187)*v(274)+v(190)*v(277)+v(193)*v(279)+v(201)*v(443)+v(206)*v(444)+v(208)*v(445))
v(447)=(-1d0)+v(330)
v(341)=-(v(383)*(v(206)*v(329)+v(281)*v(447)))
v(339)=-(v(383)*(v(201)*v(329)+v(280)*v(447)))
v(337)=-(v(383)*(v(208)*v(329)+v(282)*v(447)))
v(311)=props(8)*v(380)*(v(161)*v(274)+v(164)*v(277)+v(167)*v(279)+v(175)*v(443)+v(180)*v(444)+v(182)*v(445))
v(448)=(-1d0)+v(311)
v(322)=-(v(383)*(v(180)*v(381)+v(281)*v(448)))
v(320)=-(v(383)*(v(175)*v(381)+v(280)*v(448)))
v(318)=-(v(383)*(v(182)*v(381)+v(282)*v(448)))
v(300)=-(v(282)*v(383))
R(1)=-props(3)-v(268)+sqrt(0.15d1*(2d0*v(264)**2+2d0*v(265)**2+2d0*v(266)**2+v(273)**2+v(275)**2+v(278)**2))
R(2)=1d0-v(274)*v(383)+v(121)*(-(statev(5)*v(388))-statev(8)*v(391)-statev(2)*v(392))
R(3)=1d0-v(277)*v(383)+v(121)*(-(statev(3)*v(385))-statev(9)*v(389)-statev(6)*v(390))
R(4)=1d0-v(279)*v(383)+v(121)*(-(statev(7)*v(384))-statev(10)*v(386)-statev(4)*v(387))
R(5)=v(302)+v(121)*(-(statev(5)*v(385))-statev(2)*v(389)-statev(8)*v(390))
R(6)=v(304)+v(121)*(-(statev(9)*v(384))-statev(3)*v(386)-statev(6)*v(387))
R(7)=v(300)+v(121)*(-(statev(10)*v(388))-statev(4)*v(391)-statev(7)*v(392))
R(8)=v(300)+v(121)*(-(statev(2)*v(384))-statev(5)*v(386)-statev(8)*v(387))
R(9)=v(302)+v(121)*(-(statev(3)*v(388))-statev(6)*v(391)-statev(9)*v(392))
R(10)=v(304)+v(121)*(-(statev(10)*v(385))-statev(7)*v(389)-statev(4)*v(390))
R(11)=1d0-v(383)*(v(161)*v(381)+v(274)*v(448))+statev(14)*v(449)+statev(17)*v(450)+statev(11)*v(451)
R(12)=1d0-v(383)*(v(164)*v(381)+v(277)*v(448))+statev(12)*v(452)+statev(15)*v(453)+statev(18)*v(454)
R(13)=1d0-v(383)*(v(167)*v(381)+v(279)*v(448))+statev(16)*v(455)+statev(13)*v(456)+statev(19)*v(457)
R(14)=v(320)+statev(14)*v(452)+statev(17)*v(453)+statev(11)*v(454)
R(15)=v(322)+statev(18)*v(455)+statev(15)*v(456)+statev(12)*v(457)
R(16)=v(318)+statev(19)*v(449)+statev(13)*v(450)+statev(16)*v(451)
R(17)=v(318)+statev(11)*v(455)+statev(17)*v(456)+statev(14)*v(457)
R(18)=v(320)+statev(12)*v(449)+statev(15)*v(450)+statev(18)*v(451)
R(19)=v(322)+statev(19)*v(452)+statev(13)*v(453)+statev(16)*v(454)
R(20)=1d0-v(383)*(v(187)*v(329)+v(274)*v(447))+statev(23)*v(458)+statev(26)*v(459)+statev(20)*v(460)
R(21)=1d0-v(383)*(v(190)*v(329)+v(277)*v(447))+statev(21)*v(461)+statev(24)*v(462)+statev(27)*v(463)
R(22)=1d0-v(383)*(v(193)*v(329)+v(279)*v(447))+statev(25)*v(464)+statev(22)*v(465)+statev(28)*v(466)
R(23)=v(339)+statev(23)*v(461)+statev(26)*v(462)+statev(20)*v(463)
R(24)=v(341)+statev(27)*v(464)+statev(24)*v(465)+statev(21)*v(466)
R(25)=v(337)+statev(28)*v(458)+statev(22)*v(459)+statev(25)*v(460)
R(26)=v(337)+statev(20)*v(464)+statev(26)*v(465)+statev(23)*v(466)
R(27)=v(339)+statev(21)*v(458)+statev(24)*v(459)+statev(27)*v(460)
R(28)=v(341)+statev(28)*v(461)+statev(22)*v(462)+statev(25)*v(463)
R(29)=1d0-v(383)*(v(213)*v(347)+v(274)*v(446))+statev(32)*v(467)+statev(35)*v(468)+statev(29)*v(469)
R(30)=1d0-v(383)*(v(216)*v(347)+v(277)*v(446))+statev(30)*v(470)+statev(33)*v(471)+statev(36)*v(472)
R(31)=1d0-v(383)*(v(219)*v(347)+v(279)*v(446))+statev(34)*v(473)+statev(31)*v(474)+statev(37)*v(475)
R(32)=v(357)+statev(32)*v(470)+statev(35)*v(471)+statev(29)*v(472)
R(33)=v(359)+statev(36)*v(473)+statev(33)*v(474)+statev(30)*v(475)
R(34)=v(355)+statev(37)*v(467)+statev(31)*v(468)+statev(34)*v(469)
R(35)=v(355)+statev(29)*v(473)+statev(35)*v(474)+statev(32)*v(475)
R(36)=v(357)+statev(30)*v(467)+statev(33)*v(468)+statev(36)*v(469)
R(37)=v(359)+statev(37)*v(470)+statev(31)*v(471)+statev(34)*v(472)
R(38)=1d0-v(383)*(v(239)*v(365)+v(274)*v(442))+statev(41)*v(476)+statev(44)*v(477)+statev(38)*v(478)
R(39)=1d0-v(383)*(v(242)*v(365)+v(277)*v(442))+statev(39)*v(479)+statev(42)*v(480)+statev(45)*v(481)
R(40)=1d0-v(383)*(v(245)*v(365)+v(279)*v(442))+statev(43)*v(482)+statev(40)*v(483)+statev(46)*v(484)
R(41)=v(375)+statev(41)*v(479)+statev(44)*v(480)+statev(38)*v(481)
R(42)=v(377)+statev(45)*v(482)+statev(42)*v(483)+statev(39)*v(484)
R(43)=v(373)+statev(46)*v(476)+statev(40)*v(477)+statev(43)*v(478)
R(44)=v(373)+statev(38)*v(482)+statev(44)*v(483)+statev(41)*v(484)
R(45)=v(375)+statev(39)*v(476)+statev(42)*v(477)+statev(45)*v(478)
R(46)=v(377)+statev(46)*v(479)+statev(40)*v(480)+statev(43)*v(481)
END

max_stars_repo_name: kengwit/MaterialModels

max_issues_repo_name: kengwit/MaterialModels

max_forks_repo_name: kengwit/MaterialModels


module type_GaussianGrid

    use, intrinsic :: ISO_Fortran_env, only: &
        stderr => ERROR_UNIT

    use spherepack_precision, only: &
        wp, & ! working precision
        ip ! integer precision

    use type_SphericalGrid, only: &
        SphericalGrid

    use gaussian_latitudes_and_weights_routines, only: &
        compute_gaussian_latitudes_and_weights

    ! Explicit typing only
    implicit none

    ! Everything is private unless stated otherwise
    private
    public :: GaussianGrid

    type, public, extends(SphericalGrid) :: GaussianGrid
        ! Type components
        real(wp), allocatable, public :: gaussian_weights(:)
    contains
        ! Type-bound procedures
        procedure, public  :: assert_initialized => assert_init_gaussian_grid
        procedure, public  :: create => create_gaussian_grid
        procedure, public  :: destroy => destroy_gaussian_grid
        procedure, public, nopass :: get_latitudes_and_gaussian_weights
        procedure, public  :: unformatted_print
        generic,   public  :: assignment (=) => copy_gaussian_grid
        procedure, private :: copy_gaussian_grid
    end type GaussianGrid

    ! Declare user-defined constructor
    interface GaussianGrid
        module procedure gaussian_grid_constructor
    end interface

contains

    function gaussian_grid_constructor(nlat, nlon) &
        result (return_value)

        ! Dummy arguments
        integer(ip), intent(in) :: nlat ! number of latitudinal points 0 <= theta <= pi
        integer(ip), intent(in) :: nlon ! number of longitudinal points 0 <= phi <= 2*pi
        type(GaussianGrid)      :: return_value

        call return_value%create(nlat, nlon)

    end function gaussian_grid_constructor

    subroutine assert_init_gaussian_grid(self, routine)

        ! Dummy arguments
        class(GaussianGrid), intent(inout) :: self
        character(len=*),    intent(in)    :: routine

        ! Check if object is usable
        select type(self)
            class is (GaussianGrid)
            if (.not.self%initialized) then
                write (stderr, '(a)') &
                    'Uninitialized object of class(GaussianGrid) in '//routine
            end if
        end select

    end subroutine assert_init_gaussian_grid

    subroutine copy_gaussian_grid(self, other)

        ! Dummy arguments
        class(GaussianGrid), intent(out) :: self
        class(GaussianGrid), intent(in)  :: other

        ! Check if object is usable
        if (.not.other%initialized) then
            error stop 'Uninitialized object of class(GaussianGrid): '&
                //'in assignment (=) '
        end if

        !  Make copies
        self%initialized = other%initialized
        self%NUMBER_OF_LONGITUDES = other%NUMBER_OF_LONGITUDES
        self%NUMBER_OF_LATITUDES = other%NUMBER_OF_LATITUDES
        self%LONGITUDINAL_MESH = other%LONGITUDINAL_MESH
        self%latitudes = other%latitudes
        self%longitudes = other%longitudes
        self%gaussian_weights = other%gaussian_weights
        self%grid_type = other%grid_type

    end subroutine copy_gaussian_grid

    subroutine create_gaussian_grid(self, nlat, nlon)

        ! Dummy arguments
        class(GaussianGrid), intent(inout) :: self
        integer(ip),         intent(in)    :: nlat ! number of latitudinal points 0 <= theta <= pi
        integer(ip),         intent(in)    :: nlon ! number of longitudinal points 0 <= phi <= 2*pi

        ! Ensure that object is usable
        call self%destroy()

        ! Set contants
        self%NUMBER_OF_LATITUDES = nlat
        self%NUMBER_OF_LONGITUDES = nlon

        ! Set the gaussian grid type
        allocate (self%grid_type, source='gaussian')

        ! Set longitudinal grid: 0 <= phi <= 2*pi
        call self%get_equally_spaced_longitudes(nlon, self%longitudes)

        ! Set latitudinal grid: 0 <= theta <= pi
        call self%get_latitudes_and_gaussian_weights(nlat, self%latitudes, self%gaussian_weights)

        ! Set initialization flag
        self%initialized = .true.

    end subroutine create_gaussian_grid

    subroutine destroy_gaussian_grid(self)

        ! Dummy arguments
        class(GaussianGrid), intent(inout) :: self

        ! Check initialization flag
        if (.not.self%initialized) return

        !  Release memory
        if (allocated(self%gaussian_weights)) deallocate (self%gaussian_weights)

        call self%destroy_grid()

        ! Reset initialization flag
        self%initialized = .false.

    end subroutine destroy_gaussian_grid

    subroutine get_latitudes_and_gaussian_weights(nlat, gaussian_latitudes, gaussian_weights)

        ! Dummy arguments
        integer(ip),           intent(in)    :: nlat ! number of latitudinal points
        real(wp), allocatable, intent(out)   :: gaussian_latitudes(:) ! latitudinal points: 0 <= theta <= pi
        real(wp), allocatable, intent(out)   :: gaussian_weights(:)

        ! Local variables
        integer(ip)  :: error_flag

        ! Check input argument
        if (nlat <= 0) then
            error stop 'Object of class(GaussianGrid): '&
                //'invalid argument nlat <= 0 '&
                //'in get_equally_spaced_latitudes'
        end if

        ! Allocate memory
        allocate (gaussian_latitudes(nlat))
        allocate (gaussian_weights(nlat))

        ! Compute gaussian weights and latitudes
        call compute_gaussian_latitudes_and_weights( &
            nlat, gaussian_latitudes, gaussian_weights, error_flag)

        ! Address error flag
        select case (error_flag)
            case(0)
                return
            case(1)
                error stop 'Object of class(GaussianGrid) '&
                    //'fails to satisfy nlat >= 0 '&
                    //'in get_latitudes_and_gaussian_weights'
            case default
                error stop 'Object of class(GaussianGrid): '&
                    //'undetermined error in '&
                    //'get_latitudes_and_gaussian_weights'
        end select

    end subroutine get_latitudes_and_gaussian_weights

    subroutine unformatted_print(self, header)

        ! Dummy arguments
        class(GaussianGrid), intent(inout) :: self
        character(len=*),    intent(in)    :: header

        ! Local variables
        integer(ip)  :: file_unit

        ! Check if object is usable
        call self%assert_initialized('unformatted_print')

        ! Write latitudes and longitudes
        call self%print_to_unformatted_binary_files(header)

        ! Write gaussian weights
        associate (wts => self%gaussian_weights)
            open( newunit=file_unit, &
                file=header//'gaussian_weights.dat', &
                status='replace', form='unformatted', &
                action='write', access='stream')
            write (file_unit) wts
            close( file_unit)
        end associate

    end subroutine unformatted_print

end module type_GaussianGrid


text:
While Davis itself may be flat, it is quite close to some killer areas for backpacking and other outdoor activities. Newcomers to the area should be aware that poison oak is common in areas below 5,000 feet elevation and hikers can easily contract contact dermatitis even when the shrubs are leafless. Also see Camping and Climbing.

 Trails and Hiking Areas

 Hiking Trails in Yolo County
  Blue Ridge/Rocky Ridge Blue Ridge can be accessed only via trailheads at the north end of the ridge.
  Cache Creek Nature Preserve
  http://www.yolocounty.org/Index.aspx?page379 Cache Creek Regional Park
  Davis Wetlands  At the site of the Waste Water Treatment facility, the wetlands provide a scenic, if somewhat flat, opportunity to view local Davis Wildlife wildlife.
  http://www.ecoanimal.com/dogfun/ca/north.html#yol Dog Park Directory  Yolo County
  Fremont Weir State Wildlife Area
  Frog Pond Trail
  Yolo County Grasslands Regional Park Grasslands Regional Park
  http://www.summitpost.org/mountain/rock/154977/littlebluepeak.html Little Blue Peak. Highest point in Yolo County (elevation 3120 ft).
  Pierce Canyon Falls Trail
  Putah Creek Riparian Reserve
  Putah Creek South Fork Preserve
  Sacramento Bypass Wildlife Area
  Vic Fazio Yolo Wildlife Area  Yolo Bypass Wildlife Area
  http://www.summitpost.org/object_list.php?parent_id171146&object_type1 SummitPost, Yolo County. Trails to the peaks: Fiske, Glascock Mountain, Little Blue Peak, Lowrey.
  http://www.yolohiker.org/sitemap/index.html YoloHiker  Hiking Trails in Yolo County. Yolo Bypass Wildlife Area, Pierce Canyon Falls Trail, Frog Pond Trail, Blue Ridge Trail, Rayhouse Road, Fiske Creek Trail, Putah Creek Reserve, Fremont Weir, Valley Vista Regional Park Trails.
  Urban hiking in City of Davis: The Greenbelt, Arboretum, Fitness Trails, Davis Bike Loop, Covell Drainage Channel Trail, Wildhorse Agricultural Buffer, Willowbank Ditch, Julie Partansky Pond, West Pond.

Within two hours one can be at many excellent backpacking and hiking areas. If youre willing to drive a little more the wilderness opportunities are endless. Its a little impossible to list all of the local destinations but here are a couple of ideas. See the Adventures Outside of Davis, Camping, or Climbing page for more places.  Dog owners should visit our Dogs page for petfriendly hikes.

 Hiking Trails in Colusa County
  Bear Valley (Colusa County) Bear Valley has one of the most spectacular displays of spring wildflowers in Northern California.
  Cache Creek Ridge Trail
  http://www.dfg.ca.gov/lands/wa/region2/cachecreek.html Cache Creek Wildlife Area
  http://www.dfg.ca.gov/lands/wa/region2/colusabypass.html Colusa Bypass Wildlife Area
  Colusa National Wildlife Refuge
  ColusaSacramento River State Recreation Area
  Colusa Weir Recreation Area
  Delevan National Wildlife Refuge
  Sacramento National Wildlife Refuge
  wiki:napa:Snow_Mountain_Wilderness Snow Mountain Wilderness Snow Mountain has the highest point of elevation in Colusa County (7,056 ft).
  http://www.yolohiker.org/sitemap/index.html Hiking YoloHiker  Trails in Colusa County. Billys Hill, Brophy Canyon, Cache Creek Ridge Trail, Cowboy Camp, Craig CanyonThompson Canyon Trail, High Bridge Trail, Lynch Canyon Trail, Signal/Eagle Rocks Loop on Walker Ridge.

 Hiking Trails in Other Counties
  Hiking Trails in California Statewide lists of trails, parks and hiking areas organized by county.
  wiki:sacramento:Hiking Sacramento County Hiking Trails  Sacramento Wiki. Includes: American River Bike Trail & more.
  wiki:vacaville:Hiking Solano County Hiking Trails  Vacaville Wiki. Includes: Stebbins Cold Canyon Reserve, Lake Solano County Park, Jepson Prairie Reserve & more.
  wiki:napa:Hiking Napa County & Lake County Hiking Trails  Napa Wiki. Includes: Long Canyon, Lake Berryessa trails & more.
  wiki:lincoln:Hiking Placer County Hiking Trails  Lincoln Wiki.
  wiki:sonomacounty:Hiking Sonoma County Hiking Trails  Sonoma County Wiki Includes: Mount Saint Helena & more.
  wiki:yubasutter:Hiking Yuba, Sutter, Sierra & Nevada County Hiking Trails  YubaSutter Wiki
  http://brtinsights.blogspot.com/2010/02/californiahikingtrailsdirectory.html California Hiking Trails Directory National, State & County Parks, Hiking Trails List, Hiking Groups, & Photos.
  wiki:PCT:Hiking Pacific Crest Trail Wiki describes the Pacific Crest Trail, which stretches from Mexico to Canada.

 Resources
Listed here are some great free resources for hiking and backpacking
  http://brtinsights.blogspot.com/2010/02/californiahikingtrailsdirectory.html California Hiking Trails Directory
  http://www.ndparks.com/camping/topcampingtoolspeopleforgetbringtrip.php Camping Tools and Equipment that People Forget to Bring on Their First Trip
  http://www.medicalofficer.net/wildernesssurvival/aboutwildernessfirstaid.php Basics of Wilderness First Aid
  http://www.medicalofficer.net/campingsafety/survivaltipsforcampinginthewilderness.php Survival Tips for Camping in the Wilderness

 Maps
One really doesnt need to buy maps with these available, however, common sense dictates that you should bring a printed copy of the map of your hiking destination with you on the hike.
  http://brtinsights.blogspot.com/2013/04/hikinginyolocountycamapto.html Yolo County Trailheads Map, http://brtinsights.blogspot.com/2013/04/hikingincolusacountycaliforniamap.html Colusa County Trailheads Map
  http://www.mapsofworld.com/usa/countymaps/california/yolocountymap.html Yolo County Map http://www.mapsofworld.com/usa/states/california/californiacountymap.html California County Map
  http://brtinsights.blogspot.com/2008/01/usingtopographicmapsfromacmemapper.html Acme Mapper provides dynamic online road, topographic and satellite maps. It shows coordinates for any point on the map. Enter coordinates or a landmark name and jump to that location on the map.
  Outdoor Adventures has a massive collection of maps, ranging from common topo maps, to guide books to a cool wallsized popup map.
  http://casil.ucdavis.edu/mapsurfer/ Mapsurfer   Hosted by our very own UC Davis, this map tool is truly awesome for the hiker.
  http://www.topozone.com TopoZone  Lets you search by location names which is very helpful.
  Shields Library Shields Library Map Room  Located in the basement is a complete set of California topographic maps to look at and photocopy. You dont even have to put them away once you are done looking at them!
  http://earth.google.com/ Google Earth  Satellite imagery of the entire planet.

 Libraries & Websites
  Outdoor Adventures has an entire bookshelf devoted to local, regional and national recreation areas.
  Shields Library  sections F859.3, F868 and GV199 contain an amazing array of hiking guides to the Sierra, the coast and beyond
  Davis Branch Library  has a large selection of recent trail guides
  http://brtinsights.blogspot.com/2010/02/californiahikingtrailsdirectory.html California Hiking Trails Directory
  http://kevingong.com/Hiking/ Kevin Gongs Hiking Page  trip reports of many local backpacking destinations
  http://www.waterfallswest.com/index.shtml Waterfalls West  good hiking directions to local (and not local) waterfalls
  http://www.sierranevadawild.gov Sierra Nevada Wild  current information and recent news about the 28 designated Wilderness areas in the Sierra Nevada
  http://www.highsierratopix.com HighSierraTopix.com  excellent internet community for everything Sierra Nevada
  http://mgordonphotography.com/sierra/sierra_new.htm Gordies Ramblings in the High Sierra  contains a good set of links to follow for more information

 Places to rent gear
  Outdoor Adventures

 Places to buy gear
  Kens Bike & Ski
  Outdoor Davis
  Outdoor Adventures has a biannual Equipment Swap in October and April. OA sells their used gear, vendors sell sample items, and the general public is allowed to sell or trade their goodquality used equipment and clothing.
  Big 5 Sporting Goods
  REI in Sacramento, Roseville, Folsom or Berkeley

 Clubs and Organizations
  Capay Valley Hiking Club (aka http://yolohiker.org/ Yolo Hiker)
  DOERS  Devoted Outdoor ExploreRS, sometimes also Davis Outdoor ExploreRS
  EEB Club
  Outdoor Adventures
     Users/MarieHuynh Marie Huynh and Users/RevChad Chad Van Schoelandt went on the weekend Outdoor Adventures OA backpacking trip to the Sierra over the summer of 2004. They thought it was awesome. http://pg.photos.yahoo.com/ph/dreameroflight/album?.dir3bfa&.srcph&store&prodid&.donehttp%3a//pg.photos.yahoo.com/ph//my_photos Photos
     Outdoor Adventures also offers a Wilderness Navigation course, which can be very helpful for the navigationally impaired hiker
  http://groups.yahoo.com/group/Sacramento_Hiking_Group/ Sacramento Hiking Group
  http://hiking.meetup.com/20/ Sacramento Hiking Meetup Group
  Sierra Club has many organized outdoor activities.

20080413 13:59:14 nbsp Can anyone recommend a good backpacking location near Davis (an hour or two away) that has some sort of water nearby to frolic in? Users/Ash
 You might want to check out the Sierra Club http://motherlode.sierraclub.org/deltasierra/GSPAGES/gspage0.htm scrambles. They are a combination hike and air mattress float trip.  My wife has done a couple of them and had a great time users/jimstewart

20090625 01:00:07 nbsp Where was the picture at the top of the page taken? Users/IDoNotExist

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      DOUBLE PRECISION FUNCTION MYMAX(A, B)
      DOUBLE PRECISION A, B
      MYMAX = DMAX1(A,B)
      RETURN
      END

max_stars_repo_name: blackrez/gonum

max_issues_repo_name: mingrammer/gonum

max_forks_repo_name: mingrammer/gonum


text:
C++*******************************************************************
C
C MEDIAN1
C
C To compute the median filtered version of the input image
C Uses NAG routine M01ANF
C
C Input parameters:
C  INPUT: input image
C  OUTPUT: output image
C  IBOX: Size of the box for which the median is comuted,
C        This should be an odd, positive integer greater than 3
C
C Syntax :   MEDIAN1 IN_FMT,OUT_FMT INPUT OUTPUT IBOX
C Example:   MEDIAN1 5,5 TEST TEST_MEDIAN 5
C
C JLP
C  Version of 25/03/90
C--**************************************************************************
	PROGRAM MEDIAN1
	INTEGER*4 PNTR_IN,PNTR_OUT
	CHARACTER NAMEIN*40,NAMEOUT*40,COMMENTS*80
	INTEGER*4 MADRID(1)
	COMMON /VMR/MADRID

10	FORMAT(A)

C Possibility of direct input of the parameters in the command line:
	CALL JLP_BEGIN

C Inquire the format of the files :
	CALL JLP_INQUIFMT

C Read the input image :
	NAMEIN=' '
	CALL JLP_VM_READIMAG(PNTR_IN,NX,NY,NAMEIN,COMMENTS)

C Obtain a name for the output :
	PRINT *,' OUTPUT IMAGE ?'
	READ(5,10) NAMEOUT

C Get dynamical memory for the output:
	I=NX*NY*4
	CALL JLP_GETVM(PNTR_OUT,I)
C
C NOW PICK UP THE BOX SIZE
C
	PRINT *,' Box size? (odd number > 3)'
	READ(5,*) IBOX
	IBOX=2*(IBOX/2)+1
	IF(IBOX.LT.3) THEN
	 PRINT *,' Box too small, I take Box size = 3'
	 IBOX=3
	ENDIF
	NB=NX/IBOX
	MB=NY/IBOX
	IF (NB.LE.3.OR.MB.LE.3) THEN
	 PRINT *,' FATAL ERROR : Box too large'
	 GOTO 999
	ENDIF
C
C HAVING GOT THE PARAMETERS , DO THE WORK
C
	CALL MEDIAN(MADRID(PNTR_IN),MADRID(PNTR_OUT),NX,NY,NX,IBOX)

C Ouput the results :
	LEN=INDEX(NAMEIN,'  ')-1
	WRITE(COMMENTS,11) NAMEIN(1:LEN),IBOX
11	FORMAT(' Median of: ',A,' Box:',I5)
	CALL JLP_WRITEIMAG(MADRID(PNTR_OUT),NX,NY,NX,NAMEOUT,COMMENTS)

999	CALL JLP_END
	STOP
	END
C**************************************************************************
      SUBROUTINE MEDIAN(IN,OUT,N,M,IDIM,IBOX)
C
C      IN IS THE INPUT ARRAY ; OUT THE OUTPUT ARRAY BOTH OF SIZE N BY M
C
C      THE ROUTINE APPLIES A MEDIAN FILTER TO THE INPUT ARRAY , YIELDING
C      THE OUTPUT ARRAY.
C      IBOX DEFINES THE SIZE OF THE BOX OVER WHICH THE MEDIAN IS COMPUTED
C
	REAL*4 IN(IDIM,*),OUT(IDIM,*),MEDN
C
C      SORT IS USED TO SORT THE PIXELS IN EACH BOX FOR THE MEDIAN
C
      DOUBLE PRECISION SORT(2500)
C
C      BOXSQ IS THE NUMBER OF PIXELS IN THE BOX (REAL) IBOXSC IS THE
C      INTEGER VERSION.
C
      IBOXSQ=IBOX*IBOX
      BOXSQ=IBOXSQ
C
C      FIRST ZEROIZE THE OUTPUT ARRAY
C
      DO 20 J=1,M
         DO 10 I=1,N
           OUT(I,J)=0
   10    CONTINUE
   20 CONTINUE
C
C      NOW FIGURE WHICH PIXELS CAN HAVE A BOX OF THIS SIZE ROUND THEM
C
      ISIDE=(IBOX-1)/2
      IS=1+ISIDE
      IE=N-ISIDE
      JS=1+ISIDE
      JE=M-ISIDE
      LOOP=0
C
C      NOW LOOP THROUGH ALL POSSIBLE PIXELS
C
      DO 500 J=JS,JE
         JHS=J-ISIDE
         JHE=J+ISIDE
         DO 400 I=IS,IE
            IHS=I-ISIDE
            IHE=I+ISIDE
C
C      SET THE PIXEL VALUES INTO A 1D ARRAY FOR SORTING
C
            K=1
            DO 200 JH=JHS,JHE
               DO 180 IH=IHS,IHE
                  SORT(K)=IN(IH,JH)
                  K=K+1
  180          CONTINUE
  200       CONTINUE
C
C Now call NAG M01ANF routine to do the sort
C
            IFAIL=0
            CALL M01ANF(SORT,1,IBOXSQ,IFAIL)
C
C and test that it is ok
C
            IF (IFAIL.NE.0) GO TO 400
C
C Now pick central value
C
            MEDN=SORT((IBOXSQ+1)/2)
C
C and store it in the output array
C
            OUT(I,J)=MEDN
  400    CONTINUE
  500 CONTINUE
C
C and return to calling segment
C
      RETURN
      END

ext: for


SUBROUTINE LA_TEST_ZGBSVX( FACT, TRANS, N, KL, KU, NRHS, AB, LDAB, AFB,&
     &  LDAFB, IPIV, EQUED, R, C, B, LDB, X, LDX, RCOND, FERR, BERR, &
     &  WORK, RWORK, INFO )
!
!  -- LAPACK95 interface driver routine (version 0.0) --
!     UNI-C, Denmark; Univ. of Tennessee, USA; NAG Ltd., UK
!     September 29, 1999
!
!  .. Use Statements ..
      USE LA_PRECISION, ONLY: WP => DP
      USE LA_AUXMOD, ONLY: ERINFO
      USE F95_LAPACK, ONLY: LA_GBSVX
!  .. Implicit Statement ..
      IMPLICIT NONE
!  .. Scalar Arguments ..
      CHARACTER(LEN=1), INTENT(INOUT) :: EQUED
      CHARACTER(LEN=1), INTENT(IN) :: FACT, TRANS
      INTEGER, INTENT(IN) :: N, KL, KU, NRHS, LDAB, LDB, LDAFB, LDX
      INTEGER, INTENT(INOUT) :: INFO
      REAL(WP), INTENT(OUT) :: RCOND
!  .. Array Arguments ..
      INTEGER, INTENT(INOUT) :: IPIV(1:N)
      REAL(WP), INTENT(OUT) :: RWORK(1:N)
      COMPLEX(WP), INTENT(INOUT) :: AB(1:LDAB,1:N), AFB(1:LDAFB,1:N),&
     &  B(1:LDB,1:N)
      REAL(WP), INTENT(OUT) :: BERR(1:NRHS), FERR(1:NRHS)
      COMPLEX(WP), INTENT(OUT) :: WORK(*), X(1:LDX,1:NRHS)
      REAL(WP), INTENT(INOUT) :: R(1:N), C(1:N)
!  .. Parameters ..
      CHARACTER(LEN=8),  PARAMETER :: SRNAME = 'LA_GBSVX'
      CHARACTER(LEN=14), PARAMETER :: SRNAMT = 'LA_TEST_ZGBSVX'
!  .. Common blocks ..
      INTEGER :: INFOTC
      COMMON /LINFO95/ INFOTC
!  .. Local Scalars ..
      CHARACTER(LEN=1) :: IFACT, ITRANS, IEQUED
      INTEGER :: I, J, IAB1, IAB2, IB1, IB2, IX1, IX2, IAFB1, IAFB2, IIPIV, &
     &  IR, IC, IFERR, IBERR, IKL, IKU
      LOGICAL, SAVE :: CTEST = .TRUE., ETEST = .TRUE.
!  .. Executable Statements ..

      IAB1 = KL+KU+1; IAB2 = N
      IB1 = N; IB2 = NRHS; IX1 = N; IX2 = NRHS
      IAFB1 =  2*KL+KU+1; IAFB2 = N; IIPIV = N; IFACT = FACT; ITRANS = TRANS
      IEQUED = EQUED; IR = N; IC = N; IFERR = NRHS; IBERR = NRHS
      IKL = KL; IKU = KU
      I = INFO / 100; J = INFO - I*100
      SELECT CASE(I)
        CASE(0)
          IF (NRHS == 1) THEN
            CALL LA_GBSVX( AB(1:IAB1,1:IAB2), B(1:IB1,1), X(1:IX1,1), KL, &
            AFB(1:IAFB1,1:IAFB2), IPIV(1:IIPIV), IFACT, ITRANS, EQUED, &
            R(1:IR), C(1:IC), FERR(1), BERR(1), RCOND, RPVGRW=RWORK(1) )
            INFO = INFOTC
          ELSE
            CALL LA_GBSVX( AB(1:IAB1,1:IAB2), B(1:IB1,1:IB2), X(1:IX1,1:IX2),&
            KL, AFB(1:IAFB1,1:IAFB2), IPIV(1:IIPIV), IFACT, ITRANS, EQUED, &
            R(1:IR), C(1:IC), FERR(1:IFERR), BERR(1:IBERR), RCOND, &
            RWORK(1) )
            INFO = INFOTC
          ENDIF
        CASE (1)
          IAB2 = 0
        CASE (2)
          IB1 = IAB2-1
        CASE (3)
    I3:   SELECT CASE(J)
            CASE(1)
              IX1 = IAB2-1
            CASE (2)
              IX2 = IB2-1
            CASE (3:)
              CALL UESTOP(SRNAMT)
          END SELECT I3
        CASE (4)
    I4:   SELECT CASE(J)
            CASE (1)
              IKL = - 1
            CASE (2)
              IKL = - 1
            CASE (3:)
              CALL UESTOP(SRNAMT)
          END SELECT I4
        CASE (5)
    I5:   SELECT CASE(J)
            CASE(1)
              iafb1 = iab1-1
            CASE (2)
              IAFB2 = iab2-1
            CASE (3:)
              CALL UESTOP(SRNAMT)
          END SELECT I5
        CASE (6)
          IIPIV = IAB2-1
        CASE (7)
    I7:   SELECT CASE(J)
            CASE(1)
              IFACT = 'T'
            CASE (2)
              SELECT CASE(NRHS)
                CASE (2:)
                  CALL LA_GBSVX(AB(1:IAB1,1:IAB2), B(1:IB1,1:IB2), X(1:IX1,1:IX2), &
                  KL, AFB(1:IAFB1,1:IAFB2), FACT='F', INFO = INFO)
                CASE (1)
                  CALL LA_GBSVX(AB(1:IAB1,1:IAB2), B(1:IB1,1), X(1:IX1,1), &
                  KL, AFB(1:IAFB1,1:IAFB2), FACT='F', INFO=INFO)
                CASE DEFAULT
                  CALL UESTOP(SRNAMT)
              END SELECT
            CASE (3)
              SELECT CASE(NRHS)
                CASE (2:)
                  CALL LA_GBSVX(AB(1:IAB1,1:IAB2), B(1:IB1,1:IB2), X(1:IX1,1:IX2), &
                  KL, IPIV=IPIV(1:IIPIV), FACT='F', INFO=INFO)
                CASE (1)
                  CALL LA_GBSVX(AB(1:IAB1,1:IAB2), B(1:IB1,1), X(1:IX1,1), &
                  KL, IPIV=IPIV(1:IIPIV), FACT='F', INFO=INFO)
                CASE DEFAULT
                  CALL UESTOP(SRNAMT)
              END SELECT
            CASE (4)
              SELECT CASE(NRHS)
                CASE (2:)
                  CALL LA_GBSVX(AB(1:IAB1,1:IAB2), B(1:IB1,1:IB2), X(1:IX1,1:IX2), &
                  KL, FACT='F', INFO=INFO)
                CASE (1)
                  CALL LA_GBSVX(AB(1:IAB1,1:IAB2), B(1:IB1,1), X(1:IX1,1), &
                  KL, FACT='F', INFO=INFO)
                CASE DEFAULT
                  CALL UESTOP(SRNAMT)
              END SELECT
            CASE (5:)
              CALL UESTOP(SRNAMT)
          END SELECT I7

        CASE (8)
          ITRANS = 'x'; IFACT = 'N'; IEQUED = 'N'
        CASE (9)
          IFACT = 'F'; ITRANS = 'N'
        I9:   SELECT CASE (J)
                CASE (1)
                  IEQUED = 'T'
                CASE (2)
                  IEQUED = 'R'
                  SELECT CASE (NRHS)
                    CASE (2:)
                      CALL LA_GBSVX(AB(1:IAB1,1:IAB2), B(1:IB1,1:IB2), X(1:IX1,1:IX2), &
                      KL, AFB(1:IAFB1,1:IAFB2), IPIV(1:IIPIV), IFACT, ITRANS, IEQUED, &
                      INFO=INFO)
                    CASE (1)
                      CALL LA_GBSVX(AB(1:IAB1,1:IAB2), B(1:IB1,1), X(1:IX1,1), &
                      AFB=AFB(1:IAFB1,1:IAFB2), IPIV=IPIV(1:IIPIV), FACT=IFACT, &
&                       TRANS=ITRANS, EQUED=IEQUED, &
                      INFO=INFO)
                    CASE DEFAULT
                      CALL UESTOP(SRNAMT)
                  END SELECT
                CASE (3)
                  IEQUED = 'B'
                  SELECT CASE (NRHS)
                    CASE (2:)
                      CALL LA_GBSVX(AB(1:IAB1,1:IAB2), B(1:IB1,1:IB2), X(1:IX1,1:IX2), &
                      AFB=AFB(1:IAFB1,1:IAFB2), IPIV=IPIV(1:IIPIV), FACT=IFACT, &
&                       TRANS=ITRANS, EQUED=IEQUED, &
                      INFO=INFO)
                    CASE (1)
                      CALL LA_GBSVX(AB(1:IAB1,1:IAB2), B(1:IB1,1), X(1:IX1,1), &
                      AFB=AFB(1:IAFB1,1:IAFB2), IPIV=IPIV(1:IIPIV), FACT=IFACT, &
&                       TRANS=ITRANS, EQUED=IEQUED, &
                      INFO=INFO)
                    CASE DEFAULT
                      CALL UESTOP(SRNAMT)
                  END SELECT
                CASE (4)
                  IEQUED = 'C'
                  SELECT CASE (NRHS)
                    CASE (2:)
                      CALL LA_GBSVX(AB(1:IAB1,1:IAB2), B(1:IB1,1:IB2), X(1:IX1,1:IX2), &
                      AFB=AFB(1:IAFB1,1:IAFB2), IPIV=IPIV(1:IIPIV), FACT=IFACT, &
&                       TRANS=ITRANS, EQUED=IEQUED, &
                      INFO=INFO)
                    CASE (1)
                      CALL LA_GBSVX(AB(1:IAB1,1:IAB2), B(1:IB1,1), X(1:IX1,1), &
                      AFB=AFB(1:IAFB1,1:IAFB2), IPIV=IPIV(1:IIPIV), FACT=IFACT, &
&                       TRANS=ITRANS, EQUED=IEQUED, &
                      INFO=INFO)
                    CASE DEFAULT
                      CALL UESTOP(SRNAMT)
                  END SELECT
                CASE (5:)
                  CALL UESTOP(SRNAMT)
              END SELECT I9
            CASE (10)
              IFACT = 'F'; R(1) = -1.0_WP; C = +1.0_WP
              I10:   SELECT CASE (J)
                       CASE(1)
                         IR = IAB2-1
                       CASE(2)
                         IEQUED = 'R'
                       CASE(3)
                         IEQUED = 'B'
                       CASE (4:)
                         CALL UESTOP(SRNAMT)
                     END SELECT I10

            CASE (11)
              IFACT = 'F'; R = +1.0_WP; C(1) = -1.0_WP
              I11:  SELECT CASE (J)
                      CASE(1)
                        IC = IAB2-1
                      CASE(2)
                        IEQUED = 'C'
                      CASE(3)
                        IEQUED = 'B'
                      CASE (4:)
                        CALL UESTOP(SRNAMT)
                    END SELECT I11
            CASE (12)
            IFERR = IB2-1; R = +1.0_WP; C = +1.0_WP; IC = IAB2; IR = IAB2
            CASE (13)
            IBERR = IB2-1; R = +1.0_WP; C = +1.0_WP; IC = IAB2; IR = IAB2
            CASE (14:)
              CALL UESTOP(SRNAMT)
          END SELECT

          IF( .NOT. (( I == 7 .AND. ( J == 2 .OR. J == 3 .OR. J == 4 ) ) &
     &      .OR. ( I == 9 .AND. ( J == 2 .OR. J == 3 .OR. J == 4 ) )) &
     &      .AND. I/=0) THEN
            SELECT CASE (NRHS)
              CASE (2:)
                CALL LA_GBSVX(AB(1:IAB1,1:IAB2), B(1:IB1,1:IB2), X(1:IX1,1:IX2), &
                KL=IKL, AFB=AFB(1:IAFB1,1:IAFB2), IPIV=IPIV(1:IIPIV), FACT=IFACT, &
     &            TRANS=ITRANS, EQUED=IEQUED, &
     &            R=R(1:IR), C=C(1:IC), FERR=FERR(1:IFERR), BERR=BERR(1:IBERR), &
     &            RCOND=RCOND, &
                RPVGRW=RWORK(1), INFO=INFO)
              CASE (1)

                IF( INFO /= 302 .AND. I /= 12 .AND. I /= 13 )THEN
                  CALL LA_GBSVX(AB(1:IAB1,1:IAB2), B(1:IB1,1), X(1:IX1,1), &
                  KL=IKL,  AFB=AFB(1:IAFB1,1:IAFB2), IPIV=IPIV(1:IIPIV), FACT=IFACT, &
     &              TRANS=ITRANS, EQUED=IEQUED, &
     &              R=R(1:IR), C=C(1:IC), FERR=FERR(1), BERR=BERR(1), RCOND=RCOND, &
     &              RPVGRW=RWORK(1), INFO=INFO)
                ELSE
                  CALL ERINFO( -I, SRNAME, INFO )
                END IF
            END SELECT
          END IF
          EQUED = IEQUED

          CALL LA_AUX_AA01( I, CTEST, ETEST, SRNAMT )
        END SUBROUTINE LA_TEST_ZGBSVX


text:
Tomato Hornworms are the biggest caterpillars you ever did see. Theyll eat your plants and your parenting mother too. Planting basil near your tomatoes will help discourage hornworms from horning in on your tomatoes tomato crop. They grow up and become the 5Spotted Hawk Moth. To do more learnin about this bug, see wiki:WikiPedia:Tomato_hornworm wikipedias article.

post a picture! Users/MichelleAccurso

I knew these things really existed!  When I was a kid, my dad grew tomato plants and there were always lots of giant green caterpillars on em.  I havent seen one in years, and was starting to credit their existence to an overactive childhood imagination.  Yeahpost a photo, Users/ArlenAbraham Arlen! Users/SummerSong ss

Not the best picture  my camera kept wanting to focus on the background instead of the hornworm.  This guy was about 2 inches long and about as big around as 2 pencils side by side.  I squashed about twenty or thirty of these guys this year.  Ive never heard of the basil trick, but I had a volunteer basil come up this year right next to my tomatos, and the hornworms are gone.  Coincidence?  Seasonal?  Mystery solved?  Hornworms also come in a tobacco variety.  Users/AlexPomeranz

Growing up in North Davis my dad used to (and still does) grow Tomatoes in his seasonal veggie plots, as a kid I would find and remove these worms, they would get pretty damn big  Users/StevenDaubert


! { dg-do compile }
! { dg-options "-Wno-all" }
!
! By default, for free-form source code: Error out
! But due to the explicit -Wno-all, compile w/o warning
!
      print *, 1                                                                                + 2
      end


text:
subroutine a1(a,b,n)
    integer i,n
    real a(n),b(n)
!$omp parallel
!$omp do
  do i=2,n
    b(i)= a(i)
  enddo
!$omp end do nowait
!$omp end parallel
  print *,'end'
end subroutine a1

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


module elast

public :: strain2stress
public :: stress2traction
public :: traction_components
public :: max_shear_stress

public :: read_halfspace_file

!--------------------------------------------------------------------------------------------------!
contains
!--------------------------------------------------------------------------------------------------!

subroutine strain2stress(strain,lame_param,shear_modulus,stress)
!----
! Calculate the stress tensor from the strain tensor and elastic moduli, assuming a linear,
! isotropic medium.
!----

implicit none

! Arguments
double precision :: stress(3,3), strain(3,3), lame_param, shear_modulus

! Local variables
double precision :: diag

diag = strain(1,1) + strain(2,2) + strain(3,3)

! Diagonal terms
stress(1,1) = lame_param*diag + 2.0d0*shear_modulus*strain(1,1)
stress(2,2) = lame_param*diag + 2.0d0*shear_modulus*strain(2,2)
stress(3,3) = lame_param*diag + 2.0d0*shear_modulus*strain(3,3)

! Off-diagonal terms
stress(1,2) = 2.0d0*shear_modulus*strain(1,2)
stress(1,3) = 2.0d0*shear_modulus*strain(1,3)
stress(2,3) = 2.0d0*shear_modulus*strain(2,3)
stress(2,1) = stress(1,2)
stress(3,1) = stress(1,3)
stress(3,2) = stress(2,3)

return
end subroutine

!--------------------------------------------------------------------------------------------------!

subroutine stress2traction(stress_tensor,normal_vector,traction_vector)
!----
! Calculate the traction vector resolved onto a plane from the stress tensor and the normal vector
! to the plane.
!----

use io, only: stderr
implicit none

! Arguments
double precision :: stress_tensor(3,3), normal_vector(3), traction_vector(3)

! Local variables
integer :: i, j
double precision :: magnitude

! Initialize traction vector
traction_vector = 0.0d0

! Normalize the normal vector
magnitude = normal_vector(1)*normal_vector(1) + &
            normal_vector(2)*normal_vector(2) + &
            normal_vector(3)*normal_vector(3)

if (magnitude.lt.1.0d-8) then
    write(stderr,*) 'stress2traction: normal vector magnitude is zero, setting traction to zero'
    return
endif

normal_vector = normal_vector/magnitude

! traction_vector = stress_tensor*normal_vector
do i = 1,3
    do j = 1,3
        traction_vector(i) = traction_vector(i) + stress_tensor(i,j)*normal_vector(j)
    enddo
enddo

return
end subroutine

!--------------------------------------------------------------------------------------------------!

subroutine traction_components(trac_vec,nor_vec,trac_nor,trac_str,trac_dip)
!----
! Calculate the normal, strike-parallel, and up-dip components of the traction vector on a plane
! defined by its normal vector.
!----

use geom, only: normal2strike, normal2updip
use algebra, only: dot_product

implicit none

! Arguments
double precision :: trac_vec(3), nor_vec(3), trac_nor, trac_str, trac_dip

! Local variables
double precision :: str_vec(3), dip_vec(3)

! Calculate unit vectors parallel to strike and up-dip
call normal2strike(nor_vec,str_vec)
call normal2updip(nor_vec,dip_vec)

! Project the traction vector onto normal, along-strike, and horizontal up-dip directions
call dot_product(trac_vec,nor_vec,trac_nor)
call dot_product(trac_vec,str_vec,trac_str)
call dot_product(trac_vec,dip_vec,trac_dip)

return
end subroutine

!--------------------------------------------------------------------------------------------------!

subroutine max_shear_stress(stress,max_shear)
!----
! Calculate maximum shear stress from a stress tensor
!----

implicit none

! Arguments
double precision :: stress(3,3), max_shear

! Local variables
double precision :: s11_s22, s11_s33, s22_s33, s12s12, s13s13, s23s23

s11_s22 = stress(1,1)-stress(2,2)
s11_s33 = stress(1,1)-stress(3,3)
s22_s33 = stress(2,2)-stress(3,3)

s12s12 = stress(1,2)*stress(1,2)
s13s13 = stress(1,3)*stress(1,3)
s23s23 = stress(2,3)*stress(2,3)

max_shear = (1.0d0/6.0d0)*(s11_s22*s11_s22+s11_s33*s11_s33+s22_s33*s22_s33) + s12s12+s13s13+s23s23

! Make sure it has units of Pa
max_shear = dsqrt(max_shear)

return
end subroutine

!--------------------------------------------------------------------------------------------------!

subroutine read_halfspace_file(halfspace_file,poisson,shearmod,lame,ierr)
!----
! Read half-space elastic moduli from file (if provided), in format:
!     lbl1 value1 lbl2 value2 [lbl3 value3]
!
! For example, if you have the shear modulus and Young's modulus, the file would be:
!     shearmod 40e9 young 8e10
!----

use io, only: stderr, line_count, fileExists

implicit none

! Arguments
character(len=*) :: halfspace_file
double precision :: poisson, shearmod, lame

! Local variables
integer :: i, j, ios, ierr
character(len=512) :: input_line
character(len=32) :: lbl(4)
character(len=8) :: halfspace_file_mode
double precision :: val(3)

ierr = 0

! If the file is not defined, use defaults set in gcmdln
if (halfspace_file.eq.'') then
    return
endif

! Check whether the half-space file exists
if (.not.fileExists(halfspace_file)) then
    write(stderr,*) 'read_halfspace: no half-space file named '//trim(halfspace_file)//' found'
    write(stderr,*) 'using default elastic moduli values'
    return
endif

! Initialize modulus variables
lbl = ''
val = 0.0d0

! Open the file and read the half-space elastic moduli
halfspace_file_mode = 'modern'
open(unit=15,file=halfspace_file,status='old')
read(15,'(A)') input_line

! For back-compatibility, can read parameters in these formats (print warning):
!     vp vs dens
!     "lame" lame shear_mod
read(input_line,*,end=1501,err=1500,iostat=ios) lbl(1:4)
1501 if (ios.ne.0) then
    halfspace_file_mode = 'legacy'
    ! Found fewer than 4 inputs
    if (lbl(1).eq.'lame'.or.lbl(1).eq.'Lame') then
        read(lbl(2),*,err=1500,iostat=ios) val(1)
        read(lbl(3),*,err=1500,iostat=ios) val(2)
        lbl(1) = 'lame'
        lbl(2) = 'shearmod'
        lbl(3:4) = ''
        write(stderr,*) 'read_halfspace: WARNING: you are using a legacy half-space format ("lame" lame shear)'
        write(stderr,*) 'In the future, you should use labeled inputs, e.g.: lame 40e9 shearmod 40e9'
    else
        read(lbl(1),*,err=1500,iostat=ios) val(1)
        read(lbl(2),*,err=1500,iostat=ios) val(2)
        read(lbl(3),*,err=1500,iostat=ios) val(3)
        lbl(1) = 'vp'
        lbl(2) = 'vs'
        lbl(3) = 'dens'
        write(stderr,*) 'read_halfspace: WARNING: you are using a legacy half-space format (vp vs dens)'
        write(stderr,*) 'In the future, you should use labeled inputs, e.g.: vp 6800 vs 3926 dens 3000'
    endif
endif

! Read in updated format, which can take three (labeled) inputs for a minimum of 4 items on line
! Examples: shear 40e9 lame 35e9
!           vp 6500 vs 3700 dens 3000
if (halfspace_file_mode.eq.'modern') then
    read(input_line,*,iostat=ios) lbl(1), val(1), lbl(2), val(2), lbl(3), val(3)

    ! Then check for two inputs
    if (ios.ne.0) then
        read(input_line,*,iostat=ios) lbl(1), val(1), lbl(2), val(2)
    endif
endif

1500 if (ios.ne.0) then
    write(stderr,*) 'read_halfspace: error reading modulus labels and values'
    write(stderr,*) 'offending line: '//trim(input_line)
    ierr = 1
    return
endif

close(15)

! Check labels make sense
do i = 1,2
    if (lbl(i).eq.'mu' .or. lbl(i).eq.'shear_modulus' .or. lbl(i).eq.'shear' .or. &
                                                                     lbl(i).eq.'shearmod') then
        lbl(i) = 'shearmod'
    elseif (lbl(i).eq.'lame' .or. lbl(i).eq.'lambda') then
        lbl(i) = 'lame'
    elseif (lbl(i).eq.'poisson' .or. lbl(i).eq.'nu') then
        lbl(i) = 'poisson'
    elseif (lbl(i).eq.'young') then
        lbl(i) = 'young'
    elseif (lbl(i).eq.'vp'.or.lbl(i).eq.'vs'.or.lbl(i).eq.'dens') then
        ! Okay
    else
        write(stderr,*) 'read_halfspace: elastic modulus "',trim(lbl(i)),'" not implemented'
        write(stderr,*) 'specify two of: shear_modulus, lame, poisson, or young OR vp vs dens'
        ierr = 2
        return
    endif
enddo

! Calculate Poisson's ratio, Lame parameter, and shear modulus
do i = 1,2

    j = mod(i,2)+1

    if (lbl(i).eq.'shearmod'.and.lbl(j).eq.'lame') then
        ! Calculate Poisson's ratio from Lame's parameter and shear modulus
        shearmod = val(i)
        lame = val(j)
        poisson = lame/(2.0d0*(lame+shearmod))

    elseif (lbl(i).eq.'shearmod'.and.lbl(j).eq.'poisson') then
        ! Calculate Lame's parameter from Poisson's ratio and shear modulus
        shearmod = val(i)
        poisson = val(j)
        lame = 2.0d0*shearmod*poisson/(1.0d0-2.0d0*poisson)

    elseif (lbl(i).eq.'shearmod'.and.lbl(j).eq.'young') then
        ! Calculate Lame's parameter and Poisson's ratio from Young's modulus and shear modulus
        shearmod = val(i)
        lame = shearmod*(val(j)-2.0d0*shearmod)/(3.0d0*shearmod-val(j))
        poisson = lame/(2.0d0*(lame+shearmod))

    elseif (lbl(i).eq.'lame'.and.lbl(j).eq.'poisson') then
        ! Calculate shear modulus from Lame's parameter and Poisson's ratio
        lame = val(i)
        poisson = val(j)
        shearmod = lame*(1.0d0-2.0d0*poisson)/(2.0d0*poisson)

    elseif (lbl(i).eq.'lame'.and.lbl(j).eq.'young') then
        ! Calculate shear modulus and Poisson's ratio from Lame's parameter and Young's modulus
        lame = val(i)
        poisson = 2.0d0*lame/(val(j)+lame+sqrt(val(j)*val(j)+9.0d0*lame*lame+2.0d0*lame*val(j)))
        shearmod = lame*(1.0d0-2.0d0*poisson)/(2.0d0*poisson)

    elseif (lbl(i).eq.'poisson'.and.lbl(j).eq.'young') then
        ! Calculate shear modulus and Lame's parameter from Poisson's ratio and Young's modulus
        poisson = val(i)
        shearmod = val(j)/(2.0d0*(1.0d0+poisson))
        lame = 2.0d0*shearmod*poisson/(1.0d0-2.0d0*poisson)

    endif
enddo

if (lbl(1).eq.'vp'.and.lbl(2).eq.'vs'.and.lbl(3).eq.'dens') then
    ! Calculate shear modulus, Lame's parameter, and Poisson's ratio from seismic velocities and density
    shearmod = val(2)*val(2)*val(3)
    lame = val(1)*val(1)*val(3) - 2.0d0*shearmod
    poisson = lame/(2.0d0*(lame+shearmod))
endif

return
end subroutine

!--------------------------------------------------------------------------------------------------!

! !--------------------------------------------------------------------------------------------------!
!
!     subroutine calc_traction_components(traction,normal,strike,updip,traction_components)
!     implicit none
!     ! I/O variables
!     double precision :: traction(3), normal(3), strike(3), updip(3), traction_components(3)
!     ! Local variables
!     integer :: i
!
!     ! Project the traction vector on normal, along-strike, and horizontal up-dip directions
!     traction_components = 0.0d0
!     do i = 1,3
!         traction_components(1) = traction_components(1) + traction(i)*normal(i)
!         traction_components(2) = traction_components(2) + traction(i)*strike(i)
!         traction_components(3) = traction_components(3) + traction(i)*updip(i)
!     enddo
!     return
!     end subroutine calc_traction_components
!
!
! !--------------------------------------------------------------------------------------------------!
!
!     subroutine calc_plane_unit_vectors(strike,dip,unit_normal,unit_strike,unit_updip)
!     implicit none
!     ! I/O variables
!     double precision :: strike, dip, unit_normal(3), unit_strike(3), unit_updip(3)
!     ! Local variables
!     double precision, parameter :: pi=datan(1.0d0)*4.0d0, d2r=pi/180.0d0, pi2=pi/2.0d0
!     double precision :: str_rad, dip_rad
!
!     str_rad = strike*d2r
!     dip_rad = dip*d2r
!
!     unit_normal(1) = dsin(dip_rad)*dsin(str_rad+pi2)
!     unit_normal(2) = dsin(dip_rad)*dcos(str_rad+pi2)
!     unit_normal(3) = dcos(dip_rad)
!
!     unit_strike(1) = dsin(str_rad)
!     unit_strike(2) = dcos(str_rad)
!     unit_strike(3) = 0.0d0
!
!     unit_updip(1) = dcos(dip_rad)*dsin(str_rad-pi2)
!     unit_updip(2) = dcos(dip_rad)*dcos(str_rad-pi2)
!     unit_updip(3) = dsin(dip_rad)
!
!     return
!     end
!
! !--------------------------------------------------------------------------------------------------!
!
!     subroutine rotate_strain(strain,strike)
!     !----
!     ! Rotate strain matrix from (x=str, y=updip horizontal, z=up) to (x=E, y=N, z=up)
!     !----
!     implicit none
!     ! I/O variables
!     double precision, dimension(3,3) :: strain
!     double precision :: strike
!     ! Local variables
!     double precision, dimension(3,3) :: rot, rot_transpose, tmp
!     double precision, parameter :: pi=datan(1.0d0)*4.0d0, d2r=pi/180.0d0
!     rot(1,1) = dcos(d2r*strike-pi/2.0d0)
!     rot(2,2) = rot(1,1)
!     rot(1,2) = dsin(d2r*strike-pi/2.0d0)
!     rot(2,1) = -rot(1,2)
!     rot(1,3) = 0.0d0
!     rot(2,3) = 0.0d0
!     rot(3,1) = 0.0d0
!     rot(3,2) = 0.0d0
!     rot(3,3) = 1.0d0
!     call mattr(rot_transpose,rot)
!     call matmult(tmp,rot,strain)
!     call matmult(strain,tmp,rot_transpose)
!     return
!     end
!
! !--------------------------------------------------------------------------------------------------!
!
!     subroutine mattr(matout,matin)
!     !----
!     ! Transpose a 3x3 matrix
!     !----
!     implicit none
!     ! I/O variables
!     double precision, dimension(3,3) :: matout, matin
!     ! Local variables
!     integer :: i, j
!     do i = 1,3
!         do j = 1,3
!             matout(i,j) = matin(j,i)
!         enddo
!     enddo
!     return
!     end
!
! !--------------------------------------------------------------------------------------------------!
!
!     subroutine matmult(matout,mat1,mat2)
!     !----
!     ! Multiply two 3x3 matrices
!     !----
!     implicit none
!     ! I/O variables
!     double precision, dimension(3,3) :: mat1, mat2, matout
!     ! Local variables
!     integer :: i, j, k
!     do i = 1,3
!         do j = 1,3
!             matout(i,j) = 0.0d0
!             do k = 1,3
!                 matout(i,j) = matout(i,j) + mat1(i,k)*mat2(k,j)
!             enddo
!         enddo
!     enddo
!     return
!     end
! !--------------------------------------------------------------------------------------------------!
!
!
!
! subroutine calc_tractions(trac_vector,shear,shearmax,normal,coulomb,stress,trg)
! !----
! ! Calculate various tractions resolved onto a fault geometry from a stress tensor
! !----
!
! use trig, only: pi, d2r
! implicit none
!
! ! Arguments
! double precision :: trac_vector(3), shear, shearmax, normal, coulomb, stress(3,3), trg(4)
!
! ! Local variables
! integer :: i
! double precision :: str, dip, rak, fric
! double precision :: n(3), s(3), r(3)
!
! str = trg(1)*d2r
! dip = trg(2)*d2r
! rak = trg(3)*d2r
! fric = trg(4)
!
! ! Unit normal vector to target plane
! n(1) = dsin(dip)*dsin(str+pi/2.0d0)
! n(2) = dsin(dip)*dcos(str+pi/2.0d0)
! n(3) = dcos(dip)
!
! ! Traction vector is stress matrix times normal vector: t = S*n
! do i = 1,3
!     trac_vector(i) = stress(i,1)*n(1) + stress(i,2)*n(2) + stress(i,3)*n(3)
! enddo
!
! ! Normal component of traction is parallel to unit normal vector; take dot product
! normal = 0.0d0
! do i = 1,3
!     normal = normal + trac_vector(i)*n(i)
! enddo
!
! ! Shear component of traction is difference between total traction vector and normal traction
! s(1) = trac_vector(1) - normal*n(1)
! s(2) = trac_vector(2) - normal*n(2)
! s(3) = trac_vector(3) - normal*n(3)
! shearmax = dsqrt(s(1)*s(1) + s(2)*s(2) + s(3)*s(3))
!
! ! Compute unit slip vector (parallel to rake)
! r(1) =  dcos(str-pi/2.0d0)*dcos(rak) + dsin(str-pi/2.0d0)*dsin(rak)*dcos(dip)
! r(2) = -dsin(str-pi/2.0d0)*dcos(rak) + dcos(str-pi/2.0d0)*dsin(rak)*dcos(dip)
! r(3) =                                                    dsin(rak)*dsin(dip)
!
! ! Project shear component on slip vector
! shear = s(1)*r(1) + s(2)*r(2) + s(3)*r(3)
!
! ! Coulomb stress (recall sign convention: pos = dilation)
! coulomb = shear + fric*normal
!
! return
! end subroutine calc_tractions

end module elast


text:
module Mod_plcd_PostprocessStrain
   use typre
   use Mod_plcd_BaseElmope
   implicit none
   private
   public SetPointersPostprocessStrain

   type, extends(PointerSetter) :: SPPostprocessStrain
contains
      procedure :: SpecificSet => SpecificSetPostprocessStrain
   end type
   type(SPPostprocessStrain) :: SetPointersPostprocessStrain

contains

   !----------------------------------------------------------------------------
   !Setting Pointers

   subroutine SpecificSetPostprocessStrain(d)
      implicit none
      class(SPPostprocessStrain) :: d

      !Post Process Strain
      if(a%PostprocessStrain) then
         call ConcatenateProcedures(ProcHook%InGaussElmats,SaveStrainToPostprocess)
      endif

   end subroutine

   !Save Strain to PostProcess
   subroutine SaveStrainToPostprocess
      implicit none

      class(PLCDMaterial), pointer :: Mater => NULL()

      call a%ElementMaterialsData(ielem)%p%GetMaterialPointer(Mater)
      call Mater%CT%GetStrain(gradDisp,a%Strain(ielem)%a(:,igaus))
   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