!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE GETARPS ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE getarps(nx_ext,ny_ext,nz_ext,nzsoil_ext, & 2,21
dir_extd,extdname,extdopt,extdfmt, &
extdinit,extdfcst,julfname,nstyps, &
iproj_ext,scale_ext, &
trlon_ext,latnot_ext,x0_ext,y0_ext, &
lat_ext,lon_ext,latu_ext,lonu_ext,latv_ext,lonv_ext, &
p_ext,hgt_ext,zp_ext,zpsoil_ext,t_ext,qv_ext, &
u_ext,vatu_ext,v_ext,uatv_ext,w_ext, &
qc_ext,qr_ext,qi_ext,qs_ext,qh_ext, &
soiltyp_ext,stypfrct_ext,vegtyp_ext, &
lai_ext,roufns_ext,veg_ext, &
tsoil_ext,qsoil_ext,wetcanp_ext, &
tem1_ext,tem2_ext,istatus)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! ARPS version.
!
! Reads an ARPS file for processing by ext2arps, a program
! that converts external files to ARPS variables and format.
! This version is useful when you want to use an ARPS file
! with a different orientation or terrain than your final
! ARPS product so arpsr2h does not work.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Keith Brewster
! November, 1995
!
! MODIFICATION HISTORY:
! 01/16/1996 (Yuhe Liu)
! Added three arguments to specify the "external" ARPS data file
! names and format.
!
! 2000/08/14 (Gene Bassett)
! Added multiple soil types, sfcdata variables and grid staggering
! for use with arps history format of external model data.
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! dir_extd Directory name for external file
! extdname Prefix string of external file name
! extdopt Option of external data sources
! extdfmt Option of external data format
!
! extdinit Initialized time in mm-dd-yyyy:hh:mm:ss format
! extdfcst Forecast hour in HHH:MM:SS format
! julfname File name in yyjjjhhmm format
!
! OUTPUT:
!
! iproj_ext Map projection number of external data
! scale_ext Scale factor of external data
! trlon_ext True longitude of external data (degrees E)
! latnot_ext(2) True latitude(s) of external data (degrees N)
! x0_ext x coordinate of origin of external data
! y0_ext y coordinate of origin of external data
! lat_ext latitude of external data points (degrees N)
! lon_ext longitude of external data points (degrees E)
! latu_ext latitude of external u points (degrees N)
! lonu_ext longitude of external u points (degrees E)
! latv_ext latitude of external v points (degrees N)
! lonv_ext longitude of external v points (degrees E)
! p_ext pressure (Pascal)
! hgt_ext height (m)
! zp_ext height (m) (on arps grid)
! zpsoil_ext height (m) (soil model on arps grid)
! t_ext temperature (K)
! qv_ext specific humidity (kg/kg)
! u_ext u wind component (m/s) (staggered grid, true n-s component)
! vatu_ext v wind component at u location (m/s)
! v_ext v wind component (m/s) (staggered grid, true e-w component)
! uatv_ext u wind component at v location (m/s)
! qc_ext Cloud water mixing ratio (kg/kg)
! qr_ext Rain water mixing ratio (kg/kg)
! qi_ext Ice mixing ratio (kg/kg)
! qs_ext Snow mixing ratio (kg/kg)
! qh_ext Hail mixing ratio (kg/kg)
!
! soiltyp_ext Soil type
! stypfrct_ext Soil type fraction
! vegtyp_ext Vegetation type
! lai_ext Leaf Area Index
! roufns_ext Surface roughness
! veg_ext Vegetation fraction
!
! tsoil_ext Soil temperature (K)
! qsoil_ext Soil moisture
! wetcanp_ext Water content on canopy
!
! istatus status indicator
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
CHARACTER (LEN=*) :: dir_extd
CHARACTER (LEN=*) :: extdname
INTEGER :: extdopt
INTEGER :: extdfmt
INTEGER :: nstyps
CHARACTER (LEN=19) :: extdinit
CHARACTER (LEN=9) :: extdfcst
CHARACTER (LEN=9) :: julfname
!
!-----------------------------------------------------------------------
!
! External grid variables
!
!-----------------------------------------------------------------------
!
INTEGER :: iproj_ext
REAL :: scale_ext,trlon_ext
REAL :: latnot_ext(2)
REAL :: x0_ext,y0_ext
!
!-----------------------------------------------------------------------
!
! Output external variable arrays
!
!-----------------------------------------------------------------------
!
INTEGER :: nx_ext,ny_ext,nz_ext,nzsoil_ext
REAL :: lat_ext(nx_ext,ny_ext)
REAL :: lon_ext(nx_ext,ny_ext)
REAL :: latu_ext(nx_ext,ny_ext)
REAL :: lonu_ext(nx_ext,ny_ext)
REAL :: latv_ext(nx_ext,ny_ext)
REAL :: lonv_ext(nx_ext,ny_ext)
REAL :: p_ext(nx_ext,ny_ext,nz_ext) ! Pressure (Pascals)
REAL :: hgt_ext(nx_ext,ny_ext,nz_ext) ! Height (m)
REAL :: zp_ext(nx_ext,ny_ext,nz_ext) ! Height (m) (on arps grid)
REAL :: zpsoil_ext(nx_ext,ny_ext,nz_ext) ! Height (m) (soil model)
REAL :: t_ext(nx_ext,ny_ext,nz_ext) ! Temperature (K)
REAL :: qv_ext(nx_ext,ny_ext,nz_ext) ! Specific humidity (kg/kg)
REAL :: u_ext(nx_ext,ny_ext,nz_ext) ! Eastward wind component (m/s)
REAL :: v_ext(nx_ext,ny_ext,nz_ext) ! Northward wind component (m/s)
REAL :: uatv_ext(nx_ext,ny_ext,nz_ext)
REAL :: vatu_ext(nx_ext,ny_ext,nz_ext)
REAL :: w_ext(nx_ext,ny_ext,nz_ext) ! Vertical velocity component (m/s)
REAL :: qc_ext(nx_ext,ny_ext,nz_ext) ! Cloud H2O mixing ratio (kg/kg)
REAL :: qr_ext(nx_ext,ny_ext,nz_ext) ! Rain H2O mixing ratio (kg/kg)
REAL :: qi_ext(nx_ext,ny_ext,nz_ext) ! Ice mixing ratio (kg/kg)
REAL :: qs_ext(nx_ext,ny_ext,nz_ext) ! Snow mixing ratio (kg/kg)
REAL :: qh_ext(nx_ext,ny_ext,nz_ext) ! Hail mixing ratio (kg/kg)
INTEGER soiltyp_ext (nx_ext,ny_ext,nstyps) ! Soil type
REAL stypfrct_ext(nx_ext,ny_ext,nstyps) ! Soil type
INTEGER vegtyp_ext (nx_ext,ny_ext) ! Vegetation type
REAL lai_ext (nx_ext,ny_ext) ! Leaf Area Index
REAL roufns_ext (nx_ext,ny_ext) ! Surface roughness
REAL veg_ext (nx_ext,ny_ext) ! Vegetation fraction
REAL :: tsoil_ext(nx_ext,ny_ext,nzsoil_ext,0:nstyps) ! Soil temperature(K)
REAL :: qsoil_ext(nx_ext,ny_ext,nzsoil_ext,0:nstyps) ! Soil moisture (m3/m3)
REAL :: wetcanp_ext(nx_ext,ny_ext,0:nstyps) ! Canopy water amount
REAL :: snowdpth_ext(nx_ext,ny_ext) ! Snow depth (m)
REAL :: tem1_ext(nx_ext,ny_ext,nz_ext) ! Temporary work array
REAL :: tem2_ext(nx_ext,ny_ext,nz_ext) ! Temporary work array
!
!-----------------------------------------------------------------------
!
! Other external variable arrays
!
!-----------------------------------------------------------------------
!
REAL :: x_ext(nx_ext)
REAL :: y_ext(ny_ext)
REAL :: z_ext(nz_ext)
REAL :: xsc(nx_ext)
REAL :: ysc(ny_ext)
!
! REAL :: uprt_ext(nx_ext,ny_ext,nz_ext) ! x velocity component (m/s)
! REAL :: vprt_ext(nx_ext,ny_ext,nz_ext) ! y velocity component (m/s)
REAL :: ubar_ext(nx_ext,ny_ext,nz_ext) ! Base state x velocity
! component (m/s)
REAL :: vbar_ext(nx_ext,ny_ext,nz_ext) ! Base state y velocity
! component (m/s)
REAL :: wbar_ext(nx_ext,ny_ext,nz_ext) ! Base state z velocity
! component (m/s)
REAL :: ptbar_ext(nx_ext,ny_ext,nz_ext) ! Base state potential
! temperature (K)
REAL :: pbar_ext(nx_ext,ny_ext,nz_ext) ! Base state pressure (Pascal)
REAL :: qvbar_ext(nx_ext,ny_ext,nz_ext) ! Base state water vapor
! mixing ratio (kg/kg)
! real raing_ext (nx_ext,ny_ext) ! Grid supersaturation rain
! real rainc_ext (nx_ext,ny_ext) ! Cumulus convective rain
! real prcrate_ext(nx_ext,ny_ext,4)
! precipitation rate (kg/(m**2*s))
! prcrate(1,1,1) = total precip. rate
! prcrate(1,1,2) = grid scale precip. rate
! prcrate(1,1,3) = cumulus precip. rate
! prcrate(1,1,4) = microphysics precip. rate
! real radfrc_ext(nx,ny,nz) ! Radiation forcing (K/s)
! real radsw_ext (nx,ny) ! Solar radiation reaching the surface
! real rnflx_ext (nx,ny) ! Net radiation flux absorbed by surface
! real usflx_ext (nx,ny) ! Surface flux of u-momentum (kg/(m*s**2))
! real vsflx_ext (nx,ny) ! Surface flux of v-momentum (kg/(m*s**2))
! real ptsflx_ext(nx,ny) ! Surface heat flux (K*kg/(m*s**2))
! real qvsflx_ext(nx,ny) ! Surface moisture flux (kg/(m**2*s))
INTEGER :: istatus
!
!-----------------------------------------------------------------------
!
! Misc internal variables
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=3) :: fmtn
CHARACTER (LEN=80) :: timsnd
INTEGER :: lenrun, tmstrln
INTEGER :: i,j,k,ldir,ireturn
INTEGER :: ihr,imin,isec
REAL :: govrd
REAL :: xctr,yctr,dx_ext,dy_ext,tvbot,tvtop,tvbar
CHARACTER (LEN=80) :: runname_org
CHARACTER (LEN=80) :: cmnt_org(50)
INTEGER :: nocmnt_org,mapproj_org
INTEGER :: month_org,day_org,year_org
INTEGER :: hour_org,minute_org,second_org
REAL :: latnot(2)
REAL :: umove_org,vmove_org,xgrdorg_org,ygrdorg_org
REAL :: trulat1_org,trulat2_org,trulon_org,sclfct_org
REAL :: latitud_org,ctrlat_org,ctrlon_org
REAL :: dx_org,dy_org
REAL :: lat_org,lon_org
REAL :: dz_org,dzmin_org,zrefsfc_org,dlayer1_org,dlayer2_org
REAL :: zflat_org,strhtune_org
INTEGER :: strhopt_org
CHARACTER (LEN=80) :: grdbasfn_ext
CHARACTER (LEN=80) :: datafn_ext
INTEGER :: nchanl_ext,lengbf_ext
INTEGER :: lendtf_ext
REAL :: time_ext
!
!-----------------------------------------------------------------------
!
! Include files
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'grid.inc'
INCLUDE 'phycst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
! Build file names
!
!-----------------------------------------------------------------------
!
IF ( extdfcst == ' ') extdfcst='000:00:00'
lenrun=LEN(dir_extd)
ldir=lenrun
CALL strlnth
( dir_extd, ldir )
IF ( ldir == 0 .OR. dir_extd(1:ldir) == ' ' ) THEN
dir_extd = '.'
ldir = 1
END IF
IF( dir_extd(ldir:ldir) /= '/' .AND. ldir < lenrun ) THEN
ldir = ldir + 1
dir_extd(ldir:ldir) = '/'
END IF
lenrun = LEN( extdname )
CALL strlnth( extdname, lenrun )
IF( extdfmt == 1 ) THEN
fmtn = 'bin'
ELSE IF ( extdfmt == 2 ) THEN
fmtn = 'asc'
ELSE IF ( extdfmt == 3 ) THEN
fmtn = 'hdf'
ELSE IF ( extdfmt == 4 ) THEN
fmtn = 'pak'
ELSE IF ( extdfmt == 6 ) THEN
fmtn = 'bn2'
ELSE IF ( extdfmt == 7 ) THEN
fmtn = 'net'
ELSE IF ( extdfmt == 8 ) THEN
fmtn = 'npk'
ELSE IF ( extdfmt == 9 ) THEN
fmtn = 'gad'
ELSE IF ( extdfmt == 10 ) THEN
fmtn = 'grb'
ELSE
WRITE(6,'(a,a,a)') &
'Unknown format, ', extdfmt, '. Program stopped in GETARPS.'
STOP
END IF
READ(extdfcst,'(i3,1x,i2,1x,i2)') ihr,imin,isec
time_ext = FLOAT( (ihr*3600)+(imin*60)+isec )
CALL cvttsnd( time_ext, timsnd, tmstrln )
grdbasfn_ext = dir_extd(1:ldir)//extdname(1:lenrun) &
//'.'//fmtn//'grdbas'
lengbf_ext = ldir + lenrun + 10
datafn_ext = dir_extd(1:ldir)//extdname(1:lenrun) &
//'.'//fmtn//timsnd(1:tmstrln)
lendtf_ext = ldir + lenrun + 4 + tmstrln
WRITE(6,*) 'The external grid and base file, grdbasfn = ', &
grdbasfn_ext(1:lendtf_ext)
WRITE(6,*) 'The external time dependent file, datafn = ', &
datafn_ext(1:lengbf_ext)
!
!-----------------------------------------------------------------------
!
! Since the data reader will change certain parameters stored
! in common, they need to be saved and restored in common
! after reading is done.
!
!-----------------------------------------------------------------------
!
runname_org=runname
nocmnt_org=nocmnt
IF( nocmnt > 0 ) THEN
DO i=1,nocmnt
cmnt_org(i)=cmnt(i)
END DO
END IF
mapproj_org=mapproj
month_org=month
day_org=day
year_org=year
hour_org=hour
minute_org=minute
second_org=second
!
umove_org=umove
vmove_org=vmove
trulat1_org=trulat1
trulat2_org=trulat2
trulon_org=trulon
sclfct_org=sclfct
latitud_org=latitud
ctrlat_org=ctrlat
ctrlon_org=ctrlon
dx_org=dx
dy_org=dy
xgrdorg_org=xgrdorg
ygrdorg_org=ygrdorg
dz_org=dz
dzmin_org=dzmin
zrefsfc_org=zrefsfc
dlayer1_org=dlayer1
dlayer2_org=dlayer2
zflat_org=zflat
strhtune_org=strhtune
strhopt_org=strhopt
CALL xytoll(1,1,0.,0.,lat_org,lon_org)
!
!-----------------------------------------------------------------------
!
! Read ARPS data file
!
!-----------------------------------------------------------------------
!
! uatv_ext & vatu_ext used as temporary variables here
!
CALL dtaread(nx_ext,ny_ext,nz_ext,nzsoil_ext,nstyps, &
extdfmt, nchanl_ext, grdbasfn_ext,lengbf_ext, &
datafn_ext, lendtf_ext, time_ext, &
x_ext,y_ext,z_ext,zp_ext,zpsoil_ext, &
u_ext,v_ext,w_ext,t_ext,p_ext, &
qv_ext, qc_ext, qr_ext, &
qi_ext, qs_ext, qh_ext, vatu_ext,vatu_ext,vatu_ext, &
ubar_ext, vbar_ext, wbar_ext, &
ptbar_ext, pbar_ext, vatu_ext, qvbar_ext, &
soiltyp_ext,stypfrct_ext,vegtyp_ext, &
lai_ext,roufns_ext,veg_ext, &
tsoil_ext,qsoil_ext,wetcanp_ext, &
snowdpth_ext, &
tem1_ext(1,1,1),tem1_ext(1,1,2),tem2_ext, &
tem2_ext,tem1_ext(1,1,1),tem1_ext(1,1,2), &
tem2_ext,tem2_ext, &
tem1_ext(1,1,1),tem1_ext(1,1,2), &
tem1_ext(1,1,3),tem1_ext(1,1,4), &
ireturn,tem1_ext,tem2_ext,uatv_ext)
!
!-----------------------------------------------------------------------
!
! Set maprojection parameters
!
!-----------------------------------------------------------------------
!
iproj_ext=mapproj
scale_ext=sclfct
trlon_ext=trulon
latnot_ext(1)=trulat1
latnot_ext(2)=trulat2
CALL setmapr(iproj_ext,scale_ext,latnot_ext,trlon_ext)
CALL lltoxy(1,1,ctrlat,ctrlon,xctr,yctr)
dx_ext=x_ext(2)-x_ext(1)
x0_ext=xctr - 0.5*(nx_ext-3)*dx_ext
dy_ext=y_ext(2)-y_ext(1)
y0_ext=yctr - 0.5*(ny_ext-3)*dy_ext
CALL setorig(1,x0_ext,y0_ext)
!
!-----------------------------------------------------------------------
!
! Find lat,lon of scalar points
!
!-----------------------------------------------------------------------
!
DO i=1,nx_ext-1
xsc(i)=0.5*(x_ext(i)+x_ext(i+1))
END DO
xsc(nx_ext)=2.*xsc(nx_ext-1)-xsc(nx_ext-2)
DO j=1,ny_ext-1
ysc(j)=0.5*(y_ext(j)+y_ext(j+1))
END DO
ysc(ny_ext)=2.*ysc(ny_ext-1)-ysc(ny_ext-2)
CALL xytoll(nx_ext,ny_ext,xsc,ysc,lat_ext,lon_ext)
CALL xytoll(nx_ext,ny_ext,x_ext,ysc,latu_ext,lonu_ext)
CALL xytoll(nx_ext,ny_ext,xsc,y_ext,latv_ext,lonv_ext)
!
!-----------------------------------------------------------------------
!
! Move z field onto the scalar grid.
!
!-----------------------------------------------------------------------
!
DO k=1,nz_ext-1
DO j=1,ny_ext-1
DO i=1,nx_ext-1
hgt_ext(i,j,k)=0.5*(zp_ext(i,j,k)+zp_ext(i,j,k+1))
END DO
END DO
END DO
DO j=1,ny_ext-1
DO i=1,nx_ext-1
hgt_ext(i,j,nz_ext)=(2.*hgt_ext(i,j,nz_ext-1)) &
-hgt_ext(i,j,nz_ext-2)
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Combine wind perturbations and mean fields.
!
!-----------------------------------------------------------------------
!
u_ext = u_ext + ubar_ext
v_ext = v_ext + vbar_ext
w_ext = w_ext + wbar_ext
!
!-----------------------------------------------------------------------
!
! Orient u & v to true north.
!
!-----------------------------------------------------------------------
!
DO k=1,nz_ext
! get u at v grid point locations
DO j=2,ny_ext
DO i=1,nx_ext-1
uatv_ext(i,j,k) = 0.25*(u_ext(i,j-1,k)+u_ext(i+1,j-1,k) &
+u_ext(i,j,k)+u_ext(i+1,j,k))
END DO
END DO
DO i=1,nx_ext-1
uatv_ext(i,1,k) = 0.5*(u_ext(i,1,k)+u_ext(i+1,1,k))
END DO
DO j=2,ny_ext
uatv_ext(nx_ext,j,k) = 0.5*(u_ext(nx_ext,j-1,k)+u_ext(nx_ext,j,k))
END DO
uatv_ext(nx_ext,1,k) = u_ext(nx_ext,1,k)
! get v at u grid point locations
DO j=1,ny_ext-1
DO i=2,nx_ext
vatu_ext(i,j,k) = 0.25*(v_ext(i-1,j,k)+v_ext(i,j,k) &
+v_ext(i-1,j+1,k)+v_ext(i,j+1,k))
END DO
END DO
DO j=1,ny_ext-1
vatu_ext(1,j,k) = 0.5*(v_ext(1,j,k)+v_ext(1,j+1,k))
END DO
DO i=2,nx_ext
vatu_ext(i,ny_ext,k) = 0.5*(v_ext(i-1,ny_ext,k)+v_ext(i,ny_ext,k))
END DO
vatu_ext(1,ny_ext,k) = v_ext(1,ny_ext,k)
CALL uvmptoe(nx_ext,ny_ext,u_ext(1,1,k),vatu_ext(1,1,k),lonu_ext, &
tem1_ext(1,1,k),tem2_ext(1,1,k))
u_ext(1:nx_ext,1:ny_ext,k) = tem1_ext(1:nx_ext,1:ny_ext,k)
vatu_ext(1:nx_ext,1:ny_ext,k) = tem2_ext(1:nx_ext,1:ny_ext,k)
CALL uvmptoe(nx_ext,ny_ext,uatv_ext(1,1,k),v_ext(1,1,k),lonv_ext, &
tem1_ext(1,1,k),tem2_ext(1,1,k))
uatv_ext(1:nx_ext,1:ny_ext,k) = tem1_ext(1:nx_ext,1:ny_ext,k)
v_ext(1:nx_ext,1:ny_ext,k) = tem2_ext(1:nx_ext,1:ny_ext,k)
END DO
!
!-----------------------------------------------------------------------
!
! Combine perturbations and mean fields of scalars
! Convert potential temperature to potential temperature
!
!-----------------------------------------------------------------------
!
DO k=2,nz_ext-1
DO j=2,ny_ext-1
DO i=2,nx_ext-1
p_ext(i,j,k)=p_ext(i,j,k)+pbar_ext(i,j,k)
t_ext(i,j,k)=t_ext(i,j,k)+ptbar_ext(i,j,k)
t_ext(i,j,k)=t_ext(i,j,k)*((p_ext(i,j,k)/p0)**rddcp)
qv_ext(i,j,k)=qv_ext(i,j,k)+qvbar_ext(i,j,k)
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Supply data at the edge points (zero gradient, where missing)
!
!-----------------------------------------------------------------------
!
CALL edgfill(hgt_ext,1,nx_ext,2,nx_ext-1,1,ny_ext,2,ny_ext-1, &
1,nz_ext,1,nz_ext)
CALL edgfill(p_ext,1,nx_ext,2,nx_ext-1,1,ny_ext,2,ny_ext-1, &
1,nz_ext,1,nz_ext)
CALL edgfill(t_ext,1,nx_ext,2,nx_ext-1,1,ny_ext,2,ny_ext-1, &
1,nz_ext,1,nz_ext)
CALL edgfill(qv_ext,1,nx_ext,2,nx_ext-1,1,ny_ext,2,ny_ext-1, &
1,nz_ext,2,nz_ext-1)
! u & v now not on scalar points, so don't edgfill
! CALL edgfill(u_ext,1,nx_ext,2,nx_ext-1,1,ny_ext,2,ny_ext-1, &
! 1,nz_ext,2,nz_ext-1)
! CALL edgfill(v_ext,1,nx_ext,2,nx_ext-1,1,ny_ext,2,ny_ext-1, &
! 1,nz_ext,2,nz_ext-1)
!
!-----------------------------------------------------------------------
!
! Make top and bottom mass fields via hydrostatic extrapolation.
!
!-----------------------------------------------------------------------
!
govrd=g/rd
DO j=1,ny_ext-1
DO i=1,nx_ext-1
!
t_ext(i,j,1)=2.*t_ext(i,j,2)-t_ext(i,j,3)
tvbot=t_ext(i,j,1) * ( 1.0 + 0.608*qv_ext(i,j,1) )
tvtop=t_ext(i,j,2) * ( 1.0 + 0.608*qv_ext(i,j,2) )
tvbar=0.5*(tvtop+tvbot)
p_ext(i,j,1)=p_ext(i,j,2) &
*EXP(govrd*(hgt_ext(i,j,2)-hgt_ext(i,j,1))/tvbar)
!
t_ext(i,j,nz_ext)=2.*t_ext(i,j,nz_ext-1)-t_ext(i,j,nz_ext-2)
tvbot=t_ext(i,j,nz_ext-1)*(1.0 + 0.608*qv_ext(i,j,nz_ext-1))
tvtop=t_ext(i,j,nz_ext) *(1.0 + 0.608*qv_ext(i,j,nz_ext))
tvbar=0.5*(tvtop+tvbot)
p_ext(i,j,nz_ext)=p_ext(i,j,nz_ext-1) &
*EXP(govrd*(hgt_ext(i,j,nz_ext-1)- &
hgt_ext(i,j,nz_ext))/tvbar)
END DO
END DO
CALL edgfill(p_ext,1,nx_ext,1,nx_ext-1,1,ny_ext,1,ny_ext-1, &
1,nz_ext,1,nz_ext)
CALL edgfill(t_ext,1,nx_ext,1,nx_ext-1,1,ny_ext,1,ny_ext-1,1 &
,nz_ext,1,nz_ext)
!
!-----------------------------------------------------------------------
!
! Reset info in common to original values
!
!-----------------------------------------------------------------------
!
runname=runname_org
nocmnt=nocmnt_org
IF( nocmnt > 0 ) THEN
DO i=1,nocmnt
cmnt(i)=cmnt_org(i)
END DO
END IF
mapproj=mapproj_org
month=month_org
day=day_org
year=year_org
hour=hour_org
minute=minute_org
second=second_org
!
umove=umove_org
vmove=vmove_org
xgrdorg=xgrdorg_org
ygrdorg=ygrdorg_org
trulat1=trulat1_org
trulat2=trulat2_org
trulon=trulon_org
sclfct=sclfct_org
latitud=latitud_org
ctrlat=ctrlat_org
ctrlon=ctrlon_org
dx=dx_org
dy=dy_org
dz=dz_org
dzmin=dzmin_org
zrefsfc=zrefsfc_org
dlayer1=dlayer1_org
dlayer2=dlayer2_org
zflat=zflat_org
strhtune=strhtune_org
strhopt=strhopt_org
xgrdorg=xgrdorg_org
ygrdorg=ygrdorg_org
!
!-----------------------------------------------------------------------
!
! Reset map projection to previous values
!
!-----------------------------------------------------------------------
!
latnot(1)=trulat1
latnot(2)=trulat2
CALL setmapr(mapproj,sclfct,latnot,trulon)
CALL setorig(2,lat_org,lon_org)
istatus=1
RETURN
!
!-----------------------------------------------------------------------
!
! Error destination
!
!-----------------------------------------------------------------------
!
! 598 CONTINUE
! WRITE(6,'(a)') ' Error reading last field, returning'
! RETURN
END SUBROUTINE getarps
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE GETNMCRUC87 ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE getnmcruc87(nx_ext,ny_ext,nz_ext, & 1,12
dir_extd,extdname,extdopt,extdfmt, &
extdinit,extdfcst,julfname, &
iproj_ext,scale_ext, &
trlon_ext,latnot_ext,x0_ext,y0_ext, &
lat_ext,lon_ext, &
p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext, &
qc_ext,qr_ext,qi_ext,qs_ext,qh_ext, &
tsfc_ext,tdeep_ext,wetsfc_ext,wetdp_ext,wetcanp_ext, &
trn_ext,psfc_ext, &
istatus)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! ARPS version.
!
! Reads an ARPS file for processing by ext2arps, a program
! that converts external files to ARPS variables and format.
! This version is useful when you want to use an ARPS file
! with a different orientation or terrain than your final
! ARPS product so arpsr2h does not work.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Yuhe Liu and Jinxing Zong
! 01/18/1996
!
! MODIFICATION HISTORY:
! 6/5/1997 Jinxing Zong and Yuhe Liu
! Virtualization of temperature is accounted for in variable
! conversion. Values of constants cp, rd/cp and g used in RUC are
! adopted in making the variable conversion. Subroutine TV2TQ and
! function ESW are added to facilitate the conversion.
!
! 01/26/1998 (Yuhe Liu)
! Removed function esw, instead, used unified function f_es in
! file thermolib3d.f
!
! 1999/11/30 (Gene Bassett)
! Changed deep soil temperature and moisture to be an average from
! 0-100cm.
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! dir_extd Directory name for external file
! extdname Prefix string of external file name
! extdopt Option of external data sources
! extdfmt Option of external data format
! extdinit Initialized time in mm-dd-yyyy:hh:mm:ss format
! extdfcst Forecast hour in HHH:MM:SS format
! julfname File name in yyjjjhhmm format
!
! OUTPUT:
!
! iproj_ext Map projection number of external data
! scale_ext Scale factor of external data
! trlon_ext True longitude of external data (degrees E)
! latnot_ext(2) True latitude(s) of external data (degrees N)
! x0_ext x coordinate of origin of external data
! y0_ext y coordinate of origin of external data
! lat_ext latitude of external data points (degrees N)
! lon_ext longitude of external data points (degrees E)
! p_ext pressure (Pascal)
! hgt_ext height (m)
! t_ext temperature (K)
! qv_ext specific humidity (kg/kg)
! u_ext u wind component (m/s)
! v_ext v wind component (m/s)
! qc_ext Cloud water mixing ratio (kg/kg)
! qr_ext Rain water mixing ratio (kg/kg)
! qi_ext Ice mixing ratio (kg/kg)
! qs_ext Snow mixing ratio (kg/kg)
! qh_ext Hail mixing ratio (kg/kg)
!
! tsfc_ext Surface temperature
! tdeep_ext Soil temperature
! wetsfc_ext Top layer soil moisture (fraction)
! wetdp_ext Deep soil moisture (fraction)
! wetcanp_ext Water content on canopy
!
! trn_ext External terrain (m)
! psfc_ext Surface pressure (Pa)
!
! istatus status indicator
!
! WORK ARRAYS:
!
! var_grb3d Arrays to store the GRIB 3-D variables:
! var_grb3d(nxgrb,nygrb,nzgrb,1,1) - pressure (Pa)
! var_grb3d(nxgrb,nygrb,nzgrb,2,1) - Montgomery stream
! function (m**2/s**2)
! var_grb3d(nxgrb,nygrb,nzgrb,3,1) - Potential
! temperature (K)
! var_grb3d(nxgrb,nygrb,nzgrb,4,1) - Condensation pressure
! of lifted parcel (Pa)
! var_grb3d(nxgrb,nygrb,nzgrb,5,1) - u wind (m/s)
! var_grb3d(nxgrb,nygrb,nzgrb,6,1) - v wind (m/s)
! var_grb3d(nxgrb,nygrb,nzgrb,1,2) - soil temp. (K)
! var_grb3d(nxgrb,nygrb,nzgrb,2,2) - soil moist.
! (fraction)
!
! var_grb2d Arrays to store the GRIB 2-D variables:
! var_grb2d(nxgrb,nygrb,1) - Ground surface temperature (K)
! var_grb2d(nxgrb,nygrb,2) - Geometric height (m)
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INCLUDE 'gribcst.inc'
CHARACTER (LEN=*) :: dir_extd
CHARACTER (LEN=*) :: extdname
INTEGER :: extdopt
INTEGER :: extdfmt
CHARACTER (LEN=19) :: extdinit
CHARACTER (LEN=9) :: extdfcst
CHARACTER (LEN=9) :: julfname
!
!-----------------------------------------------------------------------
!
! External grid variables
!
!-----------------------------------------------------------------------
!
INTEGER :: iproj_ext
REAL :: scale_ext,trlon_ext
REAL :: latnot_ext(2)
REAL :: x0_ext,y0_ext
REAL :: dx_ext,dy_ext
!
!-----------------------------------------------------------------------
!
! Output external variable arrays
!
!-----------------------------------------------------------------------
!
INTEGER :: nx_ext, ny_ext, nz_ext
REAL :: lat_ext(nx_ext,ny_ext)
REAL :: lon_ext(nx_ext,ny_ext)
REAL :: p_ext (nx_ext,ny_ext,nz_ext) ! Pressure (Pascals)
REAL :: hgt_ext(nx_ext,ny_ext,nz_ext) ! Height (m)
REAL :: t_ext (nx_ext,ny_ext,nz_ext) ! Temperature (K)
REAL :: qv_ext (nx_ext,ny_ext,nz_ext) ! Specific humidity (kg/kg)
REAL :: u_ext (nx_ext,ny_ext,nz_ext) ! Eastward wind component
REAL :: v_ext (nx_ext,ny_ext,nz_ext) ! Northward wind component
REAL :: qc_ext (nx_ext,ny_ext,nz_ext) ! Cloud H2O mixing ratio (kg/kg)
REAL :: qr_ext (nx_ext,ny_ext,nz_ext) ! Rain H2O mixing ratio (kg/kg)
REAL :: qi_ext (nx_ext,ny_ext,nz_ext) ! Ice mixing ratio (kg/kg)
REAL :: qs_ext (nx_ext,ny_ext,nz_ext) ! Snow mixing ratio (kg/kg)
REAL :: qh_ext (nx_ext,ny_ext,nz_ext) ! Hail mixing ratio (kg/kg)
REAL :: tsfc_ext (nx_ext,ny_ext) ! Temperature at surface (K)
REAL :: tdeep_ext (nx_ext,ny_ext) ! Deep soil temperature (K)
REAL :: wetsfc_ext (nx_ext,ny_ext) ! Surface soil moisture (fraction)
REAL :: wetdp_ext (nx_ext,ny_ext) ! Deep soil moisture (fraction)
REAL :: wetcanp_ext(nx_ext,ny_ext) ! Canopy water amount
REAL :: trn_ext (nx_ext,ny_ext) ! Geometrical heights
REAL :: psfc_ext (nx_ext,ny_ext) ! Surface pressure (Pa)
!
!-----------------------------------------------------------------------
!
! Work arrays for storing grib data
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: var_grb2d(:,:,:,:) ! GRIB variables
REAL, ALLOCATABLE :: var_grb3d(:,:,:,:,:) ! GRIB 3-D variables
INTEGER, ALLOCATABLE :: var_lev3d(:,:,:) ! Levels (hybrid) for
! each 3-D variable
REAL, ALLOCATABLE :: rcdata(:) ! temporary data array
!
!-----------------------------------------------------------------------
!
! Other external variable arrays
!
!-----------------------------------------------------------------------
!
REAL :: x_ext(nx_ext)
REAL :: y_ext(ny_ext)
INTEGER :: istatus
!
!-----------------------------------------------------------------------
!
! Original grid variables
!
!-----------------------------------------------------------------------
!
INTEGER :: iproj
REAL :: scale,trlon,x0,y0
REAL :: latnot(2)
!
!-----------------------------------------------------------------------
!
! Misc internal variables
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=80) :: gribfile
CHARACTER (LEN=13) :: gribtime
INTEGER :: i,j,k
INTEGER :: iyr,imo,iday,myr, jldy
INTEGER :: ihr,imin,isec
INTEGER :: ifhr,ifmin,ifsec
INTEGER :: grbflen, len1, lenrun
INTEGER :: m,n,nz1,max_nr2d,max_nr3d
REAL :: pilev
REAL :: tv_ext, tvc_ext
REAL :: rovcp_p, cpd_p, g0_p, rd_p
INTEGER :: chklev, lvscan, kk, jj
REAL :: tema, temb
INTEGER :: iret ! Return flag
REAL, ALLOCATABLE :: utmp(:,:), vtmp(:,:)
!
!-----------------------------------------------------------------------
!
! Include files
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
! GRIB grid information
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=42) :: gridesc ! Grid description
INTEGER :: iproj_grb ! Map projection indicator
INTEGER :: gthin ! Indicator of whether the grid is "thinned"
INTEGER :: ni_grb ! Number of points along x-axis
INTEGER :: nj_grb ! Number of points along y-axis
INTEGER :: np_grb ! Total number of horizontal grid points
INTEGER :: nk_grb ! Number of vertical parameters
REAL :: zk_grb(nz_ext) ! Vertical coordinate parameters
INTEGER :: npeq ! Number of lat circles from pole to equator
INTEGER :: nit(nz_ext) ! Number of x-points for thinned grid
REAL :: pi_grb ! x-coordinate of pole point
REAL :: pj_grb ! y-coordinate of pole point
INTEGER :: ipole ! Projection center flag
REAL :: di_grb ! x-direction increment or grid length
REAL :: dj_grb ! y-direction increment or grid length
REAL :: latsw ! Latitude of South West corner point
REAL :: lonsw ! Longitude of South West corner point
REAL :: latne ! Latitude of North East corner point
REAL :: lonne ! Longitude of North East corner point
REAL :: lattru1 ! Latitude (1st) at which projection is true
REAL :: lattru2 ! Latitude (2nd) at which projection is true
REAL :: lontrue ! Longitude at which projection is true
REAL :: latrot ! Latitude of southern pole of rotation
REAL :: lonrot ! Longitude of southern pole of rotation
REAL :: angrot ! Angle of rotation
REAL :: latstr ! Latitude of the pole of stretching
REAL :: lonstr ! Longitude of the pole of stretching
REAL :: facstr ! Stretching factor
INTEGER :: scanmode ! Scanning indicator
INTEGER :: iscan ! x-direction scanning indicator
INTEGER :: jscan ! y-direction scanning indicator
INTEGER :: kscan ! FORTRAN index scanning indicator
INTEGER :: ires ! Resolution direction increments indicator
INTEGER :: iearth ! Earth shape indicator: spherical or oblate?
INTEGER :: icomp ! (u,v) components decomposition indicator
INTEGER :: jpenta ! J-Pentagonal resolution parameter
INTEGER :: kpenta ! K-Pentagonal resolution parameter
INTEGER :: mpenta ! M-Pentagonal resolution parameter
INTEGER :: ispect ! Spectral representation type
INTEGER :: icoeff ! Spectral coefficient storage mode
REAL :: xp_grb ! X coordinate of sub-satellite point
REAL :: yp_grb ! Y coordinate of sub-satellite point
REAL :: xo_grb ! X coordinate of image sector origin
REAL :: yo_grb ! Y coordinate of image sector origin
REAL :: zo_grb ! Camera altitude from center of Earth
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
ALLOCATE(var_grb2d(nx_ext,ny_ext,n2dvs,n2dlvt))
ALLOCATE(var_grb3d(nx_ext,ny_ext,nz_ext,n3dvs,n3dlvt))
ALLOCATE(rcdata(nx_ext*ny_ext))
ALLOCATE(var_lev3d(nz_ext,n3dvs,n3dlvt))
ALLOCATE(utmp(nx_ext,ny_ext))
ALLOCATE(vtmp(nx_ext,ny_ext))
IF ( extdfcst == ' ') extdfcst='000:00:00'
READ (extdinit,'(i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,i2)') &
iyr,imo,iday,ihr,imin,isec
CALL julday(iyr,imo,iday,jldy)
myr=MOD(iyr,100)
ifhr=0
ifmin=0
ifsec=0
READ (extdfcst,'(i3)',ERR=4,END=4) ifhr
4 CONTINUE
WRITE (gribtime,'(i4.4,i2.2,i2.2,i2.2,a1,i2.2)') &
iyr,imo,iday,ihr,'f',ifhr
len1=LEN(dir_extd)
grbflen=len1
CALL strlnth( dir_extd, grbflen )
IF( grbflen == 0 .OR. dir_extd(1:grbflen) == ' ' ) THEN
dir_extd = '.'
grbflen=1
END IF
IF( dir_extd(grbflen:grbflen) /= '/'.AND.grbflen < len1 ) THEN
grbflen=grbflen+1
dir_extd(grbflen:grbflen)='/'
END IF
lenrun = LEN( extdname )
CALL strlnth( extdname, lenrun )
gribfile = dir_extd(1:grbflen)//extdname(1:lenrun) &
//'.'//gribtime(1:13)
grbflen = grbflen + lenrun + 14
!
!-----------------------------------------------------------------------
!
! RDNMCGRB reads NMC GRIB data
!
!-----------------------------------------------------------------------
!
gridtyp = ruc87grid
mproj_grb = ruc87proj
n2dvars = ruc87nvs2d
n2dlvtps = ruc87nlvt2d
DO m=1,n2dlvtps
DO n=1,n2dvars
var_id2d(n,m) = ruc87var_id2d(n,m)
END DO
levtyp2d(m) = ruc87levs2d(m)
END DO
n3dvars = ruc87nvs3d
n3dlvtps = ruc87nlvt3d
DO m=1,n3dlvtps
DO n=1,n3dvars
var_id3d(n,m) = ruc87var_id3d(n,m)
END DO
levtyp3d(m) = ruc87levs3d(m)
END DO
CALL rdnmcgrb(nx_ext,ny_ext,nz_ext,gribfile,grbflen, gribtime, &
gridesc, iproj_grb, gthin, &
ni_grb,nj_grb,np_grb, nk_grb,zk_grb, npeq,nit, &
pi_grb,pj_grb,ipole, di_grb,dj_grb, &
latsw,lonsw, latne,lonne, &
latrot,lonrot,angrot, &
latstr,lonstr,facstr, &
lattru1,lattru2,lontrue, &
scanmode, iscan,jscan,kscan, &
ires,iearth,icomp, &
jpenta,kpenta,mpenta,ispect,icoeff, &
xp_grb,yp_grb, xo_grb,yo_grb,zo_grb, &
rcdata,var_grb2d,var_grb3d,var_lev3d,iret)
max_nr2d = 0
DO n=1,n2dvars
DO m=1,n2dlvtps
max_nr2d = MAX( max_nr2d, var_nr2d(n,m) )
END DO
END DO
max_nr3d = 0
DO n=1,n3dvars
max_nr3d = MAX( max_nr3d, var_nr3d(n,1) )
END DO
IF ( max_nr3d == 0 ) THEN
WRITE (6,'(a)') &
'No 3-D variables was found in the GRIB file', &
'Program stopped in GETNMCRUC.'
STOP
END IF
IF ( max_nr2d == 0 ) THEN
WRITE (6,'(a)') &
'No 2-D variables was found in the GRIB file'
END IF
! write (6,'(/a7,2x,6(i7))')
! : 'Lev\\VID',(var_id3d(n,1),n=1,n3dvars)
! DO 60 k=1,max_nr3d
! write (6,'(i5,4x,6(i7))')
! : k,(var_lev3d(k,n,1),n=1,n3dvars)
60 CONTINUE
DO k=1,max_nr3d
DO n=2,n3dvars
IF ( var_lev3d(k,1,1) /= var_lev3d(k,n,1) ) THEN
WRITE (6,'(a)') &
'Variables were not at the same level.', &
'Program stopped in GETNMCRUC.'
STOP
END IF
END DO
END DO
IF ( iproj_grb == 5 .AND. ipole == 0 ) THEN ! Center in NP
iproj_ext = 1
ELSE IF ( iproj_grb == 5 .AND. ipole == 1 ) THEN ! Center in SP
iproj_ext = -1
ELSE IF ( iproj_grb == 3 .AND. ipole == 0 ) THEN ! Center in NP
iproj_ext = 2
ELSE IF ( iproj_grb == 3 .AND. ipole == 1 ) THEN ! Center in SP
iproj_ext = -2
ELSE IF ( iproj_grb == 1 ) THEN
iproj_ext = 3
ELSE IF ( iproj_grb == 0 ) THEN
iproj_ext = 4
ELSE
WRITE (6,'(a)') &
'Unknown map projection. Set to non-projection.'
iproj_ext = 0
END IF
scale_ext = 1.0
latnot_ext(1) = lattru1
latnot_ext(2) = lattru2
trlon_ext = lontrue
dx_ext = di_grb
dy_ext = dj_grb
CALL getmapr(iproj,scale,latnot,trlon,x0,y0)
CALL setmapr(iproj_ext,scale_ext,latnot_ext,trlon_ext)
CALL lltoxy(1,1,latsw,lonsw,x0_ext,y0_ext)
DO i=1,nx_ext
x_ext(i)=x0_ext+(i-1)*dx_ext
END DO
DO j=1,ny_ext
y_ext(j)=y0_ext+(j-1)*dy_ext
END DO
CALL xytoll(nx_ext,ny_ext,x_ext,y_ext,lat_ext,lon_ext)
!
!-----------------------------------------------------------------------
!
! Retrieve 2-D variables
!
!-----------------------------------------------------------------------
!
DO j=1,ny_ext
DO i=1,nx_ext
IF ( var_nr2d(1,1) == 0 ) THEN
tsfc_ext (i,j) = -999.0
ELSE
tsfc_ext (i,j) = var_grb2d(i,j,1,1)
END IF
IF ( var_nr2d(2,1) == 0 ) THEN
trn_ext (i,j) = -999.0
ELSE
trn_ext (i,j) = var_grb2d(i,j,2,1)
END IF
IF ( var_nr3d(1,2) == 0 ) THEN
tsfc_ext (i,j) = var_grb3d(i,j,1,1,2)
! note that this is the 5cm value and not the surface value
tdeep_ext (i,j) = 0.1 * var_grb3d(i,j,1,1,2) & ! 5cm
+ 0.2 * var_grb3d(i,j,2,1,2) & ! 20cm
+ 0.4 * var_grb3d(i,j,3,1,2) & ! 40cm
+ 0.3 * var_grb3d(i,j,4,1,2) ! 160cm
wetsfc_ext (i,j) = var_grb3d(i,j,1,2,2)
! note that this is the 5cm value and not the surface value
wetdp_ext (i,j) = 0.1 * var_grb3d(i,j,1,2,2) & ! 5cm
+ 0.2 * var_grb3d(i,j,2,2,2) & ! 20cm
+ 0.4 * var_grb3d(i,j,3,2,2) & ! 40cm
+ 0.3 * var_grb3d(i,j,4,2,2) ! 160cm
wetcanp_ext(i,j) = var_grb2d(i,j,3,1)*1.e-3 ! in meters
ELSE
tsfc_ext (i,j) = -999.0
tdeep_ext (i,j) = -999.0
wetsfc_ext (i,j) = -999.0
wetdp_ext (i,j) = -999.0
wetcanp_ext(i,j) = -999.0
END IF
psfc_ext (i,j) = -999.0
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Retrieve 3-D variables
!
!-----------------------------------------------------------------------
!
cpd_p = 1004.686 ! cp in RUC
rovcp_p = 0.285714 ! rd/cp used in RUC
g0_p = 9.80665 ! gravity in RUC
nz1 = MIN(var_nr3d(1,1),nz_ext)
IF ( var_lev3d(1,1,1) < var_lev3d(nz1,1,1) ) THEN ! 1st level at sfc
chklev = 1
lvscan = 0
ELSE
chklev = -1
lvscan = nz1+1
END IF
DO k=1,nz1
kk = chklev * k + lvscan
DO j=1,ny_ext
DO i=1,nx_ext
p_ext(i,j,kk) = var_grb3d(i,j,k,1,1) ! Pressure (Pa)
pilev = (p_ext(i,j,kk)/100000.)**rovcp_p
tv_ext = var_grb3d(i,j,k,3,1)*pilev ! Virtual Temperature (K)
hgt_ext(i,j,kk) = (var_grb3d(i,j,k,2,1)-cpd_p*tv_ext)/g0_p
! Height (m) from Mongmery function with
! M = CpTv + gz
u_ext(i,j,kk) = var_grb3d(i,j,k,5,1) ! u wind (m/s)
v_ext(i,j,kk) = var_grb3d(i,j,k,6,1) ! v wind (m/s)
tvc_ext = var_grb3d(i,j,k,3,1) &
* (var_grb3d(i,j,k,4,1)/100000.)**rovcp_p
! Virtual Temperature (K) at LCL
CALL tv2tq( tvc_ext, var_grb3d(i,j,k,4,1), tema, temb )
t_ext(i,j,kk) = tema &
* (p_ext(i,j,kk)/var_grb3d(i,j,k,4,1))**rovcp_p ! Temperature (K)
qv_ext(i,j,kk) = temb ! Specific humidity
qc_ext(i,j,kk) = -999.
qr_ext(i,j,kk) = -999.
qi_ext(i,j,kk) = -999.
qs_ext(i,j,kk) = -999.
qh_ext(i,j,kk) = -999.
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Rotate winds to be relative to true north.
! The RUC data are sent as grid-relative.
!
!-----------------------------------------------------------------------
!
DO k=1,nz1
CALL uvmptoe(nx_ext,ny_ext,u_ext(1,1,k),v_ext(1,1,k), &
lon_ext,utmp,vtmp)
u_ext(:,:,k) = utmp(:,:)
v_ext(:,:,k) = vtmp(:,:)
END DO
!
!-----------------------------------------------------------------------
!
! Reset map projection to previous values
!
!-----------------------------------------------------------------------
!
CALL setmapr(iproj,scale,latnot,trlon)
CALL setorig(1,x0,y0)
istatus = 1
DEALLOCATE(var_grb2d,var_grb3d,rcdata,var_lev3d)
DEALLOCATE(utmp)
DEALLOCATE(vtmp)
RETURN
END SUBROUTINE getnmcruc87
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE GETNMCRUC211 ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE getnmcruc211(nx_ext,ny_ext,nz_ext, & 1,12
dir_extd,extdname,extdopt,extdfmt, &
extdinit,extdfcst,julfname, &
iproj_ext,scale_ext, &
trlon_ext,latnot_ext,x0_ext,y0_ext, &
lat_ext,lon_ext, &
p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext, &
qc_ext,qr_ext,qi_ext,qs_ext,qh_ext, &
tsfc_ext,tdeep_ext,wetsfc_ext,wetdp_ext,wetcanp_ext, &
trn_ext,psfc_ext, t_2m_ext, rh_2m_ext, u_10m_ext, &
v_10m_ext, istatus)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! ARPS version.
!
! Reads an ARPS file for processing by ext2arps, a program
! that converts external files to ARPS variables and format.
! This version is useful when you want to use an ARPS file
! with a different orientation or terrain than your final
! ARPS product so arpsr2h does not work.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Dennis Moon, SSESCO
! 06/19/1996
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! dir_extd Directory name for external file
! extdname Prefix string of external file name
! extdopt Option of external data sources
! extdfmt Option of external data format
! extdinit Initialized time in mm-dd-yyyy:hh:mm:ss format
! extdfcst Forecast hour in HHH:MM:SS format
! julfname File name in yyjjjhhmm format
!
! OUTPUT:
!
! iproj_ext Map projection number of external data
! scale_ext Scale factor of external data
! trlon_ext True longitude of external data (degrees E)
! latnot_ext(2) True latitude(s) of external data (degrees N)
! x0_ext x coordinate of origin of external data
! y0_ext y coordinate of origin of external data
! lat_ext latitude of external data points (degrees N)
! lon_ext longitude of external data points (degrees E)
! p_ext pressure (Pascal)
! hgt_ext height (m)
! t_ext temperature (K)
! qv_ext specific humidity (kg/kg)
! u_ext u wind component (m/s)
! v_ext v wind component (m/s)
! qc_ext Cloud water mixing ratio (kg/kg)
! qr_ext Rain water mixing ratio (kg/kg)
! qi_ext Ice mixing ratio (kg/kg)
! qs_ext Snow mixing ratio (kg/kg)
! qh_ext Hail mixing ratio (kg/kg)
!
! tsfc_ext Surface temperature
! tdeep_ext Soil temperature
! wetsfc_ext Top layer soil moisture (fraction)
! wetdp_ext Deep soil moisture (fraction)
! wetcanp_ext Water content on canopy
!
! trn_ext External terrain (m)
! psfc_ext Surface pressure (Pa)
!
! T_2m_ext Temperature at 2m AGL
! rh_2m_ext Specific Humidity at 2m AGL
! U_10m_ext U at 10m AGL
! V_10m_ext V at 10m AGL
!
! istatus status indicator
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INCLUDE 'gribcst.inc'
CHARACTER (LEN=*) :: dir_extd
CHARACTER (LEN=*) :: extdname
INTEGER :: extdopt
INTEGER :: extdfmt
CHARACTER (LEN=19) :: extdinit
CHARACTER (LEN=9) :: extdfcst
CHARACTER (LEN=9) :: julfname
!
!-----------------------------------------------------------------------
!
! External grid variables
!
!-----------------------------------------------------------------------
!
INTEGER :: iproj_ext
REAL :: scale_ext,trlon_ext
REAL :: latnot_ext(2)
REAL :: x0_ext,y0_ext
REAL :: dx_ext,dy_ext
!
!-----------------------------------------------------------------------
!
! Output external variable arrays
!
!-----------------------------------------------------------------------
!
INTEGER :: nx_ext,ny_ext,nz_ext
REAL :: lat_ext(nx_ext,ny_ext)
REAL :: lon_ext(nx_ext,ny_ext)
REAL :: p_ext (nx_ext,ny_ext,nz_ext) ! Pressure (Pascals)
REAL :: hgt_ext(nx_ext,ny_ext,nz_ext) ! Height (m)
REAL :: t_ext (nx_ext,ny_ext,nz_ext) ! Temperature (K)
REAL :: qv_ext (nx_ext,ny_ext,nz_ext) ! Specific humidity (kg/kg)
REAL :: u_ext (nx_ext,ny_ext,nz_ext) ! Eastward wind component
REAL :: v_ext (nx_ext,ny_ext,nz_ext) ! Northward wind component
REAL :: qc_ext (nx_ext,ny_ext,nz_ext) ! Cloud H2O mixing ratio (kg/kg)
REAL :: qr_ext (nx_ext,ny_ext,nz_ext) ! Rain H2O mixing ratio (kg/kg)
REAL :: qi_ext (nx_ext,ny_ext,nz_ext) ! Ice mixing ratio (kg/kg)
REAL :: qs_ext (nx_ext,ny_ext,nz_ext) ! Snow mixing ratio (kg/kg)
REAL :: qh_ext (nx_ext,ny_ext,nz_ext) ! Hail mixing ratio (kg/kg)
REAL :: tsfc_ext (nx_ext,ny_ext) ! Temperature at surface (K)
REAL :: tdeep_ext (nx_ext,ny_ext) ! Deep soil temperature (K)
REAL :: wetsfc_ext (nx_ext,ny_ext) ! Surface soil moisture (fraction)
REAL :: wetdp_ext (nx_ext,ny_ext) ! Deep soil moisture (fraction)
REAL :: wetcanp_ext(nx_ext,ny_ext) ! Canopy water amount
REAL :: trn_ext (nx_ext,ny_ext) ! Geometrical heights
REAL :: psfc_ext (nx_ext,ny_ext) ! Surface pressure (Pa)
REAL :: t_2m_ext (nx_ext,ny_ext)
REAL :: rh_2m_ext(nx_ext,ny_ext)
REAL :: u_10m_ext(nx_ext,ny_ext)
REAL :: v_10m_ext(nx_ext,ny_ext)
!
!----------------------------------------------------------------------
!
! Other external variable arrays
!
!----------------------------------------------------------------------
!
REAL :: x_ext(nx_ext)
REAL :: y_ext(ny_ext)
REAL :: z_ext(nz_ext)
INTEGER :: istatus
!
!-----------------------------------------------------------------------
!
! Work arrays for storing grib data
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: var_grb2d(:,:,:,:) ! GRIB variables
REAL, ALLOCATABLE :: var_grb3d(:,:,:,:,:) ! GRIB 3-D variables
INTEGER, ALLOCATABLE :: var_lev3d(:,:,:) ! Levels (hybrid) for
! each 3-D variable
REAL, ALLOCATABLE :: rcdata(:) ! temporary data array
!
!-----------------------------------------------------------------------
!
! Original grid variables
!
!-----------------------------------------------------------------------
!
INTEGER :: iproj
REAL :: scale,trlon,x0,y0
REAL :: latnot(2)
!
!-----------------------------------------------------------------------
!
! Misc internal variables
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=80) :: gribfile
CHARACTER (LEN=13) :: gribtime
CHARACTER (LEN=10) :: runstr
CHARACTER (LEN=3) :: fmtn
INTEGER :: ichr,bchar,echar
INTEGER :: i,j,k,ldir,ireturn
INTEGER :: iyr,imo,iday,myr, jldy
INTEGER :: ihr,imin,isec
INTEGER :: ifhr,ifmin,ifsec
INTEGER :: grbflen, len1, lenrun
INTEGER :: m,n,nz1,max_nr2d,max_nr3d
REAL :: qvsat, pilev, qvsatice
REAL :: tema, temb
INTEGER :: chklev, lvscan, kk, jj
INTEGER :: iret ! Return flag
REAL, ALLOCATABLE :: utmp(:,:), vtmp(:,:)
!
!-----------------------------------------------------------------------
!
! Include files
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
! GRIB grid information
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=42) :: gridesc ! Grid description
INTEGER :: iproj_grb ! Map projection indicator
INTEGER :: gthin ! Indicator of whether the grid is "thinned"
INTEGER :: ni_grb ! Number of points along x-axis
INTEGER :: nj_grb ! Number of points along y-axis
INTEGER :: np_grb ! Total number of horizontal grid points
INTEGER :: nk_grb ! Number of vertical parameters
REAL :: zk_grb(nz_ext) ! Vertical coordinate parameters
INTEGER :: npeq ! Number of lat circles from pole to equator
INTEGER :: nit(nz_ext) ! Number of x-points for thinned grid
REAL :: pi_grb ! x-coordinate of pole point
REAL :: pj_grb ! y-coordinate of pole point
INTEGER :: ipole ! Projection center flag
REAL :: di_grb ! x-direction increment or grid length
REAL :: dj_grb ! y-direction increment or grid length
REAL :: lrb ! y-direction increment or grid length
REAL :: latsw ! Latitude of South West corner point
REAL :: lonsw ! Longitude of South West corner point
REAL :: latne ! Latitude of North East corner point
REAL :: lonne ! Longitude of North East corner point
REAL :: lattru1 ! Latitude (1st) at which projection is true
REAL :: lattru2 ! Latitude (2nd) at which projection is true
REAL :: lontrue ! Longitude at which projection is true
REAL :: latrot ! Latitude of southern pole of rotation
REAL :: lonrot ! Longitude of southern pole of rotation
REAL :: angrot ! Angle of rotation
REAL :: latstr ! Latitude of the pole of stretching
REAL :: lonstr ! Longitude of the pole of stretching
REAL :: facstr ! Stretching factor
INTEGER :: scanmode ! Scanning indicator
INTEGER :: iscan ! x-direction scanning indicator
INTEGER :: jscan ! y-direction scanning indicator
INTEGER :: kscan ! FORTRAN index scanning indicator
INTEGER :: ires ! Resolution direction increments indicator
INTEGER :: iearth ! Earth shape indicator: spherical or oblate?
INTEGER :: icomp ! (u,v) components decomposition indicator
INTEGER :: jpenta ! J-Pentagonal resolution parameter
INTEGER :: kpenta ! K-Pentagonal resolution parameter
INTEGER :: mpenta ! M-Pentagonal resolution parameter
INTEGER :: ispect ! Spectral representation type
INTEGER :: icoeff ! Spectral coefficient storage mode
REAL :: xp_grb ! X coordinate of sub-satellite point
REAL :: yp_grb ! Y coordinate of sub-satellite point
REAL :: xo_grb ! X coordinate of image sector origin
REAL :: yo_grb ! Y coordinate of image sector origin
REAL :: zo_grb ! Camera altitude from center of Earth
!
!-----------------------------------------------------------------------
!
! Function f_qvsatl and inline directive for Cray PVP
!
!-----------------------------------------------------------------------
!
REAL :: f_qvsatl
!fpp$ expand (f_desdt)
!fpp$ expand (f_qvsatl)
!dir$ inline always f_desdt, f_qvsatl
!*$* inline routine (f_desdt, f_qvsatl)
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
ALLOCATE(var_grb2d(nx_ext,ny_ext,n2dvs,n2dlvt))
ALLOCATE(var_grb3d(nx_ext,ny_ext,nz_ext,n3dvs,n3dlvt))
ALLOCATE(rcdata(nx_ext*ny_ext))
ALLOCATE(var_lev3d(nz_ext,n3dvs,n3dlvt))
ALLOCATE(utmp(nx_ext,ny_ext))
ALLOCATE(vtmp(nx_ext,ny_ext))
IF ( extdfcst == ' ') extdfcst='000:00:00'
READ (extdinit,'(i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,i2)') &
iyr,imo,iday,ihr,imin,isec
CALL julday(iyr,imo,iday,jldy)
myr=MOD(iyr,100)
ifhr=0
ifmin=0
ifsec=0
READ (extdfcst,'(i3)',ERR=4,END=4) ifhr
4 CONTINUE
WRITE (gribtime,'(i4.4,i2.2,i2.2,i2.2,a1,i2.2)') &
iyr,imo,iday,ihr,'f',ifhr
len1=LEN(dir_extd)
grbflen=len1
CALL strlnth( dir_extd, grbflen )
IF( grbflen == 0 .OR. dir_extd(1:grbflen) == ' ' ) THEN
dir_extd = '.'
grbflen=1
END IF
IF( dir_extd(grbflen:grbflen) /= '/'.AND.grbflen < len1 ) THEN
grbflen=grbflen+1
dir_extd(grbflen:grbflen)='/'
END IF
lenrun = LEN( extdname )
CALL strlnth( extdname, lenrun )
gribfile = &
dir_extd(1:grbflen)//extdname(1:lenrun)//'.'//gribtime(1:13)
grbflen = grbflen + lenrun + 14
!
!-----------------------------------------------------------------------
!
! RDNMCGRB reads NMC GRIB data
!
!-----------------------------------------------------------------------
!
gridtyp = ruc211grid
mproj_grb = ruc211proj
n2dvars = ruc211nvs2d
n2dlvtps = ruc211nlvt2d
DO k=1,n2dlvtps
DO n=1,n2dvars
var_id2d(n,k) = ruc211var_id2d(n,k)
END DO
levtyp2d(k) = ruc211levs2d(k)
END DO
n3dvars = ruc211nvs3d
n3dlvtps = ruc211nlvt3d
DO m=1,n3dlvtps
DO n=1,n3dvars
var_id3d(n,m) = ruc211var_id3d(n,m)
END DO
levtyp3d(m) = ruc211levs3d(m)
END DO
CALL rdnmcgrb(nx_ext,ny_ext,nz_ext,gribfile,grbflen, gribtime, &
gridesc, iproj_grb, gthin, &
ni_grb,nj_grb,np_grb, nk_grb,zk_grb, npeq,nit, &
pi_grb,pj_grb,ipole, di_grb,dj_grb, &
latsw,lonsw, latne,lonne, &
latrot,lonrot,angrot, &
latstr,lonstr,facstr, &
lattru1,lattru2,lontrue, &
scanmode, iscan,jscan,kscan, &
ires,iearth,icomp, &
jpenta,kpenta,mpenta,ispect,icoeff, &
xp_grb,yp_grb, xo_grb,yo_grb,zo_grb, &
rcdata,var_grb2d,var_grb3d,var_lev3d,iret)
max_nr2d = 0
DO n=1,n2dvars
DO m=1,n2dlvtps
max_nr2d = MAX( max_nr2d, var_nr2d(n,m) )
END DO
END DO
max_nr3d = 0
DO n=1,n3dvars
max_nr3d = MAX( max_nr3d, var_nr3d(n,1))
END DO
IF ( max_nr3d == 0 ) THEN
WRITE (6,'(a)') &
'No 3-D variables was found in the GRIB file', &
'Program stopped in GETNMCRUC211.'
STOP
END IF
IF ( max_nr2d == 0 ) THEN
WRITE (6,'(a)') &
'No 2-D variables was found in the GRIB file'
END IF
! write (6,'(/a7,2x,6(i7))')
! : 'Lev\\VID',(var_id3d(n,1),n=1,n3dvars)
! DO 60 k=1,max_nr3d
! write (6,'(i5,4x,6(i7))')
! : k,(var_lev3d(k,n,1),n=1,n3dvars)
60 CONTINUE
DO k=1,max_nr3d
DO n=2,n3dvars
IF ( var_lev3d(k,1,1) /= var_lev3d(k,n,1) ) THEN
WRITE (6,'(a)') &
'Variables were not at the same level.', &
'Program stopped in GETNMCRUC211.'
STOP
END IF
END DO
END DO
IF ( iproj_grb == 5 .AND. ipole == 0 ) THEN ! Center in NP
iproj_ext = 1
ELSE IF ( iproj_grb == 5 .AND. ipole == 1 ) THEN ! Center in SP
iproj_ext = -1
ELSE IF ( iproj_grb == 3 .AND. ipole == 0 ) THEN ! Center in NP
iproj_ext = 2
ELSE IF ( iproj_grb == 3 .AND. ipole == 1 ) THEN ! Center in SP
iproj_ext = -2
ELSE IF ( iproj_grb == 1 ) THEN
iproj_ext = 3
ELSE IF ( iproj_grb == 0 ) THEN
iproj_ext = 4
ELSE
WRITE (6,'(a)') &
'Unknown map projection. Set to non-projection.'
iproj_ext = 0
END IF
scale_ext = 1.0
latnot_ext(1) = lattru1
latnot_ext(2) = lattru2
trlon_ext = lontrue
dx_ext = di_grb
dy_ext = dj_grb
CALL getmapr(iproj,scale,latnot,trlon,x0,y0)
CALL setmapr(iproj_ext,scale_ext,latnot_ext,trlon_ext)
CALL lltoxy(1,1,latsw,lonsw,x0_ext,y0_ext)
DO i=1,nx_ext
x_ext(i)=x0_ext+(i-1)*dx_ext
END DO
DO j=1,ny_ext
y_ext(j)=y0_ext+(j-1)*dy_ext
END DO
CALL xytoll(nx_ext,ny_ext,x_ext,y_ext,lat_ext,lon_ext)
!
!-----------------------------------------------------------------------
!
! Retrieve 2-D variables
!
!-----------------------------------------------------------------------
!
DO j=1,ny_ext
DO i=1,nx_ext
IF ( var_nr2d(1,1) == 0 ) THEN
psfc_ext (i,j) = -999.0
ELSE
psfc_ext (i,j) = var_grb2d(i,j,1,1)
END IF
IF ( var_nr2d(2,1) == 0 ) THEN
trn_ext (i,j) = -999.0
ELSE
trn_ext (i,j) = var_grb2d(i,j,2,1)
END IF
IF( var_nr2d(1,2) == 0.) THEN
t_2m_ext(i,j)= -999.0
ELSE
t_2m_ext(i,j)= var_grb2d(i,j,1,2)
END IF
! at this point we are reading in the relative humidity
! later we'll convert to specific humidity
IF( var_nr2d(2,2) == 0.) THEN
rh_2m_ext(i,j)= -999.0
ELSE
rh_2m_ext(i,j)= var_grb2d(i,j,2,2)
END IF
IF( var_nr2d(3,2) == 0.) THEN
u_10m_ext(i,j)= -999.0
ELSE
u_10m_ext(i,j)= var_grb2d(i,j,3,2)
END IF
IF( var_nr2d(4,2) == 0.) THEN
v_10m_ext(i,j)= -999.0
ELSE
v_10m_ext(i,j)= var_grb2d(i,j,4,2)
END IF
tsfc_ext (i,j) = -999.0
tdeep_ext (i,j) = -999.0
wetsfc_ext (i,j) = -999.0
wetdp_ext (i,j) = -999.0
wetcanp_ext(i,j) = -999.0
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Retrieve 3-D variables
!
!-----------------------------------------------------------------------
!
nz1 = MIN(var_nr3d(1,1),nz_ext)
IF ( var_lev3d(1,1,1) > var_lev3d(nz1,1,1) ) THEN ! 1st level at sfc
chklev = 1
lvscan = 0
ELSE
chklev = -1
lvscan = nz1+1
END IF
DO k=1,nz1
kk = chklev * k + lvscan
DO j=1,ny_ext
DO i=1,nx_ext
p_ext(i,j,kk) = 100.0 * FLOAT(var_lev3d(k,1,1)) ! Pressure
t_ext(i,j,kk) = var_grb3d(i,j,k,2,1) ! Temperature (K)
u_ext(i,j,kk) = var_grb3d(i,j,k,4,1) ! u wind (m/s)
v_ext(i,j,kk) = var_grb3d(i,j,k,5,1) ! v wind (m/s)
hgt_ext(i,j,kk) = var_grb3d(i,j,k,1,1) ! height (m)
! check for portions of constant pressure grids that are below the surface
! IF(((kk == 1) .OR. (p_ext(i,j,kk) > psfc_ext(i,j))) .AND. &
! (psfc_ext(i,j) > 0.) ) THEN
! p_ext(i,j,kk)= psfc_ext(i,j) - 1. * (kk - 1)
! t_ext(i,j,kk)= t_2m_ext(i,j)
! u_ext(i,j,kk)= u_10m_ext(i,j)
! v_ext(i,j,kk)= v_10m_ext(i,j)
! hgt_ext(i,j,kk)= trn_ext(i,j) + 0.1 * (kk - 1)
! END IF
IF( (p_ext(i,j,kk) > 0.0) .AND. (t_ext(i,j,kk) > 0.0) ) THEN
qvsat = f_qvsatl( p_ext(i,j,kk), t_ext(i,j,kk) )
qv_ext(i,j,kk)= var_grb3d(i,j,k,3,1) * qvsat * 0.01
IF((kk == 1) .OR. (p_ext(i,j,kk) > psfc_ext(i,j)))THEN
qv_ext(i,j,kk)= rh_2m_ext(i,j) * qvsat * 0.01
END IF
qc_ext(i,j,kk)= 0.0
qi_ext(i,j,kk)= 0.0
ELSE
qv_ext(i,j,kk)= 0.0
qi_ext(i,j,kk)= 0.0
qc_ext(i,j,kk)= 0.0
END IF
qr_ext (i,j,kk) = -999.
qs_ext (i,j,kk) = -999.
qh_ext (i,j,kk) = -999.
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Rotate winds to be relative to true north.
! The RUCawips data are sent as grid-relative.
!
!-----------------------------------------------------------------------
!
DO k=1,nz1
CALL uvmptoe(nx_ext,ny_ext,u_ext(1,1,k),v_ext(1,1,k), &
lon_ext,utmp,vtmp)
u_ext(:,:,k) = utmp(:,:)
v_ext(:,:,k) = vtmp(:,:)
END DO
!
!-----------------------------------------------------------------------
!
! Reset map projection to previous values
!
!-----------------------------------------------------------------------
!
CALL setmapr(iproj,scale,latnot,trlon)
CALL setorig(1,x0,y0)
istatus = 1
DEALLOCATE(var_grb2d,var_grb3d,rcdata,var_lev3d)
DEALLOCATE(utmp)
DEALLOCATE(vtmp)
RETURN
END SUBROUTINE getnmcruc211
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE GETREANALT62 ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE getreanalt62(nx_ext,ny_ext,nz_ext, & 1,4
dir_extd,extdname,extdopt,extdfmt, &
extdinit,extdfcst,julfname, &
iproj_ext,scale_ext, &
trlon_ext,latnot_ext,x0_ext,y0_ext, &
lat_ext,lon_ext, &
p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext, &
qc_ext,qr_ext,qi_ext,qs_ext,qh_ext, &
tsfc_ext,tdeep_ext,wetsfc_ext,wetdp_ext,wetcanp_ext, &
trn_ext,psfc_ext, istatus)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
!
! Reads a Global ReAnalaysis file for processing by ext2arps,
! a program that converts external files to ARPS variables and format.
! This version is useful when you want to use an ARPS file
! with a different orientation or terrain than your final
! ARPS product so arpsr2h does not work.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Dennis Moon, SSESCO
! 06/19/1998
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! dir_extd Directory name for external file
! extdname Prefix string of external file name
! extdopt Option of external data sources
! extdfmt Option of external data format
!
! extdinit Initialized time in mm-dd-yyyy:hh:mm:ss format
! extdfcst Forecast hour in HHH:MM:SS format
! julfname File name in yyjjjhhmm format
!
! OUTPUT:
!
! iproj_ext Map projection number of external data
! scale_ext Scale factor of external data
! trlon_ext True longitude of external data (degrees E)
! latnot_ext(2) True latitude(s) of external data (degrees N)
! x0_ext x coordinate of origin of external data
! y0_ext y coordinate of origin of external data
! lat_ext latitude of external data points (degrees N)
! lon_ext longitude of external data points (degrees E)
! p_ext pressure (Pascal)
! hgt_ext height (m)
! t_ext temperature (K)
! qv_ext specific humidity (kg/kg)
! u_ext u wind component (m/s)
! v_ext v wind component (m/s)
! qc_ext Cloud water mixing ratio (kg/kg)
! qr_ext Rain water mixing ratio (kg/kg)
! qi_ext Ice mixing ratio (kg/kg)
! qs_ext Snow mixing ratio (kg/kg)
! qh_ext Hail mixing ratio (kg/kg)
!
! tsfc_ext Surface temperature
! tdeep_ext Soil temperature
! wetsfc_ext Top layer soil moisture (fraction)
! wetdp_ext Deep soil moisture (fraction)
! wetcanp_ext Water content on canopy
!
! trn_ext External terrain (m)
! psfc_ext Surface pressure (Pa)
!
! istatus status indicator
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INCLUDE 'gribcst.inc'
CHARACTER (LEN=*) :: dir_extd
CHARACTER (LEN=*) :: extdname
INTEGER :: extdopt
INTEGER :: extdfmt
CHARACTER (LEN=19) :: extdinit
CHARACTER (LEN=9) :: extdfcst
CHARACTER (LEN=9) :: julfname
!
!-----------------------------------------------------------------------
!
! External grid variables
!
!-----------------------------------------------------------------------
!
INTEGER :: iproj_ext
REAL :: scale_ext,trlon_ext
REAL :: latnot_ext(2)
REAL :: x0_ext,y0_ext
REAL :: dx_ext,dy_ext
!
!-----------------------------------------------------------------------
!
! Output external variable arrays
!
!-----------------------------------------------------------------------
!
INTEGER :: nx_ext,ny_ext,nz_ext
REAL :: lat_ext(nx_ext,ny_ext)
REAL :: lon_ext(nx_ext,ny_ext)
REAL :: p_ext (nx_ext,ny_ext,nz_ext) ! Pressure (Pascals)
REAL :: hgt_ext(nx_ext,ny_ext,nz_ext) ! Height (m)
REAL :: t_ext (nx_ext,ny_ext,nz_ext) ! Temperature (K)
REAL :: qv_ext (nx_ext,ny_ext,nz_ext) ! Specific humidity (kg/kg)
REAL :: u_ext (nx_ext,ny_ext,nz_ext) ! Eastward wind component
REAL :: v_ext (nx_ext,ny_ext,nz_ext) ! Northward wind component
REAL :: qc_ext (nx_ext,ny_ext,nz_ext) ! Cloud H2O mixing ratio (kg/kg)
REAL :: qr_ext (nx_ext,ny_ext,nz_ext) ! Rain H2O mixing ratio (kg/kg)
REAL :: qi_ext (nx_ext,ny_ext,nz_ext) ! Ice mixing ratio (kg/kg)
REAL :: qs_ext (nx_ext,ny_ext,nz_ext) ! Snow mixing ratio (kg/kg)
REAL :: qh_ext (nx_ext,ny_ext,nz_ext) ! Hail mixing ratio (kg/kg)
REAL :: tsfc_ext (nx_ext,ny_ext) ! Temperature at surface (K)
REAL :: tdeep_ext (nx_ext,ny_ext) ! Deep soil temperature (K)
REAL :: wetsfc_ext (nx_ext,ny_ext) ! Surface soil moisture (fraction)
REAL :: wetdp_ext (nx_ext,ny_ext) ! Deep soil moisture (fraction)
REAL :: wetcanp_ext(nx_ext,ny_ext) ! Canopy water amount
REAL :: trn_ext (nx_ext,ny_ext) ! Geometrical heights
REAL :: psfc_ext (nx_ext,ny_ext) ! Surface pressure (Pa)
!
!----------------------------------------------------------------------
!
! Other external variable arrays
!
!----------------------------------------------------------------------
!
REAL :: x_ext(nx_ext)
REAL :: y_ext(ny_ext)
REAL :: z_ext(nz_ext)
INTEGER :: istatus
!
!-----------------------------------------------------------------------
!
! Work arrays for storing grib data
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: var_grb2d(:,:,:,:) ! GRIB variables
REAL, ALLOCATABLE :: var_grb3d(:,:,:,:,:) ! GRIB 3-D variables
INTEGER, ALLOCATABLE :: var_lev3d(:,:,:) ! Levels (hybrid) for
! each 3-D variable
REAL, ALLOCATABLE :: rcdata(:) ! temporary data array
!
!-----------------------------------------------------------------------
!
! Original grid variables
!
!-----------------------------------------------------------------------
!
INTEGER :: iproj
REAL :: scale,trlon,x0,y0
REAL :: latnot(2)
!
!-----------------------------------------------------------------------
!
! Misc internal variables
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=80) :: gribfile
CHARACTER (LEN=13) :: gribtime
CHARACTER (LEN=10) :: runstr
CHARACTER (LEN=3) :: fmtn
INTEGER :: ichr,bchar,echar
INTEGER :: i,j,k,ldir,ireturn
INTEGER :: iyr,imo,iday,myr, jldy
INTEGER :: ihr,imin,isec
INTEGER :: ifhr,ifmin,ifsec
INTEGER :: grbflen, len1, lenrun
INTEGER :: m,n,nz1,max_nr2d,max_nr3d
REAL :: qvsat, pilev, qvsatice
REAL :: tema, temb, qvavg, mixrat, tvirtbar, tmpval
REAL :: gausslat(200)
INTEGER :: chklev, lvscan, kk, jj
INTEGER :: iret ! Return flag
!
!-----------------------------------------------------------------------
!
! Include files
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
! GRIB grid information
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=42) :: gridesc ! Grid description
INTEGER :: iproj_grb ! Map projection indicator
INTEGER :: gthin ! Indicator of whether the grid is "thinned"
INTEGER :: ni_grb ! Number of points along x-axis
INTEGER :: nj_grb ! Number of points along y-axis
INTEGER :: np_grb ! Total number of horizontal grid points
INTEGER :: nk_grb ! Number of vertical parameters
REAL :: zk_grb(nz_ext) ! Vertical coordinate parameters
INTEGER :: npeq ! Number of lat circles from pole to equator
INTEGER :: nit(nz_ext) ! Number of x-points for thinned grid
REAL :: pi_grb ! x-coordinate of pole point
REAL :: pj_grb ! y-coordinate of pole point
INTEGER :: ipole ! Projection center flag
REAL :: di_grb ! x-direction increment or grid length
REAL :: dj_grb ! y-direction increment or grid length
REAL :: lrb ! y-direction increment or grid length
REAL :: latsw ! Latitude of South West corner point
REAL :: lonsw ! Longitude of South West corner point
REAL :: latne ! Latitude of North East corner point
REAL :: lonne ! Longitude of North East corner point
REAL :: lattru1 ! Latitude (1st) at which projection is true
REAL :: lattru2 ! Latitude (2nd) at which projection is true
REAL :: lontrue ! Longitude at which projection is true
REAL :: latrot ! Latitude of southern pole of rotation
REAL :: lonrot ! Longitude of southern pole of rotation
REAL :: angrot ! Angle of rotation
REAL :: latstr ! Latitude of the pole of stretching
REAL :: lonstr ! Longitude of the pole of stretching
REAL :: facstr ! Stretching factor
INTEGER :: scanmode ! Scanning indicator
INTEGER :: iscan ! x-direction scanning indicator
INTEGER :: jscan ! y-direction scanning indicator
INTEGER :: kscan ! FORTRAN index scanning indicator
INTEGER :: ires ! Resolution direction increments indicator
INTEGER :: iearth ! Earth shape indicator: spherical or oblate?
INTEGER :: icomp ! (u,v) components decomposition indicator
INTEGER :: jpenta ! J-Pentagonal resolution parameter
INTEGER :: kpenta ! K-Pentagonal resolution parameter
INTEGER :: mpenta ! M-Pentagonal resolution parameter
INTEGER :: ispect ! Spectral representation type
INTEGER :: icoeff ! Spectral coefficient storage mode
REAL :: xp_grb ! X coordinate of sub-satellite point
REAL :: yp_grb ! Y coordinate of sub-satellite point
REAL :: xo_grb ! X coordinate of image sector origin
REAL :: yo_grb ! Y coordinate of image sector origin
REAL :: zo_grb ! Camera altitude from center of Earth
!
!-----------------------------------------------------------------------
!
! T62 Gaussian grid latitudes
!
!-----------------------------------------------------------------------
!
DATA gausslat / 88.5420, 86.6532, 84.7532, 82.8508, 80.9474, &
79.0435, 77.1393, 75.2351, 73.3307, 71.4262, &
69.5217, 67.6171, 65.7125, 63.8079, 61.9033, &
59.9986, 58.0940, 56.1893, 54.2846, 52.3799, &
50.4752, 48.5705, 46.6658, 44.7611, 42.8564, &
40.9517, 39.0470, 37.1422, 35.2375, 33.3328, &
31.4281, 29.5234, 27.6186, 25.7139, 23.8092, &
21.9044, 19.9997, 18.0950, 16.1902, 14.2855, &
12.3808, 10.4760, 8.5713, 6.6666, 4.7618, &
2.8571, 0.9524, -0.9524, -2.8571, -4.7618, &
-6.6666, -8.5713, -10.4760, -12.3808, -14.2855, &
-16.1902, -18.0950, -19.9997, -21.9044, -23.8092, &
-25.7139, -27.6186, -29.5234, -31.4281, -33.3328, &
-35.2375, -37.1422, -39.0470, -40.9517, -42.8564, &
-44.7611, -46.6658, -48.5705, -50.4752, -52.3799, &
-54.2846, -56.1893, -58.0940, -59.9986, -61.9033, &
-63.8079, -65.7125, -67.6171, -69.5217, -71.4262, &
-73.3307, -75.2351, -77.1393, -79.0435, -80.9474, &
-82.8508, -84.7532, -86.6532, -88.5420, 106*0.0/
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
ALLOCATE(var_grb2d(nx_ext,ny_ext,n2dvs,n2dlvt))
ALLOCATE(var_grb3d(nx_ext,ny_ext,nz_ext,n3dvs,n3dlvt))
ALLOCATE(rcdata(nx_ext*ny_ext))
ALLOCATE(var_lev3d(nz_ext,n3dvs,n3dlvt))
IF ( extdfcst == ' ') extdfcst='000:00:00'
READ (extdinit,'(i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,i2)') &
iyr,imo,iday,ihr,imin,isec
CALL julday(iyr,imo,iday,jldy)
myr=MOD(iyr,100)
ifhr=0
ifmin=0
ifsec=0
READ (extdfcst,'(i3)',ERR=4,END=4) ifhr
4 CONTINUE
WRITE (gribtime,'(i4.4,i2.2,i2.2,i2.2,a1,i2.2)') &
iyr,imo,iday,ihr,'f',ifhr
len1=LEN(dir_extd)
grbflen=len1
CALL strlnth( dir_extd, grbflen )
IF( grbflen == 0 .OR. dir_extd(1:grbflen) == ' ' ) THEN
dir_extd = '.'
grbflen=1
END IF
IF( dir_extd(grbflen:grbflen) /= '/'.AND.grbflen < len1 ) THEN
grbflen=grbflen+1
dir_extd(grbflen:grbflen)='/'
END IF
lenrun = LEN( extdname )
CALL strlnth( extdname, lenrun )
gribfile = &
dir_extd(1:grbflen)//extdname(1:lenrun)//'.'//gribtime(1:13)
grbflen = grbflen + lenrun + 14
!
!-----------------------------------------------------------------------
!
! RDNMCGRB reads NMC GRIB data
!
!-----------------------------------------------------------------------
!
gridtyp = reanalt62grid
mproj_grb = reanalt62proj
n2dvars = reanalt62nvs2d
n2dlvtps = reanalt62nlvt2d
DO k=1,n2dlvtps
DO n=1,n2dvars
var_id2d(n,k) = reanalt62var_id2d(n,k)
END DO
levtyp2d(k) = reanalt62levs2d(k)
END DO
n3dvars = reanalt62nvs3d
n3dlvtps = reanalt62nlvt3d
DO m=1,n3dlvtps
DO n=1,n3dvars
var_id3d(n,m) = reanalt62var_id3d(n,m)
END DO
levtyp3d(m) = reanalt62levs3d(m)
END DO
CALL rdnmcgrb(nx_ext,ny_ext,nz_ext,gribfile,grbflen, gribtime, &
gridesc, iproj_grb, gthin, &
ni_grb,nj_grb,np_grb, nk_grb,zk_grb, npeq,nit, &
pi_grb,pj_grb,ipole, di_grb,dj_grb, &
latsw,lonsw, latne,lonne, &
latrot,lonrot,angrot, &
latstr,lonstr,facstr, &
lattru1,lattru2,lontrue, &
scanmode, iscan,jscan,kscan, &
ires,iearth,icomp, &
jpenta,kpenta,mpenta,ispect,icoeff, &
xp_grb,yp_grb, xo_grb,yo_grb,zo_grb, &
rcdata,var_grb2d,var_grb3d,var_lev3d,iret)
max_nr2d = 0
DO n=1,n2dvars
DO m=1,n2dlvtps
max_nr2d = MAX( max_nr2d, var_nr2d(n,m) )
END DO
END DO
max_nr3d = 0
DO n=1,n3dvars
max_nr3d = MAX( max_nr3d, var_nr3d(n,1))
END DO
IF ( max_nr3d == 0 ) THEN
WRITE (6,'(a)') &
'No 3-D variables was found in the GRIB file', &
'Program stopped in GETREANALT62.'
STOP
END IF
IF ( max_nr2d == 0 ) THEN
WRITE (6,'(a)') &
'No 2-D variables was found in the GRIB file'
END IF
! write (6,'(/a7,2x,6(i7))')
! : 'Lev\\VID',(var_id3d(n,1),n=1,n3dvars)
! DO 60 k=1,max_nr3d
! write (6,'(i5,4x,6(i7))')
! : k,(var_lev3d(k,n,1),n=1,n3dvars)
60 CONTINUE
DO k=1,max_nr3d
DO n=2,n3dvars
IF ( var_lev3d(k,1,1) /= var_lev3d(k,n,1) ) THEN
WRITE (6,'(a)') &
'Variables were not at the same level.', &
'Program stopped in GETREANALT62.'
STOP
END IF
END DO
END DO
IF ( iproj_grb == 5 .AND. ipole == 0 ) THEN ! Center in NP
iproj_ext = 1
ELSE IF ( iproj_grb == 5 .AND. ipole == 1 ) THEN ! Center in SP
iproj_ext = -1
ELSE IF ( iproj_grb == 3 .AND. ipole == 0 ) THEN ! Center in NP
iproj_ext = 2
ELSE IF ( iproj_grb == 3 .AND. ipole == 1 ) THEN ! Center in SP
iproj_ext = -2
ELSE IF ( iproj_grb == 1 ) THEN
iproj_ext = 3
ELSE IF ( (iproj_grb == 0) .OR. (iproj_grb == 4) ) THEN
iproj_ext = 4
ELSE
WRITE (6,'(a)') &
'Unknown map projection. Set to non-projection.'
iproj_ext = 0
END IF
scale_ext = 1.0
latnot_ext(1) = lattru1
latnot_ext(2) = lattru2
trlon_ext = lontrue
latnot_ext(1)= 0.
trlon_ext= 0.
dx_ext = di_grb
dy_ext = dj_grb
DO i=1, nx_ext
DO j=1, ny_ext
lat_ext(i,j)= gausslat(j)
lon_ext(i,j)= lonsw + (i-1) * dx_ext
END DO
END DO
! swap the latitude values so that +j moves to the north
DO i=1, nx_ext
DO j=1, ny_ext/2
tmpval= lat_ext(i,ny_ext-j+1)
lat_ext(i,ny_ext-j+1)= lat_ext(i,j)
lat_ext(i,j)= tmpval
END DO
END DO
! call getmapr(iproj,scale,latnot,trlon,x0,y0)
! call setmapr(iproj_ext,scale_ext,latnot_ext,trlon_ext)
! call lltoxy(1,1,LatSW,LonSW,x0_ext,y0_ext)
! DO 100 i=1,nx_ext
! x_ext(i)=x0_ext+(i-1)*dx_ext
!100 CONTINUE
! DO 110 j=1,ny_ext
! y_ext(j)=y0_ext+(j-1)*dy_ext
!110 CONTINUE
! CALL xytoll(nx_ext,ny_ext,x_ext,y_ext,lat_ext,lon_ext)
!
!-----------------------------------------------------------------------
!
! Retrieve 2-D variables
!
!-----------------------------------------------------------------------
!
DO j=1,ny_ext
DO i=1,nx_ext
IF ( var_nr2d(1,1) == 0 ) THEN
psfc_ext (i,j) = -999.0
ELSE
psfc_ext (i,j) = var_grb2d(i,j,1,1)
END IF
IF ( var_nr2d(2,1) == 0 ) THEN
trn_ext (i,j) = -999.0
ELSE
trn_ext (i,j) = var_grb2d(i,j,2,1)
END IF
tsfc_ext (i,j) = -999.0
tdeep_ext (i,j) = -999.0
wetsfc_ext (i,j) = -999.0
wetdp_ext (i,j) = -999.0
wetcanp_ext(i,j) = -999.0
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Retrieve 3-D variables
!
!-----------------------------------------------------------------------
!
nz1 = MIN(var_nr3d(1,1),nz_ext)
IF ( var_lev3d(1,1,1) > var_lev3d(nz1,1,1) ) THEN ! 1st level at sfc
chklev = 1
lvscan = 0
ELSE
chklev = -1
lvscan = nz1+1
END IF
DO k=1,nz1
kk = chklev * k + lvscan
DO j=1,ny_ext
DO i=1,nx_ext
p_ext(i,j,kk) = 1.e-4 * FLOAT(var_lev3d(k,1,1)) &
* psfc_ext(i,j) ! Pressure
t_ext(i,j,kk) = var_grb3d(i,j,k,1,1) ! Temperature (K)
qv_ext(i,j,kk)= var_grb3d(i,j,k,2,1) ! Specific Humidity
u_ext(i,j,kk) = var_grb3d(i,j,k,3,1) ! u wind (m/s)
v_ext(i,j,kk) = var_grb3d(i,j,k,4,1) ! v wind (m/s)
! calculate height from sigma-P values
IF( k == 1) THEN
mixrat= 1.0 / ( 1./qv_ext(i,j,kk) - 1.)
tvirtbar= t_ext(i,j,kk) * ( 1.0 + .61* mixrat)
hgt_ext(i,j,kk)= trn_ext(i,j) + 29.3 * tvirtbar * &
LOG(psfc_ext(i,j)/p_ext(i,j,kk))
ELSE
qvavg= 0.5 * (qv_ext(i,j,kk) + qv_ext(i,j,kk-1))
mixrat= 1.0 / (1./qvavg - .1)
tvirtbar= 0.5 * (t_ext(i,j,kk) + t_ext(i,j,kk-1)) * &
( 1.0 + .61 * mixrat)
hgt_ext(i,j,kk)= hgt_ext(i,j,kk-1) + 29.3 * tvirtbar * &
LOG(p_ext(i,j,kk-1)/p_ext(i,j,kk))
END IF
qc_ext(i,j,kk)= 0.0
qi_ext(i,j,kk)= 0.0
qr_ext (i,j,kk) = -999.
qs_ext (i,j,kk) = -999.
qh_ext (i,j,kk) = -999.
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! No need to rotate winds to be relative to true north.
! The Re-analysis grid winds are true N-S, E-W
!
!-----------------------------------------------------------------------
!
! DO 300 k=1,nz1
! CALL uvmptoe(nx_ext,ny_ext,u_ext(1,1,k),v_ext(1,1,k),
! : lon_ext,u_ext(1,1,k),v_ext(1,1,k))
!300 CONTINUE
!
!-----------------------------------------------------------------------
!
! Reset map projection to previous values
!
!-----------------------------------------------------------------------
!
! call setmapr(iproj,scale,latnot,trlon)
! call setorig(1,x0,y0)
istatus = 1
DEALLOCATE(var_grb2d,var_grb3d,rcdata,var_lev3d)
RETURN
END SUBROUTINE getreanalt62
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE TV2TQ ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE tv2tq(tv,p,t,q) 2,4
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
REAL :: tv ! virtual temperature (K)
REAL :: p ! pressure (pa)
REAL :: t ! temperature (K)
REAL :: q ! specific humidity (kg/kg)
REAL :: tg,qg,fac,tvg,dtvdt,tgnew
REAL :: es, es1
INTEGER :: it
!
!-----------------------------------------------------------------------
!
! Function and inline directive for Cray PVP
!
!-----------------------------------------------------------------------
!
REAL :: f_es
!fpp$ expand (f_es)
!dir$ inline always f_es
!*$* inline routine (f_es)
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
it = 0
!
! -- TG - guess at non-virtual temperature
!
es = f_es(p,tv)
tg = tv/( 1.0 + 0.6078 * 0.62197*es &
/ ( p - es ) )
10 CONTINUE
it = it + 1
!
! -- QG - guess at specific humidity
!
es = f_es(p,tg)
qg = 0.62197*es/( p - es )
fac = 1.+0.6078*qg
!
! -- TVG - guess at virtual temperature
!
tvg = tg*fac
!
! -- DTVDT - d(Tv)/dT estimated over 0.1 K interval
! from 1st term Taylor series expansion
!
es1 = f_es(p,tg+0.1)
dtvdt = fac + tg * 0.6078 * 0.62197*p/(p-es)**2 &
* ( es1 - es ) / 0.1
tgnew = tg + (tv-tvg)/dtvdt
!
! write(12,*) ' '
! write(12,*) 'IT=',IT,' TV=',TV,' TVG=',TVG,
! : ' TV-TVG=',TV-TVG,' DTVDT=',DTVDT
!
IF (ABS(tv-tvg) < 1.0E-4) GO TO 20
tg = tgnew
GO TO 10
20 CONTINUE
t = tgnew
es = f_es(p,t)
q = 0.62197*es / ( p - es )
RETURN
END SUBROUTINE tv2tq
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE GETNMCETA212 ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE getnmceta212(nx_ext,ny_ext,nz_ext,nzsoil_ext, & 1,11
nstyps_ext,dir_extd,extdname,extdopt,extdfmt, &
extdinit,extdfcst,julfname, &
iproj_ext,scale_ext, &
trlon_ext,latnot_ext,x0_ext,y0_ext, &
lat_ext,lon_ext,zpsoil_ext, &
p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext, &
qc_ext,qr_ext,qi_ext,qs_ext,qh_ext, &
tsoil_ext,qsoil_ext,wetcanp_ext, &
snowdpth_ext,trn_ext,psfc_ext,soiltyp_ext, &
istatus)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Reads a NMC GRIB (Grid #212, 40km) ETA file for processing by
! ext2arps, a program that converts external files to ARPS variables
! and format.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Yuhe Liu
! 02/22/1996
!
! MODIFICATION HISTORY:
!
! 1999/08/03 (Gene Bassett)
! Modified the deep soil moisture and temperature to be an average of
! the three soil layers covering 0 to 100 cm in depth.
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! dir_extd Directory name for external file
! extdname Prefix string of external file name
! extdopt Option of external data sources
! extdfmt Option of external data format
! extdinit Initialized time in mm-dd-yyyy:hh:mm:ss format
! extdfcst Forecast hour in HHH:MM:SS format
! julfname File name in yyjjjhhmm format
!
! OUTPUT:
!
! iproj_ext Map projection number of external data
! scale_ext Scale factor of external data
! trlon_ext True longitude of external data (degrees E)
! latnot_ext(2) True latitude(s) of external data (degrees N)
! x0_ext x coordinate of origin of external data
! y0_ext y coordinate of origin of external data
! lat_ext latitude of external data points (degrees N)
! lon_ext longitude of external data points (degrees E)
! p_ext pressure (Pascal)
! hgt_ext height (m)
! t_ext temperature (K)
! qv_ext specific humidity (kg/kg)
! u_ext u wind component (m/s)
! v_ext v wind component (m/s)
! qc_ext Cloud water mixing ratio (kg/kg)
! qr_ext Rain water mixing ratio (kg/kg)
! qi_ext Ice mixing ratio (kg/kg)
! qs_ext Snow mixing ratio (kg/kg)
! qh_ext Hail mixing ratio (kg/kg)
!
! tsoil_ext Soil temperature (K)
! qsoil_ext Soil moisture (m**3/m**3)
! wetcanp_ext Water content on canopy
!
! trn_ext External terrain (m)
! psfc_ext Surface pressure (Pa)
!
! istatus status indicator
!
! WORK ARRAYS:
!
! var_grb3d Arrays to store the GRIB 3-D variables:
! var_grb3d(nxgrb,nygrb,nzgrb,1,1) - Temperature (K)
! var_grb3d(nxgrb,nygrb,nzgrb,2,1) - Specific humidity
! (kg/kg)
! var_grb3d(nxgrb,nygrb,nzgrb,3,1) - u wind (m/s)
! var_grb3d(nxgrb,nygrb,nzgrb,4,1) - v wind (m/s)
! var_grb3d(nxgrb,nygrb,nzgrb,5,1) - Geopotential
! height (gpm)
! var_grb3d(nxgrb,nygrb,nzgrb,6,1) - Pressure vertical
! velocity (Pa/s)
! (if applied)
! var_grb3d(nxgrb,nygrb,nzgrb,1,2) - soil temp. (K)
! var_grb3d(nxgrb,nygrb,nzgrb,2,2) - vol. soil moist.
! (m**3/m**3)
!
! var_grb2d Arrays to store the GRIB 2-D variables:
! var_grb2d(nxgrb,nygrb,1) - Surface pressure (Pa)
! var_grb2d(nxgrb,nygrb,2) - Geopotential height (gpm)
! var_grb2d(nxgrb,nygrb,3) - Surface temperature (K)
! var_grb2d(nxgrb,nygrb,4) - Plant canopy surface
! water (kg/m**2)
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INCLUDE 'gribcst.inc'
CHARACTER (LEN=*) :: dir_extd
CHARACTER (LEN=*) :: extdname
INTEGER :: extdopt
INTEGER :: extdfmt
CHARACTER (LEN=19) :: extdinit
CHARACTER (LEN=9) :: extdfcst
CHARACTER (LEN=9) :: julfname
!
!-----------------------------------------------------------------------
!
! External grid variables
!
!-----------------------------------------------------------------------
!
INTEGER :: iproj_ext
REAL :: scale_ext,trlon_ext
REAL :: latnot_ext(2)
REAL :: x0_ext,y0_ext
REAL :: dx_ext,dy_ext
!
!-----------------------------------------------------------------------
!
! Output external variable arrays
!
!-----------------------------------------------------------------------
!
INTEGER :: nx_ext,ny_ext,nz_ext
INTEGER :: nzsoil_ext, nstyps_ext
REAL :: lat_ext(nx_ext,ny_ext)
REAL :: lon_ext(nx_ext,ny_ext)
REAL :: p_ext (nx_ext,ny_ext,nz_ext) ! Pressure (Pascals)
REAL :: hgt_ext(nx_ext,ny_ext,nz_ext) ! Height (m)
REAL :: t_ext (nx_ext,ny_ext,nz_ext) ! Temperature (K)
REAL :: qv_ext (nx_ext,ny_ext,nz_ext) ! Specific humidity (kg/kg)
REAL :: u_ext (nx_ext,ny_ext,nz_ext) ! Eastward wind component
REAL :: v_ext (nx_ext,ny_ext,nz_ext) ! Northward wind component
REAL :: qc_ext (nx_ext,ny_ext,nz_ext) ! Cloud H2O mixing ratio (kg/kg)
REAL :: qr_ext (nx_ext,ny_ext,nz_ext) ! Rain H2O mixing ratio (kg/kg)
REAL :: qi_ext (nx_ext,ny_ext,nz_ext) ! Ice mixing ratio (kg/kg)
REAL :: qs_ext (nx_ext,ny_ext,nz_ext) ! Snow mixing ratio (kg/kg)
REAL :: qh_ext (nx_ext,ny_ext,nz_ext) ! Hail mixing ratio (kg/kg)
REAL :: zpsoil_ext(nx_ext,ny_ext,nzsoil_ext) !Soil depths (m)
REAL :: tsoil_ext (nx_ext,ny_ext,nzsoil_ext) ! Soil temperature (K)
REAL :: qsoil_ext (nx_ext,ny_ext,nzsoil_ext) ! Soil moisture (m**3/m**3)
REAL :: wetcanp_ext(nx_ext,ny_ext) ! Canopy water amount
REAL :: snowdpth_ext(nx_ext,ny_ext) ! Snow depth (m)
REAL :: trn_ext (nx_ext,ny_ext) ! External terrain (m)
REAL :: psfc_ext (nx_ext,ny_ext) ! Surface pressure (Pa)
! INTEGER soiltyp_ext (nx_ext,ny_ext,nstyps_ext) ! Soil type
INTEGER soiltyp_ext (nx_ext,ny_ext) ! Soil type
!
!-----------------------------------------------------------------------
!
! Other external variable arrays
!
!-----------------------------------------------------------------------
!
REAL :: x_ext(nx_ext)
REAL :: y_ext(ny_ext)
INTEGER :: istatus
!
!-----------------------------------------------------------------------
!
! Work arrays for storing grib data
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: var_grb2d(:,:,:,:) ! GRIB variables
REAL, ALLOCATABLE :: var_grb3d(:,:,:,:,:) ! GRIB 3-D variables
INTEGER, ALLOCATABLE :: var_lev3d(:,:,:) ! Levels (hybrid) for
! each 3-D variable
REAL, ALLOCATABLE :: rcdata(:) ! temporary data array
!
!-----------------------------------------------------------------------
!
! Original grid variables
!
!-----------------------------------------------------------------------
!
INTEGER :: iproj
REAL :: scale,trlon,x0,y0
REAL :: latnot(2)
!
!-----------------------------------------------------------------------
!
! Misc internal variables
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=80) :: gribfile
CHARACTER (LEN=13) :: gribtime
INTEGER :: i,j,k,kk
INTEGER :: iyr,imo,iday,myr, jldy
INTEGER :: ihr,imin,isec
INTEGER :: ifhr,ifmin,ifsec
INTEGER :: grbflen, len1, lenrun
INTEGER :: m,n,nz1,max_nr2d,max_nr3d
REAL :: govrd
! REAL :: soildepth_ext(4)
REAL :: soildepth_ext(5) ! EMK 15 June 2002
INTEGER :: chklev, lvscan
INTEGER :: iret ! Return flag
REAL, ALLOCATABLE :: utmp(:,:), vtmp(:,:)
!
!-----------------------------------------------------------------------
!
! Include files
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'phycst.inc'
INCLUDE 'grid.inc'
!
!-----------------------------------------------------------------------
!
! GRIB grid information
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=42) :: gridesc ! Grid description
INTEGER :: iproj_grb ! Map projection indicator
INTEGER :: gthin ! Indicator of whether the grid is "thinned"
INTEGER :: ni_grb ! Number of points along x-axis
INTEGER :: nj_grb ! Number of points along y-axis
INTEGER :: np_grb ! Total number of horizontal grid points
INTEGER :: nk_grb ! Number of vertical parameters
REAL :: zk_grb(nz_ext) ! Vertical coordinate parameters
INTEGER :: npeq ! Number of lat circles from pole to equator
INTEGER :: nit(nz_ext) ! Number of x-points for thinned grid
REAL :: pi_grb ! x-coordinate of pole point
REAL :: pj_grb ! y-coordinate of pole point
INTEGER :: ipole ! Projection center flag
REAL :: di_grb ! x-direction increment or grid length
REAL :: dj_grb ! y-direction increment or grid length
REAL :: latsw ! Latitude of South West corner point
REAL :: lonsw ! Longitude of South West corner point
REAL :: latne ! Latitude of North East corner point
REAL :: lonne ! Longitude of North East corner point
REAL :: lattru1 ! Latitude (1st) at which projection is true
REAL :: lattru2 ! Latitude (2nd) at which projection is true
REAL :: lontrue ! Longitude at which projection is true
REAL :: latrot ! Latitude of southern pole of rotation
REAL :: lonrot ! Longitude of southern pole of rotation
REAL :: angrot ! Angle of rotation
REAL :: latstr ! Latitude of the pole of stretching
REAL :: lonstr ! Longitude of the pole of stretching
REAL :: facstr ! Stretching factor
INTEGER :: scanmode ! Scanning indicator
INTEGER :: iscan ! x-direction scanning indicator
INTEGER :: jscan ! y-direction scanning indicator
INTEGER :: kscan ! FORTRAN index scanning indicator
INTEGER :: ires ! Resolution direction increments indicator
INTEGER :: iearth ! Earth shape indicator: spherical or oblate?
INTEGER :: icomp ! (u,v) components decomposition indicator
INTEGER :: jpenta ! J-Pentagonal resolution parameter
INTEGER :: kpenta ! K-Pentagonal resolution parameter
INTEGER :: mpenta ! M-Pentagonal resolution parameter
INTEGER :: ispect ! Spectral representation type
INTEGER :: icoeff ! Spectral coefficient storage mode
REAL :: xp_grb ! X coordinate of sub-satellite point
REAL :: yp_grb ! Y coordinate of sub-satellite point
REAL :: xo_grb ! X coordinate of image sector origin
REAL :: yo_grb ! Y coordinate of image sector origin
REAL :: zo_grb ! Camera altitude from center of Earth
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
ALLOCATE(var_grb2d(nx_ext,ny_ext,n2dvs,n2dlvt))
ALLOCATE(var_grb3d(nx_ext,ny_ext,nz_ext,n3dvs,n3dlvt))
ALLOCATE(rcdata(nx_ext*ny_ext))
ALLOCATE(var_lev3d(nz_ext,n3dvs,n3dlvt))
ALLOCATE(utmp(nx_ext,ny_ext))
ALLOCATE(vtmp(nx_ext,ny_ext))
! DATA soildepth_ext/0.0,-0.05,-0.25,-0.70/
DATA soildepth_ext/0.0,-0.05,-0.25,-0.70,-1.50/ ! EMK 15 June 2002
IF(extdfcst == ' ') extdfcst='000:00:00'
READ (extdinit,'(i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,i2)') &
iyr,imo,iday,ihr,imin,isec
CALL julday(iyr,imo,iday,jldy)
myr=MOD(iyr,100)
ifhr=0
ifmin=0
ifsec=0
READ(extdfcst,'(i3)',ERR=4,END=4) ifhr
4 CONTINUE
WRITE (gribtime,'(i4.4,i2.2,i2.2,i2.2,a1,i2.2)') &
iyr,imo,iday,ihr,'f',ifhr
len1=LEN(dir_extd)
grbflen=len1
CALL strlnth( dir_extd, grbflen )
IF( grbflen == 0 .OR. dir_extd(1:grbflen) == ' ' ) THEN
dir_extd = '.'
grbflen=1
END IF
IF( dir_extd(grbflen:grbflen) /= '/'.AND.grbflen < len1 ) THEN
grbflen=grbflen+1
dir_extd(grbflen:grbflen)='/'
END IF
lenrun = LEN( extdname )
CALL strlnth( extdname, lenrun )
gribfile = dir_extd(1:grbflen)//extdname(1:lenrun) &
//'.'//gribtime(1:13)
grbflen = grbflen + lenrun + 14
!
!-----------------------------------------------------------------------
!
! RDNMCGRB reads NMC GRIB data
!
!-----------------------------------------------------------------------
!
govrd = g/rd
gridtyp = eta212grid
mproj_grb = eta212proj
n2dvars = eta212nvs2d
n2dlvtps = eta212nlvt2d
DO k=1,n2dlvtps
DO n=1,n2dvars
var_id2d(n,k) = eta212var_id2d(n,k)
END DO
levtyp2d(k) = eta212levs2d(k)
END DO
n3dvars = eta212nvs3d
n3dlvtps = eta212nlvt3d
DO m=1,n3dlvtps
DO n=1,n3dvars
var_id3d(n,m) = eta212var_id3d(n,m)
END DO
levtyp3d(m) = eta212levs3d(m)
END DO
CALL rdnmcgrb(nx_ext,ny_ext,nz_ext,gribfile,grbflen, gribtime, &
gridesc, iproj_grb, gthin, &
ni_grb,nj_grb,np_grb, nk_grb,zk_grb, npeq,nit, &
pi_grb,pj_grb,ipole, di_grb,dj_grb, &
latsw,lonsw, latne,lonne, &
latrot,lonrot,angrot, &
latstr,lonstr,facstr, &
lattru1,lattru2,lontrue, &
scanmode, iscan,jscan,kscan, &
ires,iearth,icomp, &
jpenta,kpenta,mpenta,ispect,icoeff, &
xp_grb,yp_grb, xo_grb,yo_grb,zo_grb, &
rcdata,var_grb2d,var_grb3d,var_lev3d,iret)
max_nr2d = 0
DO n=1,n2dvars
DO m=1,n2dlvtps
max_nr2d = MAX( max_nr2d, var_nr2d(n,m) )
END DO
END DO
max_nr3d = 0
DO n=1,n3dvars
max_nr3d = MAX( max_nr3d, var_nr3d(n,1) )
END DO
IF ( max_nr3d == 0 ) THEN
WRITE (6,'(a)') &
'No 3-D variable was found in the GRIB file', &
'Program stopped in GETNMCETA212.'
STOP
END IF
IF ( max_nr2d == 0 ) THEN
WRITE (6,'(a)') &
'No 2-D variables was found in the GRIB file'
END IF
! write (6,'(/a7,2x,6(i7))')
! : 'Lev\\VID',(var_id3d(n,1),n=1,n3dvars)
! DO 60 k=1,max_nr3d
! write (6,'(/i5,4x,6(i7))')
! : k,(var_lev3d(k,n,1),n=1,n3dvars)
60 CONTINUE
DO k=1,max_nr3d
DO n=2,n3dvars
IF ( var_lev3d(k,1,1) /= var_lev3d(k,n,1) ) THEN
WRITE (6,'(a)') &
'Variables were not at the same level.', &
'Program stopped in GETNMCETA212.'
STOP
END IF
END DO
END DO
IF ( iproj_grb == 5 .AND. ipole == 0 ) THEN ! Center in NP
iproj_ext = 1
ELSE IF ( iproj_grb == 5 .AND. ipole == 1 ) THEN ! Center in SP
iproj_ext = -1
ELSE IF ( iproj_grb == 3 .AND. ipole == 0 ) THEN ! Center in NP
iproj_ext = 2
ELSE IF ( iproj_grb == 3 .AND. ipole == 1 ) THEN ! Center in SP
iproj_ext = -2
ELSE IF ( iproj_grb == 1 ) THEN
iproj_ext = 3
ELSE IF ( iproj_grb == 0 ) THEN
iproj_ext = 4
ELSE
WRITE (6,'(a)') &
'Unknown map projection. Set to non-projection.'
iproj_ext = 0
END IF
scale_ext = 1.0
latnot_ext(1) = lattru1
latnot_ext(2) = lattru2
trlon_ext = lontrue
dx_ext = di_grb
dy_ext = dj_grb
CALL getmapr(iproj,scale,latnot,trlon,x0,y0)
CALL setmapr(iproj_ext,scale_ext,latnot_ext,trlon_ext)
CALL lltoxy(1,1,latsw,lonsw,x0_ext,y0_ext)
DO i=1,nx_ext
x_ext(i)=x0_ext+(i-1)*dx_ext
END DO
DO j=1,ny_ext
y_ext(j)=y0_ext+(j-1)*dy_ext
END DO
CALL xytoll(nx_ext,ny_ext,x_ext,y_ext,lat_ext,lon_ext)
!---------------------------------------------------------------------
! Define soil depths
!---------------------------------------------------------------------
DO k=1,nzsoil_ext
DO j=1,ny_ext
DO i=1,nx_ext
zpsoil_ext(i,j,k) = soildepth_ext(k)
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Retrieve 2-D variables
!
!-----------------------------------------------------------------------
!
DO j=1,ny_ext
DO i=1,nx_ext
IF ( var_nr2d(1,1) == 0 ) THEN
psfc_ext (i,j) = -999.0
ELSE
psfc_ext (i,j) = var_grb2d(i,j,1,1) * 100.0
END IF
IF ( var_nr2d(2,1) == 0 ) THEN
trn_ext (i,j) = -999.0
ELSE
trn_ext (i,j) = var_grb2d(i,j,2,1)
END IF
IF ( var_nr3d(1,2) == 0 ) THEN
DO k=1,nzsoil_ext
tsoil_ext (i,j,k) = -999.0
qsoil_ext (i,j,k) = -999.0
END DO
ELSE
IF (soilscheme == 0) THEN ! Old ARPS Force-Restore Soil Model
tsoil_ext(i,j,1) = var_grb2d(i,j,3,1)
IF ( nint(var_grb2d(i,j,5,1)) == 1 ) THEN ! soil temp over land
tsoil_ext(i,j,2) = 0.1 * var_grb3d(i,j,1,1,2) & ! 0-10cm
+ 0.3 * var_grb3d(i,j,2,1,2) & ! 10-40cm
+ 0.6 * var_grb3d(i,j,3,1,2) ! 40-100cm
soiltyp_ext(i,j) = 0
ELSE ! sfc temp over sea
tsoil_ext(i,j,2) = var_grb2d(i,j,3,1)
soiltyp_ext(i,j) = 13 ! Set soil type to water
END IF
qsoil_ext(i,j,1) = var_grb3d(i,j,1,2,2)
qsoil_ext(i,j,2) = 0.1 * var_grb3d(i,j,1,2,2) & ! 0-10cm
+ 0.3 * var_grb3d(i,j,2,2,2) & ! 10-40cm
+ 0.6 * var_grb3d(i,j,3,2,2) ! 40-100cm
ELSE ! OU Soil model
!EMK 15 June 2002...Reorganized, added another soil level, and added
!water temperature.
IF ( nint(var_grb2d(i,j,5,1)) == 1 ) THEN ! Land
tsoil_ext (i,j,1) = var_grb2d(i,j,3,1) ! Ground temperature
qsoil_ext (i,j,1) = var_grb3d(i,j,1,2,2) ! Assumed to be same
! as first below ground
! level.
DO k=2,nzsoil_ext
! "TSOIL" in GRIB is below ground, treated as separate
! variable from ground temperature.
tsoil_ext (i,j,k) = var_grb3d(i,j,k-1,1,2) ! Not a bug;
qsoil_ext (i,j,k) = var_grb3d(i,j,k-1,2,2)
END DO
ELSE ! Water
DO k=1,nzsoil_ext
tsoil_ext (i,j,k) = var_grb2d(i,j,7,1) ! Water temperature
qsoil_ext (i,j,k) = 1. ! 100% water
END DO
END IF ! Land or water?
! qsoil_ext (i,j,1) = var_grb3d(i,j,1,2,2) ! 0 cm
! DO k=2,nzsoil_ext
! qsoil_ext (i,j,k) = var_grb3d(i,j,k-1,2,2)
! END DO
!EMK END 15 June 2002
END IF ! soilscheme
END IF
IF ( var_nr2d(4,1) == 0 ) THEN
wetcanp_ext(i,j) = -999.0
ELSE
wetcanp_ext(i,j) = var_grb2d(i,j,4,1)*1.e-3 ! in meter
END IF
IF ( var_nr2d(6,1) == 0 ) THEN
snowdpth_ext(i,j) = -999
ELSE
! Convert water equiv. of accum. snow depth (kg/m**2) to meters
! (where 1 meter liquid water is set equivqlent to 10 meters snow).
! 0.01 = 10. (m snow/m liquid) / (1000 kg/m**3)
snowdpth_ext(i,j) = 0.01 * var_grb2d(i,j,6,1) ! in meters
END IF
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Retrieve 3-D variables
!
!-----------------------------------------------------------------------
!
nz1 = MIN(var_nr3d(1,1),nz_ext)
IF ( var_lev3d(1,1,1) > var_lev3d(nz1,1,1) ) THEN ! 1st level at sfc
chklev = 1
lvscan = 0
ELSE
chklev = -1
lvscan = nz1+1
END IF
DO k=1,nz1
kk = chklev * k + lvscan
DO j=1,ny_ext
DO i=1,nx_ext
! p_ext (i,j,kk) = psfc_ext(i,j)
! : * float(var_lev3d(k,1))/10000.
! Pressure derived from sigma ccordinates
p_ext (i,j,kk) = 100.0 * FLOAT(var_lev3d(k,1,1)) ! Pressure
t_ext (i,j,kk) = var_grb3d(i,j,k,1,1) ! Temperature (K)
qv_ext (i,j,kk) = var_grb3d(i,j,k,2,1) ! Spec. humidity
! (kg/kg)
u_ext (i,j,kk) = var_grb3d(i,j,k,3,1) ! u wind (m/s)
v_ext (i,j,kk) = var_grb3d(i,j,k,4,1) ! v wind (m/s)
hgt_ext(i,j,kk) = var_grb3d(i,j,k,5,1) ! height (m)
qc_ext (i,j,kk) = -999.
qr_ext (i,j,kk) = -999.
qi_ext (i,j,kk) = -999.
qs_ext (i,j,kk) = -999.
qh_ext (i,j,kk) = -999.
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Calculate the height by using hytrostatical equation
!
!-----------------------------------------------------------------------
!
! DO 220 j=1,ny_ext
! DO 220 i=1,nx_ext
! tema = 0.5 * ( tsfc_ext (i,j) + t_ext (i,j,1) )
! temb = 0.5 * ( qvsfc_ext(i,j) + qv_ext(i,j,1) )
! temc = tema * (1.0 + 0.608*temb) ! Virtual temperature
! tema = 0.5 * ( psfc_ext (i,j) + p_ext (i,j,1) )
! temb = temc / tema
!
! hgt_ext(i,j,1) = trn_ext(i,j) + govrd * temb ! Height (m)
! : * ( psfc_ext(i,j) - p_ext(i,j,1) )
!220 CONTINUE
!
! DO 230 k=2,nz1
! DO 230 j=1,ny_ext
! DO 230 i=1,nx_ext
! tema = 0.5 * ( t_ext (i,j,k-1) + t_ext (i,j,k) )
! temb = 0.5 * ( qv_ext(i,j,k-1) + qv_ext(i,j,k) )
! temc = tema * (1.0 + 0.608*temb) ! Virtual temperature
! tema = 0.5 * ( p_ext (i,j,k-1) + p_ext (i,j,k) )
! temb = temc / tema
!
! hgt_ext(i,j,k) = hgt_ext(i,j,k-1) + govrd * temb ! Height (m)
! : * ( p_ext (i,j,k-1) - p_ext (i,j,k) )
!230 CONTINUE
!
!
!-----------------------------------------------------------------------
!
! Rotate winds to be relative to true north.
! The RUC data are sent as grid-relative.
!
!-----------------------------------------------------------------------
!
DO k=1,nz1
CALL uvmptoe(nx_ext,ny_ext,u_ext(1,1,k),v_ext(1,1,k), &
lon_ext,utmp,vtmp)
u_ext(:,:,k) = utmp(:,:)
v_ext(:,:,k) = vtmp(:,:)
END DO
!
!-----------------------------------------------------------------------
!
! Reset map projection to previous values
!
!-----------------------------------------------------------------------
!
CALL setmapr(iproj,scale,latnot,trlon)
CALL setorig(1,x0,y0)
istatus = 1
DEALLOCATE(var_grb2d,var_grb3d,rcdata,var_lev3d)
DEALLOCATE(utmp)
DEALLOCATE(vtmp)
RETURN
END SUBROUTINE getnmceta212
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE SWAPJ ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE swapj(nx,ny,nz,varval)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! This routine swaps values in the y direction so that the
! reanalysis grid goes S to N in the +j direction
!
!-----------------------------------------------------------------------
!
! AUTHOR: Dennis Moon, SSESCO
! 06/19/1998
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT:
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER*4 nx,ny,nz,i,j,k
REAL*4 varval(nx,ny,nz), tmpval
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
DO k=1,nz
DO i=1,nx
DO j=1,ny/2
tmpval = varval(i,ny-j+1,k)
varval(i,ny-j+1,k) = varval(i,j,k)
varval(i,j,k) = tmpval
END DO
END DO
END DO
RETURN
END SUBROUTINE swapj
!
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE GETNMCRUCN236 ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE getnmcrucn236(nx_ext,ny_ext,nz_ext, & 1,11
dir_extd,extdname,extdopt,extdfmt, &
extdinit,extdfcst,julfname, &
iproj_ext,scale_ext, &
trlon_ext,latnot_ext,x0_ext,y0_ext, &
lat_ext,lon_ext, &
p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext, &
qc_ext,qr_ext,qi_ext,qs_ext,qh_ext, &
tsfc_ext,tdeep_ext,wetsfc_ext,wetdp_ext,wetcanp_ext, &
trn_ext,psfc_ext,snowdpth_ext, &
istatus)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Reads in a GRIB file containing native coordinate RUC2 data
! (Grid 236) and extracts/converts selected variables for use by
! EXT2ARPS.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Eric Kemp
! 09/17/1999
! Based on subroutine GETNMCRUC87.
!
! MODIFICATION HISTORY:
!
! 11/01/1999 Eric Kemp
! Corrected several bugs in variable retrievals, and added
! retrieval of RUC2 hydrometeor fields.
!
! 1999/11/30 Gene Bassett
! Corrected RUC2 soil moisture and changed deep soil moisture &
! temperature to be an average from 0 to 100cm.
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! dir_extd Directory name for external file
! extdname Prefix string of external file name
! extdopt Option of external data sources
! extdfmt Option of external data format
! extdinit Initialized time in mm-dd-yyyy:hh:mm:ss format
! extdfcst Forecast hour in HHH:MM:SS format
! julfname File name in yyjjjhhmm format
!
! OUTPUT:
!
! iproj_ext Map projection number of external data
! scale_ext Scale factor of external data
! trlon_ext True longitude of external data (degrees E)
! latnot_ext(2) True latitude(s) of external data (degrees N)
! x0_ext x coordinate of origin of external data
! y0_ext y coordinate of origin of external data
! lat_ext latitude of external data points (degrees N)
! lon_ext longitude of external data points (degrees E)
! p_ext pressure (Pascal)
! hgt_ext height (m)
! t_ext temperature (K)
! qv_ext specific humidity (kg/kg)
! u_ext u wind component (m/s)
! v_ext v wind component (m/s)
! qc_ext Cloud water mixing ratio (kg/kg)
! qr_ext Rain water mixing ratio (kg/kg)
! qi_ext Ice mixing ratio (kg/kg)
! qs_ext Snow mixing ratio (kg/kg)
! qh_ext Hail mixing ratio (kg/kg)
!
! tsfc_ext Surface temperature
! tdeep_ext Soil temperature
! wetsfc_ext Top layer soil moisture (fraction)
! wetdp_ext Deep soil moisture (fraction)
! wetcanp_ext Water content on canopy
!
! trn_ext External terrain (m)
! psfc_ext Surface pressure (Pa)
! snowdpth_ext Snow depth (m)
!
! istatus status indicator
!
! WORK ARRAYS:
!
! var_grb3d Arrays to store the GRIB 3-D variables:
! var_grb3d(nxgrb,nygrb,nzgrb,1,1) - pressure (Pa)
! var_grb3d(nxgrb,nygrb,nzgrb,2,1) - height (m)
! var_grb3d(nxgrb,nygrb,nzgrb,3,1) - Virtual Potential
! temperature (K)
! var_grb3d(nxgrb,nygrb,nzgrb,4,1) - Water vapor
! mixing ratio
! (kg/kg)
! var_grb3d(nxgrb,nygrb,nzgrb,5,1) - u wind (m/s)
! var_grb3d(nxgrb,nygrb,nzgrb,6,1) - v wind (m/s)
! var_grb3d(nxgrb,nygrb,nzgrb,7,1)
! - cloud water mixing ratio (kg/kg)
! var_grb3d(nxgrb,nygrb,nzgrb,8,1)
! - rain water mixing ratio (kg/kg)
! var_grb3d(nxgrb,nygrb,nzgrb,9,1)
! - ice mixing ratio (kg/kg)
! var_grb3d(nxgrb,nygrb,nzgrb,10,1)
! - snow mixing ratio (kg/kg)
! var_grb3d(nxgrb,nygrb,nzgrb,11,1)
! - graupel mixing ratio (kg/kg)
! var_grb3d(nxgrb,nygrb,nzgrb,1,2) - soil temp. (K)
! var_grb3d(nxgrb,nygrb,nzgrb,2,2) - soil moist.
! (fraction)
!
! var_grb2d Arrays to store the GRIB 2-D variables:
! var_grb2d(nxgrb,nygrb,1) - Canopy Water
! (kg/m**2)
! var_grb2d(nxgrb,nygrb,2) - Snow depth (m)
! var_grb2d(nxgrb,nygrb,3) - Surface soil temperature (K)
! var_grb2d(nxgrb,nygrb,4) - Surface soil moisture
! (fraction)
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INCLUDE 'gribcst.inc'
CHARACTER (LEN=*) :: dir_extd
CHARACTER (LEN=*) :: extdname
INTEGER :: extdopt
INTEGER :: extdfmt
CHARACTER (LEN=19) :: extdinit
CHARACTER (LEN=9) :: extdfcst
CHARACTER (LEN=9) :: julfname
!
!
!-----------------------------------------------------------------------
!
! External grid variables
!
!-----------------------------------------------------------------------
!
INTEGER :: iproj_ext
REAL :: scale_ext,trlon_ext
REAL :: latnot_ext(2)
REAL :: x0_ext,y0_ext
REAL :: dx_ext,dy_ext
!
!-----------------------------------------------------------------------
!
! Output external variable arrays
!
!-----------------------------------------------------------------------
!
INTEGER :: nx_ext, ny_ext, nz_ext
REAL :: lat_ext(nx_ext,ny_ext)
REAL :: lon_ext(nx_ext,ny_ext)
REAL :: p_ext (nx_ext,ny_ext,nz_ext) ! Pressure (Pascals)
REAL :: hgt_ext(nx_ext,ny_ext,nz_ext) ! Height (m)
REAL :: t_ext (nx_ext,ny_ext,nz_ext) ! Temperature (K)
REAL :: qv_ext (nx_ext,ny_ext,nz_ext) ! Specific humidity (kg/kg)
REAL :: u_ext (nx_ext,ny_ext,nz_ext) ! Eastward wind component
REAL :: v_ext (nx_ext,ny_ext,nz_ext) ! Northward wind component
REAL :: qc_ext (nx_ext,ny_ext,nz_ext) ! Cloud H2O mixing ratio (kg/kg)
REAL :: qr_ext (nx_ext,ny_ext,nz_ext) ! Rain H2O mixing ratio (kg/kg)
REAL :: qi_ext (nx_ext,ny_ext,nz_ext) ! Ice mixing ratio (kg/kg)
REAL :: qs_ext (nx_ext,ny_ext,nz_ext) ! Snow mixing ratio (kg/kg)
REAL :: qh_ext (nx_ext,ny_ext,nz_ext) ! Hail mixing ratio (kg/kg)
REAL :: tsfc_ext (nx_ext,ny_ext) ! Temperature at surface (K)
REAL :: tdeep_ext (nx_ext,ny_ext) ! Deep soil temperature (K)
REAL :: wetsfc_ext (nx_ext,ny_ext) ! Surface soil moisture (fraction)
REAL :: wetdp_ext (nx_ext,ny_ext) ! Deep soil moisture (fraction)
REAL :: wetcanp_ext(nx_ext,ny_ext) ! Canopy water amount
REAL :: trn_ext (nx_ext,ny_ext) ! Geometrical heights
REAL :: psfc_ext (nx_ext,ny_ext) ! Surface pressure (Pa)
REAL :: snowdpth_ext (nx_ext,ny_ext) ! Snow depth (m)
!
!-----------------------------------------------------------------------
!
! Other external variable arrays
!
!-----------------------------------------------------------------------
!
REAL :: x_ext(nx_ext)
REAL :: y_ext(ny_ext)
INTEGER :: istatus
!
!-----------------------------------------------------------------------
!
! Work arrays for storing grib data
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: var_grb2d(:,:,:,:) ! GRIB variables
REAL, ALLOCATABLE :: var_grb3d(:,:,:,:,:) ! GRIB 3-D variables
INTEGER, ALLOCATABLE :: var_lev3d(:,:,:) ! Levels (hybrid) for
! each 3-D variable
REAL, ALLOCATABLE :: rcdata(:) ! temporary data array
!
!-----------------------------------------------------------------------
!
! Original grid variables
!
!-----------------------------------------------------------------------
!
INTEGER :: iproj
REAL :: scale,trlon,x0,y0
REAL :: latnot(2)
!
!-----------------------------------------------------------------------
!
! Misc internal variables
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=80) :: gribfile
CHARACTER (LEN=13) :: gribtime
INTEGER :: i,j,k
INTEGER :: iyr,imo,iday,myr, jldy
INTEGER :: ihr,imin,isec
INTEGER :: ifhr,ifmin,ifsec
INTEGER :: grbflen, len1, lenrun
INTEGER :: m,n,nz1,max_nr2d,max_nr3d
REAL :: pilev
REAL :: tv_ext, tvc_ext
REAL :: rovcp_p, cpd_p, g0_p, rd_p
INTEGER :: chklev, lvscan, kk, jj
REAL :: tema, temb
REAL :: a,b
INTEGER :: iret ! Return flag
REAL, ALLOCATABLE :: utmp(:,:), vtmp(:,:)
!
!-----------------------------------------------------------------------
!
! Include files
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
! GRIB grid information
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=42) :: gridesc ! Grid description
INTEGER :: iproj_grb ! Map projection indicator
INTEGER :: gthin ! Indicator of whether the grid is "thinned"
INTEGER :: ni_grb ! Number of points along x-axis
INTEGER :: nj_grb ! Number of points along y-axis
INTEGER :: np_grb ! Total number of horizontal grid points
INTEGER :: nk_grb ! Number of vertical parameters
REAL :: zk_grb(nz_ext) ! Vertical coordinate parameters
INTEGER :: npeq ! Number of lat circles from pole to equator
INTEGER :: nit(nz_ext) ! Number of x-points for thinned grid
REAL :: pi_grb ! x-coordinate of pole point
REAL :: pj_grb ! y-coordinate of pole point
INTEGER :: ipole ! Projection center flag
REAL :: di_grb ! x-direction increment or grid length
REAL :: dj_grb ! y-direction increment or grid length
REAL :: latsw ! Latitude of South West corner point
REAL :: lonsw ! Longitude of South West corner point
REAL :: latne ! Latitude of North East corner point
REAL :: lonne ! Longitude of North East corner point
REAL :: lattru1 ! Latitude (1st) at which projection is true
REAL :: lattru2 ! Latitude (2nd) at which projection is true
REAL :: lontrue ! Longitude at which projection is true
REAL :: latrot ! Latitude of southern pole of rotation
REAL :: lonrot ! Longitude of southern pole of rotation
REAL :: angrot ! Angle of rotation
REAL :: latstr ! Latitude of the pole of stretching
REAL :: lonstr ! Longitude of the pole of stretching
REAL :: facstr ! Stretching factor
INTEGER :: scanmode ! Scanning indicator
INTEGER :: iscan ! x-direction scanning indicator
INTEGER :: jscan ! y-direction scanning indicator
INTEGER :: kscan ! FORTRAN index scanning indicator
INTEGER :: ires ! Resolution direction increments indicator
INTEGER :: iearth ! Earth shape indicator: spherical or oblate?
INTEGER :: icomp ! (u,v) components decomposition indicator
INTEGER :: jpenta ! J-Pentagonal resolution parameter
INTEGER :: kpenta ! K-Pentagonal resolution parameter
INTEGER :: mpenta ! M-Pentagonal resolution parameter
INTEGER :: ispect ! Spectral representation type
INTEGER :: icoeff ! Spectral coefficient storage mode
REAL :: xp_grb ! X coordinate of sub-satellite point
REAL :: yp_grb ! Y coordinate of sub-satellite point
REAL :: xo_grb ! X coordinate of image sector origin
REAL :: yo_grb ! Y coordinate of image sector origin
REAL :: zo_grb ! Camera altitude from center of Earth
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
ALLOCATE(var_grb2d(nx_ext,ny_ext,n2dvs,n2dlvt))
ALLOCATE(var_grb3d(nx_ext,ny_ext,nz_ext,n3dvs,n3dlvt))
ALLOCATE(rcdata(nx_ext*ny_ext))
ALLOCATE(var_lev3d(nz_ext,n3dvs,n3dlvt))
ALLOCATE(utmp(nx_ext,ny_ext))
ALLOCATE(vtmp(nx_ext,ny_ext))
IF ( extdfcst == ' ') extdfcst='000:00:00'
READ (extdinit,'(i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,i2)') &
iyr,imo,iday,ihr,imin,isec
CALL julday(iyr,imo,iday,jldy)
myr=MOD(iyr,100)
ifhr=0
ifmin=0
ifsec=0
READ (extdfcst,'(i3)',ERR=4,END=4) ifhr
4 CONTINUE
WRITE (gribtime,'(i4.4,i2.2,i2.2,i2.2,a1,i2.2)') &
iyr,imo,iday,ihr,'f',ifhr
len1=LEN(dir_extd)
grbflen=len1
CALL strlnth( dir_extd, grbflen )
IF( grbflen == 0 .OR. dir_extd(1:grbflen) == ' ' ) THEN
dir_extd = '.'
grbflen=1
END IF
IF( dir_extd(grbflen:grbflen) /= '/'.AND.grbflen < len1 ) THEN
grbflen=grbflen+1
dir_extd(grbflen:grbflen)='/'
END IF
lenrun = LEN( extdname )
CALL strlnth( extdname, lenrun )
gribfile = dir_extd(1:grbflen)//extdname(1:lenrun) &
//'.'//gribtime(1:13)
grbflen = grbflen + lenrun + 14
!
!-----------------------------------------------------------------------
!
! RDNMCGRB reads NMC GRIB data
!
!-----------------------------------------------------------------------
!
gridtyp = rucn236grid
mproj_grb = rucn236proj
n2dvars = rucn236nvs2d
n2dlvtps = rucn236nlvt2d
DO m=1,n2dlvtps
DO n=1,n2dvars
var_id2d(n,m) = rucn236var_id2d(n,m)
END DO
levtyp2d(m) = rucn236levs2d(m)
END DO
n3dvars = rucn236nvs3d
n3dlvtps = rucn236nlvt3d
DO m=1,n3dlvtps
DO n=1,n3dvars
var_id3d(n,m) = rucn236var_id3d(n,m)
END DO
levtyp3d(m) = rucn236levs3d(m)
END DO
CALL rdnmcgrb(nx_ext,ny_ext,nz_ext,gribfile,grbflen, gribtime, &
gridesc, iproj_grb, gthin, &
ni_grb,nj_grb,np_grb, nk_grb,zk_grb, npeq,nit, &
pi_grb,pj_grb,ipole, di_grb,dj_grb, &
latsw,lonsw, latne,lonne, &
latrot,lonrot,angrot, &
latstr,lonstr,facstr, &
lattru1,lattru2,lontrue, &
scanmode, iscan,jscan,kscan, &
ires,iearth,icomp, &
jpenta,kpenta,mpenta,ispect,icoeff, &
xp_grb,yp_grb, xo_grb,yo_grb,zo_grb, &
rcdata,var_grb2d,var_grb3d,var_lev3d,iret)
max_nr2d = 0
DO n=1,n2dvars
DO m=1,n2dlvtps
max_nr2d = MAX( max_nr2d, var_nr2d(n,m) )
END DO
END DO
max_nr3d = 0
DO n=1,n3dvars
max_nr3d = MAX( max_nr3d, var_nr3d(n,1) )
END DO
IF ( max_nr3d == 0 ) THEN
WRITE (6,'(a)') &
'No 3-D variables was found in the GRIB file', &
'Program stopped in GETNMCRUC.'
STOP
END IF
IF ( max_nr2d == 0 ) THEN
WRITE (6,'(a)') &
'No 2-D variables was found in the GRIB file'
END IF
! write (6,'(/a7,2x,6(i7))')
! : 'Lev\\VID',(var_id3d(n,1),n=1,n3dvars)
! DO 60 k=1,max_nr3d
! write (6,'(i5,4x,6(i7))')
! : k,(var_lev3d(k,n,1),n=1,n3dvars)
60 CONTINUE
DO k=1,max_nr3d
DO n=2,n3dvars
IF ( var_lev3d(k,1,1) /= var_lev3d(k,n,1) ) THEN
WRITE (6,'(a)') &
'Variables were not at the same level.', &
'Program stopped in GETNMCRUC.'
WRITE(6,*) &
'var_lev3d(k,1,1) = ',var_lev3d(k,1,1), &
'var_lev3d(k,n,1) = ',var_lev3d(k,n,1)
STOP
END IF
END DO
END DO
IF ( iproj_grb == 5 .AND. ipole == 0 ) THEN ! Center in NP
iproj_ext = 1
ELSE IF ( iproj_grb == 5 .AND. ipole == 1 ) THEN ! Center in SP
iproj_ext = -1
ELSE IF ( iproj_grb == 3 .AND. ipole == 0 ) THEN ! Center in NP
iproj_ext = 2
ELSE IF ( iproj_grb == 3 .AND. ipole == 1 ) THEN ! Center in SP
iproj_ext = -2
ELSE IF ( iproj_grb == 1 ) THEN
iproj_ext = 3
ELSE IF ( iproj_grb == 0 ) THEN
iproj_ext = 4
ELSE
WRITE (6,'(a)') &
'Unknown map projection. Set to non-projection.'
iproj_ext = 0
END IF
scale_ext = 1.0
latnot_ext(1) = lattru1
latnot_ext(2) = lattru2
trlon_ext = lontrue
dx_ext = di_grb
dy_ext = dj_grb
CALL getmapr(iproj,scale,latnot,trlon,x0,y0)
CALL setmapr(iproj_ext,scale_ext,latnot_ext,trlon_ext)
CALL lltoxy(1,1,latsw,lonsw,x0_ext,y0_ext)
DO i=1,nx_ext
x_ext(i)=x0_ext+(i-1)*dx_ext
END DO
DO j=1,ny_ext
y_ext(j)=y0_ext+(j-1)*dy_ext
END DO
CALL xytoll(nx_ext,ny_ext,x_ext,y_ext,lat_ext,lon_ext)
!
!-----------------------------------------------------------------------
!
! Retrieve 2-D variables
!
!-----------------------------------------------------------------------
!
DO j=1,ny_ext
DO i=1,nx_ext
tsfc_ext(i,j) = -999.0
trn_ext(i,j) = -999.0
IF ( var_nr3d(1,2) == 0 ) THEN
tsfc_ext (i,j) = -999.0
tdeep_ext (i,j) = -999.0
wetsfc_ext (i,j) = -999.0
wetdp_ext (i,j) = -999.0
wetcanp_ext(i,j) = -999.0
ELSE
! tsfc_ext (i,j) = var_grb3d(i,j,1,1,2)
tsfc_ext (i,j) = var_grb2d(i,j,3,1)
tdeep_ext (i,j) = 0.1 * var_grb3d(i,j,1,1,2) & ! 5cm
+ 0.2 * var_grb3d(i,j,2,1,2) & ! 20cm
+ 0.4 * var_grb3d(i,j,3,1,2) & ! 40cm
+ 0.3 * var_grb3d(i,j,4,1,2) ! 160cm
! wetsfc_ext (i,j) = var_grb3d(i,j,1,2,2)
wetsfc_ext (i,j) = var_grb2d(i,j,4,1)
wetdp_ext (i,j) = 0.1 * var_grb3d(i,j,1,2,2) & ! 5cm
+ 0.2 * var_grb3d(i,j,2,2,2) & ! 20cm
+ 0.4 * var_grb3d(i,j,3,2,2) & ! 40cm
+ 0.3 * var_grb3d(i,j,4,2,2) ! 160cm
wetcanp_ext(i,j) = var_grb2d(i,j,1,1)*1.e-3 ! in meters
END IF
psfc_ext (i,j) = -999.0
IF ( var_nr2d(2,1) == 0 ) THEN
snowdpth_ext(i,j) = -999
ELSE
snowdpth_ext(i,j) = var_grb2d(i,j,2,1) ! in meters
END IF
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Retrieve 3-D variables
!
!-----------------------------------------------------------------------
!
cpd_p = 1004.686 ! cp in RUC
rovcp_p = 0.285714 ! rd/cp used in RUC
g0_p = 9.80665 ! gravity in RUC
nz1 = MIN(var_nr3d(1,1),nz_ext)
IF ( var_lev3d(1,1,1) < var_lev3d(nz1,1,1) ) THEN ! 1st level at sfc
chklev = 1
lvscan = 0
ELSE
chklev = -1
lvscan = nz1+1
END IF
DO k=1,nz1
kk = chklev * k + lvscan
DO j=1,ny_ext
DO i=1,nx_ext
p_ext(i,j,kk) = var_grb3d(i,j,k,1,1) ! Pressure (Pa)
hgt_ext(i,j,kk) = var_grb3d(i,j,k,2,1) ! Height (m)
u_ext(i,j,kk) = var_grb3d(i,j,k,5,1) ! u wind (m/s)
v_ext(i,j,kk) = var_grb3d(i,j,k,6,1) ! v wind (m/s)
a = REAL(100000)/var_grb3d(i,j,k,1,1)
a = a**rovcp_p
tvc_ext = var_grb3d(i,j,k,3,1)/a ! Virtual Temperature
b = 0.61*var_grb3d(i,j,k,4,1)
b = REAL(1) + b
t_ext(i,j,kk) = tvc_ext/b ! Temperature (K)
a = var_grb3d(i,j,k,4,1)*var_grb3d(i,j,k,1,1)
a = a/(0.622 - var_grb3d(i,j,k,4,1))
qv_ext(i,j,kk) = a*0.622/var_grb3d(i,j,k,1,1) ! SpecificHumidity
!
!-----------------------------------------------------------------------
!
! Retrieve hydrometeor data.
!
!-----------------------------------------------------------------------
!
IF (var_nr3d(7,1) == 0) THEN
qc_ext(i,j,kk) = -999.
ELSE
qc_ext(i,j,kk) = var_grb3d(i,j,k,7,1)
END IF
IF (var_nr3d(8,1) == 0) THEN
qr_ext(i,j,kk) = -999.
ELSE
qr_ext(i,j,kk) = var_grb3d(i,j,k,8,1)
END IF
IF (var_nr3d(9,1) == 0) THEN
qi_ext(i,j,kk) = -999.
ELSE
qi_ext(i,j,kk) = var_grb3d(i,j,k,9,1)
END IF
IF (var_nr3d(10,1) == 0) THEN
qs_ext(i,j,kk) = -999.
ELSE
qs_ext(i,j,kk) = var_grb3d(i,j,k,10,1)
END IF
IF (var_nr3d(11,1) == 0) THEN
qh_ext(i,j,kk) = -999.
ELSE
qh_ext(i,j,kk) = var_grb3d(i,j,k,11,1)
END IF
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Rotate winds to be relative to true north.
! The RUC data are sent as grid-relative.
!
!-----------------------------------------------------------------------
!
DO k=1,nz1
CALL uvmptoe(nx_ext,ny_ext,u_ext(1,1,k),v_ext(1,1,k), &
lon_ext,utmp,vtmp)
u_ext(:,:,k) = utmp(:,:)
v_ext(:,:,k) = vtmp(:,:)
END DO
!
!-----------------------------------------------------------------------
!
! Reset map projection to previous values
!
!-----------------------------------------------------------------------
!
CALL setmapr(iproj,scale,latnot,trlon)
CALL setorig(1,x0,y0)
istatus = 1
DEALLOCATE(var_grb2d,var_grb3d,rcdata,var_lev3d)
DEALLOCATE(utmp)
DEALLOCATE(vtmp)
RETURN
END SUBROUTINE getnmcrucn236
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE GETNMCRUCP236 ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE getnmcrucp236(nx_ext,ny_ext,nz_ext, & 1,12
dir_extd,extdname,extdopt,extdfmt, &
extdinit,extdfcst,julfname, &
iproj_ext,scale_ext, &
trlon_ext,latnot_ext,x0_ext,y0_ext, &
lat_ext,lon_ext, &
p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext, &
qc_ext,qr_ext,qi_ext,qs_ext,qh_ext, &
tsfc_ext,tdeep_ext,wetsfc_ext,wetdp_ext,wetcanp_ext, &
trn_ext,psfc_ext, t_2m_ext, rh_2m_ext, &
u_10m_ext, v_10m_ext, snowdpth_ext, &
istatus)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Reads in a GRIB file containing isobaric RUC2 data (Grid 236)
! and extracts/converts selected variables for use by EXT2ARPS.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Eric Kemp
! 09/17/1999
! Based on subroutine GETNMCRUC211.
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! dir_extd Directory name for external file
! extdname Prefix string of external file name
! extdopt Option of external data sources
! extdfmt Option of external data format
! extdinit Initialized time in mm-dd-yyyy:hh:mm:ss format
! extdfcst Forecast hour in HHH:MM:SS format
! julfname File name in yyjjjhhmm format
!
! OUTPUT:
!
! iproj_ext Map projection number of external data
! scale_ext Scale factor of external data
! trlon_ext True longitude of external data (degrees E)
! latnot_ext(2) True latitude(s) of external data (degrees N)
! x0_ext x coordinate of origin of external data
! y0_ext y coordinate of origin of external data
! lat_ext latitude of external data points (degrees N)
! lon_ext longitude of external data points (degrees E)
! p_ext pressure (Pascal)
! hgt_ext height (m)
! t_ext temperature (K)
! qv_ext specific humidity (kg/kg)
! u_ext u wind component (m/s)
! v_ext v wind component (m/s)
! qc_ext Cloud water mixing ratio (kg/kg)
! qr_ext Rain water mixing ratio (kg/kg)
! qi_ext Ice mixing ratio (kg/kg)
! qs_ext Snow mixing ratio (kg/kg)
! qh_ext Hail mixing ratio (kg/kg)
!
! tsfc_ext Surface temperature
! tdeep_ext Soil temperature
! wetsfc_ext Top layer soil moisture (fraction)
! wetdp_ext Deep soil moisture (fraction)
! wetcanp_ext Water content on canopy
!
! trn_ext External terrain (m)
! psfc_ext Surface pressure (Pa)
!
! T_2m_ext Temperature at 2m AGL
! rh_2m_ext Specific Humidity at 2m AGL
! U_10m_ext U at 10m AGL
! V_10m_ext V at 10m AGL
!
! snowdpth_ext Snow depth (m)
!
! istatus status indicator
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INCLUDE 'gribcst.inc'
CHARACTER (LEN=*) :: dir_extd
CHARACTER (LEN=*) :: extdname
INTEGER :: extdopt
INTEGER :: extdfmt
CHARACTER (LEN=19) :: extdinit
CHARACTER (LEN=9) :: extdfcst
CHARACTER (LEN=9) :: julfname
!
!-----------------------------------------------------------------------
!
! External grid variables
!
!-----------------------------------------------------------------------
!
INTEGER :: iproj_ext
REAL :: scale_ext,trlon_ext
REAL :: latnot_ext(2)
REAL :: x0_ext,y0_ext
REAL :: dx_ext,dy_ext
!
!-----------------------------------------------------------------------
!
! Output external variable arrays
!
!-----------------------------------------------------------------------
!
INTEGER :: nx_ext,ny_ext,nz_ext
REAL :: lat_ext(nx_ext,ny_ext)
REAL :: lon_ext(nx_ext,ny_ext)
REAL :: p_ext (nx_ext,ny_ext,nz_ext) ! Pressure (Pascals)
REAL :: hgt_ext(nx_ext,ny_ext,nz_ext) ! Height (m)
REAL :: t_ext (nx_ext,ny_ext,nz_ext) ! Temperature (K)
REAL :: qv_ext (nx_ext,ny_ext,nz_ext) ! Specific humidity (kg/kg)
REAL :: u_ext (nx_ext,ny_ext,nz_ext) ! Eastward wind component
REAL :: v_ext (nx_ext,ny_ext,nz_ext) ! Northward wind component
REAL :: qc_ext (nx_ext,ny_ext,nz_ext) ! Cloud H2O mixing ratio
! (kg/kg)
REAL :: qr_ext (nx_ext,ny_ext,nz_ext) ! Rain H2O mixing ratio
! (kg/kg)
REAL :: qi_ext (nx_ext,ny_ext,nz_ext) ! Ice mixing ratio (kg/kg)
REAL :: qs_ext (nx_ext,ny_ext,nz_ext) ! Snow mixing ratio (kg/kg)
REAL :: qh_ext (nx_ext,ny_ext,nz_ext) ! Hail mixing ratio (kg/kg)
REAL :: tsfc_ext (nx_ext,ny_ext) ! Temperature at surface (K)
REAL :: tdeep_ext (nx_ext,ny_ext) ! Deep soil temperature (K)
REAL :: wetsfc_ext (nx_ext,ny_ext) ! Surface soil moisture (fraction)
REAL :: wetdp_ext (nx_ext,ny_ext) ! Deep soil moisture (fraction)
REAL :: wetcanp_ext(nx_ext,ny_ext) ! Canopy water amount
REAL :: trn_ext (nx_ext,ny_ext) ! Geometrical heights
REAL :: psfc_ext (nx_ext,ny_ext) ! Surface pressure (Pa)
REAL :: t_2m_ext (nx_ext,ny_ext)
REAL :: rh_2m_ext(nx_ext,ny_ext)
REAL :: u_10m_ext(nx_ext,ny_ext)
REAL :: v_10m_ext(nx_ext,ny_ext)
REAL :: snowdpth_ext(nx_ext,ny_ext) ! Snow depth (m)
!
!----------------------------------------------------------------------
!
! Other external variable arrays
!
!----------------------------------------------------------------------
!
REAL :: x_ext(nx_ext)
REAL :: y_ext(ny_ext)
REAL :: z_ext(nz_ext)
INTEGER :: istatus
!
!-----------------------------------------------------------------------
!
! Work arrays for storing grib data
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: var_grb2d(:,:,:,:) ! GRIB variables
REAL, ALLOCATABLE :: var_grb3d(:,:,:,:,:) ! GRIB 3-D variables
INTEGER, ALLOCATABLE :: var_lev3d(:,:,:) ! Levels (hybrid) for
! each 3-D variable
REAL, ALLOCATABLE :: rcdata(:) ! temporary data array
!
!-----------------------------------------------------------------------
!
! Original grid variables
!
!-----------------------------------------------------------------------
!
INTEGER :: iproj
REAL :: scale,trlon,x0,y0
REAL :: latnot(2)
!
!-----------------------------------------------------------------------
!
! Misc internal variables
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=80) :: gribfile
CHARACTER (LEN=13) :: gribtime
CHARACTER (LEN=10) :: runstr
CHARACTER (LEN=3) :: fmtn
INTEGER :: ichr,bchar,echar
INTEGER :: i,j,k,ldir,ireturn
INTEGER :: iyr,imo,iday,myr, jldy
INTEGER :: ihr,imin,isec
INTEGER :: ifhr,ifmin,ifsec
INTEGER :: grbflen, len1, lenrun
INTEGER :: m,n,nz1,max_nr2d,max_nr3d
REAL :: qvsat, pilev, qvsatice
REAL :: tema, temb
INTEGER :: chklev, lvscan, kk, jj
INTEGER :: iret ! Return flag
REAL, ALLOCATABLE :: utmp(:,:), vtmp(:,:)
!
!-----------------------------------------------------------------------
!
! Include files
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'phycst.inc'
!
!
!-----------------------------------------------------------------------
!
! GRIB grid information
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=42) :: gridesc ! Grid description
INTEGER :: iproj_grb ! Map projection indicator
INTEGER :: gthin ! Indicator of whether the grid is "thinned"
INTEGER :: ni_grb ! Number of points along x-axis
INTEGER :: nj_grb ! Number of points along y-axis
INTEGER :: np_grb ! Total number of horizontal grid points
INTEGER :: nk_grb ! Number of vertical parameters
REAL :: zk_grb(nz_ext) ! Vertical coordinate parameters
INTEGER :: npeq ! Number of lat circles from pole to equator
INTEGER :: nit(nz_ext) ! Number of x-points for thinned grid
REAL :: pi_grb ! x-coordinate of pole point
REAL :: pj_grb ! y-coordinate of pole point
INTEGER :: ipole ! Projection center flag
REAL :: di_grb ! x-direction increment or grid length
REAL :: dj_grb ! y-direction increment or grid length
REAL :: lrb ! y-direction increment or grid length
REAL :: latsw ! Latitude of South West corner point
REAL :: lonsw ! Longitude of South West corner point
REAL :: latne ! Latitude of North East corner point
REAL :: lonne ! Longitude of North East corner point
REAL :: lattru1 ! Latitude (1st) at which projection is true
REAL :: lattru2 ! Latitude (2nd) at which projection is true
REAL :: lontrue ! Longitude at which projection is true
REAL :: latrot ! Latitude of southern pole of rotation
REAL :: lonrot ! Longitude of southern pole of rotation
REAL :: angrot ! Angle of rotation
REAL :: latstr ! Latitude of the pole of stretching
REAL :: lonstr ! Longitude of the pole of stretching
REAL :: facstr ! Stretching factor
INTEGER :: scanmode ! Scanning indicator
INTEGER :: iscan ! x-direction scanning indicator
INTEGER :: jscan ! y-direction scanning indicator
INTEGER :: kscan ! FORTRAN index scanning indicator
INTEGER :: ires ! Resolution direction increments indicator
INTEGER :: iearth ! Earth shape indicator: spherical or oblate?
INTEGER :: icomp ! (u,v) components decomposition indicator
INTEGER :: jpenta ! J-Pentagonal resolution parameter
INTEGER :: kpenta ! K-Pentagonal resolution parameter
INTEGER :: mpenta ! M-Pentagonal resolution parameter
INTEGER :: ispect ! Spectral representation type
INTEGER :: icoeff ! Spectral coefficient storage mode
REAL :: xp_grb ! X coordinate of sub-satellite point
REAL :: yp_grb ! Y coordinate of sub-satellite point
REAL :: xo_grb ! X coordinate of image sector origin
REAL :: yo_grb ! Y coordinate of image sector origin
REAL :: zo_grb ! Camera altitude from center of Earth
!
!
!-----------------------------------------------------------------------
!
! Function f_qvsatl and inline directive for Cray PVP
!
!-----------------------------------------------------------------------
!
REAL :: f_qvsatl
!fpp$ expand (f_desdt)
!fpp$ expand (f_qvsatl)
!dir$ inline always f_desdt, f_qvsatl
!*$* inline routine (f_desdt, f_qvsatl)
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
ALLOCATE(var_grb2d(nx_ext,ny_ext,n2dvs,n2dlvt))
ALLOCATE(var_grb3d(nx_ext,ny_ext,nz_ext,n3dvs,n3dlvt))
ALLOCATE(rcdata(nx_ext*ny_ext))
ALLOCATE(var_lev3d(nz_ext,n3dvs,n3dlvt))
ALLOCATE(utmp(nx_ext,ny_ext))
ALLOCATE(vtmp(nx_ext,ny_ext))
IF ( extdfcst == ' ') extdfcst='000:00:00'
READ (extdinit,'(i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,i2)') &
iyr,imo,iday,ihr,imin,isec
CALL julday(iyr,imo,iday,jldy)
myr=MOD(iyr,100)
ifhr=0
ifmin=0
ifsec=0
READ (extdfcst,'(i3)',ERR=4,END=4) ifhr
4 CONTINUE
WRITE (gribtime,'(i4.4,i2.2,i2.2,i2.2,a1,i2.2)') &
iyr,imo,iday,ihr,'f',ifhr
len1=LEN(dir_extd)
grbflen=len1
CALL strlnth( dir_extd, grbflen )
IF( grbflen == 0 .OR. dir_extd(1:grbflen) == ' ' ) THEN
dir_extd = '.'
grbflen=1
END IF
IF( dir_extd(grbflen:grbflen) /= '/'.AND.grbflen < len1 ) THEN
grbflen=grbflen+1
dir_extd(grbflen:grbflen)='/'
END IF
lenrun = LEN( extdname )
CALL strlnth( extdname, lenrun )
gribfile = &
dir_extd(1:grbflen)//extdname(1:lenrun)//'.'//gribtime(1:13)
grbflen = grbflen + lenrun + 14
!
!-----------------------------------------------------------------------
!
! RDNMCGRB reads NMC GRIB data
!
!-----------------------------------------------------------------------
!
gridtyp = rucp236grid
mproj_grb = rucp236proj
n2dvars = rucp236nvs2d
n2dlvtps = rucp236nlvt2d
DO k=1,n2dlvtps
DO n=1,n2dvars
var_id2d(n,k) = rucp236var_id2d(n,k)
END DO
levtyp2d(k) = rucp236levs2d(k)
END DO
n3dvars = rucp236nvs3d
n3dlvtps = rucp236nlvt3d
DO m=1,n3dlvtps
DO n=1,n3dvars
var_id3d(n,m) = rucp236var_id3d(n,m)
END DO
levtyp3d(m) = rucp236levs3d(m)
END DO
CALL rdnmcgrb(nx_ext,ny_ext,nz_ext,gribfile,grbflen, gribtime, &
gridesc, iproj_grb, gthin, &
ni_grb,nj_grb,np_grb, nk_grb,zk_grb, npeq,nit, &
pi_grb,pj_grb,ipole, di_grb,dj_grb, &
latsw,lonsw, latne,lonne, &
latrot,lonrot,angrot, &
latstr,lonstr,facstr, &
lattru1,lattru2,lontrue, &
scanmode, iscan,jscan,kscan, &
ires,iearth,icomp, &
jpenta,kpenta,mpenta,ispect,icoeff, &
xp_grb,yp_grb, xo_grb,yo_grb,zo_grb, &
rcdata,var_grb2d,var_grb3d,var_lev3d,iret)
max_nr2d = 0
DO n=1,n2dvars
DO m=1,n2dlvtps
max_nr2d = MAX( max_nr2d, var_nr2d(n,m) )
END DO
END DO
max_nr3d = 0
DO n=1,n3dvars
max_nr3d = MAX( max_nr3d, var_nr3d(n,1))
END DO
IF ( max_nr3d == 0 ) THEN
WRITE (6,'(a)') &
'No 3-D variables was found in the GRIB file', &
'Program stopped in GETNMCRUCP236.'
STOP
END IF
IF ( max_nr2d == 0 ) THEN
WRITE (6,'(a)') &
'No 2-D variables was found in the GRIB file'
END IF
! write (6,'(/a7,2x,6(i7))')
! : 'Lev\\VID',(var_id3d(n,1),n=1,n3dvars)
! DO 60 k=1,max_nr3d
! write (6,'(i5,4x,6(i7))')
! : k,(var_lev3d(k,n,1),n=1,n3dvars)
60 CONTINUE
DO k=1,max_nr3d
DO n=2,n3dvars
IF ( var_lev3d(k,1,1) /= var_lev3d(k,n,1) ) THEN
WRITE (6,'(a)') &
'Variables were not at the same level.', &
'Program stopped in GETNMCRUCP236.'
STOP
END IF
END DO
END DO
IF ( iproj_grb == 5 .AND. ipole == 0 ) THEN ! Center in NP
iproj_ext = 1
ELSE IF ( iproj_grb == 5 .AND. ipole == 1 ) THEN ! Center in SP
iproj_ext = -1
ELSE IF ( iproj_grb == 3 .AND. ipole == 0 ) THEN ! Center in NP
iproj_ext = 2
ELSE IF ( iproj_grb == 3 .AND. ipole == 1 ) THEN ! Center in SP
iproj_ext = -2
ELSE IF ( iproj_grb == 1 ) THEN
iproj_ext = 3
ELSE IF ( iproj_grb == 0 ) THEN
iproj_ext = 4
ELSE
WRITE (6,'(a)') &
'Unknown map projection. Set to non-projection.'
iproj_ext = 0
END IF
scale_ext = 1.0
latnot_ext(1) = lattru1
latnot_ext(2) = lattru2
trlon_ext = lontrue
dx_ext = di_grb
dy_ext = dj_grb
CALL getmapr(iproj,scale,latnot,trlon,x0,y0)
CALL setmapr(iproj_ext,scale_ext,latnot_ext,trlon_ext)
CALL lltoxy(1,1,latsw,lonsw,x0_ext,y0_ext)
DO i=1,nx_ext
x_ext(i)=x0_ext+(i-1)*dx_ext
END DO
DO j=1,ny_ext
y_ext(j)=y0_ext+(j-1)*dy_ext
END DO
CALL xytoll(nx_ext,ny_ext,x_ext,y_ext,lat_ext,lon_ext)
!
!-----------------------------------------------------------------------
!
! Retrieve 2-D variables
!
!-----------------------------------------------------------------------
!
DO j=1,ny_ext
DO i=1,nx_ext
IF ( var_nr2d(1,1) == 0 ) THEN
psfc_ext (i,j) = -999.0
ELSE
psfc_ext (i,j) = var_grb2d(i,j,1,1)
END IF
IF ( var_nr2d(2,1) == 0 ) THEN
trn_ext (i,j) = -999.0
ELSE
trn_ext (i,j) = var_grb2d(i,j,2,1)
END IF
IF( var_nr2d(1,2) == 0.) THEN
t_2m_ext(i,j)= -999.0
ELSE
t_2m_ext(i,j)= var_grb2d(i,j,1,2)
END IF
! at this point we are reading in the relative humidity
! later we'll convert to specific humidity
IF( var_nr2d(2,2) == 0.) THEN
rh_2m_ext(i,j)= -999.0
ELSE
rh_2m_ext(i,j)= var_grb2d(i,j,2,2)
END IF
IF( var_nr2d(3,2) == 0.) THEN
u_10m_ext(i,j)= -999.0
ELSE
u_10m_ext(i,j)= var_grb2d(i,j,3,2)
END IF
IF( var_nr2d(4,2) == 0.) THEN
v_10m_ext(i,j)= -999.0
ELSE
v_10m_ext(i,j)= var_grb2d(i,j,4,2)
END IF
tsfc_ext (i,j) = -999.0
tdeep_ext (i,j) = -999.0
wetsfc_ext (i,j) = -999.0
wetdp_ext (i,j) = -999.0
wetcanp_ext(i,j) = -999.0
IF ( var_nr2d(3,1) == 0 ) THEN
snowdpth_ext(i,j) = -999
ELSE
snowdpth_ext(i,j) = var_grb2d(i,j,3,1) ! in meters
END IF
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Retrieve 3-D variables
!
!-----------------------------------------------------------------------
!
nz1 = MIN(var_nr3d(1,1),nz_ext)
IF ( var_lev3d(1,1,1) > var_lev3d(nz1,1,1) ) THEN ! 1st level at
!sfc
chklev = 1
lvscan = 0
ELSE
chklev = -1
lvscan = nz1+1
END IF
DO k=1,nz1
kk = chklev * k + lvscan
DO j=1,ny_ext
DO i=1,nx_ext
p_ext(i,j,kk) = 100.0 * FLOAT(var_lev3d(k,1,1)) ! Pressure
t_ext(i,j,kk) = var_grb3d(i,j,k,2,1) ! Temperature (K)
u_ext(i,j,kk) = var_grb3d(i,j,k,4,1) ! u wind (m/s)
v_ext(i,j,kk) = var_grb3d(i,j,k,5,1) ! v wind (m/s)
hgt_ext(i,j,kk) = var_grb3d(i,j,k,1,1) ! height (m)
! check for portions of constant pressure grids that are below the
! surface
!
! IF(((kk == 1) .OR. (p_ext(i,j,kk) > psfc_ext(i,j))) .AND. &
! (psfc_ext(i,j) > 0.) ) THEN
! p_ext(i,j,kk)= psfc_ext(i,j) - 1. * (kk - 1)
! t_ext(i,j,kk)= t_2m_ext(i,j)
! u_ext(i,j,kk)= u_10m_ext(i,j)
! v_ext(i,j,kk)= v_10m_ext(i,j)
! hgt_ext(i,j,kk)= trn_ext(i,j) + 0.1 * (kk - 1)
! END IF
IF( (p_ext(i,j,kk) > 0.0) .AND. (t_ext(i,j,kk) > 0.0) ) THEN
qvsat = f_qvsatl( p_ext(i,j,kk), t_ext(i,j,kk) )
qv_ext(i,j,kk)= var_grb3d(i,j,k,3,1) * qvsat * 0.01
IF((kk == 1) .OR. (p_ext(i,j,kk) > psfc_ext(i,j)))THEN
qv_ext(i,j,kk)= rh_2m_ext(i,j) * qvsat * 0.01
END IF
qc_ext(i,j,kk)= 0.0
qi_ext(i,j,kk)= 0.0
ELSE
qv_ext(i,j,kk)= 0.0
qi_ext(i,j,kk)= 0.0
qc_ext(i,j,kk)= 0.0
END IF
qr_ext (i,j,kk) = -999.
qs_ext (i,j,kk) = -999.
qh_ext (i,j,kk) = -999.
END DO
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Rotate winds to be relative to true north.
! The RUCawips data are sent as grid-relative.
!
!-----------------------------------------------------------------------
!
DO k=1,nz1
CALL uvmptoe(nx_ext,ny_ext,u_ext(1,1,k),v_ext(1,1,k), &
lon_ext,utmp,vtmp)
u_ext(:,:,k) = utmp(:,:)
v_ext(:,:,k) = vtmp(:,:)
END DO
!
!-----------------------------------------------------------------------
!
! Reset map projection to previous values
!
!-----------------------------------------------------------------------
!
CALL setmapr(iproj,scale,latnot,trlon)
CALL setorig(1,x0,y0)
istatus = 1
DEALLOCATE(var_grb2d,var_grb3d,rcdata,var_lev3d)
DEALLOCATE(utmp)
DEALLOCATE(vtmp)
RETURN
END SUBROUTINE getnmcrucp236
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE GETNCEPAVN3 ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE getncepavn3(nx_ext,ny_ext,nz_ext, & 1,5
dir_extd,extdname,extdopt,extdfmt, &
extdinit,extdfcst,julfname, &
iproj_ext,scale_ext, &
trlon_ext,latnot_ext,x0_ext,y0_ext, &
lat_ext,lon_ext, &
p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext, &
qc_ext,qr_ext,qi_ext,qs_ext,qh_ext, &
tsfc_ext,tdeep_ext,wetsfc_ext,wetdp_ext,wetcanp_ext, &
snowdpth_ext,trn_ext,psfc_ext, &
istatus)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Reads a NCEP AVN GRIB (Grid #3, 1x1 degree) data file for
! processing by ext2arps.
!
!-----------------------------------------------------------------------
!
! AUTHOR: Donghai Wang and Yuhe Liu
! 05/19/1999
!
! MODIFICATION HISTORY:
! 03/23/2000 (Donghai Wang)
! Fixed some bugs and modified the code for ARPS4.5.1 version.
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! dir_extd Directory name for external file
! extdname Prefix string of external file name
! extdopt Option of external data sources
! extdfmt Option of external data format
! extdinit Initialized time in mm-dd-yyyy:hh:mm:ss format
! extdfcst Forecast hour in HHH:MM:SS format
! julfname File name in yyjjjhhmm format
!
! OUTPUT:
!
! iproj_ext Map projection number of external data
! scale_ext Scale factor of external data
! trlon_ext True longitude of external data (degrees E)
! latnot_ext(2) True latitude(s) of external data (degrees N)
! x0_ext x coordinate of origin of external data
! y0_ext y coordinate of origin of external data
! lat_ext latitude of external data points (degrees N)
! lon_ext longitude of external data points (degrees E)
! p_ext pressure (Pascal)
! hgt_ext height (m)
! t_ext temperature (K)
! qv_ext specific humidity (kg/kg)
! u_ext u wind component (m/s)
! v_ext v wind component (m/s)
! qc_ext Cloud water mixing ratio (kg/kg)
! qr_ext Rain water mixing ratio (kg/kg)
! qi_ext Ice mixing ratio (kg/kg)
! qs_ext Snow mixing ratio (kg/kg)
! qh_ext Hail mixing ratio (kg/kg)
!
! tsfc_ext Surface temperature
! tdeep_ext Soil temperature
! wetsfc_ext Top layer soil moisture
! wetdp_ext Deep soil moisture
! wetcanp_ext Water content on canopy
!
! trn_ext External terrain (m)
! psfc_ext Surface pressure (Pa)
!
! istatus status indicator
!
! WORK ARRAYS:
!
! var_grb3d Arrays to store the GRIB 3-D variables:
! var_grb3d(nxgrb,nygrb,nzgrb,1,1) - Temperature (K)
! var_grb3d(nxgrb,nygrb,nzgrb,2,1) - Specific humidity
! (kg/kg)
! var_grb3d(nxgrb,nygrb,nzgrb,3,1) - u wind (m/s)
! var_grb3d(nxgrb,nygrb,nzgrb,4,1) - v wind (m/s)
! var_grb3d(nxgrb,nygrb,nzgrb,5,1) - Geopotential
! height (gpm)
! var_grb3d(nxgrb,nygrb,nzgrb,6,1) - Pressure vertical
! velocity (Pa/s)
! (if applied)
! var_grb3d(nxgrb,nygrb,nzgrb,1,2) - soil temp. (K)
! var_grb3d(nxgrb,nygrb,nzgrb,2,2) - vol. soil moist.
! (m**3/m**3)
!
! var_grb2d Arrays to store the GRIB 2-D variables:
! var_grb2d(nxgrb,nygrb,1) - Surface pressure (Pa)
! var_grb2d(nxgrb,nygrb,2) - Geopotential height (gpm)
! var_grb2d(nxgrb,nygrb,3) - Surface temperature (K)
! var_grb2d(nxgrb,nygrb,4) - Plant canopy surface
! water (kg/m**2)
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INCLUDE 'gribcst.inc'
CHARACTER (LEN=*) :: dir_extd
CHARACTER (LEN=*) :: extdname
INTEGER :: extdopt
INTEGER :: extdfmt
CHARACTER (LEN=19) :: extdinit
CHARACTER (LEN=9) :: extdfcst
CHARACTER (LEN=9) :: julfname
!
!-----------------------------------------------------------------------
!
! External grid variables
!
!-----------------------------------------------------------------------
!
INTEGER :: iproj_ext
REAL :: scale_ext,trlon_ext
REAL :: latnot_ext(2)
REAL :: x0_ext,y0_ext
REAL :: dx_ext,dy_ext
!
!-----------------------------------------------------------------------
!
! Output external variable arrays
!
!-----------------------------------------------------------------------
!
INTEGER :: nx_ext,ny_ext,nz_ext
REAL :: lat_ext(nx_ext,ny_ext)
REAL :: lon_ext(nx_ext,ny_ext)
REAL :: p_ext (nx_ext,ny_ext,nz_ext) ! Pressure (Pascals)
REAL :: hgt_ext(nx_ext,ny_ext,nz_ext) ! Height (m)
REAL :: t_ext (nx_ext,ny_ext,nz_ext) ! Temperature (K)
REAL :: qv_ext (nx_ext,ny_ext,nz_ext) ! Specific humidity (kg/kg)
REAL :: u_ext (nx_ext,ny_ext,nz_ext) ! Eastward wind component
REAL :: v_ext (nx_ext,ny_ext,nz_ext) ! Northward wind component
REAL :: qc_ext (nx_ext,ny_ext,nz_ext) ! Cloud H2O mixing ratio (kg/kg)
REAL :: qr_ext (nx_ext,ny_ext,nz_ext) ! Rain H2O mixing ratio (kg/kg)
REAL :: qi_ext (nx_ext,ny_ext,nz_ext) ! Ice mixing ratio (kg/kg)
REAL :: qs_ext (nx_ext,ny_ext,nz_ext) ! Snow mixing ratio (kg/kg)
REAL :: qh_ext (nx_ext,ny_ext,nz_ext) ! Hail mixing ratio (kg/kg)
REAL :: tsfc_ext (nx_ext,ny_ext) ! Temperature at surface (K)
REAL :: tdeep_ext (nx_ext,ny_ext) ! Deep soil temperature (K)
REAL :: wetsfc_ext (nx_ext,ny_ext) ! Surface soil moisture
REAL :: wetdp_ext (nx_ext,ny_ext) ! Deep soil moisture
REAL :: wetcanp_ext(nx_ext,ny_ext) ! Canopy water amount
REAL :: snowdpth_ext(nx_ext,ny_ext) ! Snow depth (m)
REAL :: trn_ext (nx_ext,ny_ext) ! External terrain (m)
REAL :: psfc_ext (nx_ext,ny_ext) ! Surface pressure (Pa)
!
!-----------------------------------------------------------------------
!
! Other external variable arrays
!
!-----------------------------------------------------------------------
!
REAL :: x_ext(nx_ext)
REAL :: y_ext(ny_ext)
INTEGER :: istatus
!
!-----------------------------------------------------------------------
!
! Work arrays for storing grib data
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: var_grb2d(:,:,:,:) ! GRIB variables
REAL, ALLOCATABLE :: var_grb3d(:,:,:,:,:) ! GRIB 3-D variables
INTEGER, ALLOCATABLE :: var_lev3d(:,:,:) ! Levels (hybrid) for
! each 3-D variable
REAL, ALLOCATABLE :: rcdata(:) ! temporary data array
!
!-----------------------------------------------------------------------
!
! Original grid variables
!
!-----------------------------------------------------------------------
!
INTEGER :: iproj
REAL :: scale,trlon,x0,y0
REAL :: latnot(2)
!
!-----------------------------------------------------------------------
!
! Misc internal variables
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=80) :: gribfile
CHARACTER (LEN=13) :: gribtime
INTEGER :: i,j,k,kk
INTEGER :: iyr,imo,iday,myr, jldy
INTEGER :: ihr,imin,isec
INTEGER :: ifhr,ifmin,ifsec
INTEGER :: grbflen, len1, lenrun
INTEGER :: m,n,nz1,max_nr2d,max_nr3d,min_nr3d,nz2
REAL :: govrd
INTEGER :: chklev, lvscan
INTEGER :: iret ! Return flag
!
!-----------------------------------------------------------------------
!
! Include files
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
! GRIB grid information
!
!-----------------------------------------------------------------------
!
CHARACTER (LEN=42) :: gridesc ! Grid description
INTEGER :: iproj_grb ! Map projection indicator
INTEGER :: gthin ! Indicator of whether the grid is "thinned"
INTEGER :: ni_grb ! Number of points along x-axis
INTEGER :: nj_grb ! Number of points along y-axis
INTEGER :: np_grb ! Total number of horizontal grid points
INTEGER :: nk_grb ! Number of vertical parameters
REAL :: zk_grb(nz_ext) ! Vertical coordinate parameters
INTEGER :: npeq ! Number of lat circles from pole to equator
INTEGER :: nit(nz_ext) ! Number of x-points for thinned grid
REAL :: pi_grb ! x-coordinate of pole point
REAL :: pj_grb ! y-coordinate of pole point
INTEGER :: ipole ! Projection center flag
REAL :: di_grb ! x-direction increment or grid length
REAL :: dj_grb ! y-direction increment or grid length
REAL :: latsw ! Latitude of South West corner point
REAL :: lonsw ! Longitude of South West corner point
REAL :: latne ! Latitude of North East corner point
REAL :: lonne ! Longitude of North East corner point
REAL :: lattru1 ! Latitude (1st) at which projection is true
REAL :: lattru2 ! Latitude (2nd) at which projection is true
REAL :: lontrue ! Longitude at which projection is true
REAL :: latrot ! Latitude of southern pole of rotation
REAL :: lonrot ! Longitude of southern pole of rotation
REAL :: angrot ! Angle of rotation
REAL :: latstr ! Latitude of the pole of stretching
REAL :: lonstr ! Longitude of the pole of stretching
REAL :: facstr ! Stretching factor
INTEGER :: scanmode ! Scanning indicator
INTEGER :: iscan ! x-direction scanning indicator
INTEGER :: jscan ! y-direction scanning indicator
INTEGER :: kscan ! FORTRAN index scanning indicator
INTEGER :: ires ! Resolution direction increments indicator
INTEGER :: iearth ! Earth shape indicator: spherical or oblate?
INTEGER :: icomp ! (u,v) components decomposition indicator
INTEGER :: jpenta ! J-Pentagonal resolution parameter
INTEGER :: kpenta ! K-Pentagonal resolution parameter
INTEGER :: mpenta ! M-Pentagonal resolution parameter
INTEGER :: ispect ! Spectral representation type
INTEGER :: icoeff ! Spectral coefficient storage mode
REAL :: xp_grb ! X coordinate of sub-satellite point
REAL :: yp_grb ! Y coordinate of sub-satellite point
REAL :: xo_grb ! X coordinate of image sector origin
REAL :: yo_grb ! Y coordinate of image sector origin
REAL :: zo_grb ! Camera altitude from center of Earth
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
ALLOCATE(var_grb2d(nx_ext,ny_ext,n2dvs,n2dlvt))
ALLOCATE(var_grb3d(nx_ext,ny_ext,nz_ext,n3dvs,n3dlvt))
ALLOCATE(rcdata(nx_ext*ny_ext))
ALLOCATE(var_lev3d(nz_ext,n3dvs,n3dlvt))
IF(extdfcst == ' ') extdfcst='000:00:00'
READ (extdinit,'(i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,i2)') &
iyr,imo,iday,ihr,imin,isec
CALL julday(iyr,imo,iday,jldy)
myr=MOD(iyr,100)
ifhr=0
ifmin=0
ifsec=0
READ(extdfcst,'(i3)',ERR=4,END=4) ifhr
4 CONTINUE
WRITE (gribtime,'(i4.4,i2.2,i2.2,i2.2,a1,i2.2)') &
iyr,imo,iday,ihr,'f',ifhr
len1=LEN(dir_extd)
grbflen=len1
CALL strlnth( dir_extd, grbflen )
IF( grbflen == 0 .OR. dir_extd(1:grbflen) == ' ' ) THEN
dir_extd = '.'
grbflen=1
END IF
IF( dir_extd(grbflen:grbflen) /= '/'.AND.grbflen < len1 ) THEN
grbflen=grbflen+1
dir_extd(grbflen:grbflen)='/'
END IF
lenrun = LEN( extdname )
CALL strlnth( extdname, lenrun )
gribfile = dir_extd(1:grbflen)//extdname(1:lenrun) &
//'.'//gribtime(1:13)
grbflen = grbflen + lenrun + 14
!
!-----------------------------------------------------------------------
!
! RDNMCGRB reads NMC GRIB data
!
!-----------------------------------------------------------------------
!
gridtyp = avn3grid
mproj_grb = avn3proj
n2dvars = avn3nvs2d
n2dlvtps = avn3nlvt2d
DO k=1,n2dlvtps
DO n=1,n2dvars
var_id2d(n,k) = avn3var_id2d(n,k)
END DO
levtyp2d(k) = avn3levs2d(k)
END DO
n3dvars = avn3nvs3d
n3dlvtps = avn3nlvt3d
DO m=1,n3dlvtps
DO n=1,n3dvars
var_id3d(n,m) = avn3var_id3d(n,m)
END DO
levtyp3d(m) = avn3levs3d(m)
END DO
CALL rdnmcgrb(nx_ext,ny_ext,nz_ext,gribfile,grbflen, gribtime, &
gridesc, iproj_grb, gthin, &
ni_grb,nj_grb,np_grb, nk_grb,zk_grb, npeq,nit, &
pi_grb,pj_grb,ipole, di_grb,dj_grb, &
latsw,lonsw, latne,lonne, &
latrot,lonrot,angrot, &
latstr,lonstr,facstr, &
lattru1,lattru2,lontrue, &
scanmode, iscan,jscan,kscan, &
ires,iearth,icomp, &
jpenta,kpenta,mpenta,ispect,icoeff, &
xp_grb,yp_grb, xo_grb,yo_grb,zo_grb, &
rcdata,var_grb2d,var_grb3d,var_lev3d,iret)
max_nr2d = 0
DO n=1,n2dvars
DO m=1,n2dlvtps
max_nr2d = MAX( max_nr2d, var_nr2d(n,m) )
END DO
END DO
max_nr3d = 0
min_nr3d = nz_ext
DO n=1,n3dvars
max_nr3d = MAX( max_nr3d, var_nr3d(n,1) )
min_nr3d = MIN( min_nr3d, var_nr3d(n,1) )
END DO
IF ( max_nr3d == 0 ) THEN
WRITE (6,'(a)') &
'No 3-D variable was found in the GRIB file', &
'Program stopped in GETNCEPAVN3.'
STOP
END IF
IF ( max_nr2d == 0 ) THEN
WRITE (6,'(a)') &
'No 2-D variables was found in the GRIB file'
END IF
! write (6,'(/a7,2x,6(i7))')
! : 'Lev\\VID',(var_id3d(n,1),n=1,n3dvars)
! DO 60 k=1,max_nr3d
! var_lev3d(k,5,1) = var_lev3d(k,1,1)
! var_lev3d(k,6,1) = var_lev3d(k,1,1)
! write (6,'(/i5,4x,6(i7))')
! : k,(var_lev3d(k,n,1),n=1,n3dvars)
60 CONTINUE
DO k=1,min_nr3d
DO n=2,n3dvars
IF ( var_lev3d(k,1,1) /= var_lev3d(k,n,1) ) THEN
WRITE (6,'(a)') &
'Variables were not at the same level.', &
'Program stopped in GETNCEPAVN3.'
STOP
END IF
END DO
END DO
IF ( iproj_grb == 5 .AND. ipole == 0 ) THEN ! Center in NP
iproj_ext = 1
ELSE IF ( iproj_grb == 5 .AND. ipole == 1 ) THEN ! Center in SP
iproj_ext = -1
ELSE IF ( iproj_grb == 3 .AND. ipole == 0 ) THEN ! Center in NP
iproj_ext = 2
ELSE IF ( iproj_grb == 3 .AND. ipole == 1 ) THEN ! Center in SP
iproj_ext = -2
ELSE IF ( iproj_grb == 1 ) THEN
iproj_ext = 3
ELSE IF ( iproj_grb == 0 ) THEN
iproj_ext = 4
ELSE
WRITE (6,'(a)') &
'Unknown map projection. Set to non-projection.'
iproj_ext = 0
END IF
scale_ext = 1.0
latnot_ext(1) = lattru1
latnot_ext(2) = lattru2
trlon_ext = lontrue
dx_ext = di_grb
dy_ext = dj_grb
DO i=1, nx_ext
DO j=1, ny_ext
lon_ext(i,j)= lonsw + (i-1) * dx_ext
lat_ext(i,j)= latsw + (j-1) * dy_ext
END DO
END DO
PRINT *,'LatSW = ',latsw,' LonSW = ',lonsw
!
!-----------------------------------------------------------------------
!
! Retrieve 2-D variables
!
!-----------------------------------------------------------------------
!
DO j=1,ny_ext
DO i=1,nx_ext
IF ( var_nr2d(1,1) == 0 ) THEN
psfc_ext (i,j) = -999.0
ELSE
psfc_ext (i,j) = var_grb2d(i,j,1,1) * 100.0
END IF
IF ( var_nr2d(2,1) == 0 ) THEN
trn_ext (i,j) = -999.0
ELSE
trn_ext (i,j) = var_grb2d(i,j,2,1)/g
END IF
IF ( var_nr3d(1,2) == 0 ) THEN
tsfc_ext (i,j) = -999.0
tdeep_ext (i,j) = -999.0
wetsfc_ext(i,j) = -999.0
wetdp_ext (i,j) = -999.0
ELSE
tsfc_ext (i,j) = var_grb2d(i,j,3,1) ! sfc temp.
IF ( nint(var_grb2d(i,j,5,1)) == 1 ) THEN ! soil temp over land
tdeep_ext (i,j) = var_grb3d(i,j,1,1,2)
IF ( tdeep_ext (i,j) <= 200. ) THEN
tdeep_ext (i,j) = tsfc_ext(i,j)
END IF
wetsfc_ext(i,j) = var_grb3d(i,j,2,2,2)
wetdp_ext(i,j) = var_grb3d(i,j,1,2,2)
ELSE ! sfc temp over sea
tdeep_ext (i,j) = tsfc_ext(i,j)
wetsfc_ext(i,j) = 1.0
wetdp_ext(i,j) = 1.0
END IF
END IF
IF ( var_nr2d(4,1) == 0 ) THEN
wetcanp_ext(i,j) = -999.0
ELSE
wetcanp_ext(i,j) = var_grb2d(i,j,4,1)*1.e-3 ! in meter
END IF
IF ( var_nr2d(6,1) == 0 ) THEN
snowdpth_ext(i,j) = -999.
ELSE
! Convert water equiv. of accum. snow depth (kg/m**2) to meters
! (where 1 meter liquid water is set equivqlent to 10 meters snow).
! 0.01 = 10. (m snow/m liquid) / (1000 kg/m**3)
snowdpth_ext(i,j) = 0.01 * var_grb2d(i,j,6,1) ! in meters
END IF
END DO
END DO
!
!-----------------------------------------------------------------------
!
! Retrieve 3-D variables
!
!-----------------------------------------------------------------------
!
nz1 = MIN(var_nr3d(1,1),nz_ext)
IF ( var_lev3d(1,1,1) > var_lev3d(nz1,1,1) ) THEN ! 1st level at sfc
chklev = 1
lvscan = 0
ELSE
chklev = -1
lvscan = nz1+1
END IF
DO k=1,nz1
kk = chklev * k + lvscan
DO j=1,ny_ext
DO i=1,nx_ext
p_ext (i,j,kk) = 100.0 * FLOAT(var_lev3d(k,1,1)) ! Pressure
hgt_ext(i,j,kk) = var_grb3d(i,j,k,1,1)
u_ext (i,j,kk) = var_grb3d(i,j,k,2,1) ! u wind (m/s)
v_ext (i,j,kk) = var_grb3d(i,j,k,3,1) ! v wind (m/s)
t_ext (i,j,kk) = var_grb3d(i,j,k,4,1) ! Temperature (K)
qc_ext (i,j,kk) = -999.
qr_ext (i,j,kk) = -999.
qi_ext (i,j,kk) = -999.
qs_ext (i,j,kk) = -999.
qh_ext (i,j,kk) = -999.
END DO
END DO
END DO
CALL getqvs(nx_ext,ny_ext,nz1, 1,nx_ext,1,ny_ext,1,nz1, &
p_ext, t_ext, qv_ext )
nz2 = MIN( nz1, min_nr3d )
DO k=1,nz2
kk = chklev * k + lvscan
DO j=1,ny_ext
DO i=1,nx_ext
qv_ext(i,j,kk) = 0.01*var_grb3d(i,j,k,5,1)*qv_ext(i,j,kk)
END DO
END DO
END DO
IF ( nz2 < nz1 ) THEN
DO k=nz2+1,nz1
kk = chklev * k + lvscan
DO j=1,ny_ext
DO i=1,nx_ext
qv_ext(i,j,kk) = 0.0
var_lev3d(k,5,1) = var_lev3d(k,1,1)
END DO
END DO
END DO
END IF
istatus = 1
DEALLOCATE(var_grb2d,var_grb3d,rcdata,var_lev3d)
RETURN
END SUBROUTINE getncepavn3
!
!##################################################################
!##################################################################
!###### ######
!###### SUBROUTINE GET_AVN_BIN ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
SUBROUTINE get_avn_bin(nx_ext,ny_ext,nz_ext, & 1,4
dir_extd,extdname,extdinit,extdfcst,julfname, &
iproj_ext,scale_ext,trlon_ext,latnot_ext,x0_ext,y0_ext, &
lat_ext,lon_ext,p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext, &
qc_ext,qr_ext,qi_ext,qs_ext,qh_ext, &
tsfc_ext,tdeep_ext,wetsfc_ext,wetdp_ext,wetcanp_ext, &
snowdpth_ext,trn_ext,psfc_ext, &
istatus)
!
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Reads and pass out a section of NCEP AVN GRIB
! (Grid #3, 1x1 degree) data file (created by program EXTRACT_AVN).
!
!-----------------------------------------------------------------------
!
! AUTHOR: M. Xue
! 07/25/2000
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
! INPUT:
!
! dir_extd Directory name for external file
! extdname Prefix string of external file name
! extdinit Initialized time in mm-dd-yyyy:hh:mm:ss format
! extdfcst Forecast hour in HHH:MM:SS format
! julfname File name in yyjjjhhmm format
!
! OUTPUT:
!
! iproj_ext Map projection number of external data
! scale_ext Scale factor of external data
! trlon_ext True longitude of external data (degrees E)
! latnot_ext(2) True latitude(s) of external data (degrees N)
! x0_ext x coordinate of origin of external data
! y0_ext y coordinate of origin of external data
! lat_ext latitude of external data points (degrees N)
! lon_ext longitude of external data points (degrees E)
! p_ext pressure (Pascal)
! hgt_ext height (m)
! t_ext temperature (K)
! qv_ext specific humidity (kg/kg)
! u_ext u wind component (m/s)
! v_ext v wind component (m/s)
! qc_ext Cloud water mixing ratio (kg/kg)
! qr_ext Rain water mixing ratio (kg/kg)
! qi_ext Ice mixing ratio (kg/kg)
! qs_ext Snow mixing ratio (kg/kg)
! qh_ext Hail mixing ratio (kg/kg)
!
! tsfc_ext Surface temperature
! tdeep_ext Soil temperature
! wetsfc_ext Top layer soil moisture
! wetdp_ext Deep soil moisture
! wetcanp_ext Water content on canopy
!
! trn_ext External terrain (m)
! psfc_ext Surface pressure (Pa)
!
! istatus status indicator
!
! WORK ARRAYS:
!
!-----------------------------------------------------------------------
!
! Variable Declarations.
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
INTEGER :: nx_ext,ny_ext,nz_ext
CHARACTER (LEN=*) :: dir_extd
CHARACTER (LEN=*) :: extdname
CHARACTER (LEN=19) :: extdinit
CHARACTER (LEN=9) :: extdfcst
CHARACTER (LEN=9) :: julfname
!
!-----------------------------------------------------------------------
!
! External grid variables
!
!-----------------------------------------------------------------------
!
INTEGER :: iproj_ext
REAL :: scale_ext,trlon_ext
REAL :: latnot_ext(2)
REAL :: x0_ext,y0_ext
REAL :: dx_ext,dy_ext
!
!-----------------------------------------------------------------------
!
! Output external variable arrays
!
!-----------------------------------------------------------------------
!
REAL :: lat_ext(nx_ext,ny_ext)
REAL :: lon_ext(nx_ext,ny_ext)
REAL :: p_ext (nx_ext,ny_ext,nz_ext) ! Pressure (Pascals)
REAL :: hgt_ext(nx_ext,ny_ext,nz_ext) ! Height (m)
REAL :: t_ext (nx_ext,ny_ext,nz_ext) ! Temperature (K)
REAL :: qv_ext (nx_ext,ny_ext,nz_ext) ! Specific humidity (kg/kg)
REAL :: u_ext (nx_ext,ny_ext,nz_ext) ! Eastward wind component
REAL :: v_ext (nx_ext,ny_ext,nz_ext) ! Northward wind component
REAL :: qc_ext (nx_ext,ny_ext,nz_ext) ! Cloud H2O mixing ratio (kg/kg)
REAL :: qr_ext (nx_ext,ny_ext,nz_ext) ! Rain H2O mixing ratio (kg/kg)
REAL :: qi_ext (nx_ext,ny_ext,nz_ext) ! Ice mixing ratio (kg/kg)
REAL :: qs_ext (nx_ext,ny_ext,nz_ext) ! Snow mixing ratio (kg/kg)
REAL :: qh_ext (nx_ext,ny_ext,nz_ext) ! Hail mixing ratio (kg/kg)
REAL :: tsfc_ext (nx_ext,ny_ext) ! Temperature at surface (K)
REAL :: tdeep_ext (nx_ext,ny_ext) ! Deep soil temperature (K)
REAL :: wetsfc_ext (nx_ext,ny_ext) ! Surface soil moisture
REAL :: wetdp_ext (nx_ext,ny_ext) ! Deep soil moisture
REAL :: wetcanp_ext(nx_ext,ny_ext) ! Canopy water amount
REAL :: snowdpth_ext(nx_ext,ny_ext) ! Snow depth (m)
REAL :: trn_ext (nx_ext,ny_ext) ! External terrain (m)
REAL :: psfc_ext (nx_ext,ny_ext) ! Surface pressure (Pa)
REAL :: temp_ext (nx_ext,ny_ext) ! Automatic work array
!-----------------------------------------------------------------------
!
! Local variables
!
!-----------------------------------------------------------------------
INTEGER :: istatus,ierr
INTEGER :: nunit, idummy
REAL :: rdummy
CHARACTER (LEN=10) :: var_label
CHARACTER (LEN=80) :: gribfile
CHARACTER (LEN=13) :: gribtime
INTEGER :: i,j,k
INTEGER :: iyr,imo,iday,myr, jldy
INTEGER :: ihr,imin,isec
INTEGER :: ifhr,ifmin,ifsec
INTEGER :: grbflen, len1, lenrun
INTEGER :: iuout, ivout, ipout, ihout,itout, &
iqvout, itsfcout,itsoilout,iwsfcout,iwdpout, &
iwcnpout,isnowout,itrnout,ipsfcout, &
iugrdout,ivgrdout,itgrdout,iptgrdout,irhgrdout,ipmslout
INTEGER :: nx_ext_in, ny_ext_in, nz_ext_in
REAL :: latbgn,latend,lonbgn,lonend, del_lat, del_lon
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
IF(extdfcst == ' ') extdfcst='000:00:00'
READ (extdinit,'(i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,i2)') &
iyr,imo,iday,ihr,imin,isec
CALL julday(iyr,imo,iday,jldy)
myr=MOD(iyr,100)
ifhr=0
ifmin=0
ifsec=0
READ(extdfcst,'(i3)',ERR=4,END=4) ifhr
4 CONTINUE
WRITE (gribtime,'(i4.4,i2.2,i2.2,i2.2,a1,i2.2)') &
iyr,imo,iday,ihr,'f',ifhr
len1=len_trim(dir_extd)
grbflen=len1
IF( grbflen == 0 .OR. dir_extd(1:grbflen) == ' ' ) THEN
dir_extd = '.'
grbflen=1
END IF
IF( dir_extd(grbflen:grbflen) /= '/'.AND.grbflen < len1 ) THEN
grbflen=grbflen+1
dir_extd(grbflen:grbflen)='/'
END IF
lenrun = len_trim( extdname )
gribfile = dir_extd(1:grbflen)//extdname(1:lenrun)//'.'//gribtime(1:13)
grbflen = grbflen + lenrun + 14
psfc_ext = -999.0
trn_ext = -999.0
tsfc_ext = -999.0
tdeep_ext = -999.0
wetsfc_ext = -999.0
wetdp_ext = -999.0
wetcanp_ext= -999.0
snowdpth_ext=-999.0
u_ext = -999.0
v_ext = -999.0
p_ext = -999.0
hgt_ext = -999.0
t_ext = -999.0
qv_ext = -999.0
CALL getunit( nunit)
CALL asnctl ('NEWLOCAL', 1, ierr)
CALL asnfile(trim(gribfile), '-F f77 -N ieee', ierr)
OPEN(UNIT=nunit,FILE=trim(gribfile), &
STATUS='old',FORM='unformatted',IOSTAT=istatus)
IF( istatus /= 0 ) THEN
WRITE(6,'(1x,a,a,/1x,i3,a)') &
'Error occured when opening file ',trim(gribfile), &
'using FORTRAN unit ',nunit,' Program stopped.'
STOP
END IF
PRINT*,'To read file ',trim(gribfile)
READ (nunit,ERR=910,END=920) nx_ext_in,ny_ext_in, nz_ext_in
IF( nx_ext_in /= nx_ext.OR.ny_ext_in /= ny_ext .OR. nz_ext_in /= nz_ext ) THEN
WRITE(6,'(1x,a)') &
' Dimensions in GET_AVN_BIN inconsistent with data.'
WRITE(6,'(1x,a,3I15)') ' Read were: ', &
nx_ext_in, ny_ext_in, nz_ext_in
WRITE(6,'(1x,a,3I15)') ' Expected: ', &
nx_ext, ny_ext, nz_ext
WRITE(6,'(1x,a)') &
' Program aborted in GET_AVN_BIN.'
STOP
END IF
READ (nunit,ERR=910,END=920) lonbgn,lonend,latbgn,latend
READ (nunit,ERR=910,END=920) del_lon, del_lat
READ (nunit,ERR=910,END=920) &
iuout, ivout, ipout, ihout,itout, &
iqvout, itsfcout,itsoilout,iwsfcout,iwdpout, &
iwcnpout,isnowout,itrnout,ipsfcout,iugrdout, &
ivgrdout,itgrdout,iptgrdout,irhgrdout,ipmslout, &
iproj_ext,idummy,idummy, idummy, idummy, &
idummy, idummy, idummy, idummy, idummy
READ (nunit,ERR=910,END=920) &
scale_ext,trlon_ext,latnot_ext(1),latnot_ext(2), &
x0_ext, y0_ext, rdummy, rdummy, rdummy, rdummy, &
rdummy, rdummy, rdummy, rdummy, rdummy, &
rdummy, rdummy, rdummy, rdummy, rdummy, &
rdummy, rdummy, rdummy, rdummy, rdummy, &
rdummy, rdummy, rdummy, rdummy, rdummy
IF( iuout == 1 ) THEN
READ (nunit,ERR=910,END=920) var_label
READ (nunit,ERR=910,END=920) u_ext
WRITE(6,'(a,a,a)')'Read ',var_label,' into array u_ext.'
END IF
IF( ivout == 1 ) THEN
READ (nunit,ERR=910,END=920) var_label
READ (nunit,ERR=910,END=920) v_ext
WRITE(6,'(a,a,a)')'Read ',var_label,' into array v_ext.'
END IF
IF( ipout == 1 ) THEN
READ (nunit,ERR=910,END=920) var_label
READ (nunit,ERR=910,END=920) p_ext
WRITE(6,'(a,a,a)')'Read ',var_label,' into array p_ext.'
END IF
IF( ihout == 1 ) THEN
READ (nunit,ERR=910,END=920) var_label
READ (nunit,ERR=910,END=920) hgt_ext
WRITE(6,'(a,a,a)')'Read ',var_label,' into array hgt_ext.'
END IF
IF( itout == 1 ) THEN
READ (nunit,ERR=910,END=920) var_label
READ (nunit,ERR=910,END=920) t_ext
WRITE(6,'(a,a,a)')'Read ',var_label,' into array t_ext.'
END IF
IF( iqvout== 1 ) THEN
READ (nunit,ERR=910,END=920) var_label
READ (nunit,ERR=910,END=920) qv_ext
WRITE(6,'(a,a,a)')'Read ',var_label,' into array qv_ext.'
END IF
IF( itsfcout==1 ) THEN
READ (nunit,ERR=910,END=920) var_label
READ (nunit,ERR=910,END=920) tsfc_ext
WRITE(6,'(a,a,a)')'Read ',var_label,' into array tsfc_ext.'
END IF
IF( itsoilout==1) THEN
READ (nunit,ERR=910,END=920) var_label
READ (nunit,ERR=910,END=920) tdeep_ext
WRITE(6,'(a,a,a)')'Read ',var_label,' into array tdeep_ext.'
END IF
IF( iwsfcout ==1) THEN
READ (nunit,ERR=910,END=920) var_label
READ (nunit,ERR=910,END=920) wetsfc_ext
WRITE(6,'(a,a,a)')'Read ',var_label,' into array wetsfc_ext.'
END IF
IF( iwdpout ==1 ) THEN
READ (nunit,ERR=910,END=920) var_label
READ (nunit,ERR=910,END=920) wetdp_ext
WRITE(6,'(a,a,a)')'Read ',var_label,' into array wetdp_ext.'
END IF
IF( iwcnpout ==1) THEN
READ (nunit,ERR=910,END=920) var_label
READ (nunit,ERR=910,END=920) wetcanp_ext
WRITE(6,'(a,a,a)')'Read ',var_label,' into array wetcanp_ext.'
END IF
IF( isnowout ==1) THEN
READ (nunit,ERR=910,END=920) var_label
READ (nunit,ERR=910,END=920) snowdpth_ext
WRITE(6,'(a,a,a)')'Read ',var_label,' into array snowdpth_ext.'
END IF
IF( itrnout ==1 ) THEN
READ (nunit,ERR=910,END=920) var_label
READ (nunit,ERR=910,END=920) trn_ext
WRITE(6,'(a,a,a)')'Read ',var_label,' into array trn_ext.'
END IF
IF( ipsfcout ==1) THEN
READ (nunit,ERR=910,END=920) var_label
READ (nunit,ERR=910,END=920) psfc_ext
WRITE(6,'(a,a,a)')'Read ',var_label,' into array psfc_ext.'
END IF
IF( iugrdout ==1) THEN
READ (nunit,ERR=910,END=920) var_label
READ (nunit,ERR=910,END=920) temp_ext ! discarded
WRITE(6,'(a,a,a)')'Read ',var_label,' into array temp_ext.'
END IF
IF( ivgrdout ==1) THEN
READ (nunit,ERR=910,END=920) var_label
READ (nunit,ERR=910,END=920) temp_ext ! discarded
WRITE(6,'(a,a,a)')'Read ',var_label,' into array temp_ext.'
END IF
IF( itgrdout ==1) THEN
READ (nunit,ERR=910,END=920) var_label
READ (nunit,ERR=910,END=920) temp_ext ! discarded
WRITE(6,'(a,a,a)')'Read ',var_label,' into array temp_ext.'
END IF
IF( irhgrdout ==1) THEN
READ (nunit,ERR=910,END=920) var_label
READ (nunit,ERR=910,END=920) temp_ext ! discarded
WRITE(6,'(a,a,a)')'Read ',var_label,' into array temp_ext.'
END IF
IF( iptgrdout ==1) THEN
READ (nunit,ERR=910,END=920) var_label
READ (nunit,ERR=910,END=920) temp_ext ! discarded
WRITE(6,'(a,a,a)')'Read ',var_label,' into array temp_ext.'
END IF
IF( ipmslout ==1) THEN
READ (nunit,ERR=910,END=920) var_label
READ (nunit,ERR=910,END=920) temp_ext ! discarded
WRITE(6,'(a,a,a)')'Read ',var_label,' into array temp_ext.'
END IF
CLOSE (UNIT=nunit)
CALL retunit( nunit)
PRINT*,'Finished reading file ',trim(gribfile)
PRINT*,' '
900 CONTINUE
dx_ext = del_lon
dy_ext = del_lat
IF( lonend > 180.0 ) THEN
lonend = lonend - 360.0
END IF
IF( lonbgn > 180.0 ) THEN
lonbgn = lonbgn - 360.0
END IF
IF( lonbgn > lonend ) lonbgn = lonbgn - 360.0
DO i=1,nx_ext
DO j=1,ny_ext
lon_ext(i,j)= lonbgn + (i-1) * dx_ext
lat_ext(i,j)= latbgn + (j-1) * dy_ext - 90.0
END DO
END DO
istatus = 1
RETURN
!
!-----------------------------------------------------------------------
!
! Error during read
!
!----------------------------------------------------------------------
!
910 CONTINUE
WRITE(6,'(/a/)') ' Error reading data in GET_AVN_BIN.'
istatus =2
RETURN
!
!-----------------------------------------------------------------------
!
! End-of-file during read
!
!----------------------------------------------------------------------
!
920 CONTINUE
WRITE(6,'(/a/)') ' End of file reached in GET_AVN_BIN.'
istatus =3
RETURN
END SUBROUTINE get_avn_bin