text:
c
c -------------------------------------------------------------
c
      subroutine outgmsh(qin,nvar,mptr,irr,
     1                  mitot,mjtot,qx,qy,lwidth,lstgrd,
     2                  dx,dy,xlow,ylow,time)
c
      implicit double precision (a-h,o-z)
      include "cirr.i"
      dimension irr(mitot,mjtot), q(mitot,mjtot,nvar)
      dimension qin(mitot,mjtot,nvar)
      dimension qx(mitot,mjtot,nvar),qy(mitot,mjtot,nvar)
      dimension qxx(mitot,mjtot,nvar),qyy(mitot,mjtot,nvar)
      dimension qxy(mitot,mjtot,nvar)
      dimension qxxx(mitot,mjtot,nvar),qxyy(mitot,mjtot,nvar)
      dimension qxxy(mitot,mjtot,nvar),qyyy(mitot,mjtot,nvar)
      dimension valprim(4)
      dimension exactsoln(1)
      dimension exactval(mitot,mjtot,nvar)
      dimension dlimit(mitot,mjtot)
      dimension recon(4,nvar), corners(4,2)
      dimension numvert(irrsize)
      dimension dnumsol(nvar)
      common /userdt/ cfl,gamma,gamma1,xprob,yprob,alpha,Re,iprob,
     .                ismp,gradThreshold,ilts,ibc,imodal
      common /order2/ ssw, quad, nolimiter
      logical  eout, quad, nolimiter
      logical pwconst

      data      spval/-1.e10/, eout/.true./
c      character*23  filename
      common /dump/ uvel,vvel
      dimension uvel(849,848),vvel(848,849)

      q = qin

!      open(231, file = 'out.pos')
      ibuff = 4

      xlowb = xlow - lwidth*dx
      ylowb = ylow - lwidth*dy

      qx = 0.d0
      qy = 0.d0
      qxx = 0.d0
      qxy = 0.d0
      qyy = 0.d0
      qxxx = 0.d0
      qxxy = 0.d0
      qxyy = 0.d0
      qyyy = 0.d0
!
!      call vctoprm(q,q,mitot,mjtot)
c
c     ### call for exterior bcs at each stage so can use slopes
c    ## NOTE THAT BNDRY CELLS FROM OTHER GRIDS NOT SET
            xhigh = xlowb + mitot*dx
            yhigh = ylowb + mjtot*dy
            call pphysbdlin(xlowb,xhigh,ylowb,yhigh,level,mitot,mjtot,
     &                   nvar,q,time,dx,dy,qx,qy,irr,lstgrd)

      if(quad .and. (nolimiter .eqv. .true.) ) then ! yes in primitive variables, and without limiter (this is an accurate projection between the two
      call slopes(q,qx,qy,qxx,qxy,qyy,
     &               qxxx, qxxy, qxyy, qyyy,
     &               mitot,mjtot,irr,lstgrd,lwidth,dx,dy,
     &               xlowb,ylowb,nvar)

      call vctoprm(mitot,mjtot,nvar,
     . dx, dy, lstgrd, xlowb, ylowb, irr,
     . q,q,qx,qy,qxx, qxy, qyy, qxxx, qxxy, qxyy, qyyy)

      call slopes(q,qx,qy,qxx,qxy,qyy,
     &               qxxx, qxxy, qxyy, qyyy,
     &               mitot,mjtot,irr,lstgrd,lwidth,dx,dy,
     &               xlowb,ylowb,nvar)

      elseif(quad .and. (nolimiter .eqv. .false.) ) then ! yes in primitive variables, and with limiter
                                                      ! only works for p = 1
      call slopes(q,qx,qy,qxx,qxy,qyy,
     &               qxxx, qxxy, qxyy, qyyy,
     &               mitot,mjtot,irr,lstgrd,lwidth,dx,dy,
     &               xlowb,ylowb,nvar)

      call vctoprm(mitot,mjtot,nvar,
     . dx, dy, lstgrd, xlowb, ylowb, irr,
     . q,q,qx,qy,qxx, qxy, qyy, qxxx, qxxy, qxyy, qyyy)

      call slopes(q,qx,qy,qxx,qxy,qyy,
     &               qxxx, qxxy, qxyy, qyyy,
     &               mitot,mjtot,irr,lstgrd,lwidth,dx,dy,
     &               xlowb,ylowb,nvar)

      ! limiter is only applied to the primitive variables.
!      if(imodal .eq. 0) then
      call limiter(q,  qx, qy,  mitot, mjtot, irr, nvar,dx, dy,
     .               lwidth, lstgrd, xlowb, ylowb)
!      else
!      call limiter_modal(q,  qx, qy,  qxx,qxy,qyy,
!     .             qxxx, qxxy, qxyy, qyyy,
!     .               mitot, mjtot, irr, nvar,dx, dy,
!     .               lwidth, lstgrd, xlowb, ylowb)
!      endif

      elseif (ssw .ne. 0.d0) then   ! conserved variables
         call slopes(q,qx,qy,qxx,qxy,qyy,
     &               qxxx, qxxy, qxyy, qyyy,
     &               mitot,mjtot,irr,lstgrd,lwidth,dx,dy,
     &               xlowb,ylowb,nvar)
      endif

      ivar = 1
      write(231,*) 'View "U0 - T = ',time,'" {'
      do 15 i = lwidth+1 - ibuff, mitot-lwidth+ibuff
      do 15 j = lwidth+1-ibuff, mjtot-lwidth+ibuff
         kirr = irr(i,j)
         if (kirr .eq. -1) go to 15
         if (kirr .eq. lstgrd ) then
            xcen = xlowb + (dfloat(i)-.5d0)* dx
            ycen = ylowb + (dfloat(j)-.5d0)* dy

             xcorner = xlowb + (dfloat(i)-1.)* dx
             ycorner = ylowb + (dfloat(j)-1.)* dy
             do itimes = 1,4
                if (itimes .eq. 1) then ! get all 4 corners of mesh, pw constant sol
                   xc = xcorner
                   yc = ycorner
                else if (itimes .eq. 2) then
                   xc = xcorner
                   yc = ycorner + dy
                else if (itimes .eq. 3) then
                   xc = xcorner + dx
                   yc = ycorner + dy
                else if (itimes .eq. 4) then
                   xc = xcorner + dx
                   yc = ycorner
                endif

            corners(itimes,1) = xc
            corners(itimes,2) = yc

       deltax = xc-xcen
       deltay = yc-ycen

       do mmm = 1,nvar
       call evalU(dnumsol,deltax, deltay, dx, dy,
     .            q, qx, qy,
     .            qxx, qxy, qyy,
     .            qxxx, qxxy, qxyy, qyyy,
     .            i,j, kirr, mitot, mjtot, nvar, mmm)
       end do
       recon(:,itimes) = dnumsol(:)

!            recon(:,itimes) = q(i,j,:)
!     .   +(xc-xcen)*qx(i,j,:)
!     .   +(yc-ycen)*qy(i,j,:)
!     .   + 0.5d0*qxx(i,j,:)*( (xc-xcen)**2-(dx**2)/12.d0)
!     .   + (xc-xcen)*(yc-ycen)*qxy(i,j,:)
!     .   + 0.5d0*qyy(i,j,:)*( (yc-ycen)**2-(dy**2)/12.d0)
              end do

       write(231,*)  'SQ ( ', corners(1,1),",",
     .                        corners(1,2),",0,",
     .                        corners(2,1),",",
     .                        corners(2,2),",0,",
     .                        corners(3,1),",",
     .                        corners(3,2),",0,",
     .                        corners(4,1),",",
     .                        corners(4,2),",0)",
     .                        "{",
     .                        recon(ivar,1),",",
     .                        recon(ivar,2),",",
     .                        recon(ivar,3),",",
     .                        recon(ivar,4),
     .                        "};"

         else
            ivert = 1
            do 290 while (poly(ivert+1,1,kirr) .ne. -11.)
               ivert = ivert + 1
  290       continue
             if(ihob .eq. 0 ) then
                itri = ivert - 3
             else
                itri = ivert - 4 ! don't do the last one
             endif
                  idx1 = 1
                  do 291 it = 1, itri ! for each  triangle
                    idx2 = it + 1
                    idx3 = it + 2

                    x1 = poly(idx1,1,kirr)
                    y1 = poly(idx1,2,kirr)

                    x2 = poly(idx2,1,kirr)
                    y2 = poly(idx2,2,kirr)

                    x3 = poly(idx3,1,kirr)
                    y3 = poly(idx3,2,kirr)

                xcen = xcirr(kirr)
                ycen = ycirr(kirr)

            xc = x1
            yc = y1
            recon(:,1) = q(i,j,:)
     .   +(xc-xcen)*qx(i,j,:)
     .   +(yc-ycen)*qy(i,j,:)
     .   + 0.5d0*qxx(i,j,:)*( (xc-xcen)**2-(dx**2)/12.d0)
     .   + (xc-xcen)*(yc-ycen)*qxy(i,j,:)
     .   + 0.5d0*qyy(i,j,:)*( (yc-ycen)**2-(dy**2)/12.d0)

            xc = x2
            yc = y2
            recon(:,2) = q(i,j,:)
     .   +(xc-xcen)*qx(i,j,:)
     .   +(yc-ycen)*qy(i,j,:)
     .   + 0.5d0*qxx(i,j,:)*( (xc-xcen)**2-(dx**2)/12.d0)
     .   + (xc-xcen)*(yc-ycen)*qxy(i,j,:)
     .   + 0.5d0*qyy(i,j,:)*( (yc-ycen)**2-(dy**2)/12.d0)

            xc = x3
            yc = y3
            recon(:,3) = q(i,j,:)
     .   +(xc-xcen)*qx(i,j,:)
     .   +(yc-ycen)*qy(i,j,:)
     .   + 0.5d0*qxx(i,j,:)*( (xc-xcen)**2-(dx**2)/12.d0)
     .   + (xc-xcen)*(yc-ycen)*qxy(i,j,:)
     .   + 0.5d0*qyy(i,j,:)*( (yc-ycen)**2-(dy**2)/12.d0)

       write(231,*)  'ST ( ', x1,",",
     .                        y1,",0,",
     .                        x2,",",
     .                        y2,",0,",
     .                        x3,",",
     .                        y3,",0)",
     .                        "{",
     .                        recon(ivar,1),",",
     .                        recon(ivar,2),",",
     .                        recon(ivar,3),
     .                        "};"
291   continue

                  idx1 = ivert-2
                  x1 = poly(idx1,1,kirr)
                  y1 = poly(idx1,2,kirr)
                  do 11 it = 1, nbdry-1 ! for each  triangle
                    idx2 = it + 0
                    idx3 = it + 1

                    x2 = bdry(idx2,1,kirr)
                    y2 = bdry(idx2,2,kirr)

                    x3 = bdry(idx3,1,kirr)
                    y3 = bdry(idx3,2,kirr)

                xcen = xcirr(kirr)
                ycen = ycirr(kirr)

            xc = x1
            yc = y1
            recon(:,1) = q(i,j,:)
     .   +(xc-xcen)*qx(i,j,:)
     .   +(yc-ycen)*qy(i,j,:)
     .   + 0.5d0*qxx(i,j,:)*( (xc-xcen)**2-(dx**2)/12.d0)
     .   + (xc-xcen)*(yc-ycen)*qxy(i,j,:)
     .   + 0.5d0*qyy(i,j,:)*( (yc-ycen)**2-(dy**2)/12.d0)

            xc = x2
            yc = y2
            recon(:,2) = q(i,j,:)
     .   +(xc-xcen)*qx(i,j,:)
     .   +(yc-ycen)*qy(i,j,:)
     .   + 0.5d0*qxx(i,j,:)*( (xc-xcen)**2-(dx**2)/12.d0)
     .   + (xc-xcen)*(yc-ycen)*qxy(i,j,:)
     .   + 0.5d0*qyy(i,j,:)*( (yc-ycen)**2-(dy**2)/12.d0)

            xc = x3
            yc = y3
            recon(:,3) = q(i,j,:)
     .   +(xc-xcen)*qx(i,j,:)
     .   +(yc-ycen)*qy(i,j,:)
     .   + 0.5d0*qxx(i,j,:)*( (xc-xcen)**2-(dx**2)/12.d0)
     .   + (xc-xcen)*(yc-ycen)*qxy(i,j,:)
     .   + 0.5d0*qyy(i,j,:)*( (yc-ycen)**2-(dy**2)/12.d0)

       write(231,*)  'ST ( ', x1,",",
     .                        y1,",0,",
     .                        x2,",",
     .                        y2,",0,",
     .                        x3,",",
     .                        y3,",0)",
     .                        "{",
     .                        recon(ivar,1),",",
     .                        recon(ivar,2),",",
     .                        recon(ivar,3),
     .                        "};"
11    continue

         endif
 15   continue
      write(231,*) "};"

!      write(231,*) 'View "error" {'
!      do 915 i = lwidth+1 - ibuff, mitot-lwidth+ibuff
!      do 915 j = lwidth+1-ibuff, mjtot-lwidth+ibuff
!         kirr = irr(i,j)
!         if (kirr .eq. -1) go to 915
!         if (kirr .eq. lstgrd ) then
!            xcen = xlowb + (dfloat(i)-.5d0)* dx
!            ycen = ylowb + (dfloat(j)-.5d0)* dy
!
!             xcorner = xlowb + (dfloat(i)-1.)* dx
!             ycorner = ylowb + (dfloat(j)-1.)* dy
!             do itimes = 1,4
!                if (itimes .eq. 1) then ! get all 4 corners of mesh, pw constant sol
!                   xc = xcorner
!                   yc = ycorner
!                else if (itimes .eq. 2) then
!                   xc = xcorner
!                   yc = ycorner + dy
!                else if (itimes .eq. 3) then
!                   xc = xcorner + dx
!                   yc = ycorner + dy
!                else if (itimes .eq. 4) then
!                   xc = xcorner + dx
!                   yc = ycorner
!                endif
!
!            corners(itimes,1) = xc
!            corners(itimes,2) = yc
!
!              end do
!
!
!                call computed_error(i,j,dx,dy,xlow,ylow,
!     .          time, irr,
!     .          q, qx, qy, qxx, qxy, qyy,
!     .          qxxx, qxxy, qxyy, qyyy,
!     .          mitot, mjtot, nvar,
!     .          ivar, iprob, lstgrd,lwidth,
!     .          douterr, doutmaxerr, doutemag, doutmaxemag,
!     .          douterr2, doutmaxerr2, doutemag2, doutmaxemag2)
!
!       write(231,*)  'SQ ( ', corners(1,1),",",
!     .                        corners(1,2),",0,",
!     .                        corners(2,1),",",
!     .                        corners(2,2),",0,",
!     .                        corners(3,1),",",
!     .                        corners(3,2),",0,",
!     .                        corners(4,1),",",
!     .                        corners(4,2),",0)",
!     .                        "{",
!     .                        douterr,",",
!     .                        douterr,",",
!     .                        douterr,",",
!     .                        douterr,
!     .                        "};"
!
!
!         else
!            ivert = 1
!            do 9290 while (poly(ivert+1,1,kirr) .ne. -11.)
!               ivert = ivert + 1
! 9290       continue
!            itri = ivert - 3
!
!                  idx1 = 1
!                  do 991 it = 1, itri ! for each  triangle
!                    idx2 = it + 1
!                    idx3 = it + 2
!
!                    x1 = poly(idx1,1,kirr)
!                    y1 = poly(idx1,2,kirr)
!
!                    x2 = poly(idx2,1,kirr)
!                    y2 = poly(idx2,2,kirr)
!
!                    x3 = poly(idx3,1,kirr)
!                    y3 = poly(idx3,2,kirr)
!
!
!                xcen = xcirr(kirr)
!                ycen = ycirr(kirr)
!
!                call computed_error(i,j,dx,dy,xlow,ylow,
!     .          time, irr,
!     .          q, qx, qy, qxx, qxy, qyy,
!     .          qxxx, qxxy, qxyy, qyyy,
!     .          mitot, mjtot, nvar,
!     .          ivar, iprob, lstgrd,lwidth,
!     .          douterr, doutmaxerr, doutemag, doutmaxemag,
!     .          douterr2, doutmaxerr2, doutemag2, doutmaxemag2)
!
!       write(231,*)  'ST ( ', x1,",",
!     .                        y1,",0,",
!     .                        x2,",",
!     .                        y2,",0,",
!     .                        x3,",",
!     .                        y3,",0)",
!     .                        "{",
!     .                        douterr,",",
!     .                        douterr,",",
!     .                        douterr,
!     .                        "};"
!991   continue
!
!
!
!         endif
!915   continue
!      write(231,*) "};"
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!      write(231,*) 'View "numhoods" {'
!      do 215 i = lwidth+1 - ibuff, mitot-lwidth+ibuff
!      do 215 j = lwidth+1-ibuff, mjtot-lwidth+ibuff
!         kirr = irr(i,j)
!         if (kirr .eq. -1) go to 215
!         if (kirr .eq. lstgrd ) then
!            xcen = xlowb + (dfloat(i)-.5d0)* dx
!            ycen = ylowb + (dfloat(j)-.5d0)* dy
!
!             xcorner = xlowb + (dfloat(i)-1.)* dx
!             ycorner = ylowb + (dfloat(j)-1.)* dy
!             do itimes = 1,4
!                if (itimes .eq. 1) then ! get all 4 corners of mesh, pw constant sol
!                   xc = xcorner
!                   yc = ycorner
!                else if (itimes .eq. 2) then
!                   xc = xcorner
!                   yc = ycorner + dy
!                else if (itimes .eq. 3) then
!                   xc = xcorner + dx
!                   yc = ycorner + dy
!                else if (itimes .eq. 4) then
!                   xc = xcorner + dx
!                   yc = ycorner
!                endif
!
!            corners(itimes,1) = xc
!            corners(itimes,2) = yc
!
!            recon(:,itimes) = q(i,j,:)
!     .   +(xc-xcen)*qx(i,j,:)
!     .   +(yc-ycen)*qy(i,j,:)
!     .   + 0.5d0*qxx(i,j,:)*( (xc-xcen)**2-(dx**2)/12.d0)
!     .   + (xc-xcen)*(yc-ycen)*qxy(i,j,:)
!     .   + 0.5d0*qyy(i,j,:)*( (yc-ycen)**2-(dy**2)/12.d0)
!              end do
!
!       write(231,*)  'SQ ( ', corners(1,1),",",
!     .                        corners(1,2),",0,",
!     .                        corners(2,1),",",
!     .                        corners(2,2),",0,",
!     .                        corners(3,1),",",
!     .                        corners(3,2),",0,",
!     .                        corners(4,1),",",
!     .                        corners(4,2),",0)",
!     .                        "{",
!     .                        numhoods(i,j),",",
!     .                        numhoods(i,j),",",
!     .                        numhoods(i,j),",",
!     .                        numhoods(i,j),
!     .                        "};"
!
!
!         else
!            ivert = 1
!            do 2290 while (poly(ivert+1,1,kirr) .ne. -11.)
!               ivert = ivert + 1
! 2290       continue
!            itri = ivert - 3
!
!                  idx1 = 1
!                  do 2291 it = 1, itri ! for each  triangle
!                    idx2 = it + 1
!                    idx3 = it + 2
!
!                    x1 = poly(idx1,1,kirr)
!                    y1 = poly(idx1,2,kirr)
!
!                    x2 = poly(idx2,1,kirr)
!                    y2 = poly(idx2,2,kirr)
!
!                    x3 = poly(idx3,1,kirr)
!                    y3 = poly(idx3,2,kirr)
!
!
!                xcen = xcirr(kirr)
!                ycen = ycirr(kirr)
!
!            xc = x1
!            yc = y1
!            recon(:,1) = q(i,j,:)
!     .   +(xc-xcen)*qx(i,j,:)
!     .   +(yc-ycen)*qy(i,j,:)
!     .   + 0.5d0*qxx(i,j,:)*( (xc-xcen)**2-(dx**2)/12.d0)
!     .   + (xc-xcen)*(yc-ycen)*qxy(i,j,:)
!     .   + 0.5d0*qyy(i,j,:)*( (yc-ycen)**2-(dy**2)/12.d0)
!
!            xc = x2
!            yc = y2
!            recon(:,2) = q(i,j,:)
!     .   +(xc-xcen)*qx(i,j,:)
!     .   +(yc-ycen)*qy(i,j,:)
!     .   + 0.5d0*qxx(i,j,:)*( (xc-xcen)**2-(dx**2)/12.d0)
!     .   + (xc-xcen)*(yc-ycen)*qxy(i,j,:)
!     .   + 0.5d0*qyy(i,j,:)*( (yc-ycen)**2-(dy**2)/12.d0)
!
!            xc = x3
!            yc = y3
!            recon(:,3) = q(i,j,:)
!     .   +(xc-xcen)*qx(i,j,:)
!     .   +(yc-ycen)*qy(i,j,:)
!     .   + 0.5d0*qxx(i,j,:)*( (xc-xcen)**2-(dx**2)/12.d0)
!     .   + (xc-xcen)*(yc-ycen)*qxy(i,j,:)
!     .   + 0.5d0*qyy(i,j,:)*( (yc-ycen)**2-(dy**2)/12.d0)
!
!       write(231,*)  'ST ( ', x1,",",
!     .                        y1,",0,",
!     .                        x2,",",
!     .                        y2,",0,",
!     .                        x3,",",
!     .                        y3,",0)",
!     .                        "{",
!     .                        numhoods(i,j),",",
!     .                        numhoods(i,j),",",
!     .                        numhoods(i,j),
!     .                        "};"
!2291   continue
!
!
!
!
!         endif
! 215   continue
!      write(231,*) "};"
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!      write(231,*) 'View "ncount" {'
!      do 2159 i = lwidth+1 - ibuff, mitot-lwidth+ibuff
!      do 2159 j = lwidth+1-ibuff, mjtot-lwidth+ibuff
!         kirr = irr(i,j)
!         if (kirr .eq. -1) go to 2159
!         if (kirr .eq. lstgrd ) then
!            xcen = xlowb + (dfloat(i)-.5d0)* dx
!            ycen = ylowb + (dfloat(j)-.5d0)* dy
!
!             xcorner = xlowb + (dfloat(i)-1.)* dx
!             ycorner = ylowb + (dfloat(j)-1.)* dy
!             do itimes = 1,4
!                if (itimes .eq. 1) then ! get all 4 corners of mesh, pw constant sol
!                   xc = xcorner
!                   yc = ycorner
!                else if (itimes .eq. 2) then
!                   xc = xcorner
!                   yc = ycorner + dy
!                else if (itimes .eq. 3) then
!                   xc = xcorner + dx
!                   yc = ycorner + dy
!                else if (itimes .eq. 4) then
!                   xc = xcorner + dx
!                   yc = ycorner
!                endif
!
!            corners(itimes,1) = xc
!            corners(itimes,2) = yc
!
!            recon(:,itimes) = q(i,j,:)
!     .   +(xc-xcen)*qx(i,j,:)
!     .   +(yc-ycen)*qy(i,j,:)
!     .   + 0.5d0*qxx(i,j,:)*( (xc-xcen)**2-(dx**2)/12.d0)
!     .   + (xc-xcen)*(yc-ycen)*qxy(i,j,:)
!     .   + 0.5d0*qyy(i,j,:)*( (yc-ycen)**2-(dy**2)/12.d0)
!              end do
!
!       write(231,*)  'SQ ( ', corners(1,1),",",
!     .                        corners(1,2),",0,",
!     .                        corners(2,1),",",
!     .                        corners(2,2),",0,",
!     .                        corners(3,1),",",
!     .                        corners(3,2),",0,",
!     .                        corners(4,1),",",
!     .                        corners(4,2),",0)",
!     .                        "{",
!     .                        0,",",
!     .                        0,",",
!     .                        0,",",
!     .                        0,
!     .                        "};"
!
!
!         else
!            ivert = 1
!            do 22909 while (poly(ivert+1,1,kirr) .ne. -11.)
!               ivert = ivert + 1
!22909       continue
!            itri = ivert - 3
!
!                  idx1 = 1
!                  do 2299 it = 1, itri ! for each  triangle
!                    idx2 = it + 1
!                    idx3 = it + 2
!
!                    x1 = poly(idx1,1,kirr)
!                    y1 = poly(idx1,2,kirr)
!
!                    x2 = poly(idx2,1,kirr)
!                    y2 = poly(idx2,2,kirr)
!
!                    x3 = poly(idx3,1,kirr)
!                    y3 = poly(idx3,2,kirr)
!
!
!                xcen = xcirr(kirr)
!                ycen = ycirr(kirr)
!
!            xc = x1
!            yc = y1
!            recon(:,1) = q(i,j,:)
!     .   +(xc-xcen)*qx(i,j,:)
!     .   +(yc-ycen)*qy(i,j,:)
!     .   + 0.5d0*qxx(i,j,:)*( (xc-xcen)**2-(dx**2)/12.d0)
!     .   + (xc-xcen)*(yc-ycen)*qxy(i,j,:)
!     .   + 0.5d0*qyy(i,j,:)*( (yc-ycen)**2-(dy**2)/12.d0)
!
!            xc = x2
!            yc = y2
!            recon(:,2) = q(i,j,:)
!     .   +(xc-xcen)*qx(i,j,:)
!     .   +(yc-ycen)*qy(i,j,:)
!     .   + 0.5d0*qxx(i,j,:)*( (xc-xcen)**2-(dx**2)/12.d0)
!     .   + (xc-xcen)*(yc-ycen)*qxy(i,j,:)
!     .   + 0.5d0*qyy(i,j,:)*( (yc-ycen)**2-(dy**2)/12.d0)
!
!            xc = x3
!            yc = y3
!            recon(:,3) = q(i,j,:)
!     .   +(xc-xcen)*qx(i,j,:)
!     .   +(yc-ycen)*qy(i,j,:)
!     .   + 0.5d0*qxx(i,j,:)*( (xc-xcen)**2-(dx**2)/12.d0)
!     .   + (xc-xcen)*(yc-ycen)*qxy(i,j,:)
!     .   + 0.5d0*qyy(i,j,:)*( (yc-ycen)**2-(dy**2)/12.d0)
!
!       write(231,*)  'ST ( ', x1,",",
!     .                        y1,",0,",
!     .                        x2,",",
!     .                        y2,",0,",
!     .                        x3,",",
!     .                        y3,",0)",
!     .                        "{",
!     .                        ncount(kirr),",",
!     .                        ncount(kirr),",",
!     .                        ncount(kirr),
!     .                        "};"
!2299   continue
!
!
!
!
!         endif
! 2159   continue
!      write(231,*) "};"

      close(231)

!      if(iprob .eq. 11) then
!      call p11plot(q, qx, qy,lwidth, mitot, mjtot, nvar, lstgrd, irr)
!      endif

      return
      end


text:
      SUBROUTINE SB04MY( N, M, IND, A, LDA, B, LDB, C, LDC, D, IPR,
     $                   INFO )
C
C     SLICOT RELEASE 5.7.
C
C     Copyright (c) 2002-2020 NICONET e.V.
C
C     PURPOSE
C
C     To construct and solve a linear algebraic system of order M whose
C     coefficient matrix is in upper Hessenberg form. Such systems
C     appear when solving Sylvester equations using the Hessenberg-Schur
C     method.
C
C     ARGUMENTS
C
C     Input/Output Parameters
C
C     N       (input) INTEGER
C             The order of the matrix B.  N >= 0.
C
C     M       (input) INTEGER
C             The order of the matrix A.  M >= 0.
C
C     IND     (input) INTEGER
C             The index of the column in C to be computed.  IND >= 1.
C
C     A       (input) DOUBLE PRECISION array, dimension (LDA,M)
C             The leading M-by-M part of this array must contain an
C             upper Hessenberg matrix.
C
C     LDA     INTEGER
C             The leading dimension of array A.  LDA >= MAX(1,M).
C
C     B       (input) DOUBLE PRECISION array, dimension (LDB,N)
C             The leading N-by-N part of this array must contain a
C             matrix in real Schur form.
C
C     LDB     INTEGER
C             The leading dimension of array B.  LDB >= MAX(1,N).
C
C     C       (input/output) DOUBLE PRECISION array, dimension (LDC,N)
C             On entry, the leading M-by-N part of this array must
C             contain the coefficient matrix C of the equation.
C             On exit, the leading M-by-N part of this array contains
C             the matrix C with column IND updated.
C
C     LDC     INTEGER
C             The leading dimension of array C.  LDC >= MAX(1,M).
C
C     Workspace
C
C     D       DOUBLE PRECISION array, dimension (M*(M+1)/2+2*M)
C
C     IPR     INTEGER array, dimension (2*M)
C
C     Error Indicator
C
C     INFO    INTEGER
C             = 0:  successful exit;
C             > 0:  if INFO = IND, a singular matrix was encountered.
C
C     METHOD
C
C     A special linear algebraic system of order M, with coefficient
C     matrix in upper Hessenberg form is constructed and solved. The
C     coefficient matrix is stored compactly, row-wise.
C
C     REFERENCES
C
C     [1] Golub, G.H., Nash, S. and Van Loan, C.F.
C         A Hessenberg-Schur method for the problem AX + XB = C.
C         IEEE Trans. Auto. Contr., AC-24, pp. 909-913, 1979.
C
C     NUMERICAL ASPECTS
C
C     None.
C
C     CONTRIBUTORS
C
C     Release 3.0: V. Sima, Katholieke Univ. Leuven, Belgium, Sep. 1997.
C     Supersedes Release 2.0 routine SB04AY by G. Golub, S. Nash, and
C     C. Van Loan, Stanford University, California, United States of
C     America, January 1982.
C
C     REVISIONS
C
C     -
C
C     KEYWORDS
C
C     Hessenberg form, orthogonal transformation, real Schur form,
C     Sylvester equation.
C
C     ******************************************************************
C
C     .. Scalar Arguments ..
      INTEGER           INFO, IND, LDA, LDB, LDC, M, N
C     .. Array Arguments ..
      INTEGER           IPR(*)
      DOUBLE PRECISION  A(LDA,*), B(LDB,*), C(LDC,*), D(*)
C     .. Local Scalars ..
      INTEGER           I, I2, J, K, K1, K2, M1
C     .. External Subroutines ..
      EXTERNAL          DAXPY, DCOPY, SB04MW
C     .. Intrinsic Functions ..
      INTRINSIC         MAX
C     .. Executable Statements ..
C
      DO 20 I = IND + 1, N
         CALL DAXPY( M, -B(IND,I), C(1,I), 1, C(1,IND), 1 )
   20 CONTINUE
C
      M1 = M + 1
      I2 = ( M*M1 )/2 + M1
      K2 = 1
      K  = M
C
C     Construct the linear algebraic system of order M.
C
      DO 40 I = 1, M
         J = M1 - K
         CALL DCOPY ( K, A(I,J), LDA, D(K2), 1 )
         K1 = K2
         K2 = K2 + K
         IF ( I.GT.1 ) THEN
            K1 = K1 + 1
            K  = K - 1
         END IF
         D(K1) = D(K1) + B(IND,IND)
C
C        Store the right hand side.
C
         D(I2) = C(I,IND)
         I2 = I2 + 1
   40 CONTINUE
C
C     Solve the linear algebraic system and store the solution in C.
C
      CALL SB04MW( M, D, IPR, INFO )
C
      IF ( INFO.NE.0 ) THEN
         INFO = IND
      ELSE
C
         DO 60 I = 1, M
            C(I,IND) = D(IPR(I))
   60    CONTINUE
C
      END IF
C
      RETURN
C *** Last line of SB04MY ***
      END

ext: f

max_stars_repo_name: bnavigator/SLICOT-Reference

max_issues_repo_name: RJHKnight/slicotr

max_forks_repo_name: RJHKnight/slicotr


text:
program list06p

! CAP372 - exercise 06 - parallel version - 2019-09-21
! integration of pi : 4.0/(1+x*) dx, interval: 0 to hs
! module load intel_psxe/2019
! ifort -g -check all -fpe0 -warn -traceback -debug extended -qopenmp &
!       -o ilist06p list06p.f90
! ./ilist06s
! or if using gnu:
! gfortran -Og -Wall -fcheck=all -fopenmp -o list06p list06p.f90
! ./list06s

  use omp_lib
  implicit none

  integer, parameter :: n=2**26, t=2   ! 1
!  integer, parameter :: n=2**26, t=12  ! 2
!  integer, parameter :: n=2**26, t=24  ! 3
!  integer, parameter :: n=2**28, t=2   ! 4
!  integer, parameter :: n=2**28, t=12  ! 5
!  integer, parameter :: n=2**28, t=24  ! 6
!  integer, parameter :: n=2**30, t=2   ! 7
!  integer, parameter :: n=2**30, t=12  ! 8
!  integer, parameter :: n=2**30, t=24  ! 9
  double precision, parameter  :: a=0.0, b=1.0
  double precision :: pi, t1, t2, x, h, integral=0.0
  integer :: i

  call cpu_time(t1)
  integral = ( f(a) + f(b) ) / 2.0
  x = a
  h = ( b - a ) / n

  !$omp parallel do private(x) reduction(+:integral) num_threads(t)
  do i = 1, n -1
    x = a + h * i
    integral = integral + 4.0 / (1.0 + x * x)
  end do
  !$omp end parallel do

  pi = integral * h
  call cpu_time(t2)
  print*, "Result:", pi, " Error:", dacos(-1.d0) - pi
  print*, "Partitions:", n,"    Threads:", t, "    Elapsed time:", t2 - t1

contains

  function f(x)
    implicit none

    double precision :: f
    double precision, intent(in) :: x

    f = ( 4.0 / ( 1.0 + x * x ) )

  end function f

end program list06p

max_stars_repo_name: efurlanm/hpc

max_issues_repo_name: efurlanm/hpc

max_forks_repo_name: efurlanm/hpc


text:
program doloop_02
implicit none
integer :: i, j, a, b
j = 0
a = 1
b = 10
do i = a, b
    j = j + i
end do
if (j /= 55) error stop

a = 0
do i = 1, 10
    do j = 1, 10
        a = a + (i-1)*10+j
    end do
end do
if (a /= 100*101/2) error stop

a = 0
do i = 1, 10
    do j = 1, i
        a = a + j
    end do
end do
if (a /= 220) error stop
end

max_stars_repo_name: Thirumalai-Shaktivel/lfortran

max_issues_repo_name: Thirumalai-Shaktivel/lfortran

max_forks_repo_name: Thirumalai-Shaktivel/lfortran


text:
!----------------------------------------------------------------------------
module nemsio_read
!
!$$$ documentation clock
!
! module: nemsio_read      read data fields from a nemsio file
!  Programmer: J. Wang          date: 2011-01-13
!
! abstract: this module provides subroutines to read data fields out of a
!           nemsio file. The data file could be 'bin4','bin8', or 'grib',
!           The dat field could be read out by the record number
!           or by given the data field name, level type and level. Overload
!           interfaces are provided to handle different data type (real(4)
!           or real(8) ) of the array that holds the data.
!
! Possible return code
!          0   Successful call
!         -31  get dimension ffrom gfile
!         -32  read data field by record number using w3d
!         -33  read data field by given data field name,level type and level using w3d
!         -34  read data field by record number using w34
!         -35  read data field by given data field name,level type and level using w34
!         -41  read data field by record number from 4 byte real binary file
!         -42  read data field by record number from 8 byte real binary file
!         -43  read data field by field name,levtyp and lev from 4 byte real binary file
!         -44  read data field by field name,levtyp and lev from 8 byte real binary file
!         -45  read data field by record number using w34 from grib data
!         -46  read data field by field name,level type and level using w34 from grib data
!         -47  read data field by record number using w3d from grib data
!         -48  read data field by field name,level type and level using w3d from grib data
!------------------------------------------------------------------------------
!
  use nemsio_openclose
!
  implicit none
!
  private
!------------------------------------------------------------------------------
!----- interface
!
  interface nemsio_readrec
    module procedure nemsio_readrec4
    module procedure nemsio_readrec8
  end interface nemsio_readrec
!
  interface nemsio_readrecv
    module procedure nemsio_readrecv4
    module procedure nemsio_readrecv8
  end interface nemsio_readrecv
!
  interface nemsio_readrecw34
    module procedure nemsio_readrec4w34
    module procedure nemsio_readrec8w34
  end interface nemsio_readrecw34
!
  interface nemsio_readrecvw34
    module procedure nemsio_readrecv4w34
    module procedure nemsio_readrecv8w34
  end interface nemsio_readrecvw34
!
!public mehtods
  public nemsio_readrec,nemsio_readrecv,nemsio_readrecw34,nemsio_readrecvw34
!
!---------------------------------------------------------
! local data
!
  character(8) :: mygdatatype
  character(255) :: mygfname
  integer mydimx,mydimy,mydimz,mynframe,myfieldsize,mytlmeta,myflunit
  character(255),save :: mygfnamep=''
  integer,save :: mymbuf,mynnum,mynlen,mymnum
  character,allocatable,save  :: mycbuf(:)
  logical do_byteswap
!
contains
!
!------------------------------------------------------------------------------
  subroutine nemsio_getgfile(gfile,iret)
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - - - - - -
! abstract: read nemsio data by record number into a 2D 32 bits array
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - - - - - -
    implicit none
    type(nemsio_gfile),intent(in)                 :: gfile
    integer(nemsio_intkind),optional,intent(out)  :: iret
    integer ios
    character(8) :: tmpgdatatype
!
    if(present(iret)) iret= -31
!
    call nemsio_getfilehead(gfile,iret=ios,gdatatype=tmpgdatatype,dimx=mydimx,   &
           dimy=mydimy,dimz=mydimz,nframe=mynframe,tlmeta=mytlmeta,              &
           flunit=myflunit,gfname=mygfname,do_byteswap=do_byteswap )
    if(ios/=0) then
       if(present(iret)) then
         iret=ios
         return
       else
         print *,'ERROR: NEMSIO readrec in getting file head'
         stop
       endif
    endif

    myfieldsize=(mydimx+2*mynframe)*(mydimy+2*mynframe)
    mygdatatype=tmpgdatatype(1:4)
    if(trim(mygfnamep)/=trim(mygfname)) then
       mygfnamep=mygfname
       if(trim(mygdatatype)=='grib') then
         mymbuf=256*1024
         mynnum=0
         mynlen=0
         mymnum=-1
         if(allocated(mycbuf)) deallocate(mycbuf)
         allocate(mycbuf(mymbuf))
       endif
     endif
!     print *,'in read,mygdatatype=',mygdatatype,'do_byteswap=',do_byteswap
!
     if(present(iret)) iret=0
!
  end subroutine nemsio_getgfile
!------------------------------------------------------------------------------
!
!------------------------------------------------------------------------------
  subroutine nemsio_readrec4(gfile,jrec,data,nframe,iret)
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - - - - - -
! abstract: read nemsio data by record number into a 2D 32 bits array
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - - - - - -
    implicit none
    type(nemsio_gfile),intent(inout)              :: gfile
    integer(nemsio_intkind),intent(in)            :: jrec
    real(nemsio_realkind),intent(inout)           :: data(:)
    integer(nemsio_intkind),optional,intent(out)  :: iret
    integer(nemsio_intkind),optional,intent(in)  :: nframe
    real(nemsio_realkind),allocatable             :: datatmp(:)
    real(nemsio_dblekind),allocatable            :: datatmp8(:)
    integer :: i,j,ios
!------------------------------------------------------------
! read 4 byte rec
!------------------------------------------------------------
   if(present(iret)) iret=-32
!---
   call nemsio_getgfile(gfile,iret)
!---
   if ( mygdatatype .eq. 'bin4') then
     if(.not.present(nframe) ) then
       call nemsio_readrecbin4d4(gfile,jrec,data,ios)
     else
      allocate(datatmp(myfieldsize) )
      call nemsio_readrecbin4d4(gfile,jrec,datatmp,ios)
     endif
   else if ( mygdatatype .eq. 'bin8') then
     allocate(datatmp8(myfieldsize) )
     call nemsio_readrecbin8d8(gfile,jrec,datatmp8,ios)
   else
     allocate(datatmp8(myfieldsize) )
     call nemsio_readrecgrb8(gfile,jrec,datatmp8,ios)
   endif
   if ( ios .ne.0 ) then
     if(present(iret)) then
       iret=ios
       return
     else
       call nemsio_stop
     endif
   endif
!---
   if ( present(nframe) ) then
     if(mygdatatype .eq. 'bin4') then
       do j=1,mydimy+2*mynframe-2*nframe
        do i=1,mydimx+2*mynframe -2*nframe
         data(i+(j-1)*(mydimx+2*mynframe-2*nframe))=datatmp(i+nframe        &
           +(j-1+nframe)*(mydimx+2*mynframe))
        enddo
       enddo
       deallocate(datatmp)
     elseif(mygdatatype=='bin8'.or.mygdatatype=='grib') then
       do j=1,mydimy+2*mynframe-2*nframe
        do i=1,mydimx+2*mynframe -2*nframe
         data(i+(j-1)*(mydimx+2*mynframe-2*nframe))=datatmp8(i+nframe        &
           +(j-1+nframe)*(mydimx+2*mynframe))
        enddo
       enddo
       deallocate(datatmp8)
     endif
   else
     if(mygdatatype=='bin8'.or.mygdatatype=='grib') then
       data=datatmp8
       deallocate(datatmp8)
     endif
   endif
!---
   if(present(iret)) iret=0
   return
  end subroutine nemsio_readrec4
!------------------------------------------------------------------------------
  subroutine nemsio_readrec8(gfile,jrec,data,nframe,iret)
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - - - - - -
! abstract: read nemsio data (bin) by record number into a 2D 32 bits array
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - - - - - -
    implicit none
    type(nemsio_gfile),intent(inout)              :: gfile
    integer(nemsio_intkind),intent(in)            :: jrec
    real(nemsio_dblekind),intent(inout)           :: data(:)
    integer(nemsio_intkind),optional,intent(out)  :: iret
    integer(nemsio_intkind),optional,intent(in)   :: nframe
    real(nemsio_realkind),allocatable             :: datatmp4(:)
    real(nemsio_dblekind),allocatable             :: datatmp(:)
    integer :: i,j,ios
!------------------------------------------------------------
! read 4 byte rec
!------------------------------------------------------------
   if(present(iret)) iret=-32
!---
   call nemsio_getgfile(gfile,iret)
!
   if ( mygdatatype .eq. 'bin4') then
     allocate(datatmp4(myfieldsize))
     call nemsio_readrecbin4d4(gfile,jrec,datatmp4,ios)
   else if ( mygdatatype .eq. 'bin8') then
     if(.not.present(nframe)) then
       call nemsio_readrecbin8d8(gfile,jrec,data,ios)
     else
       allocate(datatmp(myfieldsize))
       call nemsio_readrecbin8d8(gfile,jrec,datatmp,ios)
     endif
   else
     if(.not.present(nframe)) then
       call nemsio_readrecgrb8(gfile,jrec,data,ios)
     else
       allocate(datatmp(myfieldsize))
       call nemsio_readrecgrb8(gfile,jrec,datatmp,ios)
     endif
   endif
   if ( ios .ne.0 ) then
     if(present(iret)) then
       iret=ios
       return
     else
       call nemsio_stop
     endif
   endif
!---
   if ( present(nframe) ) then
     if(mygdatatype=='bin4') then
        do j=1,mydimy+2*mynframe-2*nframe
         do i=1,mydimx+2*mynframe -2*nframe
          data(i+(j-1)*(mydimx+2*mynframe-2*nframe))=datatmp4(i+nframe        &
            +(j-1+nframe)*(mydimx+2*mynframe))
         enddo
        enddo
        deallocate(datatmp4)
     elseif(mygdatatype=='bin8'.or.mygdatatype=='grib') then
        do j=1,mydimy+2*mynframe-2*nframe
         do i=1,mydimx+2*mynframe -2*nframe
          data(i+(j-1)*(mydimx+2*mynframe-2*nframe))=datatmp(i+nframe        &
            +(j-1+nframe)*(mydimx+2*mynframe))
         enddo
        enddo
        deallocate(datatmp)
     endif
   else
     if(mygdatatype=='bin4') then
       data=datatmp4
       deallocate(datatmp4)
     endif
   endif
!
   if(present(iret)) iret=0
   return
  end subroutine nemsio_readrec8
!------------------------------------------------------------------------------
  subroutine nemsio_readrecv4(gfile,name,levtyp,lev,data,nframe,iret)
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - - - - - -
! abstract: read nemsio data by record number into a 2D 32 bits array
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - - - - - -
    implicit none
    type(nemsio_gfile),intent(inout)              :: gfile
    character(*),intent(in)                       :: name
    character(*),intent(in),optional              :: levtyp
    integer(nemsio_intkind),optional,intent(in)   :: lev
    real(nemsio_realkind),intent(inout)           :: data(:)
    integer(nemsio_intkind),optional,intent(out)  :: iret
    integer(nemsio_intkind),optional,intent(in)   :: nframe
    real(nemsio_realkind),allocatable             :: datatmp(:)
    real(nemsio_dblekind),allocatable             :: datatmp8(:)
    integer :: i,j,ios
!------------------------------------------------------------
! read 4 byte rec
!------------------------------------------------------------
   if(present(iret)) iret=-33
!---
   call nemsio_getgfile(gfile,iret)
!---
   if ( mygdatatype .eq. 'bin4') then
     if(.not.present(nframe) ) then
       call nemsio_readrecvbin4d4(gfile,name,levtyp,lev,data,ios)
     else
       allocate(datatmp(myfieldsize) )
       call nemsio_readrecvbin4d4(gfile,name,levtyp,lev,datatmp,ios)
     endif
   else if ( mygdatatype .eq. 'bin8') then
     allocate(datatmp8(myfieldsize) )
     call nemsio_readrecvbin8d8(gfile,name,levtyp,lev,datatmp8,ios)
   else
     allocate(datatmp8(myfieldsize) )
     call nemsio_readrecvgrb8(gfile,name,levtyp,lev,datatmp8,ios)
   endif
   if ( ios .ne.0 ) then
     if(present(iret)) then
       iret=ios
       return
     else
       call nemsio_stop
     endif
   endif
!---
   if ( present(nframe) ) then
     if(mygdatatype=='bin4') then
        do j=1,mydimy+2*mynframe-2*nframe
         do i=1,mydimx+2*mynframe -2*nframe
          data(i+(j-1)*(mydimx+2*mynframe-2*nframe))=datatmp(i+nframe        &
            +(j-1+nframe)*(mydimx+2*mynframe))
         enddo
        enddo
        deallocate(datatmp)
     elseif(mygdatatype=='bin8'.or.mygdatatype=='grib' ) then
        do j=1,mydimy+2*mynframe-2*nframe
         do i=1,mydimx+2*mynframe -2*nframe
          data(i+(j-1)*(mydimx+2*mynframe-2*nframe))=datatmp8(i+nframe        &
            +(j-1+nframe)*(mydimx+2*mynframe))
         enddo
        enddo
        deallocate(datatmp8)
     endif
   else
     if(mygdatatype=='bin8'.or.mygdatatype=='grib' ) then
       data=datatmp8
       deallocate(datatmp8)
     endif
   endif
!---
   if(present(iret)) iret=0
   return
  end subroutine nemsio_readrecv4
!------------------------------------------------------------------------------
  subroutine nemsio_readrecv8(gfile,name,levtyp,lev,data,nframe,iret)
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - - - - - -
! abstract: read nemsio data by record number into a 2D 32 bits array
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - - - - - -
    implicit none
    type(nemsio_gfile),intent(inout)              :: gfile
    character(*),intent(in)                       :: name
    character(*),intent(in),optional              :: levtyp
    integer(nemsio_intkind),optional,intent(in)   :: lev
    real(nemsio_dblekind),intent(inout)           :: data(:)
    integer(nemsio_intkind),optional,intent(out)  :: iret
    integer(nemsio_intkind),optional,intent(in)   :: nframe
    real(nemsio_realkind),allocatable             :: datatmp4(:)
    real(nemsio_dblekind),allocatable             :: datatmp(:)
    integer :: i,j,ios
!------------------------------------------------------------
! read 8 byte rec
!------------------------------------------------------------
   if(present(iret)) iret=-33
!---
   call nemsio_getgfile(gfile,iret)
!---
   if ( mygdatatype .eq. 'bin4') then
     allocate(datatmp4(myfieldsize) )
     call nemsio_readrecvbin4d4(gfile,name,levtyp,lev,datatmp4,ios)
   else if ( mygdatatype .eq. 'bin8') then
     if(.not.present(nframe) ) then
       call nemsio_readrecvbin8d8(gfile,name,levtyp,lev,data,ios)
     else
       allocate(datatmp(myfieldsize) )
       call nemsio_readrecvbin8d8(gfile,name,levtyp,lev,datatmp,ios)
     endif
   else
     if(.not.present(nframe) ) then
       call nemsio_readrecvgrb8(gfile,name,levtyp,lev,data,ios)
     else
       allocate(datatmp(myfieldsize) )
       call nemsio_readrecvgrb8(gfile,name,levtyp,lev,datatmp,ios)
     endif
   endif
   if ( ios .ne.0 ) then
     if(present(iret)) then
       iret=ios
       return
     else
       call nemsio_stop
     endif
   endif
!---
   if ( present(nframe) ) then
      if(mygdatatype=='bin4') then
        do j=1,mydimy+2*mynframe-2*nframe
         do i=1,mydimx+2*mynframe -2*nframe
          data(i+(j-1)*(mydimx+2*mynframe-2*nframe))=datatmp4(i+nframe        &
           +(j-1+nframe)*(mydimx+2*mynframe))
         enddo
        enddo
        deallocate(datatmp4)
      elseif(mygdatatype=='bin8'.or.mygdatatype=='grib') then
        do j=1,mydimy+2*mynframe-2*nframe
         do i=1,mydimx+2*mynframe -2*nframe
          data(i+(j-1)*(mydimx+2*mynframe-2*nframe))=datatmp(i+nframe        &
           +(j-1+nframe)*(mydimx+2*mynframe))
         enddo
        enddo
        deallocate(datatmp)
      endif
   else
     if(mygdatatype=='bin4') then
       data=datatmp4
       deallocate(datatmp4)
     endif
   endif
!
   if(present(iret)) iret=0
   return
  end subroutine nemsio_readrecv8
!------------------------------------------------------------------------------
  subroutine nemsio_readrec4w34(gfile,jrec,data,nframe,iret)
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - - - - - -
! abstract: read nemsio data by record number into a 2D 32 bits array,
!           using w3_4 library to compile
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - - - - - -
    implicit none
    type(nemsio_gfile),intent(inout)              :: gfile
    integer(nemsio_intkind),intent(in)            :: jrec
    real(nemsio_realkind),intent(inout)           :: data(:)
    integer(nemsio_intkind),optional,intent(out)  :: iret
    integer(nemsio_intkind),optional,intent(in)   :: nframe
    real(nemsio_realkind),allocatable             :: datatmp(:)
    real(nemsio_dblekind),allocatable             :: datatmp8(:)
    integer :: i,j,ios
!------------------------------------------------------------
! read 4 byte rec
!------------------------------------------------------------
   if(present(iret)) iret=-34
!---
   call nemsio_getgfile(gfile,iret)
!---
   if ( mygdatatype .eq. 'bin4') then
     if(.not.present(nframe)) then
       call nemsio_readrecbin4d4(gfile,jrec,data,ios)
     else
       allocate(datatmp(myfieldsize) )
       call nemsio_readrecbin4d4(gfile,jrec,datatmp,ios)
     endif
   else if ( mygdatatype .eq. 'bin8') then
      allocate(datatmp8(myfieldsize) )
      call nemsio_readrecbin8d8(gfile,jrec,datatmp8,ios)
   else
     if(.not.present(nframe)) then
       call nemsio_readrecgrb4w34(gfile,jrec,data,ios)
     else
       allocate(datatmp(myfieldsize) )
       call nemsio_readrecgrb4w34(gfile,jrec,datatmp,ios)
     endif
   endif
   if ( ios .ne.0 ) then
     if(present(iret)) then
       iret=ios
       return
     else
       call nemsio_stop
     endif
   endif
!---
   if ( present(nframe) ) then
      if(mygdatatype=='bin4'.or.mygdatatype=='grib') then
        do j=1,mydimy+2*mynframe-2*nframe
         do i=1,mydimx+2*mynframe -2*nframe
          data(i+(j-1)*(mydimx+2*mynframe-2*nframe))=datatmp(i+nframe        &
            +(j-1+nframe)*(mydimx+2*mynframe))
         enddo
        enddo
        deallocate(datatmp)
      elseif(mygdatatype=='bin8') then
        do j=1,mydimy+2*mynframe-2*nframe
         do i=1,mydimx+2*mynframe -2*nframe
          data(i+(j-1)*(mydimx+2*mynframe-2*nframe))=datatmp8(i+nframe        &
            +(j-1+nframe)*(mydimx+2*mynframe))
         enddo
        enddo
        deallocate(datatmp8)
      endif
   else
     if(mygdatatype=='bin8') then
       data=datatmp8
       deallocate(datatmp8)
     endif
   endif
!---
   if(present(iret)) iret=0
   return
  end subroutine nemsio_readrec4w34
!------------------------------------------------------------------------------
!------------------------------------------------------------------------------
  subroutine nemsio_readrec8w34(gfile,jrec,data,nframe,iret)
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - - - - - -
! abstract: read nemsio data by record number into a 2D 32 bits array,
!           using w3_4 library to compile
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - - - - - -
    implicit none
    type(nemsio_gfile),intent(inout)              :: gfile
    integer(nemsio_intkind),intent(in)            :: jrec
    real(nemsio_dblekind),intent(inout)           :: data(:)
    integer(nemsio_intkind),optional,intent(out)  :: iret
    integer(nemsio_intkind),optional,intent(in)   :: nframe
    real(nemsio_realkind),allocatable             :: datatmp4(:)
    real(nemsio_dblekind),allocatable             :: datatmp(:)
    integer :: i,j,ios
!------------------------------------------------------------
! read 4 byte rec
!------------------------------------------------------------
   if(present(iret)) iret=-34
!---
   call nemsio_getgfile(gfile,iret)
!---
   if ( mygdatatype .eq. 'bin4') then
     allocate(datatmp4(myfieldsize) )
     call nemsio_readrecbin4d4(gfile,jrec,datatmp4,ios)
   else if ( mygdatatype .eq. 'bin8') then
     if(.not.present(nframe) ) then
      call nemsio_readrecbin8d8(gfile,jrec,data,ios)
     else
      allocate(datatmp(myfieldsize) )
      call nemsio_readrecbin8d8(gfile,jrec,datatmp,ios)
     endif
   else
     allocate(datatmp4(myfieldsize) )
     call nemsio_readrecgrb4w34(gfile,jrec,datatmp4,ios)
   endif
   if ( ios .ne.0 ) then
     if(present(iret)) then
       iret=ios
       return
     else
       call nemsio_stop
     endif
   endif
!---
   if ( present(nframe) ) then
     if(mygdatatype .eq. 'bin4'.or.mygdatatype .eq. 'grib' ) then
       do j=1,mydimy+2*mynframe-2*nframe
       do i=1,mydimx+2*mynframe -2*nframe
        data(i+(j-1)*(mydimx+2*mynframe-2*nframe))=datatmp4(i+nframe        &
          +(j-1+nframe)*(mydimx+2*mynframe))
       enddo
       enddo
       deallocate(datatmp4)
     else if(mygdatatype .eq. 'bin8') then
       do j=1,mydimy+2*mynframe-2*nframe
       do i=1,mydimx+2*mynframe -2*nframe
        data(i+(j-1)*(mydimx+2*mynframe-2*nframe))=datatmp(i+nframe       &
          +(j-1+nframe)*(mydimx+2*mynframe))
       enddo
       enddo
       deallocate(datatmp)
     endif
   else
     if(mygdatatype .eq. 'bin4'.or.mygdatatype .eq. 'grib' ) then
       data=datatmp4
       deallocate(datatmp4)
     endif
   endif
!---
   if(present(iret)) iret=0
   return
  end subroutine nemsio_readrec8w34
!------------------------------------------------------------------------------
  subroutine nemsio_readrecv4w34(gfile,name,levtyp,lev,data,nframe,iret)
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - - - - - -
! abstract: read nemsio data by record number into a 2D 32 bits array
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - - - - - -
    implicit none
    type(nemsio_gfile),intent(inout)              :: gfile
    character(*),intent(in)                       :: name
    character(*),intent(in),optional              :: levtyp
    integer(nemsio_intkind),optional,intent(in)   :: lev
    real(nemsio_realkind),intent(inout)           :: data(:)
    integer(nemsio_intkind),optional,intent(out)  :: iret
    integer(nemsio_intkind),optional,intent(in)   :: nframe
    real(nemsio_realkind),allocatable             :: datatmp(:)
    real(nemsio_dblekind),allocatable             :: datatmp8(:)
    integer :: i,j,ios
!------------------------------------------------------------
! read 4 byte rec
!------------------------------------------------------------
   if(present(iret)) iret=-35
!---
   call nemsio_getgfile(gfile,iret)
!---
   if ( mygdatatype .eq. 'bin4') then
     if(.not.present(nframe)) then
       call nemsio_readrecvbin4d4(gfile,name,levtyp,lev,data,ios)
     else
       allocate(datatmp(myfieldsize) )
       call nemsio_readrecvbin4d4(gfile,name,levtyp,lev,datatmp,ios)
     endif
   else if ( mygdatatype .eq. 'bin8') then
     allocate(datatmp8(myfieldsize) )
     call nemsio_readrecvbin8d8(gfile,name,levtyp,lev,datatmp8,ios)
   else
     if(.not.present(nframe)) then
       call nemsio_readrecvgrb4w34(gfile,name,levtyp,lev,data,ios)
     else
       allocate(datatmp(myfieldsize) )
       call nemsio_readrecvgrb4w34(gfile,name,levtyp,lev,datatmp,ios)
     endif
   endif
   if ( ios .ne.0 ) then
     if(present(iret)) then
       iret=ios
       return
     else
       call nemsio_stop
     endif
   endif
!---
   if ( present(nframe) ) then
      if(mygdatatype=='bin4'.or.mygdatatype=='grib') then
        do j=1,mydimy+2*mynframe-2*nframe
         do i=1,mydimx+2*mynframe -2*nframe
           data(i+(j-1)*(mydimx+2*mynframe-2*nframe))=datatmp(i+nframe        &
            +(j-1+nframe)*(mydimx+2*mynframe))
         enddo
        enddo
        deallocate(datatmp)
      elseif(mygdatatype=='grib8') then
        do j=1,mydimy+2*mynframe-2*nframe
         do i=1,mydimx+2*mynframe -2*nframe
           data(i+(j-1)*(mydimx+2*mynframe-2*nframe))=datatmp8(i+nframe        &
            +(j-1+nframe)*(mydimx+2*mynframe))
         enddo
        enddo
        deallocate(datatmp8)
      endif
   else
      if(mygdatatype=='grib8') then
        data=datatmp8
        deallocate(datatmp8)
      endif
   endif
!---
   if(present(iret)) iret=0
   return
  end subroutine nemsio_readrecv4w34
!
!------------------------------------------------------------------------------
  subroutine nemsio_readrecv8w34(gfile,name,levtyp,lev,data,nframe,iret)
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - - - - - -
! abstract: read nemsio data by record number into a 2D 32 bits array
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - - - - - -
    implicit none
    type(nemsio_gfile),intent(inout)              :: gfile
    character(*),intent(in)                       :: name
    character(*),intent(in),optional              :: levtyp
    integer(nemsio_intkind),optional,intent(in)   :: lev
    real(nemsio_dblekind),intent(inout)           :: data(:)
    integer(nemsio_intkind),optional,intent(out)  :: iret
    integer(nemsio_intkind),optional,intent(in)   :: nframe
    real(nemsio_dblekind),allocatable             :: datatmp(:)
    real(nemsio_realkind),allocatable             :: datatmp4(:)
    integer :: i,j,ios
!------------------------------------------------------------
! read 8 byte rec
!------------------------------------------------------------
   if(present(iret)) iret=-35
!---
   call nemsio_getgfile(gfile,iret)
!---
   if ( mygdatatype .eq. 'bin4') then
     allocate(datatmp4(myfieldsize) )
     call nemsio_readrecvbin4d4(gfile,name,levtyp,lev,datatmp4,ios)
   else if ( mygdatatype .eq. 'bin8') then
     if(.not.present(nframe)) then
       call nemsio_readrecvbin8d8(gfile,name,levtyp,lev,data,ios)
     else
       allocate(datatmp(myfieldsize) )
       call nemsio_readrecvbin8d8(gfile,name,levtyp,lev,datatmp,ios)
     endif
   else
     allocate(datatmp4(myfieldsize) )
     call nemsio_readrecvgrb4w34(gfile,name,levtyp,lev,datatmp4,ios)
   endif
   if ( ios .ne.0 ) then
     if(present(iret)) then
       iret=ios
       return
     else
       call nemsio_stop
     endif
   endif
!---
   if ( present(nframe) ) then
      if(mygdatatype .eq. 'bin4'.or.mygdatatype .eq. 'grib') then
       do j=1,mydimy+2*mynframe-2*nframe
       do i=1,mydimx+2*mynframe -2*nframe
        data(i+(j-1)*(mydimx+2*mynframe-2*nframe))=datatmp4(i+nframe        &
          +(j-1+nframe)*(mydimx+2*mynframe))
       enddo
       enddo
       deallocate(datatmp4)
      elseif(mygdatatype .eq. 'bin8') then
       do j=1,mydimy+2*mynframe-2*nframe
       do i=1,mydimx+2*mynframe -2*nframe
        data(i+(j-1)*(mydimx+2*mynframe-2*nframe))=datatmp(i+nframe        &
          +(j-1+nframe)*(mydimx+2*mynframe))
       enddo
       enddo
       deallocate(datatmp)
     endif
   else
     if(mygdatatype .eq. 'bin4'.or.mygdatatype .eq. 'grib') then
       data=datatmp4
       deallocate(datatmp4)
     endif
   endif
!
   if(present(iret)) iret=0
   return
  end subroutine nemsio_readrecv8w34
!------------------------------------------------------------------------------
!
!*****************   read bin data set :  ********************************
!
!------------------------------------------------------------------------------
  subroutine nemsio_readrecbin4d4(gfile,jrec,data,iret)
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - - - - - -
! abstract: read nemsio data (bin) by record number into a 2D 32 bits array
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - - - - - -
    implicit none
    type(nemsio_gfile),intent(in)                 :: gfile
    integer(nemsio_intkind),intent(in)            :: jrec
    real(nemsio_realkind),intent(inout)           :: data(:)
    integer(nemsio_intkind),optional,intent(out)  :: iret
    integer(nemsio_intkind8) :: iskip,iread,nread

    if(present(iret)) iret=-41
    iskip=mytlmeta+int(jrec-1,8)*int(kind(data)*myfieldsize+8,8)
    iread=int(nemsio_realkind,8)*int(size(data),8)
    call bafrreadl(myflunit,iskip,iread,nread,data)
    if(nread.lt.iread) return
    if(do_byteswap) call byteswap(data,nemsio_realkind,size(data))
    if(present(iret)) iret=0

    return
  end subroutine nemsio_readrecbin4d4
!------------------------------------------------------------------------------
  subroutine nemsio_readrecvbin4d4(gfile,name,levtyp,lev,data,iret)
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - - - - - -
! abstract: read nemsio data (bin) by record number into a 2D 32 bits array
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - - - - - -
    implicit none
    type(nemsio_gfile),intent(in)                 :: gfile
    character(*),intent(in)                      :: name
    character(*),intent(in),optional             :: levtyp
    integer(nemsio_intkind),optional,intent(in)  :: lev
    real(nemsio_realkind),intent(out)             :: data(:)
    integer(nemsio_intkind),optional,intent(out)  :: iret
    integer(nemsio_intkind8) :: iskip,iread,nread
    integer :: jrec, ierr

    if(present(iret)) iret=-42
    call nemsio_searchrecv(gfile,jrec,name,levtyp,lev,ierr)
    if ( ierr .ne. 0)  return
    iskip=mytlmeta+int(jrec-1,8)*int(nemsio_realkind*myfieldsize+8,8)
    iread=int(kind(data),8)*int(size(data),8)
    call bafrreadl(myflunit,iskip,iread,nread,data)
    if(nread.lt.iread) return
    if(do_byteswap) call byteswap(data,nemsio_realkind,size(data))
    if(present(iret)) iret=0

    return
  end subroutine nemsio_readrecvbin4d4
!------------------------------------------------------------------------------
  subroutine nemsio_readrecbin8d8(gfile,jrec,data,iret)
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - - - - - -
! abstract: read nemsio data (bin) by record number into a 2D 32 bits array
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - - - - - -
    implicit none
    type(nemsio_gfile),intent(in)                 :: gfile
    integer(nemsio_intkind),intent(in)            :: jrec
    real(nemsio_dblekind),intent(out)             :: data(:)
    integer(nemsio_intkind),optional,intent(out)  :: iret
    integer(nemsio_intkind8) :: iskip,iread,nread

    if(present(iret)) iret=-42
    iskip=mytlmeta+int(jrec-1,8)*int(nemsio_dblekind*myfieldsize+8,8)
    iread=int(nemsio_dblekind,8)*int(size(data),8)
    call bafrreadl(myflunit,iskip,iread,nread,data)
    if(nread.lt.iread) return
    if(do_byteswap) call byteswap(data,nemsio_dblekind,size(data))
    if(present(iret)) iret=0

    return
  end subroutine nemsio_readrecbin8d8
!------------------------------------------------------------------------------
  subroutine nemsio_readrecvbin8d8(gfile,name,levtyp,lev,data,iret)
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - - - - - -
! abstract: read nemsio data (bin) by record number into a 2D 32 bits array
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - - - - - -
    implicit none
    type(nemsio_gfile),intent(in)                 :: gfile
    character(*),intent(in)                      :: name
    character(*),intent(in),optional             :: levtyp
    integer(nemsio_intkind),optional,intent(in)  :: lev
    real(nemsio_dblekind),intent(out)             :: data(:)
    integer(nemsio_intkind),optional,intent(out)  :: iret
    integer(nemsio_intkind8) :: iskip,iread,nread
    integer :: jrec, ierr

    if(present(iret)) iret=-44
    call nemsio_searchrecv(gfile,jrec,name,levtyp,lev,ierr)
    if ( ierr .ne. 0) return
    iskip=mytlmeta+int(jrec-1,8)*int(nemsio_dblekind*myfieldsize+8,8)
    iread=int(nemsio_dblekind,8)*int(size(data),8)
    call bafrreadl(myflunit,iskip,iread,nread,data)
    if(nread.lt.iread) return
    if(do_byteswap) call byteswap(data,nemsio_dblekind,size(data))
    if(present(iret)) iret=0

    return
  end subroutine nemsio_readrecvbin8d8
!------------------------------------------------------------------------------
!
!*****************   read w34 data set :  *************************************
!
!------------------------------------------------------------------------------
  subroutine  nemsio_readrecgrb4w34(gfile,jrec,data,iret)
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - - - - - -
! abstract: read nemsio data by record number into a 2D 32 bits array,
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - - - - - -
    implicit none
    type(nemsio_gfile),intent(inout)              :: gfile
    integer(nemsio_intkind),intent(in)            :: jrec
    real(nemsio_realkind),intent(out)             :: data(:)
    integer(nemsio_intkind),optional,intent(out)  :: iret
    type(nemsio_grbmeta)         :: grbmeta
    integer(nemsio_intkind)      :: luidx
    integer(nemsio_intkind)      :: kf,k,kpds(200),kgds(200)
    logical*1,allocatable       :: lbms(:)
    integer(nemsio_intkind)      :: N=nemsio_kpds_intfill
    integer(nemsio_intkind)      :: ios,i,w34
!
!------------------------------------------------------------
! set up grib meta
!------------------------------------------------------------
    luidx=0
    if ( present(iret)) iret=-45
    w34=1
    call nemsio_setrqst(gfile,grbmeta,ios,jrec=jrec,w34=w34)
    if (ios.ne.0) then
       if ( present(iret))  then
         iret=ios
         return
       else
         call nemsio_stop
       endif
    endif
    allocate(lbms(grbmeta%jf))
    N=0
!------------------------------------------------------------
! get data from getgb
!------------------------------------------------------------
    call getgbm(myflunit,luidx,grbmeta%jf,N,grbmeta%jpds,grbmeta%jgds,&
      mymbuf,mycbuf,mynlen,mynnum,mymnum, &
      kf,k,kpds,kgds,lbms,data,ios)
    deallocate(lbms,grbmeta%lbms)
    if(ios.ne.0) then
       if ( present(iret))  then
         print *,'getgb_ios=',ios
         return
       else
         call nemsio_stop
       endif
    endif
    if (present(iret)) iret=0
  end subroutine nemsio_readrecgrb4w34
!------------------------------------------------------------------------------
  subroutine nemsio_readrecvgrb4w34(gfile,vname,vlevtyp,vlev,data,iret)
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - - - - - -
! abstract: read nemsio data by field name into 32 bits array,
!           using w3_4 library to compile
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - - - - - -
    implicit none
    type(nemsio_gfile),intent(inout)              :: gfile
    character*(*),intent(in)                      :: vname,vlevtyp
    integer(nemsio_intkind),intent(in)            :: vlev
    real(nemsio_realkind),intent(out)             :: data(:)
    integer(nemsio_intkind),optional,intent(out)  :: iret
    type(nemsio_grbmeta)         :: grbmeta
    integer(nemsio_intkind)      :: luidx
    integer(nemsio_intkind)      :: kf,k,kpds(200),kgds(200)
    logical*1,allocatable       :: lbms(:)
    integer(nemsio_intkind)      :: N=nemsio_kpds_intfill
    integer(nemsio_intkind)      :: ios,i,w34
!
!------------------------------------------------------------
! set up grib meta
!------------------------------------------------------------
    luidx=0
    if ( present(iret)) iret=-45
    w34=1
    call nemsio_setrqst(gfile,grbmeta,ios,vname=vname, &
      vlevtyp=vlevtyp, vlev=vlev ,w34=w34)
    if (ios.ne.0) then
       if ( present(iret))  then
         iret=ios
         return
       else
         call nemsio_stop
       endif
    endif
!------------------------------------------------------------
! get data from getgb _w34
!------------------------------------------------------------
    allocate(lbms(grbmeta%jf))
    N=0
    call getgbm(myflunit,luidx,grbmeta%jf,N,grbmeta%jpds,grbmeta%jgds,&
      mymbuf,mycbuf,mynlen,mynnum,mymnum, &
      kf,k,kpds,kgds,lbms,data,ios)
    deallocate(lbms,grbmeta%lbms)
    if(ios.ne.0) then
       if ( present(iret))  then
          print *,'getgb_ios=',ios
         return
       else
         call nemsio_stop
       endif
    endif
    if ( present(iret)) iret=0
  end subroutine nemsio_readrecvgrb4w34
!------------------------------------------------------------------------------
!
!*****************   read grb data set w3d:  *************************************
!
!------------------------------------------------------------------------------
  subroutine nemsio_readrecgrb8(gfile,jrec,data,iret)
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - - - - - -
! abstract: read nemsio data by record number into a 2D 64 bits array,
!           using w3_d library to compile
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - - - - - -
    implicit none
    type(nemsio_gfile),intent(inout)              :: gfile
    integer(nemsio_intkind),intent(in)            :: jrec
    real(nemsio_dblekind),intent(out)             :: data(:)
    integer(nemsio_intkind),optional,intent(out)  :: iret
    type(nemsio_grbmeta)         :: grbmeta
    integer(nemsio_intkind)      :: luidx
    integer(nemsio_intkind)      :: kf,k,kpds(200),kgds(200)
    logical*1,allocatable       :: lbms(:)
    integer(nemsio_intkind)      :: N=nemsio_kpds_intfill
    integer(nemsio_intkind)      :: ios,i
!
!------------------------------------------------------------
! set up grib meta
!------------------------------------------------------------
    luidx=0
    if ( present(iret)) iret=-46
    call nemsio_setrqst(gfile,grbmeta,ios,jrec=jrec)
    if (ios.ne.0) then
       if ( present(iret))  then
         iret=ios
         return
       else
         call nemsio_stop
       endif
    endif
!------------------------------------------------------------
! get data from getgb _w3d
!------------------------------------------------------------
    allocate(lbms(grbmeta%jf))
    N=0
    call getgbm(myflunit,luidx,grbmeta%jf,N,grbmeta%jpds,grbmeta%jgds,&
      mymbuf,mycbuf,mynlen,mynnum,mymnum, &
      kf,k,kpds,kgds,lbms,data,ios)
    deallocate(lbms,grbmeta%lbms)
    if(ios.ne.0) then
       if ( present(iret))  then
         print *,'getgb_ios=',ios
         return
       else
         call nemsio_stop
       endif
    endif
    if (present(iret)) iret=0
  end subroutine nemsio_readrecgrb8
!------------------------------------------------------------------------------
  subroutine nemsio_readrecvgrb8(gfile,vname,vlevtyp,vlev,data,iret)
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - - - - - -
! abstract: read nemsio data by field name into a 2D 64bits array,
!           using w3_d library to compile
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - - - - - -
    implicit none
    type(nemsio_gfile),intent(inout)              :: gfile
    character*(*),intent(in)                     :: vname,vlevtyp
    integer(nemsio_intkind),intent(in)            :: vlev
    real(nemsio_dblekind),intent(out)             :: data(:)
    integer(nemsio_intkind),optional,intent(out)  :: iret
    type(nemsio_grbmeta)         :: grbmeta
    integer(nemsio_intkind)      :: luidx
    integer(nemsio_intkind)      :: kf,k,kpds(200),kgds(200)
    logical*1,allocatable       :: lbms(:)
    integer(nemsio_intkind)      :: N=nemsio_kpds_intfill
    integer(nemsio_intkind)      :: ios,i
!
!------------------------------------------------------------
! set up grib meta
!------------------------------------------------------------
    luidx=0
    if ( present(iret)) iret=-47
    call nemsio_setrqst(gfile,grbmeta,ios,vname=vname, &
      vlevtyp=vlevtyp, vlev=vlev )
    if (ios.ne.0) then
       if ( present(iret))  then
         iret=ios
         return
       else
         call nemsio_stop
       endif
    endif
!------------------------------------------------------------
! get data from getgb _w3d
!------------------------------------------------------------
    allocate(lbms(grbmeta%jf))
    N=0
    call getgbm(myflunit,luidx,grbmeta%jf,N,grbmeta%jpds,grbmeta%jgds,&
      mymbuf,mycbuf,mynlen,mynnum,mymnum, &
      kf,k,kpds,kgds,lbms,data,ios)
    deallocate(lbms,grbmeta%lbms)
    if(ios.ne.0) then
       if ( present(iret))  then
          print *,'getgb_ios=',ios
         return
       else
         call nemsio_stop
       endif
    endif
    if ( present(iret)) iret=0
  end subroutine nemsio_readrecvgrb8
!------------------------------------------------------------------------------
end module nemsio_read

max_stars_repo_name: GMAO-SI-Team/jedi-stack

max_issues_repo_name: GMAO-SI-Team/jedi-stack

max_forks_repo_name: GMAO-SI-Team/jedi-stack


text:
      MODULE kopp2_I
      INTERFACE
!...Generated by Pacific-Sierra Research 77to90  4.3E  14:44:54  12/27/06
      SUBROUTINE kopp2 (POS, RAD3, J, JPRIM, PAR, ANTKO)
      integer, INTENT(IN) :: POS
      character (LEN = 200), INTENT(OUT) :: RAD3
      integer, DIMENSION(20), INTENT(IN) :: J
      integer, DIMENSION(20), INTENT(IN) :: JPRIM
      integer, INTENT(IN) :: PAR
      integer, DIMENSION(20), INTENT(IN) :: ANTKO
      END SUBROUTINE
      END INTERFACE
      END MODULE

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: sylas/grasp-continuum

max_forks_repo_name: sylas/grasp-continuum


text:
Chicken wings are often spicy foods spicy, and widely loved bits of chicken fried or BBQd up at sordid dumps (and the occasional quality Restaurants restaurant!) around the nation, you can also make them yourself. They come in many varieties, the most common being Buffalo style (commonly dubbed hot wings, although there are plenty of other varieties of spicy wings that fit the bill) and BBQ.  Theyre also naturally a Super Bowl staple.

Wings come in the bonein form and boneless (which are not technically wings).  Historically, chicken wings (as in the actual part of the chicken) were dirt cheap compared to the meat of the chicken, running less than 25% the cost of breast meat in the mid 1990s.  By 2010, the wholesale cost of wings had jumped to 140% the price of boneless, skinless breast meat, not only causing a price inversion, but making it much more cost effective to push boneless wings cut out of breast meat.

The best place for getting wings that are high quality, local, and cruelty/cage free is Nugget! They offer a rocky Jr wings that are billed as fryer wings but are perfect for grilling. Its also entirely the first cut and not the 2nd (1st cut is basically shoulder to elbow, 2nd is elbow to hand someone correct me if Im wrong) they make a solid base to start from when you are going to BBQ 	 Daubert

Should you feel bad about eating chicken wings? Maybe. Check out the Ethical Eating page and decide for yourself.

 Davis (and Nearby) Restaurants Serving Wings
Add your favorite wings destinations!
 Applebees  A dozen boneless or 10 bonein wings for $8.99.  They have Classic Buffalo, Hot Buffalo, Honey BBQ, Southern BBQ, and Sweet & Spicy sauces to choose from.
   7/10  Cheap and surprisingly delicious.  The Hot Buffalo isnt particularly hotpretty average buffalo wings, meaning the classic would probably be a bit on the mild side.  The sauce is delicious, though.  More peppery than most.  Wings came out just slightly on the greasylooking side, but it didnt detract from the taste or the texture.
 Cornish Kitchen and Catering
 Dominos Pizza  $5.99/10pcs.  Available in Buffalo, hot Buffalo, and BBQ, both bonein and boneless.  Hot Buffalo is your typical hot wings.  Dominos delivers, of course, with a delivery fee of $1.99.
   4/10  These wings are greasy and have very little sauce.  The sauce itself is decent, and the price is pretty low, but unless you dont like wings enough to care, its not worth it.  Good amount of meat. Users/TomGarberson
  The Graduate  The hot wings are nice and spicy, and quite tasty.  Its a fairly standard Buffalo style sauce.  They run about $9/dozen.  They also have BBQ wings... anyone whos tried them want to comment?
   67.5/10  The sauce is a little thin but tastes great; texture isnt quite as good as Sudwerk, but still very respectable.  Decent amount of meat on the bones; more than youll get from the storebought party wings, but not as much as some of the primo places like Sudwerk.
    June 2010 update: the last two times I got these, they had a lot less hot sauce than usual and werent nearly as spicy.  Last week, I didnt even get the usual cup of blue cheese dressing with my order.  Not sure if theres a new chef, if theyre just off their game, or if theyre getting stingier with the sauces.  Id drop the rating to 6.5 or 7 if the last couple of times become the standard. Users/TomGarberson
  Lamppost Pizza  They sell both Buffalo Wings and Barbeque Wings, that come with ranch sauce.    6 wings is a half order, 12 is a full order ($11), and they also do 20 wings for sharing.
   68/10  Variable quality, ranging from good to great hot wings.  I got a dozen to go during the world cup.  Huge amount of meat on the wingsa dozen easily filled me up.  They dont come heavily sauced, but they include plenty of sauce on the side (buffalo and ranch were both included, and I didnt even have to ask).  Price is reasonable, slightly below average in town.  They were very slightly greasy looking, but that may have been because they were cooked and sat in foil for 15 or 20 minutes before I partook.  I want to try them there next time, fresh out of the fryer.  They might well bump up to 9/10.  Update: on subsequent visits they havent been as good.  Consistently a little on the greasy side.  Still very tasty sauce with decent meat on the bones.  Even fresh out of the fryer, they arent quite as crispy as Id like. Still reasonably good, just not all Id hoped for. Users/TomGarberson
  9/10  Have gotten wings twice from here recently. As noted above, there is much sauce on the wings themselves, but there is a lot of sauce on the side. The skin is slightly crispy. One of the times the wings were sort of small but they through in a few extras, so no complaints here.  mikey crews
    Little update: the next time I tried them, they really werent that good.  Not much sauce and still greasy, even fresh from the oven.  Second try earns them a 6.  Ill have to give them a few more chances to see how things average out. tg
  Little Caesars  $5 for 8 wings.
     4 out of 10 stars. 1 star for the price, 2 stars for how fast I got them after I ordered them (15 seconds), and 1 star for the taste (didnt taste too well). Users/NikhilDahal
     4/10  I agree with Nikhil on this one.  Fast, easy, cheap, with a flavor and texture texture to match.  Users/TomGarberson
  Original Steves  $7.99 for a dozen.
  Plainfield Station  $5.35 for a small order of hot wings as of 2007.
  Pyramid Alehouse  Not in Davis, but its a great destination in Sacramento, Berkeley, or Walnut Creek if you want some interesting hot wings.
   10/10  truly the best hot wings Ive ever had.  Theyre atypical: they have a nice, crispy crust after being fried, and arent sauced.  Instead, theyre served with a habanero hot sauce on the side.  Massive amounts of meat on the bone, and if I recall correctly, theyre all drumettes.  The crispy coating is some sort of breading, and it gives a great texture.  It also means that the meat stays incredibly moist.  Be warned, though: the crust traps the heat inside, as well as the moisture, so they stay extremely hot.  Its easy to burn your tongue before you even get to the deliciously spicy hot sauce. Users/TomGarberson
   As of May 2011, Pyramids wings are available tossed or sauce on the side.  Tossed gets you slightly less crispy but ohsodelicious wings as the breading soaks up a ton of that yummy habanero sauce.  They also apparently have a house Louisianastyle sauce as an alternative to the habanero sauce.  Havent tried it out yet. tg
  Sudwerk  Another Buffalo style sauce.  You can get the wings by the halfdozen or dozen as an appetizer ($6.25 for 6, $11.25 for 12), or get the hot wings entree which lands you a plate of a dozen wings, a pile of their thick fries, and a big heap of onion rings, for $12.75.  You can also substitute in a fruit cup or cole slaw or whatever, if you happen not to hate your arteries.
   9/10  quite possibly the best wings in Davis.  Great sauce, nice and thick, good breading, tons of meat on the bones! Users/TomGarberson
  Thai Canteen  $5 for 1/2 dozen at happy hour.
   4/10  Not much flavor here.  The wings have little or no seasoning and no sauce when they come out.  Theres a mild, sweet and sourtype dipping sauce on the side.  The wings are well cooked and have a little bit of crispiness, which is nice, but they just arent very interesting.
  Tommy Js, the kitchen in Froggys  $5.75 1/2 dozen, $10.25 dozen for hot wings.
   7/10  They use Franks Red Hot with a decent ratio of Franks to butter.  Not super spicy, but theres a good bit of zing to them.  Average amount of meat on the bones (more than you get if you buy wings at the supermarket, but less than, e.g., Lamppost, Pyramid, or Sudwerk).  Good price and quick service certainly dont hurt.  You might try a half dozen with one of their tasty burgers.
  Uncle Vitos Slice of N.Y.  They have wings in Buffalo (traditional hot wing), BBQ, and Thai Chili sauces.  A Small order is 1 pound of wings, 810 pcs., for (I think) $7.95; a Large is 2 pounds, 1620 pcs., for $13.95.
   6/10  Very crispy skin, which I like a lot.  The buffalo sauce isnt that spicy, though, nor is it that flavorful.  A large order still left me a bit hungry; theres really not all that much meat on the bones, despite getting 1620 pieces (Id guess no more than 16).  While the buffalo sauce isnt particularly inspiring, the great crispiness would probably lend itself well to any sauce, so if they have a good BBQ or Thai sauce, Ill bet itd be delicious.  Not bad, but if youre looking for good wings, theres better available in town. Users/TomGarberson
  Village Pizza & Grill  $7.79 for 1/2 dozen; $12.79 for 1 dozen.
   7/10  The wings were pricey, but quite good.  Smallish wings, but very tasty sauce.  Its fairly spicy and a bit sweet.  Great color  a few shades darker than the usual Buffalo sauce.  The gorgonzola sauce that comes with the wings is not good.  I dont know whats in it, but it looked and tasted like a subpar buttermilk ranch dressing or something of that sort.  Since I dont dip my wings it wasnt really an issue, but if you are a dipper, look elsewhere.  Presentation of the wings was really nice  the dark sauce was very pretty, and they were laid out on two very fresh leaves of butter lettuce.  If they werent so expensive and if the wings themselves were a little larger, theyd be 8/10.  Users/TomGarberson
  Wingstop  coming soon
  Woodstocks Pizza  Half(6) for $5.29, full(12) for $7.99, and platter(18) for $11.99.
  3rd & U Cafe  Happy hour price half for $3.25
  8/10  Got em at happy hour and super cheap! Lots of sauce, decent spice level (I love spicy foods, though, so it might come off as quite a bit spicier to others). Decent sized, not very crunchy but again, its because they are welldrenched in the sauce. Came with a side of what looks like ranch with bits of blue cheese in it? The half order, with the happy hour fries, comes out to a pretty good meal at under $5. So, with taste, size and cost playing a huge factor, an 8. May be different had they not been happy hour priced! Users/HannahToru

 Sauces
Hot
Hot wings are typically served with some celery and either ranch or bleu cheese dressing.  With good wings, though, none of that is remotely necessary.  Live the heat, love the heat!

  Davis has an extraordinary selection of hot sauce hot sauces.  If the selection on the hot sauce hot sauce entry spins your head, just try Franks Red Hot for a start.  It is one of the more common sauces with a classic Buffalo taste.  Mix with melted butter for a basic Buffalo sauce.

BBQ
 Sweet Baby Rays is fantastic, but there are countless options in your local grocery store.  You can also http://allrecipes.com/Search/Recipes.aspx?WithTermbarbecue%20sauce&SortByRating&DirectionDescending make your own!
 Stubbs can give sweet baby rays a run for its money and is also pretty good
 WingTime Buffalo sauce was created in Davis, its good and has roughly 5 levels of spice
 Everett and Jones BBQ sauce is a bay area tradition sine the 60s The yellow lid medium spice is great
 Roadhouse thick & chunky can be order online (used to be able to get it at savemart) get the spicy

Other
Have a look around for recipes.  http://www.allrecipes.com Allrecipes is a great resource with a huge number of recipes and reviews.  There are some great sauces out there, and they dont have to be associated with wings.  Spicy Thai and Caribbean sauces go great.  You can even do teriyaki and other mild, sweet sauces.  Wing purists will likely mock you for it, but thats OK.  Have a little selfesteem and do it your way.  Who cares what anyone else things?  Theyre your wings!

 Other Ways to Get Your Fix
 Frozen  Supermarkets & Grocery Stores Supermarkets typically have a few different sorts of bagged wings in the frozen foods aisles.  Safeway has TGI Fridays wings as well as a couple of other options.  These arent nearly as good as fresh wings, and they arent that much cheaper.  Dont settle!  Costco has some very good frozen wings for a good price.
 Deli wings  Safeway sometimes (always?) has wings in the deli section.  Not sure on the price.

Do It Yourself!
Theres a reason wings are often found at fairly basic eateries: theyre easy to make.  Raw wings are available any place with a butcher (except the Meat Lab, which doesnt do chicken).  Theyre pretty cheap  typically $2/lb. or under.  You can get the whole wings and break them down yourself; you can have the butcher cut them up; or you can often find packs that are already cut into drumettes and wings.  If youre breaking them down yourself, you may want to discard the tips, as theres very little meat on them.  The wings you get in at the store tend to have a little less meat on them than the restaurant wings.  Presumably the stores simply sell a slightly lower grade of bird, as with beef.

__Frying__

 1. Flour: tossing the wings in seasoned flour gives you the crispy crust after cooking.  Use allpurpose flour, bread flour, whatever you prefer.  For seasonings, throw in a little of whatever you like.  Paprika, cayenne, chili powder, cumin, salt, pepper, garlic powder, onion powder, dried leaf herbs, and so on are all fair game.  Any good dry rub for grilling will go great, too.
 2. Refrigerate: its tempting to skip this step, since it adds a fair amount of time.  Dont.  Refrigerating after tossing the wings in the flour will give you the great texture you know and love at the restaurants.  If you dont do this, chances are youre going to wind up with a gummy mess on the outside of the wings.
 3. Cook: 1015 minutes in 375 degree oil will do the trick, or ~40 minutes in the oven at 375.  Naturally, always check for doneness.  When you pull them out, drain any excess grease and toss them into your sauce!
 4. Sauce: There are a lot of premixed sauces out there that have everything you need.  Likewise, you can use BBQ sauces asis.  If you get the bottles of Franks hot sauce and you dont want to go too spicy, mix it with some melted butter.  Check out the extensive entry about hot sauce here on the wiki for a wide range of flavor and heat options.

__Grilling__

 1. Marinate (Optional): Just about anything works.  If you want to go quick and easy, just toss them in a little Italian dressing and soy sauce and let them sit for a while.
 2. Grill: Medium heat, 20ish minutes, make sure you dont burn them.
 3. Sauce: After they cool for a few minutes, you can just toss them in a large ziplock bag with your sauce, seal, and shake to coat beautifully. Alternatively you can sauce them for the last 5 minutes they are on the grill, Its recommended to do 23 solid coats with no turning once sauced.

20100128 09:43:05 nbsp Any other wing fans out there interested in going around and sampling the wings at the various places offering them?  Itd be fun to put together some reviews and rankings.  And my wife has a much lower threshold for that delicious burning, so I dont think shes in. Users/TomGarberson

20110621 19:57:32 nbsp Ding How has the best, their not on the menu. I think they still might make them, it has been awhile since ive been there, their called Sauteed chicken wings. And their killer, let me know. Users/Theangrydiner

20111216 10:34:55 nbsp ExBuffalonian here...nothing is even close to back home. Ive tried a bunch in Davis and Sac. Looks like Ill have to invest in a deep fryer :) Users/sba

20121120 01:07:25 nbsp If youre in Sacramento, Wingstop is quite good, if a little greasy. Its a chain, so dont expect an immaculate dining experience. But if youre deciding to eat wings you probably dont care about that anyway. Users/CecilioPadilla

And here to the content discussion: http://daviswiki.org/Chicken_Wings/Talk

20130207 16:49:17 nbsp Could someone do a review of the wings at 3rd and U? :D Users/HannahToru
 Nevermind, did it myself xD ht

20130304 09:37:41 nbsp Growing up in Buffalo I came to think that wings were as ubiquitous as burgers, fries and pizza. In fact every pizza place in the Buffalo/Rochester area has wings on the menu as do the bars with a kitchen.  Leaving Buffalo made me appreciate how good those chicken wings were. There are places that serve pretty bad wings there but on the whole most are pretty good. Here in the Davis/Sacramento area there are fewer places to get decent wings but there is one place in Sacramento that is as close to the average Buffalo NY wing place as you’ll get anywhere called Sam’s Sub shop on Fulton Ave. They are just called wings here as they should be and come in 4 different heats and sizes from 8 to a bucket of 50.  This is the norm back home. Here in Davis I will usually make my own (I’d put my sauce recipe against the best anywhere) but if I go out I go to Sudwerk. Not great but will do if I need a quick fix. Users/pedalpower

20130309 12:06:11 nbsp PYRAMID ALEHOUSE HAS CLOSED MARCH 4, 2013 :(((((( Users/HannahToru

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      SUBROUTINE CMCKCD
C
C     THIS SUBROUTINE DETERMINES WHETHER MANUALLY SPECIFIED CONNECTION
C     ENTRIES ARE ALLOWABLE BASED ON THE PRESCRIBED GEOMETRIC TOLERANCE.
C
      INTEGER         SCSFIL,COMBO,SCORE,IST(7),SCCONN,CE(9),AAA(2),
     1                BUF2,OUTT
      DIMENSION       IPNUM(7),COORD(7,3),DIFF2(3)
      CHARACTER       UFM*23,UWM*25,UIM*29
      COMMON /XMSSG / UFM,UWM,UIM
      COMMON /CMB001/ SCR1,SCR2,SCBDAT,SCSFIL,SCCONN,SCMCON,SCTOC,
     1                GEOM4,CASECC,SCCSTM,SCR3
      COMMON /CMB002/ BUF1,BUF2,BUF3,BUF4,BUF5,SCORE,LCORE,INTP,OUTT
      COMMON /CMB003/ COMBO(7,5),CONSET,IAUTO,TOLER,NPSUB,CONECT,TRAN,
     1                MCON,RESTCT(7,7),ISORT,ORIGIN(7,3),IPRINT
      COMMON /ZZZZZZ/ Z(1)
      COMMON /BLANK / STEP,IDRY
      DATA    AAA   / 4HCMCK,4HCD   /
C
C     READ ALL BGSS INTO OPEN CORE
C
      IT    = 2
      IERR  = 0
      LLCO  = LCORE
      J     = 0
      IFILE = SCSFIL
      CALL OPEN (*200,SCSFIL,Z(BUF2),0)
      DO 30 I = 1,NPSUB
      NREC  = COMBO(I,5) + 1
      DO 10 JJ = 1,NREC
      CALL FWDREC (*210,SCSFIL)
   10 CONTINUE
      CALL READ (*210,*20,SCSFIL,Z(SCORE+J),LLCO,1,NNN)
      GO TO 220
   20 IST(I) = SCORE + J
      J = J + NNN
      LLCO  = LLCO - NNN
      CALL SKPFIL (SCSFIL,1)
   30 CONTINUE
      CALL CLOSE (SCSFIL,1)
C
C     READ CONNECTION ENTRIES AND LOAD INTO COORD ARRAY
C
      IFILE = SCCONN
      CALL OPEN (*200,SCCONN,Z(BUF2),0)
   40 CALL READ (*180,*50,SCCONN,CE,10,1,NNN)
C
C     LOAD COORD ARRAY
C     CE(3)... UP TO CE(9) ARE INTERNAL POINT NO.
C     IZ(IADD) IS THE COORD (CSTM) ID OF THE INTERNAL PTS.
C     Z(IADD+1,+2,+3) ARE THE COORD. ORIGINS
C
   50 NPT  = 0
      DO 80 I = 1,NPSUB
      IF (CE(I+2)) 80,80,60
   60 NPT  = NPT + 1
      IADD = 4*(CE(I+2)-1) + IST(I)
      IPNUM(NPT) = CE(I+2)
      DO 70 J = 1,3
      COORD(NPT,J) = Z(IADD+J)
   70 CONTINUE
   80 CONTINUE
C
C     COMPARE ALL PAIRS OF COORDINATES AGAINST TOLER.
C
      NPTM1 = NPT - 1
      DO 170 I = 1,NPTM1
      IT = IT - 1
      JJ = I  + 1
      DO 160 J = JJ,NPT
      DO 90 KK = 1,3
      DIFF2(KK) = (COORD(J,KK)-COORD(I,KK))**2
   90 CONTINUE
      SUM  = 0.0
      DO 100 KK = 1,3
      SUM  = SUM + DIFF2(KK)
  100 CONTINUE
      DIST = SQRT(SUM)
      IF (DIST .LE. TOLER) GO TO 160
      IF (IT .GT. 1) GO TO 120
      WRITE  (OUTT,110) UFM
  110 FORMAT (A23,' 6514, ERRORS HAVE BEEN FOUND IN MANUALLY SPECIFIED',
     1       ' CONNECTION ENTRIES. SUMMARY FOLLOWS')
      IERR = 1
      IDRY =-2
      IT   = 2
  120 IF (IT .GT. 2) GO TO 140
      WRITE  (OUTT,130) (CE(KDH),KDH=1,NNN)
  130 FORMAT ('0*** GEOMETRIC ERRORS HAVE BEEN FOUND IN THE FOLLOWING',
     1        ' CONNECTION ENTRY', /5X,9I10)
      IT   = 3
  140 WRITE (OUTT,150) IPNUM(I),(COORD(I,MM),MM=1,3),
     1                 IPNUM(J),(COORD(J,MM),MM=1,3)
  150 FORMAT ('0*** IP NUMBER',I10,13H COORDINATES  ,3E16.6,4H AND, /,
     1        '     IP NUMBER',I10,13H COORDINATES  ,3E16.6,
     2        ' ARE NOT WITHIN TOLER UNITS.')
  160 CONTINUE
  170 CONTINUE
      GO TO 40
C
  180 IF (IERR .EQ. 0) WRITE (OUTT,190) UIM
  190 FORMAT (A29,' 6516, ALL MANUAL CONNECTIONS SPECIFIED ARE ',
     1       'ALLOWABLE WITH RESPECT TO TOLERANCE')
      CALL CLOSE (SCCONN,1)
      GO TO 250
C
  200 IMSG = -1
      GO TO 230
  210 IMSG = -2
      GO TO 230
  220 IMSG = -8
  230 CALL MESAGE (IMSG,IFILE,AAA)
C
  250 RETURN
      END


text:
      program mtrinit
      external trinit
      common /ctrinit/ i

      print *, i

      end

      blockdata trinit
      common /ctrinit/ i
      data i /3/
      end


text:
      SUBROUTINE DTRC2C (B,NB,C1,YTAB, YPTAB)
c Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
c ALL RIGHTS RESERVED.
c Based on Government Sponsored Research NAS7-03001.
C>> 1994-11-11 DTRC2C Krogh   Declared all vars.
C>> 1994-10-20 DTRC2C Krogh  Changes to use M77CON
C>> 1987-12-09 DTRC2C Lawson  Initial code.
c--D replaces "?": ?TRC2C, ?C2BAS
c          Given C1() containing coefficients relative to a
c     cubic B-spline basis, this subr computes and stores into
c     YTAB() and YPTAB()  respectively, the value and first derivative
c     of the spline curve at each nodal point given in B().
C     INPUT..
C     (B(I),I=1,NB)     PARTITION POINTS,INCLUDING ENDPOINTS
C     NB                NO. OF POINTS IN B().
C                       NO. OF SEGMENTS IS NSEG = NB-1
c     (C1(I),I=1,NB+2)   Coefficients relative to a cubic B-spline basis
C
C     OUTPUT..
c     (YTAB(i),i=1,NB) and (YPTAB(i),i=1,NB)
c                    Values of spline function and its first
c                    derivative at the points given in B().
c                    YTAB(2*j-1) = f(B(j),   YPTAB(2*j) = fprime(B(j)).
c     ------------------------------------------------------------------
c     Based on code written by Lawson, Hanson, Lang, and Campbell, JPL,
c     1968-1974.
c     Modified July 1984 and July 1987 by C. L. Lawson, JPL, for
c     inclusion in the MATH 77 library.
c     ------------------------------------------------------------------
c     Subprogram referenced:  DC2BAS
c     ------------------------------------------------------------------
      integer NB, MB, NSEG, IB, ISEG, ISEGM1, IDERIV, I
      double precision B(NB), C1(NB+2), YTAB(NB), YPTAB(NB), P(4)
      double precision SUM, ZERO
      logical NEWSEG
      parameter( ZERO = 0.0D0)
c     ------------------------------------------------------------------
      MB=NB
      NSEG=MB-1
      do 30  IB=1,MB
         ISEG=min(NSEG,IB)
         ISEGM1=ISEG-1
         NEWSEG = .TRUE.
         do 20  IDERIV=0,1
            call DC2BAS(B(IB), IDERIV, ISEG, NEWSEG, B, NB, P)
            SUM=ZERO
            do 10  I=1,4
               SUM = SUM + C1(ISEGM1+I)*P(I)
   10       continue
            if(IDERIV .eq. 0) then
               YTAB(IB) = SUM
            else
               YPTAB(IB) = SUM
            endif
   20    continue
   30 continue
      return
      end


text:
Alive and Kicking was a local media periodical focused on the music scene of the area. It was published out of Sacramento by local music promoter, Jerry Perry, and ran from August 1991 to September 2008. Alive and Kicking Presents still holds shows every Thursday, Friday, and Saturday at Old Ironsides.

Place to find Alive and Kicking in Davis:
  Delta of Venus


! filename = sum.f90
! Purpose: To be used as a shared library

! Function to sum two real valued numbers
FUNCTION sum(x,y)
    double precision :: sum
    double precision, intent(in) :: x, y

    sum = x+y
END FUNCTION sum

! Subroutine to add two real valued numbers
SUBROUTINE sum_sub(x,y,output)
    double precision, intent(in) :: x, y
    double precision, intent(out) :: output

    output = x+y
END SUBROUTINE sum_sub

! function to calculate average of array of doubles
FUNCTION average(length, values)
    integer, intent(in) :: length
    double precision, intent(in) :: values(length)
    double precision :: average

    average = 0.0
    DO i = 1, length
        average = average + values(i)
    ENDDO
    average = average/length
END FUNCTION average