23 akz,bkz,hf,tt2,td2,h,wst,hmixplus)
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
85 if (psurf.lt.akz(i)) llev=i
89 if (llev.eq.1) llev = 2
90 if (llev.gt.nuvz) llev = nuvz-1
101 tvold=tt2*(1.+0.378*
ew(td2)/psurf)
104 rhold=
ew(td2)/
ew(tt2)
106 thetaref=tvold*(100000./pold)**(r_air/cpa)+excess
113 pint=akz(k)+bkz(k)*psurf
114 tv=ttlev(k)*(1.+0.608*qvlev(k))
116 if (abs(tv-tvold).gt.0.2)
then
117 z=zold+const*log(pold/pint)*(tv-tvold)/log(tv/tvold)
119 z=zold+const*log(pold/pint)*tv
122 theta=tv*(100000./pint)**(r_air/cpa)
124 rh = qvlev(k) /
f_qvsat( pint, ttlev(k) )
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)
135 if (ri.gt.ric .and. thetaold.lt.theta) goto 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)
161 if (ril.gt.ric) goto 25
168 thetam=0.5*(theta1+theta2)
169 wspeed=sqrt(ul**2+vl**2)
170 bvfsq=(ga/thetam)*(theta2-theta1)/(zl2-zl1)
181 hmixplus=wspeed/bvf*convke
189 wst=(-h*ga/thetaref*hf/cpa)**0.333
190 excess=-bs*hf/cpa/wst
191 if (iter.lt.itmax) goto 30
subroutine richardson_gfs(psurf, ust, ttlev, qvlev, ulev, vlev, nuvz, akz, bkz, hf, tt2, td2, h, wst, hmixplus)
real function f_qvsat(p, t)