Index: IOIPSL/trunk/src/histcom.f90
===================================================================
--- IOIPSL/trunk/src/histcom.f90	(revision 957)
+++ IOIPSL/trunk/src/histcom.f90	(revision 962)
@@ -37,9 +37,9 @@
 !-
   INTERFACE histbeg
-    MODULE PROCEDURE histbeg_totreg,histbeg_regular,histbeg_irregular
+    MODULE PROCEDURE histb_reg1d,histb_reg2d,histb_irreg
   END INTERFACE
 !-
   INTERFACE histhori
-    MODULE PROCEDURE histhori_regular,histhori_irregular
+    MODULE PROCEDURE histh_reg1d,histh_reg2d,histh_irreg
   END INTERFACE
 !-
@@ -129,5 +129,5 @@
   INTEGER :: tid,xid,yid
 !-General definitions in the NETCDF file
-  INTEGER,DIMENSION(2) :: full_size=0,slab_ori,slab_sz
+  INTEGER,DIMENSION(2) :: full_size=0,slab_ori,slab_siz
 !-The horizontal axes
   CHARACTER(LEN=25),DIMENSION(nb_hax_max,2) :: hax_name
@@ -145,8 +145,4 @@
 TYPE(T_D_F),DIMENSION(nb_files_max),SAVE :: W_F
 !-
-! Counter of elements
-!-
-  INTEGER,SAVE :: nb_files=0
-!-
 ! A list of functions which require special action
 ! (Needs to be updated when functions are added
@@ -162,55 +158,10 @@
 !===
 !-
-SUBROUTINE histbeg_totreg                 &
- & (pfilename,pim,plon,pjm,plat,          &
- &  par_orix,par_szx,par_oriy,par_szy,    &
- &  pitau0,pdate0,pdeltat,phoriid,idf,domain_id,nb_bits)
-!---------------------------------------------------------------------
-!- This is just an interface for histbeg_regular in case when
-!- the user provides plon and plat as vectors.
-!- Obviously this can only be used for very regular grids.
-!-
-!- INPUT
-!-
-!- pfilename : Name of the netcdf file to be created
-!- pim       : Size of arrays in longitude direction
-!- plon      : Coordinates of points in longitude
-!- pjm       : Size of arrays in latitude direction
-!- plat      : Coordinates of points in latitude
-!-
-!- The next 4 arguments allow to define a horizontal zoom
-!- for this file. It is assumed that all variables to come
-!- have the same index space. This can not be assumed for
-!- the z axis and thus we define the zoom in histdef.
-!-
-!- par_orix  : Origin of the slab of data within the X axis (pim)
-!- par_szx   : Size of the slab of data in X
-!- par_oriy  : Origin of the slab of data within the Y axis (pjm)
-!- par_szy   : Size of the slab of data in Y
-!-
-!- pitau0    : time step at which the history tape starts
-!- pdate0    : The Julian date at which the itau was equal to 0
-!- pdeltat   : Time step, in seconds, of the counter itau
-!-             used in histwrite for instance
-!-
-!- OUTPUT
-!-
-!- phoriid   : ID of the horizontal grid
-!- idf       : ID of the netcdf file
-!-
-!- Optional INPUT arguments
-!-
-!- domain_id       : Domain identifier
-!-
-!- TO DO
-!-
-!- This package should be written in f90
-!- and use the following features :
-!-    - structures for the meta-data of the files and variables
-!-    - memory allocation as needed
-!-    - Pointers
-!-
-!- VERSION
-!-
+SUBROUTINE histb_reg1d                 &
+ & (pfilename,pim,plon,pjm,plat,       &
+ &  par_orix,par_szx,par_oriy,par_szy, &
+ &  pitau0,pdate0,pdeltat,phoriid,idf,domain_id,mode)
+!---------------------------------------------------------------------
+!- histbeg for 1D regular horizontal coordinates (see histb_all)
 !---------------------------------------------------------------------
   IMPLICIT NONE
@@ -224,82 +175,21 @@
   REAL,INTENT(IN) :: pdate0,pdeltat
   INTEGER,INTENT(OUT) :: idf,phoriid
-  INTEGER,INTENT(IN),OPTIONAL :: domain_id,nb_bits
-!-
-  REAL,ALLOCATABLE,DIMENSION(:,:) :: lon_tmp,lat_tmp
-  LOGICAL :: l_dbg
-!---------------------------------------------------------------------
-  CALL ipsldbg (old_status=l_dbg)
-!-
-  IF (l_dbg) WRITE(*,*) "histbeg_totreg"
-!-
-  ALLOCATE(lon_tmp(pim,pjm),lat_tmp(pim,pjm))
-!-
-  lon_tmp(:,:) = SPREAD(plon(:),2,pjm)
-  lat_tmp(:,:) = SPREAD(plat(:),1,pim)
-!-
-  CALL histbeg_regular &
- &  (pfilename,pim,lon_tmp,pjm,lat_tmp, &
- &   par_orix,par_szx,par_oriy,par_szy, &
- &   pitau0,pdate0,pdeltat,phoriid,idf, &
- &   .TRUE.,domain_id,nb_bits)
-!-
-  DEALLOCATE(lon_tmp,lat_tmp)
-!----------------------------
-END SUBROUTINE histbeg_totreg
+  INTEGER,INTENT(IN),OPTIONAL :: domain_id
+  CHARACTER(LEN=*),OPTIONAL,INTENT(IN) :: mode
+!---------------------------------------------------------------------
+  CALL histb_all &
+ & (1,pfilename,pim,pjm,pitau0,pdate0,pdeltat,phoriid,idf, &
+ &  x_1d=plon,y_1d=plat, &
+ &  k_orx=par_orix,k_szx=par_szx,k_ory=par_oriy,k_szy=par_szy, &
+ &  domain_id=domain_id,mode=mode)
+!-------------------------
+END SUBROUTINE histb_reg1d
 !===
-SUBROUTINE histbeg_regular &
- & (pfilename,pim,plon,pjm,plat,           &
- &  par_orix,par_szx,par_oriy,par_szy,     &
- &  pitau0,pdate0,pdeltat,phoriid,idf, &
- &  opt_rectilinear,domain_id,nb_bits)
-!---------------------------------------------------------------------
-!- This subroutine initializes a netcdf file and returns the ID.
-!- It will set up the geographical space on which the data will be
-!- stored and offers the possibility of seting a zoom.
-!- It also gets the global parameters into the I/O subsystem.
-!-
-!- INPUT
-!-
-!- pfilename : Name of the netcdf file to be created
-!- pim       : Size of arrays in longitude direction
-!- plon      : Coordinates of points in longitude
-!- pjm       : Size of arrays in latitude direction
-!- plat      : Coordinates of points in latitude
-!-
-!- The next 4 arguments allow to define a horizontal zoom
-!- for this file. It is assumed that all variables to come
-!- have the same index space. This can not be assumed for
-!- the z axis and thus we define the zoom in histdef.
-!-
-!- par_orix  : Origin of the slab of data within the X axis (pim)
-!- par_szx   : Size of the slab of data in X
-!- par_oriy  : Origin of the slab of data within the Y axis (pjm)
-!- par_szy   : Size of the slab of data in Y
-!-
-!- pitau0    : time step at which the history tape starts
-!- pdate0    : The Julian date at which the itau was equal to 0
-!- pdeltat   : Time step, in seconds, of the counter itau
-!-             used in histwrite for instance
-!-
-!- OUTPUT
-!-
-!- phoriid   : ID of the horizontal grid
-!- idf       : ID of the netcdf file
-!-
-!- Optional INPUT arguments
-!-
-!- opt_rectilinear : If true we know the grid is rectilinear
-!- domain_id       : Domain identifier
-!-
-!- TO DO
-!-
-!- This package should be written in F90 and use the following
-!- feature :
-!-   - structures for the meta-data of the files and variables
-!-   - memory allocation as needed
-!-   - Pointers
-!-
-!- VERSION
-!-
+SUBROUTINE histb_reg2d &
+ & (pfilename,pim,plon,pjm,plat,       &
+ &  par_orix,par_szx,par_oriy,par_szy, &
+ &  pitau0,pdate0,pdeltat,phoriid,idf,domain_id,mode)
+!---------------------------------------------------------------------
+!- histbeg for 2D regular horizontal coordinates (see histb_all)
 !---------------------------------------------------------------------
   IMPLICIT NONE
@@ -312,162 +202,20 @@
   REAL,INTENT(IN) :: pdate0,pdeltat
   INTEGER,INTENT(OUT) :: idf,phoriid
-  LOGICAL,INTENT(IN),OPTIONAL :: opt_rectilinear
-  INTEGER,INTENT(IN),OPTIONAL :: domain_id,nb_bits
-!-
-  INTEGER :: nfid,iret,m_c
-  CHARACTER(LEN=120) :: file
-  CHARACTER(LEN=30) :: timenow
-  LOGICAL :: rectilinear
-  LOGICAL :: l_dbg
-!---------------------------------------------------------------------
-  CALL ipsldbg (old_status=l_dbg)
-!-
-  IF (l_dbg) WRITE(*,*) "histbeg_regular 0.0"
-!-
-  nb_files = nb_files+1
-  IF (nb_files > nb_files_max) THEN
-    CALL ipslerr (3,"histbeg", &
-   &  'Table of files too small. You should increase nb_files_max', &
-   &  'in histcom.f90 in order to accomodate all these files',' ')
-  ENDIF
-  idf = nb_files
-!-
-! 1.0 Transfering into the common for future use
-!-
-  IF (l_dbg) WRITE(*,*) "histbeg_regular 1.0"
-!-
-  W_F(idf)%itau0  = pitau0
-  W_F(idf)%date0  = pdate0
-  W_F(idf)%deltat = pdeltat
-!-
-  IF (PRESENT(opt_rectilinear)) THEN
-    rectilinear = opt_rectilinear
-  ELSE
-    rectilinear = .FALSE.
-  ENDIF
-!-
-! 2.0 Initializes all variables for this file
-!-
-  IF (l_dbg) WRITE(*,*) "histbeg_regular 2.0"
-!-
-  W_F(idf)%n_var = 0
-  W_F(idf)%n_tax = 0
-  W_F(idf)%n_hax = 0
-  W_F(idf)%n_zax = 0
-!-
-  W_F(idf)%slab_ori(1:2) = (/ par_orix,par_oriy /)
-  W_F(idf)%slab_sz(1:2)  = (/ par_szx, par_szy  /)
-!-
-! 3.0 Opening netcdf file and defining dimensions
-!-
-  IF (l_dbg) WRITE(*,*) "histbeg_regular 3.0"
-!-
-! Add DOMAIN number and ".nc" suffix in file name if needed
-!-
-  file = pfilename
-  CALL flio_dom_file (file,domain_id)
-!-
-! Check the mode
-  IF (PRESENT(nb_bits)) THEN
-    SELECT CASE (nb_bits)
-    CASE(32)
-      m_c = NF90_CLOBBER
-    CASE(64)
-      m_c = IOR(NF90_CLOBBER,NF90_64BIT_OFFSET)
-    CASE DEFAULT
-      CALL ipslerr (3,"histbeg", &
- &      'Invalid argument nb_bits for file :',TRIM(file), &
- &      'Supported values are 32 or 64')
-    END SELECT
-  ELSE
-    m_c = IOR(NF90_CLOBBER,NF90_64BIT_OFFSET)
-  ENDIF
-!-
-! Create file
-  iret = NF90_CREATE (file,m_c,nfid)
-!-
-  IF (rectilinear) THEN
-    iret = NF90_DEF_DIM (nfid,'lon',par_szx,W_F(idf)%xid)
-    iret = NF90_DEF_DIM (nfid,'lat',par_szy,W_F(idf)%yid)
-  ELSE
-    iret = NF90_DEF_DIM (nfid,'x',par_szx,W_F(idf)%xid)
-    iret = NF90_DEF_DIM (nfid,'y',par_szy,W_F(idf)%yid)
-  ENDIF
-!-
-! 4.0 Declaring the geographical coordinates and other attributes
-!-
-  IF (l_dbg) WRITE(*,*) "histbeg_regular 4.0"
-!-
-! 4.3 Global attributes
-!-
-  iret = NF90_PUT_ATT (nfid,NF90_GLOBAL,'Conventions','CF-1.1')
-  iret = NF90_PUT_ATT (nfid,NF90_GLOBAL,'file_name',TRIM(file))
-  iret = NF90_PUT_ATT (nfid,NF90_GLOBAL,'production',TRIM(model_name))
-  lock_modname = .TRUE.
-  CALL ioget_timestamp (timenow)
-  iret = NF90_PUT_ATT (nfid,NF90_GLOBAL,'TimeStamp',TRIM(timenow))
-!-
-! 5.0 Saving some important information on this file in the common
-!-
-  IF (l_dbg) WRITE(*,*) "histbeg_regular 5.0"
-!-
-  IF (PRESENT(domain_id)) THEN
-    W_F(idf)%dom_id_svg = domain_id
-  ENDIF
-  W_F(idf)%ncfid = nfid
-  W_F(idf)%full_size(1:2) = (/ pim,pjm /)
-!-
-! 6.0 storing the geographical coordinates
-!-
-  W_F(idf)%regular=.TRUE.
-!-
-  CALL histhori_regular (idf,pim,plon,pjm,plat, &
- &  ' ','Default grid',phoriid,rectilinear)
-!-----------------------------
-END SUBROUTINE histbeg_regular
+  INTEGER,INTENT(IN),OPTIONAL :: domain_id
+  CHARACTER(LEN=*),OPTIONAL,INTENT(IN) :: mode
+!---------------------------------------------------------------------
+  CALL histb_all &
+ & (2,pfilename,pim,pjm,pitau0,pdate0,pdeltat,phoriid,idf,  &
+ &  x_2d=plon,y_2d=plat, &
+ &  k_orx=par_orix,k_szx=par_szx,k_ory=par_oriy,k_szy=par_szy,    &
+ &  domain_id=domain_id,mode=mode)
+!-------------------------
+END SUBROUTINE histb_reg2d
 !===
-SUBROUTINE histbeg_irregular &
- &  (pfilename,pim,plon,plon_bounds,plat,plat_bounds,   &
- &   pitau0,pdate0,pdeltat,phoriid,idf,domain_id,nb_bits)
-!---------------------------------------------------------------------
-!- This subroutine initializes a netcdf file and returns the ID.
-!- This version is for totaly irregular grids. In this case all
-!- all the data comes in as vectors and for the grid we have
-!- the coordinates of the 4 corners.
-!- It also gets the global parameters into the I/O subsystem.
-!-
-!- INPUT
-!-
-!- pfilename   : Name of the netcdf file to be created
-!- pim         : Size of arrays in longitude direction
-!- plon        : Coordinates of points in longitude
-!- plon_bounds : The 2 corners of the grid in longitude
-!- plat        : Coordinates of points in latitude
-!- plat_bounds : The 2 corners of the grid in latitude
-!-
-!- pitau0    : time step at which the history tape starts
-!- pdate0    : The Julian date at which the itau was equal to 0
-!- pdeltat   : Time step, in seconds, of the counter itau
-!-             used in histwrite for instance
-!-
-!- OUTPUT
-!-
-!- phoriid   : ID of the horizontal grid
-!- idf       : ID of the netcdf file
-!-
-!- Optional INPUT arguments
-!-
-!- domain_id       : Domain identifier
-!-
-!- TO DO
-!-
-!- This package should be written in F90 and use the following
-!- feature :
-!-   - structures for the meta-data of the files and variables
-!-   - memory allocation as needed
-!-   - Pointers
-!-
-!- VERSION
-!-
+SUBROUTINE histb_irreg &
+ &  (pfilename,pim,plon,plon_bounds,plat,plat_bounds, &
+ &   pitau0,pdate0,pdeltat,phoriid,idf,domain_id,mode)
+!---------------------------------------------------------------------
+!- histbeg for irregular horizontal coordinates (see histb_all)
 !---------------------------------------------------------------------
   IMPLICIT NONE
@@ -480,26 +228,133 @@
   REAL,INTENT(IN) :: pdate0,pdeltat
   INTEGER,INTENT(OUT) :: idf,phoriid
-  INTEGER,INTENT(IN),OPTIONAL :: domain_id,nb_bits
+  INTEGER,INTENT(IN),OPTIONAL :: domain_id
+  CHARACTER(LEN=*),OPTIONAL,INTENT(IN) :: mode
+!---------------------------------------------------------------------
+  CALL histb_all &
+ & (3,pfilename,pim,pim,pitau0,pdate0,pdeltat,phoriid,idf,  &
+ &  x_1d=plon,y_1d=plat,x_bnds=plon_bounds,y_bnds=plat_bounds, &
+ &  domain_id=domain_id,mode=mode)
+!-------------------------
+END SUBROUTINE histb_irreg
+!===
+SUBROUTINE histb_all &
+ & (k_typ,nc_name,pim,pjm,pitau0,pdate0,pdeltat,phoriid,idf, &
+ &  x_1d,y_1d,x_2d,y_2d,k_orx,k_szx,k_ory,k_szy, &
+ &  x_bnds,y_bnds,domain_id,mode)
+!---------------------------------------------------------------------
+!- General interface for horizontal grids.
+!- This subroutine initializes a netcdf file and returns the ID.
+!- It will set up the geographical space on which the data will be
+!- stored and offers the possibility of seting a zoom.
+!- In the case of irregular grids, all the data comes in as vectors
+!- and for the grid we have the coordinates of the 4 corners.
+!- It also gets the global parameters into the I/O subsystem.
+!-
+!- INPUT
+!-
+!- k_typ    : Type of the grid (1 rectilinear, 2 regular, 3 irregular)
+!- nc_name  : Name of the netcdf file to be created
+!- pim      : Size of arrays in longitude direction
+!- pjm      : Size of arrays in latitude direction (pjm=pim for type 3)
+!-
+!- pitau0   : time step at which the history tape starts
+!- pdate0   : The Julian date at which the itau was equal to 0
+!- pdeltat  : Time step, in seconds, of the counter itau
+!-            used in histwrite for instance
+!-
+!- OUTPUT
+!-
+!- phoriid  : Identifier of the horizontal grid
+!- idf      : Identifier of the file
+!-
+!- Optional INPUT arguments
+!-
+!- For rectilinear or irregular grid
+!- x_1d  : The longitudes
+!- y_1d  : The latitudes
+!- For regular grid
+!- x_2d  : The longitudes
+!- y_2d  : The latitudes
+!- One pair (x_1d,y_1d) or (x_2d,y_2d) must be supplied.
+!-
+!- For regular grid (reg1d or reg2d),                   
+!- the next 4 arguments allow to define a horizontal zoom
+!- for this file. It is assumed that all variables to come
+!- have the same index space. This can not be assumed for
+!- the z axis and thus we define the zoom in histdef.
+!- k_orx  : Origin of the slab of data within the X axis (pim)
+!- k_szx  : Size of the slab of data in X
+!- k_ory  : Origin of the slab of data within the Y axis (pjm)
+!- k_szy  : Size of the slab of data in Y
+!-
+!- For irregular grid.
+!- x_bnds : The boundaries of the grid in longitude
+!- y_bnds : The boundaries of the grid in latitude
+!-
+!- For all grids.
+!-
+!- domain_id  : Domain identifier
+!-
+!- mode       : String of (case insensitive) blank-separated words
+!-              defining the mode used to create the file.
+!-              Supported keywords : 32, 64
+!-              "32/64" defines the offset mode.
+!-              The default offset mode is 64 bits.
+!-              Keywords "NETCDF4" and "CLASSIC" are reserved
+!-              for future use.
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+  INTEGER,INTENT(IN) :: k_typ
+  CHARACTER(LEN=*),INTENT(IN) :: nc_name
+  INTEGER,INTENT(IN) :: pim,pjm
+  INTEGER,INTENT(IN) :: pitau0
+  REAL,INTENT(IN) :: pdate0,pdeltat
+  INTEGER,INTENT(OUT) :: idf,phoriid
+  REAL,DIMENSION(:),INTENT(IN),OPTIONAL :: x_1d,y_1d
+  REAL,DIMENSION(:,:),INTENT(IN),OPTIONAL :: x_2d,y_2d
+  INTEGER,INTENT(IN),OPTIONAL :: k_orx,k_szx,k_ory,k_szy
+  REAL,DIMENSION(:,:),INTENT(IN),OPTIONAL :: x_bnds,y_bnds
+  INTEGER,INTENT(IN),OPTIONAL :: domain_id
+  CHARACTER(LEN=*),OPTIONAL,INTENT(IN) :: mode
 !-
   INTEGER :: nfid,iret,m_c
   CHARACTER(LEN=120) :: file
   CHARACTER(LEN=30) :: timenow
+  CHARACTER(LEN=11) :: c_nam
   LOGICAL :: l_dbg
 !---------------------------------------------------------------------
   CALL ipsldbg (old_status=l_dbg)
 !-
-  IF (l_dbg) WRITE(*,*) "histbeg_irregular 0.0"
-!-
-  nb_files = nb_files+1
-  IF (nb_files > nb_files_max) THEN
+  IF     (k_typ == 1) THEN
+    c_nam = 'histb_reg1d'
+  ELSEIF (k_typ == 2) THEN
+    c_nam = 'histb_reg2d'
+  ELSEIF (k_typ == 3) THEN
+    c_nam = 'histb_irreg'
+  ELSE
     CALL ipslerr (3,"histbeg", &
- &    'Table of files too small. You should increase nb_files_max', &
- &    'in histcom.f90 in order to accomodate all these files',' ')
-  ENDIF
-  idf = nb_files
+ &    'Illegal value of k_typ argument','in internal interface','?')
+  ENDIF
+!-
+  IF (l_dbg) WRITE(*,*) c_nam//" 0.0"
+!-
+! Search for a free index
+!-
+  idf = -1
+  DO nfid=1,nb_files_max
+    IF (W_F(nfid)%ncfid < 0) THEN
+      idf = nfid; EXIT;
+    ENDIF
+  ENDDO
+  IF (idf < 0) THEN
+    CALL ipslerr (3,"histbeg", &
+   &  'Table of files too small. You should increase nb_files_max', &
+   &  'in histcom.f90 in order to accomodate all these files',' ')
+  ENDIF
 !-
 ! 1.0 Transfering into the common for future use
 !-
-  IF (l_dbg) WRITE(*,*) "histbeg_irregular 1.0"
+  IF (l_dbg) WRITE(*,*) c_nam//" 1.0"
 !-
   W_F(idf)%itau0  = pitau0
@@ -509,5 +364,5 @@
 ! 2.0 Initializes all variables for this file
 !-
-  IF (l_dbg) WRITE(*,*) "histbeg_irregular 2.0"
+  IF (l_dbg) WRITE(*,*) c_nam//" 2.0"
 !-
   W_F(idf)%n_var = 0
@@ -516,26 +371,33 @@
   W_F(idf)%n_zax = 0
 !-
-  W_F(idf)%slab_ori(1:2) = (/ 1,1 /)
-  W_F(idf)%slab_sz(1:2)  = (/ pim,1 /)
+  IF ( (k_typ == 1).OR.(k_typ == 2) ) THEN
+    W_F(idf)%slab_ori(1:2) = (/ k_orx,k_ory /)
+    W_F(idf)%slab_siz(1:2)  = (/ k_szx,k_szy /)
+  ELSE
+    W_F(idf)%slab_ori(1:2) = (/ 1,1 /)
+    W_F(idf)%slab_siz(1:2) = (/ pim,1 /)
+  ENDIF
 !-
 ! 3.0 Opening netcdf file and defining dimensions
 !-
-  IF (l_dbg) WRITE(*,*) "histbeg_irregular 3.0"
+  IF (l_dbg) WRITE(*,*) c_nam//" 3.0"
 !-
 ! Add DOMAIN number and ".nc" suffix in file name if needed
 !-
-  file  = pfilename
+  file = nc_name
   CALL flio_dom_file (file,domain_id)
 !-
 ! Check the mode
-  IF (PRESENT(nb_bits)) THEN
-    SELECT CASE (nb_bits)
-    CASE(32)
+!? See fliocom for HDF4 ????????????????????????????????????????????????
+!-
+  IF (PRESENT(mode)) THEN
+    SELECT CASE (TRIM(mode))
+    CASE('32')
       m_c = NF90_CLOBBER
-    CASE(64)
+    CASE('64')
       m_c = IOR(NF90_CLOBBER,NF90_64BIT_OFFSET)
     CASE DEFAULT
       CALL ipslerr (3,"histbeg", &
- &      'Invalid argument nb_bits for file :',TRIM(file), &
+ &      'Invalid argument mode for file :',TRIM(file), &
  &      'Supported values are 32 or 64')
     END SELECT
@@ -545,25 +407,32 @@
 !-
 ! Create file
-  iret = NF90_CREATE (file,m_c,nfid)
-!-
-  iret = NF90_DEF_DIM (nfid,'x',pim,W_F(idf)%xid)
-  W_F(idf)%yid = 0
+!-
+  iret = NF90_CREATE(file,m_c,nfid)
+!-
+  IF     (k_typ == 1) THEN
+    iret = NF90_DEF_DIM(nfid,'lon',k_szx,W_F(idf)%xid)
+    iret = NF90_DEF_DIM(nfid,'lat',k_szy,W_F(idf)%yid)
+  ELSEIF (k_typ == 2) THEN
+    iret = NF90_DEF_DIM(nfid,'x',k_szx,W_F(idf)%xid)
+    iret = NF90_DEF_DIM(nfid,'y',k_szy,W_F(idf)%yid)
+  ELSEIF (k_typ == 3) THEN
+    iret = NF90_DEF_DIM(nfid,'x',pim,W_F(idf)%xid)
+    W_F(idf)%yid = W_F(idf)%xid
+  ENDIF
 !-
 ! 4.0 Declaring the geographical coordinates and other attributes
 !-
-  IF (l_dbg) WRITE(*,*) "histbeg_irregular 4.0"
-!-
-! 4.3 Global attributes
-!-
-  iret = NF90_PUT_ATT (nfid,NF90_GLOBAL,'Conventions','CF-1.1')
-  iret = NF90_PUT_ATT (nfid,NF90_GLOBAL,'file_name',TRIM(file))
-  iret = NF90_PUT_ATT (nfid,NF90_GLOBAL,'production',TRIM(model_name))
+  IF (l_dbg) WRITE(*,*) c_nam//" 4.0"
+!-
+  iret = NF90_PUT_ATT(nfid,NF90_GLOBAL,'Conventions','CF-1.1')
+  iret = NF90_PUT_ATT(nfid,NF90_GLOBAL,'file_name',TRIM(file))
+  iret = NF90_PUT_ATT(nfid,NF90_GLOBAL,'production',TRIM(model_name))
   lock_modname = .TRUE.
   CALL ioget_timestamp (timenow)
-  iret = NF90_PUT_ATT (nfid,NF90_GLOBAL,'TimeStamp',TRIM(timenow))
+  iret = NF90_PUT_ATT(nfid,NF90_GLOBAL,'TimeStamp',TRIM(timenow))
 !-
 ! 5.0 Saving some important information on this file in the common
 !-
-  IF (l_dbg) WRITE(*,*) "histbeg_irregular 5.0"
+  IF (l_dbg) WRITE(*,*) c_nam//" 5.0"
 !-
   IF (PRESENT(domain_id)) THEN
@@ -571,20 +440,89 @@
   ENDIF
   W_F(idf)%ncfid = nfid
-  W_F(idf)%full_size(1:2) = (/ pim,1 /)
+  IF ( (k_typ == 1).OR.(k_typ == 2) ) THEN
+    W_F(idf)%full_size(1:2) = (/ pim,pjm /)
+    W_F(idf)%regular=.TRUE.
+  ELSEIF (k_typ == 3) THEN
+    W_F(idf)%full_size(1:2) = (/ pim,1 /)
+    W_F(idf)%regular=.FALSE.
+  ENDIF
 !-
 ! 6.0 storing the geographical coordinates
 !-
-  W_F(idf)%regular=.FALSE.
-!-
-  CALL histhori_irregular &
- &  (idf,pim,plon,plon_bounds,plat,plat_bounds, &
- &   ' ','Default grid',phoriid)
-!-------------------------------
-END SUBROUTINE histbeg_irregular
+  IF     (k_typ == 1) THEN
+    CALL histh_all &
+ &   (k_typ,idf,pim,pjm,' ','Default grid',phoriid, &
+ &    x_1d=x_1d,y_1d=y_1d)
+  ELSEIF (k_typ == 2) THEN
+    CALL histh_all &
+ &   (k_typ,idf,pim,pjm,' ','Default grid',phoriid, &
+ &    x_2d=x_2d,y_2d=y_2d)
+  ELSEIF (k_typ == 3) THEN
+    CALL histh_all &
+ &   (k_typ,idf,pim,pim,' ','Default grid',phoriid, &
+ &    x_1d=x_1d,y_1d=y_1d,x_bnds=x_bnds,y_bnds=y_bnds)
+  ENDIF
+!-----------------------
+END SUBROUTINE histb_all
 !===
-SUBROUTINE histhori_regular &
- &  (idf,pim,plon,pjm,plat,phname,phtitle,phid,opt_rectilinear)
-!---------------------------------------------------------------------
-!- This subroutine is made to declare a new horizontale grid.
+SUBROUTINE histh_reg1d &
+ &  (idf,pim,plon,pjm,plat,phname,phtitle,phid)
+!---------------------------------------------------------------------
+!- histhori for 1d regular grid (see histh_all)
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+  INTEGER,INTENT(IN) :: idf,pim,pjm
+  REAL,INTENT(IN),DIMENSION(:) :: plon,plat
+  CHARACTER(LEN=*),INTENT(IN) :: phname,phtitle
+  INTEGER,INTENT(OUT) :: phid
+!---------------------------------------------------------------------
+  CALL histh_all &
+ & (1,idf,pim,pjm,phname,phtitle,phid,x_1d=plon,y_1d=plat)
+!-------------------------
+END SUBROUTINE histh_reg1d
+!===
+SUBROUTINE histh_reg2d &
+ & (idf,pim,plon,pjm,plat,phname,phtitle,phid)
+!---------------------------------------------------------------------
+!- histhori for 2d regular grid (see histh_all)
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+  INTEGER,INTENT(IN) :: idf,pim,pjm
+  REAL,INTENT(IN),DIMENSION(:,:) :: plon,plat
+  CHARACTER(LEN=*),INTENT(IN) :: phname,phtitle
+  INTEGER,INTENT(OUT) :: phid
+!---------------------------------------------------------------------
+  CALL histh_all &
+ & (2,idf,pim,pjm,phname,phtitle,phid,x_2d=plon,y_2d=plat)
+!-------------------------
+END SUBROUTINE histh_reg2d
+!===
+SUBROUTINE histh_irreg &
+ & (idf,pim,plon,plon_bounds,plat,plat_bounds,phname,phtitle,phid)
+!---------------------------------------------------------------------
+!- histhori for irregular grid (see histh_all)
+!---------------------------------------------------------------------
+  IMPLICIT NONE
+!-
+  INTEGER,INTENT(IN) :: idf,pim
+  REAL,DIMENSION(:),INTENT(IN) :: plon,plat
+  REAL,DIMENSION(:,:),INTENT(IN) :: plon_bounds,plat_bounds
+  CHARACTER(LEN=*),INTENT(IN) :: phname,phtitle
+  INTEGER,INTENT(OUT) :: phid
+!---------------------------------------------------------------------
+  CALL histh_all &
+ & (3,idf,pim,pim,phname,phtitle,phid, &
+ &  x_1d=plon,y_1d=plat,x_bnds=plon_bounds,y_bnds=plat_bounds)
+!-------------------------
+END SUBROUTINE histh_irreg
+!===
+SUBROUTINE histh_all &
+ & (k_typ,idf,pim,pjm,phname,phtitle,phid, &
+ &  x_1d,y_1d,x_2d,y_2d,x_bnds,y_bnds)
+!---------------------------------------------------------------------
+!- General interface for horizontal grids.
+!- This subroutine is made to declare a new horizontal grid.
 !- It has to have the same number of points as
 !- the original and thus in this routine we will only
@@ -596,9 +534,8 @@
 !- INPUT
 !-
+!- k_typ   : Type of the grid (1 rectilinear, 2 regular, 3 irregular)
 !- idf     : The id of the file to which the grid should be added
 !- pim     : Size in the longitude direction
-!- plon    : The longitudes
-!- pjm     : Size in the latitude direction
-!- plat    : The latitudes
+!- pjm     : Size in the latitude direction (pjm=pim for type 3)
 !- phname  : The name of grid
 !- phtitle : The title of the grid
@@ -608,33 +545,61 @@
 !- phid    : Id of the created grid
 !-
-!- OPTIONAL
-!-
-!- opt_rectilinear : If true we know the grid is rectilinear.
-!-
+!- Optional INPUT arguments
+!-
+!- For rectilinear or irregular grid
+!- x_1d  : The longitudes
+!- y_1d  : The latitudes
+!- For regular grid
+!- x_2d  : The longitudes
+!- y_2d  : The latitudes
+!- One pair (x_1d,y_1d) or (x_2d,y_2d) must be supplied.
+!-
+!- For irregular grid.
+!- x_bnds : The boundaries of the grid in longitude
+!- y_bnds : The boundaries of the grid in latitude
 !---------------------------------------------------------------------
   IMPLICIT NONE
 !-
+  INTEGER,INTENT(IN) :: k_typ
   INTEGER,INTENT(IN) :: idf,pim,pjm
-  REAL,INTENT(IN),DIMENSION(pim,pjm) :: plon,plat
   CHARACTER(LEN=*),INTENT(IN) :: phname,phtitle
   INTEGER,INTENT(OUT) :: phid
-  LOGICAL,INTENT(IN),OPTIONAL :: opt_rectilinear
+  REAL,DIMENSION(:),INTENT(IN),OPTIONAL :: x_1d,y_1d
+  REAL,DIMENSION(:,:),INTENT(IN),OPTIONAL :: x_2d,y_2d
+  REAL,DIMENSION(:,:),INTENT(IN),OPTIONAL :: x_bnds,y_bnds
 !-
   CHARACTER(LEN=25) :: lon_name,lat_name
-  CHARACTER(LEN=80) :: tmp_title,tmp_name
-  INTEGER :: ndim
-  INTEGER,DIMENSION(2) :: dims
+  CHARACTER(LEN=30) :: lonbound_name,latbound_name
+  INTEGER :: i_s,i_e
+  INTEGER,DIMENSION(2) :: dims,dims_b
+  INTEGER :: nbbounds
+  INTEGER :: nlonidb,nlatidb,twoid
+  LOGICAL :: transp = .FALSE.
+  REAL,ALLOCATABLE,DIMENSION(:,:) :: bounds_trans
+  REAL :: wmn,wmx
   INTEGER :: nlonid,nlatid
-  INTEGER :: orix,oriy,par_szx,par_szy
+  INTEGER :: o_x,o_y,s_x,s_y
   INTEGER :: iret,nfid
-  LOGICAL :: rectilinear
+  CHARACTER(LEN=11) :: c_nam
   LOGICAL :: l_dbg
 !---------------------------------------------------------------------
   CALL ipsldbg (old_status=l_dbg)
 !-
+  IF     (k_typ == 1) THEN
+    c_nam = 'histh_reg1d'
+  ELSEIF (k_typ == 2) THEN
+    c_nam = 'histh_reg2d'
+  ELSEIF (k_typ == 3) THEN
+    c_nam = 'histh_irreg'
+  ELSE
+    CALL ipslerr (3,"histhori", &
+ &    'Illegal value of k_typ argument','in internal interface','?')
+  ENDIF
+!-
 ! 1.0 Check that all fits in the buffers
 !-
   IF (    (pim /= W_F(idf)%full_size(1)) &
- &    .OR.(pjm /= W_F(idf)%full_size(2)) ) THEN
+ &    .OR.(W_F(idf)%regular.AND.(pjm /= W_F(idf)%full_size(2)))  &
+ &    .OR.(.NOT.W_F(idf)%regular.AND.(W_F(idf)%full_size(2) /= 1)) ) THEN
     CALL ipslerr (3,"histhori", &
  &   'The new horizontal grid does not have the same size', &
@@ -643,36 +608,55 @@
   ENDIF
 !-
-  IF (PRESENT(opt_rectilinear)) THEN
-    rectilinear = opt_rectilinear
-  ELSE
-    rectilinear = .FALSE.
-  ENDIF
-!-
 ! 1.1 Create all the variables needed
 !-
-  IF (l_dbg) WRITE(*,*) "histhori_regular 1.0"
+  IF (l_dbg) WRITE(*,*) c_nam//" 1.0"
 !-
   nfid = W_F(idf)%ncfid
 !-
-  ndim = 2
+  IF (k_typ == 3) THEN
+    IF     (SIZE(x_bnds,DIM=1) == pim) THEN
+      nbbounds = SIZE(x_bnds,DIM=2)
+      transp = .TRUE.
+    ELSEIF (SIZE(x_bnds,DIM=2) == pim) THEN
+      nbbounds = SIZE(x_bnds,DIM=1)
+      transp = .FALSE.
+    ELSE
+      CALL ipslerr (3,"histhori", &
+ &     'The boundary variable does not have any axis corresponding', &
+ &     'to the size of the longitude or latitude variable','.')
+    ENDIF
+    ALLOCATE(bounds_trans(nbbounds,pim))
+    iret = NF90_DEF_DIM(nfid,'nbnd',nbbounds,twoid)
+    dims_b(1:2) = (/ twoid,W_F(idf)%xid /)
+  ENDIF
+!-
   dims(1:2) = (/ W_F(idf)%xid,W_F(idf)%yid /)
 !-
-  tmp_name = phname
-  IF (rectilinear) THEN
+  IF     (k_typ == 1) THEN
     IF (W_F(idf)%n_hax == 0) THEN
       lon_name = 'lon'
       lat_name = 'lat'
     ELSE
-      lon_name = 'lon_'//TRIM(tmp_name)
-      lat_name = 'lat_'//TRIM(tmp_name)
-    ENDIF
-  ELSE
+      lon_name = 'lon_'//TRIM(phname)
+      lat_name = 'lat_'//TRIM(phname)
+    ENDIF
+  ELSEIF (k_typ == 2) THEN
     IF (W_F(idf)%n_hax == 0) THEN
       lon_name = 'nav_lon'
       lat_name = 'nav_lat'
     ELSE
-      lon_name = 'nav_lon_'//TRIM(tmp_name)
-      lat_name = 'nav_lat_'//TRIM(tmp_name)
-    ENDIF
+      lon_name = 'nav_lon_'//TRIM(phname)
+      lat_name = 'nav_lat_'//TRIM(phname)
+    ENDIF
+  ELSEIF (k_typ == 3) THEN
+    IF (W_F(idf)%n_hax == 0) THEN
+      lon_name = 'nav_lon'
+      lat_name = 'nav_lat'
+    ELSE
+      lon_name = 'nav_lon_'//TRIM(phname)
+      lat_name = 'nav_lat_'//TRIM(phname)
+    ENDIF
+    lonbound_name = TRIM(lon_name)//'_bounds'
+    latbound_name = TRIM(lat_name)//'_bounds'
   ENDIF
 !-
@@ -682,266 +666,121 @@
   W_F(idf)%n_hax = phid
   W_F(idf)%hax_name(phid,1:2) = (/ lon_name,lat_name /)
-  tmp_title = phtitle
 !-
 ! 2.0 Longitude
 !-
-  IF (l_dbg) WRITE(*,*) "histhori_regular 2.0"
-!-
-  IF (rectilinear) THEN
-    ndim = 1
-    dims(1:1) = (/ W_F(idf)%xid /)
-  ENDIF
-  iret = NF90_DEF_VAR (nfid,lon_name,NF90_REAL4,dims(1:ndim),nlonid)
-  IF (rectilinear) THEN
-    iret = NF90_PUT_ATT (nfid,nlonid,'axis',"X")
-  ENDIF
-  iret = NF90_PUT_ATT (nfid,nlonid,'standard_name',"longitude")
-  iret = NF90_PUT_ATT (nfid,nlonid,'units',"degrees_east")
-  iret = NF90_PUT_ATT (nfid,nlonid,'valid_min', &
- &                     REAL(MINVAL(plon),KIND=4))
-  iret = NF90_PUT_ATT (nfid,nlonid,'valid_max', &
- &                     REAL(MAXVAL(plon),KIND=4))
-  iret = NF90_PUT_ATT (nfid,nlonid,'long_name',"Longitude")
-  iret = NF90_PUT_ATT (nfid,nlonid,'nav_model',TRIM(tmp_title))
+  IF (l_dbg) WRITE(*,*) c_nam//" 2.0"
+!-
+  i_s = 1;
+  IF ( (k_typ == 1).OR.(k_typ == 3) ) THEN
+    i_e = 1; wmn = MINVAL(x_1d); wmx = MAXVAL(x_1d);
+  ELSEIF (k_typ == 2) THEN
+    i_e = 2; wmn = MINVAL(x_2d); wmx = MAXVAL(x_2d);
+  ENDIF
+  iret = NF90_DEF_VAR(nfid,lon_name,NF90_REAL4,dims(i_s:i_e),nlonid)
+  IF (k_typ == 1) THEN
+    iret = NF90_PUT_ATT(nfid,nlonid,'axis',"X")
+  ENDIF
+  iret = NF90_PUT_ATT(nfid,nlonid,'standard_name',"longitude")
+  iret = NF90_PUT_ATT(nfid,nlonid,'units',"degrees_east")
+  iret = NF90_PUT_ATT(nfid,nlonid,'valid_min',REAL(wmn,KIND=4))
+  iret = NF90_PUT_ATT(nfid,nlonid,'valid_max',REAL(wmx,KIND=4))
+  iret = NF90_PUT_ATT(nfid,nlonid,'long_name',"Longitude")
+  iret = NF90_PUT_ATT(nfid,nlonid,'nav_model',TRIM(phtitle))
+!-
+  IF (k_typ == 3) THEN
+!---
+!-- 2.1 Longitude bounds
+!---
+    iret = NF90_PUT_ATT(nfid,nlonid,'bounds',TRIM(lonbound_name))
+    iret = NF90_DEF_VAR(nfid,lonbound_name,NF90_REAL4,dims_b(1:2),nlonidb)
+    iret = NF90_PUT_ATT(nfid,nlonidb,'long_name', &
+ &          'Boundaries for coordinate variable '//TRIM(lon_name))
+  ENDIF
 !-
 ! 3.0 Latitude
 !-
-  IF (l_dbg) WRITE(*,*) "histhori_regular 3.0"
-!-
-  IF (rectilinear) THEN
-    ndim = 1
-    dims(1:1) = (/ W_F(idf)%yid /)
-  ENDIF
-  iret = NF90_DEF_VAR (nfid,lat_name,NF90_REAL4,dims(1:ndim),nlatid)
-  IF (rectilinear) THEN
-    iret = NF90_PUT_ATT (nfid,nlatid,'axis',"Y")
-  ENDIF
-  iret = NF90_PUT_ATT (nfid,nlatid,'standard_name',"latitude")
-  iret = NF90_PUT_ATT (nfid,nlatid,'units',"degrees_north")
-  iret = NF90_PUT_ATT (nfid,nlatid,'valid_min', &
- &                     REAL(MINVAL(plat),KIND=4))
-  iret = NF90_PUT_ATT (nfid,nlatid,'valid_max', &
- &                     REAL(MAXVAL(plat),KIND=4))
-  iret = NF90_PUT_ATT (nfid,nlatid,'long_name',"Latitude")
-  iret = NF90_PUT_ATT (nfid,nlatid,'nav_model',TRIM(tmp_title))
-!-
-  iret = NF90_ENDDEF (nfid)
+  IF (l_dbg) WRITE(*,*) c_nam//" 3.0"
+!-
+  i_e = 2;
+  IF ( (k_typ == 1).OR.(k_typ == 3) ) THEN
+    i_s = 2; wmn = MINVAL(y_1d); wmx = MAXVAL(y_1d);
+  ELSEIF (k_typ == 2) THEN
+    i_s = 1; wmn = MINVAL(y_2d); wmx = MAXVAL(y_2d);
+  ENDIF
+  iret = NF90_DEF_VAR(nfid,lat_name,NF90_REAL4,dims(i_s:i_e),nlatid)
+  IF (k_typ == 1) THEN
+    iret = NF90_PUT_ATT(nfid,nlatid,'axis',"Y")
+  ENDIF
+!-
+  iret = NF90_PUT_ATT(nfid,nlatid,'standard_name',"latitude")
+  iret = NF90_PUT_ATT(nfid,nlatid,'units',"degrees_north")
+  iret = NF90_PUT_ATT(nfid,nlatid,'valid_min',REAL(wmn,KIND=4))
+  iret = NF90_PUT_ATT(nfid,nlatid,'valid_max',REAL(wmx,KIND=4))
+  iret = NF90_PUT_ATT(nfid,nlatid,'long_name',"Latitude")
+  iret = NF90_PUT_ATT(nfid,nlatid,'nav_model',TRIM(phtitle))
+!-
+  IF (k_typ == 3) THEN
+!---
+!-- 3.1 Latitude bounds
+!---
+    iret = NF90_PUT_ATT(nfid,nlatid,'bounds',TRIM(latbound_name))
+    iret = NF90_DEF_VAR(nfid,latbound_name,NF90_REAL4,dims_b(1:2),nlatidb)
+    iret = NF90_PUT_ATT(nfid,nlatidb,'long_name', &
+ &          'Boundaries for coordinate variable '//TRIM(lat_name))
+  ENDIF
+!-
+  iret = NF90_ENDDEF(nfid)
 !-
 ! 4.0 storing the geographical coordinates
 !-
-  IF (l_dbg) WRITE(*,*) "histhori_regular 4.0"
-!-
-  orix = W_F(idf)%slab_ori(1)
-  oriy = W_F(idf)%slab_ori(2)
-  par_szx = W_F(idf)%slab_sz(1)
-  par_szy = W_F(idf)%slab_sz(2)
-!-
-! Transfer the longitude
-!-
-  IF (rectilinear) THEN
-    iret = NF90_PUT_VAR (nfid,nlonid,plon(orix:orix+par_szx-1,1))
-  ELSE
-    iret = NF90_PUT_VAR (nfid,nlonid, &
- &           plon(orix:orix+par_szx-1,oriy:oriy+par_szy-1))
-  ENDIF
-!-
-! Transfer the latitude
-!-
-  IF (rectilinear) THEN
-    iret = NF90_PUT_VAR (nfid,nlatid,plat(1,oriy:oriy+par_szy-1))
-  ELSE
-    iret = NF90_PUT_VAR (nfid,nlatid, &
- &           plat(orix:orix+par_szx-1,oriy:oriy+par_szy-1))
-  ENDIF
-!-
-  iret = NF90_REDEF (nfid)
-!------------------------------
-END SUBROUTINE histhori_regular
-!===
-SUBROUTINE histhori_irregular &
- &  (idf,pim,plon,plon_bounds,plat,plat_bounds, &
- &   phname,phtitle,phid)
-!---------------------------------------------------------------------
-!- This subroutine is made to declare a new horizontale grid.
-!- It has to have the same number of points as
-!- the original and thus in this routine we will only
-!- add two variable (longitude and latitude).
-!- Any variable in the file can thus point to this pair
-!- through an attribute. This routine is very usefull
-!- to allow staggered grids.
-!-
-!- INPUT
-!-
-!- idf         : The id of the file to which the grid should be added
-!- pim         : Size in the longitude direction
-!- plon        : The longitudes
-!- plon_bounds : The boundaries of the grid in longitude
-!- plat        : The latitudes
-!- plat_bounds : Boundaries of the grid in latitude
-!- phname      : The name of grid
-!- phtitle     : The title of the grid
-!-
-!- OUTPUT
-!-
-!- phid    : Id of the created grid
-!---------------------------------------------------------------------
-  IMPLICIT NONE
-!-
-  INTEGER,INTENT(IN) :: idf,pim
-  REAL,DIMENSION(pim),INTENT(IN) :: plon,plat
-  REAL,DIMENSION(:,:),INTENT(IN) :: plon_bounds,plat_bounds
-  CHARACTER(LEN=*),INTENT(IN) :: phname,phtitle
-  INTEGER,INTENT(OUT) :: phid
-!-
-  CHARACTER(LEN=25) :: lon_name,lat_name
-  CHARACTER(LEN=30) :: lonbound_name,latbound_name
-  CHARACTER(LEN=80) :: tmp_title,tmp_name,longname
-  INTEGER :: ndim,dims(2)
-  INTEGER :: ndimb,dimsb(2)
-  INTEGER :: nbbounds
-  INTEGER :: nlonid,nlatid,nlonidb,nlatidb
-  INTEGER :: iret,nfid,twoid
-  LOGICAL :: transp = .FALSE.
-  REAL,ALLOCATABLE,DIMENSION(:,:) :: bounds_trans
-  LOGICAL :: l_dbg
-!---------------------------------------------------------------------
-  CALL ipsldbg (old_status=l_dbg)
-!-
-! 1.0 Check that all fits in the buffers
-!-
-  IF (    (pim /= W_F(idf)%full_size(1)) &
- &    .OR.(W_F(idf)%full_size(2) /= 1) ) THEN
-    CALL ipslerr (3,"histhori", &
- &   'The new horizontal grid does not have the same size', &
- &   'as the one provided to histbeg. This is not yet ', &
- &   'possible in the hist package.')
-  ENDIF
-!-
-! 1.1 Create all the variables needed
-!-
-  IF (l_dbg) WRITE(*,*) 'histhori_irregular 1.0'
-!-
-  nfid = W_F(idf)%ncfid
-!-
-  IF     (SIZE(plon_bounds,DIM=1) == pim) THEN
-    nbbounds = SIZE(plon_bounds,DIM=2)
-    transp = .TRUE.
-  ELSEIF (SIZE(plon_bounds,DIM=2) == pim) THEN
-    nbbounds = SIZE(plon_bounds,DIM=1)
-    transp = .FALSE.
-  ELSE
-    CALL ipslerr (3,"histhori", &
-   &  'The boundary variable does not have any axis corresponding', &
-   &  'to the size of the longitude or latitude variable', &
-   &  '.')
-  ENDIF
-!-
-  ALLOCATE(bounds_trans(nbbounds,pim))
-!-
-  iret = NF90_DEF_DIM (nfid,'nbnd',nbbounds,twoid)
-  ndim = 1
-  dims(1) = W_F(idf)%xid
-  ndimb = 2
-  dimsb(1:2) = (/ twoid,W_F(idf)%xid /)
-!-
-  tmp_name = phname
-  IF (W_F(idf)%n_hax == 0) THEN
-    lon_name = 'nav_lon'
-    lat_name = 'nav_lat'
-  ELSE
-    lon_name = 'nav_lon_'//TRIM(tmp_name)
-    lat_name = 'nav_lat_'//TRIM(tmp_name)
-  ENDIF
-  lonbound_name = TRIM(lon_name)//'_bounds'
-  latbound_name = TRIM(lat_name)//'_bounds'
-!-
-! 1.2 Save the informations
-!-
-  phid = W_F(idf)%n_hax+1
-  W_F(idf)%n_hax = phid
-  W_F(idf)%hax_name(phid,1:2) = (/ lon_name,lat_name /)
-  tmp_title = phtitle
-!-
-! 2.0 Longitude
-!-
-  IF (l_dbg) WRITE(*,*) "histhori_irregular 2.0"
-!-
-  iret = NF90_DEF_VAR (nfid,lon_name,NF90_REAL4,dims(1:ndim),nlonid)
-  iret = NF90_PUT_ATT (nfid,nlonid,'standard_name',"longitude")
-  iret = NF90_PUT_ATT (nfid,nlonid,'units',"degrees_east")
-  iret = NF90_PUT_ATT (nfid,nlonid,'valid_min', &
- &                     REAL(MINVAL(plon),KIND=4))
-  iret = NF90_PUT_ATT (nfid,nlonid,'valid_max', &
- &                     REAL(MAXVAL(plon),KIND=4))
-  iret = NF90_PUT_ATT (nfid,nlonid,'long_name',"Longitude")
-  iret = NF90_PUT_ATT (nfid,nlonid,'nav_model',TRIM(tmp_title))
-!-
-! 2.1 Longitude bounds
-!-
-  iret = NF90_PUT_ATT (nfid,nlonid,'bounds',TRIM(lonbound_name))
-  iret = NF90_DEF_VAR (nfid,lonbound_name,NF90_REAL4, &
- &                     dimsb(1:ndimb),nlonidb)
-  longname = 'Boundaries for coordinate variable '//TRIM(lon_name)
-  iret = NF90_PUT_ATT (nfid,nlonidb,'long_name',TRIM(longname))
-!-
-! 3.0 Latitude
-!-
-  IF (l_dbg) WRITE(*,*) "histhori_irregular 3.0"
-!-
-  iret = NF90_DEF_VAR (nfid,lat_name,NF90_REAL4,dims(1:ndim),nlatid)
-  iret = NF90_PUT_ATT (nfid,nlatid,'standard_name',"latitude")
-  iret = NF90_PUT_ATT (nfid,nlatid,'units',"degrees_north")
-  iret = NF90_PUT_ATT (nfid,nlatid,'valid_min', &
- &                     REAL(MINVAL(plat),KIND=4))
-  iret = NF90_PUT_ATT (nfid,nlatid,'valid_max', &
- &                     REAL(MAXVAL(plat),KIND=4))
-  iret = NF90_PUT_ATT (nfid,nlatid,'long_name',"Latitude")
-  iret = NF90_PUT_ATT (nfid,nlatid,'nav_model',TRIM(tmp_title))
-!-
-! 3.1 Latitude bounds
-!-
-  iret = NF90_PUT_ATT (nfid,nlatid,'bounds',TRIM(latbound_name))
-  iret = NF90_DEF_VAR (nfid,latbound_name,NF90_REAL4, &
- &                     dimsb(1:ndimb),nlatidb)
-  longname = 'Boundaries for coordinate variable '//TRIM(lat_name)
-  iret = NF90_PUT_ATT (nfid,nlatidb,'long_name',TRIM(longname))
-!-
-  iret = NF90_ENDDEF (nfid)
-!-
-! 4.0 storing the geographical coordinates
-!-
-  IF (l_dbg) WRITE(*,*) "histhori_irregular 4.0"
-!-
-! 4.1 Write the longitude
-!-
-  iret = NF90_PUT_VAR (nfid,nlonid,plon(1:pim))
-!-
-! 4.2 Write the longitude bounds
-!-
-  IF (transp) THEN
-    bounds_trans = TRANSPOSE(plon_bounds)
-  ELSE
-    bounds_trans = plon_bounds
-  ENDIF
-  iret = NF90_PUT_VAR (nfid,nlonidb,bounds_trans(1:nbbounds,1:pim))
-!-
-! 4.3 Write the latitude
-!-
-  iret = NF90_PUT_VAR (nfid,nlatid,plat(1:pim))
-!-
-! 4.4 Write the latitude bounds
-!-
-  IF (transp) THEN
-    bounds_trans = TRANSPOSE(plat_bounds)
-  ELSE
-    bounds_trans = plat_bounds
-  ENDIF
-  iret = NF90_PUT_VAR (nfid,nlatidb,bounds_trans(1:nbbounds,1:pim))
-!-
-  DEALLOCATE(bounds_trans)
-!-
-  iret = NF90_REDEF (nfid)
-!--------------------------------
-END SUBROUTINE histhori_irregular
+  IF (l_dbg) WRITE(*,*) c_nam//" 4.0"
+!-
+  IF ( (k_typ == 1).OR.(k_typ == 2) ) THEN
+    o_x = W_F(idf)%slab_ori(1)
+    o_y = W_F(idf)%slab_ori(2)
+    s_x = W_F(idf)%slab_siz(1)
+    s_y = W_F(idf)%slab_siz(2)
+!---
+!-- Transfer the longitude and  the latitude
+!---
+    IF     (k_typ == 1) THEN
+      iret = NF90_PUT_VAR(nfid,nlonid,x_1d(o_x:o_x+s_x-1))
+      iret = NF90_PUT_VAR(nfid,nlatid,y_1d(o_y:o_y+s_y-1))
+    ELSEIF (k_typ == 2) THEN
+      iret = NF90_PUT_VAR(nfid,nlonid, &
+ &            x_2d(o_x:o_x+s_x-1,o_y:o_y+s_y-1))
+      iret = NF90_PUT_VAR(nfid,nlatid, &
+ &            y_2d(o_x:o_x+s_x-1,o_y:o_y+s_y-1))
+    ENDIF
+  ELSEIF (k_typ == 3) THEN
+!---
+!-- Transfer the longitude and the longitude bounds
+!---
+    iret = NF90_PUT_VAR(nfid,nlonid,x_1d(1:pim))
+!---
+    IF (transp) THEN
+      bounds_trans = TRANSPOSE(x_bnds)
+    ELSE
+      bounds_trans = x_bnds
+    ENDIF
+    iret = NF90_PUT_VAR(nfid,nlonidb,bounds_trans(1:nbbounds,1:pim))
+!---
+!-- Transfer the latitude and the latitude bounds
+!---
+    iret = NF90_PUT_VAR(nfid,nlatid,y_1d(1:pim))
+!---
+    IF (transp) THEN
+      bounds_trans = TRANSPOSE(y_bnds)
+    ELSE
+      bounds_trans = y_bnds
+    ENDIF
+    iret = NF90_PUT_VAR(nfid,nlatidb,bounds_trans(1:nbbounds,1:pim))
+!---
+    DEALLOCATE(bounds_trans)
+  ENDIF
+!-
+  iret = NF90_REDEF(nfid)
+!-----------------------
+END SUBROUTINE histh_all
 !===
 SUBROUTINE histvert (idf,pzaxname,pzaxtitle,pzaxunit, &
@@ -1239,5 +1078,5 @@
  &  (/ W_F(idf)%slab_ori(1),W_F(idf)%slab_ori(2),par_oriz /)
   W_F(idf)%W_V(iv)%zsize(1:3) = &
- &  (/ W_F(idf)%slab_sz(1),W_F(idf)%slab_sz(2),par_szz /)
+ &  (/ W_F(idf)%slab_siz(1),W_F(idf)%slab_siz(2),par_szz /)
 !-
 ! Is the size of the full array the same as that of the coordinates ?
@@ -1257,10 +1096,10 @@
 ! Is the size of the zoom smaller than the coordinates ?
 !-
-  IF (    (W_F(idf)%full_size(1) < W_F(idf)%slab_sz(1)) &
- &    .OR.(W_F(idf)%full_size(2) < W_F(idf)%slab_sz(2)) ) THEN
+  IF (    (W_F(idf)%full_size(1) < W_F(idf)%slab_siz(1)) &
+ &    .OR.(W_F(idf)%full_size(2) < W_F(idf)%slab_siz(2)) ) THEN
     str70 = &
  &   "Size of variable should be greater or equal to those of the zoom"
     WRITE(str71,'("Size of XY zoom :",2I4)') &
- &   W_F(idf)%slab_sz(1),W_F(idf)%slab_sz(2)
+ &   W_F(idf)%slab_siz(1),W_F(idf)%slab_siz(2)
     WRITE(str72,'("Size declared for variable ",A," :",2I4)') &
  &   TRIM(tmp_name),pxsize,pysize
@@ -1820,5 +1659,5 @@
 !     Thanks Marine for showing us what errors users can make !
 !-
-  IF ( (idf < 1).OR.(idf > nb_files) ) THEN
+  IF ( (idf < 1).OR.(idf > nb_files_max) ) THEN
     CALL ipslerr (3,"histwrite", &
  &    'Illegal file ID in the histwrite of variable',pvarname,' ')
@@ -2415,5 +2254,4 @@
   INTEGER :: ifile,iret,i_s,i_e
 !-
-  LOGICAL :: file_exists
   LOGICAL :: l_dbg
 !---------------------------------------------------------------------
@@ -2425,5 +2263,5 @@
 !-
   IF (PRESENT(file)) THEN
-    IF ( (file >= 1).AND.(file <= nb_files) ) THEN
+    IF ( (file >= 1).AND.(file <= nb_files_max) ) THEN
       IF (W_F(ifile)%ncfid > 0) THEN
         i_s = file
@@ -2440,5 +2278,5 @@
   ELSE
     i_s = 1
-    i_e = nb_files
+    i_e = nb_files_max
   ENDIF
 !-
@@ -2480,5 +2318,5 @@
 !-
   IF (PRESENT(idf)) THEN
-    IF ( (idf >= 1).AND.(idf <= nb_files) ) THEN
+    IF ( (idf >= 1).AND.(idf <= nb_files_max) ) THEN
       IF (W_F(ifile)%ncfid > 0) THEN
         i_s = idf
@@ -2495,5 +2333,5 @@
   ELSE
     i_s = 1
-    i_e = nb_files
+    i_e = nb_files_max
   ENDIF
 !-
