program calibrate2_hrly_rgn3
c
c Just like calibrate2_hrly, except reads the gridded output for "regional" coefficients. DRB, 8/18/2004.
c Just like calibrate.f, except reads the output of the combined probability table.
parameter (igemy=129, igemx=185) ! GEMPAK grid size
c
real prct(igemy,igemx,2),prob,pa(121),pb(121),rk,
& grd(igemx,igemy),dum1,dum2,bin(igemx,igemy),
& hit(121,2,igemx,igemy),bias(121,igemx,igemy)
integer igemy,igemx,itcnt,jstop,jdel
character*150 pctfila,pctfilb
character dummy*20,date*8,fhour*4,frun*2,dateetc*72,namcomb*50
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...
namcomb = ''
namcomb(1:18) = 'combine_pops_rgn3_'
namcomb(19:20) = frun(1:2)
namcomb(21:21) = '_'
namcomb(22:25) = fhour(1:4)
namcomb(26:29) = '.out'
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:67) = 'CALTRWHRRGN3'
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
ccc open(unit=11,file=namcomb(1:len_trim(namcomb)),status='old')
open(unit=11,file=namcomb(1:len_trim(namcomb)),status='old',
& form='unformatted')
ccc do ii=1,itcnt
read(11) pa
read(11) pb
read(11) bin
read(11) hit
close(11)
c read(11,*) pa(ii),pb(ii)
c read(11,*) dum1
c read(11,*) dum1,dum2
c read(11,*) dum1,dum2
do ii = 1,igemx
do jj = 1,igemy
do kk = 1,itcnt
if(hit(kk,1,ii,jj).lt.-9000.0.or.
& hit(kk,2,ii,jj).lt.-9000.0.or.
& bin(ii,jj).lt.-9000.0) then
bias(kk,ii,jj) = -9999.0
elseif(hit(kk,1,ii,jj).lt.0.0e-15.and.
& hit(kk,2,ii,jj).lt.0.0e-15) then
iiii = -9999
jjjj = -9999
c find closest grid point in the area with good data...
do iii=ii-20,ii+20
do jjj=jj-20,jj+20
if(iii.lt.1.or.jjj.lt.1.or.iii.gt.igemx.or.jjj.gt.igemy)
& goto 801
if(hit(kk,1,iii,jjj).ge.0.0e-15.and.
& hit(kk,2,iii,jjj).ge.0.0e-15)then
if(abs(ii-iii).lt.abs(ii-iiii).and.
& abs(jj-jjj).lt.abs(jj-jjjj)) then
iiii = iii
jjjj = jjj
endif
endif
801 continue
enddo
enddo
if(iiii.gt.-9000.and.jjjj.gt.-9000) then
dum1= (hit(kk,1,iiii,jjjj)/bin(iiii,jjjj))*100.0
dum2= (hit(kk,2,iiii,jjjj)/bin(iiii,jjjj))*100.0
bias(kk,ii,jj) = (dum2/(dum1+dum2))*100.
else
bias(kk,ii,jj) = -9999.0
endif
else
dum1= (hit(kk,1,ii,jj)/bin(ii,jj))*100.0
dum2= (hit(kk,2,ii,jj)/bin(ii,jj))*100.0
bias(kk,ii,jj) = (dum2/(dum1+dum2))*100.
endif
c bias(ii) = (dum2/(dum1+dum2))*100.
enddo
enddo
enddo
598 format(2f10.2)
600 format(f10.2)
601 format(8f10.2)
ccc enddo
ccc 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
c grd(i,j) = bias(ii)
grd(i,j) = bias(ii,i,j)
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