SUBROUTINE moment (rect,badpts,npt,usage,rotate,buffer,hx,hy) 1
!
  REAL :: rect(28),badpts(2,npt),evec1(2),evec2(2),buffer,              &
           rotate
!
! input parameters:
!  badpts      = x,y coords of flagged badpts grouped into clusters
!                are in the first two rows
!  npt         = num. of badpts. in the cluster.
!  buffer      = expand new grids with buffer zone of size buffer
!                this prevents degenerate point rectangles.
! output parameters:
!  usage       = ratio of flagged to unflagged badpts. in new grid
!                measures goodness of fit and clustering
!    rect( )      = stores some info. for grid created herein.
!                sometimes rect = rnode, sometimes = temp. array.
!                depending on calling prog. (grdfit or expand)
!
!
!  moment = compute 2nd moments about the mean of flagged badpts.
! fit an ellipse. the eigenvectors determine rectangle orientation.
! same some info., even tho. usage might be low and rect. scrapped.
!
  rn = FLOAT(npt)
  xmean = 0.0
  ymean = 0.0
  cross = 0.0
  xsq   = 0.0
  ysq   = 0.0
  DO ipt = 1, npt
    xmean = xmean + badpts(1,ipt)
    ymean = ymean + badpts(2,ipt)
    xsq   = xsq   + badpts(1,ipt)**2
    ysq   = ysq   + badpts(2,ipt)**2
    cross = cross + badpts(1,ipt)*badpts(2,ipt)
  END DO
!
  xmean = xmean / rn
  ymean = ymean / rn
  xm2   = xmean ** 2
  ym2   = ymean ** 2
  xsq   = xsq / rn
  ysq   = ysq / rn
  cross = cross / rn - xmean * ymean
!
  tl    = xsq - xm2
  br    = ysq - ym2
!
  eval1 = ( (tl+br) + SQRT((tl-br)**2 + 4.0*cross**2) ) / 2.0
  eval2 = ( (tl+br) - SQRT((tl-br)**2 + 4.0*cross**2) ) / 2.0
!
! matrix = <tl, cross;  cross, br>.
! determine eigenvectors.  normalize for ease in later computation.
! make 1st evector lie in southeast quadrant.
! make 2nd evector lie in northeast quadrant.
!
  IF (ABS(cross) <= rotate) GO TO 30
  evec1(2) = (eval1-tl) / cross
  evec2(2) = (eval2-tl) / cross
  IF (evec1(2) <= evec2(2)) GO TO 20
  temp     = evec1(2)
  evec1(2) = evec2(2)
  evec2(2) = temp
  20       CONTINUE
  fac1     = SQRT(1.0+evec1(2)**2)
  evec1(1) = 1.0 / fac1
  evec1(2) = evec1(2) / fac1
  fac2     = SQRT(1.0 + evec2(2)**2)
  evec2(1) = 1.0 / fac2
  evec2(2) = evec2(2) / fac2
  GO TO 40
  30   CONTINUE
  evec1(1) = 1.0
  evec2(1) = 0.0
  evec1(2) = 0.0
  evec2(2) = 1.0
  40   CONTINUE
!
! compute length of enclosing rectangles to include all flagged badpts.
! take projection (dotproduct) of bad pt. on evectors to do this.
! expand on all sides by 'buffer' zone.  this also prevents
!  case of degenerate rectangle (line).
!
  emx1 = badpts(1,1)*evec1(1) + badpts(2,1)*evec1(2)
  emn1 = emx1
  emx2 = badpts(1,1)*evec2(1) + badpts(2,1)*evec2(2)
  emn2 = emx2
  IF (npt == 1) GO TO 90
  DO ipt = 2, npt
    dot1  = badpts(1,ipt)*evec1(1) + badpts(2,ipt)*evec1(2)
    dot2  = badpts(1,ipt)*evec2(1) + badpts(2,ipt)*evec2(2)
    IF (dot1 <= emx1) GO TO 50
    emx1 = dot1
    50       IF (dot1 >= emn1) GO TO 60
    emn1 = dot1
    60       IF (dot2 <= emx2) GO TO 70
    emx2 = dot2
    70       IF (dot2 >= emn2) CYCLE
    emn2 = dot2
  END DO
  90   emx1  = emx1 + buffer
  emx2  = emx2 + buffer
  emn1  = emn1 - buffer
  emn2  = emn2 - buffer
!
! from length of the sides, determine the 4 rect. corners
! number the corners clockwise
!
  rect(1) = emn1*evec1(1)+emn2*evec2(1)
  rect(2) = emn1*evec1(2)+emn2*evec2(2)
  rect(3) = emn1*evec1(1)+emx2*evec2(1)
  rect(4) = emn1*evec1(2)+emx2*evec2(2)
  rect(5) = emx1*evec1(1)+emx2*evec2(1)
  rect(6) = emx1*evec1(2)+emx2*evec2(2)
  rect(7) = emx1*evec1(1)+emn2*evec2(1)
  rect(8) = emx1*evec1(2)+emn2*evec2(2)
!
!  set some other fields here that would be difficult elsewhere.
!
  rect(16)  = emx1  +.00001
  rect(17)  = emn1  -.00001
  rect(18)  = emx2  +.00001
  rect(19)  = emn2  -.00001
  rect(12)  = evec1(1)
  rect(13)  = evec1(2)
  rect(14) = evec2(1)
  rect(15) = evec2(2)
  dist = SQRT( (rect(7)-rect(1))**2                                     &
              +(rect(8)-rect(2))**2 )
  rect( 21)=(rect(7)-rect(1))/dist
  rect( 22)=(rect(8)-rect(2))/dist
!
! compute volume cutoff ratio
!
  iside1 = (hx+emx1-emn1-2.0*buffer+.00001) / hx
  iside2 = (hy+emx2-emn2-2.0*buffer+.00001) / hy
  gpall  = iside1 * iside2
  usage  = rn / gpall
!
  RETURN
END SUBROUTINE moment
!
! --------------------------------------------------------------------
!


SUBROUTINE outtre(mlev,ddsc,dmpc) 1,3
!
  INCLUDE 'nodal.inc'
  INCLUDE 'agricpu.inc'
  INCLUDE 'agricst.inc'
  LOGICAL :: ddsc,dmpc
!
! ******************************************************************
! OUTTRE - OUTPUT SUBTREE
! INPUT PARAMETERS:
!    MLEV   - PTR TO SUBTREE TO OUTPUT I.E., START AT LEVEL(MLEV)
! TREE IS OUTPUT FROM 'LEVEL' TO FINEST LEVEL.
! ******************************************************************
!
!
!  dump out generic grid data if asked for
!
  IF(ddsc) THEN
    CALL dmpdsc
    CALL dmpstr
  END IF
!
! dump out tree if needed
!
  IF( .NOT. verbose1 ) RETURN
!
  time = rnode(20,mlev)/60.

  WRITE(6,1) time
  1     FORMAT(//,1X,'  GRID STRUCTURE AT  ',f10.3,' MINUTES ')
!
  level = node(4, mlev)
!
  IF( verbose6 ) WRITE(6,'(''  LEVEL, MLEV,LFINE '',3I5)')level,mlev,lfine
!
  23000 IF(.NOT. (level <= lfine))GO TO 23001
  mptr    = lstart(level)
  23002     IF(.NOT. (mptr /= 0))GO TO 23003
  CALL outmsh(mptr)
!c               if (dmpc) call dmpcnst(mptr)
  mptr = node(10, mptr)
  GO TO 23002
  23003     CONTINUE
  level = level + 1
  GO TO 23000
  23001 CONTINUE
!
  RETURN
END SUBROUTINE outtre
!
! --------------------------------------------------------------------
!


SUBROUTINE outmsh(mptr) 1
  INCLUDE 'nodal.inc'
!
!  this routine prints out the grid information
!
  IF( (node(4,mptr) < 1) ) RETURN
  time = rnode(20,mptr)/60.
  WRITE(6,1)mptr,node(4,mptr),time
  1     FORMAT(/,1X,' ******** GRID ',i3,'  LEVEL ',i3,                 &
               ' AT ',f10.3,' MINUTES ')
  WRITE(6,2) (rnode(i,mptr),i=1,8),rnode(29,mptr),rnode(30,mptr)
  2     FORMAT(/,10X,'corners        x              y   ',              &
            /,13X,'1',3X,2(2X,e13.6),/,13X,'2',3X,2(2X,e13.6),          &
            /,13X,'3',3X,2(2X,e13.6),/,13X,'4',3X,2(2X,e13.6),          &
            //,' z at bottom = ',e13.6,' z at top = ',e13.6,/)
  WRITE(6,3) (rnode(i,mptr),i=9,11),rnode(31,mptr)
  3     FORMAT(3X,'  dx =  ',e12.5,'   dy =  ',e12.5,                   &
               /,3X,'  dt =  ',e12.5,'   dz =  ',e12.5)
  WRITE(6,4)rnode(21,mptr),rnode(22,mptr),                              &
            (rnode(i,mptr),i=12,19),                                    &
            (rnode(i,mptr),i=23,28)
  4     FORMAT(3X,'  cos = ',e12.5,'   sin = ',e12.5,                   &
               //,'  misc grid description information  ',              &
               /,4(/,2X,4(1X,e12.5)) )
  WRITE(6,5) node(5,mptr),node(6,mptr),node(14,mptr),                   &
             node(7,mptr),node(8,mptr),node(16,mptr)
  5     FORMAT(/,'  NODAL INFORMATION ',/,                              &
                 1X,'   nx = ',i4,'  ny = ',i4,'  nz = ',i4,/,          &
                 1X,'   bigstep file = ',i3,' smlstp file = ',i3,       &
                    ' constants file = ',i3)
  WRITE(6,6) node(10,mptr),node(12,mptr),node(13,mptr)
  6     FORMAT(1X,'   next grid on level = ',i3,                        &
                  ' previous grid on level = ',i3,/,                    &
                1X,'   number of timesteps taken = ',i4)
  RETURN
END SUBROUTINE outmsh
!
! --------------------------------------------------------------------
!


SUBROUTINE outch( a,tmp,m,n,kk )
  DIMENSION a(m,n,kk),tmp(1)
!
!  plot out a few levels
!
  kmid = kk/2
  jmid = n/2
  imid = m/2
!  call parray(a(1,1,kmid),m,n)
  DO k=1,kk
    DO j=1,n
      tmp(j+(k-1)*n) = a(imid,j,k)
    END DO
  END DO
!  call parray(tmp,n,kk)
  DO k=1,kk
    DO i=1,m
      tmp(i+(k-1)*m) = a(i,jmid,k)
    END DO
  END DO
!  call parray(tmp,m,kk)
  RETURN
END SUBROUTINE outch
!
! --------------------------------------------------------------------
!


SUBROUTINE dmpdsc 1
  INCLUDE 'agrigrid.inc'
  INCLUDE 'grddsc.inc'
  INCLUDE 'agricst.inc'
!
!  print the generic grid description information if wanted
!
  IF( .NOT. verbose2 ) RETURN
!
  WRITE(6,'(''  GENERIC GRID INFORMATION '')')
  WRITE(6,1) nsint,nsreal,nx1d,ny1d,nz1d,                               &
             nxy2d,nxz2d,nyz2d,nxyz3d
!
  IF (nxy2d > 0) THEN
    DO i=1,nxy2d
      WRITE(6,2) i
      WRITE(6,3) stgxy(1,i),stgxy(2,i)
      WRITE(6,4) idmxy(1,i),idmxy(2,i)
      WRITE(6,5) ipkxy(i)
      WRITE(6,6) iupxy(1,i),iupxy(2,i)
      WRITE(6,7) ibdxy(1,i),ibdxy(2,i),ibdxy(3,i)
    END DO
  END IF
!
  IF (nxz2d > 0) THEN
    DO i=1,nxz2d
      WRITE(6,12) i
      WRITE(6,3) stgxz(1,i),stgxz(2,i)
      WRITE(6,4) idmxz(1,i),idmxz(2,i)
      WRITE(6,5) ipkxz(i)
      WRITE(6,6) iupxz(1,i),iupxz(2,i)
      WRITE(6,7) ibdxz(1,i),ibdxz(2,i)
    END DO
  END IF
!
  IF (nyz2d > 0) THEN
    DO i=1,nyz2d
      WRITE(6,22) i
      WRITE(6,3) stgyz(1,i),stgyz(2,i)
      WRITE(6,4) idmyz(1,i),idmyz(2,i)
      WRITE(6,5) ipkyz(i)
      WRITE(6,6) iupyz(1,i),iupyz(2,i)
      WRITE(6,7) ibdyz(1,i),ibdyz(2,i)
    END DO
  END IF
!
  IF (nxyz3d > 0) THEN
    DO i=1,nxyz3d
      WRITE(6,32) i
      WRITE(6,33) stgxyz(1,i),stgxyz(2,i),stgxyz(3,i)
      WRITE(6,34) idmxyz(1,i),idmxyz(2,i),idmxyz(3,i)
      WRITE(6,5)  ipkxyz(i)
      WRITE(6,6)  iupxyz(1,i),iupxyz(2,i)
      WRITE(6,7)  ibdxyz(1,i),ibdxyz(2,i),ibdxyz(3,i)
    END DO
  END IF
!
  1   FORMAT( /,2X,' nsint  = ',i5,                                     &
              /,2X,' nsreal = ',i5,                                     &
              /,2X,' nx1d   = ',i5,                                     &
              /,2X,' ny1d   = ',i5,                                     &
              /,2X,' nz1d   = ',i5,                                     &
              /,2X,' nxy2d  = ',i5,                                     &
              /,2X,' nxz2d  = ',i5,                                     &
              /,2X,' nyz2d  = ',i5,                                     &
              /,2X,' nxyz3d = ',i5  )
!
  2   FORMAT( /,'  data for 2-D xy variable ',i4 )
!
  3   FORMAT( '  x staggering ',f6.3,'  y staggering ',f6.3 )
  4   FORMAT( '  x max dimension ' ,i4,'  y max dimension ',i4 )
  5   FORMAT( '  pack variable - ',i2)
  6   FORMAT( '  update - ',i2,'  other vector ',i2)
  7   FORMAT( '  boundi - ',i2,'  other vector ',i2,                    &
              '  t-dt field - ',i2)
!
  12   FORMAT( /,'  data for 2-D xz variable ',i4 )
  22   FORMAT( /,'  data for 2-D yz variable ',i4 )
  32   FORMAT( /,'  data for 3-D xyz variable ',i4 )
  33   FORMAT( '  x staggering ',f6.3,'  y staggering ',f6.3,           &
               '  z staggering ',f6.3                         )
  34   FORMAT( '  x max dimension ' ,i4,'  y max dimension ',i4,        &
               '  z max dimension ',i4 )
  RETURN
END SUBROUTINE dmpdsc
!


SUBROUTINE settmp 4
  INCLUDE 'agrialloc.inc'
  INCLUDE 'manage.inc'
  INCLUDE 'agricst.inc'
!
  ntemp = lfree(2,1)
  IF( verbose4 ) THEN
    WRITE(6,'(''  IN SETTMP '')')
    WRITE(6,'(''  NTEMP = '',I7)') ntemp
  END IF
  RETURN
END SUBROUTINE settmp
!
! --------------------------------------------------------------------
!


SUBROUTINE dmpstr 1
  INCLUDE 'nodal.inc'
  INCLUDE 'agrialloc.inc'
  INCLUDE 'agrigrid.inc'
  INCLUDE 'grddsc.inc'
  INCLUDE 'agricst.inc'
!
!  output storage locations for the grids if necessary
!
  IF( .NOT. verbose4 ) RETURN
!
  level = 1
  10    IF(level > lfine) GO TO 999
  mptr = lstart(level)
  20    IF(mptr == 0) GO TO 100
!
  WRITE(6,'(/,''  STORAGE LOCATIONS FOR GRID '',I4)') mptr
  WRITE(6,1001) ipint(mptr),ipreal(mptr),ips1d(mptr),ipx(mptr),         &
                ipy(mptr),ipz(mptr)
!
  IF(nxy2d > 0) THEN
    WRITE(6,'('' STORAGE LOCATIONS  FOR XY VARIABLES '')')
    DO i=1,nxy2d
      WRITE(6,1002) i,ipxy(i,mptr)
    END DO
  END IF
  IF(nxz2d > 0) THEN
    WRITE(6,'('' STORAGE LOCATIONS  FOR XZ VARIABLES '')')
    DO i=1,nxz2d
      WRITE(6,1002) i,ipxz(i,mptr)
    END DO
  END IF
  IF(nyz2d > 0) THEN
    WRITE(6,'('' STORAGE LOCATIONS  FOR YZ VARIABLES '')')
    DO i=1,nyz2d
      WRITE(6,1002) i,ipyz(i,mptr)
    END DO
  END IF
  IF(nxyz3d > 0) THEN
    WRITE(6,'('' STORAGE LOCATIONS  FOR XYZ VARIABLES '')')
    DO i=1,nxyz3d
      WRITE(6,1002) i,ipxyz(i,mptr)
    END DO
  END IF
  1001  FORMAT(2X,5(2X,i9))
  1002  FORMAT(2X,i5,2X,i9)
!
  mptr = node(10,mptr)
  GO TO 20
  100   level = level+1
  GO TO 10
  999   RETURN
END SUBROUTINE dmpstr
!


SUBROUTINE setrot(mlevel) 1
  INCLUDE 'nodal.inc'
!
!
!  SETROT:  SET THE ROTATION FIELDS IN RNODE FOR ALL GRIDS AT
!        THE SAME LEVEL AS MPTR.  THESE FIELDS ARE
!        COS AND SIN OF THE ANGLE OF ROTATION OF THE GRID.
!        (THESE ARE SET ELSEWHERE).
!        ALSO SET SLOPE/INTERCEPT FORM FOR ALL 4 LINES OF THE RECT.
!       SET THE MAX/MIN PROJECTION ON EACH EIGENVECTOR FIELD.
!
  mptr = mlevel
  23000 IF(.NOT. (mptr /= 0))GO TO 23001
  dist = SQRT( (rnode(7,mptr)-rnode(1,mptr))**2                         &
              +(rnode(8,mptr)-rnode(2,mptr))**2 )
  rnode( 21,mptr)=(rnode(7,mptr)-rnode(1,mptr))/dist
  rnode( 22,mptr)=(rnode(8,mptr)-rnode(2,mptr))/dist
!
  rnode( 12,mptr) =  rnode(21,mptr)
  rnode( 13,mptr) =  rnode(22,mptr)
  rnode( 14,mptr) = -rnode(22,mptr)
  rnode( 15,mptr) =  rnode(21,mptr)
!
  delx = rnode(3,mptr) - rnode(1,mptr)
  dely = rnode(4,mptr) - rnode(2,mptr)
  IF(.NOT. (ABS(delx) > .0001))GO TO 23002
  rnode(23,mptr) = dely / delx
  GO TO 23003
  23002     CONTINUE
  rnode(23,mptr) = 1000000000
  23003     CONTINUE
  IF(.NOT. (ABS(rnode(23,mptr)) > .0001))GO TO 23004
  rnode(24,mptr) = -1.0 / rnode(23,mptr)
  GO TO 23005
  23004     CONTINUE
  rnode(24,mptr) = - 1000000000
  23005     CONTINUE
  rnode(25,mptr) = rnode(2,mptr)-rnode(23,mptr)*rnode(1,mptr)
  rnode(26,mptr) = rnode(2,mptr)-rnode(24,mptr)*rnode(1,mptr)
  rnode(27,mptr) = rnode(6,mptr) - rnode(23,mptr)*rnode(5,mptr)
  rnode(28,mptr) = rnode(6,mptr) - rnode(24,mptr)*rnode(5,mptr)
!
  rnode( 16,mptr) = rnode(5,mptr)*rnode(12,mptr) +rnode(6,mptr)         &
                       *rnode(13,mptr)+.0001
  rnode( 17,mptr) = rnode(1,mptr)*rnode(12,mptr) +rnode(2,mptr)         &
                       *rnode(13,mptr)-.0001
  rnode( 18,mptr) = rnode(5,mptr)*rnode(14,mptr) +rnode(6,mptr)         &
                       *rnode(15,mptr)+.0001
  rnode( 19,mptr) = rnode(1,mptr)*rnode(14,mptr) +rnode(2,mptr)         &
                       *rnode(15,mptr)-.0001
  mptr = node(10, mptr)
  GO TO 23000
  23001 CONTINUE
!
  RETURN
END SUBROUTINE setrot
!


SUBROUTINE cpyfld( mptr,ifldin,ifldout,idim ) 6,2
  INCLUDE 'nodal.inc'
  INCLUDE 'agrialloc.inc'
!
! this routine copies the data in xy or xyz field
! ifldin onto ifldout and replaces ifldout.
!
  nx = node(5,mptr)
  ny = node(6,mptr)
  nz = node(14,mptr)
!
  IF (idim == 2) THEN
    inptrs = igtnxy( mptr,  ifldin, 1 )
    inptrr = igtnxy( mptr, ifldout, 1 )
!
    DO ij=1,nx*ny
      a(inptrr+ij-1) = a(inptrs+ij-1)
    END DO
!
    CALL retnxy( mptr, ifldout, 1, inptrr, .true. )
!
  ELSE IF (idim == 3) THEN
!
    inptrs = igtxyz( mptr,  ifldin, 1 )
    inptrr = igtxyz( mptr, ifldout, 1 )
    DO ijk=1,nx*ny*nz
      a(inptrr+ijk-1) = a(inptrs+ijk-1)
    END DO
    CALL retxyz( mptr, ifldout, 1, inptrr, .true. )
!
  ELSE
    WRITE(6,'(''  ERROR, NOT 2-D OR 3-D IN CPYFLD '',                   &
    &             4I9)') mptr,ifldin,ifldout,idim
  END IF
  RETURN
END SUBROUTINE cpyfld


SUBROUTINE addgrd( mptr ) 2,7

  INCLUDE 'agrialloc.inc'
  INCLUDE 'agrigrid.inc'
  INCLUDE 'grddsc.inc'
  INCLUDE 'agricst.inc'
  CHARACTER (LEN=80) :: filenm
!
!  this routine takes a grid stored at unit iunit
!  and adds it to the end of the grid storage
!
  CALL resett
  CALL setstr( mptr )
  CALL settmp
  in = ipint( mptr )
!
!  get the proper file
!
  CALL mkflnm( 'tmpdmp',' ','.z',0.,mptr,filenm,length )
  CALL fexist( filenm(1:length),length,.true.,.false. )
  iunit = igtunit( dum )
  OPEN( UNIT=iunit,FILE=filenm(1:length),FORM='unformatted',            &
        STATUS='old'                                        )
  IF( verbose3 ) WRITE(6,*) ' read grid ',mptr,' from tmp file ',       &
                filenm(1:length),', io through unit ',iunit
  REWIND( iunit )
!
! Yuhe: when EXBC turned on, we must add the EXBC arrays for grid 1.
!
  IF ( lexbc == 1 .AND. mptr == 1 ) THEN
    nwords = ipexbc(nexbc3d+1,mptr) - ipint(mptr)
  ELSE
    nwords = ipxyz(nxyz3d+1,mptr) - ipint(mptr)
  END IF

  IF( verbose3 ) WRITE(6,*) ' calling rdngrd ',mptr,iunit,nwords
  CALL rdngrd( iunit,a(in),nwords )
!
!  close and delete file, return unit number
!
  CLOSE( UNIT=iunit,STATUS='keep' )
  CALL retnunit( iunit )
!
! we're finished here
!
  RETURN
END SUBROUTINE addgrd
!


SUBROUTINE dmpgrd( mptr,saved ) 3,7

  LOGICAL :: saved
  INCLUDE 'agrialloc.inc'
  INCLUDE 'agrigrid.inc'
  INCLUDE 'grddsc.inc'
  INCLUDE 'agricst.inc'
  CHARACTER (LEN=80) :: filenm

!
!  this routine removes grid storage from the incore
!  array and, if desired, writes it out to file iunit
!  which, at present, will be the grid number plus 39
!
  IF ( saved ) THEN
!
! dump data to disk file
!
    CALL mkflnm( 'tmpdmp',' ','.z',0.,mptr,filenm,length )
    iunit = igtunit( dum )
    CALL fexist( filenm(1:length),length,.false.,.false. )
    OPEN( UNIT=iunit,FILE=filenm(1:length),FORM='unformatted',          &
          STATUS='new'                                        )
    IF( verbose3 ) WRITE(6,*) ' dump grid ',mptr,' to tmp file ',       &
                 filenm(1:length),', io through unit ',iunit
    REWIND( iunit )
!
! Yuhe: when EXBC turned on, we must add the EXBC arrays for grid 1.
!
    IF ( lexbc == 1 .AND. mptr == 1 ) THEN
      nwords = ipexbc(nexbc3d+1,mptr) - ipint(mptr)
    ELSE
      nwords = ipxyz(nxyz3d+1,mptr) - ipint(mptr)
    END IF

    IF( verbose3 ) WRITE(6,*) ' calling wrtgrd ',mptr,iunit,nwords
    CALL wrtgrd( iunit, a(ipint(mptr)), nwords )
!
!  close and delete file, return unit number
!
    CLOSE( UNIT=iunit,STATUS='keep' )
    CALL retnunit( iunit )
!
  END IF
!
! free up this space
!
  IF( verbose4 ) WRITE(6,*) ' reclaiming space for grid ',mptr
  CALL resett
  CALL reclam( ipint(mptr),nwords )
  CALL settmp
!
! we're finished here
!
  RETURN
END SUBROUTINE dmpgrd
!


SUBROUTINE rdngrd( iunit,a,nwords ) 2
  DIMENSION a(nwords)
  READ(iunit) a
  RETURN
END SUBROUTINE rdngrd
!


SUBROUTINE wrtgrd( iunit,a,nwords ) 2
  DIMENSION a(nwords)
  WRITE(iunit) a
  RETURN
END SUBROUTINE wrtgrd
!


SUBROUTINE dmplvl( level ) 1,1
  INCLUDE 'nodal.inc'
  INCLUDE 'agricst.inc'
!
!  this dumps out all the grids at level 'level'
!  into individual files using the dmpgrd subroutine
!
!  grids in lstart are dumped, starting with the last one.
!  it is assumed that the present level is the highest
!  existing in memory
!
  mptr = lback(level)
  IF(mptr == 0) THEN
    WRITE(6,*) ' error in dmplvl, no grids on level ',level
    STOP
  END IF
  5     IF( mptr == 0) GO TO 99
  IF( verbose3 ) WRITE(6,*) ' dumping grid - level ',mptr,level
  CALL dmpgrd( mptr,.true. )
  mptr = node(12,mptr)
  GO TO 5
!
  99    RETURN
END SUBROUTINE dmplvl
!


SUBROUTINE exchng( level ) 1,3
  INCLUDE 'nodal.inc'
  INCLUDE 'agricst.inc'
!
!  this subroutine calls the boundary condition exchange
!  for overlapping fine grids
!
!  this version doesn't do an overlap check first.
!  later version should, we just need to put in overlap
!  test to do so
!
  cpu0 = f_cputime()
  mptr = lstart( level )
  5     IF(mptr == 0) GO TO 999
!
  mptre = node(10,mptr)
  10    IF(mptre == 0) GO TO 20
!
!  do boundary condition exchange between mptr and mptre
!
!  here's where the overlap test will be
!
  WRITE(6,*) ' b.c. exchange between ',mptr,mptre
  CALL exchbc( mptr,mptre )
  mptre = node(10,mptre)
  GO TO 10
!
  20    CONTINUE
!
!  we've gone through the list for grid mptr, now go
!  to the next grid and go through the rest of the list
!
  mptr = node(10,mptr)
  GO TO 5
!
  999   CONTINUE
  cpu_bndexch = cpu_bndexch + f_cputime() - cpu0
!
!  we're finished here
!
  RETURN
END SUBROUTINE exchng
!


SUBROUTINE dmpall( time ) 1,1
  INCLUDE 'nodal.inc'
  INCLUDE 'agricst.inc'
!
!  this dumps out all the grids for a restart
!  into individual files using the dumpgrid subroutine.
!  this routine does not clear internal storage, thus integrations
!  can be continued after this dump, as opposed to the
!  dumps performed by dmplvl, which does clear the internal storage
!
!  grids in lback are dumped, starting with the last one.
!
  level = 1
  mptr = lback(level)
!
  IF(mptr == 0) THEN
    WRITE(6,*) ' error in dmpall, no grids on level 1 '
    STOP
  END IF
!
  5    IF( mptr == 0) GO TO 10
!
  IF( verbose3 ) WRITE(6,*) ' dumping grid - level ',mptr,level
  CALL dumpgrid( mptr,time )
  mptr = node(12,mptr)
  GO TO 5
  10    level = level + 1
  mptr = lback(level)
  IF(level <= lfine) GO TO 5
!
  99    RETURN
END SUBROUTINE dmpall
!


SUBROUTINE readall( time ) 1,1
  INCLUDE 'nodal.inc'
  INCLUDE 'agricst.inc'
!
!  this reads in all the grids for a restart
!  from individual files using the readgrid subroutine.
!  this routine assumes that the internal storage is
!  already configured, thus no storage management needs to be
!  done (as opposed to reading in a new grid using addgrd,
!  which does allocate new mwmory for the grid)
!
  level = 1
  mptr = lback(level)
!
  IF(mptr == 0) THEN
    WRITE(6,*) ' error in readall, no grids on level 1 '
    STOP
  END IF
!
  5   IF( mptr == 0) GO TO 10
!
  IF( verbose3 ) WRITE(6,*) ' reading grid - level ',mptr,level
  CALL readgrid( mptr,time )
  mptr = node(12,mptr)
  GO TO 5
  10    level = level + 1
  mptr = lback(level)
  IF(level <= lfine) GO TO 5
!
  99    RETURN
END SUBROUTINE readall
!
!


SUBROUTINE dumpgrid( mptr,time ) 1,4

  IMPLICIT NONE

  INCLUDE 'agrialloc.inc'
  INCLUDE 'agrigrid.inc'
  INCLUDE 'grddsc.inc'
  INCLUDE 'agricst.inc'
  INCLUDE 'globcst.inc'

  CHARACTER (LEN=80) :: filenm
  INTEGER :: length

  INTEGER :: mptr
  INTEGER :: nwords
  INTEGER :: iunit, igtunit
  INTEGER :: machsti

  REAL :: dum
  REAL :: time
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  this routine dumps out grid data for an individual
!  grid for later restart.
!
  CALL mkflnm( runname(1:lfnkey),' ','.t',time,mptr,filenm,length )

  iunit = igtunit( dum )
  WRITE(6,*) ' restart dump for run ',runname(1:lfnkey),                &
             ', grid ',mptr,' at time ',time

  IF( verbose3 ) WRITE(6,*) ' dumping to file ',filenm(1:length),       &
                ' through unit ',iunit
  CALL fexist( filenm(1:length),length,.false.,.true. )

  OPEN( UNIT=iunit,FILE=filenm(1:length),FORM='unformatted',            &
        STATUS='new' )
  REWIND( iunit )
  WRITE(iunit) runname,filenm,machst,lfnkey,mptr,time

!
! Yuhe: when EXBC turned on, we must add the EXBC arrays for grid 1.
!
  IF ( lexbc == 1 .AND. mptr == 1 ) THEN
    nwords = ipexbc(nexbc3d+1,mptr) - ipint(mptr)
  ELSE
    nwords = ipxyz(nxyz3d+1,mptr) - ipint(mptr)
  END IF

  IF( verbose3 ) WRITE(6,*) ' calling wrtgrd ',mptr,iunit,nwords
  CALL wrtgrd( iunit, a(ipint(mptr)), nwords )
  CLOSE( UNIT=iunit,STATUS='keep' )

  CALL retnunit( iunit )
  IF( verbose3 ) WRITE(6,*) ' closed ',filenm(1:length),                &
                ' on unit ',iunit
!
! we're finished dumping this grid data
!
  RETURN
END SUBROUTINE dumpgrid
!


SUBROUTINE readgrid( mptr,time ) 1,4

!  implicit none

  INCLUDE 'agrialloc.inc'
  INCLUDE 'agrigrid.inc'
  INCLUDE 'grddsc.inc'
  INCLUDE 'agricst.inc'
!  include 'globcst.inc'

  CHARACTER (LEN=80) :: filenm

!
!  this routine reads in grid data for an individual
!  grid for restart.
!
  CALL mkflnm( runold(1:nmlntho),' ','.t',time,mptr,filenm,length )
  iunit = igtunit( dum )
  WRITE(6,*) ' restart read from run ',runold(1:nmlntho),               &
             ', grid ',mptr,' at time ',time
  IF( verbose3 ) WRITE(6,*) ' reading from file ',filenm(1:length),     &
                ' through unit ',iunit
  CALL fexist( filenm(1:length),length,.true.,.false. )

  OPEN( UNIT=iunit,FILE=filenm(1:length),FORM='unformatted',            &
        STATUS='old' )
  REWIND( iunit )
  READ( iunit ) runnami,filenmi,machsti,nmlnthi,mptri,timei
!
!  here we assume everything is ok, we could check the
!  runname and filenames etc... that come from the file
!  if that is deemed necessary
!
!
! Yuhe: when EXBC turned on, we must add the EXBC arrays for grid 1.
!
  IF ( lexbc == 1 .AND. mptr == 1 ) THEN
    nwords = ipexbc(nexbc3d+1,mptr) - ipint(mptr)
  ELSE
    nwords = ipxyz(nxyz3d+1,mptr) - ipint(mptr)
  END IF

  IF( verbose3 ) WRITE(6,*) ' calling rdngrd ',mptr,iunit,nwords
  CALL rdngrd( iunit, a(ipint(mptr)), nwords )

  CLOSE( UNIT=iunit,STATUS='keep' )
  CALL retnunit( iunit )
  IF( verbose3 ) WRITE(6,*) ' closed ',filenm(1:length),                &
                ' on unit ',iunit
!
! we're finished reading this grid data
!
  RETURN
END SUBROUTINE readgrid
!
!


SUBROUTINE iosetup 1
  COMMON / unitnmb / iunitn(10),nleft,ntot
  DO i=1,10
    iunitn(i) = 40+i
  END DO
  nleft = 10
  ntot = 10

  WRITE(6,'(/2(/1x,a)/)')                                               &
      'Fortran I/O units 40-50 reserved for AGR interface.',            &
      'Do not use them for other purposes.'

  RETURN
END SUBROUTINE iosetup
!
!


  INTEGER FUNCTION igtunit( dum )
  COMMON / unitnmb / iunitn(10),nleft,ntot
  INCLUDE 'agricst.inc'
!
  IF(nleft == 0) THEN
    WRITE(6,*) ' no more units left for io in iunitn list '
    WRITE(6,*) ' error stop '
  END IF
!
  igtunit = iunitn( nleft )
  IF( verbose3 ) WRITE(6,*) ' igtunit getting unit= ',igtunit
  nleft = nleft-1
  RETURN
  END FUNCTION igtunit
!
!


SUBROUTINE retnunit( iunit ) 8
  COMMON / unitnmb / iunitn(10),nleft,ntot
!
  nleft = nleft+1
  IF(nleft > ntot) THEN
    WRITE(6,*) ' something is screwed, there are '
    WRITE(6,*) ' more units in iunitn than when we started!! '
    WRITE(6,*) ' returning unit ',iunit
    WRITE(6,*) ' ntot, nleft = ',ntot,nleft
    DO i=1,ntot
      WRITE(6,2) i,iunitn(i)
    END DO
  END IF
  2    FORMAT(2X,i5,1X,i5)
!
  iunitn(nleft) = iunit
  RETURN
END SUBROUTINE retnunit
!
!-----------------------------------------------------------------------
!


SUBROUTINE mkflnm( instr1,instr2,instr3,                                & 7
           time,grid,outstr,length )

  IMPLICIT NONE
  REAL :: time
  INTEGER :: i, j, k, length, grid, tlength,glength
  CHARACTER (LEN=*) :: instr1, instr2, instr3, outstr
  CHARACTER (LEN=80) :: tmp, chmake
  EXTERNAL chmake
  DATA tlength,glength/ 5, 2 /

!-----------------------------------------------------------------------
! Find lengths of the strings
!-----------------------------------------------------------------------
  i = INDEX( instr1,' ' ) - 1
  IF(i == -1) i = LEN( instr1 )
  j = INDEX( instr2,' ' ) - 1
  IF(j == -1) j = LEN( instr2 )
  k = INDEX( instr3,' ' ) - 1
  IF(k == -1) k = LEN( instr3 )
  length = 0
!
  IF( i > 0 ) THEN
    outstr(1:i) = instr1(1:i)
    length = i
  END IF
!
  IF( j > 0 ) THEN
    outstr(length+1:length+j) = instr2(1:j)
    length = length+j
  END IF
!
!-----------------------------------------------------------------------
! if time is < 0, then all we want is this string
!-----------------------------------------------------------------------
!
  IF( time < 0 ) RETURN
!
!-----------------------------------------------------------------------
! If TIME > 0 add time to the string
!-----------------------------------------------------------------------
!
  tmp = chmake( IFIX(time),tlength )
!
  IF(k > 0) THEN
    outstr(i+j+1:) = instr3(1:k)//tmp(1:tlength)
  ELSE
    outstr(i+j+1:) = tmp(1:tlength)
  END IF
!
  length = length + tlength + k
!
!-----------------------------------------------------------------------
! Create the grid no string if grid number greater than zero
!-----------------------------------------------------------------------
!
  IF( grid > 0 ) THEN
    tmp = chmake( grid,glength )
!      outstr(length+1:) = '.g'//tmp(1:glength)
!      length = length + glength + 2
    outstr(length+1:) = '.'//tmp(1:glength)
    length = length + glength + 1
  END IF
!
!  we're finished
!
  WRITE(6,'(a,a)') 'The file name created is ',outstr(1:length)

  999    RETURN
END SUBROUTINE mkflnm
!
!


  FUNCTION chmake( numb,length )

  IMPLICIT NONE
  INTEGER :: numb,length,numbc,i
  CHARACTER (LEN=*) :: chmake
  CHARACTER (LEN=16) :: tmp

!
!  check for degenerate cases
!
  IF(length <= 0)  THEN
    WRITE(6,*) ' length is <= zero in chmake ',length
    WRITE(6,*) ' error exit '
    STOP
  END IF

!
  IF(numb == 0) THEN
!
!  make the string
!
    DO i=1,length
      chmake(i:i) = '0'
    END DO
    RETURN
  END IF
!
  numbc = INT(ALOG10(ABS(FLOAT(numb)))) + 1
  IF( numb < 0 ) numbc = numbc+1
!
  IF( numbc > length ) THEN
    WRITE(6,*) ' number too large in chmake ',                          &
               ' for given length ',length,numb,numbc
    WRITE(6,*) ' error exit '
    STOP
  END IF
!
  IF( numbc > 16 ) THEN
    WRITE(6,*) ' number has more than 16 chars ',                       &
               ' in chmake, you must be kidding!!! '
    WRITE(6,*) ' error exit '
    STOP
  END IF
!
!  make the string
!
  DO i=1,length
    chmake(i:i) = '0'
  END DO
!
  WRITE(tmp,2) numb
  1      FORMAT(i1)
  2      FORMAT(i16)
  chmake(length-numbc+1:length) = tmp(16-numbc+1:16)
!
!  we're finished here
!
  RETURN
  END FUNCTION chmake
!
!


SUBROUTINE cleanf( mptrs ) 2,2
  INCLUDE 'nodal.inc'
  INCLUDE 'agricst.inc'
  CHARACTER (LEN=80) :: filenm
  LOGICAL :: existf
!
!  this routine closes the files for the temp grid
!  data dumps used in regrid, mptr is the first grid
!  on  a level to be dumped
!
  iunit = igtunit( dum )
  mptr = mptrs
  10   IF( mptr == 0) GO TO 99
  CALL mkflnm( 'tmpdmp',' ','.z',0.,mptr,filenm,length )
  INQUIRE(FILE=filenm(1:length),EXIST=existf)
!
  IF(existf) THEN
    OPEN( UNIT=iunit,FILE=filenm(1:length),                             &
          FORM='unformatted',STATUS='old'   )
    CLOSE(UNIT=iunit,STATUS='delete')
  ELSE
    WRITE(6,*) ' in cleanf, expected to find file ',                    &
                filenm(1:length),', yet does not exist '
  END IF
  mptr = node(10,mptr)
  GO TO 10
  99   CONTINUE
  CALL retnunit( iunit )
  RETURN
END SUBROUTINE cleanf
!


SUBROUTINE fexist( filenm,length,exst,incr ) 6
  CHARACTER (LEN=*) :: filenm
  INTEGER :: length
  LOGICAL :: exst,incr,existf
!
!  this subroutine checks on the existence of the file
!  filenm(1:length).
!
!  exst: true if the file should exist, false if it shouldn't
!     exist
!  incr: if true and exst is false but the file exists, then
!     return a new name for a file that doesn't exist.
!     here we just patch on a version number to the file
!
!  if the file should exist but doesn't, error exit
!   if the file shouldn't exist but does and no increment
!   is wanted, then we error exit
!   later we'll build in an erase option taht will allow
!   removal of existing files that we do not want.
!
!  first we check for file that should exist
!
  IF( exst ) THEN
    INQUIRE( FILE=filenm(1:length),EXIST=existf )
    IF(.NOT.existf ) THEN
      WRITE(6,*) ' exected to find file ',filenm(1:length)
      WRITE(6,*) ' file does not exist, error exit '
      STOP
    ELSE
      RETURN
    END IF
  END IF
!
!  now for file that shouldn't exist
!
  inc = 0
  10   INQUIRE( FILE=filenm(1:length),EXIST=existf )
!
  IF (.NOT. existf ) RETURN
!
  IF (.NOT. incr ) THEN
    WRITE(6,*) ' cannot overwrite pre-existing file ',                  &
                filenm(1:length),', error exit '
    STOP
  END IF
!
!  here we put a number onto the end of the file
!
  inc = inc+1
  IF(inc == 1) THEN
    WRITE(filenm(length+1:length+2),15) inc
    length = length+2
  ELSE IF(inc < 10) THEN
    WRITE(filenm(length-1:length  ),15) inc
  ELSE IF(inc > 10) THEN
    WRITE(6,*) '  to many version numbers of file ',                    &
               filenm(1:length-2),', error exit '
    STOP
  END IF
  GO TO 10
!
  999  RETURN
  15   FORMAT('.',i1)
END SUBROUTINE fexist


SUBROUTINE filzero(a,n) 4
  REAL :: a(n)

  DO i=1,n
    a(i) = 0.0
  END DO

  RETURN
END SUBROUTINE filzero
!
!-----------------------------------------------------------------------
!
! Created by Louis Wicker, November 1992
! Writes out the cpu usage for the model
!
! Modified for ARPS by Ming Xue, 12/10/1992.
!
!-----------------------------------------------------------------------
!


SUBROUTINE prtcpu(flag,UNIT) 2

  INCLUDE 'agricpu.inc'
  INTEGER :: flag, UNIT

!-----------------------------------------------------------------------
!  Flag = 0, initialize values to zero
!-----------------------------------------------------------------------

  IF( flag == 0 ) THEN

    cpu_main       = 0.0
    cpu_advect     = 0.0
    cpu_smlstep    = 0.0
    cpu_mixing     = 0.0
    cpu_filter     = 0.0
    cpu_cloud      = 0.0
    cpu_boundary   = 0.0
    cpu_copy       = 0.0
    cpu_init0      = 0.0
    cpu_usrout     = 0.0
    cpu_regrid     = 0.0
    cpu_bndexch    = 0.0
    cpu_bndcint    = 0.0
    cpu_update     = 0.0
    cpu_pack       = 0.0
    cpu_unpack     = 0.0
    cpu_simulation = 0.0

!-----------------------------------------------------------------------
!  Flag=1, write out cpu information
!-----------------------------------------------------------------------

  ELSE

    cpu_solver = cpu_advect+cpu_smlstep+cpu_mixing+cpu_cloud            &
               + cpu_boundary+cpu_copy+cpu_main+cpu_filter

    cpu_interface = cpu_simulation - cpu_solver - cpu_usrout            &
                  - cpu_init0

    WRITE(UNIT,999)
    WRITE(UNIT,101) 'ADAPTIVE MESH CLOUD MODEL CPU STATISTICS'
    WRITE(UNIT,999)

    WRITE(UNIT,101) 'Total CPU for SIMULATION  is ',cpu_simulation
    WRITE(UNIT,999)

    WRITE(UNIT,101) 'Total CPU for INTERFACE   is ',cpu_interface
    WRITE(UNIT,101) 'Total CPU for REGRID      is ',cpu_regrid
    WRITE(UNIT,101) 'Total CPU for EXCHBC      is ',cpu_bndexch
    WRITE(UNIT,101) 'Total CPU for UPDBC       is ',cpu_bndcint
    WRITE(UNIT,101) 'Total CPU for UPDATE      is ',cpu_update
    WRITE(UNIT,101) 'Total CPU for PACKING     is ',cpu_pack
    WRITE(UNIT,101) 'Total CPU for UNPACKING   is ',cpu_unpack
    WRITE(UNIT,101) 'Total CPU for INIT        is ',cpu_init0
    WRITE(UNIT,101) 'Total CPU for USROUT      is ',cpu_usrout
    WRITE(UNIT,999)

    WRITE(UNIT,101) 'Total CPU for SOLVER      is ',cpu_solver
    WRITE(UNIT,101) 'Total CPU for CLOUD3D     is ',cpu_main
    WRITE(UNIT,101) 'Total CPU for SUB SMLSTEP is ',cpu_smlstep
    WRITE(UNIT,101) 'Total CPU for SUB ADVECT  is ',cpu_advect
    WRITE(UNIT,101) 'Total CPU for SUB CLOUD   is ',cpu_cloud
    WRITE(UNIT,101) 'Total CPU for SUB MIXING  is ',cpu_mixing
    WRITE(UNIT,101) 'Total CPU for SUB FILTER  is ',cpu_filter
    WRITE(UNIT,101) 'Total CPU for SUB BOUNDS  is ',cpu_boundary
    WRITE(UNIT,101) 'Total CPU for SUB COPY    is ',cpu_copy

    WRITE(UNIT,999)
    101     FORMAT(1X,a,f15.5)
    999     FORMAT(/79('-')/)

  END IF

  RETURN
END SUBROUTINE prtcpu



SUBROUTINE writint( icnst,ncnst , UNIT)

  INTEGER :: ncnst , UNIT
  INTEGER :: icnst(ncnst)

  DO i=1,ncnst
    WRITE(UNIT, '(1x,2i10)') i, icnst(i)
  END DO

  RETURN
END SUBROUTINE writint


SUBROUTINE readint( icnst,ncnst , UNIT)

  INTEGER :: ncnst , UNIT
  INTEGER :: icnst(ncnst)

  DO i=1,ncnst
    READ(UNIT, '(1x,2i10)') ii, icnst(i)
    WRITE(6, '(1x,2(a,i10))') 'Integer No. ',i,'=',icnst(i)
  END DO

  RETURN
END SUBROUTINE readint