23 akz,bkz,hf,tt2,td2,h,wst,hmixplus)
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
86 tvold=tt2*(1.+0.378*
ew(td2)/psurf)
91 thetaref=tvold*(100000./pold)**(r_air/cpa)+excess
99 pint=akz(k)+bkz(k)*psurf
100 tv=ttlev(k)*(1.+0.608*qvlev(k))
102 if (abs(tv-tvold).gt.0.2)
then
103 z=zold+const*log(pold/pint)*(tv-tvold)/log(tv/tvold)
105 z=zold+const*log(pold/pint)*tv
108 theta=tv*(100000./pint)**(r_air/cpa)
110 rh = qvlev(k) /
f_qvsat( pint, ttlev(k) )
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)
121 if (ri.gt.ric .and. thetaold.lt.theta) goto 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)
147 if (ril.gt.ric) goto 25
154 thetam=0.5*(theta1+theta2)
155 wspeed=sqrt(ul**2+vl**2)
156 bvfsq=(ga/thetam)*(theta2-theta1)/(zl2-zl1)
167 hmixplus=wspeed/bvf*convke
175 wst=(-h*ga/thetaref*hf/cpa)**0.333
176 excess=-bs*hf/cpa/wst
177 if (iter.lt.itmax) goto 30
subroutine richardson_ecmwf(psurf, ust, ttlev, qvlev, ulev, vlev, nuvz, akz, bkz, hf, tt2, td2, h, wst, hmixplus)
real function f_qvsat(p, t)