wetdepo.f90 18.1 KB
Newer Older
Matthias Langer's avatar
 
Matthias Langer committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
!                                                                     *
! This file is part of FLEXPART.                                      *
!                                                                     *
! FLEXPART is free software: you can redistribute it and/or modify    *
! it under the terms of the GNU General Public License as published by*
! the Free Software Foundation, either version 3 of the License, or   *
! (at your option) any later version.                                 *
!                                                                     *
! FLEXPART is distributed in the hope that it will be useful,         *
! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
! GNU General Public License for more details.                        *
!                                                                     *
! You should have received a copy of the GNU General Public License   *
! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
!**********************************************************************

subroutine wetdepo(itime,ltsample,loutnext)
23
  !                  i      i        i
Matthias Langer's avatar
 
Matthias Langer committed
24
25
26
27
28
29
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
  !*****************************************************************************
  !                                                                            *
  ! Calculation of wet deposition using the concept of scavenging coefficients.*
  ! For lack of detailed information, washout and rainout are jointly treated. *
  ! It is assumed that precipitation does not occur uniformly within the whole *
  ! grid cell, but that only a fraction of the grid cell experiences rainfall. *
  ! This fraction is parameterized from total cloud cover and rates of large   *
  ! scale and convective precipitation.                                        *
  !                                                                            *
  !    Author: A. Stohl                                                        *
  !                                                                            *
  !    1 December 1996                                                         *
  !                                                                            *
  ! Correction by Petra Seibert, Sept 2002:                                    *
  ! use centred precipitation data for integration                             *
  ! Code may not be correct for decay of deposition!                           *
  !                                                                            *
  !*****************************************************************************
  !                                                                            *
  ! Variables:                                                                 *
  ! cc [0-1]           total cloud cover                                       *
  ! convp [mm/h]       convective precipitation rate                           *
  ! grfraction [0-1]   fraction of grid, for which precipitation occurs        *
  ! ix,jy              indices of output grid cell for each particle           *
  ! itime [s]          actual simulation time [s]                              *
  ! jpart              particle index                                          *
  ! ldeltat [s]        interval since radioactive decay was computed           *
  ! lfr, cfr           area fraction covered by precipitation for large scale  *
  !                    and convective precipitation (dependent on prec. rate)  *
  ! loutnext [s]       time for which gridded deposition is next output        *
  ! loutstep [s]       interval at which gridded deposition is output          *
  ! lsp [mm/h]         large scale precipitation rate                          *
  ! ltsample [s]       interval over which mass is deposited                   *
  ! prec [mm/h]        precipitation rate in subgrid, where precipitation occurs*
  ! wetdeposit         mass that is wet deposited                              *
  ! wetgrid            accumulated deposited mass on output grid               *
  ! wetscav            scavenging coefficient                                  *
  !                                                                            *
  ! Constants:                                                                 *
  !                                                                            *
  !*****************************************************************************

  use point_mod
  use par_mod
  use com_mod

  implicit none

  integer :: jpart,itime,ltsample,loutnext,ldeltat,i,j,ix,jy
73
74
  integer :: ngrid,itage,nage,hz,il,interp_time, n, clouds_v
  integer :: ks, kp
75
!  integer :: n1,n2, icbot,ictop, indcloud !TEST
Matthias Langer's avatar
 
Matthias Langer committed
76
77
  real :: S_i, act_temp, cl, cle ! in cloud scavenging
  real :: clouds_h ! cloud height for the specific grid point
78
  real :: xtn,ytn,lsp,convp,cc,grfraction(3),prec(3),wetscav,totprec
Matthias Langer's avatar
 
Matthias Langer committed
79
80
  real :: wetdeposit(maxspec),restmass
  real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
Espen Sollum's avatar
Espen Sollum committed
81
  !save lfr,cfr
Matthias Langer's avatar
 
Matthias Langer committed
82

Espen Sollum's avatar
Espen Sollum committed
83
84
  real, parameter :: lfr(5) = (/ 0.5,0.65,0.8,0.9,0.95/)
  real, parameter :: cfr(5) = (/ 0.4,0.55,0.7,0.8,0.9 /)
Matthias Langer's avatar
 
Matthias Langer committed
85

86
!ZHG aerosol below-cloud scavenging removal polynomial constants for rain and snow
Espen Sollum's avatar
Espen Sollum committed
87
88
  real, parameter :: bclr(6) = (/274.35758, 332839.59273, 226656.57259, 58005.91340, 6588.38582, 0.244984/) !rain (Laakso et al 2003)
  real, parameter :: bcls(6) = (/22.7, 0.0, 0.0, 1321.0, 381.0, 0.0/) !now (Kyro et al 2009)
89
90
91
  real :: frac_act, liq_frac, dquer_m

  integer :: blc_count, inc_count
92
  real    :: Si_dummy, wetscav_dummy
93
94
95



Matthias Langer's avatar
 
Matthias Langer committed
96
97
98
99
100
101
102
103
104
105
106
107
  ! Compute interval since radioactive decay of deposited mass was computed
  !************************************************************************

  if (itime.le.loutnext) then
    ldeltat=itime-(loutnext-loutstep)
  else                                  ! first half of next interval
    ldeltat=itime-loutnext
  endif

  ! Loop over all particles
  !************************

108
109
110
  blc_count=0
  inc_count=0

Matthias Langer's avatar
 
Matthias Langer committed
111
  do jpart=1,numpart
112

Matthias Langer's avatar
 
Matthias Langer committed
113
114
115
116
117
118
    if (itra1(jpart).eq.-999999999) goto 20
    if(ldirect.eq.1)then
      if (itra1(jpart).gt.itime) goto 20
    else
      if (itra1(jpart).lt.itime) goto 20
    endif
119

Matthias Langer's avatar
 
Matthias Langer committed
120
121
122
123
124
  ! Determine age class of the particle
    itage=abs(itra1(jpart)-itramem(jpart))
    do nage=1,nageclass
      if (itage.lt.lage(nage)) goto 33
    end do
125
33  continue
Matthias Langer's avatar
 
Matthias Langer committed
126
127
128
129
130
131
132
133


  ! Determine which nesting level to be used
  !*****************************************

    ngrid=0
    do j=numbnests,1,-1
      if ((xtra1(jpart).gt.xln(j)).and.(xtra1(jpart).lt.xrn(j)).and. &
134
      (ytra1(jpart).gt.yln(j)).and.(ytra1(jpart).lt.yrn(j))) then
Matthias Langer's avatar
 
Matthias Langer committed
135
136
137
138
        ngrid=j
        goto 23
      endif
    end do
139
23  continue
Matthias Langer's avatar
 
Matthias Langer committed
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160


  ! Determine nested grid coordinates
  !**********************************

    if (ngrid.gt.0) then
      xtn=(xtra1(jpart)-xln(ngrid))*xresoln(ngrid)
      ytn=(ytra1(jpart)-yln(ngrid))*yresoln(ngrid)
      ix=int(xtn)
      jy=int(ytn)
    else
      ix=int(xtra1(jpart))
      jy=int(ytra1(jpart))
    endif


  ! Interpolate large scale precipitation, convective precipitation and
  ! total cloud cover
  ! Note that interpolated time refers to itime-0.5*ltsample [PS]
  !********************************************************************
    interp_time=nint(itime-0.5*ltsample)
Espen Sollum's avatar
Espen Sollum committed
161
    
162
163
164
165
166
167
168
169
170
    if (ngrid.eq.0) then
      call interpol_rain(lsprec,convprec,tcc,nxmax,nymax, &
      1,nx,ny,memind,real(xtra1(jpart)),real(ytra1(jpart)),1, &
      memtime(1),memtime(2),interp_time,lsp,convp,cc)
    else
      call interpol_rain_nests(lsprecn,convprecn,tccn, &
      nxmaxn,nymaxn,1,maxnests,ngrid,nxn,nyn,memind,xtn,ytn,1, &
      memtime(1),memtime(2),interp_time,lsp,convp,cc)
    endif
Matthias Langer's avatar
 
Matthias Langer committed
171

172
!  If total precipitation is less than 0.01 mm/h - no scavenging occurs
Espen Sollum's avatar
Espen Sollum committed
173
    if ((lsp.lt.0.01).and.(convp.lt.0.01)) goto 20
174

Matthias Langer's avatar
 
Matthias Langer committed
175
  ! get the level were the actual particle is in
176
177
178
179
180
181
182
    do il=2,nz
      if (height(il).gt.ztra1(jpart)) then
        hz=il-1
        goto 26
      endif
    end do
26  continue
Matthias Langer's avatar
 
Matthias Langer committed
183

184
185
    n=memind(2)
    if (abs(memtime(1)-interp_time).lt.abs(memtime(2)-interp_time)) &
186
         n=memind(1)
Matthias Langer's avatar
 
Matthias Langer committed
187

188
189
190
191
    if (ngrid.eq.0) then
      clouds_v=clouds(ix,jy,hz,n)
      clouds_h=cloudsh(ix,jy,n)
    else
192
      ! new removal not implemented for nests yet 
193
194
195
      clouds_v=cloudsn(ix,jy,hz,n,ngrid)
      clouds_h=cloudsnh(ix,jy,n,ngrid)
    endif
Matthias Langer's avatar
 
Matthias Langer committed
196

197
198
199
  ! if there is no precipitation or the particle is above the clouds no
  ! scavenging is done

200
    if (clouds_v.le.1) goto 20
201
  
Matthias Langer's avatar
 
Matthias Langer committed
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
  ! 1) Parameterization of the the area fraction of the grid cell where the
  !    precipitation occurs: the absolute limit is the total cloud cover, but
  !    for low precipitation rates, an even smaller fraction of the grid cell
  !    is used. Large scale precipitation occurs over larger areas than
  !    convective precipitation.
  !**************************************************************************

    if (lsp.gt.20.) then
      i=5
    else if (lsp.gt.8.) then
      i=4
    else if (lsp.gt.3.) then
      i=3
    else if (lsp.gt.1.) then
      i=2
    else
      i=1
    endif

    if (convp.gt.20.) then
      j=5
    else if (convp.gt.8.) then
      j=4
    else if (convp.gt.3.) then
      j=3
    else if (convp.gt.1.) then
      j=2
    else
      j=1
    endif

233

234
235
236
  !ZHG oct 2014 : Calculated for 1) both 2) lsp 3) convp 
  ! Tentatively differentiate the grfraction for lsp and convp for treating differently the two forms
  ! for now they are treated the same
237
238
239
240
    grfraction(1)=max(0.05,cc*(lsp*lfr(i)+convp*cfr(j))/(lsp+convp))
    grfraction(2)=max(0.05,cc*(lfr(i)))
    grfraction(3)=max(0.05,cc*(cfr(j)))

Matthias Langer's avatar
 
Matthias Langer committed
241
242
243

  ! 2) Computation of precipitation rate in sub-grid cell
  !******************************************************
244
245
246
    prec(1)=(lsp+convp)/grfraction(1)
    prec(2)=(lsp)/grfraction(2)
    prec(3)=(convp)/grfraction(3)
Matthias Langer's avatar
 
Matthias Langer committed
247
248
249
250
251
252


  ! 3) Computation of scavenging coefficients for all species
  !    Computation of wet deposition
  !**********************************************************

253
254
255
256
    do ks=1,nspec      ! loop over species
      wetdeposit(ks)=0. 
      wetscav=0.   

257
      !ZHG test if it nested?
Espen Sollum's avatar
Espen Sollum committed
258
259
260
261
262
      if (ngrid.gt.0) then
        act_temp=ttn(ix,jy,hz,n,ngrid)
      else
        act_temp=tt(ix,jy,hz,n)
      endif
263
     
264

265
  !***********************
266
267
  ! BELOW CLOUD SCAVENGING
  !***********************  
268
269
      if (clouds_v.ge.4) then !below cloud

flexpart's avatar
flexpart committed
270
        if (weta(ks).gt.0. .or. wetb(ks).gt.0.) then !if positive below-cloud parameters given in SPECIES file (either A or B)
Espen Sollum's avatar
Espen Sollum committed
271
          blc_count=blc_count+1
272
273

!GAS
Espen Sollum's avatar
Espen Sollum committed
274
          if (dquer(ks) .le. 0.) then  !is gas
275
276
277
            wetscav=weta(ks)*prec(1)**wetb(ks)
            
!AEROSOL
Espen Sollum's avatar
Espen Sollum committed
278
          else !is particle
279
280
281
282
!NIK 17.02.2015
! For the calculation here particle size needs to be in meter and not um as dquer is changed to in readreleases
! for particles larger than 10 um use the largest size defined in the parameterizations (10um)
            dquer_m=min(10.,dquer(ks))/1000000. !conversion from um to m
flexpart's avatar
flexpart committed
283
            if (act_temp .ge. 273 .and. weta(ks).gt.0.)  then !Rain 
284
285
286
              ! ZHG 2014 : Particle RAIN scavenging coefficient based on Laakso et al 2003, 
              ! the below-cloud scavenging (rain efficienty) 
              ! parameter A (=weta) from SPECIES file
Espen Sollum's avatar
Espen Sollum committed
287
288
              wetscav= weta(ks)*10**(bclr(1)+ (bclr(2)*(log10(dquer_m))**(-4))+(bclr(3)*(log10(dquer_m))**(-3))+ (bclr(4)* &
                   (log10(dquer_m))**(-2))+ (bclr(5)*(log10(dquer_m))**(-1))+ bclr(6)* (prec(1))**(0.5))
289

290
291
292
293
            elseif (act_temp .lt. 273 .and. wetb(ks).gt.0.)  then ! Snow
              ! ZHG 2014 : Particle SNOW scavenging coefficient based on Kyro et al 2009, 
              ! the below-cloud scavenging (Snow efficiency) 
              ! parameter B (=wetb) from SPECIES file
Espen Sollum's avatar
Espen Sollum committed
294
295
              wetscav= wetb(ks)*10**(bcls(1)+ (bcls(2)*(log10(dquer_m))**(-4))+(bcls(3)*(log10(dquer_m))**(-3))+ (bcls(4)* &
                   (log10(dquer_m))**(-2))+ (bcls(5)*(log10(dquer_m))**(-1))+ bcls(6)* (prec(1))**(0.5))
296

Espen Sollum's avatar
Espen Sollum committed
297
            endif
298
299
300

!             write(*,*) 'bl-cloud, act_temp=',act_temp, ',prec=',prec(1),',wetscav=', wetscav, ', jpart=',jpart

Espen Sollum's avatar
Espen Sollum committed
301
302
          endif !gas or particle
        endif ! positive below-cloud scavenging parameters given in Species file
303
      endif !end BELOW
304
305
306
307


  !********************
  ! IN CLOUD SCAVENGING
308
309
      !******************************************************
      if (clouds_v.lt.4) then ! In-cloud
flexpart's avatar
flexpart committed
310
! NIK 13 may 2015: only do incloud if positive in-cloud scavenging parameters are given in species file
311
        if (weta_in(ks).gt.0. .or. wetb_in(ks).gt.0.) then 
312
          inc_count=inc_count+1
flexpart's avatar
flexpart committed
313
314
315
! if negative coefficients (turned off) set to zero for use in equation
          if (weta_in(ks).lt.0.) weta_in(ks)=0.
          if (wetb_in(ks).lt.0.) wetb_in(ks)=0.
Espen Sollum's avatar
Espen Sollum committed
316

317
318
          !ZHG 2015 Cloud liquid & ice water (CLWC+CIWC) from ECMWF
          if (readclouds) then                  !icloud_stats(ix,jy,4,n) has units kg/m2
319
320
321
!            cl =icloud_stats(ix,jy,4,n)*(grfraction(1)/cc)
! ESO: stop using icloud_stats
            cl =clw4(ix,jy,n)*(grfraction(1)/cc)
322
323
324
          else                                  !parameterize cloudwater m2/m3
            !ZHG updated parameterization of cloud water to better reproduce the values coming from ECMWF
            cl=1.6E-6*prec(1)**0.36
325
326
          endif

327
            !ZHG: Calculate the partition between liquid and water phase water. 
flexpart's avatar
flexpart committed
328
329
330
331
332
333
334
            if (act_temp .le. 253) then
              liq_frac=0
            else if (act_temp .ge. 273) then
              liq_frac=1
            else
              liq_frac =((act_temp-273)/(273-253))**2
            endif
335
           ! ZHG: Calculate the aerosol partition based on cloud phase and Ai and Bi
flexpart's avatar
flexpart committed
336
            frac_act = liq_frac*weta_in(ks) +(1-liq_frac)*wetb_in(ks)
337
 
Espen Sollum's avatar
Espen Sollum committed
338
  !ZHG Use the activated fraction and the liqid water to calculate the washout
flexpart's avatar
flexpart committed
339
340

  ! AEROSOL
341
          !**************************************************
flexpart's avatar
flexpart committed
342
343
344
          if (dquer(ks).gt. 0.) then ! is particle

            S_i= frac_act/cl
345

346
          !*********************
Espen Sollum's avatar
Espen Sollum committed
347
  ! GAS
flexpart's avatar
flexpart committed
348
349
350
          else ! is gas
               
            cle=(1-cl)/(henry(ks)*(r_air/3500.)*act_temp)+cl
351
352
353
            !REPLACE to switch old/ new scheme 
            ! S_i=frac_act/cle
            S_i=1/cle
flexpart's avatar
flexpart committed
354
          endif ! gas or particle
355

Espen Sollum's avatar
Espen Sollum committed
356
  ! scavenging coefficient based on Hertel et al 1995 - using the S_i for either gas or aerosol
357
358
           !OLD 
          if (readclouds) then
359
            wetscav=S_i*(prec(1)/3.6E6)
360
          else
361
            wetscav=S_i*(prec(1)/3.6E6)/clouds_h
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
          endif

!ZHG 2015 TEST
!          Si_dummy=frac_act/2E-7*prec(1)**0.36
!           wetscav_dummy=Si_dummy*(prec(1)/3.6E6)/clouds_h
!           if (clouds_v.lt.4) then
!           talltest=talltest+1
!if(talltest .eq. 1) OPEN(UNIT=199, FILE=utfil,FORM='FORMATTED',STATUS = 'UNKNOWN')
!if(talltest .lt. 100001)  write(199,*) prec(1)/3.6E6, cl, clouds_h*2E-7*prec(1)**0.36,clouds_v,ytra1(jpart)-90
!if(talltest .lt. 100001)  write(199,*) wetscav, wetscav_dummy,prec(1),ytra1(jpart)-90,clouds_v,cl
!if(talltest .eq. 100001) CLOSE(199)
!if(talltest .eq. 100001) STOP
!
!write(*,*)  'PREC kg/m2s CLOUD kg/m2', (prec(1)/3.6E6), cl !, '2E-7*prec(1)**0.36',  2E-7*prec(1)**0.36,'2E-7*prec(1)**0.36*clouds_h',2E-7*prec(1)**0.36*clouds_h
!write(*,*)  'PREC kg/m2s LSP+convp kg/m2', prec(1), convp+lsp
!write(*,*)  wetscav, wetscav_dummy
!write(*,*) cc, grfraction(1), cc/grfraction(1)

!write(*,*)  'Lmbda_old', (prec(1)/3.6E6)/(clouds_h*2E-7*prec(1)**0.36)


!write(*,*) '**************************************************'
!write(*,*)  'clouds_h', clouds_h, 'clouds_v',clouds_v,'abs(ltsample)', abs(ltsample)
!write(*,*)  'readclouds', readclouds, 'wetscav',wetscav, 'wetscav_dummy', wetscav_dummy
!write(*,*)  'S_i', S_i , 'Si_dummy', Si_dummy, 'prec(1)', prec(1) 


!           write(*,*) 'PRECIPITATION ,cl  ECMWF , cl PARAMETIZED, clouds_v, lat' &
!                      ,prec(1)/3.6E6, cl, clouds_h*2E-7*prec(1)**0.36,clouds_v,ytra1(jpart)-90

!endif 
393

flexpart's avatar
flexpart committed
394
        endif ! positive in-cloud scavenging parameters given in Species file
395
      endif !incloud
396
!END ZHG TEST
Espen Sollum's avatar
Espen Sollum committed
397
     
398
399
400
401
402
  !**************************************************
  ! CALCULATE DEPOSITION 
  !**************************************************
  !     if (wetscav.le.0) write (*,*) 'neg, or 0 wetscav!'
  !     +          ,wetscav,cle,cl,act_temp,prec,clouds_h,clouds_v
Matthias Langer's avatar
 
Matthias Langer committed
403

404
      if (wetscav.gt.0.) then
Matthias Langer's avatar
 
Matthias Langer committed
405
        wetdeposit(ks)=xmass1(jpart,ks)* &
Espen Sollum's avatar
Espen Sollum committed
406
             (1.-exp(-wetscav*abs(ltsample)))*grfraction(1)  ! wet deposition
407
408
409
410
!write(*,*) 'MASS DEPOSITED: PREC, WETSCAV, WETSCAVP', prec(1), wetdeposit(ks), xmass1(jpart,ks)* & 
!             (1.-exp(-wetscav_dummy*abs(ltsample)))*grfraction(1), clouds_v


411
412
413
      else ! if no scavenging
        wetdeposit(ks)=0.
      endif
414

415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
      restmass = xmass1(jpart,ks)-wetdeposit(ks)
      if (ioutputforeachrelease.eq.1) then
        kp=npoint(jpart)
      else
        kp=1
      endif
      if (restmass .gt. smallnum) then
        xmass1(jpart,ks)=restmass
  !   depostatistic
  !   wetdepo_sum(ks,kp)=wetdepo_sum(ks,kp)+wetdeposit(ks)
  !   depostatistic
      else
        xmass1(jpart,ks)=0.
      endif
  !   Correct deposited mass to the last time step when radioactive decay of
  !   gridded deposited mass was calculated
      if (decay(ks).gt.0.) then
        wetdeposit(ks)=wetdeposit(ks)*exp(abs(ldeltat)*decay(ks))
      endif
434
435


436
    end do !all species
Matthias Langer's avatar
 
Matthias Langer committed
437

438
  ! Sabine Eckhardt, June 2008 create deposition runs only for forward runs
Matthias Langer's avatar
 
Matthias Langer committed
439
440
441
  ! Add the wet deposition to accumulated amount on output grid and nested output grid
  !*****************************************************************************

442
443
    if (ldirect.eq.1) then
      call wetdepokernel(nclass(jpart),wetdeposit,real(xtra1(jpart)), &
Espen Sollum's avatar
Espen Sollum committed
444
           real(ytra1(jpart)),nage,kp)
445
      if (nested_output.eq.1) call wetdepokernel_nest(nclass(jpart), &
Espen Sollum's avatar
Espen Sollum committed
446
           wetdeposit,real(xtra1(jpart)),real(ytra1(jpart)),nage,kp)
447
    endif
Matthias Langer's avatar
 
Matthias Langer committed
448
449

20  continue
450
  end do ! all particles
Matthias Langer's avatar
 
Matthias Langer committed
451

Espen Sollum's avatar
Espen Sollum committed
452
453
454
  ! count the total number of below-cloud and in-cloud occurences:
  tot_blc_count=tot_blc_count+blc_count
  tot_inc_count=tot_inc_count+inc_count
455

Matthias Langer's avatar
 
Matthias Langer committed
456
end subroutine wetdepo