From c280a48e83c29e12f7d0ebb2bc3365b78c84507a Mon Sep 17 00:00:00 2001
From: Philippe WAUTELET <philippe.wautelet@aero.obs-mip.fr>
Date: Wed, 7 Sep 2022 10:39:59 +0200
Subject: [PATCH] Philippe 07/09/2022: correction for array boundaries +
 cleaning

---
 src/MNH/ch_init_fieldn.f90 | 337 +++++++++++++++++++------------------
 1 file changed, 170 insertions(+), 167 deletions(-)

diff --git a/src/MNH/ch_init_fieldn.f90 b/src/MNH/ch_init_fieldn.f90
index 4c9853b05..28331efcd 100644
--- a/src/MNH/ch_init_fieldn.f90
+++ b/src/MNH/ch_init_fieldn.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1995-2019 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1995-2022 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.
@@ -71,46 +71,35 @@ END MODULE MODI_CH_INIT_FIELD_n
 !!    04/06/07 (M. Leriche & JP Pinty) add pH initialization
 !!    20/04/10 (M. Leriche) remove pH initialization to ini_modeln
 !  P. Wautelet 20/05/2019: add name argument to ADDnFIELD_ll + new ADD4DFIELD_ll subroutine
-!  P. Tulet  20/05/2021:  correction for CON to MIX transformation unit (aerosols only)
-!!
-!!    EXTERNAL
-!!    --------
-!!     GET_DIM_EXT_ll : get extended sub-domain sizes
-!!     GET_INDICE_ll  : get physical sub-domain bounds
-!!
-!!------------------------------------------------------------------------------
-!!
+!  P. Tulet    20/05/2021: correction for CON to MIX transformation unit (aerosols only)
+!  P. Wautelet 07/09/2022: correction for array boundaries + cleaning
+!------------------------------------------------------------------------------
+!
+USE MODD_ARGSLIST_ll, ONLY: LIST_ll  ! for update_halo
+USE MODD_CONF,        ONLY: CPROGRAM, L1D, L2D
+USE MODD_CH_AERO_n
+USE MODD_CH_AEROSOL
+USE MODD_CH_M9_n,     ONLY: CNAMES, NEQ
+USE MODD_CH_MNHC_n,   ONLY: CCH_SCHEME, LCH_INIT_FIELD
+USE MODD_CST,         ONLY: XMD, XAVOGADRO
+USE MODD_FIELD_n,     ONLY: XSVT       ! scalar variable at t
+USE MODD_GRID_n,      ONLY: XZZ   ! height z
+USE MODD_LSFIELD_n,   ONLY: XLBXSVM, XLBYSVM
+USE MODD_NSV,         ONLY: NSV_CHEM, NSV_CHEMBEG,NSV_CHEMEND, &
+                            NSV_AER, NSV_AERBEG,NSV_AEREND
+USE MODD_PARAMETERS,  ONLY: JPVEXT, JPHEXT  ! number of External points
+USE MODD_REF_n,       ONLY: XRHODREF ! dry density of ref. state
+
+USE MODE_AERO_PSD,    ONLY: CON2MIX
+USE MODE_ARGSLIST_ll, ONLY: ADD4DFIELD_ll, CLEANLIST_ll
+USE MODE_EXCHANGE_ll, ONLY: UPDATE_HALO_ll
+USE MODE_MSG
+USE MODE_TOOLS_ll, ONLY: GET_INDICE_ll, LWEST_ll, LEAST_ll, LNORTH_ll, LSOUTH_ll
+
+USE MODI_CH_AER_EQM_INIT_n
 USE MODI_CH_FIELD_VALUE_n ! returns value of chemical species at each grid point
 USE MODI_CH_INIT_CONST_n
-USE MODI_CH_AER_EQM_INIT_n
-USE MODE_ll
-USE MODE_AERO_PSD
-!!
-!!    IMPLICIT ARGUMENTS
-!!    ------------------
-USE MODD_GRID_n,      ONLY : XZZ,       &! height z
-                             XLAT,XLON   ! latitude and longitude
-USE MODD_REF_n,       ONLY : XRHODREF,  &! dry density of ref. state
-                             XRHODJ      ! ( rhod J ) = dry density
-USE MODD_LBC_n
-USE MODD_NSV,         ONLY : NSV_CHEM, NSV_CHEMBEG,NSV_CHEMEND, &
-                             NSV_AER, NSV_AERBEG,NSV_AEREND
-USE MODD_CST,         ONLY : XMD, XAVOGADRO
 
-USE MODD_FIELD_n,     ONLY : XSVT       ! scalar variable at t
-USE MODD_PARAMETERS,  ONLY : JPVEXT, JPHEXT  ! number of External points
-USE MODD_ARGSLIST_ll, ONLY : LIST_ll  ! for update_halo
-USE MODD_CH_CONST_n                   ! for Chemical constants
-USE MODD_CONF,        ONLY : CPROGRAM, L1D, L2D 
-USE MODD_CONF_n,      ONLY : NRRL
-USE MODD_CH_MNHC_n
-USE MODD_CH_M9_n,     ONLY : CNAMES, NEQ
-USE MODD_CH_AEROSOL
-USE MODD_CH_AERO_n
-USE MODD_LSFIELD_n,   ONLY : XLBXSVM, XLBYSVM
-USE MODD_DYN_n,       ONLY : NRIMX,NRIMY
-!!
-!!
 !-------------------------------------------------------------------------------
 !
 !*       0.     DECLARATIONS
@@ -153,7 +142,7 @@ INTEGER             :: IKE  ! indice K End       in z direction
 !
 TYPE(LIST_ll), POINTER :: TZFIELDS_ll ! pointer for the list of 3D fields
 INTEGER :: IINFO_ll  ! Return code of //routines
-INTEGER :: IOR, JOR, IEND, JEND, KINFO, NIU,NJU, ILBX, ILBY, IRIMX, IRIMY
+INTEGER :: ILBX, ILBY, IRIMX, IRIMY
 !
 !-------------------------------------------------------------------------------
 !
@@ -165,39 +154,28 @@ NULLIFY(TZFIELDS_ll)
 !*       1.    PREPARE INITIALIZATION
 !              ----------------------
 
-IF (CORGANIC == TRIM("MPMPO") .OR. CORGANIC == TRIM("PUN") .OR. CORGANIC == TRIM("EQSAM2")) THEN
-  IF ((CCH_SCHEME .EQ. TRIM("NONE")) .OR. (CCH_SCHEME .EQ. TRIM("RELACS"))&
-   .OR. (CCH_SCHEME .EQ. TRIM("RACM"))) THEN
-  WRITE(KLUOUT,FMT=*) '**********************************************'
-  WRITE(KLUOUT,FMT=*) 'WARNING : NO SOA !!!!'
-  WRITE(KLUOUT,FMT=*) 'YOU WANT TO USE SOA GAS PARTICLE BALANCE'
-  WRITE(KLUOUT,FMT=*) 'BUT THE SCHEME NEED TO BE CACM or RELACS 2'
-  WRITE(KLUOUT,FMT=*) 'CORGANIC HAS BEEN SET TO NONE'
-  WRITE(KLUOUT,FMT=*) 'OTHERWISE COMPILE THE CORRECT SCHEME BEFORE'
-  WRITE(KLUOUT,FMT=*) '**********************************************'
+IF ( ( CORGANIC == "MPMPO" .OR. CORGANIC == "PUN" .OR. CORGANIC == "EQSAM2" ) .AND. &
+     ( CCH_SCHEME == "NONE" .OR. CCH_SCHEME == "RELACS" .OR. CCH_SCHEME == "RACM" ) ) THEN
+  CMNHMSG(1) = 'no SOA'
+  CMNHMSG(2) = 'you want to use SOA gas particle balance'
+  CMNHMSG(3) = 'but the scheme need to be CACM or RELACS 2'
+  CMNHMSG(4) = 'CORGANIC has been set to NONE'
+  CMNHMSG(5) = 'otherwise compile the correct scheme before'
+  CALL PRINT_MSG( NVERB_WARNING, 'GEN', 'CH_INIT_FIELD_n' )
   CORGANIC = "NONE"
-  END IF
 END IF
 !
 !*       1.1   compute dimensions of arrays
 !
-
+IIU = SIZE(XSVT,1)
+IJU = SIZE(XSVT,2)
+IKU = SIZE(XSVT,3)
 CALL GET_INDICE_ll (IIB,IJB,IIE,IJE)
-CALL GET_GLOBALDIMS_ll(IIE,IJE)
-IIU = IIE + 2 * JPHEXT
-IJU = IJE + 2 * JPHEXT
-CALL GET_INTERSECTION_ll(1+JPHEXT, 1+JPHEXT, IIU-JPHEXT , IJU-JPHEXT, IOR, JOR, IEND, JEND, "EXTE", KINFO)
 IKB = 1 + JPVEXT
-IKU = SIZE(XSVT,3)
 IKE = IKU - JPVEXT
-CALL GET_DIM_EXT_ll('B',NIU,NJU)
 !
 !        1.1.1 find maximum height level
 ILEVMAX=INT(MAXVAL(XZZ(:,:,IKE)/10.))+1
-! the following print serves to break compiler optimization with MAXVAL
-! (pb. on OS2000 with option -O3 for example, Peter Bechtold had
-! similar surprises with MAXVAL on Fuji VPP700 in the convection scheme)
-WRITE(KLUOUT,*) "CH_INIT_FIELD_n: ILEVMAX =",ILEVMAX 
 ALLOCATE(ZHEIGHT(ILEVMAX))
 ALLOCATE(ZSVINIT(ILEVMAX,NEQ))
 ALLOCATE(ZSVINITA(ILEVMAX,NSV_AER))
@@ -212,44 +190,42 @@ ZDEN2MOL = 1E-6 * XAVOGADRO / XMD
 !*       2.    INITIALIZE T FIELDS AND CONVERT CONC. TO MIXING RATIO
 !              -----------------------
 !
-
 YUNIT="MIX"
 
 IF (LORILAM) THEN
-  IF (.NOT.(ASSOCIATED(XN3D)))      ALLOCATE(XN3D(SIZE(XSVT,1),SIZE(XSVT,2),IKU,JPMODE))
-  IF (.NOT.(ASSOCIATED(XRG3D)))     ALLOCATE(XRG3D(SIZE(XSVT,1),SIZE(XSVT,2),IKU,JPMODE))
-  IF (.NOT.(ASSOCIATED(XSIG3D)))    ALLOCATE(XSIG3D(SIZE(XSVT,1),SIZE(XSVT,2),IKU,JPMODE))
-  IF (.NOT.(ASSOCIATED(XRHOP3D)))   ALLOCATE(XRHOP3D(SIZE(XSVT,1),SIZE(XSVT,2),IKU,JPMODE))
-  IF (.NOT.(ASSOCIATED(XM3D)))      ALLOCATE(XM3D(SIZE(XSVT,1),SIZE(XSVT,2),IKU,JPMODE*3))
-  IF (.NOT.(ASSOCIATED(XSEDA)))     ALLOCATE(XSEDA(SIZE(XSVT,1),SIZE(XSVT,2),IKU,JPMODE*3))
-  IF (.NOT.(ASSOCIATED(XCTOTA3D))) &
-    ALLOCATE(XCTOTA3D(SIZE(XSVT,1),SIZE(XSVT,2),IKU,NSP+NCARB+NSOA,JPMODE))
-  IF (.NOT.(ASSOCIATED(XVDEPAERO))) ALLOCATE(XVDEPAERO(SIZE(XSVT,1),SIZE(XSVT,2),JPIN))
+  IF (.NOT.(ASSOCIATED(XN3D)))      ALLOCATE(XN3D     (IIU,IJU,IKU,JPMODE))
+  IF (.NOT.(ASSOCIATED(XRG3D)))     ALLOCATE(XRG3D    (IIU,IJU,IKU,JPMODE))
+  IF (.NOT.(ASSOCIATED(XSIG3D)))    ALLOCATE(XSIG3D   (IIU,IJU,IKU,JPMODE))
+  IF (.NOT.(ASSOCIATED(XRHOP3D)))   ALLOCATE(XRHOP3D  (IIU,IJU,IKU,JPMODE))
+  IF (.NOT.(ASSOCIATED(XM3D)))      ALLOCATE(XM3D     (IIU,IJU,IKU,JPMODE*3))
+  IF (.NOT.(ASSOCIATED(XSEDA)))     ALLOCATE(XSEDA    (IIU,IJU,IKU,JPMODE*3))
+  IF (.NOT.(ASSOCIATED(XCTOTA3D)))  ALLOCATE(XCTOTA3D (IIU,IJU,IKU,NSP+NCARB+NSOA,JPMODE))
+  IF (.NOT.(ASSOCIATED(XVDEPAERO))) ALLOCATE(XVDEPAERO(IIU,IJU,JPIN))
   IF (.NOT.(ALLOCATED(XFAC)))       ALLOCATE(XFAC(NSP+NSOA+NCARB))
   IF (.NOT.(ALLOCATED(XRHOI)))      ALLOCATE(XRHOI(NSP+NSOA+NCARB))
   IF (.NOT.(ASSOCIATED(XFRAC))) THEN
-    ALLOCATE(XFRAC(SIZE(XSVT,1),SIZE(XSVT,2),IKU,NEQ))
+    ALLOCATE(XFRAC(IIU,IJU,IKU,NEQ))
     XFRAC(:,:,:,:) = 0.
   END IF
   IF (.NOT.(ASSOCIATED(XMI))) THEN
-    ALLOCATE(XMI(SIZE(XSVT,1),SIZE(XSVT,2),IKU,NSP+NCARB+NSOA))
+    ALLOCATE(XMI(IIU,IJU,IKU,NSP+NCARB+NSOA))
   END IF
-  IF (.NOT.(ASSOCIATED(XJNUC)))        ALLOCATE(XJNUC(SIZE(XSVT,1),SIZE(XSVT,2),IKU))
-  IF (.NOT.(ASSOCIATED(XJ2RAT)))       ALLOCATE(XJ2RAT(SIZE(XSVT,1),SIZE(XSVT,2),IKU))
-  IF (.NOT.(ASSOCIATED(XCONC_MASS)))   ALLOCATE(XCONC_MASS(SIZE(XSVT,1),SIZE(XSVT,2),IKU))
-  IF (.NOT.(ASSOCIATED(XCOND_MASS_I))) ALLOCATE(XCOND_MASS_I(SIZE(XSVT,1),SIZE(XSVT,2),IKU))
-  IF (.NOT.(ASSOCIATED(XCOND_MASS_J))) ALLOCATE(XCOND_MASS_J(SIZE(XSVT,1),SIZE(XSVT,2),IKU))
-  IF (.NOT.(ASSOCIATED(XNUCL_MASS)))   ALLOCATE(XNUCL_MASS(SIZE(XSVT,1),SIZE(XSVT,2),IKU))
+  IF (.NOT.(ASSOCIATED(XJNUC)))        ALLOCATE(XJNUC       (IIU,IJU,IKU))
+  IF (.NOT.(ASSOCIATED(XJ2RAT)))       ALLOCATE(XJ2RAT      (IIU,IJU,IKU))
+  IF (.NOT.(ASSOCIATED(XCONC_MASS)))   ALLOCATE(XCONC_MASS  (IIU,IJU,IKU))
+  IF (.NOT.(ASSOCIATED(XCOND_MASS_I))) ALLOCATE(XCOND_MASS_I(IIU,IJU,IKU))
+  IF (.NOT.(ASSOCIATED(XCOND_MASS_J))) ALLOCATE(XCOND_MASS_J(IIU,IJU,IKU))
+  IF (.NOT.(ASSOCIATED(XNUCL_MASS)))   ALLOCATE(XNUCL_MASS  (IIU,IJU,IKU))
 
-  IF (.NOT.(ASSOCIATED(XMBEG)))   ALLOCATE(XMBEG(SIZE(XSVT,1),SIZE(XSVT,2),IKU,JPIN))
-  IF (.NOT.(ASSOCIATED(XMINT)))   ALLOCATE(XMINT(SIZE(XSVT,1),SIZE(XSVT,2),IKU,JPIN))
-  IF (.NOT.(ASSOCIATED(XMEND)))   ALLOCATE(XMEND(SIZE(XSVT,1),SIZE(XSVT,2),IKU,JPIN))
+  IF (.NOT.(ASSOCIATED(XMBEG)))   ALLOCATE(XMBEG(IIU,IJU,IKU,JPIN))
+  IF (.NOT.(ASSOCIATED(XMINT)))   ALLOCATE(XMINT(IIU,IJU,IKU,JPIN))
+  IF (.NOT.(ASSOCIATED(XMEND)))   ALLOCATE(XMEND(IIU,IJU,IKU,JPIN))
 
-  IF (.NOT.(ASSOCIATED(XDMINTRA)))  ALLOCATE(XDMINTRA(SIZE(XSVT,1),SIZE(XSVT,2),IKU,JPIN))
-  IF (.NOT.(ASSOCIATED(XDMINTER)))  ALLOCATE(XDMINTER(SIZE(XSVT,1),SIZE(XSVT,2),IKU,JPIN))
-  IF (.NOT.(ASSOCIATED(XDMCOND)))   ALLOCATE(XDMCOND(SIZE(XSVT,1),SIZE(XSVT,2),IKU,JPIN))
-  IF (.NOT.(ASSOCIATED(XDMNUCL)))   ALLOCATE(XDMNUCL(SIZE(XSVT,1),SIZE(XSVT,2),IKU,JPIN))
-  IF (.NOT.(ASSOCIATED(XDMMERG)))   ALLOCATE(XDMMERG(SIZE(XSVT,1),SIZE(XSVT,2),IKU,JPIN))
+  IF (.NOT.(ASSOCIATED(XDMINTRA)))  ALLOCATE(XDMINTRA(IIU,IJU,IKU,JPIN))
+  IF (.NOT.(ASSOCIATED(XDMINTER)))  ALLOCATE(XDMINTER(IIU,IJU,IKU,JPIN))
+  IF (.NOT.(ASSOCIATED(XDMCOND)))   ALLOCATE(XDMCOND (IIU,IJU,IKU,JPIN))
+  IF (.NOT.(ASSOCIATED(XDMNUCL)))   ALLOCATE(XDMNUCL (IIU,IJU,IKU,JPIN))
+  IF (.NOT.(ASSOCIATED(XDMMERG)))   ALLOCATE(XDMMERG (IIU,IJU,IKU,JPIN))
   !
   XJNUC(:,:,:)        = 1.0E-7
   XJ2RAT(:,:,:)       = 0.
@@ -267,18 +243,12 @@ IF (LORILAM) THEN
   XDMCOND(:,:,:,:)      = 0.
   XDMNUCL(:,:,:,:)      = 0.
   XDMMERG(:,:,:,:)      = 0.
-
 END IF
 !
-!*          print info for user
-IF ((LCH_INIT_FIELD).AND.(CPROGRAM/='DIAG  ')) THEN
-!
-  WRITE(KLUOUT,*) "CH_INIT_FIELD_n will now initialize XSVT fields"
-!
-!
-  jlev_loop : DO JLEV=1,ILEVMAX
+IF ( LCH_INIT_FIELD .AND. CPROGRAM /= 'DIAG  ' ) THEN
+  DO JLEV=1,ILEVMAX
     ZHEIGHT=REAL(JLEV-1)*10.
-    jn_loop : DO JN = 1, NEQ
+    DO JN = 1, NEQ
       ZSVINIT(JLEV,JN) = &
            CH_FIELD_VALUE_n(ZHEIGHT(JLEV), "LLZ", &
            CNAMES(JN), YUNIT, KLUOUT, KVERB)
@@ -286,92 +256,123 @@ IF ((LCH_INIT_FIELD).AND.(CPROGRAM/='DIAG  ')) THEN
       ! CH_FIELD_VALUE_n ("LLZ"=lon-lat-Z)
       ! in future developpements, "IJK" may be used in order
       ! to pass the grid indices rather than coordinates
-    END DO jn_loop
-  END DO jlev_loop
-
-  XSVT(:,:,:,NSV_CHEMBEG:NSV_CHEMEND) = 0.
-  jk_loop : DO JK = IKB, IKE
-    jj_loop : DO JJ = JOR, JEND
-      ji_loop : DO JI = IOR, IEND
+    END DO
+  END DO
 
+  DO JK = IKB, IKE
+    DO JJ = IJB, IJE
+      DO JI = IIB, IIE
         JLEV=INT(MAX(XZZ(JI,JJ,JK),0.)/10.)+1
         XSVT(JI,JJ,JK,NSV_CHEMBEG:NSV_CHEMEND) = ZSVINIT(JLEV,:)
+      END DO
+    END DO
+  END DO
 
-      END DO ji_loop
-    END DO jj_loop
-  END DO jk_loop
-  DO JN = NSV_CHEMBEG,NSV_CHEMEND
-    DO JK=1,JPVEXT
-      XSVT(:,:,IKB-JPVEXT,JN) = XSVT(:,:,IKB,JN)
-      XSVT(:,:,IKE+JPVEXT,JN) = XSVT(:,:,IKE,JN)
+  DO JN = NSV_CHEMBEG, NSV_CHEMEND
+    DO JK = 1, JPVEXT
+      XSVT(:,:,IKB-JK,JN) = XSVT(:,:,IKB,JN)
+      XSVT(:,:,IKE+JK,JN) = XSVT(:,:,IKE,JN)
+    END DO
 
-      XSVT(IIB-JPHEXT,:,:,JN) = XSVT(IIB,:,:,JN)
-      XSVT(IIU,:,:,JN) = XSVT(IIU-JPHEXT,:,:,JN)
+    IF ( LWEST_ll() )  THEN
+      DO JI = 1, JPHEXT
+        XSVT(JI,:,:,JN) = XSVT(IIB,:,:,JN)
+      END DO
+    END IF
 
-      XSVT(:,IJB-JPHEXT,:,JN) = XSVT(:,IJB,:,JN)
-      XSVT(:,IJU,:,JN) = XSVT(:,IJU-JPHEXT,:,JN)
-    END DO
+    IF ( LEAST_ll() )  THEN
+      DO JI = 1, JPHEXT
+        XSVT(IIE+JI,:,:,JN) = XSVT(IIE,:,:,JN)
+      END DO
+    END IF
+
+    IF ( LSOUTH_ll() ) THEN
+      DO JJ = 1, JPHEXT
+        XSVT(:,JJ,:,JN) = XSVT(:,IJB,:,JN)
+      END DO
+    END IF
+
+    IF ( LNORTH_ll() ) THEN
+      DO JJ = 1, JPHEXT
+        XSVT(:,IJE+JJ,:,JN) = XSVT(:,IJE,:,JN)
+      END DO
+    END IF
   END DO
   !
-  IF (YUNIT .EQ. "CON") THEN
-    WRITE(KLUOUT,*) "CH_INIT_FIELD_n: converting initial values to mixing ratio"
+  IF ( YUNIT == "CON" ) THEN
+    CALL PRINT_MSG( NVERB_INFO, 'GEN', 'CH_INIT_FIELD_n', 'converting initial values to mixing ratio' )
     DO JN = NSV_CHEMBEG,NSV_CHEMEND
       XSVT(:,:,:,JN) = XSVT(:,:,:,JN)/(XRHODREF(:,:,:)*ZDEN2MOL)
     ENDDO
   ELSE
-    WRITE(KLUOUT,*)"CH_INIT_FIELD_n: initial values are used as is (mixing ratio)"
+    CALL PRINT_MSG( NVERB_INFO, 'GEN', 'CH_INIT_FIELD_n', 'initial values are used as is (mixing ratio)' )
   ENDIF
 !
 !
-  IF (LORILAM) THEN
-  jlev_loop2 : DO JLEV=1,ILEVMAX
-    ZHEIGHT=REAL(JLEV-1)*10.
-    jn_loop2 : DO JN = 1, NSV_AER
-      ZSVINITA(JLEV,JN) = &
-           CH_FIELD_VALUE_n(ZHEIGHT(JLEV), "LLZ", &
-           CAERONAMES(JN), YUNIT, KLUOUT, KVERB)
-      ! "LLZ" identifies the type of x-y-z values passed on to
-      ! CH_FIELD_VALUE_n ("LLZ"=lon-lat-Z)
-      ! in future developpements, "IJK" may be used in order
-      ! to pass the grid indices rather than coordinates
-    END DO jn_loop2
-  END DO jlev_loop2
-  !
-  XSVT(:,:,:,NSV_AERBEG:NSV_AEREND) = 0.
-  jk_loop2 : DO JK = IKB, IKE
-    jj_loop2 : DO JJ = JOR, JEND
-      ji_loop2 : DO JI = IOR, IEND
+  ORILAM: IF (LORILAM) THEN
+    DO JLEV=1,ILEVMAX
+      ZHEIGHT=REAL(JLEV-1)*10.
+      DO JN = 1, NSV_AER
+        ZSVINITA(JLEV,JN) = &
+             CH_FIELD_VALUE_n(ZHEIGHT(JLEV), "LLZ", &
+             CAERONAMES(JN), YUNIT, KLUOUT, KVERB)
+        ! "LLZ" identifies the type of x-y-z values passed on to
+        ! CH_FIELD_VALUE_n ("LLZ"=lon-lat-Z)
+        ! in future developpements, "IJK" may be used in order
+        ! to pass the grid indices rather than coordinates
+      END DO
+    END DO
+    !
+    DO JK = IKB, IKE
+      DO JJ = IJB, IJE
+        DO JI = IIB, IIE
+          JLEV=INT(MAX(XZZ(JI,JJ,JK),0.)/10.)+1
+          XSVT(JI,JJ,JK,NSV_AERBEG:NSV_AEREND) = ZSVINITA(JLEV,:)
+        END DO
+      END DO
+    END DO
 
-        JLEV=INT(MAX(XZZ(JI,JJ,JK),0.)/10.)+1
-        XSVT(JI,JJ,JK,NSV_AERBEG:NSV_AEREND) = ZSVINITA(JLEV,:)
+    DO JN = NSV_AERBEG,NSV_AEREND
+      DO JK = 1, JPVEXT
+        XSVT(:,:,IKB-JK,JN) = XSVT(:,:,IKB,JN)
+        XSVT(:,:,IKE+JK,JN) = XSVT(:,:,IKE,JN)
+      END DO
 
-      END DO ji_loop2
-    END DO jj_loop2
-  END DO jk_loop2
-  DO JN = NSV_AERBEG,NSV_AEREND
-    DO JK=1,JPVEXT
-      XSVT(:,:,IKB-JPVEXT,JN) = XSVT(:,:,IKB,JN)
-      XSVT(:,:,IKE+JPVEXT,JN) = XSVT(:,:,IKE,JN)
+      IF ( LWEST_ll() )  THEN
+        DO JI = 1, JPHEXT
+          XSVT(JI,:,:,JN) = XSVT(IIB,:,:,JN)
+        END DO
+      END IF
 
-      XSVT(IIB-JPHEXT,:,:,JN) = XSVT(IIB,:,:,JN)
-      XSVT(IIU,:,:,JN) = XSVT(IIU-JPHEXT,:,:,JN)
+      IF ( LEAST_ll() )  THEN
+        DO JI = 1, JPHEXT
+          XSVT(IIE+JI,:,:,JN) = XSVT(IIE,:,:,JN)
+        END DO
+      END IF
 
-      XSVT(:,IJB-JPHEXT,:,JN) = XSVT(:,IJB,:,JN)
-      XSVT(:,IJU,:,JN) = XSVT(:,IJU-JPHEXT,:,JN)
-    END DO
-  END DO
-  !
-  IF (YUNIT .EQ. "CON") THEN
-    WRITE(KLUOUT,*) "CH_INIT_FIELD_n (ORILAM): converting initial values µg/m3 to mixing ratio"
-    CALL CON2MIX (XSVT(:,:,:,NSV_AERBEG:NSV_AEREND), XRHODREF)
-  ELSE
-    WRITE(KLUOUT,*)"CH_INIT_FIELD_n (ORILAM): initial values are used as is (mixing ratio)"
-  ENDIF
+      IF ( LSOUTH_ll() ) THEN
+        DO JJ = 1, JPHEXT
+          XSVT(:,JJ,:,JN) = XSVT(:,IJB,:,JN)
+        END DO
+      END IF
 
+      IF ( LNORTH_ll() ) THEN
+        DO JJ = 1, JPHEXT
+          XSVT(:,IJE+JJ,:,JN) = XSVT(:,IJE,:,JN)
+        END DO
+      END IF
+    END DO
+    !
+    IF ( YUNIT == "CON" ) THEN
+      CALL PRINT_MSG( NVERB_INFO, 'GEN', 'CH_INIT_FIELD_n', 'ORILAM: converting initial values ug/m3 to mixing ratio' )
+      CALL CON2MIX (XSVT(:,:,:,NSV_AERBEG:NSV_AEREND), XRHODREF)
+    ELSE
+      CALL PRINT_MSG( NVERB_INFO, 'GEN', 'CH_INIT_FIELD_n', 'ORILAM: initial values are used as is (mixing ratio)' )
+    ENDIF
+    !
+  END IF ORILAM
   !
-  ENDIF !LORILAM
-  !
-ENDIF
+END IF
 !
 !
 CALL ADD4DFIELD_ll(TZFIELDS_ll, XSVT(:,:,:,NSV_CHEMBEG:NSV_CHEMEND), 'CH_INIT_FIELD_n::XSVT(:,:,:,NSV_CHEMBEG:NSV_CHEMEND)' )
@@ -397,8 +398,10 @@ IF (LORILAM) THEN
                          XM3D,XRHOP3D,XSIG3D,& 
                          XRG3D,XN3D, XRHODREF, XCTOTA3D) 
   DO JN = 1,JPIN
-      XM3D(:,:,IKB-JPVEXT,JN) =  XM3D(:,:,IKB,JN)
-      XM3D(:,:,IKE+JPVEXT,JN) =  XM3D(:,:,IKE,JN)
+    DO JK = 1, JPVEXT
+      XM3D(:,:,IKB-JK,JN) =  XM3D(:,:,IKB,JN)
+      XM3D(:,:,IKE+JK,JN) =  XM3D(:,:,IKE,JN)
+    END DO
   END DO
   !
   CALL ADD4DFIELD_ll(TZFIELDS_ll, XSVT(:,:,:,NSV_CHEMBEG:NSV_CHEMEND), 'CH_INIT_FIELD_n::XSVT(:,:,:,NSV_CHEMBEG,NSV_CHEMEND)' )
@@ -422,22 +425,22 @@ IF ((LCH_INIT_FIELD).AND.(CPROGRAM/='DIAG  ').AND.(KMI .EQ. 1)) THEN
     IF(LWEST_ll() .AND. .NOT. L1D) &
       XLBXSVM(1:IRIMX+1,        :,:,JN) = XSVT(1:IRIMX+1,        :,:,JN)
     IF(LEAST_ll() .AND. .NOT. L1D) &
-      XLBXSVM(ILBX-IRIMX:ILBX,:,:,JN)   = XSVT(NIU-IRIMX:NIU,    :,:,JN)
+      XLBXSVM(ILBX-IRIMX:ILBX,:,:,JN)   = XSVT(IIU-IRIMX:IIU,    :,:,JN)
     IF(LSOUTH_ll() .AND. .NOT. L1D .AND. .NOT. L2D) &
       XLBYSVM(:,1:IRIMY+1,        :,JN) = XSVT(:,1:IRIMY+1,      :,JN)
     IF(LNORTH_ll() .AND. .NOT. L1D .AND. .NOT. L2D) &
-      XLBYSVM(:,ILBY-IRIMY:ILBY,:,JN)   = XSVT(:,NJU-IRIMY:NJU,  :,JN)
+      XLBYSVM(:,ILBY-IRIMY:ILBY,:,JN)   = XSVT(:,IJU-IRIMY:IJU,  :,JN)
   END DO
   IF (LORILAM) THEN
   DO JN = NSV_AERBEG,NSV_AEREND
     IF(LWEST_ll() .AND. .NOT. L1D) &
       XLBXSVM(1:IRIMX+1,        :,:,JN) = XSVT(1:IRIMX+1,        :,:,JN)
     IF(LEAST_ll() .AND. .NOT. L1D) &
-      XLBXSVM(ILBX-IRIMX:ILBX,:,:,JN)   = XSVT(NIU-IRIMX:NIU,    :,:,JN)
+      XLBXSVM(ILBX-IRIMX:ILBX,:,:,JN)   = XSVT(IIU-IRIMX:IIU,    :,:,JN)
     IF(LSOUTH_ll() .AND. .NOT. L1D .AND. .NOT. L2D) &
       XLBYSVM(:,1:IRIMY+1,        :,JN) = XSVT(:,1:IRIMY+1,      :,JN)
     IF(LNORTH_ll() .AND. .NOT. L1D .AND. .NOT. L2D) &
-      XLBYSVM(:,ILBY-IRIMY:ILBY,:,JN)   = XSVT(:,NJU-IRIMY:NJU,  :,JN)
+      XLBYSVM(:,ILBY-IRIMY:ILBY,:,JN)   = XSVT(:,IJU-IRIMY:IJU,  :,JN)
   END DO
   ENDIF
 !
-- 
GitLab