Index: trunk/LMDZ.GENERIC/libf/phystd/callkeys_mod.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/callkeys_mod.F90	(revision 3294)
+++ trunk/LMDZ.GENERIC/libf/phystd/callkeys_mod.F90	(revision 3299)
@@ -6,4 +6,9 @@
       logical,save :: calladj,calltherm,co2cond,callsoil
 !$OMP THREADPRIVATE(calladj,calltherm,co2cond,callsoil)
+
+      ! do we have volcanoes (injecting tracers in atmosphere) ?
+      logical,save :: callvolcano
+!$OMP THREADPRIVATE(callvolcano)
+
       logical,save :: season,diurnal,tlocked,rings_shadow,lwrite
 !$OMP THREADPRIVATE(season,diurnal,tlocked,rings_shadow,lwrite)
Index: trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90	(revision 3294)
+++ trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90	(revision 3299)
@@ -503,4 +503,16 @@
 !       ": ok_slab_heat_transp = ",ok_slab_heat_transp
 
+! Volcanoes (as point sources of tracers)
+      if (is_master) write(*,*) trim(rname)//": Erupt volcano ?"
+      callvolcano=.false. ! default value
+      call getin_p("callvolcano",callvolcano)
+      if (is_master) write(*,*) trim(rname)//": callvolcano = ",callvolcano
+      
+      ! Sanity check: volcano requires tracers
+      if (callvolcano.and.(.not.tracer)) then
+        call abort_physic(rname,&
+             " Error: volcano option requires using tracers",1)
+      endif
+
 ! Photochemistry and chemistry in the thermosphere
 
Index: trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90	(revision 3294)
+++ trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90	(revision 3299)
@@ -66,5 +66,5 @@
                               ok_slab_ocean, photochem, rings_shadow, rmean, &
                               season, sedimentation,generic_condensation, &
-                              specOLR, callthermos, &
+                              specOLR, callthermos, callvolcano, &
                               startphy_file, testradtimes, tlocked, &
                               tracer, UseTurbDiff, water, watercond, &
@@ -77,5 +77,6 @@
       use phys_state_var_mod
       use callcorrk_mod, only: callcorrk
-      use conduction_mod
+      use conduction_mod, only: conduction
+      use volcano_mod, only: volcano
       use pindex_mod, only: pindex
       use vdifc_mod, only: vdifc
@@ -358,4 +359,5 @@
       REAL,allocatable,save :: zdqschim(:,:)  ! Calchim_asis routine
 !$OMP THREADPRIVATE(zdqchim,zdqschim)
+      real zdqvolc(ngrid,nlayer,nq) ! injection by volcanoes (kg/kg_of_air/s)
 
       REAL array_zero1(ngrid)
@@ -1481,4 +1483,14 @@
       if (tracer) then
 
+  ! -----------------------------------------------------
+  !   VI.0. Volcanoes injecting tracers in the atmosphere
+  ! -----------------------------------------------------
+        if (callvolcano) then
+          call volcano(ngrid,nlayer,pplev,pu,pv,pt,zzlev,zdqvolc,nq,massarea,&
+                       zday,ptimestep)
+          pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + &
+                                       zdqvolc(1:ngrid,1:nlayer,1:nq)
+        endif
+
   ! ---------------------
   !   VI.1. Water and ice
Index: trunk/LMDZ.GENERIC/libf/phystd/tracer_h.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/tracer_h.F90	(revision 3294)
+++ trunk/LMDZ.GENERIC/libf/phystd/tracer_h.F90	(revision 3299)
@@ -115,4 +115,40 @@
        !$OMP igcm_c2h3,igcm_c2h4,igcm_c2h6,igcm_ch2co,igcm_ch3co,igcm_hcaer)
 
+       contains
+       
+       ! some basic auxiliary functions for tracers
+       
+       logical function is_known_tracer(tracer_name)
+         ! checks if tracer_name is a know tracer
+         ! returns .true. if yes and .false. otherwise
+         implicit none
+         character(len=*),intent(in) :: tracer_name
+         integer :: iq
+         ! initialization
+         is_known_tracer=.false.
+         do iq=1,nqtot
+           if( trim(tracer_name) == trim(noms(iq)) ) then
+             is_known_tracer=.true.
+             exit
+           endif
+         enddo
+       end function is_known_tracer
+       
+       integer function tracer_index(tracer_name)
+         ! gets the index of the tracer tracer_name
+         ! returns index or 0 if unknown tracer
+         implicit none
+         character(len=*),intent(in) :: tracer_name
+         integer :: iq
+         ! initialization
+         tracer_index=0
+         do iq=1,nqtot
+           if( trim(tracer_name) == trim(noms(iq)) ) then
+             tracer_index=iq
+             exit
+           endif
+         enddo
+       end function tracer_index
+        
        end module tracer_h
 
Index: trunk/LMDZ.GENERIC/libf/phystd/volcano.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/volcano.F90	(revision 3299)
+++ trunk/LMDZ.GENERIC/libf/phystd/volcano.F90	(revision 3299)
@@ -0,0 +1,414 @@
+MODULE volcano_mod
+
+IMPLICIT NONE
+  ! saved,protected variables: (local to the core)
+  
+  integer,save,protected :: nvolc ! number of active volcanoes 
+!$OMP THREADPRIVATE(nvolc)
+
+  real,save,allocatable,protected :: latlonvals(:,:) ! dimension (nvolc,2)
+                                     ! lat, lon position of each volcano (deg)
+!$OMP THREADPRIVATE(latlonvals)
+  
+  integer,save,allocatable,protected :: heightvals(:)   ! dimension (nvolc)
+                                        ! altitude grid point index 
+                                        ! of eruption for each volcano
+!$OMP THREADPRIVATE(heightvals)
+
+  integer,save,allocatable,protected :: eruptfreqs(:,:) ! dimension (nvolc,2)
+                                        ! for each volcano: interval 
+                                        ! between each eruption in days, 
+                                        ! day offset from 0
+!$OMP THREADPRIVATE(eruptfreqs)
+  
+  real,save,allocatable,protected :: volfluxes(:,:) ! dimension (nvolc,nq)
+                                     ! outgassing flux of each tracer for 
+                                     ! each volcano, in kg/s
+!$OMP THREADPRIVATE(volfluxes)
+
+  integer,save,allocatable,protected :: ivolc(:) ! dimension (nvolc)
+                                        ! horizontal grid index of volcano
+!$OMP THREADPRIVATE(ivolc)
+
+CONTAINS
+
+  SUBROUTINE volcano(ngrid,nlay,pplev,wu,wv,pt,zzlev,ssource,nq, &
+                     massarea,zday,ptimestep)
+
+      use time_phylmdz_mod, only: daysec
+
+      IMPLICIT NONE
+
+!=======================================================================
+!   Volcanic Activity (updated to run in parallel)
+!
+!   Description :
+!     Author : Lionel Wilson
+!     Adapted for the LMD/GCM by Laura Kerber, J.-B. Madeleine, Saira Hamid
+!     Adapted to be more generalisable by Ashwin Braude (02.11.2023)
+!
+!   Reference :
+!       @ARTICLE{2007JVGR..163...83W,
+!       author = {{Wilson}, L. and {Head}, J.~W.},
+!       title = "{Explosive volcanic eruptions on Mars: Tephra and
+!       accretionary lapilli formation, dispersal and recognition in the
+!       geologic record}",
+!       journal = {Journal of Volcanology and Geothermal Research},
+!       year = 2007,
+!       month = jun,
+!       volume = 163}
+!=======================================================================
+
+!     Variable statement
+!     __________________________________________________________________
+
+
+! Parameters
+! ----------
+
+! PLEASE DEFINE THE LOCATION OF THE ERUPTION BELOW :
+! =============================================
+!     Source flux (kg/s)
+!       REAL, parameter :: mmsource = 1E8  !Mastin et al. 2009
+!       REAL, parameter :: wsource = 1.E6 !1wt% magma water content(Greeley 1987)
+!       REAL, parameter :: h2so4source = 1.E8
+! =============================================
+! Local variables
+! ---------------
+
+      INTEGER :: l
+      INTEGER :: iq                     ! Tracer identifier
+      INTEGER :: ig
+      INTEGER :: volci
+
+      LOGICAL,SAVE :: firstcall=.true.
+!$OMP THREADPRIVATE(firstcall)
+
+! Inputs
+! ------
+
+      INTEGER, intent(in) :: nq               ! Number of tracers
+      INTEGER, intent(in) :: ngrid            ! Number of grid points
+      INTEGER,intent(in) ::  nlay            ! Number of vertical levels
+      REAL,intent(in) :: pplev(ngrid,nlay+1) ! Pressure between the layers (Pa)
+      REAL,intent(in) :: zzlev(ngrid,nlay+1) ! height between the layers (km)
+      REAL,intent(in) :: wv(ngrid,nlay+1)    ! wind
+      REAL,intent(in) :: wu(ngrid,nlay+1)    ! wind
+      REAL,intent(in) :: pt(ngrid,nlay+1)    ! temp
+      REAL,intent(in) :: massarea(ngrid,nlay)
+      REAL,intent(in) :: zday ! "Universal time": in sols (and fraction thereof). 
+                              ! since the North. Spring equinox (Ls=0)
+      REAL,intent(in) :: ptimestep             ! Physics timestep (s).
+
+! Outputs
+! -------
+      REAL, intent(out) :: ssource(ngrid,nlay,nq)    ! Source tendency (kg.kg-1.s-1)
+
+!     Initialization
+!     __________________________________________________________________
+   
+      if (firstcall) then
+        ! initialize everything the very first time this routine is called
+        call read_volcano(nq,ngrid)
+        firstcall=.false.
+      endif
+      
+      ssource(1:ngrid,1:nlay,1:nq)=0 !all arrays=zero since it's a local variable within routine that need to be initialized
+      
+      do volci=1,nvolc
+        
+        !If time to erupt volcano
+        if(modulo(nint(zday)-eruptfreqs(volci,2),eruptfreqs(volci,1)).eq.0)then
+
+          !and then emit source flux at that point
+          WRITE(*,*) 'Volcano ',volci,' Time: ', nint(zday), eruptfreqs(volci,:)
+          WRITE(*,*) 'Erupt volcano at ivolc=',ivolc(volci),' height=',heightvals(volci)
+          l = heightvals(volci)
+          ig=ivolc(volci)
+          do iq=1,nq
+            ssource(ig,l,iq) = ssource(ig,l,iq)+&
+                                          volfluxes(volci,iq)/massarea(ig,l)
+          enddo
+
+        endif
+
+      enddo ! of do volci=1,nvolc
+
+  END SUBROUTINE volcano
+  
+  !=======================================================================
+  
+  SUBROUTINE read_volcano(nq,ngrid)
+  ! this routine reads in volcano.def and initializes module variables
+    USE mod_phys_lmdz_para, only: is_master, bcast
+    use mod_phys_lmdz_transfert_para, only: reduce_sum
+    use comcstfi_mod, only: pi
+    use tracer_h, only: is_known_tracer, tracer_index
+    use mod_grid_phy_lmdz,  only: nvertex ! number of grid cell "corners"
+    use geometry_mod, only: boundslon, boundslat ! grid cell boundaries (rad)
+    use geometry_mod, only: longitude_deg, latitude_deg ! grid cell coordinates (deg)
+    IMPLICIT NONE
+
+    integer, intent(in) :: nq ! number of tracers
+    INTEGER, intent(in) :: ngrid ! number of atmospheric columns
+    real, allocatable :: latlonvals_total(:,:) !lat, lon position of each volcano
+    integer, allocatable :: heightvals_total(:) ! alt index of eruption
+    real, allocatable :: eruptfreqs_total(:,:) ! date and freq of eruption
+    real, allocatable :: volfluxes_total(:,:) ! injected tracer rate 
+    integer :: ierr ! error return code
+    integer :: iq,l,volci,ivoltrac,nvoltrac,blank
+    integer :: eruptfreq,eruptoffset
+    integer :: nvolctrac ! number of outgased tracers
+    integer :: ig,i
+    integer :: nvolc_total ! total (planet-wide) number of volcanoes
+    integer :: nvolc_check
+    real    :: lat_volc,lon_volc,volflux
+    character(len=500) :: voltrac(nq)
+    character(len=500) :: tracline   ! to read volcano.def line
+    real :: boundslon_deg(nvertex) ! cell longitude boundaries (deg)
+    real :: boundslat_deg(nvertex) ! cell latitude boundaries (deg)
+    real :: minlat, minlon, maxlat, maxlon
+    logical,allocatable :: belongs_here(:)
+    integer,allocatable :: volc_index(:)
+    character(len=*),parameter :: rname="read_volcano"
+
+    ! 1. Only the master reads and processes the volcano.def file
+    if (is_master) then
+      ! 1.1 Open file and get total number of volcanoes
+        open(407, form = 'formatted', status = 'old', &
+                         file = 'volcano.def', iostat=ierr)
+        if (ierr /=0) then
+          print*,trim(rname)//': Problem in opening volcano.def'
+          stop
+        end if
+        
+        write(*,*)trim(rname)//' Reading file volcano.def'
+        DO
+          READ(407,'(A)',iostat=ierr) tracline
+          IF (ierr==0) THEN
+            IF (index(tracline,'#').ne.1) THEN ! Allows arbitary number of comments lines in the header
+              !Read in the number of point sources/volcanoes in volcano.def
+              READ(tracline,*) nvolc_total
+              EXIT
+            ENDIF
+          ELSE ! If pb, or if reached EOF without having found nqtot
+          !        call abort_physic('initracer','Unable to read numbers '// &
+          !  'of tracers in traceur.def',1)
+              print*,trim(rname)//": unable to read numbers of tracer in volcano.def"
+              stop
+          ENDIF
+        ENDDO
+    endif ! of if (is_master)
+    ! tell the world about nvolc_total
+    CALL bcast(nvolc_total)
+    
+    ! 1.2. allocate arrays (all cores) to store all the information
+    ALLOCATE(latlonvals_total(nvolc_total,2),stat=ierr)
+    if (ierr/=0) then
+      write(*,*) trim(rname)//" failed allocation for latlonvals_total()"
+      call abort_physic(rname,"failled allocation",1)
+    endif
+    ALLOCATE(heightvals_total(nvolc_total),stat=ierr)
+    if (ierr/=0) then
+      write(*,*) trim(rname)//" failed allocation for heightvals_total()"
+      call abort_physic(rname,"failled allocation",1)
+    endif
+    ALLOCATE(eruptfreqs_total(nvolc_total,2),stat=ierr)
+    if (ierr/=0) then
+      write(*,*) trim(rname)//" failed allocation for eruptfreqs_total()"
+      call abort_physic(rname,"failled allocation",1)
+    endif
+    ALLOCATE(volfluxes_total(nvolc_total,nq),stat=ierr)
+    if (ierr/=0) then
+      write(*,*) trim(rname)//" failed allocation for volfluxes_total()"
+      call abort_physic(rname,"failled allocation",1)
+    endif
+    ! initialize volfluxes_total
+    volfluxes_total(:,:)=0.0
+    
+    ! 1.3. continue (master only) to read and process volcano.def
+    if (is_master) then
+      volci = 0
+      !Iterate for each volcano
+      do while (volci.lt.nvolc_total)
+        DO
+          READ(407,'(A)',iostat=ierr) tracline
+          IF (ierr==0) THEN
+            IF (index(tracline,'#').ne.1) THEN 
+              READ(tracline,*) volci
+              EXIT
+            ENDIF
+          ELSE ! If pb, or if reached EOF without having found nqtot
+        !        call abort_physic('initracer','Unable to read numbers '// &
+        !  'of tracers in traceur.def',1)
+            print*,trim(rname)//": unable to read line in volcano.def"
+            stop
+          ENDIF
+        ENDDO
+        write(*,*)trim(rname)//': volcano number:',volci
+        READ(407,*,iostat=ierr) lat_volc,lon_volc,l
+        write(*,*)trim(rname)//' latitude:',lat_volc
+        write(*,*)trim(rname)//' longitude:',lon_volc
+        latlonvals_total(volci,1) = lat_volc
+        latlonvals_total(volci,2) = lon_volc
+        heightvals_total(volci) = l
+        !Read in number of active tracers being outgassed
+        READ(407,*,iostat=ierr) nvoltrac
+        write(*,*)trim(rname)//' number of outgased tracers:',nvoltrac
+        !Read in eruption frequency (in multiples of iphysiq) and offset of eruption
+        READ(407,'(A)',iostat=ierr) tracline
+        READ(tracline,*) eruptfreq, eruptoffset
+        write(*,*)trim(rname)//' eruption frequency (sols):',eruptfreq
+        write(*,*)trim(rname)//' eruption offset (sols):',eruptoffset
+        eruptfreqs_total(volci,1) = eruptfreq
+        eruptfreqs_total(volci,2) = eruptoffset
+        
+        do ivoltrac=1,nvoltrac
+
+          !Read in tracer name and check it's in traceur.def
+          read(407,'(A)') tracline
+          blank = index(tracline,' ') ! Find position of 1st blank in tracline
+          voltrac(ivoltrac) = tracline(1:blank)
+          ! check it is a valid tracer:
+          if (is_known_tracer(voltrac(ivoltrac))) then
+            !Read in flux of tracer out of volcano (kg/s)
+            READ(407,'(A)',iostat=ierr) tracline
+            READ(tracline,*) volflux
+            iq=tracer_index(voltrac(ivoltrac))
+            volfluxes_total(volci,iq) = volflux
+          else
+            write(*,*) trim(rname)//": problem with ",trim(voltrac(ivoltrac))
+            call abort_physic(rname, 'tracer '//&
+                              trim(voltrac(ivoltrac))//' not found', 1)
+          endif
+          
+        enddo ! of do ivoltrac=1,nvoltrac
+        
+      enddo ! of do while (volci.lt.nvolc_total)
+    
+      close(407)
+    endif ! of if (is_master)
+    
+    ! 2. Share data with all cores
+    call bcast(latlonvals_total)
+    call bcast(heightvals_total)
+    call bcast(eruptfreqs_total)
+    call bcast(volfluxes_total)
+  
+    ! 3. Each core extracts and stores data relevant to its domain only
+    ! 3.1. find out how many volcanoes for this core
+    nvolc=0
+    allocate(belongs_here(nvolc_total),stat=ierr)
+    if (ierr/=0) then
+      write(*,*) trim(rname)//" failed allocation for belong_here()"
+      call abort_physic(rname,"failled allocation",1)
+    endif
+    allocate(volc_index(nvolc_total),stat=ierr)
+    if (ierr/=0) then
+      write(*,*) trim(rname)//" failed allocation for volc_index()"
+      call abort_physic(rname,"failled allocation",1)
+    endif
+    
+    do volci=1,nvolc_total
+      belongs_here(volci)=.false.
+      volc_index(volci)=0
+      
+      lat_volc=latlonvals_total(volci,1)
+      lon_volc=latlonvals_total(volci,2)
+
+      ! find out if this volcano is in the ig cell
+      do ig=1,ngrid
+        boundslon_deg(:)=boundslon(ig,:)/pi*180.
+        boundslat_deg(:)=boundslat(ig,:)/pi*180.
+        
+        minlon=minval(boundslon_deg)
+        maxlon=maxval(boundslon_deg)
+        minlat=minval(boundslat_deg)
+        maxlat=maxval(boundslat_deg)
+        
+        if ((minlat<=lat_volc).and.(lat_volc<=maxlat).and. &
+            (((minlon<=lon_volc).and.(lon_volc<=maxlon)).or. &
+             ((minlon+360.<=lon_volc).and.(lon_volc<=maxlon+360.)).or. &
+             ((minlon-360.<=lon_volc).and.(lon_volc<=maxlon-360.)))) then
+          nvolc=nvolc+1
+          belongs_here(volci)=.true.
+          volc_index(volci)=ig
+          ! the job is done move on to next volcano
+          exit
+        endif
+      enddo ! of do ig=1,ngrid
+    enddo ! of volci=1,nvolc_total
+    
+    ! 3.2. allocate arrays, now that nvolc is known
+    if (nvolc>0) then
+      allocate(latlonvals(nvolc,2),stat=ierr)
+      if (ierr/=0) then
+        write(*,*) trim(rname)//" failed allocation for latlonvals()"
+        call abort_physic(rname,"failled allocation",1)
+      endif
+      allocate(heightvals(nvolc),stat=ierr)
+      if (ierr/=0) then
+        write(*,*) trim(rname)//" failed allocation for heightvals()"
+        call abort_physic(rname,"failled allocation",1)
+      endif
+      allocate(eruptfreqs(nvolc,2),stat=ierr)
+      if (ierr/=0) then
+        write(*,*) trim(rname)//" failed allocation for eruptfreqs()"
+        call abort_physic(rname,"failled allocation",1)
+      endif
+      allocate(volfluxes(nvolc,nq),stat=ierr)
+      if (ierr/=0) then
+        write(*,*) trim(rname)//" failed allocation for volfluxes()"
+        call abort_physic(rname,"failled allocation",1)
+      endif
+      allocate(ivolc(nvolc),stat=ierr)
+      if (ierr/=0) then
+        write(*,*) trim(rname)//" failed allocation for ivolc()"
+        call abort_physic(rname,"failled allocation",1)
+      endif
+    endif ! of if (nvolc>0) 
+    
+    ! 3.3. copy over global data to local arrays
+    do volci=1,nvolc
+      do i=1,nvolc_total
+        if (belongs_here(i)) then
+          ! copy data to local arrays
+          latlonvals(volci,1:2)=latlonvals_total(i,1:2)
+          heightvals(volci)=heightvals_total(i)
+          eruptfreqs(volci,1:2)=eruptfreqs_total(i,1:2)
+          volfluxes(volci,1:nq)=volfluxes_total(i,1:nq)
+          ivolc(volci)=volc_index(i)
+          ! the job is done, move on to next one
+          belongs_here(i)=.false.
+          exit
+        endif ! of if (belongs_here(i))
+      enddo ! of do i=1,nvolc_total
+    enddo ! of do volci=1,nvolc
+    
+    ! 4. Sanity check
+    ! verify that the sum of volcanoes over all cores is indeed nvolc_total
+    if (is_master) nvolc_check=0
+    call bcast(nvolc_check)
+    call reduce_sum(nvolc,nvolc_check)
+    call bcast(nvolc_check)
+    if (nvolc_check/=nvolc_total) then
+      write(*,*) trim(rname)//" error. Sum of volcanoes over all cores ", &
+       nvolc_check," differs from nvolc_total=",nvolc_total
+      call abort_physic(rname," error. Check and fix volcano.F90",1)
+    endif
+    
+    ! give a summary, for each core:
+    write(*,*) trim(rname)//": nvolc=",nvolc," /",nvolc_total
+    do volci=1,nvolc
+      write(*,*) trim(rname)//": latlonvals(1:2)=",latlonvals(volci,1:2)
+      write(*,*) trim(rname)//": heightvals=",heightvals(volci)
+      write(*,*) trim(rname)//": eruptfreqs(1:2)=",eruptfreqs(volci,1:2)
+      write(*,*) trim(rname)//": volfluxes(1:nq)=",volfluxes(volci,1:nq)
+      write(*,*) trim(rname)//": ivolc=",ivolc(volci)
+      write(*,*) trim(rname)//": longitude(ivolc)=",longitude_deg(ivolc(volci))
+      write(*,*) trim(rname)//": latitude(ivolc)=",latitude_deg(ivolc(volci))
+    enddo
+    
+  END SUBROUTINE read_volcano
+
+END MODULE volcano_mod
