;==================================================================
; IDL procedure history
; ---------------------
;
; January 2008  : Author : Eric Tromeur/CCNY/NOAA/CREST
;                 Numbers of Clusters, min lat, max lat, number of latitude in data
;                 and legend names as input of the procedure
;==================================================================

;==================================================================
; ISSUE
; -----
; - This procedure makes histograms showing the frequency of distribution
;   of each cluster for each lag and MJO events
;==================================================================

;==================================================================
; USAGE
; -----
; clushist_MJO_sublat_subIndex_bars_pentad,nclus,nlat,sublatmin,sublatmax,nyearmin,nyearmax,region,loc,directory
;==================================================================

;==================================================================
; INPUT
; -----
; nclus      : number of clusters (ex: 8)
; nlat       : number of latitude in the data (maxlat-minlat)/2.5 (ex: 28)
; sublatmin  : minimum of latitude in the data (ex: 18)
; sublatmax  : maximum of latitude in the data (ex: 28)
; nyearmin   : starts at year min (ex: 1997)
; nyearmax   : ends at year max (ex: 2000)
; region     : legend name related to the region (ex: 'Subtropics')
; loc        : legend name related to the location (ex: 'ST')
; directory  : path of "maps" directory (ex: 'Subtropics/maps')

; nota: loc
; ----
; Tropics    : loc='TR'
; SubTropics : loc='ST'
; MidlatN    : loc='MN'
; MidlatS    : loc='MS'

; nota: region
; ----
; Tropics     : region='Tropics'
; SubTropics  : region='Subtropics'
; MidlatN     : region='MidlatN'
; MidlatS     : region='MidlatS'
;==================================================================

;*****************************************************************************

pro clushist_MJO_sublat_subIndex_bars_pentad,nclus,nlat,sublatmin,sublatmax,nyearmin,nyearmax,region,loc,directory

;*****************************************************************************

; ------------------------------------
; Open files'name
; ------------------------------------
;namefile='../MJO_CHI200_index_ISCCPD1_JanDec8304.dat'
;namefile='../MJO_CHI200_index_ISCCPD1_MayOct8304.dat'
namefile='../MJO_CHI200_index_ISCCPD1_NovApr8304.dat'
;namefile='../MJO_CHI200_index_DelGenio.dat'
; ------------------------------------

;--------------------
; Initialize
;--------------------
maxlon=144
nlag=7
;--------------------
;rows=8
;rows=1568	;; column number for Jan-Dec period
;rows=801	;; column number for May-Oct period
rows=756	;; column number for Nov-Apr period
min_val=-2.2
max_val=-10.
;--------------------

; ------------------------------------
; Open files'name
; ------------------------------------
openr,1,namefile
; ------------------------------------

; ---------------------
; Define graph variable
; ---------------------
!p.multi=0
!p.charsize=1
!p.charthick=4
; ---------------------

;--------------------
; Initialize
;--------------------
index = FLTARR(rows,10)
index_tmp = FLTARR(rows,6)
strdate = ULONARR(1)
i1 = FLTARR(1)
i2 = FLTARR(1)
i3 = FLTARR(1)
i4 = FLTARR(1)
i5 = FLTARR(1)
i6 = FLTARR(1)
i7 = FLTARR(1)
i8 = FLTARR(1)
i9 = FLTARR(1)
i10 = FLTARR(1)
;--------------------
sz=FLTARR(nclus+1,nlag)
tsz=FLTARR(nlag)
rfo=FLTARR(nclus+1,nlag)
;--------------------

; ------------------------------------
; Use for input and output files'names
; ------------------------------------
if (nclus lt 10) then begin
 nclust=string(nclus,format='(i1)')
 nclust2='0'+string(nclus,format='(i1)')
endif else begin
 nclust=string(nclus,format='(i2)')
 nclust2=string(nclus,format='(i2)')
endelse
; ------------------------------------
ymin=string(nyearmin,format='(i4)')
ymax=string(nyearmax,format='(i4)')
; ------------------------------------

days=[31,28,31,30,31,30,31,31,30,31,30,31]
; ------------------------------------

;=============================================================================================
;                                       LOOP BEGINNING
;=============================================================================================

; -----------------
;  code for "lag"
; -----------------
; lag 3  : rfo(n,6)
; lag 2  : rfo(n,5)
; lag 1  : rfo(n,4)
; lag 0  : rfo(n,3)
; lag -1 : rfo(n,2)
; lag -2 : rfo(n,1)
; lag -3 : rfo(n,0)
; -----------------

mjocount=0

;#####################################################################################

FOR iloop=0,rows-1 DO BEGIN  	;; ---------------------------------------------------

;#####################################################################################

print,'************************'
print,'data row # :',iloop
print,'************************'
print,' '

 READF, 1, strdate, i9, i10, i1, i2, i3, i4, i5, i6, i7, i8
; print,'iloop=',iloop,strdate, i9, i10, i1, i2, i3, i4, i5, i6, i7, i8

index(iloop,0)=i9
index(iloop,1)=i10
index(iloop,2)=i1
index(iloop,3)=i2
index(iloop,4)=i3
index(iloop,5)=i4
index(iloop,6)=i5
index(iloop,7)=i6
index(iloop,8)=i7
index(iloop,9)=i8

 FOR i=1,6 DO BEGIN
     index_tmp(iloop,i-1)=index(iloop,i)
 ENDFOR

min=min_val
 max=max_val
 ind=-999
 IF(MIN(index_tmp(iloop,*)) gt max) THEN BEGIN
   FOR i=1,6 DO BEGIN
    IF(index(iloop,i) le min) THEN BEGIN
      min=index(iloop,i)
      ind=i
    ENDIF
   ENDFOR
 ENDIF

IF(ind ne -999) THEN BEGIN		;; -------------------------------------------

;=====================================================================================
strdate2 = strcompress(strdate,/remove)
yr=strmid(strdate2,0,4)
m=strmid(strdate2,4,2)-1
d=strmid(strdate2,6,2)-1

if(yr eq 1984) then days[1]=29
if(yr eq 1988) then days[1]=29
if(yr eq 1992) then days[1]=29
if(yr eq 1996) then days[1]=29
if(yr eq 2000) then days[1]=29
if(yr eq 2004) then days[1]=29

;=====================
; 160E : pos=64, ind=6
; 140E : pos=56, ind=5
; 120E : pos=48, ind=4
; 100E : pos=40, ind=3
; 80E  : pos=32, ind=2
; 70E  : pos=28, ind=1
;=====================
CASE ind OF
1: i=28-1
2: i=32-1
3: i=40-1
4: i=48-1
5: i=56-1
6: i=64-1
ENDCASE

clusarr_MJO=FLTARR(nlat,280)

dayref=0
; ----------------------------------
; 1 - 136 : D-17 to D (lag-30 to lag0)
; ----------------------------------
FOR idate=17,0,-1 DO BEGIN		;; idate loop --------------------------------------------

dlag=d+1
mlag=m
year=yr

   dlag=dlag-idate
   IF(dlag lt 1) THEN BEGIN
     mlag=mlag-1
     IF(mlag lt 0) THEN BEGIN
       year=year-1
       if(year eq 1984) then days[1]=29
       if(year eq 1988) then days[1]=29
       if(year eq 1992) then days[1]=29
       if(year eq 1996) then days[1]=29
       if(year eq 2000) then days[1]=29
       if(year eq 2004) then days[1]=29

       mlag=11
     ENDIF
     dlag=dlag+days[mlag]
   ENDIF

   day=dlag-1

   yyyy=string(year,format='(i4)')
   inputfile='../../../classify/'+region+'/'+nclust+'clus/clus'+nclust2+'.'+loc+'.'+yyyy+'.bin'
   clusarr=read_binary(inputfile,data_dims=[maxlon,nlat,12,31,8],data_type=3)

   FOR j=sublatmin-1,sublatmax-1 DO BEGIN           ;; latitude loop --------------------
    FOR h=0,7 DO BEGIN 
      clusarr_MJO[j,dayref+h]=clusarr[i,j,mlag,day,h] 
    ENDFOR
   ENDFOR
 
dayref=dayref+8

ENDFOR 		;; idate loop ----------------------------------------------------------


dayref=18*8
; --------------------------------------
; 137 - 280 : D+1 to D+17 (lag1 to lag30)
; --------------------------------------
FOR idate=1,17 DO BEGIN		;; idate loop --------------------------------------------

dlag=d+1
mlag=m
year=yr

   dlag=dlag+idate
   IF(dlag gt days[mlag]) THEN BEGIN
     dlag=dlag-days[mlag]
     mlag=mlag+1
     IF(mlag gt 11) THEN BEGIN
       year=year+1
       if(year eq 1984) then days[1]=29
       if(year eq 1988) then days[1]=29
       if(year eq 1992) then days[1]=29
       if(year eq 1996) then days[1]=29
       if(year eq 2000) then days[1]=29
       if(year eq 2004) then days[1]=29

       mlag=mlag-11-1
     ENDIF
   ENDIF

   day=dlag-1

   yyyy=string(year,format='(i4)')
   inputfile='../../../classify/'+region+'/'+nclust+'clus/clus'+nclust2+'.'+loc+'.'+yyyy+'.bin'
   clusarr=read_binary(inputfile,data_dims=[maxlon,nlat,12,31,8],data_type=3)

   FOR j=sublatmin-1,sublatmax-1 DO BEGIN           ;; latitude loop --------------------
    FOR h=0,7 DO BEGIN 
      clusarr_MJO[j,dayref+h]=clusarr[i,j,mlag,day,h] 
    ENDFOR
   ENDFOR
 
dayref=dayref+8

ENDFOR 		;; idate loop ----------------------------------------------------------


;=====================================================================================
;		LAG -1 / LAG -2 / LAG -3/ LAG 0 / LAG 1 / LAG 2 / LAG 3
;=====================================================================================

;================
; Case 0 : lag-3
; Case 1 : lag-2
; Case 2 : lag-1
; Case 3 : lag 0
; Case 4 : lag 1
; Case 5 : lag 2
; Case 6 : lag 3
;================

FOR ilag=0,nlag-1 DO BEGIN           ;; ilag loop -----------------------------------------

 CASE ilag OF
 0: hmin=0
 1: hmin=40
 2: hmin=80
 3: hmin=120
 4: hmin=160
 5: hmin=200
 6: hmin=240
 ENDCASE

 CASE ilag OF
 0: hmax=39
 1: hmax=79
 2: hmax=119
 3: hmax=159
 4: hmax=199
 5: hmax=239
 6: hmax=279
 ENDCASE

  FOR j=sublatmin-1,sublatmax-1 DO BEGIN           ;; latitude loop --------------------
    FOR h=hmin,hmax DO BEGIN 

      w=where(clusarr_MJO[j,h] eq 0)
      nb1=n_elements(w)
      w=w(0)
      if(w eq -1) then nb1=0
      sz(nclus,ilag)=sz(nclus,ilag)+nb1

      FOR n=0,nclus-1 DO BEGIN
        w=where(clusarr_MJO[j,h] eq n+1)
        nb1=n_elements(w)
        w=w(0)
        if(w eq -1) then nb1=0
        sz(n,ilag)=sz(n,ilag)+nb1
      ENDFOR

    ENDFOR
  ENDFOR            ;; latitude loop (j) -------------------------------

ENDFOR 		;; ilag loop --------------------------------------------------------

;=====================================================================================
   

;#####################################################################################

mjocount=mjocount+1

ENDIF		;; ind loop --------------------------------------------

ENDFOR 		;; iloop loop --------------------------------------------------------

CLOSE,1
;#####################################################################################


;=====================================================================================
;			RELATIVE FREQUENCY OF OCCURENCE
;=====================================================================================

for ilag=0,nlag-1 do begin
 for n=0,nclus do begin
    tsz(ilag)=tsz(ilag)+sz(n,ilag)
 endfor
endfor

for ilag=0,nlag-1 do begin
 for n=0,nclus do begin
    rfo(n,ilag)=sz(n,ilag)/tsz(ilag)*100
    if(sz(n,ilag) eq 0) then rfo(n,ilag)=0
    print,sz(n,ilag),tsz(ilag),rfo(n,ilag)
 endfor
endfor

;=====================================================================================


;=====================================================================================
;               	Plot cluster historgram per lag
;=====================================================================================
!P.CHARSIZE =2
!P.THICK =5
!X.THICK =5
!Y.THICK =5

set_plot,'PS'

fname=directory+'/Hist.'+loc+'.'+string(nclus,format='(i1)')+'.MJO.ps'
;device,/color,filename=fname,bits_per_pixel=8,/portrait
device,/color,filename=fname,bits_per_pixel=8,/landscape

;Make axes black:
;---------------
;tvlct, [0], [0], [0],0
;!P.COLOR=0

;Create a 2D array, equal in size to array, that has identical
;color index values across each row to ensure that the same item is
;represented by the same color in all bars:
;colors = INTARR(nclus+1,lag)
;FOR I = 0, 6 DO colors[*,I]=(20*I)+20

colors = INTARR(nclus+1,7)
;------------------------------
;       Resulting colors
;------------------------------
;     100:    Red
;     101:    Orange
;     102:    Yellow
;     103:    Green
;     104:    Cyan
;     105:    Blue
;       7:    Dark Blue
;     254:    Black
;     255:    White
;------------------------------
;Red
;------
TVLCT, 255, 0, 0, 100
;Orange
;------
TVLCT, 255, 150, 0, 101
;Yellow
;------
TVLCT, 255, 255, 0, 102
;Green
;------
TVLCT, 0, 255, 0, 103
;Cyan
;------
TVLCT, 0, 255, 255, 104
;Blue
;------
TVLCT, 0, 0, 255, 105

FOR I = 0, nlag-1 DO BEGIN
   colors[0,I]=100
   colors[1,I]=101
   colors[2,I]=102
   colors[3,I]=103
   colors[4,I]=104
   colors[5,I]=105
   colors[6,I]=7
ENDFOR
;------------------------------

;Scale range to accommodate the total bar lengths:
;-------------------------------------------------
;!Y.RANGE = [0, MAX(rfo)]
!Y.RANGE = [0, 60]
nrows = N_ELEMENTS(rfo[0,*])
ncols = N_ELEMENTS(rfo[*,0])

FOR I = 0,nrows-1 DO BEGIN
   lagcount=I-3
   lag=[' ',' ',string(lagcount,format='(i3)'),' ',' ',' ',' ']
;   BAR_PLOT, rfo[*,I], BACKGROUND=255, $
;      BARWIDTH=0.99, BARSPACE=0.2, BAROFFSET=I*(1.3*ncols), $
;      OVER=(I GT 0), BASERANGE=0.112,/OUTLINE,XTITLE='Lag in pentads',YTITLE='RFO (%)',TITLE='RFO of each cloud regime in 60E-180E region (MJO events: index<-1)',COLORS=colors[*,I],BARNAMES=lag
   BAR_PLOT, rfo[*,I], BACKGROUND=255, $
      BARWIDTH=0.99, BARSPACE=0.2, BAROFFSET=I*(1.3*ncols), $
      OVER=(I GT 0), BASERANGE=0.112,/OUTLINE,XTITLE='Lag in pentads',YTITLE='RFO (%)',TITLE='b/',COLORS=colors[*,I],BARNAMES=lag
ENDFOR
;=====================================================================================


;=====================================================================================
;                               	LEGEND
;=====================================================================================
;Load color table:
;----------------
LOADCT, 33, ncolors=7

;Make axes black:
;---------------
tvlct, [0], [0], [0],0
!P.COLOR=0

; vertical
;colorbar,range=[0.5,7.5],vertical=1,top=1,divisions=7,ncolors=254,position=[0.93, 0.20, 0.95, 0.90],format='(f4.1)'

;horizontal
;colorbar,range=[0.5,7.5],top=1,ncolors=7,position=[0.87, 0.35, 0.89, 0.75],divisions=7,format='(f4.1)',/reverse
colorbar,range=[0.5,7.5],top=1,ncolors=7,position=[0.73, 0.35, 0.75, 0.75],divisions=7,format='(f4.1)',/reverse,charsize=1.3,font=0
;colorbar,range=[0.5,7.5],top=1,ncolors=254,position=[0.87, 0.35, 0.89, 0.75],divisions=7,format='(f4.1)'

;=====================================================================================
text=string(mjocount,format='("MJO events : ",i3)')
print,text
text2=+ymin+'-'+ymax+' Period'

legend,[text],charsize=1.6,charthick=6,box=0,/right
legend,[text2],charsize=1.6,charthick=6,box=0,/left
;legend,[text],textcolors=!d.table_size-1,charsize=1,charthick=6,box=0,/right
; ==================================================================

device,/close

;=====================================================================================


END