source: branches/flexpart91_hasod/src_parallel/richardson.f90 @ 8

Last change on this file since 8 was 8, checked in by hasod, 11 years ago

Added parallel version of Flexpart91

File size: 7.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 richardson(psurf,ust,ttlev,qvlev,ulev,vlev,nuvz, &
23       akz,bkz,hf,tt2,td2,h,wst,hmixplus)
24  !                        i    i    i     i    i    i    i
25  ! i   i  i   i   i  o  o     o
26  !****************************************************************************
27  !                                                                           *
28  !     Calculation of mixing height based on the critical Richardson number. *
29  !     Calculation of convective time scale.                                 *
30  !     For unstable conditions, one iteration is performed. An excess        *
31  !     temperature (dependent on hf and wst) is calculated, added to the     *
32  !     temperature at the lowest model level. Then the procedure is repeated.*
33  !                                                                           *
34  !     Author: A. Stohl                                                      *
35  !                                                                           *
36  !     22 August 1996                                                        *
37  !                                                                           *
38  !     Literature:                                                           *
39  !     Vogelezang DHP and Holtslag AAM (1996): Evaluation and model impacts  *
40  !     of alternative boundary-layer height formulations. Boundary-Layer     *
41  !     Meteor. 81, 245-269.                                                  *
42  !                                                                           *
43  !     Update: 1999-02-01 by G. Wotawa                                       *
44  !                                                                           *
45  !     Two meter level (temperature, humidity) is taken as reference level   *
46  !     instead of first model level.                                         *
47  !     New input variables tt2, td2 introduced.                              *
48  !                                                                           *
49  !****************************************************************************
50  !                                                                           *
51  ! Variables:                                                                *
52  ! h                          mixing height [m]                              *
53  ! hf                         sensible heat flux                             *
54  ! psurf                      surface pressure at point (xt,yt) [Pa]         *
55  ! tv                         virtual temperature                            *
56  ! wst                        convective velocity scale                      *
57  !                                                                           *
58  ! Constants:                                                                *
59  ! ric                        critical Richardson number                     *
60  !                                                                           *
61  !****************************************************************************
62
63  use par_mod
64
65  implicit none
66
67  integer :: i,k,nuvz,iter
68  real :: tv,tvold,zref,z,zold,pint,pold,theta,thetaref,ri
69  real :: akz(nuvz),bkz(nuvz),ulev(nuvz),vlev(nuvz),hf,wst,tt2,td2,ew
70  real :: psurf,ust,ttlev(nuvz),qvlev(nuvz),h,excess
71  real :: thetaold,zl,ul,vl,thetal,ril,hmixplus,wspeed,bvfsq,bvf
72  real :: f_qvsat,rh,rhold,rhl,theta1,theta2,zl1,zl2,thetam
73  real,parameter    :: const=r_air/ga, ric=0.25, b=100., bs=8.5
74  integer,parameter :: itmax=3
75
76  excess=0.0
77  iter=0
78
79  ! Compute virtual temperature and virtual potential temperature at
80  ! reference level (2 m)
81  !*****************************************************************
82
8330   iter=iter+1
84
85  pold=psurf
86  tvold=tt2*(1.+0.378*ew(td2)/psurf)
87  zold=2.0
88  zref=zold
89  rhold=ew(td2)/ew(tt2)
90
91  thetaref=tvold*(100000./pold)**(r_air/cpa)+excess
92  thetaold=thetaref
93
94
95  ! Integrate z up to one level above zt
96  !*************************************
97
98  do k=2,nuvz
99    pint=akz(k)+bkz(k)*psurf  ! pressure on model layers
100    tv=ttlev(k)*(1.+0.608*qvlev(k))
101
102    if (abs(tv-tvold).gt.0.2) then
103      z=zold+const*log(pold/pint)*(tv-tvold)/log(tv/tvold)
104    else
105      z=zold+const*log(pold/pint)*tv
106    endif
107
108    theta=tv*(100000./pint)**(r_air/cpa)
109  ! Petra
110    rh = qvlev(k) / f_qvsat( pint, ttlev(k) )
111
112
113  !alculate Richardson number at each level
114  !****************************************
115
116    ri=ga/thetaref*(theta-thetaref)*(z-zref)/ &
117         max(((ulev(k)-ulev(2))**2+(vlev(k)-vlev(2))**2+b*ust**2),0.1)
118
119  !  addition of second condition: MH should not be placed in an
120  !  unstable layer (PS / Feb 2000)
121    if (ri.gt.ric .and. thetaold.lt.theta) goto 20
122
123    tvold=tv
124    pold=pint
125    rhold=rh
126    thetaold=theta
127    zold=z
128  end do
129
13020   continue
131
132  ! Determine Richardson number between the critical levels
133  !********************************************************
134
135  zl1=zold
136  theta1=thetaold
137  do i=1,20
138    zl=zold+real(i)/20.*(z-zold)
139    ul=ulev(k-1)+real(i)/20.*(ulev(k)-ulev(k-1))
140    vl=vlev(k-1)+real(i)/20.*(vlev(k)-vlev(k-1))
141    thetal=thetaold+real(i)/20.*(theta-thetaold)
142    rhl=rhold+real(i)/20.*(rh-rhold)
143    ril=ga/thetaref*(thetal-thetaref)*(zl-zref)/ &
144         max(((ul-ulev(2))**2+(vl-vlev(2))**2+b*ust**2),0.1)
145    zl2=zl
146    theta2=thetal
147    if (ril.gt.ric) goto 25
148    zl1=zl
149    theta1=thetal
150  end do
151
15225   continue
153  h=zl
154  thetam=0.5*(theta1+theta2)
155  wspeed=sqrt(ul**2+vl**2)                    ! Wind speed at z=hmix
156  bvfsq=(ga/thetam)*(theta2-theta1)/(zl2-zl1) ! Brunt-Vaisala frequency
157                                              ! at z=hmix
158
159  ! Under stable conditions, limit the maximum effect of the subgrid-scale topography
160  ! by the maximum lifting possible from the available kinetic energy
161  !*****************************************************************************
162
163  if(bvfsq.le.0.) then
164    hmixplus=9999.
165  else
166    bvf=sqrt(bvfsq)
167    hmixplus=wspeed/bvf*convke
168  endif
169
170
171  ! Calculate convective velocity scale
172  !************************************
173
174  if (hf.lt.0.) then
175    wst=(-h*ga/thetaref*hf/cpa)**0.333
176    excess=-bs*hf/cpa/wst
177    if (iter.lt.itmax) goto 30
178  else
179    wst=0.
180  endif
181
182end subroutine richardson
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG