Index: trunk/LMDZ.MARS/libf/aeronomars/calchim.F90
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/calchim.F90	(revision 2024)
+++ trunk/LMDZ.MARS/libf/aeronomars/calchim.F90	(revision 2027)
@@ -91,5 +91,5 @@
       real :: pdq(ngrid,nlayer,nq)   ! previous tendencies
       real :: zday                   ! date (time since Ls=0, in martian days)
-      real :: tauref(ngrid)          ! optical depth at 7 hPa
+      real :: tauref(ngrid)          ! Reference column optical depth
       real :: co2ice(ngrid)          ! co2 ice surface layer (kg.m-2)
       real :: surfdust(ngrid,nlayer) ! dust surface area (m2/m3)
@@ -150,5 +150,5 @@
       real    :: zu(ngrid,nlayer)      ! u component of the wind
       real    :: zv(ngrid,nlayer)      ! v component of the wind
-      real    :: taucol                ! optical depth at 7 hPa
+      real    :: taucol                ! dust optical depth at the surface
 
       logical,save :: firstcall = .true.
@@ -659,7 +659,7 @@
 
             call photochemistry(nlayer,nq,                            &
-                                ig,lswitch,zycol,szacol,ptimestep,    &
-                                zpress,ztemp,zdens,zmmean,dist_sol,   &
-                                surfdust1d,surfice1d,jo3,taucol,iter)
+                           ig,lswitch,zycol,szacol,ptimestep,         &
+                           zpress,zlocal,ztemp,zdens,zmmean,dist_sol, &
+                           surfdust1d,surfice1d,jo3,taucol,iter)
 
 !        ozone photolysis, for output
Index: trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F90
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F90	(revision 2024)
+++ trunk/LMDZ.MARS/libf/aeronomars/photochemistry.F90	(revision 2027)
@@ -13,7 +13,7 @@
 !*****************************************************************
 
-subroutine photochemistry(nlayer, nq,                                 &
-                          ig, lswitch, zycol, sza, ptimestep, press,  &
-                          temp, dens, zmmean, dist_sol, surfdust1d,   &
+subroutine photochemistry(nlayer, nq,                                   &
+                          ig, lswitch, zycol, sza, ptimestep, press,    &
+                          alt, temp, dens, zmmean, dist_sol, surfdust1d,&
                           surfice1d, jo3, tau, iter)
 
@@ -34,4 +34,5 @@
 real :: ptimestep             ! physics timestep (s)
 real :: press(nlayer)         ! pressure (hpa)
+real :: alt(nlayer)           ! altitude (km)
 real :: temp(nlayer)          ! temperature (k)
 real :: dens(nlayer)          ! density (cm-3)
@@ -40,5 +41,5 @@
 real :: surfdust1d(nlayer)    ! dust surface area (cm2/cm3)
 real :: surfice1d(nlayer)     ! ice surface area (cm2/cm3)
-real :: tau                   ! optical depth at 7 hpa
+real :: tau                   ! reference dust optical depth
 
 !===================================================================
@@ -65,4 +66,5 @@
 
 logical, save :: firstcall = .true.
+logical :: jonline
 
 parameter (nesp = 17)         ! number of species in the chemistry
@@ -141,6 +143,17 @@
 !===================================================================
 
-call photolysis(nlayer, lswitch, press, temp, sza, tau, zmmean, dist_sol, & 
-                rm(:,i_co2), rm(:,i_o3), v_phot)
+jonline = .false.  ! .false. recommended until the end of testing phase
+
+if (jonline) then
+   tau = tau*press(1)/6.1  ! temporary
+   call photolysis_online(nlayer,lswitch, alt, press, temp,                  &
+                          zmmean, i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,  &
+                          i_h2, i_oh, i_ho2, i_h2o2, i_h2o,                  &
+                          i_n, i_n2d, i_no, i_no2, i_n2, nesp, rm,           &
+                          tau, sza, dist_sol, v_phot)
+else
+   call photolysis(nlayer, lswitch, press, temp, sza, tau, zmmean, dist_sol, & 
+                   rm(:,i_co2), rm(:,i_o3), v_phot)
+end if
 
 ! save o3 photolysis for output
Index: trunk/LMDZ.MARS/libf/aeronomars/photolysis_online.F
===================================================================
--- trunk/LMDZ.MARS/libf/aeronomars/photolysis_online.F	(revision 2027)
+++ trunk/LMDZ.MARS/libf/aeronomars/photolysis_online.F	(revision 2027)
@@ -0,0 +1,3391 @@
+!==============================================================================
+
+      subroutine photolysis_online(nlayer,lswitch, alt, press, temp,
+     $           zmmean, i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,
+     $           i_h2, i_oh, i_ho2, i_h2o2, i_h2o,               
+     $           i_n, i_n2d, i_no, i_no2, i_n2, nesp, rm,
+     $           tau, sza, dist_sol, v_phot)
+
+      implicit none
+
+#include "chimiedata.h"
+
+!     input 
+
+      integer :: nesp                                                ! total number of chemical species
+      integer :: nlayer, lswitch
+      integer :: i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,           
+     $           i_h2, i_oh, i_ho2, i_h2o2, i_h2o,                 
+     $           i_n, i_n2d, i_no, i_no2, i_n2
+
+      real, dimension(nlayer) :: press, temp, zmmean                 ! pressure (hpa)/temperature (k)/mean molecular mass (g.mol-1)
+      real, dimension(nlayer) :: alt                                 ! altitude (km)
+      real, dimension(nlayer,nesp) :: rm                             ! mixing ratios 
+      real :: tau                                                    ! integrated aerosol optical depth at the surface
+      real :: sza                                                    ! solar zenith angle (degrees)
+      real :: dist_sol                                               ! solar distance (au)
+
+!     output
+
+      real (kind = 8), dimension(nlayer,nb_phot_max) :: v_phot       ! photolysis rates (s-1)
+ 
+!     spectral grid
+
+!     integer, parameter :: nw = 3789                                ! number of spectral intervals (high-res)
+      integer, parameter :: nw = 152                                 ! number of spectral intervals (low-res)
+      real, dimension(nw) :: wl, wc, wu                              ! lower, center, upper wavelength for each interval
+      save wl, wc, wu
+
+!     solar flux
+
+      real, dimension(nw) :: f                                       !  solar flux (w.m-2.nm-1) at 1 au
+      real :: factor
+      save f
+
+!     cross-sections and quantum yields
+
+      integer, parameter  :: nabs = 10                               ! number of absorbing gases
+      real, dimension(nlayer,nw,nd) :: sj                            ! general cross-section array (cm2)
+      real, dimension(nw) :: xsco2_195, xsco2_295, xsco2_370         ! co2 absorption cross-section at 195-295-370 k (cm2)
+      real, dimension(nw) :: yieldco2                                ! co2 photodissociation yield
+      real, dimension(nw) :: xso2_150, xso2_200, xso2_250, xso2_300  ! o2 absorption cross-section at 150-200-250-300 k (cm2)
+      real, dimension(nw) :: yieldo2                                 ! o2 photodissociation yield
+      real, dimension(nw) :: xso3_218, xso3_298                      ! o3 absorption cross-section at 218-298 k (cm2)
+      real, dimension(nw) :: xsh2o                                   ! h2o absorption cross-section (cm2)
+      real, dimension(nw) :: xsh2o2                                  ! h2o2 absorption cross-section (cm2)
+      real, dimension(nw) :: xsho2                                   ! ho2 absorption cross-section (cm2)
+      real, dimension(nw) :: xsh2                                    ! h2 absorption cross-section (cm2)
+      real, dimension(nw) :: yieldh2                                 ! h2 photodissociation yield
+      real, dimension(nw) :: xsno2, xsno2_220, xsno2_294             ! no2 absorption cross-section at 220-294 k (cm2)
+      real, dimension(nw) :: yldno2_248, yldno2_298                  ! no2 quantum yield at 248-298 k
+      real, dimension(nw) :: xsno                                    ! no absorption cross-section (cm2)
+      real, dimension(nw) :: yieldno                                 ! no photodissociation yield
+      real, dimension(nw) :: xsn2                                    ! n2 absorption cross-section (cm2)
+      real, dimension(nw) :: yieldn2                                 ! n2 photodissociation yield
+      real, dimension(nw) :: albedo                                  ! surface albedo
+
+!     atmosphere
+
+      real, dimension(nlayer) :: colinc                              ! air column increment (molecule.cm-2) 
+      real, dimension(nlayer,nw) :: dtrl                             ! rayleigh optical depth 
+      real, dimension(nlayer,nw) :: dtaer                            ! aerosol optical depth 
+      real, dimension(nlayer,nw) :: omaer                            ! aerosol single scattering albedo 
+      real, dimension(nlayer,nw) :: gaer                             ! aerosol asymmetry parameter
+      real, dimension(nlayer,nw) :: dtcld                            ! cloud optical depth 
+      real, dimension(nlayer,nw) :: omcld                            ! cloud single scattering albedo 
+      real, dimension(nlayer,nw) :: gcld                             ! cloud asymmetry parameter
+      real, dimension(nlayer,nw,nabs) :: dtgas                       ! optical depth for each gas 
+      real, dimension(nlayer,nw) :: dagas                            ! total gas optical depth
+      real, dimension(nlayer) :: edir, edn, eup                      ! normalised irradiances
+      real, dimension(nlayer) :: fdir, fdn, fup                      ! normalised actinic fluxes
+      real, dimension(nlayer) :: saflux                              ! total actinic flux
+
+      integer, dimension(0:nlayer) :: nid
+      real, dimension(0:nlayer,nlayer) :: dsdh
+
+      integer :: j_o2_o, j_o2_o1d, j_co2_o, j_co2_o1d, j_o3_o1d, j_o3_o,
+     $           j_h2o, j_h2o2, j_ho2, j_h2, j_no, j_no2, j_n2
+
+      integer :: a_o2, a_co2, a_o3, a_h2o, a_h2o2, a_ho2, a_h2, a_no, 
+     $           a_no2, a_n2
+      integer :: i, ilay, iw, ialt, mopt
+      real    :: deltaj
+
+      logical, save :: firstcall
+      data firstcall / .true. /
+
+!     high-res/low-res switch
+!
+!     mopt = 1 high-resolution
+!     mopt = 2 low-resolution (recommended for gcm use)
+
+      if (nw >= 3789) then
+         mopt = 1
+      else
+         mopt = 2
+      end if
+
+!     absorbing gas numbering
+
+      a_o2      = 1      ! o2
+      a_co2     = 2      ! co2
+      a_o3      = 3      ! o3
+      a_h2o     = 4      ! h2o
+      a_h2o2    = 5      ! h2o2
+      a_ho2     = 6      ! ho2
+      a_h2      = 7      ! h2
+      a_no      = 8      ! no
+      a_no2     = 9      ! no2
+      a_n2      = 10     ! no2
+
+!     photodissociation rates numbering
+
+      j_o2_o    = 1      ! o2 + hv     -> o + o
+      j_o2_o1d  = 2      ! o2 + hv     -> o + o(1d)
+      j_co2_o   = 3      ! co2 + hv    -> co + o
+      j_co2_o1d = 4      ! co2 + hv    -> co + o(1d)
+      j_o3_o1d  = 5      ! o3 + hv     -> o2 + o(1d)
+      j_o3_o    = 6      ! o3 + hv     -> o2 + o
+      j_h2o     = 7      ! h2o + hv    -> h + oh
+      j_h2o2    = 8      ! h2o2 + hv   -> oh + oh
+      j_ho2     = 9      ! ho2 + hv    -> oh + o
+      j_h2      = 10     ! h2 + hv     -> h + h
+      j_no      = 11     ! no + hv     -> n + o
+      j_no2     = 12     ! no2 + hv    -> no + o
+      j_n2      = 13     ! n2 + hv     -> n + n
+
+!===== day/night criterion
+
+      if (sza <= 120.) then
+
+!===== read once absorption cross-sections, quantum yields, and albedo
+
+      if (firstcall) then
+
+!        set wavelength grid
+
+         call gridw(nw,wl,wc,wu,mopt)
+
+!        read and grid solar flux data
+ 
+         call rdsolarflux(nw,wl,wc,f)
+
+!        read and grid o2 cross-sections
+
+         call rdxso2(nw,wl,xso2_150,xso2_200,xso2_250,xso2_300,yieldo2)
+
+!        read and grid co2 cross-sections
+ 
+         call rdxsco2(nw,wl,xsco2_195,xsco2_295,xsco2_370,yieldco2)
+
+!        read and grid o3 cross-sections
+
+         call rdxso3(nw,wl,xso3_218,xso3_298)
+
+!        read and grid h2o cross-sections
+
+         call rdxsh2o(nw,wl,xsh2o)
+
+!        read and grid h2o2 cross-sections
+
+         call rdxsh2o2(nw,wl,xsh2o2)
+
+!        read and grid ho2 cross-sections
+
+         call rdxsho2(nw,wl,xsho2)
+
+!        read and grid h2 cross-sections
+
+         call rdxsh2(nw,wl,wc,xsh2,yieldh2)
+
+!        read and grid no cross-sections
+
+         call rdxsno(nw,wl,xsno,yieldno)
+
+!        read and grid no2 cross-sections
+
+         call rdxsno2(nw,wl,xsno2,xsno2_220,xsno2_294,
+     $                yldno2_248, yldno2_298)
+
+!        read and grid n2 cross-sections
+
+         call rdxsn2(nw,wl,xsn2,yieldn2)
+
+!        set surface albedo
+
+         call setalb(nw,wl,albedo)
+
+         firstcall = .false.
+
+      end if
+
+!==== air column increments and rayleigh optical depth
+
+      call setair(nlayer, nw, wl, wc, press, temp, zmmean, colinc, dtrl)
+
+!==== set temperature-dependent cross-sections and optical depths
+
+!     o2:
+
+      call seto2(nd, nlayer, nw, wc, mopt, temp, xso2_150, xso2_200, 
+     $           xso2_250, xso2_300, yieldo2, j_o2_o, j_o2_o1d,
+     $           colinc, rm(:,i_o2), dtgas(:,:,a_o2), sj)
+
+!     co2:
+
+      call setco2(nd, nlayer, nw, wc, temp, xsco2_195, xsco2_295, 
+     $            xsco2_370, yieldco2, j_co2_o, j_co2_o1d,
+     $            colinc, rm(:,i_co2), dtgas(:,:,a_co2), sj)
+
+!     o3:
+
+      call seto3(nd, nlayer, nw, wc, temp, xso3_218, xso3_298, 
+     $           j_o3_o, j_o3_o1d,
+     $           colinc, rm(:,i_o3), dtgas(:,:,a_o3), sj)
+
+!     h2o2:
+
+      call seth2o2(nd, nlayer, nw, wc, temp, xsh2o2, j_h2o2,
+     $             colinc, rm(:,i_h2o2), dtgas(:,:,a_h2o2), sj)
+
+!     no2:
+
+      call setno2(nd, nlayer, nw, wc, temp, xsno2, xsno2_220,
+     $            xsno2_294, yldno2_248, yldno2_298, j_no2,
+     $            colinc, rm(:,i_no2), dtgas(:,:,a_no2), sj)
+
+!==== temperature independent optical depths and cross-sections
+
+!     optical depths
+
+      do ilay = 1,nlayer
+         do iw = 1,nw-1
+            dtgas(ilay,iw,a_h2o) = colinc(ilay)*rm(ilay,i_h2o)*xsh2o(iw)
+            dtgas(ilay,iw,a_ho2) = colinc(ilay)*rm(ilay,i_ho2)*xsho2(iw)
+            dtgas(ilay,iw,a_h2)  = colinc(ilay)*rm(ilay,i_h2)*xsh2(iw)
+            dtgas(ilay,iw,a_no)  = colinc(ilay)*rm(ilay,i_no)*xsno(iw)
+            dtgas(ilay,iw,a_n2)  = colinc(ilay)*rm(ilay,i_n2)*xsn2(iw)
+         end do
+      end do
+
+!     total gas optical depth
+
+      dagas(:,:) = 0.
+
+      do ilay = 1,nlayer
+         do iw = 1,nw-1
+            do i = 1,nabs
+               dagas(ilay,iw) = dagas(ilay,iw) + dtgas(ilay,iw,i)
+            end do
+         end do
+      end do
+
+!     cross-sections
+
+      do ilay = 1,nlayer
+         do iw = 1,nw-1
+            sj(ilay,iw,j_h2o) = xsh2o(iw)                ! h2o
+            sj(ilay,iw,j_ho2) = xsho2(iw)                ! ho2
+            sj(ilay,iw,j_h2)  = xsh2(iw)*yieldh2(iw)     ! h2
+            sj(ilay,iw,j_no)  = xsno(iw)*yieldno(iw)     ! no
+            sj(ilay,iw,j_n2)  = xsn2(iw)*yieldn2(iw)     ! n2
+         end do
+      end do
+
+!     ialt = 49
+!     print*, 'altitude = ', alt(ialt), ' km'
+!     do iw = 1,nw
+!        write(35,'(f8.4,12e13.5)') wc(iw), dtrl(ialt,iw), 
+!    $                      (max(dtgas(ialt,iw,i),1.e-30),i = 1,nabs)
+!     end do
+
+!==== set aerosol properties and optical depth
+
+      call setaer(nlayer,alt,tau,nw,dtaer,omaer,gaer)
+
+!==== set cloud properties and optical depth
+
+      call setcld(nlayer,nw,dtcld,omcld,gcld)
+
+!==== slant path lengths in spherical geometry
+
+      call sphers(nlayer, alt, sza, dsdh, nid)
+
+!==== solar flux at mars
+
+      factor = (1./dist_sol)**2.
+      do iw = 1,nw-1
+         f(iw) = f(iw)*factor
+      end do
+
+!==== initialise photolysis rates
+
+      v_phot(:,1:nd) = 0.
+
+!==== start of wavelength lopp
+
+      do iw = 1,nw-1
+
+!     monochromatic radiative transfer. outputs are:
+!     normalized irradiances     edir(nlayer), edn(nlayer), eup(nlayer) 
+!     normalized actinic fluxes  fdir(nlayer), fdn(nlayer), fup(nlayer)
+!     where 
+!     dir = direct beam, dn = down-welling diffuse, up = up-welling diffuse
+
+         call rtlink(nlayer, nw, iw, albedo(iw), sza, dsdh, nid, dtrl,
+     $               dagas, dtcld, omcld, gcld, dtaer, omaer, gaer,
+     $               edir, edn, eup, fdir, fdn, fup)
+
+
+!     spherical actinic flux
+
+         do ilay = 1,lswitch-1
+            saflux(ilay) = f(iw)*(fdir(ilay) + fdn(ilay) + fup(ilay))
+         end do
+
+!     photolysis rate integration
+
+         do i = 1,nd
+            do ilay = 1, lswitch-1
+               deltaj = saflux(ilay)*sj(ilay,iw,i)
+               v_phot(ilay,i) = v_phot(ilay,i) + deltaj*(wu(iw)-wl(iw))
+            end do
+         end do
+
+!     eliminate small values
+
+         where (v_phot(:,1:nd) < 1.e-30)
+            v_phot(:,1:nd) = 0.
+         end where
+
+      end do ! iw
+
+      else   ! night
+
+         v_phot(:,1:nd) = 0.
+
+      end if ! sza
+
+      end subroutine photolysis_online
+
+!==============================================================================
+
+      subroutine gridw(nw,wl,wc,wu,mopt)
+
+!     Create the wavelength grid for all interpolations and radiative transfer 
+!     calculations.  Grid may be irregularly spaced.  Wavelengths are in nm.  
+!     No gaps are allowed within the wavelength grid.                        
+
+      implicit none
+
+!     input
+ 
+      integer :: nw      ! number of wavelength grid points
+      integer :: mopt    ! high-res/low-res switch
+
+!     output
+
+      real, dimension(nw) :: wl, wc, wu   ! lower, center, upper wavelength for each interval
+
+!     local
+
+      real :: wincr    ! wavelength increment
+      integer :: iw, kw
+
+!     mopt = 1    high-resolution mode (3789 intervals)
+!
+!                   0-108 nm :  1.0  nm
+!                 108-124 nm :  0.1  nm
+!                 124-175 nm :  0.5  nm
+!                 175-205 nm :  0.01 nm
+!                 205-365 nm :  0.5  nm
+!                 365-850 nm :  5.0  nm  
+!
+!     mopt = 2    low-resolution mode
+!
+!                    0-60 nm :  6.0 nm
+!                   60-80 nm :  2.0 nm
+!                  80-120 nm :  5.0 nm
+!                 120-123 nm :  0.2 nm
+!                 123-163 nm :  5.0 nm
+!                 163-175 nm :  2.0 nm
+!                 175-205 nm :  0.5 nm
+!                 205-245 nm :  5.0 nm
+!                 245-415 nm : 10.0 nm
+!                 415-815 nm : 50.0 nm
+
+      if (mopt == 1) then   ! high-res
+
+! define wavelength intervals of width 1.0 nm from 0 to 108 nm:
+
+      kw = 0
+      wincr = 1.0
+      do iw = 0, 107
+        kw = kw + 1
+        wl(kw) = real(iw)
+        wu(kw) = wl(kw) + wincr
+        wc(kw) = (wl(kw) + wu(kw))/2.
+      end do
+
+! define wavelength intervals of width 0.1 nm from 108 to 124 nm:
+
+      wincr = 0.1
+      do iw = 1080, 1239, 1
+        kw = kw + 1
+        wl(kw) = real(iw)/10.
+        wu(kw) = wl(kw) + wincr
+        wc(kw) = (wl(kw) + wu(kw))/2.
+      end do
+
+! define wavelength intervals of width 0.5 nm from 124 to 175 nm:
+
+      wincr = 0.5
+      do iw = 1240, 1745, 5
+        kw = kw + 1
+        wl(kw) = real(iw)/10.
+        wu(kw) = wl(kw) + wincr
+        wc(kw) = (wl(kw) + wu(kw))/2.
+      end do
+
+! define wavelength intervals of width 0.01 nm from 175 to 205 nm:
+
+      wincr = 0.01
+      do iw = 17500, 20499, 1
+        kw = kw + 1
+        wl(kw) = real(iw)/100.
+        wu(kw) = wl(kw) + wincr
+        wc(kw) = (wl(kw) + wu(kw))/2.
+      end do
+
+! define wavelength intervals of width 0.5 nm from 205 to 365 nm:
+
+      wincr = 0.5
+      do iw = 2050, 3645, 5
+        kw = kw + 1
+        wl(kw) = real(iw)/10.
+        wu(kw) = wl(kw) + wincr
+        wc(kw) = (wl(kw) + wu(kw))/2.
+      end do
+
+! define wavelength intervals of width 5.0 nm from 365 to 855 nm:
+
+      wincr = 5.0
+      do iw = 365, 850, 5
+        kw = kw + 1
+        wl(kw) = real(iw)
+        wu(kw) = wl(kw) + wincr
+        wc(kw) = (wl(kw) + wu(kw))/2.
+      end do
+      wl(kw+1) = wu(kw)
+
+      else if (mopt == 2) then   ! low-res
+
+* define wavelength intervals of width 6.0 nm from 0 to 60 nm:
+
+      kw = 0
+      wincr = 6.0
+      DO iw = 0, 54, 6
+        kw = kw + 1
+        wl(kw) = real(iw)
+        wu(kw) = wl(kw) + wincr
+        wc(kw) = (wl(kw) + wu(kw))/2.
+      END DO
+
+* define wavelength intervals of width 2.0 nm from 60 to 80 nm:
+
+      wincr = 2.0
+      DO iw = 60, 78, 2
+        kw = kw + 1
+        wl(kw) = real(iw)
+        wu(kw) = wl(kw) + wincr
+        wc(kw) = (wl(kw) + wu(kw))/2.
+      END DO
+
+* define wavelength intervals of width 5.0 nm from 80 to 120 nm:
+
+      wincr = 5.0
+      DO iw = 80, 115, 5
+        kw = kw + 1
+        wl(kw) = real(iw)
+        wu(kw) = wl(kw) + wincr
+        wc(kw) = (wl(kw) + wu(kw))/2.
+      END DO
+
+
+* define wavelength intervals of width 0.2 nm from 120 to 123 nm:
+
+      wincr = 0.2
+      DO iw = 1200, 1228, 2
+        kw = kw + 1
+        wl(kw) = real(iw)/10.
+        wu(kw) = wl(kw) + wincr
+        wc(kw) = (wl(kw) + wu(kw))/2.
+      ENDDO
+
+* define wavelength intervals of width 5.0 nm from 123 to 163 nm:
+
+      wincr = 5.0
+      DO iw = 123, 158, 5
+        kw = kw + 1
+        wl(kw) = real(iw)
+        wu(kw) = wl(kw) + wincr
+        wc(kw) = (wl(kw) + wu(kw))/2.
+      ENDDO
+
+* define wavelength intervals of width 2.0 nm from 163 to 175 nm:
+
+      wincr = 2.0
+      DO iw = 163, 173, 2
+        kw = kw + 1
+        wl(kw) = real(iw)
+        wu(kw) = wl(kw) + wincr
+        wc(kw) = (wl(kw) + wu(kw))/2.
+      ENDDO
+
+* define wavelength intervals of width 0.5 nm from 175 to 205 nm:
+
+      wincr = 0.5
+      DO iw = 1750, 2045, 5
+        kw = kw + 1
+        wl(kw) = real(iw)/10.
+        wu(kw) = wl(kw) + wincr
+        wc(kw) = (wl(kw) + wu(kw))/2.
+      ENDDO
+
+* define wavelength intervals of width 5.0 nm from 205 to 245 nm:
+
+      wincr = 5.
+      DO iw = 205, 240, 5
+        kw = kw + 1
+        wl(kw) = real(iw)
+        wu(kw) = wl(kw) + wincr
+        wc(kw) = (wl(kw) + wu(kw))/2.
+      ENDDO
+
+* define wavelength intervals of width 10.0 nm from 245 to 415 nm:
+
+      wincr = 10.0
+      DO iw = 245, 405, 10
+        kw = kw + 1
+        wl(kw) = real(iw)
+        wu(kw) = wl(kw) + wincr
+        wc(kw) = (wl(kw) + wu(kw))/2.
+      ENDDO
+
+* define wavelength intervals of width 50.0 nm from 415 to 815 nm:
+
+      wincr = 50.0
+      DO iw = 415, 815, 50
+        kw = kw + 1
+        wl(kw) = real(iw)
+        wu(kw) = wl(kw) + wincr
+        wc(kw) = (wl(kw) + wu(kw))/2.
+      ENDDO
+
+      wl(kw+1) = wu(kw)
+
+      end if  ! mopt
+      
+!     do iw = 1,nw
+!        write(20,*) iw, wl(iw), wu(iw)
+!     end do
+      print*, 'number of spectral intervals : ', kw+1
+
+      return
+      end
+
+!==============================================================================
+ 
+      subroutine rdsolarflux(nw,wl,wc,f)
+
+!     Read and re-grid solar flux data.            
+
+      use datafile_mod, only: datadir
+
+      implicit none
+
+!     input
+ 
+      integer :: nw      ! number of wavelength grid points
+      real, dimension(nw) :: wl, wc   ! lower and central wavelength for each interval
+
+!     output
+
+      real, dimension(nw) :: f  ! solar flux (w.m-2.nm-1)
+
+!     local 
+
+      integer, parameter :: kdata = 20000    ! max dimension of input solar flux
+      integer :: msun   ! choice of solar flux
+      integer :: iw, nhead, ihead, n, i, ierr, kin
+
+      real, parameter :: deltax = 1.e-4
+      real, dimension(kdata) :: x1, y1      ! input solar flux
+      real, dimension(nw)    :: yg1         ! gridded solar flux
+
+      character(len=100) :: fil
+
+      kin = 10    ! input logical unit
+
+! select desired extra-terrestrial solar irradiance, using msun:
+
+! 18 = atlas3_thuillier_tuv.txt  0-900 nm  November 1994
+!      Thuillier et al., Adv. Space. Res., 34, 256-261, 2004 
+ 
+      msun = 18
+
+      if (msun == 18) THEN
+
+         fil = trim(datadir)//'/solar_fluxes/atlas3_thuillier_tuv.txt'
+         print*, 'solar flux : ', fil
+         OPEN(UNIT=kin,FILE=fil,STATUS='old')
+         nhead = 9
+         n = 19193
+         DO ihead = 1, nhead
+            READ(kin,*)
+         ENDDO
+         DO i = 1, n
+            READ(kin,*) x1(i), y1(i)
+            y1(i) = y1(i)*1.e-3    ! mw -> w
+         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,      1.e+38,0.)
+         CALL inter2(nw,wl,yg1,n,x1,y1,ierr)
+
+         IF (ierr .NE. 0) THEN
+            WRITE(*,*) ierr, fil
+            STOP
+         ENDIF
+
+!     convert to photon.s-1.nm-1.cm-2
+!     5.039e11 = 1.e-4*1e-9/(hc = 6.62e-34*2.998e8)
+
+         DO iw = 1, nw-1
+            f(iw) = yg1(iw)*wc(iw)*5.039e11
+!           write(25,*) iw, wc(iw), f(iw)
+         ENDDO
+
+      end if
+
+      end subroutine rdsolarflux
+
+!==============================================================================
+
+      subroutine addpnt ( x, y, ld, n, xnew, ynew )
+
+*-----------------------------------------------------------------------------*
+*=  PURPOSE:                                                                 =*
+*=  Add a point <xnew,ynew> to a set of data pairs <x,y>.  x must be in      =*
+*=  ascending order                                                          =*
+*-----------------------------------------------------------------------------*
+*=  PARAMETERS:                                                              =*
+*=  X    - REAL vector of length LD, x-coordinates                       (IO)=*
+*=  Y    - REAL vector of length LD, y-values                            (IO)=*
+*=  LD   - INTEGER, dimension of X, Y exactly as declared in the calling  (I)=*
+*=         program                                                           =*
+*=  N    - INTEGER, number of elements in X, Y.  On entry, it must be:   (IO)=*
+*=         N < LD.  On exit, N is incremented by 1.                          =*
+*=  XNEW - REAL, x-coordinate at which point is to be added               (I)=*
+*=  YNEW - REAL, y-value of point to be added                             (I)=*
+*-----------------------------------------------------------------------------*
+
+      IMPLICIT NONE
+
+C calling parameters
+
+      INTEGER ld, n
+      REAL x(ld), y(ld)
+      REAL xnew, ynew
+      INTEGER ierr
+
+C local variables
+
+      INTEGER insert
+      INTEGER i
+
+C-----------------------------------------------------------------------
+
+* initialize error flag
+
+      ierr = 0
+
+* check n<ld to make sure x will hold another point
+
+      IF (n .GE. ld) THEN
+         WRITE(0,*) '>>> ERROR (ADDPNT) <<<  Cannot expand array '
+         WRITE(0,*) '                        All elements used.'
+         STOP
+      ENDIF
+
+      insert = 1
+      i = 2
+
+* check, whether x is already sorted.
+* also, use this loop to find the point at which xnew needs to be inserted
+* into vector x, if x is sorted.
+
+ 10   CONTINUE
+      IF (i .LT. n) THEN
+        IF (x(i) .LT. x(i-1)) THEN
+           print*, x(i-1), x(i)
+           WRITE(0,*) '>>> ERROR (ADDPNT) <<<  x-data must be '//
+     >                'in ascending order!'
+           STOP
+        ELSE
+           IF (xnew .GT. x(i)) insert = i + 1
+        ENDIF
+        i = i+1
+        GOTO 10
+      ENDIF
+
+* if <xnew,ynew> needs to be appended at the end, just do so,
+* otherwise, insert <xnew,ynew> at position INSERT
+
+      IF ( xnew .GT. x(n) ) THEN
+ 
+         x(n+1) = xnew
+         y(n+1) = ynew
+  
+      ELSE
+
+* shift all existing points one index up
+
+         DO i = n, insert, -1
+           x(i+1) = x(i)
+           y(i+1) = y(i)
+         ENDDO
+
+* insert new point
+
+         x(insert) = xnew
+         y(insert) = ynew
+  
+      ENDIF
+
+* increase total number of elements in x, y
+
+      n = n+1
+
+      END
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      subroutine inter2(ng,xg,yg,n,x,y,ierr)
+
+*-----------------------------------------------------------------------------*
+*=  PURPOSE:                                                                 =*
+*=  Map input data given on single, discrete points onto a set of target     =*
+*=  bins.                                                                    =*
+*=  The original input data are given on single, discrete points of an       =*
+*=  arbitrary grid and are being linearly interpolated onto a specified set  =*
+*=  of target bins.  In general, this is the case for most of the weighting  =*
+*=  functions (action spectra, molecular cross section, and quantum yield    =*
+*=  data), which have to be matched onto the specified wavelength intervals. =*
+*=  The average value in each target bin is found by averaging the trapezoi- =*
+*=  dal area underneath the input data curve (constructed by linearly connec-=*
+*=  ting the discrete input values).                                         =*
+*=  Some caution should be used near the endpoints of the grids.  If the     =*
+*=  input data set does not span the range of the target grid, an error      =*
+*=  message is printed and the execution is stopped, as extrapolation of the =*
+*=  data is not permitted.                                                   =*
+*=  If the input data does not encompass the target grid, use ADDPNT to      =*
+*=  expand the input array.                                                  =*
+*-----------------------------------------------------------------------------*
+*=  PARAMETERS:                                                              =*
+*=  NG  - INTEGER, number of bins + 1 in the target grid                  (I)=*
+*=  XG  - REAL, target grid (e.g., wavelength grid);  bin i is defined    (I)=*
+*=        as [XG(i),XG(i+1)] (i = 1..NG-1)                                   =*
+*=  YG  - REAL, y-data re-gridded onto XG, YG(i) specifies the value for  (O)=*
+*=        bin i (i = 1..NG-1)                                                =*
+*=  N   - INTEGER, number of points in input grid                         (I)=*
+*=  X   - REAL, grid on which input data are defined                      (I)=*
+*=  Y   - REAL, input y-data                                              (I)=*
+*-----------------------------------------------------------------------------*
+
+      IMPLICIT NONE
+
+* input:
+      INTEGER ng, n
+      REAL x(n), y(n), xg(ng)
+
+* output:
+      REAL yg(ng)
+
+* local:
+      REAL area, xgl, xgu
+      REAL darea, slope
+      REAL a1, a2, b1, b2
+      INTEGER ngintv
+      INTEGER i, k, jstart
+      INTEGER ierr
+*_______________________________________________________________________
+
+      ierr = 0
+
+*  test for correct ordering of data, by increasing value of x
+
+      DO 10, i = 2, n
+         IF (x(i) .LE. x(i-1)) THEN
+            ierr = 1
+            WRITE(*,*)'data not sorted'
+            WRITE(*,*) x(i), x(i-1)
+            RETURN
+         ENDIF
+   10 CONTINUE     
+
+      DO i = 2, ng
+        IF (xg(i) .LE. xg(i-1)) THEN
+           ierr = 2
+          WRITE(0,*) '>>> ERROR (inter2) <<<  xg-grid not sorted!'
+          RETURN
+        ENDIF
+      ENDDO
+
+* check for xg-values outside the x-range
+
+      IF ( (x(1) .GT. xg(1)) .OR. (x(n) .LT. xg(ng)) ) THEN
+          WRITE(0,*) '>>> ERROR (inter2) <<<  Data do not span '//
+     >               'grid.  '
+          WRITE(0,*) '                        Use ADDPNT to '//
+     >               'expand data and re-run.'
+          STOP
+      ENDIF
+
+*  find the integral of each grid interval and use this to 
+*  calculate the average y value for the interval      
+*  xgl and xgu are the lower and upper limits of the grid interval
+
+      jstart = 1
+      ngintv = ng - 1
+      DO 50, i = 1,ngintv
+
+* initalize:
+
+            area = 0.0
+            xgl = xg(i)
+            xgu = xg(i+1)
+
+*  discard data before the first grid interval and after the 
+*  last grid interval
+*  for internal grid intervals, start calculating area by interpolating
+*  between the last point which lies in the previous interval and the
+*  first point inside the current interval
+
+            k = jstart
+            IF (k .LE. n-1) THEN
+
+*  if both points are before the first grid, go to the next point
+   30         CONTINUE
+                IF (x(k+1) .LE. xgl) THEN
+                   jstart = k - 1
+                   k = k+1
+                   IF (k .LE. n-1) GO TO 30
+                ENDIF
+
+
+*  if the last point is beyond the end of the grid, complete and go to the next
+*  grid
+   40         CONTINUE
+                 IF ((k .LE. n-1) .AND. (x(k) .LT. xgu)) THEN          
+
+                    jstart = k-1
+
+* compute x-coordinates of increment
+
+                    a1 = MAX(x(k),xgl)
+                    a2 = MIN(x(k+1),xgu)
+
+*  if points coincide, contribution is zero
+
+                    IF (x(k+1).EQ.x(k)) THEN
+                       darea = 0.e0
+                    ELSE
+                       slope = (y(k+1) - y(k))/(x(k+1) - x(k))
+                       b1 = y(k) + slope*(a1 - x(k))
+                       b2 = y(k) + slope*(a2 - x(k))
+                       darea = (a2 - a1)*(b2 + b1)/2.
+                    ENDIF
+
+
+*  find the area under the trapezoid from a1 to a2
+
+                    area = area + darea
+
+* go to next point
+              
+                    k = k+1
+                    GO TO 40
+
+                ENDIF
+
+            ENDIF
+
+*  calculate the average y after summing the areas in the interval
+            yg(i) = area/(xgu - xgl)
+
+   50 CONTINUE
+
+      RETURN
+      END
+
+!==============================================================================
+ 
+      subroutine rdxsco2(nw,wl,xsco2_195,xsco2_295,xsco2_370,yieldco2)
+
+*-----------------------------------------------------------------------------*
+*=  PURPOSE:                                                                 =*
+*=  Read and grid CO2 absorption cross-sections and photodissociation yield   =*
+*-----------------------------------------------------------------------------*
+*=  PARAMETERS:                                                              =*
+*=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
+*=           wavelength grid                                                 =*
+*=  WL     - REAL, vector of lower limits of wavelength intervals in      (I)=*
+*=           working wavelength grid                                         =*
+*=  XSCO2  - REAL, molecular absoprtion cross section (cm^2) of CO2 at    (O)=*
+*=           each specified wavelength                                       =*
+*-----------------------------------------------------------------------------*
+
+      use datafile_mod, only: datadir
+
+      implicit none
+
+!     input
+
+      integer :: nw               ! number of wavelength grid points
+      real, dimension(nw) :: wl   ! lower and central wavelength for each interval
+
+!     output
+
+      real, dimension(nw) :: xsco2_195, xsco2_295, xsco2_370 ! co2 cross-sections (cm2)
+      real, dimension(nw) :: yieldco2                        ! co2 photodissociation yield
+
+!     local
+
+      integer, parameter :: kdata = 36000
+      real, parameter :: deltax = 1.e-4
+      real, dimension(kdata) :: x1, y1, y2, y3, xion, ion
+      real, dimension(nw) :: yg
+      real :: xl, xu
+      integer :: ierr, i, l, n, n1, n2, n3, n4
+      CHARACTER*100 fil
+
+      integer :: kin, kout ! input/ouput logical units
+
+      kin  = 10
+      kout = 30
+
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+c     CO2 absorption cross-sections
+c
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+c     195K: Chan et al. (1993) + Starck et al. (2006) 
+c           + Yoshino et al. (1996) + Parkinson et al. (2003) + extrapolation
+c           11 lignes d'en-tete et n1 = 34537 
+c           fichier: co2_euv_uv_195k.txt 
+c
+c     295K: Chan et al. (1993) + Starck et al. (2006) 
+c           + Yoshino et al. (1996) + Parkinson et al. (2003) + extrapolation
+c           11 lignes d'en-tete et n2 = 35360 
+c           fichier: co2_euv_uv_295k.txt 
+c
+c     370K: Chan et al. (1993) + Starck et al. (2006) 
+c           + Lewis and Carver (1983) + extrapolation
+c           11 lignes d'en-tete et n3 = 3884 
+c           fichier: co2_euv_uv_370k.txt 
+c
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+      n1 = 34537
+      n2 = 35360
+      n3 = 3884
+c
+c     195K:
+c
+      fil = trim(datadir)//'/cross_sections/co2_euv_uv_195k.txt'
+      print*, 'section efficace CO2 195K: ', fil
+
+      OPEN(UNIT=kin,FILE=fil,STATUS='old')
+      DO i = 1,11 
+         read(kin,*)
+      END DO
+
+      DO i = 1, n1
+         READ(kin,*) x1(i), y1(i)
+      END DO
+      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,yg,n1,x1,y1,ierr)
+      IF (ierr .NE. 0) THEN
+         WRITE(*,*) ierr, fil
+         STOP
+      ENDIF
+      
+      DO l = 1, nw-1
+         xsco2_195(l) = yg(l)
+      END DO
+c
+c     295K:
+c
+      fil = trim(datadir)//'/cross_sections/co2_euv_uv_295k.txt'
+      print*, 'section efficace CO2 295K: ', fil
+c
+      OPEN(UNIT=kin,FILE=fil,STATUS='old')
+      DO i = 1,11
+         read(kin,*)
+      END DO
+c
+      DO i = 1, n2
+         READ(kin,*) x1(i), y1(i)
+      END DO
+      CLOSE (kin)
+
+      CALL addpnt(x1,y1,kdata,n2,x1(1)*(1.-deltax),0.)
+      CALL addpnt(x1,y1,kdata,n2,          0.,0.)
+      CALL addpnt(x1,y1,kdata,n2,x1(n2)*(1.+deltax),0.)
+      CALL addpnt(x1,y1,kdata,n2,      1.e+38,0.)
+      CALL inter2(nw,wl,yg,n2,x1,y1,ierr)
+      IF (ierr .NE. 0) THEN
+         WRITE(*,*) ierr, fil
+         STOP
+      ENDIF
+c
+      DO l = 1, nw-1
+         xsco2_295(l) = yg(l)
+      END DO
+c
+c     370K:
+c
+      fil = trim(datadir)//'/cross_sections/co2_euv_uv_370k.txt'
+      print*, 'section efficace CO2 370K: ', fil
+c
+      OPEN(UNIT=kin,FILE=fil,STATUS='old')
+      DO i = 1,11
+         read(kin,*)
+      END DO
+c
+      DO i = 1, n3
+         READ(kin,*) x1(i), y1(i)
+      END DO
+      CLOSE (kin)
+c
+      CALL addpnt(x1,y1,kdata,n3,x1(1)*(1.-deltax),0.)
+      CALL addpnt(x1,y1,kdata,n3,          0.,0.)
+      CALL addpnt(x1,y1,kdata,n3,x1(n3)*(1.+deltax),0.)
+      CALL addpnt(x1,y1,kdata,n3,      1.e+38,0.)
+      CALL inter2(nw,wl,yg,n3,x1,y1,ierr)
+      IF (ierr .NE. 0) THEN
+         WRITE(*,*) ierr, fil
+         STOP
+      ENDIF
+
+      DO l = 1, nw-1
+         xsco2_370(l) = yg(l)
+      END DO
+
+!     photodissociation yield:
+
+      fil = trim(datadir)//
+     $'/cross_sections/efdis_co2-o2_schunkandnagy2000.txt'
+      print*, 'photodissociation yield CO2: ', fil
+
+      OPEN(UNIT=kin,FILE=fil,STATUS='old')
+
+      do i = 1,3
+         read(kin,*)
+      end do
+
+      n4 = 17
+      do i = 1, n4
+         read(kin,*) xl, xu, ion(i)
+         xion(i) = (xl + xu)/2.
+         ion(i) = max(ion(i), 0.)
+      end do
+      close(kin)
+
+      CALL addpnt(xion,ion,kdata,n4,xion(1)*(1.-deltax),0.)
+      CALL addpnt(xion,ion,kdata,n4,          0.,0.)
+      CALL addpnt(xion,ion,kdata,n4,xion(n4)*(1.+deltax),1.)
+      CALL addpnt(xion,ion,kdata,n4,      1.e+38,1.)
+      CALL inter2(nw,wl,yieldco2,n4,xion,ion,ierr)
+      IF (ierr .NE. 0) THEN
+         WRITE(*,*) ierr, fil
+         STOP
+      ENDIF
+
+!     DO l = 1, nw-1
+!        write(kout,*) wl(l), xsco2_195(l),
+!    $                        xsco2_295(l),
+!    $                        xsco2_370(l),
+!    $                        yieldco2(l)
+!     END DO
+
+      end subroutine rdxsco2
+
+!==============================================================================
+ 
+      subroutine rdxso2(nw,wl,xso2_150,xso2_200,xso2_250,xso2_300,
+     $                  yieldo2)
+
+*-----------------------------------------------------------------------------*
+*=  PURPOSE:                                                                 =*
+*=  Read and grid O2 cross-sections and photodissociation yield              =*
+*-----------------------------------------------------------------------------*
+*=  PARAMETERS:                                                              =*
+*=  NW      - INTEGER, number of specified intervals + 1 in working       (I)=*
+*=            wavelength grid                                                =*
+*=  WL      - REAL, vector of lower limits of wavelength intervals in     (I)=*
+*=            working wavelength grid                                        =*
+*=  XSO2    - REAL, molecular absorption cross section                       =*
+*-----------------------------------------------------------------------------*
+
+      use datafile_mod, only: datadir
+
+      implicit none
+
+!     input
+
+      integer :: nw               ! number of wavelength grid points
+      real, dimension(nw) :: wl   ! lower and central wavelength for each interval
+
+!     output
+
+      real, dimension(nw) :: xso2_150, xso2_200, xso2_250, xso2_300 ! o2 cross-sections (cm2)
+      real, dimension(nw) :: yieldo2                                ! o2 photodissociation yield
+
+!     local
+
+      integer, parameter :: kdata = 18000
+      real, parameter :: deltax = 1.e-4
+      real, dimension(kdata) :: x1, y1, x2, y2, x3, y3, x4, y4
+      real, dimension(kdata) :: xion, ion
+      real    :: factor, xl, xu, dummy
+      integer :: i, ierr, n, n1, n2, n3, n4, nhead
+      integer :: kin, kout ! input/output logical units
+      character*100 fil
+
+      kin  = 10
+      kout = 30
+
+!     read o2 cross section data
+
+      nhead = 22
+      n     = 17434 
+
+      fil = trim(datadir)//'/cross_sections/o2_composite_2018_150K.txt'
+      print*, 'section efficace O2 150K: ', fil
+      OPEN(UNIT=kin,FILE=fil,STATUS='old')
+
+      DO i = 1,nhead 
+         read(kin,*)
+      END DO
+
+      n1 = n
+      DO i = 1, n1
+         READ(kin,*) x1(i), y1(i)
+      END DO
+      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,xso2_150,n1,x1,y1,ierr)
+      IF (ierr .NE. 0) THEN
+         WRITE(*,*) ierr, fil
+         STOP
+      ENDIF
+
+      fil = trim(datadir)//'/cross_sections/o2_composite_2018_200K.txt'
+      print*, 'section efficace O2 200K: ', fil
+      OPEN(UNIT=kin,FILE=fil,STATUS='old')
+
+      DO i = 1,nhead 
+         read(kin,*)
+      END DO
+
+      n2 = n
+      DO i = 1, n2
+         READ(kin,*) x2(i), y2(i)
+      END DO
+      CLOSE (kin)
+
+      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,xso2_200,n2,x2,y2,ierr)
+      IF (ierr .NE. 0) THEN
+         WRITE(*,*) ierr, fil
+         STOP
+      ENDIF
+
+      fil = trim(datadir)//'/cross_sections/o2_composite_2018_250K.txt'
+      print*, 'section efficace O2 250K: ', fil
+      OPEN(UNIT=kin,FILE=fil,STATUS='old')
+
+      DO i = 1,nhead 
+         read(kin,*)
+      END DO
+
+      n3 = n
+      DO i = 1, n3
+         READ(kin,*) x3(i), y3(i)
+      END DO
+      CLOSE (kin)
+
+      CALL addpnt(x3,y3,kdata,n3,x3(1)*(1.-deltax),0.)
+      CALL addpnt(x3,y3,kdata,n3,               0.,0.)
+      CALL addpnt(x3,y3,kdata,n3,x3(n3)*(1.+deltax),0.)
+      CALL addpnt(x3,y3,kdata,n3,           1.e+38,0.)
+      CALL inter2(nw,wl,xso2_250,n3,x3,y3,ierr)
+      IF (ierr .NE. 0) THEN
+         WRITE(*,*) ierr, fil
+         STOP
+      ENDIF
+
+      fil = trim(datadir)//'/cross_sections/o2_composite_2018_300K.txt'
+      print*, 'section efficace O2 300K: ', fil
+      OPEN(UNIT=kin,FILE=fil,STATUS='old')
+
+      DO i = 1,nhead 
+         read(kin,*)
+      END DO
+
+      n4 = n
+      DO i = 1, n4
+         READ(kin,*) x4(i), y4(i)
+      END DO
+      CLOSE (kin)
+
+      CALL addpnt(x4,y4,kdata,n4,x4(1)*(1.-deltax),0.)
+      CALL addpnt(x4,y4,kdata,n4,               0.,0.)
+      CALL addpnt(x4,y4,kdata,n4,x4(n4)*(1.+deltax),0.)
+      CALL addpnt(x4,y4,kdata,n4,           1.e+38,0.)
+      CALL inter2(nw,wl,xso2_300,n4,x4,y4,ierr)
+      IF (ierr .NE. 0) THEN
+         WRITE(*,*) ierr, fil
+         STOP
+      ENDIF
+
+!     photodissociation yield
+
+      fil = trim(datadir)//
+     $     '/cross_sections/efdis_co2-o2_schunkandnagy2000.txt'
+      OPEN(UNIT=kin,FILE=fil,STATUS='old')
+
+      do i = 1,11
+         read(kin,*)
+      end do
+
+      n = 9
+      DO i = 1, n
+         READ(kin,*) xl, xu, dummy, ion(i)
+         xion(i) = (xl + xu)/2.
+         ion(i) = max(ion(i), 0.)
+      END DO
+      CLOSE (kin)
+
+      CALL addpnt(xion,ion,kdata,n,xion(1)*(1.-deltax),0.)
+      CALL addpnt(xion,ion,kdata,n,          0.,0.)
+      CALL addpnt(xion,ion,kdata,n,xion(n)*(1.+deltax),1.)
+      CALL addpnt(xion,ion,kdata,n,      1.e+38,1.)
+      CALL inter2(nw,wl,yieldo2,n,xion,ion,ierr)
+      IF (ierr .NE. 0) THEN
+         WRITE(*,*) ierr, fil
+         STOP
+      ENDIF
+
+      end subroutine rdxso2
+
+!==============================================================================
+ 
+      subroutine rdxso3(nw,wl,xso3_218,xso3_298)
+
+*-----------------------------------------------------------------------------*
+*=  PURPOSE:                                                                 =*
+*=  Read ozone molecular absorption cross section.  Re-grid data to match    =*
+*=  specified wavelength working grid.                                       =*
+*-----------------------------------------------------------------------------*
+*=  PARAMETERS:                                                              =*
+*=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
+*=           wavelength grid                                                 =*
+*=  WL     - REAL, vector of lower limits of wavelength intervals in      (I)=*
+*=           working wavelength grid                                         =*
+*=  XSO3_218 REAL, molecular absoprtion cross section (cm^2) of O3 at     (O)=*
+*=           each specified wavelength (JPL 2006)  218 K                     =*
+*=  XSO3_298 REAL, molecular absoprtion cross section (cm^2) of O3 at     (O)=*
+*=           each specified wavelength (JPL 2006)  298 K                     =*
+*-----------------------------------------------------------------------------*
+
+      use datafile_mod, only: datadir
+
+      implicit none
+
+!     input
+
+      integer :: nw               ! number of wavelength grid points
+      real, dimension(nw) :: wl   ! lower and central wavelength for each interval
+
+!     output
+
+      real, dimension(nw) :: xso3_218, xso3_298 ! o3 cross-sections (cm2)
+
+!     local
+
+      integer, parameter :: kdata = 200
+      real, parameter :: deltax = 1.e-4
+      real, dimension(kdata) :: x1, x2, y1, y2
+      real, dimension(nw) :: yg
+      real :: a1, a2
+
+      integer :: i, ierr, iw, n, n1, n2
+      integer :: kin, kout ! input/output logical units
+
+      character*100 fil
+
+!     JPL 2006 218 K
+
+      fil = trim(datadir)//
+     $     '/cross_sections/o3_cross-sections_jpl_2006_218K.txt' 
+      print*, 'section efficace O3 218K: ', fil
+
+      OPEN(UNIT=kin,FILE=fil,STATUS='old')
+      n1 = 167
+      DO i = 1, n1
+         READ(kin,*) a1, a2, y1(i)
+         x1(i) = (a1+a2)/2.
+      END DO
+      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,yg,n1,x1,y1,ierr)
+      IF (ierr .NE. 0) THEN
+         WRITE(*,*) ierr, fil
+         STOP
+      ENDIF
+c            
+      DO iw = 1, nw-1
+         xso3_218(iw) = yg(iw)
+      END DO
+
+!     JPL 2006 298 K
+
+      fil = trim(datadir)//
+     $      '/cross_sections/o3_cross-sections_jpl_2006_298K.txt' 
+      print*, 'section efficace O3 298K: ', fil
+
+      OPEN(UNIT=kin,FILE=fil,STATUS='old')
+      n2 = 167
+      DO i = 1, n2
+         READ(kin,*) a1, a2, y2(i)
+         x2(i) = (a1+a2)/2.
+      END DO
+      CLOSE (kin)
+
+      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,yg,n2,x2,y2,ierr)
+      IF (ierr .NE. 0) THEN
+         WRITE(*,*) ierr, fil
+         STOP
+      ENDIF
+c
+      DO iw = 1, nw-1
+         xso3_298(iw) = yg(iw)
+      END DO
+
+      end subroutine rdxso3
+
+!==============================================================================
+
+      subroutine rdxsh2o(nw, wl, yg)
+
+*-----------------------------------------------------------------------------*
+*=  PURPOSE:                                                                 =*
+*=  Read H2O molecular absorption cross section.  Re-grid data to match      =*
+*=  specified wavelength working grid.                                       =*
+*-----------------------------------------------------------------------------*
+*=  PARAMETERS:                                                              =*
+*=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
+*=           wavelength grid                                                 =*
+*=  WL     - REAL, vector of lower limits of wavelength intervals in      (I)=*
+*=           working wavelength grid                                         =*
+*=  YG     - REAL, molecular absoprtion cross section (cm^2) of H2O at    (O)=*
+*=           each specified wavelength                                       =*
+*-----------------------------------------------------------------------------*
+
+      use datafile_mod, only: datadir
+
+      IMPLICIT NONE
+
+!     input
+
+      integer :: nw               ! number of wavelength grid points
+      real, dimension(nw) :: wl   ! lower and central wavelength for each interval
+
+!     output
+
+      real, dimension(nw) :: yg   ! h2o cross-sections (cm2)
+
+!     local
+
+      integer, parameter :: kdata = 500
+      real, parameter :: deltax = 1.e-4
+      REAL x1(kdata)
+      REAL y1(kdata)
+      INTEGER ierr
+      INTEGER i, n
+      CHARACTER*100 fil
+      integer :: kin, kout ! input/output logical units
+
+      kin = 10
+
+      fil = trim(datadir)//'/cross_sections/h2o_composite_250K.txt'
+      print*, 'section efficace H2O: ', fil
+
+      OPEN(UNIT=kin,FILE=fil,STATUS='old') 
+
+      DO i = 1,26 
+         read(kin,*)
+      END DO
+
+      n = 420 
+      DO i = 1, n
+         READ(kin,*) x1(i), y1(i)
+      END DO
+      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,      1.e+38,0.)
+      CALL inter2(nw,wl,yg,n,x1,y1,ierr)
+      IF (ierr .NE. 0) THEN
+         WRITE(*,*) ierr, fil
+         STOP
+      ENDIF
+      
+      end subroutine rdxsh2o
+
+!==============================================================================
+
+      subroutine rdxsh2o2(nw, wl, xsh2o2)
+
+*-----------------------------------------------------------------------------*
+*=  PURPOSE:                                                                 =*
+*=  Read and grid H2O2 cross-sections
+*=         H2O2 + hv -> 2 OH                                                 =*
+*=  Cross section:  Schuergers and Welge, Z. Naturforsch. 23a (1968) 1508    =*
+*=                  from 125 to 185 nm, then JPL97 from 190 to 350 nm.       =*
+*-----------------------------------------------------------------------------*
+*=  PARAMETERS:                                                              =*
+*=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
+*=           wavelength grid                                                 =*
+*=  WL     - REAL, vector of lower limits of wavelength intervals in      (I)=*
+*=           working wavelength grid                                         =*
+*-----------------------------------------------------------------------------*
+
+      use datafile_mod, only: datadir
+
+      implicit none
+
+!     input
+
+      integer :: nw               ! number of wavelength grid points
+      real, dimension(nw) :: wl   ! lower wavelength for each interval
+
+!     output
+
+      real, dimension(nw) :: xsh2o2   ! h2o2 cross-sections (cm2)
+
+!     local
+
+      real, parameter :: deltax = 1.e-4
+      integer, parameter :: kdata = 100
+      real, dimension(kdata) :: x1, y1
+      real, dimension(nw)    :: yg
+      integer :: i, ierr, iw, n, idum
+      integer :: kin, kout ! input/output logical units
+      character*100 fil
+
+      kin = 10
+
+!     read cross-sections
+
+      fil = trim(datadir)//'/cross_sections/h2o2_composite.txt'
+      print*, 'section efficace H2O2: ', fil
+
+      OPEN(kin,FILE=fil,STATUS='OLD')
+      READ(kin,*) idum,n
+      DO i = 1, idum-2
+         READ(kin,*)
+      ENDDO
+      DO i = 1, n
+         READ(kin,*) x1(i), y1(i)
+      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,           1.e+38,0.)
+      CALL inter2(nw,wl,yg,n,x1,y1,ierr)
+      IF (ierr .NE. 0) THEN
+         WRITE(*,*) ierr, fil
+         STOP
+      ENDIF
+
+      DO iw = 1, nw - 1
+         xsh2o2(iw) = yg(iw)
+      END DO
+
+      end subroutine rdxsh2o2
+
+!==============================================================================
+
+      subroutine rdxsho2(nw, wl, yg)
+
+*-----------------------------------------------------------------------------*
+*=  PURPOSE:                                                                 =*
+*=  Read ho2 cross-sections                                                  =*
+*=  JPL 2006 recommendation                                                  =*
+*-----------------------------------------------------------------------------*
+*=  PARAMETERS:                                                              =*
+*=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
+*=           wavelength grid                                                 =*
+*-----------------------------------------------------------------------------*
+
+      use datafile_mod, only: datadir
+
+      IMPLICIT NONE
+
+!     input
+
+      integer :: nw               ! number of wavelength grid points
+      real, dimension(nw) :: wl   ! lower wavelength for each interval
+
+!     output
+
+      real, dimension(nw) :: yg   ! ho2 cross-sections (cm2)
+
+!     local
+
+      real, parameter :: deltax = 1.e-4
+      integer, parameter :: kdata = 100
+      real, dimension(kdata) :: x1, y1
+      integer :: i, n, ierr
+      character*100 fil
+      integer :: kin, kout ! input/output logical units
+
+      kin = 10
+
+**** cross sections from Sander et al. [2003]
+
+      fil = trim(datadir)//'/cross_sections/ho2_jpl2003.txt'
+      print*, 'section efficace HO2: ', fil
+
+      OPEN(kin,FILE=fil,STATUS='OLD')
+      READ(kin,*) n
+      DO i = 1, n
+        READ(kin,*) x1(i), y1(i)
+      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)
+  
+      IF (ierr .NE. 0) THEN
+        WRITE(*,*) ierr, fil
+        STOP
+      ENDIF
+ 
+      end subroutine rdxsho2
+
+!==============================================================================
+
+      subroutine rdxsh2(nw, wl, wc, yg, yieldh2)
+
+*-----------------------------------------------------------------------------*
+*=  PURPOSE:                                                                 =*
+*=  Read h2 cross-sections and photodissociation yield                       =*
+*-----------------------------------------------------------------------------*
+*=  PARAMETERS:                                                              =*
+*=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
+*=           wavelength grid                                                 =*
+*-----------------------------------------------------------------------------*
+
+      use datafile_mod, only: datadir
+
+      implicit none
+
+!     input
+
+      integer :: nw                   ! number of wavelength grid points
+      real, dimension(nw) :: wl, wc   ! lower and central wavelength for each interval
+
+!     output
+
+      real, dimension(nw) :: yg        ! h2 cross-sections (cm2)
+      real, dimension(nw) :: yieldh2   ! photodissociation yield
+
+!     local
+
+      integer, parameter :: kdata = 1000
+      real, parameter :: deltax = 1.e-4
+      real, dimension(kdata) :: x1, y1, x2, y2
+      real :: xl, xu
+      integer :: i, iw, n, ierr
+      integer :: kin, kout ! input/output logical units
+      character*100 fil
+
+      kin = 10
+
+!     h2 cross sections
+
+      fil = trim(datadir)//'/cross_sections/h2secef.txt'
+      print*, 'section efficace H2: ', fil
+
+      OPEN(kin,FILE=fil,STATUS='OLD')
+
+      n = 792 
+      read(kin,*) ! avoid first line with wavelength = 0.
+      DO i = 1, n
+        READ(kin,*) x1(i), y1(i)
+      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)
+  
+      IF (ierr .NE. 0) THEN
+        WRITE(*,*) ierr, fil
+        STOP
+      ENDIF
+ 
+!     photodissociation yield
+
+      fil = trim(datadir)//'/cross_sections/h2_ionef_schunknagy2000.txt'
+      OPEN(UNIT=kin,FILE=fil,STATUS='old')
+
+      n = 19
+      read(kin,*)
+      DO i = 1, n
+         READ(kin,*) xl, xu, y2(i)
+         x2(i) = (xl + xu)/2.
+         y2(i) = max(1. - y2(i),0.)
+      END DO
+      CLOSE (kin)
+
+      CALL addpnt(x2,y2,kdata,n,x2(1)*(1.-deltax),0.)
+      CALL addpnt(x2,y2,kdata,n,          0.,0.)
+      CALL addpnt(x2,y2,kdata,n,x2(n)*(1.+deltax),1.)
+      CALL addpnt(x2,y2,kdata,n,      1.e+38,1.)
+      CALL inter2(nw,wl,yieldh2,n,x2,y2,ierr)
+      IF (ierr .NE. 0) THEN
+        WRITE(*,*) ierr, fil
+        STOP
+      ENDIF
+
+      end subroutine rdxsh2
+
+!==============================================================================
+
+      subroutine rdxsno2(nw,wl,xsno2,xsno2_220,xsno2_294,
+     $                   yldno2_248, yldno2_298)
+
+*-----------------------------------------------------------------------------*
+*=  PURPOSE:                                                                 =*
+*=  read and grid cross section + quantum yield for NO2                      =*
+*=  photolysis                                                               =*
+*=  Jenouvrier et al., 1996  200-238 nm
+*=  Vandaele et al., 1998    238-666 nm 220K and 294K 
+*=  quantum yield from jpl 2006
+*-----------------------------------------------------------------------------*
+*=  PARAMETERS:                                                              =*
+*=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
+*=           wavelength grid                                                 =*
+*=  WL     - REAL, vector of lower limits of wavelength intervals in      (I)=*
+*=           working wavelength grid                                         =*
+*=  SQ     - REAL, cross section x quantum yield (cm^2) for each          (O)=*
+*=           photolysis reaction defined, at each defined wavelength and     =*
+*=           at each defined altitude level                                  =*
+*-----------------------------------------------------------------------------*
+
+      use datafile_mod, only: datadir
+
+      implicit none
+
+!     input
+
+      integer :: nw               ! number of wavelength grid points
+      real, dimension(nw) :: wl   ! lower and central wavelength for each interval
+      
+!     output
+
+      real, dimension(nw) :: xsno2, xsno2_220, xsno2_294 ! no2 cross-sections (cm2)
+      real, dimension(nw) :: yldno2_248, yldno2_298      ! quantum yields at 248-298 k
+
+!     local
+
+      integer, parameter :: kdata = 28000
+      real, parameter :: deltax = 1.e-4
+      real, dimension(kdata) :: x1, x2, x3, x4, x5, y1, y2, y3, y4, y5
+      real, dimension(nw) :: yg1, yg2, yg3, yg4, yg5
+      real :: dum, qy
+      integer :: i, iw, n, n1, n2, n3, n4, n5, ierr
+      character*100 fil
+      integer :: kin, kout ! input/output logical units
+
+      kin = 10
+
+!*************** NO2 photodissociation
+
+*  Jenouvrier 1996 + Vandaele 1998 (JPL 2006)  
+
+      fil = trim(datadir)//'/cross_sections/no2_xs_jenouvrier.txt'
+      print*, 'section efficace NO2: ', fil
+
+      OPEN(UNIT=kin,FILE=fil,status='old')
+      DO i = 1, 3
+         READ(kin,*)
+      ENDDO
+      n1 = 10001
+      DO i = 1, n1
+         READ(kin,*) x1(i), y1(i)
+      end do
+
+      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)
+
+      fil = trim(datadir)//'/cross_sections/no2_xs_vandaele_294K.txt'
+      print*, 'section efficace NO2: ', fil
+
+      OPEN(UNIT=kin,FILE=fil,status='old')
+      DO i = 1, 3
+         READ(kin,*)
+      ENDDO
+      n2 = 27993
+      DO i = 1, n2
+         READ(kin,*) x2(i), y2(i)
+      end do
+
+      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)
+
+      fil = trim(datadir)//'/cross_sections/no2_xs_vandaele_220K.txt'
+      print*, 'section efficace NO2: ', fil
+
+      OPEN(UNIT=kin,FILE=fil,status='old')
+      DO i = 1, 3
+         READ(kin,*)
+      ENDDO
+      n3 = 27993
+      DO i = 1, n3
+         READ(kin,*) x3(i), y3(i)
+      end do
+
+      CALL addpnt(x3,y3,kdata,n3,x3(1)*(1.-deltax), 0.)
+      CALL addpnt(x3,y3,kdata,n3,               0., 0.)
+      CALL addpnt(x3,y3,kdata,n3,x3(n3)*(1.+deltax),0.)
+      CALL addpnt(x3,y3,kdata,n3,           1.e+38, 0.)
+      CALL inter2(nw,wl,yg3,n3,x3,y3,ierr)
+
+      do iw = 1, nw - 1
+         xsno2(iw)     = yg1(iw)
+         xsno2_294(iw) = yg2(iw)
+         xsno2_220(iw) = yg3(iw)
+      end do
+
+!     photodissociation efficiency from jpl 2006
+
+      fil = trim(datadir)//'/cross_sections/no2_yield_jpl2006.txt'
+      print*, 'quantum yield NO2: ', fil
+
+      OPEN(UNIT=kin,FILE=fil,STATUS='old')
+      DO i = 1, 5
+         READ(kin,*) 
+      ENDDO
+      n = 25
+      n4 = n
+      n5 = n
+      DO i = 1, n
+         READ(kin,*) x4(i), y4(i), y5(i)
+         x5(i) = x4(i)
+      ENDDO
+      CLOSE(kin)
+
+      CALL addpnt(x4,y4,kdata,n4,x4(1)*(1.-deltax),y4(1))
+      CALL addpnt(x4,y4,kdata,n4,               0.,y4(1))
+      CALL addpnt(x4,y4,kdata,n4,x4(n4)*(1.+deltax),  0.)
+      CALL addpnt(x4,y4,kdata,n4,           1.e+38,   0.)
+      CALL inter2(nw,wl,yg4,n4,x4,y4,ierr)
+      IF (ierr .NE. 0) THEN
+         WRITE(*,*) ierr, fil
+         STOP
+      ENDIF
+
+      CALL addpnt(x5,y5,kdata,n5,x5(1)*(1.-deltax),y5(1))
+      CALL addpnt(x5,y5,kdata,n5,               0.,y5(1))
+      CALL addpnt(x5,y5,kdata,n5,x5(n5)*(1.+deltax),  0.)
+      CALL addpnt(x5,y5,kdata,n5,           1.e+38,   0.)
+      CALL inter2(nw,wl,yg5,n5,x5,y5,ierr)
+      IF (ierr .NE. 0) THEN
+         WRITE(*,*) ierr, fil
+         STOP
+      ENDIF
+
+      do iw = 1, nw - 1
+         yldno2_298(iw) = yg4(iw)
+         yldno2_248(iw) = yg5(iw)
+      end do
+      
+      end subroutine rdxsno2
+
+!==============================================================================
+
+      subroutine rdxsno(nw, wl, yg, yieldno)
+
+*-----------------------------------------------------------------------------*
+*=  PURPOSE:                                                                 =*
+*=  Read NO cross-sections  and photodissociation efficiency                 =*
+*=  Lida et al 1986 (provided by Francisco Gonzalez-Galindo)                 =*
+*-----------------------------------------------------------------------------*
+*=  PARAMETERS:                                                              =*
+*=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
+*=           wavelength grid                                                 =*
+*-----------------------------------------------------------------------------*
+
+      use datafile_mod, only: datadir
+
+      implicit none
+
+!     input
+
+      integer :: nw               ! number of wavelength grid points
+      real, dimension(nw) :: wl   ! lower wavelength for each interval
+
+!     output
+
+      real, dimension(nw) :: yg        ! no cross-sections (cm2)
+      real, dimension(nw) :: yieldno   ! no photodissociation efficiency
+
+!     local
+
+      integer, parameter :: kdata = 110
+      real, parameter :: deltax = 1.e-4
+      real, dimension(kdata) :: x1, y1, x2, y2
+      integer :: i, iw, n, ierr
+      character*100 fil
+      integer :: kin, kout ! input/output logical units
+
+      kin = 10
+
+!     no cross-sections
+
+      fil = trim(datadir)//'/cross_sections/no_xs_francisco.txt'
+      print*, 'section efficace NO: ', fil
+      OPEN(kin,FILE=fil,STATUS='OLD')
+
+      n = 99
+      DO i = 1, n
+        READ(kin,*) x1(i), y1(i)
+      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)
+      IF (ierr .NE. 0) THEN
+         WRITE(*,*) ierr, fil
+         STOP
+      ENDIF
+  
+!     photodissociation yield
+
+      fil = trim(datadir)//'/cross_sections/noefdis.txt'
+      OPEN(UNIT=kin,FILE=fil,STATUS='old')
+
+      n = 33
+      DO i = 1, n
+         READ(kin,*) x2(n-i+1), y2(n-i+1)
+      END DO
+      CLOSE (kin)
+
+      CALL addpnt(x2,y2,kdata,n,x2(1)*(1.-deltax),0.)
+      CALL addpnt(x2,y2,kdata,n,          0.,0.)
+      CALL addpnt(x2,y2,kdata,n,x2(n)*(1.+deltax),1.)
+      CALL addpnt(x2,y2,kdata,n,      1.e+38,1.)
+      CALL inter2(nw,wl,yieldno,n,x2,y2,ierr)
+      IF (ierr .NE. 0) THEN
+        WRITE(*,*) ierr, fil
+        STOP
+      ENDIF
+
+      end subroutine rdxsno
+
+!==============================================================================
+
+      subroutine rdxsn2(nw, wl, yg, yieldn2)
+
+*-----------------------------------------------------------------------------*
+*=  PURPOSE:                                                                 =*
+*=  Read n2 cross-sections and photodissociation yield                       =*
+*-----------------------------------------------------------------------------*
+*=  PARAMETERS:                                                              =*
+*=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
+*=           wavelength grid                                                 =*
+*-----------------------------------------------------------------------------*
+
+      use datafile_mod, only: datadir
+
+      implicit none
+
+!     input
+
+      integer :: nw               ! number of wavelength grid points
+      real, dimension(nw) :: wl   ! lower wavelength for each interval
+
+!     output
+
+      real, dimension(nw) :: yg        ! n2 cross-sections (cm2)
+      real, dimension(nw) :: yieldn2   ! n2 photodissociation yield
+
+!     local
+
+      integer, parameter :: kdata = 1100
+      real, parameter :: deltax = 1.e-4
+      real, dimension(kdata) :: x1, y1, x2, y2
+      real :: xl, xu
+      integer :: i, iw, n, ierr
+      integer :: kin, kout ! input/output logical units
+      character*100 fil
+
+      kin = 10
+
+!     n2 cross sections
+
+      fil = trim(datadir)//'/cross_sections/n2secef_01nm.txt'
+      print*, 'section efficace N2: ', fil
+      OPEN(kin,FILE=fil,STATUS='OLD')
+
+      n = 1020
+      DO i = 1, n
+        READ(kin,*) x1(i), y1(i)
+      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)
+  
+      IF (ierr .NE. 0) THEN
+        WRITE(*,*) ierr, fil
+        STOP
+      ENDIF
+ 
+!     photodissociation yield
+
+      fil = trim(datadir)//'/cross_sections/n2_ionef_schunknagy2000.txt'
+      OPEN(UNIT=kin,FILE=fil,STATUS='old')
+
+      n = 19
+      read(kin,*)
+      DO i = 1, n
+         READ(kin,*) xl, xu, y2(i)
+         x2(i) = (xl + xu)/2.
+         y2(i) = 1. - y2(i)
+      END DO
+      CLOSE (kin)
+
+      CALL addpnt(x2,y2,kdata,n,x2(1)*(1.-deltax),0.)
+      CALL addpnt(x2,y2,kdata,n,          0.,0.)
+      CALL addpnt(x2,y2,kdata,n,x2(n)*(1.+deltax),1.)
+      CALL addpnt(x2,y2,kdata,n,      1.e+38,1.)
+      CALL inter2(nw,wl,yieldn2,n,x2,y2,ierr)
+      IF (ierr .NE. 0) THEN
+        WRITE(*,*) ierr, fil
+        STOP
+      ENDIF
+
+      end subroutine rdxsn2
+
+!==============================================================================
+ 
+      subroutine setalb(nw,wl,albedo)
+ 
+*-----------------------------------------------------------------------------*
+*=  PURPOSE:                                                                 =*
+*=  Set the albedo of the surface.  The albedo is assumed to be Lambertian,  =*
+*=  i.e., the reflected light is isotropic, and idependt of the direction    =*
+*=  of incidence of light.  Albedo can be chosen to be wavelength dependent. =*
+*-----------------------------------------------------------------------------*
+*=  PARAMETERS:                                                              =*
+*=  NW      - INTEGER, number of specified intervals + 1 in working       (I)=*
+*=            wavelength grid                                                =*
+*=  WL     - REAL, vector of lower limits of wavelength intervals in      (I)=*
+*=           working wavelength grid                                         =*
+*=  ALBEDO  - REAL, surface albedo at each specified wavelength           (O)=*
+*-----------------------------------------------------------------------------*
+
+      implicit none
+
+! input: (wavelength working grid data)
+
+      INTEGER nw
+      REAL wl(nw)
+
+! output:
+
+      REAL albedo(nw)
+
+! local:
+
+      INTEGER iw
+      REAL alb
+
+!     0.015: mean value from clancy et al., icarus, 49-63, 1999.
+
+      alb = 0.015
+
+      do iw = 1, nw - 1
+         albedo(iw) = alb
+      end do
+
+      return
+      end
+ 
+!==============================================================================
+
+      subroutine setair(nlev, nw, wl, wc, press, temp, zmmean, 
+     $                  colinc, dtrl)
+
+*-----------------------------------------------------------------------------*
+*=  PURPOSE:                                                                 =*
+*=  computes air column increments and rayleigh optical depth                =*
+*-----------------------------------------------------------------------------*
+
+      implicit none
+
+!     input:
+
+      integer :: nlev, nw
+
+      real, dimension(nw)   :: wl, wc              ! lower and central wavelength grid (nm)
+      real, dimension(nlev) :: press, temp, zmmean ! pressure (hpa), temperature (k), molecular mass (g.mol-1)
+
+!     output:
+
+      real, dimension(nlev) :: colinc    ! air column increments (molecule.cm-2)
+      real, dimension(nlev,nw) :: dtrl   ! rayleigh optical depth
+
+!     local:
+
+      real, parameter :: avo = 6.022e23
+      real, parameter :: g = 3.72
+      real :: dp
+      real :: alpha, rho, df, pi
+      real, dimension(nw) :: srayl
+      integer :: ilev, iw
+
+c     compute column increments 
+
+      do ilev = 1, nlev-1
+         dp = (press(ilev) - press(ilev+1))*100.
+         colinc(ilev) = avo*0.1*dp/(zmmean(ilev)*g)
+      end do
+      colinc(nlev) = avo*0.1*press(nlev)*100./(zmmean(nlev)*g)
+
+!     do ilev = nlev,1,-1
+!        print*, ilev, press(ilev), temp(ilev), zmmean(ilev), 
+!    $           colinc(ilev)
+!     end do
+
+      do iw = 1, nw - 1
+
+c        calcul de section efficace Rayleigh: voir Atreya and Gu, JGR, 
+c        13133-13145, 1994.
+c
+c        rho   : normal depolarization ratio       = 0.0774        for CO2
+c        df    : depolarization factor
+c        alpha : average dielectric polarizability = 2.911e-24 cm3 for CO2
+
+         pi    = acos(-1.0)
+         rho   = 0.0774
+         alpha = 2.911e-24*1.e21  ! nm3
+
+         df = (6. + 3.*rho)/(6. - 7.*rho)
+
+         srayl(iw) = 8.*pi/3.*((2.*pi/wc(iw))**4)*(alpha**2)*df   ! nm2
+
+         do ilev = 1, nlev
+            dtrl(ilev,iw) = colinc(ilev)*srayl(iw)*1.e-14     ! cm2
+         end do
+      end do
+
+      return
+      end
+
+!==============================================================================
+
+      subroutine setco2(nd, nlayer, nw, wc, tlay, xsco2_195, xsco2_295, 
+     $                  xsco2_370, yieldco2, j_co2_o, j_co2_o1d,
+     $                  colinc, rm, dt, sj)
+
+*-----------------------------------------------------------------------------*
+*=  PURPOSE:                                                                 =*
+*=  Set up the CO2 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_co2_o, j_co2_o1d                              ! photolysis indexes
+      real, dimension(nw)     :: wc                              ! central wavelength for each interval
+      real, dimension(nw)     :: xsco2_195, xsco2_295, xsco2_370 ! co2 cross-sections (cm2)
+      real, dimension(nw)     :: yieldco2                        ! co2 photodissociation yield
+      real, dimension(nlayer) :: tlay                            ! temperature (k)
+      real, dimension(nlayer) :: rm                              ! co2 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 :: extrapol
+      integer :: i, l
+      real    :: temp, sco2
+
+!     extrapol  = 0         no extrapolation below 195 k
+!     extrapol  = 1         extrapolation below 195 k
+
+      extrapol  = 0
+
+      do i = 1, nlayer
+         if (extrapol == 1) then
+            temp = tlay(i)             
+         else
+            temp = max(tlay(i), 195.)
+         end if
+         temp = min(temp, 370.)
+         do l = 1, nw-1
+            if (temp <= 295.) then
+               if (xsco2_195(l) /= 0. .and. xsco2_295(l) /= 0.) then
+                  sco2 =  alog(xsco2_195(l))
+     $                 + (alog(xsco2_295(l)) - alog(xsco2_195(l)))
+     $                   /(295. - 195.)*(temp - 195.)
+                  sco2 = exp(sco2)
+               else
+                  sco2 = 0.
+               end if
+            else
+               if (xsco2_295(l) /= 0. .and. xsco2_370(l) /= 0.) then
+                  sco2 =  alog(xsco2_295(l))
+     $                 + (alog(xsco2_370(l)) - alog(xsco2_295(l)))
+     $                   /(370. - 295.)*(temp - 295.)
+                  sco2 = exp(sco2)
+               else
+                  sco2 = 0.
+               end if
+            end if
+
+!           optical depth
+
+            dt(i,l) = colinc(i)*rm(i)*sco2
+   
+!           production of o(1d) for wavelengths shorter than 167 nm
+
+            if (wc(l) >= 167.) then
+               sj(i,l,j_co2_o) = sco2*yieldco2(l)
+               sj(i,l,j_co2_o1d) = 0.
+            else
+               sj(i,l,j_co2_o) = 0.
+               sj(i,l,j_co2_o1d) = sco2*yieldco2(l)
+            end if
+         end do
+      end do
+
+!     do l = 1,nw-1
+!        write(35,*) wc(l), sj(1,l,j_co2_o1d),
+!    $                      sj(1,l,j_co2_o), tlay(1)
+!     end do
+
+      end subroutine setco2
+
+!==============================================================================
+
+      subroutine seto2(nd, nlayer, nw, wc, mopt, tlay, xso2_150, 
+     $                 xso2_200, xso2_250, xso2_300, yieldo2, j_o2_o, 
+     $                 j_o2_o1d, colinc, rm, dt, sj)
+
+*-----------------------------------------------------------------------------*
+*=  PURPOSE:                                                                 =*
+*=  Set up the O2 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 :: mopt                                          ! high-res/low-res switch
+      integer :: j_o2_o, j_o2_o1d                              ! photolysis indexes
+      real, dimension(nw)     :: wc                            ! central wavelength for each interval
+      real, dimension(nw)     :: xso2_150, xso2_200, xso2_250, ! o2 cross-sections (cm2)
+     $                           xso2_300
+      real, dimension(nw)     :: yieldo2                       ! o2 photodissociation yield
+      real, dimension(nlayer) :: tlay                          ! temperature (k)
+      real, dimension(nlayer) :: rm                            ! o2 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
+      real    :: xso2, factor
+
+!     correction by factor if low-resolution in schumann-runge bands
+
+      if (mopt == 1) then
+         factor = 1.
+      else if (mopt == 2) then
+         factor = 0.8
+      end if
+
+!     calculate temperature dependance
+
+      do ilev = 1,nlayer
+         temp = max(tlay(ilev),150.)
+         temp = min(temp, 300.)
+         do iw = 1, nw-1
+            if (tlay(ilev) > 250.) then
+               xso2 = xso2_250(iw) + (xso2_300(iw) - xso2_250(iw))
+     $                              /(300. - 250.)*(temp - 250.) 
+            else if (tlay(ilev) > 200.) then
+               xso2 = xso2_200(iw) + (xso2_250(iw) - xso2_200(iw))
+     $                              /(250. - 200.)*(temp - 200.) 
+            else
+               xso2 = xso2_150(iw) + (xso2_200(iw) - xso2_150(iw))
+     $                              /(200. - 150.)*(temp - 150.) 
+            end if
+
+            if (wc(iw) > 180. .and. wc(iw) < 200.) then
+               xso2 = xso2*factor
+            end if
+
+!     optical depth
+
+            dt(ilev,iw) = colinc(ilev)*rm(ilev)*xso2
+
+!     production of o(1d) for wavelengths shorter than 175 nm
+
+            if (wc(iw) >= 175.) then
+               sj(ilev,iw,j_o2_o)   = xso2*yieldo2(iw)
+               sj(ilev,iw,j_o2_o1d) = 0.
+            else
+               sj(ilev,iw,j_o2_o)   = 0.
+               sj(ilev,iw,j_o2_o1d) = xso2*yieldo2(iw)
+            end if
+
+         end do
+      end do
+
+      end subroutine seto2
+
+!==============================================================================
+
+      subroutine seto3(nd, nlayer, nw, wc, tlay, xso3_218, xso3_298, 
+     $                 j_o3_o, j_o3_o1d,
+     $                 colinc, rm, dt, sj)
+
+*-----------------------------------------------------------------------------*
+*=  PURPOSE:                                                                 =*
+*=  Set up the O3 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_o3_o, j_o3_o1d                              ! photolysis indexes
+      real, dimension(nw)     :: wc                            ! central wavelength for each interval
+      real, dimension(nw)     :: xso3_218, xso3_298            ! o3 cross-sections (cm2)
+      real, dimension(nlayer) :: tlay                          ! temperature (k)
+      real, dimension(nlayer) :: rm                            ! o3 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
+      real, dimension(nw) :: xso3(nw)
+      real, dimension(nw) :: qy1d                 ! quantum yield for o(1d) production
+      real :: q1, q2, a1, a2, a3
+
+      do ilev = 1, nlayer
+         temp = max(tlay(ilev), 218.)
+         temp = min(temp,298.)
+         do iw = 1, nw-1
+            xso3(iw) = xso3_218(iw) + (xso3_298(iw) - xso3_218(iw))
+     $                                /(298. - 218.) *(temp - 218.) 
+
+!     optical depth
+
+            dt(ilev,iw) = colinc(ilev)*rm(ilev)*xso3(iw)
+
+         end do
+
+!     calculate quantum yield for o(1d) production (jpl 2006)
+
+         temp = max(tlay(ilev),200.)
+         temp = min(temp,320.)
+         do iw = 1, nw-1
+            if (wc(iw) <= 306.) then
+               qy1d(iw) = 0.90
+            else if (wc(iw) > 306. .and. wc(iw) < 328.) then
+               q1 = 1.
+               q2 = exp(-825.518/(0.695*temp))
+               a1 = (304.225 - wc(iw))/5.576
+               a2 = (314.957 - wc(iw))/6.601
+               a3 = (310.737 - wc(iw))/2.187
+               qy1d(iw) = (q1/(q1 + q2))*0.8036*exp(-(a1*a1*a1*a1))
+     $                  + (q2/(q1 + q2))*8.9061*(temp/300.)**2.
+     $                   *exp(-(a2*a2))
+     $                  + 0.1192*(temp/300.)**1.5*exp(-(a3*a3))
+     $                  + 0.0765
+            else if (wc(iw) >= 328. .and. wc(iw) <= 340.) then
+               qy1d(iw) = 0.08
+            else
+               qy1d(iw) = 0.
+            endif
+         end do
+         do iw = 1, nw-1
+            sj(ilev,iw,j_o3_o)   = xso3(iw)*(1. - qy1d(iw)) 
+            sj(ilev,iw,j_o3_o1d) = xso3(iw)*qy1d(iw)
+         end do
+      end do
+
+      end subroutine seto3
+
+!==============================================================================
+
+      subroutine seth2o2(nd, nlayer, nw, wc, tlay, xsh2o2, j_h2o2,
+     $                   colinc, rm, dt, sj)
+
+*-----------------------------------------------------------------------------*
+*=  PURPOSE:                                                                 =*
+*=  Set up the h2o2 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_h2o2                                        ! photolysis index
+      real, dimension(nw)     :: wc                            ! central wavelength for each interval
+      real, dimension(nw)     :: xsh2o2                        ! h2o2 cross-sections (cm2)
+      real, dimension(nlayer) :: tlay                          ! temperature (k)
+      real, dimension(nlayer) :: rm                            ! h2o2 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    :: a0, a1, a2, a3, a4, a5, a6, a7
+      real    :: b0, b1, b2, b3, b4
+      real    :: lambda, suma, sumb, chi, temp, xs
+
+      A0 = 6.4761E+04            
+      A1 = -9.2170972E+02        
+      A2 = 4.535649              
+      A3 = -4.4589016E-03        
+      A4 = -4.035101E-05         
+      A5 = 1.6878206E-07
+      A6 = -2.652014E-10
+      A7 = 1.5534675E-13
+
+      B0 = 6.8123E+03
+      B1 = -5.1351E+01
+      B2 = 1.1522E-01
+      B3 = -3.0493E-05
+      B4 = -1.0924E-07
+
+!     temperature dependance: jpl 2006
+
+      do ilev = 1,nlayer
+         temp = min(max(tlay(ilev),200.),400.)            
+         chi = 1./(1. + exp(-1265./temp))
+         do iw = 1, nw-1
+            if ((wc(iw) >= 260.) .and. (wc(iw) < 350.)) then
+               lambda = wc(iw)
+               sumA = ((((((A7*lambda + A6)*lambda + A5)*lambda + 
+     $                      A4)*lambda +A3)*lambda + A2)*lambda + 
+     $                      A1)*lambda + A0
+               sumB = (((B4*lambda + B3)*lambda + B2)*lambda + 
+     $                   B1)*lambda + B0
+               xs = (chi*sumA + (1. - chi)*sumB)*1.e-21
+               sj(ilev,iw,j_h2o2) = xs
+            else
+               sj(ilev,iw,j_h2o2) = xsh2o2(iw)
+            end if
+
+!     optical depth
+
+            dt(ilev,iw) = colinc(ilev)*rm(ilev)*sj(ilev,iw,j_h2o2)
+         end do
+      end do
+
+      end subroutine seth2o2
+
+!==============================================================================
+
+      subroutine setno2(nd, nlayer, nw, wc, tlay, xsno2, xsno2_220,
+     $                  xsno2_294, yldno2_248, yldno2_298, j_no2,
+     $                  colinc, rm, dt, sj)
+
+*-----------------------------------------------------------------------------*
+*=  PURPOSE:                                                                 =*
+*=  Set up the no2 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_no2                                         ! photolysis index
+      real, dimension(nw)     :: wc                            ! central wavelength for each interval
+      real, dimension(nw) :: xsno2, xsno2_220, xsno2_294       ! no2 absorption cross-section at 220-294 k (cm2)
+      real, dimension(nw) :: yldno2_248, yldno2_298            ! no2 quantum yield at 248-298 k
+      real, dimension(nlayer) :: tlay                          ! temperature (k)
+      real, dimension(nlayer) :: rm                            ! no2 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, qy
+
+!     temperature dependance: jpl 2006
+
+      do ilev = 1,nlayer
+         temp = max(220.,min(tlay(ilev),294.))
+         do iw = 1, nw - 1
+            if (wc(iw) < 238.) then
+               sj(ilev,iw,j_no2) = xsno2(iw)
+            else
+               sj(ilev,iw,j_no2) = xsno2_220(iw) 
+     $                           + (xsno2_294(iw) - xsno2_220(iw))
+     $                            /(294. - 220.)*(temp - 220.)
+            end if
+
+!     optical depth
+
+            dt(ilev,iw) = colinc(ilev)*rm(ilev)*sj(ilev,iw,j_no2)
+         end do
+      end do
+
+!     quantum yield: jpl 2006
+
+      do ilev = 1,nlayer
+         temp = max(248.,min(tlay(ilev),298.))
+         do iw = 1, nw - 1
+            qy = yldno2_248(iw) + (yldno2_298(iw) - yldno2_248(iw))
+     $                           /(298. - 248.)*(temp - 248.)
+            sj(ilev,iw,j_no2) = sj(ilev,iw,j_no2)*qy 
+         end do
+      end do
+
+!     do iw = 1, nw-1
+!        write(35,*) wc(iw), sj(1,iw,j_no2), tlay(1)
+!     end do
+!     stop
+
+      end subroutine setno2
+
+!==============================================================================
+
+      subroutine setaer(nlayer,alt,tau,nw,dtaer,omaer,gaer)
+
+*-----------------------------------------------------------------------------*
+*=  PURPOSE:                                                                 =*
+*=  Set aerosol properties for each specified altitude layer.  Properties    =*
+*=  may be wavelength dependent.                                             =*
+*-----------------------------------------------------------------------------*
+
+      implicit none
+
+!     input
+
+      integer :: nlayer                              ! number of vertical layers
+      integer :: nw                                  ! number of wavelength grid points
+      real, dimension(nlayer) :: alt                 ! altitude (km)
+      real :: tau                                    ! integrated aerosol optical depth at the surface
+
+!     output 
+
+      real, dimension(nlayer,nw) :: dtaer            ! aerosol optical depth 
+      real, dimension(nlayer,nw) :: omaer            ! aerosol single scattering albedo 
+      real, dimension(nlayer,nw) :: gaer             ! aerosol asymmetry parameter
+
+!     local
+
+      integer :: ilay, iw
+      real, dimension(nlayer) :: aer                                 ! dust extinction
+      real :: omega, g, scaleh, gamma
+      real :: dz, tautot, q0
+
+      omega  = 0.622 ! single scattering albedo : wolff et al.(2010) at 258 nm
+      g      = 0.88  ! asymmetry factor : mateshvili et al. (2007) at 210 nm
+      scaleh = 10.   ! scale height (km)
+      gamma  = 0.03  ! conrath parameter
+
+      dtaer(:,:) = 0.
+      omaer(:,:) = 0.
+      gaer(:,:)  = 0.
+
+!     optical depth profile:
+
+      tautot = 0.
+      do ilay = 1, nlayer-1
+         dz = alt(ilay+1) - alt(ilay)
+         tautot = tautot + exp(gamma*(1. - exp(alt(ilay)/scaleh)))*dz
+      end do
+      
+!     print*, 'tau = ', tau
+!     print*, 'q0  = ', q0
+      q0 = tau/tautot
+      do ilay = 1, nlayer-1
+         dz = alt(ilay+1) - alt(ilay)
+         dtaer(ilay,:) = q0*exp(gamma*(1. - exp(alt(ilay)/scaleh)))*dz
+         omaer(ilay,:) = omega
+         gaer(ilay,:)  = g
+      end do
+
+!     do ilay = nlayer,1,-1
+!        print*, alt(ilay), q0*exp(gamma*(1. - exp(alt(ilay)/scaleh))), 
+!    $           dtaer(ilay,1)
+!     end do
+
+      end subroutine setaer
+
+!==============================================================================
+
+      subroutine setcld(nlayer,nw,dtcld,omcld,gcld)
+
+*-----------------------------------------------------------------------------*
+*=  PURPOSE:                                                                 =*
+*=  Set cloud properties for each specified altitude layer.  Properties      =*
+*=  may be wavelength dependent.                                             =*
+*-----------------------------------------------------------------------------*
+
+      implicit none
+
+!     input
+
+      integer :: nlayer                              ! number of vertical layers
+      integer :: nw                                  ! number of wavelength grid points
+
+!     output 
+
+      real, dimension(nlayer,nw) :: dtcld            ! cloud optical depth 
+      real, dimension(nlayer,nw) :: omcld            ! cloud single scattering albedo 
+      real, dimension(nlayer,nw) :: gcld             ! cloud asymmetry parameter
+
+!     local
+
+      integer :: ilay, iw
+
+c     dtcld : optical depth
+c     omcld : single scattering albedo
+c     gcld  : asymmetry factor
+
+      do ilay = 1, nlayer - 1
+         do iw = 1, nw - 1
+            dtcld(ilay,iw) = 0.       ! no clouds for the moment
+            omcld(ilay,iw) = 0.99
+            gcld(ilay,iw)  = 0.85
+         end do
+      end do
+
+      end
+
+!==============================================================================
+
+      subroutine sphers(nlev, z, zen, dsdh, nid)
+
+*-----------------------------------------------------------------------------*
+*=  PURPOSE:                                                                 =*
+*=  Calculate slant path over vertical depth ds/dh in spherical geometry.    =*
+*=  Calculation is based on:  A.Dahlback, and K.Stamnes, A new spheric model =*
+*=  for computing the radiation field available for photolysis and heating   =*
+*=  at twilight, Planet.Space Sci., v39, n5, pp. 671-683, 1991 (Appendix B)  =*
+*-----------------------------------------------------------------------------*
+*=  PARAMETERS:                                                              =*
+*=  NZ      - INTEGER, number of specified altitude levels in the working (I)=*
+*=            grid                                                           =*
+*=  Z       - REAL, specified altitude working grid (km)                  (I)=*
+*=  ZEN     - REAL, solar zenith angle (degrees)                          (I)=*
+*=  DSDH    - REAL, slant path of direct beam through each layer crossed  (O)=*
+*=            when travelling from the top of the atmosphere to layer i;     =*
+*=            DSDH(i,j), i = 0..NZ-1, j = 1..NZ-1                            =*
+*=  NID     - INTEGER, number of layers crossed by the direct beam when   (O)=*
+*=            travelling from the top of the atmosphere to layer i;          =*
+*=            NID(i), i = 0..NZ-1                                            =*
+*-----------------------------------------------------------------------------*
+c
+      implicit none
+
+* input
+      INTEGER nlev
+      REAL zen, z(nlev)
+
+* output
+      INTEGER nid(0:nlev)
+      REAL dsdh(0:nlev,nlev)
+
+* more program constants
+      REAL re, ze(nlev)
+      REAL  dr
+      real radius
+      parameter (radius = 3393.)
+
+* local 
+
+      REAL pi, zenrad, rpsinz, rj, rjp1, dsj, dhj, ga, gb, sm
+      INTEGER i, j, k
+      INTEGER id
+
+      INTEGER nlay
+      REAL zd(0:nlev-1)
+
+*-----------------------------------------------------------------------------
+
+      pi = acos(-1.0)
+      dr = pi/180.
+      zenrad = zen*dr
+
+* number of layers:
+      nlay = nlev - 1
+
+* include the elevation above sea level to the radius of Mars:
+      re = radius + z(1)
+* correspondingly z changed to the elevation above Mars surface:
+      DO k = 1, nlev
+         ze(k) = z(k) - z(1)
+      END DO
+
+* inverse coordinate of z
+      zd(0) = ze(nlev)
+      DO k = 1, nlay
+        zd(k) = ze(nlev - k)
+      END DO
+
+* initialize dsdh(i,j), nid(i)
+      DO i = 0, nlev
+       nid(i) = 0
+       DO j = 1, nlev
+        dsdh(i,j) = 0.
+       END DO
+      END DO
+
+* calculate ds/dh of every layer
+      DO 100 i = 0, nlay
+
+        rpsinz = (re + zd(i)) * SIN(zenrad)
+ 
+        IF ( (zen .GT. 90.0) .AND. (rpsinz .LT. re) ) THEN
+           nid(i) = -1
+        ELSE
+
+*
+* Find index of layer in which the screening height lies
+*
+           id = i 
+           IF( zen .GT. 90.0 ) THEN
+              DO 10 j = 1, nlay
+                 IF( (rpsinz .LT. ( zd(j-1) + re ) ) .AND.
+     $               (rpsinz .GE. ( zd(j) + re )) ) id = j
+ 10           CONTINUE
+           END IF
+ 
+           DO 20 j = 1, id
+
+             sm = 1.0
+             IF(j .EQ. id .AND. id .EQ. i .AND. zen .GT. 90.0)
+     $          sm = -1.0
+ 
+             rj = re + zd(j-1)
+             rjp1 = re + zd(j)
+ 
+             dhj = zd(j-1) - zd(j)
+ 
+             ga = rj*rj - rpsinz*rpsinz
+             gb = rjp1*rjp1 - rpsinz*rpsinz
+             IF (ga .LT. 0.0) ga = 0.0
+             IF (gb .LT. 0.0) gb = 0.0
+ 
+             IF(id.GT.i .AND. j.EQ.id) THEN
+                dsj = SQRT( ga )
+             ELSE
+                dsj = SQRT( ga ) - sm*SQRT( gb )
+             END IF
+             dsdh(i,j) = dsj / dhj
+
+ 20        CONTINUE
+ 
+           nid(i) = id
+ 
+        END IF
+
+ 100  CONTINUE
+c
+      RETURN
+      END
+
+!==============================================================================
+
+      SUBROUTINE rtlink(nlev, nw, iw, ag, zen, dsdh, nid, dtrl, 
+     $                  dagas, dtcld, omcld, gcld, dtaer, omaer, gaer,
+     $                  edir, edn, eup, fdir, fdn, fup)
+
+      implicit none
+
+* input
+
+      integer :: nw                                  ! number of wavelength grid points
+      INTEGER nlev, iw
+      REAL ag
+      REAL zen
+      REAL dsdh(0:nlev,nlev)
+      INTEGER nid(0:nlev)
+
+      REAL dtrl(nlev,nw)
+      REAL dagas(nlev,nw)
+      REAL dtcld(nlev,nw), omcld(nlev,nw), gcld(nlev,nw)
+      REAL dtaer(nlev,nw), omaer(nlev,nw), gaer(nlev,nw)
+
+* output
+
+      REAL edir(nlev), edn(nlev), eup(nlev)
+      REAL fdir(nlev), fdn(nlev), fup(nlev)
+
+* local:
+
+      REAL dt(nlev), om(nlev), g(nlev)
+      REAL dtabs,dtsct,dscld,dsaer,dacld,daaer
+      INTEGER i, ii
+      real, parameter :: largest = 1.e+36
+
+* specific two ps2str
+
+      REAL ediri(nlev), edni(nlev), eupi(nlev)
+      REAL fdiri(nlev), fdni(nlev), fupi(nlev)
+
+      logical, save :: delta = .true.
+
+*_______________________________________________________________________
+
+* initialize:
+
+      do i = 1, nlev
+         fdir(i) = 0.
+         fup(i)  = 0.
+         fdn(i)  = 0.
+         edir(i) = 0.
+         eup(i)  = 0.
+         edn(i)  = 0.
+      end do
+c
+      do i = 1, nlev - 1
+         dscld = dtcld(i,iw)*omcld(i,iw)
+         dacld = dtcld(i,iw)*(1.-omcld(i,iw))
+c
+         dsaer = dtaer(i,iw)*omaer(i,iw)
+         daaer = dtaer(i,iw)*(1.-omaer(i,iw))
+c
+         dtsct = dtrl(i,iw) + dscld + dsaer
+         dtabs = dagas(i,iw) + dacld + daaer
+c
+         dtabs = amax1(dtabs,1./largest)
+         dtsct = amax1(dtsct,1./largest)
+c
+c     invert z-coordinate:
+c
+         ii = nlev - i
+         dt(ii) = dtsct + dtabs
+         om(ii) = dtsct/(dtsct + dtabs)
+         IF(dtsct .EQ. 1./largest) om(ii) = 1./largest
+         g(ii) = (gcld(i,iw)*dscld + 
+     $            gaer(i,iw)*dsaer)/dtsct
+      end do
+c
+c     call rt routine:
+c
+      call ps2str(nlev, zen, ag, dt, om, g,
+     $            dsdh, nid, delta,
+     $            fdiri, fupi, fdni, ediri, eupi, edni)
+c
+c     output (invert z-coordinate)
+c
+      do i = 1, nlev
+         ii = nlev - i + 1
+         fdir(i) = fdiri(ii)
+         fup(i) = fupi(ii)
+         fdn(i) = fdni(ii)
+         edir(i) = ediri(ii)
+         eup(i) = eupi(ii)
+         edn(i) = edni(ii)
+      end do
+c
+      return
+      end
+
+*=============================================================================*
+
+      SUBROUTINE ps2str(nlev,zen,rsfc,tauu,omu,gu,
+     $                  dsdh, nid, delta,
+     $                  fdr, fup, fdn, edr, eup, edn)
+
+*-----------------------------------------------------------------------------*
+*=  PURPOSE:                                                                 =*
+*=  Solve two-stream equations for multiple layers.  The subroutine is based =*
+*=  on equations from:  Toon et al., J.Geophys.Res., v94 (D13), Nov 20, 1989.=*
+*=  It contains 9 two-stream methods to choose from.  A pseudo-spherical     =*
+*=  correction has also been added.                                          =*
+*-----------------------------------------------------------------------------*
+*=  PARAMETERS:                                                              =*
+*=  NLEVEL  - INTEGER, number of specified altitude levels in the working (I)=*
+*=            grid                                                           =*
+*=  ZEN     - REAL, solar zenith angle (degrees)                          (I)=*
+*=  RSFC    - REAL, surface albedo at current wavelength                  (I)=*
+*=  TAUU    - REAL, unscaled optical depth of each layer                  (I)=*
+*=  OMU     - REAL, unscaled single scattering albedo of each layer       (I)=*
+*=  GU      - REAL, unscaled asymmetry parameter of each layer            (I)=*
+*=  DSDH    - REAL, slant path of direct beam through each layer crossed  (I)=*
+*=            when travelling from the top of the atmosphere to layer i;     =*
+*=            DSDH(i,j), i = 0..NZ-1, j = 1..NZ-1                            =*
+*=  NID     - INTEGER, number of layers crossed by the direct beam when   (I)=*
+*=            travelling from the top of the atmosphere to layer i;          =*
+*=            NID(i), i = 0..NZ-1                                            =*
+*=  DELTA   - LOGICAL, switch to use delta-scaling                        (I)=*
+*=            .TRUE. -> apply delta-scaling                                  =*
+*=            .FALSE.-> do not apply delta-scaling                           =*
+*=  FDR     - REAL, contribution of the direct component to the total     (O)=*
+*=            actinic flux at each altitude level                            =*
+*=  FUP     - REAL, contribution of the diffuse upwelling component to    (O)=*
+*=            the total actinic flux at each altitude level                  =*
+*=  FDN     - REAL, contribution of the diffuse downwelling component to  (O)=*
+*=            the total actinic flux at each altitude level                  =*
+*=  EDR     - REAL, contribution of the direct component to the total     (O)=*
+*=            spectral irradiance at each altitude level                     =*
+*=  EUP     - REAL, contribution of the diffuse upwelling component to    (O)=*
+*=            the total spectral irradiance at each altitude level           =*
+*=  EDN     - REAL, contribution of the diffuse downwelling component to  (O)=*
+*=            the total spectral irradiance at each altitude level           =*
+*-----------------------------------------------------------------------------*
+c
+      implicit none
+c
+*******
+* input:
+*******
+      INTEGER nlev
+      REAL zen, rsfc
+      REAL tauu(nlev), omu(nlev), gu(nlev)
+      REAL dsdh(0:nlev,nlev)
+      INTEGER nid(0:nlev)
+      LOGICAL delta
+
+*******
+* output:
+*******
+      REAL fup(nlev),fdn(nlev),fdr(nlev)
+      REAL eup(nlev),edn(nlev),edr(nlev)
+
+*******
+* local:
+*******
+      REAL tausla(0:nlev), tauc(0:nlev)
+      REAL mu2(0:nlev), mu, sum
+
+* internal coefficients and matrix
+      REAL lam(nlev),taun(nlev),bgam(nlev)
+      REAL e1(nlev),e2(nlev),e3(nlev),e4(nlev)
+      REAL cup(nlev),cdn(nlev),cuptn(nlev),cdntn(nlev)
+      REAL mu1(nlev)
+      INTEGER row
+      REAL a(2*nlev),b(2*nlev),d(2*nlev),e(2*nlev),y(2*nlev)
+
+*******
+* other:
+*******
+      REAL pifs, fdn0
+      REAL gi(nlev), omi(nlev), tempg
+      REAL f, g, om
+      REAL gam1, gam2, gam3, gam4
+      real, parameter :: largest = 1.e+36
+      real, parameter :: precis = 1.e-7
+
+* For calculations of Associated Legendre Polynomials for GAMA1,2,3,4
+* in delta-function, modified quadrature, hemispheric constant,
+* Hybrid modified Eddington-delta function metods, p633,Table1.
+* W.E.Meador and W.R.Weaver, GAS,1980,v37,p.630
+* W.J.Wiscombe and G.W. Grams, GAS,1976,v33,p2440, 
+* uncomment the following two lines and the appropriate statements further
+* down.
+C     REAL YLM0, YLM2, YLM4, YLM6, YLM8, YLM10, YLM12, YLMS, BETA0,
+C    >     BETA1, BETAn, amu1, subd
+
+      REAL expon, expon0, expon1, divisr, temp, up, dn
+      REAL ssfc
+      INTEGER nlayer, mrows, lev
+
+      INTEGER i, j
+
+* Some additional program constants:
+
+      real pi, dr
+      REAL eps
+      PARAMETER (eps = 1.E-3)
+*_______________________________________________________________________
+
+* MU = cosine of solar zenith angle
+* RSFC = surface albedo
+* TAUU =  unscaled optical depth of each layer
+* OMU  =  unscaled single scattering albedo
+* GU   =  unscaled asymmetry factor
+* KLEV = max dimension of number of layers in atmosphere
+* NLAYER = number of layers in the atmosphere
+* NLEVEL = nlayer + 1 = number of levels
+
+* initial conditions:  pi*solar flux = 1;  diffuse incidence = 0
+
+      pifs = 1.      
+      fdn0 = 0.
+
+      nlayer = nlev - 1
+
+      pi = acos(-1.)
+      dr = pi/180.
+      mu = COS(zen*dr)
+
+************** compute coefficients for each layer:
+* GAM1 - GAM4 = 2-stream coefficients, different for different approximations
+* EXPON0 = calculation of e when TAU is zero
+* EXPON1 = calculation of e when TAU is TAUN
+* CUP and CDN = calculation when TAU is zero
+* CUPTN and CDNTN = calc. when TAU is TAUN
+* DIVISR = prevents division by zero
+
+      do j = 0, nlev
+         tauc(j) = 0.
+         tausla(j) = 0.
+         mu2(j) = 1./SQRT(largest)
+      end do
+c
+      IF (.NOT. delta) THEN
+         DO i = 1, nlayer
+            gi(i) = gu(i)
+            omi(i) = omu(i)
+            taun(i) = tauu(i)
+         END DO
+      ELSE 
+
+* delta-scaling. Have to be done for delta-Eddington approximation, 
+* delta discrete ordinate, Practical Improved Flux Method, delta function,
+* and Hybrid modified Eddington-delta function methods approximations
+
+         DO i = 1, nlayer
+            f = gu(i)*gu(i)
+            gi(i) = (gu(i) - f)/(1 - f)
+            omi(i) = (1 - f)*omu(i)/(1 - omu(i)*f)       
+            taun(i) = (1 - omu(i)*f)*tauu(i)
+         END DO
+      END IF
+*
+* calculate slant optical depth at the top of the atmosphere when zen>90.
+* in this case, higher altitude of the top layer is recommended which can 
+* be easily changed in gridz.f.
+*
+      IF (zen .GT. 90.0) THEN
+         IF (nid(0) .LT. 0) THEN
+            tausla(0) = largest
+         ELSE
+            sum = 0.0
+            DO j = 1, nid(0)
+               sum = sum + 2.*taun(j)*dsdh(0,j)
+            END DO
+            tausla(0) = sum 
+         END IF
+      END IF
+*
+      DO 11, i = 1, nlayer
+          g = gi(i)
+          om = omi(i)
+          tauc(i) = tauc(i-1) + taun(i)
+
+* stay away from 1 by precision.  For g, also stay away from -1
+
+          tempg = AMIN1(abs(g),1. - precis)
+          g = SIGN(tempg,g)
+          om = AMIN1(om,1.-precis)
+
+* calculate slant optical depth
+*              
+          IF (nid(i) .LT. 0) THEN
+             tausla(i) = largest
+          ELSE
+             sum = 0.0
+             DO j = 1, MIN(nid(i),i)
+                sum = sum + taun(j)*dsdh(i,j)
+             END DO
+             DO j = MIN(nid(i),i)+1,nid(i)
+                sum = sum + 2.*taun(j)*dsdh(i,j)
+             END DO
+             tausla(i) = sum 
+             IF (tausla(i) .EQ. tausla(i-1)) THEN
+                mu2(i) = SQRT(largest)
+             ELSE
+                mu2(i) = (tauc(i)-tauc(i-1))/(tausla(i)-tausla(i-1))
+                mu2(i) = SIGN( AMAX1(ABS(mu2(i)),1./SQRT(largest)),
+     $                     mu2(i) )
+             END IF
+          END IF
+*
+*** the following gamma equations are from pg 16,289, Table 1
+*** save mu1 for each approx. for use in converting irradiance to actinic flux
+
+* Eddington approximation(Joseph et al., 1976, JAS, 33, 2452):
+
+c       gam1 =  (7. - om*(4. + 3.*g))/4.
+c       gam2 = -(1. - om*(4. - 3.*g))/4.
+c       gam3 = (2. - 3.*g*mu)/4.
+c       gam4 = 1. - gam3
+c       mu1(i) = 0.5
+
+* quadrature (Liou, 1973, JAS, 30, 1303-1326; 1974, JAS, 31, 1473-1475):
+
+c          gam1 = 1.7320508*(2. - om*(1. + g))/2.
+c          gam2 = 1.7320508*om*(1. - g)/2.
+c          gam3 = (1. - 1.7320508*g*mu)/2.
+c          gam4 = 1. - gam3
+c          mu1(i) = 1./sqrt(3.)
+         
+* hemispheric mean (Toon et al., 1089, JGR, 94, 16287):
+
+           gam1 = 2. - om*(1. + g)
+           gam2 = om*(1. - g)
+           gam3 = (2. - g*mu)/4.
+           gam4 = 1. - gam3
+           mu1(i) = 0.5
+
+* PIFM  (Zdunkovski et al.,1980, Conrib.Atmos.Phys., 53, 147-166):
+c         GAM1 = 0.25*(8. - OM*(5. + 3.*G))
+c         GAM2 = 0.75*OM*(1.-G)
+c         GAM3 = 0.25*(2.-3.*G*MU)
+c         GAM4 = 1. - GAM3
+c         mu1(i) = 0.5
+
+* delta discrete ordinates  (Schaller, 1979, Contrib.Atmos.Phys, 52, 17-26):
+c         GAM1 = 0.5*1.7320508*(2. - OM*(1. + G))
+c         GAM2 = 0.5*1.7320508*OM*(1.-G)
+c         GAM3 = 0.5*(1.-1.7320508*G*MU)
+c         GAM4 = 1. - GAM3
+c         mu1(i) = 1./sqrt(3.)
+
+* Calculations of Associated Legendre Polynomials for GAMA1,2,3,4
+* in delta-function, modified quadrature, hemispheric constant,
+* Hybrid modified Eddington-delta function metods, p633,Table1.
+* W.E.Meador and W.R.Weaver, GAS,1980,v37,p.630
+* W.J.Wiscombe and G.W. Grams, GAS,1976,v33,p2440
+c      YLM0 = 2.
+c      YLM2 = -3.*G*MU
+c      YLM4 = 0.875*G**3*MU*(5.*MU**2-3.)
+c      YLM6=-0.171875*G**5*MU*(15.-70.*MU**2+63.*MU**4)
+c     YLM8=+0.073242*G**7*MU*(-35.+315.*MU**2-693.*MU**4
+c    *+429.*MU**6)
+c     YLM10=-0.008118*G**9*MU*(315.-4620.*MU**2+18018.*MU**4
+c    *-25740.*MU**6+12155.*MU**8)
+c     YLM12=0.003685*G**11*MU*(-693.+15015.*MU**2-90090.*MU**4
+c    *+218790.*MU**6-230945.*MU**8+88179.*MU**10)
+c      YLMS=YLM0+YLM2+YLM4+YLM6+YLM8+YLM10+YLM12
+c      YLMS=0.25*YLMS
+c      BETA0 = YLMS
+c
+c         amu1=1./1.7320508
+c      YLM0 = 2.
+c      YLM2 = -3.*G*amu1
+c      YLM4 = 0.875*G**3*amu1*(5.*amu1**2-3.)
+c      YLM6=-0.171875*G**5*amu1*(15.-70.*amu1**2+63.*amu1**4)
+c     YLM8=+0.073242*G**7*amu1*(-35.+315.*amu1**2-693.*amu1**4
+c    *+429.*amu1**6)
+c     YLM10=-0.008118*G**9*amu1*(315.-4620.*amu1**2+18018.*amu1**4
+c    *-25740.*amu1**6+12155.*amu1**8)
+c     YLM12=0.003685*G**11*amu1*(-693.+15015.*amu1**2-90090.*amu1**4
+c    *+218790.*amu1**6-230945.*amu1**8+88179.*amu1**10)
+c      YLMS=YLM0+YLM2+YLM4+YLM6+YLM8+YLM10+YLM12
+c      YLMS=0.25*YLMS
+c      BETA1 = YLMS
+c
+c         BETAn = 0.25*(2. - 1.5*G-0.21875*G**3-0.085938*G**5
+c    *-0.045776*G**7)
+
+
+* Hybrid modified Eddington-delta function(Meador and Weaver,1980,JAS,37,630):
+c         subd=4.*(1.-G*G*(1.-MU))
+c         GAM1 = (7.-3.*G*G-OM*(4.+3.*G)+OM*G*G*(4.*BETA0+3.*G))/subd
+c         GAM2 =-(1.-G*G-OM*(4.-3.*G)-OM*G*G*(4.*BETA0+3.*G-4.))/subd
+c         GAM3 = BETA0
+c         GAM4 = 1. - GAM3
+c         mu1(i) = (1. - g*g*(1.- mu) )/(2. - g*g)
+
+*****
+* delta function  (Meador, and Weaver, 1980, JAS, 37, 630):
+c         GAM1 = (1. - OM*(1. - beta0))/MU
+c         GAM2 = OM*BETA0/MU
+c         GAM3 = BETA0
+c         GAM4 = 1. - GAM3
+c         mu1(i) = mu
+*****
+* modified quadrature (Meador, and Weaver, 1980, JAS, 37, 630):
+c         GAM1 = 1.7320508*(1. - OM*(1. - beta1))
+c         GAM2 = 1.7320508*OM*beta1
+c         GAM3 = BETA0
+c         GAM4 = 1. - GAM3
+c         mu1(i) = 1./sqrt(3.)
+
+* hemispheric constant (Toon et al., 1989, JGR, 94, 16287):
+c         GAM1 = 2.*(1. - OM*(1. - betan))
+c         GAM2 = 2.*OM*BETAn
+c         GAM3 = BETA0
+c         GAM4 = 1. - GAM3
+c         mu1(i) = 0.5
+
+*****
+
+* lambda = pg 16,290 equation 21
+* big gamma = pg 16,290 equation 22
+* if gam2 = 0., then bgam = 0. 
+
+         lam(i) = sqrt(gam1*gam1 - gam2*gam2)
+
+         IF (gam2 .NE. 0.) THEN
+            bgam(i) = (gam1 - lam(i))/gam2
+         ELSE
+            bgam(i) = 0.
+         END IF
+
+         expon = EXP(-lam(i)*taun(i))
+
+* e1 - e4 = pg 16,292 equation 44
+         
+         e1(i) = 1. + bgam(i)*expon
+         e2(i) = 1. - bgam(i)*expon
+         e3(i) = bgam(i) + expon
+         e4(i) = bgam(i) - expon
+
+* the following sets up for the C equations 23, and 24
+* found on page 16,290
+* prevent division by zero (if LAMBDA=1/MU, shift 1/MU^2 by EPS = 1.E-3
+* which is approx equiv to shifting MU by 0.5*EPS* (MU)**3
+
+         expon0 = EXP(-tausla(i-1))
+         expon1 = EXP(-tausla(i))
+          
+         divisr = lam(i)*lam(i) - 1./(mu2(i)*mu2(i))
+         temp = AMAX1(eps,abs(divisr))
+         divisr = SIGN(temp,divisr)
+
+         up = om*pifs*((gam1 - 1./mu2(i))*gam3 + gam4*gam2)/divisr
+         dn = om*pifs*((gam1 + 1./mu2(i))*gam4 + gam2*gam3)/divisr
+         
+* cup and cdn are when tau is equal to zero
+* cuptn and cdntn are when tau is equal to taun
+
+         cup(i) = up*expon0
+         cdn(i) = dn*expon0
+         cuptn(i) = up*expon1
+         cdntn(i) = dn*expon1
+ 
+   11 CONTINUE
+
+***************** set up matrix ******
+* ssfc = pg 16,292 equation 37  where pi Fs is one (unity).
+
+      ssfc = rsfc*mu*EXP(-tausla(nlayer))*pifs
+
+* MROWS = the number of rows in the matrix
+
+      mrows = 2*nlayer     
+      
+* the following are from pg 16,292  equations 39 - 43.
+* set up first row of matrix:
+
+      i = 1
+      a(1) = 0.
+      b(1) = e1(i)
+      d(1) = -e2(i)
+      e(1) = fdn0 - cdn(i)
+
+      row=1
+
+* set up odd rows 3 thru (MROWS - 1):
+
+      i = 0
+      DO 20, row = 3, mrows - 1, 2
+         i = i + 1
+         a(row) = e2(i)*e3(i) - e4(i)*e1(i)
+         b(row) = e1(i)*e1(i + 1) - e3(i)*e3(i + 1)
+         d(row) = e3(i)*e4(i + 1) - e1(i)*e2(i + 1)
+         e(row) = e3(i)*(cup(i + 1) - cuptn(i)) + 
+     $        e1(i)*(cdntn(i) - cdn(i + 1))
+   20 CONTINUE
+
+* set up even rows 2 thru (MROWS - 2): 
+
+      i = 0
+      DO 30, row = 2, mrows - 2, 2
+         i = i + 1
+         a(row) = e2(i + 1)*e1(i) - e3(i)*e4(i + 1)
+         b(row) = e2(i)*e2(i + 1) - e4(i)*e4(i + 1)
+         d(row) = e1(i + 1)*e4(i + 1) - e2(i + 1)*e3(i + 1)
+         e(row) = (cup(i + 1) - cuptn(i))*e2(i + 1) - 
+     $        (cdn(i + 1) - cdntn(i))*e4(i + 1)
+   30 CONTINUE
+
+* set up last row of matrix at MROWS:
+
+      row = mrows
+      i = nlayer
+      
+      a(row) = e1(i) - rsfc*e3(i)
+      b(row) = e2(i) - rsfc*e4(i)
+      d(row) = 0.
+      e(row) = ssfc - cuptn(i) + rsfc*cdntn(i)
+
+* solve tri-diagonal matrix:
+
+      CALL tridiag(a, b, d, e, y, mrows)
+
+**** unfold solution of matrix, compute output fluxes:
+
+      row = 1 
+      lev = 1
+      j = 1
+      
+* the following equations are from pg 16,291  equations 31 & 32
+
+      fdr(lev) = EXP( -tausla(0) )
+      edr(lev) = mu * fdr(lev)
+      edn(lev) = fdn0
+      eup(lev) =  y(row)*e3(j) - y(row + 1)*e4(j) + cup(j)
+      fdn(lev) = edn(lev)/mu1(lev)
+      fup(lev) = eup(lev)/mu1(lev)
+
+      DO 60, lev = 2, nlayer + 1
+         fdr(lev) = EXP(-tausla(lev-1))
+         edr(lev) =  mu *fdr(lev)
+         edn(lev) =  y(row)*e3(j) + y(row + 1)*e4(j) + cdntn(j)
+         eup(lev) =  y(row)*e1(j) + y(row + 1)*e2(j) + cuptn(j)
+         fdn(lev) = edn(lev)/mu1(j)
+         fup(lev) = eup(lev)/mu1(j)
+
+         row = row + 2
+         j = j + 1
+   60 CONTINUE
+c
+      RETURN
+      END
+
+*=============================================================================*
+
+      SUBROUTINE tridiag(a,b,c,r,u,n)
+
+*_______________________________________________________________________
+* solves tridiagonal system.  From Numerical Recipies, p. 40
+*_______________________________________________________________________
+
+      IMPLICIT NONE
+
+* input:
+      INTEGER n
+      REAL a, b, c, r
+      DIMENSION a(n),b(n),c(n),r(n)
+
+* output:
+      REAL u
+      DIMENSION u(n)
+
+* local:
+      INTEGER j
+
+      REAL bet, gam
+      DIMENSION gam(n)
+*_______________________________________________________________________
+
+      IF (b(1) .EQ. 0.) STOP 1001
+      bet   = b(1)
+      u(1) = r(1)/bet
+      DO 11, j = 2, n   
+         gam(j) = c(j - 1)/bet
+         bet = b(j) - a(j)*gam(j)
+         IF (bet .EQ. 0.) STOP 2002 
+         u(j) = (r(j) - a(j)*u(j - 1))/bet
+   11 CONTINUE
+      DO 12, j = n - 1, 1, -1  
+         u(j) = u(j) - gam(j + 1)*u(j + 1)
+   12 CONTINUE
+*_______________________________________________________________________
+
+      RETURN
+      END
