initial_cond_calc.f90 6.05 KB
Newer Older
1
2
! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt
! SPDX-License-Identifier: GPL-3.0-or-later
3

Matthias Langer's avatar
 
Matthias Langer committed
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
subroutine initial_cond_calc(itime,i)
  !                               i   i
  !*****************************************************************************
  !                                                                            *
  !     Calculation of the sensitivity to initial conditions for BW runs       *
  !                                                                            *
  !     Author: A. Stohl                                                       *
  !                                                                            *
  !     15 January 2010                                                        *
  !                                                                            *
  !*****************************************************************************

  use unc_mod
  use outg_mod
  use par_mod
  use com_mod

  implicit none

  integer :: itime,i,ix,jy,ixp,jyp,kz,ks
  integer :: il,ind,indz,indzp,nrelpointer
  real :: rddx,rddy,p1,p2,p3,p4,dz1,dz2,dz
  real :: ddx,ddy
  real :: rhoprof(2),rhoi,xl,yl,wx,wy,w
Espen Sollum's avatar
Espen Sollum committed
28
29
  integer :: mind2
  ! mind2        eso: pointer to 2nd windfield in memory
Matthias Langer's avatar
 
Matthias Langer committed
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73


  ! For forward simulations, make a loop over the number of species;
  ! for backward simulations, make an additional loop over the release points
  !**************************************************************************


  if (itra1(i).ne.itime) return

  ! Depending on output option, calculate air density or set it to 1
  ! linit_cond: 1=mass unit, 2=mass mixing ratio unit
  !*****************************************************************


  if (linit_cond.eq.1) then     ! mass unit

    ix=int(xtra1(i))
    jy=int(ytra1(i))
    ixp=ix+1
    jyp=jy+1
    ddx=xtra1(i)-real(ix)
    ddy=ytra1(i)-real(jy)
    rddx=1.-ddx
    rddy=1.-ddy
    p1=rddx*rddy
    p2=ddx*rddy
    p3=rddx*ddy
    p4=ddx*ddy

    do il=2,nz
      if (height(il).gt.ztra1(i)) then
        indz=il-1
        indzp=il
        goto 6
      endif
    end do
6   continue

    dz1=ztra1(i)-height(indz)
    dz2=height(indzp)-ztra1(i)
    dz=1./(dz1+dz2)

  ! Take density from 2nd wind field in memory (accurate enough, no time interpolation needed)
  !*****************************************************************************
Espen Sollum's avatar
Espen Sollum committed
74
75
    mind2=memind(2)

Matthias Langer's avatar
 
Matthias Langer committed
76
    do ind=indz,indzp
Espen Sollum's avatar
Espen Sollum committed
77
78
79
80
      rhoprof(ind-indz+1)=p1*rho(ix ,jy ,ind,mind2) &
           +p2*rho(ixp,jy ,ind,mind2) &
           +p3*rho(ix ,jyp,ind,mind2) &
           +p4*rho(ixp,jyp,ind,mind2)
Matthias Langer's avatar
 
Matthias Langer committed
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
    end do
    rhoi=(dz1*rhoprof(2)+dz2*rhoprof(1))*dz
  elseif (linit_cond.eq.2) then    ! mass mixing ratio unit
    rhoi=1.
  endif


  !****************************************************************************
  ! 1. Evaluate grid concentrations using a uniform kernel of bandwidths dx, dy
  !****************************************************************************


  ! For backward simulations, look from which release point the particle comes from
  ! For domain-filling trajectory option, npoint contains a consecutive particle
  ! number, not the release point information. Therefore, nrelpointer is set to 1
  ! for the domain-filling option.
  !*****************************************************************************

  if ((ioutputforeachrelease.eq.0).or.(mdomainfill.eq.1)) then
    nrelpointer=1
  else
    nrelpointer=npoint(i)
  endif

  do kz=1,numzgrid                ! determine height of cell
    if (outheight(kz).gt.ztra1(i)) goto 21
  end do
21   continue
  if (kz.le.numzgrid) then           ! inside output domain


    xl=(xtra1(i)*dx+xoutshift)/dxout
    yl=(ytra1(i)*dy+youtshift)/dyout
    ix=int(xl)
    if (xl.lt.0.) ix=ix-1
    jy=int(yl)
    if (yl.lt.0.) jy=jy-1


  ! If a particle is close to the domain boundary, do not use the kernel either
  !****************************************************************************

    if ((xl.lt.0.5).or.(yl.lt.0.5).or. &
         (xl.gt.real(numxgrid-1)-0.5).or. &
         (yl.gt.real(numygrid-1)-0.5)) then             ! no kernel, direct attribution to grid cell
      if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
           (jy.le.numygrid-1)) then
        do ks=1,nspec
          init_cond(ix,jy,kz,ks,nrelpointer)= &
               init_cond(ix,jy,kz,ks,nrelpointer)+ &
               xmass1(i,ks)/rhoi
        end do
      endif

    else                                 ! attribution via uniform kernel

      ddx=xl-real(ix)                   ! distance to left cell border
      ddy=yl-real(jy)                   ! distance to lower cell border
      if (ddx.gt.0.5) then
        ixp=ix+1
        wx=1.5-ddx
      else
        ixp=ix-1
        wx=0.5+ddx
      endif

      if (ddy.gt.0.5) then
        jyp=jy+1
        wy=1.5-ddy
      else
        jyp=jy-1
        wy=0.5+ddy
      endif


  ! Determine mass fractions for four grid points
  !**********************************************

      if ((ix.ge.0).and.(ix.le.numxgrid-1)) then
        if ((jy.ge.0).and.(jy.le.numygrid-1)) then
          w=wx*wy
          do ks=1,nspec
            init_cond(ix,jy,kz,ks,nrelpointer)= &
                 init_cond(ix,jy,kz,ks,nrelpointer)+xmass1(i,ks)/rhoi*w
          end do
        endif

        if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then
          w=wx*(1.-wy)
          do ks=1,nspec
            init_cond(ix,jyp,kz,ks,nrelpointer)= &
                 init_cond(ix,jyp,kz,ks,nrelpointer)+xmass1(i,ks)/rhoi*w
          end do
        endif
      endif


      if ((ixp.ge.0).and.(ixp.le.numxgrid-1)) then
        if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then
          w=(1.-wx)*(1.-wy)
          do ks=1,nspec
            init_cond(ixp,jyp,kz,ks,nrelpointer)= &
                 init_cond(ixp,jyp,kz,ks,nrelpointer)+xmass1(i,ks)/rhoi*w
          end do
        endif

        if ((jy.ge.0).and.(jy.le.numygrid-1)) then
          w=(1.-wx)*wy
          do ks=1,nspec
            init_cond(ixp,jy,kz,ks,nrelpointer)= &
                 init_cond(ixp,jy,kz,ks,nrelpointer)+xmass1(i,ks)/rhoi*w
          end do
        endif
      endif
    endif

  endif

end subroutine initial_cond_calc