!########################################################################
!########################################################################
!###### ######
!###### PROGRAM ARPSVERIF ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!########################################################################
!########################################################################
PROGRAM arpsverif,120
!-----------------------------------------------------------------------
!
! PURPOSE:
!
! Calculates verification statistics for ARPS forecasts.
!
!-----------------------------------------------------------------------
!
! HISTORY:
! 2002/04/03 First written by Jason J. Levit
!
! 2005/08/05 MPI version. Kevin W. Thomas
!
!-----------------------------------------------------------------------
!
! Use modules
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
IMPLICIT NONE ! Force explicit declarations
INCLUDE 'globcst.inc'
INCLUDE 'phycst.inc'
INCLUDE 'vericst.inc'
INCLUDE 'grid.inc'
INCLUDE 'mp.inc'
!
!-----------------------------------------------------------------------
!
! Data arrays
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: x (:) ! The x-coord. of the physical and
! computational grid. Defined at u-point.
REAL, ALLOCATABLE :: y (:) ! The y-coord. of the physical and
! computational grid. Defined at v-point.
REAL, ALLOCATABLE :: z (:) ! The z-coord. of the computational grid.
! Defined at w-point on the staggered grid.
REAL, ALLOCATABLE :: zp (:,:,:) ! The physical height coordinate defined at
! w-point of the staggered grid.
REAL, ALLOCATABLE :: zpsoil(:,:,:) ! The physical height coordinate of soil
! model
!
REAL, ALLOCATABLE :: hterain (:,:) ! Terrain height
!
REAL, ALLOCATABLE :: uprt (:,:,:) ! Perturbation u-velocity (m/s)
REAL, ALLOCATABLE :: vprt (:,:,:) ! Perturbation v-velocity (m/s)
REAL, ALLOCATABLE :: wprt (:,:,:) ! Perturbation w-velocity (m/s)
REAL, ALLOCATABLE :: ptprt (:,:,:) ! Perturbation potential temperature (K)
REAL, ALLOCATABLE :: pprt (:,:,:) ! Perturbation pressure (Pascal)
REAL, ALLOCATABLE :: qvprt (:,:,:) ! Perturbation water vapor specific
! humidity (kg/kg)
REAL, ALLOCATABLE :: qc (:,:,:) ! Cloud water mixing ratio (kg/kg)
REAL, ALLOCATABLE :: qr (:,:,:) ! Rain water mixing ratio (kg/kg)
REAL, ALLOCATABLE :: qi (:,:,:) ! Cloud ice mixing ratio (kg/kg)
REAL, ALLOCATABLE :: qs (:,:,:) ! Snow mixing ratio (kg/kg)
REAL, ALLOCATABLE :: qh (:,:,:) ! Hail mixing ratio (kg/kg)
REAL, ALLOCATABLE :: tke (:,:,:) ! Turbulent Kinetic Energy ((m/s)**2)
REAL, ALLOCATABLE :: kmh (:,:,:) ! Horizontal turb. mixing coef. for
! momentum. ( m**2/s )
REAL, ALLOCATABLE :: kmv (:,:,:) ! Vertical turb. mixing coef. for
! momentum. ( m**2/s )
REAL, ALLOCATABLE :: ubar (:,:,:) ! Base state u-velocity (m/s)
REAL, ALLOCATABLE :: vbar (:,:,:) ! Base state v-velocity (m/s)
REAL, ALLOCATABLE :: wbar (:,:,:) ! Base state w-velocity (m/s)
REAL, ALLOCATABLE :: ptbar (:,:,:) ! Base state potential temperature (K)
REAL, ALLOCATABLE :: pbar (:,:,:) ! Base state pressure (Pascal)
REAL, ALLOCATABLE :: rhobar (:,:,:) ! Base state air density (kg/m**3)
REAL, ALLOCATABLE :: qvbar (:,:,:) ! Base state water vapor specific
! humidity (kg/kg)
INTEGER :: nstyps
INTEGER, ALLOCATABLE :: soiltyp (:,:,:) ! Soil type
REAL, ALLOCATABLE :: stypfrct(:,:,:) ! Soil type
INTEGER, ALLOCATABLE :: vegtyp (:,:) ! Vegetation type
REAL, ALLOCATABLE :: lai (:,:) ! Leaf Area Index
REAL, ALLOCATABLE :: roufns (:,:) ! Surface roughness
REAL, ALLOCATABLE :: veg (:,:) ! Vegetation fraction
REAL, ALLOCATABLE :: tsoil (:,:,:,:) ! soil temperature (K)
REAL, ALLOCATABLE :: qsoil (:,:,:,:) ! soil moisture
REAL, ALLOCATABLE :: wetcanp(:,:,:) ! Canopy water amount
REAL, ALLOCATABLE :: snowdpth(:,:) ! Snow depth (m)
REAL, ALLOCATABLE :: raing(:,:) ! Grid supersaturation rain
REAL, ALLOCATABLE :: rainc(:,:) ! Cumulus convective rain
REAL, ALLOCATABLE :: prcrate(:,:,:) ! 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) = cumulative precip. rate
! prcrate(1,1,4) = microphysics precip. rate
REAL, ALLOCATABLE :: radfrc(:,:,:) ! Radiation forcing (K/s)
REAL, ALLOCATABLE :: radsw (:,:) ! Solar radiation reaching the surface
REAL, ALLOCATABLE :: rnflx (:,:) ! Net radiation flux absorbed by surface
REAL, ALLOCATABLE :: radswnet(:,:) ! Net shortwave radiation
REAL, ALLOCATABLE :: radlwin(:,:) ! Incoming longwave radiation
REAL, ALLOCATABLE :: usflx (:,:) ! Surface flux of u-momentum (kg/(m*s**2))
REAL, ALLOCATABLE :: vsflx (:,:) ! Surface flux of v-momentum (kg/(m*s**2))
REAL, ALLOCATABLE :: ptsflx(:,:) ! Surface heat flux (K*kg/(m*s**2))
REAL, ALLOCATABLE :: qvsflx(:,:) ! Surface moisture flux (kg/(m**2*s)
REAL, ALLOCATABLE :: model_data(:,:,:)
REAL, ALLOCATABLE :: obsrv_data(:,:,:)
!(:,:,1) = sfc temp
!(:,:,2) = sfc dewpoint
!(:,:,3) = sfc wind dir
!(:,:,4) = sfc wind speed
!(:,:,5) = sfc pressure
REAL, ALLOCATABLE :: stats_data(:,:,:)
!(:,1,1) = RMSE data (temp)
!(:,1,2) = RMSE data (dewpoint)
!(:,1,3) = RMSE data (sfc wind dir)
!(:,1,4) = RMSE data (sfc wind speed)
!(:,1,5) = RMSE data (sfc pressure)
!(:,2,1) = Bias data (temp)
!(:,2,2) = Bias data (dewpoint)
!(:,2,3) = Bias data (sfc wind dir)
!(:,2,4) = Bias data (sfc wind speed)
!(:,2,5) = Bias data (sfc pressure)
! etc.
!
!-----------------------------------------------------------------------
!
! Namelists
!
!-----------------------------------------------------------------------
!
INTEGER :: verifopt
NAMELIST /verif_opt/ verifopt
INTEGER :: sndopt
CHARACTER (LEN=256) :: sndlist
CHARACTER (LEN=256) :: sndrunname
CHARACTER (LEN=256) :: snddomlist
INTEGER :: proopt
CHARACTER (LEN=256) :: prolist
CHARACTER (LEN=256) :: prorunname
CHARACTER (LEN=256) :: prodomlist
INTEGER :: sfcopt
CHARACTER (LEN=256) :: sfclist
CHARACTER (LEN=256) :: sfcobsdir
CHARACTER (LEN=256) :: sfcrunname
CHARACTER (LEN=256) :: blackfile
INTEGER :: precveropt
CHARACTER (LEN=256) :: preclist
CHARACTER (LEN=256) :: precrunname
CHARACTER (LEN=256) :: precdomlist
INTEGER :: mosopt
CHARACTER (LEN=256) :: moslist
CHARACTER (LEN=256) :: mosrunname
CHARACTER (LEN=256) :: mosdomlist
INTEGER :: arpsnn_opt
CHARACTER (LEN=256) :: nnrunname
CHARACTER (LEN=256) :: wtsdir
CHARACTER (LEN=256) :: nnoutputfn
NAMELIST /single/ sndopt,sndlist,sndrunname,snddomlist, &
proopt,prolist,prorunname,prodomlist, &
sfcopt,sfclist,sfcobsdir, &
sfcrunname,blackfile,nsfcuse,sfcuse, &
precveropt,preclist,precrunname,precdomlist, &
mosopt,moslist,mosrunname,mosdomlist, &
arpsnn_opt,nnrunname,wtsdir,nnoutputfn
INTEGER :: dummy
NAMELIST /gridded/ dummy
NAMELIST /message_passing/ nproc_x,nproc_y,readsplit
!
!-----------------------------------------------------------------------
!
! Work Arrays
!
!-----------------------------------------------------------------------
!
REAL, ALLOCATABLE :: tem1(:,:,:)
REAL, ALLOCATABLE :: tem2(:,:,:)
REAL, ALLOCATABLE :: tem3(:,:,:)
REAL, ALLOCATABLE :: tem4(:,:)
!
!-----------------------------------------------------------------------
!
! Misc. internal variables
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,k
INTEGER :: nx,ny,nz
INTEGER :: nzsoil
INTEGER :: ireturn, istatus
INTEGER :: hinfmt
INTEGER :: nhisfile_max,nhisfile
PARAMETER (nhisfile_max=200)
CHARACTER (LEN=256) :: grdbasfn
CHARACTER (LEN=256) :: filename
CHARACTER (LEN=256) :: hisfile(nhisfile_max)
CHARACTER (LEN=4), PARAMETER :: tail=".hdf"
INTEGER :: lengbf,nf,lenfil
INTEGER, PARAMETER :: max_dim=200
INTEGER, PARAMETER :: unum=5 ! unit number for namelist read
INTEGER :: istat,sd_id,sds_id
CHARACTER (LEN=256) :: hdfname
INTEGER :: dims(2)
INTEGER :: lname
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
nsfcuse = 0
CALL mpinit_proc
IF (myproc == 0) THEN
WRITE(6,'(/,8(5x,a,/),/)') &
'###############################################################',&
'###############################################################',&
'### ###',&
'### Welcome to ARPSVERIF ###',&
'### A program that computes verification of ARPS forecasts. ###',&
'### ###',&
'###############################################################',&
'###############################################################'
READ (unum,message_passing)
WRITE(6,'(a)')'Namelist block message_passing sucessfully read.'
READ (unum,verif_opt)
WRITE(6,'(a)')'Namelist block verif_opt sucessfully read.'
CALL get_input_file_names(hinfmt,grdbasfn,hisfile,nhisfile)
READ (unum,single)
WRITE(6,'(a)')'Namelist block single sucessfully read.'
READ (unum,gridded)
WRITE(6,'(a)')'Namelist block gridded sucessfully read.'
filename = grdbasfn
lengbf = len_trim(filename)
IF ( mp_opt > 0 .AND. readsplit <= 0 ) THEN
filename(lengbf+1:lengbf+5) = '_0101'
lengbf = lengbf + 5
END IF
CALL get_dims_from_data(hinfmt,filename(1:lengbf), &
nx,ny,nz,nzsoil,nstyps, ireturn)
IF( ireturn /= 0 ) THEN
PRINT*,'Problem occured when trying to get dimensions from data.'
PRINT*,'Program stopped.'
CALL arpsstop("Fail in get_dims_from_data",1)
END IF
WRITE(6,'(3(a,i5))') 'nx =',nx,', ny=',ny,', nz=',nz
print*,'nstyps =', nstyps
END IF
CALL mpupdatei(nproc_x,1)
CALL mpupdatei(nproc_y,1)
CALL mpupdatei(readsplit,1)
CALL mpupdatei(nx,1)
CALL mpupdatei(ny,1)
CALL mpupdatei(nz,1)
CALL mpupdatei(nzsoil,1)
CALL mpupdatei(nstyps,1)
CALL mpupdatei(nhisfile,1)
CALL mpupdatei(hinfmt,1)
CALL mpupdatec(grdbasfn,256)
CALL mpupdatec(hisfile,nhisfile*256)
CALL mpupdatei(sndopt,1)
CALL mpupdatei(proopt,1)
CALL mpupdatei(sfcopt,1) !sfcopt ignored, always enabled
CALL mpupdatei(nsfcuse,1)
CALL mpupdatec(sfcuse,nsfcuse*4)
CALL mpupdatei(precveropt,1)
CALL mpupdatei(mosopt,1)
CALL mpupdatei(arpsnn_opt,1)
CALL mpupdatei(verifopt,1)
IF ( mp_opt == 0 ) THEN
nproc_x = 1
nproc_y = 1
readsplit = 0
END IF
! call mpinit_var
!
!-----------------------------------------------------------------------
!
! Single-point verification
!
!-----------------------------------------------------------------------
!
IF (verifopt.eq.1.or.verifopt.eq.3) THEN
IF ( mp_opt > 0 .AND. readsplit > 0 ) THEN
nx = (nx - 3) / nproc_x + 3
ny = (ny - 3) / nproc_y + 3
IF ( myproc == 0 ) WRITE(6,*) 'Processor nx/ny now',nx,ny
END IF
CALL mpinit_var
!
! ALLOCATE arrays
!
ALLOCATE(x (nx),STAT=istatus)
CALL check_alloc_status(istatus, "arps:x")
x = 0.0
ALLOCATE(y (ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:y")
y = 0.0
ALLOCATE(z (nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:z")
z = 0.0
ALLOCATE(zp (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:zp")
zp = 0.0
ALLOCATE(zpsoil (nx,ny,nzsoil),STAT=istatus)
CALL check_alloc_status(istatus, "arps:zpsoil")
zpsoil = 0.0
ALLOCATE(uprt (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:uprt")
uprt = 0.0
ALLOCATE(vprt (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:vprt")
vprt = 0.0
ALLOCATE(wprt (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:wprt")
wprt = 0.0
ALLOCATE(hterain (nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:hterain")
hterain = 0.0
ALLOCATE(ptprt (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:ptprt")
ptprt = 0.0
ALLOCATE(pprt (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:pprt")
pprt = 0.0
ALLOCATE(qvprt (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:qvprt")
qvprt = 0.0
ALLOCATE(qc (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:qc")
qc = 0.0
ALLOCATE(qr (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:qr")
qr = 0.0
ALLOCATE(qi (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:qi")
qi = 0.0
ALLOCATE(qs (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:qs")
qs = 0.0
ALLOCATE(qh (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:qh")
qh = 0.0
ALLOCATE(tke (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:tke")
tke = 0.0
ALLOCATE(kmh (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:kmh")
kmh = 0.0
ALLOCATE(kmv (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:kmv")
kmv = 0.0
ALLOCATE(ubar (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:ubar")
ubar = 0.0
ALLOCATE(vbar (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:vbar")
vbar = 0.0
ALLOCATE(wbar (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:wbar")
wbar = 0.0
ALLOCATE(ptbar (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:ptbar")
ptbar = 0.0
ALLOCATE(pbar (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:pbar")
pbar = 0.0
ALLOCATE(rhobar (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:rhobar")
rhobar = 0.0
ALLOCATE(qvbar (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:qvbar")
qvbar = 0.0
ALLOCATE(soiltyp (nx,ny,nstyps),STAT=istatus)
CALL check_alloc_status(istatus, "arps:soiltyp")
soiltyp = 0.0
ALLOCATE(stypfrct (nx,ny,nstyps),STAT=istatus)
CALL check_alloc_status(istatus, "arps:stypfrct")
stypfrct = 0.0
ALLOCATE(vegtyp (nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:vegtyp")
vegtyp = 0.0
ALLOCATE(lai (nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:lai")
lai = 0.0
ALLOCATE(roufns (nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:roufns")
roufns = 0.0
ALLOCATE(veg (nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:veg")
veg = 0.0
ALLOCATE(tsoil (nx,ny,nzsoil,0:nstyps),STAT=istatus)
CALL check_alloc_status(istatus, "arps:tsoil")
tsoil = 0.0
ALLOCATE(qsoil (nx,ny,nzsoil,0:nstyps),STAT=istatus)
CALL check_alloc_status(istatus, "arps:qsoil")
qsoil = 0.0
ALLOCATE(wetcanp (nx,ny,0:nstyps),STAT=istatus)
CALL check_alloc_status(istatus, "arps:wetcanp")
wetcanp = 0.0
ALLOCATE(snowdpth (nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:snowdpth")
snowdpth = 0.0
ALLOCATE(raing (nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:raing")
raing = 0.0
ALLOCATE(rainc (nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:rainc")
rainc = 0.0
ALLOCATE(prcrate (nx,ny,4),STAT=istatus)
CALL check_alloc_status(istatus, "arps:prcrate")
prcrate = 0.0
ALLOCATE(radfrc (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:radfrc")
radfrc = 0.0
ALLOCATE(radsw (nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:radsw")
radsw = 0.0
ALLOCATE(rnflx (nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:rnflx")
rnflx = 0.0
ALLOCATE(radswnet (nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:radswnet")
radswnet = 0.0
ALLOCATE(radlwin (nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:radlwin")
radlwin = 0.0
ALLOCATE(usflx (nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:usflx")
usflx = 0.0
ALLOCATE(vsflx (nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:vsflx")
vsflx = 0.0
ALLOCATE(ptsflx (nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:ptsflx")
ptsflx = 0.0
ALLOCATE(qvsflx (nx,ny),STAT=istatus)
CALL check_alloc_status(istatus, "arps:qvsflx")
qvsflx = 0.0
ALLOCATE(tem1 (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:tem1")
tem1 = 0.0
ALLOCATE(tem2 (nx,ny,nz),STAT=istatus)
CALL check_alloc_status(istatus, "arps:tem2")
tem2 = 0.0
ALLOCATE(model_data(sfcmax,nhisfile,5),STAT=istatus)
CALL check_alloc_status(istatus, "arps:model_data")
model_data=0.0
ALLOCATE(obsrv_data(sfcmax,nhisfile,5),STAT=istatus)
CALL check_alloc_status(istatus, "arps:obsrv_data")
obsrv_data=-99.9 ! Set all obs equal to an initial missing value
ALLOCATE(stats_data(nhisfile,2,5),STAT=istatus)
CALL check_alloc_status(istatus, "arps:stats_data")
stats_data=0.0
readstns=.false.
IF (myproc == 0) THEN
IF (sndopt.eq.1) THEN
write (6,*) "Sounding verification is currently unavailable."
write (6,*) "Program will terminate."
CALL arpsstop("Configuration Error",1)
END IF
IF (proopt.eq.1) THEN
write (6,*) "Profiler verification is currently unavailable."
write (6,*) "Program will terminate."
CALL arpsstop("Configuration Error",1)
END IF
IF (precveropt.eq.1) THEN
write (6,*) "Precip. verification is currently unavailable."
write (6,*) "Program will terminate."
CALL arpsstop("Configuration Error",1)
END IF
IF (mosopt.eq.1) THEN
write (6,*) "MOS computation is currently unavailable."
write (6,*) "Program will terminate."
CALL arpsstop("Configuration Error",1)
END IF
IF (arpsnn_opt.eq.1) THEN
write (6,*) "The ARPS Neural Network option is currently unavailable."
write (6,*) "Program will terminate."
CALL arpsstop("Configuration Error",1)
END IF
ENDIF
CALL his2ver(nx,ny,nz,nzsoil,x,y,z,zp,zpsoil,hterain, &
uprt,vprt,wprt,ptprt,pprt,qvprt, &
qc,qr,qi,qs,qh,tke,kmh,kmv, &
ubar,vbar,wbar,ptbar,pbar,rhobar,qvbar, &
nstyps,soiltyp,stypfrct,vegtyp,lai,roufns, &
veg,tsoil,qsoil,wetcanp,snowdpth, &
raing,rainc,prcrate,radfrc,radsw,rnflx, &
radswnet,radlwin, &
usflx,vsflx,ptsflx,qvsflx, &
hinfmt,grdbasfn,hisfile,nhisfile, &
sndopt,sndlist,sndrunname,snddomlist, &
proopt,prolist,prorunname,prodomlist, &
sfcopt,sfclist,sfcobsdir,sfcrunname,blackfile, &
precveropt,preclist,precrunname,precdomlist, &
mosopt,moslist,mosrunname,mosdomlist, &
arpsnn_opt,nnrunname,wtsdir,nnoutputfn, &
model_data,obsrv_data)
!
! Collect output if we're MPI.
!
IF (mp_opt > 0 ) THEN
CALL verif_collect(model_data,obsrv_data,nhisfile)
END IF
IF (myproc == 0) THEN
!
! Calculate statistics on surface data.
!
DO i=1,nhisfile
! RSME calculations
tem1(1:sfcstn,1,1)=obsrv_data(1:sfcstn,i,1)
tem2(1:sfcstn,1,1)=model_data(1:sfcstn,i,1)
CALL root_mean_square(sfcstn,1,1,tem1,tem2,0,0,0,stats_data(i,1,1))
tem1=0.0
tem2=0.0
tem1(1:sfcstn,1,1)=obsrv_data(1:sfcstn,i,2)
tem2(1:sfcstn,1,1)=model_data(1:sfcstn,i,2)
CALL root_mean_square(sfcstn,1,1,tem1,tem2,0,0,0,stats_data(i,1,2))
tem1=0.0
tem2=0.0
tem1(1:sfcstn,1,1)=obsrv_data(1:sfcstn,i,3)
tem2(1:sfcstn,1,1)=model_data(1:sfcstn,i,3)
CALL root_mean_square(sfcstn,1,1,tem1,tem2,0,0,0,stats_data(i,1,3))
tem1=0.0
tem2=0.0
tem1(1:sfcstn,1,1)=obsrv_data(1:sfcstn,i,4)
tem2(1:sfcstn,1,1)=model_data(1:sfcstn,i,4)
CALL root_mean_square(sfcstn,1,1,tem1,tem2,0,0,0,stats_data(i,1,4))
tem1=0.0
tem2=0.0
tem1(1:sfcstn,1,1)=obsrv_data(1:sfcstn,i,5)
tem2(1:sfcstn,1,1)=model_data(1:sfcstn,i,5)
CALL root_mean_square(sfcstn,1,1,tem1,tem2,0,0,0,stats_data(i,1,5))
! Bias calculations
tem1=0.0
tem2=0.0
tem1(1:sfcstn,1,1)=obsrv_data(1:sfcstn,i,1)
tem2(1:sfcstn,1,1)=model_data(1:sfcstn,i,1)
CALL bias_calc(sfcstn,1,1,tem1,tem2,0,0,0,stats_data(i,2,1))
tem1=0.0
tem2=0.0
tem1(1:sfcstn,1,1)=obsrv_data(1:sfcstn,i,2)
tem2(1:sfcstn,1,1)=model_data(1:sfcstn,i,2)
CALL bias_calc(sfcstn,1,1,tem1,tem2,0,0,0,stats_data(i,2,2))
tem1=0.0
tem2=0.0
tem1(1:sfcstn,1,1)=obsrv_data(1:sfcstn,i,3)
tem2(1:sfcstn,1,1)=model_data(1:sfcstn,i,3)
CALL bias_calc(sfcstn,1,1,tem1,tem2,0,0,0,stats_data(i,2,3))
tem1=0.0
tem2=0.0
tem1(1:sfcstn,1,1)=obsrv_data(1:sfcstn,i,4)
tem2(1:sfcstn,1,1)=model_data(1:sfcstn,i,4)
CALL bias_calc(sfcstn,1,1,tem1,tem2,0,0,0,stats_data(i,2,4))
tem1=0.0
tem2=0.0
tem1(1:sfcstn,1,1)=obsrv_data(1:sfcstn,i,5)
tem2(1:sfcstn,1,1)=model_data(1:sfcstn,i,5)
CALL bias_calc(sfcstn,1,1,tem1,tem2,0,0,0,stats_data(i,2,5))
END DO
!
! Write out all verification data into a specific HDF file.
!
lname=LEN_TRIM(sfcrunname)
hdfname=sfcrunname(1:lname)//".hdf"
CALL hdfopen(hdfname,2,sd_id)
dims(2)=sfcstn
dims(1)=len(sfcstid(1))
CALL write_sds_char(sds_id,sd_id,'sfcstid',sfcstid,2,dims)
CALL hdfwrt2d(model_data(:,:,1),sfcmax,nhisfile,sd_id,0,0,'model_temp', &
'2m ARPS Temperature','F',istatus)
CALL hdfwrt2d(obsrv_data(:,:,1),sfcmax,nhisfile,sd_id,0,0,'obs_temp', &
'Observation Temperature','F',istatus)
CALL hdfwrt2d(model_data(:,:,2),sfcmax,nhisfile,sd_id,0,0,'model_dewp', &
'2m ARPS Dewpoint','F',istatus)
CALL hdfwrt2d(obsrv_data(:,:,2),sfcmax,nhisfile,sd_id,0,0,'obs_dewp', &
'Observation Dewpoint','F',istatus)
CALL hdfwrt2d(model_data(:,:,3),sfcmax,nhisfile,sd_id,0,0,'model_wdir', &
'ARPS surface wind direction','m/s',istatus)
CALL hdfwrt2d(obsrv_data(:,:,3),sfcmax,nhisfile,sd_id,0,0,'obs_wdir', &
'Observation surface wind direction','m/s',istatus)
CALL hdfwrt2d(model_data(:,:,4),sfcmax,nhisfile,sd_id,0,0,'model_wspd', &
'ARPS surface wind speed','m/s',istatus)
CALL hdfwrt2d(obsrv_data(:,:,4),sfcmax,nhisfile,sd_id,0,0,'obs_wspd', &
'Observation surface wind speed','m/s',istatus)
CALL hdfwrt2d(model_data(:,:,5),sfcmax,nhisfile,sd_id,0,0,'model_pres', &
'ARPS surface pressure','mb',istatus)
CALL hdfwrt2d(obsrv_data(:,:,5),sfcmax,nhisfile,sd_id,0,0,'obs_pres', &
'Observation surface pressure','mb',istatus)
CALL hdfwrt1d(stats_data(:,1,1),nhisfile,sd_id,'rmse_temp', &
'2m Surface Temperature RMSE','none')
CALL hdfwrt1d(stats_data(:,1,2),nhisfile,sd_id,'rmse_dewp', &
'2m Surface Dewpoint RMSE','none')
CALL hdfwrt1d(stats_data(:,1,3),nhisfile,sd_id,'rmse_wdir', &
'Surface Wind Direction RMSE','none')
CALL hdfwrt1d(stats_data(:,1,4),nhisfile,sd_id,'rmse_wspd', &
'Surface Wind Speed RMSE','none')
CALL hdfwrt1d(stats_data(:,1,5),nhisfile,sd_id,'rmse_pres', &
'Surface Pressure RMSE','none')
CALL hdfwrt1d(stats_data(:,2,1),nhisfile,sd_id,'bias_temp', &
'2m Surface Temperature Bias','none')
CALL hdfwrt1d(stats_data(:,2,2),nhisfile,sd_id,'bias_dewp', &
'2m Surface Dewpoint Bias','none')
CALL hdfwrt1d(stats_data(:,2,3),nhisfile,sd_id,'bias_wdir', &
'Surface Wind Direction Bias','none')
CALL hdfwrt1d(stats_data(:,2,4),nhisfile,sd_id,'bias_wspd', &
'Surface Wind Speed Bias','none')
CALL hdfwrt1d(stats_data(:,2,5),nhisfile,sd_id,'bias_pres', &
'Surface Pressure Bias','none')
CALL hdfclose(sd_id,istat)
END IF
END IF
IF (mp_opt > 0 ) CALL mpexit(0)
STOP
END PROGRAM arpsverif