Index: codes/icosagcm/devel/src/output/set_bounds.f90
===================================================================
--- codes/icosagcm/devel/src/output/set_bounds.f90	(revision 880)
+++ codes/icosagcm/devel/src/output/set_bounds.f90	(revision 880)
@@ -0,0 +1,206 @@
+MODULE set_bounds_mod
+  USE geometry, ONLY : swap_geometry, lon_e, lat_e
+  USE dimensions, ONLY : swap_dimensions, u_pos
+  USE math_const, ONLY : Pi
+  USE domain_mod, ONLY : t_domain, t_cellset, domloc_glo_ind
+  USE spherical_geom_mod, ONLY : xyz2lonlat
+
+  IMPLICIT NONE
+  PRIVATE
+  SAVE
+
+  PUBLIC :: set_bounds
+
+CONTAINS
+
+  SUBROUTINE set_bounds_primal(d, cells, all, own)
+    TYPE(t_domain) :: d
+    TYPE(t_cellset) :: cells
+    LOGICAL :: all, own(:,:)         ! if all is .TRUE., include halo cells
+    REAL :: lon,lat
+    INTEGER :: i,j,k, n, halo_size
+
+    halo_size = MERGE(1,0,all)
+
+    ! count primal cells
+    n=0
+    DO j=d%jj_begin-halo_size, d%jj_end+halo_size
+       DO i=d%ii_begin-halo_size, d%ii_end+halo_size
+          IF (own(i,j) .OR. all ) n=n+1
+       ENDDO
+    ENDDO
+    cells%ncell = n
+    
+    ! now set bounds
+    ALLOCATE(cells%ij(n), cells%lon(n), cells%lat(n), cells%ind_glo(n))
+    ALLOCATE(cells%bnds_lon(0:5,n), cells%bnds_lat(0:5,n))
+
+    n=0
+    DO j=d%jj_begin-halo_size, d%jj_end+halo_size
+       DO i=d%ii_begin-halo_size, d%ii_end+halo_size
+          IF (own(i,j) .OR. all) THEN
+             n=n+1
+             CALL xyz2lonlat(d%xyz(:,i,j), lon, lat)
+             cells%lon(n)=lon*180./Pi
+             cells%lat(n)=lat*180./Pi
+             DO k=0,5
+                CALL xyz2lonlat(d%vertex(:,k,i,j), lon, lat)
+                cells%bnds_lon(k,n)=lon*180./Pi
+                cells%bnds_lat(k,n)=lat*180./Pi
+             END DO
+             cells%ij(n)=d%iim*(j-1)+i
+             cells%ind_glo(n) = d%assign_cell_glo(i,j)-1
+          END IF
+       END DO
+    END DO
+
+!    PRINT *, 'set_bounds_primal', all, halo_size, cells%ncell
+  END SUBROUTINE set_bounds_primal
+
+  SUBROUTINE set_bounds_dual(d, cells)
+    USE metric, ONLY : vup, vdown
+    TYPE(t_domain) :: d
+    TYPE(t_cellset) :: cells
+    REAL :: lonc, latc, lon(0:2), lat(0:2)
+    INTEGER :: i,j,k,n
+
+    ! count dual cells
+    n=0
+    DO j=d%jj_begin+1,d%jj_end
+       DO i=d%ii_begin,d%ii_end-1
+          n=n+2
+       ENDDO
+    ENDDO
+    cells%ncell = n
+
+    ! now set bounds
+    ALLOCATE(cells%ind_glo(n)) ! not set but must be allocated
+    ALLOCATE(cells%ij(n), cells%lon(n), cells%lat(n))
+    ALLOCATE(cells%bnds_lon(0:2,n), cells%bnds_lat(0:2,n))
+
+    n=0
+    DO j=d%jj_begin+1,d%jj_end
+       DO i=d%ii_begin,d%ii_end-1
+          n=n+1
+          CALL xyz2lonlat(d%vertex(:,vdown,i,j), lonc, latc)
+          CALL xyz2lonlat(d%xyz(:,i,j),     lon(0), lat(0))
+          CALL xyz2lonlat(d%xyz(:,i,j-1),   lon(1), lat(1))
+          CALL xyz2lonlat(d%xyz(:,i+1,j-1), lon(2), lat(2))
+          cells%lon(n)=lonc*180./Pi
+          cells%lat(n)=latc*180/Pi
+          DO k=0,2
+             cells%bnds_lat(k,n)=lat(k)*180./Pi
+             cells%bnds_lon(k,n)=lon(k)*180./Pi
+          END DO
+          cells%ij(n) = d%z_down + d%iim*(j-1)+i
+       ENDDO
+    ENDDO
+    
+    DO j=d%jj_begin,d%jj_end-1
+       DO i=d%ii_begin+1,d%ii_end
+          n=n+1
+          CALL xyz2lonlat(d%vertex(:,vup,i,j), lonc, latc)
+          CALL xyz2lonlat(d%xyz(:,i,j),     lon(0), lat(0))
+          CALL xyz2lonlat(d%xyz(:,i,j+1),   lon(1), lat(1))
+          CALL xyz2lonlat(d%xyz(:,i-1,j+1), lon(2), lat(2))
+          cells%lon(n)=lonc*180./Pi
+          cells%lat(n)=latc*180/Pi
+          DO k=0,2
+             cells%bnds_lat(k,n)=lat(k)*180./Pi
+             cells%bnds_lon(k,n)=lon(k)*180./Pi
+          END DO
+          cells%ij(n) = d%z_up + d%iim*(j-1)+i
+       ENDDO
+    ENDDO
+
+  END SUBROUTINE set_bounds_dual
+
+  SUBROUTINE set_bounds_edge(ind, d, cells)
+    USE metric, ONLY : cell_glo
+    INTEGER, INTENT(IN) :: ind
+    TYPE(t_domain) :: d
+    TYPE(t_cellset) :: cells
+    REAL :: lon(2), lat(2)
+    INTEGER :: i,j,ij,k,kk,n
+    
+    ! count edges
+    n=0
+    DO j=d%jj_begin,d%jj_end
+       DO i=d%ii_begin,d%ii_end
+          DO k=0,5
+             IF (d%edge_assign_domain(k,i,j)==domloc_glo_ind(ind) &
+                  .AND. d%edge_assign_i(k,i,j)==i &
+                  .AND. d%edge_assign_j(k,i,j)==j &
+                  .AND. d%edge_assign_pos(k,i,j)==k) n=n+1
+          END DO
+       END DO
+    END DO
+    cells%ncell = n
+
+    ! now set bounds
+    ALLOCATE(cells%ij(n), cells%lon(n), cells%lat(n), cells%ind_glo(n))
+    ALLOCATE(cells%sgn(n)) ! flip sign when reading/writing
+    ALLOCATE(cells%bnds_lon(2,n), cells%bnds_lat(2,n))
+
+    CALL swap_dimensions(ind)
+    CALL swap_geometry(ind)
+    
+    n=0
+    DO j=d%jj_begin,d%jj_end
+       DO i=d%ii_begin,d%ii_end
+          DO k=0,5
+             IF (d%edge_assign_domain(k,i,j)==domloc_glo_ind(ind) &
+                  .AND. d%edge_assign_i(k,i,j)==i &
+                  .AND. d%edge_assign_j(k,i,j)==j &
+                  .AND. d%edge_assign_pos(k,i,j)==k) THEN
+                n=n+1
+                ij=(j-1)*d%iim+i+u_pos(k+1)
+                kk = MOD(k+d%delta(i,j)+6,6)
+                cells%ij(n)     = ij
+                cells%sgn(n)    = d%edge_assign_sign(k,i,j)
+                cells%ind_glo(n)= cell_glo(d%assign_cell_glo(i,j))%edge(kk)-1
+                cells%lon(n)    = lon_e(ij)*180./Pi
+                cells%lat(n)    = lat_e(ij)*180./Pi
+               
+                kk = MOD(k-1+6,6)
+                CALL xyz2lonlat(d%vertex(:,kk,i,j), lon(1),lat(1))
+                CALL xyz2lonlat(d%vertex(:,k, i,j), lon(2),lat(2))
+                cells%bnds_lon(:,n)=lon(:)*180./Pi
+                cells%bnds_lat(:,n)=lat(:)*180/Pi
+             END IF
+          END DO
+       END DO
+    END DO
+
+  END SUBROUTINE set_bounds_edge
+
+  SUBROUTINE set_bounds(domain_type, glo)
+    TYPE(t_domain), POINTER :: domain_type(:), d
+    LOGICAL :: glo
+    INTEGER :: ind
+    !$OMP BARRIER    
+    !$OMP MASTER
+    ! IF glo is .TRUE. we are dealing with the global mesh, otherwise with the local mesh
+    ! write_field uses the global mesh and may want halo cells
+    ! output_field (XIOS) uses the local mesh and uses only own cells
+    ! output_field uses edges while write_field uses only primal and dual cells
+    DO ind=1, SIZE(domain_type)
+       d=>domain_type(ind)
+       IF(glo) THEN ! global mesh / write_field
+          ! primal cell i,j is owned if d%assign_domain(i,j)==ind
+          CALL set_bounds_primal(d, d%primal_own, .FALSE., d%assign_domain==ind)
+          CALL set_bounds_primal(d, d%primal_all, .TRUE., d%assign_domain==ind)
+          CALL set_bounds_dual(d, d%dual_own)
+          CALL set_bounds_dual(d, d%dual_all)
+       ELSE  ! local mesh / XIOS
+          ! primal cell i,j is owned if d%own(i,j)==.TRUE.
+          CALL set_bounds_primal(d, d%primal_own, .FALSE., d%own)
+          CALL set_bounds_dual(d, d%dual_own)
+          CALL set_bounds_edge(ind, d, d%edge_own)
+       END IF
+    END DO
+    !$OMP END MASTER
+    !$OMP BARRIER
+  END SUBROUTINE set_bounds
+  
+END MODULE set_bounds_mod
Index: codes/icosagcm/devel/src/output/write_field.f90
===================================================================
--- codes/icosagcm/devel/src/output/write_field.f90	(revision 879)
+++ codes/icosagcm/devel/src/output/write_field.f90	(revision 880)
@@ -1,21 +1,7 @@
 module write_field_mod
   USE genmod
+  USE write_field_vars_mod
   IMPLICIT NONE
   PRIVATE  
-  INTEGER,SAVE :: ncprec
-  
-  TYPE ncvar
-    INTEGER :: size
-    INTEGER,POINTER :: nc_id(:)
-    INTEGER :: displ
-  END TYPE ncvar
-
-  INTEGER, PARAMETER :: MaxWriteField = 1000
-  INTEGER, DIMENSION(MaxWriteField),SAVE :: FieldId
-  TYPE(ncvar), dimension(MaxWriteField),SAVE :: FieldVarId
-  INTEGER, DIMENSION(MaxWriteField),SAVE :: FieldIndex
-  CHARACTER(len=255), DIMENSION(MaxWriteField) ::  FieldName 
-   
-  INTEGER,SAVE :: NbField = 0
   
   PUBLIC init_writeField, writefield, close_files
@@ -43,25 +29,5 @@
     END SUBROUTINE init_writeField
     
-    function GetFieldIndex(name)
-    implicit none
-      integer          :: GetFieldindex
-      character(len=*) :: name
-    
-      character(len=255) :: TrueName
-      integer            :: i
-       
-      
-      TrueName=TRIM(ADJUSTL(name))
-    
-      GetFieldIndex=-1
-      do i=1,NbField
-        if (TrueName==FieldName(i)) then
-          GetFieldIndex=i
-          exit
-        endif
-      enddo
-    end function GetFieldIndex
-    
-    SUBROUTINE Writefield(name_in,field,nind,once)
+   SUBROUTINE Writefield(name_in,field,nind,once)
     USE domain_mod
     USE field_mod
@@ -71,5 +37,4 @@
     USE netcdf_mod
     USE grid_param
-    IMPLICIT NONE  
      CHARACTER(LEN=*),INTENT(IN) :: name_in
       TYPE(t_field),POINTER :: field(:)
@@ -81,5 +46,4 @@
       
       IF (no_io) RETURN
-      IF (grid_type==grid_unst) RETURN
 
 !$OMP BARRIER
@@ -91,14 +55,18 @@
       END IF
 
-      IF (field(1)%ndim==2) THEN
-        CALL allocate_field_glo(field_glo,field(1)%field_type,field(1)%data_type)
-      ELSE IF (field(1)%ndim==3) THEN
-        CALL allocate_field_glo(field_glo,field(1)%field_type,field(1)%data_type,field(1)%dim3) 
-      ELSE IF (field(1)%ndim==4) THEN
-        CALL allocate_field_glo(field_glo,field(1)%field_type,field(1)%data_type,field(1)%dim3,field(1)%dim4)
-      ENDIF
-      
-      CALL gather_field(field,field_glo)
-      
+      IF (grid_type==grid_ico) THEN
+         IF (field(1)%ndim==2) THEN
+            CALL allocate_field_glo(field_glo,field(1)%field_type,field(1)%data_type)
+         ELSE IF (field(1)%ndim==3) THEN
+            CALL allocate_field_glo(field_glo,field(1)%field_type,field(1)%data_type,field(1)%dim3) 
+         ELSE IF (field(1)%ndim==4) THEN
+            CALL allocate_field_glo(field_glo,field(1)%field_type,field(1)%data_type,field(1)%dim3,field(1)%dim4)
+         ENDIF
+         
+         CALL gather_field(field,field_glo)
+      ELSE
+         field_glo => field ! FIXME unstructured + MPI
+      END IF
+
       IF (mpi_rank==0) THEN
         IF (PRESENT(nind)) THEN
@@ -109,5 +77,5 @@
       ENDIF
       
-      CALL deallocate_field_glo(field_glo)
+      IF(grid_type == grid_ico) CALL deallocate_field_glo(field_glo)
 !$OMP END MASTER
 !$OMP BARRIER
@@ -115,255 +83,27 @@
    END SUBROUTINE writefield
    
-!   SUBROUTINE Writefield(name_in,field,nind)
-!   USE netcdf
-!   USE domain_mod
-!   use field_mod
-!   USE dimensions
-!   USE geometry
-!   IMPLICIT NONE
-!     CHARACTER(LEN=*),INTENT(IN) :: name_in
-!     TYPE(t_field),POINTER :: field(:)
-!     INTEGER,OPTIONAL,INTENT(IN) :: nind
-!     REAL(r8),ALLOCATABLE :: field_val2d(:)
-!     REAL(r8),ALLOCATABLE :: field_val3d(:,:)
-!     REAL(r8),ALLOCATABLE :: field_val4d(:,:,:)
-!     TYPE(t_domain),POINTER :: d
-!     INTEGER :: Index
-!     INTEGER :: ind,i,j,k,n,ncell,q
-!     INTEGER :: iie,jje,iin,jjn
-!     INTEGER :: status
-!     CHARACTER(len=255) :: name
-!     CHARACTER(len=255) :: str_ind
-!     INTEGER :: ind_b,ind_e
-!     INTEGER :: halo_size
-!     LOGICAL :: single
-!     
-!     
-!     name=TRIM(ADJUSTL(name_in))
-
-!     IF (PRESENT(nind)) THEN
-!       name=TRIM(name)//"_"//TRIM(int2str(nind))
-!       PRINT *,"NAME",nind,int2str(nind),name
-!       ind_b=nind
-!       ind_e=nind
-!       halo_size=1
-!       single=.TRUE.
-!     ELSE
-!       ind_b=1
-!       ind_e=ndomain
-!       halo_size=0
-!       single=.FALSE.
-!     ENDIF      
-
-!     Index=GetFieldIndex(name)
-!     if (Index==-1) then
-!       call create_header(name,field,nind)
-!	Index=GetFieldIndex(name)
-!     else
-!       FieldIndex(Index)=FieldIndex(Index)+1.
-!     endif
-!     
-!     IF (Field(ind_b)%field_type==field_T) THEN
-!       ncell=1
-!       DO ind=ind_b,ind_e
-!         d=>domain(ind)
-!         IF (Field(ind)%field_type/=field_T) THEN
-!           PRINT *,"Writefield, grille non geree"
-!           RETURN
-!         ENDIF
-
-!       n=0
-!       DO j=d%jj_begin-halo_size,d%jj_end+halo_size
-!         DO i=d%ii_begin-halo_size,d%ii_end+halo_size
-!           IF (d%own(i,j) .OR. single) n=n+1
-!         ENDDO
-!       ENDDO
-
-!       IF (field(ind)%ndim==2) THEN
-!         ALLOCATE(Field_val2d(n))
-!         n=0
-!         DO j=d%jj_begin-halo_size,d%jj_end+halo_size
-!           DO i=d%ii_begin-halo_size,d%ii_end+halo_size
-!             k=d%iim*(j-1)+i
-!             IF (d%own(i,j) .OR. single) THEN
-!               n=n+1
-!               Field_val2d(n)=field(ind)%rval2d(k)
-!             ENDIF
-!           ENDDO
-!         ENDDO
-!         status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val2d,  &
-!                             start=(/ ncell,FieldIndex(Index) /),count=(/n,1 /))
-!         DEALLOCATE(field_val2d)
-!       ELSE IF (field(ind)%ndim==3) THEN
-!         ALLOCATE(Field_val3d(n,size(field(ind)%rval3d,2)))
-!         n=0
-!         DO j=d%jj_begin-halo_size,d%jj_end+halo_size
-!           DO i=d%ii_begin-halo_size,d%ii_end+halo_size
-!             k=d%iim*(j-1)+i
-!             IF (d%own(i,j) .OR. single) THEN
-!               n=n+1
-!               Field_val3d(n,:)=field(ind)%rval3d(k,:)
-!             ENDIF
-!           ENDDO
-!         ENDDO
-!          status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ ncell,1,FieldIndex(Index) /),   &
-!                              count=(/n,size(field(1)%rval3d,2),1 /))
-!         DEALLOCATE(field_val3d)
-!       ELSE IF (field(1)%ndim==4) THEN
-
-!         DO q=1,FieldVarId(index)%size
-!         
-!           ALLOCATE(Field_val3d(n,size(field(ind)%rval4d,2)))
-!           n=0
-!           DO j=d%jj_begin-halo_size,d%jj_end+halo_size
-!             DO i=d%ii_begin-halo_size,d%ii_end+halo_size
-!               k=d%iim*(j-1)+i
-!               IF (d%own(i,j) .OR. single) THEN
-!                 n=n+1
-!                 Field_val3d(n,:)=field(ind)%rval4d(k,:,q)
-!               ENDIF
-!             ENDDO
-!           ENDDO
-!           status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ ncell,1,FieldIndex(Index) /),   &
-!                              count=(/n,size(field(1)%rval4d,2),1 /))
-!           DEALLOCATE(field_val3d)
-!         ENDDO         
-!       ENDIF
-!       
-!        PRINT *,NF90_STRERROR(status)
-!       ncell=ncell+n
-
-!    ENDDO
-!    
-!    ELSE IF (Field(ind_b)%field_type==field_Z) THEN
-!       ncell=1
-!       n=0
-!       DO ind=ind_b,ind_e
-!         d=>domain(ind)
-!         CALL swap_geometry(ind)
-!         CALL swap_dimensions(ind)
-! 
-!         n=0
-!         DO j=jj_begin+1,jj_end
-!           DO i=ii_begin,ii_end-1
-!             n=n+1
-!           ENDDO
-!         ENDDO
-
-!         DO j=jj_begin,jj_end-1
-!           DO i=ii_begin+1,ii_end
-!             n=n+1
-!           ENDDO
-!         ENDDO
-
-!       IF (field(ind)%ndim==2) THEN
-!         ALLOCATE(Field_val2d(n))
-
-!         n=0
-!         DO j=jj_begin+1,jj_end
-!           DO i=ii_begin,ii_end-1
-!             n=n+1
-!             k=iim*(j-1)+i
-!             Field_val2d(n)=field(ind)%rval2d(k+z_down)
-!           ENDDO
-!         ENDDO
-
-!         DO j=jj_begin,jj_end-1
-!           DO i=ii_begin+1,ii_end
-!             n=n+1
-!             k=iim*(j-1)+i
-!             Field_val2d(n)=field(ind)%rval2d(k+z_up)
-!           ENDDO
-!         ENDDO
-
-!         status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),                       &
-!                             Field_val2d,start=(/ ncell,FieldIndex(Index) /),count=(/n,1 /))
-!         DEALLOCATE(field_val2d)
-
-!       ELSE IF (field(ind)%ndim==3) THEN
-!         ALLOCATE(Field_val3d(n,size(field(ind)%rval3d,2)))
-!         n=0
-!         DO j=jj_begin+1,jj_end
-!           DO i=ii_begin,ii_end-1
-!             n=n+1
-!             k=iim*(j-1)+i
-!             Field_val3d(n,:)=field(ind)%rval3d(k+z_down,:)
-!           ENDDO
-!         ENDDO
-
-!         DO j=jj_begin,jj_end-1
-!           DO i=ii_begin+1,ii_end
-!             n=n+1
-!             k=iim*(j-1)+i
-!             Field_val3d(n,:)=field(ind)%rval3d(k+z_up,:)
-!           ENDDO
-!         ENDDO
-!          status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ ncell,1,FieldIndex(Index) /),   &
-!                              count=(/n,size(field(1)%rval3d,2),1 /))
-!         DEALLOCATE(field_val3d)
-!       ELSE IF (field(1)%ndim==4) THEN
-
-!         DO q=1,FieldVarId(index)%size
-!           ALLOCATE(Field_val3d(n,size(field(ind)%rval4d,2)))
-!           n=0
-!           DO j=jj_begin+1,jj_end
-!             DO i=ii_begin,ii_end-1
-!               n=n+1
-!               k=iim*(j-1)+i
-!               Field_val3d(n,:)=field(ind)%rval4d(k+z_down,:,q)
-!             ENDDO
-!           ENDDO
-
-!           DO j=jj_begin,jj_end-1
-!             DO i=ii_begin+1,ii_end
-!               n=n+1
-!               k=iim*(j-1)+i
-!               Field_val3d(n,:)=field(ind)%rval4d(k+z_up,:,q)
-!             ENDDO
-!           ENDDO
-
-!           status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ ncell,1,1,FieldIndex(Index) /),   &
-!                               count=(/n,size(field(1)%rval4d,2),1 /))
-!           DEALLOCATE(field_val3d)
-!         ENDDO
-!       ENDIF
-!       
-!        PRINT *,NF90_STRERROR(status)
-!       ncell=ncell+n
-
-!    ENDDO
-!    
-!    ENDIF
-!    status=NF90_SYNC(FieldId(Index))
-!     
-!   END SUBROUTINE Writefield
-
-
-    SUBROUTINE Writefield_gen(name_in, field, domain_type, ind_b_in, ind_e_in,once )
-    USE netcdf_mod
-    USE domain_mod
-    USE field_mod
-    IMPLICIT NONE
-      CHARACTER(LEN=*),INTENT(IN) :: name_in
-      TYPE(t_field), POINTER :: field(:)
-      TYPE(t_domain),INTENT(IN),TARGET :: domain_type(:)
-      INTEGER,OPTIONAL,INTENT(IN) :: ind_b_in
-      INTEGER,OPTIONAL,INTENT(IN) :: ind_e_in
-      REAL(r8),ALLOCATABLE :: field_val2d(:)
-      REAL(r8),ALLOCATABLE :: field_val3d(:,:)
-      REAL(r8),ALLOCATABLE :: field_val4d(:,:,:)
-      LOGICAL, INTENT(IN) :: once
-      TYPE(t_domain),POINTER :: d
-      INTEGER :: Index
-      INTEGER :: ind,i,j,k,n,ncell,q
-      INTEGER :: iie,jje,iin,jjn
-      INTEGER :: status
-      CHARACTER(len=255) :: name
-      CHARACTER(len=255) :: str_ind
-      INTEGER :: ind_b,ind_e
-      INTEGER :: halo_size
-      LOGICAL :: single
-      
-      name=TRIM(ADJUSTL(name_in))
+   SUBROUTINE Writefield_gen(name_in, field, domain_type, ind_b_in, ind_e_in,once )
+     USE netcdf_mod
+     USE domain_mod
+     USE field_mod
+     CHARACTER(LEN=*),INTENT(IN) :: name_in
+     TYPE(t_field), POINTER :: field(:)
+     TYPE(t_domain),INTENT(IN),TARGET :: domain_type(:)
+     INTEGER,OPTIONAL,INTENT(IN) :: ind_b_in
+     INTEGER,OPTIONAL,INTENT(IN) :: ind_e_in
+     LOGICAL, INTENT(IN) :: once
+     ! local variables
+     TYPE(t_cellset), POINTER :: cells(:)
+     TYPE(t_domain),POINTER :: d
+     LOGICAL :: single
+     INTEGER :: ind_b, ind_e, ind
+     CHARACTER(len=255) :: name
+     INTEGER :: index, ncell, nvert, n
+     ! for embedded routines
+     REAL(r8),ALLOCATABLE :: field_val2d(:)
+     REAL(r8),ALLOCATABLE :: field_val3d(:,:)
+     INTEGER :: status, n_begin, ij, q, dim3
+
+     name=TRIM(ADJUSTL(name_in))
 
       IF (PRESENT(ind_b_in) .AND. .NOT. PRESENT(ind_e_in)) THEN
@@ -371,5 +111,4 @@
         ind_b=ind_b_in
         ind_e=ind_b_in
-        halo_size=1
         single=.TRUE.
       ELSE IF (.NOT. PRESENT(ind_b_in) .AND. PRESENT(ind_e_in)) THEN
@@ -377,1412 +116,337 @@
         ind_b=ind_e_in
         ind_e=ind_e_in
-        halo_size=1
         single=.TRUE.
       ELSE IF ( PRESENT(ind_b_in) .AND. PRESENT(ind_e_in)) THEN
         ind_b=ind_b_in
         ind_e=ind_e_in
-        halo_size=0
         single=.FALSE.
       ELSE
-        halo_size=0
         ind_b=1
         ind_e=ndomain
         single=.FALSE.
       ENDIF      
-      
-      Index=GetFieldIndex(name)
-      if (Index==-1) then
-        call create_header_gen(name_in,field,domain_type,ind_b_in,ind_e_in,once)
-	Index=GetFieldIndex(name)
+
+      index=GetFieldIndex(name)
+      if (index==-1) then
+         call create_header_gen(name_in,field,domain_type,ind_b_in,ind_e_in,once)
+         index=GetFieldIndex(name)
       else
-        FieldIndex(Index)=FieldIndex(Index)+1.
+         FieldIndex(index)=FieldIndex(index)+1.
       endif
-      
-      IF (Field(ind_b)%field_type==field_T) THEN
-
-        ncell=0
-        DO ind=ind_b,ind_e
-          d=>domain_type(ind)
-          DO j=d%jj_begin-halo_size,d%jj_end+halo_size
-            DO i=d%ii_begin-halo_size,d%ii_end+halo_size
-              IF (d%assign_domain(i,j)==ind .OR. single) ncell=ncell+1
-            ENDDO
-          ENDDO
-        ENDDO
-        
-        IF (field(1)%ndim==2) THEN
-          ALLOCATE(Field_val2d(ncell))
-          n=0
+
+      SELECT CASE(Field(ind_b)%field_type)
+      CASE(field_T)
+         IF(single) THEN ! include halos
+            cells => mesh_glo%primal_all
+         ELSE
+            cells => mesh_glo%primal_own
+         END IF
+      CASE(field_Z)
+         IF(single) THEN ! include halos
+            cells => mesh_glo%dual_all
+         ELSE
+            cells => mesh_glo%dual_own
+         END IF
+      END SELECT
+
+      ncell=0
+      DO ind=ind_b,ind_e
+         nvert=SIZE(cells(ind)%bnds_lon,1)
+         ncell = ncell + cells(ind)%ncell
+      END DO
+
+      SELECT CASE(field(1)%ndim)
+      CASE(2)
+         CALL write_2d
+      CASE(3)
+         CALL write_3d
+      CASE(4)          
+         CALL write_4d
+      END SELECT
+      
+      status=NF90_SYNC(FieldId(index))
+      
+   CONTAINS
+    
+     SUBROUTINE write_2d
+       ALLOCATE(Field_val2d(ncell))
+       n_begin=0  
+       DO ind=ind_b,ind_e
+          DO n=1, cells(ind)%ncell
+             ij=cells(ind)%ij(n)
+             PRINT *, 'write_2d :', ind, n, n_begin+n, ij 
+             field_val2d(n_begin+n) = field(ind)%rval2d(ij)
+          END DO
+          n_begin = n_begin + cells(ind)%ncell
+       END DO
+       IF (once) THEN 
+          status=NF90_PUT_VAR(FieldId(index), FieldVarId(index)%nc_id(1), &
+               Field_val2d, start=(/ 1 /),count=(/ncell /))
+       ELSE
+          status=NF90_PUT_VAR(FieldId(Index), FieldVarId(index)%nc_id(1), &
+               Field_val2d,  start=(/ 1,FieldIndex(Index) /),count=(/ncell,1 /))
+       ENDIF
+       DEALLOCATE(field_val2d)
+     END SUBROUTINE write_2d
+
+     SUBROUTINE write_3d
+       dim3 = SIZE(field(1)%rval3d,2)
+       ALLOCATE(field_val3d(ncell,dim3))
+       n_begin=0  
+       DO ind=ind_b,ind_e
+          DO n=1, cells(ind)%ncell
+             ij=cells(ind)%ij(n)
+             field_val3d(n_begin+n,:) = field(ind)%rval3d(ij,:)
+          END DO
+       END DO
+       IF (once) THEN 
+          status=NF90_PUT_VAR(FieldId(Index), FieldVarId(index)%nc_id(1), &
+               field_val3d, start=(/ 1,1 /), &
+               count=(/ncell,dim3 /))
+       ELSE
+          status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1), &
+               field_val3d, start=(/ 1,1,FieldIndex(Index) /),   &
+               count=(/ncell,size(field(1)%rval3d,2),1 /))
+       ENDIF
+       DEALLOCATE(field_val3d)
+     END SUBROUTINE write_3d
+
+     SUBROUTINE write_4d
+       dim3 = SIZE(field(1)%rval4d,2)
+       ALLOCATE(field_val3d(ncell,dim3))
+       DO q=1,FieldVarId(index)%size
+          n_begin=0  
           DO ind=ind_b,ind_e
-            d=>domain_type(ind)
-            DO j=d%jj_begin-halo_size,d%jj_end+halo_size
-              DO i=d%ii_begin-halo_size,d%ii_end+halo_size
-                k=d%iim*(j-1)+i
-                IF (d%assign_domain(i,j)==ind .OR. single) THEN
-                  n=n+1
-                  Field_val2d(n)=field(ind)%rval2d(k)
-                ENDIF
-              ENDDO
-            ENDDO
-          ENDDO
-
+             DO n=1, cells(ind)%ncell
+                ij=cells(ind)%ij(n)
+                field_val3d(n_begin+n,:) = field(ind)%rval4d(ij,:,q)
+             END DO
+          END DO
           IF (once) THEN 
-            status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val2d,  &
-                               start=(/ 1 /),count=(/ncell /))
+             status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q), &
+                  field_val3d, start=(/ 1,1 /),   &
+                  count=(/ncell,dim3,1 /))
           ELSE
-            status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val2d,  &
-                               start=(/ 1,FieldIndex(Index) /),count=(/ncell,1 /))
+             status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q), &
+                  field_val3d, start=(/ 1,1,FieldIndex(Index) /),   &
+                  count=(/ncell,dim3,1 /))
           ENDIF
-          
-          DEALLOCATE(field_val2d)
-        ELSE IF (field(1)%ndim==3) THEN
-          ALLOCATE(Field_val3d(ncell,size(field(1)%rval3d,2)))
-          n=0
-          DO ind=ind_b,ind_e
-            d=>domain_type(ind)
-            DO j=d%jj_begin-halo_size,d%jj_end+halo_size
-              DO i=d%ii_begin-halo_size,d%ii_end+halo_size
-                k=d%iim*(j-1)+i
-                IF (d%assign_domain(i,j)==ind .OR. single) THEN
-                  n=n+1
-                  Field_val3d(n,:)=field(ind)%rval3d(k,:)
-                ENDIF
-              ENDDO
-            ENDDO
-          ENDDO
-
-          PRINT *, 'Writefield ', TRIM(name), MAXVAL(Field_val3d(:,1)), MINVAL(Field_val3d(:,1)) ! FIXME
-
-          IF (once) THEN 
-            status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ 1,1 /),   &
-                                 count=(/ncell,size(field(1)%rval3d,2) /))
-          ELSE
-             status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ 1,1,FieldIndex(Index) /),   &
-                                 count=(/ncell,size(field(1)%rval3d,2),1 /))
-          ENDIF
-                                
-          DEALLOCATE(field_val3d)
-        ELSE IF (field(1)%ndim==4) THEN
-
-          DO q=1,FieldVarId(index)%size
-          
-            ALLOCATE(Field_val3d(ncell,size(field(1)%rval4d,2)))
-            n=0
-            DO ind=ind_b,ind_e
-              d=>domain_type(ind)
-              DO j=d%jj_begin-halo_size,d%jj_end+halo_size
-                DO i=d%ii_begin-halo_size,d%ii_end+halo_size
-                  k=d%iim*(j-1)+i
-                  IF (d%assign_domain(i,j)==ind .OR. single) THEN
-                    n=n+1
-                    Field_val3d(n,:)=field(ind)%rval4d(k,:,q)
-                  ENDIF
-                ENDDO
-              ENDDO
-            ENDDO
-            
-            IF (once) THEN 
-              status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ 1,1 /),   &
-                                 count=(/ncell,size(field(1)%rval4d,2),1 /))
-            ELSE
-              status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ 1,1,FieldIndex(Index) /),   &
-                                 count=(/ncell,size(field(1)%rval4d,2),1 /))
-            ENDIF
-            DEALLOCATE(field_val3d)
-          ENDDO         
-        ENDIF
-        
-     
-     ELSE IF (Field(ind_b)%field_type==field_Z) THEN
-
-        ncell=0
-        DO ind=ind_b,ind_e
-          d=>domain_type(ind)
-  
-          DO j=d%jj_begin+1,d%jj_end
-            DO i=d%ii_begin,d%ii_end-1
-              ncell=ncell+1
-            ENDDO
-          ENDDO
-
-          DO j=d%jj_begin,d%jj_end-1
-            DO i=d%ii_begin+1,d%ii_end
-              ncell=ncell+1
-            ENDDO
-          ENDDO
-        ENDDO
-
-        IF (field(1)%ndim==2) THEN
-          ALLOCATE(Field_val2d(ncell))
-
-          n=0
-          
-          DO ind=ind_b,ind_e
-            d=>domain_type(ind)
-            DO j=d%jj_begin+1,d%jj_end
-              DO i=d%ii_begin,d%ii_end-1
-                n=n+1
-                k=d%iim*(j-1)+i
-                Field_val2d(n)=field(ind)%rval2d(k+d%z_down)
-              ENDDO
-            ENDDO
-
-            DO j=d%jj_begin,d%jj_end-1
-              DO i=d%ii_begin+1,d%ii_end
-                n=n+1
-                k=d%iim*(j-1)+i
-                Field_val2d(n)=field(ind)%rval2d(k+d%z_up)
-              ENDDO
-            ENDDO
-          ENDDO
-          
-          IF (once) THEN 
-            status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),                       &
-                                Field_val2d,start=(/ 1 /),count=(/ncell /))
-           ELSE
-            status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),                       &
-                                Field_val2d,start=(/ 1,FieldIndex(Index) /),count=(/ncell,1 /))
-          ENDIF
-          DEALLOCATE(field_val2d)
-
-        ELSE IF (field(1)%ndim==3) THEN
-          ALLOCATE(Field_val3d(ncell,size(field(1)%rval3d,2)))
-          n=0
-          DO ind=ind_b,ind_e
-            d=>domain_type(ind)
-            DO j=d%jj_begin+1,d%jj_end
-              DO i=d%ii_begin,d%ii_end-1
-                n=n+1
-                k=d%iim*(j-1)+i
-                Field_val3d(n,:)=field(ind)%rval3d(k+d%z_down,:)
-              ENDDO
-            ENDDO
-
-            DO j=d%jj_begin,d%jj_end-1
-              DO i=d%ii_begin+1,d%ii_end
-                n=n+1
-                k=d%iim*(j-1)+i
-                Field_val3d(n,:)=field(ind)%rval3d(k+d%z_up,:)
-              ENDDO
-            ENDDO
-          ENDDO
-          
-          IF (once) THEN 
-            status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ 1,1 /),   &
-                                count=(/ncell,size(field(1)%rval3d,2) /))
-          ELSE
-            status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,start=(/ 1,1,FieldIndex(Index) /),   &
-                                count=(/ncell,size(field(1)%rval3d,2),1 /))
-          ENDIF
-          DEALLOCATE(field_val3d)
-        
-        ELSE IF (field(1)%ndim==4) THEN
-
-          DO q=1,FieldVarId(index)%size
-            ALLOCATE(Field_val3d(ncell,size(field(1)%rval4d,2)))
-            n=0
-            DO ind=ind_b,ind_e
-              d=>domain_type(ind)
-              DO j=d%jj_begin+1,d%jj_end
-                DO i=d%ii_begin,d%ii_end-1
-                  n=n+1
-                  k=d%iim*(j-1)+i
-                  Field_val3d(n,:)=field(ind)%rval4d(k+d%z_down,:,q)
-                ENDDO
-              ENDDO
-
-              DO j=d%jj_begin,d%jj_end-1
-                DO i=d%ii_begin+1,d%ii_end
-                  n=n+1
-                  k=d%iim*(j-1)+i
-                  Field_val3d(n,:)=field(ind)%rval4d(k+d%z_up,:,q)
-                ENDDO
-              ENDDO
-            ENDDO
-
-            IF (once) THEN 
-              status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ 1,1,1 /),   &
-                                  count=(/ncell,size(field(1)%rval4d,2) /))
-            ELSE
-              status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,start=(/ 1,1,1,FieldIndex(Index) /),   &
-                                  count=(/ncell,size(field(1)%rval4d,2),1 /))
-            ENDIF
-            DEALLOCATE(field_val3d)
-          ENDDO
-        ENDIF
-     
-     ENDIF
-     status=NF90_SYNC(FieldId(Index))
-      
-    END SUBROUTINE Writefield_gen
-
-      
-
-    SUBROUTINE Writefield_mpi(name_in,field,nind)
-    USE netcdf_mod
-    USE domain_mod
-    use field_mod
-    USE dimensions
-    USE geometry
-    IMPLICIT NONE
-      CHARACTER(LEN=*),INTENT(IN) :: name_in
-      TYPE(t_field),POINTER :: field(:)
-      INTEGER,OPTIONAL,INTENT(IN) :: nind
-      REAL(r8),ALLOCATABLE :: field_val2d(:)
-      REAL(r8),ALLOCATABLE :: field_val3d(:,:)
-      REAL(r8),ALLOCATABLE :: field_val4d(:,:,:)
-      TYPE(t_domain),POINTER :: d
-      INTEGER :: Index
-      INTEGER :: ind,i,j,l,k,n,ncell,q
-      INTEGER :: iie,jje,iin,jjn
-      INTEGER :: status
-      CHARACTER(len=255) :: name
-      CHARACTER(len=255) :: str_ind
-      INTEGER :: ind_b,ind_e
-      INTEGER :: halo_size
-      LOGICAL :: single
-      INTEGER :: displ
-      
-      
-      name=TRIM(ADJUSTL(name_in))
-
-      IF (PRESENT(nind)) THEN
-        name=TRIM(name)//"_"//TRIM(int2str(nind))
-        PRINT *,"NAME",nind,int2str(nind),name
-        ind_b=nind
-        ind_e=nind
-        halo_size=1
-        single=.TRUE.
-      ELSE
-        ind_b=1
-        ind_e=ndomain
-        halo_size=0
-        single=.FALSE.
-      ENDIF      
-
-      Index=GetFieldIndex(name)
-      if (Index==-1) then
-        call create_header_mpi(name,field,nind)
-	Index=GetFieldIndex(name)
-      else
-        FieldIndex(Index)=FieldIndex(Index)+1.
-      endif
-      
-      IF (Field(ind_b)%field_type==field_T) THEN
-        ncell=1
-        DO ind=ind_b,ind_e
-          d=>domain(ind)
-          IF (Field(ind)%field_type/=field_T) THEN
-            PRINT *,"Writefield, grille non geree"
-            RETURN
-          ENDIF
-
-        n=0
-        DO j=d%jj_begin-halo_size,d%jj_end+halo_size
-          DO i=d%ii_begin-halo_size,d%ii_end+halo_size
-            IF (d%own(i,j) .OR. single) n=n+1
-          ENDDO
-        ENDDO
-
-        displ=FieldVarId(index)%displ
-        
-        IF (field(ind)%ndim==2) THEN
-          ALLOCATE(Field_val2d(n))
-          n=0
-          DO j=d%jj_begin-halo_size,d%jj_end+halo_size
-            DO i=d%ii_begin-halo_size,d%ii_end+halo_size
-              k=d%iim*(j-1)+i
-              IF (d%own(i,j) .OR. single) THEN
-                n=n+1
-                Field_val2d(n)=field(ind)%rval2d(k)
-              ENDIF
-            ENDDO
-          ENDDO
-          status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val2d,  &
-                              start=(/ displ+ncell,FieldIndex(Index) /),count=(/n,1 /))
-          DEALLOCATE(field_val2d)
-        ELSE IF (field(ind)%ndim==3) THEN
-          ALLOCATE(Field_val3d(n,size(field(ind)%rval3d,2)))
-          n=0
-          DO j=d%jj_begin-halo_size,d%jj_end+halo_size
-            DO i=d%ii_begin-halo_size,d%ii_end+halo_size
-              k=d%iim*(j-1)+i
-              IF (d%own(i,j) .OR. single) THEN
-                n=n+1
-                Field_val3d(n,:)=field(ind)%rval3d(k,:)
-              ENDIF
-            ENDDO
-          ENDDO
-          status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,  &
-                                start=(/ displ+ncell,1,FieldIndex(Index) /), count=(/n,size(field(ind)%rval3d,2),1 /))
-          DEALLOCATE(field_val3d)
-
-        ELSE IF (field(1)%ndim==4) THEN
-
-          DO q=1,FieldVarId(index)%size
-          
-            ALLOCATE(Field_val3d(n,size(field(ind)%rval4d,2)))
-            n=0
-            DO j=d%jj_begin-halo_size,d%jj_end+halo_size
-              DO i=d%ii_begin-halo_size,d%ii_end+halo_size
-                k=d%iim*(j-1)+i
-                IF (d%own(i,j) .OR. single) THEN
-                  n=n+1
-                  Field_val3d(n,:)=field(ind)%rval4d(k,:,q)
-                ENDIF
-              ENDDO
-             ENDDO
-
-            status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d(:,l), &
-                                start=(/ displ+ncell,l,FieldIndex(Index) /), count=(/n,size(field(ind)%rval4d,2),1 /))
-            DEALLOCATE(field_val3d)
-          ENDDO         
-        ENDIF
-        
-        ncell=ncell+n
-
-     ENDDO
-     
-     ELSE IF (Field(ind_b)%field_type==field_Z) THEN
-        ncell=1
-        n=0
-        DO ind=ind_b,ind_e
-          d=>domain(ind)
-          CALL swap_geometry(ind)
-          CALL swap_dimensions(ind)
-  
-          n=0
-          DO j=jj_begin+1,jj_end
-            DO i=ii_begin,ii_end-1
-              n=n+1
-            ENDDO
-          ENDDO
-
-          DO j=jj_begin,jj_end-1
-            DO i=ii_begin+1,ii_end
-              n=n+1
-            ENDDO
-          ENDDO
-
-        displ=FieldVarId(index)%displ
-
-        IF (field(ind)%ndim==2) THEN
-          ALLOCATE(Field_val2d(n))
-
-          n=0
-          DO j=jj_begin+1,jj_end
-            DO i=ii_begin,ii_end-1
-              n=n+1
-              k=iim*(j-1)+i
-              Field_val2d(n)=field(ind)%rval2d(k+z_down)
-            ENDDO
-          ENDDO
-
-          DO j=jj_begin,jj_end-1
-            DO i=ii_begin+1,ii_end
-              n=n+1
-              k=iim*(j-1)+i
-              Field_val2d(n)=field(ind)%rval2d(k+z_up)
-            ENDDO
-          ENDDO
-
-          status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),                       &
-                              Field_val2d,start=(/ displ+ncell,FieldIndex(Index) /),count=(/n,1 /))
-          DEALLOCATE(field_val2d)
-
-        ELSE IF (field(ind)%ndim==3) THEN
-          ALLOCATE(Field_val3d(n,size(field(ind)%rval3d,2)))
-          n=0
-          DO j=jj_begin+1,jj_end
-            DO i=ii_begin,ii_end-1
-              n=n+1
-              k=iim*(j-1)+i
-              Field_val3d(n,:)=field(ind)%rval3d(k+z_down,:)
-            ENDDO
-          ENDDO
-
-          DO j=jj_begin,jj_end-1
-            DO i=ii_begin+1,ii_end
-              n=n+1
-              k=iim*(j-1)+i
-              Field_val3d(n,:)=field(ind)%rval3d(k+z_up,:)
-            ENDDO
-          ENDDO
-
-           status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,    &
-                               start=(/ displ+ncell,1,FieldIndex(Index) /), count=(/n,size(field(ind)%rval3d,2),1 /))
-          DEALLOCATE(field_val3d)
-
-        ELSE IF (field(1)%ndim==4) THEN
-
-          DO q=1,FieldVarId(index)%size
-            ALLOCATE(Field_val3d(n,size(field(ind)%rval4d,2)))
-            n=0
-            DO j=jj_begin+1,jj_end
-              DO i=ii_begin,ii_end-1
-                n=n+1
-                k=iim*(j-1)+i
-                Field_val3d(n,:)=field(ind)%rval4d(k+z_down,:,q)
-              ENDDO
-            ENDDO
-
-            DO j=jj_begin,jj_end-1
-              DO i=ii_begin+1,ii_end
-                n=n+1
-                k=iim*(j-1)+i
-                Field_val3d(n,:)=field(ind)%rval4d(k+z_up,:,q)
-              ENDDO
-            ENDDO
-
-            status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,  &
-                                start=(/ displ+ncell,1,FieldIndex(Index) /), count=(/n,size(field(ind)%rval4d,2),1 /))
-            DEALLOCATE(field_val3d)
-          ENDDO
-        ENDIF
-        
-        ncell=ncell+n
-
-     ENDDO
-     
-     ENDIF
-     status=NF90_SYNC(FieldId(Index))
-      
-    END SUBROUTINE Writefield_mpi
-    
-    
-           
-!   SUBROUTINE Create_header(name,field,nind)
-!   USE netcdf
-!   USE field_mod
-!   USE domain_mod
-!   USE spherical_geom_mod
-!   USE dimensions
-!   USE geometry
-!   IMPLICIT NONE
-!     CHARACTER(LEN=*) :: name
-!     TYPE(t_field),POINTER :: field(:)
-!     INTEGER,OPTIONAL,INTENT(IN) :: nind
-!     INTEGER :: ncell
-!     INTEGER :: nvert
-!     REAL(rstd),ALLOCATABLE :: lon(:),lat(:),bounds_lon(:,:),bounds_lat(:,:)
-!     TYPE(t_domain),POINTER :: d
-!     INTEGER :: nvertId,ncid,lonId,latId,bounds_lonId,bounds_latId,timeId,ncellId
-!     INTEGER :: dim3id,dim4id
-!     INTEGER :: status
-!     INTEGER :: ind,i,j,k,n,q
-!     INTEGER :: iie,jje,iin,jjn
-!     INTEGER :: ind_b,ind_e
-!     INTEGER :: halo_size
-!     LOGICAL :: single 
-!     INTEGER :: nij
-!         
-!     NbField=NbField+1
-!     FieldName(NbField)=TRIM(ADJUSTL(name))
-!     FieldIndex(NbField)=1
-!     
-!     IF (PRESENT(nind)) THEN
-!       ind_b=nind
-!       ind_e=nind
-!       halo_size=1
-!       single=.TRUE.
-!     ELSE
-!       ind_b=1
-!       ind_e=ndomain
-!       halo_size=0
-!       single=.FALSE.
-!     ENDIF
-!     
-!     ncell=0
-!     
-!     IF (Field(ind_b)%field_type==field_T) THEN
-!       nvert=6
-!       
-!       DO ind=ind_b,ind_e
-!         d=>domain(ind)
-!       
-!         DO j=d%jj_begin-halo_size,d%jj_end+halo_size
-!           DO i=d%ii_begin-halo_size,d%ii_end+halo_size
-!             IF (single .OR. domain(ind)%own(i,j)) ncell=ncell+1
-!           ENDDO
-!         ENDDO
-
-!       END DO
-!     
-!       status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid)
-!       FieldId(NbField)=ncid
-!       status = NF90_DEF_DIM(ncid,'cell',ncell,ncellId)
-!       status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid)
-
-!       IF (Field(ind_b)%ndim==2)  THEN
-!         FieldVarId(NbField)%size=1
-!         ALLOCATE(FieldVarId(NbField)%nc_id(1))
-!       ELSE IF (Field(ind_b)%ndim==3)  THEN
-!         FieldVarId(NbField)%size=1
-!         ALLOCATE(FieldVarId(NbField)%nc_id(1))
-!         status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval3d,2),dim3id)
-!       ELSE IF (Field(1)%ndim==4) THEN
-!         FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3)
-!         ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size))
-!         status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval4d,2),dim3id)
-!          status = NF90_DEF_DIM(ncid,'Q',size(field(ind_b)%rval4d,3),dim4id)
-!       ENDIF
-!     
-!       status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId)
-!     
-!       status = NF90_DEF_VAR(ncid,'lon',NF90_DOUBLE,(/ ncellId /),lonId)
-!       status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude")
-!       status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east")
-!       status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon")
-!       status = NF90_DEF_VAR(ncid,'lat',NF90_DOUBLE,(/ ncellId /),latId)
-!       status = NF90_PUT_ATT(ncid,latId,"long_name","latitude")
-!       status = NF90_PUT_ATT(ncid,latId,"units","degrees_north")
-!       status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat")
-!       status = NF90_DEF_VAR(ncid,'bounds_lon',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId)
-!       status = NF90_DEF_VAR(ncid,'bounds_lat',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId)
-
-!       IF (Field(ind_b)%ndim==2) THEN
-!         status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1))
-!         status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat")
-!       ELSE IF (Field(ind_b)%ndim==3) THEN
-!         status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1))
-!         status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat")
-!       ELSE IF (Field(ind_b)%ndim==4) THEN
-!         DO i=1,FieldVarId(NbField)%size
-!           status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name))//int2str(i),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),  &
-!                                 FieldVarId(NbField)%nc_id(i))
-!           status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(i),"coordinates","lon lat")
-!         ENDDO        
-!       ENDIF 
-! 
-!     
-!       status = NF90_ENDDEF(ncid)      
-
-!       ncell=1
-!       DO ind=ind_b,ind_e
-!         d=>domain(ind)
-! 
-!         n=0
-!         DO j=d%jj_begin-halo_size,d%jj_end+halo_size
-!           DO i=d%ii_begin-halo_size,d%ii_end+halo_size
-!             IF (single .OR. d%own(i,j)) n=n+1
-!           ENDDO
-!         ENDDO
-!         
-!        ALLOCATE(lon(n),lat(n),bounds_lon(0:nvert-1,n),bounds_lat(0:nvert-1,n))
-!         
-!         n=0  
-!         DO j=d%jj_begin-halo_size,d%jj_end+halo_size
-!           DO i=d%ii_begin-halo_size,d%ii_end+halo_size
-!               IF (d%own(i,j) .OR. single) THEN
-!               n=n+1
-!               CALL xyz2lonlat(d%xyz(:,i,j),lon(n),lat(n))
-!               lon(n)=lon(n)*180/Pi
-!               lat(n)=lat(n)*180/Pi
-!               DO k=0,5
-!                 CALL xyz2lonlat(d%vertex(:,k,i,j),bounds_lon(k,n), bounds_lat(k,n))
-!                 bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
-!                 bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
-!               ENDDO
-!             ENDIF
-!           ENDDO
-!         ENDDO
-!         status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ ncell /),count=(/ n /))
-!         status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ ncell /),count=(/ n /))
-!         status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,ncell /),count=(/ nvert,n /))
-!         status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,ncell /),count=(/ nvert,n /))
-! 
-!         ncell=ncell+n
-!         DEALLOCATE(lon,lat,bounds_lon,bounds_lat)
-!     END DO 
-
-!   ELSE IF (Field(ind_b)%field_type==field_Z) THEN
-!       nvert=3
-!       DO ind=ind_b,ind_e
-!         d=>domain(ind)
-!       
-!         DO j=d%jj_begin+1,d%jj_end
-!           DO i=d%ii_begin,d%ii_end-1
-!             ncell=ncell+1
-!           ENDDO
-!         ENDDO
-
-!         DO j=d%jj_begin,d%jj_end-1
-!           DO i=d%ii_begin+1,d%ii_end
-!             ncell=ncell+1
-!           ENDDO
-!         ENDDO
-
-!       END DO
-!     
-!       status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid)
-!       FieldId(NbField)=ncid
-!       status = NF90_DEF_DIM(ncid,'cell',ncell,ncellId)
-!       status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid)
-
-!       IF (Field(ind_b)%ndim==2)  THEN
-!         FieldVarId(NbField)%size=1
-!         ALLOCATE(FieldVarId(NbField)%nc_id(1))
-!       ELSE IF (Field(ind_b)%ndim==3)  THEN
-!         FieldVarId(NbField)%size=1
-!         ALLOCATE(FieldVarId(NbField)%nc_id(1))
-!         status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval3d,2),dim3id)
-!       ELSE IF (Field(1)%ndim==4) THEN
-!         FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3)
-!         ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size))
-!         status = NF90_DEF_DIM(ncid,'Z',size(field(ind_b)%rval4d,2),dim3id)
-!       ENDIF
-
-
-!     
-!       status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId)
-!     
-!       status = NF90_DEF_VAR(ncid,'lon',NF90_DOUBLE,(/ ncellId /),lonId)
-!       status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude")
-!       status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east")
-!       status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon")
-!       status = NF90_DEF_VAR(ncid,'lat',NF90_DOUBLE,(/ ncellId /),latId)
-!       status = NF90_PUT_ATT(ncid,latId,"long_name","latitude")
-!       status = NF90_PUT_ATT(ncid,latId,"units","degrees_north")
-!       status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat")
-!       status = NF90_DEF_VAR(ncid,'bounds_lon',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId)
-!       status = NF90_DEF_VAR(ncid,'bounds_lat',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId)
-
-
-!       IF (Field(ind_b)%ndim==2) THEN
-!         status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1))
-!         status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat")
-!       ELSE IF (Field(ind_b)%ndim==3) THEN
-!         status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),NF90_DOUBLE,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1))
-!         status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat")
-!       ELSE IF (Field(ind_b)%ndim==4) THEN
-!         DO q=1,FieldVarId(NbField)%size
-!           status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)//int2str(q)),NF90_DOUBLE,             &
-!                                 (/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(q))
-!           status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(q),"coordinates","lon lat")
-!         ENDDO        
-!       ENDIF 
-!       
-!       status = NF90_ENDDEF(ncid)      
-
-!       ncell=1
-!       DO ind=ind_b,ind_e
-!         d=>domain(ind)
-!         CALL swap_geometry(ind)
-!         CALL swap_dimensions(ind)
-! 
-!         n=0
-!         DO j=jj_begin+1,jj_end
-!           DO i=ii_begin,ii_end-1
-!             n=n+1
-!           ENDDO
-!         ENDDO
-
-!         DO j=jj_begin,jj_end-1
-!           DO i=ii_begin+1,ii_end
-!             n=n+1
-!           ENDDO
-!         ENDDO
-
-!        ALLOCATE(lon(n),lat(n),bounds_lon(0:nvert-1,n),bounds_lat(0:nvert-1,n))
-!         
-!         n=0  
-!      
-!         DO j=jj_begin+1,jj_end
-!           DO i=ii_begin,ii_end-1
-!             nij=(j-1)*iim+i
-!             n=n+1
-!             CALL xyz2lonlat(xyz_v(nij+z_down,:)/radius,lon(n),lat(n))
-!             lon(n)=lon(n)*180/Pi
-!!             IF (lon(n)<0) lon(n)=lon(n)+360
-!             lat(n)=lat(n)*180/Pi
-!             CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n))
-!             CALL xyz2lonlat(xyz_i(nij+t_ldown,:)/radius,bounds_lon(1,n), bounds_lat(1,n))
-!             CALL xyz2lonlat(xyz_i(nij+t_rdown,:)/radius,bounds_lon(2,n), bounds_lat(2,n))
-
-!             DO k=0,2
-!               bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
-!               bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
-!                IF (bounds_lon(k,n)<0) bounds_lon(k,n)=bounds_lon(k,n)+360
-!             ENDDO
-!           ENDDO
-!         ENDDO
-
-!         DO j=jj_begin,jj_end-1
-!           DO i=ii_begin+1,ii_end
-!             nij=(j-1)*iim+i
-!             n=n+1
-!             CALL xyz2lonlat(xyz_v(nij+z_up,:)/radius,lon(n),lat(n))
-!             lon(n)=lon(n)*180/Pi
-!              IF (lon(n)<0) lon(n)=lon(n)+360
-!             lat(n)=lat(n)*180/Pi
-!             CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n))
-!             CALL xyz2lonlat(xyz_i(nij+t_rup,:)/radius,bounds_lon(1,n), bounds_lat(1,n))
-!             CALL xyz2lonlat(xyz_i(nij+t_lup,:)/radius,bounds_lon(2,n), bounds_lat(2,n))
-
-!             DO k=0,2
-!               bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
-!               bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
-!                IF (bounds_lon(k,n)<0) bounds_lon(k,n)=bounds_lon(k,n)+360
-!             ENDDO
-!           ENDDO
-!         ENDDO
-!         
-!         
-!         status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ ncell /),count=(/ n /))
-!         status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ ncell /),count=(/ n /))
-!         status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,ncell /),count=(/ nvert,n /))
-!         status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,ncell /),count=(/ nvert,n /))
-! 
-!         ncell=ncell+n
-!         DEALLOCATE(lon,lat,bounds_lon,bounds_lat)
-!     END DO          
-!   ENDIF
-
-
-!     
-!      status = NF90_CLOSE(ncid)
-
-!   END SUBROUTINE Create_Header
-
-
-
+       END DO
+       DEALLOCATE(field_val3d)
+     END SUBROUTINE write_4d
+
+   END SUBROUTINE Writefield_gen
    
-    SUBROUTINE Create_header_gen(name_in,field,domain_type,ind_b_in,ind_e_in,once)
-    USE netcdf_mod
-    USE field_mod
-    USE domain_mod
-    USE metric
-    USE spherical_geom_mod
-    IMPLICIT NONE
-      CHARACTER(LEN=*),INTENT(IN) :: name_in
-      TYPE(t_field),POINTER :: field(:)
-      TYPE(t_domain),INTENT(IN),TARGET :: domain_type(:)
-      INTEGER,OPTIONAL,INTENT(IN) :: ind_b_in
-      INTEGER,OPTIONAL,INTENT(IN) :: ind_e_in
-      LOGICAL,INTENT(IN) :: once
-      
-      INTEGER :: ncell
-      INTEGER :: nvert
-      REAL(rstd),ALLOCATABLE :: lon(:),lat(:),bounds_lon(:,:),bounds_lat(:,:)
-      TYPE(t_domain),POINTER :: d
-      INTEGER :: nvertId,ncid,lonId,latId,bounds_lonId,bounds_latId,timeId,ncellId
-      INTEGER :: dim3id,dim4id
-      INTEGER :: status
-      INTEGER :: ind,i,j,k,n,q
-      INTEGER :: iie,jje,iin,jjn
-      INTEGER :: ind_b,ind_e
-      INTEGER :: halo_size
-      LOGICAL :: single 
-      INTEGER :: nij
-      CHARACTER(LEN=255) :: name
-      INTEGER :: l,level_size, levId, dimlevId
-            
-      name=TRIM(ADJUSTL(name_in))
-
-      IF (PRESENT(ind_b_in) .AND. .NOT. PRESENT(ind_e_in)) THEN
+   SUBROUTINE Create_header_gen(name_in,field,domain_type,ind_b_in,ind_e_in,once)
+     USE netcdf_mod
+     USE field_mod
+     USE domain_mod
+     USE metric
+     USE spherical_geom_mod
+     CHARACTER(LEN=*),INTENT(IN) :: name_in
+     TYPE(t_field),POINTER :: field(:)
+     TYPE(t_domain),INTENT(IN),TARGET :: domain_type(:)
+     INTEGER,OPTIONAL,INTENT(IN) :: ind_b_in
+     INTEGER,OPTIONAL,INTENT(IN) :: ind_e_in
+     LOGICAL,INTENT(IN) :: once
+     
+     TYPE(t_cellset), POINTER :: cells(:)
+     INTEGER :: ncell
+     INTEGER :: nvert
+     REAL(rstd),ALLOCATABLE :: lon(:),lat(:),bounds_lon(:,:),bounds_lat(:,:)
+     TYPE(t_domain),POINTER :: d
+     INTEGER :: nvertId,ncid,lonId,latId,bounds_lonId,bounds_latId,timeId,ncellId
+     INTEGER :: dim3id,dim4id
+     INTEGER :: status
+     INTEGER :: ind,i,j,k,n,q
+     INTEGER :: iie,jje,iin,jjn
+     INTEGER :: ind_b,ind_e, n_begin, n_end
+     LOGICAL :: single 
+     INTEGER :: nij
+     CHARACTER(LEN=255) :: name
+     INTEGER :: l,level_size, levId, dimlevId
+     
+     name=TRIM(ADJUSTL(name_in))
+     
+     IF (PRESENT(ind_b_in) .AND. .NOT. PRESENT(ind_e_in)) THEN
         name=TRIM(name)//"_"//TRIM(int2str(ind_b))
         ind_b=ind_b_in
         ind_e=ind_b_in
-        halo_size=1
         single=.TRUE.
-      ELSE IF (.NOT. PRESENT(ind_b_in) .AND. PRESENT(ind_e_in)) THEN
+     ELSE IF (.NOT. PRESENT(ind_b_in) .AND. PRESENT(ind_e_in)) THEN
         name=TRIM(name)//"_"//TRIM(int2str(ind_e))
         ind_b=ind_e_in
         ind_e=ind_e_in
-        halo_size=1
         single=.TRUE.
-      ELSE IF ( PRESENT(ind_b_in) .AND. PRESENT(ind_e_in)) THEN
+     ELSE IF ( PRESENT(ind_b_in) .AND. PRESENT(ind_e_in)) THEN
         ind_b=ind_b_in
         ind_e=ind_e_in
-        halo_size=0
         single=.FALSE.
-      ELSE
-        halo_size=0
+     ELSE
         ind_b=1
         ind_e=ndomain
         single=.FALSE.
-      ENDIF      
-                
-      NbField=NbField+1
-      FieldName(NbField)=TRIM(ADJUSTL(name))
-      FieldIndex(NbField)=1
-      
-      ncell=0
-      
-      IF (Field(ind_b)%field_type==field_T) THEN
-        nvert=6
-        
-        DO ind=ind_b,ind_e
-          d=>domain_type(ind)
-        
-          DO j=d%jj_begin-halo_size,d%jj_end+halo_size
-            DO i=d%ii_begin-halo_size,d%ii_end+halo_size
-              IF (single .OR. d%assign_domain(i,j)==ind) ncell=ncell+1
-            ENDDO
-          ENDDO
-
-        END DO
-      
-        status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid)
-        FieldId(NbField)=ncid
+     ENDIF
+     
+     NbField=NbField+1
+     FieldName(NbField)=TRIM(ADJUSTL(name))
+     FieldIndex(NbField)=1
+     
+     PRINT *, 'Creating header for writefield : ', TRIM(FieldName(NbField)), ndomain, ind_b, ind_e ! FIXME
+     PRINT *, mesh_glo%ndomain, SIZE(mesh_glo%primal_own)
+     SELECT CASE(Field(ind_b)%field_type)
+     CASE(field_T)
+        PRINT *, '    field_type == field_T' ! FIXME
+        IF(single) THEN ! include halos
+           cells => mesh_glo%primal_all
+        ELSE
+           cells => mesh_glo%primal_own
+        END IF
+     CASE(field_Z)
+        PRINT *, '    field_type == field_Z' ! FIXME
+        IF(single) THEN ! include halos
+           cells => mesh_glo%dual_all
+        ELSE
+           cells => mesh_glo%dual_own
+        END IF
+     END SELECT
+     
+     PRINT *, 'writefield : ind_b, ind_e :', ind_b, ind_e
+
+     ncell=0
+     DO ind=ind_b,ind_e
+        nvert=SIZE(cells(ind)%bnds_lon,1)
+        ncell = ncell + cells(ind)%ncell
+     END DO
+     PRINT *, 'writefield : found ',ncell, ' cells.' ! FIXME
+     
+     status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid)
+     FieldId(NbField)=ncid
+     status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId)
+     
+     level_size=0
+     SELECT CASE(Field(ind_b)%ndim)
+     CASE(2)
+        FieldVarId(NbField)%size=1
+        ALLOCATE(FieldVarId(NbField)%nc_id(1))
+     CASE(3)
+        FieldVarId(NbField)%size=1
+        ALLOCATE(FieldVarId(NbField)%nc_id(1))
+        status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval3d,2),dim3id)
+        level_size=size(field(ind_b)%rval3d,2)
+     CASE(4)
+        FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3)
+        ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size))
+        status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval4d,2),dim3id)
+        level_size=size(field(ind_b)%rval4d,2)
+        !          status = NF90_DEF_DIM(ncid,'Q',size(field(ind_b)%rval4d,3),dim4id)
+     END SELECT
+     
+     PRINT*,"create_header_gen : LEVEL_SIZE=",level_size
+     IF (level_size>0) THEN
+        status = NF90_DEF_VAR(ncid,'lev',NF90_DOUBLE,(/ dim3id /),levId)
+        status = NF90_PUT_ATT(ncid,levId,"axis","Z")
+     ENDIF
+     
+     SELECT CASE(Field(ind_b)%field_type)
+     CASE(field_T)
         status = NF90_DEF_DIM(ncid,'cell_i',ncell,ncellId)
         status = NF90_DEF_DIM(ncid,'nvert_i',nvert,nvertid)
-        level_size=0
-        IF (Field(ind_b)%ndim==2)  THEN
-          FieldVarId(NbField)%size=1
-          ALLOCATE(FieldVarId(NbField)%nc_id(1))
-        ELSE IF (Field(ind_b)%ndim==3)  THEN
-          FieldVarId(NbField)%size=1
-          ALLOCATE(FieldVarId(NbField)%nc_id(1))
-          status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval3d,2),dim3id)
-          level_size=size(field(ind_b)%rval3d,2)
-        ELSE IF (Field(1)%ndim==4) THEN
-          FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3)
-          ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size))
-          status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval4d,2),dim3id)
-          level_size=size(field(ind_b)%rval4d,2)
-!          status = NF90_DEF_DIM(ncid,'Q',size(field(ind_b)%rval4d,3),dim4id)
-        ENDIF
-        PRINT*,"LEVEL_SIZE=",level_size
-
-        status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId)
-        IF (level_size>0) THEN
-          status = NF90_DEF_VAR(ncid,'lev',NF90_DOUBLE,(/ dim3id /),levId)
-          status = NF90_PUT_ATT(ncid,levId,"axis","Z")
-        ENDIF
-        
         status = NF90_DEF_VAR(ncid,'lon_i',NF90_DOUBLE,(/ ncellId /),lonId)
-        status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude")
-        status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east")
+        status = NF90_DEF_VAR(ncid,'lat_i',NF90_DOUBLE,(/ ncellId /),latId)
         status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon_i")
-        status = NF90_DEF_VAR(ncid,'lat_i',NF90_DOUBLE,(/ ncellId /),latId)
-        status = NF90_PUT_ATT(ncid,latId,"long_name","latitude")
-        status = NF90_PUT_ATT(ncid,latId,"units","degrees_north")
         status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat_i")
         status = NF90_DEF_VAR(ncid,'bounds_lon_i',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId)
         status = NF90_DEF_VAR(ncid,'bounds_lat_i',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId)
-
-        IF (Field(ind_b)%ndim==2) THEN
-          IF (once) THEN 
-            status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId /),FieldVarId(NbField)%nc_id(1))
-          ELSE
-             status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1))
-          ENDIF
-          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon_i lat_i")
-        ELSE IF (Field(ind_b)%ndim==3) THEN
-          IF (once) THEN 
-            status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,dim3id /),FieldVarId(NbField)%nc_id(1))
-          ELSE
-            status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1))
-          ENDIF
-          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon_i lat_i")
-        ELSE IF (Field(ind_b)%ndim==4) THEN
-          DO i=1,FieldVarId(NbField)%size
-            IF (once) THEN 
-              status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name))//int2str(i),ncprec,(/ ncellId,dim3id /),  &
-                                    FieldVarId(NbField)%nc_id(i))
-            ELSE
-              status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name))//int2str(i),ncprec,(/ ncellId,dim3id,timeId /),  &
-                                   FieldVarId(NbField)%nc_id(i))
-            ENDIF
-            status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(i),"coordinates","lon_i lat_i")
-          ENDDO        
-        ENDIF 
-  
-      
-        status = NF90_ENDDEF(ncid)      
-
-        if (level_size>0) status = NF90_PUT_VAR(ncid,levId,(/ (l,l=1,level_size) /))
-         
-         ALLOCATE(lon(ncell),lat(ncell),bounds_lon(0:nvert-1,ncell),bounds_lat(0:nvert-1,ncell))
-          
-         n=0  
-         DO ind=ind_b,ind_e
-           d=>domain_type(ind)
-           DO j=d%jj_begin-halo_size,d%jj_end+halo_size
-             DO i=d%ii_begin-halo_size,d%ii_end+halo_size
-               IF (d%assign_domain(i,j)==ind .OR. single) THEN
-                 n=n+1
-                 CALL xyz2lonlat(d%xyz(:,i,j),lon(n),lat(n))
-                 lon(n)=lon(n)*180/Pi
-                 lat(n)=lat(n)*180/Pi
-                 DO k=0,5
-                   CALL xyz2lonlat(d%vertex(:,k,i,j),bounds_lon(k,n), bounds_lat(k,n))
-                   bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
-                   bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
-                 ENDDO
-               ENDIF
-             ENDDO
-           ENDDO
-         ENDDO
-         status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ 1 /),count=(/ ncell /))
-         status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ 1 /),count=(/ ncell /))
-         status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,1 /),count=(/ nvert,ncell /))
-         status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,1 /),count=(/ nvert,ncell /))
-  
-         DEALLOCATE(lon,lat,bounds_lon,bounds_lat)
-
-    ELSE IF (Field(ind_b)%field_type==field_Z) THEN
-        nvert=3
-        DO ind=ind_b,ind_e
-          d=>domain_type(ind)
-        
-          DO j=d%jj_begin+1,d%jj_end
-            DO i=d%ii_begin,d%ii_end-1
-              ncell=ncell+1
-            ENDDO
-          ENDDO
-
-          DO j=d%jj_begin,d%jj_end-1
-            DO i=d%ii_begin+1,d%ii_end
-              ncell=ncell+1
-            ENDDO
-          ENDDO
-
-        END DO
-      
-        status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid)
-        FieldId(NbField)=ncid
+     CASE(field_Z)
         status = NF90_DEF_DIM(ncid,'cell_v',ncell,ncellId)
         status = NF90_DEF_DIM(ncid,'nvert_v',nvert,nvertid)
-
-        IF (Field(ind_b)%ndim==2)  THEN
-          FieldVarId(NbField)%size=1
-          ALLOCATE(FieldVarId(NbField)%nc_id(1))
-        ELSE IF (Field(ind_b)%ndim==3)  THEN
-          FieldVarId(NbField)%size=1
-          ALLOCATE(FieldVarId(NbField)%nc_id(1))
-          status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval3d,2),dim3id)
-        ELSE IF (Field(1)%ndim==4) THEN
-          FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3)
-          ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size))
-          status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval4d,2),dim3id)
-        ENDIF
-
-
-      
-        status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId)
-      
         status = NF90_DEF_VAR(ncid,'lon_v',NF90_DOUBLE,(/ ncellId /),lonId)
-        status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude")
-        status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east")
+        status = NF90_DEF_VAR(ncid,'lat_v',NF90_DOUBLE,(/ ncellId /),latId)
         status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon_v")
-        status = NF90_DEF_VAR(ncid,'lat_v',NF90_DOUBLE,(/ ncellId /),latId)
-        status = NF90_PUT_ATT(ncid,latId,"long_name","latitude")
-        status = NF90_PUT_ATT(ncid,latId,"units","degrees_north")
         status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat_v")
         status = NF90_DEF_VAR(ncid,'bounds_lon_v',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId)
         status = NF90_DEF_VAR(ncid,'bounds_lat_v',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId)
-
-
-        IF (Field(ind_b)%ndim==2) THEN
-          IF (once) THEN 
-            status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId /),FieldVarId(NbField)%nc_id(1))
-          ELSE
-            status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1))
-          ENDIF
-          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon_v lat_v")
-        ELSE IF (Field(ind_b)%ndim==3) THEN
-          IF (once) THEN 
-            status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,dim3id /),FieldVarId(NbField)%nc_id(1))
-          ELSE
-            status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1))
-          ENDIF
-          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon_v lat_v")
-        ELSE IF (Field(ind_b)%ndim==4) THEN
-          DO q=1,FieldVarId(NbField)%size
-            IF (once) THEN 
-              status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)//int2str(q)),ncprec,             &
-                                    (/ ncellId,dim3id /),FieldVarId(NbField)%nc_id(q))
-            ELSE
-              status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)//int2str(q)),ncprec,             &
-                                    (/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(q))
-            ENDIF
-            status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(q),"coordinates","lon_v lat_v")
-          ENDDO        
-        ENDIF 
-        
-        status = NF90_ENDDEF(ncid)      
-
-         ALLOCATE(lon(ncell),lat(ncell),bounds_lon(0:nvert-1,ncell),bounds_lat(0:nvert-1,ncell))
-          
-         n=0  
-       
-         DO ind=ind_b,ind_e
-           d=>domain_type(ind)
-           DO j=d%jj_begin+1,d%jj_end
-             DO i=d%ii_begin,d%ii_end-1
-               nij=(j-1)*d%iim+i
-               n=n+1
-               CALL xyz2lonlat(d%vertex(:,vdown,i,j),lon(n),lat(n))
-               lon(n)=lon(n)*180/Pi
-               lat(n)=lat(n)*180/Pi
-
-               CALL xyz2lonlat(d%xyz(:,i,j),bounds_lon(0,n), bounds_lat(0,n))
-               CALL xyz2lonlat(d%xyz(:,i,j-1),bounds_lon(1,n), bounds_lat(1,n))
-               CALL xyz2lonlat(d%xyz(:,i+1,j-1),bounds_lon(2,n), bounds_lat(2,n))
-
-               DO k=0,2
-                 bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
-                 bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
-               ENDDO
-             ENDDO
-           ENDDO
-
-           DO j=d%jj_begin,d%jj_end-1
-             DO i=d%ii_begin+1,d%ii_end
-               nij=(j-1)*d%iim+i
-               n=n+1
-               CALL xyz2lonlat(d%vertex(:,vup,i,j),lon(n),lat(n))
-               lon(n)=lon(n)*180/Pi
-               lat(n)=lat(n)*180/Pi
-               CALL xyz2lonlat(d%xyz(:,i,j),bounds_lon(0,n), bounds_lat(0,n))
-               CALL xyz2lonlat(d%xyz(:,i,j+1),bounds_lon(1,n), bounds_lat(1,n))
-               CALL xyz2lonlat(d%xyz(:,i-1,j+1),bounds_lon(2,n), bounds_lat(2,n))
-
-               DO k=0,2
-                 bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
-                 bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
-               ENDDO
-             ENDDO
-           ENDDO
-         ENDDO 
-          
-         status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ 1 /),count=(/ ncell /))
-         status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ 1 /),count=(/ ncell /))
-         status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,1 /),count=(/ nvert,ncell /))
-         status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,1 /),count=(/ nvert,ncell /))
-  
-          DEALLOCATE(lon,lat,bounds_lon,bounds_lat)
-      ENDIF
-
-
-    END SUBROUTINE Create_Header_gen
-
-    SUBROUTINE Create_header_mpi(name,field,nind)
-    USE netcdf_mod
-    USE field_mod
-    USE domain_mod
-    USE spherical_geom_mod
-    USE dimensions
-    USE geometry
-    USE mpi_mod
-    USE mpipara
-    IMPLICIT NONE
-      CHARACTER(LEN=*) :: name
-      CHARACTER(LEN=LEN_TRIM(ADJUSTL(name))) :: name_adj
-      TYPE(t_field),POINTER :: field(:)
-      INTEGER,OPTIONAL,INTENT(IN) :: nind
-      INTEGER :: ncell
-      INTEGER :: nvert
-      REAL(rstd),ALLOCATABLE :: lon(:),lat(:),bounds_lon(:,:),bounds_lat(:,:)
-      TYPE(t_domain),POINTER :: d
-      INTEGER :: nvertId,ncid,lonId,latId,bounds_lonId,bounds_latId,timeId,ncellId
-      INTEGER :: dim3id,dim4id
-      INTEGER :: status
-      INTEGER :: ind,i,j,k,n,q
-      INTEGER :: iie,jje,iin,jjn
-      INTEGER :: ind_b,ind_e
-      INTEGER :: halo_size
-      LOGICAL :: single 
-      INTEGER :: nij
-      INTEGER :: ncell_glo(0:mpi_size-1)
-      INTEGER :: displ, ncell_tot
-      
-          
-      NbField=NbField+1
-      name_adj=TRIM(ADJUSTL(name))  ! work around ICE with pgf90
-      FieldName(NbField)=name_adj
-      FieldIndex(NbField)=1
-      
-      IF (PRESENT(nind)) THEN
-        ind_b=nind
-        ind_e=nind
-        halo_size=1
-        single=.TRUE.
-      ELSE
-        ind_b=1
-        ind_e=ndomain
-        halo_size=0
-        single=.FALSE.
-      ENDIF
-      
-      ncell=0
-      
-      IF (Field(ind_b)%field_type==field_T) THEN
-        nvert=6
-        
-        DO ind=ind_b,ind_e
-          d=>domain(ind)
-        
-          DO j=d%jj_begin-halo_size,d%jj_end+halo_size
-            DO i=d%ii_begin-halo_size,d%ii_end+halo_size
-              IF (single .OR. domain(ind)%own(i,j)) ncell=ncell+1
-            ENDDO
-          ENDDO
-
+     END SELECT
+     
+     status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude")
+     status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east")
+     status = NF90_PUT_ATT(ncid,latId,"long_name","latitude")
+     status = NF90_PUT_ATT(ncid,latId,"units","degrees_north")
+     
+     SELECT CASE(Field(ind_b)%ndim)
+     CASE(2)
+        IF (once) THEN 
+           status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec, &
+                (/ ncellId /),FieldVarId(NbField)%nc_id(1))
+        ELSE
+           status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec, &
+                (/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1))
+        END IF
+        CALL put_att_coordinates(1)
+     CASE(3)
+        IF (once) THEN 
+           status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec, &
+                (/ ncellId,dim3id /),FieldVarId(NbField)%nc_id(1))
+        ELSE
+           status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec, &
+                (/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1))
+        END IF
+        CALL put_att_coordinates(1)
+     CASE(4)
+        DO i=1,FieldVarId(NbField)%size
+           IF (once) THEN 
+              status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name))//int2str(i),ncprec,(/ ncellId,dim3id /),  &
+                   FieldVarId(NbField)%nc_id(i))
+           ELSE
+              status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name))//int2str(i),ncprec,(/ ncellId,dim3id,timeId /),  &
+                   FieldVarId(NbField)%nc_id(i))
+           END IF
+           CALL put_att_coordinates(i)
+           status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(i),"coordinates","lon_i lat_i")
         END DO
-      
-        CALL MPI_ALLGATHER(ncell,1,MPI_INTEGER,ncell_glo,1,MPI_INTEGER,comm_icosa,ierr)
-
-        displ=0
-        DO i=1,mpi_rank
-          displ=displ+ncell_glo(i-1)
-        ENDDO
-        FieldVarId(NbField)%displ=displ
-        ncell_tot=sum(ncell_glo(:))
-        
-!        status = NF90_CREATE_PAR(TRIM(ADJUSTL(name))//'.nc', IOR(NF90_NETCDF4, NF90_MPIIO), comm_icosa, MPI_INFO_NULL, ncid)
-        FieldId(NbField)=ncid
-        status = NF90_DEF_DIM(ncid,'cell',ncell_tot,ncellId)
-        status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid)
-
-        IF (Field(ind_b)%ndim==2)  THEN
-          FieldVarId(NbField)%size=1
-          ALLOCATE(FieldVarId(NbField)%nc_id(1))
-        ELSE IF (Field(ind_b)%ndim==3)  THEN
-          FieldVarId(NbField)%size=1
-          ALLOCATE(FieldVarId(NbField)%nc_id(1))
-          status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval3d,2),dim3id)
-        ELSE IF (Field(1)%ndim==4) THEN
-          FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3)
-          ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size))
-          status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval4d,2),dim3id)
-!          status = NF90_DEF_DIM(ncid,'Q',size(field(ind_b)%rval4d,3),dim4id)
-        ENDIF
-      
-        status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId)
-      
-        status = NF90_DEF_VAR(ncid,'lon',NF90_DOUBLE,(/ ncellId /),lonId)
-        status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude")
-        status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east")
-        status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon")
-        status = NF90_DEF_VAR(ncid,'lat',NF90_DOUBLE,(/ ncellId /),latId)
-        status = NF90_PUT_ATT(ncid,latId,"long_name","latitude")
-        status = NF90_PUT_ATT(ncid,latId,"units","degrees_north")
-        status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat")
-        status = NF90_DEF_VAR(ncid,'bounds_lon',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId)
-        status = NF90_DEF_VAR(ncid,'bounds_lat',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId)
-
-        IF (Field(ind_b)%ndim==2) THEN
-          status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1))
-          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat")
-          status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(1), NF90_CHUNKED, (/ncell_tot,1/))
-        ELSE IF (Field(ind_b)%ndim==3) THEN
-          status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1))
-          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat")
-          status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(1), NF90_CHUNKED,   &
-                                         (/ncell_tot,size(field(ind_b)%rval3d,2),1/))
-        ELSE IF (Field(ind_b)%ndim==4) THEN
-          DO i=1,FieldVarId(NbField)%size
-            status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name))//int2str(i),ncprec,(/ ncellId,dim3id,timeId /),  &
-                                  FieldVarId(NbField)%nc_id(i))
-            status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(i),"coordinates","lon lat")
-            status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(q), NF90_CHUNKED,   &
-                                           (/ncell_tot,size(field(ind_b)%rval4d,2),1/))
-          ENDDO        
-        ENDIF 
-  
-      
-        status = NF90_ENDDEF(ncid)      
-
-        ncell=1
-        DO ind=ind_b,ind_e
-          d=>domain(ind)
-  
-          n=0
-          DO j=d%jj_begin-halo_size,d%jj_end+halo_size
-            DO i=d%ii_begin-halo_size,d%ii_end+halo_size
-              IF (single .OR. d%own(i,j)) n=n+1
-            ENDDO
-          ENDDO
-          
-         ALLOCATE(lon(n),lat(n),bounds_lon(0:nvert-1,n),bounds_lat(0:nvert-1,n))
-          
-          n=0  
-          DO j=d%jj_begin-halo_size,d%jj_end+halo_size
-            DO i=d%ii_begin-halo_size,d%ii_end+halo_size
-                IF (d%own(i,j) .OR. single) THEN
-                n=n+1
-                CALL xyz2lonlat(d%xyz(:,i,j),lon(n),lat(n))
-                lon(n)=lon(n)*180/Pi
-                lat(n)=lat(n)*180/Pi
-                DO k=0,5
-                  CALL xyz2lonlat(d%vertex(:,k,i,j),bounds_lon(k,n), bounds_lat(k,n))
-                  bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
-                  bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
-                ENDDO
-              ENDIF
-            ENDDO
-          ENDDO
-          status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ displ+ncell /),count=(/ n /))
-          status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ displ+ncell /),count=(/ n /))
-          status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,displ+ncell /),count=(/ nvert,n /))
-          status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,displ+ncell /),count=(/ nvert,n /))
-  
-          ncell=ncell+n
-          DEALLOCATE(lon,lat,bounds_lon,bounds_lat)
-      END DO 
-
-    ELSE IF (Field(ind_b)%field_type==field_Z) THEN
-        nvert=3
-        DO ind=ind_b,ind_e
-          d=>domain(ind)
-        
-          DO j=d%jj_begin+1,d%jj_end
-            DO i=d%ii_begin,d%ii_end-1
-              ncell=ncell+1
-            ENDDO
-          ENDDO
-
-          DO j=d%jj_begin,d%jj_end-1
-            DO i=d%ii_begin+1,d%ii_end
-              ncell=ncell+1
-            ENDDO
-          ENDDO
-
-        END DO
-        
-        CALL MPI_ALLGATHER(ncell,1,MPI_INTEGER,ncell_glo,1,MPI_INTEGER,comm_icosa,ierr)
-
-        displ=0
-        DO i=1,mpi_rank
-          displ=displ+ncell_glo(i-1)
-        ENDDO
-        FieldVarId(NbField)%displ=displ
-        ncell_tot=sum(ncell_glo(:))
-              
-!        status = NF90_CREATE_PAR(TRIM(ADJUSTL(name))//'.nc',IOR(NF90_NETCDF4, NF90_MPIIO), comm_icosa, MPI_INFO_NULL, ncid)
-!        status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid)
-        FieldId(NbField)=ncid
-        status = NF90_DEF_DIM(ncid,'cell',ncell_tot,ncellId)
-        status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid)
-
-        IF (Field(ind_b)%ndim==2)  THEN
-          FieldVarId(NbField)%size=1
-          ALLOCATE(FieldVarId(NbField)%nc_id(1))
-        ELSE IF (Field(ind_b)%ndim==3)  THEN
-          FieldVarId(NbField)%size=1
-          ALLOCATE(FieldVarId(NbField)%nc_id(1))
-          status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval3d,2),dim3id)
-        ELSE IF (Field(1)%ndim==4) THEN
-          FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3)
-          ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size))
-          status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval4d,2),dim3id)
-        ENDIF
-
-
-      
-        status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId)
-      
-        status = NF90_DEF_VAR(ncid,'lon',NF90_DOUBLE,(/ ncellId /),lonId)
-        status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude")
-        status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east")
-        status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon")
-        status = NF90_DEF_VAR(ncid,'lat',NF90_DOUBLE,(/ ncellId /),latId)
-        status = NF90_PUT_ATT(ncid,latId,"long_name","latitude")
-        status = NF90_PUT_ATT(ncid,latId,"units","degrees_north")
-        status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat")
-        status = NF90_DEF_VAR(ncid,'bounds_lon',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId)
-        status = NF90_DEF_VAR(ncid,'bounds_lat',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId)
-
-
-        IF (Field(ind_b)%ndim==2) THEN
-          status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1))
-          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat")
-          status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(1), NF90_CHUNKED, (/ncell_tot,1/))
-        ELSE IF (Field(ind_b)%ndim==3) THEN
-          status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1))
-          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat")
-          status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(1), NF90_CHUNKED,   &
-                                         (/ncell_tot,size(field(ind_b)%rval3d,2),1/))
-        ELSE IF (Field(ind_b)%ndim==4) THEN
-          DO q=1,FieldVarId(NbField)%size
-            status = NF90_DEF_VAR(ncid,name_adj//int2str(q),ncprec,             &
-                                  (/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(q))
-            status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(q),"coordinates","lon lat")
-           status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(q), NF90_CHUNKED, &
-                                          (/ncell_tot,size(field(ind_b)%rval4d,2),1/))
-          ENDDO        
-        ENDIF 
-        
-        status = NF90_ENDDEF(ncid)      
-
-        ncell=1
-        DO ind=ind_b,ind_e
-          d=>domain(ind)
-          CALL swap_geometry(ind)
-          CALL swap_dimensions(ind)
-  
-          n=0
-          DO j=jj_begin+1,jj_end
-            DO i=ii_begin,ii_end-1
-              n=n+1
-            ENDDO
-          ENDDO
-
-          DO j=jj_begin,jj_end-1
-            DO i=ii_begin+1,ii_end
-              n=n+1
-            ENDDO
-          ENDDO
-
-         ALLOCATE(lon(n),lat(n),bounds_lon(0:nvert-1,n),bounds_lat(0:nvert-1,n))
-          
-          n=0  
-       
-          DO j=jj_begin+1,jj_end
-            DO i=ii_begin,ii_end-1
-              nij=(j-1)*iim+i
-              n=n+1
-              CALL xyz2lonlat(xyz_v(nij+z_down,:)/radius,lon(n),lat(n))
-              lon(n)=lon(n)*180/Pi
-              lat(n)=lat(n)*180/Pi
-              CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n))
-              CALL xyz2lonlat(xyz_i(nij+t_ldown,:)/radius,bounds_lon(1,n), bounds_lat(1,n))
-              CALL xyz2lonlat(xyz_i(nij+t_rdown,:)/radius,bounds_lon(2,n), bounds_lat(2,n))
-
-              DO k=0,2
-                bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
-                bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
-              ENDDO
-            ENDDO
-          ENDDO
-
-          DO j=jj_begin,jj_end-1
-            DO i=ii_begin+1,ii_end
-              nij=(j-1)*iim+i
-              n=n+1
-              CALL xyz2lonlat(xyz_v(nij+z_up,:)/radius,lon(n),lat(n))
-              lon(n)=lon(n)*180/Pi
-              lat(n)=lat(n)*180/Pi
-              CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n))
-              CALL xyz2lonlat(xyz_i(nij+t_rup,:)/radius,bounds_lon(1,n), bounds_lat(1,n))
-              CALL xyz2lonlat(xyz_i(nij+t_lup,:)/radius,bounds_lon(2,n), bounds_lat(2,n))
-
-              DO k=0,2
-                bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
-                bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
-              ENDDO
-            ENDDO
-          ENDDO
-          
-          
-          status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ displ+ncell /),count=(/ n /))
-          status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ displ+ncell /),count=(/ n /))
-          status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,displ+ncell /),count=(/ nvert,n /))
-          status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,displ+ncell /),count=(/ nvert,n /))
-  
-          ncell=ncell+n
-          DEALLOCATE(lon,lat,bounds_lon,bounds_lat)
-      END DO          
-    ENDIF
-
-
-   END SUBROUTINE Create_Header_mpi  
-   
+     END SELECT
+     
+     status = NF90_ENDDEF(ncid)      
+     
+     if (level_size>0) status = NF90_PUT_VAR(ncid,levId,(/ (l,l=1,level_size) /))
+     
+     ALLOCATE(lon(ncell),lat(ncell))
+     ALLOCATE(bounds_lon(nvert,ncell), bounds_lat(nvert,ncell))
+     n_begin=0  
+     DO ind=ind_b,ind_e
+        n_end = n_begin + cells(ind)%ncell
+        PRINT *, 'create_header_gen ', n_begin, n_end, SHAPE(cells(ind)%bnds_lon), SHAPE(bounds_lon)
+        lat(n_begin+1:n_end) = cells(ind)%lat(:) *180./Pi
+        lon(n_begin+1:n_end) = cells(ind)%lon(:) *180./Pi
+        bounds_lon(:,n_begin+1:n_end) = cells(ind)%bnds_lon(:,:) *180./Pi
+        bounds_lat(:,n_begin+1:n_end) = cells(ind)%bnds_lat(:,:) *180./Pi
+        n_begin = n_end
+     END DO
+     status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ 1 /),count=(/ ncell /))
+     status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ 1 /),count=(/ ncell /))
+     status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,1 /),count=(/ nvert,ncell /))
+     status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,1 /),count=(/ nvert,ncell /))
+     
+     DEALLOCATE(lon,lat,bounds_lon,bounds_lat)
+     
+   CONTAINS
+     
+     SUBROUTINE put_att_coordinates(ind)
+       INTEGER :: ind
+       SELECT CASE(Field(ind_b)%field_type)
+       CASE(field_T)
+          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(ind), &
+               "coordinates","lon_i lat_i")
+       CASE(field_Z)
+          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(ind), &
+               "coordinates","lon_v lat_v")      
+       END SELECT
+     END SUBROUTINE put_att_coordinates
+     
+   END SUBROUTINE Create_Header_gen
+    
    SUBROUTINE Close_files
    USE netcdf_mod
-   IMPLICIT NONE
      INTEGER :: i,k,status
 !$OMP MASTER     
@@ -1794,23 +458,3 @@
      
      
-  function int2str(int)
-    implicit none
-    integer, parameter :: MaxLen=10
-    integer,intent(in) :: int
-    character(len=MaxLen) :: int2str
-    logical :: flag
-    integer :: i
-    flag=.true.
-    
-    i=int
-    
-    int2str=''
-    do while (flag)
-      int2str=CHAR(MOD(i,10)+48)//int2str
-      i=i/10
-      if (i==0) flag=.false.
-    enddo
-  end function int2str
-
 end module write_field_mod
-  
Index: codes/icosagcm/devel/src/output/write_field_mpi.f90
===================================================================
--- codes/icosagcm/devel/src/output/write_field_mpi.f90	(revision 880)
+++ codes/icosagcm/devel/src/output/write_field_mpi.f90	(revision 880)
@@ -0,0 +1,579 @@
+module write_field_mpi_mod
+  USE genmod
+  USE write_field_vars_mod
+  IMPLICIT NONE
+  PRIVATE  
+  
+!  TYPE ncvar
+!    INTEGER :: size
+!    INTEGER,POINTER :: nc_id(:)
+!    INTEGER :: displ
+!  END TYPE ncvar
+!
+!  INTEGER, PARAMETER :: MaxWriteField = 1000
+!  INTEGER, DIMENSION(MaxWriteField),SAVE :: FieldId
+!  TYPE(ncvar), dimension(MaxWriteField),SAVE :: FieldVarId
+!  INTEGER, DIMENSION(MaxWriteField),SAVE :: FieldIndex
+!  CHARACTER(len=255), DIMENSION(MaxWriteField) ::  FieldName 
+!   
+!  INTEGER,SAVE :: NbField = 0
+!  
+!  PUBLIC init_writeField, writefield, close_files
+  
+  CONTAINS
+  
+    SUBROUTINE Writefield_mpi(name_in,field,nind)
+    USE netcdf_mod
+    USE domain_mod
+    use field_mod
+    USE dimensions
+    USE geometry
+    IMPLICIT NONE
+      CHARACTER(LEN=*),INTENT(IN) :: name_in
+      TYPE(t_field),POINTER :: field(:)
+      INTEGER,OPTIONAL,INTENT(IN) :: nind
+      REAL(r8),ALLOCATABLE :: field_val2d(:)
+      REAL(r8),ALLOCATABLE :: field_val3d(:,:)
+      REAL(r8),ALLOCATABLE :: field_val4d(:,:,:)
+      TYPE(t_domain),POINTER :: d
+      INTEGER :: Index
+      INTEGER :: ind,i,j,l,k,n,ncell,q
+      INTEGER :: iie,jje,iin,jjn
+      INTEGER :: status
+      CHARACTER(len=255) :: name
+      CHARACTER(len=255) :: str_ind
+      INTEGER :: ind_b,ind_e
+      INTEGER :: halo_size
+      LOGICAL :: single
+      INTEGER :: displ
+      
+      name=TRIM(ADJUSTL(name_in))
+
+      IF (PRESENT(nind)) THEN
+        name=TRIM(name)//"_"//TRIM(int2str(nind))
+        PRINT *,"NAME",nind,int2str(nind),name
+        ind_b=nind
+        ind_e=nind
+        halo_size=1
+        single=.TRUE.
+      ELSE
+        ind_b=1
+        ind_e=ndomain
+        halo_size=0
+        single=.FALSE.
+      ENDIF      
+
+      Index=GetFieldIndex(name)
+      if (Index==-1) then
+        call create_header_mpi(name,field,nind)
+	Index=GetFieldIndex(name)
+      else
+        FieldIndex(Index)=FieldIndex(Index)+1.
+      endif
+      
+      IF (Field(ind_b)%field_type==field_T) THEN
+        ncell=1
+        DO ind=ind_b,ind_e
+          d=>domain(ind)
+          IF (Field(ind)%field_type/=field_T) THEN
+            PRINT *,"Writefield, grille non geree"
+            RETURN
+          ENDIF
+
+        n=0
+        DO j=d%jj_begin-halo_size,d%jj_end+halo_size
+          DO i=d%ii_begin-halo_size,d%ii_end+halo_size
+            IF (d%own(i,j) .OR. single) n=n+1
+          ENDDO
+        ENDDO
+
+        displ=FieldVarId(index)%displ
+        
+        IF (field(ind)%ndim==2) THEN
+          ALLOCATE(Field_val2d(n))
+          n=0
+          DO j=d%jj_begin-halo_size,d%jj_end+halo_size
+            DO i=d%ii_begin-halo_size,d%ii_end+halo_size
+              k=d%iim*(j-1)+i
+              IF (d%own(i,j) .OR. single) THEN
+                n=n+1
+                Field_val2d(n)=field(ind)%rval2d(k)
+              ENDIF
+            ENDDO
+          ENDDO
+          status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val2d,  &
+                              start=(/ displ+ncell,FieldIndex(Index) /),count=(/n,1 /))
+          DEALLOCATE(field_val2d)
+        ELSE IF (field(ind)%ndim==3) THEN
+          ALLOCATE(Field_val3d(n,size(field(ind)%rval3d,2)))
+          n=0
+          DO j=d%jj_begin-halo_size,d%jj_end+halo_size
+            DO i=d%ii_begin-halo_size,d%ii_end+halo_size
+              k=d%iim*(j-1)+i
+              IF (d%own(i,j) .OR. single) THEN
+                n=n+1
+                Field_val3d(n,:)=field(ind)%rval3d(k,:)
+              ENDIF
+            ENDDO
+          ENDDO
+          status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,  &
+                                start=(/ displ+ncell,1,FieldIndex(Index) /), count=(/n,size(field(ind)%rval3d,2),1 /))
+          DEALLOCATE(field_val3d)
+
+        ELSE IF (field(1)%ndim==4) THEN
+
+          DO q=1,FieldVarId(index)%size
+          
+            ALLOCATE(Field_val3d(n,size(field(ind)%rval4d,2)))
+            n=0
+            DO j=d%jj_begin-halo_size,d%jj_end+halo_size
+              DO i=d%ii_begin-halo_size,d%ii_end+halo_size
+                k=d%iim*(j-1)+i
+                IF (d%own(i,j) .OR. single) THEN
+                  n=n+1
+                  Field_val3d(n,:)=field(ind)%rval4d(k,:,q)
+                ENDIF
+              ENDDO
+             ENDDO
+
+            status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d(:,l), &
+                                start=(/ displ+ncell,l,FieldIndex(Index) /), count=(/n,size(field(ind)%rval4d,2),1 /))
+            DEALLOCATE(field_val3d)
+          ENDDO         
+        ENDIF
+        
+        ncell=ncell+n
+
+     ENDDO
+     
+     ELSE IF (Field(ind_b)%field_type==field_Z) THEN
+        ncell=1
+        n=0
+        DO ind=ind_b,ind_e
+          d=>domain(ind)
+          CALL swap_geometry(ind)
+          CALL swap_dimensions(ind)
+  
+          n=0
+          DO j=jj_begin+1,jj_end
+            DO i=ii_begin,ii_end-1
+              n=n+1
+            ENDDO
+          ENDDO
+
+          DO j=jj_begin,jj_end-1
+            DO i=ii_begin+1,ii_end
+              n=n+1
+            ENDDO
+          ENDDO
+
+        displ=FieldVarId(index)%displ
+
+        IF (field(ind)%ndim==2) THEN
+          ALLOCATE(Field_val2d(n))
+
+          n=0
+          DO j=jj_begin+1,jj_end
+            DO i=ii_begin,ii_end-1
+              n=n+1
+              k=iim*(j-1)+i
+              Field_val2d(n)=field(ind)%rval2d(k+z_down)
+            ENDDO
+          ENDDO
+
+          DO j=jj_begin,jj_end-1
+            DO i=ii_begin+1,ii_end
+              n=n+1
+              k=iim*(j-1)+i
+              Field_val2d(n)=field(ind)%rval2d(k+z_up)
+            ENDDO
+          ENDDO
+
+          status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),                       &
+                              Field_val2d,start=(/ displ+ncell,FieldIndex(Index) /),count=(/n,1 /))
+          DEALLOCATE(field_val2d)
+
+        ELSE IF (field(ind)%ndim==3) THEN
+          ALLOCATE(Field_val3d(n,size(field(ind)%rval3d,2)))
+          n=0
+          DO j=jj_begin+1,jj_end
+            DO i=ii_begin,ii_end-1
+              n=n+1
+              k=iim*(j-1)+i
+              Field_val3d(n,:)=field(ind)%rval3d(k+z_down,:)
+            ENDDO
+          ENDDO
+
+          DO j=jj_begin,jj_end-1
+            DO i=ii_begin+1,ii_end
+              n=n+1
+              k=iim*(j-1)+i
+              Field_val3d(n,:)=field(ind)%rval3d(k+z_up,:)
+            ENDDO
+          ENDDO
+
+           status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,    &
+                               start=(/ displ+ncell,1,FieldIndex(Index) /), count=(/n,size(field(ind)%rval3d,2),1 /))
+          DEALLOCATE(field_val3d)
+
+        ELSE IF (field(1)%ndim==4) THEN
+
+          DO q=1,FieldVarId(index)%size
+            ALLOCATE(Field_val3d(n,size(field(ind)%rval4d,2)))
+            n=0
+            DO j=jj_begin+1,jj_end
+              DO i=ii_begin,ii_end-1
+                n=n+1
+                k=iim*(j-1)+i
+                Field_val3d(n,:)=field(ind)%rval4d(k+z_down,:,q)
+              ENDDO
+            ENDDO
+
+            DO j=jj_begin,jj_end-1
+              DO i=ii_begin+1,ii_end
+                n=n+1
+                k=iim*(j-1)+i
+                Field_val3d(n,:)=field(ind)%rval4d(k+z_up,:,q)
+              ENDDO
+            ENDDO
+
+            status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,  &
+                                start=(/ displ+ncell,1,FieldIndex(Index) /), count=(/n,size(field(ind)%rval4d,2),1 /))
+            DEALLOCATE(field_val3d)
+          ENDDO
+        ENDIF
+        
+        ncell=ncell+n
+
+     ENDDO
+     
+     ENDIF
+     status=NF90_SYNC(FieldId(Index))
+      
+    END SUBROUTINE Writefield_mpi
+    
+    
+    SUBROUTINE Create_header_mpi(name,field,nind)
+    USE netcdf_mod
+    USE field_mod
+    USE domain_mod
+    USE spherical_geom_mod
+    USE dimensions
+    USE geometry
+    USE mpi_mod
+    USE mpipara
+    IMPLICIT NONE
+      CHARACTER(LEN=*) :: name
+      CHARACTER(LEN=LEN_TRIM(ADJUSTL(name))) :: name_adj
+      TYPE(t_field),POINTER :: field(:)
+      INTEGER,OPTIONAL,INTENT(IN) :: nind
+      INTEGER :: ncell
+      INTEGER :: nvert
+      REAL(rstd),ALLOCATABLE :: lon(:),lat(:),bounds_lon(:,:),bounds_lat(:,:)
+      TYPE(t_domain),POINTER :: d
+      INTEGER :: nvertId,ncid,lonId,latId,bounds_lonId,bounds_latId,timeId,ncellId
+      INTEGER :: dim3id,dim4id
+      INTEGER :: status
+      INTEGER :: ind,i,j,k,n,q
+      INTEGER :: iie,jje,iin,jjn
+      INTEGER :: ind_b,ind_e
+      INTEGER :: halo_size
+      LOGICAL :: single 
+      INTEGER :: nij
+      INTEGER :: ncell_glo(0:mpi_size-1)
+      INTEGER :: displ, ncell_tot
+      
+          
+      NbField=NbField+1
+      name_adj=TRIM(ADJUSTL(name))  ! work around ICE with pgf90
+      FieldName(NbField)=name_adj
+      FieldIndex(NbField)=1
+      
+      IF (PRESENT(nind)) THEN
+        ind_b=nind
+        ind_e=nind
+        halo_size=1
+        single=.TRUE.
+      ELSE
+        ind_b=1
+        ind_e=ndomain
+        halo_size=0
+        single=.FALSE.
+      ENDIF
+      
+      ncell=0
+      
+      IF (Field(ind_b)%field_type==field_T) THEN
+        nvert=6
+        
+        DO ind=ind_b,ind_e
+          d=>domain(ind)
+        
+          DO j=d%jj_begin-halo_size,d%jj_end+halo_size
+            DO i=d%ii_begin-halo_size,d%ii_end+halo_size
+              IF (single .OR. domain(ind)%own(i,j)) ncell=ncell+1
+            ENDDO
+          ENDDO
+
+        END DO
+      
+        CALL MPI_ALLGATHER(ncell,1,MPI_INTEGER,ncell_glo,1,MPI_INTEGER,comm_icosa,ierr)
+
+        displ=0
+        DO i=1,mpi_rank
+          displ=displ+ncell_glo(i-1)
+        ENDDO
+        FieldVarId(NbField)%displ=displ
+        ncell_tot=sum(ncell_glo(:))
+        
+!        status = NF90_CREATE_PAR(TRIM(ADJUSTL(name))//'.nc', IOR(NF90_NETCDF4, NF90_MPIIO), comm_icosa, MPI_INFO_NULL, ncid)
+        FieldId(NbField)=ncid
+        status = NF90_DEF_DIM(ncid,'cell',ncell_tot,ncellId)
+        status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid)
+
+        IF (Field(ind_b)%ndim==2)  THEN
+          FieldVarId(NbField)%size=1
+          ALLOCATE(FieldVarId(NbField)%nc_id(1))
+        ELSE IF (Field(ind_b)%ndim==3)  THEN
+          FieldVarId(NbField)%size=1
+          ALLOCATE(FieldVarId(NbField)%nc_id(1))
+          status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval3d,2),dim3id)
+        ELSE IF (Field(1)%ndim==4) THEN
+          FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3)
+          ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size))
+          status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval4d,2),dim3id)
+!          status = NF90_DEF_DIM(ncid,'Q',size(field(ind_b)%rval4d,3),dim4id)
+        ENDIF
+      
+        status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId)
+      
+        status = NF90_DEF_VAR(ncid,'lon',NF90_DOUBLE,(/ ncellId /),lonId)
+        status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude")
+        status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east")
+        status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon")
+        status = NF90_DEF_VAR(ncid,'lat',NF90_DOUBLE,(/ ncellId /),latId)
+        status = NF90_PUT_ATT(ncid,latId,"long_name","latitude")
+        status = NF90_PUT_ATT(ncid,latId,"units","degrees_north")
+        status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat")
+        status = NF90_DEF_VAR(ncid,'bounds_lon',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId)
+        status = NF90_DEF_VAR(ncid,'bounds_lat',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId)
+
+        IF (Field(ind_b)%ndim==2) THEN
+          status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1))
+          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat")
+          status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(1), NF90_CHUNKED, (/ncell_tot,1/))
+        ELSE IF (Field(ind_b)%ndim==3) THEN
+          status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1))
+          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat")
+          status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(1), NF90_CHUNKED,   &
+                                         (/ncell_tot,size(field(ind_b)%rval3d,2),1/))
+        ELSE IF (Field(ind_b)%ndim==4) THEN
+          DO i=1,FieldVarId(NbField)%size
+            status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name))//int2str(i),ncprec,(/ ncellId,dim3id,timeId /),  &
+                                  FieldVarId(NbField)%nc_id(i))
+            status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(i),"coordinates","lon lat")
+            status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(q), NF90_CHUNKED,   &
+                                           (/ncell_tot,size(field(ind_b)%rval4d,2),1/))
+          ENDDO        
+        ENDIF 
+  
+      
+        status = NF90_ENDDEF(ncid)      
+
+        ncell=1
+        DO ind=ind_b,ind_e
+          d=>domain(ind)
+  
+          n=0
+          DO j=d%jj_begin-halo_size,d%jj_end+halo_size
+            DO i=d%ii_begin-halo_size,d%ii_end+halo_size
+              IF (single .OR. d%own(i,j)) n=n+1
+            ENDDO
+          ENDDO
+          
+         ALLOCATE(lon(n),lat(n),bounds_lon(0:nvert-1,n),bounds_lat(0:nvert-1,n))
+          
+          n=0  
+          DO j=d%jj_begin-halo_size,d%jj_end+halo_size
+            DO i=d%ii_begin-halo_size,d%ii_end+halo_size
+                IF (d%own(i,j) .OR. single) THEN
+                n=n+1
+                CALL xyz2lonlat(d%xyz(:,i,j),lon(n),lat(n))
+                lon(n)=lon(n)*180/Pi
+                lat(n)=lat(n)*180/Pi
+                DO k=0,5
+                  CALL xyz2lonlat(d%vertex(:,k,i,j),bounds_lon(k,n), bounds_lat(k,n))
+                  bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
+                  bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
+                ENDDO
+              ENDIF
+            ENDDO
+          ENDDO
+          status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ displ+ncell /),count=(/ n /))
+          status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ displ+ncell /),count=(/ n /))
+          status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,displ+ncell /),count=(/ nvert,n /))
+          status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,displ+ncell /),count=(/ nvert,n /))
+  
+          ncell=ncell+n
+          DEALLOCATE(lon,lat,bounds_lon,bounds_lat)
+      END DO 
+
+    ELSE IF (Field(ind_b)%field_type==field_Z) THEN
+        nvert=3
+        DO ind=ind_b,ind_e
+          d=>domain(ind)
+        
+          DO j=d%jj_begin+1,d%jj_end
+            DO i=d%ii_begin,d%ii_end-1
+              ncell=ncell+1
+            ENDDO
+          ENDDO
+
+          DO j=d%jj_begin,d%jj_end-1
+            DO i=d%ii_begin+1,d%ii_end
+              ncell=ncell+1
+            ENDDO
+          ENDDO
+
+        END DO
+        
+        CALL MPI_ALLGATHER(ncell,1,MPI_INTEGER,ncell_glo,1,MPI_INTEGER,comm_icosa,ierr)
+
+        displ=0
+        DO i=1,mpi_rank
+          displ=displ+ncell_glo(i-1)
+        ENDDO
+        FieldVarId(NbField)%displ=displ
+        ncell_tot=sum(ncell_glo(:))
+              
+!        status = NF90_CREATE_PAR(TRIM(ADJUSTL(name))//'.nc',IOR(NF90_NETCDF4, NF90_MPIIO), comm_icosa, MPI_INFO_NULL, ncid)
+!        status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid)
+        FieldId(NbField)=ncid
+        status = NF90_DEF_DIM(ncid,'cell',ncell_tot,ncellId)
+        status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid)
+
+        IF (Field(ind_b)%ndim==2)  THEN
+          FieldVarId(NbField)%size=1
+          ALLOCATE(FieldVarId(NbField)%nc_id(1))
+        ELSE IF (Field(ind_b)%ndim==3)  THEN
+          FieldVarId(NbField)%size=1
+          ALLOCATE(FieldVarId(NbField)%nc_id(1))
+          status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval3d,2),dim3id)
+        ELSE IF (Field(1)%ndim==4) THEN
+          FieldVarId(NbField)%size=size(field(ind_b)%rval4d,3)
+          ALLOCATE(FieldVarId(NbField)%nc_id(FieldVarId(NbField)%size))
+          status = NF90_DEF_DIM(ncid,'lev',size(field(ind_b)%rval4d,2),dim3id)
+        ENDIF
+
+
+      
+        status = NF90_DEF_DIM(ncid,'time_counter',NF90_UNLIMITED,timeId)
+      
+        status = NF90_DEF_VAR(ncid,'lon',NF90_DOUBLE,(/ ncellId /),lonId)
+        status = NF90_PUT_ATT(ncid,lonId,"long_name","longitude")
+        status = NF90_PUT_ATT(ncid,lonId,"units","degrees_east")
+        status = NF90_PUT_ATT(ncid,lonId,"bounds","bounds_lon")
+        status = NF90_DEF_VAR(ncid,'lat',NF90_DOUBLE,(/ ncellId /),latId)
+        status = NF90_PUT_ATT(ncid,latId,"long_name","latitude")
+        status = NF90_PUT_ATT(ncid,latId,"units","degrees_north")
+        status = NF90_PUT_ATT(ncid,latId,"bounds","bounds_lat")
+        status = NF90_DEF_VAR(ncid,'bounds_lon',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_lonId)
+        status = NF90_DEF_VAR(ncid,'bounds_lat',NF90_DOUBLE,(/ nvertid,ncellId /),bounds_latId)
+
+
+        IF (Field(ind_b)%ndim==2) THEN
+          status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,timeId /),FieldVarId(NbField)%nc_id(1))
+          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat")
+          status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(1), NF90_CHUNKED, (/ncell_tot,1/))
+        ELSE IF (Field(ind_b)%ndim==3) THEN
+          status = NF90_DEF_VAR(ncid,TRIM(ADJUSTL(name)),ncprec,(/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(1))
+          status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(1),"coordinates","lon lat")
+          status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(1), NF90_CHUNKED,   &
+                                         (/ncell_tot,size(field(ind_b)%rval3d,2),1/))
+        ELSE IF (Field(ind_b)%ndim==4) THEN
+          DO q=1,FieldVarId(NbField)%size
+            status = NF90_DEF_VAR(ncid,name_adj//int2str(q),ncprec,             &
+                                  (/ ncellId,dim3id,timeId /),FieldVarId(NbField)%nc_id(q))
+            status = NF90_PUT_ATT(ncid,FieldVarId(NbField)%nc_id(q),"coordinates","lon lat")
+           status = NF90_DEF_VAR_CHUNKING(ncid, FieldVarId(NbField)%nc_id(q), NF90_CHUNKED, &
+                                          (/ncell_tot,size(field(ind_b)%rval4d,2),1/))
+          ENDDO        
+        ENDIF 
+        
+        status = NF90_ENDDEF(ncid)      
+
+        ncell=1
+        DO ind=ind_b,ind_e
+          d=>domain(ind)
+          CALL swap_geometry(ind)
+          CALL swap_dimensions(ind)
+  
+          n=0
+          DO j=jj_begin+1,jj_end
+            DO i=ii_begin,ii_end-1
+              n=n+1
+            ENDDO
+          ENDDO
+
+          DO j=jj_begin,jj_end-1
+            DO i=ii_begin+1,ii_end
+              n=n+1
+            ENDDO
+          ENDDO
+
+         ALLOCATE(lon(n),lat(n),bounds_lon(0:nvert-1,n),bounds_lat(0:nvert-1,n))
+          
+          n=0  
+       
+          DO j=jj_begin+1,jj_end
+            DO i=ii_begin,ii_end-1
+              nij=(j-1)*iim+i
+              n=n+1
+              CALL xyz2lonlat(xyz_v(nij+z_down,:)/radius,lon(n),lat(n))
+              lon(n)=lon(n)*180/Pi
+              lat(n)=lat(n)*180/Pi
+              CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n))
+              CALL xyz2lonlat(xyz_i(nij+t_ldown,:)/radius,bounds_lon(1,n), bounds_lat(1,n))
+              CALL xyz2lonlat(xyz_i(nij+t_rdown,:)/radius,bounds_lon(2,n), bounds_lat(2,n))
+
+              DO k=0,2
+                bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
+                bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
+              ENDDO
+            ENDDO
+          ENDDO
+
+          DO j=jj_begin,jj_end-1
+            DO i=ii_begin+1,ii_end
+              nij=(j-1)*iim+i
+              n=n+1
+              CALL xyz2lonlat(xyz_v(nij+z_up,:)/radius,lon(n),lat(n))
+              lon(n)=lon(n)*180/Pi
+              lat(n)=lat(n)*180/Pi
+              CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n))
+              CALL xyz2lonlat(xyz_i(nij+t_rup,:)/radius,bounds_lon(1,n), bounds_lat(1,n))
+              CALL xyz2lonlat(xyz_i(nij+t_lup,:)/radius,bounds_lon(2,n), bounds_lat(2,n))
+
+              DO k=0,2
+                bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
+                bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
+              ENDDO
+            ENDDO
+          ENDDO
+          
+          
+          status=NF90_PUT_VAR(ncid,lonid,REAL(lon,r8),start=(/ displ+ncell /),count=(/ n /))
+          status=NF90_PUT_VAR(ncid,latid,REAL(lat,r8),start=(/ displ+ncell /),count=(/ n /))
+          status=NF90_PUT_VAR(ncid,bounds_lonId,REAL(bounds_lon,r8),start=(/ 1,displ+ncell /),count=(/ nvert,n /))
+          status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,displ+ncell /),count=(/ nvert,n /))
+  
+          ncell=ncell+n
+          DEALLOCATE(lon,lat,bounds_lon,bounds_lat)
+      END DO          
+    ENDIF
+
+
+   END SUBROUTINE Create_Header_mpi  
+   
+ end module write_field_mpi_mod
+  
Index: codes/icosagcm/devel/src/output/write_field_vars.f90
===================================================================
--- codes/icosagcm/devel/src/output/write_field_vars.f90	(revision 880)
+++ codes/icosagcm/devel/src/output/write_field_vars.f90	(revision 880)
@@ -0,0 +1,60 @@
+module write_field_vars_mod
+  USE genmod
+  IMPLICIT NONE
+
+  INTEGER,SAVE :: ncprec
+
+  TYPE ncvar
+     INTEGER :: size
+     INTEGER,POINTER :: nc_id(:)
+     INTEGER :: displ
+  END TYPE ncvar
+
+  INTEGER, PARAMETER :: MaxWriteField = 1000
+  INTEGER, DIMENSION(MaxWriteField),SAVE :: FieldId
+  TYPE(ncvar), dimension(MaxWriteField),SAVE :: FieldVarId
+  INTEGER, DIMENSION(MaxWriteField),SAVE :: FieldIndex
+  CHARACTER(len=255), DIMENSION(MaxWriteField) ::  FieldName 
+
+  INTEGER,SAVE :: NbField = 0
+
+CONTAINS
+
+  FUNCTION GetFieldIndex(name)
+    integer          :: GetFieldindex
+    character(len=*) :: name
+
+    character(len=255) :: TrueName
+    integer            :: i
+
+
+    TrueName=TRIM(ADJUSTL(name))
+
+    GetFieldIndex=-1
+    do i=1,NbField
+       if (TrueName==FieldName(i)) then
+          GetFieldIndex=i
+          exit
+       endif
+    enddo
+  END FUNCTION GetFieldIndex
+
+  function int2str(int)
+    integer, parameter :: MaxLen=10
+    integer,intent(in) :: int
+    character(len=MaxLen) :: int2str
+    logical :: flag
+    integer :: i
+    flag=.true.
+
+    i=int
+
+    int2str=''
+    do while (flag)
+       int2str=CHAR(MOD(i,10)+48)//int2str
+       i=i/10
+       if (i==0) flag=.false.
+    enddo
+  end function int2str
+    
+end module write_field_vars_mod
Index: codes/icosagcm/devel/src/output/xios_mod.F90
===================================================================
--- codes/icosagcm/devel/src/output/xios_mod.F90	(revision 879)
+++ codes/icosagcm/devel/src/output/xios_mod.F90	(revision 880)
@@ -1,67 +1,63 @@
 MODULE xios_mod
-
+  
 #ifdef CPP_USING_XIOS
-   USE xios
+  USE xios
 #endif
-   IMPLICIT NONE 
-
-  PUBLIC
-  LOGICAL,SAVE :: using_xios
-
-  INTEGER,SAVE :: ncell_i
-!$OMP THREADPRIVATE(ncell_i)
-  INTEGER,SAVE :: ncell_v
-!$OMP THREADPRIVATE(ncell_v)
-  INTEGER,SAVE :: ncell_e
-!$OMP THREADPRIVATE(ncell_e)
-
-  PRIVATE ncell_i,ncell_v,ncell_e
+
+  USE prec, ONLY : rstd
+  USE field_mod, ONLY : t_field, field_T, field_U, field_Z
+  USE domain_mod, ONLY : t_domain, t_cellset, domain, ndomain, mesh_loc
+  IMPLICIT NONE   
+  PRIVATE
+  SAVE
+
+  LOGICAL :: using_xios
+
+  INTEGER :: ncell_i, ncell_v, ncell_e
+!$OMP THREADPRIVATE(ncell_i, ncell_v, ncell_e)
+
+  PUBLIC :: using_xios, xios_init, &
+       xios_init_write_field,  xios_write_field_finalize, &
+       xios_write_field, xios_read_field
 
 #ifdef CPP_USING_XIOS
+
+  PUBLIC :: xios_timestep,    & 
+       xios_set_file_attr, xios_set_fieldgroup_attr, &
+       xios_set_filegroup_attr, xios_get_axis_attr, &
+       xios_send_field, xios_read_var, &
+       xios_update_calendar, xios_set_context
   
 CONTAINS
  
  SUBROUTINE xios_init
-   USE getin_mod
-   USE xios
-   USE mpipara
-   IMPLICIT NONE
-    TYPE(xios_context) :: ctx_hdl
-
-     using_xios=.TRUE.
-     CALL xios_context_initialize("icosagcm",comm_icosa)
-     CALL xios_get_handle("icosagcm",ctx_hdl)
-     CALL xios_set_current_context(ctx_hdl)
-   
- END SUBROUTINE xios_init  
- 
+   USE mpipara, ONLY : comm_icosa
+   TYPE(xios_context) :: ctx_hdl
+
+   using_xios=.TRUE.
+   CALL xios_context_initialize("icosagcm",comm_icosa)
+   CALL xios_get_handle("icosagcm",ctx_hdl)
+   CALL xios_set_current_context(ctx_hdl)   
+
+ END SUBROUTINE xios_init
+
  SUBROUTINE xios_init_write_field
- USE genmod
- USE mpipara
- USE xios
- USE grid_param
- USE domain_mod
- USE dimensions
- USE spherical_geom_mod
- USE geometry
- USE mpi_mod
- USE time_mod
- USE metric, ONLY : vup,vdown, cell_glo
- USE disvert_mod, ONLY : presnivs
- IMPLICIT NONE
-  TYPE(xios_context) :: ctx_hdl
-  TYPE(xios_duration)      :: dtime
-  REAL(rstd) :: lev_value(llm)
-  REAL(rstd) :: lev_valuep1(llm+1)
-  REAL(rstd) :: nq_value(nqtot)
-  INTEGER :: ncell, ncell_tot, ncell_glo(0:mpi_size-1), displ
-  INTEGER :: ind, i,j,k,l,ij
-  REAL(rstd),ALLOCATABLE    :: lon(:), lat(:), bounds_lon(:,:), bounds_lat(:,:)
-  INTEGER, ALLOCATABLE      :: ind_glo(:)
-  TYPE(t_domain),POINTER :: d
-
-!$OMP BARRIER
-!$OMP MASTER
-!   CALL xios_context_initialize("icosagcm",comm_icosa)
+   USE time_mod, ONLY : dt, itau_out
+   USE grid_param, ONLY : llm, nqtot
+   USE mpi_mod, ONLY : MPI_INTEGER
+
+   TYPE(xios_context) :: ctx_hdl
+   TYPE(xios_duration)      :: dtime
+   REAL(rstd) :: lev_value(llm)
+   REAL(rstd) :: lev_valuep1(llm+1)
+   REAL(rstd) :: nq_value(nqtot)
+   INTEGER :: l, ncell, ncell_tot, displ
+   REAL(rstd),ALLOCATABLE    :: lon(:), lat(:), bounds_lon(:,:), bounds_lat(:,:)
+   INTEGER, ALLOCATABLE      :: ind_glo(:)
+   TYPE(t_domain),POINTER :: d
+   
+   !$OMP BARRIER
+   !$OMP MASTER
+   !   CALL xios_context_initialize("icosagcm",comm_icosa)
    CALL xios_get_handle("icosagcm",ctx_hdl)
    CALL xios_set_current_context(ctx_hdl)
@@ -71,262 +67,102 @@
    CALL xios_set_axis_attr("lev",n_glo=llm ,value=lev_value) ;
    CALL xios_set_axis_attr("levp1",n_glo=llm+1 ,value=lev_valuep1) ;
-   CALL xios_set_axis_attr("presnivs_mb",n_glo=llm, value=presnivs/100., unit="mb") ; 
    CALL xios_set_axis_attr("nq",n_glo=nqtot, value=nq_value) ;
    
-   ncell=0
-   DO ind=1,ndomain
-     d=>domain(ind)
-        
-     DO j=d%jj_begin,d%jj_end
-       DO i=d%ii_begin,d%ii_end
-         IF (domain(ind)%own(i,j)) ncell=ncell+1
-       ENDDO
-     ENDDO
-   ENDDO      
-   ncell_i=ncell
-   
-   CALL MPI_ALLGATHER(ncell,1,MPI_INTEGER,ncell_glo,1,MPI_INTEGER,comm_icosa,ierr)
-
-   displ=0
-   DO i=1,mpi_rank
-     displ=displ+ncell_glo(i-1)
-   ENDDO
-
-   ncell_tot=sum(ncell_glo(:))
-   
-   ALLOCATE(lon(ncell), lat(ncell), bounds_lon(0:5,ncell), bounds_lat(0:5,ncell), ind_glo(ncell)) 
-   
-   ncell=0
-   DO ind=1,ndomain
-     d=>domain(ind)
-        
-     DO j=d%jj_begin,d%jj_end
-       DO i=d%ii_begin,d%ii_end
-         IF (domain(ind)%own(i,j)) THEN 
-           ncell=ncell+1
-           CALL xyz2lonlat(d%xyz(:,i,j),lon(ncell),lat(ncell))
-           lon(ncell)=lon(ncell)*180/Pi
-           lat(ncell)=lat(ncell)*180/Pi
-           DO k=0,5
-             CALL xyz2lonlat(d%vertex(:,k,i,j),bounds_lon(k,ncell), bounds_lat(k,ncell))
-             bounds_lat(k,ncell)=bounds_lat(k,ncell)*180/Pi
-             bounds_lon(k,ncell)=bounds_lon(k,ncell)*180/Pi
-           ENDDO
-           ind_glo(ncell)=domain(ind)%assign_cell_glo(i,j)-1 
-         ENDIF
-       ENDDO
-     ENDDO
-   ENDDO         
-
+   !------------------ primal cells ------------------
+   CALL collect_bounds(6, mesh_loc%primal_own)
+   ncell_i=ncell   
    CALL xios_set_domaingroup_attr("i",ni_glo=ncell_tot, ibegin=displ, ni=ncell)
    CALL xios_set_domaingroup_attr("i", data_dim=1, type='unstructured' , nvertex=6, i_index=ind_glo)
    CALL xios_set_domaingroup_attr("i",lonvalue_1d=lon, latvalue_1d=lat, bounds_lon_1d=bounds_lon, bounds_lat_1d=bounds_lat)
-   
    DEALLOCATE(lon, lat, bounds_lon, bounds_lat,ind_glo) 
-    
-
-
-   ncell=0
-   DO ind=1,ndomain
-     d=>domain(ind)
-
-     DO j=d%jj_begin,d%jj_end
-       DO i=d%ii_begin,d%ii_end
-         DO k=0,5
-            IF (d%edge_assign_domain(k,i,j)==domloc_glo_ind(ind) .AND. d%edge_assign_i(k,i,j)==i .AND. d%edge_assign_j(k,i,j)==j &
-                                                                                 .AND. d%edge_assign_pos(k,i,j)==k) THEN
-               ncell=ncell+1
-            ENDIF
-         ENDDO
-       ENDDO
-     ENDDO
-   ENDDO
+   
+   !--------------------- dual cells ------------------
+   
+   CALL collect_bounds(3, mesh_loc%dual_own)
+   ncell_v=ncell
+   CALL xios_set_domain_attr("v",ni_glo=ncell_tot, ibegin=displ, ni=ncell)
+   CALL xios_set_domain_attr("v", data_dim=1, type='unstructured' , nvertex=3)
+   CALL xios_set_domain_attr("v",lonvalue_1d=lon, latvalue_1d=lat, bounds_lon_1d=bounds_lon, bounds_lat_1d=bounds_lat)
+   DEALLOCATE(lon, lat, bounds_lon, bounds_lat,ind_glo) 
+   
+   !---------------------- edges -----------------------
+   
+   CALL collect_bounds(2, mesh_loc%edge_own)
    ncell_e=ncell
-   
-   CALL MPI_ALLGATHER(ncell_e,1,MPI_INTEGER,ncell_glo,1,MPI_INTEGER,comm_icosa,ierr)
-   displ=0
-   DO i=1,mpi_rank
-     displ=displ+ncell_glo(i-1)
-   ENDDO
-   ncell_tot=sum(ncell_glo(:))
-   
-   ALLOCATE(lon(ncell), lat(ncell), bounds_lon(0:1,ncell), bounds_lat(0:1,ncell),ind_glo(ncell)) 
-
-
-   ncell=0
-   DO ind=1,ndomain
-     d=>domain(ind)
-     CALL swap_dimensions(ind)
-     CALL swap_geometry(ind)
-
-     DO j=d%jj_begin,d%jj_end
-       DO i=d%ii_begin,d%ii_end
-         DO k=0,5
-           IF (d%edge_assign_domain(k,i,j)==domloc_glo_ind(ind) .AND. d%edge_assign_i(k,i,j)==i .AND. d%edge_assign_j(k,i,j)==j &
-                                                                                 .AND. d%edge_assign_pos(k,i,j)==k) THEN
-              ncell=ncell+1
-              ij=(j-1)*iim+i
-
-              lon(ncell)=lon_e(ij+u_pos(k+1))*180/Pi
-              lat(ncell)=lat_e(ij+u_pos(k+1))*180/Pi
-               
-              CALL xyz2lonlat(d%vertex(:,MOD((k-1)+6,6),i,j),bounds_lon(0,ncell), bounds_lat(0,ncell))
-              CALL xyz2lonlat(d%vertex(:,k,i,j),bounds_lon(1,ncell), bounds_lat(1,ncell))
-              bounds_lon(:,ncell)=bounds_lon(:,ncell)*180/Pi
-              bounds_lat(:,ncell)=bounds_lat(:,ncell)*180/Pi
-              ind_glo(ncell)=cell_glo(d%assign_cell_glo(i,j))%edge(MOD(k+d%delta(i,j)+6,6))-1 
-           ENDIF               
-         ENDDO
-       ENDDO
-     ENDDO
-   ENDDO
    CALL xios_set_domain_attr("u",ni_glo=ncell_tot, ibegin=displ, ni=ncell)
    CALL xios_set_domain_attr("u", data_dim=1, type='unstructured' , nvertex=2, i_index=ind_glo)
    CALL xios_set_domain_attr("u",lonvalue_1d=lon, latvalue_1d=lat, bounds_lon_1d=bounds_lon, bounds_lat_1d=bounds_lat)
-
    DEALLOCATE(lon, lat, bounds_lon, bounds_lat,ind_glo) 
-
-
-   ncell=0
-   DO ind=1,ndomain
-     d=>domain(ind)
-        
-     DO j=d%jj_begin+1,d%jj_end
-       DO i=d%ii_begin,d%ii_end-1
-         ncell=ncell+1
-       ENDDO
-     ENDDO
-
-     DO j=d%jj_begin,d%jj_end-1
-       DO i=d%ii_begin+1,d%ii_end
-          ncell=ncell+1
-        ENDDO
-     ENDDO
-
-   ENDDO      
-   
-   ncell_v=ncell
-   
-   CALL MPI_ALLGATHER(ncell,1,MPI_INTEGER,ncell_glo,1,MPI_INTEGER,comm_icosa,ierr)
-
-   displ=0
-   DO i=1,mpi_rank
-     displ=displ+ncell_glo(i-1)
-   ENDDO
-
-   ncell_tot=sum(ncell_glo(:))
-   
-   ALLOCATE(lon(ncell), lat(ncell), bounds_lon(0:2,ncell), bounds_lat(0:2,ncell)) 
-   
-   ncell=0
-   DO ind=1,ndomain
-     d=>domain(ind)
- 
-     DO j=d%jj_begin+1,d%jj_end
-       DO i=d%ii_begin,d%ii_end-1
-           ncell=ncell+1
-           CALL xyz2lonlat(d%vertex(:,vdown,i,j),lon(ncell),lat(ncell))
-           lon(ncell)=lon(ncell)*180/Pi
-           lat(ncell)=lat(ncell)*180/Pi
-
-           CALL xyz2lonlat(d%xyz(:,i,j),bounds_lon(0,ncell), bounds_lat(0,ncell))
-           CALL xyz2lonlat(d%xyz(:,i,j-1),bounds_lon(1,ncell), bounds_lat(1,ncell))
-           CALL xyz2lonlat(d%xyz(:,i+1,j-1),bounds_lon(2,ncell), bounds_lat(2,ncell))
-
-           DO k=0,2
-             bounds_lat(k,ncell)=bounds_lat(k,ncell)*180/Pi
-             bounds_lon(k,ncell)=bounds_lon(k,ncell)*180/Pi
-           ENDDO
-         ENDDO
-       ENDDO 
- 
-       DO j=d%jj_begin,d%jj_end-1
-         DO i=d%ii_begin+1,d%ii_end
-           ncell=ncell+1
-           CALL xyz2lonlat(d%vertex(:,vup,i,j),lon(ncell),lat(ncell))
-           lon(ncell)=lon(ncell)*180/Pi
-           lat(ncell)=lat(ncell)*180/Pi
-           CALL xyz2lonlat(d%xyz(:,i,j),bounds_lon(0,ncell), bounds_lat(0,ncell))
-           CALL xyz2lonlat(d%xyz(:,i,j+1),bounds_lon(1,ncell), bounds_lat(1,ncell))
-           CALL xyz2lonlat(d%xyz(:,i-1,j+1),bounds_lon(2,ncell), bounds_lat(2,ncell))
-
-           DO k=0,2
-             bounds_lat(k,ncell)=bounds_lat(k,ncell)*180/Pi
-             bounds_lon(k,ncell)=bounds_lon(k,ncell)*180/Pi
-           ENDDO
-         ENDDO
-       ENDDO 
-        
-   ENDDO         
-
-  
-   CALL xios_set_domain_attr("v",ni_glo=ncell_tot, ibegin=displ, ni=ncell)
-   CALL xios_set_domain_attr("v", data_dim=1, type='unstructured' , nvertex=3)
-   CALL xios_set_domain_attr("v",lonvalue_1d=lon, latvalue_1d=lat, bounds_lon_1d=bounds_lon, bounds_lat_1d=bounds_lat)
-
-
+   
    dtime%second=dt
    CALL xios_set_timestep(dtime)
    
    CALL xios_set_fieldgroup_attr("standard_output", freq_op=itau_out*xios_timestep, freq_offset=(itau_out-1)*xios_timestep)
-    
+   
    CALL xios_close_context_definition()
-!$OMP END MASTER
-!$OMP BARRIER
+   !$OMP END MASTER
+   !$OMP BARRIER
+   
+ CONTAINS
+   
+   SUBROUTINE collect_bounds(nvert, cells)
+     USE mpipara, ONLY : comm_icosa, mpi_size, mpi_rank
+     INTEGER, INTENT(IN) :: nvert
+     TYPE(t_cellset)     :: cells(:)
+     INTEGER :: i, ind, n_beg, n_end, ierr, ncell_glo(0:mpi_size-1)
+     ncell = SUM(cells%ncell)
+     CALL MPI_ALLGATHER(ncell,1,MPI_INTEGER, &
+          ncell_glo,1, MPI_INTEGER, comm_icosa, ierr)
+     displ=0
+     DO i=1,mpi_rank
+        displ=displ+ncell_glo(i-1)
+     ENDDO
+     ncell_tot=sum(ncell_glo(:))
+     
+     ALLOCATE(lon(ncell), lat(ncell), ind_glo(ncell))
+     ALLOCATE(bounds_lon(nvert,ncell), bounds_lat(nvert,ncell))
+     
+     n_beg=0
+     DO ind=1,ndomain
+        n_end = n_beg + cells(ind)%ncell
+        ind_glo(n_beg+1:n_end) = cells(ind)%ind_glo(:)
+        lon(n_beg+1:n_end) = cells(ind)%lon(:)
+        lat(n_beg+1:n_end) = cells(ind)%lat(:)
+        bounds_lon(:,n_beg+1:n_end) = cells(ind)%bnds_lon(:,:)
+        bounds_lat(:,n_beg+1:n_end) = cells(ind)%bnds_lat(:,:)
+        n_beg = n_end
+     END DO
+   END SUBROUTINE collect_bounds
    
  END SUBROUTINE xios_init_write_field
-  
-  
+ 
  SUBROUTINE xios_write_field(name,field)
- USE field_mod
- IMPLICIT NONE 
    CHARACTER(LEN=*),INTENT(IN) :: name
    TYPE(t_field), POINTER :: field(:)
-   CHARACTER(LEN=10) :: str_number
-   INTEGER :: iq
+   TYPE(t_cellset), POINTER :: cells(:)
+   INTEGER :: ncells
 
 !$OMP BARRIER
 !$OMP MASTER
-   
-   IF (Field(1)%field_type==field_T) THEN
-     IF (field(1)%ndim==2) THEN
-        CALL xios_write_field_scalar(name,field,1,1)
-     ELSE IF (field(1)%ndim==3) THEN
-        CALL xios_write_field_scalar(name,field,size(field(1)%rval3d,2),1)
-     ELSE IF (field(1)%ndim==4) THEN
-!        DO iq=1,size(field(1)%rval4d,3)
-!          WRITE(str_number,'(i10)') iq
-!          CALL xios_write_field_scalar(name//TRIM(ADJUSTL(str_number)),field,size(field(1)%rval4d,2),iq)
-          CALL xios_write_field_scalar(name,field,size(field(1)%rval4d,2),size(field(1)%rval4d,3))
-!        ENDDO
-     ELSE
-        PRINT *, "xios_write_field : dimension > 4 are not supported for now"
-     ENDIF
-
-   ELSE IF (Field(1)%field_type==field_U) THEN
-      IF (field(1)%ndim==2) THEN
-        CALL xios_write_field_U(name,field,1,1)
-      ELSE IF (field(1)%ndim==3) THEN
-        CALL xios_write_field_U(name,field,size(field(1)%rval3d,2),1)
-      ELSE IF (field(1)%ndim==4) THEN
-        CALL xios_write_field_U(name,field,size(field(1)%rval4d,2),size(field(1)%rval4d,3))
-      ELSE
-        PRINT *, "xios_write_field : dimension > 4 are not supported for now"
-      ENDIF
-
-    ELSE IF (Field(1)%field_type==field_Z) THEN
-     IF (field(1)%ndim==2) THEN
-        CALL xios_write_field_vort(name,field,1)
-      ELSE IF (field(1)%ndim==3) THEN
-        CALL xios_write_field_vort(name,field,size(field(1)%rval3d,2))
-      ELSE IF (field(1)%ndim==4) THEN
-        DO iq=1,size(field(1)%rval4d,3)
-          WRITE(str_number,'(i10)') iq
-          CALL xios_write_field_vort(name//TRIM(ADJUSTL(str_number)),field,size(field(1)%rval4d,2),iq)
-        ENDDO
-      ELSE
-        PRINT *, "xios_write_field : dimension > 4 are not supported for now"
-      ENDIF
-    ENDIF
+
+   SELECT CASE(field(1)%field_type)
+   CASE(field_T)
+      cells => mesh_loc%primal_own
+      ncells = ncell_i
+   CASE(field_U)
+      cells => mesh_loc%edge_own
+      ncells = ncell_e
+   CASE(field_Z)
+      cells => mesh_loc%dual_own
+      ncells = ncell_v
+   END SELECT
+
+   IF (field(1)%ndim>4) THEN
+      PRINT *, "xios_write_field : dimension > 4 are not supported for now"
+   ELSE
+      CALL xios_write_field_hex(name, field, cells, & 
+           ncells, field(1)%dim3, field(1)%dim4)
+   END IF
+   
 !$OMP END MASTER
 !$OMP BARRIER
@@ -335,52 +171,31 @@
 
  SUBROUTINE xios_read_field(name,field)
- USE field_mod
- IMPLICIT NONE 
    CHARACTER(LEN=*),INTENT(IN) :: name
    TYPE(t_field), POINTER :: field(:)
-   CHARACTER(LEN=10) :: str_number
-   INTEGER :: iq
+   TYPE(t_cellset), POINTER :: cells(:)
+   INTEGER :: ncells
 
 !$OMP BARRIER
 !$OMP MASTER
-   
-   IF (Field(1)%field_type==field_T) THEN
-     IF (field(1)%ndim==2) THEN
-        CALL xios_read_field_scalar(name,field,1,1)
-      ELSE IF (field(1)%ndim==3) THEN
-        CALL xios_read_field_scalar(name,field,size(field(1)%rval3d,2),1)
-      ELSE IF (field(1)%ndim==4) THEN
-!        DO iq=1,size(field(1)%rval4d,3)
-!          WRITE(str_number,'(i10)') iq
-!          CALL xios_read_field_scalar(name//TRIM(ADJUSTL(str_number)),field,size(field(1)%rval4d,2),iq)
-!        ENDDO
-          CALL xios_read_field_scalar(name,field,size(field(1)%rval4d,2),size(field(1)%rval4d,3))
-      ELSE
-        PRINT *, "xios_write_field : dimension > 4 are not supported for now"
-      ENDIF
-   ELSE IF (Field(1)%field_type==field_U) THEN
-     IF (field(1)%ndim==2) THEN
-        CALL xios_read_field_u(name,field,1,1)
-      ELSE IF (field(1)%ndim==3) THEN
-        CALL xios_read_field_u(name,field,size(field(1)%rval3d,2),1)
-      ELSE IF (field(1)%ndim==4) THEN
-          CALL xios_read_field_u(name,field,size(field(1)%rval4d,2),size(field(1)%rval4d,3))
-      ELSE
-        PRINT *, "xios_write_field : dimension > 4 are not supported for now"
-      ENDIF
-    ELSE IF (Field(1)%field_type==field_Z) THEN
-     IF (field(1)%ndim==2) THEN
-        CALL xios_read_field_vort(name,field,1)
-      ELSE IF (field(1)%ndim==3) THEN
-        CALL xios_read_field_vort(name,field,size(field(1)%rval3d,2))
-      ELSE IF (field(1)%ndim==4) THEN
-        DO iq=1,size(field(1)%rval4d,3)
-          WRITE(str_number,'(i10)') iq
-          CALL xios_read_field_vort(name//TRIM(ADJUSTL(str_number)),field,size(field(1)%rval4d,2),iq)
-        ENDDO
-      ELSE
-        PRINT *, "xios_write_field : dimension > 4 are not supported for now"
-      ENDIF
-    ENDIF
+
+   SELECT CASE(field(1)%field_type)
+   CASE(field_T)
+      cells => mesh_loc%primal_own
+      ncells = ncell_i
+   CASE(field_U)
+      cells => mesh_loc%edge_own
+      ncells = ncell_e
+   CASE(field_Z)
+      cells => mesh_loc%dual_own
+      ncells = ncell_v
+   END SELECT
+
+   IF (field(1)%ndim>4) THEN
+      PRINT *, "xios_read_field : dimension > 4 are not supported for now"
+   ELSE
+      CALL xios_read_field_hex(name, field, cells, & 
+           ncells, field(1)%dim3, field(1)%dim4)
+   END IF
+   
 !$OMP END MASTER
 !$OMP BARRIER
@@ -388,78 +203,80 @@
  END SUBROUTINE xios_read_field
 
-
- 
- SUBROUTINE xios_write_field_scalar(name,field,nlev,nq)
- USE genmod
- USE mpipara
- USE xios
- USE grid_param
- USE domain_mod
- USE dimensions
- USE spherical_geom_mod
- USE geometry
- USE mpi_mod
- IMPLICIT NONE 
-   CHARACTER(LEN=*),INTENT(IN) :: name
-   TYPE(t_field), POINTER :: field(:)
-   INTEGER,INTENT(IN) :: nlev
-   INTEGER,INTENT(IN) :: nq
-   
-   REAL(rstd) :: field_tmp(ncell_i,nlev,nq)
-   TYPE(t_domain),POINTER :: d
-   INTEGER :: n,i,j,ij,ind
-   
-   IF (field(1)%ndim==2) THEN
-     n=0
-     DO ind=1,ndomain
-       
-       d=>domain(ind)
-        
-       DO j=d%jj_begin,d%jj_end
-         DO i=d%ii_begin,d%ii_end
-           IF (d%own(i,j)) THEN
-             n=n+1
-             ij=d%iim*(j-1)+i
-             field_tmp(n,1,1)=field(ind)%rval2d(ij)
-           ENDIF
-         ENDDO
-       ENDDO
-     ENDDO
-   ELSE IF (field(1)%ndim==3) THEN
-     n=0
-     DO ind=1,ndomain
-       d=>domain(ind)
-        
-       DO j=d%jj_begin,d%jj_end
-         DO i=d%ii_begin,d%ii_end
-           IF (d%own(i,j)) THEN
-             n=n+1
-             ij=d%iim*(j-1)+i
-             field_tmp(n,:,1)=field(ind)%rval3d(ij,:)
-           ENDIF
-         ENDDO
-       ENDDO
-     ENDDO
-   ELSE IF (field(1)%ndim==4) THEN
-     n=0
-     DO ind=1,ndomain
-       d=>domain(ind)
-        
-       DO j=d%jj_begin,d%jj_end
-         DO i=d%ii_begin,d%ii_end
-           IF (d%own(i,j)) THEN
-             n=n+1
-             ij=d%iim*(j-1)+i
-             field_tmp(n,:,:)=field(ind)%rval4d(ij,:,:)
-           ENDIF
-         ENDDO
-       ENDDO
-     ENDDO     
-   ENDIF
-   
-   CALL xios_send_field(name,field_tmp)
- 
- END SUBROUTINE xios_write_field_scalar  
-
+ SUBROUTINE xios_write_field_hex(name, field, cells, ncell_tot, nlev, nq)
+   CHARACTER(LEN=*),INTENT(IN) :: name
+   TYPE(t_field) :: field(:)
+   TYPE(t_cellset), TARGET :: cells(:)
+   INTEGER,INTENT(IN) :: ncell_tot, nlev, nq
+
+   REAL(rstd) :: field_tmp(ncell_tot,nlev,nq)
+   TYPE(t_cellset), POINTER :: cellset
+   INTEGER :: ind, n_beg, n_end, n, ij, sgn
+   LOGICAL :: signed
+
+   IF(ALLOCATED(cells(1)%sgn)) THEN
+      signed=.TRUE.
+   ELSE
+      signed=.FALSE.
+      sgn=1
+   END IF
+
+   n_beg=0
+   DO ind=1,ndomain
+      cellset => cells(ind)
+      n_end = n_beg + cellset%ncell
+      DO n=1, cellset%ncell
+         ij = cellset%ij(n)
+         IF(signed) sgn = cellset%sgn(n)
+         SELECT CASE(field(1)%ndim)
+         CASE(2)
+            field_tmp(n_beg+n,1,1) = sgn*field(ind)%rval2d(ij)
+         CASE(3)
+            field_tmp(n_beg+n,:,1) = sgn*field(ind)%rval3d(ij,:)
+         CASE(4)
+            field_tmp(n_beg+n,:,:) = sgn*field(ind)%rval4d(ij,:,:)
+         END SELECT
+      END DO
+   END DO
+   CALL xios_send_field(name,field_tmp)  
+ END SUBROUTINE xios_write_field_hex
+
+SUBROUTINE xios_read_field_hex(name, field, cells, ncell_tot, nlev, nq)
+   CHARACTER(LEN=*),INTENT(IN) :: name
+   TYPE(t_field) :: field(:)
+   TYPE(t_cellset), TARGET :: cells(:)
+   INTEGER,INTENT(IN) :: ncell_tot, nlev, nq
+
+   REAL(rstd) :: field_tmp(ncell_tot,nlev,nq)
+   TYPE(t_cellset), POINTER :: cellset
+   INTEGER :: ind, n_beg, n_end, n, ij, sgn
+   LOGICAL :: signed
+
+   CALL xios_recv_field(name,field_tmp)
+
+   IF(ALLOCATED(cells(1)%sgn)) THEN
+      signed=.TRUE.
+   ELSE
+      signed=.FALSE.
+      sgn=1
+   END IF
+
+   n_beg=0
+   DO ind=1,ndomain
+      cellset => cells(ind)
+      n_end = n_beg + cellset%ncell
+      DO n=1, cellset%ncell
+         ij = cellset%ij(n)
+         IF(signed) sgn = cellset%sgn(n)
+         SELECT CASE(field(1)%ndim)
+         CASE(2)
+            field(ind)%rval2d(ij) = sgn*field_tmp(n_beg+n,1,1) 
+         CASE(3)
+            field(ind)%rval3d(ij,:) = sgn*field_tmp(n_beg+n,:,1)
+         CASE(4)
+            field(ind)%rval4d(ij,:,:) = sgn*field_tmp(n_beg+n,:,:)
+         END SELECT
+      END DO
+   END DO
+ END SUBROUTINE xios_read_field_hex
 
  SUBROUTINE xios_read_var(name,field)
@@ -473,469 +290,6 @@
    CALL bcast_omp(field)
  END SUBROUTINE
-
- SUBROUTINE xios_read_field_scalar(name,field,nlev,nq)
- USE genmod
- USE mpipara
- USE xios
- USE grid_param
- USE domain_mod
- USE dimensions
- USE spherical_geom_mod
- USE geometry
- USE mpi_mod
- IMPLICIT NONE 
-   CHARACTER(LEN=*),INTENT(IN) :: name
-   TYPE(t_field), POINTER :: field(:)
-   INTEGER,INTENT(IN) :: nlev
-   INTEGER,INTENT(IN) :: nq
-   
-   REAL(rstd) :: field_tmp(ncell_i,nlev,nq)
-   TYPE(t_domain),POINTER :: d
-   INTEGER :: n,i,j,ij,ind
-
-   CALL xios_recv_field(name,field_tmp)
-   
-   IF (field(1)%ndim==2) THEN
-     n=0
-     DO ind=1,ndomain
-       
-       d=>domain(ind)
-        
-       DO j=d%jj_begin,d%jj_end
-         DO i=d%ii_begin,d%ii_end
-           IF (d%own(i,j)) THEN
-             n=n+1
-             ij=d%iim*(j-1)+i
-             field(ind)%rval2d(ij)=field_tmp(n,1,1)
-           ENDIF
-         ENDDO
-       ENDDO
-     ENDDO
-   ELSE IF (field(1)%ndim==3) THEN
-     n=0
-     DO ind=1,ndomain
-       d=>domain(ind)
-        
-       DO j=d%jj_begin,d%jj_end
-         DO i=d%ii_begin,d%ii_end
-           IF (d%own(i,j)) THEN
-             n=n+1
-             ij=d%iim*(j-1)+i
-             field(ind)%rval3d(ij,:)=field_tmp(n,:,1)
-           ENDIF
-         ENDDO
-       ENDDO
-     ENDDO
-   ELSE IF (field(1)%ndim==4) THEN
-     n=0
-     DO ind=1,ndomain
-       d=>domain(ind)
-        
-       DO j=d%jj_begin,d%jj_end
-         DO i=d%ii_begin,d%ii_end
-           IF (d%own(i,j)) THEN
-             n=n+1
-             ij=d%iim*(j-1)+i
-             field(ind)%rval4d(ij,:,:)=field_tmp(n,:,:)
-           ENDIF
-         ENDDO
-       ENDDO
-     ENDDO     
-   ENDIF
- 
- END SUBROUTINE xios_read_field_scalar
-
- SUBROUTINE xios_write_field_U(name,field,nlev,nq)
- USE genmod
- USE mpipara
- USE xios
- USE grid_param
- USE domain_mod
- USE dimensions
- USE spherical_geom_mod
- USE geometry
- USE mpi_mod
- IMPLICIT NONE 
-   CHARACTER(LEN=*),INTENT(IN) :: name
-   TYPE(t_field), POINTER :: field(:)
-   INTEGER,INTENT(IN) :: nlev
-   INTEGER,INTENT(IN) :: nq
-   
-   REAL(rstd) :: field_tmp(ncell_e,nlev,nq)
-   TYPE(t_domain),POINTER :: d
-   INTEGER :: n,i,j,k,ij,ind
-   
-   IF (field(1)%ndim==2) THEN
-     n=0
-     DO ind=1,ndomain
-       d=>domain(ind)
-       CALL swap_dimensions(ind)
-       CALL swap_geometry(ind)
-
-       DO j=d%jj_begin,d%jj_end
-         DO i=d%ii_begin,d%ii_end
-           DO k=0,5
-             IF (d%edge_assign_domain(k,i,j)==domloc_glo_ind(ind) .AND. d%edge_assign_i(k,i,j)==i .AND. d%edge_assign_j(k,i,j)==j &
-                                                                                 .AND. d%edge_assign_pos(k,i,j)==k) THEN
-               n=n+1
-               ij=iim*(j-1)+i
-               Field_tmp(n,1,1)=d%edge_assign_sign(k,i,j)*field(ind)%rval2d(ij+d%u_pos(k+1))
-             ENDIF
-           ENDDO
-         ENDDO
-       ENDDO
-     ENDDO       
-  
-   ELSE IF (field(1)%ndim==3) THEN
-
-     n=0
-     DO ind=1,ndomain
-       d=>domain(ind)
-       CALL swap_dimensions(ind)
-       CALL swap_geometry(ind)
-
-       DO j=d%jj_begin,d%jj_end
-         DO i=d%ii_begin,d%ii_end
-           DO k=0,5
-             IF (d%edge_assign_domain(k,i,j)==domloc_glo_ind(ind) .AND. d%edge_assign_i(k,i,j)==i .AND. d%edge_assign_j(k,i,j)==j &
-                                                                                  .AND. d%edge_assign_pos(k,i,j)==k) THEN
-               n=n+1
-               ij=iim*(j-1)+i
-               Field_tmp(n,:,1)=d%edge_assign_sign(k,i,j)*field(ind)%rval3d(ij+d%u_pos(k+1),:)
-             ENDIF
-           ENDDO
-         ENDDO
-       ENDDO
-     ENDDO       
-
-   ELSE IF (field(1)%ndim==4) THEN
-
-     n=0
-     DO ind=1,ndomain
-       d=>domain(ind)
-       CALL swap_dimensions(ind)
-       CALL swap_geometry(ind)
-
-       DO j=d%jj_begin,d%jj_end
-         DO i=d%ii_begin,d%ii_end
-           DO k=0,5
-             IF (d%edge_assign_domain(k,i,j)==domloc_glo_ind(ind) .AND. d%edge_assign_i(k,i,j)==i .AND. d%edge_assign_j(k,i,j)==j &
-                                                                                  .AND. d%edge_assign_pos(k,i,j)==k) THEN
-               n=n+1
-               ij=iim*(j-1)+i
-               Field_tmp(n,:,:)=d%edge_assign_sign(k,i,j)*field(ind)%rval4d(ij+d%u_pos(k+1),:,:)
-             ENDIF
-           ENDDO
-         ENDDO
-       ENDDO
-     ENDDO       
-
-   ENDIF
-   
-   CALL xios_send_field(name,field_tmp)
- 
- END SUBROUTINE xios_write_field_u 
-
-
- SUBROUTINE xios_read_field_u(name,field,nlev,nq)
- USE genmod
- USE mpipara
- USE xios
- USE grid_param
- USE domain_mod
- USE dimensions
- USE spherical_geom_mod
- USE geometry
- USE mpi_mod
- IMPLICIT NONE 
-   CHARACTER(LEN=*),INTENT(IN) :: name
-   TYPE(t_field), POINTER :: field(:)
-   INTEGER,INTENT(IN) :: nlev
-   INTEGER,INTENT(IN) :: nq
-   
-   REAL(rstd) :: field_tmp(ncell_e,nlev,nq)
-   TYPE(t_domain),POINTER :: d
-   INTEGER :: n,i,j,k,ij,ind
-
-   CALL xios_recv_field(name,field_tmp)
-   
-   IF (field(1)%ndim==2) THEN
-     n=0
-     DO ind=1,ndomain
-       d=>domain(ind)
-       CALL swap_dimensions(ind)
-       CALL swap_geometry(ind)
-
-       DO j=d%jj_begin,d%jj_end
-         DO i=d%ii_begin,d%ii_end
-           DO k=0,5
-             IF (d%edge_assign_domain(k,i,j)==domloc_glo_ind(ind) .AND. d%edge_assign_i(k,i,j)==i .AND. d%edge_assign_j(k,i,j)==j &
-                                                                                 .AND. d%edge_assign_pos(k,i,j)==k) THEN
-               n=n+1
-               ij=iim*(j-1)+i
-               field(ind)%rval2d(ij+d%u_pos(k+1))=Field_tmp(n,1,1)*d%edge_assign_sign(k,i,j)
-             ENDIF
-           ENDDO
-         ENDDO
-       ENDDO
-     ENDDO       
-  
-   ELSE IF (field(1)%ndim==3) THEN
-
-     n=0
-     DO ind=1,ndomain
-       d=>domain(ind)
-       CALL swap_dimensions(ind)
-       CALL swap_geometry(ind)
-
-       DO j=d%jj_begin,d%jj_end
-         DO i=d%ii_begin,d%ii_end
-           DO k=0,5
-             IF (d%edge_assign_domain(k,i,j)==domloc_glo_ind(ind) .AND. d%edge_assign_i(k,i,j)==i .AND. d%edge_assign_j(k,i,j)==j &
-                                                                                  .AND. d%edge_assign_pos(k,i,j)==k) THEN
-               n=n+1
-               ij=iim*(j-1)+i
-               field(ind)%rval3d(ij+d%u_pos(k+1),:)=Field_tmp(n,:,1)*d%edge_assign_sign(k,i,j)
-             ENDIF
-           ENDDO
-         ENDDO
-       ENDDO
-     ENDDO       
-
-   ELSE IF (field(1)%ndim==4) THEN
-
-     n=0
-     DO ind=1,ndomain
-       d=>domain(ind)
-       CALL swap_dimensions(ind)
-       CALL swap_geometry(ind)
-
-       DO j=d%jj_begin,d%jj_end
-         DO i=d%ii_begin,d%ii_end
-           DO k=0,5
-             IF (d%edge_assign_domain(k,i,j)==domloc_glo_ind(ind) .AND. d%edge_assign_i(k,i,j)==i .AND. d%edge_assign_j(k,i,j)==j &
-                                                                                  .AND. d%edge_assign_pos(k,i,j)==k) THEN
-               n=n+1
-               ij=iim*(j-1)+i
-               field(ind)%rval4d(ij+d%u_pos(k+1),:,:)=Field_tmp(n,:,:)*d%edge_assign_sign(k,i,j)
-             ENDIF
-           ENDDO
-         ENDDO
-       ENDDO
-     ENDDO       
-
-   ENDIF
-   
- 
- END SUBROUTINE xios_read_field_u 
-
-
-
-      
- SUBROUTINE xios_write_field_vort(name,field,nlev,iq)
- USE genmod
- USE mpipara
- USE xios
- USE grid_param
- USE domain_mod
- USE dimensions
- USE spherical_geom_mod
- USE geometry
- USE mpi_mod
- IMPLICIT NONE 
-   CHARACTER(LEN=*),INTENT(IN) :: name
-   TYPE(t_field), POINTER :: field(:)
-   INTEGER,INTENT(IN) :: nlev
-   INTEGER,INTENT(IN),OPTIONAL :: iq
-   
-   REAL(rstd) :: field_tmp(ncell_v,nlev)
-   TYPE(t_domain),POINTER :: d
-   INTEGER :: n,i,j,ij,ind
-   
-   IF (field(1)%ndim==2) THEN
-     n=0
-     DO ind=1,ndomain
-       d=>domain(ind)
-       CALL swap_dimensions(ind)  
-       
-       DO j=d%jj_begin+1,d%jj_end
-         DO i=d%ii_begin,d%ii_end-1
-           n=n+1
-           ij=iim*(j-1)+i
-           Field_tmp(n,1)=field(ind)%rval2d(ij+z_down)
-         ENDDO
-       ENDDO
-
-       DO j=d%jj_begin,d%jj_end-1
-         DO i=d%ii_begin+1,d%ii_end
-           n=n+1
-           ij=iim*(j-1)+i
-           Field_tmp(n,1)=field(ind)%rval2d(ij+z_up)
-          ENDDO
-       ENDDO
-          
-     ENDDO
-
-   ELSE IF (field(1)%ndim==3) THEN
-     n=0
-     DO ind=1,ndomain
-       d=>domain(ind)
-       CALL swap_dimensions(ind)    
-             
-       DO j=d%jj_begin+1,d%jj_end
-         DO i=d%ii_begin,d%ii_end-1
-           n=n+1
-           ij=iim*(j-1)+i
-           Field_tmp(n,:)=field(ind)%rval3d(ij+z_down,:)
-         ENDDO
-       ENDDO
-
-       DO j=d%jj_begin,d%jj_end-1
-         DO i=d%ii_begin+1,d%ii_end
-           n=n+1
-           ij=iim*(j-1)+i
-           Field_tmp(n,:)=field(ind)%rval3d(ij+z_up,:)
-          ENDDO
-       ENDDO
-          
-     ENDDO
-
-   ELSE IF (field(1)%ndim==4) THEN
-     n=0
-     DO ind=1,ndomain
-       d=>domain(ind)
-       CALL swap_dimensions(ind) 
-               
-       DO j=d%jj_begin+1,d%jj_end
-         DO i=d%ii_begin,d%ii_end-1
-           n=n+1
-           ij=iim*(j-1)+i
-           Field_tmp(n,:)=field(ind)%rval4d(ij+z_down,:,iq)
-         ENDDO
-       ENDDO
-
-       DO j=d%jj_begin,d%jj_end-1
-         DO i=d%ii_begin+1,d%ii_end
-           n=n+1
-           ij=iim*(j-1)+i
-           Field_tmp(n,:)=field(ind)%rval4d(ij+z_up,:,iq)
-          ENDDO
-       ENDDO
-          
-     ENDDO
-
-   ENDIF
-   
-   CALL xios_send_field(name,field_tmp)
- 
- END SUBROUTINE xios_write_field_vort 
-
- SUBROUTINE xios_read_field_vort(name,field,nlev,iq)
- USE genmod
- USE mpipara
- USE xios
- USE grid_param
- USE domain_mod
- USE dimensions
- USE spherical_geom_mod
- USE geometry
- USE mpi_mod
- IMPLICIT NONE 
-   CHARACTER(LEN=*),INTENT(IN) :: name
-   TYPE(t_field), POINTER :: field(:)
-   INTEGER,INTENT(IN) :: nlev
-   INTEGER,INTENT(IN),OPTIONAL :: iq
-   
-   REAL(rstd) :: field_tmp(ncell_v,nlev)
-   TYPE(t_domain),POINTER :: d
-   INTEGER :: n,i,j,ij,ind
-
-   CALL xios_recv_field(name,field_tmp)
-
-   
-   IF (field(1)%ndim==2) THEN
-     n=0
-     DO ind=1,ndomain
-       d=>domain(ind)
-       CALL swap_dimensions(ind)  
-       
-       DO j=d%jj_begin+1,d%jj_end
-         DO i=d%ii_begin,d%ii_end-1
-           n=n+1
-           ij=iim*(j-1)+i
-           field(ind)%rval2d(ij+z_down)=Field_tmp(n,1)
-         ENDDO
-       ENDDO
-
-       DO j=d%jj_begin,d%jj_end-1
-         DO i=d%ii_begin+1,d%ii_end
-           n=n+1
-           ij=iim*(j-1)+i
-           Field_tmp(n,1)=field(ind)%rval2d(ij+z_up)
-           field(ind)%rval2d(ij+z_up)=Field_tmp(n,1)
-          ENDDO
-       ENDDO
-          
-     ENDDO
-
-   ELSE IF (field(1)%ndim==3) THEN
-     n=0
-     DO ind=1,ndomain
-       d=>domain(ind)
-       CALL swap_dimensions(ind)    
-             
-       DO j=d%jj_begin+1,d%jj_end
-         DO i=d%ii_begin,d%ii_end-1
-           n=n+1
-           ij=iim*(j-1)+i
-           field(ind)%rval3d(ij+z_down,:)=Field_tmp(n,:)
-         ENDDO
-       ENDDO
-
-       DO j=d%jj_begin,d%jj_end-1
-         DO i=d%ii_begin+1,d%ii_end
-           n=n+1
-           ij=iim*(j-1)+i
-           field(ind)%rval3d(ij+z_up,:)=Field_tmp(n,:)
-          ENDDO
-       ENDDO
-          
-     ENDDO
-
-   ELSE IF (field(1)%ndim==4) THEN
-     n=0
-     DO ind=1,ndomain
-       d=>domain(ind)
-       CALL swap_dimensions(ind) 
-               
-       DO j=d%jj_begin+1,d%jj_end
-         DO i=d%ii_begin,d%ii_end-1
-           n=n+1
-           ij=iim*(j-1)+i
-           field(ind)%rval4d(ij+z_down,:,iq)=Field_tmp(n,:)
-         ENDDO
-       ENDDO
-
-       DO j=d%jj_begin,d%jj_end-1
-         DO i=d%ii_begin+1,d%ii_end
-           n=n+1
-           ij=iim*(j-1)+i
-           field(ind)%rval4d(ij+z_up,:,iq)=Field_tmp(n,:)
-          ENDDO
-       ENDDO
-          
-     ENDDO
-
-   ENDIF
- 
- END SUBROUTINE xios_read_field_vort 
-
-
-
-
   
  SUBROUTINE xios_write_field_finalize
- IMPLICIT NONE
-
 !$OMP BARRIER
 !$OMP MASTER 
@@ -943,9 +297,7 @@
 !$OMP END MASTER
 !$OMP BARRIER
-
  END SUBROUTINE xios_write_field_finalize
 
  SUBROUTINE xios_set_context
- IMPLICIT NONE   
   TYPE(xios_context) :: ctx_hdl
 
@@ -971,12 +323,8 @@
   
   SUBROUTINE xios_init
-   IMPLICIT NONE
-         
-    using_xios=.FALSE.
-    
+    using_xios=.FALSE.    
   END SUBROUTINE xios_init
   
   SUBROUTINE xios_send_field_scalar(name,field)
-  IMPLICIT NONE
     CHARACTER(LEN=*),INTENT(IN) :: name
     REAL,INTENT(IN) :: field
@@ -984,5 +332,4 @@
 
   SUBROUTINE xios_send_field_1d(name,field)
-  IMPLICIT NONE
     CHARACTER(LEN=*),INTENT(IN) :: name
     REAL,INTENT(IN) :: field(:)
@@ -991,5 +338,4 @@
   SUBROUTINE xios_write_field(name,field)
   USE field_mod
-  IMPLICIT NONE 
    CHARACTER(LEN=*),INTENT(IN) :: name
    TYPE(t_field), POINTER :: field(:)
@@ -998,5 +344,4 @@
   SUBROUTINE xios_read_field(name,field)
   USE field_mod
-  IMPLICIT NONE 
    CHARACTER(LEN=*),INTENT(IN) :: name
    TYPE(t_field), POINTER :: field(:)
@@ -1010,5 +355,4 @@
  
   SUBROUTINE xios_update_calendar(step)
-  IMPLICIT NONE
    INTEGER, INTENT(IN):: step  
   END SUBROUTINE xios_update_calendar
Index: codes/icosagcm/devel/src/sphere/compute_geometry.f90
===================================================================
--- codes/icosagcm/devel/src/sphere/compute_geometry.f90	(revision 879)
+++ codes/icosagcm/devel/src/sphere/compute_geometry.f90	(revision 880)
@@ -1,4 +1,18 @@
 MODULE compute_geometry_mod
   USE geometry
+  USE dimensions
+
+  USE domain_mod, ONLY : t_domain, t_cellset, &
+       domain, ndomain, assigned_domain, &
+       domain_glo, ndomain_glo, domloc_glo_ind, swap_needed
+  USE omp_para, ONLY : is_omp_level_master, is_master
+  USE transfert_mod, ONLY : req_i0, req_i1, t_message, transfert_request, transfert_message, init_message
+
+  USE spherical_geom_mod, ONLY : xyz2lonlat, circumcenter, &
+       compute_centroid, centroid, &
+       surf_triangle, dist_cart, div_arc_bis, &
+       schmidt_transform
+  USE vector, ONLY : norm, cross_product2
+
   IMPLICIT NONE
   PRIVATE
@@ -10,12 +24,5 @@
 
   SUBROUTINE update_circumcenters 
-    USE domain_mod
-    USE dimensions
-    USE spherical_geom_mod
-    USE vector
-    USE transfert_mod
-    USE omp_para
-
-    IMPLICIT NONE
+
     REAL(rstd) :: x1(3),x2(3)
     REAL(rstd) :: vect(3,6)
@@ -36,5 +43,5 @@
     
     DO ind=1,ndomain
-      IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
+
       CALL swap_dimensions(ind)
       CALL swap_geometry(ind)
@@ -55,10 +62,5 @@
 
   SUBROUTINE remap_schmidt_loc
-    USE spherical_geom_mod
     USE getin_mod
-    USE omp_para
-    USE domain_mod
-    USE dimensions
-    IMPLICIT NONE
     INTEGER :: ind,i,j,n
     REAL(rstd) :: schmidt_factor, schmidt_lon, schmidt_lat
@@ -91,13 +93,5 @@
 
   SUBROUTINE optimize_geometry
-  USE metric
-  USE spherical_geom_mod
-  USE domain_mod
-  USE dimensions
-  USE transfert_mod
-  USE vector
   USE getin_mod
-  USE omp_para
-  IMPLICIT NONE
     INTEGER :: nb_it=0
     TYPE(t_domain),POINTER :: d
@@ -193,8 +187,5 @@
     ! copy position of generators and vertices back into domain(:)%xyz/vertex
     ! so that XIOS/create_header_gen gets the right values
-    USE omp_para
-    USE dimensions
-    USE domain_mod 
-    USE transfert_mpi_mod
+    USE transfert_mpi_mod, ONLY : gather_field, bcast_field
 
     INTEGER :: ind,i,j,k,n
@@ -253,12 +244,4 @@
   SUBROUTINE set_geometry
   USE metric
-  USE vector
-  USE spherical_geom_mod
-  USE domain_mod
-  USE dimensions
-  USE transfert_mod
-  USE getin_mod
-  USE omp_para
-  IMPLICIT NONE
 
     REAL(rstd) :: surf(6)
@@ -279,7 +262,4 @@
     CALL remap_schmidt_loc
     CALL update_circumcenters
-    ! copy position of generators and vertices back into domain(:)%xyz/vertex
-    ! so that XIOS gets the right values
-    CALL update_domain
 
     DO ind=1,ndomain
@@ -527,7 +507,6 @@
  
   END SUBROUTINE set_geometry
-  
+
   SUBROUTINE compute_wee(n,pos,w)
-  IMPLICIT NONE
     INTEGER,INTENT(IN) :: n
     INTEGER,INTENT(IN) :: pos
@@ -555,12 +534,9 @@
    END SUBROUTINE compute_wee
             
-
   
   SUBROUTINE compute_geometry
     USE grid_param
-    USE domain_mod, ONLY : swap_needed
     USE init_unstructured_mod, ONLY : read_local_mesh
-    IMPLICIT NONE
-
+    USE set_bounds_mod, ONLY : set_bounds
     CALL allocate_geometry
 
@@ -568,4 +544,9 @@
     CASE(grid_ico)
        CALL set_geometry
+       ! copy position of generators and vertices back into domain_glo(:)%xyz/vertex
+       ! so that write_field gets the right values
+       CALL update_domain
+       CALL set_bounds(domain_glo, .TRUE.)
+       CALL set_bounds(domain, .FALSE.)
     CASE(grid_unst)
        swap_needed=.FALSE.
