! !################################################################## !################################################################## !###### ###### !###### 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,82 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, & snowdpth_ext,ubar_ext,vbar_ext,wbar_ext, & ptbar_ext,pbar_ext,qvbar_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 ! ! snowdpth_ext Snow depth (m) ! ubar_ext Base state u-velocity (m/s) ! vbar_ext Base state v-velocity (m/s) ! wbar_ext Base state w-velocity (m/s) ! ! ptbar_ext Base state potential temperature (K) ! pbar_ext Base state pressure (Pascal) ! qvbar_ext Base state water vapor specific humidity (kg/kg) ! ! 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 :: 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 :: 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 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=256) :: grdbasfn_ext CHARACTER (LEN=256) :: datafn_ext INTEGER :: nchanl_ext,lengbf_ext INTEGER :: lendtf_ext INTEGER :: mp_opt_save REAL :: time_ext ! !----------------------------------------------------------------------- ! ! Include files ! !----------------------------------------------------------------------- ! INCLUDE 'globcst.inc' INCLUDE 'grid.inc' INCLUDE 'phycst.inc' INCLUDE 'mp.inc' ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! !----------------------------------------------------------------------- ! ! Build file names ! !----------------------------------------------------------------------- ! IF (myproc == 0 ) THEN 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 ! ! ! We need to disable MPI here, as we DON'T want the data scattered to the ! individual processors! The current code calls for the entire external ! data to be passed to all processors. The assumption is that external ! data is much lower resolution, so this can be accomplished. Eventually, ! this will have to be changed, but not at this release. ! mp_opt_save = mp_opt mp_opt = 0 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 ! ! Restore the correct "mp_opt" for processor 0. ! mp_opt = mp_opt_save END IF ! !----------------------------------------------------------------------- ! ! Reset map projection to previous values ! !----------------------------------------------------------------------- ! CALL mpupdatei(mapproj,1) CALL mpupdater(trulat1,1) CALL mpupdater(trulat2,1) CALL mpupdater(trulon, 1) CALL mpupdater(sclfct,1) CALL mpupdater(lat_org,1) CALL mpupdater(lon_org,1) latnot(1)=trulat1 latnot(2)=trulat2 CALL setmapr(mapproj,sclfct,latnot,trulon) CALL setorig(2,lat_org,lon_org) CALL xytoll(1,1,0.,0.,lat_org,lon_org) CALL mpupdatei(iproj_ext,1) CALL mpupdater(scale_ext,1) CALL mpupdater(latnot_ext,2) CALL mpupdater(trlon_ext,1) CALL mpupdater(x0_ext,1) CALL mpupdater(y0_ext,1) CALL mpupdater(lat_ext,nx_ext*ny_ext) CALL mpupdater(lon_ext,nx_ext*ny_ext) CALL mpupdater(lat_ext,nx_ext*ny_ext) CALL mpupdater(lon_ext,nx_ext*ny_ext) CALL mpupdater(latu_ext,nx_ext*ny_ext) CALL mpupdater(lonu_ext,nx_ext*ny_ext) CALL mpupdater(latv_ext,nx_ext*ny_ext) CALL mpupdater(lonv_ext,nx_ext*ny_ext) CALL mpupdater(p_ext,nx_ext*ny_ext*nz_ext) CALL mpupdater(hgt_ext,nx_ext*ny_ext*nz_ext) CALL mpupdater(zp_ext,nx_ext*ny_ext*nz_ext) CALL mpupdater(zpsoil_ext,nx_ext*ny_ext*nz_ext) CALL mpupdater(t_ext,nx_ext*ny_ext*nz_ext) CALL mpupdater(qv_ext,nx_ext*ny_ext*nz_ext) CALL mpupdater(u_ext,nx_ext*ny_ext*nz_ext) CALL mpupdater(v_ext,nx_ext*ny_ext*nz_ext) CALL mpupdater(uatv_ext,nx_ext*ny_ext*nz_ext) CALL mpupdater(vatu_ext,nx_ext*ny_ext*nz_ext) CALL mpupdater(w_ext,nx_ext*ny_ext*nz_ext) CALL mpupdater(qc_ext,nx_ext*ny_ext*nz_ext) CALL mpupdater(qr_ext,nx_ext*ny_ext*nz_ext) CALL mpupdater(qi_ext,nx_ext*ny_ext*nz_ext) CALL mpupdater(qs_ext,nx_ext*ny_ext*nz_ext) CALL mpupdater(qh_ext,nx_ext*ny_ext*nz_ext) CALL mpupdatei(soiltyp_ext,nx_ext*ny_ext*nstyps) CALL mpupdater(stypfrct_ext,nx_ext*ny_ext*nstyps) CALL mpupdatei(vegtyp_ext,nx_ext*ny_ext) CALL mpupdater(lai_ext,nx_ext*ny_ext) CALL mpupdater(roufns_ext,nx_ext*ny_ext) CALL mpupdater(veg_ext,nx_ext*ny_ext) CALL mpupdater(tsoil_ext,nx_ext*ny_ext*nzsoil_ext*(nstyps+1)) CALL mpupdater(qsoil_ext,nx_ext*ny_ext*nzsoil_ext*(nstyps+1)) CALL mpupdater(wetcanp_ext,nx_ext*ny_ext*(nstyps+1)) CALL mpupdater(snowdpth_ext,nx_ext*ny_ext) CALL mpupdater(tem1_ext,nx_ext*ny_ext*nz_ext) CALL mpupdater(tem2_ext,nx_ext*ny_ext*nz_ext) CALL mpupdater(x_ext,nx_ext) CALL mpupdater(y_ext,ny_ext) CALL mpupdater(z_ext,nz_ext) CALL mpupdater(xsc,nx_ext) CALL mpupdater(ysc,ny_ext) CALL mpupdater(ubar_ext,nx_ext*ny_ext*nz_ext) CALL mpupdater(vbar_ext,nx_ext*ny_ext*nz_ext) CALL mpupdater(wbar_ext,nx_ext*ny_ext*nz_ext) CALL mpupdater(ptbar_ext,nx_ext*ny_ext*nz_ext) CALL mpupdater(pbar_ext,nx_ext*ny_ext*nz_ext) CALL mpupdater(qvbar_ext,nx_ext*ny_ext*nz_ext) 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=256) :: 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=256) :: 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=256) :: 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.) mixrat= qv_ext(i,j,kk) / ( 1. - qv_ext(i,j,kk) ) 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) mixrat= qvavg / (1. - qvavg) 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, & t_2m_ext, qv_2m_ext, u_10m_ext, v_10m_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. ! ! 2003/10/23 (F. KONG) ! Modified to read in from grid 212 four more near surface variables: ! 2-m temperature and humidity, and 10-m wind (u, v) ! !----------------------------------------------------------------------- ! ! 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) ! ! F.KONG add four near surface variables ! T_2m_ext Temperature at 2m AGL ! qv_2m_ext Specific Humidity at 2m AGL ! U_10m_ext U at 10m AGL ! V_10m_ext V at 10m AGL ! end F.KONG mod ! ! 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) ! Soil type REAL :: t_2m_ext (nx_ext,ny_ext) REAL :: qv_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) 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=256) :: 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(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,-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 getmapr(iproj,scale,latnot,trlon,x0,y0) 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', & 'Dataset problem in GETNMCETA212.', & ' ' !STOP ISTATUS = -888 GOTO 999 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 k=1,max_nr3d ! write (6,'(/i5,4x,6(i7))') k,(var_lev3d(k,n,1),n=1,n3dvars) ! END DO DO k=1,max_nr3d DO n=2,n3dvars IF ( var_lev3d(k,1,1) /= var_lev3d(k,n,1) ) THEN WRITE (6,'(a,/,a,I3,a,I7,/,a,I7,a,I3,/,a,/)') & 'Variables were not at the same level.', & 'The K (k= ',k,')th level for variable 1 is ',var_lev3d(k,1,1), & 'but it is ',var_lev3d(k,n,1),' for variable n = ',n, & 'Dataset problem in GETNMCETA212.' !STOP ISTATUS = -888 GOTO 999 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 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 psfc_ext (i,j) = var_grb2d(i,j,1,1) ! already with Pa unit (F.KONG) 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 (soilmodel_option == 1) 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 ! soilmodel_option 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 ! F.KONG add 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 IF( var_nr2d(2,2) == 0.) THEN qv_2m_ext(i,j)= -999.0 ELSE qv_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 ! end F.KONG mod 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 istatus = 1 ! !----------------------------------------------------------------------- ! ! Reset map projection to previous values ! !----------------------------------------------------------------------- ! 999 CONTINUE CALL setmapr(iproj,scale,latnot,trlon) CALL setorig(1,x0,y0) DEALLOCATE(var_grb2d,var_grb3d,rcdata,var_lev3d) DEALLOCATE(utmp) DEALLOCATE(vtmp) RETURN END SUBROUTINE getnmceta212 ! !################################################################## !################################################################## !###### ###### !###### SUBROUTINE GETNMCETA218 ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! SUBROUTINE getnmceta218(nx_ext,ny_ext,nz_ext,nzsoil_ext, & 1,11 nstyps_ext,dir_extd,extdname,extdinit,extdfcst, & domain_tile,npr,npc, & 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, & t_2m_ext, qv_2m_ext, u_10m_ext, v_10m_ext, & istatus) ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! Reads a NMC GRIB (Grid #218, 12km) ETA file for processing by ! ext2arps, a program that converts external files to ARPS variables ! and format. This subroutine calls getnmceta218_data for each tile. ! ! NOTE: ! ! The description about Eta 12km tiled output is available at ! http://www.emc.ncep.noaa.gov/mmb/research/tiles.218.html. ! ! The data can be downloaded from ! ftpprd.ncep.noaa.gov/pub/data/nccf/com/nam/prod/nam.YYYYMMDD/tiles.tHHz ! ! This subroutine assumes that data has been downloaded to local ! disk and stored in "dir_extd"/"extdname".tCCz.awip218HH.YY ! ! where: ! ! o "dir_extd" is specified inside arps.input; ! o "extdname" supposed to be "eta" or others; ! o CC is the model cycle, extracted from "extdinit" (00,06,12,18); ! o HH is the forecast hour, extracted from "extdfcst"; ! o YY is the tile number, determined from "domain_tile". ! ! This subroutine also requires ASCII files of the grid lat/lon ! for each tile, which can be got from ! ftp://ftpprd.ncep.noaa.gov/pub/emc/mmb/mmbpll/g218tiles.latlon ! and they should be stored to local disk inside directory ! "dir_extd"/. ! !----------------------------------------------------------------------- ! ! AUTHOR: Yunheng Wang ! 07/07/2004 ! ! MODIFICATION HISTORY: ! !----------------------------------------------------------------------- ! ! INPUT: ! ! nx_ext ! ny_ext ! nz_ext ! nzsoil_ext ! nstyp_ext ! 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 ! ! domain_tile A 2D integer array which specifies the tiles to be read ! and their order. ! npr tile number per row of the above 2D array ! npc tile number per column of the above 2D array ! ! 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) ! ! T_2m_ext Temperature at 2m AGL ! qv_2m_ext Specific Humidity at 2m AGL ! U_10m_ext U at 10m AGL ! V_10m_ext V at 10m AGL ! ! 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) ! !----------------------------------------------------------------------- IMPLICIT NONE ! !----------------------------------------------------------------------- ! ! Include files ! !----------------------------------------------------------------------- ! INCLUDE 'globcst.inc' INCLUDE 'phycst.inc' INCLUDE 'grid.inc' INCLUDE 'gribcst.inc' !----------------------------------------------------------------------- ! ! Variable declarations ! !----------------------------------------------------------------------- CHARACTER (LEN=*), INTENT(INOUT) :: dir_extd CHARACTER (LEN=*), INTENT(IN) :: extdname CHARACTER (LEN=19), INTENT(IN) :: extdinit CHARACTER (LEN=9), INTENT(IN) :: extdfcst INTEGER, INTENT(IN) :: nx_ext,ny_ext,nz_ext INTEGER, INTENT(IN) :: nzsoil_ext, nstyps_ext INTEGER, INTENT(IN) :: npr,npc INTEGER, INTENT(IN) :: domain_tile(npr,npc) ! !----------------------------------------------------------------------- ! ! External grid variables ! !----------------------------------------------------------------------- ! INTEGER, INTENT(OUT) :: iproj_ext REAL, INTENT(OUT) :: scale_ext,trlon_ext REAL, INTENT(OUT) :: latnot_ext(2) REAL, INTENT(OUT) :: x0_ext,y0_ext ! !----------------------------------------------------------------------- ! ! Output external variable arrays ! !----------------------------------------------------------------------- ! REAL, INTENT(OUT) :: lat_ext(nx_ext,ny_ext) REAL, INTENT(OUT) :: lon_ext(nx_ext,ny_ext) REAL, INTENT(OUT) :: p_ext (nx_ext,ny_ext,nz_ext) ! Pressure (Pascals) REAL, INTENT(OUT) :: hgt_ext(nx_ext,ny_ext,nz_ext) ! Height (m) REAL, INTENT(OUT) :: t_ext (nx_ext,ny_ext,nz_ext) ! Temperature (K) REAL, INTENT(OUT) :: qv_ext (nx_ext,ny_ext,nz_ext) ! Specific humidity (kg/kg) REAL, INTENT(OUT) :: u_ext (nx_ext,ny_ext,nz_ext) ! Eastward wind component REAL, INTENT(OUT) :: v_ext (nx_ext,ny_ext,nz_ext) ! Northward wind component REAL, INTENT(OUT) :: qc_ext (nx_ext,ny_ext,nz_ext) ! Cloud H2O mixing ratio (kg/kg) REAL, INTENT(OUT) :: qr_ext (nx_ext,ny_ext,nz_ext) ! Rain H2O mixing ratio (kg/kg) REAL, INTENT(OUT) :: qi_ext (nx_ext,ny_ext,nz_ext) ! Ice mixing ratio (kg/kg) REAL, INTENT(OUT) :: qs_ext (nx_ext,ny_ext,nz_ext) ! Snow mixing ratio (kg/kg) REAL, INTENT(OUT) :: qh_ext (nx_ext,ny_ext,nz_ext) ! Hail mixing ratio (kg/kg) REAL, INTENT(OUT) :: zpsoil_ext(nx_ext,ny_ext,nzsoil_ext) !Soil depths (m) REAL, INTENT(OUT) :: tsoil_ext (nx_ext,ny_ext,nzsoil_ext) ! Soil temperature (K) REAL, INTENT(OUT) :: qsoil_ext (nx_ext,ny_ext,nzsoil_ext) ! Soil moisture (m**3/m**3) REAL, INTENT(OUT) :: wetcanp_ext (nx_ext,ny_ext) ! Canopy water amount REAL, INTENT(OUT) :: snowdpth_ext(nx_ext,ny_ext) ! Snow depth (m) REAL, INTENT(OUT) :: trn_ext (nx_ext,ny_ext) ! External terrain (m) REAL, INTENT(OUT) :: psfc_ext (nx_ext,ny_ext) ! Surface pressure (Pa) INTEGER, INTENT(OUT) :: soiltyp_ext (nx_ext,ny_ext) ! Soil type REAL, INTENT(OUT) :: t_2m_ext (nx_ext,ny_ext) REAL, INTENT(OUT) :: qv_2m_ext(nx_ext,ny_ext) REAL, INTENT(OUT) :: u_10m_ext(nx_ext,ny_ext) REAL, INTENT(OUT) :: v_10m_ext(nx_ext,ny_ext) ! !----------------------------------------------------------------------- ! ! Other external variable arrays ! !----------------------------------------------------------------------- ! REAL :: dx_ext,dy_ext REAL :: x_ext(nx_ext) REAL :: y_ext(ny_ext) REAL :: latsw, lonsw REAL :: latsw_t, lonsw_t INTEGER :: istatus ! !----------------------------------------------------------------------- ! ! Misc internal variables ! !----------------------------------------------------------------------- ! CHARACTER (LEN=256) :: gribfile CHARACTER (LEN=256) :: latlonfile CHARACTER (LEN=13) :: gribtime INTEGER :: it,jt INTEGER :: i,j,k,kk INTEGER :: iyr,imo,iday INTEGER :: ihr,imin,isec INTEGER :: ifhr,ifmin,ifsec INTEGER :: nxgrb,nygrb INTEGER :: offset_x, offset_y INTEGER :: inunit, dir_len REAL :: alat, alon INTEGER :: iret ! Return flag ! !----------------------------------------------------------------------- ! ! 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 REAL, ALLOCATABLE :: utmp(:,:), vtmp(:,:) INTEGER, PARAMETER :: nxmax = 69, nymax = 72 ! max size of the tiles INTEGER, PARAMETER :: nxbdy = 62, nybdy = 68 ! size of tile at boundary REAL, ALLOCATABLE :: qvs_ext(:,:,:) REAL, ALLOCATABLE :: lat(:,:), lon(:,:) ! !----------------------------------------------------------------------- ! ! Original grid variables ! !----------------------------------------------------------------------- ! INTEGER :: iproj REAL :: scale,trlon,x0,y0 REAL :: latnot(2) ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Begining of executable code ..... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ALLOCATE(var_grb2d(nxmax,nymax,n2dvs,n2dlvt), STAT = iret) ALLOCATE(var_grb3d(nxmax,nymax,nz_ext,n3dvs,n3dlvt), STAT = iret) ALLOCATE(rcdata(nxmax*nymax), STAT = iret) ALLOCATE(var_lev3d(nz_ext,n3dvs,n3dlvt), STAT = iret) ALLOCATE(utmp(nx_ext,ny_ext), STAT = iret) ALLOCATE(vtmp(nx_ext,ny_ext), STAT = iret) !----------------------------------------------------------------------- ! ! Save orignal grid map projection information ! !----------------------------------------------------------------------- CALL getmapr(iproj,scale,latnot,trlon,x0,y0) !----------------------------------------------------------------------- ! ! Construct gribfile name and gribtime ! !----------------------------------------------------------------------- ifhr=0 ifmin=0 ifsec=0 ! 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 READ(extdfcst,'(i3)') ifhr WRITE (gribtime,'(i4.4,i2.2,i2.2,i2.2,a1,i2.2)') & iyr,imo,iday,ihr,'f',ifhr dir_len=LEN_TRIM(dir_extd) IF( dir_len == 0 .OR. dir_extd(1:dir_len) == ' ' ) THEN dir_extd = '.' dir_len = 1 END IF IF( dir_extd(dir_len:dir_len) /= '/' ) THEN dir_len = dir_len + 1 dir_extd(dir_len:dir_len) = '/' END IF !----------------------------------------------------------------------- ! ! Begining loop over each tile ! !----------------------------------------------------------------------- DO jt = 1,npc DO it = 1,npr WRITE(gribfile,'(a,3(a,I2.2))') & TRIM(extdname),'.t',ihr,'z.awip218',ifhr,'.',domain_tile(it,jt) nxgrb = nxmax nygrb = nymax IF ( MOD(domain_tile(it,jt), 9) == 0) nxgrb = nxbdy ! East boudnary IF ( domain_tile(it,jt) >= 46 ) nygrb = nybdy ! North boudnary !----------------------------------------------------------------------- ! ! Call getnmceta218_data ! !----------------------------------------------------------------------- CALL getnmceta218_data(nx_ext,ny_ext,nz_ext,nzsoil_ext, & nstyps_ext,nxgrb,nygrb,it,jt, & dir_extd,gribfile,gribtime,dx_ext,dy_ext, & iproj_ext,scale_ext, & trlon_ext,latnot_ext,latsw_t,lonsw_t, & 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, & t_2m_ext, qv_2m_ext, u_10m_ext, v_10m_ext, & var_grb2d,var_grb3d,rcdata,var_lev3d,istatus) IF (it == 1 .AND. jt == 1) THEN latsw = latsw_t lonsw = lonsw_t END IF IF (istatus < 0) GOTO 999 END DO END DO !----------------------------------------------------------------------- ! ! Check for data consistence. ! !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! ! Set map projection and check for lat_ext, lon_ext etc. ! !----------------------------------------------------------------------- ALLOCATE(lat(nx_ext,ny_ext), STAT = iret) ALLOCATE(lon(nx_ext,ny_ext), STAT = iret) 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,lon) !----------------------------------------------------------------------- ! ! Get grid lat/lon from external ASCII file ! ! NOTE: This block is not necessary after we make sure everything ! is working correctly above. ! !----------------------------------------------------------------------- ! lat_ext(:,:) = lat ! lon_ext(:,:) = lon lat_ext(:,:) = 0.0 lon_ext(:,:) = 0.0 inunit = 123 DO jt = 1, npc DO it = 1, npr offset_x = (it-1)*nxmax offset_y = (jt-1)*nymax WRITE(latlonfile,'(a,a,I2.2)') dir_extd(1:dir_len), & 'latlon.grid218.',domain_tile(it,jt) OPEN(inunit, FILE=TRIM(latlonfile),FORM='FORMATTED',STATUS='OLD') DO WHILE(.TRUE.) READ(inunit,*, END=995) i,j, alat, alon lat_ext(offset_x+i,offset_y+j) = alat lon_ext(offset_x+i,offset_y+j) = -1*alon END DO 995 CONTINUE CLOSE(inunit) END DO END DO IF( ANY(lat_ext <= 12.0) .OR. ANY(lat_ext >= 62.0) .OR. & ANY(lon_ext >= -49.0) .OR. ANY(lon_ext <= -160.0) ) THEN WRITE(6,*) 'Some lat_ext or lon_ext have wrong value. ', & 'Please check the code.' STOP END IF !----------------------------------------------------------------------- ! ! Check grid latitude and longitude ! !----------------------------------------------------------------------- DO j = 1, ny_ext DO i = 1, nx_ext IF ( ABS(lat(i,j)-lat_ext(i,j)) > 0.01 .OR. & ABS(lon(i,j)-lon_ext(i,j)) > 0.01 ) THEN WRITE(6,'(2(a,I3),/,2(a,2F7.2,/))') & 'Unmatched grid at i = ',i,' j = ',j, & 'Expecting lat,lon = ',lat_ext(i,j), lon_ext(i,j), & 'Found lat,lon = ',lat(i,j),lon(i,j) STOP END IF END DO END DO !----------------------------------------------------------------------- ! ! Convert relative humidity to specific humidity ! !----------------------------------------------------------------------- ALLOCATE(qvs_ext(nx_ext,ny_ext,nz_ext), STAT = iret) CALL check_alloc_status(iret, "getnmceta218:qvs_ext") CALL getqvs(nx_ext,ny_ext,nz_ext, 1,nx_ext,1,ny_ext,1,nz_ext, & p_ext, t_ext, qvs_ext ) DO k=1,nz_ext DO j=1,ny_ext DO i=1,nx_ext qv_ext(i,j,k) = 0.01*qv_ext(i,j,k)*qvs_ext(i,j,k) END DO END DO END DO ! ! Assume 2-m pressure is the same as psfc ! CALL getqvs(nx_ext,ny_ext,1, 1,nx_ext,1,ny_ext,1,1, & psfc_ext, t_2m_ext, qvs_ext ) DO j=1,ny_ext DO i=1,nx_ext qv_2m_ext(i,j) = 0.01*qv_2m_ext(i,j)*qvs_ext(i,j,1) END DO END DO !----------------------------------------------------------------------- ! ! Rotate winds to be relative to true north. ! The Eta data are sent as grid-relative. ! !----------------------------------------------------------------------- ! DO k=1,nz_ext 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 istatus = 1 ! !----------------------------------------------------------------------- ! ! Reset map projection to previous values ! !----------------------------------------------------------------------- ! 999 CONTINUE CALL setmapr(iproj,scale,latnot,trlon) CALL setorig(1,x0,y0) DEALLOCATE(var_grb2d,var_grb3d,rcdata,var_lev3d) DEALLOCATE(utmp, vtmp) DEALLOCATE(qvs_ext) DEALLOCATE(lat,lon) RETURN END SUBROUTINE getnmceta218 SUBROUTINE getnmceta218_data(nx_ext,ny_ext,nz_ext,nzsoil_ext, & 1,1 nstyps_ext,nxgrb,nygrb,tile_x,tile_y, & dir_extd,extdname,gribtime,dx_ext,dy_ext, & iproj_ext,scale_ext, & trlon_ext,latnot_ext,latsw_t,lonsw_t, & zpsoil_ext,p_ext,hgt_ext,t_ext,rh_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, & t_2m_ext, rh_2m_ext, u_10m_ext, v_10m_ext, & var_grb2d,var_grb3d,rcdata,var_lev3d,istatus) ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! Reads a NMC GRIB (Grid #218, 12km) ETA file for processing by ! EXT2ARPS. This subroutine reads one tile data and insert the data ! into the global domain (maybe multiple tiles or one tile) arrays. ! !----------------------------------------------------------------------- ! ! AUTHOR: Yunheng Wang ! 07/06/2004 ! Based on subroutine getnmceta212. ! ! MODIFICATION HISTORY: ! !----------------------------------------------------------------------- ! ! INPUT: ! ! nx_ext X dimension size of the global domain (maybe mupltiple tiles) ! ny_ext Y dimension size of the global domain (maybe mupltiple tiles) ! nz_ext Z dimension size (39 pressure levels) ! nzsoil_ext 4 layers (0-10cm, 10-40cm, 40-100cm, 100-200cm) ! nstyps_ext =1 ! nxgrb Tile size in X direction (69/62) ! nygrb Tile size in Y direction (72/68) ! tile_x This tile location in the multiple tile domain in X direction ! tile_y This tile location in the multiple tile domain in Y direction ! dir_extd Directory name for external file ! extdname External file name ! gribtime Data valid time, YYYYMMDDHHfHH ! ! ! 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) ! latsw_t latitude of the southwest corner of this tile (degrees N) ! lonsw_t longitude of the southwest corner of this tile (degrees E) ! p_ext pressure (Pascal) ! hgt_ext height (m) ! t_ext temperature (K) ! rh_ext relative humidity (%) ! 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) ! ! T_2m_ext Temperature at 2m AGL ! Rh_2m_ext Relative Humidity at 2m AGL ! U_10m_ext U at 10m AGL ! V_10m_ext V at 10m AGL ! ! 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,1) - Surface pressure (Pa) ! var_grb2d(nxgrb,nygrb,2,1) - Geopotential height (gpm) ! var_grb2d(nxgrb,nygrb,3,1) - Surface temperature (K) ! var_grb2d(nxgrb,nygrb,4,1) - Plant canopy surface ! water (kg/m**2) ! var_grb2d(nxgrb,nygrb,1,2) - 2-m temperature (Pa) ! var_grb2d(nxgrb,nygrb,2,2) - 2-m specific humidity (kg/kg) ! var_grb2d(nxgrb,nygrb,3,2) - 10-m U velocity (m/s) ! var_grb2d(nxgrb,nygrb,4,2) - 10-m V velocity (m/s) ! ! rcdata ! ! var_lev3d ! !----------------------------------------------------------------------- ! ! Variable Declarations. ! !----------------------------------------------------------------------- ! IMPLICIT NONE INCLUDE 'gribcst.inc' CHARACTER(LEN=*), INTENT(INOUT) :: dir_extd CHARACTER(LEN=*), INTENT(IN) :: extdname CHARACTER(LEN=*), INTENT(IN) :: gribtime !----------------------------------------------------------------------- ! ! External grid variables ! !----------------------------------------------------------------------- ! INTEGER, INTENT(IN) :: nx_ext,ny_ext ! multiple tiles size INTEGER, INTENT(IN) :: nxgrb, nygrb ! one tile size INTEGER, INTENT(IN) :: nz_ext INTEGER, INTENT(IN) :: nzsoil_ext, nstyps_ext INTEGER, INTENT(IN) :: tile_x, tile_y INTEGER, INTENT(OUT) :: iproj_ext REAL, INTENT(OUT) :: scale_ext,trlon_ext REAL, INTENT(OUT) :: latnot_ext(2) REAL, INTENT(OUT) :: latsw_t REAL, INTENT(OUT) :: lonsw_t REAL, INTENT(OUT) :: dx_ext,dy_ext ! !----------------------------------------------------------------------- ! ! Output external variable arrays ! !----------------------------------------------------------------------- ! REAL, INTENT(OUT) :: p_ext (nx_ext,ny_ext,nz_ext) ! Pressure (Pascals) REAL, INTENT(OUT) :: hgt_ext(nx_ext,ny_ext,nz_ext) ! Height (m) REAL, INTENT(OUT) :: t_ext (nx_ext,ny_ext,nz_ext) ! Temperature (K) REAL, INTENT(OUT) :: rh_ext (nx_ext,ny_ext,nz_ext) ! Relative humidity (%) REAL, INTENT(OUT) :: u_ext (nx_ext,ny_ext,nz_ext) ! Eastward wind component REAL, INTENT(OUT) :: v_ext (nx_ext,ny_ext,nz_ext) ! Northward wind component REAL, INTENT(OUT) :: qc_ext (nx_ext,ny_ext,nz_ext) ! Cloud H2O mixing ratio (kg/kg) REAL, INTENT(OUT) :: qr_ext (nx_ext,ny_ext,nz_ext) ! Rain H2O mixing ratio (kg/kg) REAL, INTENT(OUT) :: qi_ext (nx_ext,ny_ext,nz_ext) ! Ice mixing ratio (kg/kg) REAL, INTENT(OUT) :: qs_ext (nx_ext,ny_ext,nz_ext) ! Snow mixing ratio (kg/kg) REAL, INTENT(OUT) :: qh_ext (nx_ext,ny_ext,nz_ext) ! Hail mixing ratio (kg/kg) REAL, INTENT(OUT) :: zpsoil_ext(nx_ext,ny_ext,nzsoil_ext) !Soil depths (m) REAL, INTENT(OUT) :: tsoil_ext (nx_ext,ny_ext,nzsoil_ext) ! Soil temperature (K) REAL, INTENT(OUT) :: qsoil_ext (nx_ext,ny_ext,nzsoil_ext) ! Soil moisture (m**3/m**3) REAL, INTENT(OUT) :: wetcanp_ext (nx_ext,ny_ext) ! Canopy water amount REAL, INTENT(OUT) :: snowdpth_ext(nx_ext,ny_ext) ! Snow depth (m) REAL, INTENT(OUT) :: trn_ext (nx_ext,ny_ext) ! External terrain (m) REAL, INTENT(OUT) :: psfc_ext (nx_ext,ny_ext) ! Surface pressure (Pa) INTEGER, INTENT(OUT) :: soiltyp_ext (nx_ext,ny_ext) ! Soil type REAL, INTENT(OUT) :: t_2m_ext (nx_ext,ny_ext) REAL, INTENT(OUT) :: rh_2m_ext (nx_ext,ny_ext) REAL, INTENT(OUT) :: u_10m_ext (nx_ext,ny_ext) REAL, INTENT(OUT) :: v_10m_ext (nx_ext,ny_ext) ! !----------------------------------------------------------------------- ! ! Include files ! !----------------------------------------------------------------------- ! INCLUDE 'globcst.inc' INCLUDE 'phycst.inc' INCLUDE 'grid.inc' ! !----------------------------------------------------------------------- ! ! Other external variable arrays ! !----------------------------------------------------------------------- ! INTEGER, INTENT(OUT) :: istatus ! !----------------------------------------------------------------------- ! ! Work arrays for storing grib data ! !----------------------------------------------------------------------- ! REAL, INTENT(OUT) :: var_grb2d(nxgrb,nygrb,n2dvs,n2dlvt) ! GRIB variables REAL, INTENT(OUT) :: var_grb3d(nxgrb,nygrb,nz_ext,n3dvs,n3dlvt) ! GRIB 3-D variables INTEGER, INTENT(OUT) :: var_lev3d(nz_ext,n3dvs,n3dlvt) ! Levels (hybrid) for ! each 3-D variable REAL, INTENT(OUT) :: rcdata(nxgrb*nygrb) ! temporary data array ! !----------------------------------------------------------------------- ! ! Misc internal variables ! !----------------------------------------------------------------------- ! CHARACTER (LEN=256) :: gribfile INTEGER :: offset_x, offset_y INTEGER :: i,j,k,kk INTEGER :: grbflen INTEGER :: m,n,nz1,max_nr2d,max_nr3d REAL :: soildepth_ext(5) INTEGER :: chklev, lvscan INTEGER :: iret ! Return flag ! !----------------------------------------------------------------------- ! ! 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 LOGICAL, PARAMETER :: verbose = .FALSE. INTEGER, PARAMETER :: nxmax = 69, nymax = 72 ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! DATA soildepth_ext/0.0,-0.05,-0.25,-0.70,-1.50/ ! ground and soil layers gribfile = TRIM(dir_extd)//TRIM(extdname) grbflen = LEN_TRIM(gribfile) ! !----------------------------------------------------------------------- ! ! RDNMCGRB reads NMC GRIB data ! !----------------------------------------------------------------------- ! gridtyp = eta218grid mproj_grb = eta218proj n2dvars = eta218nvs2d n2dlvtps = eta218nlvt2d DO k=1,n2dlvtps DO n=1,n2dvars var_id2d(n,k) = eta218var_id2d(n,k) END DO levtyp2d(k) = eta218levs2d(k) END DO n3dvars = eta218nvs3d n3dlvtps = eta218nlvt3d DO m=1,n3dlvtps DO n=1,n3dvars var_id3d(n,m) = eta218var_id3d(n,m) END DO levtyp3d(m) = eta218levs3d(m) END DO CALL rdnmcgrb(nxgrb,nygrb,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 ',TRIM(gribfile), & ' Dataset problem in GETNMCETA218.', & ' ' istatus = -888 RETURN END IF IF ( max_nr2d == 0 ) THEN WRITE (6,'(a)') & 'No 2-D variables was found in the GRIB file' END IF IF (verbose) THEN WRITE (6,'(/a7,2x,6(i7))') 'Lev\\VID',(var_id3d(n,1),n=1,n3dvars) DO k=1,max_nr3d WRITE (6,'(/i5,4x,6(i7))') k,(var_lev3d(k,n,1),n=1,n3dvars) END DO END IF 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.', & 'Dataset problem in GETNMCETA218.', & ' ' istatus = -888 RETURN 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 ! ! Do some check here, such as nxgrd, nygrd. ! IF (iproj_grb /= mproj_grb) THEN WRITE(6,*) 'Map projection unmatch, expecting ',mproj_grb,' Found ',iproj_ext STOP END IF IF(ni_grb /= nxgrb .OR. nj_grb /= nygrb) THEN WRITE(6,'(2(a,I4),/,2(a,I4))') & 'Data array size inconsistent, expecting nxgrb = ',nxgrb, & ' nygrb = ',nygrb, & ' Found: ni_grb = ',ni_grb, & ' nj_grb = ',nj_grb istatus = -888 RETURN END IF dx_ext = di_grb dy_ext = dj_grb ! 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) latsw_t = latsw lonsw_t = lonsw offset_x = (tile_x-1)*nxmax ! Offset of this tile in the multiple offset_y = (tile_y-1)*nymax ! tile domain. !--------------------------------------------------------------------- ! Define soil depths !--------------------------------------------------------------------- DO k=1,nzsoil_ext DO j=1,nygrb DO i=1,nxgrb zpsoil_ext(offset_x+i,offset_y+j,k) = soildepth_ext(k) END DO END DO END DO ! !----------------------------------------------------------------------- ! ! Retrieve 2-D variables ! !----------------------------------------------------------------------- ! DO j=1,nygrb DO i=1,nxgrb IF ( var_nr2d(1,1) == 0 ) THEN psfc_ext (offset_x+i,offset_y+j) = -999.0 ELSE psfc_ext (offset_x+i,offset_y+j) = var_grb2d(i,j,1,1) END IF IF ( var_nr2d(2,1) == 0 ) THEN trn_ext (offset_x+i,offset_y+j) = -999.0 ELSE trn_ext (offset_x+i,offset_y+j) = var_grb2d(i,j,2,1) END IF IF ( var_nr2d(4,1) == 0 ) THEN wetcanp_ext(offset_x+i,offset_y+j) = -999.0 ELSE wetcanp_ext(offset_x+i,offset_y+j) = var_grb2d(i,j,4,1)*1.e-3 ! in meter END IF IF ( var_nr2d(6,1) == 0 ) THEN snowdpth_ext(offset_x+i,offset_y+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(offset_x+i,offset_y+j) = 0.01 * var_grb2d(i,j,6,1) ! in meters END IF IF( var_nr2d(1,2) == 0.) THEN t_2m_ext(offset_x+i,offset_y+j)= -999.0 ELSE t_2m_ext(offset_x+i,offset_y+j)= var_grb2d(i,j,1,2) ! 2-m temperature END IF IF( var_nr2d(2,2) == 0.) THEN rh_2m_ext(offset_x+i,offset_y+j)= -999.0 ELSE rh_2m_ext(offset_x+i,offset_y+j)= var_grb2d(i,j,2,2) ! 2-m relative humidity END IF IF( var_nr2d(3,2) == 0.) THEN u_10m_ext(offset_x+i,offset_y+j)= -999.0 ELSE u_10m_ext(offset_x+i,offset_y+j)= var_grb2d(i,j,3,2) ! 10-m U END IF IF( var_nr2d(4,2) == 0.) THEN v_10m_ext(offset_x+i,offset_y+j)= -999.0 ELSE v_10m_ext(offset_x+i,offset_y+j)= var_grb2d(i,j,4,2) ! 10-m V 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,nygrb DO i=1,nxgrb p_ext (offset_x+i,offset_y+j,kk) = 100.0 * FLOAT(var_lev3d(k,1,1)) ! Pressure t_ext (offset_x+i,offset_y+j,kk) = var_grb3d(i,j,k,1,1) ! Temperature (K) rh_ext (offset_x+i,offset_y+j,kk) = var_grb3d(i,j,k,2,1) ! Relative humidity (%) u_ext (offset_x+i,offset_y+j,kk) = var_grb3d(i,j,k,3,1) ! u wind (m/s) v_ext (offset_x+i,offset_y+j,kk) = var_grb3d(i,j,k,4,1) ! v wind (m/s) hgt_ext(offset_x+i,offset_y+j,kk) = var_grb3d(i,j,k,5,1) ! height (m) qc_ext (offset_x+i,offset_y+j,kk) = -999. qr_ext (offset_x+i,offset_y+j,kk) = -999. qi_ext (offset_x+i,offset_y+j,kk) = -999. qs_ext (offset_x+i,offset_y+j,kk) = -999. qh_ext (offset_x+i,offset_y+j,kk) = -999. END DO END DO END DO ! !----------------------------------------------------------------------- ! ! Retrive land surface variables ! !----------------------------------------------------------------------- DO j = 1,nygrb DO i = 1, nxgrb IF ( var_nr3d(1,2) == 0 ) THEN DO k=1,nzsoil_ext tsoil_ext (offset_x+i,offset_y+j,k) = -999.0 qsoil_ext (offset_x+i,offset_y+j,k) = -999.0 END DO ELSE IF (soilmodel_option == 1) THEN ! Old ARPS Force-Restore Soil Model tsoil_ext(offset_x+i,offset_y+j,1) = var_grb2d(i,j,3,1) IF ( nint(var_grb2d(i,j,5,1)) == 1 ) THEN ! Over land tsoil_ext(offset_x+i,offset_y+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(offset_x+i,offset_y+j) = 0 ELSE ! Over sea tsoil_ext(offset_x+i,offset_y+j,2) = var_grb2d(i,j,3,1) soiltyp_ext(offset_x+i,offset_y+j) = 13 ! Set soil type to water END IF qsoil_ext(offset_x+i,offset_y+j,1) = var_grb3d(i,j,1,2,2) qsoil_ext(offset_x+i,offset_y+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 IF ( nint(var_grb2d(i,j,5,1)) == 1 ) THEN ! Over land tsoil_ext (offset_x+i,offset_y+j,1) = var_grb2d(i,j,3,1) ! Ground temperature qsoil_ext (offset_x+i,offset_y+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 (offset_x+i,offset_y+j,k) = var_grb3d(i,j,k-1,1,2) ! Not a bug; qsoil_ext (offset_x+i,offset_y+j,k) = var_grb3d(i,j,k-1,2,2) END DO ELSE ! Over water DO k=1,nzsoil_ext tsoil_ext (offset_x+i,offset_y+j,k) = var_grb2d(i,j,3,1) ! Ground temperature qsoil_ext (offset_x+i,offset_y+j,k) = 1. ! 100% water END DO END IF END IF ! soilmodel_option END IF END DO END DO istatus = 1 RETURN END SUBROUTINE getnmceta218_data ! !################################################################## !################################################################## !###### ###### !###### 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 or Grid 252) 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=256) :: 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 ! !----------------------------------------------------------------------- ! IF (extdopt == 18) THEN gridtyp = rucn252grid ELSE gridtyp = rucn236grid END IF 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=256) :: 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 k=1,max_nr3d ! WRITE (6,'(i5,4x,6(i7))') k,(var_lev3d(k,n,1),n=1,n3dvars) ! END DO 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,nzsoil_ext, & 2,6 dir_extd,extdname,extdopt,nofixdim,lon_0_360, & 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, & t_2m_ext, rh_2m_ext, u_10m_ext, v_10m_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. ! ! 2004/01/15 (F. KONG) ! Modified to read in from grid 3 four more near surface variables: ! 2-m temperature and humidity, and 10-m wind (u, v) ! ! 08/11/2005 (Y. Wang) ! Modified to handle cases with Greenwich Meridian in the domain. Note ! that a new parameter, lon_0_360 was passed in and an integer array ! ii(nx_ext) was decalared. ! ! Modified to work with data over subregion instead of global data. A ! new parameter nofixdim was passed in. ! !----------------------------------------------------------------------- ! ! 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 ! qsoil_ext Soil moisture ! 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 ! ! 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, INTENT(IN) :: extdopt INTEGER, INTENT(IN) :: nofixdim LOGICAL, INTENT(IN) :: lon_0_360 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,nzsoil_ext REAL :: lat_ext(nx_ext,ny_ext) REAL :: lon_ext(nx_ext,ny_ext) REAL :: zpsoil_ext(nx_ext,ny_ext,nzsoil_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 :: tsoil_ext (nx_ext,ny_ext,nzsoil_ext) ! Temperature at surface (K) REAL :: qsoil_ext (nx_ext,ny_ext,nzsoil_ext) ! Deep soil temperature (K) 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 :: 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) INTEGER :: soiltyp_ext(nx_ext,ny_ext) ! !----------------------------------------------------------------------- ! ! Other external variable arrays ! !----------------------------------------------------------------------- ! ! REAL :: x_ext(nx_ext) ! REAL :: y_ext(ny_ext) INTEGER :: ii(nx_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=256) :: gribfile CHARACTER (LEN=14) :: 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 REAL, PARAMETER :: soildepth_ext(3) = (/0.0, 0.05, 1.05/) ! Surface and two soil layers 0-10cm & 10-200cm REAL :: rtmp ! !----------------------------------------------------------------------- ! ! 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... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! x0_ext = 0.0 y0_ext = 0.0 !--------------------------------------------------------------------- ! ! Define soil depths ! !--------------------------------------------------------------------- IF(soilmodel_option == 1 .AND. nzsoil_ext /= 2) THEN WRITE(6,'(2a,I3,a/2a,I2,a/,a,a)') & 'NCEP GFS grid #3 only provides 2 soil layers, However,',& ' You are trying to extract ',nzsoil_ext, ' layers.', & ' for ARPS two-layer force-restore model ', & '(soilmodel_option = ',soilmodel_option,')', & ' Please check the code ext2arps.f90 for NCEPAVN3 grid.',& ' Terminating ...' STOP ELSE IF( soilmodel_option == 2 .AND. nzsoil_ext /= 3) THEN WRITE(6,'(2a,I3,a/2a,I2,a/,a,a)') & 'NCEP GFS grid #3 only provides 2 soil layers, However,',& ' You are trying to extract ',nzsoil_ext, ' layers.', & ' for ARPS multi-layer OUSoil model ', & '(soilmodel_option = ',soilmodel_option,')', & ' Please check the code ext2arps.f90 for NCEPAVN3 grid.',& ' Terminating ...' STOP END IF IF(soilmodel_option == 1) THEN DO k=1,nzsoil_ext DO j=1,ny_ext DO i=1,nx_ext zpsoil_ext(i,j,k) = soildepth_ext(k+1) END DO END DO END DO ELSE 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 END IF 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)) var_grb2d = 0.0 var_grb3d = 0.0 rcdata = 0.0 var_lev3d = 0.0 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 ! ! AVN/GFS exceeds 99 hours, so we must use a 3 digit naming scheme. ! !WRITE (gribtime,'(i4.4,i2.2,i2.2,i2.2,a1,i3.3)') & ! WDT 2003-12-23 Use 2 digits, not 3, for forecast hour RLC 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:14) grbflen = grbflen + lenrun + 15 ! !----------------------------------------------------------------------- ! ! RDNMCGRB reads NMC GRIB data ! !----------------------------------------------------------------------- ! gridtyp = avn3grid mproj_grb = avn3proj IF (extdopt == 51 .OR. nofixdim == 1) THEN WRITE(6,*) 'Special handling for extdopt ', extdopt gridtyp = 255 END IF 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 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) ! END DO 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 j=1, ny_ext DO i=1, nx_ext lon_ext(i,j)= lonsw + (i-1) * dx_ext lat_ext(i,j)= latsw + (j-1) * dy_ext END DO END DO IF ( lon_0_360 ) THEN DO i = 1,nx_ext ii(i) = i END DO ELSE IF (mod(nx_ext,2) /= 0) THEN WRITE(6,'(a/)') 'Wrong size of nx_ext in GETNCEPAVN3 for lon_0_360.' CALL arpsstop('Wrong nx_ext size.',1) END IF DO i = 1,nx_ext/2 ! map 1-180 to 181-360 ii(i) = i + nx_ext/2 END DO DO i = nx_ext/2+1, nx_ext ! map 181-360 to 1-180 ii(i) = i - nx_ext/2 END DO WHERE (lon_ext >= 180) lon_ext = lon_ext - 360 DO j = 1,ny_ext ! swap lat_ext & lon_ext DO i = 1,nx_ext/2 rtmp = lat_ext(ii(i),j) lat_ext(ii(i),j) = lat_ext(i,j) lat_ext(i, j) = rtmp rtmp = lon_ext(ii(i),j) lon_ext(ii(i),j) = lon_ext(i,j) lon_ext(i, j) = rtmp END DO END DO END IF 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 psfc_ext (i,j) = var_grb2d(ii(i),j,1,1) ! already Pa (F.KONG) 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 trn_ext (i,j) = var_grb2d(ii(i),j,2,1) ! why divided by g? (F.KONG) END IF IF ( var_nr3d(1,2) == 0 ) THEN tsoil_ext (i,j,:) = -999.0 qsoil_ext (i,j,:) = -999.0 ELSE IF(soilmodel_option == 1) THEN ! ARPS: two-layer force-restore model) tsoil_ext(i,j,1) = var_grb2d(ii(i),j,3,1) ! sfc temp. IF ( nint(var_grb2d(ii(i),j,5,1)) == 1 ) THEN ! soil temp over land tsoil_ext(i,j,2) = var_grb3d(ii(i),j,1,1,2) IF ( tsoil_ext (i,j,2) <= 200. ) THEN tsoil_ext(i,j,2) = tsoil_ext(i,j,1) END IF qsoil_ext(i,j,1) = var_grb3d(ii(i),j,2,2,2) qsoil_ext(i,j,2) = var_grb3d(ii(i),j,1,2,2) soiltyp_ext(i,j)= 0 ! We do not know the soil type ELSE ! sfc temp over sea tsoil_ext(i,j,2) = tsoil_ext(i,j,1) qsoil_ext(i,j,:) = 1.0 soiltyp_ext(i,j) = 13 ! Water END IF ELSE ! ARPS: multi-layer OUSoil model tsoil_ext(i,j,1) = var_grb2d(ii(i),j,3,1) ! sfc temp. IF ( nint(var_grb2d(ii(i),j,5,1)) == 1 ) THEN ! soil temp over land tsoil_ext(i,j,2) = var_grb3d(ii(i),j,2,1,2) ! 0 - 10 cm tsoil_ext(i,j,3) = var_grb3d(ii(i),j,1,1,2) ! 10 - 200 cm qsoil_ext(i,j,1) = var_grb3d(ii(i),j,2,2,2) ! 0 - 10 cm ! Value at surface assume equals ! the first layer below ground qsoil_ext(i,j,2) = var_grb3d(ii(i),j,2,2,2) ! 0 - 10 cm qsoil_ext(i,j,3) = var_grb3d(ii(i),j,1,2,2) ! 10 - 200 cm soiltyp_ext(i,j)= 0 ! We do not know the soil type ELSE ! sfc temp over sea tsoil_ext(i,j,2:3) = tsoil_ext(i,j,1) ! over water, All soil layers = SST qsoil_ext(i,j,1:3) = 1.0 ! All soil layers = 1.0 soiltyp_ext(i,j) = 13 ! Water END IF END IF ! soilmodel_option END IF IF ( var_nr2d(4,1) == 0 ) THEN wetcanp_ext(i,j) = -999.0 ELSE wetcanp_ext(i,j) = var_grb2d(ii(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(ii(i),j,6,1) ! in meters END IF ! F.KONG add IF( var_nr2d(1,4) == 0.) THEN t_2m_ext(i,j)= -999.0 ELSE t_2m_ext(i,j)= var_grb2d(ii(i),j,1,4) END IF IF( var_nr2d(2,4) == 0.) THEN rh_2m_ext(i,j)= -999.0 ELSE rh_2m_ext(i,j)= var_grb2d(ii(i),j,2,4) END IF IF( var_nr2d(3,4) == 0.) THEN u_10m_ext(i,j)= -999.0 ELSE u_10m_ext(i,j)= var_grb2d(ii(i),j,3,4) END IF IF( var_nr2d(4,4) == 0.) THEN v_10m_ext(i,j)= -999.0 ELSE v_10m_ext(i,j)= var_grb2d(ii(i),j,4,4) END IF ! end F.KONG mod 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(ii(i),j,k,1,1) u_ext (i,j,kk) = var_grb3d(ii(i),j,k,2,1) ! u wind (m/s) v_ext (i,j,kk) = var_grb3d(ii(i),j,k,3,1) ! v wind (m/s) t_ext (i,j,kk) = var_grb3d(ii(i),j,k,4,1) ! Temperature (K) qc_ext (i,j,kk) = var_grb3d(ii(i),j,k,7,1) 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(ii(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 GETNCEPAVN2 ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## SUBROUTINE getncepavn2(nx_ext,ny_ext,nz_ext,nzsoil_ext, & 2,6 dir_extd,extdname,extdopt,extdfmt,lon_0_360, & 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 NCEP GFS or NCAR/NCEP Reanalysis GRIB (Grid #2, 2.5 x 2.5 degree) ! data file for processing by ext2arps. ! !----------------------------------------------------------------------- ! ! AUTHOR: Richard Carpenter ! 2003-07-03 ! ! MODIFICATION HISTORY: ! ! 08/11/2005 (Y. Wang) ! Modified to handle cases with Greenwich Meridian in the domain. Note ! that a new parameter, lon_0_360 was passed in and an integer array ! ii(nx_ext) was decalared. Still not tested. ! !----------------------------------------------------------------------- ! ! 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 ! qsoil_ext 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 LOGICAL, INTENT(IN) :: lon_0_360 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,nzsoil_ext REAL :: lat_ext(nx_ext,ny_ext) REAL :: lon_ext(nx_ext,ny_ext) REAL :: zpsoil_ext(nx_ext,ny_ext,nzsoil_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 :: tsoil_ext (nx_ext,ny_ext,nzsoil_ext) ! Temperature at surface (K) REAL :: qsoil_ext (nx_ext,ny_ext,nzsoil_ext) ! Deep soil temperature (K) 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) ! !----------------------------------------------------------------------- ! ! Other external variable arrays ! !----------------------------------------------------------------------- ! ! REAL :: x_ext(nx_ext) ! REAL :: y_ext(ny_ext) INTEGER :: ii(nx_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=256) :: gribfile CHARACTER (LEN=14) :: 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 REAL :: rtmp REAL, PARAMETER :: soildepth_ext(3) = (/0.0, 0.05, 1.05/) ! Surface and two soil layers 0-10cm & 10-200cm ! !----------------------------------------------------------------------- ! ! 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... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! x0_ext = 0.0 y0_ext = 0.0 !--------------------------------------------------------------------- ! ! Define soil depths ! !--------------------------------------------------------------------- IF(soilmodel_option == 1 .AND. nzsoil_ext /= 2) THEN WRITE(6,'(2a,I3,a/2a,I2,a/,a,a)') & 'NCEP GFS grid #2 only provides 2 soil layers, However,',& ' You are trying to extract ',nzsoil_ext, ' layers.', & ' for ARPS two-layer force-restore model ', & '(soilmodel_option = ',soilmodel_option,')', & ' Please check the code ext2arps.f90 for NCEPAVN2 grid.',& ' Terminating ...' STOP ELSE IF( soilmodel_option == 2 .AND. nzsoil_ext /= 3) THEN WRITE(6,'(2a,I3,a/2a,I2,a/,a,a)') & 'NCEP GFS grid #2 only provides 2 soil layers, However,',& ' You are trying to extract ',nzsoil_ext, ' layers.', & ' for ARPS multi-layer OUSoil model ', & '(soilmodel_option = ',soilmodel_option,')', & ' Please check the code ext2arps.f90 for NCEPAVN2 grid.',& ' Terminating ...' STOP END IF IF(soilmodel_option == 1) THEN DO k=1,nzsoil_ext DO j=1,ny_ext DO i=1,nx_ext zpsoil_ext(i,j,k) = soildepth_ext(k+1) END DO END DO END DO ELSE 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 END IF 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)) var_grb2d = 0.0 var_grb3d = 0.0 rcdata = 0.0 var_lev3d = 0.0 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 ! ! AVN/GFS exceeds 99 hours, so we must use a 3 digit naming scheme. ! !WRITE (gribtime,'(i4.4,i2.2,i2.2,i2.2,a1,i3.3)') & ! WDT 2003-12-23 Use 2 digits, not 3, for forecast hour RLC 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:14) grbflen = grbflen + lenrun + 15 ! !----------------------------------------------------------------------- ! ! RDNMCGRB reads NMC GRIB data ! !----------------------------------------------------------------------- ! gridtyp = avn2grid mproj_grb = avn2proj IF (extdopt == 51) THEN WRITE(6,*) 'Special handling for extdopt ', extdopt gridtyp = 255 END IF n2dvars = avn2nvs2d n2dlvtps = avn2nlvt2d DO k=1,n2dlvtps DO n=1,n2dvars var_id2d(n,k) = avn2var_id2d(n,k) END DO levtyp2d(k) = avn2levs2d(k) END DO n3dvars = avn2nvs3d n3dlvtps = avn2nlvt3d DO m=1,n3dlvtps DO n=1,n3dvars var_id3d(n,m) = avn2var_id3d(n,m) END DO levtyp3d(m) = avn2levs3d(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 GETNCEPAVN2.' 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 GETNCEPAVN2.' 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 IF ( lon_0_360 ) THEN DO i = 1,nx_ext ii(i) = i END DO ELSE IF (mod(nx_ext,2) /= 0) THEN WRITE(6,'(a/)') 'Wrong size of nx_ext in GETNCEPAVN2 for lon_0_360.' CALL arpsstop('Wrong nx_ext size.',1) END IF DO i = 1,nx_ext/2 ! map 1-180 to 181-360 ii(i) = i + nx_ext/2 END DO DO i = nx_ext/2+1, nx_ext ! map 181-360 to 1-180 ii(i) = i - nx_ext/2 END DO WHERE (lon_ext >= 180.) lon_ext = lon_ext - 360 DO j = 1,ny_ext ! swap lat_ext & lon_ext DO i = 1,nx_ext/2 rtmp = lat_ext(ii(i),j) lat_ext(ii(i),j) = lat_ext(i,j) lat_ext(i, j) = rtmp rtmp = lon_ext(ii(i),j) lon_ext(ii(i),j) = lon_ext(i,j) lon_ext(i, j) = rtmp END DO END DO END IF 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(ii(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(ii(i),j,2,1)/g END IF IF ( var_nr3d(1,2) == 0 ) THEN tsoil_ext (i,j,:) = -999.0 qsoil_ext (i,j,:) = -999.0 ELSE IF(soilmodel_option == 1) THEN ! ARPS: two-layer force-restore model) tsoil_ext(i,j,1) = var_grb2d(ii(i),j,3,1) ! sfc temp. IF ( nint(var_grb2d(ii(i),j,5,1)) == 1 ) THEN ! soil temp over land tsoil_ext(i,j,2) = var_grb3d(ii(i),j,1,1,2) IF ( tsoil_ext (i,j,2) <= 200. ) THEN tsoil_ext(i,j,2) = tsoil_ext(i,j,1) END IF qsoil_ext(i,j,1) = var_grb3d(ii(i),j,2,2,2) qsoil_ext(i,j,2) = var_grb3d(ii(i),j,1,2,2) soiltyp_ext(i,j)= 0 ! We do not know the soil type ELSE ! sfc temp over sea tsoil_ext(i,j,2) = tsoil_ext(i,j,1) qsoil_ext(i,j,:) = 1.0 soiltyp_ext(i,j) = 13 ! Water END IF ELSE ! ARPS: multi-layer OUSoil model tsoil_ext(i,j,1) = var_grb2d(ii(i),j,3,1) ! sfc temp. IF ( nint(var_grb2d(ii(i),j,5,1)) == 1 ) THEN ! soil temp over land tsoil_ext(i,j,2) = var_grb3d(ii(i),j,2,1,2) ! 0 - 10 cm tsoil_ext(i,j,3) = var_grb3d(ii(i),j,1,1,2) ! 10 - 200 cm qsoil_ext(i,j,1) = var_grb3d(ii(i),j,2,2,2) ! 0 - 10 cm ! Value at surface assume equals ! the first layer below ground qsoil_ext(i,j,2) = var_grb3d(ii(i),j,2,2,2) ! 0 - 10 cm qsoil_ext(i,j,3) = var_grb3d(ii(i),j,1,2,2) ! 10 - 200 cm soiltyp_ext(i,j)= 0 ! We do not know the soil type ELSE ! sfc temp over sea tsoil_ext(i,j,2:3) = tsoil_ext(i,j,1) ! over water, All soil layers = SST qsoil_ext(i,j,1:3) = 1.0 ! All soil layers = 1.0 soiltyp_ext(i,j) = 13 ! Water END IF END IF ! soilmodel_option END IF IF ( var_nr2d(4,1) == 0 ) THEN wetcanp_ext(i,j) = -999.0 ELSE wetcanp_ext(i,j) = var_grb2d(ii(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(ii(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(ii(i),j,k,1,1) u_ext (i,j,kk) = var_grb3d(ii(i),j,k,2,1) ! u wind (m/s) v_ext (i,j,kk) = var_grb3d(ii(i),j,k,3,1) ! v wind (m/s) t_ext (i,j,kk) = var_grb3d(ii(i),j,k,4,1) ! Temperature (K) qc_ext (i,j,kk) = var_grb3d(ii(i),j,k,7,1) 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(ii(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 getncepavn2 !################################################################## !################################################################## !###### ###### !###### 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=256) :: 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 ! !################################################################## !################################################################## !###### ###### !###### SUBROUTINE GETNNARR221 ###### !###### ###### !###### Developed by ###### !###### Atmospheric Science Division ###### !###### Lawrence Livermore National Laboratory ###### !###### ###### !################################################################## !################################################################## ! SUBROUTINE getnarr221(nx_ext,ny_ext,nz_ext,nzsoil_ext, & 1,10 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, & t_2m_ext, qv_2m_ext, u_10m_ext, v_10m_ext, & istatus) ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! Reads a NMC GRIB (Grid #221, 32km) ETA file from the North American ! Regional Reanalysis (NARR) data for processing by ext2arps, a program ! that converts external files to ARPS variables and format. ! The data are in GRIB version 1 format. ! ! NOTE: ! ! The description about Eta 32km NARR output is available at ! http://wwwt.emc.ncep.noaa.gov/mmb/rreanl/index.html ! ! The data can be downloaded from ! http://nomads.ncdc.noaa.gov/#narr_datasets ! ! This subroutine assumes that data has been downloaded to local ! disk and stored in "dir_extd"/"extdname"_YYYYMMDD_HHHH_000.grb ! ! where: ! ! o "dir_extd" is specified inside arps.input; ! o "extdname" supposed to be "narr-a_221"; ! o CC is the model cycle, extracted from "extdinit" (00,06,12,18); ! o HHHH is the forecast hour, extracted from "extdfcst"; ! o YYYY is the year ! o MM is the month ! o DD is the day ! e.g. narr-a_221_20030701_0300_000.grb ! !----------------------------------------------------------------------- ! ! AUTHOR: Tina Katopodes Chow (04/18/2005) ! modified from getnmceta212 to create getnmceta221 ! ! MODIFICATION HISTORY: ! ! 2005/06/23 (Y. Wang) ! Incorporated into ARPS package and renamed to getnarr221 as well did ! some tests. ! !----------------------------------------------------------------------- ! ! 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, does not used here ! 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) ! ! T_2m_ext Temperature at 2m AGL ! qv_2m_ext Specific Humidity at 2m AGL ! U_10m_ext U at 10m AGL ! V_10m_ext V at 10m AGL ! ! 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 ! does not use CHARACTER (LEN=9) :: julfname ! does not use ! !----------------------------------------------------------------------- ! ! 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) ! Soil type REAL :: t_2m_ext (nx_ext,ny_ext) REAL :: qv_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) 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=256) :: gribfile CHARACTER (LEN=13) :: gribtime CHARACTER (LEN=17) :: gribtimetmp INTEGER :: i,j,k,kk INTEGER :: iyr,imo,iday 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(5) ! EMK 15 June 2002 INTEGER :: chklev, lvscan INTEGER :: iret ! Return flag INTEGER :: terrainflag ! !----------------------------------------------------------------------- ! ! 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)) ! Define depths of soil layers DATA soildepth_ext/0.0,-0.05,-0.25,-0.70,-1.50/ ! Parse date into year, month, day, etc... READ (extdinit,'(i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,i2)') & iyr,imo,iday,ihr,imin,isec ifhr=0 ! fixed for NARR reanalysis data ifmin=0 ifsec=0 ! gribtime is used to store the time in the required format for later 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 len1 = grbflen lenrun = LEN( extdname ) CALL strlnth( extdname, lenrun ) !TINA - use gribtimetmp only for setting the file name. WRITE (gribtimetmp,'(i4.4,i2.2,i2.2,a1,i2.2,a)') & iyr,imo,iday,'_',ihr,'00_000' gribfile = dir_extd(1:grbflen)//extdname(1:lenrun) & //'_'//gribtimetmp(1:17)//'.grb' grbflen = grbflen + lenrun + 22 ! !----------------------------------------------------------------------- ! ! RDNMCGRB reads NMC GRIB data ! !----------------------------------------------------------------------- ! govrd = g/rd !Get grid type and projection type gridtyp = narr221grid mproj_grb = narr221proj !Get number and names of 2d variables n2dvars = narr221nvs2d n2dlvtps = narr221nlvt2d DO k=1,n2dlvtps DO n=1,n2dvars var_id2d(n,k) = narr221var_id2d(n,k) END DO levtyp2d(k) = narr221levs2d(k) END DO !Get number and names of 3d variables n3dvars = narr221nvs3d n3dlvtps = narr221nlvt3d DO m=1,n3dlvtps DO n=1,n3dvars var_id3d(n,m) = narr221var_id3d(n,m) END DO levtyp3d(m) = narr221levs3d(m) END DO !Get original map projection for restore purpose CALL getmapr(iproj,scale,latnot,trlon,x0,y0) !Read fixed field file CALL rdnmcgrb(nx_ext,ny_ext,nz_ext,dir_extd(1:len1)//'AWIP32.fixed', & len1+12,'1979110800f00', & 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 IF ( max_nr2d == 0 ) THEN WRITE (6,'(a)') & 'No 2-D fixed field was found in the GRIB file, AWIP32.fixed.' END IF terrainflag = var_nr2d(2,1) !Read grib data print *,'Reading gribfile = ', TRIM(gribfile) 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) 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.', & 'Dataset problem in GETNARR221.', & ' ' ISTATUS = -888 GOTO 999 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 k=1,max_nr3d ! WRITE (6,'(/i5,4x,6(i7))') & ! k,(var_lev3d(k,n,1),n=1,n3dvars) ! END DO 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.', & 'Dataset problem in GETNARR221.', & ' ' ISTATUS = -888 GOTO 999 END IF END DO END DO !Set projection and grid spacing 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 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 ! first level type, first variable psfc_ext (i,j) = -999.0 ELSE !psfc_ext (i,j) = var_grb2d(i,j,1,1) * 100.0 psfc_ext (i,j) = var_grb2d(i,j,1,1) ! already with Pa unit (F.KONG) END IF IF ( terrainflag == 0 ) THEN trn_ext (i,j) = -999.0 WRITE(6,'(/a/)') 'WARNING: no terrain found in grib file .... ' ELSE trn_ext (i,j) = var_grb2d(i,j,2,1) END IF IF ( var_nr3d(1,2) == 0 ) THEN ! do not contain soil temperature DO k=1,nzsoil_ext tsoil_ext (i,j,k) = -999.0 qsoil_ext (i,j,k) = -999.0 END DO ELSE ! yes, it contain for NARR 221 IF (soilmodel_option == 1) THEN ! Old ARPS Force-Restore Soil Model tsoil_ext(i,j,1) = var_grb2d(i,j,3,1) ! surface temperature 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 ! should contain 100-200cm here? 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 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,3,1) ! Water temperature qsoil_ext (i,j,k) = 1. ! 100% water END DO END IF ! Land or water? END IF ! soilmodel_option END IF IF ( var_nr2d(4,1) == 0 ) THEN wetcanp_ext(i,j) = -999.0 ELSE ! Plant canopy surface water (kg/m**2) 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 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 IF( var_nr2d(2,2) == 0.) THEN qv_2m_ext(i,j)= -999.0 ELSE qv_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 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,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) ! vertical velocity (pressure) - skipped. qc_ext (i,j,kk) = var_grb3d(i,j,k,7,1) ! height (m) qi_ext (i,j,kk) = var_grb3d(i,j,k,8,1) ! height (m) 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 RUC data are sent as grid-relative. ! !----------------------------------------------------------------------- ! !NCEP NARR output is already earth relative - TINA istatus = 1 ! !----------------------------------------------------------------------- ! ! Reset map projection to previous values ! !----------------------------------------------------------------------- ! 999 CONTINUE CALL setmapr(iproj,scale,latnot,trlon) CALL setorig(1,x0,y0) DEALLOCATE(var_grb2d,var_grb3d,rcdata,var_lev3d) RETURN END SUBROUTINE getnarr221