forrest@0: ; ***********************************************
forrest@0: ; xy_4.ncl
forrest@0: ; ***********************************************
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
forrest@0: load "/fis/cgd/cseg/people/jeff/clamp/co2/metrics_table.ncl"
forrest@0: ;************************************************
forrest@0: begin
forrest@0: ;************************************************
forrest@0: ; read in data: observed
forrest@0: ;************************************************
forrest@0:  diri  = "/fis/cgd/cseg/people/jeff/clamp_data/co2/"
forrest@0:  fili  = "co2_globalView_98.nc"
forrest@0:  g     = addfile (diri+fili,"r")
forrest@0:  val   = g->CO2_SEAS  
forrest@0:  lon   = g->LON 
forrest@0:  lat   = g->LAT
forrest@0:  sta   = chartostring(g->STATION) 
forrest@0:  delete (g)
forrest@0:  
forrest@0: ;print (sta(0))
forrest@0: 
forrest@0:  ncase = dimsizes(lat)
forrest@0: ;print (ncase)
forrest@0: 
forrest@0: ;**************************************************************
forrest@0: ; get only the lowest level at each station 
forrest@0: ;**************************************************************
forrest@0:  lat_tmp = lat
forrest@0:  lat_tmp@_FillValue = 1.e+36
forrest@0:  
forrest@0:  do n = 0,ncase-1
forrest@0:     if (.not. ismissing(lat_tmp(n))) then 
forrest@0:        indexes = ind(lat(n) .eq. lat .and. lon(n) .eq. lon)
forrest@0:        if (dimsizes(indexes) .gt. 1) then
forrest@0:           lat_tmp(indexes(1:)) = lat_tmp@_FillValue
forrest@0:        end if
forrest@0:        delete (indexes)
forrest@0:     end if
forrest@0:  end do
forrest@0: 
forrest@0:  indexes = ind(.not. ismissing(lat_tmp))
forrest@0: ;print (dimsizes(indexes))
forrest@0: ;print (indexes)
forrest@0:  
forrest@0:  lat_ob = lat(indexes)
forrest@0:  lon_ob = lon(indexes)
forrest@0:  val_ob = val(indexes,:)
forrest@0: ;printVarSummary (val_ob)
forrest@0: ;print (lat_ob +"/"+lon_ob)
forrest@0: 
forrest@0: ;************************************************
forrest@0: ; read in model data
forrest@0: ;************************************************
forrest@0:   diri2 = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
forrest@0: ; fili2 = "b30.061m_401_425_MONS_climo_atm.nc"
forrest@0:   fili2 = "b30.061n_1995-2004_MONS_climo_atm.nc"
forrest@0: 
forrest@0:   g     = addfile(diri2+fili2,"r")
forrest@0:   x     = g->CO2
forrest@0:   xi    = g->lon
forrest@0:   yi    = g->lat
forrest@0:   xdim  = dimsizes(x)
forrest@0:   nlev  = xdim(1)
forrest@0:   y     = x(:,0,:,:)
forrest@0: ; printVarSummary (y)
forrest@0:   
forrest@0: ; get the co2 at the lowest level
forrest@0:   y     = x(:,nlev-1,:,:)
forrest@0: 
forrest@0: ; change to unit of observed (u mol/mol)
forrest@0: ; Model_units [=] kgCO2 / kgDryAir
forrest@0: ; 28.966 = molecular weight of dry air
forrest@0: ; 44.       = molecular weight of CO2
forrest@0: ; u mol = 1e-6 mol
forrest@0: 
forrest@0:   factor = (28.966/44.) * 1e6
forrest@0:   y      = y * factor
forrest@0: 
forrest@0:   y@_FillValue = 1.e36
forrest@0:   y@units      = "u mol/mol"
forrest@0: ; y = where(y0 .lt. 287.,y@_FillValue,y)
forrest@0: ; printVarSummary (y)
forrest@0: ; print (min(y)+"/"+max(y))
forrest@0: 
forrest@0: ; interpolate into observed station
forrest@0: ; note: model is 0-360E, 90S-90N
forrest@0: 
forrest@0: ; to be able to handle observation at (-89.98,-24.80)
forrest@0:   print (yi(0))
forrest@0:   yi(0) = -90.
forrest@0: 
forrest@0:   i = ind(lon_ob .lt. 0.)
forrest@0:   lon_ob(i) = lon_ob(i) + 360.  
forrest@0: 
forrest@0:   yo = linint2_points_Wrap(xi,yi,y,True,lon_ob,lat_ob,0)
forrest@0: 
forrest@0:   val_model = yo(pts|:,time|:)
forrest@0: ; printVarSummary (val_model)
forrest@0: ; print (min(val_model)+"/"+max(val_model))
forrest@0: 
forrest@0: ; remove annual mean
forrest@0:   val_model = val_model - conform(val_model,dim_avg(val_model),0)
forrest@0: ; print (min(val_model)+"/"+max(val_model))
forrest@0: 
forrest@0:   nzone = 1
forrest@0:   do z = 0,nzone-1
forrest@0: 
forrest@0:   if (z .eq. 0) then 
forrest@0: ;    maximum score for the zone, 60N-90N 
forrest@0:      score_max = 5.0
forrest@0: ;    index of stations in this zone
forrest@0:      ind_z = ind(lat_ob .ge. 60.)
forrest@0: ;    print (ind_z)
forrest@0: ;    print (lat_ob(ind_z)+"/"+lon_ob(ind_z))
forrest@0: ;    print (val_ob(ind_z,:))
forrest@0: ;    print (val_model(ind_z,:))
forrest@0:   end if
forrest@0: 
forrest@0:   if (z .eq. 1) then 
forrest@0: ;    maximum score for the zone, 30N-60N 
forrest@0:      score_max = 5.0
forrest@0: ;    index of stations in this zone
forrest@0:      ind_z = ind(lat_ob .ge. 30. .and. lat_ob .lt. 60.)
forrest@0: ;    print (ind_z)
forrest@0: ;    print (lat_ob(ind_z)+"/"+lon_ob(ind_z))
forrest@0: ;    print (val_ob(ind_z,:))
forrest@0: ;    print (val_model(ind_z,:))
forrest@0:   end if
forrest@0: 
forrest@0:   if (z .eq. 2) then 
forrest@0: ;    maximum score for the zone, EQ-30N 
forrest@0:      score_max = 5.0
forrest@0: ;    index of stations in this zone
forrest@0:      ind_z = ind(lat_ob .ge. 0. .and. lat_ob .lt. 30.)
forrest@0: ;    print (ind_z)
forrest@0: ;    print (lat_ob(ind_z)+"/"+lon_ob(ind_z))
forrest@0: ;    print (val_ob(ind_z,:))
forrest@0: ;    print (val_model(ind_z,:))
forrest@0:   end if
forrest@0: 
forrest@0:   if (z .eq. 3) then 
forrest@0: ;    maximum score for the zone, 90S-EQ 
forrest@0:      score_max = 5.0
forrest@0: ;    index of stations in this zone
forrest@0:      ind_z = ind(lat_ob .lt. 0. )
forrest@0: ;    print (ind_z)
forrest@0: ;    print (lat_ob(ind_z)+"/"+lon_ob(ind_z))
forrest@0: ;    print (val_ob(ind_z,:))
forrest@0: ;    print (val_model(ind_z,:))
forrest@0:   end if
forrest@0: 
forrest@0:  npts = dimsizes(ind_z)
forrest@0:  print (npts)
forrest@0: 
forrest@0:  amp_ob        = new((/npts/),float)
forrest@0:  amp_model     = new((/npts/),float)
forrest@0: 
forrest@0:  amp_ratio_sta = new((/npts/),float)
forrest@0:  ccr_sta       = new((/npts/),float)
forrest@0:  M_sta         = new((/npts/),float)
forrest@0:  score_sta     = new((/npts/),float)
forrest@0: 
forrest@0:  do n=0,npts-1
forrest@0:     amp_ob(n)    = max(val_ob(ind_z(n),:)) - min(val_ob(ind_z(n),:)) 
forrest@0:     amp_model(n) = max(val_model(ind_z(n),:)) - min(val_model(ind_z(n),:))
forrest@0: 
forrest@0:     amp_ratio_sta(n) = amp_model(n)/amp_ob(n)
forrest@0:     ccr_sta(n) = esccr(val_ob(ind_z(n),:),val_model(ind_z(n),:),0)
forrest@0:     M_sta(n) = 1.-abs(amp_ratio_sta(n)-1.)
forrest@0:     score_sta(n) = (ccr_sta(n)*ccr_sta(n) + M_sta(n))*0.5 * score_max
forrest@0: 
forrest@0:     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))  
forrest@0:  end do
forrest@0: 
forrest@0:  amp_ratio_zone = avg(amp_ratio_sta)
forrest@0:  ccr_zone       = avg(ccr_sta)
forrest@0:  M_zone   = 1.- (sum(abs(amp_model-amp_ob)/(amp_model+amp_ob))/npts) 
forrest@0:  score_zone   = (ccr_zone*ccr_zone + M_zone)*0.5 * score_max
forrest@0: 
forrest@0:  print (npts+"/"+amp_ratio_zone+"/"+ccr_zone+"/"+M_zone+"/"+score_zone)  
forrest@0: ;****************************************************************************
forrest@0: ; Cases [Model]
forrest@0:   case      = (/ "Lat", "Lon", "AR" /) 
forrest@0:   nCase     = dimsizes(case )                 ; # of Cases [Cases]
forrest@0: 
forrest@0: ; variables compared
forrest@0:   var       = sta(ind_z) 
forrest@0:   nVar      = dimsizes(var)                   ; # of Variables
forrest@0: 
forrest@0: ; "Case A"                        
forrest@0:   CA_ratio   = (/lat(ind_z)/)
forrest@0: 
forrest@0: ; "Case B"                        
forrest@0:   CB_ratio   = (/lon(ind_z)/)
forrest@0: 
forrest@0: ; "Case C"                        
forrest@0:   CC_ratio   = (/amp_ratio_sta/)
forrest@0: 
forrest@0: 
forrest@0: ; arrays to be passed to taylor_diagram. It will calculate the x xnd y coordinates.
forrest@0:   ratio      = new ((/nCase, nVar/),typeof(CA_ratio) )  
forrest@0: 
forrest@0:   ratio(0,:) = CA_ratio 
forrest@0:   ratio(1,:) = CB_ratio
forrest@0:   ratio(2,:) = CC_ratio
forrest@0: 
forrest@0: ;**************************************************
forrest@0: ; fill an array for a "taylor metrics table"
forrest@0: ;**************************************************
forrest@0: 
forrest@0: ; season    = (/ "ANN" /)
forrest@0: ; nSeason   = dimsizes(season)
forrest@0:   season    = (/ "" /)
forrest@0:   nSeason   = dimsizes(season)
forrest@0: 
forrest@0:   table     = new ( (/nCase,nVar/), typeof(ratio) )
forrest@0:   table(0,:) = CA_ratio
forrest@0:   table(1,:) = CB_ratio
forrest@0:   table(2,:) = CC_ratio
forrest@0: 
forrest@0:   tt_opt        = True
forrest@0:   tt_opt@pltType= "ps"                  ; "eps" [default], "pdf", "ps"
forrest@0:                                          ; "png", "gif" [if you have ImageMajik 'convert']
forrest@0: 
forrest@0:   tt_opt@tableTitle = "Station in 60N_90N"
forrest@0: 
forrest@0:   varSource =var
forrest@0: 
forrest@0:   metrics_table("metrics", varSource, case , table, tt_opt)
forrest@0: 
forrest@0:  delete (ind_z)
forrest@0:  delete (amp_model)
forrest@0:  delete (amp_ob)
forrest@0:  delete (amp_ratio_sta)
forrest@0:  delete (ccr_sta)
forrest@0:  delete (M_sta)
forrest@0:  delete (score_sta)
forrest@0:  end do
forrest@0: end