From 6210b343a7ebeba3b1fc75c4631f5bc8afb25fba Mon Sep 17 00:00:00 2001
From: JorisP <pianezze.joris@gmail.com>
Date: Thu, 11 Jul 2024 15:25:11 +0200
Subject: [PATCH] Bugfix for aerosols (salt, dust and anthropogenic) : default
 valuers, add missing outputs, correct conversion, ...

---
 src/MNH/aer2lima.f90              | 60 +++++++++------------
 src/MNH/dustcamsn.f90             |  5 ++
 src/MNH/init_salt.f90             |  5 +-
 src/MNH/modd_salt.f90             |  7 +--
 src/MNH/mode_salt_psd.f90         | 73 +++++++------------------
 src/MNH/modeln.f90                |  2 +
 src/MNH/salt_filter.f90           | 88 ++++++++++++++-----------------
 src/MNH/saltcamsn.f90             | 40 ++++++++------
 src/MNH/write_lfifm1_for_diag.f90 |  6 +++
 src/SURFEX/coupling_isban.F90     | 35 +++++++++++-
 src/SURFEX/coupling_sltn.F90      | 11 ++--
 src/SURFEX/init_slt.F90           |  4 +-
 src/SURFEX/mode_dslt_surf.F90     |  9 ++--
 13 files changed, 171 insertions(+), 174 deletions(-)

diff --git a/src/MNH/aer2lima.f90 b/src/MNH/aer2lima.f90
index 92798be19..e53c2c91d 100644
--- a/src/MNH/aer2lima.f90
+++ b/src/MNH/aer2lima.f90
@@ -103,10 +103,10 @@ REAL, DIMENSION(SIZE(PSVT,1), SIZE(PSVT,2), SIZE(PSVT,3),6) :: ZAER
 REAL, DIMENSION(SIZE(PSVT,1), SIZE(PSVT,2), SIZE(PSVT,3),NSV) :: ZTOT
 REAL, DIMENSION(SIZE(PSVT,1), SIZE(PSVT,2), SIZE(PSVT,3),JPMODE) :: ZOM
 REAL, DIMENSION(NSV) :: ZMI
-INTEGER :: JSV, JJ, JI, II, IJ, IK, JK
+INTEGER :: JSV, JJ, JI, II, IJ, IK, JK, IDX
 REAL :: ZCCNRADIUS, ZRATMASSH2O
 
-ZCCNRADIUS = 0.03 ! to suppress the aitken mode  (µm)
+ZCCNRADIUS = 0.02 ! to suppress the aitken mode  (µm)
 
 IF ((CPROGRAM=="REAL  ").OR.(CPROGRAM=="IDEAL ")) CMINERAL = "EQSAM"
 
@@ -241,7 +241,6 @@ DO JSV=1,JPMODE
 ! only one CCN_FREE mode (activation is not performed upon aerosol class but by physical paramters)
 !
    WHERE (ZRG_AER(:,:,:,JSV) .GT. ZCCNRADIUS)
-   !WHERE ((ZRG_AER(:,:,:,JSV) .GT. ZCCNRADIUS).AND.(ZRATH2O(:,:,:).GT.ZRATMASSH2O))
      ZCCN_SUM(:,:,:,1) =   ZCCN_SUM(:,:,:,1) + ZN0_AER(:,:,:,JSV) 
    END WHERE
 
@@ -263,22 +262,22 @@ DO JSV=1,JPMODE
 END IF
 
 ! IFN_FREE initialization
-  WHERE (ZRATH2O(:,:,:) .LE. ZRATMASSH2O) ! fraction of dust if low water
-     ZIFN_SUM(:,:,:,1) = ZIFN_SUM(:,:,:,1) + ZN0_AER(:,:,:,JSV) * ZRATDST(:,:,:)
-  END WHERE
+  IF (LDUST) THEN
+    IDX=MIN(NMODE_DST,NMOD_IFN)+JSV
+  ELSE
+    IDX= 1 + JSV
+  END IF 
 
-! hydrophobic aerosols water mass less than 20%
-  IF (NMOD_IFN .GE. 2) THEN
-   WHERE (ZRATH2O(:,:,:) .LE. ZRATMASSH2O) ! hydrophobic aerosols can act as IFN
-     ZIFN_SUM(:,:,:,2) = ZIFN_SUM(:,:,:,2) + ZN0_AER(:,:,:,JSV) * (1.- ZRATSO4(:,:,:))
+  IF (IDX .LE. NMOD_IFN) THEN
+   ZIFN_SUM(:,:,:,IDX) = 0.
+   WHERE (ZRATH2O(:,:,:) .LE. ZRATMASSH2O) ! fraction of dust if low water
+     ZIFN_SUM(:,:,:,IDX) = ZIFN_SUM(:,:,:,IDX) + ZN0_AER(:,:,:,JSV) * ZRATDST(:,:,:)
    END WHERE
   END IF
-
 END DO
 
-
 ELSE ! keep lima class intiatialization
-  IF (CACTCCN=="ABRK") THEN
+   IF (CACTCCN=="ABRK") THEN
 ! only one CCN_FREE mode (activation is not performed upon aerosol class but by physical paramters)
      IF (NMOD_CCN .GE. 2) THEN
         DO JI = 2, NMOD_CCN
@@ -286,22 +285,24 @@ ELSE ! keep lima class intiatialization
                         PSVT(:,:,:,NSV_LIMA_CCN_FREE+JI-1) + PSVT(:,:,:,NSV_LIMA_CCN_ACTI+JI-1)
         END DO
      END IF
-  ELSE
+   ELSE
      IF (NMOD_CCN .GE. 2) THEN
         DO JI = 2, NMOD_CCN
            ZCCN_SUM(:,:,:,JI) = PSVT(:,:,:,NSV_LIMA_CCN_FREE+JI-1) + PSVT(:,:,:,NSV_LIMA_CCN_ACTI+JI-1)
         END DO
      END IF
-  END IF 
+   END IF 
 
-  IF (.NOT.(LDUST) .AND. NMOD_IFN.GE.1) &
-       ZIFN_SUM(:,:,:,1) = PSVT(:,:,:,NSV_LIMA_IFN_FREE)   + PSVT(:,:,:,NSV_LIMA_IFN_NUCL)
+   IF (LDUST) THEN
+    IDX =MIN(NMODE_DST,NMOD_IFN)+1
+   ELSE
+    IDX = 2
+   END IF
+
+   DO JI = IDX, NMOD_IFN
+     ZIFN_SUM(:,:,:,JI) = PSVT(:,:,:,NSV_LIMA_IFN_FREE+JI-1) + PSVT(:,:,:,NSV_LIMA_IFN_NUCL+JI-1)
+   END DO
 
-  IF (NMOD_IFN .GE. 2) THEN
-     DO JI = 2, NMOD_IFN
-        ZIFN_SUM(:,:,:,JI) = PSVT(:,:,:,NSV_LIMA_IFN_FREE+JI-1) + PSVT(:,:,:,NSV_LIMA_IFN_NUCL+JI-1)
-     END DO
-  END IF
 END IF ! end if sur LORILAM
 !
 !
@@ -346,24 +347,16 @@ IF (LDUST) THEN
   CALL PPP2DUST(PSVT(:,:,:,NSV_DSTBEG:NSV_DSTEND),PRHODREF,&
                 PSIG3D=ZSIG_DST,PRG3D=ZRG_DST,PN3D=ZN0_DST)
 !
-  DO JSV=1,NMODE_DST
-
-! #/m3 --> #/kg
-    ZN0_DST(:,:,:,JSV) = ZN0_DST(:,:,:,JSV) / PRHODREF(:,:,:)
-
-! IFN_FREE initialization (all dusts)
-    ZIFN_SUM(:,:,:,1) = ZIFN_SUM(:,:,:,1) + ZN0_DST(:,:,:,JSV)
-
-  END DO
+  IDX = MIN(NMODE_DST-1,NMOD_IFN)
+  ZIFN_SUM(:,:,:,1) = ZN0_DST(:,:,:,2)  / PRHODREF(:,:,:)
+  ZIFN_SUM(:,:,:,2) = ZN0_DST(:,:,:,3)  / PRHODREF(:,:,:)
 
 ELSE ! keep lima class intiatialization
    IF (NMOD_IFN.GE.1) &
         ZIFN_SUM(:,:,:,1) = PSVT(:,:,:,NSV_LIMA_IFN_FREE) + PSVT(:,:,:,NSV_LIMA_IFN_NUCL)
-
 END IF  ! endif sur LDUST
 !
 !
-!
 IF (NMOD_CCN .GE. 1) THEN
    DO JI=1,NMOD_CCN
       PSVT(:,:,:,NSV_LIMA_CCN_FREE+JI-1)   = MAX(ZCCN_SUM(:,:,:,JI) - PSVT(:,:,:,NSV_LIMA_CCN_ACTI+JI-1), 0.)
@@ -376,5 +369,4 @@ IF (NMOD_IFN .GE. 1) THEN
    END DO
 END IF
 !
-!
 END SUBROUTINE AER2LIMA
diff --git a/src/MNH/dustcamsn.f90 b/src/MNH/dustcamsn.f90
index 94732c05c..5a49cba92 100644
--- a/src/MNH/dustcamsn.f90
+++ b/src/MNH/dustcamsn.f90
@@ -129,6 +129,11 @@ DO JN = 1, NMODE_DST
   IF (JPDUSTORDER(JN) == 2) ZMASS(:,:,:,JN) = MAX(PMASSCAMS(:,:,:,2), 1E-15) ! median mode 
   IF (JPDUSTORDER(JN) == 3) ZMASS(:,:,:,JN) = MAX(PMASSCAMS(:,:,:,3), 1E-15) ! large mode
 
+  IF ((JPDUSTORDER(JN) == 1).AND.(XINIRADIUS(JN)==0.039)) &
+                      ZMASS(:,:,:,JN) = 258./(258.+861.) * PMASSCAMS(:,:,:,1) !  fin mode  of AMMA size distribution
+  IF ((JPDUSTORDER(JN) == 2).AND.(XINIRADIUS(JN)== 0.3205)) &
+                      ZMASS(:,:,:,JN) = 861./(258.+861.) * PMASSCAMS(:,:,:,1) + PMASSCAMS(:,:,:,2) ! median mode of AMMA size distribution
+
 ENDDO
 
 !
diff --git a/src/MNH/init_salt.f90 b/src/MNH/init_salt.f90
index bb99ab50f..ab9857ac7 100644
--- a/src/MNH/init_salt.f90
+++ b/src/MNH/init_salt.f90
@@ -43,9 +43,10 @@ IF ( NMODE_SLT == 8) THEN
 !Initial dry number median radius (um) from Ova et al., 2014 + MB21 (Bruch et al., 2022).
 XINIRADIUS_SLT=  (/0.009, 0.021, 0.045, 0.115, 0.415,2.5, 7.0, 25.0/)
 !Initial, standard deviation from  Ova et al., 2014
-XINISIG_SLT =  (/ 1.37, 1.5, 1.42, 1.53, 1.85,1.7, 1.8, 2.1 /)
+XINISIG_SLT =  (/ 1.37, 1.5, 1.42, 1.53, 1.85,1.55, 1.8, 2.1 /)
 !Minimum allowed number concentration for any mode (#/m3)
-XN0MIN_SLT  = (/1.e1 , 1.e1, 1.e1, 1., 1.e-4, 1.e-5, 1.e-6, 1.e-7 /)
+XN0MIN_SLT  = (/100. , 100., 100., 10., 1., 1E-2 , 1E-3, 1E-4 /)
+
 
 ELSE IF ( NMODE_SLT == 3) THEN
 
diff --git a/src/MNH/modd_salt.f90 b/src/MNH/modd_salt.f90
index c8ebbebd0..97cfdd833 100644
--- a/src/MNH/modd_salt.f90
+++ b/src/MNH/modd_salt.f90
@@ -78,12 +78,13 @@ INTEGER :: NMODE_SLT= 8  ! number of sea salt modes (default = 8)
 !REAL, DIMENSION(5) :: XINIRADIUS_SLT,XINISIG_SLT,XN0MIN_SLT
 
 !Initial dry number median radius (um) from Ova et al., 2014
-REAL,DIMENSION(8)    :: XINIRADIUS_SLT=  (/0.009, 0.021, 0.045, 0.115,0.415,2.5, 7.0, 20.0/)
+REAL,DIMENSION(8)    :: XINIRADIUS_SLT=  (/0.009, 0.021, 0.045, 0.115,0.415, 2.5, 7.0, 25.0/)
 !Initial, standard deviation from  Ova et al., 2014
-REAL,DIMENSION(8)      :: XINISIG_SLT =  (/ 1.37, 1.5, 1.42, 1.53, 1.85,1.7,1.8, 2.9 /)
+REAL,DIMENSION(8)      :: XINISIG_SLT =  (/ 1.37, 1.5, 1.42, 1.53, 1.85, 1.55, 1.8, 2.1 /)
 
 !Minimum allowed number concentration for any mode (#/m3)
-REAL,DIMENSION(8)  :: XN0MIN_SLT  = (/1.e1 , 1.e1, 1.e1, 1., 1.e-4,1.e-5, 1.e-6,1.e-7 /)
+REAL,DIMENSION(8)  :: XN0MIN_SLT  = (/100. , 100., 100., 10., 1., 1E-2 , 1E-3, 1E-4 /)
+
 !Test Thomas
 REAL, DIMENSION(:,:,:), ALLOCATABLE :: XSLTMSS   ! [kg/m3] total mass concentration of sea salt
 !
diff --git a/src/MNH/mode_salt_psd.f90 b/src/MNH/mode_salt_psd.f90
index 28da713be..0f765c41d 100644
--- a/src/MNH/mode_salt_psd.f90
+++ b/src/MNH/mode_salt_psd.f90
@@ -158,7 +158,6 @@ DO JN=1,NMODE_SLT
   ZRGMIN         = ZINIRADIUS(JN)
   IF (LVARSIG_SLT) THEN
     ZSIGMIN = XINISIG_SLT(IMODEIDX)
- !   ZSIGMIN = XSIGMIN_SLT
   ELSE
     ZSIGMIN = XINISIG_SLT(IMODEIDX)
   ENDIF
@@ -197,8 +196,6 @@ DO JN=1,NMODE_SLT
          * XAVOGADRO                           & !==> um6/mole_{air}
          / XMD                                 & !==> um6/kg_{air}
          * PRHODREF(:,:,:)                       !==> um6/m3_{air}
-     !Limit m6 concentration to minimum value
-!     ZM(:,:,:,NM6(JN)) =  MAX(ZM(:,:,:,NM6(JN)), ZMMIN(NM6(JN)))
   !
   !Get sigma (only if sigma is allowed to vary)
     !Get intermediate values for sigma M3^2/(M0*M6) (ORILAM paper, eqn 8)
@@ -212,8 +209,6 @@ DO JN=1,NMODE_SLT
     !Finally get the real sigma the negative sign is because of 
     !The way the equation is written (M3^2/(M0*M6)) instead of (M0*M6)/M3^3
     ZSIGMA(:,:,:)= EXP(1./3.*SQRT(-ZSIGMA(:,:,:)))
-    !Limit the value to reasonable ones
-    ZSIGMA(:,:,:) =  MAX( XSIGMIN_SLT, MIN( XSIGMAX_SLT, ZSIGMA(:,:,:) ) )
 
   !
     !Put back M6 so that it fits the sigma which is possibly modified above
@@ -236,10 +231,6 @@ DO JN=1,NMODE_SLT
          * (1.d0/ZRHOI)                        & !==>m3_{aer}/m3_{air}
          * XM3TOUM3_SALT                            & !==>um3_{aer}/m3_{air}
          / (XPI * 4./3.)                         !==>um3_{aer}/m3_{air} (volume ==> 3rd moment)
-!    ZM(:,:,:,NM3(JN)) = MAX(ZM(:,:,:,NM3(JN)), ZMMIN(NM3(JN)))
-!Modif salt/dust 5.1. beg
-       PSVT(:,:,:,JN) = ZM(:,:,:,NM3(JN)) * XMD * XPI * 4./3. * ZRHOI  / &
-                              (ZMI*PRHODREF(:,:,:)*XM3TOUM3_SALT)
 !Modif salt/dust 5.1. end
 
     ZM(:,:,:,NM0(JN))=  ZM(:,:,:,NM3(JN))/&
@@ -275,22 +266,6 @@ DO JN=1,NMODE_SLT
   !
   END IF
   !   
-  !Get number median radius (eqn. 7 in Orilam manuscript)
-  ! ++ JORIS DBG ++
-  IF (NVERB ==15) THEN
-    WRITE(*,*) 'SHAPE(ZM) =', SHAPE(ZM)
-    WRITE(*,*) 'MINVAL(ZM), MAXVAL(ZM) =', MINVAL(ZM), MAXVAL(ZM)
-    WRITE(*,*) 'MINLOC(ZM), MAXLOC(ZM) =', MINLOC(ZM), MAXLOC(ZM)
-    WRITE(*,*) 'SHAPE(ZRG) =', SHAPE(ZRG)
-    WRITE(*,*) 'MINVAL(ZRG), MAXVAL(ZRG) =', MINVAL(ZRG), MAXVAL(ZRG)
-    WRITE(*,*) 'MINLOC(ZRG), MAXLOC(ZRG) =', MINLOC(ZRG), MAXLOC(ZRG)
-    WRITE(*,*) 'XSIXTH_SALT =', XSIXTH_SALT
-    WRITE(*,*) 'JN =', JN
-    WRITE(*,*) 'NM0 =', NM0
-    WRITE(*,*) 'NM3 =', NM3
-    WRITE(*,*) 'NM6 =', NM6
-  ENDIF
-  ! -- JORIS DBG --
   ZRG(:,:,:)=    &
           (      &
           ZM(:,:,:,NM3(JN))*ZM(:,:,:,NM3(JN))*ZM(:,:,:,NM3(JN))*ZM(:,:,:,NM3(JN))    &
@@ -465,7 +440,17 @@ END SUBROUTINE PPP2SALT
           * (1.d0/ZRHOI)                        & !==>m3_{aer}/m3_{air}
           * XM3TOUM3_SALT                            & !==>um3_{aer}/m3_{air}
           / (XPI * 4./3.)                         !==>um3_{aer}/m3_{air} (volume ==> 3rd moment)
-  ELSE 
+     ZM(:,:,:,NM0(JN))=                         &
+         PSVT(:,:,:,1+(JN-1)*3)                & !#/molec_{air}
+         * XAVOGADRO                           & !==>#/mole
+         / XMD                                 & !==>#/kg_{air}
+         * PRHODREF(:,:,:)                       !==>#/m3
+    ZM(:,:,:,NM6(JN)) = PSVT(:,:,:,3+(JN-1)*3)  & !um6/molec_{air}*(cm3/m3)
+         * 1.d-6                               & !==> um6/molec_{air}
+         * XAVOGADRO                           & !==> um6/mole_{air}
+         / XMD                                 & !==> um6/kg_{air}
+         * PRHODREF(:,:,:)                       !==> um6/m3_{air}
+    ELSE 
     IF ((LRGFIX_SLT)) THEN
          ZM(:,:,:,NM3(JN)) =                   &
           PSVT(:,:,:,JN)                        & !molec_{aer}/molec_{aer}
@@ -483,33 +468,14 @@ END SUBROUTINE PPP2SALT
           * (1.d0/ZRHOI)                        & !==>m3_{aer}/m3_{air}
           * XM3TOUM3_SALT                            & !==>um3_{aer}/m3_{air}
           / (XPI * 4./3.)                         !==>um3_{aer}/m3_{air} (volume ==> 3rd moment)
+     ZM(:,:,:,NM0(JN))=                         &
+         PSVT(:,:,:,1+(JN-1)*2)                 & !#/molec_{air}
+         * XAVOGADRO                           & !==>#/mole
+         / XMD                                 & !==>#/kg_{air}
+         * PRHODREF(:,:,:)                       !==>#/m3
+
+     END IF
      END IF
-   END IF
-! calculate moment 0 from dispersion and mean radius
-     ZM(:,:,:,NM0(JN))=  ZM(:,:,:,NM3(JN))/&
-       ((PRG3D(:,:,:,JN)**3)*EXP(4.5 * LOG(PSIG3D(:,:,:,JN))**2))
-
-
-! calculate moment 6 from dispersion and mean radius
-     ZM(:,:,:,NM6(JN)) = ZM(:,:,:,NM0(JN)) * (PRG3D(:,:,:,JN)**6) * &
-               EXP(18 *(LOG(PSIG3D(:,:,:,JN)))**2)
-
-!     IF (LVARSIG_SLT) THEN
-!    WHERE ((ZM(:,:,:,NM0(JN)) .LT. ZMMIN(NM0(JN))).OR.&
-!           (ZM(:,:,:,NM3(JN)) .LT. ZMMIN(NM3(JN))).OR.&
-!           (ZM(:,:,:,NM6(JN)) .LT. ZMMIN(NM6(JN))))
-!    ZM(:,:,:,NM0(JN)) = ZMMIN(NM0(JN))
-!    ZM(:,:,:,NM3(JN)) = ZMMIN(NM3(JN))
-!    ZM(:,:,:,NM6(JN)) = ZMMIN(NM6(JN))
-!    END WHERE
-!    ELSE  IF (.NOT.(LRGFIX_SLT)) THEN
-
-!    WHERE ((ZM(:,:,:,NM0(JN)) .LT. ZMMIN(NM0(JN))).OR.&
-!           (ZM(:,:,:,NM3(JN)) .LT. ZMMIN(NM3(JN))))
-!    ZM(:,:,:,NM0(JN)) = ZMMIN(NM0(JN))
-!    ZM(:,:,:,NM3(JN)) = ZMMIN(NM3(JN))
-!    END WHERE
-!    ENDIF
 
      
      ! return to concentration #/m3 =>  (#/molec_{air}
@@ -638,9 +604,6 @@ ALLOCATE (ZRG(SIZE(PSVT,1)))
 ALLOCATE (ZSV(SIZE(PSVT,1), SIZE(PSVT,2)))
 ALLOCATE (ZINIRADIUS(NMODE_SLT))
 
-!Modif salt/dust 5.1. beg
-ZSV(:,:) = MAX(PSVT(:,:), XMNH_TINY)
-!Modif salt/dust 5.1. end
 
 DO JN=1,NMODE_SLT
   IMODEIDX = JPSALTORDER(JN)
diff --git a/src/MNH/modeln.f90 b/src/MNH/modeln.f90
index e988f0afe..a438833a6 100644
--- a/src/MNH/modeln.f90
+++ b/src/MNH/modeln.f90
@@ -633,6 +633,8 @@ CALL GET_DIM_EXT_ll('B',IIU,IJU)
 IKU=NKMAX+2*JPVEXT
 CALL GET_INDICE_ll (IIB,IJB,IIE,IJE)
 !
+CALL FILL_DIMPHYEX(YLDIMPHYEX, SIZE(XWT,1), SIZE(XWT,2), SIZE(XWT,3))
+!
 IF (IMI==1) THEN
   GSTEADY_DMASS=LSTEADYLS
 ELSE
diff --git a/src/MNH/salt_filter.f90 b/src/MNH/salt_filter.f90
index bc40e8891..8fc62336b 100644
--- a/src/MNH/salt_filter.f90
+++ b/src/MNH/salt_filter.f90
@@ -57,11 +57,6 @@ END MODULE MODI_SALT_FILTER
 USE MODD_SALT
 USE MODE_SALT_PSD
 
-!+ Marine
-USE MODI_INIT_SALT
-!- Marine
-
-!
 IMPLICIT NONE
 !
 !*       0.1   Declarations of dummy arguments :
@@ -83,12 +78,9 @@ INTEGER,DIMENSION(:),    ALLOCATABLE  :: NM6                 ! [idx] indexes for
 REAL,DIMENSION(:),       ALLOCATABLE  :: ZINIRADIUS          ! initial mean radius
 REAL,DIMENSION(:),       ALLOCATABLE  :: ZINISIGMA           ! initial standard deviation
 INTEGER                               :: JN,IMODEIDX,JJ      ! [idx] loop counters
+REAL,DIMENSION(:,:,:,:),       ALLOCATABLE  :: ZSIG, ZN0, ZRG
 
 !-------------------------------------------------------------------------------
-!+ Marine
-CALL INIT_SALT
-!- Marine
-
 
 ALLOCATE (NM0(NMODE_SLT))
 ALLOCATE (NM3(NMODE_SLT))
@@ -116,14 +108,9 @@ DO JN=1,NMODE_SLT
   NM6(JN) = 3+(JN-1)*3
   !Get minimum values possible
   ZMMIN(NM0(JN)) = XN0MIN_SLT(IMODEIDX)
-  ZRGMIN         = ZINIRADIUS(JN)
-  IF (LVARSIG_SLT) THEN
-    ZSIGMIN = XSIGMIN_SLT
-  ELSE
-    ZSIGMIN = XINISIG_SLT(IMODEIDX)
-  ENDIF
-  ZMMIN(NM3(JN)) = ZMMIN(NM0(JN)) * (ZRGMIN**3)*EXP(4.5 * LOG(ZSIGMIN)**2)
-  ZMMIN(NM6(JN)) = ZMMIN(NM0(JN)) * (ZRGMIN**6)*EXP(18. * LOG(ZSIGMIN)**2)
+  ZMMIN(NM3(JN)) = ZMMIN(NM0(JN)) * (ZINIRADIUS(JN)**3)*EXP(4.5 * LOG(ZINISIGMA(JN))**2)
+  ZMMIN(NM6(JN)) = ZMMIN(NM0(JN)) * (ZINIRADIUS(JN)**6)*EXP(18. * LOG(ZINISIGMA(JN))**2)
+
 END DO
 !
 !Set density of aerosol, here dust (kg/m3)
@@ -141,8 +128,27 @@ DO JN=1,NMODE_SLT
           * (1.d0/ZRHOI)                        & !==>m3_{aer}/m3_{air}
           * XM3TOUM3_SALT                            & !==>um3_{aer}/m3_{air}
           / (XPI * 4./3.)                         !==>um3_{aer}/m3_{air} (volume ==> 3rd moment)
-    ELSE
-       IF ((LRGFIX_SLT)) THEN
+       !Get number concentration (#/molec_{air}==>#/m3)
+    ZM(:,:,:,NM0(JN))=                         &
+         PSV(:,:,:,1+(JN-1)*3)                 & !#/molec_{air}
+         * XAVOGADRO                           & !==>#/mole
+         / XMD                                 & !==>#/kg_{air}
+         * PRHODREF(:,:,:)                       !==>#/m3
+
+    ZM(:,:,:,NM6(JN)) = PSV(:,:,:,3+(JN-1)*3)  & !um6/molec_{air}*(cm3/m3)
+         * 1.d-6                               & !==> um6/molec_{air}
+         * XAVOGADRO                           & !==> um6/mole_{air}
+         / XMD                                 & !==> um6/kg_{air}
+         * PRHODREF(:,:,:)                       !==> um6/m3_{air}
+
+     WHERE ((ZM(:,:,:,NM0(JN)) .LT. ZMMIN(NM0(JN))).OR.&
+            (ZM(:,:,:,NM3(JN)) .LT. ZMMIN(NM3(JN))).OR.&
+            (ZM(:,:,:,NM6(JN)) .LT. ZMMIN(NM6(JN))))
+     ZM(:,:,:,NM0(JN)) = ZMMIN(NM0(JN))
+     ZM(:,:,:,NM3(JN)) = ZMMIN(NM3(JN))
+     ZM(:,:,:,NM6(JN)) = ZMMIN(NM6(JN))
+     END WHERE
+    ELSE IF ((LRGFIX_SLT)) THEN
           ZM(:,:,:,NM3(JN)) =                   &
           PSV(:,:,:,JN)                        & !molec_{aer}/molec_{aer}
           * (ZMI/XMD)                           & !==>kg_{aer}/kg_{aer}
@@ -150,7 +156,9 @@ DO JN=1,NMODE_SLT
           * (1.d0/ZRHOI)                        & !==>m3_{aer}/m3_{air}
           * XM3TOUM3_SALT                            & !==>um3_{aer}/m3_{air}
           / (XPI * 4./3.)                         !==>um3_{aer}/m3_{air} (volume ==> 3rd moment)
-      ELSE
+
+         ZM(:,:,:,NM3(JN)) = MAX(ZM(:,:,:,NM3(JN)),ZMMIN(NM3(JN)))
+    ELSE
           ZM(:,:,:,NM3(JN)) =                   &
           PSV(:,:,:,2+(JN-1)*2)                & !molec_{aer}/molec_{aer}
           * (ZMI/XMD)                           & !==>kg_{aer}/kg_{aer}
@@ -158,36 +166,20 @@ DO JN=1,NMODE_SLT
           * (1.d0/ZRHOI)                        & !==>m3_{aer}/m3_{air}
           * XM3TOUM3_SALT                            & !==>um3_{aer}/m3_{air}
           / (XPI * 4./3.)                         !==>um3_{aer}/m3_{air} (volume ==> 3rd moment)
-      END IF
+  
+         ZM(:,:,:,NM0(JN))=                         &
+         PSV(:,:,:,1+(JN-1)*2)                 & !#/molec_{air}
+         * XAVOGADRO                           & !==>#/mole
+         / XMD                                 & !==>#/kg_{air}
+         * PRHODREF(:,:,:)                       !==>#/m3
+
+       WHERE ((ZM(:,:,:,NM0(JN)) .LT. ZMMIN(NM0(JN))).OR.&
+              (ZM(:,:,:,NM3(JN)) .LT. ZMMIN(NM3(JN))))
+         ZM(:,:,:,NM0(JN)) = ZMMIN(NM0(JN))
+         ZM(:,:,:,NM3(JN)) = ZMMIN(NM3(JN))
+       END WHERE
     END IF
 
-! calculate moment 0 from dispersion and mean radius
-     ZM(:,:,:,NM0(JN))=  ZM(:,:,:,NM3(JN))/&
-       ((ZINIRADIUS(JN)**3)*EXP(4.5 * LOG(ZINISIGMA(JN))**2))
-
-
-! calculate moment 6 from dispersion and mean radius
-     ZM(:,:,:,NM6(JN)) = ZM(:,:,:,NM0(JN)) * (ZINIRADIUS(JN)**6) * &
-               EXP(18 *(LOG(ZINISIGMA(JN)))**2)
-
-     IF (LVARSIG_SLT) THEN
-     WHERE ((ZM(:,:,:,NM0(JN)) .LT. ZMMIN(NM0(JN))).OR.&
-            (ZM(:,:,:,NM3(JN)) .LT. ZMMIN(NM3(JN))).OR.&
-            (ZM(:,:,:,NM6(JN)) .LT. ZMMIN(NM6(JN))))
-     ZM(:,:,:,NM0(JN)) = ZMMIN(NM0(JN))
-     ZM(:,:,:,NM3(JN)) = ZMMIN(NM3(JN))
-     ZM(:,:,:,NM6(JN)) = ZMMIN(NM6(JN))
-     END WHERE
-
-     ELSE IF (.NOT.(LRGFIX_SLT)) THEN
-
-     WHERE ((ZM(:,:,:,NM0(JN)) .LT. ZMMIN(NM0(JN))).OR.&
-            (ZM(:,:,:,NM3(JN)) .LT. ZMMIN(NM3(JN))))
-     ZM(:,:,:,NM0(JN)) = ZMMIN(NM0(JN))
-     ZM(:,:,:,NM3(JN)) = ZMMIN(NM3(JN))
-     END WHERE
-     ENDIF
-
     ! return to concentration #/m3 =>  (#/molec_{air}
      IF (LVARSIG_SLT) THEN
      PSV(:,:,:,1+(JN-1)*3) = ZM(:,:,:,NM0(JN)) * XMD / &
diff --git a/src/MNH/saltcamsn.f90 b/src/MNH/saltcamsn.f90
index 96dcde95d..82f062938 100644
--- a/src/MNH/saltcamsn.f90
+++ b/src/MNH/saltcamsn.f90
@@ -79,14 +79,14 @@ REAL,DIMENSION(:),       ALLOCATABLE  :: ZINIRADIUS, ZINISIGMA
 INTEGER :: IKU, IMOMENTS
 INTEGER :: JJ, JN, JK  ! loop counter
 INTEGER :: IMODEIDX  ! index mode
-REAL    :: ZRHOMIN
+REAL    :: ZRHOMAX
 
 REAL :: DELTA_1,DELTA_2,DELTA_3,DELTA_4,DELTA_5,DELTA_6,DELTA_7
 REAL :: RATIO_1,RATIO_2,RATIO_3,RATIO_4,RATIO_5, RATIO_6,RATIO_7
 REAL :: DELTA_CAMS_1,DELTA_CAMS_2,DELTA_CAMS_3
 REAL :: RAY_CAMS_1,RAY_CAMS_2,RAY_CAMS_3,RAY_CAMS_4
 REAL :: RAY_2,RAY_3,RAY_4
-REAL,DIMENSION(:,:,:,:), ALLOCATABLE  ::  ZMASS_TEST
+REAL,DIMENSION(:,:,:,:), ALLOCATABLE  ::  ZMASS_TEST, ZN0, ZRG, ZSIG
 !
 !-------------------------------------------------------------------------------
 !
@@ -97,7 +97,7 @@ REAL,DIMENSION(:,:,:,:), ALLOCATABLE  ::  ZMASS_TEST
 !
 CALL INIT_SALT
 IKU = SIZE(PSV,3)
-ZRHOMIN=MINVAL(PRHODREF)
+ZRHOMAX=MAXVAL(PRHODREF)
 !
 ALLOCATE (IM0(NMODE_SLT))
 ALLOCATE (IM3(NMODE_SLT))
@@ -109,6 +109,9 @@ ALLOCATE (ZINIRADIUS(NMODE_SLT))
 ALLOCATE (ZINISIGMA(NMODE_SLT))
 ALLOCATE (ZMMIN(NMODE_SLT*3))
 ALLOCATE (ZMASS(SIZE(PSV,1), SIZE(PSV,2), SIZE(PSV,3),NMODE_SLT))
+ALLOCATE (ZRG(SIZE(PSV,1), SIZE(PSV,2), SIZE(PSV,3),NMODE_SLT))
+ALLOCATE (ZSIG(SIZE(PSV,1), SIZE(PSV,2), SIZE(PSV,3),NMODE_SLT))
+ALLOCATE (ZN0(SIZE(PSV,1), SIZE(PSV,2), SIZE(PSV,3),NMODE_SLT))
 !
 ! Rayons des bins CAMS
 
@@ -163,14 +166,15 @@ ZMASS(:,:,:,5) = (PMASSCAMS(:,:,:,3) * RATIO_6) + ZMASS(:,:,:,5) ! Attribution M
 
 !========================================================
 ! Adjust the mass / SSA emissions after a few hours
-ZMASS(:,:,:,1) = MAX(ZMASS(:,:,:,1) * 0.16, 1E-18)
+ZMASS(:,:,:,1) = MAX(ZMASS(:,:,:,1) * 0.3, 1E-18)
 ZMASS(:,:,:,2) = MAX(ZMASS(:,:,:,2) * 0.1, 1E-17)
 ZMASS(:,:,:,3) = MAX(ZMASS(:,:,:,3) * 0.5, 1E-16)
-ZMASS(:,:,:,4) = MAX(ZMASS(:,:,:,4) * 0.1, 1E-16)
-ZMASS(:,:,:,5) = MAX(ZMASS(:,:,:,5), 1E-17) 
-IF (NMODE_SLT >= 6) ZMASS(:,:,:,6) = MAX(ZMASS(:,:,:,5) * 0.01, 1E-16)
-IF (NMODE_SLT >= 7) ZMASS(:,:,:,7) = MAX(ZMASS(:,:,:,5) * 0.001, 1E-16)
-IF (NMODE_SLT >= 8) ZMASS(:,:,:,8) = MAX(ZMASS(:,:,:,5) * 0.0001, 1E-16)
+ZMASS(:,:,:,4) = MAX(ZMASS(:,:,:,4) * 0.3, 1E-15)
+ZMASS(:,:,:,5) = MAX(ZMASS(:,:,:,5) * 0.1, 1E-14) 
+IF (NMODE_SLT >= 6) ZMASS(:,:,:,6) = MAX(ZMASS(:,:,:,5) * 0.1, 1E-14)
+IF (NMODE_SLT >= 7) ZMASS(:,:,:,7) = MAX(ZMASS(:,:,:,5) * 0.1, 1E-14)
+IF (NMODE_SLT >= 8) ZMASS(:,:,:,8) = MAX(ZMASS(:,:,:,5) * 0.1, 1E-14)
+
 !========================================================
 
 DO JN = 1, NMODE_SLT
@@ -213,7 +217,8 @@ DO JN = 1, NMODE_SLT
                     / (2.d0 * ZINIRADIUS(JN) * 1.d-6)**3 &![particle/m_salt^{-3}]==> particle/m3
                     * EXP(-4.5*(LOG(ZINISIGMA(JN)))**2) !Take into account distribution
 
-  ZM(:,:,:,IM0(JN)) = MAX(ZMMIN(IM0(JN)), ZM(:,:,:,IM0(JN)))
+  ZM(:,:,:,IM0(JN)) = MAX(10.*ZMMIN(IM0(JN)), ZM(:,:,:,IM0(JN)))
+
 !
 !*       1.2    calculate moment 3 from m0,  RG and SIG 
 !
@@ -234,29 +239,30 @@ DO JN = 1, NMODE_SLT
   IMOMENTS = INT(NSV_SLTEND - NSV_SLTBEG+1) / NMODE_SLT
   IF (IMOMENTS == 3) THEN
     PSV(:,:,:,1+(JN-1)*3) = ZM(:,:,:,IM0(JN)) * XMD / (XAVOGADRO*PRHODREF(:,:,:))
-    XSVMIN(NSV_SLTBEG-1+1+(JN-1)*3) = ZMMIN(IM0(JN)) * XMD / (XAVOGADRO*ZRHOMIN) 
+    XSVMIN(NSV_SLTBEG-1+1+(JN-1)*3) = ZMMIN(IM0(JN)) * XMD / (XAVOGADRO*ZRHOMAX) 
 
     PSV(:,:,:,2+(JN-1)*3) = ZM(:,:,:,IM3(JN)) * XMD * XPI * 4. / 3. * ZRHOI / &
                             (ZMI*XM3TOUM3_SALT*PRHODREF(:,:,:))
     XSVMIN(NSV_SLTBEG-1+2+(JN-1)*3) = ZMMIN(IM3(JN))  * XMD * XPI * 4. / 3. * ZRHOI / &
-                            (ZMI*XM3TOUM3_SALT**ZRHOMIN)
+                            (ZMI*XM3TOUM3_SALT*ZRHOMAX)
+  
 
     PSV(:,:,:,3+(JN-1)*3) = ZM(:,:,:,IM6(JN)) *  XMD / (XAVOGADRO*1.d-6*PRHODREF(:,:,:))
-    XSVMIN(NSV_SLTBEG-1+3+(JN-1)*3) = ZMMIN(IM6(JN))  *  XMD / (XAVOGADRO*1.d-6* ZRHOMIN)
+    XSVMIN(NSV_SLTBEG-1+3+(JN-1)*3) = ZMMIN(IM6(JN))  *  XMD / (XAVOGADRO*1.d-6* ZRHOMAX)
+
   ELSE IF (IMOMENTS == 2) THEN
     PSV(:,:,:,1+(JN-1)*2) = ZM(:,:,:,IM0(JN)) * XMD / (XAVOGADRO*PRHODREF(:,:,:))
-    XSVMIN(NSV_SLTBEG-1+1+(JN-1)*2) = ZMMIN(IM0(JN)) * XMD / (XAVOGADRO*ZRHOMIN)
+    XSVMIN(NSV_SLTBEG-1+1+(JN-1)*2) = ZMMIN(IM0(JN)) * XMD / (XAVOGADRO*ZRHOMAX)
 
     PSV(:,:,:,2+(JN-1)*2) = ZM(:,:,:,IM3(JN)) * XMD * XPI * 4./3. * ZRHOI / &
                             (ZMI*XM3TOUM3_SALT*PRHODREF(:,:,:))
     XSVMIN(NSV_SLTBEG-1+2+(JN-1)*2) = ZMMIN(IM3(JN))  * XMD * XPI * 4. / 3. * ZRHOI / &
-                            (ZMI*XM3TOUM3_SALT**ZRHOMIN)
-
+                            (ZMI*XM3TOUM3_SALT*ZRHOMAX)
   ELSE 
     PSV(:,:,:,JN) = ZM(:,:,:,IM3(JN)) * XMD * XPI * 4. / 3. * ZRHOI / &
                     (ZMI * XM3TOUM3_SALT*PRHODREF(:,:,:))
     XSVMIN(NSV_SLTBEG-1+JN) = ZMMIN(IM3(JN))  * XMD * XPI * 4. / 3. * ZRHOI / &
-                            (ZMI*XM3TOUM3_SALT**ZRHOMIN)
+                            (ZMI*XM3TOUM3_SALT*ZRHOMAX)
 
   END IF
 END DO
diff --git a/src/MNH/write_lfifm1_for_diag.f90 b/src/MNH/write_lfifm1_for_diag.f90
index abea790b2..dcb367276 100644
--- a/src/MNH/write_lfifm1_for_diag.f90
+++ b/src/MNH/write_lfifm1_for_diag.f90
@@ -1496,6 +1496,12 @@ IF ((LCHEMDIAG).AND.(LORILAM).AND.(LUSECHEM)) THEN
     TZFIELD%CUNITS     = 'ug m-3'
     WRITE(TZFIELD%CCOMMENT,'(A21,I1)')'MASS BC AEROSOL MODE ',JJ
     CALL IO_Field_write(TPFILE,TZFIELD,ZPTOTA(:,:,:,JP_AER_BC,JJ))
+    !
+    WRITE(TZFIELD%CMNHNAME,'(A4,I1)')'MDST',JJ
+    TZFIELD%CLONGNAME  = TRIM(TZFIELD%CMNHNAME)
+    TZFIELD%CUNITS     = 'ug m-3'
+    WRITE(TZFIELD%CCOMMENT,'(A22,I1)')'MASS DST AEROSOL MODE ',JJ
+    CALL IO_Field_write(TPFILE,TZFIELD,ZPTOTA(:,:,:,JP_AER_DST,JJ))
   ENDDO
 END IF
 ! Dust variables
diff --git a/src/SURFEX/coupling_isban.F90 b/src/SURFEX/coupling_isban.F90
index f70ce3371..087bcc7b7 100644
--- a/src/SURFEX/coupling_isban.F90
+++ b/src/SURFEX/coupling_isban.F90
@@ -106,6 +106,8 @@ USE MODD_ISBA_n,          ONLY : ISBA_S_t, ISBA_K_t, ISBA_P_t, ISBA_PE_t, &
                                  ISBA_NK_t, ISBA_NP_t, ISBA_NPE_t
 !
 USE MODD_DST_n,           ONLY : DST_NP_t, DST_t
+USE MODD_CHS_AEROSOL,     ONLY : XEMISSIGI, XEMISSIGJ,&
+                                 XEMISRADIUSI, XEMISRADIUSJ, CRGUNIT
 USE MODD_SLT_n,           ONLY : SLT_t
 !
 USE MODD_REPROD_OPER,     ONLY : CIMPLICIT_WIND
@@ -348,6 +350,7 @@ REAL, DIMENSION(KI) :: ZRNSUNLIT
 REAL                       :: ZCONVERTFACM0_SLT, ZCONVERTFACM0_DST
 REAL                       :: ZCONVERTFACM3_SLT, ZCONVERTFACM3_DST
 REAL                       :: ZCONVERTFACM6_SLT, ZCONVERTFACM6_DST
+REAL                       :: ZEMISRADIUSI, ZEMISRADIUSJ
 !
 ! for blowing snow scheme
 !
@@ -775,6 +778,7 @@ REAL, DIMENSION(PK%NSIZE_P) :: ZP_ZS      ! atmospheric model orography
 REAL, DIMENSION(PK%NSIZE_P) :: ZP_SFTQ    ! flux of water vapor <w'q'>            (kg.m-2.s-1)
 REAL, DIMENSION(PK%NSIZE_P) :: ZP_SFTH    ! flux of temperature <w'T'>            (W/m2)
 REAL, DIMENSION(PK%NSIZE_P,KSV) :: ZP_SFTS    ! flux of scalar      <w'sv'>           (mkg/kg/s)
+REAL, DIMENSION(PK%NSIZE_P) :: ZF_DSTI, ZF_DSTJ 
 REAL, DIMENSION(PK%NSIZE_P) :: ZP_SFCO2   ! flux of CO2 positive toward the atmosphere (m/s*kg_CO2/kg_air)
 REAL, DIMENSION(PK%NSIZE_P) :: ZP_USTAR   ! friction velocity                     (m/s)
 REAL, DIMENSION(PK%NSIZE_P) :: ZP_SFU     ! zonal momentum flux                   (pa)
@@ -785,6 +789,7 @@ REAL, DIMENSION(PK%NSIZE_P) :: ZP_Z0      ! roughness length for momentum (m)
 REAL, DIMENSION(PK%NSIZE_P) :: ZP_Z0H     ! roughness length for heat     (m)
 REAL, DIMENSION(PK%NSIZE_P):: ZP_QSURF   ! specific humidity at surface  (kg/kg)
 REAL, DIMENSION(PK%NSIZE_P) :: ZP_TEMP, ZP_PAR
+
 !
 !*  other forcing variables (packed for each patch)
 !
@@ -1407,15 +1412,41 @@ IF(CHI%SVI%NDSTEQ>0)THEN
          JSV_IDX = (JMODE-1)*2
        END IF
        !
+       IF (CRGUNIT=="MASS") THEN
+        ZEMISRADIUSI = XEMISRADIUSI * EXP(-3.*(LOG(XEMISSIGI))**2)
+        ZEMISRADIUSJ = XEMISRADIUSJ * EXP(-3.*(LOG(XEMISSIGJ))**2)
+       ELSE
+        ZEMISRADIUSI = XEMISRADIUSI
+        ZEMISRADIUSJ = XEMISRADIUSJ
+       END IF
+
+
+       ZF_DSTI(:) = 0.
+       ZF_DSTJ(:) = 0.
        DO JSV=1, size(HSV)  ! dust flux for orilam chemistry model
          IF ((TRIM(HSV(JSV)) == "@DSTI").AND.(JMODE==3)) THEN 
+           ! dust mass flux in kg/m2/s
+           ZF_DSTI(:) = ZP_SFTS(:,IBEG-1+JSV_IDX+2)
            ! add dust flux and conversion kg/m2/s into molec.m2/s
-           ZP_SFTS(:,JSV) = ZP_SFTS(:,JSV) + ZP_SFTS(:,IBEG-1+JSV_IDX+2)*XAVOGADRO/XMOLARWEIGHT_DST
+           ZP_SFTS(:,JSV) = ZP_SFTS(:,JSV) + ZF_DSTI(:)*XAVOGADRO/XMOLARWEIGHT_DST
+           ! dust flux in moment 0
+           ZF_DSTI(:) = ZF_DSTI(:) * 1E9 ! µg/m2/s
+           ZF_DSTI(:) = ZF_DSTI(:) / ((4./3.)*3.14292654*XDENSITY_DST*1.e-9) !moment 3
+           ZF_DSTI(:) = ZF_DSTI(:) / ((ZEMISRADIUSI**3)*EXP(4.5 * (LOG(XEMISSIGI))**2)) ! moment 0 molecules.m-2.s-1
          END IF
+         IF (TRIM(HSV(JSV)) == "@M0I") ZP_SFTS(:,JSV) = ZP_SFTS(:,JSV) + ZF_DSTI(:)
+
          IF ( (TRIM(HSV(JSV)) == "@DSTJ").AND.(JMODE==2)) THEN 
+           ! dust mass flux in kg/m2/s
+           ZF_DSTJ(:) = ZP_SFTS(:,IBEG-1+JSV_IDX+2)
            ! add dust flux and conversion kg/m2/sec into molec.m2/s
-           ZP_SFTS(:,JSV) = ZP_SFTS(:,JSV) + ZP_SFTS(:,IBEG-1+JSV_IDX+2)*XAVOGADRO/XMOLARWEIGHT_DST
+           ZP_SFTS(:,JSV) = ZP_SFTS(:,JSV) + ZF_DSTJ(:)*XAVOGADRO/XMOLARWEIGHT_DST
+           ! dust flux in moment 0
+           ZF_DSTJ(:) = ZF_DSTJ(:) * 1E9 ! µg/m2/s
+           ZF_DSTJ(:) = ZF_DSTJ(:) / ((4./3.)*3.14292654*XDENSITY_DST*1.e-9) !moment 3
+           ZF_DSTJ(:) = ZF_DSTJ(:) / ((ZEMISRADIUSJ**3)*EXP(4.5 * (LOG(XEMISSIGJ))**2)) ! moment 0 molecules.m-2.s-1
          END IF
+         IF (TRIM(HSV(JSV)) == "@M0J") ZP_SFTS(:,JSV) = ZP_SFTS(:,JSV) + ZF_DSTJ(:)
        END DO
        !
      END DO
diff --git a/src/SURFEX/coupling_sltn.F90 b/src/SURFEX/coupling_sltn.F90
index cdb6c10ee..47e9bd14d 100644
--- a/src/SURFEX/coupling_sltn.F90
+++ b/src/SURFEX/coupling_sltn.F90
@@ -230,16 +230,15 @@ IF ((CEMISPARAM_SLT .eq. "Ova14").OR.(CEMISPARAM_SLT .eq. "OvB21a").OR.(CEMISPAR
     ZSFSLT_MDE(:,2) = 0.044  * ( ZREYNOLDS(:) - 1.E5)**1.08
     ZSFSLT_MDE(:,3) = 149.64 * ( ZREYNOLDS(:) - 1.E5)**0.545
     ZSFSLT_MDE(:,4) = 2.96   * ( ZREYNOLDS(:) - 1.E5)**0.79
-  ENDWHERE
-  WHERE (ZREYNOLDS(:) > 2.E5)
-    ZSFSLT_MDE(:,5) = 0.52   * ( ZREYNOLDS(:) - 2.E5)**0.87
-  ENDWHERE
-
-   WHERE (ZREYNOLDS(:) <= 1.E5)
+  ELSEWHERE
     ZSFSLT_MDE(:,1) = 1.E-10
     ZSFSLT_MDE(:,2) = 1.E-10
     ZSFSLT_MDE(:,3) = 1.E-10
     ZSFSLT_MDE(:,4) = 1.E-10
+  ENDWHERE
+  WHERE (ZREYNOLDS(:) > 2.E5)
+    ZSFSLT_MDE(:,5) = 0.52   * ( ZREYNOLDS(:) - 2.E5)**0.87
+  ELSEWHERE
     ZSFSLT_MDE(:,5) = 1.E-10
   ENDWHERE
 
diff --git a/src/SURFEX/init_slt.F90 b/src/SURFEX/init_slt.F90
index c74c91f12..4de2abcc7 100644
--- a/src/SURFEX/init_slt.F90
+++ b/src/SURFEX/init_slt.F90
@@ -67,13 +67,13 @@ IF (LHOOK) CALL DR_HOOK('INIT_SLT',0,ZHOOK_HANDLE)
   NSLTMDE = 8
   CRGUNITS   = 'NUMB'
   XEMISRADIUS_INI_SLT = (/0.009, 0.021, 0.045, 0.115, 0.415, 2.5, 7.0, 25.0/)
-  XEMISSIG_INI_SLT = (/1.37, 1.5, 1.42, 1.53,1.70,1.80, 1.85, 2.1/)
+  XEMISSIG_INI_SLT = (/1.37, 1.5, 1.42, 1.53,1.85, 1.55, 1.8, 2.1/)
 
 IF ((CEMISPARAM_SLT.eq."OvB21a").OR.(CEMISPARAM_SLT.eq."OvB21b")) THEN
   NSLTMDE = 8
   CRGUNITS   = 'NUMB'
   XEMISRADIUS_INI_SLT = (/0.009, 0.021, 0.045, 0.115, 0.415, 2.5, 7.0, 25.0/)
-  XEMISSIG_INI_SLT = (/1.37, 1.5, 1.42, 1.53,1.70,1.80, 1.85, 2.1/)
+  XEMISSIG_INI_SLT = (/1.37, 1.5, 1.42, 1.53,1.85,1.55, 1.8, 2.1/)
 
 ELSE IF (CEMISPARAM_SLT.eq."Ova14") THEN
   NSLTMDE = 5
diff --git a/src/SURFEX/mode_dslt_surf.F90 b/src/SURFEX/mode_dslt_surf.F90
index 742a4de2f..8e121876d 100644
--- a/src/SURFEX/mode_dslt_surf.F90
+++ b/src/SURFEX/mode_dslt_surf.F90
@@ -126,11 +126,10 @@ DO JMODE=1,KMDE
   !IF TWO MOMENTS PER MODE, MASS FLUX IS SENT AS INDEX #2, #4, #6
 
   !Get flux of number in #/m2/sec from flux of mass in kg/m2/sec
-  ZFM(:,1) = PFLUX(:,JSV_IDX+2)                     & ! kg_{dst}/m2/sec
-             / PEMISRADIUS(JMODE)**3                & ! *um^{-3} ==> #/m2/sec*(m3/um3)
-             * EXP(-4.5*(LOG(PEMISSIG(JMODE)))**2)  & ! Take into account size distribution  
-             / PCONVERTFACM3           ! /(kg_{dst}/m^{3}_{dst)} ==> m^3_{dst}/m2/sec  
- 
+   ZFM(:,1) = PFLUX(:,JSV_IDX+2)                     & ! kg_{dst}/m2/sec
+             / (PEMISRADIUS(JMODE)**3                & ! *um^{-3} ==> #/m2/sec*(m3/um3)
+             * EXP(4.5*(LOG(PEMISSIG(JMODE)))**2)  & ! Take into account size distribution
+             * PCONVERTFACM3)           ! /(kg_{dst}/m^{3}_{dst)} ==> m^3_{dst}/m2/sec
    
   ! Get flux of moment 6 consistent with the other moments
   ZFM(:,3) = ZFM(:,1)                              & ! [#/m3]
-- 
GitLab