!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE PREPRADAR ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE prepradar(nx,ny,nz,nz_tab,nvar_anx,nvar_radin, & 1,3
nvar_rad,mx_rad,nsrc_rad,nz_rdr,mx_colrad,mx_pass, &
raduvobs,radrhobs,radistride,radkstride, &
iuserad,npass,refrh,rhradobs, &
xs,ys,zs,hterain,anx,qback,hqback,nlvqback, &
nradfil,radfname,srcrad,isrcrad,qsrcrad,qcthrrad, &
stnrad,latrad,lonrad,elvrad, &
latradc,lonradc,xradc,yradc,irad,nlevrad, &
distrad,uazmrad,vazmrad,hgtradc,theradc, &
obsrad,oanxrad,odifrad,qobsrad,qualrad, &
ncolrad,istatus,tem1,tem2,tem3,tem4,tem5,tem6)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Read radar data stored as columns on a remapped grid
! and prepare data for use in the analysis.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Keith Brewster
! August, 1995
!
! MODIFICATION HISTORY:
!
! 07/15/97 (K. Brewster)
! Corrected a bug in interpolation weights.
!
! Nov 1997 (KB)
! Added rhmax,refsat and rhrad to argument list.
!
! Feb, 1998 (KB)
! Added addition option switches to argument list allowing
! selection of either radar wind obs and/or radar RH pseudo-obs.
!
! Mar, 1998 (KB)
! Added uazmlim to protect against division by zero when
! azimuth is close to perpendicular to u or v component.
!
! Dec, 1998 (KB)
! Made changes for RHstar to qv as moisture variable.
!
! August, 2001 (KB)
! Corrected call to dhdr which now returns dhdr not local
! elevation angle. Modified calculation of odifrad to speed
! convergence to vr.
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! nz Number of grid points in the vertical
! nvar_anx number of analysis variables
! nvar_radin number of variables in the obsrad array
! nvar_rad number of variables in the odifrad array
! mx_rad maximum number of radars
! nsrc_rad number of radar sources
! nz_rdr maximum number of levels in a radar column
! mx_colrad maximum number of radar columns
!
! refsat Reflectivity threshold for assigning moisture ob
! rhrad Value of relative humidity to assign
!
! xs x location of scalar pts (m)
! ys y location of scalar pts (m)
! zs z location of scalar pts (m)
! hterain height of terrain (m MSL)
! anx analyzed variables (or first guess)
! qback background (first guess) error
!
! nradfil number of radar files
! radfname file name for radar dataset
! srcrad name of radar sources
! qsrcrad radar observation error
!
! OUTPUT:
!
! isrcrad index of radar source
! stnrad radar site name character*5
! latrad latitude of radar (degrees N)
! lonrad longitude of radar (degrees E)
! elvrad elevation of feed horn of radar (m MSL)
! latradc latitude of radar column (degrees N)
! lonradc longitude of radar column (degrees E)
! xradc x location of radar column
! yradc y location of radar column
! irad radar number
! nlevrad number of levels of radar data in each column
! distrad distance of radar column from source radar
! uazmrad u-component of radar beam for each column
! vazmrad v-component of radar beam for each column
! hgtradc height (m MSL) of radar observations
! obsrad radar observations
! oanxrad analysis (first guess) value at radar data location
! odifrad difference between radar observation and analysis
! qobsrad normalized observation error
! qualrad radar data quality indicator
! ncolrad number of radar columns read-in
! istatus status indicator
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz,nz_tab,nvar_anx
INTEGER :: nvar_radin,nvar_rad,mx_rad,nsrc_rad,nz_rdr,mx_colrad
INTEGER :: mx_pass
INTEGER :: raduvobs
INTEGER :: radrhobs
INTEGER :: radistride
INTEGER :: radkstride
INTEGER :: iuserad(0:nsrc_rad,mx_pass)
INTEGER :: npass
REAL :: refrh
REAL :: rhradobs
!
!-----------------------------------------------------------------------
!
! Grid variables
!
!-----------------------------------------------------------------------
!
REAL :: xs(nx)
REAL :: ys(ny)
REAL :: zs(nx,ny,nz)
REAL :: hterain(nx,ny)
REAL :: anx(nx,ny,nz,nvar_anx)
REAL :: qback(nvar_anx,nz_tab)
REAL :: hqback(nz_tab)
INTEGER :: nlvqback
INTEGER :: nradfil
CHARACTER (LEN=132) :: radfname(mx_rad)
!
!-----------------------------------------------------------------------
!
! Radar site variables
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=8) :: srcrad(nsrc_rad)
INTEGER :: isrcrad(0:mx_rad)
REAL :: qsrcrad(nvar_radin,nsrc_rad)
REAL :: qcthrrad(nvar_radin,nsrc_rad)
REAL :: latrad(mx_rad),lonrad(mx_rad)
REAL :: elvrad(mx_rad)
CHARACTER (LEN=5) :: stnrad(mx_rad)
!
!-----------------------------------------------------------------------
!
! Radar observation variables
!
!-----------------------------------------------------------------------
!
INTEGER :: irad(mx_colrad)
INTEGER :: nlevrad(mx_colrad)
REAL :: distrad(mx_colrad),uazmrad(mx_colrad),vazmrad(mx_colrad)
REAL :: latradc(mx_colrad),lonradc(mx_colrad)
REAL :: xradc(mx_colrad),yradc(mx_colrad)
REAL :: hgtradc(nz_rdr,mx_colrad)
REAL :: theradc(nz_rdr,mx_colrad)
REAL :: obsrad(nvar_radin,nz_rdr,mx_colrad)
REAL :: oanxrad(nvar_rad,nz_rdr,mx_colrad)
REAL :: odifrad(nvar_rad,nz_rdr,mx_colrad)
REAL :: qobsrad(nvar_rad,nz_rdr,mx_colrad)
INTEGER :: qualrad(nvar_rad,nz_rdr,mx_colrad)
INTEGER :: ncolrad
INTEGER :: istatus
!
!-----------------------------------------------------------------------
!
! Temporary work arrays
!
!-----------------------------------------------------------------------
!
REAL :: tem1(nx,ny,nz)
REAL :: tem2(nx,ny,nz)
REAL :: tem3(nx,ny,nz)
REAL :: tem4(nx,ny,nz)
REAL :: tem5(nx,ny,nz)
REAL :: tem6(nx,ny,nz)
!
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
!
INTEGER :: iusechk
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!-----------------------------------------------------------------------
!
! Read remapped radar data
!
!-----------------------------------------------------------------------
!
isrcrad(0)=0
IF(raduvobs > 0 .OR. radrhobs > 0) THEN
iusechk=1
CALL rdradcol
(nx,ny,nz,nsrc_rad,nvar_radin, &
mx_rad,nz_rdr,mx_colrad,mx_pass, &
radistride,radkstride, &
iuserad,iusechk,npass,nradfil,radfname, &
srcrad,isrcrad(1),stnrad,latrad,lonrad,elvrad, &
latradc,lonradc,irad,nlevrad,hgtradc,obsrad, &
ncolrad,istatus,tem1,tem2,tem3,tem4,tem5,tem6)
WRITE(6,'(a,i6/)') ' Prepradar: ncolrad = ',ncolrad
IF(ncolrad > 0) THEN
!
!-----------------------------------------------------------------------
!
! Remove precipitation terminal velocity from the radial velocity.
!
!-----------------------------------------------------------------------
!
CALL rmvterm(nvar_radin,mx_rad,nz_rdr,mx_colrad, &
latrad,lonrad,elvrad, &
latradc,lonradc,irad,nlevrad,hgtradc,obsrad, &
ncolrad,istatus)
!
!-----------------------------------------------------------------------
!
! Unfold and QC radar data
!
!-----------------------------------------------------------------------
!
CALL prcssrad(nx,ny,nz,nz_tab,nvar_anx,nvar_radin, &
nvar_rad,mx_rad,nsrc_rad,nz_rdr,mx_colrad, &
raduvobs,radrhobs,refrh,rhradobs, &
xs,ys,zs,hterain,anx,qback,hqback,nlvqback, &
srcrad,qsrcrad,qcthrrad, &
isrcrad,stnrad,latrad,lonrad,elvrad, &
latradc,lonradc,xradc,yradc,irad,nlevrad, &
hgtradc,obsrad,oanxrad,odifrad, &
distrad,uazmrad,vazmrad, &
theradc,qobsrad,qualrad, &
ncolrad,istatus,tem1)
END IF
ELSE
ncolrad=0
END IF
!
RETURN
END SUBROUTINE prepradar
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE PRCSSRAD ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE prcssrad(nx,ny,nz,nz_tab,nvar_anx,nvar_radin, & 1,6
nvar_rad,mx_rad,nsrc_rad,nz_rdr,mx_colrad, &
raduvobs,radrhobs,refrh,rhradobs, &
xs,ys,zs,hterain,anx,qback,hqback,nlvqback, &
srcrad,qsrcrad,qcthrrad, &
isrcrad,stnrad,latrad,lonrad,elvrad, &
latradc,lonradc,xradc,yradc,irad,nlevrad, &
hgtradc,obsrad,oanxrad,odifrad, &
distrad,uazmrad,vazmrad, &
theradc,qobsrad,qualrad, &
ncolrad,istatus,tem1)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Unfolds and QC's radar data stored as columns on a remapped grid.
!
! AUTHOR: Keith Brewster
! August, 1995
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! nz Number of grid points in the vertical
! nvar_anx number of analysis variables
! nvar_radin number of variables in the obsrad array
! nvar_rad number of variables in the odifrad array
! mx_rad maximum number of radars
! nsrc_rad number of radar sources
! nz_rdr maximum number of levels in a radar column
! mx_colrad maximum number of radar columns
!
! xs x location of scalar pts (m)
! ys y location of scalar pts (m)
! zs z location of scalar pts (m)
! hterain height of terrain (m MSL)
! anx analyzed variables (or first guess)
! qback background (first guess) error
! srcrad name of radar sources
! qsrcrad radar observation error
! isrcrad index of radar source
! stnrad radar site name character*5
! latrad latitude of radar (degrees N)
! lonrad longitude of radar (degrees E)
! elvrad elevation of feed horn of radar (m MSL)
! latradc latitude of radar column (degrees N)
! lonradc longitude of radar column (degrees E)
! xradc x location of radar column (m)
! yradc y location of radar column (m)
! irad radar number
! nlevrad number of levels of radar data in each column
! distrad distance of radar column from source radar
! uazmrad u-component of radar beam for each column
! vazmrad v-component of radar beam for each column
! hgtradc height (m MSL) of radar observations
! obsrad radar observations
!
! OUTPUT:
!
! oanxrad analysis (first guess) value at radar data location
! odifrad difference between radar observation and analysis
! theradc theta (potential temperature K) at each radar ob
! qobsrad normalized observation error
! qualrad radar data quality indicator
! ncolrad number of radar columns read-in
! istatus status indicator
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz,nz_tab,nvar_anx
INTEGER :: nvar_radin,nvar_rad,mx_rad,nsrc_rad,nz_rdr,mx_colrad
INTEGER :: raduvobs
INTEGER :: radrhobs
REAL :: refrh
REAL :: rhradobs
!
!-----------------------------------------------------------------------
!
! Grid variables
!
!-----------------------------------------------------------------------
!
REAL :: xs(nx)
REAL :: ys(ny)
REAL :: zs(nx,ny,nz)
REAL :: hterain(nx,ny)
REAL :: anx(nx,ny,nz,nvar_anx)
REAL :: qback(nvar_anx,nz_tab)
REAL :: hqback(nz_tab)
INTEGER :: nlvqback
!
!-----------------------------------------------------------------------
!
! Radar site variables
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=8) :: srcrad(nsrc_rad)
INTEGER :: isrcrad(0:mx_rad)
REAL :: qsrcrad(nvar_radin,nsrc_rad)
REAL :: qcthrrad(nvar_radin,nsrc_rad)
REAL :: latrad(mx_rad),lonrad(mx_rad)
REAL :: elvrad(mx_rad)
CHARACTER (LEN=5) :: stnrad(mx_rad)
!
!-----------------------------------------------------------------------
!
! Radar observation variables
!
!-----------------------------------------------------------------------
!
INTEGER :: irad(mx_colrad)
INTEGER :: nlevrad(mx_colrad)
REAL :: latradc(mx_colrad),lonradc(mx_colrad)
REAL :: xradc(mx_colrad),yradc(mx_colrad)
REAL :: hgtradc(nz_rdr,mx_colrad)
REAL :: theradc(nz_rdr,mx_colrad)
REAL :: obsrad(nvar_radin,nz_rdr,mx_colrad)
REAL :: oanxrad(nvar_rad,nz_rdr,mx_colrad)
REAL :: odifrad(nvar_rad,nz_rdr,mx_colrad)
REAL :: qobsrad(nvar_rad,nz_rdr,mx_colrad)
REAL :: distrad(mx_colrad),uazmrad(mx_colrad),vazmrad(mx_colrad)
INTEGER :: qualrad(nvar_rad,nz_rdr,mx_colrad)
INTEGER :: ncolrad
INTEGER :: istatus
!
!-----------------------------------------------------------------------
!
! Temporary work array
!
!-----------------------------------------------------------------------
!
REAL :: tem1(nx*ny*nz)
!
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
!
INTEGER :: icol,jrad,kr
INTEGER :: i,imid,iloc
INTEGER :: j,jmid,jloc
INTEGER :: k,kmid,kloc,ktab
INTEGER :: kfold
REAL :: xmid,ymid,zmid
REAL :: wx,wy,w11,w12,w21,w22,whi,wlo,wqhi,wqlo
REAL :: ddrot,ugrid,vgrid,prgrid,tgrid,qvgrid,qvsgrid,rhgrid
REAL :: vr_miss,qv_miss
REAL :: azim,dist,vr,vrdiff,avrdif,dz,eleva,range,dhdr,dsdr,dsdrinv
!
!-----------------------------------------------------------------------
!
! Include file
!
!-----------------------------------------------------------------------
!
INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
! QC and conversion parameters
!
! vrqcthr maximum radial velocity diff from background to use
! vrcnvthr lower limit of the dot product between the radial
! direction and the analysis wind component to use data
!
!-----------------------------------------------------------------------
!
REAL :: vrmax,vrqcthr,uazmlim
PARAMETER (vrmax=80., &
vrqcthr=15., & ! maximum vrdiff from background to use
uazmlim=0.087) ! sin(5 degrees)
REAL :: refmax
PARAMETER (refmax = 80.) ! maximum non-missing reflectivity
!
!-----------------------------------------------------------------------
!
! Function f_qvsat and inline directive for Cray PVP
!
!-----------------------------------------------------------------------
!
REAL :: f_qvsat
!fpp$ expand (f_qvsat)
!dir$ inline always f_qvsat
!*$* inline routine (f_qvsat)
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!-----------------------------------------------------------------------
!
! Some initializations
!
!-----------------------------------------------------------------------
!
vr_miss=-999.
qv_miss=-999.
!
!-----------------------------------------------------------------------
!
! Get x and y locations of each radar ob
!
!-----------------------------------------------------------------------
!
CALL lltoxy(ncolrad,1,latradc,lonradc,xradc,yradc)
!
imid=nx/2
xmid=xs(imid)
jmid=ny/2
ymid=ys(jmid)
kmid=nz/2
!
DO icol=1,ncolrad
IF(irad(icol) > 0 .AND. nlevrad(icol) > 0) THEN
IF( xradc(icol) >= xs(1) .AND. xradc(icol) <= xs(nx-1) &
.AND. yradc(icol) >= ys(1) .AND. yradc(icol) <= ys(ny-1)) &
THEN
!
!-----------------------------------------------------------------------
!
! Find column in grid
!
!-----------------------------------------------------------------------
!
IF(xradc(icol) < xmid) THEN
DO i=imid,2,-1
IF(xs(i) < xradc(icol)) EXIT
END DO
iloc=i
ELSE
DO i=imid,nx-1
IF(xs(i) > xradc(icol)) EXIT
END DO
iloc=i-1
END IF
wx = (xradc(icol)-xs(iloc))/ &
(xs(iloc+1)-xs(iloc))
!
IF(yradc(icol) < ymid) THEN
DO j=jmid,2,-1
IF(ys(j) < yradc(icol)) EXIT
END DO
jloc=j
ELSE
DO j=jmid,ny-1
IF(ys(j) > yradc(icol)) EXIT
END DO
jloc=j-1
END IF
wy = (yradc(icol)-ys(jloc))/ &
(ys(jloc+1)-ys(jloc))
!
!-----------------------------------------------------------------------
!
! Determine bilinear interpolation weights
!
!-----------------------------------------------------------------------
!
w22 = wx*wy
w12 = (1.0 - wx) * wy
w21 = wx * (1.0 - wy)
w11 = (1.0 - wx) * (1.0 - wy)
!
!-----------------------------------------------------------------------
!
! Find heading and distance from radar along ground
! Store correlation of azimuth with u and v drections
! in uazmrad and vazmrad, respectively.
!
!-----------------------------------------------------------------------
!
jrad=irad(icol)
CALL disthead(latradc(icol),lonradc(icol), &
latrad(jrad),lonrad(jrad), &
azim,dist)
distrad(icol)=dist
!
CALL ddrotuv(1,lonradc(icol),azim,1.,ddrot, &
uazmrad(icol),vazmrad(icol))
!
!-----------------------------------------------------------------------
!
! Interpolate grid heights in horizontal
!
!-----------------------------------------------------------------------
!
DO k=1,nz-1
tem1(k)=w11*zs( iloc, jloc,k) + &
w21*zs(iloc+1, jloc,k) + &
w12*zs( iloc,jloc+1,k) + &
w22*zs(iloc+1,jloc+1,k)
END DO
!
!-----------------------------------------------------------------------
!
! Loop in height
!
!-----------------------------------------------------------------------
!
DO kr=1,nlevrad(icol)
IF((ABS(obsrad(2,kr,icol)) < vrmax .OR. &
ABS(obsrad(1,kr,icol)) < refmax) .AND. &
hgtradc(kr,icol) >= tem1(1) .AND. &
hgtradc(kr,icol) <= tem1(nz-1)) THEN
dz=hgtradc(kr,icol)-elvrad(irad(icol))
CALL beamelv(dz,dist,eleva,range)
CALL dhdrange(eleva,range,dhdr)
dsdr=SQRT(MAX(0.,(1.-dhdr*dhdr)))
IF(dsdr /= 0.) THEN
dsdrinv=1./dsdr
ELSE
dsdrinv=0.
END IF
!
!-----------------------------------------------------------------------
!
! Find z location in grid
!
!-----------------------------------------------------------------------
!
zmid=tem1(kmid)
IF(hgtradc(kr,icol) < zmid) THEN
DO k=kmid,2,-1
IF(tem1(k) < hgtradc(kr,icol)) EXIT
END DO
kloc=k
ELSE
DO k=kmid,nz-1
IF(tem1(k) > hgtradc(kr,icol)) EXIT
END DO
kloc=k-1
END IF
!
!-----------------------------------------------------------------------
!
! Set z weights
!
!-----------------------------------------------------------------------
!
whi = (hgtradc(kr,icol)-tem1(kloc))/ &
(tem1(kloc+1)-tem1(kloc))
wlo = 1.0 - whi
!
!-----------------------------------------------------------------------
!
! Set z weights for error data
!
!-----------------------------------------------------------------------
!
DO k=2,nlvqback-1
IF(hqback(k) > hgtradc(kr,icol)) EXIT
END DO
ktab=k-1
wqhi = (hgtradc(kr,icol)-hqback(ktab))/ &
(hqback(ktab+1)-hqback(ktab))
wqlo = 1.0 - wqhi
!
!-----------------------------------------------------------------------
!
! Interpolate background field to radar location
!
! Note this is dependent on indices set so that
! u=1,v=2,theta=4
! This can be made more general in a future version.
!
!-----------------------------------------------------------------------
!
ugrid= &
wlo* (w11*anx( iloc, jloc, kloc,1) + &
w21*anx(iloc+1, jloc, kloc,1) + &
w12*anx( iloc,jloc+1, kloc,1) + &
w22*anx(iloc+1,jloc+1, kloc,1)) &
+ whi* (w11*anx( iloc, jloc,kloc+1,1) + &
w21*anx(iloc+1, jloc,kloc+1,1) + &
w12*anx( iloc,jloc+1,kloc+1,1) + &
w22*anx(iloc+1,jloc+1,kloc+1,1))
vgrid= &
wlo* (w11*anx( iloc, jloc, kloc,2) + &
w21*anx(iloc+1, jloc, kloc,2) + &
w12*anx( iloc,jloc+1, kloc,2) + &
w22*anx(iloc+1,jloc+1, kloc,2)) &
+ whi* (w11*anx( iloc, jloc,kloc+1,2) + &
w21*anx(iloc+1, jloc,kloc+1,2) + &
w12*anx( iloc,jloc+1,kloc+1,2) + &
w22*anx(iloc+1,jloc+1,kloc+1,2))
prgrid= &
wlo* (w11*anx( iloc, jloc, kloc,3) + &
w21*anx(iloc+1, jloc, kloc,3) + &
w12*anx( iloc,jloc+1, kloc,3) + &
w22*anx(iloc+1,jloc+1, kloc,3)) &
+ whi* (w11*anx( iloc, jloc,kloc+1,3) + &
w21*anx(iloc+1, jloc,kloc+1,3) + &
w12*anx( iloc,jloc+1,kloc+1,3) + &
w22*anx(iloc+1,jloc+1,kloc+1,3))
theradc(kr,icol)= &
wlo* (w11*anx( iloc, jloc, kloc,4) + &
w21*anx(iloc+1, jloc, kloc,4) + &
w12*anx( iloc,jloc+1, kloc,4) + &
w22*anx(iloc+1,jloc+1, kloc,4)) &
+ whi* (w11*anx( iloc, jloc,kloc+1,4) + &
w21*anx(iloc+1, jloc,kloc+1,4) + &
w12*anx( iloc,jloc+1,kloc+1,4) + &
w22*anx(iloc+1,jloc+1,kloc+1,4))
qvgrid= &
wlo* (w11*anx( iloc, jloc, kloc,5) + &
w21*anx(iloc+1, jloc, kloc,5) + &
w12*anx( iloc,jloc+1, kloc,5) + &
w22*anx(iloc+1,jloc+1, kloc,5)) &
+ whi* (w11*anx( iloc, jloc,kloc+1,5) + &
w21*anx(iloc+1, jloc,kloc+1,5) + &
w12*anx( iloc,jloc+1,kloc+1,5) + &
w22*anx(iloc+1,jloc+1,kloc+1,5))
!
IF(ABS(obsrad(2,kr,icol)) < vrmax .AND. raduvobs > 0) THEN
!
!-----------------------------------------------------------------------
!
! Find earth-relative u,v
!
!-----------------------------------------------------------------------
!
oanxrad(1,kr,icol)=ugrid
oanxrad(2,kr,icol)=vgrid
!
!-----------------------------------------------------------------------
!
! Find radial velocity
! This assumes motions are quasi-horizontal.
!
!-----------------------------------------------------------------------
!
vr=(uazmrad(icol)*ugrid + vazmrad(icol)*vgrid) * dsdr
!
!-----------------------------------------------------------------------
!
! Check for folding
!
!-----------------------------------------------------------------------
!
vrdiff=obsrad(2,kr,icol)-vr
kfold=nint(vrdiff/(2.*obsrad(3,kr,icol)))
obsrad(2,kr,icol)=obsrad(2,kr,icol)- &
kfold*2.*obsrad(3,kr,icol)
vrdiff=obsrad(2,kr,icol)-vr
avrdif=ABS(vrdiff)
!
!-----------------------------------------------------------------------
!
! QC check and calculation of increments.
! QC should be more sophisticated at some point.
!
!-----------------------------------------------------------------------
!
IF(avrdif < qcthrrad(2,isrcrad(irad(icol)))) THEN
!
!-----------------------------------------------------------------------
!
! Find increment in analysis u and v needed to make vrdiff zero
!
!-----------------------------------------------------------------------
!
odifrad(1,kr,icol)=dsdrinv*(uazmrad(icol)*vrdiff)
odifrad(2,kr,icol)=dsdrinv*(vazmrad(icol)*vrdiff)
!
!-----------------------------------------------------------------------
!
! Assign obs error to u according to observation angle
!
!-----------------------------------------------------------------------
!
IF(ABS(uazmrad(icol)) > uazmlim) THEN
qobsrad(1,kr,icol)=qsrcrad(2,isrcrad(irad(icol)))/ &
((uazmrad(icol)*uazmrad(icol))* &
(wqlo*qback(1,ktab)+wqhi*qback(1,ktab+1)))
qualrad(1,kr,icol)=10
ELSE
qobsrad(1,kr,icol)=vr_miss
qualrad(1,kr,icol)=-90
END IF
!
!-----------------------------------------------------------------------
!
! Assign obs error to v according to observation angle
!
!-----------------------------------------------------------------------
!
IF(ABS(vazmrad(icol)) > uazmlim) THEN
qobsrad(2,kr,icol)=qsrcrad(2,isrcrad(irad(icol)))/ &
((vazmrad(icol)*vazmrad(icol))* &
(wqlo*qback(2,ktab)+wqhi*qback(2,ktab+1)))
qualrad(2,kr,icol)=10
ELSE
qobsrad(2,kr,icol)=vr_miss
qualrad(2,kr,icol)=-90
END IF
ELSE
odifrad(1,kr,icol)=vr_miss
odifrad(2,kr,icol)=vr_miss
qobsrad(1,kr,icol)=vr_miss
qobsrad(2,kr,icol)=vr_miss
qualrad(1,kr,icol)=-99
qualrad(2,kr,icol)=-99
END IF
ELSE
odifrad(1,kr,icol)=vr_miss
odifrad(2,kr,icol)=vr_miss
qobsrad(1,kr,icol)=vr_miss
qobsrad(2,kr,icol)=vr_miss
qualrad(1,kr,icol)=-99
qualrad(2,kr,icol)=-99
END IF
!
!-----------------------------------------------------------------------
!
! Where reflectivity is greater than that
! assumed for "saturation" assign an rhs corresponding to
! rh = rhind, but not if the background rh is greater
! than rhradobs (background more saturated).
!
!-----------------------------------------------------------------------
!
IF(radrhobs > 0 .AND. obsrad(1,kr,icol) > refrh .AND. &
obsrad(1,kr,icol) < refmax ) THEN
tgrid=theradc(kr,icol)*(prgrid/p0)**rddcp
qvsgrid=f_qvsat(prgrid,tgrid)
rhgrid=qvgrid/qvsgrid
IF( rhgrid < rhradobs ) THEN
obsrad(1,kr,icol) = rhradobs*qvsgrid
odifrad(3,kr,icol) = obsrad(1,kr,icol)-qvgrid
oanxrad(3,kr,icol) = qvgrid
qualrad(3,kr,icol) = 10
qobsrad(3,kr,icol) = qsrcrad(1,isrcrad(irad(icol)))* &
qvsgrid
ELSE
obsrad(1,kr,icol) = qv_miss
odifrad(3,kr,icol) = qv_miss
oanxrad(3,kr,icol) = qv_miss
qobsrad(3,kr,icol) = qv_miss
qualrad(3,kr,icol) =-99
END IF
ELSE
obsrad(1,kr,icol) = qv_miss
odifrad(3,kr,icol) = qv_miss
oanxrad(3,kr,icol) = qv_miss
qobsrad(3,kr,icol) = qv_miss
qualrad(3,kr,icol) =-99
END IF
ELSE
odifrad(1,kr,icol)=vr_miss
odifrad(2,kr,icol)=vr_miss
qobsrad(1,kr,icol)=vr_miss
qobsrad(2,kr,icol)=vr_miss
qualrad(1,kr,icol)=-99
qualrad(2,kr,icol)=-99
obsrad(1,kr,icol) = qv_miss
odifrad(3,kr,icol) = qv_miss
oanxrad(3,kr,icol) = qv_miss
qobsrad(3,kr,icol) = qv_miss
qualrad(3,kr,icol) =-99
END IF
!
END DO
ELSE
irad(icol)=0
END IF
!
END IF
END DO
RETURN
END SUBROUTINE prcssrad
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE RADMCRO ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE radmcro(nx,ny,nz, & 1,5
mx_rad,nsrc_rad,nvar_radin,nz_rdr,mx_colrad,mx_pass, &
radistride,radkstride,iuserad,npass, &
nradfil,radfname,srcrad,isrcrad,stnrad, &
latrad,lonrad,elvrad,latradc,lonradc,irad, &
nlevrad,xradc,yradc,hgtradc,obsrad,ncolrad, &
radqvopt,radqcopt,radqropt,radptopt, &
refsat,rhradobs, &
refcld,cldrad,ceilopt,ceilmin,dzfill, &
refrain,radsetrat,radreflim,radptgain, &
xs,ys,zs,zp,j3,pr,pt,qv, &
ptbar,qvbar,rhostr, &
qc,qr,qi,qs,qh,ceillim, &
tem1,tem2,tem3,tem4,tem5,tem6)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Assign values to microphysical variables based on observed
! radar data.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Keith Brewster
! November, 1995
!
! MODIFICATION HISTORY:
! December, 1997
! February, 1998
! September, 1998 Jingsing Zong
! Corrected logic after DO 150 to continue to next column.
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! nx Number of grid points in the x-direction (east/west)
! ny Number of grid points in the y-direction (north/south)
! nz Number of grid points in the vertical
! nvar_radin number of variables in the obsrad array
! nz_rdr maximum number of levels in a radar column
! mx_colrad maximum number of radar columns
!
! xradc x location of radar column
! yradc y location of radar column
! hgtradc height (m MSL) of radar observations
! obsrad radar observations
! ncolrad number of radar columns read-in
!
! xs x location of scalar pts (m)
! ys y location of scalar pts (m)
! zs z location of scalar pts (m)
!
! OUTPUT:
!
! qc cloud water mixing ratio
!
! WORK ARRAYS
!
! ceillim Ceiling limit (m MSL)
! tem1 Temporary array
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz
INTEGER :: mx_rad,nsrc_rad,nvar_radin,nz_rdr,mx_colrad
INTEGER :: mx_pass
INTEGER :: radistride
INTEGER :: radkstride
INTEGER :: iuserad(0:nsrc_rad,mx_pass)
INTEGER :: npass
!
!-----------------------------------------------------------------------
!
! Control parameters
!
!-----------------------------------------------------------------------
!
INTEGER :: radqvopt,radqcopt,radqropt,radptopt
REAL :: refsat,rhradobs,refcld,cldrad
INTEGER :: ceilopt
REAL :: ceilmin,dzfill
REAL :: refrain
REAL :: radsetrat,radreflim,radptgain
!
!-----------------------------------------------------------------------
!
! Radar observation variables
!
!-----------------------------------------------------------------------
!
INTEGER :: nradfil
CHARACTER (LEN=132) :: radfname(mx_rad)
CHARACTER (LEN=8) :: srcrad(nsrc_rad)
INTEGER :: isrcrad(0:mx_rad)
CHARACTER (LEN=5) :: stnrad(mx_rad)
REAL :: latrad(mx_rad)
REAL :: lonrad(mx_rad)
REAL :: elvrad(mx_rad)
REAL :: latradc(mx_colrad)
REAL :: lonradc(mx_colrad)
REAL :: xradc(mx_colrad)
REAL :: yradc(mx_colrad)
INTEGER :: irad(mx_colrad)
INTEGER :: nlevrad(mx_colrad)
REAL :: hgtradc(nz_rdr,mx_colrad)
REAL :: obsrad(nvar_radin,nz_rdr,mx_colrad)
INTEGER :: ncolrad
!
REAL :: xs(nx)
REAL :: ys(ny)
REAL :: zs(nx,ny,nz)
REAL :: zp(nx,ny,nz)
REAL :: j3(nx,ny,nz)
REAL :: pr(nx,ny,nz)
REAL :: pt(nx,ny,nz)
REAL :: qv(nx,ny,nz)
REAL :: ptbar(nx,ny,nz)
REAL :: qvbar(nx,ny,nz)
REAL :: rhostr(nx,ny,nz)
REAL :: qc(nx,ny,nz)
REAL :: qr(nx,ny,nz)
REAL :: qs(nx,ny,nz)
REAL :: qi(nx,ny,nz)
REAL :: qh(nx,ny,nz)
REAL :: ceillim(nx,ny)
REAL :: tem1(nx,ny,nz)
REAL :: tem2(nx,ny,nz)
REAL :: tem3(nx,ny,nz)
REAL :: tem4(nx,ny,nz)
REAL :: tem5(nx,ny,nz)
REAL :: tem6(nx,ny,nz)
!
REAL :: reflim
PARAMETER (reflim=90.)
!
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,k,jrad,icol,klow,knex,klast
REAL :: dx,dy,dstlim2,dist
REAL :: ref,reflow,qsh
INTEGER :: iusechk
INTEGER :: istatus
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
dx=xs(2)-xs(1)
dy=ys(2)-ys(1)
dstlim2=0.95*(dx*dx+dy*dy)
PRINT *, ' dstlim2 = ',dstlim2
!
!-----------------------------------------------------------------------
!
! Because some radar variables may have been overwritten
! in previous routines, re-read radar files.
!
!-----------------------------------------------------------------------
!
iusechk=0
CALL rdradcol(nx,ny,nz,nsrc_rad,nvar_radin, &
mx_rad,nz_rdr,mx_colrad,mx_pass, &
radistride,radkstride, &
iuserad,iusechk,npass,nradfil,radfname, &
srcrad,isrcrad(1),stnrad,latrad,lonrad,elvrad, &
latradc,lonradc,irad,nlevrad,hgtradc,obsrad, &
ncolrad,istatus,tem1,tem2,tem3,tem4,tem5,tem6)
!
!-----------------------------------------------------------------------
!
! Get x and y locations of each radar ob
!
!-----------------------------------------------------------------------
!
CALL lltoxy(ncolrad,1,latradc,lonradc,xradc,yradc)
!-----------------------------------------------------------------------
!
! Set minimum ceiling, stored in ceillim
!
!-----------------------------------------------------------------------
!
IF(ceilopt == 1) THEN
WRITE(6,'(a,f6.0,a)') &
' Setting minimun ceiling to ',ceilmin,' m AGL'
DO j=1,ny-1
DO i=1,nx-1
ceillim(i,j)=ceilmin+zp(i,j,2)
END DO
END DO
ELSE IF(ceilopt == 2) THEN
WRITE(6,'(a)') &
' Setting minimun ceiling based on surface LCL'
CALL adaslcl(nx,ny,nz,zs,pr,pt,qv,ceillim)
ELSE
WRITE(6,'(a)') &
' Setting zero AGL minimun ceiling'
DO j=1,ny-1
DO i=1,nx-1
ceillim(i,j)=zp(i,j,2)
END DO
END DO
END IF
DO j=1,ny-1
DO i=1,nx-1
!
!-----------------------------------------------------------------------
!
! For each horizontal grid point, use the radar column(s)
! within the threshold distance defined by dstlim2.
!
!-----------------------------------------------------------------------
!
DO icol=1,ncolrad
jrad=irad(icol)
IF(isrcrad(jrad) > 0) THEN
dist=((xradc(icol)-xs(i))*(xradc(icol)-xs(i)) + &
(yradc(icol)-ys(j))*(yradc(icol)-ys(j)))
IF(dist < dstlim2) THEN
DO klow=1,nlevrad(icol)
IF(ABS(obsrad(1,klow,icol)) < reflim .AND. &
obsrad(1,klow,icol) >= refcld ) GO TO 151
END DO
CYCLE
151 CONTINUE
!
!-----------------------------------------------------------------------
!
! Found a reflectivity indicating at least cloud must be
! present. Find index in ARPS grid just above the height of
! this observation.
!
!-----------------------------------------------------------------------
!
180 CONTINUE
DO k=1,nz-1
IF(zs(i,j,k) >= hgtradc(klow,icol)) GO TO 201
END DO
CYCLE
201 CONTINUE
IF(zs(i,j,k) >= ceillim(i,j)) THEN
ref=obsrad(1,klow,icol)
qsh=qs(i,j,k)+qh(i,j,k)
CALL cldset(pr(i,j,k),pt(i,j,k), &
j3(i,j,k),ptbar(i,j,k),qvbar(i,j,k),rhostr(i,j,k), &
radqvopt,radqcopt,radqropt,radptopt, &
refsat,rhradobs,refcld,cldrad, &
refrain,radsetrat,radreflim,radptgain, &
ref,qc(i,j,k),qr(i,j,k),qi(i,j,k),qsh)
END IF
reflow=obsrad(1,klow,icol)
klast=k
!
!-----------------------------------------------------------------------
!
! Find next reflectivity indicating at least cloud must be
! present.
!
!-----------------------------------------------------------------------
!
220 CONTINUE
DO knex=klow+1,nlevrad(icol)
IF(ABS(obsrad(1,knex,icol)) < reflim .AND. &
obsrad(1,knex,icol) >= refcld ) GO TO 251
END DO
CYCLE
251 CONTINUE
!
!-----------------------------------------------------------------------
!
! Found another reflectivity, so determine which
! ARPS grid points are affected. Either
! 1) Fill from the previous reflectivity level ..or..
! 2) Set a new klow and klast and treat as if this were
! the first point for this column.
!
!-----------------------------------------------------------------------
!
IF((hgtradc(knex,icol)-hgtradc(klow,icol)) < &
dzfill) THEN
DO k=klast+1,nz-1
IF(zs(i,j,k) <= hgtradc(knex,icol)) THEN
IF(zs(i,j,k) > ceillim(i,j)) THEN
ref=(obsrad(1,knex,icol)-reflow) &
/(hgtradc(knex,icol)-hgtradc(klow,icol)) &
*(zs(i,j,k)-hgtradc(klow,icol)) &
+reflow
qsh=qs(i,j,k)+qh(i,j,k)
CALL cldset(pr(i,j,k),pt(i,j,k), &
j3(i,j,k),ptbar(i,j,k), &
qvbar(i,j,k),rhostr(i,j,k), &
radqvopt,radqcopt,radqropt,radptopt, &
refsat,rhradobs,refcld,cldrad, &
refrain,radsetrat,radreflim,radptgain, &
ref,qc(i,j,k),qr(i,j,k),qi(i,j,k),qsh)
END IF
klast=k
ELSE
GO TO 301
END IF
END DO
CYCLE
301 CONTINUE
reflow=obsrad(1,knex,icol)
klow=knex
GO TO 220
ELSE
reflow=obsrad(1,knex,icol)
klow=knex
GO TO 180
END IF
END IF
END IF
END DO
END DO
END DO
RETURN
END SUBROUTINE radmcro
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE CLDSET ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE cldset(pr,pt, & 2,3
j3,ptbar,qvbar,rhostr, &
radqvopt,radqcopt,radqropt,radptopt, &
refsat,rhradobs,refcld,cldrad, &
refrain,radsetrat,radreflim,radptgain, &
ref,qc,qr,qi,qsh)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Set RH and microphysical variables at a specific point
! that is above the minimum ceiling and has a reflectivity
! greater than the threshold.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Dingchen Hou, CAPS
! May, 1997
!
! MODIFICATION HISTORY:
!
! 02/02/98 (K. Brewster)
! Addition of refcldopt and old cldrad fixed cloud setting.
! Addition of radsetrat to control amount of cloud added.
! Addition of radreflim to allow limiting the reflectivity
! used to derive cloud amounts.
! Slight streamling of code and addition of documentation.
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! OUTPUT:
!
!-----------------------------------------------------------------------
!
!
IMPLICIT NONE
!
REAL :: pr,pt,j3,ptbar,qvbar,rhostr
INTEGER :: radqvopt,radqcopt,radqropt,radptopt
REAL :: refsat,rhradobs,refcld,cldrad
REAL :: refrain,radsetrat,radreflim,radptgain
REAL :: ref
REAL :: qc,qr,qi,qsh
!
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
!
REAL :: pratio,t,qv,qvsat,reflect
REAL :: ref1,ref2,qc1,qc2
REAL :: dqr,dqc,dpt,dpt1
REAL :: aaa,bbb,pttem
INTEGER :: iter
PARAMETER (ref1=62.96, &
ref2=7.38, &
qc1=2.0E-3, &
qc2=0.1E-3)
!
!-----------------------------------------------------------------------
!
! Function f_qvsat and inline directive for Cray PVP
!
!-----------------------------------------------------------------------
!
REAL :: f_qvsat
!fpp$ expand (f_qvsat)
!dir$ inline always f_qvsat
!*$* inline routine (f_qvsat)
!
!-----------------------------------------------------------------------
!
! Include files
!
!-----------------------------------------------------------------------
!
INCLUDE 'phycst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
dqc=0.
dqr=0.
dpt=0.
!
!
!-----------------------------------------------------------------------
!
! Limit reflectivity according to radreflim
!
!-----------------------------------------------------------------------
!
reflect=AMIN1(ref,radreflim)
IF( radqcopt == 1 .AND. ref > refcld) THEN
dqc=AMAX1((cldrad-(qc+qi)),0.0)
ELSE IF( radqcopt == 2 .AND. ref > refcld) THEN
dqc=(qc1-qc2)/(LOG(ref1)-LOG(ref2))* &
(LOG(reflect)-LOG(ref2))+qc2
dqc=AMAX1(((radsetrat*dqc)-(qc+qi)),0.0)
END IF
!
!-----------------------------------------------------------------------
!
! Adjustment of pt due to qc and qr
!
!-----------------------------------------------------------------------
!
IF( radqropt > 0 .AND. ref > refrain ) THEN
dqr=(10.0**(reflect/10.0)/17300.0)**(4.0/7.0) &
/(rhostr/j3*1000.0)
dqr=AMAX1(((radsetrat*dqr)-(qr+qsh)),0.0)
END IF
!
IF(radptopt == 1 .OR. radptopt == 3) dpt=radptgain*(dqr+dqc)*ptbar/(1.0+qvbar)
pttem=pt+dpt
!
!-----------------------------------------------------------------------
!
! Adjustment of water vapor
!
!-----------------------------------------------------------------------
!
IF(radqvopt > 0 .AND. ref > refsat) THEN
!
!-----------------------------------------------------------------------
!
! Adjustment of pt due to water vapor change
! bbb is the bouyancy multiplied by ptbar/g, to be conserved.
!
!-----------------------------------------------------------------------
!
IF (radptopt == 2 .OR. radptopt == 3) THEN
pratio=(pr/p0) ** rddcp
t =pt*pratio
qvsat=f_qvsat( pr, t )
qv=qvsat*rhradobs
aaa=ptbar*(1-0.622)/((0.622+qvbar)*(1.0+qvbar))
bbb=aaa*qv+pt
DO iter=1,10
t=pttem*pratio
qvsat=f_qvsat( pr, t )
dpt1=bbb-aaa*(qvsat*rhradobs)-pttem
! print *,iter,pttem,t,p,qvsat,aaa,bbb,dpt1
IF (ABS(dpt1) < 1.0E-3) EXIT
pttem=pttem+dpt1
END DO
! 101 CONTINUE
t=pttem*pratio
qvsat=f_qvsat( pr, t )
qv=qvsat*rhradobs
END IF
END IF
qc=qc+dqc
qr=qr+dqr
pt=pttem
RETURN
END SUBROUTINE cldset
!
!##################################################################
!##################################################################
!###### ######
!###### FUNCTION ADASLCL ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!
SUBROUTINE adaslcl(nx,ny,nz,zs,pr,pt,qv, & 1,1
hgtlcl)
IMPLICIT NONE
INTEGER :: nx,ny,nz
REAL :: zs(nx,ny,nz)
REAL :: pr(nx,ny,nz)
REAL :: pt(nx,ny,nz)
REAL :: qv(nx,ny,nz)
REAL :: hgtlcl(nx,ny)
!
!-----------------------------------------------------------------------
!
! Misc local variables
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,k
INTEGER :: imid,jmid
REAL :: pmb,tc,tk,td,wmr,thepcl,plcl,tlcl
REAL :: plnhi,plnlo,plnlcl,whi,wlo
REAL :: oe,wmr2td
!
!-----------------------------------------------------------------------
!
! Include files
!
!-----------------------------------------------------------------------
!
INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
! Find pressure of lcl
! Using pressure of lcl, find height of lcl
!
!-----------------------------------------------------------------------
!
imid=nx/2
jmid=ny/2
DO j=1,ny-1
DO i=1,nx-1
pmb=0.01*pr(i,j,2)
tk=pt(i,j,2)*((1000./pmb)**rddcp)
tc=tk-273.15
wmr=1000.*(qv(i,j,2)/(1.-qv(i,j,2)))
td=wmr2td(pmb,wmr)
thepcl=oe(tc,td,pmb)
CALL ptlcl(pmb,tc,td,plcl,tlcl)
plcl=plcl*100.
DO k=3,nz-2
IF(pr(i,j,k) < plcl) EXIT
END DO
! 81 CONTINUE
plnhi=LOG(pr(i,j,k))
plnlo=LOG(pr(i,j,k-1))
plnlcl=LOG(plcl)
whi=(plnlo-plnlcl)/(plnlo-plnhi)
wlo=1.-whi
hgtlcl(i,j)=whi*zs(i,j,k)+wlo*zs(i,j,k-1)
IF(i == imid .AND. j == jmid) THEN
PRINT *, 'arpslcl: ',pmb,tc,td,wmr,plcl,hgtlcl(i,j)
END IF
END DO
END DO
!
RETURN
END SUBROUTINE adaslcl