write_diachro.f90

  if ( Any( tpfields(jp)%ndimlist(:) == NMNHDIM_UNKNOWN ) ) &
    call Print_msg( NVERB_ERROR, 'IO', 'Write_diachro_nc4', &
                    'some dimensions are unknown for variable '//trim(tpfields(jp)%cmnhname) )

  icount = Count( tpfields(jp)%ndimlist(:) /= NMNHDIM_UNUSED )

  if ( tpfields(jp)%ndims /= icount ) &
    call Print_msg( NVERB_ERROR, 'IO', 'Write_diachro_nc4', &
                    'ndims is not coherent with ndimlist for variable '//trim(tpfields(jp)%cmnhname) )

  if ( jp == 1 ) then
    idims = icount
  else
    if ( idims /= icount ) &
      call Print_msg( NVERB_ERROR, 'IO', 'Write_diachro_nc4', &
                      'number of dimensions not the same for all tpfields for variable '//trim(tpfields(jp)%cmnhname) )
  end if
end do

!The dimension list should be the same for all tpfields entries
do jp = 2, Size( tpfields )
  do ji = 1, NMNHMAXDIMS
    if ( tpfields(jp)%ndimlist(ji) /= tpfields(1)%ndimlist(ji) ) then
      !For SERIES: it is possible to have NMNHDIM_NI and NMNHDIM_NI_U in the different tpfields
      !For SERIES: it is possible to have NMNHDIM_SERIES_LEVEL and NMNHDIM_SERIES_LEVEL_W in the different tpfields
      if ( tpfields(jp)%ndimlist(ji) /= NMNHDIM_NI           .and. tpfields(jp)%ndimlist(ji) /= NMNHDIM_NI_U .and.     &
           tpfields(jp)%ndimlist(ji) /= NMNHDIM_SERIES_LEVEL .and. tpfields(jp)%ndimlist(ji) /= NMNHDIM_SERIES_LEVEL_W ) then
        call Print_msg( NVERB_ERROR, 'IO', 'Write_diachro_nc4', &
                        'some dimensions are not the same for all tpfields entries for variable '//trim(tpfields(jp)%cmnhname) )
      end if
    end if
  end do
end do

!Check that if cartesian and no horizontal compression, parameters are as expected
if ( yshape == 'Cartesian' .and. .not. tpbudiachro%licompress .and. .not. tpbudiachro%ljcompress ) then
  icorr = Merge( 1, 0, tpbudiachro%lkcompress )
  if ( ( idims + icorr ) /= 3 .and. ( idims + icorr ) /= 4 ) then
    call Print_msg( NVERB_ERROR, 'IO', 'Write_diachro_nc4',                                                            &
                    'unexpected number of dimensions for cartesian shape without horizontal compression for variable ' &
                    // Trim( tpfields(1)%cmnhname ) )
  end if

  if (      (       tpfields(1)%ndimlist(1) /= NMNHDIM_BUDGET_CART_NI              &
              .and. tpfields(1)%ndimlist(1) /= NMNHDIM_BUDGET_CART_NI_U            &
              .and. tpfields(1)%ndimlist(1) /= NMNHDIM_BUDGET_CART_NI_V          ) &
       .or. (       tpfields(1)%ndimlist(2) /= NMNHDIM_BUDGET_CART_NJ              &
              .and. tpfields(1)%ndimlist(2) /= NMNHDIM_BUDGET_CART_NJ_U            &
              .and. tpfields(1)%ndimlist(2) /= NMNHDIM_BUDGET_CART_NJ_V          ) &
       .or. (       .not. tpbudiachro%lkcompress                                                     &
              .and. tpfields(1)%ndimlist(3) /= NMNHDIM_BUDGET_CART_LEVEL           &
              .and. tpfields(1)%ndimlist(3) /= NMNHDIM_BUDGET_CART_LEVEL_W       ) &
       .or. ( ( idims + icorr ) == 4 .and. tpfields(1)%ndimlist(6) /= NMNHDIM_BUDGET_NGROUPS ) ) then
    call Print_msg( NVERB_ERROR, 'IO', 'Write_diachro_nc4',                                       &
                    'unexpected dimensions for CART without horizontal compression for variable ' &
                    // Trim( tpfields(1)%cmnhname ) )
  end if
end if


select case ( idims )
  case (0)
     !Remark: [ integer:: ] is a constructor for a zero-size array of integers, [] is not allowed (type can not be determined)
      call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(1), pvar, [ integer:: ], gsplit, gdistributed )

  case (1)

    if ( Any ( tpfields(1)%ndimlist(4) == [ NMNHDIM_BUDGET_LES_TIME, NMNHDIM_BUDGET_LES_AVG_TIME, NMNHDIM_SERIES_TIME ] ) ) then
      if ( Size( tpfields ) /= 1 ) call Print_msg( NVERB_FATAL, 'IO', 'Write_diachro_nc4', &
                                                   'wrong size of tpfields (variable '//trim(tpfields(1)%cmnhname)//')' )

      call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(1), pvar, [ 4 ], gsplit, gdistributed )
    else if ( Any( tpfields(1)%ndimlist(1) == [ NMNHDIM_NI, NMNHDIM_NI_U, NMNHDIM_NI_V, NMNHDIM_BUDGET_CART_NI,     &
                                                NMNHDIM_BUDGET_CART_NI_U, NMNHDIM_BUDGET_CART_NI_V              ] ) ) then
      if ( Size( tpfields ) /= 1 ) call Print_msg( NVERB_FATAL, 'IO', 'Write_diachro_nc4', &
                                                   'wrong size of tpfields (variable '//trim(tpfields(1)%cmnhname)//')' )

      call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(1), pvar, [ 1 ], gsplit, gdistributed )

    else if ( Any( tpfields(1)%ndimlist(2) == [ NMNHDIM_NJ, NMNHDIM_NJ_U, NMNHDIM_NJ_V, NMNHDIM_BUDGET_CART_NJ,     &
                                                NMNHDIM_BUDGET_CART_NJ_U, NMNHDIM_BUDGET_CART_NJ_V              ] ) ) then
      if ( Size( tpfields ) /= 1 ) call Print_msg( NVERB_FATAL, 'IO', 'Write_diachro_nc4', &
                                                   'wrong size of tpfields (variable '//trim(tpfields(1)%cmnhname)//')' )
      call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(1), pvar, [ 2 ], gsplit, gdistributed )
    else if ( Any( tpfields(1)%ndimlist(3) == [ NMNHDIM_LEVEL, NMNHDIM_LEVEL_W,                            &
                                                NMNHDIM_BUDGET_CART_LEVEL, NMNHDIM_BUDGET_CART_LEVEL_W ] ) ) then
      if ( Size( tpfields ) /= 1 ) call Print_msg( NVERB_FATAL, 'IO', 'Write_diachro_nc4', &
                                                   'wrong size of tpfields (variable '//trim(tpfields(1)%cmnhname)//')' )
      call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(1), pvar, [ 3 ], gsplit, gdistributed )
    else if ( tpfields(1)%ndimlist(6) == NMNHDIM_BUDGET_NGROUPS .or. tpfields(1)%ndimlist(6) == NMNHDIM_PROFILER_PROC ) then
      do ji = 1, Size( pvar, 6 )
        !Remark: [ integer:: ] is a constructor for a zero-size array of integers, [] is not allowed (type can not be determined)
        call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(ji), pvar(:,:,:,:,:,ji:ji), [ integer:: ], &
                                          gsplit, gdistributed )
      end do
    else
      call Print_msg( NVERB_FATAL, 'IO', 'Write_diachro_nc4', &
                      'case not yet implemented (1D variable '//trim(tpfields(1)%cmnhname)//')' )
    end if


  case (2)

    if (       Any( tpfields(1)%ndimlist(1) == [ NMNHDIM_NI, NMNHDIM_NI_U, NMNHDIM_NI_V, NMNHDIM_BUDGET_CART_NI, &
                                                 NMNHDIM_BUDGET_CART_NI_U, NMNHDIM_BUDGET_CART_NI_V ] )          &
         .and. Any( tpfields(1)%ndimlist(2) == [ NMNHDIM_NJ, NMNHDIM_NJ_U, NMNHDIM_NJ_V, NMNHDIM_BUDGET_CART_NJ, &
                                                 NMNHDIM_BUDGET_CART_NJ_U, NMNHDIM_BUDGET_CART_NJ_V ] )          ) then
      if ( Size( tpfields ) /= 1 ) call Print_msg( NVERB_FATAL, 'IO', 'Write_diachro_nc4', &
                                                   'wrong size of tpfields (variable '//trim(tpfields(1)%cmnhname)//')' )
      call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(1), pvar, [ 1, 2 ], gsplit, gdistributed, &
                                        iil, iih, ijl, ijh, ikl, ikh )
    else if (       Any( tpfields(1)%ndimlist(1) == [ NMNHDIM_NI, NMNHDIM_NI_U, NMNHDIM_NI_V, NMNHDIM_BUDGET_CART_NI, &
                                                      NMNHDIM_BUDGET_CART_NI_U, NMNHDIM_BUDGET_CART_NI_V ] )          &
              .and. Any( tpfields(1)%ndimlist(3) == [ NMNHDIM_LEVEL, NMNHDIM_LEVEL_W,                                 &
                                                      NMNHDIM_BUDGET_CART_LEVEL, NMNHDIM_BUDGET_CART_LEVEL_W ] )      ) then
      if ( Size( tpfields ) /= 1 ) call Print_msg( NVERB_FATAL, 'IO', 'Write_diachro_nc4', &
                                                   'wrong size of tpfields (variable '//trim(tpfields(1)%cmnhname)//')' )
      call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(1), pvar, [ 1, 3 ], gsplit, gdistributed )
    else if (       Any( tpfields(1)%ndimlist(2) == [ NMNHDIM_NJ, NMNHDIM_NJ_U, NMNHDIM_NJ_V, NMNHDIM_BUDGET_CART_NJ, &
                                                      NMNHDIM_BUDGET_CART_NJ_U, NMNHDIM_BUDGET_CART_NJ_V ] )          &
              .and. Any( tpfields(1)%ndimlist(3) == [ NMNHDIM_LEVEL, NMNHDIM_LEVEL_W,                                 &
                                                      NMNHDIM_BUDGET_CART_LEVEL, NMNHDIM_BUDGET_CART_LEVEL_W ] )      ) then
      if ( Size( tpfields ) /= 1 ) call Print_msg( NVERB_FATAL, 'IO', 'Write_diachro_nc4', &
                                                   'wrong size of tpfields (variable '//trim(tpfields(1)%cmnhname)//')' )
      call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(1), pvar, [ 2, 3 ], gsplit, gdistributed )
    else if (  Any( tpfields(1)%ndimlist(1) == [ NMNHDIM_BUDGET_CART_NI, NMNHDIM_BUDGET_CART_NI_U, NMNHDIM_BUDGET_CART_NI_V ] ) &
              .and. tpfields(1)%ndimlist(6) == NMNHDIM_BUDGET_NGROUPS ) then
      ! Loop on the processes
      do ji = 1, Size( pvar, 6 )
        call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(ji), pvar(:,:,:,:,:,ji:ji), [ 1 ], gsplit, gdistributed )
      end do
    else if (  Any( tpfields(1)%ndimlist(2) == [ NMNHDIM_BUDGET_CART_NJ, NMNHDIM_BUDGET_CART_NJ_U, NMNHDIM_BUDGET_CART_NJ_V ] ) &
              .and. tpfields(1)%ndimlist(6) == NMNHDIM_BUDGET_NGROUPS ) then
      ! Loop on the processes
      do ji = 1, Size( pvar, 6 )
        call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(ji), pvar(:,:,:,:,:,ji:ji), [ 2 ], gsplit, gdistributed )
      end do
    else if (       tpfields(1)%ndimlist(3) == NMNHDIM_BUDGET_LES_LEVEL     &
       .and. (      tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_TIME &
               .or. tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_AVG_TIME ) ) then
      if ( Size( tpfields ) /= 1 ) call Print_msg( NVERB_FATAL, 'IO', 'Write_diachro_nc4', &
                                                   'wrong size of tpfields (variable '//trim(tpfields(1)%cmnhname)//')' )
      call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(1), pvar, [ 3, 4 ], gsplit, gdistributed )
    else if ( Any( tpfields(1)%ndimlist(3) == [ NMNHDIM_BUDGET_CART_LEVEL, NMNHDIM_BUDGET_CART_LEVEL_W ] ) &
              .and. tpfields(1)%ndimlist(6) == NMNHDIM_BUDGET_NGROUPS ) then
      ! Loop on the processes
      do ji = 1, Size( pvar, 6 )
        call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(ji), pvar(:,:,:,:,:,ji:ji), [ 3 ], gsplit, gdistributed )
      end do
    else if (  tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_TIME .and. tpfields(1)%ndimlist(5) == NMNHDIM_BUDGET_MASK_NBUMASK ) then
      if ( Size( tpfields ) /= 1 ) call Print_msg( NVERB_FATAL, 'IO', 'Write_diachro_nc4', &
                                                   'wrong size of tpfields (variable '//trim(tpfields(1)%cmnhname)//')' )
      call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(1), pvar, [ 4, 5 ], gsplit, gdistributed )
    else if (  (      tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_TIME &
                 .or. tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_AVG_TIME ) &
         .and. tpfields(1)%ndimlist(5) == NMNHDIM_BUDGET_LES_SV       ) then
      if ( Size( tpfields ) /= 1 ) call Print_msg( NVERB_FATAL, 'IO', 'Write_diachro_nc4', &
                                                   'wrong size of tpfields (variable '//trim(tpfields(1)%cmnhname)//')' )
      call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(1), pvar, [ 4, 5 ], gsplit, gdistributed )
    else if (  tpfields(1)%ndimlist(4) == NMNHDIM_FLYER_TIME &
         .and. tpfields(1)%ndimlist(6) == NMNHDIM_FLYER_PROC ) then
      !Correspond to FLYER_DIACHRO
      !Create local time dimension
      if ( isp == tzfile%nmaster_rank) then
        istatus = NF90_INQ_DIMID( ilevelids(NLVL_GROUP), 'time_flyer', idimid )
        if ( istatus == NF90_NOERR ) then
          !Dimension already exists, check that it is not changed
          istatus = NF90_INQUIRE_DIMENSION( ilevelids(NLVL_GROUP), idimid, len = ilen )
          if ( ilen /= Int( Size( pvar, 4), kind = CDFINT ) ) &
            call Print_msg( NVERB_FATAL, 'IO', 'Write_diachro_nc4', 'time_flyer dimension has changed' )
        else
          !Dimension does not exist yet, create it
          istatus = NF90_DEF_DIM( ilevelids(NLVL_GROUP), 'time_flyer', Int( Size( pvar, 4), kind = CDFINT ), idimid )
          if ( istatus /= NF90_NOERR ) &
            call IO_Err_handle_nc4( istatus, 'Write_diachro_nc4', 'NF90_DEF_DIM', Trim( tpfields(1)%cmnhname ) )
        end if
      end if

      ! Loop on the processes
      do ji = 1, Size( pvar, 6 )
        call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(ji), pvar(:,:,:,:,:,ji:ji), [ 4 ], gsplit, gdistributed )
      end do
    else if (  tpfields(1)%ndimlist(4) == NMNHDIM_PROFILER_TIME &
         .and. tpfields(1)%ndimlist(6) == NMNHDIM_PROFILER_PROC ) then
      !Correspond to WRITE_PROFILER_n
      ! Loop on the processes
      do ji = 1, Size( pvar, 6 )
        call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(ji), pvar(:,:,:,:,:,ji:ji), [ 4 ], gsplit, gdistributed )
      end do
    else if (  tpfields(1)%ndimlist(4) == NMNHDIM_SERIES_TIME &
         .and. tpfields(1)%ndimlist(6) == NMNHDIM_SERIES_PROC ) then
      !Correspond to WRITE_SERIES_n
      ! Loop on the processes
      do ji = 1, Size( pvar, 6 )
        call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(ji), pvar(:,:,:,:,:,ji:ji), [ 4 ], gsplit, gdistributed )
      end do
    else if (  ( tpfields(1)%ndimlist(3) == NMNHDIM_SERIES_LEVEL .or. tpfields(1)%ndimlist(3) == NMNHDIM_SERIES_LEVEL_W ) &
         .and. tpfields(1)%ndimlist(4) == NMNHDIM_SERIES_TIME ) then
      !Correspond to WRITE_SERIES_n
      if ( Size( tpfields ) /= 1 ) call Print_msg( NVERB_FATAL, 'IO', 'Write_diachro_nc4', &
                                                   'wrong size of tpfields (variable '//trim(tpfields(1)%cmnhname)//')' )
      call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(1), pvar(:,:,:,:,:,:), [ 3, 4 ], gsplit, gdistributed )
    else if (  tpfields(1)%ndimlist(4) == NMNHDIM_STATION_TIME &
         .and. tpfields(1)%ndimlist(6) == NMNHDIM_STATION_PROC ) then
      !Correspond to WRITE_STATION_n
      ! Loop on the processes
      do ji = 1, Size( pvar, 6 )
        call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(ji), pvar(:,:,:,:,:,ji:ji), [ 4 ], gsplit, gdistributed )
      end do
    else
      call Print_msg( NVERB_FATAL, 'IO', 'Write_diachro_nc4', &
                      'case not yet implemented (2D variable '//trim(tpfields(1)%cmnhname)//')' )
    end if


  case (3)

    if (       Any( tpfields(1)%ndimlist(1) == [ NMNHDIM_NI, NMNHDIM_NI_U, NMNHDIM_NI_V, NMNHDIM_BUDGET_CART_NI, &
                                                 NMNHDIM_BUDGET_CART_NI_U, NMNHDIM_BUDGET_CART_NI_V ] )          &
         .and. Any( tpfields(1)%ndimlist(2) == [ NMNHDIM_NJ, NMNHDIM_NJ_U, NMNHDIM_NJ_V, NMNHDIM_BUDGET_CART_NJ, &
                                                 NMNHDIM_BUDGET_CART_NJ_U, NMNHDIM_BUDGET_CART_NJ_V ] )          &
         .and. Any( tpfields(1)%ndimlist(3) == [ NMNHDIM_LEVEL, NMNHDIM_LEVEL_W,                                 &
                                                 NMNHDIM_BUDGET_CART_LEVEL, NMNHDIM_BUDGET_CART_LEVEL_W ] )      ) then
      if ( Size( tpfields ) /= 1 ) call Print_msg( NVERB_FATAL, 'IO', 'Write_diachro_nc4', &
                                                   'wrong size of tpfields (variable '//trim(tpfields(1)%cmnhname)//')' )
      call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(1), pvar, [ 1, 2, 3 ], gsplit, gdistributed, &
                                        iil, iih, ijl, ijh, ikl, ikh )
    else if (       Any(tpfields(1)%ndimlist(1) == [ NMNHDIM_BUDGET_CART_NI, NMNHDIM_BUDGET_CART_NI_U, NMNHDIM_BUDGET_CART_NI_V ]) &
              .and. Any(tpfields(1)%ndimlist(2) == [ NMNHDIM_BUDGET_CART_NJ, NMNHDIM_BUDGET_CART_NJ_U, NMNHDIM_BUDGET_CART_NJ_V ]) &
              .and.     tpfields(1)%ndimlist(6) ==   NMNHDIM_BUDGET_NGROUPS                                                   ) then
      ! Loop on the processes
      do ji = 1, Size( pvar, 6 )
        call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(ji), pvar(:,:,:,:,:,ji:ji), [ 1, 2 ], &
                                          gsplit, gdistributed, iil, iih, ijl, ijh, ikl, ikh )
      end do
    else if (       Any ( tpfields(1)%ndimlist(1) == [ NMNHDIM_BUDGET_CART_NI, NMNHDIM_BUDGET_CART_NI_U,          &
                                                       NMNHDIM_BUDGET_CART_NI_V ] )                               &
              .and. Any ( tpfields(1)%ndimlist(3) == [ NMNHDIM_BUDGET_CART_LEVEL, NMNHDIM_BUDGET_CART_LEVEL_W ] ) &
              .and.       tpfields(1)%ndimlist(6) ==   NMNHDIM_BUDGET_NGROUPS                                     ) then
      ! Loop on the processes
      do ji = 1, Size( pvar, 6 )
        call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(ji), pvar(:,:,:,:,:,ji:ji), [ 1, 3 ], gsplit, gdistributed )
      end do
    else if (       Any ( tpfields(1)%ndimlist(2) == [ NMNHDIM_BUDGET_CART_NJ, NMNHDIM_BUDGET_CART_NJ_U,          &
                                                       NMNHDIM_BUDGET_CART_NJ_V ] )                               &
              .and. Any ( tpfields(1)%ndimlist(3) == [ NMNHDIM_BUDGET_CART_LEVEL, NMNHDIM_BUDGET_CART_LEVEL_W ] ) &
              .and.       tpfields(1)%ndimlist(6) ==   NMNHDIM_BUDGET_NGROUPS                                     ) then
      ! Loop on the processes
      do ji = 1, Size( pvar, 6 )
        call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(ji), pvar(:,:,:,:,:,ji:ji), [ 2, 3 ], gsplit, gdistributed )
      end do
    else if (         (      tpfields(1)%ndimlist(3) == NMNHDIM_BUDGET_MASK_LEVEL     &
                        .or. tpfields(1)%ndimlist(3) == NMNHDIM_BUDGET_MASK_LEVEL_W ) &
               .and.         tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_TIME           &
               .and.         tpfields(1)%ndimlist(5) == NMNHDIM_BUDGET_MASK_NBUMASK   ) then
      !Correspond to Store_one_budget_rho (MASK)
      if ( Size( tpfields ) /= 1 ) call Print_msg( NVERB_FATAL, 'IO', 'Write_diachro_nc4', &
                                                   'wrong size of tpfields (variable '//trim(tpfields(1)%cmnhname)//')' )
      call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(1), pvar, [ 3, 4, 5 ], gsplit, gdistributed, &
                                        iil, iih, ijl, ijh, ikl, ikh )
    else if (              tpfields(1)%ndimlist(3) == NMNHDIM_BUDGET_LES_LEVEL      &
              .and. (      tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_TIME       &
                      .or. tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_AVG_TIME ) &
              .and.        tpfields(1)%ndimlist(6) == NMNHDIM_BUDGET_LES_MASK       ) then
      if ( nles_masks /= Size( pvar, 6 ) ) &
        call Print_msg( NVERB_FATAL, 'IO', 'Write_diachro_nc4',                             &
                        'last dimension size of pvar is not equal to nles_masks (variable ' &
                        // Trim( tpfields(1)%cmnhname ) // ')' )

      ! Loop on the processes
      do ji = 1, Size( pvar, 6 )
        call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(ji), pvar(:,:,:,:,:,ji:ji), [ 3, 4 ], gsplit, gdistributed )
      end do
    else if (       tpfields(1)%ndimlist(3) == NMNHDIM_BUDGET_LES_LEVEL     &
       .and. (      tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_TIME &
               .or. tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_AVG_TIME ) &
       .and. tpfields(1)%ndimlist(6) == NMNHDIM_BUDGET_TERM      ) then
      ! Loop on the processes
      do ji = 1, Size( pvar, 6 )
        call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(ji), pvar(:,:,:,:,:,ji:ji), [ 3, 4 ], gsplit, gdistributed )
      end do
    else if (       tpfields(1)%ndimlist(3) == NMNHDIM_BUDGET_LES_LEVEL      &
       .and. (      tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_TIME       &
               .or. tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_AVG_TIME ) &
              .and. tpfields(1)%ndimlist(5) == NMNHDIM_BUDGET_LES_PDF      ) then
      if ( Size( tpfields ) /= 1 ) call Print_msg( NVERB_FATAL, 'IO', 'Write_diachro_nc4', &
                                                   'wrong size of tpfields (variable '//trim(tpfields(1)%cmnhname)//')' )
      call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(1), pvar, [ 3, 4, 5 ], gsplit, gdistributed )
    else if (              tpfields(1)%ndimlist(3) == NMNHDIM_BUDGET_LES_LEVEL      &
              .and. (      tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_TIME       &
                      .or. tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_AVG_TIME ) &
              .and.        tpfields(1)%ndimlist(5) == NMNHDIM_BUDGET_LES_SV         ) then
      if ( Size( tpfields ) /= 1 ) call Print_msg( NVERB_FATAL, 'IO', 'Write_diachro_nc4', &
                                                   'wrong size of tpfields (variable '//trim(tpfields(1)%cmnhname)//')' )
      call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(1), pvar, [ 3, 4, 5 ], gsplit, gdistributed, &
                                        iil, iih, ijl, ijh, ikl, ikh )
    else if (              tpfields(1)%ndimlist(1) == NMNHDIM_SPECTRA_2PTS_NI       &
              .and.        tpfields(1)%ndimlist(3) == NMNHDIM_SPECTRA_LEVEL         &
              .and. (      tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_TIME       &
                      .or. tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_AVG_TIME ) ) then
      if ( Size( tpfields ) /= 1 ) call Print_msg( NVERB_FATAL, 'IO', 'Write_diachro_nc4', &
                                                   'wrong size of tpfields (variable '//trim(tpfields(1)%cmnhname)//')' )
      call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(1), pvar, [ 1, 3, 4 ], gsplit, gdistributed, &
                                        iil, iih, ijl, ijh, ikl, ikh )
    else if (       tpfields(1)%ndimlist(2) == NMNHDIM_SPECTRA_2PTS_NJ                   &
              .and. tpfields(1)%ndimlist(3) == NMNHDIM_SPECTRA_LEVEL                &
              .and. (      tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_TIME       &
                      .or. tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_AVG_TIME ) ) then
      if ( Size( tpfields ) /= 1 ) call Print_msg( NVERB_FATAL, 'IO', 'Write_diachro_nc4', &
                                                   'wrong size of tpfields (variable '//trim(tpfields(1)%cmnhname)//')' )
      call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(1), pvar, [ 2, 3, 4 ], gsplit, gdistributed, &
                                        iil, iih, ijl, ijh, ikl, ikh )
    else if ( ( tpfields(1)%ndimlist(3) == NMNHDIM_LEVEL .or. tpfields(1)%ndimlist(3) == NMNHDIM_LEVEL_W ) &
         .and. tpfields(1)%ndimlist(4) == NMNHDIM_FLYER_TIME &
         .and. tpfields(1)%ndimlist(6) == NMNHDIM_FLYER_PROC ) then
      !Correspond to FLYER_DIACHRO
      !Create local time dimension
      if ( isp == tzfile%nmaster_rank) then
        istatus = NF90_INQ_DIMID( ilevelids(NLVL_GROUP), 'time_flyer', idimid )
        if ( istatus == NF90_NOERR ) then
          !Dimension already exists, check that it is not changed
          istatus = NF90_INQUIRE_DIMENSION( ilevelids(NLVL_GROUP), idimid, len = ilen )
          if ( ilen /= Int( Size( pvar, 4), kind = CDFINT ) ) &
            call Print_msg( NVERB_FATAL, 'IO', 'Write_diachro_nc4', 'time_flyer dimension has changed' )
        else
          !Dimension does not exist yet, create it
          istatus = NF90_DEF_DIM( ilevelids(NLVL_GROUP), 'time_flyer', Int( Size( pvar, 4), kind = CDFINT ), idimid )
          if ( istatus /= NF90_NOERR ) &
            call IO_Err_handle_nc4( istatus, 'Write_diachro_nc4', 'NF90_DEF_DIM', Trim( tpfields(1)%cmnhname ) )
        end if
      end if

      ! Loop on the processes
      do ji = 1, Size( pvar, 6 )
        call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(ji), pvar(:,:,:,:,:,ji:ji), [ 3, 4 ], gsplit, gdistributed )
      end do
    else if (  ( tpfields(1)%ndimlist(3) == NMNHDIM_LEVEL .or. tpfields(1)%ndimlist(3) == NMNHDIM_LEVEL_W ) &
         .and. tpfields(1)%ndimlist(4) == NMNHDIM_PROFILER_TIME &
         .and. tpfields(1)%ndimlist(6) == NMNHDIM_PROFILER_PROC ) then
      !Correspond to PROFILER_DIACHRO_n
      ! Loop on the processes
      do ji = 1, Size( pvar, 6 )
        call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(ji), pvar(:,:,:,:,:,ji:ji), [ 3, 4 ], gsplit, gdistributed )
      end do
    else if (  ( tpfields(1)%ndimlist(3) == NMNHDIM_SERIES_LEVEL .or. tpfields(1)%ndimlist(3) == NMNHDIM_SERIES_LEVEL_W ) &
         .and. tpfields(1)%ndimlist(4) == NMNHDIM_SERIES_TIME &
         .and. tpfields(1)%ndimlist(6) == NMNHDIM_SERIES_PROC ) then
      !Correspond to PROFILER_DIACHRO_n
      ! Loop on the processes
      do ji = 1, Size( pvar, 6 )
        call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(ji), pvar(:,:,:,:,:,ji:ji), [ 3, 4 ], gsplit, gdistributed )
      end do
    else if (  ( tpfields(1)%ndimlist(1) == NMNHDIM_NI .or. tpfields(1)%ndimlist(1) == NMNHDIM_NI_U )      &
         .and. tpfields(1)%ndimlist(4) == NMNHDIM_SERIES_TIME &
         .and. tpfields(1)%ndimlist(6) == NMNHDIM_SERIES_PROC ) then
      !Correspond to PROFILER_DIACHRO_n
      ! Loop on the processes
      do ji = 1, Size( pvar, 6 )
        call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(ji), pvar(:,:,:,:,:,ji:ji), [ 1, 4 ], gsplit, gdistributed )
      end do
    else if (  ( tpfields(1)%ndimlist(2) == NMNHDIM_NJ .or. tpfields(1)%ndimlist(2) == NMNHDIM_NJ_U )      &
         .and. tpfields(1)%ndimlist(4) == NMNHDIM_SERIES_TIME &
         .and. tpfields(1)%ndimlist(6) == NMNHDIM_SERIES_PROC ) then
      ! Loop on the processes
      do ji = 1, Size( pvar, 6 )
        call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(ji), pvar(:,:,:,:,:,ji:ji), [ 2, 4 ], gsplit, gdistributed )
      end do
    else if (       tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_TIME         &
              .and. tpfields(1)%ndimlist(5) == NMNHDIM_BUDGET_MASK_NBUMASK &
              .and. tpfields(1)%ndimlist(6) == NMNHDIM_BUDGET_NGROUPS      ) then
      ! Loop on the processes
      do ji = 1, Size( pvar, 6 )
        call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(ji), pvar(:,:,:,:,:,ji:ji), [ 4, 5 ], gsplit, gdistributed )
      end do
    else
      call Print_msg( NVERB_FATAL, 'IO', 'Write_diachro_nc4', &
                      'case not yet implemented (3D variable '//trim(tpfields(1)%cmnhname)//')' )
    end if

  case (4)

    if (       Any( tpfields(1)%ndimlist(1) == [ NMNHDIM_BUDGET_CART_NI, NMNHDIM_BUDGET_CART_NI_U, NMNHDIM_BUDGET_CART_NI_V ] ) &
         .and. Any( tpfields(1)%ndimlist(2) == [ NMNHDIM_BUDGET_CART_NJ, NMNHDIM_BUDGET_CART_NJ_U, NMNHDIM_BUDGET_CART_NJ_V ] ) &
         .and. Any( tpfields(1)%ndimlist(3) == [ NMNHDIM_BUDGET_CART_LEVEL,NMNHDIM_BUDGET_CART_LEVEL_W ] )                      &
         .and.      tpfields(1)%ndimlist(6) == NMNHDIM_BUDGET_NGROUPS                                                        ) then
      !Correspond to Store_one_budget (CART)
      ! Loop on the processes
      do ji = 1, Size( pvar, 6 )
        call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(ji), pvar(:,:,:,:,:,ji:ji), [ 1, 2, 3 ], &
                                          gsplit, gdistributed, iil, iih, ijl, ijh, ikl, ikh )
      end do
    elseif (  (        tpfields(1)%ndimlist(3) == NMNHDIM_BUDGET_MASK_LEVEL     &
                  .or. tpfields(1)%ndimlist(3) == NMNHDIM_BUDGET_MASK_LEVEL_W ) &
       .and.           tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_TIME           &
       .and.           tpfields(1)%ndimlist(5) == NMNHDIM_BUDGET_MASK_NBUMASK   &
       .and.           tpfields(1)%ndimlist(6) == NMNHDIM_BUDGET_NGROUPS        ) then
      !Correspond to Store_one_budget (MASK)
      ! Loop on the processes
      do ji = 1, Size( pvar, 6 )
        call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(ji), pvar(:,:,:,:,:,ji:ji), [ 3, 4, 5 ], &
                                          gsplit, gdistributed, iil, iih, ijl, ijh, ikl, ikh )
      end do
    else if (              tpfields(1)%ndimlist(3) == NMNHDIM_BUDGET_LES_LEVEL      &
              .and. (      tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_TIME       &
                      .or. tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_AVG_TIME ) &
       .and.               tpfields(1)%ndimlist(5) == NMNHDIM_BUDGET_LES_SV         &
       .and.               tpfields(1)%ndimlist(6) == NMNHDIM_BUDGET_LES_MASK       ) then
      if ( nles_masks /= Size( pvar, 6 ) ) &
        call Print_msg( NVERB_FATAL, 'IO', 'Write_diachro_nc4',                             &
                        'last dimension size of pvar is not equal to nles_masks (variable ' &
                        // Trim( tpfields(1)%cmnhname ) // ')' )

      ! Loop on the processes
      do ji = 1, Size( pvar, 6 )
        call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(ji), pvar(:,:,:,:,:,ji:ji), [ 3, 4, 5 ], &
                                          gsplit, gdistributed, iil, iih, ijl, ijh, ikl, ikh )
      end do
    else if (              tpfields(1)%ndimlist(3) == NMNHDIM_BUDGET_LES_LEVEL      &
              .and. (      tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_TIME       &
                      .or. tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_AVG_TIME ) &
       .and.               tpfields(1)%ndimlist(5) == NMNHDIM_BUDGET_LES_SV         &
       .and.               tpfields(1)%ndimlist(6) == NMNHDIM_BUDGET_TERM           ) then
      ! Loop on the processes
      do ji = 1, Size( pvar, 6 )
        call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(ji), pvar(:,:,:,:,:,ji:ji), [ 3, 4, 5 ], &
                                          gsplit, gdistributed, iil, iih, ijl, ijh, ikl, ikh )
      end do
    else if (             tpfields(1)%ndimlist(1) == NMNHDIM_SPECTRA_SPEC_NI       &
             .and.        tpfields(1)%ndimlist(3) == NMNHDIM_SPECTRA_LEVEL         &
             .and. (      tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_TIME       &
                     .or. tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_AVG_TIME ) &
             .and.        tpfields(1)%ndimlist(5) == NMNHDIM_COMPLEX               ) then
      !Correspond to LES_DIACHRO_SPEC
      if ( Size( tpfields ) /= 1 ) call Print_msg( NVERB_FATAL, 'IO', 'Write_diachro_nc4', &
                                                   'wrong size of tpfields (variable '//trim(tpfields(1)%cmnhname)//')' )
      call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(1), pvar, [ 1, 3, 4, 5 ], gsplit, gdistributed )
    else if (              tpfields(1)%ndimlist(2) == NMNHDIM_SPECTRA_SPEC_NJ       &
              .and.        tpfields(1)%ndimlist(3) == NMNHDIM_SPECTRA_LEVEL         &
              .and. (      tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_TIME       &
                      .or. tpfields(1)%ndimlist(4) == NMNHDIM_BUDGET_LES_AVG_TIME ) &
              .and.        tpfields(1)%ndimlist(5) == NMNHDIM_COMPLEX               ) then
      !Correspond to LES_DIACHRO_SPEC
      if ( Size( tpfields ) /= 1 ) call Print_msg( NVERB_FATAL, 'IO', 'Write_diachro_nc4', &
                                                   'wrong size of tpfields (variable '//trim(tpfields(1)%cmnhname)//')' )
      call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(1), pvar, [ 2, 3, 4, 5 ], gsplit, gdistributed )
    else
      call Print_msg( NVERB_FATAL, 'IO', 'Write_diachro_nc4', &
                      'case not yet implemented (4D variable '//trim(tpfields(1)%cmnhname)//')' )
    end if

!   case (5)

!   case (6)

  case default
    do ji = 1, Size( pvar, 6 )
      call Diachro_one_field_write_nc4( tzfile, tpbudiachro, tpfields(ji), pvar(:,:,:,:,:,ji:ji), [ 1, 2, 3, 4, 5 ], &
                                        gsplit, gdistributed )
    end do

end select

!Write X and Y position of the flyer
if ( Present( tpflyer ) ) then
  if ( lcartesian ) then
    ystdnameprefix = 'plane'
  else
    ystdnameprefix = 'projection'
  endif

  tzfield%cmnhname   = 'X'
  tzfield%cstdname   = Trim( ystdnameprefix ) // '_x_coordinate'
  tzfield%clongname  = 'x-position of the flyer'
  tzfield%cunits     = 'm'
  tzfield%cdir       = '--'
  tzfield%ccomment   = ''
  tzfield%ngrid      = 0
  tzfield%ntype      = TYPEREAL
  tzfield%ltimedep   = .false.
  tzfield%ndims      = 1
  tzfield%ndimlist(1)  = NMNHDIM_FLYER_TIME
  tzfield%ndimlist(2:) = NMNHDIM_UNUSED

  call IO_Field_write( tzfile, tzfield, tpflyer%x )

  tzfield%cmnhname   = 'Y'
  tzfield%cstdname   = Trim( ystdnameprefix ) // '_y_coordinate'
  tzfield%clongname  = 'y-position of the flyer'

  call IO_Field_write( tzfile, tzfield, tpflyer%y )
end if

end  subroutine Write_diachro_nc4


subroutine Diachro_one_field_write_nc4( tpfile, tpbudiachro, tpfield, pvar, kdims, osplit, odistributed, &
                                        kil, kih, kjl, kjh, kkl, kkh )
use modd_budget,      only: NLVL_CATEGORY, NLVL_GROUP, NLVL_SHAPE, nbutshift, nbusubwrite, tbudiachrometadata
use modd_field,       only: tfielddata, tfield_metadata_base
use modd_io,          only: isp, tfiledata
use modd_parameters,  only: jphext

use mode_io_field_write, only: IO_Field_create, IO_Field_write, IO_Field_write_box

type(tfiledata),                                     intent(in)  :: tpfile        !File to write
type(tbudiachrometadata),                            intent(in)  :: tpbudiachro
class(tfield_metadata_base),                         intent(in)  :: tpfield
real,                        dimension(:,:,:,:,:,:), intent(in)  :: pvar
integer, dimension(:),                               intent(in)  :: kdims        !List of indices of dimensions to use
logical,                                             intent(in)  :: osplit
logical,                                             intent(in)  :: odistributed !.T. if data is distributed among all processes
integer,                                             intent(in), optional :: kil, kih
integer,                                             intent(in), optional :: kjl, kjh
integer,                                             intent(in), optional :: kkl, kkh

integer                                                    :: idims
integer                                                    :: ibutimepos
integer                                                    :: ji
integer,          dimension(size(shape(pvar)))             :: isizes_alldims
integer,          dimension(:),                allocatable :: ioffset
integer,          dimension(:),                allocatable :: isizes
real                                                       :: zdata0d
real,             dimension(:),                allocatable :: zdata1d
real,             dimension(:,:),              allocatable :: zdata2d
real,             dimension(:,:,:),            allocatable :: zdata3d
real,             dimension(:,:,:,:),          allocatable :: zdata4d
real,             dimension(:,:,:,:,:),        allocatable :: zdata5d
type(tfielddata)                                           :: tzfield

idims = Size( kdims )

if ( odistributed ) then
  if ( idims /= 2 .and. idims /= 3 )                                                                 &
    call Print_msg( NVERB_FATAL, 'IO', 'Diachro_one_field_write_nc4',                                &
                   'odistributed=.true. not allowed for dims/=3, field: ' //Trim( tzfield%cmnhname ) )

  if ( tpbudiachro%clevels(NLVL_SHAPE) /= 'Cartesian' )                                                                    &
    call Print_msg( NVERB_FATAL, 'IO', 'Diachro_one_field_write_nc4',                                         &
                   'odistributed=.true. not allowed for shape/=cartesian, field: ' //Trim( tzfield%cmnhname ) )
end if

if ( osplit ) then
  if ( idims > 3 )                                                                                          &
    call Print_msg( NVERB_FATAL, 'IO', 'Diachro_one_field_write_nc4',                                       &
                                 'osplit=.true. not allowed for dims>3, field: ' //Trim( tzfield%cmnhname ) )

  if ( tpbudiachro%clevels(NLVL_CATEGORY) /= 'Budgets' )                                                  &
    call Print_msg( NVERB_FATAL, 'IO', 'Diachro_one_field_write_nc4',                                    &
                    'osplit=.true. not allowed for category/=budget, field: ' //Trim( tzfield%cmnhname ) )
end if

Allocate( isizes(idims) )
isizes_alldims = 1
do ji = 1, idims
  isizes(ji) = Size( pvar, kdims(ji) )
  isizes_alldims(kdims(ji)) = isizes(ji)
end do

call Prepare_diachro_write( tpfield, tzfield, kdims, osplit, odistributed, ibutimepos )

NDIMS: select case( idims )
  case ( 0 ) NDIMS
    zdata0d = pvar(1, 1, 1, 1, 1, 1)

    if ( osplit ) then
      !Create the metadata of the field (has to be done only once)
      if ( nbutshift == 1 ) call IO_Field_create( tpfile, tzfield )

      call IO_Field_write( tpfile, tzfield, [ zdata0d ], koffset= [ ( nbutshift - 1 ) * nbusubwrite ] )
    else
      call IO_Field_write( tpfile, tzfield, zdata0d )
    end if


  case ( 1 ) NDIMS
    ! Copy selected dimensions into zdata (+ auto-allocate it)
    zdata1d = Reshape ( pvar(1:isizes_alldims(1), 1:isizes_alldims(2), 1:isizes_alldims(3),  &
                             1:isizes_alldims(4), 1:isizes_alldims(5), 1:isizes_alldims(6)), &
                             isizes(1:1) )

    if ( osplit ) then
      !Create the metadata of the field (has to be done only once)
      if ( nbutshift == 1 ) call IO_Field_create( tpfile, tzfield )

      Allocate( ioffset( tzfield%ndims ) )
      ioffset(:) = 0
      ioffset(ibutimepos) = ( nbutshift - 1 ) * nbusubwrite

      if ( tzfield%ndims == idims ) then
        !No time dimension was added in Prepare_diachro_write
        call IO_Field_write( tpfile, tzfield, zdata1d(:), koffset = ioffset )
      else if ( tzfield%ndims == ( idims + 1 ) ) then
        !A time dimension was added in Prepare_diachro_write
        call IO_Field_write( tpfile, tzfield, Reshape( zdata1d, [ Size(zdata1d,1), 1 ] ), &
                              koffset = ioffset )
      else
        call Print_msg( NVERB_FATAL, 'IO', 'Diachro_one_field_write_nc4', &
                                         'probable bug for ' //Trim( tzfield%cmnhname ) )
      end if
    else !.not. osplit
      call IO_Field_write( tpfile, tzfield, zdata1d )
    end if


  case ( 2 ) NDIMS
    ! Copy selected dimensions into zdata (+ auto-allocate it)
    zdata2d = Reshape ( pvar(1:isizes_alldims(1), 1:isizes_alldims(2), 1:isizes_alldims(3),  &
                             1:isizes_alldims(4), 1:isizes_alldims(5), 1:isizes_alldims(6)), &
                             isizes(1:2) )

    if ( osplit ) then
      !Create the metadata of the field (has to be done only once)
      if ( nbutshift == 1 ) call IO_Field_create( tpfile, tzfield )

      Allocate( ioffset( tzfield%ndims ) )
      ioffset(:) = 0
      ioffset(ibutimepos) = ( nbutshift - 1 ) * nbusubwrite

      if ( odistributed ) then
        if ( tzfield%ndims == idims ) then
          !No time dimension was added in Prepare_diachro_write
          call IO_Field_write_box( tpfile, tzfield, 'BUDGET',                              &
                                   zdata2d,                                                &
                                   kil + jphext, kih + jphext, kjl + jphext, kjh + jphext, &
                                   koffset = ioffset                                       )
        else if ( tzfield%ndims == ( idims + 1 ) ) then
          !A time dimension was added in Prepare_diachro_write
          call IO_Field_write_box( tpfile, tzfield, 'BUDGET',                                   &
                                   Reshape( zdata2d, [ Size(zdata2d,1), Size(zdata2d,2), 1 ] ), &
                                   kil + jphext, kih + jphext, kjl + jphext, kjh + jphext,      &
                                   koffset = ioffset                                            )
        else
          call Print_msg( NVERB_FATAL, 'IO', 'Diachro_one_field_write_nc4', &
                                             'probable bug for ' //Trim( tzfield%cmnhname ) )
        end if
      else
        !Data is already collected on the master process
        if ( tzfield%ndims == idims ) then
          !No time dimension was added in Prepare_diachro_write
          call IO_Field_write( tpfile, tzfield, zdata2d(:,:), koffset = ioffset )
        else if ( tzfield%ndims == ( idims + 1 ) ) then
          !A time dimension was added in Prepare_diachro_write
          call IO_Field_write( tpfile, tzfield, Reshape( zdata2d, [ Size(zdata2d,1), Size(zdata2d,2), 1 ] ), &
                                koffset = ioffset )
        else
          call Print_msg( NVERB_FATAL, 'IO', 'Diachro_one_field_write_nc4', &
                                             'probable bug for ' //Trim( tzfield%cmnhname ) )
        end if
      end if
    else !.not. osplit
      if ( odistributed ) then
        call IO_Field_write_box( tpfile, tzfield, 'BUDGET', zdata2d, &
                                 kil + jphext, kih + jphext, kjl + jphext, kjh + jphext )
      else
        !Data is already collected on the master process
        call IO_Field_write( tpfile, tzfield, zdata2d )
      end if
    end if


  case ( 3 ) NDIMS
    ! Copy selected dimensions into zdata (+ auto-allocate it)
    zdata3d = Reshape ( pvar(1:isizes_alldims(1), 1:isizes_alldims(2), 1:isizes_alldims(3),  &
                             1:isizes_alldims(4), 1:isizes_alldims(5), 1:isizes_alldims(6)), &
                             isizes(1:3) )

    if ( osplit ) then
      !Create the metadata of the field (has to be done only once)
      if ( nbutshift == 1 ) call IO_Field_create( tpfile, tzfield )

      Allocate( ioffset( tzfield%ndims ) )
      ioffset(:) = 0
      ioffset(ibutimepos) = ( nbutshift - 1 ) * nbusubwrite

      if ( odistributed ) then
        if ( tzfield%ndims == idims ) then
          !No time dimension was added in Prepare_diachro_write
          call IO_Field_write_box( tpfile, tzfield, 'BUDGET',                              &
                                   zdata3d,                                                &
                                   kil + jphext, kih + jphext, kjl + jphext, kjh + jphext, &
                                   koffset = ioffset                                       )
        else if ( tzfield%ndims == ( idims + 1 ) ) then
          !A time dimension was added in Prepare_diachro_write
          call IO_Field_write_box( tpfile, tzfield, 'BUDGET',                                                    &
                                   Reshape( zdata3d, [ Size(zdata3d,1), Size(zdata3d,2), Size(zdata3d,3), 1 ] ), &
                                   kil + jphext, kih + jphext, kjl + jphext, kjh + jphext,                       &
                                   koffset = ioffset                                                             )
        else
          call Print_msg( NVERB_FATAL, 'IO', 'Diachro_one_field_write_nc4', &
                                             'probable bug for ' //Trim( tzfield%cmnhname ) )
        end if
      else
        !Data is already collected on the master process
        if ( tzfield%ndims == idims ) then
          !No time dimension was added in Prepare_diachro_write
          call IO_Field_write( tpfile, tzfield, zdata3d(:,:,:), koffset = ioffset )
        else if ( tzfield%ndims == ( idims + 1 ) ) then
          !A time dimension was added in Prepare_diachro_write
          call IO_Field_write( tpfile, tzfield, Reshape( zdata3d, [ Size(zdata3d,1), Size(zdata3d,2), Size(zdata3d,3), 1 ] ), &
                                koffset = ioffset )
        else
          call Print_msg( NVERB_FATAL, 'IO', 'Diachro_one_field_write_nc4', &
                                             'probable bug for ' //Trim( tzfield%cmnhname ) )
        end if
      end if
    else !.not. osplit
      if ( odistributed ) then
        call IO_Field_write_box( tpfile, tzfield, 'BUDGET', zdata3d, &
                                 kil + jphext, kih + jphext, kjl + jphext, kjh + jphext )
      else
        !Data is already collected on the master process
        call IO_Field_write( tpfile, tzfield, zdata3d )
      end if
    end if


  case ( 4 ) NDIMS
    ! Copy selected dimensions into zdata (+ auto-allocate it)
    zdata4d = Reshape ( pvar(1:isizes_alldims(1), 1:isizes_alldims(2), 1:isizes_alldims(3),  &
                             1:isizes_alldims(4), 1:isizes_alldims(5), 1:isizes_alldims(6)), &
                             isizes(1:4) )

    call IO_Field_write( tpfile, tzfield, zdata4d )


  case ( 5 ) NDIMS
    ! Copy selected dimensions into zdata (+ auto-allocate it)
    zdata5d = Reshape ( pvar(1:isizes_alldims(1), 1:isizes_alldims(2), 1:isizes_alldims(3),  &
                             1:isizes_alldims(4), 1:isizes_alldims(5), 1:isizes_alldims(6)), &
                             isizes(1:5) )

    call IO_Field_write( tpfile, tzfield, zdata5d )


  case default NDIMS
    call Print_msg( NVERB_ERROR, 'IO', 'Diachro_one_field_write_nc4', Trim( tpfile%cname ) // &
                          ': unsupported number of dimensions' )
    return

end select NDIMS

end subroutine Diachro_one_field_write_nc4


subroutine Prepare_diachro_write( tpfieldin, tpfieldout, kdims, osplit, odistributed, kbutimepos )
use modd_field, only: NMNHDIM_BUDGET_TIME, NMNHDIM_UNUSED, NMNHMAXDIMS, tfielddata, tfield_metadata_base

class(tfield_metadata_base), intent(in)  :: tpfieldin
type(tfielddata),            intent(out) :: tpfieldout
integer, dimension(:),       intent(in)  :: kdims ! List of indices of dimensions to use
logical,                     intent(in)  :: osplit
logical,                     intent(in)  :: odistributed ! .true. if data is distributed among all the processes
integer,                     intent(out) :: kbutimepos

integer :: idims
integer :: jdim

idims = Size( kdims )

if ( idims > NMNHMAXDIMS ) call Print_msg( NVERB_FATAL, 'IO', 'Prepare_diachro_write', &
                                           'kdims is too big for ' //Trim( tpfieldin%cmnhname ) )

tpfieldout%cmnhname   = tpfieldin%cmnhname
tpfieldout%cstdname   = tpfieldin%cstdname
tpfieldout%clongname  = tpfieldin%clongname
tpfieldout%cunits     = tpfieldin%cunits
if ( .not. odistributed ) then
  tpfieldout%cdir       = '--'
else
  tpfieldout%cdir       = 'XY'
end if
tpfieldout%ccomment   = tpfieldin%ccomment
tpfieldout%ngrid      = tpfieldin%ngrid
tpfieldout%ntype      = tpfieldin%ntype
tpfieldout%ltimedep   = .false.

tpfieldout%ndims      = idims

do jdim = 1, idims
  tpfieldout%ndimlist(jdim)  = tpfieldin%ndimlist(kdims(jdim))
end do
tpfieldout%ndimlist(idims + 1:) = NMNHDIM_UNUSED

kbutimepos = -1

!Add budget time dimension if required
if ( osplit ) then
  do jdim = 1, idims
    if ( tpfieldout%ndimlist(jdim) == NMNHDIM_BUDGET_TIME ) then
      kbutimepos = jdim
      exit
    end if
  end do

  !budget time dimension was not found => add it
  if ( kbutimepos == -1 ) then
    idims = idims + 1
    if ( idims > NMNHMAXDIMS ) call Print_msg( NVERB_FATAL, 'IO', 'Prepare_diachro_write',  &
                                               'impossible to add NMNHDIM_BUDGET_TIME dimension for ' //Trim( tpfieldin%cmnhname ) )
    kbutimepos = idims
    tpfieldout%ndims = idims
    tpfieldout%ndimlist(idims) = NMNHDIM_BUDGET_TIME
  end if
end if

end subroutine Prepare_diachro_write


subroutine Att_write_c0( hlevel, kgrpid, hattname, hdata )
use NETCDF,            only: NF90_GET_ATT, NF90_INQUIRE_ATTRIBUTE, NF90_PUT_ATT, NF90_CHAR, NF90_GLOBAL, NF90_NOERR

use modd_precision,    only: CDFINT

use mode_io_tools_nc4, only: IO_Err_handle_nc4, IO_Mnhname_clean

character(len=*),     intent(in) :: hlevel
integer(kind=CDFINT), intent(in) :: kgrpid
character(len=*),     intent(in) :: hattname
character(len=*),     intent(in) :: hdata

character(len=Len(hattname))  :: yattname
character(len=:), allocatable :: yatt
integer(kind=CDFINT)          :: ilen
integer(kind=CDFINT)          :: istatus
integer(kind=CDFINT)          :: itype

call IO_Mnhname_clean( hattname, yattname )

istatus = NF90_INQUIRE_ATTRIBUTE( kgrpid, NF90_GLOBAL, yattname, xtype = itype, len = ilen )
if (istatus == NF90_NOERR ) then
  call Print_msg( NVERB_DEBUG, 'IO', 'Write_diachro_nc4', 'attribute ' // yattname // ' already exists for ' // Trim( hlevel ) )

  if ( itype /= NF90_CHAR ) then
    call Print_msg( NVERB_ERROR, 'IO', 'Write_diachro_nc4', 'type for attribute ' // yattname // &
                    ' has changed for ' // Trim( hlevel ) )
    return
  end if

  Allocate( character(len=ilen) :: yatt )
  istatus = NF90_GET_ATT( kgrpid, NF90_GLOBAL, yattname, yatt )
  if ( yatt == Trim( hdata ) ) then
    call Print_msg( NVERB_DEBUG, 'IO', 'Write_diachro_nc4', 'attribute ' // yattname // ' is unchanged for ' // Trim( hlevel ) )
    !If unchanged, no need to write it again => return
    return
  else
    cmnhmsg(1) = 'attribute ' // yattname // ' has changed for ' // Trim( hlevel )
    cmnhmsg(2) = yatt // ' -> ' // Trim( hdata )
    call Print_msg( NVERB_WARNING, 'IO', 'Write_diachro_nc4' )
  end if

end if

istatus = NF90_PUT_ATT( kgrpid, NF90_GLOBAL, yattname, Trim( hdata ) )
if (istatus /= NF90_NOERR ) &
 call IO_Err_handle_nc4( istatus, 'Write_diachro_nc4', 'NF90_PUT_ATT', Trim( yattname ) // ' for '// Trim( hlevel ) // ' group' )

end subroutine Att_write_c0


subroutine Att_write_i0( hlevel, kgrpid, hattname, kdata )
use NETCDF,            only: NF90_GET_ATT, NF90_INQUIRE_ATTRIBUTE, NF90_PUT_ATT, NF90_GLOBAL, NF90_NOERR

use modd_precision,    only: CDFINT, MNHINT_NF90

use mode_io_tools_nc4, only: IO_Err_handle_nc4, IO_Mnhname_clean

character(len=*),     intent(in) :: hlevel
integer(kind=CDFINT), intent(in) :: kgrpid
character(len=*),     intent(in) :: hattname
integer,              intent(in) :: kdata

character(len=Len(hattname)) :: yattname
integer              :: iatt
integer(kind=CDFINT) :: ilen
integer(kind=CDFINT) :: istatus
integer(kind=CDFINT) :: itype

call IO_Mnhname_clean( hattname, yattname )

istatus = NF90_INQUIRE_ATTRIBUTE( kgrpid, NF90_GLOBAL, yattname, xtype = itype, len = ilen )
if (istatus == NF90_NOERR ) then
  call Print_msg( NVERB_DEBUG, 'IO', 'Write_diachro_nc4', 'attribute ' // yattname // ' already exists for ' // Trim( hlevel ) )

  if ( itype /= MNHINT_NF90 ) then
    call Print_msg( NVERB_ERROR, 'IO', 'Write_diachro_nc4', 'type for attribute ' // yattname // &
                    ' has changed for ' // Trim( hlevel ) )
    return
  end if

  if ( ilen /= 1 ) then
    call Print_msg( NVERB_ERROR, 'IO', 'Write_diachro_nc4', 'size of attribute ' // yattname // &
                    ' has changed for ' // Trim( hlevel ) )
    return
  end if

  istatus = NF90_GET_ATT( kgrpid, NF90_GLOBAL, yattname, iatt )
  if ( iatt == kdata ) then
    call Print_msg( NVERB_DEBUG, 'IO', 'Write_diachro_nc4', 'attribute ' // yattname // ' is unchanged for ' // Trim( hlevel ) )
    !If unchanged, no need to write it again => return
    return
  else
    cmnhmsg(1) = 'attribute ' // yattname // ' has changed for ' // Trim( hlevel )
    Write( cmnhmsg(2), '( I0, " -> ", I0 )' ) iatt, kdata
    call Print_msg( NVERB_WARNING, 'IO', 'Write_diachro_nc4' )
  end if

end if

istatus = NF90_PUT_ATT( kgrpid, NF90_GLOBAL, yattname, kdata )
if (istatus /= NF90_NOERR ) &
 call IO_Err_handle_nc4( istatus, 'Write_diachro_nc4', 'NF90_PUT_ATT', Trim( yattname ) // ' for '// Trim( hlevel ) // ' group' )

end subroutine Att_write_i0


subroutine Att_write_x0( hlevel, kgrpid, hattname, pdata )
use NETCDF,            only: NF90_GET_ATT, NF90_INQUIRE_ATTRIBUTE, NF90_PUT_ATT, NF90_GLOBAL, NF90_NOERR

use modd_precision,    only: CDFINT, MNHREAL_NF90

use mode_io_tools_nc4, only: IO_Err_handle_nc4, IO_Mnhname_clean

character(len=*),     intent(in) :: hlevel
integer(kind=CDFINT), intent(in) :: kgrpid
character(len=*),     intent(in) :: hattname
real,                 intent(in) :: pdata

character(len=Len(hattname)) :: yattname
integer(kind=CDFINT) :: ilen
integer(kind=CDFINT) :: istatus
integer(kind=CDFINT) :: itype
real                 :: zatt

call IO_Mnhname_clean( hattname, yattname )

istatus = NF90_INQUIRE_ATTRIBUTE( kgrpid, NF90_GLOBAL, yattname, xtype = itype, len = ilen )
if (istatus == NF90_NOERR ) then
  call Print_msg( NVERB_DEBUG, 'IO', 'Write_diachro_nc4', 'attribute ' // yattname // ' already exists for ' // Trim( hlevel ) )

  if ( itype /= MNHREAL_NF90 ) then
    call Print_msg( NVERB_ERROR, 'IO', 'Write_diachro_nc4', 'type for attribute ' // yattname // &
                    ' has changed for ' // Trim( hlevel ) )
    return
  end if

  if ( ilen /= 1 ) then
    call Print_msg( NVERB_ERROR, 'IO', 'Write_diachro_nc4', 'size of attribute ' // yattname // &
                    ' has changed for ' // Trim( hlevel ) )
    return
  end if

  istatus = NF90_GET_ATT( kgrpid, NF90_GLOBAL, yattname, zatt )
  if ( zatt == pdata ) then
    call Print_msg( NVERB_DEBUG, 'IO', 'Write_diachro_nc4', 'attribute ' // yattname // ' is unchanged for ' // Trim( hlevel ) )
    !If unchanged, no need to write it again => return
    return
  else
    cmnhmsg(1) = 'attribute ' // yattname // ' has changed for ' // Trim( hlevel )
    Write( cmnhmsg(2), '( F15.7, " -> ", F15.7 )' ) zatt, pdata
    call Print_msg( NVERB_WARNING, 'IO', 'Write_diachro_nc4' )
  end if

end if

istatus = NF90_PUT_ATT( kgrpid, NF90_GLOBAL, yattname, pdata )
if (istatus /= NF90_NOERR ) &
 call IO_Err_handle_nc4( istatus, 'Write_diachro_nc4', 'NF90_PUT_ATT', Trim( yattname ) // ' for '// Trim( hlevel ) // ' group' )

end subroutine Att_write_x0


subroutine Move_to_next_level( kpreviouslevelid, gpreviousleveldefined, oleveluse, hlevelname, gleveldefined, klevelid )
use NETCDF,            only: NF90_DEF_GRP, NF90_INQ_NCID, NF90_NOERR

use modd_precision,    only: CDFINT

use mode_io_tools_nc4, only: IO_Err_handle_nc4, IO_Mnhname_clean

integer(kind=CDFINT), intent(in)    :: kpreviouslevelid
logical,              intent(in)    :: gpreviousleveldefined
logical,              intent(in)    :: oleveluse
! character(len=*),     intent(inout) :: hlevelname
character(len=*),     intent(in)    :: hlevelname
logical,              intent(out)   :: gleveldefined
integer(kind=CDFINT), intent(out)   :: klevelid

character(len=Len(hlevelname)) :: ylevelname
integer(kind=CDFINT) :: istatus

call IO_Mnhname_clean( hlevelname, ylevelname )

if ( oleveluse ) then
  istatus = NF90_INQ_NCID( kpreviouslevelid, Trim( ylevelname ), klevelid )
  if ( istatus == NF90_NOERR ) then
    gleveldefined = .true.
  else
    gleveldefined = .false.
    istatus = NF90_DEF_GRP( kpreviouslevelid, Trim( ylevelname ), klevelid )
    if ( istatus /= NF90_NOERR ) &
      call IO_Err_handle_nc4( istatus, 'Move_to_next_level', 'NF90_DEF_GRP', 'for ' // Trim( ylevelname ) )