source: flexpart.git/src/initialize_cbl_vel.f90 @ 3481cc1

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bug
Last change on this file since 3481cc1 was 3481cc1, checked in by Ignacio Pisso <ip@…>, 4 years ago

move license from headers to a different file

  • Property mode set to 100644
File size: 2.2 KB
Line 
1subroutine initialize_cbl_vel(idum,zp,ust,wst,h,sigmaw,wp, ol)
2!                              i/o   i  i   i  i     i  o   i 
3
4  use par_mod, only:pi
5  use com_mod, only:ldirect
6  use random_mod, only: gasdev, ran3
7
8  implicit none
9
10!===============================================================================
11! CBL skewed vertical profiles and formulation of LHH 1996 with profile of w3
12! from LHB 2000
13! LHH formulation has been modified to account for variable density profiles and
14! backward in time or forward in time simulations
15! see Cassiani et al. BLM 2014 doi  for explanations and references
16!===============================================================================
17
18
19  real :: usurad2,usurad2p,C0,costluar4,eps
20  parameter  (usurad2=0.7071067812,usurad2p=0.3989422804,C0=2,costluar4=0.66667,eps=0.000001)
21
22  integer idum
23  real :: wp,zp,ust,wst,h,dens,ddens,sigmaw,dsigmawdz,tlw,dcas,dcas1!,ran3,gasdev
24  real :: w3,w2
25  real ::  z, &   
26       skew, &
27       skew2, &
28       radw2, &
29       fluarw,fluarw2, &
30       rluarw, &
31       xluarw, &
32       aluarw, &
33       bluarw, &
34       sigmawa, &
35       sigmawb, & 
36       ath, &
37       bth, &
38       wb,wa
39  real timedir
40  real ol, transition
41
42!--------------------------------------------------------------------------- 
43  timedir=ldirect !direction of time forward (1) or backward(-1)
44  z=zp/h
45
46
47  transition=1.
48  if (-h/ol.lt.15) transition=((sin((((-h/ol)+10.)/10.)*pi)))/2.+0.5  !see also in cbl.f90
49
50  w2=sigmaw*sigmaw
51  w3=(((1.2*z*((1.-z)**(3./2.)))+eps)*wst**3) *transition
52
53  skew=w3/(w2**1.5)
54  skew2=skew*skew
55
56  radw2=sqrt(w2) !sigmaw
57
58  fluarw=costluar4*skew**0.333333333333333
59  fluarw2=fluarw*fluarw
60  rluarw=(1.+fluarw2)**3.*skew2/((3.+fluarw2)**2.*fluarw2)  !-> r
61  xluarw=rluarw**0.5 !----> r^1/2
62
63  aluarw=0.5*(1.-xluarw/(4.+rluarw)**0.5)
64  bluarw=1.-aluarw
65
66  sigmawa=radw2*(bluarw/(aluarw*(1.+fluarw2)))**0.5
67  sigmawb=radw2*(aluarw/(bluarw*(1.+fluarw2)))**0.5
68
69  wa=(fluarw*sigmawa)
70  wb=(fluarw*sigmawb)
71
72  dcas=ran3(idum)
73
74  if (dcas.le.aluarw) then
75    dcas1=gasdev(idum)
76    wp=timedir*(dcas1*sigmawa+wa)
77  else
78    dcas1=gasdev(idum)
79    wp=timedir*(dcas1*sigmawb-wb)
80  end if
81
82  return
83end subroutine initialize_cbl_vel
84
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG