PROGRAM ens_ana,77
IMPLICIT NONE
INCLUDE 'enscv.inc'
! INCLUDE 'mp.inc'
INTEGER :: nx,ny,nz
INTEGER :: i,j,ll
INTEGER :: istat,ierr,iSTATUS,nunit
INTEGER :: nmember,nf,anaopt
CHARACTER (LEN=256) :: ensdir,outdir
CHARACTER (LEN=256) :: memberheader(15), anaheader
CHARACTER (LEN=256) :: outheader
CHARACTER (LEN= 6) :: tlevel(50)
REAL :: sp(50,15) ! spread (nf, nvar)
REAL :: mse(50,15),sde(50,15),mnbias(50,15),sdbias(50,15),disp(50,15)
REAL :: corr(50,15),sd1(50,15),sd2(50,15)
REAL :: sum1(50,15),sum2(50,15)
REAL :: amax,amin
REAL :: totlegrid,totlegrid1,totlegrid2
NAMELIST /input_data/ ensdir,nmember,memberheader,nf,tlevel,anaopt,anaheader
NAMELIST /output_data/ outdir,outheader
REAL, ALLOCATABLE :: mslp(:,:)
REAL, ALLOCATABLE :: hgt250(:,:), vort250(:,:), uwind250(:,:), vwind250(:,:)
REAL, ALLOCATABLE :: hgt500(:,:), vort500(:,:), uwind500(:,:), vwind500(:,:)
REAL, ALLOCATABLE :: hgt850(:,:), vort850(:,:), uwind850(:,:), vwind850(:,:)
REAL, ALLOCATABLE :: temp850(:,:)
REAL, ALLOCATABLE :: thck(:,:), rh700(:,:), omega700(:,:),temp2m(:,:), dewp2m(:,:)
REAL, ALLOCATABLE :: accppt(:,:), convppt(:,:)
REAL, ALLOCATABLE :: sreh(:,:),li(:,:),cape(:,:),brn(:,:),cin(:,:),vhel(:,:),pw(:,:)
REAL, ALLOCATABLE :: cmpref(:,:),refl1km(:,:),refl4km(:,:),u10m(:,:),v10m(:,:)
REAL, ALLOCATABLE :: mslpa(:,:)
REAL, ALLOCATABLE :: hgt250a(:,:), vort250a(:,:), uwind250a(:,:), vwind250a(:,:)
REAL, ALLOCATABLE :: hgt500a(:,:), vort500a(:,:), uwind500a(:,:), vwind500a(:,:)
REAL, ALLOCATABLE :: hgt850a(:,:), vort850a(:,:), uwind850a(:,:), vwind850a(:,:)
REAL, ALLOCATABLE :: temp850a(:,:)
REAL, ALLOCATABLE :: thcka(:,:),rh700a(:,:),omega700a(:,:),temp2ma(:,:),dewp2ma(:,:)
REAL, ALLOCATABLE :: accppta(:,:), convppta(:,:)
REAL, ALLOCATABLE :: sreha(:,:),lia(:,:),capea(:,:),brna(:,:),cina(:,:)
REAL, ALLOCATABLE :: vhela(:,:),pwa(:,:)
REAL, ALLOCATABLE :: cmprefa(:,:),refl1kma(:,:),refl4kma(:,:),u10ma(:,:),v10ma(:,:)
REAL, ALLOCATABLE :: mslpm(:,:),hgt250m(:,:),hgt500m(:,:),hgt850m(:,:)
REAL, ALLOCATABLE :: temp850m(:,:),temp2mm(:,:),dewp2mm(:,:),accpptm(:,:),convpptm(:,:)
REAL, ALLOCATABLE :: cmprefm(:,:),refl1kmm(:,:),refl4kmm(:,:),u10mm(:,:),v10mm(:,:)
REAL, ALLOCATABLE :: mslpt(:,:),hgt250t(:,:),hgt500t(:,:),hgt850t(:,:)
REAL, ALLOCATABLE :: temp850t(:,:),temp2mt(:,:),dewp2mt(:,:),accpptt(:,:),convpptt(:,:)
REAL, ALLOCATABLE :: cmpreft(:,:),refl1kmt(:,:),refl4kmt(:,:),u10mt(:,:),v10mt(:,:)
! F.KONG add for probability calculation
REAL, ALLOCATABLE :: pbrain(:,:,:),pbref1(:,:,:),pbref4(:,:,:),pbcref(:,:,:)
REAL :: factor
! end F.KONG
CHARACTER (LEN=256) :: filenm,outfile
CHARACTER (LEN=132) :: char
CHARACTER (LEN=15) :: varname,unitnm
CHARACTER (LEN=40) :: q3dname,q3dunit
! REAL :: qswcorn, qswcorw, qnecorn, qnecorw
INTEGER :: hour,mimute,dweek,day,month,year
INTEGER :: nt,nm
! nproc_x = 1
! nproc_y = 1
! readsplit = 0
! joindmp = 1
! nprocs = 1
! max_fopen = 1
READ(5,input_data,END=100)
WRITE(6,'(/a/)')'Namelist block input_data sucessfully read in.'
WRITE(6,input_data)
READ(5,output_data,END=101)
WRITE(6,'(/a/)')'Namelist block output_data sucessfully read in.'
WRITE(6,output_data)
GO TO 111
100 WRITE(6,'(a)') &
'Error reading NAMELIST block input_data, STOP'
STOP
101 WRITE(6,'(a)') &
'Error reading NAMELIST block output_data, STOP'
STOP
111 CONTINUE
filenm = trim(ensdir)//trim(memberheader(1))//'.mspres'//tlevel(1)
print *, trim(filenm)
CALL getunit
( nunit)
OPEN(UNIT=nunit,FILE=trim(filenm),STATUS='old', FORM='unformatted', &
ERR=9000, IOSTAT=istat)
READ(nunit, ERR=9000, END=9000, IOSTAT=istat) nx,ny,nz
CLOSE(UNIT=nunit)
CALL retunit(nunit)
print *, 'nx,ny:',nx,ny
ALLOCATE(mslp(nx,ny),hgt250(nx,ny),vort250(nx,ny),uwind250(nx,ny),vwind250(nx,ny))
ALLOCATE(hgt500(nx,ny),vort500(nx,ny),uwind500(nx,ny),vwind500(nx,ny))
ALLOCATE(hgt850(nx,ny),vort850(nx,ny),uwind850(nx,ny),vwind850(nx,ny),temp850(nx,ny))
ALLOCATE(thck(nx,ny),rh700(nx,ny),omega700(nx,ny),temp2m(nx,ny),dewp2m(nx,ny))
ALLOCATE(accppt(nx,ny),convppt(nx,ny),sreh(nx,ny),li(nx,ny),cape(nx,ny))
ALLOCATE(brn(nx,ny),cin(nx,ny),vhel(nx,ny),pw(nx,ny),cmpref(nx,ny))
ALLOCATE(refl1km(nx,ny),refl4km(nx,ny),u10m(nx,ny),v10m(nx,ny))
! CALL binread2d(nx,ny,trim(filenm),varname,unitnm,mslp)
! print *, trim(varname),' ',trim(unitnm),nx,ny
!CALL a3dmax0(mslp,1,nx,1,nx-1,1,ny,1,ny,1,1,1,1,amax,amin)
!print *,'MSLP -- amax, amin: ',amax, amin
!stop
if(anaopt == 1) then ! allocate arrays for analysis fields
ALLOCATE(mslpa(nx,ny),hgt250a(nx,ny),vort250a(nx,ny),uwind250a(nx,ny),vwind250a(nx,ny))
ALLOCATE(hgt500a(nx,ny),vort500a(nx,ny),uwind500a(nx,ny),vwind500a(nx,ny))
ALLOCATE(hgt850a(nx,ny),vort850a(nx,ny),uwind850a(nx,ny),vwind850a(nx,ny),temp850a(nx,ny))
ALLOCATE(thcka(nx,ny),rh700a(nx,ny),omega700a(nx,ny),temp2ma(nx,ny),dewp2ma(nx,ny))
ALLOCATE(accppta(nx,ny),convppta(nx,ny),sreha(nx,ny),lia(nx,ny),capea(nx,ny))
ALLOCATE(brna(nx,ny),cina(nx,ny),vhela(nx,ny),pwa(nx,ny),cmprefa(nx,ny))
ALLOCATE(refl1kma(nx,ny),refl4kma(nx,ny),u10ma(nx,ny),v10ma(nx,ny))
endif
ALLOCATE(mslpm(nx,ny),hgt250m(nx,ny),hgt500m(nx,ny),hgt850m(nx,ny))
ALLOCATE(temp850m(nx,ny),temp2mm(nx,ny),dewp2mm(nx,ny),accpptm(nx,ny),convpptm(nx,ny))
ALLOCATE(cmprefm(nx,ny),refl1kmm(nx,ny),refl4kmm(nx,ny),u10mm(nx,ny),v10mm(nx,ny))
ALLOCATE(mslpt(nx,ny),hgt250t(nx,ny),hgt500t(nx,ny),hgt850t(nx,ny))
ALLOCATE(temp850t(nx,ny),temp2mt(nx,ny),dewp2mt(nx,ny),accpptt(nx,ny),convpptt(nx,ny))
ALLOCATE(cmpreft(nx,ny),refl1kmt(nx,ny),refl4kmt(nx,ny),u10mt(nx,ny),v10mt(nx,ny))
ALLOCATE(pbrain(nx,ny,4),pbref1(nx,ny,3),pbref4(nx,ny,3),pbcref(nx,ny,3))
factor = 100./float(nmember) ! probability factor (weight) from each member
! factor = 100./float(nmember)+0.01 ! probability factor (weight) from each member
sp = 0.0
mse = 0.0
sde = 0.0
mnbias = 0.0
corr = 0.0
disp = 0.0
sd1 = 0.0
sd2 = 0.0
sum1 = 0.0
sum2 = 0.0
totlegrid = float(nx*ny)
do nt=1,nf ! begin time loop
! zero following arrays
mslpm=0
hgt250m=0
hgt500m=0
hgt850m=0
temp850m=0
temp2mm=0
dewp2mm=0
accpptm=0
convpptm=0
cmprefm=0
refl1kmm=0
refl4kmm=0
u10mm=0
v10mm=0
mslpt=0
hgt250t=0
hgt500t=0
hgt850t=0
temp850t=0
temp2mt=0
dewp2mt=0
accpptt=0
convpptt=0
cmpreft=0
refl1kmt=0
refl4kmt=0
u10mt=0
v10mt=0
pbrain=0.0
pbref1=0.0
pbref4=0.0
pbcref=0.0
if(anaopt == 1) then ! read in analysis data
filenm = trim(ensdir)//trim(anaheader)//'.mspres'//tlevel(nt)
print *, trim(filenm)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,mslpa)
print *, varname,unitnm,nx,ny
filenm = trim(ensdir)//trim(anaheader)//'.hgt250'//tlevel(nt)
print *, trim(filenm)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,hgt250a)
print *, varname,unitnm,nx,ny
filenm = trim(ensdir)//trim(anaheader)//'.u250__'//tlevel(nt)
print *, trim(filenm)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,uwind250a)
print *, varname,unitnm,nx,ny
filenm = trim(ensdir)//trim(anaheader)//'.v250__'//tlevel(nt)
print *, trim(filenm)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,vwind250a)
print *, varname,unitnm,nx,ny
filenm = trim(ensdir)//trim(anaheader)//'.hgt500'//tlevel(nt)
print *, trim(filenm)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,hgt500a)
print *, varname,unitnm,nx,ny
filenm = trim(ensdir)//trim(anaheader)//'.u500__'//tlevel(nt)
print *, trim(filenm)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,uwind500a)
print *, varname,unitnm,nx,ny
filenm = trim(ensdir)//trim(anaheader)//'.v500__'//tlevel(nt)
print *, trim(filenm)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,vwind500a)
print *, varname,unitnm,nx,ny
filenm = trim(ensdir)//trim(anaheader)//'.hgt850'//tlevel(nt)
print *, trim(filenm)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,hgt850a)
print *, varname,unitnm,nx,ny
filenm = trim(ensdir)//trim(anaheader)//'.u850__'//tlevel(nt)
print *, trim(filenm)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,uwind850a)
print *, varname,unitnm,nx,ny
filenm = trim(ensdir)//trim(anaheader)//'.v850__'//tlevel(nt)
print *, trim(filenm)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,vwind850a)
print *, varname,unitnm,nx,ny
filenm = trim(ensdir)//trim(anaheader)//'.tmp850'//tlevel(nt)
print *, trim(filenm)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,temp850a)
print *, varname,unitnm,nx,ny
filenm = trim(ensdir)//trim(anaheader)//'.temp2m'//tlevel(nt)
print *, trim(filenm)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,temp2ma)
print *, varname,unitnm,nx,ny
filenm = trim(ensdir)//trim(anaheader)//'.dewp2m'//tlevel(nt)
print *, trim(filenm)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,dewp2ma)
print *, varname,unitnm,nx,ny
filenm = trim(ensdir)//trim(anaheader)//'.accppt'//tlevel(nt)
print *, trim(filenm)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,accppta)
print *, varname,unitnm,nx,ny
filenm = trim(ensdir)//trim(anaheader)//'.srh03_'//tlevel(nt)
print *, trim(filenm)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,sreha)
print *, varname,unitnm,nx,ny
filenm = trim(ensdir)//trim(anaheader)//'.vhel__'//tlevel(nt)
print *, trim(filenm)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,vhela)
print *, varname,unitnm,nx,ny
filenm = trim(ensdir)//trim(anaheader)//'.sbcape'//tlevel(nt)
print *, trim(filenm)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,capea)
print *, varname,unitnm,nx,ny
filenm = trim(ensdir)//trim(anaheader)//'.sbcins'//tlevel(nt)
print *, trim(filenm)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,cina)
print *, varname,unitnm,nx,ny
filenm = trim(ensdir)//trim(anaheader)//'.pwat__'//tlevel(nt)
print *, trim(filenm)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,pwa)
print *, varname,unitnm,nx,ny
filenm = trim(ensdir)//trim(anaheader)//'.cmpref'//tlevel(nt)
print *, trim(filenm)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,cmprefa)
print *, varname,unitnm,nx,ny
filenm = trim(ensdir)//trim(anaheader)//'.ref1km'//tlevel(nt)
print *, trim(filenm)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,refl1kma)
print *, varname,unitnm,nx,ny
filenm = trim(ensdir)//trim(anaheader)//'.u10m__'//tlevel(nt)
print *, trim(filenm)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,u10ma)
print *, varname,unitnm,nx,ny
filenm = trim(ensdir)//trim(anaheader)//'.v10m__'//tlevel(nt)
print *, trim(filenm)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,v10ma)
print *, varname,unitnm,nx,ny
endif ! end read in analysis data
do nm=1,nmember ! begin member loop
filenm = trim(ensdir)//trim(memberheader(nm))//'.mspres'//tlevel(nt)
print *, trim(filenm)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,mslp)
print *, varname,unitnm,nx,ny
do i=1,nx
do j=1,ny
mslpm(i,j)=mslpm(i,j)+mslp(i,j)
mslpt(i,j)=mslpt(i,j)+mslp(i,j)*mslp(i,j)
enddo
enddo
filenm = trim(ensdir)//trim(memberheader(nm))//'.hgt250'//tlevel(nt)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,hgt250)
print *, varname,unitnm,nx,ny
do i=1,nx
do j=1,ny
hgt250m(i,j)=hgt250m(i,j)+hgt250(i,j)
hgt250t(i,j)=hgt250t(i,j)+hgt250(i,j)*hgt250(i,j)
enddo
enddo
!
filenm = trim(ensdir)//trim(memberheader(nm))//'.u250__'//tlevel(nt)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,uwind250)
print *, varname,unitnm,nx,ny
!
filenm = trim(ensdir)//trim(memberheader(nm))//'.v250__'//tlevel(nt)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,vwind250)
print *, varname,unitnm,nx,ny
filenm = trim(ensdir)//trim(memberheader(nm))//'.hgt500'//tlevel(nt)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,hgt500)
print *, varname,unitnm,nx,ny
do i=1,nx
do j=1,ny
hgt500m(i,j)=hgt500m(i,j)+hgt500(i,j)
hgt500t(i,j)=hgt500t(i,j)+hgt500(i,j)*hgt500(i,j)
enddo
enddo
filenm = trim(ensdir)//trim(memberheader(nm))//'.u500__'//tlevel(nt)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,uwind500)
print *, varname,unitnm,nx,ny
filenm = trim(ensdir)//trim(memberheader(nm))//'.v500__'//tlevel(nt)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,vwind500)
print *, varname,unitnm,nx,ny
filenm = trim(ensdir)//trim(memberheader(nm))//'.hgt850'//tlevel(nt)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,hgt850)
print *, varname,unitnm,nx,ny
!CALL a3dmax0(hgt850,1,nx,1,nx-1,1,ny,1,ny,1,1,1,1,amax,amin)
!print *,'hgt850 -- amax, amin: ',amax, amin
totlegrid1 = 0.0
do i=1,nx
do j=1,ny
IF(hgt850(i,j) > -500.0) THEN ! excluding -9999.0 grids
hgt850m(i,j)=hgt850m(i,j)+hgt850(i,j)
hgt850t(i,j)=hgt850t(i,j)+hgt850(i,j)*hgt850(i,j)
totlegrid1 = totlegrid1 + 1.0
END IF
enddo
enddo
filenm = trim(ensdir)//trim(memberheader(nm))//'.u850__'//tlevel(nt)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,uwind850)
print *, varname,unitnm,nx,ny
filenm = trim(ensdir)//trim(memberheader(nm))//'.v850__'//tlevel(nt)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,vwind850)
print *, varname,unitnm,nx,ny
filenm = trim(ensdir)//trim(memberheader(nm))//'.tmp850'//tlevel(nt)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,temp850)
print *, varname,unitnm,nx,ny
!CALL a3dmax0(temp850,1,nx,1,nx-1,1,ny,1,ny,1,1,1,1,amax,amin)
!print *,'temp850 -- amax, amin: ',amax, amin
totlegrid2 = 0.0
do i=1,nx
do j=1,ny
IF(temp850(i,j) > -500.0) THEN ! excluding -9999.0 grids
temp850m(i,j)=temp850m(i,j)+temp850(i,j)
temp850t(i,j)=temp850t(i,j)+temp850(i,j)*temp850(i,j)
totlegrid2 = totlegrid2 + 1.0
END IF
enddo
enddo
filenm = trim(ensdir)//trim(memberheader(nm))//'.temp2m'//tlevel(nt)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,temp2m)
print *, varname,unitnm,nx,ny
do i=1,nx
do j=1,ny
temp2mm(i,j)=temp2mm(i,j)+temp2m(i,j)
temp2mt(i,j)=temp2mt(i,j)+temp2m(i,j)*temp2m(i,j)
enddo
enddo
filenm = trim(ensdir)//trim(memberheader(nm))//'.dewp2m'//tlevel(nt)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,dewp2m)
print *, varname,unitnm,nx,ny
do i=1,nx
do j=1,ny
dewp2mm(i,j)=dewp2mm(i,j)+dewp2m(i,j)
dewp2mt(i,j)=dewp2mt(i,j)+dewp2m(i,j)*dewp2m(i,j)
enddo
enddo
filenm = trim(ensdir)//trim(memberheader(nm))//'.accppt'//tlevel(nt)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,accppt)
print *, varname,unitnm,nx,ny
do i=1,nx
do j=1,ny
accpptm(i,j)=accpptm(i,j)+accppt(i,j)
accpptt(i,j)=accpptt(i,j)+accppt(i,j)*accppt(i,j)
enddo
enddo
! calculate probability of accppt
do i=1,nx
do j=1,ny
if(accppt(i,j)/25.4.ge.0.10) pbrain(i,j,1)=pbrain(i,j,1)+factor ! >0.10 inch
if(accppt(i,j)/25.4.ge.0.25) pbrain(i,j,2)=pbrain(i,j,2)+factor ! >0.25 inch
if(accppt(i,j)/25.4.ge.0.50) pbrain(i,j,3)=pbrain(i,j,3)+factor ! >0.50 inch
if(accppt(i,j)/25.4.ge.1.00) pbrain(i,j,4)=pbrain(i,j,4)+factor ! >1.00 inch
enddo
enddo
! 5 more fields to be read in
! print *, '5 more fields to be read in.'
filenm = trim(ensdir)//trim(memberheader(nm))//'.srh03_'//tlevel(nt)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,sreh)
print *, varname,unitnm,nx,ny
filenm = trim(ensdir)//trim(memberheader(nm))//'.vhel__'//tlevel(nt)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,vhel)
print *, varname,unitnm,nx,ny
filenm = trim(ensdir)//trim(memberheader(nm))//'.sbcape'//tlevel(nt)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,cape)
print *, varname,unitnm,nx,ny
filenm = trim(ensdir)//trim(memberheader(nm))//'.sbcins'//tlevel(nt)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,cin)
print *, varname,unitnm,nx,ny
filenm = trim(ensdir)//trim(memberheader(nm))//'.pwat__'//tlevel(nt)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,pw)
print *, varname,unitnm,nx,ny
filenm = trim(ensdir)//trim(memberheader(nm))//'.cmpref'//tlevel(nt)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,cmpref)
print *, varname,unitnm,nx,ny
do i=1,nx
do j=1,ny
cmprefm(i,j)=cmprefm(i,j)+cmpref(i,j)
cmpreft(i,j)=cmpreft(i,j)+cmpref(i,j)*cmpref(i,j)
enddo
enddo
! calculate probability of cmpref
do i=1,nx
do j=1,ny
if(cmpref(i,j).ge.35.0) pbcref(i,j,1)=pbcref(i,j,1)+factor ! >35 dBZ
if(cmpref(i,j).ge.45.0) pbcref(i,j,2)=pbcref(i,j,2)+factor ! >45 dBZ
if(cmpref(i,j).ge.55.0) pbcref(i,j,3)=pbcref(i,j,3)+factor ! >55 dBZ
enddo
enddo
filenm = trim(ensdir)//trim(memberheader(nm))//'.ref1km'//tlevel(nt)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,refl1km)
print *, varname,unitnm,nx,ny
do i=1,nx
do j=1,ny
refl1kmm(i,j)=refl1kmm(i,j)+refl1km(i,j)
refl1kmt(i,j)=refl1kmt(i,j)+refl1km(i,j)*refl1km(i,j)
enddo
enddo
! calculate probability of refl
do i=1,nx
do j=1,ny
if(refl1km(i,j).ge.35.0) pbref1(i,j,1)=pbref1(i,j,1)+factor ! >35 dBZ
if(refl1km(i,j).ge.45.0) pbref1(i,j,2)=pbref1(i,j,2)+factor ! >45 dBZ
if(refl1km(i,j).ge.55.0) pbref1(i,j,3)=pbref1(i,j,3)+factor ! >55 dBZ
enddo
enddo
filenm = trim(ensdir)//trim(memberheader(nm))//'.u10m__'//tlevel(nt)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,u10m)
print *, varname,unitnm,nx,ny
do i=1,nx
do j=1,ny
u10mm(i,j)=u10mm(i,j)+u10m(i,j)
u10mt(i,j)=u10mt(i,j)+u10m(i,j)*u10m(i,j)
enddo
enddo
filenm = trim(ensdir)//trim(memberheader(nm))//'.v10m__'//tlevel(nt)
CALL binread2d(nx,ny,trim(filenm),varname,unitnm,v10m)
print *, varname,unitnm,nx,ny
do i=1,nx
do j=1,ny
v10mm(i,j)=v10mm(i,j)+v10m(i,j)
v10mt(i,j)=v10mt(i,j)+v10m(i,j)*v10m(i,j)
enddo
enddo
! if(nt == 1 .AND. nm ==1) then
! do j=1,ny
! write(15,'(5e16.8)') (temp850(i,j),i=1,nx)
! enddo
! endif
enddo ! end member loop
! calculate ensemble mean & square spread
do i=1,nx
do j=1,ny
mslpm(i,j)=mslpm(i,j)/float(nmember)
mslpt(i,j)=mslpt(i,j)/float(nmember)-mslpm(i,j)*mslpm(i,j)
hgt250m(i,j)=hgt250m(i,j)/float(nmember)
hgt250t(i,j)=hgt250t(i,j)/float(nmember)-hgt250m(i,j)*hgt250m(i,j)
hgt500m(i,j)=hgt500m(i,j)/float(nmember)
hgt500t(i,j)=hgt500t(i,j)/float(nmember)-hgt500m(i,j)*hgt500m(i,j)
hgt850m(i,j)=hgt850m(i,j)/float(nmember)
hgt850t(i,j)=hgt850t(i,j)/float(nmember)-hgt850m(i,j)*hgt850m(i,j)
temp850m(i,j)=temp850m(i,j)/float(nmember)
temp850t(i,j)=temp850t(i,j)/float(nmember)-temp850m(i,j)*temp850m(i,j)
temp2mm(i,j)=temp2mm(i,j)/float(nmember)
temp2mt(i,j)=temp2mt(i,j)/float(nmember)-temp2mm(i,j)*temp2mm(i,j)
dewp2mm(i,j)=dewp2mm(i,j)/float(nmember)
dewp2mt(i,j)=dewp2mt(i,j)/float(nmember)-dewp2mm(i,j)*dewp2mm(i,j)
accpptm(i,j)=accpptm(i,j)/float(nmember)
accpptt(i,j)=accpptt(i,j)/float(nmember)-accpptm(i,j)*accpptm(i,j)
! convpptm(i,j)=convpptm(i,j)/float(nmember)
! convpptt(i,j)=convpptt(i,j)/float(nmember)-convpptm(i,j)*convpptm(i,j)
cmprefm(i,j)=cmprefm(i,j)/float(nmember)
cmpreft(i,j)=cmpreft(i,j)/float(nmember)-cmprefm(i,j)*cmprefm(i,j)
refl1kmm(i,j)=refl1kmm(i,j)/float(nmember)
refl1kmt(i,j)=refl1kmt(i,j)/float(nmember)-refl1kmm(i,j)*refl1kmm(i,j)
u10mm(i,j)=u10mm(i,j)/float(nmember)
u10mt(i,j)=u10mt(i,j)/float(nmember)-u10mm(i,j)*u10mm(i,j)
v10mm(i,j)=v10mm(i,j)/float(nmember)
v10mt(i,j)=v10mt(i,j)/float(nmember)-v10mm(i,j)*v10mm(i,j)
enddo
enddo
do i=1,nx
do j=1,ny
mslpt(i,j)=abs(mslpt(i,j))
hgt250t(i,j)=abs(hgt250t(i,j))
hgt500t(i,j)=abs(hgt500t(i,j))
hgt850t(i,j)=abs(hgt850t(i,j))
temp850t(i,j)=abs(temp850t(i,j))
temp2mt(i,j)=abs(temp2mt(i,j))
dewp2mt(i,j)=abs(dewp2mt(i,j))
accpptt(i,j)=abs(accpptt(i,j))
! convpptt(i,j)=abs(convpptt(i,j))
cmpreft(i,j)=abs(cmpreft(i,j))
refl1kmt(i,j)=abs(refl1kmt(i,j))
u10mt(i,j)=abs(u10mt(i,j))
v10mt(i,j)=abs(v10mt(i,j))
enddo
enddo
! write out ensemble means
outfile = trim(outdir)//'mean/'//trim(outheader)//'.'
filenm=trim(outfile)//'mslp__'//tlevel(nt)
print *, trim(filenm)
q3dname='MSLP'
q3dunit='hPa'
CALL bindump2d(nx,ny,trim(filenm),q3dname,q3dunit,mslpm)
filenm=trim(outfile)//'hgt250'//tlevel(nt)
print *, trim(filenm)
q3dname='HGHT'
q3dunit='m'
CALL bindump2d(nx,ny,trim(filenm),q3dname,q3dunit,hgt250m)
filenm=trim(outfile)//'hgt500'//tlevel(nt)
print *, trim(filenm)
q3dname='HGHT'
q3dunit='m'
CALL bindump2d(nx,ny,trim(filenm),q3dname,q3dunit,hgt500m)
filenm=trim(outfile)//'hgt850'//tlevel(nt)
print *, trim(filenm)
q3dname='HGHT'
q3dunit='m'
CALL bindump2d(nx,ny,trim(filenm),q3dname,q3dunit,hgt850m)
filenm=trim(outfile)//'tmp850'//tlevel(nt)
print *, trim(filenm)
q3dname='TMPC'
q3dunit='C'
CALL bindump2d(nx,ny,trim(filenm),q3dname,q3dunit,temp850m)
filenm=trim(outfile)//'temp2m'//tlevel(nt)
print *, trim(filenm)
q3dname='temp2m'
q3dunit='F'
CALL bindump2d(nx,ny,trim(filenm),q3dname,q3dunit,temp2mm)
filenm=trim(outfile)//'dewp2m'//tlevel(nt)
print *, trim(filenm)
q3dname='dewp2m'
q3dunit='F'
CALL bindump2d(nx,ny,trim(filenm),q3dname,q3dunit,dewp2mm)
filenm=trim(outfile)//'accppt'//tlevel(nt)
print *, trim(filenm)
q3dname='RAIN'
q3dunit='mm'
CALL bindump2d(nx,ny,trim(filenm),q3dname,q3dunit,accpptm)
filenm=trim(outfile)//'cmpref'//tlevel(nt)
print *, trim(filenm)
q3dname='cmpref'
q3dunit='dBZ'
CALL bindump2d(nx,ny,trim(filenm),q3dname,q3dunit,cmprefm)
filenm=trim(outfile)//'ref1km'//tlevel(nt)
print *, trim(filenm)
q3dname='refl1km'
q3dunit='dBZ'
CALL bindump2d(nx,ny,trim(filenm),q3dname,q3dunit,refl1kmm)
filenm=trim(outfile)//'u10m__'//tlevel(nt)
print *, trim(filenm)
q3dname='UREL'
q3dunit='m/s'
CALL bindump2d(nx,ny,trim(filenm),q3dname,q3dunit,u10mm)
filenm=trim(outfile)//'v10m__'//tlevel(nt)
print *, trim(filenm)
q3dname='VREL'
q3dunit='m/s'
CALL bindump2d(nx,ny,trim(filenm),q3dname,q3dunit,v10mm)
! write out probabilities
outfile = trim(outdir)//'prob/'//trim(outheader)//'.'
filenm=trim(outfile)//'pbr010'//tlevel(nt)
print *, trim(filenm)
q3dname='PBRAIN010'
q3dunit='%'
CALL bindump2d(nx,ny,trim(filenm),q3dname,q3dunit,pbrain(1,1,1))
filenm=trim(outfile)//'pbr025'//tlevel(nt)
print *, trim(filenm)
q3dname='PBRAIN025'
q3dunit='%'
CALL bindump2d(nx,ny,trim(filenm),q3dname,q3dunit,pbrain(1,1,2))
filenm=trim(outfile)//'pbr050'//tlevel(nt)
print *, trim(filenm)
q3dname='PBRAIN050'
q3dunit='%'
CALL bindump2d(nx,ny,trim(filenm),q3dname,q3dunit,pbrain(1,1,3))
filenm=trim(outfile)//'pbr100'//tlevel(nt)
print *, trim(filenm)
q3dname='PBRAIN100'
q3dunit='%'
CALL bindump2d(nx,ny,trim(filenm),q3dname,q3dunit,pbrain(1,1,4))
filenm=trim(outfile)//'pbz135'//tlevel(nt)
print *, trim(filenm)
q3dname='PBREFLC1KM35'
q3dunit='%'
CALL bindump2d(nx,ny,trim(filenm),q3dname,q3dunit,pbref1(1,1,1))
filenm=trim(outfile)//'pbz145'//tlevel(nt)
print *, trim(filenm)
q3dname='PBREFLC1KM45'
q3dunit='%'
CALL bindump2d(nx,ny,trim(filenm),q3dname,q3dunit,pbref1(1,1,2))
filenm=trim(outfile)//'pbz155'//tlevel(nt)
print *, trim(filenm)
q3dname='PBREFLC1KM55'
q3dunit='%'
CALL bindump2d(nx,ny,trim(filenm),q3dname,q3dunit,pbref1(1,1,3))
filenm=trim(outfile)//'pbcz35'//tlevel(nt)
print *, trim(filenm)
q3dname='PBCMPREF35'
q3dunit='%'
CALL bindump2d(nx,ny,trim(filenm),q3dname,q3dunit,pbcref(1,1,1))
filenm=trim(outfile)//'pbcz45'//tlevel(nt)
print *, trim(filenm)
q3dname='PBCMPREF45'
q3dunit='%'
CALL bindump2d(nx,ny,trim(filenm),q3dname,q3dunit,pbcref(1,1,2))
filenm=trim(outfile)//'pbcz55'//tlevel(nt)
print *, trim(filenm)
q3dname='PBCMPREF55'
q3dunit='%'
CALL bindump2d(nx,ny,trim(filenm),q3dname,q3dunit,pbcref(1,1,3))
! calculate domain rms spread
do i=1,nx
do j=1,ny
sp(nt,1)=sp(nt,1)+mslpt(i,j)
sp(nt,2)=sp(nt,2)+hgt250t(i,j)
sp(nt,3)=sp(nt,3)+hgt500t(i,j)
sp(nt,4)=sp(nt,4)+hgt850t(i,j)
sp(nt,5)=sp(nt,5)+temp850t(i,j)
sp(nt,6)=sp(nt,6)+temp2mt(i,j)
sp(nt,7)=sp(nt,7)+dewp2mt(i,j)
sp(nt,8)=sp(nt,8)+accpptt(i,j)
! sp(nt,9)=sp(nt,9)+convpptt(i,j)
sp(nt,10)=sp(nt,10)+cmpreft(i,j)
sp(nt,11)=sp(nt,11)+refl1kmt(i,j)
sp(nt,12)=sp(nt,12)+u10mt(i,j)
sp(nt,13)=sp(nt,13)+v10mt(i,j)
enddo
enddo
print *,'totlegrid,totlegrid1,totlegrid2:', &
int(totlegrid),int(totlegrid1),int(totlegrid2)
sp(nt,1)=sqrt(sp(nt,1)/totlegrid) ! domain average spread for mslp
sp(nt,2)=sqrt(sp(nt,2)/totlegrid) ! domain average spread for hgt250
sp(nt,3)=sqrt(sp(nt,3)/totlegrid) ! domain average spread for hgt500
sp(nt,4)=sqrt(sp(nt,4)/totlegrid1) ! domain average spread for hgt850
sp(nt,5)=sqrt(sp(nt,5)/totlegrid2) ! domain average spread for temp850
sp(nt,6)=sqrt(sp(nt,6)/totlegrid) ! domain average spread for temp2m
sp(nt,7)=sqrt(sp(nt,7)/totlegrid) ! domain average spread for dewp2m
sp(nt,8)=sqrt(sp(nt,8)/totlegrid) ! domain average spread for accppt
! sp(nt,9)=sqrt(sp(nt,9)/totlegrid) ! domain average spread for convppt
sp(nt,10)=sqrt(sp(nt,10)/totlegrid) ! domain average spread for cmpref
sp(nt,11)=sqrt(sp(nt,11)/totlegrid) ! domain average spread for refl
sp(nt,12)=sqrt(sp(nt,12)/totlegrid) ! domain average spread for u10m
sp(nt,13)=sqrt(sp(nt,13)/totlegrid) ! domain average spread for v10m
! calculate and write out ensemble spread
do i=1,nx
do j=1,ny
mslpt(i,j)=sqrt(mslpt(i,j))
hgt250t(i,j)=sqrt(hgt250t(i,j))
hgt500t(i,j)=sqrt(hgt500t(i,j))
hgt850t(i,j)=sqrt(hgt850t(i,j))
temp850t(i,j)=sqrt(temp850t(i,j))
temp2mt(i,j)=sqrt(temp2mt(i,j))
dewp2mt(i,j)=sqrt(dewp2mt(i,j))
accpptt(i,j)=sqrt(accpptt(i,j))
! convpptt(i,j)=sqrt(convpptt(i,j))
cmpreft(i,j)=sqrt(cmpreft(i,j))
refl1kmt(i,j)=sqrt(refl1kmt(i,j))
u10mt(i,j)=sqrt(u10mt(i,j))
v10mt(i,j)=sqrt(v10mt(i,j))
enddo
enddo
! write out binary format spread fields
outfile = trim(outdir)//'sprd/'//trim(outheader)//'.'
filenm=trim(outfile)//'mslp__'//tlevel(nt)
print *, trim(filenm)
q3dname='MSLP'
q3dunit='hPa'
CALL bindump2d(nx,ny,trim(filenm),q3dname,q3dunit,mslpt)
filenm=trim(outfile)//'hgt500'//tlevel(nt)
print *, trim(filenm)
q3dname='HGHT'
q3dunit='m'
CALL bindump2d(nx,ny,trim(filenm),q3dname,q3dunit,hgt500t)
filenm=trim(outfile)//'temp2m'//tlevel(nt)
print *, trim(filenm)
q3dname='temp2m'
q3dunit='F'
CALL bindump2d(nx,ny,trim(filenm),q3dname,q3dunit,temp2mt)
filenm=trim(outfile)//'dewp2m'//tlevel(nt)
print *, trim(filenm)
q3dname='dewp2m'
q3dunit='F'
CALL bindump2d(nx,ny,trim(filenm),q3dname,q3dunit,dewp2mt)
filenm=trim(outfile)//'accppt'//tlevel(nt)
print *, trim(filenm)
q3dname='RAIN'
q3dunit='mm'
CALL bindump2d(nx,ny,trim(filenm),q3dname,q3dunit,accppt)
filenm=trim(outfile)//'ref1km'//tlevel(nt)
print *, trim(filenm)
q3dname='refl1km'
q3dunit='dBZ'
CALL bindump2d(nx,ny,trim(filenm),q3dname,q3dunit,refl1kmt)
filenm=trim(outfile)//'u10m__'//tlevel(nt)
print *, trim(filenm)
q3dname='u10m'
q3dunit='m/s'
CALL bindump2d(nx,ny,trim(filenm),q3dname,q3dunit,u10mt)
filenm=trim(outfile)//'v10m__'//tlevel(nt)
print *, trim(filenm)
q3dname='v10m'
q3dunit='m/s'
CALL bindump2d(nx,ny,trim(filenm),q3dname,q3dunit,v10mt)
if(anaopt ==1) then ! calculate error-related quantities
! calculate mse, ...
do i=1,nx
do j=1,ny
sum1(nt,1)=sum1(nt,1)+mslpm(i,j)
sum1(nt,3)=sum1(nt,3)+hgt500m(i,j)
sum1(nt,5)=sum1(nt,5)+temp850m(i,j)
sum1(nt,6)=sum1(nt,6)+temp2mm(i,j)
sum1(nt,7)=sum1(nt,7)+dewp2mm(i,j)
sum1(nt,8)=sum1(nt,8)+accpptm(i,j)
sum1(nt,11)=sum1(nt,11)+refl1kmm(i,j)
sum1(nt,12)=sum1(nt,12)+u10mm(i,j)
sum1(nt,13)=sum1(nt,13)+v10mm(i,j)
sum2(nt,1)=sum2(nt,1)+mslpa(i,j)
sum2(nt,3)=sum2(nt,3)+hgt500a(i,j)
sum2(nt,5)=sum2(nt,5)+temp850a(i,j)
sum2(nt,6)=sum2(nt,6)+temp2ma(i,j)
sum2(nt,7)=sum2(nt,7)+dewp2ma(i,j)
sum2(nt,8)=sum2(nt,8)+accppta(i,j)
sum2(nt,11)=sum2(nt,11)+refl1kma(i,j)
sum2(nt,12)=sum2(nt,12)+u10ma(i,j)
sum2(nt,13)=sum2(nt,13)+v10ma(i,j)
mse(nt,1)=mse(nt,1)+(mslpm(i,j)-mslpa(i,j))**2
mse(nt,3)=mse(nt,3)+(hgt500m(i,j)-hgt500a(i,j))**2
mse(nt,5)=mse(nt,5)+(temp850m(i,j)-temp850a(i,j))**2
mse(nt,6)=mse(nt,6)+(temp2mm(i,j)-temp2ma(i,j))**2
mse(nt,7)=mse(nt,7)+(dewp2mm(i,j)-dewp2ma(i,j))**2
mse(nt,8)=mse(nt,8)+(accpptm(i,j)-accppta(i,j))**2
mse(nt,11)=mse(nt,11)+(refl1kmm(i,j)-refl1kma(i,j))**2
mse(nt,12)=mse(nt,12)+(u10mm(i,j)-u10ma(i,j))**2
mse(nt,13)=mse(nt,13)+(v10mm(i,j)-v10ma(i,j))**2
mnbias(nt,1)=mnbias(nt,1)+(mslpm(i,j)-mslpa(i,j))
mnbias(nt,3)=mnbias(nt,3)+(hgt500m(i,j)-hgt500a(i,j))
mnbias(nt,5)=mnbias(nt,5)+(temp850m(i,j)-temp850a(i,j))
mnbias(nt,6)=mnbias(nt,6)+(temp2mm(i,j)-temp2ma(i,j))
mnbias(nt,7)=mnbias(nt,7)+(dewp2mm(i,j)-dewp2ma(i,j))
mnbias(nt,8)=mnbias(nt,8)+(accpptm(i,j)-accppta(i,j))
mnbias(nt,11)=mnbias(nt,11)+(refl1kmm(i,j)-refl1kma(i,j))
mnbias(nt,12)=mnbias(nt,12)+(u10mm(i,j)-u10ma(i,j))
mnbias(nt,13)=mnbias(nt,13)+(v10mm(i,j)-v10ma(i,j))
enddo
enddo
sum1(nt,1)=sum1(nt,1)/totlegrid
sum1(nt,3)=sum1(nt,3)/totlegrid
sum1(nt,5)=sum1(nt,5)/totlegrid2
sum1(nt,6)=sum1(nt,6)/totlegrid
sum1(nt,7)=sum1(nt,7)/totlegrid
sum1(nt,8)=sum1(nt,8)/totlegrid
sum1(nt,11)=sum1(nt,11)/totlegrid
sum1(nt,12)=sum1(nt,12)/totlegrid
sum1(nt,13)=sum1(nt,13)/totlegrid
sum2(nt,1)=sum2(nt,1)/totlegrid
sum2(nt,3)=sum2(nt,3)/totlegrid
sum2(nt,5)=sum2(nt,5)/totlegrid2
sum2(nt,6)=sum2(nt,6)/totlegrid
sum2(nt,7)=sum2(nt,7)/totlegrid
sum2(nt,8)=sum2(nt,8)/totlegrid
sum2(nt,11)=sum2(nt,11)/totlegrid
sum2(nt,12)=sum2(nt,12)/totlegrid
sum2(nt,13)=sum2(nt,13)/totlegrid
mse(nt,1)=sqrt(mse(nt,1)/totlegrid)
mse(nt,3)=sqrt(mse(nt,3)/totlegrid)
mse(nt,5)=sqrt(mse(nt,5)/totlegrid2)
mse(nt,6)=sqrt(mse(nt,6)/totlegrid)
mse(nt,7)=sqrt(mse(nt,7)/totlegrid)
mse(nt,8)=sqrt(mse(nt,8)/totlegrid)
mse(nt,11)=sqrt(mse(nt,11)/totlegrid)
mse(nt,12)=sqrt(mse(nt,12)/totlegrid)
mse(nt,13)=sqrt(mse(nt,13)/totlegrid)
mnbias(nt,1)=mnbias(nt,1)/totlegrid
mnbias(nt,3)=mnbias(nt,3)/totlegrid
mnbias(nt,5)=mnbias(nt,5)/totlegrid2
mnbias(nt,6)=mnbias(nt,6)/totlegrid
mnbias(nt,7)=mnbias(nt,7)/totlegrid
mnbias(nt,8)=mnbias(nt,8)/totlegrid
mnbias(nt,11)=mnbias(nt,11)/totlegrid
mnbias(nt,12)=mnbias(nt,12)/totlegrid
mnbias(nt,13)=mnbias(nt,13)/totlegrid
do i=1,nx
do j=1,ny
sd1(nt,1)=sd1(nt,1)+(mslpm(i,j)-sum1(nt,1))**2
sd1(nt,3)=sd1(nt,3)+(hgt500m(i,j)-sum1(nt,3))**2
sd1(nt,5)=sd1(nt,5)+(temp850m(i,j)-sum1(nt,5))**2
sd1(nt,6)=sd1(nt,6)+(temp2mm(i,j)-sum1(nt,6))**2
sd1(nt,7)=sd1(nt,7)+(dewp2mm(i,j)-sum1(nt,7))**2
sd1(nt,8)=sd1(nt,8)+(accpptm(i,j)-sum1(nt,8))**2
sd1(nt,11)=sd1(nt,11)+(refl1kmm(i,j)-sum1(nt,11))**2
sd1(nt,12)=sd1(nt,12)+(u10mm(i,j)-sum1(nt,12))**2
sd1(nt,13)=sd1(nt,13)+(v10mm(i,j)-sum1(nt,13))**2
sd2(nt,1)=sd2(nt,1)+(mslpa(i,j)-sum2(nt,1))**2
sd2(nt,3)=sd2(nt,3)+(hgt500a(i,j)-sum2(nt,3))**2
sd2(nt,5)=sd2(nt,5)+(temp850a(i,j)-sum2(nt,5))**2
sd2(nt,6)=sd2(nt,6)+(temp2ma(i,j)-sum2(nt,6))**2
sd2(nt,7)=sd2(nt,7)+(dewp2ma(i,j)-sum2(nt,7))**2
sd2(nt,8)=sd2(nt,8)+(accppta(i,j)-sum2(nt,8))**2
sd2(nt,11)=sd2(nt,11)+(refl1kma(i,j)-sum2(nt,11))**2
sd2(nt,12)=sd2(nt,12)+(u10ma(i,j)-sum2(nt,12))**2
sd2(nt,13)=sd2(nt,13)+(v10ma(i,j)-sum2(nt,13))**2
sde(nt,1)=sde(nt,1)+(mslpm(i,j)-sum1(nt,1)-mslpa(i,j)+sum2(nt,1))**2
sde(nt,3)=sde(nt,3)+(hgt500m(i,j)-sum1(nt,3)-hgt500a(i,j)+sum2(nt,3))**2
sde(nt,5)=sde(nt,5)+(temp850m(i,j)-sum1(nt,5)-temp850a(i,j)+sum2(nt,5))**2
sde(nt,6)=sde(nt,6)+(temp2mm(i,j)-sum1(nt,6)-temp2ma(i,j)+sum2(nt,6))**2
sde(nt,7)=sde(nt,7)+(dewp2mm(i,j)-sum1(nt,7)-dewp2ma(i,j)+sum2(nt,7))**2
sde(nt,8)=sde(nt,8)+(accpptm(i,j)-sum1(nt,8)-accppta(i,j)+sum2(nt,8))**2
sde(nt,11)=sde(nt,11)+(refl1kmm(i,j)-sum1(nt,11)-refl1kma(i,j)+ &
sum2(nt,11))**2
sde(nt,12)=sde(nt,12)+(u10mm(i,j)-sum1(nt,12)-u10ma(i,j)+sum2(nt,12))**2
sde(nt,13)=sde(nt,13)+(v10mm(i,j)-sum1(nt,13)-v10ma(i,j)+sum2(nt,13))**2
corr(nt,1)=corr(nt,1)+(mslpm(i,j)-sum1(nt,1))*(mslpa(i,j)-sum2(nt,1))
corr(nt,3)=corr(nt,3)+(hgt500m(i,j)-sum1(nt,3))*(hgt500a(i,j)-sum2(nt,3))
corr(nt,5)=corr(nt,5)+(temp850m(i,j)-sum1(nt,5))*(temp850a(i,j)-sum2(nt,5))
corr(nt,6)=corr(nt,6)+(temp2mm(i,j)-sum1(nt,6))*(temp2ma(i,j)-sum2(nt,6))
corr(nt,7)=corr(nt,7)+(dewp2mm(i,j)-sum1(nt,7))*(dewp2ma(i,j)-sum2(nt,7))
corr(nt,8)=corr(nt,8)+(accpptm(i,j)-sum1(nt,8))*(accppta(i,j)-sum2(nt,8))
corr(nt,11)=corr(nt,11)+(refl1kmm(i,j)-sum1(nt,11))*(refl1kma(i,j)- &
sum2(nt,11))
corr(nt,12)=corr(nt,12)+(u10mm(i,j)-sum1(nt,12))*(u10ma(i,j)-sum2(nt,12))
corr(nt,13)=corr(nt,13)+(v10mm(i,j)-sum1(nt,13))*(v10ma(i,j)-sum2(nt,13))
enddo
enddo
sd1(nt,1)=sqrt(sd1(nt,1)/totlegrid)
sd1(nt,3)=sqrt(sd1(nt,3)/totlegrid)
sd1(nt,5)=sqrt(sd1(nt,5)/totlegrid2)
sd1(nt,6)=sqrt(sd1(nt,6)/totlegrid)
sd1(nt,7)=sqrt(sd1(nt,7)/totlegrid)
sd1(nt,8)=sqrt(sd1(nt,8)/totlegrid)
sd1(nt,11)=sqrt(sd1(nt,11)/totlegrid)
sd1(nt,12)=sqrt(sd1(nt,12)/totlegrid)
sd1(nt,13)=sqrt(sd1(nt,13)/totlegrid)
sd2(nt,1)=sqrt(sd2(nt,1)/totlegrid)
sd2(nt,3)=sqrt(sd2(nt,3)/totlegrid)
sd2(nt,5)=sqrt(sd2(nt,5)/totlegrid2)
sd2(nt,6)=sqrt(sd2(nt,6)/totlegrid)
sd2(nt,7)=sqrt(sd2(nt,7)/totlegrid)
sd2(nt,8)=sqrt(sd2(nt,8)/totlegrid)
sd2(nt,11)=sqrt(sd2(nt,11)/totlegrid)
sd2(nt,12)=sqrt(sd2(nt,12)/totlegrid)
sd2(nt,13)=sqrt(sd2(nt,13)/totlegrid)
sde(nt,1)=sqrt(sde(nt,1)/totlegrid)
sde(nt,3)=sqrt(sde(nt,3)/totlegrid)
sde(nt,5)=sqrt(sde(nt,5)/totlegrid2)
sde(nt,6)=sqrt(sde(nt,6)/totlegrid)
sde(nt,7)=sqrt(sde(nt,7)/totlegrid)
sde(nt,8)=sqrt(sde(nt,8)/totlegrid)
sde(nt,11)=sqrt(sde(nt,11)/totlegrid)
sde(nt,12)=sqrt(sde(nt,12)/totlegrid)
sde(nt,13)=sqrt(sde(nt,13)/totlegrid)
sdbias(nt,1)=sd1(nt,1)-sd2(nt,1)
sdbias(nt,3)=sd1(nt,3)-sd2(nt,3)
sdbias(nt,5)=sd1(nt,5)-sd2(nt,5)
sdbias(nt,6)=sd1(nt,6)-sd2(nt,6)
sdbias(nt,7)=sd1(nt,7)-sd2(nt,7)
sdbias(nt,8)=sd1(nt,8)-sd2(nt,8)
sdbias(nt,11)=sd1(nt,11)-sd2(nt,11)
sdbias(nt,12)=sd1(nt,12)-sd2(nt,12)
sdbias(nt,13)=sd1(nt,13)-sd2(nt,13)
! corr(nt,1)=corr(nt,1)/(totlegrid*sd1(nt,1)*sd2(nt,1))
! corr(nt,3)=corr(nt,3)/(totlegrid*sd1(nt,3)*sd2(nt,3))
! corr(nt,5)=corr(nt,5)/(totlegrid*sd1(nt,5)*sd2(nt,5))
! corr(nt,6)=corr(nt,6)/(totlegrid*sd1(nt,6)*sd2(nt,6))
! corr(nt,7)=corr(nt,7)/(totlegrid*sd1(nt,7)*sd2(nt,7))
do ll=1,13
if(sd1(nt,ll)*sd2(nt,ll)<=0.0) then
corr(nt,ll)=1.0
else
if(ll == 5) then
corr(nt,ll)=amin1(corr(nt,ll)/(totlegrid2*sd1(nt,ll)*sd2(nt,ll)),1.0)
else
corr(nt,ll)=amin1(corr(nt,ll)/(totlegrid*sd1(nt,ll)*sd2(nt,ll)),1.0)
endif
endif
enddo
! corr(nt,11)=corr(nt,11)/(totlegrid*sd1(nt,11)*sd2(nt,11))
! corr(nt,12)=corr(nt,12)/(totlegrid*sd1(nt,12)*sd2(nt,12))
! corr(nt,13)=corr(nt,13)/(totlegrid*sd1(nt,13)*sd2(nt,13))
disp(nt,1)=sqrt(2.*(1.-corr(nt,1))*sd1(nt,1)*sd2(nt,1))
disp(nt,3)=sqrt(2.*(1.-corr(nt,3))*sd1(nt,3)*sd2(nt,3))
disp(nt,5)=sqrt(2.*(1.-corr(nt,5))*sd1(nt,5)*sd2(nt,5))
disp(nt,6)=sqrt(2.*(1.-corr(nt,6))*sd1(nt,6)*sd2(nt,6))
disp(nt,7)=sqrt(2.*(1.-corr(nt,7))*sd1(nt,7)*sd2(nt,7))
disp(nt,8)=sqrt(2.*(1.-corr(nt,8))*sd1(nt,8)*sd2(nt,8))
disp(nt,11)=sqrt(2.*(1.-corr(nt,11))*sd1(nt,11)*sd2(nt,11))
disp(nt,12)=sqrt(2.*(1.-corr(nt,12))*sd1(nt,12)*sd2(nt,12))
disp(nt,13)=sqrt(2.*(1.-corr(nt,13))*sd1(nt,13)*sd2(nt,13))
endif ! end calculating error-related quantities
enddo ! end time loop
! write out spread arrays
open(4, file=trim(outdir)//'ens_sp.asc',form='formatted')
write(4,*) '& ntimes =',nf
write(4,'(a)') '& mslp hgt250 hgt500 hgt850 temp850 temp2m dewp2m accppt convppt cmpref refl u10m v10m'
do j=1,nf
write(4,'(i5,13f10.4)') j-1,(sp(j,i),i=1,13)
enddo
close(4)
if(anaopt ==1) then ! write out ensemble error analysis
open(4, file=trim(outdir)//'ens_err_rmse.asc',form='formatted')
write(4,*) '& ntimes =',nf
write(4,'(a)') '& mslp hgt250 hgt500 hgt850 temp850 temp2m dewp2m accppt convppt cmpref refl u10m v10m'
do j=1,nf
write(4,'(i5,13f10.4)') j-1,(mse(j,i),i=1,13)
enddo
close(4)
open(4, file=trim(outdir)//'ens_err_sde.asc',form='formatted')
write(4,*) '& ntimes =',nf
write(4,'(a)') '& mslp hgt250 hgt500 hgt850 temp850 temp2m dewp2m accppt convppt cmpref refl u10m v10m'
do j=1,nf
write(4,'(i5,13f10.4)') j-1,(sde(j,i),i=1,13)
enddo
close(4)
open(4, file=trim(outdir)//'ens_err_mnbias.asc',form='formatted')
write(4,*) '& ntimes =',nf
write(4,'(a)') '& mslp hgt250 hgt500 hgt850 temp850 temp2m dewp2m accppt convppt cmpref refl u10m v10m'
do j=1,nf
write(4,'(i5,13f10.4)') j-1,(mnbias(j,i),i=1,13)
enddo
close(4)
open(4, file=trim(outdir)//'ens_err_sdbias.asc',form='formatted')
write(4,*) '& ntimes =',nf
write(4,'(a)') '& mslp hgt250 hgt500 hgt850 temp850 temp2m dewp2m accppt convppt cmpref refl u10m v10m'
do j=1,nf
write(4,'(i5,13f10.4)') j-1,(sdbias(j,i),i=1,13)
enddo
close(4)
open(4, file=trim(outdir)//'ens_err_corr.asc',form='formatted')
write(4,*) '& ntimes =',nf
write(4,'(a)') '& mslp hgt250 hgt500 hgt850 temp850 temp2m dewp2m accppt convppt cmpref refl u10m v10m'
do j=1,nf
write(4,'(i5,13f10.4)') j-1,(100.*corr(j,i),i=1,13)
enddo
close(4)
open(4, file=trim(outdir)//'ens_err_disp.asc',form='formatted')
write(4,*) '& ntimes =',nf
write(4,'(a)') '& mslp hgt250 hgt500 hgt850 temp850 temp2m dewp2m accppt convppt cmpref refl u10m v10m'
do j=1,nf
write(4,'(i5,13f10.4)') j-1,(disp(j,i),i=1,13)
enddo
close(4)
endif
STOP
9000 CONTINUE ! I/O errors
CLOSE(UNIT=nunit)
CALL retunit(nunit)
IF (istat < 0) THEN
iSTATUS = 2
PRINT *, 'BINREAD2D: I/O ERRORS OCCURRED ', &
'(possible end of record or file): ',istat, iSTATUS
ELSE IF (istat > 0) THEN
iSTATUS = 2
PRINT *, 'BINREAD2D: UNRECOVERABLE I/O ERRORS OCCURRED: ', &
istat,iSTATUS
END IF
END PROGRAM ens_ana