subroutine testconservfast(ngrid,nlayer,nq,ptimestep, &
          pplev,zdq,zdqs,car1,car2)

      use comgeomfi_h          
      use comcstfi_mod, only: pi, g
      use geometry_mod, only: latitude, longitude, cell_area
      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
      REAL zdq(ngrid)
      REAL zdqs(ngrid)
      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)
!-----------------------------------------------------------------------

      dWtot=0.0
      dWtots=0.0
      nconsMAX=0.0
      do ig = 1, ngrid
         vdifcncons(ig)=0.0

        ! sum of atmospheric mass : kg
         dWtot = dWtot + zdq(ig)*ptimestep*cell_area(ig)*pplev(ig,1)/g
        ! for each column, total mass lost per sec : kg(tracer) / m2 /s
         vdifcncons(ig)=vdifcncons(ig) + pplev(ig,1)/g*zdq(ig)
 
         dWtots = dWtots + zdqs(ig)*ptimestep*cell_area(ig)
         vdifcncons(ig)=vdifcncons(ig)+zdqs(ig)

             
             ! 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

      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=',latitude(myig)*180./pi,longitude(myig)*180./pi
      ENDIF
      end subroutine testconservfast