program power_mean_den
c
c Hard coded for the ligtning intesity (21 layers).  DRB.  Oct 4 2005.
c
c Compute the power mean of all the predictors.  This gives all layers equal weight but favors those
c layers with a higher probability.  A power = 2 is called a quadratic mean (or RMS mean), while 
c a power = 1 returns the arithmetic mean.  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,21),prob,pa(3000),pb(3000),rk,bias(3000),
     &       grd(igemx,igemy),dum1,dum2,pc(3000),pd(3000),power
       integer igemy,igemx,itcnt,jstop,jdel,iv,jv,iiv,jjv,ncnt
       character*150 pctfil(21),pctin
       character dummy*20,date*8,fhour*4,frun*2,dateetc*72,gfunc*6
       
       ncnt = 21
       power = 5.0
       
       do i=1,ncnt
        call getarg(i,pctfil(i))
       enddo
       call getarg(ncnt+1,date) ! YYYYMMDD
       call getarg(ncnt+2,frun) ! 09 or 21
       call getarg(ncnt+3,fhour) ! f000 ... f063...
       call getarg(ncnt+4,gfunc) ! 6 char gempak vrbl

       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:55+len_trim(gfunc)) = gfunc(1:len_trim(gfunc))
c 
 1    format(a)
c Now read them in...
c
c 1. Read in the percent file:
      icnt = 0
      write(*,*) ' '
      write(*,*) 'Calculate the power mean for ',ncnt,' input files'
      write(*,*) 'using a power of ',power
      write(*,*) ' '
      
 6543 continue      
      icnt = icnt + 1
      pctin = ''
      pctin = pctfil(icnt)
      open(unit=10,file=pctin,status='old',form='formatted')
      write(*,*) 'Reading the GEMPAK forecast file now: ',pctin
 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,icnt),i=1,igemx)
	 else
	 read(10,*) dummy(1:8),irow,(prct(j,i,icnt),i=1,igemx)
	 endif
	enddo
       else
	goto 7775
       endif
       close(unit=10)
       if(icnt.lt.ncnt) goto 6543
       
c
c Now the grids are in...calculate the power mean! 
       do i=1,igemx
       do j=1,igemy
        grd(i,j) = 0.0
       do k=1,ncnt       
         if(prct(j,i,k).lt.-9000.0) then
	  grd(i,j) = -9999.0
	  goto 5432
	 else
	  grd(i,j) = grd(i,j) + (prct(j,i,k)**power)
         endif
       enddo
       grd(i,j) = (grd(i,j)/float(ncnt))**(1.0/power)
 5432  continue       
       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='power_mean_den.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 calculating the power mean.'
       end