!! ALL MODULES

!! subroutines for sparse matrix multiplication of Slater determinant coefficient vectors
!! the matrices are represented by derived type configptr or sparseptr depending on sparseopt
!! these derived types are handled using derived type walktype - their indexing etc is stored there

#include "Definitions.INC"

!!!!!!      BYPROC SUBROUTINES   !!!!!!!

module arbitrarymultmod
  use mpisubmod

contains

  subroutine arbitrary_sparsemult_singles_byproc(firstproc,lastproc,www,mattrans, &
       rvector,insmallvector,outsmallvector,mynumr,diagflag)
    use walkmod
    use sparse_parameters
    use fileptrmod
    implicit none
    type(walktype),intent(in) :: www
    integer,intent(in) :: firstproc,lastproc
    integer,intent(in) :: mynumr,diagflag
    DATATYPE,intent(in) :: insmallvector(mynumr,www%allbotconfigs(firstproc):www%alltopconfigs(lastproc)), &
         mattrans(www%singlematsize,www%botconfig:www%topconfig)
    DATATYPE,intent(out) :: outsmallvector(mynumr,www%botconfig:www%topconfig)
    DATAECS,intent(in) :: rvector(mynumr)
    DATATYPE :: outsum(mynumr)
    DATAECS :: myrvector(mynumr)
    integer :: config1,ihop

    if (sparseconfigflag.eq.0) then
       OFLWR "error no arbitrary_sparsemult if sparseconfigflag.eq.0"; CFLST
    endif

    if (www%topconfig-www%botconfig+1.eq.0) then
       return
    endif

    outsmallvector(:,:)=0d0

    if (www%alltopconfigs(lastproc).lt.www%allbotconfigs(firstproc)) then
       return
    endif

    if (diagflag.eq.0) then

!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(config1,ihop,outsum,myrvector)
       myrvector(:)=rvector(:)
!$OMP DO SCHEDULE(DYNAMIC)
       do config1=www%botconfig,www%topconfig
          outsum(:)=0d0

          do ihop=www%firstsinglehopbyproc(firstproc,config1),www%lastsinglehopbyproc(lastproc,config1)
             outsum(:)=outsum(:)+mattrans(ihop,config1) * &
                  insmallvector(:,www%singlehop(ihop,config1)) * myrvector(:)
          enddo
          outsmallvector(:,config1)=outsum(:)
       enddo
!$OMP END DO
!$OMP END PARALLEL

    else

!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(config1,myrvector)
       myrvector(:)=rvector(:)

!$OMP DO SCHEDULE(DYNAMIC)
       do config1=www%botconfig,www%topconfig
          if (www%singlehopdiagflag(config1).ne.0) then
             outsmallvector(:,config1)=mattrans(www%singlediaghop(config1),config1) * &
                  insmallvector(:,config1) * myrvector(:)
          endif
       enddo
!$OMP END DO
!$OMP END PARALLEL

    endif

  end subroutine arbitrary_sparsemult_singles_byproc

  subroutine arbitrary_sparsemult_doubles_byproc(firstproc,lastproc,www,mattrans,&
       rvector,insmallvector,outsmallvector,mynumr,diagflag)
    use fileptrmod
    use sparse_parameters
    use walkmod
    implicit none
    type(walktype),intent(in) :: www
    integer,intent(in) :: firstproc,lastproc,mynumr,diagflag
    DATATYPE,intent(in) :: insmallvector(mynumr,www%allbotconfigs(firstproc):www%alltopconfigs(lastproc)), &
         mattrans(www%doublematsize,www%botconfig:www%topconfig)
    DATATYPE,intent(out) :: outsmallvector(mynumr,www%botconfig:www%topconfig)
    DATAECS,intent(in) :: rvector(mynumr)
    DATATYPE :: outsum(mynumr)
    DATAECS :: myrvector(mynumr)
    integer :: config1,ihop

    if (sparseconfigflag.eq.0) then
       OFLWR "error no arbitrary_sparsemult_doubles if sparseconfigflag.eq.0"; CFLST
    endif

    if (www%topconfig-www%botconfig+1.eq.0) then
       return
    endif

    outsmallvector(:,:)=0d0

    if (www%alltopconfigs(lastproc).lt.www%allbotconfigs(firstproc)) then
       return
    endif

    if (diagflag.eq.0) then

!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(config1,ihop,outsum,myrvector)
       myrvector(:)=rvector(:)

!$OMP DO SCHEDULE(DYNAMIC)
       do config1=www%botconfig,www%topconfig
          outsum(:)=0d0
          do ihop=www%firstdoublehopbyproc(firstproc,config1),www%lastdoublehopbyproc(lastproc,config1)
             outsum(:)=outsum(:)+mattrans(ihop,config1) * &
                  insmallvector(:,www%doublehop(ihop,config1)) * myrvector(:)
          enddo
          outsmallvector(:,config1)=outsum(:)
       enddo
!$OMP END DO
!$OMP END PARALLEL

    else

!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(config1,myrvector)
       myrvector(:)=rvector(:)

!$OMP DO SCHEDULE(DYNAMIC)
       do config1=www%botconfig,www%topconfig
          if (www%doublehopdiagflag(config1).ne.0) then
             outsmallvector(:,config1)=mattrans(www%doublediaghop(config1),config1) * &
                  insmallvector(:,config1) * myrvector(:)
          endif
       enddo
!$OMP END DO
!$OMP END PARALLEL

    endif

  end subroutine arbitrary_sparsemult_doubles_byproc

  subroutine arbitraryconfig_mult_singles_byproc_hhh(firstproc,lastproc,www,in_onebodymat, &
       rvector, avectorin, avectorout,inrnum,diagflag,inholeflag)
    use walkmod
    use mpimod   !! myrank
    implicit none
    type(walktype),intent(in) :: www
    integer,intent(in) :: firstproc,lastproc,inholeflag,inrnum,diagflag
    DATATYPE,intent(in) :: in_onebodymat(www%nspf,www%nspf), &
         avectorin(inrnum,www%allbotconfigs(firstproc):www%alltopconfigs(lastproc))
    DATATYPE,intent(out) :: avectorout(inrnum,www%botconfig:www%topconfig)
    DATAECS,intent(in) :: rvector(inrnum)
    DATATYPE :: onebodymat(www%nspf,www%nspf)    !! AUTOMATIC
    DATATYPE :: csum, outsum(inrnum), holetrace
    DATAECS :: myrvector(inrnum)
    integer ::    config2, config1, iwalk, idiag,ihop,ispf

    if (www%topconfig-www%botconfig+1.eq.0) then
       return
    endif

    avectorout(:,:)=0d0

    if (www%alltopconfigs(lastproc).lt.www%allbotconfigs(firstproc)) then
       return
    endif

    onebodymat(:,:)=in_onebodymat(:,:)

!! HOLEFLAG

    if (www%holeflag.ne.0.and.inholeflag.eq.1) then
       holetrace=0
       do ispf=1,www%nspf
          holetrace=holetrace + onebodymat(ispf,ispf) * 2
       enddo

       onebodymat(:,:) = (-1) * in_onebodymat(:,:)

       if (myrank.ge.firstproc.and.myrank.le.lastproc) then
          do config1=www%botconfig,www%topconfig
             avectorout(:,config1) = holetrace*avectorin(:,config1)*rvector(:)
          enddo
       endif
    endif    !! holeflag

    if (diagflag.eq.0) then

!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(config1,ihop,iwalk,config2,csum,outsum,myrvector)
       myrvector(:)=rvector(:)

!$OMP DO SCHEDULE(DYNAMIC)
       do config1=www%botconfig,www%topconfig
          outsum(:)=avectorout(:,config1)
          do ihop=www%firstsinglehopbyproc(firstproc,config1),www%lastsinglehopbyproc(lastproc,config1)
             config2=www%singlehop(ihop,config1)
             csum=0d0
             do iwalk=www%singlehopwalkstart(ihop,config1),www%singlehopwalkend(ihop,config1)
                csum=csum +  onebodymat(www%singlewalkopspf(1,iwalk+www%scol(config1)), &
                     www%singlewalkopspf(2,iwalk+www%scol(config1))) *  &
                     www%singlewalkdirphase(iwalk+www%scol(config1))
             enddo
             outsum(:)=outsum(:)+avectorin(:,config2) * csum * myrvector(:)
          enddo
          avectorout(:,config1)=outsum(:)

       enddo
!$OMP END DO
!$OMP END PARALLEL

    else

!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(config1,idiag,iwalk,csum,myrvector)
       myrvector(:)=rvector(:)

!$OMP DO SCHEDULE(DYNAMIC)
       do config1=www%botconfig,www%topconfig
          csum=0d0
          do idiag=1,www%numsinglediagwalks(config1)
             iwalk=www%singlediag(idiag,config1)
             csum=csum+&
                  onebodymat(www%singlewalkopspf(1,iwalk+www%scol(config1)), &
                  www%singlewalkopspf(2,iwalk+www%scol(config1))) *  &
                  www%singlewalkdirphase(iwalk+www%scol(config1))
          enddo
          avectorout(:,config1)=avectorout(:,config1) + &
               avectorin(:,config1) * csum * myrvector(:)
       enddo
!$OMP END DO
!$OMP END PARALLEL

    endif

  end subroutine arbitraryconfig_mult_singles_byproc_hhh

  subroutine arbitraryconfig_mult_singles_byproc(firstproc,lastproc,www,onebodymat, &
       rvector, avectorin, avectorout,inrnum,diagflag)
    use walkmod
    use mpimod   !! myrank
    implicit none
    type(walktype),intent(in) :: www
    integer,intent(in) :: firstproc,lastproc,inrnum,diagflag
    DATATYPE,intent(in) :: onebodymat(www%nspf,www%nspf), &
         avectorin(inrnum,www%allbotconfigs(firstproc):www%alltopconfigs(lastproc))
    DATATYPE,intent(out) :: avectorout(inrnum,www%botconfig:www%topconfig)
    DATAECS,intent(in) :: rvector(inrnum)

    call arbitraryconfig_mult_singles_byproc_hhh(firstproc,lastproc,www,onebodymat, &
         rvector, avectorin, avectorout,inrnum,diagflag,1)

  end subroutine arbitraryconfig_mult_singles_byproc

  subroutine arbitraryconfig_mult_doubles_byproc(firstproc,lastproc,www,in_twobodymat, &
       rvector, avectorin, avectorout,inrnum,diagflag)
    use walkmod
    use mpimod   !! myrank
    implicit none
    type(walktype),intent(in) :: www
    integer,intent(in) :: firstproc,lastproc,inrnum,diagflag
    DATATYPE,intent(in) :: in_twobodymat(www%nspf,www%nspf,www%nspf,www%nspf), &
         avectorin(inrnum,www%allbotconfigs(firstproc):www%alltopconfigs(lastproc))
    DATATYPE,intent(out) :: avectorout(inrnum,www%botconfig:www%topconfig)
    DATAECS,intent(in) :: rvector(inrnum)
    DATATYPE, allocatable :: twobodymat(:,:,:,:)
    DATATYPE :: csum, outsum(inrnum),holetrace
    DATAECS :: myrvector(inrnum)
    integer ::   config2, config1, iwalk, idiag, ihop,ispf,jspf

    if (www%topconfig-www%botconfig+1.eq.0) then
       return
    endif

    avectorout(:,:)=0.d0

    if (www%alltopconfigs(lastproc).lt.www%allbotconfigs(firstproc)) then
       return
    endif

    allocate(twobodymat(www%nspf,www%nspf,www%nspf,www%nspf))
    twobodymat(:,:,:,:) = in_twobodymat(:,:,:,:)

    if (www%holeflag.ne.0) then

       holetrace=0
       do ispf=1,www%nspf
          do jspf=1,www%nspf
             holetrace=holetrace + twobodymat(ispf,ispf,jspf,jspf) * 2 &
                  - twobodymat(ispf,jspf,jspf,ispf)
          enddo
       enddo

!!       OFLWR "HTRACE002",holetrace,www%nspf; CFL

!! NO!  Two-electron operator       twobodymat(:,:,:,:) = (-1) * in_twobodymat(:,:,:,:)

       if (myrank.ge.firstproc.and.myrank.le.lastproc) then
          do config1=www%botconfig,www%topconfig
             avectorout(:,config1) = holetrace*avectorin(:,config1) * rvector(:)
          enddo
       endif
    endif    !! holeflag

    if (diagflag.eq.0) then

!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(config1,ihop,iwalk,config2,csum,outsum,myrvector)
       myrvector(:)=rvector(:)
!$OMP DO SCHEDULE(DYNAMIC)
       do config1=www%botconfig,www%topconfig
          outsum(:)=0d0
          do ihop=www%firstdoublehopbyproc(firstproc,config1),www%lastdoublehopbyproc(lastproc,config1)
             config2=www%doublehop(ihop,config1)
             csum=0d0
             do iwalk=www%doublehopwalkstart(ihop,config1),www%doublehopwalkend(ihop,config1)
                csum=csum+&
                     twobodymat(www%doublewalkdirspf(1,iwalk+www%dcol(config1)), &
                     www%doublewalkdirspf(2,iwalk+www%dcol(config1)),   &
                     www%doublewalkdirspf(3,iwalk+www%dcol(config1)),   &
                     www%doublewalkdirspf(4,iwalk+www%dcol(config1)))* &
                     www%doublewalkdirphase(iwalk+www%dcol(config1))
             enddo
             outsum(:)=outsum(:)+avectorin(:,config2) *  csum * myrvector(:)
          enddo
          avectorout(:,config1)=avectorout(:,config1)+outsum(:)
       enddo
!$OMP END DO
!$OMP END PARALLEL

    else

!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(config1,idiag,iwalk,csum,myrvector)
       myrvector(:)=rvector(:)
!$OMP DO SCHEDULE(DYNAMIC)
       do config1=www%botconfig,www%topconfig
          csum=0d0
          do idiag=1,www%numdoublediagwalks(config1)
             iwalk=www%doublediag(idiag,config1)
             csum=csum+&
                  twobodymat(www%doublewalkdirspf(1,iwalk+www%dcol(config1)), &
                  www%doublewalkdirspf(2,iwalk+www%dcol(config1)),   &
                  www%doublewalkdirspf(3,iwalk+www%dcol(config1)),   &
                  www%doublewalkdirspf(4,iwalk+www%dcol(config1)))* &
                  www%doublewalkdirphase(iwalk+www%dcol(config1))
          enddo
          avectorout(:,config1)=avectorout(:,config1) + &
               avectorin(:,config1) * csum * myrvector(:)
       enddo
!$OMP END DO
!$OMP END PARALLEL

    endif

    deallocate(twobodymat)

  end subroutine arbitraryconfig_mult_doubles_byproc

  subroutine arbitraryconfig_mult_singles_noparcon(www,onebodymat, rvector, avectorin, avectorout,inrnum)
    use walkmod
    use mpimod   !! myrank,nprocs
    use fileptrmod
    implicit none
    type(walktype),intent(in) :: www
    integer,intent(in) :: inrnum
    DATATYPE,intent(in) :: onebodymat(www%nspf,www%nspf), avectorin(inrnum,www%firstconfig:www%lastconfig)
    DATATYPE,intent(out) :: avectorout(inrnum,www%firstconfig:www%lastconfig)
    DATAECS,intent(in) :: rvector(inrnum)

    if (www%parconsplit.ne.0) then
       OFLWR "ERROR NOPARCON ARB"; CFLST
    endif   !! now botconfig:topconfig = firstconfig:lastconfig

    call arbitraryconfig_mult_singles_byproc(1,nprocs,www,onebodymat, &
         rvector, avectorin,avectorout(:,www%botconfig:www%topconfig),inrnum,0)

    call mpiallgather(avectorout,www%numconfig*inrnum,www%configsperproc(:)*inrnum,&
         www%maxconfigsperproc*inrnum)

  end subroutine arbitraryconfig_mult_singles_noparcon

#ifdef MPIFLAG

  subroutine arbitraryconfig_mult_singles_gather(www,onebodymat, rvector, avectorin, avectorout,inrnum)
    use walkmod
    use mpimod   !! myrank,nprocs
    use fileptrmod
    implicit none
    type(walktype),intent(in) :: www
    integer,intent(in) :: inrnum
    DATATYPE,intent(in) :: onebodymat(www%nspf,www%nspf), avectorin(inrnum,www%firstconfig:www%lastconfig)
    DATATYPE,intent(out) :: avectorout(inrnum,www%firstconfig:www%lastconfig)
    DATAECS,intent(in) :: rvector(inrnum)
    DATATYPE, allocatable :: workvector(:,:)

    if (www%parconsplit.eq.0) then
       OFLWR "ERROR PARCON ARB"; CFLST
    endif   !! now botconfig:topconfig = firstconfig:lastconfig

    allocate(workvector(inrnum,www%numconfig))
    workvector(:,:)=0d0

    if (www%topconfig.ge.www%botconfig) then
       workvector(:,www%botconfig:www%topconfig) = avectorin(:,:)
    endif
    call mpiallgather(workvector,www%numconfig*inrnum,www%configsperproc(:)*inrnum,&
         www%maxconfigsperproc*inrnum)

    if (www%lastconfig.ge.www%firstconfig) then
       call arbitraryconfig_mult_singles_byproc(1,nprocs,www,onebodymat,&
            rvector, workvector, avectorout,inrnum,0)
    endif
    deallocate(workvector)

  end subroutine arbitraryconfig_mult_singles_gather

  subroutine arbitraryconfig_mult_singles_summa(www,onebodymat, rvector, avectorin, avectorout,inrnum)
    use walkmod
    use mpimod   !! myrank,nprocs
    use fileptrmod
    implicit none
    type(walktype),intent(in) :: www
    integer,intent(in) :: inrnum
    DATATYPE,intent(in) :: onebodymat(www%nspf,www%nspf), &
         avectorin(inrnum,www%firstconfig:www%lastconfig)
    DATATYPE,intent(out) :: avectorout(inrnum,www%firstconfig:www%lastconfig)
    DATAECS,intent(in) :: rvector(inrnum)
    DATATYPE :: workvector(inrnum,www%maxconfigsperproc),&
         outtemp(inrnum,www%botconfig:www%topconfig+1)         !! AUTOMATIC
    integer :: iproc

    if (www%parconsplit.eq.0) then
       OFLWR "ERROR PARCON ARB"; CFLST
    endif   !! now botconfig:topconfig = firstconfig:lastconfig

    workvector=0; outtemp=0

    if (www%lastconfig.ge.www%firstconfig) then
       avectorout(:,:)=0d0
    endif

    do iproc=1,nprocs

       if (myrank.eq.iproc.and.www%lastconfig.ge.www%firstconfig) then
          workvector(:,1:www%topconfig-www%botconfig+1) = avectorin(:,:)
       endif
       if (www%alltopconfigs(iproc).ge.www%allbotconfigs(iproc)) then
          call mympibcast(workvector,iproc,&
               (www%alltopconfigs(iproc)-www%allbotconfigs(iproc)+1)*inrnum)
          if (www%lastconfig.ge.www%firstconfig) then
             call arbitraryconfig_mult_singles_byproc(iproc,iproc,www,onebodymat, &
                  rvector, workvector, outtemp,inrnum,0)
             avectorout(:,:) = avectorout(:,:) + outtemp(:,www%botconfig:www%topconfig)
          endif
       endif
    enddo

  end subroutine arbitraryconfig_mult_singles_summa

  subroutine arbitraryconfig_mult_singles_circ(www,onebodymat, rvector, avectorin, avectorout,inrnum)
    use walkmod
    use mpimod   !! myrank,nprocs
    use fileptrmod
    implicit none
    type(walktype),intent(in) :: www
    integer,intent(in) :: inrnum
    DATATYPE,intent(in) :: onebodymat(www%nspf,www%nspf),&
         avectorin(inrnum,www%firstconfig:www%lastconfig)
    DATATYPE,intent(out) :: avectorout(inrnum,www%firstconfig:www%lastconfig)
    DATAECS,intent(in) :: rvector(inrnum)
    DATATYPE :: workvector(inrnum,www%maxconfigsperproc),workvector2(inrnum,www%maxconfigsperproc),&
         outtemp(inrnum,www%botconfig:www%topconfig+1)         !! AUTOMATIC
    integer :: iproc,prevproc,nextproc,deltaproc

    if (www%parconsplit.eq.0) then
       OFLWR "ERROR PARCON ARB"; CFLST
    endif   !! now botconfig:topconfig = firstconfig:lastconfig

!! doing circ mult slightly different than e.g. SINCDVR/coreproject.f90
!!     and ftcore.f90, holding hands in a circle, prevproc and nextproc,
!!     each chunk gets passed around the circle
    prevproc=mod(nprocs+myrank-2,nprocs)+1
    nextproc=mod(myrank,nprocs)+1

    workvector=0; workvector2=0; outtemp=0

    if (www%topconfig.ge.www%botconfig) then
       avectorout=0d0
       workvector(:,1:www%topconfig-www%botconfig+1) = avectorin(:,:)
    endif

    do deltaproc=0,nprocs-1

!! PASSING BACKWARD (plus deltaproc)
       iproc=mod(myrank-1+deltaproc,nprocs)+1

       if (www%topconfig.ge.www%botconfig.and.www%alltopconfigs(iproc).ge.www%allbotconfigs(iproc)) then
          call arbitraryconfig_mult_singles_byproc(iproc,iproc,www,onebodymat, &
               rvector, workvector, outtemp,inrnum,0)
          avectorout(:,:) = avectorout(:,:) + outtemp(:,www%botconfig:www%topconfig)
       endif

!! PASSING BACKWARD
!! mympisendrecv(sendbuf,recvbuf,dest,source,...)

       call mympisendrecv(workvector,workvector2,prevproc,nextproc,deltaproc,&
            inrnum * www%maxconfigsperproc)
       workvector(:,:)=workvector2(:,:)

    enddo

  end subroutine arbitraryconfig_mult_singles_circ

#endif

  subroutine arbitraryconfig_mult_singles(www,onebodymat, rvector, avectorin, avectorout,inrnum)
    use walkmod
    use sparse_parameters
    use fileptrmod
    implicit none
    type(walktype),intent(in) :: www
    integer,intent(in) :: inrnum
    DATATYPE,intent(in) :: onebodymat(www%nspf,www%nspf), &
         avectorin(inrnum,www%firstconfig:www%lastconfig)
    DATATYPE,intent(out) :: avectorout(inrnum,www%firstconfig:www%lastconfig)
    DATAECS,intent(in) :: rvector(inrnum)

#ifdef MPIFLAG
    if (www%parconsplit.eq.0) then
#endif
       call arbitraryconfig_mult_singles_noparcon(www,onebodymat, rvector, avectorin, avectorout,inrnum)
#ifdef MPIFLAG
    else
       select case (sparsesummaflag)
       case (0)
          call arbitraryconfig_mult_singles_gather(www,onebodymat, rvector, avectorin, avectorout,inrnum)
       case (1)
          call arbitraryconfig_mult_singles_summa(www,onebodymat, rvector, avectorin, avectorout,inrnum)
       case (2)
          call arbitraryconfig_mult_singles_circ(www,onebodymat, rvector, avectorin, avectorout,inrnum)
       case default
          OFLWR "Error sparsesummaflag ", sparsesummaflag; CFLST
       end select
    endif
#endif

  end subroutine arbitraryconfig_mult_singles

end module arbitrarymultmod

module sparsemultmod
use arbitrarymultmod
contains

!! DIRECT MATVEC (with matrix_ptr not sparse_ptr)

  subroutine direct_sparseconfigmult_byproc(firstproc,lastproc,www,invector,outvector,matrix_ptr, &
       boflag, nucflag, pulseflag, conflag,time,onlytdflag,botr,topr,diagflag,imc)
    use ham_parameters
    use r_parameters
    use configptrmod
    use opmod   !! rkemod, proderivmod
    use mpimod  !! myrank
    use walkmod
    use pulsesubmod
    implicit none
    type(walktype),intent(in) :: www
    integer,intent(in) :: firstproc,lastproc,imc
    integer,intent(in) :: conflag,boflag,nucflag,pulseflag,onlytdflag,botr,topr,diagflag
    DATATYPE,intent(in) :: invector(botr:topr,www%allbotconfigs(firstproc):www%alltopconfigs(lastproc))
    DATATYPE,intent(out) :: outvector(botr:topr,www%botconfig:www%topconfig)
    Type(CONFIGPTR),intent(in) :: matrix_ptr
    real*8, intent(in) :: time
    integer :: mynumr,ir,flag,ii
    DATATYPE :: tempvector(botr:topr,www%botconfig:www%topconfig) !! AUTOMATIC
    DATATYPE :: tempmatel(www%nspf,www%nspf),facs(3),csum
    DATAECS :: rvector(botr:topr)
    real*8 :: gg

    if (www%topconfig.lt.www%botconfig) then
       return
    endif

    outvector(:,:)=0d0

    if (www%alltopconfigs(lastproc).lt.www%allbotconfigs(firstproc)) then
       return
    endif

    if (onlytdflag.ne.0) then
       facs(:)=0;     facs(onlytdflag)=1d0
    else
       call vectdpot(time,velflag,facs,imc)
    endif

    mynumr=topr-botr+1

    if (boflag==1) then

!! easy nuclear repulsion, hardwired like this for now

       if (myrank.ge.firstproc.and.myrank.le.lastproc) then
          do ir=botr,topr
             outvector(ir,:)=outvector(ir,:) + invector(ir,www%botconfig:www%topconfig) * ( &
                  energyshift + (nucrepulsion+frozenpotdiag)/bondpoints(ir) + &
                  frozenkediag/bondpoints(ir)**2 ) * matrix_ptr%constfac
          enddo
       endif

       rvector(:)=1/bondpoints(botr:topr)
       call arbitraryconfig_mult_doubles_byproc(firstproc,lastproc,www,matrix_ptr%xtwoematel(:,:,:,:),&
            rvector,invector,tempvector,mynumr,diagflag)
       outvector(:,:)=outvector(:,:)+tempvector(:,:)

       if (www%holeflag.ne.0) then
          call arbitraryconfig_mult_singles_byproc_hhh(firstproc,lastproc,www,matrix_ptr%xtwoe_htrace(:,:),&
               rvector,invector,tempvector,mynumr,diagflag,0)
          outvector(:,:)=outvector(:,:) + tempvector(:,:)
       endif

       call arbitraryconfig_mult_singles_byproc(firstproc,lastproc,www,matrix_ptr%xpotmatel(:,:),&
            rvector,invector,tempvector,mynumr,diagflag)
       outvector(:,:)=outvector(:,:)+tempvector(:,:)

       rvector(:)=1/bondpoints(botr:topr)**2
       call arbitraryconfig_mult_singles_byproc(firstproc,lastproc,www,matrix_ptr%xopmatel(:,:),&
            rvector,invector,tempvector,mynumr,diagflag)
       outvector(:,:)=outvector(:,:)+tempvector(:,:)

    endif

    if (conflag==1.and.constraintflag.ne.0) then
       rvector(:)=1

       tempmatel(:,:)=matrix_ptr%xconmatel(:,:)
      if (tdflag.ne.0) then
          tempmatel(:,:)=tempmatel(:,:)+matrix_ptr%xconmatelxx(:,:)*facs(1)
          tempmatel(:,:)=tempmatel(:,:)+matrix_ptr%xconmatelyy(:,:)*facs(2)
          tempmatel(:,:)=tempmatel(:,:)+matrix_ptr%xconmatelzz(:,:)*facs(3)
       endif
       call arbitraryconfig_mult_singles_byproc(firstproc,lastproc,www,tempmatel,&
            rvector,invector,tempvector,mynumr,diagflag)
       outvector(:,:)=outvector(:,:)-tempvector(:,:)

    endif
    if (pulseflag.eq.1.and.tdflag.ne.0.or.onlytdflag.ne.0) then
       if (velflag.eq.0) then
          rvector(:)=bondpoints(botr:topr)
       else
          rvector(:)=1/bondpoints(botr:topr)
       endif

       tempmatel(:,:)=0d0; flag=0
       if (onlytdflag.eq.0.or.onlytdflag.eq.1) then
          tempmatel(:,:)= tempmatel(:,:) + matrix_ptr%xpulsematelxx*facs(1); flag=1
       endif
       if (onlytdflag.eq.0.or.onlytdflag.eq.2) then
          tempmatel(:,:)= tempmatel(:,:) + matrix_ptr%xpulsematelyy*facs(2); flag=1
       endif
       if (onlytdflag.eq.0.or.onlytdflag.eq.3) then
          tempmatel(:,:)= tempmatel(:,:) + matrix_ptr%xpulsematelzz*facs(3); flag=1
       endif
       if (flag.eq.1) then
          call arbitraryconfig_mult_singles_byproc(firstproc,lastproc,www,tempmatel,&
               rvector,invector,tempvector,mynumr,diagflag)
          outvector(:,:)=outvector(:,:)+tempvector(:,:)
       endif

!! STILL NEED A-VECTOR TERM IN BOND LENGTH FOR HETERO!! TO DO!!

       gg=0.25d0

       if (myrank.ge.firstproc.and.myrank.le.lastproc) then

!!  A-squared term and nuclear dipole:  nuclear dipole always
          if (velflag.eq.0) then
             csum=0
             do ii=1,3
                if (onlytdflag.eq.0.or.onlytdflag.eq.ii) then
                   csum=csum+facs(ii)*matrix_ptr%xpulsenuc(ii)
                endif
             enddo
             do ir=botr,topr
                outvector(ir,:)=outvector(ir,:)+invector(ir,www%botconfig:www%topconfig) * &
                     csum*bondpoints(ir)
             enddo
          else if (onlytdflag.eq.0) then
!! a-squared
             csum=matrix_ptr%constfac * www%numelec * (facs(1)**2 + facs(2)**2 + facs(3)**2 ) * 2  * gg
             do ir=botr,topr
!! NO R FACTOR !!
                outvector(ir,:)=outvector(ir,:)+invector(ir,www%botconfig:www%topconfig)*csum
             enddo
          endif
       endif
    endif   !! PULSE

    if (nonuc_checkflag.eq.0.and.nucflag.eq.1.and.mynumr.eq.numr) then
       if (diagflag.eq.0) then
          rvector(:)=1d0
          call arbitraryconfig_mult_singles_byproc(firstproc,lastproc,www,matrix_ptr%xymatel,&
               rvector,invector,tempvector,numr,0)
          call MYGEMM('N','N',numr,www%topconfig-www%botconfig+1,numr,DATAONE, &
               proderivmod,numr,     tempvector,           numr,DATAONE,outvector,numr)
          if (myrank.ge.firstproc.and.myrank.le.lastproc) then
             call MYGEMM('N','N',numr,www%topconfig-www%botconfig+1,numr,DATAONE*matrix_ptr%kefac,&
                  rkemod,numr,     invector(:,www%botconfig),numr,DATAONE,outvector,numr)
          endif
!!$     else
!!$        do ir=1,numr
!!$           outvector(ir,:)=outvector(ir,:)+&
!!$                invector(ir,www%botconfig:www%topconfig)*matrix_ptr%kefac*rkemod(ir,ir)
!!$        enddo
       endif
    endif

  end subroutine direct_sparseconfigmult_byproc

!! SPARSE MATVEC (with sparse_ptr not matrix_ptr)

  subroutine sparsesparsemult_byproc(firstproc,lastproc,www,invector,outvector,sparse_ptr,&
       boflag,nucflag,pulseflag,conflag,time,onlytdflag,botr,topr,diagflag,imc)
    use sparse_parameters
    use ham_parameters
    use r_parameters
    use sparseptrmod
    use walkmod
    use opmod   !! rkemod, proderivmod
    use mpimod  !! myrank
    use fileptrmod
    use pulsesubmod
    implicit none
    type(walktype),intent(in) :: www
    integer,intent(in) :: firstproc,lastproc,imc
    integer,intent(in) :: boflag,nucflag,pulseflag,conflag,onlytdflag,botr,topr,diagflag
    DATATYPE,intent(in)  :: invector(botr:topr,www%allbotconfigs(firstproc):www%alltopconfigs(lastproc))
    DATATYPE,intent(out) :: outvector(botr:topr,www%botconfig:www%topconfig)
    Type(SPARSEPTR),intent(in) :: sparse_ptr
    real*8,intent(in) :: time
    DATATYPE :: facs(3),csum
    DATATYPE :: tempvector(botr:topr,www%botconfig:www%topconfig), &
         tempsparsemattr(www%singlematsize,www%botconfig:www%topconfig) !! AUTOMATIC
    real*8 :: gg
    DATAECS :: rvector(botr:topr)
    integer :: ir,mynumr,ii,flag

    if (sparseconfigflag.eq.0) then
       OFLWR "error no sparsesparsemult if sparseconfigflag.eq.0"; CFLST
    endif

    if (www%topconfig.lt.www%botconfig) then
       return
    endif

    if (botr.ne.1.or.topr.ne.numr) then
       OFLWR "REPROGRAMMMM."; CFLST
    endif

    outvector(:,:)=0d0

    if (www%alltopconfigs(lastproc).lt.www%allbotconfigs(firstproc)) then
       return
    endif

    if (onlytdflag.ne.0) then
       facs(:)=0;     facs(onlytdflag)=1d0
    else
       call vectdpot(time,velflag,facs,imc)
    endif

    mynumr=topr-botr+1

    if (boflag==1) then

       if (myrank.ge.firstproc.and.myrank.le.lastproc) then
          do ir=botr,topr
             outvector(ir,:)=outvector(ir,:) + invector(ir,www%botconfig:www%topconfig) * ( &
                  energyshift + (nucrepulsion+frozenpotdiag)/bondpoints(ir) + &
                  frozenkediag/bondpoints(ir)**2 ) * sparse_ptr%constfac
          enddo
       endif

       rvector(:)=1/bondpoints(botr:topr)
       call arbitrary_sparsemult_doubles_byproc(firstproc,lastproc,www,sparse_ptr%xpotsparsemattr,&
            rvector,invector,tempvector,mynumr,diagflag)
       outvector(:,:)=outvector(:,:)+tempvector(:,:)

       if (www%holeflag.ne.0) then
          call arbitrary_sparsemult_singles_byproc(firstproc,lastproc,www,sparse_ptr%xpottracemattr,&
               rvector,invector,tempvector,mynumr,diagflag)
          outvector(:,:)=outvector(:,:) + tempvector(:,:)
       endif

       call arbitrary_sparsemult_singles_byproc(firstproc,lastproc,www,sparse_ptr%xonepotsparsemattr,&
            rvector,invector,tempvector,mynumr,diagflag)
       outvector(:,:)=outvector(:,:)+tempvector(:,:)

       rvector(:)=1/bondpoints(botr:topr)**2
       call arbitrary_sparsemult_singles_byproc(firstproc,lastproc,www,sparse_ptr%xopsparsemattr,&
            rvector,invector,tempvector,mynumr,diagflag)
       outvector(:,:)=outvector(:,:)+tempvector(:,:)

    endif

    if (conflag==1.and.constraintflag.ne.0) then

       rvector(:)=1

       tempsparsemattr(:,:)=sparse_ptr%xconsparsemattr(:,:)
       if (tdflag.ne.0) then
          tempsparsemattr(:,:)=tempsparsemattr(:,:)+sparse_ptr%xconsparsemattrxx(:,:)*facs(1)
          tempsparsemattr(:,:)=tempsparsemattr(:,:)+sparse_ptr%xconsparsemattryy(:,:)*facs(2)
          tempsparsemattr(:,:)=tempsparsemattr(:,:)+sparse_ptr%xconsparsemattrzz(:,:)*facs(3)
       endif
       call arbitrary_sparsemult_singles_byproc(firstproc,lastproc,www,tempsparsemattr,&
            rvector,invector,tempvector,mynumr,diagflag)
       outvector(:,:)=outvector(:,:)-tempvector(:,:)

    endif

    if (pulseflag.eq.1.and.tdflag.ne.0.or.onlytdflag.ne.0) then
       if (velflag.eq.0) then
          rvector(:)=bondpoints(botr:topr)
       else
          rvector(:)=1/bondpoints(botr:topr)
       endif

       tempsparsemattr(:,:)=0d0; flag=0
       if (onlytdflag.eq.0.or.onlytdflag.eq.1) then
          tempsparsemattr= tempsparsemattr + sparse_ptr%xpulsesparsemattrxx*facs(1); flag=1
       endif
       if (onlytdflag.eq.0.or.onlytdflag.eq.2) then
          tempsparsemattr= tempsparsemattr + sparse_ptr%xpulsesparsemattryy*facs(2); flag=1
       endif
       if (onlytdflag.eq.0.or.onlytdflag.eq.3) then
          tempsparsemattr= tempsparsemattr + sparse_ptr%xpulsesparsemattrzz*facs(3); flag=1
       endif
       if (flag.eq.1) then
          call arbitrary_sparsemult_singles_byproc(firstproc,lastproc,www,tempsparsemattr,&
               rvector,invector,tempvector,mynumr,diagflag)
          outvector(:,:)=outvector(:,:)+tempvector(:,:)
       endif

!! STILL NEED A-VECTOR TERM IN BOND LENGTH FOR HETERO !! TO DO !!

!! A-squared term and nuclear dipole : nuclear dipole always

       gg=0.25d0

       if (myrank.ge.firstproc.and.myrank.le.lastproc) then

          if (velflag.eq.0) then
             csum=0
             do ii=1,3
                if (onlytdflag.eq.0.or.onlytdflag.eq.ii) then
                   csum=csum+facs(ii)*sparse_ptr%xpulsenuc(ii)
                endif
             enddo
             do ir=botr,topr
                outvector(ir,:)=outvector(ir,:)+invector(ir,www%botconfig:www%topconfig) * &
                     csum*bondpoints(ir)
             enddo
          else if (onlytdflag.eq.0) then
!! a-squared
             csum= sparse_ptr%constfac * www%numelec * (facs(1)**2 + facs(2)**2 + facs(3)**2 ) * 2 * gg
             do ir=botr,topr
!! NO R FACTOR !!
                outvector(ir,:)=outvector(ir,:)+invector(ir,www%botconfig:www%topconfig)*csum
             enddo
          endif
       endif
    endif !! PULSE

    if (nonuc_checkflag.eq.0.and.nucflag.eq.1.and.mynumr.eq.numr) then
       if (diagflag.eq.0) then
          rvector(:)=1d0
          call arbitrary_sparsemult_singles_byproc(firstproc,lastproc,www,sparse_ptr%xysparsemattr,&
               rvector,invector,tempvector,numr,0)
          call MYGEMM('N','N',numr,www%topconfig-www%botconfig+1,numr,DATAONE, &
               proderivmod,numr,     tempvector,           numr,DATAONE,outvector,numr)
          if (myrank.ge.firstproc.and.myrank.le.lastproc) then
             call MYGEMM('N','N',numr,www%topconfig-www%botconfig+1,numr,DATAONE*sparse_ptr%kefac,&
                  rkemod,numr,     invector(:,www%botconfig),numr,DATAONE,outvector,numr)
          endif
!!$     else
!!$        do ir=1,numr
!!$           outvector(ir,:)=outvector(ir,:)+&
!!$             invector(ir,www%botconfig:www%topconfig)*sparse_ptr%kefac*rkemod(ir,ir)
!!$        enddo
       endif
    endif

  end subroutine sparsesparsemult_byproc

  subroutine sparseconfigmult_byproc(firstproc,lastproc,www,invector,outvector,matrix_ptr,sparse_ptr,&
       boflag, nucflag, pulseflag, conflag,time,onlytdflag,botr,topr,diagflag,imc)
    use sparse_parameters
    use configptrmod
    use sparseptrmod
    use walkmod
    implicit none
    type(walktype),intent(in) :: www
    integer,intent(in) :: firstproc,lastproc
    integer,intent(in) :: conflag,boflag,nucflag,pulseflag,onlytdflag,botr,topr,diagflag,imc
    DATATYPE,intent(in) :: invector(botr:topr,www%allbotconfigs(firstproc):www%alltopconfigs(lastproc))
    DATATYPE,intent(out) :: outvector(botr:topr,www%botconfig:www%topconfig)
    real*8,intent(in) :: time
    Type(CONFIGPTR),intent(in) :: matrix_ptr
    Type(SPARSEPTR),intent(in) :: sparse_ptr

    if (www%topconfig-www%botconfig+1.eq.0) then
       return
    endif

    if (sparseopt.eq.0) then
       call direct_sparseconfigmult_byproc(firstproc,lastproc,www,invector,outvector,matrix_ptr, &
            boflag, nucflag, pulseflag, conflag,time,onlytdflag,botr,topr,diagflag,imc)
    else
       call sparsesparsemult_byproc(firstproc,lastproc,www,invector,outvector,sparse_ptr, &
            boflag, nucflag, pulseflag, conflag,time,onlytdflag,botr,topr,diagflag,imc)
    endif

  end subroutine sparseconfigmult_byproc

!! sparseconfigmultone:
!! For one bond length, born oppenheimer Hamiltonian

  subroutine sparseconfigmultone_noparcon(www,invector,outvector,matrix_ptr,sparse_ptr, &
       boflag,pulseflag,conflag,isplit,time,imc)
    use configptrmod
    use sparseptrmod
    use walkmod
    use mpimod    !! nprocs,myrank
    use fileptrmod
    implicit none
    type(walktype),intent(in) :: www
    integer,intent(in) :: conflag, boflag, pulseflag, isplit,imc
    DATATYPE,intent(in) :: invector(www%firstconfig:www%lastconfig)
    DATATYPE,intent(out) :: outvector(www%firstconfig:www%lastconfig)
    Type(CONFIGPTR),intent(in) :: matrix_ptr
    Type(SPARSEPTR),intent(in) :: sparse_ptr
    real*8,intent(in) :: time

    if (www%parconsplit.ne.0) then
       OFLWR "ERROR NOPARCON", www%parconsplit; CFLST
    endif

    if (www%topconfig.ge.www%botconfig) then
       call sparseconfigmult_byproc(1,nprocs,www,invector,outvector(www%botconfig:www%topconfig),&
            matrix_ptr,sparse_ptr, boflag, 0, pulseflag, conflag,time,0,isplit,isplit,0,imc)
    endif
    call mpiallgather(outvector,www%numconfig,www%configsperproc(:),www%maxconfigsperproc)

  end subroutine sparseconfigmultone_noparcon

#ifdef MPIFLAG

  subroutine sparseconfigmultone_gather(www,invector,outvector,matrix_ptr,sparse_ptr, &
       boflag,pulseflag,conflag,isplit,time,imc)
    use configptrmod
    use sparseptrmod
    use walkmod
    use mpimod    !! nprocs,myrank
    use fileptrmod
    implicit none
    type(walktype),intent(in) :: www
    integer,intent(in) :: conflag, boflag, pulseflag, isplit,imc
    DATATYPE,intent(in) :: invector(www%firstconfig:www%lastconfig)
    DATATYPE,intent(out) :: outvector(www%firstconfig:www%lastconfig)
    Type(CONFIGPTR),intent(in) :: matrix_ptr
    Type(SPARSEPTR),intent(in) :: sparse_ptr
    real*8,intent(in) :: time
    DATATYPE,allocatable :: workvector(:)

    if (www%parconsplit.eq.0) then
       OFLWR "ERROR PARCON"; CFLST
    endif    !! now botconfig:topconfig = firstconfig:lastconfig

    allocate(workvector(www%numconfig))
    workvector(:)=0d0

    if (www%topconfig.ge.www%botconfig) then
       workvector(www%botconfig:www%topconfig)=invector(:)
    endif
    call mpiallgather(workvector,www%numconfig,www%configsperproc(:),www%maxconfigsperproc)

    if (www%topconfig.ge.www%botconfig) then
       call sparseconfigmult_byproc(1,nprocs,www,workvector,outvector,&
          matrix_ptr,sparse_ptr, boflag, 0, pulseflag, conflag,time,0,isplit,isplit,0,imc)
    endif
    deallocate(workvector)

  end subroutine sparseconfigmultone_gather

  subroutine sparseconfigmultone_summa(www,invector,outvector,matrix_ptr,sparse_ptr, &
       boflag,pulseflag,conflag,isplit,time,imc)
    use configptrmod
    use sparseptrmod
    use walkmod
    use mpimod    !! nprocs,myrank
    use fileptrmod
    implicit none
    type(walktype),intent(in) :: www
    integer,intent(in) :: conflag, boflag, pulseflag, isplit,imc
    DATATYPE,intent(in) :: invector(www%firstconfig:www%lastconfig)
    DATATYPE,intent(out) :: outvector(www%firstconfig:www%lastconfig)
    Type(CONFIGPTR),intent(in) :: matrix_ptr
    Type(SPARSEPTR),intent(in) :: sparse_ptr
    real*8,intent(in) :: time
    DATATYPE :: workvector(www%maxconfigsperproc),outtemp(www%botconfig:www%topconfig+1)  !! AUTOMATIC
    integer :: iproc

    if (www%parconsplit.eq.0) then
       OFLWR "ERROR PARCON"; CFLST
    endif    !! now botconfig:topconfig = firstconfig:lastconfig

    workvector=0; outtemp=0

    if (www%lastconfig.ge.www%firstconfig) then
       outvector(:)=0d0
    endif

    do iproc=1,nprocs
       if (myrank.eq.iproc.and.www%lastconfig.ge.www%firstconfig) then
          workvector(1:www%topconfig-www%botconfig+1)=invector(:)
       endif
       if (www%alltopconfigs(iproc).ge.www%allbotconfigs(iproc)) then
          call mympibcast(workvector,iproc,(www%alltopconfigs(iproc)-www%allbotconfigs(iproc)+1))
          if (www%lastconfig.ge.www%firstconfig) then
             call sparseconfigmult_byproc(iproc,iproc,www,workvector,outtemp,matrix_ptr,sparse_ptr, &
                  boflag, 0, pulseflag, conflag,time,0,isplit,isplit,0,imc)
             outvector(:)=outvector(:)+outtemp(www%botconfig:www%topconfig)
          endif
       endif
    enddo

  end subroutine sparseconfigmultone_summa

  subroutine sparseconfigmultone_circ(www,invector,outvector,matrix_ptr,sparse_ptr, &
       boflag,pulseflag,conflag,isplit,time,imc)
    use configptrmod
    use sparseptrmod
    use walkmod
    use mpimod    !! nprocs,myrank
    use fileptrmod
    implicit none
    type(walktype),intent(in) :: www
    integer,intent(in) :: conflag, boflag, pulseflag, isplit,imc
    DATATYPE,intent(in) :: invector(www%firstconfig:www%lastconfig)
    DATATYPE,intent(out) :: outvector(www%firstconfig:www%lastconfig)
    Type(CONFIGPTR),intent(in) :: matrix_ptr
    Type(SPARSEPTR),intent(in) :: sparse_ptr
    real*8,intent(in) :: time
    DATATYPE :: workvector(www%maxconfigsperproc), workvector2(www%maxconfigsperproc),&
         outtemp(www%botconfig:www%topconfig+1)      !! AUTOMATIC
    integer :: iproc,prevproc,nextproc,deltaproc

    if (www%parconsplit.eq.0) then
       OFLWR "ERROR PARCON"; CFLST
    endif     !! now botconfig:topconfig = firstconfig:lastconfig

!! doing circ mult slightly different than e.g. SINCDVR/coreproject.f90
!!     and ftcore.f90, holding hands in a circle, prevproc and nextproc,
!!     each chunk gets passed around the circle
    prevproc=mod(nprocs+myrank-2,nprocs)+1
    nextproc=mod(myrank,nprocs)+1

    workvector=0; workvector2=0; outtemp=0

    if (www%topconfig.ge.www%botconfig) then
       outvector(:)=0d0
       workvector(1:www%topconfig-www%botconfig+1)=invector(:)
    endif

    do deltaproc=0,nprocs-1

!! PASSING BACKWARD (plus deltaproc)
       iproc=mod(myrank-1+deltaproc,nprocs)+1

       if (www%topconfig.ge.www%botconfig.and.www%alltopconfigs(iproc).ge.www%allbotconfigs(iproc)) then
          call sparseconfigmult_byproc(iproc,iproc,www,workvector,outtemp,matrix_ptr,sparse_ptr, &
               boflag, 0, pulseflag, conflag,time,0,isplit,isplit,0,imc)
          outvector(:)=outvector(:)+outtemp(www%botconfig:www%topconfig)
       endif

!! PASSING BACKWARD
!! mympisendrecv(sendbuf,recvbuf,dest,source,...)

       call mympisendrecv(workvector,workvector2,prevproc,nextproc,deltaproc,www%maxconfigsperproc)
       workvector(:)=workvector2(:)

    enddo

  end subroutine sparseconfigmultone_circ

#endif

  subroutine sparseconfigmultone(www,invector,outvector,matrix_ptr,sparse_ptr, &
       boflag,pulseflag,conflag,isplit,time,imc)
    use sparse_parameters
    use configptrmod
    use sparseptrmod
    use walkmod
    use fileptrmod
    implicit none
    type(walktype),intent(in) :: www
    integer,intent(in) :: conflag, boflag, pulseflag, isplit,imc
    DATATYPE,intent(in) :: invector(www%firstconfig:www%lastconfig)
    DATATYPE,intent(out) :: outvector(www%firstconfig:www%lastconfig)
    Type(CONFIGPTR),intent(in) :: matrix_ptr
    Type(SPARSEPTR),intent(in) :: sparse_ptr
    real*8,intent(in) :: time

#ifdef MPIFLAG
    if (www%parconsplit.eq.0) then
#endif
       call sparseconfigmultone_noparcon(www,invector,outvector,matrix_ptr,sparse_ptr, &
            boflag,pulseflag,conflag,isplit,time,imc)
#ifdef MPIFLAG
    else
       select case(sparsesummaflag)
       case(0)
          call sparseconfigmultone_gather(www,invector,outvector,matrix_ptr,sparse_ptr, &
               boflag,pulseflag,conflag,isplit,time,imc)
       case(1)
          call sparseconfigmultone_summa(www,invector,outvector,matrix_ptr,sparse_ptr, &
               boflag,pulseflag,conflag,isplit,time,imc)
       case(2)
          call sparseconfigmultone_circ(www,invector,outvector,matrix_ptr,sparse_ptr, &
               boflag,pulseflag,conflag,isplit,time,imc)
       case default
          OFLWR "Error sparsesummaflag ",sparsesummaflag; CFLST
       end select
    endif
#endif

  end subroutine sparseconfigmultone

  subroutine sparseconfigmultxxx_noparcon(www,invector,outvector,matrix_ptr,sparse_ptr, &
       boflag, nucflag, pulseflag, conflag,time,onlytdflag,imc)
    use r_parameters
    use configptrmod
    use sparseptrmod
    use mpimod    !! myrank,nprocs
    use walkmod
    use fileptrmod
    implicit none
    type(walktype),intent(in) :: www
    integer,intent(in) :: conflag,boflag,nucflag,pulseflag,onlytdflag,imc
    DATATYPE,intent(in) :: invector(numr,www%firstconfig:www%lastconfig)
    DATATYPE,intent(out) :: outvector(numr,www%firstconfig:www%lastconfig)
    Type(CONFIGPTR),intent(in) :: matrix_ptr
    Type(SPARSEPTR),intent(in) :: sparse_ptr
    real*8,intent(in) :: time

    if (www%parconsplit.ne.0) then
       OFLWR "ERROR NOPARCON XXX", www%parconsplit; CFLST
    endif

    if (www%topconfig.ge.www%botconfig) then
       call sparseconfigmult_byproc(1,nprocs,www,invector,outvector(:,www%botconfig:www%topconfig),&
            matrix_ptr,sparse_ptr, boflag, nucflag, pulseflag, conflag,time,onlytdflag,1,numr,0,imc)
    endif
    call mpiallgather(outvector,www%numconfig*numr,www%configsperproc(:)*numr,www%maxconfigsperproc*numr)

  end subroutine sparseconfigmultxxx_noparcon

#ifdef MPIFLAG

  subroutine sparseconfigmultxxx_gather(www,invector,outvector,matrix_ptr,sparse_ptr, &
       boflag, nucflag, pulseflag, conflag,time,onlytdflag,imc)
    use r_parameters
    use configptrmod
    use sparseptrmod
    use mpimod    !! myrank,nprocs
    use walkmod
    use fileptrmod
    implicit none
    type(walktype),intent(in) :: www
    integer,intent(in) :: conflag,boflag,nucflag,pulseflag,onlytdflag,imc
    DATATYPE,intent(in) :: invector(numr,www%firstconfig:www%lastconfig)
    DATATYPE,intent(out) :: outvector(numr,www%firstconfig:www%lastconfig)
    Type(CONFIGPTR),intent(in) :: matrix_ptr
    Type(SPARSEPTR),intent(in) :: sparse_ptr
    real*8,intent(in) :: time
    DATATYPE, allocatable :: workvector(:,:)

    if (www%parconsplit.eq.0) then
       OFLWR "ERROR PARCON XXX"; CFLST
    endif   !! now botconfig:topconfig = firstconfig:lastconfig

    allocate(workvector(numr,www%numconfig))
    workvector(:,:)=0d0

    if (www%topconfig.ge.www%botconfig) then
       workvector(:,www%botconfig:www%topconfig) = invector(:,:)
    endif
    call mpiallgather(workvector,www%numconfig*numr,www%configsperproc(:)*numr,&
         www%maxconfigsperproc*numr)

    if (www%topconfig.ge.www%botconfig) then
       call sparseconfigmult_byproc(1,nprocs,www,workvector,outvector,&
            matrix_ptr,sparse_ptr, boflag, nucflag, pulseflag, conflag,time,&
            onlytdflag,1,numr,0,imc)
    endif
    deallocate(workvector)

  end subroutine sparseconfigmultxxx_gather

  subroutine sparseconfigmultxxx_summa(www,invector,outvector,matrix_ptr,sparse_ptr, &
       boflag, nucflag, pulseflag, conflag,time,onlytdflag,imc)
    use r_parameters
    use configptrmod
    use sparseptrmod
    use mpimod    !! myrank,nprocs
    use walkmod
    use fileptrmod
    implicit none
    type(walktype),intent(in) :: www
    integer,intent(in) :: conflag,boflag,nucflag,pulseflag,onlytdflag,imc
    DATATYPE,intent(in) :: invector(numr,www%firstconfig:www%lastconfig)
    DATATYPE,intent(out) :: outvector(numr,www%firstconfig:www%lastconfig)
    Type(CONFIGPTR),intent(in) :: matrix_ptr
    Type(SPARSEPTR),intent(in) :: sparse_ptr
    real*8,intent(in) :: time
    DATATYPE :: workvector(numr,www%maxconfigsperproc), &
         outtemp(numr,www%botconfig:www%topconfig+1)    !!AUTOMATIC
    integer :: iproc

    if (www%parconsplit.eq.0) then
       OFLWR "ERROR PARCON XXX"; CFLST
    endif   !! now botconfig:topconfig = firstconfig:lastconfig

    workvector=0; outtemp=0

    if (www%lastconfig.ge.www%firstconfig) then
       outvector=0d0
    endif

    do iproc=1,nprocs

       if (myrank.eq.iproc.and.www%lastconfig.ge.www%firstconfig) then
          workvector(:,1:www%topconfig-www%botconfig+1) = invector(:,:)
       endif
       if (www%alltopconfigs(iproc).ge.www%allbotconfigs(iproc)) then
          call mympibcast(workvector,iproc,&
               (www%alltopconfigs(iproc)-www%allbotconfigs(iproc)+1)*numr)
          if (www%lastconfig.ge.www%firstconfig) then
             call sparseconfigmult_byproc(iproc,iproc,www,workvector,outtemp,matrix_ptr,sparse_ptr, &
                  boflag, nucflag, pulseflag, conflag,time,onlytdflag,1,numr,0,imc)
             outvector(:,:)=outvector(:,:) + outtemp(:,www%botconfig:www%topconfig)
          endif
       endif
    enddo

  end subroutine sparseconfigmultxxx_summa

  subroutine sparseconfigmultxxx_circ(www,invector,outvector,matrix_ptr,sparse_ptr, &
       boflag, nucflag, pulseflag, conflag,time,onlytdflag,imc)
    use r_parameters
    use configptrmod
    use sparseptrmod
    use mpimod    !! myrank,nprocs
    use walkmod
    use fileptrmod
    implicit none
    type(walktype),intent(in) :: www
    integer,intent(in) :: conflag,boflag,nucflag,pulseflag,onlytdflag,imc
    DATATYPE,intent(in) :: invector(numr,www%firstconfig:www%lastconfig)
    DATATYPE,intent(out) :: outvector(numr,www%firstconfig:www%lastconfig)
    Type(CONFIGPTR),intent(in) :: matrix_ptr
    Type(SPARSEPTR),intent(in) :: sparse_ptr
    real*8,intent(in) :: time
    DATATYPE :: workvector(numr,www%maxconfigsperproc), workvector2(numr,www%maxconfigsperproc),&
         outtemp(numr,www%botconfig:www%topconfig+1)      !! AUTOMATIC
    integer :: iproc,deltaproc,prevproc,nextproc

    if (www%parconsplit.eq.0) then
       OFLWR "ERROR PARCON XXX"; CFLST
    endif   !! now botconfig:topconfig = firstconfig:lastconfig

!! doing circ mult slightly different than e.g. SINCDVR/coreproject.f90
!!     and ftcore.f90, holding hands in a circle, prevproc and nextproc,
!!     each chunk gets passed around the circle
    prevproc=mod(nprocs+myrank-2,nprocs)+1
    nextproc=mod(myrank,nprocs)+1

    workvector=0; workvector2=0; outtemp=0

    if (www%topconfig.ge.www%botconfig) then
       outvector(:,:)=0d0
       workvector(:,1:www%topconfig-www%botconfig+1)=invector(:,:)
    endif

    do deltaproc=0,nprocs-1

!! PASSING BACKWARD (plus deltaproc)
       iproc=mod(myrank-1+deltaproc,nprocs)+1

       if (www%topconfig.ge.www%botconfig.and.www%alltopconfigs(iproc).ge.www%allbotconfigs(iproc)) then
          call sparseconfigmult_byproc(iproc,iproc,www,workvector,outtemp,matrix_ptr,sparse_ptr, &
               boflag, nucflag, pulseflag, conflag,time,onlytdflag,1,numr,0,imc)
          outvector(:,:)=outvector(:,:) + outtemp(:,www%botconfig:www%topconfig)
       endif

!! PASSING BACKWARD
!! mympisendrecv(sendbuf,recvbuf,dest,source,...)

       call mympisendrecv(workvector,workvector2,prevproc,nextproc,deltaproc,&
            numr*www%maxconfigsperproc)
       workvector(:,:)=workvector2(:,:)

    enddo

  end subroutine sparseconfigmultxxx_circ

#endif

  subroutine sparseconfigmultxxx(www,invector,outvector,matrix_ptr,sparse_ptr, &
       boflag, nucflag, pulseflag, conflag,time,onlytdflag,imc)
    use r_parameters
    use configptrmod
    use sparseptrmod
    use sparse_parameters
    use walkmod
    use fileptrmod
    implicit none
    type(walktype),intent(in) :: www
    integer,intent(in) :: conflag,boflag,nucflag,pulseflag,onlytdflag,imc
    DATATYPE,intent(in) :: invector(numr,www%firstconfig:www%lastconfig)
    DATATYPE,intent(out) :: outvector(numr,www%firstconfig:www%lastconfig)
    Type(CONFIGPTR),intent(in) :: matrix_ptr
    Type(SPARSEPTR),intent(in) :: sparse_ptr
    real*8,intent(in) :: time

#ifdef MPIFLAG
    if (www%parconsplit.eq.0) then
#endif
       call sparseconfigmultxxx_noparcon(www,invector,outvector,matrix_ptr,sparse_ptr, &
            boflag, nucflag, pulseflag, conflag,time,onlytdflag,imc)
#ifdef MPIFLAG
    else
       select case (sparsesummaflag)
       case(0)
          call sparseconfigmultxxx_gather(www,invector,outvector,matrix_ptr,sparse_ptr, &
               boflag, nucflag, pulseflag, conflag,time,onlytdflag,imc)
       case(1)
          call sparseconfigmultxxx_summa(www,invector,outvector,matrix_ptr,sparse_ptr, &
               boflag, nucflag, pulseflag, conflag,time,onlytdflag,imc)
       case(2)
          call sparseconfigmultxxx_circ(www,invector,outvector,matrix_ptr,sparse_ptr, &
               boflag, nucflag, pulseflag, conflag,time,onlytdflag,imc)
       case default
          OFLWR "Error sparsesummaflag ",sparsesummaflag; CFLST
       end select
    endif
#endif

  end subroutine sparseconfigmultxxx

!! All purpose subroutine.

  subroutine sparseconfigmult(www,invector,outvector,matrix_ptr,sparse_ptr, &
       boflag, nucflag, pulseflag, conflag,time,imc)
    use r_parameters
    use configptrmod
    use sparseptrmod
    use walkmod
    implicit none
    type(walktype),intent(in) :: www
    integer,intent(in) :: conflag,boflag,nucflag,pulseflag,imc
    DATATYPE,intent(in) :: invector(www%localnconfig*numr)
    DATATYPE,intent(out) :: outvector(www%localnconfig*numr)
    Type(CONFIGPTR),intent(in) :: matrix_ptr
    Type(SPARSEPTR),intent(in) :: sparse_ptr
    real*8,intent(in) :: time

    call sparseconfigmultxxx(www,invector,outvector,matrix_ptr,sparse_ptr,&
         boflag, nucflag, pulseflag, conflag,time,0,imc)

  end subroutine sparseconfigmult

!! All purpose subroutine, just the pulse.

  subroutine sparseconfigpulsemult(www,invector,outvector,matrix_ptr, sparse_ptr,which,imc)
    use r_parameters
    use configptrmod
    use sparseptrmod
    use walkmod
    implicit none
    type(walktype),intent(in) :: www
    integer,intent(in) :: which,imc
    DATATYPE,intent(in) :: invector(www%totadim)
    DATATYPE,intent(out) :: outvector(www%totadim)
    Type(CONFIGPTR),intent(in) :: matrix_ptr
    Type(SPARSEPTR),intent(in) :: sparse_ptr

    call sparseconfigmultxxx(www,invector,outvector,matrix_ptr,sparse_ptr, 0,0,1,0,-1d0,which,imc)

  end subroutine sparseconfigpulsemult

  subroutine parsparseconfigpreconmult(www,invector,outvector,matrix_ptr, sparse_ptr, &
       boflag, nucflag, pulseflag, conflag, inenergy,time)
    use r_parameters
    use configptrmod
    use sparseptrmod
    use walkmod
    use mpimod   !! myrank
    implicit none
    type(walktype),intent(in) :: www
    DATATYPE,intent(in) :: invector(numr,www%botconfig:www%topconfig)
    DATATYPE,intent(out) :: outvector(numr,www%botconfig:www%topconfig)
    DATATYPE,intent(in) :: inenergy
    Type(CONFIGPTR),intent(in) :: matrix_ptr
    Type(SPARSEPTR),intent(in) :: sparse_ptr
    integer,intent(in) ::  conflag,boflag,nucflag,pulseflag
    real*8,intent(in) :: time
    DATATYPE :: workvector(numr,www%botconfig:www%topconfig+1),&
         tempvector(numr,www%botconfig:www%topconfig+1)     !!AUTOMATIC

    if (www%topconfig-www%botconfig+1.eq.0) then
       return
    endif

    workvector(:,:)=1d0; tempvector(:,:)=0d0
    call sparseconfigmult_byproc(myrank,myrank,www,workvector,tempvector,matrix_ptr, sparse_ptr, &
         boflag, nucflag, pulseflag, conflag,time,0,1,numr,1,-1)

    outvector(:,:)=invector(:,:)/(tempvector(:,www%botconfig:www%topconfig)-inenergy)

  end subroutine parsparseconfigpreconmult

end module sparsemultmod


text:
SUBROUTINE TMSR

  ! **  SUBROUTINE TMSR WRITES TIME SERIES FILES FOR SURFACE ELEVATON,
  ! **  VELOCITY, CONCENTRATION, AND VOLUME SOURCES AT SPECIFIED
  ! **  (I,J) POINTS

  ! CHANGE RECORD

  USE GLOBAL

  IMPLICIT NONE

  INTEGER :: MLTM,MSDIG,MTMP,MFDIG,NX,NT,NTSSS,MTSSS,L,NMD,LMDCHHT
  INTEGER :: LMDCHUT,LMDCHVT,MTSCC,I,J,LN,K,KTMP,KTMP1,NSXD,NDOC

  REAL    :: TAUW1,TAUW2,TIME,SEDBTMP,SNDBTMP,PPTMP,TMPVAL,UTMP1,VTMP1,UTMP,VTMP
  REAL    :: TBEAST,TBNORT,TAUBDYN,TAUB2,CURANG,TAUTOT,VELMAG,TMPDRAG,TAUBSEDDYN
  REAL    :: TAUBSNDDYN,QBEDLOADX,QBEDLOADY,CQBEDLOADXT,CQBEDLOADYT,RUVTMP,QRNT

  CHARACTER*80 TITLE1,TITLE2,TITLE3,TITLE4,TITLE5,TITLE6,TITLE7, &
      TITLE11,TITLE12,TITLE13,TITLE14,TITLE15,TITLE16,TITLE17, &
      TITLE18,TITLE19,TFNTXWT,TFNTXWF,TFNTXWC,TFNTXWP, &
      TFNTXBT,TFNTXBF,TFNTXBC,TFNTXBP
  CHARACTER*10 CTUNIT
  CHARACTER*1 CZTT(0:9)
  CHARACTER*1 CCHTMF,CCHTMS
  CHARACTER*2 CNTOX(25),CNSND(10),CNSBL(10)

  REAL,ALLOCATABLE,DIMENSION(:) :: ATMPP
  REAL,ALLOCATABLE,DIMENSION(:) :: PORH
  REAL,ALLOCATABLE,DIMENSION(:) :: QCHANUIJ
  REAL,ALLOCATABLE,DIMENSION(:) :: QCHANVIJ
  REAL,ALLOCATABLE,DIMENSION(:) :: SEDSND
  REAL,ALLOCATABLE,DIMENSION(:) :: SEDSNDB
  REAL,ALLOCATABLE,DIMENSION(:) :: TXBC
  REAL,ALLOCATABLE,DIMENSION(:) :: TXBF
  REAL,ALLOCATABLE,DIMENSION(:) :: TXBP
  REAL,ALLOCATABLE,DIMENSION(:) :: TXBT
  REAL,ALLOCATABLE,DIMENSION(:) :: TXWC
  REAL,ALLOCATABLE,DIMENSION(:) :: TXWF
  REAL,ALLOCATABLE,DIMENSION(:) :: TXWP

  CHARACTER* 2,SAVE,ALLOCATABLE,DIMENSION(:)     :: CNTMSR
  CHARACTER*20,SAVE,ALLOCATABLE,DIMENSION(:)     :: FNAVB
  CHARACTER*20,SAVE,ALLOCATABLE,DIMENSION(:)     :: FNAVV
  CHARACTER*20,SAVE,ALLOCATABLE,DIMENSION(:)     :: FNBED
  CHARACTER*20,SAVE,ALLOCATABLE,DIMENSION(:)     :: FNDOX
  CHARACTER*20,SAVE,ALLOCATABLE,DIMENSION(:)     :: FNDYE
  CHARACTER*20,SAVE,ALLOCATABLE,DIMENSION(:)     :: FNNHX
  CHARACTER*20,SAVE,ALLOCATABLE,DIMENSION(:)     :: FNQ3D
  CHARACTER*20,SAVE,ALLOCATABLE,DIMENSION(:)     :: FNQQE
  CHARACTER*20,SAVE,ALLOCATABLE,DIMENSION(:)     :: FNSAL
  CHARACTER*22,SAVE,ALLOCATABLE,DIMENSION(:,:)   :: FNSBL
  CHARACTER*20,SAVE,ALLOCATABLE,DIMENSION(:)     :: FNSED
  CHARACTER*20,SAVE,ALLOCATABLE,DIMENSION(:)     :: FNSEL
  CHARACTER*20,SAVE,ALLOCATABLE,DIMENSION(:)     :: FNSFL
  CHARACTER*22,SAVE,ALLOCATABLE,DIMENSION(:,:)   :: FNSND
  CHARACTER*20,SAVE,ALLOCATABLE,DIMENSION(:)     :: FNTEM
  CHARACTER*20,SAVE,ALLOCATABLE,DIMENSION(:)     :: FNTOC
  CHARACTER*22,SAVE,ALLOCATABLE,DIMENSION(:,:)   :: FNTOX
  CHARACTER*23,SAVE,ALLOCATABLE,DIMENSION(:,:)   :: FNTXBC
  CHARACTER*23,SAVE,ALLOCATABLE,DIMENSION(:,:)   :: FNTXBF
  CHARACTER*23,SAVE,ALLOCATABLE,DIMENSION(:,:)   :: FNTXBP
  CHARACTER*23,SAVE,ALLOCATABLE,DIMENSION(:,:)   :: FNTXBT
  CHARACTER*23,SAVE,ALLOCATABLE,DIMENSION(:,:)   :: FNTXWC
  CHARACTER*23,SAVE,ALLOCATABLE,DIMENSION(:,:)   :: FNTXWF
  CHARACTER*23,SAVE,ALLOCATABLE,DIMENSION(:,:)   :: FNTXWP
  CHARACTER*23,SAVE,ALLOCATABLE,DIMENSION(:,:)   :: FNTXWT
  CHARACTER*20,SAVE,ALLOCATABLE,DIMENSION(:)     :: FNU3D
  CHARACTER*20,SAVE,ALLOCATABLE,DIMENSION(:)     :: FNUVE
  CHARACTER*20,SAVE,ALLOCATABLE,DIMENSION(:)     :: FNUVT
  CHARACTER*20,SAVE,ALLOCATABLE,DIMENSION(:)     :: FNV3D

  ALLOCATE(ATMPP(KCM))
  ALLOCATE(PORH(KBM))
  ALLOCATE(QCHANUIJ(LCM))
  ALLOCATE(QCHANVIJ(LCM))
  ALLOCATE(SEDSND(KCM))
  ALLOCATE(SEDSNDB(KBM))
  ALLOCATE(TXBC(KBM))
  ALLOCATE(TXBF(KBM))
  ALLOCATE(TXBP(KBM))
  ALLOCATE(TXBT(KBM))
  ALLOCATE(TXWC(KCM))
  ALLOCATE(TXWF(KCM))
  ALLOCATE(TXWP(KCM))

  ! *** ALLOCATE LOCAL ARRAYS
  IF(  .NOT. ALLOCATED(CNTMSR) )THEN
    ALLOCATE(CNTMSR(MLTMSRM))
    ALLOCATE(FNAVB(MLTMSRM))
    ALLOCATE(FNAVV(MLTMSRM))
    ALLOCATE(FNBED(MLTMSRM))
    ALLOCATE(FNDOX(MLTMSRM))
    ALLOCATE(FNDYE(MLTMSRM))
    ALLOCATE(FNNHX(MLTMSRM))
    ALLOCATE(FNQ3D(MLTMSRM))
    ALLOCATE(FNQQE(MLTMSRM))
    ALLOCATE(FNSAL(MLTMSRM))
    ALLOCATE(FNSBL(MLTMSRM,NSNM))
    ALLOCATE(FNSED(MLTMSRM))
    ALLOCATE(FNSEL(MLTMSRM))
    ALLOCATE(FNSFL(MLTMSRM))
    ALLOCATE(FNSND(MLTMSRM,NSNM))
    ALLOCATE(FNTEM(MLTMSRM))
    ALLOCATE(FNTOC(MLTMSRM))
    ALLOCATE(FNTOX(MLTMSRM,NTXM))
    ALLOCATE(FNTXBC(MLTMSRM,NTXM))
    ALLOCATE(FNTXBF(MLTMSRM,NTXM))
    ALLOCATE(FNTXBP(MLTMSRM,NTXM))
    ALLOCATE(FNTXBT(MLTMSRM,NTXM))
    ALLOCATE(FNTXWC(MLTMSRM,NTXM))
    ALLOCATE(FNTXWF(MLTMSRM,NTXM))
    ALLOCATE(FNTXWP(MLTMSRM,NTXM))
    ALLOCATE(FNTXWT(MLTMSRM,NTXM))
    ALLOCATE(FNU3D(MLTMSRM))
    ALLOCATE(FNUVE(MLTMSRM))
    ALLOCATE(FNUVT(MLTMSRM))
    ALLOCATE(FNV3D(MLTMSRM))

    CNTMSR = ' '
    FNAVB = ' '
    FNAVV = ' '
    FNBED = ' '
    FNDOX = ' '
    FNDYE = ' '
    FNNHX = ' '
    FNQ3D = ' '
    FNQQE = ' '
    FNSAL = ' '
    FNSBL = ' '
    FNSED = ' '
    FNSEL = ' '
    FNSFL = ' '
    FNSND = ' '
    FNTEM = ' '
    FNTOC = ' '
    FNU3D = ' '
    FNUVE = ' '
    FNUVT = ' '
    FNV3D = ' '
  ENDIF

  ! *** INITIALIZE LOCAL ARRAYS
  PORH=0.
  QCHANUIJ=0.
  QCHANVIJ=0.
  SEDSND=0.
  SEDSNDB=0.
  TXBC=0.
  TXBF=0.
  TXBP=0.
  TXBT=0.
  TXWC=0.
  TXWF=0.
  TXWP=0.

  IF( JSTMSR /= 1) GOTO 300
  !
  IF( MLTMSR > MLTMSRM )THEN
    WRITE (6,600)
    STOP
  ENDIF
  IF( MLTMSR > 99 )THEN
    WRITE (6,601)
    STOP
  ENDIF
  !
    600 FORMAT(' NUMBER OF TIME SER LOC, MLTMSR, EXCEEDS DIM, MLTMSRM')
    601 FORMAT(' NUMBER OF TIME SERIES LOCATIONS EXCEED 99')
  !
        TAUW1=0.
        TAUW2=0.
  !
  CZTT(0)='0'
  CZTT(1)='1'
  CZTT(2)='2'
  CZTT(3)='3'
  CZTT(4)='4'
  CZTT(5)='5'
  CZTT(6)='6'
  CZTT(7)='7'
  CZTT(8)='8'
  CZTT(9)='9'
  !
  DO MLTM=1,MLTMSR
  MSDIG=MOD(MLTM,10)
  MTMP=MLTM-MSDIG
  MFDIG=MTMP/10
  CCHTMF=CZTT(MFDIG)
  CCHTMS=CZTT(MSDIG)
  CNTMSR(MLTM)= CCHTMF // CCHTMS
  ENDDO
  !
  IF( TCTMSR == 1.) CTUNIT='SECONDS'
  IF( TCTMSR == 60.) CTUNIT='MINUTES'
  IF( TCTMSR == 3600.) CTUNIT='HOURS'
  IF( TCTMSR == 86400.) CTUNIT='DAYS'
  !
   CNTOX( 1)= '01'
   CNTOX( 2)= '02'
   CNTOX( 3)= '03'
   CNTOX( 4)= '04'
   CNTOX( 5)= '05'
   CNTOX( 6)= '06'
   CNTOX( 7)= '07'
   CNTOX( 8)= '08'
   CNTOX( 9)= '09'
   CNTOX(10)= '10'
   CNTOX(11)= '11'
   CNTOX(12)= '12'
   CNTOX(13)= '13'
   CNTOX(14)= '14'
   CNTOX(15)= '15'
   CNTOX(16)= '16'
   CNTOX(17)= '17'
   CNTOX(18)= '18'
   CNTOX(19)= '19'
   CNTOX(20)= '20'
   CNTOX(21)= '21'
   CNTOX(22)= '22'
   CNTOX(23)= '23'
   CNTOX(24)= '24'
   CNTOX(25)= '25'
  !
   CNSND( 1)= '01'
   CNSND( 2)= '02'
   CNSND( 3)= '03'
   CNSND( 4)= '04'
   CNSND( 5)= '05'
   CNSND( 6)= '06'
   CNSND( 7)= '07'
   CNSND( 8)= '08'
   CNSND( 9)= '09'
   CNSND(10)= '10'
  !
   CNSBL( 1)= '01'
   CNSBL( 2)= '02'
   CNSBL( 3)= '03'
   CNSBL( 4)= '04'
   CNSBL( 5)= '05'
   CNSBL( 6)= '06'
  !

  ! **  WRITE HEADINGS
  !
  TITLE1=' SALINITY (PSU) TIME SERIES, K=1,KC'
  TITLE2=' TEMPERATURE (DEG C) TIME SERIES, K=1,KC'
  TITLE3=' DYE CONC (KG/M**3) TIME SERIES, K=1,KC'
  TITLE4=' SED CONC (MG/LITER) TIME SERIES, K=1,KC'
  TITLE5=' TOXIC CONC (M/TOT VOL - UG/LITER) 1-4 BED,5-8 WC'
  TITLE6=' VISCOSITY (CM**2/S) TIME SERIES, K=1,KS'
  TITLE7=' DIFFUSIVITY (CM**2/S) TIME SERIES, K=1,KS'
  TITLE11=' SURFACE ELEVATION & DEPTH (METERS) TIME SERIES'
  TITLE12=' EXT MODE E,N VEL (CM/S) TBX TBY TB TBCG TBNG (CM/S)**2'
  TITLE13=' EXT MODE U,V TRANSPORT (M**3/S) TIME SERIES'
  TITLE14=' INT MODE EAST VEL (CM/S) TIME SERIES, K=1,KC'
  TITLE15=' INT MODE NORTH VEL (CM/S) TIME SERIES, K=1,KC'
  TITLE16=' EXT MODE VOLUME S/S (M**3/S) TIME SERIES'
  TITLE17=' INT MODE VOL S/S (M**3/S) TIME SERIES, K=1,KC'
  TITLE18=' SED BED LOAD QSX QSY (GM/S) CQSX CQSY (MG/L) '
  TITLE19=' BED TOP  KBT HBED(KBT) HBED(KBT-1) VOIDR FRAC SED/SND'
  TFNTXWT=' TOTAL TOXIC CONC WATER COL (M/TOT VOL), UG/LITER'
  TFNTXWF=' FREE DIS TOXIC CONC WATER COL (M/TOT VOL), UG/LITER'
  TFNTXWC=' DOC COMP TOXIC CONC WATER COL (M/TOT VOL), UG/LITER'
  TFNTXWP=' TOT PART TOXIC CONC WATER COL (M/M), UG/GM'
  TFNTXBT=' TOTAL TOXIC CONC SED BED (M/TOT VOL), UG/LITER'
  TFNTXBF=' FREE DIS TOXIC CONC SED BED (M/PORE VOL), UG/LITER'
  TFNTXBC=' DOC COMP TOXIC CONC SED BED (M/PORE VOL), UG/LITER'
  TFNTXBP=' TOT PART TOXIC CONC SED BED (M/M), UG/GM'
  !
  IF( ISTMSR == 2 )THEN
  DO MLTM=1,MLTMSR
    IF( MTMSRC(MLTM) == 1 )THEN
      IF( ISTRAN(1) >= 1 )THEN
        FNSAL(MLTM)=OUTDIR//'SALTS' // CNTMSR(MLTM) // '.OUT'
      ENDIF
      IF( ISTRAN(2) >= 1 )THEN
        FNTEM(MLTM)=OUTDIR//'TEMTS' // CNTMSR(MLTM) // '.OUT'
      ENDIF
      IF( ISTRAN(3) >= 1 )THEN
        FNDYE(MLTM)=OUTDIR//'DYETS' // CNTMSR(MLTM) // '.OUT'
      ENDIF
      IF( ISTRAN(4) >= 1 )THEN
        FNSFL(MLTM)=OUTDIR//'SFLTS' // CNTMSR(MLTM) // '.OUT'
      ENDIF
      IF( ISTRAN(6) >= 1 )THEN
        FNSED(MLTM)=OUTDIR//'SEDTS' // CNTMSR(MLTM) // '.OUT'
      ENDIF
  !          IF( ISTRAN(7) >= 1 )THEN
  !            FNSND(MLTM)=OUTDIR//'SNDTS' // CNTMSR(MLTM) // '.OUT'
  !          ENDIF
      IF( ISTRAN(7) >= 1 )THEN
        DO NX=1,NSND
        FNSND(MLTM,NX)=OUTDIR//'SND' // CNSND(NX) // 'TS' // CNTMSR(MLTM) // '.OUT'
        FNSBL(MLTM,NX)=OUTDIR//'SBL' // CNSBL(NX) // 'TS' // CNTMSR(MLTM) // '.OUT'
        ENDDO
      ENDIF
      IF( ISTRAN(8) >= 1 )THEN
        FNDOX(MLTM)=OUTDIR//'DOXTS' // CNTMSR(MLTM) // '.OUT'
        FNTOC(MLTM)=OUTDIR//'TOCTS' // CNTMSR(MLTM) // '.OUT'
        FNNHX(MLTM)=OUTDIR//'NHXTS' // CNTMSR(MLTM) // '.OUT'
      ENDIF
      IF( ISTRAN(5) >= 1 )THEN
        DO NT=1,NTOX
        FNTXWT(MLTM,NT)=OUTDIR//'TXWT' // CNTOX(NT) // 'TS' // CNTMSR(MLTM) // '.OUT'
        FNTXWF(MLTM,NT)=OUTDIR//'TXWF' // CNTOX(NT) // 'TS' // CNTMSR(MLTM) // '.OUT'
        FNTXWC(MLTM,NT)=OUTDIR//'TXWC' // CNTOX(NT) // 'TS' // CNTMSR(MLTM) // '.OUT'
        FNTXWP(MLTM,NT)=OUTDIR//'TXWP' // CNTOX(NT) // 'TS' // CNTMSR(MLTM) // '.OUT'
        FNTXBT(MLTM,NT)=OUTDIR//'TXBT' // CNTOX(NT) // 'TS' // CNTMSR(MLTM) // '.OUT'
        FNTXBF(MLTM,NT)=OUTDIR//'TXBF' // CNTOX(NT) // 'TS' // CNTMSR(MLTM) // '.OUT'
        FNTXBC(MLTM,NT)=OUTDIR//'TXBC' // CNTOX(NT) // 'TS' // CNTMSR(MLTM) // '.OUT'
        FNTXBP(MLTM,NT)=OUTDIR//'TXBP' // CNTOX(NT) // 'TS' // CNTMSR(MLTM) // '.OUT'
        ENDDO
      ENDIF
    ENDIF
    IF( MTMSRA(MLTM) == 1 )THEN
      FNAVV(MLTM)=OUTDIR//'AVVTS' // CNTMSR(MLTM) // '.OUT'
      FNAVB(MLTM)=OUTDIR//'AVBTS' // CNTMSR(MLTM) // '.OUT'
    ENDIF
    IF( MTMSRP(MLTM) == 1 )THEN
      FNSEL(MLTM)=OUTDIR//'SELTS' // CNTMSR(MLTM) // '.OUT'
    ENDIF
    IF( MTMSRUE(MLTM) == 1 )THEN
      FNUVE(MLTM)=OUTDIR//'UVETS' // CNTMSR(MLTM) // '.OUT'
    ENDIF
    IF( MTMSRUT(MLTM) == 1 )THEN
      FNUVT(MLTM)=OUTDIR//'UVTTS' // CNTMSR(MLTM) // '.OUT'
    ENDIF
    IF( MTMSRU(MLTM) == 1 )THEN
      FNU3D(MLTM)=OUTDIR//'U3DTS' // CNTMSR(MLTM) // '.OUT'
      FNV3D(MLTM)=OUTDIR//'V3DTS' // CNTMSR(MLTM) // '.OUT'
    ENDIF
    IF( MTMSRQE(MLTM) == 1 )THEN
      FNQQE(MLTM)=OUTDIR//'QQETS' // CNTMSR(MLTM) // '.OUT'
    ENDIF
    IF( MTMSRQ(MLTM) == 1 )THEN
      FNQ3D(MLTM)=OUTDIR//'Q3DTS' // CNTMSR(MLTM) // '.OUT'
    ENDIF
  ENDDO
  JSTMSR=0
  ENDIF
  !
  IF( JSTMSR == 0 ) GOTO 300
  !
  DO MLTM=1,MLTMSR
    IF( MTMSRC(MLTM) == 1 )THEN
      IF( ISTRAN(1) >= 1 )THEN
        FNSAL(MLTM)=OUTDIR//'SALTS' // CNTMSR(MLTM) // '.OUT'
        OPEN(11,FILE=FNSAL(MLTM),STATUS='UNKNOWN')
        CLOSE(11,STATUS='DELETE')
        OPEN(11,FILE=FNSAL(MLTM),STATUS='UNKNOWN')
        WRITE (11,100) TITLE1
        WRITE (11,101) CLTMSR(MLTM)
        WRITE (11,103)ILTMSR(MLTM),JLTMSR(MLTM)
        WRITE (11,102) CTUNIT
        CLOSE(11)
      ENDIF
      IF( ISTRAN(2) >= 1 )THEN
        FNTEM(MLTM)=OUTDIR//'TEMTS' // CNTMSR(MLTM) // '.OUT'
        OPEN(21,FILE=FNTEM(MLTM),STATUS='UNKNOWN')
        CLOSE(21,STATUS='DELETE')
        OPEN(21,FILE=FNTEM(MLTM),STATUS='UNKNOWN')
        WRITE (21,100) TITLE2
        WRITE (21,101) CLTMSR(MLTM)
        WRITE (21,103)ILTMSR(MLTM),JLTMSR(MLTM)
        WRITE (21,102) CTUNIT
        CLOSE(21)
      ENDIF
      IF( ISTRAN(3) >= 1 )THEN
        FNDYE(MLTM)=OUTDIR//'DYETS' // CNTMSR(MLTM) // '.OUT'
        OPEN(31,FILE=FNDYE(MLTM),STATUS='UNKNOWN')
        CLOSE(31,STATUS='DELETE')
        OPEN(31,FILE=FNDYE(MLTM),STATUS='UNKNOWN')
        WRITE (31,100) TITLE3
        WRITE (31,101) CLTMSR(MLTM)
        WRITE (31,103)ILTMSR(MLTM),JLTMSR(MLTM)
        WRITE (31,102) CTUNIT
        CLOSE(31)
      ENDIF
      IF( ISTRAN(4) >= 1 )THEN
        FNDYE(MLTM)=OUTDIR//'SFLTS' // CNTMSR(MLTM) // '.OUT'
        OPEN(31,FILE=FNSFL(MLTM),STATUS='UNKNOWN')
        CLOSE(31,STATUS='DELETE')
        OPEN(31,FILE=FNSFL(MLTM),STATUS='UNKNOWN')
        WRITE (31,100) TITLE3
        WRITE (31,101) CLTMSR(MLTM)
        WRITE (31,103)ILTMSR(MLTM),JLTMSR(MLTM)
        WRITE (31,102) CTUNIT
        CLOSE(31)
      ENDIF
      IF( ISTRAN(6) >= 1 )THEN
        FNSED(MLTM)=OUTDIR//'SEDTS' // CNTMSR(MLTM) // '.OUT'
        OPEN(41,FILE=FNSED(MLTM),STATUS='UNKNOWN')
        CLOSE(41,STATUS='DELETE')
        OPEN(41,FILE=FNSED(MLTM),STATUS='UNKNOWN')
        WRITE (41,100) TITLE4
        WRITE (41,101) CLTMSR(MLTM)
        WRITE (41,103)ILTMSR(MLTM),JLTMSR(MLTM)
        WRITE (41,102) CTUNIT
        CLOSE(41)
      ENDIF
      IF( ISTRAN(7) >= 1 )THEN
        DO NX=1,NSND
        FNSND(MLTM,NX)=OUTDIR//'SND'// CNSND(NX) // 'TS' // CNTMSR(MLTM) // '.OUT'
        OPEN(41,FILE=FNSND(MLTM,NX),STATUS='UNKNOWN')
        CLOSE(41,STATUS='DELETE')
        OPEN(41,FILE=FNSND(MLTM,NX),STATUS='UNKNOWN')
        WRITE (41,100) TITLE4
        WRITE (41,101) CLTMSR(MLTM)
        WRITE (41,103)ILTMSR(MLTM),JLTMSR(MLTM)
        WRITE (41,102) CTUNIT
        CLOSE(41)
        ENDDO
        DO NX=1,NSND
        FNSBL(MLTM,NX)=OUTDIR//'SBL'// CNSBL(NX) // 'TS' // CNTMSR(MLTM) // '.OUT'
        OPEN(41,FILE=FNSBL(MLTM,NX),STATUS='UNKNOWN')
        CLOSE(41,STATUS='DELETE')
        OPEN(41,FILE=FNSBL(MLTM,NX),STATUS='UNKNOWN')
        WRITE (41,100) TITLE18
        WRITE (41,101) CLTMSR(MLTM)
        WRITE (41,103)ILTMSR(MLTM),JLTMSR(MLTM)
        WRITE (41,102) CTUNIT
        CLOSE(41)
        ENDDO
  !            FNSND(MLTM)=OUTDIR//'SNDTS' // CNTMSR(MLTM) // '.OUT'
  !            OPEN(41,FILE=FNSND(MLTM),STATUS='UNKNOWN')
  !            CLOSE(41,STATUS='DELETE')
  !            OPEN(41,FILE=FNSND(MLTM),STATUS='UNKNOWN')
  !            WRITE (41,100) TITLE4
  !            WRITE (41,101) CLTMSR(MLTM)
  !            WRITE (41,103)ILTMSR(MLTM),JLTMSR(MLTM)
  !            WRITE (41,102) CTUNIT
  !            CLOSE(41)
      ENDIF
      IF( ISTRAN(6) >= 1 .OR. ISTRAN(7) >= 1 )THEN
        FNBED(MLTM)=OUTDIR//'BEDTS' // CNTMSR(MLTM) // '.OUT'
        OPEN(41,FILE=FNBED(MLTM),STATUS='UNKNOWN')
        CLOSE(41,STATUS='DELETE')
        OPEN(41,FILE=FNBED(MLTM),STATUS='UNKNOWN')
        WRITE (41,100) TITLE19
        WRITE (41,101) CLTMSR(MLTM)
        WRITE (41,103)ILTMSR(MLTM),JLTMSR(MLTM)
        WRITE (41,102) CTUNIT
        CLOSE(41)
      ENDIF
      IF( ISTRAN(8) >= 1 )THEN
        FNDOX(MLTM)=OUTDIR//'DOXTS' // CNTMSR(MLTM) // '.OUT'
        OPEN(41,FILE=FNDOX(MLTM),STATUS='UNKNOWN')
        CLOSE(41,STATUS='DELETE')
        OPEN(41,FILE=FNDOX(MLTM),STATUS='UNKNOWN')
        WRITE (41,100) TITLE4
        WRITE (41,101) CLTMSR(MLTM)
        WRITE (41,103)ILTMSR(MLTM),JLTMSR(MLTM)
        WRITE (41,102) CTUNIT
        CLOSE(41)
        FNTOC(MLTM)=OUTDIR//'TOCTS' // CNTMSR(MLTM) // '.OUT'
        OPEN(42,FILE=FNTOC(MLTM),STATUS='UNKNOWN')
        CLOSE(42,STATUS='DELETE')
        OPEN(42,FILE=FNTOC(MLTM),STATUS='UNKNOWN')
        WRITE (42,100) TITLE4
        WRITE (42,101) CLTMSR(MLTM)
        WRITE (42,103)ILTMSR(MLTM),JLTMSR(MLTM)
        WRITE (42,102) CTUNIT
        CLOSE(42)
        FNNHX(MLTM)=OUTDIR//'NHXTS' // CNTMSR(MLTM) // '.OUT'
        OPEN(43,FILE=FNNHX(MLTM),STATUS='UNKNOWN')
        CLOSE(43,STATUS='DELETE')
        OPEN(43,FILE=FNNHX(MLTM),STATUS='UNKNOWN')
        WRITE (43,100) TITLE4
        WRITE (43,101) CLTMSR(MLTM)
        WRITE (43,103)ILTMSR(MLTM),JLTMSR(MLTM)
        WRITE (43,102) CTUNIT
        CLOSE(43)
      ENDIF
      IF( ISTRAN(5) >= 1 )THEN
        DO NT=1,NTOX
          FNTOX(MLTM,NT)=OUTDIR//'TOX' // CNTOX(NT) // 'TS' // CNTMSR(MLTM) // '.OUT'
          OPEN(51,FILE=FNTOX(MLTM,NT),STATUS='UNKNOWN')
          CLOSE(51,STATUS='DELETE')
          OPEN(51,FILE=FNTOX(MLTM,NT),STATUS='UNKNOWN')
          WRITE (51,100) TITLE5
          WRITE (51,101) CLTMSR(MLTM)
          WRITE (51,103)ILTMSR(MLTM),JLTMSR(MLTM)
          WRITE (51,102) CTUNIT
          CLOSE(51)
          FNTXWT(MLTM,NT)=OUTDIR//'TXWT' // CNTOX(NT) // 'TS' // CNTMSR(MLTM) // '.OUT'
          OPEN(51,FILE=FNTXWT(MLTM,NT),STATUS='UNKNOWN')
          CLOSE(51,STATUS='DELETE')
          OPEN(51,FILE=FNTXWT(MLTM,NT),STATUS='UNKNOWN')
          WRITE (51,100) TFNTXWT
          WRITE (51,101) CLTMSR(MLTM)
          WRITE (51,103)ILTMSR(MLTM),JLTMSR(MLTM)
          WRITE (51,102) CTUNIT
          CLOSE(51)
          FNTXWF(MLTM,NT)=OUTDIR//'TXWF' // CNTOX(NT) // 'TS' // CNTMSR(MLTM) // '.OUT'
          OPEN(51,FILE=FNTXWF(MLTM,NT),STATUS='UNKNOWN')
          CLOSE(51,STATUS='DELETE')
          OPEN(51,FILE=FNTXWF(MLTM,NT),STATUS='UNKNOWN')
          WRITE (51,100) TFNTXWF
          WRITE (51,101) CLTMSR(MLTM)
          WRITE (51,103)ILTMSR(MLTM),JLTMSR(MLTM)
          WRITE (51,102) CTUNIT
          CLOSE(51)
          FNTXWC(MLTM,NT)=OUTDIR//'TXWC' // CNTOX(NT) // 'TS' // CNTMSR(MLTM) // '.OUT'
          OPEN(51,FILE=FNTXWC(MLTM,NT),STATUS='UNKNOWN')
          CLOSE(51,STATUS='DELETE')
          OPEN(51,FILE=FNTXWC(MLTM,NT),STATUS='UNKNOWN')
          WRITE (51,100) TFNTXWC
          WRITE (51,101) CLTMSR(MLTM)
          WRITE (51,103)ILTMSR(MLTM),JLTMSR(MLTM)
          WRITE (51,102) CTUNIT
          CLOSE(51)
          FNTXWP(MLTM,NT)=OUTDIR//'TXWP' // CNTOX(NT) // 'TS' // CNTMSR(MLTM) // '.OUT'
          OPEN(51,FILE=FNTXWP(MLTM,NT),STATUS='UNKNOWN')
          CLOSE(51,STATUS='DELETE')
          OPEN(51,FILE=FNTXWP(MLTM,NT),STATUS='UNKNOWN')
          WRITE (51,100) TFNTXWP
          WRITE (51,101) CLTMSR(MLTM)
          WRITE (51,103)ILTMSR(MLTM),JLTMSR(MLTM)
          WRITE (51,102) CTUNIT
          CLOSE(51)
          FNTXBT(MLTM,NT)=OUTDIR//'TXBT' // CNTOX(NT) // 'TS' // CNTMSR(MLTM) // '.OUT'
          OPEN(51,FILE=FNTXBT(MLTM,NT),STATUS='UNKNOWN')
          CLOSE(51,STATUS='DELETE')
          OPEN(51,FILE=FNTXBT(MLTM,NT),STATUS='UNKNOWN')
          WRITE (51,100) TFNTXBT
          WRITE (51,101) CLTMSR(MLTM)
          WRITE (51,103)ILTMSR(MLTM),JLTMSR(MLTM)
          WRITE (51,102) CTUNIT
          CLOSE(51)
          FNTXBF(MLTM,NT)=OUTDIR//'TXBF' // CNTOX(NT) // 'TS' // CNTMSR(MLTM) // '.OUT'
          OPEN(51,FILE=FNTXBF(MLTM,NT),STATUS='UNKNOWN')
          CLOSE(51,STATUS='DELETE')
          OPEN(51,FILE=FNTXBF(MLTM,NT),STATUS='UNKNOWN')
          WRITE (51,100) TFNTXBF
          WRITE (51,101) CLTMSR(MLTM)
          WRITE (51,103)ILTMSR(MLTM),JLTMSR(MLTM)
          WRITE (51,102) CTUNIT
          CLOSE(51)
          FNTXBC(MLTM,NT)=OUTDIR//'TXBC' // CNTOX(NT) // 'TS' // CNTMSR(MLTM) // '.OUT'
          OPEN(51,FILE=FNTXBC(MLTM,NT),STATUS='UNKNOWN')
          CLOSE(51,STATUS='DELETE')
          OPEN(51,FILE=FNTXBC(MLTM,NT),STATUS='UNKNOWN')
          WRITE (51,100) TFNTXBC
          WRITE (51,101) CLTMSR(MLTM)
          WRITE (51,103)ILTMSR(MLTM),JLTMSR(MLTM)
          WRITE (51,102) CTUNIT
          CLOSE(51)
          FNTXBP(MLTM,NT)=OUTDIR//'TXBP' // CNTOX(NT) // 'TS' // CNTMSR(MLTM) // '.OUT'
          OPEN(51,FILE=FNTXBP(MLTM,NT),STATUS='UNKNOWN')
          CLOSE(51,STATUS='DELETE')
          OPEN(51,FILE=FNTXBP(MLTM,NT),STATUS='UNKNOWN')
          WRITE (51,100) TFNTXBP
          WRITE (51,101) CLTMSR(MLTM)
          WRITE (51,103)ILTMSR(MLTM),JLTMSR(MLTM)
          WRITE (51,102) CTUNIT
          CLOSE(51)
        ENDDO
      ENDIF
    ENDIF    ! ***  MTMSRC(MLTM) == 1

    IF( MTMSRA(MLTM) == 1 )THEN
      FNAVV(MLTM)=OUTDIR//'AVVTS' // CNTMSR(MLTM) // '.OUT'
      OPEN(61,FILE=FNAVV(MLTM),STATUS='UNKNOWN')
      CLOSE(61,STATUS='DELETE')
      OPEN(61,FILE=FNAVV(MLTM),STATUS='UNKNOWN')
      WRITE (61,100) TITLE6
      WRITE (61,101) CLTMSR(MLTM)
      WRITE (61,103)ILTMSR(MLTM),JLTMSR(MLTM)
      WRITE (61,102) CTUNIT
      CLOSE(61)

      FNAVB(MLTM)=OUTDIR//'AVBTS' // CNTMSR(MLTM) // '.OUT'
      OPEN(71,FILE=FNAVB(MLTM),STATUS='UNKNOWN')
      CLOSE(71,STATUS='DELETE')
      OPEN(71,FILE=FNAVB(MLTM),STATUS='UNKNOWN')
      WRITE (71,100) TITLE7
      WRITE (71,101) CLTMSR(MLTM)
      WRITE (71,103)ILTMSR(MLTM),JLTMSR(MLTM)
      WRITE (71,102) CTUNIT
      CLOSE(71)
    ENDIF

    IF( MTMSRP(MLTM) == 1 )THEN
      FNSEL(MLTM)=OUTDIR//'SELTS' // CNTMSR(MLTM) // '.OUT'
      OPEN(11,FILE=FNSEL(MLTM),STATUS='UNKNOWN')
      CLOSE(11,STATUS='DELETE')
      OPEN(11,FILE=FNSEL(MLTM),STATUS='UNKNOWN')
      WRITE (11,100) TITLE11
      WRITE (11,101) CLTMSR(MLTM)
      WRITE (11,103)ILTMSR(MLTM),JLTMSR(MLTM)
      WRITE (11,102) CTUNIT
      CLOSE(11)
    ENDIF
    IF( MTMSRUE(MLTM) == 1 )THEN
      FNUVE(MLTM)=OUTDIR//'UVETS' // CNTMSR(MLTM) // '.OUT'
      OPEN(21,FILE=FNUVE(MLTM),STATUS='UNKNOWN')
      CLOSE(21,STATUS='DELETE')
      OPEN(21,FILE=FNUVE(MLTM),STATUS='UNKNOWN')
      WRITE (21,100) TITLE12
      WRITE (21,101) CLTMSR(MLTM)
      WRITE (21,103)ILTMSR(MLTM),JLTMSR(MLTM)
      WRITE (21,102) CTUNIT
      CLOSE(21)
    ENDIF
    IF( MTMSRUT(MLTM) == 1 )THEN
      FNUVT(MLTM)=OUTDIR//'UVTTS' // CNTMSR(MLTM) // '.OUT'
      OPEN(31,FILE=FNUVT(MLTM),STATUS='UNKNOWN')
      CLOSE(31,STATUS='DELETE')
      OPEN(31,FILE=FNUVT(MLTM),STATUS='UNKNOWN')
      WRITE (31,100) TITLE13
      WRITE (31,101) CLTMSR(MLTM)
      WRITE (31,103)ILTMSR(MLTM),JLTMSR(MLTM)
      WRITE (31,102) CTUNIT
      CLOSE(31)
    ENDIF
    IF( MTMSRU(MLTM) == 1 )THEN
      FNU3D(MLTM)=OUTDIR//'U3DTS' // CNTMSR(MLTM) // '.OUT'
      OPEN(41,FILE=FNU3D(MLTM),STATUS='UNKNOWN')
      CLOSE(41,STATUS='DELETE')
      OPEN(41,FILE=FNU3D(MLTM),STATUS='UNKNOWN')
      WRITE (41,100) TITLE14
      WRITE (41,101) CLTMSR(MLTM)
      WRITE (41,103)ILTMSR(MLTM),JLTMSR(MLTM)
      WRITE (41,102) CTUNIT
      CLOSE(41)
      FNV3D(MLTM)=OUTDIR//'V3DTS' // CNTMSR(MLTM) // '.OUT'
      OPEN(51,FILE=FNV3D(MLTM),STATUS='UNKNOWN')
      CLOSE(51,STATUS='DELETE')
      OPEN(51,FILE=FNV3D(MLTM),STATUS='UNKNOWN')
      WRITE (51,100) TITLE15
      WRITE (51,101) CLTMSR(MLTM)
      WRITE (51,103)ILTMSR(MLTM),JLTMSR(MLTM)
      WRITE (51,102) CTUNIT
      CLOSE(51)
    ENDIF
    IF( MTMSRQE(MLTM) == 1 )THEN
      FNQQE(MLTM)=OUTDIR//'QQETS' // CNTMSR(MLTM) // '.OUT'
      OPEN(61,FILE=FNQQE(MLTM),STATUS='UNKNOWN')
      CLOSE(61,STATUS='DELETE')
      OPEN(61,FILE=FNQQE(MLTM),STATUS='UNKNOWN')
      WRITE (61,100) TITLE16
      WRITE (61,101) CLTMSR(MLTM)
      WRITE (61,103)ILTMSR(MLTM),JLTMSR(MLTM)
      WRITE (61,102) CTUNIT
      CLOSE(61)
    ENDIF
    IF( MTMSRQ(MLTM) == 1 )THEN
      FNQ3D(MLTM)=OUTDIR//'Q3DTS' // CNTMSR(MLTM) // '.OUT'
      OPEN(71,FILE=FNQ3D(MLTM),STATUS='UNKNOWN')
      CLOSE(71,STATUS='DELETE')
      OPEN(71,FILE=FNQ3D(MLTM),STATUS='UNKNOWN')
      WRITE (71,100) TITLE17
      WRITE (71,101) CLTMSR(MLTM)
      WRITE (71,103)ILTMSR(MLTM),JLTMSR(MLTM)
      WRITE (71,102) CTUNIT
      CLOSE(71)
    ENDIF
  ENDDO
  !
  JSTMSR=0
  !
  !----------------------------------------------------------------------C
  !
    300 CONTINUE
  !
  !----------------------------------------------------------------------C
  !
  IF( ISDYNSTP == 0 )THEN
    TIME=(DT*FLOAT(N)+TCON*TBEGIN)/TCTMSR
  ELSE
    TIME=TIMESEC/TCTMSR
  ENDIF
  !
  FOURDPI=4./PI
  !
  ! **  STEP CURRENT TIME INTERVALS FOR WRITE SCENARIOS
  !
  DO NTSSS=1,NTSSTSP
   DO MTSSS=1,MTSSTSP(NTSSS)
    IF( TIME >= TSSTRT(MTSSS,NTSSS) )THEN
    IF( TIME <= TSSTOP(MTSSS,NTSSS) )THEN
      MTSCUR(NTSSS)=MTSSS
    ENDIF
    ENDIF
   ENDDO
  ENDDO
  !
  IF( MDCHH >= 1 )THEN
    DO L=2,LA
      QCHANUIJ(L)=0.0
      QCHANVIJ(L)=0.0
    ENDDO
    DO NMD=1,MDCHH
      LMDCHHT=LMDCHH(NMD)
      LMDCHUT=LMDCHU(NMD)
      LMDCHVT=LMDCHV(NMD)
      IF( MDCHTYP(NMD) == 1 )THEN
        QCHANUIJ(LMDCHHT)=QCHANUIJ(LMDCHHT)+QCHANU(NMD)
        QCHANUIJ(LMDCHUT)=QCHANUIJ(LMDCHUT)-QCHANU(NMD)
      ENDIF
      IF( MDCHTYP(NMD) == 2 )THEN
        QCHANVIJ(LMDCHHT)=QCHANVIJ(LMDCHHT)+QCHANV(NMD)
        QCHANVIJ(LMDCHVT)=QCHANVIJ(LMDCHVT)-QCHANV(NMD)
      ENDIF
      IF( MDCHTYP(NMD) == 3 )THEN
        QCHANUIJ(LMDCHHT)=QCHANUIJ(LMDCHHT)+QCHANU(NMD)
        QCHANUIJ(LMDCHUT)=QCHANUIJ(LMDCHUT)-QCHANU(NMD)
        QCHANVIJ(LMDCHHT)=QCHANVIJ(LMDCHHT)+QCHANV(NMD)
        QCHANVIJ(LMDCHVT)=QCHANVIJ(LMDCHVT)-QCHANV(NMD)
      ENDIF
    ENDDO
  ENDIF
  !
  DO MLTM=1,MLTMSR
   NTSSS=NTSSSS(MLTM)
   MTSCC=MTSCUR(NTSSS)
   IF( TIME >= TSSTRT(MTSCC,NTSSS) )THEN
   IF( TIME <= TSSTOP(MTSCC,NTSSS) )THEN
    I=ILTMSR(MLTM)
    J=JLTMSR(MLTM)
    L=LIJ(I,J)
    LN=LNC(L)
    IF( MTMSRC(MLTM) == 1 )THEN
      IF( ISTRAN(1) >= 1 )THEN
        OPEN(11,FILE=FNSAL(MLTM),POSITION='APPEND')
        WRITE(11,201)TIME,(SAL(L,K),K=1,KC)
        CLOSE(11)
      ENDIF
      IF( ISTRAN(2) >= 1 )THEN
        OPEN(21,FILE=FNTEM(MLTM),POSITION='APPEND')
        WRITE(21,201)TIME,(TEM(L,K),K=1,KC)
        CLOSE(21)
      ENDIF
      IF( ISTRAN(3) >= 1 )THEN
        OPEN(31,FILE=FNDYE(MLTM),POSITION='APPEND')
        WRITE(31,201)TIME,(DYE(L,K),K=1,KC)
        CLOSE(31)
      ENDIF
      IF( ISTRAN(4) >= 1 )THEN
        OPEN(31,FILE=FNSFL(MLTM),POSITION='APPEND')
        WRITE(31,201)TIME,(SFL(L,K),K=1,KC)
        CLOSE(31)
      ENDIF
      IF( ISTRAN(6) >= 1 )THEN
        OPEN(41,FILE=FNSED(MLTM),POSITION='APPEND')
        IF( ISNDAL == 2 )THEN
          SEDBTMP=SEDBT(L,KBT(L))+SEDBT(L,KBT(L)-1)
        ELSE
          SEDBTMP=SEDBT(L,KBT(L))
        ENDIF
        WRITE(41,201)TIME,SEDBTMP,(SEDT(L,K),K=1,KC)
        CLOSE(41)
      ENDIF
      IF( ISTRAN(7) >= 1 )THEN
        DO NX=1,NSND
          OPEN(41,FILE=FNSND(MLTM,NX),POSITION='APPEND')
          IF( ISNDAL == 2 )THEN
            SNDBTMP=SNDB(L,KBT(L),NX)+SNDB(L,KBT(L)-1,NX)
          ELSE
            SNDBTMP=SNDB(L,KBT(L),NX)
          ENDIF
          WRITE(41,201)TIME,SNDBTMP,(SND(L,K,NX),K=1,KC), &
                                   SNDEQSAV(L,NX)
          CLOSE(41)
        ENDDO
        DO NX=1,NSND
          OPEN(41,FILE=FNSBL(MLTM,NX),POSITION='APPEND')
  !            CQBEDLOADX=0.
  !            CQBEDLOADY=0.
          IF( UHDYE(L) /= 0.0 ) CQBEDLOADX(L,NX) = QSBDLDX(L,NX)/UHDYE(L)
          IF( VHDXE(L) /= 0.0 ) CQBEDLOADY(L,NX) = QSBDLDY(L,NX)/VHDXE(L)
          WRITE(41,201)TIME,QSBDLDX(L,NX),QSBDLDY(L,NX), &
                 CQBEDLOADX(L,NX),CQBEDLOADY(L,NX),SNDFBL(L,NX)
          CLOSE(41)
        ENDDO
  !            OPEN(41,FILE=FNSND(MLTM,NX),POSITION='APPEND')
  !            WRITE(41,201)TIME,SNDBT(L,KBT(L)),(SNDT(L,K),K=1,KC)
  !            CLOSE(41)
      ENDIF
      IF( ISTRAN(6) >= 1 .OR. ISTRAN(7) >= 1 )THEN
        OPEN(41,FILE=FNBED(MLTM),POSITION='APPEND')
        KTMP=KBT(L)
        KTMP1=KBT(L)-1
        KTMP1=MAX(KTMP1,1)
        NSXD=NSED+NSND
  !            WRITE(41,221)TIME,KTMP,HBED(L,KTMP),HBED(L,KTMP1),
        WRITE(41,221)TIME,KTMP,HBED(L,KTMP),HBED(L,KTMP1),HBEDA(L), &
                VDRBED(L,KTMP),(VFRBED(L,KTMP,NX),NX=1,NSXD)
        CLOSE(41)
      ENDIF
      IF( ISTRAN(8) >= 1 )THEN
        OPEN(41,FILE=FNDOX(MLTM),POSITION='APPEND')
        WRITE(41,201)TIME,(WQV(L,K,19),K=1,KC)
        CLOSE(41)
        OPEN(42,FILE=FNTOC(MLTM),POSITION='APPEND')
        WRITE(42,201)TIME,(WQV(L,K,6),K=1,KC)
        CLOSE(42)
        OPEN(43,FILE=FNNHX(MLTM),POSITION='APPEND')
        WRITE(43,201)TIME,(WQV(L,K,14),K=1,KC)
        CLOSE(43)
      ENDIF
      IF( ISTRAN(5) >= 1 )THEN
        DO NT=1,NTOX
  !
          OPEN(51,FILE=FNTOX(MLTM,NT),POSITION='APPEND')
          OPEN(52,FILE=FNTXWT(MLTM,NT),POSITION='APPEND')
          OPEN(53,FILE=FNTXWF(MLTM,NT),POSITION='APPEND')
          OPEN(54,FILE=FNTXWC(MLTM,NT),POSITION='APPEND')
          OPEN(55,FILE=FNTXWP(MLTM,NT),POSITION='APPEND')
          OPEN(56,FILE=FNTXBT(MLTM,NT),POSITION='APPEND')
          OPEN(57,FILE=FNTXBF(MLTM,NT),POSITION='APPEND')
          OPEN(58,FILE=FNTXBC(MLTM,NT),POSITION='APPEND')
          OPEN(59,FILE=FNTXBP(MLTM,NT),POSITION='APPEND')
  !
          NDOC=NSED+NSND+1
          DO K=1,KC
            SEDSND(K)=SEDT(L,K)+SNDT(L,K)
            TXWF(K)=TOXFDFW(L,K,NT)*TOX(L,K,NT)
            TXWC(K)=TOXCDFW(L,K,NT)*TOX(L,K,NT)
            TXWP(K)=TOXPFTW(L,K,NT)*TOX(L,K,NT)
          ENDDO
          DO K=1,KB
            SEDSNDB(K)=SEDBT(L,K)+SNDBT(L,K)
            TXBF(K)=0.
            TXBC(K)=0.
            TXBP(K)=0.
            TXBT(K)=0.
          ENDDO
          DO K=1,KBT(L)
            PORH(K)=1.0/PORBED(L,K)
            TXBF(K)=TOXFDFB(L,K,NT)*TOXB(L,K,NT)/HBED(L,K)
            TXBC(K)=TOXCDFB(L,K,NT)*TOXB(L,K,NT)/HBED(L,K)
            TXBP(K)=TOXPFTB(L,K,NT)*TOXB(L,K,NT)/HBED(L,K)
            TXBT(K)=TOXB(L,K,NT)/HBED(L,K)
          ENDDO

          K=KBT(L)
          WRITE(51,201)TIME,TXBT(K),TXBF(K),TXBC(K),TXBP(K),TOX(L,KSZ(L),NT),TXWF(1),TXWC(1),TXWP(1)
          DO K=1,KC
            IF( SEDSND(K) > 0.0)TXWP(K)=1000.*TXWP(K)/SEDSND(K)
          ENDDO

          DO K=1,KBT(L)
            TXBP(K)=TXBP(K)*HBED(L,K)
            TXBF(K)=TXBF(K)*PORH(K)
            TXBC(K)=TXBC(K)*PORH(K)
            IF( SEDSNDB(K) > 0.0)TXBP(K)=1000.*TXBP(K)/SEDSNDB(K)
          ENDDO
          WRITE(52,201)TIME,(TOX(L,K,NT),K=1,KC)
          WRITE(53,201)TIME,(TXWF(K),K=1,KC)
          WRITE(54,201)TIME,(TXWC(K),K=1,KC)
          WRITE(55,201)TIME,(TXWP(K),K=1,KC)
          WRITE(56,201)TIME,(TXBT(K),K=1,KB)
          WRITE(57,201)TIME,(TXBF(K),K=1,KB)
          WRITE(58,201)TIME,(TXBC(K),K=1,KB)
          WRITE(59,201)TIME,(TXBP(K),K=1,KB)
  !
          CLOSE(51)
          CLOSE(52)
          CLOSE(53)
          CLOSE(54)
          CLOSE(55)
          CLOSE(56)
          CLOSE(57)
          CLOSE(58)
          CLOSE(59)
  !
        ENDDO
  !            DO NT=1,NTOX
  !            OPEN(51,FILE=FNTOXPF(MLTM,NT),POSITION='APPEND')
  !            WRITE(51,201)TIME,TOXPFTB(L,KBT(L),NT),
  !     &                  (TOXPFTW(L,K,NT),K=1,KC)
  !            CLOSE(51)
  !            ENDDO
  !            DO NT=1,NTOX
  !            OPEN(51,FILE=FNTOXFD(MLTM,NT),POSITION='APPEND')
  !            WRITE(51,201)TIME,TOXPFTB(L,KBT(L),NT),
  !     &                  (TOXPFTW(L,K,NT),K=1,KC)
  !            CLOSE(51)
  !            ENDDO
      ENDIF
    ENDIF    ! ***  MTMSRC(MLTM) == 1

    IF( MTMSRA(MLTM) == 1 )THEN
      OPEN(61,FILE=FNAVV(MLTM),POSITION='APPEND')
      OPEN(71,FILE=FNAVB(MLTM),POSITION='APPEND')
      DO K=1,KS
        ATMPP(K)=10000.*AV(L,K)*HP(L)
      ENDDO
      WRITE(61,201)TIME,(ATMPP(K),K=1,KS)
      DO K=1,KS
        ATMPP(K)=10000.*AB(L,K)*HP(L)
      ENDDO
      WRITE(71,201)TIME,(ATMPP(K),K=1,KS)
      CLOSE(61)
      CLOSE(71)
    ENDIF
    IF( MTMSRP(MLTM) == 1 )THEN
      OPEN(11,FILE=FNSEL(MLTM),POSITION='APPEND')
      PPTMP=HP(L)+BELV(L)
      TMPVAL=VDWASTE(L)/DXYP(L)
      IF( (ISTRAN(6)+ISTRAN(7)) > 0 )THEN
        WRITE(11,201)TIME,PPTMP,HP(L),BELV(L),HBEDA(L),ZELBEDA(L),TMPVAL,VDWASTE(L)
      ELSE
        WRITE(11,201)TIME,PPTMP,HP(L),BELV(L),0.,0.,TMPVAL,VDWASTE(L)
      ENDIF
      CLOSE(11)
    ENDIF
    IF( MTMSRUE(MLTM) == 1 )THEN
      OPEN(21,FILE=FNUVE(MLTM),POSITION='APPEND')
      UTMP1=50.*(UHDYE(LEC(L))+UHDYE(L))/(DYP(L)*HP(L))
      VTMP1=50.*(VHDXE(LN)+VHDXE(L))/(DXP(L)*HP(L))
      IF( SPB(L) == 0 )THEN
        UTMP1=2.*UTMP1
        VTMP1=2.*VTMP1
      ENDIF
      UTMP=CUE(L)*UTMP1+CVE(L)*VTMP1
      VTMP=CUN(L)*UTMP1+CVN(L)*VTMP1
      UTMP1=5000.*(TBX(LEC(L))+TBX(L))
      VTMP1=5000.*(TBY(LN)+TBY(L))
      TBEAST=CUE(L)*UTMP1+CVE(L)*VTMP1
      TBNORT=CUN(L)*UTMP1+CVN(L)*VTMP1
      TAUBDYN=10000.*TAUB(L)
      TAUB2=TAUBDYN*TAUBDYN
      IF( ISWAVE > 0 )THEN
        CURANG=ATAN2(VTMP,UTMP)
        TAUW1=10000.*QQWV1(L)
        TAUW2=10000.*QQWV2(L)
        TAUB2=TAUB2 + (TAUW2*TAUW2) + 2.*TAUBDYN*TAUW2*COS(CURANG-WV(L).DIR)
      END IF
      TAUB2=MAX(TAUB2,0.)
      TAUTOT=SQRT(TAUB2)
      VELMAG=UTMP*UTMP+VTMP*VTMP
      TMPDRAG=0.0

      IF( (ISTRAN(6)+ISTRAN(7)) > 0 )THEN
        TAUBSEDDYN=10000.*TAUBSED(L)
        TAUBSNDDYN=10000.*TAUBSND(L)
        IF( VELMAG > 0.0)TMPDRAG=TAUTOT/VELMAG
        WRITE(21,201)TIME,UTMP,VTMP,TBEAST,TBNORT,TAUBDYN,TAUBSEDDYN,TAUBSNDDYN
      ELSE
        WRITE(21,201)TIME,UTMP,VTMP,TBEAST,TBNORT,TAUBDYN,0.,0.
      ENDIF
      CLOSE(21)
    ENDIF
    IF( MTMSRUT(MLTM) == 1 )THEN
      OPEN(31,FILE=FNUVT(MLTM),POSITION='APPEND')
      QBEDLOADX=0.
      QBEDLOADY=0.
      CQBEDLOADXT=0.
      CQBEDLOADYT=0.
      DO NX=1,NSND
        QBEDLOADX=QBEDLOADX+QSBDLDX(L,NX)
        QBEDLOADY=QBEDLOADY+QSBDLDY(L,NX)
        CQBEDLOADXT=CQBEDLOADXT+CQBEDLOADX(L,NX)
        CQBEDLOADYT=CQBEDLOADYT+CQBEDLOADY(L,NX)
      ENDDO
  !          IF( UHDYE(L) /= 0.0)CQBEDLOADX=QBEDLOADX/UHDYE(L)
  !          IF( VHDXE(L) /= 0.0)CQBEDLOADY=QBEDLOADY/VHDXE(L)
      WRITE(31,201)TIME,UHDYE(L),VHDXE(L),QBEDLOADX,QBEDLOADY, &
                   CQBEDLOADXT,CQBEDLOADYT
      CLOSE(31)
    ENDIF
    IF( MTMSRU(MLTM) == 1 )THEN
      OPEN(41,FILE=FNU3D(MLTM),POSITION='APPEND')
      OPEN(51,FILE=FNV3D(MLTM),POSITION='APPEND')
      RUVTMP=50.
      IF( SPB(L) == 0 ) RUVTMP=100.
      DO K=1,KC
       UTMP1=RUVTMP*(U(LEC(L),K)+U(L,K))
       VTMP1=RUVTMP*(V(LN,K)+V(L,K))
       ATMPP(K)=CUE(L)*UTMP1+CVE(L)*VTMP1
       BTMP(K)=CUN(L)*UTMP1+CVN(L)*VTMP1
      ENDDO
      WRITE(41,201)TIME,(ATMPP(K),K=1,KC)
      WRITE(51,201)TIME,(BTMP(K),K=1,KC)
      CLOSE(41)
      CLOSE(51)
    ENDIF
    IF( MTMSRQE(MLTM) == 1 )THEN
      OPEN(61,FILE=FNQQE(MLTM),POSITION='APPEND')
      QRNT=DXYP(L)*RAINT(L)
      WRITE(61,201)TIME,QSUME(L),QRNT,EVAPSW(L),EVAPGW(L),QGW(L),QCHANUIJ(L),QCHANVIJ(L),QDWASTE(L),VDWASTE(L)
      CLOSE(61)
    ENDIF
    IF( MTMSRQ(MLTM) == 1 )THEN
      OPEN(71,FILE=FNQ3D(MLTM),POSITION='APPEND')
      WRITE (71,201)TIME,(QSUM(L,K),K=1,KC)
      CLOSE(71)
    ENDIF
   ENDIF
   ENDIF
  ENDDO
  !
  !**********************************************************************C
  !
    100 FORMAT(A80)
    101 FORMAT('  AT LOCATION  ',A20)
    102 FORMAT('  TIME IN FIRST COLUMN HAS UNITS OF ',A10)
    103 FORMAT('  CELL I,J = ',2I5)
    201 FORMAT(F12.5,12E12.4)
    221 FORMAT(F12.5,I5,14E12.4)
  !
  !**********************************************************************C
  !
  RETURN
END

max_stars_repo_name: dsi-llc/EFDCPlus

max_forks_repo_name: dsi-llc/EFDCPlus