Index: codes/icosagcm/trunk/src/write_field.f90
===================================================================
--- codes/icosagcm/trunk/src/write_field.f90	(revision 21)
+++ codes/icosagcm/trunk/src/write_field.f90	(revision 26)
@@ -6,4 +6,5 @@
     INTEGER :: size
     INTEGER,POINTER :: nc_id(:)
+    INTEGER :: displ
   END TYPE ncvar
 
@@ -37,46 +38,521 @@
       enddo
     end function GetFieldIndex
- 
-
-      
-    subroutine WriteField_gen(name,Field,dimx,dimy,dimz)
-    implicit none
-!    include 'netcdf.inc'
-      character(len=*) :: name
-      integer :: dimx,dimy,dimz
-      real,dimension(dimx,dimy,dimz) :: Field
-      integer,dimension(dimx*dimy*dimz) :: ndex
-      integer :: status
-      integer :: index
-      integer :: start(4)
-      integer :: count(4)
-      
-           
+    
+    SUBROUTINE Writefield(name_in,field,nind)
+    USE domain_mod
+    USE field_mod
+    USE transfert_mpi_mod
+    USE dimensions
+    USE mpipara
+    IMPLICIT NONE  
+     CHARACTER(LEN=*),INTENT(IN) :: name_in
+      TYPE(t_field),POINTER :: field(:)
+      INTEGER,OPTIONAL,INTENT(IN) :: nind
+      TYPE(t_field),POINTER :: field_glo(:)
+      
+      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 (mpi_rank==0) THEN
+        IF (PRESENT(nind)) THEN
+         CALL writefield_gen(name_in,field_glo,domain_glo,nind)
+        ELSE
+         CALL writefield_gen(name_in,field_glo,domain_glo,1,ndomain_glo)
+        ENDIF
+      ENDIF
+      
+      CALL deallocate_field(field_glo)
+      
+   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 )
+    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(:)
+      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(:,:,:)
+      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(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
+        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
+        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 CreateNewField(name,dimx,dimy,dimz)
-        Index=GetFieldIndex(name)
+        call create_header_gen(name_in,field,domain_type,ind_b_in,ind_e_in)
+	Index=GetFieldIndex(name)
       else
         FieldIndex(Index)=FieldIndex(Index)+1.
       endif
       
-      start(1)=1
-      start(2)=1
-      start(3)=1
-      start(4)=FieldIndex(Index)
-
-      count(1)=dimx
-      count(2)=dimy
-      count(3)=dimz
-      count(4)=1
-
-!      status = NF_PUT_VARA_DOUBLE(FieldId(Index),FieldVarId(Index),start,count,Field)
-!      status = NF_SYNC(FieldId(Index))
-      
-    end subroutine WriteField_gen
-
-
-    SUBROUTINE Writefield(name_in,field,nind)
-    USE netcdf
+      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
+          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
+          status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val2d,  &
+                              start=(/ 1,FieldIndex(Index) /),count=(/ncell,1 /))
+          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
+           
+          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 /))
+          DEALLOCATE(field_val3d)
+        ELSE IF (field(1)%ndim==4) THEN
+
+          DO q=1,FieldVarId(index)%size
+          
+            ALLOCATE(Field_val3d(n,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
+            
+            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 /))
+            DEALLOCATE(field_val3d)
+          ENDDO         
+        ENDIF
+        
+!        PRINT *,NF90_STRERROR(status)
+!        ncell=ncell+n
+
+!     ENDDO
+     
+     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
+          
+          status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),                       &
+                              Field_val2d,start=(/ 1,FieldIndex(Index) /),count=(/ncell,1 /))
+          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
+          
+          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 /))
+          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
+
+            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 /))
+            DEALLOCATE(field_val3d)
+          ENDDO
+        ENDIF
+        
+!        PRINT *,NF90_STRERROR(status)
+!        ncell=ncell+n
+!
+!     ENDDO
+     
+     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
@@ -92,5 +568,5 @@
       TYPE(t_domain),POINTER :: d
       INTEGER :: Index
-      INTEGER :: ind,i,j,k,n,ncell,q
+      INTEGER :: ind,i,j,l,k,n,ncell,q
       INTEGER :: iie,jje,iin,jjn
       INTEGER :: status
@@ -100,4 +576,5 @@
       INTEGER :: halo_size
       LOGICAL :: single
+      INTEGER :: displ
       
       
@@ -120,5 +597,5 @@
       Index=GetFieldIndex(name)
       if (Index==-1) then
-        call create_header(name,field,nind)
+        call create_header_mpi(name,field,nind)
 	Index=GetFieldIndex(name)
       else
@@ -142,4 +619,6 @@
         ENDDO
 
+        displ=FieldVarId(index)%displ
+        
         IF (field(ind)%ndim==2) THEN
           ALLOCATE(Field_val2d(n))
@@ -155,5 +634,5 @@
           ENDDO
           status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val2d,  &
-                              start=(/ ncell,FieldIndex(Index) /),count=(/n,1 /))
+                              start=(/ displ+ncell,FieldIndex(Index) /),count=(/n,1 /))
           DEALLOCATE(field_val2d)
         ELSE IF (field(ind)%ndim==3) THEN
@@ -169,6 +648,8 @@
             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 /))
+!           DO l=1,size(field(ind)%rval3d,2)
+            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 /))
+!           ENDDO
           DEALLOCATE(field_val3d)
         ELSE IF (field(1)%ndim==4) THEN
@@ -186,7 +667,9 @@
                 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 /))
+             ENDDO
+!            DO l=1,size(field(ind)%rval4d,2)
+            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 /))
+!            ENDDO
             DEALLOCATE(field_val3d)
           ENDDO         
@@ -219,4 +702,6 @@
           ENDDO
 
+        displ=FieldVarId(index)%displ
+
         IF (field(ind)%ndim==2) THEN
           ALLOCATE(Field_val2d(n))
@@ -240,5 +725,5 @@
 
           status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),                       &
-                              Field_val2d,start=(/ ncell,FieldIndex(Index) /),count=(/n,1 /))
+                              Field_val2d,start=(/ displ+ncell,FieldIndex(Index) /),count=(/n,1 /))
           DEALLOCATE(field_val2d)
 
@@ -261,6 +746,8 @@
             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 /))
+!           DO l=1,size(field(ind)%rval3d,2)
+           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 /))
+!           ENDDO
           DEALLOCATE(field_val3d)
         ELSE IF (field(1)%ndim==4) THEN
@@ -285,6 +772,9 @@
             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 /))
+!            DO l=1,size(field(ind)%rval4d,2)
+
+            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 /))
+!            ENDDO
             DEALLOCATE(field_val3d)
           ENDDO
@@ -299,9 +789,630 @@
      status=NF90_SYNC(FieldId(Index))
       
-    END SUBROUTINE Writefield
-      
+    END SUBROUTINE Writefield_mpi
+    
+    
            
-    SUBROUTINE Create_header(name,field,nind)
-    USE netcdf
+!   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
+
+
+
+   
+    SUBROUTINE Create_header_gen(name_in,field,domain_type,ind_b_in,ind_e_in)
+    USE netcdf_mod
+    USE field_mod
+    USE domain_mod
+    USE metric
+    USE spherical_geom_mod
+!    USE dimensions
+!    USE geometry
+    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
+      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
+      
+            
+      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
+        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
+        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      
+                
+      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
+        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_type(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%assign_domain(i,j)==ind) n=n+1
+!            ENDDO
+!          ENDDO
+          
+         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 /))
+  
+!         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_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
+        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_type(ind)
+!  
+!          n=0
+!          DO j=d%jj_begin+1,d%jj_end
+!            DO i=d%ii_begin,d%ii_end-1
+!              n=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
+!            ENDDO
+!          ENDDO
+
+!         ALLOCATE(lon(n),lat(n),bounds_lon(0:nvert-1,n),bounds_lat(0:nvert-1,n))
+         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
+ !              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))
+
+               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
+!                IF (bounds_lon(k,n)<0) bounds_lon(k,n)=bounds_lon(k,n)+360
+               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(xyz_v(nij+z_up,:)/radius,lon(n),lat(n))
+               CALL xyz2lonlat(d%vertex(:,vup,i,j),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))
+               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
+!                 IF (bounds_lon(k,n)<0) bounds_lon(k,n)=bounds_lon(k,n)+360
+               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 /))
+  
+!          ncell=ncell+n
+          DEALLOCATE(lon,lat,bounds_lon,bounds_lat)
+!      END DO          
+    ENDIF
+
+
+      
+!      status = NF90_CLOSE(ncid)
+
+    END SUBROUTINE Create_Header_gen
+
+    SUBROUTINE Create_header_mpi(name,field,nind)
+    USE netcdf_mod
     USE field_mod
     USE domain_mod
@@ -309,4 +1420,6 @@
     USE dimensions
     USE geometry
+    USE mpi_mod
+    USE mpipara
     IMPLICIT NONE
       CHARACTER(LEN=*) :: name
@@ -326,4 +1439,7 @@
       LOGICAL :: single 
       INTEGER :: nij
+      INTEGER :: ncell_glo(0:mpi_size-1)
+      INTEGER :: displ, ncell_tot
+      
           
       NbField=NbField+1
@@ -359,7 +1475,16 @@
         END DO
       
-        status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid)
+        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,ncellId)
+        status = NF90_DEF_DIM(ncid,'cell',ncell_tot,ncellId)
         status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid)
 
@@ -394,7 +1519,9 @@
           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")
+          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)),NF90_DOUBLE,(/ 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
@@ -402,4 +1529,5 @@
                                   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 
@@ -437,8 +1565,8 @@
             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 /))
+          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
@@ -464,8 +1592,18 @@
 
         END DO
-      
-        status = NF90_CREATE(TRIM(ADJUSTL(name))//'.nc', NF90_CLOBBER, ncid)
+        
+        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,ncellId)
+        status = NF90_DEF_DIM(ncid,'cell',ncell_tot,ncellId)
         status = NF90_DEF_DIM(ncid,'nvert',nvert,nvertid)
 
@@ -502,7 +1640,9 @@
           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")
+          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)),NF90_DOUBLE,(/ 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
@@ -510,4 +1650,5 @@
                                   (/ 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 
@@ -579,8 +1720,8 @@
           
           
-          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 /))
+          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
@@ -593,7 +1734,17 @@
 !      status = NF90_CLOSE(ncid)
 
-    end subroutine Create_Header
+    end subroutine Create_Header_mpi  
    
-  
+   SUBROUTINE Close_files
+   USE netcdf
+   IMPLICIT NONE
+     INTEGER :: i,k,status
+     
+     DO i=1,NbField
+         status=NF90_CLOSE(FieldId(i))
+     ENDDO
+   END SUBROUTINE  Close_files
+     
+     
   function int2str(int)
     implicit none
