subroutine testconserv(ngrid,nlayer,nq,igcm1,igcm2,ptimestep, &
          pplev,zdq,zdqs,car1,car2)

      use comgeomfi_h          
      implicit none

!==================================================================
!     Purpose
!     -------
!     Test conservation of tracers
!  
!     Inputs
!     ------ 
!     ngrid                 Number of vertical columns
!     nlayer                Number of layers
!     pplev(ngrid,nlayer+1) Pressure levels 
!     zdq
!     
!     Outputs
!     -------
!     
!     Both
!     ----
!
!     Authors
!     ------- 
!     Tanguy Bertrand 2015
!     
!==================================================================

#include "dimensions.h"
#include "dimphys.h"
#include "comcstfi.h"
#include "surfdat.h"
#include "comvert.h"
#include "callkeys.h"
#include "tracer.h"

!-----------------------------------------------------------------------
!     Arguments

      INTEGER ngrid, nlayer, nq
      INTEGER igcm1, igcm2
      REAL zdq(ngrid,nlayer,nq)
      REAL zdqs(ngrid,nq)
      REAL ptimestep 
      REAL pplev(ngrid,nlayer+1)
      character(len=3) :: car1      
      character(len=7) :: car2      

      ! LOCAL VARIABLES
      INTEGER l,ig,iq
      REAL masse

      ! OUTPUT
      REAL dWtot   
      REAL dWtots
      REAL nconsMAX
      INTEGER myig
      REAL vdifcncons(ngrid)
!-----------------------------------------------------------------------

      write(*,*) 'TB igcm=',igcm1,igcm2
      write(*,*) 'TB zdq1=',zdq(1,1,igcm1)
      dWtot=0.0
      dWtots=0.0
      nconsMAX=0.0
      do ig = 1, ngrid
         vdifcncons(ig)=0.0
         do l = 1, nlayer
             masse = (pplev(ig,l) - pplev(ig,l+1))/g

             iq    = igcm1
             ! sum of atmospheric mass : kg
             dWtot = dWtot + masse*zdq(ig,l,iq)*ptimestep*area(ig)
             ! for each column, total mass lost per sec : kg(tracer) / m2 /s
             vdifcncons(ig)=vdifcncons(ig) + masse*zdq(ig,l,iq)
 
             ! if clouds :
             IF ((igcm2.eq.igcm_co_ice).or.(igcm2.eq.igcm_ch4_ice)) THEN
              iq    = igcm2
              dWtot = dWtot + masse*zdq(ig,l,iq)*ptimestep*area(ig)
              vdifcncons(ig)=vdifcncons(ig) + masse*zdq(ig,l,iq)
             ENDIF
               
         enddo

             iq     = igcm1
             dWtots = dWtots + zdqs(ig,iq)*ptimestep*area(ig)
             vdifcncons(ig)=vdifcncons(ig)+zdqs(ig,iq)

             IF ((igcm2.eq.igcm_co_ice).or.(igcm2.eq.igcm_ch4_ice)) THEN
              iq     = igcm2
              dWtots = dWtots + zdqs(ig,iq)*ptimestep*area(ig)
              vdifcncons(ig)=vdifcncons(ig)+zdqs(ig,iq)
             ENDIF
             
             ! vdifcncons is the total amount of material that appear or
             ! disapear per second in the routine
             ! it is the non conservative factor 

             if(vdifcncons(ig).gt.nconsMAX)then
                nconsMAX=vdifcncons(ig)
                myig=ig
             endif

      enddo
      write(*,*) 'TB Dwtot=',dWtots,dWtot

      dWtot  = dWtot/totarea
      dWtots = dWtots/totarea
      print*,'-------------------------------------------'
      print*,'In ',car2,' atmospheric ',car1,' change=',dWtot,' kg m-2'
      print*,'In ',car2,' surface ',car1,' change  =',dWtots,' kg m-2'
      print*,'--> non-cons factor     =',dWtot+dWtots,' kg m-2'
      print*,'--> MAX non-cons factor =',nconsMAX,' kg m-2 s-1'
      IF (nconsMAX.gt.0.) then
        print*,'--> obtained at lat/lon=',lati(myig)*180./pi,long(myig)*180./pi
      ENDIF
      end subroutine testconserv