text:
C Copyright(C) 2009-2017 National Technology & Engineering Solutions of
C Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains certain rights in this software.
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions are
C met:
C
C     * Redistributions of source code must retain the above copyright
C       notice, this list of conditions and the following disclaimer.
C
C     * Redistributions in binary form must reproduce the above
C       copyright notice, this list of conditions and the following
C       disclaimer in the documentation and/or other materials provided
C       with the distribution.
C     * Neither the name of NTESS nor the names of its
C       contributors may be used to endorse or promote products derived
C       from this software without specific prior written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

C=======================================================================
      SUBROUTINE TPLABV (NPLT, IVAR, NAME, NE, LABSTR, MAPEL, MAPND)
C=======================================================================

C   --*** TPLABV *** (TPLOT) Get a plot label
C   --   Written by Amy Gilkey - revised 11/03/87
C   --
C   --TPLABV returns a plot label.
C   --
C   --Parameters:
C   --   NPLT - IN - the label type:
C   --      -1 = name and number (long form)
C   --       0 = name and number (short form)
C   --   IVAR - IN - the variable index which determines type
C   --   NAME - IN - the variable name
C   --   NE - IN - the variable number
C   --   LABSTR - OUT - the label string
C   --
C   --Common Variables:
C   --   Uses NVARNP, NVAREL of /DBNUMS/

      CHARACTER*(*) NAME
      CHARACTER*(*) LABSTR
      INTEGER MAPEL(*), MAPND(*)

      CHARACTER TYP

      CALL DBVTYP_BL (IVAR, TYP, IDUM)

      IF (NPLT .LE. -1) THEN
         IF ((TYP .EQ. 'H') .OR. (TYP .EQ. 'G')) THEN
            WRITE (LABSTR, 10000, IOSTAT=IDUM) NAME
         ELSE IF (TYP .EQ. 'N') THEN
           WRITE (LABSTR, 10000, IOSTAT=IDUM) NAME, 'at NODE',
     *       MAPND(NE)
         ELSE IF (TYP .EQ. 'E') THEN
            WRITE (LABSTR, 10000, IOSTAT=IDUM) NAME, 'at ELEMENT',
     *       MAPEL(NE)
         END IF

      ELSE

         IF ((TYP .EQ. 'H') .OR. (TYP .EQ. 'G')) THEN
            WRITE (LABSTR, 10000, IOSTAT=IDUM) NAME
         ELSE IF (TYP .EQ. 'N') THEN
            WRITE (LABSTR, 10000, IOSTAT=IDUM) NAME, 'NODE',
     *       MAPND(NE)
         ELSE IF (TYP .EQ. 'E') THEN
           WRITE (LABSTR, 10000, IOSTAT=IDUM) NAME, 'ELEM',
     *       MAPEL(NE)
         END IF
      END IF

      CALL SQZSTR (LABSTR, LSTR)

      RETURN
10000  FORMAT (A, :, ' ', A, I12)
      END

max_stars_repo_name: brian-kelley/seacas

max_issues_repo_name: brian-kelley/seacas

max_forks_repo_name: brian-kelley/seacas


text:
      MODULE set_I
      INTERFACE
!...Generated by Pacific-Sierra Research 77to90  4.4G  08:35:52  03/09/06
      SUBROUTINE set (S1, S2, NA, NB, RAB, II)
      USE vast_kind_param,ONLY: DOUBLE
      REAL(DOUBLE), INTENT(IN) :: S1
      REAL(DOUBLE), INTENT(IN) :: S2
      INTEGER, INTENT(IN) :: NA
      INTEGER, INTENT(IN) :: NB
      REAL(DOUBLE), INTENT(IN) :: RAB
      INTEGER, INTENT(IN) :: II
      END SUBROUTINE
      END INTERFACE
      END MODULE

max_stars_repo_name: openmopac/MOPAC-archive

max_issues_repo_name: openmopac/MOPAC-archive

max_forks_repo_name: openmopac/MOPAC-archive


text:
The Best Western Palm Court Hotel is a hotels hotel in downtown Downtown Davis.  It has been rated four diamonds by AAA. Rooms cost between $119 and $199 per night. Check in time is 2pm, check out time is 11am, and any cancellations to reservations must be made by 2pm 24 hours prior to scheduled arrival.

Payment can be made by all major credit cards, with discounts available for AAA and AARP members. Wifi is free and unencrypted throughout the hotel, with the best reception on the second floor.

There are slight variations in rooms of the same type and price:

  Room 207 is a King Suite, and has two floors.  It is priced the same as a normal King Suite.  It is also directly above Cafe Bernardo Bar Bernardo, so watch out on Thursdays and crazy weekends.
  Room 214 is a Double Queen Suite.  However, it has jacuzzi jets in the bath tub, unlike the other Double Queens.
  Room 101 is smaller than other Queen Suites, but has a handicapped accessible bathroom.
  Room 302 is the quietest in the hotel.  However, some guests have had to be moved because they felt it had bad http://en.wikipedia.org/wiki/Fung_shui Feng Shui.  This has happened more than once.
  The exercise room has a jacuzzi, sauna, and a big rooftop balcony that would be perfect for watching parades or other downtown events.

Davis has another Best Western called the Best Western University Lodge.  Best Westerns are independently owned and operated franchises. Best Western operates as a nonprofit membership association.

   Looking for a place to stay (or for a guest to stay)?  Davis has many hotels hotel options!

20070830 22:57:05 nbsp This hotel is more elegant than most Best Westerns.  In fact, for the price, it was one of the nicest hotels Ive ever stayed in.  University departments will put up prospective faculty here when they are trying to woo them.  The location is right in the middle of town, so convenient to just about everything, and you can practically fall out of your room into Cafe Bernardo Café Bernardo, restaurantside or counterside. Users/CovertProfessor

20091212 05:43:48 nbsp WARNING: INHOUSE THEFT! My family was subject to a room break in, with subsequent auto theft in the supposedly secure garage. The security cameras were nonfunctioning when this happened (Dec 2009)  clearly an inside job. Despite the fact that a camera, PDA and rental car were stolen, the management assumed no responsibility and did not even make the gesture of offering a refund of the cost of the stay.

 Users/Alucky

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


module global_variables
  implicit none
  real(8),parameter :: pi = 4d0*atan(1d0)
  complex(8),parameter :: zi = (0d0,1d0)
  real(8),parameter :: fs = 0.024189d0
  real(8),parameter :: ev = 27.2114d0
  real(8),parameter :: b_a = 0.529d0
  real(8),parameter :: clight = 137.035999139d0

! material parameters
  real(8),parameter :: tau_z = 1d0, delta_gap = 0d0
! Fermi velocity of graphene PRL 101, 226405 (2008).
  real(8),parameter :: velocity = clight*1.12d6/299792458d0
  integer,parameter :: nkx = 32, nky = 32
  real(8),parameter :: kx_max = 1d0, ky_max = 1d0
! Fermi-Dirac distribution
  real(8),parameter :: mu_F = 0d0/ev
  real(8),parameter :: kbT  = 80d0/11604.505d0/ev

  real(8) :: kx(nkx),ky(nky)
  real(8) :: dkx,dky
  complex(8) :: zpsi(2,2,nkx,nky)
  real(8) :: occ(2,nkx,nky), eps(2,nkx,nky)

! time propagation
  real(8),parameter :: Tprop = 2.2d3/fs, dt = 0.08d0
  integer,parameter :: nt = aint(Tprop/dt)+1

  real(8) :: jt(2,0:nt), ac(2,0:nt), ac_dt2(2,0:nt)

end module global_variables
!----------------------------------------------------------------------------------------!
program main
  use global_variables
  implicit none

  call init
  call init_ac

  call time_propagation

end program main
!----------------------------------------------------------------------------------------!
subroutine init
  use global_variables
  implicit none
  integer :: ikx,iky
  real(8) :: kxt, kyt, theta
  complex(8) :: zs

  do ikx = 1,nkx
    kx(ikx) = -kx_max + 2d0*kx_max*dble(ikx-1)/dble(nkx-1)
  end do
  dkx = 2d0*kx_max/dble(nkx-1)

  do iky = 1,nky
    ky(iky) = -ky_max + 2d0*ky_max*dble(iky-1)/dble(nky-1)
  end do
  dky = 2d0*ky_max/dble(nky-1)

  eps = 0d0

  if(delta_gap == 0d0)then ! Dirac cone
    do ikx = 1, nkx
      kxt = kx(ikx)
      do iky = 1,nky
        kyt = ky(iky)

        zs = tau_z*kxt+zI*kyt
        if(zs /= 0d0)then
          theta = -aimag(log(-zs))
        else
          theta = 0d0
        end if

        zpsi(1,1,ikx,iky) = exp(zI*theta)/sqrt(2d0)
        zpsi(2,1,ikx,iky) = 1d0/sqrt(2d0)
        eps(1,ikx,iky) = -velocity*sqrt(tau_z**2*kxt**2 + kyt**2)

        zpsi(1,2,ikx,iky) = exp(zI*(theta+pi))/sqrt(2d0)
        zpsi(2,2,ikx,iky) = 1d0/sqrt(2d0)
        eps(2,ikx,iky) = velocity*sqrt(tau_z**2*kxt**2 + kyt**2)

      end do
    end do
  else
    stop 'Kane band is not implemented yet.'
  end if

! Fermi-Dirac distribution
  do ikx = 1,nkx
    do iky = 1,nky

      occ(1,ikx,iky) = 2d0/(exp((eps(1,ikx,iky)-mu_F)/kbT)+1d0)
      occ(2,ikx,iky) = 2d0/(exp((eps(2,ikx,iky)-mu_F)/kbT)+1d0)

    end do
  end do

end subroutine init
!----------------------------------------------------------------------------------------!
subroutine time_propagation
  use global_variables
  implicit none
  real(8) :: jxy(2)
  integer :: it

  open(21,file="current.out")
  call current(jxy(:),0)
  write(21,"(999e26.16e3)")0d0,ac(:,0),jxy
  do it = 0,nt-1
    if(mod(it,max(100,nt/100)) == 0)write(*,*)"it=",it,nt
    call dt_evolve(it)
    call current(jxy(:),it)
    write(21,"(999e26.16e3)")dt*(it+1),ac(:,it+1),jxy
  end do
  close(21)

end subroutine time_propagation
!----------------------------------------------------------------------------------------!
subroutine dt_evolve(it)
  implicit none
  integer,intent(inout) :: it

  call dt_evolve_Taylor(it)
!  call dt_evolve_Magnus(it)
!  call dt_evolve_Magnus_1st(it)

end subroutine dt_evolve
!!----------------------------------------------------------------------------------------!
!subroutine dt_evolve_Magnus(it)
!  use global_variables
!  implicit none
!  integer,intent(in) :: it
!  real(8) :: acx0, acx1, acx2
!  real(8) :: acy0, acy1, acy2
!  integer :: ikx, iky
!  real(8) :: kxt, kyt,kxt_eff, kyt_eff
!  real(8) :: delta,lambda(2)
!  complex(8) :: zHeff(2,2),zalpha,zx,zvec(2,2),zc(2)
!  real(8) :: const,ss
!
!
!  if(delta_gap /= 0d0)stop 'delta_gap has to be zero in Magnus propagator.'
!
!  acx0 = ac(1,it)
!  acx2 = (ac(1,it)-2d0*ac_dt2(1,it)+ac(1,it+1))/(0.5d0*dt)**2
!  acx1 = (ac(1,it+1)-acx0-0.5d0*dt**2*acx2)/dt
!
!  acy0 = ac(2,it)
!  acy2 = (ac(2,it)-2d0*ac_dt2(2,it)+ac(2,it+1))/(0.5d0*dt)**2
!  acy1 = (ac(2,it+1)-acy0-0.5d0*dt**2*acy2)/dt
!
!  const = -dt**2/6d0*tau_z*velocity**2
!
!
!!$omp parallel default(shared), private(ikx,iky,kxt,kyt,kxt_eff,kyt_eff,ss,delta,zalpha, &
!!$omp & lambda,zx,zvec,zc)
!!$omp do collapse(2)
!  do ikx = 1,nkx
!    do iky = 1,nky
!      kxt = kx(ikx) + acx0
!      kxt_eff = kxt + 0.5d0*dt*acx1 + dt**2/6d0*acx2
!      kyt = ky(iky) + acy0
!      kyt_eff = kyt + 0.5d0*dt*acy1 + dt**2/6d0*acy2
!
!      ss = kxt*acy1-acx1*kyt
!!      zHeff(1,1) =  const*ss
!!      zHeff(2,1) =  velocity*(tau_z*kxt_eff + zI*kyt_eff)
!!      zHeff(1,2) =  velocity*(tau_z*kxt_eff - zI*kyt_eff)
!!      zHeff(2,2) = -const*ss
!
!      delta  = const*ss
!      zalpha = velocity*(tau_z*kxt_eff + zI*kyt_eff)
!!
!      if(delta >= 0d0)then
!
!! vector 1
!        lambda(1) = sqrt(delta**2 + abs(zalpha)**2)
!        zx = zalpha/(lambda(1) + delta)
!        ss = 1d0/sqrt(1d0 + abs(zx)**2)
!        zvec(1,1) = ss
!        zvec(2,1) = ss*zx
!
!! vector 2
!        lambda(2) = -lambda(1)
!        zvec(1,2) = -ss*conjg(zx)
!        zvec(2,2) =  ss
!
!      else
!
!! vector 1
!        lambda(1) = sqrt(delta**2 + abs(zalpha)**2)
!        zx = conjg(zalpha)/(lambda(1) - delta)
!        ss = 1d0/sqrt(1d0 + abs(zx)**2)
!        zvec(1,1) = ss*zx
!        zvec(2,1) = ss
!
!! vector 2
!        lambda(2) = -lambda(1)
!        zvec(1,2) =  ss
!        zvec(2,2) = -ss*conjg(zx)
!
!      end if
!
!! projection
!      zc(1) = conjg(zvec(1,1))*zpsi(1,ikx,iky) + conjg(zvec(2,1))*zpsi(2,ikx,iky)
!      zc(2) = conjg(zvec(1,2))*zpsi(1,ikx,iky) + conjg(zvec(2,2))*zpsi(2,ikx,iky)
!
!! propagation
!      zc(1) = zc(1)*exp(-zI*lambda(1)*dt)
!      zc(2) = zc(2)*exp(-zI*lambda(2)*dt)
!
!! update wavefunction
!      zpsi(:,ikx,iky) = zc(1)*zvec(:,1) + zc(2)*zvec(:,2)
!
!
!    end do
!  end do
!!$omp end do
!!$omp end parallel
!
!end subroutine dt_evolve_Magnus
!!----------------------------------------------------------------------------------------!
!subroutine dt_evolve_Magnus_1st(it)
!  use global_variables
!  implicit none
!  integer,intent(in) :: it
!  real(8) :: acx0, acx1, acx2
!  real(8) :: acy0, acy1, acy2
!  integer :: ikx, iky
!  real(8) :: kxt, kyt,kxt_eff, kyt_eff
!  real(8) :: delta,lambda(2)
!  complex(8) :: zHeff(2,2),zalpha,zx,zvec(2,2),zc(2)
!  real(8) :: const,ss
!
!
!  if(delta_gap /= 0d0)stop 'delta_gap has to be zero in Magnus propagator.'
!
!  acx0 = ac(1,it)
!  acx2 = (ac(1,it)-2d0*ac_dt2(1,it)+ac(1,it+1))/(0.5d0*dt)**2
!  acx1 = (ac(1,it+1)-acx0-0.5d0*dt**2*acx2)/dt
!
!  acy0 = ac(2,it)
!  acy2 = (ac(2,it)-2d0*ac_dt2(2,it)+ac(2,it+1))/(0.5d0*dt)**2
!  acy1 = (ac(2,it+1)-acy0-0.5d0*dt**2*acy2)/dt
!
!  delta = 0d0
!
!
!!$omp parallel default(shared), private(ikx,iky,kxt,kyt,kxt_eff,kyt_eff,ss,zalpha, &
!!$omp & lambda,zx,zvec,zc)
!!$omp do collapse(2)
!  do ikx = 1,nkx
!    do iky = 1,nky
!      kxt = kx(ikx) + acx0
!      kxt_eff = kxt + 0.5d0*dt*acx1 + dt**2/6d0*acx2
!      kyt = ky(iky) + acy0
!      kyt_eff = kyt + 0.5d0*dt*acy1 + dt**2/6d0*acy2
!
!      zalpha = velocity*(tau_z*kxt_eff + zI*kyt_eff)
!!
!! vector 1
!      lambda(1) = sqrt(delta**2 + abs(zalpha)**2)
!      zx = zalpha/(lambda(1) + delta)
!      ss = 1d0/sqrt(1d0 + abs(zx)**2)
!      zvec(1,1) = ss
!      zvec(2,1) = ss*zx
!
!! vector 2
!      lambda(2) = -lambda(1)
!      zvec(1,2) = -ss*conjg(zx)
!      zvec(2,2) =  ss
!
!
!! projection
!      zc(1) = conjg(zvec(1,1))*zpsi(1,ikx,iky) + conjg(zvec(2,1))*zpsi(2,ikx,iky)
!      zc(2) = conjg(zvec(1,2))*zpsi(1,ikx,iky) + conjg(zvec(2,2))*zpsi(2,ikx,iky)
!
!! propagation
!      zc(1) = zc(1)*exp(-zI*lambda(1)*dt)
!      zc(2) = zc(2)*exp(-zI*lambda(2)*dt)
!
!! update wavefunction
!      zpsi(:,ikx,iky) = zc(1)*zvec(:,1) + zc(2)*zvec(:,2)
!
!
!    end do
!  end do
!!$omp end do
!!$omp end parallel
!
!end subroutine dt_evolve_Magnus_1st
!!----------------------------------------------------------------------------------------!
subroutine dt_evolve_Taylor(it)
  use global_variables
  implicit none
  integer,intent(in) :: it
  integer :: ikx, iky
  real(8) :: kxt, kyt
  complex(8) :: zpsi_t(2,2,0:4)

!$omp parallel default(shared), private(ikx,iky,kxt,kyt,zpsi_t)
!$omp do collapse(2)
  do ikx = 1,nkx
    do iky = 1,nky

      kxt = kx(ikx) + ac(1,it)
      kyt = ky(iky) + ac(2,it)
! RK1
      zpsi_t(:,:,0) = zpsi(:,:,ikx,iky)
      zpsi_t(1,:,1) = 0.5d0*delta_gap*zpsi_t(1,:,0) &
        + velocity*(tau_z*kxt - zI*kyt)*zpsi_t(2,:,0)
      zpsi_t(2,:,1) = velocity*(tau_z*kxt + zI*kyt)*zpsi_t(1,:,0)  &
        - 0.5d0*delta_gap*zpsi_t(2,:,0)

      kxt = kx(ikx) + ac_dt2(1,it)
      kyt = ky(iky) + ac_dt2(2,it)
! RK2
      zpsi_t(:,:,0) = zpsi(:,:,ikx,iky) - zI*0.5d0*dt*zpsi_t(:,:,1)

      zpsi_t(1,:,2) = 0.5d0*delta_gap*zpsi_t(1,:,0) &
        + velocity*(tau_z*kxt - zI*kyt)*zpsi_t(2,:,0)
      zpsi_t(2,:,2) = velocity*(tau_z*kxt + zI*kyt)*zpsi_t(1,:,0)  &
        - 0.5d0*delta_gap*zpsi_t(2,:,0)

! RK3
      zpsi_t(:,:,0) = zpsi(:,:,ikx,iky) - zI*0.5d0*dt*zpsi_t(:,:,2)

      zpsi_t(1,:,3) = 0.5d0*delta_gap*zpsi_t(1,:,0) &
        + velocity*(tau_z*kxt - zI*kyt)*zpsi_t(2,:,0)
      zpsi_t(2,:,3) = velocity*(tau_z*kxt + zI*kyt)*zpsi_t(1,:,0)  &
        - 0.5d0*delta_gap*zpsi_t(2,:,0)

      kxt = kx(ikx) + ac(1,it+1)
      kyt = ky(iky) + ac(2,it+1)

! RK4
      zpsi_t(:,:,0) = zpsi(:,:,ikx,iky) - zI*dt*zpsi_t(:,:,3)

      zpsi_t(1,:,4) = 0.5d0*delta_gap*zpsi_t(1,:,0) &
        + velocity*(tau_z*kxt - zI*kyt)*zpsi_t(2,:,0)
      zpsi_t(2,:,4) = velocity*(tau_z*kxt + zI*kyt)*zpsi_t(1,:,0)  &
        - 0.5d0*delta_gap*zpsi_t(2,:,0)

      zpsi(:,:,ikx,iky) = zpsi(:,:,ikx,iky) -zI*dt/6d0*( &
        zpsi_t(:,:,1) + 2d0*zpsi_t(:,:,2) + 2d0*zpsi_t(:,:,3) + zpsi_t(:,:,4) )

    end do
  end do
!$omp end do
!$omp end parallel

end subroutine dt_evolve_Taylor
!----------------------------------------------------------------------------------------!
subroutine init_ac
  use global_variables
  implicit none
  integer :: it
  real(8) :: tt,tt2, xx, xx2
  real(8) :: E0,omega,tpulse
  real(8) :: E_SD, T_SD
  real(8) :: gap_Floquet

!  gap_Floquet = 0.1d0/ev
  E0 = 1d7*b_a*1d-10/ev
  omega = 0.19074d0/ev ! 6.5 micron
  tpulse = 1d3/fs

!  E0 = omega*0.5d0/velocity*sqrt( &
!    (gap_Floquet + omega)**2 - omega**2 &
!    )

  write(*,"(A,2x,e26.16e3,A)")"Field strength =",E0*ev/b_a*1d10, "eV/m"

! Source-drain
  E_SD = 1d-8
  T_SD = 10d0/fs

  ac = 0d0
  ac_dt2 = 0d0

  do it = 0,nt

    tt = dt*it
    tt2 = dt*it + 0.5d0*dt

! source-drain
    if(tt<T_SD)then
      ac(1,it) = ac(1,it) - E_SD*(tt/2d0 - T_SD/(2d0*pi)*sin(pi*tt/T_SD))
    else
      ac(1,it) = ac(1,it) - E_SD*(tt - T_SD/2d0)
    end if
    if(tt2<T_SD)then
      ac_dt2(1,it) = ac_dt2(1,it) - E_SD*(tt2/2d0 - T_SD/(2d0*pi)*sin(pi*tt2/T_SD))
    else
      ac_dt2(1,it) = ac_dt2(1,it) - E_SD*(tt2 - T_SD/2d0)
    end if

! circular laser
    xx = tt - T_SD - 0.5d0*tpulse
    xx2 = tt2 - T_SD - 0.5d0*tpulse

     if( abs(xx)<0.5d0*tpulse )then
        ac(1,it) = ac(1,it) - E0/omega*cos(pi*xx/tpulse)**4*sin(omega*xx)
        ac(2,it) = ac(2,it) - E0/omega*cos(pi*xx/tpulse)**4*cos(omega*xx)
     end if

     if( abs(xx2)<0.5d0*tpulse )then
        ac_dt2(1,it) = ac_dt2(1,it) - E0/omega*cos(pi*xx2/tpulse)**4*sin(omega*xx2)
        ac_dt2(2,it) = ac_dt2(2,it) - E0/omega*cos(pi*xx2/tpulse)**4*cos(omega*xx2)
     end if

  end do

end subroutine init_ac
!----------------------------------------------------------------------------------------!
subroutine current(jxy,it)
  use global_variables
  implicit none
  real(8),intent(out) :: jxy(2)
  integer,intent(in) :: it
  real(8) :: jx,jy
  integer :: ikx, iky
  real(8) :: kxt, kyt
  real(8) :: jxt,jyt,jx0,jy0,xx
  complex(8) :: zs

  jx = 0d0
  jy = 0d0

!$omp parallel default(shared), private(ikx,iky,kxt,kyt,jxt,jyt,zs,jx0,jy0)
!$omp do reduction(+:jx,jy) collapse(2)
  do ikx = 1,nkx
    do iky = 1,nky

       kxt = kx(ikx) + ac(1,it)
       kyt = ky(iky) + ac(2,it)

       jxt = 2d0*occ(1,ikx,iky)*real(zpsi(1,1,ikx,iky)*conjg(zpsi(2,1,ikx,iky))) &
            +2d0*occ(2,ikx,iky)*real(zpsi(1,2,ikx,iky)*conjg(zpsi(2,2,ikx,iky)))
       jyt = 2d0*occ(1,ikx,iky)*real(zI*zpsi(1,1,ikx,iky)*conjg(zpsi(2,1,ikx,iky))) &
            +2d0*occ(2,ikx,iky)*real(zI*zpsi(1,2,ikx,iky)*conjg(zpsi(2,2,ikx,iky)))

       zs =  -(tau_z*kxt-zI*kyt)&
            /sqrt(kxt**2+kyt**2)

       jx0 =  2d0*0.5d0*real(zs)   *(occ(1,ikx,iky)-occ(2,ikx,iky))
       jy0 =  2d0*0.5d0*real(zI*zs)*(occ(1,ikx,iky)-occ(2,ikx,iky))

       jx = jx + jxt - jx0
       jy = jy + jyt - jy0

     end do
  end do
!$omp end do
!$omp end parallel

  jxy(1) = jx*dkx*dky
  jxy(2) = jy*dkx*dky

end subroutine current
!----------------------------------------------------------------------------------------!
!----------------------------------------------------------------------------------------!
!----------------------------------------------------------------------------------------!


text:
## ||<#E0E0FF>Established||
## || Date, i.e. 2014MMDD ||
## ||<#E0E0FF>Payment Methods||
## || Type of payment ||

Cork it Again Wine Sellers is a discount wine outlet and tasting room.  All wine is $5.99 a bottle.

Current Wine Listing

  Cabernet Sauvignon
  Chardonnay
  Pinot Noir
  Syrah
  Valdiguie


!---------------------------------------------------------------------
!---------------------------------------------------------------------
! snh: Bifurcations of global reinjection orbits near a saddle-node
!      Hopf bifurcation: cycle-to-point connections
!
! Parameters:
!   PAR(1) : nu1
!   PAR(2) : nu2
!   PAR(3) : d
!   PAR(4) : codimension of connection (0 or 1)
!
!   PAR(5) : delta: distance from cycle to start connection
!   PAR(6) : eps: distance from end connection to point
!   PAR(7) : mu: log of the Floquet multiplier
!   PAR(8) : h: norm of eigenfunction for cycle at 0
!   PAR(9) : T^-: time for connection from cycle (U(7:9)) to section
!   PAR(10) : T^+: time for connection from section (U(10:12)) to point
!   PAR(11) : T: period of the cycle
!   PAR(21) : sigma-: U1(9)-pi (phi-distance W^u(P) from section phi=pi)
!                     ( distance to W^u(b) for codim 0 connection)
!   PAR(22) : sigma+: U0(12)-pi (phi-distance W^s(E) from section phi=pi)
!   PAR(23) : eta: gap size for Lin vector
!   PAR(24) : Z_x: Lin vector (x-coordinate)
!   PAR(25) : Z_y: Lin vector (y-coordinate)
!   PAR(26) : Z_z: Lin vector (z-coordinate)
!---------------------------------------------------------------------
!---------------------------------------------------------------------

  SUBROUTINE RHS(U,PAR,F,JAC,A)

    IMPLICIT NONE
    DOUBLE PRECISION, INTENT(IN) :: U(3), PAR(*)
    LOGICAL, INTENT(IN) :: JAC
    DOUBLE PRECISION, INTENT(OUT) :: F(3), A(3,3)

    DOUBLE PRECISION x, y, phi
    DOUBLE PRECISION nu1, nu2, alpha, beta, omega, c, d, fp, s
    DOUBLE PRECISION sphi, x2y2, cphinu, cphinu2

    nu1 = PAR(1)
    nu2 = PAR(2)
    d = PAR(3)
    fp = 4*atan(1.0d0)*d
    alpha = -1
    s = -1
    c = 0
    omega = 1
    beta = 0

    x = U(1)
    y = U(2)
    phi = U(3)

    sphi = sin(phi)
    x2y2 = x*x + y*y
    cphinu = 2*cos(phi)+nu2
    cphinu2 = cphinu*cphinu

    F(1) = nu1*x - omega*y - (alpha*x - beta*y)*sphi - x2y2*x + d*cphinu2
    F(2) = nu1*y + omega*x - (alpha*y + beta*x)*sphi - x2y2*y + fp*cphinu2
    F(3) = nu2 + s*x2y2 + 2.0*cos(phi) + c*x2y2*x2y2

    IF(JAC)THEN
       A(1,1) = nu1 - alpha*sphi - 3*x**2 - y**2
       A(1,2) = -omega + beta*sphi - 2*x*y
       A(1,3) = -(alpha*x - beta*y)*cos(phi) - d*4*(2*cos(phi)+nu2)*sphi

       A(2,1) = omega - beta*sphi - 2*x*y
       A(2,2) = nu1 - alpha*sphi - 3*y**2 - x*x
       A(2,3) = -(alpha*y + beta*x)*cos(phi) - fp*4*(2*cos(phi) + nu2)*sphi

       A(3,1) = 2*s*x + 4*c*x2y2*x
       A(3,2) = 2*s*y + 4*c*x2y2*y
       A(3,3) = -2*sphi
    ENDIF

  END SUBROUTINE RHS

  SUBROUTINE FUNC(NDIM,U,ICP,PAR,IJAC,F,DFDU,DFDP)
! ---------- ---

    IMPLICIT NONE
    INTEGER, INTENT(IN) :: NDIM, IJAC, ICP(*)
    DOUBLE PRECISION, INTENT(IN) :: U(NDIM), PAR(*)
    DOUBLE PRECISION, INTENT(OUT) :: F(NDIM)
    DOUBLE PRECISION, INTENT(INOUT) :: DFDU(NDIM,NDIM),DFDP(NDIM,*)

    DOUBLE PRECISION T,mu
    DOUBLE PRECISION A(3,3)

    CALL RHS(U,PAR,F,NDIM>3,A)
    IF(NDIM==3)RETURN

    F(4:6) = MATMUL(A,U(4:6))

    T = PAR(11)
    F(1:6) = F(1:6) * T

    ! log of Floquet multiplier in PAR(7)
    mu = PAR(7)
    F(4:6) = F(4:6) - mu*U(4:6)

    IF (NDIM==6) RETURN

    CALL RHS(U(7:9),PAR,F(7:9),.FALSE.,A)

    T = PAR(9)
    F(7:9) = F(7:9) * T

    IF (NDIM==9) RETURN

    CALL RHS(U(10:12),PAR,F(10:12),.FALSE.,A)

    T = PAR(10)
    F(10:12) = F(10:12) * T

  END SUBROUTINE FUNC

  SUBROUTINE STPNT(NDIM,U,PAR,T)
  !--------- -----

    IMPLICIT NONE
    INTEGER, INTENT(IN) :: NDIM
    DOUBLE PRECISION, INTENT(INOUT) :: U(NDIM),PAR(*)
    DOUBLE PRECISION, INTENT(IN) :: T

    INTEGER, PARAMETER :: codim = 1
    DOUBLE PRECISION, PARAMETER :: nu1 = 0, nu2 = -1.46d0, d = 0.01d0
    DOUBLE PRECISION, PARAMETER :: delta = -1d-5, eps = 1d-6
    DOUBLE PRECISION fp(3), ev(3), pi
    DOUBLE PRECISION, SAVE :: s(6)

    IF(NDIM==9)THEN
       IF(T==0)THEN
          s(1:6) = U(1:6)
       ENDIF
       U(7:9) = s(1:3) + PAR(5)*s(4:6)
       RETURN
    ELSEIF(NDIM==12)THEN
       fp = (/0d0, 0d0, acos(-nu2/2)/)
       ev = (/0d0, 0d0, 1d0/)
       U(10:12) = fp(1:3) + PAR(6)*ev(1:3)
       RETURN
    ENDIF

    PAR(1:3) = (/nu1,nu2,d/)
    PAR(4) = codim
    PAR(5) = delta
    PAR(6) = eps
    PAR(21:22) = 0.0

    ! homotopy: equilibrium b and its phase-shifted version, to find
    ! the "heteroclinic" orbit in negative time.
    pi = 4 * ATAN(1D0)
    PAR(11) = -0.1
    PAR(12) = 0
    PAR(13) = 0
    PAR(14) = acos(-PAR(2)/2)
    PAR(15) = 0
    PAR(16) = 0
    PAR(17) = acos(-PAR(2)/2)+2*pi
    PAR(18) = 0.01d0 ! epsilon_0 distance

    ! equilibrium a
    U(1:3) = (/ 0.0d0, 0.0d0, 2*pi-acos(-nu2/2) /)

  END SUBROUTINE STPNT

  SUBROUTINE PVLS(NDIM,U,PAR)
  !--------- ----

    IMPLICIT NONE
    INTEGER, INTENT(IN) :: NDIM
    DOUBLE PRECISION, INTENT(IN) :: U(NDIM)
    DOUBLE PRECISION, INTENT(INOUT) :: PAR(*)

    DOUBLE PRECISION, EXTERNAL :: GETP
    DOUBLE PRECISION d(3), normlv
    INTEGER i, NBC, codim
    DOUBLE PRECISION pi
    LOGICAL, SAVE :: FIRST = .TRUE.

    IF (FIRST) THEN
       FIRST = .FALSE.
       ! initialization for BCND
       pi = 4d0 * ATAN(1d0)
       IF (NDIM==9) THEN
          codim = NINT(PAR(4))
          IF (codim == 0) THEN
             IF (PAR(21)==0) THEN
                ! distance to W^u(b) where b is the phase-shifted equilibrium
                PAR(21) = GETP("BV1",9,U) - PAR(17)
             ENDIF
          ELSE
             PAR(21) = GETP("BV0",9,U) - pi
          ENDIF
       ELSEIF (NDIM==12) THEN
          NBC = AINT(GETP("NBC",0,U))
          IF (NBC==15) THEN
             PAR(22) = GETP("BV0",12,U) - pi
          ELSE
             ! check if Lin vector initialized:
             IF (DOT_PRODUCT(PAR(24:26),PAR(24:26)) > 0) RETURN
             DO i=1,3
                d(i) = GETP("BV1",6+i,U) - GETP("BV0",9+i,U)
             ENDDO
             normlv = sqrt(DOT_PRODUCT(d,d))
             ! gap size in PAR(23)
             PAR(23) = normlv
             ! Lin vector in PAR(24)-PAR(26)
             PAR(24:26) = d(1:3)/normlv
          ENDIF
       ENDIF
    ENDIF
  END SUBROUTINE PVLS

  SUBROUTINE BCND(NDIM,PAR,ICP,NBC,U0,U1,FB,IJAC,DBC)
  !--------- ----

    IMPLICIT NONE
    INTEGER, INTENT(IN) :: NDIM, ICP(*), NBC, IJAC
    DOUBLE PRECISION, INTENT(IN) :: PAR(*), U0(NDIM), U1(NDIM)
    DOUBLE PRECISION, INTENT(OUT) :: FB(NBC)
    DOUBLE PRECISION, INTENT(INOUT) :: DBC(NBC,*)

    DOUBLE PRECISION nu2, delta, eps, fp(3), ev(3), eta
    DOUBLE PRECISION pi
    INTEGER codim

    ! Periodicity boundary conditions on state variables
    FB(1:3) = U0(1:3) - U1(1:3)

    ! Floquet boundary condition
    FB(4:6) = U1(4:6) - U0(4:6)

    ! normalization
    FB(7) = PAR(8) - DOT_PRODUCT(U0(4:6),U0(4:6))
    IF (NBC==7) RETURN

    delta = PAR(5)
    FB(8:10) = U0(7:9) - (U0(1:3) + delta*U0(4:6))
    pi = 4d0 * ATAN(1d0)
    codim = NINT(PAR(4))
    IF(codim == 0)THEN
       ! projection boundary condition for codimension-zero connection
       nu2 = PAR(2)
       FB(11) = COS(U1(9) - PAR(21)) +nu2/2
       RETURN
    ENDIF
    FB(11) = U1(9) - pi - PAR(21)

    IF (NBC==11) RETURN

    nu2 = PAR(2)
    eps = PAR(6)
    IF (abs(nu2) > 2) THEN
       ! truncate nu2 to [-2,2] to avoid floating point exceptions with acos
       nu2 = sign(2d0, nu2)
    ENDIF

    fp = (/0d0, 0d0, acos(-nu2/2)/)
    ev = (/0d0, 0d0, 1d0/)

    FB(12:14) = U1(10:12) - (fp(1:3) + eps*ev(1:3))

    IF (NBC==15) THEN
       FB(15) = U0(12) - pi - PAR(22)
       RETURN
    ENDIF

    eta = PAR(23)
    FB(15:17) = U1(7:9) - U0(10:12) - eta*PAR(24:26)

  END SUBROUTINE BCND

  SUBROUTINE ICND(NDIM,PAR,ICP,NINTS,U,UOLD,UDOT,UPOLD,FI,IJAC,DINT)
  !--------- ----
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: NDIM, ICP(*), NINTS, IJAC
    DOUBLE PRECISION, INTENT(IN) :: PAR(*)
    DOUBLE PRECISION, INTENT(IN) :: U(NDIM), UOLD(NDIM), UDOT(NDIM), UPOLD(NDIM)
    DOUBLE PRECISION, INTENT(OUT) :: FI(NINTS)
    DOUBLE PRECISION, INTENT(INOUT) :: DINT(NINTS,*)

    INTEGER codim

    ! Integral phase condition
    FI(1) = DOT_PRODUCT(U(1:3),UPOLD(1:3))
    IF (NINTS==1) RETURN

    codim = NINT(PAR(4))
    IF (codim == 0) THEN
       ! phase condition for codimension-zero connection
       FI(2) = DOT_PRODUCT(UPOLD(7:9),U(7:9)-UOLD(7:9))
       RETURN
    ENDIF
    FI(2) = DOT_PRODUCT(UPOLD(10:12),U(10:12)-UOLD(10:12))
  END SUBROUTINE ICND

  SUBROUTINE FOPT(NDIM,U,ICP,PAR,IJAC,FS,DFDU,DFDP)
  END SUBROUTINE FOPT


text:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!   EVB-QMDFF - RPMD molecular dynamics and rate constant calculations on
!               black-box generated potential energy surfaces
!
!   Copyright (c) 2021 by Julien Steffen (steffen@pctc.uni-kiel.de)
!                         Stefan Grimme (grimme@thch.uni-bonn.de) (QMDFF code)
!
!   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.
!
!   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.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!
!     subroutine umbrella: Calculate the actual value of the
!       sampling coordinate for a structure and compare it with
!       the reference value(s) for this sampling, then apply
!       bias potential and gradient correction
!
!     part of EVB
!
subroutine umbrella(centroid,xi_ideal,int_ideal,xi_real,dxi_act,grad_xyz,mode)
use general
use evb_mod
implicit none
integer::i,j,i1,i2,j2,k ! loop indices
integer::mode  ! which type of calculation has invoked this subroutine:
               ! 0: usual umbrella sampling
               ! 1: constraint: calculate restrain xi values
               ! 2: pre equilibrate for tri-/tetramolecular reactions
               ! 3: child trajectory: calculate usual xi, no bias applied
real(kind=8)::centroid(3,natoms)  !  the centroid with all the
                                  !  com's of all beads
real(kind=8)::int_ideal(nat6),int_real(nat6),int_tmp(nat6)
real(kind=8)::xi_ideal,xi_real
real(kind=8)::act_coord(3,natoms)
real(kind=8)::v_add,delta
real(kind=8)::pi
real(kind=8)::int_grad(nat6)
real(kind=8)::grad_xyz(3,natoms,nbeads),grad_bias(3,natoms)
real(kind=8)::grad_xyz2(3*natoms),grad_bias2(3*natoms)
!     for pre equilibration addition of bias potentials
real(kind=8)::r_eds(sum_eds,sum_eds) ! the actual distane between two educts
real(kind=8)::Red(3,sum_eds,sum_eds)  ! all possible educt-educt distances
real(kind=8)::com(sum_eds,3)
real(kind=8)::Rinv ! the inverse bond length/distance (for derivatives)
integer::atom  ! the current atom index

!
real(kind=8)::dxi_act(3,natoms) ! the calculated Xi derivative
real(kind=8)::d2xi_act(3,natoms,3,natoms)  ! the calculated Xi hessian
real(kind=8)::fs, fs2, log_fs, coeff1, coeff2, dhams ! for bias potential..
pi=3.1415926535897932384d0

!
!     calculate actual value of xi and compare it with
!     desired value
!
act_coord(1,:)=centroid(1,:)
act_coord(2,:)=centroid(2,:)
act_coord(3,:)=centroid(3,:)
call xyz_2int(act_coord,int_real,natoms)

grad_bias = 0
!
!      call routine to calculate Xi value and its derivatives
!
if (mode .eq. 0) then
   call calc_xi(act_coord,xi_ideal,xi_real,dxi_act,d2xi_act,1)
!
!      if a constrained trajectory is sampled, calculate directly the relative xi value
!      and go back to verlet without applying the bias
!
else if (mode .eq. 1) then
      call calc_xi(act_coord,xi_ideal,xi_real,dxi_act,d2xi_act,2)
!
!      directly return to main routine!
!
   return
!
!      if a child trajectory is sampled, calculate usual Xi values and return directly
!      to main routine
!
else if (mode .eq. 3) then
   call calc_xi(act_coord,xi_ideal,xi_real,dxi_act,d2xi_act,1)
   return
end if
!
!      apply umbrella potential to gradient
!
delta=(xi_real-xi_ideal)
v_add=0.5*k_force*delta*delta
do i=1,nbeads
   grad_xyz(:,:,i)=grad_xyz(:,:,i)+k_force*delta*dxi_act
end do
!
!      apply bias potential (with hessian) to gradient (why???)
!      ---> temporarily deactivated!
fs2=0.0d0
do i = 1, 3
   do j = 1, Natoms
        fs2 = fs2 + dxi_act(i,j) * dxi_act(i,j) / mass(j)
   end do
end do
coeff1 = 2.0d0 * pi * beta
fs2 = fs2 / coeff1
fs = sqrt(fs2)
log_fs = log(fs)
coeff2 = -1.0d0 / beta

! Add bias term to potential
V_add = V_add + coeff2 * log_fs
! Add bias term to forces
!   if (k_force .ne. 0.d0) then
do i = 1, 3
   do j = 1, Natoms
      dhams = 0.0d0
      do i2 = 1, 3
         do j2 = 1, Natoms
            dhams = dhams + d2xi_act(i2,j2,i,j) * dxi_act(i2,j2) / mass(j2)
         end do
      end do
      dhams = dhams * coeff2 / (coeff1 * fs2)
      do k=1,nbeads
         grad_xyz(i,j,k) = grad_xyz(i,j,k) + dhams
      end do
   end do
end do

return
end subroutine umbrella

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: chrinide/EVB-QMDFF

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: chrinide/EVB-QMDFF

max_forks_repo_licenses: MIT Unlicense


text:
      SUBROUTINE PHO_CASCMOD
C*********************************************************************

C     Initializes default set of particles for cascade calcualtions

C*********************************************************************
      IMPLICIT NONE
      INTEGER i

      INCLUDE 'inc/pobeam'
      INCLUDE 'inc/poinou'
      INCLUDE 'inc/podebg'

      INTEGER projids(4)

      DATA projids/2212 , 2112 , 211 , 321/
C      3122,3112,130,111 /

      DO i = 1 , 4
         CALL PHO_SETPAR(1,projids(i),0,0.D0)
         CALL PHO_SETPAR(2,2212,0,0.D0)
         CALL PHO_SETPCOMB
         CALL PHO_SETPAR(2,2112,0,0.D0)
         CALL PHO_SETPCOMB
      END DO

      END SUBROUTINE

max_stars_repo_name: pzhristov/DPMJET

max_issues_repo_name: pzhristov/DPMJET

max_forks_repo_name: pzhristov/DPMJET


text:
!!# MODULE <<PAR_QuasiDiffusion>>
MODULE PAR_QuasiDiffusion

!!## PURPOSE
!! Provide the runtime parameters needed for QuasiDiffusion.

!!## DEFAULT IMPLICIT
IMPLICIT NONE

!!## EXTERNAL PARAMETERS

!! * system formulation options
!
INTEGER,PARAMETER :: QDF_System01 = 1
INTEGER,PARAMETER :: QDF_System02 = 2
INTEGER,PARAMETER :: QDF_System03 = 3
INTEGER,PARAMETER :: QDF_System04 = 4
INTEGER,PARAMETER :: QDF_System05 = 5
!
CHARACTER(*),PARAMETEr :: QDF_KEY_System(1:5) = (/"GGK0",& !Goldin,Goldina,Kolpakov on essential meshes only
                                                  "GGK ",& !Goldin,Goldina,Kolpakov on cell based meshes
                                                  "AK  ",& !Aristova,Kolpakov on cell based
                                                  "JM  ",& !Jim Morel on cell based
                                                  "JMAC"/) !Jim Morel on cell-based PLUS acceleration only capability
!! * system solver options, <Solver>
!
INTEGER,PARAMETER :: QDF_Direct    = 01
INTEGER,PARAMETER :: QDF_Iterative = 02
!
CHARACTER(*),PARAMETER :: QDF_KEY_Solver(1:2)=(/"Direct   ",&
                                                "Iterative"/)
!
!
!! * system subsolver options, <SubSolver>
!
INTEGER,PARAMETER :: QDF_SMLIB_GMRES       = 01
INTEGER,PARAMETER :: QDF_SMLIB_BiCGStab    = 02
INTEGER,PARAMETER :: QDF_SMLIB_CGS         = 03
INTEGER,PARAMETER :: QDF_SMLIB_BandedGauss = 04
INTEGER,PARAMETER :: QDF_LAPACK_DenseGauss = 05
INTEGER,PARAMETER :: QDF_SPARSKIT_BiCGstab = 06
INTEGER,PARAMETER :: QDF_Botchev_BiCGstab2 = 07
INTEGER,PARAMETER :: QDF_nspiv             = 08
INTEGER,PARAMETER :: QDF_SPARSKIT_CG       = 09
INTEGER,PARAMETER :: QDF_SPARSKIT_CGNR     = 10
INTEGER,PARAMETER :: QDF_SPARSKIT_BCG      = 11
INTEGER,PARAMETER :: QDF_SPARSKIT_DBCG     = 12
INTEGER,PARAMETER :: QDF_SPARSKIT_TFQMR    = 13
INTEGER,PARAMETER :: QDF_SPARSKIT_FOM      = 14
INTEGER,PARAMETER :: QDF_SPARSKIT_GMRES    = 15
INTEGER,PARAMETER :: QDF_SPARSKIT_FGMRES   = 16
INTEGER,PARAMETER :: QDF_SPARSKIT_DQGMRES  = 17
!!
CHARACTER(*),PARAMETER :: QDF_KEY_SubSolver(1:17)=(/"SMLIB-GMRES      ",&
                                                    "SMLIB-BiCGStab   ",&
                                                    "SMLIB-CGS        ",&
                                                    "SMLIB-BandedGauss",&
                                                    "LAPACK-DenseGauss",&
                                                    "SPARSKIT-BiCGstab",&
                                                    "Botchev-BiCGstab2",&
                                                    "nspiv            ",&
                                                    "SPARSKIT-CG      ",&
                                                    "SPARSKIT-CGNR    ",&
                                                    "SPARSKIT-BCG     ",&
                                                    "SPARSKIT-DBCG    ",&
                                                    "SPARSKIT-TFQMR   ",&
                                                    "SPARSKIT-FOM     ",&
                                                    "SPARSKIT-GMRES   ",&
                                                    "SPARSKIT-FGMRES  ",&
                                                    "SPARSKIT-DQGMRES "/)
!
CHARACTER(*),PARAMETER :: QDF_KEY_ReorderStrategies(1:5)=(/"standard   ",&
                                                           "RCM        ",&
                                                           "RevRCM     ",&
                                                           "MinDegree  ",&
                                                           "MinNeighbor"/)
!
INTEGER,PARAMETER :: QDF_PRECONDITIONER_ILU = 1
INTEGER,PARAMETER :: QDF_PRECONDITIONER_ApproxInv = 2
!
CHARACTER(*),PARAMETER :: QDF_KEY_Preconditioners(1:2)=(/"ILU      ",&
                                                         "ApproxInv"/)
!
!
!! * the options for evaluating areas of half-cells <HalfCellShape>
!
INTEGER,PARAMETER :: QDF_FaceBisect = 1
INTEGER,PARAMETER :: QDF_Corners    = 2
CHARACTER(*),PARAMETER :: QDF_KEY_HalfCellShape(1:2)=(/"FaceBisect",&
                                                       "Corners   "/)
!
!
!! * the options for evaluating interfaces to use:
!!   S is Strong, C is Current, W is Weak, F is Flux and
!!   EF is Eddington and Flux combined <InterfaceType>
!
INTEGER,PARAMETER :: QDF_SC_PLUS_WF  = 1
INTEGER,PARAMETER :: QDF_SF_PLUS_WC  = 2
INTEGER,PARAMETER :: QDF_SC_PLUS_WEF = 3
INTEGER,PARAMETER :: QDF_SEF_PLUS_WC = 4
!
CHARACTER(*),PARAMETER :: QDF_KEY_InterfaceType(1:4)=(/"SC+WF ",&
                                                       "SF+WC ",&
                                                       "SC+WEF",&
                                                       "SEF+WC"/)
!
!
!! * the options for evaluating boundary conditions
INTEGER,PARAMETER :: QDF_Goldin        = 1
INTEGER,PARAMETER :: QDF_Larsen        = 2
INTEGER,PARAMETER :: QDF_reflective    = 3
INTEGER,PARAMETER :: QDF_fixedPhi      = 4
INTEGER,PARAMETER :: QDF_fixedJ        = 5
INTEGER,PARAMETER :: QDF_functionPhi   = 6
INTEGER,PARAMETER :: QDF_functionJ     = 7
INTEGER,PARAMETER :: QDF_functionMixed = 8
!
CHARACTER(*),PARAMETER :: QDF_KEY_LOBC(1:8)=(/"Goldin            ",&
                                              "Larsen            ",&
                                              "Reflective        ",&
                                              "fixedPhi          ",&
                                              "fixedJ            ",&
                                              "functionPhi       ",&
                                              "functionJ         ",&
                                              "functionMixed     "/)
!
!! * scalar flux function in a cell generation options
INTEGER,PARAMETER :: QDF_FLAT         = 1
INTEGER,PARAMETER :: QDF_LINEAR_LLS   = 2
INTEGER,PARAMETER :: QDF_LINEAR_GAUSS = 3
CHARACTER(5),PARAMETER :: QDF_KEY_CellFunction(1:3)=(/"Flat ",&
                                                      "LLS  ",&
                                                      "Gauss"/)
!! * scalar flux function in a cell fixup options
INTEGER,PARAMETER :: QDF_NONE = 1
INTEGER,PARAMETER :: QDF_TAU  = 2
INTEGER,PARAMETER :: QDF_OI3  = 3
INTEGER,PARAMETER :: QDF_Min0 = 4
INTEGER,PARAMETER :: QDF_Switch = 5
INTEGER,PARAMETER :: QDF_TAUN  = 6
INTEGER,PARAMETER :: QDF_TAUDN  = 7
CHARACTER(6),PARAMETER :: QDF_KEY_NonlinearFixup(1:7)=(/"None  ",&
                                                        "Tau   ",&
                                                        "OI3   ",&
                                                        "Min0  ",&
                                                        "Switch",&
                                                        "TauN  ",&
                                                        "TauDN "/)

END MODULE

max_stars_repo_name: wawiesel/tapack

max_issues_repo_name: wawiesel/tapack

max_forks_repo_name: wawiesel/tapack


text:
Some of this discussion was moved from Identity/Talk

I put this in as a feature request on the new localwiki, but now thinking it could provide some solution to this issue as well.  Ive had some very useful offline conversations with a few editors during some of the more recent wiki blowups, (these were with people whove posted contact information on their profiles), and it was nice because it allowed for both direct personal interactions, and privacy.  While posting contact info is probably much less viable than getting everyone to use their realname,  I think a private message system on daviswiki could be helpful.

Private messaging could allow people to remain anonymous to the overall wiki, but engage in private conversation with other editors to discuss personal issues, identities, privacy concerns, etc; if so inclined.  PMing would provide a more realistic farmers market scenario, i.e. you wouldnt have to share your name with everyone you meet, but if you thought someone was cool and wanted to get to know them better, you could do so via private message.  Similarly, if a skimask wearing editor is upset about their service at the popsicle stand, you could send them a private message and try and find out if they are just being trolls or if their face was actually, legitimately cold.   Additionally, if someone really wanted to retain their privacy, but wanted to provide some legitimacy to a more critical edit, they could private message a editor they feel is trustworthy and see if they would vouch for them (this is done quite frequently on reddit).  In fact, I think Tom did something similar with the Taxi talk page a couple weeks ago.  A driver was accused of having DUIs and Art privately emailed driving/insurance records to Tom, who then backed up Art.  Art was able to protect his drivers privacy, and the accusations were knocked down a good bit when Tom backed him up.

This would of course mean more work for the localwiki team, (sorry, not trying to be unappreciative, guys.  I think youve already done a metric shitton of great work), but just throwing it out there.  Theres probably some other concerns too (fostering a behindthescenes paranoia, innercircle stuff, etc), but might be worth thinking about.  Users/JeffTolentino
  Some of us have been asking for this for some time... it would also help solve the problem of it being hard to get the attention of new editors in order to communicate with them.   Users/CovertProfessor
  One thing I like about DavisWiki is that there is no PM. The was valuable to me because I am socially inept, and only from this do I learn how real human beings interact with other real human beings to resolve conflicts. It is only through this that I establish trust that the editors are acting in good faith, and have high moral integrity. Without it, these values are lost for me. I could not speak for the experience of others, so this is just one sample point. The loss (I would face) might be justifiable for the potential benefit of being more inclusive. If there is a PM feature, may I ask how this would be resolved: Suppose A is an anonymous user, who wrote something, that B, a real name user found objective. In the old days, B would post publicly (either on the page being edited or on As user page). In the new days, B would have to send a PM to A, because by default, B cannot assume that A would not be offended. The question is this: When A receives that PM from B, could A still opt to discuss the issue in public? According to our logic, the change in the policy is to protect A. But if A chooses to waive that protection, would A opt to discuss in public just like the old days? The second question is this: Suppose it is the policy that A could opt to discuss in public, when I receive a PM, how do I tell if I am A or B? How do I know that I have the right to opt to speak in public? Suppose I ask B, but B said that he is actually A and I have to protect his privacy. At that point, what am I suppose to do if I feel that I am a victim of oppression, targeted by B? About half a year ago, I was on another forum where the moderators asked me to continue a discussion through PM, I agreed, with the condition that they need to publicly acknowledge the questions I asked them. They never made that public acknowledgement and in their PM replies they never cited any questions I made. Then effectively, if you look at history of their replies, you could claim that the list of questions were never made, thus by doing so, they have both censored information and falsified history. The chilling problem is that the issue that was discussed was not enough of great importance that would affect peoples lives, and nothing was at stake other than saving faces. However, if the tool for censorship could corrupt people at such small scale to save face, how do you ensure that people are not would keep their integrity when their lives are at stake? I am not against PM, but I think power and options could corrupt people. I want to know whether A could still opt to discuss in public and how that could work in the policy. If power could corrupt good people, before you attain that power, you need to setup traps to catch yourself in case you do get corrupted, while you are still uncorrupted. If you have seen good people becoming corrupted, and you are not arrogant enough to believe that you are a saint, then this is the kind of things you worry about. Users/EdgarWai

  EW:  I appreciate your response.  My thoughts below are in italics  Users/JeffTolentino
   Suppose A is an anonymous user, who wrote something, that B, a real name user found objective. In the old days, B would post publicly (either on the page being edited or on As user page). In the new days, B would have to send a PM to A, because by default, B cannot assume that A would not be offended. The question is this: When A receives that PM from B, could A still opt to discuss the issue in public?
    Sure.  The wiki would still be an open forum.  B wouldnt even have to send a PM to A in the first place.  PMing would just provide the option to reach out privately if B or A wanted to.
   According to our logic, the change in the policy is to protect A. But if A chooses to waive that protection, would A opt to discuss in public just like the old days?
    Yes.  It would be up to A and B to respond whatever way they felt was appropriate.
   Suppose I ask B, but B said that he is actually A and I have to protect his privacy. At that point, what am I suppose to do if I feel that I am a victim of oppression, targeted by B?
    Im actually a little confused about what your are asking here.  Are you saying if you find out about a separate anonymous account that a realname user uses, do you need to respect their privacy?  And what happens if they start harassing you on either account?  If thats what youre asking, then I guess the answer would still be up to you.  If they werent harassing you, it would look pretty lame on your part to reveal a secret made in confidence to the rest of the wiki for no apparent reason.  If they were harassing you, the fact that you knew about their alternate account would allow you some advantage if things get ugly.  If you do reveal it, I suppose they could still deny it, but in the long run, I dont think this scenario has much advantage over the current situation.  People already do talk offline here when contact information is available, and people have been both respectful with offline info and spilled the beans on each other.  The PM option would just allow people an alternative to posting private contact info online.
   About half a year ago, I was on another forum where the moderators asked me to continue a discussion through PM, I agreed, with the condition that they need to publicly acknowledge the questions I asked them. They never made that public acknowledgement and in their PM replies they never cited any questions I made. Then effectively, if you look at history of their replies, you could claim that the list of questions were never made, thus by doing so, they have both censored information and falsified history. The chilling problem is that the issue that was discussed was not enough of great importance that would affect peoples lives, and nothing was at stake other than saving faces. However, if the tool for censorship could corrupt people at such small scale to save face, how do you ensure that people are not would keep their integrity when their lives are at stake?
    If this happened on DW, I suppose you could still go back to the public forum and say I tried to resolve this by PM, but such and such happened and none of the problems Im concerned about were dealt with adequately.  Sure, there is no public record of the PM conversation, but by bringing your concerns back into the public forum, whoever you were having a problem with in the first place would be publicly confronted with not solving the problem privately and the interaction would continue from there.  It actually sounds like the key difference with that situation, however, was that the forum you were on was moderated, whereas DW is not moderated, (I guess JW and PN have some power to ban spammers, but I think both of them have pretty good records regarding that privilege).  Nobody else can censor any of the public statements on DW.  Things do get removed, but there is always a record in the edit history, and you can always revert or redit something if you still have a problem with something.  Really, Id want everyone to be accountable in public for what they say all the time, but my intent with the PM proposal is to address the privacy issues that AFB and CP are bringing up here, i.e. how can we facilitate accountability and privacy at the same time.
   I want to know whether A could still opt to discuss in public and how that could work in the policy.
    Yeah, anyone could still discuss whatever they want in public.  You couldnt force somebody else to respond in public, but you cant really do that now anyway.
   If power could corrupt good people, before you attain that power, you need to setup traps to catch yourself in case you do get corrupted, while you are still uncorrupted. If you have seen good people becoming corrupted, and you are not arrogant enough to believe that you are a saint, then this is the kind of things you worry about.
    I dont think having the ability to PM gives anyone a whole lot of power.  You can still air you grievances publicly if you want.  PMing would just let people who are worried about privacy, a chance to provide proof of accountability to editors they trust.  Its not going to solve every problem, but it could be a useful tool for some people.
  JT: Thank you for responding. One of the differences is between the forum I mentioned and DavisWiki, is that that other forum has no public edit record. If power could corrupt, the people would already have been corrupted. Therefore when I agree that the concerns above are resolved (in the way I see them). I dont have more concerns from my side against PM. Thank you for addressing it. Users/EdgarWai
   EdgarWai, again, I find everything you say to be amazing.  I don’t have answers yet but your comments are inspiring to those ends.  JeffTolentino, I like your suggestions, too.  Accountability is an important issue, and I think Davisites have the capacity to find an innovative way to achieve this while getting more people involved. Users/ActionFigureBarbie
  A private messaging system is a bit dangerous. The community is a whole lot better if everyone can be involved in a discussion. By offering a private messaging system, people will actually use it—and if it is better than the comment bar (very easy to do), they will use it by default. Yes, there is sometimes a need for discussion of things in private, but its not the common case and email will do. When to use private messages and when to use userpages... one more thing to teach newcomers. Oy. Of course, I also hate the WYSIWYG editor in XsaplingX localwiki, so feel free to dismiss me as a stuckinhisways oldtimer. And I place a higher value on openness and transparency than most. Users/WilliamLewis
   I just want to say that I have similar thoughts about whether PM would make DavisWiki really quiet. But since I dont feel strongly whether PM would be for better or for worse, I will see whether that problem is addressed according to William. Users/EdgarWai
  from a technical standpoint a PM system is also something that can be abused by spammers / people up to no good Daubert a wise woman once taught me that less is more...
   I actually get tons on spam on my personal email account, and Ive never gotten any PM spam through FB or Reddit.  Part of the reason I dont post my personal email address here, is because I dont want more spam!  A big part of me sympathizes with you less is more point, however.  In a perfect world, Id still advocate simply using your real name as the best solution, but AFB and CP have a number of good points and wanted to find a way to meet those too.Users/JeffTolentino


text:
subroutine BladeIndVel(NT,ntTerm,NBE,NB,NE,XP,YP,ZP,UP,VP,WP,DUDX,Mode,CalcDer)

    use blade
    use wake

    integer :: nt, ntTerm, nbe, nb, ne, Mode, CalcDer
    real :: XP, YP, ZP, UP, VP, WP, DUDX
    integer :: i, j, k, nei, nej, nt1
    integer :: VFlag

    ! COMPUTE THE VORTEX INDUCED VELOCITY AT POINT XP,YP,ZP FROM THE BLADE SYSTEM
    ! Mode: 0 -> Calc whole system
    !       1 -> Calc bound vorticity only
    !       2 -> Calc whole system minus bound vorticity
    ! CalcDer: 1 to calc DUDX

    NT1=NT-1
    UP=0.0
    VP=0.0
    WP=0.0
    DUDX=0.0

    if (Mode == 0) then
        ! Calc induced velocity from the whole blade system (bound, spanwise, and trailing wake vorticity)

        ! CALCULATE THE VELOCITY CONTRIBUTIONS DUE TO TRAILING VORTICIES  ( GT(1:NT-1,:) )
        ! ntTerm represents the furthest away wake elements that are to be considered. (Calculated using user input xstop)
        VFlag=1
        do i=1,ne
            do j=ntTerm,NT1
                Call VorIVel(VFlag,CalcDer,GT(j,i),X(j,i),Y(j,i),Z(j,i),X(j+1,i),Y(j+1,i),Z(j+1,i),XP,YP,ZP,UP,VP,WP,DUDX)
            end do
        end do

        ! CALCULATE THE VELOCITY CONTRIBUTIONS DUE TO SPANWISE VORTICIES, including current bound vorticity ( GS(1:NT,:) )
        do i=1,nb
            nei=(i-1)*(nbe+1)
            do j=1,nbe
                nej=nei+j
                do k=ntTerm,NT
                    if (k==NT) then
                        VFlag=0
                    else
                        VFlag=2
                    end if
                    Call VorIVel(VFlag,CalcDer,GS(k,nej),X(k,nej),Y(k,nej),Z(k,nej),X(k,nej+1),Y(k,nej+1),Z(k,nej+1),XP,YP,ZP,UP,VP,WP,DUDX)
                end do
            end do
        end do

    else if (Mode == 1) then
        ! Calc induced velocity from bound vorticity only

        VFlag=0
        do i=1,nb
            nei=(i-1)*(nbe+1)
            do j=1,nbe
                nej=nei+j
                k=NT
                Call VorIVel(VFlag,CalcDer,GS(k,nej),X(k,nej),Y(k,nej),Z(k,nej),X(k,nej+1),Y(k,nej+1),Z(k,nej+1),XP,YP,ZP,UP,VP,WP,DUDX)
            end do
        end do

    else
        ! Calc induced velocity from spanwise and trailing wake vorticity (exclude bound vorticity)

        ! CALCULATE THE VELOCITY CONTRIBUTIONS DUE TO TRAILING VORTICIES  ( GT(1:NT-1,:) )
        ! ntTerm represents the furthest away wake elements that are to be considered. (Calculated using user input xstop)
        VFlag=1
        do i=1,ne
            do j=ntTerm,NT1
                Call VorIVel(VFlag,CalcDer,GT(j,i),X(j,i),Y(j,i),Z(j,i),X(j+1,i),Y(j+1,i),Z(j+1,i),XP,YP,ZP,UP,VP,WP,DUDX)
            end do
        end do

        ! CALCULATE THE VELOCITY CONTRIBUTIONS DUE TO SPANWISE VORTICIES, NOT including current bound vorticity ( GS(1:NT-1,:) )
        VFlag=2
        do i=1,nb
            nei=(i-1)*(nbe+1)
            do j=1,nbe
                nej=nei+j
                do k=ntTerm,NT1
                    Call VorIVel(VFlag,CalcDer,GS(k,nej),X(k,nej),Y(k,nej),Z(k,nej),X(k,nej+1),Y(k,nej+1),Z(k,nej+1),XP,YP,ZP,UP,VP,WP,DUDX)
                end do
            end do
        end do

    end if

    return
end subroutine BladeIndVel

max_stars_repo_name: FloatVAWT/CACTUS

max_issues_repo_name: FloatVAWT/CACTUS

max_forks_repo_name: FloatVAWT/CACTUS


text:
Chill Tree is the creation of local artist Emily Earlenbaugh; who believes that art feeds the soul, and the people are hungry.

Emily has been making personalized art and jewelry since she was a child, and loves to empower people by showing them reflections of themselves in art.

The mission at Chill Tree is to create affordable art and magical items for your life’s journey. This is why we started offering Art by Donations: At Chill Tree we see artmaking as a community service. We want everyone to have access to the experience of personalized art, so instead of listing prices, we list the cost of the materials and the time spent. You decide the donation.

Check out www.chilltree.com for more info


text:
      SUBROUTINE DIS2D8
C
C     2-D, 8 GRID POINT ISOPARAMETRIC STRUCTURAL ELEMENT
C     DIFFERENTIAL STIFFNESS MATRIX ROUTINE
C
      REAL             KX,KY
      DOUBLE PRECISION KIJ,G,B,XI,ETA,DNXI,DNETA,XX,TB,DNL,PT,H,XJB,
     1                 XXJB,DETERM,DNC,GSUBE,DUMARG,
     2                 BT,TEMPAR,TEMP,DNX,DNY,SAVE,TSAVE,
     3                 KWD(36),CID(18),CJD(18),KMULT(18),DHH,
     4                 THICK,PSTMUL(9),PREMUL(9),E1T
      DIMENSION        SIG(3),ALPHAS(3),ST(3),SEMP(9),TTB(9),STB(9),
     1                 BB(72),DB(72),S(6),R(9),SE1T(6),DN(8),
     2                 G(9),QQ(15),XI(8),ETA(8),TB(9),XY1(3),XY2(3),
     3                 B(12),BT(12),ECPT(1),TEMP(9),TEMPAR(1),DNX(1),
     4                 DNY(1),DNXI(1),DNETA(1),SAVE(72),TSAVE(72),
     5                 VEC(3),VVEC(3),VECI(3),VECJ(3),VECK(3),E1T(9),
     6                 IWS(2,3)
      COMMON /DS1AAA/  NPVT,ICSTM,NCSTM
      COMMON /DS1AET/  NECPT(1),NGRID(8),ID1,TH,MATID1,T,ISYS1,X1,Y1,Z1,
     1                 ISYS2,X2,Y2,Z2,ISYS3,X3,Y3,Z3,ISYS4,X4,Y4,Z4,
     2                 ISYS5,X5,Y5,Z5,ISYS6,X6,Y6,Z6,ISYS7,X7,Y7,Z7,
     3                 ISYS8,X8,Y8,Z8,TTEMP,EDT,ISETNO,TGRID(8),DISP(24)
      COMMON /MATIN /  MATID,INFLAG,ELTEMP,STRESS,SINTH,COSTH
      COMMON /MATOUT/  G11,G12,G13,G22,G23,G33,RHO,ALPHA1,ALPHA2,ALP12,
     1                 TREF,GE,KX,KY,C
      COMMON /DS1ADP/  KIJ(36),XX(16),DNC(16),DNL(16),XXJB(2,2),XJB(4),
     1                 PT(3),H(3),G,B,BT,TB,DETERM,GSUBE,DUMARG,TSAVE
      EQUIVALENCE      (ALPHAS(1),ALPHA1),
     1                 (ECPT(1),NECPT(1)),(TEMP(1),B(1)),
     2                 (DNC(1),DNXI(1)),(DNC(9),DNETA(1)),
     3                 (DNL(1),DNX(1)),(DNL(9),DNY(1)),(QQ(1),G11),
     4                 (TEMPAR(1),BT(1)),(XY1(1),X1),(XY2(1),X2)
      DATA     XI   /  -1.D0, 1.D0, 1.D0,-1.D0, 0.D0, 1.D0, 0.D0,-1.D0/
      DATA     ETA  /  -1.D0,-1.D0, 1.D0, 1.D0,-1.D0, 0.D0, 1.D0, 0.D0/
C
C     ECPT LIST
C                                      IN
C                                      THIS
C     ECPT       DESCRIPTION           ROUTINE        TYPE
C     ******************************************************************
C     ECPT( 1) = ELEMENT ID            NECPT(1)       INTEGER
C     ECPT( 2) = GRID POINT 1          NGRID(1)       INTEGER
C     ECPT( 3) = GRID POINT 2          NGRID(2)       INTEGER
C     ECPT( 4) = GRID POINT 3          NGRID(3)       INTEGER
C     ECPT( 5) = GRID POINT 4          NGRID(4)       INTEGER
C     ECPT( 6) = GRID POINT 5          NGRID(5)       INTEGER
C     ECPT( 7) = GRID POINT 6          NGRID(6)       INTEGER
C     ECPT( 8) = GRID POINT 7          NGRID(7)       INTEGER
C     ECPT( 9) = GRID POINT 8          NGRID(8)       INTEGER
C     ECPT(10) = COORD SYS ID-STRESS   ID1            INTEGER
C     ECPT(11) = ANIS. MATERIAL ANGLE  TH             REAL
C     ECPT(12) = MATERIAL ID           MATID1         INTEGER
C     ECPT(13) = THICKNESS             T              REAL
C     ECPT(14) = COORD SYS ID 1        ISYS1          INTEGER
C     ECPT(15) = X1                    X1             REAL
C     ECPT(16) = Y1                    Y1             REAL
C     ECPT(17) = Z1                    Z1             REAL
C     ECPT(18) = COORD SYS ID 2        ISYS2          INTEGER
C     ECPT(19) = X2                    X2             REAL
C     ECPT(20) = Y2                    Y2             REAL
C     ECPT(21) = Z2                    Z2             REAL
C     ECPT(22) = COORD SYS ID 3        ISYS3          INTEGER
C     ECPT(23) = X3                    X3             REAL
C     ECPT(24) = Y3                    Y3             REAL
C     ECPT(25) = Z3                    Z3             REAL
C     ECPT(26) = COORD SYS ID 4        ISYS4          INTEGER
C     ECPT(27) = X4                    X4             REAL
C     ECPT(28) = Y4                    Y4             REAL
C     ECPT(29) = Z4                    Z4             REAL
C     ECPT(30) = COORD SYS ID 5        ISYS5          INTEGER
C     ECPT(31) = X5                    X5             REAL
C     ECPT(32) = Y5                    Y5             REAL
C     ECPT(33) = Z5                    Z5             REAL
C     ECPT(34) = COORD SYS ID 6        ISYS6          INTEGER
C     ECPT(35) = X6                    XL             REAL
C     ECPT(36) = Y6                    Y6             REAL
C     ECPT(37) = Z6                    Z6             REAL
C     ECPT(38) = COORD SYS ID 7        ISYS7          INTEGER
C     ECPT(39) = X7                    X7             REAL
C     ECPT(40) = Y7                    Y7             REAL
C     ECPT(41) = Z7                    Z7             REAL
C     ECPT(42) = COORD SYS ID 8        ISYS8          INTEGER
C     ECPT(43) = X8                    X8             REAL
C     ECPT(44) = Y8                    Y8             REAL
C     ECPT(45) = Z8                    Z8             REAL
C     ECPT(46) = ELEMENT TEMP          TTEMP          REAL
C     ECPT(47) = 0.                    EDT            REAL
C     ECPT(48) = TEMPERATURE SET       ISETNO         INTEGER
C     ECPT(49) = *
C     ECPT(. ) = *  GRID POINT TEMPERATURES
C     ECPT(56) = *
C     ECPT(57) = *
C     ECPT(. ) = *  TRANSLATIONAL DOF-S OF GRIDS FOR THIS ELEMENT
C     ECPT(80) = *
C
C
C     TEST FOR PIVOT POINT
C
      DO 10 KK = 1,8
      IF (NGRID(KK) .EQ. NPVT) GO TO 20
   10 CONTINUE
C
C     IF FALL HERE NO ELEMENT GRID POINT IS THE PIVOT POINT
C
      CALL MESAGE (-30,34,ECPT(1))
C
C     UNIT I VECTOR IS FROM GRID POINT 1 TO GRID POINT 2
C
   20 DO 30 I = 1,3
      VECI(I) = XY2(I) - XY1(I)
   30 CONTINUE
      VECIL = SQRT(VECI(1)**2 + VECI(2)**2 + VECI(3)**2)
      IF (VECIL .EQ. 0.0) GO TO 60
      VECI(1) = VECI(1)/VECIL
      VECI(2) = VECI(2)/VECIL
      VECI(3) = VECI(3)/VECIL
C
C     K VECTOR IS OBTAINED BY CROSSING I INTO VECTOR FROM GRID PT. 1 TO
C     GRID
C
      VECK(1) = VECI(2)*(Z4-Z1) - VECI(3)*(Y4-Y1)
      VECK(2) = VECI(3)*(X4-X1) - VECI(1)*(Z4-Z1)
      VECK(3) = VECI(1)*(Y4-Y1) - VECI(2)*(X4-X1)
      VECKL=SQRT(VECK(1)**2 + VECK(2)**2 + VECK(3)**2)
      IF (VECKL .EQ. 0.0) GO TO 60
      VECK(1) = VECK(1)/VECKL
      VECK(2) = VECK(2)/VECKL
      VECK(3) = VECK(3)/VECKL
C
C     J VECTOR IS OBTAINED BY CROSSING K INTO I
C
      VECJ(1) = VECK(2)*VECI(3) - VECK(3)*VECI(2)
      VECJ(2) = VECK(3)*VECI(1) - VECK(1)*VECI(3)
      VECJ(3) = VECK(1)*VECI(2) - VECK(2)*VECI(1)
C
      E1T(1) = VECI(1)
      E1T(2) = VECI(2)
      E1T(3) = VECI(3)
      E1T(4) = VECJ(1)
      E1T(5) = VECJ(2)
      E1T(6) = VECJ(3)
      E1T(7) = VECK(1)
      E1T(8) = VECK(2)
      E1T(9) = VECK(3)
      DO 40 I = 1,6
   40 SE1T(I) = E1T(I)
C
C     STORE ELEMENT COORDS FOR GRIDS 1 AND 2
C
      XX(1) = 0.
      XX(2) = 0.
      XX(3) = VECIL
      XX(4) = 0.
C
C     FOR GRIDS 3-8, THE X COORDINATE IS THE DOT PRODUCT OF HTE VECTOR
C     FROM GRID POINT 1 TO THE GRID POINT AND THE I VECTOR. THE Y COORD.
C     IS THE L OF THE I VECTOR CROSSED INTO THE VECTOR FROM GRID 1 TO
C     THE GRID POINT.
C
      DO 50 I = 3,8
      IXX  = 2*I - 1
      ISUB = 4*I + 11
      VEC(1)  = ECPT(ISUB  ) - X1
      VEC(2)  = ECPT(ISUB+1) - Y1
      VEC(3)  = ECPT(ISUB+2) - Z1
      XX(IXX) = VEC(1)*VECI(1) + VEC(2)*VECI(2) + VEC(3)*VECI(3)
      VVEC(1) = VECI(2)*VEC(3) - VECI(3)*VEC(2)
      VVEC(2) = VECI(3)*VEC(1) - VECI(1)*VEC(3)
      VVEC(3) = VECI(1)*VEC(2) - VECI(2)*VEC(1)
      XX(IXX+1) = SQRT(VVEC(1)**2 + VVEC(2)**2 + VVEC(3)**2)
   50 CONTINUE
      GO TO 70
C
C     INAPPROPRIATE GEOMETRY
C
   60 CALL MESAGE (30,31,ECPT(1))
      NOGO = 1
      RETURN
C
C     COMPUTE MATERIAL PROPERTIES
C
   70 TTH   = TH*3.1415927/180.
      SINTH = SIN(TTH)
      COSTH = COS(TTH)
      ELTEMP= TTEMP
      INFLAG= 2
      MATID = MATID1
      CALL MAT (ECPT(1))
      DO 80 I = 1,3
   80 G(I)  = QQ(I)
      G(4)  = QQ(2)
      G(5)  = QQ(4)
      G(6)  = QQ(5)
      G(7)  = QQ(3)
      G(8)  = QQ(5)
      G(9)  = QQ(6)
      THICK = T
      DO 90 I = 1,9
   90 R(I)  = G(I)
      IF (ISETNO .NE. 0) CALL GMMATS (R,3,3,0,ALPHAS,3,1,0,ST)
C
C     ZERO OUT THE KIJ AND SAVE MATRICES
C
      DO 100 I = 1,36
      KWD(I) = 0.D0
  100 KIJ(I) = 0.D0
      DO 110 I = 1,72
  110 SAVE(I) = 0.D0
C
      PT(1) =-0.57735027D0
      PT(2) =-PT(1)
      H(1)  = 1.D0
      H(2)  = 1.D0
      IF (ID1 .EQ. 2) GO TO 120
      PT(1) =-0.77459667D0
      PT(2) = 0.D0
      PT(3) =-PT(1)
      H(1)  = 5.D0/9.D0
      H(2)  = 8.D0/9.D0
      H(3)  = H(1)
C
C     2 OR 3 QUADRATURE POINTS
C
  120 DO 410 III = 1,ID1
      DO 410 JJJ = 1,ID1
C
C     COMPUTE GAUSS POINT STRESSES
C
C
C     COMPUTE DERIVATIVES WITH RESPECT TO XI AND ETA
C     EACH GRID POINT
C
      DO 130 N = 1,4
      DNXI(N)  = .25D0*XI(N)*(1.D0+PT(JJJ)*ETA(N))*
     1           (2.D0*PT(III)*XI(N)+PT(JJJ)*ETA(N))
      DNETA(N) = .25D0*ETA(N)*(1.D0+PT(III)*XI(N))*
     1           (PT(III)*XI(N)+2.D0*PT(JJJ)*ETA(N))
  130 CONTINUE
      DO 140 N = 5,7,2
C
      DNXI(N)  = -PT(III)*(1.D0+PT(JJJ)*ETA(N))
      DNETA(N) = .5D0*(1.D0-PT(III)*PT(III))*ETA(N)
  140 CONTINUE
C
      DO 150 N = 6,8,2
      DNXI(N)  = .5D0*XI(N)*(1.D0-PT(JJJ)*PT(JJJ))
      DNETA(N) = -PT(JJJ)*(1.D0+PT(III)*XI(N))
  150 CONTINUE
C
C     COMPUTE JACOBEAN
C
C           N1XI   N2XI   N3XI   N4XI   N5XI   N6XI   N7XI   N8XI
C     DNC = N1ETA  N2ETA  N3ETA  N4ETA  N5ETA  N6ETA  N7ETA  N8ETA
C
C          X1  Y1
C          X2  Y2
C          X3  Y3
C     XX = X4  Y4
C          X5  Y5
C          X6  Y6
C          X7  Y7
C          X8  Y8
C
      CALL GMMATD (DNC,2,8,0,XX,8,2,0,XJB)
C
C     XJB IS ROW-STORED-IT MUST BE COLUMN-STORED AND DOUBLY DIMENSIONED
C     FOR INVERSION
C
      K = 0
      DO 160 I = 1,2
      DO 160 J = 1,2
      K = K + 1
  160 XXJB(I,J) = XJB(K)
C
C     COMPUTE INVERSE AND DETERMINANT OF JACOBEAN
C
      CALL INVERD (2,XXJB,2,DUMARG,0,DETERM,ISING,IWS)
      IF (ISING .EQ. 2) CALL MESAGE (-30,143,ECPT(1))
      DHH = DETERM*H(III)*H(JJJ)
C
C     COMPUTE DERIVATIVES WITH RESPECT TO X AND Y
C
      K = 0
      DO 170 I = 1,2
      DO 170 J = 1,2
      K = K + 1
  170 XJB(K) = XXJB(I,J)
      CALL GMMATD (XJB,2,2,0,DNC,2,8,0,DNL)
C
C           N1X N2X N3X N4X N5X N6X N7X N8X
C     DNL = N1Y N2Y N3Y N4Y N5Y N6Y N7Y N8Y
C
      DO 180 I = 1,72
  180 BB(I) = 0.
C
C     SET UP INDICATOR FOR GRID POINT TEMPERATURES
C
      IDTEMP = 0
      DO 190 I = 1,8
      IF (TGRID(I) .NE. 0.) GO TO 200
  190 CONTINUE
      GO TO 210
  200 IDTEMP = 1
C
  210 DO 270 N = 1,8
C
      DO 220 I = 1,9
  220 SEMP(I) = 0.
      DO 230 I = 1,6
  230 S(I) = 0.
      S(1) = DNX(N)
      S(4) = DNY(N)
      S(5) = DNY(N)
      S(6) = DNX(N)
C
C     TRANSFORM TO ELEMENT COORDINATES
C
      IF (NECPT(4*N+10) .EQ. 0) GO TO 240
      CALL TRANSS (NECPT(4*N+10),TTB)
      CALL GMMATS (SE1T,2,3,0,TTB,3,3,0,STB)
      GO TO 260
  240 DO 250 I = 1,6
  250 STB(I) = SE1T(I)
  260 CALL GMMATS (S,3,2,0,STB,2,3,0,SEMP(1))
      N3 = 3*N
      BB(N3- 2) = SEMP(1)
      BB(N3- 1) = SEMP(2)
      BB(N3   ) = SEMP(3)
      BB(N3+22) = SEMP(4)
      BB(N3+23) = SEMP(5)
      BB(N3+24) = SEMP(6)
      BB(N3+46) = SEMP(7)
      BB(N3+47) = SEMP(8)
      BB(N3+48) = SEMP(9)
  270 CONTINUE
C
C     BRING IN G MATRIX
C
      CALL GMMATS (R,3,3,0,BB,3,24,0,DB)
C
C     COMPUTE STRESSES
C
      CALL GMMATS (DB,3,24,0,DISP,24,1,0,SIG)
C
C
C     COMPUTE GAUSS POINT  TEMPERATURES
C
      IF (ISETNO .EQ. 0) GO TO 350
      IF (IDTEMP .EQ. 1) GO TO 280
      RGTEMP = ELTEMP - TREF
      GO TO 330
C
C     ALL TEMPERATURES ARE DEFAULT VALUE
C
  280 DO 290 N = 1,4
      DN(N) = .25*(1.+PT(III)*XI(N))*(1.+PT(JJJ)*ETA(N))
     1        *(PT(III)*XI(N)+PT(JJJ)*ETA(N)-1.)
  290 CONTINUE
      DO 300 N = 5,7,2
      DN(N) = .5*(1.-PT(III)*PT(III))*(1.+PT(JJJ)*ETA(N))
  300 CONTINUE
      DO 310 N = 6,8,2
      DN(N) = .5*(1.+PT(III)*XI(N))*(1.-PT(JJJ)*PT(JJJ))
  310 CONTINUE
      GSTEMP = 0.
      DO 320 N = 1,8
      GSTEMP = GSTEMP + DN(N)*TGRID(N)
  320 CONTINUE
      RGTEMP = GSTEMP - TREF
  330 DO 340 I = 1,3
      SIG(I) = SIG(I) - ST(I)*RGTEMP
  340 CONTINUE
C
  350 CONTINUE
C
C     FORM KWD MATRIX
C
      KWD( 1) = SIG(2)
      KWD( 2) =-SIG(3)
      KWD( 7) =-SIG(3)
      KWD( 8) = SIG(1)
      KWD(15) = SIG(1) + SIG(2)
      KWD(16) =-SIG(3)
      KWD(17) = SIG(3)
      KWD(18) = SIG(1) - SIG(2)
      KWD(21) =-SIG(3)
      KWD(27) = SIG(3)
      KWD(33) = SIG(1) - SIG(2)
C
C     FORM CID FOR I = NPVT
C
      DO 360 I = 1,18
  360 CID( I) = 0.D0
      CID( 3) = DNY(KK)
      CID( 6) =-DNX(KK)
      CID( 7) =-.5*DNY(KK)
      CID( 8) = .5*DNX(KK)
      CID(10) = DNX(KK)
      CID(14) = DNY(KK)
      CID(16) =.5*DNY(KK)
      CID(17) =.5*DNX(KK)
C
      CALL GMMATD (CID,6,3,1,KWD,6,6,0,KMULT)
C
C     LOOP FOR THE 8 6X6 PARTITIONS CORRESPONDING TO THE PRESENT
C     PIVOT POINT
C
      DO 400 N = 1,8
C
      DO 370 I = 1,18
  370 CJD(I) = 0.D0
C
      CJD( 3) = DNY(N)
      CJD( 6) =-DNX(N)
      CJD( 7) =-.5*DNY(N)
      CJD( 8) = .5*DNX(N)
      CJD(10) = DNX(N)
      CJD(14) = DNY(N)
      CJD(16) =.5*DNY(N)
      CJD(17) =.5*DNX(N)
C
      CALL GMMATD (KMULT,3,6,0,CJD,6,3,0,TEMPAR(1))
C
C     THROW IN JACOBEAN DETERMINANT AND WEIGHT FACTORS
C
      DO 380 I = 1,9
      TEMPAR(I) = TEMPAR(I)*DHH
  380 CONTINUE
C
C     ADD THE RESULTS OF THIS INTEGRATION TO THE PREVIOUS RESULTS
C
      LL = 9*(N-1)
      DO 390 I = 1,9
      L = LL + I
      SAVE(L) = SAVE(L) + TEMPAR(I)
  390 CONTINUE
C
C     LOOP FOR MORE PARTITIONS
C
  400 CONTINUE
C
C     LOOP FOR MORE GAUSS POINTS
C
  410 CONTINUE
C
C     CHECK ON NECESSITY OF PRE-MULTIPLYING COORDINATE TRANSFORMATIONS
C
      ISUB = 4*KK + 10
      IF (NECPT(ISUB) .EQ. 0) GO TO 420
C
C     ELEMENT TO GLOBAL
C
      CALL TRANSD (NECPT(ISUB),TB)
      CALL GMMATD (E1T,3,3,0,TB,3,3,0,PREMUL)
      GO TO 440
  420 DO 430 I = 1,9
      PREMUL(I) = E1T(I)
  430 CONTINUE
  440 DO 460 N = 1,8
      LL = 9*N - 8
      CALL GMMATD (PREMUL,3,3,1,SAVE(LL),3,3,0,TEMP)
C
C     STORE THE 3 X 3 IN TSAVE
C
      DO 450 I = 1,9
      L = 9*N + I - 9
  450 TSAVE(L) = TEMP(I)
C
  460 CONTINUE
C
C     NOW CHECK ON THE NECESSITY FOR POST-MULTIPLYING TRANSFORMATIONS
C
      DO 500 N = 1,8
      ISUB = 4*N + 10
      LL = 9*N - 8
      IF (NECPT(ISUB) .EQ. 0) GO TO 470
C
C     GLOBAL TO ELEMENT
C
      CALL TRANSD (NECPT(ISUB),TB)
      CALL GMMATD (E1T,3,3,0,TB,3,3,0,PSTMUL)
      GO TO 490
  470 DO 480 I = 1,9
      PSTMUL(I) = E1T(I)
  480 CONTINUE
C
C     POST-MULTIPLY
C
  490 CALL GMMATD (TSAVE(LL),3,3,0,PSTMUL,3,3,0,TEMP)
C
C     FILL OUT THE 6 X 6 PARTITION
C
      KIJ( 1) = TEMP(1)*THICK
      KIJ( 2) = TEMP(2)*THICK
      KIJ( 3) = TEMP(3)*THICK
      KIJ( 7) = TEMP(4)*THICK
      KIJ( 8) = TEMP(5)*THICK
      KIJ( 9) = TEMP(6)*THICK
      KIJ(13) = TEMP(7)*THICK
      KIJ(14) = TEMP(8)*THICK
      KIJ(15) = TEMP(9)*THICK
C
C     INSERT INTO THE OVERALL STIFFNESS MATRIX
C
      CALL DS1B (KIJ,NECPT(N+1))
C
C     LOOP FOR MORE PARTITIONS
C
  500 CONTINUE
C
      RETURN
      END


text:
Address(Evenstar Lane) is one of the many Navigating Davis streets in Davis Village Homes neighborhood that are named after places or characters in J.R.R. Tolkiens books, The Hobbit and The Lord of the Rings.

Elendil Lane branches off from Evenstar Lane.

Evenstar is the last name of Arwen Evenstar, a brave Elven warriormaiden who helped the four hobbits escape from the Nazgul by calling on the river. Arwen was Elronds daughter, and lived in Rivendell and Lorien for roughly 3000 years in contentment, until she met and fell in love with Aragorn. After the War of the Ring, Arwen and Aragorn were married on midyears day, and Arwen became Queen of Gondor.
  From http://www.villagehomesdavis.org/public/about/street_names Village Homes Street Names


text:
	SUBROUTINE HSCINT  ( iret )
C************************************************************************
C* HSCINT - XW								*
C*									*
C* This subroutine initializes default set of colors on a device.	*
C* Each device has its own default colors.  The colors may be set	*
C* by calling the device subroutine DSNAME to set each color by		*
C* name.								*
C*									*
C* HSCINT  ( IRET )							*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C**									*
C* Log:									*
C* J. Whistler/SSAI	 8/91	XW device driver			*
C* M. desJardins/NMC	01/92	Add initialization of background	*
C* K. Brill/NMC		01/92	Add 16 more colors			*
C* K. Brill/NMC		02/92	Eliminate duplicate colors		*
C* K. Brill/NMC		03/92	FIREBRICK -> FIREBRIC			*
C* S. Jacobs/EAI	10/93	Changed color list			*
C* S. Jacobs/NMC	 7/94	General clean up			*
C* S. Maxwell/GSC        6/97   Documentation changes                   *
C************************************************************************
C
C*	Set the colors from the default color table.
C
	CALL XSCINT ( iret )
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
    module StringTests
    use StringUtils
    implicit none
    contains

    function RunStringTests() result (fails)
    integer fails
    character(LEN=:), allocatable :: S, S2
    fails = 0

    if (UpperCase('Low@er')/='LOW@ER' .or. LowerCase('1L oW@eR')/='1l ow@er') then
        print *, 'Error in uppercase/lowercase'
        fails = fails+1
    end if

    S = FormatString('Test %d for %f %s', 91, 3.04, 'result')
    if (S/='Test 91 for 3.0400 result') then
        fails = fails+1
        print *, 'Error in FormatString: ', S
    end if

    S2 = Join(',', S,'XP')
    if (S2 /= S//',XP') then
        fails = fails+1
        print *, 'Error in Join: ', S2
    end if

    if (fails==0) print *, 'OK string tests'

    end function

    end module

max_stars_repo_name: LBJ-Wade/forutils

max_issues_repo_name: LBJ-Wade/forutils

max_forks_repo_name: LBJ-Wade/forutils


text:
	FUNCTION TRAPEZ(Y,N,DX)
c	FINDS AREA BELOW CURVE (X,Y) ACCORDING TO TRAPEZOIDAL RULE.
c	ASSUMES EQUAL SPACING DX BETWEEN N POINTS.
	REAL Y(N)
	TRAPEZ=0.
	IF (N.LE.1) RETURN
	TRAPEZ=(Y(1) + Y(N))*DX/2.
	IF (N.EQ.2) RETURN
	A=0.
	DO I=2,N-1
	A = A + Y(I)
	ENDDO
	TRAPEZ = TRAPEZ + A*DX
	RETURN
	END

max_stars_repo_name: jchavesmontero/galaxpy

max_issues_repo_name: jchavesmontero/galaxpy

max_forks_repo_name: jchavesmontero/galaxpy


text:
! mpiexec -np 4 bin/Debug/WGFM_test2.exe
program mainWGF
    use meshread_mod
    use spherePEC_mod
    implicit none

    real(8), parameter :: pi = 4d0 * atan(1d0)
    complex(8), parameter :: IU = (0d0, 1d0)

    ! Variables
    type(Mesh) :: msh_evl
    complex(8), allocatable :: Esrc(:, :)  ! Source electric field
    complex(8), allocatable :: Eexact(:, :)  ! Exact electric field solution

    real(8) :: ExNormSrc, ExNormTot
    integer :: nEvl

    ! Parameters
    real(8), parameter :: ep0 = 1d0 / (35950207149.4727056d0 * PI)  ! vacuum permittivity [F/m]
    real(8), parameter :: mu0 = PI * 4d-7                           ! vacuum permeability [H/m]
    real(8), parameter :: k0 = 2*PI          ! free-space wavenumber [1/m]

    real(8), parameter :: ep1_r = 1d0          ! upper half-space relative permittivity
    real(8), parameter :: ep1 = ep0*ep1_r     ! upper half-space permittivity [F/m]
    real(8), parameter :: mu1 = mu0           ! upper half-space permeability [H/m]

    real(8), parameter :: w = k0 / sqrt(ep0 * mu0)   ! angular frequency [rad/s]
    real(8), parameter :: k1 = w * sqrt(ep1 * mu1)  ! upper half-space wavenumber [1/m]

    character(len=100), parameter :: file_evl_msh = 'meshes/eval_msh/dome.msh'
    real(8), parameter :: sph_radius = 1d0      ! hemisphere radius

    ! Planewave parameters
    integer, parameter :: n_terms = 100     ! number of terms in Mie series
    real(8), parameter :: alpha = PI/32d0      ! planewave angle of incidence in (0, pi)
    logical, parameter :: TE_mode = .false.  ! planewave mode: TEmode = .true.
                                            !                 TMmode = .false.

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

    ! Load eval mesh
    call load_gmsh(file_evl_msh, msh_evl)
    nEvl = size(msh_evl%POS, 1)

    ! Compute source field Esrc
    call halfsphere_on_plane_source_field(alpha, k1, msh_evl, Esrc, TE_mode)

    ! Compute exact solution Es
    print *, "Exact solution"
    call halfsphere_on_plane_scattered_field(msh_evl%POS, k1, sph_radius,  &
                                            alpha, n_terms, Eexact, TE_mode)
    print *, "done"

    ! ExNormSrc
    ExNormSrc = maxval(sqrt(abs(Eexact(:,1))**2+abs(Eexact(:,2))**2+abs(Eexact(:,3))**2))

    ! ExNormTot
    Eexact = Eexact + Esrc
    ExNormTot = maxval(sqrt(abs(Eexact(:,1))**2+abs(Eexact(:,2))**2+abs(Eexact(:,3))**2))

    ! ExNorm
    print *, "ExNormSrc", ExNormSrc
    print *, "ExNormTot", ExNormTot
end program

max_stars_repo_name: Riarrieta/WindowedGreenFunctionMethod

max_issues_repo_name: Riarrieta/WindowedGreenFunctionMethod

max_forks_repo_name: Riarrieta/WindowedGreenFunctionMethod


text:
C
C
C  Copyright (C) 2000, 2001 Silicon Graphics, Inc.  All Rights Reserved.
C
C  This program is free software; you can redistribute it and/or modify it
C  under the terms of version 2.1 of the GNU Lesser General Public License
C  as published by the Free Software Foundation.
C
C  This program is distributed in the hope that it would be useful, but
C  WITHOUT ANY WARRANTY; without even the implied warranty of
C  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
C
C  Further, this software is distributed without any warranty that it is
C  free of the rightful claim of any third person regarding infringement
C  or the like.  Any license provided herein, whether implied or
C  otherwise, applies only to this software file.  Patent licenses, if
C  any, provided herein do not apply to combinations of this program with
C  other software, or any other product whatsoever.
C
C  You should have received a copy of the GNU Lesser General Public
C  License along with this program; if not, write the Free Software
C  Foundation, Inc., 59 Temple Place - Suite 330, Boston MA 02111-1307,
C  USA.
C
C  Contact information:  Silicon Graphics, Inc., 1600 Amphitheatre Pky,
C  Mountain View, CA 94043, or:
C
C  http://www.sgi.com
C
C  For further information regarding this notice, see:
C
C  http://oss.sgi.com/projects/GenInfo/NoticeExplan
C
C

      subroutine isortb( ad, n, x, incx, id, incd )
      implicit none

*     ... Scalar arguments ...
        integer     incx, incd, n
        character*1 ad

*     ... Vector arguments ...
        integer     id( 1 )
        integer     x( * )

*     ... Local vector ...

*       For use in a call to GETTMC().  Acquire system information,
*       specifically, the number of banks.

*     ... Local variables ...
        integer ldn1, iend, iord, k, kiord
        integer i1, i2, j, kk, ll, kbd
        integer iarg, imin, imax, info
        logical asc, des
        integer amin, amax, x1, x2

*     ... External routines ...
	integer ihost_banks@

*     ... Begin Execution ...
*
CDIR$ ID "@(#) libu/sci/isortb.f	92.1	10/20/99 12:36:37"
*
        iarg = numarg()

*       Acquire banksize.  GETTMC is a system call the retreives
*       system information.  You pass it the address of the 128-word
*       array, MCTABLE, and it returns the information.  In location
*       2 of MCTABLE, is the bank size of the target system.  I
*       need this value in order to effectively govern when
*       I am going to do loop inversions.
        kbd = ihost_banks@()

*       Test input parameters
        des = ad .eq. 'D' .or. ad .eq. 'd'
        asc = ad .eq. 'A' .or. ad .eq. 'a'

        info = 0
        if( .not.des .and. .not.asc ) then
          info = 1
        elseif( n .le. 0 ) then
          info = 2
        elseif( incx .le. 0 ) then
          info = 4
        elseif( iarg .eq. 6 ) then
	  if( incd .le. 0 ) info = 6
        elseif( iarg.ne.4 ) then
          info = 7
        endif
        if( info .ne. 0 ) then
          call xersor( 'ISORTB ', info )
          goto 1500
        endif

        ldn1 = int( log( real( n ) - 0.5 ) / log( 2. ) )
*       Ascending order
        if( asc ) then
          if( iarg .eq. 4 ) goto 10
          if( iarg .eq. 6 ) goto 400

*         In-place sort
10        continue
          if( n .eq. 1 ) return

cmic$ parallel shared( incx, kbd, ldn1, n, x )
cmic$1       private ( amax, amin, i1, i2, iend, iord, j, k, kiord, kk )
cmic$2       private ( ll, x1, x2 )

*         This is the outer loop.  This governs the size of
*         the segments to be merged.
          do 300 ll = ldn1, 0, -1
            iord = 2**ll
            if( iord .ge. kbd .or. 2*iord*iord .gt. n ) goto 130
*
*           NOTE 1:
*           Suppose iord is 8.  What this doubly nested DO-loop
*           does, in effect, is that it starts out with the
*           first 8 elements.  It splits this segment in two
*           and compares each of the first 4 elements with
*           the corresponding element from the second set of
*           4 elements.  It then goes to the next segment of
*           8 elements and repeats the process until the end of
*           the array is reached.
*
cmic$ do parallel
            do 120 j = 1, iord
cdir$ ivdep
              do 110 i1 = j, n-iord, 2*iord
                i2 = i1 + iord
                x1 = x( 1+incx*(i1-1) )
                x2 = x( 1+incx*(i2-1) )
                amin = cvmgm( x2, x1, x2-x1 )
                amax = cvmgm( x1, x2, x2-x1 )

                x( 1+incx*(i1-1) ) = amin
                x( 1+incx*(i2-1) ) = amax
110           continue
120         continue
            goto 160

*
*           NOTE 2:
*           This is an inverted version of the above loop.  Loop
*           inversion helps with performance.  I always try to
*           keep the inner loops long for vectorization and
*           the outer loops short for parallel processing.
*
130         continue
cmic$ do parallel
            do 150 j = 1, n-iord, 2*iord
              iend = min( j+iord-1, n-iord )
cdir$ ivdep
              do 140 i1 = j, iend
                i2 = i1 + iord
                x1 = x( 1+incx*(i1-1) )
                x2 = x( 1+incx*(i2-1) )
                amin = cvmgm( x2, x1, x2-x1 )
                amax = cvmgm( x1, x2, x2-x1 )

                x( 1+incx*(i1-1) ) = amin
                x( 1+incx*(i2-1) ) = amax
140           continue
150         continue

*
*           NOTE 3:
*           This triply nested DO-loop performs the task of "merging"
*           the results from the previous steps.  Going back to
*           segment length of 8, it finds out how many such segments
*           exist in the array.  For the sake of convenience, suppose
*           the size of the array is 64.  That means there are 8
*           segments of length 8.  What this part of the code
*           does is to first merge the first segment and the fifth
*           segment, the second segment and the sixth, etc.  Then it
*           does the first and third, the second and fourth, etc.
*           Finally it does the first and second, the third and
*           fourth.  When doing the merge of the say the first
*           and second segments of 8, the upper 4 elements of
*           the first segment are compared-exchanged with the
*           lower 4-elements of the second segment and so it goes.
*           The illustration in Knuth is instructive.
*
160         continue
            do 230 kk = ldn1-ll, 1, -1
              k = 2**kk-1
              kiord = k*iord
              if( iord .ge. kbd .or. 2*iord*iord .gt. n ) goto 190

cmic$ do parallel
              do 180 j = 1, iord
cdir$ ivdep
                do 170 i1 = j+iord, n-kiord, 2*iord
                  i2 = i1 + kiord
                  x1 = x( 1+incx*(i1-1) )
                  x2 = x( 1+incx*(i2-1) )
                  amin = cvmgm( x2, x1, x2-x1 )
                  amax = cvmgm( x1, x2, x2-x1 )

                  x( 1+incx*(i1-1) ) = amin
                  x( 1+incx*(i2-1) ) = amax
170             continue
180           continue
              goto 230

*
*             NOTE 4:
*             Inverted version of the previous doubly nested DO-loop.
*             The rest of the code consists of variations of the
*             present code except that it may either being doing
*             the sort in descending order and/or indexed.
*
190           continue
cmic$ do parallel
              do 210 j = 1, n-kiord-iord, 2*iord
                iend = min( j+2*iord-1, n-kiord )
cdir$ ivdep
                  do 200 i1 = j+iord, iend
                  i2 = i1 + kiord
                  x1 = x( 1+incx*(i1-1) )
                  x2 = x( 1+incx*(i2-1) )
                  amin = cvmgm( x2, x1, x2-x1 )
                  amax = cvmgm( x1, x2, x2-x1 )

                  x( 1+incx*(i1-1) ) = amin
                  x( 1+incx*(i2-1) ) = amax
200             continue
210           continue
230         continue
300       continue
cmic$ end parallel
          goto 610

*         Indexed sort
400       continue
          if( n .eq. 1 ) then
            id( 1 ) = 1
            return
          endif
*
          do 405 i1 = 1, n
            id( 1+incd*(i1-1) ) = i1
405       continue
*
cmic$ parallel shared ( id, incd, incx, kbd, ldn1, n, x )
cmic$1    private ( i1, i2, iend, imax, imin, iord, j, k, kiord, kk )
cmic$2    private ( ll, x1, x2 )
          do 600 ll = ldn1, 0, -1
            iord = 2**ll
            if( iord .ge. kbd .or. 2*iord*iord .gt. n ) goto 430
*
cmic$ do parallel
            do 420 j = 1, iord
cdir$ ivdep
              do 410 i1 = j, n-iord, 2*iord
                i2 = i1 + iord
                x1 = x( 1+incx*( id( 1+incd*( i1-1 ) )-1 ) )
                x2 = x( 1+incx*( id( 1+incd*( i2-1 ) )-1 ) )

                imin = cvmgm( id( 1+incd*(i2-1) ),
     &                        id( 1+incd*(i1-1) ),
     &                        x2-x1 )

                imax = cvmgm( id( 1+incd*(i1-1) ),
     &                        id( 1+incd*(i2-1) ),
     &                        x2-x1 )

                id( 1+incd*( i1-1 ) ) = imin
                id( 1+incd*( i2-1 ) ) = imax
410           continue
420         continue
            goto 460

430         continue
cmic$ do parallel
            do 450 j = 1, n-iord, 2*iord
              iend = min( j+iord-1, n-iord )
cdir$ ivdep
              do 440 i1 = j, iend
                i2 = i1 + iord
                x1 = x( 1+incx*( id( 1+incd*( i1-1 ) )-1 ) )
                x2 = x( 1+incx*( id( 1+incd*( i2-1 ) )-1 ) )

                imin = cvmgm( id( 1+incd*(i2-1) ),
     &                        id( 1+incd*(i1-1) ),
     &                        x2-x1 )

                imax = cvmgm( id( 1+incd*(i1-1) ),
     &                        id( 1+incd*(i2-1) ),
     &                        x2-x1 )

                id( 1+incd*( i1-1 ) ) = imin
                id( 1+incd*( i2-1 ) ) = imax
440           continue
450         continue

460         continue
            do 530 kk = ldn1-ll, 1, -1
              k = 2**kk-1
              kiord = k*iord
              if( iord .ge. kbd .or. 2*iord*iord .gt. n ) goto 490

cmic$ do parallel
              do 480 j = 1, iord
cdir$ ivdep
                do 470 i1 = j+iord, n-kiord, 2*iord
                  i2 = i1 + kiord
                  x1 = x( 1+incx*( id( 1+incd*( i1-1 ) )-1 ) )
                  x2 = x( 1+incx*( id( 1+incd*( i2-1 ) )-1 ) )

                  imin = cvmgm( id( 1+incd*(i2-1) ),
     &                          id( 1+incd*(i1-1) ),
     &                          x2-x1 )

                  imax = cvmgm( id( 1+incd*(i1-1) ),
     &                          id( 1+incd*(i2-1) ),
     &                          x2-x1 )

                  id( 1+incd*( i1-1 ) ) = imin
                  id( 1+incd*( i2-1 ) ) = imax
470             continue
480           continue
              goto 530

490           continue
cmic$ do parallel
              do 510 j = 1, n-kiord-iord, 2*iord
                iend = min( j+2*iord-1, n-kiord )
cdir$ ivdep
                do 500 i1 = j+iord, iend
                  i2 = i1 + kiord
                  x1 = x( 1+incx*( id( 1+incd*( i1-1 ) )-1 ) )
                  x2 = x( 1+incx*( id( 1+incd*( i2-1 ) )-1 ) )

                  imin = cvmgm( id( 1+incd*(i2-1) ),
     &                          id( 1+incd*(i1-1) ),
     &                          x2-x1 )

                  imax = cvmgm( id( 1+incd*(i1-1) ),
     &                          id( 1+incd*(i2-1) ),
     &                          x2-x1 )

                  id( 1+incd*( i1-1 ) ) = imin
                  id( 1+incd*( i2-1 ) ) = imax
500             continue
510           continue
530         continue
600       continue
cmic$ end parallel

610       continue
        elseif( des ) then
          if( iarg .eq. 4 ) goto 710
          if( iarg .eq. 6 ) goto 1100

*         In-place sort
710       continue
cmic$ parallel shared( incx, kbd, ldn1, n, x )
cmic$1       private ( amax, amin, i1, i2, iend, iord, j, k, kiord, kk )
cmic$2       private ( ll, x1, x2 )
          do 1000 ll = ldn1, 0, -1
            iord = 2**ll
            if( iord .ge. kbd .or. 2*iord*iord .gt. n ) goto 830
*
cmic$ do parallel
            do 820 j = 1, iord
cdir$ ivdep
              do 810 i1 = j, n-iord, 2*iord
                i2 = i1 + iord
                x1 = x( 1+incx*(i1-1) )
                x2 = x( 1+incx*(i2-1) )
                amin = cvmgm( x2, x1, x2-x1 )
                amax = cvmgm( x1, x2, x2-x1 )

                x( 1+incx*(i1-1) ) = amax
                x( 1+incx*(i2-1) ) = amin
810           continue
820         continue
            goto 860

830         continue
cmic$ do parallel
            do 850 j = 1, n-iord, 2*iord
              iend = min( j+iord-1, n-iord )
cdir$ ivdep
              do 840 i1 = j, iend
                i2 = i1 + iord
                x1 = x( 1+incx*(i1-1) )
                x2 = x( 1+incx*(i2-1) )
                amin = cvmgm( x2, x1, x2-x1 )
                amax = cvmgm( x1, x2, x2-x1 )

                x( 1+incx*(i1-1) ) = amax
                x( 1+incx*(i2-1) ) = amin
840           continue
850         continue

860         continue
            do 930 kk = ldn1-ll, 1, -1
              k = 2**kk-1
              kiord = k*iord
              if( iord .ge. kbd .or. 2*iord*iord .gt. n ) goto 890

cmic$ do parallel
              do 880 j = 1, iord
cdir$ ivdep
                do 870 i1 = j+iord, n-kiord, 2*iord
                  i2 = i1 + kiord
                  x1 = x( 1+incx*(i1-1) )
                  x2 = x( 1+incx*(i2-1) )
                  amin = cvmgm( x2, x1, x2-x1 )
                  amax = cvmgm( x1, x2, x2-x1 )

                  x( 1+incx*(i1-1) ) = amax
                  x( 1+incx*(i2-1) ) = amin
870             continue
880           continue
              goto 930

890           continue
cmic$ do parallel
              do 910 j = 1, n-kiord-iord, 2*iord
                iend = min( j+2*iord-1, n-kiord )
cdir$ ivdep
                  do 900 i1 = j+iord, iend
                  i2 = i1 + kiord
                  x1 = x( 1+incx*(i1-1) )
                  x2 = x( 1+incx*(i2-1) )
                  amin = cvmgm( x2, x1, x2-x1 )
                  amax = cvmgm( x1, x2, x2-x1 )

                  x( 1+incx*(i1-1) ) = amax
                  x( 1+incx*(i2-1) ) = amin
900             continue
910           continue
930         continue
1000      continue
cmic$ end parallel
          goto 1310

*         Indexed sort
1100       continue
*
          i2 = 1
          do 1105 i1 = 1, n
            id( i2 ) = i1
            i2 = i2 + incd
1105      continue
*
cmic$ parallel shared ( id, incd, incx, kbd, ldn1, n, x )
cmic$1       private ( i1, i2, iend, imax, imin, iord, j, k, kiord, kk )
cmic$2       private ( ll, x1, x2 )
          do 1300 ll = ldn1, 0, -1
            iord = 2**ll
            if( iord .ge. kbd .or. 2*iord*iord .gt. n ) goto 1130
*
cmic$ do parallel
            do 1120 j = 1, iord
cdir$ ivdep
              do 1110 i1 = j, n-iord, 2*iord
                i2 = i1 + iord
                x1 = x( 1+incx*( id( 1+incd*( i1-1 ) )-1 ) )
                x2 = x( 1+incx*( id( 1+incd*( i2-1 ) )-1 ) )

                imin = cvmgm( id( 1+incd*(i2-1) ),
     &                        id( 1+incd*(i1-1) ),
     &                        x2-x1 )

                imax = cvmgm( id( 1+incd*(i1-1) ),
     &                        id( 1+incd*(i2-1) ),
     &                        x2-x1 )

                id( 1+incd*( i1-1 ) ) = imax
                id( 1+incd*( i2-1 ) ) = imin
1110          continue
1120        continue
            goto 1160

1130        continue
cmic$ do parallel
            do 1150 j = 1, n-iord, 2*iord
              iend = min( j+iord-1, n-iord )
cdir$ ivdep
              do 1140 i1 = j, iend
                i2 = i1 + iord
                x1 = x( 1+incx*( id( 1+incd*( i1-1 ) )-1 ) )
                x2 = x( 1+incx*( id( 1+incd*( i2-1 ) )-1 ) )

                imin = cvmgm( id( 1+incd*(i2-1) ),
     &                        id( 1+incd*(i1-1) ),
     &                        x2-x1 )

                imax = cvmgm( id( 1+incd*(i1-1) ),
     &                        id( 1+incd*(i2-1) ),
     &                        x2-x1 )

                id( 1+incd*( i1-1 ) ) = imax
                id( 1+incd*( i2-1 ) ) = imin
1140          continue
1150        continue

1160        continue
            do 1230 kk = ldn1-ll, 1, -1
              k = 2**kk-1
              kiord = k*iord
              if( iord .ge. kbd .or. 2*iord*iord .gt. n ) goto 1190

cmic$ do parallel
              do 1180 j = 1, iord
cdir$ ivdep
                do 1170 i1 = j+iord, n-kiord, 2*iord
                  i2 = i1 + kiord
                  x1 = x( 1+incx*( id( 1+incd*( i1-1 ) )-1 ) )
                  x2 = x( 1+incx*( id( 1+incd*( i2-1 ) )-1 ) )

                  imin = cvmgm( id( 1+incd*(i2-1) ),
     &                          id( 1+incd*(i1-1) ),
     &                          x2-x1 )

                  imax = cvmgm( id( 1+incd*(i1-1) ),
     &                          id( 1+incd*(i2-1) ),
     &                          x2-x1 )

                  id( 1+incd*( i1-1 ) ) = imax
                  id( 1+incd*( i2-1 ) ) = imin
1170            continue
1180          continue
              goto 1230

1190          continue
cmic$ do parallel
              do 1210 j = 1, n-kiord-iord, 2*iord
                iend = min( j+2*iord-1, n-kiord )
cdir$ ivdep
                do 1200 i1 = j+iord, iend
                  i2 = i1 + kiord
                  x1 = x( 1+incx*( id( 1+incd*( i1-1 ) )-1 ) )
                  x2 = x( 1+incx*( id( 1+incd*( i2-1 ) )-1 ) )

                  imin = cvmgm( id( 1+incd*(i2-1) ),
     &                          id( 1+incd*(i1-1) ),
     &                          x2-x1 )

                  imax = cvmgm( id( 1+incd*(i1-1) ),
     &                          id( 1+incd*(i2-1) ),
     &                          x2-x1 )

                  id( 1+incd*( i1-1 ) ) = imax
                  id( 1+incd*( i2-1 ) ) = imin
1200            continue
1210          continue
1230        continue
1300      continue
cmic$ end parallel

1310      continue
        endif

1500  continue
      return
      end

max_stars_repo_name: sharugupta/OpenUH

max_issues_repo_name: sharugupta/OpenUH

max_forks_repo_name: sharugupta/OpenUH