c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ c GRIB to GrADS MAIN PROGRAM; c By Pejanovicz Aug-98 Malta c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ c purposes: c to read grib files and to store them into whished grads format c according to the control grads files. c arguments [YYMMDDSS] c 1) yymmddhh if in corresponding activectl.LST tf=-999 c 2) no arguments if in activectl.LST tf=-1 c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ c note: c 1 run this program: grb2grads [00032500] c 2 make: g77 -c -O1 tograds.g; g77 -O1 grb2grads.f tograds.o -o grb2grads c 3 sugestion: when you make grb2grads copy this exe version to c /usr/local/grads/bin to be seen from all others directories c 4 check activectl.LST; additional sintax is c *indset /..../...../ENS/grib ( place for grib) c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ c--------------------------- include 'tograds.inc' parameter (MXRECN=20000) c--------------------------- character * 8 Cdate8 c--------------------------- character * 150 dirn character * 620 str character * 5 c5 character * 2 c2 character * 3 cdlmd,cdphd dimension id7(7),id7h(7), ar7(7) dimension IPDS(120), IBMS(130), GDS(150) dimension ALL (120) c dimension NPOS (MXRECN), VAL (400*200), U33(400*200) dimension NPOS (MXRECN), VAL (800*800), U33(800*800) dimension acp361 (800*800) dimension acp363 (800*800) call getarg (1, Cdate8) call getarg (2, cdlmd) call getarg (3, cdphd) NCOUNT = 0 read (Cdate8(1:2),'(i2)') id7(1) read (Cdate8(3:4),'(i2)') id7(2) read (Cdate8(5:6),'(i2)') id7(3) read (Cdate8(7:8),'(i2)') id7(4) c print *,Cdate8 c print *,id7 RLEV = 0 id7(5) = 0 id7(6) = 0 id7(7) = 0 do k=1, MXRECN NPOS (k) = 0 enddo c open (99,file='sst.datll',form='unformatted') c ------------------------------------------------ call getctls (1, id7) call mklist (NF) if (NF.eq.0) call fmess('NF=0') c ------------------------------------------------ open (12,file='L.out',status='old') call system ('rm S.out') 10 call clears (dirn) read (12,'(a)',end=99) dirn c ------------------------------------------------ nb = index (dirn,' ') - 1 print *,dirn(1:nb) c ------------------------------------------------ c/usr/local/grads/bin/wgrib call system c +('/usr/local/grads/bin/gribscan -gv -gd -ig '//dirn(1:nb)// c +('/usr/local/grads/binOLD/gribscan -gv -gd -ig '//dirn(1:nb)// c +('/usr/local/grads/bin/gribscan -gv -gd -ig '//dirn(1:nb)// ! +('/opt/grads-2.0.a5/bin/gribscan -gv -gd -ig '//dirn(1:nb)// ! + ' >> S.out') +('/usr/local/grads-2.0.1/bin/gribscan -gv -gd -ig '//dirn(1:nb)// + ' >> S.out') c ------------------------------------------------ goto 10 c ------------------------------------------------ 99 continue rewind (12) c ------------------------------------------------ c ------------------------------------------------ print *,' ............. ok ' c ------------------------------------------------ open (25,file='S.out') c ------------------------------------------------ IFILE = 0 c ------------------------------------------------ NSOUT = 0 20 call clears(str) c write (33,*) '*******************************' read (25,'(a)',end=88,err=88) str NSOUT = NSOUT + 1 c print *, ' ... Nrec=',NSOUT c ------------------------------------------------ ips = 1 if( index(str,"PDS #").eq.0) goto 20 call s2iass ( str, ips,'PDS #', IREC) c ------------------------------------------------ c ne file ! if (IREC.eq.1) then IFILE = IFILE + 1 call clears (dirn) read (12,'(a)') dirn NB = index (dirn,' ') - 1 endif c ------------------------------------------------ call s2i ( str, ips, IGRID) call s2i ( str, ips, IPX) call s2i ( str, ips, ITL) call s2i ( str, ips, IL) call s2i ( str, ips, I61) call s2i ( str, ips, I62) call s2i ( str, ips, I63) call s2iass ( str, ips,'BMSFL', IBMDUM) call s2i ( str, ips, ID7 (1)) call s2i ( str, ips, ID7 (2)) call s2i ( str, ips, ID7 (3)) call s2i ( str, ips, ID7 (4)) call s2i ( str, ips, IBMDUM) call s2i ( str, ips, IBMDUM) call s2i ( str, ips, IP1 ) call s2i ( str, ips, IP2 ) call s2i ( str, ips, IACC ) c call s2i ( str, ips, IACC ) c ID7(4) = id7h(3) ID7 (5) = IP1 if ( IPX.eq.71.and.IP1-IP2.lt.0) goto 20 c------------------------------------------------------- if ( IACC.eq.4) Id7(5) = IP2 c------------------------------------------------------- call s2iass ( str, ips,'GDS ', IBMDUM ) call s2i ( str, ips, NX) call s2i ( str, ips, NY) call s2r ( str, ips, WB) call s2r ( str, ips, EB) call s2r ( str, ips, DLM) call s2r ( str, ips, SB) call s2r ( str, ips, RNB) call s2r ( str, ips, DPH) c------------------------------------------------------- c print *,IBMDUM,nx,ny if ( IBMDUM.eq.203) goto 203 c------------------------------------------------------- c print *,id7, ' IP1,IP2,IACC', IP1, IP2, IACC c print *, IGRIDLL, NX, WB,EB,DLM c print *, IGRIDLL, Ny, SB,RNB,DPH IGRIDLL = 0 c gauss-case! if ( DLM.eq.937.0.or.DLM.eq.938.0) then DLM = 0.9375 DPH = DLM EB = -DLM IGRIDLL = -1 endif c---------------------------------------------------- IYREV = 0 NSFLAG = 0 c thinned-case! if ( IGRID.ge.37.and.IGRID.le.44) IGRIDLL = -2 if ( IGRIDLL.eq.-2) then NX = NY DLM = DPH ENDIF c ---------------------------------- if ( EB.lt.WB) EB = EB + 360. c ---------------------------------- if ( RNB.lt.SB) then HLP = SB SB = RNB RNB = HLP IYREV = 1 endif if (SB.lt.0) NSFLAG = -1 c ---------------------------------- ar7(1) = 0 ! warning! ar7(2) = NX ar7(3) = WB ar7(4) = DLM ar7(5) = NY ar7(6) = SB ar7(7) = DPH c wafs RLEV =il if ( itl.ne.100) RLEV = 0 c eta2d grid ! goto 1200 c etagrid 203 c o n t i n u e IGRIDLL = -999 c read (cdlmd,"(*)") irdlmd c read (cdphd,"(*)") irdphd ipos = 1 c call s2i(cdlmd,ipos,irdlmd) c call s2i(cdphd,ipos,irdphd) c print *," now --- ",cdlmd," X ",cdphd c print *," now --- ",irdlmd,irdphd c print *," dlm dph ",dlm,dph c if ( irdlmd*irdphd.ne.0) then c dlm = 1./irdlmd c dph = 1./irdphd c else if ( DLM.eq.166.) DLM= 166.66666 if ( DLM.eq.083.) DLM= 083.33333 if ( DLM.eq.041.) DLM= 041.66666 if ( DPH.eq.166.) DPH= 166.66666 if ( DPH.eq.083.) DPH= 083.33333 if ( DPH.eq.041.) DPH= 041.66666 DPH = 0.001 * DPH DLM = 0.001 * DLM c endif c print *," dlm dph ",dlm,dph WBD = -(NX/2-0)*2*DLM SBD = -(NY/2-0)*DPH TLM0d = EB TPH0d = RNB ar7(1) = 2 ! warning! ar7(2) = TLM0D ar7(3) = WBD ar7(4) = DLM ar7(5) = TPH0d ar7(6) = SBD ar7(7) = DPH c print *,' ar7=',ar7 c print *,NCTL, ' ip,itl,il ',ipx,itl,il,NX,NY c print *,NCTL, ' id7 ' ,(id7(m),m=1,5) c---------------------------------------- 1200 c o n t i n u e write (c5,'(i5)') IREC call system c +('/usr/local/grads/bin/wgrib -v '//dirn(1:nb)//' -d '//c5//' >1.0 ! +(' /opt/grads-2.0.a5/bin/wgrib -v '//dirn(1:nb)//' -d '//c5//' >1.0 ! +') +('/usr/local/grads-2.0.1/bin/wgrib '//dirn(1:nb)//' -d '//c5//' >1.0') ! stop cpj IRET is invoked! I3 = 0 if ( mod(id7(5),6).ne.0) I3=1 if (IGRIDLL.eq.-2) then call MW3FT33x(IRET, NSFLAG, val) else if (IGRIDLL.ge.-1) then call readdump (IRET,IYREV, NX,NY,val) endif EL = IL JHR = id7(5) c############################################### if (IBMDUM.ne.203) goto 1700 call readdump (IRET,IYREV, NX,NY,val) c############################################### if ( IPX.eq.33) then do j=1,NX*NY U33(j) = val(j) enddo c ======= goto 20 c ======= endif if ( IPX.eq.61.and.I3.eq.1) then do j=1,NX*NY acp361(j) = val(j) enddo endif if ( IPX.eq.63.and.I3.eq.1) then do j=1,NX*NY acp363(j) = val(j) enddo endif c===================================================== call vrfc(IPX,ITL,IL,JHR,6,IPXVRF,ITLVRF,ILVRF) ELVRF = ILVRF c===================================================== c only wind! if ( IPX.eq.34) then IWFL = 2 if(ITL.eq.105) IWFL = 1 CALL WGRADS(id7,33,34,IWFL,ITL,IL,ar7,NX,NY,1,EL,U33,VAL) IPXV1 = 33 IPXV2 = 34 if ( IPXVRF.gt.33000) then IPXV1 = IPXVRF-1000 IPXV2 = IPXVRF endif c print *,' wind IPXVRF ... ',IPXVRF,IPXV1,IPXV2 CALL WGRADS(id7,IPXV1,IPXV2,IWFL, 99999 . ,ILVRF,ar7,NX,NY,1,ELVRF,U33,VAL) c print *,' endofwind IPXVRF ... ',IPXVRF,IPXV1,IPXV2 else c main write! on 203 grid cxcxcx CALL WGRADS(id7,IPX,0,0,ITL,IL,ar7,NX,NY,1,EL,VAL,-99.) c verification if (IPXVRF.gt.0) then if (IPXVRF.eq.61.or.IPXVRF.eq.63) then c print *,NXNY,' rain verif ... ',IPXVRF,ITLVRF,ILVRF do j=1,NX*NY if( IPXVRF.eq.61)val (j) = val(j) + acp361(j) if( IPXVRF.eq.63)val (j) = val(j) + acp363(j) enddo endif c print *,id7,' p lev ',IPXVRF,ILVRF CALL WGRADS(id7,IPXVRF,0,0, 99999 . ,ILVRF,ar7,NX,NY,1,ELVRF,VAL,-99) endif endif goto 20 c############################################### 1700 c o n t i n u e ! no eta CASE c############################################### c main write CALL WGRADS (id7,IPX,0,0,ITL,IL,ar7,NX,NY,1,EL,val,-99.) cverification call vrfc(IPX,ITL,IL,JHR,6,IPXVRF,ITLVRF,ILVRF) ELVRF = ILVRF c verification if ( IPXVRF.gt.0) then CALL WGRADS . (id7,IPXVRF,0,0, 99999,ILVRF,ar7,NX,NY,1,ELVRF,VAL,-99.) endif c endofverif! c############################################### c---------------------------------------- goto 20 c---------------------------------------- 88 continue c---------------------------------------- czzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz czzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz stop 205 continue end c================================================================ subroutine vrfc(IPX,ITL,IL,IHR,IHRINC,IPXVRF,ITLVRF,ILVRF) c================================================================ ILVRF = -99 IPXVRF= -99 ITLVRF= 99999 if ( mod( IHR,IHRINC).ne.0) return ILVRF = 000 + IHR if ( ITL.eq.100) then IPXVRF = IPX*1000 + mod(IL,1000) endif if (IPX.eq.02.and.ITL.eq.102.and.IL.eq.0) IPXVRF = IPX if (IPX.eq.61.and.ITL.eq.001.and.IL.eq.0) IPXVRF = IPX if (IPX.eq.63.and.ITL.eq.001.and.IL.eq.0) IPXVRF = IPX if (IPX.eq.65.and.ITL.eq.001.and.IL.eq.0) IPXVRF = IPX if (IPX.eq.52.and.ITL.eq.105.and.IL.eq.2) IPXVRF = IPX if (IPX.eq.11.and.ITL.eq.105.and.IL.eq.2) IPXVRF = IPX if (IPX.eq.33.and.ITL.eq.105.and.IL.eq.10) IPXVRF = IPX if (IPX.eq.34.and.ITL.eq.105.and.IL.eq.10) IPXVRF = IPX if (IPX.eq.71.and.ITL.eq.200.and.IL.eq.10) IPXVRF = IPX return end c===================================================== c--------------------------- subroutine mklist (NF) include 'tograds.inc' c--------------------------- character *250 str character *250 fstr character *10 hstr j1 = 1 m = j1 call mvwhile (' ','ne',m,+1,indset) j2 = m - 1 if ( j2.lt.j1) call fmess ('indset name!') if (indset(j2:j2).ne.'/') then j2 = j2 +1 indset(j2:j2)='/' endif call system ('ls -l '//indset(j1:j2)//' > L.out') print *,' indset=',indset(j1:j2),'**' NF = 0 open (12,file='L.out',status='old',err=299) open (22,file='L.out1',status='unknown') 10 call clears (str) read (12,'(a)',end=99) str if ( index(str,':').eq.0) goto 10 NF = NF + 1 m = 250 call mvwhile (' ','eq',m,-1,str) i2 = m call mvwhile (' ','ne',m,-1,str) i1 = m +1 lend = j2-j1+1 lens = i2-i1+1 L = lend + lens fstr(1:lend) = indset(j1:j2) fstr(lend+1:L) = str(i1:i2) c print *,'***',indset(j1:j2),'****' c print *,'***',str(i1:i2),'****' c print *,'***',fstr(1:L),'****' call system c +('/usr/local/grads/bin/wgrib -v '//fstr(1:L)//' -d 1 > wc.out1') ! +(' /opt/grads-2.0.a5/bin/wgrib -v '//fstr(1:L)//' -d 1 > wc.out1') +(' /usr/local/grads-2.0.1/bin/wgrib '//fstr(1:L)//' -d 1 > wc.out1') call system ('wc -l wc.out1 > wc.out') open (62,file='wc.out',status='old') read (62,'(a)') hstr close (62) ips = 1 call s2i ( hstr, ips, NLWC ) if ( NLWC.eq.0) goto 10 write (22,'(a)') fstr(1:L) goto 10 99 continue close (12) close (22) call system('mv L.out1 L.out') call system('rm wc.out*') 299 return end subroutine check_truncation ( DLMIN,DLM) c1/6 if ( DLMIN.eq.166) DLM = 016.66666 c1/12 if ( DLMIN.eq.083) DLM = 083.33333 c1/24 if ( DLMIN.eq.041) DLM = 041.66666 c1/30 if ( DLMIN.eq.033) DLM = 033.33333 c1/36 if ( DLMIN.eq.027) DLM = 027.77777 c1/48 c if ( DLMIN.eq.020) DLM = 020.08333 c1/60 if ( DLMIN.eq.016) DLM = 016.66666 c1/72 if ( DLMIN.eq.013) DLM = 013.88888 return end