Index: LMDZ6/branches/LMDZ_ECRad/libf/phylmd/StratAer/interp_sulf_input.F90
===================================================================
--- LMDZ6/branches/LMDZ_ECRad/libf/phylmd/StratAer/interp_sulf_input.F90	(revision 4175)
+++ LMDZ6/branches/LMDZ_ECRad/libf/phylmd/StratAer/interp_sulf_input.F90	(revision 4727)
@@ -12,5 +12,5 @@
   USE mod_phys_lmdz_omp_data, ONLY :  is_omp_root
   USE phys_local_var_mod, ONLY : budg_3D_backgr_ocs, budg_3D_backgr_so2
-  USE phys_local_var_mod, ONLY : OCS_lifetime, SO2_lifetime
+  USE phys_local_var_mod, ONLY : OCS_lifetime, SO2_lifetime, H2SO4_lifetime, O3_clim
   USE mod_phys_lmdz_para 
   USE dimphy
@@ -19,4 +19,5 @@
   USE aerophys
   USE YOMCST
+  USE strataer_local_var_mod, ONLY : flag_newclim_file,flag_verbose_strataer
 
   IMPLICIT NONE
@@ -42,14 +43,14 @@
   REAL paprs_glo(klon_glo,klev+1)
 
-  REAL, POINTER:: latitude(:)
+  REAL, ALLOCATABLE :: latitude(:)
 ! (of input data sorted in strictly ascending order)
 
-  REAL, POINTER:: longitude(:)
+  REAL, ALLOCATABLE :: longitude(:)
 ! (of input data sorted in strictly ascending order)
 
-  REAL, POINTER:: time(:)
+  REAL, ALLOCATABLE :: time(:)
 ! (of input data sorted in strictly ascending order)
 
-  REAL, POINTER:: lev(:)
+  REAL, ALLOCATABLE :: lev(:)
 ! levels of input data
 
@@ -62,4 +63,11 @@
   REAL OCS_clim_glo(klon_glo,klev)
   REAL SO2_clim_glo(klon_glo,klev)
+  
+  REAL, ALLOCATABLE :: O3_clim_in(:, :, :, :)
+  REAL, ALLOCATABLE :: O3_clim_mth(:, :, :)
+  REAL, ALLOCATABLE :: O3_clim_tmp(:, :)
+  REAL O3_clim_glo(klon_glo,klev)
+  
+  ! fixed climatos
   REAL, ALLOCATABLE :: OCS_lifetime_in(:, :, :, :)
   REAL, ALLOCATABLE :: SO2_lifetime_in(:, :, :, :)
@@ -70,4 +78,10 @@
   REAL OCS_lifetime_glo(klon_glo,klev)
   REAL SO2_lifetime_glo(klon_glo,klev)
+  
+  REAL, ALLOCATABLE :: H2SO4_lifetime_in(:, :, :, :)
+  REAL, ALLOCATABLE :: H2SO4_lifetime_mth(:, :, :)
+  REAL, ALLOCATABLE :: H2SO4_lifetime_tmp(:, :)
+
+  REAL H2SO4_lifetime_glo(klon_glo,klev)
 !
   REAL, ALLOCATABLE, SAVE :: OCS_clim(:,:)
@@ -79,6 +93,8 @@
 
 ! For NetCDF:
-  INTEGER ncid_in  ! IDs for input files
-  INTEGER varid, ncerr
+  INTEGER             :: ncid_in  ! IDs for input files
+  INTEGER             :: varid, ncerr
+  CHARACTER (LEN=3)   :: nc_lat, nc_lon
+  CHARACTER (LEN=256) :: nc_fname
     
   INTEGER, PARAMETER :: lev_input=17
@@ -103,44 +119,70 @@
   IF (is_mpi_root.AND.is_omp_root) THEN
 
-!--reading emission files
-    CALL nf95_open("ocs_so2_annual_lmdz.nc", nf90_nowrite, ncid_in)
-
+    !--init ncdf variables
+    IF(flag_newclim_file) THEN
+      nc_fname = "ocs_so2_h2so4_annual_lmdz.nc"
+      nc_lat   = "LAT"
+      nc_lon   = "LON"
+    ELSE
+      ! old file for retro compatibility
+      nc_fname = "ocs_so2_annual_lmdz.nc"
+      nc_lat   = "lat"
+      nc_lon   = "lon"
+    ENDIF
+
+    !--reading emission files
+    CALL nf95_open(nc_fname, nf90_nowrite, ncid_in)
+    IF(flag_verbose_strataer) print *, "Open file ", nc_fname
+    
     CALL nf95_inq_varid(ncid_in, "LEV", varid)
     CALL nf95_gw_var(ncid_in, varid, lev)
     n_lev = size(lev)
-
-    CALL nf95_inq_varid(ncid_in, "lat", varid)
+    
+    CALL nf95_inq_varid(ncid_in, nc_lat, varid)
     CALL nf95_gw_var(ncid_in, varid, latitude)
     n_lat = size(latitude)
-
-    CALL nf95_inq_varid(ncid_in, "lon", varid)
+    
+    CALL nf95_inq_varid(ncid_in, nc_lon, varid)
     CALL nf95_gw_var(ncid_in, varid, longitude)
     n_lon = size(longitude)
-
+    
     CALL nf95_inq_varid(ncid_in, "TIME", varid)
     CALL nf95_gw_var(ncid_in, varid, time)
     n_mth = size(time)
-
+    
     IF (.NOT.ALLOCATED(OCS_clim_in))     ALLOCATE(OCS_clim_in(n_lon, n_lat, n_lev, n_mth))
     IF (.NOT.ALLOCATED(SO2_clim_in))     ALLOCATE(SO2_clim_in(n_lon, n_lat, n_lev, n_mth))
+    IF (.NOT.ALLOCATED(O3_clim_in))      ALLOCATE(O3_clim_in(n_lon, n_lat, n_lev, n_mth))
+    
     IF (.NOT.ALLOCATED(OCS_lifetime_in)) ALLOCATE(OCS_lifetime_in(n_lon, n_lat, n_lev, n_mth))
     IF (.NOT.ALLOCATED(SO2_lifetime_in)) ALLOCATE(SO2_lifetime_in(n_lon, n_lat, n_lev, n_mth))
-
+    IF (.NOT.ALLOCATED(H2SO4_lifetime_in)) ALLOCATE(H2SO4_lifetime_in(n_lon, n_lat, n_lev, n_mth))
+    
     CALL nf95_inq_varid(ncid_in, "OCS", varid)
     ncerr = nf90_get_var(ncid_in, varid, OCS_clim_in)
-    print *,'code erreur OCS=', ncerr, varid
+    IF(flag_verbose_strataer) print *,'code erreur OCS=', ncerr, varid
 
     CALL nf95_inq_varid(ncid_in, "SO2", varid)
     ncerr = nf90_get_var(ncid_in, varid, SO2_clim_in)
-    print *,'code erreur SO2=', ncerr, varid
+    IF(flag_verbose_strataer) print *,'code erreur SO2=', ncerr, varid
 
     CALL nf95_inq_varid(ncid_in, "OCS_LIFET", varid)
     ncerr = nf90_get_var(ncid_in, varid, OCS_lifetime_in)
-    print *,'code erreur OCS_lifetime_in=', ncerr, varid
+    IF(flag_verbose_strataer) print *,'code erreur OCS_lifetime_in=', ncerr, varid
 
     CALL nf95_inq_varid(ncid_in, "SO2_LIFET", varid)
     ncerr = nf90_get_var(ncid_in, varid, SO2_lifetime_in)
-    print *,'code erreur SO2_lifetime_in=', ncerr, varid
-
+    IF(flag_verbose_strataer) print *,'code erreur SO2_lifetime_in=', ncerr, varid
+    
+    IF(flag_newclim_file) THEN
+       CALL nf95_inq_varid(ncid_in, "O3", varid)
+       ncerr = nf90_get_var(ncid_in, varid, O3_clim_in)
+       IF(flag_verbose_strataer) print *,'code erreur O3=', ncerr, varid
+       
+       CALL nf95_inq_varid(ncid_in, "H2SO4_LIFET", varid)
+       ncerr = nf90_get_var(ncid_in, varid, H2SO4_lifetime_in)
+       IF(flag_verbose_strataer) print *,'code erreur H2SO4_lifetime_in=', ncerr, varid
+    ENDIF
+    
     CALL nf95_close(ncid_in)
 
@@ -149,9 +191,14 @@
     IF (.NOT.ALLOCATED(OCS_clim_tmp)) ALLOCATE(OCS_clim_tmp(klon_glo, n_lev))
     IF (.NOT.ALLOCATED(SO2_clim_tmp)) ALLOCATE(SO2_clim_tmp(klon_glo, n_lev))
+    IF (.NOT.ALLOCATED(O3_clim_mth))  ALLOCATE(O3_clim_mth(n_lon, n_lat, n_lev))
+    IF (.NOT.ALLOCATED(O3_clim_tmp))  ALLOCATE(O3_clim_tmp(klon_glo, n_lev))
+    
     IF (.NOT.ALLOCATED(OCS_lifetime_mth)) ALLOCATE(OCS_lifetime_mth(n_lon, n_lat, n_lev))
     IF (.NOT.ALLOCATED(SO2_lifetime_mth)) ALLOCATE(SO2_lifetime_mth(n_lon, n_lat, n_lev))
     IF (.NOT.ALLOCATED(OCS_lifetime_tmp)) ALLOCATE(OCS_lifetime_tmp(klon_glo, n_lev))
     IF (.NOT.ALLOCATED(SO2_lifetime_tmp)) ALLOCATE(SO2_lifetime_tmp(klon_glo, n_lev))
-
+    IF (.NOT.ALLOCATED(H2SO4_lifetime_mth)) ALLOCATE(H2SO4_lifetime_mth(n_lon, n_lat, n_lev))
+    IF (.NOT.ALLOCATED(H2SO4_lifetime_tmp)) ALLOCATE(H2SO4_lifetime_tmp(klon_glo, n_lev))
+    
 !---select the correct month, undo multiplication with 1.e12 (precision reasons)
 !---correct latitudinal order and convert input from volume mixing ratio to mass mixing ratio
@@ -161,4 +208,14 @@
       SO2_lifetime_mth(:,j,:) = SO2_lifetime_in(:,n_lat+1-j,:,mth_cur)
       OCS_lifetime_mth(:,j,:) = OCS_lifetime_in(:,n_lat+1-j,:,mth_cur)
+      
+      ! O3 from 2d model is not tracer, in VMR
+      IF(flag_newclim_file) THEN
+         H2SO4_lifetime_mth(:,j,:) = H2SO4_lifetime_in(:,n_lat+1-j,:,mth_cur)
+         ! new input files
+         O3_clim_mth(:,j,:) = 1.e-6*O3_clim_in(:,n_lat+1-j,:,mth_cur)
+      ELSE
+         H2SO4_lifetime_mth(:,j,:) = 1.e-6
+         O3_clim_mth(:,j,:) = 1.e-6
+      ENDIF
     ENDDO
 
@@ -166,10 +223,25 @@
     CALL grid2dTo1d_glo(OCS_clim_mth,OCS_clim_tmp)
     CALL grid2dTo1d_glo(SO2_clim_mth,SO2_clim_tmp)
+    CALL grid2dTo1d_glo(O3_clim_mth,O3_clim_tmp)
+    
     CALL grid2dTo1d_glo(OCS_lifetime_mth,OCS_lifetime_tmp)
     CALL grid2dTo1d_glo(SO2_lifetime_mth,SO2_lifetime_tmp)
+    CALL grid2dTo1d_glo(H2SO4_lifetime_mth,H2SO4_lifetime_tmp)
 
 !--set lifetime to very high value in uninsolated areas
     DO i=1, klon_glo
       DO kk=1, n_lev
+!--added tests on VMR of species for realism
+        IF (OCS_clim_tmp(i,kk).LT.1.e-20) THEN
+          OCS_clim_tmp(i,kk)=1.0e-20
+        ENDIF
+        IF (SO2_clim_tmp(i,kk).LT.1.e-20) THEN
+          SO2_clim_tmp(i,kk)=1.0e-20
+        ENDIF
+        !       50 ppbv min for O3
+        IF (O3_clim_tmp(i,kk).LT.50.0e-9) THEN
+           O3_clim_tmp(i,kk)=50.0e-9
+        ENDIF
+        
         IF (OCS_lifetime_tmp(i,kk)==0.0) THEN
           OCS_lifetime_tmp(i,kk)=1.0e12
@@ -178,4 +250,7 @@
           SO2_lifetime_tmp(i,kk)=1.0e12
         ENDIF
+       IF (H2SO4_lifetime_tmp(i,kk)==0.0) THEN
+          H2SO4_lifetime_tmp(i,kk)=1.0e12
+       ENDIF
       ENDDO
     ENDDO
@@ -186,6 +261,10 @@
      OCS_lifetime_glo(i,k)=0.0
      SO2_lifetime_glo(i,k)=0.0
+     O3_clim_glo(i,k)=0.0
+     
      OCS_clim_glo(i,k)=0.0
      SO2_clim_glo(i,k)=0.0
+     H2SO4_lifetime_glo(i,k)=0.0
+     
      DO kk=1, n_lev
       OCS_lifetime_glo(i,k)=OCS_lifetime_glo(i,k)+ &
@@ -195,4 +274,11 @@
            MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk))-MAX(paprs_glo(i,k+1),paprs_input(kk+1))) &
            *SO2_lifetime_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1))
+      IF(flag_newclim_file) THEN
+         H2SO4_lifetime_glo(i,k)=H2SO4_lifetime_glo(i,k)+ &
+              MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk)) &
+              -MAX(paprs_glo(i,k+1),paprs_input(kk+1))) &
+              *H2SO4_lifetime_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1))
+      ENDIF
+      
       OCS_clim_glo(i,k)=OCS_clim_glo(i,k)+ &
            MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk))-MAX(paprs_glo(i,k+1),paprs_input(kk+1))) &
@@ -201,4 +287,10 @@
            MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk))-MAX(paprs_glo(i,k+1),paprs_input(kk+1))) &
            *SO2_clim_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1))
+      IF(flag_newclim_file) THEN
+         O3_clim_glo(i,k)=O3_clim_glo(i,k)+ &
+              MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk)) &
+              -MAX(paprs_glo(i,k+1),paprs_input(kk+1))) &
+              *O3_clim_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1))
+      ENDIF
       ENDDO
     ENDDO
@@ -213,15 +305,17 @@
   CALL scatter(OCS_clim_glo, OCS_clim)
   CALL scatter(SO2_clim_glo, SO2_clim)
+  CALL scatter(O3_clim_glo, O3_clim)
   CALL scatter(OCS_lifetime_glo, OCS_lifetime)
   CALL scatter(SO2_lifetime_glo, SO2_lifetime)
-
+  CALL scatter(H2SO4_lifetime_glo, H2SO4_lifetime)
+  
   IF (is_mpi_root.AND.is_omp_root) THEN
 !
-    DEALLOCATE(OCS_clim_in,SO2_clim_in)
-    DEALLOCATE(OCS_clim_mth,SO2_clim_mth)
-    DEALLOCATE(OCS_clim_tmp,SO2_clim_tmp)
-    DEALLOCATE(OCS_lifetime_in,SO2_lifetime_in)
-    DEALLOCATE(OCS_lifetime_mth,SO2_lifetime_mth)
-    DEALLOCATE(OCS_lifetime_tmp,SO2_lifetime_tmp)
+    DEALLOCATE(OCS_clim_in,SO2_clim_in,O3_clim_in)
+    DEALLOCATE(OCS_clim_mth,SO2_clim_mth,O3_clim_mth)
+    DEALLOCATE(OCS_clim_tmp,SO2_clim_tmp,O3_clim_tmp)
+    DEALLOCATE(OCS_lifetime_in, SO2_lifetime_in, H2SO4_lifetime_in)
+    DEALLOCATE(OCS_lifetime_mth, SO2_lifetime_mth, H2SO4_lifetime_mth)
+    DEALLOCATE(OCS_lifetime_tmp, SO2_lifetime_tmp, H2SO4_lifetime_tmp)
 !
   ENDIF !--is_mpi_root and is_omp_root
