source: trunk/src/richardson_gfs.f90 @ 28

Last change on this file since 28 was 4, checked in by mlanger, 9 years ago
File size: 7.9 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  !     CHANGE: 17/11/2005 Caroline Forster NCEP GFS version                  *
45  !                                                                           *
46  !     Two meter level (temperature, humidity) is taken as reference level   *
47  !     instead of first model level.                                         *
48  !     New input variables tt2, td2 introduced.                              *
49  !                                                                           *
50  !****************************************************************************
51  !                                                                           *
52  ! Variables:                                                                *
53  ! h                          mixing height [m]                              *
54  ! hf                         sensible heat flux                             *
55  ! psurf                      surface pressure at point (xt,yt) [Pa]         *
56  ! tv                         virtual temperature                            *
57  ! wst                        convective velocity scale                      *
58  !                                                                           *
59  ! Constants:                                                                *
60  ! ric                        critical Richardson number                     *
61  !                                                                           *
62  !****************************************************************************
63
64  use par_mod
65
66  implicit none
67
68  integer :: i,k,nuvz,iter,llev
69  real :: tv,tvold,zref,z,zold,pint,pold,theta,thetaref,ri
70  real :: akz(nuvz),bkz(nuvz),ulev(nuvz),vlev(nuvz),hf,wst,tt2,td2,ew
71  real :: psurf,ust,ttlev(nuvz),qvlev(nuvz),h,excess
72  real :: thetaold,zl,ul,vl,thetal,ril,hmixplus,wspeed,bvfsq,bvf
73  real :: f_qvsat,rh,rhold,rhl,theta1,theta2,zl1,zl2,thetam
74  real,parameter :: const=r_air/ga, ric=0.25, b=100., bs=8.5
75  integer,parameter :: itmax=3
76
77  excess=0.0
78  iter=0
79
80  ! NCEP version: find first model level above ground
81  !**************************************************
82
83      llev = 0
84      do i=1,nuvz
85       if (psurf.lt.akz(i)) llev=i
86      end do
87       llev = llev+1
88  ! sec llev should not be 1!
89       if (llev.eq.1) llev = 2
90       if (llev.gt.nuvz) llev = nuvz-1
91  ! NCEP version
92
93
94  ! Compute virtual temperature and virtual potential temperature at
95  ! reference level (2 m)
96  !*****************************************************************
97
9830   iter=iter+1
99
100  pold=psurf
101  tvold=tt2*(1.+0.378*ew(td2)/psurf)
102  zold=2.0
103  zref=zold
104  rhold=ew(td2)/ew(tt2)
105
106  thetaref=tvold*(100000./pold)**(r_air/cpa)+excess
107  thetaold=thetaref
108
109  ! Integrate z up to one level above zt
110  !*************************************
111
112  do k=llev,nuvz
113    pint=akz(k)+bkz(k)*psurf  ! pressure on model layers
114    tv=ttlev(k)*(1.+0.608*qvlev(k))
115
116    if (abs(tv-tvold).gt.0.2) then
117      z=zold+const*log(pold/pint)*(tv-tvold)/log(tv/tvold)
118    else
119      z=zold+const*log(pold/pint)*tv
120    endif
121
122    theta=tv*(100000./pint)**(r_air/cpa)
123  ! Petra
124    rh = qvlev(k) / f_qvsat( pint, ttlev(k) )
125
126
127  !alculate Richardson number at each level
128  !****************************************
129
130    ri=ga/thetaref*(theta-thetaref)*(z-zref)/ &
131         max(((ulev(k)-ulev(2))**2+(vlev(k)-vlev(2))**2+b*ust**2),0.1)
132
133  !  addition of second condition: MH should not be placed in an
134  !  unstable layer (PS / Feb 2000)
135    if (ri.gt.ric .and. thetaold.lt.theta) goto 20
136
137    tvold=tv
138    pold=pint
139    rhold=rh
140    thetaold=theta
141    zold=z
142  end do
143
14420   continue
145
146  ! Determine Richardson number between the critical levels
147  !********************************************************
148
149  zl1=zold
150  theta1=thetaold
151  do i=1,20
152    zl=zold+real(i)/20.*(z-zold)
153    ul=ulev(k-1)+real(i)/20.*(ulev(k)-ulev(k-1))
154    vl=vlev(k-1)+real(i)/20.*(vlev(k)-vlev(k-1))
155    thetal=thetaold+real(i)/20.*(theta-thetaold)
156    rhl=rhold+real(i)/20.*(rh-rhold)
157    ril=ga/thetaref*(thetal-thetaref)*(zl-zref)/ &
158         max(((ul-ulev(2))**2+(vl-vlev(2))**2+b*ust**2),0.1)
159    zl2=zl
160    theta2=thetal
161    if (ril.gt.ric) goto 25
162    zl1=zl
163    theta1=thetal
164  end do
165
16625   continue
167  h=zl
168  thetam=0.5*(theta1+theta2)
169  wspeed=sqrt(ul**2+vl**2)                    ! Wind speed at z=hmix
170  bvfsq=(ga/thetam)*(theta2-theta1)/(zl2-zl1) ! Brunt-Vaisala frequency
171                                              ! at z=hmix
172
173  ! Under stable conditions, limit the maximum effect of the subgrid-scale topography
174  ! by the maximum lifting possible from the available kinetic energy
175  !*****************************************************************************
176
177  if(bvfsq.le.0.) then
178    hmixplus=9999.
179  else
180    bvf=sqrt(bvfsq)
181    hmixplus=wspeed/bvf*convke
182  endif
183
184
185  ! Calculate convective velocity scale
186  !************************************
187
188  if (hf.lt.0.) then
189    wst=(-h*ga/thetaref*hf/cpa)**0.333
190    excess=-bs*hf/cpa/wst
191    if (iter.lt.itmax) goto 30
192  else
193    wst=0.
194  endif
195
196end subroutine richardson
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG