Index: LMDZ5/trunk/libf/phylmd/phys_local_var_mod.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/phys_local_var_mod.F90	(revision 2534)
+++ LMDZ5/trunk/libf/phylmd/phys_local_var_mod.F90	(revision 2536)
@@ -137,4 +137,6 @@
       !$OMP THREADPRIVATE(ZFSDN0,ZFSUP0)
 !
+      REAL, SAVE, ALLOCATABLE :: stratomask(:,:) 
+      !$OMP THREADPRIVATE(stratomask) 
       REAL, SAVE, ALLOCATABLE :: tausum_aero(:,:,:) 
       !$OMP THREADPRIVATE(tausum_aero) 
@@ -463,4 +465,5 @@
 !      allocate(tausum_aero(klon,nwave,naero_spc))
 !      allocate(tau3d_aero(klon,klev,nwave,naero_spc)) 
+      allocate(stratomask(klon,klev))
 !--correction mini bug OB
       allocate(tausum_aero(klon,nwave,naero_tot))
@@ -670,4 +673,5 @@
       deallocate(topsw0_aero,solsw0_aero)
       deallocate(topswcf_aero,solswcf_aero)
+      deallocate(stratomask) 
       deallocate(tausum_aero) 
       deallocate(tau3d_aero) 
Index: LMDZ5/trunk/libf/phylmd/phys_output_ctrlout_mod.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 2534)
+++ LMDZ5/trunk/libf/phylmd/phys_output_ctrlout_mod.F90	(revision 2536)
@@ -850,7 +850,4 @@
        ctrl_out((/ 2, 6, 10, 10, 10, 10, 11, 11, 11 /),'OD550_'//name_aero_tau(14),       &
        "Aerosol Optical depth at 550 nm "//name_aero_tau(14),"14", (/ ('', i=1, 9) /)) /)
-  
-
-
 !
   TYPE(ctrl_out), SAVE :: o_tausumaero_lw = ctrl_out((/ 2, 6, 10, 10, 10, 10, 11, 11, 11 /), &
@@ -1032,4 +1029,6 @@
   TYPE(ctrl_out), SAVE :: o_lcc3dstra = ctrl_out((/ 2,  6, 10, 10, 10, 10, 11, 11, 11 /), &
     'lcc3dstra', 'Stratiform cloud liquid fraction', '1', (/ ('', i=1, 9) /))
+  TYPE(ctrl_out), SAVE :: o_stratomask = ctrl_out((/ 2,  6, 10, 10, 10, 10, 11, 11, 11 /), &
+    'stratomask', 'Stratospheric fraction', '1', (/ ('', i=1, 9) /))
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
Index: LMDZ5/trunk/libf/phylmd/phys_output_write_mod.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/phys_output_write_mod.F90	(revision 2534)
+++ LMDZ5/trunk/libf/phylmd/phys_output_write_mod.F90	(revision 2536)
@@ -112,5 +112,5 @@
          o_swtoacf_zero, o_swsrfcf_zero, &
          o_topswai, o_solswai, o_scdnc, &
-         o_cldncl, o_reffclws, o_reffclwc, &
+         o_cldncl, o_reffclws, o_reffclwc, o_stratomask,&
          o_cldnvi, o_lcc, o_lcc3d, o_lcc3dcon, &
          o_lcc3dstra, o_reffclwtop, o_ec550aer, &
@@ -230,5 +230,5 @@
          toplwad_aero, toplwad0_aero, sollwad_aero, &
          sollwad0_aero, toplwai_aero, sollwai_aero, &
-         scdnc, cldncl, reffclws, reffclwc, cldnvi, &
+         scdnc, cldncl, reffclws, reffclwc, cldnvi, stratomask,&
          lcc, lcc3d, lcc3dcon, lcc3dstra, reffclwtop, &
          ec550aer, flwc, fiwc, t_seri, theta, q_seri, &
@@ -1037,4 +1037,6 @@
        CALL histwrite_phy(o_zfull,zphi/RG)
 
+       CALL histwrite_phy(o_stratomask, stratomask)
+
        IF (vars_defined)  THEN
         zx_tmp_fi3d(:,1)= pphis(:)/RG
Index: LMDZ5/trunk/libf/phylmd/physiq_mod.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/physiq_mod.F90	(revision 2534)
+++ LMDZ5/trunk/libf/phylmd/physiq_mod.F90	(revision 2536)
@@ -3079,10 +3079,18 @@
              CALL readaerosolstrato(debut)
            ELSE
-             abort_message='flag_aerosol_strat must equal 1 for rrtm=0' 
+             abort_message='flag_aerosol_strat must equal 1 for rrtm=0'
              call abort_physic(modname,abort_message,1)
            ENDIF
           ELSE
 #ifdef CPP_RRTM
-             CALL readaerosolstrato_rrtm(debut)
+            IF (flag_aerosol_strat.EQ.1) THEN
+             CALL readaerosolstrato1_rrtm(debut)
+            ELSEIF (flag_aerosol_strat.EQ.2) THEN
+             CALL stratosphere_mask(t_seri, pplay, latitude_deg)
+             CALL readaerosolstrato2_rrtm(debut)
+            ELSE
+             abort_message='flag_aerosol_strat must equal 1 or 2 for rrtm=1'
+             call abort_physic(modname,abort_message,1)
+            ENDIF
 #else
              abort_message='You should compile with -rrtm if running ' &
@@ -3093,5 +3101,4 @@
        ENDIF
        !--fin STRAT AEROSOL
-
 
        !   On prend la somme des fractions nuageuses et des contenus en eau
Index: LMDZ5/trunk/libf/phylmd/rrtm/readaerosolstrato1_rrtm.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/rrtm/readaerosolstrato1_rrtm.F90	(revision 2536)
+++ LMDZ5/trunk/libf/phylmd/rrtm/readaerosolstrato1_rrtm.F90	(revision 2536)
@@ -0,0 +1,231 @@
+!
+! $Id: readaerosolstrato1_rrtm.F90 2526 2016-05-26 22:13:40Z oboucher $
+!
+subroutine readaerosolstrato1_rrtm(debut)
+
+    use netcdf95, only: nf95_close, nf95_gw_var, nf95_inq_dimid, & 
+                        nf95_inq_varid, nf95_open
+    use netcdf, only: nf90_get_var, nf90_noerr, nf90_nowrite
+
+    USE phys_cal_mod, ONLY : mth_cur
+    USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, klon_glo, grid2dTo1d_glo
+    USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
+    USE mod_phys_lmdz_omp_data, ONLY :  is_omp_root
+    USE mod_phys_lmdz_para 
+    USE phys_state_var_mod
+    USE phys_local_var_mod
+    USE aero_mod
+    USE dimphy
+    USE YOERAD   , ONLY : NLW
+
+    implicit none
+
+    include "YOMCST.h"
+
+! Variable input
+    logical debut
+
+! Variables locales
+    integer n_lat   ! number of latitudes in the input data
+    integer n_lon   ! number of longitudes in the input data
+    integer n_lev   ! number of levels in the input data
+    integer n_month ! number of months in the input data
+    real, pointer:: latitude(:)
+    real, pointer:: longitude(:)
+    real, pointer:: time(:)
+    real, pointer:: lev(:)
+    integer k, band, wave, i
+    integer, save :: mth_pre
+
+    real, allocatable, dimension(:,:), save :: tau_aer_strat
+!$OMP THREADPRIVATE(tau_aer_strat)
+
+! Champs reconstitues
+    real, allocatable:: tauaerstrat(:, :, :, :)
+    real, allocatable:: tauaerstrat_mois(:, :, :)
+    real, allocatable:: tauaerstrat_mois_glo(:, :)
+
+    real, allocatable:: sum_tau_aer_strat(:)
+
+! For NetCDF:
+    integer ncid_in  ! IDs for input files
+    integer varid, ncerr
+
+! Stratospheric aerosols optical properties
+! alpha_sw_strat over the 6 bands is normalised by the 550 nm extinction coefficient
+    real, dimension(nbands_sw_rrtm) :: alpha_sw_strat, piz_sw_strat, cg_sw_strat
+    data alpha_sw_strat/0.8545564, 0.8451642, 0.9821724, 0.8145110, 0.3073565, 7.7966176E-02/
+    data cg_sw_strat   /0.6997170, 0.6810035, 0.7403592, 0.7562674, 0.6676504, 0.3478689/
+    data piz_sw_strat  /0.9999998, 0.9999998, 1.000000000, 0.9999958, 0.9977155, 0.4510679/
+!
+!--diagnostics AOD in the SW
+! alpha_sw_strat_wave is *not* normalised by the 550 nm extinction coefficient
+    real, dimension(nwave) :: alpha_sw_strat_wave
+    data alpha_sw_strat_wave/3.708007,4.125824,4.136584,3.887478,3.507738/
+!
+!--diagnostics AOD in the LW at 10 um 
+    real :: alpha_lw_strat_wave
+    data alpha_lw_strat_wave/0.2746812/
+!
+    real, dimension(nbands_lw_rrtm) :: alpha_lw_abs_rrtm
+    data alpha_lw_abs_rrtm/   8.8340312E-02, 6.9856711E-02, 6.2652975E-02, 5.7188231E-02, &
+                              6.3157059E-02, 5.5072524E-02, 5.0571125E-02, 0.1349073, &    
+                              0.1381676, 9.6506312E-02, 5.1312990E-02, 2.4256418E-02, &
+                              2.7191756E-02, 3.3862915E-02, 1.6132960E-02, 1.4275438E-02/ ! calculated with Mie_SW_LW_RRTM_V2.4 (bimodal, corrected)
+                                                                                          ! for r_0=/0.13E-6, 0.41E-6/ m, sigma_g=/1.26, 1.30/
+                                                                                          ! order: increasing wavelength!
+!--------------------------------------------------------
+
+    IF (.not.ALLOCATED(tau_aer_strat)) ALLOCATE(tau_aer_strat(klon,klev))
+    IF (.not.ALLOCATED(sum_tau_aer_strat)) ALLOCATE(sum_tau_aer_strat(klon))
+
+!--we only read monthly strat aerosol data
+    IF (debut.OR.mth_cur.NE.mth_pre) THEN
+
+!--only root reads the data
+    IF (is_mpi_root.AND.is_omp_root) THEN
+
+    IF (nbands_sw_rrtm.NE.6) THEN 
+        print *,'nbands_sw_rrtm doit etre egal a 6 dans readaerosolstrat_rrtm'
+        STOP
+    ENDIF
+
+    CALL nf95_open("taustrat.nc", nf90_nowrite, ncid_in)
+
+    CALL nf95_inq_varid(ncid_in, "LEV", varid)
+    CALL nf95_gw_var(ncid_in, varid, lev)
+    n_lev = size(lev)
+    IF (n_lev.NE.klev) THEN 
+       print *,'Le nombre de niveaux n est pas egal a klev'
+       STOP
+    ENDIF
+
+    CALL nf95_inq_varid(ncid_in, "LAT", varid)
+    CALL nf95_gw_var(ncid_in, varid, latitude)
+    n_lat = size(latitude)
+    print *, 'LAT aerosol strato=', n_lat, latitude
+    IF (n_lat.NE.nbp_lat) THEN 
+       print *,'Le nombre de lat n est pas egal a nbp_lat'
+       STOP
+    ENDIF
+
+    CALL nf95_inq_varid(ncid_in, "LON", varid)
+    CALL nf95_gw_var(ncid_in, varid, longitude)
+    n_lon = size(longitude)
+    print *, 'LON aerosol strato=', n_lon, longitude
+    IF (n_lon.NE.nbp_lon) THEN 
+       print *,'Le nombre de lon n est pas egal a nbp_lon'
+       STOP
+    ENDIF
+
+    CALL nf95_inq_varid(ncid_in, "TIME", varid)
+    CALL nf95_gw_var(ncid_in, varid, time)
+    n_month = size(time)
+    print *, 'TIME aerosol strato=', n_month, time
+    IF (n_month.NE.12) THEN 
+       print *,'Le nombre de month n est pas egal a 12'
+       STOP
+    ENDIF
+
+    ALLOCATE(tauaerstrat(n_lon, n_lat, n_lev, n_month))
+    ALLOCATE(tauaerstrat_mois(n_lon, n_lat, n_lev))
+    ALLOCATE(tauaerstrat_mois_glo(klon_glo, n_lev))
+
+!--reading stratospheric AOD at 550 nm
+    CALL nf95_inq_varid(ncid_in, "TAUSTRAT", varid)
+    ncerr = nf90_get_var(ncid_in, varid, tauaerstrat)
+    print *,'code erreur readaerosolstrato=', ncerr, varid
+
+    CALL nf95_close(ncid_in)
+
+!---select the correct month
+    IF (mth_cur.LT.1.OR.mth_cur.GT.12) THEN
+      print *,'probleme avec le mois dans readaerosolstrat =', mth_cur
+    ENDIF
+    tauaerstrat_mois(:,:,:) = tauaerstrat(:,:,:,mth_cur)
+
+!---reduce to a klon_glo grid 
+    CALL grid2dTo1d_glo(tauaerstrat_mois,tauaerstrat_mois_glo)
+
+    ENDIF !--is_mpi_root and is_omp_root
+
+!$OMP BARRIER
+
+!--keep memory of previous month
+    mth_pre=mth_cur
+
+!--scatter on all proc
+    CALL scatter(tauaerstrat_mois_glo,tau_aer_strat)
+
+    IF (is_mpi_root.AND.is_omp_root) THEN
+!
+    DEALLOCATE(tauaerstrat)
+    DEALLOCATE(tauaerstrat_mois)
+    DEALLOCATE(tauaerstrat_mois_glo)
+!
+    ENDIF !--is_mpi_root and is_omp_root
+
+!$OMP BARRIER
+
+    ENDIF !--debut ou nouveau mois
+
+!--total vertical aod at the 5 SW wavelengths
+    DO wave=1, nwave
+    DO k=1, klev
+    tausum_aero(:,wave,id_STRAT_phy)=tausum_aero(:,wave,id_STRAT_phy)+ &
+       tau_aer_strat(:,k)*alpha_sw_strat_wave(wave)/alpha_sw_strat_wave(2)
+    ENDDO
+    ENDDO
+
+!--weighted average for cg, piz and tau, adding strat aerosols on top of tropospheric ones
+    DO band=1, nbands_sw_rrtm
+!--anthropogenic aerosols bands 1 to nbands_sw_rrtm 
+    cg_aero_sw_rrtm(:,:,2,band)  = ( cg_aero_sw_rrtm(:,:,2,band)*piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) + &
+                                  cg_sw_strat(band)*piz_sw_strat(band)*alpha_sw_strat(band)*tau_aer_strat(:,:) ) /           &
+                             MAX( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) +                                &
+                                  piz_sw_strat(band)*alpha_sw_strat(band)*tau_aer_strat(:,:), 1.e-15 )
+    piz_aero_sw_rrtm(:,:,2,band)  = ( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) +                            &
+                              piz_sw_strat(band)*alpha_sw_strat(band)*tau_aer_strat(:,:) ) /                                 &
+                              MAX( tau_aero_sw_rrtm(:,:,2,band) + alpha_sw_strat(band)*tau_aer_strat(:,:), 1.e-15 )
+    tau_aero_sw_rrtm(:,:,2,band)  = tau_aero_sw_rrtm(:,:,2,band) + alpha_sw_strat(band)*tau_aer_strat(:,:)
+!--natural aerosols bands 1 to nbands_sw_rrtm
+    cg_aero_sw_rrtm(:,:,1,band)  = ( cg_aero_sw_rrtm(:,:,1,band)*piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) + &
+                             cg_sw_strat(band)*piz_sw_strat(band)*alpha_sw_strat(band)*tau_aer_strat(:,:) ) /                &
+                             MAX( piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) +                                &
+                                  piz_sw_strat(band)*alpha_sw_strat(band)*tau_aer_strat(:,:), 1.e-15 )
+    piz_aero_sw_rrtm(:,:,1,band)  = ( piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) +                            &
+                              piz_sw_strat(band)*alpha_sw_strat(band)*tau_aer_strat(:,:) ) /                                 &
+                              MAX( tau_aero_sw_rrtm(:,:,1,band) + alpha_sw_strat(band)*tau_aer_strat(:,:),1.e-15 )
+    tau_aero_sw_rrtm(:,:,1,band)  = tau_aero_sw_rrtm(:,:,1,band) + alpha_sw_strat(band)*tau_aer_strat(:,:)
+!--no stratospheric aerosol in index 1 for these tests
+!    cg_aero_sw_rrtm(:,:,1,band)  =  cg_aero_sw_rrtm(:,:,1,band)
+!    piz_aero_sw_rrtm(:,:,1,band)  = piz_aero_sw_rrtm(:,:,1,band)
+!    tau_aero_sw_rrtm(:,:,1,band)  = tau_aero_sw_rrtm(:,:,1,band)
+    ENDDO
+
+!--stratospheric AOD in LW
+    IF (nbands_lw_rrtm .NE. NLW) then
+      print*, 'different values for NLW (=',NLW,') and nbands_lw_rrtm (=', nbands_lw_rrtm, ')'
+      STOP
+    ENDIF 
+
+    DO band=1, nbands_lw_rrtm
+    tau_aero_lw_rrtm(:,:,2,band)  = tau_aero_lw_rrtm(:,:,2,band) + alpha_lw_abs_rrtm(band)*tau_aer_strat(:,:)
+    tau_aero_lw_rrtm(:,:,1,band)  = tau_aero_lw_rrtm(:,:,1,band) + alpha_lw_abs_rrtm(band)*tau_aer_strat(:,:) 
+!--no stratospheric aerosols in index 1 for these tests
+!    tau_aero_lw_rrtm(:,:,1,band)  = tau_aero_lw_rrtm(:,:,1,band) 
+    ENDDO
+
+!--default SSA value if there is no aerosol
+!--to avoid 0 values that seems to cause some problem to RRTM
+    WHERE (tau_aero_sw_rrtm.LT.1.e-14)
+      piz_aero_sw_rrtm = 1.0
+    ENDWHERE
+
+!--in principle this should not be necessary 
+!--as these variables have min values already but just in case
+!--put 1e-15 min value to both SW and LW AOD
+    tau_aero_sw_rrtm = MAX(tau_aero_sw_rrtm,1.e-15)
+    tau_aero_lw_rrtm = MAX(tau_aero_lw_rrtm,1.e-15)
+
+end subroutine readaerosolstrato1_rrtm
Index: LMDZ5/trunk/libf/phylmd/rrtm/readaerosolstrato2_rrtm.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/rrtm/readaerosolstrato2_rrtm.F90	(revision 2536)
+++ LMDZ5/trunk/libf/phylmd/rrtm/readaerosolstrato2_rrtm.F90	(revision 2536)
@@ -0,0 +1,339 @@
+!
+! $Id: readaerosolstrato2_rrtm.F90 2526 2016-05-26 22:13:40Z oboucher $
+!
+subroutine readaerosolstrato2_rrtm(debut)
+
+    use netcdf95, only: nf95_close, nf95_gw_var, nf95_inq_dimid, & 
+                        nf95_inq_varid, nf95_open
+    use netcdf, only: nf90_get_var, nf90_noerr, nf90_nowrite
+
+    USE phys_cal_mod, ONLY : mth_cur
+    USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, klon_glo, grid2dTo1d_glo
+    USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
+    USE mod_phys_lmdz_omp_data, ONLY :  is_omp_root
+    USE mod_phys_lmdz_para 
+    USE phys_state_var_mod
+    USE phys_local_var_mod
+    USE aero_mod
+    USE dimphy
+    USE YOERAD   , ONLY : NLW
+
+    implicit none
+
+    include "YOMCST.h"
+
+    CHARACTER (len = 80) :: abort_message
+    CHARACTER (LEN=20) :: modname = 'readaerosolstrato2'
+
+! Variable input
+    logical, intent(in) ::  debut
+
+! Variables locales
+    integer n_lat   ! number of latitudes in the input data
+    integer n_lon   ! number of longitudes
+    integer n_lev   ! number of levels in the input data
+    integer n_month ! number of months in the input data
+    integer n_wav   ! number of wavelengths in the input data
+    real, pointer:: latitude(:)
+    real, pointer:: longitude(:)
+    real, pointer:: time(:)
+    real, pointer:: lev(:)
+    real, pointer:: wav(:)
+    integer i,k,wave,band
+    integer, save :: mth_pre
+
+    real, allocatable, dimension(:,:,:), save :: tau_aer_strat
+    real, allocatable, dimension(:,:,:), save :: piz_aer_strat
+    real, allocatable, dimension(:,:,:), save :: cg_aer_strat
+    real, allocatable, dimension(:,:,:), save :: taulw_aer_strat
+!$OMP THREADPRIVATE(tau_aer_strat,piz_aer_strat,cg_aer_strat,taulw_aer_strat)
+
+! Champs reconstitues
+    real, allocatable:: tauaerstrat(:, :, :, :)
+    real, allocatable:: pizaerstrat(:, :, :, :)
+    real, allocatable:: cgaerstrat(:, :, :, :)
+    real, allocatable:: taulwaerstrat(:, :, :, :)
+
+    real, allocatable:: tauaerstrat_mois(:, :, :, :)
+    real, allocatable:: pizaerstrat_mois(:, :, :, :)
+    real, allocatable:: cgaerstrat_mois(:, :, :, :)
+    real, allocatable:: taulwaerstrat_mois(:, :, :, :)
+
+    real, allocatable:: tauaerstrat_mois_glo(:, :, :)
+    real, allocatable:: pizaerstrat_mois_glo(:, :, :)
+    real, allocatable:: cgaerstrat_mois_glo(:, :, :)
+    real, allocatable:: taulwaerstrat_mois_glo(:, :, :)
+
+! For NetCDF:
+    integer ncid_in  ! IDs for input files
+    integer varid, ncerr
+
+!--------------------------------------------------------
+
+    IF (.not.ALLOCATED(tau_aer_strat)) ALLOCATE(tau_aer_strat(klon,klev,nbands_sw_rrtm))
+    IF (.not.ALLOCATED(piz_aer_strat)) ALLOCATE(piz_aer_strat(klon,klev,nbands_sw_rrtm))
+    IF (.not.ALLOCATED(cg_aer_strat))  ALLOCATE(cg_aer_strat(klon,klev,nbands_sw_rrtm))
+
+    IF (.not.ALLOCATED(taulw_aer_strat)) ALLOCATE(taulw_aer_strat(klon,klev,nbands_lw_rrtm))
+
+!--we only read monthly strat aerosol data
+    IF (debut.OR.mth_cur.NE.mth_pre) THEN
+
+!--only root reads the data
+    IF (is_mpi_root.AND.is_omp_root) THEN
+
+!--check mth_cur
+    IF (mth_cur.LT.1.OR.mth_cur.GT.12) THEN
+      print *,'probleme avec le mois dans readaerosolstrat =', mth_cur
+    ENDIF
+
+!--initialize n_lon as input data is 2D (lat-alt) only
+    n_lon = nbp_lon
+
+!--Starts with SW optical properties
+    IF (nbands_sw_rrtm.NE.6) THEN 
+       abort_message='nbands_sw_rrtm doit etre egal a 6 dans readaerosolstrat_rrtm'
+       CALL abort_physic(modname,abort_message,1)
+    ENDIF
+
+    CALL nf95_open("tauswstrat.2D.nc", nf90_nowrite, ncid_in)
+
+    CALL nf95_inq_varid(ncid_in, "LEV", varid)
+    CALL nf95_gw_var(ncid_in, varid, lev)
+    n_lev = size(lev)
+    IF (n_lev.NE.klev) THEN 
+       abort_message='Le nombre de niveaux n est pas egal a klev'
+       CALL abort_physic(modname,abort_message,1)
+    ENDIF
+
+    CALL nf95_inq_varid(ncid_in, "LAT", varid)
+    CALL nf95_gw_var(ncid_in, varid, latitude)
+    n_lat = size(latitude)
+    IF (n_lat.NE.nbp_lat) THEN 
+       print *, 'latitude=', n_lat, nbp_lat
+       abort_message='Le nombre de lat n est pas egal a nbp_lat'
+       CALL abort_physic(modname,abort_message,1)
+    ENDIF
+
+    CALL nf95_inq_varid(ncid_in, "TIME", varid)
+    CALL nf95_gw_var(ncid_in, varid, time)
+    n_month = size(time)
+    IF (n_month.NE.12) THEN 
+       abort_message='Le nombre de month n est pas egal a 12'
+       CALL abort_physic(modname,abort_message,1)
+    ENDIF
+
+    CALL nf95_inq_varid(ncid_in, "WAV", varid)
+    CALL nf95_gw_var(ncid_in, varid, wav)
+    n_wav = size(wav)
+    print *, 'WAV aerosol strato=', n_wav, wav
+    IF (n_wav.NE.nbands_sw_rrtm) THEN 
+       abort_message='Le nombre de wav n est pas egal a NSW'
+       CALL abort_physic(modname,abort_message,1)
+    ENDIF
+
+    ALLOCATE(tauaerstrat(n_lat, n_lev, n_wav, n_month))
+    ALLOCATE(pizaerstrat(n_lat, n_lev, n_wav, n_month))
+    ALLOCATE(cgaerstrat(n_lat, n_lev, n_wav, n_month))
+
+    ALLOCATE(tauaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
+    ALLOCATE(pizaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
+    ALLOCATE(cgaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
+
+    ALLOCATE(tauaerstrat_mois_glo(klon_glo, n_lev, n_wav))
+    ALLOCATE(pizaerstrat_mois_glo(klon_glo, n_lev, n_wav))
+    ALLOCATE(cgaerstrat_mois_glo(klon_glo, n_lev, n_wav))
+
+!--reading stratospheric aerosol tau per layer
+    CALL nf95_inq_varid(ncid_in, "TAU_SUN", varid)
+    ncerr = nf90_get_var(ncid_in, varid, tauaerstrat)
+    print *,'code erreur readaerosolstrato=', ncerr, varid
+
+!--reading stratospheric aerosol omega per layer
+    CALL nf95_inq_varid(ncid_in, "OME_SUN", varid)
+    ncerr = nf90_get_var(ncid_in, varid, pizaerstrat)
+    print *,'code erreur readaerosolstrato=', ncerr, varid
+
+!--reading stratospheric aerosol g per layer
+    CALL nf95_inq_varid(ncid_in, "GGG_SUN", varid)
+    ncerr = nf90_get_var(ncid_in, varid, cgaerstrat)
+    print *,'code erreur readaerosolstrato sw=', ncerr, varid
+
+    CALL nf95_close(ncid_in)
+
+!--select the correct month
+!--and copy into 1st longitude
+    tauaerstrat_mois(1,:,:,:) = tauaerstrat(:,:,:,mth_cur)
+    pizaerstrat_mois(1,:,:,:) = pizaerstrat(:,:,:,mth_cur)
+    cgaerstrat_mois(1,:,:,:)  = cgaerstrat(:,:,:,mth_cur)
+
+!--copy longitudes
+    DO i=2, n_lon
+     tauaerstrat_mois(i,:,:,:) = tauaerstrat_mois(1,:,:,:)
+     pizaerstrat_mois(i,:,:,:) = pizaerstrat_mois(1,:,:,:)
+      cgaerstrat_mois(i,:,:,:)  = cgaerstrat_mois(1,:,:,:)
+    ENDDO
+
+!---reduce to a klon_glo grid 
+    DO band=1, nbands_sw_rrtm
+      CALL grid2dTo1d_glo(tauaerstrat_mois(:,:,:,band),tauaerstrat_mois_glo(:,:,band))
+      CALL grid2dTo1d_glo(pizaerstrat_mois(:,:,:,band),pizaerstrat_mois_glo(:,:,band))
+      CALL grid2dTo1d_glo(cgaerstrat_mois(:,:,:,band),cgaerstrat_mois_glo(:,:,band))
+    ENDDO
+
+!--Now LW optical properties
+!
+    IF (nbands_lw_rrtm .NE. NLW) then
+      abort_message='different values for NLW and nbands_lw_rrtm'
+      CALL abort_physic(modname,abort_message,1)
+    ENDIF 
+
+    CALL nf95_open("taulwstrat.2D.nc", nf90_nowrite, ncid_in)
+
+    CALL nf95_inq_varid(ncid_in, "LEV", varid)
+    CALL nf95_gw_var(ncid_in, varid, lev)
+    n_lev = size(lev)
+    IF (n_lev.NE.klev) THEN 
+       abort_message='Le nombre de niveaux n est pas egal a klev'
+       CALL abort_physic(modname,abort_message,1)
+    ENDIF
+
+    CALL nf95_inq_varid(ncid_in, "LAT", varid)
+    CALL nf95_gw_var(ncid_in, varid, latitude)
+    n_lat = size(latitude)
+    IF (n_lat.NE.nbp_lat) THEN 
+       abort_message='Le nombre de lat n est pas egal a nbp_lat'
+       CALL abort_physic(modname,abort_message,1)
+    ENDIF
+
+    CALL nf95_inq_varid(ncid_in, "TIME", varid)
+    CALL nf95_gw_var(ncid_in, varid, time)
+    n_month = size(time)
+    IF (n_month.NE.12) THEN 
+       abort_message='Le nombre de month n est pas egal a 12'
+       CALL abort_physic(modname,abort_message,1)
+    ENDIF
+
+    CALL nf95_inq_varid(ncid_in, "WAV", varid)
+    CALL nf95_gw_var(ncid_in, varid, wav)
+    n_wav = size(wav)
+    print *, 'WAV aerosol strato=', n_wav, wav
+    IF (n_wav.NE.nbands_lw_rrtm) THEN 
+       abort_message='Le nombre de wav n est pas egal a NLW'
+       CALL abort_physic(modname,abort_message,1)
+    ENDIF
+
+    ALLOCATE(taulwaerstrat(n_lat, n_lev, n_wav, n_month))
+    ALLOCATE(taulwaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
+    ALLOCATE(taulwaerstrat_mois_glo(klon_glo, n_lev, n_wav))
+
+!--reading stratospheric aerosol lw tau per layer
+    CALL nf95_inq_varid(ncid_in, "TAU_EAR", varid)
+    ncerr = nf90_get_var(ncid_in, varid, taulwaerstrat)
+    print *,'code erreur readaerosolstrato lw=', ncerr, varid
+
+    CALL nf95_close(ncid_in)
+
+!--select the correct month
+!--and copy into 1st longitude
+    taulwaerstrat_mois(1,:,:,:) = taulwaerstrat(:,:,:,mth_cur)
+!--copy longitudes
+    DO i=2, n_lon
+      taulwaerstrat_mois(i,:,:,:) = taulwaerstrat_mois(1,:,:,:)
+    ENDDO
+
+!---reduce to a klon_glo grid 
+    DO band=1, nbands_lw_rrtm
+      CALL grid2dTo1d_glo(taulwaerstrat_mois(:,:,:,band),taulwaerstrat_mois_glo(:,:,band))
+    ENDDO
+
+    ENDIF !--is_mpi_root and is_omp_root
+
+!$OMP BARRIER
+
+!--keep memory of previous month
+    mth_pre=mth_cur
+
+!--scatter on all proc
+    CALL scatter(tauaerstrat_mois_glo,tau_aer_strat)
+    CALL scatter(pizaerstrat_mois_glo,piz_aer_strat)
+    CALL scatter(cgaerstrat_mois_glo,cg_aer_strat)
+    CALL scatter(taulwaerstrat_mois_glo,taulw_aer_strat)
+
+    IF (is_mpi_root.AND.is_omp_root) THEN
+!
+    DEALLOCATE(tauaerstrat, pizaerstrat, cgaerstrat)
+    DEALLOCATE(tauaerstrat_mois, pizaerstrat_mois, cgaerstrat_mois)
+    DEALLOCATE(tauaerstrat_mois_glo,pizaerstrat_mois_glo,cgaerstrat_mois_glo)
+
+    DEALLOCATE(taulwaerstrat,taulwaerstrat_mois,taulwaerstrat_mois_glo)
+!
+    ENDIF !--is_mpi_root and is_omp_root
+
+!$OMP BARRIER
+
+    ENDIF !--debut ou nouveau mois
+
+!--total vertical aod at the 5 SW wavelengths
+!--for now use band 3 AOD into all 5 wavelengths
+    band=3
+    DO i=1, klon
+    DO k=1, klev
+      IF (stratomask(i,k).GT.0.999999) THEN
+        DO wave=1, nwave
+          tausum_aero(:,wave,id_STRAT_phy)=tausum_aero(:,wave,id_STRAT_phy)+tau_aer_strat(:,k,band)
+        ENDDO
+      ENDIF
+    ENDDO
+    ENDDO
+
+!--weighted average for cg, piz and tau, adding strat aerosols on top of tropospheric ones
+    DO band=1, nbands_sw_rrtm
+    WHERE (stratomask.GT.0.999999)
+!--anthropogenic aerosols bands 1 to nbands_sw_rrtm 
+    cg_aero_sw_rrtm(:,:,2,band)  = ( cg_aero_sw_rrtm(:,:,2,band)*piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) + &
+                                     cg_aer_strat(:,:,band)*piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) /              &
+                                MAX( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) +                             &
+                                     piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band), 1.e-15 )
+    piz_aero_sw_rrtm(:,:,2,band) = ( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) +                             &
+                                     piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) /                                     &
+                                MAX( tau_aero_sw_rrtm(:,:,2,band) + tau_aer_strat(:,:,band), 1.e-15 )
+    tau_aero_sw_rrtm(:,:,2,band)  = tau_aero_sw_rrtm(:,:,2,band) + tau_aer_strat(:,:,band)
+!--natural aerosols bands 1 to nbands_sw_rrtm
+    cg_aero_sw_rrtm(:,:,1,band)  = ( cg_aero_sw_rrtm(:,:,1,band)*piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) + &
+                                     cg_aer_strat(:,:,band)*piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) /              &
+                                MAX( piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) +                             &
+                                     piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band), 1.e-15 )
+    piz_aero_sw_rrtm(:,:,1,band) = ( piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) +                             &
+                                     piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) /                                     &
+                                MAX( tau_aero_sw_rrtm(:,:,1,band) + tau_aer_strat(:,:,band), 1.e-15 )
+    tau_aero_sw_rrtm(:,:,1,band)  = tau_aero_sw_rrtm(:,:,1,band) + tau_aer_strat(:,:,band)
+!--no stratospheric aerosol in index 1 for these tests
+!    cg_aero_sw_rrtm(:,:,1,band)  =  cg_aero_sw_rrtm(:,:,1,band)
+!    piz_aero_sw_rrtm(:,:,1,band)  = piz_aero_sw_rrtm(:,:,1,band)
+!    tau_aero_sw_rrtm(:,:,1,band)  = tau_aero_sw_rrtm(:,:,1,band)
+    ENDWHERE
+    ENDDO
+
+    DO band=1, nbands_lw_rrtm
+    WHERE (stratomask.GT.0.999999)
+    tau_aero_lw_rrtm(:,:,2,band)  = tau_aero_lw_rrtm(:,:,2,band) + taulw_aer_strat(:,:,band)
+    tau_aero_lw_rrtm(:,:,1,band)  = tau_aero_lw_rrtm(:,:,1,band) + taulw_aer_strat(:,:,band)
+!--no stratospheric aerosols in index 1 for these tests
+!    tau_aero_lw_rrtm(:,:,1,band)  = tau_aero_lw_rrtm(:,:,1,band) 
+    ENDWHERE
+    ENDDO
+
+!--default SSA value if there is no aerosol
+!--to avoid 0 values that seems to cause some problem to RRTM
+    WHERE (tau_aero_sw_rrtm.LT.1.e-14)
+      piz_aero_sw_rrtm = 1.0
+    ENDWHERE
+
+!--in principle this should not be necessary 
+!--as these variables have min values already but just in case
+!--put 1e-15 min value to both SW and LW AOD
+    tau_aero_sw_rrtm = MAX(tau_aero_sw_rrtm,1.e-15)
+    tau_aero_lw_rrtm = MAX(tau_aero_lw_rrtm,1.e-15)
+
+end subroutine readaerosolstrato2_rrtm
Index: LMDZ5/trunk/libf/phylmd/stratosphere_mask.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/stratosphere_mask.F90	(revision 2534)
+++ LMDZ5/trunk/libf/phylmd/stratosphere_mask.F90	(revision 2536)
@@ -1,3 +1,3 @@
-SUBROUTINE stratosphere_mask(t_seri, pplay, xlat, is_strato)
+SUBROUTINE stratosphere_mask(t_seri, pplay, xlat)
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -27,5 +27,5 @@
 
 USE dimphy
-!!USE phys_local_var_mod, ONLY: p_tropopause
+USE phys_local_var_mod, ONLY: stratomask
 
 IMPLICIT NONE
@@ -34,6 +34,4 @@
 REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay   ! pression pour le mileu de chaque couche (en Pa)
 REAL,DIMENSION(klon),INTENT(IN)        :: xlat    ! latitudes pour chaque point 
-
-LOGICAL,DIMENSION(klon,klev),INTENT(OUT) :: is_strato
 
 REAL, PARAMETER                        :: plimu=45000.
@@ -84,7 +82,7 @@
 DO k=1, klev
   IF (pplay(i,k).LT.tp(i)) THEN
-    is_strato(i,k)=.TRUE.
+    stratomask(i,k)=1.0
   ELSE
-    is_strato(i,k)=.FALSE.
+    stratomask(i,k)=0.0
   ENDIF
 ENDDO
