program calibrate2 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,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,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:13) = 'combine_pops_' namcomb(14:15) = frun(1:2) namcomb(16:16) = '_' namcomb(17:20) = fhour(1:4) namcomb(21:24) = '.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:63) = 'CALTRWHR' 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=namcomb(1:len_trim(namcomb)),status='old') do ii=1,itcnt read(11,*) pa(ii),pb(ii) read(11,*) dum1 read(11,*) dum1,dum2 read(11,*) dum1,dum2 bias(ii) = (dum2/(dum1+dum2))*100. 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