subroutine testconservmass(ngrid,nlayer, &
          ps,surfmass)

      use comcstfi_mod, only: g
      use geometry_mod, only: cell_area
      USE comgeomfi_h, only: totarea
      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"

!-----------------------------------------------------------------------
!     Arguments
      INTEGER ngrid, nlayer
      REAL ps(ngrid)
      REAL surfmass(ngrid)

      ! LOCAL VARIABLES
      INTEGER l,ig,iq
      REAL atmmass
      REAL pstot

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

      atmmass=0.0
      atmmasstot=0.0
      surfmasstot=0.0
      pstot=0.0

      do ig = 2, ngrid-1

         ! mass atm kg/m2
         atmmass = ps(ig)/g
          
         pstot = pstot+ps(ig)*cell_area(ig)

         ! mass atm kg
         atmmasstot= atmmasstot + atmmass*cell_area(ig)

         surfmasstot= surfmasstot + surfmass(ig)*cell_area(ig)

      enddo
      atmmasstot= atmmasstot + ps(1)/g*cell_area(1)*iim
      atmmasstot= atmmasstot + ps(ngrid)/g*cell_area(ngrid)*iim

      surfmasstot= surfmasstot + surfmass(1)*cell_area(1)*iim
      surfmasstot= surfmasstot + surfmass(ngrid)*cell_area(ngrid)*iim

      pstot= pstot + ps(1)*cell_area(1)*iim
      pstot= pstot + ps(ngrid)*cell_area(ngrid)*iim
               
      atmmasstot=  atmmasstot/(totarea+(iim-1)*(cell_area(1)+cell_area(ngrid)))
      pstot=  pstot/(totarea+(iim-1)*(cell_area(1)+cell_area(ngrid)))
      surfmasstot = surfmasstot/(totarea+(iim-1)*(cell_area(1)+cell_area(ngrid)))

      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 testconservmass