;**************************************************************************************
; Script produced by Justin Glisan at Iowa State University for the purpose of cutting
; out the full SN extreme analysis box (Over North America). After this is completed,
; masking commands are used to mask the the ocean and sea ice from Pan-Arctic WRF data.
; The data that is left is strictly land-based precipitation.
;***************************************************************************************
;***********************************************
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/wrf/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
;load "/usr/local/ncl/lib/ncarg/nclscripts/csm/shea_util.ncl"
;load "/usr/local/ncl/lib/ncarg/nclscripts/csm/contributed.ncl"
;************************************************

begin

;************************************************
; Set some universal time and lat/lon. coords
;************************************************

  TimeChara = t->time
  DimTimeChara = dimsizes(TimeChara)
  nTime    = TimeChara
  nTimes = dimsizes(nTime) - 1

  lat = a->lat
  lon = a->lon

  res = True

;************************************************
; Read in precip. data
;************************************************

  wrf_mean     = a->Daily_Precip

  wrf_mean!0   = "time"
  wrf_mean!1   = "lat"
  wrf_mean!2   = "lon"

;************************************************
; Read seaice, land, and ocean flags
;************************************************

  seaice    = b->SeaIce  
  seaice!0  = "time"
  seaice!1  = "lat"
  seaice!2  = "lon"

  land      = b->LandMask
  land!0    = "time"
  land!1    = "lat"
  land!2    = "lon"

;************************************************
; Time average data
;************************************************

  land_mean   = dim_avg( land(lat|:, lon|:, time|:) )
  seaice_mean = dim_avg( seaice(lat|:, lon|:, time|:) )

;************************************************
; Here, the NA domain is selected. Within this 
; full domain, the four sub-boxes are found.
;************************************************

  latMin  =  58
  latMax  =  68

  lonMin = -168
  lonMax = -140

;************************************************
; Next, we slice the domain out of each field
;************************************************

  wrf_mean  = mask(wrf_mean,(wrf_mean@lat2d.ge.latMin .and. wrf_mean@lat2d.le.latMax  .and. \
                wrf_mean@lon2d.ge.lonMin .and. wrf_mean@lon2d.le.lonMax), True)


  seaice_mean = mask(seaice_mean,(seaice_mean@lat2d.ge.latMin .and. seaice_mean@lat2d.le.latMax  .and. \
                seaice_mean@lon2d.ge.lonMin .and. seaice_mean@lon2d.le.lonMax), True)

  land_mean = mask(land_mean,(land_mean@lat2d.ge.latMin .and. land_mean@lat2d.le.latMax  .and. \
                land_mean@lon2d.ge.lonMin .and. land_mean@lon2d.le.lonMax), True)
  
;************************************************
; use mask function to mask out seaice (seaice=0
; from ocean data and then ocean data from land
; ocean = 1
;************************************************

  seaice_only  = mask(wrf_mean,seaice_mean,0) 
  seaice_only@lat2d = lat    ; 2D coordinate arrays
  seaice_only@lon2d = lon
 
  land_mask = mask(seaice_only,land_mean,1)
  land_mask@lat2d = lat    ; 2D coordinate arrays
  land_mask@lon2d = lon

  land_mask!0 = "time"
  land_mask!1 = "lat"
  land_mask!2 = "lon"

;************************************************
; Here, print out data in a redirected file or 
; output an ASCII file. Land/seaice-masked regions
; should be denoted at missing values
;************************************************

  print(land_mask)

  format = "f6.2"
  opt = True
  opt@fout = "Alaska_North_Land_Region_Precip.txt"
  write_matrix(multi_year_var,format,opt)