Index: LMDZ5/branches/testing/libf/phylmd/cosp/cosp_output_mod.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/cosp/cosp_output_mod.F90	(revision 2787)
+++ LMDZ5/branches/testing/libf/phylmd/cosp/cosp_output_mod.F90	(revision 2839)
@@ -20,8 +20,8 @@
       INTEGER, DIMENSION(3), SAVE  :: nhoricosp,nvert,nvertmcosp,nvertcol,nvertbze, &
                                       nvertsratio,nvertisccp,nvertp,nverttemp,nvertmisr, &
-                                      nvertReffIce,nvertReffLiq
+                                      nvertReffIce,nvertReffLiq,nverttau
       REAL, DIMENSION(3), SAVE                :: zoutm_cosp
 !$OMP THREADPRIVATE(nhoricosp, nvert,nvertmcosp,nvertcol,nvertsratio,nvertbze,nvertisccp,nvertp,zoutm_cosp,nverttemp,nvertmisr)
-!$OMP THREADPRIVATE(nvertReffIce,nvertReffLiq)
+!$OMP THREADPRIVATE(nvertReffIce,nvertReffLiq,nverttau)
       REAL, SAVE                   :: zdtimemoy_cosp
 !$OMP THREADPRIVATE(zdtimemoy_cosp) 
@@ -207,6 +207,5 @@
   SUBROUTINE cosp_output_open(Nlevlmdz, Ncolumns, presnivs, dtime, freq_cosp, &
                               ok_mensuelCOSP, ok_journeCOSP, ok_hfCOSP, ok_all_xml,  &
-                              ecrit_mth, ecrit_day, ecrit_hf, use_vgrid, vgrid)
-
+                              ecrit_mth, ecrit_day, ecrit_hf, use_vgrid, vgrid, stlidar)
 
   USE iophy
@@ -229,4 +228,5 @@
   logical                  :: ok_mensuelCOSP, ok_journeCOSP, ok_hfCOSP, use_vgrid, ok_all_xml                    
   type(cosp_vgrid)   :: vgrid   ! Information on vertical grid of stats
+  type(cosp_lidarstats) :: stlidar ! Summary statistics from lidar simulator
 
 !!! Variables locales
@@ -236,4 +236,5 @@
   real,dimension(2,SR_BINS) :: sratio_bounds
   real,dimension(SR_BINS)   ::  sratio_ax
+  real,dimension(DBZE_BINS) ::  dbze_ax
   CHARACTER(LEN=20), DIMENSION(3)  :: chfreq = (/ '1day', '1d  ', '3h  ' /)            
 
@@ -257,18 +258,19 @@
     enddo
 
-!    do ii=1,DBZE_BINS
-!     dbze_ax(i) = CFAD_ZE_MIN + CFAD_ZE_WIDTH*(ii - 0.5)
-!    enddo
-
-!   sratio_bounds(2,:)=stlidar%srbval(:) ! srbval contains the upper
+    do i=1,DBZE_BINS
+     dbze_ax(i) = CFAD_ZE_MIN + CFAD_ZE_WIDTH*(i - 0.5)
+    enddo
+ 
+
+    sratio_bounds(2,:)=stlidar%srbval(:) ! srbval contains the upper
 !                                         limits from lmd_ipsl_stats.f90
-!   sratio_bounds(1,2:SR_BINS) = stlidar%srbval(1:SR_BINS-1)
-!   sratio_bounds(1,1)         = 0.0
-!   sratio_bounds(2,SR_BINS)   = 1.e5 ! This matches with Chepfer et al., JGR,
+    sratio_bounds(1,2:SR_BINS) = stlidar%srbval(1:SR_BINS-1)
+    sratio_bounds(1,1)         = 0.0
+    sratio_bounds(2,SR_BINS)   = 1.e5 ! This matches with Chepfer et al., JGR,
 !                                    ! 2009. However, it is not consistent 
                                      ! with the upper limit in
                                      ! lmd_ipsl_stats.f90, which is
                                      ! LIDAR_UNDEF-1=998.999
-!    sratio_ax(:) = (sratio_bounds(1,:)+sratio_bounds(2,:))/2.0 
+     sratio_ax(:) = (sratio_bounds(1,:)+sratio_bounds(2,:))/2.0 
 
     cosp_outfilenames(1) = 'histmthCOSP'
@@ -331,8 +333,10 @@
    CALL wxios_add_vaxis("temp", LIDAR_NTEMP, LIDAR_PHASE_TEMP)
    CALL wxios_add_vaxis("cth16", MISR_N_CTH, MISR_CTH)
-!   CALL wxios_add_vaxis("dbze", DBZE_BINS, dbze_ax) 
-!   CALL wxios_add_vaxis("scatratio", SR_BINS, sratio_ax)
+   CALL wxios_add_vaxis("dbze", DBZE_BINS, dbze_ax) 
+   CALL wxios_add_vaxis("scatratio", SR_BINS, sratio_ax)
    CALL wxios_add_vaxis("ReffIce", numMODISReffIceBins, reffICE_binCenters)
    CALL wxios_add_vaxis("ReffLiq", numMODISReffLiqBins, reffLIQ_binCenters)
+   print*,'reffICE_binCenters=',reffICE_binCenters
+   CALL wxios_add_vaxis("tau", 7, ISCCP_TAU)
 
 #endif
@@ -384,7 +388,9 @@
                     nvertReffLiq(iff))
 
-!      CALL histvert(cosp_nidfiles(iff),"dbze","equivalent_reflectivity_factor","dBZ",DBZE_BINS,dbze_ax,nvertbze(iff))
+      CALL histvert(cosp_nidfiles(iff),"dbze","equivalent_reflectivity_factor","dBZ",DBZE_BINS,dbze_ax,nvertbze(iff))
      
-!      CALL histvert(cosp_nidfiles(iff),"scatratio","backscattering_ratio","1",SR_BINS,sratio_ax,nvertsratio(iff))
+      CALL histvert(cosp_nidfiles(iff),"scatratio","backscattering_ratio","1",SR_BINS,sratio_ax,nvertsratio(iff))
+
+      CALL histvert(cosp_nidfiles(iff),"tau","cloud optical depth","1",7,ISCCP_TAU,nverttau(iff)) 
      
 !!! Valeur indefinie en cas IOIPSL
Index: LMDZ5/branches/testing/libf/phylmd/cosp/cosp_output_write_mod.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/cosp/cosp_output_write_mod.F90	(revision 2787)
+++ LMDZ5/branches/testing/libf/phylmd/cosp/cosp_output_write_mod.F90	(revision 2839)
@@ -17,5 +17,5 @@
    CONTAINS
 
-  SUBROUTINE cosp_output_write(Nlevlmdz, Npoints, Ncolumns, itap, dtime, freq_COSP, &
+  SUBROUTINE cosp_output_write(Nlevlmdz, Npoints, Ncolumns, itap, dtime, freq_COSP, missing_cosp, &
                                cfg, gbx, vgrid, sglidar, sgradar, stlidar, stradar, &
                                isccp, misr, modis)
@@ -32,5 +32,5 @@
 !!! Variables d'entree
   integer               :: itap, Nlevlmdz, Ncolumns, Npoints
-  real                  :: freq_COSP, dtime
+  real                  :: freq_COSP, dtime, missing_val, missing_cosp
   type(cosp_config)     :: cfg     ! Control outputs
   type(cosp_gridbox)    :: gbx     ! Gridbox information. Input for COSP
@@ -49,4 +49,11 @@
   integer               :: itau_wcosp
   real, dimension(Npoints,PARASOL_NREFL) :: parasolcrefl, Ncref
+
+
+#ifdef CPP_XIOS
+  missing_val=missing_cosp
+#else
+  missing_val=0.
+#endif
 
   Nlevout = vgrid%Nlvgrid
@@ -95,5 +102,5 @@
      do ip = 1,Npoints
      if(stlidar%lidarcld(ip,k).eq.R_UNDEF)then
-      stlidar%lidarcld(ip,k)=Cosp_fill_value
+      stlidar%lidarcld(ip,k)=missing_val
      endif
      enddo
@@ -102,5 +109,5 @@
       do ip = 1,Npoints
        if(stlidar%cfad_sr(ip,ii,k).eq.R_UNDEF)then
-        stlidar%cfad_sr(ip,ii,k)=Cosp_fill_value
+        stlidar%cfad_sr(ip,ii,k)=missing_val
        endif
       enddo
@@ -111,10 +118,10 @@
    do k = 1,Nlevlmdz
      if(sglidar%beta_mol(ip,k).eq.R_UNDEF)then
-      sglidar%beta_mol(ip,k)=Cosp_fill_value
+      sglidar%beta_mol(ip,k)=missing_val
      endif
 
      do ii= 1,Ncolumns
        if(sglidar%beta_tot(ip,ii,k).eq.R_UNDEF)then
-        sglidar%beta_tot(ip,ii,k)=Cosp_fill_value
+        sglidar%beta_tot(ip,ii,k)=missing_val
        endif
      enddo
@@ -126,5 +133,5 @@
     do ip = 1,Npoints
      if(stlidar%cldlayer(ip,k).eq.R_UNDEF)then
-      stlidar%cldlayer(ip,k)=Cosp_fill_value
+       stlidar%cldlayer(ip,k)=missing_val
      endif
     enddo
@@ -133,9 +140,9 @@
 ! AI 11 / 2015
 
-   where(stlidar%parasolrefl == R_UNDEF) stlidar%parasolrefl = 0.0
-   where(stlidar%lidarcldtmp == R_UNDEF) stlidar%lidarcldtmp = 0.0
-   where(stlidar%cldlayerphase == R_UNDEF) stlidar%cldlayerphase = 0.0
-   where(stlidar%lidarcldphase == R_UNDEF) stlidar%lidarcldphase = 0.0
-   where(stlidar%lidarcldtmp == R_UNDEF) stlidar%lidarcldtmp = 0.0
+   where(stlidar%parasolrefl == R_UNDEF) stlidar%parasolrefl = missing_val
+   where(stlidar%lidarcldtmp == R_UNDEF) stlidar%lidarcldtmp = missing_val
+   where(stlidar%cldlayerphase == R_UNDEF) stlidar%cldlayerphase = missing_val
+   where(stlidar%lidarcldphase == R_UNDEF) stlidar%lidarcldphase = missing_val
+   where(stlidar%lidarcldtmp == R_UNDEF) stlidar%lidarcldtmp = missing_val
    
 
@@ -169,9 +176,12 @@
    CALL histwrite3d_cosp(o_clcalipsotmpun,stlidar%lidarcldtmp(:,:,4),nverttemp)
 
+#ifdef CPP_XIOS
+   CALL histwrite4d_cosp(o_cfad_lidarsr532,stlidar%cfad_sr)
+#else
    do icl=1,SR_BINS
       CALL histwrite3d_cosp(o_cfad_lidarsr532,stlidar%cfad_sr(:,icl,:),nvert,icl)
    enddo
-
-   CALL histwrite3d_cosp(o_parasol_refl,stlidar%parasolrefl,nvertp)
+     CALL histwrite3d_cosp(o_parasol_refl,stlidar%parasolrefl,nvertp)
+#endif
 
    do k=1,PARASOL_NREFL
@@ -182,5 +192,5 @@
         Ncref(ip,k) = 1.
      else
-        parasolcrefl(ip,k)=0.
+        parasolcrefl(ip,k)=missing_val
         Ncref(ip,k) = 0.
      endif
@@ -190,11 +200,21 @@
    CALL histwrite3d_cosp(o_parasol_crefl,parasolcrefl,nvertp)
 
+#ifdef CPP_XIOS
+   CALL histwrite4d_cosp(o_atb532,sglidar%beta_tot)
+#else
    do icl=1,Ncolumns 
       CALL histwrite3d_cosp(o_atb532,sglidar%beta_tot(:,icl,:),nvertmcosp,icl)
    enddo
+#endif
+
    CALL histwrite3d_cosp(o_beta_mol532,sglidar%beta_mol,nvertmcosp) 
  endif !Lidar
 
  if (cfg%Lradar_sim) then
+
+#ifdef CPP_XIOS
+   CALL histwrite4d_cosp(o_dbze94,sgradar%Ze_tot)
+   CALL histwrite4d_cosp(o_cfadDbze94,stradar%cfad_ze)
+#else
    do icl=1,Ncolumns
       CALL histwrite3d_cosp(o_dbze94,sgradar%Ze_tot(:,icl,:),nvertmcosp,icl)
@@ -203,12 +223,13 @@
     CALL histwrite3d_cosp(o_cfadDbze94,stradar%cfad_ze(:,icl,:),nvert,icl)
    enddo
+#endif
  endif
 
  if (cfg%Llidar_sim .and. cfg%Lradar_sim) then
    where(stradar%lidar_only_freq_cloud == R_UNDEF) &
-                           stradar%lidar_only_freq_cloud = 0.0
+                           stradar%lidar_only_freq_cloud = missing_val
    CALL histwrite3d_cosp(o_clcalipso2,stradar%lidar_only_freq_cloud,nvert)
    where(stradar%radar_lidar_tcc == R_UNDEF) &
-                           stradar%radar_lidar_tcc = 0.0 
+                           stradar%radar_lidar_tcc = missing_val
    CALL histwrite2d_cosp(o_cltlidarradar,stradar%radar_lidar_tcc)
  endif
@@ -219,20 +240,20 @@
    do ip = 1,Npoints
     if(isccp%totalcldarea(ip).eq.R_UNDEF)then
-      isccp%totalcldarea(ip)=Cosp_fill_value
+      isccp%totalcldarea(ip)=missing_val
     endif
     if(isccp%meanptop(ip).eq.R_UNDEF)then
-      isccp%meanptop(ip)=Cosp_fill_value
+      isccp%meanptop(ip)=missing_val
     endif
     if(isccp%meantaucld(ip).eq.R_UNDEF)then
-      isccp%meantaucld(ip)=Cosp_fill_value
+      isccp%meantaucld(ip)=missing_val
     endif
     if(isccp%meanalbedocld(ip).eq.R_UNDEF)then
-      isccp%meanalbedocld(ip)=Cosp_fill_value
+     isccp%meanalbedocld(ip)=missing_val
     endif
     if(isccp%meantb(ip).eq.R_UNDEF)then
-      isccp%meantb(ip)=Cosp_fill_value
+     isccp%meantb(ip)=missing_val
     endif
     if(isccp%meantbclr(ip).eq.R_UNDEF)then
-      isccp%meantbclr(ip)=Cosp_fill_value
+       isccp%meantbclr(ip)=missing_val
     endif
 
@@ -240,5 +261,5 @@
      do ii=1,7
      if(isccp%fq_isccp(ip,ii,k).eq.R_UNDEF)then
-      isccp%fq_isccp(ip,ii,k)=Cosp_fill_value
+      isccp%fq_isccp(ip,ii,k)=missing_val
      endif
      enddo
@@ -247,5 +268,5 @@
     do ii=1,Ncolumns
      if(isccp%boxtau(ip,ii).eq.R_UNDEF)then
-       isccp%boxtau(ip,ii)=Cosp_fill_value
+       isccp%boxtau(ip,ii)=missing_val
      endif
     enddo
@@ -253,5 +274,5 @@
     do ii=1,Ncolumns
      if(isccp%boxptop(ip,ii).eq.R_UNDEF)then
-       isccp%boxptop(ip,ii)=Cosp_fill_value
+      isccp%boxptop(ip,ii)=missing_val
      endif
     enddo
@@ -259,7 +280,11 @@
 
    CALL histwrite2d_cosp(o_sunlit,gbx%sunlit)
+#ifdef CPP_XIOS
+   CALL histwrite4d_cosp(o_clisccp2,isccp%fq_isccp)
+#else
    do icl=1,7
    CALL histwrite3d_cosp(o_clisccp2,isccp%fq_isccp(:,icl,:),nvertisccp,icl) 
    enddo
+#endif
    CALL histwrite3d_cosp(o_boxtauisccp,isccp%boxtau,nvertcol)
    CALL histwrite3d_cosp(o_boxptopisccp,isccp%boxptop,nvertcol) 
@@ -278,5 +303,5 @@
        do k=1,MISR_N_CTH
         if(misr%fq_MISR(ip,ii,k).eq.R_UNDEF)then
-              misr%fq_MISR(ip,ii,k)=Cosp_fill_value
+            misr%fq_MISR(ip,ii,k)=missing_val
         endif
        enddo
@@ -284,7 +309,11 @@
    enddo
 
+#ifdef CPP_XIOS
+   CALL histwrite4d_cosp(o_clMISR,misr%fq_MISR)
+#else
    do icl=1,7 
       CALL histwrite3d_cosp(o_clMISR,misr%fq_MISR(:,icl,:),nvertmisr,icl)
    enddo
+#endif
  endif
 
@@ -294,47 +323,57 @@
   do ip=1,Npoints
     if(modis%Cloud_Fraction_Low_Mean(ip).eq.R_UNDEF)then
-       modis%Cloud_Fraction_Low_Mean(ip)=Cosp_fill_value
+      modis%Cloud_Fraction_Low_Mean(ip)=missing_val
     endif
     if(modis%Cloud_Fraction_High_Mean(ip).eq.R_UNDEF)then
-       modis%Cloud_Fraction_High_Mean(ip)=Cosp_fill_value
+       modis%Cloud_Fraction_High_Mean(ip)=missing_val 
     endif
     if(modis%Cloud_Fraction_Mid_Mean(ip).eq.R_UNDEF)then
-       modis%Cloud_Fraction_Mid_Mean(ip)=Cosp_fill_value
+       modis%Cloud_Fraction_Mid_Mean(ip)=missing_val
     endif
     if(modis%Cloud_Fraction_Total_Mean(ip).eq.R_UNDEF)then
-       modis%Cloud_Fraction_Total_Mean(ip)=Cosp_fill_value
+       modis%Cloud_Fraction_Total_Mean(ip)=missing_val
     endif
     if(modis%Cloud_Fraction_Water_Mean(ip).eq.R_UNDEF)then
-       modis%Cloud_Fraction_Water_Mean(ip)=Cosp_fill_value
+       modis%Cloud_Fraction_Water_Mean(ip)=missing_val
     endif
     if(modis%Cloud_Fraction_Ice_Mean(ip).eq.R_UNDEF)then
-       modis%Cloud_Fraction_Ice_Mean(ip)=Cosp_fill_value
+       modis%Cloud_Fraction_Ice_Mean(ip)=missing_val
     endif
     if(modis%Optical_Thickness_Total_Mean(ip).eq.R_UNDEF)then
-       modis%Optical_Thickness_Total_Mean(ip)=Cosp_fill_value
+       modis%Optical_Thickness_Total_Mean(ip)=missing_val
     endif
     if(modis%Optical_Thickness_Water_Mean(ip).eq.R_UNDEF)then
-       modis%Optical_Thickness_Water_Mean(ip)=Cosp_fill_value
+       modis%Optical_Thickness_Water_Mean(ip)=missing_val
     endif
     if(modis%Optical_Thickness_Ice_Mean(ip).eq.R_UNDEF)then
-       modis%Optical_Thickness_Ice_Mean(ip)=Cosp_fill_value
+       modis%Optical_Thickness_Ice_Mean(ip)=missing_val
     endif
     if(modis%Cloud_Particle_Size_Water_Mean(ip).eq.R_UNDEF)then
-       modis%Cloud_Particle_Size_Water_Mean(ip)=Cosp_fill_value
+       modis%Cloud_Particle_Size_Water_Mean(ip)=missing_val
     endif
     if(modis%Cloud_Particle_Size_Ice_Mean(ip).eq.R_UNDEF)then
-       modis%Cloud_Particle_Size_Ice_Mean(ip)=Cosp_fill_value
+       modis%Cloud_Particle_Size_Ice_Mean(ip)=missing_val
     endif
     if(modis%Cloud_Top_Pressure_Total_Mean(ip).eq.R_UNDEF)then
-       modis%Cloud_Top_Pressure_Total_Mean(ip)=Cosp_fill_value
+       modis%Cloud_Top_Pressure_Total_Mean(ip)=missing_val
     endif
     if(modis%Liquid_Water_Path_Mean(ip).eq.R_UNDEF)then
-       modis%Liquid_Water_Path_Mean(ip)=Cosp_fill_value
+       modis%Liquid_Water_Path_Mean(ip)=missing_val
     endif
     if(modis%Ice_Water_Path_Mean(ip).eq.R_UNDEF)then
-       modis%Ice_Water_Path_Mean(ip)=Cosp_fill_value
+       modis%Ice_Water_Path_Mean(ip)=missing_val
     endif
 
   enddo
+
+    where(modis%Optical_Thickness_Total_LogMean == R_UNDEF) &
+          modis%Optical_Thickness_Total_LogMean = missing_val
+           
+ 
+    where(modis%Optical_Thickness_Water_LogMean == R_UNDEF) &
+          modis%Optical_Thickness_Water_LogMean = missing_val
+
+    where(modis%Optical_Thickness_Ice_LogMean == R_UNDEF) &
+          modis%Optical_Thickness_Ice_LogMean = missing_val
     
    CALL histwrite2d_cosp(o_cllmodis,modis%Cloud_Fraction_Low_Mean)
@@ -360,5 +399,5 @@
        do k=1,7
        if(modis%Optical_Thickness_vs_Cloud_Top_Pressure(ip,ii,k).eq.R_UNDEF)then
-          modis%Optical_Thickness_vs_Cloud_Top_Pressure(ip,ii,k)=0.
+          modis%Optical_Thickness_vs_Cloud_Top_Pressure(ip,ii,k)=missing_val
         endif
        enddo
@@ -366,15 +405,25 @@
     enddo 
 
+#ifdef CPP_XIOS
+   CALL histwrite4d_cosp(o_clmodis,modis%Optical_Thickness_vs_Cloud_Top_Pressure)
+#else
    do icl=1,7
    CALL histwrite3d_cosp(o_clmodis, &
      modis%Optical_Thickness_vs_Cloud_Top_Pressure(:,icl,:),nvertisccp,icl)           
    enddo
+#endif
 
     where(modis%Optical_Thickness_vs_ReffIce == R_UNDEF) &
-          modis%Optical_Thickness_vs_ReffIce = Cosp_fill_value
+          modis%Optical_Thickness_vs_ReffIce = missing_val
 
     where(modis%Optical_Thickness_vs_ReffLiq == R_UNDEF) &
-          modis%Optical_Thickness_vs_ReffLiq = Cosp_fill_value
-
+          modis%Optical_Thickness_vs_ReffLiq = missing_val
+
+#ifdef CPP_XIOS
+!    print*,'dimension de crimodis=',size(modis%Optical_Thickness_vs_ReffIce,2),&
+!                                    size(modis%Optical_Thickness_vs_ReffIce,3)
+    CALL histwrite4d_cosp(o_crimodis,modis%Optical_Thickness_vs_ReffIce)
+    CALL histwrite4d_cosp(o_crlmodis,modis%Optical_Thickness_vs_ReffLiq)
+#else
     do icl=1,7
    CALL histwrite3d_cosp(o_crimodis, &
@@ -383,4 +432,5 @@
      modis%Optical_Thickness_vs_ReffLiq(:,icl,:),nvertReffLiq,icl)
     enddo
+#endif
  endif
 
@@ -784,4 +834,58 @@
   END SUBROUTINE histwrite3d_cosp
 
+! ug NOUVELLE VERSION DES WRITE AVEC LA BOUCLE DO RENTREE
+! AI sept 2013
+  SUBROUTINE histwrite4d_cosp(var, field)
+  USE dimphy
+  USE mod_phys_lmdz_para
+  USE ioipsl
+  use iophy
+  USE mod_grid_phy_lmdz, ONLY: nbp_lon
+  USE print_control_mod, ONLY: lunout,prt_level
+
+#ifdef CPP_XIOS
+  USE xios, only: xios_send_field
+#endif
+
+
+  IMPLICIT NONE
+  INCLUDE 'clesphys.h'
+
+    TYPE(ctrl_outcosp), INTENT(IN)    :: var
+    REAL, DIMENSION(:,:,:), INTENT(IN)  :: field ! --> field(klon,:)
+
+    INTEGER :: iff, k
+
+    REAL,DIMENSION(klon_mpi,SIZE(field,2),SIZE(field,3)) :: buffer_omp
+    REAL :: field4d(nbp_lon,jj_nb,SIZE(field,2),SIZE(field,3))
+    INTEGER :: ip, n, nlev, nlev2
+    INTEGER, ALLOCATABLE, DIMENSION(:) :: index4d
+    CHARACTER(LEN=20) ::  nomi, nom
+
+  IF (prt_level >= 9) write(lunout,*)'Begin histrwrite4d ',var%name
+
+  IF(cosp_varsdefined) THEN
+    !Et sinon on.... Ã©crit
+    IF (SIZE(field,1)/=klon) &
+   CALL abort_physic('iophy::histwrite3d','Field first DIMENSION not equal to klon',1)            
+
+    nlev=SIZE(field,2)
+    nlev2=SIZE(field,3)
+    CALL Gather_omp(field,buffer_omp)
+!$OMP MASTER
+    CALL grid1Dto2D_mpi(buffer_omp,field4d)
+
+#ifdef CPP_XIOS
+    IF (ok_all_xml) THEN
+     CALL xios_send_field(var%name, Field4d(:,:,1:nlev,1:nlev2))
+     IF (prt_level >= 1) WRITE(lunout,*)'xios_send_field ',var%name
+    ENDIF
+#endif
+
+!$OMP END MASTER   
+  ENDIF ! vars_defined
+  IF (prt_level >= 9) write(lunout,*)'End histrwrite4d_cosp ',nom
+  END SUBROUTINE histwrite4d_cosp
+
   SUBROUTINE conf_cospoutputs(nam_var,cles_var)
 !!! Lecture des noms et cles de sortie des variables dans config.def
Index: LMDZ5/branches/testing/libf/phylmd/cosp/modis_simulator.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/cosp/modis_simulator.F90	(revision 2787)
+++ LMDZ5/branches/testing/libf/phylmd/cosp/modis_simulator.F90	(revision 2839)
@@ -532,11 +532,13 @@
 
     ! ########################################################################################
-    ! Normalize pixel counts to fraction. The first three cloud fractions have been set to -1 
-    ! in cloud-free areas, so set those places to 0.
-    ! ########################################################################################
-    Cloud_Fraction_High_Mean(1:nPoints)  = Cloud_Fraction_High_Mean(1:nPoints) /nSubcols
-    Cloud_Fraction_Mid_Mean(1:nPoints)   = Cloud_Fraction_Mid_Mean(1:nPoints)  /nSubcols
-    Cloud_Fraction_Low_Mean(1:nPoints)   = Cloud_Fraction_Low_Mean(1:nPoints)  /nSubcols
-
+    ! Normalize pixel counts to fraction.
+    ! ########################################################################################
+    Cloud_Fraction_High_Mean(1:nPoints)  = Cloud_Fraction_High_Mean(1:nPoints)  /nSubcols
+    Cloud_Fraction_Mid_Mean(1:nPoints)   = Cloud_Fraction_Mid_Mean(1:nPoints)   /nSubcols
+    Cloud_Fraction_Low_Mean(1:nPoints)   = Cloud_Fraction_Low_Mean(1:nPoints)   /nSubcols
+    Cloud_Fraction_Total_Mean(1:nPoints) = Cloud_Fraction_Total_Mean(1:nPoints) /nSubcols
+    Cloud_Fraction_Ice_Mean(1:nPoints)   = Cloud_Fraction_Ice_Mean(1:nPoints)   /nSubcols
+    Cloud_Fraction_Water_Mean(1:nPoints) = Cloud_Fraction_Water_Mean(1:nPoints) /nSubcols
+    
     ! ########################################################################################
     ! Set clear-scenes to undefined
Index: LMDZ5/branches/testing/libf/phylmd/cosp/phys_cosp.F90
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/cosp/phys_cosp.F90	(revision 2787)
+++ LMDZ5/branches/testing/libf/phylmd/cosp/phys_cosp.F90	(revision 2839)
@@ -7,5 +7,5 @@
   subroutine phys_cosp( itap,dtime,freq_cosp, &
                         ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
-                        ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, &
+                        ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
                         Nptslmdz,Nlevlmdz,lon,lat, presnivs,overlaplmdz,sunlit, &
                         ref_liq,ref_ice,fracTerLic,u_wind,v_wind,phis,phi,ph,p,skt,t, &
@@ -130,5 +130,5 @@
 ! Declaration necessaires pour les sorties IOIPSL
   integer :: ii
-  real    :: ecrit_day,ecrit_hf,ecrit_mth
+  real    :: ecrit_day,ecrit_hf,ecrit_mth, missing_val
   logical :: ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, ok_all_xml
 
@@ -192,5 +192,7 @@
 ! Allocate memory for gridbox type
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-        print *, 'Allocating memory for gridbox type...'
+! AI mars 2017
+!        print *, 'Allocating memory for gridbox type...'
+
 !! AI
 !        call construct_cosp_gridbox(dble(itap),radar_freq,surface_radar,use_mie_tables,use_gas_abs,do_ray,melt_lay,k2, &
@@ -215,5 +217,5 @@
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
-        print *, 'Populating input structure...'
+!        print *, 'Populating input structure...'
         gbx%longitude = lon
         gbx%latitude = lat
@@ -299,5 +301,5 @@
         ! Define new vertical grid
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-        print *, 'Defining new vertical grid...'
+!        print *, 'Defining new vertical grid...'
         call construct_cosp_vgrid(gbx,Nlr,use_vgrid,csat_vgrid,vgrid)
 
@@ -305,5 +307,5 @@
        ! Allocate memory for other types
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-        print *, 'Allocating memory for other types...'
+!        print *, 'Allocating memory for other types...'
         call construct_cosp_subgrid(Npoints, Ncolumns, Nlevels, sgx)
         call construct_cosp_sgradar(cfg,Npoints,Ncolumns,Nlevels,N_HYDRO,sgradar)
@@ -325,5 +327,5 @@
         call cosp_output_open(Nlevlmdz, Ncolumns, presnivs, dtime, freq_cosp, &
                               ok_mensuelCOSP, ok_journeCOSP, ok_hfCOSP, ok_all_xml, &
-                              ecrit_mth, ecrit_day, ecrit_hf, use_vgrid, vgrid)
+                              ecrit_mth, ecrit_day, ecrit_hf, use_vgrid, vgrid, stlidar)
       !$OMP END MASTER
       !$OMP BARRIER
@@ -333,5 +335,5 @@
         ! Call simulator
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-        print *, 'Calling simulator...'
+!        print *, 'Calling simulator...'
 !! AI 
 !        call cosp(overlap,Ncolumns,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,stradar,stlidar)
@@ -347,6 +349,6 @@
 !!!!!!!!!!!!!!!!!! Ecreture des sorties Cosp !!!!!!!!!!!!!!r!!!!!!:!!!!!
 
-       print *, 'Calling write output'
-        call cosp_output_write(Nlevlmdz, Npoints, Ncolumns, itap, dtime, freq_COSP, &
+!       print *, 'Calling write output'
+        call cosp_output_write(Nlevlmdz, Npoints, Ncolumns, itap, dtime, freq_COSP, missing_val, &
                                cfg, gbx, vgrid, sglidar, sgradar, stlidar, stradar, & 
                                isccp, misr, modis)
@@ -355,5 +357,5 @@
         ! Deallocate memory in derived types
 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-        print *, 'Deallocating memory...'
+!        print *, 'Deallocating memory...'
         call free_cosp_gridbox(gbx)
         call free_cosp_subgrid(sgx)
