source: flexpart.git/src/richardson.f90

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

add SPDX-License-Identifier to all .f90 files

  • Property mode set to 100644
File size: 7.6 KB
Line 
1! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt
2! SPDX-License-Identifier: GPL-3.0-or-later
3
4subroutine richardson(psurf,ust,ttlev,qvlev,ulev,vlev,nuvz, &
5       akz,bkz,hf,tt2,td2,h,wst,hmixplus,metdata_format)
6  !                        i    i    i     i    i    i    i
7  ! i   i  i   i   i  o  o     o
8  !****************************************************************************
9  !                                                                           *
10  !     Calculation of mixing height based on the critical Richardson number. *
11  !     Calculation of convective time scale.                                 *
12  !     For unstable conditions, one iteration is performed. An excess        *
13  !     temperature (dependent on hf and wst) is calculated, added to the     *
14  !     temperature at the lowest model level. Then the procedure is repeated.*
15  !                                                                           *
16  !     Author: A. Stohl                                                      *
17  !                                                                           *
18  !     22 August 1996                                                        *
19  !                                                                           *
20  !     Literature:                                                           *
21  !     Vogelezang DHP and Holtslag AAM (1996): Evaluation and model impacts  *
22  !     of alternative boundary-layer height formulations. Boundary-Layer     *
23  !     Meteor. 81, 245-269.                                                  *
24  !                                                                           *
25  !****************************************************************************
26  !                                                                           *
27  !     Update: 1999-02-01 by G. Wotawa                                       *
28  !                                                                           *
29  !     Two meter level (temperature, humidity) is taken as reference level   *
30  !     instead of first model level.                                         *
31  !     New input variables tt2, td2 introduced.                              *
32  !                                                                           *
33  !     CHANGE: 17/11/2005 Caroline Forster NCEP GFS version                  *
34  !                                                                           *
35  !     Unified ECMWF and GFS builds                                          *
36  !     Marian Harustak, 12.5.2017                                            *
37  !       - Merged richardson and richardson_gfs into one routine using       *
38  !         if-then for meteo-type dependent code                             *
39  !                                                                           *
40  !****************************************************************************
41  !                                                                           *
42  ! Variables:                                                                *
43  ! h                          mixing height [m]                              *
44  ! hf                         sensible heat flux                             *
45  ! psurf                      surface pressure at point (xt,yt) [Pa]         *
46  ! tv                         virtual temperature                            *
47  ! wst                        convective velocity scale                      *
48  ! metdata_format             format of metdata (ecmwf/gfs)                  *
49  !                                                                           *
50  ! Constants:                                                                *
51  ! ric                        critical Richardson number                     *
52  !                                                                           *
53  !****************************************************************************
54
55  use par_mod
56  use class_gribfile
57
58  implicit none
59
60  integer :: metdata_format
61  integer :: i,k,nuvz,iter,llev,loop_start
62  real :: tv,tvold,zref,z,zold,pint,pold,theta,thetaref,ri
63  real :: akz(nuvz),bkz(nuvz),ulev(nuvz),vlev(nuvz),hf,wst,tt2,td2,ew
64  real :: psurf,ust,ttlev(nuvz),qvlev(nuvz),h,excess
65  real :: thetaold,zl,ul,vl,thetal,ril,hmixplus,wspeed,bvfsq,bvf
66  real :: f_qvsat,rh,rhold,rhl,theta1,theta2,zl1,zl2,thetam
67  real,parameter    :: const=r_air/ga, ric=0.25, b=100., bs=8.5
68  integer,parameter :: itmax=3
69
70  excess=0.0
71  iter=0
72
73  if (metdata_format.eq.GRIBFILE_CENTRE_NCEP) then
74    ! NCEP version: find first model level above ground
75    !**************************************************
76
77     llev = 0
78     do i=1,nuvz
79       if (psurf.lt.akz(i)) llev=i
80     end do
81     llev = llev+1
82    ! sec llev should not be 1!
83     if (llev.eq.1) llev = 2
84     if (llev.gt.nuvz) llev = nuvz-1
85    ! NCEP version
86  end if
87
88
89  ! Compute virtual temperature and virtual potential temperature at
90  ! reference level (2 m)
91  !*****************************************************************
92
9330   iter=iter+1
94
95  pold=psurf
96  tvold=tt2*(1.+0.378*ew(td2)/psurf)
97  zold=2.0
98  zref=zold
99  rhold=ew(td2)/ew(tt2)
100
101  thetaref=tvold*(100000./pold)**(r_air/cpa)+excess
102  thetaold=thetaref
103
104
105  ! Integrate z up to one level above zt
106  !*************************************
107  if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
108    loop_start=2
109  else
110    loop_start=llev
111  end if
112  do k=loop_start,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  ! Calculate 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  k=k-1 ! ESO: make sure k <= nuvz (ticket #139)
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