! $Id: stratemit.F90  2022-07-04 mmarchand $

SUBROUTINE STRATDISTRIB(altLMDz,altemiss,sigma_alt,f_lay_emiss)

  USE dimphy, ONLY: klon,klev
  USE strataer_local_var_mod
  USE lmdz_yomcst, ONLY: RPI

  IMPLICIT NONE

! Input argument
!---------------
    REAL,DIMENSION(klev+1)  :: altLMDz    ! altitude of layer interfaces in m
    REAL                    :: sigma_alt
    REAL                    :: altemiss

! Local variables
!----------------
    REAL                 :: f_lay_sum     ! sum of layer emission fractions
    REAL                 :: alt
    INTEGER              :: k,i_int
    INTEGER,PARAMETER    :: n_int_alt=10  ! subintervals nb for integration over Gaussian emission

! Output argument
!----------------
    REAL,DIMENSION(klev) :: f_lay_emiss   ! emiss fraction for every emiss for every vertical layer
    
    
! Inject distrib
!----------------
    SELECT CASE(flag_emit_distrib)
       
    CASE(0) ! Gaussian distribution
            ! Compute distribution of emission to vertical model layers
            ! (based on Gaussian peak in altitude)
       f_lay_sum=0.0
       DO k=1, klev
          f_lay_emiss(k)=0.0
          DO i_int=1, n_int_alt
             alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
             f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt)* &
   exp(-0.5*((alt-altemiss)/sigma_alt)**2.)* &
   (altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
          ENDDO
          f_lay_sum=f_lay_sum+f_lay_emiss(k)
        ENDDO
        
     CASE(1) ! Uniform distribution
             ! In this case, parameter sigma_alt is considered
             ! to be half the
             ! height of the injection, centered around altemiss)
       f_lay_sum=0.0
       DO k=1, klev
          f_lay_emiss(k)=max(min(altemiss+sigma_alt,altLMDz(k+1))- &
   max(altemiss-sigma_alt, &
   altLMDz(k)),0.)/(2.*sigma_alt)
          f_lay_sum=f_lay_sum+f_lay_emiss(k)
       ENDDO

       END SELECT        ! End CASE over flag_emit_distrib)
       
       WRITE(*,*) "IN stratdistrib f_lay_emiss=",f_lay_emiss

       !correct for step integration error
       f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum

END SUBROUTINE STRATDISTRIB