Index: LMDZ6/trunk/libf/phylmd/StratAer/strataer_local_var_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/StratAer/strataer_local_var_mod.F90	(revision 4626)
+++ LMDZ6/trunk/libf/phylmd/StratAer/strataer_local_var_mod.F90	(revision 4750)
@@ -1,3 +1,2 @@
-! $Id: strataer_local_var_mod.F90 2022-10-17 19:39:09Z nlebas $
 MODULE strataer_local_var_mod
 ! This module contains strato microphysic model parameters & variables
@@ -41,4 +40,8 @@
   LOGICAL,SAVE :: flag_newclim_file
   !$OMP THREADPRIVATE(flag_newclim_file)
+  
+  ! flag to choose new H2SO4 density and weight percent from Tabazadeh et al. (1994).
+  LOGICAL,SAVE :: flag_new_strat_compo
+  !$OMP THREADPRIVATE(flag_new_strat_compo)
   
   ! Verbose mode to get more print info
@@ -178,4 +181,5 @@
     flag_H2SO4_photolysis = .FALSE. ! H2SO4 photolysis (default: No)
     flag_min_rreduce = .TRUE.       ! Minimum lifetime=1.5 pdt phys (default: Yes)
+    flag_new_strat_compo =.TRUE.    ! H2SO4/H2O weight percent & density routine (default: S. Bekki)
     ok_qemiss = .FALSE.             ! H2O emission flag
     
@@ -208,4 +212,5 @@
     CALL getin_p('flag_H2SO4_photolysis',flag_H2SO4_photolysis)
     CALL getin_p('flag_min_rreduce',flag_min_rreduce)
+    CALL getin_p('flag_new_strat_compo',flag_new_strat_compo)
     CALL getin_p('ok_qemiss',ok_qemiss)
     
@@ -235,4 +240,5 @@
           WRITE(lunout,*) 'IN STRATAER : flag_H2SO4_photolysis = ',flag_H2SO4_photolysis
           WRITE(lunout,*) 'IN STRATAER : flag_min_rreduce = ',flag_min_rreduce
+          WRITE(lunout,*) 'IN STRATAER : flag_new_strat_compo = ',flag_new_strat_compo
           WRITE(lunout,*) 'IN STRATAER : ok_qemiss = ',ok_qemiss
        ENDIF
Index: LMDZ6/trunk/libf/phylmd/StratAer/sulfate_aer_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/StratAer/sulfate_aer_mod.F90	(revision 4626)
+++ LMDZ6/trunk/libf/phylmd/StratAer/sulfate_aer_mod.F90	(revision 4750)
@@ -1,5 +1,2 @@
-!
-! $Id$
-!
 MODULE sulfate_aer_mod
 
@@ -8,4 +5,74 @@
 
 CONTAINS
+
+!*******************************************************************
+  SUBROUTINE STRACOMP_BIN(sh,t_seri,pplay)
+!
+!   Aerosol H2SO4 weight fraction as a function of PH2O and temperature
+!   INPUT:
+!   sh: VMR of H2O
+!   t_seri: temperature (K)
+!   pplay: middle layer pression (Pa)
+!
+!   OUTPUT:
+!   R2SO4: aerosol H2SO4 weight fraction (percent)
+    
+    USE dimphy, ONLY : klon,klev ! nb of longitude and altitude bands
+    USE aerophys
+    USE phys_local_var_mod, ONLY: R2SO4
+    
+    IMPLICIT NONE
+    
+    REAL,DIMENSION(klon,klev),INTENT(IN)          :: t_seri  ! Temperature
+    REAL,DIMENSION(klon,klev),INTENT(IN)          :: pplay   ! pression in the middle of each layer (Pa)
+    REAL,DIMENSION(klon,klev),INTENT(IN)          :: sh      ! specific humidity
+    
+    REAL ks(7)
+    REAL t,qh2o,ptot,pw
+    REAL a,b,c,det
+    REAL xsb,msb
+    
+    INTEGER ilon,ilev
+    DATA ks/-21.661,2724.2,51.81,-15732.0,47.004,-6969.0,-4.6183/
+    
+!*******************************************************************
+!***     liquid aerosols process
+!*******************************************************************
+!        BINARIES LIQUID AEROROLS:
+    
+    DO ilon=1,klon
+       DO ilev=1,klev
+          
+          t = max(t_seri(ilon,ilev),185.)
+          qh2o=sh(ilon,ilev)/18.*28.9
+          ptot=pplay(ilon,ilev)/100.
+          pw = qh2o*ptot/1013.0
+          pw = min(pw,2.e-3/1013.)
+          pw = max(pw,2.e-5/1013.)
+          
+!*******************************************************************
+!***     binaries aerosols h2so4/h2o
+!*******************************************************************
+          a = ks(3) + ks(4)/t
+          b = ks(1) + ks(2)/t
+          c = ks(5) + ks(6)/t + ks(7)*log(t) - log(pw)
+          
+          det = b**2 - 4.*a*c
+          
+          IF (det > 0.) THEN
+             xsb = (-b - sqrt(det))/(2.*a)
+             msb = 55.51*xsb/(1.0 - xsb)
+          ELSE
+             msb = 0.
+          ENDIF
+          R2SO4(ilon,ilev) = 100*msb*0.098076/(1.0 + msb*0.098076)
+          
+          ! H2SO4 min dilution: 0.5%
+          R2SO4(ilon,ilev) = max( R2SO4(ilon,ilev), 0.005 )
+       ENDDO
+    ENDDO
+100 RETURN
+    
+  END SUBROUTINE STRACOMP_BIN
 
 !********************************************************************
@@ -479,4 +546,76 @@
 
 !****************************************************************
+    SUBROUTINE DENH2SA_TABA(t_seri)
+
+!   AERSOL DENSITY AS A FUNCTION OF H2SO4 WEIGHT PERCENT AND T
+!   from Tabazadeh et al. (1994) abaques
+!   ---------------------------------------------
+
+!   
+!   INPUT: 
+!   R2SO4: aerosol H2SO4 weight fraction (percent)
+!   t_seri: temperature (K)
+!   klon: number of latitude bands in the model domain
+!   klev: number of altitude bands in the model domain
+!   for IFS: perhaps add another dimension for longitude
+!
+!   OUTPUT: 
+!   DENSO4: aerosol mass density (gr/cm3 = aerosol mass/aerosol volume)
+!   
+    USE dimphy, ONLY : klon,klev
+    USE phys_local_var_mod, ONLY: R2SO4, DENSO4
+    
+    IMPLICIT NONE
+    
+    REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri  ! Temperature
+        
+    INTEGER i,j
+    
+!----------------------------------------------------------------------
+!	... Local variables
+!----------------------------------------------------------------------
+      real, parameter :: a9 = -268.2616e4, a10 = 576.4288e3
+      
+      real :: a0, a1, a2, a3, a4, a5, a6, a7 ,a8
+      real :: c1, c2, c3, c4, w
+      
+      
+!   Loop on model domain (2 dimension for UPMC model; 3 for IFS)
+    DO i=1,klon
+       DO j=1,klev
+!----------------------------------------------------------------------
+!	... Temperature variables
+!----------------------------------------------------------------------
+          c1 = t_seri(I,J)- 273.15
+          c2 = c1**2
+          c3 = c1*c2
+          c4 = c1*c3
+!----------------------------------------------------------------------
+!	Polynomial Coefficients
+!----------------------------------------------------------------------
+          a0 = 999.8426 + 334.5402e-4*c1 - 569.1304e-5*c2
+          a1 = 547.2659 - 530.0445e-2*c1 + 118.7671e-4*c2 + 599.0008e-6*c3
+          a2 = 526.295e1 + 372.0445e-1*c1 + 120.1909e-3*c2 - 414.8594e-5*c3 + 119.7973e-7*c4
+          a3 = -621.3958e2 - 287.7670*c1 - 406.4638e-3*c2 + 111.9488e-4*c3 + 360.7768e-7*c4
+          a4 = 409.0293e3 + 127.0854e1*c1 + 326.9710e-3*c2 - 137.7435e-4*c3 - 263.3585e-7*c4
+          a5 = -159.6989e4 - 306.2836e1*c1 + 136.6499e-3*c2 + 637.3031e-5*c3
+          a6 = 385.7411e4 + 408.3717e1*c1 - 192.7785e-3*c2
+          a7 = -580.8064e4 - 284.4401e1*c1
+          a8 = 530.1976e4 + 809.1053*c1
+!----------------------------------------------------------------------
+!	... Summation
+!----------------------------------------------------------------------
+!     w : H2SO4 Weight fraction
+          w=r2SO4(i,j)*0.01
+          DENSO4(i,j) = 0.001*(a0 + w*(a1 + w*(a2 + w*(a3 + w*(a4 +  & 
+               w*(a5 + w*(a6 + w*(a7 + w*(a8 + w*(a9 + w*a10))))))))))
+          DENSO4(i,j) = max (0.0, DENSO4(i,j) )
+
+       ENDDO
+    ENDDO
+
+  END SUBROUTINE DENH2SA_TABA
+  
+!****************************************************************
     SUBROUTINE DENH2SA(t_seri)
 
Index: LMDZ6/trunk/libf/phylmd/StratAer/traccoag_mod.F90
===================================================================
--- LMDZ6/trunk/libf/phylmd/StratAer/traccoag_mod.F90	(revision 4626)
+++ LMDZ6/trunk/libf/phylmd/StratAer/traccoag_mod.F90	(revision 4750)
@@ -1,5 +1,2 @@
-!
-! $Id$
-!
 MODULE traccoag_mod
 !
@@ -119,11 +116,22 @@
     WHERE (stratomask.GT.0.5) is_strato=.TRUE. 
 
-! STRACOMP (H2O, P, t_seri -> aerosol composition (R2SO4)) 
-! H2SO4 mass fraction in aerosol (%)
-    CALL stracomp(sh,t_seri,pplay)
-
-! aerosol density (gr/cm3)
-    CALL denh2sa(t_seri)
-
+    IF(flag_new_strat_compo) THEN
+       WRITE(lunout,*) 'traccoag: USE STRAT COMPO from Tabazadeh 1994', flag_new_strat_compo
+       ! STRACOMP (H2O, P, t_seri -> aerosol composition (R2SO4)) : binary routine (from reprobus)
+       ! H2SO4 mass fraction in aerosol (%) from Tabazadeh et al. (1994).
+       CALL stracomp_bin(sh,t_seri,pplay)
+       
+       ! aerosol density (gr/cm3) - from Tabazadeh
+       CALL denh2sa_taba(t_seri)
+    ELSE
+       WRITE(lunout,*) 'traccoag: USE STRAT COMPO from Bekki 2D model', flag_new_strat_compo
+       ! STRACOMP (H2O, P, t_seri -> aerosol composition (R2SO4)) 
+       ! H2SO4 mass fraction in aerosol (%)
+       CALL stracomp(sh,t_seri,pplay)
+       
+       ! aerosol density (gr/cm3)
+       CALL denh2sa(t_seri)
+    ENDIF
+    
 ! compute factor for converting dry to wet radius (for every grid box)
     f_r_wet(:,:) = (dens_aer_dry/(DENSO4(:,:)*1000.)/(R2SO4(:,:)/100.))**(1./3.)
