module write_field_mod
USE genmod
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
  
  CONTAINS
  
    SUBROUTINE init_writeField
    USE ioipsl
    use netcdf_mod
    IMPLICIT NONE
      CHARACTER(LEN=255) :: netcdf_prec
      
      netcdf_prec='float'
      CALL getin("netcdf_prec",netcdf_prec)
      
      SELECT CASE(TRIM(netcdf_prec))
        CASE('float')
          ncprec=NF90_FLOAT
        CASE('double')
          ncprec=NF90_DOUBLE
        CASE default
        PRINT*,'Bad selector for variable netcdf_prec : <', TRIM(netcdf_prec),"> options are <float>, <double>" 
        STOP   
      END SELECT
    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)
    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
      LOGICAL,OPTIONAL,INTENT(IN) :: once
      LOGICAL                     :: once_
      
      TYPE(t_field),POINTER :: field_glo(:)

!$OMP BARRIER
!$OMP MASTER      
      IF(PRESENT(once)) THEN
         once_=once
      ELSE
         once_=.FALSE.
      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 (mpi_rank==0) THEN
        IF (PRESENT(nind)) THEN
         CALL writefield_gen(name_in,field_glo,domain_glo,nind,once=once_)
        ELSE
         CALL writefield_gen(name_in,field_glo,domain_glo,1,ndomain_glo,once=once_)
        ENDIF
      ENDIF
      
      CALL deallocate_field_glo(field_glo)
!$OMP END MASTER
!$OMP BARRIER
      
   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))

      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 create_header_gen(name_in,field,domain_type,ind_b_in,ind_e_in,once)
	Index=GetFieldIndex(name)
      else
        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
          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

          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-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



   
    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
      
            
      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_i',ncell,ncellId)
        status = NF90_DEF_DIM(ncid,'nvert_i',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_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_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)      
         
         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
        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_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 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  
   
   SUBROUTINE Close_files
   USE netcdf
   IMPLICIT NONE
     INTEGER :: i,k,status
!$OMP MASTER     
     DO i=1,NbField
         status=NF90_CLOSE(FieldId(i))
     ENDDO
!$OMP END MASTER
   END SUBROUTINE  Close_files
     
     
  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