program calibrate2_svr
c
c Just like calibrate2.f, except reads the output of the combined sevr probability table. DRB 12/10/2004
parameter (igemy=129, igemx=185) ! GEMPAK grid size
c
real prct(igemy,igemx,2),prob,pa(200),pb(200),rk,bias(200),
& grd(igemx,igemy),dum1,dum2
integer igemy,igemx,itcnt,jstop,jdel
character*150 pctfila,pctfilb
character dummy*20,date*8,fhour*4,frun*2,dateetc*72
call getarg(1,pctfila) ! forecast file
call getarg(2,pctfilb) ! forecast file2
call getarg(3,date) ! YYYYMMDD
call getarg(4,frun) ! 09 or 21
call getarg(5,fhour) ! f000 ... f063...
do i=1,72
dateetc(i:i) = ' '
enddo
dateetc(1:6) = date(3:8)
dateetc(7:7) = '/'
dateetc(8:9) = frun(1:2)
dateetc(10:12) = '00F'
dateetc(13:15) = fhour(2:4)
dateetc(40:40) = '0'
dateetc(51:54) = 'NONE'
dateetc(56:61) = 'CALSVR'
c
c fill all possible forecast values...
itcnt = 0
do i=0,100,10
do j = 0,100,10
itcnt = itcnt + 1
pa(itcnt) = float(i) ! fcst1
pb(itcnt) = float(j) ! fcst2
enddo
enddo
write(*,*) 'Number of possible combinations... ',itcnt
1 format(a)
c Now read them in...
c
c 1. Read in the percent file:
open(unit=10,file=pctfila,status='old',form='formatted')
write(*,*) 'Reading the GEMPAK forecast file 1 now: ',pctfila
7775 read(10,1) dummy
if(dummy(2:4).eq.'ROW') then
backspace(unit=10)
do j=igemy,1,-1
if(j.ge.100) then
read(10,*) dummy(1:8),(prct(j,i,1),i=1,igemx)
else
read(10,*) dummy(1:8),irow,(prct(j,i,1),i=1,igemx)
endif
enddo
else
goto 7775
endif
close(unit=10)
c 2. Read in the percent2 file:
open(unit=10,file=pctfilb,status='old',form='formatted')
write(*,*) 'Reading the GEMPAK forecast file 2 now: ',pctfilb
1775 read(10,1) dummy
if(dummy(2:4).eq.'ROW') then
backspace(unit=10)
do j=igemy,1,-1
if(j.ge.100) then
read(10,*) dummy(1:8),(prct(j,i,2),i=1,igemx)
else
read(10,*) dummy(1:8),irow,(prct(j,i,2),i=1,igemx)
endif
enddo
else
goto 1775
endif
close(unit=10)
c
c Now the grids are in...
c Do a little work here. First, bin probs to 10% intervals. Then, change verification
c to 1's or 0's.
do i=1,igemx
do j=1,igemy
do k=1,2
if(prct(j,i,k).lt.-1.0) prct(j,i,k) = -9999.0
if(prct(j,i,k).ge.-1.0.and.prct(j,i,k).lt.5.0)prct(j,i,k)= 0.0
if(prct(j,i,k).ge.5.0.and.prct(j,i,k).lt.15.) prct(j,i,k)=10.0
if(prct(j,i,k).ge.15.0.and.prct(j,i,k).lt.25.)prct(j,i,k)=20.0
if(prct(j,i,k).ge.25.0.and.prct(j,i,k).lt.35.)prct(j,i,k)=30.0
if(prct(j,i,k).ge.35.0.and.prct(j,i,k).lt.45.)prct(j,i,k)=40.0
if(prct(j,i,k).ge.45.0.and.prct(j,i,k).lt.55.)prct(j,i,k)=50.0
if(prct(j,i,k).ge.55.0.and.prct(j,i,k).lt.65.)prct(j,i,k)=60.0
if(prct(j,i,k).ge.65.0.and.prct(j,i,k).lt.75.)prct(j,i,k)=70.0
if(prct(j,i,k).ge.75.0.and.prct(j,i,k).lt.85.)prct(j,i,k)=80.0
if(prct(j,i,k).ge.85.0.and.prct(j,i,k).lt.95.)prct(j,i,k)=90.0
if(prct(j,i,k).ge.95.0) prct(j,i,k)=100.0
enddo
enddo
enddo
c
c Read the file that contains the dependency coefficient and the bias...
c
open(unit=11,file='combine_probs_svr1.out',status='old')
do ii=1,itcnt
read(11,*) pa(ii),pb(ii)
read(11,*) dum1
read(11,*) dum1,dum2
read(11,*) dum1,dum2
if((dum1+dum2).lt.0.001e-10) then
bias(ii) = 0.0
else
bias(ii) = (dum2/(dum1+dum2))*100.
endif
598 format(2f10.2)
600 format(f10.2)
601 format(8f10.2)
enddo
close (11)
c Now...bin the forecasts into all psbl combinations. Record if a hit occurred
c during this combined period too.
c
do i=1,igemx
do j=1,igemy
ii = 0
do k=1,itcnt
if(nint(prct(j,i,1)).eq.nint(pa(k)).and.
& nint(prct(j,i,2)).eq.nint(pb(k))) then
ii = k
goto 1001
endif
enddo
1001 continue
if(prct(j,i,1).lt.-9990.0.or.prct(j,i,2).lt.-9990.0.or.
& ii.eq.0) then
grd(i,j) = -9999.0
else
grd(i,j) = bias(ii)
endif
enddo ! j loop
enddo ! i loop
c
c Write the gempak file to the disk. GDEDIT will put in the gempak file.
c
c generic "stuff" follows
open(unit=20,file='file3.out',status='unknown')
write(20,1) ' GRID IDENTIFIER:'
write(20,1) ' TIME1 TIME2 LEVL1 LEVL2 VC
&ORD PARM'
write(20,1) dateetc(1:72)
write(20,*) ' AREA: GRID GRID SIZE: ',igemx,igemy
write(20,*) 'COLUMNS: 1 ',igemx,' ROWS: 1 ',igemy
write(20,*) ' '
write(20,1) ' Scale factor: 10** 0'
write(20,1) ' '
write(20,1) ' '
jstop = 8
jdel = 7
20 continue
if(jstop.eq.8) then
write(20,29) ' COLUMN:',(j,j=jstop-jdel,jstop)
else
write(20,29) ' ',(j,j=jstop-jdel,jstop)
endif
jstop = jstop + 8
if(jstop.le.igemx) then
goto 20
elseif(jstop.ne.(igemx+8)) then
jdel = igemx - (jstop - 8) - 1
jstop = igemx
goto 20
else
jstop = 8
jdel = 7
endif
29 format(a8,2x,i5,7(4x,i5))
jstop = 8
jdel = 7
do i = igemy,1,-1
30 continue
c write(*,*) 'row,jstart,jstop= ',i, jstop-jdel, jstop
if(jstop.eq.8) then
write(20,31) ' ROW',i,(grd(j,i),j=jstop-jdel,jstop)
else
write(20,32) ' ',' ',(grd(j,i),j=jstop-jdel,jstop)
endif
jstop = jstop + 8
if(jstop.le.igemx) then
goto 30
elseif(jstop.ne.(igemx+8)) then
jdel = igemx - (jstop - 8) - 1
jstop = igemx
goto 30
else
jstop = 8
jdel = 7
endif
enddo
31 format(a4,i3,1x,8(1x,f8.2))
32 format(a4,a3,1x,8(1x,f8.2))
close(unit=20)
c
stop 'Done adjusting probabilities'
end