Index: src/richardson_gfs.f90
===================================================================
--- src/richardson_gfs.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ src/richardson_gfs.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
@@ -0,0 +1,196 @@
+!**********************************************************************
+! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
+! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
+! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
+!                                                                     *
+! This file is part of FLEXPART.                                      *
+!                                                                     *
+! FLEXPART is free software: you can redistribute it and/or modify    *
+! it under the terms of the GNU General Public License as published by*
+! the Free Software Foundation, either version 3 of the License, or   *
+! (at your option) any later version.                                 *
+!                                                                     *
+! FLEXPART is distributed in the hope that it will be useful,         *
+! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
+! GNU General Public License for more details.                        *
+!                                                                     *
+! You should have received a copy of the GNU General Public License   *
+! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
+!**********************************************************************
+
+subroutine richardson(psurf,ust,ttlev,qvlev,ulev,vlev,nuvz, &
+       akz,bkz,hf,tt2,td2,h,wst,hmixplus)
+  !                        i    i    i     i    i    i    i
+  ! i   i  i   i   i  o  o     o
+  !****************************************************************************
+  !                                                                           *
+  !     Calculation of mixing height based on the critical Richardson number. *
+  !     Calculation of convective time scale.                                 *
+  !     For unstable conditions, one iteration is performed. An excess        *
+  !     temperature (dependent on hf and wst) is calculated, added to the     *
+  !     temperature at the lowest model level. Then the procedure is repeated.*
+  !                                                                           *
+  !     Author: A. Stohl                                                      *
+  !                                                                           *
+  !     22 August 1996                                                        *
+  !                                                                           *
+  !     Literature:                                                           *
+  !     Vogelezang DHP and Holtslag AAM (1996): Evaluation and model impacts  *
+  !     of alternative boundary-layer height formulations. Boundary-Layer     *
+  !     Meteor. 81, 245-269.                                                  *
+  !                                                                           *
+  !     Update: 1999-02-01 by G. Wotawa                                       *
+  !     CHANGE: 17/11/2005 Caroline Forster NCEP GFS version                  *
+  !                                                                           *
+  !     Two meter level (temperature, humidity) is taken as reference level   *
+  !     instead of first model level.                                         *
+  !     New input variables tt2, td2 introduced.                              *
+  !                                                                           *
+  !****************************************************************************
+  !                                                                           *
+  ! Variables:                                                                *
+  ! h                          mixing height [m]                              *
+  ! hf                         sensible heat flux                             *
+  ! psurf                      surface pressure at point (xt,yt) [Pa]         *
+  ! tv                         virtual temperature                            *
+  ! wst                        convective velocity scale                      *
+  !                                                                           *
+  ! Constants:                                                                *
+  ! ric                        critical Richardson number                     *
+  !                                                                           *
+  !****************************************************************************
+
+  use par_mod
+
+  implicit none
+
+  integer :: i,k,nuvz,iter,llev
+  real :: tv,tvold,zref,z,zold,pint,pold,theta,thetaref,ri
+  real :: akz(nuvz),bkz(nuvz),ulev(nuvz),vlev(nuvz),hf,wst,tt2,td2,ew
+  real :: psurf,ust,ttlev(nuvz),qvlev(nuvz),h,excess
+  real :: thetaold,zl,ul,vl,thetal,ril,hmixplus,wspeed,bvfsq,bvf
+  real :: f_qvsat,rh,rhold,rhl,theta1,theta2,zl1,zl2,thetam
+  real,parameter :: const=r_air/ga, ric=0.25, b=100., bs=8.5
+  integer,parameter :: itmax=3
+
+  excess=0.0
+  iter=0
+
+  ! NCEP version: find first model level above ground
+  !**************************************************
+
+      llev = 0
+      do i=1,nuvz
+       if (psurf.lt.akz(i)) llev=i
+      end do
+       llev = llev+1
+  ! sec llev should not be 1!
+       if (llev.eq.1) llev = 2
+       if (llev.gt.nuvz) llev = nuvz-1
+  ! NCEP version
+
+
+  ! Compute virtual temperature and virtual potential temperature at
+  ! reference level (2 m)
+  !*****************************************************************
+
+30   iter=iter+1
+
+  pold=psurf
+  tvold=tt2*(1.+0.378*ew(td2)/psurf)
+  zold=2.0
+  zref=zold
+  rhold=ew(td2)/ew(tt2)
+
+  thetaref=tvold*(100000./pold)**(r_air/cpa)+excess
+  thetaold=thetaref
+
+  ! Integrate z up to one level above zt
+  !*************************************
+
+  do k=llev,nuvz
+    pint=akz(k)+bkz(k)*psurf  ! pressure on model layers
+    tv=ttlev(k)*(1.+0.608*qvlev(k))
+
+    if (abs(tv-tvold).gt.0.2) then
+      z=zold+const*log(pold/pint)*(tv-tvold)/log(tv/tvold)
+    else
+      z=zold+const*log(pold/pint)*tv
+    endif
+
+    theta=tv*(100000./pint)**(r_air/cpa)
+  ! Petra
+    rh = qvlev(k) / f_qvsat( pint, ttlev(k) )
+
+
+  !alculate Richardson number at each level
+  !****************************************
+
+    ri=ga/thetaref*(theta-thetaref)*(z-zref)/ &
+         max(((ulev(k)-ulev(2))**2+(vlev(k)-vlev(2))**2+b*ust**2),0.1)
+
+  !  addition of second condition: MH should not be placed in an
+  !  unstable layer (PS / Feb 2000)
+    if (ri.gt.ric .and. thetaold.lt.theta) goto 20
+
+    tvold=tv
+    pold=pint
+    rhold=rh
+    thetaold=theta
+    zold=z
+  end do
+
+20   continue
+
+  ! Determine Richardson number between the critical levels
+  !********************************************************
+
+  zl1=zold
+  theta1=thetaold
+  do i=1,20
+    zl=zold+real(i)/20.*(z-zold)
+    ul=ulev(k-1)+real(i)/20.*(ulev(k)-ulev(k-1))
+    vl=vlev(k-1)+real(i)/20.*(vlev(k)-vlev(k-1))
+    thetal=thetaold+real(i)/20.*(theta-thetaold)
+    rhl=rhold+real(i)/20.*(rh-rhold)
+    ril=ga/thetaref*(thetal-thetaref)*(zl-zref)/ &
+         max(((ul-ulev(2))**2+(vl-vlev(2))**2+b*ust**2),0.1)
+    zl2=zl
+    theta2=thetal
+    if (ril.gt.ric) goto 25
+    zl1=zl
+    theta1=thetal
+  end do
+
+25   continue
+  h=zl
+  thetam=0.5*(theta1+theta2)
+  wspeed=sqrt(ul**2+vl**2)                    ! Wind speed at z=hmix
+  bvfsq=(ga/thetam)*(theta2-theta1)/(zl2-zl1) ! Brunt-Vaisala frequency
+                                              ! at z=hmix
+
+  ! Under stable conditions, limit the maximum effect of the subgrid-scale topography
+  ! by the maximum lifting possible from the available kinetic energy
+  !*****************************************************************************
+
+  if(bvfsq.le.0.) then
+    hmixplus=9999.
+  else
+    bvf=sqrt(bvfsq)
+    hmixplus=wspeed/bvf*convke
+  endif
+
+
+  ! Calculate convective velocity scale
+  !************************************
+
+  if (hf.lt.0.) then
+    wst=(-h*ga/thetaref*hf/cpa)**0.333
+    excess=-bs*hf/cpa/wst
+    if (iter.lt.itmax) goto 30
+  else
+    wst=0.
+  endif
+
+end subroutine richardson
