subroutine massn2(ngrid,nlayer,pplev,qsurf,dqsurf,pdpsurf,  & 
                       ptimestep)
  
      use comgeomfi_h    
      implicit none

!==================================================================
!     Purpose
!     -------
!     Test conservation of n2 mass 
!  
!     Inputs
!     ------ 
!     ngrid                 Number of vertical columns
!     nlayer                Number of layers
!     pplev(ngrid,nlayer+1) Pressure levels 
!     surfmass    
!
!     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
      REAL pplev(ngrid,nlayer+1)
      REAL qsurf(ngrid)
      REAL dqsurf(ngrid)
      REAL pdpsurf(ngrid)
      REAL ptimestep

      ! LOCAL VARIABLES
      INTEGER l,ig,iq
      REAL atmmass
      REAL pstot
      REAL zps(ngrid)
      REAL zqsurf(ngrid)

      ! OUTPUT
      REAL atmmasstot
      REAL surfmasstot
      REAL totmass
!-----------------------------------------------------------------------


      atmmass=0.0
      atmmasstot=0.0
      surfmasstot=0.0
      pstot=0.0
      zps(:)=0.
      zqsurf(:)=0.

      do ig = 1, ngrid
         zps(ig)=pplev(ig,1)+pdpsurf(ig)*ptimestep
         zqsurf(ig)=qsurf(ig)+dqsurf(ig)*ptimestep
 
         atmmass = zps(ig)/g
         pstot = pstot+zps(ig)*area(ig)

         atmmasstot= atmmasstot + atmmass*area(ig)
         surfmasstot= surfmasstot + zqsurf(ig)*area(ig)

       enddo
               
      atmmasstot=  atmmasstot/totarea
      pstot=  pstot/totarea
      surfmasstot = surfmasstot/totarea
      totmass=surfmasstot+atmmasstot

!      print*,'-------------------------------------------'
      print*,'Total mass surface : ',surfmasstot,' kg m-2'
      print*,'Total mass atmosphere :',atmmasstot,' kg m-2'
      print*,'Total mean pressure:',pstot,' Pa'
      print*,'Total mass : ',surfmasstot+atmmasstot,' kg m-2'
 

      end subroutine massn2