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

; add another model to plot

; add panel plot to 22.lines.ncl

; add zone plot to 21.lines.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"

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

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.061n_1995-2004_MONS_climo_atm.nc"

  fili3 = "b30.061m_401_425_MONS_climo_atm.nc"

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

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

  x1   = g->CO2

  xi   = g->lon

  yi   = g->lat

  delete (g)

  xdim  = dimsizes(x1)

  nlev  = xdim(1)

  y1     = x1(:,0,:,:)

; printVarSummary (y1)

; get the co2 at the lowest level

  y1     = x1(:,nlev-1,:,:)

  delete (x1)

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

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

  x2    = g->CO2

  delete (g)

  y2     = x2(:,0,:,:)

  y2     = x2(:,nlev-1,:,:)

  delete (x2)

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

; 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

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

  y1 = y1 * factor

  y1@_FillValue = 1.e36

  y1@units      = "u mol/mol"

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

; printVarSummary (y1)

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

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

  y2 = y2 * factor

  y2@_FillValue = 1.e36

  y2@units      = "u mol/mol"

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

; interpolate into observed station

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

;       ob    is -180-180, 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,y1,True,lon_ob,lat_ob,0)

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

  val_model1_0 = val_model1

; printVarSummary (val_model1)

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

; remove annual mean

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

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

  delete (yo)

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

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

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

  val_model2_0 = val_model2

; printVarSummary (val_model2)

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

; remove annual mean

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

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

  delete (yo)

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

; for x-axis in xyplot

  mon = ispan(1,12,1)

  mon@long_name = "month"

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

  plot_type = "ps"

  plot_type_new = "png"

  res                   = True                      ; plot mods desired

  res@xyLineThicknesses = (/2.0,2.0,2.0/)           ; make 2nd lines thicker

  res@xyLineColors      = (/"red","blue","black"/)  ; change line color

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

; Add a boxed legend using the more simple method

  res@pmLegendDisplayMode    = "Always"

; res@pmLegendWidthF         = 0.1

  res@pmLegendWidthF         = 0.08

  res@pmLegendHeightF        = 0.06

; res@pmLegendOrthogonalPosF = -1.17

; res@pmLegendOrthogonalPosF = -1.00  ;(downward)

  res@pmLegendOrthogonalPosF = -0.30  ;(downward)

; res@pmLegendParallelPosF   =  0.18

  res@pmLegendParallelPosF   =  0.23  ;(rightward)

; res@lgPerimOn             = False

  res@lgLabelFontHeightF     = 0.015

  res@xyExplicitLegendLabels = (/"b30.061n","b30.061m","observed"/)

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

  nzone = 4

  do z = 0,nzone-1

  if (z .eq. 0) then 

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

     score_max = 5.0

     zone = "60N-90N"

;    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 

     score_max = 5.0

     zone = "30N-60N"

;    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 

     score_max = 5.0

     zone = "EQ-30N"

;    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 

     score_max = 5.0

     zone = "90S-EQ"     

;    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)

  npts_str = ""

  npts_str = npts

; print (npts_str)

  plot_data   = new((/3,12,npts/),float)

  plot_data_0 = new((/2,12,npts/),float)

  plot_data!0 = "case"

  plot_data!1 = "month"

  plot_data!2 = "pts"

  plot_data@long_name   = "CO2 Seasonal"

  plot_data_0!0 = "month"

  plot_data_0!1 = "pts"

  plot_data_0@long_name = "CO2"

  do n=0,npts-1

    plot_data(0,:,n) = (/val_model1(ind_z(n),:)/)

    plot_data(1,:,n) = (/val_model2(ind_z(n),:)/)

    plot_data(2,:,n) = (/val_ob(ind_z(n),:)/)

    plot_data_0(0,:,n) = (/val_model1_0(ind_z(n),:)/)

    plot_data_0(1,:,n) = (/val_model2_0(ind_z(n),:)/)

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

; create plot with 3 lines for each observed station

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

    plot_name = sta(ind_z(n))    

    title = plot_name+"("+lat(ind_z(n))+","+lon(ind_z(n))+")"

;   print (title)

;   print (plot_name)

    wks = gsn_open_wks (plot_type,plot_name)   ; open workstation

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

;   for panel plot

    plot=new(2,graphic)                        ; create graphic array

    res@gsnFrame     = False                   ; Do not draw plot 

    res@gsnDraw      = False                   ; Do not advance frame

    pres                            = True     ; panel plot mods desired

    pres@gsnPanelYWhiteSpacePercent = 5        ; increase white space around

                                               ; indiv. plots in panel

    pres@gsnMaximize                = True     ; fill the page

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

    res@tiMainString = title                           ; add title

    plot(0)=gsn_csm_xy(wks,mon,plot_data(:,:,n),res)   ; create plot 1

;   plot(1)=gsn_csm_xy(wks,mon,plot_data_0(:,:,n),res) ; create plot 2

    plot(1)=gsn_csm_xy(wks,mon,plot_data_0(0,:,n),res) ; create plot 2

    gsn_panel(wks,plot,(/2,1/),pres)                 ; create panel plot

    system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

    system("rm "+plot_name+"."+plot_type)

    clear (wks)

 end do

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

; plot zone average

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

    plot_name = "All_"+npts_str

    title = plot_name + " in "+ zone

;   print (title)

;   print (plot_name)

    wks = gsn_open_wks (plot_type,plot_name)        ; open workstation

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

;   for panel plot

    plot=new(2,graphic)                        ; create graphic array

    res@gsnFrame     = False                   ; Do not draw plot 

    res@gsnDraw      = False                   ; Do not advance frame

    pres                            = True     ; panel plot mods desired

    pres@gsnPanelYWhiteSpacePercent = 5        ; increase white space around

                                               ; indiv. plots in panel

    pres@gsnMaximize                = True     ; fill the page

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

    res@tiMainString = title                                     ; add title

    plot(0) = gsn_csm_xy (wks,mon,dim_avg_Wrap(plot_data),res)   ; create plot 1

;   plot(1) = gsn_csm_xy (wks,mon,dim_avg_Wrap(plot_data_0),res) ; create plot 2

    xx = dim_avg_Wrap(plot_data_0)

    plot(1) = gsn_csm_xy (wks,mon,xx(0,:),res) ; create plot 2

    gsn_panel(wks,plot,(/2,1/),pres)                 ; create panel plot

    system("convert "+plot_name+"."+plot_type+" "+plot_name+"."+plot_type_new)

    system("rm "+plot_name+"."+plot_type)

    clear (wks)

    delete (ind_z)

    delete (plot_data)

    delete (plot_data_0)    

 end do

end