;test_land_1.ncl
;=======================================================================
; Creat a function "unique_index" by Dennis Shea
;========================================================================
function unique_index (nMax[1], nUnique[1]:integer, nGenerate[1]:integer)
local nRandom, ratio, IR, LU, n, i, ni, iu, niu

begin
if (nUnique.gt.nMax) then
       print("-------")
       print("unique_subscripts: nUnique > nMax ; impossible")
       print("-------")
       return(-999)
   end if

; to ensure that enough random integers are generated check the ratio

   nRandom = nMax
   ratio = nRandom/nUnique
   if (ratio.lt.5 .or. nGenerate.gt.5) then
       nRandom = nRandom*max( (/5,nGenerate/) )
   end if

   IR  = round(random_uniform(0,(nMax-1),nRandom), 3)

   LU  = new( nRandom, "logical")
   LU  = True
                     ; this *will* be slow if nRandom is large
   do n=0,nRandom-1
      i   = ind(IR.eq.IR(n) .and. LU(n))
      ni  = dimsizes(i)
      if (ni.gt.1) then
         LU(i(1:)) = False
      end if
      delete(i)
   end do

   iu  = ind(LU)                       ; index of unique subscripts
   niu = dimsizes(iu)                  ; number of unique subscripts
   if (niu.ge.nUnique) then
       return(IR(iu(0:nUnique-1)))     ;
   else
       print("-------")
       print("unique_subscripts: no. of unique subscripts < NR ; increase nGenerate")
       print("-------")
       return(IR(iu(0:niu-1)))
   end if
printVarSummary(iu)
end
;
;====================================================================================
; Original Program
;====================================================================================

;*******************************************************
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
;===================================================================
; Assume variables T, PS and ORO exist and that they have coordinate
; variables time, lev, lat, lon.
;===================================================================
 f    = addfile("wrfout_july_all.nc","r")
 e    = addfile("wrfout_LANDMASK.nc", "r")
 T   = f->T(:,:,:,:)                 ;T(30,6,70,150)
 oro = e->LM(15,:,:)             ; Extraction of LANDMASK(70,150) 
 lat = f->lat                    ;latitude of WRF
 lon = f->lon                    ;longitude of WRF
 printVarSummary(T)
 ;======================================================================
T!2 ="lat"
T&lat=lat
T!3="lon"
T&lon=lon
;=======================================================================


oro!0 ="lat"
oro&lat=lat
oro!1="lon"
oro&lon=lon
;=========================================================================
T@_FillValue =  -9999.         ; sets _FillValue to -9999

TLand = T
TLand  =  mask(T,conform(T, oro, (/2,3/)), 1)
TLand@_FillValue =  -9999.
n= dimsizes(TLand)
printVarSummary(TLand)

;print(n)
T1D= new((/31,6,10500/),float)

do a= 0,30
   do b=0,5
   T1D(a,b,:)=ndtooned(TLand(a,b,:,:))
   end do
end do

T1D@_FillValue =  -9999.

N=num(ismissing(T1D(0,0,:)))

ii= ind(ismissing(T1D(0,0,:)))
print(ii)

n= dimsizes(T1D(0,0,:))

iT= unique_index (10500,100,0) ; 70X150=10500



print(iT)
TR= new((/31,6,100/),float)
do a = 0,30
   do b=0,5
   TR(a,b,:) = T1D(a,b, iT )        ; randomly sampled T
   end if
  end do
end do
TR@_FillValue =  -9999.         ; sets _FillValue to -9999
printVarSummary(TR)
TR_avg= new((/31,6/),float)

;====================================================================
ir = ind_resolve(iT,(/70,150/) )

wks = gsn_open_wks("ps","land_100")
res = True
res@gsnDraw = False
res@gsnFrame = False
res@mpMinLonF            =23.               ; select a subregion
res@mpMaxLonF            = 156.
res@mpMinLatF            = 40.
res@mpMaxLatF            =  72.
;res@gsnCenterString = "100 grid points "
plot = gsn_csm_map_ce(wks,res)

xres = True
xres@gsMarkerIndex = 1
xres@gsMarkerColor = "blue"
xres@gsMarkerThicknessF = 5.0
pmarker = gsn_add_polymarker(wks,plot,lon(ir(:,1)),lat(ir(:,0)),xres)
draw(plot)
frame(wks)
;=====================================================================
end