Index: /LMDZ6/branches/DYNAMICO-conv/DefLists/context_lmdz.xml
===================================================================
--- /LMDZ6/branches/DYNAMICO-conv/DefLists/context_lmdz.xml	(revision 3412)
+++ /LMDZ6/branches/DYNAMICO-conv/DefLists/context_lmdz.xml	(revision 3413)
@@ -12,18 +12,35 @@
   <file_definition src="./file_def_histhf_lmdz.xml"/>
   <file_definition src="./file_def_histins_lmdz.xml"/>
-<!--  <file_definition src="./file_def_histLES_lmdz.xml"/> -->
+  <file_definition src="./file_def_histLES_lmdz.xml"/>
   <file_definition src="./file_def_histmth_lmdz.xml"/>
-<!--  <file_definition src="./file_def_histstn_lmdz.xml"/> -->
+  <file_definition src="./file_def_histstn_lmdz.xml"/>
   <file_definition src="./file_def_histmthNMC_lmdz.xml"/>
   <file_definition src="./file_def_histdayNMC_lmdz.xml"/>
   <file_definition src="./file_def_histhfNMC_lmdz.xml"/>
-<!--  <file_definition src="./file_def_histmthCOSP_lmdz.xml"/> -->
-<!--  <file_definition src="./file_def_histdayCOSP_lmdz.xml"/> -->
-<!--  <file_definition src="./file_def_histhfCOSP_lmdz.xml"/> -->
+  <file_definition src="./file_def_histmthCOSP_lmdz.xml"/> -->
+  <file_definition src="./file_def_histdayCOSP_lmdz.xml"/>  -->
+  <file_definition src="./file_def_histhfCOSP_lmdz.xml"/> --> 
   <file_definition src="./file_def_histstrataer_lmdz.xml"/>
- 
+
+
+
+
+
   <!-- Define domains and groups of domains -->
   <domain_definition>
-    <domain id="dom_glo" data_dim="2" />
+
+    <domain id="dom_glo"  data_dim="1" />
+
+    <domain id="greordered"  domain_ref="dom_glo">
+      <reorder_domain invert_lat="true" shift_lon_fraction="0.5" min_lon="0" max_lon="360" />
+    </domain> 
+   
+    <domain id="dom_regular" ni_glo="144" nj_glo="142" type="rectilinear"  >
+      <generate_rectilinear_domain/>
+      <interpolate_domain order="1"/>
+    </domain>
+   
+    <domain id="dom_out" domain_ref="dom_glo"/>
+    
   </domain_definition>
   
@@ -36,5 +53,5 @@
   <grid_definition>
     <grid id="grid_scalar" >
-    <scalar/>
+      <scalar/>
     </grid>
   </grid_definition>
@@ -42,96 +59,132 @@
   <!-- Define groups of vertical axes -->
   <axis_definition>
-    <axis id="presnivs" standard_name="Vertical levels" unit="Pa">
+    <axis id="time_month" n_glo="12" value="(0,11) [1 2 3 4 5 6 7 8 9 10 11 12]"/>
+    <axis id="time_year" unit="day" />
+    <axis id="presnivs" standard_name="Vertical levels" unit="Pa"/>
+    <axis id="Ahyb" standard_name="Ahyb comp of Hyb Cord" unit="Pa"/>
+    <axis id="Bhyb" standard_name="Bhyb comp of Hyb Cord" unit=""/>
+    <axis id="Ahyb_inter" standard_name="A comp of Hyb Cord at interface" unit="Pa"/>
+    <axis id="Bhyb_inter" standard_name="B comp of Hyb Cord at interface" unit=""/>
+    <axis id="Alt" standard_name="Height approx for scale heigh of 8km at levels" unit="km"/>
+    <axis id="plev" standard_name="model_level_number" unit="Pa"/>
+    <axis id="klev"  prec="8" long_name="number of layers" standard_name="number of layers" unit="1" />
+    <axis id="klevp1"  prec="8" long_name="number of layer interfaces" standard_name="number of layer interfaces" unit="1" />
+    <axis id="bnds" standard_name="bounds" unit="1" />
+    <axis id="spectband" standard_name="Sensor Band Central Radiation Wavenumber" unit="m-1"/>
+    <axis id="axis_lat" standard_name="Latitude axis">
+      <reduce_domain operation="average" direction="iDir" />
     </axis>
-    <axis id="Ahyb" standard_name="Ahyb comp of Hyb Cord" unit="Pa">
-    </axis> 
-    <axis id="Bhyb" standard_name="Bhyb comp of Hyb Cord" unit="">
-    </axis>
-    <axis id="Ahyb_inter" standard_name="A comp of Hyb Cord at interface" unit="Pa">
-    </axis> 
-    <axis id="Bhyb_inter" standard_name="B comp of Hyb Cord at interface" unit="">
-    </axis>
-    <axis id="Alt" standard_name="Height approx for scale heigh of 8km at levels" unit="km">
-    </axis>
-    <axis id="plev" standard_name="model_level_number" unit="Pa">
-    </axis>
-    <axis id="klev"  prec="8" long_name = "hybrid sigma pressure coordinate"
-          standard_name ="atmosphere_hybrid_sigma_pressure_coordinate" unit="1">
-    </axis>
-    <axis id="bnds" standard_name="bounds" unit="1" >
-    </axis>
 
 <!-- Cosp axis definitions-->
-    <axis id="height" standard_name="Cosp levels" unit="m"> 
-    </axis>
-    <axis id="height_mlev" standard_name="height_mlev" unit="m"> 
-    </axis>
-    <axis id="sza" standard_name="solar_zenith_angle" unit="degrees">
-    </axis>
-    <axis id="pressure2" standard_name="pressure" unit="mb">
-    </axis>
-    <axis id="column" standard_name="column" unit="count"> 
-    </axis> 
-    <axis id="temp" standard_name="temperature" unit="K"> 
-    </axis> 
-    <axis id="cth16" name="cth" standard_name="altitude" unit="m"> 
-    </axis> 
-    <axis id="ReffIce" standard_name="ReffIce" unit="microne" >
-    </axis>
-    <axis id="ReffLiq" standard_name="ReffLiq" unit="microne" >
-    </axis>
-    <axis id="scatratio" standard_name="scatratio" unit="1" >
-    </axis>
-    <axis id="dbze" standard_name="dbze" unit="dBZ" >
-    </axis>
-    <axis id="tau" standard_name="tau" unit="1" >
-    </axis>
+    <axis id="height" standard_name="Cosp levels" unit="m"/> 
+    <axis id="height_mlev" standard_name="height_mlev" unit="m"/> 
+    <axis id="sza" standard_name="solar_zenith_angle" unit="degrees"/>
+    <axis id="pressure2" standard_name="pressure" unit="mb"/>
+    <axis id="column" standard_name="column" unit="count"/> 
+    <axis id="temp" standard_name="temperature" unit="K"/> 
+    <axis id="cth16" standard_name="altitude" unit="m"/> 
+    <axis id="ReffIce" standard_name="ReffIce" unit="microne" />
+    <axis id="ReffLiq" standard_name="ReffLiq" unit="microne" />
+    <axis id="scatratio" standard_name="scatratio" unit="1" />
+    <axis id="dbze" standard_name="dbze" unit="dBZ" />
+    <axis id="tau" standard_name="tau" unit="1" />              
   </axis_definition>
 
   <grid_definition>
-<!-- Define Scalar grid for GHG, orbital parameters and solar constants -->
-    <grid id="grid_scalar">
-    </grid>
 
     <grid id="klev_bnds"> <axis axis_ref="klev" /> <axis axis_ref="bnds" /> </grid>
+    <grid id="klevp1_bnds"> <axis axis_ref="klevp1" /> <axis axis_ref="bnds" /> </grid>
+
 
      <grid id="grid_glo">
-        <domain id="dom_glo" />
-     </grid>
+        <domain domain_ref="dom_glo" />
+     </grid>
+
+     <grid id="grid_out">
+        <domain domain_ref="dom_out" />
+     </grid>
+
 
      <grid id="grid_glo_presnivs">
-        <domain id="dom_glo" />
-        <axis id="presnivs" />
+        <domain domain_ref="dom_glo" />
+        <axis axis_ref="presnivs" />
+     </grid>
+
+     <grid id="grid_out_presnivs">
+        <domain domain_ref="dom_out" />
+        <axis axis_ref="presnivs" />
      </grid>
 
 
      <grid id="grid_glo_plev">
-        <domain id="dom_glo" />
-        <axis id="plev" />
-     </grid>
+        <domain domain_ref="dom_glo" />
+        <axis axis_ref="plev" />
+     </grid>
+
+     <grid id="grid_out_plev">
+        <domain domain_ref="dom_out" />
+        <axis axis_ref="plev" />
+     </grid>
+
+
+     <grid id="grid_glo_spectband">
+        <domain domain_ref="dom_glo" />
+        <axis axis_ref="spectband" />
+     </grid>
+
+     <grid id="grid_out_spectband">
+        <domain domain_ref="dom_out" />
+        <axis axis_ref="spectband" />
+     </grid>
+     
 
      <grid id="grid_glo_height">
-        <domain id="dom_glo" />
-        <axis id="height" />
-     </grid>
+        <domain domain_ref="dom_glo" />
+        <axis axis_ref="height" />
+     </grid>
+     <grid id="grid_out_height">
+        <domain domain_ref="dom_out" />
+        <axis axis_ref="height" />
+     </grid>
+
 
      <grid id="grid_glo_heightmlev">
-        <domain id="dom_glo" />
-        <axis id="height_mlev" />
-     </grid>
+        <domain domain_ref="dom_glo" />
+        <axis axis_ref="height_mlev" />
+     </grid>
+     <grid id="grid_out_heightmlev">
+        <domain domain_ref="dom_out" />
+        <axis axis_ref="height_mlev" />
+     </grid>
+
 
      <grid id="grid_glo_temp">
-        <domain id="dom_glo" />
-        <axis id="temp" />
-     </grid>
+        <domain domain_ref="dom_glo" />
+        <axis axis_ref="temp" />
+     </grid>
+     <grid id="grid_out_temp">
+        <domain domain_ref="dom_out" />
+        <axis axis_ref="temp" />
+     </grid>
+
+
 
      <grid id="grid_glo_sza">
-        <domain id="dom_glo" />
-        <axis id="sza" />
-     </grid>
+        <domain domain_ref="dom_glo" />
+        <axis axis_ref="sza" />
+     </grid>
+     <grid id="grid_out_sza">
+        <domain domain_ref="dom_out" />
+        <axis axis_ref="sza" />
+     </grid>
+
 
      <grid id="grid_glo_column">
-        <domain id="dom_glo" />
-        <axis id="column" />
+        <domain domain_ref="dom_glo" />
+        <axis axis_ref="column" />
+     </grid>
+     <grid id="grid_out_column">
+        <domain domain_ref="dom_out" />
+        <axis axis_ref="column" />
      </grid>
 
@@ -139,49 +192,113 @@
 <!-- Define 4D grids for Cosp simulator -->
      <grid id="grid4Dcol">
-        <domain id="dom_glo" />
-        <axis id="height_mlev" />
-        <axis id="column" />
+        <domain domain_ref="dom_glo" />
+        <axis axis_ref="column" />
+        <axis axis_ref="height_mlev" />
+     </grid>
+     <grid id="grid4Dcol_out">
+        <domain domain_ref="dom_out" />
+        <axis axis_ref="column" />
+        <axis axis_ref="height_mlev" />
      </grid>
 
      <grid id="grid4Dsrbin">
-        <domain id="dom_glo" />
-        <axis id="height" />
-        <axis id="scatratio" />
-     </grid>
+        <domain domain_ref="dom_glo" />
+        <axis axis_ref="height" />
+        <axis axis_ref="scatratio" />
+     </grid>
+     <grid id="grid4Dsrbin_out">
+        <domain domain_ref="dom_out" />
+        <axis axis_ref="height" />
+        <axis axis_ref="scatratio" />
+     </grid>
+
 
      <grid id="grid4Ddbze">
-        <domain id="dom_glo" />
-        <axis id="height" />
-        <axis id="dbze" />
+        <domain domain_ref="dom_glo" />
+        <axis axis_ref="height" />
+        <axis axis_ref="dbze" />
+     </grid>
+     <grid id="grid4Ddbze_out">
+        <domain domain_ref="dom_out" />
+        <axis axis_ref="height" />
+        <axis axis_ref="dbze" />
      </grid>
 
      <grid id="grid4Dtau">
-        <domain id="dom_glo" />
-        <axis id="tau" />
-        <axis id="pressure2" />
+        <domain domain_ref="dom_glo" />
+        <axis axis_ref="tau" />
+        <axis axis_ref="pressure2" />
+     </grid>
+     <grid id="grid4Dtau_out">
+        <domain domain_ref="dom_out" />
+        <axis axis_ref="tau" />
+        <axis axis_ref="pressure2" />
      </grid>
 
      <grid id="grid4Dmisr">
-        <domain id="dom_glo" />
-        <axis id="tau" />
-        <axis id="cth16" />
+        <domain domain_ref="dom_glo" />
+        <axis axis_ref="cth16" />
+        <axis axis_ref="tau" />
+     </grid>
+     <grid id="grid4Dmisr_out">
+        <domain domain_ref="dom_out" />
+        <axis axis_ref="cth16" />
+        <axis axis_ref="tau" />
      </grid>
 
      <grid id="grid4Dreffi">
-        <domain id="dom_glo" />
-        <axis id="tau" />
-        <axis id="ReffIce" />
-     </grid>
-
+        <domain domain_ref="dom_glo" />
+        <axis axis_ref="tau" />
+        <axis axis_ref="ReffIce" />
+     </grid>
+     <grid id="grid4Dreffi_out">
+        <domain domain_ref="dom_out" />
+        <axis axis_ref="tau" />
+        <axis axis_ref="ReffIce" />
+     </grid>
+
+     <grid id="grid4Dreffl_out">
+        <domain domain_ref="dom_glo" />
+        <axis axis_ref="tau" />
+        <axis axis_ref="ReffLiq" />
+     </grid>
      <grid id="grid4Dreffl">
-        <domain id="dom_glo" />
-        <axis id="tau" />
-        <axis id="ReffLiq" />
-     </grid>
+        <domain domain_ref="dom_out" />
+        <axis axis_ref="tau" />
+        <axis axis_ref="ReffLiq" />
+     </grid>
+
+
      <grid id="grid4Dcol2">
-        <domain id="dom_glo" />
-        <axis id="height" />
-        <axis id="column" />
-     </grid>
+        <domain domain_ref="dom_glo" />
+        <axis axis_ref="height" />
+        <axis axis_ref="column" />
+     </grid>
+     <grid id="grid4Dcol2_out">
+        <domain domain_ref="dom_out" />
+        <axis axis_ref="height" />
+        <axis axis_ref="column" />
+     </grid>
+
+ <!-- Grid definitions to allow summing of a 3D variable -->    
+      <grid id="grid_3D" >
+         <domain domain_ref="dom_glo" />
+         <axis axis_ref="lev" />
+       </grid>
+
+       <grid id="grid_sum_axis">
+         <domain domain_ref="dom_glo" />
+         <scalar>
+            <reduce_axis operation="sum" />
+         </scalar>
+       </grid>
+
+       <grid id="grid_sum">
+         <scalar >
+           <reduce_domain operation="sum" local="true" />
+           <reduce_scalar operation="sum" />
+         </scalar>
+         <scalar/>
+       </grid>
 
   </grid_definition>
Index: /LMDZ6/branches/DYNAMICO-conv/bld.cfg
===================================================================
--- /LMDZ6/branches/DYNAMICO-conv/bld.cfg	(revision 3412)
+++ /LMDZ6/branches/DYNAMICO-conv/bld.cfg	(revision 3413)
@@ -66,5 +66,5 @@
 # Example of how to set specific compiling options for a specific file
 # -> this can be including in the arch.opt file
-#bld::tool::fflags::phys::readaerosol         %BASE_FFLAGS %PROD_FFLAGS  %INCDIR -C hopt -pi auto
+#bld::tool::fflags::phys::physiq_mod         %BASE_FFLAGS %DEBUG_FFLAGS  %INCDIR
 
 
Index: /LMDZ6/branches/DYNAMICO-conv/build_gcm
===================================================================
--- /LMDZ6/branches/DYNAMICO-conv/build_gcm	(revision 3412)
+++ /LMDZ6/branches/DYNAMICO-conv/build_gcm	(revision 3413)
@@ -45,6 +45,11 @@
 # run "fcm build" command
 ${dirname}fcm build $*
+status=$?
 
 # cleanup
 \rm -f '.lock' 
 
+if [ $status -ne 0 ]; then
+  exit $status
+fi
+
Index: /LMDZ6/branches/DYNAMICO-conv/libf/dynphy_lonlat/phylmd/etat0dyn_netcdf.F90
===================================================================
--- /LMDZ6/branches/DYNAMICO-conv/libf/dynphy_lonlat/phylmd/etat0dyn_netcdf.F90	(revision 3412)
+++ /LMDZ6/branches/DYNAMICO-conv/libf/dynphy_lonlat/phylmd/etat0dyn_netcdf.F90	(revision 3413)
@@ -100,5 +100,6 @@
 
   deg2rad = pi/180.0
-
+  y(:,:,:)=0  !ym warning unitialized variable
+  
 ! Compute psol AND tsol, knowing phis.
 !*******************************************************************************
Index: /LMDZ6/branches/DYNAMICO-conv/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90
===================================================================
--- /LMDZ6/branches/DYNAMICO-conv/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90	(revision 3412)
+++ /LMDZ6/branches/DYNAMICO-conv/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90	(revision 3413)
@@ -44,5 +44,5 @@
     zmax0,fevap, rnebcon,falb_dir, wake_fip,    agesno,  detr_therm, pbl_tke,  &
     phys_state_var_init, ql_ancien, qs_ancien, prlw_ancien, prsw_ancien, &
-    prw_ancien
+    prw_ancien, sollwdown
   USE comconst_mod, ONLY: pi, dtvr
 
@@ -200,4 +200,6 @@
   solsw      = 165.
   sollw      = -53.
+!ym warning missing init for sollwdown => set to 0
+  sollwdown  = 0.
   t_ancien   = 273.15
   q_ancien   = 0.
@@ -314,5 +316,8 @@
   ALLOCATE(zmea0(iml,jml),zstd0(iml,jml)) !--- Mean orography and std deviation
   ALLOCATE(zsig0(iml,jml),zgam0(iml,jml)) !--- Slope and nisotropy
+  zsig0(:,:)=0   !ym uninitialized variable
+  zgam0(:,:)=0   !ym uninitialized variable
   ALLOCATE(zthe0(iml,jml))                !--- Highest slope orientation
+  zthe0(:,:)=0   !ym uninitialized variable
   ALLOCATE(zpic0(iml,jml),zval0(iml,jml)) !--- Peaks and valley heights
 
Index: /LMDZ6/branches/DYNAMICO-conv/libf/misc/wxios.F90
===================================================================
--- /LMDZ6/branches/DYNAMICO-conv/libf/misc/wxios.F90	(revision 3412)
+++ /LMDZ6/branches/DYNAMICO-conv/libf/misc/wxios.F90	(revision 3413)
@@ -152,4 +152,5 @@
 !$OMP MASTER
         !Initialisation du contexte:
+        !!CALL xios_context_initialize(g_ctx_name, g_comm)
         CALL xios_context_initialize(g_ctx_name, COMM_LMDZ_PHY)
         CALL xios_get_handle(g_ctx_name, xios_ctx)    !Récupération
@@ -301,4 +302,5 @@
         CALL xios_set_domain_attr_hdl(dom, nj_glo=nbp_lat, jbegin=jj_begin-1, nj=jj_nb, data_dim=2)
         CALL xios_set_domain_attr_hdl(dom, lonvalue_1d=io_lon(1:nbp_lon), latvalue_1d=io_lat(jj_begin:jj_end))
+        CALL xios_set_domain_attr("dom_out", domain_ref=dom_id)
 
         IF (.NOT.is_sequential) THEN
@@ -332,4 +334,5 @@
         USE mod_phys_lmdz_para 
         USE nrtype, ONLY : PI
+        USE ioipsl_getin_p_mod, ONLY : getin_p
         IMPLICIT NONE
         CHARACTER(len=*),INTENT(IN) :: dom_id ! domain identifier
@@ -339,6 +342,6 @@
         REAL :: boundslat_mpi(klon_mpi,nvertex)
         INTEGER :: ind_cell_glo_mpi(klon_mpi)
-        TYPE(xios_domaingroup) :: dom
-
+        TYPE(xios_domain) :: dom
+        LOGICAL :: remap_output
 
         CALL gather_omp(longitude*180/PI,lon_mpi)
@@ -348,7 +351,9 @@
         CALL gather_omp(ind_cell_glo,ind_cell_glo_mpi)
         
+        remap_output=.TRUE.
+        CALL getin_p("remap_output",remap_output)
 
 !$OMP MASTER
-        CALL xios_get_domaingroup_handle(dom_id, dom)
+        CALL xios_get_domain_handle(dom_id, dom)
         
         !On parametrise le domaine:
@@ -357,4 +362,9 @@
                            bounds_lon_1d=TRANSPOSE(boundslon_mpi), bounds_lat_1d=TRANSPOSE(boundslat_mpi) )
         CALL xios_set_attr(dom, i_index=ind_cell_glo_mpi(:)-1)
+        IF (remap_output) THEN
+          CALL xios_set_domain_attr("dom_out", domain_ref="dom_regular")
+          CALL xios_set_fieldgroup_attr("dom_out", domain_ref="dom_regular")
+          CALL xios_set_fieldgroup_attr("remap_expr", expr="@this_ref")
+        ENDIF
 !$OMP END MASTER
 
Index: /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/conf_phys_m.F90
===================================================================
--- /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/conf_phys_m.F90	(revision 3412)
+++ /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/conf_phys_m.F90	(revision 3413)
@@ -1009,4 +1009,5 @@
     ! - 1 = stratospheric aerosols scaled from 550 nm AOD 
     ! - 2 = stratospheric aerosol properties from CMIP6
+    !Option 2 is only available with RRTM, this is tested later on 
     !Config Def  = 0
     !Config Help = Used in physiq.F
@@ -1014,7 +1015,5 @@
     !
     flag_aerosol_strat_omp = 0
-    IF (iflag_rrtm_omp==1) THEN
-      CALL getin('flag_aerosol_strat',flag_aerosol_strat_omp)
-    ENDIF
+    CALL getin('flag_aerosol_strat',flag_aerosol_strat_omp)
 
     !
Index: /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/cosp/phys_cosp.F90
===================================================================
--- /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/cosp/phys_cosp.F90	(revision 3412)
+++ /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/cosp/phys_cosp.F90	(revision 3413)
@@ -330,5 +330,4 @@
       !$OMP END MASTER
       !$OMP BARRIER
-        debut_cosp=.false.
       endif ! debut_cosp
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -341,5 +340,5 @@
 !        call cosp(overlap,Ncolumns,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,modis,rttov,stradar,stlidar)
 !#else
-        call cosp(overlap,Ncolumns,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,modis,stradar,stlidar)
+     if (.NOT. debut_cosp) call cosp(overlap,Ncolumns,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,modis,stradar,stlidar)
 !#endif
 !!
@@ -350,5 +349,5 @@
 
 !       print *, 'Calling write output'
-        call cosp_output_write(Nlevlmdz, Npoints, Ncolumns, itap, dtime, freq_COSP, missing_val, &
+     if (.NOT. debut_cosp) call cosp_output_write(Nlevlmdz, Npoints, Ncolumns, itap, dtime, freq_COSP, missing_val, &
                                cfg, gbx, vgrid, sglidar, sgradar, stlidar, stradar, & 
                                isccp, misr, modis)
@@ -376,4 +375,7 @@
 !  call system_clock(t1,count_rate,count_max)
 !  print *,(t1-t0)*1.0/count_rate
+    if (debut_cosp) then 
+      debut_cosp=.false.
+    endif
  
   CONTAINS 
Index: /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/create_etat0_unstruct.F90
===================================================================
--- /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/create_etat0_unstruct.F90	(revision 3412)
+++ /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/create_etat0_unstruct.F90	(revision 3413)
@@ -135,4 +135,8 @@
     rain_fall = 0.; snow_fall = 0.
     solsw = 165.;   sollw = -53.
+!ym warning missing init for sollwdown => set to 0
+  sollwdown  = 0.
+    
+    
     t_ancien = 273.15
     u_ancien=0
Index: /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/cv3_routines.F90
===================================================================
--- /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/cv3_routines.F90	(revision 3412)
+++ /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/cv3_routines.F90	(revision 3413)
@@ -3756,12 +3756,5 @@
         IF (i<=inb(il) .AND. nent(il,i)==0 .AND. iflag(il)<=1) THEN              ! cld
           qcond(il, i) = qcond(il, i) + (1.-ep(il,i))*clw(il, i)                 ! cld
-!jyg<   Bug correction 20180620
-!      PROBLEM: Should not qent(il,i,i) be taken into account even if nent(il,i)/=0?
-!!        qtment(il, i) = qent(il,k,i) + qtment(il,i)
-! cld
-          qtment(il, i) = qent(il,i,i) + qtment(il,i)
-! cld
-!>jyg
-
+          qtment(il, i) = qent(il,k,i) + qtment(il,i)                          ! cld
           nqcond(il, i) = nqcond(il, i) + 1.                                     ! cld
         END IF                                                                   ! cld
Index: /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/limit_read_mod.F90
===================================================================
--- /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/limit_read_mod.F90	(revision 3412)
+++ /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/limit_read_mod.F90	(revision 3413)
@@ -258,5 +258,5 @@
           ENDIF
           ierr=NF90_INQUIRE_DIMENSION(nid, ndimid, len=nn)
-          WRITE(abort_message,'(a,2(i3,a))')'limit.nc records number (',nn,') does no'//&
+          WRITE(abort_message,'(a,2(i0,a))')'limit.nc records number (',nn,') does no'//&
             't match year length (',year_len,')'
           IF(nn/=year_len) CALL abort_physic(modname,abort_message,1)
@@ -336,5 +336,5 @@
 
 !$OMP MASTER  ! Only master thread
-       IF (is_mpi_root) THEN ! Only master processus
+       IF (is_mpi_root) THEN ! Only master processus!
 
           ierr = NF90_OPEN ('limit.nc', NF90_NOWRITE, nid)
Index: /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/phys_cal_mod.F90
===================================================================
--- /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/phys_cal_mod.F90	(revision 3412)
+++ /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/phys_cal_mod.F90	(revision 3413)
@@ -36,8 +36,10 @@
   
   SUBROUTINE phys_cal_init(annee_ref,day_ref)
-  USE IOIPSL, ONLY:  ymds2ju, ioconf_calendar
-  USE mod_phys_lmdz_para, ONLY:  is_master,is_omp_master
-  USE ioipsl_getin_p_mod, ONLY: getin_p
-  IMPLICIT NONE
+
+    USE IOIPSL, ONLY:  ymds2ju, ioconf_calendar
+    USE mod_phys_lmdz_para, ONLY:  is_master,is_omp_master
+    USE ioipsl_getin_p_mod, ONLY: getin_p
+
+    IMPLICIT NONE
     INTEGER,INTENT(IN) :: annee_ref
     INTEGER,INTENT(IN) :: day_ref
@@ -46,5 +48,5 @@
     calend = 'earth_360d' ! default
     CALL getin_p("calend",calend)
-     
+
     IF (is_omp_master) THEN
       IF (calend == 'earth_360d') THEN
@@ -59,5 +61,5 @@
     ENDIF
 !$OMP BARRIER
-
+     
     CALL ymds2ju(annee_ref, 1, day_ref, 0., jD_ref)
     jD_ref=INT(jD_ref)
Index: /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/phys_state_var_mod.F90
===================================================================
--- /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/phys_state_var_mod.F90	(revision 3412)
+++ /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/phys_state_var_mod.F90	(revision 3413)
@@ -23,5 +23,5 @@
 !$OMP THREADPRIVATE(cvpas)
 !$OMP THREADPRIVATE(wkpas)
-      REAL, SAVE :: phys_tstep, solaire_etat0
+      REAL, SAVE :: phys_tstep=0, solaire_etat0
 !$OMP THREADPRIVATE(phys_tstep, solaire_etat0)
 
@@ -607,6 +607,8 @@
         du_gwd_rando(:,:)=0.
       ENDIF
-      if (.not. ok_hines .and. ok_gwd_rando) allocate(du_gwd_front(klon, klev))
-
+      IF (.not. ok_hines .and. ok_gwd_rando) THEN
+        ALLOCATE(du_gwd_front(klon, klev))
+        du_gwd_front(:,:) = 0 !ym missing init   
+      ENDIF
 END SUBROUTINE phys_state_var_init
 
Index: /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/physiq_mod.F90
===================================================================
--- /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/physiq_mod.F90	(revision 3412)
+++ /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/physiq_mod.F90	(revision 3413)
@@ -258,4 +258,7 @@
     USE climoz_mod
     USE limit_read_mod, ONLY : init_limit_read
+    USE regr_horiz_time_climoz_m, ONLY: regr_horiz_time_climoz
+    USE readaerosol_mod, ONLY : init_aero_fromfile
+    USE readaerosolstrato_m, ONLY : init_readaerosolstrato
 
 
@@ -1195,4 +1198,6 @@
        CALL phys_state_var_init(read_climoz)
        CALL phys_output_var_init
+       IF(read_climoz>=1 .AND. create_etat0_limit) CALL regr_horiz_time_climoz(read_climoz,ok_daily_climoz)
+
        print*, '================================================='
        !
@@ -1492,5 +1497,22 @@
        !
 
-       CALL create_climoz(read_climoz)
+!       CALL create_climoz(read_climoz)
+       CALL init_aero_fromfile(flag_aerosol)  !! initialise aero from file for XIOS interpolation (unstructured_grid)
+       CALL init_readaerosolstrato(flag_aerosol_strat)  !! initialise aero strato from file for XIOS interpolation (unstructured_grid)
+
+#ifdef CPP_COSP
+       CALL phys_cosp(itap,phys_tstep,freq_cosp, &
+               ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
+               ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
+               klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
+               JrNt,ref_liq,ref_ice, &
+               pctsrf(:,is_ter)+pctsrf(:,is_lic), &
+               zu10m,zv10m,pphis, &
+               zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
+               qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
+               prfl(:,1:klev),psfl(:,1:klev), &
+               pmflxr(:,1:klev),pmflxs(:,1:klev), &
+               mr_ozone,cldtau, cldemi)
+#endif
 
        CALL phys_output_write(itap, pdtphys, paprs, pphis,                    &
@@ -1504,6 +1526,6 @@
        IF (is_omp_master) CALL xios_update_calendar(1)
 #endif
+       IF(read_climoz>=1 .AND. create_etat0_limit) CALL regr_horiz_time_climoz(read_climoz,ok_daily_climoz)
        CALL create_etat0_limit_unstruct
-
        CALL phyetat0 ("startphy.nc",clesphy0,tabcntr0)
 
@@ -1535,7 +1557,4 @@
        ENDIF
 
-       CALL printflag( tabcntr0,radpas,ok_journe, &
-            ok_instan, ok_region )
-       !
 !       IF (ABS(phys_tstep-pdtphys).GT.0.001) THEN
 !          WRITE(lunout,*) 'Pas physique n est pas correct',phys_tstep, &
@@ -4733,13 +4752,14 @@
        !         write(97) u_seri,v_seri,t_seri,q_seri
        !         close(97)
-!       !$OMP MASTER
-       IF (read_climoz >= 1) THEN
-          IF (is_mpi_root) THEN
-             CALL nf95_close(ncid_climoz)
-          ENDIF
-          DEALLOCATE(press_edg_climoz) ! pointer
-          DEALLOCATE(press_cen_climoz) ! pointer
-       ENDIF
-!       !$OMP END MASTER
+      
+       IF (is_omp_master) THEN
+       
+         IF (read_climoz >= 1) THEN
+           IF (is_mpi_root) CALL nf95_close(ncid_climoz)
+            DEALLOCATE(press_edg_climoz) ! pointer
+            DEALLOCATE(press_cen_climoz) ! pointer
+         ENDIF
+       
+       ENDIF
 #ifdef CPP_XIOS
        IF (is_omp_master) CALL xios_context_finalize
Index: /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/radlwsw_m.F90
===================================================================
--- /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/radlwsw_m.F90	(revision 3412)
+++ /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/radlwsw_m.F90	(revision 3413)
@@ -385,4 +385,8 @@
   lldebug=.FALSE.
 
+  ztopsw_aero(:,:)  = 0. !ym missing init : warning : not initialized in SW_AEROAR4
+  ztopsw0_aero(:,:) = 0. !ym missing init : warning : not initialized in SW_AEROAR4
+  zsolsw_aero(:,:)  = 0. !ym missing init : warning : not initialized in SW_AEROAR4
+  zsolsw0_aero(:,:) = 0. !ym missing init : warning : not initialized in SW_AEROAR4
   
   !
@@ -597,8 +601,4 @@
       ztopswaiaero(i)=0.
       zsolswaiaero(i)=0.
-      ztopsw_aero(i,:)  = 0. !ym missing init : warning : not initialized in SW_AEROAR4
-      ztopsw0_aero(i,:) = 0. !ym missing init : warning : not initialized in SW_AEROAR4
-      zsolsw_aero(i,:)  = 0. !ym missing init : warning : not initialized in SW_AEROAR4
-      zsolsw0_aero(i,:) = 0. !ym missing init : warning : not initialized in SW_AEROAR4
       ENDDO
 !     print *,'Avant SW_LMDAR4: PSCT zrmu0 zfract',PSCT, zrmu0, zfract
Index: /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/readaerosol.F90
===================================================================
--- /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/readaerosol.F90	(revision 3412)
+++ /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/readaerosol.F90	(revision 3413)
@@ -4,4 +4,11 @@
 
   REAL, SAVE :: not_valid=-333.
+  
+  INTEGER, SAVE :: nbp_lon_src
+!$OMP THREADPRIVATE(nbp_lon_src)  
+  INTEGER, SAVE :: nbp_lat_src
+!$OMP THREADPRIVATE(nbp_lat_src)  
+  REAL, ALLOCATABLE, SAVE    :: psurf_interp(:,:)
+!$OMP THREADPRIVATE(psurf_interp)  
 
 CONTAINS
@@ -167,4 +174,74 @@
 
 
+  SUBROUTINE init_aero_fromfile(flag_aerosol)
+  USE netcdf
+  USE mod_phys_lmdz_para
+  USE mod_grid_phy_lmdz, ONLY: grid_type, unstructured
+  USE xios
+  IMPLICIT NONE
+  INTEGER, INTENT(IN) :: flag_aerosol
+  REAL,ALLOCATABLE :: lat_src(:)
+  REAL,ALLOCATABLE :: lon_src(:)
+  CHARACTER(LEN=*),PARAMETER :: file_aerosol='aerosols.nat.nc'
+  CHARACTER(LEN=*),PARAMETER :: file_so4='so4.nat.nc'
+  INTEGER :: klev_src
+  INTEGER :: ierr ,ncid, dimID, varid
+  REAL :: null_array(0)
+
+  IF (flag_aerosol>0 .AND. grid_type==unstructured) THEN
+  
+    IF (is_omp_root) THEN
+  
+      IF (is_mpi_root) THEN
+    
+        IF (nf90_open(TRIM(file_aerosol), NF90_NOWRITE, ncid) /= NF90_NOERR) THEN
+          CALL check_err( nf90_open(TRIM(file_so4), NF90_NOWRITE, ncid), "pb open "//trim(file_so4) )
+        ENDIF
+
+        ! Read and test longitudes
+        CALL check_err( nf90_inq_dimid(ncid, "lon", dimID),"pb inq dim lon") 
+        CALL check_err( nf90_inquire_dimension(ncid, dimID, len = nbp_lon_src),"pb inq dim lon")
+        CALL check_err( nf90_inq_varid(ncid, 'lon', varid),"pb inq lon" )
+        ALLOCATE(lon_src(nbp_lon_src))
+        CALL check_err( nf90_get_var(ncid, varid, lon_src(:)),"pb get lon" )
+
+        ! Read and test latitudes
+        CALL check_err( nf90_inq_dimid(ncid, "lat", dimID),"pb inq dim lat") 
+        CALL check_err( nf90_inquire_dimension(ncid, dimID, len = nbp_lat_src),"pb inq dim lat")
+        CALL check_err( nf90_inq_varid(ncid, 'lat', varid),"pb inq lat" )
+        ALLOCATE(lat_src(nbp_lat_src))
+        CALL check_err( nf90_get_var(ncid, varid, lat_src(:)),"pb get lat" )
+        IF (nf90_inq_dimid(ncid, 'lev', dimid) /= NF90_NOERR) THEN
+          IF (nf90_inq_dimid(ncid, 'presnivs', dimid)/= NF90_NOERR) THEN
+             CALL check_err(nf90_inq_dimid(ncid, 'PRESNIVS', dimid),'dimension lev,PRESNIVS or presnivs not in file')
+          ENDIF
+        ENDIF
+        CALL check_err( nf90_inquire_dimension(ncid, dimid, len = klev_src),"pb inq dim for PRESNIVS or lev" )
+        CALL check_err( nf90_close(ncid),"pb in close" )   
+      ENDIF
+
+      CALL bcast_mpi(nbp_lat_src)
+      CALL bcast_mpi(nbp_lon_src)
+      CALL bcast_mpi(klev_src)
+
+      IF (is_mpi_root ) THEN
+        CALL xios_set_domain_attr("domain_aerosol",nj_glo=nbp_lat_src, nj=nbp_lat_src, jbegin=0, latvalue_1d=lat_src)
+        CALL xios_set_domain_attr("domain_aerosol",ni_glo=nbp_lon_src, ni=nbp_lon_src, ibegin=0, lonvalue_1d=lon_src)
+      ELSE
+        CALL xios_set_domain_attr("domain_aerosol",nj_glo=nbp_lat_src, nj=0, jbegin=0, latvalue_1d=null_array )
+        CALL xios_set_domain_attr("domain_aerosol",ni_glo=nbp_lon_src, ni=0, ibegin=0, lonvalue_1d=null_array)
+      ENDIF
+      CALL xios_set_axis_attr("axis_aerosol",n_glo=klev_src)
+      CALL xios_set_fieldgroup_attr("aerosols", enabled=.TRUE.)
+  
+    ENDIF
+    
+  ENDIF    
+ 
+  END SUBROUTINE init_aero_fromfile
+
+
+    
+
   SUBROUTINE get_aero_fromfile(varname, cyr, filename, klev_src, pt_ap, pt_b, pt_year, psurf_out, load_out)
 !****************************************************************************************
@@ -187,9 +264,10 @@
     USE netcdf
     USE dimphy
-    USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, klon_glo, &
-                                 grid2Dto1D_glo
+    USE mod_grid_phy_lmdz, ONLY: nbp_lon_=>nbp_lon, nbp_lat_=>nbp_lat, klon_glo, &
+                                 grid2Dto1D_glo, grid_type, unstructured
     USE mod_phys_lmdz_para
     USE iophy, ONLY : io_lon, io_lat
     USE print_control_mod, ONLY: lunout
+    USE xios
 
     IMPLICIT NONE
@@ -205,6 +283,9 @@
     REAL, POINTER, DIMENSION(:)           :: pt_b         ! Pointer for describing the vertical levels      
     REAL, POINTER, DIMENSION(:,:,:)       :: pt_year      ! Pointer-variabale from file, 12 month, grid : klon,klev_src
+    REAL, POINTER, DIMENSION(:,:,:)       :: pt_year_mpi  ! Pointer-variabale from file, 12 month, grid : klon,klev_src
     REAL, DIMENSION(klon,12), INTENT(OUT) :: psurf_out    ! Surface pression for 12 months
+    REAL, DIMENSION(klon_mpi,12)          :: psurf_out_mpi    ! Surface pression for 12 months
     REAL, DIMENSION(klon,12), INTENT(OUT) :: load_out     ! Aerosol mass load in each column
+    REAL, DIMENSION(klon_mpi,12)          :: load_out_mpi     ! Aerosol mass load in each column
     INTEGER                               :: nbr_tsteps   ! number of month in file read
 
@@ -220,15 +301,39 @@
     REAL, ALLOCATABLE, DIMENSION(:)       :: varktmp
 
-    REAL, DIMENSION(nbp_lon,nbp_lat,12)   :: psurf_glo2D   ! Surface pression for 12 months on dynamics global grid
+    REAL,  ALLOCATABLE                    :: psurf_glo2D(:,:,:)   ! Surface pression for 12 months on dynamics global grid
     REAL, DIMENSION(klon_glo,12)          :: psurf_glo1D   ! -"- on physical global grid
-    REAL, DIMENSION(nbp_lon,nbp_lat,12)   :: load_glo2D    ! Load for 12 months on dynamics global grid
+    REAL,  ALLOCATABLE                    :: load_glo2D(:,:,:)    ! Load for 12 months on dynamics global grid
     REAL, DIMENSION(klon_glo,12)          :: load_glo1D    ! -"- on physical global grid
-    REAL, DIMENSION(nbp_lon,nbp_lat)      :: vartmp
-    REAL, DIMENSION(nbp_lon)              :: lon_src              ! longitudes in file
-    REAL, DIMENSION(nbp_lat)              :: lat_src, lat_src_inv ! latitudes in file
+    REAL, ALLOCATABLE, DIMENSION(:,:)     :: vartmp
+    REAL, ALLOCATABLE,DIMENSION(:)        :: lon_src              ! longitudes in file
+    REAL, ALLOCATABLE,DIMENSION(:)        :: lat_src, lat_src_inv ! latitudes in file
     LOGICAL                               :: new_file             ! true if new file format detected
     LOGICAL                               :: invert_lat           ! true if the field has to be inverted for latitudes
-
-
+    INTEGER                               :: nbp_lon, nbp_lat
+    LOGICAL,SAVE                          :: first=.TRUE.
+!$OMP THREADPRIVATE(first)
+    
+    IF (grid_type==unstructured) THEN
+      nbp_lon=nbp_lon_src
+      nbp_lat=nbp_lat_src
+    ELSE
+      nbp_lon=nbp_lon_
+      nbp_lat=nbp_lat_
+    ENDIF
+    
+    IF (is_mpi_root) THEN
+    
+      ALLOCATE(psurf_glo2D(nbp_lon,nbp_lat,12))
+      ALLOCATE(load_glo2D(nbp_lon,nbp_lat,12))
+      ALLOCATE(vartmp(nbp_lon,nbp_lat))
+      ALLOCATE(lon_src(nbp_lon))
+      ALLOCATE(lat_src(nbp_lat))
+      ALLOCATE(lat_src_inv(nbp_lat))
+    ELSE
+      ALLOCATE(varyear(0,0,0,0))
+      ALLOCATE(psurf_glo2D(0,0,0))
+      ALLOCATE(load_glo2D(0,0,0))
+    ENDIF
+            
     ! Deallocate pointers
     IF (ASSOCIATED(pt_ap)) DEALLOCATE(pt_ap)
@@ -245,41 +350,44 @@
        CALL check_err( nf90_open(TRIM(fname), NF90_NOWRITE, ncid), "pb open "//trim(fname) )
 
+
+       IF (grid_type/=unstructured) THEN
+
 ! Test for equal longitudes and latitudes in file and model
 !****************************************************************************************
        ! Read and test longitudes
-       CALL check_err( nf90_inq_varid(ncid, 'lon', varid),"pb inq lon" )
-       CALL check_err( nf90_get_var(ncid, varid, lon_src(:)),"pb get lon" )
+         CALL check_err( nf90_inq_varid(ncid, 'lon', varid),"pb inq lon" )
+         CALL check_err( nf90_get_var(ncid, varid, lon_src(:)),"pb get lon" )
        
-       IF (maxval(ABS(lon_src - io_lon)) > 0.001) THEN
-          WRITE(lunout,*) 'Problem in longitudes read from file : ',TRIM(fname)
-          WRITE(lunout,*) 'longitudes in file ', TRIM(fname),' : ', lon_src
-          WRITE(lunout,*) 'longitudes in model :', io_lon
+         IF (maxval(ABS(lon_src - io_lon)) > 0.001) THEN
+            WRITE(lunout,*) 'Problem in longitudes read from file : ',TRIM(fname)
+            WRITE(lunout,*) 'longitudes in file ', TRIM(fname),' : ', lon_src
+            WRITE(lunout,*) 'longitudes in model :', io_lon
           
-          CALL abort_physic('get_aero_fromfile', 'longitudes are not the same in file and model',1)
-       END IF
-
-       ! Read and test latitudes
-       CALL check_err( nf90_inq_varid(ncid, 'lat', varid),"pb inq lat" )
-       CALL check_err( nf90_get_var(ncid, varid, lat_src(:)),"pb get lat" )
-
-       ! Invert source latitudes
-       DO j = 1, nbp_lat
-          lat_src_inv(j) = lat_src(nbp_lat +1 -j)
-       END DO
-
-       IF (maxval(ABS(lat_src - io_lat)) < 0.001) THEN
-          ! Latitudes are the same
-          invert_lat=.FALSE.
-       ELSE IF (maxval(ABS(lat_src_inv - io_lat)) < 0.001) THEN
-          ! Inverted source latitudes correspond to model latitudes
-          WRITE(lunout,*) 'latitudes will be inverted for file : ',TRIM(fname)
-          invert_lat=.TRUE.
-       ELSE
-          WRITE(lunout,*) 'Problem in latitudes read from file : ',TRIM(fname)
-          WRITE(lunout,*) 'latitudes in file ', TRIM(fname),' : ', lat_src      
-          WRITE(lunout,*) 'latitudes in model :', io_lat
-          CALL abort_physic('get_aero_fromfile', 'latitudes do not correspond between file and model',1)
-       END IF
-
+            CALL abort_physic('get_aero_fromfile', 'longitudes are not the same in file and model',1)
+         END IF
+
+         ! Read and test latitudes
+         CALL check_err( nf90_inq_varid(ncid, 'lat', varid),"pb inq lat" )
+         CALL check_err( nf90_get_var(ncid, varid, lat_src(:)),"pb get lat" )
+
+         ! Invert source latitudes
+         DO j = 1, nbp_lat
+            lat_src_inv(j) = lat_src(nbp_lat +1 -j)
+         END DO
+
+         IF (maxval(ABS(lat_src - io_lat)) < 0.001) THEN
+            ! Latitudes are the same
+            invert_lat=.FALSE.
+         ELSE IF (maxval(ABS(lat_src_inv - io_lat)) < 0.001) THEN
+            ! Inverted source latitudes correspond to model latitudes
+            WRITE(lunout,*) 'latitudes will be inverted for file : ',TRIM(fname)
+            invert_lat=.TRUE.
+         ELSE
+            WRITE(lunout,*) 'Problem in latitudes read from file : ',TRIM(fname)
+            WRITE(lunout,*) 'latitudes in file ', TRIM(fname),' : ', lat_src      
+            WRITE(lunout,*) 'latitudes in model :', io_lat
+            CALL abort_physic('get_aero_fromfile', 'latitudes do not correspond between file and model',1)
+         END IF
+       ENDIF
 
 ! 2) Check if old or new file is avalabale.
@@ -487,60 +595,64 @@
        END IF
 
-!     - Invert latitudes if necessary
-       DO imth=1, 12
-          IF (invert_lat) THEN
-
-             ! Invert latitudes for the variable
-             varmth(:,:,:) = varyear(:,:,:,imth) ! use varmth temporarly
-             DO k=1,klev_src
-                DO j=1,nbp_lat
-                   DO i=1,nbp_lon
-                      varyear(i,j,k,imth) = varmth(i,nbp_lat+1-j,k)
-                   END DO
-                END DO
-             END DO
+
+       IF (grid_type/=unstructured) THEN
+  !     - Invert latitudes if necessary
+         DO imth=1, 12
+            IF (invert_lat) THEN
+
+               ! Invert latitudes for the variable
+               varmth(:,:,:) = varyear(:,:,:,imth) ! use varmth temporarly
+               DO k=1,klev_src
+                  DO j=1,nbp_lat
+                     DO i=1,nbp_lon
+                        varyear(i,j,k,imth) = varmth(i,nbp_lat+1-j,k)
+                     END DO
+                  END DO
+               END DO
              
-             ! Invert latitudes for surface pressure
-             vartmp(:,:) = psurf_glo2D(:,:,imth)
-             DO j=1,nbp_lat
-                DO i=1,nbp_lon
-                   psurf_glo2D(i,j,imth)= vartmp(i,nbp_lat+1-j)
-                END DO
-             END DO
+               ! Invert latitudes for surface pressure
+               vartmp(:,:) = psurf_glo2D(:,:,imth)
+               DO j=1,nbp_lat
+                  DO i=1,nbp_lon
+                     psurf_glo2D(i,j,imth)= vartmp(i,nbp_lat+1-j)
+                  END DO
+               END DO
              
-             ! Invert latitudes for the load
-             vartmp(:,:) = load_glo2D(:,:,imth)
-             DO j=1,nbp_lat
-                DO i=1,nbp_lon
-                   load_glo2D(i,j,imth)= vartmp(i,nbp_lat+1-j)
-                END DO
-             END DO
-          END IF ! invert_lat
+               ! Invert latitudes for the load
+               vartmp(:,:) = load_glo2D(:,:,imth)
+               DO j=1,nbp_lat
+                  DO i=1,nbp_lon
+                     load_glo2D(i,j,imth)= vartmp(i,nbp_lat+1-j)
+                  END DO
+               END DO
+            END IF ! invert_lat
              
-          ! Do zonal mead at poles and distribut at whole first and last latitude
-          DO k=1, klev_src
-             npole=0.  ! North pole, j=1
-             spole=0.  ! South pole, j=nbp_lat        
-             DO i=1,nbp_lon
-                npole = npole + varyear(i,1,k,imth)
-                spole = spole + varyear(i,nbp_lat,k,imth)
-             END DO
-             npole = npole/REAL(nbp_lon)
-             spole = spole/REAL(nbp_lon)
-             varyear(:,1,    k,imth) = npole
-             varyear(:,nbp_lat,k,imth) = spole
-          END DO
-       END DO ! imth
+            ! Do zonal mead at poles and distribut at whole first and last latitude
+            DO k=1, klev_src
+               npole=0.  ! North pole, j=1
+               spole=0.  ! South pole, j=nbp_lat        
+               DO i=1,nbp_lon
+                  npole = npole + varyear(i,1,k,imth)
+                  spole = spole + varyear(i,nbp_lat,k,imth)
+               END DO
+               npole = npole/REAL(nbp_lon)
+               spole = spole/REAL(nbp_lon)
+               varyear(:,1,    k,imth) = npole
+               varyear(:,nbp_lat,k,imth) = spole
+            END DO
+         END DO ! imth
        
-       ALLOCATE(varyear_glo1D(klon_glo, klev_src, 12), stat=ierr)
-       IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 3',1)
+         ALLOCATE(varyear_glo1D(klon_glo, klev_src, 12), stat=ierr)
+         IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 3',1)
        
-       ! Transform from 2D to 1D field
-       CALL grid2Dto1D_glo(varyear,varyear_glo1D)
-       CALL grid2Dto1D_glo(psurf_glo2D,psurf_glo1D)
-       CALL grid2Dto1D_glo(load_glo2D,load_glo1D)
-
+         ! Transform from 2D to 1D field
+         CALL grid2Dto1D_glo(varyear,varyear_glo1D)
+         CALL grid2Dto1D_glo(psurf_glo2D,psurf_glo1D)
+         CALL grid2Dto1D_glo(load_glo2D,load_glo1D)
+      
+      ENDIF
+      
     ELSE
-      ALLOCATE(varyear_glo1D(0,0,0))        
+        ALLOCATE(varyear_glo1D(0,0,0))        
     END IF ! is_mpi_root .AND. is_omp_root
 
@@ -566,11 +678,29 @@
     IF (ASSOCIATED(pt_year)) DEALLOCATE(pt_year)
     ALLOCATE(pt_year(klon, klev_src, 12), stat=ierr)
+    ALLOCATE(pt_year_mpi(klon_mpi, klev_src, 12), stat=ierr)
     IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 5',1)
 
-    ! Scatter global field to local domain at local process
-    CALL scatter(varyear_glo1D, pt_year)
-    CALL scatter(psurf_glo1D, psurf_out)
-    CALL scatter(load_glo1D,  load_out)
-
+    IF (grid_type==unstructured) THEN
+      IF (is_omp_master) THEN
+        CALL xios_send_field(TRIM(varname)//"_in",varyear)
+        CALL xios_recv_field(TRIM(varname)//"_out",pt_year_mpi)
+        CALL xios_send_field("load_"//TRIM(varname)//"_in",load_glo2D)
+        CALL xios_recv_field("load_"//TRIM(varname)//"_out",load_out_mpi)
+        IF (first) THEN
+          ALLOCATE(psurf_interp(klon_mpi,12))
+          CALL xios_send_field("psurf_aerosol_in",psurf_glo2D)
+          CALL xios_recv_field("psurf_aerosol_out",psurf_interp)
+        ENDIF
+      ENDIF
+      CALL scatter_omp(pt_year_mpi,pt_year)
+      CALL scatter_omp(load_out_mpi,load_out)
+      CALL scatter_omp(psurf_interp,psurf_out)
+      first=.FALSE.
+    ELSE
+      ! Scatter global field to local domain at local process
+      CALL scatter(varyear_glo1D, pt_year)
+      CALL scatter(psurf_glo1D, psurf_out)
+      CALL scatter(load_glo1D,  load_out)
+    ENDIF
 ! 7) Test for negative values
 !****************************************************************************************
Index: /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/readaerosolstrato.F90
===================================================================
--- /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/readaerosolstrato.F90	(revision 3412)
+++ /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/readaerosolstrato.F90	(revision 3413)
@@ -7,5 +7,5 @@
     USE phys_cal_mod, ONLY : mth_cur
     USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, klon_glo, &
-                                 grid2dto1d_glo
+                                 grid2dto1d_glo, grid_type, unstructured
     USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
     USE mod_phys_lmdz_omp_data, ONLY :  is_omp_root
@@ -15,4 +15,5 @@
     USE aero_mod
     USE dimphy
+    USE xios
 
     implicit none
@@ -43,4 +44,5 @@
     real, allocatable:: tauaerstrat_mois(:, :, :)
     real, allocatable:: tauaerstrat_mois_glo(:, :)
+    real, allocatable:: tau_aer_strat_mpi(:, :)
 
 ! For NetCDF:
@@ -87,18 +89,22 @@
     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
-
+    IF (grid_type/=unstructured) THEN
+      IF (n_lat.NE.nbp_lat) THEN 
+         print *,'Le nombre de lat n est pas egal a nbp_lat'
+         STOP
+      ENDIF
+    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
-
+    IF (grid_type/=unstructured) THEN
+      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
+    ENDIF
+    
     CALL nf95_inq_varid(ncid_in, "TIME", varid)
     CALL nf95_gw_var(ncid_in, varid, time)
@@ -130,10 +136,23 @@
     CALL grid2dTo1d_glo(tauaerstrat_mois,tauaerstrat_mois_glo)
 
+    ELSE
+      ALLOCATE(tauaerstrat_mois(0,0,0))
     ENDIF !--is_mpi_root and is_omp_root
 
 !$OMP BARRIER
 
+    IF (grid_type==unstructured) THEN
+      IF (is_omp_master) THEN
+        CALL xios_send_field("taustrat_in",tauaerstrat_mois)
+        ALLOCATE(tau_aer_strat_mpi(klon_mpi, klev))
+        CALL xios_recv_field("taustrat_out",tau_aer_strat_mpi)
+      ELSE
+        ALLOCATE(tau_aer_strat_mpi(0,0))
+      ENDIF
+      CALL scatter_omp(tau_aer_strat_mpi,tau_aer_strat)
+    ELSE  
 !--scatter on all proc
-    CALL scatter(tauaerstrat_mois_glo,tau_aer_strat)
+      CALL scatter(tauaerstrat_mois_glo,tau_aer_strat)
+    ENDIF
 
 !--keep memory of previous month
Index: /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/regr_horiz_time_climoz_m.F90
===================================================================
--- /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/regr_horiz_time_climoz_m.F90	(revision 3412)
+++ /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/regr_horiz_time_climoz_m.F90	(revision 3413)
@@ -2,5 +2,5 @@
 
   USE interpolation,     ONLY: locate
-  USE mod_grid_phy_lmdz, ONLY: nlon_ou => nbp_lon, nlat_ou => nbp_lat
+  USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, grid_type, unstructured
   USE nrtype,            ONLY: pi
   USE netcdf,   ONLY: NF90_CLOBBER, NF90_FLOAT,     NF90_GET_VAR, NF90_OPEN,   &
@@ -11,4 +11,5 @@
           NF95_CLOSE, NF95_ENDDEF,  NF95_PUT_ATT,   NF95_PUT_VAR, NF95_COPY_ATT
   USE print_control_mod, ONLY: lunout
+  USE dimphy
   IMPLICIT NONE
   PRIVATE
@@ -16,4 +17,10 @@
   REAL, PARAMETER :: deg2rad=pi/180.
   CHARACTER(LEN=13), PARAMETER :: vars_in(2)=['tro3         ','tro3_daylight']
+
+  INTEGER :: nlat_ou, nlon_ou
+  REAL, ALLOCATABLE :: latitude_glo(:)
+!$OMP THREADPRIVATE(latitude_glo)
+  INTEGER, ALLOCATABLE :: ind_cell_glo_glo(:)
+!$OMP THREADPRIVATE(ind_cell_glo_glo)
 
 CONTAINS
@@ -52,4 +59,5 @@
   USE assert_m,           ONLY: assert
   USE cal_tools_m,        ONLY: year_len, mid_month
+!!  USE control_mod,        ONLY: anneeref
   USE time_phylmdz_mod,   ONLY: annee_ref
   USE ioipsl,             ONLY: ioget_year_len, ioget_calendar
@@ -58,4 +66,9 @@
   USE regular_lonlat_mod, ONLY: boundslon_reg, boundslat_reg, south, west, east
   USE slopes_m,           ONLY: slopes
+  USE xios
+  USE mod_phys_lmdz_para, ONLY: is_mpi_root, is_master, is_omp_master, gather, gather_mpi, bcast_mpi, klon_mpi
+  USE geometry_mod, ONLY : latitude_deg, ind_cell_glo
+  USE mod_grid_phy_lmdz, ONLY: klon_glo
+
 !-------------------------------------------------------------------------------
 ! Arguments:
@@ -83,5 +96,7 @@
   CHARACTER(LEN=20) :: cal_in              ! Calendar
   REAL, ALLOCATABLE :: o3_in3(:,:,:,:,:)   ! Ozone climatologies
+  REAL, ALLOCATABLE :: o3_in3bis(:,:,:,:,:)   ! Ozone climatologies
   REAL, ALLOCATABLE :: o3_in2  (:,:,:,:)   ! Ozone climatologies
+  REAL, ALLOCATABLE :: o3_in2bis(:,:,:,:,:)   ! Ozone climatologies
   ! last index: 1 for the day-night average, 2 for the daylight field.
   REAL :: NaN
@@ -91,6 +106,9 @@
   REAL, ALLOCATABLE :: o3_regr_lonlat(:,:,:,:,:) ! (nlon_ou,nlat_ou,:,0:13   ,:)
   REAL, ALLOCATABLE :: o3_out3       (:,:,:,:,:) ! (nlon_ou,nlat_ou,:,ntim_ou,:)
+  REAL, ALLOCATABLE :: o3_out3_glo   (:,:,:,:) !   (nbp_lat,:,ntim_ou,:)
   REAL, ALLOCATABLE :: o3_regr_lat     (:,:,:,:) !         (nlat_in,:,0:13   ,:)
   REAL, ALLOCATABLE :: o3_out2         (:,:,:,:) !         (nlat_ou,:,ntim_ou,:)
+  REAL, ALLOCATABLE :: o3_out2_glo     (:,:,:,:) !         (nbp_lat,:,ntim_ou,:)
+  REAL, ALLOCATABLE :: o3_out          (:,:,:,:) !         (nbp_lat,:,ntim_ou,:)
 ! Dimension number  | Interval                | Contains  | For variables:
 !   1 (longitude)   | [rlonu(i-1), rlonu(i)]  | rlonv(i)  | all
@@ -116,282 +134,463 @@
   INTEGER, ALLOCATABLE :: sta(:), cnt(:)
   CHARACTER(LEN=80) :: sub, dim_nam, msg
-!-------------------------------------------------------------------------------
-  sub="regr_horiz_time_climoz"
-  WRITE(lunout,*)"Call sequence information: "//TRIM(sub)
-  CALL assert(read_climoz == 1 .OR. read_climoz == 2, "regr_lat_time_climoz")
-
-  CALL  NF95_OPEN("climoz.nc"  , NF90_NOWRITE, fID_in)
-  lprev=NF90_OPEN("climoz_m.nc", NF90_NOWRITE, fID_in_m)==NF90_NOERR
-  lnext=NF90_OPEN("climoz_p.nc", NF90_NOWRITE, fID_in_p)==NF90_NOERR
-
-  !--- Get coordinates from the input file. Converts lon/lat in radians.
-  !    Few inversions because "regr_conserv" and gcm need ascending vectors.
-  CALL NF95_INQ_VARID(fID_in, vars_in(1), varid)
-  CALL NF95_INQUIRE_VARIABLE(fID_in, varid, dimids=dIDs, ndims=ndims)
-  l3D=ndims==4; l2D=ndims==3
-  IF(l3D) WRITE(lunout,*)"Input files contain full 3D ozone fields."
-  IF(l2D) WRITE(lunout,*)"Input files contain zonal 2D ozone fields."
-  DO i=1,ndims
-    CALL NF95_INQUIRE_DIMENSION(fID_in, dIDs(i), name=dim_nam, nclen=dln)
-    CALL NF95_INQ_VARID(fID_in, dim_nam, varid)
-    ii=i; IF(l2D) ii=i+1                              !--- ndims==3:NO LONGITUDE 
-    SELECT CASE(ii)
-      CASE(1)                                         !--- LONGITUDE
-        CALL NF95_GW_VAR(fID_in, varid, lon_in)
-        ldec_lon=lon_in(1)>lon_in(dln); IF(ldec_lon) lon_in=lon_in(dln:1:-1)
-        nlon_in=dln; lon_in=lon_in*deg2rad
-      CASE(2)                                         !--- LATITUDE
-        CALL NF95_GW_VAR(fID_in, varid, lat_in)
-        ldec_lat=lat_in(1)>lat_in(dln); IF(ldec_lat) lat_in=lat_in(dln:1:-1)
-        nlat_in=dln; lat_in=lat_in*deg2rad
-      CASE(3)                                         !--- PRESSURE LEVELS
-        CALL NF95_GW_VAR(fID_in, varid, lev_in)
-        ldec_lev=lev_in(1)>lev_in(dln); IF(ldec_lev) lev_in=lev_in(dln:1:-1)
-        nlev_in=dln
-        CALL NF95_GET_ATT(fID_in, varid, "units", press_unit)
-        k=LEN_TRIM(press_unit)
-        DO WHILE(ICHAR(press_unit(k:k))==0)
-          press_unit(k:k)=' '; k=LEN_TRIM(press_unit) !--- REMOVE NULL END CHAR
-        END DO
-        IF(press_unit ==  "Pa") THEN
-          lev_in = lev_in/100.                        !--- CONVERT TO hPa
-        ELSE IF(press_unit /= "hPa") THEN
-          CALL abort_physic(sub, "the only recognized units are Pa and hPa.",1)
+  REAL :: null_array(0)
+  LOGICAL,SAVE :: first=.TRUE.
+!$OMP THREADPRIVATE(first)  
+  REAL, ALLOCATABLE :: test_o3_in(:,:)
+  REAL, ALLOCATABLE :: test_o3_out(:)
+
+
+  IF (grid_type==unstructured) THEN 
+    IF (first) THEN
+      IF (is_master) THEN
+        ALLOCATE(latitude_glo(klon_glo))
+        ALLOCATE(ind_cell_glo_glo(klon_glo))
+      ELSE
+        ALLOCATE(latitude_glo(0))
+        ALLOCATE(ind_cell_glo_glo(0))
+      ENDIF
+      CALL gather(latitude_deg,  latitude_glo)
+      CALL gather(ind_cell_glo,  ind_cell_glo_glo)
+    ENDIF 
+  ENDIF
+    
+  IF (is_omp_master) THEN
+    nlat_ou=nbp_lat
+    nlon_ou=nbp_lon
+    
+   !-------------------------------------------------------------------------------
+    IF (is_mpi_root) THEN
+      sub="regr_horiz_time_climoz"
+      WRITE(lunout,*)"Call sequence information: "//TRIM(sub)
+      CALL assert(read_climoz == 1 .OR. read_climoz == 2, "regr_lat_time_climoz")
+
+      CALL  NF95_OPEN("climoz.nc"  , NF90_NOWRITE, fID_in)
+      lprev=NF90_OPEN("climoz_m.nc", NF90_NOWRITE, fID_in_m)==NF90_NOERR
+      lnext=NF90_OPEN("climoz_p.nc", NF90_NOWRITE, fID_in_p)==NF90_NOERR
+
+      !--- Get coordinates from the input file. Converts lon/lat in radians.
+      !    Few inversions because "regr_conserv" and gcm need ascending vectors.
+      CALL NF95_INQ_VARID(fID_in, vars_in(1), varid)
+      CALL NF95_INQUIRE_VARIABLE(fID_in, varid, dimids=dIDs, ndims=ndims)
+      l3D=ndims==4; l2D=ndims==3
+      IF(l3D) WRITE(lunout,*)"Input files contain full 3D ozone fields."
+      IF(l2D) WRITE(lunout,*)"Input files contain zonal 2D ozone fields."
+      DO i=1,ndims
+        CALL NF95_INQUIRE_DIMENSION(fID_in, dIDs(i), name=dim_nam, nclen=dln)
+        CALL NF95_INQ_VARID(fID_in, dim_nam, varid)
+        ii=i; IF(l2D) ii=i+1                              !--- ndims==3:NO LONGITUDE 
+        SELECT CASE(ii)
+          CASE(1)                                         !--- LONGITUDE
+            CALL NF95_GW_VAR(fID_in, varid, lon_in)
+            ldec_lon=lon_in(1)>lon_in(dln); IF(ldec_lon) lon_in=lon_in(dln:1:-1)
+            nlon_in=dln; lon_in=lon_in*deg2rad
+          CASE(2)                                         !--- LATITUDE
+            CALL NF95_GW_VAR(fID_in, varid, lat_in)
+            ldec_lat=lat_in(1)>lat_in(dln); IF(ldec_lat) lat_in=lat_in(dln:1:-1)
+            nlat_in=dln; lat_in=lat_in*deg2rad
+          CASE(3)                                         !--- PRESSURE LEVELS
+            CALL NF95_GW_VAR(fID_in, varid, lev_in)
+            ldec_lev=lev_in(1)>lev_in(dln); IF(ldec_lev) lev_in=lev_in(dln:1:-1)
+            nlev_in=dln
+            CALL NF95_GET_ATT(fID_in, varid, "units", press_unit)
+            k=LEN_TRIM(press_unit)
+            DO WHILE(ICHAR(press_unit(k:k))==0)
+              press_unit(k:k)=' '; k=LEN_TRIM(press_unit) !--- REMOVE NULL END CHAR
+            END DO
+            IF(press_unit ==  "Pa") THEN
+              lev_in = lev_in/100.                        !--- CONVERT TO hPa
+            ELSE IF(press_unit /= "hPa") THEN
+              CALL abort_physic(sub, "the only recognized units are Pa and hPa.",1)
+            END IF
+          CASE(4)                                         !--- TIME
+            CALL NF95_INQUIRE_DIMENSION(fID_in, dIDs(i), nclen=nmth_in)
+            cal_in='gregorian'
+            IF(NF90_GET_ATT(fID_in, varid, 'calendar', cal_in)/=NF90_NOERR)        &  
+            WRITE(lunout,*)'WARNING: missing "calendar" attribute for "'//       &
+            TRIM(dim_nam)//'" in "climoz.nc". Choosing default: "gregorian".'
+            k=LEN_TRIM(cal_in)
+            DO WHILE(ICHAR(cal_in(k:k))==0)
+              cal_in(k:k)=' '; k=LEN_TRIM(cal_in)         !--- REMOVE NULL END CHAR
+            END DO
+        END SELECT
+      END DO
+
+      !--- Prepare quantities for time interpolation
+      tmidmonth=mid_month(annee_ref, cal_in)
+      IF(interpt) THEN
+        ntim_ou=ioget_year_len(annee_ref)
+        ALLOCATE(tmidday(ntim_ou))
+        tmidday=[(REAL(k)-0.5,k=1,ntim_ou)]
+        CALL ioget_calendar(cal_ou)
+      ELSE
+        ntim_ou=14
+        cal_ou=cal_in
+      END IF
+    ENDIF
+
+    IF (grid_type==unstructured) THEN
+      CALL bcast_mpi(nlon_in)
+      CALL bcast_mpi(nlat_in)
+      CALL bcast_mpi(nlev_in)
+      CALL bcast_mpi(l3d)
+      CALL bcast_mpi(tmidmonth)
+      CALL bcast_mpi(tmidday)
+      CALL bcast_mpi(ntim_ou)
+    
+      IF (is_mpi_root) THEN
+        CALL xios_set_domain_attr("domain_climoz",nj_glo=nlat_in, nj=nlat_in, jbegin=0, latvalue_1d=lat_in/deg2rad)
+        IF (l3D) THEN
+          CALL xios_set_domain_attr("domain_climoz",ni_glo=nlon_in, ni=nlon_in, ibegin=0, lonvalue_1d=lon_in/deg2rad)
+        ELSE
+          CALL xios_set_domain_attr("domain_climoz",ni_glo=8, ni=8, ibegin=0, lonvalue_1d = (/ 0.,45.,90.,135.,180.,225.,270., 315. /))
+        ENDIF
+      ELSE
+        CALL xios_set_domain_attr("domain_climoz",nj_glo=nlat_in, nj=0, jbegin=0, latvalue_1d=null_array )
+        IF (l3D) THEN
+          CALL xios_set_domain_attr("domain_climoz",ni_glo=nlon_in, ni=0, ibegin=0, lonvalue_1d=null_array)
+        ELSE
+          CALL xios_set_domain_attr("domain_climoz",ni_glo=8, ni=0, ibegin=0, lonvalue_1d=null_array)
+        ENDIF
+      ENDIF
+      CALL  xios_set_axis_attr("axis_climoz", n_glo=nlev_in)
+      CALL  xios_set_axis_attr("time_axis_climoz", n_glo=ntim_ou)
+      CALL  xios_set_axis_attr("time_axis_climoz", n_glo=ntim_ou)
+      CALL  xios_set_axis_attr("tr_climoz", n_glo=read_climoz)
+      CALL  xios_set_field_attr("tro3_out", enabled=.TRUE.)
+      CALL  xios_set_field_attr("tro3_out", enabled=.TRUE.)
+
+      IF (first) THEN
+        first=.FALSE.
+        RETURN
+      ENDIF
+    ENDIF
+    
+    
+    IF (is_mpi_root) THEN     
+      !--- Longitudes management:
+      !    * Need to shift data if the origin of input file longitudes /= -pi
+      !    * Need to add some margin in longitude to ensure input interval contains
+      !      all the output intervals => at least one longitudes slice has to be
+      !      duplicated, possibly more for undersampling.
+      IF(l3D) THEN
+        IF (grid_type==unstructured) THEN
+          dx2=0
+        ELSE
+          !--- Compute input edges longitudes vector (no end point yet)
+          ALLOCATE(v1(nlon_in+1))
+          v1(1)=(lon_in(nlon_in)+lon_in(1))/2.-pi
+          FORALL(i=2:nlon_in) v1(i)=(lon_in(i-1)+lon_in(i))/2.
+          v1(nlon_in+1)=v1(1)+2.*pi
+          DEALLOCATE(lon_in)
+
+          !--- Shift input longitudes vector until it contains first output point boundslon_reg(1,west)
+          v1=v1+2*pi*REAL(FLOOR((boundslon_reg(1,west)-v1(1))/(2.*pi)))
+
+          !--- Ensure first input longitudes interval contains first output point boundslon_reg(1,west)
+          dx1=locate(v1,boundslon_reg(1,west))-1
+          v1=CSHIFT(v1,SHIFT=dx1,DIM=1) 
+          v1(nlon_in-dx1+1:)=v1(nlon_in-dx1+1:)+2.*pi
+    
+          !--- Extend input longitudes vector until last interval contains boundslon_reg(nlat_ou,east)
+          dx2=0; DO WHILE(v1(1+dx2)+2.*pi<boundslon_reg(nlon_ou,east)); dx2=dx2+1; END DO
+
+          !--- Final edges longitudes vector (with margin and end point)
+          ALLOCATE(lon_in_edge(nlon_in+dx2+1)); lon_in_edge=[v1,v1(2:1+dx2)+2.*pi]
+          DEALLOCATE(v1)
+        ENDIF
+      END IF
+
+      !--- Compute sinus of intervals edges latitudes:
+      ALLOCATE(sinlat_in_edge(nlat_in+1))
+      sinlat_in_edge(1) = -1. ; sinlat_in_edge(nlat_in+1) = 1.
+      FORALL(j=2:nlat_in) sinlat_in_edge(j)=SIN((lat_in(j-1)+lat_in(j))/2.)
+      DEALLOCATE(lat_in)
+
+
+
+      !--- Check for contiguous years:
+      ib=0; ie=13
+      IF(nmth_in == 14) THEN; lprev=.FALSE.; lnext=.FALSE.
+        WRITE(lunout,*)'Using 14 months ozone climatology "climoz.nc"...'
+      ELSE  
+        IF(     lprev) WRITE(lunout,*)'Using "climoz_m.nc" last record (previous year).'
+        IF(.NOT.lprev) WRITE(lunout,*)"No previous year file ; assuming periodicity."
+        IF(     lnext) WRITE(lunout,*)'Using "climoz_p.nc" first record (next year).'
+        IF(.NOT.lnext) WRITE(lunout,*)"No next year file ; assuming periodicity."
+        IF(.NOT.lprev) ib=1
+        IF(.NOT.lnext) ie=12
+      END IF
+      ALLOCATE(sta(ndims),cnt(ndims)); sta(:)=1  
+      IF(l3D) cnt=[nlon_in,nlat_in,nlev_in,1]
+      IF(l2D) cnt=[        nlat_in,nlev_in,1]  
+      IF(l3D) ALLOCATE(o3_in3(nlon_in+dx2,nlat_in,nlev_in,ib:ie,read_climoz))
+      IF(l2D) ALLOCATE(o3_in2(            nlat_in,nlev_in,ib:ie,read_climoz))
+
+      !--- Read full current file and one record each available contiguous file
+      DO iv=1,read_climoz
+        msg=TRIM(sub)//" NF90_GET_VAR "//TRIM(vars_in(iv))
+        CALL NF95_INQ_VARID(fID_in, vars_in(1), vID_in(iv))
+        IF(l3D) ncerr=NF90_GET_VAR(fID_in, vID_in(iv), o3_in3(1:nlon_in,:,:,1:12,iv))
+        IF(l2D) ncerr=NF90_GET_VAR(fID_in, vID_in(iv), o3_in2(          :,:,1:12,iv))
+        CALL handle_err(TRIM(msg), ncerr, fID_in)
+        IF(lprev) THEN; sta(ndims)=12  
+          CALL NF95_INQ_VARID(fID_in_m, vars_in(1), vID_in(iv))
+          IF(l3D) ncerr=NF90_GET_VAR(fID_in_m,vID_in(iv),o3_in3(1:nlon_in,:,:, 0,iv),sta,cnt)
+          IF(l2d) ncerr=NF90_GET_VAR(fID_in_m,vID_in(iv),o3_in2(          :,:, 0,iv),sta,cnt)
+          CALL handle_err(TRIM(msg)//" previous", ncerr, fID_in_m)
         END IF
-      CASE(4)                                         !--- TIME
-        CALL NF95_INQUIRE_DIMENSION(fID_in, dIDs(i), nclen=nmth_in)
-        cal_in='gregorian'
-        IF(NF90_GET_ATT(fID_in, varid, 'calendar', cal_in)/=NF90_NOERR)        &
-          WRITE(lunout,*)'WARNING: missing "calendar" attribute for "'//       &
-          TRIM(dim_nam)//'" in "climoz.nc". Choosing default: "gregorian".'
-        k=LEN_TRIM(cal_in)
-        DO WHILE(ICHAR(cal_in(k:k))==0)
-          cal_in(k:k)=' '; k=LEN_TRIM(cal_in)         !--- REMOVE NULL END CHAR
-        END DO
-    END SELECT
-  END DO
-
-  !--- Longitudes management:
-  !    * Need to shift data if the origin of input file longitudes /= -pi
-  !    * Need to add some margin in longitude to ensure input interval contains
-  !      all the output intervals => at least one longitudes slice has to be
-  !      duplicated, possibly more for undersampling.
-  IF(l3D) THEN
-    !--- Compute input edges longitudes vector (no end point yet)
-    ALLOCATE(v1(nlon_in+1))
-    v1(1)=(lon_in(nlon_in)+lon_in(1))/2.-pi
-    FORALL(i=2:nlon_in) v1(i)=(lon_in(i-1)+lon_in(i))/2.
-    v1(nlon_in+1)=v1(1)+2.*pi
-    DEALLOCATE(lon_in)
-
-    !--- Shift input longitudes vector until it contains first output point boundslon_reg(1,west)
-    v1=v1+2*pi*REAL(FLOOR((boundslon_reg(1,west)-v1(1))/(2.*pi)))
-
-    !--- Ensure first input longitudes interval contains first output point boundslon_reg(1,west)
-    dx1=locate(v1,boundslon_reg(1,west))-1
-    v1=CSHIFT(v1,SHIFT=dx1,DIM=1); v1(nlon_in-dx1+1:)=v1(nlon_in-dx1+1:)+2.*pi
-
-    !--- Extend input longitudes vector until last interval contains boundslon_reg(nlat_ou,east)
-    dx2=0; DO WHILE(v1(1+dx2)+2.*pi<boundslon_reg(nlon_ou,east)); dx2=dx2+1; END DO
-
-    !--- Final edges longitudes vector (with margin and end point)
-    ALLOCATE(lon_in_edge(nlon_in+dx2+1)); lon_in_edge=[v1,v1(2:1+dx2)+2.*pi]
-    DEALLOCATE(v1)
-  END IF
-
-  !--- Compute sinus of intervals edges latitudes:
-  ALLOCATE(sinlat_in_edge(nlat_in+1))
-  sinlat_in_edge(1) = -1. ; sinlat_in_edge(nlat_in+1) = 1.
-  FORALL(j=2:nlat_in) sinlat_in_edge(j)=SIN((lat_in(j-1)+lat_in(j))/2.)
-  DEALLOCATE(lat_in)
-
-  !--- Prepare quantities for time interpolation
-  tmidmonth=mid_month(annee_ref, cal_in)
-  IF(interpt) THEN
-    ntim_ou=ioget_year_len(annee_ref)
-    ALLOCATE(tmidday(ntim_ou))
-    tmidday=[(REAL(k)-0.5,k=1,ntim_ou)]
-    CALL ioget_calendar(cal_ou)
-  ELSE
-    ntim_ou=14
-    cal_ou=cal_in
-  END IF
-
-  !--- Create the output file and get the variable IDs:
-  CALL prepare_out(fID_in,nlev_in,ntim_ou, fID_ou,levID_ou,timID_ou,vID_ou, &
-                   ndims, cal_ou)
-
-  !--- Write remaining coordinate variables:
-  CALL NF95_PUT_VAR(fID_ou, levID_ou, lev_in); DEALLOCATE(lev_in)
-  IF(     interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidday)
-  IF(.NOT.interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidmonth)
-
-  !--- Check for contiguous years:
-  ib=0; ie=13
-  IF(nmth_in == 14) THEN; lprev=.FALSE.; lnext=.FALSE.
-    WRITE(lunout,*)'Using 14 months ozone climatology "climoz.nc"...'
-  ELSE
-    IF(     lprev) WRITE(lunout,*)'Using "climoz_m.nc" last record (previous year).'
-    IF(.NOT.lprev) WRITE(lunout,*)"No previous year file ; assuming periodicity."
-    IF(     lnext) WRITE(lunout,*)'Using "climoz_p.nc" first record (next year).'
-    IF(.NOT.lnext) WRITE(lunout,*)"No next year file ; assuming periodicity."
-    IF(.NOT.lprev) ib=1
-    IF(.NOT.lnext) ie=12
-  END IF
-  ALLOCATE(sta(ndims),cnt(ndims)); sta(:)=1
-  IF(l3D) cnt=[nlon_in,nlat_in,nlev_in,1]
-  IF(l2D) cnt=[        nlat_in,nlev_in,1]
-  IF(l3D) ALLOCATE(o3_in3(nlon_in+dx2,nlat_in,nlev_in,ib:ie,read_climoz))
-  IF(l2D) ALLOCATE(o3_in2(            nlat_in,nlev_in,ib:ie,read_climoz))
-
-  !--- Read full current file and one record each available contiguous file
-  DO iv=1,read_climoz
-    msg=TRIM(sub)//" NF90_GET_VAR "//TRIM(vars_in(iv))
-    CALL NF95_INQ_VARID(fID_in, vars_in(1), vID_in(iv))
-    IF(l3D) ncerr=NF90_GET_VAR(fID_in, vID_in(iv), o3_in3(1:nlon_in,:,:,1:12,iv))
-    IF(l2D) ncerr=NF90_GET_VAR(fID_in, vID_in(iv), o3_in2(          :,:,1:12,iv))
-    CALL handle_err(TRIM(msg), ncerr, fID_in)
-    IF(lprev) THEN; sta(ndims)=12
-      CALL NF95_INQ_VARID(fID_in_m, vars_in(1), vID_in(iv))
-      IF(l3D) ncerr=NF90_GET_VAR(fID_in_m,vID_in(iv),o3_in3(1:nlon_in,:,:, 0,iv),sta,cnt)
-      IF(l2d) ncerr=NF90_GET_VAR(fID_in_m,vID_in(iv),o3_in2(          :,:, 0,iv),sta,cnt)
-      CALL handle_err(TRIM(msg)//" previous", ncerr, fID_in_m)
+        IF(lnext) THEN; sta(ndims)=1  
+          CALL NF95_INQ_VARID(fID_in_p, vars_in(1), vID_in(iv))
+          IF(l3D) ncerr=NF90_GET_VAR(fID_in_p,vID_in(iv),o3_in3(1:nlon_in,:,:,13,iv),sta,cnt)
+          IF(l2D) ncerr=NF90_GET_VAR(fID_in_p,vID_in(iv),o3_in2(          :,:,13,iv),sta,cnt)
+          CALL handle_err(TRIM(msg)//" next", ncerr, fID_in_p)
+        END IF
+      END DO
+      IF(lprev.OR.lnext) DEALLOCATE(sta,cnt)
+      IF(lprev) CALL NF95_CLOSE(fID_in_m)
+      IF(lnext) CALL NF95_CLOSE(fID_in_p)
+
+      !--- Revert decreasing coordinates vector
+      IF(l3D) THEN
+        IF(ldec_lon) o3_in3(1:nlon_in,:,:,:,:) = o3_in3(nlon_in:1:-1,:,:,:,:)
+        IF(ldec_lat) o3_in3 = o3_in3(:,nlat_in:1:-1,:,:,:)
+        IF(ldec_lev) o3_in3 = o3_in3(:,:,nlev_in:1:-1,:,:)
+        
+        IF (grid_type /= unstructured) THEN
+          !--- Shift values for longitude and duplicate some longitudes slices
+          o3_in3(1:nlon_in,:,:,:,:)=CSHIFT(o3_in3(1:nlon_in,:,:,:,:),SHIFT=dx1,DIM=1)
+          o3_in3(nlon_in+1:nlon_in+dx2,:,:,:,:)=o3_in3(1:dx2,:,:,:,:)
+        ENDIF
+      ELSE
+        IF(ldec_lat) o3_in2 = o3_in2(  nlat_in:1:-1,:,:,:)
+        IF(ldec_lev) o3_in2 = o3_in2(  :,nlev_in:1:-1,:,:)
+      END IF
+
+     !--- Deal with missing values 
+      DO m=1, read_climoz
+        WRITE(msg,'(a,i0)')"regr_lat_time_climoz: field Nr.",m
+        IF(NF90_GET_ATT(fID_in,vID_in(m),"missing_value",NaN)/= NF90_NOERR) THEN
+          IF(NF90_GET_ATT(fID_in, vID_in(m),"_FillValue",NaN)/= NF90_NOERR) THEN
+            WRITE(lunout,*)TRIM(msg)//": no missing value attribute found."; CYCLE
+          END IF
+        END IF
+        WRITE(lunout,*)TRIM(msg)//": missing value attribute found."
+        WRITE(lunout,*)"Trying to fill in NaNs ; a full field would be better."
+
+        !--- Check top layer contains no NaNs & search NaNs from top to ground
+        msg=TRIM(sub)//": NaNs in top layer !"
+        IF(l3D) THEN
+          IF(ANY(o3_in3(:,:,1,:,m)==NaN)) CALL abort_physic(sub,msg,1)
+          DO k = 2,nlev_in
+            WHERE(o3_in3(:,:,k,:,m)==NaN) o3_in3(:,:,k,:,m)=o3_in3(:,:,k-1,:,m)
+          END DO
+        ELSE
+          IF(ANY(o3_in2(  :,1,:,m)==NaN)) THEN
+            WRITE(lunout,*)msg
+            !--- Fill in latitudes where all values are missing
+            DO l=1,nmth_in
+              !--- Next to south pole
+              j=1;       DO WHILE(o3_in2(j,1,l,m)==NaN); j=j+1; END DO
+              IF(j>1) &
+                o3_in2(:j-1,:,l,m)=SPREAD(o3_in2(j,:,l,m),DIM=1,ncopies=j-1)
+              !--- Next to north pole
+              j=nlat_in; DO WHILE(o3_in2(j,1,l,m)==NaN); j=j+1; END DO
+              IF(j<nlat_in) &
+                o3_in2(j+1:,:,l,m)=SPREAD(o3_in2(j,:,l,m),DIM=1,ncopies=nlat_in-j)
+            END DO
+          END IF
+
+          !--- Fill in high latitudes missing values
+          !--- Highest level been filled-in, so has always valid values.
+          DO k = 2,nlev_in
+            WHERE(o3_in2(:,k,:,m)==NaN) o3_in2(:,k,:,m)=o3_in2(:,k-1,:,m)
+          END DO
+        END IF
+      END DO
+
+    ENDIF
+    
+    !=============================================================================
+    IF(l3D) THEN                                                   !=== 3D FIELDS
+    !=============================================================================
+     IF (grid_type==unstructured) THEN
+       nlat_ou=klon_mpi
+       
+       IF (is_mpi_root) THEN 
+         ALLOCATE(o3_in3bis(nlon_in,nlat_in,nlev_in,0:13,read_climoz)) 
+         o3_in3bis(:,:,:,ib:ie,:)=o3_in3(1:nlon_in,:,:,ib:ie,:)
+       ELSE
+         ALLOCATE(o3_in3bis(0,0,0,0,read_climoz))
+       ENDIF
+       ALLOCATE(o3_regr_lonlat(1, nlat_ou, nlev_in, 0:13, read_climoz))
+       
+       CALL xios_send_field("tro3_in",o3_in3bis(:,:,:,:,:))
+       CALL xios_recv_field("tro3_out",o3_regr_lonlat(1,:,:,:,:))
+
+     ELSE
+          
+       !--- Regrid in longitude
+        ALLOCATE(o3_regr_lon(nlon_ou, nlat_in, nlev_in, ie-ib+1, read_climoz))
+        CALL regr_conserv(1, o3_in3, xs = lon_in_edge,                             &
+                            xt = [boundslon_reg(1,west),boundslon_reg(:,east)],    &
+                            vt = o3_regr_lon, slope = slopes(1,o3_in3, lon_in_edge))
+        DEALLOCATE(o3_in3)
+
+        !--- Regrid in latitude: averaging with respect to SIN(lat) is
+        !                        equivalent to weighting by COS(lat)
+        !--- (inverted indices in "o3_regr_lonlat" because "rlatu" is decreasing)
+        ALLOCATE(o3_regr_lonlat(nlon_ou, nlat_ou, nlev_in, 0:13, read_climoz))
+        CALL regr_conserv(2, o3_regr_lon, xs = sinlat_in_edge,                     &
+                        xt = [- 1., SIN(boundslat_reg(nlat_ou-1:1:-1,south)), 1.], &
+                        vt = o3_regr_lonlat(:,nlat_ou:1:- 1,:,ib:ie,:),            &
+                   slope = slopes(2,o3_regr_lon, sinlat_in_edge))
+        DEALLOCATE(o3_regr_lon)
+
+     ENDIF
+
+     !--- Duplicate previous/next record(s) if they are not available
+     IF(.NOT.lprev) o3_regr_lonlat(:,:,:, 0,:) = o3_regr_lonlat(:,:,:,12,:)
+     IF(.NOT.lnext) o3_regr_lonlat(:,:,:,13,:) = o3_regr_lonlat(:,:,:, 1,:)
+      
+     !--- Regrid in time by linear interpolation:
+     ALLOCATE(o3_out3(nlon_ou, nlat_ou, nlev_in, ntim_ou, read_climoz))
+     IF(     interpt) CALL regr_lint(4,o3_regr_lonlat,tmidmonth,tmidday,o3_out3)
+     IF(.NOT.interpt) o3_out3=o3_regr_lonlat
+     DEALLOCATE(o3_regr_lonlat)
+
+     nlat_ou=nbp_lat
+     IF (grid_type==unstructured) THEN
+       CALL xios_send_field('o3_out',o3_out3)
+       ndims=3 
+       ALLOCATE(o3_out3_glo(nlat_ou, nlev_in, ntim_ou, read_climoz))
+       CALL gather_mpi(o3_out3(1,:,:,:,:), o3_out3_glo)
+     ENDIF
+
+    !--- Create the output file and get the variable IDs:
+    CALL prepare_out(fID_in,nlev_in,ntim_ou, fID_ou,levID_ou,timID_ou,vID_ou, &
+                     ndims, cal_ou)
+
+    IF (is_mpi_root) THEN
+      !--- Write remaining coordinate variables:
+      CALL NF95_PUT_VAR(fID_ou, levID_ou, lev_in); DEALLOCATE(lev_in)
+      IF(     interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidday)
+      IF(.NOT.interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidmonth)
+
+      !--- Write to file (the order of "rlatu" is inverted in the output file):
+        IF (grid_type==unstructured) THEN
+
+          ALLOCATE(o3_out(nlat_ou, nlev_in, ntim_ou, read_climoz))
+          DO i=1,klon_glo
+            o3_out(ind_cell_glo_glo(i),:,:,:)=o3_out3_glo(i,:,:,:)
+          ENDDO
+
+          DO m = 1, read_climoz
+            CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out(nlat_ou:1:-1,:,:,m))
+          END DO
+          
+        ELSE
+          DO m = 1, read_climoz
+            CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out3(:,nlat_ou:1:-1,:,:,m))
+          END DO
+      ENDIF
+      CALL NF95_CLOSE(fID_ou)
+
+
+    ENDIF
+
+
+    !=============================================================================
+    ELSE                                                         !=== ZONAL FIELDS
+    !=============================================================================
+   
+     IF (grid_type==unstructured) THEN
+       nlat_ou=klon_mpi
+
+       IF (is_mpi_root) THEN 
+         ALLOCATE(o3_in2bis(8,nlat_in,nlev_in,0:13,read_climoz)) 
+         o3_in2bis(:,:,:,ib:ie,:)=SPREAD(o3_in2,1,8)
+       ELSE
+         ALLOCATE(o3_in2bis(0,0,0,0,read_climoz))
+       ENDIF
+       ALLOCATE(o3_regr_lat(nlat_ou, nlev_in, 0:13, read_climoz))
+       CALL xios_send_field("tro3_in",o3_in2bis(:,:,:,:,:))
+       CALL xios_recv_field("tro3_out",o3_regr_lat(:,:,:,:))
+       IF(.NOT.lprev) o3_regr_lat(:,:, 0, :) = o3_regr_lat(:,:,12,:)
+       IF(.NOT.lnext) o3_regr_lat(:,:,13, :) = o3_regr_lat(:,:, 1,:)
+       
+     
+     ELSE 
+        !--- Regrid in latitude: averaging with respect to SIN(lat) is
+        !                        equivalent to weighting by COS(lat)
+        !--- (inverted indices in "o3_regr_lat" because "rlatu" is decreasing)
+        ALLOCATE(o3_regr_lat(nlat_ou, nlev_in, 0:13, read_climoz))
+        CALL regr_conserv(1, o3_in2, xs = sinlat_in_edge,                          &
+                        xt = [- 1., SIN(boundslat_reg(nlat_ou-1:1:-1,south)), 1.], &
+                        vt = o3_regr_lat(nlat_ou:1:- 1,:,ib:ie,:),                 &
+                     slope = slopes(1,o3_in2, sinlat_in_edge))
+        DEALLOCATE(o3_in2)
+
+        !--- Duplicate previous/next record(s) if they are not available
+        IF(.NOT.lprev) o3_regr_lat(:,:, 0,:) = o3_regr_lat(:,:,12,:)
+        IF(.NOT.lnext) o3_regr_lat(:,:,13,:) = o3_regr_lat(:,:, 1,:)
+
+     ENDIF
+     
+      !--- Regrid in time by linear interpolation:
+      ALLOCATE(o3_out2(nlat_ou, nlev_in, ntim_ou, read_climoz))
+      IF(     interpt) CALL regr_lint(3,o3_regr_lat, tmidmonth, tmidday, o3_out2)
+      IF(.NOT.interpt) o3_out2=o3_regr_lat
+      DEALLOCATE(o3_regr_lat)
+
+      nlat_ou=nbp_lat
+    
+      IF (grid_type==unstructured) THEN
+        ndims=3 
+        ALLOCATE(o3_out2_glo(nlat_ou, nlev_in, ntim_ou, read_climoz))
+        CALL gather_mpi(o3_out2, o3_out2_glo)
+      ENDIF
+      
+      !--- Create the output file and get the variable IDs:
+      CALL prepare_out(fID_in,nlev_in,ntim_ou, fID_ou,levID_ou,timID_ou,vID_ou, &
+                         ndims, cal_ou)
+
+      IF (is_mpi_root) THEN
+      
+        !--- Write remaining coordinate variables:
+        CALL NF95_PUT_VAR(fID_ou, levID_ou, lev_in); DEALLOCATE(lev_in)
+        IF(     interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidday)
+        IF(.NOT.interpt) CALL NF95_PUT_VAR(fID_ou, timID_ou, tmidmonth)
+
+        IF (grid_type==unstructured) THEN
+
+          ALLOCATE(o3_out3_glo(nlat_ou, nlev_in, ntim_ou, read_climoz))
+          DO i=1,klon_glo
+            o3_out(ind_cell_glo_glo(i),:,:,:)=o3_out2_glo(i,:,:,:)
+          ENDDO
+
+
+          DO m = 1, read_climoz
+            CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out(nlat_ou:1:-1,:,:,m))
+          END DO
+        ELSE
+          !--- Write to file (the order of "rlatu" is inverted in the output file):
+          DO m = 1, read_climoz
+            CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out2(nlat_ou:1:-1,:,:,m))
+          END DO
+        ENDIF
+        
+        CALL NF95_CLOSE(fID_ou)
+      
+      ENDIF
+
+    !=============================================================================
     END IF
-    IF(lnext) THEN; sta(ndims)=1
-      CALL NF95_INQ_VARID(fID_in_p, vars_in(1), vID_in(iv))
-      IF(l3D) ncerr=NF90_GET_VAR(fID_in_p,vID_in(iv),o3_in3(1:nlon_in,:,:,13,iv),sta,cnt)
-      IF(l2D) ncerr=NF90_GET_VAR(fID_in_p,vID_in(iv),o3_in2(          :,:,13,iv),sta,cnt)
-      CALL handle_err(TRIM(msg)//" next", ncerr, fID_in_p)
-    END IF
-  END DO
-  IF(lprev.OR.lnext) DEALLOCATE(sta,cnt)
-  IF(lprev) CALL NF95_CLOSE(fID_in_m)
-  IF(lnext) CALL NF95_CLOSE(fID_in_p)
-
-  !--- Revert decreasing coordinates vector
-  IF(l3D) THEN
-    IF(ldec_lon) o3_in3(1:nlon_in,:,:,:,:) = o3_in3(nlon_in:1:-1,:,:,:,:)
-    IF(ldec_lat) o3_in3 = o3_in3(:,nlat_in:1:-1,:,:,:)
-    IF(ldec_lev) o3_in3 = o3_in3(:,:,nlev_in:1:-1,:,:)
-    !--- Shift values for longitude and duplicate some longitudes slices
-    o3_in3(1:nlon_in,:,:,:,:)=CSHIFT(o3_in3(1:nlon_in,:,:,:,:),SHIFT=dx1,DIM=1)
-    o3_in3(nlon_in+1:nlon_in+dx2,:,:,:,:)=o3_in3(1:dx2,:,:,:,:)
-  ELSE
-    IF(ldec_lat) o3_in2 = o3_in2(  nlat_in:1:-1,:,:,:)
-    IF(ldec_lev) o3_in2 = o3_in2(  :,nlev_in:1:-1,:,:)
-  END IF
-
- !--- Deal with missing values
-  DO m=1, read_climoz
-    WRITE(msg,'(a,i0)')"regr_lat_time_climoz: field Nr.",m
-    IF(NF90_GET_ATT(fID_in,vID_in(m),"missing_value",NaN)/= NF90_NOERR) THEN
-      IF(NF90_GET_ATT(fID_in, vID_in(m),"_FillValue",NaN)/= NF90_NOERR) THEN
-        WRITE(lunout,*)TRIM(msg)//": no missing value attribute found."; CYCLE
-      END IF
-    END IF
-    WRITE(lunout,*)TRIM(msg)//": missing value attribute found."
-    WRITE(lunout,*)"Trying to fill in NaNs ; a full field would be better."
-
-    !--- Check top layer contains no NaNs & search NaNs from top to ground
-    msg=TRIM(sub)//": NaNs in top layer !"
-    IF(l3D) THEN
-      IF(ANY(o3_in3(:,:,1,:,m)==NaN)) CALL abort_physic(sub,msg,1)
-      DO k = 2,nlev_in
-        WHERE(o3_in3(:,:,k,:,m)==NaN) o3_in3(:,:,k,:,m)=o3_in3(:,:,k-1,:,m)
-      END DO
-    ELSE
-      IF(ANY(o3_in2(  :,1,:,m)==NaN)) THEN
-        WRITE(lunout,*)msg
-        !--- Fill in latitudes where all values are missing
-        DO l=1,nmth_in
-          !--- Next to south pole
-          j=1;       DO WHILE(o3_in2(j,1,l,m)==NaN); j=j+1; END DO
-          IF(j>1) &
-            o3_in2(:j-1,:,l,m)=SPREAD(o3_in2(j,:,l,m),DIM=1,ncopies=j-1)
-          !--- Next to north pole
-          j=nlat_in; DO WHILE(o3_in2(j,1,l,m)==NaN); j=j+1; END DO
-          IF(j<nlat_in) &
-            o3_in2(j+1:,:,l,m)=SPREAD(o3_in2(j,:,l,m),DIM=1,ncopies=nlat_in-j)
-        END DO
-      END IF
-
-      !--- Fill in high latitudes missing values
-      !--- Highest level been filled-in, so has always valid values.
-      DO k = 2,nlev_in
-        WHERE(o3_in2(:,k,:,m)==NaN) o3_in2(:,k,:,m)=o3_in2(:,k-1,:,m)
-      END DO
-    END IF
-  END DO
-  CALL NF95_CLOSE(fID_in)
-
-  !=============================================================================
-  IF(l3D) THEN                                                   !=== 3D FIELDS
-  !=============================================================================
-    !--- Regrid in longitude
-    ALLOCATE(o3_regr_lon(nlon_ou, nlat_in, nlev_in, ie-ib+1, read_climoz))
-    CALL regr_conserv(1, o3_in3, xs = lon_in_edge,                             &
-                        xt = [boundslon_reg(1,west),boundslon_reg(:,east)],    &
-                        vt = o3_regr_lon, slope = slopes(1,o3_in3, lon_in_edge))
-    DEALLOCATE(o3_in3)
-
-    !--- Regrid in latitude: averaging with respect to SIN(lat) is
-    !                        equivalent to weighting by COS(lat)
-    !--- (inverted indices in "o3_regr_lonlat" because "rlatu" is decreasing)
-    ALLOCATE(o3_regr_lonlat(nlon_ou, nlat_ou, nlev_in, 0:13, read_climoz))
-    CALL regr_conserv(2, o3_regr_lon, xs = sinlat_in_edge,                     &
-                    xt = [- 1., SIN(boundslat_reg(nlat_ou-1:1:-1,south)), 1.], &
-                    vt = o3_regr_lonlat(:,nlat_ou:1:- 1,:,ib:ie,:),            &
-                 slope = slopes(2,o3_regr_lon, sinlat_in_edge))
-    DEALLOCATE(o3_regr_lon)
-
-    !--- Duplicate previous/next record(s) if they are not available
-    IF(.NOT.lprev) o3_regr_lonlat(:,:,:, 0,:) = o3_regr_lonlat(:,:,:,12,:)
-    IF(.NOT.lnext) o3_regr_lonlat(:,:,:,13,:) = o3_regr_lonlat(:,:,:, 1,:)
-
-    !--- Regrid in time by linear interpolation:
-    ALLOCATE(o3_out3(nlon_ou, nlat_ou, nlev_in, ntim_ou, read_climoz))
-    IF(     interpt) CALL regr_lint(4,o3_regr_lonlat,tmidmonth,tmidday,o3_out3)
-    IF(.NOT.interpt) o3_out3=o3_regr_lonlat
-    DEALLOCATE(o3_regr_lonlat)
-
-    !--- Write to file (the order of "rlatu" is inverted in the output file):
-    DO m = 1, read_climoz
-      CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out3(:,nlat_ou:1:-1,:,:,m))
-    END DO
-
-  !=============================================================================
-  ELSE                                                         !=== ZONAL FIELDS
-  !=============================================================================
-    !--- Regrid in latitude: averaging with respect to SIN(lat) is
-    !                        equivalent to weighting by COS(lat)
-    !--- (inverted indices in "o3_regr_lat" because "rlatu" is decreasing)
-    ALLOCATE(o3_regr_lat(nlat_ou, nlev_in, 0:13, read_climoz))
-    CALL regr_conserv(1, o3_in2, xs = sinlat_in_edge,                          &
-                    xt = [- 1., SIN(boundslat_reg(nlat_ou-1:1:-1,south)), 1.], &
-                    vt = o3_regr_lat(nlat_ou:1:- 1,:,ib:ie,:),                 &
-                 slope = slopes(1,o3_in2, sinlat_in_edge))
-    DEALLOCATE(o3_in2)
-
-    !--- Duplicate previous/next record(s) if they are not available
-    IF(.NOT.lprev) o3_regr_lat(:,:, 0,:) = o3_regr_lat(:,:,12,:)
-    IF(.NOT.lnext) o3_regr_lat(:,:,13,:) = o3_regr_lat(:,:, 1,:)
-
-    !--- Regrid in time by linear interpolation:
-    ALLOCATE(o3_out2(nlat_ou, nlev_in, ntim_ou, read_climoz))
-    IF(     interpt) CALL regr_lint(3,o3_regr_lat, tmidmonth, tmidday, o3_out2)
-    IF(.NOT.interpt) o3_out2=o3_regr_lat
-    DEALLOCATE(o3_regr_lat)
-
-    !--- Write to file (the order of "rlatu" is inverted in the output file):
-    DO m = 1, read_climoz
-      CALL NF95_PUT_VAR(fID_ou, vID_ou(m), o3_out2(nlat_ou:1:-1,:,:,m))
-    END DO
-
-  !=============================================================================
-  END IF
-  !=============================================================================
-
-  CALL NF95_CLOSE(fID_ou)
-
+    !=============================================================================
+
+    IF (is_mpi_root) CALL NF95_CLOSE(fID_in)
+
+  ENDIF ! is_omp_master
+
+  first=.FALSE.
 END SUBROUTINE regr_horiz_time_climoz
 !
@@ -408,4 +607,8 @@
 !-------------------------------------------------------------------------------
   USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
+  USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
+  USE mod_phys_lmdz_para, ONLY: is_mpi_root
+  USE mod_grid_phy_lmdz, ONLY: klon_glo
+!
 !-------------------------------------------------------------------------------
 ! Arguments:
@@ -419,68 +622,83 @@
   INTEGER :: dlonID, dlatID, dlevID, dtimID, dIDs(4)
   INTEGER :: vlonID, vlatID, ncerr,  is
+  REAL,ALLOCATABLE    :: latitude_glo_(:)
   CHARACTER(LEN=80) :: sub
-!-------------------------------------------------------------------------------
-  sub="prepare_out"
-  WRITE(lunout,*)"CALL sequence information: "//TRIM(sub)
-  CALL NF95_CREATE("climoz_LMDZ.nc", NF90_clobber, fID_ou)
+  INTEGER :: i
+
+
+!-------------------------------------------------------------------------------
+  
+  IF (is_mpi_root) THEN  
+    sub="prepare_out"
+    WRITE(lunout,*)"CALL sequence information: "//TRIM(sub)
+    CALL NF95_CREATE("climoz_LMDZ.nc", NF90_clobber, fID_ou)
 
   !--- Dimensions:
-  IF(ndims==4) &
-  CALL NF95_DEF_DIM(fID_ou, "rlonv", nlon_ou, dlonID)
-  CALL NF95_DEF_DIM(fID_ou, "rlatu", nlat_ou, dlatID)
-  CALL NF95_DEF_DIM(fID_ou, "plev",  nlev_in, dlevID)
-  CALL NF95_DEF_DIM(fID_ou, "time",  ntim_ou, dtimID)
-
-  !--- Define coordinate variables:
-  IF(ndims==4) &
-  CALL NF95_DEF_VAR(fID_ou, "rlonv", NF90_FLOAT, dlonID, vlonID)
-  CALL NF95_DEF_VAR(fID_ou, "rlatu", NF90_FLOAT, dlatID, vlatID)
-  CALL NF95_DEF_VAR(fID_ou, "plev",  NF90_FLOAT, dlevID, vlevID)
-  CALL NF95_DEF_VAR(fID_ou, "time",  NF90_FLOAT, dtimID, vtimID)
-  IF(ndims==4) &
-  CALL NF95_PUT_ATT(fID_ou, vlonID, "units", "degrees_east")
-  CALL NF95_PUT_ATT(fID_ou, vlatID, "units", "degrees_north")
-  CALL NF95_PUT_ATT(fID_ou, vlevID, "units", "millibar")
-  CALL NF95_PUT_ATT(fID_ou, vtimID, "units", "days since 2000-1-1")
-  IF(ndims==4) &
-  CALL NF95_PUT_ATT(fID_ou, vlonID, "standard_name", "longitude")
-  CALL NF95_PUT_ATT(fID_ou, vlatID, "standard_name", "latitude")
-  CALL NF95_PUT_ATT(fID_ou, vlevID, "standard_name", "air_pressure")
-  CALL NF95_PUT_ATT(fID_ou, vtimID, "standard_name", "time")
-  CALL NF95_PUT_ATT(fID_ou, vlevID, "long_name",     "air pressure")
-  CALL NF95_PUT_ATT(fID_ou, vtimID, "calendar",      cal_ou)
+    IF(ndims==4) &
+    CALL NF95_DEF_DIM(fID_ou, "rlonv", nlon_ou, dlonID)
+    CALL NF95_DEF_DIM(fID_ou, "rlatu", nlat_ou, dlatID)
+    CALL NF95_DEF_DIM(fID_ou, "plev",  nlev_in, dlevID)
+    CALL NF95_DEF_DIM(fID_ou, "time",  ntim_ou, dtimID)
+
+    !--- Define coordinate variables:
+    IF(ndims==4) &
+    CALL NF95_DEF_VAR(fID_ou, "rlonv", NF90_FLOAT, dlonID, vlonID)
+    CALL NF95_DEF_VAR(fID_ou, "rlatu", NF90_FLOAT, dlatID, vlatID)
+    CALL NF95_DEF_VAR(fID_ou, "plev",  NF90_FLOAT, dlevID, vlevID)
+    CALL NF95_DEF_VAR(fID_ou, "time",  NF90_FLOAT, dtimID, vtimID)
+    IF(ndims==4) &
+    CALL NF95_PUT_ATT(fID_ou, vlonID, "units", "degrees_east")
+    CALL NF95_PUT_ATT(fID_ou, vlatID, "units", "degrees_north")
+    CALL NF95_PUT_ATT(fID_ou, vlevID, "units", "millibar")
+    CALL NF95_PUT_ATT(fID_ou, vtimID, "units", "days since 2000-1-1")
+    IF(ndims==4) &
+    CALL NF95_PUT_ATT(fID_ou, vlonID, "standard_name", "longitude")
+    CALL NF95_PUT_ATT(fID_ou, vlatID, "standard_name", "latitude")
+    CALL NF95_PUT_ATT(fID_ou, vlevID, "standard_name", "air_pressure")
+    CALL NF95_PUT_ATT(fID_ou, vtimID, "standard_name", "time")
+    CALL NF95_PUT_ATT(fID_ou, vlevID, "long_name",     "air pressure")
+    CALL NF95_PUT_ATT(fID_ou, vtimID, "calendar",      cal_ou)
 
   !--- Define the main variables:
-  IF(ndims==3) dIDs(1:3) = [ dlatID, dlevID, dtimID]
-  IF(ndims==4) dIDs=[dlonID, dlatID, dlevID, dtimID]
-  CALL NF95_DEF_VAR(fID_ou, vars_in(1), NF90_FLOAT, dIDs(1:ndims), vID_ou(1))
-  CALL NF95_PUT_ATT(fID_ou, vID_ou(1), "long_name", "ozone mole fraction")
-  CALL NF95_PUT_ATT(fID_ou, vID_ou(1), "standard_name", "mole_fraction_of_ozone&
+    IF(ndims==3) dIDs(1:3) = [ dlatID, dlevID, dtimID]
+    IF(ndims==4) dIDs=[dlonID, dlatID, dlevID, dtimID]
+    CALL NF95_DEF_VAR(fID_ou, vars_in(1), NF90_FLOAT, dIDs(1:ndims), vID_ou(1))
+    CALL NF95_PUT_ATT(fID_ou, vID_ou(1), "long_name", "ozone mole fraction")
+    CALL NF95_PUT_ATT(fID_ou, vID_ou(1), "standard_name", "mole_fraction_of_ozone&
       &_in_air")
-  IF(SIZE(vID_ou) == 2) THEN
-    CALL NF95_DEF_VAR(fID_ou, vars_in(2), NF90_FLOAT, dIDs(1:ndims), vID_ou(2))
-    CALL NF95_PUT_ATT(fID_ou, vID_ou(2), "long_name","ozone mole fraction in da&
-      &ylight")
-  END IF
+    IF(SIZE(vID_ou) == 2) THEN
+      CALL NF95_DEF_VAR(fID_ou, vars_in(2), NF90_FLOAT, dIDs(1:ndims), vID_ou(2))
+      CALL NF95_PUT_ATT(fID_ou, vID_ou(2), "long_name","ozone mole fraction in da&
+        &ylight")
+    END IF
 
   !--- Global attributes:
   ! The following commands, copying attributes, may fail. That is OK.
   ! It should just mean that the attribute is not defined in the input file.
-  CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"Conventions",fID_ou,NF90_GLOBAL, ncerr)
-  CALL handle_err_copy_att("Conventions")
-  CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"title",      fID_ou,NF90_GLOBAL, ncerr)
-  CALL handle_err_copy_att("title")
-  CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"institution",fID_ou,NF90_GLOBAL, ncerr)
-  CALL handle_err_copy_att("institution")
-  CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"source",     fID_ou,NF90_GLOBAL, ncerr)
-  CALL handle_err_copy_att("source")
-  CALL NF95_PUT_ATT (fID_ou,NF90_GLOBAL,"comment", "Regridded for LMDZ")
-  CALL NF95_ENDDEF(fID_ou)
-
-  !--- Write one of the coordinate variables:
-  IF(ndims==4) CALL NF95_PUT_VAR(fID_ou, vlonID, lon_reg/deg2rad)
-  CALL NF95_PUT_VAR(fID_ou, vlatID, lat_reg(nlat_ou:1:-1)/deg2rad)
-  !    (convert from rad to degrees and sort in ascending order)
-
+    CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"Conventions",fID_ou,NF90_GLOBAL, ncerr)
+    CALL handle_err_copy_att("Conventions")
+    CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"title",      fID_ou,NF90_GLOBAL, ncerr)
+    CALL handle_err_copy_att("title")
+    CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"institution",fID_ou,NF90_GLOBAL, ncerr)
+    CALL handle_err_copy_att("institution")
+    CALL NF95_COPY_ATT(fID_in,NF90_GLOBAL,"source",     fID_ou,NF90_GLOBAL, ncerr)
+    CALL handle_err_copy_att("source")
+    CALL NF95_PUT_ATT (fID_ou,NF90_GLOBAL,"comment", "Regridded for LMDZ")
+    CALL NF95_ENDDEF(fID_ou)
+
+    IF (grid_type==unstructured) THEN
+      ALLOCATE(latitude_glo_(klon_glo))
+      DO i=1,klon_glo
+        latitude_glo_(ind_cell_glo_glo(i))=latitude_glo(i)
+      ENDDO
+      CALL NF95_PUT_VAR(fID_ou, vlatID, latitude_glo_)
+    ELSE
+      !--- Write one of the coordinate variables:
+      IF(ndims==4) CALL NF95_PUT_VAR(fID_ou, vlonID, lon_reg/deg2rad)
+      CALL NF95_PUT_VAR(fID_ou, vlatID, lat_reg(nlat_ou:1:-1)/deg2rad)
+    !    (convert from rad to degrees and sort in ascending order)
+    ENDIF
+  ENDIF
+  
 CONTAINS
 
Index: /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/regr_pr_time_av_m.F90
===================================================================
--- /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/regr_pr_time_av_m.F90	(revision 3412)
+++ /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/regr_pr_time_av_m.F90	(revision 3413)
@@ -106,13 +106,14 @@
   USE assert_m,       ONLY: assert
   USE assert_eq_m,    ONLY: assert_eq
-!  USE comvert_mod,    ONLY: scaleheight
+!!  USE comvert_mod,    ONLY: scaleheight
   USE interpolation,  ONLY: locate
   USE regr_conserv_m, ONLY: regr_conserv
   USE regr_lint_m,    ONLY: regr_lint
   USE slopes_m,       ONLY: slopes
-  USE mod_phys_lmdz_mpi_data,       ONLY: is_mpi_root
-  USE mod_grid_phy_lmdz,            ONLY: nbp_lon, nbp_lat, nbp_lev
-  USE mod_phys_lmdz_transfert_para, ONLY: scatter2d, scatter
+  USE mod_phys_lmdz_para,           ONLY: is_mpi_root,is_master
+  USE mod_grid_phy_lmdz,            ONLY: nbp_lon, nbp_lat, nbp_lev, klon_glo, grid_type, unstructured
+  USE mod_phys_lmdz_transfert_para, ONLY: scatter2d, scatter, gather
   USE phys_cal_mod,                 ONLY: calend, year_len, days_elapsed, jH_cur
+  USE geometry_mod,                 ONLY: ind_cell_glo
 !-------------------------------------------------------------------------------
 ! Arguments:
@@ -157,4 +158,5 @@
   REAL, DIMENSION(nbp_lon, nbp_lat)   :: ps1, pt1, ot1
   REAL, DIMENSION(klon)               :: ps2, pt2, ot2, ptropou
+  INTEGER,ALLOCATABLE :: ind_cell_glo_glo(:)
   LOGICAL :: ll
 !-------------------------------------------------------------------------------
@@ -230,10 +232,32 @@
     CALL bcast(lO3Trop); CALL bcast(linterp)
   END IF
+  
+  IF (is_master) THEN 
+    ALLOCATE(ind_cell_glo_glo(klon_glo))
+  ELSE
+    ALLOCATE(ind_cell_glo_glo(0))
+  ENDIF
+  CALL gather(ind_cell_glo,ind_cell_glo_glo)
+  IF (is_master .AND. grid_type==unstructured) v1(:,:,:,:)=v1(:,ind_cell_glo_glo(:),:,:)
+  
   CALL scatter2d(v1,v2)
+
   !--- No "ps" in input file => assumed to be equal to current LMDZ ground press
-  IF(lPrSurf) THEN; CALL scatter2d(ps1,ps2); ELSE; ps2=pint_ou(:,1); END IF
-  IF(lPrTrop) CALL scatter2d(pt1,pt2)
-  IF(lO3Trop) CALL scatter2d(ot1,ot2)
-
+  IF(lPrSurf) THEN
+    IF (is_master .AND. grid_type==unstructured) ps1(:,:)=ps1(:,ind_cell_glo_glo(:))
+    CALL scatter2d(ps1,ps2)
+  ELSE
+   ps2=pint_ou(:,1)
+  END IF
+
+  IF(lPrTrop) THEN
+    IF (is_master .AND. grid_type==unstructured) pt1(:,:)=pt1(:,ind_cell_glo_glo(:))
+    CALL scatter2d(pt1,pt2)
+  ENDIF
+
+  IF(lO3Trop) THEN
+    IF (is_master .AND. grid_type==unstructured) ot1(:,:)=ot1(:,ind_cell_glo_glo(:))
+    CALL scatter2d(ot1,ot2)
+  ENDIF
   !--- REGRID IN PRESSURE ; 3rd index inverted because "paprs" is decreasing
   IF(.NOT.lAdjTro) THEN
Index: /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/rrtm/readaerosolstrato1_rrtm.F90
===================================================================
--- /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/rrtm/readaerosolstrato1_rrtm.F90	(revision 3412)
+++ /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/rrtm/readaerosolstrato1_rrtm.F90	(revision 3413)
@@ -2,4 +2,5 @@
 ! $Id: readaerosolstrato1_rrtm.F90 2526 2016-05-26 22:13:40Z oboucher $
 !
+
 SUBROUTINE readaerosolstrato1_rrtm(debut)
 
@@ -9,7 +10,5 @@
 
     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_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, klon_glo, grid2dTo1d_glo, grid_type, unstructured
     USE mod_phys_lmdz_para 
     USE phys_state_var_mod
@@ -19,4 +18,5 @@
     USE YOERAD, ONLY : NLW
     USE YOMCST
+    USE xios
 
     IMPLICIT NONE
@@ -45,4 +45,5 @@
     REAL, ALLOCATABLE:: tauaerstrat_mois(:, :, :)
     REAL, ALLOCATABLE:: tauaerstrat_mois_glo(:, :)
+    REAL, ALLOCATABLE:: tauaerstrat_mpi(:, :)
 
 ! For NetCDF:
@@ -102,18 +103,25 @@
     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
-
+
+    IF (grid_type/=unstructured) THEN
+      IF (n_lat.NE.nbp_lat) THEN 
+         print *,'Le nombre de lat n est pas egal a nbp_lat'
+         STOP
+      ENDIF
+    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
-
+
+    IF (grid_type/=unstructured) THEN
+      IF (n_lon.NE.nbp_lon) THEN 
+         print *,'Le nombre de lon n est pas egal a nbp_lon'
+         STOP
+      ENDIF
+    ENDIF
+    
+    
     CALL nf95_inq_varid(ncid_in, "TIME", varid)
     CALL nf95_gw_var(ncid_in, varid, time)
@@ -144,5 +152,7 @@
 !---reduce to a klon_glo grid 
     CALL grid2dTo1d_glo(tauaerstrat_mois,tauaerstrat_mois_glo)
-
+    
+    ELSE
+      ALLOCATE(tauaerstrat_mois(0,0,0))
     ENDIF !--is_mpi_root and is_omp_root
 
@@ -153,6 +163,18 @@
 
 !--scatter on all proc
-    CALL scatter(tauaerstrat_mois_glo,tau_aer_strat)
-
+    
+    IF (grid_type==unstructured) THEN
+      IF (is_omp_master) THEN
+        ALLOCATE(tauaerstrat_mpi(klon_mpi,klev))
+        CALL xios_send_field("taustrat_in",tauaerstrat_mois)
+        CALL xios_recv_field("taustrat_out",tauaerstrat_mpi)
+      ELSE
+        ALLOCATE(tauaerstrat_mpi(0,0))
+      ENDIF
+      CALL scatter_omp(tauaerstrat_mpi,tau_aer_strat)
+    ELSE  
+      CALL scatter(tauaerstrat_mois_glo,tau_aer_strat)
+    ENDIF
+    
     IF (is_mpi_root.AND.is_omp_root) THEN
 !
Index: /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/rrtm/readaerosolstrato2_rrtm.F90
===================================================================
--- /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/rrtm/readaerosolstrato2_rrtm.F90	(revision 3412)
+++ /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/rrtm/readaerosolstrato2_rrtm.F90	(revision 3413)
@@ -2,4 +2,5 @@
 ! $Id: readaerosolstrato2_rrtm.F90 2526 2016-05-26 22:13:40Z oboucher $
 !
+
 SUBROUTINE readaerosolstrato2_rrtm(debut)
 
@@ -9,7 +10,7 @@
 
     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_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, klon_glo, grid2dTo1d_glo, grid_type, unstructured
+    USE mod_phys_lmdz_mpi_data
+    USE mod_phys_lmdz_omp_data
     USE mod_phys_lmdz_para 
     USE phys_state_var_mod
@@ -19,5 +20,5 @@
     USE YOERAD, ONLY : NLW
     USE YOMCST
-
+    USE xios
     IMPLICIT NONE
 
@@ -65,4 +66,8 @@
     REAL, ALLOCATABLE:: cgaerstrat_mois_glo(:, :, :)
     REAL, ALLOCATABLE:: taulwaerstrat_mois_glo(:, :, :)
+    REAL, ALLOCATABLE:: tauaerstrat_mpi(:, :, :)
+    REAL, ALLOCATABLE:: pizaerstrat_mpi(:, :, :)
+    REAL, ALLOCATABLE:: cgaerstrat_mpi(:, :, :)
+    REAL, ALLOCATABLE:: taulwaerstrat_mpi(:, :, :)
 
 ! For NetCDF:
@@ -107,8 +112,11 @@
         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)
+
+        IF (grid_type/=unstructured) THEN
+           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
         ENDIF
 
@@ -134,12 +142,4 @@
         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)
@@ -159,26 +159,36 @@
         CALL nf95_close(ncid_in)
 
+       
+        IF (grid_type/=unstructured) THEN
+          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))
 !--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)
+          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
+          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, NSW
-          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
-
+          DO band=1, NSW
+            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
+        ENDIF
 !--Now LW optical properties
 !
+
         CALL nf95_open("taulwstrat.2D.nc", nf90_nowrite, ncid_in)
 
@@ -194,9 +204,12 @@
         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
-
+
+        IF (grid_type/=unstructured) THEN
+          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
+        ENDIF
+        
         CALL nf95_inq_varid(ncid_in, "TIME", varid)
         CALL nf95_gw_var(ncid_in, varid, time)
@@ -217,6 +230,4 @@
 
         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
@@ -227,17 +238,23 @@
         CALL nf95_close(ncid_in)
 
+        IF (grid_type/=unstructured) THEN
+
+          ALLOCATE(taulwaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
+          ALLOCATE(taulwaerstrat_mois_glo(klon_glo, n_lev, n_wav))
+
 !--select the correct month
 !--and copy into 1st longitude
-        taulwaerstrat_mois(1,:,:,:) = taulwaerstrat(:,:,:,mth_cur)
+          taulwaerstrat_mois(1,:,:,:) = taulwaerstrat(:,:,:,mth_cur)
 !--copy longitudes
-        DO i=2, n_lon
-          taulwaerstrat_mois(i,:,:,:) = taulwaerstrat_mois(1,:,:,:)
-        ENDDO
+          DO i=2, n_lon
+            taulwaerstrat_mois(i,:,:,:) = taulwaerstrat_mois(1,:,:,:)
+          ENDDO
 
 !---reduce to a klon_glo grid 
-        DO band=1, NLW
-          CALL grid2dTo1d_glo(taulwaerstrat_mois(:,:,:,band),taulwaerstrat_mois_glo(:,:,band))
-        ENDDO
-
+          DO band=1, NLW
+            CALL grid2dTo1d_glo(taulwaerstrat_mois(:,:,:,band),taulwaerstrat_mois_glo(:,:,band))
+          ENDDO
+        ENDIF
+      
       ELSE !--proc other than mpi_root and omp_root
            !--dummy allocation needed for debug mode
@@ -248,4 +265,10 @@
         ALLOCATE(taulwaerstrat_mois_glo(1,1,1))
 
+        ALLOCATE(tauaerstrat(0,0,0,12))
+        ALLOCATE(pizaerstrat(0,0,0,12))
+        ALLOCATE(cgaerstrat(0,0,0,12))
+        ALLOCATE(taulwaerstrat(0,0,0,12))
+
+
       ENDIF !--is_mpi_root and is_omp_root
 
@@ -255,20 +278,47 @@
       mth_pre=mth_cur
 
+      IF (grid_type==unstructured) THEN
+
+        IF (is_omp_master) THEN
+          ALLOCATE(tauaerstrat_mpi(klon_mpi, klev, NSW))
+          ALLOCATE(pizaerstrat_mpi(klon_mpi, klev, NSW))
+          ALLOCATE(cgaerstrat_mpi(klon_mpi, klev, NSW))        
+          ALLOCATE(taulwaerstrat_mpi(klon_mpi, klev, NLW))
+         
+          CALL xios_send_field("tauaerstrat_in",SPREAD(tauaerstrat(:,:,:,mth_cur),1,8))
+          CALL xios_recv_field("tauaerstrat_out",tauaerstrat_mpi)
+          CALL xios_send_field("pizaerstrat_in",SPREAD(pizaerstrat(:,:,:,mth_cur),1,8))
+          CALL xios_recv_field("pizaerstrat_out",pizaerstrat_mpi)
+          CALL xios_send_field("cgaerstrat_in",SPREAD(cgaerstrat(:,:,:,mth_cur),1,8))
+          CALL xios_recv_field("cgaerstrat_out",cgaerstrat_mpi)
+          CALL xios_send_field("taulwaerstrat_in",SPREAD(taulwaerstrat(:,:,:,mth_cur),1,8))
+          CALL xios_recv_field("taulwaerstrat_out",taulwaerstrat_mpi)
+        ELSE
+          ALLOCATE(tauaerstrat_mpi(0, 0, 0))
+          ALLOCATE(pizaerstrat_mpi(0, 0, 0))
+          ALLOCATE(cgaerstrat_mpi(0, 0, 0))        
+          ALLOCATE(taulwaerstrat_mpi(0, 0, 0))
+        ENDIF  
+        
+        CALL scatter_omp(tauaerstrat_mpi,tau_aer_strat)
+        CALL scatter_omp(pizaerstrat_mpi,piz_aer_strat)
+        CALL scatter_omp(cgaerstrat_mpi,cg_aer_strat)
+        CALL scatter_omp(taulwaerstrat_mpi,taulw_aer_strat)
+      ELSE  
+        
 !--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)
+        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)
+
+        DEALLOCATE(tauaerstrat_mois, pizaerstrat_mois, cgaerstrat_mois)
+        DEALLOCATE(taulwaerstrat_mois)
+      ENDIF 
 
       IF (is_mpi_root.AND.is_omp_root) THEN
-!
-        DEALLOCATE(tauaerstrat, pizaerstrat, cgaerstrat)
-        DEALLOCATE(tauaerstrat_mois, pizaerstrat_mois, cgaerstrat_mois)
-        DEALLOCATE(taulwaerstrat,taulwaerstrat_mois)
-!
-      ENDIF !--is_mpi_root and is_omp_root
-
-      DEALLOCATE(tauaerstrat_mois_glo,pizaerstrat_mois_glo,cgaerstrat_mois_glo)
-      DEALLOCATE(taulwaerstrat_mois_glo)
+        DEALLOCATE(tauaerstrat, pizaerstrat, cgaerstrat,taulwaerstrat)
+      ENDIF
+      
 
 !$OMP BARRIER
Index: /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/surf_land_mod.F90
===================================================================
--- /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/surf_land_mod.F90	(revision 3412)
+++ /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/surf_land_mod.F90	(revision 3413)
@@ -38,5 +38,10 @@
     USE surf_land_orchidee_nofrein_mod
 #else
+#if ORCHIDEE_NOUNSTRUCT
+    ! Compilation with cpp key ORCHIDEE_NOUNSTRUCT
+    USE surf_land_orchidee_nounstruct_mod
+#else
     USE surf_land_orchidee_mod
+#endif
 #endif
 #endif
Index: /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/surf_land_orchidee_mod.F90
===================================================================
--- /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/surf_land_orchidee_mod.F90	(revision 3412)
+++ /LMDZ6/branches/DYNAMICO-conv/libf/phylmd/surf_land_orchidee_mod.F90	(revision 3413)
@@ -4,4 +4,5 @@
 #ifndef ORCHIDEE_NOZ0H
 #ifndef ORCHIDEE_NOFREIN
+#ifndef ORCHIDEE_NOUNSTRUCT
 !
 ! This module controles the interface towards the model ORCHIDEE.
@@ -23,8 +24,9 @@
   USE cpl_mod,      ONLY : cpl_send_land_fields
   USE surface_data, ONLY : type_ocean
-  USE geometry_mod, ONLY : dx, dy
+  USE geometry_mod, ONLY : dx, dy, boundslon, boundslat,longitude, latitude, cell_area,  ind_cell_glo
   USE mod_grid_phy_lmdz
   USE mod_phys_lmdz_para, mpi_root_rank=>mpi_master
-
+  USE nrtype, ONLY : PI
+  
   IMPLICIT NONE
 
@@ -170,4 +172,7 @@
     REAL, ALLOCATABLE, DIMENSION(:,:), SAVE   :: lalo
     !$OMP THREADPRIVATE(lalo)
+! boundaries of cells
+    REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE   :: bounds_lalo
+    !$OMP THREADPRIVATE(bounds_lalo)
 ! pts voisins
     INTEGER,ALLOCATABLE, DIMENSION(:,:), SAVE :: neighbours
@@ -183,4 +188,8 @@
     !$OMP THREADPRIVATE(lon_scat,lat_scat)
 
+! area of cells
+    REAL, ALLOCATABLE, DIMENSION (:), SAVE  :: area  
+    !$OMP THREADPRIVATE(area)
+
     LOGICAL, SAVE                             :: lrestart_read = .TRUE.
     !$OMP THREADPRIVATE(lrestart_read)
@@ -214,6 +223,17 @@
     !$OMP THREADPRIVATE(riverflow)
     
+    INTEGER :: orch_mpi_rank
+    INTEGER :: orch_mpi_size
     INTEGER :: orch_omp_rank
     INTEGER :: orch_omp_size
+
+    REAL, ALLOCATABLE, DIMENSION(:)         :: longitude_glo
+    REAL, ALLOCATABLE, DIMENSION(:)         :: latitude_glo
+    REAL, ALLOCATABLE, DIMENSION(:,:)       :: boundslon_glo
+    REAL, ALLOCATABLE, DIMENSION(:,:)       :: boundslat_glo
+    INTEGER, ALLOCATABLE, DIMENSION(:)      :: ind_cell_glo_glo
+    INTEGER, ALLOCATABLE, SAVE,DIMENSION(:) :: ind_cell
+    !$OMP THREADPRIVATE(ind_cell)
+    INTEGER :: begin, end
 !
 ! Fin definition
@@ -258,66 +278,123 @@
        jg(klon) = nbp_lat
 
-       IF ((.NOT. ALLOCATED(lalo))) THEN
-          ALLOCATE(lalo(knon,2), stat = error)
+       IF ((.NOT. ALLOCATED(area))) THEN
+          ALLOCATE(area(knon), stat = error)
           IF (error /= 0) THEN
+             abort_message='Pb allocation area'
+             CALL abort_physic(modname,abort_message,1)
+          ENDIF
+       ENDIF
+       DO igrid = 1, knon
+          area(igrid) = cell_area(knindex(igrid))
+       ENDDO
+       
+       IF (grid_type==unstructured) THEN
+
+
+         IF ((.NOT. ALLOCATED(lon_scat))) THEN
+            ALLOCATE(lon_scat(nbp_lon,nbp_lat), stat = error)
+            IF (error /= 0) THEN
+               abort_message='Pb allocation lon_scat'
+               CALL abort_physic(modname,abort_message,1)
+            ENDIF
+         ENDIF
+ 
+         IF ((.NOT. ALLOCATED(lat_scat))) THEN
+            ALLOCATE(lat_scat(nbp_lon,nbp_lat), stat = error)
+            IF (error /= 0) THEN
+               abort_message='Pb allocation lat_scat'
+               CALL abort_physic(modname,abort_message,1)
+            ENDIF
+         ENDIF
+         CALL Gather(rlon,rlon_g)
+         CALL Gather(rlat,rlat_g)
+
+         IF (is_mpi_root) THEN
+            index = 1
+            DO jj = 2, nbp_lat-1
+               DO ij = 1, nbp_lon
+                  index = index + 1
+                  lon_scat(ij,jj) = rlon_g(index)
+                  lat_scat(ij,jj) = rlat_g(index)
+               ENDDO
+            ENDDO
+            lon_scat(:,1) = lon_scat(:,2)
+            lat_scat(:,1) = rlat_g(1)
+            lon_scat(:,nbp_lat) = lon_scat(:,2)
+            lat_scat(:,nbp_lat) = rlat_g(klon_glo)
+         ENDIF
+     
+         CALL bcast(lon_scat)
+         CALL bcast(lat_scat)
+                
+       ELSE IF (grid_type==regular_lonlat) THEN
+
+         IF ((.NOT. ALLOCATED(lalo))) THEN
+            ALLOCATE(lalo(knon,2), stat = error)
+            IF (error /= 0) THEN
+               abort_message='Pb allocation lalo'
+               CALL abort_physic(modname,abort_message,1)
+            ENDIF
+         ENDIF
+       
+         IF ((.NOT. ALLOCATED(bounds_lalo))) THEN
+           ALLOCATE(bounds_lalo(knon,nvertex,2), stat = error)
+           IF (error /= 0) THEN
              abort_message='Pb allocation lalo'
              CALL abort_physic(modname,abort_message,1)
-          ENDIF
-       ENDIF
-       IF ((.NOT. ALLOCATED(lon_scat))) THEN
-          ALLOCATE(lon_scat(nbp_lon,nbp_lat), stat = error)
-          IF (error /= 0) THEN
-             abort_message='Pb allocation lon_scat'
-             CALL abort_physic(modname,abort_message,1)
-          ENDIF
-       ENDIF
-       IF ((.NOT. ALLOCATED(lat_scat))) THEN
-          ALLOCATE(lat_scat(nbp_lon,nbp_lat), stat = error)
-          IF (error /= 0) THEN
-             abort_message='Pb allocation lat_scat'
-             CALL abort_physic(modname,abort_message,1)
-          ENDIF
-       ENDIF
-       lon_scat = 0.
-       lat_scat = 0.
-       DO igrid = 1, knon
-          index = knindex(igrid)
-          lalo(igrid,2) = rlon(index)
-          lalo(igrid,1) = rlat(index)
-       ENDDO
-
-       
-       
-       CALL Gather(rlon,rlon_g)
-       CALL Gather(rlat,rlat_g)
-
-       IF (is_mpi_root) THEN
-          index = 1
-          DO jj = 2, nbp_lat-1
-             DO ij = 1, nbp_lon
-                index = index + 1
-                lon_scat(ij,jj) = rlon_g(index)
-                lat_scat(ij,jj) = rlat_g(index)
-             ENDDO
-          ENDDO
-          lon_scat(:,1) = lon_scat(:,2)
-          lat_scat(:,1) = rlat_g(1)
-          lon_scat(:,nbp_lat) = lon_scat(:,2)
-          lat_scat(:,nbp_lat) = rlat_g(klon_glo)
-       ENDIF
+           ENDIF
+         ENDIF
+       
+         IF ((.NOT. ALLOCATED(lon_scat))) THEN
+            ALLOCATE(lon_scat(nbp_lon,nbp_lat), stat = error)
+            IF (error /= 0) THEN
+               abort_message='Pb allocation lon_scat'
+               CALL abort_physic(modname,abort_message,1)
+            ENDIF
+         ENDIF
+         IF ((.NOT. ALLOCATED(lat_scat))) THEN
+            ALLOCATE(lat_scat(nbp_lon,nbp_lat), stat = error)
+            IF (error /= 0) THEN
+               abort_message='Pb allocation lat_scat'
+               CALL abort_physic(modname,abort_message,1)
+            ENDIF
+         ENDIF
+         lon_scat = 0.
+         lat_scat = 0.
+         DO igrid = 1, knon
+            index = knindex(igrid)
+            lalo(igrid,2) = rlon(index)
+            lalo(igrid,1) = rlat(index)
+            bounds_lalo(igrid,:,2)=boundslon(index,:)*180./PI
+            bounds_lalo(igrid,:,1)=boundslat(index,:)*180./PI
+         ENDDO
+
+       
+       
+         CALL Gather(rlon,rlon_g)
+         CALL Gather(rlat,rlat_g)
+
+         IF (is_mpi_root) THEN
+            index = 1
+            DO jj = 2, nbp_lat-1
+               DO ij = 1, nbp_lon
+                  index = index + 1
+                  lon_scat(ij,jj) = rlon_g(index)
+                  lat_scat(ij,jj) = rlat_g(index)
+               ENDDO
+            ENDDO
+            lon_scat(:,1) = lon_scat(:,2)
+            lat_scat(:,1) = rlat_g(1)
+            lon_scat(:,nbp_lat) = lon_scat(:,2)
+            lat_scat(:,nbp_lat) = rlat_g(klon_glo)
+         ENDIF
    
-       CALL bcast(lon_scat)
-       CALL bcast(lat_scat)
+         CALL bcast(lon_scat)
+         CALL bcast(lat_scat)
+       
+       ENDIF
 !
 ! Allouer et initialiser le tableau des voisins et des fraction de continents
 !
-       IF ( (.NOT.ALLOCATED(neighbours))) THEN
-          ALLOCATE(neighbours(knon,8), stat = error)
-          IF (error /= 0) THEN
-             abort_message='Pb allocation neighbours'
-             CALL abort_physic(modname,abort_message,1)
-          ENDIF
-       ENDIF
-       neighbours = -1.
        IF (( .NOT. ALLOCATED(contfrac))) THEN
           ALLOCATE(contfrac(knon), stat = error)
@@ -334,5 +411,29 @@
 
 
-       CALL Init_neighbours(knon,neighbours,knindex,pctsrf(:,is_ter))
+       IF (grid_type==regular_lonlat) THEN
+ 
+         IF ( (.NOT.ALLOCATED(neighbours))) THEN
+          ALLOCATE(neighbours(knon,8), stat = error)
+          IF (error /= 0) THEN
+             abort_message='Pb allocation neighbours'
+             CALL abort_physic(modname,abort_message,1)
+          ENDIF
+         ENDIF
+         neighbours = -1.
+         CALL Init_neighbours(knon,neighbours,knindex,pctsrf(:,is_ter))
+
+       ELSE IF (grid_type==unstructured) THEN
+ 
+         IF ( (.NOT.ALLOCATED(neighbours))) THEN
+          ALLOCATE(neighbours(knon,12), stat = error)
+          IF (error /= 0) THEN
+             abort_message='Pb allocation neighbours'
+             CALL abort_physic(modname,abort_message,1)
+          ENDIF
+         ENDIF
+         neighbours = -1.
+ 
+       ENDIF
+         
 
 !
@@ -345,10 +446,13 @@
           ENDIF
        ENDIF
-       DO igrid = 1, knon
-          ij = knindex(igrid)
-          resolution(igrid,1) = dx(ij)
-          resolution(igrid,2) = dy(ij)
-       ENDDO
-     
+       
+       IF (grid_type==regular_lonlat) THEN
+         DO igrid = 1, knon
+            ij = knindex(igrid)
+            resolution(igrid,1) = dx(ij)
+           resolution(igrid,2) = dy(ij)
+         ENDDO
+       ENDIF
+        
        ALLOCATE(coastalflow(klon), stat = error)
        IF (error /= 0) THEN
@@ -401,17 +505,60 @@
     IF (debut) THEN
        CALL Init_orchidee_index(knon,knindex,offset,ktindex)
-       CALL Get_orchidee_communicator(orch_comm,orch_omp_size,orch_omp_rank)
+       CALL Get_orchidee_communicator(orch_comm,orch_mpi_size,orch_mpi_rank, orch_omp_size,orch_omp_rank)
+
+       IF (grid_type==unstructured) THEN
+         IF (knon==0) THEN
+           begin=1
+           end=0
+         ELSE
+           begin=offset+1
+           end=offset+ktindex(knon)
+         ENDIF
+        
+         IF (orch_mpi_rank==orch_mpi_size-1 .AND. orch_omp_rank==orch_omp_size-1) end=nbp_lon*nbp_lat 
+          
+         ALLOCATE(lalo(end-begin+1,2))
+         ALLOCATE(bounds_lalo(end-begin+1,nvertex,2))
+         ALLOCATE(ind_cell(end-begin+1))
+         
+         ALLOCATE(longitude_glo(klon_glo))
+         CALL gather(longitude,longitude_glo)
+         CALL bcast(longitude_glo)
+         lalo(:,2)=longitude_glo(begin:end)*180./PI
+ 
+         ALLOCATE(latitude_glo(klon_glo))
+         CALL gather(latitude,latitude_glo)
+         CALL bcast(latitude_glo)
+         lalo(:,1)=latitude_glo(begin:end)*180./PI
+
+         ALLOCATE(boundslon_glo(klon_glo,nvertex))
+         CALL gather(boundslon,boundslon_glo)
+         CALL bcast(boundslon_glo)
+         bounds_lalo(:,:,2)=boundslon_glo(begin:end,:)*180./PI
+ 
+         ALLOCATE(boundslat_glo(klon_glo,nvertex))
+         CALL gather(boundslat,boundslat_glo)
+         CALL bcast(boundslat_glo)
+         bounds_lalo(:,:,1)=boundslat_glo(begin:end,:)*180./PI
+         
+         ALLOCATE(ind_cell_glo_glo(klon_glo))
+         CALL gather(ind_cell_glo,ind_cell_glo_glo)
+         CALL bcast(ind_cell_glo_glo)
+         ind_cell(:)=ind_cell_glo_glo(begin:end)
+         
+       ENDIF
        CALL Init_synchro_omp
        
        IF (knon > 0) THEN
 #ifdef CPP_VEGET
-         CALL Init_intersurf(nbp_lon,nbp_lat,knon,ktindex,offset,orch_omp_size,orch_omp_rank,orch_comm)
+         CALL Init_intersurf(nbp_lon,nbp_lat,knon,ktindex,offset,orch_omp_size,orch_omp_rank,orch_comm,grid=grid_type)
 #endif
        ENDIF
 
        
-       IF (knon > 0) THEN
+       IF (knon > 0) THEN 
 
 #ifdef CPP_VEGET
+
          CALL intersurf_initialize_gathered (itime+itau_phy-1, nbp_lon, nbp_lat, knon, ktindex, dtime, &
                lrestart_read, lrestart_write, lalo, contfrac, neighbours, resolution, date0, &
@@ -421,5 +568,6 @@
                evap, fluxsens, fluxlat, coastalflow, riverflow, &
                tsol_rad, tsurf_new, qsurf, albedo_out, emis_new, z0m_new, &    
-               lon_scat, lat_scat, q2m, t2m, z0h_new, nvm_orch)
+               lon_scat, lat_scat, q2m(1:knon), t2m(1:knon), z0h_new(1:knon), nvm_orch, &
+               grid=grid_type, bounds_latlon=bounds_lalo, cell_area=area, ind_cell_glo=ind_cell)
 #endif         
        ENDIF
@@ -450,5 +598,5 @@
             evap(1:knon), fluxsens(1:knon), fluxlat(1:knon), coastalflow(1:knon), riverflow(1:knon), &
             tsol_rad(1:knon), tsurf_new(1:knon), qsurf(1:knon), albedo_out(1:knon,:), emis_new(1:knon), z0m_new(1:knon), &
-            lon_scat, lat_scat, q2m, t2m, z0h_new(1:knon),&
+            lon_scat, lat_scat, q2m(1:knon), t2m(1:knon), z0h_new(1:knon),&
             veget(1:knon,:),lai(1:knon,:),height(1:knon,:),&
             coszang=yrmu0(1:knon))
@@ -525,5 +673,5 @@
 ! 
 
-  SUBROUTINE Get_orchidee_communicator(orch_comm,orch_omp_size,orch_omp_rank)
+  SUBROUTINE Get_orchidee_communicator(orch_comm, orch_mpi_size, orch_mpi_rank, orch_omp_size,orch_omp_rank)
   USE  mod_surf_para
       
@@ -533,4 +681,6 @@
 
     INTEGER,INTENT(OUT) :: orch_comm
+    INTEGER,INTENT(OUT) :: orch_mpi_size
+    INTEGER,INTENT(OUT) :: orch_mpi_rank
     INTEGER,INTENT(OUT) :: orch_omp_size
     INTEGER,INTENT(OUT) :: orch_omp_rank
@@ -552,4 +702,6 @@
 #ifdef CPP_MPI    
       CALL MPI_COMM_SPLIT(COMM_LMDZ_PHY,color,mpi_rank,orch_comm,ierr)
+      CALL MPI_COMM_SIZE(orch_comm,orch_mpi_size,ierr)
+      CALL MPI_COMM_RANK(orch_comm,orch_mpi_rank,ierr)
 #endif
     
@@ -683,3 +835,4 @@
 #endif
 #endif
+#endif
 END MODULE surf_land_orchidee_mod
