lai/99.all.ncl.new
changeset 1 4be95183fbcd
equal deleted inserted replaced
-1:000000000000 0:804435136736
       
     1 ;*************************************************************
       
     2 ; remove histogram plots.
       
     3 ;
       
     4 ; required command line input parameters:
       
     5 ;  ncl 'model_name="10cn" model_grid="T42" dirm="/.../ film="..."' 01.npp.ncl
       
     6 ;
       
     7 ; histogram normalized by rain and compute correleration
       
     8 ;**************************************************************
       
     9 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl.test"
       
    10 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl.test"
       
    11 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
       
    12 ;**************************************************************
       
    13 procedure set_line(lines:string,nline:integer,newlines:string) 
       
    14 begin
       
    15 ; add line to ascci/html file
       
    16     
       
    17   nnewlines = dimsizes(newlines)
       
    18   if(nline+nnewlines-1.ge.dimsizes(lines))
       
    19     print("set_line: bad index, not setting anything.") 
       
    20     return
       
    21   end if 
       
    22   lines(nline:nline+nnewlines-1) = newlines
       
    23 ;  print ("lines = " + lines(nline:nline+nnewlines-1))
       
    24   nline = nline + nnewlines
       
    25   return 
       
    26 end
       
    27 ;**************************************************************
       
    28 ; Main code.
       
    29 begin
       
    30  
       
    31   plot_type     = "ps"
       
    32   plot_type_new = "png"
       
    33 
       
    34 ;************************************************
       
    35 ; read data: model       
       
    36 ;************************************************
       
    37 
       
    38 ;film  = "b30.061n_1995-2004_MONS_climo_lnd.nc"
       
    39 ;model_name = "b30.061n"
       
    40 ;model_grid = "T31"
       
    41 
       
    42 ;film  = "newcn05_ncep_1i_MONS_climo_lnd.nc"
       
    43 ;model_name = "newcn"
       
    44 ;model_grid = "1.9"
       
    45 
       
    46 ;film  = "i01.06cn_1798-2004_MONS_climo.nc"
       
    47 ;model_name = "06cn"
       
    48 ;model_grid = "T42"
       
    49 
       
    50 ;film  = "i01.06casa_1798-2004_MONS_climo.nc"
       
    51 ;model_name = "06casa"
       
    52 ;model_grid = "T42"
       
    53 
       
    54  film  = "i01.10cn_1948-2004_MONS_climo.nc"
       
    55  model_name = "10cn"
       
    56  model_grid = "T42"
       
    57 
       
    58 ;film  = "i01.10casa_1948-2004_MONS_climo.nc"
       
    59 ;model_name = "10casa"
       
    60 ;model_grid = "T42"
       
    61 
       
    62  dirm  = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
       
    63  fm    = addfile(dirm+film,"r")
       
    64       
       
    65  laimod  = fm->TLAI
       
    66       
       
    67 ;************************************************
       
    68 ; read data: observed
       
    69 ;************************************************
       
    70 
       
    71   ob_name = "MODIS MOD 15A2 2000-2005"
       
    72 
       
    73   diro  = "/fis/cgd/cseg/people/jeff/clamp_data/lai/ob/"
       
    74   filo1  = "land_class_"+model_grid+".nc"
       
    75   filo2  = "LAI_2000-2005_MONS_"+model_grid+".nc"
       
    76 
       
    77   fo1 = addfile(diro+filo1,"r")
       
    78   fo2 = addfile(diro+filo2,"r")
       
    79  
       
    80   classob    = tofloat(fo1->LAND_CLASS)               
       
    81   laiob      = fo2->LAI
       
    82 
       
    83 ; input observed data has 20 land-type classes
       
    84   nclass = 20
       
    85 
       
    86 ;************************************************
       
    87 ; global res
       
    88 ;************************************************
       
    89   resg                      = True             ; Use plot options
       
    90   resg@cnFillOn             = True             ; Turn on color fill
       
    91   resg@gsnSpreadColors      = True             ; use full colormap
       
    92   resg@cnLinesOn            = False            ; Turn off contourn lines
       
    93   resg@mpFillOn             = False            ; Turn off map fill
       
    94   resg@cnLevelSelectionMode = "ManualLevels"   ; Manual contour invtervals
       
    95 
       
    96 ;************************************************
       
    97 ; plot global land class: observed
       
    98 ;************************************************
       
    99 
       
   100   resg@cnMinLevelValF       = 1.              ; Min level
       
   101   resg@cnMaxLevelValF       = 19.             ; Max level
       
   102   resg@cnLevelSpacingF      = 1.              ; interval
       
   103 
       
   104 ; global contour ob
       
   105   classob@_FillValue = 1.e+36
       
   106   classob = where(classob.eq.0,classob@_FillValue,classob)
       
   107   
       
   108   plot_name = "global_class_ob"
       
   109   title     = ob_name
       
   110   resg@tiMainString  = title
       
   111 
       
   112   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation
       
   113   gsn_define_colormap(wks,"gui_default")     ; choose colormap
       
   114 
       
   115   plot = gsn_csm_contour_map_ce(wks,classob,resg)   
       
   116   frame(wks)
       
   117 
       
   118   clear (wks)
       
   119  
       
   120 ;*******************************************************************
       
   121 ; for html table : all 4 components (Mean, Max, Phase, Growth)
       
   122 ;*******************************************************************
       
   123 
       
   124 ; column (not including header column)
       
   125 
       
   126   component = (/"Mean","Max","Phase","Growth"/)
       
   127   col_head2 = (/"observed",model_name,"M_score" \
       
   128                ,"observed",model_name,"M_score" \
       
   129                ,"observed",model_name,"M_score" \
       
   130                ,"observed",model_name,"M_score" \
       
   131                /)
       
   132   
       
   133   n_comp = dimsizes(component) 
       
   134   ncol   = dimsizes(col_head2)
       
   135 
       
   136 ; row (not including header row)
       
   137   row_head  = (/"Evergreen Needleleaf Forests" \
       
   138                ,"Evergreen Broadleaf Forests" \
       
   139                ,"Deciduous Needleleaf Forest" \
       
   140                ,"Deciduous Broadleaf Forests" \
       
   141                ,"Mixed Forests" \                      
       
   142                ,"Closed Bushlands" \                   
       
   143                ,"Open Bushlands" \                     
       
   144                ,"Woody Savannas (S. Hem.)" \           
       
   145                ,"Savannas (S. Hem.)" \                 
       
   146                ,"Grasslands" \                         
       
   147                ,"Permanent Wetlands" \                 
       
   148                ,"Croplands" \                                           
       
   149                ,"Cropland/Natural Vegetation Mosaic" \             
       
   150                ,"Barren or Sparsely Vegetated" \                             
       
   151                ,"Woody Savannas (N. Hem.)" \           
       
   152                ,"Savannas (N. Hem.)" \
       
   153                ,"All Biome" \                
       
   154                /)  
       
   155   nrow = dimsizes(row_head)                  
       
   156 
       
   157 ; arrays to be passed to table. 
       
   158   text4 = new ((/nrow, ncol/),string )
       
   159 
       
   160 ; M_comp
       
   161   M_comp = (/"M_lai_mean","M_lai_max","M_lai_phase","M_lai_grow"/) 
       
   162 
       
   163 ; total M_score
       
   164   M_total = 0.
       
   165 
       
   166 ;********************************************************************
       
   167 ; use land-type class to bin the data in equally spaced ranges
       
   168 ;********************************************************************
       
   169 
       
   170   nclassn     = nclass + 1
       
   171   range       = fspan(0,nclassn-1,nclassn)
       
   172 ; print (range)
       
   173 
       
   174 ; Use this range information to grab all the values in a
       
   175 ; particular range, and then take an average.
       
   176 
       
   177   nr           = dimsizes(range)
       
   178   nx           = nr-1
       
   179   xvalues      = new((/2,nx/),float)
       
   180   xvalues(0,:) = range(0:nr-2) + (range(1:)-range(0:nr-2))/2.
       
   181   dx           = xvalues(0,1) - xvalues(0,0)       ; range width
       
   182   dx4          = dx/4                              ; 1/4 of the range
       
   183   xvalues(1,:) = xvalues(0,:) - dx/5.
       
   184 
       
   185 ;************************************************************************
       
   186 ; go through all components
       
   187 ;************************************************************************
       
   188 
       
   189   do n = 0,n_comp-1
       
   190 
       
   191 ;===================
       
   192 ; get data:
       
   193 ;===================
       
   194 ; (A) Mean
       
   195  
       
   196   if (n .eq. 0) then
       
   197      data_ob  = dim_avg_Wrap(laiob (lat|:,lon|:,time|:))
       
   198      data_mod = dim_avg_Wrap(laimod(lat|:,lon|:,time|:))
       
   199   end if
       
   200 
       
   201 ; (B) Max
       
   202 
       
   203   if (n .eq. 1) then
       
   204 
       
   205 ;    observed  
       
   206      data_ob = laiob(0,:,:)
       
   207      s       = laiob(:,0,0)
       
   208      data_ob@long_name = "Leaf Area Index Max"
       
   209  
       
   210      dsizes_z = dimsizes(laiob)
       
   211      nlat     = dsizes_z(1)
       
   212      nlon     = dsizes_z(2)
       
   213   
       
   214      do j = 0,nlat-1
       
   215      do i = 0,nlon-1
       
   216         s = laiob(:,j,i) 
       
   217         data_ob(j,i) = max(s)
       
   218      end do
       
   219      end do
       
   220 
       
   221      delete (s)
       
   222      delete (dsizes_z)          
       
   223 
       
   224 ;    model  
       
   225      data_mod = laimod(0,:,:)
       
   226      s        = laimod(:,0,0)
       
   227      data_mod@long_name = "Leaf Area Index Max"
       
   228  
       
   229      dsizes_z = dimsizes(laimod)
       
   230      nlat     = dsizes_z(1)
       
   231      nlon     = dsizes_z(2)
       
   232   
       
   233      do j = 0,nlat-1
       
   234      do i = 0,nlon-1
       
   235         s = laimod(:,j,i) 
       
   236         data_mod(j,i) = max(s)
       
   237      end do
       
   238      end do
       
   239 
       
   240      delete (s)
       
   241      delete (dsizes_z)          
       
   242   end if
       
   243 
       
   244 ; (C) phase
       
   245 
       
   246   if (n .eq. 2) then  
       
   247 
       
   248 ;    observed
       
   249      data_ob = laiob(0,:,:)
       
   250      s       = laiob(:,0,0)
       
   251      data_ob@long_name = "Leaf Area Index Max Month"
       
   252  
       
   253      dsizes_z = dimsizes(laiob)
       
   254      nlat     = dsizes_z(1)
       
   255      nlon     = dsizes_z(2)
       
   256   
       
   257      do j = 0,nlat-1
       
   258      do i = 0,nlon-1
       
   259         s = laiob(:,j,i) 
       
   260         data_ob(j,i) = maxind(s) + 1
       
   261      end do
       
   262      end do
       
   263 
       
   264      delete (s)
       
   265      delete (dsizes_z)          
       
   266 
       
   267 ;    model
       
   268      data_mod = laimod(0,:,:)
       
   269      s        = laimod(:,0,0)
       
   270      data_mod@long_name = "Leaf Area Index Max Month"
       
   271  
       
   272      dsizes_z = dimsizes(laimod)
       
   273      nlat     = dsizes_z(1)
       
   274      nlon     = dsizes_z(2)
       
   275   
       
   276      do j = 0,nlat-1
       
   277      do i = 0,nlon-1
       
   278         s = laimod(:,j,i) 
       
   279         data_mod(j,i) = maxind(s) + 1
       
   280      end do
       
   281      end do
       
   282 
       
   283      delete (s)
       
   284      delete (dsizes_z)          
       
   285   end if
       
   286 
       
   287 ; (D) grow day
       
   288 
       
   289   if (n .eq. 3) then   
       
   290 
       
   291      day_of_data = (/31,28,31,30,31,30,31,31,30,31,30,31/)
       
   292 
       
   293 ;    observed
       
   294      data_ob = laiob(0,:,:)
       
   295      data_ob@long_name = "Days of Growing Season"
       
   296  
       
   297      dsizes_z = dimsizes(laiob)
       
   298      ntime    = dsizes_z(0)
       
   299      nlat     = dsizes_z(1)
       
   300      nlon     = dsizes_z(2)
       
   301   
       
   302      do j = 0,nlat-1
       
   303      do i = 0,nlon-1
       
   304         nday = 0.
       
   305         do k = 0,ntime-1
       
   306            if (.not. ismissing(laiob(k,j,i)) .and. laiob(k,j,i) .gt. 1.0) then
       
   307               nday = nday + day_of_data(k)
       
   308            end if
       
   309         end do
       
   310 
       
   311         data_ob(j,i) = nday
       
   312      end do
       
   313      end do
       
   314 
       
   315      delete (dsizes_z)     
       
   316 
       
   317 ;    model
       
   318      data_mod = laimod(0,:,:)
       
   319      data_mod@long_name = "Days of Growing Season"
       
   320  
       
   321      dsizes_z = dimsizes(laimod)
       
   322      ntime    = dsizes_z(0)
       
   323      nlat     = dsizes_z(1)
       
   324      nlon     = dsizes_z(2)
       
   325   
       
   326      do j = 0,nlat-1
       
   327      do i = 0,nlon-1
       
   328         nday = 0.
       
   329         do k = 0,ntime-1
       
   330            if (.not. ismissing(laimod(k,j,i)) .and. laimod(k,j,i) .gt. 1.0) then
       
   331               nday = nday + day_of_data(k)
       
   332            end if
       
   333         end do
       
   334 
       
   335         data_mod(j,i) = nday
       
   336      end do
       
   337      end do
       
   338 
       
   339      delete (dsizes_z)
       
   340   end if
       
   341 
       
   342 ;==============================
       
   343 ; put data into bins
       
   344 ;==============================
       
   345 
       
   346   base_1D  = ndtooned(classob)
       
   347   data1_1D = ndtooned(data_ob)
       
   348   data2_1D = ndtooned(data_mod)
       
   349 
       
   350 ; output for data in bins
       
   351 
       
   352   yvalues = new((/2,nx/),float)
       
   353   count   = new((/2,nx/),float)
       
   354 
       
   355 ; put data into bins
       
   356 
       
   357   do nd=0,1
       
   358 
       
   359 ;   See if we are doing data1 (nd=0) or data2 (nd=1).
       
   360 
       
   361     base = base_1D
       
   362 
       
   363     if(nd.eq.0) then
       
   364       data = data1_1D
       
   365     else
       
   366       data = data2_1D
       
   367     end if
       
   368 
       
   369 ;   Loop through each range, using base.
       
   370 
       
   371     do i=0,nr-2
       
   372       if (i.ne.(nr-2)) then
       
   373 ;        print("")
       
   374 ;        print("In range ["+range(i)+","+range(i+1)+")")
       
   375          idx = ind((base.ge.range(i)).and.(base.lt.range(i+1)))
       
   376       else
       
   377 ;        print("")
       
   378 ;        print("In range ["+range(i)+",)")
       
   379          idx = ind(base.ge.range(i))
       
   380       end if
       
   381 
       
   382 ;     Calculate average 
       
   383 
       
   384       if(.not.any(ismissing(idx))) then
       
   385         yvalues(nd,i) = avg(data(idx))
       
   386         count(nd,i)   = dimsizes(idx)
       
   387       else
       
   388         yvalues(nd,i) = yvalues@_FillValue
       
   389         count(nd,i)   = 0
       
   390       end if
       
   391 
       
   392 ;#############################################################
       
   393 ;     set the following 4 classes to _FillValue:
       
   394 ;     Water Bodies(0), Urban and Build-Up(13),
       
   395 ;     Permenant Snow and Ice(15), Unclassified(17)
       
   396 
       
   397       if (i.eq.0 .or. i.eq.13 .or. i.eq.15 .or. i.eq.17) then
       
   398          yvalues(nd,i) = yvalues@_FillValue
       
   399          count(nd,i)   = 0
       
   400       end if
       
   401 ;############################################################# 
       
   402 
       
   403 ;     print(nd + ": " + count(nd,i) + " points, avg = " + yvalues(nd,i))
       
   404 
       
   405 ;     Clean up for next time in loop.
       
   406 
       
   407       delete(idx)
       
   408     end do
       
   409 
       
   410     delete(data)
       
   411   end do
       
   412 
       
   413   delete (base)
       
   414   delete (base_1D)
       
   415   delete (data1_1D)
       
   416   delete (data2_1D)
       
   417 
       
   418 ;=====================================
       
   419 ; compute correlation coef and M score 
       
   420 ;=====================================
       
   421 
       
   422   u = yvalues(0,:)
       
   423   v = yvalues(1,:)
       
   424 
       
   425   good = ind(.not.ismissing(u) .and. .not.ismissing(v))
       
   426   uu = u(good)
       
   427   vv = v(good)
       
   428 
       
   429 ; compute correlation coef
       
   430   cc = esccr(uu,vv,0)
       
   431 
       
   432   if (n .eq. 2) then
       
   433      bias   = avg(abs(vv-uu))
       
   434      bias   = where((bias.gt. 6.),12.-bias,bias)
       
   435      Mscore = ((6. - bias)/6.)*5.
       
   436      M_score = sprintf("%.2f", Mscore)
       
   437   else
       
   438      bias  = sum(abs(vv-uu)/abs(vv+uu))
       
   439      Mscore = (1.- (bias/dimsizes(uu)))*5.
       
   440      M_score = sprintf("%.2f", Mscore)
       
   441   end if
       
   442 
       
   443 ; compute M_total
       
   444   
       
   445   M_total = M_total + Mscore
       
   446 
       
   447 ;==================
       
   448 ; output M_score
       
   449 ;==================
       
   450 
       
   451   print (Mscore)
       
   452 ;=======================
       
   453 ; output to html table
       
   454 ;=======================
       
   455 
       
   456   nn = n*3
       
   457 
       
   458   do i=0,nrow-2
       
   459      text4(i,nn)   = sprintf("%.2f",u(i))
       
   460      text4(i,nn+1) = sprintf("%.2f",v(i))
       
   461      text4(i,nn+2) = "-"
       
   462   end do
       
   463   text4(nrow-1,nn)   = sprintf("%.2f",avg(u))
       
   464   text4(nrow-1,nn+1) = sprintf("%.2f",avg(v))
       
   465   text4(nrow-1,nn+2) = M_score
       
   466 
       
   467   delete (u)
       
   468   delete (v)
       
   469   delete (uu)
       
   470   delete (vv)
       
   471   delete (yvalues)
       
   472   delete (good)
       
   473 
       
   474 ;======================================== 
       
   475 ; global res changes for each component
       
   476 ;========================================
       
   477   delta = 0.00001  
       
   478 
       
   479   if (n .eq. 0) then
       
   480      resg@cnMinLevelValF       = 0.             
       
   481      resg@cnMaxLevelValF       = 10.             
       
   482      resg@cnLevelSpacingF      = 1.
       
   483 
       
   484      data_ob = where(ismissing(data_ob).and.(ismissing(data_mod).or.(data_mod.lt.delta)),0.,data_ob)
       
   485   end if
       
   486 
       
   487   if (n .eq. 1) then
       
   488      resg@cnMinLevelValF       = 0.             
       
   489      resg@cnMaxLevelValF       = 10.             
       
   490      resg@cnLevelSpacingF      = 1.
       
   491 
       
   492      data_ob = where(ismissing(data_ob).and.(ismissing(data_mod).or.(data_mod.lt.delta)),0.,data_ob)
       
   493   end if
       
   494 
       
   495   if (n .eq. 2) then
       
   496      resg@cnMinLevelValF       = 1.             
       
   497      resg@cnMaxLevelValF       = 12.             
       
   498      resg@cnLevelSpacingF      = 1.
       
   499 
       
   500      data_ob = where(ismissing(data_ob).and.(ismissing(data_mod).or.(data_mod.lt.delta)),0.,data_ob)
       
   501   end if
       
   502 
       
   503   if (n .eq. 3) then
       
   504      resg@cnMinLevelValF       = 60.             
       
   505      resg@cnMaxLevelValF       = 360.             
       
   506      resg@cnLevelSpacingF      = 20.
       
   507 
       
   508      data_ob@_FillValue = 1.e+36
       
   509      data_ob = where(data_ob .lt. 10.,data_ob@_FillValue,data_ob)
       
   510 
       
   511      data_mod@_FillValue = 1.e+36
       
   512      data_mod = where(data_mod .lt. 10.,data_mod@_FillValue,data_mod)        
       
   513   end if                                                        
       
   514 
       
   515 ;=========================
       
   516 ; global contour : ob
       
   517 ;=========================
       
   518   
       
   519   plot_name = "global_"+component(n)+"_ob"
       
   520   title     = ob_name
       
   521   resg@tiMainString  = title
       
   522 
       
   523   wks = gsn_open_wks (plot_type,plot_name)   ; open workstation
       
   524   gsn_define_colormap(wks,"gui_default")     ; choose colormap
       
   525 
       
   526   plot = gsn_csm_contour_map_ce(wks,data_ob,resg)   
       
   527   frame(wks)
       
   528 
       
   529   clear (wks)
       
   530   delete (plot)
       
   531 
       
   532 ;============================
       
   533 ; global contour : model
       
   534 ;============================
       
   535 
       
   536   plot_name = "global_"+component(n)+"_model"  
       
   537   title     = "Model " + model_name
       
   538   resg@tiMainString  = title
       
   539 
       
   540   wks = gsn_open_wks (plot_type,plot_name)
       
   541   gsn_define_colormap(wks,"gui_default")     
       
   542 
       
   543   plot = gsn_csm_contour_map_ce(wks,data_mod,resg)   
       
   544   frame(wks)
       
   545 
       
   546   clear (wks)
       
   547   delete (plot)
       
   548 
       
   549 ;================================
       
   550 ; global contour: model vs ob
       
   551 ;================================
       
   552 
       
   553   plot_name = "global_"+component(n)+"_model_vs_ob"
       
   554 
       
   555   wks = gsn_open_wks (plot_type,plot_name)   
       
   556   gsn_define_colormap(wks,"gui_default")     
       
   557 
       
   558   plot=new(3,graphic)                        ; create graphic array
       
   559 
       
   560   resg@gsnFrame             = False          ; Do not draw plot 
       
   561   resg@gsnDraw              = False          ; Do not advance frame
       
   562 
       
   563 ; plot correlation coef
       
   564 
       
   565   gRes               = True
       
   566   gRes@txFontHeightF = 0.02
       
   567   gRes@txAngleF      = 90
       
   568 
       
   569   correlation_text = "(correlation coef = "+sprintf("%.2f", cc)+")"
       
   570 
       
   571   gsn_text_ndc(wks,correlation_text,0.20,0.50,gRes)
       
   572   
       
   573 ; plot ob
       
   574 
       
   575   title     = ob_name
       
   576   resg@tiMainString  = title
       
   577 
       
   578   plot(0) = gsn_csm_contour_map_ce(wks,data_ob,resg)       
       
   579 
       
   580 ; plot model
       
   581 
       
   582   title     = "Model "+ model_name
       
   583   resg@tiMainString  = title
       
   584 
       
   585   plot(1) = gsn_csm_contour_map_ce(wks,data_mod,resg) 
       
   586 
       
   587 ; plot model-ob
       
   588 
       
   589   if (n .eq. 0) then
       
   590      resg@cnMinLevelValF  = -2.           
       
   591      resg@cnMaxLevelValF  =  2.            
       
   592      resg@cnLevelSpacingF =  0.4
       
   593   end if
       
   594 
       
   595   if (n .eq. 1) then
       
   596      resg@cnMinLevelValF  = -6.           
       
   597      resg@cnMaxLevelValF  =  6.            
       
   598      resg@cnLevelSpacingF =  1.
       
   599   end if
       
   600 
       
   601   if (n .eq. 2) then
       
   602      resg@cnMinLevelValF  = -6.           
       
   603      resg@cnMaxLevelValF  =  6.            
       
   604      resg@cnLevelSpacingF =  1.
       
   605   end if
       
   606 
       
   607   if (n .eq. 3) then
       
   608      resg@cnMinLevelValF  = -100.           
       
   609      resg@cnMaxLevelValF  =  100.            
       
   610      resg@cnLevelSpacingF =  20.
       
   611   end if            
       
   612 
       
   613   zz = data_mod
       
   614   zz = data_mod - data_ob
       
   615   title = "Model_"+model_name+" - Observed"
       
   616   resg@tiMainString    = title
       
   617 
       
   618   plot(2) = gsn_csm_contour_map_ce(wks,zz,resg) 
       
   619 
       
   620 ; plot panel
       
   621 
       
   622   pres                            = True        ; panel plot mods desired
       
   623   pres@gsnMaximize                = True        ; fill the page
       
   624 
       
   625   gsn_panel(wks,plot,(/3,1/),pres)              ; create panel plot
       
   626 
       
   627   clear (wks)
       
   628   delete (plot)
       
   629 
       
   630   end do
       
   631 ;**************************************************
       
   632 ; html table
       
   633 ;**************************************************
       
   634   output_html = "table_model_vs_ob.html"
       
   635 
       
   636   header_text = "<H1>LAI: Model "+model_name+" vs Observed</H1>" 
       
   637 
       
   638   header = (/"<HTML>" \
       
   639             ,"<HEAD>" \
       
   640             ,"<TITLE>CLAMP metrics</TITLE>" \
       
   641             ,"</HEAD>" \
       
   642             ,header_text \
       
   643             /) 
       
   644   footer = "</HTML>"
       
   645 
       
   646   table_header = (/ \
       
   647         "<table border=1 cellspacing=0 cellpadding=3 width=100%>" \
       
   648        ,"<tr>" \
       
   649        ,"   <th bgcolor=DDDDDD rowspan=2>Biome Class</th>" \
       
   650        ,"   <th bgcolor=DDDDDD colspan=3>"+component(0)+"</th>" \
       
   651        ,"   <th bgcolor=DDDDDD colspan=3>"+component(1)+"</th>" \
       
   652        ,"   <th bgcolor=DDDDDD colspan=3>"+component(2)+"</th>" \
       
   653        ,"   <th bgcolor=DDDDDD colspan=3>"+component(3)+"</th>" \
       
   654        ,"</tr>" \
       
   655        ,"<tr>" \
       
   656        ,"   <th bgcolor=DDDDDD >observed</th>" \
       
   657        ,"   <th bgcolor=DDDDDD >"+model_name+"</th>" \
       
   658        ,"   <th bgcolor=DDDDDD >M_score</th>" \
       
   659        ,"   <th bgcolor=DDDDDD >observed</th>" \
       
   660        ,"   <th bgcolor=DDDDDD >"+model_name+"</th>" \
       
   661        ,"   <th bgcolor=DDDDDD >M_score</th>" \
       
   662        ,"   <th bgcolor=DDDDDD >observed</th>" \
       
   663        ,"   <th bgcolor=DDDDDD >"+model_name+"</th>" \
       
   664        ,"   <th bgcolor=DDDDDD >M_score</th>" \
       
   665        ,"   <th bgcolor=DDDDDD >observed</th>" \
       
   666        ,"   <th bgcolor=DDDDDD >"+model_name+"</th>" \
       
   667        ,"   <th bgcolor=DDDDDD >M_score</th>" \
       
   668        ,"</tr>" \
       
   669        /)
       
   670   table_footer = "</table>"
       
   671   row_header = "<tr>"
       
   672   row_footer = "</tr>"
       
   673 
       
   674   lines = new(50000,string)
       
   675   nline = 0
       
   676 
       
   677   set_line(lines,nline,header)
       
   678   set_line(lines,nline,table_header)
       
   679 ;-----------------------------------------------
       
   680 ;row of table
       
   681 
       
   682   do n = 0,nrow-1
       
   683      set_line(lines,nline,row_header)
       
   684 
       
   685      txt1  = row_head(n)
       
   686      txt2  = text4(n,0)
       
   687      txt3  = text4(n,1)
       
   688      txt4  = text4(n,2)
       
   689      txt5  = text4(n,3)
       
   690      txt6  = text4(n,4)
       
   691      txt7  = text4(n,5)
       
   692      txt8  = text4(n,6)
       
   693      txt9  = text4(n,7)
       
   694      txt10 = text4(n,8)
       
   695      txt11 = text4(n,9)
       
   696      txt12 = text4(n,10)
       
   697      txt13 = text4(n,11) 
       
   698 
       
   699      set_line(lines,nline,"<th>"+txt1+"</th>")
       
   700      set_line(lines,nline,"<th>"+txt2+"</th>")
       
   701      set_line(lines,nline,"<th>"+txt3+"</th>")
       
   702      set_line(lines,nline,"<th>"+txt4+"</th>")
       
   703      set_line(lines,nline,"<th>"+txt5+"</th>")
       
   704      set_line(lines,nline,"<th>"+txt6+"</th>")
       
   705      set_line(lines,nline,"<th>"+txt7+"</th>")
       
   706      set_line(lines,nline,"<th>"+txt8+"</th>")
       
   707      set_line(lines,nline,"<th>"+txt9+"</th>")
       
   708      set_line(lines,nline,"<th>"+txt10+"</th>")
       
   709      set_line(lines,nline,"<th>"+txt11+"</th>")
       
   710      set_line(lines,nline,"<th>"+txt12+"</th>")
       
   711      set_line(lines,nline,"<th>"+txt13+"</th>")
       
   712 
       
   713      set_line(lines,nline,row_footer)
       
   714   end do
       
   715 ;-----------------------------------------------
       
   716   set_line(lines,nline,table_footer)
       
   717   set_line(lines,nline,footer) 
       
   718 
       
   719 ; Now write to an HTML file.
       
   720   idx = ind(.not.ismissing(lines))
       
   721   if(.not.any(ismissing(idx))) then
       
   722     asciiwrite(output_html,lines(idx))
       
   723   else
       
   724    print ("error?")
       
   725   end if
       
   726 
       
   727 ;***************************************************************************
       
   728 ; write total score to file
       
   729 ;***************************************************************************
       
   730 
       
   731   asciiwrite("M_save.lai", M_total)
       
   732 
       
   733 ;***************************************************************************
       
   734 end
       
   735