pro flag_fraction_world, datelist, output_path


common cbtz_dz_variables, lat_vector, lon_vector,                 $
          cbtz_bin_mid, cbtz_dbin, dz_bin_mid, dz_dbin,           $
          cld_column_count, total_column_count,                   $
          cbz_dz1_grid, ctz_dz1_grid, cbz_dz2_grid, ctz_dz2_grid

; ------------------------------------------------------------------ ;
output_path_base = output_path+'World'+'_'+datelist[0] + '-' +  datelist[n_elements(datelist)-1]
;output_path_base = "/uufs/chpc.utah.edu/common/home/mace_grp/web/html/qzhang/ppp/S_World_200607_200902"
; ------------------------------------------------------------------ ;

max_lat=+70. & min_lat=-70. & comin_lon=0. & comax_lon=359. & dlat=3.0 & dlon=4.

; set up the first map
map_area='World'

;
max_lat=max_lat-dlat

colon_vector=lon_vector
result=where(colon_vector lt 0,count)
if count gt 0 then colon_vector[result]=360.0+colon_vector[result]

lat_vector2=min_lat
while max(lat_vector2) le max_lat do lat_vector2=[lat_vector2,(max(lat_vector2))+(2.*dlat)]

colon_vector2=comin_lon
while max(colon_vector2) le comax_lon do colon_vector2=[colon_vector2,(max(colon_vector2))+(2.*dlon)]

lon_vector2=colon_vector2
result=where(lon_vector2 gt 180., count)

if count gt 0. then lon_vector2[result]=colon_vector2[result]-360.

max_lon=comax_lon
if comax_lon gt 180. then max_lon=comax_lon-360.

min_lon=comin_lon
if comin_lon gt 180. then min_lon=comin_lon-360.

if max_lon lt min_lon then begin
    t=max_lon & max_lon=min_lon & min_lon=t
endif

layer_base_map=fltarr(n_elements(colon_vector2),n_elements(lat_vector2))

;for base_ind1 = 0, 0 do begin
;for base_ind2 = n_elements(cbtz_bin_mid)-2, n_elements(cbtz_bin_mid)-2 do begin
for base_ind1 = 0, n_elements(cbtz_bin_mid)-1 do begin
for base_ind2 = base_ind1, n_elements(cbtz_bin_mid)-1 do begin

;for dz_ind1 = 0, 0 do begin
;for dz_ind2 = n_elements(dz_bin_mid)-1, n_elements(dz_bin_mid)-1 do begin
for dz_ind1 = 0, n_elements(dz_bin_mid)-1 do begin
for dz_ind2 = n_elements(dz_bin_mid)-1, n_elements(dz_bin_mid)-1 do begin	
;for dz_ind2=dz_ind1,n_elements(dz_bin_mid)-1 do begin

   for j=0,n_elements(lat_vector2)-1 do begin

   i=0
   while lat_vector[i] le lat_vector2[j]-dlat and i lt n_elements(lat_vector)-2 do i=i+1 & lat_start_index=i
   while lat_vector[i] lt lat_vector2[j]+dlat and i lt n_elements(lat_vector)-2 do i=i+1 & lat_end_index=i

   for k=0,n_elements(colon_vector2)-1 do begin

   i=0
   if colon_vector2[k] lt 180 then begin
       while colon_vector[i] ne 0.0 do i=i+1
   endif

   while colon_vector[i] le colon_vector2[k]-dlon and i lt n_elements(colon_vector)-2 do i=i+1 & lon_start_index=i

   if colon_vector2[k]+dlon lt max(colon_vector) then begin
      while colon_vector[i] lt colon_vector2[k]+dlon and i lt n_elements(colon_vector)-2 do i=i+1 & lon_end_index=i
   endif else begin
      while colon_vector[i] lt max(colon_vector) and i lt n_elements(colon_vector)-2 do i=i+1 & lon_end_index=i
   endelse              ; if colon_vector2[k]+dlon lt max(colon_vector)

   ; ------------------------ ;
   tot_column_count = total(float(total_column_count[lat_start_index:lat_end_index, lon_start_index:lon_end_index, *]))
   
   if tot_column_count ne 0 then begin
   
      tot_ctz_dz1_grid = total(float( cbz_dz1_grid[lat_start_index:lat_end_index,                     $
                                            lon_start_index:lon_end_index, *,                              $
                                            base_ind1:base_ind2, dz_ind1:dz_ind2] ))       $
                         + total(float(ctz_dz1_grid[lat_start_index:lat_end_index,          $
                                                  lon_start_index:lon_end_index, *,       $
                                                  base_ind2:base_ind1, dz_ind1:dz_ind2]))
                                                               
      tot_ctz_dz2_grid = total(float(cbz_dz2_grid[lat_start_index:lat_end_index,                      $
                                           lon_start_index:lon_end_index, *,                               $
                                           0:base_ind2, *, base_ind1:base_ind2, dz_ind1:dz_ind2]))       $
                         + total( float(ctz_dz2_grid[lat_start_index:lat_end_index,   $
                                                   lon_start_index:lon_end_index, *,  $
                                                   base_ind2:n_elements(cbtz_bin_mid)-1, *,   $
                                                   base_ind2:base_ind1, dz_ind1:dz_ind2]) )
                                                   
      ; if_mult = 1 gives all the tops including multiple layers
      ; if_mult = 0 gives only the highest tops
         if_mult=0  
      if if_mult eq 1 then begin                                                                        
         layer_base_map[k,j] = ( tot_ctz_dz1_grid + tot_ctz_dz2_grid ) / tot_column_count
      endif else begin
         layer_base_map[k,j] = tot_ctz_dz1_grid  / tot_column_count
      endelse
      
      if k eq n_elements(colon_vector2)-1 then begin
         ;print, tot_ctz_dz1_grid, tot_column_count ,lon_start_index, lon_end_index
         layer_base_map[k,j]=(layer_base_map[k,j]+layer_base_map[0,j])/2.
      endif

   endif else begin
      layer_base_map[k,j]=0.0
   endelse         ; if tot_column_count ne 0


   endfor     ; for k=0,n_elements(colon_vector2)-1
   endfor     ; for j=0,n_elements(lat_vector2)-1

   ; now plot the map

   lat_avg_string=strtrim(string(2.*dlat,format='(f3.1)'))
   lon_avg_string=strtrim(string(2.*dlon,format='(f3.1)'))
   colorbar_title='Coverage'

   min_base = cbtz_bin_mid[base_ind1]-cbtz_dbin[base_ind1]
   max_base = cbtz_bin_mid[base_ind2]+cbtz_dbin[base_ind2]
   min_dz = dz_bin_mid[dz_ind1]-dz_dbin[dz_ind1]
   max_dz = dz_bin_mid[dz_ind2]+dz_dbin[dz_ind2]

   title_string='CloudSat/Calipso Hydrometeor Coverage (Normalization: All Profiles) . Avg Box: ' + $
	            lat_avg_string+'X'+lon_avg_string+'. For Period '+   $
	            '_'+datelist[0]+'-'+datelist[n_elements(datelist)-1]

   title_string2='Bases ' + strtrim(string(fix(min_base),format='(i6)'),2) + '-' $
	                      + strtrim(string(fix(max_base),format='(i6)'),2) +     $
	       ' Thickness: ' + strtrim(string(fix(min_dz),format='(i6)'),2)   +     $
	                  '-' + strtrim(string(fix(max_dz),format='(i6)'),2)

	if base_ind1 eq 0 and base_ind2 eq base_ind1 then begin
		spawn, 'mkdir ' + output_path_base
		;spawn, 'mkdir '+output_path+'/World'+'_'+datelist[0]+'-'+datelist[n_elements(datelist)-1]
	endif

	output_fname=output_path_base +   $
	      '/RLGeoprof_base_AbsFraction-World'$
	      +'_'+datelist[0]+'-'+datelist[n_elements(datelist)-1]+'_'+$
	      lat_avg_string+'X'+lon_avg_string+'_base_'+strtrim(string(fix(min_base),format='(i6)'),2)+'-'+$
	      strtrim(string(fix(max_base),format='(i6)'),2)+'_dz_'+$
	      strtrim(string(fix(min_dz),format='(i6)'),2)+'-'+$
	      strtrim(string(fix(max_dz),format='(i6)'),2)+'dBZ-All'

   dmax=max(layer_base_map)+0.1*(max(layer_base_map))
   dmin=min(layer_base_map)

   if dmax gt 0. and dmax ne dmin then begin

      plot_cloudsat_coverage_map_world, layer_base_map, 1000, 600, 5, 2.5, max_lat, min_lat, comax_lon, comin_lon, $
         0.1, 0.15, 0.9, 0.85, 0, 0, -70, 0, 70, 360, dmax, dmin, title_string, $
         lon_vector2,lat_vector2, title_string2, colorbar_title, 1.0,$
         output_fname+'.gif'

      write_2d_ncdf_map, output_fname+'.cdf',$
         layer_base_map,lon_vector2,lat_vector2,(comax_lon+comin_lon)/2.0,min_lat,$
         max_lon,max_lat,min_lon,comax_lon,comin_lon

   endif else begin
      print, 'no data '         ;,min_base, max_base, min_dz, max_dz
   endelse                      ; if dmax gt 0. and dmax ne dmin


endfor	; dz_ind2
endfor	; dz_ind1
endfor	; cbztz_ind2
endfor	; cbztz_ind1


return

end