From e25fb1cfcec7b92867324d9662ce9d5f0f8c2d5b Mon Sep 17 00:00:00 2001
From: Philippe WAUTELET <philippe.wautelet@aero.obs-mip.fr>
Date: Wed, 26 Jun 2024 13:56:16 +0200
Subject: [PATCH] Philippe 26/06/2024: IO_Field_read: 3D fields: ensure that
 data can not be overwritten in case of read error + other improvements

---
 src/LIB/SURCOUCHE/src/mode_allocbuff.f90     |  65 +++-
 src/LIB/SURCOUCHE/src/mode_io_field_read.f90 | 365 +++++++++----------
 2 files changed, 230 insertions(+), 200 deletions(-)

diff --git a/src/LIB/SURCOUCHE/src/mode_allocbuff.f90 b/src/LIB/SURCOUCHE/src/mode_allocbuff.f90
index 17e624971..5dc326933 100644
--- a/src/LIB/SURCOUCHE/src/mode_allocbuff.f90
+++ b/src/LIB/SURCOUCHE/src/mode_allocbuff.f90
@@ -17,11 +17,11 @@ IMPLICIT NONE
 PRIVATE
 
 INTERFACE ALLOCBUFFER_ll
-  MODULE PROCEDURE                                  &
-    ALLOCBUFFER_X1, ALLOCBUFFER_NEW_X1, ALLOCBUFFER_X2, ALLOCBUFFER_NEW_X2, ALLOCBUFFER_X3, &
+  MODULE PROCEDURE                                                                                              &
+    ALLOCBUFFER_X1, ALLOCBUFFER_NEW_X1, ALLOCBUFFER_X2, ALLOCBUFFER_NEW_X2, ALLOCBUFFER_X3, ALLOCBUFFER_NEW_X3, &
     ALLOCBUFFER_X4, ALLOCBUFFER_X5, ALLOCBUFFER_X6, &
-    ALLOCBUFFER_N1, ALLOCBUFFER_NEW_N1, ALLOCBUFFER_N2, ALLOCBUFFER_NEW_N2, ALLOCBUFFER_N3, &
-    ALLOCBUFFER_N4,                                 &
+    ALLOCBUFFER_N1, ALLOCBUFFER_NEW_N1, ALLOCBUFFER_N2, ALLOCBUFFER_NEW_N2, ALLOCBUFFER_N3, ALLOCBUFFER_NEW_N3, &
+    ALLOCBUFFER_N4,                                                                          &
     ALLOCBUFFER_L1, ALLOCBUFFER_NEW_L1
 END INTERFACE
 
@@ -171,6 +171,35 @@ CASE default
 END SELECT
 END SUBROUTINE ALLOCBUFFER_N3
 
+SUBROUTINE ALLOCBUFFER_NEW_N3( KTAB_OUT, KTAB_IN, HDIR )
+  USE MODD_IO, ONLY: LPACK, L2D
+  !
+  INTEGER, DIMENSION(:,:,:), ALLOCATABLE, INTENT(OUT) :: KTAB_OUT
+  INTEGER, DIMENSION(:,:,:),              INTENT(IN)  :: KTAB_IN
+  CHARACTER(LEN=*),                       INTENT(IN)  :: HDIR
+
+  INTEGER :: IIMAX, IJMAX
+
+  SELECT CASE(HDIR)
+  CASE('XX')
+    CALL GET_GLOBALDIMS_ll(IIMAX,IJMAX)
+    ALLOCATE(KTAB_OUT(IIMAX+2*JPHEXT,SIZE(KTAB_IN,2),SIZE(KTAB_IN,3)))
+  CASE('YY')
+    CALL GET_GLOBALDIMS_ll(IIMAX,IJMAX)
+    ALLOCATE(KTAB_OUT(IJMAX+2*JPHEXT,SIZE(KTAB_IN,2),SIZE(KTAB_IN,3)))
+  CASE('XY')
+    CALL GET_GLOBALDIMS_ll(IIMAX,IJMAX)
+    IF (LPACK .AND. L2D) THEN
+      ! 2D compact case
+      ALLOCATE(KTAB_OUT(IIMAX+2*JPHEXT,1,SIZE(KTAB_IN,3)))
+    ELSE
+      ALLOCATE(KTAB_OUT(IIMAX+2*JPHEXT,IJMAX+2*JPHEXT,SIZE(KTAB_IN,3)))
+    END IF
+  CASE default
+    ALLOCATE( KTAB_OUT, MOLD = KTAB_IN )
+  END SELECT
+END SUBROUTINE ALLOCBUFFER_NEW_N3
+
 SUBROUTINE ALLOCBUFFER_N4(KTAB_P,KTAB,HDIR,OALLOC)
 USE MODD_IO, ONLY: LPACK, L2D
 !
@@ -448,6 +477,34 @@ CASE default
 END SELECT
 END SUBROUTINE ALLOCBUFFER_X3
 
+SUBROUTINE ALLOCBUFFER_NEW_X3( PTAB_OUT, PTAB_IN, HDIR )
+  USE MODD_IO, ONLY: LPACK, L2D
+  !
+  REAL, DIMENSION(:,:,:), ALLOCATABLE, INTENT(OUT) :: PTAB_OUT
+  REAL, DIMENSION(:,:,:),              INTENT(IN)  :: PTAB_IN
+  CHARACTER(LEN=*),                    INTENT(IN)  :: HDIR
+
+  INTEGER :: IIMAX,IJMAX
+
+  SELECT CASE(HDIR)
+  CASE('XX')
+    CALL GET_GLOBALDIMS_ll(IIMAX,IJMAX)
+    ALLOCATE(PTAB_OUT(IIMAX+2*JPHEXT,SIZE(PTAB_IN,2),SIZE(PTAB_IN,3)))
+  CASE('YY')
+    CALL GET_GLOBALDIMS_ll(IIMAX,IJMAX)
+    ALLOCATE(PTAB_OUT(IJMAX+2*JPHEXT,SIZE(PTAB_IN,2),SIZE(PTAB_IN,3)))
+  CASE('XY')
+    CALL GET_GLOBALDIMS_ll(IIMAX,IJMAX)
+    IF (LPACK .AND. L2D) THEN ! 2D compact case
+      ALLOCATE(PTAB_OUT(IIMAX+2*JPHEXT,1,SIZE(PTAB_IN,3)))
+    ELSE
+      ALLOCATE(PTAB_OUT(IIMAX+2*JPHEXT,IJMAX+2*JPHEXT,SIZE(PTAB_IN,3)))
+    END IF
+  CASE default
+    ALLOCATE( PTAB_OUT, MOLD = PTAB_IN )
+  END SELECT
+END SUBROUTINE ALLOCBUFFER_NEW_X3
+
 SUBROUTINE ALLOCBUFFER_X4(PTAB_P,PTAB,HDIR,OALLOC)
 USE MODD_IO, ONLY: LPACK, L2D
 !
diff --git a/src/LIB/SURCOUCHE/src/mode_io_field_read.f90 b/src/LIB/SURCOUCHE/src/mode_io_field_read.f90
index 30a1f0987..56678d594 100644
--- a/src/LIB/SURCOUCHE/src/mode_io_field_read.f90
+++ b/src/LIB/SURCOUCHE/src/mode_io_field_read.f90
@@ -554,7 +554,6 @@ USE MODD_PARAMETERS_ll,    ONLY: JPHEXT
 USE MODD_TIMEZ,            ONLY: TIMEZ
 USE MODD_VAR_ll,           ONLY: MNH_STATUSES_IGNORE
 !
-USE MODE_ALLOCBUFFER_ll
 #ifdef MNH_GA
 USE MODE_GA
 USE MODI_GET_HALO
@@ -581,13 +580,13 @@ INTEGER                               :: JK,JKK
 INTEGER                               :: INB_REQ
 INTEGER,ALLOCATABLE,DIMENSION(:)      :: IREQ_TAB
 INTEGER, DIMENSION(MPI_STATUS_SIZE)   :: ISTATUS
-LOGICAL                               :: GALLOC, GALLOC_ll
 logical                               :: glfi, gnc4
 REAL,DIMENSION(:,:),POINTER           :: ZTX2DP
-REAL,DIMENSION(:,:),POINTER           :: ZSLICE_ll,ZSLICE
-real,dimension(:),     pointer        :: zfieldp1d
-real,dimension(:,:),   pointer        :: zfieldp2d
-REAL,DIMENSION(:,:,:), POINTER        :: ZFIELDP
+REAL, DIMENSION(:,:), ALLOCATABLE, TARGET :: ZSLICE_ll
+REAL, DIMENSION(:,:),   pointer       :: zslice
+real, dimension(:),     allocatable   :: zfield1d
+real, dimension(:,:),   allocatable   :: zfield2d
+REAL, DIMENSION(:,:,:), allocatable, target :: zfield3d
 REAL(kind=MNHTIME), DIMENSION(2)      :: ZT0, ZT1, ZT2
 REAL(kind=MNHTIME), DIMENSION(2)      :: ZT11, ZT22
 CHARACTER(LEN=2)                      :: YDIR
@@ -606,10 +605,8 @@ CALL PRINT_MSG(NVERB_DEBUG,'IO','IO_Field_read_byfield_X3',TRIM(TPFILE%CNAME)//'
 CALL SECOND_MNH2(ZT11)
 !
 TZFILE => NULL()
-GALLOC    = .FALSE.
-GALLOC_ll = .FALSE.
+zslice => null()
 IRESP = 0
-ZFIELDP => NULL()
 YDIR = TPFIELD%CDIR
 !
 IHEXTOT = 2*JPHEXT+1
@@ -619,9 +616,9 @@ CALL IO_File_read_check(TPFILE,'IO_Field_read_byfield_X3',IRESP)
 call IO_Format_read_select( tpfile, glfi, gnc4 )
 
 IF (IRESP==0) THEN
-  IF (GSMONOPROC  .AND. TPFILE%NSUBFILES_IOZ==0 ) THEN ! sequential execution
+  IF ( GSMONOPROC .AND. TPFILE%NSUBFILES_IOZ==0 ) THEN ! sequential execution
     if ( lpack .and. l1d .and. Size( pfield, 1 ) == ihextot .and. Size( pfield, 2 ) == ihextot ) then
-      Allocate( tzfield, source = tpfield )
+      allocate( tzfield, source = tpfield )
       if ( tpfile%ldimreduced ) then
         tzfield%ndims = tzfield%ndims - 2
         if ( tzfield%ndimlist(1) /= NMNHDIM_UNKNOWN ) then
@@ -629,23 +626,23 @@ IF (IRESP==0) THEN
           tzfield%ndimlist(2)  = tzfield%ndimlist(4) !Necessary if time dimension
           tzfield%ndimlist(3:) = NMNHDIM_UNUSED
         end if
-        zfieldp1d => pfield(jphext+1, jphext+1, :)
-        if ( gnc4 ) call IO_Field_read_nc4( tpfile, tzfield, zfieldp1d, iresp )
-        if ( glfi ) call IO_Field_read_lfi( tpfile, tzfield, zfieldp1d, iresp )
-        pfield(:, :, :) = Spread( Spread( pfield(jphext + 1, jphext + 1, :), dim = 1, ncopies = ihextot ), &
-                                  dim = 2, ncopies = ihextot )
+        allocate( zfield1d(size(pfield,2)) )
+        if ( gnc4 ) call IO_Field_read_nc4( tpfile, tzfield, zfield1d, iresp )
+        if ( glfi ) call IO_Field_read_lfi( tpfile, tzfield, zfield1d, iresp )
+        if ( iresp == 0 .or. iresp == -111 ) &
+          pfield(:, :, :) = Spread( Spread( zfield1d(:), dim = 1, ncopies = ihextot ), dim = 2, ncopies = ihextot )
       else
         if ( tzfield%ndimlist(1) /= NMNHDIM_UNKNOWN ) then
           tzfield%ndimlist(1:2) = NMNHDIM_ONE
         end if
-        zfieldp => pfield(jphext + 1 : jphext + 1, jphext + 1 : jphext + 1, :)
-        if ( gnc4 ) call IO_Field_read_nc4( tpfile, tzfield, zfieldp, iresp )
-        if ( glfi ) call IO_Field_read_lfi( tpfile, tzfield, zfieldp, iresp )
-        pfield(:, :, :) = Spread( Spread( pfield(jphext + 1, jphext + 1, :), dim = 1, ncopies = ihextot ), &
-                                  dim = 2, ncopies = ihextot )
+        allocate( zfield3d(1, 1, size(pfield,3)) )
+        if ( gnc4 ) call IO_Field_read_nc4( tpfile, tzfield, zfield3d, iresp )
+        if ( glfi ) call IO_Field_read_lfi( tpfile, tzfield, zfield3d, iresp )
+        if ( iresp == 0 .or. iresp == -111 ) &
+          pfield(:, :, :) = Spread( Spread( zfield3d(1, 1, :), dim = 1, ncopies = ihextot ), dim = 2, ncopies = ihextot )
       endif
     else if ( lpack .and. l2d .and. Size( pfield, 2 ) == ihextot ) then
-      Allocate( tzfield, source = tpfield )
+      allocate( tzfield, source = tpfield )
       if ( tpfile%ldimreduced ) then
         tzfield%ndims = tzfield%ndims - 1
         if ( tzfield%ndimlist(1) /= NMNHDIM_UNKNOWN ) then
@@ -653,33 +650,34 @@ IF (IRESP==0) THEN
           tzfield%ndimlist(3)  = tzfield%ndimlist(4) !Necessary if time dimension
           tzfield%ndimlist(4:) = NMNHDIM_UNUSED
         end if
-        zfieldp2d => pfield(:, jphext + 1, :)
-        if ( gnc4 ) call IO_Field_read_nc4( tpfile, tzfield, zfieldp2d, iresp )
-        if ( glfi ) call IO_Field_read_lfi( tpfile, tzfield, zfieldp2d, iresp )
-        pfield(:, :, :) = Spread( pfield(:, jphext + 1, :), dim = 2, ncopies = ihextot )
+        allocate( zfield2d(size(pfield,1), size(pfield,3)) )
+        if ( gnc4 ) call IO_Field_read_nc4( tpfile, tzfield, zfield2d, iresp )
+        if ( glfi ) call IO_Field_read_lfi( tpfile, tzfield, zfield2d, iresp )
+        if ( iresp == 0 .or. iresp == -111 ) &
+          pfield(:, :, :) = Spread( zfield2d(:, :), dim = 2, ncopies = ihextot )
       else
         if ( tzfield%ndimlist(1) /= NMNHDIM_UNKNOWN ) then
           tzfield%ndimlist(2)  = NMNHDIM_ONE
         end if
-        zfieldp => pfield(:, jphext + 1 : jphext + 1, :)
-        if ( gnc4 ) call IO_Field_read_nc4( tpfile, tzfield, zfieldp, iresp )
-        if ( glfi ) call IO_Field_read_lfi( tpfile, tzfield, zfieldp, iresp )
-        pfield(:,:, :) = Spread( pfield(:, jphext + 1, :), dim = 2, ncopies = ihextot )
+        allocate( zfield3d(size(pfield,1), 1, size(pfield,3)) )
+        if ( gnc4 ) call IO_Field_read_nc4( tpfile, tzfield, zfield3d, iresp )
+        if ( glfi ) call IO_Field_read_lfi( tpfile, tzfield, zfield3d, iresp )
+        if ( iresp == 0 .or. iresp == -111 ) &
+          pfield(:,:, :) = Spread( zfield3d(:, 1, :), dim = 2, ncopies = ihextot )
       endif
     else
       if ( gnc4 ) call IO_Field_read_nc4( tpfile, tpfield, pfield, iresp )
       if ( glfi ) call IO_Field_read_lfi( tpfile, tpfield, pfield, iresp )
     end if
-  ELSE IF ( TPFILE%NSUBFILES_IOZ==0 .OR. YDIR == '--' ) THEN ! multiprocesses execution & 1 IO proc
+  ELSE IF ( TPFILE%NSUBFILES_IOZ == 0 .OR. YDIR == '--' ) THEN ! multiprocesses execution & 1 IO proc
     IF (ISP == TPFILE%NMASTER_RANK)  THEN
       ! I/O process case
-      CALL ALLOCBUFFER_ll(ZFIELDP,PFIELD,YDIR,GALLOC)
-      if ( gnc4 ) call IO_Field_read_nc4( tpfile, tpfield, zfieldp, iresp )
-      if ( glfi ) call IO_Field_read_lfi( tpfile, tpfield, zfieldp, iresp )
+      call Allocbuffer_ll( zfield3d, pfield, ydir )
+      if ( gnc4 ) call IO_Field_read_nc4( tpfile, tpfield, zfield3d, iresp )
+      if ( glfi ) call IO_Field_read_lfi( tpfile, tpfield, zfield3d, iresp )
     ELSE
       !Not really necessary but useful to suppress alerts with Valgrind
-      ALLOCATE(ZFIELDP(0,0,0))
-      GALLOC = .TRUE.
+      allocate( zfield3d(0,0,0) )
     END IF
     !
     CALL MPI_BCAST(IRESP,1,MNHINT_MPI,TPFILE%NMASTER_RANK-1,TPFILE%NMPICOMM,IERR)
@@ -688,24 +686,30 @@ IF (IRESP==0) THEN
     !because metadata of field has been modified in IO_Field_read_xxx
     IF (IRESP==-111) CALL IO_Field_metadata_bcast(TPFILE,TPFIELD)
     !
-    IF (YDIR == 'XX' .OR. YDIR =='YY') THEN
-      ! XX or YY Scatter Field
-      CALL SCATTER_XXFIELD(YDIR,ZFIELDP,PFIELD,TPFILE%NMASTER_RANK,TPFILE%NMPICOMM)
-    ELSE IF (YDIR == 'XY') THEN
-      IF (LPACK .AND. L2D) THEN
-        ! 2D compact case
-        call Print_msg( NVERB_FATAL, 'GEN', 'IO_Field_read_byfield_X3', '2D not (yet) allowed for parallel execution' )
-        CALL SCATTER_XXFIELD('XX',ZFIELDP(:,1,:),PFIELD(:,JPHEXT+1,:),TPFILE%NMASTER_RANK,TPFILE%NMPICOMM)
-        PFIELD(:,:,:) = SPREAD(PFIELD(:,JPHEXT+1,:),DIM=2,NCOPIES=IHEXTOT)
+    !Share data only if no error
+    if ( iresp == 0 .or. iresp == -111 ) then
+      IF (YDIR == 'XX' .OR. YDIR =='YY') THEN
+        ! XX or YY Scatter Field
+        CALL SCATTER_XXFIELD(YDIR,zfield3d,PFIELD,TPFILE%NMASTER_RANK,TPFILE%NMPICOMM)
+      ELSE IF (YDIR == 'XY') THEN
+        IF (LPACK .AND. L2D) THEN
+          ! 2D compact case
+          call Print_msg( NVERB_FATAL, 'IO', 'IO_Field_read_byfield_X3', '2D not (yet) allowed for parallel execution' )
+          CALL SCATTER_XXFIELD('XX',zfield3d(:,1,:),PFIELD(:,JPHEXT+1,:),TPFILE%NMASTER_RANK,TPFILE%NMPICOMM)
+          PFIELD(:,:,:) = SPREAD(PFIELD(:,JPHEXT+1,:),DIM=2,NCOPIES=IHEXTOT)
+        ELSE
+          ! XY Scatter Field
+          CALL SCATTER_XYFIELD(zfield3d,PFIELD,TPFILE%NMASTER_RANK,TPFILE%NMPICOMM)
+        END IF
       ELSE
-        ! XY Scatter Field
-        CALL SCATTER_XYFIELD(ZFIELDP,PFIELD,TPFILE%NMASTER_RANK,TPFILE%NMPICOMM)
+        ! Broadcast Field
+        if ( isp == tpfile%nmaster_rank ) pfield(:,:,:) = zfield3d(:,:,:)
+        CALL MPI_BCAST(PFIELD,SIZE(PFIELD),MNHREAL_MPI,TPFILE%NMASTER_RANK-1,TPFILE%NMPICOMM,IERR)
       END IF
-    ELSE
-      ! Broadcast Field
-      CALL MPI_BCAST(PFIELD,SIZE(PFIELD),MNHREAL_MPI,TPFILE%NMASTER_RANK-1,TPFILE%NMPICOMM,IERR)
-    END IF
+    end if
   ELSE  ! multiprocesses execution & // IO
+    if ( ydir /= 'XY' .or. ( lpack .and. l2d ) ) &
+      call Print_msg( NVERB_FATAL, 'IO', 'IO_Field_read_byfield_X3', 'Zsplit files not supported for this case' )
 !
 !JUAN BG Z SLICE
 !
@@ -719,7 +723,6 @@ IF (IRESP==0) THEN
     ! read the data
     !
     ALLOCATE(ZSLICE_ll(0,0)) ! to avoid bug on test of size
-    GALLOC_ll = .TRUE.
     IRESP_ISP=0
     Allocate( tzfield, mold = tpfield )
     DO JKK=1,SIZE(PFIELD,3) ! IKU_ll
@@ -734,10 +737,10 @@ IF (IRESP==0) THEN
       !
       IK_RANK = TZFILE%NMASTER_RANK
       !
-      IF (ISP == IK_RANK )  THEN
+      IF (ISP == IK_RANK ) THEN
         IF ( SIZE(ZSLICE_ll) .EQ. 0 ) THEN
           DEALLOCATE(ZSLICE_ll)
-          CALL ALLOCBUFFER_ll(ZSLICE_ll,ZSLICE,YDIR,GALLOC_ll)
+          CALL ALLOCBUFFER_ll(ZSLICE_ll,ZSLICE,YDIR
         END IF
         !
         CALL SECOND_MNH2(ZT0)
@@ -766,55 +769,50 @@ IF (IRESP==0) THEN
     !
     ! get the columun data in this proc
     !
-    ! temp buf to avoid problem with none stride PFIELDS buffer  with HALO
-    ALLOCATE (ZFIELD_GA (SIZE(PFIELD,1),SIZE(PFIELD,2),SIZE(PFIELD,3)))
-    !print*,"IO_READ_FIELD_BYFIELD_X3::nga_get=",g_a, lo_col, hi_col, ld_col, TPFIELD%CMNHNAME ; call flush(6)
-    CALL NGA_GET(G_A, LO_COL, HI_COL,ZFIELD_GA(1,1,1) , LD_COL)
-    PFIELD = ZFIELD_GA
-    call ga_sync()
-    CALL GET_HALO(PFIELD)
-    DEALLOCATE(ZFIELD_GA)
+    if ( iresp == 0 .or. iresp == -111 ) then
+      ! temp buf to avoid problem with none stride PFIELDS buffer  with HALO
+      ALLOCATE (ZFIELD_GA (SIZE(PFIELD,1),SIZE(PFIELD,2),SIZE(PFIELD,3)))
+      !print*,"IO_READ_FIELD_BYFIELD_X3::nga_get=",g_a, lo_col, hi_col, ld_col, TPFIELD%CMNHNAME ; call flush(6)
+      CALL NGA_GET(G_A, LO_COL, HI_COL,ZFIELD_GA(1,1,1) , LD_COL)
+      PFIELD = ZFIELD_GA
+      call ga_sync()
+      CALL GET_HALO(PFIELD)
+    end if
 #else
-    ALLOCATE(ZSLICE_ll(0,0))
-    GALLOC_ll = .TRUE.
     IRESP_ISP=0
     INB_PROC_REAL = MIN(TPFILE%NSUBFILES_IOZ,ISNPROC)
     ALLOCATE(IREQ_TAB((ISNPROC-1)*INB_PROC_REAL))
     ALLOCATE(T_TX2DP((ISNPROC-1)*INB_PROC_REAL))
     Allocate( tzfield, mold = tpfield )
+    Allocate( zfield3d, mold = pfield )
+
     Z_SLICE: DO JK=1,SIZE(PFIELD,3),INB_PROC_REAL
-      !
-      ! read the data
-      !
+      ! Read and send the data
       JK_MAX=MIN(SIZE(PFIELD,3),JK+INB_PROC_REAL-1)
       !
       INB_REQ=0
       DO JKK=JK,JK_MAX
-        IF (TPFILE%NSUBFILES_IOZ .GT. 1 ) THEN
-          IK_FILE = IO_Level2filenumber_get(JKK,TPFILE%NSUBFILES_IOZ)
-          TZFILE => TPFILE%TFILES_IOZ(IK_FILE+1)%TFILE
-          TZFIELD = TPFIELD
-          WRITE(YSUFFIX,'(I4.4)') JKK
-          TZFIELD%CMNHNAME = TRIM(TPFIELD%CMNHNAME)//TRIM(YSUFFIX)
-          IF (LEN_TRIM(TZFIELD%CSTDNAME)>0)  TZFIELD%CSTDNAME  = TRIM(TZFIELD%CSTDNAME)//'_at_level_'//YSUFFIX
-          IF (LEN_TRIM(TZFIELD%CLONGNAME)>0) TZFIELD%CLONGNAME = TRIM(TZFIELD%CLONGNAME)//' at level '//YSUFFIX
-          TZFIELD%NDIMS = 2
-        ELSE
-          TZFILE => TPFILE
-          TZFIELD = TPFIELD
-        END IF
+        IK_FILE = IO_Level2filenumber_get(JKK,TPFILE%NSUBFILES_IOZ)
+        TZFILE => TPFILE%TFILES_IOZ(IK_FILE+1)%TFILE
+        TZFIELD = TPFIELD
+        WRITE(YSUFFIX,'(I4.4)') JKK
+        TZFIELD%CMNHNAME = TRIM(TPFIELD%CMNHNAME)//TRIM(YSUFFIX)
+        IF (LEN_TRIM(TZFIELD%CSTDNAME)>0)  TZFIELD%CSTDNAME  = TRIM(TZFIELD%CSTDNAME)//'_at_level_'//YSUFFIX
+        IF (LEN_TRIM(TZFIELD%CLONGNAME)>0) TZFIELD%CLONGNAME = TRIM(TZFIELD%CLONGNAME)//' at level '//YSUFFIX
+        TZFIELD%NDIMS = 2
         IK_RANK = TZFILE%NMASTER_RANK
-        IF (ISP == IK_RANK )  THEN
-          IF ( SIZE(ZSLICE_ll) .EQ. 0 ) THEN
-            DEALLOCATE(ZSLICE_ll)
-            CALL ALLOCBUFFER_ll(ZSLICE_ll,ZSLICE,YDIR,GALLOC_ll)
-          END IF
+        IF ( ISP == IK_RANK )  THEN
+          !Remark: zlice is not used in the following call (no problem if not allocated)
+          if ( .not.allocated( zslice_ll ) ) call Allocbuffer_ll( zslice_ll, zslice, ydir )
           CALL SECOND_MNH2(ZT0)
           WRITE(YK,'(I4.4)')  JKK
           YRECZSLICE = TRIM(TPFIELD%CMNHNAME)//YK
           call IO_Format_read_select( tzfile, glfi, gnc4 ) !Safer to do that (probably useless)
           if ( gnc4 ) call IO_Field_read_nc4( tzfile, tzfield, zslice_ll, iresp_tmp )
           if ( glfi ) call IO_Field_read_lfi( tzfile, tzfield, zslice_ll, iresp_tmp )
+          ! Do not treat special case where metadata has changed if Zsplit files
+          ! This allows also to only keep true errors
+          if ( iresp_tmp == -111 ) iresp_tmp = 0
           IF (IRESP_TMP .NE. 0 ) IRESP_ISP = IRESP_TMP
           CALL SECOND_MNH2(ZT1)
           TIMEZ%T_READ3D_READ=TIMEZ%T_READ3D_READ + ZT1 - ZT0
@@ -829,7 +827,7 @@ IF (IRESP==0) THEN
                              TZFILE%NMPICOMM,IREQ_TAB(INB_REQ),IERR)
               !CALL MPI_BSEND(ZTX2DP,SIZE(ZTX2DP),MNHREAL_MPI,JI-1,199+IK_RANK,TZFILE%NMPICOMM,IERR)
             ELSE
-              PFIELD(:,:,JKK) = ZTX2DP(:,:)
+              zfield3d(:,:,jkk) = ztx2dp(:,:)
             END IF
           END DO
           CALL SECOND_MNH2(ZT2)
@@ -838,53 +836,29 @@ IF (IRESP==0) THEN
         !
         TZFILE => NULL()
       END DO
-      !
-      ! broadcast the data
-      !
-      IF (YDIR == 'XX' .OR. YDIR =='YY') THEN
-        ! XX or YY Scatter Field
-        call Print_msg( NVERB_FATAL, 'GEN', 'IO_Field_read_byfield_X3', 'XX/YY not (yet) allowed for parallel I/O' )
-        CALL SCATTER_XXFIELD(YDIR,ZFIELDP,PFIELD,TPFILE%NMASTER_RANK,TPFILE%NMPICOMM)
-      ELSE IF (YDIR == 'XY') THEN
-        IF (LPACK .AND. L2D) THEN
-          ! 2D compact case
-          call Print_msg( NVERB_FATAL, 'GEN', 'IO_Field_read_byfield_X3', '2D not (yet) allowed for parallel execution' )
-          CALL SCATTER_XXFIELD('XX',ZFIELDP(:,1,:),PFIELD(:,JPHEXT+1,:),TPFILE%NMASTER_RANK,TPFILE%NMPICOMM)
-          PFIELD(:,:,:) = SPREAD(PFIELD(:,JPHEXT+1,:),DIM=2,NCOPIES=IHEXTOT)
-        ELSE
-          !
-          ! XY Scatter Field
-          !
-          CALL SECOND_MNH2(ZT0)
-          DO JKK=JK,JK_MAX
-            !
-            ! get the file & rank
-            !
-            IF (TPFILE%NSUBFILES_IOZ .GT. 1 ) THEN
-               IK_FILE = IO_Level2filenumber_get(JKK,TPFILE%NSUBFILES_IOZ)
-               TZFILE => TPFILE%TFILES_IOZ(IK_FILE+1)%TFILE
-            ELSE
-              TZFILE => TPFILE
-            END IF
-            !
-            IK_RANK = TZFILE%NMASTER_RANK
-            !
-            ZSLICE => PFIELD(:,:,JKK)
-            !CALL SCATTER_XYFIELD(ZSLICE_ll,ZSLICE,TZFILE%NMASTER_RANK,TZFILE%NMPICOMM)
-            IF (ISP .NE. IK_RANK) THEN
-              CALL MPI_RECV(ZSLICE,SIZE(ZSLICE),MNHREAL_MPI,IK_RANK-1,199+IK_RANK, &
-                            TZFILE%NMPICOMM,ISTATUS,IERR)
-            END IF
-            TZFILE => NULL()
-          END DO
-          CALL SECOND_MNH2(ZT1)
-          TIMEZ%T_READ3D_RECV=TIMEZ%T_READ3D_RECV + ZT1 - ZT0
+
+      ! Receive the data
+      CALL SECOND_MNH2(ZT0)
+      DO JKK=JK,JK_MAX
+        ! Get the file & rank
+        IK_FILE = IO_Level2filenumber_get(JKK,TPFILE%NSUBFILES_IOZ)
+        TZFILE => TPFILE%TFILES_IOZ(IK_FILE+1)%TFILE
+        !
+        IK_RANK = TZFILE%NMASTER_RANK
+        !
+        zslice => zfield3d(:,:,jkk)
+        !CALL SCATTER_XYFIELD(ZSLICE_ll,ZSLICE,TZFILE%NMASTER_RANK,TZFILE%NMPICOMM)
+        IF (ISP .NE. IK_RANK) THEN
+          CALL MPI_RECV(ZSLICE,SIZE(ZSLICE),MNHREAL_MPI,IK_RANK-1,199+IK_RANK, &
+                        TZFILE%NMPICOMM,ISTATUS,IERR)
         END IF
-      ELSE
-        ! Broadcast Field
-        call Print_msg( NVERB_FATAL, 'GEN', 'IO_Field_read_byfield_X3', 'broadcast field not yet planned on Blue Gene' )
-        CALL MPI_BCAST(PFIELD,SIZE(PFIELD),MNHREAL_MPI,TPFILE%NMASTER_RANK-1,TPFILE%NMPICOMM,IERR)
-      END IF
+        TZFILE => NULL()
+        zslice => null()
+      END DO
+      CALL SECOND_MNH2(ZT1)
+      TIMEZ%T_READ3D_RECV=TIMEZ%T_READ3D_RECV + ZT1 - ZT0
+
+      ! Wait for the end of the sends
       CALL SECOND_MNH2(ZT0)
       IF (INB_REQ .GT.0 ) THEN
         CALL MPI_WAITALL(INB_REQ,IREQ_TAB,MNH_STATUSES_IGNORE,IERR)
@@ -899,18 +873,20 @@ IF (IRESP==0) THEN
     !
     CALL MPI_ALLREDUCE(-ABS(IRESP_ISP),IRESP_TMP,1,MNHINT_MPI,MPI_MIN,TPFILE%NMPICOMM,IRESP)
     IF (IRESP_TMP/=0) IRESP = IRESP_TMP !Keep last "error"
+
     !Broadcast header only if IRESP==-111
     !because metadata of field has been modified in IO_Field_read_xxx
-    IF (IRESP==-111) CALL IO_Field_metadata_bcast(TPFILE,TPFIELD)
-    !
+    !DISABLED because TPFIELD metadata has not been modified if Zsplit files (we work with a copy in TZFIELD)
+    !and metadata from main file is not really used here
+    !IF (IRESP==-111) CALL IO_Field_metadata_bcast(TPFILE,TPFIELD)
+
+    !Copy data only if no error
+    if ( iresp == 0 .or. iresp == -111 ) pfield(:,:,:) = zfield3d(:,:,:)
 #endif
 !JUAN BG Z SLICE
   END IF !(GSMONOPROC)
 END IF
 !
-IF (GALLOC)    DEALLOCATE (ZFIELDP)
-IF (GALLOC_ll) DEALLOCATE (ZSLICE_ll)
-!
 IF (IRESP==-111) IRESP = 0 !-111 is not really an error (metadata has changed)
 !
 IF (PRESENT(KRESP)) KRESP = IRESP
@@ -1674,7 +1650,6 @@ USE MODD_IO,            ONLY: GSMONOPROC, ISP, LPACK, L1D, L2D
 USE MODD_PARAMETERS_ll, ONLY: JPHEXT
 USE MODD_TIMEZ,         ONLY: TIMEZ
 !
-USE MODE_ALLOCBUFFER_ll
 USE MODE_SCATTER_ll
 !
 TYPE(TFILEDATA),                   INTENT(IN)    :: TPFILE
@@ -1682,22 +1657,19 @@ CLASS(TFIELDMETADATA),             INTENT(INOUT) :: TPFIELD
 INTEGER, DIMENSION(:,:,:), TARGET, INTENT(INOUT) :: KFIELD   ! array containing the data field
 INTEGER, OPTIONAL,                 INTENT(OUT)   :: KRESP    ! return-code
 !
-INTEGER                           :: IERR
-integer, dimension(:),     pointer :: ifieldp1d
-integer, dimension(:,:),   pointer :: ifieldp2d
-INTEGER, DIMENSION(:,:,:), POINTER :: IFIELDP
-LOGICAL                           :: GALLOC
-logical                           :: glfi, gnc4
-INTEGER                           :: IRESP
-INTEGER                           :: IHEXTOT
-class(tfieldmetadata), allocatable :: tzfield
-!
+INTEGER                                :: IERR
+integer, dimension(:),     allocatable :: ifield1d
+integer, dimension(:,:),   allocatable :: ifield2d
+integer, dimension(:,:,:), allocatable :: ifield3d
+logical                                :: glfi, gnc4
+INTEGER                                :: IRESP
+INTEGER                                :: IHEXTOT
+class(tfieldmetadata),     allocatable :: tzfield
+
 CALL PRINT_MSG(NVERB_DEBUG,'IO','IO_Field_read_byfield_N3',TRIM(TPFILE%CNAME)//': reading '//TRIM(TPFIELD%CMNHNAME))
-!
-GALLOC = .FALSE.
+
 IRESP = 0
-IFIELDP => NULL()
-!
+
 IHEXTOT = 2*JPHEXT+1
 CALL IO_File_read_check(TPFILE,'IO_Field_read_byfield_N3',IRESP)
 
@@ -1706,7 +1678,7 @@ call IO_Format_read_select( tpfile, glfi, gnc4 )
 IF (IRESP==0) THEN
   IF (GSMONOPROC) THEN ! sequential execution
     if ( lpack .and. l1d .and. Size( kfield, 1 ) == ihextot .and. Size( kfield, 2 ) == ihextot ) then
-      Allocate( tzfield, source = tpfield )
+      allocate( tzfield, source = tpfield )
       if ( tpfile%ldimreduced ) then
         tzfield%ndims = tzfield%ndims - 2
         if ( tzfield%ndimlist(1) /= NMNHDIM_UNKNOWN ) then
@@ -1714,23 +1686,23 @@ IF (IRESP==0) THEN
           tzfield%ndimlist(2)  = tzfield%ndimlist(4) !Necessary if time dimension
           tzfield%ndimlist(3:) = NMNHDIM_UNUSED
         end if
-        ifieldp1d => kfield(jphext + 1, jphext + 1, :)
-        if ( gnc4 ) call IO_Field_read_nc4( tpfile, tzfield, ifieldp1d, iresp )
-        if ( glfi ) call IO_Field_read_lfi( tpfile, tzfield, ifieldp1d, iresp )
-        kfield(:, :, :) = Spread( Spread( kfield(jphext + 1, jphext + 1, :), dim = 1, ncopies = ihextot ), &
-                                  dim = 2, ncopies = ihextot )
+        allocate( ifield1d(size(kfield,2)) )
+        if ( gnc4 ) call IO_Field_read_nc4( tpfile, tzfield, ifield1d, iresp )
+        if ( glfi ) call IO_Field_read_lfi( tpfile, tzfield, ifield1d, iresp )
+        if ( iresp == 0 .or. iresp == -111 ) &
+          kfield(:, :, :) = Spread( Spread( ifield1d(:), dim = 1, ncopies = ihextot ), dim = 2, ncopies = ihextot )
       else
         if ( tzfield%ndimlist(1) /= NMNHDIM_UNKNOWN ) then
           tzfield%ndimlist(1:2) = NMNHDIM_ONE
         end if
-        ifieldp => kfield(jphext + 1 : jphext + 1, jphext + 1 : jphext + 1, :)
-        if ( gnc4 ) call IO_Field_read_nc4( tpfile, tzfield, ifieldp, iresp )
-        if ( glfi ) call IO_Field_read_lfi( tpfile, tzfield, ifieldp, iresp )
-        kfield(:, :, :) = Spread( Spread( kfield(jphext + 1, jphext + 1, :), dim = 1, ncopies = ihextot ), &
-                                  dim = 2, ncopies = ihextot )
+        allocate( ifield3d(1, 1, size(kfield,3)) )
+        if ( gnc4 ) call IO_Field_read_nc4( tpfile, tzfield, ifield3d, iresp )
+        if ( glfi ) call IO_Field_read_lfi( tpfile, tzfield, ifield3d, iresp )
+        if ( iresp == 0 .or. iresp == -111 ) &
+          kfield(:, :, :) = Spread( Spread( ifield3d(1, 1, :), dim = 1, ncopies = ihextot ), dim = 2, ncopies = ihextot )
       endif
     else if ( lpack .and. l2d .and. Size( kfield, 2 ) == ihextot ) then
-      Allocate( tzfield, source = tpfield )
+      allocate( tzfield, source = tpfield )
       if ( tpfile%ldimreduced ) then
         tzfield%ndims = tzfield%ndims - 1
         if ( tzfield%ndimlist(1) /= NMNHDIM_UNKNOWN ) then
@@ -1738,18 +1710,20 @@ IF (IRESP==0) THEN
           tzfield%ndimlist(3)  = tzfield%ndimlist(4) !Necessary if time dimension
           tzfield%ndimlist(4:) = NMNHDIM_UNUSED
         end if
-        ifieldp2d => kfield(:, jphext + 1, :)
-        if ( gnc4 ) call IO_Field_read_nc4( tpfile, tzfield, ifieldp2d, iresp )
-        if ( glfi ) call IO_Field_read_lfi( tpfile, tzfield, ifieldp2d, iresp )
-        kfield(:, :, :) = Spread( kfield(:, jphext + 1, :), dim = 2, ncopies = ihextot )
+        allocate( ifield2d(size(kfield,1), size(kfield,3)) )
+        if ( gnc4 ) call IO_Field_read_nc4( tpfile, tzfield, ifield2d, iresp )
+        if ( glfi ) call IO_Field_read_lfi( tpfile, tzfield, ifield2d, iresp )
+        if ( iresp == 0 .or. iresp == -111 ) &
+          kfield(:, :, :) = Spread( ifield2d(:, :), dim = 2, ncopies = ihextot )
       else
         if ( tzfield%ndimlist(1) /= NMNHDIM_UNKNOWN ) then
           tzfield%ndimlist(2)  = NMNHDIM_ONE
         end if
-        ifieldp => kfield(:, jphext + 1 : jphext + 1, :)
-        if ( gnc4 ) call IO_Field_read_nc4( tpfile, tzfield, ifieldp, iresp )
-        if ( glfi ) call IO_Field_read_lfi( tpfile, tzfield, ifieldp, iresp )
-        kfield(:, :, :) = Spread( kfield(:, jphext + 1, :), dim = 2, ncopies = ihextot )
+        allocate( ifield3d(size(kfield,1), 1, size(kfield,3)) )
+        if ( gnc4 ) call IO_Field_read_nc4( tpfile, tzfield, ifield3d, iresp )
+        if ( glfi ) call IO_Field_read_lfi( tpfile, tzfield, ifield3d, iresp )
+        if ( iresp == 0 .or. iresp == -111 ) &
+          kfield(:, :, :) = Spread( ifield3d(:, 1, :), dim = 2, ncopies = ihextot )
       endif
     else
       if ( gnc4 ) call IO_Field_read_nc4( tpfile, tpfield, kfield, iresp )
@@ -1758,13 +1732,12 @@ IF (IRESP==0) THEN
   ELSE
     IF (ISP == TPFILE%NMASTER_RANK)  THEN
       ! I/O process case
-      CALL ALLOCBUFFER_ll(IFIELDP,KFIELD,TPFIELD%CDIR,GALLOC)
-      if ( gnc4 ) call IO_Field_read_nc4( tpfile, tpfield, ifieldp, iresp )
-      if ( glfi ) call IO_Field_read_lfi( tpfile, tpfield, ifieldp, iresp )
+      call Allocbuffer_ll( ifield3d, kfield, tpfield%cdir )
+      if ( gnc4 ) call IO_Field_read_nc4( tpfile, tpfield, ifield3d, iresp )
+      if ( glfi ) call IO_Field_read_lfi( tpfile, tpfield, ifield3d, iresp )
     ELSE
       !Not really necessary but useful to suppress alerts with Valgrind
-      ALLOCATE(IFIELDP(0,0,0))
-      GALLOC = .TRUE.
+      allocate( ifield3d(0,0,0) )
     END IF
     !
     CALL MPI_BCAST(IRESP,1,MNHINT_MPI,TPFILE%NMASTER_RANK-1,TPFILE%NMPICOMM,IERR)
@@ -1773,30 +1746,30 @@ IF (IRESP==0) THEN
     !because metadata of field has been modified in IO_Field_read_xxx
     IF (IRESP==-111) CALL IO_Field_metadata_bcast(TPFILE,TPFIELD)
     !
-    IF (TPFIELD%CDIR == 'XX' .OR. TPFIELD%CDIR == 'YY') THEN
-      ! XX or YY Scatter Field
-      CALL SCATTER_XXFIELD(TPFIELD%CDIR,IFIELDP,KFIELD,TPFILE%NMASTER_RANK,TPFILE%NMPICOMM)
-      ! Broadcast Field
-      CALL MPI_BCAST(KFIELD,SIZE(KFIELD),MNHREAL_MPI,TPFILE%NMASTER_RANK-1,TPFILE%NMPICOMM,IERR)
-    ELSE IF (TPFIELD%CDIR == 'XY') THEN
-      IF (LPACK .AND. L2D) THEN
-        ! 2D compact case
-        call Print_msg( NVERB_FATAL, 'GEN', 'IO_Field_read_byfield_N3', '2D not (yet) allowed for parallel execution' )
-        CALL SCATTER_XXFIELD('XX',IFIELDP(:,1,:),KFIELD(:,JPHEXT+1,:),TPFILE%NMASTER_RANK,TPFILE%NMPICOMM)
-        KFIELD(:,:,:) = SPREAD(KFIELD(:,JPHEXT+1,:),DIM=2,NCOPIES=IHEXTOT)
+    !Share data only if no error
+    if ( iresp == 0 .or. iresp == -111 ) then
+      IF (TPFIELD%CDIR == 'XX' .OR. TPFIELD%CDIR == 'YY') THEN
+        ! XX or YY Scatter Field
+        CALL SCATTER_XXFIELD(TPFIELD%CDIR,ifield3d,KFIELD,TPFILE%NMASTER_RANK,TPFILE%NMPICOMM)
+      ELSE IF (TPFIELD%CDIR == 'XY') THEN
+        IF (LPACK .AND. L2D) THEN
+          ! 2D compact case
+          call Print_msg( NVERB_FATAL, 'GEN', 'IO_Field_read_byfield_N3', '2D not (yet) allowed for parallel execution' )
+          CALL SCATTER_XXFIELD('XX',ifield3d(:,1,:),KFIELD(:,JPHEXT+1,:),TPFILE%NMASTER_RANK,TPFILE%NMPICOMM)
+          KFIELD(:,:,:) = SPREAD(KFIELD(:,JPHEXT+1,:),DIM=2,NCOPIES=IHEXTOT)
+        ELSE
+          ! XY Scatter Field
+          CALL SCATTER_XYFIELD(ifield3d,KFIELD,TPFILE%NMASTER_RANK,TPFILE%NMPICOMM)
+        END IF
       ELSE
-        ! XY Scatter Field
-        CALL SCATTER_XYFIELD(IFIELDP,KFIELD,TPFILE%NMASTER_RANK,TPFILE%NMPICOMM)
+        ! Broadcast Field
+        if ( isp == tpfile%nmaster_rank ) kfield(:,:,:) = ifield3d(:,:,:)
+        CALL MPI_BCAST(KFIELD,SIZE(KFIELD),MNHINT_MPI,TPFILE%NMASTER_RANK-1,TPFILE%NMPICOMM,IERR)
       END IF
-    ELSE
-      IF (ISP == TPFILE%NMASTER_RANK) KFIELD = IFIELDP
-      CALL MPI_BCAST(KFIELD,SIZE(KFIELD),MNHINT_MPI,TPFILE%NMASTER_RANK-1,TPFILE%NMPICOMM,IERR)
-    END IF
+    end if
   END IF
 END IF
 !
-IF (GALLOC) DEALLOCATE (IFIELDP)
-!
 IF (IRESP==-111) IRESP = 0 !-111 is not really an error (metadata has changed)
 !
 IF (PRESENT(KRESP)) KRESP = IRESP
-- 
GitLab