ohreaction.f90 5.49 KB
Newer Older
Matthias Langer's avatar
 
Matthias Langer committed
1
2
3
4
5
subroutine ohreaction(itime,ltsample,loutnext)
  !                     i      i        i
  !*****************************************************************************
  !                                                                            *
  !                                                                            *
6
  !    Author: R.L. Thompson                                                   *
Matthias Langer's avatar
 
Matthias Langer committed
7
  !                                                                            *
8
  !    Nov 2014                                                                *
Matthias Langer's avatar
 
Matthias Langer committed
9
10
11
12
  !                                                                            *
  !                                                                            *
  !*****************************************************************************
  ! Variables:                                                                 *
13
14
15
16
17
18
19
20
  ! 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                 *
Matthias Langer's avatar
 
Matthias Langer committed
21
22
23
24
25
26
27
28
29
  !                                                                            *
  !*****************************************************************************

  use oh_mod
  use par_mod
  use com_mod

  implicit none

30
31
32
33
34
35
36
37
38
  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
Matthias Langer's avatar
 
Matthias Langer committed
39
40
41
42
43
44
45
46
47
48
49
  real(kind=dp) :: jul

  ! Compute interval since radioactive decay of deposited mass was computed
  !************************************************************************

  if (itime.le.loutnext) then
    ldeltat=itime-(loutnext-loutstep)
  else                                  ! first half of next interval
    ldeltat=itime-loutnext
  endif

50
51
52
53
  jul=bdate+real(itime,kind=dp)/86400.
  call caldate(jul,jjjjmmdd,hhmmss)
  m=(jjjjmmdd-(jjjjmmdd/10000)*10000)/100
  h=hhmmss/10000
Matthias Langer's avatar
 
Matthias Langer committed
54

55
  ! Loop over particles
Matthias Langer's avatar
 
Matthias Langer committed
56
57
  !*****************************************

58
59
60
  do jpart=1,numpart

    ! Determine which nesting level to be used
Matthias Langer's avatar
 
Matthias Langer committed
61
62
63
64
65
66
67
68
    ngrid=0
    do j=numbnests,1,-1
      if ((xtra1(jpart).gt.xln(j)).and.(xtra1(jpart).lt.xrn(j)).and. &
           (ytra1(jpart).gt.yln(j)).and.(ytra1(jpart).lt.yrn(j))) then
        ngrid=j
        goto 23
      endif
    end do
69
23  continue
Matthias Langer's avatar
 
Matthias Langer committed
70

71
    ! Determine nested grid coordinates
Matthias Langer's avatar
 
Matthias Langer committed
72
73
74
75
76
77
78
79
80
81
    if (ngrid.gt.0) then
      xtn=(xtra1(jpart)-xln(ngrid))*xresoln(ngrid)
      ytn=(ytra1(jpart)-yln(ngrid))*yresoln(ngrid)
      ix=int(xtn)
      jy=int(ytn)
    else
      ix=int(xtra1(jpart))
      jy=int(ytra1(jpart))
    endif

82
83
84
    interp_time=nint(itime-0.5*ltsample)
    n=2
    if(abs(memtime(1)-interp_time).lt.abs(memtime(2)-interp_time)) n=1
Matthias Langer's avatar
 
Matthias Langer committed
85

86
87
88
89
90
91
    do i=2,nz
      if (height(i).gt.ztra1(jpart)) then
        indz=i-1
        goto 6
      endif
    end do
Matthias Langer's avatar
 
Matthias Langer committed
92
93
6   continue

94
95
    ! Get OH from nearest grid-cell and specific month 
    !*************************************************
Matthias Langer's avatar
 
Matthias Langer committed
96

97
    ! world coordinates
Matthias Langer's avatar
 
Matthias Langer committed
98
99
100
101
102
103
    xlon=xtra1(jpart)*dx+xlon0
    if (xlon.gt.180) then
       xlon=xlon-360
    endif
    ylat=ytra1(jpart)*dy+ylat0

104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
    ! 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))
Matthias Langer's avatar
 
Matthias Langer committed
124
125

    if (oh_average.gt.smallnum) then
126
127
128
129
130
131
132
133

      ! Computation of the OH reaction
      !**********************************************************

      temp=tt(ix,jy,indz,n)

      do k=1,nspec                                
        if (ohcconst(k).gt.0.) then
134
          ohrate=ohcconst(k)*temp**ohnconst(k)*exp(-ohdconst(k)/temp)*oh_average
135
136
137
138
139
140
141
142
          ! 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)))
Matthias Langer's avatar
 
Matthias Langer committed
143
        else
144
          ohreacted=0.
Matthias Langer's avatar
 
Matthias Langer committed
145
        endif
146
147
148
149
150
      end do

    endif  ! oh_average.gt.smallnum 

  end do  !continue loop over all particles
Matthias Langer's avatar
 
Matthias Langer committed
151
152
153


end subroutine ohreaction
154