!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE RDPROF ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE rdprof(nvar,nzua,mxua,jsrc,proffile, & 1,2
stnua,elevua,xua,yua,hgtua,obsua, &
qualua,isrcua,nlevsua, &
rmiss,nprev,ntotal,istatus)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Read ASCII file containing wind profiler data.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Keith Brewster
! July, 1995
!
! MODIFICATION HISTORY:
! 9/3/96 Keith Brewster
! Added full documentation.
!
! 2/16/98 Keith Brewster
! Added jsrc, source number to the variable list
!
!-----------------------------------------------------------------------
!
! INPUT :
!
! nvar Number of variables in analysis (array dimension)
! nzua Maximum number of vertical levels (array dimension)
! mxua Maximum number of UA stations (array dimension)
! jsrc Source number of data set
! proffile Name of profiler data file to read
! rmiss Missing data fill value
! nprev Number of stations read previously into UA observation
! data arrays (array index)
!
! OUTPUT :
!
! stnua station name (character*4)
! elevua station elevation (m MSL)
! xua station location x-coordinate (m)
! yua station location y-coordinate (m)
! hgtua height of data (m MSL)
! obsua observation data
! qualua observation quality indicator
! isrcua data source index
! nlevsua number of levels of data
! ntotal number of stations
! istatus status indicator
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nvar,nzua,mxua
!
INTEGER :: jsrc
CHARACTER (LEN=256) :: proffile
CHARACTER (LEN=5) :: stnua(mxua)
REAL :: elevua(mxua)
REAL :: xua(mxua)
REAL :: yua(mxua)
REAL :: hgtua(nzua,mxua)
REAL :: obsua(nvar,nzua,mxua)
INTEGER :: qualua(nvar,nzua,mxua)
INTEGER :: isrcua(mxua)
INTEGER :: nlevsua(mxua)
REAL :: rmiss
INTEGER :: nprev,ntotal
INTEGER :: istatus
!
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
!
INTEGER :: ista,ivar,jlev,jend,ksta,mxprof,nprof
INTEGER :: numsta
REAL :: hgtdum,rlat,rlon,DIRECT,speed,ddrot
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
OPEN(12,FILE=trim(proffile),ERR=400,STATUS='old')
!
mxprof=mxua-nprev
!
!-----------------------------------------------------------------------
!
! Fill arrays with "missing data" indicator
!
!-----------------------------------------------------------------------
!
DO ista=1,mxprof
DO jlev=1,nzua
DO ivar=1,nvar
ksta=ista+nprev
obsua(ivar,jlev,ksta)=rmiss
qualua(ivar,jlev,ksta)=-99
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Main data-reading loop
!
!-----------------------------------------------------------------------
!
DO ista=1,mxprof
ksta=ista+nprev
READ(12,'(i12,i12,f11.0,f15.0,f15.0,5x,a5)',ERR=250,END=250) &
numsta,nlevsua(ksta),rlat,rlon,elevua(ksta),stnua(ksta)
PRINT *, 'Reading Wind Profiler <', stnua(ksta), '>'
!
isrcua(ksta)=jsrc
CALL lltoxy
(1,1,rlat,rlon,xua(ksta),yua(ksta))
!
jend=MIN(nlevsua(ksta),nzua)
DO jlev=1,jend
READ(12,*,ERR=250,END=250) hgtua(jlev,ksta),DIRECT,speed
IF(DIRECT >= 0. .AND. speed >= 0.) THEN
CALL ddrotuv(1,rlon,DIRECT,speed,ddrot, &
obsua(1,jlev,ksta),obsua(2,jlev,ksta))
! print *, ' direct,speed,u,v=',direct,speed,
! + obsua(1,jlev,ksta),obsua(2,jlev,ksta)
qualua(1,jlev,ksta)=100
qualua(2,jlev,ksta)=100
END IF
!
END DO
!
!-----------------------------------------------------------------------
!
! Read any extra levels after nzua, but discard them.
!
!-----------------------------------------------------------------------
!
IF(nlevsua(ksta) > nzua) THEN
WRITE(6,'(//a,a/,a,i4,a)') ' RDPROF: WARNING profiler: ', &
stnua(ksta),' truncated to nzua=',nzua,' levels.'
WRITE(6,'(a,i4,a)') ' Data file has ',nlevsua(ksta), &
' levels.'
WRITE(6,'(a/)') ' Increase nz_ua in adas.inc'
DO jlev=jend+1,nlevsua(ksta)
READ(12,*,ERR=250,END=250) hgtdum,DIRECT,speed
END DO
nlevsua(ksta)=nzua
END IF
END DO
!
!-----------------------------------------------------------------------
!
! End-of-file destination
!
!-----------------------------------------------------------------------
!
250 CONTINUE
nprof=ista-1
ntotal=nprev+nprof
WRITE(6,'(a,i4,a)') ' Read ',nprof,' profiler sites'
CLOSE(12)
istatus=1
RETURN
!
!-----------------------------------------------------------------------
!
! Error opening file destination
!
!-----------------------------------------------------------------------
!
400 CONTINUE
WRITE(6,'(a)') ' Error opening profiler file: ',proffile
ntotal=nprev
istatus=-1
RETURN
END SUBROUTINE rdprof
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE RDRAOB ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE rdraob(nvar,nzua,mxua,jsrc,raobfile, & 26,3
stnua,elevua,xua,yua,hgtua,obsua, &
qualua,isrcua,nlevsua, &
rmiss,nprev,ntotal,istatus)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Read ASCII file containing rawinsonde observations (RAOBs).
!
!-----------------------------------------------------------------------
!
! AUTHOR: Keith Brewster
! Jan, 1994
!
! MODIFICATION HISTORY:
! 9/3/96 Keith Brewster
! Restored proper calculation of moisture variable.
! Added full documentation.
!
! 2/16/98 Keith Brewster
! Added jsrc, source number to the variable list
!
!-----------------------------------------------------------------------
!
! INPUT :
!
! nvar Number of variables in analysis (array dimension)
! nzua Maximum number of vertical levels (array dimension)
! mxua Maximum number of UA stations (array dimension)
! jsrc Source number of RAOB data set
! raobfile Name of profiler data file to read
! rmiss Missing data fill value
! nprev Number of stations read previously into UA observation
! data arrays (array index)
!
! OUTPUT :
!
! stnua station name (character*4)
! elevua station elevation (m MSL)
! xua station location x-coordinate (m)
! yua station location y-coordinate (m)
! hgtua height of data (m MSL)
! obsua observation data
! qualua observation quality indicator
! isrcua data source index
! nlevsua number of levels of data
! ntotal number of stations
! istatus status indicator
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nvar,nzua,mxua
!
INTEGER :: jsrc
CHARACTER (LEN=256) :: raobfile
CHARACTER (LEN=5) :: stnua(mxua)
REAL :: elevua(mxua)
REAL :: xua(mxua)
REAL :: yua(mxua)
REAL :: hgtua(nzua,mxua)
REAL :: obsua(nvar,nzua,mxua)
INTEGER :: qualua(nvar,nzua,mxua)
INTEGER :: isrcua(mxua)
INTEGER :: nlevsua(mxua)
REAL :: rmiss
INTEGER :: nprev,ntotal
INTEGER :: istatus
!
REAL :: mbtopa
PARAMETER (mbtopa=100.)
!
INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
!
INTEGER :: ista,ivar,jlev,jend,ksta,nraob,mxraob
INTEGER :: numsta
INTEGER :: nlevlast
REAL :: rlat,rlon,hgtdum,press,temp,tdew,DIRECT,speed,ddrot,tdk
REAL :: u1last,v1last,pr1last,pt1last,rh1last
!
!-----------------------------------------------------------------------
!
! 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...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
OPEN(12,FILE=trim(raobfile),ERR=400,STATUS='old')
!
mxraob=mxua-nprev
nlevlast=0
u1last=-9999.
v1last=-9999.
pr1last=-9999.
pt1last=-9999.
rh1last=-9999.
!
!-----------------------------------------------------------------------
!
! Fill arrays with "missing data" indicator
!
!-----------------------------------------------------------------------
!
DO ista=1,mxraob
DO jlev=1,nzua
DO ivar=1,nvar
ksta=ista+nprev
obsua(ivar,jlev,ksta)=rmiss
qualua(ivar,jlev,ksta)=-99
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Main data-reading loop
!
!-----------------------------------------------------------------------
!
ksta=nprev+1
ntotal=nprev
DO ista=1, mxraob
READ(12,'(i12,i12,f11.4,f15.4,f15.0,5x,a5)',ERR=250,END=250) &
numsta,nlevsua(ksta),rlat,rlon,elevua(ksta),stnua(ksta)
!
isrcua(ksta)=jsrc
CALL lltoxy(1,1,rlat,rlon,xua(ksta),yua(ksta))
!
jend=MIN(nlevsua(ksta),nzua)
DO jlev=1,jend
READ(12,*,ERR=250,END=250) &
hgtua(jlev,ksta),press,temp,tdew,DIRECT,speed
!
!-----------------------------------------------------------------------
!
! observed u,v
!
!-----------------------------------------------------------------------
!
IF(DIRECT >= 0. .AND. speed >= 0.) THEN
CALL ddrotuv(1,rlon,DIRECT,speed,ddrot, &
obsua(1,jlev,ksta),obsua(2,jlev,ksta))
qualua(1,jlev,ksta)=100
qualua(2,jlev,ksta)=100
END IF
!
!-----------------------------------------------------------------------
!
! observed press and potential temperature
!
!-----------------------------------------------------------------------
!
IF(press >= 0.) THEN
obsua(3,jlev,ksta)=press*mbtopa
qualua(3,jlev,ksta)=100
IF(temp > -99.) THEN
obsua(4,jlev,ksta)= &
(temp+273.15)*((p0/obsua(3,jlev,ksta))**rddcp)
qualua(4,jlev,ksta)=100
END IF
!
!-----------------------------------------------------------------------
!
! observed specific humidity
!
!-----------------------------------------------------------------------
!
IF(temp > -99. .AND. tdew > -99.) THEN
tdk=tdew + 273.15
obsua(5,jlev,ksta)= &
MAX(1.0E-08,f_qvsat(obsua(3,jlev,ksta),tdk))
qualua(5,jlev,ksta)=100
END IF
END IF
END DO
!
!-----------------------------------------------------------------------
!
! Read any extra levels after nzua, but discard them.
!
!-----------------------------------------------------------------------
!
IF(nlevsua(ksta) > nzua) THEN
WRITE(6,'(//a,a/,a,i4,a)') ' RDRAOB: WARNING rawinsonde: ', &
stnua(ksta),' truncated to nzua=',nzua,' levels.'
WRITE(6,'(a,i4,a)') ' Data file has ',nlevsua(ksta), &
' levels.'
WRITE(6,'(a/)') ' Increase nz_ua in adas.inc'
!
DO jlev=jend+1,nlevsua(ksta)
READ(12,*,ERR=250,END=250) &
hgtdum,press,temp,tdew,DIRECT,speed
END DO
nlevsua(ksta)=nzua
END IF
!
!-----------------------------------------------------------------------
!
! Check for duplicate sounding to to error some .snd files
! If dupe is found, do not increment ksta or ntotal counter.
!
!-----------------------------------------------------------------------
!
IF(nlevsua(ksta) == nlevlast .AND. obsua(1,1,ksta) == u1last .AND. &
obsua(2,1,ksta) == v1last .AND. &
obsua(3,1,ksta) == pr1last .AND. &
obsua(4,1,ksta) == pt1last .AND. &
obsua(5,1,ksta) == rh1last ) THEN
WRITE(6,'(a,a)') &
' Discarding erroneous duplicate data stored for', &
stnua(ksta)
ELSE
!wdt update
WRITE(6,'(a,i5,2A)') ' Read', nlevsua(ksta), &
' levels of data from radiosonde station ', TRIM(stnua(ksta))
nlevlast=nlevsua(ksta)
u1last=obsua(1,1,ksta)
v1last=obsua(2,1,ksta)
pr1last=obsua(3,1,ksta)
pt1last=obsua(4,1,ksta)
rh1last=obsua(5,1,ksta)
ksta=ksta+1
ntotal=ntotal+1
END IF
END DO
!
!-----------------------------------------------------------------------
!
! End-of-file destination
!
!-----------------------------------------------------------------------
!
250 CONTINUE
nraob=ntotal-nprev
CLOSE(12)
RETURN
!
!-----------------------------------------------------------------------
!
! Error opening file destination
!
!-----------------------------------------------------------------------
!
400 CONTINUE
WRITE(6,'(a)') ' Error opening raob file: ',raobfile
ntotal=nprev
RETURN
END SUBROUTINE rdraob
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE RDGPS ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE rdgps(nx,ny,nz,nvar,nzua,mxua,nsrcua,gpsfile, & 1,6
anx,xs,ys,zp,su,sv,st,spres,shght,sqv,sqvsat, &
srcua,stnua,elevua,xua,yua,hgtua,obsua, &
qualua,isrcua,nlevsua, &
rmiss,nprev,ntotal,istatus)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Read ASCII file containing GPS Integrated Precipitable Water
! data and generate a psuedo-sounding of qv from the background
! qv field at the data location and the IPW information.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Keith Brewster
! June, 2006
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT :
!
! nvar Number of variables in analysis (array dimension)
! nzua Maximum number of vertical levels (array dimension)
! mxua Maximum number of UA stations (array dimension)
! jsrc Source number of RAOB data set
! raobfile Name of profiler data file to read
! rmiss Missing data fill value
! nprev Number of stations read previously into UA observation
! data arrays (array index)
!
! OUTPUT :
!
! stnua station name (character*4)
! elevua station elevation (m MSL)
! xua station location x-coordinate (m)
! yua station location y-coordinate (m)
! hgtua height of data (m MSL)
! obsua observation data
! qualua observation quality indicator
! isrcua data source index
! nlevsua number of levels of data
! ntotal number of stations
! istatus status indicator
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx,ny,nz
INTEGER :: nvar,nzua,mxua,nsrcua
!
CHARACTER (LEN=132) :: gpsfile
CHARACTER (LEN=8) :: srcua(nsrcua)
CHARACTER (LEN=5) :: stnua(mxua)
REAL :: anx(nvar,nx,ny,nz)
REAL :: xs(nx)
REAL :: ys(ny)
REAL :: zp(nx,ny,nz)
REAL :: su(nz)
REAL :: sv(nz)
REAL :: st(nz)
REAL :: spres(nz)
REAL :: shght(nz)
REAL :: sqv(nz)
REAL :: sqvsat(nz)
REAL :: elevua(mxua)
REAL :: xua(mxua)
REAL :: yua(mxua)
REAL :: hgtua(nzua,mxua)
REAL :: obsua(nvar,nzua,mxua)
INTEGER :: qualua(nvar,nzua,mxua)
INTEGER :: isrcua(mxua)
INTEGER :: nlevsua(mxua)
REAL :: rmiss
INTEGER :: nprev,ntotal
INTEGER :: istatus
!
REAL, PARAMETER :: mbtopa=100.
!
INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
! GPS iteration parameters
!
!-----------------------------------------------------------------------
INTEGER, PARAMETER :: maxiter = 20
REAL, PARAMETER :: epsipw = 2.0E-03
REAL, PARAMETER :: delvmax = 200.
REAL, PARAMETER :: stdlapse = -(6.5/1000.)
!
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
!
CHARACTER(LEN=5) :: stngps
CHARACTER(LEN=8) :: srcgps
CHARACTER(LEN=12) :: datetime
INTEGER :: ista,isrc,ivar,ksta,k,kk,mxgps,ngps,nlevs
INTEGER :: ipt,jpt,ireturn,iostatus,niter
REAL :: rlat,rlon,elev,sfcpr,sfct,sfcrh,gpsipw
REAL :: sfcqv,sfcqvsat,sfcqvflg
REAL :: xgps,ygps,selev,temp
LOGICAL :: found
!
!-----------------------------------------------------------------------
!
! 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...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
OPEN(12,FILE=trim(gpsfile),STATUS='old',iostat=iostatus)
IF(iostatus /= 0) THEN
WRITE(6,'(/a,a/)') ' Error opening GPS file: ',gpsfile
ntotal=nprev
RETURN
END IF
READ(12,'(a12)',iostat=iostatus) datetime
IF(iostatus /= 0) THEN
WRITE(6,'(/a,a/)') ' Error reading header in GPS file: ',gpsfile
ntotal=nprev
CLOSE(12)
RETURN
END IF
WRITE(6,'(//a,a)') ' Reading GPS data for date-time: ',datetime
!
mxgps=mxua-nprev
!
!-----------------------------------------------------------------------
!
! Main data-reading loop
!
!-----------------------------------------------------------------------
!
ksta=nprev
ntotal=nprev
DO ista=1, mxgps
READ(12,'(a5,1x,a8,9x,f10.5,f11.5,f8.2,f8.2,f6.2,f8.2,f7.1)', &
iostat=iostatus) &
stngps,srcgps,rlat,rlon,elev,sfcpr,gpsipw,sfct,sfcrh
print *, ' '
WRITE(6,'(/a,a5,1x,a8,1x,f10.5,f11.5,f8.2)') &
' stngps,srcgps,rlat,rlon,elev: ',stngps,srcgps,rlat,rlon,elev
WRITE(6,'(a,4f8.2)') &
' sfcpr,gpsipw,sfct,sfcrh: ',sfcpr,gpsipw,sfct,sfcrh
IF(iostatus /= 0) EXIT
IF(gpsipw < 0. .OR. sfcpr < 0. ) CYCLE
sfcpr=mbtopa*sfcpr
!
CALL lltoxy(1,1,rlat,rlon,xgps,ygps)
!
!-----------------------------------------------------------------------
!
! Find location in ARPS grid.
!
!-----------------------------------------------------------------------
!
CALL findlc(nx,ny,xs,ys,xgps,ygps,ipt,jpt,ireturn)
WRITE(6,'(a,2f12.1,2i8)') ' xgps,ygps,ipt,jpt:',xgps,ygps,ipt,jpt
!
!-----------------------------------------------------------------------
!
! Interpolate if extrapoltion is not indicated.
!
!-----------------------------------------------------------------------
!
IF(ireturn == 0) THEN
found=.false.
DO isrc=1,nsrcua
IF(srcua(isrc) == srcgps) THEN
found=.true.
EXIT
END IF
END DO
IF(.NOT. found) THEN
WRITE(6,'(a,a,a)') ' Could not find GPS source ',srcgps, &
' among upper air sources in input file'
STOP
END IF
write(6,'(a,a)') stngps, &
' inside ADAS domain, creating pseudo sounding'
CALL colint(nx,ny,nz,nvar, &
xs,ys,zp,xgps,ygps,ipt,jpt,anx, &
su,sv,st,spres,shght,sqv,selev, &
nlevs)
WRITE(6,'(a,f6.1,a,f6.1)') &
' GPS Elev: ',elev,' Model elev: ',selev
IF(abs(elev-selev) > delvmax) CYCLE
!
! Convert theta to temperature and compute qv_sat
!
DO k=1,nlevs
temp=st(k)*((spres(k)/p0)**rddcp)
sqvsat(k)=f_qvsat(spres(k),temp)
st(k)=temp
END DO
IF(sfct > -50.0 ) THEN
sfct=sfct+273.15
ELSE
sfct=st(1)+stdlapse*(elev-shght(1))
END IF
sfcqvsat=f_qvsat(sfcpr,sfct)
IF(sfcrh > 0.) THEN
sfcqv=(0.01*sfcrh)*sfcqvsat
sfcqvflg=1
ELSE
sfcqv=sqv(1)
sfcqvflg=0
END IF
CALL ipw2snd(nlevs,shght,spres,st,sqv,sqvsat, &
elev,sfcpr,sfct,sfcqv,sfcqvsat,sfcqvflg,gpsipw, &
maxiter,epsipw,niter,istatus)
IF(istatus == 0) THEN
ksta=ksta+1
ntotal=ntotal+1
DO k=1,nlevs
DO ivar=1,nvar
obsua(ivar,k,ksta)=rmiss
qualua(ivar,k,ksta)=-13
END DO
END DO
xua(ksta)=xgps
yua(ksta)=ygps
elevua(ksta)=elev
stnua(ksta)=stngps
isrcua(ksta)=isrc
hgtua(1,ksta)=elev
obsua(5,1,ksta)=sfcqv
qualua(5,1,ksta)=100
kk=1
DO k=1,nlevs
IF(spres(k) < sfcpr .AND. shght(k) > elev) THEN
kk=kk+1
hgtua(kk,ksta)=shght(k)
obsua(5,kk,ksta)=sqv(k)
qualua(5,kk,ksta)=100
END IF
END DO
nlevsua(ksta)=kk
END IF ! success in IPW adjustment
END IF ! location in grid
END DO
!
!-----------------------------------------------------------------------
!
! End-of-file destination
!
!-----------------------------------------------------------------------
!
ngps=ntotal-nprev
WRITE(6,'(/a,i6,a/)')' Read and processed ',ngps,' GPS stations'
CLOSE(12)
RETURN
!
!-----------------------------------------------------------------------
!
! Error opening file destination
!
!-----------------------------------------------------------------------
!
RETURN
END SUBROUTINE rdgps
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE IPW2SND ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE ipw2snd(nzsnd,htsnd,psnd,tsnd,qvsnd,qvsat, & 1,1
elev,sfcpr,sfct,sfcqv,sfcqvsat,sfcqvflg,obsipw, &
maxiter,epsipw,niter,istatus)
!
! Given a background sounding, a measurement of integrated
! precipitable water (IPW), surface pressure, temperature
! and specific humidity generate a pseudo-sounding of
! specific humidity.
!
IMPLICIT NONE
INTEGER, INTENT(IN) :: nzsnd
REAL, INTENT(IN) :: htsnd(nzsnd)
REAL, INTENT(IN) :: psnd(nzsnd)
REAL, INTENT(IN) :: tsnd(nzsnd)
REAL, INTENT(INOUT) :: qvsnd(nzsnd)
REAL, INTENT(INOUT) :: qvsat(nzsnd)
REAL, INTENT(IN) :: elev
REAL, INTENT(IN) :: sfcpr
REAL, INTENT(IN) :: sfct
REAL, INTENT(INOUT) :: sfcqv
REAL, INTENT(IN) :: sfcqvsat
INTEGER, INTENT(IN) :: sfcqvflg
REAL, INTENT(IN) :: obsipw
INTEGER, INTENT(IN) :: maxiter
REAL, INTENT(IN) :: epsipw
INTEGER, INTENT(OUT) :: niter
INTEGER, INTENT(OUT) :: istatus
!
! Misc internal variables
!
REAL, PARAMETER :: g=9.8
REAL, PARAMETER :: rhmin=0.05
REAL, PARAMETER :: rhmax=1.0
REAL :: pwcst,sfcipw,sndipw,qvbar,dipw,pwratio,qvnew
REAL :: f_qvsat
INTEGER :: k,kstart,iter
!
! Initializations
!
istatus=0
pwcst=0.1/g
DO k=1,nzsnd-1
IF(psnd(k) > 0.) THEN
qvsat(k)=f_qvsat(psnd(k),tsnd(k))
ELSE
qvsat(k)=1.0E-06
END IF
END DO
DO k=1,nzsnd-1
IF(psnd(k) < sfcpr .AND. htsnd(k) > elev) THEN
kstart=k
EXIT
END IF
END DO
!
! Iterate until integrated precipitable water difference
! falls below the threshold, epsilon ipw
!
IF(sfcqvflg == 1) THEN ! observed surface qv
DO iter=1,maxiter
qvbar=0.5*(sfcqv+qvsnd(kstart))
sndipw=pwcst*qvbar*(sfcpr-psnd(kstart))
sfcipw=0.5*sndipw
DO k=kstart,nzsnd-1
IF(psnd(k) > 0. .AND. psnd(k+1) > 0.) THEN
qvbar=0.5*(qvsnd(k)+qvsnd(k+1))
sndipw=sndipw+(pwcst*qvbar*(psnd(k)-psnd(k+1)))
END IF
END DO
dipw=abs(sndipw-obsipw)
pwratio=(obsipw-sfcipw)/(sndipw-sfcipw)
WRITE(6,'(a,i4,4(a,f9.2))') &
' Iter:',iter,' IPW:',obsipw,' Snd IPW:',sndipw, &
' Diff:',dipw,' Ratio:',pwratio
IF( dipw < epsipw) EXIT
DO k=1,nzsnd
IF(psnd(k) > 0.) THEN
qvnew=pwratio*qvsnd(k)
qvsnd(k)=min(max(qvnew,(rhmin*qvsat(k))),(rhmax*qvsat(k)))
END IF
END DO
END DO
ELSE ! sfcqv is estimated
DO iter=1,maxiter
qvbar=0.5*(sfcqv+qvsnd(kstart))
sndipw=pwcst*qvbar*(sfcpr-psnd(kstart))
DO k=kstart,nzsnd-1
IF(psnd(k) > 0. .AND. psnd(k+1) > 0.) THEN
qvbar=0.5*(qvsnd(k)+qvsnd(k+1))
sndipw=sndipw+(pwcst*qvbar*(psnd(k)-psnd(k+1)))
END IF
END DO
dipw=abs(sndipw-obsipw)
pwratio=obsipw/sndipw
WRITE(6,'(a,i4,4(a,f9.2))') &
' Iter:',iter,' IPW:',obsipw,' Snd IPW:',sndipw, &
' Diff:',dipw,' Ratio:',pwratio
IF( dipw < epsipw) EXIT
qvnew=pwratio*sfcqv
sfcqv=min(max(qvnew,(rhmin*sfcqvsat)),(rhmax*sfcqvsat))
DO k=1,nzsnd
IF(psnd(k) > 0.) THEN
qvnew=pwratio*qvsnd(k)
qvsnd(k)=min(max(qvnew,(rhmin*qvsat(k))),(rhmax*qvsat(k)))
END IF
END DO
END DO
END IF
IF(iter > maxiter) THEN
niter=maxiter
istatus=-1
ELSE
niter=iter
istatus=0
END IF
RETURN
END SUBROUTINE IPW2SND