source: flexpart.git/src/re_initialize_particle.f90 @ 02095e3

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bugunivie
Last change on this file since 02095e3 was 5f9d14a, checked in by Espen Sollum ATMOS <eso@…>, 9 years ago

Updated wet depo scheme

  • Property mode set to 100644
File size: 4.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 re_initialize_particle(zp,ust,wst,h,sigmaw,wp,nrand,ol)
23!                                      i   i  i   i  i    io  io    i
24!=============== CBL skewed vertical profiles and formulation of LHH 1996 with profile of w3 from lHB 2000                                       ======
25!=============== LHH formulation has been modified to account for variable density profiles and backward in time or forward in time simulations  ======           
26!=============== this routine re-initiaalize particle velocity if a numerical instability in the cbl scheme generated a NaN value                ======
27!=============== the particle velocity is extracted from the updraft and downdraft distribution as required                                      ======
28!=============== the re-initialization si not perfect see e.g. Cassiani et al(2015) BLM                                                          ======
29!======================================================================================================================================================   
30!======================================================================================================================================================   
31  use par_mod, only:pi
32  use com_mod, only:ldirect,rannumb
33
34  implicit none
35
36
37  real :: usurad2,usurad2p,C0,costluar4,eps
38  parameter  (usurad2=0.7071067812,usurad2p=0.3989422804,C0=2,costluar4=0.66667,eps=0.000001)
39
40  integer idum,nrand
41  real :: wp,zp,ust,wst,h,dens,ddens,sigmaw,dsigmawdz,tlw,dcas,dcas1 !,ran3,gasdev
42  real :: w3,w2
43  real ::  z, &   
44       skew, &
45       skew2, &
46       radw2, &
47       fluarw,fluarw2, &
48       rluarw, &
49       xluarw, &
50       aluarw, &
51       bluarw, &
52       sigmawa, &
53       sigmawb, & 
54       ath, &
55       bth, &
56       wb,wa
57  real timedir
58  real ol,transition
59
60!--------------------------------------------------------------------------- 
61!timedir direction of time forward (1) or backward(-1)
62  nrand=nrand+1
63  dcas1=rannumb(nrand)
64  timedir=ldirect
65  z=zp/h
66  transition=1.
67
68  if (-h/ol.lt.15) transition=((sin((((-h/ol)+10.)/10.)*pi)))/2.+0.5
69
70  w2=sigmaw*sigmaw 
71  w3=(((1.2*z*((1.-z)**(3./2.)))+eps)*wst**3)*transition
72  skew=w3/(w2**1.5)
73  skew2=skew*skew
74  radw2=sqrt(w2) !sigmaw
75
76  fluarw=costluar4*skew**0.333333333333333
77  fluarw2=fluarw*fluarw
78  rluarw=(1.+fluarw2)**3.*skew2/((3.+fluarw2)**2.*fluarw2)  !-> r
79  xluarw=rluarw**0.5 !(1.+fluarw2)**1.5*skew/((3.+fluarw2)*fluarw)    !----> r^1/2
80
81  aluarw=0.5*(1.-xluarw/(4.+rluarw)**0.5)
82  bluarw=1.-aluarw
83
84  sigmawa=radw2*(bluarw/(aluarw*(1.+fluarw2)))**0.5
85  sigmawb=radw2*(aluarw/(bluarw*(1.+fluarw2)))**0.5
86
87  wa=(fluarw*sigmawa)
88  wb=(fluarw*sigmawb)
89
90
91
92  if ((sign(1.,wp)*timedir).gt.0) then !updraft   
93100 wp=(dcas1*sigmawa+wa)
94    if (wp.lt.0)  then
95      nrand=nrand+1
96      dcas1=rannumb(nrand)
97      goto 100
98    end if
99    wp=wp*timedir
100  else if ((sign(1.,wp)*timedir).lt.0) then !downdraft   
101101 wp=(dcas1*sigmawb-wb)
102    if (wp.gt.0)  then
103      nrand=nrand+1
104      dcas1=rannumb(nrand)
105      goto 101
106    end if
107    wp=wp*timedir
108  end if
109
110  return
111end subroutine re_initialize_particle
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG