Index: trunk/LMDZ.GENERIC/README
===================================================================
--- trunk/LMDZ.GENERIC/README	(revision 2283)
+++ trunk/LMDZ.GENERIC/README	(revision 2297)
@@ -1512,5 +1512,5 @@
   Everybody should be using AU, that's all.
 
-== 24/02/2020 (r2244) == JVO
+== 24/02/2020 (r2245) == JVO
 + Add a 'versH2H2cia' int key in callphys that allows two values (2011 or 2018) to
   deal with updated HITRAN file (for interpolateH2H2.F90) from 2018 that includes the
@@ -1533,3 +1533,29 @@
   + Add a sanity check in callcorrk instead of leaving out-of-bounds planckir indexes.
 
-
+== 26/04/2020 (r2297) == JVO
+Add a generic n-layer aerosol scheme to replace the former buggy 2-layer scheme as well as the hard-coded NH3 cloud.
+ It can be called using 'aeronlay=.true.' in callphys.def, and set the number of layers (up to 4) with 'nlayaero'.
+ Then, the following parameters are read as arrays of size nlayaero in callphys.def (separated by blank space)
+ 
+*aeronlay_tauref  (Optical depth of aerosol layer at ref wavelenght)
+*aeronlay_lamref  (Ref wavelenght (m))
+*aeronlay_choice  (Choice for vertical profile - 1:tau follows atm scale height btwn top and bottom - 2:tau follows it own scale height)
+*aeronlay_pbot    (Bottom pressure (Pa))
+*aeronlay_ptop    (Top pressure (Pa) - useful only if choice=1)
+*aeronlay_sclhght (Ratio of aerosol layer scale height / atmospheric scale height  - useful only if choice=2 )
+*aeronlay_size    (Particle size (m))
+*optprop_aeronlay_vis  File for VIS opt properties.
+*optprop_aeronlay_ir   File for IR opt properties.
+
++Extra info :
+   + In addition of solving the bug from 2-layer it enables different optical properties.
+   + The former scheme are left for retrocompatibility (for now) but you should use the new one.
+   + See aeropacity.F90 for the calculations
+        + Each layer can have different optical properties, size of particle ...
+        + You have different choices for vertical profile of the aerosol layers :
+           * aeronlay_choice = 1 : Layer tau is spread between ptop and pbot following atm scale height.
+           * aeronlay_choice = 2 : Layer tau follows its own scale height above cloud deck (pbot).
+                                   In this case ptop is dummy and sclhght gives the ratio H_cl/H_atm.
+
+   + The reference wavelenght for input optical depth is now read as input (aeronlay_lamref)
+   + Following the last point some comment is added in suaer_corrk about the 'not-really-dummy'ness of IR lamref..
Index: trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90	(revision 2283)
+++ trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90	(revision 2297)
@@ -11,5 +11,7 @@
 		pres_bottom_tropo,pres_top_tropo,obs_tau_col_tropo,	&
 		pres_bottom_strato,pres_top_strato,obs_tau_col_strato,  &
-                tau_nh3_cloud, pres_nh3_cloud 
+                tau_nh3_cloud, pres_nh3_cloud,                          &
+                nlayaero, aeronlay_tauref, aeronlay_choice,             & 
+                aeronlay_pbot, aeronlay_ptop, aeronlay_sclhght
                   
        implicit none
@@ -27,4 +29,5 @@
 !     update J.-B. Madeleine (2008)
 !     dust removal, simplification by Robin Wordsworth (2009)
+!     Generic n-layer aerosol - J. Vatant d'Ollone (2020)
 !
 !     Input
@@ -62,10 +65,12 @@
       logical,intent(in) :: clearsky
 
-      real aerosol0, obs_tau_col_aurora, pm, pente_cloud
-      
+      real aerosol0, obs_tau_col_aurora, pm
+      real pcloud_deck, cloud_slope
+
       real dp_strato(ngrid)
       real dp_tropo(ngrid)
-
-      INTEGER l,ig,iq,iaer
+      real dp_layer(ngrid)
+
+      INTEGER l,ig,iq,iaer,ia
 
       LOGICAL,SAVE :: firstcall=.true.
@@ -129,4 +134,7 @@
         if (iaero_nh3.ne.0) then
           print*,'iaero_nh3= ',iaero_nh3
+        endif
+        if (iaero_nlay(1).ne.0) then
+          print*,'iaero_nlay= ',iaero_nlay(:)
         endif
         if (iaero_aurora.ne.0) then
@@ -375,4 +383,8 @@
 !    Two-layer aerosols (unknown composition)
 !    S. Guerlet (2013) - Modif by J. Vatant d'Ollone (2020)
+!    
+!    This scheme is deprecated and left for retrocompatibility
+!    You should use the n-layer scheme below !
+!
 !==================================================================
 
@@ -435,5 +447,4 @@
                  aerosol(ig,l,iaer) = obs_tau_col_strato*aerosol(ig,l,iaer)/dp_strato(ig)
                ENDIF
-               write(*,*), aerosol(ig,l,iaer)
             ENDDO
          ENDDO
@@ -445,4 +456,5 @@
 !    Saturn/Jupiter ammonia cloud = thin cloud (scale height 0.2 hard coded...)
 !    S. Guerlet (2013)
+!    JVO 20 : You should now use the generic n-layer scheme below
 !==================================================================
 
@@ -461,6 +473,6 @@
 
              ELSEIF (pplev(ig,l) .le. pres_nh3_cloud ) THEN 
-	     pente_cloud=5. !!(hard-coded, correspond to scale height 0.2)
-             aerosol(ig,l,iaer) = ((pplev(ig,l)/pres_nh3_cloud)**(pente_cloud))*tau_nh3_cloud 
+	     cloud_slope=5. !!(hard-coded, correspond to scale height 0.2)
+             aerosol(ig,l,iaer) = ((pplev(ig,l)/pres_nh3_cloud)**(cloud_slope))*tau_nh3_cloud 
 
              ENDIF
@@ -470,8 +482,8 @@
 	  
 !       3. Re-normalize to observed total column
-         tau_col(:)=0.0
+         dp_layer(:)=0.0
          DO l=1,nlayer
           DO ig=1,ngrid
-               tau_col(ig) = tau_col(ig) &
+               dp_layer(ig) = dp_layer(ig) &
                      + aerosol(ig,l,iaer)/tau_nh3_cloud
             ENDDO
@@ -480,9 +492,102 @@
          DO ig=1,ngrid
            DO l=1,nlayer-1
-                aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/tau_col(ig)
+                aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/dp_layer(ig)
            ENDDO
          ENDDO
 
      end if ! if NH3 cloud  
+
+!=========================================================================================================
+!    Generic N-layers aerosols/clouds
+!    Author : J. Vatant d'Ollone (2020)
+!    
+!    Purpose: Replaces and extents the former buggy 2-layer scheme as well as hard-coded NH3 cloud
+!    
+!    + Each layer can have different optical properties, size of particle ...
+!    + Enables up to n=4 layers as we apparently cannot run with more scatterers (could be worth checking...)
+!    + You have different choices for vertical profile of the aerosol layers :
+!           * aeronlay_choice = 1 : Layer tau is spread between ptop and pbot following atm scale height.
+!           * aeronlay_choice = 2 : Layer tau follows its own scale height above cloud deck (pbot).
+!                                   In this case ptop is dummy and sclhght gives the ratio H_cl/H_atm.
+!           * aeronlay_choice = ... feel free to add more cases  !
+!    + Layers can overlap if needed (if you want a 'transition layer' as in the 2-scheme, just add it)
+!
+!=========================================================================================================
+
+      if (iaero_nlay(1) .ne.0) then
+
+        DO ia=1,nlayaero
+           iaer=iaero_nlay(ia)
+
+!          a. Initialization
+           aerosol(1:ngrid,1:nlayer,iaer)=0.D0
+
+!          b. Opacity calculation
+           
+           ! Case 1 : Follows atmospheric scale height between boundaries pressures
+           ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+           IF (aeronlay_choice(ia).eq.1) THEN
+
+             dp_layer(:)=0.D0
+             DO ig=1,ngrid
+               DO l=1,nlayer-1
+                 !! i. Opacity follows scale height
+                 IF ( pplev(ig,l).le.aeronlay_pbot(ia)   .AND.          &
+                      pplev(ig,l).ge.aeronlay_ptop(ia) ) THEN
+                   aerosol(ig,l,iaer) = ( pplev(ig,l) - pplev(ig,l+1) )
+                   dp_layer(ig) = dp_layer(ig) + aerosol(ig,l,iaer)
+                 !! ii. Outside aerosol layer boundaries: no aerosols
+                 ELSE
+                   aerosol(ig,l,iaer) = 0.D0
+                 ENDIF
+               ENDDO
+             ENDDO
+             ! iii. Re-normalize to required total opacity
+             DO ig=1,ngrid
+               DO l=1,nlayer-1
+                 IF ( pplev(ig,l).le.aeronlay_pbot(ia)   .AND.          &
+                      pplev(ig,l).ge.aeronlay_ptop(ia) ) THEN
+                  aerosol(ig,l,iaer) = aerosol(ig,l,iaer) / dp_layer(ig) &
+                                     * aeronlay_tauref(ia)
+                 ENDIF
+               ENDDO
+             ENDDO
+
+           ! Case 2 : Follows input scale height
+           ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+           ELSE IF (aeronlay_choice(ia).eq.2) THEN
+           
+             cloud_slope  = 1.D0/aeronlay_sclhght(ia)
+             pcloud_deck  = aeronlay_pbot(ia)
+             dp_layer(:)  = 0.D0
+
+             DO ig=1,ngrid
+               DO l=1,nlayer-1
+                 !! i. Below cloud layer: no opacity
+                 IF (pplev(ig,l) .gt. pcloud_deck) THEN
+                   aerosol(ig,l,iaer) = 0.D0            
+                 !! ii. Follows scale height above cloud deck
+                 ELSEIF (pplev(ig,l) .le. pcloud_deck) THEN 
+                   aerosol(ig,l,iaer) = ((pplev(ig,l)/pcloud_deck)**(cloud_slope))
+                   dp_layer(ig) = dp_layer(ig) + aerosol(ig,l,iaer)
+                 ENDIF
+               ENDDO
+             ENDDO
+             ! iii. Re-normalize to required total opacity
+             DO ig=1,ngrid
+               DO l=1,nlayer-1
+                 IF (pplev(ig,l) .le. pcloud_deck) THEN 
+                  aerosol(ig,l,iaer) = aerosol(ig,l,iaer) / dp_layer(ig) &
+                                     * aeronlay_tauref(ia)
+                 ENDIF
+               ENDDO
+             ENDDO
+
+           ENDIF ! aeronlay_choice
+
+          ENDDO ! loop on n aerosol layers
+
+      end if ! if N-layer aerosols
+  
 !==================================================================
 !    Jovian auroral aerosols (unknown composition) NON-GENERIC: vertical and meridional profile tuned to observations
@@ -506,5 +611,5 @@
 	 
  !       3. Meridional distribution, and re-normalize to observed total column
-         tau_col(:)=0.D0
+         dp_layer(:)=0.D0
          DO ig=1,ngrid
 	  !!Jupiter
@@ -528,5 +633,5 @@
 
  	  DO l=1,nlayer  
-               tau_col(ig) = tau_col(ig) + aerosol(ig,l,iaer)/obs_tau_col_aurora
+               dp_layer(ig) = dp_layer(ig) + aerosol(ig,l,iaer)/obs_tau_col_aurora
           ENDDO
          ENDDO
@@ -534,5 +639,5 @@
          DO ig=1,ngrid
            DO l=1,nlayer-1
-                aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/tau_col(ig)
+                aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/dp_layer(ig)
            ENDDO
          ENDDO
Index: trunk/LMDZ.GENERIC/libf/phystd/aerosol_mod.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/aerosol_mod.F90	(revision 2283)
+++ trunk/LMDZ.GENERIC/libf/phystd/aerosol_mod.F90	(revision 2297)
@@ -18,7 +18,9 @@
  ! NH3 cloud
       integer :: iaero_nh3 = 0
- ! Auroral aerosols
+! N-layer aerosol model (replaces the 2-layer and hard-coded clouds)
+      integer,dimension(:),allocatable,save :: iaero_nlay
+! Auroral aerosols
       integer :: iaero_aurora = 0
-!$OMP THREADPRIVATE(iaero_co2,iaero_h2o,iaero_dust,iaero_h2so4,noaero,iaero_back2lay,iaero_nh3,iaero_aurora)
+!$OMP THREADPRIVATE(iaero_co2,iaero_h2o,iaero_dust,iaero_h2so4,noaero,iaero_back2lay,iaero_nh3,iaero_nlay,iaero_aurora)
       
 !==================================================================
Index: trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90	(revision 2283)
+++ trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90	(revision 2297)
@@ -28,5 +28,5 @@
       use gases_h, only: ngasmx
       use radii_mod, only : su_aer_radii,co2_reffrad,h2o_reffrad,dust_reffrad,h2so4_reffrad,back2lay_reffrad
-      use aerosol_mod, only : iaero_co2,iaero_h2o,iaero_dust,iaero_h2so4, iaero_back2lay, iaero_nh3, iaero_aurora
+      use aerosol_mod, only : iaero_co2,iaero_h2o,iaero_dust,iaero_h2so4, iaero_back2lay, iaero_aurora
       use tracer_h, only: igcm_h2o_vap, igcm_h2o_ice, igcm_co2_ice
       use comcstfi_mod, only: pi, mugaz, cpp
@@ -82,5 +82,5 @@
       
       ! OUTPUT
-      REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! Aerosol tau.
+      REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! Aerosol tau at reference wavelenght.
       REAL,INTENT(OUT) :: dtlw(ngrid,nlayer)             ! Heating rate (K/s) due to LW radiation.
       REAL,INTENT(OUT) :: dtsw(ngrid,nlayer)             ! Heating rate (K/s) due to SW radiation.
@@ -164,8 +164,8 @@
 
       ! Local aerosol optical properties for each column on RADIATIVE grid.
-      real*8,save,allocatable ::  QXVAER(:,:,:)
+      real*8,save,allocatable ::  QXVAER(:,:,:) ! Extinction coeff (QVISsQREF*QREFvis)
       real*8,save,allocatable ::  QSVAER(:,:,:)
       real*8,save,allocatable ::  GVAER(:,:,:)
-      real*8,save,allocatable ::  QXIAER(:,:,:)
+      real*8,save,allocatable ::  QXIAER(:,:,:) ! Extinction coeff (QIRsQREF*QREFir)
       real*8,save,allocatable ::  QSIAER(:,:,:)
       real*8,save,allocatable ::  GIAER(:,:,:)
@@ -351,7 +351,7 @@
             call back2lay_reffrad(ngrid,reffrad(1,1,iaero_back2lay),nlayer,pplev)
          endif
-!         if(iaer.eq.iaero_nh3)then
-!	    call nh3_reffrad(ngrid,nlayer,reffrad(1,1,iaero_nh3))
-!         endif
+
+         !  For n-layer aerosol size set once for all at firstcall in su_aer_radii
+
 !         if(iaer.eq.iaero_aurora)then
 !	    call aurora_reffrad(ngrid,nlayer,reffrad(1,1,iaero_aurora))
@@ -511,4 +511,6 @@
                pweight=(pplay(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))/   &
                        (pplev(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))
+               ! As 'aerosol' is at reference (visible) wavelenght we scale it as
+               ! it will be multplied by qxi/v in optci/v
                temp=aerosol(ig,L_NLAYRAD-k,iaer)/QREFvis3d(ig,L_NLAYRAD-k,iaer)
                tauaero(2*k+2,iaer)=max(temp*pweight,0.d0)
Index: trunk/LMDZ.GENERIC/libf/phystd/callkeys_mod.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/callkeys_mod.F90	(revision 2283)
+++ trunk/LMDZ.GENERIC/libf/phystd/callkeys_mod.F90	(revision 2297)
@@ -43,6 +43,6 @@
       logical,save :: aeroco2,aeroh2o,aeroh2so4,aeroback2lay
 !$OMP THREADPRIVATE(aeroco2,aeroh2o,aeroh2so4,aeroback2lay)
-      logical,save :: aeronh3, aeroaurora
-!$OMP THREADPRIVATE(aeronh3,aeroaurora)
+      logical,save :: aeronh3, aeronlay, aeroaurora
+!$OMP THREADPRIVATE(aeronh3,aeronlay,aeroaurora)
       logical,save :: aerofixco2,aerofixh2o
 !$OMP THREADPRIVATE(aerofixco2,aerofixh2o)
@@ -62,10 +62,19 @@
 !$OMP THREADPRIVATE(photochem)
 
-      integer,save :: versH2H2cia
       integer,save :: iddist
       integer,save :: iaervar
       integer,save :: iradia
       integer,save :: startype
-!$OMP THREADPRIVATE(versH2H2cia,iddist,iaervar,iradia,startype)
+      integer,save :: versH2H2cia
+      integer,save :: nlayaero
+!$OMP THREADPRIVATE(iddist,iaervar,iradia,startype,versH2H2cia,nlayaero)
+      integer,dimension(:),allocatable,save :: aeronlay_choice
+!$OMP THREADPRIVATE(aeronlay_choice)
+
+      character(64),save :: optprop_back2lay_vis
+      character(64),save :: optprop_back2lay_ir
+      character(64),dimension(:),allocatable,save :: optprop_aeronlay_vis
+      character(64),dimension(:),allocatable,save :: optprop_aeronlay_ir
+!$OMP THREADPRIVATE(optprop_back2lay_vis,optprop_back2lay_ir,optprop_aeronlay_vis,optprop_aeronlay_ir)
 
       real,save :: tplanckmin
@@ -84,6 +93,4 @@
       real,save :: obs_tau_col_strato
 !$OMP THREADPRIVATE(Tstrat,tplanet,obs_tau_col_tropo,obs_tau_col_strato)
-      character(64),save :: optprop_back2lay_vis
-      character(64),save :: optprop_back2lay_ir
       real,save :: pres_bottom_tropo
       real,save :: pres_top_tropo
@@ -101,4 +108,11 @@
       real,save :: tau_nh3_cloud
 !$OMP THREADPRIVATE(size_nh3_cloud, pres_nh3_cloud, tau_nh3_cloud)
+      real,dimension(:),allocatable,save :: aeronlay_tauref
+      real,dimension(:),allocatable,save :: aeronlay_lamref
+      real,dimension(:),allocatable,save :: aeronlay_ptop
+      real,dimension(:),allocatable,save :: aeronlay_pbot
+      real,dimension(:),allocatable,save :: aeronlay_sclhght
+      real,dimension(:),allocatable,save :: aeronlay_size
+!$OMP THREADPRIVATE(aeronlay_tauref,aeronlay_lamref,aeronlay_ptop,aeronlay_pbot,aeronlay_sclhght,aeronlay_size)
       real,save :: co2supsat
       real,save :: pceil
Index: trunk/LMDZ.GENERIC/libf/phystd/iniaerosol.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/iniaerosol.F	(revision 2283)
+++ trunk/LMDZ.GENERIC/libf/phystd/iniaerosol.F	(revision 2297)
@@ -5,5 +5,5 @@
       use aerosol_mod
       use callkeys_mod, only: aeroco2,aeroh2o,dusttau,aeroh2so4,
-     &		aeroback2lay,aeronh3, aeroaurora
+     &		aeroback2lay,aeronh3, nlayaero, aeronlay, aeroaurora
 
       IMPLICIT NONE
@@ -19,5 +19,9 @@
 c=======================================================================
 
-      integer ia
+      integer i, ia
+
+      ! Special case, dyn. allocation for n-layer depending on callphys.def
+      IF(.NOT.ALLOCATED(iaero_nlay)) ALLOCATE(iaero_nlay(nlayaero))
+      iaero_nlay(:) = 0
 
       ia=0
@@ -58,4 +62,12 @@
       write(*,*) '--- NH3 Cloud = ', iaero_nh3
 
+      if (aeronlay) then
+         do i=1,nlayaero
+           ia=ia+1
+           iaero_nlay(i)=ia
+         enddo
+      endif
+      write(*,*) '--- N-layer aerosol = ', iaero_nlay
+
       if (aeroaurora) then
          ia=ia+1
Index: trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90	(revision 2283)
+++ trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90	(revision 2297)
@@ -554,6 +554,13 @@
      call getin_p("aeroh2so4",aeroh2so4)
      write(*,*)" aeroh2so4 = ",aeroh2so4
-	 
+ 
+     write(*,*)"Radiatively active auroral aerosols?"
+     aeroaurora=.false.     ! default value
+     call getin_p("aeroaurora",aeroaurora)
+     write(*,*)" aeroaurora = ",aeroaurora
+
 !=================================
+! TWOLAY scheme and NH3 cloudare left for retrocompatibility only,
+! You should now use N-LAYER scheme (see below).
 
      write(*,*)"Radiatively active two-layer aerosols?"
@@ -562,4 +569,9 @@
      write(*,*)" aeroback2lay = ",aeroback2lay
 
+     if (aeroback2lay) then
+       print*,'Warning : The TWOLAY AEROSOL scheme is deprecated and buggy...'
+       print*,'You should use the generic n-layer scheme (see aeronlay).'
+     endif
+
      write(*,*)"Radiatively active ammonia cloud?"
      aeronh3=.false.     ! default value
@@ -567,8 +579,8 @@
      write(*,*)" aeronh3 = ",aeronh3
 
-     write(*,*)"Radiatively active auroral aerosols?"
-     aeroaurora=.false.     ! default value
-     call getin_p("aeroaurora",aeroaurora)
-     write(*,*)" aeroaurora = ",aeroaurora
+     if (aeronh3) then
+       print*,'Warning : You are using specific NH3 cloud scheme ...'
+       print*,'You should use the generic n-layer scheme (see aeronlay).'
+     endif
 
      write(*,*)"TWOLAY AEROSOL: total optical depth ", &
@@ -613,5 +625,5 @@
        print*,'Error : TWOLAY AEROSOL, Please ensure that in callphys.def'
        print*,'you have pres_top_tropo > pres_bottom_strato !'
-       stop
+       call abort
      endif
 
@@ -647,4 +659,83 @@
      call getin_p("size_nh3_cloud",size_nh3_cloud)
      write(*,*)" size_nh3_cloud = ",size_nh3_cloud
+
+!=================================
+! Generic N-LAYER aerosol scheme
+
+     write(*,*)"Radiatively active generic n-layer aerosols?"
+     aeronlay=.false.     ! default value
+     call getin_p("aeronlay",aeronlay)
+     write(*,*)" aeronlay = ",aeronlay
+
+     write(*,*)"Number of generic aerosols layers?"
+     nlayaero=1     ! default value
+     call getin_p("nlayaero",nlayaero)
+     ! Avoid to allocate arrays of size 0
+     if (aeronlay .and. nlayaero.lt.1) then
+       print*, " You are trying to set no generic aerosols..."
+       print*, " Set aeronlay=.false. instead ! I abort."
+       call abort
+     endif
+     if (.not. aeronlay) nlayaero=1
+     write(*,*)" nlayaero = ",nlayaero
+
+     ! This is necessary, we just set the number of aerosol layers
+     IF(.NOT.ALLOCATED(aeronlay_tauref))      ALLOCATE(aeronlay_tauref(nlayaero))     
+     IF(.NOT.ALLOCATED(aeronlay_lamref))      ALLOCATE(aeronlay_lamref(nlayaero))     
+     IF(.NOT.ALLOCATED(aeronlay_choice))      ALLOCATE(aeronlay_choice(nlayaero))     
+     IF(.NOT.ALLOCATED(aeronlay_pbot))        ALLOCATE(aeronlay_pbot(nlayaero))     
+     IF(.NOT.ALLOCATED(aeronlay_ptop))        ALLOCATE(aeronlay_ptop(nlayaero))     
+     IF(.NOT.ALLOCATED(aeronlay_sclhght))     ALLOCATE(aeronlay_sclhght(nlayaero))     
+     IF(.NOT.ALLOCATED(aeronlay_size))        ALLOCATE(aeronlay_size(nlayaero))     
+     IF(.NOT.ALLOCATED(optprop_aeronlay_ir))  ALLOCATE(optprop_aeronlay_ir(nlayaero))     
+     IF(.NOT.ALLOCATED(optprop_aeronlay_vis)) ALLOCATE(optprop_aeronlay_vis(nlayaero))     
+
+     write(*,*)"Generic n-layer aerosols: Optical depth at reference wavelenght"
+     aeronlay_tauref=1.0E-1
+     call getin_p("aeronlay_tauref",aeronlay_tauref)
+     write(*,*)" aeronlay_tauref = ",aeronlay_tauref
+
+     write(*,*)"Generic n-layer aerosols: Reference wavelenght for optical depths (m)"
+     aeronlay_lamref=0.6E-6
+     call getin_p("aeronlay_lamref",aeronlay_lamref)
+     write(*,*)" aeronlay_lamref = ",aeronlay_lamref
+
+     write(*,*)"Generic n-layer aerosols: Vertical profile choice : &
+                     (1) Tau btwn ptop and pbot follows atm. scale height & 
+                 or  (2) Tau above pbot follows its own scale height"
+     aeronlay_choice=1
+     call getin_p("aeronlay_choice",aeronlay_choice)
+     write(*,*)" aeronlay_choice = ",aeronlay_choice
+
+     write(*,*)"Generic n-layer aerosols: bottom pressures (Pa)"
+     aeronlay_pbot=2000.0
+     call getin_p("aeronlay_pbot",aeronlay_pbot)
+     write(*,*)" aeronlay_pbot = ",aeronlay_pbot
+
+     write(*,*)"Generic n-layer aerosols: (if choice=1) Top pressures (Pa) "
+     aeronlay_ptop=300000.0
+     call getin_p("aeronlay_ptop",aeronlay_ptop)
+     write(*,*)" aeronlay_ptop = ",aeronlay_ptop
+
+     write(*,*)"Generic n-layer aerosols: (if choice=2) Scale height / atm. scale height"
+     aeronlay_ptop=0.2
+     call getin_p("aeronlay_sclhght",aeronlay_sclhght)
+     write(*,*)" aeronlay_sclhght = ",aeronlay_sclhght
+
+     write(*,*)"Generic n-layer aerosols: particles sizes (m)"
+     aeronlay_size=1.e-6
+     call getin_p("aeronlay_size",aeronlay_size)
+     write(*,*)" aeronlay_size = ",aeronlay_size
+
+     write(*,*)"Generic n-layer aerosols: VIS optical properties file"
+     optprop_aeronlay_vis = 'optprop_saturn_vis_n20.dat'
+     call getin_p("optprop_aeronlay_vis",optprop_aeronlay_vis)
+     write(*,*)" optprop_aeronlay_vis = ",optprop_aeronlay_vis
+
+     write(*,*)"Generic n-layer aerosols: IR optical properties file"
+     optprop_aeronlay_ir = 'optprop_saturn_ir_n20.dat'
+     call getin_p("optprop_aeronlay_ir",optprop_aeronlay_ir)
+     write(*,*)" optprop_aeronlay_ir = ",optprop_aeronlay_ir
+     
 
 !=================================
Index: trunk/LMDZ.GENERIC/libf/phystd/radcommon_h.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/radcommon_h.F90	(revision 2283)
+++ trunk/LMDZ.GENERIC/libf/phystd/radcommon_h.F90	(revision 2297)
@@ -57,7 +57,13 @@
 ! 
 ! For the "naerkind" kind of aerosol radiative properties : 
-! QIRsQREF :  Qext / Qext("longrefvis")
+! QIRsQREF :  Qext / Qext("longrefir")
 ! omegaIR  :  mean single scattering albedo
 ! gIR      :  mean assymetry factor
+!
+!
+! Note - QIRsQREF in the martian model was scaled to longrefvis,
+! here it is scaled to longrefir, which is actually a dummy parameter,
+! as we do not output scaled aerosol opacity ...
+!
 
       REAL*8 BWNI(L_NSPECTI+1), WNOI(L_NSPECTI), DWNI(L_NSPECTI), WAVEI(L_NSPECTI) !BWNI read by master in setspi
Index: trunk/LMDZ.GENERIC/libf/phystd/radii_mod.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/radii_mod.F90	(revision 2283)
+++ trunk/LMDZ.GENERIC/libf/phystd/radii_mod.F90	(revision 2297)
@@ -8,7 +8,5 @@
 !     water cloud optical properties
 
-      use callkeys_mod, only: radfixed,Nmix_co2,                    &
-                pres_bottom_tropo,pres_top_tropo,size_tropo,        &
-                pres_bottom_strato,size_strato
+      use callkeys_mod, only: radfixed,Nmix_co2
       
       real, save ::  rad_h2o
@@ -39,6 +37,6 @@
       use radinc_h, only: naerkind
       use aerosol_mod, only: iaero_back2lay, iaero_co2, iaero_dust, &
-                             iaero_h2o, iaero_h2so4,iaero_nh3,iaero_aurora
-      use callkeys_mod, only: size_nh3_cloud
+                             iaero_h2o, iaero_h2so4, iaero_nh3, iaero_nlay, iaero_aurora
+      use callkeys_mod, only: size_nh3_cloud, nlayaero, aeronlay_size
 
       Implicit none
@@ -52,5 +50,5 @@
       logical, save :: firstcall=.true.
 !$OMP THREADPRIVATE(firstcall)
-      integer :: iaer   
+      integer :: iaer, ia   
       
       print*,'enter su_aer_radii'
@@ -60,4 +58,5 @@
 !     a fixed aerosol layer, and be able to define reffrad in a 
 !     .def file. To be improved!
+!                |-> Done in th n-layer aerosol case (JVO 20)
 
             if(iaer.eq.iaero_co2)then ! CO2 ice
@@ -91,4 +90,11 @@
                nueffrad(1:ngrid,1:nlayer,iaer) = 0.1 
             endif
+
+            do ia=1,nlayaero
+              if(iaer.eq.iaero_nlay(ia))then ! N-layer aerosols
+                 reffrad(1:ngrid,1:nlayer,iaer) = aeronlay_size(ia)
+                 nueffrad(1:ngrid,1:nlayer,iaer) = 0.1 
+              endif
+            enddo
 
 	    if(iaer.eq.iaero_aurora)then ! Auroral aerosols
@@ -345,7 +351,8 @@
 !
 !==================================================================
+      use callkeys_mod, only: pres_bottom_tropo,pres_top_tropo,size_tropo,  &
+                              pres_bottom_strato,size_strato
  
-      use aerosol_mod   !! Particle sizes and boundaries of aerosol layers defined there
-     Implicit none
+      Implicit none
 
       integer,intent(in) :: ngrid
Index: trunk/LMDZ.GENERIC/libf/phystd/suaer_corrk.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/suaer_corrk.F90	(revision 2283)
+++ trunk/LMDZ.GENERIC/libf/phystd/suaer_corrk.F90	(revision 2297)
@@ -11,5 +11,7 @@
       use radcommon_h, only: qrefvis,qrefir,omegarefir !,omegarefvis
       use aerosol_mod
-      use callkeys_mod, only: tplanet, optprop_back2lay_vis, optprop_back2lay_ir
+      use callkeys_mod, only: tplanet, optprop_back2lay_vis, optprop_back2lay_ir, &
+                              optprop_aeronlay_vis, optprop_aeronlay_ir,          &
+                              aeronlay_lamref, nlayaero
 
       implicit none
@@ -80,25 +82,41 @@
       REAL epref                       ! reference extinction ep
 
-!      REAL epav(L_NSPECTI)            ! Average ep (= <Qext>/Qext(lamref) if epref=1)
-!      REAL omegav(L_NSPECTI)          ! Average single scattering albedo
-!      REAL gav(L_NSPECTI)             ! Average assymetry parameter
-
-      REAL epavVI(L_NSPECTV)            ! Average ep (= <Qext>/Qext(lamref) if epref=1)
+      REAL epavVI(L_NSPECTV)            ! Average ep (= <Qext>/Qext(lamrefvis) if epref=1)
       REAL omegavVI(L_NSPECTV)          ! Average single scattering albedo
       REAL gavVI(L_NSPECTV)             ! Average assymetry parameter
 
-      REAL epavIR(L_NSPECTI)            ! Average ep (= <Qext>/Qext(lamref) if epref=1)
+      REAL epavIR(L_NSPECTI)            ! Average ep (= <Qext>/Qext(lamrefir) if epref=1)
       REAL omegavIR(L_NSPECTI)          ! Average single scattering albedo
       REAL gavIR(L_NSPECTI)             ! Average assymetry parameter
       
       logical forgetit                  ! use francois' data?
-      integer iwvl
+      integer iwvl, ia
 
 !     Local saved variables:
 
-      CHARACTER(LEN=30), DIMENSION(naerkind,2), SAVE :: file_id
+      CHARACTER(LEN=50), DIMENSION(naerkind,2), SAVE :: file_id
 !$OMP THREADPRIVATE(file_id)
 !---- Please indicate the names of the optical property files below
 !     Please also choose the reference wavelengths of each aerosol      
+
+!--------- README TO UNDERSTAND WHAT FOLLOWS (JVO 20) -------
+!     The lamref variables comes from the Martian model
+!     where the visible one is the one used for computing
+!     and the IR one is used to output scaled opacity to
+!     match instrumental data ... This is done (at least
+!     for now) in the generic, so lamrefir is dummy*!
+
+!     The important one is the VISIBLE one as it will be used
+!     to rescale the values in callcork.F90 assuming 'aerosol' is
+!     at this visible reference wavelenght.
+
+!     *Actually if you change lamrefir here there is a
+!     slight sensitvity in the outputs because of some
+!     machine precision differences that amplifys and lead
+!     up to 10-6 differences in the radiative balance...
+!     This could be good to clean the code but would require
+!     a lot of modifs and to take care !
+
+!--------------------------------------------------------------
 
       if (noaero) then
@@ -125,5 +143,5 @@
            file_id(iaer,2) = 'optprop_co2ice_ir_n50.dat'
         endif
-        lamrefir(iaer)=12.1E-6   ! 825cm-1 TES/MGS ???
+        lamrefir(iaer)=12.1E-6   ! Dummy in generic phys. (JVO 20)
        endif ! CO2 aerosols  
 ! NOTE: these lamref's are currently for dust, not CO2 ice.
@@ -139,5 +157,5 @@
 !     infrared
          file_id(iaer,2) = 'optprop_iceir_n50.dat'
-         lamrefir(iaer)=12.1E-6   ! 825cm-1 TES/MGS
+         lamrefir(iaer)=12.1E-6   ! Dummy in generic phys. (JVO 20)
        endif
 
@@ -151,6 +169,5 @@
 !     infrared
          file_id(iaer,2) = 'optprop_dustir_n50.dat'
-         !lamrefir(iaer)=12.1E-6   ! 825cm-1 TES/MGS
-         lamrefir(iaer)=9.3E-6
+         lamrefir(iaer)=9.3E-6     ! Dummy in generic phys. (JVO 20)
        endif 
 
@@ -160,10 +177,8 @@
 !     visible
          file_id(iaer,1) = 'optprop_h2so4vis_n20.dat'
-         !lamrefir(iaer)=  doesn't exist?
          lamrefvis(iaer)=1.5E-6   ! no idea, must find
 !     infrared
          file_id(iaer,2) = 'optprop_h2so4ir_n20.dat'
-         !lamrefir(iaer)=12.1E-6   ! 825cm-1 TES/MGS
-         lamrefir(iaer)=9.3E-6 ! no idea, must find
+         lamrefir(iaer)=9.3E-6 ! ! Dummy in generic phys. (JVO 20)
 ! added by LK
        endif
@@ -174,8 +189,8 @@
 !     visible
          file_id(iaer,1) = TRIM(optprop_back2lay_vis)
-         lamrefvis(iaer)=0.8E-6  ! 
+         lamrefvis(iaer)=0.8E-6  ! This is the important one.
 !     infrared
          file_id(iaer,2) = TRIM(optprop_back2lay_ir)
-         lamrefir(iaer)=6.E-6  ! 
+         lamrefir(iaer)=6.E-6    ! This is dummy.
 ! added by SG
        endif
@@ -189,8 +204,22 @@
 !     infrared
          file_id(iaer,2) = 'optprop_nh3ice_ir.dat'
-         lamrefir(iaer)=6.E-6  ! 
+         lamrefir(iaer)=6.E-6  ! dummy 
 ! added by SG
        endif
 
+       do ia=1,nlayaero
+         if (iaer.eq.iaero_nlay(ia)) then
+           print*, 'naerkind= nlay', iaer
+           
+!       visible
+           file_id(iaer,1) = TRIM(optprop_aeronlay_vis(ia))
+           lamrefvis(iaer)=aeronlay_lamref(ia)
+!       infrared
+           file_id(iaer,2) = TRIM(optprop_aeronlay_ir(ia))
+           lamrefir(iaer)=6.E-6 ! Dummy value
+         endif
+       enddo
+! added by JVO
+      
        if (iaer.eq.iaero_aurora) then
          print*, 'naerkind= aurora', iaer
@@ -201,5 +230,5 @@
 !     infrared
          file_id(iaer,2) = 'optprop_aurora_ir.dat'
-         lamrefir(iaer)=6.E-6  ! 
+         lamrefir(iaer)=6.E-6  ! dummy 
 ! added by SG
        endif
