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

SUBROUTINE STRATEMIT(pdtphys,pdt,xlat,xlon,t_seri,pplay,paprs,tr_seri,&
                    m_emiss_vol_daily,latmin,latmax,lonmin,lonmax,altemiss,sigma_alt,id_spec,&
                    stretchlong,ispeci,id_species_total) 

  USE dimphy, ONLY: klon,klev
  USE strataer_local_var_mod
  USE phys_cal_mod
  USE phys_local_var_mod, ONLY: d_q_emiss
  USE lmdz_yomcst, ONLY: RD, RPI, RG
  USE lmdz_geometry, ONLY: cell_area, boundslat
  USE aerophys
  USE infotrac_phy
  USE lmdz_phys_mpi_data, ONLY:  is_mpi_root
  USE lmdz_grid_phy, ONLY: nbp_lon
 
  IMPLICIT NONE

! Input argument
!---------------
    REAL,INTENT(IN)    :: pdtphys, pdt               ! Pas d'integration pour la physique(seconde)
    REAL,DIMENSION(klon),INTENT(IN)    :: xlat       ! latitudes pour chaque point
    REAL,DIMENSION(klon),INTENT(IN)    :: xlon       ! longitudes pour chaque point
    REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay
    ! pression pour le mileu de chaque couche (en Pa)
    REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs
    ! pression pour chaque inter-couche (en Pa)
    REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri ! Temperature
    REAL,INTENT(IN)    :: latmin,latmax,lonmin,lonmax
    REAL,INTENT(IN)    :: m_emiss_vol_daily          ! daily injection mass emission
    REAL,INTENT(IN)    :: stretchlong                ! length to stretch emission (day)
    REAL,DIMENSION(klon,klev) :: m_air_gridbox       ! mass ofair in every grid box[kg
    REAL,INTENT(IN)     :: sigma_alt, altemiss
    INTEGER,INTENT(IN)  :: id_spec, ispeci, id_species_total
    
! Output argument
!----------------
    REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT)  :: tr_seri ! Concentration Traceur [U/KgA]
! Local variables
!----------------
    INTEGER              :: ieru, i, k, i_int, ilon, ilev
    REAL                 :: dlat_loc
    REAL,DIMENSION(klev+1)  :: altLMDz    ! altitude of layer interfaces in m
    REAL                 :: zrho          ! Density of air [kg/m3]
    REAL,DIMENSION(klev) :: zdm           ! mass of atm. model layer in kg  
    REAL                 :: zdz           ! thickness of atm. model layer in m
    REAL                 :: f_lay_sum     ! sum of layer emission fractions
    REAL,DIMENSION(klev) :: f_lay_emiss   ! emiss fraction for every emiss for every vertical layer
    REAL                 :: emission      ! emission
    REAL                 :: theta_min, theta_max, theta ! for SAI computation between two latitudes
    
    IF (is_mpi_root) THEN
       WRITE(*,*) 'IN STRATEMIT: date from phys_cal_mod=',year_cur,'-',&
    mth_cur,'-',day_cur,'-',hour,' flh2o=',flh2o
    ENDIF

!--calculate mass of air in every grid box !NL: DONE ABOVE DIRECTLY IN MAIN LOOP I,K
!    DO ilon=1, klon
!       DO ilev=1, klev
!          m_air_gridbox(ilon,ilev)=(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG*cell_area(ilon)
!       ENDDO
!    ENDDO
    
    DO i=1,klon
       
       dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3))
       dlon=360./2./FLOAT(nbp_lon) 
       ! dlat = half difference of boundary latitudes
       IF(flag_verbose_strataer) WRITE(*,*) 'IN STRATEMIT, dlat=',dlat_loc,xlat(i),xlon(i)
       
       theta_min = max(xlat(i)-dlat_loc,latmin)
       theta_max = min(xlat(i)+dlat_loc,latmax)
                
       IF (  xlat(i)>=latmin-dlat_loc .AND. &
     xlat(i)<latmax+dlat_loc .AND. &
     xlon(i)>=lonmin-dlon .AND. &
     xlon(i)<lonmax+dlon ) THEN

          WRITE(*,*) 'coordinates of volcanic injection point=',&
      xlat(i),xlon(i),day_cur,mth_cur,year_cur
          WRITE(*,*) 'DD m_emiss_vol_daily=', &
      m_emiss_vol_daily
          
          !compute altLMDz
          altLMDz(:)=0.0
          DO k=1, klev
             zrho=pplay(i,k)/t_seri(i,k)/RD       !air density in kg/m3
             zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG  !mass of layer in kg
             zdz=zdm(k)/zrho                   !thickness of layer in m
             altLMDz(k+1)=altLMDz(k)+zdz       !altitude of interface
          ENDDO
          
          CALL STRATDISTRIB(altLMDz,altemiss,sigma_alt,f_lay_emiss)
          
          IF (flag_emit==3) THEN
             theta=(sin(theta_max/180.*RPI)-sin(theta_min/180.*RPI))/ &
   (sin(xlat_max_sai/180.*RPI)-sin(xlat_min_sai/180.*RPI))
          ELSE
             theta=1.
          ENDIF

          !vertically distributed emission
          DO k=1, klev
            !--calculate mass of air in every grid box
            m_air_gridbox(i,k)=(paprs(i,k)-paprs(i,k+1))/RG*cell_area(i)
            ! stretch emission over stretchlong period
            emission=m_emiss_vol_daily/m_air_gridbox(i,k)*f_lay_emiss(k)/stretchlong/ &
       (86400.-pdt)*theta
            
             IF(flag_verbose_strataer) WRITE(*,*) 'IN STRATEMIT: emission avant/apres', &
    'alt= ', altLMDz(k), &
    'flh2o= ',flh2o, &
    'id_speac= ',id_spec,id_species_total, &
    'emission= ',emission, &
    'pdtphys= ',pdtphys, &
    'rapport m_emiss/m_air*f_lay= ', m_emiss_vol_daily/m_air_gridbox(i,k)*f_lay_emiss(k), &
     'stretchlong= ', stretchlong, &
     'theta= ', theta
             
             IF(emission < 1.E-34) emission = 0.0
             
             IF (flh2o==0) THEN
                IF(flag_verbose_strataer) WRITE(*,*) 'IN STRATEMIT: tr_ser avant/apres',&
    'i= ',i,'k= ',k, 'flh2o= ',flh2o, &
    tr_seri(i,k,id_spec), &
    tr_seri(i,k,id_spec)+emission*pdtphys
              
                tr_seri(i,k,id_spec)=tr_seri(i,k,id_spec)+emission*pdtphys
                IF (id_species_total/=0) THEN
                   tr_seri(i,k,id_species_total)=tr_seri(i,k,id_species_total)+emission*pdtphys
                ENDIF
             ELSE IF(flh2o==1) THEN
                d_q_emiss(i,k)=emission*pdtphys
                IF(d_q_emiss(i,k) > 1.E34) THEN
                   WRITE(*,*) 'IN STRATEMIT: d_q_emiss(i,k) > 1.E34 = ',d_q_emiss(i,k)
                ELSE IF(d_q_emiss(i,k) < -1.E34) THEN
                   WRITE(*,*) 'IN STRATEMIT: d_q_emiss(i,k) < -1.E34 = ',d_q_emiss(i,k)
                ENDIF
                
                IF(flag_verbose_strataer) WRITE(*,*) 'IN STRATEMIT: ',&
    'i= ',i,'k= ',k, 'flh2o= ',flh2o, &
    'emission= ',emission, &
    'd_q_emiss(i,k)= ',d_q_emiss(i,k)
                
                IF(d_q_emiss(i,k) > 1.E34) THEN
                   d_q_emiss(i,k) = 1.E34
                   WRITE(*,*) 'IN STRATEMIT: d_q_emiss(i,k) > 1.E34 : ',d_q_emiss(i,k)
                ELSE IF(d_q_emiss(i,k) < -1.E34) THEN
                   d_q_emiss(i,k) = -1.E34
                   WRITE(*,*) 'IN STRATEMIT: d_q_emiss(i,k) < -1.E34 : ',d_q_emiss(i,k)
                ENDIF
             ENDIF
             budg_emi(i,ispeci)=budg_emi(i,ispeci)+emission*zdm(k)
          ENDDO
          
       ENDIF ! emission grid cell

    ENDDO ! klon loop

END SUBROUTINE STRATEMIT