source: flexpart.git/src/initial_cond_calc.f90 @ 02095e3

devrelease-10univie
Last change on this file since 02095e3 was 5f9d14a, checked in by Espen Sollum ATMOS <eso@…>, 5 years ago

Updated wet depo scheme

  • Property mode set to 100644
File size: 7.3 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 initial_cond_calc(itime,i)
23  !                               i   i
24  !*****************************************************************************
25  !                                                                            *
26  !     Calculation of the sensitivity to initial conditions for BW runs       *
27  !                                                                            *
28  !     Author: A. Stohl                                                       *
29  !                                                                            *
30  !     15 January 2010                                                        *
31  !                                                                            *
32  !*****************************************************************************
33
34  use unc_mod
35  use outg_mod
36  use par_mod
37  use com_mod
38
39  implicit none
40
41  integer :: itime,i,ix,jy,ixp,jyp,kz,ks
42  integer :: il,ind,indz,indzp,nrelpointer
43  real :: rddx,rddy,p1,p2,p3,p4,dz1,dz2,dz
44  real :: ddx,ddy
45  real :: rhoprof(2),rhoi,xl,yl,wx,wy,w
46  integer :: mind2
47  ! mind2        eso: pointer to 2nd windfield in memory
48
49
50  ! For forward simulations, make a loop over the number of species;
51  ! for backward simulations, make an additional loop over the release points
52  !**************************************************************************
53
54
55  if (itra1(i).ne.itime) return
56
57  ! Depending on output option, calculate air density or set it to 1
58  ! linit_cond: 1=mass unit, 2=mass mixing ratio unit
59  !*****************************************************************
60
61
62  if (linit_cond.eq.1) then     ! mass unit
63
64    ix=int(xtra1(i))
65    jy=int(ytra1(i))
66    ixp=ix+1
67    jyp=jy+1
68    ddx=xtra1(i)-real(ix)
69    ddy=ytra1(i)-real(jy)
70    rddx=1.-ddx
71    rddy=1.-ddy
72    p1=rddx*rddy
73    p2=ddx*rddy
74    p3=rddx*ddy
75    p4=ddx*ddy
76
77    do il=2,nz
78      if (height(il).gt.ztra1(i)) then
79        indz=il-1
80        indzp=il
81        goto 6
82      endif
83    end do
846   continue
85
86    dz1=ztra1(i)-height(indz)
87    dz2=height(indzp)-ztra1(i)
88    dz=1./(dz1+dz2)
89
90  ! Take density from 2nd wind field in memory (accurate enough, no time interpolation needed)
91  !*****************************************************************************
92    mind2=memind(2)
93
94    do ind=indz,indzp
95      rhoprof(ind-indz+1)=p1*rho(ix ,jy ,ind,mind2) &
96           +p2*rho(ixp,jy ,ind,mind2) &
97           +p3*rho(ix ,jyp,ind,mind2) &
98           +p4*rho(ixp,jyp,ind,mind2)
99    end do
100    rhoi=(dz1*rhoprof(2)+dz2*rhoprof(1))*dz
101  elseif (linit_cond.eq.2) then    ! mass mixing ratio unit
102    rhoi=1.
103  endif
104
105
106  !****************************************************************************
107  ! 1. Evaluate grid concentrations using a uniform kernel of bandwidths dx, dy
108  !****************************************************************************
109
110
111  ! For backward simulations, look from which release point the particle comes from
112  ! For domain-filling trajectory option, npoint contains a consecutive particle
113  ! number, not the release point information. Therefore, nrelpointer is set to 1
114  ! for the domain-filling option.
115  !*****************************************************************************
116
117  if ((ioutputforeachrelease.eq.0).or.(mdomainfill.eq.1)) then
118    nrelpointer=1
119  else
120    nrelpointer=npoint(i)
121  endif
122
123  do kz=1,numzgrid                ! determine height of cell
124    if (outheight(kz).gt.ztra1(i)) goto 21
125  end do
12621   continue
127  if (kz.le.numzgrid) then           ! inside output domain
128
129
130    xl=(xtra1(i)*dx+xoutshift)/dxout
131    yl=(ytra1(i)*dy+youtshift)/dyout
132    ix=int(xl)
133    if (xl.lt.0.) ix=ix-1
134    jy=int(yl)
135    if (yl.lt.0.) jy=jy-1
136
137
138  ! If a particle is close to the domain boundary, do not use the kernel either
139  !****************************************************************************
140
141    if ((xl.lt.0.5).or.(yl.lt.0.5).or. &
142         (xl.gt.real(numxgrid-1)-0.5).or. &
143         (yl.gt.real(numygrid-1)-0.5)) then             ! no kernel, direct attribution to grid cell
144      if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
145           (jy.le.numygrid-1)) then
146        do ks=1,nspec
147          init_cond(ix,jy,kz,ks,nrelpointer)= &
148               init_cond(ix,jy,kz,ks,nrelpointer)+ &
149               xmass1(i,ks)/rhoi
150        end do
151      endif
152
153    else                                 ! attribution via uniform kernel
154
155      ddx=xl-real(ix)                   ! distance to left cell border
156      ddy=yl-real(jy)                   ! distance to lower cell border
157      if (ddx.gt.0.5) then
158        ixp=ix+1
159        wx=1.5-ddx
160      else
161        ixp=ix-1
162        wx=0.5+ddx
163      endif
164
165      if (ddy.gt.0.5) then
166        jyp=jy+1
167        wy=1.5-ddy
168      else
169        jyp=jy-1
170        wy=0.5+ddy
171      endif
172
173
174  ! Determine mass fractions for four grid points
175  !**********************************************
176
177      if ((ix.ge.0).and.(ix.le.numxgrid-1)) then
178        if ((jy.ge.0).and.(jy.le.numygrid-1)) then
179          w=wx*wy
180          do ks=1,nspec
181            init_cond(ix,jy,kz,ks,nrelpointer)= &
182                 init_cond(ix,jy,kz,ks,nrelpointer)+xmass1(i,ks)/rhoi*w
183          end do
184        endif
185
186        if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then
187          w=wx*(1.-wy)
188          do ks=1,nspec
189            init_cond(ix,jyp,kz,ks,nrelpointer)= &
190                 init_cond(ix,jyp,kz,ks,nrelpointer)+xmass1(i,ks)/rhoi*w
191          end do
192        endif
193      endif
194
195
196      if ((ixp.ge.0).and.(ixp.le.numxgrid-1)) then
197        if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then
198          w=(1.-wx)*(1.-wy)
199          do ks=1,nspec
200            init_cond(ixp,jyp,kz,ks,nrelpointer)= &
201                 init_cond(ixp,jyp,kz,ks,nrelpointer)+xmass1(i,ks)/rhoi*w
202          end do
203        endif
204
205        if ((jy.ge.0).and.(jy.le.numygrid-1)) then
206          w=(1.-wx)*wy
207          do ks=1,nspec
208            init_cond(ixp,jy,kz,ks,nrelpointer)= &
209                 init_cond(ixp,jy,kz,ks,nrelpointer)+xmass1(i,ks)/rhoi*w
210          end do
211        endif
212      endif
213    endif
214
215  endif
216
217end subroutine initial_cond_calc
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG