Index: src/ohreaction.f90
===================================================================
--- src/ohreaction.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ src/ohreaction.f90	(revision 8a65cb0801b09eb0ae933a1afabbe10361a9021b)
@@ -25,19 +25,19 @@
   !                                                                            *
   !                                                                            *
-  !    Author: S. Eckhardt                                                     *
+  !    Author: R.L. Thompson                                                   *
   !                                                                            *
-  !    June 2007                                                               *
+  !    Nov 2014                                                                *
   !                                                                            *
   !                                                                            *
   !*****************************************************************************
   ! Variables:                                                                 *
-  ! ix,jy              indices of output grid cell for each particle           *
-  ! itime [s]          actual simulation time [s]                              *
-  ! jpart              particle index                                          *
-  ! ldeltat [s]        interval since radioactive decay was computed           *
-  ! loutnext [s]       time for which gridded deposition is next output        *
-  ! loutstep [s]       interval at which gridded deposition is output          *
-  ! oh_average [mol/m^3]   OH Concentration                                    *
-  ! ltsample [s]       interval over which mass is deposited                   *
+  ! ix,jy                indices of output grid cell for each particle         *
+  ! itime [s]            actual simulation time [s]                            *
+  ! jpart                particle index                                        *
+  ! ldeltat [s]          interval since radioactive decay was computed         *
+  ! loutnext [s]         time for which gridded deposition is next output      *
+  ! loutstep [s]         interval at which gridded deposition is output        *
+  ! oh_average [molecule/cm^3] OH Concentration                                *
+  ! ltsample [s]         interval over which mass is deposited                 *
   !                                                                            *
   !*****************************************************************************
@@ -49,14 +49,13 @@
   implicit none
 
-  integer :: jpart,itime,ltsample,loutnext,ldeltat,j,k,ix,jy
-  integer :: ngrid,il,interp_time,n,mm,indz,i
-  integer :: jjjjmmdd,ihmmss,OHx,OHy,dOHx,dOHy,OHz
-  real :: xtn,ytn,oh_average
-  !real oh_diurn_var,sum_ang
-  !real zenithangle, ang
-  real :: restmass,ohreacted,OHinc
-  real :: xlon, ylat, gas_const, act_energy
-  real :: ohreact_temp_corr, act_temp
-  real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
+  integer :: jpart,itime,ltsample,loutnext,ldeltat,j,k,ix,jy!,ijx,jjy
+  integer :: ngrid,interp_time,n,m,h,indz,i!,ia,il
+  integer :: jjjjmmdd,hhmmss,OHx,OHy,OHz
+  real, dimension(nzOH) :: altOHtop
+  real :: xlon,ylat 
+  real :: xtn,ytn
+  real :: restmass,ohreacted,oh_average
+  real :: ohrate,temp 
+  real, parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
   real(kind=dp) :: jul
 
@@ -64,8 +63,4 @@
   !************************************************************************
 
-  gas_const=8.314 ! define gas constant
-  act_energy=10000 ! activation energy
-
-  !write(*,*) 'OH reaction n:',n,ohreact(1)
   if (itime.le.loutnext) then
     ldeltat=itime-(loutnext-loutstep)
@@ -74,17 +69,15 @@
   endif
 
+  jul=bdate+real(itime,kind=dp)/86400.
+  call caldate(jul,jjjjmmdd,hhmmss)
+  m=(jjjjmmdd-(jjjjmmdd/10000)*10000)/100
+  h=hhmmss/10000
 
-    dOHx=360/(maxxOH-1)
-    dOHy=180/(maxyOH-1)
-
-    jul=bdate+real(itime,kind=dp)/86400._dp
-    call caldate(jul,jjjjmmdd,ihmmss)
-    mm=int((jjjjmmdd-(jjjjmmdd/10000)*10000)/100)
-
-    do jpart=1,numpart
-
-  ! Determine which nesting level to be used
+  ! Loop over particles
   !*****************************************
 
+  do jpart=1,numpart
+
+    ! Determine which nesting level to be used
     ngrid=0
     do j=numbnests,1,-1
@@ -95,10 +88,7 @@
       endif
     end do
-23   continue
+23  continue
 
-
-  ! Determine nested grid coordinates
-  !**********************************
-
+    ! Determine nested grid coordinates
     if (ngrid.gt.0) then
       xtn=(xtra1(jpart)-xln(ngrid))*xresoln(ngrid)
@@ -111,28 +101,20 @@
     endif
 
-  n=2
-  if (abs(memtime(1)-interp_time).lt.abs(memtime(2)-interp_time)) &
-       n=1
+    interp_time=nint(itime-0.5*ltsample)
+    n=2
+    if(abs(memtime(1)-interp_time).lt.abs(memtime(2)-interp_time)) n=1
 
-  do i=2,nz
-    if (height(i).gt.ztra1(jpart)) then
-      indz=i-1
-      goto 6
-    endif
-  end do
+    do i=2,nz
+      if (height(i).gt.ztra1(jpart)) then
+        indz=i-1
+        goto 6
+      endif
+    end do
 6   continue
 
-  ! The concentration from the nearest available gridcell is taken
-  ! get OH concentration for the specific month and solar angle
+    ! Get OH from nearest grid-cell and specific month 
+    !*************************************************
 
-  !  write(*,*) OH_field(1,1,1,1),OH_field(10,1,1,10)
-  !  write(*,*) OH_field(1,maxxOH-1,maxyOH-1,1)
-  !  write(*,*) OH_field(10,maxxOH-1,maxyOH-1,10)
-  !  write(*,*) OH_field_height(1,10,4,1),OH_field_height(10,4,10,10)
-  !  write(*,*) OH_field_height(1,maxxOH-1,maxyOH-1,1)
-  !  write(*,*) OH_field_height(10,maxxOH-1,maxyOH-1,10)
-    interp_time=nint(itime-0.5*ltsample)
-
-  ! World coordinates
+    ! world coordinates
     xlon=xtra1(jpart)*dx+xlon0
     if (xlon.gt.180) then
@@ -140,74 +122,54 @@
     endif
     ylat=ytra1(jpart)*dy+ylat0
-  ! get position in the OH field - assume that the OH field is global
-    OHx=(180+xlon-1)/dOHx
-    OHy=(90+ylat-1)/dOHy
-  !  sum_ang=0
-  ! get the level of the OH height field were the actual particle is in
-  ! ztra1 is the z-coordinate of the trajectory above model orography in m
-  ! OH_field_height is the heigth of the OH field above orography
-      OHz=maxzOH
-  ! assume equally distrib. OH field, OH_field_height gives the middle of
-  ! the z coordinate
-      OHinc=(OH_field_height(3)-OH_field_height(2))/2
-      do il=2,maxzOH+1
-        if ((OH_field_height(il-1)+OHinc).gt.ztra1(jpart)) goto 26
+
+    ! get position in the OH field
+    OHx=minloc(abs(lonOH-xlon),dim=1,mask=abs(lonOH-xlon).eq.minval(abs(lonOH-xlon)))
+    OHy=minloc(abs(latOH-ylat),dim=1,mask=abs(latOH-ylat).eq.minval(abs(latOH-ylat)))
+
+    ! get the level of the OH field for the particle
+    ! ztra1 is the z-coord of the trajectory above model orography in metres
+    ! altOH is the height of the centre of the level in the OH field above orography
+    do i=2,nzOH
+      altOHtop(i-1)=altOH(i)+0.5*(altOH(i)-altOH(i-1))
+    end do
+    altOHtop(nzOH)=altOH(nzOH)+0.5*(altOH(nzOH)-altOH(nzOH-1))
+    OHz=minloc(abs(altOHtop-ztra1(jpart)),dim=1,mask=abs(altOHtop-ztra1(jpart))&
+            &.eq.minval(abs(altOHtop-ztra1(jpart))))
+
+    ! Interpolate between hourly OH fields to current time
+    !*****************************************************
+
+    oh_average=OH_hourly(OHx,OHy,OHz,1)+&
+               &(OH_hourly(OHx,OHy,OHz,2)-OH_hourly(OHx,OHy,OHz,1))*&
+               &(itime-memOHtime(1))/(memOHtime(2)-memOHtime(1))
+
+    if (oh_average.gt.smallnum) then
+
+      ! Computation of the OH reaction
+      !**********************************************************
+
+      temp=tt(ix,jy,indz,n)
+
+      do k=1,nspec                                
+        if (ohcconst(k).gt.0.) then
+          ohrate=ohcconst(k)*temp**2*exp(-ohdconst(k)/temp)*oh_average
+          ! new particle mass
+          restmass = xmass1(jpart,k)*exp(-1*ohrate*abs(ltsample))
+          if (restmass .gt. smallnum) then
+            xmass1(jpart,k)=restmass
+          else
+            xmass1(jpart,k)=0.
+          endif
+          ohreacted=xmass1(jpart,k)*(1-exp(-1*ohrate*abs(ltsample)))
+        else
+          ohreacted=0.
+        endif
       end do
-26     continue
 
-     OHz=il-1
-  !   loop was not interrupted il would be 8 (9-1)
-     if (OHz.gt.maxzOH) OHz=7
-  !   write (*,*) 'OH height: '
-  !    +        ,ztra1(jpart),jpart,OHz,OH_field_height(OHz),OHinc,
-  !    +        OH_field_height
+    endif  ! oh_average.gt.smallnum 
 
-    oh_average=OH_field(mm,OHx,OHy,OHz)
-    if (oh_average.gt.smallnum) then
-  !**********************************************************
-  ! if there is noOH concentration no reaction
-  ! for performance reason take average concentration and
-  ! ignore diurnal variation
-  ! do 28 il=1,24
-  !      ang=70-zenithangle(ylat,xlon,jul+(24-il)/24.)
-  !      if (ang.lt.0) then
-  !          ang=0
-  !      endif
-  !      sum_ang=sum_ang+ang
-  !28         enddo
-  !    oh_diurn_var=(ang/sum_ang)*(oh_average*24)
-  !    oh_average=oh_diurn_var
-  !**********************************************************
+  end do  !continue loop over all particles
 
 
-  !    Computation of the OH reaction
-  !**********************************************************
-    act_temp=tt(ix,jy,indz,n)
+end subroutine ohreaction
 
-    do k=1,nspec                                  ! loop over species
-      if (ohreact(k).gt.0.) then
-        ohreact_temp_corr=ohreact(k)*oh_average* &
-             exp((act_energy/gas_const)*(1/298.15-1/act_temp))
-        ohreacted=xmass1(jpart,k)* &
-             (1.-exp(-1*ohreact_temp_corr*abs(ltsample)))
-  !      new particle mass:
-        restmass = xmass1(jpart,k)-ohreacted
-        if (restmass .gt. smallnum) then
-          xmass1(jpart,k)=restmass
-  !   write (104) xlon,ylat,ztra1(jpart),k,oh_diurn_var,jjjjmmdd,
-  !    +               ihmmss,restmass,ohreacted
-        else
-          xmass1(jpart,k)=0.
-        endif
-  !      write (*,*) 'restmass: ',restmass
-      else
-        ohreacted=0.
-      endif
-    end do
-
-  endif
-  !endif OH concentration gt 0
-    end do
-  !continue loop over all particles
-
-end subroutine ohreaction
