Changeset 8a65cb0 in flexpart.git for src/ohreaction.f90


Ignore:
Timestamp:
Mar 2, 2015, 3:11:55 PM (9 years ago)
Author:
Espen Sollum ATMOS <espen@…>
Branches:
master, 10.4.1_pesei, GFS_025, bugfixes+enhancements, dev, release-10, release-10.4.1, scaling-bug, univie
Children:
1d207bb
Parents:
60403cd
Message:

Added code, makefile for dev branch

File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/ohreaction.f90

    re200b7a r8a65cb0  
    2525  !                                                                            *
    2626  !                                                                            *
    27   !    Author: S. Eckhardt                                                     *
     27  !    Author: R.L. Thompson                                                   *
    2828  !                                                                            *
    29   !    June 2007                                                               *
     29  !    Nov 2014                                                                *
    3030  !                                                                            *
    3131  !                                                                            *
    3232  !*****************************************************************************
    3333  ! 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                   *
     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 [molecule/cm^3] OH Concentration                                *
     41  ! ltsample [s]         interval over which mass is deposited                 *
    4242  !                                                                            *
    4343  !*****************************************************************************
    … …  
    4949  implicit none
    5050
    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
     51  integer :: jpart,itime,ltsample,loutnext,ldeltat,j,k,ix,jy!,ijx,jjy
     52  integer :: ngrid,interp_time,n,m,h,indz,i!,ia,il
     53  integer :: jjjjmmdd,hhmmss,OHx,OHy,OHz
     54  real, dimension(nzOH) :: altOHtop
     55  real :: xlon,ylat
     56  real :: xtn,ytn
     57  real :: restmass,ohreacted,oh_average
     58  real :: ohrate,temp
     59  real, parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
    6160  real(kind=dp) :: jul
    6261
    … …  
    6463  !************************************************************************
    6564
    66   gas_const=8.314 ! define gas constant
    67   act_energy=10000 ! activation energy
    68 
    69   !write(*,*) 'OH reaction n:',n,ohreact(1)
    7065  if (itime.le.loutnext) then
    7166    ldeltat=itime-(loutnext-loutstep)
    … …  
    7469  endif
    7570
     71  jul=bdate+real(itime,kind=dp)/86400.
     72  call caldate(jul,jjjjmmdd,hhmmss)
     73  m=(jjjjmmdd-(jjjjmmdd/10000)*10000)/100
     74  h=hhmmss/10000
    7675
    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
     76  ! Loop over particles
    8777  !*****************************************
    8878
     79  do jpart=1,numpart
     80
     81    ! Determine which nesting level to be used
    8982    ngrid=0
    9083    do j=numbnests,1,-1
    … …  
    9588      endif
    9689    end do
    97 23   continue
     9023  continue
    9891
    99 
    100   ! Determine nested grid coordinates
    101   !**********************************
    102 
     92    ! Determine nested grid coordinates
    10393    if (ngrid.gt.0) then
    10494      xtn=(xtra1(jpart)-xln(ngrid))*xresoln(ngrid)
    … …  
    111101    endif
    112102
    113   n=2
    114   if (abs(memtime(1)-interp_time).lt.abs(memtime(2)-interp_time)) &
    115       n=1
     103    interp_time=nint(itime-0.5*ltsample)
     104    n=2
     105    if(abs(memtime(1)-interp_time).lt.abs(memtime(2)-interp_time)) n=1
    116106
    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
     107    do i=2,nz
     108      if (height(i).gt.ztra1(jpart)) then
     109        indz=i-1
     110        goto 6
     111      endif
     112    end do
    1231136   continue
    124114
    125   ! The concentration from the nearest available gridcell is taken
    126   ! get OH concentration for the specific month and solar angle
     115    ! Get OH from nearest grid-cell and specific month
     116    !*************************************************
    127117
    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
     118    ! world coordinates
    137119    xlon=xtra1(jpart)*dx+xlon0
    138120    if (xlon.gt.180) then
    … …  
    140122    endif
    141123    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
     124
     125    ! get position in the OH field
     126    OHx=minloc(abs(lonOH-xlon),dim=1,mask=abs(lonOH-xlon).eq.minval(abs(lonOH-xlon)))
     127    OHy=minloc(abs(latOH-ylat),dim=1,mask=abs(latOH-ylat).eq.minval(abs(latOH-ylat)))
     128
     129    ! get the level of the OH field for the particle
     130    ! ztra1 is the z-coord of the trajectory above model orography in metres
     131    ! altOH is the height of the centre of the level in the OH field above orography
     132    do i=2,nzOH
     133      altOHtop(i-1)=altOH(i)+0.5*(altOH(i)-altOH(i-1))
     134    end do
     135    altOHtop(nzOH)=altOH(nzOH)+0.5*(altOH(nzOH)-altOH(nzOH-1))
     136    OHz=minloc(abs(altOHtop-ztra1(jpart)),dim=1,mask=abs(altOHtop-ztra1(jpart))&
     137            &.eq.minval(abs(altOHtop-ztra1(jpart))))
     138
     139    ! Interpolate between hourly OH fields to current time
     140    !*****************************************************
     141
     142    oh_average=OH_hourly(OHx,OHy,OHz,1)+&
     143               &(OH_hourly(OHx,OHy,OHz,2)-OH_hourly(OHx,OHy,OHz,1))*&
     144               &(itime-memOHtime(1))/(memOHtime(2)-memOHtime(1))
     145
     146    if (oh_average.gt.smallnum) then
     147
     148      ! Computation of the OH reaction
     149      !**********************************************************
     150
     151      temp=tt(ix,jy,indz,n)
     152
     153      do k=1,nspec                               
     154        if (ohcconst(k).gt.0.) then
     155          ohrate=ohcconst(k)*temp**2*exp(-ohdconst(k)/temp)*oh_average
     156          ! new particle mass
     157          restmass = xmass1(jpart,k)*exp(-1*ohrate*abs(ltsample))
     158          if (restmass .gt. smallnum) then
     159            xmass1(jpart,k)=restmass
     160          else
     161            xmass1(jpart,k)=0.
     162          endif
     163          ohreacted=xmass1(jpart,k)*(1-exp(-1*ohrate*abs(ltsample)))
     164        else
     165          ohreacted=0.
     166        endif
    155167      end do
    156 26     continue
    157168
    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
     169    endif  ! oh_average.gt.smallnum
    164170
    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   !**********************************************************
     171  end do  !continue loop over all particles
    181172
    182173
    183   !    Computation of the OH reaction
    184   !**********************************************************
    185     act_temp=tt(ix,jy,indz,n)
     174end subroutine ohreaction
    186175
    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
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG