re_initialize_particle.f90 6.26 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
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
!***********************************************************************
!* Copyright 2012,2013                                                *
!* Jerome Brioude, Delia Arnold, Andreas Stohl, Wayne Angevine,       *
!* John Burkhart, Massimo Cassiani, Adam Dingwell, Richard C Easter, Sabine Eckhardt,*
!* Stephanie Evan, Jerome D Fast, Don Morton, Ignacio Pisso,          *
!* Petra Seibert, Gerard Wotawa, Caroline Forster, Harald Sodemann,   *
!*                                                                     *
!* This file is part of FLEXPART WRF                                   *
!*                                                                     *
!* 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 re_initialize_particle(zp,ust,wst,h,sigmaw,wp,nrand,ol)
!                                      i   i  i   i  i    i/o  i/o 
!zp: particle position
!ust: velocity scale
!wst: velocity scale
!sigmaw: vertical velcotiy standard deviation
!wp: particle velocity
!nrand: random number counter
!=============== CBL skewed vertical profiles and formulation of LHH 1996 with profile of w3 from lHB 2000   ================================================================
!=============== LHH formulation has been modified to account for variable density profiles and backward in time or forward in time simulations                             =                    
!=============== by  Massimo Cassiani ( mc ) , NILU,  2012-2013                                                                                                             =
!=============== this routine re-initialize particle velocity if a numerical instability in the cbl scheme met a specific condition                                         =
!=============== (see routine cbl.f90 and Cassiani et al. 2013                                                                                                              =
!=============== the particle velocity is extracted from the updraft and downdraft distribution as required                                                                 =
!=============== this re-initialization si not perfectly consistent with teh well-mixed condition see Cassiani et al. 2013 for details but the error introduced is small    =
!=============== but for the rpesent this is faste and simpler and shoudl be ok                                                                                             =         
!============================================================================================================================================================================   
    use par_mod, only:pi
    use com_mod, only:ldirect,rannumb
!    use ieee_arithmetic
 
    implicit none


    real :: usurad2,usurad2p,C0,costluar4,eps 
    parameter  (usurad2=0.7071067812,usurad2p=0.3989422804,C0=2,costluar4=0.66667,eps=0.000001)

    integer idum,nrand
    real :: wp,zp,ust,wst,h,dens,ddens,sigmaw,dsigmawdz,tlw,dcas,dcas1,ran3,gasdev
    real :: w3,w2,wpold
    real ::  z, &    
    skew, &
    skew2, &
    radw2, &
    fluarw,fluarw2, &
    rluarw, &
    xluarw, &
    aluarw, &
    bluarw, &
    sigmawa, &
    sigmawb, &  
    ath, &
    bth, &
    wb,wa 
    real timedir
    real ol,transition
71
    !---------------------------------------------------------------------------
72
73
74
75
76
77
78
79
80
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
!     print*,'INSIDE INIT',zp,ust,wst,h,sigmaw,wp,nrand
!   wpold=wp 
          nrand=nrand+1
     dcas1=rannumb(nrand)
     timedir=ldirect !direction of time forward (1) or backward(-1)
    z=zp/h
     
     transition=1.  !comment by mc: in this version added transtion fucntion see Cassiani et al. 2013 
     if (-h/ol.lt.15) transition=((sin((((-h/ol)+10.)/10.)*pi)))/2.+0.5
    
     !   w2=((1.7*(z*(1.-0.7*z)*(1.-z))**(2./3.))+1.e-2)*(wst**2)
    w2=sigmaw*sigmaw !this is correct and use hanna routine if commented it is for test reason 
    !w3=(((1.2*z*((1.-z)**(3./2.)))+eps)*wst**3) *1.5	!the 1.5 is to test with increased skeweness see also cbl.f90
     w3=(((1.2*z*((1.-z)**(3./2.)))+eps)*wst**3) *transition !note added a transition fucntion, comemnt by mc
    skew=w3/(w2**1.5)
    skew2=skew*skew
     radw2=sigmaw !sqrt(w2) 
    
     if (skew.ne.0) then   ! the  limit must be considered explicitly to avoid NaN
     fluarw=costluar4*skew**0.333333333333333
    fluarw2=fluarw*fluarw
    rluarw=(1.+fluarw2)**3.*skew2/((3.+fluarw2)**2.*fluarw2)  !-> r
     xluarw=rluarw**0.5 !(1.+fluarw2)**1.5*skew/((3.+fluarw2)*fluarw)    !----> r^1/2
     else        
     fluarw=0.
     fluarw2=0.
     rluarw=0.        
     xluarw=0.        
     end if     
     
     aluarw=0.5*(1.-xluarw/(4.+rluarw)**0.5)
    bluarw=1.-aluarw
        
     sigmawa=radw2*(bluarw/(aluarw*(1.+fluarw2)))**0.5
     sigmawb=radw2*(aluarw/(bluarw*(1.+fluarw2)))**0.5

    wa=(fluarw*sigmawa)
    wb=(fluarw*sigmawb)

    

    if ((sign(1.,wp)*timedir).gt.0) then !updraft
     
100     wp=(dcas1*sigmawa+wa)
      if (wp.lt.0)  then
          nrand=nrand+1
          dcas1=rannumb(nrand)
          goto 100
      end if
     wp=wp*timedir
    else if ((sign(1.,wp)*timedir).lt.0) then !downdraft
101     wp=(dcas1*sigmawb-wb)
      if (wp.gt.0)  then 
           nrand=nrand+1
          dcas1=rannumb(nrand)
          goto 101
      end if
      wp=wp*timedir
    end if   
!if (ieee_is_nan(wp)) print*,'PROBLEM INSIDE',dcas1,nrand,sigmawa,fluarw,sigmawb,wb,aluarw,fluarw,xluarw,zp,ust,wst,h,sigmaw,wpold,nrand
         return
         end