!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE PREPSGLOBS ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE prepsglobs(mxsng,ntime,srcsng,nsrc_sng,sngfile, & 1,6
stnsng,latsng,lonsng,xsng,ysng,hgtsng, &
obsng,qobsng,qualsng,isrcsng,qsrc,chsrc, &
nx,ny,nz,nvar,anx,xs,ys,zp, &
shght,su,sv,spres,stheta,sqv, &
rmiss,nobsng,istatus)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! This subroutine reads non-surface single-level data.
! The data are appended to surface data already collected.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Keith Brewster, CAPS
! October, 1997
!
! MODIFICATION HISTORY:
!
! 5/4/98 Developed based on RDRASS to assign height from background
! Can also be used to read satellite winds (John Manobianco).
!
! 1998/07/13 (Gene Bassett) Implemented general format for sgl obs
! and used this for MDCRS data.
!
! 1999/09/16 (Dingchen Hou) Modified the interpretion of MDCRS data.
!
! 1999/11/19 (Gene Bassett, Keith Brewster)
! Modified MDCRS data to use pressure only for deriving other
! observational quantities but not as a pressure observation in
! ADAS proper.
!
! 11/9/2001 (Keith Brewster)
! Added FORTRAN-90 features and rearranged the logic to eliminate
! GO-TOs, added processing for humidity and pressure data,
! and generally make it more suitable for additional data
! sources other than aircraft data.
!
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nvar,mxsng,ntime
!
INTEGER :: jsrc
INTEGER :: nsrc_sng
CHARACTER (LEN=132) :: sngfile, line
CHARACTER (LEN=8) :: srcsng(nsrc_sng)
CHARACTER (LEN=8) :: chsrc(mxsng,ntime)
CHARACTER (LEN=5) :: stnsng(mxsng)
REAL :: latsng(mxsng)
REAL :: lonsng(mxsng)
REAL :: xsng(mxsng)
REAL :: ysng(mxsng)
REAL :: hgtsng(mxsng)
REAL :: obsng(nvar,mxsng)
REAL :: qobsng(nvar,mxsng)
INTEGER :: qualsng(nvar,mxsng)
REAL :: qsrc(nvar,mxsng)
INTEGER :: isrcsng(mxsng)
REAL :: rmiss
INTEGER :: nprev
INTEGER :: nobsng
INTEGER :: istatus
INTEGER :: nx,ny,nz
REAL :: anx(nx,ny,nz,nvar)
REAL :: xs(nx),ys(ny)
REAL :: zp(nx,ny,nz)
REAL :: shght(nz)
REAL :: su(nz)
REAL :: sv(nz)
REAL :: spres(nz)
REAL :: stheta(nz)
REAL :: sqv(nz)
REAL :: ktstoms,mbtopa
PARAMETER (ktstoms=0.5144,mbtopa=100.)
!
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
!
INTEGER :: ivar,ksta,nsingle
INTEGER :: ipt,jpt,kgrid
INTEGER :: indom
INTEGER :: nlevs
CHARACTER (LEN=8) :: obs_id
INTEGER :: obs_num
CHARACTER (LEN=8) :: date_str
CHARACTER (LEN=6) :: time_str
INTEGER :: num_lines
CHARACTER (LEN=8) :: obs_type
REAL :: obs_lat,obs_lon,hgt_msl,hgt_agl
REAL :: obs_t,obs_td,obs_wspd,obs_wdir,obs_gust
REAL :: obs_psta,obs_pmsl,obs_palt,obs_vis
REAL :: xsta,ysta,tk,tdk,psta,qv,qvsat
CHARACTER (LEN=16) :: obs_wx
CHARACTER (LEN=6) :: fmt_in
INTEGER :: k,iline,istat
REAL :: weight
REAL :: ddrot,selev
REAL :: altpres ! presure derived from altimeter reading
!
!-----------------------------------------------------------------------
!
! Functions
!
!-----------------------------------------------------------------------
!
REAL :: alttostpr
REAL :: msltostpr
REAL :: f_qvsat
REAL :: ztopsa
!fpp$ expand (alttostpr)
!dir$ inline always alttostpr
!*$* inline routine (alttostpr)
!fpp$ expand (msltostpr)
!dir$ inline always msltostpr
!*$* inline routine (msltostpr)
!fpp$ expand (f_qvsat)
!dir$ inline always f_qvsat
!*$* inline routine (f_qvsat)
!fpp$ expand (ztopsa)
!dir$ inline always ztopsa
!*$* inline routine (ztopsa)
!
!-----------------------------------------------------------------------
!
! Include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'phycst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
nprev=nobsng
WRITE(6,'(a,a)') ' Opening: ',sngfile
WRITE(6,'(a,i6,i6)') ' In PREPSGLOBS , mxsng, nprev = ', &
mxsng,nprev
OPEN(12,FILE=trim(sngfile),STATUS='old',iostat=istat)
IF( istat /= 0 ) THEN
WRITE(6,'(a)') ' Error opening file: ',sngfile
istatus=-1
RETURN
END IF
!
!-----------------------------------------------------------------------
!
! Main data-reading loop
!
!-----------------------------------------------------------------------
!
ksta = nobsng
DO iline=1,999999
!
!-----------------------------------------------------------------------
!
! INPUT FORMAT:
!
! first line: obs_id, obs_num, date_str, time_str,
! obs_lat, obs_lon, hgt_msl, hgt_agl,
! obs_type, num_lines
!
! See MDCRS format below.
!
! obs_id station id (8 characters)
! obs_num station id number
! date_str yyyymmdd
! time_str hhmm
! obs_lat latidute (degrees)
! obs_lon longitude (degrees)
! hgt_msl height above mean sea level (m)
! value considered missing if <= rmiss
! hgt_agl height above ground level (m), optional,
! intended for use with surface obs to help with
! mis-matches in terrain
! value considered missing if <= rmiss
! obs_type ADAS observation type (8 characters)
! num_lines number of observation lines to follow
! the header line
!
! currently unused in ADAS:
! obs_id, obs_num, date_str, time_str, hgt_agl
!
!
! currently supported observation lines formatted as follows:
!
! FMT1 20.0 -10.0 215 10.0 100.0 1013.2 1013.2 100. +TSRA
! ______ _____ _____ ___ _____ _____ ______ ______ ____ ________________
!
! FMT1 obs_t, obs_td, obs_wdir, obs_wspd, obs_gust,
! obs_pmsl, obs_palt, obs_vis, obs_wx
!
! obs_t temperature (C)
! value considered missing if <= rmiss
! obs_td dew point (C)
! value considered missing if <= rmiss
! obs_wdir wind direction (degrees clockwise from north)
! value considered missing if < 0
! obs_wspd wind speed (m/s)
! value considered missing if <= rmiss
! obs_wdir wind direction (degrees clockwise from north)
! value considered missing if < 0
! obs_gust wind gust (m/s)
! value considered missing if <= rmiss
! obs_pmsl sea level pressure (mb)
! value considered missing if <= rmiss
! obs_palt altimiter setting (mb)
! value considered missing if <= rmiss
! obs_vis visibility (km)
! value considered missing if < 0
! obs_wx current weather (METAR format)
!
! currently unused in ADAS:
! obs_gust, obs_palt, obs_vis, obs_wx
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
! Initially set all non-header line values to missing
!
!-----------------------------------------------------------------------
!
obs_t = rmiss
obs_td = rmiss
obs_wdir = -1
obs_wspd = rmiss
obs_wdir = -1
obs_gust = rmiss
obs_psta = rmiss
obs_pmsl = rmiss
obs_palt = rmiss
obs_vis = rmiss
obs_wx = ' '
altpres = rmiss
psta = rmiss
!
!-----------------------------------------------------------------------
!
! Read in the info for one station
! Skip comments
! Comments begin with # or !. Also, skip lines that begin with a space.
!
! Sample format
!XXXXXXXX 00000000 19990601_153600 42.7167 -84.9517 9448 -999.00 MDCRS 1
!-------- -------- -------- ------ -------- --------- ------ ------- -------- --
! A8, 1X, I8,1X,I4,2I2,1X,2I2,2X,1X,F8.4,1X,F9.4, 1X, I6,1X,F7.2,1X,A8,1X,I2
! MDCRS -39.8 -999.0 267 28.3 YYZ ORD FSL00307
!-----------------------------------------------------------------------
!
READ (12,'(A)',iostat=istat) line
IF (istat /= 0) EXIT
IF (line(1:1) == '#' .OR. line(1:1) == '!' .OR. line(1:1) == ' ') CYCLE
READ(line, &
'(a8,1x,i8,1x,a8,1x,a6,1x,f8.4,1x,f9.4,1x,f6.0,1x,f7.2,1x,a8,1x,i2)', &
iostat=istat) obs_id, obs_num, date_str, time_str, &
obs_lat, obs_lon, hgt_msl, hgt_agl, obs_type, num_lines
IF (istat /= 0) CYCLE
READ(12, &
'(2X,a6,1X,f6.1,1X,f6.1,1X,f3.0,1X,f5.1)',iostat=istat) &
fmt_in, obs_t, obs_td, obs_wdir, obs_wspd
IF (istat /= 0) CYCLE
! WRITE(6,'(2F9.2, F7.0, 2F7.1, F5.0, F5.1)') &
! obs_lat, obs_lon, hgt_msl, obs_t, obs_td, obs_wdir, obs_wspd
CALL lltoxy
(1,1,obs_lat,obs_lon,xsta,ysta)
CALL findlc(nx,ny,xs,ys,xsta,ysta,ipt,jpt,indom)
IF(indom == 0) THEN
!
!-----------------------------------------------------------------------
!
! Station pressure (mb)
! Pressure altitude is assumed here.
!
!-----------------------------------------------------------------------
!
IF ( obs_type == 'MDCRS ' ) THEN
!
!-----------------------------------------------------------------------
!
! Need to convert from pressure altitude to pressure and set
! true height based on the background pressure field.
!
!-----------------------------------------------------------------------
!
altpres = mbtopa*ztopsa(hgt_msl)
!
!-----------------------------------------------------------------------
!
! Interpolate background field (in the horizontal)
! for the whole vertical column.
!
!-----------------------------------------------------------------------
!
CALL colint(nx,ny,nz,nvar, &
xs,ys,zp,xsta,ysta,ipt,jpt,anx, &
su,sv,stheta,spres,shght,sqv,selev, &
nlevs)
!
!-----------------------------------------------------------------------
!
! Interpolate hgts in the vertical from background.
! Interpolation is linear in log-p.
!
!-----------------------------------------------------------------------
!
kgrid = 0
hgt_msl = rmiss
DO k=3,nlevs-1
IF(spres(k) < altpres) THEN
kgrid = k
EXIT
END IF
END DO
IF (kgrid > 0) THEN
weight = (LOG(altpres)-LOG(spres(kgrid)))/ &
(LOG(spres(kgrid-1))-LOG(spres(kgrid)))
hgt_msl = (1.-weight)*shght(kgrid) &
+ weight*shght(kgrid-1)
END IF
END IF
!
!-----------------------------------------------------------------------
!
! Use this datapoint if we have a valid height for it.
!
!-----------------------------------------------------------------------
!
IF (chsrc(ksta,1) == 'MDCRS ') THEN
ksta = ksta + 1
IF(ksta > mxsng) THEN
WRITE(6,'(a,i6,a)') &
' ERROR: number of single level stations,',mxsng, &
'exceeded in PREPSGLOBS. Increase mxsng and retry.'
istatus = -2
STOP
END IF
DO ivar=1,nvar
obsng(ivar,ksta)=rmiss
qobsng(ivar,ksta)=999.
qualsng(ivar,ksta)=-99
END DO
stnsng(ksta) = obs_id(1:5)
latsng(ksta) = obs_lat
lonsng(ksta) = obs_lon
xsng(ksta) = xsta
ysng(ksta) = ysta
chsrc(ksta,1) = obs_type
DO jsrc=1,nsrc_sng
IF (chsrc(ksta,1) == srcsng(jsrc)) THEN
isrcsng(ksta)=jsrc
EXIT
END IF
END DO
jsrc=isrcsng(ksta)
!
!-----------------------------------------------------------------------
!
! Set obs height
!
!-----------------------------------------------------------------------
!
hgtsng(ksta) = hgt_msl
!
!-----------------------------------------------------------------------
!
! Station pressure (mb)
! For now we'll convert msl pressure to station pressure
! using the current temperature.
!
!-----------------------------------------------------------------------
!
IF(obs_psta > rmiss ) THEN
psta=mbtopa*obs_psta
obsng(3,ksta)=psta
qobsng(3,ksta)=qsrc(3,jsrc)
qualsng(3,ksta)=100
ELSE IF(obs_palt > rmiss ) THEN
psta=mbtopa*alttostpr(obs_palt,hgt_msl)
obsng(3,ksta)=psta
qobsng(3,ksta)=qsrc(3,jsrc)
qualsng(3,ksta)=100
ELSE IF(obs_pmsl > rmiss .AND. &
obs_t > rmiss .AND. &
hgt_msl < 500. ) THEN
tk=obs_t + 273.15
psta=mbtopa*msltostpr(obs_pmsl,tk,hgt_msl)
obsng(3,ksta)=psta
qobsng(3,ksta)=qsrc(3,jsrc)
qualsng(3,ksta)=100
ELSE IF (altpres > rmiss) THEN
psta=altpres
ELSE
psta=mbtopa*alttostpr(1013.,hgt_msl)
END IF
!
!-----------------------------------------------------------------------
!
! Set winds.
!
!-----------------------------------------------------------------------
!
IF(obs_wdir >= 0. .AND. obs_wspd >= 0.) THEN
CALL ddrotuv(1,lonsng(ksta),obs_wdir,obs_wspd,ddrot, &
obsng(1,ksta),obsng(2,ksta))
qualsng(1,ksta)=100
qualsng(2,ksta)=100
qobsng(1,ksta)=qsrc(1,jsrc)
qobsng(2,ksta)=qsrc(2,jsrc)
END IF
!
!-----------------------------------------------------------------------
!
! Potential temperature (K)
!
!-----------------------------------------------------------------------
!
IF(obs_t > rmiss .AND. psta > rmiss) THEN
tk=obs_t + 273.15
obsng(4,ksta)=tk*((p0/psta)**rddcp)
qualsng(4,ksta)=100
qobsng(4,ksta)=qsrc(4,jsrc)
END IF
!-----------------------------------------------------------------------
!
! Specific humidity (kg/kg)
!
!-----------------------------------------------------------------------
!
IF(obs_td > rmiss .AND. obs_t > rmiss .AND. &
psta > rmiss ) THEN
tk=obs_t + 273.15
tdk=obs_td + 273.15
qv = f_qvsat( psta,tdk )
qvsat = f_qvsat( psta, tk )
obsng(5,ksta)=max((0.05*qvsat),min(qv,qvsat))
qobsng(5,ksta)=qsrc(5,jsrc)
qualsng(5,ksta)=100
END IF
!
END IF
END IF
END DO
nobsng = ksta
nsingle = nobsng - nprev
WRITE(6,'(a,i6,i6,i6)') ' PREPSGLOBS: nsingle,nprev,nobsng', &
nsingle,nprev,nobsng
istatus=0
CLOSE(12)
RETURN
END SUBROUTINE prepsglobs