program add_sub_maxmin 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),grdl(igemx,igemy) integer igemy,igemx,jstop,jdel,restart,lyrval character*150 pctfila,pctfilb character dummy*20,date*8,fhour*4,frun*2,dateetc*72,gfunc*6, & 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) ! 6 character output name call getarg(7,dowhat) ! add or sub or max or min (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:3).eq.'max'.or.dowhat(1:3).eq.'MAX') then dowhat(1:3) = 'max' write(*,*) 'Find the maximum of the two files.' elseif(dowhat(1:3).eq.'min'.or.dowhat(1:3).eq.'MIN') then dowhat(1:3) = 'min' write(*,*) 'Find the minimum of the two files.' 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:61) = gfunc(1:6) 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 read in the grid that contains which layer the max probability came from. If it does not c exist, then this is a brand new run... if(dowhat(1:3).eq.'max') then open(unit=10,file='layermax.out',status='old',err=712) write(*,*) 'Reading the GEMPAK file: layermax.out' 7779 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),(grdl(i,j),i=1,igemx) else read(10,*) dummy(1:8),irow,(grdl(i,j),i=1,igemx) endif enddo else goto 7779 endif close (10) restart = 0 lyrval = -9999 if(pctfila(8:9).eq.'01') lyrval = 1 if(pctfila(8:9).eq.'02') lyrval = 2 if(pctfila(8:9).eq.'03') lyrval = 3 if(pctfila(8:9).eq.'04') lyrval = 4 if(pctfila(8:9).eq.'05') lyrval = 5 if(pctfila(8:9).eq.'06') lyrval = 6 if(pctfila(8:9).eq.'07') lyrval = 7 if(pctfila(8:9).eq.'08') lyrval = 8 if(pctfila(8:9).eq.'09') lyrval = 9 if(pctfila(8:9).eq.'10') lyrval = 10 if(pctfila(8:9).eq.'11') lyrval = 11 if(pctfila(8:9).eq.'12') lyrval = 12 if(pctfila(8:9).eq.'13') lyrval = 13 if(pctfila(8:9).eq.'14') lyrval = 14 if(pctfila(8:9).eq.'15') lyrval = 15 if(pctfila(8:9).eq.'16') lyrval = 16 if(pctfila(8:9).eq.'17') lyrval = 17 if(pctfila(8:9).eq.'18') lyrval = 18 if(pctfila(8:9).eq.'19') lyrval = 19 if(pctfila(8:9).eq.'20') lyrval = 20 if(pctfila(8:9).eq.'21') lyrval = 21 if(pctfila(8:9).eq.'22') lyrval = 22 if(pctfila(8:9).eq.'23') lyrval = 23 if(pctfila(8:9).eq.'24') lyrval = 24 if(pctfila(8:9).eq.'25') lyrval = 25 goto 713 712 continue c here due to the file not existing...set everything to missing... write(*,*) 'Setup the GEMPAK file: layermax.out' restart = 1 do i=1,igemx do j=1,igemy grdl(i,j) = -9999.0 enddo enddo 713 continue endif c c Now...add or subtract here... c do i=1,igemx do j=1,igemy if(nint(prct(j,i,1)).eq.-9999.or. & nint(prct(j,i,2)).eq.-9999) then grd(i,j) = -9999.0 goto 1001 endif if(prct(j,i,1).lt.-0.001) prct(j,i,1) = 0.0 if(prct(j,i,2).lt.-0.001) prct(j,i,2) = 0.0 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:3).eq.'max') then grd(i,j) = max(prct(j,i,2),prct(j,i,1)) if(restart.eq.1) then if(prct(j,i,2).gt.0.01.and.prct(j,i,2).ge.prct(j,i,1)) then grdl(i,j) = 2.0 endif if(prct(j,i,1).gt.0.01.and.prct(j,i,1).ge.prct(j,i,2)) then grdl(i,j) = 1.0 endif else if(prct(j,i,1).gt.0.01.and.prct(j,i,1).ge.prct(j,i,2)) then grdl(i,j) = float(lyrval) endif endif elseif(dowhat(1:3).eq.'min') then grd(i,j) = min(prct(j,i,2),prct(j,i,1)) else stop 'I AM SO CONFUSED' endif 1001 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') restart = 0 805 continue 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 c now output the max layer if this is a max application... if(dowhat(1:3).eq.'max'.and.restart.eq.0) then dateetc(62:62) = 'L' restart = 1 do i=1,igemx do j=1,igemy grd(i,j) = grdl(i,j) enddo enddo open(unit=20,file='layermax.out',status='unknown') goto 805 endif c stop 'Done adjusting probabilities' end