module write_field_mpi_mod
   USE genmod
   USE write_field_vars_mod
   IMPLICIT NONE
   PRIVATE

CONTAINS

   SUBROUTINE Writefield_mpi(name_in, field, nind)
      USE netcdf_mod
      USE domain_mod
      use field_mod
      USE dimensions
      USE geometry
      IMPLICIT NONE
      CHARACTER(LEN=*), INTENT(IN) :: name_in
      TYPE(t_field), POINTER :: field(:)
      INTEGER, OPTIONAL, INTENT(IN) :: nind
      REAL(r8), ALLOCATABLE :: field_val2d(:)
      REAL(r8), ALLOCATABLE :: field_val3d(:, :)
      TYPE(t_domain), POINTER :: d
      INTEGER :: Index
      INTEGER :: ind, i, j, l, k, n, ncell, q
      INTEGER :: status
      CHARACTER(len=255) :: name
      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.
      END IF

      Index = GetFieldIndex(name)
      if (Index == -1) then
         call create_header_mpi(name, field, nind)
         Index = GetFieldIndex(name)
      else
         FieldIndex(Index) = FieldIndex(Index) + 1.
      end if

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

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

            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)
                     END IF
                  END DO
               END DO
               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, :)
                     END IF
                  END DO
               END DO
               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)
                        END IF
                     END DO
                  END DO

                  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)
               END DO
            END IF

            ncell = ncell + n

         END DO

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

            DO j = jj_begin, jj_end - 1
               DO i = ii_begin + 1, ii_end
                  n = n + 1
               END DO
            END DO

            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)
                  END DO
               END DO

               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)
                  END DO
               END DO

               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, :)
                  END DO
               END DO

               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, :)
                  END DO
               END DO

               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)
                     END DO
                  END DO

                  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)
                     END DO
                  END DO

                  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)
               END DO
            END IF

            ncell = ncell + n

         END DO

      END IF
      status = NF90_SYNC(FieldId(Index))

   END SUBROUTINE Writefield_mpi

   SUBROUTINE Create_header_mpi(name, field, nind)
      USE netcdf_mod
      USE field_mod
      USE domain_mod
      USE spherical_geom_mod
      USE dimensions
      USE geometry
      USE mpi_mod
      USE mpipara
      IMPLICIT NONE
      CHARACTER(LEN=*) :: name
      CHARACTER(LEN=LEN_TRIM(ADJUSTL(name))) :: name_adj
      TYPE(t_field), POINTER :: field(:)
      INTEGER, OPTIONAL, INTENT(IN) :: nind
      INTEGER :: ncell
      INTEGER :: nvert
      REAL(rstd), ALLOCATABLE :: lon(:), lat(:), bounds_lon(:, :), bounds_lat(:, :)
      TYPE(t_domain), POINTER :: d
      INTEGER :: nvertId, ncid, lonId, latId, bounds_lonId, bounds_latId, timeId, ncellId
      INTEGER :: dim3id, dim4id
      INTEGER :: status
      INTEGER :: ind, i, j, k, n, q
      INTEGER :: 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.
      END IF

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

         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)
         END DO
         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)
         END IF

         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/))
            END DO
         END IF

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

            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
                     END DO
                  END IF
               END DO
            END DO
            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
               END DO
            END DO

            DO j = d%jj_begin, d%jj_end - 1
               DO i = d%ii_begin + 1, d%ii_end
                  ncell = ncell + 1
               END DO
            END DO

         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)
         END DO
         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)
         END IF

         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/))
            END DO
         END IF

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

            DO j = jj_begin, jj_end - 1
               DO i = ii_begin + 1, ii_end
                  n = n + 1
               END DO
            END DO

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

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

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

   END SUBROUTINE Create_Header_mpi

end module write_field_mpi_mod
