; Calculate average historic NPP (1981-2010) for NP, CP, and SP west and east of the 100th meridian, and create extreme thresholds 25% above and below this average per region and west/east.

; ---------------------------------------------------

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"
load "~/NIFA/NCL_6.5.0_functions/simple_legend_ndc.ncl"

; ---------------------------------------------------

begin

RCP_list = [/"rcp_45","rcp_85"/]

	model_list = [/"CCSM4","GFDL-ESM2M","HadGEM2-ES365","IPSL-CM5A-LR","MRI-CGCM3"/]

RCP_title = [/"RCP 4.5","RCP 8.5"/]


do h = 0, ListCount(RCP_list) - 1 

	do k = 0, ListCount(model_list) - 1 

	
  ; *************************************************
  ; 1 - READ IN MC2 NPP FILES FOR NP, CP, SP
  ; *************************************************

	print("")
	print("1a - Read in historic reference")
	
  NPP_hist = addfile("~/NIFA/MC2/subset_lat_zones/extreme_NPP/hist/NPP_hist_nfs_extreme_frequency_intensity_zone_0-3_west_east.nc","r")

	print("1b - Read in future MC2 NPP files for NP, CP, and SP - " + RCP_list[h](0) + " " + model_list[k](0))
	
	diri = "~/NIFA/MC2/subset_lat_zones/" + RCP_list[h](0) + "/nfs/" + model_list[k](0) + "/NPP_3_GROUPED/"

  NPP_1_f = addfile(diri + RCP_list[h](0) + "_nfs_" + model_list[k](0) + "_NPP_3_grouped_lat_zone_1.nc","r")
  NPP_2_f = addfile(diri + RCP_list[h](0) + "_nfs_" + model_list[k](0) + "_NPP_3_grouped_lat_zone_2.nc","r")
  NPP_3_f = addfile(diri + RCP_list[h](0) + "_nfs_" + model_list[k](0) + "_NPP_3_grouped_lat_zone_3.nc","r")


	;*********************************************
	; 2 - READ IN PLOT VARIABLES
	;********************************************* 

  print("2 - Read in variables")

	lat_1 = NPP_1_f->lat
	lon_1_w = NPP_1_f->lon(0:119)
	lon_1_e = NPP_1_f->lon(121:240)
	lat_2 = NPP_2_f->lat
	lon_2_w = NPP_2_f->lon(0:119)
	lon_2_e = NPP_2_f->lon(121:240)
	lat_3 = NPP_3_f->lat
	lon_3_w = NPP_3_f->lon(0:119)
	lon_3_e = NPP_3_f->lon(121:240)

	year = NPP_1_f->year(26:50)


	;************************************************************
	; 3 - CALCULATE AVERAGE NPP FOR NP, CP, SP, AND WEST AND EAST
	;************************************************************

  print("3 - Calculate average NPP per grid cell and determine exceedance")

	; extract NPP variables
	NPP_ALL_1_w = (/NPP_1_f->NPP_ALL(26:50,:,0:119)/)				; total NPP for the NP for west and east of the 100th longitude
	NPP_ALL_1_w!0 = "year"
	NPP_ALL_1_w&year = year
	NPP_ALL_1_w!1 = "lat_1"
	NPP_ALL_1_w&lat_1 = lat_1
	NPP_ALL_1_w!2 = "lon_1_w"
	NPP_ALL_1_w&lon_1_w = lon_1_w
	NPP_ALL_1_e = (/NPP_1_f->NPP_ALL(26:50,:,121:240)/)
	NPP_ALL_1_e!0 = "year"
	NPP_ALL_1_e&year = year
	NPP_ALL_1_e!1 = "lat_1"
	NPP_ALL_1_e&lat_1 = lat_1
	NPP_ALL_1_e!2 = "lon_1_e"
	NPP_ALL_1_e&lon_1_e = lon_1_e
	NPP_ALL_2_w = (/NPP_2_f->NPP_ALL(26:50,:,0:119)/)				; total NPP for the CP for west and east of the 100th longitude
	NPP_ALL_2_w!0 = "year"
	NPP_ALL_2_w&year = year
	NPP_ALL_2_w!1 = "lat_2"
	NPP_ALL_2_w&lat_2 = lat_2
	NPP_ALL_2_w!2 = "lon_2_w"
	NPP_ALL_2_w&lon_2_w = lon_2_w
	NPP_ALL_2_e = (/NPP_2_f->NPP_ALL(26:50,:,121:240)/)
	NPP_ALL_2_e!0 = "year"
	NPP_ALL_2_e&year = year
	NPP_ALL_2_e!1 = "lat_2"
	NPP_ALL_2_e&lat_2 = lat_2
	NPP_ALL_2_e!2 = "lon_2_e"
	NPP_ALL_2_e&lon_2_e = lon_2_e
	NPP_ALL_3_w = (/NPP_3_f->NPP_ALL(26:50,:,0:119)/)				; total NPP for the SP for west and east of the 100th longitude
	NPP_ALL_3_w!0 = "year"
	NPP_ALL_3_w&year = year
	NPP_ALL_3_w!1 = "lat_3"
	NPP_ALL_3_w&lat_3 = lat_3
	NPP_ALL_3_w!2 = "lon_3_w"
	NPP_ALL_3_w&lon_3_w = lon_3_w
	NPP_ALL_3_e = (/NPP_3_f->NPP_ALL(26:50,:,121:240)/)
	NPP_ALL_3_e!0 = "year"
	NPP_ALL_3_e&year = year
	NPP_ALL_3_e!1 = "lat_3"
	NPP_ALL_3_e&lat_3 = lat_3
	NPP_ALL_3_e!2 = "lon_3_e"
	NPP_ALL_3_e&lon_3_e = lon_3_e

	; historic averages (1981-2010), used to determine if future cells are extreme or not
	NPP_1_w_avg = NPP_hist->NPP_1_w_avg
	NPP_1_e_avg = NPP_hist->NPP_1_e_avg
	NPP_2_w_avg = NPP_hist->NPP_2_w_avg
	NPP_2_e_avg = NPP_hist->NPP_2_e_avg
	NPP_3_w_avg = NPP_hist->NPP_3_w_avg
	NPP_3_e_avg = NPP_hist->NPP_3_e_avg

	; Define extreme gridcells in individual years
	NPP_1_w_low = new(dimsizes(NPP_ALL_1_w),float,9.96921E36)
	NPP_1_e_low = new(dimsizes(NPP_ALL_1_e),float,9.96921E36)
	NPP_2_w_low = new(dimsizes(NPP_ALL_2_w),float,9.96921E36)
	NPP_2_e_low = new(dimsizes(NPP_ALL_2_e),float,9.96921E36)
	NPP_3_w_low = new(dimsizes(NPP_ALL_3_w),float,9.96921E36)
	NPP_3_e_low = new(dimsizes(NPP_ALL_3_e),float,9.96921E36)

	do x = 0,dimsizes(year)-1		; low extreme NPP = value, otherwise = missing value
		NPP_1_w_low(x,:,:) = where((NPP_ALL_1_w(x,:,:) .lt. NPP_1_w_avg*0.75),NPP_ALL_1_w(x,:,:),9.96921E36)
		NPP_1_e_low(x,:,:) = where((NPP_ALL_1_e(x,:,:) .lt. NPP_1_e_avg*0.75),NPP_ALL_1_e(x,:,:),9.96921E36)
		NPP_2_w_low(x,:,:) = where((NPP_ALL_2_w(x,:,:) .lt. NPP_2_w_avg*0.75),NPP_ALL_2_w(x,:,:),9.96921E36)
		NPP_2_e_low(x,:,:) = where((NPP_ALL_2_e(x,:,:) .lt. NPP_2_e_avg*0.75),NPP_ALL_2_e(x,:,:),9.96921E36)
		NPP_3_w_low(x,:,:) = where((NPP_ALL_3_w(x,:,:) .lt. NPP_3_w_avg*0.75),NPP_ALL_3_w(x,:,:),9.96921E36)
		NPP_3_e_low(x,:,:) = where((NPP_ALL_3_e(x,:,:) .lt. NPP_3_e_avg*0.75),NPP_ALL_3_e(x,:,:),9.96921E36)
	end do

	; Define near normal gridcells in individual years
	NPP_1_w_hi = new(dimsizes(NPP_ALL_1_w),float,9.96921E36)
	NPP_1_e_hi = new(dimsizes(NPP_ALL_1_e),float,9.96921E36)
	NPP_2_w_hi = new(dimsizes(NPP_ALL_2_w),float,9.96921E36)
	NPP_2_e_hi = new(dimsizes(NPP_ALL_2_e),float,9.96921E36)
	NPP_3_w_hi = new(dimsizes(NPP_ALL_3_w),float,9.96921E36)
	NPP_3_e_hi = new(dimsizes(NPP_ALL_3_e),float,9.96921E36)

	do x = 0,dimsizes(year)-1		; high extreme NPP = value, extreme NPP = missing value
		NPP_1_w_hi(x,:,:) = where((NPP_ALL_1_w(x,:,:) .gt. NPP_1_w_avg*1.25),NPP_ALL_1_w(x,:,:),9.96921E36)
		NPP_1_e_hi(x,:,:) = where((NPP_ALL_1_e(x,:,:) .gt. NPP_1_e_avg*1.25),NPP_ALL_1_e(x,:,:),9.96921E36)
		NPP_2_w_hi(x,:,:) = where((NPP_ALL_2_w(x,:,:) .gt. NPP_2_w_avg*1.25),NPP_ALL_2_w(x,:,:),9.96921E36)
		NPP_2_e_hi(x,:,:) = where((NPP_ALL_2_e(x,:,:) .gt. NPP_2_e_avg*1.25),NPP_ALL_2_e(x,:,:),9.96921E36)
		NPP_3_w_hi(x,:,:) = where((NPP_ALL_3_w(x,:,:) .gt. NPP_3_w_avg*1.25),NPP_ALL_3_w(x,:,:),9.96921E36)
		NPP_3_e_hi(x,:,:) = where((NPP_ALL_3_e(x,:,:) .gt. NPP_3_e_avg*1.25),NPP_ALL_3_e(x,:,:),9.96921E36)
	end do


	;************************************************************
	; 4 - CALCULATE FREQUENCY PER GRID CELL AND REGION
	;************************************************************

  print("4 - Calculate frequency of extreme and average years per grid cell and region")

	NPP_1_w_low_y_n = NPP_1_w_low		; NPP_1_w_low is 3D matrix with low extreme values remaining, near normal values are missing values
	NPP_1_e_low_y_n = NPP_1_e_low
	NPP_2_w_low_y_n = NPP_2_w_low
	NPP_2_e_low_y_n = NPP_2_e_low
	NPP_3_w_low_y_n = NPP_3_w_low
	NPP_3_e_low_y_n = NPP_3_e_low

	NPP_1_w_hi_y_n = NPP_1_w_hi
	NPP_1_e_hi_y_n = NPP_1_e_hi
	NPP_2_w_hi_y_n = NPP_2_w_hi
	NPP_2_e_hi_y_n = NPP_2_e_hi
	NPP_3_w_hi_y_n = NPP_3_w_hi
	NPP_3_e_hi_y_n = NPP_3_e_hi

	; delete missing value attribute for where-function to work
	delete([/NPP_1_w_low_y_n@_FillValue,NPP_1_e_low_y_n@_FillValue,NPP_2_w_low_y_n@_FillValue,NPP_2_e_low_y_n@_FillValue,NPP_3_w_low_y_n@_FillValue,NPP_3_e_low_y_n@_FillValue/])
	delete([/NPP_1_w_hi_y_n@_FillValue,NPP_1_e_hi_y_n@_FillValue,NPP_2_w_hi_y_n@_FillValue,NPP_2_e_hi_y_n@_FillValue,NPP_3_w_hi_y_n@_FillValue,NPP_3_e_hi_y_n@_FillValue/])

	; Replace extreme values with 1's and everything else with zero
	NPP_1_w_low_y_n = where(NPP_1_w_low_y_n .ne. 9.96921E36,1,0)	; 3D matrix, low extreme year = 1, else = 0 for every year
	NPP_1_e_low_y_n = where(NPP_1_e_low_y_n .ne. 9.96921E36,1,0)
	NPP_2_w_low_y_n = where(NPP_2_w_low_y_n .ne. 9.96921E36,1,0)
	NPP_2_e_low_y_n = where(NPP_2_e_low_y_n .ne. 9.96921E36,1,0)
	NPP_3_w_low_y_n = where(NPP_3_w_low_y_n .ne. 9.96921E36,1,0)
	NPP_3_e_low_y_n = where(NPP_3_e_low_y_n .ne. 9.96921E36,1,0)

	NPP_1_w_hi_y_n = where(NPP_1_w_hi_y_n .ne. 9.96921E36,1,0)	; 3D matrix, high extreme year = 1, else = 0 for every year
	NPP_1_e_hi_y_n = where(NPP_1_e_hi_y_n .ne. 9.96921E36,1,0)
	NPP_2_w_hi_y_n = where(NPP_2_w_hi_y_n .ne. 9.96921E36,1,0)
	NPP_2_e_hi_y_n = where(NPP_2_e_hi_y_n .ne. 9.96921E36,1,0)
	NPP_3_w_hi_y_n = where(NPP_3_w_hi_y_n .ne. 9.96921E36,1,0)
	NPP_3_e_hi_y_n = where(NPP_3_e_hi_y_n .ne. 9.96921E36,1,0)


	; sum up low and high extreme years per grid cell and average per decade per grid cell, zone, and east/west
	
	NPP_1_w_low_freq = dim_sum_n(NPP_1_w_low_y_n,0)/3	; 2D matrix, add up 30 years per grid cell then divide by 3 to get average per decade
	NPP_1_e_low_freq = dim_sum_n(NPP_1_e_low_y_n,0)/3
	NPP_2_w_low_freq = dim_sum_n(NPP_2_w_low_y_n,0)/3
	NPP_2_e_low_freq = dim_sum_n(NPP_2_e_low_y_n,0)/3
	NPP_3_w_low_freq = dim_sum_n(NPP_3_w_low_y_n,0)/3
	NPP_3_e_low_freq = dim_sum_n(NPP_3_e_low_y_n,0)/3

	NPP_1_w_hi_freq = dim_sum_n(NPP_1_w_hi_y_n,0)/3	; 2D matrix, add up 30 years per grid cell then divide by 3 to get average per decade
	NPP_1_e_hi_freq = dim_sum_n(NPP_1_e_hi_y_n,0)/3
	NPP_2_w_hi_freq = dim_sum_n(NPP_2_w_hi_y_n,0)/3
	NPP_2_e_hi_freq = dim_sum_n(NPP_2_e_hi_y_n,0)/3
	NPP_3_w_hi_freq = dim_sum_n(NPP_3_w_hi_y_n,0)/3
	NPP_3_e_hi_freq = dim_sum_n(NPP_3_e_hi_y_n,0)/3

	; meta data
	NPP_1_w_low_freq@long_name = "Number of extreme years per decade per grid cell, NP, west, low extreme (<25% of average historic NPP)"
	NPP_1_e_low_freq@long_name = "Number of extreme years per decade per grid cell, NP, east, low extreme (<25% of average historic NPP)"
	NPP_2_w_low_freq@long_name = "Number of extreme years per decade per grid cell, CP, west, low extreme (<25% of average historic NPP)"
	NPP_2_e_low_freq@long_name = "Number of extreme years per decade per grid cell, CP, east, low extreme (<25% of average historic NPP)"
	NPP_3_w_low_freq@long_name = "Number of extreme years per decade per grid cell, SP, west, low extreme (<25% of average historic NPP)"
	NPP_3_e_low_freq@long_name = "Number of extreme years per decade per grid cell, SP, east, low extreme (<25% of average historic NPP)"
	
	NPP_1_w_hi_freq@long_name = "Number of extreme years per decade per grid cell, NP, west, high extreme (>25% of average historic NPP)"
	NPP_1_e_hi_freq@long_name = "Number of extreme years per decade per grid cell, NP, east, high extreme (>25% of average historic NPP)"
	NPP_2_w_hi_freq@long_name = "Number of extreme years per decade per grid cell, CP, west, high extreme (>25% of average historic NPP)"
	NPP_2_e_hi_freq@long_name = "Number of extreme years per decade per grid cell, CP, east, high extreme (>25% of average historic NPP)"
	NPP_3_w_hi_freq@long_name = "Number of extreme years per decade per grid cell, SP, west, high extreme (>25% of average historic NPP)"
	NPP_3_e_hi_freq@long_name = "Number of extreme years per decade per grid cell, SP, east, high extreme (>25% of average historic NPP)"


	; Average grid cell frequency per zone, west/east, and low/high
	
	NPP_low_hi_freq_avg = new((/3,2,2/),float,9.96921E36)

	; dimensions
;	NPP_low_hi_freq_avg!0 = "zone"
;	NPP_low_hi_freq_avg&zone = (/1,2,3/)
;	NPP_low_hi_freq_avg!1 = "west_east"
;	NPP_low_hi_freq_avg&west_east = (/1,2/)
;	NPP_low_hi_freq_avg!2 = "low_high"
;	NPP_low_hi_freq_avg&low_high = (/1,2/)
	
	; Northern Plains
	NPP_low_hi_freq_avg(0,0,0) = avg(NPP_1_w_low_freq)		; one 1981-2010 average for each zone, east/west, and low/high -> 12 numbers total
	NPP_low_hi_freq_avg(0,0,1) = avg(NPP_1_w_hi_freq)
	NPP_low_hi_freq_avg(0,1,0) = avg(NPP_1_e_low_freq)
	NPP_low_hi_freq_avg(0,1,1) = avg(NPP_1_e_hi_freq)
	; Central Plains
	NPP_low_hi_freq_avg(1,0,0) = avg(NPP_2_w_low_freq)
	NPP_low_hi_freq_avg(1,0,1) = avg(NPP_2_w_hi_freq)
	NPP_low_hi_freq_avg(1,1,0) = avg(NPP_2_e_low_freq)
	NPP_low_hi_freq_avg(1,1,1) = avg(NPP_2_e_hi_freq)
	; Southern Plains
	NPP_low_hi_freq_avg(2,0,0) = avg(NPP_3_w_low_freq)
	NPP_low_hi_freq_avg(2,0,1) = avg(NPP_3_w_hi_freq)
	NPP_low_hi_freq_avg(2,1,0) = avg(NPP_3_e_low_freq)
	NPP_low_hi_freq_avg(2,1,1) = avg(NPP_3_e_hi_freq)

	; meta data
	NPP_low_hi_freq_avg@long_name = "Mid-century (2041-2065) zonal averages for number of years per decade with NPP >±25% historic average NPP, RCP: non fire suppression, for: " + RCP_title[h](0) + ", Model: " + model_list[k](0)
	NPP_low_hi_freq_avg@comment = "Dimensions: NP/CP/SP zone, west/east of 100W, low/high extreme frequency per decade"


	;************************************************************
	; 5 - CALCULATE INTENSITY PER GRID CELL AND REGION
	;************************************************************

  print("5 - Calculate intensity of extreme and average years per grid cell and region")

	; Average per grid cell (1981-2010)
	; low extreme intensity
	NPP_1_w_low_int = dim_avg_n(NPP_1_w_low,0)	; 2D matrix, low extreme NPP for 1981-2010 averaged per grid cell
	NPP_1_e_low_int = dim_avg_n(NPP_1_e_low,0)
	NPP_2_w_low_int = dim_avg_n(NPP_2_w_low,0)
	NPP_2_e_low_int = dim_avg_n(NPP_2_e_low,0)
	NPP_3_w_low_int = dim_avg_n(NPP_3_w_low,0)
	NPP_3_e_low_int = dim_avg_n(NPP_3_e_low,0)

	; high extreme intensity
	NPP_1_w_hi_int = dim_avg_n(NPP_1_w_hi,0)		; 2D matrix, high extreme NPP for 1981-2010 averaged per grid cell
	NPP_1_e_hi_int = dim_avg_n(NPP_1_e_hi,0)
	NPP_2_w_hi_int = dim_avg_n(NPP_2_w_hi,0)
	NPP_2_e_hi_int = dim_avg_n(NPP_2_e_hi,0)
	NPP_3_w_hi_int = dim_avg_n(NPP_3_w_hi,0)
	NPP_3_e_hi_int = dim_avg_n(NPP_3_e_hi,0)

	; meta data
	NPP_1_w_low_int@long_name = "NPP intensity of extreme years per grid cell, NP, west, low extreme (<25% of average NPP)"
	NPP_1_e_low_int@long_name = "NPP intensity of extreme years per grid cell, NP, east, low extreme (<25% of average NPP)"
	NPP_2_w_low_int@long_name = "NPP intensity of extreme years per grid cell, CP, west, low extreme (<25% of average NPP)"
	NPP_2_e_low_int@long_name = "NPP intensity of extreme years per grid cell, CP, east, low extreme (<25% of average NPP)"
	NPP_3_w_low_int@long_name = "NPP intensity of extreme years per grid cell, SP, west, low extreme (<25% of average NPP)"
	NPP_3_e_low_int@long_name = "NPP intensity of extreme years per grid cell, SP, east, low extreme (<25% of average NPP)"
	
	NPP_1_w_hi_int@long_name = "NPP intensity of extreme years per grid cell, NP, west, high extreme (>25% of average NPP)"
	NPP_1_e_hi_int@long_name = "NPP intensity of extreme years per grid cell, NP, east, high extreme (>25% of average NPP)"
	NPP_2_w_hi_int@long_name = "NPP intensity of extreme years per grid cell, CP, west, high extreme (>25% of average NPP)"
	NPP_2_e_hi_int@long_name = "NPP intensity of extreme years per grid cell, CP, east, high extreme (>25% of average NPP)"
	NPP_3_w_hi_int@long_name = "NPP intensity of extreme years per grid cell, SP, west, high extreme (>25% of average NPP)"
	NPP_3_e_hi_int@long_name = "NPP intensity of extreme years per grid cell, SP, east, high extreme (>25% of average NPP)"

	
	; Average grid cell intensity per zone, west/east, and low/high
	
	NPP_low_hi_int_avg = new((/3,2,2/),float,9.96921E36)

	; dimensions
;	NPP_low_hi_int_avg!0 = "zone"
;	NPP_low_hi_int_avg&zone = (/1,2,3/)
;	NPP_low_hi_int_avg!1 = "west_east"
;	NPP_low_hi_int_avg&west_east = (/1,2/)
;	NPP_low_hi_int_avg!2 = "low_high"
;	NPP_low_hi_int_avg&low_high = (/1,2/)
	

	; Northern Plains
	NPP_low_hi_int_avg(0,0,0) = avg(NPP_1_w_low_int)		; one 1981-2010 average for each zone, east/west, and low/high -> 12 numbers total
	NPP_low_hi_int_avg(0,0,1) = avg(NPP_1_w_hi_int)
	NPP_low_hi_int_avg(0,1,0) = avg(NPP_1_e_low_int)
	NPP_low_hi_int_avg(0,1,1) = avg(NPP_1_e_hi_int)
	; Central Plains
	NPP_low_hi_int_avg(1,0,0) = avg(NPP_2_w_low_int)
	NPP_low_hi_int_avg(1,0,1) = avg(NPP_2_w_hi_int)
	NPP_low_hi_int_avg(1,1,0) = avg(NPP_2_e_low_int)
	NPP_low_hi_int_avg(1,1,1) = avg(NPP_2_e_hi_int)
	; Southern Plains
	NPP_low_hi_int_avg(2,0,0) = avg(NPP_3_w_low_int)
	NPP_low_hi_int_avg(2,0,1) = avg(NPP_3_w_hi_int)
	NPP_low_hi_int_avg(2,1,0) = avg(NPP_3_e_low_int)
	NPP_low_hi_int_avg(2,1,1) = avg(NPP_3_e_hi_int)
	
	; meta data
	NPP_low_hi_int_avg@long_name = "Mid-century (2041-2065) zonal averages for average NPP for years with NPP >±25% historic average NPP, RCP: non fire suppression, for: " + RCP_title[h](0) + ", Model: " + model_list[k](0)
	NPP_low_hi_int_avg@comment = "Dimensions: NP/CP/SP zone, west/east of 100W, low/high extreme intensity"


  ; ***********************************************
  ; 6 - WRITE DATA INTO NEW netCDF FILE
  ; ***********************************************

  print("6 - Write data into new netCDF file")

    netCDF1 = True						; Output format is NetCDF

  if (netCDF1) then
    diro = "~/NIFA/MC2/subset_lat_zones/extreme_NPP/future/"
    filo = RCP_list[h](0) + "_" + model_list[k](0) + "_NPP_future_extreme_frequency_intensity_zone_0-3_west_east.nc"

  end if

  setfileoption("nc","Format","LargeFile")

  if (netCDF1) then
    system("/bin/rm -f "+ diro + filo)				; removes old files with the same name if they are present
    ncdf = addfile(diro + filo,"c")

    setfileoption(ncdf,"DefineMode",True)

    ; create attributes
    fAtt = True
    fAtt@title = "Extreme NPP frequency and intensity for mid-century (2041-2065), extreme means >±25% of zonal historic average NPP (1981-2010), historic NPP simulated from PRISM data, future NPP based on: " + RCP_title[h](0) + ", non fire suppression, Model: " + model_list[k](0)
    fAtt@Conventions = "None"
    fAtt@creation_date = systemfunc("date")

    fileattdef(ncdf,fAtt)

    dimNames = (/"lat_1","lon_1_w","lon_1_e","lat_2","lon_2_w","lon_2_e","lat_3","lon_3_w","lon_3_e"/)
    dimSizes = (/dimsizes(lat_1),dimsizes(lon_1_w),dimsizes(lon_1_e),dimsizes(lat_2),dimsizes(lon_2_w),dimsizes(lon_2_e),dimsizes(lat_3),dimsizes(lon_3_w),dimsizes(lon_3_e)/)
    dimUnlim = (/False,False,False,False,False,False,False,False,False/)

    filedimdef(ncdf,dimNames,dimSizes,dimUnlim)

    filevardef(ncdf,"NPP_1_w_low_freq",typeof(NPP_1_w_low_freq),(/"lat_1","lon_1_w"/))
    filevardef(ncdf,"NPP_1_e_low_freq",typeof(NPP_1_e_low_freq),(/"lat_1","lon_1_e"/))
    filevardef(ncdf,"NPP_2_w_low_freq",typeof(NPP_2_w_low_freq),(/"lat_2","lon_2_w"/))
    filevardef(ncdf,"NPP_2_e_low_freq",typeof(NPP_2_e_low_freq),(/"lat_2","lon_2_e"/))
    filevardef(ncdf,"NPP_3_w_low_freq",typeof(NPP_3_w_low_freq),(/"lat_3","lon_3_w"/))
    filevardef(ncdf,"NPP_3_e_low_freq",typeof(NPP_3_e_low_freq),(/"lat_3","lon_3_e"/))

    filevardef(ncdf,"NPP_1_w_hi_freq",typeof(NPP_1_w_hi_freq),(/"lat_1","lon_1_w"/))
    filevardef(ncdf,"NPP_1_e_hi_freq",typeof(NPP_1_e_hi_freq),(/"lat_1","lon_1_e"/))
    filevardef(ncdf,"NPP_2_w_hi_freq",typeof(NPP_2_w_hi_freq),(/"lat_2","lon_2_w"/))
    filevardef(ncdf,"NPP_2_e_hi_freq",typeof(NPP_2_e_hi_freq),(/"lat_2","lon_2_e"/))
    filevardef(ncdf,"NPP_3_w_hi_freq",typeof(NPP_3_w_hi_freq),(/"lat_3","lon_3_w"/))
    filevardef(ncdf,"NPP_3_e_hi_freq",typeof(NPP_3_e_hi_freq),(/"lat_3","lon_3_e"/))

    filevardef(ncdf,"NPP_1_w_low_int",typeof(NPP_1_w_low_int),(/"lat_1","lon_1_w"/))
    filevardef(ncdf,"NPP_1_e_low_int",typeof(NPP_1_e_low_int),(/"lat_1","lon_1_e"/))
    filevardef(ncdf,"NPP_2_w_low_int",typeof(NPP_2_w_low_int),(/"lat_2","lon_2_w"/))
    filevardef(ncdf,"NPP_2_e_low_int",typeof(NPP_2_e_low_int),(/"lat_2","lon_2_e"/))
    filevardef(ncdf,"NPP_3_w_low_int",typeof(NPP_3_w_low_int),(/"lat_3","lon_3_w"/))
    filevardef(ncdf,"NPP_3_e_low_int",typeof(NPP_3_e_low_int),(/"lat_3","lon_3_e"/))

    filevardef(ncdf,"NPP_1_w_hi_int",typeof(NPP_1_w_hi_int),(/"lat_1","lon_1_w"/))
    filevardef(ncdf,"NPP_1_e_hi_int",typeof(NPP_1_e_hi_int),(/"lat_1","lon_1_e"/))
    filevardef(ncdf,"NPP_2_w_hi_int",typeof(NPP_2_w_hi_int),(/"lat_2","lon_2_w"/))
    filevardef(ncdf,"NPP_2_e_hi_int",typeof(NPP_2_e_hi_int),(/"lat_2","lon_2_e"/))
    filevardef(ncdf,"NPP_3_w_hi_int",typeof(NPP_3_w_hi_int),(/"lat_3","lon_3_w"/))
    filevardef(ncdf,"NPP_3_e_hi_int",typeof(NPP_3_e_hi_int),(/"lat_3","lon_3_e"/))


		; now write all the variables to the file
    
		ncdf->lat_1 = lat_1			; Northern Plains
		ncdf->lon_1_w = lon_1_w
		ncdf->lon_1_e = lon_1_e
		ncdf->lat_2 = lat_2			; Central Plains
		ncdf->lon_2_w = lon_2_w
		ncdf->lon_2_e = lon_2_e
		ncdf->lat_3 = lat_3			; Southern Plains
		ncdf->lon_3_w = lon_3_w
		ncdf->lon_3_e = lon_3_e

		ncdf->NPP_1_w_low_freq = (NPP_1_w_low_freq)
		ncdf->NPP_1_e_low_freq = (NPP_1_e_low_freq)
		ncdf->NPP_2_w_low_freq = (NPP_2_w_low_freq)
		ncdf->NPP_2_e_low_freq = (NPP_2_e_low_freq)
		ncdf->NPP_3_w_low_freq = (NPP_3_w_low_freq)
		ncdf->NPP_3_e_low_freq = (NPP_3_e_low_freq)

		ncdf->NPP_1_w_hi_freq = (NPP_1_w_hi_freq)
		ncdf->NPP_1_e_hi_freq = (NPP_1_e_hi_freq)
		ncdf->NPP_2_w_hi_freq = (NPP_2_w_hi_freq)
		ncdf->NPP_2_e_hi_freq = (NPP_2_e_hi_freq)
		ncdf->NPP_3_w_hi_freq = (NPP_3_w_hi_freq)
		ncdf->NPP_3_e_hi_freq = (NPP_3_e_hi_freq)
		
		ncdf->NPP_low_hi_freq_avg = (NPP_low_hi_freq_avg)

		ncdf->NPP_1_w_low_int = (NPP_1_w_low_int)
		ncdf->NPP_1_e_low_int = (NPP_1_e_low_int)
		ncdf->NPP_2_w_low_int = (NPP_2_w_low_int)
		ncdf->NPP_2_e_low_int = (NPP_2_e_low_int)
		ncdf->NPP_3_w_low_int = (NPP_3_w_low_int)
		ncdf->NPP_3_e_low_int = (NPP_3_e_low_int)

		ncdf->NPP_1_w_hi_int = (NPP_1_w_hi_int)
		ncdf->NPP_1_e_hi_int = (NPP_1_e_hi_int)
		ncdf->NPP_2_w_hi_int = (NPP_2_w_hi_int)
		ncdf->NPP_2_e_hi_int = (NPP_2_e_hi_int)
		ncdf->NPP_3_w_hi_int = (NPP_3_w_hi_int)
		ncdf->NPP_3_e_hi_int = (NPP_3_e_hi_int)

		ncdf->NPP_low_hi_int_avg = (NPP_low_hi_int_avg)

	end if

  print("**** DONE ****")
  print("")

	end do	; model loop

end do	; RCP loop

print("******** ALL DONE ********")

end