From 005ca9c1bbc017aa97d8201a7a0b6c76af70a4ad Mon Sep 17 00:00:00 2001
From: Sebastien MASSON <massons@irene151.c-irene.tgcc.ccc.cea.fr>
Date: Thu, 1 Sep 2022 12:43:32 +0200
Subject: [PATCH] fix lib_fortran routines for all halo size, #68

---
 src/OCE/lib_fortran.F90         | 272 ++++----------------------------
 src/OCE/lib_fortran_generic.h90 | 245 ++++++++++++++++++----------
 2 files changed, 193 insertions(+), 324 deletions(-)

diff --git a/src/OCE/lib_fortran.F90 b/src/OCE/lib_fortran.F90
index a9857b8ae..ded29e259 100644
--- a/src/OCE/lib_fortran.F90
+++ b/src/OCE/lib_fortran.F90
@@ -88,6 +88,22 @@ CONTAINS
 #     define DIM_3d
 #        include "lib_fortran_generic.h90"
 #     undef DIM_3d
+#  define LOCALONLY
+#     define DIM_2d
+#        include "lib_fortran_generic.h90"
+#     undef DIM_2d
+#     define DIM_3d
+#        include "lib_fortran_generic.h90"
+#     undef DIM_3d
+#  undef LOCALONLY
+#  define VEC
+#     define DIM_3d
+#        include "lib_fortran_generic.h90"
+#     undef DIM_3d
+#     define DIM_4d
+#        include "lib_fortran_generic.h90"
+#     undef DIM_4d
+#  undef VEC
 #  undef GLOBSUM_CODE
 
 #  define GLOBMINMAX_CODE
@@ -107,71 +123,26 @@ CONTAINS
 #           include "lib_fortran_generic.h90"
 #        undef OPERATION_GLOBMAX
 #     undef DIM_3
+#  define VEC
+#     define DIM_3d
+#        define OPERATION_GLOBMIN
+#           include "lib_fortran_generic.h90"
+#        undef OPERATION_GLOBMIN
+#        define OPERATION_GLOBMAX
+#           include "lib_fortran_generic.h90"
+#        undef OPERATION_GLOBMAX
+#     undef DIM_3d
+#     define DIM_4d
+#        define OPERATION_GLOBMIN
+#           include "lib_fortran_generic.h90"
+#        undef OPERATION_GLOBMIN
+#        define OPERATION_GLOBMAX
+#           include "lib_fortran_generic.h90"
+#        undef OPERATION_GLOBMAX
+#     undef DIM_4d
+#  undef VEC
 #  undef GLOBMINMAX_CODE
 
-!                          ! FUNCTION local_sum !
-
-   FUNCTION local_sum_2d( ptab )
-      !!----------------------------------------------------------------------
-      REAL(wp),  INTENT(in   ) ::   ptab(:,:) ! array on which operation is applied
-      COMPLEX(dp)              ::  local_sum_2d
-      !
-      !!-----------------------------------------------------------------------
-      !
-      COMPLEX(dp)::   ctmp
-      REAL(wp)   ::   ztmp
-      INTEGER    ::   ji, jj    ! dummy loop indices
-      INTEGER    ::   ipi, ipj  ! dimensions
-      !!-----------------------------------------------------------------------
-      !
-      ipi = SIZE(ptab,1)   ! 1st dimension
-      ipj = SIZE(ptab,2)   ! 2nd dimension
-      !
-      ctmp = CMPLX( 0.e0, 0.e0, wp )   ! warning ctmp is cumulated
-
-      DO jj = 1, ipj
-         DO ji = 1, ipi
-            ztmp =  ptab(ji,jj) * tmask_i(ji,jj)
-            CALL DDPDD( CMPLX( ztmp, 0.e0, dp ), ctmp )
-         END DO
-      END DO
-      !
-      local_sum_2d = ctmp
-       
-   END FUNCTION local_sum_2d
-
-   FUNCTION local_sum_3d( ptab )
-      !!----------------------------------------------------------------------
-      REAL(wp),  INTENT(in   ) ::   ptab(:,:,:) ! array on which operation is applied
-      COMPLEX(dp)              ::  local_sum_3d
-      !
-      !!-----------------------------------------------------------------------
-      !
-      COMPLEX(dp)::   ctmp
-      REAL(wp)   ::   ztmp
-      INTEGER    ::   ji, jj, jk   ! dummy loop indices
-      INTEGER    ::   ipi, ipj, ipk    ! dimensions
-      !!-----------------------------------------------------------------------
-      !
-      ipi = SIZE(ptab,1)   ! 1st dimension
-      ipj = SIZE(ptab,2)   ! 2nd dimension
-      ipk = SIZE(ptab,3)   ! 3rd dimension
-      !
-      ctmp = CMPLX( 0.e0, 0.e0, wp )   ! warning ctmp is cumulated
-
-      DO jk = 1, ipk
-        DO jj = 1, ipj
-          DO ji = 1, ipi
-             ztmp =  ptab(ji,jj,jk) * tmask_i(ji,jj)
-             CALL DDPDD( CMPLX( ztmp, 0.e0, dp ), ctmp )
-          END DO
-        END DO
-      END DO
-      !
-      local_sum_3d = ctmp
-       
-   END FUNCTION local_sum_3d
-
 !                          ! FUNCTION sum3x3 !
 
    SUBROUTINE sum3x3_2d( p2d )
@@ -283,181 +254,6 @@ CONTAINS
 
    END SUBROUTINE sum3x3_3d
 
-
-   FUNCTION glob_sum_vec_3d( cdname, ptab ) RESULT( ptmp )
-      !!----------------------------------------------------------------------
-      CHARACTER(len=*),  INTENT(in) ::   cdname      ! name of the calling subroutine
-      REAL(wp),          INTENT(in) ::   ptab(:,:,:) ! array on which operation is applied
-      REAL(wp), DIMENSION(SIZE(ptab,3)) ::   ptmp
-      !
-      COMPLEX(dp), DIMENSION(:), ALLOCATABLE ::   ctmp
-      REAL(wp)    ::   ztmp
-      INTEGER     ::   ji , jj , jk     ! dummy loop indices
-      INTEGER     ::   ipi, ipj, ipk    ! dimensions
-      INTEGER     ::   iis, iie, ijs, ije   ! loop start and end
-      !!-----------------------------------------------------------------------
-      !
-      ipi = SIZE(ptab,1)   ! 1st dimension
-      ipj = SIZE(ptab,2)   ! 2nd dimension
-      ipk = SIZE(ptab,3)   ! 3rd dimension
-      !
-      IF( ipi == jpi .AND. ipj == jpj ) THEN   ! do 2D loop only over the inner domain (-> avoid to use undefined values)
-         iis = Nis0   ;   iie = Nie0
-         ijs = Njs0   ;   ije = Nje0
-      ELSE                                     ! if ptab is not defined over the whole domain (-> avoid out of bounds)
-         iis = 1   ;   iie = jpi-2*nn_hls
-         ijs = 1   ;   ije = jpj-2*nn_hls
-      ENDIF
-      !
-      ALLOCATE( ctmp(ipk) )
-      !
-      DO jk = 1, ipk
-         ctmp(jk) = CMPLX( 0.e0, 0.e0, dp )   ! warning ctmp is cumulated
-         DO jj = ijs, ije
-            DO ji = iis, iie
-               ztmp =  ptab(ji,jj,jk) * tmask_i(ji,jj)
-               CALL DDPDD( CMPLX( ztmp, 0.e0, dp ), ctmp(jk) )
-            END DO
-         END DO
-      END DO
-      CALL mpp_sum( cdname, ctmp(:) )   ! sum over the global domain
-      !
-      ptmp = REAL( ctmp(:), wp )
-      !
-      DEALLOCATE( ctmp )
-      !
-   END FUNCTION glob_sum_vec_3d
-
-   FUNCTION glob_sum_vec_4d( cdname, ptab ) RESULT( ptmp )
-      !!----------------------------------------------------------------------
-      CHARACTER(len=*),  INTENT(in) ::   cdname        ! name of the calling subroutine
-      REAL(wp),          INTENT(in) ::   ptab(:,:,:,:) ! array on which operation is applied
-      REAL(wp), DIMENSION(SIZE(ptab,4)) ::   ptmp
-      !
-      COMPLEX(dp), DIMENSION(:), ALLOCATABLE ::   ctmp
-      REAL(wp)    ::   ztmp
-      INTEGER     ::   ji , jj , jk , jl     ! dummy loop indices
-      INTEGER     ::   ipi, ipj, ipk, ipl    ! dimensions
-      INTEGER     ::   iis, iie, ijs, ije    ! loop start and end
-      !!-----------------------------------------------------------------------
-      !
-      ipi = SIZE(ptab,1)   ! 1st dimension
-      ipj = SIZE(ptab,2)   ! 2nd dimension
-      ipk = SIZE(ptab,3)   ! 3rd dimension
-      ipl = SIZE(ptab,4)   ! 4th dimension
-      !
-      IF( ipi == jpi .AND. ipj == jpj ) THEN   ! do 2D loop only over the inner domain (-> avoid to use undefined values)
-         iis = Nis0   ;   iie = Nie0
-         ijs = Njs0   ;   ije = Nje0
-      ELSE                                     ! if ptab is not defined over the whole domain (-> avoid out of bounds)
-         iis = 1   ;   iie = jpi-2*nn_hls
-         ijs = 1   ;   ije = jpj-2*nn_hls
-      ENDIF
-      !
-      ALLOCATE( ctmp(ipl) )
-      !
-      DO jl = 1, ipl
-         ctmp(jl) = CMPLX( 0.e0, 0.e0, dp )   ! warning ctmp is cumulated
-         DO jk = 1, ipk
-            DO jj = ijs, ije
-               DO ji = iis, iie
-                  ztmp =  ptab(ji,jj,jk,jl) * tmask_i(ji,jj)
-                  CALL DDPDD( CMPLX( ztmp, 0.e0, dp ), ctmp(jl) )
-               END DO
-            END DO
-         END DO
-      END DO
-      CALL mpp_sum( cdname, ctmp(:) )   ! sum over the global domain
-      !
-      ptmp = REAL( ctmp(:), wp )
-      !
-      DEALLOCATE( ctmp )
-      !
-   END FUNCTION glob_sum_vec_4d
-
-   FUNCTION glob_min_vec_3d( cdname, ptab ) RESULT( ptmp )
-      !!----------------------------------------------------------------------
-      CHARACTER(len=*),  INTENT(in) ::   cdname        ! name of the calling subroutine
-      REAL(wp),          INTENT(in) ::   ptab(:,:,:)   ! array on which operation is applied
-      REAL(wp), DIMENSION(SIZE(ptab,3)) ::   ptmp
-      !
-      INTEGER     ::   jk    ! dummy loop indice & dimension
-      INTEGER     ::   ipk   ! dimension
-      !!-----------------------------------------------------------------------
-      !
-      ipk = SIZE(ptab,3)
-      DO jk = 1, ipk
-         ptmp(jk) = MINVAL( ptab(:,:,jk) * tmask_i(:,:) )
-      ENDDO
-      !
-      CALL mpp_min( cdname, ptmp (:) )
-      !
-   END FUNCTION glob_min_vec_3d
-
-   FUNCTION glob_min_vec_4d( cdname, ptab ) RESULT( ptmp )
-      !!----------------------------------------------------------------------
-      CHARACTER(len=*),  INTENT(in) ::   cdname          ! name of the calling subroutine
-      REAL(wp),          INTENT(in) ::   ptab(:,:,:,:)   ! array on which operation is applied
-      REAL(wp), DIMENSION(SIZE(ptab,4)) ::   ptmp
-      !
-      INTEGER     ::   jk , jl    ! dummy loop indice & dimension
-      INTEGER     ::   ipk, ipl   ! dimension
-      !!-----------------------------------------------------------------------
-      !
-      ipk = SIZE(ptab,3)
-      ipl = SIZE(ptab,4)
-      DO jl = 1, ipl
-            ptmp(jl) = MINVAL( ptab(:,:,1,jl) * tmask_i(:,:) )         
-         DO jk = 2, ipk
-            ptmp(jl) = MIN( ptmp(jl), MINVAL( ptab(:,:,jk,jl) * tmask_i(:,:) ) )
-         ENDDO
-      ENDDO
-      !
-      CALL mpp_min( cdname, ptmp (:) )
-      !
-   END FUNCTION glob_min_vec_4d
-   
-   FUNCTION glob_max_vec_3d( cdname, ptab ) RESULT( ptmp )
-      !!----------------------------------------------------------------------
-      CHARACTER(len=*),  INTENT(in) ::   cdname        ! name of the calling subroutine
-      REAL(wp),          INTENT(in) ::   ptab(:,:,:)   ! array on which operation is applied
-      REAL(wp), DIMENSION(SIZE(ptab,3)) ::   ptmp
-      !
-      INTEGER     ::   jk    ! dummy loop indice & dimension
-      INTEGER     ::   ipk   ! dimension
-      !!-----------------------------------------------------------------------
-      !
-      ipk = SIZE(ptab,3)
-      DO jk = 1, ipk
-         ptmp(jk) = MAXVAL( ptab(:,:,jk) * tmask_i(:,:) )
-      ENDDO
-      !
-      CALL mpp_max( cdname, ptmp (:) )
-      !
-   END FUNCTION glob_max_vec_3d
-
-   FUNCTION glob_max_vec_4d( cdname, ptab ) RESULT( ptmp )
-      !!----------------------------------------------------------------------
-      CHARACTER(len=*),  INTENT(in) ::   cdname          ! name of the calling subroutine
-      REAL(wp),          INTENT(in) ::   ptab(:,:,:,:)   ! array on which operation is applied
-      REAL(wp), DIMENSION(SIZE(ptab,4)) ::   ptmp
-      !
-      INTEGER     ::   jk , jl    ! dummy loop indice & dimension
-      INTEGER     ::   ipk, ipl   ! dimension
-      !!-----------------------------------------------------------------------
-      !
-      ipk = SIZE(ptab,3)
-      ipl = SIZE(ptab,4)
-      DO jl = 1, ipl
-            ptmp(jl) = MAXVAL( ptab(:,:,1,jl) * tmask_i(:,:) )         
-         DO jk = 2, ipk
-            ptmp(jl) = MAX( ptmp(jl), MAXVAL( ptab(:,:,jk,jl) * tmask_i(:,:) ) )
-         ENDDO
-      ENDDO
-      !
-      CALL mpp_max( cdname, ptmp (:) )
-      !
-   END FUNCTION glob_max_vec_4d
    
    SUBROUTINE DDPDD( ydda, yddb )
       !!----------------------------------------------------------------------
diff --git a/src/OCE/lib_fortran_generic.h90 b/src/OCE/lib_fortran_generic.h90
index e0702190b..2b84ed31a 100644
--- a/src/OCE/lib_fortran_generic.h90
+++ b/src/OCE/lib_fortran_generic.h90
@@ -1,132 +1,205 @@
-#if defined GLOBSUM_CODE
-!                          ! FUNCTION FUNCTION_GLOBSUM !
+/**/
+/*-----------------------------*/
+/*   DEFINE COMMON VARIABLES   */
+/*-----------------------------*/
+/**/
 #   if defined DIM_1d
-#      define XD                1d
-#      define ARRAY_TYPE(i,j,k)    REAL(wp)                 , INTENT(in   ) ::   ARRAY_IN(i,j,k)
-#      define ARRAY_IN(i,j,k)   ptab(i)
-#      define K_SIZE(ptab)      1
-#      define MASK_ARRAY(i,j)   1.
+#      define XD                  1d
+#      define ARRAY_IN(i,j,k,l)   ptab(i)
 #   endif
 #   if defined DIM_2d
-#      define XD                2d
-#      define ARRAY_TYPE(i,j,k)    REAL(wp)                 , INTENT(in   ) ::   ARRAY_IN(i,j,k)
-#      define ARRAY_IN(i,j,k)   ptab(i,j)
-#      define K_SIZE(ptab)      1
-#      define MASK_ARRAY(i,j)   tmask_i(i,j)
+#      define XD                  2d
+#      define ARRAY_IN(i,j,k,l)   ptab(i,j)
+#      define K_SIZE(ptab)        1
+#      define L_SIZE(ptab)        1
+#      define LAST_SIZE           -1
 #   endif
 #   if defined DIM_3d
-#      define XD                3d
-#      define ARRAY_TYPE(i,j,k)    REAL(wp)                 , INTENT(in   ) ::   ARRAY_IN(i,j,k)
-#      define ARRAY_IN(i,j,k)   ptab(i,j,k)
-#      define K_SIZE(ptab)      SIZE(ptab,3)
-#      define MASK_ARRAY(i,j)   tmask_i(i,j)
-#   endif
-
-   FUNCTION glob_sum_/**/XD/**/( cdname, ptab )
+#      define XD                  3d
+#      define ARRAY_IN(i,j,k,l)   ptab(i,j,k)
+#      define K_SIZE(ptab)        SIZE(ptab,3)
+#      define L_SIZE(ptab)        1
+#      define LAST_SIZE           SIZE(ptab,3)
+#   endif
+#   if defined DIM_4d
+#      define XD                  4d
+#      define ARRAY_IN(i,j,k,l)   ptab(i,j,k,l)
+#      define K_SIZE(ptab)        SIZE(ptab,3)
+#      define L_SIZE(ptab)        SIZE(ptab,4)
+#      define LAST_SIZE           SIZE(ptab,4)
+#   endif
+#   if defined VEC
+#      define ISVEC               _vec
+#   else
+#      define ISVEC               
+#   endif
+#   if defined LOCALONLY
+#      define TYPENAME           local
+#   else
+#      define TYPENAME           glob
+#   endif
+/**/
+/*-------------------------------*/
+/*   FUNCTION FUNCTION_GLOBSUM   */
+/*-------------------------------*/
+/**/
+#if defined GLOBSUM_CODE
+/**/
+/*   DEFINE LOCAL VARIABLES   */
+/**/
+!
+#   if defined LOCALONLY
+FUNCTION TYPENAME/**/_sum/**/ISVEC/**/_/**/XD/**/(         ptab ) RESULT( ptmp )
       !!----------------------------------------------------------------------
-      CHARACTER(len=*),  INTENT(in   ) ::   cdname  ! name of the calling subroutine
-      ARRAY_TYPE(:,:,:)                             ! array on which operation is applied
-      REAL(wp)   ::  glob_sum_/**/XD
-      !
-      !!-----------------------------------------------------------------------
+#   else
+FUNCTION TYPENAME/**/_sum/**/ISVEC/**/_/**/XD/**/( cdname, ptab ) RESULT( ptmp )
+      !!----------------------------------------------------------------------
+      CHARACTER(len=*), INTENT(in   ) ::   cdname  ! name of the calling subroutine
+#   endif
+      REAL(wp)        , INTENT(in   ) ::   ARRAY_IN(:,:,:,:)   ! array on which operation is applied
       !
+#   if defined VEC
+      REAL(wp)   , DIMENSION(LAST_SIZE) ::   ptmp
+      COMPLEX(dp), DIMENSION(LAST_SIZE) ::   ctmp
+#   else
+      REAL(wp)   ::  ptmp
       COMPLEX(dp)::   ctmp
-      REAL(wp)   ::   ztmp
-      INTEGER    ::   ji, jj, jk           ! dummy loop indices
-      INTEGER    ::   ipi, ipj, ipk        ! dimensions
+#   endif
+      INTEGER    ::    ji,  jj,  jk,  jl        ! dummy loop indices
+      INTEGER    ::   ipi, ipj, ipk, ipl        ! dimensions
       INTEGER    ::   iisht, ijsht
       !!-----------------------------------------------------------------------
       !
+#  if defined DIM_1d
       ctmp = CMPLX( 0.e0, 0.e0, dp )   ! warning ctmp is cumulated
-      !
-#   if defined DIM_1d
       DO ji = 1, SIZE(ptab,1)
          CALL DDPDD( CMPLX( ptab(ji), 0.e0, dp ), ctmp )
       END DO
-#   else
+      !
+#  else
       ipi = SIZE(ptab,1)   ! 1st dimension
       ipj = SIZE(ptab,2)   ! 2nd dimension
       ipk = K_SIZE(ptab)   ! 3rd dimension
+      ipl = L_SIZE(ptab)   ! 4th dimension
       !
       iisht = ( jpi - ipi ) / 2
       ijsht = ( jpj - ipj ) / 2   ! should be the same as iisht...
       !
-      DO jk = 1, ipk
-         DO_2D( 0, 0, 0, 0 )
-            ztmp = ARRAY_IN(ji-iisht,jj-ijsht,jk) * MASK_ARRAY(ji,jj)   ! warning tmask_iis defined over the full MPI domain
-            CALL DDPDD( CMPLX( ztmp, 0.e0, dp ), ctmp )
-         END_2D
-      END DO
+      ctmp = CMPLX( 0.e0, 0.e0, dp )   ! warning ctmp is cumulated
+      !
+      DO jl = 1, ipl
+         DO jk = 1, ipk
+            DO_2D( 0, 0, 0, 0 )
+               ! warning tmask_i is defined over the full MPI domain but maybe not ptab
+#   define ARRAY_LOOP             ARRAY_IN(ji-iisht,jj-ijsht,jk,jl) * tmask_i(ji,jj)
+#   if   defined VEC && defined DIM_3d
+               CALL DDPDD( CMPLX( ARRAY_LOOP, 0.e0, dp ), ctmp(jk) )
+#   endif
+#   if   defined VEC && defined DIM_4d
+               CALL DDPDD( CMPLX( ARRAY_LOOP, 0.e0, dp ), ctmp(jl) )
+#   endif
+#   if ! defined VEC
+               CALL DDPDD( CMPLX( ARRAY_LOOP, 0.e0, dp ), ctmp )
 #   endif
+            END_2D
+         END DO
+      END DO
+      !
+#  endif
+#  if defined LOCALONLY
+      ptmp = ctmp
+#  else       
       CALL mpp_sum( cdname, ctmp )   ! sum over the global domain
-      glob_sum_/**/XD = REAL(ctmp,wp)
-
-   END FUNCTION glob_sum_/**/XD
-
-#undef XD
-#undef ARRAY_TYPE
-#undef ARRAY_IN
-#undef K_SIZE
-#undef MASK_ARRAY
+      ptmp = REAL(ctmp, wp)
+#  endif
+      !
+   END FUNCTION TYPENAME/**/_sum/**/ISVEC/**/_/**/XD
 !
 # endif
+/**/
+/*----------------------------------*/
+/*   FUNCTION FUNCTION_GLOBMINMAX   */
+/*----------------------------------*/
+/**/
 #if defined GLOBMINMAX_CODE
-!                          ! FUNCTION FUNCTION_GLOBMINMAX !
-#   if defined DIM_2d
-#      define XD                2d
-#      define ARRAY_TYPE(i,j,k)    REAL(wp)                 , INTENT(in   ) ::   ARRAY_IN(i,j,k)
-#      define ARRAY_IN(i,j,k)   ptab(i,j)
-#      define K_SIZE(ptab)      1
-#   endif
-#   if defined DIM_3d
-#      define XD                3d
-#      define ARRAY_TYPE(i,j,k)    REAL(wp)                 , INTENT(in   ) ::   ARRAY_IN(i,j,k)
-#      define ARRAY_IN(i,j,k)   ptab(i,j,k)
-#      define K_SIZE(ptab)      SIZE(ptab,3)
-#   endif
+/**/
+/*   DEFINE LOCAL VARIABLES   */
+/**/
 #   if defined OPERATION_GLOBMIN
-#      define OPER min
+#      define OPER      min
+#      define DEFAULT   HUGE(1._wp)
 #   endif
 #   if defined OPERATION_GLOBMAX
-#      define OPER max
+#      define OPER      max
+#      define DEFAULT   -HUGE(1._wp)
 #   endif
-
-   FUNCTION glob_/**/OPER/**/_/**/XD/**/( cdname, ptab )
+!
+#   if defined LOCALONLY
+   FUNCTION TYPENAME/**/_/**/OPER/**//**/ISVEC/**/_/**/XD/**/(         ptab ) RESULT( ptmp )
       !!----------------------------------------------------------------------
-      CHARACTER(len=*),  INTENT(in   ) ::   cdname  ! name of the calling subroutine
-      ARRAY_TYPE(:,:,:)                             ! array on which operation is applied
-      REAL(wp)   ::  glob_/**/OPER/**/_/**/XD
-      !
-      !!-----------------------------------------------------------------------
+#   else
+   FUNCTION TYPENAME/**/_/**/OPER/**//**/ISVEC/**/_/**/XD/**/( cdname, ptab ) RESULT( ptmp )
+      !!----------------------------------------------------------------------
+      CHARACTER(len=*), INTENT(in   ) ::   cdname  ! name of the calling subroutine
+#   endif
+      REAL(wp)        , INTENT(in   ) ::   ARRAY_IN(:,:,:,:)   ! array on which operation is applied
       !
-      COMPLEX(dp)::   ctmp
-      REAL(wp)   ::   ztmp
-      INTEGER    ::   jk       ! dummy loop indices
-      INTEGER    ::   ipi, ipj, ipk        ! dimensions
+#   if defined VEC
+      REAL(wp), DIMENSION(LAST_SIZE) ::   ptmp
+#   else
+      REAL(wp)   ::  ptmp
+#   endif
+      INTEGER    ::    ji,  jj,  jk,  jl        ! dummy loop indices
+      INTEGER    ::   ipi, ipj, ipk, ipl        ! dimensions
       INTEGER    ::   iisht, ijsht
       !!-----------------------------------------------------------------------
       !
       ipi = SIZE(ptab,1)   ! 1st dimension
       ipj = SIZE(ptab,2)   ! 2nd dimension
       ipk = K_SIZE(ptab)   ! 3rd dimension
+      ipl = L_SIZE(ptab)   ! 4th dimension
       !
       iisht = ( jpi - ipi ) / 2
       ijsht = ( jpj - ipj ) / 2   ! should be the same as iisht...
       !
-      ztmp = OPER/**/val( ARRAY_IN(Nis0-iisht:Nie0-iisht,Njs0-ijsht:Nje0-ijsht,1)*tmask_i(Nis0:Nie0,Njs0:Nje0) )
-      DO jk = 2, ipk
-         ztmp = OPER/**/(ztmp, OPER/**/val( ARRAY_IN(Nis0-iisht:Nie0-iisht,Njs0-ijsht:Nje0-ijsht,jk)*tmask_i(Nis0:Nie0,Njs0:Nje0) ))
-      ENDDO
-
-      CALL mpp_/**/OPER/**/( cdname, ztmp)
-
-      glob_/**/OPER/**/_/**/XD = ztmp
-
-   END FUNCTION glob_/**/OPER/**/_/**/XD
-
+      ptmp = DEFAULT
+      !
+      DO jl = 1, ipl
+         DO jk = 1, ipk
+#   define ARRAY_LOOP   ARRAY_IN(Nis0-iisht:Nie0-iisht,Njs0-ijsht:Nje0-ijsht,jk,jl)*tmask_i(Nis0:Nie0,Njs0:Nje0) 
+#   if   defined VEC && defined DIM_3d
+            ptmp(jk) = OPER/**/( ptmp(jk), OPER/**/val( ARRAY_LOOP ) )
+#   endif
+#   if   defined VEC && defined DIM_4d
+            ptmp(jl) = OPER/**/( ptmp(jl), OPER/**/val( ARRAY_LOOP ) )
+#   endif
+#   if ! defined VEC
+            ptmp     = OPER/**/( ptmp    , OPER/**/val( ARRAY_LOOP ) )
+#   endif
+         END DO
+      END DO
+      !
+#   if ! defined LOCAL
+      CALL mpp_/**/OPER/**/( cdname, ptmp )
+#   endif
+      !
+   END FUNCTION TYPENAME/**/_/**/OPER/**//**/ISVEC/**/_/**/XD
+!
+#   undef DEFAULT
+#   undef OPER
+# endif
+/**/
+/*                               */
+/*   UNDEFINE COMMON VARIABLES   */
+/*                               */
+/**/
 #undef XD
-#undef ARRAY_TYPE
 #undef ARRAY_IN
+#   if ! defined DIM_1d  
 #undef K_SIZE
-#undef OPER
-# endif
+#undef L_SIZE
+#undef LAST_SIZE
+#   endif
+#undef ISVEC
+#undef TYPENAME
+#undef ARRAY_LOOP
-- 
GitLab