source: trunk/LMDZ.PLUTO/libf/phypluto/spreadglacier_simple.F90 @ 3607

      subroutine spreadglacier_simple(ngrid,nq,pqsurf,timstep)

      use geometry_mod, only: latitude, longitude, cell_area
      use surfdat_h, only: phisfi
      use tracer_h, only: igcm_n2

      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"

!-----------------------------------------------------------------------
!     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
      DATA firstcall/.true./
      INTEGER,SAVE,ALLOCATABLE :: indglac(:)
      INTEGER, 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
         ALLOCATE(indglac(ngrid))
         nglac=0
         indglac(:)=0
         DO ig=1,ngrid
            IF (phisfi(ig).lt.phiref) THEN
               nglac=nglac+1
               indglac(nglac)=ig
            ENDIF
         ENDDO
         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=',latitude(indglac(ig2))*180/3.14,longitude(indglac(ig2))*180/3.14
             !.and.(phisfi(indglac(ig2)).le.phisfi(indglac(ig)))) THEN
             locdist=dist_pluto2(latitude(indglac(ig2)),longitude(indglac(ig2)),latitude(indglac(ig)),longitude(indglac(ig)))
             !write(*,*) 'dist=',locdist
             IF (locdist.lt.dist) THEN
              !write(*,*) 'calcul=',latitude(indglac(ig2))*180/3.14,longitude(indglac(ig2))*180/3.14,latitude(indglac(ig))*180/3.14,longitude(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*cell_area(indglac(ig))/cell_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