source: flexpart.git/src/get_vdep_prob.f90 @ 7062fab

devrelease-10univie
Last change on this file since 7062fab was 7062fab, checked in by Sabine <sabine.eckhardt@…>, 3 years ago

bug fix zpoint in timem-get_vdep_prob/advance_rec

  • Property mode set to 100644
File size: 8.0 KB
Line 
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
22subroutine advance_rec(itime,xt,yt,zt,prob)
23  !                    i     i  i  i  o
24  !*****************************************************************************
25  !                                                                            *
26  !  Calculation of turbulent particle trajectories utilizing a                *
27  !  zero-acceleration scheme, which is corrected by a numerically more        *
28  !  accurate Petterssen scheme whenever possible.                             *
29  !                                                                            *
30  !  Particle positions are read in, incremented, and returned to the calling  *
31  !  program.                                                                  *
32  !                                                                            *
33  !  In different regions of the atmosphere (PBL vs. free troposphere),        *
34  !  different parameters are needed for advection, parameterizing turbulent   *
35  !  velocities, etc. For efficiency, different interpolation routines have    *
36  !  been written for these different cases, with the disadvantage that there  *
37  !  exist several routines doing almost the same. They all share the          *
38  !  included file 'interpol_mod'. The following                               *
39  !  interpolation routines are used:                                          *
40  !                                                                            *
41  !  interpol_all(_nests)     interpolates everything (called inside the PBL)  *
42  !  interpol_misslev(_nests) if a particle moves vertically in the PBL,       *
43  !                           additional parameters are interpolated if it     *
44  !                           crosses a model level                            *
45  !  interpol_wind(_nests)    interpolates the wind and determines the         *
46  !                           standard deviation of the wind (called outside   *
47  !                           PBL) also interpolates potential vorticity       *
48  !  interpol_wind_short(_nests) only interpolates the wind (needed for the    *
49  !                           Petterssen scheme)                               *
50  !  interpol_vdep(_nests)    interpolates deposition velocities               *
51  !                                                                            *
52  !                                                                            *
53  !     Author: A. Stohl                                                       *
54  !                                                                            *
55  !     16 December 1997                                                       *
56  !                                                                            *
57  !  Changes:                                                                  *
58  !                                                                            *
59  !  8 April 2000: Deep convection parameterization                            *
60  !                                                                            *
61  !  May 2002: Petterssen scheme introduced                                    *
62  !                                                                            *
63  !*****************************************************************************
64  !                                                                            *
65  ! Variables:                                                                 *
66  ! itime [s]          time at which this subroutine is entered                *
67  ! itimec [s]         actual time, which is incremented in this subroutine    *
68  ! href [m]           height for which dry deposition velocity is calculated  *
69  ! ldirect            1 forward, -1 backward                                  *
70  ! ldt [s]            Time step for the next integration                      *
71  ! lsynctime [s]      Synchronisation interval of FLEXPART                    *
72  ! ngrid              index which grid is to be used                          *
73  ! prob               probability of absorption due to dry deposition         *
74  ! vdepo              Deposition velocities for all species                   *
75  ! xt,yt,zt           Particle position                                       *
76  !                                                                            *
77  !*****************************************************************************
78
79  use point_mod
80  use par_mod
81  use com_mod
82  use interpol_mod
83
84  implicit none
85
86  real(kind=dp) :: xt,yt
87  real :: zt,xtn,ytn
88  integer :: itime,i,j,k,memindnext
89  integer :: nix,njy,ks
90  real :: prob(maxspec),vdepo(maxspec)
91  real,parameter :: eps=nxmax/3.e5
92
93  if (DRYDEP) then    ! reset probability for deposition
94    do ks=1,nspec
95      depoindicator(ks)=.true.
96      prob(ks)=0.
97    end do
98  endif
99
100
101  ! Determine whether lat/long grid or polarstereographic projection
102  ! is to be used
103  ! Furthermore, determine which nesting level to be used
104  !*****************************************************************
105
106  if (nglobal.and.(yt.gt.switchnorthg)) then
107    ngrid=-1
108  else if (sglobal.and.(yt.lt.switchsouthg)) then
109    ngrid=-2
110  else
111    ngrid=0
112    do j=numbnests,1,-1
113      if ((xt.gt.xln(j)+eps).and.(xt.lt.xrn(j)-eps).and. &
114           (yt.gt.yln(j)+eps).and.(yt.lt.yrn(j)-eps)) then
115        ngrid=j
116        goto 23
117      endif
118    end do
11923   continue
120  endif
121
122
123  !***************************
124  ! Interpolate necessary data
125  !***************************
126
127  if (abs(itime-memtime(1)).lt.abs(itime-memtime(2))) then
128    memindnext=1
129  else
130    memindnext=2
131  endif
132
133  ! Determine nested grid coordinates
134  !**********************************
135
136  if (ngrid.gt.0) then
137    xtn=(xt-xln(ngrid))*xresoln(ngrid)
138    ytn=(yt-yln(ngrid))*yresoln(ngrid)
139    ix=int(xtn)
140    jy=int(ytn)
141    nix=nint(xtn)
142    njy=nint(ytn)
143  else
144    ix=int(xt)
145    jy=int(yt)
146    nix=nint(xt)
147    njy=nint(yt)
148  endif
149  ixp=ix+1
150  jyp=jy+1
151
152
153  ! Determine probability of deposition
154  !************************************
155
156      if ((DRYDEP).and.(zt.lt.2.*href)) then
157        do ks=1,nspec
158          if (DRYDEPSPEC(ks)) then
159            if (depoindicator(ks)) then
160              if (ngrid.le.0) then
161                call interpol_vdep(ks,vdepo(ks))
162              else
163                call interpol_vdep_nests(ks,vdepo(ks))
164              endif
165            endif
166  ! correction by Petra Seibert, 10 April 2001
167  !   this formulation means that prob(n) = 1 - f(0)*...*f(n)
168  !   where f(n) is the exponential term
169               prob(ks)=vdepo(ks)
170!               prob(ks)=vdepo(ks)/2./href
171! instead of prob - return vdepo -> result kg/m2/s
172          endif
173        end do
174      endif
175
176
177end subroutine advance_rec
178
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG