SUBROUTINE www_graph (ishow, nl, nr, nz, t, z)
!= w3g.f90 (10Jun2010) Updated (2010) gnuplot driver
 PARAMETER (l=512)
 IMPLICIT DOUBLE PRECISION (a-h, o-z)
 CHARACTER (LEN=l) :: gnuplot, workdir, arg, gnudata, gnuload, image
 LOGICAL test
 COMMON /www_graph_local/ test
 DIMENSION t(nl:nr), z(nz,nl:nr); DATA idatafile/ 321 /; test = .FALSE.
 gnuplot = '/usr/bin/gnuplot'
!! workdir = '/var/www/html/cisti/tmp/'
!! workdir = '/afs/ist.utl.pt/users/3/8/ist11038/web/tmp/'
 workdir = './tmp/'
 image = 'gnu.png'
! ... Create a temporary, work directory
 IF (IARGC() == 1) THEN
    CALL GETARG (1, arg); workdir = TRIM(workdir) // TRIM(arg) // '/'
 END IF
97100 FORMAT (' W3G ', A, ' file:', T42, '| "', A, '" opened' / T42, '|', T4, A)
! ... Open the DATA file, "gnudata"
 gnudata = TRIM(workdir) // 'gnudata.dat'
 OPEN (idatafile, FILE=TRIM(gnudata), STATUS='UNKNOWN')
 IF (test) WRITE (*, 97100) 'Data', 'gnudata', TRIM(gnudata)
 CALL data_to_file (idatafile, nl, nr, nz, t, z, ishow)
! ... Create the 'gnuplot' LOAD file, "gnuload"
 gnuload = TRIM(workdir) // 'gnuload.dat'; loadfile = idatafile + 1
 OPEN (loadfile, FILE=gnuload, STATUS='UNKNOWN')
 IF (test) WRITE (*, 97100) 'Load', 'gnuload', TRIM(gnuload)
 image = TRIM(workdir) // image
 CALL gnuload_make (loadfile, TRIM(image), TRIM(gnudata))
! ... Now the graph
 CALL SYSTEM (gnuplot // ' ' // gnuload)
 CLOSE (idatafile, STATUS='DELETE'); CLOSE (loadfile, STATUS='DELETE')
 CALL SYSTEM ('chmod a+r ' // TRIM(workdir) // 'gnu.png')
 RETURN
 END SUBROUTINE !END!

 SUBROUTINE gnuload_make (loadfile, image, gnudata)
! ... Constructs 'gnuload', the load file ..............................
 IMPLICIT DOUBLE PRECISION (a-h, o-z)
 PARAMETER (l=512)
 CHARACTER (LEN=*) :: image, gnudata, line*(l)
 COMMON /www_graph_local/ test; LOGICAL test
 WRITE (loadfile, "(A)") 'set output "' // image // '"'
 IF (test) WRITE (*, 97100) image
97100 FORMAT (' W3G set output...   File:', T42, '| image' / T42, '|', T4, A)
! ... Get-put parameters BEFORE 'plot' or 'splot'
 DO; READ (*, "(A)", END=2900) line; line = ADJUSTL(line)
    ixp = INDEX(line, 'plot '); ixs = INDEX(line, 'splot ')
    IF (ixp == 1 .OR. ixs == 1) GOTO 3100
    WRITE (loadfile, "(A)") TRIM(line)
 END DO;
2900 WRITE (*, "(' *** Error: No plotting directive ***')"); STOP 
! ... Put plotting directive
3100 WRITE (loadfile, "(A)") line(1:5)//' "'//gnudata//'" \'
! ... Get-put parameters AFTER 'plot' or 'splot'
 DO; READ (*, "(A)", END=3900) line; line = ADJUSTL(line)
    WRITE (loadfile, "(A)") TRIM(line)
 END DO
3900 RETURN
 END SUBROUTINE !END!

 SUBROUTINE data_to_file (ifile, nl, nr, nz, t, z, ishow)
! ... ..................................................................
 IMPLICIT DOUBLE PRECISION (a-h, o-z)
 DIMENSION t(nl:nr), z(nz,nl:nr)
 DATA nzmax/ 50 / ! no. of curves
! ... Supply the data to file 'ifile'
 WRITE (ifile, "('#Result vectors', I3)") nz ! comment
 DO i=nl, nr; WRITE (ifile, "(1X, 50E20.9E3)") t(i), z(:MIN(nzmax,nz),i)
 END DO
 IF (ishow == 0) RETURN
 WRITE (*, "(' n_min, _max:', 2I8)") nl, nr
 WRITE (*, "(1X, 6X, 3X, 9X, 't,', 3X, 5X, 'z(''s):')")
 DO i=nl, nr; WRITE (*, 99920) i, t(i), z(:MIN(nzmax,nz),i)
99920 FORMAT (1X, I6, 3X, 1P, 50G15.7)
 END DO
 RETURN
 END SUBROUTINE !END!