diff --git a/src/gray.f b/src/gray.f
index d742660..2e07b0d 100644
--- a/src/gray.f
+++ b/src/gray.f
@@ -1,9 +1,9 @@
-      subroutine gray(ijetto, mr, mz, r, z, psin, psiax, psibnd,
+      subroutine gray(ijetto, mr, mz, mrd, r, z, psin, psiax, psibnd,
      .    rax, zax, nbnd, rbnd, zbnd, nrho, psijet, f, te, dne,
      .    zeff, qsf, powin, dpdv, jcd, pec, icd, ierr)
 c   input arguments
       integer ijetto, mr, mz, nbnd, nrho
-      real*8 r(mr), z(mz), psin(mr,mz)
+      real*8 r(mr), z(mz), psin(mrd,mz)
       real*8 psiax, psibnd, rax, zax
       real*8 rbnd(nbnd), zbnd(nbnd)
       real*8 psijet(nrho), f(nrho), qsf(nrho), te(nrho), dne(nrho)
@@ -12,15 +12,16 @@ c   output arguments
       real*8 dpdv(nrho), jcd(nrho), powin, pec, icd
       integer ierr
 c   local variables
-      real*8 fgray(nrho), qgray(nrho), jcdgry(nrho), icdgry
+c      real*8 fgray(nrho), qgray(nrho), jcdgry(nrho), icdgry
 c
 c   input arguments
 c
 c     ijetto  Equilibrium source (1 EFIT, 2 ESCO)
 c             If IJETTO=2, then PSIN values are valid only inside
-c             plasma boudary (PSIN=1 outside)
+c             plasma boudary (PSIN=0 outside)
 c     mr      Size of flux map grid in R direction
 c     mz      Size of flux map grid in Z direction
+c     mrd     Leading dimension of the psin(:,:) array, mrd>mr
 c     r       R coordinates of flux map grid points [m]
 c     z       Z coordinates of flux map grid points [m]
 c     psin    Normalised poloidal flux psin=(psi-psiax)/(psibnd-psiax)
@@ -61,22 +62,24 @@ c   JETTO coordinate system assumes toroidal angle increasing CW
 c   in GRAY toroidal angle increases CCW --> adapt signs on input data
 c
       do i=1,nrho
-        fgray(i)=-f(i)
-        qgray(i)=-qsf(i)
+        f(i)=-f(i)
+        qsf(i)=-qsf(i)
       end do
 c
 c   call main subroutine
 c
-      call gray_main(ijetto, mr, mz, r, z, psin, psiax, psibnd,
-     .    rax, zax, nbnd, rbnd, zbnd, nrho, psijet, fgray, te, dne,
-     .    zeff, qgray, powin, dpdv, jcdgry, pec, icdgry, ierr)
+      call gray_main(ijetto, mr, mz, r, z, psin(1:mr,:), psiax, psibnd,
+     .    rax, zax, nbnd, rbnd, zbnd, nrho, psijet, f, te, dne,
+     .    zeff, qsf, powin, dpdv, jcd, pec, icd, ierr)
 c
 c   adapt output data to JETTO convention on toroidal angle
 c
       do i=1,nrho
-        jcd(i)=-jcdgry(i)
+        jcd(i)=-jcd(i)
+        f(i)=-f(i)
+        qsf(i)=-qsf(i)
       end do
-      icd=-icdgry
+      icd=-icd
 c
       return
-      end
\ No newline at end of file
+      end