Index: IOIPSL/trunk/src/restcom.f90
===================================================================
--- IOIPSL/trunk/src/restcom.f90	(revision 2350)
+++ IOIPSL/trunk/src/restcom.f90	(revision 3188)
@@ -8,4 +8,5 @@
 USE netcdf
 !-
+USE defprec
 USE errioipsl, ONLY : ipslerr,ipsldbg, ipslout
 USE stringop
@@ -26,5 +27,6 @@
   MODULE PROCEDURE &
  &  restput_r3d, restput_r2d, restput_r1d, &
- &  restput_opp_r2d, restput_opp_r1d
+ &  restput_opp_r2d, restput_opp_r1d, restput_opp_r3d, &
+ &  restput_opp_r4d, restput_opp_r5d
 END INTERFACE
 !-
@@ -32,5 +34,6 @@
   MODULE PROCEDURE &
  &  restget_r3d,restget_r2d,restget_r1d, &
- &  restget_opp_r2d,restget_opp_r1d
+ &  restget_opp_r2d, restget_opp_r1d, restget_opp_r3d, &
+ &  restget_opp_r4d, restget_opp_r5d
 END INTERFACE
 !-
@@ -40,5 +43,6 @@
 !-
   INTEGER,PARAMETER :: &
- &  max_var=500, max_file=50, max_dim=NF90_MAX_VAR_DIMS
+ &  max_var=500, max_file=50, MAX_DIM=NF90_MAX_VAR_DIMS, RESTART_MAX_DIMS=7 !
+! RESTART_MAX_DIMS = Max dimensions in 1 variable IOIPSL can manage for restart files (> 3)
 !-
   CHARACTER(LEN=9),SAVE :: calend_str='unknown'
@@ -88,8 +92,19 @@
 !   ?ax_nb(if)
 !-
-  INTEGER,DIMENSION(max_file,max_dim,2),SAVE :: &
- &  xax_infs,yax_infs,zax_infs
-  INTEGER,DIMENSION(max_file),SAVE :: &
- &  xax_nb=0,yax_nb=0,zax_nb=0
+  INTEGER,DIMENSION(max_file,RESTART_MAX_DIMS,max_dim,2),SAVE :: &
+ &  ax_infs!,yax_infs,zax_infs,wax_infs
+  INTEGER,DIMENSION(max_file,RESTART_MAX_DIMS),SAVE :: &
+ &  ax_nb=0! ,yax_nb=0,zax_nb=0,wax_nb=0
+!-
+! Dimensions names, must be equal to RESTART_MAX_DIMS:
+! - DO not repeat any char
+! - ONLY 1 single char per dimension
+  CHARACTER, DIMENSION(RESTART_MAX_DIMS), PARAMETER :: RESTART_DIMS_NAMES = &
+     (/ 'x','y','z','l','m','n', 'o'/)
+
+!- Dimensions constants for basic dimensions
+  INTEGER, PARAMETER :: RESTART_DIM_X = 1
+  INTEGER, PARAMETER :: RESTART_DIM_Y = 2
+  INTEGER, PARAMETER :: RESTART_DIM_Z = 3
 !-
 ! Description of the time axes in the input and output files
@@ -158,5 +173,5 @@
 SUBROUTINE restini &
  & (fnamein,iim,jjm,lon,lat,llm,lev, &
- &  fnameout,itau,date0,dt,fid,owrite_time_in,domain_id)
+ &  fnameout,itau,date0,dt,fid,owrite_time_in,domain_id )
 !---------------------------------------------------------------------
 !- This subroutine sets up all the restart process.
@@ -203,4 +218,5 @@
 !-                  overwrite the time in the restart file
 !- domain_id      : Domain identifier
+!-
 !---------------------------------------------------------------------
   IMPLICIT NONE
@@ -242,6 +258,4 @@
  &   'and recompiling ioipsl.')
   ENDIF
-!-
-! 0.1 Define the open flags
 !-
   l_fi = (TRIM(fnamein)  /= 'NONE')
@@ -407,5 +421,5 @@
 ! then we open it in write mode
 !-
-  IF (l_rw) THEN; id = NF90_WRITE; ELSE; id = NF90_NOWRITE; ENDIF
+  IF (l_rw) THEN; id = IOR(NF90_WRITE,NF90_NETCDF4); ELSE; id = IOR(NF90_NOWRITE,NF90_NETCDF4); ENDIF
   iret = NF90_OPEN(fname,id,ncfid)
   IF (iret /= NF90_NOERR) THEN
@@ -440,11 +454,11 @@
       WRITE (*,*) "restopenin 0.0 dimname",id,TRIM(dimname(id))
     ENDIF
-    IF      (TRIM(dimname(id)) == 'x') THEN
+    IF      (TRIM(dimname(id)) == RESTART_DIMS_NAMES(RESTART_DIM_X)) THEN
       iread = dimlen(id)
       IF (l_dbg) WRITE (*,*) "iread",iread
-    ELSE IF (TRIM(dimname(id)) == 'y') THEN
+    ELSE IF (TRIM(dimname(id)) == RESTART_DIMS_NAMES(RESTART_DIM_Y)) THEN
       jread = dimlen(id)
       IF (l_dbg) WRITE (*,*) "jread",jread
-    ELSE IF (TRIM(dimname(id)) == 'z') THEN
+    ELSE IF (TRIM(dimname(id)) == RESTART_DIMS_NAMES(RESTART_DIM_Z)) THEN
       lread = dimlen(id)
       IF (l_dbg) WRITE (*,*) "lread",lread
@@ -468,22 +482,19 @@
       ENDIF
 !-----
-      xax_nb(fid) = 0
-      yax_nb(fid) = 0
-      zax_nb(fid) = 0
+      ax_nb(fid,:) = 0
 !-----
       DO id=1,nb_dim
-        IF      (dimname(id)(1:1) == 'x') THEN
-          xax_nb(fid) = xax_nb(fid)+1
-          xax_infs(fid,xax_nb(fid),1) = dimlen(id)
-          xax_infs(fid,xax_nb(fid),2) = id
-        ELSE IF (dimname(id)(1:1) == 'y') THEN
-          yax_nb(fid) = yax_nb(fid)+1
-          yax_infs(fid,yax_nb(fid),1) = dimlen(id)
-          yax_infs(fid,yax_nb(fid),2) = id
-        ELSE IF (dimname(id)(1:1) == 'z') THEN
-          zax_nb(fid) = zax_nb(fid)+1
-          zax_infs(fid,zax_nb(fid),1) = dimlen(id)
-          zax_infs(fid,zax_nb(fid),2) = id
+        ! Order is important
+        IF      (dimname(id)(1:1) == RESTART_DIMS_NAMES(id)) THEN
+          ax_nb(fid,id) = ax_nb(fid,id)+1
+          ax_infs(fid,id,ax_nb(fid,id),1) = dimlen(id)
+          ax_infs(fid,id,ax_nb(fid,id),2) = id
         ENDIF
+      ENDDO
+!-----
+!---  Initialize non basic dimensions
+!-----
+      DO id=4,RESTART_MAX_DIMS
+        ax_nb(fid,id) = 1
       ENDDO
     ENDIF
@@ -832,4 +843,5 @@
   CHARACTER(LEN=30) :: timenow
   LOGICAL :: l_dbg
+  INTEGER :: cmode, id
 !---------------------------------------------------------------------
   CALL ipsldbg (old_status=l_dbg)
@@ -840,5 +852,5 @@
 !- we will not even call restopenout
 !-
-  iret = NF90_CREATE(fname,cmode=or(NF90_NOCLOBBER,NF90_64BIT_OFFSET),ncid=ncfid)
+  iret = NF90_CREATE(fname,cmode=or(NF90_NOCLOBBER,NF90_NETCDF4),ncid=ncfid)
   IF (iret == -35) THEN
     CALL ipslerr (3,'restopenout',&
@@ -847,19 +859,25 @@
       & ' generate the new one with another name')
   ENDIF
-!-
-  iret = NF90_DEF_DIM(ncfid,'x',iim,x_id)
-  xax_nb(fid) = xax_nb(fid)+1
-  xax_infs(fid,xax_nb(fid),1) = iim
-  xax_infs(fid,xax_nb(fid),2) = x_id
-!-
-  iret = NF90_DEF_DIM(ncfid,'y',jjm,y_id)
-  yax_nb(fid)  = yax_nb(fid)+1
-  yax_infs(fid,yax_nb(fid),1) = jjm
-  yax_infs(fid,yax_nb(fid),2) = y_id
-!-
-  iret = NF90_DEF_DIM(ncfid,'z',llm,z_id)
-  zax_nb(fid) = zax_nb(fid)+1
-  zax_infs(fid,zax_nb(fid),1) = llm
-  zax_infs(fid,zax_nb(fid),2) = z_id
+!- Define basic dimensions (x, y, z)
+  iret = NF90_DEF_DIM(ncfid,RESTART_DIMS_NAMES(RESTART_DIM_X),iim,x_id)
+  ax_nb(fid,RESTART_DIM_X) = ax_nb(fid, RESTART_DIM_X)+1
+  ax_infs(fid,RESTART_DIM_X, ax_nb(fid, RESTART_DIM_X),1) = iim
+  ax_infs(fid,RESTART_DIM_X, ax_nb(fid, RESTART_DIM_X),2) = x_id
+!-
+  iret = NF90_DEF_DIM(ncfid,RESTART_DIMS_NAMES(RESTART_DIM_Y),jjm,y_id)
+  ax_nb(fid,RESTART_DIM_Y) = ax_nb(fid, RESTART_DIM_Y)+1
+  ax_infs(fid,RESTART_DIM_Y, ax_nb(fid, RESTART_DIM_Y),1) = jjm
+  ax_infs(fid,RESTART_DIM_Y, ax_nb(fid, RESTART_DIM_Y),2) = y_id
+!-
+  iret = NF90_DEF_DIM(ncfid,RESTART_DIMS_NAMES(RESTART_DIM_Z),llm,z_id)
+  ax_nb(fid,RESTART_DIM_Z) = ax_nb(fid, RESTART_DIM_Z)+1
+  ax_infs(fid,RESTART_DIM_Z, ax_nb(fid, RESTART_DIM_Z),1) = llm
+  ax_infs(fid,RESTART_DIM_Z, ax_nb(fid, RESTART_DIM_Z),2) = z_id
+!-
+!  Initialize non basic dimensions
+!-
+  DO id=4,RESTART_MAX_DIMS
+    ax_nb(fid,id) = 1
+  ENDDO
 !-
   iret = NF90_DEF_DIM(ncfid,'time',NF90_UNLIMITED,tdimid_out(fid))
@@ -1014,4 +1032,5 @@
   CHARACTER(LEN=7) :: topp
   LOGICAL :: l_dbg
+  INTEGER :: list_dims(2) ! x and y
 !---------------------------------------------------------------------
   CALL ipsldbg (old_status=l_dbg)
@@ -1021,7 +1040,7 @@
   req_sz = 1
   IF (nbindex == iim .AND. jjm <= 1 .AND. llm <= 1) THEN
-    IF (xax_infs(fid,1,1) > 0) req_sz = req_sz*xax_infs(fid,1,1)
-    IF (yax_infs(fid,1,1) > 0) req_sz = req_sz*yax_infs(fid,1,1)
-    IF (zax_infs(fid,1,1) > 0) req_sz = req_sz*zax_infs(fid,1,1)
+    IF (ax_infs(fid,RESTART_DIM_X,1,1) > 0) req_sz = req_sz*ax_infs(fid,RESTART_DIM_X,1,1)
+    IF (ax_infs(fid,RESTART_DIM_Y,1,1) > 0) req_sz = req_sz*ax_infs(fid,RESTART_DIM_Y,1,1)
+    IF (ax_infs(fid,RESTART_DIM_Z,1,1) > 0) req_sz = req_sz*ax_infs(fid,RESTART_DIM_Z,1,1)
   ELSE
     CALL ipslerr (3,'resget_opp_r1d', &
@@ -1039,7 +1058,8 @@
 ! 2.0 Here we get the variable from the restart file
 !-
+  list_dims = (/ ax_infs(fid,RESTART_DIM_X,1,1), ax_infs(fid,RESTART_DIM_Y,1,1) /)
+!-
   CALL restget_real &
-    (fid,vname_q,xax_infs(fid,1,1),yax_infs(fid,1,1), &
-     zax_infs(fid,1,1),itau,def_beha,buff_tmp2)
+    (fid,vname_q,list_dims,itau,def_beha,buff_tmp2)
 !-
 ! 4.0 Transfer the buffer obtained from the restart file
@@ -1083,4 +1103,5 @@
   CHARACTER(LEN=7) :: topp
   LOGICAL :: l_dbg
+  INTEGER :: list_dims(3)
 !---------------------------------------------------------------------
   CALL ipsldbg (old_status=l_dbg)
@@ -1090,6 +1111,6 @@
   req_sz = 1
   IF (nbindex == iim  .AND. llm <= 1) THEN
-    IF (xax_infs(fid,1,1) > 0) req_sz = req_sz*xax_infs(fid,1,1)
-    IF (yax_infs(fid,1,1) > 0) req_sz = req_sz*yax_infs(fid,1,1)
+    IF (ax_infs(fid,RESTART_DIM_X,1,1) > 0) req_sz = req_sz*ax_infs(fid,RESTART_DIM_X,1,1)
+    IF (ax_infs(fid,RESTART_DIM_Y,1,1) > 0) req_sz = req_sz*ax_infs(fid,RESTART_DIM_Y,1,1)
   ELSE
     CALL ipslerr (3,'resget_opp_r2d', &
@@ -1113,7 +1134,8 @@
 ! 2.0 Here we get the full variable from the restart file
 !-
+  list_dims = (/ ax_infs(fid,RESTART_DIM_X,1,1), ax_infs(fid,RESTART_DIM_Y,1,1), jjm /)
+!-
   CALL restget_real &
- & (fid,vname_q,xax_infs(fid,1,1),yax_infs(fid,1,1), &
- &  jjm,itau,def_beha,buff_tmp2)
+ & (fid,vname_q,list_dims,itau,def_beha,buff_tmp2)
 !-
 ! 4.0 Transfer the buffer obtained from the restart file
@@ -1139,4 +1161,250 @@
 END SUBROUTINE restget_opp_r2d
 !===
+SUBROUTINE restget_opp_r3d &
+ & (fid,vname_q,iim,jjm,llm,itau,def_beha, &
+ &  var,MY_OPERATOR,nbindex,ijndex)
+!---------------------------------------------------------------------
+!- This subroutine serves as an interface to restget_real
+!-
+!- Should work as restput_opp_r2d but the other way around !
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+  INTEGER :: fid
+  CHARACTER(LEN=*) :: vname_q
+  INTEGER :: iim,jjm,llm,itau
+  LOGICAL def_beha
+  REAL :: var(:,:,:)
+  CHARACTER(LEN=*) :: MY_OPERATOR
+  INTEGER :: nbindex,ijndex(nbindex)
+!-
+  INTEGER :: ll,jj,req_sz,ist,var_sz,siz1
+  REAL :: scal
+  CHARACTER(LEN=7) :: topp
+  LOGICAL :: l_dbg
+  INTEGER :: list_dims(4)
+!---------------------------------------------------------------------
+  CALL ipsldbg (old_status=l_dbg)
+!-
+! 0.0 What size should be the data in the file
+!-
+  req_sz = 1
+  IF (nbindex == iim) THEN
+    IF (ax_infs(fid,RESTART_DIM_X,1,1) > 0) req_sz = req_sz*ax_infs(fid,RESTART_DIM_X,1,1)
+    IF (ax_infs(fid,RESTART_DIM_Y,1,1) > 0) req_sz = req_sz*ax_infs(fid,RESTART_DIM_Y,1,1)
+  ELSE
+    CALL ipslerr (3,'resget_opp_r3d', &
+      'Unable to performe an operation on this variable as it has', &
+      'a second and third dimension',vname_q)
+  ENDIF
+!-
+  IF (jjm < 1) THEN
+    CALL ipslerr (3,'resget_opp_r3d', &
+      'Please specify a second dimension which is the', &
+      'layer on which the operations are performed',vname_q)
+  ENDIF
+!-
+  IF (llm < 1) THEN
+    CALL ipslerr (3,'resget_opp_r3d', &
+      'Please specify a third dimension which is the', &
+      'layer on which the operations are performed',vname_q)
+  ENDIF
+!-
+! 1.0 Allocate the temporary buffer we need
+!     to put the variable in right dimension
+!-
+  siz1 = SIZE(var,1)
+  CALL rest_alloc (1,siz1,l_dbg,'restget_opp_r3d')
+  CALL rest_alloc (2,req_sz*jjm*llm,l_dbg,'restget_opp_r3d')
+!-
+! 2.0 Here we get the full variable from the restart file
+!-
+  list_dims = (/ ax_infs(fid,RESTART_DIM_X,1,1), ax_infs(fid,RESTART_DIM_Y,1,1), jjm, llm /)
+!-
+  CALL restget_real &
+ & (fid,vname_q,list_dims,itau,def_beha,buff_tmp2)
+!-
+! 4.0 Transfer the buffer obtained from the restart file
+!     into the variable the model expects
+!-
+  topp = MY_OPERATOR(1:MIN(LEN_TRIM(MY_OPERATOR),7))
+!-
+  IF (INDEX(indchfun,topp(:LEN_TRIM(topp))) > 0) THEN
+    scal = missing_val
+    var_sz = siz1
+    DO ll = 1, llm
+      DO jj = 1,jjm
+          ist = (((ll-1)*jjm) + (jj-1))*req_sz+1
+          CALL mathop (topp,req_sz,buff_tmp2(ist:ist+req_sz-1), &
+ &           missing_val,nbindex,ijndex,scal,var_sz,buff_tmp1)
+          var(:,jj,ll) = buff_tmp1(1:siz1)
+      ENDDO
+    ENDDO
+  ELSE
+    CALL ipslerr (3,'resget_opp_r3d', &
+      'The operation you wish to do on the variable for the ',&
+      'restart file is not allowed.',topp)
+  ENDIF
+!-----------------------------
+END SUBROUTINE restget_opp_r3d
+!===
+SUBROUTINE restget_opp_r4d &
+ & (fid,vname_q,iim,jjm,llm,mmm,itau,def_beha, &
+ &  var,MY_OPERATOR,nbindex,ijndex)
+!---------------------------------------------------------------------
+!- This subroutine serves as an interface to restget_real
+!-
+!- Should work as restput_opp_r2d but the other way around !
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+  INTEGER :: fid
+  CHARACTER(LEN=*) :: vname_q
+  INTEGER :: iim,jjm,llm,mmm,itau
+  LOGICAL def_beha
+  REAL :: var(:,:,:,:)
+  CHARACTER(LEN=*) :: MY_OPERATOR
+  INTEGER :: nbindex,ijndex(nbindex)
+!-
+  INTEGER :: mm,ll,jj,req_sz,ist,var_sz,siz1
+  REAL :: scal
+  CHARACTER(LEN=7) :: topp
+  LOGICAL :: l_dbg
+  INTEGER :: list_dims(5)
+!---------------------------------------------------------------------
+  CALL ipsldbg (old_status=l_dbg)
+!-
+! 0.0 What size should be the data in the file
+!-
+  req_sz = 1
+  IF (nbindex == iim) THEN
+    IF (ax_infs(fid,RESTART_DIM_X,1,1) > 0) req_sz = req_sz*ax_infs(fid,RESTART_DIM_X,1,1)
+    IF (ax_infs(fid,RESTART_DIM_Y,1,1) > 0) req_sz = req_sz*ax_infs(fid,RESTART_DIM_Y,1,1)
+  ELSE
+    CALL ipslerr (3,'restget_opp_r4d', &
+      'Unable to performe an operation on this variable as it has', &
+      'a second and third dimension',vname_q)
+  ENDIF
+!-
+! 1.0 Allocate the temporary buffer we need
+!     to put the variable in right dimension
+!-
+  siz1 = SIZE(var,1)
+  CALL rest_alloc (1,siz1,l_dbg,'restget_opp_r4d')
+  CALL rest_alloc (2,req_sz*jjm*llm*mmm,l_dbg,'restget_opp_r4d')
+!-
+! 2.0 Here we get the full variable from the restart file
+!-
+  list_dims = (/ ax_infs(fid,RESTART_DIM_X,1,1), ax_infs(fid,RESTART_DIM_Y,1,1), jjm, llm, mmm /)
+!-
+  CALL restget_real &
+ & (fid,vname_q,list_dims,itau,def_beha,buff_tmp2)
+!-
+! 4.0 Transfer the buffer obtained from the restart file
+!     into the variable the model expects
+!-
+  topp = MY_OPERATOR(1:MIN(LEN_TRIM(MY_OPERATOR),7))
+!-
+  IF (INDEX(indchfun,topp(:LEN_TRIM(topp))) > 0) THEN
+    scal = missing_val
+    var_sz = siz1
+    DO mm = 1, mmm
+      DO ll = 1, llm
+        DO jj = 1,jjm
+            ist = ((((mm-1)*llm + (ll-1))*jjm) + (jj-1))*req_sz+1
+            CALL mathop (topp,req_sz,buff_tmp2(ist:ist+req_sz-1), &
+ &             missing_val,nbindex,ijndex,scal,var_sz,buff_tmp1)
+            var(:,jj,ll,mm) = buff_tmp1(1:siz1)
+        ENDDO
+      ENDDO
+    ENDDO
+  ELSE
+    CALL ipslerr (3,'restget_opp_r4d', &
+      'The operation you wish to do on the variable for the ',&
+      'restart file is not allowed.',topp)
+  ENDIF
+!-----------------------------
+END SUBROUTINE restget_opp_r4d
+!===
+SUBROUTINE restget_opp_r5d &
+ & (fid,vname_q,iim,jjm,llm,mmm,nnm,itau,def_beha, &
+ &  var,MY_OPERATOR,nbindex,ijndex)
+!---------------------------------------------------------------------
+!- This subroutine serves as an interface to restget_real
+!-
+!- Should work as restput_opp_r2d but the other way around !
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+  INTEGER :: fid
+  CHARACTER(LEN=*) :: vname_q
+  INTEGER :: iim,jjm,llm,mmm,nnm,itau
+  LOGICAL def_beha
+  REAL :: var(:,:,:,:,:)
+  CHARACTER(LEN=*) :: MY_OPERATOR
+  INTEGER :: nbindex,ijndex(nbindex)
+!-
+  INTEGER :: mm,ll,jj,nn,req_sz,ist,var_sz,siz1
+  REAL :: scal
+  CHARACTER(LEN=7) :: topp
+  LOGICAL :: l_dbg
+  INTEGER :: list_dims(6)
+!---------------------------------------------------------------------
+  CALL ipsldbg (old_status=l_dbg)
+!-
+! 0.0 What size should be the data in the file
+!-
+  req_sz = 1
+  IF (nbindex == iim) THEN
+    IF (ax_infs(fid,RESTART_DIM_X,1,1) > 0) req_sz = req_sz*ax_infs(fid,RESTART_DIM_X,1,1)
+    IF (ax_infs(fid,RESTART_DIM_Y,1,1) > 0) req_sz = req_sz*ax_infs(fid,RESTART_DIM_Y,1,1)
+  ELSE
+    CALL ipslerr (3,'restget_opp_r5d', &
+      'Unable to performe an operation on this variable as it has', &
+      'a second and third dimension',vname_q)
+  ENDIF
+!-
+! 1.0 Allocate the temporary buffer we need
+!     to put the variable in right dimension
+!-
+  siz1 = SIZE(var,1)
+  CALL rest_alloc (1,siz1,l_dbg,'restget_opp_r5d')
+  CALL rest_alloc (2,req_sz*jjm*llm*mmm*nnm,l_dbg,'restget_opp_r5d')
+!-
+! 2.0 Here we get the full variable from the restart file
+!-
+  list_dims = (/ ax_infs(fid,RESTART_DIM_X,1,1), ax_infs(fid,RESTART_DIM_Y,1,1), jjm, llm, mmm, nnm /)
+!-
+  CALL restget_real &
+ & (fid,vname_q,list_dims,itau,def_beha,buff_tmp2)
+!-
+! 4.0 Transfer the buffer obtained from the restart file
+!     into the variable the model expects
+!-
+  topp = MY_OPERATOR(1:MIN(LEN_TRIM(MY_OPERATOR),7))
+!-
+  IF (INDEX(indchfun,topp(:LEN_TRIM(topp))) > 0) THEN
+    scal = missing_val
+    var_sz = siz1
+    DO nn=1, nnm
+      DO mm = 1, mmm
+        DO ll = 1, llm
+          DO jj = 1,jjm
+              ist = (((((nn-1)*mmm) + (mm-1)*llm + (ll-1))*jjm) + (jj-1))*req_sz+1
+              CALL mathop (topp,req_sz,buff_tmp2(ist:ist+req_sz-1), &
+ &               missing_val,nbindex,ijndex,scal,var_sz,buff_tmp1)
+              var(:,jj,ll,mm,nn) = buff_tmp1(1:siz1)
+          ENDDO
+        ENDDO
+      ENDDO
+    ENDDO
+  ELSE
+    CALL ipslerr (3,'restget_opp_r5d', &
+      'The operation you wish to do on the variable for the ',&
+      'restart file is not allowed.',topp)
+  ENDIF
+!-----------------------------
+END SUBROUTINE restget_opp_r5d
+!===
 SUBROUTINE restget_r1d &
  & (fid,vname_q,iim,jjm,llm,itau,def_beha,var)
@@ -1155,4 +1423,5 @@
   CHARACTER(LEN=70) :: str,str2
   LOGICAL :: l_dbg
+  INTEGER :: list_dims(2)
 !---------------------------------------------------------------------
   CALL ipsldbg (old_status=l_dbg)
@@ -1188,6 +1457,8 @@
   ENDIF
 !-
+  list_dims = (/ iim,jjm /)
+!-
   CALL restget_real &
- & (fid,vname_q,iim,jjm,llm,itau,def_beha,buff_tmp1)
+ & (fid,vname_q,list_dims,itau,def_beha,buff_tmp1)
 !-
 ! 4.0 Transfer the buffer obtained from the restart file
@@ -1218,4 +1489,5 @@
   CHARACTER(LEN=70) :: str,str2
   LOGICAL :: l_dbg
+  INTEGER :: list_dims(2)
 !---------------------------------------------------------------------
   CALL ipsldbg (old_status=l_dbg)
@@ -1253,6 +1525,8 @@
   ENDIF
 !-
+  list_dims = (/ iim,jjm /)
+!-
   CALL restget_real &
- & (fid,vname_q,iim,jjm,llm,itau,def_beha,buff_tmp1)
+ & (fid,vname_q,list_dims,itau,def_beha,buff_tmp1)
 !-
 ! 4.0 Transfer the buffer obtained from the restart file
@@ -1285,4 +1559,5 @@
   CHARACTER(LEN=70) :: str,str2
   LOGICAL :: l_dbg
+  INTEGER :: list_dims(3)
 !---------------------------------------------------------------------
   CALL ipsldbg (old_status=l_dbg)
@@ -1321,6 +1596,8 @@
   ENDIF
 !-
+  list_dims = (/ iim,jjm,llm /)
+!-
   CALL restget_real &
-    (fid,vname_q,iim,jjm,llm,itau,def_beha,buff_tmp1)
+    (fid,vname_q,list_dims,itau,def_beha,buff_tmp1)
 !-
 ! 4.0 Transfer the buffer obtained from the restart file
@@ -1340,5 +1617,5 @@
 !===
 SUBROUTINE restget_real &
-  (fid,vname_q,iim,jjm,llm,itau,def_beha,var)
+  (fid,vname_q,list_dims,itau,def_beha,var)
 !---------------------------------------------------------------------
 !- This subroutine is for getting a variable from the restart file.
@@ -1368,17 +1645,18 @@
   INTEGER :: fid
   CHARACTER(LEN=*) :: vname_q
-  INTEGER :: iim,jjm,llm,itau
+  INTEGER,DIMENSION(:),INTENT(in) :: list_dims
+  INTEGER :: itau
   LOGICAL :: def_beha
   REAL :: var(:)
 !-
-  INTEGER :: vid,vnb,ncfid,iret,index,it,ndim,ia
+  INTEGER :: vid,vnb,ncfid,iret,index,it,ndim,ia,itedim
   CHARACTER(LEN=70) str,str2
   CHARACTER(LEN=80) attname
-  INTEGER,DIMENSION(4) :: corner,edge
+  INTEGER,DIMENSION(RESTART_MAX_DIMS) :: corner,edge
   LOGICAL :: l_dbg
 !---------------------------------------------------------------------
   CALL ipsldbg (old_status=l_dbg)
 !---------------------------------------------------------------------
-  IF (l_dbg) WRITE(ipslout,*) 'RESTGET 0.0 : ',netcdf_name(fid,2),vname_q,iim,jjm,llm,itau,def_beha
+  IF (l_dbg) WRITE(ipslout,*) 'RESTGET 0.0 : ',netcdf_name(fid,2),vname_q,list_dims,itau,def_beha
 !-
   ncfid = netcdf_id(fid,1)
@@ -1390,4 +1668,19 @@
 !-
   IF (l_dbg) WRITE(ipslout,*) 'RESTGET 1.0 : ',vnb
+!-
+  IF (ANY(list_dims < 1)) THEN
+      CALL ipslerr (3,'restget', &
+ &      'All values in list_dims must be positive','','')
+  ENDIF
+  IF (SIZE(list_dims, 1) > RESTART_MAX_DIMS) THEN
+      CALL ipslerr (3,'restget', &
+ &      'Limit reached for dimensions','Please increase RESTART_MAX_DIMS to allow', &
+        ' more dimensions in single variable ')
+  ENDIF
+  IF (SIZE(list_dims, 1) == 0) THEN
+      CALL ipslerr (3,'restget', &
+ &      'No dimension is given','Make sure to pass a valid list_dims', &
+        ' ')
+  ENDIF
 !-
   IF (vnb < 0) THEN
@@ -1412,5 +1705,5 @@
       touched_in(fid,vnb) = .TRUE.
 !-----
-      CALL restdefv (fid,vname_q,iim,jjm,llm,.TRUE.)
+      CALL restdefv (fid,vname_q,list_dims,.TRUE.)
       IF (l_dbg) WRITE(ipslout,*) 'RESTGET 1.1 : ',vnb
 !-----
@@ -1457,39 +1750,16 @@
     str='Incorrect dimension for '//TRIM(vname_q)
     ndim = 0
-    IF (iim > 0) THEN
-      ndim = ndim+1
-      IF (vardims_in(fid,vnb,ndim) == iim) THEN
-        corner(ndim) = 1
-        edge(ndim) = iim
-      ELSE
-        WRITE (str2,'("Incompatibility for iim : ",I6,I6)') &
-             iim,vardims_in(fid,vnb,ndim)
-        CALL ipslerr (3,'restget',str,str2,' ')
-      ENDIF
-    ENDIF
-!---
-    IF (jjm > 0) THEN
-      ndim = ndim+1
-      IF (vardims_in(fid,vnb,ndim) == jjm) THEN
-        corner(ndim) = 1
-        edge(ndim) = jjm
-      ELSE
-        WRITE (str2,'("Incompatibility for jjm : ",I6,I6)') &
-             jjm,vardims_in(fid,vnb,ndim)
-        CALL ipslerr (3,'restget',str,str2,' ')
-      ENDIF
-    ENDIF
-!---
-    IF (llm > 0) THEN
-      ndim = ndim+1
-      IF (vardims_in(fid,vnb,ndim) == llm) THEN
-        corner(ndim) = 1
-        edge(ndim) = llm
-      ELSE
-        WRITE (str2,'("Incompatibility for llm : ",I6,I6)') &
-             llm,vardims_in(fid,vnb,ndim)
-        CALL ipslerr (3,'restget',str,str2,' ')
-      ENDIF
-    ENDIF
+!--
+    DO itedim=1, SIZE(list_dims,1)
+        ndim = ndim+1
+        IF (vardims_in(fid,vnb,ndim) == list_dims(itedim)) THEN
+          corner(ndim) = 1
+          edge(ndim) = list_dims(itedim)
+        ELSE
+          WRITE (str2,'("Incompatibility for I2 dimension : ",I6,I6,I6)') &
+               itedim,list_dims(itedim),vardims_in(fid,vnb,ndim)
+          CALL ipslerr (3,'restget',str,str2,' ')
+        ENDIF
+    ENDDO
 !---
 !-- Time
@@ -1497,5 +1767,4 @@
     ndim = ndim+1
     corner(ndim) = index
-!!????? edge(ndim) = index
     edge(ndim) = 1
 !---
@@ -1510,5 +1779,5 @@
  &      .AND.(netcdf_id(fid,2) > 0) ) THEN
 !-----
-      CALL restdefv (fid,vname_q,iim,jjm,llm,.FALSE.)
+      CALL restdefv (fid,vname_q,list_dims,.FALSE.)
 !-----
       DO ia = 1,varatt_in(fid,vnb)
@@ -1553,4 +1822,5 @@
   CHARACTER(LEN=7) :: topp
   LOGICAL :: l_dbg
+  INTEGER :: list_dims(2)
 !---------------------------------------------------------------------
   CALL ipsldbg (old_status=l_dbg)
@@ -1560,7 +1830,6 @@
   req_sz = 1
   IF ( nbindex == iim .AND. jjm <= 1 .AND. llm <= 1) THEN
-    IF (xax_infs(fid,1,1) > 0) req_sz = req_sz*xax_infs(fid,1,1)
-    IF (yax_infs(fid,1,1) > 0) req_sz = req_sz*yax_infs(fid,1,1)
-    IF (zax_infs(fid,1,1) > 0) req_sz = req_sz*zax_infs(fid,1,1)
+    IF (ax_infs(fid,RESTART_DIM_X,1,1) > 0) req_sz = req_sz*ax_infs(fid,RESTART_DIM_X,1,1)
+    IF (ax_infs(fid,RESTART_DIM_Y,1,1) > 0) req_sz = req_sz*ax_infs(fid,RESTART_DIM_Y,1,1)
   ELSE
     CALL ipslerr (3,'restput_opp_r1d', &
@@ -1594,7 +1863,8 @@
   ENDIF
 !-
+  list_dims = (/ ax_infs(fid,RESTART_DIM_X,1,1), ax_infs(fid,RESTART_DIM_Y,1,1) /)
+!-
   CALL restput_real &
- & (fid,vname_q,xax_infs(fid,1,1),yax_infs(fid,1,1), &
- &  zax_infs(fid,1,1),itau,buff_tmp2)
+ & (fid,vname_q,list_dims,itau,buff_tmp2)
 !-----------------------------
 END SUBROUTINE restput_opp_r1d
@@ -1629,4 +1899,5 @@
   CHARACTER(LEN=7) :: topp
   LOGICAL :: l_dbg
+  INTEGER :: list_dims(3)
 !---------------------------------------------------------------------
   CALL ipsldbg (old_status=l_dbg)
@@ -1636,6 +1907,6 @@
   req_sz = 1
   IF ( nbindex == iim .AND. llm <= 1) THEN
-    IF (xax_infs(fid,1,1) > 0) req_sz = req_sz*xax_infs(fid,1,1)
-    IF (yax_infs(fid,1,1) > 0) req_sz = req_sz*yax_infs(fid,1,1)
+    IF (ax_infs(fid,RESTART_DIM_X,1,1) > 0) req_sz = req_sz*ax_infs(fid,RESTART_DIM_X,1,1)
+    IF (ax_infs(fid,RESTART_DIM_Y,1,1) > 0) req_sz = req_sz*ax_infs(fid,RESTART_DIM_Y,1,1)
   ELSE
     CALL ipslerr (3,'restput_opp_r2d', &
@@ -1678,9 +1949,267 @@
   ENDIF
 !-
-  CALL restput_real &
- & (fid,vname_q,xax_infs(fid,1,1),yax_infs(fid,1,1), &
- &  jjm,itau,buff_tmp2)
+  list_dims = (/ ax_infs(fid,RESTART_DIM_X,1,1),ax_infs(fid,RESTART_DIM_Y,1,1),jjm /)
+!-
+  CALL restput_real (fid,vname_q, list_dims,itau,buff_tmp2)
 !-----------------------------
 END SUBROUTINE restput_opp_r2d
+!===
+SUBROUTINE restput_opp_r3d &
+ & (fid,vname_q,iim,jjm,llm,itau,var,MY_OPERATOR,nbindex,ijndex)
+!---------------------------------------------------------------------
+!- This subroutine is the interface to restput_real which allows
+!- to re-index data onto the original grid of the restart file.
+!- The logic we use is still fuzzy in my mind but that is probably
+!- only because I have not yet though through everything.
+!-
+!- In the case iim = nbindex it means that the user attempts
+!- to project the first dimension of the matrix back onto a 3D field
+!- where jjm will be the third dimension.
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+  INTEGER :: fid
+  CHARACTER(LEN=*) :: vname_q
+  INTEGER :: iim,jjm,llm,itau
+  REAL :: var(:,:,:)
+  CHARACTER(LEN=*) :: MY_OPERATOR
+  INTEGER :: nbindex,ijndex(nbindex)
+!-
+  INTEGER :: jj,ll,req_sz,ist,siz1
+  REAL :: scal
+  CHARACTER(LEN=7) :: topp
+  LOGICAL :: l_dbg
+  INTEGER :: list_dims(4)
+!---------------------------------------------------------------------
+  CALL ipsldbg (old_status=l_dbg)
+!-
+! 0.0 What size should be the data in the file
+!-
+  req_sz = 1
+  IF ( nbindex == iim) THEN
+    IF (ax_infs(fid,RESTART_DIM_X,1,1) > 0) req_sz = req_sz*ax_infs(fid,RESTART_DIM_X,1,1)
+    IF (ax_infs(fid,RESTART_DIM_Y,1,1) > 0) req_sz = req_sz*ax_infs(fid,RESTART_DIM_Y,1,1)
+  ELSE
+    CALL ipslerr (3,'restput_opp_r3d', &
+      'Unable to performe an operation on this variable as it has', &
+      'a second and third dimension',vname_q)
+  ENDIF
+!-
+  IF (jjm < 1 .OR. llm < 1) THEN
+    CALL ipslerr (3,'restput_opp_r3d', &
+      'Please specify a second dimension which is the', &
+      'layer on which the operations are performed',vname_q)
+  ENDIF
+!-
+! 1.0 Allocate the temporary buffer we need
+!     to put the variable in right dimension
+!-
+  siz1 = SIZE(var,1)
+  CALL rest_alloc (1,siz1,l_dbg,'restput_opp_r3d')
+  CALL rest_alloc (2,req_sz*jjm*llm,l_dbg,'restput_opp_r3d')
+!-
+! 2.0 We do the operation needed.
+!     It can only be a re-indexing operation.
+!     You would not want to change the values in a restart file or ?
+!-
+  topp = MY_OPERATOR(1:MIN(LEN_TRIM(MY_OPERATOR),7))
+!-
+  IF (INDEX(indchfun,topp(:LEN_TRIM(topp))) > 0) THEN
+    scal = missing_val
+    DO ll = 1,llm
+        DO jj = 1,jjm
+          buff_tmp1(1:siz1) = var(:,jj,ll)
+          ist = (((ll-1)*jjm) + (jj-1))*req_sz+1
+          CALL mathop &
+ &        (topp,siz1,buff_tmp1,missing_val,nbindex,ijndex, &
+ &           scal,req_sz,buff_tmp2(ist:ist+req_sz-1))
+        ENDDO
+    ENDDO
+  ELSE
+    CALL ipslerr (3,'restput_opp_r3d', &
+ &    'The operation you wish to do on the variable for the ', &
+ &    'restart file is not allowed.',topp)
+  ENDIF
+!-
+  list_dims = (/ ax_infs(fid,RESTART_DIM_X,1,1),ax_infs(fid,RESTART_DIM_Y,1,1),jjm,llm /)
+!-
+  CALL restput_real &
+ & (fid,vname_q,list_dims,itau,buff_tmp2)
+!-----------------------------
+END SUBROUTINE restput_opp_r3d
+!===
+SUBROUTINE restput_opp_r4d &
+ & (fid,vname_q,iim,jjm,llm,mmm,itau,var,MY_OPERATOR,nbindex,ijndex)
+!---------------------------------------------------------------------
+!- This subroutine is the interface to restput_real which allows
+!- to re-index data onto the original grid of the restart file.
+!- The logic we use is still fuzzy in my mind but that is probably
+!- only because I have not yet though through everything.
+!-
+!- In the case iim = nbindex it means that the user attempts
+!- to project the first dimension of the matrix back onto a 3D field
+!- where jjm will be the third dimension.
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+  INTEGER :: fid
+  CHARACTER(LEN=*) :: vname_q
+  INTEGER :: iim,jjm,llm,mmm,itau
+  REAL :: var(:,:,:,:)
+  CHARACTER(LEN=*) :: MY_OPERATOR
+  INTEGER :: nbindex,ijndex(nbindex)
+!-
+  INTEGER :: jj,ll,mm,req_sz,ist,siz1
+  REAL :: scal
+  CHARACTER(LEN=7) :: topp
+  LOGICAL :: l_dbg
+  INTEGER :: list_dims(5)
+!---------------------------------------------------------------------
+  CALL ipsldbg (old_status=l_dbg)
+!-
+! 0.0 What size should be the data in the file
+!-
+  req_sz = 1
+  IF ( nbindex == iim) THEN
+    IF (ax_infs(fid,RESTART_DIM_X,1,1) > 0) req_sz = req_sz*ax_infs(fid,RESTART_DIM_X,1,1)
+    IF (ax_infs(fid,RESTART_DIM_Y,1,1) > 0) req_sz = req_sz*ax_infs(fid,RESTART_DIM_Y,1,1)
+  ELSE
+    CALL ipslerr (3,'restput_opp_r4d', &
+      'Unable to performe an operation on this variable as it has', &
+      'a second and third dimension',vname_q)
+  ENDIF
+!-
+  IF (jjm < 1 .OR. llm < 1 .OR. mmm < 1) THEN
+    CALL ipslerr (3,'restput_opp_r4d', &
+      'Please specify a second dimension which is the', &
+      'layer on which the operations are performed',vname_q)
+  ENDIF
+!-
+! 1.0 Allocate the temporary buffer we need
+!     to put the variable in right dimension
+!-
+  siz1 = SIZE(var,1)
+  CALL rest_alloc (1,siz1,l_dbg,'restput_opp_r4d')
+  CALL rest_alloc (2,req_sz*jjm*llm*mmm,l_dbg,'restput_opp_r4d')
+!-
+! 2.0 We do the operation needed.
+!     It can only be a re-indexing operation.
+!     You would not want to change the values in a restart file or ?
+!-
+  topp = MY_OPERATOR(1:MIN(LEN_TRIM(MY_OPERATOR),7))
+!-
+  IF (INDEX(indchfun,topp(:LEN_TRIM(topp))) > 0) THEN
+    scal = missing_val
+    DO mm = 1,mmm
+      DO ll = 1,llm
+        DO jj = 1,jjm
+          buff_tmp1(1:siz1) = var(:,jj,ll,mm)
+          ist = ((((mm-1)*llm + (ll-1))*jjm) + (jj-1))*req_sz+1
+          CALL mathop &
+ &        (topp,siz1,buff_tmp1,missing_val,nbindex,ijndex, &
+ &           scal,req_sz,buff_tmp2(ist:ist+req_sz-1))
+        ENDDO
+      ENDDO
+    ENDDO
+  ELSE
+    CALL ipslerr (3,'restput_opp_r4d', &
+ &    'The operation you wish to do on the variable for the ', &
+ &    'restart file is not allowed.',topp)
+  ENDIF
+!-
+  list_dims = (/ ax_infs(fid,RESTART_DIM_X,1,1),ax_infs(fid,RESTART_DIM_Y,1,1),jjm,llm,mmm /)
+!-
+  CALL restput_real &
+ & (fid,vname_q,list_dims,itau,buff_tmp2)
+!-----------------------------
+END SUBROUTINE restput_opp_r4d
+!===
+SUBROUTINE restput_opp_r5d &
+ & (fid,vname_q,iim,jjm,llm,mmm,nnm,itau,var,MY_OPERATOR,nbindex,ijndex)
+!---------------------------------------------------------------------
+!- This subroutine is the interface to restput_real which allows
+!- to re-index data onto the original grid of the restart file.
+!- The logic we use is still fuzzy in my mind but that is probably
+!- only because I have not yet though through everything.
+!-
+!- In the case iim = nbindex it means that the user attempts
+!- to project the first dimension of the matrix back onto a 3D field
+!- where jjm will be the third dimension.
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+  INTEGER :: fid
+  CHARACTER(LEN=*) :: vname_q
+  INTEGER :: iim,jjm,llm,mmm,nnm,itau
+  REAL :: var(:,:,:,:,:)
+  CHARACTER(LEN=*) :: MY_OPERATOR
+  INTEGER :: nbindex,ijndex(nbindex)
+!-
+  INTEGER :: jj,ll,mm,nn,req_sz,ist,siz1
+  REAL :: scal
+  CHARACTER(LEN=7) :: topp
+  LOGICAL :: l_dbg
+  INTEGER :: list_dims(6)
+!---------------------------------------------------------------------
+  CALL ipsldbg (old_status=l_dbg)
+!-
+! 0.0 What size should be the data in the file
+!-
+  req_sz = 1
+  IF ( nbindex == iim) THEN
+    IF (ax_infs(fid,RESTART_DIM_X,1,1) > 0) req_sz = req_sz*ax_infs(fid,RESTART_DIM_X,1,1)
+    IF (ax_infs(fid,RESTART_DIM_Y,1,1) > 0) req_sz = req_sz*ax_infs(fid,RESTART_DIM_Y,1,1)
+  ELSE
+    CALL ipslerr (3,'restput_opp_r5d', &
+      'Unable to performe an operation on this variable as it has', &
+      'a second and third dimension',vname_q)
+  ENDIF
+!-
+  IF (jjm < 1 .OR. llm < 1 .OR. mmm < 1 .OR. nnm < 1) THEN
+    CALL ipslerr (3,'restput_opp_r5d', &
+      'Please make sure all dimenensions are at least 1', &
+      '',vname_q)
+  ENDIF
+!-
+! 1.0 Allocate the temporary buffer we need
+!     to put the variable in right dimension
+!-
+  siz1 = SIZE(var,1)
+  CALL rest_alloc (1,siz1,l_dbg,'restput_opp_r5d')
+  CALL rest_alloc (2,req_sz*jjm*llm*mmm*nnm,l_dbg,'restput_opp_r5d')
+!-
+! 2.0 We do the operation needed.
+!     It can only be a re-indexing operation.
+!     You would not want to change the values in a restart file or ?
+!-
+  topp = MY_OPERATOR(1:MIN(LEN_TRIM(MY_OPERATOR),7))
+!-
+  IF (INDEX(indchfun,topp(:LEN_TRIM(topp))) > 0) THEN
+    scal = missing_val
+    DO nn = 1, nnm
+      DO mm = 1,mmm
+        DO ll = 1,llm
+          DO jj = 1,jjm
+            buff_tmp1(1:siz1) = var(:,jj,ll,mm,nn)
+            ist = (((((nn-1) * mmm) + (mm-1)*llm + (ll-1))*jjm) + (jj-1))*req_sz+1
+            CALL mathop &
+ &          (topp,siz1,buff_tmp1,missing_val,nbindex,ijndex, &
+ &             scal,req_sz,buff_tmp2(ist:ist+req_sz-1))
+          ENDDO
+        ENDDO
+      ENDDO
+    ENDDO
+  ELSE
+    CALL ipslerr (3,'restput_opp_r5d', &
+ &    'The operation you wish to do on the variable for the ', &
+ &    'restart file is not allowed.',topp)
+  ENDIF
+!-
+  list_dims = (/ ax_infs(fid,RESTART_DIM_X,1,1),ax_infs(fid,RESTART_DIM_Y,1,1),jjm,llm,mmm,nnm /)
+!-
+  CALL restput_real &
+ & (fid,vname_q,list_dims,itau,buff_tmp2)
+!-----------------------------
+END SUBROUTINE restput_opp_r5d
 !===
 SUBROUTINE restput_r1d (fid,vname_q,iim,jjm,llm,itau,var)
@@ -1698,4 +2227,5 @@
   CHARACTER(LEN=70) :: str,str2
   LOGICAL :: l_dbg
+  INTEGER :: list_dims(1)
 !---------------------------------------------------------------------
   CALL ipsldbg (old_status=l_dbg)
@@ -1704,4 +2234,5 @@
 !     to put the variable in right dimension
 !-
+  list_dims(1) = iim
   siz1 = SIZE(var)
   var_sz = siz1
@@ -1738,5 +2269,5 @@
   ENDDO
 !-
-  CALL restput_real (fid,vname_q,iim,jjm,llm,itau,buff_tmp1)
+  CALL restput_real (fid,vname_q,list_dims,itau,buff_tmp1)
 !-------------------------
 END SUBROUTINE restput_r1d
@@ -1756,4 +2287,5 @@
   CHARACTER(LEN=70) :: str,str2
   LOGICAL :: l_dbg
+  INTEGER :: list_dims(2)
 !---------------------------------------------------------------------
   CALL ipsldbg (old_status=l_dbg)
@@ -1762,4 +2294,5 @@
 !     to put the variable in right dimension
 !-
+  list_dims = (/ iim,jjm /)
   siz1 = SIZE(var,1)
   siz2 = SIZE(var,2)
@@ -1798,5 +2331,5 @@
   ENDDO
 !-
-  CALL restput_real(fid,vname_q,iim,jjm,llm,itau,buff_tmp1)
+  CALL restput_real(fid,vname_q,list_dims,itau,buff_tmp1)
 !-------------------------
 END SUBROUTINE restput_r2d
@@ -1816,4 +2349,5 @@
   CHARACTER(LEN=70) :: str,str2
   LOGICAL :: l_dbg
+  INTEGER :: list_dims(3)
 !---------------------------------------------------------------------
   CALL ipsldbg (old_status=l_dbg)
@@ -1822,4 +2356,5 @@
 !     to put the variable in right dimension
 !-
+  list_dims = (/ iim, jjm, llm /)
   siz1 = SIZE(var,1)
   siz2 = SIZE(var,2)
@@ -1862,9 +2397,9 @@
   ENDDO
 !-
-  CALL restput_real(fid,vname_q,iim,jjm,llm,itau,buff_tmp1)
+  CALL restput_real(fid,vname_q,list_dims,itau,buff_tmp1)
 !-------------------------
 END SUBROUTINE restput_r3d
 !===
-SUBROUTINE restput_real (fid,vname_q,iim,jjm,llm,itau,var)
+SUBROUTINE restput_real (fid,vname_q,list_dims,itau,var)
 !---------------------------------------------------------------------
 !- This subroutine will put a variable into the restart file.
@@ -1883,5 +2418,5 @@
 !- fid         : Id of the file in which we will write the variable
 !- vname_q     : Name of the variable to be written
-!- iim,jjm,llm : Size in 3D of the variable
+!- list_dim    : Size of the variable (1 to RESTART_MAX_DIMS)
 !- itau        : Time step at which the variable is written
 !- var         : Variable
@@ -1894,12 +2429,13 @@
 !-
   CHARACTER(LEN=*) :: vname_q
-  INTEGER :: fid,iim,jjm,llm,itau
+  INTEGER :: fid,itau
+  INTEGER,DIMENSION(:),INTENT(in) :: list_dims
   REAL :: var(:)
 !-
-  INTEGER :: iret,vid,ncid,iv,vnb
+  INTEGER :: iret,vid,ncid,iv,vnb,itedim
   INTEGER :: ierr
   REAL :: secsince,one_day,one_year
   INTEGER :: ndims
-  INTEGER,DIMENSION(4) :: corner,edge
+  INTEGER,DIMENSION(RESTART_MAX_DIMS) :: corner,edge
   LOGICAL :: l_dbg
 !---------------------------------------------------------------------
@@ -1916,5 +2452,5 @@
 !-
 ! 0.0 show arguments
-  IF (l_dbg) WRITE(ipslout,*) 'RESTPUT 0.0 : ',netcdf_name(fid,2),vname_q,iim,jjm,llm,itau
+  IF (l_dbg) WRITE(ipslout,*) 'RESTPUT 0.0 : ',netcdf_name(fid,2),vname_q,list_dims,itau
 !-
 ! 1.0 Check if the variable is already present
@@ -1932,5 +2468,5 @@
 !-
   IF (vnb <= 0) THEN
-    CALL restdefv (fid,vname_q,iim,jjm,llm,.TRUE.)
+    CALL restdefv (fid,vname_q,list_dims,.TRUE.)
     vnb = nbvar_out(fid)
   ENDIF
@@ -2005,19 +2541,10 @@
 !-
   ndims = 0
-  IF (iim > 0) THEN
+  DO itedim=1, SIZE(list_dims,1)
     ndims = ndims+1
     corner(ndims) = 1
-    edge(ndims) = iim
-  ENDIF
-  IF (jjm > 0) THEN
-    ndims = ndims+1
-    corner(ndims) = 1
-    edge(ndims) = jjm
-  ENDIF
-  IF (llm > 0) THEN
-    ndims = ndims+1
-    corner(ndims) = 1
-    edge(ndims) = llm
-  ENDIF
+    edge(ndims) = list_dims(itedim)
+  ENDDO
+
   ndims = ndims+1
   corner(ndims) = tstp_out(fid)
@@ -2028,5 +2555,5 @@
 !-
   IF (iret /= NF90_NOERR) THEN
-    CALL ipslerr (2,'restput_real',NF90_STRERROR(iret), &
+    CALL ipslerr (3,'restput_real',NF90_STRERROR(iret), &
  &    'Bug in restput.',&
  &    'Please, verify compatibility between get and put commands.')
@@ -2039,5 +2566,5 @@
 END  SUBROUTINE restput_real
 !===
-SUBROUTINE restdefv (fid,varname,iim,jjm,llm,write_att)
+SUBROUTINE restdefv (fid,varname,list_dims,write_att)
 !---------------------------------------------------------------------
 ! This subroutine adds a variable to the output file.
@@ -2048,13 +2575,34 @@
   INTEGER ::fid
   CHARACTER(LEN=*) :: varname
-  INTEGER :: iim,jjm,llm
+! define its dimension size, this must respect same order as RESTART_DIMS_NAMES
+! x,y,z,....
+  INTEGER, DIMENSION(:),INTENT(in) :: list_dims  ! Dimensions input
   LOGICAL :: write_att
 !-
-  INTEGER :: dims(4),ic,xloc,ndim,ncfid
-  INTEGER :: iret,ax_id
+  INTEGER :: dims(RESTART_MAX_DIMS),ic,xloc,ndim,ncfid
+  INTEGER :: iret,ax_id,ite
   CHARACTER(LEN=3) :: str
   LOGICAL :: l_dbg
 !---------------------------------------------------------------------
   CALL ipsldbg (old_status=l_dbg)
+!- Check consistency
+!-
+  IF (ANY(list_dims < 1)) THEN
+    CALL ipslerr (3,'restdefv', &
+      'All defined dimension sizes must be bigger than 0', &
+      'Please check list_dims argument',' ')
+  ENDIF
+!-
+  IF (SIZE(list_dims,1) == 0) THEN
+    CALL ipslerr (3,'restdefv', &
+      'You must add at least one dimension', &
+      'Please check list_dims argument',' ')
+  ENDIF
+!-
+  IF (SIZE(list_dims,1) > RESTART_MAX_DIMS) THEN
+    CALL ipslerr (3,'restdefv', &
+      'Maximum limit of dimensions in 1 variable reached', &
+      'Please increase RESTART_MAX_DIM value',' ')
+  ENDIF
 !-
   ncfid = netcdf_id(fid,2)
@@ -2084,66 +2632,26 @@
   ndim = 0
 !-
-! 1.1 Work on x
-!-
-  IF (iim > 0) THEN
-    ndim = ndim+1
-    xloc = 0
-    DO ic=1,xax_nb(fid)
-      IF (xax_infs(fid,ic,1) == iim) xloc = ic
-    ENDDO
-!---
-    IF (xloc > 0) THEN
-      dims(ndim) = xax_infs(fid,xloc,2)
-    ELSE
-      str='x_'//CHAR(96+xax_nb(fid))
-      iret = NF90_DEF_DIM(ncfid,str,iim,ax_id)
-      xax_nb(fid) = xax_nb(fid)+1
-      xax_infs(fid,xax_nb(fid),1) = iim
-      xax_infs(fid,xax_nb(fid),2) = ax_id
-      dims(ndim) = ax_id
-    ENDIF
-  ENDIF
-!-
-! 1.2 Work on y
-!-
-  IF (jjm > 0) THEN
-    ndim = ndim+1
-    xloc = 0
-    DO ic=1,yax_nb(fid)
-      IF (yax_infs(fid,ic,1) == jjm) xloc = ic
-    ENDDO
-!---
-    IF (xloc > 0) THEN
-      dims(ndim) = yax_infs(fid,xloc,2)
-    ELSE
-      str='y_'//CHAR(96+yax_nb(fid))
-      iret = NF90_DEF_DIM(ncfid,str,jjm,ax_id)
-      yax_nb(fid) = yax_nb(fid)+1
-      yax_infs(fid,yax_nb(fid),1) = jjm
-      yax_infs(fid,yax_nb(fid),2) = ax_id
-      dims(ndim) = ax_id
-    ENDIF
-  ENDIF
-!-
-! 1.3 Work on z
-!-
-  IF (llm > 0) THEN
-    ndim = ndim+1
-    xloc = 0
-    DO ic=1,zax_nb(fid)
-      IF (zax_infs(fid,ic,1) == llm) xloc = ic
-    ENDDO
-!---
-    IF (xloc > 0) THEN
-      dims(ndim) = zax_infs(fid,xloc,2)
-    ELSE
-      str='z_'//CHAR(96+zax_nb(fid))
-      iret = NF90_DEF_DIM(ncfid,str,llm,ax_id)
-      zax_nb(fid) = zax_nb(fid)+1
-      zax_infs(fid,zax_nb(fid),1) = llm
-      zax_infs(fid,zax_nb(fid),2) = ax_id
-      dims(ndim) = ax_id
-    ENDIF
-  ENDIF
+! 1.1 Work on each dimension
+!-
+  DO ite=1, SIZE(list_dims, 1)  
+    IF (list_dims(ite) > 0) THEN
+      ndim = ndim+1
+      xloc = 0
+      DO ic=1,ax_nb(fid, ite)
+        IF (ax_infs(fid,ite,ic,1) == list_dims(ite)) xloc = ic
+      ENDDO
+!---
+      IF (xloc > 0) THEN
+        dims(ndim) = ax_infs(fid,ite,xloc,2)
+      ELSE
+        str=RESTART_DIMS_NAMES(ite)//'_'//CHAR(96+ax_nb(fid,ite))
+        iret = NF90_DEF_DIM(ncfid,str,list_dims(ite),ax_id)
+        ax_nb(fid,ite) = ax_nb(fid,ite)+1
+        ax_infs(fid,ite,ax_nb(fid,ite),1) = list_dims(ite)
+        ax_infs(fid,ite,ax_nb(fid,ite),2) = ax_id
+        dims(ndim) = ax_id
+      ENDIF
+    ENDIF
+  ENDDO
 !-
 ! 1.4  Time needs to be added
@@ -2183,4 +2691,7 @@
     iret = NF90_PUT_ATT(ncfid,varid_out(fid,nbvar_out(fid)), &
  &                      'missing_value',REAL(missing_val,KIND=4))
+!---
+    iret = NF90_PUT_ATT(ncfid,varid_out(fid,nbvar_out(fid)), &
+ &                      '_FillValue',REAL(val_exp,r_std))
 !---
     IF (itau_out(fid) >= 0) THEN
