
INTEGER FUNCTION istrlen(str)
  IMPLICIT REAL(a-h,o-z)
  CHARACTER str(*)

  istrlen = 0
  DO
     IF (str(istrlen+1) == ' ') EXIT
     istrlen = istrlen + 1
  END DO
END FUNCTION istrlen

SUBROUTINE append_num(str,num)
  IMPLICIT REAL (a-h,o-z)

  CHARACTER str(*)
  len = istrlen(str)

  nstp = num
  DO ipos = len+3,len+1,-1
     idigit = MOD(nstp,10)
     str(ipos:ipos) = CHAR(ICHAR('0') + idigit)
     nstp = nstp / 10
  END DO
END SUBROUTINE append_num


SUBROUTINE read_q(maxmx, maxmy, mbc, LS, fname,iunit,mxf,myf,tf ,q)
  IMPLICIT REAL(a-h,o-z)

  !! INCLUDE "parms.h"
  REAL, DIMENSION(1-mbc:maxmx+mbc,1-mbc:maxmy+mbc) :: q
  CHARACTER*LS fname, fname_tmp

  lenf = istrlen(fname)
  fname_tmp(1:LS)= ''
  fname_tmp(1:lenf) = fname(1:lenf)
  CALL append_num(fname_tmp,iunit)

  OPEN(10,file=fname_tmp)
  READ(10,*) tf, mxf, myf, mqn
  DO j = 1,myf
     DO i = 1,mxf
        READ(10,*) q(i,j)
     END DO
  END DO
  CLOSE(10)

END SUBROUTINE read_q

SUBROUTINE write_q(maxmx,maxmy, mbc, LS, mx,my,dta, t, fname,iunit)
  IMPLICIT REAL (a-h,o-z)

  DIMENSION dta(1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc)
  CHARACTER(len=LS) fname, fname_tmp
  fname_tmp(1:LS) = ''
  lenf = istrlen(fname)
  fname_tmp(1:lenf) = fname(1:lenf)

  IF (iunit >= 100) THEN
     CALL append_num(fname_tmp,iunit)
  END IF
  OPEN(iunit,file=fname_tmp)

  meqn = 1
  WRITE(iunit,1000) t,mx,my, meqn
1000 FORMAT(E16.8, 3I5)

  DO j=1,my
     DO i=1,mx
        !! # exponents with more than 2 digits cause problems
        !! reading into matlab... reset tiny values to zero:
        IF (ABS(dta(i,j)) <  1d-99) dta(i,j) = 0.d0
        WRITE(iunit,1001) dta(i,j)
     END DO
     WRITE(iunit,*) ' '
  END DO

1001 FORMAT(4e24.16)

  WRITE(iunit,*) ' '
  CLOSE(iunit)
  RETURN
END SUBROUTINE write_q
