program combine_svr_dependence
c
c Apr 18, 2005...Use the Sangster and Hughes model to combine svr probabilities based on a
c                previously computed dependency parameter. DRB.
c
c Dec 10 2004... Simply reads the 8 grids that make up a 24h period and combines the 
c probabilities as if they are independent.  DRB.       
       
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,8),prob,pa(3000),pb(3000),rk,bias(3000),
     &       grd(igemx,igemy),dum1,dum2,pc(3000),pd(3000),dep,remove
       integer igemy,igemx,itcnt,jstop,jdel,iv,jv,iiv,jjv
       character*150 pctfila,pctfilb,pctfilc,pctfild,pctfile,pctfilf,
     &             pctfilg,pctfilh
       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,pctfile) ! forecast file2
       call getarg(6,pctfilf) ! forecast file2
       call getarg(7,pctfilg) ! forecast file2
       call getarg(8,pctfilh) ! forecast file2
       call getarg(9,date) ! YYYYMMDD
       call getarg(10,frun) ! 09 or 21
       call getarg(11,fhour) ! f000 ... f063...
       
       open(unit=11,file='dep_param.input',status='old',err=45)
       read(11,*) dep
       close(11)
       goto 46 
 45    continue
       dep = 1.0
 46    continue
       write(*,*) 'Using a dependency parameter of: ',dep

       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) = 'CALSVR24'
c 
 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)
      open(unit=10,file=pctfile,status='old',form='formatted')
      write(*,*) 'Reading the GEMPAK forecast file 5 now: ',pctfile
 9645 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,5),i=1,igemx)
	 else
	 read(10,*) dummy(1:8),irow,(prct(j,i,5),i=1,igemx)
	 endif
	enddo
       else
	goto 9645
       endif
       close(unit=10)
      open(unit=10,file=pctfilf,status='old',form='formatted')
      write(*,*) 'Reading the GEMPAK forecast file 6 now: ',pctfilf
 9745 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,6),i=1,igemx)
	 else
	 read(10,*) dummy(1:8),irow,(prct(j,i,6),i=1,igemx)
	 endif
	enddo
       else
	goto 9745
       endif
       close(unit=10)
      open(unit=10,file=pctfilg,status='old',form='formatted')
      write(*,*) 'Reading the GEMPAK forecast file 7 now: ',pctfilg
 9845 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,7),i=1,igemx)
	 else
	 read(10,*) dummy(1:8),irow,(prct(j,i,7),i=1,igemx)
	 endif
	enddo
       else
	goto 9845
       endif
       close(unit=10)
      open(unit=10,file=pctfilh,status='old',form='formatted')
      write(*,*) 'Reading the GEMPAK forecast file 8 now: ',pctfilh
 9945 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,8),i=1,igemx)
	 else
	 read(10,*) dummy(1:8),irow,(prct(j,i,8),i=1,igemx)
	 endif
	enddo
       else
	goto 9945
       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
         if(prct(j,i,1).lt.-0.1.or.prct(j,i,2).lt.-0.1) then
	  grd(i,j) = -9999.0
	 else
c	  grd(i,j) = ( (prct(j,i,1)*.01 + prct(j,i,2)*.01) -
c     &                   (prct(j,i,1)*.01*prct(j,i,2)*.01) )*100.0
          if(prct(j,i,1).ge.prct(j,i,2)) then
	   remove = ((prct(j,i,1)*.01)**dep)*(prct(j,i,2)*.01)
	   grd(i,j) = ( (prct(j,i,1)*.01 + prct(j,i,2)*.01) -
     &                   remove )*100.0	   
	  else
	   remove = ((prct(j,i,2)*.01)**dep)*(prct(j,i,1)*.01)
	   grd(i,j) = ( (prct(j,i,1)*.01 + prct(j,i,2)*.01) -
     &                   remove )*100.0	   
	  endif
         endif
       enddo
       enddo
c
       do k=3,8
       do i=1,igemx
       do j=1,igemy
         if(prct(j,i,k).lt.-0.1.or.grd(i,j).lt.-0.1) then
	  grd(i,j) = -9999.0
	 else
c	  grd(i,j) = ( (grd(i,j)*.01 + prct(j,i,k)*.01) -
c     &                   (grd(i,j)*.01*prct(j,i,k)*.01) )*100.0
          if(prct(j,i,k).ge.grd(i,j)) then
	   remove = ((prct(j,i,k)*.01)**dep)*(grd(i,j)*.01)
	   grd(i,j) = ( (prct(j,i,k)*.01 + grd(i,j)*.01) -
     &                   remove )*100.0	   
	  else
	   remove = ((grd(i,j)*.01)**dep)*(prct(j,i,k)*.01)
	   grd(i,j) = ( (prct(j,i,k)*.01 + grd(i,j)*.01) -
     &                   remove )*100.0	   
	  endif
         endif
       enddo
       enddo
       enddo       
       
c
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='file_svr_24h.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 making 24 hour severe probabilities'
       end