text:
c======================================================================
      subroutine fc2d_geoclaw_rpn2(ixy,maxm,meqn,mwaves,maux,mbc,mx,
     &      ql,qr,auxl,auxr,fwave,s,amdq,apdq)
c======================================================================
c
c Solves normal Riemann problems for the 2D SHALLOW WATER equations
c     with topography:
c     #        h_t + (hu)_x + (hv)_y = 0                           #
c     #        (hu)_t + (hu^2 + 0.5gh^2)_x + (huv)_y = -ghb_x      #
c     #        (hv)_t + (huv)_x + (hv^2 + 0.5gh^2)_y = -ghb_y      #

c On input, ql contains the state vector at the left edge of each cell
c     qr contains the state vector at the right edge of each cell
c
c This data is along a slice in the x-direction if ixy=1
c     or the y-direction if ixy=2.

c  Note that the i'th Riemann problem has left state qr(i-1,:)
c     and right state ql(i,:)
c  From the basic clawpack routines, this routine is called with
c     ql = qr
c
c
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                                                                           !
!      # This Riemann solver is for the shallow water equations.            !
!                                                                           !
!       It allows the user to easily select a Riemann solver in             !
!       riemannsolvers_geo.f. this routine initializes all the variables    !
!       for the shallow water equations, accounting for wet dry boundary    !
!       dry cells, wave speeds etc.                                         !
!                                                                           !
!           David George, Vancouver WA, Feb. 2009                           !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      use geoclaw_module, only: g => grav, drytol => dry_tolerance
      use geoclaw_module, only: earth_radius, deg2rad
      use amr_module, only: mcapa

      implicit none

      !input
      integer maxm,meqn,maux,mwaves,mbc,mx,ixy

      double precision  fwave(meqn, mwaves, 1-mbc:maxm+mbc)
      double precision  s(mwaves, 1-mbc:maxm+mbc)
      double precision  ql(meqn, 1-mbc:maxm+mbc)
      double precision  qr(meqn, 1-mbc:maxm+mbc)
      double precision  apdq(meqn,1-mbc:maxm+mbc)
      double precision  amdq(meqn,1-mbc:maxm+mbc)
      double precision  auxl(maux,1-mbc:maxm+mbc)
      double precision  auxr(maux,1-mbc:maxm+mbc)

      !local only
      integer m,i,mw,maxiter,mu,nv
      double precision wall(3)
      double precision fw(3,3)
      double precision sw(3)

      double precision hR,hL,huR,huL,uR,uL,hvR,hvL,vR,vL,phiR,phiL
      double precision bR,bL,sL,sR,sRoe1,sRoe2,sE1,sE2,uhat,chat
      double precision s1m,s2m
      double precision hstar,hstartest,hstarHLL,sLtest,sRtest
      double precision tw,dxdc

      logical rare1,rare2

      !loop through Riemann problems at each grid cell
      do i=2-mbc,mx+mbc

!-----------------------Initializing-----------------------------------
         !inform of a bad riemann problem from the start
         if((qr(1,i-1).lt.0.d0).or.(ql(1,i) .lt. 0.d0)) then
            write(*,*) 'Negative input: hl,hr,i=',qr(1,i-1),ql(1,i),i
         endif

         !Initialize Riemann problem for grid interface
         do mw=1,mwaves
              s(mw,i)=0.d0
                 fwave(1,mw,i)=0.d0
                 fwave(2,mw,i)=0.d0
                 fwave(3,mw,i)=0.d0
         enddo

c        !set normal direction
         if (ixy.eq.1) then
            mu=2
            nv=3
         else
            mu=3
            nv=2
         endif

         !zero (small) negative values if they exist
         if (qr(1,i-1).lt.0.d0) then
               qr(1,i-1)=0.d0
               qr(2,i-1)=0.d0
               qr(3,i-1)=0.d0
         endif

         if (ql(1,i).lt.0.d0) then
               ql(1,i)=0.d0
               ql(2,i)=0.d0
               ql(3,i)=0.d0
         endif

         !skip problem if in a completely dry area
         if (qr(1,i-1) <= drytol .and. ql(1,i) <= drytol) then
            go to 30
         endif

         !Riemann problem variables
         hL = qr(1,i-1)
         hR = ql(1,i)
         huL = qr(mu,i-1)
         huR = ql(mu,i)
         bL = auxr(1,i-1)
         bR = auxl(1,i)

         hvL=qr(nv,i-1)
         hvR=ql(nv,i)

         !check for wet/dry boundary
         if (hR.gt.drytol) then
            uR=huR/hR
            vR=hvR/hR
            phiR = 0.5d0*g*hR**2 + huR**2/hR
         else
            hR = 0.d0
            huR = 0.d0
            hvR = 0.d0
            uR = 0.d0
            vR = 0.d0
            phiR = 0.d0
         endif

         if (hL.gt.drytol) then
            uL=huL/hL
            vL=hvL/hL
            phiL = 0.5d0*g*hL**2 + huL**2/hL
         else
            hL=0.d0
            huL=0.d0
            hvL=0.d0
            uL=0.d0
            vL=0.d0
            phiL = 0.d0
         endif

         wall(1) = 1.d0
         wall(2) = 1.d0
         wall(3) = 1.d0
         if (hR.le.drytol) then
            call riemanntype(hL,hL,uL,-uL,hstar,s1m,s2m,
     &                                  rare1,rare2,1,drytol,g)
            hstartest=max(hL,hstar)
            if (hstartest+bL.lt.bR) then !right state should become ghost values that mirror left for wall problem
c                bR=hstartest+bL
               wall(2)=0.d0
               wall(3)=0.d0
               hR=hL
               huR=-huL
               bR=bL
               phiR=phiL
               uR=-uL
               vR=vL
            elseif (hL+bL.lt.bR) then
               bR=hL+bL
            endif
         elseif (hL.le.drytol) then ! right surface is lower than left topo
            call riemanntype(hR,hR,-uR,uR,hstar,s1m,s2m,
     &                                  rare1,rare2,1,drytol,g)
            hstartest=max(hR,hstar)
            if (hstartest+bR.lt.bL) then  !left state should become ghost values that mirror right
c               bL=hstartest+bR
               wall(1)=0.d0
               wall(2)=0.d0
               hL=hR
               huL=-huR
               bL=bR
               phiL=phiR
               uL=-uR
               vL=vR
            elseif (hR+bR.lt.bL) then
               bL=hR+bR
            endif
         endif

         !determine wave speeds
         sL=uL-sqrt(g*hL) ! 1 wave speed of left state
         sR=uR+sqrt(g*hR) ! 2 wave speed of right state

         uhat=(sqrt(g*hL)*uL + sqrt(g*hR)*uR)/(sqrt(g*hR)+sqrt(g*hL)) ! Roe average
         chat=sqrt(g*0.5d0*(hR+hL)) ! Roe average
         sRoe1=uhat-chat ! Roe wave speed 1 wave
         sRoe2=uhat+chat ! Roe wave speed 2 wave

         sE1 = min(sL,sRoe1) ! Eindfeldt speed 1 wave
         sE2 = max(sR,sRoe2) ! Eindfeldt speed 2 wave

         !--------------------end initializing...finally----------
         !solve Riemann problem.

         maxiter = 1

         call riemann_aug_JCP(maxiter,3,3,hL,hR,huL,
     &        huR,hvL,hvR,bL,bR,uL,uR,vL,vR,phiL,phiR,sE1,sE2,
     &                                    drytol,g,sw,fw)

c         call riemann_ssqfwave(maxiter,meqn,mwaves,hL,hR,huL,huR,
c     &     hvL,hvR,bL,bR,uL,uR,vL,vR,phiL,phiR,sE1,sE2,drytol,g,sw,fw)

c          call riemann_fwave(meqn,mwaves,hL,hR,huL,huR,hvL,hvR,
c     &      bL,bR,uL,uR,vL,vR,phiL,phiR,sE1,sE2,drytol,g,sw,fw)

c        !eliminate ghost fluxes for wall
         do mw=1,3
            sw(mw)=sw(mw)*wall(mw)

               fw(1,mw)=fw(1,mw)*wall(mw)
               fw(2,mw)=fw(2,mw)*wall(mw)
               fw(3,mw)=fw(3,mw)*wall(mw)
         enddo

         do mw=1,mwaves
            s(mw,i)=sw(mw)
            fwave(1,mw,i)=fw(1,mw)
            fwave(mu,mw,i)=fw(2,mw)
            fwave(nv,mw,i)=fw(3,mw)
!            write(51,515) sw(mw),fw(1,mw),fw(2,mw),fw(3,mw)
!515         format("++sw",4e25.16)
         enddo

 30      continue
      enddo

c==========Capacity for mapping from latitude longitude to physical space====
        if (mcapa.gt.0) then
         do i=2-mbc,mx+mbc
          if (ixy.eq.1) then
             dxdc=(earth_radius*deg2rad)
          else
             dxdc=earth_radius*cos(auxl(3,i))*deg2rad
          endif

          do mw=1,mwaves
c             if (s(mw,i) .gt. 316.d0) then
c               # shouldn't happen unless h > 10 km!
c                write(6,*) 'speed > 316: i,mw,s(mw,i): ',i,mw,s(mw,i)
c                endif
	           s(mw,i)=dxdc*s(mw,i)
               fwave(1,mw,i)=dxdc*fwave(1,mw,i)
               fwave(2,mw,i)=dxdc*fwave(2,mw,i)
               fwave(3,mw,i)=dxdc*fwave(3,mw,i)
          enddo
         enddo
        endif

c===============================================================================

c============= compute fluctuations=============================================
         amdq(1:3,:) = 0.d0
         apdq(1:3,:) = 0.d0
         do i=2-mbc,mx+mbc
            do  mw=1,mwaves
               if (s(mw,i) < 0.d0) then
                     amdq(1:3,i) = amdq(1:3,i) + fwave(1:3,mw,i)
               else if (s(mw,i) > 0.d0) then
                  apdq(1:3,i)  = apdq(1:3,i) + fwave(1:3,mw,i)
               else
                 amdq(1:3,i) = amdq(1:3,i) + 0.5d0 * fwave(1:3,mw,i)
                 apdq(1:3,i) = apdq(1:3,i) + 0.5d0 * fwave(1:3,mw,i)
               endif
            enddo
         enddo
!--       do i=2-mbc,mx+mbc
!--            do m=1,meqn
!--                write(51,151) m,i,amdq(m,i),apdq(m,i)
!--                write(51,152) fwave(m,1,i),fwave(m,2,i),fwave(m,3,i)
!--151             format("++3 ampdq ",2i4,2e25.15)
!--152             format("++3 fwave ",8x,3e25.15)
!--            enddo
!--        enddo

      return
      end subroutine

max_stars_repo_name: ECLAIRWaveS/ForestClaw

max_issues_repo_name: ECLAIRWaveS/ForestClaw

max_forks_repo_name: ECLAIRWaveS/ForestClaw


text:
      SUBROUTINE amebsa(p,y,mp,np,ndim,pb,yb,ftol,funk,iter,temptr)
      INTEGER iter,mp,ndim,np,NMAX
      REAL ftol,temptr,yb,p(mp,np),pb(np),y(mp),funk
      PARAMETER (NMAX=200)
      EXTERNAL funk
CU    USES amotsa,funk,ran1
      INTEGER i,idum,ihi,ilo,inhi,j,m,n
      REAL rtol,sum,swap,tt,yhi,ylo,ynhi,ysave,yt,ytry,psum(NMAX),
     *amotsa,ran1
      COMMON /ambsa/ tt,idum
      tt=-temptr
1     do 12 n=1,ndim
        sum=0.
        do 11 m=1,ndim+1
          sum=sum+p(m,n)
11      continue
        psum(n)=sum
12    continue
2     ilo=1
      inhi=1
      ihi=2
      ylo=y(1)+tt*log(ran1(idum))
      ynhi=ylo
      yhi=y(2)+tt*log(ran1(idum))
      if (ylo.gt.yhi) then
        ihi=1
        inhi=2
        ilo=2
        ynhi=yhi
        yhi=ylo
        ylo=ynhi
      endif
      do 13 i=3,ndim+1
        yt=y(i)+tt*log(ran1(idum))
        if(yt.le.ylo) then
          ilo=i
          ylo=yt
        endif
        if(yt.gt.yhi) then
          inhi=ihi
          ynhi=yhi
          ihi=i
          yhi=yt
        else if(yt.gt.ynhi) then
          inhi=i
          ynhi=yt
        endif
13    continue
      rtol=2.*abs(yhi-ylo)/(abs(yhi)+abs(ylo))
      if (rtol.lt.ftol.or.iter.lt.0) then
        swap=y(1)
        y(1)=y(ilo)
        y(ilo)=swap
        do 14 n=1,ndim
          swap=p(1,n)
          p(1,n)=p(ilo,n)
          p(ilo,n)=swap
14      continue
        return
      endif
      iter=iter-2
      ytry=amotsa(p,y,psum,mp,np,ndim,pb,yb,funk,ihi,yhi,-1.0)
      if (ytry.le.ylo) then
        ytry=amotsa(p,y,psum,mp,np,ndim,pb,yb,funk,ihi,yhi,2.0)
      else if (ytry.ge.ynhi) then
        ysave=yhi
        ytry=amotsa(p,y,psum,mp,np,ndim,pb,yb,funk,ihi,yhi,0.5)
        if (ytry.ge.ysave) then
          do 16 i=1,ndim+1
            if(i.ne.ilo)then
              do 15 j=1,ndim
                psum(j)=0.5*(p(i,j)+p(ilo,j))
                p(i,j)=psum(j)
15            continue
              y(i)=funk(psum)
            endif
16        continue
          iter=iter-ndim
          goto 1
        endif
      else
        iter=iter+1
      endif
      goto 2
      END

ext: for

max_stars_repo_path: NR-Functions/Numerical Recipes- Example & Functions/Functions/amebsa.for


text:
module operators
  use utils, only: fp, ip
  use linfuncs
  implicit none
  private

  public :: operator_cont

  type :: operator_cont
     integer(ip) :: nbr_elements
     integer(ip) :: dim
     integer(ip), allocatable,dimension(:) :: el_dim
     complex(fp), allocatable,dimension(:,:,:) :: cc,aa,nn,Iphi,Qphi,sx,sy,sz
     complex(fp), allocatable,dimension(:,:) :: eye
     contains
       procedure :: initialize
    end type operator_cont
contains

  subroutine initialize(this,el_dim_in)
    use linfuncs
    implicit none
    class(operator_cont) :: this
    integer(ip),dimension(:),intent(in) :: el_dim_in
    integer(ip) :: ii,jj,csize
    integer(ip), dimension(size(el_dim_in)+1) :: edw

    !write(*,*) 'Initializing the operators now.'
    this%nbr_elements = size(el_dim_in)
    allocate(this%el_dim(this%nbr_elements))
    this%el_dim = el_dim_in
    !write(*,*) 'el_dim:', this%el_dim
    edw(1) = 1
    edw(2:this%nbr_elements+1) = this%el_dim
    this%dim = product(this%el_dim)
    !write(*,*) 'dim:', this%dim
    !write(*,*) 'nbr_elements:', this%nbr_elements

    allocate(this%cc(this%dim,this%dim,this%nbr_elements))
    allocate(this%aa(this%dim,this%dim,this%nbr_elements))
    allocate(this%nn(this%dim,this%dim,this%nbr_elements))
    allocate(this%Iphi(this%dim,this%dim,this%nbr_elements))
    allocate(this%Qphi(this%dim,this%dim,this%nbr_elements))
    allocate(this%eye(this%dim,this%dim))
    allocate(this%sx(this%dim,this%dim,this%nbr_elements))
    allocate(this%sy(this%dim,this%dim,this%nbr_elements))
    allocate(this%sz(this%dim,this%dim,this%nbr_elements))

    !write(*,*) 'Allocated them all.'

    this%cc = (0.0_fp,0.0_fp)
    this%aa = (0.0_fp,0.0_fp)
    this%nn = (0.0_fp,0.0_fp)
    this%Iphi = (0.0_fp,0.0_fp)
    this%Qphi = (0.0_fp,0.0_fp)
    this%eye = get_eye(this%dim)
    this%sx = (0.0_fp,0.0_fp)
    this%sy = (0.0_fp,0.0_fp)

    do ii=1,this%nbr_elements
       csize = 1
       this%cc(:,:,ii) = get_eye(this%dim)
       this%aa(:,:,ii) = get_eye(this%dim)
       do jj=1,this%nbr_elements
          if(jj==ii) then
             this%cc(1:(csize*edw(jj+1)),1:(csize*edw(jj+1)),ii)&
                  =kron(this%cc(1:csize,1:csize,ii),get_creation_op(edw(jj+1)))
             this%aa(1:(csize*edw(jj+1)),1:(csize*edw(jj+1)),ii)&
                  =kron(this%aa(1:csize,1:csize,ii),get_annihilation_op(edw(jj+1)))
             !write(*,*) 'The same!'
          else
             this%cc(1:(csize*edw(jj+1)),1:(csize*edw(jj+1)),ii)&
                  =kron(this%cc(1:csize,1:csize,ii),get_eye(edw(jj+1)))
             this%aa(1:(csize*edw(jj+1)),1:(csize*edw(jj+1)),ii)&
                  =kron(this%aa(1:csize,1:csize,ii),get_eye(edw(jj+1)))
             !write(*,*) 'Different!'
          endif
          !call print_matrix(cc(1:(csize*edw(jj+1)),1:(csize*edw(jj+1)),ii),.true.)

          csize = csize*this%el_dim(jj)
       end do
    end do

    do ii=1,this%nbr_elements
       this%nn(:,:,ii) = matmul(this%cc(:,:,ii),this%aa(:,:,ii))
       this%Iphi(:,:,ii) = (0.5,0.0)*(this%cc(:,:,ii)*EXP((0.0,1.0)*0.0)&
            +this%aa(:,:,ii)*EXP(-(0.0,1.0)*0.0))
       this%Qphi(:,:,ii) = (0.0,0.5)*(this%cc(:,:,ii)*EXP((0.0,1.0)*0.0)&
            -this%aa(:,:,ii)*EXP(-(0.0,1.0)*0.0))
       !Note the added factor of 2. It seems to produce correct results, but I don't know why.
       ! The correct results to where???? I removed the factor of two, because it definitely should not be there
       this%sx(:,:,ii) = 1.0*(this%cc(:,:,ii) + this%aa(:,:,ii))
       this%sy(:,:,ii) = -(0.0_fp,1.0_fp)*(this%cc(:,:,ii) - this%aa(:,:,ii))
       this%sz(:,:,ii) = -this%nn(:,:,ii)*2+this%eye ! sz now corresponds to the standard definition
    end do
    !write(*,*) 'The final results.'
    !call print_matrix(cc(:,:,1),.true.)
    !call print_matrix(aa(:,:,1),.true.)
  end subroutine initialize

end module operators

max_stars_repo_name: ittnas/qsc

max_issues_repo_name: ittnas/qsc

max_forks_repo_name: ittnas/qsc


text:
      subroutine mpiwrapper_export_fortran_constants
      implicit none
      include "mpif.h"

      include "mpiwrapper_definitions_fortran.h"

      end

ext: f

max_stars_repo_name: ocaisa/MPIwrapper

max_issues_repo_name: ocaisa/MPIwrapper

max_forks_repo_name: ocaisa/MPIwrapper


text:
      subroutine brkpt(nbrkpt,mxint)
c******************************************************************
c                                                                 *
c  This subroutine is called from SR STMGET to read in rainfall   *
c  breakpoint data.  The first column of data read in is time     *
c  after midnight in hours(f4.2) and the second column is the     *
c  accumulated rainfall in mm(f6.2).  Irrigation is not an option *
c  when breakpoint input is used.                                 *
c                                                                 *
c******************************************************************
c                                                                 *
c  Arguments                                                      *
c    nbrkpt - number of rainfall breakpoints                      *
c    mxint  - maximum rainfall intensity (m/s)                    *
c                                                                 *
c******************************************************************
c                                                                 *
c  Argument Declarations                                          *
c                                                                 *
      integer nbrkpt
      real mxint
c******************************************************************
c                                                                 *
c  Parameters                                                     *
c      mxplan : maximum number of over land flow elements         *
c      mxtlsq : maximum number of tillage sequences               *
c      mxtill : maximum number of tillage operations per tillage  *
c                                                                 *
c******************************************************************
c
      include 'pmxelm.inc'
      include 'pmxpln.inc'
      include 'pmxtil.inc'
      include 'pmxtls.inc'
      include 'pmxhil.inc'
c
c******************************************************************
c                                                                 *
c  Common Blocks                                                  *
c                                                                 *
c******************************************************************
c
      include 'cdiss1.inc'
      include 'cdiss3.inc'
      include 'chydrol.inc'
      include 'cupdate.inc'
c     Added by S. Dun, Feb 20, 2008
c     for record storm start time for breakpoint data
      include 'cclim.inc'
c     End adding
c
c******************************************************************
c
c
c     local variables
c
      real pptcum(200), drain, dtime,pktime
      integer i
c
c     read in break points and convert into seconds and meters
      do 10 i = 1, nbrkpt
        read (13,*,end=20) timem(i), pptcum(i)
        if (i.eq.1) stmstr = timem(1)
        timem(i) = (timem(i)-stmstr) * 3600.
        pptcum(i) = pptcum(i) / 1000.
   10 continue
   20 continue
c
      stmdur = 0.0
      mxint = 0.0
      pktime = 0.0
c     compute rainfall intensity in m/s
      do 30 i = 1, nbrkpt - 1
c
        drain = pptcum(i+1) - pptcum(i)
        if (drain.eq.0.0) then
c
c XXX     Commented out following line - The time value should not be
c         zero according to Jeff Stone - even for periods of zero
c         intensity within a storm.
c         dtime = 0.0
          dtime = timem(i+1) - timem(i)
          intsty(i) = 0.0
        else if (drain.lt.0.0) then
          write (6,1000) day, mon, year
          stop
        else
          dtime = timem(i+1) - timem(i)
          if (dtime.le.0.0) then
            write (6,1100) day, mon, year
            stop
          end if
          intsty(i) = drain / dtime
        end if
c
c       compute storm duration
        stmdur = stmdur + dtime
        dur = stmdur
        if (intsty(i).gt.mxint) then
          mxint = intsty(i)
          pktime = stmdur - dtime / 2.
        end if
c     write (6,*)timem(i)/3600.,intsty(i)*3.6e+06
   30 continue
c
c     compute normalized peak intensity and time to peak
c XXX Commented out following 2 lines because these parameter values
c     are not needed in breakpoint simulations - and now their values
c     will be incorrectly estimated since STMDUR computation above
c     has been corrected.    dcf  9/19/94
c     timep = pktime / stmdur
c     ip = mxint / (pptcum(nbrkpt)/stmdur)
c
c     set rainfall intensity of last interval equal to zero
      intsty(nbrkpt) = 0.0
      prcp = pptcum(nbrkpt)
      p = prcp
c
      return
 1000 format (/1x,'*** Error in breakpoint input data. ***',//1x,
     1    'Date of precipitation : ',i2,1x,i2,1x,i4,/,1x,
     1    'The cumulative precipitation volume decreased. Check',/1x,
     1    'climate input file and restart simulation.')
 1100 format (/1x,'*** Error in breakpoint input data. ***',//1x,
     1    'Date of precipitation : ',i2,1x,i2,1x,i4,/,1x,
     1    'The breakpoint time is less than the preceeding time.',/,
     1    1x,'Check climate input file and restart simulation.')
      end

max_stars_repo_name: jarad/dep

max_issues_repo_name: jarad/dep

max_forks_repo_name: akrherz/idep


text:
module quadpack_constants

    implicit none

    ! Code working precision
    integer, parameter :: wp = kind(1.0D0)

end module quadpack_constants

max_stars_repo_name: IRFM/ALOHA

max_issues_repo_name: IRFM/ALOHA

max_forks_repo_name: IRFM/ALOHA


text:
  character(26) :: alpha
  integer :: i

  do i = 1, 26
    alpha(i:i) = achar(iachar('a') + i - 1)
  end do

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


text:
      SUBROUTINE inrgcm()
C**********************************************************************
C
C     SUBROUTINE INRGCM()
C          INitialize Random number Generator CoMmon
C
C
C                              Function
C
C
C     Initializes common area  for random number  generator.  This saves
C     the  nuisance  of  a  BLOCK DATA  routine  and the  difficulty  of
C     assuring that the routine is loaded with the other routines.
C
C**********************************************************************
C     .. Parameters ..
      INTEGER numg
      PARAMETER (numg=32)
C     ..
C     .. Scalars in Common ..
      INTEGER a1,a1vw,a1w,a2,a2vw,a2w,m1,m2
C     ..
C     .. Arrays in Common ..
      INTEGER cg1(numg),cg2(numg),ig1(numg),ig2(numg),lg1(numg),
     +        lg2(numg)
      LOGICAL qanti(numg)
C     ..
C     .. Local Scalars ..
      INTEGER i
      LOGICAL qdum
C     ..
C     .. External Functions ..
      LOGICAL qrgnsn
      EXTERNAL qrgnsn
C     ..
C     .. Common blocks ..
      COMMON /globe/m1,m2,a1,a2,a1w,a2w,a1vw,a2vw,ig1,ig2,lg1,lg2,cg1,
     +       cg2,qanti
C     ..
C     .. Save statement ..
      SAVE /globe/
C     ..
C     .. Executable Statements ..
C     V=20;                            W=30;
C
C     A1W = MOD(A1**(2**W),M1)         A2W = MOD(A2**(2**W),M2)
C     A1VW = MOD(A1**(2**(V+W)),M1)    A2VW = MOD(A2**(2**(V+W)),M2)
C
C   If V or W is changed A1W, A2W, A1VW, and A2VW need to be recomputed.
C    An efficient way to precompute a**(2*j) MOD m is to start with
C    a and square it j times modulo m using the function MLTMOD.
C
      m1 = 2147483563
      m2 = 2147483399
      a1 = 40014
      a2 = 40692
      a1w = 1033780774
      a2w = 1494757890
      a1vw = 2082007225
      a2vw = 784306273
      DO 10,i = 1,numg
          qanti(i) = .FALSE.
   10 CONTINUE
C
C     Tell the world that common has been initialized
C
      qdum = qrgnsn(.TRUE.)
      RETURN

      END

max_stars_repo_name: ppernot/MC-ChemDB

max_issues_repo_name: ppernot/MC-ChemDB

max_forks_repo_name: ppernot/MC-ChemDB


text:
** Copyright (c) 1989, NVIDIA CORPORATION.  All rights reserved.
**
** Licensed under the Apache License, Version 2.0 (the "License");
** you may not use this file except in compliance with the License.
** You may obtain a copy of the License at
**
**     http://www.apache.org/licenses/LICENSE-2.0
**
** Unless required by applicable law or agreed to in writing, software
** distributed under the License is distributed on an "AS IS" BASIS,
** WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
** See the License for the specific language governing permissions and
** limitations under the License.

*   SAVE statements in conjunction with EQUIVALENCES.

C   NONSTANDARD:
C      Use of character strings to initialize numeric variable.

	program p
	parameter(n = 10)
	save
	integer rslts(n), expect(n)

c --- tests 1 - 6:

	call sub1(.true.,  rslts(1))
	call sub1(.false., rslts(4))

c --- tests 7 - 10:

	call sub2(rslts(7))
	call sub2(rslts(9))

c --- check results:

	call check(rslts, expect, n)

	data expect / 12, 'abcd', 2, 13, 'axcd', 5,
     +                'abcd', 7, 23, 'xyzw'        /
	end

c ccccccccccccccccccccccccccccccccccccccc

	subroutine sub1(flag, ir)
	save
	equivalence (a, b), (i, j), (c(2:2), d), (c, k)
	logical flag
	integer ir(*)
	character*4 c, d*1
	integer*2 i, j

	if (flag) then
	    m = 2
	    i = 10
	    c = 'abcd'
	    b = 2.6
	else
	    i = j + 1
	    d = 'x'
	    a = a * 2.0
	endif

	ir(1) = j + m
	ir(2) = k
	ir(3) = int(a)
	end

ccccccccccccccccccccccccccccccccccccc

	subroutine sub2(ir)
	integer b, ir(*)
	equivalence (a, b), (c, d), (i, j), (m, n)
	save b, i
	character*5 j, a, c, n

	data a / 'abcde'/, i / 7 /

	ir(1) = b
	ir(2) = i
	j = 'xyzwv'
	b = 23
	end

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang


text:
SUBROUTINE fixtem(iwrit,npoin,iftmp,label)

  !***  APPLY fixities

  USE c_input
  USE ift_db
  USE curv_db, ONLY: getcun
  IMPLICIT NONE
  INTEGER (kind=4), INTENT(IN) :: iwrit,npoin,label(:)
  INTEGER (kind=4), INTENT(IN OUT) :: iftmp(:,:)

  INTEGER (kind=4) :: ifix(3),g,i,j,l,n,ipoin,chnode,np,nposn,nrve
  INTEGER (kind=4), PARAMETER :: nn = 1000000, nn1 = 1000001
  TYPE (ift_nod), POINTER :: ift
  REAL (kind=8) xg(ndoft),xf(ndoft)
  REAL (kind=8), ALLOCATABLE :: auxil(:,:)
  TYPE (rpt_set), POINTER :: rves
  TYPE (rpt_nod), POINTER :: rven

  ntfix = 0  !initializes number of fixed values

  IF(iwrit == 1) THEN
    IF(ndoft == 1) WRITE(lures,"(/,' Boundary Conditions', &
                               & /,'   Node    T')",ERR=9999)
    IF(ndoft == 2) WRITE(lures,"(/,' Boundary Conditions', &
                               & /,'   Node    TI TS')",ERR=9999)
    IF(ndoft == 3) WRITE(lures,"(/,' Boundary Conditions', &
                               & /,'   Node    TN TI TS')",ERR=9999)
  END IF

  ift => head                !point to first node in the list
  DO n=1,nift                !for each node in the list

    np    = ift%ifix(1)                 !node label
    ipoin = chnode(np)   !node internal number
    ifix(1:ndoft) = ift%ifix(2:ndoft+1) !restriction codes
    DO i=1,ndoft                        !for each DOF
      SELECT CASE (iftmp(i,ipoin))      !according to previous codes
      CASE (0)                          ! Active DOF
        IF(ifix(i) == 1) THEN           !if a restriction is included
          ntfix = ntfix+1               !increase number of fixed values
          iftmp(i,ipoin) = -ntfix       !assign a position
        END IF

      CASE (1)                      ! Not an active DOF
        IF(ifix(i) == 0) THEN           !if a release code is included
          WRITE(lures,"(' WARNING, Node ',i5,' DOF ',i2, &
                      & ' Has been released')",ERR=9999)label(ipoin),i
          iftmp(i,ipoin) = 0            !release DOF
        ELSE IF( ift%val(i) /= 0d0 )THEN
          ntfix = ntfix+1               !increase number of fixed values
          iftmp(i,ipoin) = -ntfix       !assign a position
        END IF

      CASE (:-1)                        !for a fixed or prescribed DOF
        IF(ifix(i) == 0) THEN           !if a release code included, ERROR
          WRITE(lures,"(' ERROR, Node ',i5,' DOF ',i2, &
                      & ' was previously constrained')",ERR=9999)label(ipoin),i
          CALL runend('FIXTEM: Inconsistent input data    ')
        ELSE IF(ABS(ifix(i)) == 1) THEN       !constrained twice, WARNING
          WRITE(lures,"(' WARNING, Node ',i5,' DOF ',i2, &
                      & ' was previously constrained')",ERR=9999)label(ipoin),i
        END IF
      END SELECT
    END DO
    ift => ift%next                     !point to next node
  END DO

  ALLOCATE ( fxtem(ntfix) )

  ! second loop to gather fixed temperatures

  ift => head                !point to first node in the list
  DO n=1,nift                !for each node in the list

    np  = ift%ifix(1)        !node label
    ipoin = chnode(np)       !node internal number
    ifix(1:ndoft) = ift%ifix(2:ndoft+1) !restriction codes
    DO i=1,ndoft                        !for each DOF
      IF( iftmp(i,ipoin) < 0 )THEN
        np = -iftmp(i,ipoin)
        fxtem(np) = ift%val(i)
      END IF
    END DO
    IF(iwrit == 1) THEN
      IF(ndoft == 1)WRITE(lures,"(i7,3x,1i2,1f12.4)",ERR=9999)ift%ifix(1:2),ift%val(1)
      IF(ndoft == 2)WRITE(lures,"(i7,3x,2i3,2f12.4)",ERR=9999)ift%ifix(1:3),ift%val(1:2)
      IF(ndoft == 3)WRITE(lures,"(i7,3x,3i3,3f12.4)",ERR=9999)ift%ifix(1:4),ift%val(1:3)
    END IF
    ift => ift%next                     !point to next node
  END DO

  !***  Applies prescribed temperatures

  nprev = 0    !initializes

  IF (npret > 0) THEN  !if prescribed temperatures exist

    ALLOCATE( lctmp(npret+1) )                !space for associated curves
    ALLOCATE( auxil(0:ndoft*npoin,npret) )    !auxiliar space
    auxil = 0d0                               !initializes

    rves => headv                             !point to first set

    DO l=1,npret                   !loop on each set

      lctmp(l)  = getcun (rves%lc)  ! get associated curve number for this set
      auxil(0,l)= rves%factor       !associated factor for this set

      IF(iwrit == 1) THEN
        WRITE(lures, "(//, &
        & 5X,'Curve scaling this temperature set .',i10/ &
        & 5X,'Scaling factor .....................',e14.7,/)",ERR=9999) &
          lctmp(l),auxil(0,l)

        WRITE(lures, "(//,5X,'Prescribed temperatures',//)",ERR=9999)
        IF (ndoft == 1 ) WRITE(lures, &
           "(6X,'NODE',6X,'Temp.')",ERR=9999)
        IF (ndoft == 2 ) WRITE(lures, &
           "(6X,'NODE',6X,'Temp-S',8X,'Temp-I')",ERR=9999)
        IF (ndoft == 3 ) WRITE(lures, &
           "(6X,'NODE',6X,'Temp-M',8X,'Temp-I',8X,'Temp-S')",ERR=9999)
      END IF

      nrve = rves%nrv         !number of nodes in this set
      rven => rves%head       !point to first node in the set
      DO j=1,nrve             !loop for each node in the set
        g = rven%node                  !node label
        n = chnode(g)         !internal number

        xg(1:ndoft) = rven%v(1:ndoft)  !prescribed temperatures components
        xf = xg                        !copy onto XF
        DO i=1,ndoft                   !for each DOF
          nposn = iftmp(i,n)           !restriction code
          SELECT CASE (nposn)
          CASE (-nn:-1)                       !fixed values
            ! this is an error
            WRITE (lures, "(' Warning. At node',i8,' previously fixed ', &
                          & 'boundary conditions overwritten.')",ERR=9999)  g
          CASE (:-nn1)                        !prescribed values
            auxil(-nposn+nn,l) = xg(i)        !assign value
          CASE (0,1)                          !active DOF or non existent
            IF(xg(i) /= 0d0)THEN              !if a non-zero value
              nprev = nprev+1                 !increase number of fixed values
              iftmp(i,n) = -nprev-nn          !modify restriction code
              auxil(nprev,l) = xg(i)          !assign value
            END IF
          END SELECT
        END DO
        ! echo effectively assigned values
        IF(iwrit == 1)  WRITE(lures,"(i10,3e14.5)",ERR=9999) label(n), xf(1:ndoft)
        rven => rven%next        !point to next node
      END DO
      rves => rves%next       !point to next set
    END DO

    ALLOCATE( prtmp(nprev+1,npret+1) )     !reserve space for definitive array
    prtmp(1:nprev,1:npret) = auxil(1:nprev,1:npret)  !transfer temperature data
    prtmp(nprev+1,1:npret) = auxil(0,1:npret)        !transfer factors
    prtmp(1:nprev,npret+1) = 0.d0  !initializes compounded prescribed temperatures
    DEALLOCATE ( auxil )           !release auxiliar array

  ELSE  !no prescribed temperatures

    ALLOCATE( prtmp(1,1) ) !to avoid null pointers only
    ALLOCATE( lctmp(npret+1) )
    prtmp = 0d0

  END IF

RETURN
 9999 CALL runen2('')
END SUBROUTINE fixtem

max_stars_repo_name: jerebenitez/IFE-simpact-openfoam

max_issues_repo_name: jerebenitez/IFE-simpact-openfoam

max_forks_repo_name: jerebenitez/IFE-simpact-openfoam


module mod_misc
  !-----------------------------------------------------------------------------
  !
  !++ Used modules
  !
  use prec
  use mod_precision
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public, including parameters
  !
  include 'problem_size.inc'

  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: DEBUG_rapstart
  public :: DEBUG_rapend
  public :: DEBUG_rapreport
  public :: DEBUG_valuecheck

  interface DEBUG_valuecheck
     module procedure DEBUG_valuecheck_1D
     module procedure DEBUG_valuecheck_2D
     module procedure DEBUG_valuecheck_3D
     module procedure DEBUG_valuecheck_4D
     module procedure DEBUG_valuecheck_5D
     module procedure DEBUG_valuecheck_6D
  end interface DEBUG_valuecheck

  public :: MISC_make_idstr        !--- make file name with a number
  public :: MISC_get_available_fid !--- get an available file ID

  public :: ADM_proc_stop

  public :: CONST_setup

  public :: trace_start, trace_end
  public :: wait_message, test_message

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  integer,             public :: EX_CSTEP_diffusion = 0
  integer,             public :: EX_TSTEP_diffusion = 60
  integer,             public :: EX_CSTEP_divdamp3d = 0
  integer,             public :: EX_TSTEP_divdamp3d = 140
  integer,             public :: EX_fid
  integer,             public :: EX_err
  character(len=1024), public :: EX_fname
  character(len=16),   public :: EX_item
  real(RP),            public :: EX_max
  real(RP),            public :: EX_min
  real(RP),            public :: EX_sum

  type t_message
    !! dummy here
  end type t_message
  public t_message

  REAL(rstd), public :: ptop ! pressure at top of atmosphere l=llm+1
!!$  REAL(rstd), public, allocatable :: mass_ak(:) ! llm
!!$  REAL(rstd), public, allocatable :: mass_bk(:) ! llm

!!$  REAL(rstd), public, allocatable  :: le(:)          ! lenght of a edge
!!$  REAL(rstd), public, allocatable  :: Ai(:)          ! area of a cell
!!$
!!$  REAL(rstd), public, allocatable :: Av(:) ! area of dual mesk cell
!!$  real(rstd), public, allocatable :: de(:) ! distance from a neighbour == lenght of an edge of the dual mesh
!!$  real(rstd), public, allocatable :: Riv2(:,:) ! weigt
!!$  REAL(rstd), public, allocatable :: Wee(:,:,:)     ! weight
!!$  real(rstd), public, allocatable :: fv(:) ! coriolis (evaluted on a vertex)

  REAL(rstd), public, allocatable  :: bp(:)

  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  private :: DEBUG_rapid

  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  integer,                 private, parameter :: DEBUG_rapnlimit = 100
  integer,                 private            :: DEBUG_rapnmax   = 0
  character(len=ADM_NSYS), private            :: DEBUG_rapname(DEBUG_rapnlimit)
  real(8),                 private            :: DEBUG_raptstr(DEBUG_rapnlimit)
  real(8),                 private            :: DEBUG_rapttot(DEBUG_rapnlimit)
  integer,                 private            :: DEBUG_rapnstr(DEBUG_rapnlimit)
  integer,                 private            :: DEBUG_rapnend(DEBUG_rapnlimit)

!!$#ifdef _FIXEDINDEX_
!!$  real(DP), public              :: GRD_gz   (ADM_kall)
!!$  real(DP), public              :: GRD_gzh  (ADM_kall)
!!$  real(DP), public              :: GRD_dgz  (ADM_kall)
!!$  real(DP), public              :: GRD_dgzh (ADM_kall)
!!$  real(DP), public              :: GRD_rdgz (ADM_kall)
!!$  real(DP), public              :: GRD_rdgzh(ADM_kall)
!!$
!!$  real(DP), public              :: GRD_afact(ADM_kall)
!!$  real(DP), public              :: GRD_bfact(ADM_kall)
!!$  real(DP), public              :: GRD_cfact(ADM_kall)
!!$  real(DP), public              :: GRD_dfact(ADM_kall)
!!$
!!$  real(RP), public              :: GRD_xr   (ADM_gall   ,ADM_kall,ADM_lall   ,ADM_AI:ADM_AJ,XDIR:ZDIR)
!!$  real(RP), public              :: GRD_xr_pl(ADM_gall_pl,ADM_kall,ADM_lall_pl,              XDIR:ZDIR)
!!$
!!$  real(RP), public              :: GMTR_p   (ADM_gall   ,ADM_KNONE,ADM_lall   ,              GMTR_p_nmax   )
!!$  real(RP), public              :: GMTR_p_pl(ADM_gall_pl,ADM_KNONE,ADM_lall_pl,              GMTR_p_nmax   )
!!$  real(RP), public              :: GMTR_t   (ADM_gall   ,ADM_KNONE,ADM_lall   ,ADM_TI:ADM_TJ,GMTR_t_nmax   )
!!$  real(RP), public              :: GMTR_t_pl(ADM_gall_pl,ADM_KNONE,ADM_lall_pl,              GMTR_t_nmax   )
!!$  real(RP), public              :: GMTR_a   (ADM_gall   ,ADM_KNONE,ADM_lall   ,ADM_AI:ADM_AJ,GMTR_a_nmax   )
!!$  real(RP), public              :: GMTR_a_pl(ADM_gall_pl,ADM_KNONE,ADM_lall_pl,              GMTR_a_nmax_pl)
!!$
!!$  real(RP), public              :: VMTR_RGSQRTH     (ADM_gall   ,ADM_kall,ADM_lall   )
!!$  real(RP), public              :: VMTR_RGSQRTH_pl  (ADM_gall_pl,ADM_kall,ADM_lall_pl)
!!$  real(RP), public              :: VMTR_RGAM        (ADM_gall   ,ADM_kall,ADM_lall   )
!!$  real(RP), public              :: VMTR_RGAM_pl     (ADM_gall_pl,ADM_kall,ADM_lall_pl)
!!$  real(RP), public              :: VMTR_RGAMH       (ADM_gall   ,ADM_kall,ADM_lall   )
!!$  real(RP), public              :: VMTR_RGAMH_pl    (ADM_gall_pl,ADM_kall,ADM_lall_pl)
!!$  real(RP), public              :: VMTR_RGSGAM2     (ADM_gall   ,ADM_kall,ADM_lall   )
!!$  real(RP), public              :: VMTR_RGSGAM2_pl  (ADM_gall_pl,ADM_kall,ADM_lall_pl)
!!$  real(RP), public              :: VMTR_RGSGAM2H    (ADM_gall   ,ADM_kall,ADM_lall   )
!!$  real(RP), public              :: VMTR_RGSGAM2H_pl (ADM_gall_pl,ADM_kall,ADM_lall_pl)
!!$  real(RP), public              :: VMTR_GSGAM2H     (ADM_gall   ,ADM_kall,ADM_lall   )
!!$  real(RP), public              :: VMTR_GSGAM2H_pl  (ADM_gall_pl,ADM_kall,ADM_lall_pl)
!!$  real(RP), public              :: VMTR_C2WfactGz   (ADM_gall   ,ADM_kall,6,ADM_lall   )
!!$  real(RP), public              :: VMTR_C2WfactGz_pl(ADM_gall_pl,ADM_kall,6,ADM_lall_pl)
!!$#else
!!$  real(DP), public, allocatable :: GRD_gz   (:) ! gsi (z-star) coordinate
!!$  real(DP), public, allocatable :: GRD_gzh  (:) ! gsi (z-star) coordinate at the half point
!!$  real(DP), public, allocatable :: GRD_dgz  (:) ! d(gsi)
!!$  real(DP), public, allocatable :: GRD_dgzh (:) ! d(gsi) at the half point
!!$  real(DP), public, allocatable :: GRD_rdgz (:)
!!$  real(DP), public, allocatable :: GRD_rdgzh(:)
!!$
!!$  real(DP), public, allocatable :: GRD_afact(:) ! From the cell center value to the cell wall value
!!$  real(DP), public, allocatable :: GRD_bfact(:) !    A(k-1/2) = ( afac(k) A(k) + bfac(k) * A(k-1) ) / 2
!!$  real(DP), public, allocatable :: GRD_cfact(:) ! From the cell wall value to the cell center value
!!$  real(DP), public, allocatable :: GRD_dfact(:) !    A(k) = ( cfac(k) A(k+1/2) + dfac(k) * A(k-1/2) ) / 2
!!$
!!$  real(RP), public, allocatable :: GRD_xr   (:,:,:,:,:) ! mass centroid position
!!$  real(RP), public, allocatable :: GRD_xr_pl(:,:,:,:)
!!$
!!$  real(RP), public, allocatable :: GMTR_p   (:,:,:,:)   ! geometrics for the cell point
!!$  real(RP), public, allocatable :: GMTR_p_pl(:,:,:,:)
!!$  real(RP), public, allocatable :: GMTR_t   (:,:,:,:,:) ! geometrics for the cell vertex
!!$  real(RP), public, allocatable :: GMTR_t_pl(:,:,:,:)
!!$  real(RP), public, allocatable :: GMTR_a   (:,:,:,:,:) ! geometrics for the cell arc
!!$  real(RP), public, allocatable :: GMTR_a_pl(:,:,:,:)
!!$
!!$  real(RP), public, allocatable :: VMTR_RGSQRTH     (:,:,:)
!!$  real(RP), public, allocatable :: VMTR_RGSQRTH_pl  (:,:,:)
!!$  real(RP), public, allocatable :: VMTR_RGAM        (:,:,:)
!!$  real(RP), public, allocatable :: VMTR_RGAM_pl     (:,:,:)
!!$  real(RP), public, allocatable :: VMTR_RGAMH       (:,:,:)
!!$  real(RP), public, allocatable :: VMTR_RGAMH_pl    (:,:,:)
!!$  real(RP), public, allocatable :: VMTR_RGSGAM2     (:,:,:)
!!$  real(RP), public, allocatable :: VMTR_RGSGAM2_pl  (:,:,:)
!!$  real(RP), public, allocatable :: VMTR_RGSGAM2H    (:,:,:)
!!$  real(RP), public, allocatable :: VMTR_RGSGAM2H_pl (:,:,:)
!!$  real(RP), public, allocatable :: VMTR_GSGAM2H     (:,:,:)
!!$  real(RP), public, allocatable :: VMTR_GSGAM2H_pl  (:,:,:)
!!$  real(RP), public, allocatable :: VMTR_C2WfactGz   (:,:,:,:)
!!$  real(RP), public, allocatable :: VMTR_C2WfactGz_pl(:,:,:,:)
!!$#endif

  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  function DEBUG_rapid( rapname ) result(id)
    implicit none

    character(len=*), intent(in) :: rapname

    integer :: id
    !---------------------------------------------------------------------------

    if ( DEBUG_rapnmax >= 1 ) then
       do id = 1, DEBUG_rapnmax
          if( trim(rapname) == trim(DEBUG_rapname(id)) ) return
       enddo
    endif

    DEBUG_rapnmax     = DEBUG_rapnmax + 1
    id                = DEBUG_rapnmax
    DEBUG_rapname(id) = trim(rapname)
    DEBUG_raptstr(id) = 0.D0
    DEBUG_rapttot(id) = 0.D0
    DEBUG_rapnstr(id) = 0
    DEBUG_rapnend(id) = 0

  end function DEBUG_rapid

  !-----------------------------------------------------------------------------
  subroutine DEBUG_rapstart( rapname )
    implicit none

    character(len=*), intent(in) :: rapname

    real(8) :: time
    integer :: id
    !---------------------------------------------------------------------------

    id = DEBUG_rapid( rapname )

    call CPU_TIME(time)

    DEBUG_raptstr(id) = time
    DEBUG_rapnstr(id) = DEBUG_rapnstr(id) + 1

    !write(ADM_LOG_FID,*) rapname, DEBUG_rapnstr(id)

#ifdef _FAPP_
    call fapp_start( rapname, id, 1 )
#endif

    return
  end subroutine DEBUG_rapstart

  !-----------------------------------------------------------------------------
  subroutine DEBUG_rapend( rapname )
    implicit none

    character(len=*), intent(in) :: rapname

    real(8) :: time
    integer :: id
    !---------------------------------------------------------------------------

    id = DEBUG_rapid( rapname )

    call CPU_TIME(time)

    DEBUG_rapttot(id) = DEBUG_rapttot(id) + ( time-DEBUG_raptstr(id) )
    DEBUG_rapnend(id) = DEBUG_rapnend(id) + 1

#ifdef _FAPP_
    call fapp_stop( rapname, id, 1 )
#endif

    return
  end subroutine DEBUG_rapend

  !-----------------------------------------------------------------------------
  subroutine DEBUG_rapreport
    implicit none

    integer :: id
    !---------------------------------------------------------------------------

    if ( DEBUG_rapnmax >= 1 ) then

       do id = 1, DEBUG_rapnmax
          if ( DEBUG_rapnstr(id) /= DEBUG_rapnend(id) ) then
              write(*,*) '*** Mismatch Report',id,DEBUG_rapname(id),DEBUG_rapnstr(id),DEBUG_rapnend(id)
          endif
       enddo

       write(ADM_LOG_FID,*)
       write(ADM_LOG_FID,*) '*** Computational Time Report'

       do id = 1, DEBUG_rapnmax
          write(ADM_LOG_FID,'(1x,A,I3.3,A,A,A,F10.3,A,I7)') &
          '*** ID=',id,' : ',DEBUG_rapname(id),' T=',DEBUG_rapttot(id),' N=',DEBUG_rapnstr(id)
       enddo

    else
       write(ADM_LOG_FID,*)
       write(ADM_LOG_FID,*) '*** Computational Time Report: NO item.'
    endif

    return
  end subroutine DEBUG_rapreport

  !-----------------------------------------------------------------------------
  subroutine DEBUG_valuecheck_1D( &
       varname, &
       var      )
    implicit none

    character(len=*),  intent(in)  :: varname
    real(RP),          intent(in)  :: var(:)
    !---------------------------------------------------------------------------

     EX_item = trim  (varname)
     EX_max  = maxval(var)
     EX_min  = minval(var)
     EX_sum  = sum   (var)
     write(ADM_LOG_FID,'(1x,A,A16,3(A,ES24.16))') '+check[',EX_item,'] max=',EX_max,',min=',EX_min,',sum=',EX_sum

    return
  end subroutine DEBUG_valuecheck_1D

  !-----------------------------------------------------------------------------
  subroutine DEBUG_valuecheck_2D( &
       varname, &
       var      )
    implicit none

    character(len=*),  intent(in)  :: varname
    real(RP),          intent(in)  :: var(:,:)
    !---------------------------------------------------------------------------

     EX_item = trim  (varname)
     EX_max  = maxval(var)
     EX_min  = minval(var)
     EX_sum  = sum   (var)
     write(ADM_LOG_FID,'(1x,A,A16,3(A,ES24.16))') '+check[',EX_item,'] max=',EX_max,',min=',EX_min,',sum=',EX_sum

    return
  end subroutine DEBUG_valuecheck_2D

  !-----------------------------------------------------------------------------
  subroutine DEBUG_valuecheck_3D( &
       varname, &
       var      )
    implicit none

    character(len=*),  intent(in)  :: varname
    real(RP),          intent(in)  :: var(:,:,:)
    !---------------------------------------------------------------------------

     EX_item = trim  (varname)
     EX_max  = maxval(var)
     EX_min  = minval(var)
     EX_sum  = sum   (var)
     write(ADM_LOG_FID,'(1x,A,A16,3(A,ES24.16))') '+check[',EX_item,'] max=',EX_max,',min=',EX_min,',sum=',EX_sum

    return
  end subroutine DEBUG_valuecheck_3D

  !-----------------------------------------------------------------------------
  subroutine DEBUG_valuecheck_4D( &
       varname, &
       var      )
    implicit none

    character(len=*),  intent(in)  :: varname
    real(RP),          intent(in)  :: var(:,:,:,:)
    !---------------------------------------------------------------------------

     EX_item = trim  (varname)
     EX_max  = maxval(var)
     EX_min  = minval(var)
     EX_sum  = sum   (var)
     write(ADM_LOG_FID,'(1x,A,A16,3(A,ES24.16))') '+check[',EX_item,'] max=',EX_max,',min=',EX_min,',sum=',EX_sum

    return
  end subroutine DEBUG_valuecheck_4D

  !-----------------------------------------------------------------------------
  subroutine DEBUG_valuecheck_5D( &
       varname, &
       var      )
    implicit none

    character(len=*),  intent(in)  :: varname
    real(RP),          intent(in)  :: var(:,:,:,:,:)
    !---------------------------------------------------------------------------

     EX_item = trim  (varname)
     EX_max  = maxval(var)
     EX_min  = minval(var)
     EX_sum  = sum   (var)
     write(ADM_LOG_FID,'(1x,A,A16,3(A,ES24.16))') '+check[',EX_item,'] max=',EX_max,',min=',EX_min,',sum=',EX_sum

    return
  end subroutine DEBUG_valuecheck_5D

  !-----------------------------------------------------------------------------
  subroutine DEBUG_valuecheck_6D( &
       varname, &
       var      )
    implicit none

    character(len=*),  intent(in)  :: varname
    real(RP),          intent(in)  :: var(:,:,:,:,:,:)
    !---------------------------------------------------------------------------

     EX_item = trim  (varname)
     EX_max  = maxval(var)
     EX_min  = minval(var)
     EX_sum  = sum   (var)
     write(ADM_LOG_FID,'(1x,A,A16,3(A,ES24.16))') '+check[',EX_item,'] max=',EX_max,',min=',EX_min,',sum=',EX_sum

    return
  end subroutine DEBUG_valuecheck_6D

  !-----------------------------------------------------------------------------
  subroutine MISC_make_idstr( &
       str,    &
       prefix, &
       ext,    &
       numID,  &
       digit   )
    implicit none

    character(len=*),  intent(out) :: str    !< combined extention string
    character(len=*),  intent(in)  :: prefix !< prefix
    character(len=*),  intent(in)  :: ext    !< extention ( e.g. .rgn )
    integer,           intent(in)  :: numID  !< number
    integer, optional, intent(in)  :: digit  !< digit

    logical, parameter            :: NSTR_ZERO_START = .true. ! number of separated file starts from 0 ?
    integer, parameter            :: NSTR_MAX_DIGIT  = 6      ! digit of separated file

    character(len=128) :: rankstr
    integer            :: setdigit
    !---------------------------------------------------------------------------

    if ( NSTR_ZERO_START ) then
       write(rankstr,'(I128.128)') numID-1
    else
       write(rankstr,'(I128.128)') numID
    endif

    if ( present(digit) ) then
       setdigit = digit
    else
       setdigit = NSTR_MAX_DIGIT
    endif

    rankstr(1:setdigit) = rankstr(128-(setdigit-1):128)
    rankstr(setdigit+1:128) = ' '

    str = trim(prefix)//'.'//trim(ext)//trim(rankstr) ! -> prefix.ext00000

    return
  end subroutine MISC_make_idstr

  !-----------------------------------------------------------------------------
  !> Search and get available machine id
  !> @return fid
  function MISC_get_available_fid() result(fid)
    implicit none

    integer :: fid

    integer, parameter :: min_fid =  7 !< minimum available fid
    integer, parameter :: max_fid = 99 !< maximum available fid

    logical :: i_opened
    !---------------------------------------------------------------------------

    do fid = min_fid, max_fid
       inquire(fid,opened=i_opened)
       if( .NOT. i_opened ) return
    enddo

  end function MISC_get_available_fid

  !-----------------------------------------------------------------------------
  subroutine ADM_proc_stop

    stop

  end subroutine ADM_proc_stop

  !-----------------------------------------------------------------------------
  subroutine CONST_setup
    implicit none

    write(ADM_LOG_FID,*) '*** setup constants'

    CONST_PI   = 4.E0_RP * atan( 1.0_RP )
    CONST_EPS  = epsilon(0.0_RP)
    CONST_HUGE =    huge(0.0_RP)

    write(ADM_LOG_FID,*) '*** PI   = ', CONST_PI
    write(ADM_LOG_FID,*) '*** EPS  = ', CONST_EPS
    write(ADM_LOG_FID,*) '*** HUGE = ', CONST_HUGE

    return
  end subroutine CONST_setup

  !! from trace.F90
  SUBROUTINE trace_start(name)
  IMPLICIT NONE
    CHARACTER(LEN=*),INTENT(IN) :: name
#ifdef VTRACE
#include <vt_user.inc>
#endif

!$OMP MASTER
#ifdef VTRACE
     VT_USER_START(name)
#endif
!$OMP END MASTER

  END SUBROUTINE trace_start

  !! from trace.F90
  SUBROUTINE trace_end(name)
  IMPLICIT NONE
#ifdef VTRACE
#include <vt_user.inc>
#endif

    CHARACTER(LEN=*),INTENT(IN) :: name

!$OMP MASTER
#ifdef VTRACE
     VT_USER_END(name)
#endif
!$OMP END MASTER

  END SUBROUTINE trace_end

  !! from transfert_mpi.f90 with _seq version.
  SUBROUTINE wait_message(message)
  IMPLICIT NONE
    TYPE(t_message) :: message

  END SUBROUTINE wait_message

  !! from transfert_mpi.f90 with _seq version.
  SUBROUTINE test_message(message)
  IMPLICIT NONE
    TYPE(t_message) :: message
  END SUBROUTINE  test_message

end module mod_misc

!!$ Local Variables:
!!$ f90-do-indent: 2
!!$ f90-if-indent: 2
!!$ f90-program-indent: 2
!!$ f90-type-indent: 2
!!$ f90-directive-comment-re: "!\$omp "
!!$ end:


text:
      program transp_mn
      implicit none
c
c --- dummy version.
c
      write(6,'(a,a)')
     &  'use transp_mn3, or transp_mn3_lm, ',
     &  'with flnm_p and flnm_n set to NONE'
      stop
      end

max_stars_repo_name: TillRasmussen/HYCOM-tools

max_issues_repo_name: TillRasmussen/HYCOM-tools

max_forks_repo_name: TillRasmussen/HYCOM-tools


text:
      module basis
        ! Atomic basis set and integral matrices as a wrapper of "eval".
!$      use omp_lib
        use eval
        implicit none
        private
        integer,parameter::MXB=1000,MXC=20
        integer::nb_,ls_(3,MXB),is_(MXB),ns_(MXB)
        real(8)::cs_(MXC,MXB),es_(MXC,MXB)
        ! --------------------------------------------------------------
        public ::
     &    basis_init,
     &    basis_add,
     &    basis_set,
     &    basis_dim,
     &    basis_sij,
     &    basis_hij,
     &    basis_dijk,
     &    basis_ijkl,
     &    basis_enuc
        ! --------------------------------------------------------------

      contains

        subroutine basis_init()
          implicit none
          nb_=0
          ls_(:,:)=0
          is_(:)=0
          ns_(:)=0
          cs_(:,:)=0d0
          es_(:,:)=0d0
        end subroutine

        subroutine basis_dim(num)
          ! The number of basis set.
          implicit none
          integer,intent(out)::num
          num=nb_
        end subroutine

        subroutine basis_set(key,ia,i1)
          ! Add an atomic basis set.
          implicit none
          character(*),intent(in)::key
          integer,intent(in)::ia,i1
          if (key=='STO3G') then
            call bs_sto3g_(ia,i1)
          else if (key=='321G') then
            call bs_321g_(ia,i1)
          else
            write(*,'("error: no such basis set: ",a)') TRIM(key)
            stop 1
          end if
        end subroutine

        subroutine basis_add(i1,n1,l1,c1,e1)
          ! Add Gaussian basis set to the list.
          implicit none
          integer,intent(in)::i1,n1,l1(3)
          real(8),intent(in)::c1(n1),e1(n1)
          nb_=nb_+1
          if (nb_>MXB) then
            write(*,'("error: out of range: basis_add()")')
            stop 1
          end if
          ns_(nb_)=n1
          is_(nb_)=i1
          ls_(:3,nb_)=l1(:3)
          cs_(:n1,nb_)=c1(:n1)
          es_(:n1,nb_)=e1(:n1)
        end subroutine

        subroutine basis_enuc(r,ian,val)
          ! Nuclear repulsion energy.
          integer,intent(in)::ian(:)
          real(8),intent(in)::r(:,:)
          real(8),intent(out)::val
          real(8)::rij
          integer::i,j
          val=0d0
          do i=1,SIZE(ian)
            do j=1,i-1
              rij=SQRT(SUM((r(:,i)-r(:,j))**2))
              val=val+ian(i)*ian(j)/rij
            end do
          end do
        end subroutine

        subroutine basis_sij(r,val)
          ! AO overlap matrix (i|j).
          implicit none
          real(8),intent(in)::r(:,:)
          real(8),intent(out),allocatable::val(:)
          integer::ij,i,j,nb2
          real(8)::vt
          ! Get memory.
          nb2=ip_(nb_,nb_)
          allocate(val(nb2))
          ! Evaluate all the symmetry-unique matrix elements.
          do ij=1,SIZE(val)
            call iup_(ij,i,j)
            if (i==j) then
              val(ij)=1d0
            else
              call eval_overlap(
     &          ns_(i),ns_(j),cs_(:,i),cs_(:,j),es_(:,i),es_(:,j),
     &          r(:,is_(i)),r(:,is_(j)),ls_(:,i),ls_(:,j),vt)
              val(ij)=vt
            end if
          end do
        end subroutine

        subroutine basis_hij(r,ian,val)
          ! AO core Hamiltonian matrix (i|h|j).
          implicit none
          integer,intent(in)::ian(:)
          real(8),intent(in)::r(:,:)
          real(8),intent(out),allocatable::val(:)
          integer::ij,i,j,ic,nb2
          real(8)::vt
          ! Get memory.
          nb2=ip_(nb_,nb_)
          allocate(val(nb2))
          ! Evaluate all the symmetry-unique matrix elements.
          do ij=1,SIZE(val)
            call iup_(ij,i,j)
            call eval_kinetic(
     &        ns_(i),ns_(j),cs_(:,i),cs_(:,j),es_(:,i),es_(:,j),
     &        r(:,is_(i)),r(:,is_(j)),ls_(:,i),ls_(:,j),vt)
            val(ij)=vt
            do ic=1,SIZE(ian)
              call eval_attract(
     &          ns_(i),ns_(j),cs_(:,i),cs_(:,j),es_(:,i),es_(:,j),
     &          r(:,is_(i)),r(:,is_(j)),ls_(:,i),ls_(:,j),r(:,ic),vt)
              val(ij)=val(ij)-ian(ic)*vt
            end do
          end do
        end subroutine

        subroutine basis_dijk(r,val)
          ! AO derivative coupling matrix (i|d/dRk|j)
          implicit none
          real(8),intent(in)::r(:,:)
          real(8),intent(out),allocatable::val(:,:,:,:)
          real(8),allocatable::vb(:,:,:)
          integer::i,j,ix,ia,natm
          real(8)::vt
          natm=SIZE(r(1,:))
          allocate(val(nb_,nb_,3,natm))
          allocate(vb(nb_,nb_,3))
          ! AO velocity integral (i|d/dr|j)
          do ix=1,3
            do i=1,nb_
              do j=1,nb_
                if (i==j) then
                  vt=0d0
                else
                  call eval_velocity(
     &              ns_(i),ns_(j),cs_(:,i),cs_(:,j),es_(:,i),es_(:,j),
     &              r(:,is_(i)),r(:,is_(j)),ls_(:,i),ls_(:,j),ix,vt)
                end if
                vb(i,j,ix)=vt
              end do
            end do
          end do
          val(:,:,:,:)=0d0
          do ia=1,natm
            do ix=1,3
              do i=1,nb_
                do j=1,nb_
                  if (is_(j)==ia) then
                    val(i,j,ix,ia)=(-1d0)*vb(i,j,ix)
                  end if
                end do
              end do
            end do
          end do
        end subroutine

        subroutine basis_ijkl(r,val)
          ! AO 2e integral matrix (ij|kl).
          implicit none
          real(8),intent(in)::r(:,:)
          real(8),intent(out),allocatable::val(:)
          integer::ijkl,i,j,k,l,nb3
          real(8)::vt
          ! Get memory.
          nb3=ip2_(nb_,nb_,nb_,nb_)
          allocate(val(nb3))
          ! Evaluate all the symmetry-unique matrix elements.
!$omp     parallel private(ijkl,i,j,k,l,vt)
!$omp     do
          do ijkl=1,SIZE(val)
            call iup2_(ijkl,i,j,k,l)
            call eval_repulse(
     &        ns_(i),ns_(j),ns_(k),ns_(l),
     &        cs_(:,i),cs_(:,j),cs_(:,k),cs_(:,l),
     &        es_(:,i),es_(:,j),es_(:,k),es_(:,l),
     &        r(:,is_(i)),r(:,is_(j)),r(:,is_(k)),r(:,is_(l)),
     &        ls_(:,i),ls_(:,j),ls_(:,k),ls_(:,l),vt)
            val(ijkl)=vt
          end do
!$omp     end do
!$omp     end parallel
        end subroutine

        subroutine add_s_(i1,cf,ex)
          implicit none
          integer,intent(in)::i1
          real(8),intent(in)::cf(:),ex(:)
          integer::l1(3)
          l1(:)=0
          call basis_add(i1,SIZE(ex),l1,cf,ex)
        end subroutine

        subroutine add_l_(i1,cfs,cfp,ex)
          implicit none
          integer,intent(in)::i1
          real(8),intent(in)::cfs(:),cfp(:),ex(:)
          integer::l1(3),i
          l1(:)=0
          call basis_add(i1,SIZE(ex),l1,cfs,ex)
          do i=1,3
            l1(:)=0
            l1(i)=1
            call basis_add(i1,SIZE(ex),l1,cfp,ex)
          end do
        end subroutine

        subroutine bs_sto3g_(ia,i1)
          ! STO-3G basis set (H-Ne).
          implicit none
          integer,intent(in)::ia,i1
          integer,parameter::MXNA=10
          integer,parameter::NG=3
          real(8),parameter::
     &      zt1(MXNA)=(/
     &        1.24d0,1.69d0,2.69d0,3.68d0,4.68d0,
     &        5.67d0,6.67d0,7.66d0,8.65d0,9.64d0
     &        /),
     &      zt2(MXNA)=(/
     &        0.00d0,0.00d0,0.80d0,1.15d0,1.45d0,
     &        1.72d0,1.95d0,2.25d0,2.55d0,2.88d0
     &        /),
     &      ex1s(NG)=(/
     &        2.22766d0,0.405771d0,0.109818d0
     &        /),
     &      cf1s(NG)=(/
     &        0.154328967295d0,0.535328142282d0,0.444634542185d0
     &        /),
     &      ex2sp(NG)=(/
     &        0.994203d0,0.231031d0,0.0751386d0
     &        /),
     &      cf2s(NG)=(/
     &        -0.099967229187d0,0.399512826089d0,0.700115468880d0
     &        /),
     &      cf2p(NG)=(/
     &        0.155916274999d0,0.607683718598d0,0.391957393099d0
     &        /)
          real(8)::ex(NG)
          ! H - He.
          if (1<=ia.and.ia<=2) then
            ex(:)=ex1s(:)*zt1(ia)**2
            call add_s_(i1,cf1s,ex)
          ! Li - Ne.
          else if (3<=ia.and.ia<=10) then
            ex(:)=ex1s(:)*zt1(ia)**2
            call add_s_(i1,cf1s,ex)
            ex(:)=ex2sp(:)*zt2(ia)**2
            call add_l_(i1,cf2s,cf2p,ex)
          else
            write(*,'("error: out of range: basis.bs_sto3g_()")')
            stop 1
          end if
        end subroutine

        subroutine bs_321g_(ia,i1)
          ! 3-21G basis set (H-Ne).
          implicit none
          integer,intent(in)::ia,i1
          ! H:
          if (ia==1) then
            call add_s_(i1,
     &        (/0.156284978695d0,0.904690876670d0/),
     &        (/5.4471780d0,0.8245472d0/)
     &        )
            call add_s_(i1,
     &        (/1.0d0/),(/0.1831916d0/)
     &        )
          ! He:
          else if (ia==2) then
            call add_s_(i1,
     &        (/0.175229871839d0,0.893482346517d0/),
     &        (/13.6267000d0,1.9993500d0/)
     &        )
            call add_s_(i1,
     &        (/1.0d0/),(/0.3829930d0/)
     &        )
          ! Li:
          else if (ia==3) then
            call add_s_(i1,
     &        (/0.069668663812d0,0.381346349289d0,0.681702624397d0/),
     &        (/36.8382000d0,5.4817200d0,1.1132700d0/)
     &        )
            call add_l_(i1,
     &        (/-0.263126405761d0,1.143387417797d0/),
     &        (/0.161545970837d0,0.915662834700d0/),
     &        (/0.5402050d0,0.1022550d0/)
     &        )
            call add_l_(i1,
     &        (/1.0d0/),(/1.0d0/),(/0.0285645d0/)
     &        )
          ! Be:
          else if (ia==4) then
            call add_s_(i1,
     &        (/0.064426309746d0,0.366096055378d0,0.695934105272d0/),
     &        (/71.8876000d0,10.7289000d0,2.2220500d0/)
     &        )
            call add_l_(i1,
     &        (/-0.421064065879d0,1.224070191516d0/),
     &        (/0.205131923736d0,0.882527671894d0/),
     &        (/1.2954800d0,0.2688810d0/)
     &        )
            call add_l_(i1,
     &        (/1.0d0/),(/1.0d0/),(/0.0773501d0/)
     &        )
          ! B:
          else if (ia==5) then
            call add_s_(i1,
     &        (/0.062960465890d0,0.363303803175d0,0.697254622252d0/),
     &        (/116.4340000d0,17.4314000d0,3.6801600d0/)
     &        )
            call add_l_(i1,
     &        (/-0.368663477319d0,1.199444806451d0/),
     &        (/0.231151902315d0,0.866763633704d0/),
     &        (/2.2818700d0,0.4652480d0/)
     &        )
            call add_l_(i1,
     &         (/1.0d0/),(/1.0d0/),(/0.1243280d0/)
     &        )
          ! C:
          else if (ia==6) then
            call add_s_(i1,
     &        (/0.061766907377d0,0.358794042852d0,0.700713083689d0/),
     &        (/172.2560000d0,25.9109000d0,5.5333500d0/)
     &        )
            call add_l_(i1,
     &        (/-0.395895162119d0,1.215834355681d0/),
     &        (/0.236459946619d0,0.860618805716d0/),
     &        (/3.6649800d0,0.7705450d0/)
     &        )
            call add_l_(i1,
     &         (/1.0d0/),(/1.0d0/),(/0.1958570d0/)
     &        )
          ! N:
          else if (ia==7) then
            call add_s_(i1,
     &        (/0.059865700508d0,0.352955002994d0,0.706513005993d0/),
     &        (/242.7660000d0,36.4851000d0,7.8144900d0/)
     &        )
            call add_l_(i1,
     &        (/-0.413300077430d0,1.224417266851d0/),
     &        (/0.237972016222d0,0.858953058551d0/),
     &        (/5.4252200d0,1.1491500d0/)
     &        )
            call add_l_(i1,
     &         (/1.0d0/),(/1.0d0/),(/0.2832050d0/)
     &        )
          ! O:
          else if (ia==8) then
            call add_s_(i1,
     &        (/0.059239393389d0,0.351499960776d0,0.707657921031d0/),
     &        (/322.0370000d0,48.4308000d0,10.4206000d0/)
     &        )
            call add_l_(i1,
     &        (/-0.404453583190d0,1.221561761397d0/),
     &        (/0.244586106967d0,0.853955373466d0/),
     &        (/7.4029400d0,1.5762000d0/)
     &        )
            call add_l_(i1,
     &         (/1.0d0/),(/1.0d0/),(/0.3736840d0/)
     &        )
          ! F:
          else if (ia==9) then
            call add_s_(i1,
     &        (/0.058548302929d0,0.349308017477d0,0.709632035505d0/),
     &        (/413.8010000d0,62.2446000d0,13.4340000d0/)
     &        )
            call add_l_(i1,
     &        (/-0.407326277682d0,1.223137830990d0/),
     &        (/0.246680003198d0,0.852321011049d0/),
     &        (/9.7775900d0,2.0861700d0/)
     &        )
            call add_l_(i1,
     &         (/1.0d0/),(/1.0d0/),(/0.4823830d0/)
     &        )
          ! Ne:
          else if (ia==10) then
            call add_s_(i1,
     &        (/0.058143030441d0,0.347951182169d0,0.710714372092d0/),
     &        (/515.7240000d0,77.6538000d0,16.8136000d0/)
     &        )
            call add_l_(i1,
     &        (/-0.409922320837d0,1.224310958240d0/),
     &        (/0.247459983573d0,0.851742943461d0/),
     &        (/12.4830000d0,2.6645100d0/)
     &        )
            call add_l_(i1,
     &         (/1.0d0/),(/1.0d0/),(/0.6062500d0/)
     &        )
          else
            write(*,'("error: out of range: basis.bs_sto3g_()")')
            stop 1
          end if
        end subroutine

        integer function ip_(i,j) result(ij)
          ! Canonical index packing for a symmetric matrix.
          implicit none
          integer,intent(in)::i,j
          integer::ii,jj
          ii=max(i,j)
          jj=min(i,j)
          ij=(ii-1)*ii/2+jj
        end function

        integer function ip2_(i,j,k,l) result(ijkl)
          ! Canonical index packing for a symmetric matrix.
          implicit none
          integer,intent(in)::i,j,k,l
          integer::ij,kl
          ij=ip_(i,j)
          kl=ip_(k,l)
          ijkl=ip_(ij,kl)
        end function

        subroutine iup_(ij,i,j)
          ! Unpack the canonical index for a symmetric matrix (1e).
          implicit none
          integer,intent(in)::ij
          integer,intent(out)::i,j
          real(8)::vt
          vt=0.5d0*(-1d0+SQRT(1d0+8*ij))
          i=CEILING(vt)
          j=ij-(i-1)*i/2
        end subroutine

        subroutine iup2_(ijkl,i,j,k,l)
          ! Unpack the canonical index for a symmetric matrix (2e).
          implicit none
          integer,intent(in)::ijkl
          integer,intent(out)::i,j,k,l
          integer::ij,kl
          call iup_(ijkl,ij,kl)
          call iup_(ij,i,j)
          call iup_(kl,k,l)
        end subroutine

      end module

max_stars_repo_name: kyasci/qcref

max_issues_repo_name: kyasci/qcref

max_forks_repo_name: kyasci/qcref


text:
*DECK XSETUA
      SUBROUTINE XSETUA (IUNITA, N)
C***BEGIN PROLOGUE  XSETUA
C***PURPOSE  Set logical unit numbers (up to 5) to which error
C            messages are to be sent.
C***LIBRARY   SLATEC (XERROR)
C***CATEGORY  R3B
C***TYPE      ALL (XSETUA-A)
C***KEYWORDS  ERROR, XERROR
C***AUTHOR  Jones, R. E., (SNLA)
C***DESCRIPTION
C
C     Abstract
C        XSETUA may be called to declare a list of up to five
C        logical units, each of which is to receive a copy of
C        each error message processed by this package.
C        The purpose of XSETUA is to allow simultaneous printing
C        of each error message on, say, a main output file,
C        an interactive terminal, and other files such as graphics
C        communication files.
C
C     Description of Parameters
C      --Input--
C        IUNIT - an array of up to five unit numbers.
C                Normally these numbers should all be different
C                (but duplicates are not prohibited.)
C        N     - the number of unit numbers provided in IUNIT
C                must have 1 .LE. N .LE. 5.
C
C***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
C                 Error-handling Package, SAND82-0800, Sandia
C                 Laboratories, 1982.
C***ROUTINES CALLED  J4SAVE, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   790801  DATE WRITTEN
C   861211  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900510  Change call to XERRWV to XERMSG.  (RWC)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  XSETUA
      DIMENSION IUNITA(5)
      CHARACTER *8 XERN1
C***FIRST EXECUTABLE STATEMENT  XSETUA
C
      IF (N.LT.1 .OR. N.GT.5) THEN
         WRITE (XERN1, '(I8)') N
         CALL XERMSG ('SLATEC', 'XSETUA',
     *      'INVALID NUMBER OF UNITS, N = ' // XERN1, 1, 2)
         RETURN
      ENDIF
C
      DO 10 I=1,N
         INDEX = I+4
         IF (I.EQ.1) INDEX = 3
         JUNK = J4SAVE(INDEX,IUNITA(I),.TRUE.)
   10 CONTINUE
      JUNK = J4SAVE(5,N,.TRUE.)
      RETURN
      END


text:
****************************************************************
*   Virtual matrix element for
*   0 -> q(-p1) + qb(-p4) + a(p2) + a(p3) + lb(p5) + l(p6)
*   l = charged lepton/neutrino
****************************************************************
      subroutine qqb_zaa_v(p,msqv)
      implicit none
      include 'constants.f'
      include 'zprods_decl.f'
      include 'sprods_com.f'
      include 'scheme.f'

      include 'scale.f'
      include 'epinv2.f'
      include 'epinv.f'

      double precision p(mxpart,4),msqv(-nf:nf,-nf:nf)
      double complex a60h(4,16),a6vh(4,16)
      double precision qqb(2),qbq(2),pbdk(mxpart,4)
      integer i,j,k
c-----set scheme
      scheme='dred'
c-----swap momenta
c-----implemented (ala bdk)
c----- 0 -> q(p1) + qb(p4) + gam(p2) + gam(p3) + lb(p5) + l(p6)
c-----mcfm input:
c----- 0 -> q(p1) + qb(p2) + l(p3) + lb(p4) + gam(p5) + gam(p6)
      do i=1,4
         pbdk(1,i)=p(2,i)
         pbdk(4,i)=p(1,i)
         pbdk(2,i)=p(6,i)
         pbdk(3,i)=p(5,i)
         pbdk(5,i)=p(4,i)
         pbdk(6,i)=p(3,i)
      enddo
c-----evaluate spinor products,invariants
      call spinoru(6,pbdk,za,zb)
c-----initialize matelem
      do i=1,2
         qqb(i)=zip
         qbq(i)=zip
      enddo
c-----compute color ordered amplitudes
c-----then compute squared matrix element
c-----qqbar subprocess
      call zaa_a60h(1,2,3,4,5,6,za,zb,a60h)
      call zaa_a6vh(1,2,3,4,5,6,za,zb,a6vh)
      call zaa_m6vsq(2,a60h,a6vh,qqb(2))
      call zaa_m6vsq(1,a60h,a6vh,qqb(1))
c-----qbarq subprocess
      call zaa_a60h(4,2,3,1,5,6,za,zb,a60h)
      call zaa_a6vh(4,2,3,1,5,6,za,zb,a6vh)
      call zaa_m6vsq(2,a60h,a6vh,qbq(2))
      call zaa_m6vsq(1,a60h,a6vh,qbq(1))
c-----initialize msq
      do i=-nf,nf
      do j=-nf,nf
         msqv(i,j)=zip
      enddo
      enddo
c-----fill msq
      msqv(1,-1)=qqb(1)
      msqv(3,-3)=qqb(1)
      msqv(5,-5)=qqb(1)
      msqv(2,-2)=qqb(2)
      msqv(4,-4)=qqb(2)
      msqv(-1,1)=qbq(1)
      msqv(-3,3)=qbq(1)
      msqv(-5,5)=qbq(1)
      msqv(-2,2)=qbq(2)
      msqv(-4,4)=qbq(2)
c-----donehere
      return
      end

***********************************************
* virtual matrix element squared
* qqb_zaa, for initial state flavor qi
***********************************************
      subroutine zaa_m6vsq(qi,a60h,a6vh,msqv)
      implicit none
      include 'constants.f'
      include 'zprods_decl.f'
      include 'sprods_com.f'
      include 'ewcouple.f'
      include 'ewcharge.f'
      include 'qcdcouple.f'
      include 'zcouple.f'
      include 'masses.f'
      include 'ipsgen.f'
      double precision t
      integer j1,j2,j3,j4,j5,j6
      double complex a60h(4,16),m60hA(16),m60hB(16),m60hC(16),m60hD(16)
      double complex a6vh(4,16),m6vhA(16),m6vhB(16),m6vhC(16),m6vhD(16)
      double complex propzQQ,propzLL,propzQL(3:4)
      double precision m6sqvh(16),msqv
      double precision m6sqvhAA(16),m6sqvhBB(16),m6sqvhCC(16),
     .m6sqvhDD(16),m6sqvhAB(16),m6sqvhAC(16),m6sqvhAD(16),m6sqvhBC(16),
     .m6sqvhBD(16),m6sqvhCD(16)
      double precision iwdth,iph,izz,qq1
      integer qi,i,j,k
c-----Z propagator
      iwdth = 1D0
c-----QQ propagator
      propzQQ=s(5,6)/Dcmplx(s(5,6)-zmass**2,iwdth*zwidth*zmass)
c-----LL propagator
      propzLL=s(1,4)/Dcmplx(s(1,4)-zmass**2,iwdth*zwidth*zmass)
c-----QL propagator (2,3)
      propzQL(3)=t(1,4,2)/Dcmplx(t(1,4,2)-zmass**2,iwdth*zwidth*zmass)
c-----QL propagator (3,2)
      propzQL(4)=t(1,4,3)/Dcmplx(t(1,4,3)-zmass**2,iwdth*zwidth*zmass)
c-----dress the amplitude with couplings etc.
c-----iph,izz -> switches for photon/z contribution
      iph=one
      izz=one
c-----qq1 -> choose electron/neutrino contributions
      qq1=dabs(q1)
c-----LO amplitude
c-----M(+-xx-+)
      do i=1,4
         if (ipsgen .ne. 2) then
         m60hA(i)=four*esq**2*
     .   (+Q(qi)**2*(q1*Q(qi)+r1*r(qi)*propzQQ)*a60h(1,i))
         endif
         if (ipsgen .ne. 3) then
         m60hB(i)=four*esq**2*
     .   (+qq1*(q1*Q(qi)*iph+izz*r1*r(qi)*propzLL)*a60h(2,i))
         endif
         if ((ipsgen .eq. 2) .or. (ipsgen .eq. 3)) then
         m60hC(i)=four*esq**2*
     .   (-qq1*Q(qi)*(q1*Q(qi)+r1*r(qi)*propzQL(3))*a60h(3,i))
         endif
         if ((ipsgen .eq. 3) .or. (ipsgen .eq. 4)) then
         m60hD(i)=four*esq**2*
     .   (-qq1*Q(qi)*(q1*Q(qi)+r1*r(qi)*propzQL(4))*a60h(4,i))
         endif
      enddo
c-----M(+-xx+-)
      do i=5,8
         if (ipsgen .ne. 2) then
         m60hA(i)=four*esq**2*
     .   (+Q(qi)**2*(q1*Q(qi)+l1*r(qi)*propzQQ)*a60h(1,i))
         endif
         if (ipsgen .ne. 3) then
         m60hB(i)=four*esq**2*
     .   (+qq1*(q1*Q(qi)+l1*r(qi)*propzLL)*a60h(2,i))
         endif
         if ((ipsgen .eq. 2) .or. (ipsgen .eq. 3)) then
         m60hC(i)=four*esq**2*
     .   (+qq1*Q(qi)*(q1*Q(qi)+l1*r(qi)*propzQL(3))*a60h(3,i))
         endif
         if ((ipsgen .eq. 3) .or. (ipsgen .eq. 4)) then
         m60hD(i)=four*esq**2*
     .   (+qq1*Q(qi)*(q1*Q(qi)+l1*r(qi)*propzQL(4))*a60h(4,i))
         endif
      enddo
c-----(-+xx+-)
      do i=9,12
         if (ipsgen .ne. 2) then
         m60hA(i)=four*esq**2*
     .   (+Q(qi)**2*(q1*Q(qi)+l1*l(qi)*propzQQ)*a60h(1,i))
         endif
         if (ipsgen .ne. 3) then
         m60hB(i)=four*esq**2*
     .   (+qq1*(q1*Q(qi)+l1*l(qi)*propzLL)*a60h(2,i))
         endif
         if ((ipsgen .eq. 2) .or. (ipsgen .eq. 3)) then
         m60hC(i)=four*esq**2*
     .   (-qq1*Q(qi)*(q1*Q(qi)+l1*l(qi)*propzQL(3))*a60h(3,i))
         endif
         if ((ipsgen .eq. 3) .or. (ipsgen .eq. 4)) then
         m60hD(i)=four*esq**2*
     .   (-qq1*Q(qi)*(q1*Q(qi)+l1*l(qi)*propzQL(4))*a60h(4,i))
         endif
      enddo
c-----(-+xx-+)
      do i=13,16
         if (ipsgen .ne. 2) then
         m60hA(i)=four*esq**2*
     .   (+Q(qi)**2*(q1*Q(qi)+r1*l(qi)*propzQQ)*a60h(1,i))
         endif
         if (ipsgen .ne. 3) then
         m60hB(i)=four*esq**2*
     .   (+qq1*(q1*Q(qi)+r1*l(qi)*propzLL)*a60h(2,i))
         endif
         if ((ipsgen .eq. 2) .or. (ipsgen .eq. 3)) then
         m60hC(i)=four*esq**2*
     .   (+qq1*Q(qi)*(q1*Q(qi)+r1*l(qi)*propzQL(3))*a60h(3,i))
         endif
         if ((ipsgen .eq. 3) .or. (ipsgen .eq. 4)) then
         m60hD(i)=four*esq**2*
     .   (+qq1*Q(qi)*(q1*Q(qi)+r1*l(qi)*propzQL(4))*a60h(4,i))
         endif
      enddo
c-----Virtual amplitude
c-----M(+-xx-+)
      do i=1,4
         if (ipsgen .ne. 2) then
         m6vhA(i)=four*esq**2*
     .   (+Q(qi)**2*(q1*Q(qi)+r1*r(qi)*propzQQ)*a6vh(1,i))
         endif
         if (ipsgen .ne. 3) then
         m6vhB(i)=four*esq**2*
     .   (+qq1*(q1*Q(qi)*iph+izz*r1*r(qi)*propzLL)*a6vh(2,i))
         endif
         if ((ipsgen .eq. 2) .or. (ipsgen .eq. 3)) then
         m6vhC(i)=four*esq**2*
     .   (-qq1*Q(qi)*(q1*Q(qi)+r1*r(qi)*propzQL(3))*a6vh(3,i))
         endif
         if ((ipsgen .eq. 3) .or. (ipsgen .eq. 4)) then
         m6vhD(i)=four*esq**2*
     .   (-qq1*Q(qi)*(q1*Q(qi)+r1*r(qi)*propzQL(4))*a6vh(4,i))
         endif
      enddo
c-----M(+-xx+-)
      do i=5,8
         if (ipsgen .ne. 2) then
         m6vhA(i)=four*esq**2*
     .   (+Q(qi)**2*(q1*Q(qi)+l1*r(qi)*propzQQ)*a6vh(1,i))
         endif
         if (ipsgen .ne. 3) then
         m6vhB(i)=four*esq**2*
     .   (+qq1*(q1*Q(qi)+l1*r(qi)*propzLL)*a6vh(2,i))
          endif
         if ((ipsgen .eq. 2) .or. (ipsgen .eq. 3)) then
        m6vhC(i)=four*esq**2*
     .   (+qq1*Q(qi)*(q1*Q(qi)+l1*r(qi)*propzQL(3))*a6vh(3,i))
         endif
         if ((ipsgen .eq. 3) .or. (ipsgen .eq. 4)) then
         m6vhD(i)=four*esq**2*
     .   (+qq1*Q(qi)*(q1*Q(qi)+l1*r(qi)*propzQL(4))*a6vh(4,i))
         endif
      enddo
c-----(-+xx+-)
      do i=9,12
         if (ipsgen .ne. 2) then
         m6vhA(i)=four*esq**2*
     .   (+Q(qi)**2*(q1*Q(qi)+l1*l(qi)*propzQQ)*a6vh(1,i))
         endif
         if (ipsgen .ne. 3) then
         m6vhB(i)=four*esq**2*
     .   (+qq1*(q1*Q(qi)+l1*l(qi)*propzLL)*a6vh(2,i))
         endif
         if ((ipsgen .eq. 2) .or. (ipsgen .eq. 3)) then
         m6vhC(i)=four*esq**2*
     .   (-qq1*Q(qi)*(q1*Q(qi)+l1*l(qi)*propzQL(3))*a6vh(3,i))
         endif
         if ((ipsgen .eq. 3) .or. (ipsgen .eq. 4)) then
         m6vhD(i)=four*esq**2*
     .   (-qq1*Q(qi)*(q1*Q(qi)+l1*l(qi)*propzQL(4))*a6vh(4,i))
         endif
      enddo
c-----(-+xx-+)
      do i=13,16
         if (ipsgen .ne. 2) then
         m6vhA(i)=four*esq**2*
     .   (+Q(qi)**2*(q1*Q(qi)+r1*l(qi)*propzQQ)*a6vh(1,i))
         endif
         if (ipsgen .ne. 3) then
         m6vhB(i)=four*esq**2*
     .   (+qq1*(q1*Q(qi)+r1*l(qi)*propzLL)*a6vh(2,i))
         endif
         if ((ipsgen .eq. 2) .or. (ipsgen .eq. 3)) then
         m6vhC(i)=four*esq**2*
     .   (+qq1*Q(qi)*(q1*Q(qi)+r1*l(qi)*propzQL(3))*a6vh(3,i))
         endif
         if ((ipsgen .eq. 3) .or. (ipsgen .eq. 4)) then
         m6vhD(i)=four*esq**2*
     .   (+qq1*Q(qi)*(q1*Q(qi)+r1*l(qi)*propzQL(4))*a6vh(4,i))
         endif
      enddo
c-----square the helicity amplitudes
      do i=1,16
         if (ipsgen .eq. 1) then
         m6sqvhAA(i)=two*dble(dconjg(m60hA(i))*m6vhA(i))
         m6sqvhAB(i)=two*dble( dconjg(m60hA(i))*m6vhB(i)
     .                        +dconjg(m60hB(i))*m6vhA(i) )
         endif
         if (ipsgen .eq. 2) then
         m6sqvhBB(i)=two*dble(dconjg(m60hB(i))*m6vhB(i))
         m6sqvhBC(i)=two*dble( dconjg(m60hB(i))*m6vhC(i)
     .                        +dconjg(m60hC(i))*m6vhB(i) )
         endif
         if (ipsgen .eq. 3) then
         m6sqvhCC(i)=two*dble(dconjg(m60hC(i))*m6vhC(i))
         m6sqvhAC(i)=two*dble( dconjg(m60hA(i))*m6vhC(i)
     .                        +dconjg(m60hC(i))*m6vhA(i) )
         m6sqvhCD(i)=two*dble( dconjg(m60hC(i))*m6vhD(i)
     .                        +dconjg(m60hD(i))*m6vhC(i) )
         endif
         if (ipsgen .eq. 4) then
         m6sqvhDD(i)=two*dble(dconjg(m60hD(i))*m6vhD(i))
         m6sqvhAD(i)=two*dble( dconjg(m60hA(i))*m6vhD(i)
     .                        +dconjg(m60hD(i))*m6vhA(i) )
         m6sqvhBD(i)=two*dble( dconjg(m60hB(i))*m6vhD(i)
     .                        +dconjg(m60hD(i))*m6vhB(i) )
         endif
      enddo
      do i=1,16
c--------commented pieces below only for cross checking with ols PSgen
c         m6sqvh(i)= m6sqvhAA(i)+m6sqvhBB(i)+m6sqvhCC(i)+m6sqvhDD(i)
c     .             +m6sqvhAB(i)+m6sqvhAC(i)+m6sqvhAD(i)+m6sqvhBC(i)
c     .             +m6sqvhBD(i)+m6sqvhCD(i)
         if (ipsgen.eq.1) then
             m6sqvh(i)= m6sqvhAA(i)+m6sqvhAB(i) !AA+AB
         elseif (ipsgen.eq.2) then
             m6sqvh(i)= m6sqvhBB(i)+m6sqvhBC(i) !BB+BC
         elseif (ipsgen.eq.3) then
             m6sqvh(i)= m6sqvhCC(i)+m6sqvhAC(i)+m6sqvhCD(i) !CC+AC+CD
         elseif (ipsgen.eq.4) then
             m6sqvh(i)= m6sqvhDD(i)+m6sqvhAD(i)+m6sqvhBD(i) !DD+AD+BD
         endif
      enddo
c-----sum them up
      msqv = zip
      do i=1,16
         msqv=msqv+m6sqvh(i)
      enddo
c-----multiply by color factors,identical particle factors
c-----average over colors,spins
      msqv=msqv*half*aveqq*(Nc-one/Nc)*three*gsq/16D0/pi**2
c-----donehere
      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
( Simple tilemap collision )
depend sample/platformer/lib/array2d.f

( what sides the object collided )
0 value lwall?
0 value rwall?
0 value floor?
0 value ceiling?

[undefined] tstep@ [if] : tstep@ 16 16 ; [then]

extend: _actor
    var mbw  var mbh   \ object collision box width,height
;class

defer on-tilemap-collide  ' drop is on-tilemap-collide  ( tilecell - )
defer tileprops@   :noname drop 0 ; is tileprops@  ( tilecell - bitmask )

#1
    bit BIT_CEL
    bit BIT_FLR
    bit BIT_WLT
    bit BIT_WRT
value tile-bits

define collisioning

    0 value map
    : map@  ( col row - tile )  map loc @ ;

    : cel? BIT_CEL and ; \ ' ceiling '
    : flr? BIT_FLR and ; \ ' floor '
    : wlt? BIT_WLT and ; \ ' wall left '
    : wrt? BIT_WRT and ; \ ' wall right '

    : vector   create 0 , here 0 , constant ;
    vector nx ny

    : gap  ( - n ) tstep@ drop ;  \ just square tiles supported for now

    : px x @ ;
    : py y @ ;

    variable t
    : xy>cr  ( x y tilesize - ) dup  2/  2pfloor ;
    : pt  gap xy>cr  map@ dup t !  tileprops@ ;          \ point

    ( increment coordinates )
    : ve+  swap gap +  mbw @ #1 - px +  min  swap ;
    : he+  gap +  mbh @ #1 - ny @ +  min ;

    : +vy ny +! ny @ py - vy ! ;
    : +vx nx +! nx @ px - vx ! ;

    ( push up/down )
    : pu ( xy ) nip gap mod negate +vy  true to floor?  t @ on-tilemap-collide  ;
    : pd ( xy ) nip gap mod negate gap + +vy  true to ceiling?  t @ on-tilemap-collide ;

    ( check up/down )
    : cu mbw @ gap / 2 + for 2dup pt cel? if pd unloop exit then ve+ loop 2drop ;
    : cd mbw @ gap / 2 + for 2dup pt flr? if pu unloop exit then ve+ loop 2drop ;

    ( push left/right )
    : pl ( xy ) drop gap mod negate ( -1.0 + ) +vx  true to rwall?  t @ on-tilemap-collide ;
    : pr ( xy ) drop gap mod negate gap + +vx  true to lwall?  t @ on-tilemap-collide ;

    ( check left/right )
    : cl mbh @ gap /  2 + for 2dup pt wrt? if pr unloop exit then he+ loop 2drop ;
    : crt mbh @ gap / 2 + for 2dup pt wlt? if pl unloop exit then he+ loop 2drop ;

    : ud vy @ -exit vy @ 0 < if px ny @ cu exit then px ny @ mbh @ + cd ;
    : lr vx @ -exit vx @ 0 < if nx 2@ cl exit then   nx @ mbw @ + ny @ crt ;

    : init   px py  vx 2@  2+  nx 2!  0 to lwall? 0 to rwall? 0 to floor? 0 to ceiling? ;

only forth definitions fixed
also collisioning

: collide-tilemap ( array2d - ) to map  init ud lr ;

only forth definitions

max_stars_repo_name: RogerLevy/RamenEngine

max_issues_repo_name: Alim-Oezdemir/RamenEngine

max_forks_repo_name: Alim-Oezdemir/RamenEngine


text:
        subroutine msgmerge(a,b,c)
        save
        character *1 a(1),b(1),c(1),ast
        data ast/'*'/
c
        do 1200 i=1,1000
c
        if(a(i) .eq. ast) goto 1400
        c(i)=a(i)
        iadd=i
 1200 continue
c
 1400 continue
c
        do 1800 i=1,1000
c
        c(iadd+i)=b(i)
        if(b(i) .eq. ast) return
 1800 continue
        return
        end
c
c
c
c
c

max_stars_repo_name: ZhekehZ/catboost

max_issues_repo_name: birichie/catboost

max_forks_repo_name: birichie/catboost


text:
PROGRAM rwfnplot
   IMPLICIT NONE

   INTEGER, PARAMETER:: NPTS0 = 5000

   DOUBLE PRECISION pg(NPTS0), qg(NPTS0), rg(NPTS0), pgg(NPTS0)
   DOUBLE PRECISION rg2(NPTS0)
   DOUBLE PRECISION energy, a0
   CHARACTER title*6, orbl*4, nnstr*2, new*3, xa*3, el*5
   CHARACTER leg*50, format_string*50, name*100
   DIMENSION el(30)
   INTEGER np, lp, jp, nn, laky, ll, jj, npts, j, i, nwf

   write (*, *) 'RWFNPLOT'
   write (*, *) 'Program to generate Matlab/GNU Octave and'
   write (*, *) 'Xmgrace files that plot radial orbitals'
   write (*, *) 'Input file:  name.w'
   write (*, *) 'Output files: octave_name.m, xmgrace_name.agr'
   write (*, *)
   write (*, *) 'To plot orbital: press enter'
   write (*, *) 'To remove orbital: type "d" or "D" and press enter'
   write (*, *)
   write (*, *) '                             Jorgen Ekman Jun 2015'
   write (*, *)

   WRITE (*, *) 'Name of state:'
   READ (*, *) name
   OPEN (3, FILE=trim(name)//'.w', STATUS='OLD', FORM='UNFORMATTED')
   OPEN (4, FILE='octave_'//trim(name)//'.m', STATUS='UNKNOWN')
   OPEN (8, FILE='xmgrace_'//trim(name)//'.agr', STATUS='UNKNOWN')

   write (*, *)
   write (*, *) 'To have r on x-axis: type "y" otherwise "n" for sqrt(r)'
   read (*, *) xa
   write (*, *)
   i = 1

   write (4, *) '%     r                   P(r)                   Q(r)'

   ! Xmgrace stuff
   if (xa .eq. 'y  ') then
      write (8, *) '@    xaxis  label "r"'
   else
      write (8, *) '@    xaxis  label "sqrt(r)"'
   end if
   write (8, *) '@    yaxis  label "P(r)"'
   write (8, *) '#     r                   P(r)                   Q(r)'

   write (4, *) 'clf'
   READ (3) title
   IF (title .NE. 'G92RWF') THEN   ! Extra safety
      PRINT *, 'title = ', title, 'does not match G92RWF'
      STOP
   ENDIF

   DO
      READ (3, END=20) nn, laky, energy, npts
      write (nnstr, '(I2)') nn
      IF (laky .GT. 0) THEN
         ll = laky
         jj = -1
      ELSEIF (laky .LE. -1) THEN
         ll = -laky - 1
         jj = 1
      ELSE
         WRITE (*, *) 'Unexpected case in reading mcdf.w'
         STOP
      ENDIF

      if (jj .eq. 1) then
         if (ll .eq. 0) orbl = nnstr//'s '
         if (ll .eq. 1) orbl = nnstr//'p '
         if (ll .eq. 2) orbl = nnstr//'d '
         if (ll .eq. 3) orbl = nnstr//'f '
         if (ll .eq. 4) orbl = nnstr//'g '
         if (ll .eq. 5) orbl = nnstr//'h '
         if (ll .eq. 6) orbl = nnstr//'i '
         if (ll .eq. 7) orbl = nnstr//'k '
         if (ll .eq. 8) orbl = nnstr//'l '
         if (ll .eq. 9) orbl = nnstr//'m '
      else
         if (ll .eq. 1) orbl = nnstr//'p-'
         if (ll .eq. 2) orbl = nnstr//'d-'
         if (ll .eq. 3) orbl = nnstr//'f-'
         if (ll .eq. 4) orbl = nnstr//'g-'
         if (ll .eq. 5) orbl = nnstr//'h-'
         if (ll .eq. 6) orbl = nnstr//'i-'
         if (ll .eq. 7) orbl = nnstr//'k-'
         if (ll .eq. 8) orbl = nnstr//'l-'
         if (ll .eq. 9) orbl = nnstr//'m-'
      end if

      IF (npts .GT. NPTS0) THEN
         WRITE (*, *) 'df2hf: npts .GT. NPTS0'
         STOP
      ENDIF

      READ (3) a0, (pg(j), j=1, npts), (qg(j), j=1, npts)
      READ (3) (rg(j), j=1, npts)

      WRITE (*, '(2x,A,A,$)') orbl, ' = '
      READ (*, '(A)') new
      IF (NEW .NE. 'd  ' .AND. NEW .NE. 'D  ') THEN
         el(i) = orbl

         ! Xmgrace stuff
         if (i - 1 .lt. 10) then
            format_string = "(A,I1,3A)"
         else
            format_string = "(A,I2,3A)"
         end if
         write (leg, format_string) &
            '@    s', i - 1, '  legend  "', el(i), '"'
         write (8, *) trim(leg)
         write (8, *) '# ', el(i)

         i = i + 1
         write (4, *) 'P = ['
         DO j = 1, npts
            rg2(j) = sqrt(rg(j))
            if (xa .eq. 'y  ') then
               if (abs(pg(j)) .gt. 0.0005 .OR. j .EQ. 1) &
                  WRITE (4, '(3D20.10)') rg(j), pg(j), qg(j)
               WRITE (8, '(3D20.10)') rg(j), pg(j), qg(j)
            else
               if (abs(pg(j)) .gt. 0.0005 .OR. j .EQ. 1) &
                  WRITE (4, '(3D20.10)') rg2(j), pg(j), qg(j)
               WRITE (8, '(3D20.10)') rg2(j), pg(j), qg(j)
            end if
         ENDDO
         write (8, *)
         write (4, *) '];'
         write (4, *) 'plot(P(:,1), P(:,2))'
         write (4, *) 'hold all'
      ENDIF
   ENDDO

20 CONTINUE

   nwf = i - 1

   if (xa .eq. 'y  ') then
      write (4, *) 'xlabel (''r'')'
   else
      write (4, *) 'xlabel (''sqrt(r)'')'
   end if

   write (4, *) 'ylabel (''P(r)'')'
   write (4, *) 'grid on'
   write (4, 101, advance='no') 'legend('
   DO i = 1, nwf
      IF (i .LT. nwf) THEN
         write (4, 102, advance='no') '''', el(i), ''','
      ELSE
         write (4, 102, advance='no') '''', el(i), ''')'
      END IF
   ENDDO
101 format(a)
102 format(3a)
   PRINT *, ' Finished.'
   STOP
END

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: sylas/grasp-continuum

max_forks_repo_name: sylas/grasp-continuum


! PR debug/71642
! { dg-do compile }
! { dg-options "-g" }

MODULE gauss_colloc
  INTEGER, PARAMETER :: dp=8
CONTAINS
SUBROUTINE collocGauss(h,h_inv,grid,poly,alphai,posi,max_r2,&
        periodic,gdim,local_bounds,local_shift,poly_shift,scale,lgrid,error)
    REAL(dp), DIMENSION(0:, 0:, 0:), &
      INTENT(inout)                          :: grid
    INTEGER, INTENT(inout), OPTIONAL :: lgrid
    CONTAINS
    SUBROUTINE kloop6
    IF (kJump/=1 .AND. (ikstart+kmax-kstart>=ndim(2)+l_shift(2) .OR.&
        ikstart2+kmin-kstart2<=l_ub(2)-ndim(2))) THEN
        DO
            DO k=kstart2,kend2,-1
                IF ( PRESENT ( lgrid ) ) THEN
                  grid(ik,ij,ii) = grid(ik,ij,ii) + p_v*res_k
                END IF
            END DO
        END DO
    END IF
    END SUBROUTINE
END SUBROUTINE
END MODULE gauss_colloc


text:
c>> 2006-03-11 DRDSPGE  Krogh  Initial code.
c--D replaces "?": DR?SPGE, ?SPGE
c Demo driver for DRSPGE
      integer N, LA, LIA
      parameter (N=6, LA=300, LIA=LA+5*N)
      integer IA(LA+5*N), ISPEC(7), I, J, K
      double precision A(LA),OPT(3),SIZE,B(N),BB(N),X(N),ERRX,ERRRES
      data IA(1),A(1), IA(2),A(2), IA(3),A(3) /2,0.D0, 1,-2.D0, 4,1.D0/
      data IA(4),A(4), IA(5),A(5) /1,0.D0, 3,-2.E0/
      data IA(6),A(6), IA(7),A(7), IA(8),A(8) /2,0.D0, 2,-2.D0, 5,-2.D0/
      data IA( 9),A( 9), IA(10),A(10), IA(11),A(11), IA(12),A(12)  /
     1  3,0.D0, 1,-1.D0, 5,-1.D0, 6,2.D0 /
      data IA(13),A(13), IA(14),A(14), IA(15),A(15), IA(16),A(16),
     1  IA(17),A(17), IA(18),A(18) / 5,0.D0, 1,-1.D0, 2,-2.D0, 3,2.D0,
     2  4,-1.D0, 6,-2.D0 /
      data IA(19),A(19), IA(20),A(20), IA(21),A(21), IA(22),A(22) /
     1  3,0.D0,  1,-2.D0, 3,2.D0, 6,1.D0 /
      data X / -2.0D0, -1.0D0, -2.0D0, -1.0D0, -1.0D0, -1.0D0 /
      data B /  8.0D0,  6.0D0, -2.0D0, -1.0D0,  5.0D0, -1.0D0 /
      data ISPEC / 0, LIA, LA, 0, 3, 4, 0 /

      do 100 I = 1, N
         BB(I) = B(I)
 100  continue
c              Call DSPGE to get the solution
      call DSPGE(N, ISPEC, IA, A, B, OPT)
c            Get maximal errors in the residual and the solution.
      ERRRES = 0.D0
      SIZE = 0.D0
      K = 1
      do 280 J = 1, N
        do 270 I = K+1, K + IA(K)
          SIZE = max(SIZE, abs(A(I) * X(J)))
          BB(IA(I)) = BB(IA(I)) - A(I) * B(J)
 270    continue
        K = I
 280  continue
      do 300 J = 1, N
        ERRRES = max(ERRRES, abs(BB(J)))
 300  continue
      ERRRES = ERRRES / SIZE
      ERRX = 0.D0
      do 320 I = 1, N
        ERRX = max(ERRX, abs(X(I) - B(I)))
 320  continue

      I = int(log10(abs(OPT(2))))
      if (I .gt. 0) I = I + 1
      OPT(2) = OPT(2)*10.D0**(-I)
      I = I + nint(OPT(3))
      print '('' Problem     N   RESERR     XERR    Unused    Used'',
     1  ''    RCOND      DET   x10^?'')'
      print '(''DEMOTEST     6'',1P,2E10.2,2I8,E10.2,0PF9.5,I6)',
     1  ERRRES,ERRX,ISPEC(4),ISPEC(3)-ISPEC(4),OPT(1),OPT(2),I
      stop
      end


text:
C$Procedure                  CRTPTR (Create pointer)

      CHARACTER*(*) FUNCTION CRTPTR (BASE, INDEX, PNTER)

C$ Abstract
C
C     Returns the symbol 'BASE~INDEX~PNTER'.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C     None.
C
C$ Keywords
C
C     None.
C
C$ Declarations

      CHARACTER*(*)         BASE
      INTEGER               INDEX
      CHARACTER*(*)         PNTER

      INTEGER               SYMLEN
      PARAMETER           ( SYMLEN = 1024 )

C$ Brief_I/O
C
C     Variable  I/O  Description
C     --------  ---  --------------------------------------------------
C      BASE,
C      INDEX,
C      PNTER     I    Components of the symbol 'BASE~INDEX~PNTER'.
C
C      SYMLEN    P    Maximum length of the symbol.
C
C$ Detailed_Input
C
C     BASE,
C     INDEX,
C     PNTER      are components of the symbol 'BASE~INDEX~PNTER'.
C
C$ Detailed_Output
C
C     CRTPTR     is the symbol 'BASE~INDEX~PNTER'.
C
C$ Parameters
C
C     SYMLEN     is the maximum length of the symbol 'BASE~INDEX~PNTER'.
C
C$ Exceptions
C
C     1) If the length of the symbol 'BASE~INDEX~PNTER' exceeds SYMLEN,
C        the error SPICE(BUFFERTOOSMALL) is signalled.
C
C     2) If the length of the symbol 'BASE~INDEX~PNTER' exceeds
C        LEN(CRTPTR), the error SPICE(DIMENSIONTOOSMALL) is signalled.
C
C$ Files
C
C     None.
C
C$ Particulars
C
C     This routine creates a symbol that may be used to look up
C     nodes in the symbol table created by CPARSE.
C
C$ Examples
C
C     None.
C
C$ Restrictions
C
C     None.
C
C$ Literature_References
C
C     None.
C
C$ Author_and_Institution
C
C     M.J. Spencer   (JPL)
C
C$ Version
C
C-    Version 1.1.1, 13-JAN-2007, (EDW)
C
C        Corrected typo in the previous version string;
C        from:
C
C           09-DEC-203
C
C        to
C
C           09-DEC-2003
C
C-    Version 1.1.0, 09-DEC-2003, (EDW)
C
C       Set the SYMLEN value to 1024 to match the same
C       value in niospk.
C
C-    Beta Version 1.0.0, 11-AUG-1992 (MJS)
C
C-&

C
C     SPICELIB functions
C
      INTEGER               RTRIM
      LOGICAL               RETURN

C
C     Local variables
C
      CHARACTER*(10)        CNUM
      CHARACTER*(SYMLEN)    SYM

      INTEGER               BLEN
      INTEGER               CLEN
      INTEGER               PLEN
      INTEGER               TOTAL

C
C     This routine will use discovery check-in.
C
      IF (RETURN()) THEN
         RETURN
      END IF

C
C     Compute the lengths of the strings involved.
C
      CALL INTSTR (INDEX, CNUM)

      CLEN  = RTRIM (CNUM)
      BLEN  = RTRIM (BASE)
      PLEN  = RTRIM (PNTER)
      TOTAL = CLEN+BLEN+PLEN+2

C
C     TOTAL must be SYMLEN characters, or fewer.
C
      IF (TOTAL .GT. SYMLEN) THEN
         CALL CHKIN  ('CRTPTR'                                     )
         CALL SETMSG ('Symbol exceeds # characters. Increase the '  //
     .                'value of SYMLEN.'                           )
         CALL ERRINT ('#', SYMLEN                                  )
         CALL SIGERR ('SPICE(BUFFERTOOSMALL)'                      )
         CALL CHKOUT ('CRTPTR'                                     )
         RETURN
      END IF

C
C     And TOTAL must be LEN(CRTPTR) characters, or fewer.
C
      IF (TOTAL .GT. LEN(CRTPTR)) THEN
         CALL CHKIN  ('CRTPTR'                                 )
         CALL SETMSG ('Symbol exceeds the dimension of CRTPTR.')
         CALL SIGERR ('SPICE(DIMENSIONTOOSMALL)'               )
         CALL CHKOUT ('CRTPTR'                                 )
         RETURN
      END IF

C
C     Form the symbol 'BASE~INDEX~PNTER'.
C
      SYM = ' '
      SYM = BASE(1:BLEN) //'~'// CNUM(1:CLEN) //'~'// PNTER(1:PLEN)

      CRTPTR = SYM

      RETURN
      END

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
\
\ Math
\

:0> 0< 0= ;
:0<> 0= 0= ;

:1+ 1 + ;
:1- 1 - ;
:2* dup + ;

:*/ >r * r> /mod drop ;
:/ /mod drop ;
:mod /mod nip ;

:< - 0< ;
:> - 0> ;

\
\ Stack operators
\

:2drop  drop drop ;
:swap over >r nip r> ;
:rot >r swap r> swap ;

\
\ Memory access
\

:@ >a a@+ ;
:! >a a!+ ;

:2@ >a a@+ a@+ ;
:2! >a a!+ a!+ ;

:2, over , nip , ;

\
\ Combinators
\
\ largely from https://gist.github.com/crcx/8060687#file-combinators-forth

:dip swap >r call r> ;
:sip over >r call r> ;
:bi [ sip ] dip call ;
:bi* [ dip ] dip call ;
:bi@ dup bi* ;

\ Call the quotes f, g, and h with x, leaving the result of each on the stack
\
\ ( x f g h -- )
:tri swap >r [ sip ] dip sip r> call swap ;

\
\ Variables
\

\ Binds a name to a function that pushes the given value to the stack
\
\ ( value 'name -- )
:constant here &lit , rot , &exit , swap bind ;

0 'false constant
1 'true constant

\
\ Control flow
\

:?swap  ?exit swap ;

\ Calls false or true execution token depending on the flag
\
\ ( flag false true -- )
:choose  rot ?swap nip call ;

\ Like choose, but does a goto instead of a call
\
\ ( flag false true -- )
:whereto rot ?swap nip goto ;

\ Call the token if the flag is true
\
\ ( flag token -- )
:if  [ ] swap choose ;

\ Opposite of -if
\
\ ( flag token -- )
:-if [ ] choose ;

\ Decorator to execute the rest of a function n times. The top item of the
\ stack counts to zero
\
\ ( n -- n' )
:times>  1 -  r@  over [ drop ] [ r> drop call &times> goto ] whereto ;

\
\ Types
\

:number? type 'number = ;
:string? type 'string = ;

max_stars_repo_name: aotmr/evef

max_issues_repo_name: aotmr/elen

max_forks_repo_name: aotmr/elen