MODULE MP2M_MICROPHYSICS
    !============================================================================
    !
    !     Purpose
    !     -------
    !     Interface to main microphysics subroutine.
    !     The interface computes all aerosols microphysics processes in a single call.
    !
    !     The module contains two methods:
    !        - mm_muphys(m3as_prod,dm0a_s,dm3a_s,dm0a_f,dm3a_f)
    !        - mm_diagnostics(dt,aer_prec,aer_s_w,aer_f_w,aer_s_flux,aer_f_flux,rc_sph,rc_fra)
    !
    !     Authors
    !     -------
    !     B. de Batz de Trenquelléon, J. Burgalat (11/2024)
    !
    !============================================================================

    USE MP2M_MPREC
    USE MP2M_GLOBALS
    USE MP2M_HAZE
    USE MP2M_METHODS
    IMPLICIT NONE
    
    PUBLIC :: mm_muphys, mm_diagnostics
  
    CONTAINS
  
    FUNCTION mm_muphys(m3as_prod,dm0a_s,dm3a_s,dm0a_f,dm3a_f) RESULT(ret)
        !! Compute the evolution of moments tracers through haze microphysics processes.
        !!
        !! This method computes the evolution of all the microphysics tracers, given under the form
        !! of moments during a time step.
        !!
        !! The method requires that global variables of the model (i.e. variables declared in mm_globals
        !! module) are initialized/updated correctly (see mm_global_init, mm_column_init, and mm_aerosols_init).
        !!
        !! The tendencies returned by the method are defined on the vertical __layers__ of the model from the __GROUND__ to
        !! the __TOP__ of the atmosphere. They should be added to the input variables used in the initialization methods
        !! before the latter are called to initialize a new step.
        !!

        ! Production of the 3rd order moment of the spherical mode distribution (m3.m-2).
        REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: m3as_prod
        ! Tendency of the 0th order moment of the spherical mode distribution (m-2).
        REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm0a_s
        ! Tendency of the 3rd order moment of the spherical mode distribution (m3.m-2).
        REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm3a_s
        ! Tendency of the 0th order moment of the fractal mode distribution (m-2).
        REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm0a_f
        ! Tendency of the 3rd order moment of the fractal mode distribution (m3.m-2).
        REAL(kind=mm_wp), INTENT(out), DIMENSION(:) :: dm3a_f

        ! .true. on succes (i.e. model has been initialized at least once previously), .false. otherwise.
        LOGICAL :: ret
        
        ! Local variables.
        ! Production of the spherical aerosols (m3.m-3).
        REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE :: m3a_s_prod
        ALLOCATE(m3a_s_prod(mm_nla))
        
        ret = (mm_ini_col.AND.mm_ini_aer)
        
        IF (.NOT.ret) RETURN
      
        ! Reverse vectors so they go from top to ground
        ! @note: mm_dzlev is already from top to ground
        m3a_s_prod = m3as_prod(mm_nla:1:-1) / mm_dzlev(:)

        ! Calls haze microphysics
        call mm_haze_microphysics(m3a_s_prod,dm0a_s,dm3a_s,dm0a_f,dm3a_f)
        
        ! Reverse vectors so they go from ground to top
        dm0a_s = dm0a_s(mm_nla:1:-1) * mm_dzlev(mm_nla:1:-1)
        dm3a_s = dm3a_s(mm_nla:1:-1) * mm_dzlev(mm_nla:1:-1)
        dm0a_f = dm0a_f(mm_nla:1:-1) * mm_dzlev(mm_nla:1:-1)
        dm3a_f = dm3a_f(mm_nla:1:-1) * mm_dzlev(mm_nla:1:-1)
        
        RETURN
    END FUNCTION mm_muphys


    SUBROUTINE mm_diagnostics(dt,aer_s_prec,aer_f_prec,aer_s_w,aer_f_w,aer_s_flux,aer_f_flux,rc_sph,rc_fra)
        !! Get various diagnostic fields of the microphysics.
        !!
        !! @note
        !! Fluxes values are always negative as they account for sedimentation fluxes. They are set as
        !! vector and are ordered from __GROUND__ to __TOP__.
        !!
        !! @note
        !! Precipitations are always positive and defined in kg.m-2.s-1.
        !!

        ! Physics timestep (s).
        REAL(kind=8), INTENT(IN) :: dt

        ! Aerosol precipitation (kg.m-2.s-1).
        REAL(kind=mm_wp), INTENT(out), OPTIONAL :: aer_s_prec
        REAL(kind=mm_wp), INTENT(out), OPTIONAL :: aer_f_prec

        ! Aerosol settling velocity (m.s-1).
        REAL(kind=mm_wp), INTENT(out), DIMENSION(:), OPTIONAL :: aer_s_w
        REAL(kind=mm_wp), INTENT(out), DIMENSION(:), OPTIONAL :: aer_f_w

        ! Aerosol mass flux (kg.m-2.s-1).
        REAL(kind=mm_wp), INTENT(out), DIMENSION(:), OPTIONAL :: aer_s_flux
        REAL(kind=mm_wp), INTENT(out), DIMENSION(:), OPTIONAL :: aer_f_flux

        ! Aerosol characteristic radius (m).
        REAL(kind=mm_wp), INTENT(out), DIMENSION(:), OPTIONAL :: rc_sph
        REAL(kind=mm_wp), INTENT(out), DIMENSION(:), OPTIONAL :: rc_fra

        IF (PRESENT(aer_s_prec)) aer_s_prec = ABS(mm_aers_prec) / dt
        IF (PRESENT(aer_f_prec)) aer_f_prec = ABS(mm_aerf_prec) / dt
        IF (PRESENT(aer_s_w))    aer_s_w    = -mm_m3as_vsed(mm_nla:1:-1)
        IF (PRESENT(aer_f_w))    aer_f_w    = -mm_m3af_vsed(mm_nla:1:-1)
        IF (PRESENT(aer_s_flux)) aer_s_flux = -mm_aer_s_flux(mm_nla:1:-1)
        IF (PRESENT(aer_f_flux)) aer_f_flux = -mm_aer_f_flux(mm_nla:1:-1)

        IF (mm_ini_aer) THEN
            IF (PRESENT(rc_sph)) rc_sph = mm_rcs(mm_nla:1:-1)
            IF (PRESENT(rc_fra)) rc_fra = mm_rcf(mm_nla:1:-1)
        ELSE
            IF (PRESENT(rc_sph)) rc_sph = 0._mm_wp
            IF (PRESENT(rc_fra)) rc_fra = 0._mm_wp
        ENDIF    
    END SUBROUTINE mm_diagnostics

END MODULE MP2M_MICROPHYSICS