text:
c     ----------------------------------------------------------------------
c     Subroutines for Frontal Detection Tool
c     ----------------------------------------------------------------------

c     -----------------------------------------------
c     label connected fronts
c     -----------------------------------------------

      subroutine clustering(outar,ntot,inar,nx,ny,grid)

c     Given the input array <inar> with 0/1 entries (for example the
c     NCOFF array), the individual events are counted and each event
c     gets an individual 'stamp' which is specific to this event. An
c     event is characterised by 1-entries in <inar> which are connected.
c     <nx>,<ny> give the dimension of the arrays, <ntot> the total
c     number of events i.e. objects.

      implicit none

c     Declaration of subroutine parameters
      integer, intent(in)	:: nx,ny
      integer, intent(in)	:: grid
      integer, intent(inout) 	:: outar(nx,ny)
      integer, intent(inout)  	:: ntot
      real, intent(in) 		:: inar(nx,ny)

c     Auxiliary variables
      integer			:: i,j

c     Copy inar to outar
      outar=nint(inar)

c     -------- Do the clustering -----------
c     total number of cluster
      ntot=1

c     Find an element which does not yet belong to a cluster
       do j=1,ny
        do i=1,nx
          if (outar(i,j).eq.1) then
            ntot=ntot+1
            outar(i,j)=ntot
            call connect(outar,ntot,i,j,nx,ny,grid)
          endif
        enddo
      enddo

c     Correct the output array (so far, the index is 1 unit to large)
      where (outar.gt.1) outar=outar-1
      ntot=ntot-1

c     -----------------------------------------------
      end subroutine clustering
c     -----------------------------------------------

c     -----------------------------------------------
c     mark connected gridpoints with fronts
c     -----------------------------------------------

      subroutine connect (outar,label,i,j,nx,ny,grid)

      implicit none

c     Mark all elements in outar which are connected to element i,j

c     Declaration of subroutine parameters
      integer, intent(in)	:: nx,ny
      integer, intent(in)	:: grid
      integer, intent(in)	:: label
      integer, intent(in) 	:: i,j
      integer, intent(inout)	:: outar(nx,ny)

c     Auxiliary variables
      integer  :: il,ir,ju,jd,im,jm
      integer  :: k
      integer  :: stack
      integer  :: indx(100000),indy(100000)
      integer  :: workar(0:nx+1,0:ny+1)

c     Initialize working array
      workar=0
      workar(1:nx,1:ny)=outar

c     Push the first element on the stack
      stack=1
      indx(stack)=i
      indy(stack)=j

c     Define the indices of the neighboring elements
      wstack: do

        il=indx(stack)-1
        ir=indx(stack)+1
        ju=indy(stack)+1
        jd=indy(stack)-1
        im=indx(stack)
        jm=indy(stack)
        stack=stack-1

c      Grid handling
c       periodic closed, cyclic point
        if ( grid.eq.2 ) then
           if ( il.eq.0      ) il = nx-1
           if ( ir.eq.(nx+1) ) ir = 2
c       periodic no cyclic point (not closed)
        elseif ( grid.eq.1 ) then
           if ( il.eq.0      ) il = nx
           if ( ir.eq.(nx+1) ) ir = 1
c       not periodic
        else
           if ( il.eq.0      ) il = 1
           if ( ir.eq.(nx+1) ) ir = nx
        endif

c       Mark all connected elements (build up the stack)
        if (workar(il,jm).eq.1) then
          workar(il,jm)=label
          stack=stack+1
          indx(stack)=il
          indy(stack)=jm
        endif
        if (workar(ir,jm).eq.1) then
          workar(ir,jm)=label
          stack=stack+1
          indx(stack)=ir
          indy(stack)=jm
        endif
        if (workar(im,ju).eq.1) then
          workar(im,ju)=label
          stack=stack+1
          indx(stack)=im
          indy(stack)=ju
        endif
        if (workar(im,jd).eq.1) then
          workar(im,jd)=label
          stack=stack+1
          indx(stack)=im
          indy(stack)=jd
        endif
        if (workar(il,jd).eq.1) then
          workar(il,jd)=label
          stack=stack+1
          indx(stack)=il
          indy(stack)=jd
        endif
        if (workar(ir,jd).eq.1) then
          workar(ir,jd)=label
          stack=stack+1
          indx(stack)=ir
          indy(stack)=jd
        endif
        if (workar(il,ju).eq.1) then
          workar(il,ju)=label
          stack=stack+1
          indx(stack)=il
          indy(stack)=ju
        endif
        if (workar(ir,ju).eq.1) then
          workar(ir,ju)=label
          stack=stack+1
          indx(stack)=ir
          indy(stack)=ju
        endif

        if (stack.eq.0) exit wstack

      enddo wstack

      outar=workar(1:nx,1:ny)

c     -----------------------------------------------
      end subroutine connect
c     -----------------------------------------------

c     -----------------------------------------------------------------
c     Gradient
c     -----------------------------------------------------------------

      subroutine grad(gradx,grady,scalar,
     >                vert,xmin,ymin,dx,dy,nx,ny,nz,mdv)

c     Calculate the two vector components <gradx> and <grady> of the
c     scalar field <scalar>.  The vertical coordinate is specified in
c     <vert>, the grid in <xmin,ymin,dx,dy,nx,ny,nz,mdv>.

      implicit none

c     Declaration of subroutine parameters
      integer nx,ny,nz
      real    gradx(nx,ny,nz)
      real    grady(nx,ny,nz)
      real    scalar(nx,ny,nz)
      real    vert(nx,ny,nz)
      real    xmin,ymin,dx,dy
      real    mdv

c     Calculate the derivatives in x and y direction
      call deriv (gradx,scalar,'x',vert,xmin,ymin,dx,dy,nx,ny,nz,mdv)
      call deriv (grady,scalar,'y',vert,xmin,ymin,dx,dy,nx,ny,nz,mdv)

      end

c     -----------------------------------------------------------------
c     Create the aura around a 3d field
c     -----------------------------------------------------------------

      subroutine aura (gri,aur,dir,nx,ny,nz,xmin,ymin,dx,dy,mdv)

c     Create a one-point aura around the grid, in order to avoid nasty
c     problems when calculating fields at boundary points

      implicit none

c     Declaration of subroutine parameters
      integer nx,ny,nz
      real    gri(nx,ny,nz)
      real    aur(0:nx+1,0:ny+1,0:nz+1)
      integer dir
      real    xmin,ymin,dx,dy
      real    mdv

c     Numerical and physical parameters
      real       eps
      parameter  (eps=0.01)

c     Auxiliary variables
      integer i,j,k
      real    xmax,ymax
      integer domx,domy
      real    mean,count
      integer kmin,kmax

c     Set the z dimension for the output array
      if (nz.gt.1) then
         kmin=0
         kmax=nz+1
      elseif (nz.eq.1) then
         kmin=0
         kmax=2
      endif

c     Determine the x topology of the grid
c     1: periodic, not closed;
c     2: periodic, closed;
c     0: not periodic (and therefore not closed)
      xmax=xmin+real(nx-1)*dx
      ymax=ymin+real(ny-1)*dy
      if (abs(xmax-xmin-360.).lt.eps) then
         domx=2
      elseif (abs(xmax-xmin-360.+dx).lt.eps) then
         domx=1
      else
         domx=0
      endif

c     Determine the y topology of the grid
c     1    : neither north, nor south pole;
c     mod 2: exactly at south pole (closed)
c     mod 3: exactly at north pole (closed)
c     mod 5: one grid point north of south pole (not closed)
c     mod 7: one grid point south of north pole (not closed)
      domy=1
      if (abs(ymin+90.).lt.eps) then
         domy=2
      endif
      if (abs(ymax-90.).lt.eps) then
         domy=domy*3
      endif
      if (abs(ymin+90.-dy).lt.eps) then
         domy=domy*5
      endif
      if (abs(ymin-90.+dy).lt.eps) then
         domy=domy*7
      endif

c     Forward transformation (create aura)
      if (dir.eq.1) then

c        Copy the interior part
         aur(1:nx,1:ny,1:nz)=gri(1:nx,1:ny,1:nz)

c        Upper boundary
         if (nz.gt.1) then
            do i=1,nx
               do j=1,ny
                  if ((abs(aur(i,j,  nz)-mdv).gt.eps).and.
     >                (abs(aur(i,j,nz-1)-mdv).gt.eps)) then
                     aur(i,j,nz+1) = 2.*aur(i,j,nz) - aur(i,j,nz-1)
                  else
                     aur(i,j,nz+1) = mdv
                  endif
               enddo
            enddo
         else
            do i=1,nx
               do j=1,ny
                  aur(i,j,2)=aur(i,j,1)
               enddo
            enddo
         endif

c        Lower boundary
         if (nz.gt.1) then
            do i=1,nx
               do j=1,ny
                  if ((abs(aur(i,j,1)-mdv).gt.eps).and.
     >                (abs(aur(i,j,2)-mdv).gt.eps)) then
                     aur(i,j,0) = 2.*aur(i,j,1) - aur(i,j,2)
                  else
                     aur(i,j,0) = mdv
                  endif
               enddo
            enddo
         else
            do i=1,nx
               do j=1,ny
                  aur(i,j,0)=aur(i,j,1)
               enddo
            enddo
         endif

c        Northern and southern boundary, not near the poles
         if (mod(domy,1).eq.0) then
            do i=1,nx
               do k=kmin,kmax
                  if ((abs(aur(i,  ny,k)-mdv).gt.eps).and.
     >                (abs(aur(i,ny-1,k)-mdv).gt.eps)) then
                     aur(i,ny+1,k) = 2.*aur(i,ny,k)-aur(i,ny-1,k)
                  else
                     aur(i,ny+1,k) = mdv
                  endif
                  if  ((abs(aur(i,1,k)-mdv).gt.eps).and.
     >                 (abs(aur(i,2,k)-mdv).gt.eps)) then
                     aur(i,0,k) = 2.*aur(i,1,k)-aur(i,2,k)
                  else
                     aur(i,0,k) = mdv
                  endif
               enddo
            enddo
         endif

c        Southern boundary, one grid point north of pole
c        Set the pole point to the mean of the nearest points
         if (mod(domy,5).eq.0) then
            do k=kmin,kmax
               mean=0.
               count=0.
               do i=1,nx
                  if (abs(aur(i,1,k)-mdv).gt.eps) then
                     mean=mean+aur(i,1,k)
                     count=count+1.
                  endif
               enddo
               if (count.gt.0.) then
                  mean=mean/count
               else
                  mean=mdv
               endif
               do i=1,nx
                  aur(i,0,k) = mean
               enddo
            enddo
         endif

c        Northern boundary, one grid point south of pole
c        Set the pole point to the mean of the nearest points
         if (mod(domy,7).eq.0) then
            do k=kmin,kmax
               mean=0.
               count=0.
               do i=1,nx
                  if (abs(aur(i,ny,k)-mdv).gt.eps) then
                     mean=mean+aur(i,ny,k)
                     count=count+1.
                  endif
               enddo
               if (count.gt.0.) then
                  mean=mean/count
               else
                  mean=mdv
               endif
               do i=1,nx
                  aur(i,ny+1,k) = mean
               enddo
            enddo
         endif

c        Southern boundary, exactly at south pole
         if (mod(domy,2).eq.0) then
            do i=1,nx
               do k=kmin,kmax
                  aur(i,0,k)=mdv
               enddo
            enddo
         endif

c        Northern boundary, exactly at north pole
         if (mod(domy,3).eq.0) then
            do i=1,nx
               do k=kmin,kmax
                  aur(i,ny+1,k)=mdv
               enddo
            enddo
         endif

c        The domain is periodic in x, but not closed (no cyclic point)
         if (domx.eq.1) then
            do j=0,ny+1
               do k=kmin,kmax
                  aur(   0,j,k) = aur(nx,j,k)
                  aur(nx+1,j,k) = aur( 1,j,k)
               enddo
            enddo
         endif

c        The domain is periodic in x and closed (cyclic point)
         if (domx.eq.2) then
            do j=0,ny+1
               do k=kmin,kmax
                  aur(   0,j,k) = aur(nx-1,j,k)
                  aur(nx+1,j,k) = aur(   2,j,k)
               enddo
            enddo
         endif

c        The domain is not periodic in x
         if (domx.eq.0) then
            do j=0,ny+1
               do k=kmin,kmax
                  if ((abs(aur(1,j,k)-mdv).gt.eps).and.
     >                (abs(aur(2,j,k)-mdv).gt.eps)) then
                     aur(0,j,k) = 2.*aur(1,j,k) - aur(2,j,k)
                  else
                     aur(0,j,k) = mdv
                  endif
                  if ((abs(aur(  nx,j,k)-mdv).gt.eps).and.
     >                (abs(aur(nx-1,j,k)-mdv).gt.eps)) then
                     aur(nx+1,j,k) = 2.*aur(nx,j,k) - aur(nx-1,j,k)
                  else
                     aur(nx+1,j,k) = mdv
                  endif
               enddo
            enddo
         endif

      endif

c     Backward transformation
      if (dir.eq.-1) then

         if (nz.gt.1) then
            gri(1:nx,1:ny,1:nz)=aur(1:nx,1:ny,1:nz)
         elseif (nz.eq.1) then
            gri(1:nx,1:ny,1)=aur(1:nx,1:ny,1)
         endif

      endif

      end

c     -----------------------------------------------------------------
c     Horizontal and vertical derivatives for 3d fields
c     -----------------------------------------------------------------

      subroutine deriv (dfield,field,direction,
     >                  vert,xmin,ymin,dx,dy,nx,ny,nz,mdv)

c     Calculate horizontal and vertical derivatives of the 3d field <field>.
c     The direction of the derivative is specified in <direction>
c         'x','y'          : Horizontal derivative in x and y direction
c         'p','z','t','m'  : Vertical derivative (pressure, height, theta, model)
c     The 3d field <vert> specifies the isosurfaces along which the horizontal
c     derivatives are calculated or the levels for the vertical derivatives. If
c     horizontal derivatives along model levels should be calculated, pass an
c     index array <vert(i,j,k)=k>.

      implicit none

c     Input and output parameters
      integer    nx,ny,nz
      real       dfield(nx,ny,nz)
      real       field(nx,ny,nz)
      real       vert(nx,ny,nz)
      character  direction
      real       xmin,ymin,dx,dy
      real       mdv

c     Numerical and physical parameters
      real       pi180
      parameter  (pi180=3.141592654/180.)
      real       deltay
      parameter  (deltay=1.111775) ! transformation of units in 1/100km
      real       zerodiv
      parameter  (zerodiv=0.00000001)
      real       eps
      parameter  (eps=0.01)

c     Auxiliary variables
      integer    i,j,k
      real       vmin,vmax
      real       scale,lat
      real       vu,vl,vuvl,vlvu
      real       df(0:nx+1,0:ny+1,0:nz+1)
      real       f(0:nx+1,0:ny+1,0:nz+1)
      real       v(0:nx+1,0:ny+1,0:nz+1)

c     Create the aura around the grid for fast boundary handling
      call aura (field,f,1,nx,ny,nz,xmin,ymin,dx,dy,mdv)
      call aura (vert, v,1,nx,ny,nz,xmin,ymin,dx,dy,mdv)

c     Vertical derivative
      if ((direction.eq.'z').or.
     >    (direction.eq.'th').or.
     >    (direction.eq.'p').or.
     >    (direction.eq.'m').and.
     >    (nz.gt.1)) then

c        Finite differencing
         do i=1,nx
            do j=1,ny
               do k=1,nz

                  if ((abs(f(i,j,k+1)-mdv).gt.eps).and.
     >                (abs(f(i,j,k-1)-mdv).gt.eps).and.
     >                (abs(v(i,j,k+1)-mdv).gt.eps).and.
     >                (abs(v(i,j,k-1)-mdv).gt.eps).and.
     >                (abs(v(i,j,  k)-mdv).gt.eps)) then

                     vu = v(i,j,  k)-v(i,j,k+1)
                     vl = v(i,j,k-1)-v(i,j,  k)
                     vuvl = vu/(vl+zerodiv)
                     vlvu = 1./(vuvl+zerodiv)

                     df(i,j,k) = 1./(vu+vl)
     >                           * (vuvl*(f(i,j,k-1)-f(i,j,  k))
     >                           +  vlvu*(f(i,j,  k)-f(i,j,k+1)))

                  else
                     df(i,j,k) = mdv
                  endif

               enddo
            enddo
         enddo

c     Horizontal derivative in the y direction: 3d
      elseif ((direction.eq.'y').and.(nz.gt.1)) then

c        Scale factor for derivative in 100km
         scale=1./(2.*dy*deltay)

c        Finite differencing
         do i=1,nx
            do j=1,ny
               do k=1,nz

                 if ((abs(f(i,j+1,k)-mdv).gt.eps).and.
     >               (abs(f(i,j-1,k)-mdv).gt.eps).and.
     >               (abs(f(i,j,k+1)-mdv).gt.eps).and.
     >               (abs(f(i,j,k-1)-mdv).gt.eps).and.
     >               (abs(v(i,j,k+1)-mdv).gt.eps).and.
     >               (abs(v(i,j,k+1)-mdv).gt.eps).and.
     >               (abs(v(i,j,k-1)-mdv).gt.eps).and.
     >               (abs(v(i,j,k-1)-mdv).gt.eps).and.
     >               (abs(v(i,j+1,k)-mdv).gt.eps).and.
     >               (abs(v(i,j+1,k)-mdv).gt.eps).and.
     >               (abs(v(i,j-1,k)-mdv).gt.eps).and.
     >               (abs(v(i,j-1,k)-mdv).gt.eps)) then

                     df(i,j,k) =
     >                   scale*(f(i,j+1,k)-f(i,j-1,k))
     >                  -(f(i,j,k+1)-f(i,j,k-1))/(v(i,j,k+1)-v(i,j,k-1))
     >                  *scale*(v(i,j+1,k)-v(i,j-1,k))

                  else
                     df(i,j,k) = mdv
                  endif

               enddo
            enddo
         enddo

c     Horizontal derivative in the x direction: 3d
      elseif ((direction.eq.'x').and.(nz.gt.1)) then

c        Finite differencing
         do j=1,ny

c           Scale factor for derivatives in 100km (latitude dependent)
            lat=ymin+real(j-1)*dy
            scale=1./(2.*dx*deltay*cos(pi180*lat)+zerodiv)

            do i=1,nx
               do k=1,nz

                 if ((abs(f(i+1,j,k)-mdv).gt.eps).and.
     >               (abs(f(i-1,j,k)-mdv).gt.eps).and.
     >               (abs(f(i,j,k+1)-mdv).gt.eps).and.
     >               (abs(f(i,j,k-1)-mdv).gt.eps).and.
     >               (abs(v(i,j,k+1)-mdv).gt.eps).and.
     >               (abs(v(i,j,k+1)-mdv).gt.eps).and.
     >               (abs(v(i,j,k-1)-mdv).gt.eps).and.
     >               (abs(v(i,j,k-1)-mdv).gt.eps).and.
     >               (abs(v(i+1,j,k)-mdv).gt.eps).and.
     >               (abs(v(i+1,j,k)-mdv).gt.eps).and.
     >               (abs(v(i-1,j,k)-mdv).gt.eps).and.
     >               (abs(v(i-1,j,k)-mdv).gt.eps)) then

                     df(i,j,k) =
     >                   scale*(f(i+1,j,k)-f(i-1,j,k))
     >                  -(f(i,j,k+1)-f(i,j,k-1))/(v(i,j,k+1)-v(i,j,k-1))
     >                  *scale*(v(i+1,j,k)-v(i-1,j,k))

                  else
                     df(i,j,k) = mdv
                  endif

               enddo
            enddo
         enddo

c     Horizontal derivative in the y direction: 2d
      elseif ((direction.eq.'y').and.(nz.eq.1)) then

c        Scale factor for derivative in 100km
         scale=1./(2.*dy*deltay)

c        Finite differencing
         do i=1,nx
            do j=1,ny

               if ((abs(f(i,j+1,1)-mdv).gt.eps).and.
     >             (abs(f(i,j-1,1)-mdv).gt.eps)) then

                  df(i,j,1) = scale*(f(i,j+1,1)-f(i,j-1,1))

               else
                  df(i,j,1) = mdv
               endif

            enddo
         enddo

c     Horizontal derivative in the x direction: 2d
      elseif ((direction.eq.'x').and.(nz.eq.1)) then

c        Finite differencing
         do j=1,ny

c           Scale factor for derivatives in 100km (latitude dependent)
            lat=ymin+real(j-1)*dy
            scale=1./(2.*dx*deltay*cos(pi180*lat)+zerodiv)

            do i=1,nx

               if ((abs(f(i+1,j,1)-mdv).gt.eps).and.
     >             (abs(f(i-1,j,1)-mdv).gt.eps)) then

                  df(i,j,1) = scale*(f(i+1,j,1)-f(i-1,j,1))

               else
                  df(i,j,1) = mdv
               endif

            enddo
         enddo

c     Undefined direction for derivative
      else

         print*,'Invalid direction of derivative... Stop'
         stop

      endif

c     Get rid of the aura
      call aura (dfield,df,-1,nx,ny,nz,xmin,ymin,dx,dy,mdv)

      end

c     -----------------------------------------------------------------
c     Divergence
c     -----------------------------------------------------------------

      subroutine div(divgce,comp1,comp2,
     >               vert,xmin,ymin,dx,dy,nx,ny,nz,mdv)

c     Calculate the divergence <divgce> of a vector field with the
c     components <comp1> and <comp2>.  The vertical coordinate is
c     specified in <vert>, the grid in <xmin,ymin,dx,dy,nx,ny,nz,mdv>.

      implicit none

c     Declaration of subroutine parameters
      integer nx,ny,nz
      real    divgce(nx,ny,nz)
      real    comp1(nx,ny,nz),comp2(nx,ny,nz)
      real    tmp(nx,ny,nz)
      real    vert(nx,ny,nz)
      real    xmin,ymin,dx,dy
      real    mdv

c     Calculate the derivatives in x and y direction
      call deriv (tmp,comp1,'x',vert,xmin,ymin,dx,dy,nx,ny,nz,mdv)
      call deriv (divgce,comp2,'y',vert,xmin,ymin,dx,dy,nx,ny,nz,mdv)

      divgce=divgce+tmp

      end

c     -----------------------------------------------------------------
c     Computation of absolute value of 2D cartesian vector field
c     -----------------------------------------------------------------

      subroutine abs2d (xfield,yfield,absfield,nx,ny,nz,mdv)

      implicit none

      integer,intent(in) :: nx,ny,nz
      real,intent(in) :: mdv
      real,dimension(nx,ny,nz),intent(in) :: xfield,yfield
      real,dimension(nx,ny,nz),intent(out) :: absfield

      real,parameter :: eps=1.e-5

      absfield=sqrt(xfield**2+yfield**2)

      where ((abs(xfield-mdv).lt.eps).or.(abs(yfield-mdv).lt.eps))
     >  absfield=mdv

      end

c     -----------------------------------------------------------------
c     Computation of scalar product of 2d cartesian vector field
c     -----------------------------------------------------------------

      subroutine scal2d (xfield1,yfield1,xfield2,yfield2,
     >                   scalfield,nx,ny,nz,mdv)

      implicit none

      integer,intent(in) :: nx,ny,nz
      real,intent(in) :: mdv
      real,dimension(nx,ny,nz),intent(in) :: xfield1,yfield1
      real,dimension(nx,ny,nz),intent(in) :: xfield2,yfield2
      real,dimension(nx,ny,nz),intent(out) :: scalfield

      real,parameter :: eps=1.e-5

      scalfield=xfield1*xfield2+yfield1*yfield2

      where ((abs(xfield1-mdv).lt.eps).or.(abs(xfield2-mdv).lt.eps).or.
     >       (abs(yfield1-mdv).lt.eps).or.(abs(yfield2-mdv).lt.eps))
     >  scalfield=mdv

      end

c     -----------------------------------------------------------------
c     Standardization of 2d cartesian vector field
c     -----------------------------------------------------------------

      subroutine norm2d (xfield,yfield,nx,ny,nz,mdv)

c     Computation of unit 2d vectors / division by the absolute value

      implicit none

      integer,intent(in) :: nx,ny,nz
      real,intent(in) :: mdv
      real,dimension(nx,ny,nz),intent(inout) :: xfield,yfield

      real,dimension(nx,ny,nz) :: absfield

      real,parameter :: eps=1.e-6

      call abs2d(xfield,yfield,absfield,nx,ny,nz,mdv)

      where ((abs(xfield-mdv).gt.eps).and.(abs(yfield-mdv).gt.eps))
        xfield=xfield/absfield
        yfield=yfield/absfield
      endwhere

      end

c     -----------------------------------------------------------------
c     Check environment of gridpoint for zero crossings
c     -----------------------------------------------------------------

      function zcross(f33,mdv) result(zout)

c     Check adjacent gridpoints for sign
c     Return true if zero crossing present and center gridpoint
c     closest to zero compared to the others

      implicit none

      real,dimension(3,3),intent(in) :: f33
      real,intent(in) :: mdv
      real,dimension(1,9) :: fcalc
      real,parameter :: eps=1.e-6
      integer :: i
      logical :: zout

      fcalc=reshape(f33,(/ 1,9 /))

      zout=(abs(fcalc(5)-mdv).gt.eps)

      if (.not.zout) return

      where (abs(fcalc-mdv).lt.eps) fcalc=fcalc(5)

      do i=1,4
        zout=((sign(1.,fcalc(1,i)).ne.sign(1.,fcalc(1,5))).and.
     >      (abs(fcalc(1,5)).le.abs(fcalc(1,i))))
        if (zout) return
      enddo
      do i=6,9
        zout=((sign(1.,fcalc(1,i)).ne.sign(1.,fcalc(1,5))).and.
     >      (abs(fcalc(1,5)).le.abs(fcalc(1,i))))
        if (zout) return
      enddo

      end

c     -----------------------------------------------------------------
c     Determine zero contours in gridded data
c     -----------------------------------------------------------------

      subroutine zcontour (field,nx,ny,nz,xmin,ymin,
     >                     dx,dy,mdv)

c     Locate zero contour lines in a gridded scalar field
c     input/output : field
c     (contains ones for contour line and zeros elsewhere)

      implicit none

      integer,intent(in) :: nx,ny,nz
      real,intent(in) :: xmin,ymin,dx,dy,mdv
      real,dimension(nx,ny,nz),intent(inout) :: field

      integer :: i,j,k
      real,dimension(0:nx+1,0:ny+1,0:nz+1) :: f
      real,parameter :: eps=1.e-6

      logical,external :: zcross

      call aura(field,f,1,nx,ny,nz,xmin,ymin,dx,dy,mdv)

      field=0.

      do k=1,nz
        do j=1,ny
          do i=1,nx
            if (zcross(f(i-1:i+1,j-1:j+1,k),mdv)) field(i,j,k)=1.
          enddo
        enddo
      enddo

      end

c     -----------------------------------------------------------------
c     Calculate wind direction in degrees
c     -----------------------------------------------------------------

      subroutine winddir(nx,ny,nz,mdv,u3d,v3d,dir3d)

      ! calculation of wind direction of u and v field

      integer,intent(in) :: nx,ny,nz

      real,intent(in) :: mdv

      real,dimension(nx,ny,nz),intent(in) :: u3d,v3d

      real,dimension(nx,ny,nz) :: help

      real,dimension(nx,ny,nz),intent(out) :: dir3d

      real :: pi,twopi,temp

      integer :: i,j,k

      real,parameter :: eps=1.e-5

      pi=4.*atan(1.)
      twopi=8.*atan(1.)

      do k=1,nz
        do j=1,ny
          do i=1,nx
            help(i,j,k)=atan(u3d(i,j,k)/v3d(i,j,k))
          enddo
        enddo
      enddo

      where (v3d>0.) help=pi+help
      where ((u3d>0.).and.(v3d<0.)) help=twopi+help

      dir3d=nint(360.*help/twopi)

      where (dir3d==360) dir3d=0

      where ((v3d==0.).and.(u3d<0.)) dir3d=90
      where ((v3d==0.).and.(u3d>0.)) dir3d=270
      where ((u3d==0.).and.(v3d<0.)) dir3d=0
      where ((u3d==0.).and.(v3d>0.)) dir3d=180
      where ((u3d==0.).and.(v3d==0.)) dir3d=mdv

      where ((abs(u3d-mdv).lt.eps).or.(abs(v3d-mdv).lt.eps))
     >      dir3d=mdv

      end

c     -----------------------------------------------------------------
      SUBROUTINE T_PHUV2LLUV(PHU,PHV,LLU,LLV,PHX,PHY)
c     -----------------------------------------------------------------

c     declaration of input/output vars:
      REAL        PHU,PHV,LLU,LLV,PHX,PHY

c     declaration of internal vars:
      real        phip,lamp,phipol,cosphip,sinphip,coslamp

      real      g2r
      parameter(g2r=0.0174533)

      real      pollon,pollat
      common /pole/ pollon,pollat

c     convert to radian:
      phip=PHY*g2r
      lamp=PHX*g2r
      phipol=pollat*g2r

c     abbrevations in rotpol.icl:
      zcospol=Cos(phipol)
      zsinpol=Sin(phipol)

c     other abbrevations:
      cosphip=Cos(phip)
      sinphip=Sin(phip)
      coslamp=Cos(lamp)

c     calculate the lat/lon u component:

      LLU=8*Cos(ASin(Coslamp*zcospol*Cosphip +
     -     zsinpol*Sinphip))*
     -   (-2*PHV*Sin(lamp - phipol) - 2*PHV*Sin(lamp + phipol) -
     -     PHU*Sin(lamp - phipol - phip) -
     -     2*PHU*Sin(phipol - phip) -
     -     PHU*Sin(lamp + phipol - phip) + PHU*Sin(lamp -
     -     phipol + phip) -
     -     2*PHU*Sin(phipol + phip) + PHU*Sin(lamp + phipol + phip))/
     -  (-20 + 4*Cos(2*lamp) + 2*Cos(2*(lamp - phipol)) -
     -     4*Cos(2*phipol) +
     -    2*Cos(2*(lamp + phipol)) - 4*Cos(lamp - 2*phipol - 2*phip) +
     -    4*Cos(lamp + 2*phipol - 2*phip) + 2*Cos(2*(lamp - phip)) +
     -    Cos(2*(lamp - phipol - phip)) + 6*Cos(2*(phipol - phip)) +
     -    Cos(2*(lamp + phipol - phip)) - 4*Cos(2*phip) +
     -    2*Cos(2*(lamp + phip)) + Cos(2*(lamp - phipol + phip)) +
     -    6*Cos(2*(phipol + phip)) + Cos(2*(lamp + phipol + phip)) +
     -    4*Cos(lamp - 2*phipol + 2*phip) -
     -    4*Cos(lamp + 2*phipol + 2*phip))

c     calculate the lat/lon v component:

      LLV=(PHV*Cosphip*zsinpol +
     -     zcospol*(-(PHU*Sin(lamp)) -
     -     PHV*Coslamp*Sinphip))/
     -     Sqrt(1 - (Coslamp*zcospol*Cosphip +
     -     zsinpol*Sinphip)**2)

      RETURN
      END

c     -----------------------------------------------------------------
c     Apply a special connecting filter for extraction of lines
c     -----------------------------------------------------------------

      subroutine linesplus(field,nx,ny,nz,xmin,ymin,dx,dy,mdv)

      ! connecting lines in 1/0 field

      implicit none

      integer,intent(in) :: nx,ny,nz
      real,intent(in) :: xmin,ymin,dx,dy,mdv
      real,dimension(nx,ny,nz),intent(inout) :: field

      integer :: i,j,k
      real,dimension(0:nx+1,0:ny+1,0:nz+1) :: f
      real,parameter :: eps=1.e-6

      call aura(field,f,1,nx,ny,nz,xmin,ymin,dx,dy,mdv)

      do k=1,nz
        do j=1,ny
          do i=1,nx
            if (f(i,j,k).gt.eps) cycle
            if ((f(i-1,j,k).gt.eps).and.(f(i+1,j,k).gt.eps)) then
              field(i,j,k)=1
              cycle
            endif
            if ((f(i,j-1,k).gt.eps).and.(f(i,j+1,k).gt.eps)) then
              field(i,j,k)=1
              cycle
            endif
            if ((f(i+1,j+1,k).gt.eps).and.(f(i-1,j-1,k).gt.eps)) then
              field(i,j,k)=1
              cycle
            endif
            if ((f(i+1,j-1,k).gt.eps).and.(f(i-1,j+1,k).gt.eps)) then
              field(i,j,k)=1
              cycle
            endif
            if ((f(i+1,j+1,k).gt.eps).and.
     >          ((f(i,j-1,k).gt.eps).or.(f(i-1,j,k).gt.eps))) then
              field(i,j,k)=1
              cycle
            endif
            if ((f(i+1,j-1,k).gt.eps).and.
     >          ((f(i,j+1,k).gt.eps).or.(f(i-1,j,k).gt.eps))) then
              field(i,j,k)=1
              cycle
            endif
            if ((f(i-1,j+1,k).gt.eps).and.
     >          ((f(i,j-1,k).gt.eps).or.(f(i+1,j,k).gt.eps))) then
              field(i,j,k)=1
              cycle
            endif
            if ((f(i-1,j-1,k).gt.eps).and.
     >          ((f(i,j+1,k).gt.eps).or.(f(i+1,j,k).gt.eps))) then
              field(i,j,k)=1
            endif
          enddo
        enddo
      enddo

      end

max_stars_repo_name: ruestefa/stormtrack

max_issues_repo_name: ruestefa/stormtrack

max_forks_repo_name: ruestefa/stormtrack


text:
      REAL FUNCTION RNDRL (LO, HI)

C  RNDRL returns a random real number, X, uniformly distributed,
C  such that LO .LE. X .LT. HI.

      REAL       RND01, LO, HI

      RNDRL = RND01()*(HI-LO) + LO

      END

ext: f


text:
subroutine invertdielf(pola,polb,ipwndx,npwndx,ntpwndx,ipwc,npwc,nwpt,nqpt,lossfn)
implicit none
integer :: npwndx,ntpwndx,npwc,nwpt,nqpt
integer :: ipwndx(2,ntpwndx),ipwc(3,npwc)
double complex :: dielf(npwndx),pola(npwndx),polb(npwndx),lossfn(ntpwndx)
integer :: iqpt,iwpt,ipw1,ipw2,iipw
double complex :: A(npwc,npwc),B(npwc,npwc),AA(npwc,npwc)

  do ipw1=1,npwc
  do ipw2=1,npwc
    A(ipw1,ipw2)=0.d0
  enddo
  enddo
  do iipw=1,npwndx
    ipw1=ipwndx(1,iipw)
    ipw2=ipwndx(2,iipw)
    if (ipw1.ne.ipw2) then
      A(ipw1,ipw2)=-pola(iipw)-polb(iipw)
      A(ipw2,ipw1)=-conjg(pola(iipw))+conjg(polb(iipw))
    else
      A(ipw1,ipw2)=1.d0-dble(pola(iipw))-dimag(polb(iipw))*(0.d0,1.d0)
    endif
    dielf(iipw)=A(ipw1,ipw2)
  enddo
  do ipw1=1,npwc
  do ipw2=1,npwc
    AA(ipw1,ipw2)=A(ipw1,ipw2)
  enddo
  enddo
  call inversec(A,B,npwc,npwc)
  do iipw=1,ntpwndx
    ipw1=ipwndx(1,iipw)
    ipw2=ipwndx(2,iipw)
    lossfn(iipw)=B(ipw1,ipw2)
  enddo

return
end subroutine invertdielf

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

subroutine invertdielfT(pola,polb,ipwndx,npwndx,ntpwndx,ipwc,npwc,nwpt,nqpt,lossfn)
use geometry
implicit none
integer :: npwndx,ntpwndx,npwc,nwpt,nqpt
integer :: npwct
integer :: ipwndx(2,ntpwndx),ipwc(3,npwc)
double complex :: pola(9,npwndx),polb(9,npwndx),lossfn(9,ntpwndx)
integer :: iqpt,iqt,iwpt,ipw1,ipw2,ipw3,iipw,ii,jj,kk,ll,ix,jx
double precision :: blatinv(3,3)
double complex :: dummy(3,3), dummy2(3,3), det
double complex :: A(3*npwc,3*npwc),B(3*npwc,3*npwc),AA(npwc,npwc),BB(npwc,npwc)
double complex :: rllossfn(3,3,npwc,npwc),rlDielf(3,3,npwc,npwc),rlAlpha(3,3,npwc,npwc) ! expressed in reciprocal lattice
double complex :: cartlossfn(3,3,npwc,npwc),cartDielf(3,3,npwc,npwc),cartAlpha(3,3,npwc,npwc) ! expressed in cartesian coordinates
double precision :: prefact
double complex :: testlossfn(3,3,npwc,npwc),testD(3,3,npwc,npwc),testdummy(3,3,npwc,npwc),testB(3,3,npwc,npwc)

! Polarizability in reciprocal lattice
  do ipw1=1,npwc
  do ipw2=1,npwc
    do ii=1,3
    do jj=1,3
      rlAlpha(ii,jj,ipw1,ipw2)=(0.d0,0.d0)
    enddo
    enddo
  enddo
  enddo
  do iipw=1,npwndx
    ipw1=ipwndx(1,iipw)
    ipw2=ipwndx(2,iipw)
    do ii=1,3
    do jj=1,3
      iqt=ii+3*(jj-1)
      if (ipw1.ne.ipw1) then
        rlAlpha(ii,jj,ipw1,ipw2)=-pola(iqt,iipw)-polb(iqt,iipw)
        rlAlpha(ii,jj,ipw2,ipw1)=-conjg(pola(iqt,iipw))+conjg(polb(iqt,iipw))
      else
        rlAlpha(ii,jj,ipw1,ipw2)=-dble(pola(iqt,iipw))-dimag(polb(iqt,iipw))*(0.d0,1.d0)
      endif
    enddo
    enddo
  enddo
! Polarizability and dielectric function in cartesian coordiantes
  do ipw1=1,npwc
  do ipw2=1,npwc
    do ii=1,3
    do jj=1,3
      cartAlpha(ii,jj,ipw1,ipw2)=(0.d0,0.d0)
      do kk=1,3
      do ll=1,3
        cartAlpha(ii,jj,ipw1,ipw2)=cartAlpha(ii,jj,ipw1,ipw2)+blat(kk,ii)*rlAlpha(kk,ll,ipw1,ipw2)*blat(ll,jj)
      enddo
      enddo
      cartDielf(ii,jj,ipw1,ipw2)=cartAlpha(ii,jj,ipw1,ipw2)
    enddo
    enddo
    do ii=1,3
      if (ipw1.eq.ipw2) cartDielf(ii,ii,ipw1,ipw2)=cartDielf(ii,ii,ipw1,ipw2)+1.d0
    enddo
  enddo
  enddo
! Invert, lossfn in cartesian coordinates
  do ii=1,3
  do jj=1,3
    do ipw1=1,npwc
    do ipw2=1,npwc
      ix=(ipw1-1)*3+ii
      jx=(ipw2-1)*3+jj
      A(ix,jx)=cartDielf(ii,jj,ipw1,ipw2)
!      AA(ipw1,ipw2)=cartDielf(ii,jj,ipw1,ipw2)
    enddo
    enddo
!    call inversec(AA,BB,npwc,npwc)
  enddo
  enddo
  call inversec(A,B,3*npwc,3*npwc)
  do ii=1,3
  do jj=1,3
    do ipw1=1,npwc
    do ipw2=1,npwc
      ix=(ipw1-1)*3+ii
      jx=(ipw2-1)*3+jj
      cartlossfn(ii,jj,ipw1,ipw2)=B(ix,jx)
!      cartlossfn(ii,jj,ipw1,ipw2)=BB(ipw1,ipw2)
    enddo
    enddo
  enddo
  enddo
! lossfn in reciprocal lattice
  prefact = 1.d0/(2.d0*acos(-1.d0))**2
  do ipw1=1,npwc
  do ipw2=1,npwc
    do ii=1,3
    do jj=1,3
      rllossfn(ii,jj,ipw1,ipw2)=(0.d0,0.d0)
      do kk=1,3
      do ll=1,3
        rllossfn(ii,jj,ipw1,ipw2)=rllossfn(ii,jj,ipw1,ipw2)+rprimd(ii,kk)*cartlossfn(kk,ll,ipw1,ipw2)*rprimd(jj,ll)
      enddo
      enddo
      rllossfn(ii,jj,ipw1,ipw2)=rllossfn(ii,jj,ipw1,ipw2)*prefact
    enddo
    enddo
  enddo
  enddo
  do iipw=1,npwndx
    ipw1=ipwndx(1,iipw)
    ipw2=ipwndx(2,iipw)
    do ii=1,3
    do jj=1,3
      iqt=ii+3*(jj-1)
      lossfn(iqt,iipw)=rllossfn(ii,jj,ipw1,ipw2)
    enddo
    enddo
  enddo

return
end subroutine invertdielfT

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

subroutine locateelement(iqq,ipw1,ipw2,ngkpt,iqsymndx,npwc,npwx,nsym,pwsymndx,invpw2ndx,iipw)
implicit none
integer :: nsym,npwc,npwx,ngkpt(3),pwcmin(3),pwcmax(3)
integer :: iqq(3),ipw1,ipw2,isym,ipw3,ipw4,iipw
integer :: iqsymndx(ngkpt(1),ngkpt(2),ngkpt(3))
integer :: pwsymndx(npwc,2*nsym)
integer :: invpw2ndx(npwx,npwx)

isym=iqsymndx(iqq(1),iqq(2),iqq(3))
ipw3=pwsymndx(ipw1,isym)
ipw4=pwsymndx(ipw2,isym)
if (ipw3.ne.0.and.ipw4.ne.0) then
  iipw=invpw2ndx(ipw3,ipw4)
else
  iipw=0
endif

return
end subroutine locateelement

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

subroutine mkpwsymndx(ipwc,npwc,syminv,nsym,invpwndx,pwcmin,pwcmax,pwsymndx)
implicit none
integer :: nsym,npwc,pwcmin(3),pwcmax(3)
integer :: syminv(3,3,nsym),ipwc(3,npwc)
integer :: invpwndx(pwcmin(1):pwcmax(1),pwcmin(2):pwcmax(2),pwcmin(3):pwcmax(3))
integer :: pwsymndx(npwc,2*nsym)
integer :: igg(3),iggr(3)
integer :: ii,jj,isym,ipw,isign,iisym

do ipw=1,npwc
  igg=ipwc(:,ipw)
  do isym=1,2*nsym
    iisym=mod(isym-1,nsym)+1
    isign=-(2*((isym-1)/nsym)-1)
    iggr=(/0,0,0/)
    do ii=1,3
    do jj=1,3
      iggr(ii)=iggr(ii)+syminv(jj,ii,iisym)*isign*igg(jj)
    enddo
    enddo
    if (iggr(1).ge.pwcmin(1).and.iggr(1).le.pwcmax(1).and. &
&       iggr(2).ge.pwcmin(2).and.iggr(2).le.pwcmax(2).and. &
&       iggr(3).ge.pwcmin(3).and.iggr(3).le.pwcmax(3)) then
      pwsymndx(ipw,isym)=invpwndx(iggr(1),iggr(2),iggr(3))
    else
      pwsymndx(ipw,isym)=0
    endif
  enddo
enddo

return
end subroutine mkpwsymndx

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

subroutine mkinvpwndx(ipwc,npwc,invpwndx,pwcmax,pwcmin)
implicit none
integer :: npwc,pwcmax(3),pwcmin(3)
integer :: ipwc(3,npwc)
integer :: invpwndx(pwcmin(1):pwcmax(1),pwcmin(2):pwcmax(2),pwcmin(3):pwcmax(3))
integer :: ix,iy,iz,ii

do ix=pwcmin(1),pwcmax(1)
do iy=pwcmin(2),pwcmax(2)
do iz=pwcmin(3),pwcmax(3)
  invpwndx(ix,iy,iz)=0
enddo
enddo
enddo

do ii=1,npwc
  invpwndx(ipwc(1,ii),ipwc(2,ii),ipwc(3,ii))=ii
enddo

return
end subroutine mkinvpwndx

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

subroutine mkinvpw2ndx(ipwndx,napwndx,npwx,invpw2ndx)
implicit none
integer :: napwndx,npwx
integer :: ipwndx(2,napwndx)
integer :: invpw2ndx(npwx,npwx)
integer :: ipw1,ipw2,ii

do ipw1=1,npwx
do ipw2=1,npwx
  invpw2ndx(ipw1,ipw2)=0
enddo
enddo

do ii=1,napwndx
  invpw2ndx(ipwndx(1,ii),ipwndx(2,ii))=ii
enddo

return
end subroutine mkinvpw2ndx

max_stars_repo_name: dinisAbranches/nwchem

max_issues_repo_name: dinisAbranches/nwchem

max_forks_repo_name: dinisAbranches/nwchem


text:
module a
  type :: dot
    integer :: x, y, z
      contains
        procedure:: sub => sub1
  end type dot

  type :: cut
    type(dot) :: a, b
  end type cut

  type :: cutarray
    type(cut), dimension(10) :: dat
  end type cutarray

  type :: object
    integer :: id
    type(cutarray) :: ar
  end type object

  contains

  real function l(o)
    call t1%sub()
    return o%ar%dat(5)%b%x
  end function
end module a

max_stars_repo_name: LChernigovskaya/fortran-plugin

max_issues_repo_name: LChernigovskaya/fortran-plugin

max_forks_repo_name: LChernigovskaya/fortran-plugin


text:
C
C	$Id: drawpv.f,v 1.6 2008-07-27 00:23:04 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE DRAWPV (IX,IY,IND)
C
C DRAWPV INTERCEPTS THE CALL TO PLOTIT.
C FOR IND=1 (PEN DOWN) IT CHECKS IF A LINE (PART OF IT) HAS TO BE DRAWN
C OR REMOVED.
C FOR IND=0 (PEN UP) IT CHECKS IF THE PEN HAS TO BE MOVED OR IF IT IS
C ALREADY CLOSE ENOUGH TO THE WANTED POSITION.
C FOR IND=2 DO NOT MOVE PEN.JUST SET POSITION OF PEN TO (IMPOS,IMPOS).
C
      COMMON/INTPR/IPAU,FPART,TENSN,NP,SMALL,L1,ADDLR,ADDTB,MLLINE,
     1    ICLOSE
C IOFFS IS AN INTERNAL PARAMETER. IT IS REFERENCED IN  DRAWPV AND
C FDVDLD AND INITIALIZED IN DASHBDX.
C
      COMMON /SMFLAG/ IOFFS
C
C IN DSAVE3 THE STARTING POINT OF THE NEXT LINE SEGMENT TO BE DRAWN
C OR REMOVED IS SAVED.
C
      COMMON /DSAVE3/ IXSTOR,IYSTOR
C
C IX1,IY1,IX2 AND IY2 ARE USED TO STORE THE NEXT LINE SEGMENT TO BE
C MARKED, WHICH IS THE SEGMENT THAT WAS DRAWN LAST.
C
      COMMON /DSAVE4/ IX1,IY1,IX2,IY2
C
C ISKIP IS USED TO ADJUST DASHSUPR TO THE SIZE OF THE MODEL PICTURE AS
C DEFINED IN SUBROUTINE REMOVE. ISKIP IS INITIALIZED IN DASHBDX.
C
      COMMON /DSUP1/ ISKIP
C
      LOGICAL LI, LM, HID
      LOGICAL LDUMMY
      SAVE
C
      DATA LDUMMY /.TRUE./
C
C BRANCH DEPENDING ON FUNCTION TO BE PERFORMED.
C
      IIND = IND + 1
      GOTO (2,1,3), IIND
C
    1 CONTINUE
C
C *************************************
C
C INTERCEPT PLOTIT(IX,IY,1) TO CHECK IF THE LINE TO BE DRAWN (PART OF
C IT) HAS TO BE REMOVED. DRAW AND MARK LINE IF FEASIBLE.
C
C
      IF (IOFFS .EQ. 0) GO TO 10
C
C DO NOT MARK LINES IF NOT REQUIRED.
C
      IXSTOR = IX
      IYSTOR = IY
      CALL PLOTIT (IX,IY,1)
      RETURN
C
C LINE HAS TO BE DRAWN AND MARKED.
C
   10 CONTINUE
      MX = IXSTOR
      MY = IYSTOR
      IXSTOR = IX
      IYSTOR = IY
      DX = IXSTOR - MX
      DY = IYSTOR - MY
C
C CHECK BEGINNING AND END POINT OF LINE, IF THEY ARE ALREADY MARKED.
C
      CALL REMOVE (IXSTOR,IYSTOR,LI,3)
      CALL REMOVE (MX,MY,LM,3)
C
C NOTHING TO DO , IF BEGINNING AND END POINT ARE ALREADY MARKED.
C
      IF (LI .AND. LM) RETURN
C
C DRAW A POINT IF LINE IS OF LENGTH 0.
C
      IF (DX .NE. 0. .OR. DY .NE. 0.) GO TO 15
      CALL PLOTIT (IXSTOR,IYSTOR,1)
      RETURN
   15 CONTINUE
C
      IF (LI .OR. LM) GO TO 20
C
C IF NEITHER STARTING NOR END POINT IS MARKED, DRAW LINE AND MARK IT.
C
      CALL PLOTIT (IXSTOR,IYSTOR,1)
      GO TO 80
C
C EITHER BEGINNING OR END POINT OF LINE WERE ALREADY MARKED.
C
   20 D = MAX(ABS(DX),ABS(DY))
      ID = D/REAL(ISKIP)
      DX = (DX/D)*REAL(ISKIP)
      DY = (DY/D)*REAL(ISKIP)
      IF (LM) GO TO 50
C
C STARTING POINT WAS NOT MARKED. FIND THE FIRST MARKED POINT.
C
      X = MX
      Y = MY
      DO 30 I=1,ID
         X = X+DX
         Y = Y+DY
         CALL REMOVE (INT(X+.5),INT(Y+.5),HID,3)
         IF (HID) GO TO 40
   30 CONTINUE
      X = REAL(IXSTOR)+DX
      Y = REAL(IYSTOR)+DY
   40 IXP = X+.5-DX
      IYP = Y+.5-DY
C FIRST MARKED POINT FOUND.
C DRAW LINE, SET PEN, MARK LINE.
      CALL PLOTIT (IXP,IYP,1)
      CALL PLOTIT (IXSTOR,IYSTOR,0)
      GO TO 80
C
C     FIRST PART IS HIDDEN
C       FIND THE FIRST VISIBLE POINT
C
   50 X = MX
      Y = MY
      DO 60 I=1,ID
         X = X+DX
         Y = Y+DY
         CALL REMOVE (INT(X+.5),INT(Y+.5),HID,3)
         IF (.NOT.HID) GO TO 70
   60 CONTINUE
      X = IXSTOR
      Y = IYSTOR
   70 IXP = X+.5
      IYP = Y+.5
C FIRST VISIBLE POINT FOUND.
C SET PEN, DRAW LINE, MARK LINE.
      CALL PLOTIT (IXP,IYP,0)
      CALL PLOTIT (IXSTOR,IYSTOR,1)
C
C
   80 CONTINUE
      CALL MARKL (IX1,IY1,IX2,IY2)
      IX1 = MX
      IY1 = MY
      IX2 = IXSTOR
      IY2 = IYSTOR
      RETURN
C
    2 CONTINUE
C
C *************************************
C
C INTERCEPT PLOTIT(IX,IY,0) TO SET THE STARTING POINT FOR THE NEXT
C LINE TO BE DRAWN OR REMOVED.
C
C CHECK IF PEN IS ALREADY CLOSE ENOUGH TO THE WANTED POSITION.
      DIFF = REAL(ABS(IXSTOR-IX)+ABS(IYSTOR-IY))
      IF (DIFF .LE. REAL(ICLOSE)) GO TO 110
C
      IXSTOR = IX
      IYSTOR = IY
      CALL PLOTIT (IXSTOR,IYSTOR,0)
      RETURN
C
    3 CONTINUE
C
C ****************************
C
C DO NOT MOVE PEN. JUST SET POSITION OF PEN TO (IMPOS,IMPOS).
C
      IXSTOR = IX
      IYSTOR = IY
C
  110 CONTINUE
C
      RETURN
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
C  =====================================================================
C  pgm: SHHRFX .. Constrain the hour range to 0-23 as opposed to 0-24
C
C  use:     CALL SHHRFX(IYR,IMO,IDA,IHR)
C
C  i/o: IYR .......... zulu date year number - INT
C  i/o: IMO .......... zulu date month number - INT
C  i/o: IDA .......... zulu date day number - INT
C  i/o: IHR .......... zulu date hour number (0-23 for output) - INT
C
C  rqd: SHGJUL,SHIJUL,SHGCAL
C
C  cmt: The original SHEF Code allowed hours in the range 0-24.  This
C  cmt:  filter adjusts an output date that has an hour number of 24 to
C  cmt:  a date with an hour number of 0 by incrementing it by one day.
C  cmt:  It does this for output only, without altering the logic of the
C  cmt:  parsing code.  Note that there are conflicts in the comments
C  cmt:  throughout the parsing code about whether the range is
C  cmt:  0-24 or 0-23.
C  =====================================================================
      SUBROUTINE SHHRFX(IYR,IMO,IDA,IHR)

      EXTERNAL           SHGJUL,SHIJUL,SHGCAL

      INTEGER            IYR,IMO,IDA,IHR,JUL
C
C    ================================= RCS keyword statements ==========
      CHARACTER*68     RCSKW1,RCSKW2
      DATA             RCSKW1,RCSKW2 /                                 '
     .$Source: /fs/hseb/ob5/rfc/ofs/src/shefpars_driv/RCS/shhrfx.f,v $
     . $',                                                             '
     .$Id: shhrfx.f,v 1.2 2000/03/01 16:08:00 dws Exp $
     . $' /
C    ===================================================================
C

C        If the hour is 24, convert the yy-mm-dd to julda (SHGJUL), add
C         one to the julda (SHIJUL), and recover the yy-mm-dd (SHGCAL)

        IF (IHR .EQ. 24) THEN
          IHR = 0
          CALL SHGJUL(JUL,IYR,IMO,IDA)
          CALL SHIJUL(JUL,IYR,1)
          CALL SHGCAL(JUL,IYR,IMO,IDA)
        ENDIF

      RETURN
      END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
!
! tests bug 2554037 regarding include path
!   - this should run correctly without giving an include path
!
include "bug-2554037.h"
end

max_stars_repo_name: OpenFortranProject/ofp-sdf

max_issues_repo_name: OpenFortranProject/ofp-sdf

max_forks_repo_name: OpenFortranProject/ofp-sdf


text:
!
! The Laboratory of Algorithms
!
! The MIT License
!
! Copyright 2011-2015 Andrey Pudov.
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the 'Software'), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in
! all copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
! THE SOFTWARE.
!

! Longest Collatz sequence
!
! The following iterative sequence is defined for the set of positive integers:
!
!    n  > n/2 (n is even)
!    n  > 3n + 1 (n is odd)
!
! Using the rule above and starting with 13, we generate the following sequence:
!
!    13 > 40 > 20 > 10 > 5 > 16 > 8 > 4 > 2 > 1
!
! It can be seen that this sequence (starting at 13 and finishing at 1) contains
! 10  terms. Although it has not been proved yet (Collatz Problem), it is
! thought that all starting numbers finish at 1.
!
! Which starting number, under one million, produces the longest chain?
!
! NOTE: Once the chain starts the terms are allowed to go above one million.
module MPEProblem14

    implicit none
    private

    type, public :: TPEProblem14
    contains
        procedure, nopass :: present
    end type
contains
    subroutine present
        write (*, '(A)') 'Problem 14. Longest Collatz sequence.'

        write (*, '(A, I)') 'Chain 1: ', chain1()
    end subroutine

    function chain1
        integer :: chain1
        integer :: index
        integer :: value
        integer :: max

        ! initial values
        chain1 = 0
        index  = 0
        value  = 0
        max    = 0

        do index = 1, 1000000
            value = get_length(index)

            if (value .gt. max) then
                max    = value
                chain1 = index
            end if
        end do
    end function

    function get_length(number)
        integer, intent(in) :: number

        integer :: get_length
        integer*8 :: value

        ! initial value
        get_length = 1
        value      = number

        do while (value .gt. 1)
            if (mod(value, 2) .eq. 0) then
                value = value / 2
            else
                value = (value * 3) + 1
            end if

            get_length = get_length + 1
        end do
    end function
end module

max_stars_repo_name: andreypudov/Algorithms

max_issues_repo_name: jlokimlin/Algorithms

max_forks_repo_name: jlokimlin/Algorithms


text:
subroutine lsf_gh(x,n,xcenter,par,npar,lsf)

!replicates APOGEE's lsf_gh.pro

use share, only: dp

implicit none

!input/output
real(dp), intent(in)    :: x(n)	   ! pixel location array
real(dp), intent(in)    :: par(npar)!param. array
real(dp), intent(in)    :: xcenter ! position of the lsf center
integer, intent(in)     :: n,npar
real(dp), intent(out)   :: lsf(n)  ! lsf array

!locals
integer,parameter       :: maxporder=5
integer					:: i,j
real(dp)				:: binsize = 0.0_dp    ! width of a pixel in x units  -- par(1)
real(dp)				:: xoffset = 0.0_dp, xcenter1 = 0.0_dp
integer, parameter	    :: Horder  = 5         ! highest allowed hermite order
integer					:: AHorder             ! actual highest hermite order -- par(3)
integer					:: Porder(Horder+1)    ! polynomial order for variation with x
								               ! of each parameter -- par(4:Horder+4)
integer					:: cstart(Horder+2)
real(dp)				:: GHcoefs(Horder+1)   ! There are Porder[i]+1
								               ! coefficients for parameter i.
								               ! The Hermite parameters start with H2
								               ! since we fix H1=1.
real(dp)				:: coefarr(Horder+1,maxporder+1)
real(dp)                :: cpar(npar-(Horder+4)+1)
real(dp)				:: inpar(Horder+4)

binsize=par(1)
xoffset=par(2)
AHorder = floor(par(3))
Porder(1:AHorder+1) = floor(par(4:AHorder+4))

!write(*,*)'binsize,xoffset,AHorder=',binsize,xoffset,AHorder
!write(*,*)'Porder=',Porder(1:AHorder+1)

!checks
if (binsize > 0.1) then
        write(*,*) 'lsf_gh: WARNING'
        write(*,*) 'binsize =',binsize,'> 0  '
        write(*,*)' -- but no binning will be used in building the LSF kernel'
endif
if (maxval(Porder(1:AHorder+1)) > maxporder) then
        write(*,*) 'lsf_gh: ERROR'
        write(*,*) 'maxval(Porder) = ',maxval(Porder),' > maxporder=',maxporder
        stop
endif

cpar = par(AHorder+4+1:npar+1)
cstart(:)=0
do i=1,AHorder+1
	cstart(i+1)=cstart(i)+Porder(i)+1
enddo

!write(*,*)'cpar=',cpar
!write(*,*)'cstart=',cstart

coefarr(:,:)=0._dp
do i=1,AHorder+1
	coefarr(i,:)=0._dp
	coefarr(i,1:Porder(i)+1)= cpar(cstart(i)+1:cstart(i)+Porder(i)+1)
	!write(*,*)i,'coefarr(i,1:porder(i)+1)',coefarr(i,1:Porder(i)+1)
enddo

Xcenter1 = xcenter+xoffset

GHcoefs(:)=0._dp
do i=1,AHorder+1
	!write(*,*)'i=',i
	GHcoefs(i) = 0.0_dp
	do j=1,Porder(i)+1
		GHcoefs(i)=GHcoefs(i)+coefarr(i,j)*Xcenter1**(j-1)
	enddo
	!write(*,*)'GHcoefs(',i,')=',GHcoefs(i)
enddo

inpar(:) = 0.0_dp
inpar(1) = 1.0_dp
inpar(2) = Xcenter
inpar(3) = GHcoefs(1)
inpar(4) = 1.0_dp

!write(*,*)'xcenter,xoffset,xcenter1=',xcenter,xoffset,xcenter1

if (Horder > 0) then
	do i=1,AHorder
		inpar(4+i)=GHcoefs(i+1)
	enddo
endif

!write(*,*)'n,AHorder+4=',n,AHorder+4
!write(*,*)'inpar=',inpar(1:AHorder+4)

call gausshermite(x,n,inpar(1:AHorder+4),AHorder+4,lsf)

end subroutine lsf_gh

max_stars_repo_name: dnidever/apogee

max_issues_repo_name: dnidever/apogee

max_forks_repo_name: dnidever/apogee


text:
      PROGRAM xbeta
C     driver for routine beta
      INTEGER i,nval
      REAL beta,value,w,z
      CHARACTER text*13
      open(7,file='FNCVAL.DAT',status='OLD')
10    read(7,'(a)') text
      if (text.ne.'Beta Function') goto 10
      read(7,*) nval
      write(*,*) text
      write(*,'(1x,t5,a1,t11,a1,t24,a6,t43,a9)')
     *     'W','Z','Actual','BETA(W,Z)'
      do 11 i=1,nval
        read(7,*) w,z,value
        write(*,'(2f6.2,2e20.6)') w,z,value,beta(w,z)
11    continue
      close(7)
      END

ext: for

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


text:
      function area(r,z,ix)
c     implicit double precision (a-h,o-z)                                    dp
      dimension r(*),z(*),ix(*),xdr(4),xdz(4)
      do 10 i=1,4
      xdr(i)=r(ix(i))
   10 xdz(i)=z(ix(i))
      zd42=xdz(4)-xdz(2)
      ac1=xdr(2)*(xdz(3)-xdz(4))+xdr(3)*zd42+xdr(4)*(xdz(2)-xdz(3))
      ac2=xdr(2)*(xdz(4)-xdz(1))+xdr(4)*(xdz(1)-xdz(2))-xdr(1)*zd42
      area=ac1+ac2
      return
      end

max_stars_repo_name: imohame/LabCode

max_issues_repo_name: imohame/LabCode

max_forks_repo_name: imohame/LabCode


text:
STA Travel closed its doors on campus Dec. 13, 2013, but according to a http://dateline.ucdavis.edu/dl_detail.lasso?id14677 Dateline Newsbrief, will continue to offer phone and email service to the campus.  STA Travel is a full service travel agency catering to students, teachers and the general public.  They offer everything from flexible/refundable tickets to tours, cruises, hotels, rail passes, travel insurance, and even volunteer projects everywhere in the world.

They are the only place locally that sells the International Student ID Card (ISIC Card) which cost $26 (or about $22 if you bring a passport picture). The ISIC Card gives you both local and global discounts.  Some of our favorite discounts in Davis are 15% off total purchase at Outdoor Davis, 10% off all cash purchases at Cultive Frozen Yogurt, $1 off any burger or sandwich at Tommy Js, free small drink with purchase of an entree at El Mariachi, 10% off dinner at Rajas, 10% off order at 3rd and U, and 15% off Amtrak tickets if you buy them 3 days in advance.  To check out more discounts, check out their website http://www.myisic.com/MyISIC/index.aspx ISIC Passport Photos Passport Pictures are $10.

As a travel agency, of course it costs a little more than reserving a trip oneself. However, if youre traveling on the Universitys dime, using STA is very convenient.
  Wiki world travelers Travelling? Davis offers many Travel Agencies.

I used STA Travel to book my hostels AND flight on my short European trip. It seemed like a pretty good deal. I saved a bunch of money by going through them. Also it was extremely easy. However, they kind of overcharge on little things like International Student IDs (you dont really need them. They claim discounts with them, but a regular UCD ID will do overseas!). And my biggest gripe is the fact that their international offices DO NOT communicate with each other. So if you are in London, you cant go into their London office and inquire about your bookings. You have to call specifically the location that booked you (quite a hassle!). But this isnt really a problem if you have no qwells/questions with your booking, really. Just be aware!  Users/YawenChen Yawen

I went to the STA in Berkeley and found it to be a load of shit. I had never really traveled before, so I thought I needed reservations, but I ended up paying through the nose for a hostel that was not very nice and a student ID card that I didnt need. They also knew more or less nothing about where I was going (Paris) and proceeded to waste my time. I will say that they did help me get my HI card, but my guide book told me that I needed it as well.  Users/ArlenAbraham

20090602 16:33:17 nbsp a regular UCD ID will do overseas
: Of course such a general statement is bound to be wrong sometimes. I had my student discount denied at least once because the UCD ID does not have an expiration date printed on it. The ISIC ID does have a date on it, so that probably would have worked. You have to weigh whether youre likely to take advantage of many student discounts (e.g. by going to a lot of museums in Europe) or not. Users/JoFeuerstein

20090602 18:13:48 nbsp An international student ID can save you a LOT of money, depending on where you are, and what mode of transportation you are using.  Jo is correct about the expiration date issue.  You can get significant discounts on trains, performances, and other things you may do in Europe if you have an international student ID.  Alternatively, you may also be able to use a letter on the schools letterhead explaining that you are a student.  Get this before you travel. Users/IDoNotExist

20090627 08:17:04 nbsp STA is nice if youre going to change your ticket, you can do it over email and often the fare difference hasnt been too big of a deal. HOWEVER, if you miss a flight with STA, you are totally screwed. They are ticket consolidators which means their tickets are the first to get bumped. Also, airline people cant rebook you, you have to call STA. I missed a flight to europe and had to pay $500 to rebook it. That doesnt happen when you go through online deals.  Users/AmyGoogenspa

20120514 22:53:34 nbsp I suggest using the STA website instead of the agency. Numerous times, I caught my agent charging me more than the offered price on the website. Instead of double checking their work, it is MUCH easier to book everything yourself, including the flight, international student id card, and insurance. I tried to use an agent to save myself time, but instead spent more time realizing that I was overcharged. You have to use an agent for the packaged tours, but even that I would try to do over email so that everything is documented.  Users/IshmaelSmith

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!       NASA/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
!-----------------------------------------------------------------------
!
!
! !PROGRAM:      zeit_usage
!
! !DESCRIPTION:  This program returns usage information about the particular
!                program that was called
!
! !USAGE:        zeit_usage(progname)
!                where:  progname    Name of called program
!
      subroutine zeit_usage(progname)
      implicit none
      character*7 progname
      if(progname.eq.'zeit_pr')then
      print*,'Usage:     zeit_pr.x [-i fname] [-o out_fname] [-d] [-h]'
         print*,'-i         Indicate use of a users registry file'
         print*,'fname      Input registry file. Default is .zeit'
         print*,'-o         Indicate name of output file'
         print*,'out_fname  Output file to write current timings'
         print*,'           Default is standard output'
         print*,'-d         Option to remove registry when finished'
         print*,'-h         Option to print usage message'
      elseif(progname.eq.'zeit_ci')then
         print*,'Usage is:   zeit_ci.x [-r fname] [-v] [name]'
         print*,'-r fname    Time register file name,'
         print*,'            default is .zeit'
         print*,'-v          Verbose mode: prints wallclock'
         print*,'            time'
         print*,'name        Name of program to be timed'
      elseif(progname.eq.'zeit_co')then
         print*,'Usage is:   zeit_co.x [-r fname] [-v] [name]'
         print*,'-r fname    Time register file name,'
         print*,'            default is .zeit'
         print*,'-v          Verbose mode: prints wallclock'
         print*,'            time'
         print*,'name        Name of program to be timed'
      endif
      return
      stop
      end subroutine zeit_usage


text:
C-------GF 08.11.1980------------------------------ P A R S
      IMPLICIT NONE
      INTEGER*2 XXBUCH,XXCODH,XXITEH,XXLINH,XXMEMH,XXPARH
      PARAMETER (XXBUCH=127)
      PARAMETER (XXCODH=128)
      PARAMETER (XXITEH=600)
      PARAMETER (XXLINH=41)
      PARAMETER (XXMEMH=1200)
      PARAMETER (XXPARH=30)
      INTEGER*2 XXPREH,XXPRIH,XXPROH,XXSEMH,XXSETH,XXSLOH
      PARAMETER (XXPREH=1400)
      PARAMETER (XXPRIH=61)
      PARAMETER (XXPROH=350)
      PARAMETER (XXSEMH=XXPREH)
      PARAMETER (XXSETH=XXMEMH)
      PARAMETER (XXSLOH=6)
      INTEGER*2 XXSPAH,XXSTAH,XXSTKH,XXSTRH,XXSYMH,XXTARH
      PARAMETER (XXSPAH=40)
      PARAMETER (XXSTAH=500)
      PARAMETER (XXSTKH=80)
      PARAMETER (XXSTRH=4*XXSTAH)
      PARAMETER (XXSYMH=500)
      PARAMETER (XXTARH=35)
      INTEGER*2 XXTRAH,XXTREH,XXDIRH,XXUDIR,XXUTRE,XXUPRI
      PARAMETER (XXTRAH=XXSETH/3)
      PARAMETER (XXTREH=40)
      PARAMETER (XXDIRH=(XXTREH+1)*(XXSLOH+1))
      PARAMETER (XXUDIR=7)
      PARAMETER (XXUTRE=6)
      PARAMETER (XXUPRI=6)
      INTEGER*2 XXULIN,XXUTAR,XXUPAR,XXUASS
      PARAMETER (XXULIN=4)
      PARAMETER (XXUTAR=3)
      PARAMETER (XXUPAR=2)
      PARAMETER (XXUASS=1)
      INTEGER*2 UPRI,ULIN,UTAR
      INTEGER*2 UPAR,UASS,UTRE
      INTEGER*2 UDIR
      INTEGER*2 PARVAL(XXPARH),PARME(XXPARH)
      INTEGER*2 PARM(3,XXPARH)
      INTEGER*2 PARLAS,PARHIB
      COMMON /AAAAAA/
     =  UPRI      ! XXUPRI: TEST-PRINTOUT TO THIS UNIT
     = ,ULIN      ! XXULIN: SOURCE-PROGRAM FROM HERE
     = ,UTAR      ! XXUTAR: TARGET-PROGRAM GOES HERE
     = ,UPAR      ! XXUPAR: READ PARAMETERS FROM THIS FILE
     = ,UASS      ! XXUASS: WRITE ASSERT-TEXTS TO THIS DIRECT-FILE
     = ,UTRE    ! XXUTRE: EXTERNAL TREE ON THIS UNIT
     = ,UDIR    ! XXUDIR: DIRECT-FILE FOR PARSING TABLE
     = ,PARLAS  ! ACTUAL NUMBER OF PARAMETERS
     = ,PARHIB  ! MAXIMAL ...
     = ,PARVAL  ! VALUES OF THE PARAMETERS
     = ,PARME   ! ENDING POSITION OF THE TEXTS
     = ,PARM      ! NAMES OF THE PARAMETERS

max_stars_repo_name: gfis/jextra

max_issues_repo_name: gfis/jextra

max_forks_repo_name: gfis/jextra


text:
      subroutine volcorr(x_solid,nsd,nen,ien,nn_solid,ne_solid,x_solid_init)
! x_solid: the potision of the solid domain at current time step
! nsd: dimension of the problem
! nen: number of nodes per element of solid mesh
! g: correction of the updated solid points
! ien: connectiviy of the solid mesh
! nn_solid: number of nodes of solid mesh
! ne_solid: number of elements of solid mesh
! e_area: the area of every element
! rhs: (det(F_{ij}))^{-1}-1 right hand side of the equations

      integer nsd
	  integer nn_solid
	  integer ne_solid
	  integer ien(ne_solid,nen)
      real(8) x_solid(nsd,nn_solid)
	  real(8) g(nsd,nn_solid)

	  real(8) x_solid_init(nsd,nn_solid)
!=============================================================
! shape funciton variable
	  integer iquad
	  integer nquad
	  real(8) xq(nsd,1)
	  real(8) wq
	  integer iq
	  real(8) sq(0:nsd,nen,1)
	  real(8) sh(0:nsd,nen)
	  real(8) xr(nsd,nsd)
	  real(8) cf(nsd,nsd)
          real(8) sx(3,3)
	  real(8) det
! iquad: number of integration points always 1 here
! nquad: always 1
! sq: the same meaning of in the main program but only 1 in 3nd dimension here
! iq: always 1 here just for convenience
!===============================================================
! variables to calculate rhs & e_area
      real(8) e_area(ne_solid)
	  real(8) rhs(ne_solid)
	  real(8) y(nsd,nen)
	  real(8) rs(nsd)
	  real(8) xj(nsd,nsd)
	  real(8) xji(nsd,nsd)
	  real(8) detJ
	  real(8) toxj(nsd,nsd)
	  real(8) toxji(nsd,nsd)
	  real(8) todet
	  real(8) xto(nsd,nsd)
	  real(8) xot(nsd,nsd)
	  real(8) toc(nsd,nsd)
	  real(8) detf
! det: determinant of Jacobian at current time step
! detf: det(F_{ij}^{-1})
!=======================================
! Seting up left and right hand side variables
      integer ine
	  integer inen
	  integer index
	  integer jne
	  integer jnen
	  integer indey
	  integer inn
	  real(8) x(nsd,nen)
!	  real(8) y(nsd,nen)
	  real(8) crrm(nsd*nn_solid,nsd*nn_solid)
	  real(8) crra(nsd*nn_solid)
	  real(8) s_area(nn_solid)
! crrm: the whole matirx of the left hand side
! crra: the whole matrix of the right hand side
! s_area: the sum area of the elements sounding every point
!====================================================
! varialbe to solve the linear system
      integer np
      integer mp
      integer error
      real(8) crra_cg(nsd*nn_solid)
!=============================================================
! First step: calculate the shape funciton in the parent domain, based on nsd & nen
	  iquad=1
	  iq=1
          xq(:,:)=0.0d0
	  if (nsd==2) then
	     if (nen==3) then
 	!              write(*,*) 'I am here in volcorr 84', xq(:,:)
		     call quad2d3n(iquad,nquad,xq,wq,nsd,1)
         !             write(*,*) 'I am after call quad', xq(:,:)
		 else if (nen==4) then
		     call quad2d4n(iquad,nquad,xq,wq,nsd,1)
		 end if
		 if(nen==3) then
		  sq(0,1,iq) = xq(1,iq)
		  sq(0,2,iq) = xq(2,iq)
		  sq(0,3,iq) = 1 - xq(1,iq) - xq(2,iq)
        else if (nen==4) then
		  sq(0,1,iq) = (1 - xq(1,iq)) * (1 - xq(2,iq)) / 4
		  sq(0,2,iq) = (1 + xq(1,iq)) * (1 - xq(2,iq)) / 4
		  sq(0,3,iq) = (1 + xq(1,iq)) * (1 + xq(2,iq)) / 4
		  sq(0,4,iq) = (1 - xq(1,iq)) * (1 + xq(2,iq)) / 4

		  sq(1,1,iq) = - (1 - xq(2,iq)) / 4
		  sq(1,2,iq) = + (1 - xq(2,iq)) / 4
		  sq(1,3,iq) = + (1 + xq(2,iq)) / 4
		  sq(1,4,iq) = - (1 + xq(2,iq)) / 4

		  sq(2,1,iq) = - (1 - xq(1,iq)) / 4
		  sq(2,2,iq) = - (1 + xq(1,iq)) / 4
		  sq(2,3,iq) = + (1 + xq(1,iq)) / 4
		  sq(2,4,iq) = + (1 - xq(1,iq)) / 4

        endif
	  else if (nsd==3) then
	     if (nen==4) then
		     call quad3d4n(iquad,nquad,xq,wq,nsd,1)
		 else if  (nen==8) then
		     call quad3d8n(iquad,nquad,xq,wq,nsd,1)
		 end if

		 if(nen.eq.4) then
		  sq(0,1,iq) = xq(1,iq)
		  sq(0,2,iq) = xq(2,iq)
		  sq(0,3,iq) = xq(3,iq)
		  sq(0,4,iq) = 1 - xq(1,iq) - xq(2,iq) - xq(3,iq)
        else
		  sq(0,1,iq) = (1 - xq(1,iq))	&
			   * (1 - xq(2,iq)) * (1 - xq(3,iq)) / 8
		  sq(0,2,iq) = (1 + xq(1,iq))	&
			   * (1 - xq(2,iq)) * (1 - xq(3,iq)) / 8
		  sq(0,3,iq) = (1 + xq(1,iq))	&
			   * (1 + xq(2,iq)) * (1 - xq(3,iq)) / 8
		  sq(0,4,iq) = (1 - xq(1,iq))	&
			   * (1 + xq(2,iq)) * (1 - xq(3,iq)) / 8
		  sq(0,5,iq) = (1 - xq(1,iq))	&
			   * (1 - xq(2,iq)) * (1 + xq(3,iq)) / 8
		  sq(0,6,iq) = (1 + xq(1,iq))	&
			   * (1 - xq(2,iq)) * (1 + xq(3,iq)) / 8
		  sq(0,7,iq) = (1 + xq(1,iq))	&
			   * (1 + xq(2,iq)) * (1 + xq(3,iq)) / 8
		  sq(0,8,iq) = (1 - xq(1,iq))	&
			   * (1 + xq(2,iq)) * (1 + xq(3,iq)) / 8
		  sq(1,1,iq) = - (1 - xq(2,iq)) * (1 - xq(3,iq)) / 8
		  sq(1,2,iq) = + (1 - xq(2,iq)) * (1 - xq(3,iq)) / 8
		  sq(1,3,iq) = + (1 + xq(2,iq)) * (1 - xq(3,iq)) / 8
		  sq(1,4,iq) = - (1 + xq(2,iq)) * (1 - xq(3,iq)) / 8
		  sq(1,5,iq) = - (1 - xq(2,iq)) * (1 + xq(3,iq)) / 8
		  sq(1,6,iq) = + (1 - xq(2,iq)) * (1 + xq(3,iq)) / 8
		  sq(1,7,iq) = + (1 + xq(2,iq)) * (1 + xq(3,iq)) / 8
		  sq(1,8,iq) = - (1 + xq(2,iq)) * (1 + xq(3,iq)) / 8
		  sq(2,1,iq) = - (1 - xq(1,iq)) * (1 - xq(3,iq)) / 8
		  sq(2,2,iq) = - (1 + xq(1,iq)) * (1 - xq(3,iq)) / 8
		  sq(2,3,iq) = + (1 + xq(1,iq)) * (1 - xq(3,iq)) / 8
		  sq(2,4,iq) = + (1 - xq(1,iq)) * (1 - xq(3,iq)) / 8
		  sq(2,5,iq) = - (1 - xq(1,iq)) * (1 + xq(3,iq)) / 8
		  sq(2,6,iq) = - (1 + xq(1,iq)) * (1 + xq(3,iq)) / 8
		  sq(2,7,iq) = + (1 + xq(1,iq)) * (1 + xq(3,iq)) / 8
		  sq(2,8,iq) = + (1 - xq(1,iq)) * (1 + xq(3,iq)) / 8
		  sq(3,1,iq) = - (1 - xq(1,iq)) * (1 - xq(2,iq)) / 8
		  sq(3,2,iq) = - (1 + xq(1,iq)) * (1 - xq(2,iq)) / 8
		  sq(3,3,iq) = - (1 + xq(1,iq)) * (1 + xq(2,iq)) / 8
		  sq(3,4,iq) = - (1 - xq(1,iq)) * (1 + xq(2,iq)) / 8
		  sq(3,5,iq) = + (1 - xq(1,iq)) * (1 - xq(2,iq)) / 8
		  sq(3,6,iq) = + (1 + xq(1,iq)) * (1 - xq(2,iq)) / 8
		  sq(3,7,iq) = + (1 + xq(1,iq)) * (1 + xq(2,iq)) / 8
		  sq(3,8,iq) = + (1 - xq(1,iq)) * (1 + xq(2,iq)) / 8
        endif
	  end if

!==================================================================
! Second step: Set up the left and right handside matrix
      crra(:)=0
	  crrm(:,:)=0
	  s_area(:)=0

! 2.0 calculate deformation matrix at 1 gau point and area per element
      do ine=1,ne_solid
	     do inen=1,nen
		    y(1:nsd,inen)=x_solid(1:nsd,ien(ine,inen)) ! current position
			x(1:nsd,inen)=x_solid_init(1:nsd,ien(ine,inen)) ! initial position
		 end do
		 rs(1:nsd)=xq(1:nsd,1)
		 call r_element(rs)
		 call r_jacob(y,xj,xji,detJ)
		 call r_jacob(x,toxj,toxji,todet)
		 call r_bdpd_curr(xji)
		 call r_bdpd_init(toxji)
		 call r_stoxc(xto,xot,xj,xji,toxj,toxji,toc,ine)
		 e_area(ine)=wq*detJ
		 call determinant(xot,nsd,nsd,detf)
		 rhs(ine)=detf-1
	  end do
       continue

!loop over element
!==========================================================
! NOTE
! Alougth the system of linear equations is derived from per nodes,
! it is more convenient to calculate the terms in both right and left handside
! matrix from element to elment.
! The formation may be tricky and need to be tested. Also, how I break the coefficients
! of this linear system should be write down in the notes.
!==========================================================

	  do ine=1,ne_solid
! 2.1: get N_{j,i}
         do inen=1,nen
		    x(:,inen)=x_solid(1:nsd,ien(ine,inen))
		 end do
      !...  calculate the shape function and the weight at quad point
		if (nsd==2) then
		    if (nen.eq.3) then !calculate shape function at quad point
			   include "sh2d3n.h"
			elseif (nen.eq.4) then
				include "sh2d4n.h"
			endif
		elseif (nsd==3) then
		    if (nen.eq.4) then !calculate shape function at quad point
			   include "sh3d4n.h"
			elseif (nen.eq.8) then
				include "sh3d8n.h"
			endif
		endif
! 2.2: loop over nsd and elements to set up the left and rigt handside matrix
		do isd=1,nsd
		do inen=1,nen
		   index=ien(ine,inen)+(isd-1)*nn_solid
		   do jsd=1,nsd
		   do jnen=1,nen
		   indey=ien(ine,jnen)+(jsd-1)*nn_solid
           crrm(index,indey)=crrm(index,indey)+sh(isd,inen)*sh(jsd,jnen)
! left handside matrix but missing 1/A_{\alpha} for the diagonal terms

		   end do
		   end do
		   crra(index)=crra(index)+rhs(ine)*sh(isd,inen)

		end do
		end do
		end do
		do ine=1,ne_solid
		   do inen=1,nen
! right handside matrix
		   s_area(ien(ine,inen))=s_area(ien(ine,inen))+e_area(ine)
! s_area is the area sum of the sounding elements
           end do
		end do

          continue
        do inn=1,nn_solid
		   do isd=1,nsd
		      index=inn+(isd-1)*nn_solid
			  crrm(index,index)=crrm(index,index)+1.0/s_area(inn)
! adding 1/A_{\alpha} to the diagonal terms
		   end do
		end do
! Finish setting up left and right handside matrix
!====================================================================
! Step 3 solve the linear equations & update the correction
! For testing just use gausian
!        open (unit=100, file='a.txt', status='unknown')
!        open (unit=101, file='b.txt', status='unknown')
!        do ine=1,nsd*nn_solid
!           do inen=1,nsd*nn_solid
!               write(100,*) crrm(ine,inen)
!           end do
!           write(101,*) crra(ine)
!        end do
        np=nsd*nn_solid
        mp=1
        continue
!=======================================================
!We have two options to solve linear system eqautions
! Gausian
!        call gaussj(crrm,np,np,crra,np,mp)
! CG
         error=0
         crra_cg(:)=0.0d0
         call cg_method(np,crrm,crra,crra_cg,error)
		do isd=1,nsd
		   do inn=1,nn_solid
		      index=inn+(isd-1)*nn_solid
		      !g(isd,inn)=crra(index) !For guasian
                      g(isd,inn)=crra_cg(index) !For CG
			  x_solid(isd,inn)=x_solid(isd,inn)+g(isd,inn)
			  end do
		end do
		return
		end

max_stars_repo_name: biofluids/IFEM-archive

max_issues_repo_name: biofluids/IFEM-archive

max_forks_repo_name: biofluids/IFEM-archive


text:
./tb.v
../mapped/pipeline_riscv_mapped.v
../rtl/instmem.v
../rtl/datamem.v
../lib/saed90nm.v


text:
c******************************************************************************
c
c      Routines for reading old unformatted files. This must be compiled
c      with special compiler switches.
c
c      For Digital Visual Fortran, compiler settings for this file should be
c      for PowerPC 4.0 i/o format and library compatibility, unformatted file
c      conversion = none.
c
c******************************************************************************
c
       SUBROUTINE OUFOVV(IU,FNAME,IERR)
c
c      Opens unformatted file FNAME on unit IU. Returns IERR=0 if OK, else 1.
c------------------------------------------------------------------------------
       CHARACTER*32  FNAME
c------------------------------------------------------------------------------
       OPEN (IU,FILE=FNAME,FORM='UNFORMATTED',ERR=10)
       IERR = 0
       RETURN
10     IERR = 1
       RETURN
       END
c******************************************************************************
c
       SUBROUTINE RSEOVV(IUF,IERR)
c
c      Reads V1 or V1.5 Microsoft compiled .ESS file opened on unit IUF
c      and loads for V2.
c------------------------------------------------------------------------------
c    designators
       CHARACTER*48    NEGM,NEOP,NEMP,NAMB,NICC,NECC,NFLM,NTBM,NVVC,
     ;                 NPCC,NPDS,NFGM,
     ;                 NEGR,NFGT,NFPM,NHTM,NIMM,NEMM,NVRA,NRUN,
     ;                 NXXX,NOLDS,UNSET,NGPM,NIVP,NEVP,NPMP
       CHARACTER*1     NGPMX(48),STRING(80),CU
       CHARACTER*32    NAMEF
c------------------------------------------------------------------------------
c    cylinder work arrays
       DIMENSION   IW(50),RW(400)
       COMMON /IWVV/ IW
       COMMON /RWVV/ RW
c    old work array
       INTEGER*2   IW2(50)
c-----------------------------------------------------------------------------
       EQUIVALENCE (KIVP,IW(18))
       EQUIVALENCE (KEVP,IW(19))
c    Note: in V1 Microsoft-compiled the valve program type 2 was user, but
c          in V2 type 3 is user.
c-----------------------------------------------------------------------------
c    run setup
       COMMON /NRUNVV/ NRUN
c    engine geometry model
       COMMON /NEGMVV/ NEGM
c    operating parameters
       COMMON /NEOPVV/ NEOP
c    model parameters
       COMMON /NEMPVV/ NEMP
c    ambient conditions
       COMMON /NAMBVV/ NAMB
c    intake valve cosine-flat-cosine
       COMMON /NICCVV /NICC
c    exhaust valve cosine-flat-cosine
       COMMON /NECCVV /NECC
c    valve reference areas
       COMMON /NVRAVV /NVRA
c    valve control
       COMMON /NVVCVV /NVVC
c    EGR
       COMMON /NEGRVV /NEGR
c    piston/cylinder conventional crankshaft
       COMMON /NPCCVV/ NPCC
c    piston/cylinder dual stroke
       COMMON /NPDSVV/ NPDS
c    flow model
       COMMON /NFLMVV /NFLM
c    turbulence model
       COMMON /NTBMVV /NTBM
c    flame geometry table
       COMMON /NFGTVV /NFGT
c    flame propagation model
       COMMON /NFPMVV /NFPM
c    gas properties model
       COMMON /NGPMVV /NGPM
c    heat transfer model
       COMMON /NHTMVV /NHTM
c    intake manifold model
       COMMON /NIMMVV /NIMM
c    exhaust manifold model
       COMMON /NEMMVV /NEMM
c    inlet valve profile name
       COMMON /NIVPVV/NIVP
c    exhaust valve profile name
       COMMON /NEVPVV/NEVP
c-----------------------------------------------------------------------------
c    user valve programs
       DIMENSION IVDI(48),FVDI(48)
       DIMENSION IVDE(48),FVDE(48)
       COMMON  /UVIDVV/ NDI,FMXI,IVDI,FVDI
       COMMON  /UVEDVV/ NDE,FMXE,IVDE,FVDE
c-----------------------------------------------------------------------------
       DATA  NOLDS/'Loaded from old-style setup file                '/
c----------------------------------------------------------------------------
c    unset designators for WCR's interface system
       IU = -22222
       RU = -1.1111111E+11
       CU = CHAR(239)
c
c     construct an unset string
        DO I=1,48
           STRING(I) = CU
           ENDDO
        CALL ASETVV(CU,UNSET)
c
c *   Read the old file
c
c    assume style 1
       KESS = 1
c
c    set work array length
10     IF (KESS.EQ.1)  THEN
c            style 1 with 300 RW fields
               NDRW = 300
               DO I=301,400
                   RW(I) = RU
                   ENDDO
           ELSE
c            style 2 with 400 RW fields
               NDRW = 400
           ENDIF
c
c    read the work arrays
        READ(IUF,END=800,ERR=800) IW2,(RW(I),I=1,NDRW)
c    convert integers
       DO I=1,50
           IW(I) = IW2(I)
           ENDDO
c    read the properties data
       CALL RTPOVV(IUF,IERR)
       IF (IERR.NE.0) GOTO 800
c    read the flame geometry table
       CALL RFGOVV(IUF,IERR)
       IF (IERR.NE.0) GOTO 800
c    read the names of parameters and property sets
       IF (KESS.EQ.1)  THEN
c            style 1
               READ (IUF,END=800,ERR=800)  NEGM,NEOP,
     ;         NEVP,NFGM,NFLM,NFPM,NGPMX,NHTM,NIVP,NPMP,NTBM
c           check for user-specified inlet valve program
              IF (KIVP.EQ.2)  THEN
                  CALL RVPOVV(IUF,NIVP,NDI,FMXI,IVDI,FVDI,IERR)
                  IF (IERR.NE.0) GOTO 800
c               change code for V2+
                  KIVP = 4
                  ENDIF
c           check for user-specified exhaust valve program
              IF (KEVP.EQ.2)  THEN
                  CALL RVPOVV(IUF,NEVP,NDE,FMAE,IVDE,FVDE,IERR)
                  IF (IERR.NE.0) GOTO 800
c               change code for V2+
                  KEVP = 4
                  ENDIF
c            set defaults for style 1 given options in style 2
               IW(6) = 1
               DO I=17,20
                   IW(I) = 1
                   ENDDO
               RW(217) = 0
c            set null data not in style 1 files
               DO I=24,29
                   IW(I) = IU
                   ENDDO
c            check ignition
               IF ((IW(11).LT.0).OR.(IW(11).GT.710)) THEN
c                set for motoring
                   IW(11) = IU
                   IW(6) = 2
                   ENDIF
c            set for loaded table
               IW(17) = 2
c            set for no setup save
               IW(35) = 2
c            set for loaded flame geometry table not saved
               IW(36) = 2
               IW(37) = 1
               IW(38) = 2
c            set for no manifold
               IW(39) = 1
               IW(40) = 1
c            convert NGPMX
                CALL ACONVV(NGPMX,NGPM)
c            set new designators
                NRUN = NOLDS
                NEMP = NOLDS
                NAMB = NOLDS
                NICC = UNSET
                NECC = UNSET
                NVRA = NOLDS
                NVVC = NOLDS
                NEGR = NOLDS
                NPCC = NOLDS
                NPDS = UNSET
                NFGT = NFGM
                NIMM = UNSET
                NEMM = UNSET
c            convert NGPMX
                CALL ACONVV(NGPMX,NGPM)
           ELSE
c            style 2
               READ (IUF,END=800,ERR=800)  NEGM,NEOP,
     ;         NEVP,NFGM,NFLM,NFPM,NGPM,NHTM,NIVP,NPMP,NTBM,NIMM,NEMM
c            set for no setup save
               IW(35) = 2
c            set for flame geometry table not saved
               IW(36) = 2
c            set new designators
               NRUN = NOLDS
               NEMP = NOLDS
               NAMB = NOLDS
               NICC = UNSET
               NECC = UNSET
               NVRA = NOLDS
               NVVC = NOLDS
               NEGR = NOLDS
               NPCC = NOLDS
               NPDS = UNSET
               NFGT = NFGM
               NIMM = UNSET
               NEMM = UNSET
           ENDIF
c    check ignition
       IF ((IW(11).LT.0).OR.(IW(11).GT.719)) THEN
c        set for motoring
           IW(11) = IU
           IW(6) = 2
           ENDIF
c    check old designators
       CALL NFIXVV(NEGM)
       CALL NFIXVV(NEOP)
       CALL NFIXVV(NGPM)
       CALL NFIXVV(NHTM)
       CALL NFIXVV(NTBM)
       CALL NFIXVV(NFLM)
       CALL NFIXVV(NFPM)

c    wrap it up
       CLOSE(IUF)
       IERR = 0
       RETURN
c
c   file read failure
800    IF (KESS.EQ.1) THEN
c        try style 2
           KESS = 2
           REWIND(IUF)
           GOTO 10
           ENDIF
c    unresolveable error
       CLOSE(IUF)
       IERR = 1
       RETURN
c
       END
c******************************************************************************
c
       SUBROUTINE OMUFVV(NAME,IU,IERR)
c
c      Opens Microsoft-compiled unformatted file NAME on unit IU
c------------------------------------------------------------------------------
       CHARACTER*32    NAME
c------------------------------------------------------------------------------
       OPEN(IU,FILE=NAME,FORM='UNFORMATTED',ERR=90)
         IERR=0
         RETURN
c    opening error
90     IERR = 1
       RETURN
       END
c******************************************************************************
c
       SUBROUTINE  RTPOVV(IU,IERR)
c
c      Reads old V1 thermodynamic property file opened on unit IU. File left
c      open. Sets up for interpolations.
c
c      Returns IERR = 0 if ok, 1 if no file read.
c------------------------------------------------------------------------------
       CHARACTER*8     CHEM
       CHARACTER*48    NGPM
c------------------------------------------------------------------------------
c    Array dimensions
       PARAMETER (IMAX = 48)
       PARAMETER (NSMAX = 50)
c------------------------------------------------------------------------------
       DIMENSION   CHEM(NSMAX),RMOLS(NSMAX)
       DIMENSION   HP(IMAX),HR(IMAX),PVP(IMAX),PVR(IMAX)
c------------------------------------------------------------------------------
c    Properties file data
       COMMON  /NGPMVV/ NGPM
       COMMON  /PROCVV/ CHEM
       COMMON  /PRONVV/ NF,NO,NP
       COMMON  /PRODVV/ HP,HR,PVP,PVR
       COMMON  /PROPVV/ PP,PR,RFT,RMOLS
c------------------------------------------------------------------------------
       INTEGER*2        NF2,NO2,NP2
c------------------------------------------------------------------------------
c   read the file
      READ (IU,ERR=20,END=20)  NGPM,RFT,NF2,NO2,PR
c   check name
       CALL NFIXVV(NGPM)
c    convert integers
       NF = NF2
       NO = NO2
       NR = NF + NO
       READ (IU,ERR=20,END=20)  (CHEM(J),J=1,NR),(RMOLS(J),J=1,NR)
       READ (IU,ERR=20,END=20)  (HR(J),J=1,IMAX),(PVR(J),J=1,IMAX)
       READ (IU,ERR=20,END=20)  NP2,PP
       NP = NP2
       J1 = NR + 1
       J2 = NR + NP
       READ (IU,ERR=20,END=20)  (CHEM(J),J=J1,J2),
     ;             (HP(J),J=1,IMAX),(PVP(J),J=1,IMAX)
       IERR = 0
       RETURN
c
c    error reading file
20     IERR = 1
       RETURN
       END
c******************************************************************************
c
       SUBROUTINE  RFGOVV(IU,IERR)
c
c      Reads old V1 flame geometry table file opened on unit IU. Leaves open.
c      Returns IERR = 0 read, 1 if not.
c------------------------------------------------------------------------------
       CHARACTER*48    NFGT
       DIMENSION       FAHB(11),RFAP(11)
c------------------------------------------------------------------------------
c    Flame geometry table
       COMMON /FGMTVV/ FAHB,RFAP
       COMMON /NFGTVV/ NFGT
c------------------------------------------------------------------------------
c    read the file
       READ (IU,ERR=20,END=20) NFGT,(FAHB(I),I=1,11),(RFAP(I),I=1,11)
       CALL NFIXVV(NFGT)
       IERR = 0
       RETURN
c    error
20     IERR = 1
       RETURN
       END
c******************************************************************************
c
       SUBROUTINE  RVPOVV(IU,NUVP,ND,FMAX,IVD,FVD,IERR)
c
c      Reads an old-style valve program from opened unit IU.
c      Returns IERR = 0 ok, 1 if not.
c------------------------------------------------------------------------------
c      Nomeclature:
c          NUVP    user valve program name
c          ND      number of data points
c          FMAX    maximum value of relative lift
c          IVD     list of specified angles from valv open angle, deg (1,ND)
c          FVD     list of corresponding lifts (arbitrary units)
c------------------------------------------------------------------------------
       CHARACTER*48    NUVP
       INTEGER*2       ND2,IVD2(48)
       DIMENSION       FVD(48),IVD(48)
c------------------------------------------------------------------------------
c    read the file
       READ (IU,ERR=20,END=20) NUVP,ND2,FMAX
       ND = ND2
       READ (IU,ERR=20,END=20) (IVD2(I),I=1,ND),(FVD(I),I=1,ND)
       CALL NFIXVV(NUVP)
c    convert integers
       DO I=1,ND
           IVD(I) = IVD2(I)
           ENDDO
       IERR = 0
       RETURN
c    error
20     IERR = 1
       RETURN
       END
c******************************************************************************
c
       SUBROUTINE NFIXVV(NAME)
c
c      Replaces blank designator
c------------------------------------------------------------------------------
       CHARACTER*48    NAME,BLANK,NONE
       DATA BLANK/'                                                '/
       DATA NONE /'Undesignated data from old style ESP file       '/
c------------------------------------------------------------------------------
       IF (NAME.EQ.BLANK) NAME = NONE
       RETURN
       END
c******************************************************************************

max_stars_repo_name: kjfergu/esp

max_issues_repo_name: fergu/esp

max_forks_repo_name: fergu/esp


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

!    Date:       7 Januari 1994
!    Program:    BLOUTC.FOR
!    Version:    0.1
!    Programmer: Jos van Gils
!
!    Set output controls for BLOOM
!
!    Called by: BLOOMC
!    Calls    : -

      SUBROUTINE BLOUTC (HISTOR,LPRINO,LDUMPO)
!
!     Arguments
!
!     Name    Type  Length   I/O  Description
!
!     HISTOR  L     1        I    Flag to activate output
!     LPRINO  I     1        I    Saves original value of LPRINT
!     LDUMPO  I     1        I    Saves original value of IDUMP

      LOGICAL         HISTOR
      INTEGER         LPRINO, LDUMPO

!     Common block variables used
!
!     Name    Type  Length   I/O  Inc-file  Description
!
!     LPRINT  I     1        I    sumout    Print flag
!     IDUMP   I     1        I    phyt2     Print flag

      INCLUDE 'blmdim.inc'
      INCLUDE 'phyt2.inc'
      INCLUDE 'sumout.inc'
!
!     Local variables
!
!     Name    Type  Length   I/O  Description
!
!     -
!
!     LPRINT and IDUMP are output control flags of Bloom
!     They have been read from the input and their value has been
!     saved in LPRINO and LDUMPU
!     Here we set them to their original values, ONLY if HISTOR is
!     true, that is for history elements at history times.
!     This is to avoid excessively sized output files of Bloom

      LPRINT = 0
      IDUMP  = 0
      IF (HISTOR) THEN
          LPRINT = LPRINO
          IDUMP  = LDUMPO
      ENDIF

      RETURN
      END

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


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

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

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

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

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

! End MIT license text.

      SUBROUTINE GET_ELEM_AGRID_BGRID ( INT_ELEM_ID, CHECK_AGRID )

! Gets element actual and internal grid numbers given the element's internal ID

      USE PENTIUM_II_KIND, ONLY       :  BYTE, LONG
      USE IOUNT1, ONLY                :  WRT_ERR, WRT_LOG, ERR, F04, F06
      USE SCONTR, ONLY                :  BLNK_SUB_NAM, FATAL_ERR, medat0_cuserin, MELGP, NGRID
      USE TIMDAT, ONLY                :  TSEC
      USE SUBR_BEGEND_LEVELS, ONLY    :  GET_ELEM_AGRID_BGRID_BEGEND
      USE MODEL_STUF, ONLY            :  AGRID, BGRID, EDAT, EID, ELGP, EPNT, ETYPE, GRID, GRID_ID, TYPE

      USE GET_ELEM_AGRID_BGRID_USE_IFs

      IMPLICIT NONE

      CHARACTER(LEN=LEN(BLNK_SUB_NAM)):: SUBR_NAME = 'GET_ELEM_AGRID_BGRID'
      CHARACTER(LEN=*), INTENT(IN)    :: CHECK_AGRID       ! If 'Y' perform check on AGRID's to see if appropriate type

      INTEGER(LONG), INTENT(IN)       :: INT_ELEM_ID       ! Internal element ID for which
      INTEGER(LONG)                   :: IERR        = 0   ! Local error count for BGRID not defined
      INTEGER(LONG)                   :: EPNTK             ! Value from array EPNT at the row for this internal elem ID. It is the
!                                                            row number in array EDAT where data begins for this element.
      INTEGER(LONG)                   :: GRID_ID_ROW_NUM   ! Row num in GRID_ID where AGRID(I) exists
      INTEGER(LONG)                   :: I                 ! DO loop index
      INTEGER(LONG)                   :: DELTA             ! Offset in EDAT (from 1st record for an elem) where grid no's begin
      INTEGER(LONG), PARAMETER        :: SUBR_BEGEND = GET_ELEM_AGRID_BGRID_BEGEND

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

! **********************************************************************************************************************************
      EPNTK = EPNT(INT_ELEM_ID)
      TYPE  = ETYPE(INT_ELEM_ID)
      EID   = EDAT(EPNTK)

! AGRID/BGRID contain the G.P. no's (actual/internal) for points that the elem connects to (not for the v vector)

      DO I=1,MELGP+1
         AGRID(I) = 0
         BGRID(I) = 0
      ENDDO

      CALL GET_ELGP ( INT_ELEM_ID )

      DO I=1,ELGP
         DELTA = 1
         IF (TYPE == 'USERIN  ') THEN
            DELTA = MEDAT0_CUSERIN - 1
         ENDIF
         AGRID(I) = EDAT(EPNTK+I+DELTA)
         CALL GET_ARRAY_ROW_NUM ( 'GRID_ID', SUBR_NAME, NGRID, GRID_ID, AGRID(I), BGRID(I) )
         IF (BGRID(I) == -1) THEN
            WRITE(ERR,1900) AGRID(I), EID, TYPE
            WRITE(F06,1900) AGRID(I), EID, TYPE
            IERR = IERR + 1
            FATAL_ERR = FATAL_ERR + 1
            CYCLE
         ENDIF
      ENDDO

! Test to determine if AGRID's are appropriate for the elem TYPE (do test only if grid exists)

      IF (CHECK_AGRID == 'Y') THEN
         IF ((TYPE(1:4) /= 'ELAS') .AND. (TYPE /= 'USERIN  ')) THEN
            DO I=1,ELGP
               CALL GET_ARRAY_ROW_NUM ( 'GRID_ID', SUBR_NAME, NGRID, GRID_ID, AGRID(I), GRID_ID_ROW_NUM )
               IF (GRID_ID_ROW_NUM > 0) THEN
                  IF (GRID(GRID_ID_ROW_NUM,6) /= 6) THEN
                     IERR = IERR + 1
                     FATAL_ERR = FATAL_ERR + 1
                     WRITE(ERR,1951) TYPE, EID, AGRID(I)
                     WRITE(F06,1951) TYPE, EID, AGRID(I)
                  ENDIF
               ENDIF
            ENDDO
         ENDIF
      ENDIF

      IF (IERR > 0) THEN
         CALL OUTA_HERE ( 'Y' )
      ENDIF

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

      RETURN

! **********************************************************************************************************************************
 1900 FORMAT(' *ERROR  1900: GRID ',I8,' ON ELEMENT ',I8,' TYPE ',A,' NOT  DEFINED')

 1951 FORMAT(' *ERROR  1951: ',A,I8,' USES GRID ',I8,' WHICH IS A SCALAR POINT. SCALAR POINTS NOT ALLOWED FOR THIS ELEM TYPE')

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

      END SUBROUTINE GET_ELEM_AGRID_BGRID

max_stars_repo_name: luzpaz/MYSTRAN

max_issues_repo_name: luzpaz/MYSTRAN

max_forks_repo_name: luzpaz/MYSTRAN


text:
      module rrlw_ref

      use parkind, only : im => kind_im, rb => kind_rb

      implicit none
      save

!------------------------------------------------------------------
! rrtmg_lw reference atmosphere
! Based on standard mid-latitude summer profile
!
! Initial version:  JJMorcrette, ECMWF, jul1998
! Revised: MJIacono, AER, jun2006
! Revised: MJIacono, AER, aug2008
!------------------------------------------------------------------

!  name     type     purpose
! -----  :  ----   : ----------------------------------------------
! pref   :  real   : Reference pressure levels
! preflog:  real   : Reference pressure levels, ln(pref)
! tref   :  real   : Reference temperature levels for MLS profile
! chi_mls:  real   :
!------------------------------------------------------------------

      real(kind=rb) , dimension(59) :: pref
      real(kind=rb) , dimension(59) :: preflog
      real(kind=rb) , dimension(59) :: tref
      real(kind=rb) :: chi_mls(7,59)

      end module rrlw_ref

max_stars_repo_name: nfeldl/climlab

max_issues_repo_name: nfeldl/climlab

max_forks_repo_name: nfeldl/climlab


text:
c-----------------------------------------------------------------------
      subroutine prtrlt(v,m)
      implicit real*8(a-h,o-z)
c
c     ----- print out the lower triangle of a symmetric matrix (stored
c           in packed canonical form (actually an upper triangle) !) -----
c
      include 'stdunit.h'
      dimension v(m*(m+1)/2)

      max=5
      imax = 0
      do while(imax.lt.m)
        imin = imax+1
        imax = min(imax+max,m)
        write(lulog,'(/,5x,10(6x,i4,5x)/)') (i,i = imin,imax)
        do i=1,m
          ii = i*(i-1)/2
          mm = imin + ii
          kk = min(i,imax) + ii
          if(mm.le.kk) then
            write(lulog,'(i4,1x,10e15.7)') i,(v(j),j=mm,kk)
          end if
        end do
      end do
      write(*,*)
      return
      end

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public


text:
The Attendance Office is a lil house on L Street that hosts some pretty serious events.

In 2009, the Attendance Office was referred to as the The Hospice of French Ruud Solace but was later retired due to trying too hard to be clever.

The houses current name comes from a sign from one of the occupants high school attendance offices that now hangs over the door to the backyard. Image(Picture 8.jpg, Ash Reiter and her band., right, thumbnail, 350)Image(Picture 4.jpg, right, 350, thumbnail, The Shants.)

The Attendance Office has hosted, among other shows, the following:

• March 25, 2011

http://www.facebook.com/pages/Youth/136910316350802?skinfo Youth
http://cole.bandcamp.com/album/pyt Cole & the Rats
http://www.myspace.com/gloomprarie Pregnant

• October 10th, 2009
http://www.myspace.com/heycapybara Capybara
http://www.myspace.com/djpoppetsadventures Poppet

• January 10th, 2010
http://www.myspace.com/youryellowdress Your Yellow Dress
http://www.myspace.com/brotherbird All Blood
http://www.myspace.com/nrfitzgerald Nellie Fitzgerald
http://www.myspace.com/neilruud Neil Ruud

• February 27, 2010 (show relocated from The Greenhouse)
http://www.youtube.com/watch?vmfd4LsSsewY MountShout
http://www.youtube.com/watch?vtvBmSraGCY The Shants
http://www.youtube.com/watch?vhnKwza3S36o&featurerelated Ash Reiter

• May 1, 2010
Lions, Bats, and Bears
http://www.youtube.com/watch?vQPU3BnAkGjc&featurerelated Retrospect
http://www.myspace.com/nrfitzgerald Nellie Fitzgerald
R&S

The Attendance Office has also hosted numerous themed parties to breathe some life into Californias Central Valley (e.g. themed Lady Gaga, Iron Chef, etc.).

The house is also part of the KDVS house show network, though maybe not as active as its brothers and sisters: Robot Rocket Residence, Dam Haus, Villanova House, Bomb Shelter.

20110318 14:53:50 nbsp Are there currently any plans for future events? Users/WesP