Index: /trunk/LMDZ.COMMON/libf/dyn3d/caladvtrac.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d/caladvtrac.F	(revision 2125)
+++ /trunk/LMDZ.COMMON/libf/dyn3d/caladvtrac.F	(revision 2126)
@@ -9,7 +9,8 @@
 c
       USE infotrac, ONLY : nqtot
-      USE control_mod, ONLY : iapp_tracvl,planet_type
+      USE control_mod, ONLY : iapp_tracvl,planet_type,
+     &                        force_conserv_tracer
       USE comconst_mod, ONLY: dtvr
- 
+      USE planetary_operations, ONLY: planetary_tracer_amount_from_mass
       IMPLICIT NONE
 c
@@ -47,5 +48,14 @@
       INTEGER ij,l, iq, iapptrac
       REAL finmasse(ip1jmp1,llm), dtvrtrac
+      REAL :: totaltracer_old(nqtot),totaltracer_new(nqtot)
+      REAL :: ratio
 
+! Ehouarn : try to fix tracer conservation issues:
+      if (force_conserv_tracer) then
+        do iq=1,nqtot
+          call planetary_tracer_amount_from_mass(masse,q(:,:,iq),
+     &                                totaltracer_old(iq))
+        enddo
+      endif
 cc
 c
@@ -107,5 +117,16 @@
 c
         endif ! of if (planet_type.eq."earth")
-      ELSE
+        
+        ! Ehouarn : try to fix tracer conservation after tracer advection
+        if (force_conserv_tracer) then
+          do iq=1,nqtot
+            call planetary_tracer_amount_from_mass(masse,q(:,:,iq),
+     &                                  totaltracer_new(iq))
+            ratio=totaltracer_old(iq)/totaltracer_new(iq)
+            q(:,:,iq)=q(:,:,iq)*ratio
+          enddo
+        endif !of if (force_conserv_tracer)
+        
+      ELSE ! i.e. iapptrac.NE.iapp_tracvl
         if (planet_type.eq."earth") then
 ! Earth-specific treatment for the first 2 tracers (water)
Index: /trunk/LMDZ.COMMON/libf/dyn3d/conf_gcm.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d/conf_gcm.F90	(revision 2125)
+++ /trunk/LMDZ.COMMON/libf/dyn3d/conf_gcm.F90	(revision 2126)
@@ -17,5 +17,5 @@
                          output_grads_dyn, periodav, planet_type, &
                          raz_date, resetvarc, starttime, timestart, &
-                         ecritstart
+                         ecritstart,force_conserv_tracer
   USE infotrac, ONLY : type_trac
   use assert_m, only: assert
@@ -233,4 +233,8 @@
   iapp_tracvl = iperiod
   CALL getin('iapp_tracvl',iapp_tracvl)
+
+! Enforce tracer conservation in integrd & caladvtrac ?
+  force_conserv_tracer = .false. 
+  CALL getin('force_conserv_tracer',force_conserv_tracer)
 
 !Config  Key  = iconser
@@ -949,4 +953,6 @@
  write(lunout,*)' iphysiq = ', iphysiq
  write(lunout,*)' iflag_trac = ', iflag_trac
+ write(lunout,*)' iapp_tracvl = ', iapp_tracvl
+ write(lunout,*)' force_conserv_tracer = ', force_conserv_tracer
  write(lunout,*)' clon = ', clon
  write(lunout,*)' clat = ', clat
Index: /trunk/LMDZ.COMMON/libf/dyn3d/integrd.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d/integrd.F	(revision 2125)
+++ /trunk/LMDZ.COMMON/libf/dyn3d/integrd.F	(revision 2126)
@@ -7,5 +7,5 @@
      &  )
 
-      use control_mod, only : planet_type
+      use control_mod, only : planet_type,force_conserv_tracer
       USE comvert_mod, ONLY: ap,bp
       USE comconst_mod, ONLY: pi
@@ -66,4 +66,5 @@
       REAL vscr( ip1jm ),uscr( ip1jmp1 ),hscr( ip1jmp1 ),pscr(ip1jmp1)
       REAL massescr( ip1jmp1,llm )
+      REAL :: massratio(ip1jmp1,llm)
 !      REAL finvmasse(ip1jmp1,llm)
       REAL p(ip1jmp1,llmp1)
@@ -242,4 +243,15 @@
 
       endif ! of if (planet_type.eq."earth")
+
+      if (force_conserv_tracer) then
+        ! Ehouarn: try to keep total amont of tracers fixed 
+        ! by acounting for mass change in each cell
+        massratio(1:ip1jmp1,1:llm)=massescr(1:ip1jmp1,1:llm)
+     &                             /masse(1:ip1jmp1,1:llm)
+        do iq=1,nq
+        q(1:ip1jmp1,1:llm,iq)=q(1:ip1jmp1,1:llm,iq)
+     &                        *massratio(1:ip1jmp1,1:llm)
+        enddo
+      endif ! of if (force_conserv_tracer)
 c
 c
Index: /trunk/LMDZ.COMMON/libf/dyn3d/planetary_operations.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d/planetary_operations.F90	(revision 2126)
+++ /trunk/LMDZ.COMMON/libf/dyn3d/planetary_operations.F90	(revision 2126)
@@ -0,0 +1,133 @@
+!
+! $Id: $
+!
+module planetary_operations
+! module which contains functions to obtain total values over the
+! entire globe (trivial in serial mode, but not in parallel)
+
+implicit none
+
+contains
+
+  subroutine planetary_atm_mass_from_ps(ps,value)
+  implicit none
+  include "dimensions.h"
+  include "paramet.h"
+  include "comgeom.h"
+  real,intent(in) :: ps(ip1jmp1)
+  real,intent(out) :: value
+  integer :: i,j,ij
+
+  ! compute total atmospheric mass over the whole planet
+  ! North Pole
+  value=ps(1)*airesurg(1)*iim
+  do j=2,jjm
+    do i=1,iim
+      ij=i+(j-1)*iip1
+      value=value+ps(ij)*airesurg(ij)
+    enddo
+  enddo
+  ! South pole
+  value=value+ps(ip1jmp1-iim)*airesurg(ip1jmp1-iim)*iim
+  
+  end subroutine planetary_atm_mass_from_ps
+
+  subroutine planetary_atm_mass_from_mass(mass,value)
+  implicit none
+  include "dimensions.h"
+  include "paramet.h"
+  real,intent(in) :: mass(ip1jmp1,llm) ! air mass (kg)
+  real,intent(out) :: value
+  integer :: i,j,ij,l
+  real :: column(ip1jmp1) ! columns amount of air (kg)
+
+  ! compute total atmospheric mass over the whole planet
+  ! 1. build column amout of tracer (kg)
+  column(:)=0
+  do ij=1,ip1jmp1
+    do l=1,llm
+      column(ij)=column(ij)+mass(ij,l)
+    enddo
+  enddo
+  
+  ! 2. Compute total tracer mass over the planet
+  !North pole
+  value=column(1)*iim
+  do j=2,jjm
+    do i=1,iim
+      ij=i+(j-1)*iip1
+      value=value+column(ij)
+    enddo
+  enddo
+  ! South pole
+  value=value+column(ip1jmp1-iim)*iim
+  
+  end subroutine planetary_atm_mass_from_mass
+  
+  subroutine planetary_tracer_amount_from_p(p,q,amount)
+  use comconst_mod, ONLY: g
+  implicit none
+  include "dimensions.h"
+  include "paramet.h"
+  include "comgeom.h"
+  real,intent(in) :: p(ip1jmp1,llmp1) ! pressure at vertical mesh interfaces
+  real,intent(in) :: q(ip1jmp1,llm) ! 3D tracer (kg/kg_air)
+  real,intent(out) :: amount
+  integer :: i,j,ij,l
+  real :: column(ip1jmp1) ! columns amount of tracer (kg.m-2)
+  
+  ! 1. build column amout of tracer (kg.m-2)
+  column(:)=0
+  do ij=1,ip1jmp1
+    do l=1,llm
+      column(ij)=column(ij)+q(ij,l)*(p(ij,l)-p(ij,l+1))/g
+    enddo
+  enddo
+  
+  ! 2. Compute total tracer mass over the planet
+  !North pole
+  amount=column(1)*aire(1)*iim
+  do j=2,jjm
+    do i=1,iim
+      ij=i+(j-1)*iip1
+      amount=amount+column(ij)*aire(ij)
+    enddo
+  enddo
+  ! South pole
+  amount=amount+column(ip1jmp1-iim)*aire(ip1jmp1-iim)*iim
+
+  end subroutine planetary_tracer_amount_from_p
+
+  subroutine planetary_tracer_amount_from_mass(mass,q,amount)
+  implicit none
+  include "dimensions.h"
+  include "paramet.h"
+  real,intent(in) :: mass(ip1jmp1,llm) ! air mass (kg)
+  real,intent(in) :: q(ip1jmp1,llm) ! 3D tracer (kg/kg_air)
+  real,intent(out) :: amount
+  integer :: i,j,ij,l
+  real :: column(ip1jmp1) ! columns amount of tracer (kg)
+  
+  ! 1. build column amout of tracer (kg)
+  column(:)=0
+  do ij=1,ip1jmp1
+    do l=1,llm
+      column(ij)=column(ij)+q(ij,l)*mass(ij,l)
+    enddo
+  enddo
+  
+  ! 2. Compute total tracer mass over the planet
+  !North pole
+  amount=column(1)*iim
+  do j=2,jjm
+    do i=1,iim
+      ij=i+(j-1)*iip1
+      amount=amount+column(ij)
+    enddo
+  enddo
+  ! South pole
+  amount=amount+column(ip1jmp1-iim)*iim
+
+  end subroutine planetary_tracer_amount_from_mass
+
+end module planetary_operations
Index: /trunk/LMDZ.COMMON/libf/dyn3d_common/control_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3d_common/control_mod.F90	(revision 2125)
+++ /trunk/LMDZ.COMMON/libf/dyn3d_common/control_mod.F90	(revision 2126)
@@ -29,4 +29,5 @@
   integer,save :: ip_ebil_dyn
   logical,save :: offline
+  logical,save :: force_conserv_tracer ! enforce conservation of tracer mass
   character(len=4),save :: config_inca
   character(len=10),save :: planet_type ! planet type ('earth','mars',...)
Index: /trunk/LMDZ.COMMON/libf/dyn3dpar/advtrac_p.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3dpar/advtrac_p.F90	(revision 2125)
+++ /trunk/LMDZ.COMMON/libf/dyn3dpar/advtrac_p.F90	(revision 2126)
@@ -18,6 +18,8 @@
   USE times
   USE infotrac, ONLY: nqtot, iadv
-  USE control_mod, ONLY: iapp_tracvl, day_step, planet_type
+  USE control_mod, ONLY: iapp_tracvl, day_step, planet_type, &
+                         force_conserv_tracer
   USE comconst_mod, ONLY: dtvr
+  USE planetary_operations_p, ONLY: planetary_tracer_amount_from_mass_p
   IMPLICIT NONE
   !
@@ -77,4 +79,14 @@
   REAL,SAVE :: teta_tmp(ip1jmp1,llm)
   REAL,SAVE :: pk_tmp(ip1jmp1,llm)
+  REAL :: totaltracer_old(nqtot),totaltracer_new(nqtot)
+  REAL :: ratio
+
+! Ehouarn : try to fix tracer conservation issues:
+  if (force_conserv_tracer) then
+    do iq=1,nqtot
+       call planetary_tracer_amount_from_mass_p(masse,q(:,:,iq), &
+                                          totaltracer_old(iq))
+    enddo
+  endif
 
   ijb_u=ij_begin
@@ -467,4 +479,16 @@
      endif ! of if (planet_type=="earth")
 
+     ! Ehouarn : try to fix tracer conservation after tracer advection
+     if (force_conserv_tracer) then
+       ijb=ij_begin
+       ije=ij_end
+       do iq=1,nqtot
+         call planetary_tracer_amount_from_mass_p(masse,q(:,:,iq), &
+                                     totaltracer_new(iq))
+         ratio=totaltracer_old(iq)/totaltracer_new(iq)
+         q(ijb:ije,1:llm,iq)=q(ijb:ije,1:llm,iq)*ratio
+       enddo
+     endif !of if (force_conserv_tracer)
+
         !------------------------------------------------------------------
         !   on reinitialise a zero les flux de masse cumules
Index: /trunk/LMDZ.COMMON/libf/dyn3dpar/conf_gcm.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3dpar/conf_gcm.F90	(revision 2125)
+++ /trunk/LMDZ.COMMON/libf/dyn3dpar/conf_gcm.F90	(revision 2126)
@@ -255,4 +255,8 @@
   iapp_tracvl = iperiod
   CALL getin('iapp_tracvl',iapp_tracvl)
+
+! Enforce tracer conservation in integrd & caladvtrac ?
+  force_conserv_tracer = .false. 
+  CALL getin('force_conserv_tracer',force_conserv_tracer)
 
 !Config  Key  = iconser
@@ -987,4 +991,6 @@
       write(lunout,*)' iphysiq = ', iphysiq
       write(lunout,*)' iflag_trac = ', iflag_trac
+      write(lunout,*)' iapp_tracvl = ', iapp_tracvl
+      write(lunout,*)' force_conserv_tracer = ', force_conserv_tracer
       write(lunout,*)' clon = ', clon
       write(lunout,*)' clat = ', clat
Index: /trunk/LMDZ.COMMON/libf/dyn3dpar/integrd_p.F
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3dpar/integrd_p.F	(revision 2125)
+++ /trunk/LMDZ.COMMON/libf/dyn3dpar/integrd_p.F	(revision 2126)
@@ -7,5 +7,5 @@
       USE parallel_lmdz, ONLY: ij_begin, ij_end, pole_nord, pole_sud,
      &                         omp_chunk
-      USE control_mod, only : planet_type
+      USE control_mod, only : planet_type,force_conserv_tracer
       USE comvert_mod, ONLY: ap,bp
       USE comconst_mod, ONLY: pi
@@ -65,4 +65,5 @@
       REAL vscr( ip1jm ),uscr( ip1jmp1 ),hscr( ip1jmp1 ),pscr(ip1jmp1)
       REAL massescr( ip1jmp1,llm )
+      REAL :: massratio(ip1jmp1,llm)
 !      REAL finvmasse(ip1jmp1,llm)
       REAL,SAVE :: p(ip1jmp1,llmp1)
@@ -361,4 +362,15 @@
       endif ! of if (planet_type.eq."earth")
 
+
+      if (force_conserv_tracer) then
+        ! Ehouarn: try to keep total amont of tracers fixed 
+        ! by acounting for mass change in each cell
+        massratio(ijb:ije,1:llm)=massescr(ijb:ije,1:llm)
+     &                             /masse(ijb:ije,1:llm)
+        do iq=1,nq
+        q(ijb:ije,1:llm,iq)=q(ijb:ije,1:llm,iq)
+     &                        *massratio(ijb:ije,1:llm)
+        enddo
+      endif ! of if (force_conserv_tracer)
 c
 c
Index: /trunk/LMDZ.COMMON/libf/dyn3dpar/planetary_operations_p.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/dyn3dpar/planetary_operations_p.F90	(revision 2126)
+++ /trunk/LMDZ.COMMON/libf/dyn3dpar/planetary_operations_p.F90	(revision 2126)
@@ -0,0 +1,108 @@
+!
+! $Id: $
+!
+module planetary_operations_p
+! module which contains functions to obtain total values over the
+! entire globe (in parallel)
+
+implicit none
+
+contains
+
+  subroutine planetary_atm_mass_from_ps_p(ps,value)
+  USE parallel_lmdz, ONLY: Gather_Field, mpi_rank
+  USE mod_const_mpi, ONLY: COMM_LMDZ
+  implicit none
+  include "dimensions.h"
+  include "paramet.h"
+  include "comgeom.h"
+  include 'mpif.h'
+  real,intent(in) :: ps(ip1jmp1)
+  real,intent(out) :: value
+  integer :: i,j,ij
+  integer :: ierr
+
+  ! compute total atmospheric mass over the whole planet
+  ! 1. Gather ps on master to do the computations
+  call Gather_Field(ps,ip1jmp1,1,0)
+
+  ! 2. compute on master
+  if (mpi_rank==0) then
+    ! North Pole
+    value=ps(1)*airesurg(1)*iim
+    do j=2,jjm
+      do i=1,iim
+        ij=i+(j-1)*iip1
+        value=value+ps(ij)*airesurg(ij)
+      enddo
+    enddo
+    ! South pole
+    value=value+ps(ip1jmp1-iim)*airesurg(ip1jmp1-iim)*iim
+  endif ! of if (mpi_rank==0)
+
+  ! broadcast the result to all procs
+!$OMP CRITICAL (MPI_BCAST)
+#ifdef CPP_MPI
+  call MPI_BCAST(value,1,MPI_REAL8,0,COMM_LMDZ,ierr)
+#endif
+!$OMP END CRITICAL (MPI_BCAST)
+  
+  end subroutine planetary_atm_mass_from_ps_p
+
+  subroutine planetary_tracer_amount_from_mass_p(mass,q,amount)
+  USE parallel_lmdz, ONLY: ij_begin,ij_end, &
+                           pole_nord,pole_sud, &
+                           Gather_Field, mpi_rank
+  USE mod_const_mpi, ONLY: COMM_LMDZ
+  implicit none
+  include "dimensions.h"
+  include "paramet.h"
+  include 'mpif.h'
+  real,intent(in) :: mass(ip1jmp1,llm) ! air mass (kg)
+  real,intent(in) :: q(ip1jmp1,llm) ! 3D tracer (kg/kg_air)
+  real,intent(out) :: amount
+  integer :: i,j,ij,l
+  integer :: ijb,ije
+  integer :: ierr
+  real :: column(ip1jmp1) ! columns amount of tracer (kg)
+  
+  ! 1. build column amout of tracer (kg) on each proc
+  ijb=ij_begin-iip1
+  ije=ij_end+2*iip1
+  if (pole_nord) ijb=ij_begin
+  if (pole_sud)  ije=ij_end
+  
+  column(ijb:ije)=0
+  do ij=ijb,ije
+    do l=1,llm
+      column(ij)=column(ij)+q(ij,l)*mass(ij,l)
+    enddo
+  enddo
+  
+  ! 2 Gather "column" to do the upcoming computations on master
+  call Gather_Field(column,ip1jmp1,1,0)
+  
+  ! 3. Compute total tracer mass over the planet on master
+  if (mpi_rank==0) then
+    !North pole
+    amount=column(1)*iim
+    do j=2,jjm
+      do i=1,iim
+        ij=i+(j-1)*iip1
+        amount=amount+column(ij)
+      enddo
+    enddo
+    ! South pole
+    amount=amount+column(ip1jmp1-iim)*iim
+  endif ! of if (mpi_rank==0)
+  
+  ! broadcast the result to all procs
+!$OMP CRITICAL (MPI_BCAST)
+#ifdef CPP_MPI
+  call MPI_BCAST(amount,1,MPI_REAL8,0,COMM_LMDZ,ierr)
+#endif
+!$OMP END CRITICAL (MPI_BCAST)
+
+  end subroutine planetary_tracer_amount_from_mass_p
+
+end module planetary_operations_p
