source: trunk/src/initial_cond_calc.f90 @ 28

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