source: branches/flexpart91_hasod/src_parallel/initial_cond_calc.f90 @ 8

Last change on this file since 8 was 8, checked in by hasod, 9 years ago

Added parallel version of Flexpart91

File size: 7.5 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, only: init_cond
35  use outg_mod, only: outheight
36  use com_mod, only: linit_cond, dx, dy, xmass1, ztra1, xtra1, ytra1, npoint, &
37      height, rho, itra1, xoutshift, youtshift, nz, numzgrid, numygrid, numxgrid, &
38      mdomainfill, nspec, dxout, dyout, ioutputforeachrelease
39
40  implicit none
41
42  integer :: itime,i,ix,jy,ixp,jyp,kz,ks
43  integer :: il,ind,indz,indzp,nrelpointer
44  real :: rddx,rddy,p1,p2,p3,p4,dz1,dz2,dz
45  real :: ddx,ddy
46  real :: rhoprof(2),rhoi,xl,yl,wx,wy,w
47
48
49  ! For forward simulations, make a loop over the number of species;
50  ! for backward simulations, make an additional loop over the release points
51  !**************************************************************************
52
53
54  if (itra1(i).ne.itime) return
55
56  ! Depending on output option, calculate air density or set it to 1
57  ! linit_cond: 1=mass unit, 2=mass mixing ratio unit
58  !*****************************************************************
59
60
61  if (linit_cond.eq.1) then     ! mass unit
62
63    ix=int(xtra1(i))
64    jy=int(ytra1(i))
65    ixp=ix+1
66    jyp=jy+1
67    ddx=xtra1(i)-real(ix)
68    ddy=ytra1(i)-real(jy)
69    rddx=1.-ddx
70    rddy=1.-ddy
71    p1=rddx*rddy
72    p2=ddx*rddy
73    p3=rddx*ddy
74    p4=ddx*ddy
75
76    do il=2,nz
77      if (height(il).gt.ztra1(i)) then
78        indz=il-1
79        indzp=il
80        goto 6
81      endif
82    end do
836   continue
84
85    dz1=ztra1(i)-height(indz)
86    dz2=height(indzp)-ztra1(i)
87    dz=1./(dz1+dz2)
88
89  ! Take density from 2nd wind field in memory (accurate enough, no time interpolation needed)
90  !*****************************************************************************
91    do ind=indz,indzp
92      rhoprof(ind-indz+1)=p1*rho(ix ,jy ,ind,2) &
93           +p2*rho(ixp,jy ,ind,2) &
94           +p3*rho(ix ,jyp,ind,2) &
95           +p4*rho(ixp,jyp,ind,2)
96    end do
97    rhoi=(dz1*rhoprof(2)+dz2*rhoprof(1))*dz
98  elseif (linit_cond.eq.2) then    ! mass mixing ratio unit
99    rhoi=1.
100  endif
101
102
103  !****************************************************************************
104  ! 1. Evaluate grid concentrations using a uniform kernel of bandwidths dx, dy
105  !****************************************************************************
106
107
108  ! For backward simulations, look from which release point the particle comes from
109  ! For domain-filling trajectory option, npoint contains a consecutive particle
110  ! number, not the release point information. Therefore, nrelpointer is set to 1
111  ! for the domain-filling option.
112  !*****************************************************************************
113
114  if ((ioutputforeachrelease.eq.0).or.(mdomainfill.eq.1)) then
115    nrelpointer=1
116  else
117    nrelpointer=npoint(i)
118  endif
119
120  do kz=1,numzgrid                ! determine height of cell
121    if (outheight(kz).gt.ztra1(i)) goto 21
122  end do
12321   continue
124  if (kz.le.numzgrid) then           ! inside output domain
125
126
127    xl=(xtra1(i)*dx+xoutshift)/dxout
128    yl=(ytra1(i)*dy+youtshift)/dyout
129    ix=int(xl)
130    if (xl.lt.0.) ix=ix-1
131    jy=int(yl)
132    if (yl.lt.0.) jy=jy-1
133
134
135  ! If a particle is close to the domain boundary, do not use the kernel either
136  !****************************************************************************
137
138    if ((xl.lt.0.5).or.(yl.lt.0.5).or. &
139         (xl.gt.real(numxgrid-1)-0.5).or. &
140         (yl.gt.real(numygrid-1)-0.5)) then             ! no kernel, direct attribution to grid cell
141      if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
142           (jy.le.numygrid-1)) then
143!$OMP CRITICAL
144        do ks=1,nspec
145          init_cond(ix,jy,kz,ks,nrelpointer)= &
146               init_cond(ix,jy,kz,ks,nrelpointer)+ &
147               xmass1(i,ks)/rhoi
148        end do
149!$OMP END CRITICAL
150      endif
151
152    else                                 ! attribution via uniform kernel
153
154      ddx=xl-real(ix)                   ! distance to left cell border
155      ddy=yl-real(jy)                   ! distance to lower cell border
156      if (ddx.gt.0.5) then
157        ixp=ix+1
158        wx=1.5-ddx
159      else
160        ixp=ix-1
161        wx=0.5+ddx
162      endif
163
164      if (ddy.gt.0.5) then
165        jyp=jy+1
166        wy=1.5-ddy
167      else
168        jyp=jy-1
169        wy=0.5+ddy
170      endif
171
172
173  ! Determine mass fractions for four grid points
174  !**********************************************
175
176!$OMP CRITICAL
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!$OMP END CRITICAL
214
215    endif
216  endif
217
218end subroutine initial_cond_calc
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG