From 015bd30e4dcdefc7125a610c4d5fbad375ffe2cc Mon Sep 17 00:00:00 2001
From: Philippe WAUTELET <philippe.wautelet@cnrs.fr>
Date: Wed, 9 Oct 2024 14:31:38 +0200
Subject: [PATCH] Philippe 09/10/2024: ibm_prep_ls: secure array sizes (check
 if too small and enlarge them if necessary)

---
 src/MNH/ibm_prep_ls.f90 | 93 +++++++++++++++++++++++++++--------------
 1 file changed, 62 insertions(+), 31 deletions(-)

diff --git a/src/MNH/ibm_prep_ls.f90 b/src/MNH/ibm_prep_ls.f90
index f4929d7c3..b79914c16 100644
--- a/src/MNH/ibm_prep_ls.f90
+++ b/src/MNH/ibm_prep_ls.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 2021-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 2021-2024 CNRS, Meteo-France and Universite Paul Sabatier
 !MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
 !MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
 !MNH_LIC for details. version 1.
@@ -82,6 +82,7 @@ SUBROUTINE IBM_PREP_LS(OIBM,HIBM_TYPE,PPHI)
   !    -------------
   !      Original         01/06/2021
   !  P. Wautelet 23/07/2021: replace non-standard FLOAT function by REAL function
+  !  P. Wautelet 09/10/2024: secure array sizes (check if too small and enlarge them if necessary)
   !------------------------------------------------------------------------------
   !       
   !**** 0. DECLARATIONS
@@ -131,6 +132,9 @@ SUBROUTINE IBM_PREP_LS(OIBM,HIBM_TYPE,PPHI)
   !
   !       0.2  declaration of local variables
   !
+  INTEGER, PARAMETER :: ISIZE_INIT_XYZ   = 1000000
+  INTEGER, PARAMETER :: ISIZE_INIT_FACES = 1000000
+  !
   INTEGER :: JN,JM,JNM,JL,JMM,JI,JJ,JK,JF,JV                             ! loop index
   INTEGER :: IIU,IJU
   REAL    :: ZX_MIN,ZX_MAX,ZY_MIN,ZY_MAX,DX_LOW,DY_LOW,DX_HIGH,DY_HIGH
@@ -141,10 +145,10 @@ SUBROUTINE IBM_PREP_LS(OIBM,HIBM_TYPE,PPHI)
   INTEGER :: IIBM_NUMB_SURF                                              ! number of surfaces in each type 
   REAL    :: ZIBM_X1,ZIBM_X2,ZIBM_Y1,ZIBM_Y2,ZIBM_Z1,ZIBM_Z2             ! location of surface points for one object
   REAL    :: ZIBM_TYPE_SURF
-  REAL, DIMENSION(:,:), ALLOCATABLE :: ZIBM_XYZ1,ZIBM_XYZ2               ! location of surface points for all object
-  REAL, DIMENSION(:,:), ALLOCATABLE :: ZV1,ZV1_2,ZV2,ZV2_2,ZV3,ZV3_2     ! face vectors
-  REAL, DIMENSION(:,:), ALLOCATABLE :: NORM_FACES,NORM_FACES2                 ! norm of the faces
-  REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZIBM_FACES,ZIBM_FACES2,ZIBM_FACES2b  ! extremities of triangle faces for all object
+  REAL, DIMENSION(:,:),   ALLOCATABLE :: ZIBM_XYZ1,ZIBM_XYZ2             ! location of surface points for all object
+  REAL, DIMENSION(:,:),   ALLOCATABLE :: ZV1, ZV2, ZV3                   ! face vectors
+  REAL, DIMENSION(:,:),   ALLOCATABLE :: ZNORM_FACES                     ! norm of the faces
+  REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZIBM_FACES, ZIBM_FACES2b        ! extremities of triangle faces for all object
   TYPE(LIST_ll), POINTER :: TZFIELDS_ll
   INTEGER                :: IINFO_ll
   CHARACTER(LEN=12)      :: HFILEGENE, HFILEIDEA
@@ -197,13 +201,13 @@ SUBROUTINE IBM_PREP_LS(OIBM,HIBM_TYPE,PPHI)
      !
      !Allocate the tables containing the vortices, the faces locations,
      !the norms, which are needed to calculate the LSF
-     ALLOCATE(ZIBM_XYZ1(5400000,3))
-     ALLOCATE(ZIBM_FACES(3150000,3,3))
-     ALLOCATE(ZIBM_FACES2b(3150000,3,3))
-     ALLOCATE(NORM_FACES(3150000,3))
-     ALLOCATE(ZV1(3150000,3))
-     ALLOCATE(ZV2(3150000,3))
-     ALLOCATE(ZV3(3150000,3))
+     ALLOCATE( ZIBM_XYZ1   (ISIZE_INIT_XYZ,   3 )   )
+     ALLOCATE( ZIBM_FACES  (ISIZE_INIT_FACES, 3, 3) )
+     ALLOCATE( ZIBM_FACES2b(ISIZE_INIT_FACES, 3, 3) )
+     ALLOCATE( ZNORM_FACES (ISIZE_INIT_FACES, 3)    )
+     ALLOCATE( ZV1         (ISIZE_INIT_FACES, 3)    )
+     ALLOCATE( ZV2         (ISIZE_INIT_FACES, 3)    )
+     ALLOCATE( ZV3         (ISIZE_INIT_FACES, 3)    )
      !
      OPEN(NEWUNIT=ILUIBMGENE , FILE=HFILEGENE , IOSTAT=IRESPIBMGENE , STATUS='OLD')
      !
@@ -237,6 +241,8 @@ SUBROUTINE IBM_PREP_LS(OIBM,HIBM_TYPE,PPHI)
            ZIBM_XYZ1(JV,1) = ZIBM_XYZ1(JV,1) +200.
            ZIBM_XYZ1(JV,2) = ZIBM_XYZ1(JV,2) +200.
            JV=JV+1
+           ! Increase size of array if necessary
+           IF ( JV > SIZE(ZIBM_XYZ1,1) ) CALL ENLARGE_2D( ZIBM_XYZ1 )
         ENDIF
         IF (JN==1.AND.TRIM(YSTRING(1:2))=='f') THEN
            NS2=INDEX(TRIM(YSTRING)," ",back=.true.)
@@ -253,6 +259,8 @@ SUBROUTINE IBM_PREP_LS(OIBM,HIBM_TYPE,PPHI)
            IF (ZIBM_XYZ1(ZN1,2)<ZY_MIN.AND.ZIBM_XYZ1(ZN2,2)<ZY_MIN.AND.ZIBM_XYZ1(ZN3,2)<ZY_MIN) CYCLE
            IF (ZIBM_XYZ1(ZN1,2)>ZY_MAX.AND.ZIBM_XYZ1(ZN2,2)>ZY_MAX.AND.ZIBM_XYZ1(ZN3,2)>ZY_MAX) CYCLE
            JF=JF+1
+           ! Increase size of array if necessary
+           IF ( JF > SIZE(ZIBM_FACES,1) ) CALL ENLARGE_3D( ZIBM_FACES )
            !
            ZIBM_FACES(JF,1,1) = ZIBM_XYZ1(ZN1,1)
            ZIBM_FACES(JF,1,2) = ZIBM_XYZ1(ZN1,2)
@@ -275,8 +283,18 @@ SUBROUTINE IBM_PREP_LS(OIBM,HIBM_TYPE,PPHI)
            IF (ZNB(1)==0..AND.ZNB(2)==0..AND.ZNB(3)==0.) CYCLE
            IF (ZNA(2)==0..AND.ZNA(3)==0..AND.ZNB(2)==0..AND.ZNB(3)==0.) CYCLE
            IF (ZNA(1)==ZNB(1).AND.ZNA(2)==ZNB(2).AND.ZNA(3)==ZNB(3)) CYCLE
+           !
            JCOUNT=JCOUNT+1
-           NORM_FACES(JCOUNT,:)= CROSSPRODUCT(ZNA,ZNB)
+           ! Increase size of arrays if necessary
+           IF ( JCOUNT > SIZE(ZV1,1) ) THEN
+            CALL ENLARGE_2D( ZNORM_FACES )
+            CALL ENLARGE_3D( ZIBM_FACES2b )
+            CALL ENLARGE_2D( ZV1 )
+            CALL ENLARGE_2D( ZV2 )
+            CALL ENLARGE_2D( ZV3 )
+         END IF
+           !
+           ZNORM_FACES(JCOUNT,:)= CROSSPRODUCT(ZNA,ZNB)
            ZIBM_FACES2b(JCOUNT,:,:)=ZIBM_FACES(JF,:,:)
            !
            !Equation (6) of Jones (1995)
@@ -337,24 +355,11 @@ SUBROUTINE IBM_PREP_LS(OIBM,HIBM_TYPE,PPHI)
         !
      END DO
      !
-     ALLOCATE(ZIBM_FACES2(JCOUNT,3,3))
-     ALLOCATE(NORM_FACES2(JCOUNT,3))
-     ALLOCATE(ZV1_2(JCOUNT,3))
-     ALLOCATE(ZV2_2(JCOUNT,3))
-     ALLOCATE(ZV3_2(JCOUNT,3))
-     !
-     NORM_FACES2 = NORM_FACES(:JCOUNT,:)
-     ZV1_2 = ZV1(:JCOUNT,:)
-     ZV2_2 = ZV2(:JCOUNT,:)
-     ZV3_2 = ZV3(:JCOUNT,:)
-     ZIBM_FACES2 = ZIBM_FACES2b(:JCOUNT,:,:)
-     !
-  ENDIF
   !
   !----------------------------
   !---------Idealized case-----
   !----------------------------
-  IF (HIBM_TYPE=='IDEA') THEN
+  ELSE IF (HIBM_TYPE=='IDEA') THEN
      !
      OPEN(NEWUNIT=ILUIBMIDEA , FILE= HFILEIDEA , IOSTAT=IRESPIBMIDEA , FORM='FORMATTED' , &
           STATUS='OLD', ACCESS='SEQUENTIAL', ACTION='READ')
@@ -398,10 +403,9 @@ SUBROUTINE IBM_PREP_LS(OIBM,HIBM_TYPE,PPHI)
   !                            locations of interface (ellipsoidal/parallelepipedic shapes)
   !
   IF (HIBM_TYPE=='GENE') THEN
-     CALL IBM_GENERLS(ZIBM_FACES2,NORM_FACES2,ZV1_2,ZV2_2,ZV3_2,ZX_MIN,ZY_MIN,ZX_MAX,ZY_MAX,PPHI)
-  ENDIF
-  !
-  IF (HIBM_TYPE=='IDEA') then
+     CALL IBM_GENERLS(ZIBM_FACES2b(:JCOUNT,:,:),ZNORM_FACES(:JCOUNT,:),ZV1(:JCOUNT,:),ZV2(:JCOUNT,:),ZV3(:JCOUNT,:), &
+                      ZX_MIN,ZY_MIN,ZX_MAX,ZY_MAX,PPHI)
+  ELSEIF (HIBM_TYPE=='IDEA') then
      DO JN=1,JNM
         IF (abs(ZIBM_XYZ2(JN,7)-1.).lt.XIBM_EPSI) CALL IBM_IDEALRP(JN,ZIBM_XYZ2,PPHI)
         IF (abs(ZIBM_XYZ2(JN,7)-2.).lt.XIBM_EPSI) CALL IBM_IDEALEE(JN,ZIBM_XYZ2,PPHI)
@@ -428,4 +432,31 @@ CONTAINS
     !
   END FUNCTION CROSSPRODUCT
   !
+  SUBROUTINE ENLARGE_2D( PA )
+    REAL, DIMENSION(:,:), ALLOCATABLE, INTENT(INOUT) :: PA
+    !
+    INTEGER, PARAMETER :: ISTEP = 1000000 ! step of size increment
+    !
+    REAL, DIMENSION(:,:), ALLOCATABLE :: ZTMP2D ! temporary array
+    !
+    ALLOCATE( ZTMP2D( SIZE(PA,1)+ISTEP, SIZE(PA,2) ) )
+    !
+    ZTMP2D(1:SIZE(PA,1)+ISTEP,:) = PA(:,:)
+    !
+    CALL MOVE_ALLOC( FROM = ZTMP2D, TO = PA )
+   END SUBROUTINE ENLARGE_2D
+  !
+  SUBROUTINE ENLARGE_3D( PA )
+    REAL, DIMENSION(:,:,:), ALLOCATABLE, INTENT(INOUT) :: PA
+    !
+    INTEGER, PARAMETER :: ISTEP = 1000000 ! step of size increment
+    !
+    REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZTMP3D ! temporary array
+    !
+    ALLOCATE( ZTMP3D( SIZE(PA,1)+ISTEP, SIZE(PA,2), SIZE(PA,3) ) )
+    !
+    ZTMP3D(1:SIZE(PA,1)+ISTEP,:,:) = PA(:,:,:)
+    !
+    CALL MOVE_ALLOC( FROM = ZTMP3D, TO = PA )
+   END SUBROUTINE ENLARGE_3D
 END SUBROUTINE IBM_PREP_LS
-- 
GitLab