! source file: /sfs/fs6/home-geomar/smomw258/UVOK_1.1/Kiel_Feb_2019/source/mom/matrix.F
      SUBROUTINE MATRIX_INIT( )

      IMPLICIT NONE

      RETURN
      END

      SUBROUTINE MATRIX_STORE_EXPLICIT(
     I                  ntile,tracer,twodt )

      IMPLICIT NONE

      include "size.h"
      include "matrix.h"

      integer ntile
      real tracer(imt,km,jmw), twodt(km)

      character(120) :: fname, new_file_name
      integer k

!     compute and accumulate explicit matrix
	  do k=1,km
		Aexp(:,k,:,ntile) = Aexp(:,k,:,ntile) +
     &       ((tracer(:,k,:)-tile(:,k,:,ntile))/twodt(k))
      enddo
	  if (ntile.eq.numtiles) then
	    stepCount = stepCount + 1
	    expMatrixCounter = expMatrixCounter + 1
	    write(*,*)'DEBUG3: accumulating explicit matrix ',
     &   expMatrixCounter, stepCount, twodt(1),twodt(2)
      endif

!     reset tracer to tiles/initial condition
	  tracer(:,:,:)=tile(:,:,:,ntile)

      RETURN
      END

      SUBROUTINE MATRIX_STORE_IMPLICIT(
     I                  ntile,tracer )

      IMPLICIT NONE

      include "size.h"
      include "matrix.h"

      integer ntile
      real tracer(imt,km,jmw)

      character(120) :: fname, new_file_name

!     compute and accumulate implicit matrix
	  Aimp(:,:,:,ntile) = Aimp(:,:,:,ntile) + tracer(:,:,:)

	  if (ntile.eq.numtiles) then
	    impMatrixCounter = impMatrixCounter + 1
	    write(*,*)'DEBUG4: accumulating implicit matrix ',
     &   impMatrixCounter
      endif

      RETURN
      END

      SUBROUTINE MATRIX_WRITE(writeExpMatrix, writeImpMatrix)

      IMPLICIT NONE

      include "size.h"
      include "pconst.h"
      include "matrix.h"

      logical writeExpMatrix, writeImpMatrix

      integer ntile
      character(120) :: fn, fname, new_file_name
      logical exists
      real recipExpMatrixCounter, recipImpMatrixCounter

C      writeExpMatrix = .TRUE.
C      writeImpMatrix = .TRUE.

      IF (writeExpMatrix) THEN
		recipExpMatrixCounter = c1/expMatrixCounter
		write(*,*)'DEBUG1 ',recipExpMatrixCounter,expMatrixCounter
        Aexp(:,:,:,:) = Aexp(:,:,:,:)*recipExpMatrixCounter
        do ntile=1,numtiles
		  write(*,*) 'Writing explicit matrix :',ntile,
     &        expMatrixWriteCount, expMatrixCounter
		  write(fn,'(A9,I2.2,A4)') 'MATRIXEXP',ntile,'.bin'
		  fname = new_file_name (fn)
		  call write_tracer_array(fname,Aexp(1,1,1,ntile),
     &	                          expMatrixWriteCount)
        enddo
		expMatrixCounter=0
		expMatrixWriteCount = expMatrixWriteCount + 1

C       Reset explicit matrix
        Aexp(:,:,:,:) = c0
      ENDIF

      IF (writeImpMatrix) THEN
		recipImpMatrixCounter = c1/impMatrixCounter
		write(*,*)'DEBUG2 ',recipImpMatrixCounter,impMatrixCounter
        Aimp(:,:,:,:) = Aimp(:,:,:,:)*recipImpMatrixCounter
        do ntile=1,numtiles
		  write(*,*) 'Writing implicit matrix :',ntile,
     &        impMatrixWriteCount, impMatrixCounter
		  write(fn,'(A9,I2.2,A4)') 'MATRIXIMP',ntile,'.bin'
		  fname = new_file_name (fn)
		  call write_tracer_array(fname,Aimp(1,1,1,ntile),
     &	                          impMatrixWriteCount)
        enddo
		impMatrixCounter=0
		impMatrixWriteCount = impMatrixWriteCount + 1

C       Reset implicit matrix
        Aimp(:,:,:,:) = c0
      ENDIF

      RETURN
      END

C======================================================
C This subroutine reads a model tracer array from file
C======================================================
      SUBROUTINE READ_TRACER_ARRAY(filename,fld,iSlice)

      implicit none

      include "size.h"
      include "param.h"

      integer iSlice
      real*8 fld(1:imt,1:km,1:jmt)
      character*(*) filename

      real*8 tmp(imtm2,jmtm2,km)
      integer j, k

      call read_r8_field(imtm2*jmtm2,km,iSlice,tmp,filename)

      do j=2,jmtm1
        do k=1,km
          fld(2:imtm1,k,j)=tmp(1:imtm2,j-1,k)
        enddo
      enddo
      fld(2:imtm1,1:km,1)=fld(2:imtm1,1:km,2)
      fld(2:imtm1,1:km,jmt)=fld(2:imtm1,1:km,jmtm1)
	  do j=1,jmt
		call setbcx (fld(1:imt,1:km,j),imt,km)
	  enddo

      END

C======================================================
C This subroutine writes a model tracer array from file
C======================================================
      SUBROUTINE WRITE_TRACER_ARRAY(filename,fld,iSlice)

      implicit none

      include "size.h"
      include "param.h"

      integer iSlice
      real*8 fld(1:imt,1:km,1:jmt)
      character*(*) filename

      real*8 tmp(imtm2,jmtm2,km)
      integer j, k

      do j=2,jmtm1
        do k=1,km
          tmp(1:imtm2,j-1,k)=fld(2:imtm1,k,j)
        enddo
      enddo

      call write_r8_field(imtm2*jmtm2,km,iSlice,tmp,filename)

      END