program add_sub_pcpn
c
c         Just reads the ascii gempak files and adds or subtracts grids.  DRB

       parameter (igemy=129, igemx=185) ! GEMPAK grid size
c
       real prct(igemy,igemx,2),grd(igemx,igemy)
       integer igemy,igemx,jstop,jdel
       character*150 pctfila,pctfilb
       character dummy*20,date*8,fhour*4,frun*2,dateetc*72,gfunc*4,
     &          dowhat*3
       
       call getarg(1,pctfila) ! forecast file
       call getarg(2,pctfilb) ! forecast file2
       call getarg(3,date) ! YYYYMMDD
       call getarg(4,frun) ! 00 or 12
       call getarg(5,fhour) ! f000 ... f063...
       call getarg(6,gfunc) ! 4 character output name
       call getarg(7,dowhat) ! add or sub or mul (add or subtract grids)
       if(dowhat(1:1).eq.'a'.or.dowhat(1:1).eq.'A') then
        dowhat(1:1) = 'a'
	write(*,*) 'Simply adding the two files together.'
       elseif(dowhat(1:1).eq.'s'.or.dowhat(1:1).eq.'S') then
        dowhat(1:1) = 's'
	write(*,*) 'Subtract the first file FROM the second file.'
       elseif(dowhat(1:1).eq.'m'.or.dowhat(1:1).eq.'M') then
        dowhat(1:1) = 'm'
	write(*,*) 'Multiply the first file and the second file.'
       endif

       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:59) = gfunc(1:4)
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
c Now...add or subtract here...
c
       do i=1,igemx
       do j=1,igemy
        if(prct(j,i,1).lt.-0.001.or.prct(j,i,2).lt.-0.001) then
	 grd(i,j) = -9999.0
	 goto 777
	endif
        if(dowhat(1:1).eq.'a') then
	 grd(i,j) = prct(j,i,1)+prct(j,i,2)
	elseif(dowhat(1:1).eq.'s') then
	 grd(i,j) = prct(j,i,2)-prct(j,i,1)
	elseif(dowhat(1:1).eq.'m') then
	 grd(i,j) = prct(j,i,2)*prct(j,i,1)*.01 ! assume they are probabilities...
	else
	 stop 'I AM SO CONFUSED'
	endif
 777    continue	
       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