text:
! ###################################################################
! Copyright (c) 2015-2020, Marc De Graef Research Group/Carnegie Mellon University
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without modification, are
! permitted provided that the following conditions are met:
!
!     - Redistributions of source code must retain the above copyright notice, this list
!        of conditions and the following disclaimer.
!     - Redistributions in binary form must reproduce the above copyright notice, this
!        list of conditions and the following disclaimer in the documentation and/or
!        other materials provided with the distribution.
!     - Neither the names of Marc De Graef, Carnegie Mellon University nor the names
!        of its contributors may be used to endorse or promote products derived from
!        this software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
! ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
! LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
! USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
! ###################################################################

!--------------------------------------------------------------------------
! EMsoft:EMPEDkin.f90
!--------------------------------------------------------------------------
!
! PROGRAM: EMPEDkin
!
!> @author Marc De Graef, Carnegie Mellon University
!
!> @brief Kinematical precession electron diffraction dictionary creation
!
!> @date 03/01/15 MDG 1.0 original
!> @date 03/03/15 MDG 1.1 first tests with realistic parameters; reasonable results
!> @date 11/14/15 MDG 1.2 minor name change; verification; HDF5 output
!> @date 05/21/16 MDG 1.3 change for HDF internal file reorganization
!> @date 08/15/19 MDG 1.4 program moved into public repository
!--------------------------------------------------------------------------
program EMPEDkin

use local
use NameListTypedefs
use NameListHandlers
use files
use io
use stringconstants

IMPLICIT NONE

character(fnlen)                  :: nmldeffile, progname, progdesc
type(PEDkinNameListType)          :: pednl

nmldeffile = 'EMPEDkin.nml'
progname = 'EMPEDkin.f90'
progdesc = 'Kinematical Precession Electron Diffraction Dictionary Generation'

! print some information
call EMsoft(progname, progdesc)

! deal with the command line arguments, if any
call Interpret_Program_Arguments(nmldeffile,1,(/ 15 /), progname)

! deal with the namelist stuff
call GetPEDkinNameList(nmldeffile,pednl)

! generate a set of kinematical PED patterns
 call PEDkin_dictionary(pednl,progname, nmldeffile)

end program EMPEDkin

!--------------------------------------------------------------------------
!
! SUBROUTINE:PEDkin_dictionary
!
!> @author Marc De Graef, Carnegie Mellon University
!
!> @brief compute a kinematical PED dictionary
!
!> @param nmlfile namelist file name
!
!> @date 03/02/15 MDG 1.0 original
!> @date 03/06/15 MDG 1.1 testing complete; added EulerAngles.txt output
!> @date 11/14/15 MDG 1.2 added HDF5 support; write patterns using hyperslabbing
!--------------------------------------------------------------------------
subroutine PEDkin_dictionary(pednl,progname, nmldeffile)

use local
use typedefs
use dictmod
use crystal
use initializersHDF
use initializers
use gvectors
use io
use timing
use diffraction
use symmetry
use quaternions
use NameListTypedefs
use constants
use rotations
use so3
use math
use HDF5
use NameListHDFwriters
use HDFsupport
use stringconstants
use error

IMPLICIT NONE

type(PEDkinNameListType),INTENT(INOUT)  :: pednl
character(fnlen),INTENT(IN)             :: progname
character(fnlen),INTENT(IN)             :: nmldeffile

integer(kind=irg)               :: FZcnt, pgnum, ival, i, j, l
type(FZpointd),pointer          :: FZlist, FZtmp, FZtmp2
real(kind=sgl)                  :: la, dval, dmin, glen, gmax, io_real(3), om(3,3), k(3), sgmax, FN(3), xgmin, Ig, Igmax, &
                                   maxint, w, ku(3), kp(3), rnmpp, dx, dy, eu(3), tstart, tstop, x, y, ma, mi
integer(kind=irg)               :: gp(3), imh, imk, iml, nref, gg(3), ix, iy, iz, io_int(5), ww, nsize, tdp, sx, sy, hdferr, &
                                   ninbatch, nbatches, nremainder,ibatch,istat, gridtype, tickstart
integer(HSIZE_T)                :: dims3(3), hdims(3), offset(3)
logical                         :: verbose, insert=.TRUE., overwrite=.TRUE., exists
character(fnlen)                :: groupname, dataset, outname
character(11)                   :: dstr
character(15)                   :: tstrb
character(15)                   :: tstre
character(fnlen, KIND=c_char),allocatable,TARGET :: stringarray(:)
character(fnlen,kind=c_char)                     :: line2(1)
character(fnlen)                :: fname, atype

real(kind=sgl),allocatable      :: pedpattern(:,:), xx(:,:), yy(:,:), line(:), dot(:,:), eulerarray(:,:), sglread(:)
character(len=1),allocatable    :: pedp(:,:,:), pedpat(:,:)
integer(kind=irg),allocatable   :: patinbatch(:)
integer(kind=irg)               :: totnumberbatch, nextra
type(unitcell)                  :: cell
type(DynType),save              :: Dyn
type(gnode),save                :: rlp
type(reflisttype),pointer       :: reflist, nexts, rltmpa

type(HDFobjectStackType)          :: HDF_head

sgmax = 0.50

call timestamp(datestring=dstr, timestring=tstrb)
call Time_tick(tickstart)

!=============================================
!=============================================
! crystallography section
!nullify(cell)
!allocate(cell)

verbose = .TRUE.

call Initialize_Cell(cell,Dyn,rlp,pednl%xtalname, pednl%dmin, pednl%voltage, verbose)

! determine the point group number
j=0
do i=1,32
 if (SGPG(i).le.cell % SYM_SGnum) j=i
end do
pgnum = j
io_int(1) = pgnum
call WriteValue(' Point group number                       : ',io_int, 1, "(I3)")

!=============================================
!=============================================
! generation of all potential reflections inside a reciprocal space sphere
! computed from the camera length and the detector size ...

! first set the maximum |g| value that can possibly give rise to a diffracted beam on the detector (diagonal)
  gmax = sqrt(2.0) * float(pednl%npix) * pednl%rnmpp
  io_real(1) = gmax
  call WriteValue(' Length of longest g-vector               : ', io_real, 1, "(F8.4)")

! this code is taken from the Initialize_ReflectionList routine, but we do not
! need everything from that routine; first get the size of the lookup table
  gp = shape(cell%LUT)
  imh = (gp(1)-1)/4
  imk = (gp(2)-1)/4
  iml = (gp(3)-1)/4

! initialize the reflection list
  nullify(reflist)
  nullify(rltmpa)
  nref = 0

! transmitted beam always has excitation error zero
  gg = (/ 0,0,0 /)
  call AddReflection(rltmpa, reflist, cell, nref, gg)   ! this guarantees that 000 is always the first reflection
  rltmpa%xg = 0.0
  xgmin = 100000.0
  Igmax = 0.0

! now compute |U_g|^2 for all allowed reflections;
ixl: do ix=-imh,imh
iyl:  do iy=-imk,imk
izl:   do iz=-iml,iml
        if ((abs(ix)+abs(iy)+abs(iz)).ne.0) then  ! avoid double counting the origin
         gg = (/ ix, iy, iz /)
         glen = CalcLength(cell, float(gg), 'r' )

! find all reflections, ignoring double diffraction spots
         if ((IsGAllowed(cell,gg)).and.(glen.le.gmax)) then ! allowed by the lattice centering, if any
            call AddReflection(rltmpa, reflist, cell, nref, gg )
! we'll use the sangle field of the rltail structure to store |Ug|^2; we will also need the extinction distance
            rltmpa%sangle = abs(cell%LUT(ix, iy, iz))**2
            if (rltmpa%sangle.gt.Igmax) Igmax = rltmpa%sangle
            rltmpa%xg = 1.0/(abs(cell%LUT(ix,iy,iz))*cell%mLambda)
            if (rltmpa%xg.lt.xgmin) xgmin = rltmpa%xg
         end if ! IsGAllowed
        end if
       end do izl
      end do iyl
    end do ixl

io_int(1) = nref
call WriteValue(' Length of the master list of reflections : ', io_int, 1, "(I8)")

!=============================================
!=============================================
! create the coordinate arrays for the Gaussian peaks that will represent the diffraction spots
rnmpp = 1.0/pednl%rnmpp
ww = 6
tdp = 2*ww+1
allocate(xx(-ww:ww,-ww:ww), yy(-ww:ww,-ww:ww), line(-ww:ww), dot(-ww:ww,-ww:ww))
line = (/ (float(i),i=-ww,ww) /) * rnmpp
xx = spread(line,dim=1,ncopies=2*ww+1)
yy = transpose(xx)

!=============================================
!=============================================
! create the output array
nsize = pednl%npix/2 + ww
allocate(pedpattern(-nsize:nsize,-nsize:nsize))
maxint = Igmax
write (*,*) 'Maximum diffracted intensity : ',maxint

!=============================================
!=============================================
! open the HDF5 output file; in the old code, the output was written into two separate binary files
! here, we do just one output file in the standard EMsoft style
  nullify(HDF_head%next)

! Initialize FORTRAN interface.
  call h5open_EMsoft(hdferr)

! Create a new file using the default properties.
  outname = trim(EMsoft_getEMdatapathname())//trim(pednl%outname)
  outname = EMsoft_toNativePath(outname)
  hdferr =  HDF_createFile(outname, HDF_head)

! write the EMheader to the file
  groupname = SC_PEDkin
  call HDF_writeEMheader(HDF_head, dstr, tstrb, tstre, progname, groupname)

! create a namelist group to write all the namelist files into
  groupname = SC_NMLfiles
  hdferr = HDF_createGroup(groupname, HDF_head)

! read the text file and write the array to the file
  dataset = SC_PEDkinNameList
  hdferr = HDF_writeDatasetTextFile(dataset, nmldeffile, HDF_head)

! leave this group
  call HDF_pop(HDF_head)

! create a namelist group to write all the namelist files into
  groupname = SC_NMLparameters
  hdferr = HDF_createGroup(groupname, HDF_head)
  call HDFwritePEDkinNameList(HDF_head, pednl)

! leave this group
  call HDF_pop(HDF_head)

! then the remainder of the data in a EMData group
  groupname = SC_EMData
  hdferr = HDF_createGroup(groupname, HDF_head)

!=============================================
!=============================================
! rotation sampling section
! we need to get a sampling of orientation space, starting from the
! cubochoric representation or read from file depending on input

if(pednl%sampling .eq. 'dict') then
    nullify(FZlist)
    FZcnt = 0
    gridtype = 0
    call sampleRFZ(pednl%ncubochoric, pgnum, gridtype, FZcnt, FZlist)

! convert the Rodrigues vectors to Euler angle triplets and delete the entries
    allocate(eulerarray(3,FZcnt))
    FZtmp => FZlist                        ! point to the top of the list
    do i=1,FZcnt
        eulerarray(1:3,i) = ro2eu(FZtmp%rod)
        FZtmp => FZtmp%next                  ! point to the next entry
    end do
! delete list
    FZtmp => FZlist%next
    FZtmp2 => FZlist
    do
        deallocate(FZtmp2)
        if (.not. associated(FZtmp) ) EXIT
        FZtmp2 => FZtmp
        FZtmp => FZtmp%next
    end do
    nullify(FZlist)
else if(pednl%sampling .eq. 'file') then
    fname = trim(EMsoft_getEMdatapathname())//trim(pednl%eulerfile)
    inquire(file=fname,exist=exists)
    if(.not.exists) call FatalError('PEDkin:','eulerfile not found : '//trim(fname))
    open(unit=dataunit,file=fname,action='read',status='unknown')

    read(dataunit,*) atype
    read(dataunit,*)FZcnt
    allocate(eulerarray(3,FZcnt))

    if(atype .eq. 'eu') then
        allocate(sglread(3))
        do i = 1,FZcnt
            read(dataunit,*) sglread
            eulerarray(1:3,i) = sglread*cPi/180.0
        end do
    else if(atype .eq. 'qu') then
        allocate(sglread(4))
        do i = 1,FZcnt
            read(dataunit,*) sglread
            eulerarray(1:3,i) = qu2eu(sglread)
        end do
    else
        call FatalError('PEDkin:','Unknown orientation format in eulerfile :'//trim(fname)//&
             '. (only euler angles or quaternions allowed)')
    end if
else
    call FatalError('PEDkin:','unknown sampling type in namelist file :'//trim(nmldeffile))
end if

call Message(' Done creating Euler angle array.')

io_int(1) = FZcnt
call WriteValue(' Number of orientation in dictionary       : ', io_int, 1, "(I8)")

! Euler angle array size
dataset = SC_FZcnt
hdferr = HDF_writeDatasetInteger(dataset, FZcnt, HDF_head)

! and write them to the HDF5 file
dataset = SC_EulerAngles
hdferr = HDF_writeDatasetFloatArray2D(dataset, eulerarray*180.0/sngl(cPi), 3, FZcnt, HDF_head)

! keep the datafile open for writing; we'll collect 1024 patterns in an array,
! and then use hyperslab writing to put them into the HDF5 file.

!=============================================
!=============================================
! and loop over all orientations

! allocate the hyperslab array that will hold the batches of PED patterns
ninbatch = 1024
nbatches = floor(float(FZcnt)/float(ninbatch))
nremainder = mod(FZcnt,ninbatch)

nextra = 0
if(nremainder .gt. 0) nextra = 1
totnumberbatch = nbatches + nextra

allocate(patinbatch(totnumberbatch))
patinbatch = ninbatch
if(nextra .eq. 1) patinbatch(totnumberbatch) = nremainder
allocate(pedpat(-nsize:nsize,-nsize:nsize),stat=istat)

call Message(' Starting main computation loop.')
call Message(' ')
batchloop: do i = 1, totnumberbatch       ! loop over all batches in the eulerarray
! allocate and set the individual output pattern to zero
    if(allocated(pedp)) deallocate(pedp)
    allocate(pedp(pednl%npix,pednl%npix,patinbatch(i)),stat=istat)
    pedp = ' '
    pedpattern = 0.0
    pedpat = ' '
    orientationloop: do l = 1,patinbatch(i)
        pedpattern = 0.0
        pedpat = ' '
        ival = (i-1)*patinbatch(i) + l

! convert the rodrigues vector to a passive rotation matrix.
        om = eu2om(eulerarray(1:3,ival))

! multiplication with (0,0,1) produces the normalized beam direction in a
! cartesian reference frame; so now we can compute the excitation errors
! for every reflection and keep only the ones that are sufficiently small
        k = (/ 0.0, 0.0, 1.0 /)
        ku = matmul(om,k)
        FN = ku
        k = ku/sngl(cell%mLambda)

! first we go through the entire reflection list and compute the excitation errors
! those points that satisfy the cutoff are linked via the nexts pointers
        rltmpa => reflist%next
        nexts => rltmpa
        do j=1,nref
          gg = rltmpa%hkl
          rltmpa%sg = Calcsg(cell,float(gg),k,FN)
! should we consider this point any further ? If so, add it to the strong reflection linked list
          if (abs(rltmpa%sg).le.sgmax) then
            nexts%nexts => rltmpa
            nexts => rltmpa
          end if
          rltmpa => rltmpa%next
        end do

! then, for each point in the nexts list, we compute the components of k' = k+g+s
! and place them in the proper reference frame; we skip the incident beam since it is
! meaningless in the kinematical approximation
      nexts => reflist%next%nexts
      do
! determine the vector k'
        kp = k + float(nexts%hkl) + nexts%sg*ku
        kp = matmul(transpose(om),kp)

! get the intensity for each point
        w = sngl(cPi)*nexts%sg*pednl%thickness
        if (abs(w).lt.1.0e-6) then
          Ig = nexts%sangle  ! * (sngl(cPi)*pednl%thickness/nexts%xg)**2
        else
          Ig = nexts%sangle * (sin(w)/w)**2 ! * (sngl(cPi)*pednl%thickness/nexts%xg)**2
        end if

! determine the spot coordinates on the detector
        x = rnmpp * kp(1)
        y = rnmpp * kp(2)

! and plot that spot as a small Gaussian in the pedpattern array, assuming it falls on the detector.
        if ((abs(x).le.nsize-ww).and.(abs(y).le.nsize-ww)) then
          sx = nint(x)
          sy = nint(y)
          dx = x-sx
          dy = y-sy
          dot = (Ig/Igmax)**0.2 * exp(-((xx-dx)**2+(yy-dy)**2)*0.0025)
          pedpattern(sx-ww:sx+ww,sy-ww:sy+ww) = pedpattern(sx-ww:sx+ww,sy-ww:sy+ww) + dot(-ww:ww,-ww:ww)
        end if

! and repeat this until the end of the list
        if (.not. associated(nexts%nexts)) EXIT
        nexts => nexts%nexts
      end do

      ma = maxval(pedpattern)
      mi = minval(pedpattern)
      pedpattern = ((pedpattern - mi)/ (ma-mi))
      pedpat = char(nint(255.0*pedpattern))

! save the pedpattern to the pedp array
      pedp(1:pednl%npix,1:pednl%npix,l) = pedpat(-nsize+ww+1:nsize-ww,-nsize+ww+1:nsize-ww)

! reset the nexts linked list and start over
      nexts => reflist%next
      rltmpa => nexts%nexts
      do
        nullify(nexts%nexts)
        if (.not. associated(rltmpa%nexts)) EXIT
        nexts => rltmpa
        rltmpa => rltmpa%nexts
      end do
    end do orientationloop

    offset = (/ 0, 0, (i-1)*patinbatch(i) /)
    hdims = (/ pednl%npix, pednl%npix, FZcnt /)
    dims3 = (/ pednl%npix, pednl%npix, patinbatch(i) /)

    if (i.eq.1) then
        dataset = SC_PEDpatterns
        hdferr = HDF_writeHyperslabCharArray3D(dataset, pedp, hdims, offset, dims3, HDF_head)
    else
        dataset = SC_PEDpatterns
        hdferr = HDF_writeHyperslabCharArray3D(dataset, pedp, hdims, offset, dims3, HDF_head, insert)
    end if

    io_real(1) = float(sum(patinbatch(1:i)))*100/float(sum(patinbatch))
    call WriteValue(' Completed ',io_real,1,'(F10.2," % ")')

end do batchloop

call HDF_pop(HDF_head)

! and update the end time
call timestamp(datestring=dstr, timestring=tstre)
groupname = SC_EMheader
hdferr = HDF_openGroup(groupname, HDF_head)

groupname = SC_PEDkin
hdferr = HDF_openGroup(groupname, HDF_head)

! stop time /EMheader/StopTime 'character'
dataset = SC_StopTime
line2(1) = dstr//', '//tstre
hdferr = HDF_writeDatasetStringArray(dataset, line2, 1, HDF_head, overwrite)

tstop = Time_tock(tickstart)
dataset = SC_Duration
hdferr = HDF_writeDatasetFloat(dataset, tstop, HDF_head)

! close the datafile
call HDF_pop(HDF_head,.TRUE.)

! close the Fortran interface
call h5close_EMsoft(hdferr)

end subroutine PEDkin_dictionary

max_stars_repo_name: hwagit/EMsoft

max_stars_repo_licenses: Intel BSD-3-Clause MIT-CMU Unlicense

max_issues_repo_name: hwagit/EMsoft

max_issues_repo_licenses: Intel BSD-3-Clause MIT-CMU Unlicense

max_forks_repo_name: hwagit/EMsoft

max_forks_repo_licenses: Intel BSD-3-Clause MIT-CMU Unlicense


text:
c ********************************************************************
c *                                                                  *
c *  routines to support crack growth by cohesive criterion          *
c *                                                                  *
c ********************************************************************
c
c     ****************************************************************
c     *                                                              *
c     *                subroutine dam_param_cohes                    *
c     *                                                              *
c     *                       written by : rhd                       *
c     *                                                              *
c     *                   last modified : 7/31/2013 rhd              *
c     *                                                              *
c     *     for a killable cohesive element not yet killed,          *
c     *     determine if the element should be killed now.           *
c     *     return values of cohesive parameters to support optional *
c     *     printing of values at start of each step                 *
c     *                                                              *
c     ****************************************************************
c
c
c
      subroutine dam_param_cohes( elem, kill_now, debug, values,
     &                            compute_type )
      use global_data ! old common.main
c
      use main_data,       only : elems_to_blocks
      use elem_block_data, only : history_blocks, eps_n_blocks,
     &                            urcs_n_blocks, history_blk_list
      use damage_data
c
      implicit integer (a-z)
c
c              parameter declarations
c
      logical kill_now, debug
      double precision
     &     values(*)  ! see caller for size
c
c              local declarations
c
      double precision
     &     zero, one, half, t1, t2, tn, d1, d2, dn, fpngp, beta,
     &     deff, deff_peak, teff, normalized_teff, teff_peak,
     &     normalized_deff, peak_normal_stress,
     &     peak_shear_stress, dn_max_hist, ds_max_hist,
     &     dn_limit, ds_limit, dn_at_peak, ds_at_peak,
     &     ratio_normal, ratio_shear, dn_half_peak, ds_half_peak,
     &     kill_fract, deff_ratio, deff_gp, deff_peak_gp
      real chk_value
      double precision,
     &  dimension(:), pointer :: history, urcs_n, eps_n
      logical beyond_peak, fgm_cohes, option_ppr,
     &        option_exponential, kill_criterion_element,
     &        kill_gp, local_debug, option_cavit
      data zero, one, half
     &       / 0.0d00, 1.0d00, 0.5d00 /
c
c
c               elem -- absolute number for interface element
c               kill_now -- output, logical
c               debug -- input, logical
c               values -- output vector, float
c               compute_type -- input integer.
c                    = 1 compute/return the "values" vector
c                        but not kill_now flag (for checking,
c                        printing, etc.)
c                    = 2 compute/return the "values" vector
c                        and the kill_now flag
c
c
c               get basic information about element and the
c               cohesive material model for the element. set pointers
c               into block data structures to simplify access
c               to stress, strain, history values for this specific
c               element.
c
c               cohesive types:
c                  1 -- linear elastic; 2 -- bilinear; 3 -- ramp
c                  4 -- exponential_1; 5 -- exponential_2
c                  6 -- PPR; 7 -- cavit
c                  Only types 1, 4, 6, 7 implemented.
c
      local_debug = .false.
      if( local_debug ) write(out,9030)
      ngp         = iprops(6,elem)
      blk         = elems_to_blocks(elem,1)
      rel_elem    = elems_to_blocks(elem,2)
      hist_size   = history_blk_list(blk)
      hoffset     = (rel_elem-1)*hist_size*ngp
      epsoffset   = (rel_elem-1)*nstr*ngp
      sigoffset   = (rel_elem-1)*nstrs*ngp
      urcs_n      => urcs_n_blocks(blk)%ptr
      eps_n       => eps_n_blocks(blk)%ptr
      history     => history_blocks(blk)%ptr
      cohes_type  = iprops(27,elem)
      option_exponential = cohes_type .eq. 4
      option_ppr         = cohes_type .eq. 6
      option_cavit       = cohes_type .eq. 7
c
c               loop over element integration pts. compute the simple
c               average of the normal traction and (relative)
c               displacement. Also compute average of the two shear
c               tractions and (total) sliding displacement
c
      t1 = zero; t2 = zero; tn = zero; d1 = zero; d2 = zero; dn = zero
c
      do gp = 1, ngp
        t1 = t1 + urcs_n(sigoffset+1)
        t2 = t2 + urcs_n(sigoffset+2)
        tn = tn + urcs_n(sigoffset+3)
        d1 = d1 + eps_n(epsoffset+1)
        d2 = d2 + eps_n(epsoffset+2)
        dn = dn + eps_n(epsoffset+3)
        if( local_debug) then
         i = sigoffset; j = epsoffset
          write(out,9020) gp, urcs_n(i+1), urcs_n(i+2), urcs_n(i+3),
     &                    eps_n(j+1), eps_n(j+2), eps_n(j+3)
        end if
        epsoffset = epsoffset + nstr
        sigoffset = sigoffset + nstrs
      end do
c
c    finish computing the average value for each of the
c    3 tractions and 3 interface displacements. compute
c    resultant shear traction and resultant sliding displacement.
c
c       values        content
c        1         averaged normal traction
c        2         averaged total shear traction (1+2)
c        3         averaged normal opening
c        4         averaged (total) sliding displacement (1+2)
c
c                    for exponential option (=4)
c                    ===========================
c
c        5         averaged single "effective" traction (reflects a
c                  user-specified beta value)
c        6         averaged single "effective" interface displacement (reflects
c                  user-specified beta value)
c        7         user-specified "equivalent" interface displacement
c                  at peak stress
c        8         averaged equivalent stress / user-specified peak stress
c        9         averaged equivalent interface displacement / user-specified
c                  equivalent displacement at peak stress
c
c                    for PPR option (=6)
c                    ===================
c
c        5         averaged normal traction / peak normal traction
c                  on curve
c        6         averaged shear traction / peak shear traction
c                  on curve
c        7         the number of Gauss points as follows
c                   for each GP: get maximum normal displacement
c                   over loading history.
c                   if it exceeds the normal displacement at peak
c                   normal stress on curve, increment counter
c                   of GPs
c        8         the number of Gauss points as follows
c                   for each GP: get maximum sliding displacement
c                   (1+2) over loading history.
c                   if it exceeds the sliding displacement at peak shear
c                   stress on curve, increment counter of GPs
c        9         the number of Gauss points as follows
c                   for each GP: get maximum normal displacement
c                   over loading history.
c                   if it exceeds the 50% of normal displacement at
c                   peak normal stress on curve, increment GP counter
c       10         the number of Gauss points as follows
c                   for each GP: get maximum shear displacement
c                   over loading history.
c                   if it exceeds the sliding displacement at peak shear
c                   stress on curve, increment counter
c       11          maximum value of an "effective" interface
c                   displacement at Gauss points
c                    = sqrt( dnormal**2 + dsliding**2 )
c       12          maximum value of effective peak displacement at Gauss points
c                   corresponding to above value. = sqrt( dn_peak**2 +
c                   ds_peak**2 )
c
c                    for cavit option (=7)
c                    ===================
c
c        5         averaged normal traction / peak normal traction
c                  on curve
c        6         averaged shear traction / peak shear traction
c                  on curve
c
c
      fpngp = one / dble(ngp)
      t1 = t1 * fpngp
      t2 = t2 * fpngp
      tn = tn * fpngp
      d1 = d1 * fpngp
      d2 = d2 * fpngp
      dn = dn * fpngp
      values(1) = tn
      values(2) = sqrt( t1*t1 + t2*t2 )
      values(3) = dn
      values(4) = sqrt( d1*d1 + d2*d2 )
c
c              "exponential_1" option
c
      if( option_exponential ) then
        call dam_param_cohes_exp
        return
      end if
c
c              "PPR" option
c
      if( option_ppr ) then
        call dam_param_cohes_ppr
        return
      end if
c
c              "cavit" option
c
      if( option_cavit ) then
        call dam_param_cohes_cavit
        return
      end if
c
c              unknown cohesive option
c
      write(out,*) ' '
      write(out,*) '>>> FATAL ERROR: dam_param_cohes (1)'
      write(out,*) '    unknown cohesive option. job terminated'
      call die_abort
c
 9020 format(
     &   5x,'gp:                  ',i2,
     & /,5x,'t1:                  ',f12.4,
     & /,5x,'t2:                  ',f12.4,
     & /,5x,'tn:                  ',f12.4,
     & /,5x,'d1:                  ',f12.8,
     & /,5x,'d2:                  ',f12.8,
     & /,5x,'dn:                  ',f12.8 )
 9030 format(/,"..... inside dam_param_cohes ....")
c
c
      contains
c     ========
c     ****************************************************************
c     *                                                              *
c     *                subroutine dam_param_cohes_exp                *
c     *                                                              *
c     *                       written by : rhd                       *
c     *                                                              *
c     *                   last modified : 7/31/2013 rhd              *
c     *                                                              *
c     ****************************************************************
c
      subroutine dam_param_cohes_exp
      implicit integer (a-z)
c
      beta        = props(23,elem)
      fgm_cohes   = abs(props(34,elem)) .gt. 0.1
      if( beta .eq. zero ) then
        values(4) = zero
        values(5) = tn
        values(6) = dn
      else
        if( dn .lt. zero ) then
           tn = zero
           dn = zero
        end if
        values(5) = sqrt( (t1*t1 + t2*t2)/beta/beta + tn*tn )
        values(6) = sqrt( beta*beta*(d1*d1+d2*d2) + dn*dn )
      end if
      teff      = values(5)
      deff      = values(6)
c
      values(7) = props(22,elem)  ! equivalent separation distance
                                  ! at peak stress
      teff_peak = props(13,elem)  ! peak normal stress on expon curve
      deff_peak = props(22,elem)
      if( fgm_cohes ) then
         values(7) = props(35,elem) ! peak normal stress on curve (fgm)
         teff_peak = props(37,elem) ! critical stress ductile (fgm)
         deff_peak = props(35,elem)
      end if
c
      normalized_teff = teff / teff_peak
      normalized_deff = deff / deff_peak
      values(8)       = normalized_teff
      values(9)       = normalized_deff
c
c              does element need to be killed?
c
      kill_now    = .false.
      if( compute_type .ne. 2 ) return
      deff_peak   = props(22,elem)   ! separation distance at peak on
                                     ! exponential curve
      if( fgm_cohes )  deff_peak = props(35,elem)
      deff        = values(6)
      beyond_peak =  deff .gt. deff_peak
      if ( .not. beyond_peak ) return
      if ( normalized_deff .gt. critical_cohes_deff_fract )
     &      kill_now = .true.
c
      return
      end subroutine dam_param_cohes_exp
c
c     ****************************************************************
c     *                                                              *
c     *                subroutine dam_param_cohes_ppr                *
c     *                                                              *
c     *                       written by : rhd                       *
c     *                                                              *
c     *                   last modified : 7/31/2013 rhd              *
c     *                                                              *
c     ****************************************************************
c
      subroutine dam_param_cohes_ppr
      implicit integer (a-z)
c
c              loop over integration pts of this interface element.
c               o pull from history the maximum normal and sliding
c                 displacements ever reached over loading history
c               o pull from history the normal and sliding displacements
c                 at which tractions degrade to zero
c               o pull from history the fractions of displacements
c                 at peak values of normal and shear tractions
c               o count number of gauss points at which maximum
c                 interface displacements over history exceed 50% of
c                 values at peak tractions
c               o count number of gauss points at which maximum
c                 interface displacements over history exceed values
c                 at peak tractions
c
c              Does the integration pt status meet criterion for element
c              extinction? (All GP must meet criterion for extincting
c              the element).
c
c              After much discussion, Kyoungsoo recommends using the
c              "or" condition on normal and shear loading. If either
c              the normal opening or sliding displacement exceeds a
c              fraction (e.g., 90%) of the corresponding displacement
c              when the traction degrades fully to zero, kill the element
c              (must be satisfied at all interface element Gauss points)
c
c              The reason is that the remaining shear capacity
c              for example, degrades severely as the normal
c              traction capacity degrades. Same for normal capacity
c              under shear loading.
c
      epsoffset        = (rel_elem-1)*nstr*ngp
      num_dn_past_peak = 0
      num_ds_past_peak = 0
      num_dn_half_peak = 0
      num_ds_half_peak = 0
      kill_fract = ppr_kill_displ_fraction
      kill_criterion_element = .true.
      deff_ratio = -1.0 ! initialize to small value
      values(11) = zero
      values(12) = one
      if( local_debug ) write(out,9000) elem
c
      do gp = 1, ngp
c
c                  values for possible printing
c
        dn_limit    = history(hoffset+1)
        ds_limit    = history(hoffset+2)
        dn_max_hist = history(hoffset+3)
        ds_max_hist = abs(history(hoffset+4))
        ratio_normal = props(37,elem)
        ratio_shear  = props(39,elem)
        dn_at_peak = dn_limit * ratio_normal
        ds_at_peak = ds_limit * ratio_shear
        dn_half_peak = half * dn_at_peak
        ds_half_peak = half * ds_at_peak
        if( dn_max_hist .gt. dn_at_peak )
     &                       num_dn_past_peak = num_dn_past_peak + 1
        if( ds_max_hist .gt. ds_at_peak )
     &                       num_ds_past_peak = num_ds_past_peak + 1
        if( dn_max_hist .gt. dn_half_peak )
     &                       num_dn_half_peak = num_dn_half_peak + 1
        if( ds_max_hist .gt. ds_half_peak )
     &                       num_ds_half_peak = num_ds_half_peak + 1
c
c                  does this Gauss point satisfy criterion to
c                  kill element? If not, element is not killable
c                  at this load step
c
        kill_gp = ( dn_max_hist .gt. kill_fract * dn_limit )  .or.
     &            ( ds_max_hist .gt. kill_fract * ds_limit )
        if( .not. kill_gp ) kill_criterion_element = .false.
c
c                  values that may be used by adapative load step
c                  processor
c
        d1 = eps_n(epsoffset+1)
        d2 = eps_n(epsoffset+2)
        dn = eps_n(epsoffset+3)
        epsoffset = epsoffset + nstr
        deff_gp = sqrt( dn**2 + d1**2 + d2**2 )
        deff_peak_gp = sqrt( dn_at_peak**2 + ds_at_peak**2 )
        normalized_deff = deff_gp / deff_peak_gp
c        write(*,*) '...gp, d1, d2, dn, deff_gp, dn_at_peak,ds_at_peak,
c     &     deff_peak_gp'
c        write(*,*) gp, d1, d2, dn, deff_gp, dn_at_peak,ds_at_peak,
c     &       deff_peak_gp
        if( normalized_deff .gt. deff_ratio ) then
            values(11) = deff_gp
            values(12) = deff_peak_gp
            deff_ratio = normalized_deff
        end if
c
        if( local_debug .and. gp .eq. 1 )
     &       write(out,9010) gp, dn_limit, ds_limit,
     &       dn_max_hist, ds_max_hist, ratio_normal, ratio_shear,
     &       dn_at_peak, ds_at_peak, dn_half_peak, ds_half_peak,
     &       kill_gp, deff_gp, normalized_deff
        hoffset   = hoffset + hist_size
      end do
c
      peak_normal_stress = props(13,elem)
      peak_shear_stress  = props(14,elem)
      values(5) = zero; values(6) = zero
      if( peak_normal_stress .gt. zero )
     &    values(5) = values(1) / peak_normal_stress
      if( peak_shear_stress .gt. zero )
     &    values(6) = values(2) / peak_shear_stress
c
      values(7) = zero; values(8) = zero
      values(9) = zero; values(10) = zero
      if( num_dn_past_peak .gt. 0 ) then
         values(7) = dble(num_dn_past_peak)
      end if
      if( num_ds_past_peak .gt. 0 ) then
         values(8) = dble(num_ds_past_peak)
      end if
c
      if( num_dn_half_peak .gt. 0 ) then
         values(9) = dble(num_dn_half_peak)
      end if
      if( num_ds_half_peak .gt. 0 ) then
         values(10) = dble(num_ds_half_peak)
      end if
c
c              does element need to be killed? needs fixing ....!
c
      kill_now    = .false.
      if ( compute_type .ne. 2 ) return
      kill_now =  kill_criterion_element
      return
c
 9000 format('**** debug ppr in dam_param_cohes element: ',i6)
 9010 format(
     &   5x,'gp:                  ',i2,
     & /,5x,'dn_limit:            ',f12.8,
     & /,5x,'ds_limit:            ',f12.8,
     & /,5x,'dn_max_hist:         ',f12.8,
     & /,5x,'ds_max_hist:         ',f12.8,
     & /,5x,'ratio_normal:        ',f12.8,
     & /,5x,'ratio_shear:         ',f12.8,
     & /,5x,'dn_at_peak:          ',f12.8,
     & /,5x,'ds_at_peak:          ',f12.8,
     & /,5x,'dn_half_peak:        ',f12.8,
     & /,5x,'ds_half_peak:        ',f12.8,
     & /,5x,'kill_gp:             ',l1,
     & /,5x,'deff_gp:             ',f12.8,
     & /,5x,'normalized_deff:     ',f12.8)

c
      end subroutine dam_param_cohes_ppr
c
c     ****************************************************************
c     *                                                              *
c     *                subroutine dam_param_cohes_cavit              *
c     *                                                              *
c     *                       written by : rhd                       *
c     *                                                              *
c     *                   last modified : 7/31/2013 rhd              *
c     *                                                              *
c     ****************************************************************
c
      subroutine dam_param_cohes_cavit
      implicit integer (a-z)
      double precision
     &     tn_max_over_history, dn_max_at_tn_max
c
c              loop over integration pts of this interface element.
c
c               o pull from history the maximum normal traction
c                 reached over loading history
c               o count number of integration pts at which the
c                 current normal tractions is less than 20% of the
c                 peak value
c
c              For element death, all integration points must have
c              have a normal traction below 20% of peak value.
c
      kill_fract = 0.20   !  make input parameter later
      kill_criterion_element = .true.
      if( local_debug ) write(out,9100) elem
      hoffset     = (rel_elem-1)*hist_size*ngp
      epsoffset   = (rel_elem-1)*nstr*ngp
      sigoffset   = (rel_elem-1)*nstrs*ngp
c
      count_past_peak = 0
c
      do gp = 1, ngp
c
c                  values for possible printing
c
        tn = urcs_n(sigoffset+3)
        dn = eps_n(epsoffset+3)
        tn_max_over_history =  history(hoffset+12)
        dn_max_at_tn_max    =  history(hoffset+13)
        epsoffset = epsoffset + nstr
        sigoffset = sigoffset + nstrs
        hoffset   = hoffset + hist_size
c
c                  does this point satisfy criterion to
c                  kill element? If not, element is not killable
c                  at this load step
c
        kill_gp =  dn .gt. dn_max_at_tn_max  .and.
     &             tn .le. kill_fract * tn_max_over_history
        if( dn .gt. dn_max_at_tn_max )
     &      count_past_peak = count_past_peak + 1
        if( .not. kill_gp ) kill_criterion_element = .false.
c
        if( local_debug .and. gp .eq. 1 )
     &       write(out,9010) gp, tn, dn, tn_max_over_history,
     &                       dn_max_at_tn_max
      end do
c
      values(5) = dble( count_past_peak )
c
c              does element need to be killed? needs fixing ....!
c
      kill_now    = .false.
      if ( compute_type .ne. 2 ) return
      kill_now =  kill_criterion_element
      return
c
 9100 format('**** debug cavit in dam_param_cohes element: ',i6)
 9010 format(
     &   5x,'gp:                  ',i2,
     & /,5x,'tn:                  ',f12.8,
     & /,5x,'dn:                  ',f12.8,
     & /,5x,'tn_max_over_history: ',f12.8,
     & /,5x,'dn_max_at_tn_max:    ',f12.8 )
c
      end subroutine dam_param_cohes_cavit
c
      end subroutine dam_param_cohes
c     contains completed
c     ==================
c
c     ****************************************************************
c     *                                                              *
c     *                 subroutine dam_print_elem4                   *
c     *                                                              *
c     *                       written by : rhd                       *
c     *                                                              *
c     *                   last modified : 11/28/2010 rhd             *
c     *                                                              *
c     *     This routine prints out the status of cohesive elements  *
c     *     marked as killable at the beginning of a load step       *
c     *                                                              *
c     ****************************************************************
c
      subroutine dam_print_elem4( step, iter )
      use global_data ! old common.main
      use elem_extinct_data, only : dam_state, dam_print_list
      use main_data, only : output_packets, packet_file_no
      use damage_data
      implicit integer (a-z)
c
c                    declare local variables
c
      double precision
     &    dummy, eps_plas, eps_crit, sig_mean, sig_mises, eps_plas_tol,
     &    d_eps_plas, max_d_eps_plas, values(20),half
      character(len=10) :: special_char
      logical ldummy, debug, all_killed, beyond_peak,
     & found_fgm_cohes, found_homog_cohes, found_ppr,
     & found_exponential, option_ppr, option_exponential,
     & write_header, show_results, option_cavit, found_cavit
      data zero, half, debug, eps_plas_tol / 0.0, 0.5, .false., 1.0e-9 /
c
c          1.0  print the header
c
      write(out,*) ' '
      write(out,*) ' ********************************************* '
      write(out,*) ' ***       Cohesive element status         *** '
      write(out,*) ' ********************************************* '
      write(out,*) ' '
c
c         2.0  check to see if all elements in print list have been killed.
c              if so, print a message and return
c
      all_killed = .true.
      do elem_loop = 1, num_print_list
         element  = dam_print_list(elem_loop)
         elem_ptr = dam_ptr(element)
         if( dam_ptr(element) .eq. 0 ) cycle
         if( dam_state(elem_ptr) .ne. 0 ) cycle
         all_killed = .false.
      end do
c
      if( all_killed ) then
         write(out,9040)
         return
      end if
c
c         3.0  process each element in the list of elements specifed
c              by the user. skipped already killed eleemnts.
c
      local_count       = 0
      found_homog_cohes = .false.
      found_fgm_cohes   = .false.
      found_ppr         = .false.
      found_exponential = .false.
      found_cavit       = .false.
      write_header      = .true.
c
c          loop over all model elements
c
      do elem_loop = 1, num_print_list
         element  = dam_print_list(elem_loop)
         elem_ptr = dam_ptr(element)
c
c             check if element is a killable and/or if it has
c             already been killed.
c
         if ( dam_ptr(element) .eq. 0 ) cycle
         if ( dam_state(elem_ptr) .ne. 0 ) cycle
c
c             cohesive element is still active in the model. call a
c             routine to compute a vector of values to print for
c             this element. these are average values computed from
c             the element gauss point values.
c
         call dam_param_cohes( element, ldummy, debug, values, 1 )
         cohes_type         = iprops(27,element)
         option_exponential = cohes_type .eq. 4
         option_ppr         = cohes_type .eq. 6
         option_cavit       = cohes_type .eq. 7
         if( option_exponential ) found_exponential = .true.
         if( option_ppr )         found_ppr = .true.
         if( option_cavit )       found_cavit = .true.
c
c              write results for exponential model
c
         if( option_exponential ) then
           if( write_header ) then
              write(out,9030)
              write_header = .false.
           end if
           beyond_peak = values(6) .gt.  values(7)
           special_char(1:10) = ' '
           if ( beyond_peak ) special_char(1:3) = '(*)'
           if ( values(6) .gt. half*values(7) ) then
             write(out,9004) element, (values(k),k=1,6),
     &                       special_char(1:3)
             local_count = local_count + 1
             if( abs(props(34,element)) .gt. 0.1 ) then
                 found_fgm_cohes = .true.
             else
                 found_exponential = .true.
             end if
           end if
           cycle
         end if
c
c              write results for ppr model
c
         if( option_ppr ) then
           if( write_header ) then
              write(out,9032)
              write_header = .false.
           end if
           num_gp_post_peak_normal = int(values(7))
           num_gp_post_peak_shear  = int(values(8))
           num_gp_half_peak_normal = int(values(9))
           num_gp_half_peak_shear  = int(values(10))
           show_results =
     &       (num_gp_half_peak_normal + num_gp_half_peak_shear) .gt. 0
           if( show_results ) then
              write(out,9005) element, (values(k),k=1,6),
     &                num_gp_post_peak_normal, num_gp_post_peak_shear
              local_count = local_count + 1
           end if
           cycle
         end if
c
c              write results for cavit model
c
         if( option_cavit ) then
           if( write_header ) then
              write(out,9034)
              write_header = .false.
           end if
           write(out,9008) element, (values(k),k=1,4),
     &                     int(values(5))
           cycle
         end if
c
         write(out,9100)
         call die_gracefully
      end do
c
c         4.0  write nice ending explanations for table of status
c              values just printed
c
      write(out,*) ' '
      if( found_exponential ) then
           if( local_count .gt. 0 ) write(out,9050)
           if( local_count .eq. 0 ) write(out,9052)
      end if
      if( found_ppr ) then
         if( local_count .gt. 0 ) write(out,9055)
         if( local_count .eq. 0 ) write(out,9058)
      end if
      if( found_cavit ) then
         write(out,9060)
      end if
      write(out,*) ' '
c
c         5.0 for packet output, we now know the number of data lines
c             (elements) to be output. run thru the above loop again to
c             write packets.
c
      if ( local_count .eq. 0 ) go to 9000
      if ( .not. output_packets ) go to 9000
      write(packet_file_no) 7, local_count, step, iter
      do elem_loop = 1, num_print_list
         element  = dam_print_list(elem_loop)
         elem_ptr = dam_ptr(element)
c
c             check if element is a killable element and/or if it has
c             already been killed.
c
         if ( dam_ptr(element) .eq. 0 ) cycle
         if ( dam_state(elem_ptr) .ne. 0 ) cycle
c
c             cohesive element is still active in the model. call a
c             routine to compute a vector of values to issue values into
c             binary file
c
         call dam_param_cohes( element, ldummy, debug, values, 1 )
         cohes_type  = iprops(27,element)
         option_exponential = cohes_type .eq. 4
         option_ppr         = cohes_type .eq. 6
         if( option_exponential ) then
          beyond_peak = values(6) .gt. values(7)
          special_char(1:3) = ' '
          if ( beyond_peak ) special_char(1:3) = '(*)'
          if ( values(6) .gt. 0.5*values(7) )
     &     write(packet_file_no) element, cohes_type, (values(k),k=1,6),
     &                           special_char
         end if
         if( option_ppr ) then
           num_gp_post_peak_normal = int(values(7))
           num_gp_post_peak_shear  = int(values(8))
           num_gp_half_peak_normal = int(values(9))
           num_gp_half_peak_shear  = int(values(10))
           show_results =
     &       (num_gp_half_peak_normal + num_gp_half_peak_shear) .gt. 0
           if( show_results )
     &        write(packet_file_no) element, cohes_type,
     &           (values(k),k=1,6),
     &           num_gp_post_peak_normal, num_gp_post_peak_shear
         end if
      end do
c
 9000 continue
      return
c
 9040 format(1x,'*** NOTE: all elements in the killable element',
     & ' status list have been killed.',/)
 9030 format(1x,'element       Tn             Ts            Dn    ',
     &       '        Ds           Teff       ',
     &       '  Deff    ',/,1x,
     &       '-------   --------       ---------      --------',
     &       '     ---------    ----------- ',
     &       '  -----------')
 9032 format(1x,'element    T-normal       T-shear        D-normal ',
     &       '    D-shear     Tn/Tn-peak   ',
     &       ' Ts/Ts-peak   ',/,1x,
     &       '-------   --------       ---------      --------',
     &       '     ---------    ----------- ',
     &       '  -----------')
 9034 format(1x,'element    T-normal       T-shear        D-normal ',
     &       '    D-shear     (1)',/,1x,
     &       '-------   --------       ---------      --------',
     &       '     ---------    ---- ' )
 9050 format(
     &   1x,'Exponential option: **** only elements with Deff > 0.5 x ',
     &      'Dpeak listed ****',
     & /,1x,'(*) denotes deformation beyond peak stress',
     & /,1x,'**** T and D values denote independent averages',
     &      ' computed over element integration points ****',
     & /,6x,'(T & D values may thus not exactly satisfy cohesive',
     &      ' constitutive model)' )
 9052 format(
     &   1x,'Exponential option: **** no elements with Deff > 0.5 x ',
     &      'Dpeak ****' )
 9055 format(
     &    1x,'PPR option: ',
     &       'T -> traction values, D -> relative displacement values',
     &  /,1x,'**** only elements with at least 50% peak loading ',
     &      'listed ****',
     &  /,1x,'(#,#) denotes number of Gauss points for (normal,shear)',
     &       ' displacement that have exceeded the coresponding',
     &  /,1x,'      values at peak (normal,shear) stress over the',
     &       ' loading history',
     &  /,1x,'**** T and D values shown denote independent averages',
     &       ' computed over element integration points ****',
     &  /,6x,
     & '(these T & D values may thus not exactly satisfy cohesive',
     &       ' constitutive model)' )
 9058 format(
     &    1x,'PPR option: ',
     &       '**** no elements with at least 50% peak loading **** ')
 9060 format(
     &   1x,'Cavit option: integration point average tractions & ',
     &      'displacement jumps',
     &  /,1x,'(1) Number of integration points with post-peak state ',
     &         'based on T-normal, D-normal')
 9004 format(1x,i7,1x,e14.6,1x,e14.6,
     &       e14.6, e14.6, e14.6, e14.6, 1x,a3)
 9005 format(1x,i7,1x,e14.6,1x,e14.6,
     &       e14.6, e14.6, f11.5, f13.5, 3x,'(',i1,',',i1,')')
 9008 format(1x,i7,1x,e14.6,1x,e14.6,
     &       e14.6, e14.6, 3x, I1)
 9100 format(//,'>>> FATAL ERROR: dam_print_elem4'
     &        /,'                 job terminated' )
c
      end

ext: f


text:
${ARBITRATION}/rtl/majority3.v
${MISC}/rtl/dsync.v
${STORAGE}/rtl/fifo.v
${STORAGE}/rtl/fifo_fwft.v
${UART}/rtl/uart_rx.v
${UART}/rtl/uart_rx_fifo.v
${UART}/rtl/uart_tx.v
${UART}/rtl/uart_tx_fifo.v
${UART}/rtl/uart.v
${UART}/test/testbench/uart_fifo_loopback.sv


text:
      DOUBLE PRECISION FUNCTION MYSIN(A)
      DOUBLE PRECISION A
      MYSIN = DSIN(A)
      RETURN
      END

max_stars_repo_name: blackrez/gonum

max_issues_repo_name: mingrammer/gonum

max_forks_repo_name: mingrammer/gonum


text:
fl
\
\
\ 2* ( n1 -- n1<<1 ) n2 is shifted logically left 1 bit
[ifndef 2*
: 2* _xasm2>1IMM h0001 _cnip h05F _cnip ;
]

: sym
	cr
	lastnfa
	begin
\ (nfa -- )
		c" $C_" over swap npfx
\ (nfa flag -- )
		over c" $V_" npfx or

		if
			." [ifndef " dup .strname cr
			base W@ ."     h" over nfa>pfa 2+ W@ hex . base W! ." wconstant " dup .strname cr
			." ]" cr cr
		then
		nfa>next dup 0=
	until
	drop cr
;

\ _wordArray ( n cstr -- addr)
: _wordArray
	lockdict ccreate $C_a_dovarw w, here W@ 0 w, swap 0 max h1000 min dup w, 2* allot forthentry freedict
	dup dup W@ swap 2+ swap 2* 0 fill
;

\ _waAppend( addr val -- )
: _waAppend
	over dup W@ swap 2+ W@ <
	if
		over dup W@ 2+ 2* +  W!
		1 swap W+!

	else
		hDD ERR
	then
;

\ _string ( cstr name -- addr)
: _string
	lockdict ccreate here W@ 2+ swap
	$H_cq w, dup here W@ ccopy C@ 1+ allot herewal
	$C_a_exit w,
	forthentry freedict
;

\ cstr> ( addr1 addr2 -- )
: cstr>
	0 rot2
\ ( 0 a1 a2 -- )
	2dup C@ swap C@ swap
\ ( 0 a1 a2 l1 l2 -- )
	2dup
	>
	>r
	min  rot
\ ( 0 a2 length  a1 -- )
	1+ rot 1+ rot
\ ( 0 a1+1 a2+1 length -- )
	bounds
	do
\ ( 0 a1+n -- )
		dup C@ i C@
		-
		rot + tuck
		0<>
		if
			leave
		then
		1+
	loop
	drop
	dup
	if
		r> drop 0>
	else
\ strings are the same for length compared, if a1 is longer than a2 then true
		drop r>
	then
;

\ _strBsort ( addr -- addr t/f )
: _strBsort
	-1 swap
	dup W@ 1 >
	if
		dup W@ 1- 2* over 4+ swap bounds
		do
			i W@ i 2+ W@
			cstr>
			if
				i W@ i 2+ W@
				i W! i 2+ W!
				nip 0 swap
			then
		2 +loop
	then
	swap
;

\ strBsort ( addr -- addr ) bubble sort of string array
: strBsort
	begin
		_strBsort
	until
;

\ symgen ( -- )
: symgen
	lockdict
	h400 c" _wlist" _wordArray
	lastnfa
	begin
\ (array nfa -- )
\ (nfa -- )
		c" $C_" over swap npfx
\ (nfa flag -- )
		over c" $V_" npfx or
		over c" $S_" npfx or

		if
			dup tbuf namecopy
			tbuf C@ h1F and tbuf C!
			over tbuf c" _$$" _string
			_waAppend
		then
		nfa>next dup 0=

	until
	drop
	strBsort
	dup W@ 0<>
	if

		cr cr cr
		dup W@ 2* over 4+ swap bounds
		do
			i W@
			." [ifndef " dup .cstr cr
			base W@ ."     h"
			over find
			if
				execute
			else
				drop 0
			then
			hex . base W! ." wconstant "  .cstr cr
			." ]" cr cr

		2 +loop
	then
	drop
\	c" _wlist" (forget)
	freedict
;


text:
The Davis Creamery is an ice cream parlor.  Until early 2011, it was located at Oakshade Town Center, a few doors down from Dos Coyotes. That location used to be a Marble Slab Creamery. The Davis Creamery now resides inside of Sugar Daddies in downtown.

Former News Headlines for Davis Creamery lovers, from Wendy Weitzels Comings and Goings column in the Enterprise:

 Davis Creamery, also home to Cupcake Craving and Cookie Connection, is moving from South Davis to downtown, effective in May. The owner has leased a space at 113 E St., in The Lofts, last occupied by Yogurt Shack.

 It will serve ice cream, cookies and cupcakes. Once it moves, plans are to be open seven days a week, from 11 a.m. until 2 a.m., according to a news release from Davis Commercial Properties, which handled the lease.

Regular ice cream flavors included chocolate, blueberry, vanilla brown sugar, peanut butter, salted caramel, chai tea, chocolate mint, banana, and honey.  They also offered also seasonal flavors, such as pumpkin and egg nog.  In December 2007, liver and onions was listed as a flavor, although it was conspicuously absent from the display case.  Mixins could be added to the basic flavors, with the mixing being done on the nowcommon cold slab.

The Davis Creamery made its own ice cream in the shop, and handmade cookies and cupcakes are also available.  Their motto is Udderly Delicious, and their logo features a cows cow riding a bicycles bicycle.

As of February 2009, they also offered cupcakes from http://www.cupcakecraving.com Cupcake Craving. There were a wide variety of cupcake flavors such as Bavarian Dream, Monkey Mania, Caramel Comfort, Raspberry Reward just to name a few. Different combo deals for ice cream and cupcakes are available like the single combo which is $5 for any cupcake and a scoop of ice cream. Starting in July 2009, the Davis Creamery began carrying cookies from http://www.cookieconnection.com Cookie Connection out of Sacramento, which are also wiki:woodland:The Cookie Company cookies from Woodland. They now offer about 24 different varieties of gourmet cookies, and the cowpies (ice cream sandwiches) are available both premade with these cookies and kept frozen, or made to order to fit to your specific tastes.

They have also been spotted at the Farmers Market Farmers Market, serving up some of their specialty flavors, such as avocado coconut, honey vanilla, peach brown sugar, lemon basil, and many more. The farmers market flavors are made from fresh ingredients purchased from local farmers, so the flavors vary based on what is available for purchase at the market. These flavors are available in the store, but only upon request as they are kept in the back for the market.

They want you to join the http://www.thedaviscreamery.com Davis Creamery online or follow their http://www.daviscreamery.blogspot.com/ blog, http://twitter.com/daviscreamery Twitter, and http://www.facebook.com/pages/DavisCA/DavisCreamery/22741963831 Facebook to hear about the newest flavors, biggest deals, and special sales.

Stacy Daniels took over the Davis Creamery from David Robert (former owner of the Marble Slab Creamery in this location) on May 1, 2010.

 Former Daily Specials

The also give out blue coupons for $1 off a regularly priced quart of ice cream.

  Gotta love ice cream, gelato, or frozen yogurt on a hot Summer day!

20071223 18:45:15 nbsp You can tell youre getting the real stuff here.  For example, the chocolate peanut butter ice cream has chunks of peanuts in it, and the banana ice cream  a flavor I dont normally like  actually tastes like bananas.  There isnt a huge supply of mixins, but they have the standard ones; dont overlook the homemade brownies off to the side.  Now I just need to save up some calories for a Cow Pie. Users/CovertProfessor

20080127 16:03:08 nbsp Ive tried their Cow Pie. The Vanilla ice cream was delicious and had good flavor and consistency (not too hard, not too soft). However, the Cow Pie was more than I could handle because I wasnt able to finish the second cookie! I do wish they would make their cookies chewy, instead of crunchy, so that it would be easier to eat along with the ice cream. I am looking forward to going back and trying their other flavors. Users/at86

20080201 10:49:55 nbsp how are the vegan options here? mixed berry sorbet sounds cool... Users/PxlAted

20080306 22:50:12 nbsp Saw a coupon in the California Aggie for $1 off... also saw that students get a 10% discount with ID? I guess it wouldnt hurt to ask! Users/at86

20080313 16:46:59 nbsp The Cupcake ice cream was great, but the blueberries were bad.  It was obvious they were frozen, not fresh, and were hard in the center.  Im also not a big fan of frozen strawberries in syrup which they had in lieu of fresh strawberries.  Next time Ill pass on the fruit toppings.  Users/GraceAn

20080318 15:34:09 nbsp Had just cones today, a chocolate and a raspberry.  Nothing totally unique but nearby and good quality.  Will be back for sure and likely going to give the Cow pie a try. Users/loneshark

20080326 17:52:08 nbsp Have to say this place has the best ice cream hands down. In response to at86, Ive heard that theyve changed their cookie making processes and it shows. Every time Ive had one, the cookies have been perfect. I recommend the cowpie to absolutely everyone! Oh... and be sure to ask for the chocolate chip peanut butter cookies, best kept secret of the store!    Users/blitzcatt

20080428 18:11:12 nbsp They had great customer service and fantastic ice cream! Love the Chocolate Peanut Butter! Definitely recommend it for summer time blues  Users/LizDarko

20080501 10:58:20 nbsp I had a good experience with Davis Creamery a couple months ago: good ice cream and good service. Then I went there with friends yesterday and was very disappointed at the service. Even though we had a coupon that explicitly said $1 from the Aggie with no fine print on it, they wouldnt let us use it on single and junior scoop, frozen yogurt (basically everything we wanted to get!). I dont know why they would put a coupon in the newspaper if theyre just going to add in all these stipulations last minute. Its a shame because theyre ice cream is pretty good quality. Well, they lost our business that day and Im sure we wont be going back very soon in the future either.  Users/AmberHMa

20080504 19:27:41 nbsp Best ice cream in Davis without a doubt.  Ice cream is by far my favorite food (how could anything else possibly be better!!) and after my first taste of DC ice cream I was hooked.  Try it yourself just once, like CovertProfessor said, you will instantly be able to tell you are getting the real stuff.  The ice cream here is supposedly made in Davis and they have a huge selection of things to order.  You can get Froyo (the natural kind), ice cream + mixins, Cow Pies (ice cream sandwiched between two cookies), milkshakes, and other menu items I neglect because I always get vanilla ice cream + strawberry/raspberry ice cream + strawberry toppings.  Try the milkshakes though they are excellent!  Their milkshakes have that heavy ice cream flavor and yummy ice cream chunks you might be familiar with if youve tried some milkshakes from upscale restaurants.  Another thing I like about Davis Creamery is that it has some interesting flavors... for example yesterday they had blueberry flavored ice cream.  Never seen that flavor before and it was a nice rich purple color.  Cool.
I made this daviswiki account just so I could comment on Davis Creamery.  Yes, I know it is a bit sad.  Cold Stone just closed down and I hope my beloved Davis Creamery doesnt follow the same fate.
Oh yea by the way, Davis Creamery has very reasonable prices.  3.95 for two ice creams and a mixin (thats without coupons/student id also) and every time I go my cup is overflowing from the top.  Cant give this place a better recommendation.  Go out there and try it, guys!

P.S. About GraceAns comments on DCs use of frozen strawberries as toppings  I agree.  The frozen toppings are no dice and my only gripe about Davis Creamery.  They taste alright on the outside but the inside of the fruit is still frozen and it doesnt go too well with the ice cream.  Whats up with that, DC!   Users/DavidLiu

Hi Guys!
Thank you for taking the time to provide comments and feedback about DC.  Were glad to see many of you enjoying our ice cream. We feel weve got the best stuff in town! :D
We also appreciate your patience as we continue to figure some things out, like balancing fresh fruit vs. frozen fruit and training our team members to provide consistent and excellent customer service every single time.
CP Thanks for getting the ball rolling!
at86 & PixlAted  We were playing around with our cookie recipe when we opened, (you wouldnt believe how challenging it was to find a quality cookie that could be baked consistently by whomever was working that day), I think weve found the right one now!
The only vegan or nondairy option we have at the moment is the mixedberry sorbet.  Give it a try with a granola mixin... I think youll enjoy it.
GraceAn I totally agree with you on the fresh vs. frozen strawberries. Were committed to fresh strawberries as the choice for our topping since we believe our ice cream deserves it. We dont want to blow it for the sake of convenience with frozen berries.
blitzcatt, LizDarko  Thanks for the great comments!
AmberHMa Im glad that we managed to not disappoint you with our ice cream!  We appreciate your feedback and hope youll give us another shot. Weve had a chat with our team members about honoring our coupons and have worked on revising our Ads to make them clear and enticing.  How about you come by and try out Spicy Chai Tea flavor?
DavidLiu  Thanks for the great enthusiasm!  If we can make everyone feel the way you do, well be sticking around for a long time to come.
Were committed to keep on trying our best and bringing you locally produced, innovative ice cream flavors that are delicious!  Check us out on FB and join the DC Fan Club!

20080611 22:19:54 nbsp I LOVE DAVIS CREAMERY! Users/davisrox

20080703 20:44:40 nbsp I just had the Spicy Chai Tea ice cream and it is excellent!  The other flavors are all great too, but the Chai is perfect.  Every time Ive been into DC, the customer service has been great too! Users/RMC

20080809 20:38:38 nbsp Service and ice cream were great. Cookie was still a little crisp, though not enough to deter my next first. Great place. Will be back. Users/Sar58

20080904 18:07:55 nbsp Its listed for sale on craigslist. Users/BradBenedict

20090201 16:23:06 nbsp Went by today during business hours, and it was closed. Is this place still in business? Users/jsg718
  I read in the Enterprise that they were closing temporarily for the winter, but that they would be reopening very soon with a slightly more focused business.  I think they really want to emphasize their specialty flavors, which makes me very, very happy.  Their blueberry ice cream is out of this world  keep in mind that this remark comes from someone who usually doesnt bother with dessert that isnt chocolate.   Users/CovertProfessor

20090211 15:26:40 nbsp are they opening a cupcake place in this spot? Users/Yesenia
  I read that the Creamery will offer cupcakes in addition to the ice cream.   Users/CovertProfessor

20090211 20:14:17 nbsp They have reopened and showcase off a nice variety of the cupcakes by Cupcakes Cravings.  All looked incredibly yummy and it was hard to make a choice.  Ended up with a Koo Koo for Macaroons and a Dulce de leche and they were both really good.  The Dulce was super yummy with very silky frosting.  You can see the cupcake flavors at http://www.cupcakecraving.com.  Nice edition to the Creamery and a very satisfying desert option. Users/loneshark

20090212 17:19:15 nbsp The ice cream is great, but a bit pricey. Now, for 1 scoop and a mixin, my total came out to $4.04!!! The red velvet cupcake isnt very good, either. The banana ice cream is amazing.  Users/TY3

20090214 17:21:09 nbsp DAVIS CREAMERY IS NOW REOPENED!

We closed the store for a few weeks for renovation, and to bring in our new partner, Cupcake Cravings. They are a privately owned and operated business out of Sacramento, and we are pleased to have their handmade and handdecorated cupcakes to offer along with our ice cream. The ice cream is the same and prepared fresh here in the shop, with a few new flavor offerings, such as the salted caramel and 100% organic flavors. We have changed sizing and pricing to hopefully make it simpler for both our customers and ourselves, and are open to any feedback. The prices have raised on both the banana and honey ice cream, since they are now 100% organic, but we have simplified everything to make it a better experience overall for every customer that walks through our door.

Thanks to all of you for your comments, keep them coming! We rely on feedback to continue changing and improving to offer you the best service and experience possible every time. Be sure to come and check out everything new we have to offer you, we feel that both our continued products and our new additions are delicious and we will do everything possible to make sure you are not disappointed!    Users/hjermstad

20090221 13:39:54 nbsp What are the hours? Im about to go bike all they way over there to check it out but I dont know if its open. Users/KellyM

20090221 22:39:02 nbsp Theyre open until 11 PM on Saturday Users/OrofinJackson

20090314 15:58:32 nbsp Had to post after I just had an incredibly disappointing experience with this place... Went in wanting a cowpie and first off the cookies they used were out of a plastic package (like the big stale inedible ones you buy at Safeway near the bakery section).  Maybe Im imagining this but in my memory these cowpies used to be made with fresh baked cookies and Im pretty sure the cookies were warmed up or something when I got my cowpie.  So this kid sticks these sick cookies into a cup after informing me they only have chocolate chip, then puts literally the smallest scoop of ice cream Ive ever gotten in between them.  I was too stunned to even say anything... I will say the person I was with got about 2.5 times as much ice cream I did, so either the scooper hated me or was just being really lazy.  I cant believe I just paid $3.50 for a thimbleful of vanilla ice cream and cookies that I wouldnt even feed to a pigeon.  :(  I will definitely never return to this place.   Users/soledad101

20090318 11:52:26 nbsp really cute locally owned place.  love that they make their ice cream instore.  pints and quarts are a little pricey, and some of the cupcakes dont really taste like a huge step up from safeway cupcakes that you could get across the parking lot, but some of them are awesome.  bavarian cream and red velvet rock!  they have giant snickerdoodle cookies that are REALLY good here... soft and chewy.  the ice cream is great too.  its so hard to pick because they all taste so good. theres cute little asian girl that works there knows her stuff and has really good recommendations when you dont know what to pick.  service was great.  some asshole in front of us was giving her a really hard time but the whole situation was handled really professionally and i was actually surprised. usually places like this have barelyoldenough punks who look like theyd rather be somewhere else. Users/houndstooth
   Out of curiosity, how can you give somebody a hard time over ice cream? jw

20090325 22:09:17 nbsp I finally made it back to check out this place after the reopening.  Some comments, in no particular order, partly in response to comments above:  We also had the experience of very different portion sizes for the same price, from the same server, one which was very small and one which was reasonable (for the price).  Very odd.  The cookies we saw were definitely homemade; perhaps they run out occasionally?  I am consistently drawn to their fruitflavored ice cream, in spite of the fact that I dont generally like fruitflavored ice cream.  It is that good.  (I think I had raspberry this time).  We split a cupcake, and it was yummy  they have some sort of one scoop/one cupcake deal (sorry, I dont remember the details).  Last but not least:  Please bring back brownie as a mixin option!!   Users/CovertProfessor

20090328 17:12:37 nbsp third time to try cupcakes and my partner and I were both very bummed that the cake wasnt moist and seemed a bit stale or air exposed.  Case full of great options but doesnt mean much if they arent really fresh.  Made my own brownies today. Users/loneshark

20090425 17:22:48 nbsp Always been very pleased with the ice cream/milk shakes!  Since the remodel, the wait is much shorter when theres a line and the addition of cupcakes is great (although some of the flavors are more consistent in quality than others).  I highly recommend it! Users/cfb

20090519 08:08:08 nbsp Bought a few cupcakes here for the first time last week.  Very disappointing, very cheap Safewaylike taste, especially the frosting.  The cake didnt taste much like chocolate, either.  I couldnt even finish my chocolate cupcakea first ever for me.  Try the cupcakes at the Nugget in North Davis instead. Users/LeslieMadsenBrooks

20090601 14:56:49 nbsp I love this place, and the staff is so friendly. Now that the weather is getting hotter I know Ill be trying more of their flavors.  Users/nkristis

20090601 22:53:41 nbsp YAY for dollar Mondays! I went there tonight for the first time, and got a scoop of peanut butter icecream mixed with chocolate ice cream. They were both delicious, and I finished my meal wanting more. My friend had the pomegranate raspberry and enjoyed it. Im going to try to chai next time I go in.

And their cupcakes looks absolutely divine, but I havent tried them yet. My friend said they looked better than they tasted, but Ill let you know when I try one myself.
 Users/BrookeB

20090706 20:26:35 nbsp I love Davis Creamery!!! They have unique flavors and I love how its so creamy... most of all, Im so thankful for dollar scoop Mondays! For once I have a reason to not hate Mondays! Users/JanetHung

20090831 16:11:29 nbsp Best ice cream in Davis. Especially their blueberry, raspberry and banana is extremely good.
 Users/Amanda09

20090914 15:19:27 nbsp Weekly specials are the best  find myself there on Mondays, Tuesdays and Wednesdays. See above for the specials. Users/DerrickUm

20090917 01:57:26 nbsp I went here on Monday and it was really good. I plan to come back next Monday, and maybe get a cupcake. Users/NickSchmalenberger
 Im with Nick. The sorbet I tried (mango) was awesome. For a buck a scoop, you really cant go wrong. ) wl

20090918 19:16:32 nbsp Theyve killed Marble Slab!  We will never be going here again. $4.75 for a Dixie cup full of ice cream (and they call that two scoops?) with one mixin.  They were out of the mixin my wife wanted and also out of waffle cones (how does an ice cream place run out of cones?)  I miss Marble Slab.  Users/Chuckles

20091017 21:27:08 nbsp I have an idea for you: Ice cream cupcakes. They would be more gourmet than cake cupcakes, could be sold for a reasonable price, etc. Users/SteveDavison

20091111 11:21:44 nbsp I love this place. It is wellpriced, the customer service is always great, and their hours are convenient. They often have a couple of unique and interesting flavors of ice cream that I try when I am feeling adventurous. Yay for a local business! Users/GJC02

20091111 20:57:12 nbsp I absolutely love the ice cream and cupcakes, etc. Such a pleasant staff and really local feel. The daily specials are fantastic. I love the ice cream sandwich, especially with oatmeal raisin cookies. Also, I like the music they play. I will keep coming back. Users/GarySaylin

20091118 23:00:58 nbsp The ice cream was really really good.  The texture is appealing and consistent while the ice cream flavors are true; very bright and intense without any hints of manufactured chemicals.  If you like pumpkin ice cream, Davis Creamerys version is sublime.   Users/ZachStednick

20091229 13:39:04 nbsp I dont often go to South Davis, but this place is worth the trip! The salted caramel ice cream is amazing! Users/EmilyHughes

20100120 12:00:28 nbsp I love this place. The daily deals are amazing. The cupcakes are super tasty and the ice cream is delicious. I had a cowpie last night and it was good, but i recommend going early so you can get a fresh cookie. They had run out of them, so we picked from the dayold cookies. Even then it was tasty. I think everything here is good and for the most part, the staff has been friendly and helpful. Ill definitely be back for the deals.  Users/emmaleanne

20100128 17:50:28 nbsp My favorite ice cream in town and the best peppermint ice cream Ive ever had!  Also love the cowpies and their nightly deals, however, the cupcakes leave something to be desired... Users/5thAve

20100203 14:42:21 nbsp Recently tried Davis Creamery for myself  I love it! Ive gone twice in the last 2 weeks. The employees are really nice and helpful, which is great especially because I have a terrible time deciding which flavors I want to try. The cupcakes are alright, but something about them seems a little off to me. Really love the daily specials! Users/AliceW

20100205 14:56:06 nbsp The ice cream here is very, very good.  It has a smooth creamy homemade texture.  The vanilla brown sugar is a nice unique twist on a classic flavor.  The cookies are also real good.  I dont care much for the cupcakes but I dont like cake much at all.  The caramelo cookies and smores cookies are awesome.  The best cowpies are the caramelo with the vanilla brown sugar.  The best deal is by far Tuesday night $2 cowpies it is such a great price for the quality that you get.  I often feel guilty only paying $2 for those things. As of June 2010 Tuesday deal is $2.50 cowpies, still a great bargain. Users/DagonJones

20100222 16:07:56 nbsp The ice cream is very good, and homemade. The yogurt is also very good. With this said it is very expensive.  Users/GiulianaWasHere

20100318 18:28:56 nbsp this place is BOMB! such good ice cream, tons of flavors, really cheap and open pretty late. No complaints! Users/cj2012aggie

20100714 22:06:25 nbsp Crackle  enough said. Users/BamCam

20100816 00:18:08 nbsp No more $1 scoop nite. (
Its $1.50 now, which I guess is still good for the taste/quality, but still....damn. Users/TheShah
  Also, their new Pineapple Blast is awesome! Users/TheShah SS
  20100818 00:23 Had a cow pie today, went up to $3. I had a macadamia nut cookie and chocolate chip/M&M cookie, and they werent that great. Ill stick to the ice cream and cupcakes, I guess. SS

20100819 14:37:18 nbsp Ive had their avocado coconut ice cream at the farmers market and it is delicious.  This is my favorite place for ice cream in Davis.  I went into the shop today and sampled the bubblegum and salted caramel flavors, and they were both tasty.  I got a sugarfree chocolate cupcake and it was very meh, but my guess is the regular cupcakes will better, so Ill be back another time to try one of those.  Service was friendly and polite. Users/DanaSullivan

20110331 15:49:01 nbsp Drove by this morning and the place is GONE!  Sign has been removed and the store has been cleared out.  The ice cream was good the cupcakes not so much.   Users/loneshark
  As noted above, they were supposed to move to 113 E St. in The Lofts in May, so maybe theyre not really gone.   Users/CovertProfessor
  Yes, they are on E St. now, two doors north of the John Natsoulas Gallery (on the same side of the street). Users/eganio

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
cdis
cdis    Open Source License/Disclaimer, Forecast Systems Laboratory
cdis    NOAA/OAR/FSL, 325 Broadway Boulder, CO 80305
cdis
cdis    This software is distributed under the Open Source Definition,
cdis    which may be found at http://www.opensource.org/osd.html.
cdis
cdis    In particular, redistribution and use in source and binary forms,
cdis    with or without modification, are permitted provided that the
cdis    following conditions are met:
cdis
cdis    - Redistributions of source code must retain this notice, this
cdis    list of conditions and the following disclaimer.
cdis
cdis    - Redistributions in binary form must provide access to this
cdis    notice, this list of conditions and the following disclaimer, and
cdis    the underlying source code.
cdis
cdis    - All modifications to this software must be clearly documented,
cdis    and are solely the responsibility of the agent making the
cdis    modifications.
cdis
cdis    - If significant modifications or enhancements are made to this
cdis    software, the FSL Software Policy Manager
cdis    (softwaremgr@fsl.noaa.gov) should be notified.
cdis
cdis    THIS SOFTWARE AND ITS DOCUMENTATION ARE IN THE PUBLIC DOMAIN
cdis    AND ARE FURNISHED "AS IS."  THE AUTHORS, THE UNITED STATES
cdis    GOVERNMENT, ITS INSTRUMENTALITIES, OFFICERS, EMPLOYEES, AND
cdis    AGENTS MAKE NO WARRANTY, EXPRESS OR IMPLIED, AS TO THE USEFULNESS
cdis    OF THE SOFTWARE AND DOCUMENTATION FOR ANY PURPOSE.  THEY ASSUME
cdis    NO RESPONSIBILITY (1) FOR THE USE OF THE SOFTWARE AND
cdis    DOCUMENTATION; OR (2) TO PROVIDE TECHNICAL SUPPORT TO USERS.
cdis
cdis
cdis
cdis
cdis

        subroutine interp_prof(ob_pr_ht_obs,ob_pr_u_obs,ob_pr_v_obs, ! I
     1                             u_diff       , v_diff,            ! I
     1                             u_interp     , v_interp,          ! O
     1                             di_interp    , sp_interp,         ! O
     1                             i_pr,ht,level,nlevels_obs_pr,     ! I
     1                             lat_pr,lon_pr,i_ob,j_ob,          ! I
     1                             r_missing_data,                   ! I
     1                             heights_3d,ni,nj,nk,              ! I
     1                             MAX_PR,MAX_PR_LEVELS,             ! I
     1                             n_vel_grids,istatus)              ! I/O

!************************ARRAYS.FOR******************************************

!       Profiler Observations

        integer nlevels_obs_pr(MAX_PR)
        real ob_pr_ht_obs(MAX_PR,MAX_PR_LEVELS)
        real ob_pr_u_obs(MAX_PR,MAX_PR_LEVELS)
        real ob_pr_v_obs(MAX_PR,MAX_PR_LEVELS)

!**************************** END ARRAYS.FOR ********************************
        real heights_3d(ni,nj,nk)

        u_interp = r_missing_data
        v_interp = r_missing_data
        di_interp = r_missing_data
        sp_interp = r_missing_data

!  ***  Interpolate the profiler observations to the input height *******
        do i_obs = 1,nlevels_obs_pr(i_pr)

          if(i_obs .gt. 1)then

            if(ob_pr_ht_obs(i_pr,i_obs-1) .le. ht .and.
     1         ob_pr_ht_obs(i_pr,i_obs  ) .ge. ht       )then

                h_lower  = ob_pr_ht_obs(i_pr,i_obs-1)
                h_upper =  ob_pr_ht_obs(i_pr,i_obs  )

                height_diff = h_upper - h_lower

                fracl = (h_upper - ht) / height_diff
                frach = 1.0 - fracl

                u_interp = ob_pr_u_obs(i_pr,i_obs) * frach
     1                    +       ob_pr_u_obs(i_pr,i_obs-1) * fracl

                v_interp = ob_pr_v_obs(i_pr,i_obs) * frach
     1                    +       ob_pr_v_obs(i_pr,i_obs-1) * fracl

!               Correct for the time lag
!               u_diff      = du/dt * [t(anal) - t(ob)]
                u_interp = u_interp + u_diff

!               v_diff      = dv/dt * [t(anal) - t(ob)]
                v_interp = v_interp + v_diff

!               Calculate direction and speed
                call uv_to_disp( u_interp,
     1                           v_interp,
     1                           di_interp,
     1                           sp_interp)

             endif
          endif

        enddo ! level

        if(.true.)return

!       Lower Tail

        if( float(level)
     1    .lt. height_to_zcoord2(ob_pr_ht_obs(i_pr,1),heights_3d
     1                          ,ni,nj,nk,i_ob,j_ob,istatus)
     1                          .and.
     1  (height_to_zcoord2(ob_pr_ht_obs(i_pr,1),heights_3d
     1                          ,ni,nj,nk,i_ob,j_ob,istatus)
     1       - float(level)) .le. 0.5)then

                u_interp  = ob_pr_u_obs(i_pr,1)
                v_interp  = ob_pr_v_obs(i_pr,1)

!               Correct for the time lag
                u_interp = u_interp + u_diff
                v_interp = v_interp + v_diff

!               Calculate direction and speed
                call uv_to_disp( u_interp,
     1                           v_interp,
     1                           di_interp,
     1                           sp_interp)

        endif

        if(istatus .ne. 1)then
            write(6,*)
     1        ' Warning, out of bounds in height_to_zcoord2/interp_prof'
        endif

!       Upper Tail

        if( float(level) .gt.
     1  height_to_zcoord2(ob_pr_ht_obs(i_pr,nlevels_obs_pr(i_pr))
     1                  ,heights_3d,ni,nj,nk,i_ob,j_ob,istatus)
     1                           .and.
     1    (float(level) -
     1     height_to_zcoord2(ob_pr_ht_obs(i_pr,nlevels_obs_pr(i_pr))
     1                  ,heights_3d,ni,nj,nk,i_ob,j_ob,istatus))
     1                                          .le. 0.5)then

                u_interp  = ob_pr_u_obs(i_pr,nlevels_obs_pr(i_pr))
                v_interp  = ob_pr_v_obs(i_pr,nlevels_obs_pr(i_pr))

!               Correct for the time lag
                u_interp = u_interp + u_diff
                v_interp = v_interp + v_diff

!               Calculate direction and speed
                call uv_to_disp( u_interp,
     1                           v_interp,
     1                           di_interp,
     1                           sp_interp)

        endif

        if(istatus .ne. 1)then
            write(6,*)' Warning, out of bounds in height_to_zcoord2/inte
     1rp_prof'
        endif

        return

        end

max_stars_repo_name: maxinye/laps-mirror

max_stars_repo_licenses: Intel Unlicense OLDAP-2.2.1 NetCDF

max_issues_repo_name: longwosion/laps-mirror

max_issues_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense

max_forks_repo_name: longwosion/laps-mirror

max_forks_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense


Module hsay_structure
  Use math_constants
  Use hs_structure
  Implicit None
  Contains

    !Function that calculates the direct correlation function of a Hard Sphere +
    !Attractive Yukawa System using Sharma-Sharma approximation Scheme and Verlet-Weiss
    !for the hard core part
    ! c(k,ϕ,T,z)
    !wave-vector "k" for "k>=0"
    !density "ϕ" for "ϕ>0"
    !temperature "T" for "T>0"
    !Yukawa potential length parameter "z>0"
    Function c_hsay_shvw(phi,T,z,k) result(c)
      Implicit None
      Real * 8, Intent(in) :: phi !Volume fraction value
      Real * 8, Intent(in) :: T !Temperature
      Real * 8, Intent(in) :: z !Yukawa potential length parameter
      Real * 8, Intent(in) :: k !Value of the wave vector
      Real * 8 :: c
      Real * 8 :: k_func,ckdum,chs
      Real * 8 :: k2,z2,k4
      chs = c_hs_vw(phi,k)
      !chs = c_hs_py(phi,k)
      k2 = k * k
      z2 = z ** 2
      if (k > 0.075d0) then
        c = ( k * cos(k) ) + ( z * sin(k) )
        c = c / k
      else
        k4 = k2 * k2
        c = (1.d0 + z) - ((3.d0 + z) / 6.d0) * k2 + ((5.d0 + z) / 120.d0) * k4
      end if
      c = 4.d0 * pi * c / (T * (k2 + z2))
      c = chs + c
      return
    End Function

    !Function that calculates the inverse structure factor of a Hard Sphere +
    !Attractive Yukawa System using Sharma-Sharma approximation Scheme and Verlet-Weiss
    ! 1/S(k,ϕ,T,z)
    !wave-vector "k" for "k>=0"
    !density "ϕ" for "ϕ>0"
    !temperature "T" for "T>0"
    !Yukawa potential length parameter "z>0"
    Function is_hsay_shvw(phi,T,z,k) result(is)
      Implicit None
      Real * 8, Intent(in) :: phi !Volume fraction value
      Real * 8, Intent(in) :: T !Temperature
      Real * 8, Intent(in) :: z !Yukawa potential length parameter
      Real * 8, Intent (in) :: k !Value of the wave vector
      Real * 8 :: is
      Real * 8 :: c
      Real * 8 :: chs
      Real * 8 :: k2,z2,k4, phi_vw
      phi_vw = phi*(1.d0 - (phi / 16.d0))
      chs = c_hs_vw(phi,k)
      !chs = c_hs_py(phi,k)
      k2 = k * k
      z2 = z ** 2
      if (k > 0.075d0) then
        c = ( k * cos(k) ) + ( z * sin(k) )
        c = c / k
      else
        k4 = k2 * k2
        c = (1.d0 + z) - ((3.d0 + z) / 6.d0) * k2 + ((5.d0 + z) / 120.d0) * k4
      end if
      c = 4.d0 * pi * c / (T * (k2 + z2))
      is =  (phi_vw * chs) + (phi * c)
      is = 1.d0 - 6.d0 * is / pi
      return
    End Function

    !Function that calculates the structure factor of a Hard Sphere +
    !Attractive Yukawa System using Sharma-Sharma approximation Scheme and Verlet-Weiss
    ! S(k,ϕ,z)
    !wave-vector "k" for "k>=0"
    !density "ϕ" for "ϕ>0"
    !temperature "T" for "T>0"
    !Yukawa potential length parameter "z>0"
    Function s_hsay_shvw(phi,T,z,k) result(s)
      Implicit None
      Real * 8, Intent(in) :: phi !Volume fraction value
      Real * 8, Intent(in) :: T !Temperature
      Real * 8, Intent(in) :: z !Yukawa potential length parameter
      Real * 8, Intent (in) :: k !Value of the wave vector
      Real * 8 :: s
      s = 1.d0 / is_hsay_shvw(phi,T,z,k)
      return
    End Function

    !Function that calculates the structure factor of a Hard Sphere +
    !Attractive Yukawa System using Generalized Sharma-Sharma approximation Scheme and Verlet-Weiss
    ! S(k,ϕ,z)
    !wave-vector "k" for "k>=0"
    !density "ϕ" for "ϕ>0"
    !temperature "T" for "T>0"
    !Yukawa potential length parameter "z>0"
    Function c_hsay_gsvw(phi,T,z,k,qn,qw) result(c)
      Implicit None
      Real * 8, Intent(in) :: phi !Volume fraction value
      Real * 8, Intent(in) :: T !Temperature
      Real * 8, Intent(in) :: z !Yukawa potential length parameter
      Real * 8, Intent(in) :: k !Value of the wave vector
      Real * 8, Intent(in), dimension(:) :: qn  !Quadrature nodes
      Real * 8, Intent(in), dimension(:) :: qw  !Quadrature weights
      Real * 8 :: c
      Real * 8 :: cay,u,chs,r,w, dumsin,kr
      Real * 8 :: k2,z2,k4
      Integer :: i1
      cay = 0.d0
      do i1=1, size(qn)
        r = qn(i1)
        kr = k * r
        if (kr > 0.075d0 ) then
          w = qw(i1) * r * sin(k * r) / k
        else
          w = qw(i1) * r**2 * (1.d0 + ( (k**2) * r/6.d0) )
        end if
        u = - exp(z * (r-1.d0) )/r
        cay = cay + ( ( exp(-u/T) - 1.d0 ) * w )
      end do
      chs = c_hs_vw(phi,k)
      c = chs + cay
      return
    End Function

    !Function that calculates the inverse structure factor of a Hard Sphere +
    !Attractive Yukawa System using Generalized Sharma-Sharma approximation Scheme and Verlet-Weiss
    ! 1/S(k,ϕ,T,z)
    !wave-vector "k" for "k>=0"
    !density "ϕ" for "ϕ>0"
    !temperature "T" for "T>0"
    !Yukawa potential length parameter "z>0"
    Function is_hsay_gsvw(phi,T,z,k,qn,qw) result(is)
      Implicit None
      Real * 8, Intent(in) :: phi !Volume fraction value
      Real * 8, Intent(in) :: T !Temperature
      Real * 8, Intent(in) :: z !Yukawa potential length parameter
      Real * 8, Intent(in) :: k !Value of the wave vector
      Real * 8, Intent(in), dimension(:) :: qn  !Quadrature nodes
      Real * 8, Intent(in), dimension(:) :: qw  !Quadrature weights
      Real * 8 :: is
      is = c_hsay_gsvw(phi,T,z,k,qn,qw)
      is = 1.d0 - (6.d0 * phi * is / pi)
      return
    End Function

    !Function that calculates the structure factor of a Hard Sphere +
    !Attractive Yukawa System using Generalized Sharma-Sharma approximation Scheme and Verlet-Weiss
    ! S(k,ϕ,z)
    !wave-vector "k" for "k>=0"
    !density "ϕ" for "ϕ>0"
    !temperature "T" for "T>0"
    !Yukawa potential length parameter "z>0"
    Function s_hsay_gsvw(phi,T,z,k,qn,qw) result(s)
      Implicit None
      Real * 8, Intent(in) :: phi !Volume fraction value
      Real * 8, Intent(in) :: T !Temperature
      Real * 8, Intent(in) :: z !Yukawa potential length parameter
      Real * 8, Intent(in) :: k !Value of the wave vector
      Real * 8, Intent(in), dimension(:) :: qn  !Quadrature nodes
      Real * 8, Intent(in), dimension(:) :: qw  !Quadrature weights
      Real * 8 :: s
      s = 1.d0 / is_hsay_gsvw(phi,T,z,k,qn,qw)
      return
    End Function

    Function Ts_hsay_shvw(phi,z) result(Ts)
      Implicit None
      Real * 8, Intent(in) :: phi !Volume fraction value
      Real * 8 :: Ts !Temperature
      Real * 8, Intent(in) :: z !Yukawa potential length parameter
      Real * 8 :: S_HS
      S_HS = s_hs_vw(phi,0.d0)
      Ts = 24.d0 * phi * (z + 1.d0)  * S_HS / (z*z)
      return
    End Function

End Module