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

; xy_4.ncl

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

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 "/fis/cgd/cseg/people/jeff/clamp/co2/metrics_table.ncl"

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

begin

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

; read in data: observed

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

 diri  = "/fis/cgd/cseg/people/jeff/clamp_data/co2/"

 fili  = "co2_globalView_98.nc"

 g     = addfile (diri+fili,"r")

 val   = g->CO2_SEAS  

 lon   = g->LON 

 lat   = g->LAT

 sta   = chartostring(g->STATION) 

 delete (g)

;print (sta(0))

 ncase = dimsizes(lat)

;print (ncase)

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

; get only the lowest level at each station 

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

 lat_tmp = lat

 lat_tmp@_FillValue = 1.e+36

 do n = 0,ncase-1

    if (.not. ismissing(lat_tmp(n))) then 

       indexes = ind(lat(n) .eq. lat .and. lon(n) .eq. lon)

       if (dimsizes(indexes) .gt. 1) then

          lat_tmp(indexes(1:)) = lat_tmp@_FillValue

       end if

       delete (indexes)

    end if

 end do

 indexes = ind(.not. ismissing(lat_tmp))

;print (dimsizes(indexes))

;print (indexes)

 lat_ob = lat(indexes)

 lon_ob = lon(indexes)

 val_ob = val(indexes,:)

;printVarSummary (val_ob)

;print (lat_ob +"/"+lon_ob)

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

; read in model data

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

  diri2 = "/fis/cgd/cseg/people/jeff/clamp_data/model/"

; fili2 = "b30.061m_401_425_MONS_climo_atm.nc"

  fili2 = "b30.061n_1995-2004_MONS_climo_atm.nc"

  g     = addfile(diri2+fili2,"r")

  x     = g->CO2

  xi    = g->lon

  yi    = g->lat

  xdim  = dimsizes(x)

  nlev  = xdim(1)

  y     = x(:,0,:,:)

; printVarSummary (y)

; get the co2 at the lowest level

  y     = x(:,nlev-1,:,:)

; change to unit of observed (u mol/mol)

; Model_units [=] kgCO2 / kgDryAir

; 28.966 = molecular weight of dry air

; 44.       = molecular weight of CO2

; u mol = 1e-6 mol

  factor = (28.966/44.) * 1e6

  y      = y * factor

  y@_FillValue = 1.e36

  y@units      = "u mol/mol"

; y = where(y0 .lt. 287.,y@_FillValue,y)

; printVarSummary (y)

; print (min(y)+"/"+max(y))

; interpolate into observed station

; note: model is 0-360E, 90S-90N

; to be able to handle observation at (-89.98,-24.80)

  print (yi(0))

  yi(0) = -90.

  i = ind(lon_ob .lt. 0.)

  lon_ob(i) = lon_ob(i) + 360.  

  yo = linint2_points_Wrap(xi,yi,y,True,lon_ob,lat_ob,0)

  val_model = yo(pts|:,time|:)

; printVarSummary (val_model)

; print (min(val_model)+"/"+max(val_model))

; remove annual mean

  val_model = val_model - conform(val_model,dim_avg(val_model),0)

; print (min(val_model)+"/"+max(val_model))

  nzone = 4

; column

  case_zone = (/"Stations","Amplitude Ratio", \

                "Correlation Coef","M score","Combined Score"/)

  nCase_zone = dimsizes(case_zone ) 

; row

  var_zone = (/"60N-90N","30N-60N","EQ-30N","90S-EQ","Total"/)  

  nVar_zone = dimsizes(var_zone)                  

; arrays to be passed to diagram. 

  case_value_zone = new ((/nCase_zone, nVar_zone/),float )  

  do z = 0,nzone-1

  if (z .eq. 0) then 

;    maximum score for the zone, 60N-90N

     zone = "60N-90N" 

     score_max = 5.0

;    index of stations in this zone

     ind_z = ind(lat_ob .ge. 60.)

;    print (ind_z)

;    print (lat_ob(ind_z)+"/"+lon_ob(ind_z))

;    print (val_ob(ind_z,:))

;    print (val_model(ind_z,:))

  end if

  if (z .eq. 1) then 

;    maximum score for the zone, 30N-60N

     zone = "30N-60N" 

     score_max = 5.0

;    index of stations in this zone

     ind_z = ind(lat_ob .ge. 30. .and. lat_ob .lt. 60.)

;    print (ind_z)

;    print (lat_ob(ind_z)+"/"+lon_ob(ind_z))

;    print (val_ob(ind_z,:))

;    print (val_model(ind_z,:))

  end if

  if (z .eq. 2) then 

;    maximum score for the zone, EQ-30N 

     zone = "EQ-30N"

     score_max = 5.0

;    index of stations in this zone

     ind_z = ind(lat_ob .ge. 0. .and. lat_ob .lt. 30.)

;    print (ind_z)

;    print (lat_ob(ind_z)+"/"+lon_ob(ind_z))

;    print (val_ob(ind_z,:))

;    print (val_model(ind_z,:))

  end if

  if (z .eq. 3) then 

;    maximum score for the zone, 90S-EQ

     zone = "90S-EQ" 

     score_max = 5.0

;    index of stations in this zone

     ind_z = ind(lat_ob .lt. 0. )

;    print (ind_z)

;    print (lat_ob(ind_z)+"/"+lon_ob(ind_z))

;    print (val_ob(ind_z,:))

;    print (val_model(ind_z,:))

  end if

 npts = dimsizes(ind_z)

 print (npts)

 amp_ob        = new((/npts/),float)

 amp_model     = new((/npts/),float)

 amp_ratio_sta = new((/npts/),float)

 ccr_sta       = new((/npts/),float)

 M_sta         = new((/npts/),float)

 score_sta     = new((/npts/),float)

 do n=0,npts-1

    amp_ob(n)    = max(val_ob(ind_z(n),:)) - min(val_ob(ind_z(n),:)) 

    amp_model(n) = max(val_model(ind_z(n),:)) - min(val_model(ind_z(n),:))

    amp_ratio_sta(n) = amp_model(n)/amp_ob(n)

    ccr_sta(n) = esccr(val_ob(ind_z(n),:),val_model(ind_z(n),:),0)

    M_sta(n) = 1.-abs(amp_ratio_sta(n)-1.)

    score_sta(n) = (ccr_sta(n)*ccr_sta(n) + M_sta(n))*0.5 * score_max

    print (sta(ind_z(n))+"/"+lat(ind_z(n))+"/"+lon(ind_z(n))+"/"+amp_ratio_sta(n)+"/"+ccr_sta(n)+"/"+M_sta(n)+"/"+score_sta(n))  

 end do

 amp_ratio_zone = avg(amp_ratio_sta)

 ccr_zone       = avg(ccr_sta)

 M_zone   = 1.- (sum(abs(amp_model-amp_ob)/(amp_model+amp_ob))/npts) 

 score_zone   = (ccr_zone*ccr_zone + M_zone)*0.5 * score_max

 print (npts+"/"+amp_ratio_zone+"/"+ccr_zone+"/"+M_zone+"/"+score_zone)

  case_value_zone(0,z) = npts

  case_value_zone(1,z) = (/amp_ratio_zone/)

  case_value_zone(2,z) = (/ccr_zone/)

  case_value_zone(3,z) = (/M_zone/)

  case_value_zone(4,z) = (/score_zone/)  

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

; plot station table

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

; column for station table

  case_sta   = (/"Latitude","Longitude","Amplitude Ratio", \

                 "Correlation Coef","M score","Combined Score"/) 

  nCase_sta  = dimsizes(case_sta )                 

; row for station table

  var_sta  = sta(ind_z) 

  nVar_sta = dimsizes(var_sta)                 

; arrays to be passed to diagram. 

  case_value_sta = new ((/nCase_sta, nVar_sta/),float )  

  case_value_sta(0,:) = (/lat(ind_z)/)

  case_value_sta(1,:) = (/lon(ind_z)/)

  case_value_sta(2,:) = (/amp_ratio_sta/)

  case_value_sta(3,:) = (/ccr_sta/)

  case_value_sta(4,:) = (/M_sta/)

  case_value_sta(5,:) = (/score_sta/)

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

; plot station table

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

  tt_opt        = True

  tt_opt@pltType= "png" ; "eps" [default], "pdf", "ps"

                        ; "png", "gif" [if you have ImageMajik 'convert']

  tt_opt@tableTitle = "Zone "+ zone

  plot_name         = "table_sta." + zone

  metrics_table(plot_name, var_sta, case_sta , case_value_sta, tt_opt)

 delete (ind_z)

 delete (amp_model)

 delete (amp_ob)

 delete (amp_ratio_sta)

 delete (ccr_sta)

 delete (M_sta)

 delete (score_sta)

 delete (var_sta)

 delete (case_value_sta)

 end do

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

; plot zone table

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

  case_value_zone(0,4) = case_value_zone(0,0)+case_value_zone(0,1)+case_value_zone(0,2)+case_value_zone(0,3) 

  case_value_zone(1,4) = 0.

  case_value_zone(2,4) = 0.

  case_value_zone(3,4) = 0.

  case_value_zone(4,4) = case_value_zone(4,0)+case_value_zone(4,1)+case_value_zone(4,2)+case_value_zone(4,3) 

  tt_opt        = True

  tt_opt@pltType= "png" ; "eps" [default], "pdf", "ps"

                        ; "png", "gif" [if you have ImageMajik 'convert']

  tt_opt@tableTitle = "Zone"

  plot_name         = "table_zone" 

  metrics_table(plot_name, var_zone, case_zone , case_value_zone, tt_opt)

end