SUBROUTINE calmufi(plev, zlev, play, zlay, temp, pq, zdq)
  !! Interface subroutine to YAMMS model for Titan LMDZ GCM.
  !!
  !! 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.
  !! 
  !! Authors : J.Burgalat, J.Vatant d'Ollone - 2017
  !!
  USE MMP_GCM
  USE tracer_h
  USE comcstfi_mod, only : g
  USE callkeys_mod, only : callclouds
  IMPLICIT NONE

  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) :: temp  !! Temperature at the center of each layer (K).

  REAL(kind=8), DIMENSION(:,:,:), INTENT(IN)  :: pq   !! Tracers (\(kg.kg^{-1}}\)).
  REAL(kind=8), DIMENSION(:,:,:), INTENT(OUT) :: zdq  !! Tendency for tracers (\(kg.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 (\(m^{3}.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 (\(m^{3}.m^{-2}\)).

  REAL(kind=8), DIMENSION(:), ALLOCATABLE   :: m0n  !! 0th order moment of the CCN distribution (\(m^{-2}\)).
  REAL(kind=8), DIMENSION(:), ALLOCATABLE   :: m3n  !! 3rd order moment of the CCN distribution (\(m^{3}.m^{-2}\)).
  REAL(kind=8), DIMENSION(:,:), ALLOCATABLE :: m3i  !! 3rd order moments of the ice components (\(m^{3}.m^{-2}\)).
  REAL(kind=8), DIMENSION(:,:), ALLOCATABLE :: gazs !! Condensible species gazs molar fraction (\(mol.mol^{-1}\)).

  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 (\(m^{3}.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 (\(m^{3}.m^{-2}\)).
  REAL(kind=8), DIMENSION(:), ALLOCATABLE   :: dm0n  !! Tendency of the 0th order moment of the _CCN_ distribution (\(m^{-2}\)).
  REAL(kind=8), DIMENSION(:), ALLOCATABLE   :: dm3n  !! Tendency of the 3rd order moment of the _CCN_ distribution (\(m^{3}.m^{-2}\)).
  REAL(kind=8), DIMENSION(:,:), ALLOCATABLE :: dm3i  !! Tendencies of the 3rd order moments of each ice components (\(m^{3}.m^{-2}\)). 
  REAL(kind=8), DIMENSION(:,:), ALLOCATABLE :: dgazs !! Tendencies of each condensible gaz species !(\(mol.mol^{-1}\)). 

  REAL(kind=8), DIMENSION(:), ALLOCATABLE ::  int2ext
  TYPE(error) :: err

  INTEGER :: ilon, i,nices

  INTEGER :: nlon,nlay

  ! Read size of arrays
  nlon  = size(play,DIM=1)
  nlay  = size(play,DIM=2)
  nices = size(ices_indx)
  ! Conversion intensive to extensive
  ALLOCATE( int2ext(nlay) )  

  ! Loop on horizontal grid points

  ! Allocate arrays
  ALLOCATE( m0as(nlay) )
  ALLOCATE( m3as(nlay) )
  ALLOCATE( m0af(nlay) )
  ALLOCATE( m3af(nlay) )
  ALLOCATE( m0n(nlay) )
  ALLOCATE( m3n(nlay) )
  ALLOCATE( m3i(nlay,nices) )
  ALLOCATE( gazs(nlay,nices) )

  ALLOCATE( dm0as(nlay) )
  ALLOCATE( dm3as(nlay) )
  ALLOCATE( dm0af(nlay) )
  ALLOCATE( dm3af(nlay) )
  ALLOCATE( dm0n(nlay) )
  ALLOCATE( dm3n(nlay) )
  ALLOCATE( dm3i(nlay,nices) )
  ALLOCATE( dgazs(nlay,nices) )

  ! Initialization of zdq here since intent=out and no action performed on every tracers
  zdq(:,:,:) = 0.0

  DO ilon = 1, nlon
    ! Convert tracers to extensive ( except for gazs where we work with molar mass ratio )
    ! We suppose a given order of tracers !
    int2ext(:) = ( plev(ilon,1:nlay) - plev(ilon,2:nlay+1) ) / g

    m0as(:) = pq(ilon,:,1) * int2ext(:)
    m3as(:) = pq(ilon,:,2) * int2ext(:)
    m0af(:) = pq(ilon,:,3) * int2ext(:)
    m3af(:) = pq(ilon,:,4) * int2ext(:)
    
    if (callclouds) then ! if call clouds
      dm0n(:) = pq(ilon,:,5) * int2ext(:)
      dm3n(:) = pq(ilon,:,6) * int2ext(:)
      do i=1,nices
        dm3i(:,nices) = pq(ilon,:,6+i) * int2ext(:)
        dgazs(:,i)    = pq(ilon,:,ices_indx(i)) * rat_mmol(ices_indx(i)) ! For gazs we work on the full tracer array !!
        ! We use the molar mass ratio from GCM in case there is discrepancy with the mm one
      enddo
    endif


    ! Initialize YAMMS atmospheric column
    err = mm_column_init(plev(ilon,:),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)
    IF (callclouds) THEN ! call clouds
      err = mm_clouds_init(m0n,m3n,m3i,gazs) ; IF (err /= 0) call abort_program(err)
    ENDIF

    ! Check on size (???)

    ! initializes tendencies:
    !dm0as(:) = 0._mm_wp ; dm3as(:) = 0._mm_wp ; dm0af(:) = 0._mm_wp ; dm3af(:) = 0._mm_wp
    !dm0n(:) = 0._mm_wp ; dm3n(:) = 0._mm_wp ; dm3i(:,:) = 0._mm_wp ; dgazs(:,:) = 0._mm_wp

    dm0as(:) = 0.0 ; dm3as(:) = 0.0 ; dm0af(:) = 0.0 ; dm3af(:) = 0.0
    dm0n(:) = 0.0 ; dm3n(:) = 0.0 ; dm3i(:,:) = 0.0 ; dgazs(:,:) = 0.0
    ! call microphysics

    IF (callclouds) THEN ! call clouds
      IF(.NOT.mm_muphys(dm0as,dm3as,dm0af,dm3af,dm0n,dm3n,dm3i,dgazs)) &
        call abort_program(error("mm_muphys aborted -> initialization not done !",-1))
    ELSE
      IF (.NOT.mm_muphys(dm0as,dm3as,dm0af,dm3af)) & 
        call abort_program(error("mm_muphys aborted -> initialization not done !",-1))
    ENDIF

    ! Convert tracers back to intensives ( except for gazs where we work with molar mass ratio )
    ! We suppose a given order of tracers !

    zdq(ilon,:,1) = dm0as(:) / int2ext(:)
    zdq(ilon,:,2) = dm3as(:) / int2ext(:)
    zdq(ilon,:,3) = dm0af(:) / int2ext(:)
    zdq(ilon,:,4) = dm3af(:) / int2ext(:)
    
    if (callclouds) then ! if call clouds
      zdq(ilon,:,5) = dm0n(:) / int2ext(:)
      zdq(ilon,:,6) = dm3n(:) / int2ext(:)
      do i=1,nices
        zdq(ilon,:,6+i) = dm3i(:,nices) / int2ext(:)
        zdq(ilon,:,ices_indx(i)) = dgazs(:,i)  / rat_mmol(ices_indx(i)) ! For gazs we work on the full tracer array !!
        ! We use the molar mass ratio from GCM in case there is discrepancy with the mm one
      enddo
    endif

  END DO ! loop on ilon

END SUBROUTINE calmufi