/*-------------------------------------------------------------------\
|  Thermal Index Calculator (ti)                       Version 1.70  |
|  written by Kevin Ford                               2-21-96       |
|  version 1 written Feb. 1992                                       |
|								     |
|  VAX/VMS Compatibility added 06/30/93 by Kim Reniska               |
|  Text graphics added 12/13/95 by Gary Helmstetter (GH)             |
|                                                                    |
|  Six modes of operation:				             |
|  Interactive:   						     |
|     ti                read from .uac file (text)                   |
|     ti -1		read from .uac file			     |
|     ti -3	        read from archived file			     |
|  Batch							     |
|     ti    filename station hightemp maxz  ; read from .uac (text)  |
|     ti -1 filename station hightemp maxz  ;  read from .uac file   |
|     ti -3 filename date hightemp maxz  ;  read from archive file   |
\-------------------------------------------------------------------*/

/* must use cc ti.c -lm to compile */

#include <stdio.h>
#include <math.h>
#include <ctype.h>
#include <string.h>
#ifdef VMS
#define O_RDONLY 0
#include <unixio.h>
#else
#include <fcntl.h>
#endif

#define t0  (273.15)      /* convert from Kelvin to Celcius */
#define RG (29.286)
#define GD (-0.00974)      /* dry adiabatic lapse rate, deg. C per meter */
#define DALR (-.0029687)   /*GH same rate, deg. C per foot */
#define GSCA (40)          /*GH width of graph in characters */
#define G0 (-0.0065)       /* standard lapse rate */
#define CFM  (0.3048)        /* convert from feet to meters GH:by multication */
#define spdcvt (1.9458)   /* convert from m/s to knots */

/* stored data structure for .uac files */

typedef struct {
  char ident[12], loc[12], loc2[12];
  long idcode;
  float lat,lon;
  long ml,sl,wl;   /* num. of mand. levels, sig. levels, wind levels */
  float mp[20], mh[20], mt[20], mdp[20], mdir[20], mspd[20];
  float sp[50], st[50], sdp[50];
  float wp[50], wdir[50], wspd[50];
} filerec;
filerec rec;

/* raw data from the file, for each level */

float p[100],           /* pressure in millibars */
      t[100],           /* temperature in degrees Kelvin */
      dp[100];          /* dewpoint in degrees Kelvin    */
int   ddir[100],        /* wind direction, degrees */
      dspd[100];        /* wind speed in knots */

/* auxilliary 1000 foot interval wind data from the file */

int winddir[160],       /* wind direction, degrees */
    windspd[160];       /* wind speed, knots       */

/* calculated data for each pressure level */

float h[100],           /* height of level in meters MSL  */
      tv[100];          /* virtual temperature in degrees Celcius  */

/* These arrays contain calculated data at 500 foot intervals */

float ti[300],          /* thermal index     */
      ta[300],		/*GH temp at altitude (replaces single variable) */
      alt[300];         /* altimeter reading */
int   dir[300],         /* wind direction, degrees */
      spd[300];         /* wind speed in knots */

float Gmin,Gmax,Glow,Gsca;/*GH lowest/highest tv - to get base & scale for: */
int   Gtmp,Glr;		/*GH Glow-based vtemp offsets, in scaled degrees C */
int   Ghh,Gdd,Gmm,Gyy;	/*GH to parse time&date from input filename */
char *Gptr,*Gptr1;	/*GH parse pointers */
char *Gmon[] = 		/*GH pretty date */
    {"JAN","FEB","MAR","APR","MAY","JUN","JUL","AUG","SEP","OCT","NOV","DEC"};
char  Gch;		/*GH char for drawing graph */
char  Graph[GSCA+2];	/*GH current line of graph */
char  Graphend=0;	/*GH null to end full graph lines */
float Gx;		/*GH scratch float */

int levels;             /* total number of pressure levels */

/* calculate the virtual temp given the temp, dewpoint and pressure */

float vtemp(t,dp,p)
  float t,dp,p;
{ float e;
  e=6.11 * exp(5423.0*(1.0/t0 - 1.0/dp));
  return(t/(1.0-.378*e/p));
}

/* altitude for a given pressure in the standard atmosphere */

float stdz(p)
  float p;
{
  return((-288.0/G0) * (1.0- pow(p/1013.25,-RG*G0)));
}

/* convert a string to an integer */

int strtoint(buf)
  char *buf;
{
  int d; char w;
  d=0;
  while ((w= *buf++)>0) d=10*d+w-'0';
  return(d);
}

/* decode a coded temperature/dewpoint from an text data file */

tempcvt(n,t,dp)
  int n;         /* coded form */
  float *t, *dp; /* decoded temp & dewpoint */
{
  int q; float h;
  q=n/100;
  h=((float) q)/10.0;
  if (q%2==0) *t=t0+h; else *t=t0-h;
  q=n%100;
  h=(float) q;
  if (q<56) *dp = *t - h/10.0;
    else if (q<99) *dp = *t - h+50.0;
      else *dp = t0 - 99.9;
}

/* virtual temperature at altitude a */

float tempalt(a)
  float a;
{
  int i;
  if (a<=h[0]) return(tv[0]);
  if (a >= h[levels-1]) return(tv[levels-1]);
  for (i=1; h[i]<a; i++);
  return(tv[i-1]+(tv[i]-tv[i-1])*(a-h[i-1])/(h[i]-h[i-1]));
}

main(argc, argv)
  int argc;
  char *argv[];
 {
  int idx,itemp,found,pr,p0,n,m,WND,hi,lo,ht,hv,i,j,
      fd, nid, stat, ti1, ti3, batch, trig[200], cb; /* fscanf() */
  long lseek(), read(); /* open() */
  FILE *fopen(),*fp;
  char buf[200],fname[100], id[10], date[20], date2[20],
     block[300][10], station[20];
  float salt,statalt,ftemp,htt,hd,k,z,a,maxz;

/*--------------------- Process command line arguments ----------------*/

  ti1 = 0; ti3=0; batch=0; i=1;
  if (argc>1) {
     if (strcmp(argv[1],"-3")==0) { ti3=1; i++; }
     if (strcmp(argv[1],"-1")==0) { ti1=1; i++; }
  }
  if (argc>i) {
    if (argc<i+4) { printf("ti error: Too few arguments\n"); exit(-1); }
    strncpy(fname, argv[i], 99); fname[99]=0; i++;
    if (ti3) { strncpy(date,argv[i],10); date[10]=0; }
    else { strncpy(station,argv[i],6); station[6]=0; }
    i++;
    sscanf(argv[i],"%d",&ht); i++;
    sscanf(argv[i],"%f",&maxz);
    batch=1;
  }

/*------------ Print headers and open file ------------- */

  if (!batch) {
    printf("Thermal Index Calculator\n");
    printf("     by Kevin Ford\n\n");
    printf("Reading files in the wxp ");
    if (ti3) printf("text archive format (one station per file)\n");
    if (ti1+ti3==0) printf("text .uac format (one time per file)\n");
    if (ti1) printf(".uac format (non-text)\n");
      printf("\nEnter filename : ");
      scanf("%s",fname);
  }
  if (!ti1) {
    if ((fp=fopen(fname,"r")) == NULL) {
        printf("ti error: cannot open file %s.\n",fname);
        printf("no data available.\n"); exit(-1);
    }
  }
  else if ((fd = open(fname, O_RDONLY,0)) == -1) {
     printf("ti error: cannot open file %s.\n",fname);
     printf("no data available.\n"); exit(-1); exit(-1);
  }

/* get month/day/year from filename */

    Gdd=0;
    if (!ti3) {
#ifdef VMS							/*GH*/
    if (Gptr=strchr(fname,']')) ++Gptr;				/*GH*/
	else Gptr=fname;					/*GH*/
    if (Gptr1=strchr(Gptr,'.')) Gptr= ++Gptr1;			/*GH*/
    sscanf(Gptr,"%4d%2d%2d%2d",&Gyy,&Gmm,&Gdd,&Ghh);		/*GH*/
/*  printf("yy=%4d mm=%3d dd=%3d hh=%3d\n",Gyy,Gmm,Gdd,Ghh);	/*GH debug*/
#else /* UNIX */
    for (Gptr=fname-1+strlen(fname); Gptr>=fname; Gptr--)
      if (*Gptr=='/') break;
    Gptr++;
    sscanf(Gptr,"%4d%2d%2d%2d",&Gyy,&Gmm,&Gdd,&Ghh);		/*GH*/
#endif /* def VMS */						/*GH*/
    } /* if not ti3 */

/*--------- Get date/station and read data for text data file --------------*/

  if (!ti1) {
    if (!batch) {
       if (ti3) { printf("Enter date (YYMMDDHH) : "); scanf("%s",date); }
       else { printf("Enter station : "); scanf("%s",station); }
    }
    found=0;
    while (!found) {
      if (fscanf(fp,"%s",buf)==EOF) found=2;
      if ((strncmp(buf,"DATE",4)==0) && (ti3)) {
        fscanf(fp,"%s",date2);
        if (strncmp(date2,date,8)==0) found=1;
      }
      if ((!ti3) && (strncmp(buf,station,3) == 0)) found=1;
    }
    if (found!=1) {
      if (ti3) printf("ti error: no data for %s.\n",date);
      else printf("ti error: no data for %s.\n",station);
      exit(-1);
    }
    if (ti3) fscanf(fp,"%s",station);
    for (i=0; block[i-1][0]!='$'; i++) fscanf(fp,"%s",block[i]);
    fclose(fp);
  }   /* if !ti1 */

/*--------- Get station id and read upper air data from .uac file -------*/

  else {
  lseek(fd,4284L,0); nid=0;
  while ((read(fd,&rec,sizeof(filerec))) == sizeof(filerec))
    strcpy(block[nid++],rec.ident);
  if (nid==0) {
     printf("ti error: file %s contains no sounding data.\n",fname);
     exit(-1);
  }

  if (!batch) {
    printf("Stations:\n\n");
    for (i=0; i<nid; i++) {
       block[i][4]=0;
       printf("%5s",block[i]); if (i%15==14) printf("\n");
    }
    for (stat= -1; stat<0;) {
      printf("\n\nEnter station id: "); scanf("%s",station);
      for (i=0; station[i]!=0; i++)
        if ((station[i]>='a')&&(station[i]<='z')) station[i] -= ('a'-'A');
      for (j=0; j<nid; j++) if (strcmp(station,block[j])==0) { stat=j;break; }
    }
  }
  else {
    stat= -1;
    for (i=0; station[i]!=0; i++)
      if ((station[i]>='a')&&(station[i]<='z')) station[i] -= ('a'-'A');
    for (j=0; j<nid; j++) if (strcmp(station,block[j])==0) { stat=j;break; }
    if (stat<0) {
      printf("ti error: station %s data not found.\n",station);
      exit(-1);
    }
  }

  lseek(fd,4284L+stat*sizeof(filerec),0);
  read(fd,&rec,sizeof(filerec));
  close(fd);
  }  /* if ti1 */

/*------------ Read in station altitude from data file ------------------*/

  statalt = 0.0;
  if ((fp = fopen("stations","r")) != NULL) {
    while (fscanf(fp,"%s %f %d\n",id,&salt,&k) != EOF)
      if (strcmp(id,station)==0) { statalt = salt; break; }
    fclose(fp);
  }

/*--------------- Parse data from text data file ---------------------*/

  if (!ti1) {

              /* read in significant levels */

  p0= -999;
  levels=0; j=53;
  while ((block[j][0]!='X')&&(block[j][0]!='-')) {
    pr = strtoint(block[j]);  if (block[j][0]=='0') pr+=1000;
    if (j==53) p0=pr;  /* sfc pressure */
    p[levels] = (float) pr;
    dspd[levels] = -999; ddir[levels] = -999;
    j++;
    if (block[j][0]!='X') {
       tempcvt(strtoint(block[j]),&t[levels],&dp[levels]); j++;
    }
    else { t[levels]= -999.0; dp[levels] = -999.0; }
    levels++;
  }
  j++; while (block[j][0]=='X') j++;
  idx=j;
  if (levels<1) {printf("ti error: data for %s not found.\n",station);
                 exit(-1); }

            /* read in mandatory levels */

  WND=0;
  for (j=0; j<53; j+=3) {
    if (block[j][0]!='X') {
      n = strtoint(block[j]);
      pr = n/1000; pr*=10; if (pr==0) pr=1000;
      if (pr==880) pr=n%1000;                   /* level of tropopause */
      else if (pr==770) { pr=n%1000; WND=1; }   /* level of max. wind  */
      if (pr<p0) {
        p[levels] = (float) pr;
        if (!WND) {
          if (block[j+1][0]!='X')
	    tempcvt(strtoint(block[j+1]),&t[levels],&dp[levels]);
          else { t[levels]= -999.0; dp[levels] = -999.0; }
        }
	else {t[levels]= -999.0; dp[levels] = -999.0; }
        if (block[j+2-WND][0]!='X') {
          n=strtoint(block[j+2-WND]);
          ddir[levels] = (n/500)*5;
          dspd[levels] = n%500;
        }
        else {
          ddir[levels] = -999; dspd[levels] = -999;
        }
        levels++;
      }
    }
  }

               /* read in winds aloft data */

  j=idx;
  for (i=0; i<=150; i++) { windspd[i] = -999; winddir[i] = -999; }
  while (block[j][0]!='$') {
    n=strtoint(block[j]);
    if ((block[j+1][0]!='X') && (block[j+1][0]!='$')) {
      m=strtoint(block[j+1]);
      windspd[n] = m%500;
      winddir[n] = (m/500)*5;
    }
    j += 2;
    if (block[j-1][0]=='$') j--;
  }

  } /* if !ti1 */

/*------------------ Organize data from .uac file --------------------*/

  else {

  for (i=0;i<rec.sl;i++) {
    p[i]=rec.sp[i];
    t[i]=rec.st[i] + t0;
    dp[i]=rec.sdp[i] + t0;
    if (dp[i] < 100.0) dp[i] = t0 - 99.9;
    dspd[i]= -999;
    ddir[i]= -999;
  }
  for (i=0,j=rec.sl;i<rec.ml;i++,j++) {
    p[j]=rec.mp[i];
    t[j]=rec.mt[i] + t0;
    dp[j]=rec.mdp[i] + t0;
    if (dp[j] < 100.0) dp[j] = t0 - 99.9;
    dspd[j]= (int) (rec.mspd[i] * spdcvt + 0.5);
    ddir[j]= ((int) ((rec.mdir[i] + 2.5)/5.0)) * 5;
  }
  for (i=0;i<100;i++) { windspd[i]= -999.0; winddir[i]= -999.0; }
  for (i=0; i<rec.wl; i++) {
    j = (int) ((rec.wp[i]/CFM + 500.0)/1000.0);
    windspd[j] = (int) (rec.wspd[i] * spdcvt + 0.5);
    winddir[j] = ((int) ((rec.wdir[i] + 2.5)/5.0)) * 5;
  }
  levels=rec.ml+rec.sl;
  }  /* if ti1 */

/*-------------------- Sort data by pressure ------------------------*/

/* printf("raw data:\n");
   for (i=0; i<levels;i++) printf("%8.1f%8.1f%8.1f\n",p[i],t[i],dp[i]);
*/

  for (i=0;i<levels-1;i++) for (j=i+1;j<levels;j++)
    if ((p[i]<p[j])||((p[i]==p[j])&&(dspd[i]>dspd[j]))) {
        ftemp=p[i];p[i]=p[j];p[j]=ftemp;
        ftemp=t[i];t[i]=t[j];t[j]=ftemp;
        ftemp=dp[i];dp[i]=dp[j];dp[j]=ftemp;
        itemp=dspd[i];dspd[i]=dspd[j];dspd[j]=itemp;
        itemp=ddir[i];ddir[i]=ddir[j];ddir[j]=itemp;
    }

  p[levels]=9999.0;
  for (i=0; i<levels; i++)
    if ((t[i]<0.0)||(p[i]<9.9)||(p[i]-p[i+1]<.1)) {
      levels--;
      for (j=i; j<=levels; j++) {
        p[j]=p[j+1]; t[j]=t[j+1]; dp[j]=dp[j+1]; dspd[j]=dspd[j+1];
        ddir[j]=ddir[j+1];
      }
      i--;
    }

/* printf("\sorted data:\n");
   for (i=0; i<levels;i++) printf("%8.1f%8.1f%8.1f%8d%8d\n",
   p[i],t[i],dp[i], dspd[i],ddir[i]);  */

  if ((levels<2) || ((ti1) && (rec.mh[0] < -100.0))) {
    printf("ti error: insufficient upper air data for ");
    if (ti3) printf("%s\n",date);
    else printf("%s\n",station);
    exit(-1);
  }

/*------------------ calculate heights of levels ------------------------*/

   /* get surface elevation */

  h[0]=statalt;

/*  printf("station elevation: %f\n",statalt);
  if (ti1) printf("lowest ml altitude: %f\n",rec.mh[0]);  */

  if (ti1) if ((h[0]-rec.mh[0])*(h[0]-rec.mh[0]) > 10.0) {
     printf("ti error: surface data missing\n");
     exit(-1);
  }

  for (i=0; i<levels;i++) tv[i] = vtemp(t[i],dp[i],p[i]);
  for (i=1; i<levels;i++) {
    ftemp = log(p[i-1]) - log(p[i]);
    if (tv[i]==tv[i-1])  h[i]=h[i-1]+RG*tv[i]*ftemp;
    else h[i]=h[i-1]+RG*(tv[i]-tv[i-1]) * ftemp/(log(tv[i])-log(tv[i-1]));
  }

/*--------------------- Get forecast hi, dp, max alt to calc TI ---------*/

    /* set dewpoint = average dewpoint in lowest 500 meters */

  z = 0.0; i=1; while (h[i] < h[0] + 500.0) {
    z += ((h[i]-h[i-1])*(dp[i]+dp[i-1])/2.0);
    i++;
  }
  if (i >= levels) hd=dp[0];
  else {
    ftemp = h[i] - h[i-1];
    a = dp[i-1]*(h[i]-(h[0]+500.0))/ftemp + dp[i]*(500.0+h[0]-h[i-1])/ftemp;
    z += (500.0+h[0]-h[i-1]) * (a + dp[i-1]) / 2.0;
    hd = z / 500.0;
  }

  if (!batch) {
    printf("\nEnter forecast high temp (F) : "); scanf("%d",&ht);
    printf("Enter maximum altitude to calculate TI : "); scanf("%f",&maxz);
  }
  ftemp = (float) (ht-32.0)/1.8+t0;
  if (ftemp<t[0])
    { ftemp=t[0]; ht = (int) ((ftemp-t0)*1.8+32.5); }
  if (ftemp > 333.0) ftemp=t[0];
  if (maxz*CFM<h[0]+305.0) maxz=(h[0]+305.0)/CFM;
  if (hd<200.0) hd = t[0]-5.0;
  if (hd>ftemp) hd=ftemp;
  cb = (int) ((ftemp - hd)*4.0 + .5) * 100;  /* est. cloud base */
  htt = vtemp(ftemp,hd,p[0]);

  if (Gdd)
    printf("  %2d-%3s-%4d %2d UTC  ", Gdd, Gmon[Gmm-1], Gyy, Ghh);
  else
    printf("           ");
  printf("Soaring report from %s upper air data.\n", station);
  if (ht < 150) {
    printf("   Forecast high:%3d F; est. base of any clouds:%5d feet AGL.",
           ht,cb);
/*    printf("Surface virtual temp: %4.1f C\n", htt-t0+.05); */
  }
  else
    printf("No forecast high or estimated cloud base.");
  printf("\n\n");


  for (i=0; i<199; i++) {dir[i]= -999.0; spd[i]= -999.0; }

/*------------------ calculate thermal index ---------------------------*/

  lo=(int) (h[0]/CFM/500.0) + 1;
  hi=(int) ((maxz+250.0)/500.0);
  if (hi>120) hi=120;  /* don't go over 60,000 feet MSL */
  if (hi<lo+9) hi=lo+9;   /* do at least 5000 feet AGL */

  for (i=0; i<levels; i++)
    if ((ddir[i]>=0) && (dspd[i]>=0)) {
        j = (int) ((h[i]/CFM+250.0)/500.0);
        dir[j] = ddir[i]; spd[j] = dspd[i];
  }

  for (Gmin=999.0, Gmax=0.0, i=lo; i<=hi; i++) {
    k=500.0 * CFM * i;
    if ((i%2==0)&&(winddir[i/2]>=0)) {
      dir[i]=winddir[i/2]; spd[i]=windspd[i/2];
    }
    itemp = ht;
    if (ht >= 150) itemp = (int) ((t[0]-t0)*1.8+32.5);
    for (j=1; j<levels;j++) if ((h[j-1]<k)&&(k<=h[j])) {
        ta[i] = tv[j-1] + (tv[j]-tv[j-1]) * (k-h[j-1])/ (h[j]-h[j-1]);
        ti[i] = ta[i] - (htt + GD * (k-h[0]));
        trig[i] = itemp + (int) (ti[i]*1.8+.5);
        break;
    }
    if (ta[i] < Gmin) Gmin = ta[i];
    if (ta[i] > Gmax) Gmax = ta[i];
  }

/*GH Gsca = chars per deg. K of extremes of overlapped (actual,DALR) ranges. */
/*GH Glow = lowest temp in ranges in deg. K, for setting min point of graph. */

  Glow = ta[hi] - ti[hi];		/*GH min of DALR range */
  if (Glow > Gmin) Glow = Gmin;		/*GH min of actual range */
  Gx = ta[lo] - ti[lo];			/*GH max of DALR range */
  if (Gx < Gmax) Gx = Gmax;		/*GH max of actual range */
  Gsca = (float) GSCA / (Gx - Glow);

/*  printf(
    "Gmax %5.1f  Gmin %5.1f   hiLR %5.1f @%5d  loLR %5.1f @%5d\n",
    Gmax, Gmin,
    ta[lo],				lo*500,
    ta[lo] + (float)(500*DALR*(hi-lo)),	hi*500); GH debug*/
/*  printf("Glow %5.1f K %6.1f F  Range %5.1f C  Gsca %5.3f\n",
    Glow, ((Glow-t0)*1.8)+32., Gx - Glow, Gsca); GH debug*/


/*---------------------- print out tables ----------------------------*/

  printf("                    === Raw Upper-Air Data ===");
  j = (levels < 10) ? levels : 10;
  printf("\nFeet MSL:");
  for (i=0; i < j; i++) printf("%7d",	(int) (h[i]/CFM+.5)	);
  printf("\nPres mb: ");
  for (i=0; i < j; i++) printf("%7.0f",	p[i]			);
  printf("\nTemp C:  ");
  for (i=0; i < j; i++) printf("%7.1f",	t[i]-t0+.01		);
  printf("\nVirT C:  ");
  for (i=0; i < j; i++) printf("%7.1f",	tv[i]-t0+.01		);
  printf("\nDewPt C: ");
  for (i=0; i < j; i++)
    if (dp[i] > 180.0)	printf("%7.1f", dp[i]-t0+.01		);
    else		printf("       ");
  printf("\nWdir@kts:");
  for (i=0; i < j; i++)
    if (dspd[i] >= 0)	printf("%4d%3d", ddir[i], dspd[i]	);
    else		printf("       ");

  if (p[0]<500.0) {
     printf("Insufficient upper air data for thermal index calculations.\n");
     exit(-1);
  }

  for (i=0; i<=GSCA; i++) Graph[i] = '-'; Graph[GSCA+1]=0;
  printf("\n\n       ");
  printf("=== Interpolations (temps in deg. F, altitudes in feet MSL) ===\n");
  printf("  MSL  *TI* Wdir@kts trig ");
  printf(" VirT  %2.1f degrees/division (\"`\": Dry Adiabatic)\n",1.8/Gsca);
  printf("-----  ---- -------- ---- . ---- %s\n", Graph);

  for (i=hi; i>=lo; i--) {
    printf("%5d", i*500);
    if (ht<150) 		printf("%6.1f", ti[i]+.05);
      else			printf("      ");
    if ((spd[i]<0)||(dir[i]<0)) printf("         ");
      else			printf("%5d%4d",dir[i],spd[i]);
    printf("%5d |%5.1f ", trig[i], ((ta[i]-t0)*1.8)+32.0);

    for (j=0; j<GSCA; j++) Graph[j]=' ';
    Glr =  ( ta[i]-ti[i] - Glow) * Gsca;
    if ((Glr >=0) && (Glr <= GSCA))
      Graph[Glr] = '`';
    Gtmp = ( ta[i] - Glow) * Gsca;
/*  printf("\nGlr:%5d   Gtmp:%5d          ",Glr,Gtmp); GH debug*/
    Gx = ta[i+1] - ta [i];
/*    if (Gx < -1.5)	Gch = '\\';
     else if (Gx > 1.5)	Gch = '/';
     else		Gch = ':'; */
    Gch=':';
    if ((Gtmp >=0) && (Gtmp <= GSCA))
      Graph[Gtmp] = Gch;
    if (Glr > Gtmp) Gtmp=Glr;
    Graph[(int) ++Gtmp] = '\0';
    printf("%s\n",Graph);
  }
  printf("\n\n   === Complete Upper Air Data ===\n\n");
  printf("   P(mb)     H(ft)    Tv(C)    T(C)    DP(C)   ");
  printf("wind dir  wind spd\n\n");
  for (i=0; i<levels;i++) {
    printf(" %8.1f%9d%9.1f%9.1f",p[i],(int) (h[i]/CFM+.5), tv[i]-t0+.01,
       t[i]-t0+.01);
    if (dp[i]>180.0) printf("%9.1f",dp[i]-t0+.01);
        else printf("         ");
    if (dspd[i]>=0) printf("%9d%9d",ddir[i], dspd[i]);
    printf("\n");
  }
  printf("Additional wind data:\n\n");
  printf("  H(ft)   dir   spd\n");
  printf(" ------   ---   ---\n");
  for (i=1; i<= 150; i++) if (windspd[i]>0) 
    printf(" %2d000    %3d   %3d\n",i,winddir[i],windspd[i]);
  printf("\n");
}