SUBROUTINE macv2sp(pphis,pplay,paprs,xlon,xlat,tau_allaer,piz_allaer,cg_allaer,dNovrN)
  !
  !--routine to read the MACv2SP plume and compute optical properties
  !--requires flag_aerosol = 7
  !--feeds into aerosol optical properties and newmicro cloud droplet size if ok_cdnc activated
  !--for this one needs to feed natural (pre-industrial) aerosols twice for nat and 1980 files
  !--pre-ind aerosols (index=1) are not changed, present-day aerosols (index=2) are incremented
  !--uses model year so year_cur needs to be correct in the model simulation
  !
  !--aod_prof = AOD per layer 
  !--ssa_prof = SSA 
  !--asy_prof = asymetry parameter
  !--dNovrN   = enhancement factor for CDNC
  !
  !ASima : *_allaer arguments and non-SP "aer" diagnostics arrive from rrtm/aeropt_5wv_rrtm.F90 with values for inca1850 aerosols

  USE mo_simple_plumes, ONLY: sp_aod550_profile, sp_aop_profile, sp_setup       ! subroutines
  USE mo_simple_plumes, ONLY: nplumes, sp_initialized, aerosols_SP_forcing_year ! variables
  USE phys_cal_mod, ONLY : year_cur, year_len, days_elapsed
  USE dimphy
  USE aero_mod
  USE phys_local_var_mod, ONLY: t_seri, od443aer, od550aer, od865aer, ec550aer, dryod550aer, od550lt1aer
  ! SP-specific diagnostics (an SP extinction coefficient could be also defined: ec550SPaer = od550SPaer_prof/dz)
  USE phys_local_var_mod, ONLY: od550SPaer, od550SPaer_prof 
  !
  USE yomcst_mod_h
IMPLICIT NONE
  !

  ! 
  REAL,DIMENSION(klon),INTENT(IN)        :: pphis   ! Geopotentiel de surface
  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 les interfaces de chaque couche (en Pa)
  REAL,DIMENSION(klon),INTENT(IN)        :: xlat    ! latitudes pour chaque point
  REAL,DIMENSION(klon),INTENT(IN)        :: xlon    ! longitudes pour chaque point
  !
  REAL, DIMENSION(klon,klev,2,nbands_sw_rrtm), INTENT(INOUT) :: tau_allaer !  epaisseur optique aerosol
  REAL, DIMENSION(klon,klev,2,nbands_sw_rrtm), INTENT(INOUT) :: piz_allaer !  single scattering albedo aerosol
  REAL, DIMENSION(klon,klev,2,nbands_sw_rrtm), INTENT(INOUT) :: cg_allaer  !  asymmetry parameter aerosol
  !
  REAL,DIMENSION(klon),INTENT(OUT)       :: dNovrN   !< anthropogenic increase in cloud drop number concentration (factor)
  !
  REAL,DIMENSION(klon,klev,nplumes) :: aod550
  REAL,DIMENSION(klon,klev) :: aod_prof, ssa_prof, asy_prof
  REAL,DIMENSION(klon,klev) :: z, dz
  REAL,DIMENSION(klon)      :: oro, zrho, zt
  !
  REAL, PARAMETER    :: l443 = 443.0, l865 = 865.0 !--wavelengths in nm for diagnostics (other than 550 nm calculated by default)
  !
  INTEGER, PARAMETER :: Nwvmax=25
  REAL, DIMENSION(0:Nwvmax), PARAMETER :: lambda=(/ 240.0, &  !--this one is for band 1
                  280.0,  300.0,  330.0,  360.0,  400.0,   &  !--these are bounds of Streamer bands
                  440.0,  480.0,  520.0,  570.0,  640.0,   &
                  690.0,  750.0,  780.0,  870.0, 1000.0,   &
                 1100.0, 1190.0, 1280.0, 1530.0, 1640.0,   &
                 2130.0, 2380.0, 2910.0, 3420.0, 4000.0   /)
  !
  REAL, DIMENSION(1:Nwvmax-1), PARAMETER :: weight =(/    &   !--and the weights to be given to the bands
                 0.01,  4.05,  9.51, 15.99, 26.07, 33.10, &   !--corresponding to a typical solar spectrum 
                33.07, 39.91, 52.67, 27.89, 43.60, 13.67, &
                42.22, 40.12, 32.70, 14.44, 19.48, 14.23, &
                13.43, 16.42,  8.33,  0.95,  0.65,  2.76  /)
  !
  REAL :: zlambda, zweight
  REAL :: decimal_year 
  !
  INTEGER band, k, Nwv
  !
  ! define the height and dheight arrays
  !
  oro(:)  = pphis(:)/RG                             ! surface height in m
  ! ASima : the following calculations of z and dz are also done in phys_output_write_mod.F90
  !      for wind100m, under "IF (vars_defined)" ; "z" therein does not include "oro" 
  !    Layers' top and middle altitudes zhalf et zfull are also calculated for diagnostics
  !      in phys_output_write_mod.F90, under "IF (vars_defined)", using the geopotential "phi" ...   
  DO k = 1, klev
    zrho(:) = pplay(:,k)/t_seri(:,k)/RD                         ! air density in kg/m3
    dz(:,k) = (paprs(:,k)-paprs(:,k+1))/zrho(:)/RG              ! layer thickness in m
    IF (k==1) THEN
       z(:,1) = oro(:) + (paprs(:,1)-pplay(:,1))/zrho(:)/RG     ! altitude middle of first layer in m
       zt(:)  = oro(:) + dz(:,1)                                ! altitude top of first layer in m
    ELSE
      z(:,k) = zt(:) + (paprs(:,k)-pplay(:,k))/zrho(:)/RG       ! altitude middle of layer k in m
      zt(:)  = zt(:) + dz(:,k)                                  ! altitude top of layer k in m
    ENDIF
  ENDDO
  !

  ! Call sp_setup to retrieve 'aerosol_SP_forcing_year' value used in decimal_year calculation below
  IF (.NOT.sp_initialized) CALL sp_setup

  !--fractional year : 
  !   Original version, bugged :       year_fr = FLOAT(year_cur) + (FLOAT(day_cur)-0.5) / FLOAT(year_len)
  !   Correction 2026-03, FH & ASima : year_fr= float(year_cur) + float(days_elapsed)/float(year_len)
  ! Name changed in decimal_year ; choice between yearly vs fixed_year forcing, depending on 'aerosols_SP_forcing_year'
  IF (aerosols_SP_forcing_year.EQ.-9999) THEN
     decimal_year= float(year_cur) + float(days_elapsed)/float(year_len)
  ELSE
     decimal_year= float(aerosols_SP_forcing_year) + float(days_elapsed)/float(year_len)
  ENDIF

  IF (decimal_year.LT.1850.0.OR.decimal_year.GE.2024.0) THEN
     CALL abort_physic ('macv2sp','year not supported by plume model',1)
  ENDIF

  !
  !--call SUBROUTINE sp_aod550_profile once ; all the wavelength-dependent profiles and 2D diagnostics will use aod550 profile
  !
  CALL sp_aod550_profile                                            ( &
       klev           ,klon       ,xlon       ,xlat      , &
       decimal_year   ,z          ,dz         ,dNovrN    , aod550   )

  !--AOD550 calculations for diagnostics
  !
  ! aod550 profile = sum over the 'nplumes' dimension
  od550SPaer_prof(:,:)=SUM(aod550(:,:,:),dim=3)
  !
  !--a/ AOD of SP anthropic aerosols = vertical sum of SP aod profile
  od550SPaer(:)=SUM(od550SPaer_prof(:,:),dim=2)
  !
  !--b/ total AOD = Inca1850 + SP
  od550aer(:)=od550aer(:)+od550SPaer(:)
  !
  !--c/ fine-mode AOD = Inca1850(fine mode) + SP
  od550lt1aer(:)=od550lt1aer(:)+od550SPaer(:)
  !
  !--d/ dry AOD 
  dryod550aer(:)=dryod550aer(:)+od550SPaer(:)
  !
  !--extinction coefficient for diagnostic
  ec550aer(:,:)=ec550aer(:,:)+od550SPaer_prof(:,:)/dz(:,:)
  !

  ! Diagnostic AOD calculations for other wavelengths : 443 and 865 nm  
  !--call to sp routine -- 443 nm
  !
  CALL sp_aop_profile                                    ( &
       klev     ,klon ,  l443         , &
       aod550   ,aod_prof ,ssa_prof  , &
       asy_prof )

  od443aer(:)= od443aer(:)+SUM(aod_prof(:,:),dim=2)

  !--call to sp routine -- 865 nm
  !
  CALL sp_aop_profile                                    ( &
       klev     ,klon ,  l865         , &
       aod550   ,aod_prof ,ssa_prof  , &
       asy_prof )

  od865aer(:)=od865aer(:)+SUM(aod_prof(:,:),dim=2)

  !
  !--re-weighting of piz and cg arrays before adding the anthropogenic aerosols
  !--index 2 = all natural + anthropogenic aerosols
  piz_allaer(:,:,2,:)=piz_allaer(:,:,2,:)*tau_allaer(:,:,2,:)
  cg_allaer(:,:,2,:) =cg_allaer(:,:,2,:)*piz_allaer(:,:,2,:)
  !
  !--now computing the same at many wavelengths to fill the model bands
  !
  DO Nwv=0,Nwvmax-1

    IF (Nwv.EQ.0) THEN          !--RRTM spectral band 1
      zlambda=lambda(Nwv)
      zweight=1.0
      band=1
    ELSEIF (Nwv.LE.5) THEN      !--RRTM spectral band 2
      zlambda=0.5*(lambda(Nwv)+lambda(Nwv+1))
      zweight=weight(Nwv)/SUM(weight(1:5))
      band=2
    ELSEIF (Nwv.LE.10) THEN     !--RRTM spectral band 3
      zlambda=0.5*(lambda(Nwv)+lambda(Nwv+1))
      zweight=weight(Nwv)/SUM(weight(6:10))
      band=3
    ELSEIF (Nwv.LE.16) THEN     !--RRTM spectral band 4
      zlambda=0.5*(lambda(Nwv)+lambda(Nwv+1))
      zweight=weight(Nwv)/SUM(weight(11:16))
      band=4
    ELSEIF (Nwv.LE.21) THEN     !--RRTM spectral band 5
      zlambda=0.5*(lambda(Nwv)+lambda(Nwv+1))
      zweight=weight(Nwv)/SUM(weight(17:21))
      band=5
    ELSE                        !--RRTM spectral band 6
      zlambda=0.5*(lambda(Nwv)+lambda(Nwv+1))
      zweight=weight(Nwv)/SUM(weight(22:Nwvmax-1))
      band=6
    ENDIF
    !
  CALL sp_aop_profile                                    ( &
       klev     ,klon ,  zlambda         , &
       aod550   ,aod_prof ,ssa_prof  , &
       asy_prof )


    !
    !--adding up the quantities tau, piz*tau and cg*piz*tau
    tau_allaer(:,:,2,band)=tau_allaer(:,:,2,band)+zweight*MAX(aod_prof(:,:),1.e-15)
    piz_allaer(:,:,2,band)=piz_allaer(:,:,2,band)+zweight*MAX(aod_prof(:,:),1.e-15)*ssa_prof(:,:)
    cg_allaer(:,:,2,band) =cg_allaer(:,:,2,band) +zweight*MAX(aod_prof(:,:),1.e-15)*ssa_prof(:,:)*asy_prof(:,:)
    !
  ENDDO
  !
  !--renormalizing cg and piz now that MACv2SP increments have been added
  cg_allaer(:,:,2,:) =cg_allaer(:,:,2,:) /piz_allaer(:,:,2,:)
  piz_allaer(:,:,2,:)=piz_allaer(:,:,2,:)/tau_allaer(:,:,2,:)
  !
END SUBROUTINE macv2sp 
