c**************************************************************
c     subroutine output2d()
c     Purpose
c     =======
c     Subroutine to output field variables at nodes and at
c     Gauss points, and for error norm computations for the 
c     Natural Element Method (NEM)
c
c**************************************************************
c
      subroutine output2d()
      implicit none
c
c     include files for the output phase (NEM program)
c
#include "include/params2d.h"
#include "include/node.h"
#include "include/neighbors.h"
#include "include/triangle.h"
#include "include/quadrature.h"
#include "include/gausspt_x.h"
#include "include/gausspt_ref.h"
#include "include/solution.h"
c
c     determinant of the Jacobian and the volume
c
      real*8 detj, dv
c
c     nodes and nodal coordinates of the vertices of a triangle
c
      integer nd1, nd2, nd3
      real*8 vertices(NSD,NODES_PER_TRIANGLE)
c
      real*8 sumux, sumuy, epsx(1,NSD+1), sig(1,NSD+1)
c
      real*8 l2norm, h1norm, exact_l2norm, num_l2norm
      real*8 energynorm, exact_energynorm, num_energynorm
c
      real*8 xstart, ystart, xend, yend
      integer no_of_points
c
c     integer counters
c
      integer i, j, tri_no, ngqp
c
      do 10 i = 1,no_of_nodes
c
c     output displacements at nodes
c
         write(50,197)u(NSD*i-1),u(NSD*i)
         write(80,198)coord(1,i),coord(2,i),u(NSD*i-1),u(NSD*i)
10    continue
c
      do 20 tri_no = 1,no_of_triangles
c 
c     get pgauss(3,nint) and wgauss(nint) - natural coordinates and
c     associated weights
c
         nint = tri_quad(tri_no)
         call gaussq_triangle(nint,pgauss,wgauss)
c
c     get nodal coordinates of the nodes at the vertex of triangle
c     "tri_no"
c
         nd1 = tri_connect(1,tri_no)
         nd2 = tri_connect(2,tri_no)
         nd3 = tri_connect(3,tri_no)
         vertices(1,1) = coord(1,nd1)
         vertices(2,1) = coord(2,nd1)
         vertices(1,2) = coord(1,nd2)
         vertices(2,2) = coord(2,nd2)
         vertices(1,3) = coord(1,nd3)
         vertices(2,3) = coord(2,nd3)
c
         do 30 ngqp = 1,nint
            xi1 = pgauss(1,ngqp)
            xi2 = pgauss(2,ngqp)
            xi3 = pgauss(3,ngqp)
c 
#ifdef FE
	    if (ngqp .eq. 1) then
	       xi1 = 1.d0
	       xi2 = 0.d0
	       xi3 = 0.d0
            elseif (ngqp .eq. 2) then
	       xi1 = 0.d0
	       xi2 = 1.d0
	       xi3 = 0.d0
            elseif (ngqp .eq. 3) then
	       xi1 = 0.d0
	       xi2 = 0.d0
	       xi3 = 1.d0
            endif
#endif
c
c     get Jacobian detj and the global coordinates (xbar,ybar) 
c     associated with the quadrature point xi1,xi2,xi3
c
            call compute_xgJ_triangle(vertices,detj)
            dv = 0.5*detj*wgauss(ngqp)
c
#ifdef FE
c
c     compute FE shape functions and derivatives
c
            call feshape2d_tri(tri_no,n,ngbrs_list,phi)
#else
c
c     compute number of natural neighbors for the Gauss point,
c     the neighbor list in ngbrs_list(), the NEM shape functions
c     (n-n coordinates) phi_I and its derivatives phi_I,j
c 
  	    call nemshape2d(n,ngbrs_list,phi)
#endif
	    sumux = 0.d0
	    sumuy = 0.d0
	    do 40 i = 1,n
               j = ngbrs_list(i)
               sumux = sumux + phi(1,i)*u(NSD*j-1)
               sumuy = sumuy + phi(1,i)*u(NSD*j)
40          continue
c
c     get strains at the Gauss points
c
            call strain2d(n,ngbrs_list,phi,u,epsx)
c
c     get stresses at the Gauss points
c
            call stress2d(epsx,sig)
c
c     output displacements, strains and stresses
c
            write(81,198)xbar,ybar,sumux,sumuy
            write(82,199)xbar,ybar,(epsx(1,j),j=1,3)
            write(83,199)xbar,ybar,(sig(1,j),j=1,3)
30       continue
20    continue
c
c     format statements
c
197   format(2E18.10)
198   format(4E18.10)
199   format(5E12.5)
c
c     norm calculations
c
#ifdef PLATE_HOLE
c
c     circular hole in an infinite plate under tension
c
      call l2energy_norm(l2norm,exact_l2norm,num_l2norm,
     &                   energynorm,exact_energynorm,num_energynorm)
c
      write(6,*)'Error norm calculations for plate with a hole problem'
      write(6,*)'_____________________________________________________'
      write(6,*)
      write(6,*)'L2-error norm is:',l2norm
      write(6,*)'Exact [!] L2-norm is:',exact_l2norm
      write(6,*)'Numerical L2-norm is:',num_l2norm
      write(6,*)'Error in energy-norm is:',energynorm
      write(6,*)'Exact [!] energy-norm is:',exact_energynorm
      write(6,*)'Numerical energy-norm is:',num_energynorm
#elif PATCH_EBC
c
c     displacement patch test
c
      call norms(l2norm,h1norm)
      write(6,*)
      write(6,*)'Error norm calculations for the patch test'
      write(6,*)'__________________________________________'
      write(6,*)
      write(6,*)'L2-error norm for the displacements: ',l2norm
      write(6,*)'H^1-error norm for the displacements: ',h1norm
#endif
c
#ifdef FE
#else
#ifdef PLATE_HOLE
c
c     Settings for the plate with a hole problem
c     ==========================================
c     xstart = 1.d-4
c     ystart = 1.001
c     xend = 1.d-4
c     yend = 4.9999
c     no_of_points = 240 ! 1 < R < 5
c
c************************************************
c
      write(*,*)'Enter xstart and ystart'
      read(*,*)xstart,ystart
      write(*,*)'Enter xend and yend'
      read(*,*)xend,yend
      write(*,*)'Enter no_of_points (> 1)'
      read(*,*)no_of_points
      call fieldout(xstart,ystart,xend,yend,no_of_points)
#endif
#endif
c
      return
      end
c
c**************************************************************
c     subroutine norms(l2norm,h1norm)
c     Purpose
c     =======
c     Compute the L2- and H1-error norms for the patch test
c
c**************************************************************
c
      subroutine norms(l2norm,h1norm)
      implicit none
c
c     include files for the norm calculations (NEM program)
c
#include "include/params2d.h"
#include "include/node.h"
#include "include/neighbors.h"
#include "include/triangle.h"
#include "include/quadrature.h"
#include "include/gausspt_x.h"
#include "include/gausspt_ref.h"
#include "include/solution.h"
c
      real*8 l2norm, h1norm
c
c     determinant of the Jacobian and the volume
c
      real*8 detj, dv
c
c     nodes and nodal coordinates of the vertices of a triangle
c
      integer nd1, nd2, nd3
      real*8 vertices(NSD,NODES_PER_TRIANGLE)
c
c     norm calculations
c
      external get_uexact, get_graduexact
      real*8 get_uexact, get_graduexact
      real*8 uerr1, uerr2
      real*8 grad_ux(NSD), grad_uy(NSD), graduerr(NSD,NSD)
c      
c     local variables and integer counters
c
      real*8 sumux, sumuy
      integer i, j, tri_no, ngqp
c
      l2norm = 0.0d0
      h1norm = 0.0d0
c
      do 10 tri_no = 1,no_of_triangles
c
c     get pgauss(3,nint) and wgauss(nint) - natural coordinates and
c     associated weights
c
#ifdef QUAD
         nint = QUAD
#else
         nint = tri_quad(tri_no)
#endif
         call gaussq_triangle(nint,pgauss,wgauss)
c
c     get nodal coordinates of the nodes at the vertex of triangle
c     "tri_no"
c
         nd1 = tri_connect(1,tri_no)
         nd2 = tri_connect(2,tri_no)
         nd3 = tri_connect(3,tri_no)
         vertices(1,1) = coord(1,nd1)
         vertices(2,1) = coord(2,nd1)
         vertices(1,2) = coord(1,nd2)
         vertices(2,2) = coord(2,nd2)
         vertices(1,3) = coord(1,nd3)
         vertices(2,3) = coord(2,nd3)
c
         do 20 ngqp = 1,nint
            xi1 = pgauss(1,ngqp)
            xi2 = pgauss(2,ngqp)
            xi3 = pgauss(3,ngqp)
c
c     get Jacobian detj and the global coordinates (xbar,ybar) 
c     associated with the quadrature point xi1,xi2,xi3
c
            call compute_xgJ_triangle(vertices,detj)
            dv = 0.5*detj*wgauss(ngqp)
c
#ifdef FE
c
c     compute FE shape functions and derivatives
c
            call feshape2d_tri(tri_no,n,ngbrs_list,phi)
#else
c
c     compute number of natural neighbors for the Gauss point,
c     the neighbor list in ngbrs_list(), the NEM shape functions
c     (n-n coordinates) phi_I and its derivatives phi_I,j
c
  	    call nemshape2d(n,ngbrs_list,phi)
#endif
            sumux = 0.0d0
            sumuy = 0.0d0
	    do 30 i = 1,n
               j = ngbrs_list(i)
               sumux = sumux + phi(1,i)*u(NSD*j-1)
               sumuy = sumuy + phi(1,i)*u(NSD*j)
30          continue
c
            uerr1 = sumux - get_uexact(xbar,ybar,1)
            uerr2 = sumuy - get_uexact(xbar,ybar,2)
            l2norm = l2norm + (uerr1**2 + uerr2**2)*dv
            h1norm = h1norm + (uerr1**2 + uerr2**2)*dv
c
c     get displacement gradients at the Gauss points [H^1 norm
c     calculation]
c
            do 40 i = 1,NSD
	       grad_ux(i) = 0.0d0
	       grad_uy(i) = 0.0d0
40          continue
c
            do 50 i = 1,n
	       j = ngbrs_list(i)
    	       grad_ux(1) = grad_ux(1) + phi(2,i)*u(NSD*j-1)
	       grad_ux(2) = grad_ux(2) + phi(3,i)*u(NSD*j-1)
	       grad_uy(1) = grad_uy(1) + phi(2,i)*u(NSD*j)
	       grad_uy(2) = grad_uy(2) + phi(3,i)*u(NSD*j)
50          continue
c
            do 60 j = 1,NSD
	       graduerr(1,j) = grad_ux(j) - get_graduexact(1,j,
     &            	       xbar,ybar)
	       graduerr(2,j) = grad_uy(j) - get_graduexact(2,j,
     &                         xbar,ybar)
60          continue
c
            do 70 i = 1,NSD
	       do 80 j = 1,NSD
		  h1norm = h1norm + graduerr(i,j)*graduerr(i,j)*dv
80             continue
70          continue
c
20       continue
10    continue
c
      l2norm = dsqrt(l2norm)
      h1norm = dsqrt(h1norm)
c
      return
      end
c
c**************************************************************
c     real*8 function get_uexact(x,y,dir)
c     Purpose
c     =======
c     Exact solution (displacements) for the patch test
c
c**************************************************************
c
      real*8 function get_uexact(x,y,dir)
      implicit none
c
      integer dir
      real*8 x, y
c
c     Exact (Analytical) solution
c
      if (dir .eq. 1) then
	 get_uexact = x
      else if (dir .eq. 2) then
	 get_uexact = y
      else
	 write(6,*)'Error in calling function get_uexact()'
	 stop
      endif
c
      return
      end
c
c**************************************************************
c     real*8 function get_graduexact(component,dir,x,y)
c     Purpose
c     =======
c     Compute the exact solution (displacement gradients) for
c     the patch test
c
c**************************************************************
c
      real*8 function get_graduexact(component,dir,x,y)
      implicit none
c
      integer component, dir
      real*8 x, y
c
c     Exact (Analytical) solution for the displacement gradients
c
      if (component .eq. 1) then
	 if (dir .eq. 1) then
	    get_graduexact = 1.d0
	 elseif (dir .eq. 2) then
	    get_graduexact = 0.d0
	 else
	   write(6,*)'Error in calling function get_graduexact()'
	   stop
         endif
      elseif (component .eq. 2) then
	 if (dir .eq. 1) then
	    get_graduexact = 0.d0
	 elseif (dir .eq. 2) then
	    get_graduexact = 1.d0
	 else
	   write(6,*)'Error in calling function get_graduexact()'
	   stop
         endif
      else
	 write(6,*)'Error in calling function get_epsexact()'
	 stop
      endif
c
      return
      end
c
c**************************************************************
c     subroutine l2energy_norm(l2norm,exact_l2norm,num_l2norm,
c    &                         energynorm,exact_energynorm,
c    &                         num_energynorm))
c     Purpose
c     =======
c     Compute the error in L2- and energy-norm, the exact value
c     (numerically computed) of the error norms, and also the
c     numerical [FE/NEM] error norms.
c
c**************************************************************
c
      subroutine l2energy_norm(l2norm,exact_l2norm,num_l2norm,
     &                         energynorm,exact_energynorm,
     &                         num_energynorm)
      implicit none
c
c     include files for the energy-norm calculations
c
#include "include/params2d.h"
#include "include/node.h"
#include "include/neighbors.h"
#include "include/triangle.h"
#include "include/quadrature.h"
#include "include/gausspt_x.h"
#include "include/gausspt_ref.h"
#include "include/solution.h"
c
      real*8 l2norm, exact_l2norm, num_l2norm
      real*8 energynorm, exact_energynorm, num_energynorm
c
c     determinant of the Jacobian and the volume
c
      real*8 detj, dv
c
c     nodes and nodal coordinates of the vertices of a triangle
c
      integer nd1, nd2, nd3
      real*8 vertices(NSD,NODES_PER_TRIANGLE)
c
c     norm calculations
c
      real*8 ux_e, uy_e
      real*8 sig(1,NSD+1), epsx(1,NSD+1)
      real*8 sig_e(1,NSD+1), epsx_e(1,NSD+1)
c      
c     local variables and integer counters
c
      real*8 sumux, sumuy
      integer i, j, tri_no, ngqp
c
      l2norm = 0.0d0
      exact_l2norm = 0.0d0
      num_l2norm = 0.0d0
      energynorm = 0.0d0
      exact_energynorm = 0.0d0
      num_energynorm = 0.0d0
c
      do 10 tri_no = 1,no_of_triangles
c
c     get pgauss(3,nint) and wgauss(nint) - natural coordinates and
c     associated weights
c
#ifdef QUAD
         nint = QUAD
#else
         nint = tri_quad(tri_no)
#endif
         call gaussq_triangle(nint,pgauss,wgauss)
c
c     get nodal coordinates of the nodes at the vertex of triangle
c     "tri_no"
c
         nd1 = tri_connect(1,tri_no)
         nd2 = tri_connect(2,tri_no)
         nd3 = tri_connect(3,tri_no)
         vertices(1,1) = coord(1,nd1)
         vertices(2,1) = coord(2,nd1)
         vertices(1,2) = coord(1,nd2)
         vertices(2,2) = coord(2,nd2)
         vertices(1,3) = coord(1,nd3)
         vertices(2,3) = coord(2,nd3)
c
         do 20 ngqp = 1,nint
            xi1 = pgauss(1,ngqp)
            xi2 = pgauss(2,ngqp)
            xi3 = pgauss(3,ngqp)
c
c     get Jacobian detj and the global coordinates (xbar,ybar) 
c     associated with the quadrature point xi1,xi2,xi3
c
            call compute_xgJ_triangle(vertices,detj)
            dv = 0.5*detj*wgauss(ngqp)
c
#ifdef FE
c
c     compute FE shape functions and derivatives
c
            call feshape2d_tri(tri_no,n,ngbrs_list,phi)
#else
c
c     compute number of natural neighbors for the Gauss point,
c     the neighbor list in ngbrs_list(), the NEM shape functions
c     (n-n coordinates) phi_I and its derivatives phi_I,j
c
  	    call nemshape2d(n,ngbrs_list,phi)
#endif
            sumux = 0.0d0
            sumuy = 0.0d0
	    do 30 i = 1,n
               j = ngbrs_list(i)
               sumux = sumux + phi(1,i)*u(NSD*j-1)
               sumuy = sumuy + phi(1,i)*u(NSD*j)
30          continue
c
            call strain2d(n,ngbrs_list,phi,u,epsx)
c
            call stress2d(epsx,sig)
c
            call exact_disp_hole(ux_e,uy_e)
            call exact_stress_strain_hole(sig_e,epsx_e)
c
            l2norm = l2norm + dv*((sumux - ux_e)**2 + (sumuy - uy_e)**2)
            exact_l2norm = exact_l2norm + dv*(ux_e**2 + uy_e**2)
            num_l2norm = num_l2norm + dv*(sumux**2 + sumuy**2)
c
            energynorm = energynorm + dv*((epsx(1,1) - epsx_e(1,1))*
     &                                    (sig(1,1) - sig_e(1,1)) + 
     &                                    (epsx(1,2) - epsx_e(1,2))*
     &                                    (sig(1,2) - sig_e(1,2)) + 
     &                                    (epsx(1,3) - epsx_e(1,3))*
     &                                    (sig(1,3) - sig_e(1,3)))
            exact_energynorm = exact_energynorm + 
     &                         (epsx_e(1,1)*sig_e(1,1) + 
     &                          epsx_e(1,2)*sig_e(1,2) +
     &                          epsx_e(1,3)*sig_e(1,3))*dv
            num_energynorm = num_energynorm + 
     &                       (epsx(1,1)*sig(1,1) + 
     &                        epsx(1,2)*sig(1,2) +
     &                        epsx(1,3)*sig(1,3))*dv
20       continue          
10    continue          
c
      l2norm = dsqrt(l2norm)
      exact_l2norm = dsqrt(exact_l2norm)
      num_l2norm = dsqrt(num_l2norm)
      energynorm = dsqrt(energynorm/2.d0)
      exact_energynorm = dsqrt(exact_energynorm/2.d0)
      num_energynorm = dsqrt(num_energynorm/2.d0)
c
      return
      end
c
c*****************************************************************
c     subroutine exact_disp_hole(ux_e,uy_e)
c     Purpose
c     =======
c     Subroutine to compute the exact displacements for the plate
c     with a hole problem
c
c*****************************************************************
c 
      subroutine exact_disp_hole(ux_e,uy_e)
      implicit none
c
#include "include/params2d.h"
#include "include/gausspt_x.h"
#include "include/plate_hole.h"
c
      real*8 ux_e, uy_e
      real*8 sig(1,NSD+1), epsx(1,NSD+1)
      real*8 r, theta, aovr
      real*8 fac
      real*8 ymbar, nubar, kappa
      common /equiv_material/ ymbar, nubar, kappa
c
      r = dsqrt(xbar**2 + ybar**2)
      theta = datan2(ybar,xbar)
      aovr = a/r
      fac = 0.25d0*a*(1.d0 + nubar)/ymbar
      ux_e = fac*((kappa + 1.0)*dcos(theta)/aovr +
     &            2.0*aovr*((1.0 + kappa)*dcos(theta) +
     &            dcos(3.0*theta)) - 2.0*(aovr**3)*dcos(3.0*theta))
      uy_e = fac*((kappa - 3.0)*dsin(theta)/aovr +
     &            2.0*aovr*((1.0 - kappa)*dsin(theta) +
     &            dsin(3.0*theta)) - 2.0*(aovr**3)*dsin(3.0*theta))
c
      return
      end
c
c*****************************************************************
c     subroutine exact_stress_strain_hole(sig,epsx)
c     Purpose
c     =======
c     Subroutine to compute the exact stress and strain fields
c     from the exact stress field solution for the plate with
c     a hole problem
c
c     sig(k,j)	: the (k-1)th derivative of the stress component j 
c     epsx(k,j)	: the (k-1)th derivative of the strain component j 
c
c     sig(1,1)  = sige11 sig(1,2)  = sig22  sig(1,3)  = sig12
c     epsx(1,1) = eps11  epsx(1,2) = eps22  epsx(1,3) = 2*epsx12
c
c*****************************************************************
c 
      subroutine exact_stress_strain_hole(sig,epsx)
      implicit none
c
#include "include/params2d.h"
#include "include/gausspt_x.h"
#include "include/plate_hole.h"
c
      real*8 sig(1,NSD+1), epsx(1,NSD+1)
      real*8 r, theta, aovr
c
c     2D Cauchy stresses
c
      r = dsqrt(xbar**2 + ybar**2)
      theta = datan2(ybar,xbar)
      aovr = a/r
      sig(1,1) = 1. - aovr**2 * (1.5 * dcos(2.*theta)
     &           + dcos(4.*theta)) + 1.5*aovr**4 * dcos(4.*theta)
      sig(1,2) = -aovr**2 * (0.5 * dcos(2.*theta) - dcos(4.*theta))
     &           - 1.5*aovr**4 * dcos(4.*theta)
      sig(1,3) = -aovr**2 * (0.5 * dsin(2.*theta) + dsin(4.*theta))
     &           + 1.5*aovr**4 * dsin(4.*theta)
c
c     compute 2D strains
c
      call strain2d_from_stress2d(sig,epsx)
c
      return
      end
c
c**************************************************************
c     subroutine strain2d_from_stress2d(sig,epsx)
c     Purpose
c     =======
c     Subroutine to compute the 2D strains from the 2D stresses
c
c**************************************************************
c
      subroutine strain2d_from_stress2d(sig,epsx)
      implicit none
c
#include "include/params2d.h"
#include "include/material.h"
c
      real*8 d
      real*8 sig(1,NSD+1), epsx(1,NSD+1)
c
c     2D small strains [engineering shear strain gamma_12 is
c     stored in epsx(1,3)]
c
      d = d1*d1 - d2*d2
      epsx(1,1) = (sig(1,1)*d1 - sig(1,2)*d2)/d
      epsx(1,2) = (sig(1,2)*d1 - sig(1,1)*d2)/d
      epsx(1,3) = sig(1,3)/d3
c      
      return
      end
c
c****************************************************************
c     subroutine fieldout(xstart,ystart,xend,yend,points)
c     Purpose
c     =======
c     Subroutine to output field variables along a line between
c     two points (x1,y1) and (x2,y2)
c
c****************************************************************
c
      subroutine fieldout(xstart,ystart,xend,yend,points)
      implicit none
c
#include "include/params2d.h"
#include "include/node.h"
#include "include/neighbors.h"
#include "include/triangle.h"
#include "include/gausspt_x.h"
#include "include/material.h"
#include "include/solution.h"
c
      real*8 xstart, ystart, xend, yend
      integer points
      real*8 xn1, yn1, xn2, yn2
c
      integer point_no
      real*8 dt, tht, r, c, s, fac
      real*8 sumux, sumuy, ur_e
      real*8 epsx(1,3), sig(1,3), epsx_e(1,3), sig_e(1,3)
      real*8 ymbar, nubar, kappa
      common /equiv_material/ ymbar, nubar, kappa 
      integer consistency_flag
      common /c_flag/ consistency_flag
c
      integer i, j
c
      open(2,file='PATH_LINE.dat')
      open(7,file='u_out.dat')
      open(8,file='e_out.dat')
      open(9,file='s_out.dat')
      dt = 1.d0/dble(points - 1)
      tht = datan2(yend-ystart,xend-xstart)
      c = dcos(tht)
      s = dsin(tht)
c
      do 10 point_no = 1,points 
 	 xbar = xstart + (point_no-1)*dt*(xend - xstart)
	 ybar = ystart + (point_no-1)*dt*(yend - ystart)
c
c     compute number of natural neighbors for the Gauss point,
c     the neighbor list in ngbrs_list(), the NEM shape functions
c     (n-n coordinates) phi_I and its derivatives phi_I,j
c
         call nemshape2d(n,ngbrs_list,phi)
	 sumux = 0.0d0
	 sumuy = 0.0d0
	 do 20 i = 1,n
            j = ngbrs_list(i)
            sumux = sumux + phi(1,i)*u(NSD*j-1)
            sumuy = sumuy + phi(1,i)*u(NSD*j)
20       continue
c
c     get strains at the output points
c
         call strain2d(n,ngbrs_list,phi,u,epsx)
c
c     get stresses at the output points
c
         call stress2d(epsx,sig)
c
c     compute exact soln for the stress
c
      r = dsqrt(xbar*xbar + ybar*ybar)
c
         call exact_stress_strain_hole(sig_e,epsx_e)
c
c     compute 2D strains
c
         call strain2d_from_stress2d(sig_e,epsx_e)
c
c     output exact and numerical stresses along the line
c     joining (x1,y1) to (x2,y2)
c
      if (consistency_flag .eq. 0) then
         write(2,199)r,(sig(1,j),j=1,3),(sig_e(1,j),j=1,3)
      endif
c
10    continue
c
c     format statement 
c
197   format(4E15.7,I3)
198   format(6E15.7)
199   format(7E15.7)
299   format(8E15.7)
c
      close(2)
      return
      end
