subroutine spreadglacier_simple(ngrid,nq,pqsurf,timstep)
          
      use radinc_h, only : naerkind
      use comgeomfi_h, only : long,lati,area
      implicit none

!==================================================================
!     Purpose
!     -------
!     Spreading the glacier in a circular crater using a speed of flow
!  
!     Inputs
!     ------ 
!     ngrid                 Number of vertical columns
!     pqsurf(ngrid,nq)      N2 ice tracer on surface kg/m2 
!     
!     Outputs
!     -------
!     pqsurf(ngrid,nq)      new value for N2 ice tracer on surface kg/m2 
!     
!     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, nq, ig, ig2
      REAL :: pqsurf(ngrid,nq)
      REAL timstep 
!-----------------------------------------------------------------------
!     Local variables

      REAL :: zdqsurf(ngrid)
      REAL phiref
      REAL dist
      REAL locdist
      REAL tau
      REAL speed
      REAL dist_pluto2  ! function
      LOGICAL firstcall
      SAVE firstcall
      DATA firstcall/.true./
      INTEGER indglac(ngridmx)  !
      SAVE indglac
      INTEGER nglac
      SAVE nglac

!---------------- INPUT ------------------------------------------------
! Defined for  Tombaugh crater : 

      dist=300.  ! reference radius km
      phiref=-1000. ! geopotential reference
!      tau=1.e10
      speed=6.4e-7   ! glacier speed km/s
!      tau=1. ! direct redistribution glacier

!--------------- FIRSTCALL : define crater index -----------------------
      IF (firstcall) then
         indglac(:)=0
         DO ig=1,ngrid
            IF (phisfi(ig).lt.phiref) THEN
               nglac=nglac+1
               indglac(nglac)=ig
            ENDIF
         ENDDO 
         !write(*,*) 'inglac=',indglac
         write(*,*) 'nglac=',nglac
 
         firstcall=.false. 

      ENDIF

!--------------- redistribution -------------------------------------
      DO ig=1,nglac
        ! redistribution if large amount of ice only : detection of reservoirs
        IF (pqsurf(indglac(ig),igcm_n2).gt.50.) THEN
         ! looking where to redistribute : if less ice, lower altitude
         ! and close to the reservoir
         DO ig2=1,nglac
           IF ((pqsurf(indglac(ig2),igcm_n2).lt.pqsurf(indglac(ig),igcm_n2)*0.8)) THEN    
             !write(*,*) 'loop=',lati(indglac(ig2))*180/3.14,long(indglac(ig2))*180/3.14
             !.and.(phisfi(indglac(ig2)).le.phisfi(indglac(ig)))) THEN
             locdist=dist_pluto2(lati(indglac(ig2)),long(indglac(ig2)),lati(indglac(ig)),long(indglac(ig)))
             !write(*,*) 'dist=',locdist
             IF (locdist.lt.dist) THEN
              !write(*,*) 'calcul=',lati(indglac(ig2))*180/3.14,long(indglac(ig2))*180/3.14,lati(indglac(ig))*180/3.14,long(indglac(ig))*180/3.14
              tau=locdist/speed  ! characteristic time for redistribution
              zdqsurf(indglac(ig2))=pqsurf(indglac(ig),igcm_n2)*(1.-exp(-timstep/tau))/timstep*area(indglac(ig))/area(indglac(ig2))
              zdqsurf(indglac(ig))=pqsurf(indglac(ig),igcm_n2)*(exp(-timstep/tau)-1.)/timstep
              pqsurf(indglac(ig2),igcm_n2)=pqsurf(indglac(ig2),igcm_n2)+zdqsurf(indglac(ig2))*timstep
              pqsurf(indglac(ig),igcm_n2)=pqsurf(indglac(ig),igcm_n2)+zdqsurf(indglac(ig))*timstep
             ENDIF
           ENDIF
         ENDDO
        ENDIF 
      ENDDO

      end subroutine spreadglacier_simple

      real function dist_pluto2(lat1,lon1,lat2,lon2)
      implicit none
      !#include "comcstfi.h"
      real, intent(in) ::  lat1,lon1,lat2,lon2
      real dlon,dlat
      real a,c,d
      real rad
      rad=1187. ! planetary radius (km) 1187 km sur Pluton

      dlon = lon2 - lon1
      dlat = lat2 - lat1
      a = (sin(dlat/2))**2 + cos(lat1) * cos(lat2) * (sin(dlon/2))**2
      c = 2 * atan2( sqrt(a), sqrt(1-a) )
      dist_pluto2 = rad * c
      return 
      end function dist_pluto2