source: branches/jerome/src_flexwrf_v3.1/initialize_cbl_vel.f90 @ 16

Last change on this file since 16 was 16, checked in by jebri, 10 years ago

sources for flexwrf v3.1

File size: 4.8 KB
Line 
1!***********************************************************************
2!* Copyright 2012,2013                                                *
3!* Jerome Brioude, Delia Arnold, Andreas Stohl, Wayne Angevine,       *
4!* John Burkhart, Massimo Cassiani, Adam Dingwell, Richard C Easter, Sabine Eckhardt,*
5!* Stephanie Evan, Jerome D Fast, Don Morton, Ignacio Pisso,          *
6!* Petra Seibert, Gerard Wotawa, Caroline Forster, Harald Sodemann,   *
7!* This file is part of FLEXPART WRF                                   *
8!*                                                                     *
9!* FLEXPART is free software: you can redistribute it and/or modify    *
10!* it under the terms of the GNU General Public License as published by*
11!* the Free Software Foundation, either version 3 of the License, or   *
12!* (at your option) any later version.                                 *
13!*                                                                     *
14!* FLEXPART is distributed in the hope that it will be useful,         *
15!* but WITHOUT ANY WARRANTY; without even the implied warranty of      *
16!* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
17!* GNU General Public License for more details.                        *
18!*                                                                     *
19!* You should have received a copy of the GNU General Public License   *
20!* along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
21!***********************************************************************
22   
23   
24    subroutine initialize_cbl_vel(idum,zp,ust,wst,h,sigmaw,wp,dcas,dcas1,ol)
25!                                  i/o   i  i   i  i     i  o   i    i       i
26! idum: for random number but not usednot used
27! zp: particle psition
28! ust: velocity scale, not used
29! wst: ocnvective velcotiy scale
30! sigmaW: standard deviaiton of vertical velocity
31! wp: particle velocity
32! dcas: for random number
33! dcas1: for random number
34! ol: Obukhov lenght
35!=============== initilization of particle velcoity based on CBL skewed vertical profiles and formulation of LHH 1996 with profile of w3 from lHB 2000         ==================
36!=============== by  Massimo Cassiani ( mc ), NILU,  2012-2013, reference to Cassiani et al. 2013..                                                            ==================
37!================================================================================================================================================================================   
38    use par_mod, only:pi
39    use com_mod, only:ldirect
40!    use ieee_arithmetic
41   
42    implicit none
43
44
45    real :: usurad2,usurad2p,C0,costluar4,eps
46    parameter  (usurad2=0.7071067812,usurad2p=0.3989422804,C0=2,costluar4=0.66667,eps=0.000001)
47
48    integer idum
49    real :: wp,zp,ust,wst,h,dens,ddens,sigmaw,dsigmawdz,tlw,dcas,dcas1,ran3,gasdev
50    real :: w3,w2
51    real ::  z, &   
52    skew, &
53    skew2, &
54    radw2, &
55    fluarw,fluarw2, &
56    rluarw, &
57    xluarw, &
58    aluarw, &
59    bluarw, &
60    sigmawa, &
61    sigmawb, & 
62    ath, &
63    bth, &
64    wb,wa
65    real timedir
66    real ol, transition
67   
68    !--------------------------------------------------------------------------- 
69         !timedir=dble(ldirect) !direction of time forward (1) or backward(-1)
70     timedir=ldirect !time direction forward (1) or backward(-1)
71         z=zp/h  !hn is the boundarylayer top
72         
73     transition=1.
74     if (-h/ol.lt.15) transition=((sin((((-h/ol)+10.)/10.)*pi)))/2.+0.5 !transtion from near neutral to unstable
75     
76     
77     !w2=(1.7*(z*(1.-0.7*z)*(1.-z))**(2./3.))*(wst**2)
78         w2=sigmaw*sigmaw
79         !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
80     w3=(((1.2*z*((1.-z)**(3./2.)))+eps)*wst**3)*transition
81         skew=w3/(w2**1.5)
82         skew2=skew*skew
83     radw2=sqrt(w2) !sigmaw
84     
85     if (skew.ne.0) then     
86     fluarw=costluar4*skew**0.333333333333333
87         fluarw2=fluarw*fluarw
88         rluarw=(1.+fluarw2)**3.*skew2/((3.+fluarw2)**2.*fluarw2)  !-> r
89     xluarw=rluarw**0.5 !(1.+fluarw2)**1.5*skew/((3.+fluarw2)*fluarw)    !----> r^1/2
90     else       
91     fluarw=0.
92     fluarw2=0.
93     rluarw=0.       
94     xluarw=0.       
95     end if   
96
97     aluarw=0.5*(1.-xluarw/(4.+rluarw)**0.5)
98         bluarw=1.-aluarw
99       
100     sigmawa=radw2*(bluarw/(aluarw*(1.+fluarw2)))**0.5
101     sigmawb=radw2*(aluarw/(bluarw*(1.+fluarw2)))**0.5
102
103         wa=(fluarw*sigmawa)
104         wb=(fluarw*sigmawb)           
105
106         !dcas=ran3(idum) !pass from outside
107
108         if (dcas.le.aluarw) then
109          !dcas1=gasdev(idum) !pass from outside
110          wp=timedir*(dcas1*sigmawa+wa)
111         else
112          !dcas1=gasdev(idum) !pass from outside
113          wp=timedir*(dcas1*sigmawb-wb)
114         end if   
115!    if (ieee_is_nan(wp)) print*,'PROBLEM INIT',wp,timedir, &
116!      dcas1,sigmawa,wa,sigmawb,wb,idum,zp,ust,wst,h,sigmaw,wp,dcas,dcas1,ol   
117    return
118    end
119         
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG