Index: trunk/LMDZ.VENUS/libf/phyvenus/photolysis_mod.F90
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/photolysis_mod.F90	(revision 3529)
+++ trunk/LMDZ.VENUS/libf/phyvenus/photolysis_mod.F90	(revision 3530)
@@ -40,5 +40,5 @@
   real, dimension(nw), save :: xshocl                                 ! hocl absorption cross-section (cm2)
   real, dimension(nw), save :: xsso2_200, xsso2_298, xsso2_360        ! so2 absorption cross-section (cm2)
-  real, dimension(nw), save :: xsso                                   ! so absorption cross-section (cm2)
+  real, dimension(nw), save :: xsso_150, xsso_250                     ! so absorption cross-section (cm2)
   real, dimension(nw), save :: xsso3                                  ! so3 absorption cross-section (cm2)
   real, dimension(nw), save :: xsclo                                  ! clo absorption cross-section (cm2)
@@ -58,6 +58,6 @@
   real, dimension(nw), save :: albedo                                 ! surface albedo
 
-!$OMP THREADPRIVATE(f,xsco2_195, xsco2_295, xsco2_370,yieldco2,xso2_150, xso2_200, xso2_250, xso2_300,yieldo2,xso3_218, xso3_298,xsh2o,xsh2o2,xsho2,xsh2,yieldh2,xsno2, xsno2_220, xsno2_294,yldno2_248,yldno2_298)
-!$OMP THREADPRIVATE(xsno,yieldno,xsn2,yieldn2,xshdo,albedo)
+!$OMP THREADPRIVATE(f, xsco2_195, xsco2_295, xsco2_370, yieldco2, xso2_150, xso2_200, xso2_250, xso2_300, yieldo2, xso3_218, xso3_298, xsh2o, xsh2o2, xsho2, xsh2, yieldh2, xsno2, xsno2_220, xsno2_294, yldno2_248, yldno2_298)
+!$OMP THREADPRIVATE(xsno, yieldno, xsn2, yieldn2, xshdo, albedo)
 
 contains
@@ -123,5 +123,5 @@
 ! read and grid so cross-sections
 
-  call rdxsso(nw,wl,xsso)
+  call rdxsso(nw,wl,xsso_150,xsso_250)
 
 ! read and grid so3 cross-sections
@@ -2327,10 +2327,10 @@
 !==============================================================================
 
-      subroutine rdxsso(nw, wl, yg)
+      subroutine rdxsso(nw,wl,xsso_150,xsso_250)
 
 !-----------------------------------------------------------------------------*
 !=  PURPOSE:                                                                 =*
-!=  Read SO cross-sections                                                  =*
-!=  JPL 2006 recommendation                                                  =*
+!=  Read SO cross-sections                                                   =*
+!=  Phillips (1981) or Heays et al. (2023)                                   =*
 !-----------------------------------------------------------------------------*
 !=  PARAMETERS:                                                              =*
@@ -2342,53 +2342,132 @@
       USE mod_phys_lmdz_transfert_para, ONLY: bcast
 
-      IMPLICIT NONE
+      implicit none
 
 !     input
+
       integer :: nw               ! number of wavelength grid points
-      real, dimension(nw) :: wl  ! lower wavelength for each interval
+      real, dimension(nw) :: wl   ! lower wavelength for each interval
 
 !     output
 
-      real, dimension(nw) :: yg   ! SO cross-sections (cm2)
+      real, dimension(nw) :: xsso_150, xsso_250 ! so cross-sections (cm2)
 
 !     local
+
+      integer, parameter :: kdata = 1000
+      integer :: i, l, n, n1, n2, ierr, iopt
+      integer :: kin, kout ! input/output logical units
+      real, dimension(nw) :: yg, yg1, yg2
+      real, dimension(kdata) :: x1, y1, x2, y2
       real :: dummy
       real, parameter :: deltax = 1.e-4
-      integer, parameter :: kdata = 1000
-      real, dimension(kdata) :: x1, y1
-      integer :: i, n, ierr
       character*100 fil
-      integer :: kin, kout ! input/output logical units
 
       kin = 10
 
-!*** cross sections from JPL [2006]
-
-      fil = 'cross_sections/so_marcq.txt'
-      print*, 'section efficace SO: ', fil
-
-      if(is_master) then
+!     iopt = 1 : Phillips (1981)                        
+!     iopt = 2 : Heays et al., Molecular Physics, 2023 
+
+      iopt = 2
+
+      if (iopt == 1) then
+
+!        cross sections from Phillips (1981)
+
+         fil = 'cross_sections/so_marcq.txt'
+
+         if (is_master) then
+
+         print*, 'section efficace SO: ', fil
 
          n = 851
-      OPEN(kin,FILE=fil,STATUS='OLD') 
-      DO i = 1,n
-        READ(kin,*) x1(i), dummy, dummy, dummy, dummy, y1(i)
-        y1(i) = y1(i) * 0.88
-      ENDDO
-      CLOSE(kin)
-
-      CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
-      CALL addpnt(x1,y1,kdata,n,          0.,0.)
-      CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
-      CALL addpnt(x1,y1,kdata,n,        1E38,0.)
-
-      CALL inter2(nw,wl,yg,n,x1,y1,ierr)
+         OPEN(kin,FILE=fil,STATUS='OLD') 
+         DO i = 1,n
+            READ(kin,*) x1(i), dummy, dummy, dummy, dummy, y1(i)
+            y1(i) = y1(i)*0.88  ! from Mills, PhD, 1998
+         ENDDO
+         CLOSE(kin)
+
+         CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
+         CALL addpnt(x1,y1,kdata,n,          0.,0.)
+         CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
+         CALL addpnt(x1,y1,kdata,n,        1E38,0.)
+
+         do i = 1,n
+            print*, x1(i), y1(i)
+         end do
+
+         CALL inter2(nw,wl,yg,n,x1,y1,ierr)
   
-      IF (ierr .NE. 0) THEN
-        WRITE(*,*) ierr, fil
-        STOP
-      ENDIF
-
-      endif !is_master
+         IF (ierr .NE. 0) THEN
+           WRITE(*,*) ierr, fil
+           STOP
+         ENDIF
+
+         DO l = 1, nw-1
+            xsso_150(l) = yg(l)
+            xsso_250(l) = yg(l)
+         END DO
+
+         endif !is_master
+
+      else if (iopt == 2) then
+
+!        cross sections from Heays et al. (2023)
+!        values given at 150 K and 250 K
+
+         fil = 'cross_sections/SO_heays_2023.txt'
+
+         if (is_master) then
+
+         print*, 'section efficace SO: ', fil
+         
+         n1 = 525
+         n2 = n1
+         OPEN(kin,FILE=fil,STATUS='OLD') 
+
+         DO i = 1, 3
+            READ(kin,*)
+         ENDDO
+
+         DO i = 1,n1
+            READ(kin,*) x1(i), dummy, dummy, dummy, dummy, dummy, &
+                        dummy, dummy, dummy, y1(i), y2(i)
+            x2(i) = x1(i)
+         ENDDO
+         CLOSE(kin)
+
+         CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax),0.)
+         CALL addpnt(x1,y1,kdata,n1,               0.,0.)
+         CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.)
+         CALL addpnt(x1,y1,kdata,n1,           1.e+38,0.)
+
+         CALL inter2(nw,wl,yg1,n1,x1,y1,ierr)
+
+         IF (ierr .NE. 0) THEN
+            WRITE(*,*) ierr, fil
+            STOP
+         ENDIF
+
+         CALL addpnt(x2,y2,kdata,n2,x2(1)*(1.-deltax),0.)
+         CALL addpnt(x2,y2,kdata,n2,               0.,0.)
+         CALL addpnt(x2,y2,kdata,n2,x2(n2)*(1.+deltax),0.)
+         CALL addpnt(x2,y2,kdata,n2,           1.e+38,0.)
+
+         CALL inter2(nw,wl,yg2,n2,x2,y2,ierr)
+
+         IF (ierr .NE. 0) THEN
+            WRITE(*,*) ierr, fil
+            STOP
+         ENDIF
+
+         DO l = 1, nw-1
+            xsso_150(l) = yg1(l)
+            xsso_250(l) = yg2(l)
+         END DO
+
+         endif !is_master
+
+      end if ! iopt
 
       call bcast(yg)
Index: trunk/LMDZ.VENUS/libf/phyvenus/photolysis_online.F
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/photolysis_online.F	(revision 3529)
+++ trunk/LMDZ.VENUS/libf/phyvenus/photolysis_online.F	(revision 3530)
@@ -183,8 +183,14 @@
      $             dtgas(1:nlayer,:,a_h2o2), sj)
 
+!     so:
+
+      call setso(nphot, nlayer, nw, wc, temp, xsso_150, xsso_250,
+     $           j_so, colinc(1:nlayer), rm(:,i_so),
+     $           dtgas(1:nlayer,:,a_so), sj)
+
 !     so2:
 
-      call setso2(nphot, nlayer, nw, wc, temp, xsso2_200,xsso2_298,
-     $            xsso2_360, j_so2,colinc(1:nlayer), rm(:,i_so2),
+      call setso2(nphot, nlayer, nw, wc, temp, xsso2_200, xsso2_298,
+     $            xsso2_360, j_so2, colinc(1:nlayer), rm(:,i_so2),
      $            dtgas(1:nlayer,:,a_so2), sj)
 
@@ -209,5 +215,4 @@
             dtgas(ilay,iw,a_hocl) = colinc(ilay)*rm(ilay,i_hocl)
      $      *xshocl(iw)
-            dtgas(ilay,iw,a_so) = colinc(ilay)*rm(ilay,i_so)*xsso(iw)
             dtgas(ilay,iw,a_so3) = colinc(ilay)*rm(ilay,i_so3)*xsso3(iw)
             dtgas(ilay,iw,a_s2) = colinc(ilay)*rm(ilay,i_s2)*xss2(iw)
@@ -245,5 +250,4 @@
             sj(ilay,iw,j_cl2) = xscl2(iw)                ! cl2
             sj(ilay,iw,j_hocl) = xshocl(iw)              ! hocl
-            sj(ilay,iw,j_so) = xsso(iw)                  ! so 
             sj(ilay,iw,j_s2) = xss2(iw)                  ! s2
             sj(ilay,iw,j_so3) = xsso3(iw)                ! so3
@@ -647,4 +651,63 @@
 
       end subroutine seto3
+!==============================================================================
+
+      subroutine setso(nd, nlayer, nw, wc, tlay, xsso_150, xsso_250, 
+     $                 j_so, colinc, rm, dt, sj)
+
+!-----------------------------------------------------------------------------*
+!=  PURPOSE:                                                                 =*
+!=  Set up the so temperature dependent cross-sections and optical depth     =*
+!-----------------------------------------------------------------------------*
+
+      implicit none
+
+!     input:
+
+      integer :: nd                                               ! number of photolysis rates
+      integer :: nlayer                                           ! number of vertical layers
+      integer :: nw                                               ! number of wavelength grid points
+      integer :: j_so                                             ! photolysis indexe
+      real, dimension(nw)     :: wc                               ! central wavelength for each interval
+      real, dimension(nw)     :: xsso_150, xsso_250               ! so cross-sections (cm2)
+      real, dimension(nlayer) :: tlay                             ! temperature (k)
+      real, dimension(nlayer) :: rm                               ! so mixing ratio
+      real, dimension(nlayer) :: colinc                           ! air column increment (molecule.cm-2) 
+
+!     output:
+
+      real, dimension(nlayer,nw)    :: dt                         !  optical depth
+      real, dimension(nlayer,nw,nd) :: sj                         !  cross-section array (cm2)
+
+!     local:
+!
+      integer :: ilev, iw
+      real :: temp, xsso
+
+!     calculate temperature dependance
+
+      do ilev = 1,nlayer
+         temp = max(tlay(ilev),150.)
+         temp = min(temp, 250.)
+         do iw = 1, nw-1
+            if (xsso_150(iw) /= 0. .and. xsso_250(iw) /= 0.) then
+               xsso = log(xsso_150(iw)) 
+     $              + (log(xsso_250(iw)) - log(xsso_150(iw)))
+     $               /(250. - 150.)*(temp - 150.) 
+
+               sj(ilev,iw,j_so) = exp(xsso)
+            else
+               sj(ilev,iw,j_so) = 0.
+            end if
+
+!     optical depth
+
+            dt(ilev,iw) = colinc(ilev)*rm(ilev)*sj(ilev,iw,j_so)
+
+         end do
+      end do
+
+      end subroutine setso
+
 !==============================================================================
 
