MODULE strataer_emiss_mod
  ! This module contains information about strato microphysic model emission parameters
  IMPLICIT NONE; PRIVATE
  PUBLIC strataer_emiss_init

CONTAINS

  SUBROUTINE strataer_emiss_init()

    USE strataer_local_var_mod
    USE lmdz_ioipsl_getin_p, ONLY: getin_p
    USE lmdz_print_control, ONLY: lunout
    USE lmdz_phys_para, ONLY: is_master
    USE lmdz_abort_physic, ONLY: abort_physic
    USE dimphy, ONLY: klon

    IMPLICIT NONE

    ! Local variables
    INTEGER :: ispec

    WRITE(lunout, *) 'IN STRATAER_EMISS INIT WELCOME!'

    IF (flag_emit==1 .OR. flag_emit==4) THEN ! Volcano
      CALL getin_p('nErupt', nErupt) !eruption nb
      CALL getin_p('injdur', injdur) !injection duration

      IF (flag_emit==1) THEN
        CALL getin_p('nAerErupt', nAerErupt) !sulfur aer nb
      ELSEIF (flag_emit==4) THEN
        CALL getin_p('nSpeciesErupt', nSpeciesErupt) !chimical species nb
      ENDIF

      IF (nErupt>0) THEN
        ALLOCATE(year_emit_vol(nErupt), mth_emit_vol(nErupt), day_emit_vol(nErupt))
        year_emit_vol = 0 ; mth_emit_vol = 0 ; day_emit_vol = 0
        ALLOCATE(altemiss_vol(nErupt), sigma_alt_vol(nErupt))
        ALLOCATE(xlat_min_vol(nErupt), xlon_min_vol(nErupt))
        ALLOCATE(xlat_max_vol(nErupt), xlon_max_vol(nErupt))
        altemiss_vol = 0. ; sigma_alt_vol = 0.
        xlon_min_vol = 0. ; xlon_max_vol = 0.
        xlat_min_vol = 0. ; xlat_max_vol = 0.
        ! injection params (dates, loc, injections params)
        CALL getin_p('year_emit_vol', year_emit_vol)
        CALL getin_p('mth_emit_vol', mth_emit_vol)
        CALL getin_p('day_emit_vol', day_emit_vol)
        CALL getin_p('altemiss_vol', altemiss_vol)
        CALL getin_p('sigma_alt_vol', sigma_alt_vol)
        CALL getin_p('xlon_min_vol', xlon_min_vol)
        CALL getin_p('xlon_max_vol', xlon_max_vol)
        CALL getin_p('xlat_min_vol', xlat_min_vol)
        CALL getin_p('xlat_max_vol', xlat_max_vol)
        IF (flag_emit==1) THEN
          ALLOCATE(m_sulf_emiss_vol(nErupt))
          ALLOCATE(m_aer_emiss_vol(nErupt, nAerErupt))
          m_aer_emiss_vol = 0. ; m_sulf_emiss_vol = 0.
          IF (ok_qemiss) THEN
            ALLOCATE(m_H2O_emiss_vol(nErupt))
            ALLOCATE(m_H2O_emiss_vol_daily(nErupt))
            !                ALLOCATE(d_q_emiss(klon,klev))
            ALLOCATE(budg_emi(klon, nAerErupt + 1))
            m_H2O_emiss_vol(:) = 0.
            m_H2O_emiss_vol_daily(:) = 0.
            !                d_q_emiss(:,:)=0.
          ELSE
            ALLOCATE(budg_emi(klon, nAerErupt))
          ENDIF
          budg_emi(:, :) = 0.
        ELSEIF (flag_emit==4) THEN
          ALLOCATE(m_Chlore_emiss_vol(nErupt))
          ALLOCATE(id_HCl)
          ALLOCATE(m_Brome_emiss_vol(nErupt))
          ALLOCATE(id_HBr)
          ALLOCATE(id_species(nSpeciesErupt))
          id_species = 0
          ALLOCATE(m_species_emiss_vol(nErupt, nSpeciesErupt))
          m_species_emiss_vol = 0.
          ALLOCATE(m_NOx_emiss_vol(nErupt))
          ALLOCATE(m_H2O_emiss_vol(nErupt))
          m_Chlore_emiss_vol = 0. ; m_Brome_emiss_vol = 0.
          m_NOx_emiss_vol = 0. ; m_H2O_emiss_vol = 0.
          ALLOCATE(id_NOx)
          ALLOCATE(id_H2O)
          ALLOCATE(budg_emi(klon, nSpeciesErupt))
          budg_emi(:, :) = 0.
        ENDIF
      ELSE
        WRITE(lunout, *) 'ERROR : Using flag_emit=1 or 4 (ie Volcanic eruption) but nErupt (', nErupt, ') <=0 !'
        CALL abort_physic('strataer_emiss_mod', &
                'No eruption define in physiq_def (nErupt=0). Add at one eruption or use background condition.', 1)
      ENDIF ! fin if nerupt

      IF (flag_emit==1) THEN
        CALL getin_p('m_sulf_emiss_vol', m_sulf_emiss_vol)
        IF (ok_qemiss) THEN
          CALL getin_p('m_H2O_emiss_vol', m_H2O_emiss_vol)
        endif
      ELSEIF (flag_emit==4) THEN
        CALL getin_p('id_species', id_species)
        CALL getin_p('m_Chlore_emiss_vol', m_Chlore_emiss_vol)
        CALL getin_p('id_HCl', id_HCl)
        CALL getin_p('m_Brome_emiss_vol', m_Brome_emiss_vol)
        CALL getin_p('id_HBr', id_HBr)
        CALL getin_p('m_NOx_emiss_vol', m_NOx_emiss_vol)
        CALL getin_p('id_NOx', id_NOx)
        CALL getin_p('m_H2O_emiss_vol', m_H2O_emiss_vol)
        CALL getin_p('id_H2O', id_H2O)
      ENDIF

    ELSEIF (flag_emit == 2) THEN ! SAI
      CALL getin_p('m_aer_emiss_sai', m_aer_emiss_sai)
      CALL getin_p('altemiss_sai', altemiss_sai)
      CALL getin_p('sigma_alt_sai', sigma_alt_sai)
      CALL getin_p('xlat_sai', xlat_sai)
      CALL getin_p('xlon_sai', xlon_sai)
      CALL getin_p('year_emit_sai_start', year_emit_sai_start)
      CALL getin_p('year_emit_sai_end', year_emit_sai_end)
      CALL getin_p('mth_emit_sai_start', mth_emit_sai_start)
      CALL getin_p('mth_emit_sai_end', mth_emit_sai_end)
      CALL getin_p('day_emit_sai_start', day_emit_sai_start)
      CALL getin_p('day_emit_sai_end', day_emit_sai_end)

    ELSEIF (flag_emit == 3) THEN ! SAI between latitudes
      CALL getin_p('m_aer_emiss_sai', m_aer_emiss_sai)
      CALL getin_p('altemiss_sai', altemiss_sai)
      CALL getin_p('sigma_alt_sai', sigma_alt_sai)
      CALL getin_p('xlon_sai', xlon_sai)
      CALL getin_p('xlat_max_sai', xlat_max_sai)
      CALL getin_p('xlat_min_sai', xlat_min_sai)
    ENDIF

    IF (flag_emit == 1) THEN
      DO ispec = 1, nAerErupt
        m_aer_emiss_vol(:, ispec) = m_sulf_emiss_vol(:)
      ENDDO
    ELSEIF (flag_emit== 4) THEN
      DO ispec = 1, nSpeciesErupt
        IF(id_species(ispec) == id_HCl) THEN
          m_species_emiss_vol(:, ispec) = m_Chlore_emiss_vol(:)
        ENDIF
        IF (id_species(ispec) == id_HBr) THEN
          m_species_emiss_vol(:, ispec) = m_Brome_emiss_vol(:)
        ENDIF
        IF (id_species(ispec) == id_NOx) THEN
          m_species_emiss_vol(:, ispec) = m_NOx_emiss_vol(:)
        ENDIF
        IF (id_species(ispec) == id_H2O) THEN
          m_species_emiss_vol(:, ispec) = m_H2O_emiss_vol(:)
        ENDIF
      ENDDO
    ENDIF

    !============= Injection ponderation =============
    IF (flag_emit > 0) THEN
      CALL strataer_ponde_init
      WRITE(lunout, *) 'IN STRATAER INIT : ponde_lonlat_vol', ponde_lonlat_vol
    ENDIF

    !============= Print params =============
    IF (is_master) THEN
      IF (nErupt > 0) THEN
        IF (flag_emit == 1 .OR. flag_emit == 4) THEN
          WRITE(lunout, *) 'IN STRATAER nErupt: ', nErupt
          WRITE(lunout, *) 'IN STRATAER injdur: ', injdur
          WRITE(lunout, *) 'IN STRATAER nAerErupt: ', nAerErupt

          WRITE(lunout, *) 'IN STRATAER : year_emit_vol', year_emit_vol
          WRITE(lunout, *) 'IN STRATAER : mth_emit_vol', mth_emit_vol
          WRITE(lunout, *) 'IN STRATAER : day_emit_vol', day_emit_vol
          WRITE(lunout, *) 'IN STRATAER : altemiss_vol', altemiss_vol
          WRITE(lunout, *) 'IN STRATAER : sigma_alt_vol', sigma_alt_vol
          WRITE(lunout, *) 'IN STRATAER : xlat_min_vol', xlat_min_vol
          WRITE(lunout, *) 'IN STRATAER : xlat_max_vol', xlat_max_vol
          WRITE(lunout, *) 'IN STRATAER : xlon_min_vol', xlon_min_vol
          WRITE(lunout, *) 'IN STRATAER : xlon_max_vol', xlon_max_vol
          IF (flag_emit==1) THEN
            WRITE(lunout, *) 'IN STRATAEREMISS : m_sulf_emiss_vol', m_sulf_emiss_vol
            WRITE(lunout, *) 'IN STRATAER : m_aer_emiss_vol', m_aer_emiss_vol
            IF (ok_qemiss) THEN
              WRITE(lunout, *) 'IN STRATAEREMISS : m_H2O_emiss_vol', m_H2O_emiss_vol
            ENDIF
          ENDIF
        ELSEIF (flag_emit == 2) THEN
          WRITE(lunout, *) 'IN STRATAER : m_aer_emiss_sai', m_aer_emiss_sai
          WRITE(lunout, *) 'IN STRATAER : altemiss_sai', altemiss_sai
          WRITE(lunout, *) 'IN STRATAER : sigma_alt_sai', sigma_alt_sai
          WRITE(lunout, *) 'IN STRATAER : xlat_sai', xlat_sai
          WRITE(lunout, *) 'IN STRATAER : xlon_sai', xlon_sai
        ELSEIF (flag_emit == 3) THEN
          WRITE(lunout, *) 'IN STRATAER : m_aer_emiss_sai', m_aer_emiss_sai
          WRITE(lunout, *) 'IN STRATAER : altemiss_sai', altemiss_sai
          WRITE(lunout, *) 'IN STRATAER : sigma_alt_sai', sigma_alt_sai
          WRITE(lunout, *) 'IN STRATAER : year_emit_sai start/end', year_emit_sai_start, year_emit_sai_end
          WRITE(lunout, *) 'IN STRATAER : mth_emit_sai start/end', mth_emit_sai_start, mth_emit_sai_end
          WRITE(lunout, *) 'IN STRATAER : day_emit_sai start/end', day_emit_sai_start, day_emit_sai_end
          WRITE(lunout, *) 'IN STRATAER : xlat_min_sai', xlat_min_sai
          WRITE(lunout, *) 'IN STRATAER : xlat_max_sai', xlat_max_sai
          WRITE(lunout, *) 'IN STRATAER : xlon_sai', xlon_sai
        ENDIF
        IF(flag_emit == 4) THEN
          WRITE(lunout, *) 'IN STRATAER : nSpeciesErupt: ', nSpeciesErupt
          WRITE(lunout, *) 'IN STRATAER : id_species = ', id_species
          WRITE(lunout, *) 'IN STRATAER : id_HCl = ', id_HCl
          WRITE(lunout, *) 'IN STRATAER : id_HBr = ', id_HBr
          WRITE(lunout, *) 'IN STRATAER : id_NOx = ', id_NOx
          WRITE(lunout, *) 'IN STRATAER : id_H2O = ', id_H2O
          WRITE(lunout, *) 'IN STRATAER : m_Chlore_emiss_vol = ', m_Chlore_emiss_vol
          WRITE(lunout, *) 'IN STRATAER : m_Brome_emiss_vol = ', m_Brome_emiss_vol
          WRITE(lunout, *) 'IN STRATAER : m_NOx_emiss_vol = ', m_NOx_emiss_vol
          WRITE(lunout, *) 'IN STRATAER : m_H2O_emiss_vol = ', m_H2O_emiss_vol
        ENDIF
      endif
    ENDIF ! if master

    WRITE(lunout, *) 'IN STRATAER_EMISS END'
  END SUBROUTINE strataer_emiss_init

  ! Compute the ponderation to applicate in each grid point for all eruptions and init
  ! dlat & dlon variables
  SUBROUTINE strataer_ponde_init()

    USE lmdz_regular_lonlat, ONLY: lon_reg, lat_reg
    USE dimphy, ONLY: klon
    USE lmdz_grid_phy, ONLY: nbp_lat, nbp_lon
    USE lmdz_print_control, ONLY: lunout
    USE strataer_local_var_mod
    USE lmdz_yomcst, ONLY: RPI

    ! local var
    REAL :: lat_reg_deg, lon_reg_deg ! lat and lon of grid points in degree
    INTEGER :: ieru, i, j

    ALLOCATE(ponde_lonlat_vol(nErupt))

    !Compute lon/lat ponderation for injection
    dlat = 180. / 2. / FLOAT(nbp_lat)   ! d latitude in degree
    dlon = 360. / 2. / FLOAT(nbp_lon)   ! d longitude in degree
    WRITE(lunout, *) 'IN STRATAER_INIT dlat=', dlat, 'dlon=', dlon
    WRITE(lunout, *) 'IN STRATAER_INIT nErupt=', nErupt
    WRITE(lunout, *) 'IN STRATAER_INIT xlat_min=', xlat_min_vol, 'xlat_max=', xlat_max_vol
    WRITE(lunout, *) 'IN STRATAER_INIT xlon_min=', xlon_min_vol, 'xlon_max=', xlon_max_vol

    DO ieru = 1, nErupt
      ponde_lonlat_vol(ieru) = 0
      DO i = 1, nbp_lon
        lon_reg_deg = lon_reg(i) * 180. / RPI
        DO j = 1, nbp_lat
          lat_reg_deg = lat_reg(j) * 180. / RPI
          IF  (lat_reg_deg>=xlat_min_vol(ieru) - dlat .AND. lat_reg_deg<xlat_max_vol(ieru) + dlat .AND. &
                  lon_reg_deg>=xlon_min_vol(ieru) - dlon .AND. lon_reg_deg<xlon_max_vol(ieru) + dlon) THEN
            ponde_lonlat_vol(ieru) = ponde_lonlat_vol(ieru) + 1
          ENDIF
        ENDDO
      ENDDO
      IF(ponde_lonlat_vol(ieru) == 0) THEN
        WRITE(lunout, *) 'STRATAER_INIT ERROR: no grid point found for eruption ieru=', ieru
      ENDIF
    ENDDO !ieru

    WRITE(lunout, *) 'IN STRATAER_PONDE_INIT ponde_lonlat: ', ponde_lonlat_vol

  END SUBROUTINE strataer_ponde_init

END MODULE strataer_emiss_mod