text:
!##############################################################################
!# Tutorial 004i: Datastructures - quadtrees
!##############################################################################

module tutorial004i

  ! Include basic Feat-2 modules
  use fsystem
  use genoutput
  use quadtreebase
  use quadtreeDP

  implicit none
  private

  public :: start_tutorial004i

contains

  ! ***************************************************************************

  subroutine start_tutorial004i

    ! Declare some variables
    type(t_quadtreeDP) :: rquadtree
    integer :: iresult

    ! Print a message
    call output_lbrk()
    call output_separator (OU_SEP_STAR)
    call output_line ("This is FEAT-2. Tutorial 004i")
    call output_separator (OU_SEP_MINUS)

    ! =================================
    ! Create quadtree with bounding box (0,1)x(0,2)
    ! which can store 20 vertices in 10 nodes before
    ! new memory is allocated
    ! =================================

    call qtree_create(rquadtree, 20, 10, 0.0_DP, 0.0_DP, 1.0_DP, 2.0_DP)

    ! =================================
    ! Insert new vertices into quadtree
    ! =================================

    iresult = qtree_insert(rquadtree, (/0.2_DP, 0.6_DP/))
    iresult = qtree_insert(rquadtree, (/0.2_DP, 0.5_DP/))
    iresult = qtree_insert(rquadtree, (/0.3_DP, 0.6_DP/))
    iresult = qtree_insert(rquadtree, (/0.5_DP, 0.6_DP/))
    iresult = qtree_insert(rquadtree, (/0.6_DP, 0.7_DP/))
    iresult = qtree_insert(rquadtree, (/0.8_DP, 2.0_DP/))
    iresult = qtree_insert(rquadtree, (/0.7_DP, 0.8_DP/))
    iresult = qtree_insert(rquadtree, (/0.4_DP, 1.6_DP/))

    select case(iresult)
    case (QTREE_FAILED)
      call output_line ("Failed to insert vertex into quadtree!")
    case (QTREE_FOUND)
      call output_line ("Vertex is already present in quadtree!")
    case (QTREE_INSERTED)
      call output_line ("Inserted vertex into quadtree!")
    end select

    iresult = qtree_insert(rquadtree, (/0.4_DP, 1.6_DP/))

    select case(iresult)
    case (QTREE_FAILED)
      call output_line ("Failed to insert vertex into quadtree!")
    case (QTREE_FOUND)
      call output_line ("Vertex is already present in quadtree!")
    case (QTREE_INSERTED)
      call output_line ("Inserted vertex into quadtree!")
    end select

    ! =================================
    ! Delete vertex from quadtree
    ! =================================

    iresult = qtree_delete(rquadtree, (/0.2_DP, 0.6_DP/))

    select case(iresult)
    case (QTREE_FAILED)
      call output_line ("Failed to delete vertex from quadtree!")
    case (QTREE_DELETED)
      call output_line ("Deleted vertex from quadtree!")
    end select

    iresult = qtree_delete(rquadtree, (/0.1_DP, 1.3_DP/))

    select case(iresult)
    case (QTREE_FAILED)
      call output_line ("Failed to delete vertex from quadtree!")
    case (QTREE_DELETED)
      call output_line ("Deleted vertex from quadtree!")
    end select

    ! =================================
    ! Find vertex in quadtree
    ! =================================

    iresult = qtree_find(rquadtree, (/0.501_DP, 0.601_DP/))

    select case(iresult)
    case (QTREE_NOT_FOUND)
      call output_line ("Unable to find vertex in quadtree!")
    case (QTREE_FOUND)
      call output_line ("Found vertex in quadtree!")
    end select

    ! =================================
    ! Find vertex in quadtree using user-defined isEqual
    ! function (with less restrictive proximity definition)
    ! =================================

    iresult = qtree_find(rquadtree, (/0.501_DP, 0.601_DP/),&
                         fcb_isEqual=my_isEqual)

    select case(iresult)
    case (QTREE_NOT_FOUND)
      call output_line ("Unable to find vertex in quadtree!")
    case (QTREE_FOUND)
      call output_line ("Found vertex in quadtree!")
    end select

    ! =================================
    ! Release quadtree
    ! =================================

    call qtree_release(rquadtree)

  end subroutine

  ! ***************************************************************************

  pure logical function my_isEqual(data1,data2) result(bisEqual)
    real(DP), dimension(2), intent(in) :: data1,data2

    bisEqual = (maxval(abs(data1-data2)) .lt. 1e-1)

  end function my_isEqual

end module

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


text:
      subroutine phdrw(ix,np,word,title,case,exist,lfno)
      use tfstk
      use tmacro
      implicit real*8(a-h,o-z)
      integer*4 ax(2)
      character*(*) word,title,case
      character*8 coord(6)
      logical*4 abbrev,exist,fconv,lod,append
      real*8 s(8)
      data coord  /'X       ','PX      ','Y       ',
     z             'PY      ','Z       ','PZ      '/
      fconv=.false.
      lfn=lfno
c .... parsing word
      if(.not.exist) goto 2
    1 call getwdl(word)
    2 continue
      ax(1)=0
      ax(2)=0
      do 11 j=1,3
        if(word(j:j).eq.'-') then
          do 12 k=1,6
            if(word(:j-1).eq.coord(k)(:j-1)) then
              ax(1)=k
            endif
            if(word(j+1:j+2).eq.coord(k)(:2)) then
              ax(2)=k
            endif
   12     continue
          goto 19
        endif
   11 continue
      do 10 i=1,6
        if(word.eq.coord(i)) then
          ax(1)=i
          ax(2)=i+2*(i-2*(i/2))-1
          goto 19
        endif
   10 continue
      if(word.eq.'OUT'.or.abbrev(word,'APP_END','_')) then
        append=abbrev(word,'APP_END','_')
        lfn=int(getva(exist))
        if(exist)then
          if(lfn .le. 0)then
            lfn=6
          endif
          goto 1
        else
          lfn=98
          inquire(unit=lfn,opened=lod,iostat=ios)
          if(lod) close(lfn)
          call getwrd(word)
          call texpfn(word)
          if(append) then
            open(lfn,file=word,status='UNKNOWN',
     1           access='APPEND',ERR=6101)
          else
            open(lfn,file=word,status='UNKNOWN',
     1           ERR=6101)
          endif
          goto 1
        endif
6101    lfn=6
        call termes(lfno,'?File open error ',word)
        goto 999
      endif
   19 continue
c
      if(ax(1)*ax(2).eq.0) then
        exist=.false.
        goto 999
      endif
c .... convert to canonical variables
      if(.not.fconv) then
        ix1=italoc(np0*8)
        fconv=.true.
      endif
      do 100 i=1,np
        do 101 j=1,8
  101     s(j)=rlist(ix+np0*(j-1)+i-1)
        call tconv(s,s,1)
        do 102 j=1,8
  102     rlist(ix1+np0*(j-1)+i-1)=s(j)
  100 continue
      call phdrwa(rlist(ix1+np0*(ax(2)-1)),rlist(ix1+np0*(ax(1)-1)),
     z            np,ax,title,case,lfn)
      goto 1
c
  999 continue
      if(fconv) then
        call tfree(int8(ix1))
      endif
      if(lfn .ne. lfno) close(lfn)
      return
      end

ext: f

max_stars_repo_name: noboruatkek/SAD

max_issues_repo_name: noboruatkek/SAD

max_forks_repo_name: noboruatkek/SAD


      subroutine ch_rtpest

!!     ~ ~ ~ PURPOSE ~ ~ ~
!!     this subroutine computes the daily stream pesticide balance
!!     (soluble and sorbed)

!!    ~ ~ ~ INCOMING VARIABLES ~ ~ ~
!!    name          |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    chpst_conc(:) |mg/(m**3)     |initial pesticide concentration in reach
!!    chpst_koc(:)  |m**3/g        |pesticide partition coefficient between
!!                                 |water and sediment in reach
!!    chpst_mix(:)  |m/day         |mixing velocity (diffusion/dispersion) for
!!                                 |pesticide in reach
!!    chpst_rea(:)  |1/day         |pesticide reaction coefficient in reach
!!    chpst_rsp(:)  |m/day         |resuspension velocity in reach for pesticide
!!                                 |sorbed to sediment
!!    chpst_stl(:)  |m/day         |settling velocity in reach for pesticide
!!                                 |sorbed to sediment
!!    chpst_vol(:)  |m/day         |pesticide volatilization coefficient in
!!                                 |reach
!!    drift(:)      |kg            |amount of pesticide drifting onto main
!!                                 |channel in subbasin
!!    rchdep        |m             |depth of flow on day
!!    rchwtr        |m^3 H2O       |water stored in reach at beginning of day
!!    sedpst_rea(:) |1/day         |pesticide reaction coefficient in river bed
!!                                 |sediment
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    ~ ~ ~ OUTGOING VARIABLES ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    bury        |mg pst        |loss of pesticide from active sediment layer
!!                               |by burial
!!    difus       |mg pst        |diffusion of pesticide from sediment to reach
!!    reactb      |mg pst        |amount of pesticide in sediment that is lost
!!                               |through reactions
!!    reactw      |mg pst        |amount of pesticide in reach that is lost
!!                               |through reactions
!!    resuspst    |mg pst        |amount of pesticide moving from sediment to
!!                               |reach due to resuspension
!!    setlpst     |mg pst        |amount of pesticide moving from water to
!!                               |sediment due to settling
!!    solpesto    |mg pst/m^3    |soluble pesticide concentration in outflow
!!                               |on day
!!    sorpesto    |mg pst/m^3    |sorbed pesticide concentration in outflow
!!                               |on day
!!    volatpst    |mg pst        |amount of pesticide in reach lost by
!!                               |volatilization
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ LOCAL DEFINITIONS ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    chpstmass   |mg pst        |mass of pesticide in reach
!!    depth       |m             |depth of water in reach
!!    fd2         |
!!    frsol       |none          |fraction of pesticide in reach that is soluble
!!    frsrb       |none          |fraction of pesticide in reach that is sorbed
!!    jrch        |none          |reach number
!!    pstin       |mg pst        |total pesticide transported into reach
!!                               |during time step
!!    sedcon      |g/m^3         |sediment concentration
!!    sedpstmass  |mg pst        |mass of pesticide in bed sediment
!!    solpstin    |mg pst        |soluble pesticide entering reach during
!!                               |time step
!!    sorpstin    |mg pst        |sorbed pesticide entering reach during
!!                               |time step
!!    tday        |days          |flow duration
!!    wtrin       |m^3 H2O       |volume of water entering reach during time
!!                               |step
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ SUBROUTINES/FUNCTIONS CALLED ~ ~ ~
!!    Intrinsic: Abs

!!    ~ ~ ~ ~ ~ ~ END SPECIFICATIONS ~ ~ ~ ~ ~ ~

      use channel_data_module
      use channel_module
      use sd_channel_module
      use ch_pesticide_module
      use hydrograph_module, only : ob, jrch, ht1, ch_stor
      use constituent_mass_module
      use pesticide_data_module

      implicit none

      integer :: ipest          !none          |pesticide counter - sequential
      integer :: jpst           !none          |pesticide counter from data base
      real :: pstin             !mg pst        |total pesticide transported into reach during time step
      real :: kd                !(mg/kg)/(mg/L) |koc * carbon
      real :: depth             !m             |depth of water in reach
      real :: chpstmass         !mg pst        |mass of pesticide in reach
      real :: sedpstmass        !mg pst        |mass of pesticide in bed sediment
      real :: fd2               !units         |description
      real :: solmax            !units         |description
      real :: sedcon            !g/m^3         |sediment concentration
      real :: tday              !none          |flow duration (fraction of 24 hr)
      real :: rchwtr            !m^3 H2O       |water stored in reach at beginning of day
      real :: por               !none          |porosity of bottom sediments
      real :: pest_init         !mg            |amount of pesticide before decay
      real :: pest_end          !mg            |amount of pesticide after decay

      !! zero outputs
      chpst_d(jrch) = chpstz

      !! initialize depth of water for pesticide calculations
      if (rchdep < 0.01) then
        depth = .01
      else
        depth = rchdep
      endif

      do ipest = 1, cs_db%num_pests
        jpst = cs_db%pest_num(ipest)

        !! volume of water entering reach and stored in reach
        wtrin = ht1%flo + ch_stor(jrch)%flo

        !! pesticide transported into reach during day
        pstin = hcs1%pest(ipest)

        !! calculate mass of pesticide in reach
        chpstmass = pstin + ch_water(jrch)%pest(ipest)

        !! calculate mass of pesticide in bed sediment
        sedpstmass = ch_benthic(jrch)%pest(ipest)

        if (chpstmass + sedpstmass < 1.e-12) then
          ch_water(jrch)%pest(ipest) = 0.
          ch_benthic(jrch)%pest(ipest) = 0.
        end if
        if (chpstmass + sedpstmass < 1.e-12) return

        !!in-stream processes
        if (wtrin / 86400. > 1.e-9) then
          !! calculate sediment concentration
          sedcon = ht1%sed / wtrin * 1.e6

          !! set kd
          kd = pestdb(jpst)%koc * sd_chd(jrch)%carbon / 100.

          !! calculate fraction of soluble and sorbed pesticide
          frsol = 1. / (1. + kd * sedcon)
          frsrb = 1. - frsol

          !! ASSUME DENSITY=2.6E6; KD2=KD1
          por = 1. - sd_chd(jrch)%bd / 2.65
          fd2 = 1. / (por + kd)

          !! calculate flow duration
          tday = rttime / 24.0
          if (tday > 1.0) tday = 1.0
          tday = 1.0

          !! calculate amount of pesticide that undergoes chemical or biological degradation on day in reach
          pest_init = chpstmass
          if (pest_init > 1.e-12) then
            pest_end = chpstmass * pestcp(jpst)%decay_a
            chpstmass = pest_end
            chpst%pest(ipest)%react = pest_init - pest_end
          end if

          !! calculate amount of pesticide that volatilizes from reach
          chpst%pest(ipest)%volat = pestdb(jpst)%aq_volat * frsol * chpstmass * tday / depth
          if (chpst%pest(ipest)%volat > frsol * chpstmass) then
            chpst%pest(ipest)%volat = frsol * chpstmass
            chpstmass = chpstmass - chpst%pest(ipest)%volat
          else
            chpstmass = chpstmass - chpst%pest(ipest)%volat
          end if

          !! calculate amount of pesticide removed from reach by settling
          chpst%pest(ipest)%settle = pestdb(jpst)%aq_settle * frsrb * chpstmass * tday / depth
          if (chpst%pest(ipest)%settle >  frsrb * chpstmass) then
            chpst%pest(ipest)%settle = frsrb * chpstmass
            chpstmass = chpstmass - chpst%pest(ipest)%settle
          else
            chpstmass = chpstmass - chpst%pest(ipest)%settle
          end if
          sedpstmass = sedpstmass + chpst%pest(ipest)%settle

          !! calculate resuspension of pesticide in reach
          chpst%pest(ipest)%resus = pestdb(jpst)%aq_resus * sedpstmass * tday / depth
          if (chpst%pest(ipest)%resus > sedpstmass) then
            chpst%pest(ipest)%resus = sedpstmass
            sedpstmass = 0.
          else
            sedpstmass = sedpstmass - chpst%pest(ipest)%resus
          end if
          chpstmass = chpstmass + chpst%pest(ipest)%resus

          !! calculate diffusion of pesticide between reach and sediment
          chpst%pest(ipest)%difus = sd_ch(jrch)%aq_mix(ipest) * (fd2 * sedpstmass - frsol * chpstmass) * tday / depth
          if (chpst%pest(ipest)%difus > 0.) then
            if (chpst%pest(ipest)%difus > sedpstmass) then
              chpst%pest(ipest)%difus = sedpstmass
              sedpstmass = 0.
            else
              sedpstmass = sedpstmass - Abs(chpst%pest(ipest)%difus)
            end if
            chpstmass = chpstmass + Abs(chpst%pest(ipest)%difus)
          else
            if (Abs(chpst%pest(ipest)%difus) > chpstmass) then
              chpst%pest(ipest)%difus = -chpstmass
              chpstmass = 0.
            else
              chpstmass = chpstmass - Abs(chpst%pest(ipest)%difus)
            end if
            sedpstmass = sedpstmass + Abs(chpst%pest(ipest)%difus)
          end if

          !! calculate removal of pesticide from active sediment layer by burial
          chpst%pest(ipest)%bury = pestdb(jpst)%ben_bury * sedpstmass / pestdb(jpst)%ben_act_dep
          if (chpst%pest(ipest)%bury > sedpstmass) then
            chpst%pest(ipest)%bury = sedpstmass
            sedpstmass = 0.
          else
            sedpstmass = sedpstmass - chpst%pest(ipest)%bury
          end if

          !! verify that water concentration is at or below solubility
          solmax = pestdb(jpst)%solub * wtrin
          if (solmax < chpstmass * frsol) then
            sedpstmass = sedpstmass + (chpstmass * frsol - solmax)
            chpstmass = chpstmass - (chpstmass * frsol - solmax)
          end if

        else
          !!insignificant flow
          sedpstmass = sedpstmass + chpstmass
          chpstmass = 0.
        end if

        !! benthic processes
        !! calculate loss of pesticide from bed sediments by reaction
        pest_init = sedpstmass
        if (pest_init > 1.e-12) then
          pest_end = sedpstmass * pestcp(jpst)%decay_b
          sedpstmass = pest_end
          chpst%pest(ipest)%react_bot = pest_init - pest_end
        end if

        !! set new pesticide mass of (in + store) after processes
        if (wtrin > 1.e-6) then
          hcs1%pest(ipest) = chpstmass
        else
          sedpstmass = sedpstmass + chpstmass
        end if
        ch_benthic(jrch)%pest(ipest) = sedpstmass

      end do

      return
      end subroutine ch_rtpest


text:
!**************************************************************
!* AceGen    7.205 MacOSX (15 Jan 21)                         *
!*           Co. J. Korelc  2020           11 Jul 21 20:29:04 *
!**************************************************************
! User     : Full professional version
! Notebook : linelasticDynamic
! Evaluation time                 : 9 s     Mode  : Optimal
! Number of formulae              : 265     Method: Automatic
! Subroutine                      : elmt16_ISW01 size: 57
! Subroutine                      : elmt16_ISW03 size: 3530
! Subroutine                      : elmt16_ISW05 size: 893
! Subroutine                      : elmt16_ISW06 size: 1065
! Subroutine                      : elmt16_ISW09 size: 450
! Total size of Mathematica  code : 5995 subexpressions
! Total size of Fortran code      : 19903 bytes

      subroutine elmt16 (d,ul,xl,ix,tl,s,p,ndf,ndm,nst,isw)

!      * * F E A P * * A Finite Element Analysis Program

!....  Copyright (c) 1984-2021: Regents of the University of California
!                               All rights reserved

!-----[--.----+----.----+----.-----------------------------------------]
!     Modification log                                Date (dd/mm/year)
!-----[--.----+----.----+----.-----------------------------------------]
!     Purpose:

!     Inputs:
!       d(*)          - Material set parameters
!       ul(ndf,nen,*) - Solution parameters
!       xl(ndm,nen)   - Element nodal coordinates
!       ix(nen1)      - Element global node numbers
!       tl(nen)       - Element vector (e.g., nodal temperatures)
!       ndf           - Maximum no dof's/node
!       ndm           - Mesh spatial dimension
!       nst           - Element matrix/vector dimension
!       isw           - Switch parameter (set by input commands)

!     Outputs:
!       s(nst,nst,2)  - Element matrix (stiffness, mass, etc.)
!       p(nst,2)      - Element vector (residual, lump mass, etc.)
!-----[--.----+----.----+----.-----------------------------------------]

      implicit none

      include 'debugs.h'
      include 'bdata.h'
      include 'cdata.h'
      include 'eldata.h'
      include 'eltran.h'
      include 'hdata.h'
      include 'iofile.h'
      include 'umac1.h'                 ! utx(1)
      include 'comblk.h'

      integer       ::  ndf             ! Max DOF's/node
      integer       ::  ndm             ! Mesh spatial dimention
      integer       ::  nst             ! Matrix 's' dimension
      integer       ::  isw             ! Switch for solution steps
      integer       ::  ix(*)           ! Element global nodes
      real (kind=8) ::  d(*)            ! Material set parameters
      real (kind=8) ::  ul(ndf,nen,*)   ! Local nodal solutions
      real (kind=8) ::  xl(ndm,*)       ! Local nodal coordinates
      real (kind=8) ::  tl(*)           ! Local nodal load array
      real (kind=8) ::  s(nst,nst,*)    ! Element matrix
      real (kind=8) ::  p(ndf,*)        ! Element vector
      logical       ::  symmetric       ! Is matrix symmetric
      logical       ::  errck           ! For error in inputs

      character (len=50) :: datades(3)
      logical       ::  pinput
      integer       ::  dofacegen,ngpo
      real (kind=8) ::  td(10)
      real (kind=8) ::  v(1125)   ! AceGen Storage
      real (kind=8), allocatable :: ua(:,:)    ! AceGen Storage
      real (kind=8), allocatable :: k(:,:)     ! AceGen stiffness
      real (kind=8), allocatable :: c(:,:)     ! AceGen damping
      real (kind=8), allocatable :: m(:,:)     ! AceGen mass
      real (kind=8), allocatable :: r(:)       ! AceGen residual

      integer       :: i,j,kk,n                 ! Loop variables
      real (kind=8) :: gp(4,25)

!     Limits: 10 PDE; 30 DOF total
      integer       :: npde, nelu(10),du(10),dofu(10,30)
      integer       :: hist1, hist3

      save

      symmetric=.false.

      if(isw.ne.1) then

!       Read number of PDEs
        npde = nint(d(100))

      endif

      if(isw.lt.0) then
        utx(1) = "elmt16"

      elseif(isw.eq.1) then                ! Input material set data
        pstyp = ndm                        ! Sets plot dimension (1,2,3)

!       Read number pdes
        nelu(:)   = 0
        du(:)     = 0
        dofu(:,:) = 0
        call elmt16_ISW01(v,npde,du,nelu,hist1,hist3)

!       Build 'dofu'
        kk = 0
        do i = 1,npde
          do j = 1, du(i)
            kk = kk + 1
            dofu(i,j) = kk
          end do ! j
        end do ! i

!       Initialize number of history variables
        nh1 = hist1
        nh3 = hist3

!       Save number of pde
        d(100) = dble(npde)

!       Build size of AceGen arrays
        dofacegen = 0
        do n = 1,npde
          dofacegen = dofacegen + nelu(n)*du(n)
        end do ! n

        if( .not.allocated(r) ) then
          if(debug) write(*,*) ' --> Allocate: DOFACEGEN =',dofacegen
          allocate( r(dofacegen) )
          allocate( m(dofacegen,dofacegen) )
          allocate( c(dofacegen,dofacegen) )
          allocate( k(dofacegen,dofacegen) )
          allocate( ua(6,dofacegen) )
        endif

!       Material properties
        errck = pinput(td,3)
        d(1:3) = td(1:3)

!       Output material set data

        datades(1)="Em-Youngs modulus"

        datades(2)="[Nu]-Poisson ratio"

        datades(3)="[Rho]0-Density"

        write(iow,"(10x,f15.5,A3,A)")
     #     (d(i)," = ",datades(i),i=1,3)

!       Set active dof's
        call SB_activedofs(du,dofu,nelu,npde,ix,ndf)

!       Set quadrature points
        nel = maxval(nelu(1:npde))
        call SB_Int(ndm, nel, ngpo,gp)

      elseif(isw.eq.3 .or.isw.eq.6) then   ! Compute residual/tangent

!       Map dofs
        call SB_ua_set(ul,ndf,nen,ua,du,dofu,nelu,npde)

        if(isw.eq.3) then

!         Initialize element arrays
          k(:,:) = 0.0d0
          c(:,:) = 0.0d0
          m(:,:) = 0.0d0
          r(:)   = 0.0d0

!         Call the routine for ISW=3
          call elmt16_ISW03(v,d,xl,ua,k,c,m,r,hr(nh1),hr(nh2),gp,ngpo)

!         Print values for debug mode
          if(debug) then
            call mprint(k,dofacegen,dofacegen,dofacegen,'ACE_K')
            call mprint(c,dofacegen,dofacegen,dofacegen,'ACE_C')
            call mprint(m,dofacegen,dofacegen,dofacegen,'ACE_M')
            call mprint(r,1,dofacegen,1,'ACE_R')
          endif

!         Combine tangent
          k(:,:) = k(:,:)*ctan(1) + c(:,:)*ctan(2) + m(:,:)*ctan(3)

!         Fill symmetric part
          if(symmetric) then
            do j = 1,dofacegen
              do i = j+1,dofacegen
                k(i,j) = k(j,i)
              end do ! i
            end do ! j
          end if

!         Map to FEAP storage
          call SB_ace2nst(k,r,dofacegen,s,p,nst,du,dofu,nelu,npde,isw)

        elseif(isw.eq.6) then

!         Initialize element arrays
          r(:)   = 0.0d0

!         Call the routine for ISW=6
          call elmt16_ISW06(v,d,xl,ua,r,hr(nh1),hr(nh2),gp,ngpo)

!         Map to FEAP storage
          call SB_ace2nst(k,r,dofacegen,s,p,nst,du,dofu,nelu,npde,isw)

        endif

      elseif(isw.eq.4 .or.isw.eq.8) then   ! Output/plot element data

      elseif(isw.eq.5) then                ! Compute mass matrix

!       Map dofs
        call SB_ua_set(ul,ndf,nen,ua,du,dofu,nelu,npde)

!       Initialize element arrays
        m(:,:) = 0.0d0
        r(:)   = 0.0d0

!       Call the routine for ISW=5
        call elmt16_ISW05(v,d,xl,ua,m,r,hr(nh1),hr(nh2),gp,ngpo)

      elseif(isw.eq.9) then                ! Compute damping matrix

!       Map dofs
        call SB_ua_set(ul,ndf,nen,ua,du,dofu,nelu,npde)

!       Initialize element arrays
        c(:,:) = 0.0d0
        r(:)   = 0.0d0

!       Call the routine for ISW=9
        call elmt16_ISW09(v,d,xl,ua,c,r,hr(nh1),hr(nh2),gp,ngpo)

      endif

      end subroutine elmt16

!******************* S U B R O U T I N E **********************
      SUBROUTINE elmt16_ISW01(v,npde,du,nelu,hist1,hist3)
      IMPLICIT NONE
      include 'sms.h'
      INTEGER npde,du(10),nelu(10),hist1,hist3
      DOUBLE PRECISION v(1125)
      npde=(1)
      du(1)=(2)
      nelu(1)=(4)
      hist1=(0)
      hist3=(0)
      END

!******************* S U B R O U T I N E **********************
      SUBROUTINE elmt16_ISW03(v,d,xl,ul,s,c,m,p,ht,hp,gp,ngpo)
      IMPLICIT NONE
      include 'sms.h'
      INTEGER ngpo,i6
      DOUBLE PRECISION v(1125),d(3),xl(2,4),ul(6,8),s(8,8),c
     &(8,8),m(8,8),p(8),ht(0),hp(0),gp(4,4)
      v(73)=ul(1,8)
      v(72)=ul(1,7)
      v(71)=ul(1,6)
      v(70)=ul(1,5)
      v(69)=ul(1,4)
      v(68)=ul(1,3)
      v(67)=ul(1,2)
      v(66)=ul(1,1)
      v(17)=xl(2,4)
      v(16)=xl(1,4)
      v(15)=xl(2,3)
      v(984)=v(15)-v(17)
      v(14)=xl(1,3)
      v(986)=v(14)-v(16)
      v(13)=xl(2,2)
      v(980)=v(13)-v(15)
      v(12)=xl(1,2)
      v(982)=v(12)-v(14)
      v(11)=xl(2,1)
      v(983)=v(11)-v(13)
      v(979)=v(11)-v(17)
      v(10)=xl(1,1)
      v(985)=v(10)-v(12)
      v(981)=v(10)-v(16)
      v(9)=d(3)
      v(8)=d(2)
      v(965)=d(1)/(1d0+v(8))
      v(164)=(v(8)*v(965))/(1d0-2d0*v(8))
      v(162)=v(965)/2d0
      v(976)=2d0*v(162)
      v(987)=v(164)+v(976)
      DO i6=1,ngpo
       v(90)=gp(1,i6)
       v(99)=1d0-v(90)
       v(105)=(-0.25d0)*v(99)
       v(97)=1d0+v(90)
       v(106)=(-0.25d0)*v(97)
       v(111)=v(105)*v(979)+v(106)*v(980)
       v(109)=v(105)*v(981)+v(106)*v(982)
       v(91)=gp(2,i6)
       v(107)=(1d0+v(91))/4d0
       v(104)=((-1d0)+v(91))/4d0
       v(110)=v(104)*v(983)+v(107)*v(984)
       v(108)=v(104)*v(985)+v(107)*v(986)
       v(94)=-(v(104)*v(99))
       v(96)=-(v(104)*v(97))
       v(98)=v(107)*v(97)
       v(101)=v(107)*v(99)
       v(112)=-(v(109)*v(110))+v(108)*v(111)
       v(969)=gp(4,i6)*v(112)
       v(966)=v(9)*v(969)
       v(968)=v(966)*v(98)
       v(967)=v(94)*v(966)
       v(325)=(v(101)*v(101))*v(966)
       v(323)=v(966)*(v(98)*v(98))
       v(321)=v(966)*(v(96)*v(96))
       v(317)=(v(94)*v(94))*v(966)
       v(977)=v(9)*(ul(5,7)*v(101)+ul(5,1)*v(94)+ul(5,3)*v(96)+ul(5,5
     & )*v(98))
       v(978)=v(9)*(ul(5,8)*v(101)+ul(5,2)*v(94)+ul(5,4)*v(96)+ul(5,6
     & )*v(98))
       v(320)=v(101)*v(967)
       v(319)=v(967)*v(98)
       v(318)=v(96)*v(967)
       v(322)=v(96)*v(968)
       v(324)=v(101)*v(968)
       v(135)=-(v(111)/v(112))
       v(149)=-(v(107)*v(135))
       v(141)=-(v(104)*v(135))
       v(136)=v(109)/v(112)
       v(152)=-(v(107)*v(136))
       v(143)=-(v(104)*v(136))
       v(137)=-(v(110)/v(112))
       v(150)=v(105)*v(137)
       v(145)=v(106)*v(137)
       v(138)=v(108)/v(112)
       v(153)=v(105)*v(138)
       v(147)=v(106)*v(138)
       v(139)=v(141)+v(150)
       v(225)=v(139)*v(162)
       v(140)=v(143)+v(153)
       v(213)=v(140)*v(164)
       v(273)=(v(139)*v(213)+v(140)*v(225))*v(969)
       v(142)=-v(141)+v(145)
       v(970)=v(142)*v(162)
       v(214)=v(142)*v(164)
       v(221)=v(214)+2d0*v(970)
       v(144)=-v(143)+v(147)
       v(971)=v(144)*v(162)
       v(215)=v(144)*v(164)
       v(288)=v(969)*(v(142)*v(215)+v(144)*v(970))
       v(275)=v(969)*(v(139)*v(215)+v(140)*v(970))
       v(233)=v(215)+2d0*v(971)
       v(282)=v(969)*(v(140)*v(233)+v(139)*v(970))
       v(281)=v(969)*(v(140)*v(214)+v(139)*v(971))
       v(274)=v(969)*(v(139)*v(221)+v(140)*v(971))
       v(146)=-v(145)+v(149)
       v(972)=v(146)*v(162)
       v(216)=v(146)*v(164)
       v(222)=v(216)+2d0*v(972)
       v(148)=-v(147)+v(152)
       v(973)=v(148)*v(162)
       v(217)=v(148)*v(164)
       v(299)=v(969)*(v(146)*v(217)+v(148)*v(972))
       v(290)=v(969)*(v(142)*v(217)+v(144)*v(972))
       v(277)=v(969)*(v(139)*v(217)+v(140)*v(972))
       v(234)=v(217)+2d0*v(973)
       v(295)=v(969)*(v(144)*v(234)+v(142)*v(972))
       v(284)=v(969)*(v(140)*v(234)+v(139)*v(972))
       v(294)=v(969)*(v(144)*v(216)+v(142)*v(973))
       v(289)=v(969)*(v(142)*v(222)+v(144)*v(973))
       v(283)=v(969)*(v(140)*v(216)+v(139)*v(973))
       v(276)=v(969)*(v(139)*v(222)+v(140)*v(973))
       v(151)=-v(149)-v(150)
       v(974)=v(151)*v(162)
       v(218)=v(151)*v(164)
       v(223)=v(218)+2d0*v(974)
       v(154)=-v(152)-v(153)
       v(975)=v(154)*v(162)
       v(219)=v(154)*v(164)
       v(306)=v(969)*(v(151)*v(219)+v(154)*v(974))
       v(301)=v(969)*(v(146)*v(219)+v(148)*v(974))
       v(292)=v(969)*(v(142)*v(219)+v(144)*v(974))
       v(279)=v(969)*(v(139)*v(219)+v(140)*v(974))
       v(235)=v(219)+2d0*v(975)
       v(304)=v(969)*(v(148)*v(235)+v(146)*v(974))
       v(297)=v(969)*(v(144)*v(235)+v(142)*v(974))
       v(286)=v(969)*(v(140)*v(235)+v(139)*v(974))
       v(303)=v(969)*(v(148)*v(218)+v(146)*v(975))
       v(300)=v(969)*(v(146)*v(223)+v(148)*v(975))
       v(296)=v(969)*(v(144)*v(218)+v(142)*v(975))
       v(291)=v(969)*(v(142)*v(223)+v(144)*v(975))
       v(285)=v(969)*(v(140)*v(218)+v(139)*v(975))
       v(278)=v(969)*(v(139)*v(223)+v(140)*v(975))
       v(159)=v(139)*v(66)+v(142)*v(68)+v(146)*v(70)+v(151)*v(72)
       v(161)=v(140)*v(67)+v(144)*v(69)+v(148)*v(71)+v(154)*v(73)
       v(166)=(v(159)+v(161))*v(164)
       v(167)=v(166)+v(159)*v(976)
       v(168)=v(162)*(v(140)*v(66)+v(139)*v(67)+v(144)*v(68)+v(142)*v
     & (69)+v(148)*v(70)+v(146)*v(71)+v(154)*v(72)+v(151)*v(73))
       v(169)=v(166)+v(161)*v(976)
       p(1)=p(1)+v(969)*(v(139)*v(167)+v(140)*v(168)+v(94)*v(977))
       p(2)=p(2)+v(969)*(v(139)*v(168)+v(140)*v(169)+v(94)*v(978))
       p(3)=p(3)+v(969)*(v(142)*v(167)+v(144)*v(168)+v(96)*v(977))
       p(4)=p(4)+v(969)*(v(142)*v(168)+v(144)*v(169)+v(96)*v(978))
       p(5)=p(5)+v(969)*(v(146)*v(167)+v(148)*v(168)+v(977)*v(98))
       p(6)=p(6)+v(969)*(v(146)*v(168)+v(148)*v(169)+v(978)*v(98))
       p(7)=p(7)+v(969)*(v(151)*v(167)+v(154)*v(168)+v(101)*v(977))
       p(8)=p(8)+v(969)*(v(151)*v(168)+v(154)*v(169)+v(101)*v(978))
       s(1,1)=s(1,1)+v(969)*((v(140)*v(140))*v(162)+(v(139)*v(139))*v
     & (987))
       s(1,2)=s(1,2)+v(273)
       s(1,3)=s(1,3)+v(274)
       s(1,4)=s(1,4)+v(275)
       s(1,5)=s(1,5)+v(276)
       s(1,6)=s(1,6)+v(277)
       s(1,7)=s(1,7)+v(278)
       s(1,8)=s(1,8)+v(279)
       s(2,1)=s(2,1)+v(273)
       s(2,2)=s(2,2)+v(969)*(v(139)*v(225)+v(140)*(v(213)+v(140)*v
     & (976)))
       s(2,3)=s(2,3)+v(281)
       s(2,4)=s(2,4)+v(282)
       s(2,5)=s(2,5)+v(283)
       s(2,6)=s(2,6)+v(284)
       s(2,7)=s(2,7)+v(285)
       s(2,8)=s(2,8)+v(286)
       s(3,1)=s(3,1)+v(274)
       s(3,2)=s(3,2)+v(281)
       s(3,3)=s(3,3)+v(969)*(v(142)*v(221)+v(144)*v(971))
       s(3,4)=s(3,4)+v(288)
       s(3,5)=s(3,5)+v(289)
       s(3,6)=s(3,6)+v(290)
       s(3,7)=s(3,7)+v(291)
       s(3,8)=s(3,8)+v(292)
       s(4,1)=s(4,1)+v(275)
       s(4,2)=s(4,2)+v(282)
       s(4,3)=s(4,3)+v(288)
       s(4,4)=s(4,4)+v(969)*(v(144)*v(233)+v(142)*v(970))
       s(4,5)=s(4,5)+v(294)
       s(4,6)=s(4,6)+v(295)
       s(4,7)=s(4,7)+v(296)
       s(4,8)=s(4,8)+v(297)
       s(5,1)=s(5,1)+v(276)
       s(5,2)=s(5,2)+v(283)
       s(5,3)=s(5,3)+v(289)
       s(5,4)=s(5,4)+v(294)
       s(5,5)=s(5,5)+v(969)*(v(146)*v(222)+v(148)*v(973))
       s(5,6)=s(5,6)+v(299)
       s(5,7)=s(5,7)+v(300)
       s(5,8)=s(5,8)+v(301)
       s(6,1)=s(6,1)+v(277)
       s(6,2)=s(6,2)+v(284)
       s(6,3)=s(6,3)+v(290)
       s(6,4)=s(6,4)+v(295)
       s(6,5)=s(6,5)+v(299)
       s(6,6)=s(6,6)+v(969)*(v(148)*v(234)+v(146)*v(972))
       s(6,7)=s(6,7)+v(303)
       s(6,8)=s(6,8)+v(304)
       s(7,1)=s(7,1)+v(278)
       s(7,2)=s(7,2)+v(285)
       s(7,3)=s(7,3)+v(291)
       s(7,4)=s(7,4)+v(296)
       s(7,5)=s(7,5)+v(300)
       s(7,6)=s(7,6)+v(303)
       s(7,7)=s(7,7)+v(969)*(v(151)*v(223)+v(154)*v(975))
       s(7,8)=s(7,8)+v(306)
       s(8,1)=s(8,1)+v(279)
       s(8,2)=s(8,2)+v(286)
       s(8,3)=s(8,3)+v(292)
       s(8,4)=s(8,4)+v(297)
       s(8,5)=s(8,5)+v(301)
       s(8,6)=s(8,6)+v(304)
       s(8,7)=s(8,7)+v(306)
       s(8,8)=s(8,8)+v(969)*(v(154)*v(235)+v(151)*v(974))
       c(1,1)=0d0+c(1,1)
       c(1,2)=0d0+c(1,2)
       c(1,3)=0d0+c(1,3)
       c(1,4)=0d0+c(1,4)
       c(1,5)=0d0+c(1,5)
       c(1,6)=0d0+c(1,6)
       c(1,7)=0d0+c(1,7)
       c(1,8)=0d0+c(1,8)
       c(2,1)=0d0+c(2,1)
       c(2,2)=0d0+c(2,2)
       c(2,3)=0d0+c(2,3)
       c(2,4)=0d0+c(2,4)
       c(2,5)=0d0+c(2,5)
       c(2,6)=0d0+c(2,6)
       c(2,7)=0d0+c(2,7)
       c(2,8)=0d0+c(2,8)
       c(3,1)=0d0+c(3,1)
       c(3,2)=0d0+c(3,2)
       c(3,3)=0d0+c(3,3)
       c(3,4)=0d0+c(3,4)
       c(3,5)=0d0+c(3,5)
       c(3,6)=0d0+c(3,6)
       c(3,7)=0d0+c(3,7)
       c(3,8)=0d0+c(3,8)
       c(4,1)=0d0+c(4,1)
       c(4,2)=0d0+c(4,2)
       c(4,3)=0d0+c(4,3)
       c(4,4)=0d0+c(4,4)
       c(4,5)=0d0+c(4,5)
       c(4,6)=0d0+c(4,6)
       c(4,7)=0d0+c(4,7)
       c(4,8)=0d0+c(4,8)
       c(5,1)=0d0+c(5,1)
       c(5,2)=0d0+c(5,2)
       c(5,3)=0d0+c(5,3)
       c(5,4)=0d0+c(5,4)
       c(5,5)=0d0+c(5,5)
       c(5,6)=0d0+c(5,6)
       c(5,7)=0d0+c(5,7)
       c(5,8)=0d0+c(5,8)
       c(6,1)=0d0+c(6,1)
       c(6,2)=0d0+c(6,2)
       c(6,3)=0d0+c(6,3)
       c(6,4)=0d0+c(6,4)
       c(6,5)=0d0+c(6,5)
       c(6,6)=0d0+c(6,6)
       c(6,7)=0d0+c(6,7)
       c(6,8)=0d0+c(6,8)
       c(7,1)=0d0+c(7,1)
       c(7,2)=0d0+c(7,2)
       c(7,3)=0d0+c(7,3)
       c(7,4)=0d0+c(7,4)
       c(7,5)=0d0+c(7,5)
       c(7,6)=0d0+c(7,6)
       c(7,7)=0d0+c(7,7)
       c(7,8)=0d0+c(7,8)
       c(8,1)=0d0+c(8,1)
       c(8,2)=0d0+c(8,2)
       c(8,3)=0d0+c(8,3)
       c(8,4)=0d0+c(8,4)
       c(8,5)=0d0+c(8,5)
       c(8,6)=0d0+c(8,6)
       c(8,7)=0d0+c(8,7)
       c(8,8)=0d0+c(8,8)
       m(1,1)=m(1,1)+v(317)
       m(1,2)=0d0+m(1,2)
       m(1,3)=m(1,3)+v(318)
       m(1,4)=0d0+m(1,4)
       m(1,5)=m(1,5)+v(319)
       m(1,6)=0d0+m(1,6)
       m(1,7)=m(1,7)+v(320)
       m(1,8)=0d0+m(1,8)
       m(2,1)=0d0+m(2,1)
       m(2,2)=m(2,2)+v(317)
       m(2,3)=0d0+m(2,3)
       m(2,4)=m(2,4)+v(318)
       m(2,5)=0d0+m(2,5)
       m(2,6)=m(2,6)+v(319)
       m(2,7)=0d0+m(2,7)
       m(2,8)=m(2,8)+v(320)
       m(3,1)=m(3,1)+v(318)
       m(3,2)=0d0+m(3,2)
       m(3,3)=m(3,3)+v(321)
       m(3,4)=0d0+m(3,4)
       m(3,5)=m(3,5)+v(322)
       m(3,6)=0d0+m(3,6)
       m(3,7)=m(3,7)+v(319)
       m(3,8)=0d0+m(3,8)
       m(4,1)=0d0+m(4,1)
       m(4,2)=m(4,2)+v(318)
       m(4,3)=0d0+m(4,3)
       m(4,4)=m(4,4)+v(321)
       m(4,5)=0d0+m(4,5)
       m(4,6)=m(4,6)+v(322)
       m(4,7)=0d0+m(4,7)
       m(4,8)=m(4,8)+v(319)
       m(5,1)=m(5,1)+v(319)
       m(5,2)=0d0+m(5,2)
       m(5,3)=m(5,3)+v(322)
       m(5,4)=0d0+m(5,4)
       m(5,5)=m(5,5)+v(323)
       m(5,6)=0d0+m(5,6)
       m(5,7)=m(5,7)+v(324)
       m(5,8)=0d0+m(5,8)
       m(6,1)=0d0+m(6,1)
       m(6,2)=m(6,2)+v(319)
       m(6,3)=0d0+m(6,3)
       m(6,4)=m(6,4)+v(322)
       m(6,5)=0d0+m(6,5)
       m(6,6)=m(6,6)+v(323)
       m(6,7)=0d0+m(6,7)
       m(6,8)=m(6,8)+v(324)
       m(7,1)=m(7,1)+v(320)
       m(7,2)=0d0+m(7,2)
       m(7,3)=m(7,3)+v(319)
       m(7,4)=0d0+m(7,4)
       m(7,5)=m(7,5)+v(324)
       m(7,6)=0d0+m(7,6)
       m(7,7)=m(7,7)+v(325)
       m(7,8)=0d0+m(7,8)
       m(8,1)=0d0+m(8,1)
       m(8,2)=m(8,2)+v(320)
       m(8,3)=0d0+m(8,3)
       m(8,4)=m(8,4)+v(319)
       m(8,5)=0d0+m(8,5)
       m(8,6)=m(8,6)+v(324)
       m(8,7)=0d0+m(8,7)
       m(8,8)=m(8,8)+v(325)
      ENDDO
      END

!******************* S U B R O U T I N E **********************
      SUBROUTINE elmt16_ISW05(v,d,xl,ul,m,p,ht,hp,gp,ngpo)
      IMPLICIT NONE
      include 'sms.h'
      INTEGER ngpo,i329
      DOUBLE PRECISION v(1125),d(3),xl(2,4),ul(6,8),m(8,8),p
     &(8),ht(0),hp(0),gp(4,4)
      v(340)=xl(2,4)
      v(339)=xl(1,4)
      v(338)=xl(2,3)
      v(998)=v(338)-v(340)
      v(337)=xl(1,3)
      v(994)=v(337)-v(339)
      v(336)=xl(2,2)
      v(992)=v(336)-v(338)
      v(335)=xl(1,2)
      v(996)=v(335)-v(337)
      v(334)=xl(2,1)
      v(997)=v(334)-v(336)
      v(991)=v(334)-v(340)
      v(333)=xl(1,1)
      v(995)=v(333)-v(339)
      v(993)=v(333)-v(335)
      DO i329=1,ngpo
       v(413)=gp(1,i329)
       v(422)=1d0-v(413)
       v(428)=(-0.25d0)*v(422)
       v(420)=1d0+v(413)
       v(429)=(-0.25d0)*v(420)
       v(414)=gp(2,i329)
       v(430)=(1d0+v(414))/4d0
       v(427)=((-1d0)+v(414))/4d0
       v(417)=-(v(422)*v(427))
       v(419)=-(v(420)*v(427))
       v(421)=v(420)*v(430)
       v(424)=v(422)*v(430)
       v(988)=d(3)*gp(4,i329)*((v(428)*v(991)+v(429)*v(992))*(v(427
     & )*v(993)+v(430)*v(994))-(v(428)*v(995)+v(429)*v(996))*(v(427
     & )*v(997)+v(430)*v(998)))
       v(990)=v(421)*v(988)
       v(989)=v(417)*v(988)
       v(551)=(v(424)*v(424))*v(988)
       v(549)=(v(421)*v(421))*v(988)
       v(547)=(v(419)*v(419))*v(988)
       v(543)=(v(417)*v(417))*v(988)
       v(546)=v(424)*v(989)
       v(545)=v(421)*v(989)
       v(544)=v(419)*v(989)
       v(548)=v(419)*v(990)
       v(550)=v(424)*v(990)
       m(1,1)=m(1,1)+v(543)
       m(1,2)=0d0+m(1,2)
       m(1,3)=m(1,3)+v(544)
       m(1,4)=0d0+m(1,4)
       m(1,5)=m(1,5)+v(545)
       m(1,6)=0d0+m(1,6)
       m(1,7)=m(1,7)+v(546)
       m(1,8)=0d0+m(1,8)
       m(2,1)=0d0+m(2,1)
       m(2,2)=m(2,2)+v(543)
       m(2,3)=0d0+m(2,3)
       m(2,4)=m(2,4)+v(544)
       m(2,5)=0d0+m(2,5)
       m(2,6)=m(2,6)+v(545)
       m(2,7)=0d0+m(2,7)
       m(2,8)=m(2,8)+v(546)
       m(3,1)=m(3,1)+v(544)
       m(3,2)=0d0+m(3,2)
       m(3,3)=m(3,3)+v(547)
       m(3,4)=0d0+m(3,4)
       m(3,5)=m(3,5)+v(548)
       m(3,6)=0d0+m(3,6)
       m(3,7)=m(3,7)+v(545)
       m(3,8)=0d0+m(3,8)
       m(4,1)=0d0+m(4,1)
       m(4,2)=m(4,2)+v(544)
       m(4,3)=0d0+m(4,3)
       m(4,4)=m(4,4)+v(547)
       m(4,5)=0d0+m(4,5)
       m(4,6)=m(4,6)+v(548)
       m(4,7)=0d0+m(4,7)
       m(4,8)=m(4,8)+v(545)
       m(5,1)=m(5,1)+v(545)
       m(5,2)=0d0+m(5,2)
       m(5,3)=m(5,3)+v(548)
       m(5,4)=0d0+m(5,4)
       m(5,5)=m(5,5)+v(549)
       m(5,6)=0d0+m(5,6)
       m(5,7)=m(5,7)+v(550)
       m(5,8)=0d0+m(5,8)
       m(6,1)=0d0+m(6,1)
       m(6,2)=m(6,2)+v(545)
       m(6,3)=0d0+m(6,3)
       m(6,4)=m(6,4)+v(548)
       m(6,5)=0d0+m(6,5)
       m(6,6)=m(6,6)+v(549)
       m(6,7)=0d0+m(6,7)
       m(6,8)=m(6,8)+v(550)
       m(7,1)=m(7,1)+v(546)
       m(7,2)=0d0+m(7,2)
       m(7,3)=m(7,3)+v(545)
       m(7,4)=0d0+m(7,4)
       m(7,5)=m(7,5)+v(550)
       m(7,6)=0d0+m(7,6)
       m(7,7)=m(7,7)+v(551)
       m(7,8)=0d0+m(7,8)
       m(8,1)=0d0+m(8,1)
       m(8,2)=m(8,2)+v(546)
       m(8,3)=0d0+m(8,3)
       m(8,4)=m(8,4)+v(545)
       m(8,5)=0d0+m(8,5)
       m(8,6)=m(8,6)+v(550)
       m(8,7)=0d0+m(8,7)
       m(8,8)=m(8,8)+v(551)
      ENDDO
      END

!******************* S U B R O U T I N E **********************
      SUBROUTINE elmt16_ISW06(v,d,xl,ul,p,ht,hp,gp,ngpo)
      IMPLICIT NONE
      include 'sms.h'
      INTEGER ngpo,i553
      DOUBLE PRECISION v(1125),d(3),xl(2,4),ul(6,8),p(8),ht
     &(0),hp(0),gp(4,4)
      v(620)=ul(1,8)
      v(619)=ul(1,7)
      v(618)=ul(1,6)
      v(617)=ul(1,5)
      v(616)=ul(1,4)
      v(615)=ul(1,3)
      v(614)=ul(1,2)
      v(613)=ul(1,1)
      v(564)=xl(2,4)
      v(563)=xl(1,4)
      v(562)=xl(2,3)
      v(1009)=v(562)-v(564)
      v(561)=xl(1,3)
      v(1011)=v(561)-v(563)
      v(560)=xl(2,2)
      v(1005)=v(560)-v(562)
      v(559)=xl(1,2)
      v(1007)=v(559)-v(561)
      v(558)=xl(2,1)
      v(1008)=v(558)-v(560)
      v(1004)=v(558)-v(564)
      v(557)=xl(1,1)
      v(1010)=v(557)-v(559)
      v(1006)=v(557)-v(563)
      v(556)=d(3)
      v(555)=d(2)
      v(999)=d(1)/(1d0+v(555))
      v(711)=(v(555)*v(999))/(1d0-2d0*v(555))
      v(709)=v(999)/2d0
      v(1000)=2d0*v(709)
      DO i553=1,ngpo
       v(637)=gp(1,i553)
       v(646)=1d0-v(637)
       v(652)=(-0.25d0)*v(646)
       v(644)=1d0+v(637)
       v(653)=(-0.25d0)*v(644)
       v(658)=v(1004)*v(652)+v(1005)*v(653)
       v(656)=v(1006)*v(652)+v(1007)*v(653)
       v(638)=gp(2,i553)
       v(654)=(1d0+v(638))/4d0
       v(651)=((-1d0)+v(638))/4d0
       v(657)=v(1008)*v(651)+v(1009)*v(654)
       v(655)=v(1010)*v(651)+v(1011)*v(654)
       v(641)=-(v(646)*v(651))
       v(643)=-(v(644)*v(651))
       v(645)=v(644)*v(654)
       v(648)=v(646)*v(654)
       v(659)=-(v(656)*v(657))+v(655)*v(658)
       v(1001)=gp(4,i553)*v(659)
       v(1002)=v(556)*(ul(5,1)*v(641)+ul(5,3)*v(643)+ul(5,5)*v(645)
     & +ul(5,7)*v(648))
       v(1003)=v(556)*(ul(5,2)*v(641)+ul(5,4)*v(643)+ul(5,6)*v(645)
     & +ul(5,8)*v(648))
       v(682)=-(v(658)/v(659))
       v(696)=-(v(654)*v(682))
       v(688)=-(v(651)*v(682))
       v(683)=v(656)/v(659)
       v(699)=-(v(654)*v(683))
       v(690)=-(v(651)*v(683))
       v(684)=-(v(657)/v(659))
       v(697)=v(652)*v(684)
       v(692)=v(653)*v(684)
       v(685)=v(655)/v(659)
       v(700)=v(652)*v(685)
       v(694)=v(653)*v(685)
       v(686)=v(688)+v(697)
       v(687)=v(690)+v(700)
       v(689)=-v(688)+v(692)
       v(691)=-v(690)+v(694)
       v(693)=-v(692)+v(696)
       v(695)=-v(694)+v(699)
       v(698)=-v(696)-v(697)
       v(701)=-v(699)-v(700)
       v(706)=v(613)*v(686)+v(615)*v(689)+v(617)*v(693)+v(619)*v(698)
       v(708)=v(614)*v(687)+v(616)*v(691)+v(618)*v(695)+v(620)*v(701)
       v(713)=(v(706)+v(708))*v(711)
       v(714)=v(1000)*v(706)+v(713)
       v(715)=(v(614)*v(686)+v(613)*v(687)+v(616)*v(689)+v(615)*v(691
     & )+v(618)*v(693)+v(617)*v(695)+v(620)*v(698)+v(619)*v(701))*v
     & (709)
       v(716)=v(1000)*v(708)+v(713)
       p(1)=p(1)+v(1001)*(v(1002)*v(641)+v(686)*v(714)+v(687)*v(715))
       p(2)=p(2)+v(1001)*(v(1003)*v(641)+v(686)*v(715)+v(687)*v(716))
       p(3)=p(3)+v(1001)*(v(1002)*v(643)+v(689)*v(714)+v(691)*v(715))
       p(4)=p(4)+v(1001)*(v(1003)*v(643)+v(689)*v(715)+v(691)*v(716))
       p(5)=p(5)+v(1001)*(v(1002)*v(645)+v(693)*v(714)+v(695)*v(715))
       p(6)=p(6)+v(1001)*(v(1003)*v(645)+v(693)*v(715)+v(695)*v(716))
       p(7)=p(7)+v(1001)*(v(1002)*v(648)+v(698)*v(714)+v(701)*v(715))
       p(8)=p(8)+v(1001)*(v(1003)*v(648)+v(698)*v(715)+v(701)*v(716))
      ENDDO
      END

!******************* S U B R O U T I N E **********************
      SUBROUTINE elmt16_ISW09(v,d,xl,ul,c,p,ht,hp,gp,ngpo)
      IMPLICIT NONE
      include 'sms.h'
      INTEGER ngpo,i759
      DOUBLE PRECISION v(1125),d(3),xl(2,4),ul(6,8),c(8,8),p
     &(8),ht(0),hp(0),gp(4,4)
      DO i759=1,ngpo
       c(1,1)=0d0+c(1,1)
       c(1,2)=0d0+c(1,2)
       c(1,3)=0d0+c(1,3)
       c(1,4)=0d0+c(1,4)
       c(1,5)=0d0+c(1,5)
       c(1,6)=0d0+c(1,6)
       c(1,7)=0d0+c(1,7)
       c(1,8)=0d0+c(1,8)
       c(2,1)=0d0+c(2,1)
       c(2,2)=0d0+c(2,2)
       c(2,3)=0d0+c(2,3)
       c(2,4)=0d0+c(2,4)
       c(2,5)=0d0+c(2,5)
       c(2,6)=0d0+c(2,6)
       c(2,7)=0d0+c(2,7)
       c(2,8)=0d0+c(2,8)
       c(3,1)=0d0+c(3,1)
       c(3,2)=0d0+c(3,2)
       c(3,3)=0d0+c(3,3)
       c(3,4)=0d0+c(3,4)
       c(3,5)=0d0+c(3,5)
       c(3,6)=0d0+c(3,6)
       c(3,7)=0d0+c(3,7)
       c(3,8)=0d0+c(3,8)
       c(4,1)=0d0+c(4,1)
       c(4,2)=0d0+c(4,2)
       c(4,3)=0d0+c(4,3)
       c(4,4)=0d0+c(4,4)
       c(4,5)=0d0+c(4,5)
       c(4,6)=0d0+c(4,6)
       c(4,7)=0d0+c(4,7)
       c(4,8)=0d0+c(4,8)
       c(5,1)=0d0+c(5,1)
       c(5,2)=0d0+c(5,2)
       c(5,3)=0d0+c(5,3)
       c(5,4)=0d0+c(5,4)
       c(5,5)=0d0+c(5,5)
       c(5,6)=0d0+c(5,6)
       c(5,7)=0d0+c(5,7)
       c(5,8)=0d0+c(5,8)
       c(6,1)=0d0+c(6,1)
       c(6,2)=0d0+c(6,2)
       c(6,3)=0d0+c(6,3)
       c(6,4)=0d0+c(6,4)
       c(6,5)=0d0+c(6,5)
       c(6,6)=0d0+c(6,6)
       c(6,7)=0d0+c(6,7)
       c(6,8)=0d0+c(6,8)
       c(7,1)=0d0+c(7,1)
       c(7,2)=0d0+c(7,2)
       c(7,3)=0d0+c(7,3)
       c(7,4)=0d0+c(7,4)
       c(7,5)=0d0+c(7,5)
       c(7,6)=0d0+c(7,6)
       c(7,7)=0d0+c(7,7)
       c(7,8)=0d0+c(7,8)
       c(8,1)=0d0+c(8,1)
       c(8,2)=0d0+c(8,2)
       c(8,3)=0d0+c(8,3)
       c(8,4)=0d0+c(8,4)
       c(8,5)=0d0+c(8,5)
       c(8,6)=0d0+c(8,6)
       c(8,7)=0d0+c(8,7)
       c(8,8)=0d0+c(8,8)
      ENDDO
      END

max_stars_repo_name: bhajay/Benchmark-CoupledProblems

max_issues_repo_name: bhajay/Benchmark-CoupledProblems

max_forks_repo_name: bhajay/Benchmark-CoupledProblems


text:
Name: Benjamin Rosenstein

Personality: Life long Townies Townie who graduated from UCD with a degree in Political Science in March of 2008.  It seems very likely that he will finally leave Davis after 22.5 years in August.

I am going to hopefully write up a review of every restaurant (I really have eaten at most of them), before I leave town.  I like the 4 star system, even though with the half stars its more like the 8 star system.  Oh, and I do think cost is big factor in my food reviews (ie I expect a lot less for $5 than $20, so thats why a fancy restaurant might get less stars than a cheap one).

American

Applebees  1.5 Stars

Its bland greesy food.  I wonder if I was put into a Applebees if I could point out anything differnant than a TGIFs or Shaniganz.  The only thing here I have every really liked are the acholic milkshakes.  These things can be really good, I like the modified mudsilde in particular.

Bakers Square

Burgers and Brew

Cindys

Froggys

The Graduate

Hot Dogger

IHOP

Plainfield Station

Sudwerk

Tommy Js

Wildhorse Grill

20060123 08:59:57 nbsp Whoa, when did your hair get so... floofy?  :D Users/ErinManley

20060130 00:24:09 nbsp that hair makes you look like a troll doll.  i mean that in the best way possible.  i love you ben. Users/JacobRoche

20060827 12:10:40 nbsp Hiya Ben.  You still working at Steves? Users/ColinFreeman

20061106 14:38:51 nbsp YAYA Users/StevenDaubert

20070104 13:25:02 nbsp pppsst Ben, no one cares about ASUCD! Users/StevenDaubert

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
c-----------------------------------------------------------------------
      program test

      include 'ceedf.h'

      integer ceed,err
      integer x,y
      integer r
      integer i,n
      integer*8 offset

      integer ne
      parameter(ne=3)

      integer*4 ind(2*ne)
      real*8 a(2*(ne+1))
      real*8 yy(4*ne)
      real*8 diff

      character arg*32

      call getarg(1,arg)
      call ceedinit(trim(arg)//char(0),ceed,err)

      call ceedvectorcreate(ceed,2*(ne+1),x,err)

      do i=1,ne+1
        a(2*i-1)=10+i-1
        a(2*i)=20+i-1
      enddo

      call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,a,err)

      do i=1,ne
        ind(2*i-1)=i-1
        ind(2*i  )=i
      enddo

      call ceedelemrestrictioncreate(ceed,ne,2,ne+1,2,ceed_mem_host,
     $  ceed_use_pointer,ind,r,err)

      call ceedvectorcreate(ceed,2*(2*ne),y,err);
      call ceedvectorsetvalue(y,0.d0,err);
      call ceedelemrestrictionapply(r,ceed_notranspose,
     $  ceed_transpose,x,y,ceed_request_immediate,err)

      call ceedvectorgetarrayread(y,ceed_mem_host,yy,offset,err)
      do i=0,ne-1
        do n=1,2
          diff=10+(2*i+n)/2-yy(i*4+n+offset)
          if (abs(diff) > 1.0D-15) then
            write(*,*) 'Error in restricted array y(',i*4+n,')=',
     $  yy(i*4+n+offset),'!=',10+(2*i+n)/2
          endif
          diff=20+(2*i+n)/2-yy(i*4+n+2+offset)
          if (abs(diff) > 1.0D-15) then
            write(*,*) 'Error in restricted array y(',i*4+n+2,')=',
     $  yy(i*4+n+2+offset),'!=',20+(2*i+n)/2
          endif
        enddo
      enddo
      call ceedvectorrestorearrayread(y,yy,offset,err)

      call ceedvectordestroy(x,err)
      call ceedvectordestroy(y,err)
      call ceedelemrestrictiondestroy(r,err)
      call ceeddestroy(ceed,err)

      end
c-----------------------------------------------------------------------

max_stars_repo_name: YohannDudouit/libCEED

max_issues_repo_name: YohannDudouit/libCEED

max_forks_repo_name: YohannDudouit/libCEED


text:
C Read two ints and two floats from a file using one FORMAT,
C then write them out to another file using a different FORMAT.
C The expected output written to file 'outfile3' are the
C following four lines:
C
C F =  23.5; I =  456
C
C F =  67.9; I =  123

      PROGRAM MAIN

      REAL X,Y
      OPEN (10, FILE="infile3")
      OPEN (20, FILE="outfile3", STATUS="REPLACE")
 30   FORMAT(/,'F = ', F5.1, '; I = ', I4)
      READ (10,10) I,X,J,Y
      WRITE (20,30) X,J
      WRITE (20,30) Y,I
 10   FORMAT (2(I3,X,F5.2,X))
      STOP
      END PROGRAM MAIN

ext: for


!---------------------------------------------------!
! Copyright (c) 2016 Shunsuke A. Sato               !
! Released under the MIT license                    !
! https://opensource.org/licenses/mit-license.php   !
!---------------------------------------------------!
!-------10--------20--------30--------40--------50--------60--------70--------80--------90
subroutine energy(Eex)
  use global_variables
  implicit none
  real(8) :: Eex
  integer :: ikr,ikz

  Eex = 0d0
!$omp parallel
!$omp do private(ikz, ikr) reduction(+:Eex)
  do ikz = -NKz,NKz
  do ikr = 1,NKr
    Eex = Eex+ deps(ikr,ikz)*abs(zCt(2,ikr,ikz))**2*kr(ikr)
  end do
  end do
!$omp end parallel
  Eex = Eex*2d0/((2d0*pi)**3)*(2d0*pi*dkr*dkz)

  return
end subroutine energy


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

!
! **********************************************************************
! * Begin subroutine to calculate solar intensity and efficiency       *
! * matrix for the phytoplankton model bloom ii                        *
! * the program has been updated most recently in december 2015        *
! * to be part of delwaq/bloom as a subroutine                         *
! * the user may specify the photosynthetic curves tablelized          *
! **********************************************************************
!
      subroutine bleffpro(lunrep, lunblm, nuecog, nz, zvec, fun, der)

      implicit none
      dimension power(51),solvec(51),rho(51,30),fun(51,30),der(51,30),
     1          sfirst(51),solar(51),time(51),cdf(51),freq(51),
     2          dens(51),tsol(51),tden(51),domf(51),rfirst(51,30),
     3          gfun(51),gder(51),zvec(51),daymul(24,30),dl(24)
      common/effcny/rho,power
      common/solrad/tsol,tden,freq,nval
      common/solval/delsol,solmax,day
      character*5 table
      character*265 inputfile
      character*265 outputfile
      integer lunblm, lunrep
      integer irc, npoint, i, j, k, nsp, nval, nz
      real*8 power, rho, solvec, time, solar, cdf, dens, freq, tsol, tden, domf
      real*8 rfirst, zvec, done, dneg, sfirst, gfun, gder, fun, der
      real*8 daymul, dl, delsol, solmax, day
!
      integer nuecog ! Number of groups
!
!  Read number of light intensity points
!  Read efficiency data for each ecogroup
!
      read (lunblm,*) npoint
!
!  Read in tables
!
   15 continue
      do 20 i=1,npoint
   20 read (lunblm,*) power(i),(rho(i,j),j=1,nuecog)
!
!  Read, integrate, and transform diurnal intensity distribution
!  Read number of points in solar radiation distribution
!
      read (lunblm,*) nsp
!
!  Read: solvec--solar radiation level; time--time of day
!
      read (lunblm,*) (solvec(k),time(k),k=1,nsp)
!
!  nval is the number of points desired in the transformed solar
!   radiation distribution function (max 50)
!
      nval=50
      call solcdf(solvec,time,nsp,nval,solar,cdf)
      call convrt(solar,cdf,nval,dens)
      call indist(nval,solar,dens,freq,tsol,tden)
!
!  Transform efficiency data
!
      call inteff(npoint,nuecog,domf,rfirst)
!
!  Determine appropriate tabulation points for convolution
!  nz is the number of points desired in the convolution (max 51)
!  Write convolution points into data set
!
      nz=51
      call zval(domf,npoint,tsol,nval,zvec,nz)
!
!  Compute convolutions and their derivatives
!  Calculate constants and determine roots for each order
!
      done = 1.0d0
      dneg = -1.0d0
      do 80 j=1,nuecog
      do 60 k=1,npoint
   60 sfirst(k)=rfirst(k,j)
      call cvolve(tsol,freq,nval,done,domf,sfirst,npoint,dneg,zvec,
     1            gfun,nz)
      call cvolve(domf,sfirst,npoint,dneg,tsol,tden,nval,done,zvec,
     1            gder,nz)
      do 70 k=1,nz
      fun(k,j)=gfun(k)
   70 der(k,j)=gder(k)
   80 continue
      return
      end subroutine
!
!  Subroutine to adjust lower end of efficiency data and transform it
!  to a logarithmic form
!
      subroutine inteff(npoint,nuecog,domf,rfirst)
      implicit none
      dimension effic(51,30),power(51),domf(51),rfirst(51,30),e(30)
      common/effcny/effic,power
      integer npoint, nuecog, i, j, i1
      real*8 power, p, effic, e, rfirst, domf
!
!  Move lower end of curve away from zero.
!
    1 if (power(1) .gt. 0.001*power(2)) go to 15
      p=0.001*power(2)
      call interm(power,effic,npoint,nuecog,p,e)
      do 12 j=1,nuecog
   12 effic(1,j)=e(j)
      power(1)=p
!
!  Transform intensity into its logarithmic form
!
   15 do 20 i=1,npoint
      i1=npoint-i+1
      do 18 j=1,nuecog
   18 rfirst(i,j)=effic(i1,j)
   20 domf(i)=-dlog(power(i1))
      return
      end
!
!  Subroutine to find cumulative distribution function
!
      subroutine solcdf(solvec,time,nsp,nval,solar,cdf)
      implicit none
      dimension solvec(51),time(51),cdf(51),solar(51),solna(51),
     1          solnb(51),timea(51),timeb(51),value(2)
      common/solval/delsol,solmax,day
      integer i, nsp, it, na, nb, j, jk, jkn, n, nval
      real*8 solmax, solvec, rn, delsol, day, time, solna, timea, solnb
      real*8 timeb, solar, cdf, rj, sol, value, dens, ri
!
!  Find maximum intensity
!
    1 solmax=solvec(1)
      do 5 i=1,nsp
    5 if (solvec(i) .gt. solmax) solmax=solvec(i)
!
!  SET CONSTANTS
!
      rn=nval-1
      delsol=solmax/rn
      day=time(nsp)-time(1)
!
!   SPLIT SOLAR RADIATION VECTOR INTO ASCENDING AND DESCENDING PARTS
!
      it=0
      i=1
   12 if (solvec(i) .eq. solmax) it=1
      solna(i)=solvec(i)
      timea(i)=time(i)
      if (it .eq. 1) go to 15
      i=i+1
      go to 12
   15 na=i
      nb=nsp-na+1
      do 18 j=1,nb
      jk=j-1
      jkn=nsp-jk
      solnb(j)=solvec(jkn)
   18 timeb(j)=time(jkn)
!
!  Calculate cdf for radiation function
!  Set endpoints of arrays
!
      solar(1)=0.0
      cdf(1)=0.0
!
!  Begin calculation of cdf
!
      do 20 j=2,nval
      rj=j-1
      sol=delsol*rj
      call interp(solna,timea,na,sol,value(1))
      call interp(solnb,timeb,nb,sol,value(2))
      cdf(j)=(value(1)-time(1)+time(nsp)-value(2))/day
   20 solar(j)=sol
      return
      end
!
!  Subroutine to convert cdf to probability density function
!
      subroutine convrt(solar,cdf,nval,dens)
      implicit none
      dimension solar(51),cdf(51),dens(51),value(2)
      common/solval/delsol,solmax,day
      integer i, n, nval, imin
      real*8 rj, sol, value, dens, ri, solc, delsol, soll, solu
      real*8 solar, cdf, solmax, day
!
!  Set constants and endpoints
!
    1 dens(1)=0.0
!
!  Iterate through intensity function
!  Determine density values
!
      n=nval-1
      do 10 i=2,n
      ri=i-1
      solc=delsol*ri
      if (i .eq. 2) go to 6
      value(1)=value(2)
      soll=solu
    4 solu=solc+0.5*delsol
      call interp(solar,cdf,nval,solu,value(2))
      go to 8
    6 soll=0.5*delsol
      call interp(solar,cdf,nval,soll,value(1))
      go to 4
    8 dens(i)=(value(2)-value(1))/delsol
   10 continue
      dens(nval)=(1.0-value(2))*2.0/delsol
      return
      end
!
!  Subroutine to integrate and normalize diurnal intensity
!  distribution, and compute mean intensity
!
      subroutine indist(nval,solar,dens,freq,tsol,tden)
      implicit none
      dimension solar(51),freq(51),dens(51),tsol(51),tden(51)
      integer i, nval
      real*8 solar, s, dens, d, freq, sbar, del, amult, tsol, tden
!
!  Move lower end of curve away from zero
!
      if (solar(1) .gt. 0.001*solar(2)) go to 15
      s=0.001*solar(2)
      call interp(solar,dens,nval,s,d)
      dens(1)=d
      solar(1)=s
!
!  Integrate to get cumulative distribution
!
   15 freq(1)=0.0
      sbar=0.0
      do 20 i=2,nval
      del=0.5*(solar(i)-solar(i-1))*(dens(i-1)+dens(i))
      freq(i)=freq(i-1)+del
   20 sbar=sbar+solar(i)*del
      amult=1.0/freq(nval)
      sbar=sbar*amult
!
!  Normalize distribution
!
      do 30 i=1,nval
      dens(i)=dens(i)*amult
      freq(i)=freq(i)*amult
      tsol(i)=dlog(solar(i)/sbar)
   30 tden(i)=solar(i)*dens(i)
      return
      end
!
!  Subroutine to determine z-values for convolutions
!
      subroutine zval(xvec,nx,yvec,ny,zvec,nz)
      implicit none
      dimension xvec(51),yvec(51),zvec(51),ix(51)
      integer i, nx, ix, nz1, nz, ny, n1, n2, j, k, imin, ixk
      real*8 rat, crat, xvec, yvec, zvec, smin, del

      do 10 i=1,nx
   10 ix(i)=1
      nz1=nz-1
      rat=(nx*ny-1)/nz1
      crat=0.0
      zvec(1)=xvec(1)+yvec(1)
      ix(1)=2
      zvec(nz)=xvec(nx)+yvec(ny)
      n2=0
!
!  Loop through desired number of z-values
!
      do 20 i=2,nz1
      n1=n2+1
      crat=crat+rat
      n2=crat+0.5
!
!  Loop through next "rat" potential z-values in ascending order
!
      do 15 j=n1,n2
      smin=zvec(nz)+1.0
      imin=0
      do 12 k=1,nx
      if (ix(k) .gt. ny) go to 12
      ixk=ix(k)
      if (xvec(k)+yvec(ixk) .ge. smin) go to 12
      smin=xvec(k)+yvec(ixk)
      imin=k
   12 continue
   15 ix(imin)=ix(imin)+1
!
!  Fill in next actual z-value
!
   20 zvec(i)=smin
!
!  Adjust for duplicates
!
      do 30 i=2,nz
      if (zvec(i) .gt. zvec(i-1)) go to 30
      do 25 j=i,nz
      if (zvec(j) .gt. zvec(i-1)) go to 26
   25 continue
   26 del =(zvec(j)-zvec(i-1))/(j-i)
      do 27 k=i,j
   27 zvec(k)=zvec(k-1)+del
   30 continue
      return
      end
!
!  Subroutine to convolve the functions f(x) and g(y)
!
      subroutine cvolve(xvec,fofx,nx,ax,yvec,gofy,ny,ay,zvec,
     1                  fstarg,nz)
      implicit none
      dimension xvec(51),fofx(51),yvec(51),gofy(51),zvec(51),
     1          fstarg(51)
      integer i, j, nx, ny, nz, ix, iy
      real*8 xvec, yvec, zvec, fofx, fstarg, bot, top
      real*8 ex1, ax, ex2, ey1, ay, ey2, f2, f1, g2, gofy, g1, d, s
      real*8 tmp, x
!
!  Add a convenience point to f(x).
!
      xvec(nx+1)=xvec(nx)+1.0
      fofx(nx+1)=fofx(nx)
      do 50 i=1,nz
      fstarg(i)=0.0
      iy=2
      bot=yvec(1)
      ix=nx
      do 1 j=1,nx
      if (zvec(i)-xvec(ix) .gt. yvec(1)) go to 20
    1 ix=ix-1
!
!  G-domain lies entirely to the right of inverted f-domain
!  integral must be zero
!
      go to 50
!
!  Integrate over the overlapping parts of f- and g-domains
!
   20 top=dmin1(yvec(iy),zvec(i)-xvec(ix))
      ex1=dexp(ax*xvec(ix))
      ex2=dexp(ax*xvec(ix+1))
      ey1=dexp(ay*yvec(iy))
      ey2=dexp(ay*yvec(iy-1))
      f2=(fofx(ix+1)-fofx(ix))/(ex2-ex1)
      f1=fofx(ix)-f2*ex1
      g2=(gofy(iy-1)-gofy(iy))/(ey2-ey1)
      g1=gofy(iy)-g2*ey1
      d=ay-ax
      s=f1*g1*(top-bot)+f1*g2*(dexp(ay*top)-dexp(ay*bot))/ay
     2 -f2*g1*(dexp(ax*(zvec(i)-top))-dexp(ax*(zvec(i)-bot)))/ax
     3 +f2*g2*dexp(ax*zvec(i))*(dexp(d*top)-dexp(d*bot))/d
      fstarg(i)=fstarg(i)+s
!
!  Update intervals and stop if x(1) or y(ny) has been reached
!
      if (top .ge. yvec(iy)) iy=iy+1
!
      tmp = zvec(i) - xvec(ix) - 1.0d-60
      if (top .ge. tmp) ix=ix-1
      if (iy .gt. ny) go to 50
      if (ix .lt. 1) go to 50
      bot=top
      go to 20
   50 continue
      return
      end
!
!  Subroutine to perform linear interpolations
!
      subroutine interp(xvec,fofx,n,x,f)
!
!  Performs a single linear interpolation
!
      implicit none
      dimension xvec(51),fofx(51),fofxm(51,30),fm(30)
      integer i, n
      real*8 x, xvec, f, fofx, fofxm, alam, fm
!
!  Check whether x is too low or too high
!
      if (x .gt. xvec(1)) go to 20
      f=fofx(1)
      return
   20 if (x .lt. xvec(n)) go to 40
      f=fofx(n)
      return
   40 do 60 i=2,n
      if (x .le. xvec(i)) go to 80
   60 continue
   80 alam=(x-xvec(i-1))/(xvec(i)-xvec(i-1))
      f=alam*fofx(i)+(1.0-alam)*fofx(i-1)
      return
      end
!
!  Entry for multiple value functions
!
      subroutine interm(xvec,fofxm,n,nuecog,x,fm)
      implicit none
      dimension xvec(51),fofx(51),fofxm(51,30),fm(30)
      integer i, j, n, nuecog
      real*8 x, xvec, fm, fofxm, alam, fofx
!
!  Check whether x is too low or too high
!
      if (x .gt. xvec(1)) go to 120
      do 110 j=1,nuecog
  110 fm(j)=fofxm(1,j)
      return
  120 if (x .lt. xvec(n)) go to 140
      do 130 j=1,nuecog
  130 fm(j)=fofxm(n,j)
      return
  140 do 160 i=2,n
      if (x .le. xvec(i)) go to 180
  160 continue
  180 alam=(x-xvec(i-1))/(xvec(i)-xvec(i-1))
      do 190 j=1,nuecog
  190 fm(j)=alam*fofxm(i,j)+(1.0-alam)*fofxm(i-1,j)
      return
      end

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
The http://www.yolotaxpayers.com/home.htm Yolo County Taxpayers Association was featured heavily in antimeasure November 2006 Election/Measure H H and November 2006 Election/Measure I I mailings and TV ads during the November 2006 Election. Some of their contributions to the campaign are enshrined and PG&E/History preserved.

They meet the Third Wednesday of each month at 7:30 PM in the Leake Room of the Woodland Public Library.

They oppose:
  Increased Property Taxes Disguised as Assessment District, Fees or Charges
  City and County Budget Excesses
  Tax Increase Where a NO Vote Required Writing a Letter of Discontent
  County Supervisors Voting Their Own Pay Raises WITHOUT voter approval
  Inappropriate Flood Control Planning Where Water, One of Our Most Precious Natural Resources is Regarded Instead as a Dangerous Nuisance.
  Government acquisition of private property by eminent domain.


text:
! Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
! See https://llvm.org/LICENSE.txt for license information.
! SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
!
! Test for IO with the whole implied-shape array.

program test
   implicit none
   integer :: i, length
   character(len = 6) :: str
   character(len = 12) :: str2
   character(len = 40) :: str3
   integer, parameter :: arr1(1:*) = (/(i, i = 1, 6)/)
   character, parameter :: arr2(1:*) = (/'a', 'b', 'c', 'd', 'e', 'f'/)
   integer, parameter :: arr3(1:*) = [10, 11, 12, 13, 14, 15]
   character, parameter :: arr4(1:*) = ['A', 'B', 'C', 'D', 'E', 'F']
   integer, parameter :: arr5(1:*,1:*) = reshape((/(i, i = 1, 20)/), (/ 4, 5 /))

   write(str, 10) arr1
   if (str .ne. '123456') stop 1
   write(str, 20) arr2
   if (str .ne. 'abcdef') stop 2
   write(str2, 30) arr3
   if (str2 .ne. '101112131415') stop 3
   write(str, 20) arr4
   if (str .ne. 'ABCDEF') stop 4
   write(str3, 40) arr5
   if (str3 .ne. ' 1 2 3 4 5 6 7 8 91011121314151617181920') stop 5

   inquire(IOLENGTH = length) arr1
   if (length .ne. 24) stop 6
   inquire(IOLENGTH = length) arr2
   if (length .ne. 6) stop 7
   inquire(IOLENGTH = length) arr3
   if (length .ne. 24) stop 8
   inquire(IOLENGTH = length) arr4
   if (length .ne. 6) stop 9
   inquire(IOLENGTH = length) arr5
   if (length .ne. 80) stop 10

   print *, 'PASS'
10 FORMAT (6I1)
20 FORMAT (6A1)
30 FORMAT (6I2)
40 FORMAT (20I2)
end

max_stars_repo_name: abrahamtovarmob/flang

max_issues_repo_name: abrahamtovarmob/flang

max_forks_repo_name: abrahamtovarmob/flang


! ##################################################################################################################################
! Begin MIT license text.
! _______________________________________________________________________________________________________

! Copyright 2019 Dr William R Case, Jr (dbcase29@gmail.com)

! Permission is hereby granted, free of charge, to any person obtaining a copy of this software and
! associated documentation files (the "Software"), to deal in the Software without restriction, including
! without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to
! the following conditions:

! The above copyright notice and this permission notice shall be included in all copies or substantial
! portions of the Software and documentation.

! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
! OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
! THE SOFTWARE.
! _______________________________________________________________________________________________________

! End MIT license text.

      SUBROUTINE CALC_ELEM_NODE_FORCES

! Calculates elem nodal forces in local elem coord system for one elem and one subcase for all element types.

      USE PENTIUM_II_KIND, ONLY       :  BYTE, LONG, DOUBLE
      USE IOUNT1, ONLY                :  WRT_ERR, WRT_LOG, ERR, F04, F06
      USE SCONTR, ONLY                :  BLNK_SUB_NAM, INT_SC_NUM, JTSUB, NCORD, NGRID, WARN_ERR
      USE TIMDAT, ONLY                :  TSEC
      USE SUBR_BEGEND_LEVELS, ONLY    :  CALC_ELEM_NODE_FORCES_BEGEND
      USE CONSTANTS_1, ONLY           :  ZERO
      USE DEBUG_PARAMETERS
      USE MODEL_STUF, ONLY            :  AGRID, BGRID, CORD, EID, ELAS_COMP, ELDOF, ELGP, GRID, KE, KEG, KEO_BUSH,                &
                                         PEB, PE_GA_GB, PEG, PEL, PTE, RCORD, SCNUM, SUBLOD, TYPE, UEB, UEG, UEL, TE, TE_GA_GB

      USE CALC_ELEM_NODE_FORCES_USE_IFs

      IMPLICIT NONE

      CHARACTER(LEN=LEN(BLNK_SUB_NAM)):: SUBR_NAME = 'CALC_ELEM_NODE_FORCES'
      CHARACTER( 1*BYTE)              :: TEMP_OPT(6)       ! Array of EMG option indicators

      INTEGER(LONG)                   :: ACID_G            ! Actual coordinate system ID
      INTEGER(LONG)                   :: I,J,K             ! DO loop indices
      INTEGER(LONG)                   :: ICID              ! Internal coord sys no. corresponding to an actual coord sys no.
      INTEGER(LONG)                   :: I1,I2             ! Calculated displ component no's for ELAS elems
      INTEGER(LONG)                   :: NCOLS             ! Number of rows in element stiffness matrix
      INTEGER(LONG)                   :: NROWS             ! Number of cols in element stiffness matrix
      INTEGER(LONG)                   :: NUM_COMPS_GRID_1  ! No. displ components for 1st grid on ELAS elems
      INTEGER(LONG), PARAMETER        :: SUBR_BEGEND = CALC_ELEM_NODE_FORCES_BEGEND

      REAL(DOUBLE)                    :: DUM1(3),DUM2(3)   ! Intermediate variables
      REAL(DOUBLE)                    :: PHID, THETAD      ! Outputs from subr GEN_T0L
      REAL(DOUBLE)                    :: T0G(3,3)          ! Matrix to transform offsets from global to basic  coords
      REAL(DOUBLE)                    :: TET(3,3)          ! Transpose of TE
      REAL(DOUBLE)                    :: TET_GA_GB(3,3)    ! Transpose of TE
      REAL(DOUBLE)                    :: TR(12,12)         ! Matrix with 4 TE matrices on the diagonal

! **********************************************************************************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9001) SUBR_NAME,TSEC
 9001    FORMAT(1X,A,' BEGN ',F10.3)
      ENDIF

! **********************************************************************************************************************************
      NROWS = ELDOF
      NCOLS = ELDOF

      DO I=1,NCOLS
         PEL(I) = ZERO
      ENDDO

! **********************************************************************************************************************************
! Calc forces for one element. The ELAS and ROD1 elem have very sparse stiffness matrices, so an explicit form is used
! for them. All other element forces are calculated by multiplication of complete stiffness matrix with the displ's.

      IF (TYPE(1:4) == 'ELAS') THEN                        ! Calculate forces for ELAS1-4 elems

         I1 = ELAS_COMP(1)
         CALL GET_GRID_NUM_COMPS ( AGRID(1), NUM_COMPS_GRID_1, SUBR_NAME )
         I2 = NUM_COMPS_GRID_1 + ELAS_COMP(2)
         PEL(I1) = KE(I1,I1)*UEL(I1) + KE(I1,I2)*UEL(I2)   ! Note: KE is global and local for the ELAS elems
         PEL(I2) = KE(I2,I1)*UEL(I1) + KE(I2,I2)*UEL(I2)

      ELSE IF (TYPE == 'ROD     ') THEN                    ! Calculate forces for ROD1 elem

         IF (SUBLOD(INT_SC_NUM,2) > 0) THEN
            PEL(1) = -PTE(1,JTSUB)
            PEL(7) = -PTE(7,JTSUB)
         ENDIF

         PEL( 1) = PEL(1) + KE( 1, 1)*UEL( 1) + KE( 1, 7)*UEL( 7)
         PEL( 4) =          KE( 4, 4)*UEL( 4) + KE( 4,10)*UEL(10)
         PEL( 7) = PEL(7) + KE( 7, 1)*UEL( 1) + KE( 7, 7)*UEL( 7)
         PEL(10) =          KE(10, 4)*UEL( 4) + KE(10,10)*UEL(10)

      ELSE IF (TYPE == 'USERIN  ') THEN

         WRITE(F06,9991) TYPE

      ELSE IF (TYPE == 'BUSH    ') THEN

         CALL ELEM_TRANSFORM_LBG ( 'KE', KE, PTE )

         DO I=1,6
            TEMP_OPT(I) ='N'
         ENDDO
         TEMP_OPT(4) = 'Y'
         CALL ELMOFF ( TEMP_OPT, 'Y' )

         DO I=1,NROWS
            DO J=1,NCOLS
               PEG(I) = PEG(I) + KEG(I,J)*UEG(J)
            ENDDO
         ENDDO

         DO I=1,2
            ACID_G = GRID(BGRID(I),3)                      ! Get global coord sys for this grid
            IF (ACID_G /= 0) THEN                          ! Global is not basic so need to transform offset from basic to global
               ICID = 0
               DO J=1,NCORD
                  IF (ACID_G == CORD(J,2)) THEN
                     ICID = J
                     EXIT
                  ENDIF
               ENDDO
               CALL GEN_T0L ( BGRID(I), ICID, THETAD, PHID, T0G )
               IF (I == 1) THEN
                  DO J=1,3
                     PEB(J)   = T0G(J,1)*PEG(J)   + T0G(J,2)*PEG(J)   + T0G(J,3)*PEG(J)
                  ENDDO
                  DO J=1,3
                     PEB(J+3) = T0G(J,1)*PEG(J+3) + T0G(J,2)*PEG(J+3) + T0G(J,3)*PEG(J+3)
                  ENDDO
               ELSE
                  DO J=1,3
                     PEB(J+6) = T0G(J,1)*PEG(J+6) + T0G(J,2)*PEG(J+6) + T0G(J,3)*PEG(J+6)
                  ENDDO
                  DO J=1,3
                     PEB(J+9) = T0G(J,1)*PEG(J+9) + T0G(J,2)*PEG(J+9) + T0G(J,3)*PEG(J+9)
                  ENDDO
               ENDIF
            ELSE                                           ! Global was basic so no transformation of coords needed
               DO J=1,12
                  PEB(J) = PEG(J)
               ENDDO
            ENDIF
         ENDDO

         DO I=1,3
            DO J=1,3
               TET(I,J) = TE(J,I)
            ENDDO
         ENDDO

         DO I=1,12
            DO J=1,12
               TR(I,J) = ZERO
            ENDDO
         ENDDO

         DO I=1,3                                          ! TR is a 12x12 matrix with 4 TET matrices on its diagonal
            DO J=1,3
               TR(I  ,J  ) = TET(I,J)
               TR(I+3,J+3) = TET(I,J)
               TR(I+6,J+6) = TET(I,J)
               TR(I+9,J+9) = TET(I,J)
            ENDDO
         ENDDO

         CALL MATMULT_FFF_T ( TR, PEB, 12, 12, 1, PEL )
         DO I=1,3
            DO J=1,3
               TET_GA_GB(I,J) = TE_GA_GB(J,I)
            ENDDO
         ENDDO

         DO I=1,12
            DO J=1,12
               TR(I,J) = ZERO
            ENDDO
         ENDDO

         DO I=1,3                                          ! TR is a 12x12 matrix with 4 TET matrices on its diagonal
            DO J=1,3
               TR(I  ,J  ) = TET_GA_GB(I,J)
               TR(I+3,J+3) = TET_GA_GB(I,J)
               TR(I+6,J+6) = TET_GA_GB(I,J)
               TR(I+9,J+9) = TET_GA_GB(I,J)
            ENDDO
         ENDDO

         CALL MATMULT_FFF_T ( TR, PEB, 12, 12, 1, PE_GA_GB )

      ELSE                                                 ! Calculate forces for any other type of elem

         DO I=1,NROWS

            IF (SUBLOD(INT_SC_NUM,2) > 0) THEN
               PEL(I) = -PTE(I,JTSUB)
            ENDIF

            DO J=1,NCOLS
               PEL(I) = PEL(I) + KE(I,J)*UEL(J)
            ENDDO
         ENDDO

      ENDIF

      IF (DEBUG(56) > 0) THEN                              ! Print UEL, PEL to f06 for debug
         WRITE(F06,5000)
         WRITE(F06,5004) TRIM(TYPE), EID, SCNUM(INT_SC_NUM)
         J = 0 ; K = 0
         DO I=1,ELDOF
            J = J+1
            IF (J == 1) THEN
               K = K+1
               WRITE(F06,5006) '  GRID, COMP, UEL, PEL = ', AGRID(K), J, UEL(I), PEL(I)
            ELSE
               WRITE(F06,5007) '        COMP, UEL, PEL = ', J, UEL(I), PEL(I)
            ENDIF
            IF (J == 6) THEN
               WRITE(F06,*)
               J = 0
            ENDIF
         ENDDO
         WRITE(F06,5000)
      ENDIF

! **********************************************************************************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9002) SUBR_NAME,TSEC
 9002    FORMAT(1X,A,' END  ',F10.3)
      ENDIF

      RETURN

! **********************************************************************************************************************************

 5000 FORMAT('=============================================================================================================',/)

 5004 FORMAT(16X, 'S U B R O U T I N E   CALC_ELEM_NODE_FORCES   F O R   ',A,I8,', S/C',I8,/,                                      &

             28X,'(displacements and node forces in element coordinates)',/)

 5006 FORMAT(21X,A,I8,1X,I3,2ES15.6)

 5007 FORMAT(21X,A,9X,I3,2ES15.6)

 9991 FORMAT(' *INFORMATION: ELEMENT NODE FORCE CALCULATION NOT PROGRAMMED FOR ',A,' ELEMENTS')

! **********************************************************************************************************************************

      END SUBROUTINE CALC_ELEM_NODE_FORCES


text:
      SUBROUTINE renum2(adjnt,start,level,lpntn,npoin,ns)
!***********************************************************************
!*
!****  computes a set of posibles starting npoin
!*
!***********************************************************************
      IMPLICIT none
!     routine parameters
      INTEGER (kind=4),INTENT(IN) :: adjnt(1:),npoin
      INTEGER (kind=4),INTENT(OUT):: lpntn(:),level(1:),start(1:),ns
!     local variables
      LOGICAL (kind=4) :: better
      INTEGER (kind=4) :: i,iroot,depth,width,nnoded,ndpth,nwdth,lhw

      INTERFACE
        INCLUDE 'renum4.h'
      END INTERFACE

!.... begin iteration
!.... select intial root node arbitrarily and generate its level
!.... structure

      iroot = 1
      DO

        CALL renum4(start,level,depth,adjnt,width,npoin,iroot,lhw)

        !..  create a list of npoin which are at maximum distance from root
        !..  node and store the root
        ns = iroot

        !..  loop over npoin at maximum distance from root node
        !..  generate level structure for each node
        !..  set switch IF a level structure of greater depth occurs
        better = .FALSE.

        DO i = 1,lhw
          nnoded = start(i)
          CALL renum4(lpntn,level,ndpth,adjnt,nwdth,npoin,nnoded,lhw)
          IF(ndpth >= depth) THEN
            IF((ndpth /= depth) .OR. (nwdth < width)) THEN
              iroot = nnoded
              depth = ndpth
              width = nwdth
              better = .TRUE.
            END IF
          END IF
        END DO
        IF (.NOT.better) EXIT
      END DO

      RETURN

      END SUBROUTINE renum2

max_stars_repo_name: jerebenitez/IFE-simpact-openfoam

max_issues_repo_name: jerebenitez/IFE-simpact-openfoam

max_forks_repo_name: jerebenitez/IFE-simpact-openfoam


text:
\ hf.f -- hiForth "protolude", implementing non-microcoded ops
\ Copyright (C) 2017-2019 Austin Seipp. See Copyright Notice in COPYING

\ Commentary:

  \ NOTE: See the corresponding hf.S for more information, including changes.
  \
  \ TODO FIXME: Lorem ipsum...

\ ------------------------------------------------------------------------------

: TRUE  1 ;
: FALSE 0 ;

: TRUE?  ( bool -- bool ) TRUE = ;
: FALSE? ( bool -- bool ) FALSE = ;

: NOT ( bool -- bool ) 0 = ;
: <> ( x y -- bool ) = NOT ;

: B>C ( bool -- char )
  0BRANCH 20
  84
  BRANCH 12
  70
  ;

: '\n' 10 ;

\ TODO FIXME: use this implementation?
: >CFA2 ( a -- codeword )
  4+ 4+        \ skip length chain
  1+           \ skip flags byte
  DUP @        \ grab length         STACK: addr len
  SWAP 1+      \ skip length byte    STACK: len addr
  +            \ add                 STACK: addr

  3 +          \ round to 4-byte boundary
  3 INVERT AND
  ;

: DFA> ( codeword -- body ) >CFA 4+ ;

: FIND ( a len -- w ) ;

: PROTOMAIN ( -- result )
  TRUE FALSE <>   \ dumb test

  0BRANCH 16      \ conditional jump
  TRUE
  BRANCH 8        \ unconditional
  FALSE

  B>C EMIT '\n' EMIT \ spit out result
  0 \ result
  ;

\ ------------------------------------------------------------------------------
\ -- El fin --------------------------------------------------------------------

\ Local Variables:
\ mode: forth
\ fill-column: 80
\ indent-tabs-mode: nil
\ c-basic-offset: 2
\ buffer-file-coding-system: utf-8-unix
\ End:

\ ------------------------------- hf.f ends here -------------------------------


text:
C
C       $Id: vvmpxy.f,v 1.13 2008-07-27 00:17:35 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
C
C-----------------------------------------------------------------------
C
      SUBROUTINE VVMPXY (X,Y,U,V,UVM,XB,YB,XE,YE,IST)
C
C This routine incorporates the functions previously known as FX, FX,
C MXF, and MYF into a single subroutine. However, a compatibility
C mode is retained so that old codes that depend on custom versions of
C these functions will still work. The recommended approach for new
C code that requires a custom transformation routine is to write your
C own version of VVUMXY (Velocity Vector -- User Map X,Y). Whenever
C IMAP has a value other than 0,1, or 2 VVUMXY will be called.
C
C The most important thing to note is that the input vector position
C (a single point -- input parameters X and Y) is in user coordinate
C space, while the output (two points defining the beginning and
C end of the vector, XB,YB and XE,YE) is in the fractional coordinate
C space.
C
C Note that each vector consists of three potentially tranformable
C components: position, direction, and magnitude. Depending on the
C user's requirement, any or all of these components may need to be
C transformed. As an example the map transformations provided by
C VVMPXY transform the position and direction of the vector, but keep
C the vector magnitude based on a uniform coordinate system. This
C is the motivation for having the input vector position in user
C coordinates, but the output in the fractional coordinate system.
C The fractional system is known to be uniform -- a unit in X always
C equals a unit in Y.
C
C Input parameters:
C
C X,Y   -- vector position in the user coordinate system
C U,V   -- vector components from the U,V arrays for this position
C UVM   -- vector magnitude,  (SQRT(U*U+V*V))
C
C Output parameters:
C
C XB,YB -- starting point of the vector in the fractional
C          coordinate system
C XE,YE -- ending point of the vector in the fractional
C          coordinate system
C IST   -- status results of the mapping: 0 indicates success
C
C
C The mapping common block: made available to user mapping routines
C
      COMMON /VVMAP/
     +                IMAP       ,
     +                XVPL       ,XVPR       ,YVPB       ,YVPT       ,
     +                WXMN       ,WXMX       ,WYMN       ,WYMX       ,
     +                XLOV       ,XHIV       ,YLOV       ,YHIV       ,
     +                SXDC       ,SYDC       ,NXCT       ,NYCT       ,
     +                RLEN       ,LNLG       ,INVX       ,INVY       ,
     +                ITRT       ,IWCT       ,FW2W       ,FH2H       ,
     +                DVMN       ,DVMX       ,RBIG       ,IBIG
C
      SAVE /VVMAP/
C
C Math constants
C
      PARAMETER (PDTOR  = 0.017453292519943,
     +           PRTOD  = 57.2957795130823,
     +           P1XPI  = 3.14159265358979,
     +           P2XPI  = 6.28318530717959,
     +           P1D2PI = 1.57079632679489,
     +           P5D2PI = 7.85398163397448)
C
C --------------------------------------------------------------------
C
C --------------------------------------------------------------------
C
C Local parameters:
C
C PRCFAC - Precision factor used to resolve float equality within
C            the precision of a 4 byte REAL
C PVFRAC - Initial fraction of the vector magnitude used to
C            determine the differential increment
C PFOVFL - Floating point overflow value
C IPMXCT - Number of times to allow the differential to increase
C PDUVML - Multiplier when the differential is increased
C PCSTST - Test value for closeness to 90 degree latitude
C
      PARAMETER (PRCFAC=1E5,
     +           PVFRAC=0.001,
     +           PFOVFL=1E12,
     +           IPMXCT=40,
     +           PDUVML=2.0,
     +           IPCTST=PRCFAC*90)
C
C Local variables:
C
C DUV   - incremental UV magnitude difference
C CLT   - cosine of the latitude
C ICT   - repetition count for the difference loop
C SGN   - sign factor, for switching from pos to neg difference
C XT,YT - temporary projected X,Y position values
C XTF,YTF - temporary X,Y fractional coordinate values
C DV1,DV2 - vector magnitudes
C T     - temporary value
C IZO   - zero magnitude flag
C
C Each end of the vector is calculated in the user coordinate space.
C
C Zero the status flag
C
      IST=0
C
C Check for zero magnitude
C
      IF (INT(UVM*PRCFAC+0.5) .EQ. 0) THEN
         IZO = 1
      ELSE
         IZO = 0
      END IF
C
      If (IMAP .EQ. 0) THEN
C
C Linear transformation
C
         IF (X .LT. WXMN .OR. X .GT. WXMX .OR.
     +        Y .LT. WYMN .OR. Y .GT. WYMX) THEN
            IST=-5
            RETURN
         ENDIF
         XB=CUFX(X)
         YB=CUFY(Y)
C
C Return special status if zero magnitude
C
         IF (IZO .EQ. 1) THEN
            XE = XB
            YE = YB
            IST = -999
            RETURN
         END IF
C
         XE=XB+U*SXDC
         YE=YB+V*SYDC
C
C The following code accounts for non-uniform window coordinate
C systems, such as log scaled coordinates. The technique used
C is similar to the one used in the EZMAP case
C
         IF (ITRT .EQ. 1) THEN
C
            DV1=SQRT((XE-XB)*(XE-XB)+(YE-YB)*(YE-YB))
C
C Set up an initial increment factor
C
            DUV=0.1/UVM
            ICT=0
            SGN=1.0
C
 10         CONTINUE
C
C Bail out if it's not working
C
            IF (ICT .GT. IPMXCT) THEN
               IST = -3
               RETURN
            END IF
C
C Calculate the incremental end points, then check to see if
C they take us out of the user coordinate boundaries. If they
C do, try incrementing in the other direction
C
            XT=X+SGN*U*DUV
            YT=Y+SGN*V*DUV
            IF (XT .LT. WXMN .OR. XT .GT. WXMX .OR.
     +           YT .LT. WYMN .OR. YT .GT. WYMX) THEN
               IF (SGN .EQ. 1.0) THEN
                  SGN = -1.0
                  GO TO 10
               ELSE
                  IST=-4
                  RETURN
               ENDIF
            END IF
C
C Convert to fractional coordinates and find the incremental
C distance in the fractional system. To ensure that this distance
C is meaningful, we require that it be between 1E3 and 1E4
C times smaller than the maximum vector length.
C
            XTF=CUFX(XT)
            YTF=CUFY(YT)
            DV2=SQRT((XTF-XB)*(XTF-XB)+(YTF-YB)*(YTF-YB))
            IF (DV2*1E3 .GT. DVMX) THEN
               ICT=ICT+1
               DUV=DUV/2.0
               GO TO 10
            ELSE IF (DV2*1E4 .LT. DVMX) THEN
               ICT=ICT+1
               DUV=DUV*2.0
               GO TO 10
            END IF
C
C The actual endpoints are found using the ratio of the incremental
C distance to the actual distance times the fractional component
C length
            XE=XB+SGN*(XTF-XB)*DV1/DV2
            YE=YB+SGN*(YTF-YB)*DV1/DV2
C
         END IF
C
C EZMAP overlay
C
      ELSE IF (IMAP .EQ. 1) THEN
C
C X is longitude, Y is latitude
C If Y is 90 degrees, a directional vector is meaningless
C Since EZMAP uses uniform user coordinates, don't convert to
C fractional coordinates until the end
C The tranformation type parameter is not used.
C
         IF (INT(ABS(Y)*PRCFAC+0.5).EQ.IPCTST) THEN
            IST=-1
            RETURN
         END IF
C
C Project the starting value: bail out if outside the window
C
         CALL MAPTRA (Y,X,XB,YB)
         IF (XB .LT. WXMN .OR. XB .GT. WXMX .OR.
     +       YB .LT. WYMN .OR. YB .GT. WYMX) THEN
            IST=-5
            RETURN
         END IF
C
C Return special status if zero magnitude
C
         IF (IZO .EQ. 1) THEN
            XE = XB
            YE = YB
            IST = -999
            RETURN
         END IF
C
C The incremental distance is proportional to a small fraction
C of the vector magnitude
C
         DUV=PVFRAC/UVM
         CLT=COS(Y*PDTOR)
C
C Project the incremental distance. If the positive difference doesn't
C work, try the negative difference. If the difference results in a
C zero length vector, try a number of progressively larger increments.
C
         ICT=0
         SGN=1.0
 20      CONTINUE

         CALL MAPTRA(Y+SGN*V*DUV,X+SGN*U*DUV/CLT,XT,YT)

         DV1=SQRT((XT-XB)*(XT-XB)+(YT-YB)*(YT-YB))
         IF (DV1 .GT. RLEN) THEN
            IF (SGN .EQ. -1.0) THEN
               IST=-4
               RETURN
            ELSE
               SGN=-1.0
               GO TO 20
            END IF
         END IF
C
         IF (INT(DV1*PRCFAC) .EQ. 0) THEN
            IF (ICT .LT. IPMXCT) THEN
               ICT = ICT + 1
               DUV=DUV*PDUVML
               GO TO 20
            ELSE
               IST=-3
               RETURN
            END IF
         END IF
C
         IF (ABS(XT) .GE. PFOVFL .OR. ABS(YT) .GE. PFOVFL) THEN
            IST=-6
            RETURN
         END IF
C
         T=SGN*((XT-XB)/DV1)*UVM
         XB=CUFX(XB)
         XE=XB+T*SXDC
         T=SGN*((YT-YB)/DV1)*UVM
         YB=CUFY(YB)
         YE=YB+T*SYDC
C
      ELSE IF (IMAP .EQ. 2) THEN
C
C ... the polar coordinate case ...
C
         TH=PDTOR*Y
         XB=CUFX(X*COS(TH))
         YB=CUFY(X*SIN(TH))
C
C Return special status if zero magnitude
C
         IF (IZO .EQ. 1) THEN
            XE = XB
            YE = YB
            IST = -999
            RETURN
         END IF
C
         XE=XB+U*SXDC
         YE=YB+V*SYDC
C
         IF (ITRT.GE.1) THEN
            DV1=SQRT((XE-XB)*(XE-XB)+(YE-YB)*(YE-YB))
            PHI=ATAN2(V,U)
            ANG=TH+PHI
            XE=XB + DV1*COS(ANG)
            YE=YB + DV1*SIN(ANG)
         END IF
C
C ... and everything else.
C
      ELSE
C
         CALL VVUMXY(X,Y,U,V,UVM,XB,YB,XE,YE,IST)
C
      END IF
C
C Done.
C
      RETURN
C
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
!----- LGPL --------------------------------------------------------------------
!
!  Copyright (C)  Stichting Deltares, 2011-2016.
!
!  This library is free software; you can redistribute it and/or
!  modify it under the terms of the GNU Lesser General Public
!  License as published by the Free Software Foundation version 2.1.
!
!  This library is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!  Lesser General Public License for more details.
!
!  You should have received a copy of the GNU Lesser General Public
!  License along with this library; if not, see <http://www.gnu.org/licenses/>.
!
!  contact: delft3d.support@deltares.nl
!  Stichting Deltares
!  P.O. Box 177
!  2600 MH Delft, The Netherlands
!
!  All indications and logos of, and references to, "Delft3D" and "Deltares"
!  are registered trademarks of Stichting Deltares, and remain the property of
!  Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
!  $Id: test_01.f90 5717 2016-01-12 11:35:24Z mourits $
!  $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/6686/src/utils_lgpl/deltares_common/tests/01_getenv_test/test_01.f90 $
subroutine test_01
!
!   Program to test some functions which are available in the util_mf libarray
!
!   util_mf: UTILities for Mor and Flow
!
!   Jan.Mooiman@deltares.nl
!   22 sep 05
!
    implicit none
    logical            :: error
    character(len=256) :: value
    character(len=20)  :: env (2)

    call util_getenv('D3D_HOME', value)
    write(*,*) 'Result D3D_HOME: <',trim(value),'>'

    env(1) = 'D3D_HOME'
    env(2) = 'd3d_home'
    value  = ' '
    call util_getenv(env(1), value)
    write(*,*) 'Result D3D_HOME: <',trim(value),'>'

!    call getmp(error, value)
!    write(*,*) 'Result getmp : ',error,'<',trim(value),'>'

end subroutine test_01


text:
C
C $Id: nnuser_f.f,v 1.5 2008-07-27 03:10:13 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE NNSETC(PNAM,CVAL)
C
      CHARACTER*(*) PNAM,CVAL
C
      ILEN = LEN(CVAL)
      CALL FNNSETC(PNAM,CVAL,ILEN)
C
      RETURN
      END
      SUBROUTINE NNGETC(PNAM,CVAL)
C
      CHARACTER*(*) PNAM,CVAL
C
      ILEN = LEN(CVAL)
      CALL FNNGETC(PNAM,CVAL,ILEN)
C
      RETURN
      END