Skip to content
Snippets Groups Projects
write_diachro.f90 67.9 KiB
Newer Older
  • Learn to ignore specific revisions
  •       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, tpfields(1), htype, pvar, [ 3, 4, 5 ], gsplit, gdistributed, &
                                            kil, kih, kjl, kjh, kkl, kkh )
    
        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, tpfields(ji), htype, 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, tpfields(ji), htype, pvar(:,:,:,:,:,ji:ji), [ 3, 4 ], 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, tpfields(1), htype, pvar, [ 3, 4, 5 ], gsplit, gdistributed, &
                                            kil, kih, kjl, kjh, kkl, kkh )
    
        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, tpfields(1), htype, pvar, [ 1, 3, 4 ], gsplit, gdistributed, &
                                            kil, kih, kjl, kjh, kkl, kkh )
    
        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, tpfields(1), htype, pvar, [ 2, 3, 4 ], gsplit, gdistributed, &
                                            kil, kih, kjl, kjh, kkl, kkh )
    
        else if (  tpfields(1)%ndimlist(3) == NMNHDIM_LEVEL      &
             .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_DEF_DIM( igrpid, 'time_flyer', Size( pvar, 4), idimid )
            if ( istatus /= NF90_NOERR ) &
              call IO_Err_handle_nc4( istatus, 'Write_diachro_nc4', 'NF90_DEF_DIM', Trim( tpfields(1)%cmnhname ) )
          end if
    
    
          ! Loop on the processes
          do ji = 1, Size( pvar, 6 )
            call Diachro_one_field_write_nc4( tzfile, tpfields(ji), htype, pvar(:,:,:,:,:,ji:ji), [ 3, 4 ], gsplit, gdistributed )
    
          end do
        else if (  tpfields(1)%ndimlist(3) == NMNHDIM_LEVEL      &
             .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, tpfields(ji), htype, pvar(:,:,:,:,:,ji:ji), [ 3, 4 ], gsplit, gdistributed )
    
        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, tpfields(ji), htype, 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, tpfields(ji), htype, pvar(:,:,:,:,:,ji:ji), [ 1, 4 ], gsplit, gdistributed )
    
        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
    
            call Diachro_one_field_write_nc4( tzfile, tpfields(ji), htype, pvar(:,:,:,:,:,ji:ji), [ 4, 5 ], gsplit, gdistributed )
    
        else
          call Print_msg( NVERB_FATAL, 'IO', 'Write_diachro_nc4', &
                          'case not yet implemented (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)
    
          do ji = 1, Size( pvar, 6 )
            call Diachro_one_field_write_nc4( tzfile, tpfields(ji), htype, pvar(:,:,:,:,:,ji:ji), [ 1, 2, 3 ], gsplit, gdistributed, &
                                              kil, kih, kjl, kjh, kkl, kkh )
    
        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, tpfields(ji), htype, pvar(:,:,:,:,:,ji:ji), [ 3, 4, 5 ], gsplit, gdistributed, &
                                              kil, kih, kjl, kjh, kkl, kkh )
    
        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, tpfields(ji), htype, pvar(:,:,:,:,:,ji:ji), [ 3, 4, 5 ], gsplit, gdistributed, &
                                              kil, kih, kjl, kjh, kkl, kkh )
    
        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, tpfields(ji), htype, pvar(:,:,:,:,:,ji:ji), [ 3, 4, 5 ], gsplit, gdistributed, &
                                              kil, kih, kjl, kjh, kkl, kkh )
    
        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, tpfields(1), htype, 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, tpfields(1), htype, pvar, [ 2, 3, 4, 5 ], gsplit, gdistributed )
    
        else
          call Print_msg( NVERB_FATAL, 'IO', 'Write_diachro_nc4', &
                          'case not yet implemented (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, tpfields(ji), htype, pvar(:,:,:,:,:,ji:ji), [ 1, 2, 3, 4, 5 ], &
                                            gsplit, gdistributed )
    
    
    !Restore id of the file root group ('/' group)
    tzfile%nncid = isavencid
    
    end  subroutine Write_diachro_nc4
    
    subroutine Diachro_one_field_write_nc4( tpfile, tpfield, htype, pvar, kdims, osplit, odistributed, kil, kih, kjl, kjh, kkl, kkh )
    use modd_budget,      only: nbutshift, nbusubwrite
    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
    class(tfield_metadata_base),                         intent(in)  :: tpfield
    character(len=*),                                    intent(in)  :: htype
    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 /= 3 )                                                                                  &
        call Print_msg( NVERB_FATAL, 'IO', 'Diachro_one_field_write_nc4',                                &
                       'odistributed=.true. not allowed for dims/=3, field: ' //Trim( tzfield%cmnhname ) )
    
      if ( htype /= 'CART' )                                                                                 &
        call Print_msg( NVERB_FATAL, 'IO', 'Diachro_one_field_write_nc4',                                    &
                       'odistributed=.true. not allowed for htype/=CART, 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 ( htype /= 'CART' .and. htype /= 'MASK' )                                                                 &
        call Print_msg( NVERB_FATAL, 'IO', 'Diachro_one_field_write_nc4',                                          &
                        'osplit=.true. not allowed for htype/=CART and /=MASK, 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 ( 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
        else !.not. osplit
          call IO_Field_write( tpfile, tzfield, zdata2d )
        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
    
    
    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
    
    
    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