!--------------------------------------------------------------------------------------------------
! PROGRAM: cis_overlap_prog
!> @author Marin Sapunar, Ruđer Bošković Institute
!> @date August, 2018
!
!> @brief Calculate overlap between two sets of CIS type wave functions.
!--------------------------------------------------------------------------------------------------
program cis_overlap_prog
    ! General
    use global_defs
    ! Variables
    use overlap_variables
    ! I/O
    use overlap_output_mod
    use overlap_input_mod
    ! Calculation
    use overlap_ao_mod
    use overlap_mo_mod
    use overlap_cis_mod
    implicit none

    ! All program variables are defined in overlap_variables_mod

    ! Begin run.
    time00 = omp_get_wtime()
    call set_defaults()
    call command_line_interface()
    if (print_level >= 1) then
        write(stdout, '(5x,a)') '-------------------------------------------------------------'
        write(stdout, '(5x,a)') '                     cis_overlap program                     '
        write(stdout, '(5x,a)') '                         version 1.1                         '
        write(stdout, '(5x,a)') '                                                             '
        write(stdout, '(5x,a)') ' Program compiled on '//__DATE__//' '//__TIME__//'.          '
        write(stdout, '(5x,a)') '-------------------------------------------------------------'
        write(stdout, *)
        write(stdout, '(1x,a,i0,a)') 'Using ', omp_get_max_threads(), ' threads.'
        write(stdout, *)
        write(stdout, '(1x,a)') 'Input for bra orbitals/states: '
        write(stdout, '(5x,a,a)') 'Format: ', input_format_1
        write(stdout, '(5x,a,a)') 'Path: ', path1
        write(stdout, '(1x,a)') 'Input for ket orbitals/states: '
        write(stdout, '(5x,a,a)') 'Format: ', input_format_2
        write(stdout, '(5x,a,a)') 'Path: ', path2
    end if

    ! Calculate AO overlaps.
    call overlap_ao()
    call output_ao(outfile_ao, s_ao)
    if (ao_stop) goto 999

    ! Calculate MO overlaps.
    call overlap_mo()
    call output_mo(outfile_mo, s_mo_a, s_mo_b)
    if (mo_stop) goto 999

    ! Calculate CIS overlaps/dyson orbitals.
    call overlap_cis()
    if (dyson_c) then
        call output_dyson()
    else
        call output_orth(orth_overlap, s_wf)
        call output_phase(match_phase, s_wf)
        call output_wf(outfile_wf, s_wf)
    end if

    ! Timings
    999 time_tot = omp_get_wtime() - time00
    if (print_level >= 2) then
        write(stdout, *)
        write(stdout,'(1x,a)') 'Program time:'
        write(stdout, '(5x, a40, f14.4)') 'Input                     - time (sec):', time_in
        write(stdout, '(5x, a40, f14.4)') 'AO overlap                - time (sec):', time_ao
        write(stdout, '(5x, a40, f14.4)') 'MO overlap                - time (sec):', time_mo
        write(stdout, '(5x, a40, f14.4)') 'WF overlap                - time (sec):', time_wf
        write(stdout, '(5x, 40x, a14)') '--------------'
        write(stdout, '(5x, a40, f14.4)') 'Total                     - time (sec):', time_tot
    end if

    if (print_level >= 1) then
        write(stdout, *)
        write(stdout, '(1x,a)') 'cis_overlap done                                             '
        write(stdout, '(1x,a)') '-------------------------------------------------------------'
    end if

end program cis_overlap_prog


text:
      subroutine carrow(u,v,dx,dy,ndx,nx,ny,
     &                  incx,incy,xedge,yedge,vrefm,vrefp)
      implicit none
c
      integer   ndx,nx,ny,incx,incy
      real      u(ndx,ny),v(ndx,ny), dx(ndx,ny),dy(ndx,ny),
     &          xedge,yedge, vrefm,vrefp
c
c**********
c*
c 1)  a gks based routine for ploting curved vector arrows on a
c      curvi-linear grid in logical array space.
c
c     the arrow head and the length of the arrow tail indicate the
c      current direction and current speed, respectively, at the
c      point marked by the head of the arrow.  the arrow tail
c      follows a streamline.  current speeds larger than 'vrefm'
c      are truncated to 'vrefm'.
c
c 2)  arrows are drawn using the gks polyline routine, 'gpl', with
c      the current atributes for linetype, linecolor, and linewidth.
c
c 3)  argument list:
c
c        u      - x-ward component of vectors to be plotted
c        v      - y-ward component of vectors to be plotted
c        dx     - x-ward grid spacing, see (4)
c        dy     - y-ward grid spacing, see (4)
c        ndx    - 1st dimension of u and v (.ge.nx)
c        nx,ny  - plot ( u(1:nx,1:ny), v(1:nx,1:ny) ), see (5)
c        incx   - plot every incx point in 1st dimension
c        incy   - plot every incy point in 2nd dimension
c        xedge  - minimum x edge (typically 1.0)
c        yedge  - minimum y edge (typically 1.0)
c        vrefm  - reference vector magnitude,   see (5)
c        vrefp  - reference vector plot length, see (5)
c
c        all arguments are unchanged on exit.
c
c 4)  the arrays 'dx' and 'dy' define the grid aspect ratio, with
c      dy(i,j)/dx(i,j) representing the ratio at node (i,j), i.e.
c      the distance between (i,j-eps) and (i,j+eps) divided by
c      the distance between (i-eps,j) and (i+eps,j).
c
c 5)  the vector (u(i,j),v(i,j)) will be plotted with the arrow head
c      at world coordinate (float(i),float(j)).  the current world
c      coordinate transformation will be used, so set the world
c      coordinates before calling 'carrow'.
c     a vector of magnitude 'vm' will have a plotted tail length,
c      in  normalized device coordinates, of 'vrefp' * 'vm'/'vrefm'.
c     all vectors larger in magnitude than 'vrefm' are plotted as
c      if they were 'vrefm' long.
c     the tail will be curved to indicate the path of a particle
c      that ends up at the head, assuming the current field does
c      not change with time, i.e. it will trace a streamline.  the
c      path will allow for changes in the grid aspect ratio, as
c      defined by 'dx' and 'dy'.
c     since 'vrefp' is in normalized device coordinates it must
c      have a value between 0.0 and 1.0, typically 0.01.
c     all plotted arrow heads are the same size.
c     no vector is plotted if it is outside xedge:nx-xedge,yedge:ny-yedge.
c
c 6)  nothing is plotted at (float(i),float(j)) whenever
c      both u(i,j) and v(i,j) are zero.  This should always
c      be the case over "land".
c     note that current values are typically required at all
c      "sea" nodes, not just as plotted nodes.
c     The u and v velocity are co-located, but if they originally
c      came from a C-grid the land boundary should be half way
c      between grid points.  For this reason, the tail will be
c      terminated if the nearest grid point contains zero speed.
c
c 8)  alan j. wallcraft,  march 1990 and may 1992.
c     updated for curvi-linear grids and HYCOM november 2008.
c*
c**********
c
      external gpl,gqcntn,gqnt
c
      integer    mxlen
      parameter (mxlen=500)
c
      logical ldebug
      integer i,ii,j,jj,k,lbest,len,len1,lentmp,lenx, ierr,ntrn
      real    di,dj,scale,sp,spp,up,uw,vp,vw,
     +        xh,xo,xq,xscale,xyscal,xw,yh,yo,yq,yscale,yw
      real    xp(mxlen+7),yp(mxlen+7), window(4),viewpt(4)
c
c     st = sin of 22.5 degrees, and ct = cos of 22.5 degrees.
c
      real       st,ct
      parameter (st=0.382683432365090, ct=0.923879532511287)
c
c     head = size of arrow head in normalized device coordinates,
c
      real       head
      parameter (head=0.006)
c
      real       zero,eps,p1,one
      parameter (zero=0.0, eps=1.e-5, p1=0.1, one=1.0)
c
c     conversion factor between n.d.c. and world coordinates.
c
      call gqcntn(ierr, ntrn)
      call gqnt(  ntrn, ierr, window, viewpt )
      xscale = (window(2) - window(1)) / (viewpt(2) - viewpt(1))
      yscale = (window(4) - window(3)) / (viewpt(4) - viewpt(3))
      xyscal = max( xscale, yscale )
c
c     loop through all nodes, drawing arrows if required.
c
      ldebug = .false. !default
c
      do 11 j=max(incy,2),ny-1,incy
        do 12 i=max(incx,2),nx-1,incx
*         ldebug = mod(i,nx/3).lt.incx .and. mod(j,ny/3).lt.incy
c
          if     (u(i,j).ne.zero .or. v(i,j).ne.zero) then
c
c           arrow in normalized device coordinates.
c
            sp  = vrefp * min( one,
     +                         sqrt(u(i,j)**2 + v(i,j)**2) / vrefm )
            up  = dy(i,j)/dx(i,j) * u(i,j)/xscale
            vp  =                   v(i,j)/yscale
            spp = sqrt(up**2 + vp**2)
            up  = up*sp/spp
            vp  = vp*sp/spp
                if     (ldebug) then
                  write(6,*) 'i,j =',i,j
                  write(6,*) 'x,yscale =',xscale,yscale
                  write(6,*) 'sp,up,vp =',sp,up,vp
                endif !ldebug
c
c           define the arrow head (in world coordinates).
c
            xh = xscale*p1*head*up/sp
            yh = yscale*p1*head*vp/sp
c
            xp(1) = i
            yp(1) = j
            xp(2) = i - xscale*head*(ct*up+st*vp)/sp
            yp(2) = j - yscale*head*(ct*vp-st*up)/sp
            xp(3) = xp(2) - xh
            yp(3) = yp(2) - yh
            xp(4) = xp(1) - xh
            yp(4) = yp(1) - yh
            xp(5) = i - xscale*head*(ct*up-st*vp)/sp - xh
            yp(5) = j - yscale*head*(ct*vp+st*up)/sp - yh
            xp(6) = xp(5) + xh
            yp(6) = yp(5) + yh
            xp(7) = i
            yp(7) = j
c
c           find the end point of the 'best' curved tail.
c
            xq    = i
            yq    = j
            lbest = max( 5, min( mxlen,
     +                           nint(sp/0.001),         !small on the plot
     +                           nint(2.0*sp*xscale),    !< 1/2 grid spacing
     +                           nint(2.0*sp*yscale) ) ) !< 1/2 grid spacing
            scale = sp/lbest
                if     (ldebug) then
                  write(6,*) '*****'
                  write(6,*) 'lbest,scale =',lbest,scale
                  write(6,*) 'k,xq,yq =',0,xq,yq
                endif !ldebug
            do 21 k= 1,lbest
              xo = xq !old value
              yo = yq !old value
              ii = min( nx-1, int(xq) )
              jj = min( ny-1, int(yq) )
              di = xq - ii
              dj = yq - jj
              xw =  (one-di)*(one-dj)*dx(ii,  jj)   +
     +              (one-di)*     dj *dx(ii,  jj+1) +
     +                   di *(one-dj)*dx(ii+1,jj)   +
     +                   di *     dj *dx(ii+1,jj+1)
              yw =  (one-di)*(one-dj)*dy(ii,  jj)   +
     +              (one-di)*     dj *dy(ii,  jj+1) +
     +                   di *(one-dj)*dy(ii+1,jj)   +
     +                   di *     dj *dy(ii+1,jj+1)
              uw = yw/xw*
     +             ( (one-di)*(one-dj)*u(ii,  jj)   +
     +               (one-di)*     dj *u(ii,  jj+1) +
     +                    di *(one-dj)*u(ii+1,jj)   +
     +                    di *     dj *u(ii+1,jj+1)  ) / xscale
              vw = ( (one-di)*(one-dj)*v(ii,  jj)   +
     +               (one-di)*     dj *v(ii,  jj+1) +
     +                    di *(one-dj)*v(ii+1,jj)   +
     +                    di *     dj *v(ii+1,jj+1)  ) / yscale
              spp = sqrt( uw**2 + vw**2 )
                  if     (ldebug) then
                    write(6,*) 'k,ii,jj =',k,ii,jj
                    write(6,*) 'k,di,dj =',k,di,dj
                    write(6,*) 'k,yw/xw,uw,vw =',k,yw/xw,uw,vw
                    write(6,*) 'k,spp =',k,spp
                  endif !ldebug
              if     (spp.gt.eps*vrefm/xyscal) then  !non-zero velocity
c               direction is from velocity in WC, i.e. in array index space
c               length is scale in NDC
                xq = xq - (xscale*uw)*(scale/spp)
                yq = yq - (yscale*vw)*(scale/spp)
                    if     (ldebug) then
                      write(6,*) 'k,xq,yq =',k,xq,yq
                    endif !ldebug
                if     (xq.ge.nx .or. xq.lt.1 .or.
     +                  yq.ge.ny .or. yq.lt.1     ) then
                  xq = xo
                  yq = yo
                      if     (ldebug) then
                        write(6,*) 'k,xo,yo =',k,xq,yq
                      endif !ldebug
                  exit !21
                elseif (u(nint(xq),nint(yq)).eq.zero .and.
     +                  v(nint(xq),nint(yq)).eq.zero       ) then
                  xq = xo
                  yq = yo
                      if     (ldebug) then
                        write(6,*) 'k,xo,yo =',k,xq,yq
                      endif !ldebug
                  exit !21
                endif
              else  !zero velocity
                exit !21
              endif
   21       continue !k
c
c           find the least numbr of points on the curve such that
c            the differnce from the best curve is too small to notice.
c
            len1 = max( 2, min( 5, lbest/3 ) )
                if     (ldebug) then
                  write(6,*) '*****'
                  write(6,*) 'lbest,len1  =',lbest,len1
                  write(6,*) 'k,xq,yq =',0,xp(7),yp(7)
                endif !ldebug
            len  = len1
            do 31 lentmp= len1,lbest
                  if     (ldebug) then
                    write(6,*) '*****'
                    write(6,*) 'lbest,lentmp  =',lbest,lentmp
                    write(6,*) 'k,xq,yq =',0,xp(7),yp(7)
                  endif !ldebug
              lenx  = 7+len
              scale = sp/len
              do 32 k= 8,7+len
                ii = min( nx-1, int(xp(k-1)) )
                jj = min( ny-1, int(yp(k-1)) )
                di = xp(k-1) - ii
                dj = yp(k-1) - jj
                xw =  (one-di)*(one-dj)*dx(ii,  jj)   +
     +                (one-di)*     dj *dx(ii,  jj+1) +
     +                     di *(one-dj)*dx(ii+1,jj)   +
     +                     di *     dj *dx(ii+1,jj+1)
                yw =  (one-di)*(one-dj)*dy(ii,  jj)   +
     +                (one-di)*     dj *dy(ii,  jj+1) +
     +                     di *(one-dj)*dy(ii+1,jj)   +
     +                     di *     dj *dy(ii+1,jj+1)
                uw = yw/xw*
     +               ( (one-di)*(one-dj)*u(ii,  jj)   +
     +                 (one-di)*     dj *u(ii,  jj+1) +
     +                      di *(one-dj)*u(ii+1,jj)   +
     +                      di *     dj *u(ii+1,jj+1)  ) / xscale
                vw = ( (one-di)*(one-dj)*v(ii,  jj)   +
     +                 (one-di)*     dj *v(ii,  jj+1) +
     +                      di *(one-dj)*v(ii+1,jj)   +
     +                      di *     dj *v(ii+1,jj+1)  ) / yscale
                spp = sqrt( uw**2 + vw**2 )
                  if     (ldebug) then
                    write(6,*) 'k,ii,jj =',k-7,ii,jj
                    write(6,*) 'k,di,dj =',k-7,di,dj
                    write(6,*) 'k,yw/xw,uw,vw =',k-7,yw/xw,uw,vw
                    write(6,*) 'k,spp =',k-7,spp
                  endif !ldebug
                if     (spp.gt.eps*vrefm/xyscal) then !non-zero velocity
                  xp(k) = xp(k-1) - (xscale*uw)*(scale/spp)
                  yp(k) = yp(k-1) - (yscale*vw)*(scale/spp)
                      if     (ldebug) then
                        write(6,*) 'k,xq,yq =',k-7,xp(k),yp(k)
                      endif !ldebug
                  if     (xp(k).ge.nx .or. xp(k).lt.1 .or.
     +                    yp(k).ge.ny .or. yp(k).lt.1     ) then
                    lenx = k-1
                    exit !32
                  elseif (u(nint(xp(k)),nint(yp(k))).eq.zero .and.
     +                    v(nint(xp(k)),nint(yp(k))).eq.zero) then
                    lenx = k-1
                    exit !32
                  endif
                else  !zero velocity
                  lenx = k-1
                  exit !32
                endif
   32         continue !k
c
              if     (abs(xq-xp(lenx))/xscale.lt.0.0005 .and.
     +                abs(yq-yp(lenx))/yscale.lt.0.0005      ) then
                exit !31
              elseif (len.eq.lbest) then
                exit !31
              else
                len = min( lbest, len + 1 + len/10 )
              endif
   31       continue !lentmp
c
            if     (minval(xp(1:lenx)).gt.   xedge .and.
     &              maxval(xp(1:lenx)).lt.nx-xedge .and.
     &              minval(yp(1:lenx)).gt.   yedge .and.
     &              maxval(yp(1:lenx)).lt.ny-yedge      ) then
              call gpl(lenx, xp,yp)
            endif !safe to draw vector
          endif !non-zero current
   12   continue !i
   11 continue !j
      return
c     end of carrow.
      end
      subroutine clegend(xc,yc,yl, vrefm,vrefp,cunits)
      implicit none
c
      character*(*) cunits
      real          xc,yc,yl, vrefm,vrefp
c
c**********
c*
c 1)  a gks based routine for ploting a key for routine carrow.
c
c     the arrow head and the length of the arrow tail indicate the
c      current direction and current speed, respectively.
c
c 2)  arrows are drawn using the gks polyline routine, 'gpl', with
c      the current atributes for linetype, linecolor, and linewidth.
c
c 3)  argument list:
c
c        xc     - x location of the arrow center in world coordinates
c        yc     - y location of the arrow center in world coordinates
c        yl     - y location of the label center in world coordinates
c        vrefm  - reference vector magnitude,   see carrow
c        vrefp  - reference vector plot length, see carrow
c        cunits - units of the vector (e.g. "cm/s")
c
c        all arguments are unchanged on exit.
c
c 8)  alan j. wallcraft, NRL, november 2008.
c*
c**********
c
c --- fieldproc specific
      integer          nperfr,locbar
      common/perframe/ nperfr,locbar
      save  /perframe/
c
      external gpl,gqcntn,gqnt
c
      character text*6
      integer   ierr,ntrn
      real      csn,xscale,yscale
      real      sp,up,vp,x,xh,y,yh
      real      xp(8),yp(8), window(4),viewpt(4)
c
c     st = sin of 22.5 degrees, and ct = cos of 22.5 degrees.
c
      real       st,ct
      parameter (st=0.382683432365090, ct=0.923879532511287)
c
c     head = size of arrow head in normalized device coordinates,
c
      real       head
      parameter (head=0.006)
c
      real       p1
      parameter (p1=0.1)
c
c     conversion factor between n.d.c. and world coordinates.
c
      call gqcntn(ierr, ntrn)
      call gqnt(  ntrn, ierr, window, viewpt )
      xscale = (window(2) - window(1)) / (viewpt(2) - viewpt(1))
      yscale = (window(4) - window(3)) / (viewpt(4) - viewpt(3))
c
c     arrow in normalized device coordinates.
c
      sp  = vrefp
      up  = sp
      vp  = 0.0
c
c     define the arrow head (in world coordinates).
c
      x   = xc + 0.5*xscale*vrefp
      y   = yc
c
      xh = xscale*p1*head*up/sp
      yh = yscale*p1*head*vp/sp
c
      xp(1) = x
      yp(1) = y
      xp(2) = x - xscale*head*(ct*up+st*vp)/sp
      yp(2) = y - yscale*head*(ct*vp-st*up)/sp
      xp(3) = xp(2) - xh
      yp(3) = yp(2) - yh
      xp(4) = xp(1) - xh
      yp(4) = yp(1) - yh
      xp(5) = x - xscale*head*(ct*up-st*vp)/sp - xh
      yp(5) = y - yscale*head*(ct*vp+st*up)/sp - yh
      xp(6) = xp(5) + xh
      yp(6) = yp(5) + yh
      xp(7) = x
      yp(7) = y
      xp(8) = x - xscale*vrefp
      yp(8) = y
c
      call gpl(8, xp,yp)
c
      if     (nperfr.ne.4) then
        csn=-0.9
      else
        csn=-0.9*0.7
      endif
      write (text,"(f6.2)") vrefm
      if     (text(4:6).eq.'.00') then
         text(4:6) = '   '
      elseif (text(6:6).eq.'0') then
         text(6:6) = ' '
      endif
      call pcloqu(xc,yl,trim(text)//' '//trim(cunits),csn,0.,0.)
      return
c     end of clegend.
      end

ext: f

max_stars_repo_name: TillRasmussen/HYCOM-tools

max_issues_repo_name: TillRasmussen/HYCOM-tools

max_forks_repo_name: TillRasmussen/HYCOM-tools


text:
	SUBROUTINE GH_TC2VG ( device, tnam, wtype, timstr, wname, origc,
     +                       wocen, wsped, wadnm, wpres, sylat,sylon,
     +			     sytim, tau0, rlat, rlon,
     +                       nc, nstrm, nknt, fdate,wwnd, wgst, lext,
     +                       idisp, mxwd, gust,wdir, itrack, szmrk,
     +                       ltype, lwidth, lcolor, reset, iret )
C************************************************************************
C* GH_TC2VG								*
C*									*
C* This subroutine writes the current tropical storm or hurricane       *
C* track, error cone  and watches/warnings to three VG files.		*
C*									*
C* GH_TC2VG ( DEVCIE, TNAM, WTYPE, TIMSTR, WNAME, ORIGC, WOCEN, WSPED,  *
C*	     WADNM, WPRES,RLAT,RLON,NC,NSTRM,NKNT,FDATE,WWND,WGST,      *
C*           LEXT,IDISP, MXWD,WGST,WDIR,ITRACK,SZMRK,LTYPE,LWIDTH,      *
C*	     LCOLOR, RESET, IRET )	         			*
C*									*
C* Input parameters:							*
C*      DEVICE          CHAR*           Current device                  *
C*	TNAM		CHAR*		Output filename			*
C* 	WTYPE    	CHAR*		Tropical storm type		*
C*	TIMSTR    	CHAR*		Advisory valid time string	*
C*	WNAME    	CHAR*		Tropical storm name		*
C*	ORIGC    	CHAR*		Issuing center                  *
C*     	WOCEN    	CHAR*		Tropical storm identifier       *
C*	WSPED    	CHAR*		Movement speed (kts)		*
C*	WADNM    	CHAR*		Storm advisory number		*
C*	WPRES    	CHAR*		Storm center pressure		*
C*	RLAT (NC,*) 	REAL		Current/forecasted Latitudes	*
C*	RLON (NC,*) 	REAL		Current/forecasted Longitudes	*
C*	NC 		INTEGER		Maximum number of lat/lon pairs *
C*	NSTRM 		INTEGER		Number of decoded files		*
C*      NKNT		INTEGER         Number of forecast points	*
C*      FDATE (NC,*)	CHAR*		Forecast dates			*
C*	WWND		CHAR*		Current storm speed             *
C*	WGST		CHAR*		Current storm wind gust         *
C*      LEXT (NC,*)     CHAR*           Forecast Post-Trop indicators   *
C*	IDISP		INTEGER		Storm dissipating flag		*
C*      MXWD (NC,*)	INTEGER		Forecast max wind speeds        *
C*      GUST (NC,*)	INTEGER		Forecast wind gusts             *
C*	WDIR		CHAR*		Direction of present movement   *
C*      ITRACK		INTEGER		Track plot flag			*
C*	SZMRK (2)	REAL		Marker size multiplier		*
C*	LTYPE (2)	INTEGER		Forecast track line type	*
C*	LWIDTH (2)	INTEGER		Forecast track line width	*
C*	LCOLOR (2)	INTEGER		Forecast track line color	*
C*									*
C* Input/Output parameters:						*
C*	RESET		LOGICAL		Color reset flag		*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C*									*
C**									*
C* Log:									*
C* m.gamazaychikov/SAIC	06/07	From GH_TCWW				*
C* m.gamazaychikov/SAIC	10/07	Fixed the misspelling of "Storm"	*
C* m.gamazaychikov/SAIC	10/08	Cleaned up, added stsrc strings to 	*
C*				gh_wtct CS, updated cyint and cylong	*
C* A. Krautkramer/NHC	 5/09	Properly set storm type and time zone	*
C* X. Guo/CWS		03/10	Removed argument rmnt and ilegend  	*
C* A. Krautkramer/NHC	 6/09	Add variable definition, replace the	*
C*				first element of the tau array with the	*
C*				passed in tau				*
C* A. Krautkramer	10/10	Corrected initial point inconsistency	*
C* 				with the calculation of dvlbl and tcdvlp*
C* A. Krautkramer/NHC	09/11	Added 'Post-Tropical Cyclone'		*
C* S. Jacobs/NCEP	 4/13	Added smoothing for the track and cone	*
C* S. Jacobs/NCEP	 9/13	Fixed problem with VG output for track	*
C*				less than 3 days			*
C* M. Sardi/NHC		 7/14	Fixed problem with 3 day cone by	*
C*				reverting first GH_KGAT call.		*
C* S. Jacobs/NCEP	 8/14	Removed previous change - did not work	*
C* B. Hebbard/NCEP	 5/20	Update for new 60hr forecast - SCN20-20	*
C*                              by increasing PARAMETER NUMHR 8 -> 9    *
C*                              and adding '60' value to DATA tau       *
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
	INCLUDE		'ghcmn.cmn'
C*
	PARAMETER	( NUMHR = 9, MXPTS = 500 )
	PARAMETER	( NUMH2 = NUMHR - 2 )
C*
	CHARACTER*(*)	device, wname, timstr, wtype, wsped, tnam,
     +                  wocen, fdate(nc,*), wwnd, wadnm, wdir, origc,
     +                  wpres, wgst, lext (nc,*)
	CHARACTER	sytim*20, ctau*4
	INTEGER		mxwd (nc,*), ltype (*), lwidth (*), lcolor (*),
     +                  gust (nc,*)
        REAL		rlat(nc,*), rlon(nc,*), szmrk (*)
        REAL		sylat, sylon
	LOGICAL		reset
C*
	CHARACTER	table*128, tbtyp*128,
     +                  ddev*12,
     +                  ddate(NUMHR)*20, tzone, filnam*40,
     +                  tdev*30,
     +			adnum*4, cyint (5), bpnam(100)*33, stnum*4,
     +			yyyy*4, stname*32, vtime*10, bpstr*3300,
     +			zone*4, status, cadvnm*5,
     +			fcstpd*10, wbas*4,
     +			wtyp*4, dummy*30, dummy2*30,
     +			cylong(5)*20, yymm*4, tdate*20, yyy*2, mm*2,
     +			fmonth*9, fmon*9, fday*9, ctime*2,
     +                  zone2*4, ampm*2, bklntp*4
	CHARACTER*30	tau (NUMHR), maxwnd(nc), wgust(nc), mslp(nc),
     +			tcdvlp(nc), dvlbl(nc), tcdir(nc), tcspd(nc),
     +			dtlbl(nc), advdat(nc), stsrc (nc)
	CHARACTER*240	advstr, taustr, wndstr, gststr, prsstr,
     +			dvpstr, dvlstr, dirstr, spdstr, dtlstr,
     +			srcstr
	CHARACTER*240	advstr6, taustr6, wndstr6, gststr6, prsstr6,
     +			dvpstr6, dvlstr6, dirstr6, spdstr6, dtlstr6,
     +			srcstr6
C*
	INTEGER		numbkp(50),
     +			numb(4), ibkpts(4,100), iarea(4,100),
     +                  isev(50), iadvtp(50), igeog(50), icnt (4, 50),
     +			jtarr(5),ibkcur(MAXBK), iarcur(MAXBK),
     +			nppart(20), kpos(MXPTS)
	CHARACTER	cmin*2,cyear*5,cmonth*3,cday*3
        INTEGER		tau0, len, leny, lenm, lend, lenmn
C Start Krautkramer 10/10 modification
	INTEGER         wwnd_int (1), num
C End Krautkramer
C*
        REAL            slat(NUMHR), slon(NUMHR), alat(NUMHR),
     +			alon(NUMHR), slnlat(MXPTS), slnlon(MXPTS),
     +			slnlt6(MXPTS), slnln6(MXPTS), dist(NUMHR),
     +			bplat(100), bplon(100),
     +			lbklat(20,MAXBK), lbklon(20,MAXBK),
     +			bkplat(MAXBK), bkplon(MAXBK),
     +			blat(MXPTS), blon(MXPTS), bdst(MXPTS),
     +			bx(MXPTS), by(MXPTS), dx(NUMHR), dy(NUMHR)
C*
	LOGICAL		hvbkpt, finish, plotit, first, extra
C*
	INCLUDE		'ERMISS.FNC'
C*
	DATA		cyint / 'M','H', 'S', 'D', 'L'/
	DATA		cylong / 'Major Hurricane',
     +                           'Hurricane',
     +				 'Tropical Storm',
     +				 'Tropical Depression',
     +                           'Remnant Low' /
	DATA		tau / '0','12','24','36','48',
     +                        '60','72','96','120'/
C-----------------------------------------------------------------------
	iret = 0
C
C*	Query current settings.
C
        CALL GH_SAVE ( ier )
C
	finish = .false.
	extra  = .false.
        ddev = 'XW'
        filnam = ' '
        yyyy = ' '
        stname = ' '
        stnum = ' '
        wbas = ' '
        vtime = ' '
        bpstr = ' '
        srcstr = ' '
        srcstr6= ' '
        dummy = '9999'
        dummy2 = ' '
        iwht = 31
        iblk = 32
        npline = 0
        DO ibpnam = 1, 100
            bpnam(ibpnam) = ' '
        END DO
        DO ii = 1, MAXBK
            bkplat (ii) = RMISSD
            bkplon (ii) = RMISSD
        END DO
C
C*	convert tau0 from int to char
C
	WRITE  ( ctau, '(I2)', IOSTAT = iostat ) tau0
C
C*	append null character to the string
C
	CALL ST_NULL(ctau,ctau,len,iret)
C
C* 	replace the first tau in the array with the passed in tau0
C
	tau(1) = ctau
C
C*      Set the flag for extratropical line in legend box
C
        ii = 1
        iextra = 0
        DO WHILE ( (.not. extra) .and. (ii .le. nknt) )
           IF ( lext(ii,nstrm) .eq. 'E' ) THEN
              iextra = 1
              extra = .true.
           END IF
           ii = ii + 1
        END DO
        CALL ST_LSTR ( tnam, lenm, ier )
        tdev= ddev // '| TRACK-W/W'
        CALL GG_SDEV ( tdev, ier )
        IF ( .not. reset ) THEN
            CALL GH_SVCL ( ier )
            reset = .true.
          ELSE
            CALL GH_RSCL ( ier )
        END IF
        CALL GSTANM ( ier )
	DO WHILE ( .not. finish )
C
C*          Reverse black and white if device is for postscript.
C
            CALL GQDEV  ( ddev, iunit, iatyp, ier )
            IF  ( ddev(:2) .eq. 'PS' )  THEN
                iwht = 32
                iblk = 1
              ELSE
                iwht = 31
                iblk = 32
            END IF
            IF  ( ddev ( :2) .ne. 'XW' ) THEN
                finish = .true.
            END IF
C
C*	    Set symbol and text attributes.
C
	    CALL GSTEXT ( 21, 2, 1.0, 1, 111, 1, 1, ier )
	    isytyp = 17
            ssize  = 1.0
            iwidth = 2
C
C*          Read in average track error file.
C
	    iertb = 0
            IF ( wocen(1:2)  .eq. 'AL' ) THEN
                table = 'altker.tbl'
                IF ( rlon(1,nstrm) .le. -85.0 ) THEN
                    tzone = 'C'
                  ELSE
                    tzone = 'E'
                END IF
              ELSE IF ( wocen(1:2)  .eq. 'EP' ) THEN
                table = 'eptker.tbl'
                tzone = 'P'
              ELSE IF ( wocen(1:2)  .eq. 'CP' ) THEN
                table = 'cptker.tbl'
                tzone = 'H'
	      ELSE
		iertb = -10
            END IF
            IF ( origc(1:4)  .eq. 'CPHC' ) THEN
                table = 'cptker.tbl'
                tzone = 'H'
            END IF
C
	    IF ( iertb .eq. 0 ) THEN
                tbtyp = 'hcnadv'
C
C*              Read the average track error file.
C
                CALL GH_KGTB ( table, tbtyp, NUMHR, dist, iertb )
	    END IF
	    IF ( iertb .ne. 0 ) THEN
		CALL ER_WMSG ( 'GPTCWW', iertb, ' ', ierr )
	    END IF
C
            im    = nstrm
            nltln = 0
            jmax  = NUMHR
	    IF ( nknt .le. NUMH2 ) jmax = NUMH2
C
            IF ( ERMISS ( rlat(1,im) ) .or. ERMISS ( rlon(1,im) ) ) THEN
                im = im - 1
                jmax = 1
            END IF
C
C*          Find the average track error lat. lon. points.
C
            DO jj = 1, jmax
                IF ( ( .not. ERMISS ( rlat (jj,im) ) ) .and.
     +               ( .not. ERMISS ( rlon (jj,im) ) ) ) THEN
                    nltln = nltln + 1
                    alat (nltln) = rlat(jj,im)
                    alon (nltln) = rlon(jj,im)
		    IF ( jj .gt. 1 ) ddate (nltln) = fdate(jj-1,im)
                END IF
            END DO
C
C*	    Smooth the central track line and the track error distances.
C
	    CALL GTRANS ( 'M', 'D', nltln, alat, alon, dx, dy, iergt )
	    dens = 5.0
	    cscl = 30.0
	    is = 0
	    ie = nltln+1
	    CALL CV_PRM3 ( nltln, dx, dy, dist, dens, MXPTS, cscl,
     +			   is, ie, nout, bx, by, bdst, iercv )
	    CALL GTRANS ( 'D', 'M', nout, bx, by, blat, blon, iergt )
C
C*	    Find the original points in the smoothed line.
C
C	    knumh2 = nout
	    DO  kk = 1, nout
	        kpos(kk) = 0
		DO  nn = 1, nltln
		    IF  ( ( ABS(blat(kk)-alat(nn)) .lt. 0.0005 ) .and.
     +			  ( ABS(blon(kk)-alon(nn)) .lt. 0.0005 ) )  THEN
			kpos(kk) = nn
			IF ( nn .eq. NUMH2 ) knumh2 = kk
		    END IF
		END DO
	    END DO
C
C*	    Calculate the area under the cone.
C
            IF ( iertb .ge. 0 ) THEN
		plotit = .true.
		inner  = MIN0 ( nltln, NUMH2 )
		minpt  = MIN0 ( nout, knumh2 )
                CALL GH_KGAT ( blat, blon, nout, bdst, MXPTS, slnlt6,
     +                         slnln6, numln6, ier )
		IF ( nout .gt. knumh2 .and. knumh2 .ne. 0 ) THEN
                    CALL GH_KGAT ( blat, blon, nout, bdst, MXPTS,
     +                             slnlat, slnlon, numln, ier )
		  ELSE
		    numln = numln6
		    DO jj = 1, numln
			slnlat ( jj ) = slnlt6 ( jj )
			slnlon ( jj ) = slnln6 ( jj )
		    END DO
		END IF
C
C*              Check for breakpoints.
C
                hvbkpt = .true.
C
C*              Get the name of VGF file containing breakpoints.
C
                CALL ST_LSTR ( wocen, lenocn, ier )
                CALL ST_NULL ( wadnm, wadnm, lena, ier )
                CALL GH_ADVN ( wadnm, numadv, iadflag, ier )
                CALL TI_C2I  (timstr, iyyyy, immdd, ihhmm, ier )
                IF ( numadv .gt. 0 ) THEN
                    numadv = numadv + 1000
                    CALL ST_INCH ( numadv, adnum, ier )
                    CALL ST_INCH ( iyyyy, yyyy, ier )
                    CALL ST_UCLC ( wocen (:4), stnum, ier)
                    IF ( iadflag .eq.0 ) THEN
                        filnam = 'tca_' // stnum // yyyy //'_'//
     +                            adnum ( 2:4 ) //'.vgf'
                    ELSE IF ( iadflag .eq.1 ) THEN
                        filnam = 'tca_' // stnum // yyyy //'_'//
     +                            adnum ( 2:4 ) //'a.vgf'
                    ELSE IF ( iadflag .eq.2 ) THEN
                        filnam = 'tca_' // stnum // yyyy //'_'//
     +                            adnum ( 2:4 ) //'b.vgf'
                    ELSE
                        filnam = ' '
                    END IF
                  ELSE
                    filnam = ' '
                END IF
                CALL ST_NULL (filnam, filnam, ilen, iret)
C
C*              Get the breakpoint information from the VGF file
C
                CALL GH_BKRV ( filnam, istnum, iyear, ibasin, cadvnm,
     +                         isstyp, stname, vtime, icount,
     +                         isev,  iadvtp, igeog,
     +                         bplat, bplon, bpstr, numbkp, zone,
     +                         status, iret)
                IF ( iret .eq. -1) THEN
                    CALL ER_WMSG ( 'GPTCWW', 4, ' ', ierr )
                    hvbkpt = .false.
                ELSEIF ( icount .EQ. 0 ) THEN
                    hvbkpt = .false.
                ELSE
C
C*                  Separate name string into an array of strings
C
                    CALL ST_CLSL ( bpstr, ';', ' ', 100,
     +                             bpnam, nbpn, iret)
C
C*                  Get the information to plot breakpoints
C
                    CALL GH_BKGB ( icount, isev, iadvtp, bpnam, numbkp,
     +                             ibkpts, numb, iarea, icnt, iret)
                END IF
                CALL ST_LSTR ( wocen, lenocn, ier )
                CALL ST_UCLC ( wocen (3:4), stnum, ier)
                CALL ST_UCLC ( wocen (:2), wbas, ier)
C
C*	        Plot the average track error line.
C
                IF ( idisp .ge. 0 ) THEN
                  CALL ST_NULL ( tnam, tnam, lena, ier )
                  CALL ST_NULL ( wadnm, wadnm, lena, ier )
                  CALL ST_NULL ( wtype, wtyp, lena, ier )
                  CALL ST_NULL ( timstr,timstr, lena, ier )
                  CALL ST_NULL ( stnum, stnum, lena, ier )
                  CALL ST_NULL ( wname,wname, lena, ier )
                  CALL ST_NULL ( wbas, wbas, lena, ier )
C
C*                3day track error cone.
C
                  fcstpd = '072'
                  iclr = iblk
                  iltyp = 1
                  ifclr = iwht
                  iftyp = 1
                  CALL ST_NULL ( fcstpd, fcstpd, lena, ier )
                  CALL GH_WTCE ( tnam, wtyp, timstr, stnum,
     +                            wname, wbas, wadnm,
     +                            fcstpd,
     +                            ilclr, iltyp, ifclr, iftyp,
     +                            numln6, slnlt6, slnln6,
     +                            iret)
                  IF ( nout .gt. knumh2 .and. knumh2 .ne. 0 ) THEN
C
C*                    5day track error cone.
C
                      fcstpd = '120'
                      iclr = iwht
                      iltyp = 1
                      ifclr = iwht
                      iftyp = 2
                      CALL ST_NULL ( fcstpd, fcstpd, lena, ier )
                      CALL GH_WTCE ( tnam, wtyp, timstr, stnum,
     +                            wname, wbas, wadnm,
     +                            fcstpd,
     +                            ilclr, iltyp, ifclr, iftyp,
     +                            numln, slnlat, slnlon,
     +                            iret)
                  END IF
C
C*	    Plot the forecasted track.
C*	    Set the color and symbol for the current location.
C
                tnam = tnam (:8)//'TCt'//tnam(12:)
                CALL ST_NULL ( tnam, tnam, lena, ier )
                iclr = iblk
                iltyp = 1
C
C*              Lat/Lons and other TC Track parameters

                DO kk = 1, nknt
                    stsrc(kk) = 'Tropical Cyclone'
                    slat(kk) = rlat(kk,im)
                    slon(kk) = rlon(kk,im)
                    IF ( kk .eq. 1 ) THEN
                      advdat(kk) = timstr(5:)
                      maxwnd(kk) = wwnd
                      wgust(kk)  = wgst
                      mslp(kk)   = wpres
                      tcdir(kk)  = wdir
                      tcspd(kk)  = wsped

		      IF (wtype .eq. 'PT') THEN
		           stsrc(kk) = 'Post-Tropical Cyclone'
		      END IF

C START Krautkramer Sept 2, 2010
C Changed due to inconsitency between the calculation of dvlbl and tcdvlp between
C the initial point (k=1) and k > 1

C	Convert the string stored in maxwnd into a integer for comparison

		      CALL ST_C2I (maxwnd(kk), 1, wwnd_int, num, ier)

		      IF (wwnd_int(1) .ge. 96 ) THEN
                         ii = 1
                      ELSE IF ( wwnd_int(1) .ge. 64 ) THEN
                         ii = 2
                      ELSE IF ( wwnd_int(1) .ge. 34 ) THEN
                         ii = 3
                      ELSE
                         ii = 4
                      END IF

		      dvlbl(kk) = cyint (ii)
                      tcdvlp(kk) = cylong (ii)

C                      IF ( wtype (1:1) .eq. 'T') THEN
C                         dvlbl(kk)  = wtype(2:2)
C                         IF ( wtype (2:2) .eq. 'D') THEN
C                            tcdvlp(kk) = "Tropical Depression"
C                          ELSE
C                            tcdvlp(kk) = "Tropical Storm"
C                         END IF
C		      ELSE IF ( wtype (1:1) .eq. 'S') THEN
C		      	 dvlbl(kk) = wtype(2:2)
C			 IF  (wtype (2:2) .eq. 'S') THEN
C			    tcdvlp(kk) = "Tropical Storm"
C			 ELSE
C			    tcdvlp(kk) = "Tropical Depression"
C			 END IF
C                      ELSE
C                         dvlbl(kk)  = wtype (1:1)
C                         tcdvlp(kk) = "Hurricane"
C                      END IF
C END Krautkramer Sept 2, 2010

                    ELSE

                      mslp(kk)  = dummy
                      tcdir(kk) = dummy
                      tcspd(kk) = dummy

                      advdat(kk) = fdate(kk-1, im)
                      CALL ST_INLN (mxwd (kk-1,im),maxwnd(kk),lena,ier)
                      CALL ST_INLN (gust (kk-1,im),wgust(kk), lena,ier)

                      IF ( mxwd ( kk - 1, im ) .ge. 96 ) THEN
                         ii = 1
                      ELSE IF ( mxwd ( kk - 1, im ) .ge. 64 ) THEN
                         ii = 2
                      ELSE IF ( mxwd ( kk - 1, im ) .ge. 34 ) THEN
                         ii = 3
                      ELSE
                         ii = 4
                      END IF

                      dvlbl(kk) = cyint (ii)
                      tcdvlp(kk) = cylong (ii)

                      IF ( lext(kk-1, im) .eq. 'P' ) THEN
                          stsrc(kk) = 'Post-Tropical Cyclone'
                      ELSE
                          stsrc(kk) = 'Tropical Cyclone'
                      ENDIF
                    ENDIF
C
C*                  Create the day/time label string.
C
                    IF ( kk .eq. 1 ) THEN
                        yymm = timstr ( 1:4 )
                        tdate = timstr
                        CALL ST_INTG ( timstr ( 5:6 ), iday, ier )
                        first = .true.
                      ELSE
                        CALL ST_INTG ( fdate (kk-1,im) (1:2), nday,ier)
                        IF ( ( nday .lt. iday ) .and. first ) THEN
                            first = .false.
                            CALL ST_INTG ( timstr(3:4), month, ier )
                            month = month + 1
                            IF ( month .gt. 12 ) THEN
                                month = 1
                                CALL ST_INTG ( timstr(1:2), iyear, ier)
                                iyear = iyear + 1
                                CALL ST_INCH ( iyear, yyy, ier )
                                IF ( iyear .lt. 10 )
     +                               yyy = '0' // yyy ( 1:1 )
                                yymm ( 1:2 ) = yyy
                            END IF
                            CALL ST_INCH ( month, mm, ier )
                            IF ( month .lt. 10 ) mm = '0' // mm ( 1:1 )
                            yymm ( 3:4 ) = mm
                        END IF
                        tdate = yymm // fdate ( kk-1,im )
                    END IF
                    CALL GH_TIME ( tdate, tzone, jtarr, hours, zone2,
     +                             fmonth, fmon, fday, ampm, ier )
		    CALL ST_NULL ( zone2, zone2, lena, ier )
		    CALL ST_INCH ( jtarr (5), cmin, ier )
		    if ( jtarr (5) .eq. 0) THEN
		    	cmin = '00'
		    ENDIF
		    CALL ST_LSTR ( cmin, lenmn, ier )
  		    CALL ST_INCH ( jtarr (4), ctime, ier )
                    CALL ST_LSTR ( ctime, lent, ier )
                    CALL ST_LSTR ( ampm, lena, ier )
		    CALL ST_INCH ( jtarr (1), cyear, ier )
		    CALL ST_LSTR ( cyear, leny, ier )
		    CALL ST_INCH ( jtarr (2), cmonth, ier )
		    IF ( jtarr(2) .lt. 10 ) THEN
		    	cmonth = '0' // cmonth
		    END IF
		    CALL ST_LSTR ( cmonth, lenm, ier )
		    CALL ST_INCH ( jtarr (3), cday, ier )
		    IF ( jtarr(3) .lt. 10 ) THEN
		    	cday = '0' // cday
		    END IF
		    CALL ST_LSTR ( cday, lend, ier )
		    dtlbl(kk)=cyear(:leny) // cmonth(:lenm)
     +   	    // cday(:lend) // ' ' // ctime(:lent)
     +              // cmin(:lenmn) // ' ' // ampm (:lena) //
     +		    ' ' // fday (1:3)
                END DO
                CALL ST_LSTC  ( advdat, inner, '|', advstr6, ier )
                CALL ST_LSTC  ( tau,    inner, '|', taustr6, ier )
                CALL ST_LSTC  ( maxwnd, inner, '|', wndstr6, ier )
                CALL ST_LSTC  ( wgust,  inner, '|', gststr6, ier )
                CALL ST_LSTC  ( mslp,   inner, '|', prsstr6, ier )
                CALL ST_LSTC  ( tcdvlp, inner, '|', dvpstr6, ier )
                CALL ST_LSTC  ( dvlbl,  inner, '|', dvlstr6, ier )
                CALL ST_LSTC  ( tcdir,  inner, '|', dirstr6, ier )
                CALL ST_LSTC  ( tcspd,  inner, '|', spdstr6, ier )
                CALL ST_LSTC  ( dtlbl,  inner, '|', dtlstr6, ier )
                CALL ST_LSTC  ( stsrc,  inner, '|', srcstr6, ier )
                CALL ST_NULL ( advstr6, advstr6,lena, ier )
                CALL ST_NULL ( taustr6, taustr6,lena, ier )
                CALL ST_NULL ( wndstr6, wndstr6,lena, ier )
                CALL ST_NULL ( gststr6, gststr6,lena, ier )
                CALL ST_NULL ( prsstr6, prsstr6,lena, ier )
                CALL ST_NULL ( dvpstr6, dvpstr6,lena, ier )
                CALL ST_NULL ( dvlstr6, dvlstr6,lena, ier )
                CALL ST_NULL ( dirstr6, dirstr6,lena, ier )
                CALL ST_NULL ( spdstr6, spdstr6,lena, ier )
                CALL ST_NULL ( dtlstr6, dtlstr6,lena, ier )
                CALL ST_NULL ( srcstr6, srcstr6,lena, ier )
C
C*              3 day track.
C
                 fcstpd = '072'
                 iclr = iblk
                 iltyp = 1
                 CALL ST_NULL ( fcstpd, fcstpd, lena, ier )
		 CALL GH_WTCT ( tnam, zone2, wtyp, timstr, stnum,
     +                         wname, wbas, wadnm,
     +                         fcstpd,
     +                         ilclr, iltyp, inner,
     +                         slat, slon,
     +                         advstr6, taustr6, wndstr6,
     +                         gststr6, prsstr6,
     +                         dvpstr6, dvlstr6, dirstr6,
     +                         spdstr6, dtlstr6, srcstr6,
     +                         iret)

                 IF ( nout .gt. knumh2 .and. knumh2 .ne. 0 ) THEN
C
C*                  5day track.
C
                    fcstpd = '120'
                    iclr = iblk
                    iltyp = 3
                    CALL ST_LSTC  ( advdat, nknt, '|', advstr, ier )
                    CALL ST_LSTC  ( tau,    nknt, '|', taustr, ier )
                    CALL ST_LSTC  ( maxwnd, nknt, '|', wndstr, ier )
                    CALL ST_LSTC  ( wgust,  nknt, '|', gststr, ier )
                    CALL ST_LSTC  ( mslp,   nknt, '|', prsstr, ier )
                    CALL ST_LSTC  ( tcdvlp, nknt, '|', dvpstr, ier )
                    CALL ST_LSTC  ( dvlbl,  nknt, '|', dvlstr, ier )
                    CALL ST_LSTC  ( tcdir,  nknt, '|', dirstr, ier )
                    CALL ST_LSTC  ( tcspd,  nknt, '|', spdstr, ier )
                    CALL ST_LSTC  ( dtlbl,  nknt, '|', dtlstr, ier )
                    CALL ST_LSTC  ( stsrc,  nknt, '|', srcstr, ier )
                    CALL ST_NULL ( advstr, advstr,lena, ier )
                    CALL ST_NULL ( taustr, taustr,lena, ier )
                    CALL ST_NULL ( wndstr, wndstr,lena, ier )
                    CALL ST_NULL ( gststr, gststr,lena, ier )
                    CALL ST_NULL ( prsstr, prsstr,lena, ier )
                    CALL ST_NULL ( dvpstr, dvpstr,lena, ier )
                    CALL ST_NULL ( dvlstr, dvlstr,lena, ier )
                    CALL ST_NULL ( dirstr, dirstr,lena, ier )
                    CALL ST_NULL ( spdstr, spdstr,lena, ier )
                    CALL ST_NULL ( dtlstr, dtlstr,lena, ier )
                    CALL ST_NULL ( fcstpd, fcstpd, lena, ier )
                    CALL ST_NULL ( srcstr, srcstr, lena, ier )
                    CALL GH_WTCT ( tnam, zone2, wtyp, timstr, stnum,
     +                           wname, wbas, wadnm,
     +                           fcstpd,
     +                           ilclr, iltyp, nknt, slat, slon,
     +                           advstr, taustr, wndstr, gststr, prsstr,
     +                           dvpstr, dvlstr, dirstr, spdstr, dtlstr,
     +                           srcstr, iret)
                END IF
C
           IF ( hvbkpt ) THEN
C
C*         Breakpoint output
C
                tnam = tnam (:8)//'TCb'//tnam(12:)
                CALL ST_NULL ( tnam, tnam, lena, ier )
                    DO ii = 1, 4
                      IF ( numb (ii) .gt. 0 ) THEN
                        IF ( ii .eq. 1 ) THEN
                           itcww = 0
                           ilclr = 2
                           iltyp = 12
                        END IF
                        IF ( ii .eq. 2 ) THEN
                           itcww = 1
                           ilclr = 11
                           iltyp = 5
                        END IF
                        IF ( ii .eq. 3 ) THEN
                           itcww = 2
                           ilclr = 4
                           iltyp = 8
                        END IF
                        IF ( ii .eq. 4 ) THEN
                           itcww = 3
                           ilclr = 5
                           iltyp = 2
                        END IF
                         indbkc = 1
                         indbkp = 1
C
C*                       Start loop over the number of segments
C
                         DO ilines = 1, numb(ii)
                          inpseg = icnt (ii,indbkc)
C
C*                        Zero out the brakpoint and area arrays
C*                        for current segment
C
                          DO ipts = 1, inpseg
                             ibkcur(ipts) = 0
                             iarcur(ipts) = 0
                          END DO
                          indb2=1
C
C*                        Start loop over the segments number of points
C
                          DO ipts = 1, icnt (ii,indbkc)
                             ibkcur(indb2) = ibkpts(ii, indbkp)
                             iarcur(indb2) = iarea(ii, indbkp)
                             IF ( iarcur(indb2) .eq. IUSGEC .or.
     +                            iarcur(indb2) .eq. IMXCSA .or.
     +                            iarcur(indb2) .eq. IPACIF .or.
     +                            iarcur(indb2) .eq. IKEYS ) THEN
                                 bklntp = "pair"
                              ELSE IF ( iarcur(indb2) .eq. ICUBA .or.
     +                                  iarcur(indb2) .eq. IHISP.or.
     +                                  iarcur(indb2) .eq. IPRICO ) THEN
                                 bklntp = "ordr"
                              ELSE IF ( iarcur(indb2) .eq. IOTHER .or.
     +                                  iarcur(indb2) .eq. IWATER ) THEN
                                 bklntp = "list"
                                 CALL GH_GBKL ( ibkcur(indb2),
     +                                          iarcur(indb2),
     +                                          1, bklntp,
     +				                nparts, nppart,
     +				                lbklat, lbklon,
     +					        iret)
                                 fcstpd = '120'
                                 CALL ST_NULL ( tnam, tnam, lena, ier )
                                 CALL ST_NULL ( wadnm, wadnm, lena,ier)
                                 CALL ST_NULL ( wtype, wtyp, lena, ier)
                                 CALL ST_NULL ( timstr,timstr,lena,ier)
                                 CALL ST_NULL ( stnum, stnum, lena,ier)
                                 CALL ST_NULL ( wname,wname, lena, ier)
                                 CALL ST_NULL ( wbas, wbas, lena, ier )
                                 CALL ST_NULL ( fcstpd,fcstpd,lena,ier)
                                 DO inparts = 1, nparts
                                   npline = nppart(inparts)
                                   DO ibk= 1, npline
                                     bkplat(ibk) = lbklat(inparts,ibk)
                                     bkplon(ibk) = lbklon(inparts,ibk)
                                   END DO
                                   CALL GH_WTCB ( tnam, wtyp, timstr,
     +                                         stnum, wname,wbas,wadnm,
     +                                         fcstpd, ilclr, iltyp,
     +                                         itcww,
     +                                         npline,
     +                                         bkplat, bkplon,
     +                                         iret)
                                END DO
                             END IF
                          indb2 = indb2+1
                          indbkp = indbkp + 1
C
C*                        End of POINTS loop
C
                          END DO

                            IF  ( bklntp .eq. "ordr" .or.
     +                            bklntp .eq. "pair") THEN
                                  CALL GH_GBKL ( ibkcur, iarcur,
     +				                 inpseg, bklntp,
     +				                 nparts, nppart,
     +				  		 lbklat, lbklon,
     +					         iret)
                              fcstpd = '120'
                              CALL ST_NULL ( tnam, tnam, lena, ier )
                              CALL ST_NULL ( wadnm, wadnm, lena, ier )
                              CALL ST_NULL ( wtype, wtyp, lena, ier )
                              CALL ST_NULL ( timstr,timstr, lena, ier )
                              CALL ST_NULL ( stnum, stnum, lena, ier )
                              CALL ST_NULL ( wname,wname, lena, ier )
                              CALL ST_NULL ( wbas, wbas, lena, ier )
                              CALL ST_NULL ( fcstpd, fcstpd, lena,ier)
                              npline = nppart(1)
                              DO ibk= 1, npline
                                     bkplat(ibk) = lbklat(1,ibk)
                                     bkplon(ibk) = lbklon(1,ibk)
                              END DO
                              CALL GH_WTCB ( tnam, wtyp, timstr,stnum,
     +                                     wname, wbas, wadnm,
     +                                     fcstpd, ilclr, iltyp,
     +                                     itcww,
     +		                           npline, bkplat, bkplon,
     +                                     iret)
                            END IF

                            DO ipts = 1, inpseg
                              ibkcur(ipts) = 0
                              iarcur(ipts) = 0
                            END DO
                            indbkc = indbkc + 1
C
C*                       End of SEGMENTS loop
C
                         END DO
                      END IF
C
C*                  End of NUMB loop
C
                    END DO
                END IF
		plotit = .false.
		finish = .true.
                END IF
	      ELSE
		plotit = .false.
		hvbkpt = .false.
		finish = .true.
            END IF
C
C*	    Flush the graphics buffer.
C
	    CALL GEPLOT ( ier )
            CALL GENANM ( ier )
C
C*          Check for an acceptable display.
C
            IF ( .not. finish ) THEN
                IF ( device(:2) .eq. 'XW' ) THEN
                    finish = .true.
                ELSE
                    ddev = device
                    CALL ST_LSTR(tnam, lenm, ier )
                    CALL ST_LSTR(ddev, lend, ier )
                    device = ddev(:lend) // '|'
     +                    // tnam(:lenm)
                    CALL GG_SDEV ( device, ier )
                    CALL GSTANM ( ier )
                END IF
            END IF
	END DO
C
C*	Reset the saved attributes.
C
	CALL GH_REST ( ier )
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
Fire hydrants are used by the Fire Department and the UC Davis Fire Department. Davis has a total of 1,699 hydrants (IN TOWN, UCD has more). This is a photographic catalog of them.  In case your fingers and toes dont count that high, it doesnt catalog all 1,699 hydrants.  Yet.

 City Hydrants

The City of Davis uses dry hydrants.  Dry hydrants do not have any water in them when they are not in use.  These hydrants are filled with water by turning the special bolt at the top of the hydrant which opens a valve, allowing the hydrant to pressurize from water in the main below.  This feature protects the hydrants from breaking in freezing temperatures.  This feature, however, makes use of the hydrants cumbersome in fire fighting operations because all outlets of the hydrants are charged at once.

 Freestanding

 Attached to Buildings

 Not attached to anything

The hydrant on the roof of the Chi Phi house.

 Campus Hydrants

Unlike the City of Davis, the UC Davis university uses wet hydrants.  The outlets in wet hydrants open individually, expanding their usability in fighting fire. Wet hydrants can freeze and break in long periods of cold weather. Sometimes contractors hit them and generate spontaneous fountains.

 Geysers

 Freestanding

The header image for this page is of Hydrant H20.

 Attached to Buildings

See also: Dogs

This page totally reminds me of the Wiki/April Fools Day, 2005 Davis Manhole Project. Users/BrentLaabs

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
c     ---------------------------------------------------------------
      program fronts
c     ---------------------------------------------------------------

c     Atmospheric Dynamics Group IAC, ETH Zurich

c     Computation of location of fronts on THE charts
c     Advanced algorithm for front locating function used
c     Output with frontal strength

c     For advection speed: Output with negative sign for warm fronts
c     and information of velocity and direction of motion

c     Optional output for frontal areas and removal of small objects
c     Optional removal of quasi-stationary fronts by min advection speeds

c     Literature:
c     -----------
c     Code is partly based on ideas of:
c     Hewson TD. 1998. Objective fronts. Met. Apps 5: 37–65, doi:10.1017/ S1350482798000553.

c     First development and application of our tool to COSMO Reanalysis:
c     Jenkner J. et al. 2010. Detection and climatology of fronts in a high-resolution model reanalysis over the alps. Met. Apps 17: 1–18, doi:10.1002/met.142.

c     Introduction of frontal areas, mobile fronts, and first global climatology by our tool:
c     Schemm S. et al. 2014. Extratropical Fronts in the Lower Troposphere – Global Perspectives obtained from two Automated Methods. Quart. J. Roy. Meteor. Soc.

c     History:
c     ---------
c     Johannes Jenkner April 2006
c                      September 2006
c                      September 2008
c     Sebastian Schemm October 2013
c     		           June 2014

c     Development of front2line tool
c     Michael Sprenger April 2014

      implicit none

c     ---------------------------------------------------------------
c     Declaration of variables
c     ---------------------------------------------------------------

      character(len=100) :: innam,outnam
      character(len=80)  :: invar,cfn

      integer :: iargc,ierr,grid
      integer :: idin,idout,idcst
      integer :: ndim,ntimes,nmin
      integer :: i,j,k,n,nz,nlev,ntot

      integer, dimension(4) :: vardim
      integer, dimension(5) :: stdate

      real :: thld
      real :: misdat
      real :: pollon,pollat
      real :: dx,dy
      real :: phx,phy
      real :: minadv

      real,dimension(    4) :: varmin,varmax,stag
      real,dimension(10000) :: time

      real,dimension(:),allocatable :: aklev,bklev,aklay,bklay

      real,dimension(:,:),allocatable :: clusterin

      real,dimension(:,:,:),allocatable :: gradx,grady,vert,help,
     >                                     help1,help2,condi1,
     >                                     condi2,tfield,farea

c     MSP
      integer,dimension(:,:),allocatable :: belowc,abovec
      integer                               nbelow,nabove
      integer                               belowl(5000)
      integer                               il,ir,ju,jd
      integer                               ind,ival

      integer,dimension(:,:),allocatable :: clusterout,onecluster

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

      logical :: winds,plus,withgap,area,mobile,objects

c      common /pole/ pollon,pollat

c     -----------------------------------------------
c     Check number of input arguments
c     -----------------------------------------------

      print*,""
      print*," FRONT DETECTION "
      print*,""

      if ( (iargc().lt.4).or.(iargc().gt.10) ) then
        write(*,*) 'USAGE: ./fronts [-winds][-plus][-withgap]'//
     >   '[-area][-mobile thld][-objects minpixel]'//
     >	 ' INFILE OUTFILE VARIABLE THRESHOLD '

   	  print*," "
	  write(*,*) 'Optional parameters:'
	  write(*,*) ' -winds: output of advection velocity and direction'
	  write(*,*) ' -plus: enlarge the frontal lines'
	  write(*,*) ' -withgap: gaps between mobile sections'
	  write(*,*) ' -area: output frontal areas'
	  write(*,*) ' -mobile threshold: remove quasi-stationary fronts'//
     >		     ' by minimum advection criterion in m/s'
	  write(*,*) ' -objects minpixel: remove small objects with size'//
     >  	     ' lower or equal of "minpixel" grid points'

	  print*, " "
	  stop

      end if

c     -----------------------------------------------
c     Read input parameters and input file
c     -----------------------------------------------

      i=1
      call getarg(i,invar)
      winds=(invar=='-winds')
      if (winds) i=i+1

      call getarg(i,invar)
      plus=(invar=='-plus')
      if (plus) i=i+1

      call getarg(i,invar)
      withgap=(invar=='-withgap')
      if (withgap) i=i+1

      call getarg(i,invar)
      area=(invar=='-area')
      if (area) i=i+1

      call getarg(i,invar)
      mobile=(invar=='-mobile')
      if (mobile) then
       i=i+1
        call getarg(i,invar)
        read(invar,*) minadv
        i=i+1
      endif

      call getarg(i,invar)
      objects=(invar=='-objects')
      if (objects) then
       i=i+1
        call getarg(i,invar)
        read(invar,*) nmin
        i=i+1
      endif

      call getarg(i,innam)
      call getarg(i+1,outnam)
      call getarg(i+3,invar)
      read(invar,*) thld
      call getarg(i+2,invar)

c     -----------------------------------------------
c     Print Info
c     -----------------------------------------------

      if (winds)   write(*,"(A)") "-winds option selected"
      if (plus)    write(*,"(A)") "-plus option selected"
      if (withgap) write(*,"(A)") "-withgap option selected"
      if (area)    write(*,"(A)") "-area option selected"
      if (mobile) write(*,"(A,F4.1)") "-mobile option selected: ",minadv
      if (objects) write(*,"(A,I4.1)") "-objects option selected:",nmin

c     -----------------------------------------------
c     Access to input netcdf file
c     -----------------------------------------------

      call cdfopn(trim(innam),idin,ierr)
      call getcfn(idin,cfn,ierr)
      call getdef(idin,trim(invar),
     >            ndim,misdat,vardim,varmin,varmax,stag,ierr)
      nz=vardim(3)

c     -----------------------------------------------
c     Allocate required arrays
c     -----------------------------------------------

      allocate(help(vardim(1),vardim(2),nz))
      allocate(gradx(vardim(1),vardim(2),nz))
      allocate(grady(vardim(1),vardim(2),nz))
      allocate(help1(vardim(1),vardim(2),nz))
      allocate(help2(vardim(1),vardim(2),nz))
      allocate(vert(vardim(1),vardim(2),nz))
      allocate(condi1(vardim(1),vardim(2),nz))
      allocate(condi2(vardim(1),vardim(2),nz))
      allocate(tfield(vardim(1),vardim(2),nz))
      allocate(farea(vardim(1),vardim(2),nz))
      allocate(clusterout(vardim(1),vardim(2)))
      allocate(clusterin(vardim(1),vardim(2)))
      allocate(onecluster(vardim(1),vardim(2)))

c     MSP
      allocate(belowc(vardim(1),vardim(2)))
      allocate(abovec(vardim(1),vardim(2)))

c     -----------------------------------------------
c     Read in definitions of input netcdf file
c     -----------------------------------------------

      call cdfopn(cfn,idcst,ierr)

c     Get grid dx / dy
      call getgrid(idcst,dx,dy,ierr)
      call clscdf(idcst,ierr)

c     Fill level reference array
      do i=1,nz
        vert(:,:,i)=real(i)
      enddo

c     Determine the x topology of the grid
c     2: periodic, cyclic point;
c     1: periodic, not cyclic;
c     0: not periodic (and therefore not closed)

c     dx   = (xmax-xmin) / real(nx-1)

      if (abs(varmax(1)-varmin(2)-360.).lt.eps) then
         grid=2
      elseif (abs(varmax(1)-varmin(2)-360.+dx).lt.eps) then
         grid=1
      else
         grid=0
      endif

c     Get time axis
      call gettimes(idin,time,ntimes,ierr)

c     -----------------------------------------------
c     Create output netcdf file
c     -----------------------------------------------

      call crecdf(trim(outnam),idout,varmin,varmax,3,
     >            trim(cfn),ierr)

c     Define output variables
      call putdef(idout,'fronts_'//trim(invar),ndim,misdat,
     >            vardim,varmin,varmax,stag,ierr)

      if (winds) then
        call putdef(idout,'fvel_'//trim(invar),ndim,
     >              misdat,vardim,varmin,varmax,stag,ierr)
        call putdef(idout,'fdir_'//trim(invar),ndim,
     >              misdat,vardim,varmin,varmax,stag,ierr)
      endif

      if (area) then
        call putdef(idout,'farea_'//trim(invar),ndim,
     >              misdat,vardim,varmin,varmax,stag,ierr)
      endif

      if (mobile) then
        call putdef(idout,'fmobile_'//trim(invar),ndim,
     >              misdat,vardim,varmin,varmax,stag,ierr)
      endif

      call putdef(idout,invar,ndim,
     >              misdat,vardim,varmin,varmax,stag,ierr)

c     -----------------------------------------------
c     Loop over times
c     -----------------------------------------------

      do n=1,ntimes

        call getdat(idin,trim(invar),time(n),0,help,ierr)
        tfield = help

c     -----------------------------------------------
c       Calculation of gradients and derived fields
c     -----------------------------------------------

c	    1) Gradient of THE

        call grad(gradx,grady,help,vert,varmin(1),varmin(2),
     >            dx,dy,vardim(1),vardim(2),nz,misdat)

c       2) Abs. value of gradient of THE (condi1)

        call abs2d(gradx,grady,condi1,vardim(1),vardim(2),
     >            nz,misdat)

c       3) Gradient of abs.value of gradient of THE (help1/help2)

        call grad(help1,help2,condi1,vert,varmin(1),varmin(2),
     >            dx,dy,vardim(1),vardim(2),nz,misdat)

c 	    4) Scalar product of 3) and 1) -> Gradient in direction of 1)

        call scal2d(help1,help2,gradx,grady,condi2,vardim(1),
     >            vardim(2),nz,misdat)

c       The TFP is condi2, not normalized by abs. value of 1)

c     -------------------------------------------------------
c       Calculate advection speed in direction of grad(TFP)
c     -------------------------------------------------------

        if (winds) then

          where (abs(condi2-misdat).gt.eps) condi2=-condi2

c        1) Gradient of TFP
          call grad(gradx,grady,condi2,vert,varmin(1),varmin(2),
     >            dx,dy,vardim(1),vardim(2),nz,misdat)

c        2) Normalized by abs. value of 1)
          call norm2d(gradx,grady,vardim(1),vardim(2),nz,misdat)

        endif

c       --------------------------------------------------------------
c       Calculation of masking condition: abs(grad TFP).gt. threshold
c       --------------------------------------------------------------

        where ((abs(condi1-misdat).gt.eps).and.(condi1.gt.thld))
          help=1
        elsewhere
          help=0
        endwhere

c       --------------------------------------------------------------
c       Save this area as frontal area
c       --------------------------------------------------------------

        where ((abs(condi1-misdat).gt.eps).and.(condi1.gt.thld))
          farea=1
        elsewhere
          farea=0
        endwhere

c       --------------------------------------------------------------
c       Calculation of locating condition; zero TFP
c       --------------------------------------------------------------

        call zcontour(condi2,vardim(1),vardim(2),nz,
     >                varmin(1),varmin(2),dx,dy,misdat)

        help = help*condi2

c       --------------------------------------------------------------
c       Avoid detection of minimum gradient
c       --------------------------------------------------------------

        call div(condi2,help1,help2,vert,varmin(1),
     >           varmin(2),dx,dy,vardim(1),vardim(2),nz,misdat)

        where (condi2.gt.0) help=0
        where (condi2.gt.0) farea=0

c       -----------------------------------------------
c       Apply special emphasizing filter for lines
c       -----------------------------------------------
        if (plus) then

         call linesplus(help,vardim(1),vardim(2),nz,
     >                   varmin(1),varmin(2),dx,dy,misdat)
        endif

c       --------------------------------------------------------------
c       Scale output with frontal intensity
c       --------------------------------------------------------------

        help  = help *condi1
        farea = farea*condi1

c       --------------------------------------------------------------
c       Set vacant pixels to missing
c       --------------------------------------------------------------

        where (help ==0)  help = misdat
        where (farea==0) farea = misdat

c       --------------------------------------------------------------
c       Further specification of frontal characteristics by winds
c       --------------------------------------------------------------

        if (winds.or.mobile) then

          call getdat(idin,'U',time(n),0,help1,ierr)
          call getdat(idin,'V',time(n),0,help2,ierr)

c       --------------------------------------------------------------
c       Advection speed in direction of normalized TFP
c       --------------------------------------------------------------

          call scal2d(help1,help2,gradx,grady,condi1,vardim(1),
     >              vardim(2),nz,misdat)

c         Rotate wind vector if necessary (COSMO only)
c          if (abs(pollat-90.).gt.eps) then
c            do k=1,nz
c              do j=1,vardim(2)
c                phy=varmin(2)+(j-1)*dy
c                do i=1,vardim(1)
c                  phx=varmin(1)+(i-1)*dx
c                  call t_phuv2lluv(gradx(i,j,k),grady(i,j,k),
c     >                             help1(i,j,k),help2(i,j,k),phx,phy)
c                enddo
c              enddo
c            enddo
c          else
            help1=gradx
            help2=grady
c          endif

c       --------------------------------------------------------------
c        Correction of frontal direction for warm fronts
c       --------------------------------------------------------------

          where ((abs(condi1-misdat).gt.eps).and.(condi1.lt.0.))
            help1 = -help1
            help2 = -help2
          end where

          call winddir(vardim(1),vardim(2),nz,misdat,help1,help2,
     >                 condi2)

          where (abs(help-misdat).lt.eps)
            condi1 = misdat
            condi2 = misdat
          end where

c       --------------------------------------------------------------
c       end further specifications in terms of winds
c       --------------------------------------------------------------

      endif

c     --------------------------------------------------------------
c     Remove quasi-stationary fronts
c     --------------------------------------------------------------

      if (mobile.and.(withgap.eq..false.)) then

c       MSP BEGIN
        do k=1,vardim(3)

c        Cluster all stationary frontal sections
         do i=1,vardim(1)
           do j=1,vardim(2)

              if ( ( abs(condi1(i,j,k))       .lt.minadv ).and.
     >             ( abs(condi1(i,j,k)-misdat).gt.eps    ).and.
     >             ( abs(condi2(i,j,k)-misdat).gt.eps    ) )
     >        then
                 clusterin(i,j) = 1
              else
                 clusterin(i,j) = 0
              endif
           enddo
         enddo
         call clustering(belowc,nbelow,clusterin,
     >                   vardim(1),vardim(2),grid)

c        Cluster all mobile frontal sections
         do i=1,vardim(1)
           do j=1,vardim(2)

              if ( ( abs(condi1(i,j,k))       .ge.minadv ).and.
     >             ( abs(condi1(i,j,k)-misdat).gt.eps    ).and.
     >             ( abs(condi2(i,j,k)-misdat).gt.eps    ) )
     >        then
                 clusterin(i,j) = 1
              else
                 clusterin(i,j) = 0
              endif
           enddo
         enddo
         call clustering(abovec,nabove,clusterin,
     >                   vardim(1),vardim(2),grid)

c       Init the removal flag for each stationary section
        do j=1,nbelow
          belowl(i) = 0
        enddo

c       Check wether the stationary clusters are connected to at least two
c       different mobilr clusters -> belowl=-1
        do i=1,vardim(1)
          do j=1,vardim(2)

             ind = belowc(i,j)
             if ( ind.ne.0 ) then
               il = i-1
               ir = i+1
               ju = j+1
               jd = j-1
               if (il.lt.1         ) il = 1
               if (ir.gt.vardim(1) ) ir = vardim(1)
               if (jd.lt.1         ) jd = 1
               if (ju.gt.vardim(2) ) ju = vardim(2)

               ival = abovec(il,j )
               if ( (ival.ne.0).and.(belowl(ind).eq.0) ) then
                  belowl(ind) = ival
               elseif ( (ival.ne.0).and.(belowl(ind).ne.ival) ) then
                  belowl(ind) = -1
               endif

               ival = abovec(ir,j )
               if ( (ival.ne.0).and.(belowl(ind).eq.0) ) then
                  belowl(ind) = ival
               elseif ( (ival.ne.0).and.(belowl(ind).ne.ival) ) then
                  belowl(ind) = -1
               endif

               ival = abovec(i ,jd)
               if ( (ival.ne.0).and.(belowl(ind).eq.0) ) then
                  belowl(ind) = ival
               elseif ( (ival.ne.0).and.(belowl(ind).ne.ival) ) then
                  belowl(ind) = -1
               endif

               ival = abovec(i ,ju)
               if ( (ival.ne.0).and.(belowl(ind).eq.0) ) then
                  belowl(ind) = ival
               elseif ( (ival.ne.0).and.(belowl(ind).ne.ival) ) then
                  belowl(ind) = -1
               endif

               ival = abovec(il,jd)
               if ( (ival.ne.0).and.(belowl(ind).eq.0) ) then
                  belowl(ind) = ival
               elseif ( (ival.ne.0).and.(belowl(ind).ne.ival) ) then
                  belowl(ind) = -1
               endif

               ival = abovec(ir,jd)
               if ( (ival.ne.0).and.(belowl(ind).eq.0) ) then
                  belowl(ind) = ival
               elseif ( (ival.ne.0).and.(belowl(ind).ne.ival) ) then
                  belowl(ind) = -1
               endif

               ival = abovec(il,ju)
               if ( (ival.ne.0).and.(belowl(ind).eq.0) ) then
                  belowl(ind) = ival
               elseif ( (ival.ne.0).and.(belowl(ind).ne.ival) ) then
                  belowl(ind) = -1
               endif

               ival = abovec(ir,ju)
               if ( (ival.ne.0).and.(belowl(ind).eq.0) ) then
                  belowl(ind) = ival
               elseif ( (ival.ne.0).and.(belowl(ind).ne.ival) ) then
                  belowl(ind) = -1
               endif

             endif

            enddo
          enddo

c         Remove all stationary sections which have not two mobile neighbors
          do i=1,vardim(1)
            do j=1,vardim(2)

               ind = belowc(i,j)
               if ( ind.gt.0 ) then
                 if ( belowl(ind).ne.-1 ) then
                    help(i,j,k)   = misdat
                    condi1(i,j,k) = misdat
                    condi2(i,j,k) = misdat
                 endif
               endif

            enddo
          enddo

        enddo

      elseif ( mobile.and.(withgap.eq..true.) ) then

        print*,'Applying advection criterion'
c 	    check if advection velocity is below threshold
        where(abs(condi1).lt.minadv)          condi1 = misdat
        where(abs(condi1-misdat).lt.eps)      condi2 = misdat
        where(abs(condi2-misdat).lt.eps)      help   = misdat

      endif

c     --------------------------------------------------------------
c     Remove very short objects
c     --------------------------------------------------------------

      if (objects) then

       do k=1,vardim(3)

          where ((abs(help(:,:,k)-misdat).lt.eps))
            clusterin(:,:)=0.
          elsewhere
            clusterin(:,:)=1.
          endwhere

     	 call clustering(clusterout,ntot,clusterin,
     >			        vardim(1),vardim(2),grid)

          do j=1,ntot

            if (count(clusterout(:,:)==j).lt.nmin) then
               where (clusterout(:,:)==j) help(:,:,k) = misdat
            endif

          enddo

       enddo

c     --------------------------------------------------------------
c     correction to fvel_var/fdir_var after min objects removal
c     --------------------------------------------------------------

       if (winds.or.mobile) then

        where(abs(help-misdat).lt.eps) condi1 = misdat
        where(abs(help-misdat).lt.eps) condi2 = misdat

       end if

c     --------------------------------------------------------------
c      correction to frontal areas after min objects removal
c     --------------------------------------------------------------

       if (area) then

        clusterout = 0

        do k=1,vardim(3)

          where ((abs(farea(:,:,k)-misdat).lt.eps))
            clusterin=0.
          elsewhere
            clusterin=1.
          endwhere

     	  call clustering(clusterout,ntot,clusterin,
     >                      vardim(1),vardim(2),grid)

        do  j = 1, ntot
            onecluster=0
            where (clusterout(:,:)==j) onecluster=clusterout

            if ( .not. any( (onecluster(:,:)*help(:,:,k)).gt.0) ) then
                where (clusterout(:,:)==j) farea(:,:,k)=misdat
            endif

          enddo

        enddo
c     --------------------------------------------------------------
c     end if area
c     --------------------------------------------------------------

       endif

c     --------------------------------------------------------------
c     end if objects
c     --------------------------------------------------------------

      endif

c     --------------------------------------------------------------
c      Output
c     --------------------------------------------------------------

        if (winds) then

          call putdat(idout,'fvel_'//trim(invar),time(n),0,condi1,
     >                ierr)
          call putdat(idout,'fdir_'//trim(invar),time(n),0,condi2,
     >                ierr)
        endif

        if (area) then
          call putdat(idout,'farea_'//trim(invar),time(n),0,farea,
     >                ierr)

        endif

        call putdat(idout,'fronts_'//trim(invar),time(n),0,
     >               help,ierr)

        call putdat(idout,invar,time(n),0,tfield,ierr)

        if (mobile) then
           help = 0.
           where( abs(condi1).lt.minadv     ) help = 1.
           where( abs(condi1).ge.minadv     ) help = 2.
           where( abs(condi1-misdat).le.eps ) help = 0.

           call putdat(idout,'fmobile_'//trim(invar),time(n),0,
     >               help,ierr)

        endif

c     --------------------------------------------------------------
      enddo ! time loop
c     --------------------------------------------------------------

c     --------------------------------------------------------------
c     Close netcdf files
c     --------------------------------------------------------------

      call clscdf(idin,ierr)
      call clscdf(idout,ierr)

c     --------------------------------------------------------------
c     Deallocate required arrays
c     --------------------------------------------------------------

       if (allocated(help))  deallocate(help)
       if (allocated(gradx)) deallocate(gradx)
       if (allocated(grady)) deallocate(grady)
       if (allocated(help1)) deallocate(help1)
       if (allocated(help2)) deallocate(help2)
       if (allocated(condi1)) deallocate(condi1)
       if (allocated(condi2)) deallocate(condi2)
       if (allocated(vert))   deallocate(vert)
       if (allocated(tfield)) deallocate(tfield)
       if (allocated(farea))  deallocate(farea)
       if (allocated(clusterout)) deallocate(clusterout)
       if (allocated(clusterin))  deallocate(clusterin)
       if (allocated(onecluster)) deallocate(onecluster)

        print*," "
        print*," DONE"
        print*," "

c     -----------------------------------------------
      end program fronts
c     -----------------------------------------------

max_stars_repo_name: ruestefa/stormtrack

max_issues_repo_name: ruestefa/stormtrack

max_forks_repo_name: ruestefa/stormtrack


text:
\ https://rosettacode.org/wiki/Sieve_of_Eratosthenes#Forth

: PRIME? ( n -- ? ) HERE + C@ 0= ;
: COMPOSITE! ( n -- ) HERE + 1 SWAP C! ;

: SIEVE ( n -- )
  HERE OVER ERASE
  2
  BEGIN
    2DUP DUP * >
  WHILE
    DUP PRIME? IF
      2DUP DUP * DO
        I COMPOSITE!
      DUP +LOOP
    THEN
    1+
  REPEAT
  DROP
  ." Primes: " 2 DO I PRIME? IF I . THEN LOOP CR ;

DECIMAL 100 SIEVE

max_stars_repo_name: happz/ducky

max_issues_repo_name: happz/ducky-legacy

max_forks_repo_name: happz/ducky-legacy


text:
      SUBROUTINE SET3SCAN(PARMNAME,LUPROMPT,LUCNTL,
     *      VALUE1,DELVALUE,NVALUES,FACTOR,IEND)
      IMPLICIT REAL*8 (A-H,O-Z)
C
C SET3SCAN IS A HANDY WAY TO HAVE PROGRAMS PROMPT THE USER FOR PARAMETER
C CONTROLS WHERE A PARAMETER IS TO VARY IN EQUAL INCREMENTS OVER A RANGE
C
C EXAMPLE:
C
C    A PROGRAM NEEDS TO DO COMPUTATIONS INVOLVING A SET OF EQUALLY
C    SPACED LONGITUDES. INSTEAD OF WRITING CODE TO PROMPT FOR AND READ
C    THE LONGITUDE RANGE, THE CALLER DOES --
C
C   CALL SET3SCAN('LONGITUDE RANGE',6,5,ALONG1,DLONG,NLONG,DEGRAD,IEND)
C
C    THIS ROUTINE WILL PROMPT FOR THE FIRST VALUE, INCREMENT, AND NUMBER
C    OF VALUES WANTED. THESE WILL BE OUTPUT TO THE CALLER THROUGH THE
C    CALLING SEQUENCE VIA VARIABLES ALONG1, DLONG, NLONG AFTER
C    CONVERTING FROM DEGREES(THE USER'S INPUT UNITS) TO RADIANS(THE
C    CALLER'S WORKING UNITS).
C
C
C VARIABLE DIM TYPE I/O DESCRIPTION
C -------- --- ---- --- -----------
C
C PARMNAME  1 CH*(*) I  TEXT TO BE PART OF THE PROMPT MESSAGE. USEFUL
C                       WHEN MORE THAN ONE SET OF CONTROLS ARE USED BY
C                       THE CALLER. THE TEXT WILL BE INSERTED IN A
C                       LINE OF PROMPT TEXT AS FOLLOWS:
C
C                              <<< parmname CONTROLS >>>
C
C LUPROMPT  1  I*4   I  THE FORTRAN UNIT NUMBER TO WHICH THE PROMPT
C                       MESSAGE IS WRITTEN. IF ZERO OR NEGATIVE, NO
C                       PROMPT IS WRITTEN.
C
C LUCNTL    1  I*4   I  THE FORTRAN UNIT NUMBER FROM WHICH THE DATA IS
C                       TO BE READ.
C
C VALUE1    1  R*8  I/O THE FIRST VALUE IN THE RANGE OF PARAMETER VALUES
C
C                       INPUT: THE VALUE THAT WILL BE OUTPUT UNLESS THE
C                              USER OVERRIDES IT. INPUT UNITS ARE THE
C                              SAME AS OUTPUT UNITS.
C
C                       OUTPUT: THE VALUE RETURNED TO THE CALLER.
C
C DELVALUE  1  R*8  I/O THE INCREMENT BETWEEN PARAMETER VALUES.
C
C                       INPUT/OUTPUT IS AS DESCRIBED UNDER VALUE1.
C
C NVALUES   1  I*4  I/O THE NUMBER OF PARAMETER VALUES.
C
C                       INPUT/OUTPUT IS AS DESCRIBED UNDER VALUE1.
C
C FACTOR    1  R*8   I  A FACTOR BY WHICH THE USER'S INPUT VALUE1 AND
C                       DELVALUE ARE TO BE DIVIDED BEFORE THEY ARE
C                       RETURNED TO THE CALLER. IF FACTOR = 0.D0, THEN
C                       1.D0 IS USED.
C
C IEND      1  I*4   O  AN END-OF-FILE RETURN FLAG.
C
C                       = 0, NO END OF FILE WAS SENSED. OUTPUT ITEMS
C                            ARE AS THE USER HAS SPECIFIED.
C                       = 1, END OF FILE WAS SENSED FOR UNIT LUCNTL. IN
C                            THIS CASE, VALUE1, DELVALUE, AND NVALUES
C                            ARE UNCHANGED.
C
C***********************************************************************
C
C BY C PETRUZZO, 6/84.
C    MODIFIED....
C
C***********************************************************************
C
      CHARACTER*(*) PARMNAME
C
C INITIALIZATION
C
      IEND = 1   ! WILL BE RESET TO ZERO IF NO END-FILE IS ENCOUNTERED
      FAC = FACTOR
      IF(FAC.EQ.0.D0) FAC=1.D0
      V1 = VALUE1*FAC
      DELV = DELVALUE*FAC
      NVAL = NVALUES
C
C PROMPT
C
      IF(LUPROMPT.GT.0) WRITE(LUPROMPT,100) PARMNAME,V1,DELV,NVAL
  100 FORMAT(/,
     * ' <<< ',A,' CONTROLS >>>'//,
     * '   FIRST VALUE = ',G19.12/,
     * '   INCREMENT =   ',G19.12/,
     * '   NUM VALUES =  ',I10//,
     * '   ENTER NEW CONTROLS -----> ',$)
C
C READ OVERRIDES
C
      READ(LUCNTL,*,END=999) V1,DELV,NVAL
      IEND = 0
C
C REPORT ON VALUES ENTERED
C
      IF(LUPROMPT.GT.0) WRITE(LUPROMPT,101) V1,DELV,NVAL
  101 FORMAT(/,'   ENTERED --> ',2G20.12,1X,I6/)
C
C SET OUTPUT VARIABLES
C
      NVALUES = NVAL
      VALUE1 = V1/FAC
      DELVALUE = DELV/FAC
C
  999 CONTINUE
      RETURN
      END

ext: for


text:
      SUBROUTINE QUIKVIS1(IDTARG,TARGNAMES,KTARGTYP,TARGPARM,IERR)
      IMPLICIT REAL*8 (A-H,O-Z)
C
C THIS ROUTINE IS PART OF THE QUIKVIS PROGRAM.  IT LOADS THE
C IDTARG, TARGNAMES, AND TARGPARM ARRAYS VIA AN INTERNAL ALGORITHM(WHEN
C THE SURVEY OPTION IS ON) OR VIA THE TOSS TARGET CATALOGUE READER.
C
C
C VARIABLE      DIM       TYPE I/O DESCRIPTION
C --------      ---       ---- --- -----------
C
C IDTARG      MAXTARGS     I*4  I  DESCRIBED IN QUIKVIS(=MAIN) PROLOGUE.
C
C TARGNAMES   MAXTARGS    CH*16 I  DESCRIBED IN QUIKVIS(=MAIN) PROLOGUE.
C
C KTARGTYP    MAXTARGS     I*4  I  DESCRIBED IN QUIKVIS(=MAIN) PROLOGUE.
C
C TARGPARM NPARMS,MAXTARGS R*8  I  DESCRIBED IN QUIKVIS(=MAIN) PROLOGUE.
C
C IERR           1         I*4  O  ERROR RETURN FLAG
C                                  =0, NO ERROR
C                                  =OTHERWISE, ERROR.
C
C***********************************************************************
C
C BY C PETRUZZO/GFSC/742.   2/86.
C       MODIFIED....
C
C***********************************************************************
C
      INCLUDE 'QUIKVIS.INC'
C
      CHARACTER*16 TARGNAMES(MAXTARGS)
      INTEGER*4 IDTARG(MAXTARGS)
      INTEGER*4 KTARGTYP(MAXTARGS)
C
C    TARGPARM ARRAY DOES DOUBLE DUTY.  (1) ON THE CALL TO RDCAT, IT IS
C    THE ARRAY INTO WHICH ALL TARGET PARAMETERS ARE LOADED; HENCE, THE
C    FIRST DIMENSION MUST BE LARGE ENOUGH TO HOLD ALL PARAMETER VALUES
C    PRESENT IN THE CATALOGUE; (2) AFTER RDCAT RETURNS, THE 3RD, 4TH,
C    AND 5TH PARAMETER POSITIONS ARE LOADED WITH UNIT VECTORS(ZEROS IF
C    TARGET IS MOVING).
      REAL*8 TARGPARM(NPARMS,MAXTARGS)
C
      CHARACTER*18 DATETIME
C
C
      IERR = 0
C
C
      CALL QUIKVIS999(-1,R8DATA,I4DATA,L4DATA)
C
C
C ******************************************************
C * FOR SURVEY OPTION, DEFINE STAR NAMES AND LOCATIONS *
C ******************************************************
C
C
C IF THE SURVEY OPTION IS ON, THE TARGET CATALOGUE IS NOT READ.
C TARGNAMES, TARGPARM, AND KTARGTYP ARE SET HERE.
C
      IF(DOSURVEY) THEN
        DELRA = TWOPI/(NRASURVEY-1)
        DELDEC = PI/(NDECSURVEY-1)
        ITARG = 0
        DO IDEC=1,NDECSURVEY
          DEC = HALFPI - (IDEC-1)*DELDEC
          DO IRA=1,NRASURVEY
            RA = (IRA-1)*DELRA
            ITARG = ITARG + 1
            TARGNAMES(ITARG) = 'SURVEY OPTION'
            TARGPARM(1,ITARG) = RA
            TARGPARM(2,ITARG) = DEC
            CALL XYZSPH(1,RA,DEC,1.D0,TARGPARM(3,ITARG),1.D0)  !UNIT VEC
            KTARGTYP(ITARG) = 3  ! FIXED TARGET
            END DO
          END DO
        GO TO 9999
        END IF
C
C
C
C *******************************************************************
C *  IF SURVEY OPTION IS OFF, READ TARGET CATLG FOR SPECIFIC TARGS  *
C *******************************************************************
C
C
C  READ TARGET INFORMATION FROM THE CATALOGUE:
C
      WRITE(LUPROMPT,5) DATETIME(0)
    5 FORMAT(/,' READING CATALOGUE INFORMATION.   ',A/)
      CALL RDCAT(LUCATLG,MAXTARGS,IDTARG,TARGNAMES,KTARGTYP,NPARMS,
     *              TARGPARM,NFOUND,0,LUERR,IER1)
      WRITE(LUPROMPT,6) DATETIME(0)
    6 FORMAT(/,' END READING CATALOGUE INFORMATION.   ',A/)
C
      IF(IER1 .NE. 0) THEN
        IERR = 1
        WRITE(LUERR,*) ' QUIKVIS1. TARGET CATALOGUE ERROR.'
        GO TO 9999
        END IF
C
C ERROR CHECK. ALL TARGETS FOUND ?
C
      IF(NUMTARGS.NE.NFOUND) THEN
        IERR = 1
        WRITE(LUERR,*)
     *     ' QUIKVIS1. SOME TARGETS ARE ABSENT FROM THE CATALOGUE.'
        GO TO 9999
        END IF
C
C ERROR CHECK. ARE ALL ENTRIES FIXED CELESTIAL OR MOVING CELESTIAL ?
C
      DO I=1,NUMTARGS
        IF(KTARGTYP(I).NE.3 .AND. KTARGTYP(I).NE.1) THEN
          IERR = 1
          WRITE(LUERR,8001) IDTARG(I)
 8001     FORMAT(/,
     *      ' QUIKVIS1 ERROR. NON-CELESTIAL TARGET SPECIFIED.   ID=',I5)
          END IF
        END DO
      IF(IERR.NE.0) GO TO 9999
C
C ALL IS OK. PUT UNIT VECTOR IN TARGPARM IF THE TARGET IS FIXED, ZEROS
C IF MOVING. IF MOVING, SET THE NEEDTIME VARIABLE TO .TRUE.
C
        DO ITARG=1,NUMTARGS
          IF(KTARGTYP(ITARG).EQ.3) THEN
            CALL XYZSPH(1,TARGPARM(1,ITARG),TARGPARM(2,ITARG),1.D0,
     *           TARGPARM(3,ITARG),1.D0)
          ELSE
            NEEDTIME = .TRUE.
            DO I=1,5
              TARGPARM(I,ITARG) = 0.D0
              END DO
            END IF
          END DO
C
C
C WRAP UP. STORE NEEDTIME IN QUIKVIS999 ARRAYS FOR RETRIEVAL ELSEWHERE.
C
      CALL QUIKVIS999(1,R8DATA,I4DATA,L4DATA)
C
C
 9999 CONTINUE
      RETURN
      END


text:
      SUBROUTINE ZTRMV(UPLO,TRANS,DIAG,N,A,LDA,X,INCX)
*     .. Scalar Arguments ..
      INTEGER INCX,LDA,N
      CHARACTER DIAG,TRANS,UPLO
*     ..
*     .. Array Arguments ..
      DOUBLE COMPLEX A(LDA,*),X(*)
*     ..
*
*  Purpose
*  =======
*
*  ZTRMV  performs one of the matrix-vector operations
*
*     x := A*x,   or   x := A**T*x,   or   x := A**H*x,
*
*  where x is an n element vector and  A is an n by n unit, or non-unit,
*  upper or lower triangular matrix.
*
*  Arguments
*  ==========
*
*  UPLO   - CHARACTER*1.
*           On entry, UPLO specifies whether the matrix is an upper or
*           lower triangular matrix as follows:
*
*              UPLO = 'U' or 'u'   A is an upper triangular matrix.
*
*              UPLO = 'L' or 'l'   A is a lower triangular matrix.
*
*           Unchanged on exit.
*
*  TRANS  - CHARACTER*1.
*           On entry, TRANS specifies the operation to be performed as
*           follows:
*
*              TRANS = 'N' or 'n'   x := A*x.
*
*              TRANS = 'T' or 't'   x := A**T*x.
*
*              TRANS = 'C' or 'c'   x := A**H*x.
*
*           Unchanged on exit.
*
*  DIAG   - CHARACTER*1.
*           On entry, DIAG specifies whether or not A is unit
*           triangular as follows:
*
*              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
*
*              DIAG = 'N' or 'n'   A is not assumed to be unit
*                                  triangular.
*
*           Unchanged on exit.
*
*  N      - INTEGER.
*           On entry, N specifies the order of the matrix A.
*           N must be at least zero.
*           Unchanged on exit.
*
*  A      - COMPLEX*16       array of DIMENSION ( LDA, n ).
*           Before entry with  UPLO = 'U' or 'u', the leading n by n
*           upper triangular part of the array A must contain the upper
*           triangular matrix and the strictly lower triangular part of
*           A is not referenced.
*           Before entry with UPLO = 'L' or 'l', the leading n by n
*           lower triangular part of the array A must contain the lower
*           triangular matrix and the strictly upper triangular part of
*           A is not referenced.
*           Note that when  DIAG = 'U' or 'u', the diagonal elements of
*           A are not referenced either, but are assumed to be unity.
*           Unchanged on exit.
*
*  LDA    - INTEGER.
*           On entry, LDA specifies the first dimension of A as declared
*           in the calling (sub) program. LDA must be at least
*           max( 1, n ).
*           Unchanged on exit.
*
*  X      - COMPLEX*16       array of dimension at least
*           ( 1 + ( n - 1 )*abs( INCX ) ).
*           Before entry, the incremented array X must contain the n
*           element vector x. On exit, X is overwritten with the
*           tranformed vector x.
*
*  INCX   - INTEGER.
*           On entry, INCX specifies the increment for the elements of
*           X. INCX must not be zero.
*           Unchanged on exit.
*
*  Further Details
*  ===============
*
*  Level 2 Blas routine.
*  The vector and matrix arguments are not referenced when N = 0, or M = 0
*
*  -- Written on 22-October-1986.
*     Jack Dongarra, Argonne National Lab.
*     Jeremy Du Croz, Nag Central Office.
*     Sven Hammarling, Nag Central Office.
*     Richard Hanson, Sandia National Labs.
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE COMPLEX ZERO
      PARAMETER (ZERO= (0.0D+0,0.0D+0))
*     ..
*     .. Local Scalars ..
      DOUBLE COMPLEX TEMP
      INTEGER I,INFO,IX,J,JX,KX
      LOGICAL NOCONJ,NOUNIT
*     ..
*     .. External Functions ..
      LOGICAL LSAME
      EXTERNAL LSAME
*     ..
*     .. External Subroutines ..
      EXTERNAL XERBLA
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC DCONJG,MAX
*     ..
*
*     Test the input parameters.
*
      INFO = 0
      IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
          INFO = 1
      ELSE IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND.
     +         .NOT.LSAME(TRANS,'C')) THEN
          INFO = 2
      ELSE IF (.NOT.LSAME(DIAG,'U') .AND. .NOT.LSAME(DIAG,'N')) THEN
          INFO = 3
      ELSE IF (N.LT.0) THEN
          INFO = 4
      ELSE IF (LDA.LT.MAX(1,N)) THEN
          INFO = 6
      ELSE IF (INCX.EQ.0) THEN
          INFO = 8
      END IF
      IF (INFO.NE.0) THEN
          CALL XERBLA('ZTRMV ',INFO)
          RETURN
      END IF
*
*     Quick return if possible.
*
      IF (N.EQ.0) RETURN
*
      NOCONJ = LSAME(TRANS,'T')
      NOUNIT = LSAME(DIAG,'N')
*
*     Set up the start point in X if the increment is not unity. This
*     will be  ( N - 1 )*INCX  too small for descending loops.
*
      IF (INCX.LE.0) THEN
          KX = 1 - (N-1)*INCX
      ELSE IF (INCX.NE.1) THEN
          KX = 1
      END IF
*
*     Start the operations. In this version the elements of A are
*     accessed sequentially with one pass through A.
*
      IF (LSAME(TRANS,'N')) THEN
*
*        Form  x := A*x.
*
          IF (LSAME(UPLO,'U')) THEN
              IF (INCX.EQ.1) THEN
                  DO 20 J = 1,N
                      IF (X(J).NE.ZERO) THEN
                          TEMP = X(J)
                          DO 10 I = 1,J - 1
                              X(I) = X(I) + TEMP*A(I,J)
   10                     CONTINUE
                          IF (NOUNIT) X(J) = X(J)*A(J,J)
                      END IF
   20             CONTINUE
              ELSE
                  JX = KX
                  DO 40 J = 1,N
                      IF (X(JX).NE.ZERO) THEN
                          TEMP = X(JX)
                          IX = KX
                          DO 30 I = 1,J - 1
                              X(IX) = X(IX) + TEMP*A(I,J)
                              IX = IX + INCX
   30                     CONTINUE
                          IF (NOUNIT) X(JX) = X(JX)*A(J,J)
                      END IF
                      JX = JX + INCX
   40             CONTINUE
              END IF
          ELSE
              IF (INCX.EQ.1) THEN
                  DO 60 J = N,1,-1
                      IF (X(J).NE.ZERO) THEN
                          TEMP = X(J)
                          DO 50 I = N,J + 1,-1
                              X(I) = X(I) + TEMP*A(I,J)
   50                     CONTINUE
                          IF (NOUNIT) X(J) = X(J)*A(J,J)
                      END IF
   60             CONTINUE
              ELSE
                  KX = KX + (N-1)*INCX
                  JX = KX
                  DO 80 J = N,1,-1
                      IF (X(JX).NE.ZERO) THEN
                          TEMP = X(JX)
                          IX = KX
                          DO 70 I = N,J + 1,-1
                              X(IX) = X(IX) + TEMP*A(I,J)
                              IX = IX - INCX
   70                     CONTINUE
                          IF (NOUNIT) X(JX) = X(JX)*A(J,J)
                      END IF
                      JX = JX - INCX
   80             CONTINUE
              END IF
          END IF
      ELSE
*
*        Form  x := A**T*x  or  x := A**H*x.
*
          IF (LSAME(UPLO,'U')) THEN
              IF (INCX.EQ.1) THEN
                  DO 110 J = N,1,-1
                      TEMP = X(J)
                      IF (NOCONJ) THEN
                          IF (NOUNIT) TEMP = TEMP*A(J,J)
                          DO 90 I = J - 1,1,-1
                              TEMP = TEMP + A(I,J)*X(I)
   90                     CONTINUE
                      ELSE
                          IF (NOUNIT) TEMP = TEMP*DCONJG(A(J,J))
                          DO 100 I = J - 1,1,-1
                              TEMP = TEMP + DCONJG(A(I,J))*X(I)
  100                     CONTINUE
                      END IF
                      X(J) = TEMP
  110             CONTINUE
              ELSE
                  JX = KX + (N-1)*INCX
                  DO 140 J = N,1,-1
                      TEMP = X(JX)
                      IX = JX
                      IF (NOCONJ) THEN
                          IF (NOUNIT) TEMP = TEMP*A(J,J)
                          DO 120 I = J - 1,1,-1
                              IX = IX - INCX
                              TEMP = TEMP + A(I,J)*X(IX)
  120                     CONTINUE
                      ELSE
                          IF (NOUNIT) TEMP = TEMP*DCONJG(A(J,J))
                          DO 130 I = J - 1,1,-1
                              IX = IX - INCX
                              TEMP = TEMP + DCONJG(A(I,J))*X(IX)
  130                     CONTINUE
                      END IF
                      X(JX) = TEMP
                      JX = JX - INCX
  140             CONTINUE
              END IF
          ELSE
              IF (INCX.EQ.1) THEN
                  DO 170 J = 1,N
                      TEMP = X(J)
                      IF (NOCONJ) THEN
                          IF (NOUNIT) TEMP = TEMP*A(J,J)
                          DO 150 I = J + 1,N
                              TEMP = TEMP + A(I,J)*X(I)
  150                     CONTINUE
                      ELSE
                          IF (NOUNIT) TEMP = TEMP*DCONJG(A(J,J))
                          DO 160 I = J + 1,N
                              TEMP = TEMP + DCONJG(A(I,J))*X(I)
  160                     CONTINUE
                      END IF
                      X(J) = TEMP
  170             CONTINUE
              ELSE
                  JX = KX
                  DO 200 J = 1,N
                      TEMP = X(JX)
                      IX = JX
                      IF (NOCONJ) THEN
                          IF (NOUNIT) TEMP = TEMP*A(J,J)
                          DO 180 I = J + 1,N
                              IX = IX + INCX
                              TEMP = TEMP + A(I,J)*X(IX)
  180                     CONTINUE
                      ELSE
                          IF (NOUNIT) TEMP = TEMP*DCONJG(A(J,J))
                          DO 190 I = J + 1,N
                              IX = IX + INCX
                              TEMP = TEMP + DCONJG(A(I,J))*X(IX)
  190                     CONTINUE
                      END IF
                      X(JX) = TEMP
                      JX = JX + INCX
  200             CONTINUE
              END IF
          END IF
      END IF
*
      RETURN
*
*     End of ZTRMV .
*
      END

max_stars_repo_name: NetlibRepository/lapack-www

max_issues_repo_name: NetlibRepository/lapack-www

max_forks_repo_name: NetlibRepository/lapack-www


text:
20071010 11:29:26 nbsp Apparently here is the closest, unless you hold one ;)

http://slashdot.org/anniversary.pl?view_id269 Users/WesHardaker

20091019 18:42:42 nbsp pretty sure RBU already has a page Users/StevenDaubert

20091019 18:44:17 nbsp yeah its rockband no space Users/StevenDaubert


text:
      subroutine writenumber (xnum)
c******************************************************************************
c     This subroutine decides on the order-of-magnitude of a number and
c     writes it to a character string in a reasonable format
c******************************************************************************

      include 'Pstuff.com'
      real*4 xnum, lognum
      integer numdec

      lognum = alog10(abs(xnum))
      if     (lognum .ge. 6.) then
         write (errmess,1002)
      elseif (lognum .ge. 0.) then
         numdec = 5 - nint(lognum)
         write (errmess,1001) numdec
      else
         write (errmess,1003)
      endif
      write (array,errmess) xnum
      return

c*****format statements
1001  format ('(f8.',i1,'$)')
1002  format ('(1pe8.1$)')
1003  format ('(f8.5$)')

      end