program calibrate2_12hr
c
c Reads the 12 hr tables and takes calenh to calenh12. DRB. 6/29/2004.
c
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,4),prob,pa(3000),pb(3000),rk,bias(3000),
& grd(igemx,igemy),dum1,dum2,dumtest,pc(3000),pd(3000)
integer igemy,igemx,itcnt,jstop,jdel,iv,jv,iiv,jjv
character*150 pctfila,pctfilb,pctfilc,pctfild
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,pctfilc) ! forecast file2
call getarg(4,pctfild) ! forecast file2
call getarg(5,date) ! YYYYMMDD
call getarg(6,frun) ! 09 or 21
call getarg(7,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:63) = 'CALENH12'
c
c fill all possible forecast values...
itcnt = 0
c do i=0,50,10
c do j = 0,50,10
c do ii = 0,50,10
c do jj = 0,50,10
do i=1,7
if(i.eq.1) iv = 0
if(i.eq.2) iv = 5
if(i.eq.3) iv = 10
if(i.eq.4) iv = 20
if(i.eq.5) iv = 30
if(i.eq.6) iv = 40
if(i.eq.7) iv = 50
do j =1,7
if(j.eq.1) jv = 0
if(j.eq.2) jv = 5
if(j.eq.3) jv = 10
if(j.eq.4) jv = 20
if(j.eq.5) jv = 30
if(j.eq.6) jv = 40
if(j.eq.7) jv = 50
do ii=1,7
if(ii.eq.1) iiv = 0
if(ii.eq.2) iiv = 5
if(ii.eq.3) iiv = 10
if(ii.eq.4) iiv = 20
if(ii.eq.5) iiv = 30
if(ii.eq.6) iiv = 40
if(ii.eq.7) iiv = 50
do jj = 1,7
if(jj.eq.1) jjv = 0
if(jj.eq.2) jjv = 5
if(jj.eq.3) jjv = 10
if(jj.eq.4) jjv = 20
if(jj.eq.5) jjv = 30
if(jj.eq.6) jjv = 40
if(jj.eq.7) jjv = 50
itcnt = itcnt + 1
pa(itcnt) = float(iv) ! fcst1
pb(itcnt) = float(jv) ! fcst2
pc(itcnt) = float(iiv) ! fcst1
pd(itcnt) = float(jjv) ! fcst2
enddo
enddo
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 2. Read in the percent2 file:
open(unit=10,file=pctfilc,status='old',form='formatted')
write(*,*) 'Reading the GEMPAK forecast file 3 now: ',pctfilc
1335 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,3),i=1,igemx)
else
read(10,*) dummy(1:8),irow,(prct(j,i,3),i=1,igemx)
endif
enddo
else
goto 1335
endif
close(unit=10)
c 2. Read in the percent2 file:
open(unit=10,file=pctfild,status='old',form='formatted')
write(*,*) 'Reading the GEMPAK forecast file 4 now: ',pctfild
1445 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,4),i=1,igemx)
else
read(10,*) dummy(1:8),irow,(prct(j,i,4),i=1,igemx)
endif
enddo
else
goto 1445
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,4
c if(prct(j,i,k).lt.-1.0) prct(j,i,k) = -9999.0
c if(prct(j,i,k).ge.-1.0.and.prct(j,i,k).lt.5.0)prct(j,i,k)= 0.0
c if(prct(j,i,k).ge.5.0.and.prct(j,i,k).lt.15.) prct(j,i,k)=10.0
c if(prct(j,i,k).ge.15.0.and.prct(j,i,k).lt.25.)prct(j,i,k)=20.0
c if(prct(j,i,k).ge.25.0.and.prct(j,i,k).lt.35.)prct(j,i,k)=30.0
c if(prct(j,i,k).ge.35.0.and.prct(j,i,k).lt.45.)prct(j,i,k)=40.0
c if(prct(j,i,k).ge.45.0.and.prct(j,i,k).lt.55.)prct(j,i,k)=50.0
c if(prct(j,i,k).ge.55.0.and.prct(j,i,k).lt.65.)prct(j,i,k)=60.0
c if(prct(j,i,k).ge.65.0.and.prct(j,i,k).lt.75.)prct(j,i,k)=70.0
c if(prct(j,i,k).ge.75.0.and.prct(j,i,k).lt.85.)prct(j,i,k)=80.0
c if(prct(j,i,k).ge.85.0.and.prct(j,i,k).lt.95.)prct(j,i,k)=90.0
c if(prct(j,i,k).ge.95.0) prct(j,i,k)=100.0
if(prct(j,i,k).lt.2.5.and.prct(j,i,k).gt.-1.0) prct(j,i,k)=0.0
if(prct(j,i,k).ge.2.5.and.prct(j,i,k).lt.7.5) prct(j,i,k)=5.0
if(prct(j,i,k).ge.7.5.and.prct(j,i,k).lt.15.0) prct(j,i,k)=10.0
if(prct(j,i,k).ge.15.0.and.prct(j,i,k).lt.25.0)prct(j,i,k)=20.0
if(prct(j,i,k).ge.25.0.and.prct(j,i,k).lt.35.0)prct(j,i,k)=30.0
if(prct(j,i,k).ge.35.0.and.prct(j,i,k).lt.45.0)prct(j,i,k)=40.0
if(prct(j,i,k).ge.45.0) prct(j,i,k)=50.0
enddo
enddo
enddo
c
c Read the file that contains the dependency coefficient and the bias...
c
open(unit=11,file='combine_pops_12hr.out',status='old')
do ii=1,itcnt
read(11,*) pa(ii),pb(ii),pc(ii),pd(ii)
read(11,*) dum1
read(11,*) dum1,dum2
read(11,*) dum1,dum2
dumtest = dum1 + dum2
if(dumtest.le.0.0e-8) then
bias(ii) = -9999.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)).and.
& nint(prct(j,i,3)).eq.nint(pc(k)).and.
& nint(prct(j,i,4)).eq.nint(pd(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.
& prct(j,i,3).lt.-9990.0.or.prct(j,i,4).lt.-9990.0.or.
& ii.eq.0) then
grd(i,j) = -9999.0
else
if(bias(ii).gt.-1.0) then
grd(i,j) = bias(ii)
else
c grd(i,j) = prct(j,i,1)+prct(j,i,2)+prct(j,i,3)+prct(j,i,4)
c grd(i,j) = grd(i,j)*0.25
grd(i,j)=max(prct(j,i,1),prct(j,i,2),prct(j,i,3),prct(j,i,4))
endif
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='file5.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 12 hour probabilities'
end