#!/bin/csh -f
# data convert from GRIB to direct-access BINARY or GRIB
#       ftp:ncardata.ucar.edu /pub/reanalysis/monthly
#       ftp:ncardata.ucar.edu /datasets/ds090.2
#
# wgrib options
#       -nh  : no-header
#       -grib : GRIB
#       -ieee : BINARY, iy=1(NP) to 73(SP) ; 
#               read(11,rec=1) ((data(ix,iy),ix=1,nx),iy=1,ny)
cd /home/viper/gpv
#set date="`/bin/date -v-1d +%Y%m%d`"12
set date=2000080312
cd $date
foreach var ( h u v )
  set ifile=h${var}xa50
  set ofile=NH500${var}.${date}
if ( -e $ofile ) \rm $ofile
  echo $ofile
# *** All records
  wgrib -s $ifile | wgrib $ifile -i -ieee -nh -o $ofile > /dev/null
end
#
# +++ Fortran +++
cat >! prog.f << EOF
      program main
      parameter(nx=144,ny=37)
      real*4 data(nx,ny,3)

      nlen=ny*nx*4
c *** input data (Z)
      open(21,file='fort.21',
     & form='unformatted',access='direct',recl=nlen)
c *** input data (u)
      open(22,file='fort.22',
     & form='unformatted',access='direct',recl=nlen)
c *** input data (v)
      open(23,file='fort.23',
     & form='unformatted',access='direct',recl=nlen)
c *** output data
c     open(51,file='fort.51',
c    & form='unformatted',access='direct',recl=nlen)

      read(21,rec=1) ((data(ix,iy,1),ix=1,nx),iy=1,ny)
      read(22,rec=1) ((data(ix,iy,2),ix=1,nx),iy=1,ny)
      read(23,rec=1) ((data(ix,iy,3),ix=1,nx),iy=1,ny)

        do iy=1,ny
        ylat=90-2.5*(iy-1)
        do ix=1,nx
        xlon=2.5*(ix-1)
          write(51,*) xlon,ylat,(data(ix,iy,iv),iv=1,3)
        enddo	! END of ix loop
          write(51,*) 360.0,ylat,(data(1,iy,iv),iv=1,3)
        enddo	! END of iy loop

c         write(51,rec=1) ((data(ix,iy),ix=1,nx),iy=1,ny)

      stop
      end
EOF
f77 -o TMP prog.f
set ifile1=NH500h.$date
set ifile2=NH500u.$date
set ifile3=NH500v.$date
#if( -e fort.2\? ) rm fort.2?
\rm fort.2?
  ln -s $ifile1 fort.21
  ln -s $ifile2 fort.22
  ln -s $ifile3 fort.23
  TMP
#  mv -f fort.51 $ofile
\rm prog.f TMP
\rm $ifile1 $ifile2 $ifile3 gala.f
\rm fort.2?
#
# +++ GMT +++
# *** NH polar stereographic mapping ***
###### input & output file path setting ######
set ifile=fort.51
set ofile=NH500Z.${date}.ps
echo $ifile
echo $ofile
#############################  GMT program  ################################
set rr=0/360/20/90		# mapping range
set jj=180/90/3.0/20		# -Js option parameter
set cint=50			# Contour interval
@ aint = $cint * 4		# Annotation interval
set ll=0/10000			# limitation of contour lines
cat >! text.in << EOF
-0.5  7.4 12 0 4 1 !! 20N-NP version !!
-0.5  7.2 12 0 4 1 JMA GPV ; ${date}+00 (analysis)
-0.5  7.0 12 0 4 1 contour=${cint}m ; annotation=${aint}m
 6.0  6.0 15 0 5 3 500hPa Z & V
 0.0  6.0 15 0 5 1 ${date}+00
EOF
gmtset UNIX_TIME_POS 0.0/-0.5
psxy -R0/8.5/0/11 -Jx1 -L -W3 -U'M. HAYASAKI ; Tsukuba Univ.' -P -N -K -X0.7 -Y2.5 << EOF >! $ofile
 -0.7 -1.0
  7.0 -1.0
  7.0  7.0
 -0.7  7.0
EOF
# +++ scale vector +++
psxy -R -Jx -Sv0.0005/0.03/0.015 -L -G0 -N -W1 -O -K << EOF >> $ofile
  0.45  6.8  180   0.2
EOF
pstext -R -Jx -N -O -K << EOF >> $ofile
  0.5 6.8  12  0  4  5  : 20 m/s
EOF
# +++ header +++
pstext text.in -R0/8.5/0/11 -Jx1 -O -N -K -X0.5 -Y0.5 >> $ofile
pscoast -R$rr -Js$jj -B60g30/30g30 -Dc -A3000 -G225 -W1 -O -K >> $ofile
# +++ Z (contour) +++
gawk '{print $1,$2,$3}' $ifile | xyz2grd -Gz.cdf -R$rr -I2.5 -L
grdcontour z.cdf -R -Js -C$cint -A$aint -L$ll -Wc2 -Wa5 -O -K >> $ofile
# +++ u,v (vector) +++
gawk '{print $1,$2,$4}' $ifile | xyz2grd -Gu.cdf -R0/360/0/60 -I5.0
gawk '{print $1,$2,$5}' $ifile | xyz2grd -Gv.cdf -R0/360/0/60 -I5.0
grdvector u.cdf v.cdf -R$rr -Js -G0 -Q0.0005/0.03/0.015 -S100 -T -W -O -K >> $ofile
gawk '{print $1,$2,$4}' $ifile | xyz2grd -Gu.cdf -R0/360/60/85 -I10.0/5.0
gawk '{print $1,$2,$5}' $ifile | xyz2grd -Gv.cdf -R0/360/60/85 -I10.0/5.0
grdvector u.cdf v.cdf -R$rr -Js -G0 -Q0.0005/0.03/0.015 -S100 -T -W -O >> $ofile
\rm -f *.cdf text.in
\rm fort.51
find . -name .gmtcommands -exec rm {} \;
find . -name .gmtdefaults -exec rm {} \;
exit

