INTEGER, PARAMETER :: IMT=1238,JMT=1046,KM=100

REAL*4, DIMENSION(IMT,JMT,2) :: psixy
REAL*4, DIMENSION(JMT,KM,2) :: psiyz
REAL*4, DIMENSION(JMT,MR,4,2) :: psiyr

There is a "2" because I was analysing two differenty runs and comparing them
So remove it or set nr=1 further down

file=trim(dir(nr))//trim(run((nr)))//'_run.csv'
print *,trim(file)
close(56)
open(56,file=trim(file),status='old')
goto 141
143 continue
print *,ntrac,rlon,rlat,zz,subvol,tt,lbox,salt,temp,dens
stop 395672

141 continue
read(56,*,end=142,err=143) ntrac,rlon,rlat,zz,subvol,tt,lbox,salt,temp!,dens

call rhoc(1,1,1,0., temp,salt,dens)


!if(nselect(ntrac,nr)<2) goto 141
i=int(rlon)+1 ; j=int(rlat)+1 ; k=int(zz)+1

mz=INT( ( depth(INT(zz)) + (zz - INT(zz))*dz(INT(zz)+1) )*float(MR)/dmax +1 )! zz -> m
if(mz< 1) mz=1
if(mz>MR) mz=MR
mt=int((temp-tmin)/dtemp)+1
if(mt< 1) mt=1
if(mt>MR) mt=MR
ms=int((salt-smin)/dsalt)+1
if(ms<1 ) ms=1
if(ms>MR) ms=MR
md=int((dens-rmin)/drho)+1
if(md<1 ) md=1
if(md>MR) md=MR

tt=tt/(24.*3600.)

In case you want to compute the residence times you can use this but you will also then need the <file>_ini.csv and <file>_out.csv

age   (i,j,k,nr)=age   (i,j,k,nr) + sngl((tt- trin(ntrac,nr,4))               *subvol )
volage(i,j,k,nr)=volage(i,j,k,nr) + sngl(subvol)
res   (i,j,k,nr)=res   (i,j,k,nr) + sngl((trut(ntrac,nr,4) - tt)              *subvol )
tra   (i,j,k,nr)=tra   (i,j,k,nr) + sngl((trut(ntrac,nr,4) - trin(ntrac,nr,4))*subvol )
volres(i,j,k,nr)=volres(i,j,k,nr) + sngl(subvol)

nt=int( tt - trin(ntrac,nr,4) )+1
if(1<=nt .and. nt<=NTIME) restime(nt,2,nr)=restime(nt,2,nr)+sngl(subvol)
nt=int( trut(ntrac,nr,4) - tt )+1
if(1<=nt .and. nt<=NTIME) restime(nt,3,nr)=restime(nt,3,nr)+sngl(subvol)
nt=int( trut(ntrac,nr,4) - trin(ntrac,nr,4) )+1
if(1<=nt .and. nt<=NTIME) restime(nt,4,nr)=restime(nt,4,nr)+sngl(subvol)

! stream functions ___________________________________________________________________
 box:if(lbox==3 .or. lbox==4) then


  if(lbox==3) then
   voltran=sngl(subvol)
  elseif(lbox==4) then
   voltran=-sngl(subvol)
  endif
  psixy(i,j,nr)=psixy(i,j,nr)+voltran
  summa(j,nr)=summa(j,nr)+voltran

  if(i>=ICEN .and. j>=JCEN) then
   psiyz(j,k ,  nr)=psiyz(j,k ,  nr)+voltran ! model levels
   psiyr(j,mz,1,nr)=psiyr(j,mz,1,nr)+voltran ! depth metres
   psiyr(j,mt,2,nr)=psiyr(j,mt,2,nr)+voltran ! temperature
   psiyr(j,ms,3,nr)=psiyr(j,ms,3,nr)+voltran ! salinity
   psiyr(j,md,4,nr)=psiyr(j,md,4,nr)+voltran ! density




goto 141
142 continue
close(56)



psixy=-psixy/365.
do i=IMT-1,2,-1
 psixy(i,:,:)=psixy(i,:,:)+psixy(i+1,:,:)
enddo

psiwmin=0. ; psismin=0. ; psiwmax=0. ; psismax=0.
do j=1,JYP
 do i=1,IXP
  if(psiwmin>psixy(i+ICW+1,j+JCS+1,1)) psiwmin=psixy(i+ICW+1,j+JCS+1,1)
  if(psiwmax<psixy(i+ICW+1,j+JCS+1,1)) psiwmax=psixy(i+ICW+1,j+JCS+1,1)
  if(psismin>psixy(i+ICW+1,j+JCS+1,2)) psismin=psixy(i+ICW+1,j+JCS+1,2)
  if(psismax<psixy(i+ICW+1,j+JCS+1,2)) psismax=psixy(i+ICW+1,j+JCS+1,2)
 enddo
enddo

print *,'maxmin of psixy water',psiwmin,psiwmax,' m3/s'
print *,'maxmin of psixy salt',psismin,psismax,' kg/s'

psixy(:,:,1)=psixy(:,:,1)*1.e-2 ! m3/s -> 100 m3/s eller cSv centiSverdrup
psixy(:,:,2)=psixy(:,:,2)*1.e-4 ! kg/s -> 10^4 kg/s


psiyz=-psiyz/365. ! *2.e-2
psiyr=-psiyr/365. ! *2.e-2


do k=KM-1,1,-1
 psiyz(:,k,  :)=psiyz(:,k,  :)+psiyz(:,k+1,  :)
 psiyr(:,k,:,:)=psiyr(:,k,:,:)+psiyr(:,k+1,:,:)
enddo



psiyz(:,:,1)=psiyz(:,:,1)*1.e-2 ! m3/s -> 100 m3/s

psiyr(:,:,:,1)=psiyr(:,:,:,1)*1.e-2 ! m3/s -> 100 m3/s