text: Alison Stevenson is currently a student at UC Davis. She is a third year English (emphasis on creative writing) and Film Studies double major, with a minor in African American Studies. She has acted in the UC Davis Department of Theater & Dances production of Ghost, Bathtub, Windmill written by Elise Kane, as well as Studio 301s productions, Comedy On The Rocks With A Twist and A Midsummer Nights Dream. Currently she is playing the role of Gertie in the play Fuddy Meers by David LindsayAbaire which will open on February 17, 2010. Last year she directed the full length comedy, Beyond Therapy by Christopher Durang and is also a board member for Studio 301 in charge of publicity. Aside from theater, Alison is also a member of Birdstrike Theater, the improv comedy troupe at Davis as well as a member of the mixed media board for Nameless Magazine and has some of her photography on the website. Her latest business venture is a tshirt website through cafepress that her and one of her good friends, Elizabeth Vazquez, created. The site is http://www.cafepress.com/alizabeth.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!$omp target map(tofrom: y) map(to:x) map(from:is_cpu)
is_cpu=omp_is_initial_device();
do i=1,ARRAY_SIZE
        y(i) = a*x(i) + y(i)
end do
!$omp end target

max_forks_repo_name: junxnone/oneAPI-samples


text:
! { dg-do run }
! { dg-options "-std=gnu" }
character(25) :: sround, ssign, sasynchronous, sdecimal, sencoding
integer :: vsize, vid
logical :: vpending

open(10, file='mydata', asynchronous="yes", blank="null", &
& decimal="comma", encoding="utf-8", sign="plus")

inquire(unit=10, round=sround, sign=ssign, size=vsize, id=vid, &
& pending=vpending, asynchronous=sasynchronous, decimal=sdecimal, &
& encoding=sencoding)

if (ssign.ne."PLUS") call abort
if (sasynchronous.ne."YES") call abort
if (sdecimal.ne."COMMA") call abort
if (sencoding.ne."UTF-8") call abort
if (vpending) call abort

close(10, status="delete")
end

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
subroutine  corrente

!****************************************************************!
! COMPUTE: current density from J = curl(B)                      !
!          [2D version]                                          !
! -------------------------------------------------------------- !
! original version:             F. Califano, 2006                !
! MPI parallel version:         M. Faganello/F. Valentini, 2008  !
! Anisotropic/FLR-MHD version:  S. S. Cerri, 2011                !
!****************************************************************!

use parameter_mod
use parallel_mod
use dom_distr_mod
use box_mod
use deriv_mod
use fields_UJ_mod
use fields_EB_mod

IMPLICIT NONE

integer :: ix, iy, iz

REAL(dp),ALLOCATABLE :: zx(:), zy(:), zz(:)
REAL(dp), ALLOCATABLE :: Bt(:,:,:)

allocate(zx(nx))
allocate(zy(ny))
allocate(zz(nz))

ALLOCATE( Bt( nx, nyt, nz ) )  ! nyt ~ ny/nproc

!--compute Jx
do iz = 1, nz
  do ix = 1, nxl
    call  dery_1( Bz(ix,:,iz), Jx(ix,:,iz) )
  enddo
enddo

!--compute Jy
call traspdist( Bz, Bt, 1 )
do iz = 1, nz
  do iy = 1, nyt
    call  derx_1(Bt(:,iy,iz), zx )
    Bt(:,iy,iz) = - zx
  enddo
enddo
call traspdist( Jy, Bt, -1 )

!--compute Jz
call traspdist( By, Bt, 1 )
do iz = 1, nz
  do iy = 1, nyt
    call  derx_1( Bt(:,iy,iz), zx )
    Bt(:,iy,iz) = zx
  enddo
enddo
call traspdist( Jz, Bt, -1 )
!
do iz = 1, nz
  do ix = 1, nxl
    call dery_1( Bx(ix,:,iz), zy )
    Jz(ix,:,iz) = Jz(ix,:,iz) - zy
  enddo
enddo

deallocate(zx)
deallocate(zy)
deallocate(zz)
deallocate(Bt)

end subroutine

max_stars_repo_name: sscerr/eTF

max_issues_repo_name: sscerr/eTF

max_forks_repo_name: sscerr/eTF


text:
       subroutine timeDependBCElas(x, iBC, BC, BC_flow, umeshold)
c
c-----------------------------------------------------------------
c
c-----------------------------------------------------------------
c
      include "common.h"
      include "mpif.h"

      real*8    x(numnp,nsd)
      real*8    umeshold(numnp,nsd)
      dimension iBC(nshg),        BC(nshg,3), BC_flow(nshg,3)
c
      if (numrbs .gt. 0) then
        call rigidBodyBCElas(x, iBC, BC, BC_flow)
      endif
c
      if (elasFDC .gt. 0) then
        call prescribedBCElas(x, iBC, BC, BC_flow, umeshold)
      endif
c
      return
      end
c
c-----------------------------------------------------------------
c
c-----------------------------------------------------------------
c
      subroutine rigidBodyBCElas(x,  iBC,  BC,  BC_flow)
c
        use rigidBodyReadData
        use rigidBodyForce
c
        include "common.h"
        include "mpif.h"
c
        real*8    x(numnp,nsd)
c        real*8    umeshold(numnp,nsd)
        dimension iBC(nshg), BC(nshg,3), BC_flow(nshg,3)
        real*8    t
        real*8, dimension(nsd) :: ra, rc, r_disp
c
        call calc_rbMotion
c
c.... loop over mesh vertices
        do i = 1,numnp
          if ( (ibits(iBC(i),14,3) .eq. 7) .and.
     &         (rbFlags(i) .gt. 0) ) then
c.... calculate the displacement due to rotation
            t = rbAng(rbFlags(i))/180.0*pi
            ra = rb_prop(rbFlags(i),5:7)
            rc = rb_prop(rbFlags(i),8:10)
            r_disp(1) =
     &               (cos(t)+ra(1)*ra(1)*(1-cos(t))-1)*(x(i,1)-rc(1))
     &         + (ra(1)*ra(2)*(1-cos(t))-ra(3)*sin(t))*(x(i,2)-rc(2))
     &         + (ra(1)*ra(3)*(1-cos(t))+ra(2)*sin(t))*(x(i,3)-rc(3))
            r_disp(2) =
     &           (ra(2)*ra(1)*(1-cos(t))+ra(3)*sin(t))*(x(i,1)-rc(1))
     &             + (cos(t)+ra(2)*ra(2)*(1-cos(t))-1)*(x(i,2)-rc(2))
     &         + (ra(2)*ra(3)*(1-cos(t))-ra(1)*sin(t))*(x(i,3)-rc(3))
            r_disp(3) =
     &           (ra(3)*ra(1)*(1-cos(t))-ra(2)*sin(t))*(x(i,1)-rc(1))
     &         + (ra(3)*ra(2)*(1-cos(t))+ra(1)*sin(t))*(x(i,2)-rc(2))
     &             + (cos(t)+ra(3)*ra(3)*(1-cos(t))-1)*(x(i,3)-rc(3))
c.... set corresponding mesh elas BC
            BC(i,1:3) = rbDisp(rbFlags(i),1:3) + r_disp(1:3)
c.... update flow BC
            BC_flow(i,1:3) = rbVel(rbFlags(i),1:3) + r_disp(1:3)/Delt(1)
          endif
        enddo ! end loop numnp
c
      return
      end
c
c-----------------------------------------------------------------
c
c-----------------------------------------------------------------
c
       subroutine prescribedBCElas(x, iBC, BC, BC_flow, umeshold)
c
c-----------------------------------------------------------------
c
c-----------------------------------------------------------------
c
      use m2gfields ! read m2g fields
      use core_snap
      use interfaceflag
      use core_rigid_body
c
      include "common.h"
      include "mpif.h"
c
c
      real*8    x(numnp,nsd)
      real*8    disp(numnp,nsd)
      real*8    umeshold(numnp,nsd)
      dimension iBC(nshg),        BC(nshg,3), BC_flow(nshg,3)
      integer   casenumber
      real*8    acc
      real*8    totalForce(3),    objMass
      dimension Forin(3),         Forout(3)
c
c.... dynamic origin, x translation, rotation frequence
      real*8    dyn_org,   xtsl,  rotf
      real*8    dyn_lnt,   shrk,  shrkfactor
      integer   answer
      real*8    cent(numrbs, 3)
c
      casenumber = 0
      if (elasFDC .gt. 0) then
        if (myrank .eq. master) then
          write(*,*) "use Force-driven case:", elasFDC
        endif
        casenumber = elasFDC
      endif
c
c.... Update BC value based on geom and iBC
c
      if ( casenumber .eq. 1 ) then
        acc = 60000.0
        do i = 1,numnp
          if ( (ibits(iBC(i),14,3) .eq. 7) .and.
     &         (abs(x(i,1)-0.8) .lt. 0.79) .and.
     &         ((x(i,2)*x(i,2) + x(i,3)*x(i,3)) .lt. 0.004225)) then
c.... 0.004225 = 0.065^2
            BC(i,1)   = acc * lstep * Delt(1) * Delt(1)
            BC(i,2)   = zero
            BC(i,3)   = zero
          endif ! end if inside box
        enddo ! end loop numnp
      endif  ! end case 1
c
c.... Update BC value based on total force on the object
c
      if ( casenumber .eq. 2 ) then
c.... debugging {
c        do j = 1,numrbs
c          call core_get_centroid(rbsTags(j), cent(j,:))
c        enddo
c.... debugging }
        totalForce(:) = zero
        objMass = 15.0 ! kg
c.... bcast force to all processors
        if (numpe > 1) then
          do j = 1, nsrfCM
            Forin  = (/ Force(1,j), Force(2,j), Force(3,j) /)
            call MPI_BCAST (Forin(1), 3, MPI_DOUBLE_PRECISION,
     &                      master,      MPI_COMM_WORLD,ierr)
            Force(1:3,j) = Forin(1:3)
          enddo
        endif
c.... collect total force
        do j = 1,nsrfCM
          totalForce(:) = totalForce(:) + Force(:,j)
        enddo ! end collect total force

        acc = totalForce(1) / objMass
        if (myrank .eq. master) then
          write(*,*) "projectile acceleration is:", acc
        endif

        do i = 1,numnp
          if ( (ibits(iBC(i),14,3) .eq. 7) .and.
     &         (abs(x(i,1)-0.8) .lt. 0.79) .and.
     &         ((x(i,2)*x(i,2) + x(i,3)*x(i,3)) .lt. 0.004225)) then
c.... 0.004225 = 0.065^2
            BC(i,1)   = (umeshold(i,1) + acc * Delt(1)) * Delt(1)
            BC(i,2)   = zero
            BC(i,3)   = zero

c.... >>> hard-coding update flow velocity boundry condition on mortar surfaces
            BC_flow(i,1:3) = BC(i,1:3) / Delt(1)
c.... <<< hard-coding

          endif ! end if inside box
        enddo ! end loop numnp
      endif  ! end case 2
c
c.... Update BC value based on geom and iBC
c.... rotate + trans + shrink
c
      if ( casenumber .eq. 3 ) then
        do i = 1,numnp
c          if ( (ibits(iBC(i),14,3) .eq. 7) .and.
          if ( (ibits(iBC(i),3,3) .eq. 7) .and.
     &         (abs(x(i,1)) .lt. 0.2) .and.
     &         (abs(x(i,2)) .lt. 0.2) .and.
     &         (abs(x(i,3)) .lt. 0.2) ) then
c
            disp(i,1) = -0.005 * (x(i,1)) !+ 0.02
     &                + (x(i,1)) * (cos(pi/600) - 1.0)
     &                - (x(i,2)) *  sin(pi/600)
            disp(i,2) = !-0.005 * (x(i,2)) !+ 0.01
     &                + (x(i,2)) * (cos(pi/600) - 1.0)
     &                + (x(i,1)) *  sin(pi/600)
            disp(i,3) = -0.005 * x(i,3)
c
            BC(i,:)   = disp(i,:)
          endif ! end if inside box
        enddo ! end loop numnp
      endif  ! end case 3
c
c.... Update BC value based on geom and iBC
c.... twist one end of a cylinder
c
      if ( casenumber .eq. 4 ) then
        do i = 1,numnp
          if ( (ibits(iBC(i),14,3) .eq. 7) .and.
     &         (x(i,1) .lt. 0.501) .and.
     &         (x(i,1) .gt. 0.499) ) then
c
            disp(i,1) = 0.0
            disp(i,2) = x(i,2) * (cos(pi/100) - 1.0)
     &                - x(i,3) *  sin(pi/100)
            disp(i,3) = x(i,3) * (cos(pi/100) - 1.0)
     &                + x(i,2) *  sin(pi/100)
c
            BC(i,:)   = disp(i,:)
          endif ! end if inside box
        enddo ! end loop numnp
      endif  ! end case 4
c
c.... test case 5
c.... original mov+rot+shrk
c
      if ( casenumber .eq. 5 ) then
        if (lstep .eq. 0) then
        dyn_org    = 2.625e-3
        do i = 1,numnp
          if ( (ibits(iBC(i),14,3) .eq. 7) .and.
     &         (x(i,1) .lt. 6.65e-3) .and. (x(i,1) .gt. -1.4e-3) .and.
     &         ((x(i,2)*x(i,2) + x(i,3)*x(i,3)) .lt. 2.89e-6) ) then
            if( x(i,2) .ge. 0.6e-3 ) then ! top
              if ( x(i,1) .le. (dyn_org-1.025e-3) ) then ! top tail
                disp(i,1) = -0.05 * (x(i,1)-0.5e-3) !+ 0.02
                disp(i,2) = -0.05 * (x(i,2)-1.125e-3) !+ 0.01
                disp(i,3) = -0.05 * x(i,3)
              else if ( x(i,1) .ge. (dyn_org+1.025e-3) ) then ! top head
                disp(i,1) = -0.05 * (x(i,1)-4.75e-3) !+ 0.02
                disp(i,2) = -0.05 * (x(i,2)-1.125e-3) !+ 0.01
                disp(i,3) = -0.05 * x(i,3)
              else ! top middle
                disp(i,1) = -0.05 * (x(i,1)-2.625e-3) !+ 0.02
                disp(i,2) = -0.05 * (x(i,2)-1.125e-3) !+ 0.01
                disp(i,3) = -0.05 * x(i,3)
              endif ! end top--switch head middle tail
            else if( x(i,2) .le. -0.6e-3 ) then ! bottom
              if ( x(i,1) .le. (dyn_org-1.025e-3) ) then ! bottom tail
                disp(i,1) = -0.05 * (x(i,1)-0.5e-3) !+ 0.02
                disp(i,2) = -0.05 * (x(i,2)+1.125e-3) !- 0.01
                disp(i,3) = -0.05 * x(i,3)
              else if ( x(i,1) .ge. (dyn_org+1.025e-3) ) then ! bottom head
                disp(i,1) = -0.05 * (x(i,1)-4.75e-3) !+ 0.02
                disp(i,2) = -0.05 * (x(i,2)+1.125e-3) !- 0.01
                disp(i,3) = -0.05 * x(i,3)
              else ! bottom middle
                disp(i,1) = -0.05 * (x(i,1)-2.625e-3) !+ 0.02
                disp(i,2) = -0.05 * (x(i,2)+1.125e-3) !- 0.01
                disp(i,3) = -0.05 * x(i,3)
              endif ! end bottom--switch head middle tail
            else ! middle
              if ( x(i,1) .le. (dyn_org-1.025e-3) ) then ! middle tail
                disp(i,1) = -0.05 * (x(i,1)-0.5e-3) !+ 0.02
                disp(i,2) = -0.05 * (x(i,2)) ! + 0.01
                disp(i,3) = -0.05 * x(i,3)
              else if ( x(i,1) .ge. (dyn_org+1.025e-3) ) then ! middle head
                disp(i,1) = -0.05 * (x(i,1)-4.75e-3) !+ 0.02
                disp(i,2) = -0.05 * (x(i,2)) ! + 0.01
                disp(i,3) = -0.05 * x(i,3)
              else ! middle middle
                disp(i,1) = -0.05 * (x(i,1)-2.625e-3) !+ 0.02
                disp(i,2) = -0.05 * (x(i,2)) ! + 0.01
                disp(i,3) = -0.05 * x(i,3)
              endif ! end middle--switch head middle tail
            endif ! end if switch top middle bottom
            BC(i,1:3)   = disp(i,1:3)
          endif ! end if inside cylinder
        enddo ! end loop numnp
        endif ! first time step
      endif ! end if case 5
c
c.... end test case 5
c
c
c.... test case 6
c.... original mov+rot+shrk
c
      if ( casenumber .eq. 6 ) then
        xtsl     = 2.0000000000000e-4
        dyn_org  = 2.6250000000000e-3 + DBLE(lstep) * xtsl
        if (myrank .eq. master) then
          write(*,*) "current lstep:", lstep, "dyn_org:", dyn_org
        endif
        rotf     = 90.0
        do i = 1,numnp
          if ( (ibits(iBC(i),14,3) .eq. 7) .and.
     &         (x(i,1) .lt. 6.65e-3) .and. (x(i,1) .gt. -1.4e-3) .and.
     &         ((x(i,2)*x(i,2) + x(i,3)*x(i,3)) .lt. 2.89e-6) ) then
            if( x(i,2) .ge. 0.6e-3 ) then ! top
              if ( x(i,1) .le. (dyn_org-1.025e-3) ) then ! top tail
                disp(i,1) = -0.1 * (x(i,1)-dyn_org+2.125e-3) + xtsl
     &                    + (x(i,1)-dyn_org+2.125e-3) * (cos(pi/rotf) - 1.0)
     &                    - (x(i,2)-1.125e-3) *  sin(pi/rotf)
                disp(i,2) = -0.1 * (x(i,2)-1.125e-3) !+ 0.01
     &                    + (x(i,2)-1.125e-3) * (cos(pi/rotf) - 1.0)
     &                    + (x(i,1)-dyn_org+2.125e-3) *  sin(pi/rotf)
                disp(i,3) = -0.1 * x(i,3)

              else if ( x(i,1) .ge. (dyn_org+1.025e-3) ) then ! top head
                disp(i,1) = -0.1 * (x(i,1)-dyn_org-2.125e-3) + xtsl
     &                    + (x(i,1)-dyn_org-2.125e-3) * (cos(pi/rotf) - 1.0)
     &                    - (x(i,2)-1.125e-3) *  sin(pi/rotf)
                disp(i,2) = -0.1 * (x(i,2)-1.125e-3) !+ 0.01
     &                    + (x(i,2)-1.125e-3) * (cos(pi/rotf) - 1.0)
     &                    + (x(i,1)-dyn_org-2.125e-3) *  sin(pi/rotf)
                disp(i,3) = -0.1 * x(i,3)

              else ! top middle
                disp(i,1) = -0.1 * (x(i,1)-dyn_org) + xtsl
     &                    + (x(i,1)-dyn_org) * (cos(pi/rotf) - 1.0)
     &                    - (x(i,2)-1.125e-3) *  sin(pi/rotf)
                disp(i,2) = -0.1 * (x(i,2)-1.125e-3) !+ 0.01
     &                    + (x(i,2)-1.125e-3) * (cos(pi/rotf) - 1.0)
     &                    + (x(i,1)-dyn_org) *  sin(pi/rotf)
                disp(i,3) = -0.1 * x(i,3)

              endif ! end top--switch head middle tail

            else if( x(i,2) .le. -0.6e-3 ) then ! bottom
              if ( x(i,1) .le. (dyn_org-1.025e-3) ) then ! bottom tail
                disp(i,1) = -0.1 * (x(i,1)-dyn_org+2.125e-3) + xtsl
     &                    - (x(i,1)-dyn_org+2.125e-3) * (cos(pi/rotf) - 1.0)
     &                    + (x(i,2)+1.125e-3) *  sin(pi/rotf)
                disp(i,2) = -0.1 * (x(i,2)+1.125e-3) !- 0.01
     &                    - (x(i,2)+1.125e-3) * (cos(pi/rotf) - 1.0)
     &                    - (x(i,1)-dyn_org+2.125e-3) *  sin(pi/rotf)
                disp(i,3) = -0.1 * x(i,3)

              else if ( x(i,1) .ge. (dyn_org+1.025e-3) ) then ! bottom head
                disp(i,1) = -0.1 * (x(i,1)-dyn_org-2.125e-3) + xtsl
     &                    - (x(i,1)-dyn_org-2.125e-3) * (cos(pi/rotf) - 1.0)
     &                    + (x(i,2)+1.125e-3) *  sin(pi/rotf)
                disp(i,2) = -0.1 * (x(i,2)+1.125e-3) !- 0.01
     &                    - (x(i,2)+1.125e-3) * (cos(pi/rotf) - 1.0)
     &                    - (x(i,1)-dyn_org-2.125e-3) *  sin(pi/rotf)
                disp(i,3) = -0.1 * x(i,3)

              else ! bottom middle
                disp(i,1) = -0.1 * (x(i,1)-dyn_org) + xtsl
     &                    - (x(i,1)-dyn_org) * (cos(pi/rotf) - 1.0)
     &                    + (x(i,2)+1.125e-3) *  sin(pi/rotf)
                disp(i,2) = -0.1 * (x(i,2)+1.125e-3) !- 0.01
     &                    - (x(i,2)+1.125e-3) * (cos(pi/rotf) - 1.0)
     &                    - (x(i,1)-dyn_org) *  sin(pi/rotf)
                disp(i,3) = -0.1 * x(i,3)

              endif ! end bottom--switch head middle tail

            else ! middle
              if ( x(i,1) .le. (dyn_org-1.025e-3) ) then ! middle tail
                disp(i,1) = -0.1 * (x(i,1)-dyn_org+2.125e-3) + xtsl
c     &                    + (x(i,1)-0.5e-3) * (cos(pi/rotf) - 1.0)
c     &                    - (x(i,2)) *  sin(pi/rotf)
                disp(i,2) = -0.1 * (x(i,2)) ! + 0.01
c     &                    + (x(i,2)) * (cos(pi/rotf) - 1.0)
c     &                    + (x(i,1)-0.5e-3) *  sin(pi/rotf)
                disp(i,3) = -0.1 * x(i,3)

              else if ( x(i,1) .ge. (dyn_org+1.025e-3) ) then ! middle head
                disp(i,1) = -0.1 * (x(i,1)-dyn_org-2.125e-3) + xtsl
c     &                    + (x(i,1)-4.75e-3) * (cos(pi/rotf) - 1.0)
c     &                    - (x(i,2)) *  sin(pi/rotf)
                disp(i,2) = -0.1 * (x(i,2)) ! + 0.01
c     &                    + (x(i,2)) * (cos(pi/rotf) - 1.0)
c     &                    + (x(i,1)-4.75e-3) *  sin(pi/rotf)
                disp(i,3) = -0.1 * x(i,3)

              else ! middle middle
                disp(i,1) = -0.1 * (x(i,1)-dyn_org) + xtsl
c     &                    + (x(i,1)-2.625e-3) * (cos(pi/rotf) - 1.0)
c     &                    - (x(i,2)) *  sin(pi/rotf)
                disp(i,2) = -0.1 * (x(i,2)) ! + 0.01
c     &                    + (x(i,2)) * (cos(pi/rotf) - 1.0)
c     &                    + (x(i,1)-2.625e-3) *  sin(pi/rotf)
                disp(i,3) = -0.1 * x(i,3)

              endif ! end middle--switch head middle tail

            endif ! end if switch top middle bottom
            BC(i,1:3)   = disp(i,1:3)
          endif ! end if inside cylinder
        enddo ! end loop numnp
      endif ! end if case 6
c
c.... end test case 6
c
c
c.... test case 7
c
      if ( casenumber .eq. 7 ) then
        do i = 1,numnp
          if ( (ibits(iBC(i),14,3) .eq. 7) .and.
     &        (x(i,1) .lt. 0.029) .and.
     &        (x(i,1) .gt.-0.009) ) then
c
            disp(i,1) = 0.001
            disp(i,2) = 0.0
            disp(i,3) = 0.0
c
            BC(i,:)   = disp(i,:)
          endif
        enddo ! end loop numnp
      endif ! end if case 7
c
c.... end test case 7
c
c
c.... test case 8
c.... mov+rot+shrk non-uniformly
c
      if ( casenumber .eq. 8 ) then
        xtsl     = 2.0000000000000e-4
        dyn_org  = 2.6250000000000e-3 + DBLE(lstep) * xtsl
        dyn_lnt  = 2.0000000000000e-3 * 0.9**lstep
        if (myrank .eq. master) then
          write(*,*) "current lstep:", lstep, "dyn_org:", dyn_org
        endif
        rotf     = 90.0
        shrkfactor = 2.0
        do i = 1,numnp
          if ( (ibits(iBC(i),14,3) .eq. 7) .and.
     &         (x(i,1) .lt. 6.65e-3) .and. (x(i,1) .gt. -1.4e-3) .and.
     &         ((x(i,2)*x(i,2) + x(i,3)*x(i,3)) .lt. 2.89e-6) ) then
            if( x(i,2) .ge. 0.6e-3 ) then ! top
              if ( x(i,1) .le. (dyn_org-1.025e-3) ) then ! top tail
                shrk = 0.5 - abs(x(i,1)-dyn_org+2.125e-3)/dyn_lnt
                shrk = shrkfactor * shrk
c
                disp(i,1) = -0.1 * (x(i,1)-dyn_org+2.125e-3) + xtsl
     &                    + (x(i,1)-dyn_org+2.125e-3) * (cos(pi/rotf) - 1.0)
     &                    - (x(i,2)-1.125e-3) *  sin(pi/rotf)
                disp(i,2) = -0.1 * (1.0 + shrk) * (x(i,2)-1.125e-3) !+ 0.01
     &                    + (x(i,2)-1.125e-3) * (cos(pi/rotf) - 1.0)
     &                    + (x(i,1)-dyn_org+2.125e-3) *  sin(pi/rotf)
                disp(i,3) = -0.1 * (1.0 + shrk) * x(i,3)

              else if ( x(i,1) .ge. (dyn_org+1.025e-3) ) then ! top head
                shrk = 0.5 - abs(x(i,1)-dyn_org-2.125e-3)/dyn_lnt
                shrk = shrkfactor * shrk
c
                disp(i,1) = -0.1 * (x(i,1)-dyn_org-2.125e-3) + xtsl
     &                    + (x(i,1)-dyn_org-2.125e-3) * (cos(pi/rotf) - 1.0)
     &                    - (x(i,2)-1.125e-3) *  sin(pi/rotf)
                disp(i,2) = -0.1 * (1.0 + shrk) * (x(i,2)-1.125e-3) !+ 0.01
     &                    + (x(i,2)-1.125e-3) * (cos(pi/rotf) - 1.0)
     &                    + (x(i,1)-dyn_org-2.125e-3) *  sin(pi/rotf)
                disp(i,3) = -0.1 * (1.0 + shrk) * x(i,3)

              else ! top middle
                shrk = 0.5 - abs(x(i,1)-dyn_org)/dyn_lnt
                shrk = shrkfactor * shrk
c
                disp(i,1) = -0.1 * (x(i,1)-dyn_org) + xtsl
     &                    + (x(i,1)-dyn_org) * (cos(pi/rotf) - 1.0)
     &                    - (x(i,2)-1.125e-3) *  sin(pi/rotf)
                disp(i,2) = -0.1 * (1.0 + shrk) * (x(i,2)-1.125e-3) !+ 0.01
     &                    + (x(i,2)-1.125e-3) * (cos(pi/rotf) - 1.0)
     &                    + (x(i,1)-dyn_org) *  sin(pi/rotf)
                disp(i,3) = -0.1 * (1.0 + shrk) * x(i,3)

              endif ! end top--switch head middle tail

            else if( x(i,2) .le. -0.6e-3 ) then ! bottom
              if ( x(i,1) .le. (dyn_org-1.025e-3) ) then ! bottom tail
                shrk = 0.5 - abs(x(i,1)-dyn_org+2.125e-3)/dyn_lnt
                shrk = shrkfactor * shrk
c
                disp(i,1) = -0.1 * (x(i,1)-dyn_org+2.125e-3) + xtsl
     &                    - (x(i,1)-dyn_org+2.125e-3) * (cos(pi/rotf) - 1.0)
     &                    + (x(i,2)+1.125e-3) *  sin(pi/rotf)
                disp(i,2) = -0.1 * (1.0 + shrk) * (x(i,2)+1.125e-3) !- 0.01
     &                    - (x(i,2)+1.125e-3) * (cos(pi/rotf) - 1.0)
     &                    - (x(i,1)-dyn_org+2.125e-3) *  sin(pi/rotf)
                disp(i,3) = -0.1 * (1.0 + shrk) * x(i,3)

              else if ( x(i,1) .ge. (dyn_org+1.025e-3) ) then ! bottom head
                shrk = 0.5 - abs(x(i,1)-dyn_org-2.125e-3)/dyn_lnt
                shrk = shrkfactor * shrk
c
                disp(i,1) = -0.1 * (x(i,1)-dyn_org-2.125e-3) + xtsl
     &                    - (x(i,1)-dyn_org-2.125e-3) * (cos(pi/rotf) - 1.0)
     &                    + (x(i,2)+1.125e-3) *  sin(pi/rotf)
                disp(i,2) = -0.1 * (1.0 + shrk) * (x(i,2)+1.125e-3) !- 0.01
     &                    - (x(i,2)+1.125e-3) * (cos(pi/rotf) - 1.0)
     &                    - (x(i,1)-dyn_org-2.125e-3) *  sin(pi/rotf)
                disp(i,3) = -0.1 * (1.0 + shrk) * x(i,3)

              else ! bottom middle
                shrk = 0.5 - abs(x(i,1)-dyn_org)/dyn_lnt
                shrk = shrkfactor * shrk
c
                disp(i,1) = -0.1 * (x(i,1)-dyn_org) + xtsl
     &                    - (x(i,1)-dyn_org) * (cos(pi/rotf) - 1.0)
     &                    + (x(i,2)+1.125e-3) *  sin(pi/rotf)
                disp(i,2) = -0.1 * (1.0 + shrk) * (x(i,2)+1.125e-3) !- 0.01
     &                    - (x(i,2)+1.125e-3) * (cos(pi/rotf) - 1.0)
     &                    - (x(i,1)-dyn_org) *  sin(pi/rotf)
                disp(i,3) = -0.1 * (1.0 + shrk) * x(i,3)

              endif ! end bottom--switch head middle tail

            else ! middle
              if ( x(i,1) .le. (dyn_org-1.025e-3) ) then ! middle tail
                shrk = 0.5 - abs(x(i,1)-dyn_org+2.125e-3)/dyn_lnt
                shrk = shrkfactor * shrk
c
                disp(i,1) = -0.1 * (x(i,1)-dyn_org+2.125e-3) + xtsl
                disp(i,2) = -0.1 * (1.0 + shrk) * (x(i,2)) ! + 0.01
                disp(i,3) = -0.1 * (1.0 + shrk) * x(i,3)

              else if ( x(i,1) .ge. (dyn_org+1.025e-3) ) then ! middle head
                shrk = 0.5 - abs(x(i,1)-dyn_org-2.125e-3)/dyn_lnt
                shrk = shrkfactor * shrk
c
                disp(i,1) = -0.1 * (x(i,1)-dyn_org-2.125e-3) + xtsl
                disp(i,2) = -0.1 * (1.0 + shrk) * (x(i,2)) ! + 0.01
                disp(i,3) = -0.1 * (1.0 + shrk) * x(i,3)

              else ! middle middle
                shrk = 0.5 - abs(x(i,1)-dyn_org)/dyn_lnt
                shrk = shrkfactor * shrk
c
                disp(i,1) = -0.1 * (x(i,1)-dyn_org) + xtsl
                disp(i,2) = -0.1 * (1.0 + shrk) * (x(i,2)) ! + 0.01
                disp(i,3) = -0.1 * (1.0 + shrk) * x(i,3)

              endif ! end middle--switch head middle tail

            endif ! end if switch top middle bottom
            BC(i,1:3)   = disp(i,1:3)
          endif ! end if inside cylinder
        enddo ! end loop numnp
      endif ! end if case 8
c
c.... end test case 8
c
c.... test case 9
c.... Update mesh2geom to identify vertex
c.... twist one end of a cylinder
c
      if ( casenumber .eq. 9 ) then
        do i = 1,numnp
c
          call core_is_in_closure(m2gClsfcn(i,1), m2gClsfcn(i,2),
     &                            3,              108,
     &                            answer)
c
          if (answer .ne. 0) then
c            disp(i,1) = 0.0015
c            disp(i,2) = 0.0
c            disp(i,3) = 0.0
            disp(i,1) = 0.0015
            disp(i,2) = x(i,2) * (cos(pi/100) - 1.0)
     &                - x(i,3) *  sin(pi/100)
            disp(i,3) = x(i,3) * (cos(pi/100) - 1.0)
     &                + x(i,2) *  sin(pi/100)
c
            BC(i,:)   = disp(i,:)
          endif ! end if inside box
        enddo ! end loop numnp
      endif  ! end case 9
c
c.... test case 10
c.... drive interface toward -x direction in stefan case
c
      if ( casenumber .eq. 10 ) then
        do i = 1,numnp
          if ( ifFlag(i) .eq. 1 ) then ! interface node
            disp(i,1) = -1.0 * Delt(1)
            disp(i,2) =  0.0 * Delt(1)
            disp(i,3) =  0.0 * Delt(1)
            BC(i,:) = disp(i,:)
          endif
        enddo ! end loop numnp
      endif  ! end case 10
c
c.... test case 11
c.... prescribe 6-grain displacement
c     mimic real simulation deformation
c     non-uniformly shrink; no rotation
c
      if ( casenumber .eq. 11 ) then
        xtsl    = 5.3333333333333e-3 * 0.98**lstep
        dyn_org = (5.3333333333333e-3 - 5.3333333333333e-3*0.98**lstep)/
     &            (1.0-0.98)
        dyn_lnt = 2.0000000000000e1 * 0.98**lstep
        shrkfactor = 2.0
c
        if (myrank .eq. master) then
          write(*,*) "current lstep:", lstep, "xtsl:", xtsl
        endif
c
        do i = 1,numnp
          if ( ifFlag(i) .eq. 1 ) then ! interface node
            if( x(i,2) .gt. 0.0 ) then ! top
              if ( x(i,1) .le. -1.25 ) then ! top tail
                shrk = 0.5 - abs(x(i,1)+2.5+dyn_org)/dyn_lnt
                shrk = shrkfactor * shrk
c
                disp(i,1) = -0.02 * (x(i,1)+2.5+dyn_org) - xtsl
                disp(i,2) = -0.02 * (1.0 + shrk) * (x(i,2)-0.75) !+ 0.01
                disp(i,3) = -0.02 * (1.0 + shrk) * x(i,3)
              else if ( x(i,1) .ge. 1.25 ) then ! top head
                shrk = 0.5 - abs(x(i,1)-2.5-dyn_org)/dyn_lnt
                shrk = shrkfactor * shrk
c
                disp(i,1) = -0.02 * (x(i,1)-2.5-dyn_org) + xtsl
                disp(i,2) = -0.02 * (1.0 + shrk) * (x(i,2)-0.75) !+ 0.01
                disp(i,3) = -0.02 * (1.0 + shrk) * x(i,3)
              else ! top middle
                shrk = 0.5 - abs(x(i,1))/dyn_lnt
                shrk = shrkfactor * shrk
c
                disp(i,1) = -0.02 * (x(i,1))
                disp(i,2) = -0.02 * (1.0 + shrk) * (x(i,2)-0.75) !+ 0.01
                disp(i,3) = -0.02 * (1.0 + shrk) * x(i,3)
              endif ! end top--switch head middle tail
            else ! bottom
              if ( x(i,1) .le. -1.25 ) then ! bottom tail
                shrk = 0.5 - abs(x(i,1)+2.5+dyn_org)/dyn_lnt
                shrk = shrkfactor * shrk
c
                disp(i,1) = -0.02 * (x(i,1)+2.5+dyn_org) - xtsl
                disp(i,2) = -0.02 * (1.0 + shrk) * (x(i,2)+0.75) !- 0.01
                disp(i,3) = -0.02 * (1.0 + shrk) * x(i,3)
              else if ( x(i,1) .ge. 1.25 ) then ! bottom head
                shrk = 0.5 - abs(x(i,1)-2.5-dyn_org)/dyn_lnt
                shrk = shrkfactor * shrk
c
                disp(i,1) = -0.02 * (x(i,1)-2.5-dyn_org) + xtsl
                disp(i,2) = -0.02 * (1.0 + shrk) * (x(i,2)+0.75) !- 0.01
                disp(i,3) = -0.02 * (1.0 + shrk) * x(i,3)
              else ! bottom middle
                shrk = 0.5 - abs(x(i,1))/dyn_lnt
                shrk = shrkfactor * shrk
c
                disp(i,1) = -0.02 * (x(i,1))
                disp(i,2) = -0.02 * (1.0 + shrk) * (x(i,2)+0.75) !- 0.01
                disp(i,3) = -0.02 * (1.0 + shrk) * x(i,3)
              endif ! end bottom--switch head middle tail
            endif ! end if switch top bottom
            BC(i,1:3)   = disp(i,1:3)
          endif
        enddo ! end loop numnp
      endif  ! end case 11
c
c.... test case 12
c.... shrink non-uniformly
c
      if ( casenumber .eq. 12 ) then
        xtsl     = 2.0000000000000e-4
        dyn_org  = 2.6250000000000e-3 ! + DBLE(lstep) * xtsl
        dyn_lnt  = 2.0000000000000e-3 * 0.9**lstep
        if (myrank .eq. master) then
          write(*,*) "current lstep:", lstep, "dyn_org:", dyn_org
        endif
        rotf     = 90.0
        shrkfactor = 2.0
        do i = 1,numnp
          if ( (ibits(iBC(i),14,3) .eq. 7) .and.
     &         (x(i,1) .lt. 6.65e-3) .and. (x(i,1) .gt. -1.4e-3) .and.
     &         ((x(i,2)*x(i,2) + x(i,3)*x(i,3)) .lt. 2.89e-6) ) then
            if( x(i,2) .ge. 0.6e-3 ) then ! top
              if ( x(i,1) .le. (dyn_org-1.025e-3) ) then ! top tail
                shrk = 0.5 - abs(x(i,1)-dyn_org+2.125e-3)/dyn_lnt
                shrk = shrkfactor * shrk
c
                disp(i,1) = -0.1 * (x(i,1)-dyn_org+2.125e-3) !+ xtsl
                disp(i,2) = -0.1 * (1.0 + shrk) * (x(i,2)-1.125e-3) !+ 0.01
                disp(i,3) = -0.1 * (1.0 + shrk) * x(i,3)

              else if ( x(i,1) .ge. (dyn_org+1.025e-3) ) then ! top head
                shrk = 0.5 - abs(x(i,1)-dyn_org-2.125e-3)/dyn_lnt
                shrk = shrkfactor * shrk
c
                disp(i,1) = -0.1 * (x(i,1)-dyn_org-2.125e-3) !+ xtsl
                disp(i,2) = -0.1 * (1.0 + shrk) * (x(i,2)-1.125e-3) !+ 0.01
                disp(i,3) = -0.1 * (1.0 + shrk) * x(i,3)

              else ! top middle
                shrk = 0.5 - abs(x(i,1)-dyn_org)/dyn_lnt
                shrk = shrkfactor * shrk
c
                disp(i,1) = -0.1 * (x(i,1)-dyn_org) !+ xtsl
                disp(i,2) = -0.1 * (1.0 + shrk) * (x(i,2)-1.125e-3) !+ 0.01
                disp(i,3) = -0.1 * (1.0 + shrk) * x(i,3)

              endif ! end top--switch head middle tail

            else if( x(i,2) .le. -0.6e-3 ) then ! bottom
              if ( x(i,1) .le. (dyn_org-1.025e-3) ) then ! bottom tail
                shrk = 0.5 - abs(x(i,1)-dyn_org+2.125e-3)/dyn_lnt
                shrk = shrkfactor * shrk
c
                disp(i,1) = -0.1 * (x(i,1)-dyn_org+2.125e-3) !+ xtsl
                disp(i,2) = -0.1 * (1.0 + shrk) * (x(i,2)+1.125e-3) !- 0.01
                disp(i,3) = -0.1 * (1.0 + shrk) * x(i,3)

              else if ( x(i,1) .ge. (dyn_org+1.025e-3) ) then ! bottom head
                shrk = 0.5 - abs(x(i,1)-dyn_org-2.125e-3)/dyn_lnt
                shrk = shrkfactor * shrk
c
                disp(i,1) = -0.1 * (x(i,1)-dyn_org-2.125e-3) !+ xtsl
                disp(i,2) = -0.1 * (1.0 + shrk) * (x(i,2)+1.125e-3) !- 0.01
                disp(i,3) = -0.1 * (1.0 + shrk) * x(i,3)

              else ! bottom middle
                shrk = 0.5 - abs(x(i,1)-dyn_org)/dyn_lnt
                shrk = shrkfactor * shrk
c
                disp(i,1) = -0.1 * (x(i,1)-dyn_org) !+ xtsl
                disp(i,2) = -0.1 * (1.0 + shrk) * (x(i,2)+1.125e-3) !- 0.01
                disp(i,3) = -0.1 * (1.0 + shrk) * x(i,3)

              endif ! end bottom--switch head middle tail

            else ! middle
              if ( x(i,1) .le. (dyn_org-1.025e-3) ) then ! middle tail
                shrk = 0.5 - abs(x(i,1)-dyn_org+2.125e-3)/dyn_lnt
                shrk = shrkfactor * shrk
c
                disp(i,1) = -0.1 * (x(i,1)-dyn_org+2.125e-3) !+ xtsl
                disp(i,2) = -0.1 * (1.0 + shrk) * (x(i,2)) ! + 0.01
                disp(i,3) = -0.1 * (1.0 + shrk) * x(i,3)

              else if ( x(i,1) .ge. (dyn_org+1.025e-3) ) then ! middle head
                shrk = 0.5 - abs(x(i,1)-dyn_org-2.125e-3)/dyn_lnt
                shrk = shrkfactor * shrk
c
                disp(i,1) = -0.1 * (x(i,1)-dyn_org-2.125e-3) !+ xtsl
                disp(i,2) = -0.1 * (1.0 + shrk) * (x(i,2)) ! + 0.01
                disp(i,3) = -0.1 * (1.0 + shrk) * x(i,3)

              else ! middle middle
                shrk = 0.5 - abs(x(i,1)-dyn_org)/dyn_lnt
                shrk = shrkfactor * shrk
c
                disp(i,1) = -0.1 * (x(i,1)-dyn_org) !+ xtsl
                disp(i,2) = -0.1 * (1.0 + shrk) * (x(i,2)) ! + 0.01
                disp(i,3) = -0.1 * (1.0 + shrk) * x(i,3)

              endif ! end middle--switch head middle tail

            endif ! end if switch top middle bottom
            BC(i,1:3)   = disp(i,1:3)
          endif ! end if inside cylinder
        enddo ! end loop numnp
      endif ! end if case 12
c
c.... end test case 12
c
c.... test case 13
c.... prescribe 6-grain displacement
c     non-uniformly shrink; no rotation
c
      if ( casenumber .eq. 13 ) then
        xtsl    = 5.3333333333333e-3 * 0.98**lstep
        dyn_org = (5.3333333333333e-3 - 5.3333333333333e-3*0.98**lstep)/
     &            (1.0-0.98)
        dyn_lnt = 2.0000000000000e1 * 0.85**lstep
        shrkfactor = 2.0
c
        if (myrank .eq. master) then
          write(*,*) "current lstep:", lstep, "xtsl:", xtsl
        endif
c
        do i = 1,numnp
          if ( ifFlag(i) .eq. 1 ) then ! interface node
            if( x(i,2) .gt. 0.0 ) then ! top
              if ( x(i,1) .le. -1.25 ) then ! top tail
                shrk = 0.5 - abs(x(i,1)+2.5+dyn_org)/dyn_lnt
                shrk = shrkfactor * shrk
c
                disp(i,1) = -0.02 * (x(i,1)+2.5+dyn_org) - xtsl
                disp(i,2) = -0.02 * (1.0 + shrk) * (x(i,2)-0.75) !+ 0.01
                disp(i,3) = -0.02 * (1.0 + shrk) * x(i,3)
              else if ( x(i,1) .ge. 1.25 ) then ! top head
                shrk = 0.5 - abs(x(i,1)-2.5-dyn_org)/dyn_lnt
                shrk = shrkfactor * shrk
c
                disp(i,1) = -0.02 * (x(i,1)-2.5-dyn_org) + xtsl
                disp(i,2) = -0.02 * (1.0 + shrk) * (x(i,2)-0.75) !+ 0.01
                disp(i,3) = -0.02 * (1.0 + shrk) * x(i,3)
              else ! top middle
                shrk = 0.5 - abs(x(i,1))/dyn_lnt
                shrk = shrkfactor * shrk
c
                disp(i,1) = -0.02 * (x(i,1))
                disp(i,2) = -0.02 * (1.0 + shrk) * (x(i,2)-0.75) !+ 0.01
                disp(i,3) = -0.02 * (1.0 + shrk) * x(i,3)
              endif ! end top--switch head middle tail
            else ! bottom
              if ( x(i,1) .le. -1.25 ) then ! bottom tail
                shrk = 0.5 - abs(x(i,1)+2.5+dyn_org)/dyn_lnt
                shrk = shrkfactor * shrk
c
                disp(i,1) = -0.02 * (x(i,1)+2.5+dyn_org) - xtsl
                disp(i,2) = -0.02 * (1.0 + shrk) * (x(i,2)+0.75) !- 0.01
                disp(i,3) = -0.02 * (1.0 + shrk) * x(i,3)
              else if ( x(i,1) .ge. 1.25 ) then ! bottom head
                shrk = 0.5 - abs(x(i,1)-2.5-dyn_org)/dyn_lnt
                shrk = shrkfactor * shrk
c
                disp(i,1) = -0.02 * (x(i,1)-2.5-dyn_org) + xtsl
                disp(i,2) = -0.02 * (1.0 + shrk) * (x(i,2)+0.75) !- 0.01
                disp(i,3) = -0.02 * (1.0 + shrk) * x(i,3)
              else ! bottom middle
                shrk = 0.5 - abs(x(i,1))/dyn_lnt
                shrk = shrkfactor * shrk
c
                disp(i,1) = -0.02 * (x(i,1))
                disp(i,2) = -0.02 * (1.0 + shrk) * (x(i,2)+0.75) !- 0.01
                disp(i,3) = -0.02 * (1.0 + shrk) * x(i,3)
              endif ! end bottom--switch head middle tail
            endif ! end if switch top bottom
            BC(i,1:3)   = disp(i,1:3)
          endif
        enddo ! end loop numnp
      endif  ! end case 13
c
      return
      end
c

max_stars_repo_name: SCOREC/phasta

max_issues_repo_name: SCOREC/phasta

max_forks_repo_name: SCOREC/phasta


text:
!
!     CalculiX - A 3-dimensional finite element program
!     Copyright (C) 1998-2015 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     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, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
!     construction of the B matrix
!
!     author: Yannick Muller
!
      subroutine flowoutput(itg,ieg,ntg,nteq,bc,lakon,ntmat_,
     &     v,shcon,nshcon,ipkon,kon,co,nflow, dtime,ttime,time,
     &     ielmat,prop,ielprop,nactdog,nacteq,iin,physcon,
     &     camt,camf,camp,uamt,uamf,uamp,rhcon,nrhcon,
     &     vold,jobnamef,set,istartset,iendset,ialset,nset,mi)
!
      implicit none
!
      character*8 lakon(*)
      character*132 jobnamef(*)
      character*81 set(*)
!
      integer mi(*),itg(*),ieg(*),ntg,nflow,ielmat(mi(3),*),i,
     &     nrhcon(*),node,
     &     ntmat_,nteq,nshcon(*),nelem,index,ipkon(*),kon(*),iin,
     &     nactdog(0:3,*),nacteq(0:3,*),ielprop(*),
     &     istartset(*),iendset(*),ialset(*),nset
!
      real*8 physcon(3),v(0:mi(2),*),shcon(0:3,ntmat_,*),co(3,*),
     &     prop(*),dtime,ttime,time,xflow,camt(*),camf(*),camp(*),
     &     rhcon(0:1,ntmat_,*),vold(0:mi(2),*),uamt,uamf,uamp,eta,
     &     bc(*)
!
      do i=1,nflow
         if(lakon(ieg(i))(2:5).eq.'LICH') then
            if((lakon(ieg(i))(6:7).eq.'SG').or.
     &           (lakon(ieg(i))(6:7).eq.'WE').or.
     &           (lakon(ieg(i))(6:7).eq.'DS')) then
               index=ipkon(ieg(i))
               node=kon(index+2)
               if(nactdog(3,node).eq.0) cycle
               index=ielprop(ieg(i))
               if(lakon(ieg(i))(6:7).eq.'SG') then
                  eta=prop(index+4)
                  nelem=int(prop(index+7))
               elseif(lakon(ieg(i))(6:7).eq.'WE') then
                  eta=prop(index+4)
                  nelem=int(prop(index+7))
               elseif(lakon(ieg(i))(6:7).eq.'DS') then
                  eta=prop(index+7)
                  nelem=int(prop(index+9))
               endif
               if(nelem.ne.0) then
                  write(*,*) '     *INFO in flowoutput: hydraulic jump'
                  write(*,*) '           in element ',nelem,'.'
                  write(*,*) '           relative location:',eta
                  write(*,*)
               endif
            endif
         endif
      enddo
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


	! Copyright 2005-2019 ECMWF. 中文注释: chinagod
	!
	! This software is licensed under the terms of the Apache Licence Version 2.0
	! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
	!
	! In applying this licence, ECMWF does not waive the privileges and immunities
	! granted to it by virtue of its status as an intergovernmental organisation
	! nor does it submit to any jurisdiction.
	!
	!
	!  Description: How to create and use and index to access messages from a file
	!               Creation of grib message, cloning one input message..
	!
	! 描述: 怎样从一个文件创建和使用index来访问消息
	!       创建grib消息，复制一个输入消息
	!
	program index
	use eccodes
	implicit none

	integer              :: iret
	integer,dimension(:),allocatable :: paramId,number
	integer                          :: paramIdSize,numberSize
	integer              :: i,j
	integer              :: idx,igrib,igribout,count,numberOfValues
	integer              :: outfile
	real (KIND=8),dimension(:), allocatable   :: values,result

	! create an index for the file 'eps' file using the keys number and paramId
	! use codes_index_create
	! 使用codes_index_createh函数，根据关键字number和paramId，创建文件'eps'的索引
	call codes_index_create(idx,'eps','paramId,number')

	! get the number of distinct values of paramId in the index
	! use codes_index_get_size
	! 使用codes_index_get_size函数，获得索引中paramID的值的数目
	call codes_index_get_size(idx,'paramId',paramIdSize)

	! allocate the array to contain the list of distinct paramId
	! 分配数组，用以存放paramId
	allocate(paramId(paramIdSize))

	! get the list of distinct paramId from the index
	! use codes_index_get
	!使用codes_index_get函数，根据索引读取paramId的值
	call codes_index_get(idx,'paramId',paramId)

	print*, 'grib contains ',paramIdSize, 'different parameters'

	! get the number of distinct values of number in the index
	! use codes_index_get_size
	!使用codes_index_get_size函数，从索引中获得number的值的数目
	call codes_index_get_size(idx,'number',numberSize)

	! allocate the array to contain the list of distinct numbers
	!分配数组，用以存放number
	allocate(number(numberSize))

	! get the list of distinct numbers from the index
	! use codes_index_get
	!使用codes_index_get函数，根据索引读取numbers的值
	call codes_index_get(idx,'number',number)

	print*, 'grib contains ',numberSize, 'different EPS members'

	count=0

	! select paramId 130 - Temperature
	! use codes_index_select
	!使用codes_index_select，选定paramID为130——对应变量为温度
	call codes_index_select(idx,'paramId','130')

	! loop over the different ensemble members
	!在不同集合成员直接循环
	do j=1,numberSize ! loop on number !number循环

		! select ensemble number=number(j)
		! use codes_index_select
		!使用codes_index_select，选定集合序号=number(j)
		call codes_index_select(idx,'number',number(j))

		! get one grib message for the above values of the index
		! use codes_new_from_index
		!使用codes_new_from_index，根据上述索引的值，得到一个grib的消息
		call codes_new_from_index(idx,igrib)

		! for the first filed allocate array for values and result
		! use codes_get_size
		!使用use codes_get_size函数获取数组大小，对第一个场，为values和results分配数组
		if ( j == 1 ) then
			call codes_get_size(igrib,'values',numberOfValues)
			allocate(values(numberOfValues), stat = iret)
			if (iret /= 0)  STOP 'Failed to allocate values'
			allocate(result(numberOfValues), stat = iret)
			if (iret /= 0)  STOP 'Failed to allocate result'

			! clone grib message
			! use codes_clone:
			!   codes_clone (gribid_src, gribid_dest [, status])
			!使用codes_clone (gribid_src, gribid_dest [, status])函数，复制grib消息
			call codes_clone( igrib, igribout)

		end if

		! get data values
		! use codes_get
		!使用use codes_get，获得data的值
		call codes_get(igrib,"values", values)

		count = count + 1

		result=result+values

		! release the grib message
		! use grib_release
		!使用grib_release，关闭grib消息
		call codes_release(igrib)

	end do ! loop on number

	! release the index
	! use codes_index_release
	!使用grib_release，关闭索引
	call codes_index_release(idx)

	print*,'We considered ',count,' members'

	result=result/count

	print*,'=============================================================================='
	print*, 'Stats for ensemble mean of T850'
	print*, 'Min: ', minval(result), 'Max: ', maxval(result), 'Avg: ', sum(result)/numberOfValues
	print*,'=============================================================================='

	! open output grib file eps_mean.grib
	! use codes_open_file:
	!   codes_open_file (ifile, filename, mode [, status])
	!使用codes_open_file，打开输出文件eps_mean.grib
	call codes_open_file(outfile, 'eps_mean.grib','w')

	! change grib headers
	! See  http://apps.ecmwf.int/codes/grib/format/grib1/local/ (definition 1 or
	! 30) and
	!      http://apps.ecmwf.int/codes/grib/format/mars/type/ for the datatype
	! Keys to change are perturbationNumber and dataType.

	! use codes_set:
	!   codes_set (gribid, key, value [, status])
	!使用codes_set，更改grib头

	call codes_set(igribout,'perturbationNumber',0)
	call codes_set(igribout,'dataType',"em")

	! change data values
	!更改数据值
	call codes_set(igribout,'values',result)

	! write message out
	! use codes_write:
	!   codes_write (gribid, ifile [, status])
	!使用codes_write，写入更改
	call codes_write(igribout,outfile)

	end program index


text:
program problem22
  implicit none
  integer, parameter         :: n=5163,m=20
  character(m), dimension(n) :: names
  character(m)               :: temp
  integer                    :: i,j,answer

  open(33,file='./data/problem22.dat')
  read(33,*) names

  ! sort
  do i=1,n
     do j=i,n
        if (names(i) > names(j)) then
           temp=names(i)
           names(i)=names(j)
           names(j)=temp
        end if
     end do
  end do

  answer=0
  do i=1,n
     answer=answer+(i*score(names(i)))
  end do

  print *, answer

contains
  ! score string
  pure function score(name)
    implicit none
    character(m), intent(in) :: name
    integer                  :: score,i

    score=0
    do i=1,len_trim(name)
       score=score+ichar(name(i:i))-ichar('@')
    end do

  end function score

end program problem22

max_stars_repo_name: jamesmcclain/Dioskouroi

max_issues_repo_name: jamesmcclain/Dioskouroi

max_forks_repo_name: jamesmcclain/Dioskouroi


text:
      subroutine bmpinit

!!    ~ ~ ~ PURPOSE ~ ~ ~
!!    This subroutine sets default values for urban bmp parameters

!!    ~ ~ ~ INCOMING VARIABLES ~ ~ ~ ~ ~ ~ ~
!!    name           |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    dtp_onoff(:)   |none          |sub-basin detention pond is associated with
!!    dtp_iy(:)      |none          |year of the simulation that the reservoir
!!                                  |becomes operational
!!    dtp_imo(:)     |none          |month the reservoir becomes operational
!!    dtp_evrsv      |none          |detention pond evaporation coefficient
!!    dtp_numweir(:) |none          |Total number of weirs in the BMP
!!    dtp_numstage(:)|none          |Total number of stages in the weir
!!    dtp_lwratio(:)   |none          |Ratio of length to width of water back up
!!    dtp_totwrwid(:)|m             |Total constructed width of the detention wall across
!!                                  |the creek
!!    dtp_stagdis(:) |none          |0=use weir/orifice discharge equation to calculate
!!                                  |outflow, 1=use stage-dicharge relationship
!!    dtp_reltype(:) |none          |Equations for Stage-Discharge relationship,1=exponential
!!                                  |function, 2=linear, 3=logarithmic, 4=cubic, 5=power
!!    dtp_intcept(:) |none          |Intercept used in regression equations
!!    dtp_expont(:)  |none          |Exponent used in the exponential equation
!!    dtp_coef1(:)   |none          |Coefficient of 3rd degree in the polynomial equation
!!    dtp_coef2(:)   |none          |Coefficient of 2nd degree in the polynomial equation
!!    dtp_coef3(:)   |none          |Coefficient of 1st degree in the polynomial equation
!!    dtp_dummy1(:)   |none         |Dummy variable, backs up space (Rename for use)
!!    dtp_dummy2(:)   |none         |Dummy variable, backs up space (Rename for use)
!!    dtp_dummy3(:)   |none         |Dummy variable, backs up space (Rename for use)
!!    dtp_weirtype(:,:)|none        |Type of weir: 1=rectangular and 2=circular
!!    dtp_weirdim(:,:)|none         |Weir dimensions, 1=read user input, 0=use model calculation
!!    dtp_wdratio(:,:)|none         |Width depth ratio of rectangular weirs
!!    dtp_depweir(:,:)|m            |Depth of rectangular wier at different stages
!!    dtp_diaweir(:,:)|m            |Diameter of orifice hole at different stages
!!    dtp_addon(:,:)  |m            |The distance between spillway levels
!!    dtp_flowrate(:,:)|m3/sec      |Maximum discharge from each stage of the weir/hole
!!    dtp_cdis(:,:)  |none          |Discharge coeffieicne for weir/orifice flow
!!    dtp_retperd(:,:)|years        |Return period at different stages
!!    dtp_pcpret(:,:)|mm            |precipitation for different return periods (not used)
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ OUTGOING VARIABLES ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    name           |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    res_out(:,:,:) |m**3/day      |measured average daily outflow from the
!!                                  |reservoir for the month (needed if IRESCO=1)
!!                                  |(read in as m**3/s and converted to m**3/day)
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ LOCAL DEFINITIONS ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    eof         |none          |end of file flag
!!    k           |none          |counter
!!    lnvol       |none          |variable to hold denominator value
!!    mon         |none          |counter
!!    titldum     |NA            |title line in .det file (not used in program)
!!    resmto      |NA            |name of reservoir outflow file
!!                               |(needed if IDETCO = 2)
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

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

      use parm
      implicit none
      integer :: j, k, eof,kk
      real*8 :: hwq,wqv,sub_ha,bmpfr_sf,bmpfr_ri, qstg, hstg

      eof = 0; bmpfr_sf=0.; bmpfr_ri=0.; hstg = 0.; qstg=0.
      sub_ha = sub_km(i) * 100.

      !! Detention pond
      !!----------------
      if (dtp_onoff(i)==1) then
         if (dtp_imo(i)<=0)      dtp_imo(i) = 1
         if (dtp_iyr(i)<=1000)   dtp_iyr(i) = iyr
   	   if (dtp_evrsv(i)<=0)    dtp_evrsv(i) = 0.1
   	   if (dtp_lwratio(i)>1)   dtp_lwratio(i) = 1
	   if (dtp_numweir(i)<=0)  dtp_numweir(i) = 1
	   if (dtp_numstage(i)<=0) dtp_numstage(i) = 1
	   if (dtp_numstage(i)>1) then
	      do k=2,dtp_numstage(i)
	          if (dtp_weirtype(i,k)==1) dtp_addon(i,k) = 0.
	      end do
	   endif

         !!	Estimating design flow rate if not entered by user
         do k=1,dtp_numstage(i)
            if (dtp_flowrate(i,k)<=0.0) then
                dtp_flowrate(i,k) = 0.5 * 1000.0 * dtp_pcpret(i,k)
     &            * subdr_km(i) / (idt*60.)
	      end if

	      if (dtp_flowrate(i,k)<0) then
	         print *,"Error.. Could not estimate emergency spillway volume"
               print *,"Please enter necessary data in *.pnd input file"
               print *,"for subbasin : ",i
   !!            stop
	      end if
         end do

         !!	Separate cumulative flow information to individual weir
         do k=2,dtp_numstage(i)
            dtp_flowrate(i,k) = dtp_flowrate(i,k) / dtp_numweir(i)
         end do

         !!Estimate weir dimensions based on existing data
         do k=1,dtp_numstage(i)
            if (dtp_weirdim(i,k)==0) then  !! Estimating weir dimensions
               if (dtp_weirtype(i,k)==2) then	!! choosing weir type
                dtp_diaweir(i,k)=(0.479081 * dtp_flowrate(i,k)
     &                              / dtp_cdis(i,k))**0.4
               else !rectangular weir
			      if (k==1) then
				     hstg = dtp_depweir(i,k)
				     dtp_wrwid(i,k) = dtp_flowrate(i,k)
     &                   / (1.84*dtp_cdis(i,k)*hstg**1.5)
                  else
				     qstg = 0.
				     do j=1,k-1
					    hstg = sum(dtp_depweir(i,j:k))
					    qstg = dtp_cdis(i,k) * 1.84
     &                   * dtp_wrwid(i,j) * hstg ** 1.5 !m3/s
					    dtp_flowrate(i,k) = max(0.,dtp_flowrate(i,k)-qstg)

					 end do
                     dtp_wrwid(i,k) = dtp_flowrate(i,k)
     &                  / (1.84*dtp_cdis(i,k)*dtp_depweir(i,k)**1.5)
                  endif
               end if
            else  !! read user-entered data
               if (dtp_weirtype(i,k)==1) then
                  dtp_wrwid(i,k) = dtp_wdratio(i,k) * dtp_depweir(i,k)
               end if
            end if
         end do
          !! divide rectangular weirs into multiple single stage weirs
         do k = 2, dtp_numstage(i)
            dtp_addon(i,k) = dtp_addon(i,k-1) + dtp_depweir(i,k-1)
         end do

         ! weir depth from the top to the bottom of each stage
		 do k = dtp_numstage(i), 2, -1
            dtp_depweir(i,k-1) = dtp_depweir(i,k) + dtp_depweir(i,k-1)
         end do
      end if

      !! Wet pond
      !!----------
      if (wtp_onoff(i)==1) then
         if (wtp_imo(i)<=0)      wtp_imo(i) = 1
         if (wtp_iyr(i)<=1000)   wtp_iyr(i) = iyr
         if (wtp_k(i) <=0)       wtp_k(i) = 0.01
         if (wtp_evrsv(i)<=0)    wtp_evrsv(i) = 0.01
         if (wtp_evrsv(i)>1)     wtp_evrsv(i) = 0.99
         if (wtp_sdslope(i)<3)   wtp_sdslope(i) = 3. !COA Manual 1.6.6.C
         if (wtp_lenwdth(i)<2)   wtp_lenwdth(i) = 2. !COA Manual 1.6.6.C
         if (wtp_pdia(i)<=0)     wtp_pdia(i) = 0.1524 !meter(=6inches), COA Manual 1.6.6
         if (wtp_plen(i)<=2)     wtp_plen(i) = 2 !meter
         if (wtp_pmann(i)<=0)    wtp_pmann(i) = 0.012 ! concrete surface
         if (wtp_ploss(i)<=0)    wtp_ploss(i) = 0
      end if

      !fraction runoff that directly enters the channel
      if(num_sf(i)>=1) then
         do kk=1,num_sf(i)
             bmpfr_sf = bmpfr_sf + sf_fr(i,kk)
         end do
      endif
      if(num_ri(i)>=1) then
         do kk=1,num_ri(i)
             bmpfr_ri = bmpfr_ri + ri_fr(i,kk)
         end do
      endif
      if (bmpfr_sf>1.or.bmpfr_ri>1) then
         write (*,*) " "
         write (*,*) "Urban BMP Warning!!"
         write (*,*) "In subbasin ", i
         write (*,*) "The fraction runoff draining to urban BMPs"
         write (*,*) " are larger than one, so the fraction values"
         write (*,*) " were automatically reassigned"
         write (*,*) " "
         do kk=1,num_sf(i)
             sf_fr(i,kk) = sf_fr(i,kk) / bmpfr_sf
         end do
          do kk=1,num_sf(i)
             ri_fr(i,kk) = ri_fr(i,kk) / bmpfr_ri
         end do
         bmpfr_sf = 1.; bmpfr_ri=1.
      endif

      !!Retention-Irrigation
      !!---------------------
      do k=1,num_ri(i)
         ! skip the pond that has zero inflow
         if (ri_fr(i,k)==0) cycle

         ! determine water quality volume for defult pond sizes
         !City of Austin Design Guideline 1.6.9.2 Table 1-12
         !Retention-Irrigation for Barton Springs Zone

         hwq = (1.8 * sub_ha_imp(i) / sub_ha_urb(i) + 0.6) !inches
         wqv = hwq / 12. * sub_ha_urb(i) * ri_fr(i,k) * 107639.104167 !ft3

         if (ri_dim(i,k)==0) then
           !Determine pond size automatically based on City of Austin's Design Guideline 1.6
            ri_vol(i,k) = wqv * 0.028317 !m3
		      ri_dep(i,k)=1.5 !assume 1.5m as default retention pond depth
            ri_sa(i,k) = ri_vol(i,k) / ri_dep(i,k)
            ri_dd (i,k)=60.0 !drawdown time, hr
            ri_k(i,k)=2.5
            ri_evrsv(i,k)=0.6
         else
            !Evaluate unit variables provided by user
            if (ri_sa(i,k)<1.or.ri_vol(i,k)<1) then
               !City of Austin Design Guideline 1.6.5
               ri_vol(i,k) = wqv * 0.028317 !m3
               ri_sa(i,k) = ri_vol(i,k) / ri_dep(i,k)
            end if

		   ri_dep(i,k) = ri_vol(i,k) / ri_sa(i,k)
		   if (ri_dd(i,k)<=0) ri_dd (i,k)=72.0
		   if (ri_k(i,k) <=0 ) ri_k(i,k)=2.5
		   if (ri_evrsv(i,k) <=0 ) ri_evrsv(i,k)=0.1
		   if (ri_dep(i,k)<=1) ri_dep(i,k)=1.
         end if

	   ! draw down time [number of time step]
	   ri_ndt(i,k) = (ri_dd(i,k) - 12.) * 60 / idt  !minus the first dry 12 hours no-pumping

	   ! pumping rate that empties the basin in 72 hours with initial 12 hour of no operatrion
	   ri_pumpv(i,k) = ri_vol(i,k) / ri_ndt(i,k) !m3/dt

         if (ri_im(i,k)<0.or.ri_im(i,k)>12) ri_im(i,k) = 0
         if (ri_iy(i,k)<1000.and.ri_iy(i,k)>0) ri_iy(i,k) = 0
         if (ri_iy(i,k)==0) ri_iy(i,k) = iyr
         if (ri_im(i,k)==0) ri_im(i,k) = 1

       write(77779,'(a11,i5)') 'Subbasin #:', i   ! bmp_sedfil.out
       write(77779,'(a46)') ''
       write(77779,'(a10,i5)') 'RI #:', K   ! bmp_sedfil.out
       write(77779,'(a17,f10.1,a4)') 'Total volume =', ri_vol(i,k),'m^3'
       write(77779,'(a17,f10.1,a4)') 'Surface area =', ri_sa(i,k),'m^2'
       write(77779,'(a17,f10.1,a3)') 'Drawdown time =', ri_dd (i,k),'hr'
       write(77779,'(a17)') ''

      end do

      !!Sedimentation-Filtration
      !!---------------------

      do k=1,num_sf(i)
         write(77778,'(a11,i5)') 'Subbasin #:', i   ! bmp_sedfil.out
         write(77778,'(a46)') ''
         write(77778,'(a10,i5)') 'SED-FIL #:', K   ! bmp_sedfil.out
         !determine water quality volume for defult pond sizes
         !City of Austin Design Guideline 1.6.2
         hwq = (0.5 + sub_ha_imp(i) / sub_ha_urb(i) - 0.2) !inches
         wqv = hwq / 12. * sub_ha_urb(i) * sf_fr(i,k) * 107639.104167 !ft3

         if (sf_dim(i,k)==0) then
            write(77778,'(a46)') 'This SED-FIL size is automatically'
     &         // ' estimated based on WQV.'
            !Determine pond size automatically based on City of Austin's Design Guideline 1.6
            if (sf_typ(i,k)==1.or.sf_typ(i,k)==3) then
               ! full scale or sedimentation basin only
               sp_pvol(i,k) = wqv * 0.028317 !m3
               sp_sa(i,k) = sp_pvol(i,k) / 1.5 !assume 1.5m depth for sed pond
               sp_pd(i,k) = sqrt(4. * 2. * sp_sa(i,k) * 1.5**0.5
     &                 / (0.6 * 172800. * 19.6**0.5) / 3.14159) * 1000. !mm
               ft_sa(i,k) = wqv/(7.+2.33*4.) * 0.093 !m2
               ft_fsa(i,k) = 1.

            else
               ! partial scale
               ft_sa(i,k) = wqv * 0.028317 !m3
!               wqv/(4.+1.33*4.) * 0.093

            end if
            sp_bpw(i,k) = 10. !m overflow weir width
            ft_pd(i,k) = 1524. !mm
            ft_dep(i,k) = 420. !mm
            ft_h(i,k) = 1200. !mm

         else
         !Evaluate unit variables given by users
            if (sf_typ(i,k)>3) sf_typ(i,k) = 1
            if (sp_sa(i,k)<1) then
               !City of Austin Design Guideline 1.6.5
               sp_pvol(i,k) = wqv * 0.028317 !m3
               sp_sa(i,k) = sp_pvol(i,k) / 1.5 !assume 1.5m depth for sed pond
            end if
            if (sp_pd(i,k)<0.1) sp_pd(i,k) = 152.4 !mm diameter, equivalent to 6in
            if (ft_h(i,k)<5) ft_h(i,k) = 1200. !mm
            if (ft_sa(i,k)<1) then
               if (sf_typ(i,k)==1) then
                  ft_sa(i,k) = wqv/(7.+2.33*ft_h(i,k)/304.8) * 0.093 !m2
                  ft_fsa(i,k) = 1.
               else
                  ft_sa(i,k) = wqv * 0.028317 !m3
                  ft_fsa(i,k) = 1. /(4.+1.33*ft_h(i,k)/304.8)
               end if
            end if
         end if

         !Outflow control
         if (sp_qfg(i,k)==0) then
            sp_pd(i,k) = sqrt(4. * 2. * sp_sa(i,k) * 1.5**0.5
     &           / (0.6 * 172800. * 19.6**0.5) / 3.14159) * 1000. !mm
         end if
         if (ft_qfg(i,k)==0) then
            ft_pd(i,k) = sqrt(4. * 2. * ft_sa(i,k) * 1.5**0.5
     &          / (0.6 * 172800. * 19.6**0.5) / 3.14159) * 1000. !mm
         end if

         !Orifice pipe for sand filter should be equal or larger than
         !sedimentation pond outlet pipe for full-type SedFils
         if (ft_pd(i,k)<sp_pd(i,k)) ft_pd(i,k) = sp_pd(i,k)

         if (ft_dep(i,k)<100) ft_dep(i,k) = 100.
         if (sf_ptp(i,k)>1) sf_ptp(i,k) = 1
   !      if (sf_typ(i,k)==1) sf_ptp(i,k) = 0 removed by Jaehak 2014
   !      if (sp_pd(i,k)>254) sp_pd(i,k) = 254. ! max 10inches dia
         if (sp_pd(i,k)<10) sp_pd(i,k) = 10. ! min 10mm dia
   !     if (ft_pd(i,k)>254) ft_pd(i,k) = 254. ! max 10inches
         if (ft_k(i,k)<1) ft_k(i,k) = 1.
         if (ft_dp(i,k)<0.0001) ft_dp(i,k) = 0.02
         if (ft_dc(i,k)<0.01) ft_dc(i,k) = 0.762
         if (ft_por(i,k)<0.1) ft_por(i,k) = 0.45
         if (tss_den(i,k)<0.5) tss_den(i,k) = 0.5
         if (ft_alp(i,k)<0.1) ft_alp(i,k) = 0.1
         if (sf_im(i,k)<0.or.sf_im(i,k)>12) sf_im(i,k) = 0
         if (sf_iy(i,k)<1000.and.sf_iy(i,k)>0) sf_iy(i,k) = 0
         if (sf_iy(i,k)==0) sf_iy(i,k) = iyr
         if (sf_im(i,k)==0) sf_im(i,k) = 1
         if (ft_fsa(i,k)==0) ft_fsa(i,k) = 0.85
         if (sf_typ(i,k)==1) ft_fsa(i,k) = 1

         if (sf_typ(i,k)==1) then
           write(77778,'(a37)') 'Full type sed-fil selected'
         elseif (sf_typ(i,k)==2) then
           write(77778,'(a40)') 'Partial type sed-fil selected'
         else
           write(77778,'(a43)') 'Sedimentation pond only selected'
         endif
         if (sf_typ(i,k)==1.or.sf_typ(i,k)==3) then
      write(77778,'(a18)') 'Sedimentation pond'
      write(77778,'(a17,f10.1,a4)') 'Total volume =', sp_pvol(i,k),'m^3'
      write(77778,'(a17,f10.1,a4)') 'Surface area =', sp_sa(i,k),'m^2'
      write(77778,'(a17,f10.1,a3)') 'Drain Pipe Dia =', sp_pd(i,k),'mm'
      write(77778,'(a17)') ''
         endif
         if (sf_typ(i,k)==1.or.sf_typ(i,k)==2) then
      write(77778,'(a11)') 'Sand Filter'
      write(77778,'(a17,f10.1,a4)') 'Surface area =', ft_sa(i,k),'m^2'
      write(77778,'(a17,f10.1,a3)') 'Max ponding =', ft_h(i,k),'mm'
      write(77778,'(a17,f10.1,a3)') 'Filter depth =', ft_dep(i,k),'mm'
      write(77778,'(a17,f10.1,a3)') 'Drain Pipe Dia =', ft_pd(i,k),'mm'
      write(77778,'(a17)') ''
        end if
        end do

      return
      end

max_stars_repo_name: GISWAT/erosion-sediment

max_issues_repo_name: GISWAT/erosion-sediment

max_forks_repo_name: GISWAT/erosion-sediment


text:
!#####################################################################
module error_handler
   ! Generic error message handler for both console and non-console
   ! applications. The routines in this module do not halt program
   ! execution; they merely store error messages for subsequent
   ! retrieval.
   ! Written by J.L. Schafer, 4/19/02; revised 6/03/03
   implicit none
   private ! by default
   ! declare public types
   public :: error_type
   ! declare public subroutines and functions
   public :: err_reset, err_handle, err_msg_present, err_get_msgs
   ! Parameters private to this module
   integer, parameter :: &
        ! max width of any single error message line
        err_msg_width = 70
   !##################################################################
   type :: msg_line_type
      ! Private type for a single node in the linked list
      sequence
      character (len=err_msg_width) :: line = ""
      type(msg_line_type), pointer :: next => null()
   end type msg_line_type
   !##################################################################
   type :: error_type
      ! Public type for holding a linked list of messages
      sequence
      private ! contents of this type are private to this module
      logical :: msg_present = .false.
      type(msg_line_type), pointer :: head => null(), tail => null()
   end type error_type
   !##################################################################
contains
   !##################################################################
   subroutine insert_msg_line( text_line, err )
      ! inserts a new message line at the end of the list
      implicit none
      ! declare arguments
      character(len=*), intent(in) :: text_line
      type(error_type), intent(inout) :: err
      ! begin
      if( .not. err%msg_present ) then
         ! begin a linked list
         allocate( err%head )
         err%tail => err%head
         err%head%line = text_line
         err%msg_present = .true.
      else
         ! add a node to the list; point tail to the new node
         allocate( err%tail%next )
         err%tail => err%tail%next
         err%tail%line = text_line
      end if
   end subroutine insert_msg_line
   !##################################################################
   subroutine err_reset( err )
      ! Public: deletes all messages from the list
      implicit none
      type(error_type), intent(inout) :: err
      type(msg_line_type), pointer :: current_line
      if( .not. err%msg_present ) return
      do
         current_line => err%head
         err%head => err%head%next
         deallocate( current_line )
         if( .not.associated(err%head) ) exit
      end do
      nullify(err%tail)
      err%msg_present = .false.
   end subroutine err_reset
   !##################################################################
   subroutine err_handle(err, err_code, called_from, file_name, &
        line_no, object_name, custom_1, custom_2, custom_3)
      ! Public: Stores a message in the error handler
      ! Meaning of err_code
      !     0 = no error
      ! 1-99: I/O errors
      !     1 = file could not be opened for read-access
      !     2 = file could not be opened for write-access
      !     3 = error in reading file
      !     4 = error in writing to file
      !     5 = error in reading file: EOR/EOF encountered
      !     6 = file open for write-access could not be closed
      !     7 = file open for read-access could not be closed
      ! 100-199: numerical errors
      !   100 = matrix apparently singular
      !   101 = matrix not positive definite
      !   102 = attempted division by zero
      !   103 = attempted logarithm of non-positive number
      !   104 = argument to exp function too large
      !   105 = attempted square root of negative number
      ! 200-299: memory errors
      !   200 = unable to dynamically allocate memory for object
      !   201 = unable to deallocate memory for object
      ! 300-399: array dimension errors
      !   300 = non-square matrix encountered where square
      !     matrix expected
      !   301 = dimensions of matrix arguments not conformable
      ! 1000: other error
      !  1000 = reserved for custom error messages
      implicit none
      ! declare required arguments
      type(error_type), intent(inout) :: err
      integer, intent(in) :: err_code
      ! declare optional arguments
      character (len=*), optional :: called_from
      character (len=*), optional :: file_name
      character (len=*), optional :: object_name
      character (len=*), optional :: custom_1
      character (len=*), optional :: custom_2
      character (len=*), optional :: custom_3
      integer, optional :: line_no
      ! local variables
      character(len=12) :: ichar
      ! begin
      select case(err_code)
         case(0)
            call insert_msg_line( &
                 "No errors", err)
         ! I/O errors
         case(1)
            call insert_msg_line( &
                 "File could not be opened for read-access", err)
         case(2)
            call insert_msg_line( &
                 "File could not be opened for write-access", err)
         case(3)
            call insert_msg_line( &
                 "Error in reading file", err)
         case(4)
            call insert_msg_line( &
                 "Error in writing to file", err)
         case(5)
            call insert_msg_line( &
                 "Error in reading file: EOR or EOF encountered", err)
         case(6)
            call insert_msg_line( &
                 "File open for write-access could not be closed", &
                 err)
         case(7)
            call insert_msg_line( &
                 "File open for read-access could not be closed", err)
         ! numerical errors
         case(100)
            call insert_msg_line( &
                 "Matrix apparently singular", err)
         case(101)
            call insert_msg_line(  &
                 "Matrix not positive definite", err)
         case(102)
            call insert_msg_line(  &
                 "Attempted division by zero", err)
         case(103)
            call insert_msg_line( &
                 "Attempted logarithm of non-positive number", err)
         case(104)
            call insert_msg_line( &
                 "Argument to exp function too large", err)
         case(105)
            call insert_msg_line( &
                 "Attempted square root of negative number", err)
         ! memory errors
         case(200)
            call insert_msg_line(  &
                 "Unable to allocate memory for object", err)
         case(201)
            call insert_msg_line( &
                 "Unable to deallocate memory for object", err)
         ! array dimension errors
         case(300)
            call insert_msg_line( &
                 "Non-square matrix encountered; square expected", &
                 err)
         case(301)
            call insert_msg_line( &
                 "Dimensions of matrix arguments not conformable", &
                 err)
         ! custom error message
         case(1000)
            ! don't do anything yet
         ! anything else
         case default
            call insert_msg_line("Unknown error code.", err)
      end select
      ! append other optional information if present
      if( present(custom_1) ) then
         call insert_msg_line(custom_1, err)
      end if
      if( present(custom_2) ) then
         call insert_msg_line(custom_2, err)
      end if
      if( present(custom_3) ) then
         call insert_msg_line(custom_3, err)
      end if
      if(present(file_name)) then
         call insert_msg_line("FILE: " // trim(file_name), err)
      end if
      if(present(line_no)) then
         write(ichar,"(I12)") line_no
         ichar = adjustl(ichar)
         call insert_msg_line("LINE: " // trim(ichar), err)
      end if
      if(present(object_name)) then
         call insert_msg_line(trim(object_name), err)
      end if
      if(present(called_from)) then
         call insert_msg_line("OCCURRED IN: " // &
              trim(called_from), err)
      end if
   end subroutine err_handle
   !##################################################################
   logical function err_msg_present(err)
      ! Public: Queries the error_type to see if a message is present
      implicit none
      type(error_type), intent(inout) :: err
      err_msg_present = err%msg_present
   end function err_msg_present
   !##################################################################
   subroutine err_get_msgs(err, msg_string, platform)
      ! Public: Retrieves all stored messages as a single character
      ! string, with message lines separated by platform-appropriate
      ! ASCII carriage control characters.
      ! Values for platform may be "UNIX", "MAC" or "PC"
      implicit none
      ! required arguments
      type(error_type), intent(inout) :: err
      character(len=*), intent(out) :: msg_string
      ! optional arguments
      character(len=*), intent(in), optional :: platform
      ! local variables
      character(len=4) :: plat
      integer :: posn
      logical :: first_time
      type(msg_line_type), pointer :: cur_line
      ! determine platform
      if( present(platform) ) then
         plat = platform
      else
         plat = "PC"
      end if
      ! clean out msg_string
      msg_string = ""
      ! step through the linked list, appending the lines
      first_time = .true.
      cur_line => err%head
      do
         if( .not.associated(cur_line) ) exit
         posn = len_trim(msg_string)
         if( (posn+3) >= len(msg_string) ) exit ! out of space
         posn = posn + 1
         if( .not.first_time) then
            select case(plat)
               case("UNIX")
                  ! Separate lines with LF
                  msg_string(posn:) = achar(10)
                  posn = posn + 1
               case("MAC")
                  ! Separate lines with CR
                  msg_string(posn:) = achar(13)
                  posn = posn + 1
               case default
                  msg_string(posn:) = achar(13) // achar(10)
                  posn = posn + 2
            end select
         end if
         msg_string(posn:) = trim(cur_line%line)
         first_time = .false.
         cur_line => cur_line%next
      end do
   end subroutine err_get_msgs
   !##################################################################
end module error_handler
!#####################################################################

max_stars_repo_name: physicsfan/Object-Oriented-GRASP

max_issues_repo_name: physicsfan/Object-Oriented-GRASP

max_forks_repo_name: physicsfan/Object-Oriented-GRASP


text: #redirect Think Tank Creative Dept.


text:
module call_stack_test
    use erloff_call_stack_m, only: call_stack_t
    use iso_varying_string, only: operator(//)
    use erloff_module_m, only: module_t
    use erloff_procedure_m, only: procedure_t
    use veggies, only: &
            Result_t, &
            Test_Item_t, &
            assert_Includes, &
            assert_Not, &
            assert_That, &
            describe, &
            it

    implicit none
    private

    public :: test_call_stack
contains
    function test_call_stack() result(tests)
        type(Test_Item_t) :: tests

        type(Test_Item_t) :: individual_tests(5)

        individual_tests(1) = it( &
                "Can tell if it originated from a module", &
                checkOriginatedFromModule)
        individual_tests(2) = it( &
                "Can tell if it originated from a procedure", &
                checkOriginatedFromProcedure)
        individual_tests(3) = it( &
                "Can tell if it contains a module", &
                checkContainsModule)
        individual_tests(4) = it( &
                "Can tell if it contains a procedure", &
                checkContainsProcedure)
        individual_tests(5) = it( &
                "String includes the given module and procedure names", &
                checkStringIncludes)
        tests = describe("CallStack_t", individual_tests)
    end function test_call_stack

    pure function checkOriginatedFromModule() result(result_)
        type(Result_t) :: result_

        type(Module_t) :: another_module
        type(procedure_t) :: another_procedure
        type(Module_t) :: the_module
        type(procedure_t) :: the_procedure
        type(call_stack_t) :: stack

        another_module = module_t("Another_m")
        another_procedure = procedure_t("another")
        the_module = module_t("Some_m")
        the_procedure = procedure_t("some")
        stack = call_stack_t(the_module, the_procedure)
        stack = stack%with_names_prepended(another_module, another_procedure)

        result_ = &
                assert_That( &
                        stack.originatedFrom.the_module, &
                        stack%repr() // '.originatedFrom.' // the_module%repr()) &
                .and.assert_Not( &
                        stack.originatedFrom.another_module, &
                        stack%repr() // '.originatedFrom.' // another_module%repr())
    end function checkOriginatedFromModule

    pure function checkOriginatedFromProcedure() result(result_)
        type(Result_t) :: result_

        type(Module_t) :: another_module
        type(procedure_t) :: another_procedure
        type(Module_t) :: the_module
        type(procedure_t) :: the_procedure
        type(call_stack_t) :: stack

        another_module = module_t("Another_m")
        another_procedure = procedure_t("another")
        the_module = module_t("Some_m")
        the_procedure = procedure_t("some")
        stack = call_stack_t(the_module, the_procedure)
        stack = stack%with_names_prepended(another_module, another_procedure)

        result_ = &
                assert_That( &
                        stack.originatedFrom.the_procedure, &
                        stack%repr() // '.originatedFrom.' // the_procedure%repr()) &
                .and.assert_Not( &
                        stack.originatedFrom.another_procedure, &
                        stack%repr() // '.originatedFrom.' // another_procedure%repr())
    end function checkOriginatedFromProcedure

    pure function checkContainsModule() result(result_)
        type(Result_t) :: result_

        type(Module_t) :: another_module
        type(procedure_t) :: another_procedure
        type(Module_t) :: other_module
        type(Module_t) :: the_module
        type(procedure_t) :: the_procedure
        type(call_stack_t) :: stack

        another_module = module_t("Another_m")
        another_procedure = procedure_t("another")
        other_module = module_t("Other_m")
        the_module = module_t("Some_m")
        the_procedure = procedure_t("some")
        stack = call_stack_t(the_module, the_procedure)
        stack = stack%with_names_prepended(another_module, another_procedure)

        result_ = &
                assert_That( &
                        stack.includes.the_module, &
                        stack%repr() // '.includes.' // the_module%repr()) &
                .and.assert_Not( &
                        stack.includes.other_module, &
                        stack%repr() // '.includes.' // other_module%repr())
    end function checkContainsModule

    pure function checkContainsProcedure() result(result_)
        type(Result_t) :: result_

        type(Module_t) :: another_module
        type(procedure_t) :: another_procedure
        type(procedure_t) :: other_procedure
        type(Module_t) :: the_module
        type(procedure_t) :: the_procedure
        type(call_stack_t) :: stack

        another_module = module_t("Another_m")
        another_procedure = procedure_t("another")
        other_procedure = procedure_t("other")
        the_module = module_t("Some_m")
        the_procedure = procedure_t("some")
        stack = call_stack_t(the_module, the_procedure)
        stack = stack%with_names_prepended(another_module, another_procedure)

        result_ = &
                assert_That( &
                        stack.includes.the_procedure, &
                        stack%repr() // '.includes.' // the_procedure%repr()) &
                .and.assert_Not( &
                        stack.includes.other_procedure, &
                        stack%repr() // '.includes.' // other_procedure%repr())
    end function checkContainsProcedure

    pure function checkStringIncludes() result(result_)
        type(Result_t) :: result_

        type(Module_t) :: another_module
        type(procedure_t) :: another_procedure
        type(Module_t) :: the_module
        type(procedure_t) :: the_procedure
        type(call_stack_t) :: stack

        another_module = module_t("Another_m")
        another_procedure = procedure_t("another")
        the_module = module_t("Some_m")
        the_procedure = procedure_t("some")
        stack = call_stack_t(the_module, the_procedure)
        stack = stack%with_names_prepended(another_module, another_procedure)

        result_ = &
                assert_Includes(the_module%to_string(), stack%to_string()) &
                .and.assert_Includes(the_procedure%to_string(), stack%to_string()) &
                .and.assert_Includes(another_module%to_string(), stack%to_string()) &
                .and.assert_Includes(another_procedure%to_string(), stack%to_string())
    end function checkStringIncludes
end module call_stack_test

max_stars_repo_name: everythingfunctional/erloff

max_issues_repo_name: everythingfunctional/erloff

max_forks_repo_name: everythingfunctional/erloff


text:
! Copyright (c) 2006-2011 IDRIS/CNRS
! Author: Philippe Wautelet (IDRIS/CNRS), wautelet@idris.fr
! Distributed under the CeCILL 2.0 license. For full terms see the file LICENSE.

#ifndef WITHOUTMPI

! This subroutine is executed by every process
subroutine init_comm
  use amr_commons
  use mpi
  implicit none

  integer::color_iogroup,color_comp,i,ierr
  integer::nbgbl,szbgbl,szsmbl

  ! Check if ncpus are coherent
  if(ncpu_io==0)then
     if(ncpu==0)then
        ncpu = ncpu_world
     else
        ncpu_io = ncpu_world-ncpu
     end if
  else
     if(ncpu==0)then
        ncpu = ncpu_world-ncpu_io
     end if
  end if

#ifndef IOPROC
  if(ncpu_io/=0) then
     if(myid_world==1)write(*,*)'Error:'
     if(myid_world==1)write(*,*)'Compiled without specialized I/O processes support'
     call clean_stop
  end if
#endif
  if(ncpu_io/=0.AND.(fileformat_in/='std'.OR.fileformat_out/='std')) then
    if(myid_world==1)write(*,*)'Error:'
    if(myid_world==1)write(*,*)'Specialized I/O processes only with std fileformats'
    call clean_stop
  end if

  if(ncpu+ncpu_io/=ncpu_world.OR.ncpu<0.OR.ncpu_io<0)then
     if(myid_world==1)write(*,*)'Error in the namelist:'
     if(myid_world==1)write(*,*)'ncpu/ncpu_io distribution is incoherent'
     call clean_stop
  end if

  if(ncpu_io>ncpu)then
     if(myid_world==1)write(*,*)'Error in the namelist:'
     if(myid_world==1)write(*,*)'ncpu should be at least ncpu_io'
     call clean_stop
  end if

  ! Distribute CPUs
  if(ncpu_io/=0)then
     if(mod(ncpu_world,ncpu_io)==0)then
        nbgbl  = ncpu_io
        szbgbl = ncpu_world/ncpu_io
        szsmbl = 0
     else
        nbgbl  = mod(ncpu_world,ncpu_io)
        szbgbl = ncpu_world/ncpu_io+1
        szsmbl = ncpu_world/ncpu_io
     end if
     if(myid_world<=nbgbl*szbgbl)then
        color_iogroup = (myid_world-1)/szbgbl
        if(mod(myid_world,szbgbl)==1)then
           color_comp = 0
        else
           color_comp = 1
        end if
     else
        color_iogroup = (myid_world-1-nbgbl*szbgbl)/szsmbl+nbgbl
        if(mod(myid_world-nbgbl*szbgbl,szsmbl)==1)then
           color_comp = 0
        else
           color_comp = 1
        end if
     end if
  else
     ! The consequence of this is that MPI_COMM_COMP=MPI_COMM_IOGROUP=MPI_COMM_WORLD
     color_iogroup = 0
     color_comp    = 1
  end if

  ! Initialize communicators
  call MPI_COMM_SPLIT(MPI_COMM_WORLD,color_iogroup,myid_world,MPI_COMM_IOGROUP,ierr)
  call MPI_COMM_RANK(MPI_COMM_IOGROUP,myid_iogroup,ierr)
  call MPI_COMM_SIZE(MPI_COMM_IOGROUP,ncpu_iogroup,ierr)
  if(ncpu_io>0)then
     myid_iogroup = myid_iogroup+1
  else
     ! Necessary to prevent problems later (a wrong branch can be taken by a process)
     myid_iogroup=MPI_UNDEFINED
  end if

  call MPI_COMM_SPLIT(MPI_COMM_WORLD,color_comp,myid_world,MPI_COMM_COMP,ierr)
  call MPI_COMM_RANK(MPI_COMM_COMP,myid,ierr)
  if(color_comp==1)then
     myid    = myid+1
     myid_io = MPI_UNDEFINED
  else
     myid_io = myid+1
     myid    = MPI_UNDEFINED
     allocate(iogroup2world(0:ncpu_iogroup-1),iogroup2comp(0:ncpu_iogroup-1))
     do i=0,ncpu_iogroup-1
        iogroup2world(i)=myid_world-1+i
     end do
     iogroup2comp(0)=MPI_UNDEFINED
     do i=1,ncpu_iogroup-1
        iogroup2comp(i)=myid_world-myid_io+i
     end do
  end if

  ncoarse = nx*ny*nz

  my_iogroup=color_iogroup+1

  ncpu = ncpu
end subroutine init_comm

#endif

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rknop/amuse


!
! Copyright (C) by Argonne National Laboratory
!     See COPYRIGHT in top-level directory
!

subroutine MPI_Add_error_string_f08(errorcode, string, ierror)
    use, intrinsic :: iso_c_binding, only : c_int, c_char
    use :: mpi_c_interface, only : c_Datatype, c_Comm, c_Request
    use :: mpi_c_interface, only : MPIR_Add_error_string_c
    use :: mpi_c_interface, only : MPIR_Fortran_string_f2c

    implicit none

    integer, intent(in) :: errorcode
    character(len=*), intent(in) :: string
    integer, optional, intent(out) :: ierror

    integer(c_int) :: errorcode_c
    character(kind=c_char) :: string_c(len_trim(string)+1)
    integer(c_int) :: ierror_c

    call MPIR_Fortran_string_f2c(string, string_c)

    if (c_int == kind(0)) then
        ierror_c = MPIR_Add_error_string_c(errorcode, string_c)
    else
        errorcode_c = errorcode
        ierror_c = MPIR_Add_error_string_c(errorcode_c, string_c)
    end if

    if (present(ierror)) ierror = ierror_c

end subroutine MPI_Add_error_string_f08