diff --git a/src/MNH/mean_field.f90 b/src/MNH/mean_field.f90
index 573046c23188f3d2f7edd9804864543ceebb933b..a8a5cd98c508dc08f14cfb3f405ae686c3d5a230 100644
--- a/src/MNH/mean_field.f90
+++ b/src/MNH/mean_field.f90
@@ -67,6 +67,8 @@ USE MODD_CST
 USE MODD_PASPOL
 USE MODE_THERMO
 USE MODI_SHUMAN
+! To play with MPI
+USE MODD_ARGSLIST_ll, ONLY: LIST_ll
 !
 ! * EOL
 USE MODD_EOL_MAIN, ONLY: LMAIN_EOL, CMETH_EOL, NMODEL_EOL
@@ -77,6 +79,8 @@ USE MODD_EOL_ADR, ONLY: XFAERO_BLEQ_RA_GLB, XFAERO_BLEQ_RA_SUM
 USE MODD_EOL_ADR, ONLY: XAOA_BLEQ_GLB, XAOA_BLEQ_SUM
 USE MODD_EOL_ALM, ONLY: XFAERO_RE_SUM, XFAERO_RE_GLB
 !
+USE MODI_SHUMAN,  ONLY: MXF, MYF, MZF
+!
 USE MODE_MODELN_HANDLER
 USE MODI_UPDATE_WELFORD
 !
@@ -97,7 +101,6 @@ INTEGER           :: IIU,IJU,IKU,IIB,IJB,IKB,IIE,IJE,IKE ! Arrays bounds
 INTEGER           :: JI,JJ,JK   ! Loop indexes
 REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZUMEAN_OLD
 REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZWMEAN_OLD
-REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZTHMEAN_OLD
 !
 REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZRT
 REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZQSAT_W
@@ -115,6 +118,9 @@ REAL, DIMENSION(:,:),   ALLOCATABLE :: ZRH_P_MAXCOL
 !
 INTEGER :: IMI !Current model index
 !
+! For communication
+TYPE(LIST_ll), POINTER :: TZFIELDS_Cov_ll  ! Field list of Wind for exchange beofre covariances
+INTEGER                :: IINFO            ! Info integer
 !
 !-----------------------------------------------------------------------
 !
@@ -228,10 +234,8 @@ END IF
    IF (LCOV_FIELD) THEN
       ALLOCATE( ZUMEAN_OLD (IIU, IJU, IKU) )
       ALLOCATE( ZWMEAN_OLD (IIU, IJU, IKU) )
-      ALLOCATE( ZTHMEAN_OLD(IIU, IJU, IKU) )
       ZUMEAN_OLD (:,:,:) = XUM_MEAN (:,:,:)
       ZWMEAN_OLD (:,:,:) = XWM_MEAN (:,:,:)
-      ZTHMEAN_OLD(:,:,:) = XTHM_MEAN(:,:,:)
    END IF
 !  Welford method for variables whom we compute variances
 !
@@ -244,12 +248,24 @@ END IF
 !
 !  Welford method for covariances
 !
-   IF (LCOV_FIELD) THEN
-     XUV_MEAN  = XUV_MEAN + (PUT-ZUMEAN_OLD)*(PVT-XVM_MEAN)
-     XUW_MEAN  = XUW_MEAN + (PUT-ZUMEAN_OLD)*(PWT-XWM_MEAN)
-     XVW_MEAN  = XVW_MEAN + (PWT-ZWMEAN_OLD)*(PVT-XVM_MEAN)
-     XWTH_MEAN = XVW_MEAN + (PWT-ZWMEAN_OLD)*(PTHT-XTHM_MEAN)
-     DEALLOCATE( ZUMEAN_OLD, ZWMEAN_OLD, ZTHMEAN_OLD )
+   IF (LCOV_FIELD) THEN 
+! 
+    NULLIFY(TZFIELDS_Cov_ll)
+    CALL ADD3DFIELD_ll( TZFIELDS_Cov_ll,PUT,        'MEAN_FIELD::PUT')
+    CALL ADD3DFIELD_ll( TZFIELDS_Cov_ll,PWT,        'MEAN_FIELD::PWT')
+    CALL ADD3DFIELD_ll( TZFIELDS_Cov_ll,PVT,        'MEAN_FIELD::PVT')
+    CALL ADD3DFIELD_ll( TZFIELDS_Cov_ll,ZUMEAN_OLD, 'MEAN_FIELD::ZUMEAN_OLD')
+    CALL ADD3DFIELD_ll( TZFIELDS_Cov_ll,ZWMEAN_OLD, 'MEAN_FIELD::ZWMEAN_OLD')
+    CALL ADD3DFIELD_ll( TZFIELDS_Cov_ll,XVM_MEAN,   'MEAN_FIELD::XVM_MEAN')
+    CALL ADD3DFIELD_ll( TZFIELDS_Cov_ll,XWM_MEAN,   'MEAN_FIELD::XWM_MEAN')
+    CALL UPDATE_HALO_ll(TZFIELDS_Cov_ll,IINFO)
+    CALL CLEANLIST_ll(  TZFIELDS_Cov_ll)
+!
+    XUV_MEAN  = XUV_MEAN  + (MXF(PUT)-MXF(ZUMEAN_OLD))*(MYF(PVT)-MYF(XVM_MEAN))
+    XUW_MEAN  = XUW_MEAN  + (MXF(PUT)-MXF(ZUMEAN_OLD))*(MZF(PWT)-MZF(XWM_MEAN))
+    XVW_MEAN  = XVW_MEAN  + (MZF(PWT)-MZF(ZWMEAN_OLD))*(MYF(PVT)-MYF(XVM_MEAN))
+    XWTH_MEAN = XWTH_MEAN + (MZF(PWT)-MZF(ZWMEAN_OLD))*(PTHT-XTHM_MEAN)
+    DEALLOCATE( ZUMEAN_OLD, ZWMEAN_OLD )
    END IF
 !
 !-----------------------------------------------------------------------
diff --git a/src/MNH/write_lfin.f90 b/src/MNH/write_lfin.f90
index 43e3f1e22a929e658dfeab027eb2e00883750b7f..9a023ff617c2f6358b4d152a9b9ec3971ab3d7a5 100644
--- a/src/MNH/write_lfin.f90
+++ b/src/MNH/write_lfin.f90
@@ -757,38 +757,6 @@ IF (MEAN_COUNT /= 0) THEN
   TZFIELD%CUNITS     = 'm s-1'
   TZFIELD%CCOMMENT   = 'X_Y_Z_U component of max wind'
   CALL IO_Field_write(TPFILE,TZFIELD,XUM_MAX)
-!
-  IF (LCOV_FIELD) THEN
-    !
-    TZFIELD%CMNHNAME   = 'UVME'
-    TZFIELD%CLONGNAME  = 'UVME'
-    TZFIELD%CUNITS     = 'm2 s-2'
-    TZFIELD%CCOMMENT   = 'X_Y_Z_UV component of mean wind variance'
-    ZWORK3D = XUV_MEAN/MEAN_COUNT
-    CALL IO_Field_write(TPFILE,TZFIELD,ZWORK3D)
-    !
-    TZFIELD%CMNHNAME   = 'UWME'
-    TZFIELD%CLONGNAME  = 'UWME'
-    TZFIELD%CUNITS     = 'm2 s-2'
-    TZFIELD%CCOMMENT   = 'X_Y_Z_UW component of mean wind variance'
-    ZWORK3D = XUW_MEAN/MEAN_COUNT
-    CALL IO_Field_write(TPFILE,TZFIELD,ZWORK3D)
-    !
-    TZFIELD%CMNHNAME   = 'VWME'
-    TZFIELD%CLONGNAME  = 'VWME'
-    TZFIELD%CUNITS     = 'm2 s-2'
-    TZFIELD%CCOMMENT   = 'X_Y_Z_VW component of mean wind variance'
-    ZWORK3D = XVW_MEAN/MEAN_COUNT
-    CALL IO_Field_write(TPFILE,TZFIELD,ZWORK3D)
-    !
-    TZFIELD%CMNHNAME   = 'WTHME'
-    TZFIELD%CLONGNAME  = 'WTHME'
-    TZFIELD%CUNITS     = 'm2 s-2'
-    TZFIELD%CCOMMENT   = 'X_Y_Z_WTH component of mean wind variance'
-    ZWORK3D = XWTH_MEAN/MEAN_COUNT
-    CALL IO_Field_write(TPFILE,TZFIELD,ZWORK3D)
-    !
-  END IF
 !
   TZFIELD = TFIELDMETADATA(                          &
     CMNHNAME   = 'generic for mean_count variables', & !Temporary name to ease identification
@@ -1028,6 +996,38 @@ IF (MEAN_COUNT /= 0) THEN
     TZFIELD%CCOMMENT   = 'X_Y_Z_max kinetic energy'
     CALL IO_Field_write(TPFILE,TZFIELD,XTKEM_MAX)
   END IF
+!
+  IF (LCOV_FIELD) THEN
+    !
+    TZFIELD%CMNHNAME   = 'UVME'
+    TZFIELD%CLONGNAME  = 'UVME'
+    TZFIELD%CUNITS     = 'm2 s-2'
+    TZFIELD%CCOMMENT   = 'X_Y_Z_UV component of mean wind variance'
+    ZWORK3D = XUV_MEAN/MEAN_COUNT
+    CALL IO_Field_write(TPFILE,TZFIELD,ZWORK3D)
+    !
+    TZFIELD%CMNHNAME   = 'UWME'
+    TZFIELD%CLONGNAME  = 'UWME'
+    TZFIELD%CUNITS     = 'm2 s-2'
+    TZFIELD%CCOMMENT   = 'X_Y_Z_UW component of mean wind variance'
+    ZWORK3D = XUW_MEAN/MEAN_COUNT
+    CALL IO_Field_write(TPFILE,TZFIELD,ZWORK3D)
+    !
+    TZFIELD%CMNHNAME   = 'VWME'
+    TZFIELD%CLONGNAME  = 'VWME'
+    TZFIELD%CUNITS     = 'm2 s-2'
+    TZFIELD%CCOMMENT   = 'X_Y_Z_VW component of mean wind variance'
+    ZWORK3D = XVW_MEAN/MEAN_COUNT
+    CALL IO_Field_write(TPFILE,TZFIELD,ZWORK3D)
+    !
+    TZFIELD%CMNHNAME   = 'WTHME'
+    TZFIELD%CLONGNAME  = 'WTHME'
+    TZFIELD%CUNITS     = 'm2 s-2'
+    TZFIELD%CCOMMENT   = 'X_Y_Z_WTH component of mean wind variance'
+    ZWORK3D = XWTH_MEAN/MEAN_COUNT
+    CALL IO_Field_write(TPFILE,TZFIELD,ZWORK3D)
+    !
+  END IF
 !
 END IF
 !