pro compare_profiles_sub,x1,v1i,x2,v2i,ivt,pct,uni_2,statopt,$ ; in
                            x1ave,x1std,x2ave,x2std,meandif,rmsdif ; out

;---------------------------------------------------------------------------
; Routine computes profile comparison statistics for 2 sets of coincident 
; profile data.  The input data must be coincident, that is, measurement 
; x1(i) corresponds to x2(i).
;
; Differences are X2 - X1
;
; Input:
;   x1.....profile data, instrument 1, fltarr(altitude,record)
;   v1i....vertical scale for x1, pressure or altitude
;   x2.....profile data, instrument 2, fltarr(altitude,record)
;   v2i....vertical scale for x2, pressure or altitude
;   ivt....character, set to 'p' if v1 & v2 are pressure, 'z' if altitude
;   pct....character, set to 'd' if differences are to be in data units,
;          otherwise they are in percent
;
;   uni_2..if dataset 2 is on a uniform vertical grid (and v2i = fltarr(n)), then uni_2 = 1
;          if dataset 2 is not on a uniform vertical grid (and v2i = fltarr(n,m)), then uni_2 = 2
;
;   statopt...(1) differences are (x2-x1) / mean([x1,x2]), (2) differences are 2*(x2-x1)/(x2+x1)
;
; Output:
;   x1ave.....average profile for instrument 1
;   x1std.....standard deviation of x1ave
;   x2ave.....average profile for instrument 2
;   x2std.....standard deviation of x2ave
;   meandif...profile of mean differences
;   rmsdif....profile of RMS differences
;
; Source: Mark Hervig, GATS Inc.
;----------------------------------------------------------------------------

   miss = 1e20

   v1 = v1i
   v2 = v2i

   if (ivt eq 'p') then v1 = alog10(v1i)
   if (ivt eq 'p') then v2 = alog10(v2i)

;- Loop over records, put data set 2 on data set 1's vertical grid
;  set extended range to missing

   nrec1 = n_elements( x1(0,*) )
   nrec2 = n_elements( x2(0,*) )

   x2g   = x1*0.

   for ir = 0,nrec1-1 do begin

;-   put data set 2 on data set 1's vertical grid

     ok = where(x2(*,ir) lt miss and x2(*,ir) gt 0.0, nok)

     if (uni_2 eq 1) then begin
     
      if (nok gt 1) then x2g(*,ir) = interpol(x2(ok,ir),v2(ok),v1)
      bd = where(v1 gt max(v2) or v1 lt min(v2),nbd)
      if (nbd gt 0) then x2g(bd,ir) = 0.
     endif
     
     if (uni_2 eq 2) then begin
     
      if (nok gt 1) then x2g(*,ir) = interpol(x2(ok,ir),v2(ok,ir),v1)
      bd = where(v1 gt max(v2(*,ir)) or v1 lt min(v2(*,ir)),nbd)
      if (nbd gt 0) then x2g(bd,ir) = 0.
     endif
     
   endfor

;- do the stats

   nmin = nrec1 * 0.5

   nl1 = n_elements(v1)

   x1ave = fltarr(nl1)
   x1std = fltarr(nl1)

   x2ave = fltarr(nl1)
   x2std = fltarr(nl1)

   meandif = fltarr(nl1)
   rmsdif  = fltarr(nl1)
   stddif  = fltarr(nl1)

   for i = 0,nl1-1 do begin    ; loop over levels in profile

     ok = where(x1(i,*)  gt 0.0 and  x1(i,*) lt miss and $
                x2g(i,*) gt 0.0 and x2g(i,*) lt miss, nok)

     if (nok gt nmin) then begin

      r1 = moment( x1(i,ok))
      r2 = moment(x2g(i,ok)) 

      x1ave(i) = r1(0)        ; mean value of x1
      x2ave(i) = r2(0)        ; mean value of x2

      x1std(i) = sqrt(r1(1))  ; std dev of x1
      x2std(i) = sqrt(r2(1))  ; std dev of x2

      if statopt eq 1 then begin
      
       xbar = ( x1ave(i)+x2ave(i) ) *0.5    ; average of means

       dx = x2g(i,ok) - x1(i,ok)

       r  = moment(dx)     

       meandif(i) = 100.*r(0)/xbar                  ; mean difference %
       if (pct eq 'd') then meandif(i) = r(0)       ; mean difference 

       r  = moment(dx^2.) 
       rmsdif(i) = 100.*sqrt(r(0))/xbar             ; RMS difference % 
       if (pct eq 'd') then rmsdif(i) = sqrt(r(0))  ; RMS difference 

       stddif(i) = sqrt(rmsdif(i)^2 - meandif(i)^2)  ; difference std. dev.

      endif
      
;-    new formulation

      if statopt eq 2 then begin

       dx = 200.*(x2g(i,ok) - x1(i,ok)) / (x2g(i,ok) + x1(i,ok))

       if (pct eq 'd') then dx = x2g(i,ok) - x1(i,ok)

       r = moment(dx)     

       meandif(i) = r(0)                ; mean diff
       rmsdif(i)  = sqrt(r(1) + r(0)^2) ; RMS diff
       stddif(i)  = sqrt(r(1))          ; std. dev. of the % diffs

      endif

     endif

   endfor

;- done

   return
   end