source: flexpart.git/src/initial_cond_calc.f90

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bug
Last change on this file was 92fab65, checked in by Ignacio Pisso <ip@…>, 4 years ago

add SPDX-License-Identifier to all .f90 files

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