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

      use comcstfi_mod, only: pi, g
      use geometry_mod, only: latitude, longitude, cell_area
      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

      REAL pq(ngrid)
      REAL pdq(ngrid)
      REAL zdq(ngrid)
      REAL pqs(ngrid)
      REAL pdqs(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
      REAL pqc(ngrid)
      REAL pqcs(ngrid)

      ! 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

        ! sum of atmospheric mass : kg
        Wtot = Wtot + pqc(ig)*cell_area(ig)*pplev(ig,1)/g
        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)
        Wtots = Wtots + pqcs(ig)*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(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)/totarea_planet,' kg m-2'
      end subroutine testconservfast