| 1 | !********************************************************************** |
|---|
| 2 | ! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 * |
|---|
| 3 | ! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, * |
|---|
| 4 | ! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann * |
|---|
| 5 | ! * |
|---|
| 6 | ! This file is part of FLEXPART. * |
|---|
| 7 | ! * |
|---|
| 8 | ! FLEXPART is free software: you can redistribute it and/or modify * |
|---|
| 9 | ! it under the terms of the GNU General Public License as published by* |
|---|
| 10 | ! the Free Software Foundation, either version 3 of the License, or * |
|---|
| 11 | ! (at your option) any later version. * |
|---|
| 12 | ! * |
|---|
| 13 | ! FLEXPART is distributed in the hope that it will be useful, * |
|---|
| 14 | ! but WITHOUT ANY WARRANTY; without even the implied warranty of * |
|---|
| 15 | ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * |
|---|
| 16 | ! GNU General Public License for more details. * |
|---|
| 17 | ! * |
|---|
| 18 | ! You should have received a copy of the GNU General Public License * |
|---|
| 19 | ! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. * |
|---|
| 20 | !********************************************************************** |
|---|
| 21 | |
|---|
| 22 | subroutine ohreaction(itime,ltsample,loutnext) |
|---|
| 23 | ! i i i |
|---|
| 24 | !***************************************************************************** |
|---|
| 25 | ! * |
|---|
| 26 | ! * |
|---|
| 27 | ! Author: S. Eckhardt * |
|---|
| 28 | ! * |
|---|
| 29 | ! June 2007 * |
|---|
| 30 | ! * |
|---|
| 31 | ! * |
|---|
| 32 | !***************************************************************************** |
|---|
| 33 | ! Variables: * |
|---|
| 34 | ! ix,jy indices of output grid cell for each particle * |
|---|
| 35 | ! itime [s] actual simulation time [s] * |
|---|
| 36 | ! jpart particle index * |
|---|
| 37 | ! ldeltat [s] interval since radioactive decay was computed * |
|---|
| 38 | ! loutnext [s] time for which gridded deposition is next output * |
|---|
| 39 | ! loutstep [s] interval at which gridded deposition is output * |
|---|
| 40 | ! oh_average [mol/m^3] OH Concentration * |
|---|
| 41 | ! ltsample [s] interval over which mass is deposited * |
|---|
| 42 | ! * |
|---|
| 43 | !***************************************************************************** |
|---|
| 44 | |
|---|
| 45 | use oh_mod |
|---|
| 46 | use par_mod |
|---|
| 47 | use com_mod |
|---|
| 48 | |
|---|
| 49 | implicit none |
|---|
| 50 | |
|---|
| 51 | integer :: jpart,itime,ltsample,loutnext,ldeltat,j,k,ix,jy |
|---|
| 52 | integer :: ngrid,il,interp_time,n,mm,indz,i |
|---|
| 53 | integer :: jjjjmmdd,ihmmss,OHx,OHy,dOHx,dOHy,OHz |
|---|
| 54 | real :: xtn,ytn,oh_average |
|---|
| 55 | !real oh_diurn_var,sum_ang |
|---|
| 56 | !real zenithangle, ang |
|---|
| 57 | real :: restmass,ohreacted,OHinc |
|---|
| 58 | real :: xlon, ylat, gas_const, act_energy |
|---|
| 59 | real :: ohreact_temp_corr, act_temp |
|---|
| 60 | real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled |
|---|
| 61 | real(kind=dp) :: jul |
|---|
| 62 | |
|---|
| 63 | ! Compute interval since radioactive decay of deposited mass was computed |
|---|
| 64 | !************************************************************************ |
|---|
| 65 | |
|---|
| 66 | gas_const=8.314 ! define gas constant |
|---|
| 67 | act_energy=10000 ! activation energy |
|---|
| 68 | |
|---|
| 69 | !write(*,*) 'OH reaction n:',n,ohreact(1) |
|---|
| 70 | if (itime.le.loutnext) then |
|---|
| 71 | ldeltat=itime-(loutnext-loutstep) |
|---|
| 72 | else ! first half of next interval |
|---|
| 73 | ldeltat=itime-loutnext |
|---|
| 74 | endif |
|---|
| 75 | |
|---|
| 76 | |
|---|
| 77 | dOHx=360/(maxxOH-1) |
|---|
| 78 | dOHy=180/(maxyOH-1) |
|---|
| 79 | |
|---|
| 80 | jul=bdate+real(itime,kind=dp)/86400._dp |
|---|
| 81 | call caldate(jul,jjjjmmdd,ihmmss) |
|---|
| 82 | mm=int((jjjjmmdd-(jjjjmmdd/10000)*10000)/100) |
|---|
| 83 | |
|---|
| 84 | do jpart=1,numpart |
|---|
| 85 | |
|---|
| 86 | ! Determine which nesting level to be used |
|---|
| 87 | !***************************************** |
|---|
| 88 | |
|---|
| 89 | ngrid=0 |
|---|
| 90 | do j=numbnests,1,-1 |
|---|
| 91 | if ((xtra1(jpart).gt.xln(j)).and.(xtra1(jpart).lt.xrn(j)).and. & |
|---|
| 92 | (ytra1(jpart).gt.yln(j)).and.(ytra1(jpart).lt.yrn(j))) then |
|---|
| 93 | ngrid=j |
|---|
| 94 | goto 23 |
|---|
| 95 | endif |
|---|
| 96 | end do |
|---|
| 97 | 23 continue |
|---|
| 98 | |
|---|
| 99 | |
|---|
| 100 | ! Determine nested grid coordinates |
|---|
| 101 | !********************************** |
|---|
| 102 | |
|---|
| 103 | if (ngrid.gt.0) then |
|---|
| 104 | xtn=(xtra1(jpart)-xln(ngrid))*xresoln(ngrid) |
|---|
| 105 | ytn=(ytra1(jpart)-yln(ngrid))*yresoln(ngrid) |
|---|
| 106 | ix=int(xtn) |
|---|
| 107 | jy=int(ytn) |
|---|
| 108 | else |
|---|
| 109 | ix=int(xtra1(jpart)) |
|---|
| 110 | jy=int(ytra1(jpart)) |
|---|
| 111 | endif |
|---|
| 112 | |
|---|
| 113 | n=2 |
|---|
| 114 | if (abs(memtime(1)-interp_time).lt.abs(memtime(2)-interp_time)) & |
|---|
| 115 | n=1 |
|---|
| 116 | |
|---|
| 117 | do i=2,nz |
|---|
| 118 | if (height(i).gt.ztra1(jpart)) then |
|---|
| 119 | indz=i-1 |
|---|
| 120 | goto 6 |
|---|
| 121 | endif |
|---|
| 122 | end do |
|---|
| 123 | 6 continue |
|---|
| 124 | |
|---|
| 125 | ! The concentration from the nearest available gridcell is taken |
|---|
| 126 | ! get OH concentration for the specific month and solar angle |
|---|
| 127 | |
|---|
| 128 | ! write(*,*) OH_field(1,1,1,1),OH_field(10,1,1,10) |
|---|
| 129 | ! write(*,*) OH_field(1,maxxOH-1,maxyOH-1,1) |
|---|
| 130 | ! write(*,*) OH_field(10,maxxOH-1,maxyOH-1,10) |
|---|
| 131 | ! write(*,*) OH_field_height(1,10,4,1),OH_field_height(10,4,10,10) |
|---|
| 132 | ! write(*,*) OH_field_height(1,maxxOH-1,maxyOH-1,1) |
|---|
| 133 | ! write(*,*) OH_field_height(10,maxxOH-1,maxyOH-1,10) |
|---|
| 134 | interp_time=nint(itime-0.5*ltsample) |
|---|
| 135 | |
|---|
| 136 | ! World coordinates |
|---|
| 137 | xlon=xtra1(jpart)*dx+xlon0 |
|---|
| 138 | if (xlon.gt.180) then |
|---|
| 139 | xlon=xlon-360 |
|---|
| 140 | endif |
|---|
| 141 | ylat=ytra1(jpart)*dy+ylat0 |
|---|
| 142 | ! get position in the OH field - assume that the OH field is global |
|---|
| 143 | OHx=(180+xlon-1)/dOHx |
|---|
| 144 | OHy=(90+ylat-1)/dOHy |
|---|
| 145 | ! sum_ang=0 |
|---|
| 146 | ! get the level of the OH height field were the actual particle is in |
|---|
| 147 | ! ztra1 is the z-coordinate of the trajectory above model orography in m |
|---|
| 148 | ! OH_field_height is the heigth of the OH field above orography |
|---|
| 149 | OHz=maxzOH |
|---|
| 150 | ! assume equally distrib. OH field, OH_field_height gives the middle of |
|---|
| 151 | ! the z coordinate |
|---|
| 152 | OHinc=(OH_field_height(3)-OH_field_height(2))/2 |
|---|
| 153 | do il=2,maxzOH+1 |
|---|
| 154 | if ((OH_field_height(il-1)+OHinc).gt.ztra1(jpart)) goto 26 |
|---|
| 155 | end do |
|---|
| 156 | 26 continue |
|---|
| 157 | |
|---|
| 158 | OHz=il-1 |
|---|
| 159 | ! loop was not interrupted il would be 8 (9-1) |
|---|
| 160 | if (OHz.gt.maxzOH) OHz=7 |
|---|
| 161 | ! write (*,*) 'OH height: ' |
|---|
| 162 | ! + ,ztra1(jpart),jpart,OHz,OH_field_height(OHz),OHinc, |
|---|
| 163 | ! + OH_field_height |
|---|
| 164 | |
|---|
| 165 | oh_average=OH_field(mm,OHx,OHy,OHz) |
|---|
| 166 | if (oh_average.gt.smallnum) then |
|---|
| 167 | !********************************************************** |
|---|
| 168 | ! if there is noOH concentration no reaction |
|---|
| 169 | ! for performance reason take average concentration and |
|---|
| 170 | ! ignore diurnal variation |
|---|
| 171 | ! do 28 il=1,24 |
|---|
| 172 | ! ang=70-zenithangle(ylat,xlon,jul+(24-il)/24.) |
|---|
| 173 | ! if (ang.lt.0) then |
|---|
| 174 | ! ang=0 |
|---|
| 175 | ! endif |
|---|
| 176 | ! sum_ang=sum_ang+ang |
|---|
| 177 | !28 enddo |
|---|
| 178 | ! oh_diurn_var=(ang/sum_ang)*(oh_average*24) |
|---|
| 179 | ! oh_average=oh_diurn_var |
|---|
| 180 | !********************************************************** |
|---|
| 181 | |
|---|
| 182 | |
|---|
| 183 | ! Computation of the OH reaction |
|---|
| 184 | !********************************************************** |
|---|
| 185 | act_temp=tt(ix,jy,indz,n) |
|---|
| 186 | |
|---|
| 187 | do k=1,nspec ! loop over species |
|---|
| 188 | if (ohreact(k).gt.0.) then |
|---|
| 189 | ohreact_temp_corr=ohreact(k)*oh_average* & |
|---|
| 190 | exp((act_energy/gas_const)*(1/298.15-1/act_temp)) |
|---|
| 191 | ohreacted=xmass1(jpart,k)* & |
|---|
| 192 | (1.-exp(-1*ohreact_temp_corr*abs(ltsample))) |
|---|
| 193 | ! new particle mass: |
|---|
| 194 | restmass = xmass1(jpart,k)-ohreacted |
|---|
| 195 | if (restmass .gt. smallnum) then |
|---|
| 196 | xmass1(jpart,k)=restmass |
|---|
| 197 | ! write (104) xlon,ylat,ztra1(jpart),k,oh_diurn_var,jjjjmmdd, |
|---|
| 198 | ! + ihmmss,restmass,ohreacted |
|---|
| 199 | else |
|---|
| 200 | xmass1(jpart,k)=0. |
|---|
| 201 | endif |
|---|
| 202 | ! write (*,*) 'restmass: ',restmass |
|---|
| 203 | else |
|---|
| 204 | ohreacted=0. |
|---|
| 205 | endif |
|---|
| 206 | end do |
|---|
| 207 | |
|---|
| 208 | endif |
|---|
| 209 | !endif OH concentration gt 0 |
|---|
| 210 | end do |
|---|
| 211 | !continue loop over all particles |
|---|
| 212 | |
|---|
| 213 | end subroutine ohreaction |
|---|
