load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" 
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" 
                   
  begin 
    
;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
;                    Read in a series of daily SST netCDF files
;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      infile = addfile("/media/amadou/COULIBALY/Toshiba/wind_seabreeze/netcdf/all_new.nc", "r")
         
      temp = infile->air_temperature
      time = infile->time

      ; read in the station IDs
      stations = tostring(infile->station_id)
      nstat = dimsizes(stations)

      time@units = "seconds since 1970-01-01 00:00:00"
      dates = cd_calendar(time, 0)

      iSize = dimsizes(time)

      T1 = new((/12,nstat/), float)

      do i = 0, nstat-1

           T = temp(:,i)

           index=ind(dates(:,3) .ge. 6 .and. dates(:,3) .le. 18 .and. T(:) .ne. T@_FillValue)

                 if (.not.any(ismissing(index)) .and. dimsizes(index) .ge. 2) then
                 
                    T0=T(index)

                    dates0=dates(index,0:3) 

                    do j = 0, 11
                        
                         m=j+1
  
                         index1 = ind(dates0(:,1) .eq. m .and. dates0(:,0) .eq. 2006)

                         if (.not.any(ismissing(index1))) then

		                 T1(j,i)=dim_avg(T0(index1))
                         end if
		       
		                 delete(index1)                  

                    end do
                    
                    delete(dates0)                 
                    delete(T0)
                  
                end if
               
                delete(index)
                delete(T)
         end do
    
        T2 = dim_avg(T1)             ; mean monthly variation of temperature for each station
      ; print(T1)

;%%%%%%%%%%%%%%%%%%%%%%%%%%  Monthly variations of  daily Sea_Surface Temperature %%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      infile1 = addfile("/media/amadou/COULIBALY/Toshiba/SST_DATA/sst_2008/amsr-avhrr-v2.2008_mean.nc", "r")

          SST = infile1->sst
       
          SST0 = SST(:,0,0,0)*0.01
          
          SST1 = SST0 + 273.15
    
          Time = infile1->time
          time@units = "days since 1978-01-01 00:00:00"
          date = cd_calendar(Time, 0)

          SST01 = new((/12/), float)         
      
                     do k = 0, 11
                        
                         n=k+1
  
                         index2 = ind(date(:,1) .eq. n)

		         SST01(k)=dim_avg(SST1(index2))
                 
		         delete(index2)                  

                    end do
;print(SST01)
    Arr = new((/6,12/),float)
    
   ;Arr(0,:) = T1(:,0)
   ; Arr(1,:) = T1(:,7)
   ; Arr(2,:) = T1(:,12)
    ;Arr(3,:) = T1(:,23)
    Arr(4,:) = T1(:,25)
    Arr(5,:) = SST01

;%%%%%%%%%%%%%%%%%%%%%%%%%%  Monthly variations of temperature %%%%%%%%%%%%%%%%%%%%%%%%%%%%%

          x = ispan(0,11,1)                       
          n = 12
;******************************************************************************************************
;                 create the plot
;******************************************************************************************************  
	  wks = gsn_open_wks("pdf","Monthly_Var_Daytime_Temp_Ikeja_2008")  ; open pdf file

	  res                       = True                    ; plot mods desired         
	  res@gsnFrame              = False                   ; don't advance frame yet
          res@gsnDraw               = False 
          res@xyLineThicknesses     = (/14.5,16.0/)             ; make 2nd lines thicker
          res@xyLineColors          = (/"red","black"/)        ; change line color

          res@xyDashPattern         = 0                       ; Make curves all solid
	  res@xyMarkLineMode        = "MarkLines"             ; Markers *and* lines
	  res@xyMarkers             = (/0,1/)                ; 2 different markers
	  res@xyMarkerColors        = (/"blue","green"/)      ; 2 different colors
          res@tiMainString          = "Monthly variation Air_temp .vs SST(Ikeja_2008)"  ; add title
          
	  res@gsnYRefLine           = 0.0                    ; draw a reference line
	  res@tmXBOn                = True                    ; turn off tickmarks at bot
	  res@trYMinF               = 297                    
          res@trYMaxF               = 307
	  res@trXMinF               = -1                      ; adds space on either end
	  res@trXMaxF               = 12                      ; of the 1st and last bars
	  res@tmXBMode              = "Explicit"              ; explicit labels
	  res@tmXBValues            = fspan(0,11,n)
	  res@tmXBLabels            = (/"J","F","M","A","M","J","J","A","S","O","N","D"/)
          res@tiYAxisString         = "Temperature(K)"
          res@tiXAxisString         = "Month"
          res@tiXAxisFontHeightF    = 0.015 
          res@tiYAxisFontHeightF    = 0.015 

         plot = gsn_csm_xy(wks,x,Arr,res) 

	; Set up resources for a customized legend.

	  lgres                    = True
	  lgres@lgLineColors     = (/"red","black"/)
	  lgres@lgItemType         = "Lines"        ; show lines only (default)
	  lgres@lgLabelFontHeightF = .08            ; legend label font thickness
	  lgres@vpWidthF           = 0.13           ; width of legend (NDC)
	  lgres@vpHeightF          = 0.10           ; height of legend (NDC)
	  lgres@lgPerimThicknessF  = 2.0            ; thicken the box perimeter
	  lgres@lgMonoDashIndex    = True
	  lgres@lgDashIndex        = 0

	  labels = (/"T_Ikeja","SST"/)
	  nitems = dimsizes(labels)                 ; number of legend items

	; Create legend
	  lbid = gsn_create_legend(wks,nitems,labels,lgres)

	; Use gsn_add_annotation to attach this legend to our existing plot.
	; This way, if we resize the plot, the legend will stay with the
	; plot and be resized automatically.
	
	  amres                  = True
	
	; Point (0,0) is the dead center of the plot. Point (0,.5) is center,
	; flush bottom. Point (0.5,0.5) is flush bottom, flush right.

	  amres@amJust           = "BottomRight"    ; Use bottom right corner of box
		                                    ; for determining its location.
	  amres@amParallelPosF   = 0.5              ; Move legend to right
	  amres@amOrthogonalPosF = 0.5              ; Move legend down.

	  annoid = gsn_add_annotation(plot,lbid,amres)  ; add legend to plot

	;**********************************
	; Now maximize the plot
	;**********************************   
	; Maximize these plots for PS or PDF output. If output is                    
	; to an X11 window or an NCGM file, then this routine doesn't                
	; do anything.                                                               
	  
	  psres = True                                                               
	  maximize_output(wks,psres)  ; calls draw and frame for you 

  end