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 (600*600), U33(600*600)
dimension acp361 (600*600)
dimension acp363 (600*600)
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)
IHRNO = 0
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
if (IP1-IP2.eq.IHRNO) 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