* Sample program for gtool_history/gtool4
* 
* Solving diffusion equation 
*     du/dt = kappa d^2 u/dx^2
* for giving values of u at x=0 and 1. 
*
      program diffeq

      integer nx, nt, ndisp
        parameter( nx=30 )                       ! $B%0%j%C%I?t(B
        parameter( nt=200 )                      ! $B;~4V%9%F%C%W?t(B
        parameter( ndisp=10 )                    ! $B=PNO4V3V(B
      double precision  dx, dt, kappa
        parameter( dx=1.0/(nx-1))                ! $B%0%j%C%I4V3V(B
        parameter( dt=0.0005 )                   ! $B;~4V4V3V(B
        parameter( kappa=1.0D0 )                 ! $BG.3H;678?t(B
      double precision x(nx)                     ! $B:BI8JQ?t(B
      double precision temp(nx), temp1(nx)       ! $B29EY(B

      integer i,it                               ! $B:n6HJQ?t(B

      character*1  dims(2), units(2)
      character*6  xtypes(2)
      character*12 lnames(2)
      integer      dimsiz(2)
      data      dims    /'x','t'/
      data      dimsiz  /nx,0/
      data      xtypes  /'double','double'/
      data      lnames  /'X-coordinate','time        '/
      data      units   /'m','s'/

      do i=1,nx
        x(i)=dx*(i-1)                            ! $B:BI8CM@_Dj(B
      enddo

      do i=1,nx
        temp(i) = exp(-((x(i)-0.5)/0.1)**2)      ! $B=i4|CM@_Dj(B
      enddo

      call hscrea( 'diff77.nc', 'Diffusion equation',      ! $B%R%9%H%j!<:n@.(B 
     &             'Sample Fortran77 progam of gtool_history/gtool4', 
     &             'GFD_Dennou Club davis project', 
     &             2, dims, dimsiz, xtypes, lnames, units, 
     &             0.0, real(ndisp*dt)          )

      call hsp1d('x',x,nx)                                     ! $B<!85JQ?t=PNO(B
      call hsavar('temp',2,dims,'temperature','K','double')    ! $BJQ?tDj5A(B

      call hsp1d('temp',temp,nx)                               ! $BJQ?t=PNO(B
      call hsaatd('temp','+gt_user_davis_kappa',kappa,1)       ! $BB0@-=PNO(B

      do it=1,nt                                 ! $B;~4V@QJ,(B
         do i=2,nx-1
            temp1(i) = temp(i) 
     &               + kappa*(temp(i+1)-2*temp(i)+temp(i-1))
     &                   /dx**2*dt
         enddo
         do i=1,nx
            temp(i) = temp1(i)
         enddo

         if ( mod(nt,ndisp) .eq. 0 ) then
            call hsp1d('temp',temp,nx)                      ! $BJQ?t=PNO(B
         endif
      enddo

      call hsclse

      stop
      end
