subroutine testchange(ngrid,nlayer,nq,igcm1,igcm2,ptimestep, &
          pplev,zdq,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 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 nconsMAX
      REAL vdifcncons(ngrid)
!-----------------------------------------------------------------------


      dWtot=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
               dWtot = dWtot + masse*zdq(ig,l,iq)*ptimestep*area(ig)
               vdifcncons(ig)=vdifcncons(ig) + masse*zdq(ig,l,iq)

               iq    = igcm2
               dWtot = dWtot + masse*zdq(ig,l,iq)*ptimestep*area(ig)
               vdifcncons(ig)=vdifcncons(ig) + masse*zdq(ig,l,iq)
               
         enddo

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

      enddo

      dWtot  = dWtot/totarea
      print*,'-------------------------------------------'
      print*,'In ',car2,' atmospheric ',car1,' change=',dWtot,' kg m-2'
      print*,'--> MAX non-cons factor =',nconsMAX,' kg m-2 s-1'
 
      end subroutine testchange