initial added

CONSTANTS.pro

+  C_NO_DATA   =  -1.
+
+  C_WATER     =   0.
+  C_MIXED     =  50.
+  C_ICE       = 100.
+
+  C_CLOUDFREE =   1
+  C_CLOUDY    =   2

LIST_COLORS.pro

+
+
+  colList=['Green',      'Red',         'Blue',           'Cyan',   'Orange',$
+           'Magenta',    'Dark Green',  'Cornflowerblue', 'Brown',  'GreenYellow'     ,$ 
+           'Dark Orchid','Tan',         'LightSeaGreen',  'Violet', 'Khaki',$
+           'LightCoral', 'Hot Pink',    'Yellow',     'Olive Drab', 'Dark Gray', 'Black']
+
+
+; see fsc_color
+;                Almond   Antique White      Aquamarine           Beige          Bisque           Black
+;                  Blue     Blue Violet           Brown       Burlywood        Charcoal      Chartreuse
+;             Chocolate           Coral        Cornsilk            Cyan  Dark Goldenrod       Dark Gray
+;            Dark Green      Dark Khaki     Dark Orchid     Dark Salmon       Deep Pink     Dodger Blue
+;             Firebrick    Forest Green            Gold       Goldenrod            Gray           Green
+;          Green Yellow        Honeydew        Hot Pink      Indian Red           Ivory           Khaki
+;              Lavender      Lawn Green     Light Coral      Light Cyan      Light Gray    Light Salmon
+;          Light Yellow      Lime Green           Linen         Magenta          Maroon     Medium Gray
+;         Medium Orchid        Moccasin            Navy           Olive      Olive Drab          Orange
+;            Orange Red          Orchid  Pale Goldenrod      Pale Green                            Peru
+;                  Pink            Plum     Powder Blue          Purple             Red            Rose
+;            Rosy Brown      Royal Blue    Saddle Brown          Salmon     Sandy Brown       Sea Green
+;              Seashell          Sienna        Sky Blue      Slate Gray            Snow    Spring Green
+;            Steel Blue             Tan         Thistle          Tomato       Turquoise          Violet
+;            Violet Red           Wheat           White          Yellow
+
+
+; not enough contrast
+; Papaya
+; 'Indian Red' and 'Brown'

LIST_GROUPS.pro

+
+  ; group list
+  grpList =['CMS','EUM','OCA','MPF','FUB','DLR','MFR','RMB','AWG','UKM','GSF','LAR','COX','OCA2','KNM']
+  ;grpList  =['CMS','EUM','OCA','MPF'      ,'DLR','MFR','RMB','AWG','UKM','GSF','LAR','COX','OCA2']
+
+  ; cannot be sorted in alphabetic order as guteInd depent on the order of the groups 
+  ; guteInd should be replaced by an automatic check
+  ; grpList = ['AWG','CMS','DLR','EUM','FUB','GSF','LAR','MFR','MPF','OCA','RMB','UKM'] ;  groups
+
+  grpLong = ['Climate Monitoring - SAF', $
+             'Eumetsat M?', $
+             'Eumetsat OCA', $
+             'Eumetsat MPEF', $
+             'Free University Berlin', $
+             'German Aerospace Center, APICS', $
+             'Meteo France', $
+             'Belgian Institute for Space Aeronomy', $
+             'Algorithm Working Group, University Madison/Wisconsin', $
+             'UK Met-Office', $
+             'NASA, Goddard Space Flight Center', $
+             'NASA, Langley Research Center', $
+             'German Aerospace Center, COCS', $
+             'Eumetsat OCA second cloud layer' $
+             , 'KNMI'] ;  groups
+
+
+;  grpLong = ['Algorithm Working Group, University Madison/Wisconsin', $
+;             'Climate Monitoring - SAF', $
+;             'German Aerospace Center', $
+;             'Eumetsat Mef', $
+;             'Free University Berlin', $
+;             'NASA, Goddard Space Flight Center', $
+;             'NASA, Langley Research Center', $
+;             'Meteo France', $
+;             'Eumetat MPF', $
+;             'Eumetsat OCA', $
+;             'Belgian Institute for Space Aeronomy', $
+;             'UK Met-Office'] ;  groups
+
+  ngrp = n_elements(grpList)
+  if verbose then print, '... defined ',ngrp,' data groups'
+
+  if (n_elements(grpList) ne n_elements(grpLong)) then print, '*** Warning: group-list and group-list-long dont have the same size.'
+
+  ; availability of products by the groups is defined with:
+  ; guteInd, see cws_stats_params.pro

LIST_PRODUCTS.pro

+  ; product listm, short string	
+  prdList     = ['rgb','ovw',$
+                  'cmb',  'ctp',  'cth',  'ctt',  'cod',  'ref',  'cph',  'lwp', $
+                  'cmb2', 'ctp2', 'cth2', 'ctt2', 'cod2', 'ref2', 'cph2', 'lwp2',$  ; second cloud layer (at the moment only OCA)
+                 'ucmb', 'uctp', 'ucth', 'uctt', 'ucod', 'uref', 'ucph', 'ulwp', $  ; uncertainty
+                 'ucmb2','uctp2','ucth2','uctt2','ucod2','uref2','ucph2','ulwp2']   ; uncertainty, second cloud layer (at the moment only OCA)
+              ;,'iwp','cty','cfr' $   ; next params exclusively DWD products!
+              ;  ,'hbl','hhl','hml','hto','hts','rgb','lst','ovw']
+
+  ; product listm, long string	
+  prdListLong = ['False Color Image', 'Atrain Overpass'$
+                , 'Cloud Mask','Cloud Top Pressure','Cloud Top Height','Cloud Top Temperature' $
+                ,'Cloud Optical Depth','Effective Radius','Cloud Phase','Cloud Water Path' $
+                ,'Cloud Mask (2)','Cloud Top Pressure (2)','Cloud Top Height (2)','Cloud Top Temperature (2)' $
+                ,'Cloud Optical Depth (2)','Effective Radius (2)','Cloud Phase (2)','Cloud Water Path (2)' $
+                ,'Uncertainty of Cloud Mask','Uncertainty of Cloud Top Pressure'$
+                ,'Uncertainty of Cloud Top Height','Uncertainty of Cloud Top Temperature' $
+                ,'Uncertainty of Cloud Optical Depth','Uncertainty of Effective Radius'$
+                ,'Uncertainty of Cloud Phase','Uncertainty of Cloud Water Path' $
+                ,'Uncertainty of Cloud Mask (2)','Uncertainty of Cloud Top Pressure (2)'$
+                ,'Uncertainty of Cloud Top Height (2)','Uncertainty of Cloud Top Temperature (2)' $
+                ,'Uncertainty of Cloud Optical Depth (2)','Uncertainty of Effective Radius (2)'$
+                ,'Uncertainty of Cloud Phase (2)','Uncertainty of Cloud Water Path (2)']
+                ;,'Ice Water Path','Cloud Type','Cloud Fraction' $   
+                ; next params exclusively DWD products!
+               ; ,'Humidity Boundary layer','hhl','hml','hto','hts','rgb','Land Surface Temperature','ovw']
+
+  nprd = n_elements(prdList)
+  if verbose then print, '... defined ',nprd,' product types'
+
+  if (n_elements(prdList) ne n_elements(prdListLong)) then print, '*** Warning: Product-list and Product-list-long dont have the same size.'
+
+  ; availability of products by the groups is defined with:
+  ; guteInd, see cws_stats_params.pro

LIST_REGIONS.pro

+; this is a list of predefined regions
+; the coordinates along x and y are in MSG pixels
+; it is possible to define more regions
+; or to use the keywords xrange and yrange, when you define 
+; your objects: cws_read or cws_stats
+
+case cut of
+   0: begin & xrange = [   0,3711] & yrange = [   0,3711] & scale = n_elements(scale) gt 0 ? scale : 0.2 & end
+   1: begin & xrange = [ 800,1200] & yrange = [2800,3200] & end
+   2: begin & xrange = [1980,2180] & yrange = [3300,3420] & end ; cpr 11317 : 1230
+   3: begin & xrange = [1600,1800] & yrange = [2060,2200] & end ; cpr 11318 1345
+   4: begin & xrange = [1600,1800] & yrange = [1800,2020] & end ; cpr 11318 1345
+   5: begin & xrange = [1888,2388] & yrange = [ 500,1000] & end ; cpr 11318 1345
+   6: begin & xrange = [1000,1500] & yrange = [3000,3400] & end ; cpr 11318 1400
+;  7: BEGIN & xrange = [3150,3600] & yrange = [2200,2700] & end ; arabian peninsula and see, cmb issue
+   7: BEGIN & xrange = [2950,3600] & yrange = [2100,2800] & end ; arabian peninsula and see, cmb issue
+   8: BEGIN & xrange = [1340,1800] & yrange = [2100,2550] & end ; west africa cmb
+   9: BEGIN & xrange = [1400,1700] & yrange = [1200,1800] & end ; cyclone
+  10: BEGIN & xrange = [2450,2749] & yrange = [1500,1799] & end ; tropical convection over land
+  11: BEGIN & xrange = [1016,1109] & yrange = [2348,2469] & end ; multi-layer
+  12: BEGIN & xrange = [2080,2110] & yrange = [2630,2660] & end ; Niamey, Nigeria
+  13: BEGIN & xrange = [2000,3110] & yrange = [2030,3660] & end ; Arabien peninsula 
+  14: begin & xrange = [   0,3711] & yrange = [   0,1855] & end ; southern hemisphere
+  15: begin & xrange = [2650,3250] & yrange = [ 300, 900] & end ; south-east begin of track 11317
+  16: begin & xrange = [1250,2120] & yrange = [1840,2640] & end ; central africa, thin cirrus
+  17: begin & xrange = [1688,2488] & yrange = [ 900,1900] & end ; marine stratocumulus, south west africa
+  18: begin & xrange = [ 450,1550] & yrange = [1800,2700] & end ; sunglint for track 11319, 2008-06-13 15:30UTC
+  19: begin & xrange = [1790,2070] & yrange = [ 943,1200] & end ; marine stratocumulus, south west africa
+                                                                ; for fellowship report, time 20080613_1345
+  20: begin & xrange = [1935,1962] & yrange = [ 945, 970] & end ; detail 
+                                  ; homogenious marine stratocumulus, south west africa
+                                  ; with CPR and CALIOP measurement                          
+                                  ; for fellowship report, time 20080613_1345
+  21: begin & xrange = [1906,1927] & yrange = [1066,1088] & end ; detail 
+                                  ; homogenious marine stratocumulus, south west africa
+                                  ; only CALIOP measurement
+                                  ; for fellowship report, time 20080613_1345
+  else: begin
+    print, '*** Error in LIST_REGIONS.pro'
+    print, '    unknown cut: ', cut
+    stop
+  end
+endcase

LIST_SENSORS.pro

+  ; sensor list
+  sensorList = ['MSG','AVHRR','MODIS','MERIS']
+
+  nsensor = n_elements(sensorList)
+  
+  if verbose then print, '... defined ', nsensor,' satellite sensors'

advanced_keyword_set.pro

+;---------------------------------------------------------------------------------------------
+;advanced_keyword_set(keyword) ermittelt ob keyword bergeben wurde so wie auch keyword_set()
+;allerdings mit dem vorteil, dass auch die moeglichkeit besteht eine null als keyword zu
+;uebergeben.
+;
+;Author: Max Reuter, 2000
+;
+;---------------------------------------------------------------------------------------------
+function advanced_keyword_set, keyword
+	temp = size(keyword)
+	temp = temp[1]
+	return, (temp gt 0)
+end

advanced_ptr_free.pro

+pro advanced_ptr_free, p, parents = parents, circle = circle
+	type = size(p, /type)
+		case type of
+			8	:	begin								;struct
+						for j = 0l, n_tags(p) - 1 do begin
+							advanced_ptr_free, p.(j), parents = keyword_set(circle) ? parents : 0, circle = circle
+						endfor
+					end
+			10	:	begin								;pointer
+						anz = ulong(n_elements(p))
+						for i = 0l, anz - 1 do begin
+							if ptr_valid(p[i]) then begin
+								if keyword_set(parents) then begin
+									if not is_element(p[i], parents) then begin
+										advanced_ptr_free, *(p[i]), parents = keyword_set(circle) ? [parents, p[i]] : 0, circle = circle
+										ptr_free, p[i]
+									endif
+								endif else begin
+									advanced_ptr_free, *(p[i]), parents = keyword_set(circle) ? p[i] : 0, circle = circle
+									ptr_free, p[i]
+								endelse
+							endif
+						endfor
+					end
+			11	:	begin								;object
+						print, 'caution: object could contain pointers, try to destroy object...'
+						obj_destroy, p
+					end
+			else:
+		endcase
+end
+
+;a = {q:findgen(10000), w:ptr_new(  {q:findgen(10000),w:ptr_new() }  ) }
+;(*a.w).w = ptr_new(a)
+;advanced_ptr_free, a

all_average.pro

+pro all_average
+
+  o=obj_new('cws_stats',scale=1.0)
+
+; arr_pro = ['cmb', 'cph', 'ctt', 'cth', 'ctp', 'cod', 'ref', 'lwp']   ; OK
+;  arr_pro = ['cmb', 'cph', 'ctt', 'cth', 'ctp', 'cod', 'ref', 'lwp']
+  arr_pro = ['cth']
+  npro = n_elements(arr_pro)
+
+  day   = [13,17,18,22,03]
+  month = [ 6, 6, 6, 6, 7]
+
+  ; for d = 0, 4 do begin 
+  for d = 0, 0 do begin 
+  ;  for h=0,23 do begin 
+  for h=12,15 do begin 
+      for m=00,45,15 do begin 
+        print, '=== TIME: ', h, m
+        o->set_date,2008,6,13,h,m
+        for p = 0, npro-1 do begin 
+          product = arr_pro[p]
+          o->set_product, product 
+        ;  avg=o->algorithm_average(/calc,/save,/log,/no_plot) ; 
+          avg=o->algorithm_average(/calc,/save,/no_plot) ; 
+        ;  std=o->algorithm_stddev(/calc,/save,/jpg)
+        ;  std=o->algorithm_stddev(/calc,/save,/jpg,/rel)
+        endfor
+      endfor
+    endfor
+  endfor 
+
+end

all_crossscatter.pro

+pro all_crossscatter
+
+  o=obj_new('cws_stats') ; ,scale=1.0
+
+;  arr_pro = ['cmb', 'cph'] ; does not work 
+;  arr_pro = ['ctt', 'cth', 'ctp', 'cod', 'ref', 'lwp']   ; OK
+
+  arr_pro = ['cod']
+  npro = n_elements(arr_pro)
+ 
+  h=12
+  m=0
+
+  for h=12,14 do begin
+    for m=00,45,15 do begin 
+      print, 'time: ', h, m
+      o->set_date,2008,6,13,h,m
+      for p = 0, npro-1 do begin 
+        product = arr_pro[p]
+        o->set_product, product
+;        o->make_corSingle                  
+;        o->make_corSingle, phase='ice'
+;        o->make_corSingle, phase='water'
+;        o->make_corSingle, surf='land'
+;        o->make_corSingle, surf='sea'      
+        o->make_corthumb                   ;;
+      endfor
+    endfor
+  endfor
+
+end

all_hist1d.pro

+pro all_hist1d
+
+  o=obj_new('cws_stats') ;,scale=1.0
+
+;  arr_pro = ['cmb', 'cph'] ; does not work 
+  arr_pro = ['ctt', 'cth', 'ctp', 'cod', 'ref', 'lwp']   ; OK
+;  arr_pro = ['lwp']
+  npro = n_elements(arr_pro)
+; mask='IND'
+
+  h=12
+  m=0
+
+  for h=12,14 do begin
+    for m=0,45,15 do begin 
+      print, 'time: ', h, m
+      o->set_date,2008,6,13,h,m
+      for p = 0, npro-1 do begin 
+        product = arr_pro[p]
+        o->set_product, product
+        o->hist1d, mask=mask
+        o->hist1d, phase='ice', mask=mask
+        o->hist1d, phase='water', mask=mask
+        o->hist1d, surf='land', mask=mask
+        o->hist1d, surf='sea', mask=mask
+      endfor 
+    endfor 
+  endfor
+
+end

all_make_image.pro

+pro all_make_image, phase=phase
+
+  o=obj_new('cws_stats') ; ,scale=1.0
+;  arr_pro = [] ; does not work 
+;  arr_pro = ['cmb', 'cph', 'ctt', 'cth', 'ctp', 'cod', 'ref', 'lwp']   ; OK
+  arr_pro = ['ctp']
+  npro = n_elements(arr_pro)
+; single=1
+  o->set_group, 'GSF'
+
+  h=12
+  m=0
+  day   = [13,17,18,22,03]
+  month = [ 6, 6, 6, 6, 7]
+; phase='water'
+
+  for d = 2, 2 do begin 
+;    for h=12,15 do begin
+    for h=0,2 do begin
+      for m=00,45,15 do begin 
+;      for m=00,00,15 do begin 
+        print, '=== TIME:', h, m
+        o->set_date,2008,month[d],day[d],h,m
+        for p = 0, npro-1 do begin 
+          product = arr_pro[p]
+          o->set_product, product  
+          ;o->make_image, product=product, /single
+          o->make_image, product=product
+        endfor
+      endfor
+    endfor
+ endfor
+
+end

all_rgb.pro

+pro all_rgb
+
+  o=obj_new('cws_stats')
+
+; arr_pro = ['cmb', 'cph', 'ctt', 'cth', 'ctp', 'cod', 'ref', 'lwp']   ; OK
+  arr_pro = ['rgb']
+  npro = n_elements(arr_pro)
+
+  for h=0,23 do begin 
+    for m=00,45,15 do begin 
+      print, '=== TIME: ', h, m
+      o->set_date,2008,6,13,h,m
+      for p = 0, npro-1 do begin 
+        product = arr_pro[p]
+        o->set_product, product 
+        o->make_image
+      endfor
+    endfor
+  endfor
+
+end

all_transform.pro

+pro all_transform, group=group
+
+  default, group, 'GSF'
+  o=obj_new('cws_read',scale=1.0)
+  o->set_group, group
+
+; arr_pro = ['cmb', 'cph', 'ctt', 'cth', 'ctp', 'cod', 'ref', 'lwp']   ; OK
+; arr_pro = ['cmb', 'cph', 'ctt', 'cth', 'ctp', 'cod', 'ref', 'lwp']
+  arr_pro = ['cmb', 'cph', 'cod', 'ref', 'lwp']
+;  arr_pro = ['cmb2', 'cph2','ctp2', 'ctt2', 'cod2', 'ref2']
+  npro = n_elements(arr_pro)
+
+  day   = [13,17,18,22,03]
+  month = [ 6, 6, 6, 6, 7]
+
+  for d = 0, 4 do begin 
+;  for d = 0, 0 do begin 
+    for p = 0, npro-1 do begin 
+      product = arr_pro[p]
+      o->set_product, product
+      for h=0,23 do begin 
+;      for h=12,12 do begin 
+        for m=00,45,15 do begin 
+;        for m=00,00,15 do begin 
+          print, '=== TIME: ', h, m
+          o->set_date,2008,month[d],day[d],h,m
+            o->transformIt, /force
+        endfor
+      endfor
+    endfor
+  endfor 
+
+end

calipso_cth.pro

+pro calipso_cth, start, nr , overpass= overpass,noPlot = noPlot
+
+
+if n_params() ne 2 then nr = 0
+
+overpass =  n_elements(overpass) eq 1 ? overpass : 11318
+
+
+outFile = '/data4/cws2/results/atrain/calipso_cth/cth_calipso_' $
+          +string(overpass,format='(i5.5)')+'_'+string(nr,format='(i2.2)')
+file_mkdir,file_dirname(outFile)
+
+o = obj_new('cws2',overpass=overpass)
+
+o->atrain::set_overpass,overpass,/allow_nonexi
+grpList = ['CMS','EUM','OCA','MPF','FUB','DLR','MFR','UWM','UKM','GSF','LAR'] ;  groups
+ind = indgen(400)+start
+
+
+rfl = o->get_data(pro='TAB_1064',sens='CALIOP',/data)
+height = o->get_data(sen='CALIOP'     ,pro='HGT',/DATA)
+
+rfl = rfl[*,ind]
+height = height[*,ind]
+
+time = o->get_data(     pro='YDT',/DATA)
+
+time = time[ind]
+
+lon  = o->get_data( pro='lon',/DATA)
+lat  = o->get_data( pro='lat',/DATA)
+
+
+oMSG = obj_new('msg_data_cl')
+caldat,time[0]/3600./24.+julday(5,24,1968,0,0),month,day,year,hour,minute
+minute = 15*(minute/15)
+oMSG -> set_date,year,month,day,hour,minute
+rgb=oMSG->get_product('rgb')
+
+latMSG = oMSG->lat()
+
+
+ !p.multi=[0,4,3]
+;win,1,xs=1100,ys=800
+;win,2,xs=1100,ys=800
+if ~keyword_set(noPlot) then startps,_extra=_extra,outname=outFile,xs=4,ys=2.
+
+view2d,rgb ,/asp,/no_ax,space_idx = where(latMSG lt -100),space_col = !p.background,pos=[-0.1,0.2,0.3,0.8]
+
+
+
+pxl = geo_2_msg(lon,lat)
+;plots,pxl.(0),pxl.(1),color=fsc_color('pink')
+
+
+pxl = geo_2_msg(lon[ind],lat[ind])
+
+plots,pxl.(0),pxl.(1),color=fsc_color('red')
+
+!p.multi=[12,4,3]
+
+xyouts,0.0,0.95,'overpass '+string(overpass),/norm,charsize=0.6
+
+view2d,rgb[*,(min(pxl.(0)) - 200):(max(pxl.(0)) + 200),(min(pxl.(1)) - 200):(max(pxl.(1)) + 200) ],/asp,/no_ax,position =[0.45,0.1,0.99,0.9 ]
+
+plots,pxl.(0)-(min(pxl.(0)) - 200),pxl.(1)-(min(pxl.(1)) - 200),color=fsc_color('red')
+stop
+
+for u=0,n_elements(sfcbin)-1 do  radar[sfcbin[u]:*,u]= -999.
+
+
+ ; cth_cpr  = cth_cpr[ind]
+  ; radar  = radar [*,ind]
+    height = height[*,ind]
+    sfcbin = sfcbin[ind]
+    
+ h = REVERSE(FINDGEN(99)*0.25) + 0.5 ; this just for quicklook height bins starting 500 m above surfcae
+    l = LINDGEN(N_ELEMENTS(ind))
+   
+zind =(99-FINDGEN(99))  #  (FLTARR(N_ELEMENTS(sfcbin))+1)
+    zind =(FLTARR(99)+1) #  (sfcbin-2)  - 2 - zind
+    ;horizontal index
+    lind =(FLTARR(99)+1) # FINDGEN(N_ELEMENTS(sfcbin)) 
+
+    r    = TRANSPOSE(radar[zind,lind])
+   
+  
+   
+ 
+
+loadct,39
+for i= 0 , n_elements(grpList) -1 do begin
+  print,grpList[i]
+  data = o->_get_data(sensor='MSG',affi=grpList[i],product='CTP')
+  cmb  = o->_get_data(sensor='MSG',affi=grpList[i],product='CMB')
+  cmb = cmb[ind]
+  
+  if n_elements(data) le 1 then continue
+  
+  idx = where(data le 50,c)
+  if c gt 0 then data[idx] = !VALUES.F_NAN
+  
+  cth_affi = o->atrain::p2z(data[ind,*],ind)
+
+     
+  loadct,11
+  for p=2,60 do tvlct,fsc_color('gray',/triple),p
+    contour,r,l,h,/FILL,/color,NLEV=10,_EXTRA=_EXTRA,YRANGE=[0,20],$
+            YSTY=5,XStyle=5,XRANGE=[MIN(l),MAX(l)] $
+            ,title=group,charsize=1., pos =[0.1,0.68,0.9,0.95],/norm
+loadct,39
+
+;
+;    contour,r,l,h,/FILL,/color,NLEV=20,_EXTRA=_EXTRA,YRANGE=[0,20],$
+;            YSTY=5,XStyle=5,XRANGE=[MIN(l),MAX(l)] $
+;            ,title=grpList[i],charsize=1.
+    oplot,cth_affi/1000.,color=fsc_color('white'),psym=1,symsize=0.14
+	oplot,height[sfcbin,indgen(400)]/1000.
+   
+    axis,yax=0,yrange=[0,20]  ,yticks=yticks,YSTY=1,ytitle='km'
+    axis,yax=1,yticks=yticks,ytickform='(A1)',yminor=1 ; draw right y axis w/o labels.
+    axis,xax=0,xticks=1,xtickform='(A1)',xminor=1 ; draw lower x axis w/o labels.
+    axis,xax=1,xticks=1,xtickform='(A1)',xminor=1 ; draw lower x axis w/o labels.
+	
+	
+	cloud_idx = where(cmb gt 1.5,c,compl=compl)
+	if c gt 0 then plots,cloud_idx,make_array(c,value=-0.5),color=fsc_color('blue'),psym=3
+	c_compl = n_elements(compl)
+	if  c_compl gt 0 then plots,compl,make_array(c_compl,value=-0.8),color=fsc_color('gold'),psym=1,symsize=0.12
+	
+	;cth_direct = o->_get_data(sensor='MSG',affi=grpList[i],product='CTH')
+	
+	; wset,2
+; 	 !p.multi=[12-i,4,3]
+; 	plot,cth_affi/1000.,cth_cpr/1000.,psym=3,title=grpList[i],xrange=[0,12],yrange=[0,12]
+; 	oplot,[0,10],[0,10],col=fsc_color('red')
+	;if n_elements(cth_direct) le 1 then continue
+    ; oplot,cth_direct[ind]/1000.,color=fsc_color('gold'),psym=5,symsize=0.3
+	;  xyouts,80,13.5,'cth_direct',/data,charsize=14.
+
+endfor
+
+
+
+
+
+; contour,r,l,h,/FILL,/color,NLEV=20,_EXTRA=_EXTRA,YRANGE=[0,15],$
+;             YSTY=5,XStyle=5,XRANGE=[MIN(l),MAX(l)] $
+;             ,title='CPR',charsize=1.
+;     oplot,cth_cpr/1000.,color=fsc_color('red'),psym=1,symsize=0.3
+
+
+spawn,'convert -density 160 '+outfile $
+               +'.eps '+outfile+'.jpg'
+if ~keyword_set(noPlot)  then endps
+stop
+
+
+
+
+
+end

cloudsat_all_params.pro

+pro cloudsat_all_params, start , nr $
+    ,overpass=overpass $ 
+    ,cut = cut $
+    ,noPS = noPS $
+    , group = group $
+    , idx = idx $
+    , sensor = sensor 
+
+
+default,sensor,'MSG'
+
+if keyword_set(idx) then begin
+  start = idx[0]
+  nr_idx = max(idx)-min(idx)
+
+endif else nr_idx = 400
+
+if n_params() ne 2 then nr = 0
+
+;colList=['brown','orange','Royal Blue','Cyan','Green','Gold','Violet','Dark Green',$
+;    'Navy','Sandy Brown','Red','Dark Gray']
+@LIST_COLORS
+
+default,overpass,11317
+default,group,'UWM'
+default,start,400
+
+outFile = !PROJECTS.cws_path+'/results/atrain/cloudsat/cloudsat_all_' $
+        +group+'_' +string(overpass,format='(i5.5)')+'_'+string(nr,format='(i4.4)')
+file_mkdir,file_dirname(outFile)
+
+
+o = obj_new('cws2',overpass=overpass)
+o->atrain::set_overpass,overpass,/allow_nonexi
+;grpList = ['CMS','EUM','OCA','MPF','FUB','DLR','MFR','UWM','UKM','GSF','LAR'] ;  groups
+
+ind = indgen(nr_idx)+start
+
+rfl = o->get_data(pro='RFL',sens='CPR',/data)
+radar  = o->get_data(sen='CPR'     ,pro='RFL',/DATA)
+height = o->get_data(sen='CPR'     ,pro='HGT',/DATA)
+sfcbin = o->get_data(sen='CPR'     ,pro='SHB',/DATA)
+time = o->get_data(pro='YDT',/DATA)
+cth_cpr =  o->get_data(sen='CPR', prod='CTH',/data)
+
+time = time[ind]
+
+lon  = o->get_data( pro='lon',/DATA)
+lat  = o->get_data( pro='lat',/DATA)
+
+cod_modis = o->get_data(sen='MODIS', prod='COD',/data)
+cod_modis = cod_modis[ind]
+
+ref_modis = o->get_data(sen='MODIS', prod='CPS_REF',/data)
+ref_modis = ref_modis[ind]
+
+lwp_amsr = o->get_data(sensor='AMSR', prod = 'LWP',/data)
+lwp_amsr = lwp_amsr[ind]
+
+for u=0,n_elements(sfcbin)-1 do  radar[sfcbin[u]:*,u]= -999.
+
+
+cth_cpr  = cth_cpr[ind]
+radar  = radar [*,ind]
+height = height[*,ind]
+sfcbin = sfcbin[ind]
+    
+h = REVERSE(FINDGEN(99)*0.25) + 0.5 ; this just for quicklook height bins starting 500 m above surfcae
+l = LINDGEN(N_ELEMENTS(ind))
+   
+zind =(99-FINDGEN(99))  #  (FLTARR(N_ELEMENTS(sfcbin))+1)
+zind =(FLTARR(99)+1) #  (sfcbin-2)  - 2 - zind
+    ;horizontal index
+lind =(FLTARR(99)+1) # FINDGEN(N_ELEMENTS(sfcbin)) 
+
+r    = TRANSPOSE(radar[zind,lind])
+   
+  
+   
+ 
+prdList = ['ctp','cod','ref','ctt']
+loadct,39
+
+
+  data = o->_get_data(sensor=sensor,affi=group,product='CTP')
+  cmb  = o->_get_data(sensor=sensor,affi=group,product='CMB')
+  cmb = round(cmb[ind])
+  cph  = o->_get_data(sensor=sensor,affi=group,product='CPH')
+  cph = round(cph[ind])
+  cod = o->_get_data(sensor=sensor,affi=group,product='COD')
+  cod = cod[ind]
+  ref = o->_get_data(sensor=sensor,affi=group,product='REF')
+  ref = ref[ind]
+  ctt = o->_get_data(sensor=sensor,affi=group,product='CTT')
+  ctt = ctt[ind]
+  lwp = o->_get_data(sensor=sensor ,affi=group,product='LWP')
+  lwp = lwp[ind]
+  
+  
+  idx = where(cmb eq 1,cCmb)
+  if cCmb gt 0 then cph[idx] = 4  
+  
+  if n_elements(data) le 1 then begin
+    print,'no data'
+    return
+  endif
+  
+  idx = where(data le 50,c)
+  if c gt 0 then data[idx] = !VALUES.F_NAN
+  
+  cth_affi = o->atrain::p2z(data[ind,*],ind)
+  thisP = !p
+  !p.charsize=0.8
+  if ~keyword_set(noPS) then startps,outname=outfile,xs=2,ys=2
+    contour,r,l,h,/FILL,/color,NLEV=20,_EXTRA=_EXTRA,YRANGE=[0,20],$
+            YSTY=5,XStyle=5,XRANGE=[MIN(l),MAX(l)] $
+            ,title=group+(sensor eq 'MSG' ? '' : ' ('+sensor+')'),charsize=1. $
+            , pos =[0.1,0.68,0.9,0.95],/norm
+    oplot,cth_affi/1000.,color=fsc_color('white'),psym=1,symsize=0.44
+	  oplot,height[sfcbin,indgen(nr_idx)]/1000.
+   
+    axis,yax=0,yrange=[0,20]  ,yticks=yticks,YSTY=1,ytitle='km'
+    axis,yax=1,yticks=yticks,ytickform='(A1)',yminor=1 ; draw right y axis w/o labels.
+    axis,xax=0,xticks=1,xtickform='(A1)',xminor=1 ; draw lower x axis w/o labels.
+    axis,xax=1,xticks=1,xtickform='(A1)',xminor=1 ; draw lower x axis w/o labels.
+	
+	 
+	cloud_idx = where(cmb gt 1.5,c,compl=compl)
+	if c gt 0 then plots,cloud_idx,make_array(c,value=-0.5),color=fsc_color('blue'),psym=3
+	
+	no_cloud_idx = where(cmb eq 1,cNoCloud) 
+	if  cNoCloud gt 0 then plots,compl,make_array(cNoCloud,value=-0.8),color=fsc_color('gold'),psym=1,symsize=0.12
+	
+	phase_idx_wat = where(cph eq 1,cWat )
+	if cWat gt 1 then plots,phase_idx_wat,make_array(cWat,value=-1.3),color=fsc_color('green'),psym=3
+	phase_idx_ice = where(cph eq 2,cIce )
+	if  cIce gt 1 then plots,phase_idx_ice,make_array(cIce,value=-1.8),color=fsc_color('red'),psym=1,symsize=0.12
+	
+	phase_idx_mix = where(cph eq 3,cMix )
+  if  cMix gt 1 then plots,phase_idx_mix,make_array(cMix,value=-2.3),color=fsc_color('purple'),psym=1,symsize=0.12
+  
+  thisP = !P
+!p.charsize /= 2.4
+xyouts,nr_idx+10, -0., 'Cl',color=fsc_color('blue')
+xyouts,nr_idx+10,-0.8, 'NCl',color=fsc_color('gold')
+xyouts,nr_idx+10,-1.6, 'Wat',color=fsc_color('green')
+xyouts,nr_idx+10,-2.4, 'Ice',color=fsc_color('red')
+xyouts,nr_idx+10,-3.2, 'Mix',color=fsc_color('purple')
+!p = thisP
+  
+  
+  cod_wat = cod * (cph eq 1) - 999.*(cph ne 1)
+;  lwp = (0.833*0.65 * ref*(cod<50)) * (cph eq 1) + (((cod<50)^(1/0.84))/0.065 )
+; ; lwp = (0.833*0.65 * ref*cod)
+;  if size(lwp_amsr,/tname) eq 'STRING' then lwp_amsr=make_array(size(lwp_amsr,/dim),val=0)
+;  plot,lat[ind],(lwp_amsr>0)*1000. $
+;               ,pos=[0.1,0.12,0.9,0.28],/noerase $
+;              , yRange=[00,1500] ,symsize=0.3 $
+;              , yTitle = 'lwp [mm]',/xstyle,xtitle='lat',/nodata
+;    xyouts,0.5,0.3,'Liquid water path',charsize=0.8
+;  oplot,lat[ind],(lwp>0) $
+;          , color=fsc_color('blue')
+;       
+;  oplot,lat[ind],(lwp_amsr>0)*1000. $
+;          , color=fsc_color('green')        
+;         
+;  legend,[group, 'AMSR-E' ] $
+;             ,linestyle = [0,0] $
+;              , charsize = 0.6    $
+;              , box = 1          $
+;              , spacing = 0.3 $
+;              ,color=[fsc_color('blue')$
+;                 ,fsc_color('green')] 
+;  
+  
+  plot,lat[ind],cod,pos=[0.1,0.48,0.9,0.64],/noerase $
+              , yRange=[0,50],psym=1 $
+              , xTickformat='NOLABEL' $
+              , yTitle = 'cod[]' $
+              , symsize=0.3,/nodata,/xstyle
+              
+  oplot,lat[ind],cod,psym=1 $
+                 ,symsize=0.6,col=fsc_color('red')   
+  
+ oplot,lat[ind],cod_wat,psym=1 $
+                 ,symsize=0.6,col=fsc_color('green')  
+                 
+   oplot,lat[ind],cod_modis,psym=1 $
+                 ,symsize=0.4,col=fsc_color('blue')              
+                 
+   legend,['Water', ' Ice' $
+             ,' MODIS'] $
+             ,psym = [1,1,1],color=[fsc_color('green')$
+                 ,fsc_color('red'),0] $
+              , charsize = 0.6    $
+              , box = 0          $
+              , spacing = 0.3 ,/Fill
+
+  xyouts,0.3,0.97,group,align=0.5,charsize=3.$
+                ,col=fsc_color('white'),/norm
+
+;  plot,indgen(400)+start,ctt $
+;               ,pos=[0.1,0.1,0.9,0.36],/noerase $
+;              , yRange=[200,300],psym=1 ,symsize=0.3 $
+;              , yTitle = 'ctt [K]'
+ ;  axis,yax=1,yrange=[0,300], yTitle= 'lwp [gm2]'
+ ; oplot,indgen(400)+start,lwp/3.+200,color=fsc_color('green')
+  ref_wat = ref * (cph eq 1) - 999.*(cph ne 1)
+   plot,lat[ind],ref $
+               ,pos=[0.1,0.3,0.9,0.46],/noerase $
+              , yRange=[00,50],psym=1 ,symsize=0.3 $
+              , yTitle = 'ref['+textoidl('\mu')+'m]',/nodata $
+               , xTickformat='NOLABEL' ,/xstyle
+              
+    oplot,lat[ind],ref,psym=1 $
+                 ,symsize=0.6,col=fsc_color('red')
+   oplot,lat[ind],ref_wat,psym=1 $
+                 ,symsize=0.6,col=fsc_color('green')
+              
+   oplot,lat[ind],ref_modis,psym=1 $
+                 ,symsize=0.4,col=fsc_color('blue')   
+ 
+ legend,[' Water', ' Ice' $
+             ,' MODIS'] $
+             ,psym = [1,1,1],color=[fsc_color('green')$
+                 ,fsc_color('red'),0] $
+              , charsize = 0.6    $
+              , box = 0          $
+              , spacing = 0.3 
+ 
+ 
+ 
+                 
+if ~keyword_set(noPS)  then begin
+  spawn,'convert -density 160 '+outfile $
+               +'.eps '+outfile+'.jpg'
+  spawn, 'rm -f '+outfile+'.eps'
+  endps
+endif
+
+
+end
+
+
+
+

cmpr_modis.pro

+pro cmpr_modis,day=day,product = product
+
+default,product,'COD'
+default,day,'165'
+
+case product of
+  'COD': begin
+    modis_hdf_name = 'Cloud_Optical_Thickness'
+  
+  end
+  else:
+endcase
+
+
+
+
+oCWS = obj_new('cws_read')
+seviri_lon = oCWS->lon()
+seviri_lat = oCWS->lat()
+
+modis_path = '/data4/cws2/data/source/MODIS/NOM/'
+index_path = '/data4/cws2/data/indexfile_'
+
+for i=12,12 do begin
+  print,i
+  nominal_time = string(i,form='(i2.2)')+'00'
+  nominal_sec = 3600. * fix(strMid(nominal_time,0,2)) $
+      + 60. * fix(strMid(nominal_time,2,2))
+  sec_seviri   = (nominal_sec $
+        + indgen(3712)* (760./3712.)) $
+		   ## make_array(3712,/double,value=1D)
+		   
+  for j = 0 , 2 do begin
+    time_modis = nominal_time + ( j * 5)
+	modisFile = file_search(modis_path,'/MOD06_L2.A2008' $
+	  + string(day,format='(i3.3)') + $
+    '.'+ string(time_modis, format='(i4.4)')+'*.hdf',c=nModisfile)
+	
+	if nModisFile eq 0 then continue
+	
+	modisFile = modisFile[0]
+ 
+	
+	modis_ref = get_hdf4_sd(modisFile,'Cloud_Effective_Radius',$
+		scale_name = modis_scale_name, offset_name = modis_offset_name,$
+		missing_name = modis_missing_name)
+		
+	; Full resolution dimensions
+	
+	full_resolution_dim = (size(modis_ref))[1:2]
+    fu_x = full_resolution_dim[0] 
+    fu_y = full_resolution_dim[1] 
+	
+	
+	; modis time seconds of the day 
+    sec_modis_low = get_hdf4_sd(modisFile,'Scan_Start_Time') mod (3600D *24) 
+    sec_modis_vec = congrid(reform(sec_modis_low[4,*]),(size(modis_ref,/DIM))[1],/interp)
+    sec_modis = sec_modis_vec ## make_array((size(modis_ref,/DIM))[0],value=1)
+	
+	;geo-coordinates MODIS
+    lat_modis = get_hdf4_sd(modisFile,'Latitude',scale_name='scale_factor')
+    lon_modis = get_hdf4_sd(modisFile,'Longitude',scale_name='scale_factor')
+    lat_modis = congrid(lat_modis,fu_x,fu_y,/interp)
+    lon_modis = congrid(lon_modis,fu_x,fu_y,/interp)
+	
+	
+	;stupid, but it works
+	if max(lon_modis) lt -70 then continue
+	if min(lon_modis) gt 70 then continue
+	if max(lat_modis) lt -70 then continue
+	if min(lat_modis) gt 70 then continue
+	
+	idx= where(between(lat_sev,min(lat_modis),max(lat_modis)) $
+	           AND between(lon_sev,min(lon_modis),max(lon_modis))  $
+			   AND between(sec_seviri,min(sec_modis_vec),max(sec_modis_vec)) $
+	           AND between(sec_seviri,min(sec_modis_vec),max(sec_modis_vec)),cIdx)
+	
+	if cIdx eq 0 then begin
+	  print,'not in the spatial or time window'
+	  continue
+	endif
+	
+	lat_sev_scene = lat_sev[idx]
+    lon_sev_scene = lon_sev[idx]
+	sec_seviri_scene = sec_seviri[idx]
+	
+	index_file =index_path+'SEVIRI'+string(day,format='(i3.3)') $
+	   +string(i,format='(i2.2)')+'_MODIS'+file_basename(modisFile)+'.sav'
+	 if file_test(index_file) eq 0 then begin
+    ; find nearest neighbor and the distance to each seviri
+      idx_seviri = nearest_neighbor(lat_sev_scene,lon_sev_scene,lat_modis,lon_modis,distance=distance,/debug)
+	  save,idx_seviri,filename=index_file
+	endif else begin
+	  restore,index_file
+	 
+	endelse  
+	save,idx,idx_seviri,filename=index_file+'new.sav'
+	 
+	
+  endfor
+
+
+endfor
+
+
+
+end

cmpr_modis_2.pro

+pro cmpr_modis_2 $
+ , day = day $
+ , product = product $
+ , single = single $
+ , phase = phase $
+ , knmi = knmi
+
+
+added_name ='A' 
+if keyword_set(knmi) then added_name += '_KNMI'
+single =1
+cws2_startup
+
+default,product,'ref'
+default,day,'165'
+
+product = strLowCase(product)
+
+modis_path = !CWS2.data_path+'source/MODIS/NOM/20080613/'
+index_path = !CWS2.data_path+'index_files/'
+
+case product of
+  'cod': begin
+    modis_hdf_name = 'Cloud_Optical_Thickness'
+    binH = 0.25
+    minH=0
+    maxH=50
+    formH = '(f5.2)'
+    unit =''
+  end
+  
+  'ref' : begin
+    modis_hdf_name = 'Cloud_Effective_Radius'
+    binH = 0.25
+    minH=0
+    maxH=50
+    formH = '(f5.2)'
+    unit = '\mum'
+  end
+  
+  'ctp' : begin
+    modis_hdf_name = 'Cloud_Top_Pressure'
+    binH = 25
+    minH=0
+    maxH=1000
+  formH = '(f5.0)'
+    unit='hPa'
+  
+  end
+  
+  else:
+endcase
+
+modis_scale_name='scale_factor'
+modis_offset_name = 'add_offset'
+modis_missing_name='_FillValue'
+
+
+
+oCWS = obj_new('cws_read')
+
+seviri_lon = oCWS->lon()
+seviri_lat = oCWS->lat()
+oCWS ->set_date, 2008,6,13,12,0
+groups_data = oCWS->get_groups(product)
+
+histEle = 1+((maxH-minH)/binH)
+
+erg = replicate({d2Hist : dblarr(histEle,histEle)},n_elements(groups_data) )
+
+for i=12,16 do begin
+  for k=0,3 do begin
+    print,i,k
+  
+    nominal_time = string(i,form='(i2.2)')+string(k*15,format='(i2.2)')
+    nominal_sec = 3600. * fix(strMid(nominal_time,0,2)) $
+                  + 60. * fix(strMid(nominal_time,2,2)) 
+    sec_seviri = (nominal_sec $
+                  + indgen(3712)* (760./3712.)) $
+		              ## make_array(3712,/double,value=1D)
+		   
+    for j = 0 , 2 do begin
+  
+      time_modis = nominal_time + ( j * 5)
+	
+	    modisFile = file_search(modis_path,'/MYD06_L2.A2008' $
+	                + string(day,format='(i3.3)') + $
+                  '.'+ string(time_modis, format='(i4.4)')+'*.hdf',c=nModisfile)
+	
+	    print,time_modis,file_basename(modisFile[0]),file_test(modisFile[0])
+
+	    if nModisFile eq 0 then continue
+	
+	    modisFile = modisFile[0]
+	    
+	    
+	    if keyword_set(knmi) then begin
+	      knmi_file = file_search(!CWS2.data_path+'source/MODIS/KNMI_2/20080613/' $
+	                      ,'/CPP4MOD.AQUA.A2008' $
+                  + oCWS->get_date(/string,/day_of_year) + $
+                  '.'+ string(time_modis, format='(i4.4)')+'.h5',c=nKnmiFile) 
+                  
+	      if nKnmiFile eq 0 then continue
+	      
+	      oKNMI = h5_read_odata(knmi_file[0],'CPP_Results/'+modis_hdf_name)
+      
+        modis_data_str = oKnmi->get('all')
+      
+        knmi_data = modis_data_str.data
+	    
+	    endif else begin
+ 
+	    modis_data = get_hdf4_sd(modisFile, modis_hdf_name,$
+		               scale_name = modis_scale_name, offset_name = modis_offset_name,$
+		               missing_name = modis_missing_name)
+		 endelse
+		 
+		 modis_data = get_hdf4_sd(modisFile, modis_hdf_name,$
+                   scale_name = modis_scale_name, offset_name = modis_offset_name,$
+                   missing_name = modis_missing_name)
+;      win,4
+;      plot,knmi_data,modis_data,xRange=[0,60],yRange=[0,60],psym=3
+;      win,5
+;      view2d,knmi_data,/cool,/colo,min=0,no_data_idx=where(knmi_data lt 0)
+      win,6
+      view2d,modis_data,/cool,/colo,min=0,no_data_idx=where(modis_data lt 0)
+      win,3
+		 
+	    modis_full = get_hdf4_sd(modisFile, 'Cloud_Optical_Thickness' ,$
+                   scale_name = modis_scale_name, offset_name = modis_offset_name,$
+                   missing_name = modis_missing_name)
+		
+	; Full resolution dimensions
+	
+	    full_resolution_dim = size(modis_full,/dim)
+      fu_x = full_resolution_dim[0] 
+      fu_y = full_resolution_dim[1] 
+    
+      if fu_x ne (size(modis_data,/dim))[0] then begin
+        modis_data = congrid(modis_data,fu_x,fu_y)
+      endif
+    
+	
+	
+	; modis time seconds of the day 
+      sec_modis_low = get_hdf4_sd(modisFile,'Scan_Start_Time') mod (3600D *24) 
+      sec_modis_vec = congrid(reform(sec_modis_low[4,*]),(size(modis_data,/DIM))[1],/interp)
+      modisSeconds = sec_modis_vec ## make_array((size(modis_data,/DIM))[0],value=1)
+	
+	;geo-coordinates MODIS
+      modis_lat = get_hdf4_sd(modisFile,'Latitude',scale_name='scale_factor')
+      modis_lon = get_hdf4_sd(modisFile,'Longitude',scale_name='scale_factor')
+      modis_lat = congrid(modis_lat,fu_x,fu_y,/interp)
+      modis_lon = congrid(modis_lon,fu_x,fu_y,/interp)
+	
+	
+	
+	    if max(modis_lon) lt -70 then continue
+	    if min(modis_lon) gt 70 then continue
+	    if max(modis_lat) lt -70 then continue
+	    if min(modis_lat) gt 70 then continue
+	
+	    idx= where(between(seviri_lat,min(modis_lat),max(modis_lat)) $
+	           AND between(seviri_lon,min(modis_lon),max(modis_lon))  $
+			       AND between(sec_seviri,min(sec_modis_vec),max(sec_modis_vec)) $
+	           AND between(sec_seviri,min(sec_modis_vec),max(sec_modis_vec)),cIdx)
+	
+	    if cIdx eq 0 then begin
+	      print,'not in the spatial or time window'
+	      continue
+	    endif
+	
+	    seviri_lat_scene = seviri_lat[idx]
+      seviri_lon_scene = seviri_lon[idx]
+	    seviriSecond_scene = sec_seviri[idx]
+	
+	    index_file = index_path+'SEVIRI_'+string(day,format='(i3.3)') $
+	                 +string(i,format='(i2.2)')+string(k*15,format='(i2.2)')+'__MODIS_' $
+	                 +file_basename(modisFile)+'.sav'
+	
+	     if file_test(index_file) eq 0 then begin
+       ; find nearest neighbor and the distance to each seviri
+         idx_seviri = nearest_neighbor(seviri_lat_scene,seviri_lon_scene $
+	                     ,modis_lat,modis_lon,distance=distance,/debug)
+	       save,idx,idx_seviri,filename=index_file
+	     endif else begin
+	       restore,index_file	 
+	     endelse  
+	
+	     dist_meter = great_circle_dist(seviri_lon_scene , seviri_lat_scene $
+	                  ,modis_lon[idx_seviri] ,modis_lat[idx_seviri] )
+	
+	     time_shift = abs(( modisSeconds[idx_seviri])-seviriSecond_scene)
+	
+	     idx_match = where( time_shift lt 600 and dist_meter lt 2000,dum ) 
+	     if dum eq 0 then continue
+	     seviri_lat_match = seviri_lat_scene[idx_match]
+	     seviri_lon_match = seviri_lon_scene[idx_match]
+	 
+	     modis_lon_match = modis_lon[idx_seviri[idx_match]]
+	     modis_lat_match = modis_lat[idx_seviri[idx_match]]
+	 
+	     modis_data_match = modis_data[idx_seviri[idx_match]]
+	 
+	     if n_elements(phase) eq 1 then begin
+	     
+	       modis_phase = get_hdf4_sd(modisFile,'Cloud_Phase_Infrared',$
+             scale_name = modis_scale_name, offset_name = modis_offset_name,$
+               missing_name = modis_missing_name)
+               
+         if fu_x ne (size(modis_phase,/dim))[0] then begin
+           modis_phase = congrid(modis_phase,fu_x,fu_y)
+         endif
+	       modis_phase_match_this = modis_phase[idx_seviri[idx_match]]
+	       modis_data_match_this = modis_data_match
+	     endif
+	 
+	; window,2,xsize=600,ySize=900
+	!p.multi = [0,2,1+n_elements(groups_data)/2]
+	oCWS ->set_date, 2008,6,13,i,k*15
+	 for gg = 0 , n_elements(groups_data) -1 do begin
+	   
+	   print,gg,groups_data[gg]
+	   seviri_data = oCWS -> get_data(gr=groups_data[gg],product=product)
+	   seviri_data_match = seviri_data[idx[idx_match]]
+	   
+	   if n_elements(phase) eq 1 then begin
+	     seviri_cph =  oCWS -> get_data(gr=groups_data[gg],product='cph')
+	     seviri_cph_match= seviri_cph[idx[idx_match]]
+	     added_name=phase
+	    
+	    case phase of
+         'W':begin
+            iWat = where((round(seviri_cph_match) eq 1) and (round(modis_phase_match_this) eq 1 ),cWat)  
+            if cWat ne 0 then begin
+              seviri_data_match = seviri_data_match[iWat]
+              modis_data_match = modis_data_match_this[iWat]
+            endif
+         end
+         'I': begin
+            iIce = where((round(seviri_cph_match) eq 2) and (round(modis_phase_match_this) eq 2 ),cIce) 
+           if cIce ne 0 then begin
+              seviri_data_match = seviri_data_match[iIce]
+              modis_data_match = modis_data_match_this[iIce]
+            endif
+          end
+         else:stop
+         endcase
+	   endif
+	   
+	   d2Hist = hist_2d(seviri_data_match,modis_data_match $
+	            ,min1=minH,min2=minH,max1=maxH,max2=maxH,bin1=binH,bin2=binH)
+	   erg[gg].(0)  += d2Hist
+	  
+	  statsH = get_stats_from2d(erg[gg].(0),max=maxH)
+	  print,minmax(erg[gg].(0))
+	  
+	   view2d,alog10(erg[gg].(0)),/cool $
+	                   ,no_data_idx=where(erg[gg].(0) eq 0) $
+	                   ,xrange=[minH,maxH]$
+	                   ,yrange=[minH,maxH] $
+	                   ,/asp $
+	                   ,title=groups_data[gg]$
+	                     ,charsize=1.6
+	   
+	   oplot,[0,maxH],[0,maxH],color=fsc_color('white')
+	   xyouts,0,maxH+1,string(statsH.bias,format=formH),/data	   
+	 endfor
+	 
+	
+  endfor
+
+endfor
+endfor
+
+
+
+outname = !CWS2.data_path+'/../results/modis/modis_'+product+added_name+'_scatterplot_2008'+day
+startps, out=outname, xs=2.,ys=2.
+
+thisP = !p 
+thisX = !x
+thisY = !y
+!p.multi = [0,3,1+n_elements(groups_data)/3]
+!y.omargin=[2,6]
+
+if product eq 'ctp' then begin
+  
+  xrng = [maxH,minH]
+  yrng = [maxH,minH]
+endif else begin
+  xrng = [minH,maxH]
+  yrng = [minH,maxH]
+endelse
+
+for gg=0,n_elements(groups_data) -1 do begin
+  if product eq 'ctp' then erg[gg].(0) = rotate(erg[gg].(0),2)
+  statsH = get_stats_from2d(erg[gg].(0),max=maxH)
+  
+                       
+                        view2d,erg[gg].(0),/cool $
+                 ,no_data_idx=where(erg[gg].(0) eq 0) $
+                     ,xrange= xrng $
+                     ,yrange= yrng $
+                     ,/asp $
+                     ,title=groups_data[gg]$
+                       ,charsize=0.8 $
+                       ,ytitle=textoidl(unit)+' [MODIS]' $
+                       , xTitle = textoidl(unit)+ ' [SEVIRI]' $
+                       ,/logarith,/colo $
+                       ,barformat = '(i5)' $
+                       ,bardist=0.2
+                       
+  oplot,[0,maxH],[0,maxH],color=fsc_color('white')
+  xyouts,xrng[0],yrng[1] + (product eq 'ctp' ? (-10) : 1) $
+        ,'bias: '+string(statsH.bias,format=formH),/data ,charsize= 0.48   
+    xyouts,xrng[1],yrng[1] + (product eq 'ctp' ? (-10) : 1),'cr: '+string(statsH.corr,format='(f4.2)'),/data $
+    ,charsize= 0.48 ,align=1
+endfor
+xyOuts,0.5,0.97,'SEVIRI / MODIS  '$
+    +strUpcase(product eq 'ctp' ? 'Cloud Top Pressure ': product) $
+    +' Phase '+added_name+' !C Day 20080613',align=0.5,/normal
+endps
+spawn,'convert -density 200 '+outname+'.eps '+outname+'.jpg'
+spawn,'rm -f '+outname+'.eps'
+!p=thisP
+!x = thisX
+!y = thisY
+; single
+
+if keyword_set(single) then begin
+
+  
+
+  thisP = !p 
+  !p.multi = 0
+  
+
+  if product eq 'ctp' then begin
+  
+    xrng = [maxH,minH]
+    yrng = [maxH,minH]
+  endif else begin
+    xrng = [minH,maxH]
+    yrng = [minH,maxH]
+  endelse
+
+  for gg=0,n_elements(groups_data) -1 do begin
+    outname = !CWS2.data_path+'/../results/modis/single/modis_'+product+added_name+'_'+groups_data[gg]+'_scatterplot_2008'+day
+    startps, out=outname, xs=1.,ys=1.
+    if product eq 'ctp' then erg[gg].(0) = rotate(erg[gg].(0),2)
+    statsH = get_stats_from2d(erg[gg].(0),max=maxH)
+    view2d,erg[gg].(0),/cool $
+                 ,no_data_idx=where(erg[gg].(0) eq 0) $
+                     ,xrange= xrng $
+                     ,yrange= yrng $
+                     ,/asp $
+                     ,title=product + ' '+added_name $
+                       ,charsize=0.6 $
+                       ,ytitle=textoidl(unit)+' [MODIS]' $
+                       , xTitle = textoidl(unit)+ ' [SEVIRI]' $
+                       ,/logarith,/colo $
+                       ,barformat = '(i5)'
+     
+    oplot,[0,maxH],[0,maxH],color=fsc_color('white')
+    xyouts,xrng[0],yrng[1] + (product eq 'ctp' ? (-10) : 1) $
+        ,'bias: '+string(statsH.bias,format=formH),/data ,charsize= 0.4   
+    xyouts,xrng[1],yrng[1] + (product eq 'ctp' ? (-10) : 1),'correlation: '+string(statsH.corr,format=formH),/data $
+    ,charsize= 0.4 ,align=1
+    endps
+    spawn,'convert -density 200 '+outname+'.eps '+outname+'.jpg'
+    spawn,'rm -f '+outname+'.eps'
+  endfor
+endif
+
+end
+
+
+pro make_everything
+cmpr_modis_2,product='cod'
+cmpr_modis_2,phase='I',product='cod'
+;
+cmpr_modis_2,phase='W',product='cod'
+;
+cmpr_modis_2,phase='I',product='ref'
+cmpr_modis_2,product='ref'
+cmpr_modis_2,phase='W',product='ref'
+
+cmpr_modis_2,phase='I',product='ctp'
+cmpr_modis_2,product='ctp'
+cmpr_modis_2,phase='W',product='ctp'
+
+
+
+end

cmpr_modis_angular.pro

+pro cmpr_modis_angular $
+ , day = day $
+ , product = product
+
+cws2_startup
+
+default,product,'cod'
+default,day,'165'
+
+product = strlowcase(product)
+
+modis_path = !CWS2.data_path+'source/MODIS/NOM/20080613/'
+index_path = !CWS2.data_path+'index_files/'
+
+case product of
+  'cod': begin
+    modis_hdf_name = 'Cloud_Optical_Thickness'
+    binH = 0.25
+    minH=0
+    maxH=50
+    formH = '(f5.2)'
+  end
+  
+  'ref' : begin
+    modis_hdf_name = 'Cloud_Effective_Radius'
+    binH = 0.25
+    minH=0
+    maxH=50
+  formH = '(f5.2)'
+  
+  end
+  
+  'ctp' : begin
+    modis_hdf_name = 'Cloud_Top_Pressure'
+    binH = 25
+    minH=0
+    maxH=1000
+  formH = '(f5.0)'
+  
+  end
+  
+  else:
+endcase
+
+modis_scale_name='scale_factor'
+modis_offset_name = 'add_offset'
+modis_missing_name='_FillValue'
+
+
+
+oCWS = obj_new('cws_read')
+
+seviri_lon = oCWS->lon()
+seviri_lat = oCWS->lat()
+oCWS ->set_date, 2008,6,13,12,0
+groups_data = oCWS->get_groups(product)
+
+result_struct = { $
+           dat_obs : make_array(71,201,value=0.) $
+          ,dat_sun : make_array(71,201,value=0.) $
+          ,dat_azi : make_array(91,201,value=0.) $
+      ,info:ptr_new() }
+
+
+erg = replicate(result_struct,n_elements(groups_data) )
+
+for i=6,17 do begin
+  for k=0,3 do begin
+    print,i,k
+  
+    nominal_time = string(i,form='(i2.2)')+string(k*15,format='(i2.2)')
+    nominal_sec = 3600. * fix(strMid(nominal_time,0,2)) $
+                  + 60. * fix(strMid(nominal_time,2,2)) 
+    sec_seviri = (nominal_sec $
+                  + indgen(3712)* (760./3712.)) $
+		              ## make_array(3712,/double,value=1D)
+		   
+    for j = 0 , 2 do begin
+  
+      time_modis = nominal_time + ( j * 5)
+	
+	    modisFile = file_search(modis_path,'/MOD06_L2.A2008' $
+	                + string(day,format='(i3.3)') + $
+                  '.'+ string(time_modis, format='(i4.4)')+'*.hdf',c=nModisfile)
+	
+	    print,time_modis,file_basename(modisFile[0]),file_test(modisFile[0])
+
+	    if nModisFile eq 0 then continue
+	
+	    modisFile = modisFile[0]
+ 
+	    modis_data = get_hdf4_sd(modisFile, modis_hdf_name,$
+		               scale_name = modis_scale_name, offset_name = modis_offset_name,$
+		               missing_name = modis_missing_name)
+		
+	    modis_full = get_hdf4_sd(modisFile, 'Cloud_Optical_Thickness' ,$
+                   scale_name = modis_scale_name, offset_name = modis_offset_name,$
+                   missing_name = modis_missing_name)
+		
+	; Full resolution dimensions
+	
+	    full_resolution_dim = size(modis_full,/dim)
+      fu_x = full_resolution_dim[0] 
+      fu_y = full_resolution_dim[1] 
+    
+      if fu_x ne (size(modis_data,/dim))[0] then begin
+        modis_data = congrid(modis_data,fu_x,fu_y)
+      endif
+    
+	
+	
+	; modis time seconds of the day 
+      sec_modis_low = get_hdf4_sd(modisFile,'Scan_Start_Time') mod (3600D *24) 
+      sec_modis_vec = congrid(reform(sec_modis_low[4,*]),(size(modis_data,/DIM))[1],/interp)
+      modisSeconds = sec_modis_vec ## make_array((size(modis_data,/DIM))[0],value=1)
+	
+	;geo-coordinates MODIS
+      modis_lat = get_hdf4_sd(modisFile,'Latitude',scale_name='scale_factor')
+      modis_lon = get_hdf4_sd(modisFile,'Longitude',scale_name='scale_factor')
+      modis_lat = congrid(modis_lat,fu_x,fu_y,/interp)
+      modis_lon = congrid(modis_lon,fu_x,fu_y,/interp)
+	
+	
+	
+	    if max(modis_lon) lt -70 then continue
+	    if min(modis_lon) gt 70 then continue
+	    if max(modis_lat) lt -70 then continue
+	    if min(modis_lat) gt 70 then continue
+	
+	    idx= where(between(seviri_lat,min(modis_lat),max(modis_lat)) $
+	           AND between(seviri_lon,min(modis_lon),max(modis_lon))  $
+			       AND between(sec_seviri,min(sec_modis_vec),max(sec_modis_vec)) $
+	           AND between(sec_seviri,min(sec_modis_vec),max(sec_modis_vec)),cIdx)
+	
+	    if cIdx eq 0 then begin
+	      print,'not in the spatial or time window'
+	      continue
+	    endif
+	    
+	    
+	
+	    seviri_lat_scene = seviri_lat[idx]
+      seviri_lon_scene = seviri_lon[idx]
+	    seviriSecond_scene = sec_seviri[idx]
+	
+	    index_file = index_path+'SEVIRI_'+string(day,format='(i3.3)') $
+	                 +string(i,format='(i2.2)')+string(k*15,format='(i2.2)')+'__MODIS_' $
+	                 +file_basename(modisFile)+'.sav'
+	
+	     if file_test(index_file) eq 0 then begin
+       ; find nearest neighbor and the distance to each seviri
+         idx_seviri = nearest_neighbor(seviri_lat_scene,seviri_lon_scene $
+	                     ,modis_lat,modis_lon,distance=distance,/debug)
+	       save,idx,idx_seviri,filename=index_file
+	     endif else begin
+	       restore,index_file	 
+	     endelse  
+	
+	     dist_meter = great_circle_dist(seviri_lon_scene , seviri_lat_scene $
+	                  ,modis_lon[idx_seviri] ,modis_lat[idx_seviri] )
+	
+	     time_shift = abs(( modisSeconds[idx_seviri])-seviriSecond_scene)
+	
+	     idx_match = where( time_shift lt 600 and dist_meter lt 2000,dum ) 
+	     if dum eq 0 then continue
+	     
+	     seviri_path = !cws2.data_path+'/source/SEVIRI/UWM/v3/'
+	     
+	     seviriFile = file_search(seviri_path,'*.2008' $
+                     +(oCWS->get_date(/day_of_year,/string)) $
+                     +'.'+string(i,form='(i2.2)') + string(k*15,format='(i2.2)')+'00.hdf',c=nSeviriFile)
+	     
+	     seviri_sun = read_geocat(seviriFile[0],'SUN')
+       seviri_obs = read_geocat(seviriFile[0],'OBS')
+       seviri_azi = read_geocat(seviriFile[0],'AZI')
+	     
+	     seviri_sun_match  = seviri_sun[idx[idx_match]]
+	     seviri_obs_match  = seviri_obs[idx[idx_match]]
+	     seviri_azi_match  = seviri_azi[idx[idx_match]]
+	     
+	     undefine,seviri_sun
+	     undefine,seviri_obs
+	     undefine,seviri_azi
+	     
+	     seviri_lat_match = seviri_lat_scene[idx_match]
+	     seviri_lon_match = seviri_lon_scene[idx_match]
+	 
+	     modis_lon_match = modis_lon[idx_seviri[idx_match]]
+	     modis_lat_match = modis_lat[idx_seviri[idx_match]]
+	 
+	     modis_data_match = modis_data[idx_seviri[idx_match]]
+	 
+	  
+	 
+	; window,2,xsize=600,ySize=900
+	!p.multi = [0,2,1+n_elements(groups_data)/2]
+	!p.multi=0
+	!p.charsize=1.6
+	oCWS ->set_date, 2008,6,13,i,k*15
+	; for gg = 4 , n_elements(groups_data) -1 do begin
+	   for gg = 4 , 4 do begin
+	   print,gg,groups_data[gg]
+	   seviri_data = oCWS -> get_data(gr=groups_data[gg],product=product)
+	   seviri_data_match = seviri_data[idx[idx_match]]
+	   
+	   idxM = where(between(modis_data_match,0,30) AND between(seviri_data_match,0,30),cidxM)
+	   if cidxM eq 0 then continue
+	   modis_data_match =  modis_data_match[idxM]
+	   seviri_data_match =  seviri_data_match[idxM]
+	   
+	   error_data = 100.*(modis_data_match - seviri_data_match)/modis_data_match
+  
+    data_sun = hist_2d(seviri_sun_match,error_data,min1=0,max1=70,min2=-100,max2=100,bin1=1,bin2=1)
+    data_obs = hist_2d(seviri_obs_match,error_data,min1=0,max1=70,min2=-100,max2=100,bin1=1,bin2=1)
+    data_azi = hist_2d(seviri_azi_match,error_data,min1=0,max1=180,min2=-100,max2=100,bin1=2,bin2=1.)
+  
+    erg[gg].dat_sun += data_sun
+    erg[gg].dat_obs += data_obs
+    erg[gg].dat_azi += data_azi
+	   
+	 view2d,(dum=(erg[gg].dat_obs)[*,2:198]),/cool,/colo $
+	      ,yRange=[-98,98],xRange=[0,70] $
+	      , title=groups_data[gg],/logarith $
+                       ,barformat = '(i5)' $
+                       ,no_data_idx=where(dum eq 0)
+	  
+	 endfor
+	 
+	
+  endfor
+
+endfor
+endfor
+stop
+
+
+end

cmpr_modis_knmi.pro

+pro cmpr_modis_knmi $
+ , day = day $
+ , product = product $
+ , single = single
+
+
+
+if keyword_set(knmi) then added_name += '_KNMI'
+single =1
+cws2_startup
+
+default,product,'ref'
+default,day,'165'
+
+product = strLowCase(product)
+
+modis_path = !CWS2.data_path+'source/MODIS/NOM/20080613/'
+index_path = !CWS2.data_path+'index_files/'
+
+case product of
+  'cod': begin
+    modis_hdf_name = 'Cloud_Optical_Thickness'
+    binH = 0.25
+    minH=0
+    maxH=100
+    formH = '(f5.2)'
+    unit =''
+  end
+  
+  'ref' : begin
+    modis_hdf_name = 'Cloud_Effective_Radius'
+    binH = 0.25
+    minH=0
+    maxH=60
+    formH = '(f5.2)'
+    unit = '\mum'
+  end
+  
+  'lwp' : begin
+    modis_hdf_name = 'Cloud_Water_Path'
+    binH = 1.
+    minH=0
+    maxH=800
+    formH = '(f5.2)'
+    unit = 'g/km2'
+  
+  
+  
+   end 
+  
+  
+  else:
+endcase
+
+modis_scale_name='scale_factor'
+modis_offset_name = 'add_offset'
+modis_missing_name='_FillValue'
+
+
+
+
+
+histEle = 1+((maxH-minH)/binH)
+
+erg = {d2Hist : dblarr(histEle,histEle)} 
+
+for i=12,16 do begin
+  for k=0,3 do begin
+    print,i,k
+  
+    nominal_time = string(i,form='(i2.2)')+string(k*15,format='(i2.2)')
+    nominal_sec = 3600. * fix(strMid(nominal_time,0,2)) $
+                  + 60. * fix(strMid(nominal_time,2,2)) 
+    sec_seviri = (nominal_sec $
+                  + indgen(3712)* (760./3712.)) $
+		              ## make_array(3712,/double,value=1D)
+		   
+    for j = 0 , 2 do begin
+  
+      time_modis = nominal_time + ( j * 5)
+	
+	    modisFile = file_search(modis_path,'/MYD06_L2.A2008' $
+	                + string(day,format='(i3.3)') + $
+                  '.'+ string(time_modis, format='(i4.4)')+'*.hdf',c=nModisfile)
+	
+	    print,time_modis,file_basename(modisFile[0]),file_test(modisFile[0])
+
+	    if nModisFile eq 0 then continue
+	
+	    modisFile = modisFile[0]
+	    
+	    
+	   
+	    knmi_file = file_search(!CWS2.data_path+'source/MODIS/KNMI_2/20080613/' $
+	                      ,'/CPP4MOD.AQUA.A2008' $
+                  + day + $
+                  '.'+ string(time_modis, format='(i4.4)')+'.h5',c=nKnmiFile) 
+                  
+	    if nKnmiFile eq 0 then continue 
+	      
+	    oKNMI = h5_read_odata(knmi_file[0],'CPP_Results/'+modis_hdf_name)
+	    if ~obj_valid(oKNMI) then continue
+      modis_data_str = oKnmi->get('all')
+      knmi_data = modis_data_str.data
+	    
+	    modis_data = get_hdf4_sd(modisFile, modis_hdf_name,$
+                   scale_name = modis_scale_name, offset_name = modis_offset_name,$
+                   missing_name = modis_missing_name)
+                   
+                   
+      modis_full = get_hdf4_sd(modisFile, 'Cloud_Optical_Thickness' ,$
+                   scale_name = modis_scale_name, offset_name = modis_offset_name,$
+                   missing_name = modis_missing_name)             
+       full_resolution_dim = size(modis_full,/dim)
+      fu_x = full_resolution_dim[0] 
+      fu_y = full_resolution_dim[1] 
+    
+      if fu_x ne (size(modis_data,/dim))[0] then begin
+        modis_data = congrid(modis_data,fu_x,fu_y)
+      endif
+    
+;
+      win,6
+      view2d,modis_data,/cool,/colo,min=0,no_data_idx=where(modis_data lt 0)
+      win,5
+		  view2d,knmi_data,/cool,/colo,min=0,no_data_idx=where(knmi_data lt 0)
+		   
+		      
+	    if n_elements(phase) eq 1 then begin
+	     
+	       modis_phase = get_hdf4_sd(modisFile,'Cloud_Phase_Infrared',$
+         scale_name = modis_scale_name, offset_name = modis_offset_name,$
+         missing_name = modis_missing_name)
+               
+         if fu_x ne (size(modis_phase,/dim))[0] then begin
+           modis_phase = congrid(modis_phase,fu_x,fu_y)
+         endif
+	      
+	    endif
+	 
+
+	 
+	   
+	   
+	   
+	   d2Hist = hist_2d(knmi_data,modis_data $
+	            ,min1=minH,min2=minH,max1=maxH,max2=maxH,bin1=binH,bin2=binH)
+	   erg.(0)  += d2Hist
+	  
+	  statsH = get_stats_from2d(erg.(0),max=maxH)
+	  print,minmax(erg.(0))
+	  
+	   view2d,erg.(0),/cool $
+	                   ,no_data_idx=where(erg.(0) eq 0) $
+	                   ,xrange=[minH,maxH]$
+	                   ,yrange=[minH,maxH] $
+	                   ,/asp $
+	                   ,title=product $
+	                     ,charsize=1.6 $
+	                     ,ytitle=textoidl(unit)+' [MYD06]' $
+                       , xTitle = textoidl(unit)+ ' [KNMI]' $
+	                     ,/colo $
+	                     ,/logar
+	   
+	   oplot,[0,maxH],[0,maxH],color=fsc_color('white')
+	   xyouts,0,maxH+1,string(statsH.bias,format=formH),/data	   
+	 endfor
+	 
+	
+  endfor
+
+endfor
+
+
+
+
+outname = !CWS2.data_path+'/../results/modis/KNMI_'+product+'_scatterplot_2008'+day
+startps, out=outname, xs=1.6,ys=1.2
+
+thisP = !p 
+thisX = !x
+thisY = !y
+!p.multi = 0
+!y.omargin=[2,6]
+
+
+  xrng = [minH,maxH]
+  yrng = [minH,maxH]
+
+
+
+  statsH = get_stats_from2d(erg.(0),max=maxH)
+  
+                       
+                        view2d,erg.(0),/cool $
+                 ,no_data_idx=where(erg.(0) eq 0) $
+                     ,xrange= xrng $
+                     ,yrange= yrng $
+                     ,/asp $
+                     ,title=product $
+                       ,charsize=0.8 $
+                       ,ytitle=textoidl(unit)+' [MYD06]' $
+                       , xTitle = textoidl(unit)+ ' [KNMI]' $
+                       ,/logarith,/colo $
+                       ,barformat = '(i5)' $
+                       ,bardist=0.2
+                       
+  oplot,[0,maxH],[0,maxH],color=fsc_color('white')
+  xyouts,xrng[0],yrng[1] + (product eq 'lwp' ? (10) : 1) $
+        ,'bias: '+string(statsH.bias,format=formH),/data ,charsize= 0.48   
+    xyouts,xrng[1],yrng[1] + (product eq 'lwp' ? (10) : 1),'cr: '+string(statsH.corr,format='(f4.2)'),/data $
+    ,charsize= 0.48 ,align=1
+
+
+endps
+spawn,'convert -density 200 '+outname+'.eps '+outname+'.jpg'
+spawn,'rm -f '+outname+'.eps'
+!p=thisP
+!x = thisX
+!y = thisY
+; single
+
+
+
+end
+
+
+pro make_everything
+;cmpr_modis_2,product='cod'
+;cmpr_modis_2,phase='I',product='cod'
+;
+;cmpr_modis_2,phase='W',product='cod'
+;
+;cmpr_modis_2,phase='I',product='ref'
+;cmpr_modis_2,product='ref'
+;cmpr_modis_2,phase='W',product='ref'
+
+cmpr_modis_knmi,phase='I',product='ctp'
+cmpr_modis_knmi,product='ctp'
+cmpr_modis_knmi,phase='W',product='ctp'
+
+
+
+end