!### macro's #####################################################
!
#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"
!
!#################################################################

module photolysis

  use GO, only : gol, goPr, goErr
  
!  use binas,           only : pi
!  use photolysis_data, only : phot_dat
!  use dims

  implicit none
  
  
  ! --- in/out -------------------------------------

  private

  !public :: get_theta, nbands, fparam, photolysis_rates, phot_dat
  public :: init_photolysis, end_photolysis
  !public :: cloud_info, ozone_info

  ! --- const -----------------------
  
  character(len=*), parameter  ::  mname = 'photolysis'



!  ! vars
!
!  integer,parameter      :: nbands=7
!  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   
!
!  real,dimension(nbands),parameter :: crs_o3 = &
!       (/ 3.6448e-19, 1.7115e-18, 2.8218e-19, 1.1091e-19, &
!          2.8145e-20, 1.9000e-23, 4.5500e-21/)
!  ! crs_o2: cross-section o2 205 nm
!  real,parameter                        :: crs_o2= 7.63e-24
!
!  ! reference tables (pressure,mu0,band)
!  real,dimension(0:102,0:19,nbands)     :: fint
!  ! extraterrestial fluxes at selected wavelengths
!  real,dimension(nbands)                :: fe


contains


  subroutine init_photolysis !(iunit,o3du)
    
    use dims           , only : nregions, im, jm
    use dims           , only : newsrun
    use photolysis_data, only : phot_dat
!#ifdef MPI
!    use mpi_const, only: root,localComm, ierr, my_real, myid, lmloc
!#endif
!    use chem_param, only : d23_data

    ! --- in/out --------------------------------------
    
!    ! input
!    integer,intent(in) :: iunit
!    type(d23_data),dimension(nregions),intent(in) :: o3du

    ! --- local ----------------------------------------
    
    integer   ::  region, imr, jmr
    !integer   :: i,m,lmr
    !real      :: lat
    
    ! --- begin ------------------------------------------

!    ! routine to read some reference tables for photolysis calculations
!    ! to be called at the start of any photolysis calculation
!
!    if(newsrun) then
!#ifdef MPI
!       if(myid == root) then
!#endif
!          open(unit=iunit,file='fintwb.data',status='old')
!
!          read(iunit,*)
!          read(iunit,*)
!          read(iunit,*) (fe(i),i=1,7)
!          read(iunit,*) fint
!
!          close(iunit)
!#ifdef MPI
!       end if
!       call mpi_bcast(fint, 103*20*nbands , my_real, root, localComm, ierr)
!       call mpi_bcast(fe,   7 , my_real, root, localComm, ierr)
!#endif
!    end if
!    ! allocate memory photolysis and initialise:

    do region=1,nregions
       imr = im(region) ; jmr = jm(region) 
!#ifdef MPI   
!       lmr = lmloc
!#else
!       lmr = lm(region)
!#endif
!       if(newsrun) allocate(phot_dat(region)%o3klim_top(jmr))
!       phot_dat(region)%o3klim_top(:) = o3du(region)%d23(:,lm(region)) 
       if ( newsrun ) then
!          allocate(phot_dat(region)%vo3_surf  (imr,jmr))
!          allocate(phot_dat(region)%taus      (imr,jmr))
!          allocate(phot_dat(region)%cloudfrac (imr,jmr))
          allocate(phot_dat(region)%albedo    (imr,jmr))
!          allocate(phot_dat(region)%vo3s_av   (imr,jmr))
!          allocate(phot_dat(region)%taus_av   (imr,jmr))
!JEW:      allocate(phot_dat(region)%ags_av    (imr,jmr))
!          allocate(phot_dat(region)%cfr_av    (imr,jmr))
!          phot_dat(region)%albedo  = 0.05
!          phot_dat(region)%cfr_av  = 0.0
!          phot_dat(region)%taus_av = 0.0
!          phot_dat(region)%ags_av  = 0.0
!          phot_dat(region)%ncloud_av = 0
!          phot_dat(region)%vo3s_av = 0.0
!          phot_dat(region)%nvo3s_av= 0
!JEW:      phot_dat(region)%nalb_av= 0
!
!          allocate(phot_dat(region)%vo3 (imr,jmr,lmr))     
!          allocate(phot_dat(region)%tauc(imr,jmr,lmr))
!          allocate(phot_dat(region)%jo3_av(imr,jmr,lm(region)))
!          allocate(phot_dat(region)%jno2_av(imr,jmr,lm(region)))
!          phot_dat(region)%jo3_av  = 0.0
!          phot_dat(region)%jno2_av = 0.0
!          phot_dat(region)%nj_av = 0
       end if
    end do

  end subroutine init_photolysis
  
  
!  real function get_theta(daynr,lon_deg,lat_deg,time)
!    implicit none
!
!    ! input
!    integer,intent(in)       :: daynr
!    real,intent(in)          :: lat_deg,lon_deg
!    real,intent(in)          :: time
!
!    ! local
!    real           :: lon_rad,lat_rad,obliq, deday, delta,houra
!
!    lat_rad=lat_deg*pi/180.
!    lon_rad=lon_deg*pi/180.
!    obliq = 23.45 * pi/180.
!    deday = 4.88 + 2.*pi/365.  * daynr
!    delta = asin(sin(obliq)*sin(deday))
!    houra = lon_rad -pi + time * (2.*pi/24.)
!
!    get_theta= acos(sin(delta)*sin(lat_rad) + &
!         cos(delta)*cos(lat_rad)*cos(houra))*180./pi
!
!  end function get_theta
!
!
!
!  real function sundis(imonth,iday)
!    !-----------------------------------------------------------------------------*
!    !=  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 =*
!    !-----------------------------------------------------------------------------*
!
!    implicit none
!    integer, intent(in)  ::iday,imonth
!
!    ! 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
!
!
!
!  subroutine ozone_info(region)
!    ! ---------------------------------------------------------------------
!    ! subroutine to calculate, from an ozone field,
!    ! the overhead ozone at all grid points.
!    ! upper values are given by the ozone climatology above the highest 
!    ! model mid-level
!    ! MK 01/2003 adapted for new TM5 memory structure.
!    ! ---------------------------------------------------------------------
!    use binas,     only: Avog
!    use dims
!    use photolysis_data, only: phot_dat
!    use global_data,     only: region_dat, mass_dat
!    use chem_param,      only: xmo3,io3
!#ifdef MPI
!    use mpi_const,       only: tracer_active, tracer_loc, my_real
!    use mpi_const,       only: localComm, ierr, proc_tracer
!    use mpi_comm,        only: scatter_after_read_k
!#endif
!    implicit none
!
!    ! input
!    integer, intent(in) :: region
!
!    ! local
!    real,dimension(:,:,:),allocatable     :: o3_overhead ! in #cm-2
!    real,dimension(:,:),pointer           :: o3_overhead_surf !in #cm-2
!    real,dimension(:),pointer             :: ozone_klimtop !in #cm-2
!    real,dimension(:,:,:),pointer         :: o3    !in kg/gridbox
!    real,dimension(:),pointer             :: area  !in m2
!    integer                               :: i,j,l
!    real,parameter                        :: conv = 0.1*Avog/xmo3
!    !  ! conv: from kg/m2 --> #/cm2
!
!    ! start
!
!    o3_overhead_surf => phot_dat(region)%vo3_surf
!    allocate(o3_overhead(im(region),jm(region),lm(region)))
!
!#ifdef MPI
!    if (tracer_active(io3) ) then    
!       o3 => mass_dat(region)%rm_t(1:im(region),1:jm(region),1:lm(region),tracer_loc(io3))
!#else
!       o3 => mass_dat(region)%rm_t(1:im(region),1:jm(region),1:lm(region),io3)
!#endif
!       ozone_klimtop => phot_dat(region)%o3klim_top
!       area => region_dat(region)%dxyp
!
!       do j=1,jm(region)
!          do i=1,im(region)
!             o3_overhead(i,j,lm(region)) = ozone_klimtop(j)  
!          end do
!       end do
!       do l=lm(region)-1,1,-1
!          do j=1,jm(region)
!             do i=1,im(region)
!                o3_overhead(i,j,l) = o3_overhead(i,j,l+1) + &
!                     0.5*conv*(o3(i,j,l)+o3(i,j,l+1))/area(j)
!             end do
!          end do
!       end do
!
!       do j=1,jm(region)
!          do i=1,im(region)
!             o3_overhead_surf(i,j) = o3_overhead(i,j,1) + &
!                  0.5*conv*o3(i,j,1)/area(j)
!          end do
!       end do
!
!#ifdef MPI
!    end if  ! tracer_active(io3)
!#endif
!
!#ifdef MPI
!    call mpi_bcast(o3_overhead_surf, im(region)*jm(region), &
!         my_real, proc_tracer(io3) , localComm, ierr)
!    call scatter_after_read_k(o3_overhead,&
!         im(region),jm(region),lm(region),0,0,0,1,&
!         phot_dat(region)%vo3,proc_tracer(io3))
!#else
!    phot_dat(region)%vo3 = o3_overhead  
!#endif
!
!    phot_dat(region)%vo3s_av  = phot_dat(region)%vo3s_av + o3_overhead_surf
!    phot_dat(region)%nvo3s_av = phot_dat(region)%nvo3s_av + 1
!
!#ifdef MPI
!    if (tracer_active(io3) ) then    
!#endif
!       nullify(o3)
!       nullify(ozone_klimtop)
!       nullify(area)
!#ifdef MPI
!    end if
!#endif
!    deallocate (o3_overhead)
!    nullify(o3_overhead_surf)
!
!  end subroutine ozone_info
!
!
!
!  subroutine cloud_info(region)
!    !
!    ! routine to convert lwc, iwc, and cf to tauc and taus, and cfr
!    ! to be called each time fresh cloud info becomes available
!    ! ADAPTED for the new data structure, MK, 01/2003
!    !
!    use binas,  only : xmair, Avog
!    use photolysis_data, only : phot_dat
!    use global_data,     only : cloud_dat, mass_dat
!#ifdef MPI
!    use mpi_const,   only : lmar, myid, lmloc
!#endif
!    implicit none
!
!    ! input
!    integer, intent(in) :: region
!
!    ! local
!    real,dimension(:,:,:),pointer     :: cloudliqw,cloudicew,cloudcov
!    real,dimension(:,:,:),pointer     :: pressure
!    real,dimension(:,:,:),allocatable :: tau_c 
!    ! tau_c: optical thickness clouds above level
!    real,dimension(:,:),  pointer     :: tau_surf 
!    ! tau_surf: optical thickness clouds at surface
!    real,dimension(:,:),  pointer     :: cloud_frac  !cloud fraction at the surface
!    real,parameter                    :: reff = 8.0   !fd effective cloud radius
!    real,parameter                    :: factor =  7.24e16*1e6*xmair*29.2605/Avog
!    real,dimension(:,:,:),allocatable :: wc      !ice + lwc
!    real                              :: airp,rhodz,maxtau
!    integer                           :: i,j,l,lsave,offsetl
!
!    ! start
!
!    ! calculates cloud optical thickness and cloud fraction
!    ! each time a new cloud field is read in (6 hours)
!
!    cloudliqw => cloud_dat(region)%lwc
!    cloudicew => cloud_dat(region)%iwc
!    cloudcov => cloud_dat(region)%cc
!    pressure => mass_dat(region)%phlb_t
!    tau_surf => phot_dat(region)%taus
!    cloud_frac => phot_dat(region)%cloudfrac
!
!    allocate(tau_c(im(region),jm(region),lm(region)))
!    allocate(wc(im(region),jm(region),lm(region)))
!
!    wc = cloudliqw + cloudicew    ! no differentiation between water and ice....
!
!    do l=1,lm(region)
!       do j=1,jm(region)
!          do i=1,im(region)
!             ! threshold to account for the cloud
!             if ( cloudcov(i,j,l) > 0.02) then 
!                airp = (pressure(i,j,l)+pressure(i,j,l+1))/2.
!                !   aird = 7.24e16*airp/t(i,j,k)           !#/cm3
!                !   rho  = aird*1e6*xmair/avo              !g/m3
!                !   dz   = 29.2605*t(i,j,k)* alog(phlb(i,j,l)/(1.0+phlb(i,j,l+1)))
!                rhodz = factor*airp*alog(pressure(i,j,l)/(1.0+pressure(i,j,l+1)))
!                !CMK for security 
!                tau_c(i,j,l)  = max(0.0,wc(i,j,l)*rhodz*1.5/reff)
!             else
!                tau_c(i,j,l)  = 0.0
!             end if
!          end do
!       end do
!    end do
!    !
!    ! use optical depth to find maximum cloudcover
!    !do j=1,jm(region)
!    !do i=1,im(region) 
!    !  cfr(i,j) = maxval(cf(i,j,:))  !maximum cloud cover in the column
!    !end do
!    !end do
!    !mktest use cloudcover from box with maximum tau
!    !
!    do j=1,jm(region)
!       do i=1,im(region) 
!          maxtau = 0.0
!          lsave = -1
!          do l=1,lm(region) 
!             if(tau_c(i,j,l) > maxtau) then 
!                maxtau = tau_c(i,j,l)
!                lsave = l
!             end if
!          end do
!          if (lsave.ne.-1) then
!             cloud_frac(i,j) = cloudcov(i,j,lsave)
!          else
!             cloud_frac(i,j) = 0.0
!          end if
!       end do
!    end do
!
!
!    cloud_frac=max(1e-7,cloud_frac)! avoid division by zero
!    do l=1,lm(region)
!       do j=1,jm(region)
!          do i=1,im(region)
!             !scale to cloudy fraction (temp in wc)
!             wc(i,j,l) = tau_c(i,j,l)/cloud_frac(i,j)
!          end do
!       end do
!    end do
!
!    ! for radiative transfer: calculate clouds above the level!!
!
!    do j=1,jm(region)
!       do i=1,im(region)
!          tau_c(i,j,lm(region)) = 0.5*wc(i,j,lm(region))
!       end do
!    end do
!
!    !cmk check.....should air mass be used?
!    do l=lm(region)-1,1,-1
!       do j=1,jm(region)
!          do i=1,im(region)
!             tau_c(i,j,l) = tau_c(i,j,l+1) + 0.5*(wc(i,j,l)+wc(i,j,l+1))
!          end do
!       end do
!    end do
!
!
!    ! cloud optical thickness above surface
!
!    do j=1,jm(region)
!       do i=1,im(region)
!          tau_surf(i,j) = tau_c(i,j,1) + 0.5*wc(i,j,1)
!       end do
!    end do
!
!    tau_c = max(tau_c,0.0)
!    tau_surf = max(tau_surf,0.0)
!
!    !...copy tau_c to tau_c which is parallel over the levels (for photolysis)
!    offsetl=0
!#ifdef MPI
!    if(myid>0)offsetl=sum(lmar(0:myid-1) )  ! offset for global value of l
!#endif
!#ifdef MPI
!    phot_dat(region)%tauc = tau_c(:,:,offsetl+1: offsetl+lmloc)
!#else
!    phot_dat(region)%tauc = tau_c(:,:,1:lm(region))
!#endif
!
!
!
!    phot_dat(region)%taus_av = phot_dat(region)%taus_av + tau_surf
!    phot_dat(region)%cfr_av  = phot_dat(region)%cfr_av  + cloud_frac
!    phot_dat(region)%ncloud_av  = phot_dat(region)%ncloud_av  + 1
!
!
!    deallocate(wc)
!    deallocate(tau_c)
!
!    nullify(cloudliqw)
!    nullify(cloudicew)
!    nullify(cloudcov)
!    nullify(pressure)
!    nullify(tau_surf)
!    nullify(cloud_frac)
!
!  end subroutine cloud_info
!
!
!
!  subroutine fparam(region,level,th,mu,fact)
!    !
!    use dims
!    use photolysis_data, only: phot_dat
!    use global_data, only: region_dat, mass_dat
!#ifdef MPI
!    use advect_tools, only : calc_phlb_k
!#endif
!
!    implicit none
!
!    ! input/output
!    integer, intent(in)                            :: region
!    integer, intent(in)                            :: level   ! the level 
!    real,dimension(im(region),jm(region))          :: th,mu  ! sza, cos(sza)
!    real,dimension(im(region),jm(region),nbands),intent(out)       :: fact
!
!    ! local
!    real,dimension(:,:,:),pointer  :: phlb
!    real,dimension(:,:),pointer    :: ps      ! surface pressure
!    real,dimension(:,:,:),pointer  :: vo3     ! overhead ozone at surface
!    real,dimension(:,:),pointer    :: vo3s    ! overhead ozone at the surface
!    real,dimension(:,:,:),pointer  :: tauc    ! optical depth
!    real,dimension(:,:),pointer    :: taus    ! optical depth at surface
!    real,dimension(:,:),pointer    :: cfr     ! cloud fraction in column
!    real,dimension(:,:),pointer    :: ags     ! surface albedo of column
!    integer,dimension(:,:),pointer :: zoomed
!
!    real :: tu,tuc,tl,tlc !optical thicknesses lower/upper
!    real :: etu,etuc,etl,etlc          !exp(tau/mu)
!    real :: eto3 !for ozon or tau
!    real :: ala0,alad
!    real :: ala0c,aladc,al0,ald,al0c,aldc
!    real :: eto3u
!    real :: au,au0,au0c,auc,ag
!    real :: f0,f0c,fs,fsc
!    real :: fu,fuc,fd,fdc
!    real :: cref
!    real,parameter      :: one_g = 1.-0.87
!    real,parameter      :: one_g2 = 1.-0.87*0.87
!    integer  :: rd=2
!
!    integer,dimension(4),parameter  :: band= (/ 3,4,5,7 /)
!    real,dimension(3:7),parameter   :: taur= &
!         (/ 1.177e-5,1.067e-5,0.919e-5,0.498e-5,0.0782e-5 /)
!
!    integer  :: i,j,ilon,ilat,indi,ip,ith ,ipres,k,jb
!    real     :: o3c_surf,mu0,pp,dp,dth,tth,taui,vo3i,vo2,p
!
!    ! start
!
!#ifdef MPI
!    phlb => mass_dat(region)%phlb_k
!#else
!    phlb => mass_dat(region)%phlb_t
!#endif
!    ps   => mass_dat(region)%p
!    vo3  => phot_dat(region)%vo3
!    vo3s => phot_dat(region)%vo3_surf
!    tauc => phot_dat(region)%tauc
!    taus => phot_dat(region)%taus
!    cfr  => phot_dat(region)%cloudfrac
!    ags  => phot_dat(region)%albedo
!    zoomed => region_dat(region)%zoomed
!
!#ifdef MPI
!    call calc_phlb_k(region)
!#endif
!
!    jb = 1   ! first Wavelength band
!    do j=jsr(region),jer(region)
!       do i=isr(region),ier(region)
!          if(zoomed(i,j)/=region) cycle
!          p = (phlb(i,j,level+1) + phlb(i,j,level))*0.5   !pressure
!          if(th(i,j).lt.90.0.and.p.lt.2e4) then
!
!             vo2  = p*sp
!             !ozone and oxygen
!             eto3 = exp(-(crs_o2*vo2 + crs_o3(jb)*vo3(i,j,level))/mu(i,j))
!
!             pp = max(102-1e-3*p,0.0)
!             ip = int(pp)
!             dp = pp-ip
!             tth = min(th(i,j),89.9999)*0.2
!             ith = int(tth)
!             dth = tth-ith
!             cref =( (1-dp)*((1-dth)*fint(ip  ,ith  ,jb)      &
!                  + dth *fint(ip  ,ith+1,jb))     &
!                  +  dp *((1-dth)*fint(ip+1,ith  ,jb)      &
!                  + dth *fint(ip+1,ith+1,jb)) )
!
!             ! actinic flux only depends on ozone absorption and ozygen absorption
!
!             fact(i,j,jb) =  cref*eto3
!
!          else
!
!             fact(i,j,jb) =  0.0 
!
!          end if
!
!       end do
!    end do
!
!    jb=2
!    do j=jsr(region),jer(region)
!       do i=isr(region),ier(region)
!          if(zoomed(i,j)/=region) cycle
!
!          p = (phlb(i,j,level+1) + phlb(i,j,level))*0.5   !pressure
!          if(th(i,j).lt.90.0.and.p.lt.2e4) then
!
!
!             eto3 = exp(-crs_o3(jb)*vo3(i,j,level)/mu(i,j))
!
!             pp = max(102-1e-3*p,0.0)
!             ip = int(pp)
!             dp = pp-ip
!             tth = min(th(i,j),89.9999)*0.2
!             ith = int(tth)
!             dth = tth-ith
!             cref =( (1-dp)*((1-dth)*fint(ip  ,ith  ,jb)      &
!                  + dth *fint(ip  ,ith+1,jb))     &
!                  +  dp *((1-dth)*fint(ip+1,ith  ,jb)      &
!                  + dth *fint(ip+1,ith+1,jb)) )
!
!             ! actinic flux only depends on ozone absorption
!
!             fact(i,j,jb) =  cref*eto3
!
!          else
!
!             fact(i,j,jb) = 0.0
!
!          end if
!       end do
!    end do !_____jb = 2
!
!    do k=1,4   !over wavelength intervals
!       jb = band(k)
!       do j=jsr(region),jer(region)
!          do i=isr(region),ier(region)
!             if(zoomed(i,j)/=region) cycle
!
!             if(th(i,j).lt.90.0) then 
!
!
!                taui = tauc(i,j,level) 
!                vo3i = vo3(i,j,level)
!                mu0  = mu(i,j)
!                ag   = ags(i,j)
!                p = (phlb(i,j,level+1) + phlb(i,j,level))*0.5   !pressure
!                tu   = taur(jb)*p     !optical thickness upper atmosphere clear
!                tuc  = tu + one_g*taui
!                tl   = taur(jb)*(ps(i,j) - p)             !lower atmosphere
!                tlc  = tl + one_g*(taus(i,j)-taui)           !clouds included
!                eto3      = exp(-crs_o3(jb)*vo3i/mu0)
!                eto3u     = exp(-crs_o3(jb)*(vo3s(i,j)-vo3i)*2.5)
!                etu       = exp(-tu/mu0)
!                etuc      = etu*exp(-taui*one_g2/mu0)   !includes cloud effects
!                etl       = exp(-tl/mu0)
!                etlc      = etl*exp(-(taus(i,j)-taui)*one_g2/mu0)
!                ala0      = 1.-(2.+3.*mu0+(2.-3.*mu0)*etl)/(4.+3.*tl)
!                alad      = 1.-4./(4.+3.*tl)
!                ala0c     = 1.-(2.+3.*mu0+(2.-3.*mu0)*etlc)/(4.+3.*tlc)
!                aladc     = 1.-4./(4.+3.*tlc)
!                al0       = ala0  + (1.-ala0)  *ag*  (1.-alad)     /(1.-ag*alad)
!                ald       = alad  + (1.-alad)  *ag*  (1.-alad)     /(1.-ag*alad)
!                al0c      = ala0c + (1.-ala0c) *ag*  (1.-aladc)    /(1.-ag*aladc)
!                aldc      = aladc + (1.-aladc) *ag*  (1.-aladc)    /(1.-ag*aladc)
!                au        = 1.-4./(4.+3*tu)
!                au0       = 1.-(2.+3.*mu0+(2.-3.*mu0)*etu) /(4.+3.*tu)
!                au0c      = 1.-(2.+3.*mu0+(2.-3.*mu0)*etuc)/(4.+3.*tuc)
!                auc       = 1.-4./(4.+3*tuc)
!                !
!                ! direct radiation
!                !
!                f0        = eto3*etu
!                f0c       = eto3*etuc
!                !
!                ! downward diffuse, only from events above p(i)
!                !
!
!                fs        = rd *mu0*(1.-au0-etu)  *eto3
!                fsc       = rd *mu0*(1.-au0c-etuc)*eto3
!                !
!                !correct lower albedo for ozone absorbtion
!                !
!
!                al0       = al0     *eto3u      
!                ald       = ald     *eto3u
!                al0c      = al0c    *eto3u
!                aldc      = aldc    *eto3u
!                !
!                ! upward diffuse radiation 
!                !
!                fu       = (rd *mu0*f0  *al0   + ald   *fs)     /(1.-ald   *au)
!                fuc      = (rd *mu0*f0c *al0c  + aldc  *fsc)    /(1.-aldc  *auc)
!                !
!                ! downward direct + diffuse radiation
!                !
!                fd       = f0     + fs     + au  *fu
!                fdc      = f0c    + fsc    + auc *fuc
!                !
!                !  interpolate actinic flux table
!                !
!                pp = max(102-1e-3*p,0.0)
!                ip = int(pp)
!                dp = pp-ip
!                tth = min(th(i,j),89.9999)*0.2
!                ith = int(tth)
!                dth = tth-ith
!                cref =( (1-dp)*((1-dth)*fint(ip  ,ith  ,jb)      &
!                     + dth *fint(ip  ,ith+1,jb))     &
!                     +  dp *((1-dth)*fint(ip+1,ith  ,jb)      &
!                     + dth *fint(ip+1,ith+1,jb)) )
!                !
!                ! calculate actinic flux, weighted with the cloud fraction
!                ! 
!                fact(i,j,jb) = (cfr(i,j)*(fdc + fuc) + &
!                     (1.-cfr(i,j))*(fu+fd)) *cref
!             else
!
!                fact(i,j,jb) = 0.0
!
!             end if
!          end do
!       end do
!    end do
!
!    jb = 6    !no ozone absorption
!    do j=jsr(region),jer(region)
!       do i=isr(region),ier(region)
!          if(zoomed(i,j)/=region) cycle
!
!          p = (phlb(i,j,level+1) + phlb(i,j,level))*0.5   !pressure
!          if(th(i,j).lt.90.0) then 
!
!             taui = tauc(i,j,level)
!             mu0  = mu(i,j)
!             ag   = ags(i,j)
!             tu   = taur(jb)*p       !optical thickness upper atmosphere clear
!             tuc  = tu + one_g*taui
!             tl   = taur(jb)*(ps(i,j) - p)             !lower atmosphere
!             tlc  = tl + one_g*(taus(i,j)-taui)           !clouds included
!             etu       = exp(-tu/mu0)
!             etuc      = etu*exp(-taui*one_g2/mu0)     !includes cloud effects
!             etl       = exp(-tl/mu0)
!             etlc      = etl*exp(-(taus(i,j)-taui)*one_g2/mu0)  
!             ala0      = 1.-(2.+3.*mu0+(2.-3.*mu0)*etl)/(4.+3.*tl)
!             alad      = 1.-4./(4.+3.*tl)
!             ala0c     = 1.-(2.+3.*mu0+(2.-3.*mu0)*etlc)/(4.+3.*tlc)
!             aladc     = 1.-4./(4.+3.*tlc)
!             al0       = ala0  + (1.-ala0)  *ag*  (1.-alad)   /(1.-ag*alad)
!             ald       = alad  + (1.-alad)  *ag*  (1.-alad)   /(1.-ag*alad)
!             al0c      = ala0c + (1.-ala0c) *ag*  (1.-aladc)  /(1.-ag*aladc)
!             aldc      = aladc + (1.-aladc) *ag*  (1.-aladc)  /(1.-ag*aladc)
!             au        = 1.-4./(4.+3*tu)
!             au0       = 1.-(2.+3.*mu0+(2.-3.*mu0)*etu) /(4.+3.*tu)
!             au0c      = 1.-(2.+3.*mu0+(2.-3.*mu0)*etuc)/(4.+3.*tuc)
!             auc       = 1.-4./(4.+3*tuc)
!             !
!             ! direct radiation
!             !
!             f0        = etu
!             f0c       = etuc
!             !
!             ! downward diffuse, only from events above p(i)
!             !
!             fs        = rd *mu0*(1.-au0-etu)  
!             fsc       = rd *mu0*(1.-au0c-etuc)
!             !
!             ! upward diffuse radiation 
!             !
!             fu       = (rd *mu0*f0   *al0   + ald  *fs)  /(1.-ald  *au)
!             fuc      = (rd *mu0*f0c  *al0c  + aldc *fsc) /(1.-aldc *auc)
!             !
!             ! downward direct + diffuse radiation
!             !
!             fd       = f0     + fs     + au  *fu
!             fdc      = f0c    + fsc    + auc *fuc
!             !
!             !  interpolate actinic flux table
!             !
!             pp = max(102-1e-3*p,0.0)
!             ip = int(pp)
!             dp = pp-ip
!             tth = min(th(i,j),89.9999)*0.2
!             ith = int(tth)
!             dth = tth-ith
!             cref =( (1-dp)*((1-dth)*fint(ip  ,ith  ,jb)      &
!                  + dth *fint(ip  ,ith+1,jb))     &
!                  +  dp *((1-dth)*fint(ip+1,ith  ,jb)      &
!                  + dth *fint(ip+1,ith+1,jb)) )
!             !
!             ! calculate actinic flux, weighted with the cloud fraction
!             ! 
!             fact(i,j,jb) = (cfr(i,j)*(fdc + fuc) + (1.-cfr(i,j))*(fu+fd)) *cref
!             !cmk   fact(i,j) = (fu+fd) *cref
!
!          else
!
!             fact(i,j,jb) = 0.0
!
!          end if
!
!       end do
!    end do   !___jb=6
!
!
!    nullify(phlb)
!    nullify(ps)
!    nullify(vo3)
!    nullify(vo3s)
!    nullify(tauc)
!    nullify(taus)
!    nullify(cfr)
!    nullify(ags)
!    nullify(zoomed)
!
!  end subroutine fparam
!
!
!
!  subroutine photolysis_rates(region,level,th,mu,fact,rjx)
!    !
!    use dims
!#ifdef MPI
!    use mpi_const,   only : lmar, myid, lmloc
!#endif
!    use global_data, only: mass_dat, region_dat, meteo_dat
!    use photolysis_data, only: nj   ,jo3d,jno2,jh2o2,jhno3,jhno4,jn2o5,jach2o, &
!         jbch2o,jmepe,jano3,jbno3,jpan,jorgn,j45,j74,jrooh, phot_dat
!
!    implicit none
!
!    ! input/output
!    !region & model level
!    integer,intent(in)                                      :: region
!    integer,intent(in)                                      :: level
!    ! fact: actinic flux calculated by fparam...
!    real,dimension(im(region),jm(region),nbands),intent(in) :: fact
!    ! th: sza (deg) mu: cos(sza)
!    real,dimension(im(region),jm(region)),intent(in)        :: th
!    real,dimension(im(region),jm(region)),intent(in)        :: mu
!    ! rjx: photolysis rates
!    real,dimension(im(region),jm(region),nj   ),intent(out) :: rjx
!
!    ! local
!    real,dimension(:,:,:), pointer                 :: phlb
!    integer,dimension(:,:)  , pointer              :: zoomed
!    ! t: temprature (K)    NOTE: this if in parallel full 3D array
!    real,dimension(:,:,:), pointer                 :: t
!    ! vo3: overhead ozone (#/cm2)
!    real,dimension(:,:,:), pointer                 :: vo3
!
!    ! ratio  incoming f0 (band)/ incoming f0(li)  
!    ! depends on extraterrestial flux chosen
!
!    ! test works with fint2.data real,dimension(1:7),parameter ::f0 = &
!    !   (/ 70.61, 4.850, 19.81, 9.303, 24.86, 18.28, 48.04/)
!    ! f0: bruhl wavelength works with fintwb.data
!    real,dimension(nbands),parameter  :: f0 = &
!         (/ 73.44, 4.977, 16.13, 8.017, 27.9,  19.00, 48.12/)
!
!    real, parameter    :: densref = 1.3e19
!
!    real    :: dens,vo2
!    real    :: tau_2,tau_3,tau_4,tau_5,tau_7             !optical depth
!    real    :: tau_1,sig_1,fint_1
!    real    :: sig_2,sig_3,sig_4,sig_5,sig_7,sig_6       !effective cross-sections
!    real    :: fint_2,fint_3,fint_4   !integrated actinic fluxes 
!    real    :: fint_5,fint_6,fint_7   !integrated actinic fluxes 
!    real    :: vo3s,v3_du2,v3_du1,v3_du,v3s1,v3s2        !ozone column variables
!    real    :: vo2s,v2s_m,v2s2                       !oxygen column variables
!    real    :: tscale,tscale2,pscale      !temperature and pressure scales
!    real    :: p2s,p3s,p4,p5,p6,p7        !tempvar for crosssection calculations
!    real    :: cc0,cc1,cc2
!    real    :: pres,esrm2
!    integer :: i1,i2,i3,ind                            !indices lookup tables
!    integer :: i,j                                     !count
!    integer :: ilon,ilat,ipres,offsetl
!
!    include 'photolysis.i'
!
!    ! start
!
!#ifdef MPI
!    phlb => mass_dat(region)%phlb_k
!#else
!    phlb => mass_dat(region)%phlb_t
!#endif
!    zoomed => region_dat(region)%zoomed
!
!    offsetl=0
!#ifdef MPI
!    if(myid>0)offsetl=sum(lmar(0:myid-1) )  ! offset for global value of l
!#endif
!
!#ifdef MPI
!    t => meteo_dat(region)%t(1:im(region),1:jm(region),offsetl+1:offsetl+lmloc)
!#else
!    t => meteo_dat(region)%t(1:im(region),1:jm(region),1:lm(region))
!#endif
!    !print *, 'CHECK in photolysis:', myid, lbound(t),ubound(t)
!    vo3 => phot_dat(region)%vo3
!
!    rjx  = 0.0
!
!    do j=jsr(region),jer(region)
!       do i=isr(region),ier(region)
!          if(zoomed(i,j)/=region) cycle
!
!          if(th(i,j).lt.90.) then
!
!             pres  = 0.5*(phlb(i,j,level+1)+phlb(i,j,level))
!             vo2= pres*sp  !make oxygen conc. in #/cm3
!             vo3s  = vo3(i,j,level)/mu(i,j)
!
!             vo2s  = vo2/mu(i,j)
!             v3_du = vo3s*todu
!             v2s_m = vo2s*to_m
!
!             !   interval: 1
!
!             if (v2s_m.ge.500.) then
!                v2s_m = 500.
!                v2s2  = 1.3602e+24
!             else
!                v2s2  = vo2s
!             end if
!
!             !   interval: 1+2
!             if (v3_du.ge.300.) then
!                v3_du1 = 300.
!                v3s1   = 8.161e+18
!             else
!                v3_du1 = v3_du
!                v3s1   = vo3s
!             end if
!
!             !   interval: 3 - 7
!             if (v3_du.ge.3000.) then
!                v3_du2   = 3000.
!                v3s2     = 8.161e+19
!             else
!                v3_du2   = v3_du
!                v3s2     = vo3s
!             end if
!
!             !   get indices for lookup-tables
!             ind  = min(int((v3_du-0.5)/2.5) +1,115)
!             i2   = ifil(ind)
!             i1   = i2
!             ind  = min(int((v3_du-5)/25) +1,115)
!             i3   = ifil(ind)
!
!             !---calculate optical depth due to ozone absorbtion
!
!             tau_1  = p1(taub1(i1,1),taua1(i1,1),v3_du1) + &
!                  p1(taub1(i1,2),taua1(i1,2),v3_du1) * v2s_m + &
!                  p1(taub1(i1,3),taua1(i1,3),v3_du1) * v2s_m**2 
!             tau_2  = p1(taub2(i2),taua2(i2),v3_du1) 
!             tau_3  = p1(taub3(i3),taua3(i3),v3_du2) 
!             tau_4  = p3(c4_tau(1),c4_tau(2),c4_tau(3),c4_tau(4),v3_du2)
!             tau_5  = p2(c5_tau(1),c5_tau(2),c5_tau(3),v3_du2)
!             tau_7  = p1(c7_tau(1),c7_tau(2),v3_du2)
!
!             !---get effective actinic fluxes in units quanta/s/cm2
!
!             fint_1 = exp(-tau_1 + crs_o3(1)*v3s1 + &
!                  crs_o2   *v2s2)*f0(1)*fact(i,j,1)*fe(1) 
!             fint_2 = exp(-tau_2 + crs_o3(2)*v3s1)*f0(2)*fact(i,j,2)*fe(2) 
!             fint_3 = exp(-tau_3 + crs_o3(3)*v3s2)*f0(3)*fact(i,j,3)*fe(3)
!             fint_4 = exp(-tau_4 + crs_o3(4)*v3s2)*f0(4)*fact(i,j,4)*fe(4)
!             fint_5 = exp(-tau_5 + crs_o3(5)*v3s2)*f0(5)*fact(i,j,5)*fe(5)
!             fint_6 =                              f0(6)*fact(i,j,6)*fe(6)
!             fint_7 = exp(-tau_7 + crs_o3(7)*v3s2)*f0(7)*fact(i,j,7)*fe(7)
!
!             tscale = (t(i,j,level)-250.)/250.
!             dens   = 7.24e16*pres/t(i,j,level)
!             pscale = (dens-densref)/densref
!
!             !--------------------
!             ! o3 -> o(1d) + o2
!             !--------------------
!             p2s     = p1(b2_o1d(i2),a2_o1d(i2),v3_du1)
!             p3s     = p1(b3_o1d(i3),a3_o1d(i3),v3_du2)
!             p4      = p3(c4_o1d(1),c4_o1d(2),c4_o1d(3),c4_o1d(4),v3_du2)
!             p5      = p3(c5_o1d(1),c5_o1d(2),c5_o1d(3),c5_o1d(4),v3_du2)
!
!             sig_2 = p3(tj_o1d(1,2),tj_o1d(2,2),tj_o1d(3,2),tj_o1d(4,2),tscale) &
!                  * p2s
!             sig_3 = p3(tj_o1d(1,3),tj_o1d(2,3),tj_o1d(3,3),tj_o1d(4,3),tscale) &
!                  * p3s
!             sig_4 = p3(tj_o1d(1,4),tj_o1d(2,4),tj_o1d(3,4),tj_o1d(4,4),tscale) &
!                  * p4
!             sig_5 = p3(tj_o1d(1,5),tj_o1d(2,5),tj_o1d(3,5),tj_o1d(4,5),tscale) &
!                  * p5
!
!             rjx(i,j,jo3d) =  sig_2 * fint_2 &
!                  +  sig_3 * fint_3 &
!                  +  sig_4 * fint_4 &
!                  +  sig_5 * fint_5 
!             !--------------------
!             ! no2
!             !--------------------
!
!             sig_2 =p1(tj_no2(1,2),tj_no2(2,2),tscale) * &
!                  p1(b2_no2(i2),a2_no2(i2),v3_du1)
!             sig_3 =p1(tj_no2(1,3),tj_no2(2,3),tscale) * &
!                  p1(b3_no2(i3),a3_no2(i3),v3_du2)
!             sig_4 =p1(tj_no2(1,4),tj_no2(2,4),tscale) * &
!                  p2(c4_no2(1) ,c4_no2(2),c4_no2(3),v3_du2)
!             sig_5 =p1(tj_no2(1,5),tj_no2(2,5),tscale) * &
!                  p2(c5_no2(1) ,c5_no2(2),c5_no2(3),v3_du2)
!             sig_6 =p1(tj_no2(1,6),tj_no2(2,6),tscale) * &
!                  p1(c6_no2(1) ,c6_no2(2),v3_du2)
!
!             rjx(i,j,jno2) =  sig_2 * fint_2 &
!                  +  sig_3 * fint_3 &
!                  +  sig_4 * fint_4 &
!                  +  sig_5 * fint_5 &
!                  +  sig_6 * fint_6 
!
!             !--------------------
!             ! h2o2
!             !--------------------
!
!             sig_1 = p1(b1_h2o2(i1,1),a1_h2o2(i1,1),v3_du1) + &
!                  p1(b1_h2o2(i1,2),a1_h2o2(i1,2),v3_du1)* v2s_m
!             sig_2 = p2(tj_h2o2(1,2),tj_h2o2(2,2),tj_h2o2(3,2),tscale) * &
!                  p1(b2_h2o2(i2),a2_h2o2(i2),v3_du1)
!             sig_3 = p2(tj_h2o2(1,3),tj_h2o2(2,3),tj_h2o2(3,3),tscale) * &
!                  p1(b3_h2o2(i3),a3_h2o2(i3),v3_du2)
!             sig_4 = p2(tj_h2o2(1,4),tj_h2o2(2,4),tj_h2o2(3,4),tscale) * &
!                  p2(c4_h2o2(1),c4_h2o2(2),c4_h2o2(3),v3_du2)
!             sig_5 = p2(tj_h2o2(1,5),tj_h2o2(2,5),tj_h2o2(3,5),tscale) * &
!                  p2(c5_h2o2(1),c5_h2o2(2),c5_h2o2(3),v3_du2)
!             sig_6 = p2(tj_h2o2(1,6),tj_h2o2(2,6),tj_h2o2(3,6),tscale) * &
!                  p1(c6_h2o2(1),c6_h2o2(2),v3_du2)
!
!             rjx(i,j,jh2o2) =  sig_1 * fint_1 &
!                  +  sig_2 * fint_2 &
!                  +  sig_3 * fint_3 &
!                  +  sig_4 * fint_4 &
!                  +  sig_5 * fint_5 &
!                  +  sig_6 * fint_6 
!
!             !--------------------
!             ! hno3
!             !--------------------
!
!             sig_1 =(p1(b1_hno3(i1,1),a1_hno3(i1,1), v3_du1) + &
!                  p1(b1_hno3(i1,2),a1_hno3(i1,2), v3_du1)*v2s_m)* &
!                  p3(tj_hno3(1,1),tj_hno3(2,1),tj_hno3(3,1),tj_hno3(4,1),tscale)
!             sig_2 = &
!                  p3(tj_hno3(1,2),tj_hno3(2,2),tj_hno3(3,2),tj_hno3(4,2),tscale) &
!                  * p1(b2_hno3(i2),a2_hno3(i2),v3_du1)
!             sig_3 = &
!                  p3(tj_hno3(1,3),tj_hno3(2,3),tj_hno3(3,3),tj_hno3(4,3),tscale) &
!                  * p1(b3_hno3(i3),a3_hno3(i3),v3_du2)
!             sig_4 = &
!                  p3(tj_hno3(1,4),tj_hno3(2,4),tj_hno3(3,4),tj_hno3(4,4),tscale) &
!                  * p2(c4_hno3(1),c4_hno3(2),c4_hno3(3),v3_du2)
!             sig_5 = &
!                  p3(tj_hno3(1,5),tj_hno3(2,5),tj_hno3(3,5),tj_hno3(4,5),tscale) &
!                  * p2(c5_hno3(1),c5_hno3(2),c5_hno3(3),v3_du2)
!             sig_6 = &
!                  p3(tj_hno3(1,6),tj_hno3(2,6),tj_hno3(3,6),tj_hno3(4,6),tscale) &
!                  * p1(c6_hno3(1),c6_hno3(2),v3_du2)
!
!             rjx(i,j,jhno3) =  sig_1 * fint_1 &
!                  +  sig_2 * fint_2 &
!                  +  sig_3 * fint_3 &
!                  +  sig_4 * fint_4 &
!                  +  sig_5 * fint_5 &
!                  +  sig_6 * fint_6 
!
!             !--------------------
!             ! hno4
!             !--------------------
!
!             sig_1 = p1(b1_hno4(i1,1),a1_hno4(i1,1), v3_du1) + &
!                  p1(b1_hno4(i1,2),a1_hno4(i1,2), v3_du1) *v2s_m
!             sig_2 = p1(b2_hno4(i2), a2_hno4(i2),v3_du1)
!             sig_3 = p1(b3_hno4(i3), a3_hno4(i3),v3_du2)
!             sig_4 = p2(c4_hno4(1),c4_hno4(2),c4_hno4(3),v3_du2)
!             sig_5 = p2(c5_hno4(1),c5_hno4(2),c5_hno4(3),v3_du2)
!
!             rjx(i,j,jhno4) =  sig_1 * fint_1 &
!                  +  sig_2 * fint_2 &
!                  +  sig_3 * fint_3 &
!                  +  sig_4 * fint_4 &
!                  +  sig_5 * fint_5
!
!             !--------------------
!             ! n2o5
!             !--------------------
!
!             tscale2 = (max(min(t(i,j,level),300.0),225.0)-250.0)/250.0
!             sig_1 = p1(b1_n2o5(i1,1),a1_n2o5(i1,1), v3_du1) + &
!                  p1(b1_n2o5(i1,2),a1_n2o5(i1,2), v3_du1) *v2s_m
!             sig_2 = p3(tj_n2o5(1,2),tj_n2o5(2,2),tj_n2o5(3,2), &
!                  tj_n2o5(4,2),tscale2) * p1(b2_n2o5(i2),a2_n2o5(i2),v3_du1)
!             sig_3 = p3(tj_n2o5(1,3),tj_n2o5(2,3),tj_n2o5(3,3), &
!                  tj_n2o5(4,3),tscale2) * p1(b3_n2o5(i3),a3_n2o5(i3),v3_du2)
!             sig_4 = p3(tj_n2o5(1,4),tj_n2o5(2,4),tj_n2o5(3,4), &
!                  tj_n2o5(4,4),tscale2) * &
!                  p2(c4_n2o5(1),c4_n2o5(2),c4_n2o5(3),v3_du2)
!             sig_5 = p3(tj_n2o5(1,5),tj_n2o5(2,5),tj_n2o5(3,5), &
!                  tj_n2o5(4,5),tscale2) * &
!                  p2(c5_n2o5(1),c5_n2o5(2),c5_n2o5(3),v3_du2)
!             sig_6 = p3(tj_n2o5(1,6),tj_n2o5(2,6),tj_n2o5(3,6), &
!                  tj_n2o5(4,6),tscale2) * p1(c6_n2o5(1),c6_n2o5(2),v3_du2)
!
!             rjx(i,j,jn2o5) =  sig_1 * fint_1 &
!                  +  sig_2 * fint_2 &
!                  +  sig_3 * fint_3 &
!                  +  sig_4 * fint_4 &
!                  +  sig_5 * fint_5 &
!                  +  sig_6 * fint_6 
!
!             !--------------------
!             ! ch2o -> co + h2
!             !--------------------
!
!             sig_3 = p1(tj_coh2(1,3),tj_coh2(2,3),tscale) * &
!                  p1(b3_coh2(i3),a3_coh2(i3),v3_du2)
!             sig_4 = p3(c4_coh2(1),c4_coh2(2),c4_coh2(3),c4_coh2(4),v3_du2)
!             sig_5 = p1(tj_coh2(1,5),tj_coh2(2,5),tscale) * &
!                  p2(c5_coh2(1),c5_coh2(2),c5_coh2(3),v3_du2)
!             cc0 = p2(      1.00997 ,  -0.398144 ,  -0.265157, tscale)
!             cc1 = p2(     0.387158 ,  -0.434852 ,  -0.221940, tscale)
!             cc2 = p2(   -0.0388784 ,  0.0571797 ,  0.0876116, tscale)
!             sig_6 = p1(c6_coh2(1),c6_coh2(2),v3_du2)/ p2(cc0,cc1,cc2,pscale)
!
!             rjx(i,j,jach2o) =  sig_3 * fint_3 &
!                  +  sig_4 * fint_4 &
!                  +  sig_5 * fint_5 &
!                  +  sig_6 * fint_6 
!
!             !--------------------
!             ! ch2o -> cho + h
!             !--------------------
!             sig_3 = p1(tj_choh(1,3),tj_choh(2,3),tscale) * &
!                  p1(b3_choh(i3),a3_choh(i3),v3_du2)
!             sig_4 = p3(c4_choh(1),c4_choh(2),c4_choh(3), c4_choh(4),v3_du2)
!             sig_5 = p1(tj_choh(1,5),tj_choh(2,5),tscale) * &
!                  p2(c5_choh(1),c5_choh(2),c5_choh(3),v3_du2)
!             sig_6 = p1(tj_choh(1,6),tj_choh(2,6),tscale) * &
!                  p1(c6_choh(1),c6_choh(2),v3_du2)
!
!             rjx(i,j,jbch2o) =  sig_3 * fint_3 &
!                  +  sig_4 * fint_4 &
!                  +  sig_5 * fint_5 &
!                  +  sig_6 * fint_6 
!
!             !--------------------
!             ! ch3ooh --> oh + ch3o
!             !--------------------
!             sig_1 =p1(b1_ch3ooh(i1,1),a1_ch3ooh(i1,1),v3_du1) + &
!                  p1(b1_ch3ooh(i1,2),a1_ch3ooh(i1,2),v3_du1) * v2s_m
!             sig_2 =p1(b2_ch3ooh(i2),a2_ch3ooh(i2),v3_du1)
!             sig_3 =p1(b3_ch3ooh(i3),a3_ch3ooh(i3),v3_du2)
!             sig_4 =p2(c4_ch3ooh(1),c4_ch3ooh(2),c4_ch3ooh(3),v3_du2)
!             sig_5 =p2(c5_ch3ooh(1),c5_ch3ooh(2),c5_ch3ooh(3),v3_du2)
!             sig_6 =p1(c6_ch3ooh(1),c6_ch3ooh(2),v3_du2)
!
!             rjx(i,j,jmepe) =  sig_1 * fint_1 &
!                  +  sig_2 * fint_2 &
!                  +  sig_3 * fint_3 &
!                  +  sig_4 * fint_4 &
!                  +  sig_5 * fint_5 &
!                  +  sig_6 * fint_6 
!
!
!             !--------------------
!             ! no3 -> no2 + o
!             !--------------------
!
!
!             sig_6 = p1(tj_no2o(1,6),tj_no2o(2,6),tscale) * &
!                  p1(c6_no2o(1),c6_no2o(2),v3_du2)
!             sig_7 = p1(tj_no2o(1,7),tj_no2o(2,7),tscale) * &
!                  p1(c7_no2o(1),c7_no2o(2),v3_du2)
!
!             rjx(i,j,jano3)  = sig_6 * fint_6 + &
!                  sig_7 * fint_7 
!
!             !--------------------
!             ! no3 -> no + o2
!             !--------------------
!
!             sig_7 =p1(tj_noo2(1,7),tj_noo2(2,7),tscale)* &
!                  p1(c7_noo2(1),c7_noo2(2),v3_du2)
!
!             rjx(i,j,jbno3) = sig_7 * fint_7 
!
!             ! ---
!             ! pan
!             ! ---
!
!             sig_1 = p3(tj_pan(1,1),tj_pan(2,1),tj_pan(3,1),tj_pan(4,1),tscale)* &
!                  (p1(b1_pan(i1,1),a1_pan(i1,1), v3_du1) + &
!                  p1(b1_pan(i1,2),a1_pan(i1,2), v3_du1) *v2s_m)
!             sig_2 = p3(tj_pan(1,2),tj_pan(2,2),tj_pan(3,2),tj_pan(4,2),tscale)* &
!                  p1(b2_pan(i2),a2_pan(i2),v3_du1)
!             sig_3 = p3(tj_pan(1,3),tj_pan(2,3),tj_pan(3,3),tj_pan(4,3),tscale)* &
!                  p1(b3_pan(i3),a3_pan(i3),v3_du2)
!             sig_4 = p3(tj_pan(1,4),tj_pan(2,4),tj_pan(3,4),tj_pan(4,4),tscale)* &
!                  p2(c4_pan(1),c4_pan(2),c4_pan(3),v3_du2)
!             sig_5 = p3(tj_pan(1,5),tj_pan(2,5),tj_pan(3,5),tj_pan(4,5),tscale)* &
!                  p2(c5_pan(1),c5_pan(2),c5_pan(3),v3_du2)
!             sig_6 = p3(tj_pan(1,6),tj_pan(2,6),tj_pan(3,6),tj_pan(4,6),tscale)* &
!                  p1(c6_pan(1),c6_pan(2),v3_du2)
!
!             rjx(i,j,jpan) =  sig_1 * fint_1 &
!                  +  sig_2 * fint_2 &
!                  +  sig_3 * fint_3 &
!                  +  sig_4 * fint_4 &
!                  +  sig_5 * fint_5 &
!                  +  sig_6 * fint_6 
!
!             ! ---
!             ! ch3ono2
!             ! ---
!
!             sig_2 = p1(b2_ch3ono2(i2),a2_ch3ono2(i2),v3_du1)
!             sig_3 = p1(b3_ch3ono2(i2),a3_ch3ono2(i2),v3_du2)
!             sig_4 = p2(c4_ch3ono2(1), c4_ch3ono2(2), c4_ch3ono2(3), v3_du2)
!             sig_5 = p2(c5_ch3ono2(1), c5_ch3ono2(2), c5_ch3ono2(3), v3_du2)
!
!             rjx(i,j,jorgn) = sig_2 * fint_2 + &
!                  sig_3 * fint_3 + &
!                  sig_4 * fint_4 + &
!                  sig_5 * fint_5
!
!             ! ---
!             ! ch3cho
!             ! ---
!
!             sig_2 = p1(b2_ch3cho(i2),a2_ch3cho(i2),v3_du1)               &
!                  / p1(tp_ch3cho(1,2),tp_ch3cho(2,2),pscale)
!             sig_3 = p1(b3_ch3cho(i2),a3_ch3cho(i2),v3_du2)               &
!                  / p1(tp_ch3cho(1,3),tp_ch3cho(2,3),pscale)
!             sig_4 = p2(c4_ch3cho(1), c4_ch3cho(2), c4_ch3cho(3), v3_du2) &
!                  / p1(tp_ch3cho(1,4),tp_ch3cho(2,4),pscale)
!             sig_5 = p2(c5_ch3cho(1), c5_ch3cho(2), c5_ch3cho(3), v3_du2) &
!                  / p1(tp_ch3cho(1,5),tp_ch3cho(2,5),pscale)
!
!             rjx(i,j,j45) = sig_2 * fint_2 + &
!                  sig_3 * fint_3 + &
!                  sig_4 * fint_4 + &
!                  sig_5 * fint_5
!
!             ! multiply methyl-nitrate to get organic nitrate
!             rjx(i,j,jorgn) = 2.5 * rjx(i,j,jorgn)
!             rjx(i,j,j74)   = 50.0* rjx(i,j,jbch2o)  !methylglyoxal
!             rjx(i,j,jrooh) = rjx(i,j,jmepe)         !other peroxides...
!          end if
!
!       end do
!    end do   ! i,j loop
!
!    ! correction photolysis for dsitance sun-earth
!    esrm2 = sundis(idate(2),idate(3))
!    ! multiply with seasonal correction factor
!    rjx  = rjx*esrm2  
!
!    phot_dat(region)%jo3_av (:,:,offsetl+level) = &
!         phot_dat(region)%jo3_av (:,:,offsetl+level) + rjx(:,:,jo3d)
!    phot_dat(region)%jno2_av(:,:,offsetl+level) = &
!         phot_dat(region)%jno2_av(:,:,offsetl+level) + rjx(:,:,jno2)
!#ifdef MPI
!    if ( level == lmloc ) phot_dat(region)%nj_av = &
!         phot_dat(region)%nj_av + 1   !proceed only after last level....
!#else
!    if ( level == lm(region) ) phot_dat(region)%nj_av = &
!         phot_dat(region)%nj_av + 1   !proceed only after last level....
!#endif
!
!    nullify(phlb)
!    nullify(zoomed)
!    nullify(t)
!    nullify(vo3)
!
!  contains
!
!    real function p1(c0,c1,x)
!      implicit none
!      real, intent(in) :: c0, c1, x
!      p1 = c0 + c1*x
!    end function p1
!
!    real function p2(c0,c1,c2,x)
!      implicit none
!      real, intent(in) :: c0, c1, c2, x
!      p2 = c0 + (c1 + c2*x)*x
!    end function p2
!
!    real function p3(c0,c1,c2,c3,x)
!      implicit none
!      real, intent(in) :: c0, c1, c2, c3, x
!      p3 = c0 + (c1 + (c2 + c3*x)*x)*x
!    end function p3
!
!  end subroutine photolysis_rates


  !
  ! subroutine that deallocates photolysis memory 
  ! and writes photolysis statistics
  ! called from trace_end
  !

  subroutine end_photolysis

    use dims           , only : nregions
    use photolysis_data, only : phot_dat
!    use global_data, only : outdir
!    use io_hdf
!#ifdef MPI
!    use mpi_const,   only : lmloc, my_real, myid, root_k, mpi_sum, com_lev, ierr
!#endif
!
!    ! called by trace_end
!    ! output of (e.g. monthly) averaged photolysis jo3 jno2
!    ! and surface data of: albedo, cloud-cover, cloud optical thickness
!    ! total ozone
!    ! Deallocation of memory
!
!    integer   :: io,n ,sfstart,sfend,sfsattr
!    integer   :: k,iret,imr,jmr,lmr
!    character(len=22) :: file_name
!    character(len=7) :: rg
    integer  :: region
!    real,dimension(:,:,:),allocatable    :: average_field
!
!    file_name = trim(outdir)//'/j_statistics.hdf'
!#ifdef MPI
!    if (myid == root_k) io = sfstart(file_name,DFACC_CREATE)
!#else
!    io = sfstart(file_name,DFACC_CREATE)
!#endif
!
!    do region=1,nregions
!       rg = '_'//trim(region_name(region))
!       imr=im(region)
!       jmr=jm(region)
!       lmr=lm(region)          
!
!       print *, 'end_photolysis: j_statistics ', &
!            '(region, nj_av, nvo3s_av, ncloud_av)', region,   &
!            phot_dat(region)%nj_av, phot_dat(region)%nvo3s_av, &
!            phot_dat(region)%ncloud_av
!       if ( phot_dat(region)%nj_av > 0 ) then
!          allocate(average_field(imr,jmr,lmr))
!#ifdef MPI          
!          if ( lmloc /= 0 ) then
!             call mpi_allreduce(phot_dat(region)%jo3_av,average_field, &
!                  imr*jmr*lmr,my_real,mpi_sum,com_lev,ierr)
!          end if
!#else
!          average_field = phot_dat(region)%jo3_av
!#endif
!#ifdef MPI          
!          if ( myid == root_k ) then
!#endif
!             average_field = average_field/phot_dat(region)%nj_av
!             call io_write(io,imr,'LON'//rg,jmr,'LAT'//rg, &
!                  lmr,'HYBRID'//rg,average_field,'jo3_av')
!#ifdef MPI          
!          end if
!#endif
!#ifdef MPI          
!          if ( lmloc /= 0 ) then
!             call mpi_allreduce(phot_dat(region)%jno2_av,average_field, &
!                  imr*jmr*lmr,my_real,mpi_sum,com_lev,ierr)
!          end if
!#else
!          average_field = phot_dat(region)%jno2_av
!
!#endif
!#ifdef MPI          
!          if ( myid == root_k ) then
!#endif
!             average_field = average_field/phot_dat(region)%nj_av
!             call io_write(io,imr,'LON'//rg,jmr,'LAT'//rg, &
!                  lmr,'HYBRID'//rg,average_field,'jno2_av')
!#ifdef MPI          
!          end if
!#endif
!          deallocate(average_field)
!       end if
!#ifdef MPI     
!       if(myid == root_k) then
!#endif
!          if (phot_dat(region)%nvo3s_av > 0) then
!             phot_dat(region)%vo3s_av = &
!                  todu*phot_dat(region)%vo3s_av/phot_dat(region)%nvo3s_av
!             call io_write(io,imr,'LON'//rg,jmr,'LAT'//rg, &
!                  phot_dat(region)%vo3s_av,'vo3s_av')
!          end if
!          if ( phot_dat(region)%nalb_av > 0 ) then
!             phot_dat(region)%ags_av = &
!                  phot_dat(region)%ags_av/phot_dat(region)%nalb_av
!             call io_write(io,imr,'LON'//rg,jmr,'LAT'//rg, &
!                  phot_dat(region)%ags_av,'ags_av')
!          end if
!#ifdef MPI     
!       end if
!#endif
!
!#ifdef MPI     
!       if ( myid == root_k ) then
!#endif
!          if ( phot_dat(region)%ncloud_av > 0 ) then
!             phot_dat(region)%taus_av = &
!                  phot_dat(region)%taus_av/phot_dat(region)%ncloud_av
!             call io_write(io,imr,'LON'//rg,jmr,'LAT'//rg, &
!                  phot_dat(region)%taus_av,'taus_av')
!             phot_dat(region)%cfr_av = &
!                  phot_dat(region)%cfr_av/phot_dat(region)%ncloud_av
!             call io_write(io,imr,'LON'//rg,jmr,'LAT'//rg, &
!                  phot_dat(region)%cfr_av,'cfr_av')
!          end if
!#ifdef MPI     
!       end if
!#endif
!    end do
!
!#ifdef MPI
!    if ( myid == root_k ) iret = sfend(io)
!#else
!    iret = sfend(io)
!#endif

    do region = 1,nregions
!       deallocate(phot_dat(region)%taus      )
!       deallocate(phot_dat(region)%cloudfrac )
       deallocate(phot_dat(region)%albedo    )
!       deallocate(phot_dat(region)%vo3s_av   )
!       deallocate(phot_dat(region)%taus_av   )
!       deallocate(phot_dat(region)%cfr_av    )
!       deallocate(phot_dat(region)%vo3       )
!       deallocate(phot_dat(region)%tauc      )
!       deallocate(phot_dat(region)%jo3_av    )
!       deallocate(phot_dat(region)%jno2_av   )
!JEW:   deallocate(phot_dat(region)%ags_av    )
    end do 

  end subroutine  end_photolysis


end module photolysis
