#define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr
#define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if
#define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if
!
#include "tm5.inc"
#include "output.inc"  
!
!----------------------------------------------------------------------------
!                    TM5                                                    !
!----------------------------------------------------------------------------
!BOP
!
! !MODULE:      PHOTOLYSIS
!
! !DESCRIPTION: <add description>
!\\
!\\
! !INTERFACE: 
!
MODULE PHOTOLYSIS
  !
  ! !USES:
  !
  use GO,             only : gol, goPr, goErr
  use tm5_distgrid,   only : dgrid, Get_DistGrid, gather
  use photolysis_data
#ifdef with_optics
  use optics,         only : wavelendep
#endif

  IMPLICIT NONE
  PRIVATE
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  public :: photolysis_init, photolysis_done ! life cycle
  public :: ozone_info_online ! calculate the overhead ozone
  public :: aerosol_info      ! aerosol optical depths
  public :: slingo            ! cloud optical properties
  public :: update_csqy       ! update ...
  public :: daysim            ! calculate solar zenith angles
  public :: get_sza           ! return solar zenith angle (function)
  public :: photo_flux        ! calculates radiation fields
  public :: photorates_tropo  ! calculates photolysis and heating rates
#ifdef with_optics
  public :: aerosol_info_m7   ! aerosol optical depths
#endif
  !
  ! !PUBLIC DATA MEMBERS:
  !
  public                     :: phot_dat   ! type defined in photolysis_data
  !
  ! !PRIVATE DATA MEMBERS:
  !
  character(len=*), parameter  ::  mname = 'photolysis'

  real, parameter      :: todu = 3.767e-17                                ! from #/cm2 --> du
  real, parameter      :: to_m = 3.767e-22                                ! from #/cm2 --> m (for o2)
  real, parameter      :: sp   = 6.022e23*1.e-4*0.2095/(28.964*1e-3*9.81) ! sp from pa ---> #o2/cm2   

  logical                :: with_o3du

  real,dimension(122,34) :: xs_o3_look
  real,dimension( 65,34) :: xs_hno3_look,xs_pan_look,qy_o3_look 
  real,dimension( 45,34) :: xs_h2o2_look
  real,dimension( 55,34) :: xs_n2o5_look
  real,dimension( 89,34) :: xs_no2_look,qy_no2_look,qy_co_look
  real,dimension( 62,34) :: xs_no3_look
  real,dimension(105,34) :: xs_ch2o_look  ! local
  ! constants for acetone qy
  real,dimension(89,34)    ::   a1_acet,a2_acet,a3_acet,a4_acet
    
  integer                 :: l  
  real                    :: wl,x, ww

#ifdef with_optics
  type(wavelendep),dimension(:),allocatable :: wdep 
#endif

  !
  ! !REVISION HISTORY: 
  !   16 Jun 2011 - P. Le Sager - new version from JEW
  !   21 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
  !    3 Oct 2012 - P. Le Sager - separate Init into true init and monthly update
  !    4 Mar 2013 - P. Le Sager - changes for optics feedback from NB
  !
  ! !REMARKS:
  !
  !EOP
  !--------------------------------------------------------------------------

CONTAINS

  ! Trap overflow
  function safe_div(n,d,altv) result(q)
    real, intent(in) :: n, d, altv
    real :: q

    if ( exponent(n)-exponent(d) >= maxexponent(n) .OR. d==0) then
       q=altv
    else
       q=n/d
    end if
  end function safe_div

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    INIT_PHOTOLYSIS
  !
  ! !DESCRIPTION:  read reference tables for photolysis, allocate photolysis
  !                data. Called from sources_sinks/Trace0 at the run start and
  !                at beginning of every month. 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE PHOTOLYSIS_INIT( first, iunit, o3du )
    !
    ! !USES:
    !
    use dims,         only : im, jm, lm, nregions, newsrun
    use ParTools,     only : isRoot, par_broadcast
    use global_types, only : d23_data
    use global_data,  only : rcfile
    use MeteoData,    only : Set
    use MeteoData,    only : sp_dat, temper_dat, lsmask_dat, phlb_dat
    use MeteoData,    only : lwc_dat, iwc_dat, cc_dat, humid_dat, gph_dat
    use GO,           only : TrcFile, Init, Done, ReadRc
#ifdef with_optics
    use optics,       only : Optics_Init
#endif
    !
    ! !INPUT PARAMETERS:
    !
    logical,                                       intent(in) :: first   ! true init for a new run?
    integer,                                       intent(in) :: iunit
    type(d23_data), dimension(nregions), optional, intent(in) :: o3du 
    !
    ! !REVISION HISTORY: 
    !   21 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    integer            :: i,j,l,m,n,region,imr,jmr,lmr,i1, i2, j1, j2
    real               :: lat
    integer            :: status, wav
    type(TrcFile)      :: rcF
    character(len=256) :: photodir

    character(len=*), parameter  ::  rname = mname//'/photolysis_init'

    ! --- begin -------------------------

    with_o3du = present(o3du)

    !--------------------------------------------------
    ! ** TRUE INIT : only once
    !--------------------------------------------------        
    if (first) then

       REG: do region = 1, nregions

          call Set(     sp_dat(region), status, used=.true. ) 
          call Set( temper_dat(region), status, used=.true. )
          call Set(  humid_dat(region), status, used=.true. )
          call Set(    gph_dat(region), status, used=.true. )
          call Set(   phlb_dat(region), status, used=.true. )
          call Set(    lwc_dat(region), status, used=.true. ) 
          call Set(    iwc_dat(region), status, used=.true. ) 
          call Set(     cc_dat(region), status, used=.true. ) 
          call Set( lsmask_dat(region), status, used=.true. ) 

          ! allocate memory photolysis and initialise:
          call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
          lmr=lm(region)

          if (with_o3du) allocate( phot_dat(region)%o3klim_top(j1:j2) )

          allocate(phot_dat(region)%albedo    (i1:i2, j1:j2))     ; phot_dat(region)%albedo  = 0.05
          allocate(phot_dat(region)%ags_av    (i1:i2, j1:j2))     ; phot_dat(region)%ags_av  = 0.0
          allocate(phot_dat(region)%sza_av    (i1:i2, j1:j2))     ; phot_dat(region)%sza_av  = 0.0
          allocate(phot_dat(region)%lwp_av    (i1:i2, j1:j2,lmr)) ; phot_dat(region)%lwp_av  = 0.0
          allocate(phot_dat(region)%cfr_av    (i1:i2, j1:j2,lmr)) ; phot_dat(region)%cfr_av  = 0.0
          allocate(phot_dat(region)%reff_av   (i1:i2, j1:j2,lmr)) ; phot_dat(region)%reff_av = 0.0

          allocate(phot_dat(region)%aero_surface_clim(i1:i2, j1:j2)) ; phot_dat(region)%aero_surface_clim = 0

          phot_dat(region)%nalb_av   = 0.
          phot_dat(region)%ncloud_av = 0.
          phot_dat(region)%naer_av   = 0.

          allocate(phot_dat(region)%cloud_reff(i1:i2, j1:j2, lmr )) ; phot_dat(region)%cloud_reff = 0.0
          allocate(phot_dat(region)%pmcld     (i1:i2, j1:j2, lmr )) ; phot_dat(region)%pmcld      = 0.0
          allocate(phot_dat(region)%taus_cld  (i1:i2, j1:j2, lmr )) ; phot_dat(region)%taus_cld   = 0.0
          allocate(phot_dat(region)%taua_cld  (i1:i2, j1:j2, lmr )) ; phot_dat(region)%taua_cld   = 0.0 

          allocate(phot_dat(region)%pmaer     (i1:i2, j1:j2, lmr, nbands_trop, ngrid )) ; phot_dat(region)%pmaer    = 0.0
          allocate(phot_dat(region)%taus_aer  (i1:i2, j1:j2, lmr, nbands_trop, ngrid )) ; phot_dat(region)%taus_aer = 0.0
          allocate(phot_dat(region)%taua_aer  (i1:i2, j1:j2, lmr, nbands_trop, ngrid )) ; phot_dat(region)%taua_aer = 0.0         

          allocate(phot_dat(region)%pmcld_av    (i1:i2, j1:j2, lmr )) ; phot_dat(region)%pmcld_av    = 0.0
          allocate(phot_dat(region)%taus_cld_av (i1:i2, j1:j2, lmr )) ; phot_dat(region)%taus_cld_av = 0.0
          allocate(phot_dat(region)%taua_cld_av (i1:i2, j1:j2, lmr )) ; phot_dat(region)%taua_cld_av = 0.0 

          allocate(phot_dat(region)%pmaer_av    (i1:i2, j1:j2, lmr, nbands_trop,ngrid )) ; phot_dat(region)%pmaer_av    = 0.0
          allocate(phot_dat(region)%taus_aer_av (i1:i2, j1:j2, lmr, nbands_trop,ngrid )) ; phot_dat(region)%taus_aer_av = 0.0
          allocate(phot_dat(region)%taua_aer_av (i1:i2, j1:j2, lmr, nbands_trop,ngrid )) ; phot_dat(region)%taua_aer_av = 0.0   

          allocate(phot_dat(region)%jo3      (i1:i2, j1:j2, lmr )) ; phot_dat(region)%jo3       = 0.0
          allocate(phot_dat(region)%jno2     (i1:i2, j1:j2, lmr )) ; phot_dat(region)%jno2      = 0.0
          allocate(phot_dat(region)%jo3_av   (i1:i2, j1:j2, lmr )) ; phot_dat(region)%jo3_av    = 0.0
          allocate(phot_dat(region)%jno2_av  (i1:i2, j1:j2, lmr )) ; phot_dat(region)%jno2_av   = 0.0
          allocate(phot_dat(region)%jhno2    (i1:i2, j1:j2, lmr )) ; phot_dat(region)%jhno2     = 0.0
          allocate(phot_dat(region)%jhno2_av (i1:i2, j1:j2, lmr )) ; phot_dat(region)%jhno2_av  = 0.0
          allocate(phot_dat(region)%jch2oa   (i1:i2, j1:j2, lmr )) ; phot_dat(region)%jch2oa    = 0.0
          allocate(phot_dat(region)%jch2oa_av(i1:i2, j1:j2, lmr )) ; phot_dat(region)%jch2oa_av = 0.0
          allocate(phot_dat(region)%jch2ob   (i1:i2, j1:j2, lmr )) ; phot_dat(region)%jch2ob    = 0.0
          allocate(phot_dat(region)%jch2ob_av(i1:i2, j1:j2, lmr )) ; phot_dat(region)%jch2ob_av = 0.0             

          phot_dat(region)%nj_av = 0

          ! absorption cross-sections
          allocate(phot_dat(region)%v2         (i1:i2, j1:j2, lmr+1  )) ; phot_dat(region)%v2          = 0.0
          allocate(phot_dat(region)%v3         (i1:i2, j1:j2, lmr+1  )) ; phot_dat(region)%v3          = 0.0
          allocate(phot_dat(region)%v3_surf    (i1:i2, j1:j2         )) ; phot_dat(region)%v3_surf     = 0.0
          allocate(phot_dat(region)%qy_ch3cocho(i1:i2, j1:j2, lmr, 52)) ; phot_dat(region)%qy_ch3cocho = 0.0
          allocate(phot_dat(region)%qy_c2o3    (i1:i2, j1:j2, lmr, 89)) ; phot_dat(region)%qy_c2o3     = 0.0

          ! SAD fields
          phot_dat(region)%nsad_av   = 0
          allocate(phot_dat(region)%sad_cld    (i1:i2, j1:j2, lmr)) ; phot_dat(region)%sad_cld  = 0.0
          allocate(phot_dat(region)%sad_ice    (i1:i2, j1:j2, lmr)) ; phot_dat(region)%sad_ice  = 0.0
          allocate(phot_dat(region)%sad_aer    (i1:i2, j1:j2, lmr)) ; phot_dat(region)%sad_aer  = 0.0
          allocate(phot_dat(region)%sad_cld_av (i1:i2, j1:j2, lmr)) ; phot_dat(region)%sad_cld_av  = 0.0
          allocate(phot_dat(region)%sad_ice_av (i1:i2, j1:j2, lmr)) ; phot_dat(region)%sad_ice_av  = 0.0
          allocate(phot_dat(region)%sad_aer_av (i1:i2, j1:j2, lmr)) ; phot_dat(region)%sad_aer_av  = 0.0
          
       end do REG

       ! read settings from rcfile:
       call Init( rcF, rcfile, status )
       call ReadRc( rcF, 'input.photo', photodir, status )

       !--------------------------------------------------------------------   
       !     wavelengths in the middle of the intervals                        
       !--------------------------------------------------------------------   
       !     see: C. Bruehl, P.J. Crutzen, Scenarios of possible changes in ...
       !          Climate Dynamics (1988)2: 173-203 
       ! ONLINE TUNING : first 13 wavelength bins are now ignored. Grid also shortened above 690nm.

       do l=1,13 
          wave_full(l)=1./(56250.-500.*l) 
       end do

       do l=14,45 
          wave_full(l)=1./(49750.-(l-13)*500.) 
       end do

       do l=46,68 
          wave_full(l)=(266.+(l-13))*1.e-7 
       end do

       do l=69,71 
          wave_full(l)=(320.5+2.*(l-68))*1.e-7 
       end do
       do l=72,176 
          wave_full(l)=(325.+5.*(l-71))*1.e-7 
       end do

       do l=2,maxwav-1 
          dwave_full(l)=0.5*(wave_full(l+1)-wave_full(l-1))
       end do
       dwave_full(1)      = dwave_full(2) 
       dwave_full(maxwav) = dwave_full(maxwav-1) 
       ! select a subset from the full spctral grid to remove the first band      
       wave=wave_full(14:135)
       dwave=dwave_full(14:135)

       !---------------------------------------------------------------------  
       !     Rayleigh scattering cross sections                                
       !     emp. formula of Nicolet   plan. space sci.32, 1467f (1984)        
       !---------------------------------------------------------------------  

       do  l=1,maxwav                         
          wl=wave(l)*1.e4                   
          x=0.389*wl+0.09426/wl-0.3228      
          cs_ray(l)=4.02e-28/wl**(4.+x)     
       end do

       !--------------------------------------------------------------         
       !    Read and store temperature dependent cross-section values
       !--------------------------------------------------------------         

       open(unit=7,  file=trim(photodir)//'/tropo_look_up_modified_cb05.dat',     status='old')
       read(7,*)
       read(7,*)        
       read(7,*)
       do i=1,34
          read(7,'(176(1Pe10.3))') (xs_o3_look(wav,i),wav=1,122)
       enddo
       read(7,*)
       do i=1,34
          read(7,*) (xs_no2_look(wav,i),wav=1,89)
       enddo
       read(7,*)
       do i=1,34
          read(7,*) (xs_hno3_look(wav,i),wav=1,65)
       enddo
       read(7,*)
       do i=1,34
          read(7,*) (xs_h2o2_look(wav,i),wav=1,45)
       enddo
       read(7,*)
       do i=1,34
          read(7,*) (xs_n2o5_look(wav,i),wav=1,55)
       enddo
       read(7,*)
       do i=1,34
          read(7,*) (xs_ch2o_look(wav,i),wav=1,90)
       enddo
       read(7,*)
       do i=1,34
          read(7,*) (xs_pan_look(wav,i),wav=1,65)
       enddo
       read(7,*)
       do i=1,34
          read(7,*) (xs_no3_look(wav,i),wav=1,62)
       enddo
       read(7,*)
       do i=1,34
          read(7,*) (qy_o3_look(wav,i),wav=1,65)
       enddo
       read(7,*)
       do i=1,34
          read(7,*)(qy_no2_look(wav,i),wav=1,89)
       enddo
       read(7,*)
       do i=1,34
          read(7,*)(qy_co_look(wav,i),wav=1,89)
       enddo
       read(7,*)
       do i=1,34
          read(7,*)(a1_acet(wav,i),wav=1,89)
       enddo
       read(7,*)
       do i=1,34
          read(7,*)(a2_acet(wav,i),wav=1,89)
       enddo
       read(7,*)
       do i=1,34
          read(7,*)(a3_acet(wav,i),wav=1,89)
       enddo
       read(7,*)
       do i=1,34
          read(7,*)(a4_acet(wav,i),wav=1,89)
       enddo
       close(7)
       
       ! put constants in correct units
       a1_acet=a1_acet*1e-18
       a2_acet=a2_acet*1e-17
       a4_acet=a4_acet*1e-15        

       !--------------------------------------------------------------         
       !    Read and store temperature dependent cross-section values

       !------------------ extraterrestrial flux from input file ------------------

!        open(unit=7,  file=trim(photodir)//'/OMI.data.reduce',status='old')     
!        read(7,*)
!        read(7,'(1p7e10.2)') flux 
!        close(7)
       
       ! Cannot be read for some unknown reason on cca@ecmwf with cray compiler.
       ! Just set it here:

       flux=(/ 2.351E+12, 2.414E+12, 2.360E+12, 3.004E+12, 8.861E+12, 4.475E+12, 9.005E+12, & 
               1.607E+13, 1.556E+13, 1.221E+13, 2.061E+13, 2.070E+13, 1.047E+13, 7.497E+12, & 
               1.729E+13, 9.523E+12, 2.018E+13, 2.265E+13, 1.275E+13, 1.750E+13, 2.590E+13, & 
               8.777E+13, 2.193E+13, 5.449E+13, 1.074E+14, 6.856E+13, 7.875E+13, 2.275E+13, & 
               1.488E+14, 2.166E+14, 2.518E+14, 3.504E+14, 1.493E+14, 5.366E+13, 4.045E+13, & 
               2.054E+13, 6.274E+13, 9.589E+13, 1.326E+14, 1.307E+13, 1.332E+14, 1.729E+14, & 
               1.349E+14, 5.685E+13, 1.342E+14, 1.391E+14, 8.107E+13, 1.459E+14, 1.965E+14, & 
               5.681E+13, 1.349E+14, 6.819E+13, 2.137E+14, 1.564E+14, 1.334E+14, 3.559E+14, & 
               2.593E+14, 4.187E+14, 7.893E+14, 1.610E+14, 1.021E+15, 9.119E+14, 1.056E+15, & 
               8.424E+14, 6.622E+14, 9.387E+14, 1.463E+15, 3.382E+14, 1.076E+15, 9.586E+14, & 
               1.783E+15, 8.718E+14, 1.864E+15, 1.661E+15, 1.938E+15, 2.101E+15, 2.058E+15, & 
               1.738E+15, 9.227E+14, 2.404E+15, 2.176E+15, 2.625E+15, 2.167E+15, 1.741E+15, & 
               2.424E+15, 1.531E+15, 1.819E+15, 2.625E+15, 2.293E+15, 2.570E+15, 2.619E+15, & 
               2.663E+15, 2.705E+15, 2.542E+15, 1.331E+15, 2.633E+15, 2.575E+15, 2.744E+15, & 
               1.958E+15, 2.642E+15, 2.691E+15, 2.646E+15, 2.720E+15, 2.703E+15, 1.188E+15, & 
               2.720E+15, 2.271E+15, 2.405E+15, 2.723E+15, 2.690E+15, 2.550E+15, 2.670E+15, & 
               2.665E+15, 2.702E+15, 2.653E+15, 2.658E+15, 2.691E+15, 2.593E+15, 2.646E+15, & 
               2.653E+15, 2.636E+15, 2.648E+15 /)
       

#ifdef with_optics

       ! define wavelengths for optics calculations
       nwdep = nbands_trop + count(lmid.ne.lmid_gridA)
       wav_grid  = 0
       wav_gridA = 0
       allocate(photo_wavelengths(nwdep))

       l=1
       do i=1,nbands_trop
          if (lmid(i)==lmid_gridA(i)) then
             photo_wavelengths(l) = wave(lmid(i))*1.e4 ! cm to um
             wav_grid(i) = l
             wav_gridA(i) = l
             l=l+1   
          else
             photo_wavelengths(l) = wave(lmid(i))*1.e4 ! cm to um
             photo_wavelengths(l+1) = wave(lmid_gridA(i))*1.e4 ! cm to um
             wav_grid(i) = l
             wav_gridA(i) = l+1
             l=l+2
          endif
       enddo
              
       allocate(wdep(nwdep))
       wdep(:)%wl = photo_wavelengths
       wdep(:)%split = .false.
       wdep(:)%insitu = .false.
       
       call Done( rcF, status )
       IF_NOTOK_RETURN(status=1)
       
       CALL OPTICS_INIT( nwdep, wdep, status )
       IF_NOTOK_RETURN(status=1)
       
       deallocate(photo_wavelengths)
       
#else  
       !**************************************************************************
       !     aerosol data from: "Models for Aerosols of the Lower Atmosphere   
       !     and the Effects of Humidity Variations on their Optical Properties
       !     E. P. Shettle and R. W. Fenn (1979), Environmental Research Paper 
       !     No. 676                                                           
       !************************************************************************** 

       open(unit=7,  file=trim(photodir)//'/aerosol_reduce.dat',status='old')
       do l = 1,4
          read(7,*)
          do i = 1,8
             read(7,*)
             read(7,20)(ext(wav,i,l),wav=1,maxwav)        
             read(7,20)(abs_eff(wav,i,l),wav=1,maxwav)
             read(7,20)(g(wav,i,l),wav=1,maxwav)
             do wav=1,maxwav
                sca(wav,i,l) = ext(wav,i,l) - abs_eff(wav,i,l)
             enddo
          enddo
       enddo
       close(7) 

20     format(6(1X,F7.5))

       do l=1,maxwav
          do i=1,8
             do j=1,4
                sca(l,i,j) = sca(l,i,j)/pn_ref(j)*1.E-5
                abs_eff(l,i,j) = abs_eff(l,i,j)/pn_ref(j)*1.E-5
             end do
          end do
       end do

       ! close rc file
       call Done( rcF, status )
#endif
       
    ELSE ! NEWSRUN


       !--------------------------------------------------
       ! ** MONTHLY UPDATE 
       !--------------------------------------------------    

       if (with_o3du) then
          do region = 1, nregions             
             call Get_DistGrid( dgrid(region), J_STRT=j1, J_STOP=j2 )
             phot_dat(region)%o3klim_top(:) = o3du(region)%d23(j1:j2,lm(region))
          end do
       end if

    ENDIF

  END SUBROUTINE PHOTOLYSIS_INIT
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:   DAYSIM
  !
  ! !DESCRIPTION: Get all solar zenith angles (THA) for a given Julian day
  !                (DAYNR) at 1/6, 3/6(=center), and 5/6 of each grid box.
  !                THA is three times oversampled. 
  !                The longitudinal variation is equal to the simulation of
  !                one day, and covers the [-180,180] range for all regions.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE DAYSIM(region, daynr, is, i1,i2,j1, j2, tha)
    !
    ! !USES:
    !
    use dims,      only : im, jm, xbeg, xend
    use binas,     only : pi
    use meteodata, only : lli
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)   :: daynr, region, is, j1, j2,i1,i2
    !
    ! !OUTPUT PARAMETERS:
    !
    real,    intent(out)  :: tha( im(region)*3, j1:j2)
    !
    ! !REVISION HISTORY: 
    !   26 Apr 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:
    !   (1) The regions are artificially expanded to cover all longitudes, in
    !       order to keep the correspondance b/w time and longitude.
    !   (2) Need to be called once a day only. Hence separated from GET_SZA,
    !       and the need to have 'im_alt' and 'shift' defined module wide.
    !   (3) Input 'is' is not i1=I_START from the distributed grid, but: (i1-1)*3+1
    !       It is not used yet. To delete?
    !   (4) Considering the cshift functions used on THA in GET_SZA, it is
    !       assumed that THA covers the entire longitude range.
    !EOP
    !------------------------------------------------------------------------
    !BOC

    real, dimension(j1:j2)  :: lat
    real                    :: houra, rdecl, sinhh, eqr, ut
    real                    :: piby, lon, tz, sintz, costz
    real                    :: sin2tz, cos2tz, sin3tz, cos3tz
    integer                 :: i, j, k
    !
    !call Get_DistGrid( dgrid(region), I_STOP=i2, J_STOP=j2 )
    
    piby = pi/180.
    !
    ! latitude range in radian
    lat  = lli(region)%lat
    !
    tz = 2.*pi*daynr/365.
    !
    ! Solar declination approximation 
    !(p.97, Brasseur et al, Atmospheric Chemistry and Global Change, 1999)
    !
    rdecl = -0.4*cos(2.*pi*real(daynr+10)/365.)
    !
    sintz = sin(tz)
    costz = cos(tz)
    sin2tz = 2.*sintz*costz
    cos2tz = costz*costz-sintz*sintz

    eqr = 0.000075 + 0.001868*costz  - 0.032077*sintz &
         - 0.014615*cos2tz - 0.040849*sin2tz
    !
    ! the sza variation over longitude [-180,180] at ut=12 equals
    ! a with variation of ut over [0,24] at longitude 0:
    ! houra = 15.*(ut-12. + eqr*24./(2.*pi) + lon*24./360)*piby
    !
    
    ! THA (sampled at 1/6, 3/6, and 5/6 of each box). As follow for 6 degree resolution:
    !
    !     -180---+---+---177---+---+---174---+---+---171---+---+---168  .......etc...
    !       |                           |                           |
    !       |        box #1             |          box #2           |
    !       |                           |                           |
    !       +---------------------------+---------------------------+   .......etc.....
    !       0    1   2    3    4   5    0    1   2    3    4   5
    !            *        *        *         *        *        *         <---- samples (lon,THA)    
    do j=j1,j2
       do i=1,im(region)*3
          lon      = -180. + (real(i)-0.5)*360./real(im(region)*3)
          !lon      = xbeg(region) + (real(i)-0.5)*(xend(region)-xbeg(region))/real(im(region)*3)
          houra    = -pi + lon*piby + eqr
          sinhh    = sin(rdecl)*sin(lat(j)) &
               + cos(rdecl)*cos(lat(j))*cos(houra)
          tha(i,j) = acos(sinhh)/piby
       enddo
    enddo
    
  END SUBROUTINE DAYSIM
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    GET_SZA
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE GET_SZA(region, i1, j1, i2, j2, tstart, deltat, tha, sza)
    !
    ! !USES:
    !
    use dims, only : im, jm
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)  :: region, i1, j1, i2, j2
    real,    intent(in)  :: tstart                  ! start time of chemistry timestep     
    real,    intent(in)  :: deltat                  ! chemistry timestep in hours          
    real,    intent(in)  :: tha(im(region)*3,j1:j2) ! three times oversampled zenith angles
    !
    ! !OUTPUT PARAMETERS:
    !
    real,    intent(out) :: sza(i1:i2,j1:j2)     ! effective solar zenith angles for this timestep
    !
    ! !REVISION HISTORY: 
    !   26 Apr 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    real, dimension(im(region)*3,j1:j2) :: mu3
    real                                :: dt
    integer                             :: it, i, ii
    !
    mu3 = 0.0
    !
    ! split chemistry timestep in intervals and calculate mu at 1/6, 3/6
    ! and 5/6 of timestep chemistry. Interpolate in time using the
    ! longitudinal dependence of tha and take the average.
    !
    do ii=1,5,2
       dt = ( tstart + ii*deltat/6. ) * (real(im(region)*3)/24.)
       it = int(dt)   
       dt = dt-it
       mu3  = mu3 + (1-dt)*cshift(tha,it,1)+dt*cshift(tha,it+1,1)
    enddo
    !
    ! average over three angles at three times resolution
    !
    do i=i1,i2
       sza(i,:) = (mu3(i*3-2,:) + mu3(i*3-1,:) + mu3(i*3,:))/9.0
    enddo

  END SUBROUTINE GET_SZA
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    PHOTO_FLUX
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE PHOTO_FLUX( region, is, js, sza, rj_new)
    !********************************************************************** 
    !                                                                     * 
    !     contact:                                                        * 
    !                                                                     *
    !     Jason Williams                                                  * 
    !     Koninklijk Nederlands Meteorologisch instituut (KNMI)           *
    !     Postbus 201                                                     *
    !     3730 AE De Bilt                                                 *
    !     tel +31 (0)30 2206748                                           * 
    !     e-mail : williams@knmi.nl                                       *
    !                                                                     *
    !********************************************************************** 
    ! purpose:
    !
    ! The calculation of the radiation fields using input from TM5 and ECMWF meteo.
    !
    ! method: 
    !
    ! Optical properties (optical thickness,cross sections, quantum yields)
    ! are calculated for lm(region) layers. Radiation fluxes are derived for 0:lm(region)
    ! levels, including the surface and top of atmosphere. Thus, level l overlays layer l.
    ! Photolysis rates in a layer are evaluated by taking the average of the
    ! photolysis rates evaluated at the upper and lower boundaries of the layer.
    !
    ! First, spectral calculations are performed for an atmosphere with absorption
    ! only, including gaseous absorption by O2 and O3 and aerosol absorption.
    ! Second, a correction is applied for scattering and surface reflection, per
    ! scattering band and using radiative transfer calculations at representative
    ! wavelengths 
    !
    ! This photolysis module is based on:
    !
    !     Landgraf and Crutzen, 1998, J. Atmos. Sci, 55, 863-878.
    !     Williams et al., 2006, Atmos.Chem.Phys., 6, 4137-4161.   
    !
    !------------------------------------------------------------------
    !
    ! Albedo field
    ! ------------
    !
    ! In older TM photolysis code, an adhoc uv-vis albedo was computed
    ! from land cover fields:
    !
    !   module dry_deposition
    !     dd_surface_fields
    !       if(root_k)then
    !         dd_calc_inisurf
    !           compute ags on glb1x1 from land cover, on root_k only
    !         dd_coarsen_vd
    !           coarsen ags to vd_temp, copy into phot_dat(region)%albedo
    !       endif
    !       broadcast phot_dat(region)%albedo to all other processors
    !       
    ! In the first stratospheric codes, the same adhoc albedo computed in 
    ! drydepos was written over the albedo meteo field, and then this one was used.
    ! To avoid this obscure destruction of ecmwf data, we use the phot_dat version here.
    ! In future, this should be replaced by ECMWF field:
    !     UV visible albedo for direct radiation  ALUVP         (0 - 1)         15      128 
    !
    !------------------------------------------------------------------ 
    ! subroutine to calculate direct flux and actinic flux
    !*******************************************************************        
    !
    ! !USES:
    !
    use dims,      only : lm, idate, ndyn, ndyn_max
    use MeteoData, only : temper_dat, gph_dat, cc_dat
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in) :: region
    integer, intent(in) :: is, js
    real,    intent(in) :: sza(is:,js:)
    !
    ! !OUTPUT PARAMETERS:
    !
    !  photolysis rates for this time-step     
    real, dimension(is:,js:,:,:), intent(out) :: rj_new ! (i1:i2,j1:j2,lm(region),nj)
    !
    ! !REVISION HISTORY: 
    !   26 Apr 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    integer :: i, j, l, k, n, table_pos
    real    :: zangle, alb, esrm2
    real, dimension(lm(region)+1)   :: levels 
    logical                         :: clear, cloudy, debug_flag
    ! temperature index array for the lookup table
    integer, parameter              :: n_ind =34 ! number of increments in the look-up table
    integer, parameter              :: jindex=9
    integer, parameter              :: latindex=21
    real, parameter                 :: incr = 5. ! 5deg C intervals in look-up table (optimised value)
    real, dimension(n_ind)          :: temp_ind   

    real, dimension(:),    allocatable  :: v2_col, v3_col, dv2_col, dv3_col, t_lay, ly_a, column_cloud
    real, dimension(:,:),  allocatable  :: cst_o3_col
    real, dimension(:,:),  allocatable  :: cst_no2_col, qy_no2_col
    real, dimension(:,:),  allocatable  :: cst_hno3_col
    real, dimension(:,:),  allocatable  :: cst_h2o2_col
    real, dimension(:,:),  allocatable  :: cst_n2o5_col
    real, dimension(:,:),  allocatable  :: cst_ch2o_col 
    real, dimension(:,:),  allocatable  :: cst_pan_col
    real, dimension(:,:),  allocatable  :: cst_no3_col
    real, dimension(:,:),  allocatable  :: qy_o1d_col
    real, dimension(:,:),  allocatable  :: qy_ch3cocho_col
    real, dimension(:,:),  allocatable  :: qy_co_col,qy_c2o3_col     
    real, dimension(:,:),  allocatable  :: fdir, fact
    real, dimension(:),    allocatable  :: taua_cld_col, taus_cld_col
    real, dimension(:,:,:),allocatable  :: taus_aer_col, taua_aer_col
    real, dimension(:,:),  allocatable  :: ts_pi_clr_col, ta_pi_clr_col, g_pi_clr_col, g_pi_cld_col
    real, dimension(:),    allocatable  :: ts_pi_cld, ta_pi_cld 
    real, dimension(:,:,:),allocatable  :: pm_aer_col
    real, dimension(:),    allocatable  :: pm_cld_col
    real, dimension(:,:),  allocatable  :: ds
    real, dimension(:,:),  allocatable  :: rj_column   
    real, dimension(:,:,:),pointer      :: temperature
    ! additional diagnostic array to compare old/new jvalue input
    real, dimension(:), allocatable      :: vo3s_col
    real                                 :: zangle_in, rdif_szalims
    integer :: i1, j1, i2, j2

    call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )

    allocate( v2_col  (lm(region)+1) )
    allocate( v3_col  (lm(region)+1) )
    allocate( dv2_col (lm(region)  ) )
    allocate( dv3_col (lm(region)  ) )
    allocate( t_lay   (lm(region)  ) )
    allocate( cst_o3_col     (lm(region),maxwav) )
    allocate( cst_no2_col    (lm(region),89)     )
    allocate( qy_no2_col     (lm(region),89)     )
    allocate( qy_co_col      (lm(region),89)     )
    allocate( qy_c2o3_col    (lm(region),89)     )    
    allocate( cst_hno3_col   (lm(region),65)     )
    allocate( cst_h2o2_col   (lm(region),45)     )
    allocate( cst_n2o5_col   (lm(region),55)     )
    allocate( cst_ch2o_col   (lm(region),105)    )
    allocate( cst_pan_col    (lm(region),65)     )
    allocate( cst_no3_col    (lm(region),62)     )
    allocate( qy_o1d_col     (lm(region),65)     )
    allocate( qy_ch3cocho_col(lm(region),52)     )
    allocate( fdir           (lm(region),maxw)   ) ; fdir(:,:) = 0.
    allocate( taus_aer_col  (lm(region),nbands_trop,ngrid))
    allocate( taua_aer_col  (lm(region),nbands_trop,ngrid))
    allocate( taus_cld_col  (lm(region)))
    allocate( taua_cld_col  (lm(region)))
    allocate( ts_pi_clr_col (lm(region),nbands_trop))
    allocate( ta_pi_clr_col (lm(region),nbands_trop))
    allocate( g_pi_clr_col  (lm(region),nbands_trop))
    allocate( g_pi_cld_col  (lm(region),nbands_trop))
    allocate( pm_aer_col    (lm(region),nbands_trop,ngrid))
    allocate( pm_cld_col    (lm(region)))
    allocate( fact          (lm(region),nbands_trop)) ; fact(:,:) = 0.
    allocate( column_cloud  (lm(region)))
    allocate( ds            (lm(region)+1,lm(region)))
    allocate( rj_column     (lm(region),nj)) ; rj_column = 0.

    allocate(vo3s_col(lm(region)+1))

    temperature => temper_dat(region)%data

    clear = .false.
    cloudy = .true.
    debug_flag = .false.

    ! define temperature index array from the lookup table
    do i = 1, 34
       temp_ind(i) = 182.5 + (i-1) * 5.0
    enddo

    ! Inverse of SZA difference in limits,i.e. 1./(94-85)
    rdif_szalims = 1.0 / ( sza_widelimit - sza_limit )

    do j = j1, j2
       do i = i1, i2

          zangle = sza(i,j)
          ! Limit value used by photolysis computation to 'sza_limit' ...
          zangle_in = min( sza_limit, zangle )
          
          t_lay = temperature(i,j,:)

          rj_new(i,j,:,:) = 0.
          
          if ( zangle <= sza_widelimit ) then
             
             ! initialise
             fdir = 0. ; dv2_col = 0. ; dv3_col = 0. ; fact = 0.

             v2_col(:) = phot_dat(region)%v2(i,j,:)
             v3_col(:) = phot_dat(region)%v3(i,j,:)
             qy_ch3cocho_col(:,:) = phot_dat(region)%qy_ch3cocho(i,j,:,:)
             qy_c2o3_col(:,:) = phot_dat(region)%qy_c2o3(i,j,:,:)
             column_cloud(:) = cc_dat(region)%data(i,j,:)
             levels=gph_dat(region)%data(i,j,:)

             ! set optical depths parameters
             taus_cld_col =  phot_dat(region)%taus_cld(i,j,:) 
             taua_cld_col =  phot_dat(region)%taua_cld(i,j,:) 
             taus_aer_col =  phot_dat(region)%taus_aer(i,j,:,:,:) 
             taua_aer_col =  phot_dat(region)%taua_aer(i,j,:,:,:)
             pm_cld_col = phot_dat(region)%pmcld(i,j,:)  
             pm_aer_col = phot_dat(region)%pmaer(i,j,:,:,:) 

             !
             ! assign the cross-sections here to avoid the exporting of large arrays
             !   
             do k=1,lm(region)

                table_pos=int((t_lay(k) - (temp_ind(1)-0.5*5)) / 5) + 1
                table_pos=min(max(1,table_pos),34)  ! bug fix for temperatures below 175.5
                cst_o3_col(k,:) = xs_o3_look(:,table_pos)
                cst_no2_col(k,:) = xs_no2_look(:,table_pos)
                cst_hno3_col(k,:) = xs_hno3_look(:,table_pos)
                cst_h2o2_col(k,:) = xs_h2o2_look(:,table_pos)
                cst_ch2o_col(k,:) = xs_ch2o_look(:,table_pos)
                cst_n2o5_col(k,:) = xs_n2o5_look(:,table_pos)
                cst_pan_col(k,:) = xs_pan_look(:,table_pos)
                cst_no3_col(k,:) = xs_no3_look(:,table_pos)
                qy_o1d_col(k,:) = qy_o3_look(:,table_pos)
                qy_no2_col(k,:) = qy_no2_look(:,table_pos)
                qy_co_col(k,:) = qy_co_look(:,table_pos)                
             enddo

             ! slant column, lyman-alpha and direct flux ; all are zenith angle dependent

             call directflux(region,zangle_in,v2_col,v3_col,cst_o3_col,t_lay,fdir,dv2_col,dv3_col,taua_cld_col, &
                  taua_aer_col,levels,ds,vo3s_col,debug_flag)

             ! now do radiative transfer calculation, calculate actinic flux
             ! set albedo
             alb = phot_dat(region)%albedo(i,j)

             ! Now call the PIFM RT solver for the online calculation ; both clear and cloudy conditions can be used

             if  (clear) call pifm(region,zangle_in,alb,cst_o3_col,dv2_col,dv3_col, &
                  taua_aer_col,taus_aer_col,pm_aer_col,fact)
             if (cloudy) call pifm_ran( region,zangle_in,alb,cst_o3_col,dv2_col,dv3_col, &
                  taua_cld_col,taus_cld_col,pm_cld_col,taua_aer_col,taus_aer_col,pm_aer_col,fact,column_cloud)

             rj_column=0.

             call photorates_tropo(region,zangle_in,cst_o3_col,cst_no2_col,cst_hno3_col,cst_h2o2_col,cst_ch2o_col,&
                  cst_n2o5_col,cst_pan_col,cst_no3_col,qy_no2_col,qy_o1d_col,qy_ch3cocho_col,qy_co_col,qy_c2o3_col, &
                  fact,fdir,rj_column,t_lay,debug_flag) 

             if ( zangle <= sza_limit) then
                ! If actual SZA is below 85 deg then use values as such...
                rj_new(i,j,:,:) = rj_column
             else
                ! Otherwise, with SZA is between 85 deg and 95 then use interpolated values 
                rj_new(i,j,:,:) = rj_column * (sza_widelimit - zangle) * rdif_szalims
             end if

          elseif(zangle > sza_widelimit) then

             vo3s_col = 0.

             do k=1,lm(region)
                table_pos=int((t_lay(k) - (temp_ind(1)-0.5*5)) / 5) + 1  
                table_pos=min(max(1,table_pos),34)  ! max(...) == temperatures below 175.5 used like 180.5
                cst_o3_col(k,:) = xs_o3_look(:,table_pos)
             enddo

             debug_flag=.false.

          endif ! sza_limit

       enddo ! i
    enddo  ! j 


    ! Correction photolysis for distance sun-earth
    esrm2 = sundis(idate(2),idate(3))
    rj_new = rj_new *esrm2

    ! Store
    phot_dat(region)%jo3 (:,:,:) = rj_new(:,:,:,jo3d) 
    phot_dat(region)%jno2 (:,:,:) = rj_new(:,:,:,jno2)
    phot_dat(region)%jhno2 (:,:,:) = rj_new(:,:,:,jhono)
    phot_dat(region)%jch2oa (:,:,:) = rj_new(:,:,:,jach2o)
    phot_dat(region)%jch2ob (:,:,:) = rj_new(:,:,:,jbch2o)
    
    ! Averages
    phot_dat(region)%nj_av     = phot_dat(region)%nj_av   + float(ndyn)/float(ndyn_max)
    phot_dat(region)%nalb_av   = phot_dat(region)%nalb_av + float(ndyn)/float(ndyn_max)
    phot_dat(region)%jo3_av    = phot_dat(region)%jo3_av  + float(ndyn)/float(ndyn_max) * rj_new(:,:,:,jo3d) 
    phot_dat(region)%jno2_av   = phot_dat(region)%jno2_av + float(ndyn)/float(ndyn_max) * rj_new(:,:,:,jno2)
    phot_dat(region)%jhno2_av  = phot_dat(region)%jhno2_av + float(ndyn)/float(ndyn_max) * rj_new(:,:,:,jhono)
    phot_dat(region)%jch2oa_av = phot_dat(region)%jch2oa_av + float(ndyn)/float(ndyn_max) * rj_new(:,:,:,jach2o)
    phot_dat(region)%jch2ob_av = phot_dat(region)%jch2ob_av + float(ndyn)/float(ndyn_max) * rj_new(:,:,:,jbch2o)
      
    phot_dat(region)%ags_av  = phot_dat(region)%ags_av  + float(ndyn)/float(ndyn_max) * phot_dat(region)%albedo
    phot_dat(region)%sza_av  = phot_dat(region)%sza_av  + float(ndyn)/float(ndyn_max) * sza

    ! Done
    nullify(temperature)

    deallocate(v2_col)
    deallocate(v3_col)
    deallocate(dv2_col)
    deallocate(dv3_col)
    deallocate(cst_o3_col)
    deallocate(cst_no2_col)
    deallocate(qy_no2_col)
    deallocate(qy_co_col)
    deallocate(qy_c2o3_col)
    deallocate(cst_hno3_col)
    deallocate(cst_h2o2_col)
    deallocate(cst_n2o5_col)
    deallocate(cst_ch2o_col)
    deallocate(cst_pan_col)
    deallocate(cst_no3_col)
    deallocate(qy_o1d_col)
    deallocate(qy_ch3cocho_col)
    deallocate(fdir)
    deallocate(t_lay)
    deallocate(taus_aer_col)
    deallocate(taua_aer_col)
    deallocate(pm_aer_col)
    deallocate(taus_cld_col)
    deallocate(taua_cld_col)
    deallocate(pm_cld_col)
    deallocate(ts_pi_clr_col)
    deallocate(ta_pi_clr_col)
    deallocate(g_pi_clr_col)
    deallocate(g_pi_cld_col)

    deallocate(fact)
    deallocate(column_cloud)
    deallocate(ds)
    deallocate(rj_column)

    deallocate(vo3s_col)


  END SUBROUTINE PHOTO_FLUX
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    DIRECTFLUX
  !
  ! !DESCRIPTION: 
  !         Schumann-Runge parameterization
  !         Koppers GAA, Murtagh DP, Ann. Geophysicae 14, 68-79
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine directflux(region,zangle,v2_col,v3_col,cst_o3_col,t_lay,&
       fdir,dv2_col,dv3_col,taua_cld_col,taua_aer_col,&
       levels,ds,vo3s_col,debug_flag)
    !
    ! !USES:
    !
    use dims,  only : lm
    use binas, only : pi, avog, grav, xmair, ae
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)    :: region 
    real,intent(in)        :: zangle
    ! oxygen and ozone column
    real,dimension(lm(region)+1),intent(in)                   :: v2_col,v3_col
    real,dimension(lm(region)),intent(in)                     :: taua_cld_col
    real,dimension(lm(region),nbands_trop,ngrid),intent(in)   :: taua_aer_col
    !
    ! !OUTPUT PARAMETERS:
    !
    real,dimension(lm(region)+1),intent(out)                  :: vo3s_col
    ! temperature
    real,dimension(lm(region)),intent(in)                     :: t_lay
    ! temperature dependent cs of ozone
    real,dimension(lm(region),maxwav),intent(in)              :: cst_o3_col
    ! flux through pure absorbing atmosphere
    real,dimension(lm(region),maxw),intent(out)               :: fdir

    real,dimension(lm(region)),intent(out)                    :: dv2_col,dv3_col
    ! levels for ds calculation
    real,dimension(lm(region)+1)                              :: levels
    real,dimension(lm(region)+1,lm(region)),intent(out)       :: ds
    logical,intent(in)                                        :: debug_flag
    !
    ! !REVISION HISTORY: 
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    real,dimension(:),allocatable   :: t_lev  
    real,dimension(:),allocatable   :: v2s,v3s,dv2s, dv3s  
    integer                         :: l, li, i, j, k, k1, k2, id, n, h, bandno 
    real                            :: dl_o2,u0     
    real,dimension(20)              :: coeff_a, coeff_b 
    real                            :: a,b,c, gamma
    real,dimension(maxw)            :: fdir_top
    real                            :: fdir_bot
    real                            :: ta_o2, ta_o3, tau
    real                            :: sp, const, press, alp, alv2, sln, ck, sf  
    !-----------------------------------------------------------------     
    ! declarations for the calculation of DS
    real,dimension(0:lm(region))     :: ze
    real,dimension(0:lm(region),lm(region)) :: ds_tmp
    real                             :: rpsinz,rj,rjp1,diffj,height1,height2
    real                             :: zenrad,sm,re,dsj,dummy

    allocate(t_lev(lm(region)+1)) ; t_lev=0.
    allocate(dv2s(lm(region))) ; dv2s = 0.
    allocate(dv3s(lm(region))) ; dv3s = 0.
    allocate(v2s(1:lm(region)+1)) ; v2s = 0.
    allocate(v3s(1:lm(region)+1)) ; v3s = 0.

    ! initialise all output
    fdir = 0. ; dv2_col = 0 ; dv3_col = 0. ; ds = 0. ; ds_tmp = 0.

    ! define temperature on grid levels. Note temperature on layers now has reversed vertical grid
    do k  = 2,lm(region)
       t_lev(k) = (t_lay(k) + t_lay(k-1) ) * 0.5
    enddo

    ! boundaries
    t_lev(1) = t_lay(1) 
    t_lev(lm(region)+1) = t_lay(lm(region))

    a = 0.50572 ; b = 6.07995 ; c = 1.63640

    dv2_col(lm(region))  = v2_col(lm(region))      
    dv3_col(lm(region))  = v3_col(lm(region)) 

    do l=lm(region)-1,1,-1
       dv2_col(l)  = v2_col(l)-v2_col(l+1)
       dv3_col(l)  = v3_col(l)-v3_col(l+1) 
    end do

    !-----correction of the air mass factor---------------------------------                                                                       
    !     F. Kasten and T. Young, Revised optical air mass tabels and       
    !     approximation formula (1989) Aplied Optics Vol 28, No. 22 P. 4735 
    !     and J. Lenoble, atmospheric radiative transfer (1993), P. 236     
    !    define air mass factor in mu

    gamma = 90. - zangle
    u0 = sin(gamma*pi/180.) + a *(gamma+b)**(-c)        
    u0 = min(1.,u0)

    if(zangle <= 75.) then
       ds = 1/u0
    elseif(zangle >75. .and. zangle <=sza_limit) then
       levels(1:lm(region)+1)=levels(lm(region)+1:1:-1)

       re = (ae+levels(lm(region)+1))*1.e-3

       do k=lm(region)+1,1,-1
          ze(k-1) = (levels(k)/1.e3)-(levels(lm(region)+1)/1.e3)
       enddo

       zenrad=acos(u0)

       do k=0,lm(region)
          rpsinz = (re+ze(k))*sin(zenrad)

          id = k
          do n=1,id
             sm = 1.0
             rj = re + ze(n-1)
             rjp1 = re + ze(n)

             diffj = ze(n-1) - ze(n)

             height1 = rj**2 - rpsinz**2
             height2 = rjp1**2 - rpsinz**2

             height1=max(0.0,height1)
             height2=max(0.0,height2)

             if(id>k .and. n==id) then
                dsj=sqrt(height1)
             else
                dsj=sqrt(height1)-sm*sqrt(height2)
             endif

             ds_tmp(k,n) = dsj/diffj

          enddo ! n
       enddo !k
       ! invert matrix to be compatible with lm(region)+1 being tOA
       ds(1:lm(region)+1,1:lm(region))=ds_tmp(lm(region):0:-1,lm(region):1:-1)
    ENDIF ! sza_rad   

    ! slant column : calculate in different way depending on zangle
    if(zangle <=75.) then 
       v2s(lm(region)+1)=ds(lm(region),lm(region))*v2_col(lm(region)+1)  
       v3s(lm(region)+1)=ds(lm(region),lm(region))*v3_col(lm(region)+1) 
       do k1=lm(region)+1,1,-1 
          v2s(k1) = v2s(lm(region)+1)
          v3s(k1) = v3s(lm(region)+1)

          do k2=lm(region),k1,-1
             v2s(k1) = v2s(k1) + ds(k1,k2)*dv2_col(k2)    
             v3s(k1) = v3s(k1) + ds(k1,k2)*dv3_col(k2)    
          end do
       end do
    elseif(zangle>75. .and. zangle<=85.) then
       v2s(lm(region)+1) = ds(lm(region),lm(region))*v2_col(lm(region)+1)
       v3s(lm(region)+1) = ds(lm(region),lm(region))*v3_col(lm(region)+1)
       do k1=lm(region)+1,1,-1
          v2s(k1) = v2s(lm(region)+1)
          v3s(k1) = v3s(lm(region)+1)
          do k2=lm(region),1,-1
             v2s(k1) = v2s(k1) + ds(k1,k2)*dv2_col(k2)    
             v3s(k1) = v3s(k1) + ds(k1,k2)*dv3_col(k2)
          enddo
       enddo
    endif

    do k1=lm(region)+1,1,-1
       vo3s_col(k1)=v3s(k1)
    enddo

    ! differential slant column
    dv2s(lm(region)) = v2s(lm(region))-v2s(lm(region)+1)
    dv3s(lm(region)) = v3s(lm(region))-v3s(lm(region)+1)    
    do k = lm(region)-1,1,-1
       dv2s(k) = v2s(k)-v2s(k+1)     
       dv3s(k) = v3s(k)-v3s(k+1)
       if(dv2s(k)<0.) dv2s(k)=-1.0*dv2s(k) 
       if(dv3s(k)<0.) dv3s(k)=-1.0*dv3s(k)   
    enddo

    ! intialise direct flux
    fdir_top(1:maxw) = flux(1:maxw)     
    ! direct flux = actinic flux for a purely abs. atmosphere
    ! here, layer averaged quantity 

    bandno=1

    if(zangle <=85.) then
       do l = 1,maxw
          do k = lm(region),1,-1 
             if(l<18) then
                ta_o2 = cross_o2(l)*dv2s(k)
             else
                ta_o2 = 0.
             endif
             ta_o3 = cst_o3_col(k,l)*dv3s(k)
             fdir_bot  = fdir_top(l) * exp(-ta_o2-ta_o3-taua_cld_col(k)-taua_aer_col(k,bandno,1)) 
             fdir(k,l) = (fdir_top(l) + fdir_bot)/2.
             fdir_top(l) = fdir_bot
             if(l>lfin(bandno)) bandno=bandno+1 
          enddo
       enddo
    endif

    deallocate(t_lev)
    deallocate(v2s)
    deallocate(v3s)
    deallocate(dv2s)
    deallocate(dv3s)

  END SUBROUTINE DIRECTFLUX
  !EOC
  
#ifdef with_optics
  ! only define when coupling with photolysis is turned on
  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    AEROSOL_INFO_M7
  !
  ! !DESCRIPTION: assignment of aerosol optical depths calculated in M7  
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE AEROSOL_INFO_M7(region, status)
    !
    ! !USES:
    !
    USE TM5_DISTGRID, ONLY : dgrid, Get_DistGrid
    USE DIMS,         ONLY : isr,ier,jsr,jer, ndyn_max, ndyn
    USE METEODATA,    ONLY : gph_dat
    USE DIMS,         ONLY : im, jm, lm
    USE OPTICS,       ONLY : optics_aop_get
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)  :: region
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out) ::  status
    !
    ! !REVISION HISTORY: 
    !      Aug 2012 - NB
    !    5 Mar 2013 - Ph. Le Sager - TM5-MP version
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter  ::  rname = mname//'/aerosol_info_m7'

    integer                               :: i,j,l,lvec, nsend, igrid, i1, i2, j1, j2, lmr
    real, dimension(:,:,:),   pointer     :: gph
    real, dimension(:,:,:,:), allocatable :: taus_aer, taua_aer, pmaer

    ! Optics output fields (to be used and allocated by methods using the optics)
    real, dimension(:,:,:),   allocatable :: aop_out_ext ! extinctions
    real, dimension(:,:),     allocatable :: aop_out_a   ! single scattering albedo
    real, dimension(:,:),     allocatable :: aop_out_g   ! assymetry factor

    ! ------------ begin -----------------
    
    ! count lvec, tile grid size
    CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
    lmr=lm(region)
    lvec = (i2-i1+1)*(j2-j1+1)*lmr

    allocate( aop_out_ext(lvec, nwdep, 1))
    allocate( aop_out_a  (lvec, nwdep)   )
    allocate( aop_out_g  (lvec, nwdep)   )

    CALL OPTICS_AOP_GET(lvec, region, nwdep, wdep, 1, .false., &
                        aop_out_ext, aop_out_a, aop_out_g, status)
    IF_NOTOK_RETURN(status=1)

    gph => gph_dat(region)%data

    allocate(taus_aer (i1:i2, j1:j2, lmr,nwdep)); taus_aer = 0.0
    allocate(taua_aer (i1:i2, j1:j2, lmr,nwdep)); taua_aer = 0.0  
    allocate(pmaer    (i1:i2, j1:j2, lmr,nwdep)); pmaer    = 0.0

    ! ---------------------------------
    ! unpack results from calculate_aop   
    ! ---------------------------------
    lvec = 0
    do l = 1, lmr
    do j = j1,j2
    do i = i1,i2

       lvec = lvec + 1

       taus_aer(i,j,l,1:nwdep) = aop_out_ext(lvec,1:nwdep,1) * aop_out_a(lvec,1:nwdep) * (gph(i,j,l+1) - gph(i,j,l))       ! from 1/m to nondim
       taua_aer(i,j,l,1:nwdep) = aop_out_ext(lvec,1:nwdep,1) * (1.-aop_out_a(lvec,1:nwdep)) * (gph(i,j,l+1) - gph(i,j,l))
       pmaer   (i,j,l,1:nwdep) = aop_out_g  (lvec,1:nwdep)
       
    enddo
    enddo
    enddo

    nullify(gph)

    do i=1, nbands_trop
       phot_dat(region)%taus_aer(:,:,:,i,1) = taus_aer(:,:,:,wav_grid(i))
       phot_dat(region)%taua_aer(:,:,:,i,1) = taua_aer(:,:,:,wav_grid(i))
       phot_dat(region)%pmaer   (:,:,:,i,1) = pmaer   (:,:,:,wav_grid(i))

       phot_dat(region)%taus_aer(:,:,:,i,2) = taus_aer(:,:,:,wav_gridA(i))
       phot_dat(region)%taua_aer(:,:,:,i,2) = taua_aer(:,:,:,wav_gridA(i))
       phot_dat(region)%pmaer   (:,:,:,i,2) = pmaer   (:,:,:,wav_gridA(i))
    enddo

    phot_dat(region)%naer_av     = phot_dat(region)%naer_av     + float(ndyn)/float(ndyn_max)
    phot_dat(region)%pmaer_av    = phot_dat(region)%pmaer_av    + float(ndyn)/float(ndyn_max) * phot_dat(region)%pmaer 
    phot_dat(region)%taus_aer_av = phot_dat(region)%taus_aer_av + float(ndyn)/float(ndyn_max) * phot_dat(region)%taus_aer 
    phot_dat(region)%taua_aer_av = phot_dat(region)%taua_aer_av + float(ndyn)/float(ndyn_max) * phot_dat(region)%taua_aer 

    deallocate(taus_aer)
    deallocate(taua_aer)
    deallocate(pmaer)

    Deallocate(aop_out_ext)
    Deallocate(aop_out_a  )
    Deallocate(aop_out_g  )

    status = 0

  END SUBROUTINE AEROSOL_INFO_M7
  !EOC
#endif

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    AEROSOL_INFO
  !
  ! !DESCRIPTION: assignment of aerosol optical depths
  ! 
  ! This is a crude method to provide some average values for the absorption and scattering
  ! terms by background aerosol choice of values defined according to the relative humidity.
  ! Absorption component can be set to zero throughout (M.van Weele, private comm.,2005).
  ! Scattering component chosen to be representative of background aerosol
  ! Moreover there is a choice as the whether the values defined by shettle and fenn are used
  ! This will require a look up table on 1x1,60 layers with indexes 1->4 with respect to
  ! location. At the moment background aerosol is chosen throughout
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE AEROSOL_INFO( region )
    !
    ! !USES:
    !
    use dims,         only : im, jm, lm, newsrun, ndyn, ndyn_max
    use binas,        only : xm_air, xm_h2o,grav
    use global_data,  only : mass_dat
    use MeteoData,    only : temper_dat, humid_dat, gph_dat, cc_dat, iwc_dat, lwc_dat, phlb_dat
    use MeteoData,    only : lsmask_dat
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)       :: region
    !
    ! !REVISION HISTORY: 
    !      Jun 2005 - JEW - 
    !   26 Apr 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    real, dimension(:,:,:,:,:), allocatable  :: taus_aer, taua_aer
    real, dimension(:,:,:,:,:), allocatable  :: pmaer
    real, dimension(:,:,:),     allocatable  :: rhum, dens, part, dz, press_lay
    real, dimension(:,:,:),     pointer      :: q, phlb, temp, gph, frac, xland
    real                                     :: a_sc, b_sc, a_ab, b_ab, a_g, b_g
    real                                     :: bsa, baa, ga
    integer                                  :: i_type
    integer, dimension(:,:),    pointer      :: aero_clim

    real      :: wv, tr, scale_aero, lay1, lay2
    integer   :: i, j, l, k, kboun, i_ref, b, wav, ll, n 
    integer   :: i1, j1, i2, j2, lmr
    logical   :: aerosol, shettle_and_fenn

    real, dimension(8) :: rh_ref = (/0., 50., 70., 80., 90., 95., 98., 99./)

    real, dimension(4) :: pn_ref = (/ 15000., 4000., 20000., 5000./)
    !
    !     different aerosol types:                                          
    !        1 = rural aerosol                                              
    !        2 = maritime aerosol                                           
    !        3 = urban aerosol                                              
    !        4 = free troposphere aerosol               
    !-----------------------------------------------------------------------

    aerosol = .false.
    shettle_and_fenn = .true.

    call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
    lmr = lm(region)
    
    allocate(taus_aer(i1:i2, j1:j2, lmr, nbands_trop, ngrid))
    allocate(taua_aer(i1:i2, j1:j2, lmr, nbands_trop, ngrid))  
    allocate(pmaer   (i1:i2, j1:j2, lmr, nbands_trop, ngrid))
    allocate(rhum    (i1:i2, j1:j2, lmr))
    allocate(dens    (i1:i2, j1:j2, lmr))
    allocate(part    (i1:i2, j1:j2, lmr))
    allocate(dz      (i1:i2, j1:j2, lmr))
    allocate(press_lay(i1:i2, j1:j2, lmr))

    ! Initialize values

    taus_aer=0.
    taua_aer=0.
    pmaer   =0.
    
    ! read in met. data
    temp      => temper_dat(region)%data  ! (i1:i2, j1:j2, 1:lmr)
    q         =>  humid_dat(region)%data  ! (i1:i2, j1:j2, 1:lmr)
    phlb      =>   phlb_dat(region)%data  ! 
    gph       =>    gph_dat(region)%data  ! (i1:i2, j1:j2, 1:lmr+1)
    frac      =>     cc_dat(region)%data  ! (i1:i2, j1:j2, 1:lmr)
    xland     => lsmask_dat(region)%data  ! (i1:i2, j1:j2, 1)
    
    aero_clim =>   phot_dat(region)%aero_surface_clim  ! (i1:i2, j1:j2)


    ! initialize values

    scale_aero = 1.0

    ! the land-sea mask is used to ascribe either marine or rural aerosol for the bottom layers

    dz = 0. ; dens = 0. ; rhum = 0. 

    phot_dat(region)%taus_aer = 0.
    phot_dat(region)%taua_aer = 0.
    phot_dat(region)%pmaer = 0. 

    if (shettle_and_fenn) then

       if(newsrun) then
          do j=j1,j2 
             do i=i1,i2
                if(j<=10 .or. j>=80) then
                   aero_clim(i,j)=4
                else    
                   if(xland(i,j,1) <  30.) aero_clim(i,j)=2
                   if(xland(i,j,1) >= 30.) aero_clim(i,j)=1
                endif
             enddo
          enddo
       endif

       do l=1,lm(region)
          do j=j1,j2 
             do i=i1,i2
                dz(i,j,l) = (gph(i,j,l+1) - gph(i,j,l))
                press_lay(i,j,l) =  (phlb(i,j,l) + phlb(i,j,l+1)) * 0.01 * 0.5 ! hPa
                dens(i,j,l) = 7.2427e18 * press_lay(i,j,l)/temp(i,j,l)
                tr=1.-373.15/temp(i,j,l)
                wv=exp((((-.1299*tr-.6445)*tr-1.976)*tr+13.3185)*tr)  
                if(frac(i,j,l)<0.95) rhum(i,j,l)=q(i,j,l)*dens(i,j,l)*xm_air/xm_h2o*temp(i,j,l)/(1013.25*wv*7.24e16)
                ! JEW assume saturation when cloud fraction is very high.
                if(frac(i,j,l)>=0.95)  rhum(i,j,l)=98.0
                rhum(i,j,l)=min(98.,rhum(i,j,l))
             enddo
          enddo
       enddo

       !
       ! use land-sea mask to ascribe either marine or rural aerosol for different grid cells
       ! in the lower few KM
       !
       do l=1,lm(region)
          do j=j1,j2 
             do i=i1,i2
                if(l<5) then
                   i_type=aero_clim(i,j)
                elseif(l>=5) then
                   i_type=4
                endif
                kboun = 1
                lay1 = (phlb(i,j,l)/phlb(i,j,kboun))**3
                lay2 = (phlb(i,j,l+1)/phlb(i,j,kboun))**3
                if(i_type<4) then
                   part(i,j,l) = scale_aero*pn_ref(i_type)*dz(i,j,l)*100.
                   kboun = l
                else
                   part(i,j,l) = scale_aero*pn_ref(i_type)*0.5*(lay1+lay2)*dz(i,j,l)*100.
                endif

                i_ref = 8

                do k = 1,8
                   if(rh_ref(k) < rhum(i,j,l)) i_ref = k
                enddo

                do n=1,ngrid
                   do b=1,nbands_trop

                      if (n==1) wav=lmid(b)
                      if (n==2) wav=lmid_gridA(b) 

                      a_sc = (sca(wav,i_ref,i_type)-sca(wav,i_ref+1,i_type))/(rh_ref(i_ref)-rh_ref(i_ref+1))

                      b_sc = sca(wav,i_ref,i_type)- a_sc*rh_ref(i_ref)

                      a_ab = (abs_eff(wav,i_ref,i_type)-abs_eff(wav,i_ref+1,i_type))/(rh_ref(i_ref)-rh_ref(i_ref+1))

                      b_ab = abs_eff(wav,i_ref,i_type)- a_ab*rh_ref(i_ref)

                      a_g =(g(wav,i_ref,i_type)-g(wav,i_ref+1,i_type))/(rh_ref(i_ref)-rh_ref(i_ref+1))

                      b_g =g(wav,i_ref,i_type) - a_g*rh_ref(i_ref)           

                      bsa = a_sc*rhum(i,j,l) + b_sc
                      baa = a_ab*rhum(i,j,l) + b_ab
                      ga = a_g*rhum(i,j,l) + b_g

                      taus_aer(i,j,l,b,n) = bsa*part(i,j,l)
                      taua_aer(i,j,l,b,n) = baa*part(i,j,l)

                      if(taus_aer(i,j,l,n,1)>0.) then
                         ! do k = 1,nmom
                         pmaer(i,j,l,b,n)=ga
                         ! enddo
                      else
                         pmaer(i,j,l,b,n) = 0.
                      endif

                   enddo !nbands_trop
                enddo !ngrid

             enddo ! l
          enddo ! j
       enddo ! i

       phot_dat(region)%taus_aer = taus_aer
       phot_dat(region)%taua_aer = taua_aer
       phot_dat(region)%pmaer = pmaer 

    endif ! shettle and fenn switch

    ! JEW switch for the aerosol absorption and scattering properties  
    if (aerosol) then

       do l=1,lm(region) 
          do j=j1,j2 
             do i=i1,i2
                dens(i,j,l) = 7.2427e18 * press_lay(i,j,l)/temp(i,j,l)
                tr=1.-373.15/temp(i,j,l)
                wv=exp((((-.1299*tr-.6445)*tr-1.976)*tr+13.3185)*tr)  
                rhum(i,j,l)=q(i,j,l)*dens(i,j,l)*xm_air/xm_h2o*temp(i,j,l)/(1013.25*wv*7.24e16)
                if(rhum(i,j,l)>40. .and. rhum(i,j,l)<80.) pmaer(i,j,l,:,:) = 0.65
                if(rhum(i,j,l)>=80.)  pmaer(i,j,l,:,:) = 0.85
             enddo
          enddo
       enddo

       phot_dat(region)%taus_aer = 3.0e-3
       phot_dat(region)%taua_aer = 1.5e-4
       phot_dat(region)%pmaer = 0.7

    endif

    ! Averages
    phot_dat(region)%naer_av     = phot_dat(region)%naer_av     + (float(ndyn)/float(ndyn_max))
    phot_dat(region)%pmaer_av    = phot_dat(region)%pmaer_av    + (float(ndyn)/float(ndyn_max)) * phot_dat(region)%pmaer 
    phot_dat(region)%taus_aer_av = phot_dat(region)%taus_aer_av + (float(ndyn)/float(ndyn_max)) * phot_dat(region)%taus_aer 
    phot_dat(region)%taua_aer_av = phot_dat(region)%taua_aer_av + (float(ndyn)/float(ndyn_max)) * phot_dat(region)%taua_aer 
    
    ! Done
    deallocate(taus_aer)
    deallocate(taua_aer)
    deallocate(pmaer)
    deallocate(rhum)
    deallocate(dens)
    deallocate(part)
    deallocate(dz)
    deallocate(press_lay)

    nullify(q)
    nullify(temp)
    nullify(phlb)
    nullify(gph)
    nullify(frac)

    nullify(xland)
    nullify(aero_clim)

  END SUBROUTINE AEROSOL_INFO
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    SLINGO
  !
  ! !DESCRIPTION: A. Slingo's data for cloud particle radiative properties
  !               (from 'A GCM Parameterization for the Shortwave Properties of Water Clouds'
  !               JAS vol. 46 may 1989 pp 1419-1427)
  !
  !               Called everytime the meteo data is updated to calculate new
  !               cloud optical properties
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE SLINGO( region )                                                       
    !
    ! !USES:
    !
    use binas,         only : xmair, Avog
    use dims,          only : im, jm, lm, lmax_conv
    use global_data,   only : mass_dat
    use MeteoData,     only : gph_dat, lwc_dat, iwc_dat, cc_dat, phlb_dat, temper_dat, lsmask_dat
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in) :: region
    !
    ! !REVISION HISTORY: 
    !   26 Apr 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    real                               :: eff_rad ! effective radius no longer fixed but linked to LWP
    real                               :: rhodz, tau_c, xsa, press, airn, D_ge
    real                               :: rclf    ! inverse cloud fraction
    real                               :: clwc, ciwc ! [g/m3]
    real                               :: lfrac, sfrac ! scaled components due to land and sea fraction
    real,dimension(:,:,:),allocatable  :: taua_cld
    real,dimension(:,:,:),allocatable  :: taus_cld
    !total scattering optical depth for cloud layer (liquid+cirrus)
    !total absorption optical depth for cloud layer (liquid_cirrus)
    real,dimension(:,:,:),allocatable  :: pmcld !phase function (HG)

    !----------------------------Local------------------------------------------ 
    real,dimension(:,:,:),allocatable :: lwp !cloud liquid water path [g/m^2]
    real,dimension(:,:,:),pointer     :: frac, lwc, iwc, gph, phlb, temp, xland
        
    real, parameter  :: factor = 7.24e16*1.e6*xmair*29.2605/avog

    real,dimension(:,:,:), allocatable :: dz, totalwc, global_cloud_reff

    real, dimension(4)    :: abar, bbar, cbar, dbar, ebar, fbar
    ! A coefficient for extinction optical depth        
    ! B coefficiant for extinction optical depth        
    ! C coefficiant for single particle scat albedo     
    ! D coefficiant for single particle scat albedo     
    ! E coefficiant for asymmetry parameter             
    ! F coefficiant for asymmetry parameter             

    data abar/ 2.817e-02, 2.682e-02,2.264e-02,1.281e-02/ 
    data bbar/ 1.305    , 1.346    ,1.454    ,1.641    / 
    data cbar/-5.62e-08 ,-6.94e-06 ,4.64e-04 ,0.201    / 
    data dbar/ 1.63e-07 , 2.35e-05 ,1.24e-03 ,7.56e-03 / 
    data ebar/ 0.829    , 0.794    ,0.754    ,0.826    / 
    data fbar/ 2.482e-03, 4.226e-03,6.560e-03,4.353e-03/ 

    real   :: abari, bbari, cbari, dbari, ebari, fbari                                                        
    ! A coefficiant for current spectral interval       
    ! B coefficiant for current spectral interval       
    ! C coefficiant for current spectral interval       
    ! D coefficiant for current spectral interval       
    ! E coefficiant for current spectral interval       
    ! F coefficiant for current spectral interval       

    real    :: gc, wc, tot, tmp1 ,tmp2, tmp3
    !asymmetry factor                                    
    !single scattering albedo                            
    !total optical depth of cloud layer

    ! constants for scattering parameterization of Fu (1996)
    real, dimension(7) :: a0, a1

    data a0/-.236447E-03,-.236447E-03,-.266955E-03,-.266955E-03,-.266955E-03,-.258858E-03,.982244E-04/              
    data a1/.253817E+01 ,.253817E+01 ,.254719E+01 ,.254719E+01 ,.254719E+01 ,.253815E+01,.250875E+01/  

    ! Parameters for computing cloud effective radius according to Martin et al. (JAS 1994)
    real            :: denom,numerator
    real, parameter :: d_land=0.43
    real, parameter :: d_sea=0.33

    integer :: i,j,l,k, indxsl, n, i1, i2, j1, j2, lmr

    !----------- START----------

    call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
    lmr = lm(region)
    
    !-------------------------------------------------------------------------           
    !     Set index for cloud particle properties based on the wavelength,  
    !     according to A. Slingo (1989) equations 1-3:                      
    !     Use index 1 (0.25 to 0.69 micrometers) for visible                
    !     Use index 2 (0.69 - 1.19 micrometers) for near-infrared           
    !     Use index 3 (1.19 to 2.38 micrometers) for near-infrared          
    !     Use index 4 (2.38 to 4.00 micrometers) for near-infrared          

    indxsl = 1 

    !     Set cloud extinction optical depth, single scatter albedo,        
    !     asymmetry parameter, and forward scattered fraction:              

    abari = abar(indxsl)
    bbari = bbar(indxsl)
    cbari = cbar(indxsl)
    dbari = dbar(indxsl)
    ebari = ebar(indxsl)
    fbari = fbar(indxsl)
    !
    !--------------------------------------------------------------------------
    !     In the parameterization of Fu(1996) wavelength dependant extinction maybe
    !     calculated using a set of pre-defined constants:
    !     
    !      nm               a0                a1
    !    250-300         -.236447E-03      .253817E+01
    !    300-330         -.266955E-03      .254719E+01
    !    330-360         -.293599E-03      .254540E+01
    !    360-400         -.258858E-03      .253815E+01
    !    400-440         -.106451E-03      .252684E+01
    !    440-480          .129121E-03      .250410E+01
    !    480-520         -.954458E-04      .252061E+01
    !    520-570         -.303108E-04      .251805E+01
    !    570-640          .982244E-04      .250875E+01     
    !
    !    used in Beta=IWC(a0+a1/Dg)
    !--------------------------------------------------------------------------     
    allocate(taua_cld(i1:i2, j1:j2, lm(region)))
    allocate(taus_cld(i1:i2, j1:j2, lm(region)))
    allocate(pmcld   (i1:i2, j1:j2, lm(region))) 
    allocate(lwp     (i1:i2, j1:j2, lm(region)))

    allocate(dz               (i1:i2, j1:j2, lm(region)  ))
    allocate(totalwc          (i1:i2, j1:j2, lm(region)  ))
    allocate(global_cloud_reff(i1:i2, j1:j2, lm(region)  ))

    ! Initialize values

    taua_cld = 0. ; taus_cld = 0. ; pmcld = 0. ; lwp = 0.

    ! assume a baseline cloud radius, required for heterogeneous chemistry
    global_cloud_reff = 6.

    !JEW The ice and water particles are now treated seperately. For cloud the values are taken from slingo 
    !JEW the refractive index for ICE is very low below 750nm therefore T~T(scatt).

    iwc   => iwc_dat   (region)%data  ! (i1:i2, j1:j2, 1:lmr)
    frac  => cc_dat    (region)%data  ! (i1:i2, j1:j2, 1:lmr)
    gph   => gph_dat   (region)%data  ! (i1:i2, j1:j2, 1:lmr+1)
    lwc   => lwc_dat   (region)%data  ! (i1:i2, j1:j2, 1:lmr)
    phlb  => phlb_dat  (region)%data  ! 
    temp  => temper_dat(region)%data  ! (i1:i2, j1:j2, 1:lmr)
    xland => lsmask_dat(region)%data  ! (i1:i2, j1:j2, 1)

    ! No negative input
    lwc=max(0.,lwc)
    iwc=max(0.,iwc)  

    ! read in new cloud data to feed into the slingo routine. The values of lwc are zero in top levels so limit layer loop
   
    do l=1,lmax_conv
       do j=j1,j2 
          do i=i1,i2
             press = (phlb(i,j,l) + phlb(i,j,l+1)) * 0.5 ! Pa
             dz(i,j,l) = ( gph(i,j,l+1) - gph(i,j,l) )

             ! We have replaced the parameterization of McFarlane et al. (1992) with that of Martin et al. (1994).
             ! In the IFS a crude filter of 0.5 for the land fraction is used with the CNN(Marine) = 40.
             ! and CNN(Land) = 900. Here we weight the final CCN with the actual land fraction so as to introduce
             ! more variability. 

             ! JEW : there is a potential problem in that nonzero cloud fractions may occur with no associated lwc value
             ! on TM5 vertical resolutions, therefore a filter w.r.t. both fraction and cloud liquid path are included. 
             if( lwc(i,j,l) > 1.e-10 .and. frac(i,j,l) > 0.02 ) then
                ! JEW calculate total water path locally : convert to g/m(2) for slingo input
                ! following the conversion procedure for LWC from ECMWF input from old cloud subroutine.

                ! Included comments from old code to explain the prefactor in rhodz:
                ! rho  = airn*1e6*xmair/avo                   ! g/m3
                ! dz   = 29.2605*temp(i,j,l)* alog(phlb(i,j,l)/(1.0+phlb(i,j,l+1)))
                rclf=1./frac(i,j,l)
                rhodz = factor*press*alog(phlb(i,j,l)/(1.0+phlb(i,j,l+1)))

                ! VH - scale lwc and lwp with cloud fraction, to compute cloud optical properties
                ! and droplet radius representative for cloudy part only
                lwp(i,j,l) = rhodz*lwc(i,j,l)*rclf

                airn=7.24e16*press/temp(i,j,l) !#/cm3
                clwc=lwc(i,j,l)*xmair*airn*1.e6/avog ! g/m3

                ! Effective radius: Martin et al. (JAS 1994) parametrization
                !if (xland(i,j,l) > 50.) then
                ! Above land use ccn of ~900
                ! D_land = 0.43
                ! CCNLND = 900
                ! NTOT=-2.10E-04*CCNLAND*CCNLAND+0.568*CCNLAND-27.9
                ! DENOM=4.0*pi*rho_w*NTOT*(1.0+D_land*D_land)**3 
                ! with rho_w in g/m3.
                denom = 6547.52*1.e6
                numerator=3.0*clwc*rclf*(1.0+3.0*D_land*D_land)**2
                lfrac=0.
                if((numerator/denom) > 1.e-20) lfrac=1.e4*exp(0.333*LOG(numerator/denom))
                !else
                ! Maritime air mass
                ! D_sea = 0.333
                ! CCNSEA = 40
                ! NTOT=-1.15E-03*CCNSEA*CCNSEA+0.963*CCNSEA+5.30
                ! DENOM=4.0*pi*rho_w*NTOT*(1.0+D_sea*D_sea)**3
                denom = 723.210*1.e6
                numerator=3.0*clwc*rclf*(1.0+3.0*D_sea*D_sea)**2
                sfrac=0.
                if((numerator/denom) > 1.e-20) sfrac=1.e4*exp(0.333*LOG(numerator/denom))
                !endif

                ! TvN: Is this the best way to average land and sea?
                eff_rad=(xland(i,j,1)/100.*lfrac)+(1.0-xland(i,j,1)/100.)*sfrac

                ! prevent the radius of non-precipitation clouds being too big
                ! these size limits are equal to those chosen in the IFS
                eff_rad=min(16.0,max(eff_rad, 4.0))

                tmp1 = abari + bbari/eff_rad
                tmp2 = 1. - cbari - dbari*eff_rad
                tmp3 = fbari*eff_rad

                !     Do not let single scatter albedo be 1; delta-eddington solution
                !     for non-conservative case:
                wc   = min(tmp2,.999999)               
                tot  = lwp(i,j,l)*tmp1      
                gc   = ebari + tmp3 
                ! JEW : no wavelength dependence for the absorption or scattering effects due to liquid clouds !!!!!
                taua_cld(i,j,l) = (1.-wc)*tot 
                taus_cld(i,j,l) = wc*tot
                ! avoid possible negatives due to input data            
                taua_cld(i,j,l)=max(0.0,taua_cld(i,j,l))

                global_cloud_reff(i,j,l) = eff_rad

                ! JEW : for calculating the scattering component due to ice 
                if(iwc(i,j,l) > 1.e-10) then
                   airn=7.24e16*press/temp(i,j,l)
                   ! iwc in g/m3 from TM5 definition       
                   ciwc=iwc(i,j,l)*xmair*airn*1.e6/avog
                   ! Following Eq. (5) of Heymsfield and McFarquhar (JAS, 1996), 
                   ! the cross-sectional surface area A (km^-1) 
                   ! can be approximated by:
                   ! 10*IWC^0.9, where IWC in g/m3.
                   ! Thus, A in m^2/m^3 = m^-1 is given by 0.01 IWC^0.9,
                   ! which implies that xsa is in cm^-1.
                   xsa=1.0e-4*ciwc**0.9
                   !
                   ! calculate D_ge using the relationship in Fu (1996) where Beta=extinction co-efficient (m-1)
                   !
                   !  D_ge = 2(3)**0.5/(3 Rho)*(IWC/xsa)
                   ! 
                   ! The above relationship is for instance given in Table 1
                   ! in McFarquhar and Heymsfield (JAS, 1997).
                   ! Below an ice density in units g/cm^3 and
                   ! IWC in g/m^3 are used.
                   ! The conversion of IWC from g/m^3 to g/cm^3
                   ! gives a factor 10^6, which together with
                   ! the division by 100 converts from cm to um.
                   !
                   D_ge=(2.*1.73205/(3.*0.917))*(ciwc/xsa)
                   ! convert to uM
                   D_ge=D_ge/100.
                   !
                   ! Cirrus scattering has a wavelength dependancy
                   !
                   taus_cld(i,j,l)=taus_cld(i,j,l)+((ciwc*(a0(3)+(a1(3)/D_ge)))*dz(i,j,l))        
                endif

                if (taus_cld(i,j,l) > 0.) then 
                   pmcld(i,j,l)=gc 
                else 
                   pmcld(i,j,l)=0. 
                end if
             end if
          enddo
       enddo
    enddo

    ! Store optical depths and phase functions
    phot_dat(region)%pmcld      = pmcld
    phot_dat(region)%taus_cld   = taus_cld  
    phot_dat(region)%taua_cld   = taua_cld
    phot_dat(region)%cloud_reff = global_cloud_reff

    ! Averages
    phot_dat(region)%lwp_av      = lwp 
    phot_dat(region)%cfr_av      = phot_dat(region)%cfr_av      + frac
    phot_dat(region)%reff_av     = global_cloud_reff
    phot_dat(region)%ncloud_av   = phot_dat(region)%ncloud_av   + 1
    phot_dat(region)%pmcld_av    = phot_dat(region)%pmcld_av    + pmcld 
    phot_dat(region)%taus_cld_av = phot_dat(region)%taus_cld_av + taus_cld 
    phot_dat(region)%taua_cld_av = phot_dat(region)%taua_cld_av + taua_cld 

    nullify(lwc)
    nullify(frac)
    nullify(iwc)
    nullify(gph)
    nullify(phlb)
    nullify(temp,xland)

    deallocate(taua_cld)
    deallocate(taus_cld)
    deallocate(pmcld)
    deallocate(lwp)

    deallocate(dz)
    deallocate(totalwc)
    deallocate(global_cloud_reff)

  END SUBROUTINE SLINGO
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    UPDATE_CSQY
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE UPDATE_CSQY( region )
    !
    ! !USES:
    !
    use dims
    use binas,           only: Avog, xmair, grav
    use global_data,     only: mass_dat
    use MeteoData,       only: phlb_dat, temper_dat, m_dat
    use chem_param,      only: xmo3, fscale, iacet
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)                    :: region   
    !
    ! !REVISION HISTORY: 
    !   26 Apr 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    integer                                :: i,j,l 
    real,dimension(:,:,:),pointer          :: temperature   
    real,dimension(:,:,:),pointer          :: phlb
    real,dimension(:,:,:),pointer          :: mass
    real,dimension(:,:,:,:),pointer        :: qyacet    
    real,dimension(:,:,:),allocatable      :: density
    real,dimension(:,:,:),allocatable      :: pressure

    integer                      :: k, m, table_pos, i1, j1, i2, j2, lmr
    real                         :: xa, xb, te, chix, ww, temper, ptorr
    real                         :: phi,  qy, tt, alpha, delta_t
    real,dimension(maxwav)       :: rd, phi0
    real                         :: a0, b0, a1, b1, a2, b2, a3, b3, c3, a4, b4
    real                         :: term,term1,term2,t_scale            
    real                         :: airn,aird 
    
    integer,parameter            :: n_ind =34 ! number of increments in the look-up table
    real,dimension(n_ind)        :: temp_ind       
    
    phlb        => phlb_dat(region)%data
    temperature => temper_dat(region)%data
    qyacet      => phot_dat(region)%qy_c2o3
    mass        => m_dat(region)%data
    
    call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
    lmr = lm(region)
    
    allocate(density (i1:i2, j1:j2, lmr))
    allocate(pressure(i1:i2, j1:j2, lmr))

    !---------------------------------------------------------------------------------   
    ! Only update the temperature sensitive regions of the spectral
    ! grids. That is, initialize values values to either zero or cs_ values for selected
    ! species.                                  
    !--------------------------------------------------------------------------------- 

    do k  = 1,lm(region)
       do j  = j1,j2
          do i  = i1,i2
             ! define air pressure
             pressure(i,j,k) = (phlb(i,j,k) + phlb(i,j,k+1))/2
          enddo
       enddo
    enddo

    ! define air density
    density = 7.24e16*pressure/temperature(i1:i2, j1:j2, 1:lmr) 

    ! define temperature index array from the lookup table
    do i = 1, 34
       temp_ind(i) = 182.5 + (i-1) * 5.0
    enddo


    !---------------------------------------------------------------------
    !        QUANTUM YIELD OF METHYGLYOXAL S. KOCH and G. K. MOORTGAT       
    !        J Phys Chem, 102, 9142-9153, 1998.
    !
    ! JEW: the qy_ch3cocho array has been reduced to 52 bins to represent the pressure.
    !      dep. spectral bins ONLY !!   
    !---------------------------------------------------------------------

    ! Calculate co-efficients for METHYGLYOXAL qy outside loop
    phi0 = 1. ; rd = 0.      
    ! for wave < 380nm qy is essentially = 1.0.

    do l = 1,39
       if(l<32) then
          phi0(l) = 1
       elseif(l>=32) then
          phi0(l) = 8.15e-9 * exp(7131./wave(l+37)*1.e-7)
       endif
       rd(l)   = 7.34e-9 * exp(8793./wave(l+37)*1.e-7)
    enddo

    do l = 40,52
       phi0(l) = 3.63e-7 * exp(5693./wave(l+37)*1.e-7) 
       rd(l)   = 1.93e4*exp(-5639./wave(l+37)*1.e-7)
    enddo

    do k  = 1,lm(region) 
       do j = j1,j2
          do i = i1,i2

             !        QUANTUM YIELD OF METHYGLYOXAL JPL 2006 P4-71       
             !        CH3COCHO -> CH3C(O)O2 + HO2 + CO     (1)
             !                 -> CH4 + 2CO                (2)
             !                 -> CH3CHO + CO              (3) 
             ! the predominant products are from channel(1) for wav 240-480nm                                    

             ! Pa -> Torr          
             !ptorr= (pressure(i,j,k)/100.)*760./1.013*1e-3
             ptorr= (pressure(i,j,k))*0.0075006
             !second absorption band     
             do l = 1,39 
                qy   = (phi0(l) * rd(l))/(rd(l) + ptorr*phi0(l)) 
                phot_dat(region)%qy_ch3cocho(i,j,k,l) = min(1.0,qy)
             end do
             ! switch to the formulation by Chen et al, JPC, 104, 11126-11131, 2000 for wav > 420nm
             do l = 40,52
                qy   = phi0(l)/(1.+(phi0(l)*rd(l)*ptorr))
                phot_dat(region)%qy_ch3cocho(i,j,k,l) = min(1.0,qy)
             enddo
             !
             ! calculate the second photolysis channel for acetone
             !
             do l=26,35
                table_pos=int((temperature(i,j,k)-(temp_ind(1)-0.0*5))/5)+1
                table_pos=min(max(1,table_pos),34)  ! bug fix for temperatures below 175.5
                qy=(1-qy_co_look(l,table_pos))/(1.0+a1_acet(l,table_pos)*density(i,j,k))
                qyacet(i,j,k,l)=max(0.0,qy)   
             enddo
             do l=36,89
                ! determine constants for calculating qy
                table_pos=int((temperature(i,j,k)-(temp_ind(1)-0.0*5))/5)+1
                table_pos=min(max(1,table_pos),34)  ! bug fix for temperatures below 175.5
                a1=1.0+(density(i,j,k)*a4_acet(l,table_pos))+a3_acet(l,table_pos)
                a2=1.0+(density(i,j,k)*a2_acet(l,table_pos))+a3_acet(l,table_pos)
                a4=1.0+(density(i,j,k)*a4_acet(l,table_pos))
                qy=(1-qy_co_look(l,table_pos))*(a1/(a2*a4))
                qyacet(i,j,k,l)=max(0.0,qy)     
             enddo
          end do  ! longitude
       end do ! latitude
    end do ! level

    deallocate(density)
    deallocate(pressure)
    nullify(temperature)
    nullify(mass)    
    nullify(phlb)
    nullify(qyacet)
    
  END SUBROUTINE UPDATE_CSQY
  !EOC


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    PIFM
  !
  ! !DESCRIPTION: 
  !                                                                      *
  !            PRACTICAL IMPROVED FLUX METHOD (PIFM)                     *
  !                 to calculate actinic fluxes                          *
  !     Zdunkowski,Welch,Korb: Beitr. Phys. Atmosph. vol. 53, p. 147 ff  *
  !                                                                      *
  !     This version is not suitable for calculation for conserving      *
  !     scattering (W0=1). W0 is limited to W0 <= 1. - 1.E-15.           *
  !     For W0 = 1, AL(4) and AL(5) has to be calculated differently.    *
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE PIFM(region, zangle, alb, cst_o3_col, dv2, dv3, taua_aer_col, taus_aer_col, paer_col, fact)
    !
    ! !USES:
    !
    use dims,  only : lm
    use binas, only : pi
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in) :: region
    real,    intent(in) :: zangle    ! zenith angle
    real,    intent(in) :: alb       ! albedo 

    real, dimension(lm(region),maxwav)  :: cst_o3_col         ! temp dep o3 cross sections
    real, dimension(lm(region))         :: dv2, dv3             ! differential column info
    real, dimension(lm(region),nbands_trop,ngrid) :: taua_aer_col, taus_aer_col   ! optical depth aerosols     
    real, dimension(lm(region),nbands_trop,ngrid) :: paer_col
    real, dimension(lm(region),nbands_trop) :: fact             !actinic flux 
    !
    ! !REVISION HISTORY: 
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    real, dimension(3*(lm(region)+1))   :: rw,  & !flux array for cloudy sky
                                           rf     !flux array for clear sky 
    real, dimension(lm(region))     :: tu1,  &   !parallel solar flux   
                                       tu2       !matrix coefficient 

    real, dimension(lm(region),5)     :: al
    real, dimension(lm(region)+1,nbands_trop):: sd, fd, fu
    real, dimension(0:nmom)         :: pray  ! phase functions
    real                            :: taus, taua, tautot, g
    real, dimension(lm(region),nbands_trop)  :: taua_clr, taus_clr  ! optical depth clear sky    

    !     diffusivity factor                                                
    real, parameter ::  u  = 2., delu0 = 1.e-3, resonc = 1.e-6
    integer             :: grid, l, k, n, nl, k3,  maxlev,  maxlay2, maxlev3       
    real                        :: p1, w0, f, b0, bu0                              
    real                        :: eps, factor, ueps2, smooth1, smooth2    
    real                        :: alph1, alph2, alph3, alph4                              
    real                        :: gam1, gam2, e, rm, tds1, td1, td2               
    real                        :: fact_bot, fact_top, arg                                         
    !downward diffuse flux            
    !upward diffuse flux 
    real                        :: a, b, c, gamma, u0
    real                        :: cs_o2

    !     internal integer variable                                         
    maxlev  = lm(region) + 1     
    maxlay2 = 2 * lm(region)     
    maxlev3 = 3 * maxlev

    ! intitialise arrays
    al = 0. ; fd = 0. ; sd = 0. ; fu = 0. ; tu1 = 0. ; tu2 = 0. ; rw = 0. ; rf = 0.

    !initialise output (actinic flux)
    fact = 0.

    !-----correction of the air mass factor---------------------------------  

    !   F. Kasten and T. Young, Revised optical air mass tabels and
    !   approximation formula (1989) Aplied Optics Vol 28, No. 22 P. 4735   
    !   and J. Lenoble, atmospheric radiative transfer (1993), P. 236

    a = 0.50572 ; b = 6.07995 ; c = 1.63640

    ! define air mass factor in mu
    gamma = 90. - zangle
    u0 = sin(gamma*pi/180.) + a *(gamma+b)**(-c)
    u0 = min(1.,u0)

    !===================================================================    
    !     calculation of the matrix coefficients a1,...,a5                  
    !=================================================================== 

    ! determine phase functions
    pray(0) = 1. ; pray(1) = 0. ; pray(2)=0.1 ; pray(3:nmom) = 0. 

    ! reverse order for input parameters

    dv2(1:lm(region)) = dv2(lm(region):1:-1)
    dv3(1:lm(region)) = dv3(lm(region):1:-1)

    ! do not calculate for band #1  
    do l  = 1,nbands_trop                                                           
       if (zangle <= 71.) then
          nl = lmid(l)
          grid = 1
          cs_o2 = 7.322e-24
       elseif (zangle > 71. .and. zangle<=sza_limit) then
          nl = lmid_gridA(l)
          grid = 2
          cs_o2 = 6.608e-24
       endif

       ! reverse order for input parameters
       cst_o3_col(1:lm(region),nl) = cst_o3_col(lm(region):1:-1,nl)
       taus_aer_col(1:lm(region),l,grid) = taus_aer_col(lm(region):1:-1,l,grid)
       taua_aer_col(1:lm(region),l,grid) = taua_aer_col(lm(region):1:-1,l,grid)
       paer_col(1:lm(region),l,grid) = paer_col(lm(region):1:-1,l,grid)

       do k  = 1,lm(region)

          ! Calculate optical depth (clear-sky)
          if(nl<18) then 
             taua_clr(k,l) = cs_o2*dv2(k) + cst_o3_col(k,nl)*dv3(k)
          else
             taua_clr(k,l) = cst_o3_col(k,nl)*dv3(k)
          endif
          taus_clr(k,l) = cs_ray(nl)*1./0.2095*dv2(k)

          taus = taus_clr(k,l) + taus_aer_col(k,l,grid)
          taua = taua_clr(k,l) + taua_aer_col(k,l,grid) 

          tautot = taus + taua

          g = 0.
          if (taus > 0.) &
               g = (paer_col(k,l,grid)* taus_aer_col(k,l,grid) + pray(1)* taus_clr(k,l) )/taus

          if (tautot .ne. 0.) then 
             w0=taus / tautot 
          else 
             w0=0. 
          end if

          w0 = min(w0,1.-1.e-15)

          !      P1: first expansion coefficient of the phase function    

          p1=3.*g 

          !      F: fraction of radiation contained in diffraction peak    
          f=g**2 

          !      B0:  fractional mean backward scattering coefficient    
          !           of diffuse light                             
          !      BU0: backward scattering coefficient of primary scattered   
          !           parallel solar light                         
          !      for small P1 SMOOTH1,SMOOTH2 manage the smooth change of B0 
          !      BU0 to 0                                          
          if (p1 <= 0.1) then 
             smooth1=1.33333333-p1*3.3333333 
             smooth2=10.*p1 
          else 
             smooth1=1. 
             smooth2=1. 
          end if

          b0=(3.-p1)/8.  *smooth1 
          bu0=0.5-u0/4.*(p1-3.*f)/(1.-f)  *smooth2 

          !      alpha coefficient                                 
          alph1=u*(1.-(1.-b0)*w0)
          alph2=u*b0*w0
          alph3=w0*bu0*(1-f)
          alph4=w0*(1.-bu0)*(1-f)

          !      epsilon and gamma coefficient                     
          eps=sqrt(alph1**2-alph2**2)

          factor=1.-w0*f 

          !      check for resonance condition in GAM1 and GAM2, if fulfil th
          !      chance U0(J) and calculate UEPS2, BU0, ALPH3, ALPH4 again.  

          ueps2=(u0 *eps)**2
          arg = ueps2-factor**2

          if (arg < 0.) arg = -1. * arg

          if (arg < resonc) then
             if (ueps2.lt.factor**2) then
                u0 = u0 - delu0
             else
                u0 = u0 + delu0
             end if
             ueps2=(u0 * eps)**2
             bu0=0.5-u0/4.*(p1-3.*f)/(1.-f)  *smooth2
             alph3=w0*bu0*(1-f)
             alph4=w0*(1.-bu0)*(1-f)
          end if

          gam1=( factor*alph3-u0 * (alph1*alph3+alph2*alph4) ) * &
               1/(factor**2-ueps2)
          gam2=(-factor*alph4-u0 * (alph1*alph4+alph2*alph3) ) * &
               1/(factor**2-ueps2)

          e=exp(-eps*tautot)
          rm=alph2/(alph1+eps)

          al(k,4)=e*(1.-rm**2.)/(1.-e**2. * rm**2.)
          al(k,5)=rm*(1.-e**2.)/(1.-e**2. * rm**2.)
          al(k,1)=exp(-factor*tautot/u0)
          al(k,2)=-al(k,4)*gam2-al(k,5)*gam1*al(k,1) + gam2*al(k,1)
          al(k,3)=-al(k,5)*gam2-al(k,4)*gam1*al(k,1) + gam1
       enddo

       !====================================================================   
       ! matrix inversion                                                      
       !====================================================================   

       !        direct solution of the first four equations                    
       rf(1) = u0 * flux(nl) 
       rf(2) = 0.

       !       5th to 10th equation: bring matrix elements on the left of the m
       !       diagonal to the rhs:  save elements on the right of the main    
       !                             diagonal in array -tu(l,1)                
       rf(3) = al(1,3) * rf(1) 
       rf(4) = al(1,1) * rf(1) 
       rf(5) = al(1,2) * rf(1) 

       tu1(1) = al(1,4) 
       tu2(1) = al(1,5)

       ! blocks of 6 equations: eliminate left matrix elements, save righ
       !                              matrix elements in array -tu(l,i),       
       !                              calculate rhs.                           
       do k=2,lm(region)
          k3=3*k
          td1 = 1./(1.-al(k,5)*tu2(k-1))
          td2 = al(k,4)*tu2(k-1)
          tu1(k) = td1*al(k,4)
          tu2(k) = al(k,5) + td2*tu1(k)
          rf(k3) = td1 * (al(k,3)*rf(k3-2) + al(k,5)*rf(k3-1))
          rf(k3+1)= al(k,1)*rf(k3-2)
          rf(k3+2)= al(k,2)*rf(k3-2) + al(k,4)*rf(k3-1)+td2*rf(k3)
       end do

       !    last two equations: the same as before                         
       tds1       = 1. / (1.-alb*tu2(lm(region)))
       rf(maxlev3) = tds1 * (alb * rf(maxlev3-2)+   &
            alb * rf(maxlev3-1))

       !        now we have created an upper triangular matrix the elements of 
       !        are -tu(l,i), 0, or 1 (in the main diagonal). the 0 and 1 eleme
       !        are not stored in an array. let us solve the system now and sto
       !        results in the arrays rf (fluxes clear sky) and rw (fluxes clou

       do k=lm(region),1,-1
          k3=3*k

          rf(k3+2) = rf(k3+2) + tu2(k)*rf(k3+3)  
          rf(k3) = rf(k3) + tu1(k)*rf(k3+3)    

          sd(k+1,l)  = rf(k3+1)                                    
          fd(k+1,l)  = rf(k3+2)                                    
          fu(k+1,l)  = rf(k3+3)                                    
       end do

       sd(1,l)  = rf(1)    
       fd(1,l)  = rf(2)    
       fu(1,l)  = rf(3)     

       !        calculate the actinic flux                                     

       fact_top = sd(1,l)/u0 + u * fd(1,l)   + u * fu(1,l)

       do k = 1,lm(region)

          fact_bot = sd(k+1,l)/u0 +     &                
               u * fd(k+1,l) + u * fu(k+1,l)                                             
          fact(k,l) = max(0. ,(fact_top + fact_bot)/2.)
          fact_top = fact_bot
       end do

       ! swap vertical levels to the default order of the model
       fact(1:lm(region),l) = fact(lm(region):1:-1,l)
       cst_o3_col(1:lm(region),nl) = cst_o3_col(lm(region):1:-1,nl)

    enddo    ! spectral bands (nbands)

  END SUBROUTINE PIFM
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:   PIFM_RAN
  !
  ! !DESCRIPTION: Pifm method including random overlap method for clouds
  ! 
  !************************************************************************
  !*            PRACTICAL IMPROVED FLUX METHOD (PIFM)                     * 
  !*                 to calculate actinic fluxes                          *
  !*     Zdunkowski,Welch,Korb: Beitr. Phys. Atmosph. vol. 53, p. 147 ff  *
  !*                                                                      *
  !*           Cloud effects added using the method of :                  *
  !*     Geleyn and Hollingworth, Contribs. Atms.Phys,52(1),1979          *
  !*                                                                      *
  !************************************************************************      
  !*     This version is not suitable for calculation for conserving      *
  !*     scattering (W0=1). W0 is limited to W0 .le. 1. - 1.E-15.         *
  !*     For W0 = 1, AL(4) and AL(5) has to be calculated differently.    *
  !************************************************************************
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE PIFM_RAN(region, zangle, alb, cst_o3_col, dv2, dv3,       &
                      taua_cld_col, taus_cld_col, pcld_col,            &
                      taua_aer_col, taus_aer_col, paer_col, fact, frac )
    !
    ! !USES:
    !
    use dims,  only : lm
    use binas, only : pi
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)               :: region
    real,    intent(in)               :: zangle     ! zenith angle
    real,    intent(in)               :: alb        ! albedo
    
    real,dimension(lm(region))        :: frac       ! cloud fraction per layer (0->1)
    real,dimension(lm(region),maxwav) :: cst_o3_col ! o3 cross sections
    real,dimension(lm(region))        :: dv2, dv3   ! differential column info

    real,dimension(lm(region))                   :: taua_cld_col, taus_cld_col ! optical depth clouds              
    real,dimension(lm(region),nbands_trop)       :: taua_clr_col, taus_clr_col ! optical depth clear sky
    real,dimension(lm(region),nbands_trop,ngrid) :: taua_aer_col, taus_aer_col ! optical depth aerosols   
    real,dimension(lm(region),nbands_trop)       :: fact     ! actinic flux 
    real,dimension(lm(region),nbands_trop,ngrid) :: paer_col ! aerosol phase functions
    real,dimension(lm(region))                   :: pcld_col ! cloud phase functions
    real,dimension(0:nmom)                       :: pray     ! rayleigh phase function
    !
    ! !REVISION HISTORY: 
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    real,dimension(3*(lm(region)+1)) :: rw,  & ! flux array for cloudy sky
                                        rf     ! flux array for clear sky 
    real,dimension(lm(region)) :: tu1, & ! parallel solar flux   
                                  tu2    ! matrix coefficient 

    real,dimension(lm(region),5) :: al
    real,dimension(lm(region)+1,nbands_trop):: sd, fd, fu
    real,dimension(lm(region),nbands_trop)  :: TS_PI_CLR,TA_PI_CLR, G_PI_CLR, TS_PI_CLD
    real,dimension(lm(region))              :: TA_PI_CLD, G_PI_CLD

    real      ::  U,DELU0,RESONC,a,b,c,gamma,u0,cs_o2 

    integer   :: grid,j,k,l,maxlev,maxlay2,maxlev3,nl,k3 

    real :: w0,p1,F,smooth1,smooth2,b0,bu0,alph1,alph2,alph3,alph4,tautot,eps,factor
    real :: rm,gam1,gam2,e,tauscat,td1,td2,tds1,ueps2,g,arg,fact_bot, fact_top

    !     diffusivity factor
    DATA   U / 2./

    DATA    DELU0 /1.E-3/ 

    DATA    RESONC /1.E-6/

    real    :: AL_CLR_1,AL_CLR_2,AL_CLR_3,AL_CLR_4,AL_CLR_5   !matrix coefficient for clear sky
    real    :: AL_CLD_1,AL_CLD_2,AL_CLD_3,AL_CLD_4,AL_CLD_5   !matrix coefficient for cloudy sky


    !-----correction of the air mass factor---------------------------------
    !   F. Kasten and T. Young, Revised optical air mass tables and               
    !   approximation formula (1989) Aplied Optics Vol 28, No. 22 P. 4735   
    !   and J. Lenoble, atmospheric radiative transfer (1993), P. 236     

    a = 0.50572 ; b = 6.07995 ; c = 1.63640

    ! define air mass factor in mu
    gamma = 90. - zangle                                                                         
    u0 = sin(gamma*pi/180.) + a *(gamma+b)**(-c)         
    u0 = min(1.,u0)    

    !---------------------------------------------------------------------     
    !     internal integer variable 

    MAXLEV  = lm(region) + 1
    MAXLAY2 = 2 * lm(region)
    MAXLEV3 = 3 * MAXLEV

    ! initialize

    fact = 0. ; sd = 0 ; fd = 0. ;  fu = 0. ; td1 = 0.

    rf = 0. ; rw = 0. ; al = 0. ;   tu1 = 0. ; tu2 = 0.

    ! determine phase functions
    pray(0) = 1. ; pray(1) = 0. ; pray(2)=0.1 ; pray(3:nmom) = 0. 

    dv2(1:lm(region)) = dv2(lm(region):1:-1)
    dv3(1:lm(region)) = dv3(lm(region):1:-1)
    frac(1:lm(region)) = frac(lm(region):1:-1)

    taua_cld_col(1:lm(region)) = taua_cld_col(lm(region):1:-1)

    taus_cld_col(1:lm(region)) = taus_cld_col(lm(region):1:-1)

    pcld_col(1:lm(region)) = pcld_col(lm(region):1:-1)

    do l  = 1,nbands_trop 

       if (zangle <= 71.) then
          nl = lmid(l)
          grid = 1
          cs_o2 = 7.322e-24
       elseif (zangle > 71. .and. zangle<=sza_limit) then
          nl = lmid_gridA(l)
          grid = 2
          cs_o2 = 6.608e-24
       endif

       ! reverse order for input parameters
       cst_o3_col(1:lm(region),nl) = cst_o3_col(lm(region):1:-1,nl)
       taus_aer_col(1:lm(region),l,grid) = taus_aer_col(lm(region):1:-1,l,grid)
       taua_aer_col(1:lm(region),l,grid) = taua_aer_col(lm(region):1:-1,l,grid)
       paer_col(1:lm(region),l,grid) = paer_col(lm(region):1:-1,l,grid)

       ! fill array for absorption and scattering components before performing
       ! calculations.   
       do K = 1,lm(region)

          taus_clr_col(k,l) = cs_ray(nl)*1./0.2095*dv2(k)

          if(nl<18) then
             taua_clr_col(k,l) = cs_o2*dv2(k) + cst_o3_col(k,nl)*dv3(k)
          else
             taua_clr_col(k,l) = cst_o3_col(k,nl)*dv3(k)
          endif

          TS_PI_CLR(K,L) = TAUS_CLR_col(K,L)+TAUS_AER_col(K,L,grid)
          TA_PI_CLR(K,L) = TAUA_CLR_col(K,L)+TAUA_AER_col(K,L,grid)
          TS_PI_CLD(K,L) = TAUS_CLD_col(K)
          TA_PI_CLD(K) = TAUA_CLD_col(K)

          IF (TAUS_CLR_col(K,L)+TAUS_AER_col(K,L,grid)>0.) THEN
             G_PI_CLR(K,L)  = (PRAY(1)*TAUS_CLR_col(K,L) +   &
                  PAER_col(k,l,grid)*TAUS_AER_col(K,L,grid))/  &
                  (TAUS_CLR_col(K,L)+TAUS_AER_col(K,L,grid))                
          ELSE
             G_PI_CLR(K,L)  = 0.                 

          ENDIF
          G_PI_CLD(K)  = PCLD_col(k) 

       enddo

       do K = 1,lm(region)       

          !*****  first: clear sky *******************************************

          TAUTOT = TS_PI_CLR(K,L)+TA_PI_CLR(K,L)


          if (tautot /= 0.) then  
             W0=TS_PI_CLR(K,L) / TAUTOT
          ELSE
             W0=0.
          ENDIF

          W0 = MIN(W0,1.-1e-15)

          !           P1: first expansion coefficient of the phase function

          P1=3.*G_PI_CLR(K,L)                            

          !          F: fraction of radiation contained in diffraction peak

          F=G_PI_CLR(K,L)**2                        

          !           B0:  fractional mean backward scattering coefficient of diffuse light
          !           BU0: backward scattering coefficient of primary scattered parallel solar light

          !           for small P1 SMOOTH1,SMOOTH2 manage the smooth change of B0 and
          !           BU0 to 0

          IF (P1<=0.1) THEN
             SMOOTH1=1.33333333-P1*3.3333333
             SMOOTH2=10.*P1
          ELSE
             SMOOTH1=1
             SMOOTH2=1
          ENDIF

          B0=(3.-P1)/8.  *SMOOTH1                                 
          BU0=0.5-U0/4.*(P1-3.*F)/(1.-F)  *SMOOTH2    

          !           alpha coefficient

          ALPH1=U*(1.-(1.-B0)*W0)     
          ALPH2=U*B0*W0               
          ALPH3=W0*BU0*(1-F)
          ALPH4=W0*(1.-BU0)*(1-F)

          !           epsilon and gamma coefficient

          EPS=SQRT(ALPH1**2-ALPH2**2)          

          FACTOR=1.-W0*F

          !           check for resonance condition in GAM1 and GAM2, if fulfil then
          !           chance U0 and calculate UEPS2, BU0, ALPH3, ALPH4 again.

          UEPS2=(U0*EPS)**2
          arg = ueps2-factor**2
          if (arg < 0.) then
             arg = -1. * arg
          endif
          if (arg < resonc) then        
             IF(UEPS2.LT.FACTOR**2) THEN
                U0=U0-DELU0
             ELSE
                U0=U0+DELU0
             ENDIF
             UEPS2=(U0*EPS)**2
             BU0=0.5-U0/4.*(P1-3.*F)/(1.-F)  *SMOOTH2    
             ALPH3=W0*BU0*(1-F)
             ALPH4=W0*(1.-BU0)*(1-F)
          ENDIF

          GAM1=( FACTOR*ALPH3-U0*(ALPH1*ALPH3+ALPH2*ALPH4) ) * &
               &           1/(FACTOR**2-UEPS2)
          GAM2=(-FACTOR*ALPH4-U0*(ALPH1*ALPH4+ALPH2*ALPH3) ) * &
               &           1/(FACTOR**2-UEPS2)

          E=EXP(-EPS*TAUTOT)
          RM=ALPH2/(ALPH1+EPS)

          AL_CLR_4= E*(1-RM**2)/(1-E**2 * RM**2)
          AL_CLR_5= RM*(1-E**2)/(1-E**2 * RM**2)
          AL_CLR_1= EXP(-FACTOR*TAUTOT/U0)

          AL_CLR_2=-AL_CLR_4*GAM2-AL_CLR_5*GAM1*     &
               &                AL_CLR_1+GAM2*AL_CLR_1  
          AL_CLR_3=-AL_CLR_5*GAM2-AL_CLR_4*GAM1*     &
               &                AL_CLR_1+GAM1    

          !******************* second: cloudy sky *****************************
          ! For ECMWF input there is the possibility that cloud fraction occurs without a 
          ! corresponding value for lwc

          IF(  FRAC(K) > 0.02 .and. TS_PI_CLD(K,L) > 1.e-5 ) THEN


             TAUSCAT = TS_PI_CLR(K,L) + TAUS_CLD_col(K)
             TAUTOT  = TA_PI_CLR(K,L) + TAUA_CLD_col(K)+TAUSCAT
             G       = G_PI_CLD(K)*TAUS_CLD_col(K)/TAUSCAT 

             IF (TAUTOT > 0.) THEN
                W0=TAUSCAT/TAUTOT    
             ELSE
                W0=0.
             ENDIF

             W0 = MIN(W0,1.-1e-15)

             P1=3.*G
             F=G**2                        

             IF (P1<0.1) THEN
                SMOOTH1=1.33333333-P1*3.3333333
                SMOOTH2=10.*P1
             ELSE
                SMOOTH1=1
                SMOOTH2=1
             ENDIF

             B0=(3.-P1)/8.  *SMOOTH1                                 
             BU0=0.5-U0/4.*(P1-3.*F)/(1.-F)  *SMOOTH2    

             ALPH1=U*(1.-(1.-B0)*W0)     
             ALPH2=U*B0*W0               
             ALPH3=W0*BU0*(1-F)
             ALPH4=W0*(1.-BU0)*(1-F)

             EPS=SQRT(ALPH1**2-ALPH2**2)  



             FACTOR=1.-W0*F

             UEPS2=(U0*EPS)**2
             arg = ueps2-factor**2
             if (arg < 0.) then
                arg = -1. * arg
             endif

             if (arg < resonc) then 
                IF(UEPS2.LT.FACTOR**2) THEN
                   U0=U0-DELU0
                ELSE
                   U0=U0+DELU0
                ENDIF
                UEPS2=(U0*EPS)**2
                BU0=0.5-U0/4.*(P1-3.*F)/(1.-F)  *SMOOTH2    
                ALPH3=W0*BU0*(1-F)
                ALPH4=W0*(1.-BU0)*(1-F)
             ENDIF

             GAM1=( FACTOR*ALPH3-U0*(ALPH1*ALPH3+ALPH2*ALPH4) ) * 1/(FACTOR**2-UEPS2)
             GAM2=(-FACTOR*ALPH4-U0*(ALPH1*ALPH4+ALPH2*ALPH3) ) * 1/(FACTOR**2-UEPS2)

             E=EXP(-EPS*TAUTOT)
             RM=ALPH2/(ALPH1+EPS)

             AL_CLD_4=E*(1-RM**2)/(1-E**2 * RM**2)
             AL_CLD_5=RM*(1-E**2)/(1-E**2 * RM**2)
             AL_CLD_1=EXP(-FACTOR*TAUTOT/U0) 

             AL_CLD_2=-AL_CLD_4*GAM2-AL_CLD_5*GAM1*AL_CLD_1+GAM2*AL_CLD_1  
             AL_CLD_3=-AL_CLD_5*GAM2-AL_CLD_4*GAM1*AL_CLD_1+GAM1   

             AL(K,1) =(1-FRAC(K))*AL_CLR_1 + FRAC(K)*AL_CLD_1
             AL(K,2) =(1-FRAC(K))*AL_CLR_2 + FRAC(K)*AL_CLD_2
             AL(K,3) =(1-FRAC(K))*AL_CLR_3 + FRAC(K)*AL_CLD_3
             AL(K,4) =(1-FRAC(K))*AL_CLR_4 + FRAC(K)*AL_CLD_4
             AL(K,5) =(1-FRAC(K))*AL_CLR_5 + FRAC(K)*AL_CLD_5
          ELSE

             AL(K,1) =  AL_CLR_1
             AL(K,2) =  AL_CLR_2
             AL(K,3) =  AL_CLR_3
             AL(K,4) =  AL_CLR_4
             AL(K,5) =  AL_CLR_5

          ENDIF
       enddo !k

       !====================================================================
       ! matrix inversion
       !====================================================================

       !        direct solution of the first four equations       

       RF(1) = U0*FLUX(NL)
       RF(2) = 0.                              

       !       5th to 10th equation: bring matrix elements on the left of the main
       !       diagonal to the rhs:  save elements on the right of the main 
       !                             diagonal in array -tu(l,1)    

       RF(3) = AL(1,3) * RF(1)         
       RF(4) = AL(1,1) * RF(1)         
       RF(5) = AL(1,2) * RF(1)         

       TU1(1) = AL(1,4)        
       TU2(1) = AL(1,5)        

       !       blocks of 6 equations: eliminate left matrix elements, save right 
       !                              matrix elements in array -tu(l,i), 
       !                              calculate rhs.

       DO  K=2,lm(region)      
          K3=3*K              
          TD1       = 1./(1.-AL(K,5)*TU2(K-1))                  
          TD2       = AL(K,4)*TU2(K-1)        
          TU1(K)  = TD1*AL(K,4)
          TU2(K)  = AL(K,5) + TD2*TU1(K) 
          RF(K3)  = TD1 * (AL(K,3)*RF(K3-2) + AL(K,5)*RF(K3-1))     
          RF(K3+1)= AL(K,1)*RF(K3-2)      
          RF(K3+2)= AL(K,2)*RF(K3-2) + AL(K,4)*RF(K3-1)+TD2*RF(K3)
       enddo

       !        last two equations: the same as before 

       TDS1          = 1. / (1.-ALB*TU2(lm(region)))  
       rf(maxlev3) = tds1 * (alb * rf(maxlev3-2)+   &
            alb * rf(maxlev3-1))

       !        now we have created an upper triangular matrix the elements of which
       !        are -tu(l,i), 0, or 1 (in the main diagonal). the 0 and 1 elements 
       !        are not stored in an array. let us solve the system now and store the
       !        results in the arrays rf (fluxes clear sky) and rw (fluxes cloudy sky)

       DO  K=lm(region),1,-1         
          K3=3*K                  
          RF(K3+2) = RF(K3+2) + TU2(K)*RF(K3+3) 
          RF(K3)   = RF(K3) + TU1(K)*RF(K3+3)   
          SD(K+1,L)  = RF(K3+1)
          FD(K+1,L)  = RF(K3+2) 
          FU(K+1,L)  = RF(K3+3)
       enddo ! K  

       SD(1,L)  = RF(1) 
       FD(1,L)  = RF(2) 
       FU(1,L)  = RF(3)

       !        calculate the actinic flux                                     

       fact_top = sd(1,l)/u0 + u * fd(1,l)   + u * fu(1,l)

       do k = 1,lm(region)
          fact_bot = sd(k+1,l)/u0 + u * fd(k+1,l) + u * fu(k+1,l)
          fact(k,l) = max(0. ,(fact_top + fact_bot)/2.)
          fact_top = fact_bot   
       end do ! K

       ! swap vertical levels to the default order of the model
       fact(1:lm(region),l) = fact(lm(region):1:-1,l)
       cst_o3_col(1:lm(region),nl) = cst_o3_col(lm(region):1:-1,nl)

    enddo ! wavelength


    return                                                  

  END SUBROUTINE PIFM_RAN
  !EOC

  !------------------------------------------------------------------------------
  !                    TM5                                                      !
  !------------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:   SUNDIS
  !
  ! !DESCRIPTION: 
  !-----------------------------------------------------------------------------*
  !=  purpose:                                                                 =*
  !=  calculate earth-sun distance variation for a given date.  based on       =*
  !=  fourier coefficients originally from:  spencer, j.w., 1971, fourier      =*
  !=  series representation of the position of the sun, search, 2:172          =*
  !-----------------------------------------------------------------------------*
  !=  parameters:                                                              =*
  !=  idate  - integer, specification of the date, from yymmdd              (i)=*
  !=  esrm2  - real, variation of the earth-sun distance                    (o)=*
  !=           esrm2 = (average e/s dist)^2 / (e/s dist on day idate)^2        =*
  !-----------------------------------------------------------------------------*
  !=  edit history:                                                            =*
  !=  01/95  changed computation of trig function values                       =*
  !-----------------------------------------------------------------------------*
  != this program is free software;  you can redistribute it and/or modify     =*
  != it under the terms of the gnu general public license as published by the  =*
  != free software foundation;  either version 2 of the license, or (at your   =*
  != option) any later version.                                                =*
  != the tuv package is distributed in the hope that it will be useful, but    =*
  != without any warranty;  without even the implied warranty of merchantibi-  =*
  != lity or fitness for a particular purpose.  see the gnu general public     =*
  != license for more details.                                                 =*
  != to obtain a copy of the gnu general public license, write to:             =*
  != free software foundation, inc., 675 mass ave, cambridge, ma 02139, usa.   =*
  !-----------------------------------------------------------------------------*
  != to contact the authors, please mail to:                                   =*
  != sasha madronich, ncar/acd, p.o.box 3000, boulder, co, 80307-3000, usa  or =*
  != send email to:  sasha@ucar.edu                                            =*
  !-----------------------------------------------------------------------------*
  != copyright (c) 1994,95,96  university corporation for atmospheric research =*
  !-----------------------------------------------------------------------------*
  !\\
  !\\
  ! !INTERFACE:
  !
  REAL FUNCTION SUNDIS( imonth, iday)
    !
    ! !USES:
    !
    use binas, only : pi
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in) :: iday, imonth
    !
    ! !REVISION HISTORY: 
    !
    !EOP
    !------------------------------------------------------------------------------
    !BOC

    ! internal:
    integer mday, month, jday
    real dayn, thet0
    real sinth, costh, sin2th, cos2th

    integer,dimension(12) ::imn=(/ 31,28,31,30,31,30,31,31,30,31,30,31 /) 

    ! start

    ! parse date to find day number (julian day)

    mday = 0
    do  month = 1, imonth-1
       mday = mday + imn(month)               
    end do
    jday = mday + iday
    dayn = float(jday - 1) + 0.5

    ! define angular day number and compute esrm2:

    thet0 = 2.*pi*dayn/365.

    ! calculate sin(2*thet0), cos(2*thet0) from
    ! addition theoremes for trig functions for better
    ! performance;  the computation of sin2th, cos2th
    ! is about 5-6 times faster than the evaluation
    ! of the intrinsic functions sin and cos
    !
    sinth = sin(thet0)
    costh = cos(thet0)
    sin2th = 2.*sinth*costh
    cos2th = costh*costh - sinth*sinth
    sundis  = 1.000110 +   0.034221*costh  +  0.001280*sinth + &
         0.000719*cos2th +  0.000077*sin2th
    !

  END FUNCTION SUNDIS
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:   OZONE_INFO_ONLINE
  !
  ! !DESCRIPTION: calculate, from an ozone field, the overhead ozone at all
  !               grid points. 
  !               Optical depth clear sky conditions
  !               upper values are given by the ozone climatology above the
  !               highest model mid-level
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE OZONE_INFO_ONLINE( region )
    !
    ! !USES:
    !
    use binas,           only : Avog, xm_air, grav
    use dims
!    use photolysis_data, only : phot_dat
    use global_data,     only : region_dat, mass_dat
    use chem_param,      only : xmo3, io3
    use meteodata,       only : phlb_dat
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)  :: region
    !
    ! !REVISION HISTORY: 
    !      Jan 2003 - MK  - adapted for new TM5 memory structure.
    !      Jul 2008 - JEW - new ozone info unit coupled for online calculations
    !   26 Apr 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    real, parameter  :: conv = 0.1*Avog/xmo3   !  from kg/m2 --> #/cm2
    real, parameter  :: sp = Avog*1.e-4*0.2095/(xm_air*grav)

    integer :: i, j, l, lmr, i1, i2, j1, j2
    real,dimension(:,:,:),allocatable :: o3_overhead_online, o2_overhead   ! #/cm2
    real,dimension(:),    pointer     :: ozone_klimtop !in #cm-2
    real,dimension(:,:),  pointer     :: v3_surf       ! #/cm2
    real,dimension(:,:,:),pointer     :: v2, v3        ! #/cm2 
    real,dimension(:,:,:,:),pointer   :: o3            ! #/cm2 
    real,dimension(:,:,:),pointer     :: phlb
    real,dimension(:),    pointer     :: area

    
    call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
    lmr = lm(region)
    
    ! allocate on all prcessors so result maybe scattered   
    allocate(o3_overhead_online(i1:i2, j1:j2, lm(region))) ; o3_overhead_online = 0.0
    allocate(o2_overhead       (i1:i2, j1:j2, lm(region))) ; o2_overhead = 0.0

    v2 => phot_dat(region)%v2
    v3 => phot_dat(region)%v3
    v3_surf => phot_dat(region)%v3_surf

    area => region_dat(region)%dxyp
    phlb =>   phlb_dat(region)%data
    o3   =>   mass_dat(region)%rm    ! (i1:i2, j1:j2, 1:lmr, io3)

    ! CALCULATE COLUMNS
    ! -----------------

    ! (1) top
    ! --------
    ! o3du is not used anymore, except in cases where stratosphere is not resolved (EC-EARTH)

    if (with_o3du) then  ! use fortuin-kelder
       ozone_klimtop => phot_dat(region)%o3klim_top

       do j=j1,j2
       do i=i1,i2
             o3_overhead_online(i,j,lm(region)) = ozone_klimtop(j)  
       end do
       end do
       
    else ! from TM4-routine
       
       do j=j1,j2
       do i=i1,i2
             o3_overhead_online(i,j,lm(region)) =  0.5*conv*o3(i,j,lm(region),io3)/area(j)
       end do
       end do
       
    endif
    
    ! (2) rest
    ! --------
    do l = lm(region)-1,1,-1
       do j = j1,j2
       do i = i1,i2
             o3_overhead_online(i,j,l) = o3_overhead_online(i,j,l+1) + &
                  0.5*conv*(o3(i,j,l,io3)+o3(i,j,l+1,io3))/area(j)         
       enddo
       enddo
    enddo

    ! Define surface ozone column field for diagnostic purposes
    do j=j1,j2
    do i=i1,i2
          v3_surf(i,j) = o3_overhead_online(i,j,1) + 0.5*conv*o3(i,j,1,io3)/area(j)
    end do
    end do

    ! Boundaries
    v3(:,:,1) = o3_overhead_online(:,:,1)  

    ! now transform o3 column from layers to levels
    do l = 2,lm(region)
       v3(:,:,l) = ( o3_overhead_online(:,:,l) + o3_overhead_online(:,:,l-1) ) * 0.5 
    enddo

    ! TOA
    do j=j1,j2
    do i=i1,i2
          v3(i,j,lm(region)+1) = 0.5*v3(i,j,lm(region))
    enddo
    enddo

    ! determine oxygen columns (can directly be calculated on levels)
    do l = 1,lm(region)
    do j = j1,j2
    do i = i1,i2
             v2(i,j,l) = phlb(i,j,l)*sp
    enddo
    enddo
    enddo

    ! top boundary for v2 
    v2(:,:,lm(region)+1) = 0.5*v2(:,:,lm(region))

    nullify(o3)
    nullify(ozone_klimtop)
    nullify(area)
    nullify(phlb)

    nullify(v2)
    nullify(v3)
    nullify(v3_surf)

    deallocate(o3_overhead_online) 
    deallocate(o2_overhead) 

  END SUBROUTINE OZONE_INFO_ONLINE
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    PHOTORATES_TROPO
  !
  ! !DESCRIPTION: calculation of photolysis and heating rates
  ! 
  !     References:                                                        
  !       J. Landgraf and P.J. Crutzen, 1998:                               
  !       An Efficient Methode for Online Calculation of Photolysis and     
  !       Heating Rates, J. Atmos. Sci., 55, 863-878
  ! 
  !     Expanded for high angles and online:
  !       J.E.Williams, J. Landgraf, B. Bregman and H. H. Walter, A modified band
  !       approach for the accurate calculation of online photolysis rates in
  !       stratospheric-tropospheric Chemistry Transport Models, Atms. Chem. Phys., 
  !       6, 4137-4161, 2006.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE PHOTORATES_TROPO(region, zangle, cst_o3_col,   cst_no2_col, cst_hno3_col, cst_h2o2_col, &
                              cst_ch2o_col,   cst_n2o5_col, cst_pan_col, cst_no3_col, qy_no2_col, qy_o1d_col, &
                              qy_ch3cocho_col, qy_co_col,   qy_c2o3_col,  fact, fdir, rj, t_lay, debug_flag)
    !
    ! !USES:
    !
    use Dims,                only : im, jm, lm, idate
    use global_data,         only : mass_dat
    !
    ! !INPUT PARAMETERS:
    !
    integer,                                intent(in) :: region
    real,                                   intent(in) :: zangle
    real,dimension(lm(region),nbands_trop), intent(in) :: fact  !actinic flux
    real,dimension(lm(region),maxw),        intent(in) :: fdir  !flux o2-o3 abs.
    logical,                                intent(in) :: debug_flag
    real                                               :: tscale
    
    real,dimension(lm(region)),intent(in)              :: t_lay
    real,dimension(lm(region),maxwav),      intent(in) :: cst_o3_col
    real,dimension(lm(region),89),          intent(in) :: cst_no2_col,  qy_no2_col, qy_co_col, qy_c2o3_col
    real,dimension(lm(region),65),          intent(in) :: cst_hno3_col, qy_o1d_col
    real,dimension(lm(region),55),          intent(in) :: cst_n2o5_col
    real,dimension(lm(region),45),          intent(in) :: cst_h2o2_col
    real,dimension(lm(region),105),         intent(in) :: cst_ch2o_col
    real,dimension(lm(region),65),          intent(in) :: cst_pan_col
    real,dimension(lm(region),62),          intent(in) :: cst_no3_col
    real,dimension(lm(region),52),          intent(in) :: qy_ch3cocho_col   
    !
    ! !OUTPUT PARAMETERS:
    !
    real,dimension(lm(region),nj),         intent(out) :: rj  !photolysis rates 
    !
    ! !REVISION HISTORY: 
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    real,dimension(nbands_trop)       :: bjo1d,     bjno2,    bjhno3, bjcoh2,  bjchoh,   bjn2o5a,  bjn2o5b,   &
                                         bjhno4,    bjpana,   bjpanb,    bjno2o,   bjnoo2, bjh2o2,  bjch3ooh, bjald2,   bjch3o2co, &
                                         bjch3cocho,bjch3ono2,bjo2  , bjispd,  bj_a_acet,bj_b_acet,bjhono

    !=================================================================================================
    !First tropospheric band is temperature independent for H2O2.N2O5,HCHO and NO2 therefore remove from large temp dep. arrays.
    real,dimension(17)      :: cs_h2o2, cs_n2o5,cs_ch2o,cs_no2, cs_o2

    data cs_h2o2 /4.34E-19,  4.07E-19,  3.85E-19,  3.63E-19,  3.41E-19,  3.18E-19,  2.95E-19, &  
         2.72E-19,  2.50E-19,  2.30E-19,  2.10E-19,  1.92E-19,  1.74E-19,  1.57E-19, &  
         1.41E-19,  1.26E-19,  1.13E-19/ 

    data cs_n2o5 /8.049E-18, 7.208E-18, 6.212E-18, 4.853E-18, 3.925E-18, 3.248E-18, 2.619E-18, &
         2.087E-18, 1.661E-18, 1.349E-18, 1.131E-18, 9.549E-19, 8.353E-18, 7.427E-19, &
         6.762E-19, 6.095E-19, 5.229E-19/ 

    data  cs_ch2o /0.00E+00, 0.00E+00, 0.00E+00, 0.00E+00, 0.00E+00, 0.00E+00, 0.00E+00, &
         0.00E+00, 0.00E+00, 0.00E+00, 0.00E+00, 1.844E-22, 3.311E-22, 3.198E-22, &
         6.982E-22, 7.341E-22, 1.383E-21/

    data cs_o2 /7.448E-24, 7.322E-24, 7.002E-24, 6.608E-24, 6.118E-24, 5.736E-24, 5.302E-24,&
         4.741E-24, 4.261E-24, 3.788E-24, 3.213E-24, 2.69E-24, 2.218E-24, 1.793E-24,&
         1.384E-24, 1.054E-24, 6.318E-25/              

    real,dimension(62)      :: cs_hno4 ! 
    ! Reference : JPL 2006, 4-40 
    data cs_hno4 /4.43E-18,  3.64E-18,  3.09E-18,  2.54E-18,  2.13E-18,  1.78E-18,  1.51E-18, &  
         1.30E-18,  1.13E-18,  1.01E-18,  9.07E-19,  8.33E-19,  7.65E-19,  7.06E-19, &  
         6.48E-19,  5.91E-19,  5.36E-19,  4.83E-19,  4.36E-19,  3.93E-19,  3.53E-19, &   
         3.10E-19,  2.69E-19,  2.31E-19,  1.96E-19,  1.61E-19,  1.27E-19,  9.53E-20, &  
         7.01E-20,  4.93E-20,  3.31E-20,  2.14E-20,  1.60E-20,  1.40E-20,  1.30E-20, &  
         1.20E-20,  1.10E-20,  1.00E-20,  9.00E-21,  8.20E-21,  7.40E-21,  6.60E-21, &  
         5.80E-21,  5.00E-21,  4.60E-21,  4.20E-21,  3.80E-21,  3.40E-21,  3.00E-21, &  
         2.80E-21,  2.60E-21,  2.40E-21,  2.20E-21,  2.00E-21,  1.80E-21,  1.50E-21, &  
         1.10E-21,  7.00E-22,  0.00E+00,  0.00E+00,  0.00E+00,  0.00E+00/
         
    ! Reference : IUPAC, P21
    real,dimension(65)       :: qy_pan  
    data qy_pan /  7.589E-01, 7.589E-01, 7.589E-01, 7.589E-01, 7.589E-01, 7.589E-01, 7.589E-01, &
        7.589E-01, 7.589E-01, 7.589E-01, 7.589E-01, 7.589E-01, 7.589E-01, 7.589E-01, &
        7.589E-01, 7.589E-01, 7.589E-01, 7.589E-01, 7.589E-01, 7.511E-01, 7.431E-01, &
        7.348E-01, 7.264E-01, 7.177E-01, 7.089E-01, 6.997E-01, 6.903E-01, 6.807E-01, &
        6.708E-01, 6.606E-01, 6.501E-01, 6.393E-01, 6.325E-01, 6.300E-01, 6.275E-01, &
        6.250E-01, 6.225E-01, 6.200E-01, 6.175E-01, 6.150E-01, 6.125E-01, 6.100E-01, &
        5.995E-01, 5.890E-01, 5.785E-01, 5.680E-01, 5.575E-01, 5.470E-01, 5.365E-01, &
        5.260E-01, 5.155E-01, 5.050E-01, 4.945E-01, 4.840E-01, 4.735E-01, 4.577E-01, & 
        4.367E-01, 4.157E-01, 3.800E-01, 3.300E-01, 2.800E-01, 2.300E-01, 0.000E+00, &
        0.000E+00, 0.000E+00/     
    real, dimension(89)     :: cs_hono
    ! Reference : IUPAC 2001, datasheet PNOx1
    data cs_hono /2.110E-18, 2.198E-18, 2.173E-18, 2.147E-18, 2.025E-18, 1.867E-18, 1.710E-18, &
       1.554E-18, 1.408E-18, 1.280E-18, 1.133E-18, 9.571E-19, 8.147E-19, 7.137E-19, &
       6.104E-19, 5.046E-19, 3.962E-19, 2.908E-19, 2.207E-19, 1.685E-19, 1.348E-19, &
       1.003E-19, 7.195E-20, 5.256E-20, 3.956E-20, 3.020E-20, 2.069E-20, 1.399E-21, &
       0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 1.400E-21, &
       2.800E-21, 4.200E-21, 5.600E-21, 7.000E-21, 8.800E-21, 1.060E-20, 1.240E-20, &
       1.420E-20, 1.600E-20, 1.780E-20, 1.960E-20, 2.140E-20, 2.320E-20, 2.500E-20, &
       2.880E-20, 3.260E-20, 3.640E-20, 4.020E-20, 4.400E-20, 4.520E-20, 4.700E-20, &
       4.940E-20, 6.290E-20, 9.300E-20, 1.305E-19, 1.680E-19, 9.600E-20, 1.150E-19, &
       2.360E-19, 8.000E-20, 1.610E-19, 2.050E-19, 4.900E-20, 9.200E-20, 1.450E-19, &
       2.400E-20, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
       0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
       0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/

    real, dimension(89)     :: cs_ch3ooh
    ! Reference : JPL 2011, D9  
    data cs_ch3ooh / 2.782E-19, 2.782E-19, 2.782E-19, 2.782E-19, 2.782E-19, 2.316E-19, 1.956E-19, & 
         1.696E-19, 1.476E-19, 1.318E-19, 1.169E-19, 1.036E-19, 9.132E-20, 8.045E-20, &
         7.084E-20, 6.199E-20, 5.483E-20, 4.808E-20, 4.260E-20, 3.744E-20, 3.263E-20, & 
         2.819E-20, 2.431E-20, 2.119E-20, 1.827E-20, 1.569E-20, 1.338E-20, 1.107E-20, & 
         9.226E-21, 7.677E-21, 6.358E-21, 5.152E-21, 4.406E-21, 4.130E-21, 3.930E-21, & 
         3.730E-21, 3.530E-21, 3.330E-21, 3.130E-21, 2.982E-21, 2.834E-21, 2.686E-21, &  
         2.538E-21, 2.390E-21, 2.276E-21, 2.162E-21, 2.048E-21, 1.934E-21, 1.820E-21, & 
         1.730E-21, 1.640E-21, 1.550E-21, 1.460E-21, 1.370E-21, 1.306E-21, 1.210E-21, & 
         1.082E-21, 9.720E-22, 7.900E-22, 6.100E-22, 4.700E-22, 3.500E-22, 2.700E-22, & 
         2.100E-22, 1.600E-22, 1.200E-22, 7.500E-23, 5.200E-23, 4.000E-23, 4.050E-23, &
         4.100E-23, 1.000E-23, 6.000E-24, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
         0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
         0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/

    real, dimension(120)     :: cs_ch3cocho
    ! Reference : JPL 2006, P 4-71  
    data cs_ch3cocho /2.280E-19, 1.489E-19, 9.624E-20, 6.968E-20, 5.036E-20, 3.627E-20, 2.581E-20, & 
         1.729E-20, 1.440E-20, 1.435E-20, 1.460E-20, 1.521E-20, 1.619E-20, 1.743E-20, & 
         1.908E-20, 2.036E-20, 2.207E-20, 2.312E-20, 2.509E-20, 2.697E-20, 2.807E-20, & 
         3.175E-20, 3.343E-20, 3.580E-20, 4.070E-20, 4.234E-20, 4.473E-20, 4.906E-20, & 
         4.782E-20, 4.712E-20, 4.813E-20, 4.120E-20, 3.760E-20, 3.690E-20, 3.700E-20, & 
         3.740E-20, 3.740E-20, 3.620E-20, 3.380E-20, 3.150E-20, 2.920E-20, 2.710E-20, & 
         2.520E-20, 2.340E-20, 2.180E-20, 2.060E-20, 1.970E-20, 1.900E-20, 1.860E-20, & 
         1.860E-20, 1.870E-20, 1.820E-20, 1.680E-20, 1.500E-20, 1.340E-20, 1.180E-20, & 
         9.670E-21, 8.110E-21, 6.470E-21, 4.950E-21, 3.220E-21, 2.950E-21, 3.850E-21, & 
         5.560E-21, 7.650E-21, 1.080E-20, 1.470E-20, 1.900E-20, 2.420E-20, 3.200E-20, & 
         4.030E-20, 4.670E-20, 5.620E-20, 6.920E-20, 8.520E-20, 9.690E-20, 1.020E-19, & 
         1.030E-19, 1.040E-19, 1.080E-19, 9.940E-20, 9.620E-20, 8.680E-20, 3.690E-20, & 
         9.140E-21, 2.680E-21, 1.080E-21, 6.270E-22, 3.920E-22, 2.430E-22, 1.740E-22, & 
         0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
         0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
         0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
         0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, &
         0.000E+00/

    real,dimension(62)      :: cs_orgn
    ! Reference : average 1-C4H9ONO2 & 2-C4H9ONO2 (IUPAC data sheet P18 and P19)
    ! JEW: The absorption spectrum of 2-C4H9ONO2 is mirrored around the maxiumum.
    data cs_orgn / 5.250E-018, 4.398E-018, 3.609E-018, 2.804E-018, 2.172E-018, 1.595E-018, 1.130E-018, &
         7.710E-019, 5.025E-019, 3.714E-019, 2.623E-019, 1.900E-019, 1.342E-019, 9.905E-020, &
         7.285E-020, 5.245E-020, 4.052E-020, 3.110E-020, 2.806E-020, 5.649E-020, 5.136E-020, &
         4.886E-020, 4.635E-020, 4.358E-020, 3.970E-020, 3.610E-020, 3.247E-020, 2.784E-020, &
         2.308E-020, 1.846E-020, 1.401E-020, 9.810E-021, 7.430E-021, 6.550E-021, 6.080E-021, &
         5.610E-021, 5.140E-021, 4.670E-021, 4.200E-021, 3.840E-021, 3.480E-021, 3.120E-021, &
         2.760E-021, 2.400E-021, 2.170E-021, 1.940E-021, 1.710E-021, 1.480E-021, 1.250E-021, &
         1.131E-021, 1.012E-021, 8.930E-022, 7.740E-022, 2.550E-022, 2.350E-022, 2.050E-022, &
         1.650E-022, 1.400E-022, 1.050E-022, 8.000E-023, 0.000E+000, 0.000E+000/

    real, dimension(89)     :: cs_ald2         
    ! Reference : average CH3CHO & C2H5CHO (JPL and IUPAC data sheet, respectively)
    data cs_ald2 /5.211E-022, 5.133E-022, 5.163E-022, 5.272E-022, 5.526E-022, 5.837E-022, 6.266E-022, &
         6.774E-022, 7.498E-022, 8.806E-022, 1.054E-021, 1.387E-021, 1.849E-021, 2.472E-021, &
         3.444E-021, 4.679E-021, 6.217E-021, 8.229E-021, 1.074E-020, 1.376E-020, 1.729E-020, &
         2.131E-020, 2.584E-020, 3.082E-020, 3.565E-020, 4.055E-020, 4.482E-020, 4.809E-020, &
         5.184E-020, 5.155E-020, 5.245E-020, 4.795E-020, 4.640E-020, 4.600E-020, 4.540E-020, &
         4.465E-020, 4.330E-020, 4.085E-020, 3.870E-020, 3.730E-020, 3.585E-020, 3.490E-020, &
         3.380E-020, 3.265E-020, 3.145E-020, 3.015E-020, 2.895E-020, 2.750E-020, 2.485E-020, &
         2.235E-020, 2.125E-020, 1.990E-020, 1.869E-020, 1.777E-020, 1.631E-020, 1.471E-020, &
         1.254E-020, 1.014E-020, 6.315E-021, 3.375E-021, 1.525E-021, 2.300E-022, 9.000E-023, &
         3.000E-023, 1.500E-023, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, &
         0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, &
         0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, &
         0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000/
                 
    real, dimension(89)     :: cs_ispd
    ! Reference : average MACR and MVK (IUPAC data sheets)      1st two bins are estimates
    data cs_ispd /0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, &
                  0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, &
                  0.000E+000, 0.000E+000, 0.000E+000, 1.700E-021, 1.900E-021, 2.107E-021, 2.157E-021, &
                  2.403E-021, 2.928E-021, 3.713E-021, 4.854E-021, 6.439E-021, 8.562E-021, 1.147E-020, &
                  1.496E-020, 1.942E-020, 2.471E-020, 3.088E-020, 3.480E-020, 3.655E-020, 3.825E-020, &
                  3.980E-020, 4.130E-020, 4.275E-020, 4.425E-020, 4.610E-020, 4.770E-020, 4.920E-020, &
                  5.055E-020, 5.180E-020, 5.355E-020, 5.540E-020, 5.685E-020, 5.815E-020, 5.920E-020, &
                  6.075E-020, 6.230E-020, 6.365E-020, 6.455E-020, 6.485E-020, 6.470E-020, 6.558E-020, &
                  6.788E-020, 6.880E-020, 7.215E-020, 6.510E-020, 6.340E-020, 6.380E-020, 4.680E-020, &
                  4.145E-020, 4.365E-020, 2.680E-020, 1.650E-020, 1.144E-020, 1.175E-020, 5.440E-021, &
                  2.455E-021, 1.525E-021, 2.300E-022, 9.000E-023, 0.000E+000, 0.000E+000, 0.000E+000, &
                  0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, &
                  0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000/
    ! Reference : 5 abs. spectra at 5 different temperatures for acetone (235,254,263,280,298)
    real, dimension(77)     :: cs_ch3coch3_235
    data cs_ch3coch3_235 /1.185E-018, 9.785E-019, 7.682E-019, 5.534E-019, 3.341E-019, 1.102E-019, 1.875E-21, &
                          2.253E-021, 2.749E-021, 3.370E-021, 4.256E-021, 5.246E-021, 6.592E-021, 8.226E-21, &
                          1.025E-020, 1.277E-020, 1.562E-020, 1.898E-020, 2.278E-020, 2.694E-020, 3.129E-20, &
                          3.568E-020, 3.973E-020, 4.352E-020, 4.592E-020, 4.791E-020, 4.753E-020, 4.708E-20, &
                          4.391E-020, 4.018E-020, 3.524E-020, 2.906E-020, 2.510E-020, 2.370E-020, 2.280E-20, &
                          2.160E-020, 2.020E-020, 1.910E-020, 1.790E-020, 1.640E-020, 1.500E-020, 1.360E-20, &
                          1.250E-020, 1.130E-020, 1.020E-020, 9.320E-021, 8.620E-021, 7.630E-021, 6.710E-21, &
                          6.030E-021, 5.370E-021, 4.610E-021, 3.940E-021, 3.330E-021, 2.890E-021, 2.130E-21, &
                          1.390E-021, 8.720E-022, 3.550E-022, 1.200E-022, 5.110E-023, 1.590E-023, 0.000E+00, &
                          0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+00, &
                          0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+00/     
    real, dimension(77)     :: cs_ch3coch3_254
    data cs_ch3coch3_254 /1.177E-018, 9.719E-019, 7.630E-019, 5.497E-019, 3.319E-019, 1.095E-019, 1.892E-021, &
                          2.257E-021, 2.729E-021, 3.321E-021, 4.176E-021, 5.134E-021, 6.450E-021, 8.067E-021, &
                          1.005E-020, 1.252E-020, 1.532E-020, 1.868E-020, 2.248E-020, 2.660E-020, 3.089E-020, &
                          3.537E-020, 3.937E-020, 4.321E-020, 4.567E-020, 4.772E-020, 4.761E-020, 4.723E-020, &
                          4.428E-020, 4.068E-020, 3.584E-020, 2.986E-020, 2.590E-020, 2.450E-020, 2.360E-020, &
                          2.240E-020, 2.100E-020, 1.990E-020, 1.860E-020, 1.710E-020, 1.570E-020, 1.440E-020, &
                          1.320E-020, 1.200E-020, 1.090E-020, 9.940E-021, 9.190E-021, 8.150E-021, 7.190E-021, &
                          6.480E-021, 5.790E-021, 4.990E-021, 4.300E-021, 3.670E-021, 3.240E-021, 2.430E-021, &
                          1.630E-021, 1.056E-021, 4.500E-022, 1.450E-022, 5.400E-023, 1.810E-023, 0.000E+000, &
                          0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, &
                          0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000/    
    real, dimension(77)     :: cs_ch3coch3_263
    data cs_ch3coch3_263 /1.177E-018, 9.719E-019, 7.630E-019, 5.497E-019, 3.319E-019, 1.095E-019, 1.892E-021, & 
                          2.257E-021, 2.729E-021, 3.321E-021, 4.176E-021, 5.134E-021, 6.450E-021, 8.067E-021, & 
                          1.005E-020, 1.252E-020, 1.532E-020, 1.868E-020, 2.248E-020, 2.660E-020, 3.089E-020, & 
                          3.537E-020, 3.937E-020, 4.321E-020, 4.567E-020, 4.772E-020, 4.761E-020, 4.723E-020, & 
                          4.428E-020, 4.068E-020, 3.584E-020, 2.986E-020, 2.590E-020, 2.450E-020, 2.360E-020, &
                          2.240E-020, 2.100E-020, 1.990E-020, 1.860E-020, 1.710E-020, 1.570E-020, 1.440E-020, &
                          1.320E-020, 1.200E-020, 1.090E-020, 9.940E-021, 9.190E-021, 8.150E-021, 7.190E-021, &
                          6.480E-021, 5.790E-021, 4.990E-021, 4.300E-021, 3.670E-021, 3.240E-021, 2.430E-021, &
                          1.630E-021, 1.056E-021, 4.500E-022, 1.450E-022, 5.400E-023, 1.810E-023, 0.000E+000, &
                          0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, &
                          0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000/
    real, dimension(77)     :: cs_ch3coch3_280                       
    data cs_ch3coch3_280 /1.196E-018, 9.883E-019, 7.759E-019, 5.590E-019, 3.375E-019, 1.112E-019, 1.855E-021, &
                          2.233E-021, 2.719E-021, 3.340E-021, 4.225E-021, 5.200E-021, 6.540E-021, 8.166E-021, &
                          1.015E-020, 1.262E-020, 1.542E-020, 1.878E-020, 2.258E-020, 2.674E-020, 3.109E-020, &
                          3.558E-020, 3.967E-020, 4.361E-020, 4.627E-020, 4.843E-020, 4.851E-020, 4.823E-020, &
                          4.541E-020, 4.188E-020, 3.714E-020, 3.106E-020, 2.710E-020, 2.570E-020, 2.480E-020, &
                          2.350E-020, 2.200E-020, 2.080E-020, 1.960E-020, 1.800E-020, 1.660E-020, 1.530E-020, &
                          1.410E-020, 1.280E-020, 1.170E-020, 1.070E-020, 9.930E-021, 8.820E-021, 7.790E-021, &
                          7.040E-021, 6.310E-021, 5.480E-021, 4.760E-021, 4.100E-021, 3.650E-021, 2.795E-021, &
                          1.940E-021, 1.305E-021, 5.880E-022, 1.930E-022, 7.200E-023, 2.380E-023, 0.000E+000, &
                          0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, &
                          0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000/
    real, dimension(77)     :: cs_ch3coch3_298                       
    data cs_ch3coch3_298 /1.246E-019, 1.032E-019, 8.135E-020, 5.903E-020, 3.623E-020, 1.295E-020, 1.835E-021, &
                          2.213E-021, 2.699E-021, 3.310E-021, 4.186E-021, 5.166E-021, 6.490E-021, 8.097E-021, &
                          1.007E-020, 1.257E-020, 1.542E-020, 1.878E-020, 2.258E-020, 2.680E-020, 3.125E-020, &
                          3.578E-020, 3.997E-020, 4.401E-020, 4.677E-020, 4.913E-020, 4.931E-020, 4.913E-020, &
                          4.648E-020, 4.298E-020, 3.824E-020, 3.216E-020, 2.820E-020, 2.670E-020, 2.580E-020, &
                          2.450E-020, 2.300E-020, 2.180E-020, 2.050E-020, 1.890E-020, 1.750E-020, 1.610E-020, &
                          1.490E-020, 1.360E-020, 1.240E-020, 1.140E-020, 1.060E-020, 9.440E-021, 8.370E-021, &
                          7.600E-021, 6.840E-021, 5.980E-021, 5.230E-021, 4.550E-021, 4.110E-021, 3.210E-021, &
                          2.290E-021, 1.575E-021, 7.400E-022, 2.480E-022, 9.120E-023, 3.010E-023, 0.000E+000, &
                          0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, &
                          0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000, 0.000E+000/
                       
    real, dimension(77)     :: cs_ch3coch3                                       
    !==========================================================================================================
    ! Similar procedure for quantum yields for : HCHO and NO3. A quantum yield of 0.8 is adopted for JN2O5 below 305nm
    ! as recommended in JPL 2006 - JEW 2009        
    real,dimension(90)       :: qyh_hco,qyh2_co
    real,dimension(62)       :: qyno2_o,qyno_o2

    data qyh_hco  /0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,  & 
                   0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,  &
                   0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 3.087E-01, 3.087E-01, 3.030E-01,  &   
                   3.113E-01, 3.325E-01, 3.649E-01, 4.059E-01, 4.535E-01, 5.061E-01, 5.611E-01,  & 
                   6.159E-01, 6.662E-01, 7.097E-01, 7.418E-01, 7.550E-01, 7.580E-01, 7.610E-01,  & 
                   7.620E-01, 7.620E-01, 7.620E-01, 7.600E-01, 7.580E-01, 7.540E-01, 7.490E-01,  & 
                   7.440E-01, 7.370E-01, 7.290E-01, 7.200E-01, 7.090E-01, 6.980E-01, 6.850E-01,  & 
                   6.710E-01, 6.560E-01, 6.390E-01, 6.220E-01, 6.030E-01, 5.830E-01, 5.500E-01,  & 
                   5.020E-01, 4.490E-01, 3.430E-01, 1.650E-01, 0.000E+00, 0.000E+00, 0.000E+00,  &
                   0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,  & 
                   0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,  &
                   0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,  &  
                   0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/ !(JPL 2006, p4-47) 
                           
    data qyh2_co  /0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,   &   
                   0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,   &
                   0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 4.913E-01, 4.913E-01, 4.970E-01,   & 
                   4.887E-01, 4.697E-01, 4.546E-01, 4.313E-01, 4.020E-01, 3.682E-01, 3.440E-01,   & 
                   3.152E-01, 2.922E-01, 2.662E-01, 2.547E-01, 2.450E-01, 2.420E-01, 2.390E-01,   & 
                   2.380E-01, 2.380E-01, 2.380E-01, 2.400E-01, 2.420E-01, 2.460E-01, 2.510E-01,   & 
                   2.560E-01, 2.630E-01, 2.710E-01, 2.800E-01, 2.910E-01, 3.020E-01, 3.150E-01,   & 
                   3.290E-01, 3.440E-01, 3.610E-01, 3.780E-01, 3.970E-01, 4.170E-01, 4.500E-01,   &  
                   4.980E-01, 5.510E-01, 6.570E-01, 7.350E-01, 6.500E-01, 5.000E-01, 3.800E-01,   & 
                   2.200E-01, 4.000E-03, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,   &  
                   0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,   & 
                   0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,   &  
                   0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00/
                                
    data qyno2_o /0.9,  0.9,       0.9,      0.9,      0.9,     0.9,     0.9,               & 
                  0.9,  0.9,       0.9,      0.9,      0.9,     0.9,     0.9,               &
                  0.9,  0.9,       0.9,      0.9,      0.9,     0.9,     1.00,              &
                  1.00,  1.00,     1.00,     1.00,     1.00,    1.00,    1.00,              &
                  1.00,  1.00,     1.00,     1.00,     1.00,    1.00,    1.00,              & 
                  1.00,  1.00,     1.00,     1.00,     1.00,    1.00,    1.00,              &
                  1.00,  1.00,     1.00,     1.00,     1.00,    1.00,    1.00,              &
                  1.00,  0.83,     0.65,     0.58,     0.51,    0.43,    0.36,              &
                  0.29,  0.22,     0.14,     0.07,     0.0,     0.0/

    data qyno_o2 /0.0,  0.0,       0.0,      0.0,      0.0,     0.0,     0.0,               & 
                  0.0,  0.0,       0.0,      0.0,      0.0,     0.0,     0.0,               &
                  0.0,  0.0,       0.0,      0.0,      0.0,     0.0,     0.0,               &
                  0.0,  0.0,       0.0,      0.0,      0.0,     0.0,     0.0,               &
                  0.0,  0.0,       0.0,      0.0,      0.0,     0.0,     0.0,               & 
                  0.0,  0.0,       0.0,      0.0,      0.0,     0.0,     0.0,               &
                  0.0,  0.0,       0.0,      0.0,      0.0,     0.0,     0.0,               &
                  0.0,  0.18,      0.35,     0.31,     0.27,    0.23,    0.19,              &
                  0.16,  0.12,     0.08,     0.04,     0.0,     0.0/

    !========================================================================================================

    integer               :: i, j, l, k, c
    real                  :: delta1, delta2, delta3, delta4, delta5, delta6, delta7, delta1_spec, delta3_spec
    integer, dimension(7) :: band_start, band_middle, band_end
    logical               :: daylight

    ! ==============================================================        
    ! CALCULATION PHOTOLYSIS RATES FOR SCATTERING ATMOSPHERE             
    ! ============================================================== 

    ! the check for daylight is skipped as the calculation of rj is only called within the photo_flux loop
    ! where a filter already applies

    rj(:,:) = 0.
                       
    ! choose the grid settings dependent on the SZA for the column       
    if(zangle<=71.) then
       band_start(:)=lini(:)
       band_middle(:)=lmid(:)
       band_end(:)=lfin(:)
    elseif(zangle>71. .and. zangle<=sza_limit) then
       band_start(:)=lini_gridA(:)
       band_middle(:)=lmid_gridA(:)
       band_end(:)=lfin_gridA(:)
    endif

    do k = 1,lm(region)         
       ! re-initialize bj values               
       do l = 1,nbands_trop    
          bjo1d(l)      = 0.0     
          bjno2(l)      = 0.0     
          bjhno3(l)     = 0.0     
          bjcoh2(l)     = 0.0     
          bjchoh(l)     = 0.0
          bjn2o5a(l)    = 0.0
          bjn2o5b(l)    = 0.0     
          bjhno4(l)     = 0.0
          bjhono(l)     = 0.0    
          bjno2o(l)     = 0.0     
          bjnoo2(l)     = 0.0     
          bjh2o2(l)     = 0.0     
          bjch3ooh(l)   = 0.0  
          bjpana(l)     = 0.0
          bjpanb(l)     = 0.0     
          bjald2(l)     = 0.0
          bjch3cocho(l) = 0.0  
          bjch3ono2(l) = 0.0
          bjo2(l)      = 0.0
          bjispd(l)    = 0.0
          bj_a_acet(l) = 0.0
          bj_b_acet(l) = 0.0
       end do
       !                                      
       ! assign cs_ch3coch3  w/ a good value
       !
       ! linearly interpolate to try and improve on J values
       !
       if(t_lay(k)<=235.) cs_ch3coch3=cs_ch3coch3_235
       if(t_lay(k)>235. .and. t_lay(k)<=254.) then
          tscale=(t_lay(k)-235.)/19.0
          cs_ch3coch3=((1.0-tscale)*cs_ch3coch3_235)+(tscale*cs_ch3coch3_254)
       endif
       if(t_lay(k)>254. .and. t_lay(k)<=263.) then
          tscale=(t_lay(k)-254.)/9.0
          cs_ch3coch3=((1.0-tscale)*cs_ch3coch3_254)+(tscale*cs_ch3coch3_263)
       endif
       if(t_lay(k)>263. .and. t_lay(k)<=280.) then
          tscale=(t_lay(k)-264.)/16.0
          cs_ch3coch3=((1.0-tscale)*cs_ch3coch3_263)+(tscale*cs_ch3coch3_280)
       endif
       if(t_lay(k)>280.) then
          tscale=(t_lay(k)-280.)/18.0 
          cs_ch3coch3=((1.0-tscale)*cs_ch3coch3_280)+(tscale*cs_ch3coch3_298)
       endif
       !==============================================================================================
       ! re-initialize the temporary arrays for the next layer
       !=============================================================================================== 

       delta1 = 0. ; delta2 = 0. ;delta3 = 0. ; delta4 = 0. ; delta5 = 0. ; delta6 = 0. ; delta7 = 0.
       delta1_spec = 0. ; delta3_spec = 0.

       ! =====================================================================================
       !       Tropo band 1
       ! =====================================================================================
       ! temp indep species for this band: no2,h2o2,n2o5,hno4,ch2o,ch3ooh,ald2,orgntr
       do l = band_start(1),band_end(1) 
          bjo1d(1)  = bjo1d(1) + cst_o3_col(k,l)*qy_o1d_col(k,l)*fdir(k,l)
          bjno2(1)  = bjno2(1) + cst_no2_col(k,l)*qy_no2_col(k,l)*fdir(k,l)
          bjh2o2(1) = bjh2o2(1) + cs_h2o2(l)*fdir(k,l)
          bjhno3(1) = bjhno3(1) + cst_hno3_col(k,l)*fdir(k,l)
          bjhno4(1) = bjhno4(1) + cs_hno4(l)*fdir(k,l)
          bjhono(1) = bjhono(1) + cs_hono(l)*fdir(k,l)        
          bjn2o5a(1)= bjn2o5a(1) + 0.8*cs_n2o5(l)*fdir(k,l)
          bjn2o5b(1)= bjn2o5b(1) + 0.2*cs_n2o5(l)*fdir(k,l)
          bjchoh(1) = bjchoh(1) + cs_ch2o(l)*qyh_hco(l)*fdir(k,l)    
          bjch3ooh(1) = bjch3ooh(1) + cs_ch3ooh(l)*fdir(k,l)
          bjald2(1)   = bjald2(1) + cs_ald2(l)*fdir(k,l)
          bjpana(1)    = bjpana(1) + cst_pan_col(k,l)*qy_pan(l)*fdir(k,l)
          bjpanb(1)    = bjpanb(1) + cst_pan_col(k,l)*(1.0-qy_pan(l))*fdir(k,l)  
          bjch3ono2(1) = bjch3ono2(1) + cs_orgn(l) * fdir(k,l)
          bjo2(1) = bjo2(1) + cs_o2(l) * fdir(k,l)   
       end do

       ! ===================================================================================== 
       !           Tropo band 2                                                      
       ! =====================================================================================
       ! temp indep species for this band: hno4,ch3ooh,ald2,orgntr
       do l = band_start(2),band_end(2) ! temp indep for : no2,hno4,ch2o,ch3ooh,orgntr
          bjo1d(2)       = bjo1d(2) + cst_o3_col(k,l)*qy_o1d_col(k,l)*fdir(k,l)
          bjno2(2)       = bjno2(2) + cst_no2_col(k,l)*qy_no2_col(k,l)*fdir(k,l)
          bjh2o2(2)      = bjh2o2(2) + cst_h2o2_col(k,l-17)*fdir(k,l)
          bjhno3(2)      = bjhno3(2) + cst_hno3_col(k,l)*fdir(k,l)
          bjhno4(2)      = bjhno4(2) + cs_hno4(l)*fdir(k,l)
          bjhono(2)      = bjhono(2) + cs_hono(l)*fdir(k,l)
          bjn2o5a(2)     = bjn2o5a(2) + 0.8*cst_n2o5_col(k,l-17)*fdir(k,l)
          bjn2o5b(2)     = bjn2o5b(2) + 0.2*cst_n2o5_col(k,l-17)*fdir(k,l)
          bjcoh2(2)      = bjcoh2(2) + cst_ch2o_col(k,l-17)*qyh2_co(l)*fdir(k,l)    
          bjchoh(2)      = bjchoh(2) + cst_ch2o_col(k,l-17)*qyh_hco(l)*fdir(k,l)                             
          bjch3ooh(2)    = bjch3ooh(5) + cs_ch3ooh(l)*fdir(k,l)
          bjpana(2)      = bjpana(2) + cst_pan_col(k,l)*qy_pan(l)*fdir(k,l)
          bjpanb(2)      = bjpanb(2) + cst_pan_col(k,l)*(1.0-qy_pan(l))*fdir(k,l)
          bjald2(2)      = bjald2(2) + cs_ald2(l)*fdir(k,l)
          bjch3cocho(2)  = bjch3cocho(2) + cs_ch3cocho(l)*fdir(k,l)
          bjch3ono2(2)   = bjch3ono2(2) + cs_orgn(l) * fdir(k,l)
          bjispd(2)      = bjispd(2) + cs_ispd(l) * fdir(k,l)
          bj_a_acet(2)   = bj_a_acet(2) + cs_ch3coch3(l)*qy_co_col(k,l)*fdir(k,l)              
          bj_b_acet(2)   = bj_b_acet(2) + cs_ch3coch3(l)*qy_c2o3_col(k,l)*fdir(k,l)       
       end do

       ! ======================================================================================= 
       !           Tropo band 3                                                      
       ! ======================================================================================= 
       ! temp indep species for this band: hno4,ch3ooh,ald2,mgly,orgntr      
       do l = band_start(3),band_end(3)
          bjo1d(3)       = bjo1d(3) + cst_o3_col(k,l)*qy_o1d_col(k,l)*fdir(k,l)
          bjno2(3)       = bjno2(3) + cst_no2_col(k,l)*qy_no2_col(k,l)*fdir(k,l)
          bjh2o2(3)      = bjh2o2(3) + cst_h2o2_col(k,l-17)*fdir(k,l)
          bjhno3(3)      = bjhno3(3) + cst_hno3_col(k,l)*fdir(k,l)
          bjhno4(3)      = bjhno4(3) + cs_hno4(l)*fdir(k,l)
          bjhono(3)      = bjhono(3) + cs_hono(l)*fdir(k,l)
          bjn2o5a(3)     = bjn2o5a(3) + 0.9*cst_n2o5_col(k,l-17)*fdir(k,l)
          bjn2o5b(3)     = bjn2o5b(3) + 0.1*cst_n2o5_col(k,l-17)*fdir(k,l)
          bjcoh2(3)      = bjcoh2(3) + cst_ch2o_col(k,l-17)*qyh2_co(l)*fdir(k,l)                      
          bjchoh(3)      = bjchoh(3) + cst_ch2o_col(k,l-17)*qyh_hco(l)*fdir(k,l)                                 
          bjch3ooh(3)    = bjch3ooh(3) + cs_ch3ooh(l)*fdir(k,l)
          bjpana(3)      = bjpana(3) + cst_pan_col(k,l)*qy_pan(l)*fdir(k,l)
          bjpanb(3)      = bjpanb(3) + cst_pan_col(k,l)*(1.0-qy_pan(l))*fdir(k,l)
          bjald2(3)      = bjald2(3) + cs_ald2(l)*fdir(k,l)    
          bjch3cocho(3)  = bjch3cocho(3) + cs_ch3cocho(l)*fdir(k,l)
          bjch3ono2(3)   = bjch3ono2(3) + cs_orgn(l)*fdir(k,l)
          bjispd(3)      = bjispd(3) + cs_ispd(l) * fdir(k,l)
          bj_a_acet(3)   = bj_a_acet(3) + cs_ch3coch3(l)*qy_co_col(k,l)*fdir(k,l)
          bj_b_acet(3)   = bj_b_acet(3) + cs_ch3coch3(l)*qy_c2o3_col(k,l)*fdir(k,l)                      
       end do

       ! ================================================================================ 
       !           Tropo band 4                                                      
       ! ================================================================================
       ! temp indep species for this band: hno4,ch3ooh,ald2,orgntr  
       do l = band_start(4),band_end(4)
          bjo1d(4)      = bjo1d(4)+cst_o3_col(k,l)*qy_o1d_col(k,l)*fdir(k,l)
          bjno2(4)      = bjno2(4) + cst_no2_col(k,l)*qy_no2_col(k,l)*fdir(k,l)
          bjh2o2(4)     = bjh2o2(4) + cst_h2o2_col(k,l-17)*fdir(k,l)
          bjhno3(4)     = bjhno3(4) + cst_hno3_col(k,l)*fdir(k,l)
          bjhno4(4)     = bjhno4(4) + cs_hno4(l)*fdir(k,l)
          bjhono(4)     = bjhono(4) + cs_hono(l)*fdir(k,l)
          bjn2o5a(4)    = bjn2o5a(4) + cst_n2o5_col(k,l-17)*fdir(k,l)
          bjcoh2(4)     = bjcoh2(4) + cst_ch2o_col(k,l-17)*qyh2_co(l)*fdir(k,l)
          bjchoh(4)     = bjchoh(4) + cst_ch2o_col(k,l-17)*qyh_hco(l)*fdir(k,l)           
          bjch3ooh(4)   = bjch3ooh(4) + cs_ch3ooh(l)*fdir(k,l)
          bjpana(4)     = bjpana(4) + cst_pan_col(k,l)*qy_pan(l)*fdir(k,l)
          bjpanb(4)     = bjpanb(4) + cst_pan_col(k,l)*(1.0-qy_pan(l))*fdir(k,l)
          bjald2(4)     = bjald2(4) + cs_ald2(l)*fdir(k,l)        
          bjch3cocho(4) = bjch3cocho(4) + cs_ch3cocho(l)*qy_ch3cocho_col(k,l-37)*fdir(k,l)
          bjch3ono2(4)  = bjch3ono2(4) + cs_orgn(l)*fdir(k,l)
          bjispd(4)     = bjispd(4) + cs_ispd(l) * fdir(k,l)
          bj_a_acet(4)  = bj_a_acet(4) + cs_ch3coch3(l)*qy_co_col(k,l)*fdir(k,l)
          bj_b_acet(4)  = bj_b_acet(4) + cs_ch3coch3(l)*qy_c2o3_col(k,l)*fdir(k,l)        
       end do

       ! ======================================================================================== 
       !           Tropo band 5                                                      
       ! ======================================================================================== 
       ! temp indep species for this band: hno4,ch3ooh,ald2,orgntr
       do l = band_start(5),band_end(5)
          bjo1d(5)      = bjo1d(5)+cst_o3_col(k,l)*qy_o1d_col(k,l)*fdir(k,l)
          bjno2(5)      = bjno2(5)+cst_no2_col(k,l)*qy_no2_col(k,l)*fdir(k,l)
          bjh2o2(5)     = bjh2o2(5)+cst_h2o2_col(k,l-17)*fdir(k,l)
          bjhno3(5)     = bjhno3(5)+cst_hno3_col(k,l)*fdir(k,l)
          bjhno4(5)     = bjhno4(5)+cs_hno4(l)*fdir(k,l)
          bjhono(5)     = bjhono(5)+cs_hono(l)*fdir(k,l)
          bjn2o5a(5)    = bjn2o5a(5)+cst_n2o5_col(k,l-17)*fdir(k,l)
          bjcoh2(5)     = bjcoh2(5)+cst_ch2o_col(k,l-17)*qyh2_co(l)*fdir(k,l)
          bjchoh(5)     = bjchoh(5)+cst_ch2o_col(k,l-17)*qyh_hco(l)*fdir(k,l)
          bjch3ooh(5)   = bjch3ooh(5)+cs_ch3ooh(l)*fdir(k,l)
          bjpana(5)     = bjpana(5) + cst_pan_col(k,l)*qy_pan(l)*fdir(k,l)
          bjpanb(5)     = bjpanb(5) + cst_pan_col(k,l)*(1.0-qy_pan(l))*fdir(k,l)
          bjald2(5)     = bjald2(5) + cs_ald2(l)*fdir(k,l)        
          bjch3cocho(5) = bjch3cocho(5)+cs_ch3cocho(l)*qy_ch3cocho_col(k,l-37)*fdir(k,l)
          bjch3ono2(5)  = bjch3ono2(5)+cs_orgn(l)*fdir(k,l)
          bjispd(5)     = bjispd(5) + cs_ispd(l) * fdir(k,l)      
       end do
       !=========================================================================================
       !           Tropo band 6
       !=========================================================================================

       do l = band_start(6),band_end(6)         
          bjno2(6)      = bjno2(6) +cst_no2_col(k,l)*qy_no2_col(k,l)*fdir(k,l)
          bjcoh2(6)     = bjcoh2(6)+cst_ch2o_col(k,l-17)*qyh2_co(l)*fdir(k,l)
          bjch3ooh(6)   = bjch3ooh(6)+cs_ch3ooh(l)*fdir(k,l)
          bjno2o(6)     = bjno2o(6)+cst_no3_col(k,l-60)*qyno2_o(l-60)*fdir(k,l)
          bjch3cocho(6) = bjch3cocho(6)+cs_ch3cocho(l)*qy_ch3cocho_col(k,l-37)*fdir(k,l)
          bjispd(6)     = bjispd(6) + cs_ispd(l) * fdir(k,l)      
       end do

       ! Add the contribution from band 6 separately for H2O2, N2O5 and HNO3 for high angle grid to allow reducion in cst arrays
       ! This has been tested using a box model to ensure that no errors are introduced compared to using main array

       if(zangle>71.) then

          bjh2o2(6)     = bjh2o2(6)+cst_h2o2_col(k,44)*fdir(k,61)
          bjh2o2(6)     = bjh2o2(6)+cst_h2o2_col(k,45)*fdir(k,62)

          do c=63,68
             bjn2o5a(6)    = bjn2o5a(6)+cst_n2o5_col(k,c-17)*fdir(k,c)
          enddo

          bjhno3(6)     = bjhno3(6)+cst_hno3_col(k,63)*fdir(k,63)

          bjpana(6)      = 0.
          bjpanb(6)      = 0.

       elseif(zangle<=71.) then

          do c=61,63
             bjhno3(6)     = bjhno3(6)+cst_hno3_col(k,c)*fdir(k,c)
          enddo

          do c=61,69
             bjn2o5a(6)    = bjn2o5a(6)+cst_n2o5_col(k,c-17)*fdir(k,c)
          enddo

          do c=61,62
             bjpana(6)        = bjpana(6)+cst_pan_col(k,c)*qy_pan(c)*fdir(k,c)
             bjpanb(6)        = bjpanb(6)+cst_pan_col(k,c)*(1.0-qy_pan(c))*fdir(k,c)
          enddo

       endif
       !======================================================================================== 
       !           Tropo band 7                                                      
       !========================================================================================= 

       do l = band_start(7),band_end(7)-2
          bjno2o(7)     = bjno2o(7)+cst_no3_col(k,l-60)*qyno2_o(l-60)*fdir(k,l)
          bjnoo2(7)     = bjnoo2(7)+cst_no3_col(k,l-60)*qyno_o2(l-60)*fdir(k,l)
       end do

       ! only set a scaling ratio for the first band once the sza > 71.
       ! only calculate limits to bands 2 and 4 for high sza         

       if(zangle<=71.) then
          delta1 = fact(k,1)/fdir(k,band_middle(1))
          delta2 = fact(k,2)/fdir(k,band_middle(2))
          delta3 = fact(k,3)/fdir(k,band_middle(3))
          delta4 = fact(k,4)/fdir(k,band_middle(4))
          delta5 = fact(k,5)/fdir(k,band_middle(5))
          delta6 = fact(k,6)/fdir(k,band_middle(6))
          delta7 = fact(k,7)/fdir(k,band_middle(7))
       elseif(zangle>71. .and. zangle<=sza_limit) then
          delta1 = safe_div( fact(k,1), fdir(k,band_middle(1)), huge(1.))
          delta2 = safe_div( fact(k,2), fdir(k,band_middle(2)), huge(1.))
          delta3 = safe_div( fact(k,3), fdir(k,band_middle(3)), huge(1.))
          delta4 = safe_div( fact(k,4), fdir(k,band_middle(4)), huge(1.))
          delta5 = safe_div( fact(k,5), fdir(k,band_middle(5)), huge(1.))
          delta6 = safe_div( fact(k,6), fdir(k,band_middle(6)), huge(1.))
          delta7 = safe_div( fact(k,7), fdir(k,band_middle(7)), huge(1.))
          if(fdir(k,band_middle(1))>5.e9) delta1_spec = fact(k,1)/fdir(k,band_middle(1))           
          if(fdir(k,band_middle(3))>1.e9) delta3_spec = fact(k,3)/fdir(k,band_middle(3))
       endif

       !        ============================================================
       !        CALCULATION PHOTOLYSIS RATES FOR SCATTERING ATMOSPHERE    
       !        ============================================================

       ! O2 -> O3P

       rj(k,jo2)   = bjo2(1)*delta1

       ! O3 -> 2OH 

       rj(k,jo3d)  =  bjo1d(1)*delta1 + bjo1d(2)*delta2 &
            + bjo1d(3)*delta3 + bjo1d(4)*delta4 &
            + bjo1d(5)*delta5

       ! NO2 -> O3

       rj(k,jno2)  = bjno2(1)*delta1 + bjno2(2)*delta2 + &
            bjno2(3)*delta3 + bjno2(4)*delta4 + &
            bjno2(5)*delta5 + bjno2(6)*delta6   

       ! HNO3 -> OH + NO2 

       rj(k,jhno3)  = bjhno3(1)*delta1 + bjhno3(2)*delta2 + &
            bjhno3(3)*delta3 + bjhno3(4)*delta4 + &
            bjhno3(5)*delta5 + bjhno3(6)*delta6

       ! HNO4 -> HO2 + NO2 

       rj(k,jhno4)  =  bjhno4(1)*delta1 + bjhno4(2)*delta2 + &
            bjhno4(3)*delta3 + bjhno4(4)*delta4 + &
            bjhno4(5)*delta5 
            
       ! HONO -> OH + NO
       
       rj(k,jhono)  =  bjhono(1)*delta1 + bjhono(2)*delta2 + &
            bjhono(3)*delta3 + bjhono(4)*delta4 + &
            bjhono(5)*delta5        

       ! N2O5 -> NO3 + NO2                                  
       rj(k,jn2o5a)  =  bjn2o5a(1)*delta1 + bjn2o5a(2)*delta2 +  &
            bjn2o5a(3)*delta3 + bjn2o5a(4)*delta4 +  &
            bjn2o5a(5)*delta5 + bjn2o5a(6)*delta6 
            
       ! N2O5 -> NO3 + NO                                  
       rj(k,jn2o5b)  =  bjn2o5b(1)*delta1 + bjn2o5b(2)*delta2 +  &
            bjn2o5b(3)*delta3       

       ! PAN   -> NO2 + C2O3

       rj(k,jpana) =   bjpana(1)*delta1  + bjpana(2)*delta2 +   &
            bjpana(3)*delta3  + bjpana(4)*delta4 +   &
            bjpana(5)*delta5  + bjpana(6)*delta6
            
       ! PAN   -> NO3 + CH3O2
       
        rj(k,jpanb) =   bjpanb(1)*delta1  + bjpanb(2)*delta2 +   &
            bjpanb(3)*delta3  + bjpanb(4)*delta4 +   &
            bjpanb(5)*delta5  + bjpanb(6)*delta6            

       ! CH2O -> CO

       rj(k,jach2o)  =  bjcoh2(2)*delta2  +  bjcoh2(3)*delta3 &
            + bjcoh2(4)*delta4 +  bjcoh2(5)*delta5 &
            + bjcoh2(6)*delta6   

       ! CH2O -> COH + H

       rj(k,jbch2o) = bjchoh(1)*delta1 + bjchoh(2)*delta2 &
                    + bjchoh(3)*delta3 + bjchoh(4)*delta4 &
                    + bjchoh(5)*delta5 + bjchoh(6)*delta6

       ! H2O2 -> 2OH

       rj(k,jh2o2) = bjh2o2(1)*delta1 + bjh2o2(2)*delta2 &
                   + bjh2o2(3)*delta3 + bjh2o2(4)*delta4 &
                   + bjh2o2(5)*delta5 + bjh2o2(6)*delta6 

       ! NO3 -> NO2 + O

       rj(k,jano3)  = bjno2o(6)*delta6 + bjno2o(7)*delta7

       ! NO3 -> NO

       rj(k,jbno3) = bjnoo2(7)*delta7

       ! CH3OOH -> HCHO + HO2 + OH                                                                        
       rj(k,jmepe) = bjch3ooh(1)*delta1 + bjch3ooh(2)*delta2 + &
                     bjch3ooh(3)*delta3 + bjch3ooh(4)*delta4 + &
                     bjch3ooh(5)*delta5 + bjch3ooh(6)*delta6 

         ! ALD2 -> 2HO2 + CO + CH3O2 
         
       rj(k,jald2) = bjald2(1)*delta1 + bjald2(2)*delta2     + &
                     bjald2(3)*delta3 + bjald2(4)*delta4     + &
                     bjald2(5)*delta5
                 

       ! CH3ONO2 -> HO2 + NO2

       rj(k,jorgn) = bjch3ono2(1)*delta1 + bjch3ono2(2)*delta2 + &
                     bjch3ono2(3)*delta3 + bjch3ono2(4)*delta4 + &
                     bjch3ono2(5)*delta5


       ! CH3COCHO -> C2O3 + HO2 + CO                              
       rj(k,jmgly) = bjch3cocho(2)*delta2 + bjch3cocho(3)*delta3 +  &
                   bjch3cocho(4)*delta4 + bjch3cocho(5)*delta5 + &
                   bjch3cocho(6)*delta6 + bjch3cocho(7)*delta7 
            
       ! ISPD -> 0.333CO + 0.067ALD2 + 0.9CH2O + 0.832PAR + 1.033HO2 + 0.7XO2 + 0.967C2O3
          
       rj(k,jispd) = bjispd(2)*delta2 + bjispd(3)*delta3 + &
                     bjispd(4)*delta4 + bjispd(5)*delta5 + &
                     bjispd(6)*delta6
                       
       !ch3coch3 -> 2CH3O2 + CO
           
       rj(k,ja_acet) =  bj_a_acet(2)*delta2 + bj_a_acet(3)*delta3 + bj_a_acet(4)*delta4
       
       rj(k,jb_acet) = bj_b_acet(2)*delta2 + bj_b_acet(3)*delta3 + bj_b_acet(4)*delta4 

       if(zangle>80. .and. zangle<=85.) then
          
          rj(k,jo3d) = 0. ; rj(k,jhno3) = 0. ; rj(k,jhno4) = 0. ; rj(k,jn2o5a) = 0. ; rj(k,jn2o5b) = 0.
          rj(k,jhono) = 0. ; rj(k,jh2o2) = 0. ; rj(k,jach2o) = 0. ; rj(k,jpana) = 0. ; rj(k,jpanb) = 0.
          rj(k,jald2) = 0. ; rj(k,jorgn) = 0. ; rj(k,jbch2o) = 0. ; rj(k,jmgly) = 0.
          rj(k,jo2) = 0. ; rj(k,jispd) = 0. ; rj(k,ja_acet) = 0. ; rj(k,jb_acet) = 0.

          ! O2 -> 2O3P

          rj(k,jo2) = bjo2(1)*delta1_spec

          ! O3 -> O(1D) + O2                                      
          rj(k,jo3d)  =  bjo1d(1)*delta1_spec + bjo1d(2)*delta2 +   &
                         bjo1d(3)*delta3_spec + bjo1d(4)*delta4 +   &
                         bjo1d(5)*delta5

          ! HNO3 -> OH + NO2                                                                      
          rj(k,jhno3)  = bjhno3(1)*delta1_spec + bjhno3(2)*delta2 + &
                         bjhno3(3)*delta3_spec + bjhno3(4)*delta4 + &
                         bjhno3(5)*delta5      + bjhno3(6)*delta6                 

          ! N2O5 -> NO3 + NO2                                  
          rj(k,jn2o5a)  =  bjn2o5a(1)*delta1_spec + bjn2o5a(2)*delta2 +  &
               bjn2o5a(3)*delta3_spec + bjn2o5a(4)*delta4 +  &
               bjn2o5a(5)*delta5  + bjn2o5a(6)*delta6
               
          ! N2O5 -> NO3 + NO                                  
          rj(k,jn2o5b)  =  bjn2o5b(1)*delta1_spec + bjn2o5b(2)*delta2 +  &
               bjn2o5b(3)*delta3_spec        

          ! PAN 

          rj(k,jpana) = bjpana(1)*delta1_spec + bjpana(2)*delta2 +   &
               bjpana(3)*delta3_spec + bjpana(4)*delta4 +   &
               bjpana(5)*delta5  + bjpana(6)*delta6 
               
          rj(k,jpanb) = bjpanb(1)*delta1_spec + bjpanb(2)*delta2 +   &
               bjpanb(3)*delta3_spec + bjpanb(4)*delta4 +   &
               bjpanb(5)*delta5  + bjpanb(6)*delta6             

          ! ORGNTR -> HO2 + NO2

          rj(k,jorgn) = bjch3ono2(1)*delta1_spec + bjch3ono2(2)*delta2 +   &
               bjch3ono2(3)*delta3_spec + bjch3ono2(4)*delta4 +   &
               bjch3ono2(5)*delta5 

         ! ALD2 -> 2HO2 + CO + CH3O2 

          rj(k,jald2) = bjald2(1)*delta1_spec + bjald2(2)*delta2 + &
                        bjald2(3)*delta3_spec + bjald2(4)*delta4 + &
                        bjald2(5)*delta5
                            
          ! NO2 -> O3

          rj(k,jno2)  =  bjno2(1)*delta1_spec + bjno2(2)*delta2 + &
               bjno2(3)*delta3_spec + bjno2(4)*delta4 + &
               bjno2(5)*delta5 + bjno2(6)*delta6

          ! CH2O -> CO

          rj(k,jach2o)  =  bjcoh2(2)*delta2 + bjcoh2(3)*delta3_spec + &
               bjcoh2(4)*delta4 + bjcoh2(5)*delta5 + &
               bjcoh2(6)*delta6

          ! CH2O -> COH + H

          rj(k,jbch2o)  =  bjchoh(1)*delta1_spec + bjchoh(2)*delta2 + &
               bjchoh(3)*delta3_spec + bjchoh(4)*delta4 + &
               bjchoh(5)*delta5 + bjchoh(6)*delta6

          rj(k,jhno4)  =  bjhno4(1)*delta1_spec + bjhno4(2)*delta2 + &
               bjhno4(3)*delta3_spec + bjhno4(4)*delta4 + &
               bjhno4(5)*delta5
               
          rj(k,jhono)  =  bjhono(1)*delta1_spec + bjhono(2)*delta2 + &
               bjhono(3)*delta3_spec + bjhono(4)*delta4 + &
               bjhono(5)*delta5        

          ! H2O2 -> 2OH                                                                           
          rj(k,jh2o2) = bjh2o2(1)*delta1_spec + bjh2o2(2)*delta2 + &
               bjh2o2(3)*delta3_spec + bjh2o2(4)*delta4 + &
               bjh2o2(5)*delta5 + bjh2o2(6)*delta6

          ! CH3OOH -> HCHO + HO2 + OH                                                                     
          rj(k,jmepe) = bjch3ooh(1)*delta1_spec + bjch3ooh(2)*delta2 + &
               bjch3ooh(3)*delta3_spec + bjch3ooh(4)*delta4 + &
               bjch3ooh(5)*delta5 + bjch3ooh(6)*delta6 

          ! CH3COCHO -> C2O3 + HO2 + CO                                   
          rj(k,jmgly) = bjch3cocho(2)*delta2 + bjch3cocho(3)*delta3_spec +  &
               bjch3cocho(4)*delta4 + bjch3cocho(5)*delta5 +  &
               bjch3cocho(6)*delta6 + bjch3cocho(7)*delta7
                      
          ! ispd -> 0.333CO + 0.067ALD2 + 0.9CH2O + 0.832PAR + 1.033HO2 + 0.7XO2 + 0.967C2O3
          
          rj(k,jispd) = bjispd(2)*delta2 + bjispd(3)*delta3_spec + &
                        bjispd(4)*delta4 + bjispd(5)*delta5 + &
                        bjispd(6)*delta6
                       
          rj(k,ja_acet) = bj_a_acet(2)*delta2 + bj_a_acet(3)*delta3_spec + bj_a_acet(4)*delta4 
                         
          rj(k,jb_acet) = bj_b_acet(2)*delta2 + bj_b_acet(3)*delta3_spec + bj_b_acet(4)*delta4
                
       endif ! sza > 80. 

       rj(k,jrooh) = rj(k,jmepe)
       !
       !  take the upper limit for the photolysis rate i.e. qy=0.05
       !                
       rj(k,jispd) = 0.05*rj(k,jispd)
       !
       ! Two channels of CH3O2NO2 photo-dissociation equal to HNO4 photolysis
       ! Taken from Browne et al, ACP, 2011 
       ! 
        
       rj(k,jmo2no2a) = rj(k,jhno4)
       rj(k,jmo2no2b) = rj(k,jhno4)       
    enddo ! layers          
    

  END SUBROUTINE PHOTORATES_TROPO
  !EOC

  !------------------------------------------------------------------------------
  !                    TM5                                                      !
  !------------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    PHOTOLYSIS_DONE
  !
  ! !DESCRIPTION: Deallocates photolysis data, and writes photolysis statistics:
  !               output (e.g. monthly) averaged photolysis jo3, jno2 and
  !               surface data of: albedo, cloud-cover, cloud optical thickness 
  !               total ozone
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE PHOTOLYSIS_DONE( status )
    !
    ! !USES:
    !
    use dims,            only : im, jm, lm, nregions
    use dims,            only : idatei, idatee
    use global_data,     only : outdir
#ifdef with_hdf4
    use file_hdf,        only : THdfFile, TSds
    use file_hdf,        only : Init, Done, WriteAttribute, Init, Done, WriteData
#endif
    use partools,        only : isRoot
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)       ::  status
    !
    ! !REVISION HISTORY: 
    !   21 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:
    !   (1) Called from SOURCES_SINKS/TRACE_END
    !
    !EOP
    !------------------------------------------------------------------------------
    !BOC

    character(len=*), parameter  ::  rname = mname//'/end_photolysis'

    integer              ::  region
    integer              ::  n 
    integer              ::  k,iret,imr,jmr,lmr
    real, allocatable    ::  average_field(:,:,:), average_field_2d(:,:), average_field_5d(:,:,:,:,:)
#ifdef with_hdf4
    type(THdfFile)       ::  io
    type(TSds)           ::  sds
#endif
    character(len=256)  ::  fname

    ! --- begin --------------------------------

    if (isRoot) then
       write (fname,'(a,"/j_statistics_",i4.4,3i2.2,"_",i4.4,3i2.2,".hdf")') &
            trim(outdir), idatei(1:4), idatee(1:4)
#ifdef with_hdf4
       call Init(io,fname,'create',status)
       IF_ERROR_RETURN(status=1)
#endif
    end if

    
    REG: do region=1,nregions
    
       if(isRoot) then
          imr=im(region)
          jmr=jm(region)
          lmr=lm(region)
       else
          imr=1; jmr=1; lmr=1
       end if

       allocate(average_field(imr,jmr,lmr))

       ! Assume that the nav is the same on all proc
       write (gol,*) 'end_photolysis: j_statistics ', &
            '(region, nj_av, nvo3s_av, ncloud_av)', region,   &
            phot_dat(region)%nj_av, phot_dat(region)%nvo3s_av; call goPr

#ifdef with_hdf4
       !---- 3D fields
       if ( phot_dat(region)%nj_av > 0 ) then

         call gather(dgrid(region), phot_dat(region)%jo3_av, average_field, 0, status)

          if ( isRoot ) then
             average_field = average_field/phot_dat(region)%nj_av
             call Init(Sds,io,'jo3_av',(/imr,jmr,lmr/),'real(4)',status)
             IF_ERROR_RETURN(status=1)
             call WriteAttribute(Sds,'Unit','1/s',status)
             IF_ERROR_RETURN(status=1)
             call WriteData(Sds,average_field,status)
             IF_ERROR_RETURN(status=1)
             call Done(Sds,status)
             IF_ERROR_RETURN(status=1)
          end if

          call gather(dgrid(region), phot_dat(region)%jno2_av, average_field, 0, status)

          if ( isRoot ) then
             average_field = average_field/phot_dat(region)%nj_av
             call Init(Sds,io,'jno2_av',(/imr,jmr,lmr/),'real(4)',status)
             IF_ERROR_RETURN(status=1)
             call WriteAttribute(Sds,'Unit','1/s',status)
             IF_ERROR_RETURN(status=1)
             call WriteData(Sds,average_field,status)
             IF_ERROR_RETURN(status=1)
             call Done(Sds,status)
             IF_ERROR_RETURN(status=1)
          end if
           
          call gather(dgrid(region), phot_dat(region)%jhno2_av, average_field, 0, status)

          if ( isRoot ) then
             average_field = average_field/phot_dat(region)%nj_av
             call Init(Sds,io,'jhno2_av',(/imr,jmr,lmr/),'real(4)',status)
             IF_ERROR_RETURN(status=1)
             call WriteAttribute(Sds,'Unit','1/s',status)
             IF_ERROR_RETURN(status=1)
             call WriteData(Sds,average_field,status)
             IF_ERROR_RETURN(status=1)
             call Done(Sds,status)
             IF_ERROR_RETURN(status=1)
          end if
          
          call gather(dgrid(region), phot_dat(region)%jch2oa_av, average_field, 0, status)

          if ( isRoot ) then
             average_field = average_field/phot_dat(region)%nj_av
             call Init(Sds,io,'jch2oa_av',(/imr,jmr,lmr/),'real(4)',status)
             IF_ERROR_RETURN(status=1)
             call WriteAttribute(Sds,'Unit','1/s',status)
             IF_ERROR_RETURN(status=1)
             call WriteData(Sds,average_field,status)
             IF_ERROR_RETURN(status=1)
             call Done(Sds,status)
             IF_ERROR_RETURN(status=1)
          end if
          
          call gather(dgrid(region), phot_dat(region)%jch2ob_av, average_field, 0, status)

          if ( isRoot ) then
             average_field = average_field/phot_dat(region)%nj_av
             call Init(Sds,io,'jch2ob_av',(/imr,jmr,lmr/),'real(4)',status)
             IF_ERROR_RETURN(status=1)
             call WriteAttribute(Sds,'Unit','1/s',status)
             IF_ERROR_RETURN(status=1)
             call WriteData(Sds,average_field,status)
             IF_ERROR_RETURN(status=1)
             call Done(Sds,status)
             IF_ERROR_RETURN(status=1)
          end if                  

          call gather(dgrid(region), phot_dat(region)%reff_av, average_field, 0, status)

          if ( isRoot ) then
             call Init(Sds,io,'reff_av',(/imr,jmr,lmr/),'real(4)',status)
             IF_ERROR_RETURN(status=1)
             call WriteAttribute(Sds,'Unit','um',status)
             IF_ERROR_RETURN(status=1)
             call WriteData(Sds,average_field,status)
             IF_ERROR_RETURN(status=1)
             call Done(Sds,status)
             IF_ERROR_RETURN(status=1)
          end if

          call gather(dgrid(region), phot_dat(region)%lwp_av, average_field, 0, status)

          if ( isRoot ) then
             call Init(Sds,io,'lwp_av',(/imr,jmr,lmr/),'real(4)',status)
             IF_ERROR_RETURN(status=1)
             call WriteAttribute(Sds,'Unit','g/m2',status)
             IF_ERROR_RETURN(status=1)
             call WriteData(Sds,average_field,status)
             IF_ERROR_RETURN(status=1)
             call Done(Sds,status)
             IF_ERROR_RETURN(status=1)
          end if
       end if
       
       ! ---- 2D fields
       if ( phot_dat(region)%nalb_av > 0 ) then

          allocate(average_field_2d(imr,jmr))

          call gather(dgrid(region), phot_dat(region)%ags_av, average_field_2d, 0, status)

          if(isRoot) then
             average_field_2d = average_field_2d / phot_dat(region)%nalb_av
             call Init(Sds,io,'ags_av',(/imr,jmr/),'real(4)',status)
             IF_ERROR_RETURN(status=1)
             call WriteAttribute(Sds,'Unit','Fraction',status)
             IF_ERROR_RETURN(status=1)
             call WriteData(Sds,average_field_2d,status)
             IF_ERROR_RETURN(status=1)
             call Done(Sds,status)
             IF_ERROR_RETURN(status=1)
          end if


          call gather(dgrid(region), phot_dat(region)%sza_av, average_field_2d, 0, status)

          if ( isRoot ) then
             average_field_2d = average_field_2d / phot_dat(region)%nalb_av
             call Init(Sds,io,'sza_av',(/imr,jmr/),'real(4)',status)
             IF_ERROR_RETURN(status=1)
             call WriteAttribute(Sds,'Unit','Zenith Angle',status)
             IF_ERROR_RETURN(status=1)
             call WriteData(Sds,average_field_2d,status)
             IF_ERROR_RETURN(status=1)
             call Done(Sds,status)
             IF_ERROR_RETURN(status=1)
          end if

          deallocate(average_field_2d)

       end if ! any 2D field
       

       ! Clouds data (3D)
       if ( phot_dat(region)%ncloud_av > 0 ) then

          call gather(dgrid(region), phot_dat(region)%pmcld_av, average_field, 0, status)

          if ( isRoot ) then
             average_field = average_field/phot_dat(region)%ncloud_av
             call Init(Sds,io,'pmcld_av',(/imr,jmr,lmr/),'real(4)',status)
             IF_ERROR_RETURN(status=1)
             call WriteAttribute(Sds,'Unit',' ',status)
             IF_ERROR_RETURN(status=1)
             call WriteData(Sds,average_field,status)
             IF_ERROR_RETURN(status=1)
             call Done(Sds,status)
             IF_ERROR_RETURN(status=1)
          endif
          
          call gather(dgrid(region), phot_dat(region)%taus_cld_av, average_field, 0, status)

          if ( isRoot ) then
             average_field = average_field/phot_dat(region)%ncloud_av
             call Init(Sds,io,'taus_cld_av',(/imr,jmr,lmr/),'real(4)',status)
             IF_ERROR_RETURN(status=1)
             call WriteAttribute(Sds,'Unit',' ',status)
             IF_ERROR_RETURN(status=1)
             call WriteData(Sds,average_field,status)
             IF_ERROR_RETURN(status=1)
             call Done(Sds,status)
             IF_ERROR_RETURN(status=1)
          end if

          call gather(dgrid(region), phot_dat(region)%taua_cld_av, average_field, 0, status)

          if ( isRoot ) then
             average_field = average_field/phot_dat(region)%ncloud_av
             call Init(Sds,io,'taua_cld_av',(/imr,jmr,lmr/),'real(4)',status)
             IF_ERROR_RETURN(status=1)
             call WriteAttribute(Sds,'Unit',' ',status)
             IF_ERROR_RETURN(status=1)
             call WriteData(Sds,average_field,status)
             IF_ERROR_RETURN(status=1)
             call Done(Sds,status)
             IF_ERROR_RETURN(status=1)
          end if

       end if

       ! Aerosols info (5D)
       if ( phot_dat(region)%naer_av > 0 ) then

          allocate(average_field_5d(imr,jmr,lmr,nbands_trop,ngrid))

          call gather(dgrid(region), phot_dat(region)%pmaer_av, average_field_5d, 0, status)

          if ( isRoot ) then
             average_field_5d = average_field_5d/phot_dat(region)%naer_av
             call Init(Sds,io,'pmaer_av',(/imr,jmr,lmr,nbands_trop,ngrid/),'real(4)',status)
             IF_ERROR_RETURN(status=1)
             call WriteAttribute(Sds,'Unit',' ',status)
             IF_ERROR_RETURN(status=1)
             call WriteData(Sds,average_field_5d,status)
             IF_ERROR_RETURN(status=1)
             call Done(Sds,status)
             IF_ERROR_RETURN(status=1)
          end if

          call gather(dgrid(region), phot_dat(region)%taus_aer_av, average_field_5d, 0, status)

          if ( isRoot ) then
             average_field_5d = average_field_5d/phot_dat(region)%naer_av
             call Init(Sds,io,'taus_aer_av',(/imr,jmr,lmr,nbands_trop,ngrid/),'real(4)',status)
             IF_ERROR_RETURN(status=1)
             call WriteAttribute(Sds,'Unit',' ',status)
             IF_ERROR_RETURN(status=1)
             call WriteData(Sds,average_field_5d,status)
             IF_ERROR_RETURN(status=1)
             call Done(Sds,status)
             IF_ERROR_RETURN(status=1)
          end if

          call gather(dgrid(region), phot_dat(region)%taua_aer_av, average_field_5d, 0, status)

          if ( isRoot ) then
             average_field_5d = average_field_5d/phot_dat(region)%naer_av
             call Init(Sds,io,'taua_aer_av',(/imr,jmr,lmr,nbands_trop,ngrid/),'real(4)',status)
             IF_ERROR_RETURN(status=1)
             call WriteAttribute(Sds,'Unit',' ',status)
             IF_ERROR_RETURN(status=1)
             call WriteData(Sds,average_field_5d,status)
             IF_ERROR_RETURN(status=1)
             call Done(Sds,status)
             IF_ERROR_RETURN(status=1)
          end if
          
          deallocate(average_field_5d)
       end if



       if ( phot_dat(region)%nsad_av >0 ) then
  
          call gather(dgrid(region), phot_dat(region)%sad_cld_av, average_field, 0, status)
          
          if ( isRoot ) then
             average_field = average_field/phot_dat(region)%nsad_av
             call Init(Sds,io,'sad_cld',(/imr,jmr,lmr/),'real(4)',status)
             IF_ERROR_RETURN(status=1)
             call WriteAttribute(Sds,'Unit','cm^2/cm^3 ',status)
             call WriteAttribute(Sds,'long_name','average water cloud surface area density',status)
             IF_ERROR_RETURN(status=1)
             call WriteData(Sds,average_field,status)
             IF_ERROR_RETURN(status=1)
             call Done(Sds,status)
             IF_ERROR_RETURN(status=1)
          endif
          
          call gather(dgrid(region), phot_dat(region)%sad_ice_av, average_field, 0, status)
          
          if ( isRoot ) then
             average_field = average_field/phot_dat(region)%nsad_av
             call Init(Sds,io,'sad_ice',(/imr,jmr,lmr/),'real(4)',status)
             IF_ERROR_RETURN(status=1)
             call WriteAttribute(Sds,'Unit','cm^2/cm^3 ',status)
             call WriteAttribute(Sds,'long_name','average ice cloud surface area density',status)
             IF_ERROR_RETURN(status=1)
             call WriteData(Sds,average_field,status)
             IF_ERROR_RETURN(status=1)
             call Done(Sds,status)
             IF_ERROR_RETURN(status=1)
          endif
          
          call gather(dgrid(region), phot_dat(region)%sad_aer_av, average_field, 0, status)

          if ( isRoot ) then
             average_field = average_field/phot_dat(region)%nsad_av
             call Init(Sds,io,'sad_aer',(/imr,jmr,lmr/),'real(4)',status)
             IF_ERROR_RETURN(status=1)
             call WriteAttribute(Sds,'Unit','cm^2/cm^3 ',status)
             call WriteAttribute(Sds,'long_name','average aerosol surface area density',status)
             IF_ERROR_RETURN(status=1)
             call WriteData(Sds,average_field,status)
             IF_ERROR_RETURN(status=1)
             call Done(Sds,status)
             IF_ERROR_RETURN(status=1)
          endif
       end if
#endif
       
       deallocate(average_field)

    end do REG

#ifdef with_hdf4    
    if ( isRoot ) then
       call Done(io, status)
       IF_ERROR_RETURN(status=1)
    endif
#endif
    
    do region = 1,nregions

       if (associated( phot_dat(region)%o3klim_top)) deallocate(phot_dat(region)%o3klim_top )
       deallocate(phot_dat(region)%albedo     )       
       deallocate(phot_dat(region)%aero_surface_clim)            
       deallocate(phot_dat(region)%cloud_reff )
       deallocate(phot_dat(region)%pmcld      )
       deallocate(phot_dat(region)%taus_cld   )
       deallocate(phot_dat(region)%taua_cld   )
       deallocate(phot_dat(region)%pmaer      )
       deallocate(phot_dat(region)%taus_aer   )
       deallocate(phot_dat(region)%taua_aer   )
       deallocate(phot_dat(region)%lwp_av     )
       deallocate(phot_dat(region)%cfr_av     )
       deallocate(phot_dat(region)%reff_av    )
       deallocate(phot_dat(region)%ags_av     )
       deallocate(phot_dat(region)%sza_av     )
       deallocate(phot_dat(region)%jo3        )
       deallocate(phot_dat(region)%jno2       )
       deallocate(phot_dat(region)%jhno2      )
       deallocate(phot_dat(region)%jch2oa     )
       deallocate(phot_dat(region)%jch2ob     )
       deallocate(phot_dat(region)%jo3_av     )
       deallocate(phot_dat(region)%jno2_av    )
       deallocate(phot_dat(region)%jhno2_av   )
       deallocate(phot_dat(region)%jch2oa_av  )
       deallocate(phot_dat(region)%jch2ob_av  )
       
       deallocate(phot_dat(region)%taua_aer_av)
       deallocate(phot_dat(region)%taus_aer_av)
       deallocate(phot_dat(region)%pmaer_av   )
       deallocate(phot_dat(region)%taua_cld_av)
       deallocate(phot_dat(region)%taus_cld_av)
       deallocate(phot_dat(region)%pmcld_av   )
       deallocate(phot_dat(region)%v2         )
       deallocate(phot_dat(region)%v3         )
       deallocate(phot_dat(region)%v3_surf    )
       deallocate(phot_dat(region)%qy_ch3cocho)
       deallocate(phot_dat(region)%qy_c2o3    )
       deallocate(phot_dat(region)%sad_cld    )
       deallocate(phot_dat(region)%sad_ice    )
       deallocate(phot_dat(region)%sad_aer    )
       deallocate(phot_dat(region)%sad_cld_av )
       deallocate(phot_dat(region)%sad_ice_av )
       deallocate(phot_dat(region)%sad_aer_av )

    end do

#ifdef with_optics
     deallocate( wdep )  !NB optics_done
#endif
    ! ok
    status = 0 

  END SUBROUTINE PHOTOLYSIS_DONE
  !EOC

END MODULE PHOTOLYSIS