Index: /trunk/LMDZ.PLUTO/libf/phypluto/mp2m_calmufi.F90
===================================================================
--- /trunk/LMDZ.PLUTO/libf/phypluto/mp2m_calmufi.F90	(revision 3559)
+++ /trunk/LMDZ.PLUTO/libf/phypluto/mp2m_calmufi.F90	(revision 3559)
@@ -0,0 +1,187 @@
+MODULE mp2m_calmufi
+    use tracer_h
+    use callkeys_mod, only : call_haze_prod_pCH4, haze_rho
+    
+    ! Microphysical model MP2M
+    use mp2m_intgcm
+    use mp2m_diagnostics
+    
+
+    implicit none
+    
+    !============================================================================
+    !
+    !     Purpose
+    !     -------
+    !     Interface subroutine to YAMMS model for LMD PCM.
+    !
+    !     The subroutine computes the microphysics processes for a single vertical column.
+    !       - All input vectors are assumed to be defined from GROUND to TOP of the atmosphere.
+    !       - All output vectors are defined from GROUND to TOP of the atmosphere.
+    !       - Only tendencies are returned.
+    !
+    !     @important
+    !     The method assumes global initialization of YAMMS model (and extras) has been already
+    !     done elsewhere.
+    !
+    !     @warning
+    !     Microphysical tracers from physics must be in X/kg_of_air and convert into X/m2 for microphysics.
+    !
+    !     @warning
+    !     We suppose a given order of tracers (1. mu_m0as, 2. mu_m3as, 3. mu_m0af, 4. mu_m3af)!
+    !
+    !     Authors
+    !     -------
+    !     B. de Batz de Trenquelléon (11/2024)
+    !
+    !============================================================================
+
+    CONTAINS
+    
+    SUBROUTINE calmufi(dt, plev, zlev, play, zlay, g3d, temp, pq, zdqfi, zdqmufi_prod, zdqmufi)
+        !! Interface subroutine to YAMMS model for LMD PCM.
+        !!
+        !! The subroutine computes the microphysics processes for a single vertical column.
+        !!
+        REAL(kind=8), INTENT(IN) :: dt  !! Physics timestep (s).
+        
+        REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: plev            ! Pressure levels (Pa).
+        REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: zlev            ! Altitude levels (m).
+        REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: play            ! Pressure layers (Pa).
+        REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: zlay            ! Altitude at the center of each layer (m).
+        REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: g3d             ! Latitude-Altitude depending gravitational acceleration (m.s-2).
+        REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: temp            ! Temperature at the center of each layer (K).
+        
+        REAL(kind=8), DIMENSION(:,:,:), INTENT(IN)  :: pq           ! Tracers (X.kg-1).
+        REAL(kind=8), DIMENSION(:,:,:), INTENT(IN)  :: zdqfi        ! Tendency from former processes for tracers (X.kg-1.s-1).
+        REAL(kind=8), DIMENSION(:,:,:), INTENT(IN)  :: zdqmufi_prod ! Aerosols production tendency (kg/kg_of_air/s).
+        REAL(kind=8), DIMENSION(:,:,:), INTENT(OUT) :: zdqmufi      ! Microphysical tendency for tracers (X.m-2 --> X.kg-1.s-1).
+        
+        ! Local tracers:
+        !~~~~~~~~~~~~~~~
+        REAL(kind=8), DIMENSION(:,:,:), ALLOCATABLE :: zq    ! Local tracers updated from former processes (X.kg-1).
+        
+        REAL(kind=8), DIMENSION(:), ALLOCATABLE :: m0as      ! 0th order moment of the spherical mode (m-2).
+        REAL(kind=8), DIMENSION(:), ALLOCATABLE :: m3as      ! 3rd order moment of the spherical mode (m3.m-2).
+        REAL(kind=8), DIMENSION(:), ALLOCATABLE :: m0af      ! 0th order moment of the fractal mode (m-2).
+        REAL(kind=8), DIMENSION(:), ALLOCATABLE :: m3af      ! 3rd order moment of the fractal mode (m3.m-2).
+
+        REAL(kind=8), DIMENSION(:), ALLOCATABLE :: m3as_prod ! Production of 3rd order moment of the spherical mode (m3.m-2).
+        
+        REAL(kind=8), DIMENSION(:), ALLOCATABLE   :: dm0as   ! Tendency of the 0th order moment of the spherical mode distribution (m-2).
+        REAL(kind=8), DIMENSION(:), ALLOCATABLE   :: dm3as   ! Tendency of the 3rd order moment of the spherical mode distribution (m3.m-2).
+        REAL(kind=8), DIMENSION(:), ALLOCATABLE   :: dm0af   ! Tendency of the 0th order moment of the fractal mode distribution (m-2).
+        REAL(kind=8), DIMENSION(:), ALLOCATABLE   :: dm3af   ! Tendency of the 3rd order moment of the fractal mode distribution (m3.m-2).
+        
+        ! Local variables:
+        !~~~~~~~~~~~~~~~~~
+        REAL(kind=8), DIMENSION(:,:), ALLOCATABLE :: int2ext ! Conversion intensive to extensive (kg.m-2).
+        REAL(kind=8), DIMENSION(:), ALLOCATABLE   :: tmp 
+        TYPE(error) :: err
+        
+        INTEGER           :: ilon,iq
+        INTEGER           :: nq,nlon,nlay
+        CHARACTER(len=10) :: tname
+        
+        ! Read size of arrays:
+        !~~~~~~~~~~~~~~~~~~~~~
+        nq    = size(pq,DIM=3)
+        nlon  = size(play,DIM=1)
+        nlay  = size(play,DIM=2)
+        
+        ! Allocate arrays:
+        !~~~~~~~~~~~~~~~~~
+        ALLOCATE(zq(nlon,nlay,nq))
+        
+        ALLOCATE(m0as(nlay))
+        ALLOCATE(m3as(nlay))
+        ALLOCATE(m0af(nlay))
+        ALLOCATE(m3af(nlay))
+
+        ALLOCATE(m3as_prod(nlay))
+        
+        ALLOCATE(dm0as(nlay))
+        ALLOCATE(dm3as(nlay))
+        ALLOCATE(dm0af(nlay))
+        ALLOCATE(dm3af(nlay))
+
+        ALLOCATE(int2ext(nlon,nlay))  
+        
+        !------------------
+        ! 1. Initialization
+        !------------------
+
+        ! Initialization of zdqmufi here since intent=out and no action performed on every tracers
+        zdqmufi(:,:,:) = 0.D0
+        
+        ! Initialize tracers updated with former processes from physics
+        zq(:,:,:) = pq(:,:,:) + zdqfi(:,:,:)*dt
+        
+        ! Loop on horizontal grid points
+        DO ilon = 1, nlon
+
+            ! Convert tracers to extensive
+            int2ext(ilon,:) = (plev(ilon,1:nlay)-plev(ilon,2:nlay+1)) / g3d(ilon,1:nlay)
+            
+            m0as(:) = zq(ilon,:,micro_indx(1)) * int2ext(ilon,:)
+            m3as(:) = zq(ilon,:,micro_indx(2)) * int2ext(ilon,:)
+            
+            m0af(:) = zq(ilon,:,micro_indx(3)) * int2ext(ilon,:)
+            m3af(:) = zq(ilon,:,micro_indx(4)) * int2ext(ilon,:)
+
+            ! Production of haze in the atmosphere by photolysis of CH4
+            if (call_haze_prod_pCH4) then
+                do iq = 1, nq
+                    tname = noms(iq)
+                    if (tname(1:4).eq."haze") then
+                        m3as_prod(:) = zdqmufi_prod(ilon,:,iq) * (int2ext(ilon,:) / haze_rho) * dt
+                    endif
+                enddo
+            else
+                m3as_prod(:) = 0.D0
+            endif
+            
+            ! Hackin the pressure level
+            tmp = plev(ilon,:)
+            if (tmp(nlay+1) == 0.0) then
+            tmp(nlay+1) = 2*tmp(nlay) - tmp(nlay-1)
+            endif
+        
+            ! Initialize YAMMS atmospheric column
+            err = mm_column_init(tmp,zlev(ilon,:),play(ilon,:),zlay(ilon,:),temp(ilon,:)) ; IF (err /= 0) call abort_program(err)
+            ! Initialize YAMMS aerosols moments column
+            err = mm_aerosols_init(m0as,m3as,m0af,m3af) ; IF (err /= 0) call abort_program(err)
+            
+            ! Initializes tendencies
+            dm0as(:) = 0._mm_wp ; dm3as(:) = 0._mm_wp ; dm0af(:) = 0._mm_wp ; dm3af(:) = 0._mm_wp
+            
+        !----------------------------
+        ! 2. Call microphysical model
+        !----------------------------
+            
+            ! Call microphysics
+            IF (.NOT.mm_muphys(m3as_prod,dm0as,dm3as,dm0af,dm3af)) THEN
+                call abort_program(error("mm_muphys aborted -> initialization not done !",-1))
+            ENDIF
+            
+            ! Save diagnostics
+            call mm_diagnostics(dt,mp2m_aer_s_prec(ilon),mp2m_aer_f_prec(ilon),  &
+                                mp2m_aer_s_w(ilon,:),mp2m_aer_f_w(ilon,:),       &
+                                mp2m_aer_s_flux(ilon,:),mp2m_aer_f_flux(ilon,:), &
+                                mp2m_rc_sph(ilon,:),mp2m_rc_fra(ilon,:))
+        
+            ! Convert tracers back to intensives
+            zdqmufi(ilon,:,micro_indx(1)) = dm0as(:) / int2ext(ilon,:)
+            zdqmufi(ilon,:,micro_indx(2)) = dm3as(:) / int2ext(ilon,:)
+            zdqmufi(ilon,:,micro_indx(3)) = dm0af(:) / int2ext(ilon,:)
+            zdqmufi(ilon,:,micro_indx(4)) = dm3af(:) / int2ext(ilon,:)
+            
+        END DO ! End loop on ilon
+
+        ! YAMMS gives a tendency which is integrated for all the timestep but in the GCM 
+        ! we want to have routines spitting tendencies in s-1 -> let's divide !
+        zdqmufi(:,:,:) = zdqmufi(:,:,:) / dt
+        
+    END SUBROUTINE calmufi
+
+END MODULE mp2m_calmufi
Index: /trunk/LMDZ.PLUTO/libf/phypluto/mp2m_diagnostics.F90
===================================================================
--- /trunk/LMDZ.PLUTO/libf/phypluto/mp2m_diagnostics.F90	(revision 3559)
+++ /trunk/LMDZ.PLUTO/libf/phypluto/mp2m_diagnostics.F90	(revision 3559)
@@ -0,0 +1,133 @@
+MODULE mp2m_diagnostics
+
+!=======================================================================
+!   Purpose
+!   -------
+!   Microphysics diagnostics variables module.
+!
+!   The module contains all the variables related to output diagnostics of the microphysics two moments.
+!   Such variables are not (and should not be) used as input in the GCM except for output data writing.
+!
+!   The module also contains three methods:
+!      - inimuphy_diag(ngrid,nlayer,nq,pq,int2ext)
+!      - alpha_s(k)
+!      - alpha_f(k)
+!
+!   Authors
+!   -------
+!      Bruno de Batz de Trenquelléon (11/2024)
+!=======================================================================
+
+  USE tracer_h
+  
+  IMPLICIT NONE
+
+  REAL(kind=8), ALLOCATABLE, DIMENSION(:)     :: mp2m_aer_s_prec ! Spherical aerosols precipitations (kg.m-2.s-1).
+  REAL(kind=8), ALLOCATABLE, DIMENSION(:)     :: mp2m_aer_f_prec ! Fractal aerosols precipitations (kg.m-2.s-1).
+  REAL(kind=8), ALLOCATABLE, DIMENSION(:,:)   :: mp2m_aer_s_w    ! Spherical aerosol settling velocity (m.s-1).
+  REAL(kind=8), ALLOCATABLE, DIMENSION(:,:)   :: mp2m_aer_f_w    ! Fractal aerosol settling velocity (m.s-1).
+  REAL(kind=8), ALLOCATABLE, DIMENSION(:,:)   :: mp2m_aer_s_flux ! Spherical aerosol mass flux (kg.m-2.s-1).
+  REAL(kind=8), ALLOCATABLE, DIMENSION(:,:)   :: mp2m_aer_f_flux ! Fractal aerosol mass flux (kg.m-2.s-1).
+  REAL(kind=8), ALLOCATABLE, DIMENSION(:,:)   :: mp2m_rc_sph     ! Spherical mode characteristic radius, i.e. bulk radius (m).
+  REAL(kind=8), ALLOCATABLE, DIMENSION(:,:)   :: mp2m_rc_fra     ! Fractal mode characteristic radius, i.e. bulk radius (m).
+  
+  !$OMP THREADPRIVATE(mp2m_aer_s_prec,mp2m_aer_f_prec,mp2m_aer_s_w,mp2m_aer_f_w,mp2m_aer_s_flux,mp2m_aer_f_flux)
+  !$OMP THREADPRIVATE(mp2m_rc_sph,mp2m_rc_fra)
+
+  CONTAINS
+
+  SUBROUTINE inimufi_diag(ngrid,nlayer,nq,pq,int2ext)
+    !! Initialize the variables associated to microphysics diagnostics.
+    !!
+    use callkeys_mod, only : haze_rc_prod, haze_rm
+
+    INTEGER, INTENT(in) :: ngrid  ! Number of points of the horizontal grid.
+    INTEGER, INTENT(in) :: nlayer ! Number of points of the vertical grid (layers).
+    INTEGER, INTENT(in) :: nq     ! Number of advected fields.
+
+    REAL,INTENT(in) :: pq(ngrid,nlayer,nq)   ! Tracers (X/kg_of_air).
+    REAL,INTENT(in) :: int2ext(ngrid,nlayer) ! Intensive to extensive factor (kg_air/m3).
+
+    ! Local variables:
+    !~~~~~~~~~~~~~~~~~
+    REAL, DIMENSION(:,:), ALLOCATABLE :: m0as ! 0th order moment of the spherical mode (m-3).
+    REAL, DIMENSION(:,:), ALLOCATABLE :: m3as ! 3rd order moment of the spherical mode (m3.m-3).
+    REAL, DIMENSION(:,:), ALLOCATABLE :: m0af ! 0th order moment of the fractal mode (m-3).
+    REAL, DIMENSION(:,:), ALLOCATABLE :: m3af ! 3rd order moment of the fractal mode (m3.m-3).
+    
+    IF (.NOT.ALLOCATED(mp2m_aer_s_prec)) ALLOCATE(mp2m_aer_s_prec(ngrid))
+    IF (.NOT.ALLOCATED(mp2m_aer_f_prec)) ALLOCATE(mp2m_aer_f_prec(ngrid))
+    IF (.NOT.ALLOCATED(mp2m_aer_s_w))    ALLOCATE(mp2m_aer_s_w(ngrid,nlayer))
+    IF (.NOT.ALLOCATED(mp2m_aer_f_w))    ALLOCATE(mp2m_aer_f_w(ngrid,nlayer))
+    IF (.NOT.ALLOCATED(mp2m_aer_s_flux)) ALLOCATE(mp2m_aer_s_flux(ngrid,nlayer))
+    IF (.NOT.ALLOCATED(mp2m_aer_f_flux)) ALLOCATE(mp2m_aer_f_flux(ngrid,nlayer))
+    IF (.NOT.ALLOCATED(mp2m_rc_sph))     ALLOCATE(mp2m_rc_sph(ngrid,nlayer))
+    IF (.NOT.ALLOCATED(mp2m_rc_fra))     ALLOCATE(mp2m_rc_fra(ngrid,nlayer))
+    
+    mp2m_aer_s_prec(:)   = 0d0
+    mp2m_aer_f_prec(:)   = 0d0
+    mp2m_aer_s_w(:,:)    = 0d0
+    mp2m_aer_f_w(:,:)    = 0d0
+    mp2m_aer_s_flux(:,:) = 0d0
+    mp2m_aer_f_flux(:,:) = 0d0
+    
+    ! Initializes the radius to their first value (useful for radiative transfer).
+    m0as = pq(:,:,micro_indx(1)) * int2ext
+    m3as = pq(:,:,micro_indx(2)) * int2ext
+    m0af = pq(:,:,micro_indx(3)) * int2ext
+    m3af = pq(:,:,micro_indx(4)) * int2ext
+
+    WHERE(m0as > 1e-10 .AND. m3as > (1e-10*alpha_s(3.)*haze_rc_prod**3))
+      mp2m_rc_sph = (m3as / (m0as*alpha_s(3.)))**(1./3.)
+    ELSEWHERE
+      mp2m_rc_sph = 0d0
+    ENDWHERE
+
+    WHERE(m0af > 1e-10 .AND. m3af > (1e-10*alpha_f(3.)*haze_rm**3))
+      mp2m_rc_fra = (m3af / (m0af*alpha_f(3.)))**(1./3.)
+    ELSEWHERE
+      mp2m_rc_fra = 0d0
+    ENDWHERE
+  
+    CONTAINS
+
+    FUNCTION alpha_s(k) RESULT (res)
+      !! Inter-moment relation for spherical aerosols size distribution law.
+      !! Mk / M0 = rc^k . alpha(k)
+      !!
+      use mp2m_intgcm
+      REAL, INTENT(in) :: k ! k order of the moment.
+      REAL :: sigma         ! Standard deviation.
+      REAL :: res           ! Alpha value.
+    
+      ! Titan's case
+      !~~~~~~~~~~~~~
+      ! res = SUM(dexp(mm_asp%a*k**2 + mm_asp%b*k + mm_asp%c))
+
+      ! Pluto's case
+      !~~~~~~~~~~~~~
+      sigma = 0.2
+      res = exp(k**2 * sigma**2 / 2.)
+    END FUNCTION alpha_s
+
+    FUNCTION alpha_f(k) RESULT (res)
+      !! Inter-moment relation for fractal aerosols size distribution law.
+      !! Mk / M0 = rc^k . alpha(k)
+      !!
+      use mp2m_intgcm
+      REAL, INTENT(in) :: k ! k order of the moment.
+      REAL :: sigma         ! Standard deviation.
+      REAL :: res           ! Alpha value.
+    
+      ! Titan's case
+      !~~~~~~~~~~~~~
+      ! res = SUM(dexp(mm_asp%a*k**2 + mm_asp%b*k + mm_asp%c))
+
+      ! Pluto's case
+      !~~~~~~~~~~~~~
+      sigma = 0.35
+      res = exp(k**2 * sigma**2 / 2.)
+    END FUNCTION alpha_f
+  END SUBROUTINE inimufi_diag
+
+END MODULE mp2m_diagnostics
Index: /trunk/LMDZ.PLUTO/libf/phypluto/mp2m_inimufi.F90
===================================================================
--- /trunk/LMDZ.PLUTO/libf/phypluto/mp2m_inimufi.F90	(revision 3559)
+++ /trunk/LMDZ.PLUTO/libf/phypluto/mp2m_inimufi.F90	(revision 3559)
@@ -0,0 +1,67 @@
+subroutine inimufi(ptimestep)
+  use callkeys_mod, only : call_haze_prod_pCH4, haze_p_prod, haze_tx_prod, haze_rc_prod, haze_rm, haze_df, haze_rho, air_rad
+  use tracer_h
+  use comcstfi_mod, only : g, rad, mugaz
+  use datafile_mod
+
+  ! Microphysical model MP2M
+  use mp2m_intgcm
+
+  implicit none
+
+  !============================================================================
+  !
+  !     Purpose
+  !     -------
+  !     This routine call mm_initialize which perform the global initialization
+  !     for the microphysical YAMMS model.
+  !     It also performs some sanity checks on microphysical tracer names. 
+  !
+  !     Authors
+  !     -------
+  !     B. de Batz de Trenquelléon (11/2024)
+  !
+  !============================================================================
+
+
+  !----------------
+  ! 0. Declarations
+  !----------------
+
+  integer :: i,idx
+  real, intent(in) :: ptimestep ! Timestep (s)
+
+  character(len=30), dimension(4), parameter :: aernames = &
+     (/"mu_m0as              ", "mu_m3as              ",   &
+       "mu_m0af              ", "mu_m3af              "/)
+  
+  logical :: err
+
+  !--------------------------
+  ! 1. Check names of tracers
+  !--------------------------
+  err = .false.
+
+  DO i=1,size(aernames)
+    idx = indexOfTracer(TRIM(aernames(i)),.false.)
+    IF (idx <= 0) THEN
+      WRITE(*,*) "inimufi: '"//TRIM(aernames(i))//"' not found in tracers table."
+      err = .true.
+    ENDIF
+  ENDDO
+
+  IF (err) THEN
+    WRITE(*,*) "You loose in inimufi.F90.... Try again"
+    STOP
+  ENDIF
+
+  !--------------------------------------------------
+  ! 2. Call microphysical model global initialization
+  !--------------------------------------------------
+  
+  ! /!\ mugaz [g/mol] --> [kg/mol] !!
+  call mm_initialize(ptimestep,call_haze_prod_pCH4,haze_p_prod,haze_tx_prod,haze_rc_prod, &
+                      haze_rm,haze_df,haze_rho,rad,g,air_rad,mugaz*1e-3,                  &
+                      config_mufi)
+
+end subroutine inimufi
