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

      use comcstfi_mod, only: pi, g
      use geometry_mod, only: latitude, longitude, cell_area
      USE tracer_h, only: igcm_co_ice, igcm_ch4_ice
      USE comgeomfi_h, only: totarea,totarea_planet
      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"

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

      INTEGER ngrid, nlayer, nq
      INTEGER igcm1, igcm2
      REAL pq(ngrid,nlayer,nq)
      REAL pdq(ngrid,nlayer,nq)
      REAL zdq(ngrid,nlayer,nq)
      REAL pqs(ngrid,nq)
      REAL pdqs(ngrid,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
      REAL pqc(ngrid,nlayer,nq)
      REAL pqcs(ngrid,nq)

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

      dWtot=0.0
      dWtots=0.0
      Wtot=0.0
      Wtots=0.0
      nconsMAX=0.0
      pqc=pq+pdq*ptimestep
      pqcs=pqs+pdqs*ptimestep
      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
             Wtot = Wtot + masse*pqc(ig,l,iq)*cell_area(ig)
             dWtot = dWtot + masse*zdq(ig,l,iq)*ptimestep*cell_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*cell_area(ig)
              Wtot = Wtot + masse*pqc(ig,l,iq)*cell_area(ig)
              vdifcncons(ig)=vdifcncons(ig) + masse*zdq(ig,l,iq)
             ENDIF
               
         enddo

         iq     = igcm1
         dWtots = dWtots + zdqs(ig,iq)*ptimestep*cell_area(ig)
         Wtots = Wtots + pqcs(ig,iq)*cell_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*cell_area(ig)
              Wtots = Wtots + pqcs(ig,iq)*cell_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(abs(vdifcncons(ig)).gt.abs(nconsMAX))then
              nconsMAX=vdifcncons(ig)
              myig=ig
        endif

      enddo

      dWtot  = dWtot/totarea_planet
      dWtots = dWtots/totarea_planet
      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=',latitude(myig)*180./pi,longitude(myig)*180./pi
      ENDIF
      print*,'--> Total Mass ',car1,' =',Wtot+Wtots,' kg m-2'
      end subroutine testconserv