******* Phan loai cap TC cua JMA:
* 2 : Tropical Depression (TD)
* 3 : Tropical Storm (TS)
* 4 : Severe Tropical Storm (STS)
* 5 : Typhoon (TY)
* 6 : Extra-tropical Cyclone (L)
* 7 : Just entering into the responsible area of Japan Meteorological Agency (JMA)
* 8 : Not used
* 9 : Tropical Cyclone of TS intensity or higher
***********************

   "reinit"
   "set display white"; "c"
   "set mpdraw on"; "set frame on"; "set grid on";"set xlab on";"set ylab on";

    mark.2=0; mark.3=0; mark.4=0; mark.5=0; mark.6=0; mark.7=0; mark.9=0
*    colo.2=4; colo.3=10; colo.4=12; colo.5=2; colo.6=1; colo.7=12; colo.9=12
    colo.2=5; colo.3=7; colo.4=8; colo.5=6; colo.6=0; colo.7=12; colo.9=12
    thic.2=4; thic.3=4; thic.4=4; thic.5=4; thic.6=4; thic.7=4; thic.9=4
    styl.2=1; styl.3=1; styl.4=1; styl.5=1; styl.6=1; styl.7=1; styl.9=1;

    NInt_tc=4
    inten.1="2"; inten.2="3"; inten.3="4"; inten.4="5"; inten.5="6"; inten.6="7"; inten.7="9"
    cdn.1="(TD)"; cdn.2="(TS)"; cdn.3="(STS)"; cdn.4="(TY)"

    x00=1; x01=8; y00=1; y01=10
*   "set vpage "x00" "x01" "y00" "y01 
   
*   "sdfopen pca_ntc.nc"
   "sdfopen 4TC_paper.nc"
*   "open lsmask.ctl" 
   lat00=0; lat01=35; lon00=90; lon01=130; dlat=lat01-lat00; dlon=lon01-lon00
   "set lon "lon00" "lon01; "set lat "lat00" "lat01
*   "set map auto"; "draw map"
   "set ccolor 0"; "set clab off"; "set grads off"
   "set mpdset hires"
*   "d pc"
   "d tc"
*   "basemap L 7 1 M "
*   "basemap L 15 1 M"

*   dx=x01-x00; dy=y01-y00
*   xscal=dx/dlon; yscal=dy/dlat

*  Define something needed:
   'query gxinfo'
   rec3 = sublin(result,3) 
   rec4 = sublin(result,4) 
   xlim1 = subwrd(rec3,4)
   xlim2 = subwrd(rec3,6)
   ylim1 = subwrd(rec4,4)
   ylim2 = subwrd(rec4,6)
   xd = xlim2 - xlim1
   yd = ylim2 - ylim1
*   say 'xd = 'xd'  yd = 'yd
*   say 'xlim2 = 'xlim2'  ylim2 = 'ylim2
   'query dims'
   rec3 = sublin(result,2) 
   rec4 = sublin(result,3) 
   xl1 = subwrd(rec3,6)
   xl2 = subwrd(rec3,8)
      if (xl1 < 0)
      xl1 = xl1  + 360.0
      endif
      if (xl2 < 0)
      xl2 = xl2  + 360.0
      endif
   yl1 = subwrd(rec4,6)
   yl2 = subwrd(rec4,8)
   lonw = xl1
   lats = yl1
   xld = xl2 - xl1
   yld = yl2 - yl1
*   say 'xld = 'xld'  yld = 'yld
*   say 'lats = 'lats'  lonw = 'lonw

*  Legend first:
   x00=1; x01=8; dx=(x01-x00)/NInt_tc
   y00=1.2; icd=0
   While (icd<NInt_tc)
      icd=icd+1; x0=x00+(icd-1)*dx; xx0=x0+0.2
      CDtc=inten.icd
      draw_tit(x0+0.9, y00, inten.icd " "cdn.icd)
      "set line "colo.CDtc" 1 4"
      "draw line "x0" "y00-0.1" "x0+dx-0.4" "y00-0.1
      draw_tit((x00+x01)/2, y00-0.3,"TC intensity")
   Endwhile
   draw_hsts((x00+x01)/2, (ylim1+ylim2)/2-0.2,"Hoang Sa")
   draw_hsts((x00+x01)/2+0.2, (ylim1+ylim2)/2-1.5,"Truong Sa")

*   "drawrec 105 120 5 25 1"
*   "drawrec 110 120 5 15 5"
*   pull dd

********************************************************
*  Years to be plotted: Case=3
   yy0=1997; yy1=1997
   Nyear=yy1-yy0+1
*  ENSO years: Case=1 for El Nino, Case=2 for La Nina
   NEL=15; NLA=18
   E.1.1=1951; E.1.2=1957; E.1.3=1963; E.1.4=1965; E.1.5=1972; E.1.6=1977; E.1.7=1982; E.1.8=1987; 
   E.1.9=1991; E.1.10=1994; E.1.11=1997; E.1.12=2002; E.1.13=2004; E.1.14=2009; E.1.15=2015

   E.2.1=1954; E.2.2=1955; E.2.3=1956; E.2.4=1964; E.2.5=1970; E.2.6=1971; E.2.7=1973; E.2.8=1974; 
   E.2.9=1975; E.2.10=1984; E.2.11=1985; E.2.12=1988; E.2.13=1998; E.2.14=1999; E.2.15=2000; 
   E.2.16=2007; E.2.17=2010; E.2.18=2011

*  ------ Chose Case for plotting: -------
   Case=3
   if (Case=1); Ncase=NEL; ica=1; Fout="EL"; endif
   if (Case=2); Ncase=NLA; ica=2; Fout="LA"; endif
   if (Case=3); Ncase=Nyear; Fout=yy0"-"yy1; endif

*  Reading and plotting:
*   filename = TC_tr.txt

*   filename = "/home/tanpv/Progs/TCs/TC_Dat/TC_B12_for_GrADS.txt"
   OutDir="/home/tanpv/Progs/TCs/TC_Statics/PDF/"
 Nmon=12; imon=0
* While (imon<Nmon)
   imon=10
*   imon=imon+1
   if (imon<10); filename="/home/tanpv/Progs/TCs/TC_Dat/TC_tr_for_GrADS_0"imon".txt"; endif
   if (imon>=10); filename="/home/tanpv/Progs/TCs/TC_Dat/TC_tr_for_GrADS_"imon".txt"; endif
* --------------------------------------
   code=0
   res=read (filename)
   code=sublin(res,1)
*   if (code !=0); say Error !!!; pull dd; endif
   Ntc=sublin(res,2)

   Main_title ((x00+x01)/2, 9.7, "TC tracks of month "imon" for period "Fout)

   While (code=0)
      res=read (filename)
      code=sublin(res,1)
      if (code=2); "draw map"; "gxprint "OutDir"Tracks_"Fout"_"imon".pdf pdf"; "close 1"; break; endif
      if (code=9); say "something wrong !! "code; "close 1"; break; endif
      xx=sublin(res,2)
      Nlf=subwrd(xx,1)
      res=read (filename)
      code=sublin(res,1)
      xx = sublin(res,2)
      id = subwrd(xx,1); yr = subwrd(xx,2); mm = subwrd(xx,3); dd = subwrd(xx,4); hh = subwrd(xx,5);
      lon1 = subwrd(xx,6); lat1 = subwrd(xx,7); CD1 = subwrd(xx,8) 
      if (CD1=9); CD1=3; endif
      lon1 = lon1*0.10; lat1 = lat1*0.10
      x1 = (lon1 - lonw)/ xld * xd; y1 = (lat1 - lats)/ yld * yd
      x1 = x1 + xlim1; y1 = y1 + ylim1

      ilf=1
      While (ilf<Nlf)
         ilf=ilf+1
         res=read (filename)
         code=sublin(res,1)
         xx = sublin(res,2)
         id = subwrd(xx,1)
         yr = subwrd(xx,2); mm = subwrd(xx,3); dd = subwrd(xx,4); hh = subwrd(xx,5);
         lon2 = subwrd(xx,6); lat2 = subwrd(xx,7)
         lon2 = lon2/10; lat2 = lat2/10; CD2 = subwrd(xx,8); if (CD2=9); CD2=3; endif

         x2 = (lon2 - lonw)/ xld * xd; y2 = (lat2 - lats)/ yld * yd
         x2 = x2 + xlim1; y2 = y2 + ylim1

         yyyy=0
         While (yyyy<Ncase)
           yyyy=yyyy+1; 
           if (Case=3);
              if (yr>=yy0 & yr<=yy1);
*                 say yr" "CD2
                 "set line "colo.CD2" 1 4"
                 if (x1>xlim1 & x1<xlim2 & x2>xlim1 & x2<xlim2 & y1>ylim1 & y1<ylim2 & y2<ylim2);
                    "draw line "x1" "y1" "x2" "y2;
                 endif
                 x1 = x2; y1 = y2;
              endif
              yyyy=Ncase
           endif

           if (Case=1 | Case=2);
              if (yr=E.ica.yyyy);
                 say yyyy" "yr" "E.ica.yyyy" "CD2
                 "set line "colo.CD2" 1 4"
                 if (x1>xlim1 & x1<xlim2 & x2>xlim1 & x2<xlim2 & y1>ylim1 & y1<ylim2 & y2<ylim2);
                    "draw line "x1" "y1" "x2" "y2;
                 endif
                 x1 = x2; y1 = y2;
              endif
           endif 
         Endwhile
         if (code=2); "draw map"; "gxprint "OutDir"Tracks_"Fout"_"imon".pdf pdf"; "close 1"; break; endif
         if (code=9); say "something wrong !! "code; "close 1"; break; endif
      Endwhile
      "draw map";
   Endwhile
   pull dd
   "draw map";"gxprint "OutDir"Tracks_"Fout"_"imon".pdf pdf"
   say "next mon: "imon
   "c"
* Endwhile
* say "Done !!!"
*   pull dd
*  ----------------------------------------------------
*******************************************************
   function Main_title (x, y, title)
     "set string 1 c 5.5 0"
     "set strsiz 0.11 0.13"
     "draw string "x" "y" "title
   return
*******************************************************

   function draw_tit(x, y, title)
      "set string 1 c 5.5"; "set strsiz 0.12 0.12"
      "draw string "x" "y" "title
   return
*******************************************************

   function draw_hsts(x, y, title)
      "set string 1 c 4.0 45"; "set strsiz 0.10 0.10"
      "draw string "x" "y" "title
   return
*******************************************************