subroutine hazesource(ngrid,nlayer,nq,ptimestep, &
          pplev,flusurf,mu0,zdq_source)

      use radinc_h, only : naerkind
      use comgeomfi_h
      implicit none

!==================================================================
!     Purpose
!     -------
!     Surface source of haze particle
!  
!     Inputs
!     ------ 
!     ngrid                 Number of vertical columns
!     nlayer                Number of layers
!     pplay(ngrid,nlayer)   Pressure layers
!     pplev(ngrid,nlayer+1) Pressure levels 
!     
!     Outputs
!     -------
!     zdq_source
!
!     Authors
!     ------- 
!     Tanguy Bertrand
!     
!==================================================================

#include "dimensions.h"
#include "dimphys.h"
#include "comcstfi.h"
#include "surfdat.h"
#include "comvert.h"
#include "callkeys.h"
#include "tracer.h"

!-----------------------------------------------------------------------
!     Arguments

      INTEGER,INTENT(IN) :: ngrid, nlayer, nq
      REAL,INTENT(IN) :: ptimestep 
      REAL,INTENT(IN) :: pplev(ngrid,nlayer+1)
      REAL,INTENT(IN) :: flusurf(ngrid,nq)     ! flux cond/sub kg.m-2.s-1
      REAL,INTENT(IN) :: mu0(ngrid)  ! cosinus of solar incident flux
      REAL,INTENT(OUT) :: zdq_source(ngrid,nlayer,nq)
!-----------------------------------------------------------------------
!     Local variables

      INTEGER l,ig,iq
      CHARACTER(len=10) :: tname
      real massfix(ngrid) !
      real tholflux(ngrid) !
!-----------------------------------------------------------------------

      !! From callphys.def : kfix, fracsource, latsource, mode_hs

      !! 1) Get mass of the layer (kg/m2) where haze is injected
      massfix(:)=0.0
      zdq_source(:,:,:)=0.0
      DO ig=1,ngrid
         DO l= 1,kfix
            massfix(ig)= massfix(ig) + (pplev(ig,l)-pplev(ig,l+1))/g
         ENDDO
      ENDDO


      SELECT CASE (mode_hs) ! mode haze source
         !!! Source haze in SP: 0.02 pourcent when n2 sublimes
         CASE (0) 
          DO ig=1,ngrid
            IF (phisfi(ig).le.-1500.) THEN  ! SP
               IF (lati(ig)*180./pi.ge.latsource.and.lati(ig)*180./pi.lt.65.) THEN

                   IF (flusurf(ig,igcm_n2).lt.0) THEN
                      IF (mu0(ig).gt.0.6) THEN
                       !! flux tholins kg/m2/s
                       tholflux(ig)=-flusurf(ig,igcm_n2)*fracsource
                       DO iq=1,nq
                        tname=noms(iq)
                        if (tname(1:4).eq."haze") then
                          DO l= 1,kfix
                            zdq_source(ig,l,iq)=tholflux(ig)/massfix(ig)
                          ENDDO
                        endif
                       ENDDO
                      ENDIF ! mu0
                   ENDIF ! flusurf

               ENDIF
            ENDIF
          ENDDO

         !!! Source haze at the BTD
         CASE (1) 
          DO ig=1,ngrid
            IF (phisfi(ig).ge.0.) THEN  
               !print*,'lat/lon=',lati(ig)*180./pi,long(ig)*180./pi
               IF (lati(ig)*180./pi.ge.0..and.lati(ig)*180./pi.lt.5.) THEN
                 IF (long(ig)*180./pi.ge.-50..and.long(ig)*180./pi.lt.-45.) THEN
                   !print*,'YES',flusurf(ig,igcm_ch4_ice)

                   IF (flusurf(ig,igcm_ch4_ice).lt.0) THEN
                       !! flux tholins kg/m2/s
                       tholflux(ig)=-flusurf(ig,igcm_ch4_ice)*fracsource
                       DO iq=1,nq
                        tname=noms(iq)
                        if (tname(1:4).eq."haze") then
                          DO l= 1,kfix
                            zdq_source(ig,l,iq)=tholflux(ig)/massfix(ig)
                          ENDDO
                        endif
                       ENDDO
                   ENDIF ! flusurf

                 ENDIF
               ENDIF
            ENDIF
          ENDDO

         !!! Basic source(geysers)
         CASE (2) 
          DO ig=1,ngrid
              IF ((abs(lati(ig)*180./pi-latsource).le.abs(lati(1)-lati(2))*180./pi)  &
             .and.(abs(long(ig)*180./pi-lonsource).le.abs(long(3)-long(4))*180./pi)) THEN

                       !! flux tholins kg/m2/s
                       tholflux(ig)=fracsource
                       DO iq=1,nq
                        tname=noms(iq)
                        if (tname(1:4).eq."haze") then
                          DO l= 1,kfix
                            zdq_source(ig,l,iq)=tholflux(ig)/massfix(ig)
                          ENDDO
                        endif
                       ENDDO
              ENDIF
          ENDDO


         CASE DEFAULT
            write(*,*) 'STOP in hazesource.F90: mod_hs not found'
            stop
      END SELECT


      end subroutine hazesource