Index: LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/conf_phys_m.F90
===================================================================
--- LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/conf_phys_m.F90	(revision 3314)
+++ LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/conf_phys_m.F90	(revision 3318)
@@ -414,4 +414,5 @@
     ! - flag_aerosol=5 => dust only
     ! - flag_aerosol=6 => all aerosol
+    ! - flag_aerosol=7 => natural aerosol + MACv2SP
 
     flag_aerosol_omp = 0
@@ -2448,10 +2449,10 @@
 
     ! Flag_aerosol cannot be to zero if we are in coupled mode for aerosol
-    IF (aerosol_couple .AND. flag_aerosol .eq. 0 ) THEN
+    IF (aerosol_couple .AND. flag_aerosol .EQ. 0 ) THEN
        CALL abort_physic('conf_phys', 'flag_aerosol cannot be to zero if aerosol_couple=y ', 1)
     ENDIF
 
     ! flag_aerosol need to be different to zero if ok_cdnc is activated
-    IF (ok_cdnc .AND. flag_aerosol .eq. 0) THEN
+    IF (ok_cdnc .AND. flag_aerosol .EQ. 0) THEN
        CALL abort_physic('conf_phys', 'flag_aerosol cannot be to zero if ok_cdnc is activated ', 1)
     ENDIF
@@ -2460,4 +2461,12 @@
     IF (ok_aie .AND. .NOT. ok_cdnc) THEN
        CALL abort_physic('conf_phys', 'ok_cdnc must be set to y if ok_aie is activated',1)
+    ENDIF
+
+    ! flag_aerosol=7 => MACv2SP climatology 
+    IF (flag_aerosol.EQ.7.AND. iflag_rrtm.NE.1) THEN
+       CALL abort_physic('conf_phys', 'flag_aerosol=7 (MACv2SP) can only be activated with RRTM',1)
+    ENDIF
+    IF (flag_aerosol.EQ.7.AND. NSW.NE.6) THEN
+       CALL abort_physic('conf_phys', 'flag_aerosol=7 (MACv2SP) can only be activated with NSW=6',1)
     ENDIF
 
Index: LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/macv2sp.F90
===================================================================
--- LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/macv2sp.F90	(revision 3318)
+++ LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/macv2sp.F90	(revision 3318)
@@ -0,0 +1,174 @@
+SUBROUTINE MACv2SP(pphis,pplay,paprs,xlon,xlat,tau_allaer,piz_allaer,cg_allaer)
+  !
+  !--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
+  !
+  USE mo_simple_plumes, ONLY: sp_aop_profile
+  USE phys_cal_mod, ONLY : year_cur, day_cur, year_len
+  USE dimphy
+  USE aero_mod
+  USE phys_local_var_mod, ONLY: t_seri, od443aer, od550aer, od865aer, ec550aer, dryod550aer, od550lt1aer, dNovrN
+  !!USE YOMCST, ONLY : RD, RG
+  !
+  IMPLICIT NONE
+  !
+  include "YOMCST.h"
+  ! 
+  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(OUT) :: tau_allaer !  epaisseur optique aerosol
+  REAL, DIMENSION(klon,klev,2,nbands_sw_rrtm), INTENT(OUT) :: piz_allaer !  single scattering albedo aerosol
+  REAL, DIMENSION(klon,klev,2,nbands_sw_rrtm), INTENT(OUT) :: cg_allaer  !  asymmetry parameter aerosol
+  !
+  REAL,DIMENSION(klon,klev) :: aod_prof, ssa_prof, asy_prof
+  REAL,DIMENSION(klon,klev) :: z, dz
+  REAL,DIMENSION(klon)      :: oro, zrho, zt
+  !
+  INTEGER, PARAMETER :: nmon = 12
+  !
+  REAL, PARAMETER    :: l443 = 443.0, l550 = 550.0, l865 = 865.0 !--wavelengths in nm
+  !
+  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 :: year_fr
+  !
+  INTEGER band, i, k, Nwv
+  !
+  ! define the height and dheight arrays
+  !
+  oro(:)  = pphis(:)/RG                             ! surface height in m
+  !
+  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
+  !
+  !--fractional year
+  !
+  year_fr = FLOAT(year_cur) + (FLOAT(day_cur)-0.5) / FLOAT(year_len)
+  print *,'year_fr=',year_fr
+  !
+  !--call to sp routine -- 443 nm
+  !
+  CALL sp_aop_profile                                    ( &
+       klev     ,klon ,l443 ,oro    ,xlon     ,xlat      , &
+       year_fr  ,z    ,dz   ,dNovrN ,aod_prof ,ssa_prof  , &
+       asy_prof )
+  !
+  !--AOD calculations for diagnostics
+  od443aer(:)= od443aer(:)+SUM(aod_prof(:,:),dim=2)
+  !
+  !--call to sp routine -- 550 nm
+  !
+  CALL sp_aop_profile                                    ( &
+       klev     ,klon ,l550 ,oro    ,xlon     ,xlat      , &
+       year_fr  ,z    ,dz   ,dNovrN ,aod_prof ,ssa_prof  , &
+       asy_prof )
+  !
+  !--AOD calculations for diagnostics
+  od550aer(:)=od550aer(:)+SUM(aod_prof(:,:),dim=2)
+  !
+  !--dry AOD calculation for diagnostics
+  dryod550aer(:)=dryod550aer(:)+od550aer(:)
+  !
+  !--fine-mode AOD calculation for diagnostics
+  od550lt1aer(:)=od550lt1aer(:)+od550aer(:)
+  !
+  !--extinction coefficient for diagnostic
+  ec550aer(:,:)=ec550aer(:,:)+aod_prof(:,:)/dz(:,:)
+  !
+  !--call to sp routine -- 865 nm
+  !
+  CALL sp_aop_profile                                    ( &
+       klev     ,klon ,l865 ,oro    ,xlon     ,xlat      , &
+       year_fr  ,z    ,dz   ,dNovrN ,aod_prof ,ssa_prof  , &
+       asy_prof )
+  !
+  !--AOD calculations for diagnostics
+  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 ,oro    ,xlon     ,xlat      , &
+         year_fr  ,z    ,dz      ,dNovrN ,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
+  !
+  !--renpomalizing 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
Index: LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/mo_simple_plumes.F90
===================================================================
--- LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/mo_simple_plumes.F90	(revision 3318)
+++ LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/mo_simple_plumes.F90	(revision 3318)
@@ -0,0 +1,413 @@
+!>
+!!
+!! @brief Module MO_SIMPLE_PLUMES: provides anthropogenic aerosol optical properties as a function of lat, lon
+!!   height, time, and wavelength
+!!
+!! @remarks
+!!
+!! @author Bjorn Stevens, Stephanie Fiedler and Karsten Peters MPI-Met, Hamburg (v1 release 2016-11-10)
+!!
+!! @change-log:
+!!          - 2016-12-05: beta release (BS, SF and KP, MPI-Met)
+!!          - 2016-09-28: revised representation of Twomey effect (SF, MPI-Met)
+!!          - 2015-09-28: bug fixes  (SF, MPI-Met)
+!!          - 2016-10-12: revised maximum longitudinal extent of European plume (KP, SF, MPI-Met)
+!! $ID: n/a$
+!!
+!! @par Origin
+!!   Based on code originally developed at the MPI-Met by Karsten Peters, Bjorn Stevens, Stephanie Fiedler
+!!   and Stefan Kinne with input from Thorsten Mauritsen and Robert Pincus
+!!
+!! @par Copyright
+!! 
+!
+MODULE MO_SIMPLE_PLUMES
+
+  USE netcdf
+
+  IMPLICIT NONE
+
+  INTEGER, PARAMETER ::                        &
+       nplumes   = 9                          ,& !< Number of plumes
+       nfeatures = 2                          ,& !< Number of features per plume
+       ntimes    = 52                         ,& !< Number of times resolved per year (52 => weekly resolution)
+       nyears    = 251                           !< Number of years of available forcing
+
+  LOGICAL, SAVE ::                             &
+       sp_initialized = .FALSE.                  !< parameter determining whether input needs to be read
+!$OMP THREADPRIVATE(sp_initialized)
+
+  REAL ::                                      &
+       plume_lat      (nplumes)               ,& !< latitude of plume center (AOD maximum)
+       plume_lon      (nplumes)               ,& !< longitude of plume center (AOD maximum)
+       beta_a         (nplumes)               ,& !< parameter a for beta function vertical profile
+       beta_b         (nplumes)               ,& !< parameter b for beta function vertical profile
+       aod_spmx       (nplumes)               ,& !< anthropogenic AOD maximum at 550 for plumes 
+       aod_fmbg       (nplumes)               ,& !< anthropogenic AOD at 550 for fine-mode natural background (idealized to mimic Twomey effect)
+       asy550         (nplumes)               ,& !< asymmetry parameter at 550nm for plume
+       ssa550         (nplumes)               ,& !< single scattering albedo at 550nm for plume
+       angstrom       (nplumes)               ,& !< Angstrom parameter for plume 
+       sig_lon_E      (nfeatures,nplumes)     ,& !< Eastward extent of plume feature
+       sig_lon_W      (nfeatures,nplumes)     ,& !< Westward extent of plume feature
+       sig_lat_E      (nfeatures,nplumes)     ,& !< Southward extent of plume feature
+       sig_lat_W      (nfeatures,nplumes)     ,& !< Northward extent of plume feature
+       theta          (nfeatures,nplumes)     ,& !< Rotation angle of plume feature
+       ftr_weight     (nfeatures,nplumes)     ,& !< Feature weights 
+       time_weight    (nfeatures,nplumes)     ,& !< Time weights 
+       time_weight_bg (nfeatures,nplumes)     ,& !< as time_weight but for natural background in Twomey effect 
+       year_weight    (nyears,nplumes)        ,& !< Yearly weight for plume
+       ann_cycle      (nfeatures,ntimes,nplumes) !< annual cycle for plume feature
+!$OMP THREADPRIVATE(plume_lat,plume_lon,beta_a,beta_b,aod_spmx,aod_fmbg,asy550,ssa550,angstrom)
+!$OMP THREADPRIVATE(sig_lon_E,sig_lon_W,sig_lat_E,sig_lat_W,theta,ftr_weight,year_weight,ann_cycle)
+
+  PUBLIC sp_aop_profile
+
+CONTAINS
+  !
+  ! ------------------------------------------------------------------------------------------------------------------------
+  ! SP_SETUP:  This subroutine should be called at initialization to read the netcdf data that describes the simple plume
+  ! climatology.  The information needs to be either read by each processor or distributed to processors.
+  !
+  SUBROUTINE sp_setup
+    !
+    USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root
+    USE mod_phys_lmdz_omp_data, ONLY: is_omp_root
+    USE mod_phys_lmdz_transfert_para, ONLY: bcast
+    !
+    ! ---------- 
+    !
+    INTEGER :: iret, ncid, DimID, VarID, xdmy
+    !
+    ! ---------- 
+    !--only one processor reads the input data
+    IF (is_mpi_root.AND.is_omp_root) THEN
+       !    
+       iret = nf90_open("MACv2.0-SP_v1.nc", NF90_NOWRITE, ncid)
+       IF (iret /= NF90_NOERR) STOP 'NetCDF File not opened'
+       !
+       ! read dimensions and make sure file conforms to expected size
+       !
+       iret = nf90_inq_dimid(ncid, "plume_number"  , DimId)
+       iret = nf90_inquire_dimension(ncid, DimId, len = xdmy)
+       IF (xdmy /= nplumes) STOP 'NetCDF improperly dimensioned -- plume_number'
+       !
+       iret = nf90_inq_dimid(ncid, "plume_feature", DimId)
+       iret = nf90_inquire_dimension(ncid, DimId, len = xdmy)
+       IF (xdmy /= nfeatures) STOP 'NetCDF improperly dimensioned -- plume_feature'
+       !
+       iret = nf90_inq_dimid(ncid, "year_fr"   , DimId)
+       iret = nf90_inquire_dimension(ncid, DimID, len = xdmy)
+       IF (xdmy /= ntimes) STOP 'NetCDF improperly dimensioned -- year_fr'
+       !
+       iret = nf90_inq_dimid(ncid, "years"   , DimId)
+       iret = nf90_inquire_dimension(ncid, DimID, len = xdmy)
+       IF (xdmy /= nyears) STOP 'NetCDF improperly dimensioned -- years'
+       !
+       ! read variables that define the simple plume climatology
+       !
+       iret = nf90_inq_varid(ncid, "plume_lat", VarId)
+       iret = nf90_get_var(ncid, VarID, plume_lat(:), start=(/1/),count=(/nplumes/))
+       IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading plume_lat'
+       iret = nf90_inq_varid(ncid, "plume_lon", VarId)
+       iret = nf90_get_var(ncid, VarID, plume_lon(:), start=(/1/),count=(/nplumes/))
+       IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading plume_lon'
+       iret = nf90_inq_varid(ncid, "beta_a"   , VarId)
+       iret = nf90_get_var(ncid, VarID, beta_a(:)   , start=(/1/),count=(/nplumes/))
+       IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading beta_a'
+       iret = nf90_inq_varid(ncid, "beta_b"   , VarId)
+       iret = nf90_get_var(ncid, VarID, beta_b(:)   , start=(/1/),count=(/nplumes/))
+       IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading beta_b'
+       iret = nf90_inq_varid(ncid, "aod_spmx" , VarId)
+       iret = nf90_get_var(ncid, VarID, aod_spmx(:)  , start=(/1/),count=(/nplumes/))
+       IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading aod_spmx'
+       iret = nf90_inq_varid(ncid, "aod_fmbg" , VarId)
+       iret = nf90_get_var(ncid, VarID, aod_fmbg(:)  , start=(/1/),count=(/nplumes/))
+       IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading aod_fmbg'
+       iret = nf90_inq_varid(ncid, "ssa550"   , VarId)
+       iret = nf90_get_var(ncid, VarID, ssa550(:)  , start=(/1/),count=(/nplumes/))
+       IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading ssa550'
+       iret = nf90_inq_varid(ncid, "asy550"   , VarId)
+       iret = nf90_get_var(ncid, VarID, asy550(:)  , start=(/1/),count=(/nplumes/))
+       IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading asy550'
+       iret = nf90_inq_varid(ncid, "angstrom" , VarId)
+       iret = nf90_get_var(ncid, VarID, angstrom(:), start=(/1/),count=(/nplumes/))
+       IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading angstrom'
+       !
+       iret = nf90_inq_varid(ncid, "sig_lat_W"     , VarId)
+       iret = nf90_get_var(ncid, VarID, sig_lat_W(:,:)    , start=(/1,1/),count=(/nfeatures,nplumes/))
+       IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading sig_lat_W'
+       iret = nf90_inq_varid(ncid, "sig_lat_E"     , VarId)
+       iret = nf90_get_var(ncid, VarID, sig_lat_E(:,:)    , start=(/1,1/),count=(/nfeatures,nplumes/))
+       IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading sig_lat_E'
+       iret = nf90_inq_varid(ncid, "sig_lon_E"     , VarId)
+       iret = nf90_get_var(ncid, VarID, sig_lon_E(:,:)    , start=(/1,1/),count=(/nfeatures,nplumes/))
+       IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading sig_lon_E'
+       iret = nf90_inq_varid(ncid, "sig_lon_W"     , VarId)
+       iret = nf90_get_var(ncid, VarID, sig_lon_W(:,:)    , start=(/1,1/),count=(/nfeatures,nplumes/))
+       IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading sig_lon_W'
+       iret = nf90_inq_varid(ncid, "theta"         , VarId)
+       iret = nf90_get_var(ncid, VarID, theta(:,:)        , start=(/1,1/),count=(/nfeatures,nplumes/))
+       IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading theta'
+       iret = nf90_inq_varid(ncid, "ftr_weight"    , VarId)
+       iret = nf90_get_var(ncid, VarID, ftr_weight(:,:)   , start=(/1,1/),count=(/nfeatures,nplumes/))
+       IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading plume_lat'
+       iret = nf90_inq_varid(ncid, "year_weight"   , VarId)
+       iret = nf90_get_var(ncid, VarID, year_weight(:,:)  , start=(/1,1/),count=(/nyears,nplumes   /))
+       IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading year_weight'
+       iret = nf90_inq_varid(ncid, "ann_cycle"     , VarId)
+       iret = nf90_get_var(ncid, VarID, ann_cycle(:,:,:)  , start=(/1,1,1/),count=(/nfeatures,ntimes,nplumes/))
+       IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading ann_cycle'
+       !
+       iret = nf90_close(ncid)
+       !
+    ENDIF !--root processor
+    !
+    CALL bcast(plume_lat)
+    CALL bcast(plume_lon)
+    CALL bcast(beta_a)
+    CALL bcast(beta_b)
+    CALL bcast(aod_spmx)
+    CALL bcast(aod_fmbg)
+    CALL bcast(asy550)
+    CALL bcast(ssa550)
+    CALL bcast(angstrom)
+    CALL bcast(sig_lon_E)
+    CALL bcast(sig_lon_W)
+    CALL bcast(sig_lat_E)
+    CALL bcast(sig_lat_W)
+    CALL bcast(theta)
+    CALL bcast(ftr_weight)
+    CALL bcast(year_weight)
+    CALL bcast(ann_cycle)
+    !
+    sp_initialized = .TRUE.
+    !
+    RETURN
+    !
+  END SUBROUTINE sp_setup
+  !
+  ! ------------------------------------------------------------------------------------------------------------------------
+  ! SET_TIME_WEIGHT:  The simple plume model assumes that meteorology constrains plume shape and that only source strength
+  ! influences the amplitude of a plume associated with a given source region.   This routine retrieves the temporal weights
+  ! for the plumes.  Each plume feature has its own temporal weights which varies yearly.  The annual cycle is indexed by
+  ! week in the year and superimposed on the yearly mean value of the weight. 
+  !
+  SUBROUTINE set_time_weight(year_fr)
+    !
+    ! ---------- 
+    !
+    REAL, INTENT(IN) ::  &
+         year_fr           !< Fractional Year (1850.0 - 2100.99)
+
+    INTEGER          ::  &
+         iyear          ,& !< Integer year values between 1 and 156 (1850-2100) 
+         iweek          ,& !< Integer index (between 1 and ntimes); for ntimes=52 this corresponds to weeks (roughly)
+         iplume            ! plume number
+    !
+    ! ---------- 
+    !
+    iyear = FLOOR(year_fr) - 1849
+    iweek = FLOOR((year_fr - FLOOR(year_fr)) * ntimes) + 1
+
+    IF ((iweek > ntimes) .OR. (iweek < 1) .OR. (iyear > nyears) .OR. (iyear < 1)) THEN 
+      CALL abort_physic('set_time_weight','Time out of bounds')
+    ENDIF
+
+    DO iplume=1,nplumes
+      time_weight(1,iplume) = year_weight(iyear,iplume) * ann_cycle(1,iweek,iplume)
+      time_weight(2,iplume) = year_weight(iyear,iplume) * ann_cycle(2,iweek,iplume)
+      time_weight_bg(1,iplume) = ann_cycle(1,iweek,iplume)
+      time_weight_bg(2,iplume) = ann_cycle(2,iweek,iplume) 
+    END DO
+    
+    RETURN
+  END SUBROUTINE set_time_weight
+  !
+  ! ------------------------------------------------------------------------------------------------------------------------
+  ! SP_AOP_PROFILE:  This subroutine calculates the simple plume aerosol and cloud active optical properties based on the
+  ! the simple plume fit to the MPI Aerosol Climatology (Version 2).  It sums over nplumes to provide a profile of aerosol
+  ! optical properties on a host models vertical grid. 
+  !
+  SUBROUTINE sp_aop_profile                                                                           ( &
+       nlevels        ,ncol           ,lambda         ,oro            ,lon            ,lat            , &
+       year_fr        ,z              ,dz             ,dNovrN         ,aod_prof       ,ssa_prof       , &
+       asy_prof       )
+    !
+    ! ---------- 
+    !
+    INTEGER, INTENT(IN)        :: &
+         nlevels,                 & !< number of levels
+         ncol                       !< number of columns
+
+    REAL, INTENT(IN)           :: &
+         lambda,                  & !< wavelength
+         year_fr,                 & !< Fractional Year (1903.0 is the 0Z on the first of January 1903, Gregorian)
+         oro(ncol),               & !< orographic height (m)
+         lon(ncol),               & !< longitude 
+         lat(ncol),               & !< latitude
+         z (ncol,nlevels),        & !< height above sea-level (m)
+         dz(ncol,nlevels)           !< level thickness (difference between half levels) (m)
+
+    REAL, INTENT(OUT)          :: &
+         dNovrN(ncol)           , & !< anthropogenic increase in cloud drop number concentration (factor)
+         aod_prof(ncol,nlevels) , & !< profile of aerosol optical depth
+         ssa_prof(ncol,nlevels) , & !< profile of single scattering albedo
+         asy_prof(ncol,nlevels)     !< profile of asymmetry parameter
+
+    INTEGER                    :: iplume, icol, k
+
+    REAL                       ::  &
+         eta(ncol,nlevels),        & !< normalized height (by 15 km)
+         z_beta(ncol,nlevels),     & !< profile for scaling column optical depth
+         prof(ncol,nlevels),       & !< scaled profile (by beta function)
+         beta_sum(ncol),           & !< vertical sum of beta function
+         ssa(ncol),                & !< single scattering albedo 
+         asy(ncol),                & !< asymmetry parameter
+         cw_an(ncol),              & !< column weight for simple plume (anthropogenic) AOD at 550 nm
+         cw_bg(ncol),              & !< column weight for fine-mode natural background AOD at 550 nm
+         caod_sp(ncol),            & !< column simple plume anthropogenic AOD at 550 nm
+         caod_bg(ncol),            & !< column fine-mode natural background AOD at 550 nm
+         a_plume1,                 & !< gaussian longitude factor for feature 1
+         a_plume2,                 & !< gaussian longitude factor for feature 2
+         b_plume1,                 & !< gaussian latitude factor for feature 1
+         b_plume2,                 & !< gaussian latitude factor for feature 2
+         delta_lat,                & !< latitude offset
+         delta_lon,                & !< longitude offset
+         delta_lon_t,              & !< threshold for maximum longitudinal plume extent used in transition from 360 to 0 degrees
+         lon1,                     & !< rotated longitude for feature 1
+         lat1,                     & !< rotated latitude for feature 2
+         lon2,                     & !< rotated longitude for feature 1
+         lat2,                     & !< rotated latitude for feature 2
+         f1,                       & !< contribution from feature 1
+         f2,                       & !< contribution from feature 2
+         f3,                       & !< contribution from feature 1 in natural background of Twomey effect
+         f4,                       & !< contribution from feature 2 in natural background of Twomey effect
+         aod_550,                  & !< aerosol optical depth at 550nm
+         aod_lmd,                  & !< aerosol optical depth at input wavelength
+         lfactor                     !< factor to compute wavelength dependence of optical properties
+    !
+    ! ---------- 
+    !
+    ! initialize input data (by calling setup at first instance) 
+    !
+    IF (.NOT.sp_initialized) CALL sp_setup
+    !
+    ! get time weights
+    !
+    CALL set_time_weight(year_fr)
+    !
+    ! initialize variables, including output
+    !
+    DO k=1,nlevels
+      DO icol=1,ncol
+        aod_prof(icol,k) = 0.0
+        ssa_prof(icol,k) = 0.0
+        asy_prof(icol,k) = 0.0
+        z_beta(icol,k)   = MERGE(1.0, 0.0, z(icol,k) >= oro(icol))
+        eta(icol,k)      = MAX(0.0,MIN(1.0,z(icol,k)/15000.))
+      END DO
+    END DO
+    DO icol=1,ncol
+      dNovrN(icol)   = 1.0
+      caod_sp(icol)  = 0.0
+      caod_bg(icol)  = 0.02
+    END DO
+    !
+    ! sum contribution from plumes to construct composite profiles of aerosol optical properties
+    !
+    DO iplume=1,nplumes
+      !
+      ! calculate vertical distribution function from parameters of beta distribution
+      !
+      DO icol=1,ncol
+        beta_sum(icol) = 0.
+      END DO
+      DO k=1,nlevels
+        DO icol=1,ncol
+          prof(icol,k)   = (eta(icol,k)**(beta_a(iplume)-1.) * (1.-eta(icol,k))**(beta_b(iplume)-1.)) * dz(icol,k)
+          beta_sum(icol) = beta_sum(icol) + prof(icol,k)
+        END DO
+      END DO
+      DO k=1,nlevels
+        DO icol=1,ncol
+          prof(icol,k)   = ( prof(icol,k) / beta_sum(icol) ) * z_beta(icol,k)
+        END DO
+      END DO
+      !
+      ! calculate plume weights
+      !
+      DO icol=1,ncol
+        !
+        ! get plume-center relative spatial parameters for specifying amplitude of plume at given lat and lon
+        !
+        delta_lat   = lat(icol) - plume_lat(iplume)
+        delta_lon   = lon(icol) - plume_lon(iplume)
+        delta_lon_t = MERGE (260., 180., iplume == 1)
+        delta_lon   = MERGE ( delta_lon-SIGN(360.,delta_lon) , delta_lon , ABS(delta_lon) > delta_lon_t)
+
+        a_plume1  = 0.5 / (MERGE(sig_lon_E(1,iplume), sig_lon_W(1,iplume), delta_lon > 0)**2)
+        b_plume1  = 0.5 / (MERGE(sig_lat_E(1,iplume), sig_lat_W(1,iplume), delta_lon > 0)**2)
+        a_plume2  = 0.5 / (MERGE(sig_lon_E(2,iplume), sig_lon_W(2,iplume), delta_lon > 0)**2)
+        b_plume2  = 0.5 / (MERGE(sig_lat_E(2,iplume), sig_lat_W(2,iplume), delta_lon > 0)**2)
+        !
+        ! adjust for a plume specific rotation which helps match plume state to climatology.
+        !
+        lon1 =   COS(theta(1,iplume))*(delta_lon) + SIN(theta(1,iplume))*(delta_lat)
+        lat1 = - SIN(theta(1,iplume))*(delta_lon) + COS(theta(1,iplume))*(delta_lat)
+        lon2 =   COS(theta(2,iplume))*(delta_lon) + SIN(theta(2,iplume))*(delta_lat)
+        lat2 = - SIN(theta(2,iplume))*(delta_lon) + COS(theta(2,iplume))*(delta_lat)
+        !
+        ! calculate contribution to plume from its different features, to get a column weight for the anthropogenic
+        ! (cw_an) and the fine-mode natural background aerosol (cw_bg)
+        !
+        f1 = time_weight(1,iplume) * ftr_weight(1,iplume) * EXP(-1.* (a_plume1 * ((lon1)**2) + (b_plume1 * ((lat1)**2)))) 
+        f2 = time_weight(2,iplume) * ftr_weight(2,iplume) * EXP(-1.* (a_plume2 * ((lon2)**2) + (b_plume2 * ((lat2)**2)))) 
+        f3 = time_weight_bg(1,iplume) * ftr_weight(1,iplume) * EXP(-1.* (a_plume1 * ((lon1)**2) + (b_plume1 * ((lat1)**2)))) 
+        f4 = time_weight_bg(2,iplume) * ftr_weight(2,iplume) * EXP(-1.* (a_plume2 * ((lon2)**2) + (b_plume2 * ((lat2)**2))))
+
+        cw_an(icol) = f1 * aod_spmx(iplume) + f2 * aod_spmx(iplume)  
+        cw_bg(icol) = f3 * aod_fmbg(iplume) + f4 * aod_fmbg(iplume) 
+        !
+        ! calculate wavelength-dependent scattering properties
+        !
+        lfactor   = MIN(1.0,700.0/lambda)
+        ssa(icol) = (ssa550(iplume) * lfactor**4) / ((ssa550(iplume) * lfactor**4) + ((1-ssa550(iplume)) * lfactor))
+        asy(icol) =  asy550(iplume) * SQRT(lfactor)
+      END DO
+      !
+      ! distribute plume optical properties across its vertical profile weighting by optical depth and scaling for
+      ! wavelength using the angstrom parameter. 
+      !      
+      lfactor = EXP(-angstrom(iplume) * LOG(lambda/550.0))
+      DO k=1,nlevels
+        DO icol = 1,ncol
+          aod_550          = prof(icol,k)     * cw_an(icol)
+          aod_lmd          = aod_550          * lfactor
+          caod_sp(icol)    = caod_sp(icol)    + aod_550
+          caod_bg(icol)    = caod_bg(icol)    + prof(icol,k) * cw_bg(icol)
+          asy_prof(icol,k) = asy_prof(icol,k) + aod_lmd * ssa(icol) * asy(icol)
+          ssa_prof(icol,k) = ssa_prof(icol,k) + aod_lmd * ssa(icol)
+          aod_prof(icol,k) = aod_prof(icol,k) + aod_lmd
+        END DO
+      END DO
+    END DO
+    !
+    ! complete optical depth weighting
+    !
+    DO k=1,nlevels
+      DO icol = 1,ncol
+        asy_prof(icol,k) = MERGE(asy_prof(icol,k)/ssa_prof(icol,k), 0.0, ssa_prof(icol,k) > TINY(1.))
+        ssa_prof(icol,k) = MERGE(ssa_prof(icol,k)/aod_prof(icol,k), 1.0, aod_prof(icol,k) > TINY(1.))
+      END DO
+    END DO
+    !
+    ! calculate effective radius normalization (divisor) factor
+    !
+    DO icol=1,ncol
+      dNovrN(icol) = LOG((1000.0 * (caod_sp(icol) + caod_bg(icol))) + 1.0)/LOG((1000.0 * caod_bg(icol)) + 1.0)
+    END DO
+
+    RETURN
+  END SUBROUTINE sp_aop_profile
+  
+END MODULE MO_SIMPLE_PLUMES
Index: LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/newmicro.F90
===================================================================
--- LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/newmicro.F90	(revision 3314)
+++ LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/newmicro.F90	(revision 3318)
@@ -3,5 +3,5 @@
 
 
-SUBROUTINE newmicro(ok_cdnc, bl95_b0, bl95_b1, paprs, pplay, t, pqlwp, pclc, &
+SUBROUTINE newmicro(flag_aerosol, ok_cdnc, bl95_b0, bl95_b1, paprs, pplay, t, pqlwp, pclc, &
     pcltau, pclemi, pch, pcl, pcm, pct, pctlwp, xflwp, xfiwp, xflwc, xfiwc, &
     mass_solu_aero, mass_solu_aero_pi, pcldtaupi, re, fl, reliq, reice, &
@@ -10,5 +10,6 @@
   USE dimphy
   USE phys_local_var_mod, ONLY: scdnc, cldncl, reffclwtop, lcc, reffclws, &
-    reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra, icc3dcon, icc3dstra, zfice
+    reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra, icc3dcon, icc3dstra, &
+    zfice, dNovrN
   USE phys_state_var_mod, ONLY: rnebcon, clwcon
   USE icefrac_lsc_mod ! computes ice fraction (JBM 3/14)
@@ -142,4 +143,5 @@
   ! within the grid cell)
 
+  INTEGER flag_aerosol
   LOGICAL ok_cdnc
   REAL bl95_b0, bl95_b1 ! Parameter in B&L 95-Formula
@@ -257,8 +259,4 @@
         ! to be within [20, 1000 cm^3]
 
-        ! --present-day case
-        cdnc(i, k) = 10.**(bl95_b0+bl95_b1*log(max(mass_solu_aero(i,k), &
-          1.E-4))/log(10.))*1.E6 !-m-3
-        cdnc(i, k) = min(cdnc_max_m3, max(cdnc_min_m3,cdnc(i,k)))
 
         ! --pre-industrial case
@@ -266,4 +264,18 @@
           1.E-4))/log(10.))*1.E6 !-m-3
         cdnc_pi(i, k) = min(cdnc_max_m3, max(cdnc_min_m3,cdnc_pi(i,k)))
+
+        ! --present-day case
+        ! --flag_aerosol=7 => MACv2SP climatology  
+        ! in this case there is an enhancement factor 
+        IF (flag_aerosol .EQ. 7) THEN
+           cdnc(i, k) = cdnc_pi(i,k)*dNovrN(i)
+        ELSE
+           !--standard case, present day
+           cdnc(i, k) = 10.**(bl95_b0+bl95_b1*log(max(mass_solu_aero(i,k), &
+                1.E-4))/log(10.))*1.E6 !-m-3
+           cdnc(i, k) = min(cdnc_max_m3, max(cdnc_min_m3,cdnc(i,k)))
+        END IF
+
+
 
         ! --present-day case
Index: LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/phys_local_var_mod.F90
===================================================================
--- LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/phys_local_var_mod.F90	(revision 3314)
+++ LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/phys_local_var_mod.F90	(revision 3318)
@@ -151,4 +151,6 @@
       REAL, SAVE, ALLOCATABLE :: scdnc(:,:)
       !$OMP THREADPRIVATE(scdnc)
+      REAL, SAVE, ALLOCATABLE :: dNovrN(:) 
+      !$OMP THREADPRIVATE(dNovrN) 
       REAL, SAVE, ALLOCATABLE :: cldncl(:)
       !$OMP THREADPRIVATE(cldncl)
@@ -615,4 +617,5 @@
       ALLOCATE(tau3d_aero(klon,klev,nwave,naero_tot)) 
       ALLOCATE(scdnc(klon, klev))
+      ALLOCATE(dNovrN(klon))
       ALLOCATE(cldncl(klon))
       ALLOCATE(reffclwtop(klon))
@@ -909,4 +912,5 @@
       DEALLOCATE(tau3d_aero) 
       DEALLOCATE(scdnc)
+      DEALLOCATE(dNovrN)
       DEALLOCATE(cldncl)
       DEALLOCATE(reffclwtop)
Index: LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/physiq_mod.F90
===================================================================
--- LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/physiq_mod.F90	(revision 3314)
+++ LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/physiq_mod.F90	(revision 3318)
@@ -3564,4 +3564,10 @@
                         tausum_aero, drytausum_aero, tau3d_aero)
 #endif
+
+                   IF (flag_aerosol .EQ. 7) THEN
+                      CALL MACv2SP(pphis,pplay,paprs,longitude_deg,latitude_deg,  &
+                                   tau_aero_sw_rrtm,piz_aero_sw_rrtm,cg_aero_sw_rrtm)
+                   ENDIF
+
                    !
                 ELSE IF (NSW.EQ.2) THEN 
@@ -3688,5 +3694,5 @@
 #endif
           ENDIF
-          CALL newmicro (ok_cdnc, bl95_b0, bl95_b1, &
+          CALL newmicro (flag_aerosol, ok_cdnc, bl95_b0, bl95_b1, &
                paprs, pplay, t_seri, cldliq, cldfra, &
                cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
Index: LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/rrtm/aeropt_5wv_rrtm.F90
===================================================================
--- LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/rrtm/aeropt_5wv_rrtm.F90	(revision 3314)
+++ LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/rrtm/aeropt_5wv_rrtm.F90	(revision 3318)
@@ -381,5 +381,5 @@
      ALLOCATE (aerosol_name(nb_aer)) 
      aerosol_name(1) = id_CIDUSTM_phy
-  ELSEIF (flag_aerosol .EQ. 6) THEN 
+  ELSEIF (flag_aerosol .EQ. 6 .OR. flag_aerosol .EQ. 7 ) THEN 
      nb_aer = 13
      ALLOCATE (aerosol_name(nb_aer)) 
Index: LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/rrtm/aeropt_6bands_rrtm.F90
===================================================================
--- LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/rrtm/aeropt_6bands_rrtm.F90	(revision 3314)
+++ LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/rrtm/aeropt_6bands_rrtm.F90	(revision 3318)
@@ -507,5 +507,5 @@
      ALLOCATE (aerosol_name(nb_aer)) 
      aerosol_name(1) = id_CIDUSTM_phy
-  ELSEIF (flag_aerosol .EQ. 6) THEN 
+  ELSEIF (flag_aerosol .EQ. 6 .OR. flag_aerosol .EQ. 7 ) THEN 
      nb_aer = 13
      ALLOCATE (aerosol_name(nb_aer)) 
Index: LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/rrtm/aeropt_lw_rrtm.F90
===================================================================
--- LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/rrtm/aeropt_lw_rrtm.F90	(revision 3314)
+++ LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/rrtm/aeropt_lw_rrtm.F90	(revision 3318)
@@ -46,5 +46,5 @@
   IF (ok_alw) THEN                                   !--aerosol LW effects
    !
-   IF (flag_aerosol.EQ.5.OR.flag_aerosol.EQ.6) THEN  !-Dust
+   IF (flag_aerosol.EQ.5.OR.flag_aerosol.EQ.6.OR.flag_aerosol.EQ.7) THEN  !-Dust
     !
     zdh(:,:)=pdel(:,:)/(RG*zrho(:,:))      ! m
